CN102608205A

CN102608205A - Multilayer-object nondestructive testing ultrasonic imaging method based on variable wave speed phase shifting

Info

Publication number: CN102608205A
Application number: CN2012100461511A
Authority: CN
Inventors: 秦开怀; 杨春; 孙丰
Original assignee: Tsinghua University
Current assignee: Tsinghua University
Priority date: 2012-02-24
Filing date: 2012-02-24
Publication date: 2012-07-25
Anticipated expiration: 2032-02-24
Also published as: CN102608205B

Abstract

The invention discloses a multilayer-object nondestructive testing ultrasonic imaging method based on variable wave speed phase shifting, which belongs to the technical field of fast imaging of an inner structure of a multilayer heterogeneous object. The method is characterized in that under a condition that equal-interval moving detection can be carried out in the length direction of the multilayer heterogeneous object based on an ultrasonic transducer, in a homogeneous area, images in each row in the homogeneous area are calculated by a constant wave speed phase shifting ultrasonic imaging technology and a fast Fourier transform technology, in a heterogeneous area, a bounding box bounding medium boundaries is built, and for the heterogeneous areas positioned above and below the boundaries in the bounding box, images in each row in the heterogeneous area are calculated by a variable wave speed phase shifting method and an N-time single output fast Fourier transform technology. The multilayer-object nondestructive testing ultrasonic imaging method is suitable for the fast ultrasonic imaging of the multilayer heterogeneous object with heterogeneous media in the depth direction and the horizontal direction, and a respective treating method adopting variable wave speed as a characteristic can be adopted close to bonding surfaces of the media, so that the accuracy is high.

Description

Multilayer object nondestructive testing ultrasonic imaging method based on variable wave speed phase shift

Technical Field

The invention relates to an ultrasonic nondestructive testing technology, a synthetic aperture focusing imaging technology and a phase migration technology, and realizes the rapid imaging of the internal structure of a complex irregular layered object.

Background

Currently, in the field of ultrasonic nondestructive testing, for ultrasonic imaging of a layered object having an interface (including contact ultrasonic imaging and liquid immersion ultrasonic imaging, if a couplant layer is regarded as a part of a tested object, liquid immersion ultrasonic imaging can also be regarded as contact ultrasonic imaging of the layered object, so that both are collectively called ultrasonic imaging of the layered object in the present invention), there are two main approaches at present: the method adopts a method of combining a time domain synthetic aperture focused ultrasound imaging technology and a Ray Tracing (Ray Tracing) technology, and adopts a Phase Shift transfer (Phase Shift transfer) technology and a frequency domain synthetic aperture focused ultrasound imaging technology.

The fundamental idea of synthetic aperture focused ultrasound imaging (SAFT) is to use a single transducer to simulate an aperture array, to scan the object to be measured in an orderly manner, and to use a delay-and-add (DAS) method (time delay or phase delay) to perform focused imaging on the scanned pulse echo signals. The SAFT model of operation is shown in FIG. 1, where the Z direction is the depth direction and the X direction is the transducer scan direction, when the transducer is located at u_nAt a distance r from the target reflector (x ', z') to the transducer_n. The SAFT technique has a division between the time and frequency domains:

based on DAS principle and multi-point dynamic focusing technology, the time domain SAFT calculates different time delays for focusing points at different depths of a target imaging region. As in fig. 1, to focus at the target reflector (x ', z'), the time-domain SAFT technique subjects the echo signals obtained by the transducer at each scanning position within the synthetic aperture length L to a time-delayed superposition process: let s (u)_nT) transducer at u_nAt the received echo signal, t is the sampling time, u_nWith a delay of (x ', z') with respect to the target reflection object

Where c is the propagation velocity of ultrasound in the medium, then the imaging at (x ', z') is

Wherein, ω is_nFor apodization functions, the delta function is a unit sample sequence.

The frequency domain SAFT is called wave number algorithm (wavenumber algorithm), and is mainly a method for generating beams based on phase delay, which is derived from the inverse theory of wave equation. The basic principle is to perform a fourier decomposition of the Green function. The algorithm process comprises the steps of firstly carrying out two-dimensional Fourier transform on ultrasonic data to obtain a two-dimensional frequency spectrum, then realizing coordinate transformation by utilizing nonlinear mapping (stopmapping) and interpolation, and finally carrying out filtering processing on the transformed frequency spectrum and carrying out two-dimensional Fourier inverse transformation to generate a reconstructed image in a space-time domain.

The Phase Shift (Phase Shift Migration) ultrasonic imaging method is an ultrasonic imaging method obtained by combining a Migration Technique (Migration Technique) in Reflection Seismology (Reflection Seismology) with a frequency domain SAFT. The phase shift technology was originally proposed by Gazdag in 1978 and has been used in seismic/radar imaging, while the recently proposed phase shift ultrasonic imaging method has applied the technology to the field of ultrasonic nondestructive testing.

The phase transition ultrasonic imaging method regards an ultrasonic detection model as an explosive reflection model, the situation that the reflectors in the object to be detected explode simultaneously at the moment when t is equal to 0 is assumed, the explosion intensity of each reflector is proportional to the reflectivity of the reflector, and the whole field intensity is measured by a group of transducers. The main idea is to extrapolate from the sound field observed at the horizontal position (i.e. the first line in the depth direction) to calculate the sound field at other positions in the depth direction. The specific algorithm comprises two main steps: the first step is to perform two-dimensional Fourier transform on the time domain data to obtain a two-dimensional frequency spectrum,

S(k，ω)＝2D-FT(s(u_n，t))

the second step is a cycle in the depth direction: firstly, phase shift is carried out on the two-dimensional frequency spectrum obtained by the last circulation, then two-dimensional Fourier inversion is carried out, t is taken as 0, a line of image is obtained,

s(u_n，t＝0)＝2D-IFT(S(k，ω)·α(Δz，k，ω))

wherein,

the phase migration amount corresponding to the step length of Δ z in the depth direction, and c is the propagation speed of the ultrasound in the medium, and the value of c is constant, so that the ultrasound imaging method uses a phase migration technique with a constant wave speed.

From the DAS principle, the time-domain SAFT algorithm needs to determine the delay according to the relative distance between the transducer and the reflector and the wave velocity. For a multilayer object, acoustic properties such as refraction and reflection of sound waves occur at the interfaces between layers, so that the propagation path changes. In addition, the media of each layer of the object are different, and the propagation speeds of the sound waves in each layer are usually also different, so that the calculation of the time delay cannot adopt a simple calculation method of dividing the relative distance by the wave speed, but the propagation path of the wave beam is obtained first, and then the propagation time is obtained section by section. The important and difficult point is to find the propagation path of the beam quickly and accurately, and the ray tracing technology is the best method for achieving the purpose. The principle of ray tracing is mainly based on the Snell theorem or the Fermat principle, and the shortest time-consuming sound wave propagation path is found out through iterative calculation. The advantage of the ray tracing technology is that the technology has no special requirements on the shape of the boundary surface of the multilayer object and can be suitable for ultrasonic imaging of objects with any shapes. However, the serious disadvantage is that the ray tracking algorithm includes iterative operation, the time complexity is high, and meanwhile, the DAS principle is actually a convolution operation, the calculated amount is large, and the imaging time is long. As shown in fig. 2(b), the protocol takes 2 hours to generate the image (5cm by 5 cm). This disadvantage severely limits the spread of the technology.

For imaging of a multilayer object, the phase shift ultrasonic imaging method does not require ray paths and convolution operations, and the imaging speed can be greatly improved, for example, only 4 minutes is needed to generate an image shown in fig. 2 (c). And the phase shift technology is an image reconstruction method based on a frequency domain, has higher distance resolution, higher azimuth resolution, lower side lobe and higher imaging quality than the SAFT technology in a time domain. However, since the wave velocity of the acoustic wave is assumed to be constant in the phase shift technology, this method can only be applied to the regular case where there are heterogeneous media in the depth direction and the same media in the horizontal direction, i.e. the interfaces between the layer media are horizontal or parallel to each other, and for the object where there are heterogeneous media in both the depth direction and the horizontal direction, i.e. the object including the irregular interfaces (curves/curved surfaces), the detection result is not accurate enough due to the different wave velocities of the acoustic wave in the media, as shown in fig. 2(c), the imaging result of the second layer interface has been seriously wrong.

Thus, there is currently a lack of rapid, accurate and efficient methods for ultrasound imaging of complex layered objects containing irregular interfaces.

Disclosure of Invention

The invention aims to provide an ultrasonic imaging method, which can realize quick, accurate and effective imaging of a complex layered object containing an irregular interface.

The invention is characterized in that the invention can image the regular layered object or the object only containing single medium, and can image the irregular layered object containing non-horizontal and non-parallel interfaces, and the imaging speed is high.

The invention is characterized in that a phase migration ultrasonic imaging technology is adopted, the phase migration technology is improved, the condition that different media exist is considered in the depth direction and the horizontal direction of a measured object respectively, different wave speeds are introduced into a phase migration formula, and the phase migration technology with variable wave speed is provided, so that the ultrasonic imaging method is suitable for irregular layered objects. The whole structure is as shown in fig. 3, the boundary lines are processed separately in the imaging process, the phase shift formula and the fast fourier transform technique of constant wave velocity are used outside the depth range of the boundary lines, and the phase shift technique and the N-time single output fast fourier transform technique of variable wave velocity are used in the range, so as to reduce the calculation amount.

The invention is characterized by comprising the following steps in sequence:

step (1): constructing a system which is composed of a computer, an ultrasonic transducer, a set of positioning controllers and an analog-to-digital converter and is used for carrying out nondestructive ultrasonic imaging on a vertical section formed in two directions of depth and horizontal of a three-layer heterogeneous object with heterogeneous media in the transverse direction and the longitudinal direction based on variable wave speed phase shift, wherein:

the ultrasonic transducer is provided with: the pulse signal input end is connected with the output end of the positioning controller, the input end of the positioning controller is connected with the corresponding positioning control signal output end of the computer, and the ultrasonic transducer is also provided with: and the output end of the analog-to-digital converter is connected with the echo sampling signal input end of the computer. The transducer is controlled by the positioning controller to move at a fixed speed of 1 step/ms on the surface of a three-layer heterogeneous object, the positioning controller is a transmission device for controlling the moving position of the ultrasonic transducer, and parameters of the transmission device are input by the computer,

the horizontal length of the three-layer heterogeneous object along the x-axis direction is L, the three-layer heterogeneous object is divided into L/delta x intervals, delta x is the interval length, and the interval length is from a point 0 to an end point N of the ultrasonic transducer on the x-axis _x1 step size per move, N_xThe point reached by the ultrasonic transducer in each movement is called probe point, and the total number of the probe points is N_xNumber n_x＝0，1，...，N_x-1, said positioning controller generating a TTL pulse at each detection point triggering said ultrasound transducer to transmit an excitation pulse in a depth direction z perpendicular to the x-axis of the three-layered heterogeneous object, then the transducer switches to a receiving mode and starts timing to receive the echo signal reflected from the object to be measured, said analog-to-digital converter couples said transducer at a detection point n_xProcessing the received echo signal by N_tSub-sampling and storing in computer, sampling number n_t＝0，1，...，N_t-1, sampling frequency f_sThe value is predetermined for the analog-to-digital converter, and s (z, x, t) is represented by x ═ n on the abscissa_xThe echo signal received at Δ x is n at t_t/f_sSampling values of the moments;

step (2): the computer starts from n_tStarting to read the probe point n sequentially when the probe point n is 0_xThe sampling value at 0 and the sampling sequence number of the first fluctuation of the sampling value is recorded

Subscript 0 indicates the number of probe points, and then n is sequentially read by repeating the above steps_x＝1，...，N_x-1 sampling values at the probe points and recording the number of first fluctuations of the sampling values at the points

To obtain N_xN in total on each probe point_xNumber of samples fluctuating for the first time

From

Select the smallest serial number

And calculating the depth coordinate value corresponding to the sequence number in the z direction of the three-layer heterogeneous object

And from

Select the largest serial number

Wherein v is₁The wave velocity in the first layer medium is determined by two horizontal lines parallel to the x-axis by the two depth coordinate values, respectively, and the boundary between the first layer medium and the second layer medium is enclosed to form a three-layer heterogeneous objectA first bounding box spanning both the first and second layers of media in the x-z longitudinal section of the volume;

and (3): carrying out one-dimensional fast Fourier transform (1D-FFT) on sampling values s (z, x, t) at all detection points one by one along a time axis t_t(s (z, x, t)), and taking the initial value z of z₀Obtaining an x-omega domain two-dimensional signal:

S(z₀，x，ω)＝1D-FFT_t(s(z₀，x，t))

wherein x is the serial number n of the detection point_xAt horizontal coordinate value, where ω is sampling angular frequency, and then, for the above x- ω domain two-dimensional signal S (z)₀X, ω) are sequentially subjected to one-dimensional fast fourier transform 1D-FFT along a horizontal axis x_x(S(z₀X, ω)) to yield z₀K-omega domain two-dimensional spectrum S (z)₀，k，ω)：

S(z₀，k，ω)＝1D-FFT_x(1D-FFT_t(s(z₀，x，t)))

Wherein k is the x-component of the wavenumber vector with the sequence number n_kThe wavenumber vector represents a spatial gradient of the wave phase,

0≤n_k≤N_x-1, and n_kIs an integer;

0≤n_ω≤N_t-1, and n_ωIs an integer;

wherein n is_ωThe value of the sampling angular frequency omega is 0 to (N)_t-1), take t ═ t₀One-dimensional sampled data s (z) of 0₀，x，t₀) I.e. the z-th ═ z on the vertical section of the three-layer heterogeneous object₀Each in 0 rowPixel values of the pixels;

and (4): generating a z-direction homogeneous region z by the following steps₀＜z≤z_minThe inner images of all rows formed by pixel points of all rows respectively:

step (4.1): calculating the profile in z using the following phase shift equation_minTwo-dimensional spectrum S (z) of_min，k，ω)：

Wherein,

wherein i is an imaginary unit, z ₀0, ω is the sampling angular frequency, v_zThe wave velocity of the sound wave in the z-th row of medium of the measured object in the depth direction is equal to or less than z_minWhen, v_z＝v₁，v₁Is the wave velocity of the acoustic wave in the first layer of medium, and Δ z is the separation distance between two adjacent pixels in the depth direction z, and α (Δ z, k, ω, v)_z) Is the amount of phase migration at a constant wave velocity,

step (4.2): the following steps are performed in sequence:

step (4.2.1): get t ═ t₀(iv) 0, for S (z) obtained in the step (4.1)_minK, ω) is taken as the gaussian integral of the independent variable ω,

step (4.2.2): conjugating the Gaussian integral result of the step (4.2.1), performing one-dimensional fast Fourier transform on k, conjugating the calculation result and multiplying by 1/N_xObtaining the coordinate value z in the depth direction in the first layer medium_minThe image formed by the row of pixel points:

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>min</mi> </msub> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msub> <mi>N</mi> <mi>x</mi> </msub> </mfrac> <mo>{</mo> <mn>1</mn> <mi>D</mi> <mo>-</mo> <msub> <mi>FFT</mi> <mi>k</mi> </msub> <mo>{</mo> <mo>[</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>π</mi> </mrow> </mfrac> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>π</mi> </mrow> <mi>π</mi> </msubsup> <mi>S</mi> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>min</mi> </msub> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <mi>dω</mi> <msup> <mo>]</mo> <mo>*</mo> </msup> <mo>)</mo> <msup> <mo>}</mo> <mo>*</mo> </msup> </mrow> </math>

in the z-direction the first layer homogeneous zone z₀＜z≤z_minAnd in the inner step, the images of all rows are the same:

s(z，x，t₀)＝s(z_min，x，t₀)，z₀＜z≤z_min；

and (5): taking v in step (4)_z＝v₁And circularly performing the following steps (5.1) and (5.2) by taking the step delta z as the step size for the z value until z is more than or equal to z_maxGenerating a depth direction z in the first bounding box_min＜z≤z_maxEach line of image formed by each line of pixel points in the interval respectively comprises the following steps:

step (5.1): calculating a two-dimensional frequency spectrum at z + Δ z using the phase shift formula at the constant wave velocity, z_min＜z+Δz≤z_max：

S(z+Δz，k，ω)＝S(z，k，ω)·α(Δz，k，ω，v_z)，

Step (5.2): the following steps are performed in sequence:

step (5.2.1): get t ═ t₀(ii) the gaussian integral of the argument ω is taken for S (z + Δ z, k, ω) obtained in the step (5.1) described above as 0,

step (5.2.2): conjugating the Gaussian integral result of the step (5.2.1), performing one-dimensional fast Fourier transform on k, conjugating the calculation result and multiplying by 1/N_xObtaining an image formed by a row of pixel points with the coordinate value of z + delta z in the depth direction, z_min＜z+Δz≤z_max：

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>+</mo> <mi>Δz</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msub> <mi>N</mi> <mi>x</mi> </msub> </mfrac> <mo>{</mo> <mn>1</mn> <mi>D</mi> <mo>-</mo> <msub> <mi>FFT</mi> <mi>k</mi> </msub> <mo>{</mo> <mo>[</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>π</mi> </mrow> </mfrac> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>π</mi> </mrow> <mi>π</mi> </msubsup> <mi>S</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>+</mo> <mi>Δz</mi> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <mi>dω</mi> <msup> <mo>]</mo> <mo>*</mo> </msup> <mo>)</mo> <msup> <mo>}</mo> <mo>*</mo> </msup> <mo>;</mo> </mrow> </math>

And (6): correcting the depth direction z in the first bounding box as follows_minTo z_maxInterval z_min＜z≤z_maxTo eliminate errors caused by inhomogeneities between the second layer and the first layer of media:

step (6.1): z obtained in step (5)_minTo z_maxTaking an interval image block as an input quantity, and extracting z by using a Canny operator edge extraction algorithm_min、z_maxThe boundary c between the first layer medium and the second layer medium₁(x，z)，

Step (6.2): from z to z_minInitially until z ≧ z_maxZ obtained in said step (4.1)_minTwo-dimensional spectrum S (z) of_minK, ω) as a starting value, and step (6.2.1) and step (6.2.2) are cyclically executed with Δ z as a step size:

step (6.2.1): and calculating to obtain an x-omega domain two-dimensional signal at a z + delta z position in the first bounding box by using an N-time single-output inverse fast Fourier transform method according to the following formula:

wherein the value of N is equal to N_x，n_x＝0，1，...，N_x-1，n_k＝0，1，...，N_x-1，α(Δz，k，ω，v_z(x) Phase shift amount as a wave velocity:

wherein, the values of k and omega are the same as the step (3), v_z(x) Is the wave velocity of the acoustic wave in the medium at coordinates (x, z):

v₁to be located at the boundary line c₁Speed of propagation, v, of acoustic waves in a first layer of medium above (x, z)₂To be located at the boundary line c₁(x, z) the propagation velocity of the acoustic wave in the second layer medium below (x, z),

the calculation steps are as follows:

step (a): respectively for different n_x，n_k(n_x＝0，1，...，N_x-1；n_k＝0，1，...，N_x-1) calculating

Step (b): for different n_x(n_x＝0，1，...，N_x-1) with G (n)_x，n_k)，n_k＝0，1，...，N_x-1 sequence as input to the clipped output FFT algorithm, for n_kAs G (n)_x，n_k) Inverse discrete Fourier transform, but taking only n_xThe value of this time-domain point is taken as an x- ω domain two-dimensional signal at a coordinate value z + Δ z in the depth direction:

step (6.2.2): get t ═ t₀(iv) when S (z + Δ z, x, ω) obtained in the step (6.2.1) is subjected to gaussian integration of the argument ω to obtain an image formed by one line of pixel points whose coordinate value in the depth direction is z + Δ z, z being equal to 0_min＜z+Δz＜z_max：

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>+</mo> <mi>Δz</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>π</mi> </mrow> </mfrac> <msubsup> <mo>&Integral;</mo> <mrow> <mo>-</mo> <mi>π</mi> </mrow> <mi>π</mi> </msubsup> <mi>S</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>+</mo> <mi>Δz</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <mi>dω</mi> </mrow> </math>

Meanwhile, storing S (z + delta z, x, omega) and performing one-dimensional Fourier transform to obtain S (z + delta z, k, omega) which is used as an initial value of the next iteration;

and (7): according to the method described in step (5), v is taken_z＝v₂And circularly executing the step (5.1) and the step (5.2) by taking the step deltaz as a step size for the z value until z is more than or equal to z_depthTo generate a depth direction z_max＜z≤z_depthEach row image v formed by each row of pixel points in interval₂Is the propagation velocity of the acoustic wave in the second layer medium, z_depthThe maximum coordinate value of the three-layer heterogeneous object in the depth direction z is represented as a system preset value;

and (8): z obtained in step (7)_maxTo z_depthTaking an interval image block as an input quantity, and extracting z by using a Canny operator edge extraction algorithm_max、z_depthThe boundary c between the second layer medium and the third layer medium₂(x, z) determining the boundary line c₂(x, z) minimum coordinate value z 'in depth direction z'_minAnd maximum coordinate value z'_maxRespectively through said z'_minAnd z'_maxTwo depth coordinate values are used for making two horizontal lines parallel to the x axis, and a boundary between the second layer of medium and the third layer of medium is surrounded to form a second surrounding box which is positioned on the x-z vertical section of the three-layer heterogeneous object and simultaneously spans the second layer of medium and the third layer of medium;

and (9): correcting the depth direction z 'in the second surrounding box according to the method of the step (6.2)'_minTo z'_maxZone z'_min＜z≤z′_maxThe images of the rows in the image table are used for eliminating errors caused by the heterogeneity between the third layer medium and the second layer medium;

step (10): according to the method described in step (5), v is taken_z＝v₃Z is z'_maxThe step (5.1) and the step (5.2) are circularly executed by taking the step delta z as a step size until the step z is more than or equal to z_depthGenerating images, v, of all rows formed by pixel points of all rows except the second bounding box in the third layer of medium₃Is the propagation speed of the acoustic wave in the third layer medium.

The invention has the advantages that the imaging precision is high, and for the measured object shown in figure 2(a), if the diameter of the transducer is 0.5mm, the moving step length of the transducer is 0.5mm, the central frequency of the ultrasonic wave emitted by the transducer is 5MHz, the sampling frequency is 100MHz, and the step length Delta z in the depth direction is 0.05mm, the error of the interface of the first layer can be controlled within 0.4mm, and the error of the interface of the second layer can be controlled within 1.0mm (as shown in figure 4).

Compared with the prior art, the invention has the following advantages:

1. the fast ultrasonic imaging can be carried out on the object with the irregular interface;

2. the imaging resolution and precision are high.

Drawings

Fig. 1 is a working model of SAFT: x-direction is the transducer scanning direction, Z-direction is the depth direction, u_nFor the position of the transducer when it is moved in equal steps, i.e. the position of the probe point, (x ', z') are the coordinates of the target reflector, r_nThe distance of the target reflector to the transducer.

FIG. 2 is a comparison of ultrasound imaging of a non-regular layered object under the same experimental environment: 2(a) is a sectional view of the object to be measured; 2(b) is an image generated by combining a time domain SAFT and a ray tracing technology with the 2(a), the imaging time is more than 2 hours, and the image is very sensitive to the parameters of an algorithm and not very good in robustness; 2(c) is an image generated by a phase shift ultrasonic imaging technology with constant wave speed, the imaging time is 4 minutes, and the error is large; 2(d) is an image generated by the imaging method, and the imaging time is 6 minutes.

Fig. 3 is a schematic view of the overall process of the present imaging method.

FIG. 4 is an analysis chart of experimental data of FIG. 2 (d): 4(a) is a first layer boundary error graph, a real curve is a standard value of a first layer boundary of a measured object, a virtual point curve is the first layer boundary obtained by the imaging method, and a dashed curve is an error curve of the first layer boundary obtained by the imaging method; and 4(b) is an error graph of the boundary of the second layer, the real curve is a standard value of the boundary of the second layer of the object to be measured, the imaginary point curve is the boundary of the second layer obtained by the imaging method, and the dashed curve is the error curve of the boundary of the second layer obtained by the imaging method.

Fig. 5 is a flow chart of the ultrasonic imaging system.

Figure 6 is a block diagram of the present ultrasound imaging hardware system.

Fig. 7 is an ultrasound transducer operational schematic.

FIG. 8 is a flow chart of a multi-layer object nondestructive testing ultrasonic imaging algorithm based on a variable wave speed phase shift technique.

Fig. 9 is a schematic diagram of a bounding box in the present ultrasound imaging method.

Fig. 10 is a schematic diagram of FFT butterfly operation when N is 8.

Detailed Description

The implementation process of the invention comprises three parts (as shown in figure 5): ultrasound data acquisition, imaging calculations, and image display. The structure diagram of the hardware platform system is shown in fig. 6, the ultrasonic imaging system is composed of a computer, an ultrasonic transducer, a set of positioning controller and an analog-to-digital converter, the pulse signal input end of the ultrasonic transducer is connected with the output end of the positioning controller, and the input end of the positioning controller is connected with the positioning control signal output end of the computer. The echo signal output end of the ultrasonic transducer is connected with the input end of the analog-to-digital converter, and the output end of the analog-to-digital converter is connected with the echo sampling signal input end of the computer.

The mode of ultrasonic data acquisition can be selected as follows: contact type, that is, the ultrasonic transducer directly contacts the surface of the object to be measured; or liquid immersion type, that is, the object to be measured is placed in liquid, and the ultrasonic transducer measures on the surface of the liquid. The transducer may use a contact probe or a liquid immersion probe depending on the detection mode. The system uses a single transmit/receive transducer that is moved by a positioning controller in the x direction (as in fig. 7) at a fixed rate of about 1 step/ms at a uniform step size deltax over the surface of the object being measured or in the liquid. The controller generates a TTL (transistor-transistor logic level) pulse at the moment that each target position is stable, the TTL pulse is used for triggering the ultrasonic transducer to transmit an excitation pulse to the depth direction of the measured object, the depth direction of the measured object is perpendicular to the x direction, then the ultrasonic transducer is switched to a receiving mode, timing is started, and echo waves reflected from the measured object are received. The position of each transducer transmitting pulse and receiving echo is a probe point. Echo signals received by the transducer are acquired by an analog-to-digital converter and stored in a memory. The Δ x needs to be determined according to the actual size of the object to be measured and the imaging precision requirement, and the smaller the value of Δ x, the more accurate the generated image, but the longer the calculation time.

The imaging calculation is that sampling data of a measured object at each detection point on a longitudinal section is used as computer input, and then a longitudinal section image of the measured object is calculated according to the following imaging steps (see a flow chart in figure 8), wherein s (z, x, t) represents a sampling value of a signal received by a transducer in the horizontal direction x at the time t; z represents a coordinate value in the depth direction; Δ z represents a step length in the depth direction (related to precision), that is, the distance between two adjacent pixel points in the longitudinal direction (i.e., the depth direction) of the generated image; z is a radical of_depthMaximum depth value in the longitudinal direction for a preset image:

step (1): the first layer of medium is homogeneous medium, and the ultrasonic reflection signal is the same, and the ultrasonic reflection signal generates fluctuation at the boundary of the first layer and the second layer. Searching the first signal fluctuation moment at each detection point according to the sampling data of the ultrasonic echo signal, and finding out the minimum moment value t from the fluctuation moments_minAnd a maximum time value t_maxThen according to the wave speed v in the first layer medium₁Finding a bounding box of a boundary between the first layer medium and the second layer medium to obtain two boundary lines of the bounding box in the depth direction (namely the z direction): z is a radical of_min＝v₁·t_minAnd z_max＝v₁·t_max(see FIG. 9);

step (2): taking the initial value z in the depth direction as z₀Performing two-dimensional Fourier transform on the parameters x and t to transform the original ultrasonic echo signal into a k-omega domain signal:

S(z₀，k，ω)＝∫∫s(z₀，x，t)e^-ikxe^-iωtdxdt (1)

wherein, ω is the sampling angular frequency, and k is the x-direction component (wave number vector refers to the spatial gradient of the wave phase) of the wave number vector (wave number vector);

and (3): generating a depth direction z₀To z_minInterval (z)₀≤z≤z_min) The image of (2). Because the media of the measured object in the section of interval are completely the same, the section of image is also completely the same:

step (3.1): calculating z using a phase shift formula_minK- ω domain signal of (c):

wherein,

i is an imaginary unit, v_zThe wave velocity of the sound wave in the z-th row of medium of the measured object in the depth direction is obtained, and when z is not less than z_minWhen, v_z＝v₁，v₁Is the wave velocity of the sound wave in the first layer of medium;

step (3.2): making two-dimensional Fourier inverse transformation, making t equal to t₀When the value is 0, a line of image with the longitudinal coordinate value z is obtained:

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>s</mi> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>min</mi> </msub> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mrow> <mn>4</mn> <mi>π</mi> </mrow> <mn>2</mn> </msup> </mfrac> <mo>&Integral;</mo> <mo>&Integral;</mo> <mi>S</mi> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>min</mi> </msub> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mi>ikx</mi> </msup> <msup> <mi>e</mi> <mi>iωt</mi> </msup> <mi>dkdω</mi> <mo>,</mo> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mn>0</mn> </msub> <mo>≤</mo> <mi>z</mi> <mo>≤</mo> <msub> <mi>z</mi> <mi>min</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow> </math>

and (4): generating a depth direction z_minTo z_maxInterval (z)_min＜z≤z_max) The image of (2). Get v_z＝v₁Step (4.1) and step (4.2) are performed cyclically in Δ z steps for z values until z ≧ z_max：

Step (4.1): the k- ω domain signal at z + Δ z is calculated using the phase shift formula:

S(z+Δz，k，ω)＝S(z，k，ω)α(Δz，k，ω，v_z) (4)

step (4.2): let z ← z + Δ z, t ═ t₀And (5) performing two-dimensional inverse Fourier transform to obtain a line of image with a longitudinal coordinate value z:

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mrow> <mn>4</mn> <mi>π</mi> </mrow> <mn>2</mn> </msup> </mfrac> <mo>&Integral;</mo> <mo>&Integral;</mo> <mi>S</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mi>ikx</mi> </msup> <msup> <mi>e</mi> <mi>iωt</mi> </msup> <mi>dkdω</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow> </math>

and (5): correction of z_minTo z_maxInterval (z)_min＜z≤z_max) The image of (2). Because the calculation in the step (4) only usesThe wave velocity of the first layer, i.e. the interval is still regarded as the same medium, and the obtained calculation result necessarily has an error, so that the step needs to be corrected. The method specifically comprises two steps:

step (5.1): for z obtained in step (4)_minTo z_maxExtracting z from the image of the interval by using a Canny edge extraction operator_minAnd z_maxThe boundary line c (x, z), z between_min≤z≤z_max. The calculation process of the Canny edge extraction operator comprises four steps:

step (a): and (6) smoothing the image. Smoothing and filtering the original image by using a two-dimensional Gaussian function to obtain an image I (x, z);

step (b): the gradient M and direction Q of each pixel of the image I (x, z) are calculated. Using a 2 x 2 template as a first order approximation of partial differentiation in the x-and z-directions, i.e.

p = \frac{1}{2} (\begin{matrix} - 1 & 1 \\ - 1 & 1 \end{matrix}),

q = \frac{1}{2} (\begin{matrix} 1 & 1 \\ - 1 & - 1 \end{matrix})

The magnitude M and direction Q of the gradient are then

Q = \arctan (\frac{q}{p});

Step (c): non-maxima suppression of the gradient image. For each pixel point I (x, z), if the gradient value M (x, z) of each pixel point is smaller than the gradient amplitude of two adjacent points along the gradient direction Q (x, z), which indicates that the point is not an edge point, setting the gray level of I (x, z) as 0;

step (d): and (4) carrying out double thresholding processing. And setting a Low threshold Low and a High threshold High required by the double-threshold method for detecting and connecting edges, and performing double-thresholding processing on the gradient image. Edges where the gradient magnitude is greater than the High threshold High; if the gradient magnitude is less than the Low threshold Low, it is not an edge; if the gradient amplitude is between the two gradient amplitudes, judging whether an edge pixel with a higher threshold value High exists in eight adjacent pixels of the pixel, if so, judging that the edge pixel exists, otherwise, judging that the edge pixel does not exist;

step (5.2): for z value from z to z_minInitially, steps (5.2.1) and (5.2.2) are performed cyclically in Δ z steps until z ≧ z_max：

Step (5.2.1): calculating the k-omega domain signal at z + delta z by using a phase migration formula of variable wave speed:

S(z+Δz，k，ω)＝S(z，k，ω)α(Δz，k，ω，v_z(x)) (6)

wherein, α (Δ z, k, ω, v)_z(x) Phase shift amount as a wave velocity:

vz (x) takes the following values:

v₁，v₂the propagation velocities of the acoustic waves in the upper and lower layers of media where the boundary c (x, z) is located are respectively equal;

step (5.2.2): let z ← z + Δ z, t ═ t₀And (5) performing two-dimensional inverse Fourier transform to obtain a line of image with a longitudinal coordinate value z:

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msup> <mrow> <mn>4</mn> <mi>π</mi> </mrow> <mn>2</mn> </msup> </mfrac> <mo>&Integral;</mo> <mo>&Integral;</mo> <mi>S</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mi>ikx</mi> </msup> <msup> <mi>e</mi> <mi>iωt</mi> </msup> <mi>dkdω</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow> </math>

and (6): generating a depth direction z_maxTo z_depthInterval (z)_max＜z≤z_depth) The image of (2). Let v_z＝v₂Calculating z according to the method of step (4)_maxTo z_depthEach line image of the interval;

and (7): the subsequent demarcation lines are processed. For the image part (z) generated in step (6)_max＜z≤z_depth) Extract the new boundary c (x, z) using the Canny edge extraction operator, let v₁，v₂The propagation velocities of acoustic waves in the upper and lower layers of the medium, where the boundary c (x, z) is located, respectively, are set to z_minAnd z_maxA new bounding box is formed for the upper and lower boundaries of the new boundary line c (x, z) in the depth direction, respectively, if z is_min＝＝z_max＝＝z_depthIf yes, ending; otherwise, executing step (5.2) to correct z_minTo z_maxAn image of the section;

and (8): if z is_max＝＝z_depthIf yes, ending; otherwise, repeating the step (6) and the step (7).

In order to implement the two-dimensional fourier transform (as formula (1)) in step (2), the invention adopts a method of performing two-dimensional fast fourier transforms (1D-FFT) on the parameters x and t in sequence to improve the calculation speed, that is, the formula (1) is implemented as:

S(z₀，k，ω)＝1D-FFT_x(1D-FFT_t(s(z₀，x，t)))

in the step (3.2) and the step (4.2), two-dimensional inverse Fourier transform is required and t is taken as t ═ t₀If the value of 0 is used as the image data, t may be set to t first₀And integrate ω and then k to reduce the amount of data for the second dimension integration operation. Then equations (3) (5) are:

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>π</mi> </mrow> </mfrac> <mo>&Integral;</mo> <mo>{</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>π</mi> </mrow> </mfrac> <mo>&Integral;</mo> <mi>S</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <mi>dω</mi> <mo>}</mo> <msup> <mi>e</mi> <mi>ikx</mi> </msup> <mi>dk</mi> </mrow> </math>

wherein the integral of ω

Can be quickly obtained by using a numerical integration method (such as Gaussian integration). The inverse fourier transform of k is still implemented by one-dimensional fast fourier transform (1D-FFT), but the data obtained by numerical Integration is first conjugated ([ Integration ]_ω(S(z，k，ω))]Directly operating by using a 1D-FFT subprogram, conjugating the operation result once, and multiplying by 1/N_xAnd obtaining the product. N is a radical of_xFor transducer movement steps, i.e. N_x＝1+L_x/. DELTA.x, wherein L_xIs the length of the object to be measured in the x direction. Therefore, the specific implementation process of the equations (3) (5) is expressed as:

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <msub> <mi>N</mi> <mi>x</mi> </msub> </mfrac> <mo>{</mo> <mn>1</mn> <mi>D</mi> <mo>-</mo> <msub> <mi>FFT</mi> <mi>k</mi> </msub> <mrow> <mo>(</mo> <mo>[</mo> <mi>Integratio</mi> <msub> <mi>n</mi> <mi>ω</mi> </msub> <mrow> <mo>(</mo> <mi>S</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <msup> <mo>]</mo> <mo>*</mo> </msup> <mo>)</mo> </mrow> <msup> <mo>}</mo> <mo>*</mo> </msup> </mrow> </math>

the phase shift amount α (Δ z, k, ω, v)_z) Needs to be addressed for each velocity v_zThe calculation is performed for each of k and ω, and the results are stored in a three-dimensional array. Wherein, k and omega take the values as follows:

0≤n_k≤N_x-1, and n_kIs an integer;

0≤n_ω≤N_t-1, and n_ωIs an integer;

where Δ x is the step size of the transducer movement in the x-direction, f_sFor the sampling frequency (the value of which is predetermined by the system), N_x＝1+L_xΔ x is the number of probing points of the transducer in the x-direction, L_xΔ x is the number of transducer moving steps, L_xIs the length (x direction) of the object to be measured, N_tThe number of times the transducer samples the echo signal at each probe point. Due to phaseThe bit migration ultrasonic imaging technology regards the working model as an explosive reflection model, and only the part of the sound wave propagating upwards (along the negative direction of the z axis) from the reflector is considered, so the phase migration quantity alpha (delta z, k, omega, v)_z) In the calculation of (d), only the value of the part ω < 0 is taken, and the phase shift amount of the part ω ≧ 0 is made 0, and if b < 0, the value of (d) is taken

According to the periodicity of the value of omega, N can be taken in the calculation when omega is required to be less than 0_t/2≤n_ω＜N_tAnd (4) partial.

In step (5.2), the phase shift amount α (Δ z, k, ω, v) is varied in wave velocity_z(x) Is a function of x, the wave velocity v_zNot constant, so step (5.2.1) can no longer be calculated directly, but needs to be incorporated into the inverse fourier transform of step (5.2.2), i.e.

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>+</mo> <mi>Δz</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <msub> <mi>t</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>π</mi> </mrow> </mfrac> <mo>&Integral;</mo> <mo>{</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>π</mi> </mrow> </mfrac> <mo>&Integral;</mo> <mo>[</mo> <mi>S</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <mi>α</mi> <mrow> <mo>(</mo> <mi>Δz</mi> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>,</mo> <msub> <mi>v</mi> <mi>z</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mo>]</mo> <msup> <mi>e</mi> <mi>ikx</mi> </msup> <mi>dk</mi> <mo>}</mo> <msup> <mi>e</mi> <mi>iωt</mi> </msup> <mi>dω</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow> </math>

The outer layer integral in the formula (8) can be obtained by using a numerical integration method (such as Gaussian integral), and the inner layer integral operation (such as the formula (9)) needs to be obtained by modifying a discrete Fourier inverse transformation calculation formula (such as the formula (10)),

<math> <mrow> <mi>s</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>+</mo> <mi>Δz</mi> <mo>,</mo> <mi>x</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <mi>π</mi> </mrow> </mfrac> <mo>&Integral;</mo> <mo>[</mo> <mi>S</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>)</mo> </mrow> <mi>α</mi> <mrow> <mo>(</mo> <mi>Δz</mi> <mo>,</mo> <mi>k</mi> <mo>,</mo> <mi>ω</mi> <mo>,</mo> <msub> <mi>v</mi> <mi>z</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mo>]</mo> <msup> <mi>e</mi> <mi>ikx</mi> </msup> <mi>dk</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow> </math>

in the formula (10), n_xThe number of the detection point where the transducer is positioned in the x direction is equal to n_kAnd α (Δ z, k, ω, v: (Δ z, k, ω), v: (v:)_zx)) is an argument of n_kAnd n_xThe Inverse Discrete Fourier Transform (IDFT) or Inverse Fast Fourier Transform (IFFT) method cannot be used directly. Therefore, the present invention uses the N-times single output inverse fast fourier transform method (nsiffcoating). The method is based on the idea of butterfly operation in fast Fourier transform, but each transform only outputs a value in one transform domain. The specific algorithm of nsifftplanning is as follows:

For n_x＝0 to N_x-1

For n_k＝0 to N_x-1

end

S(z+Δz，x＝n_x·Δx，ω)＝IFFTPruning(G，N_x，n_x)

end

wherein when n is_xFixed time, function IFFTsounding (G, N)_x，n_x) Calculating G (n)_x，n_k)，n_k＝0，1，…，N_xInverse discrete Fourier transform of-1 sequence, but taking only n_xThe value of this time-domain point, i.e.

Equation (11) can be calculated using IDFT or IFFT algorithm, but unlike the general IDFT problem: equation (11) requires only the calculation of G (n)_x，n_k) At n_xIf the IDFT or IFFT algorithm is used directly, there must be a large number of redundant calculations, for example, in fig. 10, if only n needs to be calculated_xAt 2, the butterfly operation in the dotted line is redundant. To reduce the amount of computation, the present invention uses an FFT algorithm (ifftrounding) that clips the output. The algorithm is still based on the FFT algorithm, a butterfly operation method is adopted, and the principle is as follows: according to the butterfly grouping method in FFT algorithm, the calculation is divided into H steps (Stage), wherein H satisfies 2^H＝N_x. In step H (H ═ 1.., H), N is reacted with N_xPoint input sequence G (n)_x，n_k)，n_k＝0，1，…，N_x-1 average division into 2^H-hGroups, each group being incremented by 1 number b from 0_h(b_h＝0，1，...，2^h-1) in each set of butterflies except for the number b_h＝b_h+1mod 2^h-1The remaining points are redundant. In each calculation step, only the number b of each butterfly operation is calculated_hPoint (2) of (c). As shown in FIG. 10, N_x＝8，n_xWhen 2, the number b in each butterfly calculation is calculated in step (1)₁The point 0, numbered b in each set of butterflies, is calculated in step (2)₂At point 2, the number b in each butterfly calculation is calculated in step (3)₃A point of 2.

The image width calculated according to the imaging step is N_xRather than the width L of the object to be measured_xTherefore, the x-direction interpolation operation is required to be performed in the last step of the imaging calculation stage. The specific process is as follows: and linearly inserting inst- Δ x/Δ z-1 points between every two adjacent points in the x direction.

The image display is to display the two-dimensional image data obtained in the imaging calculation stage on a display, and a gray image or a color image can be displayed according to the requirement.

In the case that the object to be measured is a regular layered object, the boundary c (x, z) extracted by using the Canny edge extraction operator in the imaging step (5.1) is a straight line, and the imaging calculation step is unchanged.

For the case that the object to be measured only contains a single medium (described as a contact ultrasonic measurement mode, otherwise, the object is regarded as a regular layered object), z is_min＝＝z_max＝＝z_depthAnd (4) stopping the imaging calculation after the step (3) is executed.

Claims

1. The multilayer object nondestructive testing ultrasonic imaging method based on variable wave speed phase shift is characterized by sequentially comprising the following steps of:

the horizontal length of the three-layer heterogeneous object along the x-axis direction is L, the three-layer heterogeneous object is divided into L/delta x intervals, delta x is the interval length, and the interval length is from a point 0 to an end point N of the ultrasonic transducer on the x-axis_x1 step size per move, N_xThe point reached by the ultrasonic transducer in each movement is called probe point, and the total number of the probe points is N_xNumber n_x＝0，1，...，N_x-1, said positioning controller generating a TTL pulse at each detection point triggering said ultrasound transducer to transmit an excitation pulse in a depth direction z perpendicular to the x-axis of the three-layered heterogeneous object, then the transducer switches to a receiving mode and starts timing to receive the echo signal reflected from the object to be measured, said analog-to-digital converter couples said transducer at a detection point n_xProcessing the received echo signal by N_tSub-sampling and storing in computer, sampling number n_t＝0，1，...，N_t-1, sampling frequency f_sThe value is predetermined for the analog-to-digital converter, and s (z, x, t) is represented by x ═ n on the abscissa_xThe echo signal received at Δ x is n at t_t/f_sSampling values of the moments;

Subscript 0 denotes the probe number, and then repeatsThe above steps read n in sequence_x＝1，...，N_x-1 sampling values at the probe points and recording the number of first fluctuations of the sampling values at the points

From

Select the smallest serial number

And from

Select the largest serial number

Wherein v is₁Respectively making two horizontal lines parallel to the x axis according to the two depth coordinate values to surround a boundary between the first layer of medium and the second layer of medium to form a first bounding box which is positioned on the x-z longitudinal section of the three-layer heterogeneous object and simultaneously spans the first layer of medium and the second layer of medium;

and (3): the sampled values s (z, x, t) at the detection points are entered along the time axis t one by oneRow one-dimensional fast Fourier transform 1D-FFT_t(s (z, x, t)), and taking the initial value z of z₀Obtaining an x-omega domain two-dimensional signal:

S(z₀，x，ω)＝1D-FFT_t(s(z₀，x，t))

S(z₀，k，ω)＝1D-FFT_x(1D-FFT_t(s(z₀，x，t)))

0≤n_k≤N_x-1, and n_kIs an integer;

0≤n_ω≤N_t-1, and n_ωIs an integer;

wherein n is_ωThe value of the sampling angular frequency omega is 0 to (N)_t-1), take t ═ t₀One-dimensional sampled data s (z) of 0₀，x，t₀) I.e. the z-th ═ z on the vertical section of the three-layer heterogeneous object₀Pixel values of pixel points in 0 row;

Wherein,

wherein i is an imaginary unit, z₀0, ω is the sampling angular frequency, v_zThe wave velocity of the sound wave in the z-th row of medium of the measured object in the depth direction is equal to or less than z_minWhen, v_z＝v₁，v₁Is the wave velocity of the acoustic wave in the first layer of medium, and Δ z is the separation distance between two adjacent pixels in the depth direction z, and α (Δ z, k, ω, v)_z) Is the amount of phase migration at a constant wave velocity,

step (4.2): the following steps are performed in sequence:

step (4.2.2): to pairConjugating the Gaussian integral result of the step (4.2.1), performing one-dimensional fast Fourier transform on k, conjugating the calculation result, and multiplying the result by 1/N_xObtaining the coordinate value z in the depth direction in the first layer medium_minThe image formed by the row of pixel points:

s(z，x，t₀)＝s(z_min，x，t₀)，z₀＜z≤z_min；

and (5): taking v in step (4)_z＝v₁And circularly performing the following steps (5.1) and (5.2) by taking the step delta z as the step size for the z value until z is more than or equal to z_maxGenerating a first enclosureDirection of internal depth z_min＜z≤z_maxEach line of image formed by each line of pixel points in the interval respectively comprises the following steps:

S(z+Δz，k，ω)＝S(z，k，ω)·α(Δz，k，ω，v_z)，

Step (5.2): the following steps are performed in sequence:

the calculation steps are as follows:

step (10): according to the method described in step (5), v is taken_z＝v₃Z is z'_maxThe step (5.1) and the step (5.2) are circularly executed by taking the step delta z as a step size until the step z is more than or equal to z_depthGenerating the parts except the second bounding box in the third layer of mediumEach row image formed by row pixel points, v₃Is the propagation speed of the acoustic wave in the third layer medium.