Numerical Investigations of Ultrasonic Reverse Time Migration for Complex Cracks Near the Surface

Metal materials will produce fatigue cracks under long-term cyclic loading. It is difficult to fully image the topography of defects with the current mainstream ultrasound imaging technology. The Reverse time migration (RTM), which has extensive application in geophysical exploration, can image complex geological structures. This research introduces the RTM into the field of ultrasonic non-destructive testing. Through numerical simulation, it is concluded that the RTM based on no absorbing boundary can make the defect contour clearer, but the RTM based on the absorbing boundary can make the image contrast metric higher. In order to obtain the RTM results with high contrast metric and clear defect contours, the solution given in this study is to increase the migration aperture. Then this research discusses the influence of different cross-correlation imaging condition on the imaging quality. And compared with the prevailing total focusing methods (TFM), it shows the superiority of the RTM to the bottom opening crack (BOC) imaging. Finally, the RTM was used to complete the imaging of the curved and bifurcated cracks with complex structures. And after reducing the number of ultrasonic transducers, there are still good imaging results, which can reduce the waste of resources and energy.


I. INTRODUCTION
In engineering structures, metal materials under long-term cyclic loading will produce fatigue cracks, which will increase the possibility of failure of metal materials [1], [2]. Ultrasonic array imaging technology can image defects in the workpiece, so that the inspection results are presented in a more intuitive image [3]- [6]. However, the traditional ultrasonic imaging method can only obtain the image of the upper surface or the upper end of the defect. Because most of the echo information on the side and underside of the defect is reflected multiple times. The traditional ultrasonic imaging methods cannot handle multiple reflection waves and are difficult to effectively image defects with longitudinal The associate editor coordinating the review of this manuscript and approving it for publication was Ananya Sen Gupta . or steep complex geometric contours. Only the upper end position and lateral size of the defect can be determined. Therefore, it is necessary to develop a structural health detection method to evaluate the contour and location of the crack.
The RTM has extensive applications in geophysical exploration. The RTM can image the complex geological structures by using reflected waves [7]- [9]. Some scholars have imaged the defects in concrete based on the RTM principle. Müller et al. imaged vertical interfaces and circular voids in concrete [10]. Based on the RTM results, Beniwal et al. locate the rebar and distinguish the layering and bonding state of the rebar [11]. Qi et al. inspected the grouting quality in the connecting casing in the prefabricated building, and can image the casing in the semi-grouted state [12]. In industrial non-destructive testing, Yang et al. imaged multiple side VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ holes in the aluminium block based on the principle of the RTM, and compared it with the TFM [13]. These indicate the RTM's ability to detect internal defects, while this study is to image the BOC near the surface. Based on the RTM principle, this research perform complete imaging of curved and bifurcated BOC with complex geometries. And after reducing the number of transducers, the RTM can still perform complete imaging of complex defects. This can reduce the waste of resources and energy. The paper is organized as follows. In the next section, the principle of the RTM is introduced, and then MATLAB is used to build a 2D model of BOC in steel, and the influence of the absorbing boundary conditions and the migration aperture on the RTM result is discussed. Based on the influence of the absorbing boundary conditions and the migration aperture on the RTM results, we used COMSOL finite element software to build a 2D model of BOC with increased migration aperture in the steel, discussed the influence of different imaging conditions on the RTM results, and the RTM results with the traditional TFM results are compared. Based on the above results, we build a curved and bifurcated BOC model with complex structures in steel. Using the RTM to complete imaging of complex defects, the ability of the RTM to image complex cracks near the surface is discussed. This paper ends with a conclusion section and description of future work.

II. PRINCIPLE THE RTM
The RTM mainly includes three steps: (1) The forward propagation of the wave source field; (2) The back propagation of the receiving wave field; (3) Reconstructing the image using imaging conditions [14].
The basis of the research on RTM is the forward propagation simulation of the wave field. In the forward propagation step of the wave source field, this paper uses the finite difference time domain method to extrapolate the wave field [15]- [18]. The finite difference method uses a Cartesian coordinate grid with constant spacing to approximate the partial derivatives of the wave equation in time and space. This method has fast calculation speed and small memory usage. It can be applied not only to homogeneous media, but also to non-homogeneous media and anisotropic media. And it is based on the wave equation, so it can not only get direct waves, first and multiple reflection waves, but also converted waves, diffracted waves, etc., which can fully simulate the kinematics and dynamics characteristics of the wave field [19].

A. WAVE FIELD EXTRAPOLATION SCHEME
In 2D, the x − z plane is regarded as the vertical plane, z is the depth coordinate, and x as the horizontal coordinate. The constant density acoustic wave equation is shown in Eq. (1).
where c(x, z) is the speed of sound at point (x, z), p(x, z, t) is the sound pressure at point (x, z) at time t, s(x, z, t) is the source excitation function of time t at point (x, z). We discretize space-time with, where j and k are the grid codes of x and z in Cartesian coordinate system respectively; n is the time code after discretization.
Using the three-point operator for the second derivatives in time leads us to the extrapolation scheme, the expression is shown as, where on the right side of the equation the space and time dependencies are implicitly assumed and the partial derivatives are approximated by, The stability condition of the above difference equation is Eq. (6).
where t is the time step of calculation; d is the size of the divided grid.

B. ABSORBING BOUNDARY CONDITIONS
Limited by calculation and memory capacity and calculation time, in the forward simulation wave field of a certain medium, only a limited target area can be artificially intercepted for calculation. When intercepting a limited target area, the ultrasonic wave will reflect when it reaches the boundary. In order to eliminate the influence of boundary reflection, the reflected wave needs to be suppressed. In modern migration technology, in order to avoid the boundary reflection waves in the useful profile range, the edges are often expanded on both sides of the profile, such as perfect matching layer. As the number of edge expansion layers increases, the amount of calculation and equipment resource usage will increase. For this reason, we use absorbing boundary conditions to absorb and attenuate the wave field [20]. If the wave field value on the boundary of the stage can be accurately calculated and used as the boundary condition on the truncated boundary, no reflection will naturally occur. Similarly, if the approximate wave field value on the truncated boundary is used as the boundary condition, the reflected wave field on the truncated boundary will also be effectively suppressed.
To this end, the left absorbing boundary condition is set as shown in Eq. (7).
where x = 0 is the abscissa of the left boundary of the discussion area. The right absorbing boundary condition is set as shown in the Eq. (8).
where x = c is the abscissa of the right boundary of the discussion area. The upper absorbing boundary condition is set as shown in Eq. (9).
where y = 0 is the depth coordinate of the upper boundary of the discussion area. The lower absorbing boundary condition is set as shown in Eq. (10).
where y = d is the depth coordinate of the lower boundary of the discussion area.

C. REVERSE TIME WAVEFIELD EXTRAPOLATION
The RTM can be considered as a multi-source synthetic focusing imaging technology. As shown in Fig. 1, only three of the propagation paths of ultrasonic trapping defect information are drawn. When one of the ultrasonic transducers emits ultrasonic waves, they are received by the other three ultrasonic transducers through the three paths 1, 2, and 3 respectively. Among them, after the ultrasonic waves of path 1 and path 2 are reflected by the defect, they are respectively received by the ultrasonic transducer through the reflection of the bottom surface and the side surface, so that the information of the side surface of the defect is contained in the ultrasonic wave. The path 3 is through one reflection, so that the top information of the defect can be obtained. When performing reverse time wavefield extrapolation, the ultrasonic wave carrying defect information in the paths 1, 2, and 3 are propagated in the reverse time, and finally focus is generated at the launch position, and the defect can be imaged after using the imaging conditions. The RTM is the inverse operation of the wave equation, which is to recurse the wave field in the opposite direction of the time axis through the wave equation, starting from the maximum recorded time and ending at the zero time. In actual operation, the received echo time of the scattered signal including defects or boundaries is reversed and used as an input signal of the sound source propagating in the reverse time to make it return to the original transmitting position. This is almost the same as the wavefield extrapolation scheme in Section II.A, except that the original receiving position is transformed into the transmitting position, and the input signal of the transmitting source is the reverse time wavefield after time inversion. The position and the number of transmitting and receiving in the RTM can be set arbitrarily. But in order to match the existing ultrasonic phased array transducers, the actual transmitting and receiving mode is that one of the transducers transmits ultrasonic signals, and all ultrasonic transducers including itself receive signals. Then the ultrasonic transducer at the next position continues the above-mentioned steps until the transducers at all positions are emitted.

D. CROSS-CORRELATION IMAGING CONDITIONS
The cross-correlation imaging conditions make full use of the imaging information and can perform multiple imaging on multiple excitations. While enhancing the imaging signal, it also effectively suppresses imaging noise. The RTM for selecting cross-correlation imaging conditions includes the following three steps: First, forward the wave source field along the time direction, and save the wave field information at each time step; Secondly, reverse the received wave field along the time direction; Finally, the source wavefield and the received wavefield are cross-correlated and superimposed at each time step. The RTM result of all single wave sources are superimposed, and the RTM result I (x, z) of the area covered by the transducer can be obtained.
The traditional cross-correlation imaging condition is, where S s (x, z, t) represents the source wave field; R s (x, z, t) represents the wave field after the time reversal. VOLUME 10, 2022 The source illumination cross-correlation and the receiving wavefield illumination cross-correlation imaging conditions can suppress noise [21], [22]. The source illumination crosscorrelation and the receiving wavefield illumination crosscorrelation imaging equations are Eq. (12) and Eq. (13), respectively.

III. MIGRATION QUALITY EVALUATION
There are four main types of commonly used methods for migration quality evaluation: higher other metric, entropy metric, contrast metric, and intensity metric [23], [24]. In this research, we mainly image the BOC, only need to compare the BOC image with the whole image, the image with the largest contrast metric value corresponds to the BOC image with the highest definition. Therefore, this paper uses the contrast metric for evaluation. The contrast metric evaluation index formula is, where m, n represent the number of rows and columns of image pixels, α x i , z j = (1/mn) m i=1 n j=1 I x i , z j k and k = 2 in this study.

A. MODEL ESTABLISHMENT
In order to simplify the model and reduce the calculation time, The MATLAB is used to create a 2-D geometric model for numerical simulation. The model is shown in Fig. 1. The model to be established is the BOC in the steel. The depth of the steel is set to 20mm and the width is set to 35mm. At the bottom of the model, there is an air with a width a = 1mm and length b = 5mm to simulate the longitudinal BOC in the steel. There are 32 ultrasonic transducers arrays above the model, and the center distance between two adjacent ultrasonic transducers is 1mm.  In the model, 5MHz ultrasound is used. When 9 grids are placed in one wavelength, the sampling frequency can be set to 70MHz according to Eq. (6), and the length of time to collect is 15µs.

B. IMAGING RESULT
If the absorbing boundary conditions are not set when collecting the ultrasonic signal, input the collected data into the RTM algorithm without setting the absorption boundary, and when the imaging condition is the receiving wavefield illumination cross-correlation, the RTM result as shown in Fig. 2(a) can be obtained. Input the collected data without setting the absorbing boundary conditions into the RTM algorithm with the left and right absorption boundary to obtain the RTM result as shown in Fig. 2(b). Setting the upper absorbing boundary condition in the RTM algorithm can get the RTM result as shown in Fig. 2 If the left and right absorbing boundary conditions are set when collecting ultrasonic signals, input the collected data into the same four RTM imaging algorithms above, and the RTM results shown in Fig. 3 can be obtained. And their contrast metrics are 71.7757, 81.7765, 43.0290, and 99.1059, respectively.
In Fig. 2 and Fig. 3, the migrated results can show the image of BOC, and the top contour of the defect is relatively clear. However, the imaging of Fig. 2(a) and Fig. 2(c) is better because they can show the side profile of the defect more clearly. The reasons can be divided into two points: (1) When collecting ultrasonic signals, the ultrasonic waves with partial defect information to the left and right borders are not absorbed. This path can refer to path 2 in Fig. 1, so that more ultrasonic signals carrying side information of the defect can be received; (2) Similarly, the left and right absorbing boundary conditions are not set in the RTM algorithm, and the reflected wave information of the left and right boundaries can be fully utilized. In other images, because part of the ultrasonic waves that are reflected to the left and right borders and carry defect side information are absorbed when collecting ultrasonic signals or in the RTM algorithm, the side profile is blurred.
Although Fig. 2(a) without absorbing boundary conditions can better present the defect profile, the contrast metric of Fig. 2 and Fig. 3 with absorbing boundary conditions is greater than that of Fig. 2(a). And through comparison, it is found that there are fewer upper shadows of defects in Fig. 2(c), Fig. 2(d), Fig. 3(c) and Fig. 3(d) where there is an upper absorbing boundary condition. This is because the ultrasonic wave reflected on the upper surface of the workpiece is absorbed in the RTM imaging algorithm, and the influence of multiple reflections of the ultrasonic wave on the upper surface of the workpiece is eliminated. In addition, although the effect of Fig. 2 (a) is good, it has certain limitations. It can only detect smaller workpieces because it needs to receive the reflected waves passing through the side boundary of the workpiece.
Our goal is to obtain a higher contrast metric image and be able to identify the contour of the defect. At this point, we will focus on Fig. 3(d). The contrast metric of Fig. 3(d) is the highest in Fig. 2 and Fig. 3. However, due to the lack of information reflected on the left and right borders, some of the ultrasound signals carrying the side information of the defect will be missing, resulting in a fuzzy side of the defect. However, there is a compensation method. As shown in Fig. 1, the information of path 2 is missing, which can be compensated by adding information similar to path 1.
To this end, we established a model to verify. Compared with the model in Section IV.A, we change the model width to 70mm, and the transducers are set to the following three types: (1) The number of transducers is 32, and the center distance between two adjacent ultrasonic transducers is 1mm; (2) The number of transducers is 32, and the center distance between two adjacent ultrasonic transducers is 2mm; (3) The number of transducers is 64, and the center distance between two adjacent ultrasonic transducers is 1mm. Only the left and  Table 1. The imaging calculation time is obtained on a computer with a CPU model of 10-Core Intel Core i9.
As shown in Fig. 4, from Fig. 4(a) to Fig. 4(b) and then to Fig. 4(c), the morphological clarity of the defect outline is gradually increasing. This provides two detection schemes for workpieces that cannot obtain reflection information on the left and right sides. Solution A: Increase the center spacing of the transducers while the number of transducers remains the same. Solution B: Increase the number of transducers under the condition that the center spacing of the transducers remains unchanged. As shown in Table. 1, in the case of the same velocity model, the imaging calculation time only increases as the number of transducers increases. When parallel computing is turned on, this increase is not linear. But in terms of image quality, Solution B is better than VOLUME 10, 2022 Solution A. In Fig. 4, the fake layer can be seen in the red dotted rectangular frame. When the inverse time wavefield is extrapolated, the BOC is not set in the velocity model. The ultrasonic wave reflected by BOC is inverted and stays at the position of the fake layer. After the cross-correlation calculation, the fake layer will appear.

V. INFLUENCE OF DIFFERENT CROSS-CORRELATION IMAGING CONDITIONS ON RTM A. MODEL ESTABLISHMENT
It takes less time to acquire ultrasonic signals when using the MATLAB for numerical simulation. The seismic source is an ideal point source and the actual ultrasonic transducer has a certain size. The COMSOL can more conveniently set the size of the ultrasonic transducer and the defects with complex morphology, and the software is a finite element numerical simulation software widely recognized by the outside world, which can simulate the actual situation more realistically. Based on the above research, in this section we use the COMSOL to perform numerical simulation to obtain the ultrasonic signals.
As shown in Fig. 1, the depth of the steel is set to 20mm, the width is set to 70mm, and an air with a width a = 1mm and a length b = 5mm is set in the middle of the bottom of the model to simulate the longitudinal BOC in the steel. Both steel and air use COMSOL built-in material properties. There are 64 ultrasonic transmitting and receiving arrays above the model, the center spacing of two adjacent ultrasonic transducers is 1mm, and the width of the ultrasonic transducer is set to 0.9mm. Set the absorbing boundary conditions on the left and right boundaries of the model to simulate steel with infinite width. In the model, 5MHz ultrasound is used, the sampling frequency of the ultrasound is set to 70MHz, and the length of time to collect data is 15µs.

B. IMAGING RESULT
After increasing the migration aperture, use different imaging conditions for imaging. According to the discussion of the absorbing boundary conditions in the previous section, we set the upper, left, and right absorbing boundary conditions in the RTM algorithm. The traditional cross-correlation imaging conditions, the receiving wavefield illumination crosscorrelation imaging conditions, and the source illumination cross-correlation imaging conditions are used for imaging. The RTM results are shown in Fig. 5. And their contrast metrics are 0.1009, 12.3910, and 37139, respectively.
In Fig. 5, Fig. 5(a) has the lowest contrast metric. Fig. 5(c) has a higher contrast metric, but compared with Fig. 5(b), Fig. 5(b) has better observability. In order to improve the observability of the picture, we have made some adjustments to the picture. We denote the maximum absolute value of the pixel in the image as MI . We can adjust the color range of the three pictures in Fig. 5   As shown in Fig. 6, the observability of the adjusted Fig.  6(a) and Fig. 6(c) has been significantly improved. However, as shown in Fig. 6(b), the observability is worse than before the adjustment. And the adjusted fig. 6(c) still has a higher contrast metric than fig. 5(b). In summary, in this study, if no post-processing is performed on the picture, the receiving wavefield illumination cross-correlation imaging conditions can be selected. But in order to pursue higher imaging quality, the imaging method in Fig. 6(c) can be used. The source illumination cross-correlation imaging conditions are selected and the color range is adjusted. The color range given in this article is [0.2MI , −0.2MI ].
In addition, we use the traditional TFM technology to image the defects in this section [25], the data used is the same as the RTM, and the imaging results are shown in Fig. 7.
As shown in Fig. 7, compared with Fig. 6(c), the imaging result of the RTM is significantly better than that of the TFM. The top and bottom positions of the defect can be roughly judged from the TFM results, but the side of the defect cannot be displayed. And compared with the RTM results, the top and bottom of the defect diverge more severely, and the contours of the top and bottom of the defect cannot be accurately obtained. In the RTM results, whether it is the top and bottom of the defect, or the side of the defect, the RTM results can be well presented, and the defect morphology can  be judged more intuitively. This shows that the RTM has better detection capabilities for longitudinal BOC.

VI. IMAGING OF THE CURVED AND BIFURCATED BOC A. MODEL ESTABLISHMENT
Based on the above research in this article, we use the COMSOL to establish a curved and bifurcated BOC model as shown in Fig. 8. The morphology of crack refers to [26], [27]. The upper part of the defect in the model consists of two identical hollow arcs rotated at different angles around the point (35, 17.5). The thickness of the defect is 0.1mm, the radius of the circle where the arc is located is 2.5mm, and the central angle is 60 • , So the chord length of the arc is  also 2.5mm. The lower part of the defect is a rectangle of 0.1mm × 2.5mm.

B. IMAGING RESULT
The imaging result of RTM is shown in Fig. 9(a). Using the same data as Fig. 9(a), the imaging result of TFM is shown in Fig. 9(b).
As shown in Fig. 9, in the RTM results, the curved and bifurcated structure can be clearly seen, and the morphology of the defect can be accurately judged. The TFM results can only roughly see the approximate location of the defect. And the TFM results cannot evaluate the morphology of the defect.
In addition, the number of transducers is reduced from 64 to 32. In order to keep the migration aperture constant, we increase the center distance of adjacent transducers from 1mm to 2mm. After re-collecting the data, the imaging results of RTM and TFM are shown in Fig. 10(a) and Fig. 10(b), respectively.
As shown in Fig. 10, after reducing the number of transducers, the RTM imaging results can still clearly see the curved and bifurcated structure, while the TFM results become less observable as the number of transducers decreases. This shows that the RTM has a better detection ability for the BOC with complex structure than the TFM, and the RTM can make more accurate assessments of the location and morphology of defects. And after reducing the number of transducers by half, the RTM still has a good detection ability, which can also play a role in saving resources, which can greatly reduce calculation time. And this can use energy more effectively and rationally.

VII. CONCLUSION
In this research, an RTM based ultrasound imaging method is proposed. Based on the principle of RTM in geophysical exploration, we have successfully applied it to the field of ultrasonic phased array nondestructive testing. This article uses RTM to complete the imaging of BOC with complex structure, which is difficult to solve by TFM. And after reducing the number of transducers, there is still effective detection capability, which can reduce the waste of resources and reduce the calculation time.
Based on the application of the RTM in non-destructive testing in the research, we give the following suggestions: (1) For small-sized workpieces that can obtain reflected waves from the side of the workpiece, in order to obtain more ultrasonic waves that carry defect information, it is not suitable to add a side absorption boundary in the RTM algorithm.
(2) For large-size workpieces that cannot obtain the reflected wave from the side of the workpiece, an absorption boundary should be added to the algorithm in the RTM to simulate an infinite medium. Use limited computing resources to calculate it. By increasing the migration aperture, the ultrasonic transducer can obtain more ultrasonic waves carrying defect information, and the imaging results can be more accurate.
(3) Adding the top absorption boundary in the RTM algorithm can eliminate the influence of multiple reflection waves passing through the bottom of the workpiece and the top of the defect on the imaging, and can improve the imaging quality of the RTM.
(4) If you do not perform any post-processing on the imaging results, the cross-correlation imaging conditions can be the receiving wavefield illumination cross-correlation imaging conditions; but in order to pursue higher imaging quality, the source illumination cross-correlation imaging conditions can be selected, and the color range should be adjusted. The color range given is [0.2MI , −0.2MI ].
Future work will be oriented towards field-testing of the technique proposed in this paper. It will also include further improvement of imaging quality and imaging research on multilayer media with complex surfaces.