Elastic least-squares reverse-time migration with density variation for flaw imaging in heterogeneous structures

Ultrasonic techniques are able to accurately detect and characterize flaws in homogeneous structures. Elastic reverse time migration (ERTM) is a powerful tool to reconstruct high-resolution images of flaws. To achieve images with better quality, the solution can be obtained by iteratively finding an image generating the modeled data which can best match the measured data in a least-squares sense, i.e. least-squares migration (LSM). Combing ERTM and LSM, conventional elastic least-squares reverse time migration (ELSRTM) methods are based on the assumption of a constant density, which can lead to inaccurate amplitudes and parameter crosstalk artifacts in the reconstructed images. In this paper, an ultrasonic imaging technique based on the ELSRTM which considers density as well as longitudinal-(L-) and shear-wave (S-wave) velocity variations is explored for imaging flaws in heterogeneous structures. The ELSRTM with density variations can simultaneously reconstruct density and L- and S-wave velocity images, which can provide amplitude-preserving images and mitigate crosstalk artifacts. This method is applied to numerical as well as physical laboratory experiments and the results appear promising for flaw identification in heterogeneous structures.


Introduction
Monitoring of the state of engineered structures is important for mitigating the risk and reducing the repair cost. To provide the early warning of structural flaws, such as the deterioration in inhomogeneous austenitic welds [1] or in concrete bridge decks [2], ultrasonic non-destructive evaluation (NDE) techniques are increasingly useful tools. Ultrasonic flaw sizing techniques can be broadly grouped into the following categories: amplitude, temporal, imaging and inversion [3]. In amplitude techniques, the amplitude of a scattered signal from the flaw and other knowledge are used to infer the flaw size. Temporal techniques use the arrival time of signals from the flaw to obtain the flaw size, e.g. time-offlight diffraction [4,5] and relative arrival time technique [6]. Imaging techniques use signals from the flaw to obtain the representation of the region of interest and its size is inferred from this representation in different ways. Ultrasonic images are often obtained by using arrays, which can capture the full data set and post-process them for B-scan [7], sectorial scan [8] and focused [9,10] images. In inversion techniques, signals from the flaw are inputted into various algorithms which determine the physical properties of the flaw, such as scattering matrix techniques [11] and full waveform inversion [12].
In many cases, a structure is simple enough for arrival times of signals to be accurately computed by traveltimebased methods. However, these traveltime-based techniques are not easy to deal with structures of complex geometry and characterize mechanical properties of the anomalous feature [13]. Advanced ultrasonic techniques usually involve full acoustic or elastic wave equation modeling, which can offer solutions to these limitations. The time reversal technique based on self-focusing of the reversed wavefield has been applied to image acoustic emission sources [14][15][16]. Based on the temporal and spatial cross-correlation of the forward and backward wavefields, the reverse time migration (RTM) algorithm has shown the ability to locate and detect interfaces and flaws [17]. The acoustic and elastic RTM were first developed in geophysics for seismic waves imaging [18,19] and have been applied to characterize flaws in NDE applications [20][21][22][23]. Compared with the acoustic RTM, elastic reverse time migration (ERTM) provides a better description of wave-propagation phenomena and achieves more accurate images [19]. However, conventional ERTM methods for detection of flaws are based on the assumption of a constant density [24][25][26][27]. Although the density does not play an important role as longitudinal (L-) and shear-wave (S-wave) velocities on signal waveforms, its contribution to wave amplitudes cannot be ignored [28,29]. For example, scatterers with density variations can generate L-and S-wave scatterings in the elastic case [30]. Therefore, the neglection of density variations in the ERTM algorithm can lead to inaccuracies in L-and S-wave velocity images.
In this paper, an ultrasonic imaging algorithm based on the ERTM in the presence of density variations for characterizing flaws in heterogeneous structures can simultaneously consider density as well as L-and S-wave velocity variations, and analyze the influence of density variations on velocity reconstruction results. This approach can also provide L-and S-wave impedance (i.e. the product of velocity and density) images by summing estimated images of L-and S-wave velocities with density, respectively [29,31]. To reduce the artifacts caused by limited transducer coverages and improve the quality of the image, least-squares migration (LSM) is proposed to match the amplitudes of the modeled data with the measured data via an iterative inversion scheme [32]. Combining this ERTM and LSM, the elastic leastsquares reverse time migration (ELSRTM) with density variations is investigated in this paper, which can provide amplitude-preserving images and more effectively mitigate density-related crosstalk compared with conventional methods. Therefore, this method can be expected to generate more accurate characterization of flaws in heterogeneous structures. This paper is organized as follows. The theory of ELSRTM with density variations is described. Then, the effect of density perturbations on velocity reflectivity estimations using ultrasonic bulk waves is analyzed. Numerical experiments of two randomly-distributed inhomogeneous materials and a layered structure are demonstrated to validate the effectiveness of this method, and the results are compared to those using the conventional total focusing method (TFM) algorithm. The physical laboratory experiment performed in the multi-material part fabricated using laser powder bed fusion is also presented. The discussion follows and finally conclusions are summarized.

Theoretical formulation for ELSRTM with density variations
The ERTM images flaws of a structure by numerically backprojecting waves measured at the receivers. The principle of ERTM is that the scatterers are points where the forward wavefield and the backward wavefield coincide [33,34]. Based on this, a reflectivity image (explained later) can be obtained by cross-correlating the forward and backward wavefields [26].
Two numerical solutions of the elastic wave equations for heterogeneous media, which can consider variable density as well as variable velocities, are required to produce a reflectivity image. One solution represents a wavefield propagating forward in time from the source and the other represents the calculation of the measured data back-propagated in time from the receivers. These elastic wave equations can be given as [35] r r where U x and U y denote the horizontal-and vertical-displacement wavefields, which can involve L-waves, S-waves and wave mode conversions between these two. V L and V S are L-wave velocity and the S-wave velocity, respectively. ρ is the density and t is time.
The free-surface boundary conditions are that the normal stress is zero [19] - and the tangential stress is zero In this paper, the parameter reflectivity models, which are defined as relative parameter perturbations, can be given by d where δ represents the model perturbation and the subscript 0 refers to the background model. In the context of LSM, the variables in equations (1) and (2) are related to the matrixvector operation (forward modeling): where d mea represents the horizontal and vertical components of measured displacement wavefields and F is an elastic forward modeling operator which computes the perturbed wavefield for a source term (forward wavefields). Then, the migration imaging m mig is the adjoint of the forward modeling [36], which can be expressed as where F T is the adjoint operator that computes the measured wavefield d mea back-propagated in time from the receivers (backward wavefields). This adjoint operator transforms the measured data into the image m mig . Based on equations (6) and (7), the migration image can be rewritten as where F T F denotes the Hessian matrix H in the LSM [37].
Equivalently, the reflectivity model m can be estimated by using the following equation However, the direct calculation of the inverse of the Hessian matrix H −1 is usually too expensive [38]. Yet, the solution of m can be efficiently obtained by minimizing the objective This objective function can evaluate the misfit between measured data d mea and the modeled data Fm in a leastsquares sense for all sources.
In order to solve the least-squares problem in equation (10), the conjugate gradient (CG) method is used in this paper. More detailed reviews of CG method can be found in [27,39]. The reflectivity model is updated by: where k is the iteration number, α is the step length for the model update, and q is the searching direction of the optimization process.

Influence of density variation on velocity reflectivity estimation
The aim of this section is to investigate the influence of the density variation on the velocity reflectivity estimation using ultrasonic bulk waves. In the heterogeneous elastic case, the measured data contain information of L-wave velocity, S-wave velocity and density. However, in the conventional ELSRTM, only L-and S-wave velocities are included in the Hessian matrix H [40], and thus this method cannot handle the density-related crosstalk. As a result, the neglection of density variations will result in inaccuracy in reflectivity models of L-and S-wave velocities.
To illustrate the effects of ignoring density variations on velocity images, the ELSRTM with and without the consideration of density variations is performed on a heterogeneous structure. Consider an example in figure 1. The L-and S-wave velocities are constant but the true density model includes a circular flaw with 5% negative anomaly, as shown in figures 1(a)-(c). It means that the measured data are only generated by density variations. Figures 1(d) and (e) show significant crosstalk artifacts in L-and S-wave velocities based on the conventional ELSRTM algorithm. The density image cannot be achieved by the conventional ELSRTM algorithm, as shown in figure 1(f). In figures 1(g)-(i), artifacts are almost invisible in velocity images and the density perturbation is well reconstructed by using ELSRTM with the consideration of density variations. These numerical results demonstrate that density variations have a nontrivial effect on velocity images. It needs to be noted that a circular transducer array is around the inspection area in this case in order to better show the results. In this paper, if not otherwise specified, a linear transducer array is placed on the top surface of a structure.

Estimation of impedance reflectivity
In the elastic case, the L-and S-wave impedances are expressed as [29] The relative perturbations of L-and S-wave impedances can be given as  presented in figures 2(d)-(h), respectively. It can be seen that the reflectors of L-and S-wave velocities as well as density are reconstructed at correct positions. The ellipses show interfaces that are caused by density variations. The images using the conventional ELSRTM method suffer from parameter crosstalk artifacts and the density reflector is not part of the velocity images and is wrongly presented in figures 2(d), (e). Compared with the conventional ELSRTM images, images obtained from ELSRTM with density variations (figures 2(f)-(g)) can clearly mitigate the false density reflector in the velocity reconstructions and correctly represent the density reflector in the density image ( figure 2(h)).
The true L-and S-wave impedance reflectivity images that are obtained by summing the true reflectivity models of L-wave velocity ( figure 2(a)) and S-wave velocity (figure 2(b)) with density ( figure 2(c)) are shown in figures 3(a) and (b), respectively. The corresponding images using ELSRTM in the presence of density variations of L-and S-wave impedance reflectivity estimations are shown in figures 3(c) and (d). The underestimations or overestimations in the images of velocity and density can be compensated in the impedance images. Therefore, the L-and S-wave impedance reflectivity estimations show good consistency with the true impedance images, and they are better than the individual velocity and density images.

Flaw in inhomogeneous materials
Material variations in engineered structures can occur, e.g. due to long-term chemical and mechanical degradation. In order to demonstrate the applicability of the ELSRTM in the presence of density variations, the structures with different degrees of inhomogeneous material properties and a 0.  (details in figures 4(b) and (f)). Both cases are simulated by 2D finite element models implemented in POGO [41]. A linear transducer array has 64 elements on the top surface, as shown in figure 4(a). The signal excited by the array is a 2-cycle Hann-windowed tone-burst with a central frequency of 2 MHz. Each element generates the ultrasonic signal and all elements measure reflected signals, which produce all source-receiver combinations in a full matrix capture (FMC) of data. The homogeneous background models (V L =6090 m s −1 , V S =3195 m s −1 and ρ=4000 kg m −3 ) are used as the forward models. Figures 5(a)-(c) show the ELSRTM reconstructions of L-wave, S-wave and density in the weak inhomogeneous structures (the delamination located in the area of material degradation up to 10% compared to background materials) after 30 iterations, respectively. Laplacian filtering is used to suppress low-wavenumber artifacts. It is clear that the shape and the location of the delamination are reconstructed successfully in all images. Accurate sizing of this delamination is obtained. This is because the delamination is located in the weak inhomogeneous area and the scattering from this area has negligible effects on the reconstruction of the delamination. The L-wave and S-wave impedance reflectivity estimations, as shown in figures 6(a) and (b), are obtained by summing the reconstructed velocity and density images in As a comparison, the conventional TFM algorithm is used to the same simulated data, and the result is shown in figure 6(c). It is clear that the TFM result shows a good indication of this delamination, however, it is not as clearly visible as it is in the ELSRTM result ( figure 6(a)). Also, the profiles appear 'blurred' and take up some pixels; this is due to the Point Spread Function which means that even an infinitesimal point will have an image indication that is finite in size. Furthermore, only L-wave information is considered in the conventional TFM. Figures 7(a)-(c) show the reconstructed images of L-wave velocity, S-wave velocity and density in the strong inhomogeneous area (the delamination located in the area of material degradation up to 35%) after 30 iterations, respectively. It can be seen from figure 7 that the thin delamination present in this inhomogeneous area is visible in each parameter image. However, some distortions and shorter lengths in the reflectivity focus for the delamination and artifacts especially around the delamination are observed, because the velocity and density mismatch caused by the randomly-distributed material properties there are more significant. The corresponding L-and S-wave impedance reflectivity estimations are shown in figures 8(a) and (b). It can be seen that the coupling effects between velocity and density are compensated. It is noted that in this case, the homogeneous background is far from the true model, and thus small flaws may be not very accurately imaged. To improve the image quality, a low-frequency modelbuilding step can be introduced to reconstruct the large-scale features of the unknown background, e.g. elastic full waveform inversion (EFWI) [42], and then the high-frequency ELSRTM in the presence of density variations is expected to improve the reflectivity image of small flaws. As a comparison, the TFM image based on the same simulation is shown in figure 8(c). The delamination can be differentiated in the TFM image, however, this image contains strong artifacts which mask the indication of the small flaw.

Delamination in a layered structure
Layered or composite materials are commonly used in engineered structures. Many benefits are provided, however, the risk of encountering defects, e.g. delamination, is significant. A structure consisting of three layers and a 0.1 mm thin delamination with a length of 3 mm between layer 2 and layer 3 is shown in figure 9(a) (detail in figure 9(b)). The physical experiment for this structure is replaced by a simulated experiment using the elastic finite difference method. Absorbing regions are introduced to the left and right sides [43]. A linear transducer array with 120 equally spaced elements is placed on the top surface. The signal excited by the array is a 2-cycle Hann-windowed tone-burst with a central frequency of 2 MHz. The FMC data are collected. The homogeneous backgrounds of layer 1/layer 3 (V L = 5100 m s −1 , V S =3200 m s −1 and ρ=2600 kg m −3 ) are used as the forward models. Figures 9(c)-(e) show the estimated ELSRTM images of L-and S-wave velocities and density after 30 iterations. It can be seen that the interfaces show good agreements with the true models and the delamination is clearly visible in each parameter image. However, the reconstructed density image is slightly compromised at the recovered delamination. This is due to the fact that contrasts in density change the waveforms significantly [44]. In velocity images, the estimated S-wave velocity anomalies (figure 9(d)) are somewhat sharper than the L-wave velocity anomalies (figure 9(c))-a result of the fact that L-wave propagates faster and thus has longer wavelength signals compared with the S-wave. Summation of velocity and density images provides impedance reflectivity estimations, which are in better consistency with true ones than the velocity and density images separately, as shown in figures 10(a) and (b). In figure 10(c), the TFM result is less comprehensive and less accurate than the L-wave impedance reflectivity estimation ( figure 10(a)).

Experiment
In order to further confirm the imaging capability of the ELSRTM with density variations, a multi-material part fabricated by laser powder bed fusion is presented, as shown in figure 11. This multi-material part, which consists of Ti-6Al-4V and 316L stainless steel with the interlayer of K220 copper-alloy, was fabricated using the SLM250HL machine from SLM Solutions AG with dimensions of 24 mm × 18 mm×12 mm [45]. The designed flaw has the size of 2 mm×1 mm in the relatively brittle Ti-6Al-4V/ K220 interface, as shown in figure 11(a). The material properties are shown in table 1.
The experimental setup is shown in figure 11(b). The 64element ultrasonic array with a central frequency of 10 MHz (Guangzhou Doppler Electronic Technologies Co., Ltd) was  placed on the surface of the multi-material part. The array has an element pitch of 0.3 mm. The FMC data are collected by using an array controller (Advanced OEM Solutions, USA).
The resulting ELSRTM images generated on the basis of a forward model using a homogeneous background model of Ti-6Al-4V after 30 iterations are shown in figure 12. It can be seen from figures 12(a)-(c) that the position and the geometry of the interface flaw are clearly imaged in the L-wave velocity, S-wave velocity as well as density images, despite the porous nature of the printed multi-material part. There is a slight upward shift of this flaw due to the closeness in the time of arrival of the interface reflection and the flaw. Figures 13(a) and (b) show the L-wave and S-wave impedance reflectivity estimations obtained by summing velocity and density images in figure 12, which demonstrate better estimations and fewer artifacts compared with the L-wave and S-wave velocities. In the TFM image, as shown in figure 13(c), the result is less accurate than the L-wave impedance reflectivity estimation in figure 13(a).

Discussion
The results presented in this paper have demonstrated the effectiveness of the ELSRTM with density variations, which is a useful tool to image small defects in complex heterogeneous structures. However, it should be noted that inaccuracies in the background models of wave velocities and density may lead to an incorrect image of small defects. To improve the detection of delaminations, cracks and other small defects in engineered structures, an alternative hybrid approach that combines EFWI and ELSRTM with density variations can be used. First, a low-frequency EFWI method  is used to reconstruct the large-scale background models of the engineered structures, and second, high-frequency ELSRTM with density variations aims at imaging small defects based on the improved velocity models from EFWI. It needs to be mentioned that a homogeneous background is usually used in the density model since it is difficult to obtain an accurate density model in reality. As a partial remedy, the initial density model can be obtained from velocity models based on empirical relationships between density and velocity [46]. However, different empirical equations can provide different density estimations. Therefore, the sensitivity of the reconstructed images with respect to different density models with different degrees of accuracy must be investigated.
In the multi-parameter inversion, although the CG method used in this work is robust, the convergence in general is slow. For computational efficiency, the quasi-Newton L-BFGS and truncated Newton methods can provide faster convergence rates [39,47], which account for Hessian without explicitly building this matrix in the multi-parameter inversion. These optimization methods combined with the proposed hybrid approach will be a subject for future work.

Conclusions
The conventional ELSRTM, which ignores the variable density, may cause inaccurate amplitudes and crosstalk artifacts in L-and S-wave velocity images. In this paper, an ultrasonic imaging algorithm based on the ELSRTM with density variations, which can simultaneously reconstruct L-and S-wave velocities as well as density, is investigated for imaging small flaws in complex heterogeneous structures. In order to obtain improved reflectivity estimations, the L-and S-wave impedance images by summing the reconstructed velocity and density images are provided in this paper. The ELSRTM with the consideration of variable density has been applied to virtual and physical laboratory experiments and this approach has promises for characterizing small flaws in complex heterogeneous structures compared with other array image method. The further study will combine the lowfrequency EFWI and the high-frequency ELSRTM with density variations to improve the detection of small flaws.