Warped Wigner-Hough Transform for Defect Reflection Enhancement in Ultrasonic Guided Wave Monitoring

To improve the defect detectability of Lamb wave inspection systems, the application of nonlinear signal processing was investigated. The approach is based on a Warped Frequency Transform WFT to compensate the dispersive behavior of ultrasonic guided waves, followed by a Wigner-Ville time-frequency analysis and the Hough Transform to further improve localization accuracy. As a result, an automatic detection procedure to locate defect-induced reflections was demonstrated and successfully tested by analyzing numerically simulated Lamb waves propagating in an aluminum plate. The proposed method is suitable for defect detection and can be easily implemented for real-world structural health monitoring applications.


Introduction
In recent years, ultrasonic guided waves GWs have received a great deal of attention among nondestructive tests community due mainly to the ability to travel long distances without substantial attenuation and to employ multimode/-frequency examination for defect classification and sizing.Among the various techniques based on GWs, the detection of defects in plates-like structures by means of Lamb waves has been, and is still, widely investigated 1-6 due to the variety of potential applications.Since the propagation characteristics are directly related to both the inherent structure and mechanical properties of the medium, the dispersiveness of GWs can reveal important information for structural health monitoring purposes.Unfortunately, several different modes appear simultaneously in the signal.These modes overlap in both time and frequency domains, and simple Fourier analysis techniques are not able to separate them.
Thus, identification of Lamb modes is a challenging step in the process of damage detection.However, recent works in the area of time-frequency representations TFRs 7, 8 show great potential for applications in nondestructive evaluation and material characterization.Specifically, TFRs can provide an effective tool for the interpretation of GWs, whose multimodal and dispersive nature make them an extremely complicated class of ultrasonic signals.
This work proposes a time-frequency TF energy density function approach that makes use of known dispersion characteristics for a propagating wave mode in order to compensate for the effect of dispersion and locate defects in plate-like structures.Our approach will be illustrated through a relevant case study, in which defects are to be located on an aluminum plate where Lamb waves are excited.

Numerical Simulation of Lamb Wave Propagation
Let us consider an aluminum plate of thickness is h 2.54 mm, Young's modulus E 69 GPa, and Poisson's ratio ν 0.33.The proposed processing requires the computation of the group velocity dispersion curves for the plate.For such task, today several formulations and tools are available.For instance, in uniform waveguides the group velocity c g f can be estimated by means of analytical-based formulations 9 , semianalytical finite element SAFE simulations 10 , and by using standard finite element codes 11 .Recent developments allow the computation of the dispersion curves also in the case of irregular waveguides 12 .The results shown in Figure 1 were obtained by a using free-SAFE-based tool that can be downloaded at http://www.guiguw.com/ 10 .As a second step, time waveforms related to Lamb waves propagating in the aluminum plate were obtained numerically by means of dedicated Finite Element FEM simulations using Abaqus explicit 13 .Thanks to the Lamb problem symmetry, a x-y plane strain condition was assumed, as shown in Figure 2. A notch of width b 0.25 mm and depth a was considered, such that a/h 0.3.The assumed notch location was on the top side at the center of the plate x 500 mm .Lamb waves were excited by applying an impulsive force p t to the left edge of the plate towards the positive x-direction: this mainly stimulates the symmetric S 0 mode.The force was shaped in time as a triangular window with a total duration of 2 μs see the top left of Figure 2 in order to excite consistent Lamb waves up to 500 kHz.
For instance, a similar excitation was obtained in 14 by focusing, through conventional optics, a laser beam to a straight line acting along the plate edge.It is shown in such work that the S 0 mode can be excited by means of such experimental setup.
To ensure accuracy to the time-transient finite element simulations 15 , the plate domain was discretized with elements of maximum side length L max 0.125 mm, and the time integration step was kept Δt < 1e − 8 sec.
Time-dependent out-of-plane displacements v t were recorded at three points on the top side of the plate y h/2 , namely, A, B, and C, respectively, located at x A 100 mm, x B 200 mm, and x C 300 mm.The recorded waveforms are shown in Figure 3.The leftmost peak in each signal corresponds to the passage of the excited S 0 mode through the recording position path 1, in Figure 4 , while oscillations in the central part of the waveforms are due to defect-induced reflections path 2, in Figure 4 , which also excite the slower A 0 mode.Spreading of these oscillations clearly reveals the effect of dispersion.Finally, further reflections from the plate edges path 3, in Figure 4 are responsible for the complicated behaviour observed in the rightmost part of the signals.

Mathematical Tools
The defect location procedure can be divided in three steps: i Warped Frequency Transform to remove the dispersive behaviour of the S 0 mode; ii equalization procedure to enhance weak reflections; iii Wigner-Hough Transform to distinguish S 0 reflections from other interfering waves.Such steps are detailed in the following subsections.

The Warped Frequency Transform (WFT)
The WFT is a unitary time-frequency transformation that produces a flexible sampling of the time-frequency domain 16 .
Given a generic signal s t whose Fourier Transform FT is S f , the frequency warping operator W w transforms the original waveform into a warped version s w t by reshaping the frequency domain through a properly designed warping map w f .The procedure is depicted in Figure 5.
The WFT can be used to compensate dispersion in GWs 17 by defining the warping map through the following relationship: Equation 3.1 relates the functional inverse w −1 f of the map to the inverse of the group velocity curve, that is, 1/c g f , of the wave that we want to consider.K is a normalization parameter selected so that w −1 0.5 w 0.5 0.5.As shown in Figure 1, for acoustic emission below 500 kHz, only the two fundamental waves A 0 and S 0 can propagate through the plate.It is assumed that mostly S 0 is actuated, therefore the group velocity curve of S 0 has been used to build w f .If s t is an undamped guided wave at a traveled distance D from the actuator, its FT is given by S f S 0 f • e −j2π f 0 τ α dα where S 0 f is the FT of the actuated wave incipient pulse and τ f D/c g f represents the frequency-dependent group delay of S 0 in this case.Using 3.1 , S f can be rewritten as

3.2
By applying warping and exploiting the invertibility property of the map, that is, w −1 w f f, yields to a signal s w t whose spectrum is

3.3
The linear dependence on the warped frequency in this equation shows that the dispersive effect is converted in a simple warped time-delay KD proportional to the distance.Dispersion is therefore compensated, and the resulting signal s w can be equivalently plotted Mathematical Problems in Engineering as a function of the distance from the source, thus allowing to locate defects of the excited waveguide by detecting the corresponding reflected waves.
It is worth noticing that, despite the different formalism, the processing described above is substantially analogous to the ones presented in 18, 19 .However, in many practical applications such dispersion compensation is not sufficient to ensure a reliable estimation of wave traveled distance, due to the weakness of reflections and to the interfering presence of different modes.For these reasons, necessary further processing steps are introduced in this work, as will be shown in the following sections.

Wave Equalization
By means of the warping procedure described in the previous section, a realignment of the time-frequency content of S 0 waves in vertical lines is produced.However, in general the energy of the waves scattered by defects is much lower than the energy of the incident wave, especially for small defects.To overcome this problem, the energy of incident and reflected waves in the acquired signal can be conveniently equalized.Such task can be accomplished with a simple but effective procedure based on a local averaging of the acquired signal.Indicating with LA |s w x | the local average of |s w x | in the neighborhood of x, the equalized signal s we x is obtained as where T is a given threshold, set as the 5% of the maximum value of |s w x |, which is used to avoid the undesired amplification of numerical noise when the signal is absent.In experimental data, the value of T must be set according to the SNR of the acquisition setup.
The warped and equalized versions of the signal in Figure 3

The Wigner-Hough Transform (WHT)
After equalization, defect detection can be performed automatically with a further processing of the signal s we x .It is worth noticing that the energy of spurious contribution caused by multimodal propagation or mode conversion is quite high in the equalized signal compared to the energy of compensated S 0 waves.In the example of Figure 6 c , the peak related to the S 0 wave reflected by the defect indicated as S 0 path 2 is just about twice as high as the peak in the following mode-converted wave A 0 .Therefore it is quite difficult to implement simple thresholding procedures capable of distinguishing different wave modes in this representation.
However, defect-induced reflections of the analyzed mode S 0 in this example appear in the compensated waveform as well-localized spikes, thus producing vertical maxima lines in a TF representation, whose warped time location can be directly converted to the defect position.On the contrary, spurious contributions related to different modes A 0 in this example show a peculiar frequency modulation due to a different group velocity curve from the one of S 0 .This can be clearly observed in a simple TFR of s we t , provided by the shorttime Fourier transform STFT , shown in Figure 7.Other energy distributions, such as the Wigner-Ville distribution WVD , defined as can be used to further improve the effectiveness of the representation.In fact, WVD provides optimal energy localization of linear chirp signals in the time-frequency plane 20 .It follows that applying vertical line detection algorithms to the WVD of the compensated signals provides an asymptotically optimal detector of wave propagating distances.
In particular, automatical detection of the desired lines of energy maxima can be performed by applying the Hough Transform HT 21 to the WVD, resulting in the so-called Wigner-Hough Transform WHT 22 .Generally speaking, the HT is an image processing tool that performs an integration on all the possible lines of a given image I and maps the value of each integral to a ρ, θ plane corresponding to the polar parametrization of  the lines.High-intensity pixels concentrated on straight lines on I will therefore produce peaks in the ρ, θ domain.In the Wigner-Hough Transform, the input image corresponds to the WVD of the considered signal and, in our approach, emphasis is placed in finding vertical lines, located at θ {π/2, 3π/2}.Therefore, in the WHT, the portion which corresponds to these angles is isolated, and the ρ value corresponding to detected peaks represents the distance traveled by the wave.One of the major limitations of the WVD is the presence of interference terms between different spectral components of the analyzed waveform induced by the WVD.However, it is worth noticing that this inconvenience is largely compensated through the integration performed by the Hough operator, as these undesired components appear as alternating positive peaks and negative valleys in the TF plane.

Procedure of the Method
The mathematical tools detailed in the previous section can be efficiently implemented by the following processing steps.
i The discrete WFT can be computed with the approach described in 23 .In essence, the warped signal is obtained by performing a nonuniform Fourier Transform 24 followed by an inverse Fourier Transform.Fast-Fourier algorithms can be exploited to compute both the direct nonuniform and the inverse transforms.
ii In the second step, warped signals are equalized in amplitude, according to 3.4 .
The local averaging window applied in the right hand term of the same equation is about 9 cm in length it is worth recalling that, in the warped domain, time intervals are related to distances .iii The software Time-Frequency Toolbox-TFTB which computes the Wigner-Ville distribution adopted in this study is available for academic use at http://tftb.nongnu.org/.In the same tool, also the code for Hough Transform computation is provided.However, for wave propagating distance estimation, the calculation of the Hough Transform coefficients in the whole ρ, θ plane is redundant, as described in Section 3.3.For this reason, in our approach such calculation is simplified with a simple integration of WVD coefficients across frequencies.

Numerical Results
The WVD of the equalized warped signal in Figure 6 is depicted in Figure 8.As it can be seen, two vertical lines appears in correspondence to the actual traveled distances of the incident  and scattered wave at 0.1 and 0.9 m, resp.together with spurious contribution due to different modes of propagation and interference terms.
The computation of the WHT, that is, the Hough Transform on the WVD, produces the image depicted in Figure 9, where three peaks can be associated with θ π/2 and 3π/2.The peaks ordinates represent the difference between the effective distance of propagation of a given wave and a reference distance of propagation of 0.7 m, which corresponds to the half of the maximum considered propagation distance in the analyses of 1.4 m.
By extracting the values related to θ π/2 and 3π/2 in the WHT and reordering them according to the distance from the origin, a novel signal s we−WHT is obtained.In such signal, depicted in Figure 10, the peaks related to the actual traveled distances of S 0 waves clearly emerge with respect to the spurious contribution due to mode-converted waves, thus greatly simplifying the definition of automated distance estimation procedures.In particular, the amplitude of the peak related to the reflected S 0 path 2 is about 10 times higher than   comparable to that of the reflected spurious A 0 wave.It follows that the location of defects by exploiting only warped signals could be difficult in real noisy applications.As for the first case discussed, the proposed robust and automated approach involves first an equalization of the warped signal to emphasize the amplitude of reflection-induced peaks and next the detection of vertical maxima lines in a space-frequency representation WHT of the equalized warped signals.
The Wigner-Hough Transform, in fact, appears as a suitable tool to isolate S 0 components and locate defect-induced reflections, as it can be seen form Figures 13 and 14, where the WHT of the equalized warped signals v Bwe and v Cwe , respectively, are displayed.
Local maxima at θ π/2 and θ 3π/2 can be easily detected, and the corresponding ρ coordinates provide the distance traveled by the incident and reflected S 0 waveform components, respectively.
From the extraction of peak coordinates in the waveforms s we−WHT depicted in Figures 10, 15, and 16, the defect responsible for reflections is located at x 503 mm, x 502 mm, and x 502 mm, respectively.Errors with respect to the actual defect position x 500 mm are thus below 3 mm, which roughly corresponds to the minimum wavelength associated with Mathematical Problems in Engineering the excited Lamb waves.Similar good results were found by considering different defect depths and positions, not shown here for the sake of brevity.

Conclusions
The work described the application of a Warped Wigner-Ville analysis to improve defect detectability of conventional Lamb wave inspection systems.The proposed equalization approach effectively enhances the amplitude of relevant peaks in warped signals, where dispersion for a GW mode of interest has been removed.This procedure may encounter limitations in the presence of especially noisy signals, as spurious components might be erroneously amplified.However, several alternatives are possible, including more sophisticated preprocessing algorithms under investigation and the averaging of multiple acquisitions.
The Wigner-Ville distribution of the equalized signal is then computed.The presence of interference terms is largely compensated through the integration of the time-frequency decomposition performed by the Hough operator.
In the resulting Winger-Hough Transform representation, vertical lines associated with relevant acoustic events can be detected.This allows for the separation of overlapping Lamb waves.In particular, it was shown with numerical examples that the contribution of S 0 can be highlighted and the one due to interfering terms such as A 0 wave deeply attenuated.
Finally, with simple thresholding procedures, the information about the distance traveled by the incident and reflected components of a monitored wave can be easily recognized.Thanks to its very high precision the developed tool could pave the way for a new class of procedures to locate defects in waveguides.

6 Figure 1 :
Figure 1: Group velocity dispersion curves for the Lamb waves propagating in a 2.54 mm thick aluminum plate Young modulus E 69 GPa, Poisson's coefficient ν 0.33, and density ρ 2700 kg/m 3 .

Figure 2 :
Figure 2: Schematic representation of the damaged aluminum plate used in the time-transient FEM simulations plate dimensions are in mm .The spatial distribution and time-amplitude shape of the actuation pulse is also shown.

Figure 3 :
Figure 3: Time-dependent out-of-plane displacements recorded on the top side y h/2 of the plate at a x A 100 mm, b x B 200 mm, and c x C 300 mm.

Figure 4 :
Figure 4: Schematic representation of the multiple paths traveled by the waves detected in a given acquisition point.

Figure 5 :
Figure 5: Computational flow of the frequency warping operator W w .F and F −1 are the direct and inverse Fourier Transform operators, respectively, while w f and ẇ f are the warping map and its first derivative.
a are plotted in Figures 6 a and 6 c , respectively.The equalization factor max{LA |s w x | , T } is depicted in Figure 6 b .

Figure 6 :
Figure 6: a Warped version v Aw of the signal v A represented in Figure 3 a .b Equalization factor for the same signal.c Equalized warped signal v Awe as a result of the equalization procedure.

S 0 path 3 S 0 path 2 SFigure 7 :
Figure 7: Spectrogram of the equalized waveform v Awe represented in Figure 6 c .

Figure 8 :
Figure 8: Wigner-Ville distribution WVD of the equalized waveform v Awe represented in Figure 6 c .

Figure 9 :Figure 10 :
Figure 9: Hough Transform of the Wigner-Ville distribution WHT represented in Figure 8. Peaks at θ {π/2, 3π/2} dashed lines correspond to vertical lines in Figure 7; ρ coordinates provide the distance from the center of the analyzed image.

Figure 11 :
Figure 11: a Warped version v Bw of the signal v B presented in Figure 3 b .b Equalization factor for the same signal.c Equalized warped signal v Bwe as a result of the equalization procedure.

Figure 12 :
Figure 12: a Warped version v Cw of the signal v C presented in Figure 3 c .b Equalization factor for the same signal.c Equalized warped signal v Cwe as a result of the equalization procedure.

Figure 13 :
Figure 13: Wigner-Hough Transform WHT of the equalized warped signal v Bwe represented in Figure 11 c .Peaks at θ {π/2, 3π/2} dashed lines correspond to vertical lines in the Wigner-Ville distribution; ρ coordinates provide the distance from the center of the analyzed domain.

Figure 14 :
Figure 14: Wigner-Hough Transform of the equalized warped signal v Cwe represented in Figure 12 c .Peaks at θ {π/2, 3π/2} dashed lines correspond to vertical lines in the Wigner Ville distribution; ρ coordinates provide the distance from the center of the analyzed domain.