A high-precision ultrasonic imaging method for internal defects of complex-shaped parts based on search vector imaging conditions

Abstract

The present invention provides a high-precision ultrasonic imaging method for internal defects of complex-shaped parts based on search vector imaging conditions. Starting from full waveform inversion, the search vector imaging conditions are constructed. Through the first-order derivative vector combined with the approximate Hessian matrix correction, high-quality images of defects in complex sound velocity models can be obtained, which is suitable for defect detection of complex surface components. The imaging conditions used in the imaging method of the present invention do not make assumptions about the wave field and are more rigorous than traditional imaging condition theory, which ensures the high precision of the inventive method. In addition, the imaging method of the present invention can highlight defects in the entire measurement area, which is very convenient for locating defects in large measurement structures. Compared with traditional methods, the imaging method of the present invention significantly improves imaging resolution and signal-to-noise ratio (SNR).

Description

Technical Field

[0001] The present invention belongs to the technical field of ultrasonic imaging, and in particular relates to a high-precision ultrasonic imaging method for internal defects of complex-shaped parts based on search vector imaging conditions. Background Art

[0002] With the continuous advancement of manufacturing technology, the constraints on component shapes are becoming less and less. As a result, components with complex shapes are mass-produced to meet the requirements of various working conditions. Shafts and camshafts are common transmission components with curved surfaces. At the same time, more complex surface components have emerged in the automotive, energy, and aerospace fields. However, these components are prone to holes and cracks during the tedious manufacturing process and severe service life, which greatly weakens the mechanical properties of the components and causes them to fail prematurely or even cause serious accidents. Therefore, high-quality detection of internal defects in complex-shaped components has become a research hotspot. However, the large inspection area and complex surface shape pose great challenges to detection.

[0003] Three typical nondestructive testing (NDT) technologies are widely used for imaging internal defects: infrared (IR) thermography, X-ray tomography, and ultrasonic scanning. In IR thermography, by heating the component, heat accumulates at the defect site, generating a localized temperature increase. This temperature rise is captured by an infrared camera and used to determine the defect's location. However, IR thermography can only detect internal defects close to the surface and cannot determine the defect's depth. X-ray tomography captures radiographs from different angles to reconstruct internal defects with high precision. However, this reconstruction process requires tedious and complex calculations, making it unsuitable for inspecting large components. Ultrasonic scanning has fewer limitations on components due to its strong penetrating power and low operating environment requirements. Furthermore, ultrasonic phased array imaging has gained prominence in NDT due to its portability, reliability, and relatively high accuracy, effectively detecting defects located within components. Full matrix capture (FMC) collects the complete time-domain signal between each pair of transmitting and receiving elements, improving imaging resolution and reducing distortion of the imaged target.

[0004] The surface of an ultrasonic phased array is typically flat, so the probe can be coupled to complex-shaped components through water or specialized wedges. Under these measurement conditions, the sound velocity within the measurement area varies both horizontally and vertically. Three imaging methods are used for this purpose: total focusing function (TFM), phase shift migration (PSM), and reverse time migration (RTM). These methods all employ excitation-time imaging to generate defect images. They consider the ultrasonic field as a combination of an incident wavefield and a reflected wavefield, where the reflected wavefield is actively emitted from the reflection point at a specific moment, and the excitation time is determined by the incident wavefield.

[0005] TFM-based methods use rays to approximate wavefield propagation, reconstruct the wavefield through time delays, and introduce ray tracing methods to calculate ray paths. However, the imaging accuracy of TFM-based methods depends on the accuracy of the path search, which is most suitable for two-layer models and easily fails when the sound velocity model is complex. PSM-based methods reconstruct the incident and reflected wavefields through phase extrapolation in the frequency-wavenumber domain, which is easily adaptable to variations in the sound velocity in the vertical direction. Researchers have introduced phase shift plus interpolation (PSPI) and non-stationary phase shift (NSPS) methods to compensate for variations in the sound velocity in the horizontal direction. However, the phase shift is a first-order approximation of the wave equation, and PSPI and NSPS assume that the sound velocity variations in the horizontal direction are small relative to the average sound velocity. Therefore, when the sound velocity varies sharply in the horizontal direction, the imaging accuracy of PSM-based methods will be significantly reduced. RTM-based methods use a finite difference solution of the wave equation to reconstruct the incident and reflected wavefields. Compared with the previous two methods, RTM-based methods make significantly fewer assumptions about wavefield propagation, thus accommodating any complex shape and improving imaging accuracy. However, the RTM-based method still suffers from the disadvantage of insufficient illumination at deep locations due to acoustic wave scattering, so it is usually impossible to highlight the defect image in the entire measurement area, and manual search is required in the local area to obtain the defect image. Summary of the Invention

[0006] To address the problems existing in the prior art, the present invention provides a high-precision ultrasonic imaging method for internal defects in complex-shaped parts based on search vector imaging conditions. This search vector (vector) imaging condition utilizes the first-order derivative vector of full waveform inversion and an approximate Hessian matrix to obtain highly accurate defect images, enabling high-precision detection of internal defects in complex-shaped parts. This imaging method can highlight defect images within the entire measurement area, making it very convenient for defect location in large measurement structures. The resulting image clearly highlights all pores and cracks, resulting in high imaging resolution and signal-to-noise ratio (SNR).

[0007] A high-precision ultrasonic imaging method for internal defects of complex-shaped parts based on search vector imaging conditions comprises the following steps:

[0008] (1) Input the excitation signal, sound velocity model and full matrix data of the measurement;

[0009] (2) Preprocess the excitation signal and the full matrix data respectively to obtain the sound source matrix and residual matrix in the frequency domain; discretize the sound velocity model to obtain the sound velocity vector, and calculate the impedance matrix;

[0010] (3) For any untraversed frequency, the LU decomposition of the impedance matrix is ​​used to solve the source matrix and the residual matrix respectively to obtain the forward wave field and residual wave field of the frequency;

[0011] (4) The impedance matrix is ​​partial-derived with respect to the acoustic velocity vector and then multiplied with the forward wave field to obtain the total virtual source matrix; at the same time, the forward wave field is normalized by the amplitude of the excitation value of the excitation signal to obtain the unit forward wave field;

[0012] (5) Multiply the total virtual source matrix with the unit forward wave field and the residual wave field respectively to obtain the partial derivative matrix and gradient matrix corresponding to the frequency;

[0013] (6) Repeat steps (3) to (5) until all frequencies are traversed, and obtain the partial derivative matrix and gradient matrix of all frequencies, and use the imaging conditions of the search vector to obtain the imaging results of the internal defects of the part.

[0014] In the above step (1), the excitation signal is S(x r , x s , t), the speed of sound model is represented by c(x, z), the full matrix data D(x r , x s , t) represents; where x r represents the rth receiving element, x s Represents the sth transmitting oscillator, r and s are independently taken from [1, M], M is the number of oscillators; t represents time.

[0015] Preferably, in step (2), the excitation signal and the full matrix data are preprocessed separately, that is, a temporal one-dimensional Fourier transform is performed on the excitation signal and the full matrix data respectively.

[0016] Preferably, in step (2), the calculation formula of the residual matrix δD(c) is as follows:

[0017]

[0018] Where D represents the measurement value of the full matrix data; is the estimated value of the full matrix data;

[0019] Assuming that the sound velocity model is smooth and no reflection occurs, the estimated value of the full matrix data is 0, and the residual matrix δD(c) = D;

[0020] The full matrix data in the time domain of the measurement is subjected to a one-dimensional Fourier transform in time to obtain the full matrix data in the frequency domain, that is, the residual matrix in the frequency domain is obtained, that is, δD(c, ω) = D(x r , x s ,ω), ω represents the frequency.

[0021] Preferably, in step (2), the sound velocity model is discretized using the finite difference method to obtain the sound velocity vector, and the impedance matrix is ​​calculated using the following formula:

[0022]

[0023] Where A(ω) is the N×N impedance matrix; represents the Laplace operator obtained after discretization; c is the N×1 sound speed vector.

[0024] Preferably, in step (3), the solution formula for the forward wave field is as follows:

[0025] P(ω)=[p1(ω), p2(ω),...,p M (ω)]=A(ω) -1 S(ω)

[0026] Where S(ω) represents the sound source matrix; P(ω) represents the forward wave field; A(ω) -1 Indicates the inversion of A(ω);

[0027] The solution formula for the residual wave field is as follows:

[0028] V(ω j )=A(ω j ) -1 W T δD(c,ω j )

[0029] Where, V(ω j ) represents the frequency ω j The corresponding residual wave field; δD(c, ω j ) represents the frequency ω j The corresponding residual matrix, j∈[1, N ω ],N ω represents the number of discrete frequencies; W is an M×N matrix, which extracts the sound field value at the node corresponding to the receiving array element in the entire wave field. In the kth row, only the element at the node corresponding to the kth receiving array element is 1, and the other elements are 0, where k∈[1,M] and M<N; the superscript T represents the transpose.

[0030] Preferably, in step (4), the total virtual source matrix is ​​defined as:

[0031] F(ω)=[F( 1 )(ω),...,F (N) (ω)]

[0032] in,

[0033]

[0034] Where, F (n) (ω) is regarded as the nth virtual source, which is an N×M matrix; c nRepresents the nth sound speed parameter (element) in the sound speed vector; n∈[1,N];

[0035] The forward wave field is normalized using the following formula:

[0036]

[0037] Where, is the unit forward wave field; a(ω j )=S(x s , x s ,ω j ), for S(x r , x s ,ω j ), only in x r =x s There is an incentive value of a(ω j ); when x r ≠x s When S(x r , x s ,ω j )=0.

[0038] Preferably, in step (5), the calculation formula for the partial derivative matrix obtained by multiplying the total virtual source matrix by the unit forward wave field is as follows:

[0039] Where, J l is the N×1 partial derivative matrix of the lth data, where l=[1,...,L], L=N ω ×M×M, represents the number of data points; l=(j-1)M 2 +(s-1)M+r;

[0040] f (l，n) is the virtual source vector corresponding to the lth data and the nth parameter, which is obtained from the total virtual source matrix, where: the total virtual source matrix F(ω)=[F (1) (ω),...,F (N) (ω)]

[0041] For frequency ω j exist:

[0042]

[0043] In the above formula, the right side of the equation replicates the virtual source matrix 64 times in the row direction;

[0044] By unit forward wave field get:

[0045]

[0046] The calculation formula of the gradient matrix obtained by multiplying the total virtual source matrix and the residual wave field is:

[0047]

[0048] In the formula, F(ω j ) represents the frequency ω j The corresponding total virtual source matrix.

[0049] Preferably, the imaging condition for searching the vector in step (6) is:

[0050]

[0051] Where Im represents the imaging result; I is an N-dimensional unit matrix; represents the N×1 gradient matrix; γ is the damping coefficient, which is used to stabilize the inversion;

[0052] H a =J T J * , where J represents the partial derivative matrix; superscript T and * represent transpose and complex conjugate respectively;

[0053] The nth element in is obtained by the following formula:

[0054]

[0055] Where, The subscript represents the element located in the (s-1)N+nth row and the rth column; n∈[1,N];

[0056] According to all partial derivative matrices, H is calculated by the following formula a The N diagonal elements of :

[0057]

[0058] Where, diag(H a ) represents the approximate Hessian matrix; denotes Hadamard multiplication, l = (j-1)M 2 +(s-1)M+r.

[0059] The theoretical derivation process of the imaging method of the present invention is as follows:

[0060] 1.1 Frequency Wave Field

[0061] The frequency domain wave field in inhomogeneous media satisfies the following equation:

[0062]

[0063] Among them, x and z are two spatial coordinates in two-dimensional space, ω is frequency, z is depth, p(x, z, ω) is the sound field (wave field), s(x, z, ω) is the sound source field, and c(x, z) is the sound speed model.

[0064] In engineering practice, formula (1) is usually solved by the finite difference method, where the wave field and sound velocity model are discretized as N = N x ×N z points, such as Figure 2 shown.

[0065] The discrete wave equation is as follows:

[0066]

[0067] Among them, c is the sound velocity vector, p(ω) is the wave field vector, and s(ω) is the sound source vector, and all three are N×1 vectors;

[0068] is an N×N impedance matrix, where represents the Laplace operator obtained after discretization.

[0069] Therefore, the wave field p(ω) can be obtained by the following formula:

[0070] p(ω)=A(ω) -1 s(ω) (3)

[0071] Among them, A(ω) -1 Indicates the inversion of A(ω); the calculation of the wave field can introduce the LU decomposition of the impedance matrix to improve the calculation efficiency.

[0072] The schematic diagram of phased array defect detection is as follows: Figure 3 As shown, each array element is excited in turn, and all array elements receive signals at the same time. Therefore, the FMC data (full matrix data) contains the time domain signal of each transmit-receive pair. Here, D(x r , x s , t) represents, where x r and x s They correspond to the rth receiving array element and the sth transmitting array element respectively; s∈[1,M], r∈[1,M], where M is the number of elements. Figure 2 As shown, the ultrasonic signal is transmitted from the array element (x s ,0) is emitted and reflected by the interface and defects and is located at (x r Each element in the phased array can be represented by a node corresponding to its center point in the discrete model, so the wave fields at these nodes can be known from the FMC dataset.

[0073] Using the sound source matrix S(ω)=[s1(ω),s2(ω),...,s M (ω)] represents the sound source of a phased array with M elements. The wave field corresponding to the excitation of these M elements can be solved at one time by the following formula:

[0074] P(ω)=[p1(ω), p2(ω),...,p M (ω)]=A(ω) -1 S(ω) (4)

[0075] Among them, P(ω) is the wave field matrix of the full matrix data.

[0076] 1.2 Full Waveform Inversion and Gradient Vectors

[0077] Given the sound source matrix and the sound velocity model, the wave field matrix of the FMC data can be calculated according to Equation (4), and the frequency domain FMC data can be obtained by taking the wave field values ​​at the corresponding nodes of the receiver, which is written as:

[0078]

[0079] in, Represents the estimated value of the full matrix data; W is an M×N matrix, which extracts the sound field value at the node corresponding to the receiving array element in the entire wave field. In the kth row, only the element at the node corresponding to the kth receiving array element is 1, and the other elements are 0, where k∈[1,M] and M<N; is an M×M matrix, and its element at the (r, s) position corresponds to the signal received by the rth element when the sth element is excited. By setting the Integration is obtained.

[0080] The calculated full matrix data The difference between the measured full matrix data D is defined as the residual matrix, written as:

[0081]

[0082] Among them, δD(c) and The whole is a function of the sound velocity vector c. The sound velocity vector c that minimizes the modulus of the residual matrix δD(c) is the closest to the true sound velocity, which is the basis of full waveform inversion.

[0083] Here, the l2 norm of the residual data δD(c) is used to define the objective function:

[0084]

[0085] Wherein, the superscripts T and * represent transpose and complex conjugate respectively; Nω is the number of discrete frequencies; M represents the number of oscillators; D(x r , x s ,ω) by D(x r , x s , t) is obtained by Fourier transform in the time dimension.

[0086] The negative gradient direction can be used to update the sound speed vector, written as:

[0087]

[0088] in, represents the gradient matrix; J is the L×N Frechet partial derivative matrix, and its elements are as follows:

[0089]

[0090] Where L = N ω ×M×M number of data points, l=(j-1)M 2 +(s-1)M+r;j∈[1,N ω ].

[0091] Apply both sides of formula (4) to the nth sound speed parameter c n Taking partial derivatives we can get formula (10):

[0092]

[0093] in, is an N×N partial derivative impedance matrix, which is the impedance matrix for the nth sound velocity parameter c n Find the result of partial derivative.

[0094] The form of formula (10) is similar to formula (4), F (n) (ω) can be regarded as the nth virtual source, which is an N×M matrix and written as:

[0095]

[0096] therefore, is an N×M matrix, which is the solution to a forward problem.

[0097] Then, at frequency ω j The Frechet partial derivative matrix when can be written as:

[0098]

[0099] Among them, J(ω j ) is an M×MN matrix, and the entire Frechet partial derivative matrix J can be obtained by integrating all discrete frequencies ωj The corresponding J(ω j )get.

[0100] The total virtual source matrix can be defined as F(ω)=[F (1) (ω),...,F (N) (ω)], which is an N×MN matrix.

[0101] Explicitly calculating the partial derivative matrix J is very time-consuming. Substituting formula (12) into formula (8), the gradient matrix can be directly calculated:

[0102]

[0103] Among them, δD(c,ω j ) is the frequency ω j The M×M residual data matrix when , The frequency is ω j The MN×M gradient matrix when .

[0104] Residual wave field V(ω j ) is an N×M matrix and is calculated as follows:

[0105] V(ω j )=A(ω j ) -1 W T δD(c,ω j ) (14)

[0106] Gradient Matrix The nth element of is calculated as follows:

[0107]

[0108] in, The subscript represents the element located at the (s-1)N+nth row and rth column.

[0109] In summary, to calculate each frequency ω j The gradient component of , needs 2M times of forward modeling, of which M times are used to obtain the forward wave field P(ω) according to formula (4), and the other M times are used to obtain the residual wave field V(ω), as shown in formula (14). Therefore, the gradient matrix The required computational effort is affordable and can be further improved through parallel computing methods.

[0110] Most importantly, the Frechet partial derivative matrix J is the wavefield generated when the receiving node collects all nodes as a single scattering point, while the residual matrix δD(c) is the actual data received by the sensor from the unknown scattering point. Therefore, as shown in formula (8), the gradient matrix is ​​actually the result of the zero-lag correlation between the corresponding scattering signal at each node and the actual scattering signal. Naturally, the correlation value is large at the actual scattering point but small at the non-scattering point, which occupies most of the measurement area. Therefore, the gradient matrix can be used as an image of the reflection point, and formula (8) can be used as an imaging condition.

[0111] 1.3 Imaging conditions of search vectors

[0112] The negative gradient direction can be directly used to minimize the objective function (7). However, in order to speed up the search, the Newton method is usually used to correct the gradient in the full waveform inversion, which is expressed as:

[0113]

[0114] Where dc is the search vector, H is the N×N Hessian matrix, and the elements of H are obtained by equation (17):

[0115]

[0116] Since only the defect position is unknown for ultrasonic imaging, the input sound velocity model is close to the actual sound velocity distribution, and the Hessian matrix can be approximated as:

[0117] H≈H a =J T J * (18)

[0118] From formulas (9) and (18), it can be seen that H a The main diagonal elements of are the zero-lag values ​​of the autocorrelation of one partial-guided wavefield at the receiving node, while the off-diagonal elements are the zero-lag values ​​of the cross-correlation of two partial-guided wavefields at the receiving node. In the high-frequency limit, H a It is diagonally dominant.

[0119] Therefore, H a The diagonal elements of are used to approximate the Hessian matrix to avoid the inverse operation of the impractical large N-dimensional matrix. The imaging condition of the search vector is written as:

[0120]

[0121] Among them, Im is the imaging result; diag(H a ) indicates H aTake the diagonal; I is an N-dimensional identity matrix; γ is a damping coefficient used to stabilize the inversion.

[0122] As mentioned in Section 1.2, the gradient matrix can be directly used as the image of the defect.

[0123] According to formula (8), the gradient matrix depends on the Frechet partial derivative matrix J, which is composed of all the partial derivative wave fields at the receiving nodes. Therefore, the characteristics of the partial derivative matrix J itself will affect the imaging results. On the one hand, the energy of the partial derivative wave field corresponding to the node at a deeper position is always smaller than that at a shallower position, so the partial derivative matrix J will cause insufficient illumination of the lower part of the imaging area. On the other hand, the different relative positions of the nodes and the sound source node and the receiving node will also cause different energies of the partial derivative wave field, so the partial derivative matrix J will also cause uneven illumination of the imaging area in the horizontal direction. The introduction of the Hessian matrix can eliminate the insufficient and uneven illumination caused by the Frechet partial derivative matrix J in the gradient, making the imaging results clearer and more focused, and balancing the imaging amplitudes at each position.

[0124] According to formula (18), the partial derivative matrix J needs to be calculated to calculate H a However, it is impractical to directly calculate the partial derivative matrix J using formula (12) because it requires too many forward calculations.

[0125] Starting from the definition of the partial derivative matrix J, formula (9) can be further written as:

[0126]

[0127] Among them, f (l，n) is the virtual source vector corresponding to the lth data and the nth parameter;

[0128] is an N×1 vector, which is only in the (x r ,0) nodes are 1, and all other nodes are 0. Therefore, J l，n The physical meaning of is the virtual source vector f (l，n) The excited sound field is (x r , the corresponding received signal at the node 0). In fact, f (l，n) Indicates the sound source corresponding to node n.

[0129] According to the principle of reciprocity, the nth node sends the r ,0) node receives the wave field signal equal to the signal received by the (x r , the wavefield signal transmitted by the 0)th node and received by the nth node.

[0130] Therefore, J l，n It can be obtained by the following formula:

[0131]

[0132] in,

[0133] For S(x r , x s ,ω j ), only in x r =x s There is an incentive value of a(ω j ), at this time a(ω j )=S(x s , x s ,ω j ); when x r ≠x s When S(x r , x s ,ω j )=0.

[0134]

[0135]

[0136]

[0137] Then we get: represents the unit forward wave field, in formula (21) It is obtained from the unit forward wave field.

[0138] The partial derivative matrix of a certain data (the lth data) with respect to all sound speed parameters can be directly calculated by the following formula:

[0139]

[0140] Among them, J l is the N×1 partial derivative matrix of the lth data;

[0141] f (l，n) is the virtual source vector corresponding to the lth data and the nth parameter, which is obtained from the total virtual source matrix, where: the total virtual source matrix F(ω)=[F (1) (ω),...,F (N) (ω)]

[0142] For frequency ω j exist:

[0143]

[0144] In the above formula, the right side of the equation replicates the virtual source matrix 64 times in the row direction;

[0145] In this way, in order to obtain the entire partial derivative matrix of a data, only one forward calculation is required, which reduces the computational complexity by 1 / N and makes the calculation of the partial derivative matrix J feasible;

[0146] H a The N diagonal elements of are calculated as follows:

[0147]

[0148] in, denotes Hadamard multiplication, l = (j-1)M 2 +(s-1)M+r.

[0149] The time zero imaging condition is a traditional imaging condition commonly used for self-transmitted and self-received signals. It assumes that the reflected signal is actively emitted by the reflection point at time 0, so the image of the reflection point can be obtained by reconstructing the wave field at time 0. For FMC data, the imaging condition is modified to the excitation time imaging condition, where the excitation time is the time when the incident wave reaches the reflector. Conventional imaging methods, including the total focusing method (TFM), phase shift method (PSM), and reverse time migration (RTM), all use the excitation time imaging condition to image the defect. The search vector imaging condition is constructed from the perspective of full waveform inversion, as shown in equation (19). This imaging condition is strictly derived from the formula and does not make any assumptions about the wave field. Therefore, its theory is more rigorous and can ensure higher imaging accuracy.

[0150] The high-precision imaging method of the present invention using ultrasonic FMC data sets to detect internal defects of complex-shaped parts can be summarized into three steps: preprocessing of input data, gradient matrix in frequency domain, and the approximate Hessian matrix diag(H a ) and the realization of the search vector imaging conditions.

[0151] The input data includes the excitation signal, the sound velocity model, and the FMC data (full matrix data). The excitation signal of the probe is usually unknown, so the excitation signal should be approximated by an appropriate wavelet, and its spectrum distribution should be as similar as possible to the obtained FMC data; the discretization of the sound velocity model is mainly determined by the finite difference method; the FMC data D(x r , x s , t) contains the reflection signals from various interfaces and defects in the measurement area. Assuming that the sound velocity model is smooth and no reflection occurs, the model data The residual data set δD(c) = D. The advantage of this is that all reflection interfaces can be clearly imaged, so the residual matrix δD(c, ω j ).

[0152] The second step is the main body of the whole method. This step is implemented in the frequency domain, so a specific calculation frequency range [ω min ,ω max ], the discrete step length dω is determined by the sampling frequency and the number of sampling points, and the number of frequencies is N ω At frequency ω j At, according to formulas (4) and (14), by the impedance matrix A(ω j ) is decomposed into LU, and the source matrix and residual matrix are processed respectively to calculate the forward wave field P(ω j ) and the residual wave field V(ω j ). In addition, as shown in formula (11), the forward wave field P(ω j ) is multiplied by all partial impedance matrices to obtain the total virtual source matrix F(ω j )=[F (1) (ω j ), ..., F (N) (ω j )]. In addition, P(ω j ) by the excitation value a(ω j ) is normalized to obtain the unit forward wave field Forward wave field and total virtual source matrix F(ω j ) are multiplied, as shown in formula (22), to calculate M 2 Deflection vector l=(j-1)M 2 +1. At the same time, the total virtual source matrix F(ω j ) is also related to the residual wave field V(ω j ) are multiplied, as shown in formula (13), to obtain the gradient matrix The above operation is for each discrete frequency ω j Therefore, all discrete frequencies are traversed according to the above steps. After the traversal is completed, the gradient matrix of all discrete frequencies is obtained by summing according to formulas (15) and (23). and the approximate Hessian matrix diag(H a ).

[0153] Finally, the imaging condition of the search vector is executed according to formula (19), and the imaging result Im is output, that is, a high-precision image of the internal defects of the part is obtained.

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

[0155] The present invention provides a high-precision ultrasonic imaging method for internal defects of complex-shaped parts based on search vector imaging conditions. Starting from full waveform inversion, the search vector imaging conditions are constructed. By combining the first-order derivative vector with the approximate Hessian matrix correction, high-quality images of defects in complex sound velocity models can be obtained, which is suitable for defect detection of complex surface components. The imaging conditions used in the imaging method of the present invention do not make assumptions about the wave field and are more rigorous than traditional imaging condition theories, which ensures the high precision of the inventive method. In addition, the imaging method of the present invention can highlight defects in the entire measurement area, which is very convenient for locating defects in large measurement structures. Compared with traditional methods, the imaging method of the present invention significantly improves imaging resolution and signal-to-noise ratio (SNR). BRIEF DESCRIPTION OF THE DRAWINGS

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

[0157] Figure 2 (a) is a schematic diagram of the discretization of the sound field; (b) is a schematic diagram of the discretization of the sound velocity model;

[0158] Figure 3 This is a schematic diagram of phased array detection of part defects;

[0159] Figure 4 Schematic diagram of defect detection in a curved aluminum component according to an embodiment of the present invention;

[0160] Figure 5 Comparisons of overall and local imaging results obtained using different imaging methods; (a)-(d) are the overall imaging results of TFM, FD-RTM, GDM, and the imaging method of an embodiment of the present invention, respectively; (e)-(h) are enlarged views of the local areas of (a)-(d), respectively. DETAILED DESCRIPTION

[0161] like Figure 1 As shown, a high-precision ultrasonic imaging method for internal defects of complex-shaped parts based on search vector imaging conditions includes the following steps:

[0162] (1) Input excitation signal S(x r , x s , t), the sound velocity model c(x, z) and the measured full matrix data D(x r , x s , t); where x r represents the rth receiving element, x s Represents the sth transmitting oscillator, r and s are both taken from [1, M], M is the number of oscillators; t represents time.

[0163] (2) Perform one-dimensional Fourier transform on the excitation signal and the full matrix data in time to obtain the sound source matrix S(x r , x s ,ω) and full matrix data D(x r , x s ,ω);

[0164] The calculation formula of the residual matrix δD(c) is as follows:

[0165]

[0166] Among them, the measurement value of the full matrix data D, that is, D(x r , x s ,ω); is the estimated value of the full matrix data;

[0167] Assuming that the sound velocity model is smooth and no reflection occurs, the estimated value of the full matrix data is 0, and the residual matrix δD(c, ω) = D(x r , x s , ω), ω represents the frequency.

[0168] The sound velocity model is discretized using the finite difference method to obtain the sound velocity vector, and the impedance matrix is ​​calculated using the following formula:

[0169]

[0170] Where A(ω) is the N×N impedance matrix; represents the Laplace operator obtained after discretization; c is the N×1 sound velocity vector.

[0171] (3) For any untraversed frequency ω j , using the impedance matrix A(ω j ) of the LU decomposition of the sound source matrix S(x r , x s ,ω j ) and the residual matrix D(x r , x s ,ω j ) are solved respectively to obtain the forward wave field and residual wave field of the frequency;

[0172] The solution formula for the forward wave field is as follows:

[0173] P(ω)=[p1(ω), p2(ω),...,p M (ω)]=A(ω) -1 S(ω)

[0174] Where S(ω) represents the sound source matrix, frequency ω jThe corresponding sound source matrix is ​​S(x r , x s ,ω j ); P(ω) represents the forward wave field, frequency ω j The corresponding forward wave field is P(ω j );

[0175] The solution formula for the residual wave field is as follows:

[0176] V(ω j )=A(ω j ) -1 WTδD(c,ω j )

[0177] Where, V(ω j ) represents the frequency ω j The corresponding residual wave field; δD(c, ω j ) represents the frequency ω j The corresponding residual matrix, j∈[1, N ω ],N ω represents the number of discrete frequencies; W is an M×N matrix, which extracts the sound field value at the node corresponding to the receiving array element in the entire wave field. In the kth row, only the element at the node corresponding to the kth receiving array element is 1, and the other elements are 0, where k∈[1,M] and M<N; the superscript T represents the transpose.

[0178] (4) The impedance matrix is ​​derived from the sound velocity vector and then multiplied by the forward wave field to obtain the total virtual source matrix. The total virtual source matrix is ​​defined as:

[0179] F(ω)=[F (1) (ω),,..,F (N) (ω)]

[0180] in,

[0181]

[0182] Where, F (n) (ω) is regarded as the nth virtual source, which is an N×M matrix; c n Represents the nth sound speed parameter (element) in the sound speed vector; n∈[1,N];

[0183] At the same time, the forward wave field is normalized by the amplitude of the excitation value of the excitation signal according to the following formula to obtain the unit forward wave field:

[0184]

[0185] Where, is the unit forward wave field; a(ω j )=S(x s , xs ,ω j ), for S(x r , x s ,ω j ), only in x r =x s There is an incentive value of a(ω j ); when x r ≠x s When S(x r , x s ,ω j )=0.

[0186] (5) Multiply the total virtual source matrix with the unit forward wave field and the residual wave field respectively to obtain the partial derivative matrix and gradient matrix corresponding to the frequency;

[0187] The calculation formula of the partial derivative matrix obtained by multiplying the total virtual source matrix and the unit forward wave field is as follows:

[0188]

[0189] Where, J l is the N×1 partial derivative matrix of the lth data, where l=[1,...,L], L=N ω ×M×M, represents the number of data points; l=(j-1)M 2 +(s-1)M+r;

[0190] f (l，n) is the virtual source vector corresponding to the lth data and the nth parameter, which is obtained from the total virtual source matrix, where: the total virtual source matrix F(ω)=[F (1) (ω),...,F (N) (ω)]

[0191] For frequency ω i exist:

[0192]

[0193] In the above formula, the right side of the equation replicates the virtual source matrix 64 times in the row direction;

[0194] By unit forward wave field get:

[0195]

[0196] The calculation formula of the gradient matrix obtained by multiplying the total virtual source matrix and the residual wave field is:

[0197]

[0198] In the formula, F(ωj ) represents the frequency ω j The corresponding total virtual source matrix.

[0199] (6) Repeat steps (3) to (5) until all frequencies are traversed, obtain the partial derivative matrix and gradient matrix of all frequencies, and use the imaging conditions of the search vector to obtain the imaging results of the internal defects of the part;

[0200] The imaging condition of the search vector is:

[0201]

[0202] Where Im represents the imaging result; I is an N-dimensional unit matrix; represents the N×1 gradient matrix; γ is the damping coefficient, which is used to stabilize the inversion;

[0203] H a =J T J * , where J represents the partial derivative matrix; superscript T and * represent transpose and complex conjugate respectively;

[0204] The nth element in is obtained by the following formula:

[0205]

[0206] Where, The subscript represents the element located in the (s-1)N+nth row and the rth column; n∈[1,N];

[0207] According to all partial derivative matrices, H is calculated by the following formula a The N diagonal elements of :

[0208]

[0209] Where, diag(H a ) represents the approximate Hessian matrix; denotes Hadamard multiplication, l = (j-1)M 2 +(s-1)M+r.

[0210] The theoretical derivation process of the above ultrasound imaging method is as follows:

[0211] 1.1 Frequency Wave Field

[0212] The frequency domain wave field in inhomogeneous media satisfies the following equation:

[0213]

[0214] Among them, x and z are two spatial coordinates in two-dimensional space, ω is frequency, z is depth, p(x, z, ω) is the sound field (wave field), s(x, z, ω) is the sound source field, and c(x, z) is the sound speed model.

[0215] In engineering practice, formula (1) is usually solved by the finite difference method, where the wave field and sound velocity model are discretized as N = N x ×N z points, such as Figure 2 shown.

[0216] The discrete wave equation is as follows:

[0217]

[0218] Among them, c is the sound velocity vector, p(ω) is the wave field vector, and s(ω) is the sound source vector, and all three are N×1 vectors;

[0219] is an N×N impedance matrix, where represents the Laplace operator obtained after discretization.

[0220] Therefore, the wave field p(ω) can be obtained by the following formula:

[0221] p(ω)=A(ω) -1 s(ω) (3)

[0222] Among them, A(ω) -1 Indicates the inversion of A(ω); the calculation of the wave field can introduce the LU decomposition of the impedance matrix to improve the calculation efficiency.

[0223] The schematic diagram of phased array defect detection is as follows: Figure 3 As shown, each array element is excited in turn, and all array elements receive signals at the same time. Therefore, the FMC data (full matrix data) contains the time domain signal of each transmit-receive pair. Here, D(x r , x s , t) represents, where x r and x s They correspond to the rth receiving array element and the sth transmitting array element respectively; s∈[1,M], r∈[1,M], where M is the number of elements. Figure 2 As shown, the ultrasonic signal is transmitted from the array element (x s ,0) is emitted and reflected by the interface and defects and is located at (x r Each element in the phased array can be represented by a node corresponding to its center point in the discrete model, so the wave fields at these nodes can be known from the FMC dataset.

[0224] Using the sound source matrix S(ω)=[s1(ω),s2(ω),...,s M (ω)] represents the sound source of a phased array with M elements. The wave field corresponding to the excitation of these M elements can be solved at one time by the following formula:

[0225] P(ω)=[p1(ω), p2(ω),,..,p M (ω)]=A(ω) -1 S(ω) (4)

[0226] Among them, P(ω) is the wave field matrix of the full matrix data.

[0227] 1.2 Full Waveform Inversion and Gradient Vectors

[0228] Given the sound source matrix and the sound velocity model, the wave field matrix of the FMC data can be calculated according to Equation (4), and the frequency domain FMC data can be obtained by taking the wave field values ​​at the corresponding nodes of the receiver, which is written as:

[0229]

[0230] in, Represents the estimated value of the full matrix data; W is an M×N matrix, which extracts the sound field value at the node corresponding to the receiving array element in the entire wave field. In the kth row, only the element at the node corresponding to the kth receiving array element is 1, and the other elements are 0, where k∈[1,M] and M<N; is an M×M matrix, and its element at the (r, s) position corresponds to the signal received by the rth element when the sth element is excited. By setting the Integration is obtained.

[0231] The calculated full matrix data The difference between the measured full matrix data D is defined as the residual matrix, written as:

[0232]

[0233] Among them, δD(c) and They are all functions of the sound velocity vector c. The sound velocity vector c that minimizes the modulus of the residual matrix δD(c) is the closest to the true sound velocity, which is the basis of full waveform inversion.

[0234] Here, the l2 norm of the residual data δD(c) is used to define the objective function:

[0235]

[0236] Wherein, the superscripts T and * represent transpose and complex conjugate respectively; Nω is the number of discrete frequencies; M represents the number of oscillators; D(x r , x s ,ω) by D(x r , x s , t) is obtained by Fourier transform in the time dimension.

[0237] The negative gradient direction can be used to update the sound speed vector, written as:

[0238]

[0239] in, represents the gradient matrix; J is the L×N Frechet partial derivative matrix, and its elements are as follows:

[0240]

[0241] Where L = N ω ×M×M number of data points, l=(j-1)M 2 +(s-1)M+r;j∈[1,N ω ].

[0242] Apply both sides of formula (4) to the nth sound speed parameter c n Taking partial derivatives we can get formula (10):

[0243]

[0244] in, is an N×N partial derivative impedance matrix, which is the impedance matrix for the nth sound velocity parameter c n Find the result of partial derivative.

[0245] The form of formula (10) is similar to formula (4), F (n) (ω) can be regarded as the nth virtual source, which is an N×M matrix and written as:

[0246]

[0247] therefore, is an N×M matrix, which is the solution to a forward problem.

[0248] Then, at frequency ω j The Frechet partial derivative matrix when can be written as:

[0249]

[0250] Among them, J(ω j ) is an M×MN matrix, and the entire Frechet partial derivative matrix J can be obtained by integrating all discrete frequencies ωj The corresponding J(ω j )get.

[0251] The total virtual source matrix can be defined as F(ω)=[F (1) (ω),...,F (N) (ω)], which is an N×MN matrix.

[0252] Explicitly calculating the partial derivative matrix J is very time-consuming. Substituting formula (12) into formula (8), the gradient matrix can be directly calculated:

[0253]

[0254] Among them, δD(c,ω j ) is the frequency ω j The M×M residual data matrix when , The frequency is ω j The MN×M gradient matrix when .

[0255] Residual wave field V(ω j ) is an N×M matrix and is calculated as follows:

[0256] V(ω j )=A(ω j ) -1 W T δD(c,ω j ) (14)

[0257] Gradient Matrix The nth element of is calculated as follows:

[0258]

[0259] in, The subscript represents the element located at the (s-1)N+nth row and rth column.

[0260] In summary, to calculate each frequency ω j The gradient component of , needs 2M times of forward modeling, of which M times are used to obtain the forward wave field P(ω) according to formula (4), and the other M times are used to obtain the residual wave field V(ω), as shown in formula (14). Therefore, the gradient matrix The required computational effort is affordable and can be further improved through parallel computing methods.

[0261] Most importantly, the Frechet partial derivative matrix j is the wavefield generated when the receiving node collects all nodes as a single scattering point, while the residual matrix δD(c) is the actual data received by the sensor from the unknown scattering point. Therefore, as shown in formula (8), the gradient matrix is ​​actually the result of the zero-lag correlation between the corresponding scattering signal at each node and the actual scattering signal. Naturally, the correlation value is large at the actual scattering point but small at the non-scattering point, which occupies most of the measurement area. Therefore, the gradient matrix can be used as an image of the reflection point, and formula (8) can be used as an imaging condition.

[0262] 1.3 Imaging conditions of search vectors

[0263] The negative gradient direction can be directly used to minimize the objective function (7). However, in order to speed up the search, the Newton method is usually used to correct the gradient in the full waveform inversion, which is expressed as:

[0264]

[0265] Where dc is the search vector, H is the N×N Hessian matrix, and the elements of H are obtained by equation (17):

[0266]

[0267] Since only the defect position is unknown for ultrasonic imaging, the input sound velocity model is close to the actual sound velocity distribution, and the Hessian matrix can be approximated as:

[0268] H≈H a =J T J * (18)

[0269] From formulas (9) and (18), it can be seen that H a The main diagonal elements of are the zero-lag values ​​of the autocorrelation of one partial-guided wavefield at the receiving node, while the off-diagonal elements are the zero-lag values ​​of the cross-correlation of two partial-guided wavefields at the receiving node. In the high-frequency limit, H a It is diagonally dominant.

[0270] Therefore, H a The diagonal elements of are used to approximate the Hessian matrix to avoid the inverse operation of the impractical large N-dimensional matrix. The imaging condition of the search vector is written as:

[0271]

[0272] Among them, Im is the imaging result; diag(H a ) indicates H aTake the diagonal; I is an N-dimensional identity matrix; γ is a damping coefficient used to stabilize the inversion.

[0273] As mentioned in Section 1.2, the gradient matrix can be directly used as the image of the defect.

[0274] According to formula (8), the gradient matrix depends on the Frechet partial derivative matrix J, which is composed of all the partial derivative wave fields at the receiving nodes. Therefore, the characteristics of the partial derivative matrix J itself will affect the imaging results. On the one hand, the energy of the partial derivative wave field corresponding to the node at a deeper position is always smaller than that at a shallower position, so the partial derivative matrix j will cause insufficient illumination of the lower part of the imaging area. On the other hand, the different relative positions of the nodes and the sound source node and the receiving node will also cause different energies of the partial derivative wave field, so the partial derivative matrix j will also cause uneven illumination of the imaging area in the horizontal direction. The introduction of the Hessian matrix can eliminate the insufficient and uneven illumination caused by the Frechet partial derivative matrix j in the gradient, making the imaging results clearer and more focused, and balancing the imaging amplitudes at each position.

[0275] According to formula (18), the partial derivative matrix j needs to be calculated to calculate H a , however, it is impractical to directly calculate the partial derivative matrix j through formula (12) because it requires too many forward calculations.

[0276] Starting from the definition of the partial derivative matrix j, formula (9) can be further written as:

[0277]

[0278] Among them, f (l，n) is the virtual source vector corresponding to the lth data and the nth parameter;

[0279] is an N×1 vector, which is only in the (x r ,0) nodes are 1, and all other nodes are 0. Therefore, J l，n The physical meaning is that the virtual source vector f (l，n) The excited sound field is (x r , the corresponding received signal at the node 0). In fact, f (l，n) Represents the sound source corresponding to node n.

[0280] According to the principle of reciprocity, the nth node sends the r ,0) node receives the wave field signal equal to the signal received by the (x r , the wavefield signal transmitted by the 0)th node and received by the nth node.

[0281] Therefore, J l，n It can be obtained by the following formula:

[0282]

[0283] in,

[0284] For S(x r , x s ,ω j ), only in x r =x s There is an incentive value of a(ω j ), at this time a(ω j )=S(x s , x s ,ω j ); when x r ≠x s When S(x r , x s ,ω j )=0.

[0285]

[0286]

[0287]

[0288] Then we get:

[0289] The partial derivative matrix of a certain data (the lth data) with respect to all sound speed parameters can be directly calculated by the following formula:

[0290]

[0291] Among them, J l is the N×1 partial derivative matrix of the lth data;

[0292] f (l，n) is the virtual source vector corresponding to the lth data and the nth parameter, which is obtained from the total virtual source matrix, where: the total virtual source matrix F(ω)=[F (1) (ω),...,F (N) (ω)]

[0293] For frequency ω j exist:

[0294]

[0295] In the above formula, the right side of the equation replicates the virtual source matrix 64 times in the row direction;

[0296] In this way, in order to obtain the entire partial derivative matrix of a data, only one forward calculation is required, which reduces the computational complexity by 1 / N and makes the calculation of the partial derivative matrix J feasible.

[0297] H a The N diagonal elements of are calculated as follows:

[0298]

[0299] in, denotes Hadamard multiplication, l = (j-1)M 2 +(s-1)M+r.

[0300] The time zero imaging condition is a traditional imaging condition commonly used for self-transmitted and self-received signals. It assumes that the reflected signal is actively emitted by the reflection point at time 0, so the image of the reflection point can be obtained by reconstructing the wave field at time 0. For FMC data, the imaging condition is modified to the excitation time imaging condition, where the excitation time is the time when the incident wave reaches the reflector. Conventional imaging methods, including the total focusing method (TFM), phase shift method (PSM), and reverse time migration (RTM), all use the excitation time imaging condition to image the defect. The search vector imaging condition is constructed from the perspective of full waveform inversion, as shown in equation (19). This imaging condition is strictly derived from the formula and does not make any assumptions about the wave field. Therefore, its theory is more rigorous and can ensure higher imaging accuracy.

[0301] Detection experiment

[0302] The following experiments were conducted to verify the effectiveness of the imaging method in this embodiment for imaging internal defects of components with complex shapes.

[0303] The FMC data (full matrix data) were collected using a 64 / 64OEM-PA (AOS. Inc., USA) with a sampling frequency of 50 MHz and a time range of 80 μs. An aluminum component with a circular surface was the measurement object, such as Figure 4 As shown, cracks and voids are the measurement targets. A 64-element array transducer (Guangdong Gaohua Co., Ltd., China) with a center frequency of 2.5 MHz and a pitch of 0.75 mm was used as an ultrasonic transceiver. The components were immersed in water. Since the phased array is not an immersion probe, a wedge made of acrylonitrile butadiene styrene (ABS) was used to protect the probe from water. The sound velocity of the wedge was 2200 m / s. The bottom of the wedge is a concave surface to focus more ultrasonic waves into the measured component. Water, as a coupling agent, filled the gap between the wedge and the component but did not reach the top surface of the wedge.

[0304] The experimental data (full matrix data) were processed using TFM, FD-RTM, GDM and the imaging method of this embodiment. The frequency bands (frequency ranges) selected by these methods are all from 1 MHz to 4 MHz. The parameter γ of the imaging method of this embodiment was selected as 0.05. The overall imaging results of the four methods are shown in Figure 2. Figure 5 As shown in (a)-(d), Figure 5 The results of the TFM and imaging methods of this embodiment in (a) and (d) are very easy to find defects; however, due to insufficient lighting, Figure 5 Almost no defect imaging results can be found in the FD-RTM and GDM results in (b) and (c).

[0305] also, Figure 5 (e)-(h) show the local magnified images of the imaging results of the above four imaging methods. Figure 5 The defect images in (f) and (g) FD-RTM and GMD results are still difficult to distinguish, and Figure 5 Defect contrast ratio in the GMD imaging result in (e) Figure 5 The contrast of the FD-RTM imaging result in (f) is high. Figure 5 In (e), the artifacts around the defect are severe and the crack is somewhat distorted. Figure 5 The imaging result of the imaging method of this embodiment in (h) can accurately image the crack morphology and effectively suppress the generation of background artifacts.

Claims

1. A high-precision ultrasonic imaging method for internal defects of complex-shaped parts based on search vector imaging conditions, characterized in that: The following steps are involved: (1) Input the excitation signal, sound velocity model and full matrix data of the measurement; (2) Preprocess the excitation signal and full matrix data separately to obtain the sound source matrix and residual matrix in the frequency domain; Discretize the sound velocity model to obtain the sound velocity vector, and calculate the impedance matrix; (3) For any untraversed frequency, the LU decomposition of the impedance matrix is ​​used to solve the source matrix and the residual matrix respectively to obtain the forward wave field and residual wave field of the frequency; (4) The impedance matrix is ​​partial-derived with respect to the acoustic velocity vector and then multiplied with the forward wave field to obtain the total virtual source matrix. At the same time, the forward wave field is normalized by the amplitude of the excitation value of the excitation signal to obtain the unit forward wave field. (5) Multiply the total virtual source matrix with the unit forward wave field and the residual wave field respectively to obtain the partial derivative matrix and gradient matrix corresponding to the frequency; (6) Repeat steps (3) to (5) until all frequencies are traversed, and obtain the partial derivative matrix and gradient matrix of all frequencies. Use the imaging conditions of the search vector to obtain the imaging results of the internal defects of the part; The imaging condition of the search vector in step (6) is: ； in, Indicates imaging results; For one N -dimensional identity matrix; represents the gradient matrix; is the damping coefficient, used to stabilize the inversion; ,in, Denotes the partial derivative matrix; superscript T and denote transpose and complex conjugate, respectively.

2. The high-precision ultrasonic imaging method for internal defects of complex-shaped parts based on search vector imaging conditions according to claim 1 is characterized in that: In step (2), the residual matrix The calculation formula is as follows: ； in, Measurement values ​​of full matrix data; is the estimated value of the full matrix data; Assuming that the sound velocity model is smooth and no reflection occurs, the estimated value of the full matrix data is 0, and the residual matrix ; The full matrix data in the measured time domain is subjected to a one-dimensional Fourier transform in time to obtain the full matrix data in the frequency domain, that is, to obtain the residual matrix in the frequency domain.

Images