Matrix techniques for Lamb-wave damage imaging in metal plates

The implementation of efficient maintenance strategies of thin-walled structural components require reliable damage detection and localization techniques. In particular, guided ultrasonic waves technology represent an auspicious approach when implemented in a structural health monitoring system. The method is usually based on distributed sensing with piezoelectric elements that act in turn as ultrasound transmitter and receiver. This work aims at a unifying framework for damage localization considering algorithms from different scientific disciplines, e.g. originated from radar and geophysics. Here, we systematically express those algorithms in matrix form and compare the respective damage localization performance with experimental measurements considering an isotropic specimen with a single and also multiple simultaneous defects. In addition, we evaluate the algorithms’ point spread function and propose performance metrics to quantitatively compare the imaging success.


Introduction
In the field of guided-wave Structural-Health-Monitoring (SHM) the localization and classification of structural damage is generally the main goal. Ultrasonic transducers arranged as a spatially distributed array are being used [1][2][3]. Of particular interest for targeted maintenance activities is the damage position extracted inversely from ultrasound signals. Therefore a 3 This author contributed to equally Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. wide range of model-based algorithms [4][5][6][7], like Delay-and-Sum beamforming (DAS), synthetic focusing aperture technique, or Migration operators from geophysics can be used. Despite their different scientific heritage and the variety in implementations, those algorithms have the same fundamental principles. Firstly, they perform a virtual back-propagation of the measured wave field to its point of origin in time and space. Secondly, by applying a suitable imaging condition the back-propagated scattered wave field is simplified for the end user in terms of an interpretable information like reflectivity or energy spanning the spatial domain. In figure 1 an overview which unifies several approaches is presented.
For the purpose of modelling, a suitable linear wave equation has to be chosen, which includes most observable effects that are present in the data and which is able to account for the geometry of the investigated structure as well as the potentially heterogeneous spatial distribution of elastic parameters and therefore the distribution of velocities. The back-propagation is realised by connecting the fundamental solutions G (Green's function) of the modelled problem with the measured scattered wave field. Precisely, for calculations in the frequency domain the conjugated Green's function G * is used, whereas in time domain simply inverting the data is adequate. This virtual back propagation is possible due to the symmetry of any linear wave equation regarding the time axis. However, using a signal from a single observation point generally will not result in a sufficient reconstruction of the complete former scattered wave field. Instead the data of an adequate amount of observation points has to be used. This whole process is known as Time reversal mirror in the field of Non-Destructive Testing [8] and is the basic principle behind any geophysical migration operator. In theory, there is a discrimination whether the observation points are treated as boundary conditions on the surface or as arbitrarily distributed points inside the observation domain. In this work the latter one is used, whereas a comparison between both approaches for an analogous geophysical setting can be found in [9].
If an analytic version of Green's function is unknown, for example in the case of complex geometries or non-trivial spatial variation of the elastic parameters, the wave propagation can be solved numerically, which is known under the term reverse-time-migration (RTM) [10,11]. For the underlying mathematical principles the chosen numerical method is irrelevant, but in terms of numerical accuracy, effort and parallelization capabilities it has significant influence on the process. [7] applied RTM to a SHM problem, modelling an elastic homogeneous and isotropic medium, whereas [12] generalized it for elastic anisotropic models. Both authors used a wave equation based on the Midlin plate theory, modelling the propagation of flexural waves and therefore using the A 0 Lamb mode for the damage visualization. [13] reformulated the RTM as a least-squares problem for reflectivity mapping in geophysics, while [14] used this in the context of guided-wave SHM.
Instead of a complete numerical solution of the underlying wave equation, Green's function can be approximated by a ray-theoretical problem or two hyperbolic partial differential equations (eikonal-and transport equation) [15]. Nevertheless, the approximation will reduce the propagator implicitly to a certain phase like first or direct arrivals without multiple reflections. This circumstance can be avoided by explicitly modelling certain phases and add them to the propagator [16].
By applying an imaging condition the multi-dimensional result in space-time or space-frequency is either reduced or transformed into an interpretable quantity. Autocorrelation can be used to calculate the energy of the scattered wave field over the whole observation time [17,18], whereas the cross-correlation of different wave types indicate converted energy [17]. Both correlation types are generally implemented without any time lag, assuming the velocity function and therefore Green's function is sufficiently accurate [19]. Taking into account additional lags in time and space can increase the accuracy or can be used for estimating model adaptations [20]. Cross-correlation of the scattered wave field with former incident wave field itself has a special significance, since it connects scattering with the principle of back propagation [21,22]. The imaging result can be interpreted as reflectivity or perturbation in the velocity model and is used nowadays in full waveform inversion routines.
In this work a unifying framework for damage imaging techniques formulated in matrix notation is presented. In particular, methodologies from different scientific fields have been implemented and tested with experimental measurements from an isotropic medium. The imaging results are based on a single damage but also multiple damages at the same time.Throughout this report, we here assume baseline reference measurements from the undamaged structure and furthermore process the differential signal, i.e. the residual between measurements from the damaged structure and the aforementioned reference state.
The remainder of this paper is organized in the following way: The mathematical background of the utilized imaging techniques is provided in section 2, where a DAS-type algorithm serves as the benchmark. In section 3 the experimental setup is described. In sections 4 and 5 results are presented and discussed, respectively. Conclusions are finally given in section 6.

Mathematical description of matrix imaging techniques
The isotropic linear guided wave propagation can be asymptotically described by the Helmholtz equation of the form: ( Here u(x, ω) is a simplified scalar wave field associated with a certain guided wave mode of interest, at the location x, for the circular frequency ω after the excitation of the source function f at the source location s. The simplification is necessary since conventional piezoelectric transducers only measure a scalar voltage instead of the three-component elastic wave field. As a result the following modelling excludes other wave modes, mode conversion and the treatment of polarization. Nevertheless, the frequency-dependent phase velocity c of the target wave mode implicitly includes the dispersion of the simplified amplitude, if the solution is calculated for all frequencies of interest. Additionally the background velocity model c 0 (x, ω) is disturbed by a perturbation α(x), so that: Due to the imposed linearity, the wave field can be split into an incidence u i and a scattered wave field u s : Accordingly, the incident wave field can be considered as the guided wave SHM baseline measurement, whereas the scattered wave field represents the difference between a damaged state and the baseline. In realistic applications with varying additional mechanical loads or environmental conditions this becomes more complicated [23,24], requiring a sufficient amount of observations and accurate compensation or selection strategies [25,26]. Using the assumptions equations (3) and (2) in equation (1) leads to the Lippmann-Schwinger equation: by integrating over the whole observation domain Ω, the frequency range of interest and using Green's function G [27,31]. The non-linear equation includes multiple scattering at the sensor location r and therefore all damage cases which could be asymptotically described by the chosen wave equation.
For the sake of simplicity, the disturbance α is imposed as considerably small compared to the background velocity, so that the interaction of the scattering field with the anomaly itself can be neglected. This assumption is known as Born approximation and simplifies above equation further to a linear integration: describing the following synthetic experiment: The source pulse propagates from the transmitter location s to any point in observation domain x. The resulting incidence wave field is interacting linearly with the anomaly distribution, causing a set of scattering signals, which are finally all independently propagated to the receiver location r and superimposed to the artificial measurement. In the following the equation will be referred as forward operator.
The characteristic of an operator becomes more clearly, if equation (5) is reformulated as a discrete linear equation of the form: for all K frequencies, R receivers and all N points of interest inside the domain: Since most transducers in SHM networks can usually act as source and receiver consecutively, a complete dataset consist of various measurements of the above form with S different sender positions x s . As a consequence the data space in equation (6) of size K × R can be increased to M = K × R × S by stacking the different datasets, without changing the model space of size N. In the following, this stacked artificial measurement setup is implicitly assumed. Conventional imaging of the velocity perturbation is realised by the adjoint of the forward operator [28]: which is the conjugated and transposed operator matrix F T multiplied with the data. For the chosen discretization and including all source points, this corresponds to socalled prestack migration [11], while in conventional Non-Destructive Testing this process is known as Total Focusing Method [29]. A physical interpretation of the adjoint can be derived by applying it on equation (5), the integral equation of the forward operator: Transformed into time domain; the expression would correspond to the zero-lag cross-correlation of the forward propagated second-order time derivative of the source pulse and the scattering field. The latter one is realized by introducing the foremost measured data at all receiver locations as virtual sources and by propagating it backwards in time. Above formulation is usually implemented without the computer memory intense assembling of the linear equation system itself.
The adjoint has the advantage of being a unique, stable and cost efficient operation, but it will not necessarily map the correct scattering intensities since it only re-adjusts phase or time shifts, respectively. Additionally it consists of noise, exhibits truncation due to limited aperture, is incomplete or suffers from aliasing due to coarse sensor spacings, finally leading to imaging artefacts [28,30]. This issue can be overcome by calculating the inverse of equation (6), which can also be interpreted as deconvolution of equation (5) [31] Since the operator matrix F is mostly not square nor invertible, the model vector m is calculated indirectly by solving an optimization problem for a given L p -norm. [30] used a preconditioned conjugated gradient least-squares scheme with p = 2 and a Kirchhoff-Migration operator, while [13] used the same method for a RTM operator. [32] compared the method of [30] with various optimization problems formulated in the L 1 -and L 0 -norm, giving compressed sensing solutions. Compressedsensing methods and sparse representation techniques [4,33] are promising approaches for enhancing imaging results.
In the work at hands, results for different operators and solvers are compared for a typical SHM experimental setup. Beside the adjoint or conventional solution, the least-squares method (LSM) of [30], the orthogonal matching pursuit (OMP) implemented by [34] and the SPGL1 [35] solver of [36] are being used. Precisely, the latter one corresponds to the basis pursuit denoising (BPDN) problem:   with the apriori regularization parameter τ . The OMP is implemented by solving the convex least-squares problem: ∥Fm − d∥ 2 2 subject to min ∥m∥ 0 ≤ E.
Here, the L 0 norm seeks a solution with a given number of nonzero entries E. In certain cases, this apriori enforced sparsity can be critical, because it limits the number of defects which can be resolved or it limits the spatial resolution of a formation of multiple defects.

Kirchhoff-migration
The Kirchhoff operator in equation (8) is realized by using an analytical Green's function: with the phase velocity c(ω) of a chosen wave mode and without acknowledging any boundary reflections. Therefore the operator only contains information of direct arrivals and treats other modes or multipath contributions as noise. The propagation is assumed to be two-dimensional, with an associated cylindrical geometrical spreading rate, because guided elastic waves in plate-like structures show similar characteristics. For the comparison with different methods, the imaging output I(x) is given by the squared absolute value of the calculated scattering intensities m or m l , respectively. Note that, especially in geophysical literature, it is common to map the real part of m. Also the used source wavelet is transformed into zero phase beforehand.

Reverse-time-migration
For reasons of better comparison with the other two operators and in contrast to classical numerical formulations, the RTM operator in this work is realized by a series of semi-analytical expressions. Each one acknowledging a specific side-wall reflection of certain order, minimizing numerical influences in the benchmark. For each source position s a reflection at each side wall, additional Green's functions are defined, with virtual source positions constructed by the method of images and an algebraic sign swap. By a subsequent recursion scheme multiple side-wall reflections can be added as well. In line with Kirchhoff-Migration, the imaging output is again given by the squared absolute value of the scattering intensities m or m l .

Delay-multiply-and-sum
Delay-and-Sum-type beamformers constitute an important class of techniques that have already been employed early [8] in the development of guided-wave damage imaging. An improved variant of Delay-and-Sum that is named Delay-Multiply-and-Sum (DMAS) emerged in the field of microwave-based breast-cancer detection [37] and has since been implemented for ground-penetrating radar [38], medical ultrasonics [39] as well as photoacoustic microscopy [40]. In the present article, DMAS is introduced for Lamb-wave damage imaging. It here serves as a state-of-the-art benchmark for the other presented techniques. The ability to inspect visually structural defects calls for an image formation process which provides a intensity map I(x) based on differential signals u s ij (t). In the DAS framework, the contribution originating from sender i and receiver j (out of the S senders and R receivers) reads: where temporal integration is carried out over a discrete window T and where furthermore D ij (x) indicates a round-trip delay. Here, D ij (x) = (|s i − x| + |r j − x|)/v is derived from the group velocity v of the chosen mode and excited centroid frequency, the distance |s i − x| from sender i to position x as well as the distance |r j − x| from this point to the receiver j. Moreover, pre-processed signalsũ ij (t) are considered: the differential signals u s ij (t) are detrended, smoothed (moving-average filter and zeroizing amplitudes below the noise threshold), also aligned in order to correct for phase jitter and moreover distance-gated with respect to the chosen mode. DMAS introduces usage of the products of the signals where all possible signal combinations are involved: The introduced multiplications are either interpreted as correlations [38] or geometric means [39,40]. Formally, ∑ T qũ ij · u mn represents the unnormalized Pearson correlation coefficient ofũ ij andũ mn which accordingly causes that coherent components in both signals contribute stronger to the image intensity I than incoherent signals. The productũ s ij ·ũ mn = (√ũ ij ·ũ mn ) 2 can also be interpreted as the squared geometric mean of both signals where squaring omits the problem associated with possibly negative signs of the products. Instead of S · R signals in the case of Delay-And-Sum, DMAS hence makes use of an increased number of 'mean' signals, namely (S · R) 2 . [41] suggested a non-linear aggregation of the aforementioned individual sender-receiver contributions, which is utilized also in the present work. Instead of linear superposition , a medianbased fusion is used I(x) = median ijmn (I ijmn (x)).  and a thickness of 1.5 mm was placed in a climatic chamber in which controlled temperature conditions are present. Twelve circular piezoelectric transducers with a diameter of 10 mm and a thickness of 0.2 mm were used. The exact positions of the transducers are listed in table 1. Each actuator-sensor pair is recorded at room temperature with a multiplexing unit described in [42]. The excitation signal in the experiment is a Hann-windowed tone-burst with a carrier frequency of 70, 100, 200, 300, 400 and 500 kHz. Damages in the form of through holes with a diameter of 5 mm were inserted at locations listed in table 2. This arrangement implies measurement data from a single damage as well as two and three simultaneous defects of the same kind at different locations, probing if the imaging algorithms are able to discriminate between those adjacent scatterers.

Benchmark parameters
Calculations were performed for all three damage configurations, the three operators, six excitation frequencies, up to four solving methods and various spatial sampling rates dx. The sampling of 10 mm are used for the sparse OMP and BPDN solvers, whereas 2.5 mm used in the other cases. Since for DMAS the additional pre-processing is essential, those calculations were carried out for the other two operators as well. For the OMP solver the maximum number of damages was limited E ≤ 5 and for the BPDN approach the regularization parameter was chosen as τ ≤ 2. The semi-analytical RTM operator included side wall reflections up to the second order. Finally all imaging results were normalized afterwards and the following benchmark parameters were extracted: peak-signal-tonoise-ratio (PSNR), maximum value of the background noise: , standard deviation σ(x b ) of the background noise and the maximum peak values: max ( ) of each single hole in a damage formation inside the associated so-called trust region. Therefore, a circular region with a diameter of 3.5 cm was defined around the center of each hole. All values inside the regions x d k were assigned to the respective hole k and all N b values outside were declared as background x b . A single hole k is detected, if its maximum peak value is higher than the background peak. And a whole damage formation is correctly detected ('✓'), if this is true for all holes. Otherwise the formation is only partially ('(✓)') or not detected ('') at all. For the normalized results the logarithmic scaled PSNR is defined as: with the perfect reconstruction K(x) for a given spatial sampling rate. Here, K(x) acts similar to an indicator function which equals 1 exactly at the damage position and is zero elsewhere. Note that the square-root of I(x) is used, since the expression was already squared by the associated imaging function.

Damage localization for single damage
All the results concerning the detection success are provided in table 3

Damage localization for multiple simultaneous defects
Considering here multi-damage scenarios, the following depictions are presented: In figure 7 the conventional/adjoint solutions at 200 kHz for all three operators are drawn. Keeping the same frequency, the results concerning the least-squares and the sparse solutions are given in figure 8. Additionally, in figure 9 the PSNR is depicted for all studied frequencies. Furthermore, detection quality among the methods is compared in figure 10.

Discussion
All tested methods, i.e. operator-solver combinations were in general able to detect damage, especially given a suitable inspections frequency selection (see table 3). From table 3 it can also be observed that the single defect could be detect with greater success, likely because a single defect approximately fulfills a Born-type assumption. In general, the methods were also able to clearly separate multiple defects at least at one of the studied central frequencies. For example, in the three-hole scenario RTM in conjunction with LSM was able to detect all defects only at 200 kHz, exhibiting the poorest performance in this regard among the tested methods and in terms of the studied metric. Form the extracted PSFs in figure 5 the increasing resolution by applying L p -Solvers can be obtained for all methods. The resolution of RTM and Kirchhoff are comparable, whereas DMAS showed significantly reduced oscillations epically along the vertical axis. This circumstance can also be obtained from figure 3. Remaining imaging artefacts around the damage in figures 3 and 4 are assumed mainly to cohere from the limited amount of sensors and modelling errors. Precisely an inaccurate velocity model and the simplification due to the chosen Helmholtz equation, implicitly neglecting polarization effects and mode conversion. With additional holes in figures 7 and 8 the noise level increases for all operators, most likely due to effects of multiple scattering, which are not included by the underlying Born assumption.
Comparing the metrics in figures 6, 9 and 10 RTM was generally outperformed by its Kirchhoff counterpart. This circumstance, at first, contradicts the intuition that additional information, i.e. acknowledging the further illumination from the virtual sources due to the side-wall reflections, should lead to better results. However, the noise in the data and the modelling error can outweigh this benefit. Firstly, reflections generally have longer travel paths and are therefore more sensitive to modelling errors by imperfectly assumed phase velocities or side-wall geometries. Secondly, the longer travel path inherently leads commonly to lower amplitudes due to damping, which in turn decreases the signal-to-noise ratio. Thirdly, especially for the considered sensor arrangement, there is no portion of the plate which is not significantly illuminated by first arrivals.
In contrast, the advantage of sparse solvers consists in a significantly improved signal-to-noise ratio. In figure 6 the OMP and BPDN based approaches consistently provide a higher PSNR than the other implemented approaches, particularly for 200 and 300 kHz. The fact that the performance varies with excitation frequency can rely on different factors, two important ones being the following reasons: (i) The applied operators can reproduce the experimental signals well, if the tuned excitation behaviour of Lamb waves [43] is captured well by the operator. For the given experimental setup, in the range of 200 and 300 kHz Lamb waves exhibit a pronounced S 0 mode which is reproduced successfully by all the studied operators. And (ii) the wavelength of the actuated waves changes with frequency. Assuming Lamb waves with pure S 0 mode, then the wavelength will drop with increasing frequency leading, at least in theory, to an enhanced resolution of smaller defects. The improved resolution of individual defects might be responsible for the good performance observed at 500 kHz in the multi-site damage configuration depicted in figure 9. Nonetheless, interpretation of the intricate scattering processes involving Lamb waves is complicated due to their complex propagation behaviour.
Among the sparse solutions, OMP outperforms BPDN in terms of PSNR in figures 6 and 9. DMAS shows overall similar performance like the sparse solvers, potentially in part due to the median-based filtering which facilitates clutter suppression. Finally, when looking again only at the most promising centroid frequency of 200 kHz, from figure 10 it is evident that the Kirchhoff operator in conjunction with the OMP solver delivers the minimal maximum-noise value (black lines) and the smallest variance in noise (red lines) which both together facilitate the clearest discrimination between defect and noisy background. Especially, for the three-hole scenario the reconstructed intensity of all three damages is comparatively high. The corresponding imaging result is given in figure 8 where Kirchhoff OMP hence provides the most promising support regarding the unambiguous identification of spots for targeted maintenance.

Conclusion
In this work, damage imaging using guided ultrasonic waves is studied in an experimental setting comprising an aluminum plate with single and multi-site damage in form of throughholes. Several model-based imaging algorithms are introduced in a unifying framework. Those methods are compared based on different figures-of-merit. The results are discussed in light of theory-based arguments. Due to the consistent formulation of the problem based on scattering theory, it is possible that the imaged intensity of the scatterer is correlated with e.g. damage size. Therefore, the presented methods could be utilized in the future in a gauged SHM system in order to estimate damage size or damage severity. Additionally the possibility to approximate Greens function could include more complex geometrical features like nozzles, varying thickness or bended parts. Moreover, through deployment of iterative procedures and by taking into account higher-order scattering processes, it is a formidable open question for future research, whether the resolution could be further improved in the case of multisite damage.