A phase-shifted ultrasound imaging method with spatial absorbing boundary conditions

Abstract

The present invention provides a phase-shifted ultrasonic imaging method with a spatial domain absorption boundary condition. In the process of extrapolation, the imaging method assigns different weights to elements according to their spatial positions, divides the absorption area into two left and right boundaries and a middle area, and the weights at the left and right boundaries decrease exponentially; and in each extrapolation process, a weight operator is applied to the wave field in the frequency space domain, and the energy propagated to the boundary area is suppressed through the exponentially decaying boundary, thereby reducing pseudo-reflected waves and improving the imaging quality in terms of signal-to-noise ratio and resolution. And according to whether the discrete layer is a uniform medium, it is determined whether sound velocity compensation in the horizontal direction is required to ensure the calculation accuracy and efficiency.

Description

Technical Field

[0001] The present invention belongs to the technical field of detection methods, and in particular relates to a phase-shifted ultrasonic imaging method with a spatial domain absorption boundary condition. Background Art

[0002] Ultrasonic phased array imaging occupies a prominent position in nondestructive testing (NDT) due to its portability, reliability and relatively high accuracy. It is very effective in detecting invisible defects or damage in components. The two commonly used data modes of phased array are line scan (B-scan) and full matrix (FMC) data. Line scan data has low requirements on 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 signal of each pair of excitation and receiving elements and contains more information about the measured object. With the continuous expansion of application areas, ultrasonic phased array imaging of parts in multilayer structures has become a research hotspot. In these cases, when the sound velocity of the medium changes in the horizontal and vertical directions, it poses a greater challenge to the imaging method.

[0003] Phase shifting (PSM) is an important method in ultrasonic imaging. It regards the measured data as the wave field on the surface and reconstructs the wave field of the measured area by extrapolation in the wavenumber-frequency domain. Phase shift was first used to process line scan data to image defects in layered media. Martin improved the computational efficiency of the phase shift method by introducing Stolt interpolation. Furthermore, Wu et al. extended the phase shift to process full matrix data to image hole defects in layered media, and Ji et al. proposed an efficient phase shift method for full matrix ultrasonic imaging of layered media. In the above applications, the sound velocity of the medium only changes in the vertical direction. For the lateral variation of the sound velocity, 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 variation. Chang et al. introduced a non-stationary phase shift method to process full matrix data, which divides the wave field according to the sound velocity at each depth and uses the corresponding sound velocity for wave field extrapolation. Therefore, the phase shift-based method can process B-scan and FMC data, and can also adapt to media with vertical and lateral sound velocity variations.

[0004] Reconstructing the wave field in the measurement area and then implementing imaging conditions is the basic idea of ​​obtaining high-quality defect images, among which the first step is the core of various methods. Full wave equation migration (reverse time migration) can adapt to the complex sound velocity distribution in the medium and has high accuracy, but it requires high memory and computing time. One-way wave equation migration simplifies the full wave equation by adopting an approximate dispersion equation, which improves its efficiency, but each downward extrapolation step requires a finite difference operator, which makes the time cost still very high. In order to further improve the imaging efficiency, phase migration uses a phase shift operator to realize the downward extrapolation of the wave field, which is actually the solution of the one-way wave equation assuming that the sound velocity is laterally uniform. When performing wave equation migration, pseudo reflections will appear at the edge of the calculation area, so the absorption boundary condition is crucial for these methods to reconstruct the wave field with high accuracy and ensure 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 absorption boundary conditions for wave field extrapolation in phase migration. Summary of the invention

[0005] In order to solve the problems existing in the prior art, the present invention provides a phase-shifted ultrasonic imaging method with a spatial domain absorption boundary condition. In the process of extrapolation, the imaging method assigns different weights to elements according to their spatial positions, and divides the imaging area into left and right boundaries and a middle area. The weights at the left and right boundaries decrease exponentially. In each extrapolation process, a weight operator is applied to the wave field in the frequency space domain, and the energy propagated to the boundary area is suppressed through the exponentially decaying boundary, thereby reducing pseudo-reflected waves and improving the imaging quality in terms of signal-to-noise ratio and resolution.

[0006] A phase-shifted ultrasonic imaging method with a spatial domain absorbing boundary condition comprises the following steps:

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

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

[0009] (3) Taking the surface wave field as the wave field of the first layer in the frequency space domain, any untraversed layer is selected as the current layer in order from top to bottom, a weight operator is applied to the wave field of the previous layer in the frequency space domain, and then the wave field is converted to the frequency wave number domain for wave field extrapolation to obtain the wave field of the current layer in the frequency wave number domain;

[0010] Convert the wave field of the current layer in the frequency wave number domain to the frequency space domain and apply a compensation operator to obtain the wave field of the current layer in the frequency space domain;

[0011] (4) imaging is performed according to the obtained wave field of the current layer in the frequency space domain to obtain an imaging result of the current layer;

[0012] (5) Repeat steps (3) and (4) until all discrete layers are traversed, and the imaging results of all layers are obtained, which are then integrated to obtain the imaging results of the workpiece.

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

[0014] Preferably, in step (1), the imaging area is discretized by setting the step size in the horizontal x direction to Δx, i.e., the spacing between adjacent 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; and the number of discrete layers is N.

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

[0016]

[0017] Where i represents the imaginary unit; k x represents the wave number in the horizontal direction; ω represents the frequency; and x represents the position of the oscillator.

[0018] Preferably, in step (3), the wave field of the nth layer in the frequency wave number domain is converted to the frequency space domain and a compensation operator is applied to obtain its wave field in the frequency space domain;

[0019] The calculation formula of the action compensation operator is:

[0020]

[0021] Among them, P(x,z n+1 ,ω) represents the wave field of the n+1th layer in the frequency space domain; P′(x,z n+1 ,ω) represents the wave field obtained by converting the wave field of the nth layer in the frequency wave number domain to the frequency space domain; represents the compensation operator; i represents the imaginary unit; ω represents the frequency; x represents the position of the oscillator; Δz represents the discretization step in the depth direction; Δs(x,z n ) indicates the depth is z n (i.e., the depth of the nth layer); n∈[1,N-1], N is the number of discrete layers.

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

[0023] If the sound velocity of the nth layer varies in the horizontal direction (non-uniform medium), the wave field P(x,z) of the n+1th layer in the frequency space domain is n+1 ,ω) is obtained by converting the wave field in the frequency wavenumber domain to the frequency wavenumber domain and performing horizontal sound velocity compensation (action compensation operator).

[0024] The boundaries in the spatial domain are the two sides of the weight operator model. The weight operator includes the left boundary, the right boundary and the middle area, such as Figure 5 As shown. The weights at the two boundaries decrease exponentially, and the intermediate weight is set to 1. Preferably, in step (3), the calculation formula for the wave field action weight operator of the nth layer in the frequency space domain is:

[0025] P″(x,z n ,ω)=P(x,z n ,ω)weight(x)

[0026] Among them, P″(x,z n ,ω) represents the wave field after the weight operator is applied; P(x,z n ,ω) represents the wave field of the nth layer in the frequency space domain;

[0027] Weight(x) represents the weight operator, which is a piecewise function of x, and its expression is as follows:

[0028]

[0029] x L and x R are the left and right boundary points respectively; Δx is the discrete step length in the x direction.

[0030] The absorbing boundary condition acts on the extrapolation process of each wave field. The exponentially decaying boundary suppresses the energy propagating to the boundary area and reduces the pseudo reflection wave.

[0031] Preferably, in step (3), in the frequency wave number domain, the calculation formula for extrapolating the wave field of the nth layer in the frequency space domain after applying the weight operator to the frequency wave number domain is:

[0032]

[0033] Among them, P′(k x ,z n+1 ,ω) represents the wave field of the n+1th layer in the frequency wavenumber domain; P(k x ,z n ,ω) represents the wave field converted to the frequency wavenumber domain after the weight operator is applied; definition 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 size in the depth direction; k x Indicates the wave number in the horizontal direction; Indicates the depth is z n The total vertical wave number at is expressed as:

[0034]

[0035] When the nth layer is a uniform medium, c n represents the speed of sound at the nth layer;

[0036] When the nth layer is a non-uniform medium, c n =1 / s 0 (z n ), s 0 (z n ) indicates the depth is z n The average of all slownesses on the horizontal plane at that time.

[0037] The specific derivation process is as follows:

[0038] The propagation equation of ultrasonic wave in two-dimensional medium is:

[0039]

[0040] Among them, p is the sound pressure, c is the sound velocity of the medium, x and z are the two-dimensional spatial coordinates, and t is the time. After Fourier transform at t, the wave equation in the frequency domain is obtained as follows:

[0041]

[0042] Where ω is the frequency. Since the energy of multiple echoes is small and can be ignored, the wave equation can be decomposed into a one-way equation as follows:

[0043]

[0044] Here, the phased array probe collects signals at the surface z = 0, and the measurement targets are distributed in the half space z>0, such 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 frequency domain respectively. The wave field of the entire measurement area can be extrapolated from the surface wave field according to formula (3).

[0045] For a homogeneous medium, the speed of sound 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. Wave field extrapolation can be accomplished by applying a phase shift to each harmonic component, as shown in formula (4):

[0046]

[0047] Among them, k z is the wave number in the vertical direction, expressed as follows:

[0048]

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

[0050]

[0051] In the frequency wavenumber domain, the calculation formula for extrapolating the wave field of the current layer to obtain the wave field of the next layer is:

[0052]

[0053]

[0054] In formula (8), c n represents the speed of sound at the nth layer.

[0055] For inhomogeneous media, if there is no horizontal change, it becomes a multilayered medium, and the sound speed is only related to the depth, which can be written as c(z). In this case, extrapolation can be done by changing the phase shift with depth z.

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

[0057]

[0058] Among them, s 0(z) is the reference slowness, which is the average of all slownesses on the water surface at depth z; Δs(x,z) includes all slowness changes on the horizontal surface at depth z. The square root term in formula (3) can be Taylor expanded and approximated as:

[0059]

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

[0061]

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

[0063] Select any untraversed vibration element as the current vibration element, convert the sound source wave field and the receiving wave field when the current vibration element is excited into the frequency space domain, obtain the surface sound source wave field and the surface receiving wave field when the current vibration element is excited, and replace the surface wave field in step (2) with the obtained surface sound source wave field and surface receiving wave field to enter step (3);

[0064] After processing in step (3), the sound source wave field and receiving wave field of the current layer in the frequency space domain are obtained respectively;

[0065] In step (4), the imaging conditions are applied to the sound source wave field and the receiving wave field of the current layer in the frequency space domain to obtain the imaging result of the current layer, and then proceed to step (5);

[0066] After processing in step (5), the imaging results of all discrete layers corresponding to the current oscillator are obtained, and after integration, the imaging result generated when the current oscillator is excited is obtained;

[0067] Repeat steps (2) to (5) until all vibration elements are traversed, obtain the imaging results of all vibration elements, and superimpose them to obtain the imaging result of the workpiece.

[0068] The current oscillator is recorded as the mth oscillator, m∈[1,M], M represents the number of oscillators; the sound source wave field of the current layer (nth layer, n∈[1,N], N represents the number of discrete layers) in the frequency space domain is recorded as S m (x,z n ,ω), the received wave field is recorded as P m (x,z n ,ω). When n=1, the source wave field and receiving wave field of the first layer are obtained by converting the full matrix data.

[0069] In the frequency wavenumber domain, the extrapolation calculation formula of the received wave field is:

[0070]

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

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

[0073]

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

[0075] Since the phase shift operator and the horizontal compensation operator of the wave field extrapolation of the sound source wave field are in complex conjugate relationship with the phase shift operator and the compensation operator corresponding to the receiving wave field; therefore, in the frequency wave number domain, the extrapolation calculation formula of the sound source wave field is:

[0076]

[0077] Among them, S′ m (k x ,z n+1 ,ω) represents the sound source wave field of the n+1th 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 weight operator is applied; n∈[1,N-1], N represents the number of discrete layers.

[0078] The calculation formula of the compensation operator for the acoustic source wave field effect is:

[0079]

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

[0081] As a further preferred embodiment, the imaging condition for imaging the sound source wave field and the receiving wave field of the nth layer in the frequency space domain is:

[0082]

[0083] Among them, I m (x,z n ) represents the imaging result of the nth layer when the mth oscillator is excited; P m (x,z n ,ω) represents the received wave field of the nth layer in the frequency space domain; Indicates S m (x,z n ,ω) conjugate, S m (x,z n ,ω) represents the sound source wave field of the nth layer in the frequency space 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.

[0084] Full matrix acquisition (FMC) is a widely used acquisition mode of phased array probes, in which each element in the phased array transmits one by one and all elements receive the signal. Therefore, FMC data contains the time domain signal of each excitation element-receiver 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 receives the signal at time t = τ. m (x,z) reaches the defect, and then the scattered wave is received by the phased array. Therefore, the assumption of the explosive reflection model should be modified because the scattering point is located at the defect and at t = τ m (x,z) is excited, so the defect image is t = τ m Here, the source (incident) wave field from the excitation element is extrapolated to simulate the wave incidence and compensate for the propagation time τ m (x,z). Therefore, its imaging conditions are as follows:

[0085]

[0086] Among them, I m (x,z) is the imaging result obtained when the mth oscillator is excited, Pm (x,z,ω) is the received wave field extrapolated from the phased array received signal, Indicates S m The conjugate of (x,z,ω), S m (x,z,ω) is the sound source wave field extrapolated from the excitation element. The excitation wave field of the mth element can be obtained by the following formula:

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

[0088] Where δ is the pulse function; x m is the position of the mth oscillator.

[0089] As a further preferred embodiment, the calculation formula for superimposing the imaging results of all vibration elements is:

[0090]

[0091] Where, I(x,z) represents the imaging result of the workpiece; I m (x,z) represents the imaging result of the mth oscillator.

[0092] The imaging method for FMC data of the present invention first initializes the parameters, and the FMC data and the excitation signal are respectively recorded as P m (t,0,ω) and S m (t,0,ω) input. The algorithm has two loops. The outer loop traverses all excitation elements; the inner loop traverses all discrete depths (all discrete layers) in the imaging area. The wavefield extrapolation method with boundary absorption conditions in step (3) is applied to the receiving wavefield and the sound source wavefield in the fk domain from z simultaneously. n-1 Extrapolate to z n In the process, it is transformed to the fx domain and acts on the compensation operator and then acts on the imaging condition to obtain the imaging result I of the nth layer. m (x,z n ). After the two-layer cycle is completed, the imaging results generated when all oscillators are excited are superimposed on the imaging result I(x,z). It should be noted that the speed of sound does not need to be divided by 2 during the calculation process, because in this case both the transmission and reflection processes of the wave are taken into account.

[0093] Preferably, when the measured data is line scan data, in step (4), according to the obtained wave field P(x, z n ,ω) The imaging conditions for imaging are:

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

[0095] Among them, I(x,z n ) represents the imaging result of the nth layer; ω represents the frequency; x represents the position of the oscillator; z n Represents the depth of the nth layer, n∈[1,N], and N represents the number of discrete layers.

[0096] Line scan (B-scan) data is a basic type of data collected by ultrasonic phased array probes. It is easy to process and has low requirements on 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, so 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.

[0097] Therefore, after reconstructing the wave field in the fx domain (frequency space domain), the imaging results can be obtained by achieving the following imaging conditions:

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

[0099] Among them, I(x,z) is the imaging result.

[0100] Introducing the wavefield extrapolation method with absorbing boundary conditions (step (3)) into the imaging method of B-scan data can reduce the pseudo reflection waves existing when reconstructing the wavefield, so as to improve the imaging quality. The imaging method for line scan data of the present invention discretizes the measurement area into N layers in the depth direction and calculates the defect image layer by layer. First, the parameters are initialized, and the B-scan data P(t,0,x) is used as input. Then P(t,0,x) is transformed into the fx domain, and the extrapolation method with absorbing boundary conditions in step (3) is used to realize the wavefield from P(x,z n-1 ,ω) to P(x,z n ,ω). In addition, according to the imaging conditions, 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 to finally obtain the defect image I(x,z) of the measurement area.

[0101] The phase-shifted ultrasonic imaging method with frequency-space domain absorption boundary conditions of the present invention, the absorption boundary conditions of the frequency-space domain (fx) are to give different weights to elements according to their spatial positions during the extrapolation process. First, a weight operator is applied to the wave field of the previous layer in the frequency-space (fx) domain; the wave field of the weight operator is converted to the frequency-wavenumber domain (fk) domain for extrapolation 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 converted to the frequency-space domain and then a compensation operator is applied (when the previous layer is a uniform medium, the compensation operator can be omitted) to obtain the wave field of the current layer in the frequency-space domain, and finally the imaging result of the current layer is obtained using the imaging conditions. Finally, the imaging results of all layers are integrated to obtain the imaging result of the workpiece.

[0102] Compared with the prior art, the present invention has the following beneficial effects:

[0103] The phase-shifted ultrasonic imaging method with frequency-space domain absorption boundary conditions of the present invention assigns different weights to elements according to their spatial positions during the extrapolation process, and divides the imaging area into two left and right boundaries and an intermediate area, and the weights at the left and right boundaries decrease exponentially; and in each extrapolation process, a weight operator is applied to the wave field in the frequency-space domain, and the energy propagated to the boundary area is suppressed through the exponentially decaying boundary, thereby reducing pseudo-reflected waves and improving the imaging quality in terms of signal-to-noise ratio and resolution. And according to whether the discrete layer is a uniform medium, it is determined whether sound velocity compensation in the horizontal direction is required to ensure calculation accuracy and calculation efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0104] Figure 1 (a) is a flow chart of wave field extrapolation with absorbing boundaries of a homogeneous medium in Examples 1 and 2; (b) is a flow chart of wave field extrapolation with absorbing boundaries of an inhomogeneous medium in Examples 1 and 2;

[0105] Figure 2 This is a flow chart of imaging of line scan data in Example 1;

[0106] Figure 3 This is a flowchart of imaging of full matrix data in Example 2;

[0107] Figure 4 Schematic diagram of phased array detection;

[0108] Figure 5 is the distribution diagram of the weight function in the frequency space domain;

[0109] Figure 6 (a) is the imaging result obtained by reconstructing the wave field using the one-step extrapolation method in a homogeneous medium; (b) is the imaging result obtained by reconstructing the wave field using the two-step extrapolation method in an inhomogeneous medium.

[0110] Figure 7 (a) and (b) are respectively based on Figure 1 Methods (a) and (b) are Figure 6 Imaging results obtained by forward modeling of homogeneous media (a) and inhomogeneous media (b);

[0111] Figure 8 (a) is a hole defect detection image of a copper part; (b) is a hole defect detection image of a fiber-reinforced PEEK part;

[0112] Fig. 9 (a) is the imaging result of the traditional phase shift imaging algorithm (without absorption boundary) of line scan data; (b) is the imaging result of the hole defect of aluminum parts using the imaging method in Example 1; (c) is the imaging result of the traditional phase shift imaging algorithm (without absorption boundary) of full matrix data; (d) is the imaging result of the hole defect of aluminum parts using the imaging method in Example 2;

[0113] Fig.10 (a) The imaging result of the traditional phase shift imaging algorithm (without absorption boundary) for line scanning data; (b) The imaging result of the hole defect of the fiber-reinforced PEEK part using the imaging method in Example 1; (c) The imaging result of the traditional phase shift imaging algorithm (without absorption boundary) for full matrix data; (d) The imaging result of the hole defect of the fiber-reinforced PEEK part using the imaging method in Example 2. DETAILED DESCRIPTION

[0114] 1. Imaging method of line scan data (B scan data)

[0115] like Figure 2 As shown, a phase shift ultrasound imaging method with frequency space domain absorption boundary conditions comprises the following steps:

[0116] (1) discretizing the imaging area of ​​the workpiece to form multiple discrete layers in the depth direction;

[0117] Among them, the step size in the horizontal x direction is Δx, that is, the spacing between adjacent vibrators; 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.

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

[0119]

[0120] Among them, ω represents the frequency; x represents the position of the oscillator.

[0121] (3) The surface wave field is used as the first layer of wave field in the frequency space domain, that is, (P(x,z 1 ,ω)=P(x,z=0,ω)), and traverse the discrete layers in order from top to bottom; apply a weight operator to the wave field of the nth layer in the frequency space domain, and then convert it to the frequency wave number domain for wave field extrapolation to obtain the wave field of the n+1th layer in the frequency wave number domain.

[0122] The wave field of the n+1th layer in the frequency wave number domain is converted to the frequency space domain and then a compensation operator is applied to obtain the wave field of the n+1th layer in the frequency space domain.

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

[0124] The specific process is:

[0125] Action weight operator:

[0126] The boundaries in the spatial domain are the two sides of the weight operator model. The weight operator includes the left boundary, the right boundary and the middle area, such as Figure 5 As shown in Figure 2, the weights at the two boundaries decrease exponentially, and the intermediate weights are set to 1.

[0127] For the nth layer of wave field P(x,z n ,ω) acts on the weight operator to obtain P″(x,z n ,ω):

[0128] P″(x,z n ,ω)=P(x,z n ,ω)weight (x) (2)

[0129] Among them, weight(x) represents the weight operator, which is a piecewise function about x, and the expression is as follows:

[0130]

[0131] x L and x R are the left and right boundary points respectively; Δx is the discrete step length in the x direction.

[0132] Extrapolation of the wave field in the frequency-wavenumber domain:

[0133] P″(x,z n ,ω) is converted to the frequency wavenumber domain, and P(k x ,z n ,ω); for P(k x ,z n ,ω) to obtain the wave field P′(k x ,zn+1 ,ω);

[0134] The extrapolation formula is:

[0135]

[0136] In formula (4), 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 size in the depth direction; k x Indicates the wave number in the horizontal direction; Indicates the depth is z n The total vertical wave number at is expressed as:

[0137]

[0138] When the nth layer is a uniform medium, c n represents the speed of sound at the nth layer;

[0139] When the nth layer is a non-uniform medium, c n =1 / s 0 (z n ), s 0 (z n ) indicates the depth is z n The average of all slownesses on the horizontal plane at that time.

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

[0141] The propagation equation of ultrasonic wave in two-dimensional medium is:

[0142]

[0143] Among them, p is the sound pressure, c is the sound velocity of the medium, x and z are the two-dimensional spatial coordinates, and t is the time. After Fourier transform at t, the wave equation in the frequency domain is obtained as follows:

[0144]

[0145] Where ω is the frequency. Since the energy of multiple echoes is small and can be ignored, the wave equation can be decomposed into a one-way equation as follows:

[0146]

[0147] Here, the phased array probe collects signals at z = 0 on the surface during measurement, and the measurement targets are distributed in the half space where z>0, such as Figure 4As shown. Therefore, the wave field on the surface can be directly acquired 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 frequency domain respectively. The wave field of the entire measurement area can be extrapolated from the surface wave field according to formula (8).

[0148] For homogeneous media, the speed of sound 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. Wave field extrapolation can be accomplished by applying a phase shift to each harmonic component, as shown in formula (9):

[0149]

[0150] Among them, k z is the wave number in the vertical direction, which is expressed as follows:

[0151]

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

[0153]

[0154] In the frequency wavenumber domain, the calculation formula for obtaining the wave field of the current layer by extrapolating the wave field of the previous layer in a layer-by-layer recursive manner is:

[0155]

[0156]

[0157] In formula (13), c n Represents the speed of sound in the nth layer, which is ∈[1,N-1].

[0158] For inhomogeneous media, if there is no horizontal change, it becomes a multilayered medium, and the sound speed is only related to the depth, which can be written as c(z). In this case, extrapolation can be done by changing the phase shift with depth z.

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

[0160]

[0161] Among them, s 0(z) is the reference slowness, which is the average of all slownesses on the water surface at depth z; Δs(x,z) includes all slowness changes on the horizontal surface at depth z. The square root term in formula (8) can be Taylor expanded and approximated as:

[0162]

[0163] If s 0 (z)>>Δs(x,z), then in equation (13), c n =1 / s 0 (z n ), formula (13) can be written as:

[0164]

[0165] Compensation of sound velocity in the horizontal direction:

[0166] The extrapolated P′(k x ,z n+1 ,ω) is converted to the frequency space domain to obtain P′(x,z n+1 ,ω), the post-action compensation operator, and the wave field P(x,z) of the n+1th layer in the frequency space domain is obtained. n+1 ,ω);

[0167] The calculation formula is:

[0168]

[0169] in, represents the compensation operator; i represents the imaginary unit; ω represents the frequency; x represents the position of the oscillator; Δz represents the discretization step in the depth direction; Δs(x,z n ) indicates the depth is z n All slowness changes on the horizontal plane when ; n∈[1,N-1], N is the number of discrete layers.

[0170] If the sound velocity of the nth layer does not change in the horizontal direction (homogeneous medium), that is, Δs(x,z n )=0, the horizontal compensation can be omitted; therefore, before traversing the discrete layers, it is possible to first determine whether the nth layer is a homogeneous medium. If so, the wave field P(x,z) of the n+1th layer in the frequency space domain is n+1 ,ω) can be directly obtained by transforming its wave field in the frequency-wavenumber domain, such as Figure 1 As shown in (a).

[0171] If the sound velocity of the nth layer varies in the horizontal direction (non-uniform medium), the wave field P(x,z) of the n+1th layer in the frequency space domain is n+1,ω) is obtained by converting the wave field in the frequency wave number domain to the frequency wave number domain and performing horizontal sound velocity compensation (action compensation operator), as Figure 1 As shown in (b).

[0172] (4) imaging is performed according to the obtained wave field of the nth layer in the frequency space domain to obtain an imaging result of the nth layer;

[0173] The imaging conditions are:

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

[0175] At this time, I(x,z n ) is the imaging result of the nth layer; n∈[1,N].

[0176] Line scan (B-scan) data is a basic type of data collected by ultrasonic phased array probes. It is easy to process and has low requirements on 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, so 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.

[0177] Therefore, after reconstructing the wave field in the fk domain (frequency wavenumber domain), the imaging results can be obtained by achieving the following imaging conditions:

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

[0179] Among them, I(x,z) is the imaging result.

[0180] (5) Repeat steps (3) and (4) until all discrete layers are traversed, and the imaging results of all layers are obtained, which are then integrated to obtain the imaging results of the workpiece.

[0181] Introducing the wave field extrapolation method with absorbing boundary conditions (step (3)) into the imaging method of B-scan data can reduce the pseudo reflection waves existing when reconstructing the wave field, thereby improving the imaging quality.

[0182] 2. Imaging method of full matrix data (FMC data)

[0183] like Figure 3 As shown, a phase-shifted ultrasonic imaging method with a spatial domain absorbing boundary condition comprises the following steps:

[0184] (1) discretizing the imaging area of ​​the workpiece to form multiple discrete layers in the depth direction;

[0185] Among them, the step size in the horizontal x direction is Δx, that is, the spacing between adjacent vibrators; 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.

[0186] (2) Select any untraversed vibration element as the current vibration element m, and set the sound source wave field S when the current vibration element m is excited. m (x,z=0,t) and the received wave field P m (x,z=0,t) is converted to the frequency space domain to obtain the surface acoustic source wave field S in the frequency space domain when the current vibration element is excited m (x,z=0,ω) and the surface receiving wave field P m (x,z=0,ω);

[0187] Take receiving TRON as an example:

[0188]

[0189] Among them, ω represents the frequency; x represents the position of the oscillator.

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

[0191] (3) The surface sound source wave field and the surface receiving wave field are used as the first layer of the sound source wave field and the receiving wave field in the frequency space domain, that is, (S m (x,z=0,ω)=S m (x,z 1 ,ω)、P m (x,z=0,ω)=P m (x,z 1 ,ω)), traverse the discrete layers in order from top to bottom; apply weight operators to the sound source wave field and receiving wave field of the nth layer in the frequency space domain, and then convert them to the frequency wave number domain for wave field extrapolation to obtain the sound source wave field and receiving wave field of the n+1th layer in the frequency wave number domain;

[0192] The sound source wave field and receiving wave field of the n+1th layer in the frequency wave number domain are converted to the frequency space domain respectively and then acted on by the compensation operator to obtain the sound source wave field and receiving wave field of the n+1th layer in the frequency space domain.

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

[0194] Taking the receiving wave field as an example, the specific derivation process is as follows:

[0195] Action weight operator:

[0196] The boundaries in the spatial domain are the two sides of the weight operator model. The weight operator includes the left boundary, the right boundary and the middle area, such as Figure 5 As shown in Figure 2, the weights at the two boundaries decrease exponentially, and the intermediate weights are set to 1.

[0197] The received wave field P of the nth layer in the frequency space domain m (x,z n ,ω) acts on the weight operator to obtain P″ m (x,z n ,ω):

[0198] P″ m (x,z n ,ω)=P m (x,z n ,ω)weight(x) (2)

[0199] Among them, weight(x) represents the weight operator, which is a piecewise function about x, and the expression is as follows:

[0200]

[0201] x L and x R are the left and right boundary points respectively; Δx is the discrete step length in the x direction.

[0202] Extrapolation of the wave field in the frequency-wavenumber domain:

[0203] P″ m (x,z n ,ω) is converted to the frequency wavenumber domain to obtain P m (k x ,z n ,ω); for P m (k x ,z n ,ω) to obtain the next layer of wave field P′ in the frequency wave number domain m (k x ,z n+1 ,ω); the extrapolation formula is:

[0204]

[0205] In formula (4), 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 size in the depth direction; k x Indicates the wave number in the horizontal direction; Indicates the depth is z nThe total vertical wave number at is expressed as:

[0206]

[0207] When the nth layer is a uniform medium, c n represents the speed of sound at the nth layer;

[0208] When the nth layer is a non-uniform medium, c n =1 / s 0 (z n ), s 0 (z n ) indicates the depth is z n The average of all slownesses on the horizontal plane at that time.

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

[0210] The propagation equation of ultrasonic wave in two-dimensional medium is:

[0211]

[0212] Among them, p is the sound pressure, c is the sound velocity of the medium, x and z are the two-dimensional spatial coordinates, and t is the time. After Fourier transform at t, the wave equation in the frequency domain is obtained as follows:

[0213]

[0214] Where ω is the frequency. Since the energy of multiple echoes is small and can be ignored, the wave equation can be decomposed into a one-way equation as follows:

[0215]

[0216] Here, the phased array probe collects signals at z = 0 on the surface during measurement, and the measurement targets are distributed in the half space where z>0, such as Figure 4 As shown. Therefore, the wave field on the surface can be directly acquired 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 frequency domain respectively. The wave field of the entire measurement area can be extrapolated from the surface wave field according to formula (8).

[0217] For homogeneous media, the speed of sound 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. Wave field extrapolation can be accomplished by applying a phase shift to each harmonic component, as shown in formula (9):

[0218]

[0219] Among them, k zis the wave number in the vertical direction, which is expressed as follows:

[0220]

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

[0222]

[0223] In the frequency wavenumber domain, the calculation formula for extrapolating the wave field of the current layer to obtain the wave field of the next layer is:

[0224]

[0225]

[0226] In formula (13), c n represents the speed of sound at the nth layer.

[0227] For inhomogeneous media, if there is no horizontal change, it becomes a multilayered medium, and the sound speed is only related to the depth, which can be written as c(z). In this case, extrapolation can be done by changing the phase shift with depth z.

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

[0229]

[0230] Among them, s 0 (z) is the reference slowness, which is the average of all slownesses on the water surface at depth z; Δs(x,z) includes all slowness changes on the horizontal surface at depth z. The square root term in formula (8) can be Taylor expanded and approximated as:

[0231]

[0232] If s 0 (z)>>Δs(x,z), then in equation (13), c n =1 / s 0 (z n ), formula (13) can be written as:

[0233]

[0234] Action compensation operator:

[0235] The extrapolated P′ m (k x ,z n+1 ,ω) is converted to the frequency space domain to obtain P′ m (x,z n+1 ,ω), the post-action compensation operator, obtains the wave field P of the current layer in the frequency space domain m (x,z n+1 ,ω);

[0236] The calculation formula is:

[0237]

[0238] in, represents the compensation operator; i represents the imaginary unit; ω represents the frequency; x represents the position of the oscillator; Δz represents the discretization step in the depth direction; Δs(x,z n ) indicates the depth is z n All slowness changes on the horizontal plane when ; n∈[1,N-1], N is the number of discrete layers.

[0239] If the sound velocity of the nth layer does not change in the horizontal direction (homogeneous medium), that is, Δs(x,z n )=0, the horizontal compensation can be omitted; therefore, before traversing the discrete layers, it is possible to first determine whether the nth layer is a homogeneous medium. If so, the wave field P(x,z) of the n+1th layer in the frequency space domain is n+1 ,ω) can be directly obtained by transforming its wave field in the frequency-wavenumber domain, such as Figure 1 As shown in (a).

[0240] If the sound velocity of the nth layer varies in the horizontal direction (non-uniform medium), the wave field P(x,z) of the n+1th layer in the frequency space domain is n+1 ,ω) is obtained by converting the wave field in the frequency wave number domain to the frequency wave number domain and performing horizontal sound velocity compensation (action compensation operator), as Figure 1 As shown in (b).

[0241] For the sound source wave field, the phase shift operator and the horizontal compensation operator of its wave field extrapolation are in complex conjugate relationship with the phase shift operator and compensation operator corresponding to the receiving wave field; therefore, in the frequency wave number domain, the extrapolation calculation formula of the sound source wave field is:

[0242]

[0243] At this time, S′ m (k x ,z n+1 ,ω) represents the sound source wave field of the n+1th layer in the frequency wavenumber domain; S m (kx ,z n ,ω) represents the sound source wave field converted to the frequency wavenumber domain after the weight operator is applied; n∈[1,N-1].

[0244] The calculation formula of the compensation operator for the acoustic source wave field effect is:

[0245]

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

[0247] (4) For the sound source wave field S of the nth layer in the frequency space domain m (x,z n ,ω) and the received wave field P m (x,z n ,ω), after applying the imaging conditions, the imaging result of the nth layer is obtained.

[0248] The imaging conditions are:

[0249]

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

[0251] Full matrix acquisition (FMC) is a widely used acquisition mode of phased array probes, in which each element in the phased array transmits one by one and all elements receive the signal. Therefore, FMC data contains the time domain signal of each excitation element-receiver 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 receives the signal at time t = τ. m (x,z) reaches the defect, and then the scattered wave is received by the phased array. Therefore, the assumption of the explosive reflection model should be modified because the scattering point is located at the defect and at t = τ m (x,z) is excited, so the defect image is t = τm Here, the source (incident) wave field from the excitation element is extrapolated to simulate the wave incidence and compensate for the propagation time τ m (x,z). Therefore, its imaging conditions are as follows:

[0252]

[0253] Among them, I m (x,z) is the imaging result obtained when the mth oscillator is excited, P m (x,z,ω) is the received wave field extrapolated from the phased array received signal, Indicates S m The conjugate of (x,z,ω), S m (x,z,ω) is the sound source wave field extrapolated from the excitation element. The excitation wave field of the mth element can be obtained by the following formula:

[0254] S m (x=x m ,z=0,t=0)=δ(x=x m ,t=0) (21)

[0255] Where δ is the pulse function; x m is the position of the mth oscillator.

[0256] (5) Repeat steps (3) and (4) until all discrete layers are traversed, and the imaging results of all layers are obtained, and then integrated to obtain the imaging result I when the current oscillator is excited. m (x,z).

[0257] (6) Repeat steps (2) to (5) to obtain the imaging results when all the oscillators are excited, and then superimpose them to obtain the imaging result I (x, z) of the workpiece:

[0258]

[0259] Forward Experiment

[0260] The imaging quality of ultrasound imaging methods depends largely on the accuracy of wave field reconstruction. Therefore, two models are designed to test the wave field reconstruction accuracy of the phase shift method. Figure 6 As shown in (a), the medium is uniform, the sound speed is 1500m / s, and a point source is located at the center of the upper boundary and excited at time zero. The wave field at 15us is reconstructed using the one-step extrapolation method. Figure 6 As shown in (a), the real wavefront is represented by a dotted line. It can be seen that the extrapolated wavefront is consistent with the real wavefront. However, there is a serious pseudo reflection on the side of the model. Figure 6As shown in (b), a three-layer medium is simulated, in which the middle layer is annular with a sound speed of 2500m / s, and the sound speed of the rest of the model is 1500m / s. The point source is located at the center of the ring and excited at time zero. Two-step extrapolation is used to construct the 14.2us wavefield. Here, the extrapolated wavefront is also consistent with the real wavefront, but the pseudo-reflection on the side of the model still exists. Pseudo-reflection seriously affects the accuracy of wavefield reconstruction, which in turn greatly affects the imaging quality of the imaging method based on phase shift.

[0261] use Figure 1 The wave field extrapolation method with boundary absorption conditions in (a) and (b) is used to model the above homogeneous medium and inhomogeneous medium respectively. The results are as follows: Figure 7 As shown in (a) and (b). Figure 7 It can be seen that the amplitudes of the pseudo reflection waves in the two media are greatly reduced, indicating that the wave field extrapolation method with boundary absorption conditions proposed in the present invention can effectively reduce pseudo reflection waves and improve imaging quality.

[0262] Detection experiment

[0263] Experiments were conducted to evaluate the imaging method proposed in the embodiment of the present invention. A 64-element phased array probe (Shantou Ultrasound) with a center frequency of 2.5 MHz and an element spacing of 0.75 mm was used in the experiment. Figure 8 As shown, the measuring element is immersed in water. In both cases, wedges are 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 cavity defects in an aluminum component (homogeneous medium) with a sound velocity of 6300m / s. A wedge (SC63-NL-Z20) with a sound velocity of 2337m / s is used. Figure 8 The second experiment measured a component made of 30% glass fiber reinforced PEEK (polyetheretherketone) with a sound velocity of 2600 m / s (inhomogeneous medium), and the wedge used was made of ABS (acrylonitrile butadiene styrene) with a sound velocity of 2200 m / s, as shown in (a). Figure 8 As shown in (b), the B-scan and FMC data were acquired using a data acquisition card 64 / 64OEM-PA (AOS. Ltd, America) with a sampling frequency of 50 MHz and a time range of 60 μs.

[0264] First experiment:

[0265] In the first experiment, all interfaces were plane. The B-scan data processing results of the traditional imaging method (without absorbing boundary conditions) are as follows: Fig. 9 It can image the three defects, but the speckle artifacts in the background will blur the defect image. Then use Figure 2 The imaging method in the B-scan data is used to process the Figure 2 The step "wave field extrapolation with absorbing boundaries" is replaced by Figure 1 The imaging results are shown in (a). Fig. 9 As shown in (b), its energy is more concentrated at the defect, that is, the absorbing boundary condition improves the signal-to-noise ratio and makes the defect image clearer.

[0266] In order to quantitatively analyze the imaging results, the total contrast index C is first introduced G , which reflects the measurement target area in the image (D f ) and background area (ND f ), which can quantitatively characterize the signal-to-noise ratio, and the calculation formula is as follows:

[0267]

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

[0269] API = A -6dB / λ 2 (twenty four)

[0270] Among them, A -6dB is the area whose ratio to 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 greater the resolution.

[0271] Listed in Table 1 Fig. 9 C of the holes in (a) and (b) G and average API, it can be seen that the frequency-space domain absorption boundary condition significantly improves the signal-to-noise ratio and resolution.

[0272] Table 1 C of imaging results in the first experiment G Tie API

[0273]

[0274] For the FMC (full matrix) data in the first experiment, conventional imaging methods and imaging methods with absorbing boundary conditions (such as Figure 3 ), we get Fig. 9 The imaging results shown in (c) and (d) are shown in Fig. Figure 3The steps "Extrapolation of receiving wave field with absorbing boundary conditions" and "Extrapolation of source wave field with absorbing boundary conditions" are replaced by Figure 1 The solution in (a). Fig. 9 As can be seen in Figure (c), the background artifacts are smaller than Fig. 9 There are much fewer artifacts in (a) (line scan data), and the signal-to-noise ratio is much improved. The absorbing boundary conditions in the frequency space domain suppress most of the artifacts, such as Fig. 9 As shown in (d). The C listed in Table 1 G The average API further illustrates that the absorption boundary condition in the frequency space domain greatly improves the signal-to-noise ratio.

[0275] Second experiment:

[0276] like Figure 8 As shown in (b), all interfaces in the second experiment are curved. Figure 2 In the imaging method for line scan data, replace the step "wavefield extrapolation with absorbing boundaries" with Figure 1 In the scheme (b), the imaging results are as follows Fig.10 The conventional imaging method (without absorbing boundary conditions) processes the B-scan data, and the result is as follows Fig.10 As shown in (a). Fig.10 As can be seen in (a), although the traditional imaging method can image three defects, the lowest defect will be submerged by the background artifacts. Fig.10 As shown in (b), the absorption boundary condition greatly improves the signal-to-noise ratio of the bottom defect image and suppresses the speckle artifact. As shown in Table 2, Fig.10 C of the hole defect in (a) and (b) G and average API show that the absorbing boundary conditions can improve the signal-to-noise ratio and resolution.

[0277] Table 2 C of imaging results in the second experiment G Tie API

[0278]

[0279] The results of processing FMC data using traditional imaging methods (without absorbing boundary conditions) are as follows: Fig.10 As shown in (c) of the figure, 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 result of processing the FMC data using the imaging method in Example 2 is shown in Fig.10 As shown in (d), the absorbing boundary condition suppresses most artifacts and greatly improves the signal-to-noise ratio of the defect image. As shown in Table 2, Fig.10 C of the hole defect in (c) and (d) GThe average API shows that the frequency-space domain absorption boundary condition improves the signal-to-noise ratio by a factor of 2 and improves the resolution more effectively.

Claims

1. A phase-shifted ultrasonic imaging method with spatial domain absorbing boundary conditions, characterized in that: The following steps are involved: (1) discretizing the imaging area of ​​the workpiece to form multiple discrete layers in the depth direction; (2) Convert the workpiece measurement data into the frequency space domain to obtain the surface wave field; (3) Taking the surface wave field as the wave field of the first layer in the frequency space domain, any untraversed layer is selected as the current layer in order from top to bottom, a weight operator is applied to the wave field of the previous layer in the frequency space domain, and then the wave field is converted to the frequency wave number domain for wave field extrapolation to obtain the wave field of the current layer in the frequency wave number domain; Convert the wave field of the current layer in the frequency wave number domain to the frequency space domain and apply a compensation operator to obtain the wave field of the current layer in the frequency space domain; (4) imaging is performed according to the obtained wave field of the current layer in the frequency space domain to obtain an imaging result of the current layer; (5) Repeat steps (3) and (4) until all discrete layers are traversed, the imaging results of all layers are obtained, and the imaging results of the workpiece are obtained by integration; In step (3), the calculation formula for the wave field action weight operator of the nth layer in the frequency space domain is: P″(x,z n ,ω)=P(x,z n ,ω)weight(x) Among them, P″(x,z n ,ω) represents the wave field after the weight operator is applied; P(x,z n ,ω) represents the wave field of the nth layer in the frequency space domain; Weight(x) represents the weight operator, which is a piecewise function of x, and its expression is as follows: x L and x R are the left and right boundary points respectively; Δx is the discrete step length in the x direction; In step (3), in the frequency wave number domain, the calculation formula for extrapolating the wave field of the nth layer in the frequency space domain after applying the weight operator to the frequency wave number domain is: Among them, P ′ (k x ,z n+1 ,ω) represents the wave field of the nth layer in the frequency wavenumber domain; P(k x ,z n ,ω) represents the wave field converted to the frequency wavenumber domain after the weight operator is applied; definition 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 size in the depth direction; k x Indicates the wave number in the horizontal direction; Indicates the depth is z n The vertical wave number at is expressed as: When the nth layer is a uniform medium, c n represents the speed of sound at the nth layer; When the nth layer is a non-uniform medium, c n =1 / s0(z n ), s0(z n ) indicates the depth is z n The average of all slownesses on the horizontal plane at that time.

2. The phase-shifted ultrasonic imaging method with spatial-domain absorbing boundary conditions according to claim 1, characterized in that: In step (3), the wave field of the n+1th layer in the frequency wave number domain is converted to the frequency space domain and a compensation operator is applied to obtain its wave field in the frequency space domain; The calculation formula of the action compensation operator is: Among them, P(x,z n+1 ,ω) represents the wave field of the n+1th layer in the frequency space domain; P ′ (x,z n+1 ,ω) represents the wave field obtained by converting the wave field of the n+1th 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 oscillator; Δz represents the discretization step in the depth direction; Δs(x,z n ) indicates the depth is z n All slowness changes on the horizontal plane when ; n∈[1,N-1], N is the number of discrete layers.

3. The phase-shifted ultrasonic imaging method with spatial-domain absorbing boundary conditions according to claim 1, characterized in that: When the measured data is line scan data, in step (4), according to the obtained wave field P(x, z n ,ω) The imaging conditions for imaging are: Among them, I(x,z n ) represents the imaging result of the nth layer; ω represents the frequency; x represents the position of the oscillator; z n Represents the depth of the nth layer, n∈[1,N], and N represents the number of discrete layers.

4. The phase-shifted ultrasonic imaging method with spatial-domain absorbing boundary conditions according to claim 1, characterized in that: When the measured data is full matrix data, the imaging process is as follows: Select any untraversed vibration element as the current vibration element, convert the sound source wave field and the receiving wave field when the current vibration element is excited into the frequency space domain, obtain the surface sound source wave field and the surface receiving wave field when the current vibration element is excited, and replace the surface wave field in step (2) with the obtained surface sound source wave field and surface receiving wave field to enter step (3); After processing in step (3), the sound source wave field and receiving wave field of the current layer in the frequency space domain are obtained respectively; In step (4), the imaging conditions are applied to the sound source wave field and the receiving wave field of the current layer in the frequency space domain to obtain the imaging result of the current layer, and then proceed to step (5); After processing in step (5), the imaging results of all discrete layers corresponding to the current oscillator are obtained, and after integration, the imaging result generated when the current oscillator is excited is obtained; Repeat steps (2) to (5) until all vibration elements are traversed, obtain the imaging results of all vibration elements, and superimpose them to obtain the imaging result of the workpiece.

5. The phase-shifted ultrasonic imaging method with spatial-domain absorbing boundary conditions according to claim 4, characterized in that: The imaging condition for imaging by the sound source wave field and the receiving wave field of the nth layer in the frequency space domain is: Among them, I m (x,z n ) represents the imaging result of the nth layer when the nth oscillator is excited; P m (x,z n ,ω) represents the received wave field of the nth layer in the frequency space domain; Indicates S m (x,z n ,ω) conjugate, S m (x,z n ,ω) represents the sound source wave field of the nth layer in the frequency space 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