CN114858926A - Full matrix data high-quality imaging method based on frequency domain reverse time migration - Google Patents
Full matrix data high-quality imaging method based on frequency domain reverse time migration Download PDFInfo
- Publication number
- CN114858926A CN114858926A CN202210504512.6A CN202210504512A CN114858926A CN 114858926 A CN114858926 A CN 114858926A CN 202210504512 A CN202210504512 A CN 202210504512A CN 114858926 A CN114858926 A CN 114858926A
- Authority
- CN
- China
- Prior art keywords
- matrix
- imaging
- excitation
- sound source
- full
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
- G01N29/04—Analysing solids
- G01N29/06—Visualisation of the interior, e.g. acoustic microscopy
- G01N29/0654—Imaging
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
- G01N29/04—Analysing solids
- G01N29/06—Visualisation of the interior, e.g. acoustic microscopy
- G01N29/0654—Imaging
- G01N29/069—Defect imaging, localisation and sizing using, e.g. time of flight diffraction [TOFD], synthetic aperture focusing technique [SAFT], Amplituden-Laufzeit-Ortskurven [ALOK] technique
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
- G01N29/04—Analysing solids
- G01N29/07—Analysing solids by measuring propagation velocity or propagation time of acoustic waves
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
- G01N29/44—Processing the detected response signal, e.g. electronic circuits specially adapted therefor
- G01N29/4445—Classification of defects
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2291/00—Indexing codes associated with group G01N29/00
- G01N2291/01—Indexing codes associated with the measuring variable
- G01N2291/011—Velocity or travel time
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2291/00—Indexing codes associated with group G01N29/00
- G01N2291/02—Indexing codes associated with the analysed material
- G01N2291/023—Solids
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N2291/00—Indexing codes associated with group G01N29/00
- G01N2291/02—Indexing codes associated with the analysed material
- G01N2291/028—Material parameters
- G01N2291/0289—Internal structure, e.g. defects, grain size, texture
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Optimization (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Life Sciences & Earth Sciences (AREA)
- Health & Medical Sciences (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Analysis (AREA)
- Analytical Chemistry (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Acoustics & Sound (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Signal Processing (AREA)
- Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
Abstract
The invention provides a full matrix data high-quality imaging method based on frequency domain reverse time migration, which comprises the steps of respectively preprocessing an excitation signal, sound velocity distribution and full matrix data to obtain an excitation sound source matrix, an impedance matrix and a receiving sound source matrix; and taking the preprocessing result as input, carrying out LU decomposition processing on the impedance matrix, solving to obtain a transmitting wave field and a receiving wave field, and finally obtaining an imaging result by using imaging conditions. The imaging method utilizes the characteristic that the impedance matrix is fixed under the same frequency, and can solve the sound pressure field generated by multiple times of excitation under the same frequency at one time through LU decomposition, thereby greatly saving the calculation time and improving the calculation efficiency. The wave field forward modeling is carried out on the basis of the full-wave equation, the simplification of wave field propagation is less, the method is suitable for measuring defects in any complex sound velocity distribution, and the high quality of an imaging result can be ensured.
Description
Technical Field
The invention belongs to the technical field of detection methods, and particularly relates to a full matrix data high-quality imaging method based on frequency domain reverse time migration.
Background
The curved surface part is widely applied to various fields. Among them, shafts and camshafts are core elements in drive systems, and in addition, more complex curved parts are widely present in the fields of automobiles, energy sources, aerospace, and the like due to dynamic requirements. However, due to the complex manufacturing process and the harsh service environment of such components, the components are prone to void and crack defects, which greatly affect the mechanical properties of the components, and can lead to premature failure of the components or even serious accidents. Therefore, high quality inspection of internal defects of curved parts is an urgent issue, but the presence of larger part sizes and complex surfaces also presents greater challenges to inspection techniques.
Phased array ultrasonic imaging obtains outstanding status in nondestructive testing, and can effectively represent invisible defects inside components. Full matrix acquisition (FMC) is an important acquisition mode in phased array ultrasound, where each array element acts as a receive and a transmitter, which captures the complete time domain signal of each pair of transmit and receive array elements. Existing imaging methods for full matrix data have shown their ability to detect parts with flat surfaces that are parallel or oblique to the surface of the phased array, but these methods cannot be extended directly to curved surfaces. In order to measure curved surface parts, scholars have introduced flexible array sensors to accommodate curved surfaces, but these techniques increase the cost of the measurement. On the other hand, non-contact measurement by water immersion is more widely adopted, which is mainly handled in the invention with respect to the full matrix data obtained in this case.
The water immersion measurement mode puts higher requirements on post-processing of full matrix data. Several methods of ultrasound imaging for curved interface structures exist. Typically, the full focus method (TFM) adapts to curved interfaces by introducing a ray tracing method. Furthermore, Sutcliffe proposes a virtual source method for high resolution imaging of defects in bilayer media with complex interfaces, where the point of incidence of the ray on the reflecting interface is still searched using the fermat principle. However, this approach does not take into account multiple echoes from the defect, which reduces the resolution of the imaging results and leads to artifacts in the results. The TFM method is a standard method in the full matrix imaging method, and it has been shown in previous studies to be efficient, but its efficiency is greatly reduced when the sound velocity model of the measurement region becomes complicated. In the above studies, the sound velocity models were all bi-level to ensure imaging efficiency, however, the computational complexity of ray tracing grows exponentially with the number of layers in the sound velocity model, and it becomes unacceptable when the number of layers increases to four layers.
Frequency domain wavefield extrapolation may also be used in the case of curved interfaces. Lukomski extends the Phase-shift-plus-interpolation (PSPI) method to the full matrix data to measure curved parts. In the PSPI method, the acoustic velocity field of the interface region is approximated as a combination of two uniform acoustic velocity fields, phase shift shifts are respectively performed according to the two acoustic velocities, and then the corrected extrapolated wavefield is superimposed. Yang et al considered multiple reflections from the bottom of the measurement part to improve the previous imaging method. Furthermore, Chang et al introduced a non-stationary phase shift (NSPS) method to compensate for the effect of horizontal sound velocity variations on the downward extrapolation, which can eliminate discontinuities in the reconstructed wavefield compared to PSPI. However, frequency domain wavefield extrapolation, which reduces the wave equation to a one-way form and approximates the acoustic velocity model, significantly degrades the imaging quality as the number of curved interfaces increases.
Reverse Time Migration (RTM) was first developed in the field of seismic imaging to reconstruct the wavefield of the measurement volume using full wave equations. The RTM has the advantages of high imaging precision, no inclination angle limitation, strong adaptability to any complex speed model and the like. In recent years, with the rapid development of computer technology, the huge calculation cost for restricting the development is greatly relieved, and the FMC ultrasonic imaging is introduced in partial scenes. Rao et al used Elastic Reverse Time Migration (ERTM) to image irregularly shaped grooves in a metal block. Liu et al successfully introduced ERTM into ultrasonic imaging of concrete structures. However, ERTM is calculated in the time domain and requires iteration at each sample time. Furthermore, for FMC data, forward and backward propagation of the excitation and received signals, respectively, needs to be achieved for each excitation. These make the method computationally complex and computationally long.
Disclosure of Invention
In order to solve the problems in the prior art, the invention provides a frequency domain reverse time migration imaging method for ultrasonic full matrix data, which can greatly improve the imaging quality of a curved surface part while controlling the calculation time. The invention extrapolates the wave field in the selected frequency range, and the wave field with different excitations under a certain frequency can be calculated by one iteration, thereby effectively reducing the calculation complexity compared with ERTM.
A full matrix data high-quality imaging method based on frequency domain reverse time migration comprises the following steps:
(1) inputting an excitation signal, sound velocity distribution and full matrix data;
(2) respectively preprocessing the excitation signal and the full matrix data to respectively obtain an excitation sound source matrix and a receiving sound source matrix;
(3) discretizing the sound velocity distribution, and calculating to obtain an impedance matrix;
(4) according to the obtained excitation sound source matrix and the received sound source matrix, adopting LU decomposition of an impedance matrix to respectively solve and obtain a transmitting wave field and a receiving wave field;
(5) and processing the transmitting wave field and the receiving wave field by adopting imaging conditions to obtain an initial image, and performing post-processing on the initial image to obtain an imaging result.
Wherein the excitation signal in step (1) is represented as E (x) r ,x s T), full matrix data is represented as D (x) r ,x s T), the sound velocity distribution is denoted as c (x, z);
in step (2), exciting the sound source momentThe matrix and the received sound source matrix respectively comprise a selected frequency band range [ omega ] min ,ω max ]Excitation sound source matrix [ E ] corresponding to each discrete frequency omega 1 (ω),E 2 (ω),…,E n (ω)]And receiving the sound source matrix
Wherein n represents the excitation times (i.e. the number of the vibration elements of the phased array probe, and one vibration element is excited once). Selecting a band range [ omega ] min ,ω max ]Selecting according to the characteristics of the excitation signal and the receiving signal (full matrix data); within the selected frequency band, the frequency domain discrete step size between every two adjacent discrete frequencies is determined by the sampling frequency, and is denoted as d ω.
Preferably, the full matrix data is preprocessed as follows:
the method comprises the steps of firstly filtering direct waves by applying a window function to full matrix data, carrying out Fourier transform on the full matrix data with the direct waves filtered out in a time dimension, then carrying out complex conjugation on data in a frequency domain obtained by the transform to realize an inverse time effect, and sorting to obtain a receiving sound source matrix.
Full matrix acquisition (FMC) acquires time domain signals between all transmitting-receiving pairs in a phased array probe, and for an n-element phased array probe, FMC data comprises n 2 An A-scan signal. D (x) r ,x s T) represents the transmission from the s-th element of the signal received by the r-th element, as shown in FIG. 2, of ultrasound from the element location (x) s 0) is emitted, reflected by a scattering point located at (x, z) and then by a scattering point located at (x, z) r And, 0) sensor reception. The fourier transform in the time domain yields a frequency domain representation of the full matrix data, as shown in the following equation:
D(x r ,x s ,ω)＝∫D(x r ,x s ,t)e jkωt dt
wherein j represents an imaginary unit; k represents a horizontal wave number; ω represents frequency; t represents time.
From the property of fourier transform, the fourier transform of the inverse time signal is equivalent to the original signal taking the conjugate, so the inverse time signal in the frequency domain can be written as follows:
wherein denotes a complex conjugate.
Preferably, the excitation signal is pre-processed as follows:
and carrying out Fourier transform processing on the excitation signal in a time domain to obtain an excitation sound source matrix.
Since the calculation region is limited, if unconditional conditions are not added, false reflection can be generated at the boundary of the calculation region, and preferably, when the sound velocity distribution is subjected to discretization processing, the impedance matrix is optimized by adding a perfect matching layer around the discretization model.
Therein, the attenuation γ (x) of the Perfect Matching Layer (PML) is defined as follows:
wherein d is the width of the PML layer; x' is a coordinate in a local coordinate system in the PML, and the origin of the local coordinate system is positioned outside the boundary of the model; c. C PML Is a defined scalar.
Preferably, in step (3), the calculation formula of the impedance matrix a (ω) is as follows:
wherein the content of the first and second substances,
representing a Laplace operator obtained after discretization; ω represents frequency; and c represents the speed of sound.
Preferably, the sound velocity distribution is discretized by a 9-point finite difference method.
Preferably, the specific steps of step (4) are as follows:
4.1, taking any frequency which is not traversed in the selected frequency band range as a current frequency, taking an excitation sound source matrix, a receiving sound source matrix and an impedance matrix corresponding to the current frequency as input, carrying out LU decomposition on the impedance matrix, and solving to obtain an excitation sound pressure field and a receiving sound pressure field corresponding to the current frequency;
4.2 go through all frequencies in the selected frequency band range, resulting in the transmit wavefield and the receive wavefield.
As a further preference, in step 4.1, the calculation formula for solving the sound pressure field generated by the sound source matrix is as follows:
L(ω)U(ω)[P 1 (ω),P 2 (ω),…,P n (ω)]＝[S 1 (ω),S 2 (ω),…,S n (ω)]
wherein n represents the number of excitations; l (omega) and U (omega) represent two factors obtained after LU decomposition of the impedance matrix A (omega); [ S ] 1 (ω),S 2 (ω),…,S n (ω)]Representing either an excitation source matrix or a reception source matrix; [ P ] 1 (ω),P 2 (ω),…,P n (ω)]Representing either an excitation sound pressure field or a received sound pressure field.
The derivation process of the above formula for calculating the sound pressure field is as follows:
the two-dimensional wave equation in the frequency domain is:
wherein x and z are two-dimensional space coordinates respectively; omega is frequency; p (x, z, ω) is the sound pressure field; s (x, z, ω) is the sound source field; c (x, z) is the sound velocity distribution.
Sound velocity distribution as shown in fig. 2(a), the 9-point finite difference method can introduce formula 1) to discretize the spatial domain, where the sound velocity is discretized into l ═ N x ×N y Grid points, as shown in fig. 2 (b); wherein N is x And N y The number of grid points on the x and y axes, respectively. The discrete acoustic wave equations can be written in matrix form as follows:
A(ω)P(ω)＝S(ω) 2)
wherein A (ω) is an impedance matrix; p (omega) is a discrete sound pressure field(ii) a S (ω) is a discrete sound source vector. A (ω) can also be interpreted as a forward propagation operator, as:
here, P (ω) and S (ω) are l × 1 column vectors, and a (ω) is an l × l matrix.
For forward propagation, S (ω) is given by the transmitted wave and a (ω) depends on the frequency (ω) and the acoustic properties of the medium, so only P (ω) is unknown and can be solved by equation (2). The LU decomposition of the impedance matrix a (ω) can be used to solve P (ω) in equation (2), and the factors of the LU decomposition of the impedance matrix a (ω) corresponding to the same frequency can be repeatedly used to solve the sound pressure field generated by any new sound source vector S (ω) corresponding to the frequency. Therefore, the sound pressure field excited for multiple times under the same frequency can be efficiently calculated at one time, and the forward propagation solving mode is very suitable for full matrix data.
For the case of multiple excitations at the same frequency, the formula 3 can be arranged
L(ω)U(ω)[P 1 (ω),P 2 (ω),…,P n (ω)]＝[S 1 (ω),S 2 (ω),…,S n (ω)] 3)
More preferably, the frequency is traversed in descending order.
The frequency band range [ omega ] is selected in the order of small to large min ,ω max ]Traversing the internal frequency, and judging whether the current frequency omega is less than omega or not after the traversal of the current frequency omega is finished max (ii) a If yes, the next frequency (omega + d omega) is traversed until the current judgment frequency is equal to omega max And outputting the solved transmitting wave field S when all the frequencies are traversed i (x, z, ω) and the received wavefield R i (x,z,ω)。
Preferably, in step (5), the expression of the imaging conditions is:
wherein I (x, z) represents an initial image; n represents the number of excitations; n is a radical of ω Representing the total number of discrete frequencies in the selected frequency band; s i (x, z, ω) represents the transmit wavefield; r i (x, z, ω) represents the receive wavefield; denotes taking the complex conjugate.
Preferably, the initial image is post-processed as follows:
firstly, carrying out sharpening processing on an initial image by adopting Laplace filtering, then carrying out Hilbert transform on the sharpened initial image to obtain the envelope of the image, and finally filtering salt and pepper noise by adopting median filtering to obtain the imaging result.
The full matrix data high-quality imaging method based on the frequency domain reverse time migration can be used for carrying out high-quality imaging on the internal defects of the curved surface part. The method can be summarized in three steps: forward extrapolation of the acoustic wavefield, backward extrapolation of the received wavefield, and implementation of the imaging conditions. The input of the imaging method comprises an excitation signal, sound velocity distribution and full matrix data, and the three are preprocessed and then wave field solving is carried out. E i (ω) and
respectively representing the excitation vector corresponding to the frequency omega at the i-th excitation and the receiving vector in the reverse time, i ∈ [1, n ∈](ii) a n is the number of array elements in the phased array probe and is equal to the excitation times. Because the impedance matrix A (omega) is fixed under a certain frequency, excitation vectors and reverse time receiving vectors in different excitations can be respectively arranged into a sound source matrix, namely an excited sound source matrix and a receiving sound source matrix, so sound pressure fields generated by different excitations can be simultaneously calculated in one iteration, and the calculation efficiency is greatly improved. The imaging method of the invention adopts a frequency-by-frequency iteration mode to calculate the selected frequency band range [ omega ] min ,ω max ]The corresponding sound pressure field at each discrete frequency, wherein the frequency band range is determined by the full matrix data, and the frequency domain discrete step length d omega is determined by the sampling frequency. S i (x, z, ω) and R i (x, z, ω) is the emission corresponding to the frequency ω at the i-th excitationAnd the wave field and the receiving wave field take the transmitting wave field and the receiving wave field as input, and an imaging result is obtained according to the imaging condition.
Compared with the prior art, the invention has the beneficial effects that:
the invention relates to a full matrix data high-quality imaging method based on frequency domain reverse time migration, which comprises the steps of respectively preprocessing an excitation signal, sound velocity distribution and full matrix data to obtain an excitation sound source matrix, an impedance matrix and a receiving sound source matrix; and taking the preprocessing result as input, carrying out LU decomposition processing on the impedance matrix, solving to obtain a transmitting wave field and a receiving wave field, and finally obtaining an imaging result by using imaging conditions. The imaging method utilizes the characteristic that the impedance matrix is fixed under the same frequency, and can solve the sound pressure field generated by multiple times of excitation under the same frequency at one time through LU decomposition, thereby greatly saving the calculation time and improving the calculation efficiency. The wave field forward modeling is carried out on the basis of a full-wave equation (a two-dimensional wave equation in a frequency domain), the simplification of wave field propagation is less, the method is suitable for measuring defects in any complex sound velocity distribution, and the high quality of an imaging result can be ensured.
Drawings
FIG. 1 is a flow chart of an imaging method of an embodiment of the present invention;
fig. 2(a) is a schematic diagram of a sound velocity distribution; (b) the sound velocity distribution is a schematic diagram after discretization processing;
FIG. 3 (a) is a flow chart of preprocessing of full matrix data; (b) a flow chart for obtaining imaging results from the transmit and receive wavefields;
FIG. 4 (a) is a schematic diagram illustrating the detection of hole defects in copper parts; (b) the detection diagram of the hole defect in the aluminum part is shown; (c) the detection of crack defects in the aluminum parts is shown schematically;
FIG. 5 (a) shows the TFM-based imaging of hole defects in copper parts; (b) the imaging result of the hole defect in the copper part by the PSM-based method is shown; (c) the result of FD-RTM imaging of the hole defect in the copper part is obtained; (d) the method is an imaging result of a TFM-based method on the hole defect in the aluminum part; (e) the imaging result of the hole defect in the aluminum part by the PSM-based method is shown; (f) the method is an imaging result of FD-RTM on the hole defect in the aluminum part; (g) the method is an imaging result of the TFM-based method on the crack defects in the aluminum parts; (h) the method is an imaging result of the PSM-based method on the crack defects in the aluminum part; (i) the method is an imaging result of FD-RTM on crack defects in the aluminum part;
FIG. 6 (a) is a graph of pixel amplitude on a labeled line in the imaging results of hole defects in copper parts by different imaging methods; (b) marking a pixel amplitude curve on a straight line for the imaging results of the hole defects in the aluminum parts by different imaging methods; (c) and marking a pixel amplitude curve on a straight line for the imaging results of the crack defects in the aluminum parts by different imaging methods.
Detailed Description
As shown in fig. 1, a full matrix data high quality imaging method based on frequency domain reverse time migration includes:
1. for input excitation signal E (x) r ,x s T), full matrix data D (x) r ,x s T) and the sound velocity distribution c (x, z) are respectively preprocessed to respectively obtain an excitation sound source matrix, a receiving sound source matrix and an impedance matrix.
1.1 Pre-processing of excitation signals
And carrying out Fourier transform processing on the excitation signal in a time domain to obtain an excitation sound source matrix.
1.2 Pre-processing of full matrix data
As shown in fig. 3 (a), a window function is first applied to the full-matrix data to filter out direct waves, the full-matrix data with the filtered direct waves is fourier-transformed in the time dimension, and then the data in the frequency domain obtained by the transformation is complex-conjugated to realize the reverse time effect, and the received sound source matrix is obtained by sorting.
Full matrix acquisition (FMC) acquires time domain signals between all transmitting-receiving pairs in a phased array probe, and for an n-element phased array probe, the FMC data comprises n 2 An A-scan signal. D (x) r ,x s T) represents the transmission from the s-th element of the signal received by the r-th element, as shown in FIG. 2, of ultrasound from the element location (x) s 0) emitted, reflected by a scattering point located at (x, z)In (x) r And, 0) sensor reception. The fourier transform in the time domain yields a frequency domain representation of the full matrix data, as shown in the following equation:
D(x r ,x s ,ω)＝∫D(x r ,x s ,t)e jkωt dt
wherein j represents an imaginary unit; k represents a horizontal wave number; ω represents frequency; t represents time.
From the property of fourier transform, the fourier transform of the inverse time signal is equivalent to the original signal taking the conjugate, so the inverse time signal in the frequency domain can be written as follows:
wherein denotes a complex conjugate.
The obtained excitation sound source matrix and the received sound source matrix respectively comprise a selective frequency band range [ omega ] min ,ω max ]Excitation sound source matrix [ E ] corresponding to each discrete frequency omega 1 (ω),E 2 (ω),…,E n (ω)]And receiving the sound source matrix
Where the subscript n indicates the number of firings. Selecting a band range [ omega ] min ,ω max ]Selecting according to the characteristics (full matrix data) of the excitation signal and the receiving signal; within the selected frequency band, the frequency domain discrete step size between every two adjacent discrete frequencies is determined by the sampling frequency, and is denoted as d ω.
1.3 preprocessing of the sound velocity distribution
Discretizing sound velocity distribution by adopting a 9-point finite difference method, and calculating by utilizing the following formula to obtain an impedance matrix A (omega):
wherein the content of the first and second substances,
representing a Laplace operator obtained after discretization; ω represents frequency; and c represents the speed of sound.
Because the calculation area is limited, false reflection can be generated at the boundary of the calculation area if no condition is added, and when discretization processing is carried out on sound velocity distribution, the impedance matrix is optimized by adding perfect matching layers around the discretization model.
Therein, the attenuation γ (x) of the Perfect Matching Layer (PML) is defined as follows:
where d is the width of the PML layer, x' is the coordinates in the local coordinate system in the PML, the origin of the local coordinate system is located outside the model boundary, c PML Is a defined scalar.
2. And respectively solving to obtain a transmitting wave field and a receiving wave field by adopting LU decomposition of an impedance matrix according to the obtained excitation sound source matrix and the receiving sound source matrix.
The method comprises the following specific steps:
2.1 selecting the frequency band range [ omega ] min ,ω max ]Taking any one of the unexploded frequencies omega as the current frequency, and using the excitation sound source matrix [ E ] corresponding to the current frequency 1 (ω),E 2 (ω),…,E n (ω)]Receiving sound source matrix
Taking the impedance matrix A (omega) as input, carrying out LU decomposition on the impedance matrix A (omega), and solving to obtain an excitation sound pressure field and a receiving sound pressure field corresponding to the current frequency;
the derivation process of the solving formula of the sound pressure field is as follows:
the two-dimensional wave equation in the frequency domain is:
wherein x and z are two-dimensional space coordinates respectively; omega is frequency; p (x, z, ω) is the sound pressure field; s (x, z, ω) is the sound source field; c (x, z) is the sound velocity distribution.
Sound velocity distribution as shown in fig. 2(a), the 9-point finite difference method can introduce formula 1) to discretize the spatial domain, where the sound velocity is discretized into l ═ N x ×N y Grid points, as shown in fig. 2 (b); wherein N is x And N y The number of grid points on the x and y axes, respectively. The discrete acoustic wave equations can be written in matrix form as follows:
A(ω)P(ω)＝S(ω) 2)
wherein A (ω) is an impedance matrix; p (omega) is a discrete sound pressure field; s (ω) is a discrete sound source vector. A (ω) can also be interpreted as a forward propagation operator, as:
here, P (ω) and S (ω) are l × 1 column vectors, and a (ω) is an l × l matrix.
For forward propagation, S (ω) is given by the transmitted wave and a (ω) depends on the frequency (ω) and the medium' S ascending characteristics, so only P (ω) is unknown and can be solved by equation (2). P (ω) in equation (2) can be solved using LU decomposition of the impedance matrix a (ω); since the impedance matrix a (ω) corresponding to the same frequency is fixed, the factors of LU decomposition of the impedance matrix a (ω) corresponding to the same frequency ω can be repeatedly used to solve the sound pressure field generated by any new sound source vector S (ω) corresponding to the frequency ω. Therefore, the sound pressure field excited for multiple times under the same frequency can be efficiently calculated at one time, and the forward propagation solving mode is very suitable for full matrix data.
For the case of multiple excitations at the same frequency, the formula 3 can be arranged
L(ω)U(ω)[P 1 (ω),P 2 (ω),…,P n (ω)]＝[S 1 (ω),S 2 (ω),…,S n (ω)] 3)
Wherein n represents the excitation frequency (namely the number of the phased array probe vibration elements); l (omega) and U (omega) represent two obtained by LU decomposition of impedance matrix A (omega)A factor; [ S ] 1 (ω),S 2 (ω),…,S n (ω)]Representing either an excitation source matrix or a reception source matrix; [ P ] 1 (ω),P 2 (ω),…,P n (ω)]Representing either an excitation sound pressure field or a received sound pressure field.
And 2.2, traversing and selecting all frequencies in the frequency band range to obtain a transmitting wave field and a receiving wave field.
The frequency band range [ omega ] is selected in the order of small to large min ,ω max ]Traversing the internal frequency, and judging whether the current frequency omega is less than omega or not after the traversal of the current frequency omega is finished max (ii) a If yes, the next frequency (omega + d omega) is traversed until the current judgment frequency is equal to omega max And outputting the solved transmitting wave field S when all the frequencies are traversed i (x, z, ω) and the received wavefield R i (x,z,ω)。
3. The transmit and receive wavefields are processed using imaging conditions, as shown in figure 3 (b), to obtain an initial image. Firstly, the initial image is sharpened by adopting Laplace filtering, then the sharpened initial image is subjected to Hilbert transform to obtain the envelope of the image, and finally the imaging result is obtained after the salt and pepper noise is filtered by adopting median filtering.
Wherein, the expression of the imaging condition of the initial image I (x, z) is:
in the above formula, n represents the number of times of excitation; n is a radical of ω Representing the total number of discrete frequencies in the selected frequency band;
S i (x, z, ω) represents the transmit wavefield; r i (x, z, ω) represents the receive wavefield; denotes taking the complex conjugate.
Detection experiment
Experiments were conducted to verify the effectiveness of the present embodiment method in imaging defects in curved surface parts. As shown in FIGS. 4 (a) and (b), first, void defects were measured in copper and aluminum parts having a radius of 80mm, and there were four void defects having a diameter of 1mm in the parts. Then, V-shaped crack defects in an aluminum part having a radius of 80mm were measured in (c) of FIG. 4. In the experiment, a 64-element phased array probe (Shantou ultrasonic) with the center frequency of 2.5MHz and the element spacing of 0.75mm is used, and FMC data acquisition is realized by matching with a data acquisition card 64/64OEM-PA (AOS. Ltd, America), wherein the sampling frequency is 50MHz, and the time range is 60 mu s. The above parts were immersed in water during measurement, and since the phased array was not a water immersion probe, a wedge made of polyphenylene sulfide (PPS) was used to isolate the probe from the water, the wedge sound velocity being 2200 m/s. The bottom surface of the wedge is a concave surface for focusing more ultrasonic waves into the measurement component. Water acts as a coupling agent filling the gap between the wedge and the component, but does not reach the top surface of the wedge.
Three imaging methods were used to process the experimental data, and the imaging results are shown in FIG. 5, wherein the comparative methods are labeled as TFM-based and PSM-based in the figure, and the method of this example is labeled as FD-RTM. The frequency band selected for all three experiments was from 1MHz to 4 MHz.
For hole defect detection in copper parts, the TFM-based method and FD-RTM method can clearly image all defects, as shown in (a) and (c) of FIG. 5, where the latter is more clear than the former. However, as shown in fig. 5 (b), the PSM-based method has poor imaging quality, and the two defects in the lower left corner are difficult to distinguish from the background.
For hole defect detection in aluminum parts, in fig. 5 (d), severe artifacts appear on the upper side of the imaging area, which severely disturb the identification of the top three defects before the TFM-based method. The situation gets worse for the PSM-based approach, where only the middle defect can be identified in fig. 5(e), with other defects buried in the background. FD-RTM suppresses artifacts very well, all defects are in sharp contrast to the background in fig. 5 (f).
The wave energy reflected by the crack defect is far larger than that of the void defect, so that the V-shaped crack can be imaged clearly by the three methods, and the background interference is low. The TFM-based method has higher crack image contrast than the PSM-based method, but both methods image sharp corners of real cracks as rounded corners. Although there are artifacts in the imaging results of FD-RTM, the imaging target is clearly identifiable and the boundary of FD-RTM is clearest, the shape being closest to the shape of the actual crack.
In addition, for further comparison, the amplitude values at the portions marked with the broken lines in (a) to (i) in fig. 5 are extracted and normalized for comparison, as shown in (a) to (c) in fig. 6. In the first two cases (hole defects in copper parts and hole defects in aluminum parts), the main lobe widths of the three imaging modes are almost the same, while the sidelobes of FD-RTM are the smallest in the second case (hole defects in aluminum parts). In addition, the relative difference of the two peaks of the FD-RTM is minimal, indicating that the FD-RTM still maintains high contrast at the edges of the imaged area. For the third case (crack defect in aluminum part), the main lobe of FD-RTM falls the fastest on both sides as shown in fig. 6 (c), which guarantees the clear boundary in fig. 5 (i).
In order to quantitatively analyze the imaging result, a dimensionless parameter API (array performance indicator) is introduced, which has the following formula
API＝A -6dB /λ 2 15)
Wherein A is -6dB Is the area with the ratio of the defect area image peak value exceeding-6 dB, and lambda is the wave number corresponding to the central frequency. Obviously, smaller APIs mean larger resolution. Table 1 lists the average API of the void defect imaging results in the first two experiments.
TABLE 1 average API for hole defects in imaging results obtained by different methods in the experiment
Since the two defects in the lower left corner of fig. 5(e) are completely indistinguishable from the background, only the average API for the other two defects is calculated. It can be seen that the API of the imaging result of FD-RTM was minimal in both experiments, and the resolution of FD-RTM method was higher than TFM-based method. In the three experiments (hole defect of copper part, hole defect of aluminum part and crack defect of aluminum part), FD-RTM images are 978s, 1019 and 1031s respectively, and TFM-based images are 1566s, 1536 s and 1554s respectively. Therefore, experiments also show that the FD-RTM method is very suitable for high-quality defect imaging of curved parts.
Claims (10)
1. A full matrix data high-quality imaging method based on frequency domain reverse time migration is characterized by comprising the following steps:
(1) inputting an excitation signal, sound velocity distribution and full matrix data;
(2) respectively preprocessing the excitation signal and the full matrix data to respectively obtain an excitation sound source matrix and a receiving sound source matrix;
(3) discretizing the sound velocity distribution, and calculating to obtain an impedance matrix;
(4) according to the obtained excitation sound source matrix and the received sound source matrix, adopting LU decomposition of an impedance matrix to respectively solve and obtain a transmitting wave field and a receiving wave field;
(5) and processing the transmitting wave field and the receiving wave field by adopting imaging conditions to obtain an initial image, and performing post-processing on the initial image to obtain an imaging result.
2. The method for high-quality imaging of full matrix data based on frequency domain reverse time migration according to claim 1, wherein the full matrix data is preprocessed as follows:
the method comprises the steps of firstly filtering direct waves by applying a window function to full matrix data, carrying out Fourier transform on the full matrix data with the direct waves filtered out in a time dimension, then carrying out complex conjugation on data in a frequency domain obtained by the transform to realize an inverse time effect, and sorting to obtain a receiving sound source matrix.
3. The method of claim 1, wherein the excitation signal is preprocessed according to the following steps:
and carrying out Fourier transform processing on the excitation signal in a time domain to obtain an excitation sound source matrix.
4. The method for high-quality imaging of the full matrix data based on the frequency domain reverse time migration as claimed in claim 1, wherein when discretizing the sound velocity distribution, the impedance matrix is optimized by adding a perfect matching layer around the discretization model.
5. The method for high-quality imaging of full matrix data based on frequency domain reverse time migration according to claim 1, wherein in step (3), the impedance matrix A (ω) is calculated as follows:
wherein the content of the first and second substances,
representing a Laplace operator obtained after discretization; ω represents frequency; and c represents the speed of sound.
6. The method for high-quality imaging of full matrix data based on frequency domain reverse time migration as claimed in claim 1, wherein the specific steps of step (4) are as follows:
4.1, taking any frequency which is not traversed in the selected frequency band range as a current frequency, taking an excitation sound source matrix, a receiving sound source matrix and an impedance matrix corresponding to the current frequency as input, carrying out LU decomposition on the impedance matrix, and solving to obtain an excitation sound pressure field and a receiving sound pressure field corresponding to the current frequency;
4.2 go through all frequencies in the selected frequency band range, resulting in the transmit wavefield and the receive wavefield.
7. The method for high-quality imaging of full matrix data based on frequency domain reverse time migration as claimed in claim 6, wherein in step 4.1, the calculation formula for solving the sound pressure field generated by the sound source matrix according to the sound source matrix is as follows:
L(ω)U(ω)[P 1 (ω),P 2 (ω),…,P n (ω)]＝[S 1 (ω),S 2 (ω),…,S n (ω)]
wherein n representsThe number of excitations; l (omega) and U (omega) represent two factors obtained after LU decomposition of the impedance matrix A (omega); [ S ] 1 (ω),S 2 (ω),…,S n (ω)]Representing a matrix of excitation sources or a matrix of reception sources; [ P ] 1 (ω),P 2 (ω),…,P n (ω)]Representing either an excitation sound pressure field or a received sound pressure field.
8. The method for high-quality imaging of full matrix data based on frequency domain reverse time migration according to claim 1, wherein in the step (5), the expression of the imaging condition is:
wherein I (x, z) represents an initial image; n represents the number of excitations; n is a radical of ω Representing the total number of discrete frequencies in the selected frequency band; s. the i (x, z, ω) represents the transmit wavefield; r i (x, z, ω) represents the receive wavefield; denotes taking the complex conjugate.
9. The method of claim 1, wherein the initial image is post-processed as follows:
firstly, carrying out sharpening processing on an initial image by adopting Laplace filtering, then carrying out Hilbert transform on the sharpened initial image to obtain the envelope of the image, and finally filtering salt and pepper noise by adopting median filtering to obtain the imaging result.
10. The full matrix data high-quality imaging method based on frequency domain reverse time migration according to claim 1, characterized in that a 9-point finite difference method is adopted to discretize the sound velocity distribution.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210504512.6A CN114858926A (en) | 2022-05-10 | 2022-05-10 | Full matrix data high-quality imaging method based on frequency domain reverse time migration |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210504512.6A CN114858926A (en) | 2022-05-10 | 2022-05-10 | Full matrix data high-quality imaging method based on frequency domain reverse time migration |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN114858926A true CN114858926A (en) | 2022-08-05 |
Family
ID=82637359
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210504512.6A Pending CN114858926A (en) | 2022-05-10 | 2022-05-10 | Full matrix data high-quality imaging method based on frequency domain reverse time migration |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114858926A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115791974A (en) * | 2023-02-02 | 2023-03-14 | 汕头市超声检测科技有限公司 | Fitting curve-based all-focusing imaging directivity weighting calculation method |
-
2022
- 2022-05-10 CN CN202210504512.6A patent/CN114858926A/en active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115791974A (en) * | 2023-02-02 | 2023-03-14 | 汕头市超声检测科技有限公司 | Fitting curve-based all-focusing imaging directivity weighting calculation method |
| CN115791974B (en) * | 2023-02-02 | 2023-08-22 | 汕头市超声检测科技有限公司 | Full-focus imaging directivity weighting calculation method based on fitting curve |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106770664B (en) | A method of edge defect detection is improved based on total focus imaging algorithm | |
| CN113888471B (en) | High-efficiency high-resolution defect nondestructive testing method based on convolutional neural network | |
| Olofsson | Phase shift migration for imaging layered objects and objects immersed in water | |
| Skjelvareid et al. | Internal pipeline inspection using virtual source synthetic aperture ultrasound imaging | |
| US20080092661A1 (en) | Methods and system for ultrasound inspection | |
| Merabet et al. | 2-D and 3-D reconstruction algorithms in the Fourier domain for plane-wave imaging in nondestructive testing | |
| Kažys et al. | Ultrasonic detection and characterization of delaminations in thin composite plates using signal processing techniques | |
| CN107356670A (en) | A kind of ultrasonic phase array weld defect detection method based on oblique incidence | |
| Ji et al. | Ultrasonic full-matrix imaging of curved-surface components | |
| CN114858926A (en) | Full matrix data high-quality imaging method based on frequency domain reverse time migration | |
| CN108508093A (en) | A kind of detection method and system of workpiece, defect height | |
| Wu et al. | Synthetic aperture imaging for multilayer cylindrical object using an exterior rotating transducer | |
| Jin et al. | Frequency domain synthetic aperture focusing technique for variable-diameter cylindrical components | |
| Robert et al. | Assessment of real-time techniques for ultrasonic non-destructive testing | |
| Karaojiuzt et al. | Defect detection in concrete using split spectrum processing | |
| JP5456367B2 (en) | Phased array aperture synthesis processing method | |
| Zhang et al. | Multi-defect detection based on ultrasonic Lamb wave sign phase coherence factor imaging method | |
| Lukomski | Non-stationary phase shift migration for flaw detection in objects with lateral velocity variations | |
| Duan et al. | Multiframe ultrasonic TOFD weld inspection imaging based on wavelet transform and image registration | |
| Osegueda et al. | Detection of cracks at rivet holes in thin plates using Lamb-wave scanning | |
| Chang et al. | Frequency-wavenumber migration of ultrasonic data | |
| Huthwaite et al. | Guided wave tomography performance analysis | |
| Bazulin | Obtaining flaw images that take the effect of multiple ultrasonic pulse reflections from the boundaries of a test object into account | |
| Olofsson et al. | Ultrasonic imaging of immersed objects using migration techniques | |
| CN114487115B (en) | High-resolution defect nondestructive testing method based on combination of Canny operator and ultrasonic plane wave imaging |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |