Lamb wave inspection for composite laminates using a combined method of sparse reconstruction and delay-and-sum

Composite Structures 223 (2019) 110973

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Lamb wave inspection for composite laminates using a combined method of sparse reconstruction and delay-and-sum

T

⁎

Caibin Xua, Zhibo Yanga, , Shaohua Tiana, Xuefeng Chena,b a b

School of Mechanical Engineering, Xi’an Jiaotong University, Xi’an 710049, China State Key Laboratory for Manufacturing Systems Engineering, Xi’an Jiaotong University, Xi’an 710049, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Lamb waves Sparse reconstruction Composite laminates Damage detection Structural health monitoring

A narrowband tone burst excitation is usually required in most Lamb wave inspection techniques in order to reduce the dispersion effect and maintain mode purity, especially in the case that dispersion compensation is hard to implement, such as the scenarios prior knowledge of accurate dispersion curves is not available. The duration of the excitation is relatively large because of its narrow bandwidth, which degrades the resolution for inspection. Following the emerging sparse representation theory, a combined method of sparse reconstruction and delay-and-sum (DAS) for high resolution Lamb wave inspection is proposed in this paper. The scattering signals are compressed in time domain without losing the information of the arrival time of each echo. Unlike most methods that signal decomposition and reconstruction use the same dictionary, the proposed approach utilizes two different dictionaries: one for signal decomposition and another for feature extraction. The DAS imaging method combined with the feature signals is applied to damage imaging for carbon fiber reinforced plastic (CFRP) laminates. The effectiveness of the proposed approach is validated on a quasi-isotropic laminated CFRP plate. The imaging results show that the proposed approach can achieve a higher localization resolution and lower noise floor without degrading the localization accuracy.

1. Introduction Composite laminates are widely used in different fields, such as aerospace and wind turbine. While owing to many good mechanical properties, such as high strength to weight ratio and corrosion resistance, composite laminates can be also failure when damages such as delamination [1,2] occur. Lamb wave, which is known as guided wave propagating in solid thin plates, is one of the most prominent tools for active damage detection and structural health monitoring (SHM) [3–7]. In Lamb wave based SHM, an array of permanently affixed actuators and sensors, like piezoelectric lead zirconate titanate (PZT) [8], are commonly used to transmit and record Lamb waves so as to assess the health status of structures. As the transmitter-receiver pair can cover a relatively large-scale area, only a few elements are required. In general, a set of baseline measurement signals needs to be recorded when the structure is in health condition. The baseline-subtracted signals (referred to as scattering signals) are then used to evaluate the location, type and severity of damages. In Lamb wave based SHM, damage localization resolution is closely related to the excitation. To achieve mode purity and suppress dispersion, a narrowband tone burst excitation is commonly transmitted [9].

⁎

The duration of the excitation cannot be very small because it must balance with the bandwidth as wider bandwidth can lead to more serious dispersion effect and less mode purity. Thus in many Lamb wave applications, the number of cycles is set to be five for the above reason [10–14]. However, even in the case of a specific tone burst excitation, the echoes of scattering signals may overlap together when scatterers are close in a structure, which limits the damage localization resolution. In addition, the priori information of dispersion curves of a structure is not always available, which causes some broadband excitation methods like pulse compression [15,16] inapplicable to damage detection. Sparse reconstruction, as a part of the field of the compressed sensing (CS) [17,18], is a new research field in mathematics and informatics, and has attracted much attention in ultrasonic guided wave based damage detection [19,13,20]. Some greedy algorithms, such as the orthogonal matching pursuit (OMP) [21] and the compressive sampling matching pursuit (CoSaMP) [22], can be taken to find an approximate solution for a specific sparse reconstruction problem. One of the key issues of compressed sensing based ultrasonic damage detection is to build an efficient overcomplete dictionary according to actual problems. In order to decompose the scattering signals, some dictionary-based methods have been developed. Mallat and Zhang [23]

Corresponding author. E-mail address: [email protected] (Z. Yang).

https://doi.org/10.1016/j.compstruct.2019.110973 Received 24 October 2018; Received in revised form 24 February 2019; Accepted 8 May 2019 Available online 11 May 2019 0263-8223/ © 2019 Elsevier Ltd. All rights reserved.

Composite Structures 223 (2019) 110973

C. Xu, et al.

s (t ) = w (t ) eiωc t .

developed a method called matching pursuit (MP) to decompose signals with a redundant dictionary of Gabor functions. Hong et al. used a dispersion-based chirp function dictionary [24] and a Gabor dictionary [25] for guided wave damage detection by the MP method. Tse and Wang [26] used MP decomposition with an optimized dictionary-based on two interfering reflection components to evaluate the length of defect. Demirli and Saniie [27] presented a model-based estimation pursuit method that made use of statistical estimation principle in echo matching. Lu and Michaels [28] presented a numerical implementation method of MP with a Gabor dictionary for the analysis of ultrasonic signals. Zhang et al. [29] studied the selection of overcomplete dictionaries for ultrasonic nondestructive evaluation. Marchi et al. [30] proposed a warped basis pursuit algorithm with the dictionaries composed by optimized atomic functions for Lamb wave damage detection. Xu et al. [31] presented a sparse reconstruction based method to compensate the dispersion effect of Lamb waves. All the above methods require the dispersion curves of guided wave modes to construct an overcomplete dictionary and then to sparsely decompose the response signal. To some extent, those methods get a certain degree of success in specific conditions. However, the prior information of dispersion curves is hard to obtain in Lamb wave inspection for composite laminated structures because of the uncertainty and anisotropy of composite material parameters. The lack of accurate dispersion curves makes it impossible for the above methods to be directly applied to composite laminated structures. A sparse reconstruction based approach using two dictionaries for high resolution damage inspection is presented in this paper. In contrast to Gabor or some other wavelet dictionaries, the first dictionary used here is constructed using the excitation and its time shifted and phase varied signals so as to be independence on dispersion curves. Another dictionary is constructed using the same way as the first one but each of its atom holds a shorter duration. Under a narrowband tone burst excitation, each echo of the scattering signals can be approximately treated as a time shifted and phase varied excitation with additional amplitude attenuation. The time of flight (TOF) information of each echo is extracted using the first dictionary by solving a sparse representation problem. Then feature signal is extracted using the second dictionary. In this way, the scattering signals are compressed using the two dictionaries. This approach is referred to as matching compression (MC) in the rest of the paper. Finally, the delay-and-sum (DAS) [32] imaging method combined with the MC approach is used to locate damages. The rest of the paper is organized as follows. In Section 2, Lamb wave propagation model in the linear wavenumber assumption is introduced. The proposed approach, signal matching compression by sparse reconstruction, is described in Section 3. Section 4 describes the experimental setup and procedure. Section 5 provides the imaging results and discussion. Finally, the conclusion is in Sections 6.

(1)

Here ωc is the center angular frequency of the tone burst excitation, and w (t ) is a window function. In this research, the Hanning window function is selected to be w (t ) , which can be expressed as

f n ⎤ ⎡ w (t ) = ⎢H (t ) − H ⎜⎛t − ⎟⎞ ⎥ ⎜⎛1 − cos ⎛2π c t ⎞ ⎟⎞, fc ⎠ ⎝ n ⎠⎠ ⎝ ⎝ ⎣ ⎦ ⎜

⎟

(2)

where H (t ) is the Heaviside function, n is the number of cycles, and fc is the center frequency. Under assumptions of linear system and single mode, the response of Lamb wave between any two points in frequency domain can be expressed as

Y (ω) =

r0 S (ω) e−ik (ω) r . r

(3)

Here Y (ω) is the Lamb wave response in frequency domain corresponding to a traveling distance r, r0 is a constant and referred to as reference distance, k (ω) is the wavenumber of Lamb wave and frequency-dependent, i is the imaginary unit, and S (ω) is the corresponding excitation in frequency domain. According to Eq. (1), the tone burst excitation in frequency domain can be expressed as

S (ω) =

∞

∞

∫−∞ s (t ) e−iωtdt = ∫−∞ w (t ) e−i (ω−ω ) tdt = W ⎛ω − ωc⎞, c

⎜

⎟

⎝

⎠

(4)

where W (·) is the Fourier transform of the window function w (t ) . Eq. (4) shows that, the excitation is equal to the window function with an additional frequency shift. Through a Taylor expansion, the frequency-dependent wavenumber k (ω) can be expressed as

k (ω) = α + β (ω − ωc ) + o (ω − ωc ),

(5)

where

α = k (ωc ) =

ωc , cp (ωc )

β = k ′ (ωc ) =

dk (ω) dω

(6)

= ω = ωc

1 , cg (ωc )

(7)

and o (ω − ωc ) is infinitesimal of (ω − ωc ) , which supports the relationship

lim ω → ωc

o (ω − ωc ) = 0. ω − ωc

(8)

In Eqs. (6) and (7), cp and cg are the phase velocity and the group velocity, respectively. They are defined as

cp (ω) =

ω , k (ω)

(9)

cg (ω) =

dω . dk (ω)

(10)

2. Lamb wave propagation with the linear wavenumber If the wavenumber is a linear function of frequency, it will be referred to as the linear wavenumber in the following sections. Consider that two sensors, which are affixed to the same surface of a solid thinwall plate with a fixed thickness, are employed as a transmitter and a receiver, respectively. The transmitter is regarded as a point-like source, which transmits waves in all directions. Because an approximately single zero-mode Lamb wave is possible to be transmitted by optimizing the center frequency and bandwidth of the excitation [33], Lamb wave with only a single mode is considered here. The entire system of Lamb wave propagation can be regarded as a linear system. In practice, both the part of unmodeled mode Lamb waves and the nonlinear part are treated as noise. In this research, only a narrowband tone burst excitation is used in order to reduce the effect of dispersion. In general, the tone burst excitation s (t ) in complex domain can be expressed as

Under the narrowband tone burst excitation, the frequency-dependent wavenumber k (ω) is treated to be a linear function of frequency in this study. In practice, the linear wavenumber assumption will be achieved by treating the high order components of k (ω) to be additional interference. Under this assumption, the component o (ω − ωc ) in Eq. (5) is discarded. Substituting Eqs. (4) and (5), into Eq. (1), the response of Lamb wave in time domain can be expressed as

y (t ) =

r0 −1 r F {W (ω − ωc ) e−i [α + β (ω − ωc )] r } = 0 e−iαr eiωc t w (t − βr ), r r (11)

where F−1 is the inverse Fourier transform. According to Eq. (1), Eq. (11) can be expressed as 2

Composite Structures 223 (2019) 110973

C. Xu, et al.

y (t ) =

r0 −iαr r ⎞ e s (t − βr ) eiωc βr = r0 r −0.5e−iαr eiωc βr s ⎛⎜t − ⎟. r c ( ωc ) ⎠ g ⎝

⎤ ⎡ Dj = ⎢ em − j + 1, em − j + 2, ⋯, em , 0, ⋯, 0⎥.  ⎥ ⎢  j ⎦ ⎣

(12)

Because r0 r −0.5e−iαr eiωc βr is time-independent, the envelope shape of y (t ) is the same as the excitation s (t ) except a constant amplitude factor r0 r −0.5 . Eq. (12) shows that under the linear wavenumber assumption, the signal y (t ) is non-dispersive. y (t ) is just equal to the excitation with a time shift, a phase variation, and additional amplitude attenuation. If the time shift is determined, which is referred to as TOF Δt , then the traveling distance r can be obtained by multiplying Δt and the group velocity cg (ωc ) . Note that when the wavenumber k (ω) is not linear, some dispersion compensation [11,31,34,35] techniques still can be used to get the same results as Eq. (12) because the dispersive wave packets will be compensated into the waveform of the excitation after dispersion compensation.

Here ej is a column vector with the dimension m × 1 and the its jth element is equal to 1 while others are zeros, i.e.

ej = (0, ⋯, 0, 1j , 0, ⋯, 0)T.

(19)

Then a sub-dictionary corresponding to a specific φ can be built using the m − m 0 possible responses (referred to as atoms) in the following way

A φ = [a1 a2 ⋯ am − m0].

(20)

For the convenience of the subsequent processing, the sub-dictionary needs to be ℓ2 -norm normalized, i.e.

A φ = [a1 a2 ⋯ am − m0] ·(‖a1‖2 )−1. 3. Combined method of sparse reconstruction and delay-and-sum

In Section 2, under the linear wavenumber assumption, the response can be treated as a weighted combination of the time shifted excitation with additional phase variation. Based on this assumption, this section will compress the echoes of the response into shorter ones in time domain, which will improve inspection resolution. Suppose that the recorded signal y (t ) has a length of m points when the signal is digitized. With the same sampling frequency fs , the same center frequency fc , and the same number of cycles n 0 , the excitation signal has a data length of m 0 points. According to Eq. (2), the data length of excitation m 0 is equal to

n 0 fs fc

.

A = [A φ1 A φ2 ⋯ A φp ].

y = Ax, (13)

(23)

)T

where x = (x1, x2, ⋯, x p (m − m0) is a column vector. Each element in x represents the corresponding weighted factor. For example, the element x j is the weighted factor of the jth atom in A . So each atom in A is corresponding to a unique time shift τ and a unique phase variation φ . In general, the signal y is a weighted combination of only a few atoms in A because the locations of scatterers in a plate are sparse compared with the whole plate. That is, the most elements in vector x are zeros and the number of nonzero elements is sparse compared with the dimension of the vector x . Then y can be expressed using only a few atoms in A and thus the corresponding weighted factor vector x can be obtained by solving the following problem

x ̂ = arg minx ‖x‖0 subject to y = Ax.

(15)

Here the superscript T represents the transposition operation. Then the other response signals aj (t ) (2 ⩽ j ⩽ m − m 0) with different TOFs can be expressed as a time shift τj of a1 (t ) , i.e.

aj (t ) = a1 (t − τj ),

(16)

where τj = (j − 1)/ fs . The corresponding discrete form of aj (t ) can be expressed as

aj = Dj a1,

(24)

Once the problem of Eq. (24) is solved, the TOF of each echo is obtained. Eq. (24) is a sparse representation problem, and it is hard to get the exact solution. In fact, it has been proven to be a NP-hard problem [36]. According to the compressed sensing theory, approximate solution of Eq. (24) can be obtained using pursuit algorithms [36]. In this

(14)

where φ is a phase variation, and it can be discretized in the range of 0 to 2π , so as to match different phase variations. When φ is a determined value, a1 (t ) can be discretized and denoted as a column vector a1 with the dimension m × 1

a1 = (a1, a2, ⋯, am0 , 0, ⋯0)T.

(22)

A schematic diagram of the sub-dictionary A φ1 is shown in Fig. 1. Note that the dictionary can match different phase variations of wavepackets because of the existence of the parameter φ in Eq. (14). In general, the damage scattering signal y may contain several echoes, and it can be expressed as a weighted combination of all the possible atoms of dictionary A , which is

Since m is determined, there are m − m 0 possible TOFs in total, and are corresponding to m − m 0 possible traveling distances. Note that m must be large enough to cover the whole monitoring area. The possible response signals are modeled now. When Lamb wave is transmitted and recorded in the same position, which means that the traveling distance is equal to zero, the response is equal to the excitation (ignore difference of amplitude) with zero-padding in order to reach the length of m points. Just like Eq. (12) shows, the phase of the response signal is changed because e−iαr eiωc βr is not equal to zero, then the recorded signal can be denoted as

f ⎞ ⎛ n ⎤⎛ ⎞ ⎡ a1 (t ) = ⎢H (t ) − H ⎜⎛t − 0 ⎟⎞ ⎥ ⎜1 − cos ⎜⎛2π c t ⎟⎞ ⎟ sin ⎜2πfc t + φ⎟, f n 0 ⎠⎠ c ⎠⎦⎝ ⎝ ⎝ ⎣ ⎝ ⎠

(21)

In Eq. (21), because of ‖a1‖2 ≡ ‖aj ‖2 (1 ⩽ j ⩽ m 0) , each column in A is ℓ2 -norm normalized. For different φ in the range of 0 to 2π , a corresponding sub-dictionary A φ can be obtained using the above procedure. Then the phase variation parameter φ is discretized into p values: φ1, φ2 , ···,φp . Finally the p sub-dictionaries are combined together, which can be expressed as

3.1. Signal matching compression by sparse reconstruction

m0 =

(18)

(17)

where Dj is a forward shift matrix with the dimension m × m and can be expressed as

Fig. 1. A schematic diagram of dictionary A φ1 (n 0 = 5). 3

Composite Structures 223 (2019) 110973

C. Xu, et al.

research, the OMP algorithm is used to solve the problem of Eq. (24) and thus the wave packets in y can be extracted. In order to extract the feature signal from the original signal y to achieve a shorter duration, another atom and its time shifted and phase varied signals are used to build another dictionary B in the same way as building the dictionary A . The new atom is defined as

b1 (t ) 2πf (t − τ ′) ⎞ ⎤ n ⎞⎤ ⎡ ⎛ ⎞ ⎛ = ⎢H ⎜t − τ ′⎟ − H ⎜t − τ ′ − 1 ⎟ ⎥ ⎡1 − cos ⎜⎛ c sin ⎟ ⎥ fc ⎢ n1 ⎝ ⎠⎦ ⎣ ⎣ ⎝ ⎠ ⎝ ⎠⎦ ⎛ ⎛ ⎞ ⎞ ⎜2πfc ⎜t − τ ′⎟ + φ⎟, ⎝ ⎝ ⎠ ⎠

(25)

Fig. 3. A schematic diagram of dictionary B φ1(n1 = 1) .

where τ ′ = (n 0 − n1)/ fc is a time delay, which makes the curves of a1 (t ) and b1 (t ) to center at the same position on the time axis. For a determined φ , the discrete form of b1 (t ) can be expressed as b1

and the sparse coefficient x ̂

y ̂ = Bx.̂

T

b1 = ⎡  0, 0, ⋯ , 0 , b1, b2, ⋯, bm1, 0, ⋯0⎤ ⎢ (m0 −m1)/2 ⎥ ⎣ ⎦

, (26)

where

m1 =

n1 fs fc

,

(27)

where n1 is a positive number and is smaller than n 0 , so that b1 has a shorter duration than a1. b1, b2 , ⋯, bm1 are digitization values of the imaginary part of the signal defined by Eq. (1) with parameter n = n1. As shown in Fig. 2, both the waveform of a1 and b1 within 0–60 μs are plotted. After center alignment, the maximum value of a1 and b1 are corresponding to the same time index. Similar to the process of defining aj , the other atoms can be defined as

bj = Dj b1 (1 ⩽ j ⩽ m − m 0).

Algorithm 1. Matching Compression (MC) Task: Compress duration of each packet in the signal y . Input: Dictionaries A and B , signal y , the sparsity L′, error threshold ε .

(28)

Initialize:l ← 0, x 0 ← 0, r 0 ← y ,

S 0 ← support{x 0} = ∅ . Main iteration: Whilel ⩽ L′, do l ← l + 1. 1. Find j0 :

Then a normalized sub-dictionary Bφ for a specific φ is

Bφ = [b1 b2 ⋯ bm − m0] ·(‖b1‖2 )−1.

(29)

When considering other values for φ in the same range from 0 to 2π , all the sub-dictionaries can be combined together as

B = [Bφ1 Bφ2 ⋯ Bφp ].

(31)

Because x ̂ contains the TOF information of the recorded signal y , and each atom of B have a shorter duration compared with A , each echo of y ̂ has a shorter duration than the original signal y . So y ̂ can be seen as a compression of y . The overlapped wave packets in y may possible to be separated by the above procedure. Since solving x ̂ in Eq. (24) can be treated as a matching process in time domain, and the feature signal y ̂ in Eq. (31) can be treated as a compression process, we refer to this approach as matching compression. Note that the method presented above does not consider anisotropy of guided waves. A pseudo-code of the matching compression is given in Algorithm 1.

∀ j ∉ Sl − 1, φ (j0 ) ⩾ φ (j ) , where φ (j ) = | 〈aj, r l − 1〉 |.

(30)

2. Update: Sl = Sl − 1 ∪ {j0 } . 3. Solve the least square problem:

The dictionaries A and B have the same dimension m × p (m − m 0) . A schematic diagram of sub-dictionary Bφ1 is shown in Fig. 3. Each atom in B can be seen as a compression of the corresponding atom in A . In ultrasonic imaging [32], the duration of the echo is crucial to resolution. A shorter duration of each echo is required for high resolution imaging. So the feature signal y ̂ can be obtained using the dictionary B

x l = arg minx ‖y − Ax‖22 s. t . support{x} = Sl . 4. Update: r l = y − Ax l . 5. ifr l < ε , then

x ̂ ← x l ; break end if end while Restructure: Calculate the feature signal: y ̂ = Bx .̂ Output: Feature signal y .̂

3.2. Damage imaging by the DAS method In order to locate the damages, an ultrasonic imaging method known as DAS [32] is introduced

⎛ P x, ⎜ ⎝

⎞ y = ⎟ ⎠

L

∑ j=1

y ⎛⎜ ⎝

rTj + rRj ⎞ ⎟, cg ⎠

(32)

rTj and

where P (x , y ) is the pixel value of imaging pixel in location (x , y ), rRj are the distances from the transmitter and the receiver of the jth sensor pair to location (x , y ) , and y (·) is the scattering signal. In practice, y (·) is usually replaced by its Hilbert envelope for discarding phase information. In this research, y (·) is the Hilbert envelope of the feature signal obtained from the proposed MC approach. As this process combines the DAS and the MC, it is referred to as MCDAS in this paper. The

Fig. 2. The waveforms of the two signals and the corresponding envelopes. The red one is time shifted in order to align center with the blue one. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 4

Composite Structures 223 (2019) 110973

C. Xu, et al.

the edge reflection, signals acquired from the transmitter-receiver pair i − j and j − i are treated to be identical). Each excitation is generated by an NI PXI-5412 arbitrary waveform generator and then amplified by a Piezo Systems EPA-104 linear amplifier. The magnitude of the excitation amplitude voltage is 120 V (peak-to-peak value). The signals are recorded by four NI PXI-5122 digitizers (each digitizer supports 2 recording channels). The software system is developed using LabView. In the experiment, signals are recorded in three different states of the plate: 1) no damage state (baseline state); 2) one damage state (#1 damage, case I); and 3) two damages state (# 1 and # 2 damages, case II). Scattering signals are obtained by subtracting the baseline signals from the corresponding signals acquired in case I and case II, respectively. 4.2. Group velocity determination The group velocity needs to be calibrated. CFRP laminated composite plates exhibit anisotropy nature, which means that the dispersion curves of the material are direction dependent. So the group velocities are also direction dependent. In addition, when the mechanical property parameters of the material are not accurately determined, the group velocities calculated by some theories are still not the same as the real ones. Considering the direction dependent nature, the group velocity calibration [38] is performed using the baseline signals by solving the problem

Fig. 4. Schematic diagram of the laminated composite plate with eight affixed sensors and two simulated damages. The circles denote the locations of sensors, and the asterisks denote the locations of simulated damages.

imaging results of both the DAS and the MCDAS will be present in the next section. 4. Experimental setup and procedure

cg = arg max c 4.1. Signal acquisition

Sensors

Simulated damages

1 2 3 4 5 6 7 8 #1 #2

0 100 100 100 0 −100 −100 −100 25 50

−100 −100 0 100 100 100 0 −100 −25 −75

j=1

dj sj ̂ ⎛t + ⎞ dt , c⎠ ⎝ ⎜

⎟

(33)

In order to evaluate the imaging performance, three metrics are defined: the localization error, the estimated damage spot size, and noise floor. Firstly, the localization error re is defined to be the distance from the actual damage location (x ∗, y∗ ) to the imaging damage location (x 0∗, y0∗ ) :

re =

(x ∗ − x 0∗)2 + (y∗ − y0∗ )2 .

(34)

(x 0∗, y0∗ )

is defined to be the location of the The imaging damage location maximum pixel value in the entire image when there is only one damage. When there is multiple damages, (x 0∗, y0∗ ) is the location of the local maximum pixel value. A smaller re denotes higher localization accuracy. Secondly, the localization resolution can be evaluated by defining the estimated damage spot size:

Sps = Y(mm)

1

4.3. Damage imaging and imaging performance

Table 1 Locations of sensors and simulated damages. X(mm)

P

̂ is the analytic representation of the baseline signal, dj is the where sj (·) distance from the transmitter to the receiver for the jth sensor pair, and t1, t2 are the lower and upper bound of the integration period, respectively. Note that through implementing Eq. (33), the group velocities in all directions are treated to be identical. In the experiment, the integration period is set to be the duration of the excitation, which is 0–62.5 μs when center frequency fc = 80 kHz and number of cycles n = 5. The result of the group velocity calibration is found to be 1430 m/s (A0 mode).

The experiment specimen is a [0/45/ −45/90]2s CFRP laminated composite plate with a dimension of 405 mm × 405 mm × 2 mm. An array of eight PZT disks with 0.5 mm thickness and 8 mm diameter are affixed to the plate with a square configuration. Two cylindrical neodymium magnets with 10 mm length and 10 mm diameter are attached to the plate as the simulated damages. Fig. 4 shows the schematic diagram of the laminated composite plate with the eight PZT disks and two simulated damages. The origin of coordinates is set at the center of the plate. The locations of PZT disks and the simulated damages are presented in Table 1. In this experiment, a chirp signal with the bandwidth from 1 kHz to 250 kHz and 100 μs duration is transmitted and an 80 kHz, five-cycle Hanning windowed tone burst responses are obtained by post-processing technique [37]. Using the post-processing technique, the extracted responses from chirp responses are the same with the directly measured responses to the tone burst signal (80 kHz center frequency, five-cycle number Hanning windowed tone burst). Within the bandwidth of the tone burst signal, the wavenumber of A0 mode is a linear function of frequency approximately. The A0 Lamb mode is dominant in the plate under these conditions. The recorded signals are acquired in a roundrobin fashion: one of the eight PZT disks takes turn to be the transmitter, the other seven serve as receivers. By this way, a total of L = 8 × 7/2 = 28 unique signals (baseline-subtracted signals) are obtained (under the assumption of linear system and without considering

Number

t2

∫t ∑

∑ P {(x , y)}

s. t. P (x , y ) ⩾ max{P (x , y )} × Th,

(35)

where P {(x , y )} is an operator for getting the area of imaging pixel (x , y ) , and Th is a threshold. In this research, Th is set to be 0.5, which means that the defined damage spot size is the total area where pixel value is equal or greater than the half-amplitude of P (x , y ) . A smaller Sps denotes a greater confidence in localization resolution, and a greater imaging energy focusability. A smaller Sps also means that we can differentiate the damage location and damage-free location more precisely. In this way, the first metric re evaluates the accuracy, and the second metric Sps evaluates the precision. Third, the noise floor is defined as 5

Composite Structures 223 (2019) 110973

C. Xu, et al.

Fig. 5. (a) The normalized scattering signal of sensor pair 1–2 in the state of one simulated damage. (b) The feature signal extracted from (a) by the presented MC algorithm. (c) The normalized scattering signal of sensor pair 1–2 in the state of two simulated damages. (d) The feature signal extracted from (c) by the presented MC algorithm.

Nf = 1 −

∑ P (x ′, y′) , ∑ P (x , y )

the noise floor of the DAS is about 1.6 times as large as the MCDAS’s. The results show that the MCDAS method has a more excellent imaging performance than the DAS method in this case.

(36)

where P (x ′, y′) is the pixel value that greater than the threshold value max{P (x , y )} × Th . A smaller Nf denotes a lower noise floor of an image. This means that we can find the damage location more easily in the image in the case of a smaller Nf.

5.2. Case II: two simulated damages Here the case of two simulated damages is investigated. Fig. 5(c) and (d) show the original scattering signal and the corresponding extracted feature signal of sensor pair 1–2 in case II, respectively. As the two figures show, the first two echoes (locate in 0.1 ms and 0.15 ms roughly) are overlapped with each other, while they are separated after processing by the MC algorithm. So it suggests that the MCDAS method has better damage localization performance than the DAS method. The imaging results using the methods of DAS and MCDAS are shown in Fig. 7(a) and (b), respectively. In Fig. 7(b), the imaging pixel values corresponding to damages # 1 and #2 are equal to 1 and 0.872, respectively. So the ratio of the smaller pixel value to the larger one is 87.20% for the two damages in the MCDAS method. However, in Fig. 7(a), the corresponding ratio is found to be only 71.92% in the original DAS method. It means that the smaller pixel value is not easy to be recognized by the DAS method because the smaller one would likely be treated to be noise when the two pixel values are much difference. The corresponding filtered images for the two methods using thresholds of 0.5 and 0.75 are shown in Fig. 7(c)–(f). As one can see that, when the threshold reaches to be 0.75, only one damage can be detected in the filtered image for the DAS method, as shown in Fig. 7(e). And in the corresponding result of the MCDAS method, as shown in Fig. 7(h), both of two damages are detected. So the original DAS method may lead to missed detection easily for multiple damages. The presented method can reduce the missed detection because the presented method can reduce the gap between the pixel values corresponding to the two damages. The values of the three metrics are presented in Table 2. The total localization error of the DAS and the MCDAS methods in case II are 16.1 mm and 9.2 mm, respectively. The spot size of Fig. 7(d) is only about 30.51% of that in Fig. 7(c), i.e. the localization resolution of the MCDAS method is about 3 times higher than that of the DAS method. Similar to the case I, the noise floor of the MCDAS method is also lower than of the DAS method in this case.

5. Results and discussion In this section, the scattering signals are decomposed and feature signals are extracted. Then both the MCDAS and the DAS methods are used to generate images. The imaging performance of the two methods is evaluated by the three metrics. Finally, the performances of the MC algorithm in several cases are also investigated. 5.1. Case I: one simulated damage Only one simulated damage (damage # 1) is considered here. As an example, Fig. 5(a) shows the scattering signal of sensor pair 1–2, and Fig. 5(b) shows the corresponding feature signal obtained by the MC algorithm. From Fig. 5(b), the duration of each wave packet is significantly smaller than the original. The imaging results of the DAS and the MCDAS are shown in Fig. 6(a) and (b), respectively. The actual simulated damage location is denoted by a triangle, and the locations of the PZT disks are denoted by white dots. In addition, the corresponding filtered images using a threshold (all pixel values smaller than the threshold are set to be zeros) are shown in Fig. 6(c) and (d), respectively. Note that when generating images, the plate is discretized into a 2 mm × 2 mm grid, which means that the range resolution is 2 mm. With the comparison of Fig. 6(c) and (d), one can see that both two methods can roughly locate the damage location. From the comparison one can see that, the result of the MCDAS holds a relatively smaller spot size, which means it achieves a higher resolution. The results of the three metrics for quantifying imaging performance are given in Table 2. The localization errors of the DAS and the MCDAS are almost equal. The difference between them is only 0.1 mm, which is much smaller than the grid size. However, the spot size Sps of the DAS is about 4 times as large as the MCDAS’s. Note that smaller Sps value denotes higher localization resolution and energy focusability of an image. In addition, 6

Composite Structures 223 (2019) 110973

C. Xu, et al.

Fig. 6. Normalized imaging results for one simulated damage. (a) The original imaging result using the DAS method; (b) The original imaging result using the MCDAS method; (c) The filtered image for (a) using a threshold of 0.5; (d) The filtered image for (b) using a threshold of 0.5. The white dots denote the PZT locations and the triangles denote the actual simulated damage location.

one can see that, the two overlapped wave packets are separated. However, the amplitude of the wave packet #1 becomes smaller with the number of cycles increases. This result suggests that the presented algorithm cannot achieve good performance when the two wave packets are overlapped seriously. Actually, it is a hard problem in guided wave based SHM to separate two wave packets in such seriously overlapped situation. Note that the parameter n 0 will not always be 5, and it should be set to be the corresponding number of cycles when designing the dictionary A . The parameter n1 is always set to be 1. The performance in the case of tone burst excitations with more higher center frequency is also investigated. The scattering signals to the tone burst excitations with higher center frequencies, from 80 kHz to 240 kHz with a step of 40 kHz, are obtained in the state of two simulated damages. The transmitter and receiver are also PZT 1 and PZT 2, respectively. As the center frequency increases in the range of 80 kHz to 240 kHz, the proportion of the A0 mode guided wave decreases while the S0 mode increases. As shown in Fig. 9(a), when the center frequency is 80 kHz, the scattering signal contains hardly any S0 mode guided wave, and the two wave packets are both A0 mode (the wave modes are determined by calculating the group velocity for each wave packet because the group velocity of S0 mode is much larger than the A0 mode). When the center frequency increases to 240 kHz, the S0 mode is dominant in the scattering signal. The two original scattering wave packets reflected from damages #1 and #2 are overlapped because of dispersion effect, as shown in Fig. 9(a). The corresponding result processed by the MC algorithm can be seen in Fig. 9(b) when the center frequency is 240 kHz. There are several wave packets besides the #1 and #2 wave packets. The similar result can be seen when the center frequency is 200 kHz, which both the A0 and S0 modes are dominant.

Table 2 Results of imaging performance.

One damage Two damages

Method

Sps (mm2)

Nf

re (mm)

DAS McDAS DAS MCDAS

1504 360 3356 1024

6.90% 4.35% 9.44% 5.49%

3.1 3.2 4.2(#1), 11.9(#2) 5.1(#1), 4.1(#2)

5.3. Discussion In order to study the performance of the presented MC algorithm for separating two overlapped wave packets in more other cases, several investigations are implemented. The performance in the case of tone burst excitations with more larger number of cycles is investigated. The scattering signals to tone burst excitations with different number of cycles, from 5 to 15 with a step of 2, are obtained in the state of two simulated damages. Excitations are transmitted in PZT 1, and responses are acquired in PZT 2. Fig. 8(a) shows the normalized scattering signals to the excitations with different number of cycles. In Fig. 8(a), the scattering wave packets reflected from damages #1 and #2 are marked with the notations “#1” and “#2 ” when number of cycles is 5. From Fig. 8(a) one can see that, the two scattering wave packets are gradually overlapped together as the number of cycles increases. Because the wave packets are wider in time domain and easier to overlap as the number of cycles increases. Ones can hardly distinguish the two wave packets when the number of cycles is larger than 9. Fig. 8(b) are the corresponding feature signals extracted by the presented MC algorithm. From Fig. 8(b) 7

Composite Structures 223 (2019) 110973

C. Xu, et al.

Fig. 7. Normalized imaging results for the case of two simulated damages: (a) The original imaging result using the DAS method; (b) The original imaging result using the MCDAS method; (c) The filtered image for (a) using a threshold of 0.5; (d) The filtered image for (b) using a threshold of 0.5; (e) The filtered image for (a) using a threshold of 0.75; (f) The filtered image for (b) using a threshold of 0.75. The white dots denote PZT locations and the triangles denote the actual simulated damages locations.

(scattering signals) are decomposed in the first dictionary (matching dictionary), and then the feature signals are extracted using the second dictionary (compression dictionary). The matching dictionary is made up of the excitation and the its time shifted and phase varied signals. Another signal owning a shorter duration, as well as its corresponding time shifted and phase varied signals make up the compression dictionary. The TOF information of the scattering signals is obtained after they are sparsely represented using the matching dictionary. Then the TOF information and the compression dictionary are used to extract the feature signals. The echoes of the feature signals are thinner than the originals, and the overlapped echoes can be separated, therefore the imaging

The result suggests that the presented MC algorithm cannot separate well the overlapped wave packets with strong dispersion effect and multiple modes because the dispersion and multiple modes are not accommodated in dictionary A . So it is necessary to select a suitable center frequency in which a single mode guided wave is dominant and the dispersion effect is relatively weak. 6. Conclusion A combined method of sparse reconstruction and DAS is developed for high resolution Lamb wave inspection in this paper. Under narrowband tone burst excitation, the baseline-subtracted response signals 8

Composite Structures 223 (2019) 110973

C. Xu, et al.

Fig. 8. Comparisons between the original scattering signals of sensor pair 1–2 in the state of two simulated damages to tone bursts with different number of cycles and the corresponding feature signals by the MC algorithm. (a) The original scattering signals to tone bursts with different number of cycles. (b) The corresponding feature signals of (a) by the presented MC algorithm. The notations “#1” and “#2 ” denote the scattering wave packets reflected from damages #1 and #2 , respectively.

Fig. 9. Comparisons between the original scattering signals of sensor pair 1–2 in the state of two simulated damages to excitations with different center frequencies and the corresponding feature signals by the MC algorithm. (a) The original scattering signals to excitations with different center frequencies. (b) The corresponding feature signals of (a) by the presented MC algorithm. The notations “#1” and “#2 ” denote the scattering wave packets from damages #1 and #2 , respectively.

Appendix A. Supplementary data

resolution can be improved. Comparing with the DAS method, the MCDAS method can achieve a smaller spot size and lower noise floor while keeping almost the same localization accuracy. The results suggest that the developed method can achieve a better imaging performance. Future work remains in accommodating anisotropy of the specimen and considering both the dispersion and multiple modes when designing the dictionaries.

Supplementary data associated with this article can be found, in the online version, athttps://doi.org/10.1016/j.compstruct.2019.110973. References [1] Yelve NP, Mitra M, Mujumdar PM. Detection of delamination in composite laminates using lamb wave based nonlinear method. Compos Struct 2017;159:257–66. [2] Ghrib M, Berthe L, Mechbal N, Rébillat M, Guskov M, Ecault R, et al. Generation of controlled delaminations in composites using symmetrical laser shock configuration. Compos Struct 2017;171. [3] Su Z, Ye L, Lu Y. Guided lamb waves for identification of damage in composite structures: a review. J Sound Vib 2006;295(3):753–80. [4] Xu K, Ta D, Moilanen P, Wang W. Mode separation of lamb waves based on dispersion compensation method. J Acoust Soc Am 2012;131(4):2714. [5] Zuo H, Yang Z, Chen X, Xie Y, Miao H. Analysis of laminated composite plates using wavelet finite element method and higher-order plate theory. Compos Struct 2015;131:248–58. [6] Mitra M, Gopalakrishnan S. Guided wave based structural health monitoring: a review. Smart Mater Struct 2016;25(5):053001.

Acknowledgments We thank the supports given by the National Natural Science Foundation of China (Nos. 51875433 & 51605365), the National Key Basic Research Program of China (No. 2015CB057400), the Young Talent Fund of University Association for Science and Technology in Shaanxi of China (No. 20170502), the Fundamental Research Funds for the Central Universities, and the Key Research and Development Program of Shanxi Province (Grant Nos. 2017ZDCXL-GY-02-01 and 2017ZDCXL-GY-02-02). 9

Composite Structures 223 (2019) 110973

C. Xu, et al.

[23] Mallat SG, Zhang Z. Matching pursuits with time-frequency dictionaries. IEEE Trans Signal Process 1993;41(12):3397–415. [24] Hong JC, Sun KH, Kim YY. Waveguide damage detection by the matching pursuit approach employing the dispersion-based chirp functions. IEEE Trans Ultrason Ferroelectr Freq Control 2006;53(3):592–605. [25] Hong JC, Sun KH, Kim YY. The matching pursuit approach based on the modulated gaussian pulse for efficient guided-wave damage inspection. Smart Mater Struct 2005;14(4):548. [26] Zhang J, Drinkwater BW, Wilcox PD, Hunter AJ. Defect detection using ultrasonic arrays: the multi-mode total focusing method. NDT & E Int 2010;43(2):123–33. [27] Demirli R, Saniie J. Model-based estimation pursuit for sparse decomposition of ultrasonic echoes. Signal Process IET 2012;6(4):313–25. [28] Lu Y, Michaels JE. Numerical implementation of matching pursuit for the analysis of complex ultrasonic signals. IEEE Trans Ultrason Ferroelectr Freq Control 2008;55(1):173–82. [29] Zhang GM, Zhang CZ, Harvey DM. Sparse signal representation and its applications in ultrasonic nde. Ultrasonics 2012;52(3):351–63. [30] De Marchi L, Ruzzene M, Xu B, Baravelli E. Warped basis pursuit for damage detection using lamb waves. IEEE Trans Ultrason Ferroelectr Freq Control 2010;57(12):2734–41. [31] Xu C, Yang Z, Chen X, Tian S, Xie Y. A guided wave dispersion compensation method based on compressed sensing. Mech Syst Signal Process 2018;103:89–104. [32] Wang CH, Rose JT, Chang FK. A synthetic time-reversal imaging method for structural health monitoring. Smart Mater Struct 2004;13(2):415–23. [33] Xu B, Giurgiutiu V. Single mode tuning effects on lamb wave time reversal with piezoelectric wafer active sensors for structural health monitoring. J Nondestr Eval 2007;26(2–4):123–34. [34] Hall JS, Michaels JE. Adaptive dispersion compensation for guided wave imaging. Review of progress in quantitative nondestructive evaluation. 2012. p. 623–30. [35] Kexel C, Harley JB, Moll J. Attenuation and phase compensation for guided wave based inspection using a filter approach. Ultrason Symp 2015:1–4. [36] Elad, M. Sparse and redundant representations: from theory to applications in signal and image processing 2010;02(1):1094–1097. [37] Michaels JE, Sang JL, Croxford AJ, Wilcox PD. Chirp excitation of ultrasonic guided waves. Ultrasonics 2013;53(1):265–70. [38] Hall JS, Mckeon P, Satyanarayan L, Michaels JE, Declercq NF, Berthelot YH. Minimum variance guided wave imaging in a quasi-isotropic composite plate. Smart Mater Struct 2011;20(2):025013.

[7] Yang ZB, Chen XF, Xie Y, Zuo H, Miao HH, Zhang XW. Wave motion analysis and modeling of membrane structures using the wavelet finite element method. Appl Math Model 2016;40(3):2407–20. [8] B Z, JG Y, XM Z, PM M. Complex guided waves in functionally graded piezoelectric cylindrical structures with sectorial cross-section. Appl Math Model 2018;63:288–302. [9] Wilcox P, Lowe M, Cawley P. The effect of dispersion on long-range inspection using ultrasonic guided waves. NDT & E Int 2001;34(1):1–9. [10] Park HW, Kim SB, Sohn H. Understanding a time reversal process in lamb wave propagation. Wave Motion 2009;46(7):451–67. [11] Liu L, Yuan FG. A linear mapping technique for dispersion removal of lamb waves. Struct Health Monit 2010;9(1):75–86. [12] Zhou C, Su Z, Cheng L. Quantitative evaluation of orientation-specific damage using elastic waves and probability-based diagnostic imaging. Mech Syst Signal Process 2011;25(6):2135–56. [13] Li X, Yang Z, Zhang H, Du Z, Chen X. Crack growth sparse pursuit for wind turbine blade. Smart Mater Struct 2015;24(1):15002. 15009(8). [14] Hao Z, Zhibo Y, Caibin X, Shaohua T, Xuefeng C. Damage identification for platelike structures using ultrasonic guided wave based on improved music method. Compos Struct 2018;203:164–71. [15] Yang Z, Chen X, Li X, Jiang Y, Miao H, He Z. Wave motion analysis in arch structures via wavelet finite element method. J Sound Vib 2014;333(2):446–69. [16] Lin J, Hua J, Zeng L, Luo Z. Excitation waveform design for lamb wave pulse compression. IEEE Trans Ultrason Ferroelectr Freq Control 2015;63(1):165–77. [17] Donoho DL. For most large underdetermined systems of linear equations the minimal 11-norm solution is also the sparsest solution. Commun Pure Appl Math 2010;59(6):797–829. [18] Donoho DL. Compressed sensing. IEEE Trans Inf Theory 2006;52(4):1289–306. [19] Levine RM, Michaels JE. Block-sparse reconstruction and imaging for lamb wave structural health monitoring. IEEE Trans Ultrason Ferroelectr Freq Control 2014;61(6):1006–15. [20] Mesnil O, Ruzzene M. Sparse wavefield reconstruction and source detection using compressed sensing. Ultrasonics 2016;67:94–104. [21] Pati YC, Rezaiifar R, Krishnaprasad PS. Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition. Signals, Systems and Computers, 1993. 1993 Conference Record of The Twenty-Seventh Asilomar Conference on. 2002. p. 40–4. vol. 1. [22] Needell D, Tropp JA. CoSaMP: iterative signal recovery from incomplete and inaccurate samples. ACM; 2010.

10