Singular spectrum analysis and its application in Lamb wave-based damage detection

2.1. Embedding

A trajectory matrix or a lagged version is obtained during this step. Let $y=[y_1,y_2,\dots ,y_N]$ be a time series of length $N$, $L$ is an integer (window length) and $K$ is the number of lagged vector ($K=N-L+1$). The obtained trajectory matrix is defined as:

Eq. (1) shows that the trajectory matrix is a Hankel matrix since all the elements along the ‘anti-diagonals’ are the same.

2.2. Singular value decomposition (SVD)

In this step, the SVD of the matrix $\mathbf{Y}$ is obtained. It is further decomposed into a sum of mutually orthogonal, unit-rank, elementary matrices. Let $\mathbf{S}=\mathbf{Y}{\mathbf{Y}}^T$ be a $L$×$L$ matrix. Let ${\lambda }_1$, ${\lambda }_2$,…, ${\lambda }_r$ be the non-zero eigenvalues of $\mathbf{S}$ ($r=L$ if no eigenvalue is zero) arranged in decreasing order, and $U_1$, $U_2$,…, $U_r$ be the corresponding eigenvectors. The vectors $V_j={\mathbf{Y}}^T(U_j/\sqrt{{\lambda }_j})$ for $j=$ 1, 2,…, $r$ is constructed.

The SVD of the trajectory matrix $\mathbf{Y}$ can be written as:

where ${\mathbf{Y}}_j=\sqrt{{\lambda }_j}{U}_j{V}_j^T$.

The obtained $L$ numbers of singular value are the square roots of the eigenvalues of matrix $\mathbf{S}$. The plot of the singular value in decreasing order is known as singular spectrum, and thus gives the method its name.

2.3. Diagonal averaging

Each matrix of Eq. (2) is transformed into a new time series of length $N$ by applying a linear transformation known as diagonal averaging. The diagonal averaging algorithm is as follows:

Since $\mathbf{Y}$ is an $L$×$K$ matrix with elements $y_{ij}$, 1 $\le i\le L$, 1 $\le j\le K$. Set $L^{\mathrm{*}}=\min (L,K)$, $K^{\mathrm{*}}=\max (L,K)$ and $N=L+K-1$. Let $y_{ij}^{\mathrm{*}}=y_{ij}$ if $L<K$ and $y_{ij}^{\mathrm{*}}=y_{ji}$ otherwise. By averaging the diagonal the matrix $\mathbf{Y}$ is transferred into the series $y_1$,…, $y_N$ using the formula:

The diagonal averaging procedure is applied to all the selected matrix components $({\mathbf{Y}}_1$, ${\mathbf{Y}}_2$,…, ${\mathbf{Y}}_r)$ to form Hankel matrices $({\stackrel{-}{\mathbf{Y}}}_1$, ${\stackrel{-}{\mathbf{Y}}}_2$,…, ${\stackrel{-}{\mathbf{Y}}}_r)$. Each Hankel matrix determines a time series that is recognized as the projection of the original signal in the corresponding direction determined by the corresponding eigentriple. Then the sum of all selected projected time series component represents the reconstructed signal.

2.4. Parameter selection

It should be mentioned that two parameters are required in SSA: the window length $L$ and the number of reconstructed components $r$. A detailed description of parameter selection is presented in Golyandina et al. [9]. Values for $L$ and $r$ could be defined using information provided by the time series under study or through additional indices. The general rule is to select $L=N/2$, which is also the value used in this paper. The dimension of the selected components, which is defined as $r$, is very important. In fact the “selection” of useful signals with SSA is still under studied with no unified approach. However, by the data structure of Hankel matrix it can be seen that the next row vector is just lagging one data point compare to the front row vector. Therefore the adjacent two lines of the noise-free Hankel matrix (in a morbid form) will be highly correlated, and the rank $r$ is much smaller than min ($m$, $n$). The number of non-zero singular value is equal to the rank of the matrix according to the theory of SVD, which means only the first $r$ singular values are non-zero, while the rest are substantially zero. And there will be a curve exists at the turning point.

For the noise constructed Hankel matrix, the autocorrelation function of noise sequence $\zeta \left(i\right)$ is as follows:

where $b$ is standard deviation of $\zeta \left(i\right)$, $\delta \left(\tau \right)$ is Kronecker function:

Eq. (4) and Eq. (5) show that there is no correlation between the two adjacent rows of the noise constructed Hankel matrix, the matrix is benign with rank $r=\mathrm{m}\mathrm{i}\mathrm{n}(m,n)$. There will be no turning point exists in the singular value curve for the noise-only condition. Thus, the turning point of the singular value curve can be the distinction of real signal and noise involved signal. At this turning point energy of the real signal is more concentrated, while the noise energy is relatively dispersed. It can be concluded that there will be a “peak” (not unique) at the boundary between signal and noise, and noise makes main contribution after that “peak”. Singular value differential spectrum value is then defined as:

By Eq. (6) we can select the critical point $r$ (also not unique) to be the reconstructed components.

3.1. Modeling Lamb wave for damage detection

This subsection is mainly about the introduction of active Lamb wave in structure damage detection. Lamb waves are firstly inspired by piezoelectric ceramics and narrow band excitation is applied here which has a form as bellows:

where $n$ is the peak number, $f_c$ is the center frequency, $A$ is the amplitude modulation and $H$ is the Heaviside step function.

The inspired field quantity is expressed as $u(x,t)$, where $x$ presents the propagation distance and $t$ presents the propagation time. If one point of $u(x,t)$ and the propagation characteristics of Lamb waves are known, then any point of $u(x,t)$ may be obtained by the phase shift of each frequency component. And $u(x,t)$ can be expressed as [20]:

where $k\left(\omega \right)$ is the wave number, $\omega$ is the angular frequency and $F\left(\omega \right)$ is the Fourier transform of $f\left(t\right)$.

During the propagation the inspired signal will encounter reflector (such as damage spot or boundary) and the reflected signal $g\left(t\right)$ is typically superimposed by different reflected sources. The reflected signal to the sensor which is going across the $j$th reflected source can be expressed as:

where $d_j$ is the distance between actuator and reflector, $x_j$ is the distance between reflector and sensor, $A{}_j{}^{\mathrm{}}\left(\omega \right)$ is the reflection coefficient.

The propagation distance from the actuator to the sensor is defined as $l_j$ and:

Then Eq. (9) can be expressed as:

Fig. 1Schematic diagram of detection principle

Eq. (9) and Eq. (11) show that signal from the $j$th reflected source $g_j\left(t\right)$ can be considered as the propagation signal from the actuator to another sensor (see Fig. 1).

3.2. Damage detection based on the maximum similarity analysis

In this subsection the maximum similarity analysis method is introduced. This method is based on the priori knowledge as below:

1. Dispersion characteristics of the propagated Lamb waves are obtained from the material properties;

2. Direct wave arrival time is known when structural model dimensions and actuator/sensor arrangement are determined;

3. Reflected signal is typically superimposed by different reflected sources and it can be determined by the form of the inspired signal and propagation distance.

It is assumed that there are totally $n$ reflected sources and all of them (damage scattered waves or boundary reflected waves) can be determined by Eq. (11), the synthetic signal that consists of damage scattered waves and boundary reflected waves is further obtained as:

where $a_n$ is the control factor of the amplitude.

Then the extraction of damage scattered signal can be converted to the optimization problem: Parameter optimization ($a_n$, $l_n$ and $n$) that base on the Genetic Algorithms (GA). Cosine method [21] is applied to be the similarity index to measure the degree of coincidence between the synthetic signal and the target signal, the formula is as follows:

where $x$ presents the synthetic signal and $y$ presents the target signal. The fitness function is as follows:

The two signals have a better similarity when fitness function value is closer to zero. Fig. 2 shows the main process of the maximum similarity analysis that based on GA.

Fig. 2Block diagram

3.3. Damage location

Time of Flight [22] method is used for damage location. The actuator and sensor can be regarded as two focus of an ellipse, which means the ellipse can be determined. However, damage can’t be located if only one ellipse is known. So another sensor is added for obtaining a second ellipse and the intersection of the two ellipses is the exact damage location.

3.4. Experiment analysis

Test device and specimen are shown in Fig. 3 and Fig. 4, including Sonox P5 circular piezoelectric plate (diameter is 6 mm), RIGOL DG1022 Arbitrary Waveform Generator, HVPA05 amplifier, Tektronix TDS2014C oscilloscopes and 2024-T3 aluminum specimen (dimension is 600 mm×600 mm×1 mm). The boundary-free aluminum plate is supported by foam and the material parameters are shown in Table 1. Three piezoelectric sheets are pasted with conductive silver colloid. PZT A is actuator, PZT B and PZT C is sensor. A transparent hole (diameter is 5 mm) is represented for the damage, five-cycle sine signal (modulated with Hanning window) is adopted as the inspired signal, with the center frequency of 300 KHz and the inspired voltage peak is 20 V. In addition, TDS2014C digital oscilloscope is set to 256 times average automatically processing in order to improve the stability and signal to noise ratio, and trigger interval is set to 1 second to ensure that the last inspiration and reflection disappear completely. Equipment and each piezoelectric sheet are accessed to the ground to avoid electromagnetic interference wave packet.

Fig. 3Test device

Fig. 4Test specimen

Original response signal of path AB is shown in Fig. 5. SSA technique is applied for the original signal and the singular value differential spectrum is shown in Fig. 6. In Fig. 6 “peak” exist when $r$ is equal to 2, 4, 6, 8, 9, 11, 13 and 19. The reconstructed results of these components are plotted in Fig. 6 respectively. The first 19 components reconstructed result is a very accurate approximation to the real response on comparison of Fig. 5 and Fig. 7. Since no aliasing happens in direct wave so the target signal can be simplified as the dashed frame part of Fig. 7(h).

Fig. 5Damage signal captured by PZT B

Fig. 6Singular value differential spectrum

Maximum similarity analysis that based on GA is applied for the optimization of the target signal and the optimization results are shown in Table 2. Here $n$ presents the number of reflected waves (damage scattered wave or boundary reflected waves). The accurate propagation distance $l$ and propagation time $t$ of the boundary reflected wave can be obtained by the position of the sensor and the boundary (see Table 3).

Fig. 7Reconstructed responses with different combinations of components

a) First 2 components

b) First 4 components

c) First 6 components

d) First 8 components

e) First 9 components

f) First 11 components

g) First 13 components

h) First 19 components

Table 3 shows the accurate propagation distance $l$ is equal to 0.415, and 0.415 also appears in Table 3 when $n=$ 2 and $n=$ 3. Since there will be only one damage scattered wave and the rest are boundary reflected waves, optimization results match the known conditions best when $n=$ 2. So the synthetic signal contains two waves: damage scattered wave and upper boundary reflected wave. The time domain waveform is shown in Fig. 8. Fig. 9 shows the comparison of synthetic signal and target signal. In addition, Fig. 4 shows that the accurate propagation of damage scattered wave $l_{accurate}$ is equal to 0.3795 m, while the optimization results fit it very well ($l_{calculate}=$0.380, see Table 2).

Fig. 8Reflected signal

a) Boundary reflected wave

b) Damage scattered wave

Fig. 9Comparison of synthetic signal and target signal

In order to obtain the exact damage location another path AC is analyzed. Fig. 10 shows the original response signal of path AC. SSA technique is also applied for the original signal and the singular value differential spectrum is shown in Fig. 11. “Peak” exists when $r$ is equal to 2, 4, 6, 8, 10, 11, 13, 15, 17, 19, 21 and 24. The reconstructed results of these components are plotted in Fig. 12 respectively. The first 21 components reconstructed result is a very accurate approximation to the real response on comparison of Fig. 10 and Fig. 12. The same optimization analysis is conducted and Table 4 shows the damage localization results. The identification results show that the damage location can be identified accurately with small positive errors when the best reconstructed components are selected. Fig. 13 shows the comparison of the actual position and the ellipse orientation position.

Fig. 10Damage signal captured by PZT C

Fig. 11Singular value differential spectrum

Fig. 12Reconstructed responses with different combinations of components

a) First 2 components

b) First 4 components

c) First 6 components

d) First 8 components

e) First 10 components

f) First 11 components

g) First 13 components

h) First 15 components

i) First 17 components

j) First 19 components

k) First 21 components

l) First 24 components