Introduction

Assessing dynamic material properties from scattered elastic waves is at the center of numerous remote sensing applications, including acoustic and elastic imaging, non-invasive diagnostics, and metamaterial design and characterization. In most cases, the wave equation is inverted using numerical or analytical methods that rely on limiting assumptions. For example, all the components of the mass density and stiffness tensors of unknown materials could be extracted from free-space scattered acoustic waves, but these methods require material samples large enough so that the diffraction from the sample edges does not contaminate the scattered field1,2.

In many scenarios of interest (e.g., medical diagnostics and metamaterial design), the samples are small and the aforementioned analytical free-space methods do not apply. For instance, the mechanical properties of biological tissue were shown to be modified by pathological processes. State-of-the-art non-invasive techniques, including impediography and elastography ensonify biological tissue in vivo with ultrasound beams and infer a small subset of parameters such as impedance, Young’s or shear moduli of potentially small abnormal regions from the scattered ultrasound3,4. However, biological tissue is a complex anisotropic material characterized by numerous distinct stiffness and mass density tensor components. Extracting all of these components (not just a small subset) is critical for assessing tissue health and diagnosing possible abnormalities early on while the affected area is small.

Even though diffracted waves potentially contain important information about the elastic behavior of materials, making sense of these complicated fields is still not well understood, and diffracted fields are usually avoided5,6,7,8. Diffraction-avoiding techniques have recently been employed extensively in metamaterial research. Metamaterial design validation is typically done by fabricating small prototypes whose dynamic behavior is measured in a series of experiments in which waveguides are used to mitigate the influence of edge diffraction on the scattered fields9. These methods involve placing the sample under test inside carefully designed waveguides and inverting reflection and transmission coefficients to obtain material properties such as dynamic mass density and various elastic moduli10,11,12,13,14. However, these methods are either unsuitable in applications in which the sample cannot be manipulated (e.g., in non-invasive diagnostics, imaging) or cumbersome to apply (e.g., for materials operating in aqueous environments in which the waveguide itself contaminates the scattered fields)15.

Instead of mitigating the effects of diffraction, this work shows that convolutional neural networks (CNNs) can effectively interpret the diffracted fields to extract important information about the dynamical behavior of complex media, as follows. First, we show that CNNs can efficiently estimate the mapping between scattered acoustic fields and macroscopic material parameters from relatively few examples of numerically simulated near-field spatial distributions scattered by material samples having randomly chosen macroscopic parameters. More importantly, we show that CNNs trained exclusively with self-labeled synthetic data can process acoustic field measurements and accurately estimate the effective material parameters of physical samples. Second, we introduce methods to analyze the trained neural networks to obtain quantitative measures of how sensitive the scattered fields are to each macroscopic property. These methods could be leveraged in material design to prioritize the realization of one material property over another. Another consequence of this sensitivity analysis is obtaining an upper boundary for the shear modulus below which a solid material behaves acoustically as a fluid. We will show that this finding has significant relevance to transformation acoustics applied to the design of metamaterials operating in aqueous environments16,17,18. Third, we show that further analysis of the trained network reveals regions of the scattered fields that contain the most information on the macroscopic material parameters. These regions located around the edges of the sample under test provide direct evidence that fields diffracted by the sample edges contain significant information about the sample’s behavior under elastic wave excitation.

Results and discussion

Approach

Figure 1a shows a typical scenario in which a material sample is ensonified by a source situated in the far field. The sample is small enough that the sound diffracted by the material’s edges dominates the scattered acoustic field. Our goal is to train a convolutional neural network that estimates all the elastic material parameters of the sample, including its effective dynamic stiffness tensor and mass density and to analyze the trained CNN to obtain insight into the physical behavior of matter excited by mechanical waves.

Fig. 1: Material characterization from scattered fields.
figure 1

a The unknown sample is ensonified by a source and the scattered fields are recorded on the reflection and transmission planes located at distances d1 and d2 in front and, respectively, behind the sample. Two sources labeled 1 and 2 are considered in this paper; b The scattered fields for an aluminum sample were computed numerically and show complicated patterns caused by the interference of edge-diffracted waves.

For illustration purposes and without loss of generality, we assume that the sample has a rectangular shape of width W, height H, and thickness T and is placed in a background fluid of known mass density ρ0 and bulk modulus K0. The acoustic pressure scattered by the sample is recorded in front and behind it on planes situated at known distances d1 and d2 away from the sample, which we call reflection and, respectively, transmission planes.

To highlight the difficulties imposed by edge diffraction effects, Fig. 1b shows the scattered pressure distribution (amplitude and phase) obtained for an aluminum block of dimensions W = H = 4λ and T = 0.16λ, where λ is the wavelength of the omnidirectional source ensonifying the block at 120 kHz. The fields were computed numerically on the reflection and transmission planes using Comsol Multiphysics, a finite element solver of the acoustic wave equation. The simulation assumed a water background (ρ0 = 1000 kg m−3, K0 = 2.15 GPa) and the acoustic pressure fields are reported relative to the incident spherical wave. The superposition between the diffracted, reflected, and transmitted fields give rise to complex patterns in both the scattered pressure amplitude and phase. This complicated field distribution makes inverting the scattered pressure very challenging using analytical approaches1,2. Nevertheless, the pattern details are very sensitive to changes in the macroscopic material parameters inside the material, which is critical to the success of our CNN-based approach.

Our aim is to find the mapping \(g:{{{{{{{\mathcal{F}}}}}}}}\to {{{{{{{\mathcal{M}}}}}}}}\) between the set of scattered fields \({{{{{{{\mathcal{F}}}}}}}}\) obtained on the reflection and transmission planes and the set of macroscopic material parameters \({{{{{{{\mathcal{M}}}}}}}}\). An element \(f\in {{{{{{{\mathcal{F}}}}}}}}\) may be represented by the four images shown in Fig. 1b. An element m in the N-dimensional set \({{{{{{{\mathcal{M}}}}}}}}\) comprises all the N unknown macroscopic properties of the sample. In general, m may contain components of the stiffness, mass density, and Willis parameter tensors.

The inverse mapping \({g}^{-1}:{{{{{{{\mathcal{M}}}}}}}}\to {{{{{{{\mathcal{F}}}}}}}}\) can be readily computed with numerical solvers such as Comsol Multiphysics or through measurements. We will solve the much more difficult problem of determining g by computing g−1(m) for \(m\in {{{{{{{{\mathcal{M}}}}}}}}}_{T}\subset {{{{{{{\mathcal{M}}}}}}}}\), where \({{{{{{{{\mathcal{M}}}}}}}}}_{T}\) represents a small set of examples called training set. We will show that CNNs taking as input the set of scattered field patterns (e.g., see Fig. 1b) and outputting the material parameters m will extrapolate this set of examples to find g(f) for any \(f\in {{{{{{{\mathcal{F}}}}}}}}\).

CNNs are ideal at finding patterns in images such as f and associating them to the scatterer identity encoded in the effective macroscopic material parameters m. This pattern recognition problem is very similar to a standard computer vision task in which CNNs are asked to recognize familiar objects in photos19. The resurgence of CNN research during the last decade is in part the result of the remarkable success of neural networks at solving this type of pattern recognition problem. Moreover, research on CNNs produced methods to analyze the trained algorithms to obtain important insight into the workings of the pattern recognition process20.

This suggests that modified versions of the computer vision CNNs should be able to solve our pattern recognition task and find the mapping g. More importantly, it suggests that analysis of the trained CNN may provide important insight into the physics of the sample in addition to the material properties we are seeking to find.

Application to elastic metamaterials

We illustrate the capability of the CNN in an example in which an isotropic metamaterial designed to operate underwater is composed of copper spheres 2 mm in diameter placed in a polylactide (PLA) frame, as depicted in Fig. 2. The frame is 3D-printed and the thickness of each frame wall is 0.2 mm. Unlike acoustic or electromagnetic media, elastic metamaterials with square lattices and unit cells with C4 rotational symmetry could be anisotropic, especially when the effective shear modulus is high21. However, the anisotropy level is typically much smaller than the ~5% accuracy of our method. In addition, we expect that the effective shear modulus of our metamaterial is modest since the components of the structure are not glued together but held in place by pressure fitting the copper spheres inside the plastic matrix. Therefore, the metamaterial shown in Fig. 2 is under a very good approximation isotropic.

Fig. 2: Metamaterial sample.
figure 2

The sample consists of copper spheres embedded in a 3D-printed polylactide (PLA) matrix. The inset shows a zoomed-in section of the metamaterial featuring the placement of spheres inside the plastic matrix.

Solid metamaterials functioning in dense backgrounds are notoriously difficult to design and characterize. Unlike the vast majority of cases reported in the literature involving metamaterials in air, in aqueous environments, mechanical waves propagate through the solids composing the metamaterial, and parameters such as the friction force between the spheres and the plastic matrix keeping the spheres in place significantly influence the effective macroscopic properties of the sample but are very hard to measure. In these situations, methods that predict macroscopic properties based on the material microstructure22,23,24,25,26 become very challenging to apply.

For instance, the shear modulus has been shown to be very low when the friction force between the solid inclusions composing a metamaterial is small17,27 with negative effects on the mechanical integrity of the sample. Conversely, when the friction force is large (in the limit, the spheres in Fig. 2 are fused to their frame), the shear will be in the gigapascals range. A question of interest becomes whether mechanically robust structures in which the friction force is large could be used in underwater transformation acoustics devices that rely on zero shear. Remarkably, we will show shortly that analysis of the trained CNN provides definitive answers.

The sample shown in Fig. 2 has dimensions W = H = 4λ and T = 0.16λ at the operating frequency of 120 kHz and is characterized by N = 4 unknown material properties. These are the mass density ρ, the bulk modulus K, the shear modulus G, and the stiffness loss tangent factor η. The stiffness tensor can be determined from these parameters as (1 + iη)C(K, G), where C is the stiffness tensor of lossless isotropic solids having the elasticity moduli K and G, and i is the imaginary unit.

In this work, our CNN is a modified version of AlexNet, a neural network that proved particularly effective at pattern recognition in images19. These computer vision tasks are significantly more difficult than ours for several reasons. First, the process of sound scattering is linear, and there is significant redundancy in the spatial distribution of the scattered fields. Second, there is no overlap between the sound scattered by the sample under test and other nearby objects. Third, the sample is ensonified from the same angle in all cases using the same known source. All these factors suggest that significantly less training data is needed in this work than in the computer vision image classification problem, which points to employing a simpler network having fewer unknown weights rather than the much deeper networks used in computer vision. In our case we experimented with various network depths and converged on the following architecture.

Our network consists of four convolutional layers followed by two fully connected layers, as described in Supplementary Note 1. The network is trained with maps f = g−1(m) obtained in numerical simulations performed with Comsol Multiphyics, where \(m\equiv (\rho ,K,G,\eta )\in {{{{{{{{\mathcal{M}}}}}}}}}_{T}\) and \({{{{{{{{\mathcal{M}}}}}}}}}_{T}\) contains 8000 material parameter quadruples sampled randomly in a very large space chosen as follows (see Supplementary Note 2 for training details and Supplementary Note 3 for the reason why 8000 material parameter quadruples are selected). The dynamical mass density ρ is searched in the range 400–6000 kg m−3. These are very conservative bounds because the effective dynamical ρ is likely to be lower than the static density of ~3000 kg m−3 since more sound propagates through the plastic matrix (same impedance as water) than copper. Similarly, the bulk modulus was sampled in the interval 0.5–12 GPa (K0 = 2.15 GPa for water). There is currently little knowledge and research on the dynamical shear modulus of metamaterials, which is reflected in the larger searching interval for G, namely 30 kPa–50 GPa. This search space was chosen to encompass the shear of copper (45 GPa). There is no reason to look below 30 kPa because a material with such a low shear behaves like a true fluid (G = 0) for all practical purposes. In fact, we will show that the analysis of the trained CNN reveals that shear below 50 MPa has a negligible effect on the elastic behavior of isotropic media. The loss parameter η is searched in the interval 0 (no loss) – 0.5 (high loss).

The training set was obtained by randomly sampling the values of ρ, \({\log }_{10}(K)\), \({\log }_{10}(G)\), and η. As explained above, the elastic moduli (bulk and shear) are searched in very large intervals spanning multiple orders of magnitude, therefore their logarithm is sampled uniformly instead of their linear values.

The sample is ensonified by an omnidirectional point source located at coordinates x = 0, y = −7.1 cm, and z = 18.5 cm relative to the center of the metamaterial [see source location 1 in Fig. 1a]. The simulation domain is surrounded by reflectionless, perfectly matched layers (PML), and the quasi-spherical incident wave is implemented as a background pressure domain condition. The acoustic pressure is recorded on 8λ by 8λ reflection and transmission planes situated d1 = 2.7λ and respectively d2 = 1.5λ from the metamaterial. The values of d1 and d2 were chosen based on the experimental setup, which will be discussed later. The amplitude and phase recorded on these planes are used as inputs f to the network (each image is a 51 × 51 matrix). All simulated fields are normalized to the incident pressure.

The performance of the CNN model trained with the material parameter set \({{{{{{{{\mathcal{M}}}}}}}}}_{T}\) (i.e., labels) and the corresponding scattered fields \({g}^{-1}({{{{{{{{\mathcal{M}}}}}}}}}_{T})\) (i.e., CNN inputs) is quantified using a validation set \({{{{{{{{\mathcal{M}}}}}}}}}_{V}\) and its corresponding scattered fields \({g}^{-1}({{{{{{{{\mathcal{M}}}}}}}}}_{V})\), where \({{{{{{{{\mathcal{M}}}}}}}}}_{T}\) and \({{{{{{{{\mathcal{M}}}}}}}}}_{V}\) do not have any elements in common. For each ground truth material parameter quadruple \(m\equiv (\rho ,K,G,\eta )\in {{{{{{{{\mathcal{M}}}}}}}}}_{V}\) we feed the corresponding field distribution f = g−1(m) to the CNN to obtain an estimation of these material parameters \(\hat{m}\equiv (\hat{\rho },\hat{K},\hat{G},\hat{\eta })=g(f)\).

Physical insights provided by CNN analysis

Figure 3 compares the ground truth m versus the CNN-estimated \(\hat{m}\) for all the elements in the validation set. For each material parameter, the figure provides the R228 and SMAPE metrics29. The comparisons show that all the four unknown material parameters are very well estimated by the CNN. One notable exception is the shear modulus G which deserves a longer discussion, as follows.

Fig. 3: Validation set analysis.
figure 3

The analysis shows with solid blue circles the CNN-estimated material properties, \(\hat{m}\), versus their real values (ground truth), m, where m represents \({\log }_{10}(K)\), \({\log }_{10}(G)\), ρ, and η. The \(\hat{m}=m\) ground truth line is shown with a solid magenta line. The R2 and SMAPE metrics of each parameter are listed in the lower right corner of each panel.

The validation set analysis may be interpreted in two ways. In the first interpretation, Fig. 3 provides confidence bounds on \(\hat{m}\). For example, the estimation error for ρ, η and \({\log }_{10}(K)\) is ±50 kg m−3, ±0.05, ±0.1 in all cases. Interestingly, the CNN appears to fail to evaluate the shear modulus correctly for values below ~50 MPa. To understand the causes of this behavior we need to consider the second interpretation of the plots \(\hat{m}\) versus m.

In the second interpretation, Fig. 3 quantifies the influence of each material property on the scattered near-field. The mass density has the most influence on the scattered field distribution and therefore \(\hat{\rho }\) is very close to its real value for all the examples in the validation set. At the opposite spectrum lies the shear modulus. For G < 50 MPa, the estimation \({\log }_{10}(\hat{G})\) tends to converge to a constant value.

This behavior reveals an important physical insight into the relationship between the shear modulus and the scattered field it generates. The training process minimizes the mean-square-error loss function \(\parallel \hat{m}-m{\parallel }^{2}\) to evaluate how good the CNN estimation is and to update the weights of the convolutional and fully connected layers. The neural network detects no effect of G on the scattered fields for low values of G and therefore chooses the estimations \({\log }_{10}(\hat{G})\) on a line situated halfway in between the minimum and maximum values of \({\log }_{10}(G)\) to minimize the loss function over the entire validation set. This important conclusion is further validated in Supplementary Note 4 with numerical simulations that show quantitatively how the scattered fields change as the shear modulus increases from 0 to higher values.

This interpretation is an especially remarkable result in the context of transformation acoustics (TA). Unlike TA devices operating in air, TA devices functioning in dense environments (e.g., water) require fluid-like media having G = 0, which is challenging to obtain in structures composed of solids. It has been shown that metamaterials composed of loosely connected solids30,31, soft rubbers16,18, and oiled interfaces17,27 might work but these approaches typically result in flimsy devices lacking structural robustness. The validation set analysis reveals that a metamaterial with an effective G < 50 MPa behaves acoustically indistinguishable from a fluid and thus extremely low shear modulus is not needed.

Having established the ability of the CNN to estimate all the unknown material parameters, we assess the ability of the neural network to process measured fields and extract the effective material parameters of the fabricated metamaterial shown in Fig. 2. The metamaterial is hung by a 0.1 mm diameter copper wire at the center of a 0.5 m × 0.5 m × 0.5 m water tank. A hydrophone connected to a Verasonics Vantage 64 research ultrasound machine is placed at the source location 1, shown in Fig. 1a, having coordinates (0, −7.1 cm, 18.5 cm). The hydrophone sends quasi-spherical waves in the form of five-cycle Gaussian pulses centered at 120 kHz. A second hydrophone scans the reflection and transmission planes shown in Fig. 1a and records the reflected and transmitted waves, which are then normalized to the incident field to replicate the conditions simulated in Comsol Multiphysics and used to train the CNN. The distances between the reflection and transmission planes are set to d1 = 3.4 cm and, respectively, d2 = 1.9 cm. These distances were chosen to allow the scanning hydrophone to move without touching the sample.

The measured normalized reflected and transmitted fields are shown in Fig. 4a. These scattered fields are fed to the CNN, which produced the following estimated material properties: \(\hat{K}=2.5\) GPa, \(\hat{G}=0.72\) MPa, \(\hat{\rho }=1285.1\) kg m−3, and \(\hat{\eta }=0.252\). These values are validated by simulating in Comsol Multiphysics the acoustic fields scattered from a continuous homogeneous material having these estimated properties and comparing the measured scattered fields against the simulations. Figure 4 shows the excellent agreement between measurements and simulations. Remarkably, every single feature in the measured acoustic pressure is captured by the simulation. To quantify the similarity between the measured and simulated scattered fields, we defined a figure of merit analogous to R2 to better quantify the match. In this case, we defined R2 for the transmission and reflection planes as \({R}^{2}=1-\parallel S-M{\parallel }_{F}^{2}/\parallel S-\overline{S}{\parallel }_{F}^{2}\), where S and M are the simulated and respectively measured spatial 2D fields obtained on either the reflection or the transmission planes, \(\overline{S}\) is the spatial average of the simulated field, and F is the Frobenius norm, i.e., the Euclidean norm of a matrix. We obtained R2 = 0.92 on the transmission plane and R2 = 0.85 on the reflection plane.

Fig. 4: Measured versus simulated scattered acoustic fields.
figure 4

Comparisons between the measured (a) and simulated (b) scattered fields on the transmission and reflection planes show identical features in both amplitude and phase. The scale of each field quantity (relative amplitude and phase) is shown underneath each plot column.

The accuracy of the determined macroscopic properties was further confirmed through numerical simulations in which only the mass density was modified from the estimated value by several percentages while all the other material parameters were kept unchanged (see Supplementary Note 5 and Supplementary Fig. 5). We also confirmed the correctness of the evaluated macroscopic properties by comparing the measurements and simulations for a different position of the source relative to the metamaterial sample, namely source location 2 in Fig. 1a having coordinates (0, 0, 18.5 cm). The comparison between the measured and simulated acoustic pressure for this ensonification is presented in Supplementary Note 5 and Supplementary Fig. 6 and shows the same excellent match.

All these comparisons demonstrate the robustness of the CNN at extracting macroscopic effective properties in experiments. There are small discrepancies between the measurements and simulations caused by the misalignment of the submerged sample inside the water tank. For instance, the sample was not perfectly vertical and was also rotated slightly about the normal to its surface. As a result, the acoustic fields shown in Fig. 4 are not symmetric as they would be in the ideal case. Nevertheless, the sound diffracted by the sample edges generates rich features in the scattered fields that do not change because of misalignment and other experimental errors, which explains the CNNs robustness to experimental errors.

In this work, we expressed the stiffness tensor in terms of bulk and shear moduli. Other stiffness quantities could also be estimated by the CNN, such as Young’s modulus and Lamé parameters. For instance, Supplementary Note 6 shows the CNN-estimated Young’s modulus, which is consistent with the estimated K and G presented in Fig. 3.

Furthermore, we analyze the trained CNN to extract more information on the dynamical behavior of materials. In particular, we are interested in understanding what sub-regions of the scattered field contain the most information on the identity of the sample. To address this issue, we occluded different portions of the input testing image with a 3 × 3 mask of pixels, i.e., the same size as the convolutional kernels used by the CNN. The mask is moved in all four input images by steps of 1 pixel, and the changes in the CNN output were recorded. Figure 5 shows the average change of each material parameter on a logarithmic scale for the reflection and transmission sides for the entire validation data. Larger values indicate regions targeted by the CNN that most influence the CNN outputs.

Fig. 5: Field regions targeted by CNN.
figure 5

Average change of CNN estimations for each material parameter when a 3 × 3 pixel mask moves and occludes the scattered fields (logarithmic scale). Higher changes indicate regions targeted by the CNN. The area inside the dotted rectangles is the specular projection of the material sample onto the a reflection and b transmission planes.

Figure 5 shows that the bulk modulus is determined mostly from the fields reflected and transmitted through the sample because the area that most influences this parameter is inside the specular projection of the sample in the measurements planes marked with dotted rectangles in the figure. The other three material parameters are mostly determined from the diffracted fields, as demonstrated by regions of high influence on these parameters located around the edges of the material. This constitutes direct evidence of the importance of diffracted fields on the estimated material parameters.

Conclusions

We presented a CNN-based approach to construct the direct map of scattered acoustic fields-to-elastic macroscopic material properties. The CNN is trained exclusively with self-labeled data obtained in numerical simulations but can robustly process measurements. More importantly, analysis of the trained CNN reveals important insight into the dynamic behavior of elastic materials. For example, it reveals quantitative bounds for the shear modulus for which the sample behaves acoustically like a fluid, which is very relevant to the design of underwater transformation acoustics devices. In addition, it shows the regions around the metamaterial where the acoustic pressure should be sampled to extract most information on the metamaterial effective macroscopic properties. This analysis produced direct evidence of the importance of diffracted fields to characterize unknown media from scattered fields. Furthermore, validation set analysis indicates confidence bounds on the extracted material parameters. Although our method is illustrated for an isotropic elastic metamaterial, the CNN approach shown here is general in that it could, in principle, be applied to extract unknown material properties for a wide range of complex media from biological tissue to Willis media and anisotropic materials. This work also shows that neural network-based approaches may become an important tool to characterize and analyze the dynamical behavior of complex media.

Method

Experimental apparatus

The experimental measurements are done inside a water tank composed of transparent polycarbonate plates 50 cm × 50 cm × 0.6 cm that are held together with T-slot aluminum rails. The metamaterial sample measured in the experiment is made up of a 3D-printed PLA frame with 22 × 22 slots and the matching number of copper spheres. The height and the width of the PLA frame are both 5 cm, the wall forming the slots is 0.2 mm in thickness, and the diameter of all copper spheres is 2 mm. The overall thickness of the metamaterial is the same as the copper sphere diameter. A Reson TC 4013 hydrophone connected to a Verasonics Vantage 64 research ultrasound machine serves as the acoustic source and produces quasi-spherical waves towards the metamaterial. The spatial and temporal distribution of the acoustic field is measured by a second Reson TC 4013 hydrophone connected to the same ultrasound machine that raster scans the area in front of and behind the metamaterial sample. The frequency spectrum of the measure acoustic field is obtained from the Fourier transformation of the time domain measurements.

Acoustic fields simulations

Comsol Multiphysics is employed to obtain all simulated acoustic fields. The simulation has the same setup as the experiment. Since the water tank and the metamaterial are symmetric, the simulation domain and metamaterial defined in Comsol Multiphysics only contain parts on one side of the y-z plane shown in Fig. 1a to reduce computation complexity. The y-z plane of the simulation has the symmetry boundary condition, and the rest of the domain is surrounded by reflectionless, perfectly matched layers (PML). The quasi-spherical incident wave is implemented as a background pressure domain condition. Comsol Multiphysics is integrated with MATLAB via the LiveLink for MATLAB, which allows us to control the simulation from MATLAB code. The material properties of the sample are selected randomly and automatically via the MATLAB code while running different simulation cases.