A Local TR-MUSIC Algorithm for Damage Imaging of Aircraft Structures

Lamb wave-based damage imaging is a promising technique for aircraft structural health monitoring, as enhancing the resolution of damage detection is a persistent challenge. In this paper, a damage imaging technique based on the Time Reversal-MUltiple SIgnal Classification (TR-MUSIC) algorithm is developed to detect damage in plate-type structures. In the TR-MUSIC algorithm, a transfer matrix is first established by exciting and sensing signals. A TR operator is constructed for eigenvalue decomposition to divide the data space into signal and noise subspaces. The structural space spectrum of the algorithm is calculated based on the orthogonality of the two subspaces. A local TR-MUSIC algorithm is proposed to enhance the image quality of multiple damages by using a moving time window to establish the local space spectrum at different times or different distances. The multidamage detection capability of the proposed enhanced TR-MUSIC algorithm is verified by simulations and experiments. The results reveal that the local TR-MUSIC algorithm can not only effectively detect multiple damages in plate-type structures with good image quality but also has a superresolution ability for detecting damage with distances smaller than half the wavelength.


Introduction
Structural health monitoring (SHM) has been considered a revolutionary technology for the operation and maintenance of aircraft structures [1]. At present, some SHM technologies, such as acoustic emission [2,3], eddy current [4][5][6], Lamb wave [7,8], impedance [9,10], vibration [11,12], and comparative vacuum [13] monitoring, have received greatly focused research and development for damage detection in aircraft structures. Lamb wave-based damage imaging, which is an active detection technology suitable for large-scale structures, is a promising technology for the health monitoring of aircraft plate structures [14,15].
However, due to the multiple modes, wave dispersion and boundary scattering of Lamb waves, the acquired Lamb wave signal is often difficult to interpret and needs to be further analyzed. With the help of advanced signal and information processing algorithms, Lamb wave damage imaging technology, which can intuitively provide the location and size of the damage, has become a topic of interest in SHM technology. Existing Lamb wave damage imaging technologies, such as probabilistic reconstruction algorithms [16,17], phased arrays [18,19], tomography [20,21] and delay-and-sum [22,23] methods, can realize damage location and even quantification in plate-type structures. However, due to the limitations of wave diffraction, existing damage imaging algorithms have poor resolution and difficulty distinguishing two neighboring damage within the Rayleigh limit. The deficiency of damage monitoring capability makes it impossible to acquire sufficient information to predict the residual life of structures, which limits the further development of Lamb wave-based SHM technology.
Multiple signal classification (MUSIC) is a typical method used to address the Rayleigh criterion. Schmidt et al. [24] proposed the MUSIC algorithm based on eigenvalue decomposition and subspace theory. The algorithm offers high estimation resolution and stability for estimating the direction of a target and has been widely used in electromagnetic waves, acoustic waves and other fields. For wave-based damage identification, each damage can be regarded as a secondary wave source according to the Huygens principle, which means that the MUSIC algorithm can be used to estimate the damage location. Engholm et al. [25] used the MUSIC algorithm to estimate the direction of the damage source but did not provide the damage distance. Ambrozinski et al. [26] applied the MUSIC algorithm for the damage location of an aluminum plate. Zhong et al. [27][28][29] developed a series of improved MUSIC algorithm to detect impact sources under varying temperature, vibration and deformation conditions. Yuan et al. [19,30,31] made use of the MUSIC algorithm for damage location estimation of aircraft structures.
Time reversal multiple signal classification (TR-MUSIC), which combines the adaptive focusing of the time reversal algorithm and the superresolution of direction estimation of the MUSIC algorithm, has also been used in ultrasonic nondestructive testing. Devaney et al. [32] first proposed TR-MUSIC and verified its superresolution ability in identifying two parallel copper wires with a short distance in water by numerical simulations. Simonetti [33] used ultrasonic TR-MUSIC to create defect images for the first time and successfully distinguished two holes with half wavelengths in a steel block, which broke through the limitation of acoustic diffraction and realized the ultrasonic superresolution imaging of solid media. Fan et al. [34] verified the excellent performance of TR-MUSIC in superresolution damage imaging by comparing TR-MUSIC with the full focus method (TFM). Several studies assessed TR-MUSIC for the damage detection of plate-type structures by using ultrasonic Lamb waves. He and Yuan [35] imaged multiple damages over a long distance in an aluminum plate by using a piezoelectric array and the DORT-MUSIC algorithm. However, their study did not discuss the superresolution capability of detecting multiple damages. In addition, two other problems exist in multidamage detection. One is image interference caused by pixel values of different orders of magnitude at different damage locations, which results from different damage sizes and different distances from the damage to the sensing array. Another is that when the amount of damage is large and even larger than the sensor numbers, the TR-MUSIC algorithm is prone to generating an incorrect damage image.
In this paper, an enhanced TR-MUSIC algorithm is proposed to study the capability of detecting multiple damages in plate-type structures. The superresolution ability to detect multiple damages, whose distances are smaller than a half-wavelength, is discussed. In addition, the enhanced algorithm can markedly improve the quality of multidamage images.
The structure of this paper is established as follows. In Section 2, the process and formulations of TR-MUSIC are introduced. In Section 3, the ability of the local TR-MUSIC algorithm is verified for multidamage detection by using the signals from the finite element model. The algorithm is further verified by experiments in Section 4. Finally, conclusions are drawn in Section 5.

Local TR-MUSIC Algorithm
The basic process and theoretical formulations have been introduced in the references [29][30][31][32]. In this paper, the formulations of the TR-MUSIC algorithm based on a linear piezoelectric sensor array and Lamb waves can be written as follows. Figure 1 shows an operational diagram of a piezoelectric sensor array. A piezoelectric array containing N elements is bonded to the surface of a plate, and piezoelectric elements are used not only as transmitters but also as receivers. Each element in the piezoelectric array is excited in turn, and all the array elements receive response signals. Assume that T m (m = 1, 2, · · · N) represents the transmitter location, R n (n = 1, 2, · · · , N) represents the receiver location, and D q , q = 1, 2 · · · , Q represents the damage location.  Figure 1 shows an operational diagram of a piezoelectric sensor array. A piezoelectric array containing N elements is bonded to the surface of a plate, and piezoelectric elements are used not only as transmitters but also as receivers. Each element in the piezoelectric array is excited in turn, and all the array elements receive response signals.  The mth 1 m N   ( ) piezoelectric element in the piezoelectric array is used as a transmitter to excite the Lamb wave propagating in the structure. The wave at location s r away from the transmitter can be expressed as

Transfer Matrix
is Green's function, which represents the process of the Lamb wave propagating from m T to s r . In the frequency domain, Equation (1) can be written as s m g t r T , respectively. The frequency symbol  will be omitted in the following mathematical expression.
After scattering by the damage at s r , the Lamb wave is acquired by the nth receiver n R as follows: For point-type damage, the influence of different damage can be neglected, and then the scattering coefficient of the structure can be expressed as follows: The mth (1 ≤ m ≤ N) piezoelectric element in the piezoelectric array is used as a transmitter to excite the Lamb wave propagating in the structure. The wave at location r s away from the transmitter can be expressed as where g(r s , T m , t) is Green's function, which represents the process of the Lamb wave propagating from T m to r s . In the frequency domain, Equation (1) can be written as where X(ω), Y sm (ω) and G(r s , T m , ω) are Fourier transforms of x(t), y sm (t) and g(r s , T m , t), respectively. The frequency symbol ω will be omitted in the following mathematical expression. After scattering by the damage at r s , the Lamb wave is acquired by the nth receiver R n as follows: where O(r s ) is the scattering coefficient of the damage. For point-type damage, the influence of different damage can be neglected, and then the scattering coefficient of the structure can be expressed as follows: where σ q (r q ) is the scattering coefficient of the damage at r q and δ(r − r q ) is the Dirac function. Substituting Equation (4) into Equation (3), the signal that is sensed by the nth receiver and transmitted by the mth transmitter can be written as Then, the transfer matrix element can be defined as where the subscript 'nm' indicates that the mth sensor is used as the transmitter and the nth sensor is used as the receiver. Therefore, the transfer matrix can be defined as where R, r and T represent the locations of receivers, damage and transmitters, respectively.
G(R, r) and G T (r, T) are Green's function matrices, which are associated sets of Green's functions G(R n , r q ) and G(r q , T m ), respectively.

Time-Reversal Operation
The TR process can be conducted from the transfer function matrix.
If X (0) and Y (0) are the input and output matrices, respectively, the following equation can be obtained: The time reversal of the output signal is equal to the phase conjugation of the frequency-domain signal, and it yields where the symbol '*' indicates the conjugate operation on the matrix and the superscript '1' indicates the first time-reversal operation. X (1) is considered to be the new input signal and is emitted, and the new output matrix can be written as Y (1) is reversed again as the new input matrix After the 2n iterations, the input matrix of the system is written as After the 2n + 1 iterations, the new input matrix X (2n+1) is From Equations (12) and (13), two kinds of time reversal operators are defined as where T 1 and T 2 are the time reversal operators for the transmitting and receiving modes, respectively. If the transfer matrix is a complex symmetric matrix, the following feature of the TR operation can be obtained:

Singular Value Decompositions
The singular value decompositions of the transfer matrix and its conjugate transposition can be expressed as where σ i is defined as the singular value of the transfer matrix, while µ i and ν i are the left and right singular vectors related to σ i , respectively. The eigenvalue decompositions of T 1 and T 2 can be expressed as follows: The eigenvalue decomposition of the TR operator is used to divide the whole data space into a signal subspace and noise subspace according to the distribution of eigenvalues. Under the assumption of ideal point damage, there is a one-to-one mapping relationship between each damage on the structure and the eigenvector corresponding to each nonzero eigenvalue of the TR operator. Therefore, the vector space composed of eigenvectors corresponding to nonzero eigenvalues is usually regarded as the signal subspace, while the vector space composed of eigenvectors corresponding to zero eigenvalues is regarded as the noise subspace.

TR-MUSIC Imaging Algorithm
The basis of the signal subspace is obtained by eigenvalue decomposition of the TR operator. Under the Born approximation, the signal subspace can also be calculated by Green's function. The signal subspace and noise subspace are orthogonal to each other so that the inner product of the noise subspace and the Green function vector at the real damage site should be equal to 0. The imaging function of TR-MUSIC is defined as where Green's function at each discrete point of the structure is used to calculate the inner product with the noise subspace to construct the spatial spectrum of the TR-MUSIC algorithm. The peak of the spatial spectrum is the location of the real damage.

Local TR-MUSIC Imaging Algorithm
TR-MUSIC has a poor capability of detecting multiple damages, especially when the element number of the sensor array is small. Therefore, to solve the above problem, an enhanced monitoring strategy based on the local TR-MUSIC algorithm has been developed. The process of the local TR-MUSIC algorithm is similar to the short-time Fourier transform. First, a Hanning window moving along the time axis is used to separate the original scattered signal y sm (t) into a set of short-segment signals: where w(t − t i ) denotes the Hanning window function at the ith time point and I represents that the signal is separated into I segments. Different Is denote different step sizes in which the window moves, and the larger the value of I is, the smaller the step size of the window function is. The local segments of scatter signals z i (t) are used to calculate the local transfer matrix and perform the local eigenvalue decomposition. Thus, there exists a timedependent eigenvalue function. Generally, it can be imagined that the signal segments that are scattered from the real damage will generate a local maximum point in the eigenvalue function. When the damages have different distances from the sensor array, the time of scattering signals caused by each damage is different, and the time of the local maximum point in the eigenvalue function is also different. This process can be expressed as where κ(•) denotes the transfer function, which is the quotient of the scattering signal and incident signal in the frequency domain. (•) and µ(•) denote the process of calculating the TR operator and the eigenvalue, respectively, whose formulations can be seen in Equations (15) and (16). For different i, µ( (κ(z i (t)))) makes a time-dependent eigenvalue function. Thus, 'argmax' means looking for the time when the eigenvalue function takes the local maximum. t is just a set of time points where local maximum values exist. Thus, for short-segment signals z i (t) at each time point in t, a damage image performed by Equation (22) can locate the real damage. The multidamage contour can be summed by the damage images at all time points in t, which is written by This can effectively avoid the occurrence of missing damage. In addition, since the scattered signals of different time periods are cut out for damage location processing from zero, this is equivalent to scanning the structure from near to far according to the distance from each structural point to the sensing array. At the same distance, the probability that the damage number exceeds the total number of sensing array elements is almost negligible. Therefore, this method can also solve the problem whereby the damage number identified by the TR-MUSIC algorithm cannot exceed the total number of elements in the sensing array.

Simulation Verification
A finite element simulation based on the ABAQUS Standard is used to verify the effectiveness and feasibility of the developed damage imaging algorithm. In Figure 2, an aluminum plate with dimensions of 600 mm × 600 mm × 1 mm, an elastic modulus of 70 GPa, a Poisson's ratio of 0.3, and a mass density of 2700 kg/m 3 is considered. The boundary of the plate is fixed all around. The element type is C3D8R, and the element size is 1 mm. There are a total of 360,000 elements. The time interval is set as 0.1 µs, and the recording time is 0.5 ms. The exciting signal, with a five-peaked Hanning windowed sinusoidal tone burst and a center frequency of 50 kHz, is loaded along the in-plane direction at array nodes to model the piezoelectric sensor. Taking the center point as the coordinate origin, a linear sensor array with 11 array elements spaced at 5 mm lies at the left of and 150 mm away from the coordinate origin. of the plate is fixed all around. The element type is C3D8R, and the element size is 1 mm. There are a total of 360,000 elements. The time interval is set as 0.1 μs, and the recording time is 0.5 ms. The exciting signal, with a five-peaked Hanning windowed sinusoidal tone burst and a center frequency of 50 kHz, is loaded along the in-plane direction at array nodes to model the piezoelectric sensor. Taking the center point as the coordinate origin, a linear sensor array with 11 array elements spaced at 5 mm lies at the left of and 150 mm away from the coordinate origin.

Single-and Multi-Damage Identification
In this section, two damage cases-a single damage and three damages-are considered. Based on the coordinates in Figure 2, a few damage locations are set casually, and detailed damage information is given in Table 1.  Figure 3a shows the baseline and current signals recorded by the 6th sensor. The signals before and after damage show almost no change. Then, the scattered signal, with a smaller magnitude than the incident wave in Figure 3a, is obtained by subtracting the baseline signal from the current signal, as shown in Figure 3b. The first wave packet in the

Single-and Multi-Damage Identification
In this section, two damage cases-a single damage and three damages-are considered. Based on the coordinates in Figure 2, a few damage locations are set casually, and detailed damage information is given in Table 1.  Figure 3a shows the baseline and current signals recorded by the 6th sensor. The signals before and after damage show almost no change. Then, the scattered signal, with a smaller magnitude than the incident wave in Figure 3a, is obtained by subtracting the baseline signal from the current signal, as shown in Figure 3b. The first wave packet in the signal is the scattering signal caused by damage, and the subsequent wave packet is the scattered signal from the boundary. The scattering signal caused by the damage is extracted to calculate the transfer matrix and the TR operator, and eigenvalue decomposition of the TR operator is performed. location of the spatial spectrum coincides with the real damage location of the structure, which means that the TR-MUSIC algorithm successfully realizes the damage image of a single damage. Figures 6 and 7 show the scattered signals, eigenvalue distribution and damage imaging results of three damages. Multiple damage imaging is also successfully realized by TR-MUSIC.     Figure 4 shows the normalized eigenvalue distribution of the TR operator. There is only one nonzero eigenvalue, and the other eigenvalues are all zeros. The signal subspace and the noise subspace can be obtained. The spatial spectrum of the TR-MUSIC algorithm is constructed to identify the damage. Figure 5 shows the imaging results of Case I in Table 1. The symbol 'X' indicates the location of the real damage in the structure. The peak-value location of the spatial spectrum coincides with the real damage location of the structure, which means that the TR-MUSIC algorithm successfully realizes the damage image of a single damage. Figures 6 and 7 show the scattered signals, eigenvalue distribution and damage imaging results of three damages. Multiple damage imaging is also successfully realized by TR-MUSIC.
signal is the scattering signal caused by damage, and the subsequent wave packet is the scattered signal from the boundary. The scattering signal caused by the damage is extracted to calculate the transfer matrix and the TR operator, and eigenvalue decomposition of the TR operator is performed. Figure 4 shows the normalized eigenvalue distribution of the TR operator. There is only one nonzero eigenvalue, and the other eigenvalues are all zeros. The signal subspace and the noise subspace can be obtained. The spatial spectrum of the TR-MUSIC algorithm is constructed to identify the damage. Figure 5 shows the imaging results of Case I in Table  1. The symbol 'X' indicates the location of the real damage in the structure. The peak-value location of the spatial spectrum coincides with the real damage location of the structure, which means that the TR-MUSIC algorithm successfully realizes the damage image of a single damage. Figures 6 and 7 show the scattered signals, eigenvalue distribution and damage imaging results of three damages. Multiple damage imaging is also successfully realized by TR-MUSIC.

Verification of Local TR-MUSIC Algorithm
The same model as shown in Figure 2 is considered. Three damages are set in the model. Damage sizes are all 4 mm × 4 mm × 1 mm. Damage coordinates are (−50 mm, −50 mm), (0 mm, 0 mm) and (50 mm, 50 mm). Figure 8a shows the scattered signal received by all sensing elements. Figure 8b shows the normalized eigenvalue distribution. There are three nonzero eigenvalues, which are consistent with the number of damage in the structure. Therefore, the vector space composed of the eigenvectors corresponding to the three nonzero eigenvalues is taken as the signal subspace, while the vector space composed of the other eigenvectors corresponding to the zero eigenvalues is used as the noise subspace for damage identification.

Verification of Local TR-MUSIC Algorithm
The same model as shown in Figure 2 is considered. Three damages are set in the model. Damage sizes are all 4 mm × 4 mm × 1 mm. Damage coordinates are (−50 mm, −50 mm), (0 mm, 0 mm) and (50 mm, 50 mm). Figure 8a shows the scattered signal received by all sensing elements. Figure 8b shows the normalized eigenvalue distribution. There are three nonzero eigenvalues, which are consistent with the number of damage in the structure. Therefore, the vector space composed of the eigenvectors corresponding to the three nonzero eigenvalues is taken as the signal subspace, while the vector space composed of the other eigenvectors corresponding to the zero eigenvalues is used as the noise subspace for damage identification.
(a) Scattered signal (b) Normalized eigenvalue distribution  Figure 9 shows the imaging results of the TR-MUSIC algorithm for three damage. Comparisons between the identified damage locations, which are obtained from the local peak locations in the image, and the real damage locations can be seen in Table 2. Figure  9 and Table 2 show that only the damage at the location (0 mm, 0 mm) has a high imaging accuracy. For the damage at the location (50 mm, 50 mm), although the damage direction identification is correct, the identification error of radial distance is very large. Moreover, the damage at the location (−50 mm, −50 mm) is not even visible in the damage image.   Figure 9 shows the imaging results of the TR-MUSIC algorithm for three damage. Comparisons between the identified damage locations, which are obtained from the local peak locations in the image, and the real damage locations can be seen in Table 2. Figure 9 and Table 2 show that only the damage at the location (0 mm, 0 mm) has a high imaging accuracy. For the damage at the location (50 mm, 50 mm), although the damage direction identification is correct, the identification error of radial distance is very large. Moreover, the damage at the location (−50 mm, −50 mm) is not even visible in the damage image. To solve this problem, the local TR-MUSIC algorithm is applied. The Hanning window length is one-fifth of the length of the exciting signal. Figure 10 shows the calculated time-dependent eigenvalue function. There are three extreme points in the function, corresponding to 0.222 ms, 0.282 ms and 0.368 ms. Therefore, the scattering signals at these  To solve this problem, the local TR-MUSIC algorithm is applied. The Hanning window length is one-fifth of the length of the exciting signal. Figure 10 shows the calculated time-dependent eigenvalue function. There are three extreme points in the function, corresponding to 0.222 ms, 0.282 ms and 0.368 ms. Therefore, the scattering signals at these three moments are imaged. To solve this problem, the local TR-MUSIC algorithm is applied. The Hanning window length is one-fifth of the length of the exciting signal. Figure 10 shows the calculated time-dependent eigenvalue function. There are three extreme points in the function, corresponding to 0.222 ms, 0.282 ms and 0.368 ms. Therefore, the scattering signals at these three moments are imaged.    Figure 11d shows the combined imaging results. The local TR-MUSIC algorithm is capable of identifying multiple damages, and it can realize damage identification from near to far according to the distance from each structural point to the sensing array. Even when the damage number is larger than the array element number, the local TR-MUSIC can provide good damage identification ability.

Super-Resolution Imaging
In this section, the superresolution capability of the TR-MUSIC algorithm is verified. By dispersion analysis, the wavelength of the A0 wave at the frequency-thickness product of 50 kHz·mm is 13.7 mm. In the same model, as shown in Figure 2, two damage with a size of 4 mm × 4 mm × 1 mm are set at (0 mm, −3 mm) and (0 mm, 3 mm), respectively. Thus, the central distance between these two damages is 6 mm, and the shortest distance is only 2 mm. the distance from each structural point to the sensing array. Even when the damage number is larger than the array element number, the local TR-MUSIC can provide good damage identification ability.

Super-Resolution Imaging
In this section, the superresolution capability of the TR-MUSIC algorithm is verified. By dispersion analysis, the wavelength of the A0 wave at the frequency-thickness product of 50 kHz•mm is 13.7 mm. In the same model, as shown in Figure 2, two damage with a size of 4 mm × 4 mm × 1 mm are set at (0 mm, −3 mm) and (0 mm, 3 mm), respectively. Thus, the central distance between these two damages is 6 mm, and the shortest distance is only 2 mm. Figure 12a shows the scattered signal. There is only one scattering wave packet in the signal, so it is impossible to distinguish the two damage types from the time domain signal. The scattered signal caused by the damage is substituted into the TR-MUSIC algorithm. The eigenvalues are shown in Table 3, and it can be seen that eigenvalues No. 1 and No. 2 are quite larger than the other eigenvalues. Therefore, the vector space composed of the eigenvectors corresponding to the first two eigenvalues is taken as the signal subspace, and the vector space composed of the eigenvectors corresponding to the other eigenvalues  Figure 12a shows the scattered signal. There is only one scattering wave packet in the signal, so it is impossible to distinguish the two damage types from the time domain signal. The scattered signal caused by the damage is substituted into the TR-MUSIC algorithm. The eigenvalues are shown in Table 3, and it can be seen that eigenvalues No. 1 and No. 2 are quite larger than the other eigenvalues. Therefore, the vector space composed of the eigenvectors corresponding to the first two eigenvalues is taken as the signal subspace, and the vector space composed of the eigenvectors corresponding to the other eigenvalues is taken as the noise subspace. Based on the orthogonality of the Green function and noise subspace, the spatial spectrum of the TR-MUSIC algorithm is constructed. Figure 13 shows the TR-MUSIC imaging results. It can be seen that the superresolution imaging of damage with a less than half-wavelength distance has been successfully realized.
is taken as the noise subspace. Based on the orthogonality of the Green function and noise subspace, the spatial spectrum of the TR-MUSIC algorithm is constructed. Figure 13 shows the TR-MUSIC imaging results. It can be seen that the superresolution imaging of damage with a less than half-wavelength distance has been successfully realized.    is taken as the noise subspace. Based on the orthogonality of the Green function and noise subspace, the spatial spectrum of the TR-MUSIC algorithm is constructed. Figure 13 shows the TR-MUSIC imaging results. It can be seen that the superresolution imaging of damage with a less than half-wavelength distance has been successfully realized.

Enhancing Radial Resolution Using Dual Arrays
It can be seen from the previous sections that the proposed TR-MUSIC algorithm has a good angle resolution for damage detection but a poor radial resolution. In order to eliminate this problem, dual arrays having an 'L' shape can be used for damage imaging. As shown in Figure 14, one horizontal and one vertical array containing 6 sensors are used to detect two damage in the aluminate panel. The same process is used for both horizontal and vertical arrays. Figure 15 gives comparisons of damage images between single array and dual arrays. It is obvious that although a single array has also a poor radial resolution, dual arrays have eliminated this illustration by a simple summation of two damage images.

Enhancing Radial Resolution Using Dual Arrays
It can be seen from the previous sections that the proposed TR-MUSIC algorithm has a good angle resolution for damage detection but a poor radial resolution. In order to eliminate this problem, dual arrays having an 'L' shape can be used for damage imaging. As shown in Figure 14, one horizontal and one vertical array containing 6 sensors are used to detect two damage in the aluminate panel. The same process is used for both horizontal and vertical arrays. Figure 15 gives comparisons of damage images between single array and dual arrays. It is obvious that although a single array has also a poor radial resolution, dual arrays have eliminated this illustration by a simple summation of two damage images.  It can be seen from the previous sections that the proposed TR-MUSIC algorithm has a good angle resolution for damage detection but a poor radial resolution. In order to eliminate this problem, dual arrays having an 'L' shape can be used for damage imaging. As shown in Figure 14, one horizontal and one vertical array containing 6 sensors are used to detect two damage in the aluminate panel. The same process is used for both horizontal and vertical arrays. Figure 15 gives comparisons of damage images between single array and dual arrays. It is obvious that although a single array has also a poor radial resolution, dual arrays have eliminated this illustration by a simple summation of two damage images.

Enhancing Radial Resolution Using Dual Arrays
(a) Image using the vertical arry (b) Image using the horizontal array (c) Image using the dual arrays

Experimental Section
An experiment was conducted to verify the feasibility and effectiveness of the proposed algorithm. As shown in Figure 16, the test sample was an aluminum panel with dimensions of 600 mm × 600 mm × 2 mm. Eleven PZT sensors that formed a linear piezoelectric sensor array were bonded to the surface of the aluminum panel. The diameter of the sensor is 8 mm, and the thickness is 0.48 mm. The distance of the adjacent sensors is 10 mm. Taking the panel center as the coordinate origin, the distance between the coordinate origin and the linear array is 180 mm. As shown in Figure 16c, a single-transmitting-and-multiplereceiving guided wave monitoring system based on the piezoelectric sensors developed by Xue et al. [33] was used to acquire signals. A five-peaked Hanning windowed sinusoidal tone burst with a center frequency of 70 kHz was used as the exciting signal. Holes were drilled in the aluminum panel to model the real damage. Two damage cases, including one damage and two damage, were considered in the experiment.

Experimental Section
An experiment was conducted to verify the feasibility and effectiveness of the proposed algorithm. As shown in Figure 16, the test sample was an aluminum panel with dimensions of 600 mm × 600 mm × 2 mm. Eleven PZT sensors that formed a linear piezoelectric sensor array were bonded to the surface of the aluminum panel. The diameter of the sensor is 8 mm, and the thickness is 0.48 mm. The distance of the adjacent sensors is 10 mm. Taking the panel center as the coordinate origin, the distance between the coordinate origin and the linear array is 180 mm. As shown in Figure 16c, a single-transmittingand-multiple-receiving guided wave monitoring system based on the piezoelectric sensors developed by Xue et al. [33] was used to acquire signals. A five-peaked Hanning windowed sinusoidal tone burst with a center frequency of 70 kHz was used as the exciting signal. Holes were drilled in the aluminum panel to model the real damage. Two damage cases, including one damage and two damage, were considered in the experiment.  To obtain accurate Green functions, a hybrid system combining a Nd:Yag pulse laser and a piezoelectric sensor was also developed, as shown in Figure 16c. The aluminum panel was fixed on an X-Y scanner to move to realize scanning of all the monitored area, while a Nd:Yag pulse laser was fixed on the optical platform to excite the Lamb wave To obtain accurate Green functions, a hybrid system combining a Nd:Yag pulse laser and a piezoelectric sensor was also developed, as shown in Figure 16c. The aluminum panel was fixed on an X-Y scanner to move to realize scanning of all the monitored area, while a Nd:Yag pulse laser was fixed on the optical platform to excite the Lamb wave propagating in the panel. Hence the plate can be monitored on a rectangular arrays for the positions of the laser source. The piezoelectric sensor array in Figure 16a was connected to an acquisition card to collect signals to calculate the green functions. Figure 17a,b shows a typical broadband signal acquired by the piezoelectric sensor when the pulse laser was used for excitation. To calculate the Green functions at a frequency of 70 kHz, the Gabor wavelet transform was used to extract the narrowband signal. As shown in Figure 18, the A0 wave was extracted to perform Fourier transform and calculate the Green function. For one exciting point of the pulse laser, 11 A0-wave signals acquired from the piezoelectric sensor array are labeled, x i (t), 1 ≤ i ≤ 11, and the Green function from this point to the sensor array can be written as where X i (ω) is the Fourier transform of x i (t). In the experiment, the rectangular area from (−100 mm, −100 mm) to (100 mm, 100 mm) was scanned by a pulse laser in steps of 1 mm, which generated 40,401 green function vectors. propagating in the panel. Hence the plate can be monitored on a rectangular arrays for the positions of the laser source. The piezoelectric sensor array in Figure 16a was connected to an acquisition card to collect signals to calculate the green functions. Figure 17a,b shows a typical broadband signal acquired by the piezoelectric sensor when the pulse laser was used for excitation. To calculate the Green functions at a frequency of 70 kHz, the Gabor wavelet transform was used to extract the narrowband signal. As shown in Figure 18, the A0 wave was extracted to perform Fourier transform and calculate the Green function. For one exciting point of the pulse laser, 11 A0-wave signals acquired from the piezoelectric sensor array are labeled, , and the Green function from this point to the sensor array can be written as where

Single Damage Case
One drilled hole with a diameter of 4 mm was set at the location (−57 mm, −37 mm) in the coordinates of Figure 16b. Figure 19 shows the scattering signal transmitted and received by both sensors. The first wave is a crosstalk. The A0 wave is obvious and used to establish the transfer function matrix and perform the time-reverse and eigenvalue decomposition processes. Figure 20 shows the histogram of the normalized eigenvalue. It propagating in the panel. Hence the plate can be monitored on a rectangular arrays for the positions of the laser source. The piezoelectric sensor array in Figure 16a was connected to an acquisition card to collect signals to calculate the green functions. Figure 17a,b shows a typical broadband signal acquired by the piezoelectric sensor when the pulse laser was used for excitation. To calculate the Green functions at a frequency of 70 kHz, the Gabor wavelet transform was used to extract the narrowband signal. As shown in Figure 18, the A0 wave was extracted to perform Fourier transform and calculate the Green function. For one exciting point of the pulse laser, 11 A0-wave signals acquired from the piezoelectric sensor array are labeled, , and the Green function from this point to the sensor array can be written as where

Single Damage Case
One drilled hole with a diameter of 4 mm was set at the location (−57 mm, −37 mm) in the coordinates of Figure 16b. Figure 19 shows the scattering signal transmitted and received by both sensors. The first wave is a crosstalk. The A0 wave is obvious and used to establish the transfer function matrix and perform the time-reverse and eigenvalue decomposition processes. Figure 20 shows the histogram of the normalized eigenvalue. It

Single Damage Case
One drilled hole with a diameter of 4 mm was set at the location (−57 mm, −37 mm) in the coordinates of Figure 16b. Figure 19 shows the scattering signal transmitted and received by both sensors. The first wave is a crosstalk. The A0 wave is obvious and used to establish the transfer function matrix and perform the time-reverse and eigenvalue decomposition processes. Figure 20 shows the histogram of the normalized eigenvalue. It can be clearly seen in Figure 20 that there is only one nonzero eigenvalue that is related to the damage. The noise subspace and Green functions are substituted into Equation (22), and the damage image can be obtained as shown in Figure 21. It can be seen that the damage image gives the correct location of the real damage, while it also shows the typical feature of MUSIC whereby the angle resolution is obviously larger than the radial resolution.
can be clearly seen in Figure 20 that there is only one nonzero eigenvalue th the damage. The noise subspace and Green functions are substituted into and the damage image can be obtained as shown in Figure 21. It can be damage image gives the correct location of the real damage, while it also sho feature of MUSIC whereby the angle resolution is obviously larger than the tion.   can be clearly seen in Figure 20 that there is only one nonzero eigenvalue th the damage. The noise subspace and Green functions are substituted into and the damage image can be obtained as shown in Figure 21. It can be damage image gives the correct location of the real damage, while it also sho feature of MUSIC whereby the angle resolution is obviously larger than the tion.   can be clearly seen in Figure 20 that there is only one nonzero eigenvalue that is related to the damage. The noise subspace and Green functions are substituted into Equation (22), and the damage image can be obtained as shown in Figure 21. It can be seen that the damage image gives the correct location of the real damage, while it also shows the typical feature of MUSIC whereby the angle resolution is obviously larger than the radial resolution.

Two-Damage Case
Two drilled holes with diameters of 4 mm were set at the locations (−57 mm, −37 mm) and (70 mm, 8 mm). Figure 22 shows the scattered signal transmitted and received by both sensors. The first wave is a crosstalk. The A0 wave is longer than Figure 19 for the twodamage scattering, and this is used to establish the transfer function matrix and perform

Two-Damage Case
Two drilled holes with diameters of 4 mm were set at the locations (−57 mm, −37 mm) and (70 mm, 8 mm). Figure 22 shows the scattered signal transmitted and received by both sensors. The first wave is a crosstalk. The A0 wave is longer than Figure 19 for the twodamage scattering, and this is used to establish the transfer function matrix and perform the time-reverse and eigenvalue decomposition processes. Local TR-MUSIC was used to analyze the A0 wave in Figure 22. Figure 23 gives the time-dependent eigenvalue function, in which there exist two peaks at the location of the scattering A0 wave. At 0.198 ms and 0.249 ms, damage images calculated by Equation (25) are given in Figure 24a,b, in which the locations of the pixel peaks show good consistency with the location of real damage. Combining Figure 24a-c gives the accurate locations of two real damages. However, due to the impact of noise, the TR-MUSIC damage image in the experiment has a worse radial accuracy than that in the simulation. sors 2021, 21, x FOR PROOF the time-reverse and eigenvalue decomposition processes. Local TR-M analyze the A0 wave in Figure 22. Figure 23 gives the time-depende tion, in which there exist two peaks at the location of the scattering A0 and 0.249 ms, damage images calculated by Equation (25) are given which the locations of the pixel peaks show good consistency with damage. Combining Figure 24a-c gives the accurate locations of two r ever, due to the impact of noise, the TR-MUSIC damage image in th worse radial accuracy than that in the simulation.   the time-reverse and eigenvalue decomposition processes. Local TR-MUSIC was used to analyze the A0 wave in Figure 22. Figure 23 gives the time-dependent eigenvalue function, in which there exist two peaks at the location of the scattering A0 wave. At 0.198 ms and 0.249 ms, damage images calculated by Equation (25) are given in Figures 24a,b, in which the locations of the pixel peaks show good consistency with the location of real damage. Combining Figure 24a-c gives the accurate locations of two real damages. However, due to the impact of noise, the TR-MUSIC damage image in the experiment has a worse radial accuracy than that in the simulation.

Conclusions
Damage imaging of plate-type structures based on Lamb waves can intuitively reveal the location and size of the damage, which has been a research topic of interest in the field of aircraft structural health monitoring. In this paper, an enhanced TR-MUSIC algorithm is proposed for multidamage detection in plate structures. The conclusions are as follows: (1) The TR-MUSIC algorithm is constructed by the orthogonality of signal and noise subspaces, which are divided by the eigenvalue decomposition of the TR operation of the transfer matrix. (2) An enhanced algorithm, the local TR-MUSIC algorithm, is proposed by using the moving time window to calculate the space spectrum of TR-MUSIC at different moments. (3) The TR-MUSIC algorithm can effectively detect damage with a higher angle location accuracy than the radial location accuracy. (4) The TR-MUSIC algorithm can break through the restriction of the Rayleigh criterion and realize the superresolution identification of multiple damage with distances smaller than a half-wavelength. (5) By using the local TR-MUSIC algorithm conducted by a moving time window, the proposed algorithm can detect multiple damages, even though the number of identified damage can break through the limitation of the number of sensor array elements.

Conclusions
Damage imaging of plate-type structures based on Lamb waves can intuitively reveal the location and size of the damage, which has been a research topic of interest in the field of aircraft structural health monitoring. In this paper, an enhanced TR-MUSIC algorithm is proposed for multidamage detection in plate structures. The conclusions are as follows: (1) The TR-MUSIC algorithm is constructed by the orthogonality of signal and noise subspaces, which are divided by the eigenvalue decomposition of the TR operation of the transfer matrix. (2) An enhanced algorithm, the local TR-MUSIC algorithm, is proposed by using the moving time window to calculate the space spectrum of TR-MUSIC at different moments. (3) The TR-MUSIC algorithm can effectively detect damage with a higher angle location accuracy than the radial location accuracy. (4) The TR-MUSIC algorithm can break through the restriction of the Rayleigh criterion and realize the superresolution identification of multiple damage with distances smaller than a half-wavelength. (5) By using the local TR-MUSIC algorithm conducted by a moving time window, the proposed algorithm can detect multiple damages, even though the number of identified damage can break through the limitation of the number of sensor array elements. (6) The damage image in the experiment has a worse radial accuracy than that in the simulation due to the effect of noise.
(7) Some time-varying parameters, like temperature, load and so on, have a significant influence on Lamb waves. Although there are some focusing on eliminating the effect of these parameters, especially the effect of temperature, it is necessary to develop a suitable compensation method for the proposed TR-MUSIC algorithm in the future.