A high-quality imaging method for full matrix data based on frequency-domain reverse time migration

Abstract

The present invention provides a high-quality imaging method for full matrix data based on frequency-domain reverse time migration. The excitation signal, sound velocity distribution, and full matrix data are respectively preprocessed to obtain an excitation sound source matrix, an impedance matrix, and a received sound source matrix. Using the above preprocessing results as inputs, LU decomposition processing is performed on the impedance matrix to solve for the transmitted wave field and the received wave field. Finally, an imaging result is obtained using an imaging condition. This imaging method utilizes the characteristic that the impedance matrix is fixed at the same frequency. Through LU decomposition, the sound pressure fields generated by multiple excitations at the same frequency can be obtained at one time, greatly saving calculation time and improving calculation efficiency. Based on the full wave equation for wave field forward modeling, there is little simplification of wave field propagation, which is suitable for measuring defects in any complex sound velocity distribution and can ensure high quality of the obtained imaging result.

Description

Technical Field

[0001] The present invention belongs to the technical field of detection methods, and particularly relates to a high-quality imaging method for full matrix data based on frequency-domain reverse time migration. Background Art

[0002] Curved surface components are widely used in various fields. Among them, shafts and camshafts are core components in drive systems. In addition, due to dynamic requirements, more complex curved surface components widely exist in fields such as automotive, energy, and aerospace. However, due to the complex manufacturing process and harsh service environment of such components, the components are prone to holes and crack defects, which greatly affect the mechanical properties of the components, may cause premature failure of the components, and even result in serious accidents. Therefore, high-quality detection of internal defects of curved surface components is an urgent problem, but the large component size and the existence of complex surfaces also pose greater challenges to detection technologies.

[0003] Phased array ultrasonic imaging has achieved a prominent position in non-destructive testing and can effectively characterize invisible defects inside components. Full matrix capture (FMC) is an important capture mode in phased array ultrasound, where each array element acts as both a receiver and a transmitter, and it can capture the complete time-domain signals of each pair of transmitting and receiving array elements. Existing imaging methods for full matrix data have shown their ability to detect components with flat surfaces, whose surfaces are parallel or inclined to the surface of the phased array, but these methods cannot be directly extended to curved surfaces. To measure curved surface parts, scholars have introduced flexible array sensors to adapt to curved surfaces, but these technologies increase the measurement cost. On the other hand, non-contact measurement methods through water immersion are more widely adopted, and the present invention mainly processes the full matrix data obtained in this case.

[0004] The water immersion measurement method poses higher requirements for the post-processing of full matrix data. Currently, there are several ultrasonic imaging methods for curved surface interface structures. Typically, the full focus method (TFM) adapts to curved interfaces by introducing a ray tracing method. In addition, Sutcliffe proposed a virtual source method for high-resolution imaging of defects in a double-layer medium with a complex interface, where the Fermat principle is still used to search for the incident points of rays on the reflection interface. However, this method does not consider multiple echoes from defects, which reduces the resolution of the imaging results and leads to artifacts in the results. The TFM method is a standard method in full matrix imaging methods and has shown to be efficient in previous studies, but when the sound speed model in the measurement area becomes complex, its efficiency will be greatly reduced. In the above studies, the sound speed models are all double-layered to ensure imaging efficiency. However, the computational complexity of ray tracing increases exponentially with the increase in the number of layers in the sound speed model, and it becomes unacceptable when the number of layers increases to four.

[0005] Frequency-domain wavefield extrapolation can also be used for the case of curved interfaces. Lukomski extended the Phase-shift-plus-interpolation (PSPI) method to full matrix data for measuring curved parts. In the PSPI method, the sound speed field in the interface region is approximated as a combination of two uniform sound speed fields. Phase-shift migration is performed separately according to the two sound speeds, and then superimposed to correct the extrapolated wavefield. Yang et al. also considered the multiple reflections from the bottom of the measured part to improve the previous imaging method. In addition, Chang et al. introduced the Non-stationary Phase Shift (NSPS) method to compensate for the influence of horizontal sound speed variations on downward extrapolation, which can eliminate the discontinuity of the reconstructed wavefield compared with PSPI. However, frequency-domain wavefield extrapolation simplifies the wave equation to a unidirectional form and approximates the sound speed model, and the imaging quality deteriorates significantly with the increase in the number of curved interfaces.

[0006] Reverse Time Migration (RTM) was first developed in the field of seismic imaging to reconstruct the wavefield in the measurement area using the full wave equation. RTM has the advantages of high imaging accuracy, no dip limitation, and strong adaptability to any complex velocity model. In recent years, with the rapid development of computer technology, the huge computational cost that restricted its development has been alleviated to a great extent, and it has been introduced into FMC ultrasonic imaging in some scenarios. Rao et al. used Elastic Reverse Time Migration (ERTM) to image irregularly shaped grooves in metal blocks. Liu et al. successfully introduced ERTM into the ultrasonic imaging of concrete structures. However, ERTM is calculated in the time domain and requires iteration at each sampling time. In addition, for FMC data, it is necessary to separately implement the forward and backward propagation of the excitation and received signals for each excitation. These make the method have high computational complexity and long computational time. Summary of the Invention

[0007] To solve the problems existing in the prior art, the present invention provides a frequency-domain reverse time migration imaging method for ultrasonic full matrix data, which can greatly improve the imaging quality of curved parts while controlling the computational time. The present invention extrapolates the wavefield in the selected frequency range, and the wavefields of different excitations at a certain frequency can be calculated in one iteration. Compared with ERTM, the computational complexity is effectively reduced.

[0008] A high-quality imaging method for full matrix data based on frequency-domain reverse time migration includes the following steps:

[0009] (1) Input the excitation signal, sound speed distribution, and full matrix data;

[0010] (2) Preprocess the excitation signal and full matrix data respectively to obtain the excitation sound source matrix and the received sound source matrix;

[0011] (3) Discretize the sound speed distribution and calculate the impedance matrix;

[0012] (4) According to the obtained excitation sound source matrix and receiving sound source matrix, use the LU decomposition of the impedance matrix to solve for the transmitted wave field and the received wave field respectively;

[0013] (5) Process the transmitted wave field and the received wave field using the imaging condition to obtain the initial image, and post-process the initial image to obtain the imaging result.

[0014] Among them, the excitation signal in step (1) is expressed as E(x r , x s , t), the full matrix data is expressed as D(x r , x s , t), and the sound speed distribution is expressed as c(x, z);

[0015] In step (2), the excitation sound source matrix and the receiving sound source matrix respectively include the excitation sound source matrix [E1(ω), E2(ω), …, E min , ω max corresponding to each discrete frequency ω within the selected frequency band range [ω n (ω)] and the receiving sound source matrix Among them, n represents the number of excitation times (i.e., the number of elements of the phased array probe, and one element is excited once). The selected frequency band range [ω min , ω max is selected according to the characteristics of the excitation signal and the receiving signal (full matrix data); within the selected frequency band range, the frequency domain step size between every two adjacent discrete frequencies is determined by the sampling frequency, denoted as dω.

[0016] Preferably, the following preprocessing is performed on the full matrix data:

[0017] First, apply a window function to the full matrix data to filter out the direct wave, perform a Fourier transform on the full matrix data with the direct wave filtered in the time dimension, and then take the complex conjugate of the data in the obtained frequency domain to achieve the effect of time reversal, and organize to obtain the receiving sound source matrix.

[0018] Full matrix capture (FMC) captures the time domain signals between all transmit-receive pairs in the phased array probe. For an n-element phased array probe, the FMC data contains n 2 A-scan signals. D(x r , x s , t) represents the signal transmitted from the s-th element and received by the r-th element. As Figure 2 shown, the ultrasonic wave is emitted from the element position (x s , 0), reflected by the scattering point located at (x, z), and then received by the element located at (x r,0) is received by the sensor at that location. The Fourier transform in the time domain gives the frequency-domain representation of the full matrix data, as shown in the following equation:

[0019] D(x r ,x s ,ω)=∫D(x r ,x s ,t)e jkωt dt

[0020] where j represents the imaginary unit; k represents the horizontal wavenumber; ω represents the frequency; and t represents the time.

[0021] According to the properties of the Fourier transform, it is known that the Fourier transform of the inverse-time signal is equivalent to taking the conjugate of the original signal. Therefore, the inverse-time signal in the frequency domain can be written as the following equation:

[0022]

[0023] where * represents the complex conjugate.

[0024] Preferably, the excitation signal is preprocessed as follows:

[0025] The Fourier transform processing of the excitation signal in the time domain is performed to obtain the excitation sound source matrix.

[0026] Since the calculation region is finite, if no conditions are added, false reflections will occur at its boundaries. Preferably, when discretizing the sound speed distribution, the impedance matrix is optimized by adding a perfectly matched layer around the discretization model.

[0027] where the attenuation γ(x) of the perfectly matched layer (PML) is defined as follows:

[0028]

[0029] where d is the width of the PML layer; x′ is the coordinate in the local coordinate system in the PML, and the origin of the local coordinate system is located outside the model boundary; c PML is a defined scalar.

[0030] Preferably, in step (3), the calculation formula of the impedance matrix A(ω) is as follows:

[0031]

[0032] where represents the Laplace operator obtained after discretization; ω represents the frequency; and c represents the sound speed.

[0033] Preferably, the 9-point finite difference method is used to discretize the sound speed distribution.

[0034] Preferably, the specific steps of step (4) are as follows:

[0035] 4.1 Select any un-traversed frequency within the selected frequency band as the current frequency. Using the excitation sound source matrix, receiving sound source matrix, and impedance matrix corresponding to the current frequency as inputs, perform LU decomposition on the impedance matrix to solve for the excitation sound pressure field and receiving sound pressure field corresponding to the current frequency;

[0036] 4.2 Traverse all frequencies within the selected frequency band to obtain the transmitted wave field and the received wave field.

[0037] As a further preference, in step 4.1, the calculation formula for solving the sound pressure field generated by the sound source matrix is as follows:

[0038] L(ω)U(ω)[P1(ω), P2(ω), …, P n (ω)] = [S1(ω), S2(ω), …, S n (ω)]

[0039] where n represents the number of excitation times; L(ω) and U(ω) represent the two factors obtained after performing LU decomposition on the impedance matrix A(ω); [S1(ω), S2(ω), …, S n (ω)] represents the excitation sound source matrix or the receiving sound source matrix; [P1(ω), P2(ω), …, P n (ω)] represents the excitation sound pressure field or the receiving sound pressure field.

[0040] The derivation process of the above calculation formula for the sound pressure field is as follows:

[0041] The two-dimensional wave equation in the frequency domain is:

[0042]

[0043] where x and z are the two-dimensional space coordinates; ω is the frequency; p(x, z, ω) is the sound pressure field; s(x, z, ω) is the sound source field; c(x, z) is the sound speed distribution.

[0044] The sound speed distribution is as Figure 2 (a) shows that the 9-point finite difference method can introduce formula 1) to discretize the spatial domain, and the sound speed is discretized into l = N x ×N y grid points, as Figure 2 (b) shows; where N x and N y are the number of grid points on the x and y axes respectively. The discrete sound wave equation can be written in matrix form as follows:

[0045] A(ω)P(ω) = S(ω) 2)

[0046] Wherein, A(ω) is the impedance matrix; P(ω) is the discrete sound pressure field; S(ω) is the discrete sound source vector. A(ω) can also be interpreted as the forward propagation operator, as shown in the formula:

[0047] Here, P(ω) and S(ω) are column vectors of l×1, and A(ω) is an l×l matrix.

[0048] For forward propagation, S(ω) is given by the emitted wave, A(ω) depends on the frequency (ω) and the acoustic properties of the medium. Therefore, only P(ω) is unknown and can be solved by Equation (2). The LU decomposition of the impedance matrix A(ω) can be used to solve P(ω) in Equation (2). The factors of the LU decomposition of the impedance matrix A(ω) corresponding to the same frequency can be reused to solve the sound pressure field generated by any new sound source vector S(ω) corresponding to this frequency. Therefore, the sound pressure fields excited multiple times at the same frequency can be calculated efficiently at one time. This forward propagation solution method is very suitable for full matrix data.

[0049] For the case of multiple excitations at the same frequency, it can be arranged as Equation (3)

[0050] L(ω)U(ω)[P 1 (ω), P2(ω), …, P n (ω)] = [S 1 (ω), S2(ω), …, S n (ω)] (3)

[0051] As a further preference, when traversing the frequencies, it is carried out in ascending order.

[0052] Traverse the frequencies within the selected frequency band range [ω min , ω max in ascending order. After completing the traversal of the current frequency ω, judge whether the current frequency ω is less than ω max ; if so, it means there are still un-traversed frequencies, then perform the traversal of the next frequency (ω + dω) until the current judgment frequency is equal to ω max , indicating that all frequencies have been traversed, then output the solved emitted wave field S i (x, z, ω) and the received wave field R i (x, z, ω).

[0053] As a preference, in step (5), the expression of the imaging condition is:

[0054]

[0055] Wherein, I(x, z) represents the initial image; n represents the number of excitations; N ωdenotes the total number of discrete frequencies in the selected frequency band; S i (x, z, ω) represents the transmitted wave field; R i (x, z, ω) represents the received wave field; * denotes taking the complex conjugate.

[0056] Preferably, the following post - processing is performed on the initial image:

[0057] First, the initial image is sharpened by Laplace filtering, then the Hilbert transform is performed on the sharpened initial image to obtain the envelope of the image, and finally, median filtering is used to filter out salt - and - pepper noise to obtain the imaging result.

[0058] The high - quality imaging method of full - matrix data based on frequency - domain reverse - time migration of the present invention can perform high - quality imaging on internal defects of curved parts. This method can be summarized into three steps: forward extrapolation of the source wave field, backward extrapolation of the received wave field, and implementation of the imaging condition. Among them, the inputs of the imaging method include the excitation signal, sound - speed distribution, and full - matrix data. The above three are respectively pre - processed and then the wave field is solved. E i (ω) and respectively represent the excitation vector and the reverse - time received vector corresponding to the frequency ω at the i - th excitation, i ∈ [1, n]; n is the number of array elements in the phased - array probe, equal to the number of excitations. Since the impedance matrix A(ω) is fixed at a certain frequency, the excitation vectors and reverse - time received vectors at different excitations can be respectively organized into a source matrix, that is, the excitation source matrix and the received source matrix. Therefore, the sound - pressure fields generated by different excitations can be calculated simultaneously in one iteration, which greatly improves the calculation efficiency. The imaging method of the present invention adopts a frequency - by - frequency iterative method to calculate the sound - pressure fields corresponding to each discrete frequency within the selected frequency - band range [ω min , ω max , where the frequency - band range is determined by the full - matrix data, and the frequency - domain discrete step dω is determined by the sampling frequency. S i (x, z, ω) and R i (x, z, ω) are respectively the transmitted wave field and the received wave field corresponding to the frequency ω at the i - th excitation. Taking the transmitted wave field and the received wave field as inputs, the imaging result is obtained according to the imaging condition.

[0059] Compared with the prior art, the beneficial effects of the present invention are:

[0060] The high-quality imaging method of full matrix data based on frequency domain reverse time migration of the present invention preprocesses the excitation signal, sound velocity distribution and full matrix data respectively to obtain the excitation sound source matrix, impedance matrix and receiving sound source matrix; the above preprocessing results are used as input, the impedance matrix is ​​LU decomposed, and the emission wave field and the receiving wave field are obtained, and finally the imaging result is obtained by using the imaging conditions. The imaging method uses the fixed characteristic of the impedance matrix at the same frequency, and the sound pressure field generated by multiple excitations at the same frequency can be obtained at one time through LU decomposition, which greatly saves the calculation time and improves the calculation efficiency. The full wave equation (two-dimensional wave equation in the frequency domain) is used as the basis for wave field forward modeling, which simplifies the wave field propagation less, is suitable for the measurement of defects in any complex sound velocity distribution, and can ensure the high quality of the imaging result. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 is a flow chart of an imaging method according to an embodiment of the present invention;

[0062] Figure 2 (a) is a schematic diagram of the sound velocity distribution; (b) is a schematic diagram of the sound velocity distribution after discretization;

[0063] Figure 3 (a) is the preprocessing flow chart of the full matrix data; (b) is the flow chart of obtaining the imaging results from the transmitted wave field and the received wave field;

[0064] Figure 4 (a) is a schematic diagram of hole defect detection in copper parts; (b) is a schematic diagram of hole defect detection in aluminum parts; (c) is a schematic diagram of crack defect detection in aluminum parts;

[0065] Figure 5 (a) is the imaging result of hole defects in copper parts using the TFM-based method; (b) is the imaging result of hole defects in copper parts using the PSM-based method; (c) is the imaging result of hole defects in copper parts using the FD-RTM method; (d) is the imaging result of hole defects in aluminum parts using the TFM-based method; (e) is the imaging result of hole defects in aluminum parts using the PSM-based method; (f) is the imaging result of hole defects in aluminum parts using the FD-RTM method; (g) is the imaging result of crack defects in aluminum parts using the TFM-based method; (h) is the imaging result of crack defects in aluminum parts using the PSM-based method; (i) is the imaging result of crack defects in aluminum parts using the FD-RTM method;

[0066] Figure 6(a) shows the pixel amplitude curves on the marked straight line in the imaging results of hole defects in copper parts by different imaging methods; (b) shows the pixel amplitude curves on the marked straight line in the imaging results of hole defects in aluminum parts by different imaging methods; (c) shows the pixel amplitude curves on the marked straight line in the imaging results of crack defects in aluminum parts by different imaging methods. Detailed implementation manners

[0067] As Figure 1 shown, a high-quality imaging method for full matrix data based on frequency-domain reverse time migration includes:

[0068] 1. Preprocess the input excitation signal E(x r , x s , t), full matrix data D(x r , x s , t) and sound speed distribution c(x, z) respectively to obtain an excitation sound source matrix, a receiving sound source matrix and an impedance matrix.

[0069] 1.1 Preprocessing of the excitation signal

[0070] Perform Fourier transform processing on the excitation signal in the time domain to obtain an excitation sound source matrix.

[0071] 1.2 Preprocessing of the full matrix data

[0072] As Figure 3 shown in (a), first use a window function to act on the full matrix data to filter out the direct wave, perform Fourier transform on the full matrix data after filtering out the direct wave in the time dimension, and then take the complex conjugate of the data in the obtained frequency domain to achieve the effect of reverse time, and organize to obtain a receiving sound source matrix.

[0073] Full matrix capture (FMC) captures the time-domain signals between all transmit-receive pairs in a phased array probe. For an n-element phased array probe, the FMC data contains n 2 A-scan signals. D(x r , x s , t) represents the signal transmitted from the s-th element and received by the r-th element. As Figure 2 shown, ultrasonic waves are emitted from the element position (x s , 0), reflected by the scatterer located at (x, z), and received by the sensor located at (x r , 0). The Fourier transform in the time domain gives the frequency-domain representation of the full matrix data, as shown in the following formula:

[0074] D(x r , x s , ω) = ∫D(x r , x s , t)ejkωt dt

[0075] Among them, j represents the imaginary unit; k represents the horizontal wavenumber; ω represents the frequency; t represents the time.

[0076] According to the properties of the Fourier transform, it is known that the Fourier transform of the time-reversed signal is equivalent to taking the conjugate of the original signal. Therefore, the time-reversed signal in the frequency domain can be written as the following formula:

[0077]

[0078] Among them, * represents the complex conjugate.

[0079] The obtained excitation sound source matrix and receiving sound source matrix respectively include the excitation sound source matrices [E1(ω), E2(ω), …, E min , ω max corresponding to each discrete frequency ω within the selected frequency band range [ω n (ω)] and the receiving sound source matrix where the subscript n represents the excitation times. The selected frequency band range [ω min , ω max is selected according to the characteristics of the excitation signal and the receiving signal (full matrix data); within the selected frequency band range, the frequency domain discrete step size between every two adjacent discrete frequencies is determined by the sampling frequency, denoted as dω.

[0080] 1.3. Preprocessing of the sound speed distribution

[0081] The sound speed distribution is discretized by using the 9-point finite difference method, and the impedance matrix A(ω) is calculated by using the following formula:

[0082]

[0083] Among them, represents the Laplace operator obtained after discretization; ω represents the frequency; c represents the sound speed.

[0084] Since the calculation region is finite, if no conditions are added, false reflections will occur at its boundaries. When discretizing the sound speed distribution, the impedance matrix is optimized by adding a perfectly matched layer around the discretization model.

[0085] Among them, the attenuation situation γ(x) of the perfectly matched layer (PML) is defined as follows:

[0086]

[0087] Among them, d is the width of the PML layer, x′ is the coordinate in the local coordinate system in the PML, the origin of the local coordinate system is located outside the model boundary, and c PML is the defined scalar.

[0088] 2. Based on the obtained excitation sound source matrix and reception sound source matrix, LU decomposition of the impedance matrix is adopted to separately solve for the transmitted wave field and the reception wave field.

[0089] The specific steps include:

[0090] 2.1. Select any un-traversed frequency ω within the selected frequency band range [ω min , ω max as the current frequency. Using the excitation sound source matrix [E1(ω), E2(ω), …, E n (ω)], the reception sound source matrix and the impedance matrix A(ω) corresponding to the current frequency as inputs, perform LU decomposition on the impedance matrix A(ω) to solve for the excitation sound pressure field and the reception sound pressure field corresponding to the current frequency;

[0091] Among them, the derivation process of the solution formula for the sound pressure field is as follows:

[0092] The two-dimensional wave equation in the frequency domain is:

[0093]

[0094] Among them, x and z are respectively two-dimensional space coordinates; ω is the frequency; p(x, z, ω) is the sound pressure field; s(x, z, ω) is the sound source field; c(x, z) is the sound speed distribution.

[0095] The sound speed distribution is as shown in Figure 2 (a). The 9-point finite difference method can introduce formula 1) to discretize the spatial domain, and the sound speed is discretized into l = N x ×N y grid points, as shown in Figure 2 (b); among them, N x and N y are respectively the number of grid points on the x and y axes. The discrete sound wave equation can be written in matrix form as follows:

[0096] A(ω)P(ω) = S(ω) 2)

[0097] Among them, A(ω) is the impedance matrix; P(ω) is the discrete sound pressure field; S(ω) is the discrete sound source vector. A(ω) can also be interpreted as the forward propagation operator, as shown in the formula:

[0098] Here, P(ω) and S(ω) are column vectors of l×1, and A(ω) is an l×l matrix.

[0099] For forward propagation, S(ω) is given by the transmitted wave, A(ω) depends on the frequency (ω) and the acoustic properties of the medium. So only P(ω) is unknown and can be solved from Equation (2). The LU decomposition of the impedance matrix A(ω) can be used to solve for P(ω) in Equation (2); since the impedance matrix A(ω) corresponding to the same frequency is fixed, therefore, the factors of the LU decomposition of the impedance matrix A(ω) corresponding to the same frequency ω can be reused to solve for the acoustic pressure field generated by any new sound source vector S(ω) corresponding to that frequency ω. Therefore, the acoustic pressure fields excited multiple times at the same frequency can be calculated efficiently at one time, and this forward propagation solution method is very suitable for full matrix data.

[0100] For the case of multiple excitations at the same frequency, it can be arranged as Equation (3).

[0101] L(ω)U(ω)[P 1 (ω), P2(ω), …, P n (ω)] = [S 1 (ω), S2(ω), …, S n (ω)] (3)

[0102] where n represents the number of excitations (i.e., the number of elements of the phased array probe); L(ω) and U(ω) represent the two factors obtained after the LU decomposition of the impedance matrix A(ω); [S1(ω), S2(ω), …, S n (ω)] represents the excitation sound source matrix or the receiving sound source matrix; [P1(ω), P2(ω), …, P n (ω)] represents the excitation acoustic pressure field or the receiving acoustic pressure field.

[0103] 2.2. Traverse all frequencies within the selected frequency band to obtain the transmitted wave field and the received wave field.

[0104] Traverse the frequencies within the selected frequency band [ω min , ω max in ascending order. After completing the traversal of the current frequency ω, judge whether the current frequency ω is less than ω max ; if so, it means there are still un-traversed frequencies, then perform the traversal of the next frequency (ω + dω) until the current judgment frequency is equal to ω max , indicating that all frequencies have been traversed, then output the solved transmitted wave field S i (x, z, ω) and the received wave field R i (x, z, ω).

[0105] 3. As Figure 3As shown in (b), the transmitted wavefield and the received wavefield are processed using an imaging condition to obtain an initial image. First, the initial image is sharpened using Laplace filtering, then the Hilbert transform is applied to the sharpened initial image to obtain the envelope of the image, and finally, median filtering is used to filter out salt-and-pepper noise to obtain the imaging result.

[0106] Among them, the expression of the imaging condition for the initial image I(x, z) is:

[0107]

[0108] In the above formula, n represents the number of excitation times; N ω represents the total number of discrete frequencies in the selected frequency band;

[0109] S i (x, z, ω) represents the transmitted wavefield; R i (x, z, ω) represents the received wavefield; * represents taking the complex conjugate.

[0110] Detection experiment

[0111] An experiment is conducted to verify the effectiveness of the method of this embodiment for imaging defects in curved surface components. As Figure 4 shown in (a) and (b), first, the hole defects in copper and aluminum components with a radius of 80 mm are measured. There are four hole defects with a diameter of 1 mm in the components. Then, Figure 4 the V-shaped crack defect in the aluminum component with a radius of 80 mm is measured in (c). In the experiment, a 64-element phased array probe (Shantou Ultrasonic) with a center frequency of 2.5 MHz and an element spacing of 0.75 mm is used, and an FMC data acquisition is realized in cooperation with a data acquisition card 64 / 64OEM-PA (AOS.Ltd, America). The sampling frequency is 50 MHz, and the time range is 60 μs. The above components are all immersed in water during measurement. Since the phased array is not a water-immersion probe, a wedge made of polyphenylene sulfide (PPS) is used to separate the probe from the water. The sound velocity of the wedge is 2200 m / s. The bottom surface of the wedge is a concave surface for focusing more ultrasonic waves into the measured component. Water is used as a couplant to fill the gap between the wedge and the component, but does not reach the top surface of the wedge.

[0112] Three imaging methods are respectively used to process the experimental data. The imaging results are as Figure 5 shown, where the comparison methods are marked as TFM-based and PSM-based in the figure, and the method of this embodiment is denoted as FD-RTM. The frequency bands selected for the three experiments are all from 1 MHz to 4 MHz.

[0113] For the detection of hole defects in copper parts, the TFM-based method and FD-RTM can clearly image all defects, asFigure 5 as shown in (a) and (c), where the latter is clearer than the former. However, as Figure 5 shown in (b), the imaging quality of the PSM-based method is poor, and the two defects in the lower left corner are difficult to distinguish from the background.

[0114] For the detection of hole defects in aluminum parts, in Figure 5 in (d), serious artifacts appear on the upper side of the imaging area, which seriously interferes with the identification of the three upper-front defects of the TFM-based method. For the PSM-based method, the situation becomes worse, and only the middle defect can be identified in Figure 5 in (e), and other defects are submerged in the background. FD-RTM well suppresses the artifacts, and all defects are in sharp contrast with the Figure 5 background in (f).

[0115] The wave energy reflected by crack defects is much greater than that of void defects. Therefore, all three methods can image V-shaped cracks relatively clearly with low background interference. The contrast of the crack image of the TFM-based method is higher than that of the PSM-based method, but both methods image the sharp corners of real cracks as rounded corners. Although there are artifacts in the imaging result of FD-RTM, the imaging target is clearly recognizable, and the boundary of FD-RTM is the clearest and the shape is closest to the actual crack shape.

[0116] In addition, for further comparison, the amplitude values at the parts marked as dotted lines in Figure 5 in (a)–(i) are taken out and normalized for comparison, as Figure 6 shown in (a)–(c). In the first two cases (hole defects in copper parts and hole defects in aluminum parts), the main lobe widths of the three imaging methods are almost the same, while the side lobe of FD-RTM is the smallest in the second case (hole defects in aluminum parts). In addition, the relative difference between the two peaks of FD-RTM is the smallest, which indicates that FD-RTM still maintains a high contrast at the edge of the imaging area. For the third case (crack defects in aluminum parts), as Figure 6 shown in (c), the main lobe of FD-RTM drops fastest on both sides, which ensures the Figure 5 clear boundary in (i).

[0117] To quantitatively analyze the imaging results, a dimensionless parameter API (array performance indicator) is introduced, as shown in the following formula

[0118] API = A -6dB / λ 2 15)

[0119] where A-6dB is the area where the ratio of the peak value of the defective area image exceeds -6 dB, and λ is the wave number corresponding to the center frequency. Obviously, the smaller the API, the larger the resolution. Table 1 lists the average API of the imaging results of the void defects in the first two experiments.

[0120] Table 1 Average API of void defects in the imaging results obtained by different methods in the experiment

[0121]

[0122] Since Figure 5 the two defects in the lower left corner of (e) could not be distinguished from the background at all, so the average API of the other two defects was only calculated. It can be seen that the API of the imaging results of FD-RTM is the smallest in both experiments, and the resolution of the FD-RTM method is higher than that of the TFM-based method. The imaging times of FD-RTM in the three experiments (void defects in copper parts, void defects in aluminum parts, and crack defects in aluminum parts) were 978 s, 1019 s, and 1031 s respectively, while the imaging times of the TFM-based method were 1566 s, 1536 s, and 1554 s respectively. Therefore, the experiment also shows that the FD-RTM method is very suitable for high-quality defect imaging of curved parts.

Claims

1. A high-quality imaging method for full matrix data based on frequency-domain reverse time migration, characterized in that, It includes the following steps: (1) Input the excitation signal, sound speed distribution, and full matrix data; (2) Preprocess the excitation signal and the full matrix data respectively to obtain the excitation sound source matrix and the received sound source matrix; The preprocessing of the excitation signal is as follows: Perform Fourier transform processing on the excitation signal in the time domain to obtain the excitation sound source matrix; The preprocessing of the full matrix data is as follows: First, use a window function to act on the full matrix data to filter out the direct wave. Perform Fourier transform on the full matrix data with the direct wave filtered out in the time dimension, and then take the complex conjugate of the data in the obtained frequency domain to achieve the effect of time reversal, and organize to obtain the received sound source matrix; (3) Discretize the sound speed distribution and calculate the impedance matrix; In step (3), the calculation formula of the impedance matrix A(ω) is as follows: Among them, represents the Laplace operator obtained after discretization processing; ω represents the frequency; c represents the speed of sound; (4) According to the obtained excitation sound source matrix and received sound source matrix, use LU decomposition of the impedance matrix to solve and obtain the transmitted wave field and the received wave field respectively; The specific steps of step (4) are as follows: 4.1 Take any un-traversed frequency within the selected frequency band as the current frequency. Use the excitation sound source matrix, received sound source matrix, and impedance matrix corresponding to the current frequency as inputs, perform LU decomposition on the impedance matrix, and solve to obtain the excitation sound pressure field and the received sound pressure field corresponding to the current frequency; 4.2 Traverse all frequencies within the selected frequency band to obtain the transmitted wave field and the received wave field; (5) Process the transmitted wave field and the received wave field using the imaging condition to obtain the initial image, and perform post-processing on the initial image to obtain the imaging result; The post-processing of the initial image is as follows: First, perform sharpening processing on the initial image using Laplace filtering, then perform Hilbert transform on the sharpened initial image to obtain the envelope of the image, and finally use median filtering to filter out salt-and-pepper noise to obtain the imaging result.

2. The high-quality imaging method for full matrix data based on frequency-domain reverse time migration according to claim 1, wherein When discretizing the sound speed distribution, optimize the impedance matrix by adding a perfectly matched layer around the discretization model.

3. The high-quality imaging method for full matrix data based on frequency-domain reverse time migration according to claim 1, characterized in that In step 4.1, the calculation formula for solving the sound pressure field generated by the sound source matrix is as follows: L(ω)U(ω)[P1(ω),P2(ω),…,P n (ω)]=[S1(ω),S2(ω),…,S n (ω)] where n represents the number of excitation times; L(ω) and U(ω) represent two factors obtained after performing LU decomposition on the impedance matrix A(ω); [S1(ω), S2(ω), …, S n (ω)] represents the excitation sound source matrix or the receiving sound source matrix; [P1(ω), P2(ω), …, P n (ω)] represents the excitation sound pressure field or the receiving sound pressure field.

4. The high-quality imaging method for full matrix data based on frequency-domain reverse time migration according to claim 1, wherein In step (5), the expression of the imaging condition is: Among them, I(x, z) represents the initial image; n represents the number of excitation times; N ω represents the total number of discrete frequencies in the selected frequency band; S i (x, z, ω) represents the transmitted wavefield; R i (x, z, ω) represents the received wavefield; * represents taking the complex conjugate.

5. The high-quality imaging method for full matrix data based on frequency-domain reverse time migration according to claim 1, wherein Use the 9-point finite difference method to discretize the sound speed distribution.

Images