Abstract
Dispersive propagation and overlapping wave modes are two main obstacles for guided Lamb wave SHM applications. In an effort to overcome such obstacles, a new signalprocessing technique taking advantage order tracking based on dispersion relation, is developed. In this approach, by referencing the wave numberfrequency function of specified mode, the operations of resampling and interpolating are performed on the frequencyspectral series of raw signal. The orders referenced to wave numberfrequency are calculated, according to which the individual wavepacket is identified and its corresponding propagating distance is estimated. In the order domain, the overlapping modes are readily separated by Gabor expansion on the frequencyspectral series of raw signal. Numerical and FEM simulations on strongly dispersive and multimode overlapping guided waves were carried out to evaluate the performance of the proposed approach. The results demonstrated that the proposed approach is effective in dispersion analysis, mode differentiation and overlapped wavepackets separation.
1. Introduction
Guided ultrasonic wave, with stress distributed through the thickness of the material, can propagate along the thin and large components of onedimensional structures with little loss in energy. It offers the potential to achieve the required permanent structural integrity monitoring [1, 2]. It is difficult to control the excitation and reception of selected guided wave modes. The main challenge for processing Lamb wave signals is how to identify and separate multimode overlapping wavepackets due to their different dispersion characteristics and different propagation velocities [3, 4]. A variety of pragmatic signal processing approaches have been developed and contributed to ascertainment of multimode overlapping wavepackets. These conventional techniques mostly focused on the improvement of timefrequency or time resolution and the suppression of dispersion. In the conventional methods, narrowbandwidth and singlemodal guided waves with the selected central frequency were generally used as the excitation and ignored the dispersion effect.
The timefrequency representation (TRF) were adopted to separate the damagereflected and boundaryreflected wavepackets and identify their individual wave mode, such as shorttime Fourier transform (STFT) [5], wavelet transform [6, 7], S transform [8] etc. However, because of the continuously evolving shape (increased duration and reduced amplitude) and modes overlapping in both the frequency and time domains, many of TFR techniques, limited by their timefrequency resolution, either lose effectiveness or lack sufficient accuracy. In order to produce a much more distinct timefrequency resolution, the spectrogram reassignments [9, 10] were introduced to enhance readability of TF representations for multimode dispersive Lamb wave signals. But the reassignment methods lack the noise robustness and their resolutions may be degraded with the increase of wavepacket number.
To improve the time resolution of general Lamb wave signals, a number of deconvolution techniques were applied including Wiener pulse shaping filter, twosided Wiener filter, weighted least squares filter, Mendel’s minimum variance deconvolution algorithm, Oldenburg’s frequency domain deconvolution algorithm and L1norm deconvolution [11, 12]. Since the deconvolution methods were very sensitive to noises, they generally could be inapplicable if strong noise were present in the sensor signal. Furthermore, even with high time resolution, the method of pinpointing wavepacket arrival time maybe loses its effectiveness. That is because Lamb waves are inherently dispersive, i.e., different frequency components propagate at different speeds in a wavepacket, which cause wavepackets to spread out in space and time when they propagate through a structure so as to make the concept of arrival time illdefined in time domain. To obtain welldefined onsettime of echoes, Chirplets based dispersion curve were employed to extract the onsettime of that single mode from multimode dispersive wave signal [13, 14]. However, in the case of coherent interferences (false reflection) in the waveform, the Chirplet method may produce misleading results. To eliminate interference of false reflections, Ajay et al. [15] adopted an algorithm based on correlation analysis and matching pursuit method to estimate onsettime of echoes. But this method still has a deficiency in that the wave propagation model and baseline signal must be known.
In an effort to overcome the difficulty of wavepacket identification caused by strong dispersion, researchers developed some signalprocessing techniques to compensate for the dispersive nature of Lamb waves in longrange inspections. The time reversal technique and its improved methods [1619] are used to achieve temporal recompression of Lamb waves under broadband signal excitation. Warped frequency transform [20] upwarps the dispersive signal, and then wraps TF points of the STFT of the signal to obtain a dispersionmatched warpogram. The frequency domain polynomial Chirplet transform [21] adopts a rotation operator and as shift operator to do dispersion compensation.
Although the methods above are useful to identify the overlapping wavepackets by taking dispersion nature into account, they don’t provide the algorithms to extract the individual mode wavepacket for consequent damage evaluation. Whatever is taken to improve the resolution of diagnostic waves and whatever is taken to suppress dispersion effect, to separate completely overlapping multimode wavepackets is always restricted to a certain extent. For this reason, a new approach is developed to process strong dispersive and overlapping multimode wavepackets. In this proposed method, by taking the wave numberfrequency function of different mode as referencing phase, the corresponding mode wave can be identified and reconstructed with Gabor expansion. This article is organized as follows. In Section 2, the mathematical model of Lamb wave propagation is presented. Section 3 explains the basic theorem of order tracking technique in Lamb wave analysis. The detailed algorithm for order tracking technique is shown in Section 4. The proposed method’s applications in processing multimode overlapping signal and analysing dispersion signal are described in Section 5 and Section 6, separately. Section 7 outlines the possibilities to apply the presented method to the SHM application and describes briefly ongoing studies and the path for future studies.
2. Mathematical model for Lamb wave propagation
Because of the dispersion phenomenon in Lamb wave propagation, when Lamb waves travel across the structural features like boundaries, stringers, etc., reflections from these features produce phase shifts with respect to the original excitation signal. To take into account the dispersion effect of elastic waves, the dispersion physics of single mode wave propagation in 1D dispersive media can be modelled by Fourier analysis. When the incident wave $h\left(t\right)$ is to be generated by a transducer and propagates along a waveguide over the distance $x$, the arrival time of each frequency component will delay $x/{c}_{p}(=kx/2\pi f)$, where ${c}_{p}$ is the phase velocity. After $h\left(t\right)$ propagating over $x$ distance, the Fourier transform of $h\left(t\right)$ (noted as $H\left(f\right)$) will be shifted in the phase by ${e}^{jk\left(f\right)x}$. The Fourier transform of received wavepacket which can be phaseshifted wave, is given by:
where $k\left(f\right)$ denotes a wave numberfrequency function or dispersion relation function. The coefficient $C$ is nearly independent to $f$ after the wave propagates along or is reflected by either damages or boundaries. In practical applications, there inevitably exists a multimode wave in the received signal, since the current methods such as double side excitation only can be used to strengthen one mode as well as weaken the other mode, but to eliminate completely one mode is impossible. Therefore, the resulting sensor signal in frequency domain can be written as a superposition of multiple mode wavepackets from multiple propagation paths, i.e.:
where $N$ is the total number of propagating wave paths. ${x}_{n}$ and ${C}_{n}$ are the propagation distance and the transmission coefficient of $n$th wave path, respectively. ${k}_{m}\left(f\right)$ is the wave numberfrequency function of mode $m$.
Although the dispersion relation is an obstacle in Lamb wave applications, it can also serve as a useful tool to discriminate one mode from the other modes. $S\left(f\right)$ in Eq. (2) is presented as the linear superposition of components whose phases are the multiples of the wave numberfrequency function ${k}_{m}\left(f\right)$.This relation between $S\left(f\right)$ and ${k}_{m}\left(f\right)$ is very useful to identify ${k}_{m}$ ($f$)related wavepackets and to track their propagating distances ${x}_{n}$. If ‘order’ is defined as ${k}_{m}\left(f\right){x}_{n}$_{}normalized by ${k}_{m}\left(f\right)$, more accurate interpretation of received signal can be obtained by converting $S\left(f\right)$ into socalled order spectrum. It is obvious that there is onetoone relationship between the ‘order’ and the individual mode wavepacket. The order spectrum provides more conveniences for separate multimode overlapping wavepackets.
3. Principle of order analysis for Lamb wave
As what is mentioned above, if the order spectrum of some specific mode can be obtained, the wavepackets in the specific mode can be identified. It will improve the efficiency of propagating distance estimation, and simplify the computation process of pinpointing damage position. According to Eq. (1), the frequency domain signal $S\left(f\right)$ can be regarded as a frequencyvarying or phaseshifting waveform. The dispersion characteristic of sensor received signal depends on the dispersion curve function $k\left(f\right)$. If $k\left(f\right)$ is a linear function of $f$, $S\left(f\right)$ is the Fourier transform of nondispersed wave. If $k\left(f\right)$ is a quadratic or higher order function of $f$, $S\left(f\right)$ represents the Fourier transform of dispersed wave. In order to describe the relationship between $S\left(f\right)$ and $k\left(f\right)$, the dispersion signal $S\left(f\right)$ sampled constant in frequency discrete signal is converted into the equiangular signal sampled constant in angle from. From the mathematical point of view, this task could be solved by resampling and interpolation theory, which is explained in the example given in Fig. 1. The excited signal is assumed to be Gaussian modulated wavepacket presented as:
where${f}_{0}=$ 60 kHz, $\sigma =$ 2×10^{5}. Assuming that the wave numberfrequency function is $k\left(f\right)=2\pi f(3f+1)$ and that the dispersion wave after propagating along x distance presents as $S\left(f\right)=CH\left(f\right){e}^{jk\left(f\right)x}$ ($C=$1, $k\left(f\right)=2\pi f(3f+1)$, $x=$ 2 and $H\left(f\right)$ is Fourier transform of $h\left(t\right)$). The assumed dispersion curve and the dispersive wave are shown in Fig. 1.
Because the phase shift of $S\left(f\right)$ has strict linear relationship with $k\left(f\right)$, which provides the feasibility of calculating resamplings for $S\left(f\right)\text{.}$ Taking $k\left(f\right)$ as the reference phase i.e., referencing to the equal interval ordinate value in Fig. 1(a), the assumed dispersed signal in Fig. 1(b) is resampled at constant increments of the phase angle. The dispersed signal is transformed from sequence $S\left(f\right)$, which is at constant frequency increments $\left(\mathrm{\Delta}f\right)$ in frequency domain into sequence $S\left(\theta \right)$, which is at constant angle increments ($\mathrm{\Delta}\theta )$ in angledomain. And then, the equiangular signal is interpolated. The interpolated signal is presented as the stationary wave in Fig. 2(a). The order spectrum (shown in Fig. 2(b)) is obtained by the FFT on the stationary wave. The order spectrum in Fig. 2(b) has a single frequency component at 2order, which is independent of the wave numberfrequency function $k\left(f\right)$. Since order analysis relates $S\left(f\right)$ to the reference wave numberfrequency function $k\left(f\right)$, the order is defined as the phase $xk\left(f\right)$ normalized by $k\left(f\right)$. As long as the order is obtained, the dispersion wavepacket is identified and its propagating distance $x$ is determined.
Fig. 1Scheme of resampling: a) the resamplings determined from the assumed dispersion curve function kf, b) the assumed dispersion wave in frequency domain Sf sampled using constant increments of angle
For multimode Lamb wave, there are multiple modes propagating simultaneously in the structure, and each is independent with distinctly different dispersion characteristic presented as different wave numberfrequency function. It can clearly be seen that the order spectrum only displays the orders related to the reference $k\left(f\right)$, the wave modes that have different dispersion properties will be excluded. Namely, for the multimode dispersion signal, just taking the numberfrequency function of specified mode as the reference phase, the specified mode wavepacket can be identified in the obtained order spectrum and its propagating distance can be estimated without its flight time. The wavepackets related to other modes would appear as uncorrelated noise. That is, the order tracking technique can be used to identify the mode of Lamb wave.
4. Algorithm of order tracking for Lamb wave
On the basis of the theorem aforementioned, the order tracking algorithm for Lamb wave analysis is proposed and can be summarized as the flowchart in Fig. 3.
Fig. 2The order analysis for the simulated signal: a) the equalized angle simulated reflection signal in frequency domain, b) order spectrum of reflection signal
a)
b)
Fig. 3Flowchart of order tracking algorithm for Lamb wave
4.1. Determination of resampling points
It has been clarified that the phase of $S\left(f\right)$ has the linear relationship with the wave numberfrequency function $k\left(f\right)$. In this case, the order tracking technique is converted to process $S\left(f\right)$ by taking ${\varphi}_{0}\left(f\right)=lk\left(f\right)$ (Generally, $l$ is a constant not equal to 1 so as to avoid too high orders.) as the referenced phase. Although the precise $k\left(f\right)$ can be approximated by highorder polynomials, generally cubic polynomial can describe dominant dispersion phenomena fair enough. In this case, the referenced phase presents as:
where $a$, $b$, $c$, $d$ are the cubic function fitting parameters for $k\left(f\right)$, $m=$ 1, 2,… is the index of resampling series, and $\mathrm{\Delta}\theta =2\pi /orde{r}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ ($orde{r}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the max order range). According to Eq. (4), the equiangular sampling frequencies can be determined from ${\varphi}_{0}\left(f\right)$. Because the referenced phase ${\varphi}_{0}\left(f\right)$ can be known from $k\left(f\right)$ in advance, it is possible to calculate the corresponding frequencies for each constant angular increment, $m\mathrm{\Delta}\theta $, and to solve Eq. (4) for given $f$as well. Setting$e=dm\mathrm{\Delta}\theta $, and $f=yb/3a$, then Eq. (4) becomes:
where $p=(3ac{b}^{2})/3{a}^{2}$, $q=(2{b}^{3}9abc+27{a}^{2}e)/27{a}^{3}$. The resolution for $y$ is:
${y}_{3}={r}_{2}\sqrt[3]{\frac{q}{2}+\sqrt{\mathrm{\Delta}}}+{r}_{1}\sqrt[3]{\frac{q}{2}\sqrt{\mathrm{\Delta}}},$
where ${r}_{\mathrm{1,2}}=\frac{1}{2}\pm \frac{\sqrt{3}}{2}i$, $\mathrm{\Delta}=\frac{{q}^{2}}{4}+\frac{{p}^{3}}{27}$.
If $\mathrm{\Delta}>0$, the Eq. (5) has only one solution i.e. ${y}_{1}$.
If $\mathrm{\Delta}=0$, the Eq. (5) has three same solutions i.e. ${y}_{1}={y}_{2}={y}_{3}=\sqrt[3]{q/2}$.
If $\mathrm{\Delta}<0$, the Eq. (5) has three different solutions, i.e.
${y}_{\mathrm{1,2},3}=2\times \sqrt[3]{\alpha}\mathrm{c}\mathrm{o}\mathrm{s}\left(\right(\phi +2k\pi )/3)$, $k=$ 0, 1, 2, where $\alpha =\sqrt{{p}^{3}/27}$, $\mathrm{c}\mathrm{o}\mathrm{s}\phi =q/2r$, $\mathrm{s}\mathrm{i}\mathrm{n}\phi =\sqrt{\mathrm{\Delta}}/r$. The effective solution in ${y}_{\mathrm{1,2},3}$ is the first solution whose value is bigger than the last sampling time. So, the resampling frequency corresponding to the given equiangular increment can be determined as${f}_{re}=yb/3a$.
4.2. Interpolation and determination of propagation distance
Once the resampling frequencies ${f}_{re}$ for equiangular increment are determined, the corresponding frequency amplitudes can be calculated by interpolating $S\left(f\right)$. The accuracy of the interpolation method determines the accuracy of the resampling amplitude. There are diverse interpolation methods that can be used, such as linear and cubic spline interpolations. In this paper, cubic spline interpolation is selected, which is to use a cubic curve to fit the four points (two raw frequency points before and two frequency points after the interpolation point). The formulation of the interpolation polynomial is given:
Firstly, to calculate the coefficients $\left\{{c}_{0},{c}_{1},{c}_{2},{c}_{3}\right\}$ which are used to determine the interpolated $S\left(f\right)$ at any given $f$, four sets of data points (i.e. pairs of ${f}_{i1,i,i+1,i+2}$ and ${S}_{i1,i,i+1,i+2}$) are substituted into Eq. (7) to generate four independent equations. These equations are solved to get the coefficients $\left\{c\right\}$, which are then used to determine $S\left(f\right)$ (the signal amplitude) at any given $f$ (the resample frequency). Once the amplitude at the any frequency is estimated, the equiangular interval signal would be obtained. Finally, Fourier transform is performed on the resampled $S\left(f\right)$ to get its order spectrum. According to the order tracking theorem in Section 2, the “order” in the order spectrum of $S\left(f\right)$ is defined as its phase normalized by the referenced phase ${\varphi}_{0}\left(f\right)$. The same order is defined as the multiple of ${\varphi}_{0}\left(f\right)$, i.e. the propagating distance $x=order/l$. By taking the wave number function of specified mode as the reference phase, the order spectrum corresponding to the specified mode is obtained. The propagating distance of individual wavepacket is onetoone corresponding to the order in the obtained order spectrum. According to the propagating distances, the individual wavepacket from different propagating path can be identified.
4.3. Extraction of order component
It has been mentioned that the order spectrum only displays the orders related to the reference wave number function and that the wave mode having different dispersion properties would be excluded. Moreover, the wavepackets in the reference mode can be identified one by one in the order spectrum and their propagating distances can be estimated according to the onetoone relationship between the wavepacket propagation distance and its corresponding order. In some applications, however, engineers do need some wavepacket waveform of particular order to quantify the damage degree. Although shortterm Fourier transform characterizes each individual transient and nonstationary harmonic (or order) in its spectrogram distribution, it is invalid to reconstruct target order/spectral components embedded in the analyzed signal. The Gabor order tracking technique is proposed with extremely good selectivity on the extracted wavepacket. The process to recover the wavepacket of specified order can be done with the Gabor transform and its inversethe Gabor expansion. The Gabor expansion coefficients of $S\left(f\right)$ are evaluated through the Gabor transform. The Gabor expansion of $S\left(f\right)$ can be expressed as:
where $i$, $j$, $I$, $J\in Z$, ${c}_{i,j}$ denote the Gabor expansion coefficient and ${g}_{i,j}\left(f\right)$ is the shift and modulated form of the Gabor elementary function. $I$ denotes the time sampling number, $J$ denotes the frequency sampling number or frequency bins. Gabor expansion coefficients, ${c}_{i,j}$ are expressed be using the Gabor transform, i.e.:
where ${g}_{i,j}^{*}\left(f\right)$ is the corresponding dual function of ${g}_{i,j}\left(f\right)$. The Gabor expansion coefficients associated with the Gabor transform are flexibly applied to reconstruct the order components of interest. The computation employs masking strategies. The order components of interest that are covered by designated masks can is singled out and reconstructed from the multicomponent reflection. The Gabor spectrogram of $S\left(f\right)$ is generated by mapping the Gabor coefficient array into a 2D image of Gabor spectrogram distribution. The Gabor spectrogram can characterize the energyintensive components of $S\left(f\right)$ in certain frequencies that are relevant to the referenced $lk\left(f\right)$. A mask operation is performed on the initial Gabor coefficient array ${c}_{i,j}$ using the referenced $lk\left(f\right)$ and its related order. Then, a mask array ${w}_{i,j}$ is generated with limited binary values, preserving ${c}_{i,j}$ when ${w}_{i,j}=$ 1 and removing ${c}_{i,j}$ when ${w}_{i,j}=0$. In this case, we can achieve a frequency domain wavepacket whose Gabor coefficients are exactly inside area defined by the mask function ${w}_{i,j}$. The Gabor coefficients are corrected by ${w}_{i,j}$, and then a modified Gabor coefficient array ${\stackrel{~}{c}}_{m,n}$ is generated, shown as:
The frequency domain wavepacket corresponding to the order component of interest is reconstructed from the modified coefficient array ${\stackrel{~}{c}}_{i,j}$:
5. Application in processing multimode overlapping wave
In this section, the performance of the proposed approach in processing multimode overlapping wave is investigated with a numerical simulation. The software Disperse (Nondestructive Testing Laboratory, Dept. of Mechanical Engineering, Imperial College, London, UK) [22] is used to generated the wave numberfrequency function for Lamb waves in a 1mm thick aluminium plate and these are shown in Fig. 4. The Gaussian modulated wavepacket in Eq. (3) with centre frequency ${f}_{0}=$ 150 kHz, time width parameter $\sigma =$ 2×10^{5} is assumed to be the incident wave, which is shown in Fig. 5(a). Its frequency spectrum is shown in Fig. 5(b), where its bandwidth range extends about $\left[{f}_{\mathrm{m}\mathrm{i}\mathrm{n}},{f}_{\mathrm{m}\mathrm{a}\mathrm{x}}\right]$ = [80 kHz, 220 kHz].
Fig. 4The dispersion curve of wave numberfrequency
Fig. 5The incident signal: a) the incident signal in time domain, b) the incident signal in frequency domain
a)
b)
According to the dispersion curve in Fig. 3, there exist two modes of the fundamental symmetric and antisymmetric modes (S0 and A0) in frequency band [80 kHz, 220 kHz]. The group velocity of A0 at centre frequency 150 kHz is 2.069 m/ms, and that of S0 is 2.701 m/ms. The frequency range of received wave can be changed because of its dispersion nature. So, the cubic polynomials are used to fit the dispersion curves of A0 and S0 mode in [50 kHz, 250 kHz]. The number wave functions calculated are as follows:
${k}_{S}\left(f\right)=0.0031{f}^{3}0.0002{f}^{2}+0.1839f2\times 1{0}^{6},$
where ${k}_{A}\left(f\right)$ is the wave numberfrequency function of mode A0, and ${k}_{A}\left(f\right)$ is that of mode S0. The wave number functions of both S0 and A0 modes can be estimated with theoretical calculation combined with experiment verification in practical application. Based on Eq. (12), an artificial signal is generated by Eq. (13) with the propagating distance of A0 mode wavepacket ${x}_{1}=$ 300 mm, that of S0 mode wavepacket ${x}_{2}=$ 400 mm. And the reflection coefficients $C$ are set as “1”. The numerical overlapped and multimode signal $s\left(t\right)$ can be described by this model:
The sampling frequency and sampling period for the synthetic signal $s\left(t\right)$ is 1 MHz and 400 us respectively. The simulation result is shown in Fig. 6(a). The two wave packets are not resolvable in the original time trace. The frequency spectrum of $s\left(t\right)$ is shown in Fig. 6(b), where too many supernumerary frequencies make the identification of mode wave more difficult in frequency domain. Because of the dispersion and multimode overlapping, the shape of $s\left(t\right)$ frequency spectrum has been changed greatly, compared with that of incident signal in Fig. 5(b). It is almost impossible to identify which one is mode A0/S0 wavepacket, because they are completely overlapping in time and frequency domain. Without prior knowledge, it’s more impossible to calculate the propagation distance of each wave packet due to multiple and interference of overlapping time arrivals.
Fig. 6The numerical simulated signal: a) the numerical simulated signal, b) the frequency spectrum of numerical simulated signal
a)
b)
The order tracking method is performed on the signal in Fig. 6(b) to identify and separate the overlapped multimode wave packets according to the steps in Fig. 3. In order to identify the A0 mode wavepacket, the wave number function ${\theta}_{0}\left(f\right)=20{k}_{A}\left(f\right)$ is evaluated to be a referenced phase. The role of multiplying factor 20 is to make the order of A0 mode not too high. Considering the length of analysed signal, the maximum propagating distance of A0 mode is about 827.6 mm. The maximum order of A0 mode is about 42. Its corresponding resampling interval is $\mathrm{\Delta}\theta =$ 0.1496. It should be noted that at a smaller interval $\mathrm{\Delta}\theta $, the more accurate amplitude of resampled wave occurs, and a more accurate order spectrum is obtained, however, that a smaller resampling interval requires more computer resources for storing a larger quantity of raw data. After resampling and interpolation, the order power spectrum obtained is shown in Fig. 7, where there is a clear peak in the order spectrum at 15th order.
Fig. 7Order spectrum of the numerical simulated signal
The 15th order in Fig. 7 means that the frequency domain phase of A0 mode wave is 15 times of ${\theta}_{0}\left(f\right)$. Thus, the propagation distance of A0 mode wave is 20×15 mm, which completely equals the A0 mode propagation distance preset in the simulation. In the same way, if the wave number function of S0 mode is taken as the referenced phase, its propagation distance can be estimated from its corresponding order spectrum. Therefore, the individual mode can be identified since at one time the order spectrum only displays one kind of mode wave correlated to the selected wave number function.
It is noteworthy that the stronger the dispersion is, the more it can show the superiority of the proposed approach. As we all know, most conventional methods identify the wavepackets according to their arrival times, so the more severe dispersion phenomenon always bring more negative effects for the identification. However, the proposed approach just takes advantage of the dispersion nature to calculate the propagation distance, which would result in higher accurateness in structure damage location than the conventional methods. In the realsituation, when the theoretical wave numberfrequency function is more closely matched to its real value, the performance of the proposed approach becomes more superior, as the wave numberfrequency function is more representative of the phase change trend derived from dispersive signal.
Another significant advantage of the proposed method is to extract specified wave packet for further quantitative damage assessment. For the damage detection in SHM, the amplitude or shape of wavepacket reflected by damage has the close relation to damage severity. Such as for the delamination of composite structure, the energy of reflected wavepacket is as a function of the normalized length of the delaminated. If the wave packet reflected by damage can be extracted completely, it is much beneficial to quantitatively evaluate damage. Fig. 8 depicts the result of the Gabor expansionbased order tracking, where computed Gabor expansion coefficients of 15th order is characterised in the Gabor spectrogram. The colour in the Gabor spectrogram indicated the magnitude of the Gabor coefficients.
Fig. 8Gabor coefficient spectrogram of the numerical simulated signal
The order component in Fig. 8 represents the reflection related to the reference phase $l{k}_{A}\left(f\right)$, i.e. the mode A0 wave packet. We limit ${w}_{m,n}$ to binary values (either one or zero) which behaves as a mask, preserving ${c}_{m,n}$ when ${w}_{m,n}=$ 1 and removing ${c}_{m,n}$ when ${w}_{m,n}=$ 0. The inverse Gabor transform is performed using only Gabor coefficients which are exactly inside area defined by the mask function ${w}_{m,n}$. While the modified Gabor coefficients are obtained, the Gabor expansion equation is applied to reconstruct the waveforms of 15th order. The convolution result of the reconstructed order component and incident signal yields the time history of A0 mode wavepacket. Fig. 9 shows that the reconstructed signal superimposed with the original A0 mode wave packet.
The match between the reconstructed and original wave packet is very good except the amplitude attenuation and the tail of the wave packet. To quantify this similarity, a quantity called damage index (DI) is defined and calculated as:
where ${s}_{A}^{\text{'}}$ and ${s}_{A}\left(t\right)$ denote the reconstructed A0 mode wave and the original A0 mode packet. This method compares the amplitude of two sets of data. The calculated similarity is 85.78 %, which denotes the strong similarity of the two wavepackets. It is noteworthy that longer propagation distance will improve the resolution in Gabor spectrogram and provide better mode wavepacket separation, because the difference between the dispersion characteristics of two modes gets more apparent. The reconstructed wavepacket is much closer to the original wave packet, especially if propagation distance is large enough.
Fig. 9Comparison between the reconstructed wave and the original A0 mode wave
6. Application in processing strong dispersive signal
Because the proposed method is based on the dispersion nature of Lamb wave, the incident signal would not be constrained on the narrowband frequency to avoid dispersion effects. In order to verify the proposed algorithm’s potential capabilities in processing strong dispersive signal, the FEM experiment is conducted with an aluminium beam measuring L2000×W20×TH1 mm. The Lamb wave is generated at left beam end and collected at the distance 50 mm away from the exciting position. The damage measuring L2×W8×TH1 mm is introduced at distance 800 mm away from the left end. The beam is finely meshed using eightnote 3D brick elements with eight layers of elements in the beam thickness. The dynamic simulation is conducted using ABAQUES/EXPLICIT. The propagation paths of the Lamb wave are schematically depicted in Fig. 10.
To demonstrate this advantage, in the FEM simulation a broadband frequency signal is chosen as the incident signal, which is a Gaussian pulse function in Eq. (3) with attenuation coefficient $\sigma =$ 2×10^{5} at central frequency of 60 kHz. The group velocity of S0 at centre frequency 60 kHz is 5.437 m/ms. The Lamb wave signal captured is shown in Fig. 11(a), where the serious dispersion corrupts the signal and the reflections are smeared together. It is difficult to identify how many reflected wavepackets in the time waveform, not to mention to identify which mode each wavepacket belongs to. The frequency transform spectrum of received signal is shown in Fig. 11(b), from which it is clear the frequency domain signal is frequencyvarying nonstationary signal. Since the time waveform of received signal is known, its frequency spectrum covers the frequency range [0175 kHz].
Fig. 10Schematic illustration of propagation paths of Lamb waves in the damaged beam
Fig. 11The received signal: a) the received signal in time domain, b) the received signal in frequency domain
a)
b)
Therefore, the dispersion curve in Fig. 3 between the range [0175 kHz] is chosen as the reference phase. Two modes S0 and A0 mode should be considered in this frequency range. But the A0 mode usually presents higher attenuation during propagation because of the dominant outofplane movement of particles in the mode shape, which leaks partial energy to surrounding medium. Furthermore, A0 mode provides weaker reflections from the damage and boundary and insensitivity to damage in the structural thickness. Therefore, in this FEM simulation only the S0 mode needs to be considered. Performing the proposed method on the received signal, the reference phase ${\varphi}_{0}\left(f\right)=lk{}_{S}{}^{}\left(f\right)$ ($l=$250 mm is used to avoid the high order, $k{}_{S}{}^{}\left(f\right)$ is wave number –frequency function of S0 mode in range [0175 kHz]). The frequency domain signal is resampled and interpolated, and then its order spectrum obtained is shown in Fig. 12.
Fig. 12 provides three order peaks fairly separated and identifiable. The orders in Fig. 12 are identified as 0.206th, 0.6193th, 15.77th orders corresponding to the propagating distance 250×0.206 = 51.5 mm, 1548 mm, 3942.5 mm. Considering the geometry of the beam, the reflections sequentially are direct arrived wavepacket, the reflection from damage and that from the right end. It is obvious that the propagation distances estimated agree well with the theoretic values. Therefore, even in the case of strong dispersion the damage position can be identified precisely.
Fig. 12The order spectrum of received signal
Fig. 13TRFs of received signal: a) STFT, b) WVD, c) Chirplet decomposition, d) S transform
a)
b)
c)
d)
In order to verify the effectiveness of the proposed method, Fig. 13 provides the TFRs of the received signal. From the STFT and WVD spectrogram (Fig. 13(a) and (b)), it is difficult to identify the reflected wavepackets from the boundary and the damage owing to the limited timefrequency resolution and the crossterms. Because of the strong dispersion of guided wave, the wavepackes are overlapped together and cannot been separated from each other. In Fig. 13(c), the over decomposition of Chirplet transform causes the serious distortion of wavepackets. In the S transform of Fig. 13(d), it is easy to recognize the wavepackets, but the time centers of wavepackets deviate away from the actual propagation time. In Table 1, the estimated distances with the proposed method are tabulated and compared with those obtained with other conventional methods. It is verified that the propagating distances estimated with the proposed method agree with the theoretical values. But the estimated values with the conventional methods completely cannot present the practical propagation distances.
Table 1Comparisons of estimated distances between the proposed method and conventional ones
Method  Directed wavepacket (mm)  Damage wavepacket (mm)  Boundary wavepacket (mm) 
Theoretical distance  50  1550  3950 
Amplitude Enevlope  54.5  1844.3  4281.2 
S transform  282.6  1944.9  4338.8 
Chirplet decomposition  326.2  2120.5  4512.8 
Proposed method  51.5  1548  3942.5 
If mode conversion happens in the reflection, the frequency domain phase would be in an intricate shape and its shift related to dispersion curve would need to be recalculated according to the distance propagating as A0 mode and the distance propagating as S0. It is assumed that the mode S0 wave is converted into A0 mode after propagating over distance ${x}_{1}$ when it interacts with damage, and then it continues to propagate over distance ${x}_{2}$. Consequently, the resulting wavepacket reflected by damage ${S}_{x}\left(t\right)$ can be regarded as a timevarying product of shapeshifted waveform, i.e.:
where ${C}_{T}$ and ${C}_{R}$ represent the transmission and reflection coefficient respectively under the same assumption that the amplitude change of wave propagation is of insignificant. From Eq. (15), it can be seen that both S0 and A0 mode dispersion function are not suitable to be the reference phase as the linear relationship between the shift phase of ${S}_{x}\left(f\right)$ and any of wavenumber function dose not exist. In this case, it seems that the proposed method lost its effectiveness. One solution is to use some priori knowledge about propagating distances ${x}_{1}$ or ${x}_{2}$. Supposing ${x}_{2}$ is the linear function of ${x}_{1}$, which is expressed as ${x}_{2}=\alpha {x}_{1}+\beta $, the referenced phase:
According to Eq. (17), it is clear that the propagation ${x}_{1}$ can be calculated by taking ${\varphi}_{0}\left(f\right)={k}_{s}\left(f\right)+\alpha {k}_{A}\left(f\right)$ as reference phase with the proposed method.
It should be noted that the noises in the order spectrum (shown in Fig. 12) is caused by the operation of resampling and interpolation. In practical application, there are some factors that influence the effectiveness of the proposed method including the dispersion strength and the accurateness of dispersion curve, the resampling rate and interpolation algorithm, as well as the mode conversion in propagation and the frequency band of incident signal. The calculation of wave numberfrequency function is a key point for the application of the proposed algorithm, because the reference phase must be derived from it. In the practical application, the wave numberfrequency function obtained by the theoretical calculation can be different to that of actual structure tested. Therefore, the wave number function from theoretical calculation generally needs to be corrected by experimental measurements.
7. Conclusions
The dispersive and multimode nature of Lamb wave causes some difficulties in the interpretation of signals produced by damage structures. When a Lamb wave interacts with two closelyspaced discontinuities, the echo signal will be useless for processing. To resolve this problem, a specific signal processing scheme is proposed based on order tracking related to dispersion curve. In the proposed scheme, the unavoidable dispersion characteristic of Lamb waves was employed, but not evaded. The numerical and FEM simulation results demonstrate that the proposed method works with several merits for damage identification in SHM are summarized: (1) Dispersion presents no real obstacle to the high resolution of damage pinpointing in structures. The technique can provide better propagating distance of incident wave, without time of arrival and mode velocity. (2) Constraints on the choice of frequency range, where the dispersion effects are small, are relaxed. Using different frequency ranges may enhance the damage of various types. (3) The reflection wave packet from damage can be extracted for further analysis of damage degree.
One the other hand, it needs to be noted that the proposed scheme is limited by some drawbacks in practical applications. For instance, the performance strongly relies on the preciseness of dispersion curve and dispersion strength. But for the narrowband incident signal, the phenomenon in its propagation is not very strong, even its dispersion is weak enough to be neglected. And its detection accurateness is influenced by resampling frequency and interpolation. In the ongoing studies, the proposed algorithm is being validated experimentally. Future work will be carried out with the aim of overcoming its drawbacks mentioned above. The refinement of the technique will bring much beneficial for precise location of damage and quantitative describe of damage severity.
References

Miller C. A., Hinders M. K. Classification of flaw severity using pattern recognition for guided wavebased structural health monitoring. Ultrasonics, Vol. 54, 2014, p. 247258.

Veidt M., Liew C. Nondestructive evaluation(NDE) of aerospace composites: structural health monitoring of aerospace structures using guided wave ultrasonics. NonDestructive Evaluation (NDE) of Polymer Matrix Composites, Vol. 43, 2013, p. 449479.

Liu Xiaofeng, Bo Lin, Luo Hongling Compensation method of propagating distance for guided wave. IET Science, Measurement and Technology Vol. 9, Issue 2, 15, p. 915.

Liu Xiaofeng, Veidt M., Bo Lin Order analysis for measurement of Lamb wave propagation distance. Measurement Science and Technology, Vol. 24, 2013, p. 211226.

Zhang H., Wu S., Ta D., Xu K., Wang W. Coded excitation of ultrasonic guided waves in long bone fracture assessment. Ultrasonics, Vol. 54, 2014, p. 12031209.

Farhidzadeh A., Salamone S. Referencefreecorrosion damage diagnosis in steel strands using guided ultrasonic waves. Ultrasonics, Vol. 57, 2015, p. 198208.

Mijarez R., Baltazar A., RodríguezRodríguez J., Ramírez Niño J. Damage detection in ACSR cables based on ultrasonic guided waves. Dyna, Vol. 81, 2014, p. 226233.

Assous S., Elkington P. Modified S transform for timefrequency analysis of bore hole flexural waves dispersion. Journal of the Acoustical Society of America, Vol. 131, 2012, p. 3349.

Latif R., Laaboubi M., Aassif E., Maze G. Dispersion analysis of acoustic circumferential waves using timefrequency representations. Vibration and Structural Acoustics Analysis, Springer, Heidelberg, 2011, p. 183205.

Quaegebeur N., Masson P., Micheau P., Mrad N. Broadband generation of ultrasonic guided waves using piezo ceramics and subband decomposition. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control Vol. 59, 2012, p. 928938.

Olofsson T. Deconvolution and modelbased restoration of clipped ultrasonic signal. IEEE Transactions on Instrumentation and Measurement, Vol. 54, 2005, p. 12351240.

Sin SamKit, Chen ChiHay A comparison of deconvolution techniques for the ultrasonic nondestructive evaluation of materials. IEEE Transactions on Image Processing, Vol. 1, 1992, p. 310.

Kuttig H., Niethammer M., Hurlebaus S., Jacobs L. J. Modelbased analysis of dispersion curves using chirplets. Journal of the Acoustical Society of America, Vol. 119, 2006, p. 21222130.

Zhao M., Zeng L., Lin J., Wu W. Mode identification and extraction of broad band ultrasonic guided waves. Measurement Science and Technology, Vol. 25, 2014, p. 115005.

Raghavan Ajay, Cesnik Carlos E. S. Guidedwave signal processing using chirplet matching pursuits and mode correlation for structural health monitoring. Smart Materials and Structures Vol. 16, 2007, p. 355366.

Cai J., Shi L., Qing X. P. A timedistance domain transform method for Lamb wave dispersion compensation considering signal waveform correction. Smart Materials and Structures, Vol. 22, 2013, p. 13105024.

Liu L., Yuan F. G. A linear mapping technique for dispersion removal of Lamb waves. Structural Health Monitoring, Vol. 9, 2010, p. 7586.

Agrahari Jitendra Kumar, Kapuria Santosh A refined Lamb wave timereversal method with enhanced sensitivity for damage detection in isotropic plates. Journal of Intelligent Material Systems and Structures, Vol. 27, 2016, p. 12831305.

Chan Eugene, Rose Francis L. R., Wang Chun H. A comparison and extensions of algorithms for quantitative imaging of laminar damage in plates. II. Nonmonopole scattering and noise tolerance. Wave Motion, Vol. 66, 2016, p. 220237.

Demarchi L., Marzani A., Caporale S., Speciale N. Ultrasonic guidedwaves characterization with warped frequency transforms. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, Vol. 56, 2009, p. 22322240.

Yang Y., Peng Z., Zhang W., Meng G. Frequencyvarying group delay estimation using frequency domain polynomial chirplet transform. Mechanical Systems and Signal Processing, Vol. 46, 2014, p. 146162.

Pavlakovic B., Lowe M. J. S. Disperse Software Manual Version 2.0.16B. Imperial College London, UK, 2003.
About this article
This work is supported by the National Natural Science Foundation Project (No. 51475052 and No. 51675064), the China Postdoctoral Science Foundation (No. 2015M582519 and No. 2016790833) and the Fundamental Research Funds for the Central Universities (No. 106112016CDJZR115502).