Sage Journals: Discover world-class research

Abstract

In order to reconstruct the possible defects on the plate surface with arrays, a new wave tomography method based on the method of moments is established in this paper. According to the relationship between the probe number and grid amount, two algorithms, that is, the neural network and principal component analysis (PCA), are proposed and used to solve the ill-conditioned inversion equations. The neural network makes imaging feasible even if input data are not enough, and the PCA can greatly improve the computational efficiency via reducing the matrix dimension. Both numerical simulations and experimental measurements are conducted with the algorithm’s correctness and high precision validated. After investigating the influence of probe number on imaging quality, it is demonstrated that the algorithm can exactly predict the defect location when the input scattering data is not enough or fewer probes are arranged. More probes are needed for reconstructing the specific shape and thickness, especially when multiple defects are included. The qualitative results and quantitative data are conducive to providing some reference for engineering applications in nondestructive testing and structural health monitoring.

Keywords

Guided wave defect reconstruction method of moments neural network principal component analysis

Introduction

Lamb wave, a type of guided wave in plates, has become one of the hotspots in the field of nondestructive testing (NDT) and structural health monitoring (SHM) regions because of its long propagation distance and dispersive feature.^1,2 The wave velocity is related to the plate thickness, and then can be utilized for detecting corrosion defects.^3,4 Taking the A0 mode for example, it is very sensitive to the plate thickness when the value of fh is small. Here, f and h are the working frequency and plate thickness, respectively. Therefore, it is viewed as the priority when detecting surface defects. Additionally, it is more easily generated and received than bulk acoustic waves,^5–7 and then plenty of methods based on the A0 mode have been developed for NDT.

During the method development, the algorithm plays an important role in final defect visualization, which is directly related to the quality of final imaging. The defect profile is required to be reconstructed after utilizing a proper algorithm based on the wave amplitude,^8,9 phase,^10,11 wave number,^12,13 and other related characteristics of received signals from transducers. Therefore, how to develop a new algorithm and establish an accurate mathematical model to realize the qualitative and quantitative analysis according to different conditions, is of great significance for NDT or SHM.

Generally, typical algorithms for defect reconstruction include the elliptic imaging algorithm,^14,15 rapid imaging algorithm,^16,17 intelligent evolutionary algorithm,^18,19 etc. They are mainly interested in the macro-information carried by the signal from transducers, including time-of-flight and amplitude. Normally, these algorithms converge quickly and are suitable for judging the defect location. However, the imaging accuracy is not satisfactory, especially when identifying the defect profile using sparse arrays of transducers. In order to improve the resolution, the annular array with more probes is proposed and developed, as well as the corresponding algorithms, including diffraction tomography,²⁰ HARBUT algorithm,²¹ full-waveform inversion,²² and so forth. These algorithms are derived from the basic acoustic governing equations, in which the scattering mechanism of defects is considered. Mostly, the imaging resolution is satisfactory. However, enough data must be inputted when using these algorithms, and the transducer number N in an annular array must satisfy $N > \frac{4 π r}{λ}$ ,²³ where r is the circular radius and $λ$ is the wavelength of the illuminating wavefield. If $N \leq \frac{4 π r}{λ}$ , The resolution of these algorithms mentioned above will decrease evidently. If a good imaging is achieved via fewer transducers, it will be convenient for engineering applications. Therefore, it is urgent for us to develop a new algorithm that can reconstruct the defects, especially when fewer probes are arranged, which is the origin of this contribution.

In this paper, a reasonable inversion model is established by using the method of moment (MOM).^24–26 The highlight is that two algorithms are proposed according to the probe number arranged and the grid amount in the detection region. After that, the influence of different probe numbers on the imaging results is investigated in detail. When the equation is underdetermined, the corresponding neural network is designed, which can update the weight through the feedback of errors.²⁷ After inversion, the underdetermined equation is simplified into a single-layer neural network model and the network weight can achieve a good convergence. For overdetermined problems, the principal component analysis (PCA) is developed,²⁸ which can reduce the data dimension and efficiently improve computational efficiency. Additionally, PCA can also coordinate the main components, remove noise interference, and reduce the error to the greatest extent.

Firstly, the inversion equation is established based on MOM in “Design of tomography algorithm”, and the corresponding algorithms are introduced systematically, including the design scheme, convergence criterion and data processing. After that, numerical simulations are carried out with the aid of ABAQUS in Section “Numerical model and results.” The wave signals obtained are implanted into previous algorithms, and the defect reconstruction is shown. Utilizing the two algorithms, the influence of probe number on imaging quality is illustrated visually. In order to furthermore validate the algorithms, experimental measurements are conducted in Section “Experimental model and results,” which is based on a self-designed air-coupled system. Finally, some conclusions are summarized in Section “Conclusion.” It is expected that the algorithms developed in this paper can provide technical guidance for NDT or SHM technologies.

Design of tomography algorithm

Basic model

It is assumed that a circular area surrounded by dotted lines in the middle of a blue rectangular plate in Figure 1 is the detecting region, around which eight transducers are set at the same angle. If a two-dimensional Cartesian coordinate system is established at the circle center, the coordinate positions of all sensors can be obtained when the probe number and the radius are known. One transducer is used for wave generation and the others are receivers. Then, transform the emitting probe one by one and repeat the above process, the scattered field can be captured. This is a typical arrangement of sparse arrays when detecting defects inside the circle. Generally, more probes can produce better imaging results. Several array patterns, such as 16, 32, and 64 arrays with the same angle of configuration, can be designed to meet different requirements of detection. However, how to efficiently solve the governing equations and realize the defect reconstruction for any arrays, is challenging. It is anticipated that the present contribution can provide some guidance on this issue.

Figure 1.

Schematic diagram of a basic detection array.

In this paper, a single A0 mode with a frequency f is utilized for NDT. Because of dispersion, the phase velocity of A0 mode is dependent on the plate thickness if the working frequency is fixed. Then, the defects in the plate can be expressed indirectly with the aid of phase velocity. In the following sections, the guided wave tomography will be adopted here to reconstruct the phase velocity distribution in the unknown area and then convert it to thickness, so that defects can be detected visually. It should be stressed that the scattering waves with more modes can be generated if the defect boundary changes drastically.²⁹ However, there is no mode conversion as the thickness varies slowly.²¹ Therefore, it is assumed that corrosion-type defects do not lead to mode conversions. In order to illustrate the defect boundary as much as possible, an iterative solution is proposed in this paper with the specific algorithm demonstrated in the next section.

Establishment of the algorithm

The basic idea of this algorithm is to collect the scattering signal $U_{s} (r_{2})$ of each receiving point and establish the corresponding relationship $U_{s} (r_{2}) = f (c (r))$ between $U_{s} (r_{2})$ and phase velocity $c (r)$ , based on which $c (r)$ can be converted into the thickness distribution. Here, $r$ represents a position in the imaging area and $r_{2}$ is a receiving point. Taking eight probes equally distributed in Figure 1 as an example, the propagation path can be assumed as a direct ray mode, and then 56 sets of scattering fields are obtained, which are the initial input data. For analysis, the imaging area is usually divided into N grids, such as in Figure 2. Generally, $N > 56$ , that is, the input information is less than the output, which makes $U_{s} (r_{2}) = f (c (r))$ underdetermined. To improve the resolution, more probes need to be arranged, for example, 4032 inputs can be obtained if 64 arrays are distributed. On this condition, N is mostly larger than 4032, that is, $U_{s} (r_{2}) = f (c (r))$ is overdetermined. The corresponding algorithms for solving underdetermined and overdetermined equations will be proposed and demonstrated later.

Figure 2.

Grid diagram of the imaging area.

Originally, the Lippmann Schwinger equation requires:^30–32

U (r_{2}) = U_{u} (r_{2}) + U_{s} (r_{2}) = U_{u} (r_{2}) - \int G (r - r_{2}) O (r) U (r) d r

(1)

in which

U_{s} (r_{2}) = - \int G (r - r_{2}) O (r) U (r) d r

(2)

O (r) = {k_{0}}^{2} [{(\frac{c_{0}}{c (r)})}^{2} - 1]

(3)

Here, $O (r)$ is the objective function, $U (r_{2})$ represents the total field, and $U_{u} (r_{2})$ denotes the incident field, which can be calculated by Green’s function $G (r, r_{2})$ . $G (r, r_{2}) = - j / 4 \cdot H_{0}^{(2)} (k_{0} | r - r_{2} |)$ with $k_{0}$ and $c_{0}$ referring to the wavenumber and phase velocity of the A0 mode guided wave in the non-destructive region. j stands for the imaginary unit and ${H_{p}}^{(q)} (x)$ is the pth order Hankel function of class q.

Considering Equation (2) is an integral expression, it is difficult to directly use it to solve the objective function. Based on the grids described above, the corresponding transformation of Equation (1) can be obtained by utilizing MOM:^25,26

U (r_{i}) = U_{u} (r_{i}) + U_{s} (r_{i}) = U_{u} (r_{i}) + \sum_{j} O (r_{j}) U (r_{j}) c_{ij},

(4)

with

c_{ij} = {\begin{matrix} \frac{j}{2} \cdot [π k_{0} a {H_{1}}^{(2)} (k_{0} a) - 2 j], i = j; \\ \frac{j}{2} π k_{0} a \cdot J_{1} (k_{0} a) {H_{0}}^{(2)} (k_{0} R_{ij}), i \neq j . \end{matrix}

(5)

Here $i, j = 1, 2, . . ., N$ stands for the number of the grid point, $r_{i}$ is the location of the ith grid, and $c_{ij}$ is the integral form of Green’s function.³³ If the receiving transducer is used as the observation point, $r_{i}$ is replaced by $r_{m}$ with $m = 1, 2, . . ., M$ representing the identifier of one probe, and then

U_{s} (r_{m}) = \sum_{j} O (r_{j}) U (r_{j}) d_{mj} .

(6)

Here, $d_{mj}$ is a special case when the receiving transducer is viewed as the observation point. M is the number of transducers contained in the annular array. By writing Equations (4) and (6) into matrix form, we can get:

U^{(t)} = U^{(in)} + CO U^{(t)}

(7)

U^{(s)} = DO U^{(t)}

(8)

where $U^{(t)}$ and $U^{(in)}$ represent the total field and incident field at the grid point, respectively. They are both matrices of order $N \times 1$ . $U^{(s)}$ in Equation (8) is the scattering field measured by the receiving point, whose dimension is $[M \times (M - 1)] \times 1$ . $O = diag (O (r_{j}))$ is a diagonal matrix composed of the objective function with the dimension of $N \times N$ . $C$ and $D$ contain $c_{ij}$ and $d_{mj}$ with the dimension of $N \times N$ and $[M \times (M - 1)] \times 1$ , respectively. Equations (7) and (8) can be solved iteratively by utilizing the distorted Born iteration method (DBIM).^34,35 The method works on the basis of the Born approximation by treating the entire inspection area as a flat surface, and then, the total field is equivalent to the incident field. However, this assumption is not reasonable when the damage is large and the edge reflection is intense. It is necessary to continuously correct the objective function $O_{k}$ during the iterations in order to reduce this error. After that, the coefficient matrix $D_{k}$ is recalculated via $D_{k} = D (I - O_{k} C)^{- 1}$ . More complex defects require more iterations, which can be terminated by setting a convergence criterion. At the initial iteration, $U^{(t)} = U^{(in)}$ is substituted into Equation (8), and the initial objective function $O_{1}$ is obtained, based on which the coefficient matrix $D_{1}$ is modified. The results from the latest objective function will be substituted into Equation (7). The solution gradually converges to stability after several iterations between Equations (7) and (8). In brief, according to Born approximation, the initial total field can be replaced by the incident field, and the objective function will be updated by several repeated iterations.

Figure 3 is the basic flow of the imaging algorithm in this paper. It also returns to the original problem in this section, that is, how to solve the underdetermined and overdetermined equation when using different arrays, which will be highlighted in the next section.

Figure 3.

Flowchart showing the stages of iterative imaging algorithms.

Solving algorithm

Before iteration, replace $U^{(in)}$ with $U^{(t)}$ in Equation (8), and take the column elements of $U^{(in)}$ to form a new diagonal matrix ${U^{(in)}}^{*}$ , then Equation (8) can be transformed into ${U^{(s)}}_{[M (M - 1)] \times 1} = {D {U^{(in)}}^{*}}_{[M (M - 1)] \times N} \cdot$ ${O (r_{j})}_{N \times 1}$ . Let ${U^{(s)}} = b_{0}$ , ${D {U^{(in)}}^{*}} = A_{0}$ and ${O (r_{j})} = X_{0}$ , Equation (8) can be simplified into a general form $A_{0} X_{0} = b_{0}$ .

At the beginning of the iteration, the total field ${U_{k}}^{(t)}$ is calculated by using Equation (7) and the latest objective function $O_{k}$ . Green’s function $D_{k}$ is also modified in this step.^34,35 Here k = 0,1,2,3… represents the number of iterations. After that, calculate $Δ {U_{k}}^{(s)} = {U_{k}}^{(s)} - U^{(s)}$ and $Δ O_{k} = O_{k} - O_{k - 1}$ . Take the column elements of ${U_{k}}^{(t)}$ in Equation (8) to form a new diagonal matrix $U_{k} t^{*}$ , then the iteration can be changed into ${Δ {U_{k}}^{(s)}}_{[M (M - 1)] \times 1} = {D_{k} U_{k} t^{*}}_{[M (M - 1)] \times N} \cdot {Δ O_{k} (r_{j})}_{N \times 1}$ , which can be marked as $A_{k} Δ X_{k} = Δ b_{k}$ .

Therefore, the problem of solving an iterative process can be summarized in a general form:

AX = b

(9)

When Equation (9) is a system of underdetermined equations, the single-layer neural network is developed to iteratively find the stable root. On the contrary, if Equation (9) is overdetermined, PCA is utilized to solve it. The two algorithms will be introduced in detail below.

Neural network for solving underdetermined equations

Let $f (X) = AX - b$ , in which X is the training weight of the neural network, and then $f (X^{t}) = A X^{t} - b$ is the output. Here $t = 1, 2, . . .$ is the training times of the sample. If $δ^{t} = 0 - f (X^{t}) = - f (X^{t})$ is introduced as the training error, the performance index of this network can be expressed as $T = \frac{1}{2} {‖ δ^{t} ‖_{2}}^{2}$ with ${‖ ‖_{2}}^{2}$ representing the square of the 2-norm.

The weight can be adjusted by the following formula:

Δ X^{t} = - η \frac{d T}{d X^{t}} = - η \frac{d δ^{t}}{d f (X^{t})} \frac{d f (X^{t})}{d X^{t}} \frac{d T}{d δ^{t}} = η A^{T} δ^{t}

(10)

Here $0 < η < 1$ is the learning rate, and in general can be obtained after applying Equation (10):

X^{t + 1} = X^{t} + Δ X^{t} = X^{t} + η A^{T} δ^{t}

(11)

with $A^{T}$ representing the transposition of $A$ . Equation (11) is a formula for finding roots using a single-layer neural network, which is suitable for solving the underdetermined problem like Equation (9).

PCA for solving overdetermined equations

To avoid complex matrix operations, this method does not follow the standard process of PCA. At first, singular value decomposition is performed on $A$ , that is,

A = U Λ V^{T} = C V^{T}

(12)

in which $V$ is a matrix composed of eigenvectors from $A^{T}$ and $C$ is the principal component matrix. $Λ$ and $A$ have the same dimension with its principal diagonal containing the singular value of $A$ . The column vector of $U$ is the standard orthogonal eigenvector of $A^{H} A$ , in which $A^{H}$ is the transpose conjugate.

Let $Y = V^{T} X$ , then Equation (9) can be attributed to solve $CY = b$ . Select the first p columns of $C$ , $V$ and $Y$ , that is, $C = (C_{1}, C_{2})$ , $V = (V_{1}, V_{2})$ and $Y = (Y_{1}, Y_{2})$ . After solving $C_{1} Y_{1} = b$ , $Y_{1} ~ {C_{1}}^{+} b$ can be obtained, and then $X ~ ({V_{1}}^{T})^{+} Y_{1}$ . Here + denotes generalized inverse, and ~ means taking the least square solution. When the matrix has full rank, it is equivalent to the inverse matrix. Since $V_{1}$ is a unitary orthogonal matrix, we may get $X ~ V_{1} Y_{1}$ . Utilizing this, an estimate of $X$ can be denoted as:

\tilde{X} = V_{1} Y_{1}

(13)

Equation (13) is a formula for finding roots using PCA. Correspondingly, $‖ A \tilde{X} - b ‖_{2}$ can be used to calculate the error. The principal component p can be selected by observing the cumulative percentage $s (p)$ of eigenvalue,²⁸ which also can be calculated via

s (p) = \frac{\sum_{i = 1}^{p} λ_{i}}{\sum_{i = 1}^{q} λ_{i}}

(14)

when $p \leq q$ , in which q represents the number of eigenvalues of $A^{T} A$ .

Convergence criterion

Section “Solving algorithm” only describes two methods for solving ill-posed equations but does not provide a convergence criterion for iteration. As defined by relative errors (RE), the following equation can be obtained:

RE = {‖ X_{k} - X_{exact} ‖}_{2} / {‖ X_{exact} ‖}_{2}

(15)

where $X_{k}$ represents the solution of the k^th iteration, $X_{exact}$ is the exact solution. In practice, $X_{exact}$ is usually unknown. In order to verify the approximation between the solution of the algorithm and the real solution, a relative residual error is redefined in this paper:

RRE = ‖ A_{k} X_{k} - b ‖_{2} / ‖ b ‖_{2} = ‖ D_{k} {U_{k}}^{(t) *} O_{k} - U^{(s)} ‖_{2} / ‖ U^{(s)} ‖_{2}

(16)

Usually, the iteration is stopped when $RRE \leq 0.01$ , which can efficiently ensure the calculation accuracy in engineering applications. However, this paper chooses a stricter limit. Considering that the signal is more likely to be affected by noise with the massive input information when the acquisition array possesses enough scattering data, the criterion is adjusted to $RRE \leq 0.001$ . It should be stressed that the computation of the neural network in the sparse array is small, and more iterations will not burden the computer. In addition, the PCA algorithm itself is efficient and can support more iterations to meet higher precision.

Data processing

Equations (7) and (8) of MOM are derived from the primal Lippmann-Schwinger equations, which are calculated in the frequency domain. Therefore, extracting the scattering field in the frequency domain during simulations or experiments plays a key role. The most basic scattering field can be achieved by the following expression:

U^{(s)} = U^{(t)} - U^{(in)}

(17)

Here $U^{(s)}$ , $U^{(t)}$ and $U^{(in)}$ respectively represent the scattering field, total field, and incident field in the frequency domain. $U^{(t)}$ can be obtained by extracting the displacement of z-direction in the time domain and making a fast Fourier transform (FFT) on it when simulating or measuring. Figure 4 is an example from a set of experimental data with a narrowband excitation. Measure the defective and flawless signals along the same propagation path, intercept the direct wave component using a reasonable window function and then conduct a FFT. The scattered field can be achieved by performing a subtraction between the defective field and that of a flat plate without defects. Repeat the above operation, numerous scattered fields can be obtained as the input components.

Figure 4.

A set of experimental received signals.

In general, the incident field is equivalent to the total field of the non-destructive region. Thus, one method to obtain $U^{(in)}$ is to prepare a signal in a flat plate without defects beforehand. The other method is directly calculating the Green’s function introduced above. Strictly speaking, the first method is more convincing. However, $U^{(in)}$ of each propagation is different, requiring repeated measurements of all flawless signals, which is tedious and difficult to realize in practical conditions. By contrast, the second method using Green’s function theory is more practical. Since the approach is theoretically calculated, the real incident field would have some phase deviation due to the detection environment and noise. Additionally, the simulation can be also subject to errors caused by the grid. Therefore, some processing of such theoretical calculations is also necessary. To reduce the error between simulation (or experiment) and the governing equation of MOM, a calibration factor can be proposed as:^20–22,36

Q = \frac{U^{(in, FE or EXP)}}{U^{(in)}}

(18)

in which $U^{(in)}$ represents the theoretical value of the incident field calculated by Green’s function, and $U^{(in, FE or EXP)}$ is the result of FFT of a non-destructive signal in the time domain selected by simulation (or experiment). $U^{(in, FE or EXP)}$ can be measured by prior visual observation of a set of undamaged inspection locations. Q would be calculated for a group of measurements. This prevents the dilemma of measuring flawless signals in all paths. The incident field is used to replace the total field in the first iteration when using Born approximation. Equation (8) can be rewritten as:

U^{(s, FE or EXP)} = DO U^{(in, FE or EXP)}

(19)

Here, $U^{(s, FE or EXP)}$ stands for scattering field in the frequency domain gained from simulations (or experiments). After introducing the calibration factor, Equation (19) can be rewritten as

\frac{U^{(s, FE or EXP)}}{Q} = DO U^{(in)}

(20)

The utilization of Q to calibrate the scattering field is a very important step in the iterative process of inversion, which can reduce the error of simulations (or experiments). Evidently, the above methods of data processing mainly introduce the extraction of the scattering field and the proposal of the calibration factor.

After illustrating details associated with signal processing, the imaging method will be demonstrated. As mentioned above, a combination of experimental and theoretical calculations is used to obtain $U^{(s)}$ and the algorithm DBIM using iterative iterations is eventually established to yield the objective function $O (r) = k_{0}^{2} [{(\frac{c_{0}}{c (r)})}^{2} - 1]$ , from which $c (r) = \frac{k_{0} c_{0}}{\sqrt{O (r) + k_{0}^{2}}}$ is deduced. In other words, the output from the algorithm has been transformed into the phase velocity $c (r)$ . Utilizing the dispersion curve of A0 mode Lamb wave, that is, the phase velocity versus thickness multiplied by frequency, the plate thickness value at each point in the detection area is formalized, and corrosion-type defects can be obtained. In the next section, both Finite Element Methods and experimental measurements will be carried out to verify the tomography algorithm.

Numerical model and results

Numerical model

In order to verify the algorithm, ABAQUS software is utilized to simulate several aluminum plates with different kinds of defects. The plate size is chosen as 1200 mm × 1200 mm × 3 mm, and the diameter for detection is 400 mm. For material parameters, Young’s modulus, Poisson ratio, and density are 70.8 GPa, 0.34, and 2700 kg/m³, respectively. Additionally, several layers of materials with increasing Rayleigh damping coefficients are added around the aluminum plate edges to form an absorption layer,^37–39 which are used to fully absorb the reflected wave from plate boundaries. The source excitation is a 5-cycle Hanning windowed toneburst signal at a central frequency of 200 kHz. C3D8R elements are distributed with at least three meshes adopted along the plate thickness direction. Recomputation has been carried out with smaller elements before final simulations, which can make the numerical results convergent and accurate.

Figure 5 shows the thickness distribution of four theoretical models. Both Figure 5(a) and (b) belong to a single defect model with the defect center at (0, 0). The thickness of the triangular defect in Figure 5(a) is uniform, while the elliptical defect in Figure 5(b) has a gradient change in thickness. Additionally, two defect models are also considered in Figure 5(c) and (d), in which Figure 5(c) has double circle defects with different sizes and Figure 5(d) consists of a triangle defect and an elliptical one.

Figure 5.

Theoretical models. (a) An equilateral triangle defect locates at (0, 0), and the side length is 80 mm with a remaining thickness of 1 mm. (b) The elliptical defect locates at (0, 0), and the long and short axes are respectively 80 and 40 mm with a remaining thickness of 1 mm. The defect thickness varies smoothly from its edge to the plate. (c) The large circle defect centered at (−100, 0 mm) has a radius of 40 mm and a remaining thickness of 2 mm, and the small circle defect located at (100, 0 mm) has a radius of 20 mm and a remaining thickness of 1 mm. (d) The equilateral triangle defect centered at (0, 60 mm) has a side length of 48 mm and a remaining thickness of 1 mm. The elliptical defect centered at (0, −60 mm), whose long axis and short axis are respectively 70 and 40 mm, the remaining thickness is 2 mm.

Reconstruction results

For detecting the defects in Figure 5, the detection region is divided into 1793 points. The array with 8, 16, 32, and 64 probes are used, respectively. The first three models obtain 56, 240, and 992 groups of input scattering data, respectively, and all of them select the neural network method for iterative solutions. The 64-array case can receive 4032 groups of scattering fields, and the PCA algorithm is selected to avoid noise pollution caused by redundant data. Correspondingly, the reconstruction results are shown in Figures 6 to 9. Meanwhile, to visually evaluate image quality, the thickness distributions across the defect center are extracted and illustrated in Figure 10. By contrast, the theoretical value of the actual defect is also included.

Figure 6.

Reconstruction results of the triangular defect by (a) 8-array model, (b) 16-array model, (c) 32-array model, and (d) 64-array model.

Figure 7.

Reconstruction results of the elliptical defect by (a) 8-array model, (b) 16-array model, (c) 32-array model, and (d) 64-array model.

Figure 8.

Reconstruction results of the double circular defects by (a) 8-array model, (b) 16-array model, (c) 32-array model, and (d) 64-array model.

Figure 9.

Reconstruction results of triangular and elliptical defects by (a) 8-array model, (b) 16-array model, (c) 32-array model, and (d) 64-array model.

Figure 10.

Cross-sections of the thickness reconstructions across the defect center from (a) Figure 5(a), (b)Figure 5(b), (c)Figure 5(c), and (d)Figure 5(d). The yellow, green, blue, and black lines represent the results calculated from the 8, 16, 32, and 64 probes, respectively, and the red solid line corresponds to the actual defect.

From a whole perspective, the reconstruction results of MOM can effectively distinguish the defect location. The shape and thickness information of defects are closer to reality as the probe number increases. The imaging of these four types calculated by eight probes is relatively poor. For single defect such as Figure 5(a) and (b), the 8-array model can only roughly make a distinction between noise and defect. For the two defects, it is even difficult to make a qualitative evaluation with a lot of overlap between the defect and surrounding noise. By contrast, the 16-array model can reconstruct a relatively good imaging, because it can effectively seek the damaged area. From the quantitative reconstruction results in Figure 10(a) and (b), it can be seen that the error is controlled within 0.2 mm for the single defect cases. Unfortunately, it is still not ideal for identifying the defect shape. The neural network method combined with 32 probes is able to conduct a quantitative evaluation in the single defect, with the contour and depth very close to the reality. However, when facing multiple defects with more scattering forms, there are still considerable errors in the reconstruction of thickness. As the probe number further increased to 64, the shape and depth of defects can be reconstructed accurately with the noise pollution at a very low level, and the reconstruction results are best. Totally speaking, Figures 6 to 10 efficiently validate the correctness of the wave tomography method proposed in this paper. It should be stressed that the calculation speed of PCA is fast with only partial principal components being selected to participate. For example, when inputting data from the 64-array is adopted, only less than 3 minutes are needed at a time. Therefore, this method is convenient to be applied in NDT.

Reviewing Figure 10, it also can be concluded that a defect with a small gradient variation along thickness is easier to be reconstructed. For illustrating this, the thickness value on the cross-sections in Figure 10 is extracted and a parameter R_NL is defined as

R_{NL} = 1 - \sqrt{[\sum_{i} {(T_{exa}^{i} - T_{cal}^{i})}^{2}] / \sum_{i} (T_{exa}^{i})^{2}}

(21)

in which $T_{cal}^{i}$ is the thickness value of a reconstructed point on the cross-section, and $T_{exa}^{i}$ is the real thickness. When the R_NL is closer to 1, the reconstruction error is smaller, which indicates that the defect is easier to be reconstructed. The R_NL variation versus the probe number is shown in Figure 11, from which it can be seen that the elliptical defect with a gradient variation is easiest to be depicted. By contrast, the case of double circular defects is the most difficult to be reconstructed, because the defect area of this case is the largest and the Born approximation brings greater initial error.

Figure 11.

The R_NL variation for different defects.

Experimental model and results

Experimental model

For furthermore validating the correctness of algorithms, as well as the influence of probe number on reconstructing different defects, experiments are carried out. An equilateral triangle defect and an elliptical defect that are the same as in Figure 5(d) are fabricated in an aluminum plate, such as in Figure 12. A pair of non-contact air-coupled probes (Japan Probe 0.2 K 14 × 20 N) are designed as transducers, and probes are located at the same position as those in numerical simulations. A 5-cycle Hanning windowed toneburst signal with the central frequency of 200 kHz is chosen as the wave generation, and the excitation and reception are integrated into a set of NDT system (JPR-600C). For eliminating the unnecessary reflected waves, damping glues are attached to the aluminum at plate edges. The air-coupled probes are suspended above the reference points, and the distance between them will be continuously changed to meet different paths of detection. For exciting a relatively pure A0 mode Lamb wave,^7,40 the deflection angle of the probe should be arranged according to Snell law, that is,

θ_{A_{0}} = \sin^{- 1} (c_{air} / c_{Lamb})

(22)

Here $c_{air} = 340 m / s$ is the propagation velocity of ultrasonic waves in air, and $c_{Lamb} = 2009 m / s$ is the phase velocity of A0 mode guided wave in the 3 mm aluminum plate when f = 200 kHz. Therefore, $θ_{A 0} = 9.74 ° \approx 10 °$ is chosen as the deflection angle. It should be stressed that the received signal has about 10 cycles, not 5 cycles, which can be attributed to two reasons. Firstly, the dispersion is inevitable, which makes the wave signal distorted. On the other hand, the air-coupled transducer has a built-in filter, and the narrowband signal will increase the wavepacket number. In spite of these, it has a tiny impact on practical signal processing in the frequency domain.

Figure 12.

Experimental setup of an air-coupled nondestructive testing system

Reconstruction results

8, 16, 32, and 64 arrays are respectively arranged, and the signals received are treated by the neural network or PCA after FFT, with the imaging results shown in Figure 13. Meanwhile, the thickness distribution across the defect center is extracted and illustrated in Figure 14.

Figure 13.

Reconstruction results of triangular and elliptical defects using experimental signals from (a) 8 probes, (b) 16 probes, (c) 32 probes, and (d) 64 probes.

Figure 14.

Cross-sections of the thickness reconstructions by experiments. The yellow, green, blue, and black lines represent the results calculated from the 8, 16, 32, and 64 probes, respectively, and the red solid line corresponds to the actual defect.

In summary, sparse arrays, such as 8 or 16 probes, can only identify the defect location. The thickness information is inaccurate. As the probe number increases, the defect profile can be distinguished gradually. For the four cases, the noise pollution of the imaging results is more serious than that in numerical simulations, which may be due to the deviation of the measurement system. In addition, the air-coupled detection system is not limited by the coupling agent. It is difficult to accurately control the exciting point at oblique incidence. In addition, the flat probe has a finite surface area for wave excitation and reception, which also can bring in discrepancies with simulation conditions. Despite those, the results based on experimental signals match well with the numerical simulations. Similarly, R_NL can be obtained and shown in Figure 15, which also demonstrates that the thickness distribution in experiments is slightly worse than that of simulations. The depth error from the 64-array imaging is less than 0.25 mm, which efficiently validates the algorithm’s effectiveness once again.

Figure 15.

The variation trend of R_NL in experiments

Final remarks

The experimental results demonstrate that the algorithm is effective to obtain the thickness imaging for plate-like structures in different arrays. Although the diffraction tomography based on the Born approximation can also reconstruct defects,²⁰ it still has limitations. The defect must be relative to the wavelength λ and of low contrast relative to the background. Therefore, a lower frequency is usually selected when detecting large flaws.^20,21 Accompanying that, the received signal of guided waves at low frequencies is commonly noisier and the imaging resolution will be reduced. Unlike diffraction tomography, the DBIM method adopted in this algorithm overcomes this drawback as the objective function is modified at each iteration step for better convergence.

The iterative approach, widely utilized in full waveform inversion for seismic waves,²² can also be adapted to reconstruct corrosion-type defects in plates. It usually requires 20–40 iterations, which takes 4.5 h.²² With the same 64-array conditions set in that case, only 8 iterations are needed in this paper, and each iteration takes less than 3 min when using a common computer with 11th Gen Intel Core i7-11700 @ 2.50 GHz. The Green’s function is recalculated at each step when the objective function is obtained, so that the new scattering features are properly expressed, which induces faster convergence. Additionally, the PCA algorithm decreases the matrix dimensionality, which further improves the computing efficiency.

The rule $N > \frac{4 π r}{λ}$ is usually kept when quantifying imaging using arrays,^20–23 in which: r, λ and N respectively stand for the circular radius, the wavelength, and the probe number. However, it is usually hard to ensure sufficient probes in engineering applications. Therefore, this paper attempts to obtain a good image using sparse arrays for a qualitative description of defects. Frankly speaking, although the DBIM method can compensate the shortcomings of Born approximation to some extent, the iteration still needs much computing costs if the defects are too complex, which is expected to be conquered in the future.

Conclusion

In summary, two algorithms that can be suitable for any number of probes are developed and utilized for the reconstruction of various defects. The influence of probe number is systematically investigated from perspectives of simulation and experiments, with the main conclusions drawn as follows:

(1) The MOM has the ability to transform the nonlinear Lippmann-Schwinger equation into a linear form, and the corresponding relationship between the objective function and the thickness is directly established.

(2) Two algorithms of machine learning, that is, neural network and PCA, are developed, both of which have high calculation accuracy when solving the underdetermined and overdetermined equations, respectively.

(3) The defect location can be roughly distinguished if the array is sparse. More probes are needed for reconstructing the specific shape and thickness, especially when multiple defects are included.

The research conducted in this paper is based on a single-mode guided wave. For improving the imaging resolution, a relatively short wave with the central frequency of 200 kHz is utilized. Meanwhile, in order to eliminate the interference of S0 mode, an air-coupled system is designed, which is used for wave generation and reception. A laser vibrometer may be better for improving imaging quality if experimental conditions permit, which is worthy of exploration in the future.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors gratefully acknowledge the supports of the Natural Science Foundation of China (12061131013, 12172171, and 11972276), the Opening Project of Applied Mechanics and Structure Safety Key Laboratory of Sichuan Province (SZDKF-202002), the State Key Laboratory of Mechanics and Control of Mechanical Structures at NUAA (No. MCMS-I-0522G01), National Natural Science Foundation of Jiangsu Province (BK20211176), Jiangsu High-Level Innovative and Entrepreneurial Talents Introduction Plan (Shuangchuang Doctor Program, JSSCBS20210166), the Fundamental Research Funds for the Central Universities (NS2022011 and NE2020002), Local Science and Technology Development Fund Projects Guided by the Central Government (2021Szvup061), and a project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD). Prof. Iren Kuznetsova thanks the Russian Ministry of Science and Higher Education (government task #FFWZ-2022-0002) for partial financial supports.

ORCID iDs

Zhi Qian

Yinghong Zhang

References

Luo

Zeng

Lin

. A novel time-frequency transform for broadband Lamb waves dispersion characteristics analysis. Struct Heal Monit Int J 2021; 20(6): 3056–3074.

Livadiotis

Ebrahimkhanlou

Salamone

. Monitoring internal corrosion in steel pipelines: a two-step helical guided wave approach for localization and quantification. Struct Heal Monit 2021; 20(5): 2694–2707.

Mitra

Gopalakrishnan

. Guided wave based structural health monitoring: a review. Smart Mater Struct 2016; 25(5): 53001–53027.

Grahn

. Lamb wave scattering from a circular partly through-thickness hole in a plate. Wave Motion 2003; 37(1): 63–80.

Nishino

Asano

Taniguchi

, et al. Precise measurement of pipe wall thickness in noncontact manner using a circumferential lamb wave generated and detected by a pair of air-coupled transducers. Jpn JAppl Phy 2011; 50(7): 07HC10.

Tat

Georgeson

Bossi

. Evaluation of air coupled ultrasound for composite aerospace structure. In: Thompson

Chimenti

(eds.) Review of progress in quantitative nondestructive evaluation. 2009, vol. 1096(35), pp. 912–919.

Rose

. Ultrasonic guided waves in solid media. West Nyack, NY: Cambridge University Press, 2014.

Lin

Wilkins

Rao

, et al. Corrosion detection with ray-based and full-waveform guided wave tomography. In: Yu

Shull

(eds.) Nondestructive characterization and monitoring of advanced materials, aerospace, civil infrastructure, and transportation VIX. Bellingham, WA: SPIE, 2020, p. 11380.

Koduru

Rose

. Mode controlled guided wave tomography using annular array transducers for SHM of water loaded plate like structures. Smart Mater Struct 2013; 22(12): 125021.

10.

Huthwaite

. Guided wave tomography with an improved scattering model. Proc R Soc A: Math Phys Eng Sci 2016; 472(2195): 20160643.

11.

Huthwaite

. Evaluation of inversion approaches for guided wave thickness mapping. Proc R Soc A: Math Phys Eng Sci 2014; 470(2166): 20140063.

12.

Liu

, et al. Quantitative detection of lamination defect in thin-walled metallic pipe by using circumferential Lamb waves based on wavenumber analysis method. NDT E Int 2019; 102: 56–67.

13.

Liu

Feng

, et al. Quantitative rectangular notch detection of laser-induced lamb waves in aluminium plates with wavenumber analysis. Trans Nanjing Univ Aeronaut Astronaut 2018; 35(2): 244–255.

14.

Chen

, et al. An improved discrete elliptic imaging algorithm based on guided waves for defect localisation in curved plates. Insight 2019; 61(11): 656–662.

15.

Luo

, et al. A new probabilistic ellipse imaging method based on adaptive signal truncation for ultrasonic guided wave defect localization on pressure vessels. Sensors 2022; 22(4): 1540.

16.

Bian

Rose

. Sparse array ultrasonic guided wave tomography. Mater Eval 2005; 63(10): 1035–1038.

17.

Hay

Royer

Gao

, et al. A comparison of embedded sensor Lamb wave ultrasonic tomography approaches for material loss detection. Smart Mater Struct 2006; 15(4): 946–951.

18.

Chen

Liu

, et al. An intelligent algorithm based on evolutionary strategy and clustering algorithm for Lamb wave defect location. Struct Heal Monit Int J 2021; 20(4): 2088–2109.

19.

Chen

Liu

Gong

, et al. Evolutionary strategy-based location algorithm for high-resolution lamb wave defect detection with sparse array. IEEE Trans Ultrason Ferroelectr Freq Control 2021; 68(6): 2277–2293.

20.

Belanger

Cawley

Simonetti

. Guided wave diffraction tomography within the born approximation. IEEE Trans Ultrason Ferroelectr Freq Control 2010; 57(6): 1405–1418.

21.

Huthwaite

Simonetti

. High-resolution imaging without iteration: a fast and robust method for breast ultrasound tomography. J Acoust Soc Am 2011; 130(3): 1721–1734.

22.

Rao

Ratassepp

Fan

. Guided wave tomography based on full waveform inversion. IEEE Trans Ultrason Ferroelectr Freq Control 2016; 63(5): 737–745.

23.

Simonetti

Huang

Duric

. On the spatial sampling of wave fields with circular ring apertures. J Appl Phys 2007; 101(8): 083103.

24.

Zhang

Qin

. A regularized method of moments for three-dimensional time-harmonic electromagnetic scattering. Appl Maths Lett 2021; 112: 106746.

25.

Johnson

Tracy

. Inverse scattering solutions by a sinc basis, multiple source, moment method—Part I: theory. Ultrason Imaging 1983; 5(4): 361–375.

26.

Tracy

Johnson

. Inverse scattering solutions by a sinc basis, multiple source, moment method—Part II: numerical evaluations. Ultrason Imaging 1983; 5(4): 376–392.

27.

Wang

Zeng

Chen

. Back propagation neural network with adaptive differential evolution algorithm for time series forecasting. Expert Syst Appl 2015; 42(2): 855–863.

28.

Abdi

Williams

. Principal component analysis. WIREs. Comput Stat 2010; 2(4): 433–459.

29.

Khalili

Cawley

. The choice of ultrasonic inspection method for the detection of corrosion at inaccessible locations. NDT E Int 2018; 99: 80–92.

30.

Pereira

Schmidt

AGM

. Exact solutions for Lippmann–Schwinger equation for the scattering by hyper-spherical potentials. Few-Body Syst 2022; 63(1): 25.

31.

Eikrem

Nævdal

Jakobsen

. Iterative solution of the Lippmann-Schwinger equation in strongly scattering acoustic media by randomized construction of preconditioners. Geophys J Int 2021; 224(3): 2121–2130.

32.

Kak

Slaney

Wang

. Principles of computerized tomographic imaging. Med Phys 2002; 29(1): 107.

33.

Richmond

. Scattering by a dielectric cylinder of arbitrary cross section shape. IEEE Trans Antennas Propag 1965; 13(3): 334–341.

34.

Chew

Wang

. Reconstruction of 2-dimensional permittivity distribution using the distorted born iterative method. IEEE Trans Med Imaging 1990; 9(2): 218–225.

35.

Haddadin

Ebbini

. Self-focusing arrays for imaging and therapy through inhomogeneous media. In: Levy

Schneider

McAvoy

(ed.) 1996 IEEE ultrasonics symposium, proceedings, San Antonio, TX, vol 1 and 2, pp. 1563–1566.

36.

Belanger

Cawley

. Feasibility of low frequency straight-ray guided wave tomography. NDT E Int 2009; 42(2): 113–119.

37.

Liu

Quek Jerry

. A non-reflecting boundary for analyzing wave propagation using the finite element method. Finite Elem Anal Des 2003; 39(5): 403–417.

38.

Moreau

Castaings

. The use of an orthogonality relation for reducing the size of finite element models for 3D guided waves scattering problems. Ultrasonics 2008; 48(5): 357–366.

39.

Castaings

Bacon

. Finite element modeling of torsional wave modes along pipes with absorbing materials. J Acoust Soc Am 2006; 119(6): 3741–3751.

40.

Auld

Green

. Acoustic fields and waves in solids: two volumes. Phys Today 1974; 27(10): 63–64.

A novel wave tomography method for defect reconstruction with various arrays

Abstract

Keywords

Introduction

Design of tomography algorithm

Basic model

Establishment of the algorithm

Solving algorithm

Neural network for solving underdetermined equations

PCA for solving overdetermined equations

Convergence criterion

Data processing

Numerical model and results

Numerical model

Reconstruction results

Experimental model and results

Experimental model

Reconstruction results

Final remarks

Conclusion

Footnotes

Declaration of conflicting interests

Funding

ORCID iDs

References