A Modified Wavenumber Algorithm of Multi-Layered Structures with Oblique Incidence Based on Full-Matrix Capture

Full-matrix capture (FMC)-based ultrasonic imaging provides good sensitivity to small defects in non-destructive testing and has gradually become a mainstream research topic. Many corresponding algorithms have been developed, e.g., the total focusing method (TFM). However, the efficiency of the TFM is limited, especially in multi-layered structures. Although the appearance of wavenumber algorithms, such as extended phase-shift migration (EPSM) methods, has improved imaging efficiency, these methods cannot be applied to cases with oblique incidence. Therefore, a modified wavenumber method for full-matrix imaging of multi-layered structures with oblique array incidence is proposed. This method performs a coordinate rotation in the frequency domain to adapt it to the oblique incidence. It then utilizes wave-field extrapolation to migrate the transmitting and receiving wave field to each imaging line, and a correlation imaging condition is used to reconstruct a total focused image. The proposed method can deal with any incident angle without precision loss. Moreover, it inherits the computational efficiency advantages of the wavenumber algorithms. The simulation and experimental results show that the proposed method performs better in terms of accuracy and efficiency than the TFM. Specifically, it is nearly 60 times faster than the TFM when processing an FMC dataset with a size of 4096 × 64 × 64.


Introduction
Complex structures such as welds, joints, and pipe bends are critical components in massive buildings, bridges, and oil pipelines, where internal defects pose considerable threats [1][2][3][4]. Therefore, an effective and accurate nondestructive testing (NDT) technique to detect and characterize defects is essential to ensure the safety and reliability of these components. The ultrasonic phased array, as a well-established technique, is widely used in industrial NDT applications [5][6][7][8] due to its portability, reliability, and relatively high accuracy. Conventional ultrasonic phased array techniques are based on delay-and-sum (DAS) beamforming, which generates coherent ultrasonic beams by controlling the time delay of transmitters [9,10]. However, the limited information in the data in these solutions leads to low image contrast and distortion of the target object and severely restricts the applications for complex structures. The appearance of full-matrix capture (FMC) has solved this problem [11]. FMC performs data acquisition by having all elements act as receivers while every successive element acts as a transmitter. Hence, FMC datasets contain more information about the objects, and the focusing process can be performed offline without applying delay laws during acquisition [12].
On the basis of FMC, many algorithms have been developed independently. The total focusing method (TFM) [13,14] is a typical post-processing algorithm that uses all the FMC data to focus at every pixel in the region of interest (ROI). It can significantly improve detection sensitivity and spatial resolution and is considered to be the gold 2 of 17 standard in ultrasonic imaging. However, the computational cost of TFM is substantial, as it is an extension of the DAS approach in the time domain [15]. In order to resolve the computational issue, some researchers have sought for solutions in the frequency domain. Hunter et al. [16] proposed a wavenumber algorithm for full-matrix imaging, which has superior computational performance. Nevertheless, this algorithm is only applicable for single-layered structures. To date, frequency-domain algorithms have been proposed to accommodate more complicated structures, such as the extended phase-shift migration (EPSM) method [17] and multi-layered wavenumber algorithms [18,19]. These methods are adapted to multi-layered structures and have improved the computational efficiency.
The algorithms mentioned above are only valid for a multi-layered medium with a surface parallel to the phased array. In practice, however, complex structures such as welds and joints usually require oblique incidence detection to obtain excellent inspection coverage [20,21]. In this situation, the TFM needs to use the Fermat principle, calculating the refraction points at each interface to accurately estimate the wave propagation time in the multi-layered medium [18]. However, this complicated step increases the computational cost dramatically. Several other approaches to full-matrix imaging have also been proposed for oblique incidence detection. Ray tracing methods [22,23] estimate the propagation time through iterative calculations, no matter how complicated the interface is, but the complexity of iterative calculations limits their efficiency. The root-mean-square (RMS) velocity [24,25] has been applied to reduce the calculation time. However, the approximation of RMS velocity introduces a huge error when the incidence angle of the array is overlarge. The virtual source TFM [26,27] can achieve total focusing through arbitrarily shaped interfaces. In addition, it reduces the processing complexity of time-of-flight calculations, leading to more efficient implementation of the TFM. Unfortunately, although many attempts have been made to improve the efficiency, these time-domain algorithms still cannot meet the high-resolution real-time imaging demand. To solve this, Lukomski proposed an FMC method based on phase-shift migration (PSM) [28], implemented in a wavenumber domain. Unlike EPSM, this method allows for imaging of a multi-layered medium with lateral velocity variations, such as a medium that is not perpendicular to the surface of the arrays. The main advantage of the algorithm is a much higher efficiency than time-domain algorithms. Nevertheless, this algorithm cannot achieve a spatial resolution as high as the standard TFM, and the signal-to-noise ratio (SNR) is lower than for the TFM. Thus, a high-efficiency and high-accuracy algorithm for ultrasonic full-matrix imaging in oblique incidence detection on multi-layered structures is intensely desired.
In this paper, a modified wavenumber method for full-matrix imaging of multi-layered structures with oblique array incidence is proposed. The proposed method performs a coordinate rotation in the frequency domain to extrapolate the original wave field, which is parallel to the linear array, to the virtual measurement line that is parallel to the object's surface. Then both transmission and reception wave-field extrapolation are performed in the rotated coordinate system, and an imaging condition is applied to obtain a total focused image of the multi-layered structure. The proposed method introduces an accurate coordinate transformation relation, so that it can deal with any incident angle without precision loss. It can be used in many applications, such as immersion detection of large objects with slightly curved surfaces and oblique incidence inspection of welds with an angled wedge. The simulations and experimental results show that the proposed method has high accuracy and efficiency.
The remainder of this paper is organized as follows. The method and the derivation of the proposed method are provided in Section 2. The results and discussion illustrating the method performance are presented in Section 3. Finally, the conclusions are given in Section 4.

Wave-Field Extrapolation for Multi-Layered Structures
In this subsection, the basic wave-field extrapolation is introduced, as well as the extension to multi-layered structures. Here, a linear array is assumed to detect a homogeneous medium with a sound velocity of c. The x-axis is directed along the interface, and the z-axis is perpendicular to the interface, as shown in Figure 1. In FMC, the ultrasound waves excited by the u-th element will propagate to the scatterers and the reflected waves will be received by all array elements. Let p u (x, z, t) be the ultrasound pressures reflected by the point scatterer. The propagation of the reflected pressures should satisfy the wave equation given in [18]: of layer m , respectively. Using the narrow beam assumption, the extrapolated wave fields above and below the interface m zz  are proportional: where ( , , ) By proceeding in this way for the remaining layers, the wave field at depth z within layer L can be extrapolated as: Since the focus is usually on the relative amplitudes within each layer, the amplitude scaling effect is insignificant in the image reconstruction and can reasonably be neglected [19].

Oblique Incidence Compensation
The multi-layer algorithm presented above is only applicable for horizontal multilayered structures as shown in Figure 1. However, in practical phased array inspection, such as weld inspection, a wedge between the array and the object is used to create an oblique beam and obtain excellent coverage. Consequently, a tilt angle is generated between the normal direction of the transducer surface and the z -axis. The wave-field ex- Meanwhile, p u (x, z, t) can be expanded in the wavenumber domain as in [28]: where P u (k x , z, ω) is the Fourier expansion of p u (x, z, t), k x is the wavenumber in the x direction, and ω is the angular frequency. By inserting Equation (2) into Equation (1), the Helmholtz equation can be obtained: where k z is the wavenumber component in the z direction, defined as The general solution of Equation (3) is given by [17]: where A(k x , ω)· exp(ik z z) and B(k x , ω)· exp(−ik z z) represent the upward-traveling and downward-traveling waves, respectively. We assume that the array is placed at z = 0 and all reflectors are located in the half-space z > 0. Thus, only upward-traveling waves can be recorded by the ultrasonic array. Considering the boundary condition P u (k x , z = 0, ω), which is the 2D Fourier transform of the received pressures p(x, z = 0, t) at depth z = 0, the wave field at an arbitrary depth z can be extrapolated as in [17]: where The sign function in Equation (7) guarantees that the k z value represents the waves traveling in the −z direction. Thus, the pressures p u (x, z, t) in time-space coordinates can be reconstructed by an inverse Fourier transform In the case of multi-layered structures where the sound velocity varies along the z-axis (Figure 1), Equation (6) cannot be applied directly. On the one hand, the phase factor exp(ik z z) in Equation (6) is determined by the sound velocity of the medium, and it should be modified in order to extend it to a multi-layer case. On the other hand, the wave transmissions through each interface must be considered. Theoretically, the transmission coefficient between different layers is a complex function of the incident angle and the acoustic impedances of the materials [29]. In practice, however, the transmission coefficient can be approximated as a constant based on the assumption of a narrow transducer beam [30].
A typical multi-layered structure is schematically illustrated in Figure 1. The sound velocity is different between layers, while it is homogeneous inside each layer. The layers are numbered m = 1, 2, . . . , L, where c m and d m denote the sound velocity and thickness of layer m, respectively. Using the narrow beam assumption, the extrapolated wave fields above and below the interface z = z m are proportional: where P u (k x , z − m , ω) and P u (k x , z + m , ω) represent the wave fields of the upper and lower side of the interface, respectively.
For a point scatterer within layer L, the wave field at the interface z = z L−1 is used as the boundary condition defining the solution within the layer. This gives us the solution By proceeding in this way for the remaining layers, the wave field at depth z within layer L can be extrapolated as: Since the focus is usually on the relative amplitudes within each layer, the amplitude scaling effect is insignificant in the image reconstruction and can reasonably be neglected [19].

Oblique Incidence Compensation
The multi-layer algorithm presented above is only applicable for horizontal multilayered structures as shown in Figure 1. However, in practical phased array inspection, such as weld inspection, a wedge between the array and the object is used to create an oblique beam and obtain excellent coverage. Consequently, a tilt angle is generated between the normal direction of the transducer surface and the z-axis. The wave-field extrapolation of the multi-layered structures represented by Equation (11) is no longer applicable. Therefore, we propose an oblique incidence compensation method to extrapolate the measured wave field to a virtual measurement line that is parallel to the object's surface.
Considering the case of oblique incidence, as shown in Figure 2, a wedge with angle θ is added between the linear array and the multi-layered structure. The sound velocity of the wedge is c. A classical Cartesian coordinate system is established using the center of the first array element as the origin. Another coordinate system Ox z is obtained by rotating the coordinate system Oxz by an angle θ, so that the x -axis is parallel to the surface of the multi-layered structure. The tilt of the wedge is compensated for by extrapolating the measured wave field from z = 0 to the virtual measurement line z = 0. The original element position (x, 0) is extrapolated to the new position (x new , z new ), which can be expressed in the original system as in [31]: Appl. Sci. 2021, 11, x FOR PEER REVIEW 5 of 18 trapolation of the multi-layered structures represented by Equation (11) is no longer applicable. Therefore, we propose an oblique incidence compensation method to extrapolate the measured wave field to a virtual measurement line that is parallel to the object's surface. Considering the case of oblique incidence, as shown in Figure 2, a wedge with angle  is added between the linear array and the multi-layered structure. The sound velocity of the wedge is c . A classical Cartesian coordinate system is established using the center of the first array element as the origin. Another coordinate system Ox z  is obtained by rotating the coordinate system Oxz by an angle  , so that the x -axis is parallel to the surface of the multi-layered structure. The tilt of the wedge is compensated for by extrapolating the measured wave field from 0 z  to the virtual measurement line 0 z  . The original element position ( ,0) x is extrapolated to the new position ( , ) new new xz , which can be expressed in the original system as in [31]: Thus, for the oblique incidence situation, the wave field along the virtual measurement line 0 z  can be obtained by inserting Equation (12) into Equation (8): After combining the exp( ) z ik ax and exp( ) x ik bx terms, Equation (13) is rewritten as Rearranging Equation (15), we obtain With the integrals transformation, replacing (14) is expressed as Thus, for the oblique incidence situation, the wave field along the virtual measurement line z = 0 can be obtained by inserting Equation (12) into Equation (8): After combining the exp(ik z ax) and exp(ik x bx) terms, Equation (13) is rewritten as where Rearranging Equation (15), we obtain With the integrals transformation, replacing k x with k x , Equation (14) is expressed as where and After the rotation transformation, the wavenumber-domain wave field P u (k x , z = 0, ω) given by Equation (17) can be used directly as the start point for wave-field extrapolation, giving where k z is the wavenumber along z -axis. Then, for the multi-layered structures, the wave field at depth z within layer L in an oblique incidence system can be extrapolated as:

Full-Matrix Imaging in Wavenumber Domain
In the case of FMC data, we consider the rotated coordinate system Ox z , where the linear array is rotated to z = 0. For simplicity, we assume that the u-th element emits ultrasound waves at t = 0, and the scatterer is located at (x , z ) in the L-th layer of the oblique multi-layer structure. The ultrasound waves reach the scatterer at t = τ u (x , z ), which is determined by the relative position of the scatterer and the element. At this moment, the reflected wave field from the scatterer is tightly focused and is not affected by the other scatterers. If there is an array element just above it to receive the reflected signals, an optimally focused image can be reconstructed [17]. Then, the imaging condition becomes t = τ u (x , z ) instead of t = 0. Thus, the focused image line at depth z can be reconstructed according to [28]: However, the computational burden of calculating τ u (x , z ) is very heavy for FMC data, especially in multi-layered structures. To solve this problem, we extend the received wave-field extrapolation to the transmitted wave field, to compensate for the propagation time. For the u-th active element, the transmitted pressures can be written as s u (x, z = 0, t). Note that the transmitted pressures travel along the +z direction. According to Equation (5), A(k x , ω) should be set to zero, and the extrapolated wave field from the active array element u can be expressed as Combining the wave-field extrapolation and the tilt compensation operations in oblique multi-layered structures, we can obtain the extrapolated transmitted wave field at position (x , z ) in the L-th layer, which is expressed as where Appl. Sci. 2021, 11, 10808 It should be noted that Equations (26) and (27) differ from Equations (16) and (19) due to the physical differences between forward and backward propagations. Thus, the sign function in Equations (7) and (15) should be different when it is used in the transmitted wave-field extrapolation.
After extrapolating the transmitted wave field at the point of interest, the image can be reconstructed by taking a correlation between the forward-extrapolated wave field of s u (x, z = 0, t) and the backward-extrapolated wave field of p u (x, z = 0, t), depth by depth. The correlation in the frequency domain can be written as where S u (x , z , ω) and P u (x , z , ω) are the inverse Fourier transforms of S u (k x , z , ω) and P u (k x , z , ω) over k x , respectively. The corresponding image condition is defined as After reconstructing all B-scan images, the final image is obtained by superimposing the B-scan results:f wheref (x , z ) represents the distribution of scatterers in the Ox z coordinates, and N u is the number of array elements. Figure 3 illustrates the implementation flowchart of the proposed algorithm. In summary, the algorithm consists of five main steps:

Implementation
1. Performing a 2D Fourier transform of the received pressures p u (x, z = 0, t) and transmitted pressures s u (x, z = 0, t).
2. Using Equations (17) and (25) for tilt compensation to extrapolate both received and transmitted wave fields from the coordinate system Oxz to Ox z .
3. Using Equations (21) and (24) for wave-field extrapolations to obtain both received and transmitted wave fields at depth z in layer L.
4. Image reconstruction at depth z by taking a correlation between the two wave fields using Equation (29).
5. Obtaining the full-matrix image by superimposing all B-scan images using Equation (30). It should be noted that the substitution of k x with k x is implemented by an interpolation in the k x domain. Note also that a frequency truncation is applied, as the frequency spectrum of ultrasonic signals is band-limited. After the 2D Fourier transform, the spectra P u (k x , z = 0, ω) and S u (k x , z = 0, ω) are truncated to obtain a subset corresponding to ω ∈ [2π f min , 2π f max ], where f min and f max are the lower and upper cutoff frequencies of the array, respectively. Therefore, the SNR of the image is improved due to the removal of high-frequency noise. The efficiency is also improved due to the dataset size reduction.
The proposed algorithm provides an accurate method for compensating for the oblique incidence of the linear array and achieving oblique incidence full-matrix imaging with a significant efficiency advantage compared with the time-domain TFM. Appl. Sci. 2021, 11, x FOR PEER REVIEW 8 of 18 It should be noted that the substitution of x k with x k  is implemented by an interpolation in the x k domain. Note also that a frequency truncation is applied, as the frequency spectrum of ultrasonic signals is band-limited. After the 2D Fourier transform, the spectra are the lower and upper cutoff frequencies of the array, respectively. Therefore, the SNR of the image is improved due to the removal of high-frequency noise. The efficiency is also improved due to the dataset size reduction.
The proposed algorithm provides an accurate method for compensating for the oblique incidence of the linear array and achieving oblique incidence full-matrix imaging with a significant efficiency advantage compared with the time-domain TFM.

Simulation
To verify the performance of the proposed method, a simulation for ultrasonic FMC imaging in an oblique incidence situation was conducted. The simulation was designed

Simulation
To verify the performance of the proposed method, a simulation for ultrasonic FMC imaging in an oblique incidence situation was conducted. The simulation was designed for immersion detection of a steel specimen with side-drilled holes (SDHs). The simulation was carried out on the CIVA platform (version 2020, the French Atomic Energy Commission, Paris, France). As shown in Figure 4, the specimen was made of steel, with a sound velocity of 5889 m/s, and the reflectors were a set of SDHs with diameters of 1 mm. The SDHs were arranged in a horizontal line with a step of 3 mm, and the vertical distance of the SDHs from the surface of the specimen was 10 mm. The coordinate origin was set on the upper surface of the specimen. The x-axis was parallel to the specimen's surface while the z-axis was along the depth direction. The ultrasonic transducer was set to be a linear array with the detailed parameters listed in Table 1. The array and specimen were immersed in water, with a sound velocity of 1483 m/s, and the angle between the array and the horizontal plane was 10 degrees. The distance between the center of the array and the upper surface of the specimen was 20 mm. The source signal was a Hanning windowed signal with one and a half cycles. The FMC data were sampled at 50 MHz frequency with a time range of the upper surface of the specimen. The x-axis was parallel to the specimen's surface while the z-axis was along the depth direction. The ultrasonic transducer was set to be a linear array with the detailed parameters listed in Table 1. The array and specimen were immersed in water, with a sound velocity of 1483 m/s, and the angle between the array and the horizontal plane was 10 degrees. The distance between the center of the array and the upper surface of the specimen was 20 mm. The source signal was a Hanning windowed signal with one and a half cycles. The FMC data were sampled at 50 MHz frequency with a time range of 20 us, and then imported into MATLAB (version 2018b, MathWorks, USA) for post-processing. The raw dataset size was 4000 × 64 × 64 pixels.  As well as the proposed method, the time-domain TFM was also used to process the FMC data, for comparison. The time-domain TFM utilizes the Fermat principle to calculate the refraction points at the interface, then accurately estimates the wave propagation time. The processing computer had an AMD Ryzen 5 3600 CPU and an NVIDIA GeForce GTX 970 GPU. All the algorithms used in the simulation and the experiment were accelerated by the GPU in MATLAB. Figure 5a shows the result reconstructed by the time-domain TFM. It can be seen that the tilted coordinate was well compensated. However, the shape of the SDHs in the reconstructed image was distorted, and this increased as the angle of the incident beam increased. Figure 5c shows an enlarged image of SDH A in the image reconstructed by TFM. The shape of SDH A is no longer a circle. This distortion is introduced by the heuristic approach used in the TFM, where the ray-based method used for the beam path calculation results in errors in the oblique incidence compensation [16].  As well as the proposed method, the time-domain TFM was also used to process the FMC data, for comparison. The time-domain TFM utilizes the Fermat principle to calculate the refraction points at the interface, then accurately estimates the wave propagation time. The processing computer had an AMD Ryzen 5 3600 CPU and an NVIDIA GeForce GTX 970 GPU. All the algorithms used in the simulation and the experiment were accelerated by the GPU in MATLAB. Figure 5a shows the result reconstructed by the time-domain TFM. It can be seen that the tilted coordinate was well compensated. However, the shape of the SDHs in the reconstructed image was distorted, and this increased as the angle of the incident beam increased. Figure 5c shows an enlarged image of SDH A in the image reconstructed by TFM. The shape of SDH A is no longer a circle. This distortion is introduced by the heuristic approach used in the TFM, where the ray-based method used for the beam path calculation results in errors in the oblique incidence compensation [16].
The result reconstructed by the proposed method is shown in Figures 5b,d is the local enlarged image of SDH A. Apparently, the proposed method yields a more accurate result than time-domain TFM. It can compensate well for the wave-field rotation introduced by the oblique incidence and avoid the distortion of the reflectors. The proposed method provides a more rigorous solution to the inverse problem, and the array performance exhibits less dependence on the lateral position x at any given depth z than the TFM. Furthermore, Figure 5b shows a lower background noise than the result for TFM. Specifically, the signalto-noise ratio (SNR) of Figure 5b    The result reconstructed by the proposed method is shown in Figure 5b, and Figure  5d is the local enlarged image of SDH A. Apparently, the proposed method yields a more accurate result than time-domain TFM. It can compensate well for the wave-field rotation introduced by the oblique incidence and avoid the distortion of the reflectors. The proposed method provides a more rigorous solution to the inverse problem, and the array performance exhibits less dependence on the lateral position x at any given depth z than the TFM. Furthermore, Figure 5b shows a lower background noise than the result for TFM. Specifically, the signal-to-noise ratio (SNR) of Figure 5b is 23.15, while the SNR of Figure  5a is 9.86. This means that the proposed method can achieve a better SNR than the TFM.
To quantitatively analyze the performance of the proposed method and quantify the image resolution, a dimensionless metric termed the array performance indicator (API) [11] was adopted and defined as To quantitatively analyze the performance of the proposed method and quantify the image resolution, a dimensionless metric termed the array performance indicator (API) [11] was adopted and defined as where A −6dB is the image area in which the amplitude is higher than −6 dB of its maximum value in the defined imaging area, and λ c is the wavelength at the center frequency. As well as the API, the full width at half maximum (FWHM) was introduced to represent the −6 dB lateral resolution of the results. Table 2 lists the APIs and FWHMs of five representative SDHs, where the number of the SDH represents the order of the SDH from left to right. From the table we can see that the average API of the result reconstructed by the TFM was 0.7264, while that of the result of the proposed method was 0.7070. In addition, the FWHM of the proposed method's result was 0.9738 mm, which is about 4.5% less than the result of the TFM. The table shows that the proposed method has a small advantage compared to the TFM in lateral resolution. Table 2 also lists the time costs for the image reconstruction with these two methods. The TFM took about 258.12 s, and the proposed method took only 4.31 s to reconstruct the image. The proposed method shows a huge computational efficiency advantage compared with the time-domain TFM. Therefore, the simulation results show that compared to the TFM, the proposed method can avoid the distortion of the reflectors and improve the SNR of the images while achieving a relatively low API. More importantly, it has a huge efficiency advantage over the time-domain TFM.

Experiment A
Three experiments were also conducted to validate the performance of the proposed method. The setup of experiment A is shown in Figure 6. Most of the parameters were the same as those in the simulation, except for the specimen. The specimen used in the experiment was an aluminum block, with a sound velocity of 6450 m/s. The reflectors were four SDHs aligned in an angular direction with a diameter of 1 mm. The horizontal distance of these SDHs was 3 mm, as was the vertical distance. The ultrasonic transducer used in the experiment was a linear array produced by Olympus Corporation Company in Japan. The detailed parameters of the array were the same as those used in the simulation, and are listed in Table 1. The FMC data acquisition was implemented with the EXPLORER 64/128 FMC equipment produced by The Phased Array Company (TPAC, West Chester, OH, USA), and the transmission and reception parameters were as shown in Table 3.  The experimental data were also processed by both the TFM and the proposed method. The raw dataset size was 4096 × 64 × 64 and the reconstructed image sizes were  The experimental data were also processed by both the TFM and the proposed method. The raw dataset size was 4096 × 64 × 64 and the reconstructed image sizes were 300 × 300. Figure 7a,b shows the results reconstructed by the TFM and the proposed method, respectively. Figure 7c-f shows the local enlarged images of the four SDHs in the result of the TFM from left to right, respectively. Figure 7g-j shows the same enlarged images in the result of the proposed method. The results show that both TFM and the proposed method can compensate for the wave-field rotation very well. However, as in the simulation, the experimental result of the TFM has a problem with distortion. The proposed method can reconstruct the reflectors more accurately. The normalized amplitudes of the experimental results are given in Figure 8. From this figure we can see that the background noise level of the TFM result is higher than that of the proposed method. The SNRs of the reconstruction results of the TFM and the proposed method were 5.85 and 9.19, respectively.   In order to compare the experimental results of these two methods in detail, we calculated the APIs and FWHMs of four SDHs in the reconstructed results of the two methods, as listed in Table 4. Furthermore, the time consumptions required for the two methods to reconstruct the experimental result are also listed in Table 4. It can be seen that the TFM costs about 255.89 s to reconstruct the image with an average API of 0.7302 and an FWHM of 1.0315 mm. In contrast, the proposed method only takes 4.25 s to reconstruct the same sized image with an average API of 0.6691 and an FWHM of 0.9726 mm. This shows that the proposed method can enhance the −6 dB lateral resolution by about 5.7% over the TFM and improve the computational efficiency by about 60 times.  In order to compare the experimental results of these two methods in detail, we calculated the APIs and FWHMs of four SDHs in the reconstructed results of the two methods, as listed in Table 4. Furthermore, the time consumptions required for the two methods to reconstruct the experimental result are also listed in Table 4. It can be seen that the TFM costs about 255.89 s to reconstruct the image with an average API of 0.7302 and an FWHM of 1.0315 mm. In contrast, the proposed method only takes 4.25 s to reconstruct the same sized image with an average API of 0.6691 and an FWHM of 0.9726 mm. This shows that the proposed method can enhance the −6 dB lateral resolution by about 5.7% over the TFM and improve the computational efficiency by about 60 times.

Experiment B
The second experiment was carried out on a two-layer structure, as shown in Figure 9. The upper layer was made of polymethyl methacrylate (PMMA), with a sound velocity of 2337 m/s and a thickness of 20 mm. The lower layer was the same specimen used in experiment A. The two-layer structure was immersed in water, and the distance between the array center and the upper PMMA surface was 10 mm. The tilt angle of the linear array was set to 8 degrees, and the rest of the experimental parameters remained the same as in experiment A.

Experiment B
The second experiment was carried out on a two-layer structure, as shown in Figure  9. The upper layer was made of polymethyl methacrylate (PMMA), with a sound velocity of 2337 m/s and a thickness of 20 mm. The lower layer was the same specimen used in experiment A. The two-layer structure was immersed in water, and the distance between the array center and the upper PMMA surface was 10 mm. The tilt angle of the linear array was set to 8 degrees, and the rest of the experimental parameters remained the same as in experiment A. Figure 10a shows the experimental result reconstructed by the proposed method. It can be seen that the proposed method can compensate for the tilt angle very well and can achieve a satisfactory result in the detection of two-layer structures. Figure 10b shows the normalized amplitudes of these four SDHs in the experimental result reconstructed by the proposed method. We also calculated the FWHM for every SDH and the SNR of the result, as listed in Table 5. As shown in Table 5, the proposed method achieved an average FWHM of 1.012 mm in the detection of the two-layer structure, and the SNR of the image was 26.18. This result is consistent with the results in experiment A. In fact, due to the presence of the water, this experimental object should be considered a three-layer structure. Therefore, the experimental results demonstrate the effectiveness of the proposed method for oblique incidence detection of multi-layer structures. Furthermore, Table 5 shows that the proposed method only costs 3.02 s to reconstruct an image with a size of 200 × 300 on an ordinary PC. Figure 9. The geometry of the two-layer structure in experiment B. Figure 9. The geometry of the two-layer structure in experiment B. Figure 10a shows the experimental result reconstructed by the proposed method. It can be seen that the proposed method can compensate for the tilt angle very well and can achieve a satisfactory result in the detection of two-layer structures. Figure 10b shows the normalized amplitudes of these four SDHs in the experimental result reconstructed by the proposed method. We also calculated the FWHM for every SDH and the SNR of the result, as listed in Table 5. As shown in Table 5, the proposed method achieved an average FWHM of 1.012 mm in the detection of the two-layer structure, and the SNR of the image was 26.18. This result is consistent with the results in experiment A. In fact, due to the presence of the water, this experimental object should be considered a three-layer structure. Therefore, the experimental results demonstrate the effectiveness of the proposed method for oblique incidence detection of multi-layer structures. Furthermore, Table 5 shows that the proposed method only costs 3.02 s to reconstruct an image with a size of 200 × 300 on an ordinary PC.  In order to further test the effectiveness of the proposed method in detecting defects at different tilt angles and in different locations, experiment C was carried out. The specimen used in this experiment was a standard B-type specimen for phased array ultrasonic testing, produced by Shandong Ruixiang Mould Company (China, Shandong Province, Jining City)., as shown in Figure 11a. It was made of 20# steel, with a measured sound velocity of 6012 m/s. The defects tested consisted of a set of SDHs with diameters of 1 mm arranged along a circular arc of radius 25 mm, as marked by the red rectangle in Figure  11a. In this experiment, we used a tilted wedge with an angle of 20 degrees for coupling, and the distance between the first array element center and the upper surface of the specimen was 39.28 mm, as shown in Figure 11b. The sound velocity of the wedge was 2330 m/s. The rest of the experimental parameters remained the same as in the former experiments. The raw dataset size was 4096 × 64 × 64, and the data were imported into MATLAB 2018b for post-processing.  In order to further test the effectiveness of the proposed method in detecting defects at different tilt angles and in different locations, experiment C was carried out. The specimen used in this experiment was a standard B-type specimen for phased array ultrasonic testing, produced by Shandong Ruixiang Mould Company (China, Shandong Province, Jining City)., as shown in Figure 11a. It was made of 20# steel, with a measured sound velocity of 6012 m/s. The defects tested consisted of a set of SDHs with diameters of 1 mm arranged along a circular arc of radius 25 mm, as marked by the red rectangle in Figure 11a. In this experiment, we used a tilted wedge with an angle of 20 degrees for coupling, and the distance between the first array element center and the upper surface of the specimen was 39.28 mm, as shown in Figure 11b. The sound velocity of the wedge was 2330 m/s. The rest of the experimental parameters remained the same as in the former experiments. The raw dataset size was 4096 × 64 × 64, and the data were imported into MATLAB 2018b for post-processing.
To compare the imaging results, we used the TFM and the proposed method for image reconstruction. As shown in Figure 12, both the TFM and the proposed method could achieve imaging of the measured specimen and could compensate for the tilt angle of the linear array very well. However, it can be seen from the figures that the reconstructed images from both methods have shortcomings. The result of the TFM still had a distortion problem, especially for the defects near to the upper surface of the specimen. Furthermore, the background noise of the TFM result was very large. Although the proposed method could overcome the image distortion problem and suppress the background noise, the echo energy of the defects near to the upper surface of the specimen was not as good as that of the TFM. Furthermore, neither of the two methods could reconstruct the images of the lowest defects well, because the refraction angle of the wave entering the specimen was too large, resulting in too little energy reaching the lowest defects. However, the proposed method still showed a great efficiency advantage, since the TFM cost 299.89 s to reconstruct this image with a size of 300 × 450, while the proposed method only cost 5.01 s.
In order to further test the effectiveness of the proposed method in detecting defects at different tilt angles and in different locations, experiment C was carried out. The specimen used in this experiment was a standard B-type specimen for phased array ultrasonic testing, produced by Shandong Ruixiang Mould Company (China, Shandong Province, Jining City)., as shown in Figure 11a. It was made of 20# steel, with a measured sound velocity of 6012 m/s. The defects tested consisted of a set of SDHs with diameters of 1 mm arranged along a circular arc of radius 25 mm, as marked by the red rectangle in Figure  11a. In this experiment, we used a tilted wedge with an angle of 20 degrees for coupling, and the distance between the first array element center and the upper surface of the specimen was 39.28 mm, as shown in Figure 11b. The sound velocity of the wedge was 2330 m/s. The rest of the experimental parameters remained the same as in the former experiments. The raw dataset size was 4096 × 64 × 64, and the data were imported into MATLAB 2018b for post-processing.  To compare the imaging results, we used the TFM and the proposed method for image reconstruction. As shown in Figure 12, both the TFM and the proposed method could achieve imaging of the measured specimen and could compensate for the tilt angle of the linear array very well. However, it can be seen from the figures that the reconstructed images from both methods have shortcomings. The result of the TFM still had a distortion problem, especially for the defects near to the upper surface of the specimen. Furthermore, the background noise of the TFM result was very large. Although the proposed method could overcome the image distortion problem and suppress the background noise, the echo energy of the defects near to the upper surface of the specimen was not as good as that of the TFM. Furthermore, neither of the two methods could reconstruct the images of the lowest defects well, because the refraction angle of the wave entering the specimen was too large, resulting in too little energy reaching the lowest defects. However, the proposed method still showed a great efficiency advantage, since the TFM cost 299.89 s to reconstruct this image with a size of 300 × 450, while the proposed method only cost 5.01 s.

Discussion
The experimental results show that the proposed method can achieve oblique incidence full-matrix imaging of multi-layered structures. It has a great efficiency advantage compared with the TFM and can suppress the background noise. It also shows an improvement in lateral resolution.
However, the proposed method also has some limitations. Firstly, the imaging results are not satisfactory and readable when the defect is in an area that the acoustic beam cannot easily access, as shown in Figure 12b. This is due to the existence of the refraction angle, resulting in the defect receiving very little acoustic energy. Secondly, although the proposed method can achieve full-matrix imaging at any tilt angle mathematically, the conversion of transverse and longitudinal waves in a solid medium must be considered in practical applications. When the tilt angle is too large, the longitudinal waves in a solid medium disappear. At this time, if the formula was still derived by using the velocity of longitudinal waves, incorrect results would appear. In this case, the effectiveness of imaging with other types of ultrasound waves, i.e., transverse waves, should be considered and tested in numerical simulations and laboratory experiments.

Conclusions
In this paper, a modified wavenumber method for full-matrix imaging of multi-lay-

Discussion
The experimental results show that the proposed method can achieve oblique incidence full-matrix imaging of multi-layered structures. It has a great efficiency advantage compared with the TFM and can suppress the background noise. It also shows an improvement in lateral resolution.
However, the proposed method also has some limitations. Firstly, the imaging results are not satisfactory and readable when the defect is in an area that the acoustic beam cannot easily access, as shown in Figure 12b. This is due to the existence of the refraction angle, resulting in the defect receiving very little acoustic energy. Secondly, although the proposed method can achieve full-matrix imaging at any tilt angle mathematically, the conversion of transverse and longitudinal waves in a solid medium must be considered in practical applications. When the tilt angle is too large, the longitudinal waves in a solid medium disappear. At this time, if the formula was still derived by using the velocity of longitudinal waves, incorrect results would appear. In this case, the effectiveness of imaging with other types of ultrasound waves, i.e., transverse waves, should be considered and tested in numerical simulations and laboratory experiments.

Conclusions
In this paper, a modified wavenumber method for full-matrix imaging of multi-layered structures with oblique array incidence was proposed. The proposed method utilizes a frequency-domain coordinate rotation to map the original wave field onto the virtual measurement line parallel to the object's surface, performs both transmission and reception wave-field extrapolation in the rotated coordinate system, and finally obtains a total focused image of the multi-layered structure by applying a correlation imaging condition. Since the coordinate transformation is accurately defined in the mathematics, the proposed method can deal with any incident angle without precision loss. It can be used in many applications, such as immersion detection of large objects with slightly curved surfaces and oblique incidence inspection of welds with an angled wedge. The performance of the proposed method was evaluated by FMC imaging in an oblique incidence situation conducted both in simulation and experimentally. Compared with the TFM, the proposed method provides a more mathematically rigorous solution based on the wave equation, and it can suppress the artifacts presented in the TFM. In addition, the algorithm efficiency is significantly improved; for example, it only takes about 4.25 s to reconstruct an FMC dataset with a size of 4096 × 64 × 64 on an ordinary PC. It was demonstrated that the proposed method is superior to the TFM in both accuracy and efficiency.
In the future, the efficiency of the proposed method can be further improved to meet the demands of real-time imaging and can be extended to 3D cases. Furthermore, the TFM is shown to outperform the wavenumber algorithm at large angles relative to the array [16], pointing the way for the improvement of the wavenumber algorithm.