A Phase-Shift Ultrasonic Imaging Method with Wavenumber Domain Absorbing Boundary Conditions

Abstract

The present invention provides a phase-shift ultrasonic imaging method with frequency-wavenumber domain absorbing boundary conditions. By using a designed window selection operator, the wave field is divided into left and right boundary wave fields and an intermediate region wave field in the frequency space domain. Through Fourier transform and wave field extrapolation of the three wave fields simultaneously, and superimposing the three extrapolated wave fields, the wave field of the current layer in the frequency-wavenumber domain is obtained. Filtering is performed according to the wavenumber of the incident wave in the boundary region to achieve the effect of absorbing the pseudo-reflected waves caused by the numerical boundary during the wave field extrapolation process, thereby improving the accuracy of wave field reconstruction. The wave field of the current layer in the frequency-wavenumber domain is transformed into the frequency space domain, and a compensation operator is applied, and finally an imaging condition is applied to obtain the imaging result of the current layer. All discrete layers are traversed to integrate and obtain the imaging result of the workpiece. The present invention can simultaneously adapt to line scan data and phased array data collected in homogeneous or inhomogeneous media, and improve the imaging quality by improving the accuracy of wave field reconstruction.

Description

Technical Field

[0001] The present invention belongs to the technical field of detection methods, and particularly relates to a phase-shift ultrasonic imaging method with wavenumber-domain absorbing boundary conditions. Background Art

[0002] Ultrasonic phased array imaging is prominent in non-destructive testing (NDT) due to its portability, reliability, and relatively high precision. It is very effective in detecting invisible defects or damages in components. The two commonly used data modes in phased arrays are line scan (B-scan) and full matrix (FMC) data. Line scan data has low requirements for acquisition equipment and is easy to obtain and process, but the information in the data is limited, resulting in a low signal-to-noise ratio (SNR) of the imaging results. The emergence of FMC data has improved this phenomenon. It captures the time-domain signals of each pair of excitation and receiving elements, containing more information about the measurement object. With the continuous expansion of application fields, ultrasonic phased array imaging of parts in multi-layer structures has become a research hotspot. In these cases, when the sound speed of the medium changes in the horizontal and vertical directions, greater challenges are posed to the imaging method.

[0003] Phase shift (PSM) is an important method in ultrasonic imaging. It regards the measurement data as a wave field on the surface and reconstructs the wave field of the measurement area through extrapolation in the wavenumber-frequency domain. Phase migration was first used to process line scan data for imaging defects in layered media. Martin improved the computational efficiency of the phase shift method by introducing Stolt interpolation. Further, Wu et al. extended the phase shift to process full matrix data for imaging hole defects in layered media. Ji et al. proposed an efficient phase shift method for ultrasonic full matrix imaging of layered media. In the above applications, the sound speed of the medium only changes in the vertical direction. For the lateral change of the sound speed, Jin et al. introduced a step-by-step migration method to process B-scan data, which uses a phase shift factor to compensate for the lateral change. Chang et al. introduced a non-stationary phase shift method to process full matrix data. This scheme divides the wave field according to the sound speed at each depth and performs wave field extrapolation using the corresponding sound speed. Therefore, the phase shift-based method can process B-scan and FMC data and can also adapt to media with vertical and lateral sound speed changes.

[0004] Reconstructing the wave field in the measurement area and then implementing the imaging condition is the basic idea for obtaining high-quality defect images, and the first step is the core of various methods. The full-wave equation migration (reverse-time migration) can adapt to the complex sound speed distribution in the medium and has high accuracy, but it requires high memory and computing time. The one-way wave equation migration simplifies the full-wave equation by adopting an approximate dispersion equation, and its efficiency is improved, but each downward extrapolation step requires a finite difference operator, which still makes the time cost very high. To further improve the imaging efficiency, phase migration uses a phase shift operator to achieve the downward extrapolation of the wave field, which is actually the solution of the one-way wave equation assuming a laterally uniform sound speed. When performing wave equation migration, pseudo-reflections will appear at the edges of the calculation area. Therefore, the absorbing boundary condition is crucial for these methods to reconstruct the wave field with high accuracy and ensure the imaging quality. Phase shift migration originates from wave equation migration, so pseudo-reflections will also appear at the boundary. However, there is currently a lack of research on the absorbing boundary condition for wave field extrapolation in phase migration. Summary of the Invention

[0005] To solve the problems existing in the prior art, the present invention provides a phase shift ultrasonic imaging method with absorbing boundary conditions in the wavenumber domain. This imaging method uses a window selection operator to divide the wave field in the upper-layer frequency space domain into three parts: the left boundary wave field, the middle region wave field, and the right boundary wave field, and uses filters to absorb the incident waves for the wave fields at both boundaries during the extrapolation process, reducing the pseudo-reflected waves in the imaging result and improving the imaging quality.

[0006] A phase shift ultrasonic imaging method with absorbing boundary conditions in the wavenumber domain includes the following steps:

[0007] (1) Discretize the imaging area of the workpiece to form multiple discrete layers in the depth direction;

[0008] (2) Convert the measurement data of the workpiece into the frequency space domain to obtain the surface wave field;

[0009] (3) Using the surface wave field as the wave field of the first layer in the frequency space domain, select any unvisited layer as the current layer in the order from top to bottom, and apply a window selection operator to divide the wave field of the upper layer in the frequency space domain into a left boundary wave field, a middle region wave field, and a right boundary wave field;

[0010] Perform Fourier transform on the middle region wave field to obtain its wave field in the frequency-wavenumber domain; after converting the left boundary wave field and the right boundary wave field into the frequency-wavenumber domain respectively and applying the corresponding filters, obtain the wave fields of the left boundary and the right boundary in the frequency-wavenumber domain respectively;

[0011] In the frequency-wavenumber domain, the wave fields at the left boundary, in the middle region, and at the right boundary are extrapolated separately and then superimposed to obtain the wave field of the current layer in the frequency-wavenumber domain; the wave field of the current layer in the frequency-wavenumber domain is transformed into the frequency-space domain and the compensation operator is applied to obtain the wave field of the current layer in the frequency-space domain;

[0012] (4) Image based on the obtained wave field of the current layer in the frequency-space domain to obtain the imaging result of the current layer;

[0013] (5) Repeat steps (3) and (4) until all discrete layers are traversed, obtain the imaging results of all layers, and integrate them to obtain the imaging result of the workpiece.

[0014] The measurement data in step (2) above can be line scan data (B-scan data) or full matrix data (FMC data). For different data types, corresponding specific refined operation processing is performed. The wave field extrapolation with boundary absorption conditions in step (3) can be applied to both homogeneous and inhomogeneous media.

[0015] Preferably, in step (1), for the discretization of the imaging region, the step size in the horizontal x direction is set as Δx, which is the spacing between adjacent vibration elements; the step size in the depth z direction is Δz, where the size of Δz is determined according to the resolution requirement in the vertical direction; the number of discrete layers is N.

[0016] Preferably, in step (2), the measurement data of the workpiece is expressed as P(x, z = 0, t), and P(x, z = 0, t) is transformed into the frequency-wavenumber domain through the following Fourier transform to obtain the surface wave field P(x, z = 0, ω):

[0017]

[0018] where, i represents the imaginary unit; ω represents the frequency; x represents the position of the vibration element; t represents the time.

[0019] Preferably, in step (3), the calculation formulas for using the window selection operator to divide the wave field of the nth layer in the frequency-space domain into the left boundary wave field, the middle region wave field, and the right boundary wave field are:

[0020] P L (x, z n , ω) = P(x, z n , ω)W L (x)

[0021] P M (x, z n , ω) = P(x, z n , ω)W M (x)

[0022] PR (x, z n , ω) = P(x, z n , ω)W R (x)

[0023] wherein, P L (x, z n , ω) represents the left boundary wave field; P M (x, z n , ω) represents the middle region wave field; P R (x, z n , ω) represents the right boundary wave field; P(x, z n , ω) represents the wave field of the nth layer in the frequency - space domain, n ∈ [1, N - 1], and N represents the number of discrete layers;

[0024] W L (x) is the left - boundary window - selection operator, and its expression is:

[0025]

[0026] W M (x) is the middle - region window - selection operator, and its expression is:

[0027]

[0028] W R (x) is the right - boundary window - selection operator, and its expression is:

[0029]

[0030] In the formula, x L and x R are the left and right boundary points respectively; Δx is the discrete step length in the x - direction; w() represents the window function.

[0031] The above - mentioned three window - selection operators are window - selection functions designed for window - selection at different positions of the imaging region. The distributions of the three window - selection operators are as Figure 5 shown. At the boundaries of window - selection, a half - window function is used for transition to reduce the frequency leakage generated during Fourier transform and improve the imaging accuracy.

[0032] As a further preference, the window function is one of Hanning window, Blackman window, Flat top window, and rectangular window. More preferably, it is the Blackman window.

[0033] As a preference, in step (3), in the frequency - wavenumber domain, the following formula is used to extrapolate the wave fields of the left boundary, the middle region, and the right boundary respectively:

[0034]

[0035] Among them, P′(k x , z n+1 , ω) represents the wave field in the frequency-wavenumber domain of the (n + 1)-th layer; P(k x , z n , ω) represents the wave field transformed into the frequency-wavenumber domain after acting on the weighting operator; it is defined that is the phase shift operator; n ∈ [1, N - 1], N represents the number of discrete layers; z n represents the depth of the n-th layer; i represents the imaginary unit; ω represents the frequency; Δz represents the discretization step in the depth direction; k x represents the wavenumber in the horizontal direction; represents the total vertical wavenumber at the depth of z n , and the expression is:

[0036]

[0037] When the n-th layer is a homogeneous medium, c n represents the sound speed of the n-th layer;

[0038] When the n-th layer is a non-homogeneous medium, c n = 1 / s0(z n ), s0(z n ) represents the average value of all slownesses on the horizontal plane at the depth of z n .

[0039] The specific derivation process is as follows:

[0040] The propagation equation of ultrasonic waves in a two-dimensional medium is:

[0041]

[0042] Among them, p is the sound pressure, c is the sound speed of the medium, x and z are the two-dimensional space coordinates respectively, and t is the time. After performing the Fourier transform on t, the wave equation in the frequency domain is obtained as:

[0043]

[0044] Among them, ω is the frequency. Since the energy of multiple echoes is small and can be ignored, the wave equation can be decomposed into a unidirectional equation as follows:

[0045]

[0046] Here, the phased array probe collects signals at the surface z = 0, and the measurement targets are evenly distributed in the half-space where z > 0, as Figure 4As shown. Therefore, the wave field on the surface can be collected by the phased array probe and can be expressed as P(x, z = 0, t) or P(x, z = 0, ω), corresponding to the time domain or the frequency domain respectively. The wave field of the entire measurement area can be extrapolated from the surface wave field according to Equation (3).

[0047] For a homogeneous medium, the sound speed is independent of the spatial coordinates x and z, so it can be represented by c. At this time, the Fourier transform in the x direction can be directly performed to separate the harmonic components of each plane wave. The wave field extrapolation can be completed by applying a phase shift to each harmonic component, as shown in Equation (4):

[0048]

[0049] where k z is the wave number in the vertical direction, and the expression is as follows:

[0050]

[0051] P(k x , z = 0, ω) is obtained by Fourier transforming the measurement data P(x, z = 0, ω):

[0052]

[0053] Then, in the frequency-wavenumber domain, the calculation formula for extrapolating the wave field of the next layer from the wave field of the current layer in a layer-by-layer recursive manner is:

[0054]

[0055]

[0056] In Equation (8), c n represents the sound speed of the nth layer.

[0057] For a non-homogeneous medium, if there is no horizontal variation, it becomes a multi-layer medium case. At this time, the sound speed is only related to the depth and can be written as c(z). At this time, the extrapolation can be completed by a phase shift that varies with the depth z.

[0058] However, if there is a horizontal variation in the sound speed, the sound speed is a function of x and z, which makes it infeasible to perform only the Fourier transform on the wave field. At this time, the slowness s(x, z) is introduced. It is the reciprocal of the sound speed and can be decomposed into two components, as shown in the following formula:

[0059]

[0060] Among them, s0(z) is the reference slowness, which is the average value of all slownesses on the water surface at depth z; Δs(x, z) contains all slowness variations on the horizontal plane at depth z. The square root term in formula (3) can be expanded by Taylor series and approximated as:

[0061]

[0062] If s0(z) >> Δs(x, z), then in equation (8), c n = 1 / s0(z n ), can be written as:

[0063]

[0064] As an absorbing boundary condition, the wave incident into the boundary region should be dissipated. For the left boundary, the direction of the incident wave is opposite to the positive direction of the x-axis, so the wave number k x of the incident wave is negative. On the contrary, at the right boundary, the direction of the incident wave is the same as the x-axis direction, so the wave number k x of the incident wave is positive. Therefore, the filter at the left boundary sets the wave field values with negative wave number k x in the frequency-wave number domain to zero; the filter at the right boundary sets the wave field values with positive wave number k x in the frequency-wave number domain to zero.

[0065] Therefore, in the frequency-wave number domain, the wave field extrapolation formulas for the wave fields at the left boundary, in the middle region, and at the right boundary are as follows:

[0066]

[0067]

[0068]

[0069] Among them, P L (k x , z n , ω), P R (k x , z n , ω) are the wave fields at the left boundary and the right boundary after the wave fields at the left boundary and the right boundary of the nth layer are transformed into the frequency-wave number domain and the filter is applied; P M (k x , z n , ω) is the wave field in the middle region of the nth layer in the frequency-wave number domain; P′ L (k x , z n+1 , ω), P′ L (k x , z n+1,ω), P' R (k x ,z n+1 ,ω) represent the wave fields at the left boundary, the middle region, and the right boundary of the (n + 1)-th layer in the frequency-wavenumber domain respectively.

[0070] Preferably, in step (3), in the frequency-wavenumber domain, the wave fields at the left boundary, the middle region, and the right boundary of the n-th layer are extrapolated and then superimposed respectively to obtain the wave field of the (n + 1)-th layer in the frequency-wavenumber domain;

[0071] Among them, the superimposition formula is as follows:

[0072] P'(k x ,z n+1 ,ω) = P' L (k x ,z n+1 ,ω) + P' M (k x ,z n+1 ,ω) + P' R (k x ,z n+1 ,ω)

[0073] In the formula, P'(k x ,z n+1 ,ω) represents the wave field of the (n + 1)-th layer in the frequency-wavenumber domain;

[0074] P' L (k x ,z n+1 ,ω), P' L (k x ,z n+1 ,ω), P' R (k x ,z n+1 ,ω) are the wave fields obtained by extrapolating the wave fields at the left boundary, the middle region, and the right boundary of the n-th layer respectively.

[0075] Preferably, in step (3), the wave field of the (n + 1)-th layer in the frequency-wavenumber domain is transformed into the frequency-space domain and a compensation operator is applied to obtain its wave field in the frequency-space domain;

[0076] The calculation formula for applying the compensation operator is:

[0077]

[0078] Among them, P(x, z n+1 ,ω) represents the wave field of the (n + 1)-th layer in the frequency-space domain; P'(x, z n+1, ω) represents the wave field obtained by converting the wave field of the (n + 1)-th layer in the frequency-wavenumber domain to the frequency space domain; represents the compensation operator; i represents the imaginary unit; ω represents the frequency; x represents the position of the vibration element; Δz represents the discretization step in the depth direction; Δs(x, z n ) represents all the slowness changes on the horizontal plane at a depth of z n (i.e., the depth of the (n - 1)-th layer); n ∈ [1, N - 1], and N represents the number of discrete layers.

[0079] In this technical solution, if the sound speed of the n-th layer does not change in the horizontal direction (homogeneous medium), that is, Δs(x, z n ) = 0, then the compensation in the horizontal direction can be omitted; therefore, before traversing the discrete layers, it can be judged first whether the n-th layer is a homogeneous medium. If so, the wave field P(x, z n+1 , ω) of the (n + 1)-th layer in the frequency space domain can be directly obtained by converting its wave field in the frequency-wavenumber domain.

[0080] If the sound speed of the n-th layer changes in the horizontal direction (non-homogeneous medium), then the wave field P(x, z n+1 , ω) of the (n + 1)-th layer in the frequency space domain is obtained after performing the sound speed compensation in the horizontal direction (applying the compensation operator) after converting its wave field in the frequency-wavenumber domain to the frequency-wavenumber domain.

[0081] Preferably, when the measurement data is line scan data, in step (4), the imaging condition for imaging according to the obtained wave field P(x, z n , ω) of the n-th layer is:

[0082] I(x, z n ) = P(x, z n , t = 0) = ∫dωP(x, z n , ω)

[0083] where, I(x, z n ) represents the imaging result of the n-th layer; ω represents the frequency; x represents the position of the vibration element; z n represents the depth of the n-th layer, n ∈ [1, N], and N represents the number of discrete layers.

[0084] Line scan (B-scan) data is a basic type of data collected by an ultrasonic phased array probe, which is convenient to process and has low requirements for the acquisition equipment. Wave field extrapolation can be directly used for wave field reconstruction of B-scan data. The explosive reflector model assumes that the sound source is located at the defect and is excited at t = 0, so the defect image is the wave field at t = 0. It should be noted that according to the requirements of the explosive reflector model, the sound speed should be divided by 2 in all cases.

[0085] Therefore, after reconstructing the wave field in the f-x domain (frequency space domain), the imaging result can be obtained by implementing the following imaging condition:

[0086] I(x, z) = P(x, z, t = 0) = ∫dωP(x, z, ω)

[0087] where I(x, z) is the imaging result.

[0088] Introducing the wave field extrapolation method with absorbing boundary conditions (step (3)) into the imaging method of B-scan data can reduce the pseudo-reflected waves existing in the reconstructed wave field to improve the imaging quality. For the imaging method of the present invention for line scan data, the measurement area is discretized into N layers in the depth direction, and the defect images are calculated layer by layer. First, initialize the parameters, take the B-scan data P(t, 0, x) as the input, then transform P(t, 0, x) to the f-x domain, and use the extrapolation method with absorbing boundary conditions in step (3) to realize the wave field extrapolation from P(x, z n , ω) to P(x, z n+1 , ω). In addition, according to the imaging condition, the reconstructed wave field is superimposed in the ω dimension to obtain the imaging result of the nth layer. The above steps are iterated at each depth, and finally the defect image I(x, z) of the measurement area is obtained.

[0089] Preferably, when the measurement data is full matrix data (FMC data), the imaging process is as follows:

[0090] Select any un-traversed element as the current element, convert the sound source wave field and received wave field when the current element is excited to the frequency space domain to obtain the surface sound source wave field and surface received wave field when the current element is excited, and substitute the obtained surface sound source wave field and surface received wave field into step (3) respectively to replace the surface wave field in step (2);

[0091] After being processed by step (3), the sound source wave field and received wave field of the current layer in the frequency space domain are obtained respectively;

[0092] In step (4), apply the imaging condition to the sound source wave field and received wave field of the current layer in the frequency space domain to obtain the imaging result of the current layer, and enter step (5);

[0093] After being processed by step (5), the imaging results of all discrete layers corresponding to the current element are obtained, and after integration, the imaging result generated when the current element is excited is obtained;

[0094] Repeat steps (2) to (5) until all elements are traversed, obtain the imaging results of all elements, and perform superposition to obtain the imaging result of the workpiece.

[0095] Among them, the current vibration element is denoted as the m-th vibration element, where m ∈ [1, M], and M represents the number of vibration elements; the sound source wave field of the current layer (the n-th layer, n ∈ [1, N], and N represents the number of discrete layers) in the frequency-space domain is denoted as S m (x, z n , ω), and the received wave field is denoted as P m (x, z n , ω). When n = 1, the sound source wave field and the received wave field of the first layer are both obtained by converting the full matrix data.

[0096] In the frequency-wavenumber domain, the extrapolation calculation formula for the received wave field is:

[0097]

[0098] Among them, P′ m (k x , z n+1 , ω) represents the received wave field of the (n + 1)-th layer in the frequency-wavenumber domain; P m (k x , z n , ω) represents the received wave field converted to the frequency-wavenumber domain after applying the weighting operator; n ∈ [1, N - 1], and N represents the number of discrete layers.

[0099] The calculation formula for acoustic velocity compensation of the received wave field in the horizontal direction is:

[0100]

[0101] Among them, P m (x, z n+1 , ω) represents the received wave field of the (n + 1)-th layer in the frequency-space domain; P′ m (x, z n+1 , ω) represents the received wave field after converting the received wave field of the (n + 1)-th layer in the frequency-wavenumber domain to the frequency-space domain; n ∈ [1, N - 1], and N represents the number of discrete layers.

[0102] Since the phase shift operator for wave field extrapolation of the sound source wave field and the compensation operator in the horizontal direction are complex conjugate relations with the corresponding phase shift operator and compensation operator of the received wave field; therefore, in the frequency-wavenumber domain, the extrapolation calculation formula for the sound source wave field is:

[0103]

[0104] Among them, S′ m (k x , z n+1 , ω) represents the sound source wave field of the (n + 1)-th layer in the frequency-wavenumber domain; S m (k x , z n, ω) represents the sound source wave field converted to the frequency-wavenumber domain after the action of the weighting operator; n ∈ [1, N - 1], and N represents the number of discrete layers.

[0105] The calculation formula for the action compensation operator of the sound source wave field is:

[0106]

[0107] Among them, S m (x, z n+1 , ω) represents the sound source wave field in the frequency spatial domain of the nth layer; S′ m (x, z n+1 , ω) represents the sound source wave field converted from the sound source wave field in the frequency-wavenumber domain of the (n + 1)th layer to the frequency spatial domain; n ∈ [1, N - 1], and N represents the number of discrete layers.

[0108] As a further preference, the imaging condition for imaging from the sound source wave field and the received wave field in the frequency spatial domain of the nth layer is:

[0109]

[0110] Among them, I m (x, z n ) represents the imaging result of the nth layer when the mth vibration element is excited; P m (x, z n , ω) represents the received wave field in the frequency spatial domain of the nth layer; represents the conjugate of S m (x, z n , ω), S m (x, z n , ω) represents the sound source wave field in the frequency spatial domain of the nth layer; ω represents the frequency; x represents the position of the vibration element; z n represents the depth of the nth layer, n ∈ [1, N], and N represents the number of discrete layers.

[0111] Full matrix capture (FMC) is a widely used capture mode of phased array probes, in which each vibration element in the phased array emits one by one, and all vibration elements receive signals. Therefore, FMC data contains the time-domain signals of each excitation vibration element-receiving vibration element pair in the phased array. Compared with B-scan data, FMC data contains more information, greatly improving the imaging resolution and signal-to-noise ratio (SNR). Here, the extrapolation of the wave field with absorbing boundary conditions (step (3)) is extended to FMC data for high-quality imaging. In this case, the wave is excited at one vibration element in the phased array and reaches the defect at time t = τ m (x, z) and then the scattered wave is received by the phased array. Therefore, the assumption of the explosion reflection model should be modified because the scattering point is located at the defect and at t = τm Excitation at (x, z), so the defect image is at t = τ m The wave field at (x, z). Here, the sound source (incident) wave field from the excitation element is extrapolated to simulate wave incidence and compensate for the propagation time τ m (x, z). Therefore, its imaging condition is as follows:

[0112]

[0113] where I m (x, z) is the imaging result obtained when the m-th element is excited, and P m (x, z, ω) is the received wave field extrapolated from the phased array received signal, denotes the conjugate of S m (x, z, ω), and S m (x, z, ω) is the sound source wave field extrapolated from the excitation element. Among them, the excitation wave field of the m-th element can be obtained by the following formula:

[0114] S m (x = x m , z = 0, t = 0) = δ(x = x m , t = 0) (18)

[0115] where δ is the impulse function; x m is the position of the m-th element.

[0116] As a further preference, the calculation formula for superimposing the imaging results of all elements is:

[0117]

[0118] where I(x, z) represents the imaging result of the workpiece; I m (x, z) represents the imaging result of the m-th element.

[0119] The imaging method for FMC data of the present invention first initializes parameters, and records the FMC data and the excitation signal as P m (t, 0, ω) and S m (t, 0, ω) and inputs them. This algorithm has two layers of loops. The outer loop traverses all excitation elements; the inner loop traverses all discrete depths (all discrete layers) in the imaging area. Apply the wave field extrapolation method with boundary absorption conditions in step (3) to the received wave field and the sound source wave field simultaneously in the f - k domain from z n-1 extrapolate to z n during the process, and transform to the f - x domain and apply the compensation operator and then apply the imaging condition to obtain the imaging result I m (x, z n)。After the two-layer loop ends, the imaging results generated when all vibration elements are excited are superimposed on the imaging result I(x, z). It should be noted that during the calculation process, the sound speed does not need to be divided by 2 because in this case, both the wave transmission and reflection processes are considered.

[0120] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0121] The phase-shift ultrasonic imaging method with frequency-wave number domain absorbing boundary conditions of the present invention uses a designed window selection operator to divide the wave field into two boundary wave fields on the left and right and an intermediate region wave field in the frequency space domain. At the same time, Fourier transforms and wave field extrapolations are performed on the three wave fields, and the three extrapolated wave fields are superimposed to obtain the wave field of the current layer in the frequency-wave number domain. Filtering is performed according to the wave number of the incident wave in the boundary region to achieve the effect of pseudo-reflected waves caused by numerical boundaries during the absorbing wave field extrapolation process, thereby improving the accuracy of wave field reconstruction; the wave field of the current layer in the frequency-wave number domain is converted to the frequency space domain, and a compensation operator is applied, and finally an imaging condition is applied to obtain the imaging result of the current layer. All discrete layers are traversed and integrated to obtain the imaging result of the workpiece. The present invention can simultaneously adapt to line scan data (B-scan data) and phased array data (FMC data) collected in homogeneous or inhomogeneous media, and improve the imaging quality by improving the accuracy of wave field reconstruction. Description of the Drawings

[0122] Figure 1 It is a flowchart of wave field extrapolation with absorbing boundaries for inhomogeneous media in Embodiments 1 and 2;

[0123] Figure 2 It is a flowchart of imaging for line scan data in Embodiment 1;

[0124] Figure 3 It is a flowchart of imaging for full matrix data in Embodiment 2;

[0125] Figure 4 It is a schematic diagram of phased array detection;

[0126] Figure 5 It is a distribution diagram of three window selection functions (window selection operator);

[0127] Figure 6 In (a), it is the imaging result obtained by reconstructing the wave field using the one-step extrapolation method for homogeneous media; (b) is the imaging result obtained by reconstructing the wave field using the two-step extrapolation method for inhomogeneous media;

[0128] Figure 7 In (a) and (b), they are respectively the imaging results obtained by forward modeling of homogeneous and inhomogeneous media in Figure 11 and Figure 1 using the method of Figure 6 ;

[0129] Figure 8 In (a), it is the detection diagram of the hole defect of the copper part; in (b), it is the detection diagram of the hole defect of the fiber-reinforced PEEK part.

[0130] Figure 9 In (a), it is the imaging result of the traditional phase-shift imaging algorithm (without absorption boundary) for line-scan data; in (b), it is the imaging result of the hole defect of the aluminum part by using the imaging method in Embodiment 1; in (c), it is the imaging result of the traditional phase-shift imaging algorithm (without absorption boundary) for full matrix data; in (d), it is the imaging result of the hole defect of the aluminum part by using the imaging method in Embodiment 2.

[0131] Figure 10 In (a), it is the imaging result of the traditional phase-shift imaging algorithm (without absorption boundary) for line-scan data; in (b), it is the imaging result of the hole defect of the fiber-reinforced PEEK part by using the imaging method in Embodiment 1; in (c), it is the imaging result of the traditional phase-shift imaging algorithm (without absorption boundary) for full matrix data; in (d), it is the imaging result of the hole defect of the fiber-reinforced PEEK part by using the imaging method in Embodiment 2.

[0132] Figure 11 It is the wave field extrapolation flow chart with absorption boundary for homogeneous media in Embodiments 1 and 2. Detailed implementation manners

[0133] 1. Imaging method for line-scan data (B-scan data)

[0134] As Figure 2 shown, a phase-shift ultrasonic imaging method with absorption boundary conditions in the frequency-space domain includes the following steps:

[0135] (1) Discretize the imaging area of the workpiece to form multiple discrete layers in the depth direction.

[0136] Among them, the step size in the horizontal x direction is Δx, that is, the spacing between adjacent vibration elements; the step size in the depth z direction is Δz, where the size of Δz is determined according to the resolution requirement in the vertical direction; the number of discrete layers is N.

[0137] (2) Convert the line-scan data (P(x, z = 0, t)) of the workpiece to the frequency-space domain to obtain the surface wave field P(x, z = 0, ω).

[0138]

[0139] Among them, ω represents frequency; x represents the position of the vibration element; t represents time.

[0140] (3) Use the surface wave field as the wave field of the first layer in the frequency-space domain, i.e., (P(x, z1, ω) = P(x, z = 0, ω)), and traverse the discrete layers in the order from top to bottom; apply a window selection operator to divide the wave field of the nth layer in the frequency-space domain into a left boundary wave field, a middle region wave field, and a right boundary wave field;

[0141] Perform a Fourier transform on the middle region wave field to obtain its wave field in the frequency-wavenumber domain; after converting the left boundary wave field and the right boundary wave field to the frequency-wavenumber domain respectively, apply the corresponding filters to obtain the wave fields of the left boundary and the right boundary in the frequency-wavenumber domain respectively;

[0142] In the frequency-wavenumber domain, extrapolate and superimpose the wave fields of the left boundary, the middle region, and the right boundary respectively to obtain the wave field of the (n + 1)th layer in the frequency-wavenumber domain; convert the wave field of the (n + 1)th layer in the frequency-wavenumber domain to the frequency-space domain and apply a compensation operator to obtain the wave field of the (n + 1)th layer in the frequency-space domain.

[0143] Since the data of the first layer can be directly obtained from the measurement data, the extrapolation traversal can start from the second layer.

[0144] The specific process is as follows:

[0145] Apply the window selection operator:

[0146] Apply the window selection operator to divide the wave field P(x, z n , ω) of the nth layer in the frequency-space domain into a left boundary wave field P L (x, z n , ω), a middle region wave field P M (x, z n , ω), and a right boundary wave field P R (x, z n , ω); among them, the calculation formula is as follows:

[0147] P L (x, z n , ω) = P(x, z n , ω)W L (x) (2)

[0148] P M (x, z n , ω) = P(x, z n , ω)W M (x) (3)

[0149] P R (x, z n , ω) = P(x, z n , ω)W R (x) (4)

[0150] In the above formula, n ∈ [1, N - 1];

[0151] W L (x) is the left - boundary window - selection operator, and its expression is:

[0152]

[0153] W M (x) is the middle - region window - selection operator, and its expression is:

[0154]

[0155] W R (x) is the right - boundary window - selection operator, and its expression is:

[0156]

[0157] In the formula, x L and x R are the left and right boundary points respectively; Δx is the discrete step length in the x - direction; w() represents the Blackman window function.

[0158] Respectively, transform P L (x, z n , ω) and P R (x, z n , ω) to the frequency - wavenumber domain and then apply the corresponding filters to obtain P L (k x , z n , ω) and P R (k x , z n , ω); transform P M (x, z n , ω) to the frequency - wavenumber domain to obtain P M (k x , z n , ω).

[0159] As the absorbing boundary condition, the waves incident on the boundary region should be dissipated. For the left boundary, the direction of the incident wave is opposite to the positive x - axis direction, so the wavenumbers k x of the incident waves are all negative. On the contrary, at the right boundary, the direction of the incident wave is the same as the x - axis direction, so the wavenumbers k x of the incident waves are all positive. Therefore, the filter at the left boundary sets the wave - field values with negative wavenumbers k x in the frequency - wavenumber domain to zero; the filter at the right boundary sets the wave - field values with positive wavenumbers k x in the frequency - wavenumber domain to zero.

[0160] Wave - field extrapolation in the frequency - wavenumber domain:

[0161] Extrapolate P L (k x , z n , ω), P M (k x , z n , ω) and P R (k x , z n , ω) respectively according to the following formulas to obtain P′ L (k x , z n+1 , ω), P′ M (k x , z n+1 , ω) and P′ R (k x , z n+1 , ω):

[0162]

[0163]

[0164]

[0165] In formulas (8) to (10), is the phase shift operator; n ∈ [1, N - 1], N represents the number of discrete layers; z n represents the depth of the nth layer; i represents the imaginary unit; ω represents the frequency; Δz represents the discretization step in the depth direction; k x represents the wave number in the horizontal direction; represents the total vertical wave number at depth z n , and the expression is:

[0166]

[0167] When the nth layer is a homogeneous medium, c n represents the sound speed of the nth layer;

[0168] When the nth layer is a non - homogeneous medium, c n = 1 / s0(z n ), s0(z n ) represents the average value of all slownesses on the horizontal plane at depth z n .

[0169] The derivation process of the above extrapolation formula is as follows:

[0170] The propagation equation of ultrasonic waves in two - dimensional media is:

[0171]

[0172] Wherein, p is the sound pressure, c is the sound speed of the medium, x and z are two-dimensional space coordinates respectively, and t is the time. After performing the Fourier transform on t, the wave equation in the frequency domain is obtained as follows:

[0173]

[0174] Wherein, ω is the frequency. Since the energy of multiple echoes is small and can be ignored, the wave equation can be decomposed into a unidirectional equation as follows:

[0175]

[0176] Here, the phased array probe collects signals at the surface z = 0 during measurement, and the measurement targets are evenly distributed in the half-space where z>0, as Figure 4 shown. Therefore, the wave field on the surface can be directly collected by the phased array probe and can be expressed as P(x, z = 0, t) or P(x, z = 0, ω), corresponding to the time domain or the frequency domain respectively. The wave field of the entire measurement area can be extrapolated from the surface wave field according to formula (14).

[0177] For a homogeneous medium, the sound speed is independent of the spatial coordinates x and z, so it can be represented by c. At this time, the Fourier transform in the x direction can be directly performed to separate the harmonic components of each plane wave. The wave field extrapolation can be completed by applying a phase shift to each harmonic component, as shown in formula (15):

[0178]

[0179] Wherein, k z is the wave number in the vertical direction, and its expression is as follows:

[0180]

[0181] P(k x , z = 0, ω) is obtained by performing the Fourier transform on the measurement data P(x, z = 0, ω):

[0182]

[0183] Then, in the frequency-wave number domain, the calculation formula for extrapolating the wave field of the current layer from the wave field of the upper layer in a layer-by-layer recursive manner is:

[0184]

[0185]

[0186] In formula (19), c n represents the sound speed of the nth layer, and at this time n ∈ [1, N - 1].

[0187] For a non-uniform medium, if there is no horizontal variation, it becomes a case of a multi-layer medium. In this case, the sound speed is only related to the depth and can be written as c(z). At this time, the extrapolation can be completed by the phase shift varying with the depth z.

[0188] However, if there is a horizontal variation in the sound speed, the sound speed is a function of x and z, which makes it infeasible to perform only a Fourier transform on the wave field. At this time, the slowness s(x, z) is introduced. It is the reciprocal of the sound speed and can be decomposed into two components as shown in the following equation:

[0189]

[0190] where s0(z) is the reference slowness, which is the average value of all slownesses on the water surface at depth z; Δs(x, z) contains all the slowness variations on the horizontal plane at depth z. The square root term in equation (14) can be Taylor-expanded and approximated as:

[0191]

[0192] If s0(z) >> Δs(x, z), then in equation (19), c n = 1 / s0(z n ), and equation (19) can be written as:

[0193]

[0194] Superpose the obtained P′ L (k x , z n+1 , ω), P′ M (k x , z n+1 , ω) and P′ R (k x , z n+1 , ω) to obtain P′(k x , z n+1 , ω):

[0195] P′(k x , z n+1 , ω) = P′ L (k x , z n+1 , ω) + P′ M (k x , z n+1 , ω) + P′ R (k x , z n+1 , ω) (23)

[0196] Compensation of the sound speed in the horizontal direction:

[0197] Convert \(P'(k x , z n+1 , \omega)\) to the frequency spatial domain to obtain \(P'(x, z n+1 , \omega)\), and then apply the compensation operator to obtain the wave field \(P(x, z n+1 , \omega)\) of the \((n + 1)\)-th layer in the frequency spatial domain;

[0198] Among them, the calculation formula is:

[0199]

[0200] Among them, represents the compensation operator; \(i\) represents the imaginary unit; \(\omega\) represents the frequency; \(x\) represents the position of the vibration element; \(\Delta z\) represents the discretization step in the depth direction; \(\Delta s(x, z n )\) represents all the slowness changes on the horizontal plane at a depth of \(z n \); \(n\in[1, N - 1]\), and \(N\) represents the number of discrete layers.

[0201] If the sound speed of the \(n\)-th layer does not change in the horizontal direction (homogeneous medium), that is, \(\Delta s(x, z n ) = 0\), then the compensation in the horizontal direction can be omitted; therefore, before traversing the discrete layers, it is possible to first judge whether the \(n\)-th layer is a homogeneous medium. If so, the wave field \(P(x, z n+1 , \omega)\) of the \((n + 1)\)-th layer in the frequency spatial domain can be directly obtained by converting its wave field in the frequency wavenumber domain, as shown in Figure 11 .

[0202] If the sound speed of the \(n\)-th layer changes in the horizontal direction (non - homogeneous medium), then the wave field \(P(x, z n+1 , \omega)\) of the \((n + 1)\)-th layer in the frequency spatial domain is obtained by converting its wave field in the frequency wavenumber domain to the frequency wavenumber domain and then performing horizontal - direction sound - speed compensation (applying the compensation operator), as shown in Figure 1 .

[0203] (4) Adopt the following imaging condition to perform imaging according to the obtained wave field of the \(n\)-th layer in the frequency spatial domain;

[0204] I(x, z n ) = P(x, z n , t = 0) = \int d\omega P(x, z n , \omega)\ (25)

[0205] At this time, \(I(x, z n )\) is the imaging result of the \(n\)-th layer; \(n\in[1, N]\).

[0206] The line scan (B-scan) data is a basic type of data collected by an ultrasonic phased array probe. It is convenient to process and has low requirements for the acquisition equipment. Wavefield extrapolation can be directly used for wavefield reconstruction of B-scan data. The explosive reflector model assumes that the sound source is located at the defect and is excited at t = 0. Therefore, the defect image is the wavefield at t = 0. It should be noted that according to the requirements of the explosive reflector model, the sound speed should be divided by 2 in all cases.

[0207] Therefore, after reconstructing the wavefield in the f-k domain (frequency-wavenumber domain), the imaging result can be obtained by implementing the following imaging condition:

[0208] I(x, z) = P(x, z, t = 0) = ∫dωP(x, z, ω)

[0209] where I(x, z) is the imaging result.

[0210] (5) Repeat steps (3) and (4) until all discrete layers are traversed, obtain the imaging results of all layers, and integrate them to obtain the imaging result of the workpiece.

[0211] Introducing the wavefield extrapolation method with absorbing boundary conditions (step (3)) into the imaging method of B-scan data can reduce the pseudo-reflected waves existing during wavefield reconstruction, so as to improve the imaging quality.

[0212] 2. Imaging method for full matrix data (FMC data)

[0213] As Figure 3 shown, a phase-shift ultrasonic imaging method with wavenumber-domain absorbing boundary conditions includes the following steps:

[0214] (1) Discretize the imaging area of the workpiece to form multiple discrete layers in the depth direction;

[0215] where the step size in the horizontal x direction is Δx, which is the spacing between adjacent vibration elements; the step size in the depth z direction is Δz, and the size of Δz is determined according to the resolution requirements in the vertical direction; the number of discrete layers is N.

[0216] (2) Select any un-traversed vibration element as the current vibration element m, and convert the sound source wavefield S m (x, z = 0, t) and the received wavefield P m (x, z = 0, t) when the current vibration element m is excited to the frequency space domain to obtain the surface sound source wavefield S m (x, z = 0, ω) and the surface received wavefield P m (x, z = 0, ω) in the frequency space domain when the current vibration element is excited;

[0217] Taking the received wavefield as an example:

[0218]

[0219] Among them, ω represents frequency; x represents the position of the vibration element; t represents time.

[0220] The conversion process of the sound source wave field is carried out according to the conversion process of the received wave field.

[0221] (3) Using the surface sound source wave field and the surface received wave field as the sound source wave field and the received wave field in the frequency space domain of the first layer, that is, S m (x, z = 0, ω) = S m (x, z1, ω), P m (x, z = 0, ω) = P m (x, z1, ω)), traverse the discrete layers in the order from top to bottom;

[0222] Taking the received wave field as an example:

[0223] Apply the window selection operator to divide the received wave field of the nth layer in the frequency space domain into the left boundary received wave field, the middle region received wave field, and the right boundary received wave field;

[0224] Perform Fourier transform on the middle region received wave field to obtain its received wave field in the frequency wavenumber domain; after converting the left boundary received wave field and the right boundary received wave field to the frequency wavenumber domain respectively and applying the corresponding filters, obtain the received wave fields of the left boundary and the right boundary in the frequency wavenumber domain respectively;

[0225] In the frequency wavenumber domain, extrapolate and superimpose the received wave fields of the left boundary, the middle region, and the right boundary respectively to obtain the received wave field of the (n + 1)th layer in the frequency wavenumber domain; convert the received wave field of the (n + 1)th layer in the frequency wavenumber domain to the frequency space domain and apply the compensation operator to obtain the received wave field of the (n + 1)th layer in the frequency space domain.

[0226] Since the data of the first layer can be directly obtained from the measurement data, the extrapolation traversal can start from the second layer.

[0227] The specific derivation process is as follows:

[0228] Apply the window selection operator:

[0229] Apply the window selection operator to divide the received wave field P m (x, z n , ω) of the nth layer in the frequency space domain into the left boundary wave field P mL (x, z n , ω), the middle region wave field P mM (x, z n , ω) and the right boundary wave field P mR (x, z n, ω); where the calculation formula is as follows:

[0230] P mL (x, z n , ω) = P m (x, z n , ω)W L (x) (2)

[0231] P mM (x, z n , ω) = P m (x, z n , ω)W M (x) (3)

[0232] P mR (x, z n , ω) = P m (x, z n , ω)W R (x) (4)

[0233] In the above formula, n ∈ [1, N - 1];

[0234] W L (x) is the left - boundary window - selection operator, and its expression is:

[0235]

[0236] W M (x) is the middle - region window - selection operator, and its expression is:

[0237]

[0238] W R (x) is the right - boundary window - selection operator, and its expression is:

[0239]

[0240] In the formula, x L and x R are the left and right boundary points respectively; Δx is the discrete step length in the x - direction; w() represents the Blackman window function.

[0241] After transforming P mL (x, z n , ω) and P mR (x, z n , ω) to the frequency - wavenumber domain and applying the corresponding filters, P mL (k x , z n , ω) and P mR (k x , z n, ω); Transform P mM (x, z n , ω) to the frequency-wavenumber domain to obtain P mM (k x , z n , ω).

[0242] As an absorbing boundary condition, the wave incident on the boundary region should be dissipated. At the left boundary, the direction of the incident wave is opposite to the positive x-axis direction, so the wavenumber k x of the incident wave is negative. On the contrary, at the right boundary, the direction of the incident wave is the same as the x-axis direction, so the wavenumber k x of the incident wave is positive. Therefore, the filter at the left boundary sets the wave field values with negative wavenumber k x in the frequency-wavenumber domain to zero; the filter at the right boundary sets the wave field values with positive wavenumber k x in the frequency-wavenumber domain to zero.

[0243] Wave extrapolation in the frequency-wavenumber domain:

[0244] Extrapolate P mL (k x , z n , ω), P mM (k x , z n , ω) and P mR (k x , z n , ω) respectively according to the following formulas to obtain P′ mL (k x , z n+1 , ω), P′ mM (k x , z n+1 , ω) and P′ mR (k x , z n+1 , ω):

[0245]

[0246]

[0247]

[0248] In Eqs. (8) - (10), is the phase shift operator; n ∈ [1, N - 1], N represents the number of discrete layers; z n represents the depth of the nth layer; i represents the imaginary unit; ω represents the frequency; Δz represents the discretization step in the depth direction; k x represents the wavenumber in the horizontal direction; represents at depth z nThe total vertical wavenumber at [location], with the expression:

[0249]

[0250] When the nth layer is a homogeneous medium, c n represents the sound speed of the nth layer;

[0251] When the nth layer is a non - homogeneous medium, c n = 1 / s0(z n ), s0(z n ) represents the average value of all slownesses on the horizontal plane at a depth of z n .

[0252] The derivation process of the above extrapolation formula is as follows:

[0253] The propagation equation of ultrasonic waves in a two - dimensional medium is:

[0254]

[0255] where p is the sound pressure, c is the sound speed of the medium, x and z are the two - dimensional space coordinates respectively, and t is the time. After performing the Fourier transform on t, the wave equation in the frequency domain is obtained as:

[0256]

[0257] where ω is the frequency. Since the energy of multiple echoes is small and can be ignored, the wave equation can be decomposed into a unidirectional equation as follows:

[0258]

[0259] Here, the phased array probe collects signals at the surface z = 0 during measurement, and the measurement targets are evenly distributed in the half - space of z>0, as Figure 4 shown. Therefore, the wave field on the surface can be directly collected by the phased array probe and can be expressed as P(x, z = 0, t) or P(x, z = 0, ω), corresponding to the time domain or the frequency domain respectively. The wave field of the entire measurement area can be extrapolated from the surface wave field according to formula (14).

[0260] For a homogeneous medium, the sound speed is independent of the spatial coordinates x and z, so it can be represented by c. At this time, the Fourier transform in the x - direction can be directly performed to separate the harmonic components of each plane wave. The wave field extrapolation can be completed by applying a phase shift to each harmonic component, as shown in formula (15):

[0261]

[0262] where k z is the vertical wavenumber, and its expression is as follows:

[0263]

[0264] P(k x , z = 0, ω) is obtained by Fourier transform of the measured data P(x, z = 0, ω):

[0265]

[0266] Then, in the frequency-wavenumber domain, the calculation formula for extrapolating the wave field of the next layer from the wave field of the current layer in a layer-by-layer recursive manner is:[[]]

[0267]

[0268]

[0269] In formula (19), c n represents the sound speed of the nth layer.[[]]

[0270] For a non-uniform medium, if there is no horizontal variation, it becomes the case of a multi-layer medium. At this time, the sound speed is only related to the depth and can be written as c(z). At this time, the extrapolation can be completed by the phase shift varying with the depth z.[[]]

[0271] However, if there is a horizontal variation in the sound speed, the sound speed is a function of x and z, which makes it infeasible to perform only a Fourier transform on the wave field. At this time, the slowness s(x, z) is introduced. It is the reciprocal of the sound speed and can be decomposed into two components as shown in the following formula:[[]]

[0272]

[0273] Among them, s0(z) is the reference slowness, which is the average value of all slownesses on the water surface at depth z; Δs(x, z) contains all the slowness variations on the horizontal plane at depth z. The square root term in formula (14) can be Taylor-expanded and approximated as:[[]]

[0274]

[0275] If s0(z) >> Δs(x, z), then in formula (19), c n = 1 / s0(z n ), formula (19) can be written as:[[]]

[0276]

[0277] Substitute the extrapolated P′ mL (k x , z n+1 , ω), P′ mM (k x , zn+1 , ω) and P' mR (k x , z n+1 , ω) are superimposed to obtain P' m (k x , z n+1 , ω):

[0278] P' m (k x , z n+1 , ω) =

[0279] P' mL (k x , z n+1 , ω) + P' mM (k x , z n+1 , ω) + P' mR (k x , z n+1 , ω) (23)

[0280] Compensation of sound speed in the horizontal direction:

[0281] Convert P' m (k x , z n+1 , ω) to the frequency spatial domain to obtain P' m (x, z n+1 , ω), and then apply the compensation operator to obtain the wave field P m (x, z n+1 , ω) of the (n + 1)-th layer in the frequency spatial domain;

[0282] Among them, the calculation formula is as follows:

[0283]

[0284] Among them, represents the compensation operator; i represents the imaginary unit; ω represents the frequency; x represents the position of the vibration element; Δz represents the discretization step in the depth direction; Δs(x, z n ) represents all the slowness changes on the horizontal plane at a depth of z n ; n ∈ [1, N - 1], and N represents the number of discrete layers.

[0285] If there is no change in the sound speed of the n-th layer in the horizontal direction (homogeneous medium), that is, Δs(x, z n ) = 0, then the compensation in the horizontal direction can be omitted; therefore, before traversing the discrete layers, it can be judged whether the n-th layer is a homogeneous medium. If so, the wave field P(x, z n+1 , ω) of the (n + 1)-th layer in the frequency spatial domain can be directly obtained by converting its wave field in the frequency wavenumber domain, asFigure 11 as shown

[0286] If the sound speed in the n-th layer varies in the horizontal direction (non-uniform medium), then the wave field P(x, z n+1 , ω) in the (n + 1)-th layer in the frequency space domain is obtained by compensating the sound speed in the horizontal direction (applying a compensation operator) after converting its wave field in the frequency-wavenumber domain to the frequency-wavenumber domain, as Figure 1 shown

[0287] Processing the source wave field of the n-th layer according to the above processing process of the received wave field, the wave field S m (x, z n+1 , ω) in the (n + 1)-th layer in the frequency space domain is obtained

[0288] Among them, since the phase shift operator and the horizontal compensation operator for the extrapolation of the wave field of the source wave field and the corresponding phase shift operator and compensation operator of the received wave field are complex conjugate relations; therefore, in the frequency-wavenumber domain, the extrapolation calculation formula of the source wave field is

[0289]

[0290]

[0291]

[0292] At this time, S mM (k x , z n , ω) represents the source wave field in the middle region of the n-th layer in the frequency-wavenumber domain; S mL (k x , z n , ω), S mR (k x , z n , ω) respectively represent the source wave fields on the left boundary and the right boundary of the source wave field of the n-th layer after converting to the frequency-wavenumber domain and applying the filter; S′ mL (k x , z n+1 , ω), S′ mM (k x , z n+1 , ω), S′ mR (k x , z n+1 , ω) respectively represent the source wave fields on the left boundary, the middle region, and the right boundary of the (n + 1)-th layer in the frequency-wavenumber domain; n ∈ [1, N - 1].

[0293] Taking S′ mL (k x , zn+1 , ω), S' mM (k x , z n+1 , ω), S' mR (k x , z n+1 , ω) are superimposed to obtain the sound source wave field S' of the (n + 1)-th layer in the frequency-wavenumber domain m (k x , z n+1 , ω):

[0294] S' m (k x , z n+1 , ω) =

[0295] S' mL (k x , z n+1 , ω) + S' mM (k x , z n+1 , ω) + S' mR (k x , z n+1 , ω) (28)

[0296] The calculation formula of the sound source wave field action compensation operator is as follows:

[0297]

[0298] Among them, S m (x, z n+1 , ω) represents the sound source wave field of the (n + 1)-th layer in the frequency space domain; S' m (x, z n+1 , ω) represents the sound source wave field S of the (n + 1)-th layer in the frequency-wavenumber domain m (k x , z n+1 , ω) transformed into the sound source wave field in the frequency space domain; n ∈ [1, N - 1], and N represents the number of discrete layers

[0299] (4) Apply the imaging condition to the sound source wave field S m (x, z n , ω) and the received wave field P m (x, z n , ω) to obtain the imaging result of the n-th layer

[0300] Among them, the imaging condition is:

[0301]

[0302] In the formula, I m (x, zn ) represents the imaging result of the nth layer when the mth element is excited; represents S m (x, z n , ω) conjugate; n ∈ [1, N].

[0303] Full matrix capture (FMC) is a widely used capture mode for phased array probes, where each element in the phased array emits one by one, and all elements receive signals. Therefore, FMC data contains the time domain signals of each excited element-receiving element pair in the phased array. Compared with B-scan data, FMC data contains more information, greatly improving the imaging resolution and signal-to-noise ratio (SNR). Here, the wavefield extrapolation with absorbing boundary conditions (step (3)) is extended to FMC data for high-quality imaging. In this case, the wave is excited at one element in the phased array and reaches the defect at time t = τ m (x, z), and then the scattered wave is received by the phased array. Therefore, the assumptions of the explosion reflection model should be modified because the scattering point is located at the defect and is excited at t = τ m (x, z), so the defect image is the wavefield at t = τ m (x, z). Here, the sound source (incident) wavefield from the excited element is extrapolated to simulate the wave incidence and compensate for the propagation time τ m (x, z). Therefore, its imaging condition is as follows:

[0304]

[0305] where, I m (x, z) is the imaging result obtained when the mth element is excited, P m (x, z, ω) is the received wavefield extrapolated from the phased array received signal, represents S m (x, z, ω) conjugate, S m (x, z, ω) is the sound source wavefield extrapolated from the excited element. Among them, the excited wavefield of the mth element can be obtained by the following formula:

[0306] S m (x = x m , z = 0, t = 0) = δ(x = x m , t = 0) (31)

[0307] where, δ is the impulse function; x m is the position of the mth element.

[0308] (5) Repeat steps (3) and (4) until all discrete layers are traversed, obtain the imaging results of all layers, and integrate them to obtain the imaging result I when the current element is excitedm (x, z).

[0309] (6) Repeat steps (2) - (5) to obtain the imaging results when all vibration elements are excited, and after superposition, obtain the imaging result I(x, z) of the workpiece:

[0310]

[0311] Forward experiment

[0312] The imaging quality of the ultrasonic imaging method depends to a large extent on the accuracy of wave field reconstruction. Therefore, two models were designed to test the wave field reconstruction accuracy of the phase shift method. As Figure 6 (a) shows, at this time the medium is uniform, the sound speed is 1500 m / s, and a point source is located at the center of the upper boundary and is excited at the zero moment. The wave field at 15 μs is reconstructed using the one-step extrapolation method. As Figure 6 (a) shows, the true wavefront is represented by a dotted line, and it can be seen that the extrapolated wavefront is consistent with the true wavefront. However, there are serious pseudo-reflections on the side of the model. Secondly, as Figure 6 (b) shows, a three-layer medium is simulated, where the middle layer is annular, the sound speed is 2500 m / s, the sound speed of the rest of the model is 1500 m / s, and the point source is located at the center of the ring and is excited at the zero moment. Two-step extrapolation is used to construct the wave field at 14.2 μs. Here, the extrapolated wavefront is also consistent with the true wavefront, but the pseudo-reflections on the side of the model still exist. The pseudo-reflections seriously affect the accuracy of wave field reconstruction, and thus greatly affect the imaging quality of the imaging method based on phase shift.

[0313] Using Figure 11 and Figure 1 the wave field extrapolation methods with boundary absorption conditions in Figure 7 are used to perform forward simulations on the above-mentioned homogeneous and inhomogeneous media respectively, and the obtained results are as shown in Figure 7 (a) and (b) in

[0314] It can be seen from

[0315] that there are only weak pseudo-reflection waves in the homogeneous medium, and almost no pseudo-reflection waves can be found in the inhomogeneous medium. It shows that the wave field extrapolation method with boundary absorption conditions proposed by the present invention can greatly reduce the pseudo-reflection waves and improve the imaging quality. Figure 8As shown, the measuring element is immersed in water. In both cases, a wedge is used to protect the phased array from water, so a three-layer structure is formed in both cases. The measurement target of the first experiment is three void defects in an aluminum component (homogeneous medium) with a sound velocity of 6300 m / s. A wedge (SC63-NL-Z20) with a sound velocity of 2337 m / s is used, as shown in Figure 8 (a) in Figure 8 . The measurement object of the second experiment is a component (non-homogeneous medium) made of 30% glass fiber-reinforced PEEK (polyetheretherketone) with a sound velocity of 2600 m / s. The wedge used is made of ABS (acrylonitrile butadiene styrene) with a sound velocity of 2200 m / s, as shown in (b) in

[0316] The first experiment:

[0317] All interfaces in the first experiment are planar. The B-scan data processing results of the traditional imaging method (without absorption boundary conditions) are shown in Figure 9 (a) in Figure 2 . It can image the three defects, but the speckle artifacts in the background will blur the defect images. Then, the imaging method in Figure 2 is used to process the B-scan data, and the step "wavefield extrapolation with absorption boundary" in Figure 11 is replaced with the scheme in Figure 9 (b) in (b) in

[0318] To quantitatively analyze the imaging results, first introduce the total contrast index C G , which reflects the contrast between the measured target area (D f ) and the background area (N - D f ) in the image. It can quantitatively characterize the signal-to-noise ratio, and the calculation formula is as follows:

[0319]

[0320] where and are the pixel amplitudes of the target area and the background area respectively, and mean and std represent the mean value and the standard deviation respectively. Here, the image signal-to-noise ratio is proportional to C G . At the same time, another dimensionless parameter API (array performance indicator) is introduced to characterize the imaging resolution, as follows:

[0321] API = A -6dB / λ 2 (24)

[0322] where A -6dB is the area where the ratio of the peak value of the defect area image exceeds -6dB, and λ is the wave number corresponding to the center frequency. Obviously, the smaller the API, the larger the resolution.

[0323] The C Figure 9 and average API of the holes in (a) and (b) listed in Table 1 can be seen that the absorbing boundary condition in the frequency-wave number domain doubles both the signal-to-noise ratio and the resolution. G

[0324] The C G and average API of the imaging results in the first experiment of Table 1

[0325]

[0326] For the FMC (Full Matrix) data in the first experiment, using the traditional imaging method and the imaging method with absorbing boundary conditions (such as the scheme in Figure 3 ), the imaging results shown in (c) and (d) in Figure 9 are obtained. Here, the steps "extrapolation of the received wave field with absorbing boundary conditions" and "extrapolation of the sound source wave field with absorbing boundary conditions" in Figure 3 are replaced by the scheme in Figure 11 . It can be seen from the figure (c) in Figure 9 that there are many fewer artifacts in the background than in (a) (line scan data) in Figure 9 , and the signal-to-noise ratio is improved a lot. The absorbing boundary condition in the frequency-wave number domain can make the background clearer, as shown in (d) in Figure 9 . The C G and average API listed in Table 1 further illustrate that the absorbing boundary condition in the frequency-wave number domain doubles the signal-to-noise ratio and effectively improves the resolution.

[0327] Second experiment:

[0328] As shown in (b) in Figure 8 , all the interfaces in the second experiment are curved. Therefore, in the imaging method for line scan data in Figure 2 , the step "wave field extrapolation with absorbing boundary" is replaced by the scheme in Figure 1 , and the imaging result is shown in (b) in Figure 10 . The traditional imaging method (without absorbing boundary conditions) processes the B-scan data, and the result is shown in (a) in Figure 10 . It can be seen from (a) in Figure 10 that although the traditional imaging method can image three defects, the lowest defect is about to be submerged by the background artifacts. As​Figure 10 As shown in (b), the absorption boundary condition in the frequency-wavenumber domain greatly improves the signal-to-noise ratio of the bottom defect image and suppresses speckle artifacts. As shown in Table 2, Figure 10 the C of the hole defects in (a) and (b) in G and the average API indicate that the absorption boundary condition in the frequency-wavenumber domain can improve the signal-to-noise ratio and resolution.

[0329] The C of the imaging results in the second experiment in Table 2 G and the average API

[0330]

[0331] The results of processing FMC data by the traditional imaging method (without absorption boundary condition) are as shown in Figure 10 (c). It can be seen that the background artifacts around the defect image are very serious, which will affect the recognition of the bottom defect. The results of processing FMC data by the imaging method in Example 2 are as shown in Figure 10 (d). The absorption boundary condition suppresses most of the artifacts and greatly improves the signal-to-noise ratio of the defect image. As shown in Table 2, Figure 10 the C of the hole defects in (c) and (d) in G and the average API indicate that the absorption boundary condition in the frequency-wavenumber domain increases the signal-to-noise ratio by 4 times and more effectively improves the resolution.

Claims

1. A phase-shift ultrasonic imaging method with wavenumber domain absorbing boundary conditions, characterized in that, It includes the following steps: (1) Discretize the imaging area of the workpiece to form multiple discrete layers in the depth direction; (2) Convert the measurement data of the workpiece to the frequency spatial domain to obtain the surface wave field; (3) Take the surface wave field as the wave field of the first layer in the frequency spatial domain. Select any un-traversed layer as the current layer in the order from top to bottom. Apply a window selection operator to divide the wave field of the upper layer in the frequency spatial domain into a left boundary wave field, an intermediate region wave field, and a right boundary wave field; Perform a Fourier transform on the intermediate region wave field to obtain its wave field in the frequency wavenumber domain; After converting the left boundary wave field and the right boundary wave field to the frequency wavenumber domain respectively, apply corresponding filters to obtain the wave fields of the left boundary and the right boundary in the frequency wavenumber domain respectively; In the frequency wavenumber domain, extrapolate and superimpose the wave fields of the left boundary, the intermediate region, and the right boundary respectively to obtain the wave field of the current layer in the frequency wavenumber domain; Convert the wave field of the current layer in the frequency wavenumber domain to the frequency spatial domain and apply a compensation operator to obtain the wave field of the current layer in the frequency spatial domain; (4) Perform imaging based on the obtained wave field of the current layer in the frequency spatial domain to obtain the imaging result of the current layer; (5) Repeat steps (3) and (4) until all discrete layers are traversed, obtain the imaging results of all layers, and integrate them to obtain the imaging result of the workpiece.

2. The phase-shift ultrasonic imaging method with a wavenumber-domain absorbing boundary condition according to claim 1, wherein In step (3), the calculation formula for applying the window selection operator to divide the wave field of the nth layer in the frequency spatial domain into a left boundary wave field, an intermediate region wave field, and a right boundary wave field is: P L (x,z n ,ω) = P(x,z n ,ω)W L (x) P M (x,z n ,ω) = P(x,z n ,ω)W M (x) P R (x, z n , ω) = P(x, z n , ω)W R (x) Among them, P L (x, z n , ω) represents the left - boundary wave field; P M (x, z n , ω) represents the middle - region wave field; P R (x, z n , ω) represents the right - boundary wave field; P(x, z n , ω) represents the wave field of the n - th layer in the frequency - space domain, where n ∈ [1, N - 1] and N represents the number of discrete layers; W L (x) is the left boundary window selection operator, and the expression is: W M (x) is the window selection operator for the intermediate region, and the expression is: W R (x) is the right boundary window selection operator, and the expression is: where x L and x R are the left and right boundary points respectively; Δx is the discretization step length in the x direction; w() represents the window function; ω represents the frequency; x represents the position of the vibration element; z n represents the depth of the nth layer.

3. The phase-shift ultrasonic imaging method with wavenumber domain absorbing boundary condition according to claim 2, wherein The window function is one of the Hanning window, Blackman window, Flat top window, and rectangular window.

4. The phase-shift ultrasonic imaging method with a wavenumber domain absorbing boundary condition according to claim 1, characterized in that, In step (3), in the frequency wavenumber domain, use the following formula to extrapolate the wave fields of the left boundary, the intermediate region, and the right boundary respectively: where \(P(k\) x , z\) n , \(\omega)\) represents the wave field in the frequency-wavenumber domain of the \(n\)th layer; \(P'(k\) x , z\) n+1 , \(\omega)\) represents the wave field in the frequency-wavenumber domain of the \((n + 1)\)th layer; define as the phase shift operator; \(n\in[1, N - 1]\), where \(N\) represents the number of discrete layers; \(z\) n represents the depth of the \(n\)th layer; \(i\) represents the imaginary unit; \(\omega\) represents the frequency; \(\Delta z\) represents the discretization step in the depth direction; \(k\) x represents the wavenumber in the horizontal direction; represents the wavenumber in the vertical direction at depth \(z\) n , and the expression is: When the nth layer is a homogeneous medium, c n represents the sound velocity of the nth layer; When the nth layer is a non-uniform medium, c n = 1 / s0(z n ), where s0(z n ) represents the average slowness of all slownesses on the horizontal plane at a depth of z n .

5. The phase-shift ultrasonic imaging method with wavenumber-domain absorbing boundary conditions according to claim 1, characterized in that In step (3), convert the wave field of the (n + 1)th layer in the frequency wavenumber domain to the frequency spatial domain and apply a compensation operator to obtain its wave field in the frequency spatial domain; The calculation formula for applying the compensation operator is: where, P(x,z n+1 , ω) represents the wave field of the (n + 1)-th layer in the frequency spatial domain; P′(x,z n+1 , ω) represents the wave field obtained by converting the wave field of the (n + 1)-th layer in the frequency wavenumber domain to the frequency spatial domain; represents the compensation operator; i represents the imaginary unit; ω represents the frequency; x represents the position of the vibration element; Δz represents the discretization step in the depth direction; Δs(x,z n ) represents all the slowness changes on the horizontal plane at a depth of z n ; n ∈ [1, N - 1], and N represents the number of discrete layers.

6. The phase-shift ultrasonic imaging method with wavenumber-domain absorbing boundary conditions according to claim 1, characterized in that When the measured data is line-scan data, in step (4), the imaging condition for imaging based on the wave field P(x, z n , ω) obtained for the nth layer in the frequency spatial domain is as follows: I(x,z n ) = P(x,z n , t = 0) = ∫dωP(x,z n , ω) Among them, I(x, z n ) represents the imaging result of the nth layer; ω represents the frequency; x represents the position of the vibration element; z n represents the depth of the nth layer, where n ∈ [1, N], and N represents the number of discrete layers.

7. The phase-shift ultrasonic imaging method with wavenumber domain absorption boundary conditions according to claim 1, characterized in that, When the measurement data is full matrix data, the imaging processing process is as follows: Select any un-traversed element as the current element. Convert the sound source wave field and the received wave field when the current element is excited to the frequency spatial domain to obtain the surface sound source wave field and the surface received wave field when the current element is excited, and use the obtained surface sound source wave field and surface received wave field to replace the surface wave field in step (2) and enter step (3); After being processed by step (3), obtain the sound source wave field and the received wave field of the current layer in the frequency spatial domain respectively; In step (4), apply an imaging condition to the sound source wave field and the received wave field of the current layer in the frequency spatial domain to obtain the imaging result of the current layer, and enter step (5); After being processed by step (5), obtain the imaging results of all discrete layers corresponding to the current element, and after integration, obtain the imaging result generated when the current element is excited; Repeat steps (2) to (5) until all elements are traversed, obtain the imaging results of all elements, and perform superposition to obtain the imaging result of the workpiece.

8. The phase-shift ultrasonic imaging method with wavenumber-domain absorbing boundary conditions according to claim 7, characterized in that, The imaging condition for imaging from the sound source wave field and the received wave field of the nth layer in the frequency spatial domain is: Among them, I m (x, z n ) represents the imaging result of the nth layer when the mth vibration element is excited; P m (x, z n , ω) represents the received wave field of the nth layer in the frequency spatial domain; denotes the conjugate of S m (x, z n , ω), S m (x, z n , ω) represents the sound source wave field of the nth layer in the frequency spatial domain; ω represents the frequency; x represents the position of the vibration element; z n represents the depth of the nth layer, n ∈ [1, N], and N represents the number of discrete layers.

Images