A high-efficiency full-matrix imaging method for curved parts using ultrasonic phased array

Abstract

The present invention provides a full-matrix efficient imaging method of an ultrasonic phased array for curved parts, which performs classification, extrapolation and imaging according to whether the sound velocity changes in the horizontal direction of each layer. For a uniform dielectric layer, the wave field of the layer in the frequency wavenumber domain obtained by extrapolation is directly used for imaging; for a non-uniform dielectric layer, it is necessary to convert the wave field of the layer in the frequency wavenumber domain obtained by extrapolation into the frequency space domain and then apply a compensation operator in the horizontal direction to obtain the wave field of the current layer in the frequency space domain, and then image the wave field of the current layer in the frequency space domain. The above-mentioned classification imaging method can not only ensure the quality of the imaging results of the non-uniform dielectric layer, but also improve the overall imaging efficiency. In addition, there are only two operators, superposition and Fourier transform, in the imaging conditions of the two cases, which greatly reduces the memory requirements and calculation complexity, and shortens the calculation time.

Description

Technical Field

[0001] The invention belongs to the technical field of detection methods, and in particular relates to a full-matrix efficient imaging method of an ultrasonic phased array for curved surface parts. Background Art

[0002] Curved components are common in various industries such as manufacturing, automotive, energy, and aerospace. Among them, shafts and camshafts are core components in drive systems. However, complex manufacturing processes and harsh use environments make components prone to holes and cracks, which greatly reduce the mechanical properties of components and even lead to unexpected failures and serious accidents. Therefore, efficient and accurate nondestructive testing (NDT) of defects in curved components is essential to ensure their safety and reliability in use. Ultrasonic phased array has gained a prominent position in industrial nondestructive testing due to its portability, reliability, and high precision. Conventional ultrasonic phased array technology generates coherent ultrasonic beams by controlling the time delay of the transmitting oscillator. However, the limited information in this data modality leads to low image contrast of the target object, which seriously hinders its application in curved components. Full matrix capture (FMC) circumvents this defect, which captures the complete time domain signal of each pair of transmitting and receiving oscillators. FMC datasets contain more information and show great potential for solving various problems.

[0003] Flexible phased arrays are designed to accommodate curved components, but they are expensive and in their infancy, so non-contact methods through water immersion are more widely adopted. However, the presence of curved interfaces causes sound velocity variations in the measurement area in both horizontal and vertical directions, so immersion methods place higher demands on the post-processing of FMC data sets. To overcome this problem, there are several imaging methods that are adapted to inhomogeneous media. The ray-based total focusing method (TFM) simplifies the beam into rays and uses the Fermat principle to determine the ray path of the ultrasound, thereby estimating the TFM delay. However, the computational complexity of the ray tracing process grows exponentially with the number of layers, so it is very time-consuming when imaging multi-layer media.

[0004] Reverse time migration (RTM) uses the full wave equation to reconstruct the wave field of the measurement area. It has the advantages of high imaging accuracy, no inclination angle restriction, and adaptability to arbitrarily complex velocity models. Rao et al. used elastic reverse time migration (ERTM) to image irregular grooves in metal blocks. However, ERTM is performed in the time domain and requires iteration at each time node, which requires a lot of calculations and long imaging time. Frequency domain reverse time migration (FD-RTM) can introduce ultrasound full matrix imaging. This method can greatly improve the efficiency of RTM because all excitations of the FMC data set can be extrapolated together to simultaneously obtain the source wave field and the received wave field of each excitation.

[0005] In addition, Lukomski (PSM-based method) extended the phase-shift-plus-interpolation (PSPI) method to FMC data for ultrasonic imaging of irregular geometries. Here, the sound velocity field with lateral sound velocity changes is decomposed into several uniform sound velocity fields according to the sound velocity values ​​in the horizontal direction, and phase shift migration (PSM) is implemented in each of these sound velocity fields, and then superimposed to obtain the extrapolated wave field. The PSM-based method has a great advantage over the previous two methods in imaging time, but its wave field reconstruction process and imaging conditions limit the computational efficiency. In addition, PSPI fails when the surface of the component is so complex that the sound velocity at the interface is difficult to decompose. Summary of the invention

[0006] In order to solve the problems existing in the prior art, the present invention proposes an efficient full-matrix imaging method for ultrasonic phased array of curved parts. The method considers the downward wavefield extrapolation method of FMC data as the combined effect of the downward displacement of the receiving wavefield and the sound source wavefield. The split-step extrapolation method (step-by-step extrapolation) is introduced into the extrapolation process of the receiving wavefield and the sound source wavefield to compensate for the lateral changes in the sound velocity, so that the wavefield at each depth of the imaging area can be obtained more accurately, thereby obtaining high-quality imaging results. Moreover, this method has obvious advantages in computational efficiency over other methods and has high imaging efficiency.

[0007] A method for high-efficiency full-matrix ultrasonic phased array imaging of curved parts, comprising the following steps:

[0008] (1) Discretizing the imaging area of ​​the curved part to form multiple discrete layers in the depth direction;

[0009] (2) Convert the full matrix data of the curved surface part into the frequency wavenumber domain through Fourier transform to obtain the surface wave field;

[0010] (3) Taking the surface wave field as the wave field of the first layer, select any untraversed discrete layer as the current layer in order from top to bottom, and determine whether the sound velocity of the previous layer changes in the horizontal direction:

[0011] If yes, in the frequency-wavenumber domain, the wave field of the previous layer is extrapolated downward and transformed to the frequency-space domain, and then a compensation operator is applied to obtain the wave field of the current layer in the frequency-space domain; the imaging condition and Fourier transform are applied to the wave field of the current layer in the frequency-space domain to obtain the imaging result of the current layer and its wave field in the frequency-wavenumber domain;

[0012] If not, then in the frequency-wavenumber domain, the wave field of the previous layer is extrapolated downward to obtain the wave field of the current layer, and the imaging conditions are applied to obtain the imaging result of the current layer;

[0013] (4) Repeat step (3) until all layers are traversed, obtain the imaging results of all layers, and integrate them to obtain the imaging results of the curved surface parts.

[0014] Preferably, in step (1), the imaging area is discretized by setting the number of discrete layers to N and the step length in the depth z direction to Δz, wherein the size of Δz is determined according to the resolution requirement in the vertical direction.

[0015] Preferably, in step (2), full matrix acquisition (FMC) acquires the time domain signals between all excitation-receiving pairs in the phased array probe. For an n-element phased array probe, the FMC data contains n 2 A scan signal. D(x r , x s , t) represents the signal transmitted from the sth oscillator and received by the rth oscillator at time t, such as Figure 2 As shown, the ultrasonic wave moves from the oscillator position (x s , 0), and after being reflected by the scattering point at (x, z), it is sent out by the r , received by the sensor at 0).

[0016] The traditional measurement method is to scan the surface with a probe to obtain measurement data. Here, D(x, t) represents the one-dimensional scan data, which can be directly regarded as the wave field on the surface P(x, z = 0, t). Compared with the scan data, FMC data has an extra dimension. Since both the excitation element and the receiving element have two-dimensional coordinates in space, in this case the wave field can be regarded as five-dimensional and expressed as P(x r , z r , x s , z s , t) represents, where x r and z r are the horizontal and vertical coordinates of the receiving element; x s and z s are the horizontal and vertical coordinates of the exciting element. Therefore, the wave field on the measurement surface (surface wave field) is P(x r ,0,x s , 0, t), and:

[0017] P(x r ,0,x s ,0,t)=D(x r , x s , t) 1)

[0018] P(x r ,0,x s , 0, t) can be obtained by a three-dimensional Fourier transform to obtain the surface wave field in the frequency wavenumber domain:

[0019]

[0020] Among them, k sx and k rx are the horizontal wave numbers of the exciting oscillator and the receiving oscillator respectively; x s and x r are the positions of the exciting oscillator and the receiving oscillator respectively; i is the imaginary unit; ω is the frequency; t is the time.

[0021] Preferably, in step (3), if the sound velocity of the previous layer does not change in the horizontal direction (uniform layer), the wave field of the current layer is obtained by extrapolating the wave field of the previous layer, that is, the wave field of the first layer is first extrapolated to obtain the wave field of the second layer; then the wave field of the second layer is extrapolated to obtain the wave field of the third layer; and so on, the wave fields of all layers are obtained. Among them, the wave field of the first layer is the surface wave field. By this layer-by-layer extrapolation method, a more accurate wave field of each layer can be obtained, thereby improving the final imaging quality. If the sound velocity of the previous layer changes in the horizontal direction (non-uniform layer), the wave field of the previous layer in the frequency wavenumber domain is used for extrapolation processing and imaging according to the change.

[0022] When the sound velocity c in the horizontal direction is uniform, the wave field in the measurement area can be directly obtained by downward extrapolation of the surface wave field in the frequency wavenumber domain. In order to obtain the wave field at depth z, the receiving element is first offset:

[0023]

[0024] Among them, k rz is the vertical wave number of the receiving element:

[0025]

[0026] Where c is the speed of sound; ω is the frequency; k rx The horizontal wave number of the receiving element;

[0027] The cross-sectional view of the receiving oscillator and the transmitting oscillator after the above offset is as follows: Figure 3 As shown in (a) and (b).

[0028] Then the excitation element is shifted to depth z:

[0029]

[0030] Among them, k sz is the vertical wave number of the excited oscillator:

[0031]

[0032] Where c is the speed of sound; ω is the frequency; k sxThe horizontal wave number of the exciting element;

[0033] At this time, the cross-section of the receiving oscillator and the exciting oscillator is as follows Figure 3 As shown in (c) and (d), both are offset to depth z.

[0034] In general, the wave field at depth z can be directly obtained by the following formula:

[0035]

[0036] Among them, k z is the total vertical wave number of the wave field extrapolated downward, defined as:

[0037]

[0038] As a further preferred embodiment, if the sound velocity of the previous layer does not change in the horizontal direction, in the frequency wave number domain, the wave field P(k rx , z n-1 , k sx , z n-1 ,ω) to extrapolate and obtain the wave field P(k rx , z n , k sx , z n , ω) is calculated as:

[0039]

[0040] Among them, z n represents the depth of the nth layer, n∈[2,N], N is the number of discrete layers; Δz represents the discrete step length in the depth direction; k rx represents the horizontal wave number of the receiving oscillator; k sx represents the horizontal wave number of the transmitting oscillator; ω represents the frequency; i represents the imaginary unit; Indicates the depth is z n-1 The total vertical wave number at ; From formula 8), we get:

[0041]

[0042] Among them, c n-1 represents the speed of sound, which is a known parameter; Indicates that the receiving element is at a depth of z n-1 The vertical wave number at ; Indicates that the excitation element is at a depth of z n-1 The vertical wave number at .

[0043] In the above calculation process, when n=2, the previous layer is the first layer in the discrete layer of the imaging area, and its wave field P(krx , z1, k sx , z1, ω) is the surface wave field P(k rx , 0, k sx , 0, ω).

[0044] In the above calculation of the vertical phase offset used in the wave field of the current layer in the frequency wave number domain and the compensation operator used in the calculation of the wave field of the current layer in the frequency space domain in the inhomogeneous medium, the frequency band range [ω min ,ω max ] to traverse each discrete frequency in the range, and obtain the sound pressure field (wave field) corresponding to all discrete frequencies. For the convenience of calculation, the phase shift operator and compensation operator corresponding to all frequencies are integrated respectively, and the sound pressure field (wave field) corresponding to all discrete frequencies is directly calculated by point multiplication. Among them, the frequency band range [ω min ,ω max ] is determined by the full matrix data, and the frequency domain discrete step length dω is determined by the sampling frequency.

[0045] Preferably, in step (3), if the sound velocity of the previous layer does not change in the horizontal direction, the imaging condition for imaging using the wave field of the current layer in the frequency wavenumber domain is:

[0046]

[0047] Among them, I(x, z n ) represents the imaging result of the current layer; k rx represents the horizontal wave number of the receiving oscillator; k sx Represents the horizontal wave number of the transmitting oscillator; ω represents the frequency; i represents the imaginary unit; x represents the coordinate when the receiving oscillator and the exciting oscillator are in the same position; n∈[1,N], N represents the number of discrete layers.

[0048] Preferably, in step (3), when the sound velocity of the previous layer changes in the horizontal direction, the surface wave field is used as the wave field of the first layer, and all discrete layers are traversed in order from top to bottom. In the frequency wave number domain, the wave field of the previous layer is extrapolated, and the extrapolated wave field is transformed into the frequency space domain, and then a compensation operator is applied to obtain the wave field of the current layer in the frequency space domain; the imaging condition is applied to the wave field of the current layer in the frequency space domain and Fourier transform is performed to obtain the imaging result of the current layer and its wave field in the frequency wave number domain;

[0049] According to the above steps, the imaging result of the next layer and its wave field in the frequency wave number domain are obtained from the wave field of the current layer in the frequency wave number domain, and so on, the imaging results of all layers are obtained.

[0050] In the frequency wave number domain, the wave field of the current layer is obtained by extrapolating the wave field of the previous layer in the above manner, that is, the wave field of the first layer is first extrapolated to obtain the wave field of the second layer; then the wave field of the second layer is extrapolated to obtain the wave field of the third layer; and so on, the wave fields of all layers are obtained. If the sound velocity of the previous layer does not change in the horizontal direction (uniform layer), then the wave field of the previous layer in the frequency wave number domain is used for extrapolation processing and imaging according to the situation that there is no change.

[0051] As a further preferred embodiment, if the sound velocity of the upper layer varies in the horizontal direction, in the frequency wave number domain, the wave field P(k rx , z n-1 , k sx , z n-1 ,ω) for extrapolation (acting as a phase shift operator) is calculated as:

[0052]

[0053] Among them, z n represents the depth of the nth layer, n∈[2,N], N is the number of discrete layers; Δz is the discrete step length in the depth direction; k rx represents the horizontal wave number of the receiving oscillator; k sx represents the horizontal wave number of the transmitting oscillator; ω represents the frequency; i represents the imaginary unit; Indicates the depth is z n-1 The total number of vertical waves;

[0054]

[0055] At this time, c n-1 =1 / s0(z n-1 ), s0(z n-1 ) indicates the depth is z n-1 The average value of all slownesses in the horizontal plane at ; Indicates that the receiving element is at a depth of z n-1 The vertical wave number at ; Indicates that the excitation element is at a depth of z n-1 The vertical wave number at .

[0056] In the above calculation process, when n=2, the previous layer is the first layer in the discrete layer of the imaging area, and its wave field P(k rx , z1, k sx , z1, ω) is the surface wave field P(k rx , 0, k sx , 0, ω).

[0057] The derivation process of the wave field extrapolation calculation formula 12) for the inhomogeneous medium layer is the same as that of the wave field extrapolation calculation formula 9) for the homogeneous medium, but In the expression of n-1 ) is used as a substitute; and the uniform dielectric layer corresponds to In the expression of, the speed of sound c is a known parameter.

[0058] Preferably, in step (3), when the sound velocity of the previous layer varies in the horizontal direction (non-uniform layer), the compensation operator Φ′(x r , x s , z n , Δz, ω) is obtained by the following formula:

[0059] Φ′(x r , x s , z n , Δz,ω)=Φ(x r , z n ,Δz,ω)Φ(x s , z n , Δz, ω)

[0060] Among them, Φ(x r , z n ,Δz,ω),Φ(x s , z n , Δz, ω) are each independently selected from any one of the following three forms:

[0061]

[0062]

[0063]

[0064] Where ω represents the frequency; z n represents the depth of the nth layer, n∈[1, N-1], N is the number of discrete layers; i represents the imaginary unit; Δz represents the discrete step length in the depth direction; x represents the position of the oscillator; Δs(x, z n ) indicates the depth is z n Slowness changes in the horizontal plane.

[0065] The theoretical derivation process of wave field extrapolation in the above inhomogeneous layer is as follows:

[0066] 1. Extrapolation of the individual wave fields of receiving elements in inhomogeneous media

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

[0068]

[0069] Among them, p is the sound pressure; s is the medium slowness, that is, the reciprocal of the sound speed c; x and z are the two-dimensional space coordinates respectively; t is time.

[0070] After Fourier transform at t, the wave equation in the frequency domain is obtained as follows:

[0071]

[0072] 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 shown below

[0073]

[0074] in, is an imaginary unit;

[0075] The slowness s is an equation in space coordinates and can be decomposed into two components:

[0076] s(x,z)=s0(z)+Δs(x,z) 17)

[0077] Where s0(z) is the reference slowness, which is the average of all slownesses on the water surface at depth z, and Δs(x, z) includes all slowness changes on the horizontal plane at depth z.

[0078] For inhomogeneous media where the sound velocity varies both horizontally and in depth, the imaging area is first discretized and converted into multiple discrete layers. The total number of discrete layers is recorded as N.

[0079] Through Taylor expansion, the square root in formula 16) can be approximated as:

[0080]

[0081] If s0(z)>>Δs(x, z), Equation 16) can be solved by downward extrapolation to convert the wave field from P(x, z) to n-1 ,ω) to P(x,z n ,ω).

[0082] First, the wave field P(x, z n-1 ,ω) is transformed into the frequency wavenumber domain, P(k x , z n-1 ,ω), and then apply a phase shift operator:

[0083]

[0084] Where ω represents the frequency; z nrepresents the depth of the nth layer; i represents the imaginary unit; Δz represents the discrete step length in the depth direction; k x Indicates the wave number in the horizontal direction;

[0085]

[0086] Then the wave field P′(k x , z n ,ω) is converted to the frequency space domain, and P′(x, z n , ω), and through a disturbance term Φ(x, z n-1 , Δz, ω) to compensate for the impact of the change in horizontal sound speed, as shown in the following formula:

[0087] P(x,z n ,ω)=P′(x,z n ,ω)Φ(x,z n-1 , Δz, ω) 21)

[0088] Among them, Φ(x, z n-1 , Δz, ω) is the compensation operator, which has three forms:

[0089] The first form is:

[0090]

[0091] The second form is:

[0092]

[0093] The third form is:

[0094]

[0095] Where ω represents the frequency; z n represents the depth of the nth layer, n∈[2,N], N is the number of discrete layers; i represents the imaginary unit; Δz represents the discrete step length in the depth direction; x represents the position of the oscillator; Δs(x,z n ) indicates the depth is z n Slowness changes in the horizontal plane.

[0096] 2. Global extrapolation of wave fields in inhomogeneous media

[0097] For inhomogeneous media, the speed of sound is a function of x and z, so it is not feasible to perform Fourier transform on the wave field alone. Therefore, the inverse slowness s(x, z) of the speed of sound c is introduced, which is an equation of spatial coordinates and can be decomposed into two components:

[0098] s(x,z)=s0(z)+Δs(x,z) 25)

[0099] Where s0(z) is the reference slowness, which is the average of all slownesses on the water surface at depth z, and Δs(x, z) includes all slowness changes on the horizontal surface at depth z. When s0(z)>>Δs(x, z), then c=1 / s0(z).

[0100] Therefore, for the sound speed c in formula 13) n-1 Then 1 / s0(z n-1 )express.

[0101] Since the sound velocity of inhomogeneous media varies in the horizontal direction, imaging cannot be performed directly from the wave field in the frequency wavenumber domain obtained by extrapolation. It is necessary to transform it into the frequency space domain and then apply a horizontal compensation operator to compensate for the impact of the change in the sound velocity in the horizontal direction and improve the imaging quality.

[0102] Regarding disturbance compensation (compensation operator), the offset in the formula for the downward extrapolation of the wave field in an inhomogeneous medium can be regarded as firstly extrapolating the receiving element downward and then offsetting the exciting element downward, so a corresponding phase correction term is required in each process.

[0103] Therefore, we first perform an inverse Fourier transform to the frequency space domain, and then extrapolate the wave field of the previous layer downward in the frequency wave number domain in step (3) (the same as the derivation process of the uniform medium layer) to obtain the wave field P′(k rx , z n , k sx , z n ,ω), first perform inverse Fourier transform to the frequency space domain to obtain P′(x r , z n , x s , z n ,ω):

[0104]

[0105] Then for P′(x r , z n , x s , z n ,ω) acts as a compensation operator to obtain the wave field P(x) of the current layer in the frequency space domain r , z n , x s , z n ,ω):

[0106] P(x r , z n , x s , z n ,ω)=P′(x r , z n , x s , zn ,ω)Φ′(x r , x s , z n-1 , Δz, ω) 27)

[0107] Φ′(x r , x s , z n-1 , Δz,ω)=Φ(x r , z n-1 ,Δz,ω)Φ(x s , z n-1 , Δz, ω) 28)

[0108] Among them, Φ(x r , z n-1 ,Δz,ω),Φ(x s , z n-1 , Δz, ω) are each independently selected from formula 22), 23) or 24).

[0109] Finally, the imaging conditions are used to obtain P(x r , z n , x s , z n , ω) images the current layer to obtain the imaging result of the current layer.

[0110] For the wave field of the current layer (nth layer) in the frequency wave number domain, the wave field P(x r , z n , x s , z n ,ω) can be transformed into the frequency wavenumber domain.

[0111] In the above calculation of the vertical phase offset used in the wave field of the current layer in the frequency wave number domain and the compensation operator used in the calculation of the wave field of the current layer in the frequency space domain in the inhomogeneous medium, the frequency band range [ω min ,ω max ] to traverse each discrete frequency in the range, and obtain the sound pressure field (wave field) corresponding to all discrete frequencies. For the convenience of calculation, the phase shift operator and compensation operator corresponding to all frequencies are integrated respectively, and the sound pressure field (wave field) corresponding to all discrete frequencies is directly calculated by point multiplication. Among them, the frequency band range [ω min ,ω max ] is determined by the full matrix data, and the frequency domain discrete step length dω is determined by the sampling frequency.

[0112] After obtaining the wave field of the current layer in the frequency space domain, it is necessary to use imaging conditions to obtain the defect image. The explosion reflection model is an effective imaging condition. It regards the receiving wave as the result of active excitation of the sound source at the defect location. It is excited at time 0, so the defect image is the wave field at time 0. However, there are four spatial dimensions in the wave field of the full matrix data. In this case, the positions of the receiving oscillator and the transmitting oscillator should be the same to achieve the coincidence of the upgoing wave and the downgoing wave, that is, x r =x s =x.

[0113] Preferably, in step (4), when the sound velocity of the previous layer changes in the horizontal direction, the imaging condition for imaging using the wave field of the current layer in the frequency space domain is:

[0114] I(x,z n )=∫P(x r =x,z n , x s =x,z n ,ω)dω 28)

[0115] Among them, I(x, z n ) represents the imaging result of the current layer; x r Indicates the position of the receiving oscillator; x s represents the position of the exciting element; ω represents the frequency; z n represents the depth of the nth layer; P(x r =x,z n , x s =x,z n ,ω) represents x r =x s =x; n∈[1,N], N represents the number of discrete layers.

[0116] The invention discloses a method for high-efficiency full-matrix imaging of curved surface parts using ultrasonic phased arrays. The imaging method is implemented in the frequency domain and the frequency range [ω min ,ω max ] to improve the signal-to-noise ratio and reduce the amount of calculation. The vertical discrete resolution Δz also needs to be selected according to the actual situation, balancing accuracy and efficiency. Δz has a greater impact on the extrapolation accuracy of inhomogeneous media, and Δz can be reduced in inhomogeneous media. First, the imaging area is discretized to obtain N discrete layers; the full matrix data D(x r , x s, t) After performing a three-dimensional fast Fourier transform FFT, the surface wave field in the frequency wave number domain is obtained, and then the obtained surface wave field is used to reconstruct the wave field at the depth of each discrete layer. First, the surface wave field is extrapolated downward in the frequency wave number domain to obtain the wave field of the second layer; then the wave field of the second layer is extrapolated to obtain the wave field of the third layer, and so on, to obtain the wave fields of all layers. Determine the depth z of the previous layer n-1 n-1 The sound velocity at the location does not change in the horizontal direction. If the sound velocity is uniform in the horizontal direction, the wave field of the current layer n in the frequency wavenumber domain is directly used for imaging to obtain the imaging result of the current layer. n-1 If the sound velocity at the location changes in the horizontal direction, the wave field of the current layer nth layer in the frequency wave number domain is first transformed to the frequency space domain and then a compensation operator in the horizontal direction is applied to obtain the wave field of the current layer in the frequency space domain. Then the wave field of the current layer in the frequency space domain is imaged to obtain the imaging result of the current layer. Finally, the imaging results of all layers are integrated to obtain the overall imaging result of the inhomogeneous medium.

[0117] During the measurement of curved parts, the imaging area may include non-uniform medium structures and uniform medium structures. Therefore, after the imaging area is discretized, there may be non-uniform medium layers and uniform medium layers. For this reason, the ultrasonic phased array full-matrix efficient imaging method for curved parts of the present invention first converts the full-matrix data into a surface wave field in the frequency wave number domain, and then performs classification, extrapolation and imaging according to whether the sound velocity changes in the horizontal direction of each layer. For the upper layer where the sound velocity does not change in the horizontal direction, the wave field of the current layer in the frequency wave number domain obtained by extrapolation is directly used for imaging; for the upper layer where the sound velocity changes in the horizontal direction, it is necessary to convert the wave field of the current layer in the frequency wave number domain obtained by extrapolation to the frequency space domain and then apply a compensation operator in the horizontal direction to obtain the wave field of the current layer in the frequency space domain, and then use the wave field of the current layer in the frequency space domain for imaging. The above-mentioned classification imaging method can not only ensure the quality of the imaging results of the non-uniform medium layer, but also improve the overall imaging efficiency. In addition, there are only two operators, superposition and Fourier transform, in the imaging conditions of both cases, which greatly reduces the memory requirements and computational complexity and shortens the calculation time.

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

[0119] The ultrasonic phased array full-matrix efficient imaging method for curved surface parts of the present invention performs classification, extrapolation and imaging according to whether the sound velocity changes in the horizontal direction of each discrete layer, which can not only ensure the quality of the imaging results of the non-uniform medium layer, but also improve the overall imaging efficiency. In addition, there are only two operators, superposition and Fourier transform, in the imaging conditions of both cases, which greatly reduces the memory requirements and calculation complexity, and shortens the calculation time. It is particularly suitable for imaging applications of parts with both non-uniform media and uniform media. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0121] Figure 2 Schematic diagram of wave field propagation in inhomogeneous media;

[0122] Figure 3 Middle: (a) is a cross-sectional view of the full matrix with a common excitation oscillator after the receiving oscillator is shifted downward to z; (b) is a cross-sectional view of the full matrix with a common receiving oscillator after the receiving oscillator is shifted downward to z; (c) is a cross-sectional view of the full matrix with a common excitation oscillator after both the receiving oscillator and the excitation oscillator are shifted to z; (d) is a cross-sectional view of the full matrix with a common receiving oscillator after both the receiving oscillator and the excitation oscillator are shifted to z;

[0123] Figure 4 In the figure: (a) measurement schematic diagram of the annular component; (b) structural schematic diagram of the measured part (annular component); (c) measurement results based on the TFM method; (d) measurement results based on the PSM method; (e) measurement results based on the FD-RTM method; (f) measurement results of the method of the present invention. DETAILED DESCRIPTION

[0124] A method for high-efficiency full-matrix ultrasonic phased array imaging of curved parts, comprising the following steps:

[0125] (1) The imaging area of ​​the curved part is discretized to form multiple discrete layers in the depth direction; the number of discrete layers is set to N, and the step length in the depth z direction is Δz, where the size of Δz is determined according to the resolution requirement in the vertical direction.

[0126] (2) Convert the full matrix data of the curved surface part into the frequency wavenumber domain through Fourier transform to obtain the surface wave field;

[0127] Full Matrix Capture (FMC) collects the time domain signals between all excitation-receive pairs in the phased array probe. For an n-element phased array probe, the FMC data contains n 2 A scan signal. D(x r , x s, t) represents the signal transmitted from the sth oscillator and received by the rth oscillator at time t, such as Figure 2 As shown, the ultrasonic wave moves from the oscillator position (x s , 0), and after being reflected by the scattering point at (x, z), it is sent out by the r , received by the sensor at 0).

[0128] The traditional measurement method is to scan the surface with a probe to obtain measurement data. Here, D(x, t) represents the one-dimensional scan data, which can be directly regarded as the wave field on the surface P(x, z = 0, t). Compared with the scan data, FMC data has an extra dimension. Since both the excitation element and the receiving element have two-dimensional coordinates in space, in this case the wave field can be regarded as five-dimensional and expressed as P(x r , z r , x s , z s , t) represents, where x r and z r are the horizontal and vertical coordinates of the receiving element; x s and z s are the horizontal and vertical coordinates of the exciting element. Therefore, the wave field on the measurement surface (surface wave field) is P(x r ,0,x s , 0, t), and:

[0129] P(x r ,0,x s ,0,t)=D(x r , x s , t) 1)

[0130] P(x r ,0,x s , 0, t) can be obtained by a three-dimensional Fourier transform to obtain the surface wave field in the frequency wavenumber domain:

[0131]

[0132] Among them, k sx and k rx are the horizontal wave numbers of the exciting oscillator and the receiving oscillator respectively; x s and x r are the positions of the exciting oscillator and the receiving oscillator respectively; i is the imaginary unit; ω is the frequency; t is the time.

[0133] (3) Select any untraversed discrete layer as the current layer in order from top to bottom, and determine whether the sound speed of the previous layer changes in the horizontal direction:

[0134] 3.1) The speed of sound in the upper layer does not change in the horizontal direction

[0135] Taking the surface wave field as the first layer of wave field, all discrete layers are traversed from top to bottom. In the frequency wavenumber domain, the wave field of the previous layer is extrapolated to obtain the wave field of the current layer, and the imaging conditions are applied to obtain the imaging result of the current layer.

[0136] Specifically, in the frequency wavenumber domain, the wavefield of the current layer is obtained by extrapolating the wavefield of the previous layer, that is, the wavefield of the first layer is first extrapolated to obtain the wavefield of the second layer; then the wavefield of the second layer is extrapolated to obtain the wavefield of the third layer; and so on, the wavefields of all layers are obtained. Among them, the wavefield of the first layer is the surface wavefield. This layer-by-layer extrapolation method can obtain a more accurate wavefield of each layer, thereby improving the final imaging quality. If the sound velocity of the previous layer changes in the horizontal direction (non-uniform layer), the wavefield of the previous layer in the frequency wavenumber domain is used for extrapolation processing and imaging according to the change.

[0137] The specific derivation process is:

[0138] When the sound velocity c in the horizontal direction is uniform, the wave field in the measurement area can be directly obtained by downward extrapolation of the surface wave field in the frequency wavenumber domain. In order to obtain the wave field at depth z, the receiving element is first offset:

[0139]

[0140] Among them, k rz is the vertical wave number of the receiving element:

[0141]

[0142] Where c is the speed of sound; ω is the frequency; k rx The horizontal wave number of the receiving element;

[0143] The cross-sectional view of the receiving oscillator and the transmitting oscillator after the above offset is as follows: Figure 3 As shown in (a) and (b).

[0144] Then the excitation element is shifted to depth z:

[0145]

[0146] Among them, k sz is the vertical wave number of the excited oscillator:

[0147]

[0148] Where c is the speed of sound; ω is the frequency; k sx The horizontal wave number of the exciting element;

[0149] At this time, the cross-section of the receiving oscillator and the exciting oscillator is as follows Figure 3 As shown in (c) and (d), both are offset to depth z.

[0150] In general, the wave field at depth z can be directly obtained by the following formula:

[0151]

[0152] Among them, k z is the total vertical wave number of the wave field extrapolated downward, defined as:

[0153]

[0154] Then in the frequency wavenumber domain, for the wave field P(k rx , z n-1 , k sx , z n-1 ,ω) to extrapolate and obtain the wave field P(k rx , z n , k sx , z n , ω) is calculated as:

[0155]

[0156] Among them, z n represents the depth of the nth layer, n∈[2,N], N is the number of discrete layers; Δz represents the discrete step length in the depth direction; k rx represents the horizontal wave number of the receiving oscillator; k sx represents the horizontal wave number of the transmitting oscillator; ω represents the frequency; i represents the imaginary unit; Indicates the depth is z n-1 The total number of vertical waves; From formula 8), we get:

[0157]

[0158] Among them, c n-1 represents the speed of sound, which is a known value in the full matrix data; Indicates that the receiving element is at a depth of z n-1 The vertical wave number at ; Indicates that the excitation element is at a depth of z n-1 The vertical wave number at .

[0159] In the above calculation process, when n=2, the previous layer is the first layer in the discrete layer of the imaging area, and its wave field P(k rx , z1, k sx , z1, ω) is the surface wave field P(k rx , 0, k sx, 0, ω).

[0160] In the above calculation of the vertical phase offset used in the wave field of the current layer in the frequency wave number domain and the compensation operator used in the calculation of the wave field of the current layer in the frequency space domain in the inhomogeneous medium, the frequency band range [ω min ,ω max ] to traverse each discrete frequency in the range, and obtain the sound pressure field (wave field) corresponding to all discrete frequencies. For the convenience of calculation, the phase shift operator and compensation operator corresponding to all frequencies are integrated respectively, and the sound pressure field (wave field) corresponding to all discrete frequencies is directly calculated by point multiplication. Among them, the frequency band range [ω min ,ω max ] is determined by the full matrix data, and the frequency domain discrete step length dω is determined by the sampling frequency.

[0161] The imaging condition for imaging using the wave field of the current layer obtained by extrapolation in the frequency wavenumber domain is:

[0162]

[0163] Among them, I(x, z n ) represents the imaging result of the current layer; k rx represents the horizontal wave number of the receiving oscillator; k sx Represents the horizontal wave number of the transmitting oscillator; ω represents the frequency; i represents the imaginary unit; x represents the coordinate when the receiving oscillator and the exciting oscillator are in the same position; n∈[1,N], N represents the number of discrete layers.

[0164] 3.2) The speed of sound in the upper layer varies in the horizontal direction

[0165] Taking the surface wave field as the wave field of the first layer, traverse all discrete layers in order from top to bottom. In the frequency wave number domain, the wave field of the previous layer is extrapolated, and the extrapolated wave field is transformed into the frequency space domain. Then a compensation operator is applied to obtain the wave field of the current layer in the frequency space domain. Imaging conditions and Fourier transform are applied to the wave field of the current layer in the frequency space domain to obtain the imaging result of the current layer and its wave field in the frequency wave number domain.

[0166] According to the above steps, the imaging result of the next layer and its wave field in the frequency wave number domain are obtained from the wave field of the current layer in the frequency wave number domain, and so on, the imaging results of all layers are obtained.

[0167] Specifically, in the frequency wavenumber domain, the wavefield of the current layer is obtained by extrapolating the wavefield of the previous layer, that is, the wavefield of the first layer is first extrapolated to obtain the wavefield of the second layer; then the wavefield of the second layer is extrapolated to obtain the wavefield of the third layer; and so on, the wavefields of all layers are obtained. If the sound velocity of the previous layer does not change in the horizontal direction (uniform layer), then the wavefield of the previous layer in the frequency wavenumber domain is used for extrapolation processing and imaging according to the situation that there is no change.

[0168] The specific derivation process is:

[0169] In the frequency wavenumber domain, for the wave field P(k rx , z n-1 , k sx , z n-1 ,ω) for extrapolation (acting as a phase shift operator) is calculated as:

[0170]

[0171] Among them, z n represents the depth of the nth layer, n∈[2,N], N is the number of discrete layers; Δz is the discrete step length in the depth direction; k rx represents the horizontal wave number of the receiving oscillator; k sx represents the horizontal wave number of the transmitting oscillator; ω represents the frequency; i represents the imaginary unit; Indicates the depth is z n-1 The total number of vertical waves;

[0172]

[0173] At this time, c n-1 =1 / s0(z n-1 ), s0(z n-1 ) indicates the depth is z n-1 The average value of all slownesses in the horizontal plane at ; Indicates that the receiving element is at a depth of z n-1 The vertical wave number at ; Indicates that the excitation element is at a depth of z n-1 The vertical wave number at .

[0174] In the above calculation process, when n=2, the previous layer is the first layer in the discrete layer of the imaging area, and its wave field P(k rx , z1, k sx , z1, ω) is the surface wave field P(k rx , 0, k sx , 0, ω).

[0175] The derivation process of the wave field extrapolation calculation formula 12) for the inhomogeneous medium layer is the same as that of the wave field extrapolation calculation formula 9) for the homogeneous medium, but The expression is, the speed of sound c n-1 1 / s0(z n-1 ) is used as a substitute; and the uniform dielectric layer corresponds to The expression is, the speed of sound c n-1 are known parameters.

[0176] The theoretical derivation process of wave field extrapolation in the above inhomogeneous layer is as follows:

[0177] 1. Extrapolation of the individual wave fields of receiving elements in inhomogeneous media

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

[0179]

[0180] Among them, p is the sound pressure; s is the medium slowness, that is, the reciprocal of the sound speed c; x and z are the two-dimensional space coordinates respectively; t is time.

[0181] After Fourier transform at t, the wave equation in the frequency domain is obtained as follows:

[0182]

[0183] 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 shown below

[0184]

[0185] in, is an imaginary unit;

[0186] The slowness s is an equation in space coordinates and can be decomposed into two components:

[0187] s(x,z)=s0(z)+Δs(x,z) 17)

[0188] Where s0(z) is the reference slowness, which is the average of all slownesses on the water surface at depth z, and Δs(x, z) includes all slowness changes on the horizontal plane at depth z.

[0189] For inhomogeneous media where the sound velocity varies both horizontally and in depth, the imaging area is first discretized and converted into multiple discrete layers. The total number of discrete layers is recorded as N.

[0190] Through Taylor expansion, the square root in formula 16) can be approximated as:

[0191]

[0192] If s0(z)>>Δs(x, z), Equation 16) can be solved by downward extrapolation to convert the wave field from P(x, z) to n-1 ,ω) to P(x,z n ,ω).

[0193] First, the wave field P(x, z n-1 ,ω) is transformed into the frequency wavenumber domain, and P(k x , z n-1 ,ω), and then apply a phase shift operator:

[0194]

[0195] Where ω represents the frequency; z n represents the depth of the nth layer; i represents the imaginary unit; Δz represents the discrete step length in the depth direction; k x Indicates the wave number in the horizontal direction;

[0196]

[0197] Then the wave field P′(k x , z n ,ω) is converted to the frequency space domain, and P′(x, z n , ω), and through a disturbance term (compensation operator) Φ(x, z n-1 , Δz, ω) to compensate for the impact of the change in horizontal sound speed, as shown in the following formula:

[0198] P(x,z n ,ω)=P′(x,z n ,ω)Φ(x,z n-1 , Δz, ω) 21)

[0199] Among them, Φ(x, z n-1 , Δz, ω) is the compensation term, which has three forms:

[0200] The first form is:

[0201]

[0202] The second form is:

[0203]

[0204] The third form is:

[0205]

[0206] Where ω represents the frequency; z nrepresents the depth of the nth layer; i represents the imaginary unit; Δz represents the discrete step length in the depth direction; x represents the position of the oscillator; Δs(x, z n ) indicates the depth is z n Slowness changes in the horizontal plane.

[0207] 2. Global extrapolation of wave fields in inhomogeneous media

[0208] For inhomogeneous media, the speed of sound is a function of x and z, so it is not feasible to perform Fourier transform on the wave field alone. Therefore, the inverse slowness s(x, z) of the speed of sound c is introduced, which is an equation of spatial coordinates and can be decomposed into two components:

[0209] s(x,z)=s0(z)+Δs(x,z) 25)

[0210] Where s0(z) is the reference slowness, which is the average of all slownesses on the water surface at depth z, and Δs(x, z) includes all slowness changes on the horizontal surface at depth z. When s0(z)>>Δs(x, z), then c=1 / s0(z).

[0211] Therefore, for the sound speed c in formula 13) n-1 Then 1 / s0(z n-1 )express.

[0212] Since the sound velocity of inhomogeneous media varies in the horizontal direction, imaging cannot be performed directly from the wave field in the frequency wavenumber domain obtained by extrapolation. It is necessary to transform it into the frequency space domain and then apply a horizontal compensation operator to compensate for the impact of the change in the sound velocity in the horizontal direction and improve the imaging quality.

[0213] Regarding disturbance compensation (compensation operator), the offset in the formula for the downward extrapolation of the wave field in an inhomogeneous medium can be regarded as firstly extrapolating the receiving element downward and then offsetting the exciting element downward, so a corresponding phase correction term is required in each process.

[0214] Therefore, we first perform an inverse Fourier transform to the frequency space domain, and then extrapolate the wave field of the previous layer downward in the frequency wave number domain in step (3) to obtain the wave field P′(k rx , z n , k sx , z n ,ω), first perform inverse Fourier transform to the frequency space domain to obtain P′(x r , z n , x s , z n ,ω):

[0215]

[0216] Then for P′(x r , z n , x s , z n ,ω) acts as a compensation operator to obtain the wave field P(x) of the current layer in the frequency space domain r , z n , x s , z n ,ω):

[0217] P(x r , z n , x s , z n ,ω)=P′(x r , z n , x s , z n ,ω)Φ′(x r , x s , z n-1 , Δz, ω) 27)

[0218] Among them, Φ′(x r , x s , z n-1 , Δz, ω) is as follows:

[0219] Φ′(x r , x s , z n-1 , Δz,ω)=Φ(x r , z n-1 ,Δz,ω)Φ(x s , z n-1 , Δz, ω) 28)

[0220] Φ(x r , z n-1 ,Δz,ω),Φ(x s , z n-1 , Δz, ω) are each independently selected from formula 22), 23) or 24).

[0221] In the above calculation of the vertical phase offset used in the wave field of the current layer in the frequency wave number domain and the compensation operator used in the calculation of the wave field of the current layer in the frequency space domain in the inhomogeneous medium, the frequency band range [ω min ,ω max ] to traverse each discrete frequency in the range, and obtain the sound pressure field (wave field) corresponding to all discrete frequencies. For the convenience of calculation, the phase shift operator and compensation operator corresponding to all frequencies are integrated respectively, and the sound pressure field (wave field) corresponding to all discrete frequencies is directly calculated by point multiplication. Among them, the frequency band range [ω min ,ω max] is determined by the full matrix data, and the frequency domain discrete step length dω is determined by the sampling frequency.

[0222] After obtaining the wave field of the current layer in the frequency space domain, it is necessary to use imaging conditions to obtain the defect image. The explosion reflection model is an effective imaging condition. It regards the receiving wave as the result of active excitation of the sound source at the defect location. It is excited at time 0, so the defect image is the wave field at time 0. However, there are four spatial dimensions in the wave field of the full matrix data. In this case, the positions of the receiving oscillator and the transmitting oscillator should be the same to achieve the coincidence of the upgoing wave and the downgoing wave, that is, x r =x s =x.

[0223] The imaging condition for imaging using the wave field of the current layer in the frequency space domain is:

[0224] I(x,z n )=∫P(x r =x,z n , x s =x,z n ,ω)dω 29)

[0225] Among them, I(x, z n ) represents the imaging result of the current layer; x r Indicates the position of the receiving oscillator; x s represents the position of the exciting element; ω represents the frequency; z n represents the depth of the nth layer; P(x r =x,z n , x s =x,z n ,ω) represents x r =x s =x; n∈[1,N], N represents the number of discrete layers.

[0226] For the wave field of the current layer (nth layer) in the frequency wave number domain, the wave field P(x r , z n , x s , z n ,ω) can be transformed into the frequency wavenumber domain.

[0227] (4) Repeat step (3) until all layers are traversed, obtain the imaging results of all layers, and integrate them to obtain the imaging results of the curved surface parts.

[0228] Detection experiment

[0229] This experiment is for a ring component such as Figure 4Measurements were performed in (b) to verify the performance of the imaging method of this embodiment in practice. Figure 4 As shown in (a), the element (ring component) is mounted on a triangular chuck, which drives the element to rotate, and a measurement is performed every 45° rotation. A 64-element phased array probe (Shantou Ultrasound) with a center frequency of 2.5MHz and an element spacing of 0.75mm was used in the experiment. The FMC data was obtained using a data acquisition card 64 / 64OEM-PA (AOS. Ltd, America), with a sampling frequency of 50MHz and a time range of 60μs. Therefore, a complete measurement includes 8 sets of FMC data. The material of the measurement component is PEEK (c = 2500m / s), and there are 8 through holes with a diameter of 3mm inside, and their spatial positions are as shown in the figure. Figure 4 The 8 sets of FMC data were processed by the TFM-based method, the PSM-based method, the frequency domain reverse time migration method (FD-RTM) and the imaging method of this embodiment (Proposed), and the imaging results of each method were obtained as shown in FIG. Figure 4 As shown in (c)-(f).

[0230] In order to quantitatively analyze the imaging results, a dimensionless parameter API (array performance indicator) is introduced as follows:

[0231] API = A -6dB / λ 2 30)

[0232] 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.

[0233] Table 1 lists the time consumption, average API and positioning error of various methods. r and E θ are radial and circumferential errors respectively.

[0234] Table 1. Computation time and average API of imaging results

[0235]

[0236] Depend on Figure 4As can be seen from (c) to (f) and Table 1, the TFM-based method requires two hours to reconstruct the entire imaging result, and both the TFM-based method and the PSM-based method have obvious background artifacts in their imaging results. The FD-RTM method is also time-consuming, and there are multiple echo artifacts in its results, while its radial positioning accuracy is the highest. The imaging method of this embodiment (Proposed) shows great advantages in both computational efficiency and imaging quality, with the highest resolution and lowest axial error, and the computational time is 1 / 2 of that of the PSM-based method, 1 / 54 of that of the FD-RTM method, and 1 / 44 of that of the TFM-based method.

Claims

1. A method for full-matrix ultrasonic phased array imaging of curved parts, characterized in that: The following steps are involved: (1) Discretizing the imaging area of ​​the curved part to form multiple discrete layers in the depth direction; (2) Convert the full matrix data of the curved surface part into the frequency wavenumber domain through Fourier transform to obtain the surface wave field; (3) Taking the surface wave field as the wave field of the first layer, select any untraversed discrete layer as the current layer in order from top to bottom, and determine whether the sound velocity of the previous layer changes in the horizontal direction: If yes, in the frequency-wavenumber domain, the wave field of the previous layer is extrapolated downward and transformed to the frequency-space domain, and then a compensation operator is applied to obtain the wave field of the current layer in the frequency-space domain; the imaging condition and Fourier transform are applied to the wave field of the current layer in the frequency-space domain to obtain the imaging result of the current layer and its wave field in the frequency-wavenumber domain; If not, then in the frequency-wavenumber domain, the wave field of the previous layer is extrapolated downward to obtain the wave field of the current layer, and the imaging conditions are applied to obtain the imaging result of the current layer; (4) Repeat step (3) until all layers are traversed, obtain the imaging results of all layers, and integrate them to obtain the imaging results of the curved surface part; In step (3), when the sound velocity of the previous layer changes in the horizontal direction, in the frequency wave number domain, the wave field P(k rx ,z n-1 ,k sx ,z n-1 ,ω) is calculated as follows: Among them, P ′ (k rx ,z n ,k sx ,z n ,ω) represents the wave field of the current layer obtained by extrapolating the wave field of the previous layer; z n represents the depth of the nth layer, n∈[2,N], N is the number of discrete layers; Δz represents the discrete step length in the depth direction; k rx Indicates the horizontal wave number of the receiving element; k sx represents the horizontal wave number of the transmitting oscillator; ω represents the frequency; i represents the imaginary unit; Indicates the depth is z n-1 The total vertical wave number at ; At this time, c n-1 =1 / s0(z n-1 ), s0(z n-1 ) indicates the depth is z n-1 The average value of all slownesses in the horizontal plane at ; Indicates that the receiving element is at a depth of z n-1 The vertical wave number at time ; Indicates that the excitation element is at a depth of z n-1 The vertical wave number at ; In step (3), the compensation operator Φ ′ (x r ,x s ,z n ,Δz,ω) is obtained by the following formula: F ′ (x r ,x s ,z n ,Δz,ω)=Φ(x r ,z n ,Δz,ω)Φ(x s ,z n ,Δz,ω) Among them, Φ(x r ,z n ,Δz,ω)、Φ(x s ,z n ,Δz,ω) are independently selected from any of the following three forms: Where ω represents the frequency; z n represents the depth of the nth layer, n∈[1, N-1], N is the number of discrete layers; i represents the imaginary unit; Δz represents the discrete step length in the depth direction; x represents the position of the oscillator; Δs(x,z n ) indicates the depth is z n Slowness changes in the horizontal plane.

2. The method for full-matrix ultrasonic phased array imaging of curved parts according to claim 1, characterized in that: In step (3), when the sound velocity of the previous layer changes in the horizontal direction, the imaging condition for imaging using the wave field of the current layer in the frequency space domain is: I(x,z n )=∫P(x r =x,z n ,x s =x,z n ,ω)dω Among them, I(x,z n ) represents the imaging result of the current layer; x r Indicates the position of the receiving oscillator; x s represents the position of the exciting oscillator; x represents the position of the oscillator; ω represents the frequency; z n represents the depth of the nth layer; P(x r =x,z n ,x s =x,z n ,ω) represents x r =x s = the wave field of the current layer at time x; n∈[1,N], where N represents the number of discrete layers.

3. The method for full-matrix ultrasonic phased array imaging of curved parts according to claim 1, characterized in that: In step (3), if the sound velocity of the previous layer does not change in the horizontal direction, in the frequency wavenumber domain, the wave field P(k rx ,z n-1 ,k sx ,z n-1 ,ω) to obtain the wave field P(k rx ,z n ,k sx ,z n ,ω) is calculated as: Among them, z n represents the depth of the nth layer, n∈[2,N], N is the number of discrete layers; Δz represents the discrete step length in the depth direction; k rx Indicates the horizontal wave number of the receiving element; k sx represents the horizontal wave number of the transmitting oscillator; ω represents the frequency; i represents the imaginary unit; Indicates the depth is z n-1 The total vertical wave number at ; Among them, c n-1 represents the speed of sound; Indicates that the receiving element is at a depth of z n-1 The vertical wave number at ; Indicates that the excitation element is at a depth of z n-1 The vertical wave number at .

4. The method for full-matrix ultrasonic phased array imaging of curved parts according to claim 1, characterized in that: In step (3), if the sound velocity of the previous layer does not change in the horizontal direction, the imaging condition for imaging using the wave field of the current layer in the frequency wavenumber domain is: Among them, I(x,z n ) represents the imaging result of the current layer; k rx Indicates the horizontal wave number of the receiving element; k sx Represents the horizontal wave number of the transmitting oscillator; ω represents the frequency; i represents the imaginary unit; x represents the coordinate when the receiving oscillator and the exciting oscillator are in the same position; n∈[1,N], N represents the number of discrete layers.

Images