HIGH-PRECISION ULTRASONIC IMAGING METHOD FOR INTERNAL DEFECT OF COMPLEX-SHAPED COMPONENT BASED ON SEARCH-VECTOR IMAGING CONDITION

Abstract

The present disclosure provides a high-precision ultrasonic imaging method for an internal defect of a complex-shaped component based on a search-vector imaging condition. Based on full waveform inversion, the search-vector imaging condition is constructed, and a first-order derivative vector is combined with an approximate Hessian matrix for correction so that a high-quality image of a defect in a complex sound velocity model can be obtained. The high-precision ultrasonic imaging method is suitable for defect testing on complex surface components. The imaging condition used for the imaging method of the present disclosure does not involve an assumption on a wave field and is stricter than a traditional imaging condition theory, thereby ensuring high precision of the inventive method. Moreover, the imaging method of the present disclosure can highlight a defect in an entire measurement area, which is very convenient for locating a defect in a large-size measurement structure.

Claims

1. A high-precision ultrasonic imaging method for an internal defect of a complex-shaped component based on a search-vector imaging condition, comprising the following steps: (1) inputting an excitation signal, a sound velocity model, and full matrix capture (FMC) data measured; (2) preprocessing the excitation signal and the FMC data to obtain a sound source matrix and a residual matrix within a frequency domain, respectively; discretizing the sound velocity model to obtain a sound velocity vector, and obtaining an impedance matrix by calculation; (3) for any frequency not traversed, solving the sound source matrix and the residual matrix by LU decomposition of the impedance matrix to obtain a forward wave field and a residual wave field of the frequency, respectively; (4) multiplying the impedance matrix by the forward wave field after obtaining a partial derivative of the sound velocity vector to obtain a total virtual source matrix, and normalizing the forward wave field by an amplitude of an excitation value of the excitation signal to obtain a unit forward wave field; (5) multiplying the total virtual source matrix by the unit forward wave field and the residual wave field to obtain a partial derivative matrix and a gradient matrix corresponding to the frequency, respectively; and (6) repeating steps (3) to (5) until all frequencies are traversed completely to obtain partial derivative matrices and gradient matrices of all the frequencies, and obtaining an imaging result of an internal defect of a complex-shaped component by using a search-vector imaging condition.

2. The high-precision ultrasonic imaging method for an internal defect of a complex-shaped component based on a search-vector imaging condition according to claim 1, wherein in step (2), a calculation formula for the residual matrix ?D(c) is as follows: $?D\left(c\right)=D-\widehat{D}\left(c\right)$ wherein D represents a measured value of the FMC data; and {circumflex over (D)}(c) represents an estimated value of the FMC data; assuming that the sound velocity model is smooth without reflection, the estimated value of the FMC data is 0, and the residual matrix is ?D(c)=D; and temporal one-dimensional Fourier transformation is performed on the measured FMC data within a time domain to obtain FMC data within a frequency domain, namely the residual matrix within the frequency domain.

3. The high-precision ultrasonic imaging method for an internal defect of a complex-shaped component based on a search-vector imaging condition according to claim 1, wherein the search-vector imaging condition in step (6) is as follows: $\mathrm{Im}=-{\left[\mathrm{diag}\left(H_{?}\right)+?I\right]}^{-1}{?}_{c}E\left(c\right)$ wherein Im represents the imaging result; I represents an N-dimension unit matrix; ?.sub.cE(c) represents the gradient matrix; and ? represents a damping coefficient for stable inversion; and H.sub.a=J.sup.TJ*, wherein J represents the partial derivative matrix; and superscripts T and * represent transposition and complex conjugation, respectively.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0114] FIG. 1 is a flowchart of an imaging method according to an embodiment of the present disclosure;

[0115] FIG. 2A is a schematic diagram of sound field discretization, and FIG. 2B is a schematic diagram of sound velocity model discretization;

[0116] FIG. 3 is a schematic diagram of phased array testing of a component defect;

[0117] FIG. 4 is a schematic diagram of defect testing in a curved-surface aluminum component in an embodiment of the present disclosure; and

[0118] FIGS. 5A-H are comparison diagrams of overall imaging results and local imaging results obtained by different imaging methods, where FIG. 5A to FIG. 5D are diagrams of overall imaging results of TFM, FD-RTM, GDM, and the imaging method in the embodiment of the present disclosure; and FIG. 5E to FIG. 5H are partially enlarged diagrams of FIG. 5A to FIG. 5D.

[0119] Reference numerals in FIG. 1: (1) indicates S(x.sub.r,x.sub.s,?.sub.j), (2) indicates c(x,y), (3) indicates D(x.sub.i,x.sub.s,t.sub.j), (4) indicates S(x.sub.r,x.sub.s,?.sub.j), (5) indicates A(w.sub.j), (6) indicates ?D(c,w.sub.j), (7) indicates P(w.sub.j), (8) indicates ?A(w.sub.j)/?c.sub.1, . . . , ?A(w.sub.j)/?c.sub.N, (9) indicates V(w.sub.j), (10) indicates {tilde over (P)}(w.sub.j)=P(w.sub.k)/a(w.sub.j), (11) indicates F(w.sub.j)=[F.sup.(1)(w.sub.j), . . . , F.sup.(N)(w.sub.j)], (12) indicates [J.sub.l, . . . , J.sub.l+M.sub.2.sub.?1], l=(j?1)M.sup.2+1, (13) indicates ?.sub.cE(c,w.sub.j), (14) indicates

[00032]$\begin{array}{cc}{?}_c{E\left(c\right)}_n=\underset{j=1}{\overset{N_{?}}{.Math.}}\underset{s=1}{\overset{M}{.Math.}}\underset{r=1}{\overset{M}{.Math.}}{?}_c{E\left(c,w_j\right)}_{\left(s-1\right)N+n,r},& \left(15\right)\end{array}$$\mathrm{diag}\left(H_a\right)=\stackrel{N_{?}}{\underset{j=1}{.Math.}}\underset{s=1}{\overset{M}{.Math.}}\underset{r=1}{\overset{M}{.Math.}}{J}_l?J_l,$

indicates

[00033]$Im=-{\left[\mathrm{diag}\left(H_a\right)+?I\right]}^{-1}{?}_{c}E\left(c\right).$

DETAILED DESCRIPTION OF THE EMBODIMENTS

[0120] As shown in FIG. 1, a high-precision ultrasonic imaging method for an internal defect of a complex-shaped component based on a search-vector imaging condition includes the following steps.

[0121] (1) An excitation signal S(x.sub.r,x.sub.s,t), a sound velocity model c(x,z), and FMC data D(x.sub.r,x.sub.s,t) measured, where x.sub.r represents the rth receiving array element; x.sub.s represents the sth transmitting array element; r and s are both taken from [1,M]; M represents a number of array elements; and t represents a time.

[0122] (2) Temporal one-dimensional Fourier transformation is performed on the excitation signal and the FMC data to obtain a sound source matrix S(x.sub.r,x.sub.s,?) and FMC data D(x.sub.r,x.sub.s,?) within a frequency domain, respectively.

[0123] A calculation formula for a residual matrix ?D(c) is as follows:

[00034]$?D\left(c\right)=D-\widehat{D}\left(c\right)$

[0124] where D represents a measured value of the FMC data, i.e., D(x.sub.r,x.sub.s,?); and {circumflex over (D)}(c) represents an estimated value of the FMC data.

[0125] Assuming that the sound velocity model is smooth without reflection, the estimated value of the FMC data is 0, and the residual matrix is ?D(c,?)=D(x.sub.r,x.sub.s,?), ? representing a frequency.

[0126] The sound velocity model is discretized by using a finite difference method to obtain the sound velocity vector, and the impedance matrix is obtained through calculation by the following formula:

[00035]$A\left(?\right)=\left[{?}^2+{?}^2/c^2\right]$

[0127] where A(?) represents an N?N impedance matrix; ?.sup.2 represents a laplace operator obtained after the discretization; and c represents an N?1 sound velocity vector.

[0128] (3) For any frequency ?.sub.j not traversed, the sound source matrix S(x.sub.r,x.sub.s,?.sub.j) and the residual matrix D(x.sub.r,x.sub.s,?.sub.j) are solved by LU decomposition of the impedance matrix A(?.sub.j) to obtain a forward wave field and a residual wave field of the frequency, respectively.

[0129] A solving formula for the forward wave field is as follows:

[00036]$P\left(?\right)=\left[p_1\left(?\right),p_2\left(?\right),.Math.,p_M\left(?\right)\right]={A\left(?\right)}^{-1}S\left(?\right)$

[0130] where S(?) represents the sound source matrix, and the sound source matrix corresponding to frequency ?.sub.j is S(x.sub.r,x.sub.s,?.sub.j); P(?) represents the forward wave field, and the forward wave field corresponding to frequency ?.sub.j is P(?.sub.j).

[0131] A solving formula for the residual wave field is as follows:

[00037]$V\left({?}_j\right)={A\left({?}_j\right)}^{-1}{W}^T?D\left(c,{?}_j\right)$

[0132] where V(?.sub.j) represents the residual wave field corresponding to frequency ?.sub.j; ?D(c,?.sub.j) represents the residual matrix corresponding to the frequency ?.sub.j, j?[1,N.sub.?]; N.sub.? represents a number of discrete frequencies; W represents a matrix of M?N, which extracts a sound field value at a node corresponding to a receiving array element in the entire wave field; an element at a node corresponding to the kth receiving array element in the kth row is 1, and other elements are all 0, where k?[1,M], and M<N; and the superscript T represents transposition.

[0133] (4) The impedance matrix is multiplied by the forward wave field after a partial derivative of the sound velocity vector is obtained to obtain a total virtual source matrix. The total virtual source matrix is defined as.

F(?)=[F.sup.(1)(?), . . . ,F.sup.(N)(?)]

[0134] where

[00038]$F^{\left(n\right)}\left(?\right)=-\frac{?A\left(?\right)}{?c_n}P\left(?\right)$

[0135] where F(W)(?) is regarded as the nth virtual source, which is an N?M matrix; c.sub.n represents the nth sound velocity parameter (element) in the sound velocity vector; and n?[1,N].

[0136] Meanwhile, the forward wave field is normalized by an amplitude of an excitation value of the excitation signal to obtain a unit forward wave field.

{tilde over (P)}(?.sub.j)=P(?.sub.j)/a(?.sub.j)

[0137] where {tilde over (P)}(?.sub.j) represents the unit forward wave field; a(?.sub.j)=S(x.sub.s,x.sub.s,?.sub.j), for S(x.sub.r,x.sub.s,?.sub.j), there is an excitation value a(?.sub.j) only at x.sub.r=x.sub.s; and when x.sub.r?x.sub.s, S(x.sub.r,x.sub.s,?.sub.j)=0.

[0138] (5) The total virtual source matrix is multiplied by the unit forward wave field and the residual wave field to obtain a partial derivative matrix and a gradient matrix corresponding to the frequency, respectively.

[0139] A calculation formula for multiplying the total virtual source matrix by the unit forward wave field to obtain the partial derivative matrix is as follows:

[00039]$J_l=\frac{?{\hat{D}}_l}{?c}={\left[f_{\left(l,1\right)},f_{\left(l,2\right)},.Math.,f_{\left(l,N\right)}\right]}^T{\tilde{p}}_{\left({?}_j,x_r\right)}$

[0140] where J.sub.l represents the N?1 partial derivative matrix of the lth data, where l=[1, . . . , L], L=N.sub.??M?M, representing a number of data points; and l=(j?1)M.sup.2+(s?1)M+r.

[0141] f.sub.(l,n) represents a virtual source vector corresponding to the lth data and the nth parameter and is obtained by the total virtual source matrix, where the total virtual source matrix is expressed as: F(?)=[F.sup.(1)(?), . . . , F.sup.(N)(?)].

[0142] For frequency ?.sub.j, the following can be derived:

[00040]$\left[f_{\left(l,n\right)},f_{\left(l+1,n\right)},.Math.,f_{\left(l+M^2-1,n\right)}\right]=\left[F^{\left(n\right)}\left({?}_j\right),F^{\left(n\right)}\left({?}_j\right),.Math.,F^{\left(n\right)}\left({?}_j\right)\right]$

[0143] In the above equation, the virtual source matrix is repeated for 64 times in a row direction on the right of the equation.

[0144] {tilde over (p)}.sub.(?.sub.j.sub.,x.sub.r.sub.) is derived from the unit forward wave field {tilde over (P)}(?.sub.j):

{tilde over (P)}(?.sub.j)=[{tilde over (p)}.sub.(?.sub.j.sub.,x.sub.1.sub.), . . . ,{tilde over (p)}.sub.(?.sub.j.sub.,x.sub.M.sub.)]

[0145] A calculation formula for multiplying the total virtual source matrix by the residual wave field to obtain the gradient matrix is as follows:

[00041]${?}_{c}E\left(c,{?}_j\right)={F\left({?}_j\right)}^T\left[{A\left({?}_j\right)}^{-1}{W}^T?D\left(c,{?}_j\right)\right]={F\left({?}_j\right)}^{T}V\left({?}_j\right)$

[0146] where F(?.sub.j) represents the total virtual source matrix corresponding to the frequency ?.sub.j.

[0147] (6) Steps (3) to (5) are repeated until all frequencies are traversed completely to obtain partial derivative matrices and gradient matrices of all the frequencies, and an imaging result of an internal defect of a complex-shaped component is obtained by using a search-vector imaging condition.

[0148] The search-vector imaging condition is as follows:

[00042]$Im=-{\left[\mathrm{diag}\left(H_a\right)+?I\right]}^{-1}{?}_{c}E\left(c\right)$

[0149] where Im represents the imaging result; I represents an N-dimension unit matrix; ?.sub.cE(c) represents an N?1 gradient matrix; and ? represents a damping coefficient for stable inversion.

[0150] H.sub.a=J.sup.TJ*, where J represents the partial derivative matrix; and superscripts T and * represent transposition and complex conjugation, respectively.

[0151] The nth element in ?.sub.cE(c) is obtained by the following formula:

[00043]${?}_c{E\left(c\right)}_n={.Math.}_{j=1}^{N_{?}}{.Math.}_{s=1}^M{.Math.}_{r=1}^M{?}_c{E\left(c,{?}_j\right)}_{\left(s-1\right)N+n,r}$

[0152] where the subscript of ?.sub.cE(c,?.sub.j).sub.(s-1)N+n,r represents an element in row (s?1)N+n and column r; and n?[1,N].

[0153] N diagonal elements of H.sub.a are obtained through calculation by the following formula according to all the partial derivative matrices:

[00044]$\mathrm{diag}\left(H_a\right)={.Math.}_{j=1}^{N_{?}}{.Math.}_{s=1}^M{.Math.}_{r=1}^M{J}_l?J_l$

[0154] where diag(H.sub.a) represents an approximate Hessian matrix; and a represents Hadamard multiplication, l=(j?1)M.sup.2+(s?1)M+r.

[0155] The theoretical derivation process of the ultrasonic imaging method is as follows.

1.1 Frequency Wave Field

[0156] A frequency-domain wave field in an inhomogeneous medium meets the following formula:

[00045]$\begin{array}{cc}|{?}^2+\frac{{?}^2}{{c\left(x,z\right)}^2}]p\left(x,z,?\right)=s\left(x,z,?\right)& \left(1\right)\end{array}$

[0157] where x and z represent two spatial coordinates of a two-dimensional space, respectively; ? represents a frequency; z represents a depth; p(x,z,?) represents a sound field (wave field); s(x,z,?) represents a sound source field; and c(x,z) represents a sound velocity model.

[0158] In engineering practice, the formula (1) is usually solved by the finite difference method, where the wave field and the sound velocity model are discretized into N=N.sub.x?N.sub.z points, as shown in FIG. 2.

[0159] A discretized wave equation is as follows:

[00046]$\begin{array}{cc}\left[{?}^2+\frac{{?}^2}{c^2}\right]p\left(?\right)=A\left(?\right)p\left(?\right)=s\left(?\right)& \left(2\right)\end{array}$

[0160] where c represents a sound velocity vector; p(?) represents a wave field vector; s(?) represents a sound source vector, and the three vectors are all N?1 vectors.

[0161] A(?)=[?.sup.2+?.sup.2/c.sup.2] is an N?N impedance matrix, where ?.sup.2 represents a laplace operator obtained after the discretization.

[0162] Therefore, the wave field p(?) may be obtained by the following formula:

[00047]$p\left(?\right)={A\left(?\right)}^{-1}s\left(?\right)$

[0163] where A(?).sup.?1 represents inverting A(?); and LU decomposition of the impedance matrix may be introduced into the calculation of the wave field to improve the calculation efficiency.

[0164] The schematic diagram of phased array testing of a defect is as shown in FIG. 3. Array elements are excited orderly, and meanwhile, all the array elements receive signals. Therefore, FMC data (full matrix capture data) contains a time-domain signal of each transmitting-receiving pair, which is represented by D(x.sub.r,x.sub.s,t) here, where x.sub.r and x.sub.s correspond to the rth receiving array element and the sth transmitting array element, respectively; s?[1,M], r?[1,M], and M represents a number of array elements. As shown in FIG. 2, an ultrasonic signal is emitted from array element (x.sub.s,0) and is received by an array element at position (x.sub.r,0) after being reflected by an interface and a defect. Each array element in a phased array can be represented by a node corresponding to the position of a center point in a discrete model, and wave fields at these nodes can be known from an FMC dataset.

[0165] The sound source matrix S(?)=[s.sub.1(?), s.sub.2(?), . . . , s.sub.M(?)] is used to represent a sound source of the phased array having M array elements, and corresponding wave fields excited for the M array elements may be solved through one step by the following formula:

[00048]$\begin{array}{cc}P\left(?\right)=\left[p_1\left(?\right),p_2\left(?\right),.Math.,p_M\left(?\right)\right]={A\left(?\right)}^{-1}S\left(?\right)& \left(4\right)\end{array}$

[0166] where P(?) represents a wave field matrix of the FMC data.

1.2 Full Waveform Inversion and Gradient Vector

[0167] Given the sound source matrix and the sound velocity model, the wave field matrix of the FMC data is calculated according to the formula (4); and frequency-domain FMC data is obtained by a wave field value at a corresponding node of a received, which is written as:

{circumflex over (D)}(c,?)=WP(?)(5)

[0168] where {circumflex over (D)}(c,?) represents an estimated value of the FMC data; W represents an M?N matrix, which extracts a sound field value at a node corresponding to a receiving array element in the entire wave field; an element at a node corresponding to the kth receiving array element in the kth row is 1, and other elements are all 0, where k?[1,M], and M<N; and {circumflex over (D)}(c,?) is an M?M matrix, and its element at position (r, s) corresponds to a signal received by the rth array element when the sth array element is excited. The entire FMC data {circumflex over (D)}(c) obtained through the calculation may be obtained by combining {circumflex over (D)}(c,?) of each frequency.

[0169] A difference between the calculated FMC data D(c) and the measured FMC data D is defined as a residual matrix, written as:

[00049]$\begin{array}{cc}?D\left(c\right)=D-\stackrel{?}{D}\left(c\right)& \left(6\right)\end{array}$

[0170] where ?D(c) and {circumflex over (D)}(c) are both functions of the sound velocity vector c. The sound velocity vector c that minimizes the modulus of the residual matrix ?D(c) is closest to the real sound velocity, which is the base of the full waveform inversion.

[0171] Here, an objective function is defined by using the l.sub.2 norm of residual data ?D(c):

[00050]$\begin{array}{cc}E\left(c\right)=\frac{1}{2}?{D\left(c\right)}^T?{d\left(c\right)}^{*}=\frac{1}{2}{.Math.}_{j=1}^{N_{?}}{.Math.}_{s=1}^M{{.Math.}_{r=1}^M\left[D\left(x_r,x_s,{?}_j\right)-\hat{D}\left(x_r,x_s,{?}_j\right)\right]}^2& \left(7\right)\end{array}$

[0172] where the superscripts T and * represent transposition and complex conjugation, respectively; N.sub.? represents a number of discrete frequencies; M represents the number of array elements; and D(x.sub.r,x.sub.s,?) is obtained by performing Fourier transformation on D(x.sub.r,x.sub.s,t) in time dimension.

[0173] A negative gradient direction may be used to update the sound velocity vector, written as:

?.sub.cE(c)=J.sup.T?D(c)(8)

[0174] where ?.sub.cE(c) represents a gradient matrix; J is an L?N Frechet partial derivative matrix, and an element thereof is expressed as:

[00051]$\begin{array}{cc}{J}_{l,n}=\frac{?{\hat{D}}_l}{?c_n}=\frac{?\hat{D}\left(x_r,x_s,{?}_j\right)}{?c_n},l=\left[1,.Math.,L\right];n=\left[1,.Math.,N\right]& \left(9\right)\end{array}$

[0175] where L=N.sub.??M?M represents a number of data points; l=(j?1)M.sup.2+(s?1)M+r; and j?[1,N].

[0176] Formula (10) may be derived by obtaining a partial derivative of the nth sound velocity parameter c.sub.n on both sides of the formula (4):

[00052]$\begin{array}{cc}\frac{?P\left(?\right)}{?c_n}={A\left(?\right)}^{-1}\left[-\frac{?A\left(?\right)}{?c_n}P\left(?\right)\right]={A\left(?\right)}^{-1}{F}^{\left(n\right)}\left(?\right)& \left(10\right)\end{array}$

[0177] where ?A(?)/?c.sub.n represents an N?N partial derivative impedance matrix, which is a result of obtaining the partial derivative of the nth sound velocity parameter c.sub.n by the impedance matrix.

[0178] The form of the formula (10) is similar to that of the formula (4). F.sup.(n)(?) may be regarded as the nth virtual source, which is an N?M matrix, written as:

[00053]$\begin{array}{cc}{F}^{\left(n\right)}\left(?\right)=-\frac{?A\left(?\right)}{?c_n}P\left(?\right)& \left(11\right)\end{array}$

[0179] Therefore, ?P(?)/?c.sub.n is an N?M matrix, which is a solution of a forward problem.

[0180] Then, the Frechet partial derivative matrix at frequency ?.sub.j may be written as:

[00054]$\begin{array}{cc}J\left({?}_l\right)=W\left[\frac{?P\left({?}_j\right)}{?c_1},.Math.,\frac{?P\left({?}_j\right)}{?c_n}\right]={\mathrm{WA}\left({?}_j\right)}^{-1}\left[F^{\left(1\right)}\left({?}_j\right),.Math.,F^{\left(N\right)}\left({?}_j\right)\right],& \left(12\right)\end{array}$

[0181] where J(?.sub.j) represents an M?MN matrix, and the whole Frechet partial derivative matrix J may be obtained by combining the J(?.sub.j) corresponding to all the discrete frequencies ?.sub.j.

[0182] A total virtual source matrix may be defined as F(?)=[F.sup.(1)(?), . . . , F.sup.(N)(?)], which is an N?MN matrix.

[0183] Explicitly calculating the partial derivative matrix J is highly time-consuming, and the gradient matrix may be directly calculated by substituting the formula (12) into the formula (8):

[00055]$\begin{array}{cc}{?}_{c}E\left(c,{?}_j\right)={F\left({?}_j\right)}^T\left[{A\left({?}_j\right)}^{-1}{W}^T?D\left(c,{?}_j\right)\right]={F\left({?}_j\right)}^{T}V\left({?}_j\right)& \left(13\right)\end{array}$

[0184] where ?D(c,?.sub.j) represents an M?M residual data matrix at frequency ?.sub.j; and ?.sub.cE(c,?.sub.j) represents an MN?M gradient matrix at frequency ?.sub.j.

[0185] Residual wave field V(?.sub.j) is an N?M matrix, which is calculated by the following formula:

[00056]$\begin{array}{cc}V\left({?}_j\right)={A\left({?}_j\right)}^{-1}{W}^T?D\left(c,{?}_j\right)& \left(14\right)\end{array}$

[0186] The nth element of the gradient matrix ?.sub.cE(c) is calculated by the following formula:

[00057]$\begin{array}{cc}{?}_c{E\left(c\right)}_n={.Math.}_{j=1}^{N_{?}}{.Math.}_{s=1}^M{.Math.}_{r=1}^M{?}_c{E\left(c,{?}_j\right)}_{\left(s-1\right)N+n,r}& \left(15\right)\end{array}$

[0187] where the subscript of ?.sub.cE(c,?.sub.j).sub.(s-1)N+n,r represents an element in row (s?1)N+n and column r.

[0188] In conclusion, to calculate a gradient component of each frequency ?.sub.j, forward modeling may be performed for 2M times, where the forward modeling is performed for M times to obtain a forward wave field P(?) according to the formula (4), and performed for another M times to obtain the residual wave field V(?), as shown in formula (14). Therefore, the calculation quantity required by the gradient matrix ?.sub.cE(c) is acceptable, and the efficiency can be further improved by a parallel calculation method.

[0189] Most importantly, the wave fields generated by all the nodes as single scattering points are collected by the Frechet partial derivative matrix J at a receiving node, and the residual matrix ?D(c) is real data received by a sensor from an unknown scattering point. Therefore, as shown in the formula (8), the gradient matrix is actually a result of zero-lag correlation between a corresponding scattering signal at each node and a real scattering signal. Naturally, a correlation value is large on real scattering point, but small on non-scattering points, and the non-scattering points occupy most of the measurement area. Therefore, the gradient matrix may be used as an image of a reflection point, and the formula (8) may be used as an imaging condition.

1.3 Search-Vector Imaging Condition

[0190] The negative gradient direction may be directly used to minimize the objective function (7). However, to increase a search speed, Newton method is usually used to correct a gradient in the full waveform inversion, denoted as:

[00058]$\begin{array}{cc}\mathrm{dc}=-H^{-1}{?}_{c}E\left(c\right)& \left(16\right)\end{array}$

[0191] where dc represents the search vector; H represents an N?N Hessian matrix; and an element of H is obtained by formula (17):

[00059]$\begin{array}{cc}{H}_{\mathrm{pq}}=\frac{{?}^{2}E\left(c\right)}{?c_p?c_q},p=\left(1,2,.Math.,N\right);q=\left(1,2,.Math.,N\right)& \left(17\right)\end{array}$

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

[00060]$\begin{array}{cc}H?H_a=J^T{J}^{*}& \left(18\right)\end{array}$

[0193] As can be seen from the formulas (9) and (18), a main diagonal element of H.sub.a is an auto-correlative zero-lag value of a partial derivative wave field on the receiving node, and non-diagonal elements are correlated zero-lag values of two partial derivative wave fields on the receiving node. In a high-frequency limit, H.sub.a is diagonally dominant.

[0194] Therefore, the diagonal elements of H.sub.a are used to approximate the Hessian matrix to avoid inverse operation of an unrealistic huge N-dimension matrix, and the search-vector imaging condition is written as:

[00061]$\begin{array}{cc}\mathrm{Im}=-{\left[\mathrm{diag}\left(H_a\right)-?I\right]}^{-1}{?}_{c}E\left(c\right)=-{\left[\mathrm{diag}\left(H_a\right)+?I\right]}^{-1}{J}^T?D\left(c\right)& \left(19\right)\end{array}$

[0195] where Im represents an imaging result; diag(H.sub.a) represents diagonalizing H.sub.a; I represents an N-dimension unit matrix; and ? represents a damping coefficient for stable inversion.

[0196] As described in section 1.2, the gradient matrix may be directly used as the image of the defect.

[0197] According to the formula (8), the gradient matrix depends on the Frechet partial derivative J. The partial derivative is composed of all partial derivative wave fields at the receiving node. Therefore, the characteristics of the partial derivative matrix J itself may affect the imaging result. On the one hand, the energy of the partial derivative wave field corresponding to a node at a deep position is always smaller than that at a shallow position, and therefore, the partial derivative matrix J may lead to insufficient illumination at a lower portion of an imageable area. On the other hand, different relative positions of a node to a sound source node and the receiving node may lead to different energy of the partial derivative wave fields, and therefore, the partial derivative matrix J may also lead to nonuniform illumination of the imageable area in the horizontal direction. The introduction of the Hessian matrix may eliminate the insufficient and nonuniform illumination caused by the Frechet partial derivative matrix J at gradients, enable a clearer and more concentrated imaging result, and balance imaging amplitudes at various positions.

[0198] According to the formula (18), the partial derivative matrix J needs to be solved to calculate H.sub.a. However, it is unrealistic to directly calculate the partial derivative matrix 1 by the formula (12) because too much forward calculation is required.

[0199] Based on the definition of the partial derivative matrix J, the formula (9) may be further written as:

[00062]$\begin{array}{cc}{J}_{l,n}=\frac{?{\widehat{D}}_l}{?c_n}=O_{\left(x_r\right)}^T{A\left({?}_j\right)}^{-1}{f}_{\left(l,n\right)}& \left(20\right)\end{array}$

[0200] where f.sub.(l,n) represents a virtual source vector corresponding to the lth data and the nth parameter; and.

[0201] O.sub.(x.sub.r.sub.) is an N?1 vector, which is 1 only at node (x.sub.r,0) and 0 at other nodes. Therefore, the physical significance of J.sub.l,n is a received signal at node f.sub.(l,n) corresponding to the sound field excited by the virtual source vector (x.sub.r,0). In effect, f.sub.(l,n) represents the corresponding sound source at node n.

[0202] According to the reciprocal theorem, the wave field signal emitted by the nth node and received by node (x.sub.r,0) is equal to the wave field signal emitted by node (x.sub.r,0) and received by the nth node.

[0203] Therefore, J.sub.l,n may be obtained by the following formula:

[00063]$\begin{array}{cc}{J}_{l,n}=f_{\left(l,n\right)}^T{A\left({?}_j\right)}^{-1}{O}_{\left(x_r\right)}=f_{\left(l,n\right)}^T{\tilde{p}}_{\left({?}_j,x_r\right)}& \left(21\right)\end{array}$

[0204] where {tilde over (p)}.sub.(?.sub.j.sub.,x.sub.r.sub.)=A(?.sub.j).sup.?1O.sub.(x.sub.r.sub.); and

[0205] for S(x.sub.r,x.sub.s,?.sub.j), there is an excitation value a(?.sub.j) only at x.sub.r=x.sub.s, and in this case, a(?.sub.j)=S(x.sub.s,x.sub.s,?.sub.j); and when x.sub.r?x.sub.s, S(x.sub.r,x.sub.s,?.sub.j)=0.

[00064]$\begin{array}{c}{p}_{\left({?}_j,x_r\right)}=a\left({?}_j\right){A\left({?}_j\right)}^{-1}{O}_{\left(x_r\right)}\\ P\left({?}_j\right)=\left[p_{\left({?}_j,x_1\right)},.Math.,p_{\left({?}_j,x_M\right)}\right]\\ P\left(?\right)={A\left(?\right)}^{-1}S\left(?\right)=a\left({?}_j\right){A\left({?}_j\right)}^{-1}\left[O_{\left(x_1\right)},.Math.,O_{\left(x_M\right)}\right]\end{array}$

[0206] The following formula is then derived: {tilde over (P)}(?.sub.j)=[{tilde over (p)}.sub.(?.sub.j.sub.,x.sub.1.sub.), . . . , {tilde over (p)}.sub.(?.sub.j.sub.,x.sub.M.sub.)]

[0207] The partial derivative matrix of certain data (the lth data) relative to all sound velocity parameters may be directly calculated by the following formula:

[00065]$\begin{array}{cc}{J}_l=\frac{?{\hat{D}}_l}{?c}={\left[f_{\left(l,1\right)},f_{\left(l,2\right)},.Math.,f_{\left(l,N\right)}\right]}^T{\tilde{p}}_{\left({?}_j,x_r\right)}& \left(22\right)\end{array}$

[0208] where J.sub.l represents the N?1 partial derivative matrix of the lth data, and

[0209] f.sub.(l,n) represents a virtual source vector corresponding to the lth data and the nth parameter and is obtained by the total virtual source matrix, where the total virtual source matrix is expressed as: F(?)=[F.sup.(1)(?), . . . , F.sup.(N)(?)].

[0210] For frequency ?.sub.j, the following can be derived:

[00066]$\left[f_{\left(l,n\right)},f_{\left(l+1,n\right)},.Math.,f_{\left(l+M^2-1,n\right)}\right]=\left[F^{\left(n\right)}\left({?}_j\right),F^{\left(n\right)}\left({?}_j\right),.Math.,F^{\left(n\right)}\left({?}_j\right)\right]$

[0211] In the above equation, the virtual source matrix is repeated for 64 times in a row direction on the right of the equation.

[0212] Thus, to obtain the whole partial derivative matrix of a piece of data, forward calculation only needs to be performed once, and the calculation complexity is reduced by 1/N, rendering the calculation of the partial derivative matrix J feasible.

[0213] N diagonal elements of H.sub.a are calculated by the following formula:

[00067]$\begin{array}{cc}\mathrm{diag}\left(H_a\right)={.Math.}_{j=1}^{N_{?}}{.Math.}_{s=1}^M{.Math.}_{r=1}^M{J}_l?J_l& \left(23\right)\end{array}$

[0214] where ? represents Hadamard multiplication, l=(j?1)M.sup.2+(s?1)M+r.

[0215] A time-zero imaging condition is a traditional imaging condition usually for self-transmitted and self-received signals, which involves an assumption that a reflected signal is actively emitted by a reflection point at time 0. Therefore, an 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 an excitation-time imaging time, where the excitation time is a time when an incident wave reaches a reflector. Conventional imaging methods, including total focusing method (TFM), phase shift migration (PSM), and reverse time migration (RTM), use the excitation-time imaging time to image defects. The search-vector imaging condition is constructed from the perspective of full waveform inversion, as shown in the formula (19). Such an imaging condition is strictly derived from a formula and there is no assumption regarding the wave field. Therefore, the theory of such an imaging condition is stricter, and higher imaging precision can be guaranteed.

Testing Experiment

[0216] The following experiment is conducted to verify the imaging effect of the internal defect of the complex-shaped component by the imaging method in the present embodiment.

[0217] The FMC data is collected by using 64/64 OEM-PA (AOS. limited company, the United States) with a sampling frequency of 50 MHz and a time range of 80 ?s. An aluminum component with a circular surface is used as a measurement object. As shown in FIG. 4, crack and pore defects are measurement objects. A 64-element array transducer (Guangdong Gaohua Corporation Limited, China) with a center frequency of 2.5 MHz and an interval of 0.75 mm is used as an ultrasonic transceiver in the experiment. The components are immersed in water. Since the phased array is not an immersion probe, a wedge block made of acrylonitrile-butadiene-styrene (ABS) is used to protect the probe from water, and the sound velocity of the wedge block is 2200 m/s. The bottom of the wedge is a concave surface so that more ultrasonic waves are concentrated in the measured component. Water, as a coupling agent, fills the gap between the wedge block and the component but does not reach the top surface of the wedge block.

[0218] Experimental data (FMC data) is processed by using TFM, FD-RTM, GDM, and the imaging method of the present embodiment separately. The frequency band (frequency range) selected by these methods is from 1 MHz to 4 MHz. The parameter ? of the imaging method of the present embodiment is chosen to be 0.05. The overall imaging results of the four methods are as shown in FIG. 5A to FIG. 5D, and defects are easy to find in the results of TFM and the imaging method of the present embodiment in FIG. 5A and FIG. 5D. Due to insufficient illumination, the defect imaging results can be hardly found in the results of FD-RTM and GDM in FIG. 5B and FIG. 5C.

[0219] Moreover, FIG. 5E to FIG. 5H show the partially enlarged images of the imaging results of the above four imaging methods. Here, it is still difficult to distinguish between the defect images in the results of FD-RTM and GMD in FIG. 5G and FIG. 5F, and the contrast of the defect in the imaging result of GMD in FIG. 5E is higher than that of the imaging result of FD-RTM in FIG. 5F. However, in FIG. 5E, the artifact around the defect is very serious, and cracks are slightly distorted. The imaging result of the imaging method of the present embodiment in FIG. 5H indicates that the crack form can be correctly imaged, and background artifacts can be well inhibited.