Extended Non-Stationary Phase-Shift Migration for Ultrasonic Imaging of Irregular Surface Component

This paper introduces an extended non-stationary phase-shift migration (ENSPS) method for ultrasound imaging of irregular surface components. This method extends non-stationary phase-shift (NSPS) technology based on the classical B-scan data to full-matrix capture (FMC) data. In addition, a new surface estimation method (SEM) is also proposed to prepare for ultrasonic imaging of irregular surface components. The method is implemented in the frequency domain, thereby ensuring the computational efficiency. In the numerical simulation, two components with different surface shapes were simulated to study the imaging ability of the method for hole defects and upper surface crack defects in irregular surface components. The simulation results show that the proposed ENSPS can detect hole defects with an interval of 0.5 mm and vertical surface cracks. In order to verify the simulation results, we set up a water-immersed ultrasonic experiment platform and obtained the ENSPS images. As a comparison, the extended phase-shift-plus-interpolation method (EPSPI) is also used to process the same data. The comparison results show that the algorithm proposed in this paper has the advantages of high resolution and low noise for the ultrasonic imaging of irregular surface components.


I. INTRODUCTION
Irregular surface component has been widely used in aerospace, petroleum tankers, automobiles, and electric power-related applications [1]. Therefore, non-destructive testing and evaluation of irregular components is of considerable significance in these industries [2]. In recent years, with the rapid development in electronic technology, ultrasonic array transducers have been widely used. Subsequently, many ultrasonic imaging algorithms in the time-domain have been developed. The synthetic aperture focusing technology (SAFT) is one of them [3], [4]. In this technique, several small apertures with large scattering angles are synthesised into one large aperture to obtain high-resolution images of a region of interest (ROI). Further, the working mode of the The associate editor coordinating the review of this manuscript and approving it for publication was Liangtian Wan . ultrasonic array is improved to FMC, i.e., record a complete time-domain data from each possible transmitter-receiver element pair, and thereby a total focusing method (TFM) with a high image resolution is developed [5]- [7].
For irregular surface components, the mainstream method is to use immersion type or specially designed wedge, and then obtain the wave propagation time of all grid pixels according to the ray tracing method, and finally use TFM for imaging. It is well known that although TFM can obtain ultrasonic images of irregular surface components, it costs a lot of computation time [8]- [11]. Other approach is to use flexible array transducers that can be closely attached to the surface, combined with optimization algorithms to obtain the best estimated surface shape, which can also solve the defect detection in complex components [12]- [14]. However, owing to the uncertainty of optimisation, an estimated rather than true surface shape is obtained and causes the image offset rotation [13]. In addition, the high cost of purchasing the transducer and the high time-consuming optimization process also need to be solved urgently.
Compared with ultrasonic imaging algorithm in the timedomain, the algorithm in the frequency domain has some advantages for imaging irregular components. The latter without increasing the calculation of ray tracing [15], [16] and use the fast Fourier transform for calculation, which improves the imaging efficiency to a certain extent. The original phase-shift migration (PSM), which originated from reflection seismology, is a frequency-domain algorithm for processing pulse-echo data. It was first proposed by Gazdag [17] in 1978. PSM extrapolates the backscattered wavefield from the recorded plane to ROI according to angular spectrum propagation theory [18], [19]. Later, in order to enable the PSM algorithm to deal with irregular formations, Gazdag and Sguazerro [20] proposed a multi-speed phase-shiftplus-interpolation method (PSPI) in 1984. Subsequently, Margrave [21] and Margrave and Ferguson [22] introduced the concept of the non-stationary linear filter operator in 1997, and divided non-stationary filters into non-stationary combined filters and non-stationary convolution filters. Then he summarized PSPI and proved that the limit form of PSPI is based on the theory of non-stationary combined linear filtering, which does not form a linear superposition of impulse response. Margrave also developed the NSPS based on the theory of non-stationary convolution filtering.
Subsequently, these improved phase-shift algorithms were applied to the ultrasound non-destructive testing field [23], and extended from the pulse-echo data to FMC. Qin et al. [24] proposed the generalized phase-shift migration (GPSM) of the PSPI method based on the theory of non-stationary combined linear filtering, and introduced the split-radium-fast-Fourier-transform (SRFFT) algorithm to speed up image reconstruction. Wu et al. [25] presented an extended phase-shift migration (EPSM) algorithm by splitting the matrix data, and demonstrated that the EPSM algorithm has higher-resolution than the TFM algorithm when the velocity changes vertically. Furthermore, Lukomski [26] extended PSPI to FMC data and realized ultrasound imaging of complex geometries. As far as the author knows, the NSPS algorithm has not been extended to FMC data in ultrasonic imaging of irregular geometries.
In our work, an extended non-stationary phase-shift (ENSPS) is proposed, which combines EPSM and non-stationary convolution filter theory. In addition, we proposed a new pulse-echo mode SEM to obtain the required velocity-depth model to achieve high-precision ultrasonic imaging of irregular geometries.
The remainder of this paper is organized as follows: in the Section II, the ENSPS theory and velocity-depth model acquisition are introduced. The simulations are conducted in Section III to evaluate the performance of the algorithm. The experiments in the Section IV. Finally, the conclusions and a discussion are given in Section V.

A. ACQUISITION OF VELOCITY-DEPTH MODEL
For irregular components, the original PSM is difficult to perform because of the multiple sound velocities at the same depth. Although some technologies can obtain the velocity-depth distribution in the geometry, such as rootmean-square velocity [27], tomographic inversion [28], and full waveform inversion [29]. However, most of the above techniques require iterative calculation for all grids in ROI, which is complex and time-consuming compared with the simple geometric operations of surface mapping technology. Although surface reconstruction based on ultrasonic imaging technology (such as TFM, etc.) has been proved to be a feasible method, the reconfiguration efficiency becomes low for large surface areas. For smooth irregular surfaces, the gap between the results of surface mapping technology and TFM imaging methods is subtle [30].
Therefore, we proposed a new SEM that is different from the pulse-echo mode in [30] and [31]. Fig. 1 depicts the new SEM. This method mainly calculates the delay of adjacent elements T i and T i+1 , and then the centre point M i of adjacent elements is regarded as the base point instead of the element T i , and establish a series of points m i on the irregular surface as the discrete points of the surface. In this study, we only considered wave propagation in two media, such as water and irregular components or wedge and irregular components. Based on this, the results of surface reconstruction can provide a reference for the velocity model. When the element spacing is relatively small, and the radius of curvature of the surface relatively big, the paths l i and l i+1 can be regarded as parallel. Thus, the coordinates of the mapping point, m i , of elements T i and T i+1 can be written as where x and z represent the horizontal and vertical coordinates of points, respectively. Parameter l represents the distance from the element to the surface. Note that the reconstructed surface is only described by N − 1 mapping points for the ultrasonic array with N elements, and these points may not be in the same vertical column as their corresponding elements. However, in the phase-shifting migration process, the velocity distribution at the same extrapolation depth must be represented by the position of each element in the ultrasonic array. Therefore, after mapping is completed, N − 1 points need to be interpolated to N to ensure that these points correspond to each array element one by one. In this study, we obtained the smooth surface information by cubic spline interpolation of mapping points. Need to declare that this new method only serves as a supplement to SEM to obtain the velocity-depth model required for imaging, but it cannot ensure the improvement of the reconstruction accuracy, at least it has not been proven in this study.

B. NONSTATIONARY PHASE-SHIFT IN IRREGULAR SURFACE GEOMETRIES
In order to obtain the imaging of irregular geometries, the PSM with constant velocity must be solved. The wavefield extrapolation process of component with horizontal surface is shown in Fig. 2. According to the exploding reflector model, the process of extrapolating the recorded data by spontaneous self-receiving to ROI can be described as [18]  where the boundary condition P r is the 2D Fourier transform of space-time wavefield p r over x and t; k x is the component of the wavevector in the x direction and ω is the angular frequency; and α (·) is the wavefield extrapolation operator, i.e., where z is an extension step. k z is the component of the wavevector in the z direction and k z is defined as a function of k x , ω, and c, i.e., where sign (·) is a sign function andĉ is half of the real speed c. In combination with Eqs. (3), (4), and (5), the recorded field is extended along the discrete Z axis piece by piece until the whole ROI can obtain the expected focus image. Before ultrasound imaging of irregular surface component, it is necessary to obtain the velocity values at each extrapolated depth from the known velocity-depth model, and then use the spatial window function to weight the wavefield. In fact, Margrave has demonstrated that the difference between the extreme form of PSPI and NSPS is the weighting order of the phase shift operator and the spatial window function in the mathematical expression, that is, PSPI adds windows in spatial dimension after phase-shifting operation, whereas NSPS is the opposite. At depth z p of each offset step, the wavefield P r may contain the aliasing information from two or more velocities, and the image defocus will be caused by using the continuation operator at a single velocity. Therefore, a different speed must be found, and then a window function must be applied to separate wavefield with different velocities. Subsequently, the Fourier transform and phase shift are performed for each velocity. The process of NSPS is depicted in Fig. 3 and can be simply expressed as where is the window function to distinguish the sound velocities at various horizontal positions at a certain depth z p and is defined as Subscript k is to the serial number of sound velocity; FT and IFT are the fast Fourier transform and inverse Fourier transform, respectively. Note that, unlike the constant velocity medium, α(·) in Eq. (6) represents the extrapolation operator in the kth medium, and the windowing process of the NSPS algorithm is performed in the spatial domain, it will not cause truncation in the frequency domain like PSPI.

C. EXTENDED TO FMC DATA
Generally, for SAFT, the ith element in the ultrasonic array transmits ultrasound to the ROI at time t = 0 and simultaneously receives the backscatter echo. Then, the operation is repeated for all array elements using B-scan. The process can be regarded as the element receiving the wave from the scatterer at t = 2τ i (x, z), where τ i (x, z) is the propagation time of the wave from the element to the scatterer. However, for FMC, the ith element transmits waves separately, and all elements receive the scatterer echo, not just the ith element. Fig. 4 depicts an FMC diagram. Note that the times-of-flight are different for every pair of transmitter-receiver is different, that is, there is delay between signals; therefore, the imaging condition (t = 0) cannot be the same as that of SAFT. Assuming the extrapolated wavefield P r k x , z p , ω of the ith transmit element, according to the principle of exploding reflector model, the horizontal focusing image line of depth z p can be obtained by extracting the wavefield p r x, z p , t at time t = 0. Therefore, the ENSPS image is reconstructed as where e jωt | t=0 = 1. Evidently, the solution to (8) cannot be obtained directly by an inverse two-dimensional Fourier transform, and τ i (x, z) of all scattering points (x, z) in the ROI must be obtained in advance. In fact, τ i (x, z) can also be considered as the delay between the transmitting and receiving wavefields, which can be rapidly obtained by solving the cross-correlation of wavefields in the frequency domain [25]. If wavefield p t x, z p , t of ith transmitter is known, the cross-correlation is where the superscript * denotes the complex conjugate, and P t k x , z p , ω and P r k x , z p , ω can be obtained from (6), but their extrapolation operators are a pair of conjugate complex numbers. Particularly noted that the exploding reflector model only applies to the pulse-echo data (single path) but not to FMC data (double path). Thus, in the extrapolation operator no longer appearsĉ but the real velocity c. After extrapolation of the transmitting wavefield from z = 0 to the scatterer, the scattering wavefield can be regarded as the wavefield emitted by the scatterer at t = 0. Therefore, t = 0 can again be regarded as a correction imaging condition of extrapolated wavefield of FMC data. In fact, Eq. (9) can be described as the reflection coefficient of the scatterer, i.e. R x, z p , ω = P r x, z p , ω P t x, z p , ω = P corr x, z p , ω P t x, z p , ω * P t x, z p , ω , where P t x, z p , ω and P r x, z p , ω represent respectively the inverse Fourier transform of the P t k x , z p , ω and P r k x , z p , ω over k x ; and the denominator represents the amplitude. However, owing to the unknown scatterer position and small denominator, the calculation may be unstable. Therefore, in our study, we window the denominator on the spatial axis to obtain [32] R s x, z p , ω = P t x, z p , ω * P r x, z p , ω Win nx P t x, z p , ω * P t x, z p , ω , (11) where Win[·] is a two-dimensional smooth window function, and nx is the number of pixels in the x-direction. The window function can be a rectangle, a triangle, a Gaussian function, etc. Although not introduced here, it is necessary to study the influence of the shape of the window function on the migration results. In this paper, a Gaussian function with nx = 8 is used.
To conclude, the information of the scatterer can be represented by R s x, z p , ω in the ENSPS algorithm. Accordingly, Eq. (8) becomes To reconstruct the ultrasonic image of irregular surface components, each depth z p in the ROI need to be calculated using (12). However, in order to obtain the final focused image, it is necessary to accumulate the focused images of all the transmitting elements, i.e., The algorithmic flow of ENSPS is depicted in Fig. 5.

III. SIMULATION
The finite-element analysis model has an ideal system interconnected with actual engineering. The main objective of the simulation is to verify the feasibility of the proposed algorithm for irregular surface components. To simplify the simulation and minimise the calculation time, three 2D geometric models were created using the COMSOL Multiphysics software. The ultrasonic array testing of carbon steel was simulated with the material with a sound velocity of 5900 m/s and the wedge PMMA with a sound velocity of 2500 m/s. The surface of Model 1 is a sinusoid, with a 2 mm amplitude and a 40 mm period, and the surface of Models 2 and 3 is a stepped-curve with a slope. In Models 1 and 2, there are eleven linear through-holes with diameter of 2 mm, and the centre distance of adjacent holes is gradually reduced from 6.5mm to 2.5mm by a step of 0.5mm. In Model 3, notches with different orientation angles (θ = 15 • , 30 • , 45 • , 60 • , 75 • and 90 • ) were set on the upper surface to simulate surface crack. The length and width of the notches were 5mm and 0.4mm, respectively. The detailed geometries of the models are depicted in Fig. 6. A 3.5 MHz, 64 element ultrasonic array is set directly above the model, with a 0.8 mm element width and a 1 mm element spacing. The excitation signal is a 3 cycle Hamming-windowed tone-burst with a centre frequency of 3.5 MHz. In the models, the left and right boundaries are set as absorbing boundaries to simulate the 3008 VOLUME 9, 2021  infinite domain. Tested the data in MATLAB 2018a running in Intel Core i7 9700F. In order to avoid the negative influence of SEM error on the focused image, we use the actual surface shape to realize the ultrasonic imaging of irregular components. For comparison, the data are also post processed by EPSPI. The reconstructed images of hole defects and crack defects in the simulation is shown in Figs. 7, 8, and 9, and all images have been VOLUME 9, 2021  up-sampled 4 times. The left and right columns in Fig. 7 represent the post-processing results of EPSPI and ENSPS, respectively. The first row depicts the results of Model 1, and the last row depicts the results of Model 2. As can be seen from Fig. 7, both ENSPS and EPSPI algorithms can obtain clear defect images. And although the groove of adjacent holes is 0.5mm (less than wavelength λ c = 1.68mm), neither of the two techniques will result in incomplete bottom contour due to the existence of defects. It can also be seen that EPSPI image has more noise than ENSPS image.
In order to quantify the accuracy of ENSPS technology, this paper introduces array performance indicator (API) defined by Holmes [7] to evaluate the defect resolution of ultrasound image. that is where A −6dB is the area within which the point spread function is greater than −6 dB down from its maximum value. Table 1 summarizes the APIs for images with hole defects. The APIs of hole #1 and hole #3 in the ENSPS images are significantly lower than that of the EPSPI images, and the average API is also slightly lower than that of the EPSPI. Combined with Figs. 8 and 9, at the θ = 15 • , both post-processing technologies have obvious artifact below the notch, which may be attributed to the reflection of wave between the upper surface and notch. As the orientation angle θ increases, the artifact gradually disappears. at the θ = 90 • , the amplitude of the notch is very weak, but the outline can still be distinguished.
In general, the ENSPS algorithm has superior signalto-noise ratio and can present the location and shape of crack-like defects clearly and accurately. Measurement results comparison of the length and orientation for surface   notch are summarized in Table 2. Unfortunately, compared with EPSPI, the proposed ENSPS method does not run faster. The processing time of single transmitted data of ENSPS method and EPSPI method is about 0.12s and 0.15s respectively, and the total running time under the condition of 64 elements is about 7.78s and 9.60s respectively.

IV. EXPERIMENT A. EXPERIMENTAL SETUP
To evaluate the results of the simulation, we designed and built an experimental system, as shown in Fig. 10. A multiplexer (MUX-128D-E, Japan Probe) with 128 independent channels was connected to the ultrasonic pulser-receiver (JPR-600C, Japan Probe), which was used to control the delay rule of the ultrasonic array and display the ultrasonic signal on the PC terminal. A commercial ultrasound array VOLUME 9, 2021  (Guangzhou Doppler Electronics Technology Co., Ltd.) was used for data acquisition. Detailed parameters of the array are listed in Table 3 The ultrasonic data of irregular surface components were obtained by immersion test. Fig. 11 displays the physical photo of samples.

B. EXPERIMENTAL RESULTS
In the experimental device, it is difficult to accurately determine the relative position of parts and sensors. Hence, we used the new SEM as described in section A of Part II to reconstruct the surface shape of an irregular surface component. The results of surface estimation for samples 1 and 2 are depicted in Fig. 12. The mapping results of the array and the B-scan image of the detection data are included in the figure. It can be seen that the SEM results of sample 1 show more complete surface than B-scan images, which are mainly observed on the back of the ridge (within the red dotted line). However, this phenomenon did not appear in sample 2 with relatively smooth surface.
Further, the results of the samples 1 and 2 are shown in Fig. 13. In Figs. 13 (a) and (c), the EPSPI images are clearly shown, but the hyperbolic mode appears in SDHs 3012 VOLUME 9, 2021     Table 4. It can be seen from Table 4 that the APIs for the ENSPS are lower than for the EPSPI. Fig. 14 shows results of post-processing experimental data from sample 3. In this figure, the valid elements is unable to receive the echo on the right side of the surface due to sample 3 at an inclination of 15 • , resulting in a missing surface image, but the missing parts can be recovered by increasing the scanning length. In addition, at θ = 15 • , the same artifact appears below the notch as in the simulation. It can be seen from EPSPI images that hyperbolic mode (marked in red) appears at the notch tip at θ = 30 • and 45 • , resulting in the measured results being larger than the actual values. at θ = 60 • , incomplete notch appeared in the EPSPI image but were still completely visible in the ENSPS image. When θ is greater than 75 • , the notches are not fully focused in the ENSPS and EPSPI images because the notches echo (except for the tip) are reflected by the left notch.. However, when the distance between the notches is increased (or the orientation θ is reduced), the complete focused image as in the simulation may be obtained again. The measurement results of notch parameters for sample 3 are shown in Table 5. With the same hardware parameters, the ENSPS algorithm proposed in this paper has a higher resolution and image signal-to-noise ratio than EPSPI for surface notchs and hole defects.

V. CONCLUSION AND DISCUSSION
With the rapid development in industrial technology, irregular surface components are being widely used in the aerospace, petroleum, automotive, and other industries. Therefore, non-destructive testing and evaluation methods for irregular surface components have become the focus of research. In this study, we proposed an extended non-stationary phase-shift migration (ENSPS) method for ultrasound imaging of irregular surface component. Different from PSPI, the truncation in the space of the NSPS algorithm occurs before the phase-shift migration. Moreover, to address the limitation of the unknown surface of irregular surface components, a new SEM is introduced to obtain velocitydepth model before ENSPS post-processing.
The results of the numerical simulation verified the feasibility of the proposed algorithm for ultrasonic imaging of irregular surface components. ENSPS is a high-precision imaging technology that can better evaluate hole defects with a gap of 0.5 mm and upper surface notches with a 90 • orientation angle. It has obvious advantages over EPSPI in terms of signal-to-noise ratio and resolution. In addition, we also set up an immersion type ultrasonic experimental platform for irregular surface components and obtained the expected imaging results.
However, it should be noted that although the ENSPS algorithm relies on fast Fourier transform, the process of layer-by-layer recursion and all B-scan repetition still exists. Its computational performance is not significantly superior to EPSPI, and it is still difficult to meet the requirements of real-time imaging. Therefore, it may need to be combined with GPU parallel processing for its calculation ability to be improved. The existence of these limitations leaves scope for future research.