CN102608205B - Multilayer-object nondestructive testing ultrasonic imaging method based on variable wave speed phase shifting - Google Patents

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

Based on the multi-layer body Non-Destructive Testing ultrasonic imaging method that becomes velocity of wave Phase shift

Technical field

The present invention relates to Ultrasonic Nondestructive technology, synthetic aperture focusing imaging technique and Phase shift technology, realize the inner structure fast imaging of the irregular layering object to complexity.

Background technology

At present, in Ultrasonic Nondestructive field, (comprise contact ultrasonic imaging and immersion method ultrasonic imaging for the ultrasonic imaging that contains interfacial layering object, if couplant layer is considered as to a part for testee, immersion method ultrasonic imaging also can be considered the contact ultrasonic imaging of layering object, therefore in the present invention, both unifications are called to the ultrasonic imaging of layering object), the current two schemes that mainly contains: the one, the method that adopts time domain ultrasonic imaging by synthetic aperture focusing technology to combine with ray trace (Ray Tracing) technology, the 2nd, adopt the Phase shift ultrasonic imaging method based on Phase shift (Phase Shift Migration) technology and frequency domain ultrasonic imaging by synthetic aperture focusing technology.

The basic thought of ultrasonic imaging by synthetic aperture focusing technology (SAFT) is to utilize a single transducer to simulate array of apertures, testee is carried out to orderly scanning, and the pulse echo signal that adopts time delay stack (DAS) method (time delay or phase delay) to obtain scanning carries out focal imaging.As shown in Figure 1, wherein, Z-direction is depth direction to SAFT working model, and X is to being transducer direction of scanning, when transducer is positioned at u _n, target reverberation (x', z') to the distance of transducer is r when at the place _n.SAFT technology has dividing of time domain and frequency domain:

Time domain SAFT, based on DAS principle and multiple spot dynamic focus technology, calculates different time delays to the focus point on the different depth of target imaging region.As Fig. 1, in order to focus at target reverberation place (x', z'), the time domain SAFT technology echoed signal that each the scanning position place in length of synthetic aperture L obtains by transducer is carried out time delay overlap-add procedure: establish s (u _n, be t) that transducer is at u _nthe echoed signal that place receives, t is sampling instant, u _nplace about the time delay of target reverberation (x', z') is

${\mathrm{\τ}}_n=\frac{2}{c}(r_n-z^{\′})=\frac{2}{c}(\sqrt{z^{\′2}+{(x^{\′}-u_n)}^2}-z^{\′}),n=\mathrm{0,1},...,L-1.$

Wherein, c is ultrasonic velocity of propagation in medium, and that locates at (x', z') is imaged as

$i(x^{\′},z^{\′})=\underset{n=0}{\overset{L-1}{\mathrm{\Σ}}}\frac{{\mathrm{\ω}}_{n}s(u_n,t-2z^{\′}/c)\mathrm{\δ}(t-{\mathrm{\τ}}_n)}{r_n}$

Wherein, ω _nfor apodizing function, δ function is unit sampling sequence.

Frequency domain SAFT is called as wave number algorithm (wavenumber algorithm), is mainly the method that generates wave beam based on phase delay, and it stems from wave equation inversion theory.Ultimate principle is that Green function (Green function) is carried out to Fourier decomposition.Algorithmic procedure is for first carrying out two-dimensional Fourier transform to ultrasound data, obtain 2-d spectrum, then utilize Nonlinear Mapping (Stolt mapping) and interpolation to realize coordinate transform, finally the frequency spectrum after conversion is done to filtering processing, and carry out two-dimensional inverse Fourier transform, the reconstructed image under the span-time domain.

Phase shift (Phase Shift Migration) ultrasonic imaging method is that the migrating technology (Migration Technique) in reflection seismology (Reflection Seismology) is combined with frequency domain SAFT, and a kind of ultrasonic imaging method obtaining.Phase shift technology is proposed in 1978 by Gazdag the earliest, and always for earthquake/radar imagery, and the Phase shift ultrasonic imaging method proposing recently has just applied to Ultrasonic NDT field by this technology.

Supersonic sounding model is considered as explosive reflection model by Phase shift ultrasonic imaging method, supposes that the reverberation in object under test exploded in the t=0 moment simultaneously, and the bursting strength of each reverberation is proportional to its reflectivity, and whole field intensity is measured with one group of transducer.Its main thought is according to extrapolating to calculate the sound field of other positions of depth direction from horizontal level (the being depth direction the first row) sound field observing.Specific algorithm comprises two main steps: the first step is carried out two-dimensional Fourier transform to time domain data, obtains 2-d spectrum,

S(k,ω)＝2D-FT(s(u _n,t))

Second step is the circulation at depth direction: first the last time 2-d spectrum obtaining that circulates done to phase shift, then makes two-dimensional inverse Fourier transform and get t=0, obtain a line image,

s(u _n,t＝0)＝2D-IFT(S(k,ω)·α(Δz,k,ω))

Wherein, for the corresponding Phase shift amount of Δ z step-length on depth direction, c is ultrasonic velocity of propagation in medium, and its value is constant, and therefore, what this ultrasonic imaging method used is the Phase shift technology of constant velocity of wave.

From the principle of DAS, time domain SAFT algorithm need to be determined time delay to the relative distance between reverberation and velocity of wave according to transducer.And the acoustic characteristics such as for multi-layer body, sound wave can reflect with the interphase place of layer at layer, reflection, thereby travel path can change.In addition, the medium that object is each layer is different, and sound wave is conventionally also unequal the velocity of propagation of each layer, therefore, the calculating of time delay can not adopt the simple computation method of relative distance divided by velocity of wave again, but need to first obtain the travel path of wave beam, then obtains piecemeal the travel-time.Its Focal point and difficult point is to find out quickly and accurately the travel path of wave beam, and Tracing Technology is realized the best approach of this object just.The principle of ray trace mainly, based on Snell theorem or Fermat principle, is found out the shortest time-consuming acoustic wave propagation path by iterative computation.The advantage of Tracing Technology is that this technology does not have special requirement to the interphase shape of multi-layer body, can be applicable to the ultrasonic imaging of arbitrary shaped body.But its serious shortcoming is that ray tracing algorithm comprises interative computation, and time complexity is high, meanwhile, DAS principle is actually a kind of convolution algorithm, and its calculated amount is also larger, and imaging time is long.As shown in Fig. 2 (b), this scheme generates this image (5cm*5cm) just needs 2 hours.This shortcoming has seriously limited promoting the use of this technology.

For the imaging of multi-layer body, Phase shift ultrasonic imaging method is without asking raypath and convolution algorithm, and image taking speed can significantly improve, and only needs 4 minutes as generated the image as shown in Fig. 2 (c).And Phase shift technology is a kind of image reconstructing method based on frequency domain, has the azimuth resolution of higher range resolution and Geng Gao and lower secondary lobe than the SAFT technology under time domain, and its image quality is higher.But, owing to supposing in Phase shift technology that the velocity of wave of sound wave is constant, make the method can only be applicable to have foreign medium on depth direction and horizontal direction must be medium of the same race, layer and the interphase between layer medium are the regular situation of level or straight line/plane parallel to each other and so on, and for the object that all has foreign medium in depth direction and horizontal direction, comprise the object of irregular interphase (curve/curved surface), due to sound wave, the velocity of wave in medium is different throughout, and cause testing result accurate not, as shown in Fig. 2 (c), there is grave error in the interfacial imaging results of the second layer.

Therefore,, aspect the ultrasonic imaging that contains irregular interfacial complicated layering object, also lack at present fast and accurately effective method.

Summary of the invention

The object of the invention is to propose a kind of ultrasonic imaging method, realize quick, accurate, effectively imaging to containing irregular interfacial complicated layering object.

The invention is characterized in, can be to regular layering object or only containing the object imaging of Single Medium, simultaneously also can be to containing non-level and the imaging of non-interfacial irregular layering object parallel to each other, and image taking speed is fast.

The invention is characterized in, adopt Phase shift ultrasonic imaging technique, and Phase shift technology is improved, in the depth direction of testee and horizontal direction, consider respectively to exist the situation of foreign medium, in Phase shift formula, introduce different velocities of wave, propose to become the Phase shift technology of velocity of wave, so that this ultrasonic imaging method is applicable to irregular layering object.One-piece construction as shown in Figure 3, in imaging process, separatrix is processed separately, outside the depth range at place, separatrix, use Phase shift formula and the Fast Fourier Transform (FFT) technology of constant velocity of wave, and within the scope of this, use the Phase shift technology and N the single output Fast Fourier Transform (FFT) technology that become velocity of wave, to reduce calculated amount.

The invention is characterized in, contain successively following steps:

Step (1): build of being formed by a computing machine, ultrasonic transducer, a set of register control and analog to digital converter based on become velocity of wave Phase shift for to the system of doing not damaged ultrasonic imaging on the vertical section forming at the degree of depth and horizontal both direction at horizontal and vertical three layers of heterogeneous object that all have a foreign medium, wherein:

Described ultrasonic transducer is provided with: the pulse signal input terminal being connected with the output terminal of described register control, the input end of described register control is connected with the corresponding positioning control signal output part of described computing machine, described ultrasonic transducer is also provided with: the echoed signal output terminal being connected with the input end of described analog to digital converter, the output terminal of described analog to digital converter is connected with the echo samples signal input part of described computing machine.Described transducer, by described register control control, moves with 1 step-length/ms fixed rate at three layers of heterogeneous body surface, and described register control is to control the gearing of described ultrasonic transducer shift position, and its parameter is by described computer input,

Three layers of heterogeneous object are L along the axial horizontal length of x, are divided into x of L/ Δ interval, and Δ x is burst length, be also described ultrasonic transducer from the point 0 x axle to terminal N _x-1 step-length only at every turn moving, N _x=1+L/ Δ x, the each mobile point reaching of described ultrasonic transducer is called sensing point, total N _xindividual, sequence number n _x=0,1 ..., N _x-1, described register control produces a TTL pulse at each sensing point place, trigger depth direction z perpendicular to the x axle transmitting driving pulse of described ultrasonic transducer to three layers of heterogeneous object, transducer transfers receiving mode to and starts timing subsequently, receive the echoed signal from testee reflection, described analog to digital converter to described transducer at sensing point n _xthe echoed signal that place receives is carried out N _tinferior sampling is also stored in computing machine, sampling sequence number n _t=0,1 ..., N _t-1, sample frequency is f _s, its value is default for analog to digital converter, and note s (z, x, t) is is x=n at horizontal ordinate _xthe echoed signal that Δ x place receives is at t=n _t/ f _sthe sampled value in moment;

Step (2): described computing machine is from n _t=0 starts sequentially to read sensing point n _xthe sampled value at=0 place, and record sampled value occur for the first time fluctuation sampling sequence number subscript 0 represents sensing point sequence number, then, repeats above-mentioned steps and reads successively n _x=1 ..., N _xthe sampled value at-1 each sensing point place, and record each point place sampled value occur for the first time fluctuation sequence number obtain N _xn altogether on individual sensing point _xthere is for the first time the sequence number of fluctuation in individual sampled value from in choose minimum sequence number and calculate this sequence number corresponding depth coordinate value in the z of described three layers of heterogeneous object direction with from in choose maximum sequence number and calculate this sequence number corresponding depth coordinate value in the z of described three layers of heterogeneous object direction wherein v ₁for the velocity of wave in ground floor medium, make two horizontal lines that are parallel to x axle by these two depth coordinate values respectively, separatrix between ground floor medium and second layer medium is surrounded, forms one and be positioned at first bounding box of simultaneously crossing over first, second two layer medium on described three layers of heterogeneous object x-z vertical section;

Step (3): along time shaft t one by one the sampled value s to each sensing point place (z, x, t) carry out one dimension Fast Fourier Transform (FFT) 1D-FFT _t(s (z, x, t)), gets the initial value z of z ₀=0, obtain x-ω territory 2D signal:

S(z ₀,x,ω)＝1D-FFT _t(s(z ₀,x,t))

Wherein, x is that sensing point sequence number is n _xthe level at place is to coordinate figure, and ω is sampling angular frequency, then, and to above-mentioned x-ω territory 2D signal S (z ₀, x, ω) and carry out successively one dimension Fast Fourier Transform (FFT) 1D-FFT along transverse axis x _x(S (z ₀, x, ω)), obtain z ₀the k-ω territory 2-d spectrum S (z at place ₀, k, ω):

S(z ₀,k,ω)＝1D-FFT _x(1D-FFT _t(s(z ₀,x,t)))

Wherein, the x that k is wave-number vector is to component, and sequence number is n _k, described wave-number vector represents the spatial gradient of wave phase,

$k=\frac{2\mathrm{\π}{n}_k}{N_x\·\mathrm{\Δx}},$ 0≤n _k≤ N _x-1, and n _kfor integer;

$\mathrm{\ω}=\frac{2\mathrm{\π}{n}_{\mathrm{\ω}}{f}_s}{N_t},$ 0≤n _ω≤ N _t-1, and n _ωfor integer;

Wherein, n _ωfor the sequence number of sampling angular frequency, value is at 0～(N _t-1), between, get t=t ₀=0, one dimension sampled data s (z ₀, x, t ₀) be z=z on three layers of heterogeneous object vertical section ₀the pixel value of the each pixel of=0 row;

Step (4): generate according to the following steps successively z direction homogeneity district z ₀<z≤z _minthe interior each row image being formed by each row pixel respectively:

Step (4.1): utilize following Phase shift formula to calculate described vertical section at z _min2-d spectrum S (the z at place _min, k, ω):

$S(z_{\min },k,\mathrm{\&omega;})=S(z_0,k,\mathrm{\&omega;})\underset{z=1}{\overset{z_{\min }}{\mathrm{\&Pi;}}}\mathrm{\&alpha;}(\mathrm{\&Delta;z},k,\mathrm{\&omega;},v_z)$

Wherein,

Wherein, i is imaginary unit, z ₀=0, ω is sampling angular frequency, v _zfor testee acoustic velocity in the capable medium of z on depth direction, at z≤z _mintime, v _z=v ₁, v ₁for the velocity of wave of sound wave in ground floor medium, Δ z is the spacing distance between upper adjacent two pixels of depth direction z, α (Δ z, k, ω, v _z) be the Phase shift amount under constant velocity of wave,

Step (4.2): carry out successively following steps:

Step (4.2.1): get t=t ₀=0, the S (z that described step (4.1) is obtained _min, k, ω) and make the Gauss integration of independent variable ω,

Step (4.2.2): the Gauss integration result to described step (4.2.1) is got conjugation, then k is done to one dimension Fast Fourier Transform (FFT), is multiplied by 1/N after result of calculation is got to conjugation again _x, obtaining coordinate figure on the interior depth direction of ground floor medium is z _minthe image that forms of one-row pixels point:

$s(z_{\min },x,t_0)=\frac{1}{N_x}{\{1D-{\mathrm{FFT}}_k\left({\left[\frac{1}{2\mathrm{\&pi;}{f}_s}{\&Integral;}_0^{2\mathrm{\&pi;}{f}_s}S(z_{\min },k,\mathrm{\&omega;})\mathrm{d\&omega;}\right]}^{*}\right)\}}^{*}$

At z direction ground floor homogeneity district z ₀<z≤z _minin, each row image is all identical:

s(z,x,t ₀)＝s(z _min,x,t ₀),z ₀<z≤z _min；

Step (5): get the v in step (4) _z=v ₁, z value is circulated taking Δ z as step-length and carries out following step (5.1), step (5.2), until z>=z _maxonly, generate depth direction z in the first bounding box _min<z≤z _maxthe interval each row image being formed by each row pixel respectively, its step is as follows:

Step (5.1): utilize the 2-d spectrum at the Phase shift formula calculating z+ Δ z place under described constant velocity of wave, z _min<z+ Δ z≤z _max:

S(z+Δz,k,ω)＝S(z,k,ω)·α(Δz,k,ω,v _z)，

Step (5.2): carry out successively following steps:

Step (5.2.1): get t=t0=0, the S (z+ Δ z, k, ω) that described step (5.1) is obtained makes the Gauss integration of independent variable ω,

Step (5.2.2): the Gauss integration result to described step (5.2.1) is got conjugation, then k is done to one dimension Fast Fourier Transform (FFT), is multiplied by 1/N after result of calculation is got to conjugation again _x, obtain the image that on depth direction, coordinate figure forms for the one-row pixels point of z+ Δ z, z _min<z+ Δ z≤z _max:

$s(z+\mathrm{\&Delta;z},x,t_0)=\frac{1}{N_x}{\{1D-{\mathrm{FFT}}_k\left({\left[\frac{1}{2\mathrm{\&pi;}{f}_s}{\&Integral;}_0^{2\mathrm{\&pi;}{f}_s}S(z_{\min },k,\mathrm{\&omega;})\mathrm{d\&omega;}\right]}^{*}\right)\}}^{*};$

Step (6): revise according to the following steps depth direction z in described the first bounding box _minto z _maxinterval z _min<z≤z _maxin each row image, the error causing to eliminate the not homogeneity between the second layer and ground floor medium:

Step (6.1): the z obtaining with step (5) _minto z _maxinterval image block, as input quantity, uses Canny operator Boundary extracting algorithm to extract described z _min, z _maxbetween ground floor medium and the separatrix c of second layer medium ₁(x, z),

Step (6.2): from z=z _minstart, until z>=z _maxonly, the z obtaining with described step (4.1) _min2-d spectrum S (the z at place _min, k, ω) and make initial value, circulation execution step (6.2.1) and the step (6.2.2) taking Δ z as step-length:

Step (6.2.1): use single output Fast Fourier Transform Inverse (FFTI) method N time, be calculated as follows the x-ω territory 2D signal that obtains z+ Δ z place in described first bounding box:

$S(z+\mathrm{\&Delta;z},x=n_x\&CenterDot;\mathrm{\&Delta;x},\mathrm{\&omega;})=\frac{1}{N_x}\underset{n_k=0}{\overset{N_x-1}{\mathrm{\&Sigma;}}}\left[S(z,k=\frac{{2\mathrm{\&pi;n}}_k}{N_x\mathrm{\&Delta;x}},\mathrm{\&omega;})\mathrm{\&alpha;}(\mathrm{\&Delta;z},k=\frac{2\mathrm{\&pi;}{n}_k}{N_x\mathrm{\&Delta;x}},\mathrm{\&omega;},v_z(x=n_x\&CenterDot;\mathrm{\&Delta;x}))\right]e^{i\frac{2\mathrm{\&pi;}}{N_x}{n}_k{n}_x}$

Wherein, the value of N equals N _x, n _x=0,1 ..., N _x-1, n _k=0,1 ..., N _x-1, α (Δ z, k, ω, v _z(x)) for becoming the Phase shift amount of velocity of wave:

Wherein, the same step of the value of k, ω (3), v _z(x) be the velocity of wave of sound wave in the medium located of coordinate (x, z):

V ₁for being positioned at separatrix c ₁sound velocity of wave propagation in ground floor medium on (x, z), v ₂for being positioned at separatrix c ₁sound velocity of wave propagation in second layer medium under (x, z),

Calculation procedure is as follows:

Step (a): respectively to different n _x, n _k(n _x=0,1 ..., N _x-1; n _k=0,1 ..., N _x-1) calculate

$G(n_x,n_k)=S(z,k=\frac{2\mathrm{\&pi;}{n}_k}{N_x\mathrm{\&Delta;x}},\mathrm{\&omega;})\mathrm{\&alpha;}(\mathrm{\&Delta;z},k=\frac{2\mathrm{\&pi;}{n}_k}{N_x\mathrm{\&Delta;x}},\mathrm{\&omega;},v_z(x=n_x\&CenterDot;\mathrm{\&Delta;x})),$

Step (b): to different n _x(n _x=0,1 ..., N _x-1), with G (n _x, n _k), n _k=0,1 ..., N _x-1 sequence is as the input quantity of cutting output fft algorithm, to n _kmake G (n _x, n _k) discrete fourier inverse transformation, but only get n _xthe value of this time domain point is the x-ω territory 2D signal at z+ Δ z place as coordinate figure on depth direction:

$S(z+\mathrm{\&Delta;z},x=n_x\&CenterDot;\mathrm{\&Delta;x},\mathrm{\&omega;})=\frac{1}{N_x}\underset{n_k=0}{\overset{N_x-1}{\mathrm{\&Sigma;}}}G(n_x,n_k)e^{i\frac{2\mathrm{\&pi;}}{N_x}{n}_k{n}_x},$

Step (6.2.2): get t=t ₀=0, the S (z+ Δ z, x, ω) that described step (6.2.1) is obtained makes the Gauss integration of independent variable ω, obtains the image that on depth direction, coordinate figure forms for the one-row pixels point of z+ Δ z, z _min<z+ Δ z<z _max:

$s(z+\mathrm{\&Delta;z},x,t_0)=\frac{1}{2\mathrm{\&pi;}{f}_s}{\text{\&Integral;}}_0^{2\mathrm{\&pi;}{f}_s}S(z_{\min },k,\mathrm{\&omega;})\mathrm{d\&omega;}$

Meanwhile, preserve S (z+ Δ z, x, ω) and obtain S (z+ Δ z, k, ω) do one-dimensional Fourier transform, as the initial value of next round iteration;

Step (7): by method step (5) described, get v _z=v ₂, z value is circulated taking Δ z as step-length and carries out described step (5.1), step (5.2), until z>=z _depthonly, generate depth direction z _max<z≤z _deptheach row image that interval each row pixel forms, v ₂for sound velocity of wave propagation in second layer medium, z _depthrepresenting the maximum coordinates value of three layers of heterogeneous object on depth direction z, is systemic presupposition value;

Step (8): the z obtaining with step (7) _maxto z _depthinterval image block, as input quantity, uses Canny operator Boundary extracting algorithm to extract described z _max, z _depthbetween second layer medium and the separatrix c of the 3rd layer of medium ₂(x, z), asks separatrix c ₂(x, z) min coordinates value z' on depth direction z _minwith maximum coordinates value z' _max, respectively by described z' _minand z' _maxtwo depth coordinate values are made two horizontal lines that are parallel to x axle, separatrix between second layer medium and the 3rd layer of medium is surrounded, forms one and be positioned at second bounding box simultaneously crossing over second, third two layer medium on described three layers of heterogeneous object x-z vertical section;

Step (9): according to depth direction z' in method correction the second bounding box step (6.2) Suo Shu _minto z' _maxinterval z' _min<z≤z' _maxin each row image, with the error of eliminating the not homogeneity between the 3rd layer and second layer medium and cause;

Step (10): by method step (5) described, get v _z=v ₃, to z value with z' _maxfor initial value, taking Δ z as step-length, described step (5.1), step (5.2) are carried out in circulation, until z>=z _depthonly, generate each row image that in the 3rd layer of medium, the each row pixel of part forms except described the second bounding box, v ₃it is sound velocity of wave propagation in the 3rd layer of medium.

The invention has the advantages that, imaging precision is high, for the testee shown in Fig. 2 (a), if transducer diameter is 0.5mm, transducer moving step length is 0.5mm, and the ultrasound wave centre frequency of transducer transmitting is 5MHz, sample frequency 100MHz, the step delta z of depth direction is taken as 0.05mm, and ground floor interphase error can be controlled in 0.4mm, and second layer interphase can be controlled at 1.0mm with interior (as Fig. 4).

The present invention has following advantage compared with prior art:

1. can carry out fast ultrasonic imaging to containing irregular interfacial object;

2. the resolution of imaging and precision are higher.

Brief description of the drawings

Fig. 1 is the working model of SAFT: X to being transducer direction of scanning, and Z-direction is depth direction, u _nfor transducer with etc. the position of step-length while moving, be also the position of sensing point, (x ', z ') is the coordinate of target reverberation, r _nfor target reverberation is to the distance of transducer.

Fig. 2 is the sectional view that the comparison diagram of irregular layering object ultrasonic imaging method under identical experiment environment: 2 (a) are testee; 2 (b) are that time domain SAFT is combined the image that 2 (a) generated with Tracing Technology, and imaging time is more than 2 hours, and very responsive to the parameter of algorithm, and robustness is not fine; The image that the Phase shift ultrasonic imaging technique that 2 (c) are constant velocity of wave generates, imaging time 4 minutes, its error is larger; The image that 2 (d) generate for this formation method, imaging time 6 minutes.

Fig. 3 is this formation method overall process schematic diagram.

Fig. 4 is the analysis of experimental data figure to Fig. 2 (d): 4 (a) are ground floor separatrix Error Graph, solid-line curve is the marginal standard value of testee ground floor, imaginary point curve is the ground floor separatrix that this formation method obtains, the marginal graph of errors of ground floor that dashed curve obtains for this formation method; 4 (b) are second layer separatrix Error Graph, solid-line curve is the marginal standard value of the testee second layer, imaginary point curve is the second layer separatrix that this formation method obtains, the marginal graph of errors of the second layer that dashed curve obtains for this formation method.

Fig. 5 is this ultrasonic image-forming system schematic flow sheet.

Fig. 6 is this ultrasonic imaging hardware system structure figure.

Fig. 7 is ultrasonic transducer work schematic diagram.

Fig. 8 is the multi-layer body Non-Destructive Testing ultrasonic imaging algorithm flow chart based on becoming velocity of wave Phase shift technology.

Fig. 9 asks bounding box schematic diagram in this ultrasonic imaging method.

FFT butterfly computation schematic diagram when Figure 10 is N=8.

Embodiment

Specific embodiment of the invention process comprises three parts (as Fig. 5): ultrasound data obtains, imaging calculating and image show.Hardware platform system structural drawing as shown in Figure 6, ultrasonic image-forming system is made up of a computing machine, ultrasonic transducer, a set of register control and an analog to digital converter, the pulse signal input terminal of ultrasonic transducer is connected with the output terminal of register control, and the input end of register control is connected with the positioning control signal output part of computing machine.The echoed signal output terminal of ultrasonic transducer is connected with the input end of analog to digital converter, and the output terminal of analog to digital converter is connected with the echo samples signal input part of computing machine.

The mode that ultrasound data obtains can be selected: contact, the i.e. surface of the direct contact measured object of ultrasonic transducer; Or liquid immersion type, object under test is placed in liquid, and ultrasonic transducer is measured at liquid surface.According to the difference of detection mode, transducer can use contact probe or liquid immersion type probe.System is used single transmitting/receiving transducer, transducer by register control in testee surface or liquid with uniform step delta x in the x-direction (as Fig. 7) move with about 1 step/ms fixed rate.The controller moment stable in each target location produces a TTL(transistor-transistor logic level) pulse, this pulse is used for triggering ultrasonic transducer to testee and the perpendicular driving pulse of depth direction transmitting of x direction, transducer transfers receiving mode to and starts timing subsequently, receives from the echo of testee reflection.Each position of transducer transponder pulse and reception echo is sensing point.The echoed signal that transducer receives is gathered and is stored in storer by analog to digital converter.Δ x need require determine according to the actual size of object under test and imaging precision, and its value is less, and the image of generation is more accurate, but computing time is also longer.

It is exactly that sampled data using testee each sensing point place on a vertical section is as computer input that imaging is calculated, then calculate the skiagraph picture of testee by following image-forming step (process flow diagram is referring to Fig. 8), wherein s (z, x, t) sampled value of the signal that receives in the t moment of the transducer that is illustrated in horizontal direction x place; Z represents the coordinate figure of depth direction; Δ z represents the step-length (relevant with precision) of depth direction, and the image that generated is in the spacing of longitudinal (that is depth direction) upper adjacent two pixels; z _depthfor default image maximum depth value in the vertical:

Step (1): ground floor medium is isomorphism medium, and its ultrasonic reflection signal is identical, and at ground floor and second layer boundary, ultrasonic reflection signal can produce fluctuation., according to the sampled data of ultrasound echo signal, search for each sensing point place the moment of signal fluctuation occurs for the first time, and find out minimal instant value t from these fluctuation moment _minwith maximum moment value t _max, then according to the velocity of wave v in ground floor medium ₁ask marginal bounding box between ground floor and second layer medium, obtain two boundary line: the zs of bounding box on depth direction (be z to) _min=v ₁t _minand z _max=v ₁t _max(as Fig. 9);

Step (2): get depth direction initial value z=z ₀=0, parameter x and t are carried out to two-dimensional Fourier transform, original ultrasound echo signal is transformed to k-ω territory signal:

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

Wherein, ω for sampling angular frequency, k be that the x of wave-number vector (wavenumber vector) is to component (wave-number vector refers to the spatial gradient of wave phase);

Step (3): generate depth direction z ₀to z _mininterval (z ₀≤ z≤z _min) image.Because the medium of testee in this section of interval is identical, therefore this section of image is also just the same:

Step (3.1): utilize Phase shift formula to calculate z _minthe k-ω territory signal at place:

$S(z_{\min },k,\mathrm{\&omega;})=S(z_0,k,\mathrm{\&omega;})\underset{z=1}{\overset{z_{\min }}{\mathrm{\&Pi;}}}\mathrm{\&alpha;}(\mathrm{\&Delta;z},k,\mathrm{\&omega;},v_z)---\left(2\right)$

Wherein,

I is imaginary unit, v _zfor testee acoustic velocity in the capable medium of z on depth direction, as z≤z _mintime, v _z=v ₁, v ₁for the velocity of wave of sound wave in ground floor medium;

Step (3.2): make two-dimensional inverse Fourier transform, make t=t ₀=0, obtaining along slope coordinate value is a line image of z:

$s(z,x,t_0)=s(z_{\min },x,t_0)=\frac{1}{4{\mathrm{\&pi;}}^2}\&Integral;\&Integral;S(z_{\min },k,\mathrm{\&omega;})e^{\mathrm{ikx}}{e}^{\mathrm{i\&omega;t}}\mathrm{dkd\&omega;},(z_0\&le;z\&le;z_{\min })---\left(3\right)$

Step (4): generate depth direction z _minto z _maxinterval (z _min<z≤z _max) image.Get v _z=v ₁, to z value with Δ z step-length circulation execution step (4.1) and step (4.2) until z>=z _max:

Step (4.1): utilize Phase shift formula to calculate the k-ω territory signal at z+ Δ z place:

S(z+Δz,k,ω)＝S(z,k,ω)α(Δz,k,ω,v _z) (4)

Step (4.2): make z ← z+ Δ z, t=t ₀=0, make two-dimensional inverse Fourier transform, obtaining along slope coordinate value is a line image of z:

$s(z,x,t_0)=\frac{1}{4{\mathrm{\&pi;}}^2}\&Integral;\&Integral;S(z,k,\mathrm{\&omega;})e^{\mathrm{ikx}}{e}^{\mathrm{i\&omega;t}}\mathrm{dkd\&omega;}---\left(5\right)$

Step (5): revise z _minto z _maxinterval (z _min<z≤z _max) image.Because the calculating in step (4) has only been used the velocity of wave of ground floor, be still considered as medium of the same race by this interval, the result of calculation obtaining certainly exists error, so this step needs to revise.Specifically comprise two steps:

Step (5.1): the z that step (4) is obtained _minto z _maxinterval image, uses Canny arithmetic operators to extract z _minwith z _maxbetween separatrix c (x, z), z _min≤ z≤z _max.The computation process of Canny arithmetic operators is divided into four steps:

Step (a): image smoothing.Original image is obtained to image I (x, z) with two-dimensional Gaussian function smothing filtering;

Step (b): the gradient M of the each pixel of computed image I (x, z) and direction Q.Adopt the first approximation of 2 × 2 templates as the partial differential to x direction and z direction,

$p=\frac{1}{2}\left(\begin{array}{cc}-1& 1\\ -1& 1\end{array}\right),q=\frac{1}{2}\left(\begin{array}{cc}1& 1\\ -1& -1\end{array}\right)$

Gradient magnitude M and direction Q are

$M=\sqrt{p\&times;p+q\&times;q},Q=\arctan \left(\frac{q}{p}\right);$

Step (c): the non-maximum value of gradient image suppresses.To each pixel I (x, z), if its Grad M is (x, z) be less than along the gradient magnitude of two consecutive point of its gradient direction Q (x, z), illustrate that this point is not marginal point, the gray scale of I (x, z) is made as to 0;

Step (d): dual threashold value processing.Set dual threshold method and detect and the low threshold value Low and the high threshold High that are connected edge needs, gradient image is carried out to dual threashold value processing.What gradient magnitude was greater than high threshold High is edge; What gradient magnitude was less than low threshold value Low is not edge; Gradient magnitude is marginal, judges in eight neighbors of this pixel whether have the edge pixel that is greater than high threshold High, is edge pixel, otherwise is not if exist;

Step (5.2): to z value from z=z _minstart, with Δ z step-length circulation execution step (5.2.1) and step (5.2.2) until z>=z _max:

Step (5.2.1): the k-ω territory signal that utilizes the Phase shift formula calculating z+ Δ z place that becomes velocity of wave:

S(z+Δz,k,ω)＝S(z,k,ω)α(Δz,k,ω,v _z(x)) (6)

Wherein, α (Δ z, k, ω, v _z(x)) for becoming the Phase shift amount of velocity of wave:

V _z(x) value is as follows:

V ₁, v ₂equal respectively sound velocity of wave propagation in the two layer medium up and down at place, separatrix c (x, z);

Step (5.2.2): make z ← z+ Δ z, t=t ₀=0, make two-dimensional inverse Fourier transform, obtaining along slope coordinate value is a line image of z:

$s(z,x,t_0)=\frac{1}{4{\mathrm{\&pi;}}^2}\&Integral;\&Integral;S(z,k,\mathrm{\&omega;})e^{\mathrm{ikx}}{e}^{\mathrm{i\&omega;t}}\mathrm{dkd\&omega;}---\left(7\right)$

Step (6): generate depth direction z _maxto z _depthinterval (z _max<z≤z _depth) image.Make v _z=v ₂, calculate z by the method for step (4) _maxto z _depthinterval each row image;

Step (7): process follow-up separatrix.Image section (the z that step (6) is generated _max<z≤z _depth) use Canny arithmetic operators to extract new separatrix c (x, z), make v ₁, v ₂be respectively sound velocity of wave propagation in the two layer medium up and down at place, separatrix c (x, z), and make z _minand z _maxbe respectively the up-and-down boundary of new separatrix c (x, z) on depth direction, form new bounding box, if z _min==z _max==z _depth, finish; Otherwise execution step (5.2) is revised z _minto z _maxinterval image;

Step (8): if z _max==z _depth, finish; Otherwise, repeated execution of steps (6) and step (7).

Wherein, for the two-dimensional Fourier transform in performing step (2) (as formula (1)), the present invention adopts the method that parameter x and t are carried out to twice one dimension Fast Fourier Transform (FFT) (1D-FFT) in order, and to improve computing velocity, formula (1) is embodied as:

S(z ₀,k,ω)＝1D-FFT _x(1D-FFT _t(s(z ₀,x,t)))

In step (3.2) and step (4.2), need to do two-dimensional inverse Fourier transform and get t=t ₀=0 value, can shilling t=t as view data ₀=0, and ω is made to integration, and then k is made to integration, to reduce the data volume of the second dimension integral operation.Formula (3) (5) turns to so:

$s(z,x,t_0)=\frac{1}{2\mathrm{\&pi;}}\&Integral;\{\frac{1}{2\mathrm{\&pi;}}\&Integral;S(z,k,\mathrm{\&omega;})\mathrm{d\&omega;}\}{e}^{\mathrm{ikx}}\mathrm{dk}$

Wherein, the integration to ω can utilize numerical integration method (as Gauss Gaussian integration) to try to achieve fast.And the Fourier inversion of k is still realized by one dimension Fast Fourier Transform (FFT) (1D-FFT), the data that just first numerical integration obtained are got conjugation ([Integration _ω(S (z, k, ω))] *), then directly utilize 1D-FFT subroutine to carry out computing, then operation result got to conjugation one time, and be multiplied by 1/N _x, can obtain.N _xfor transducer moves step number, i.e. N _x=1+L _x/ Δ x, wherein L _xfor object under test is in the length of x direction.Therefore, the specific implementation procedural representation of formula (3) (5) is:

$s(z,x,t_0)=\frac{1}{N_x}{\{1D-{\mathrm{FFT}}_k\left({\left[\mathrm{Integratio}{n}_{\mathrm{\&omega;}}\left(S(z,k,\mathrm{\&omega;})\right)\right]}^{*}\right)\}}^{*}$

Phase shift amount α (Δ z, k, ω, v _z) need to be for each speed v _z, respectively different k, ω to be calculated, its result is kept in three-dimensional array in advance.Wherein, the value of k, ω is respectively:

$k=\frac{2\mathrm{\&pi;}{n}_k}{N_x\&CenterDot;\mathrm{\&Delta;x}},$ 0≤n _k≤ N _x-1, and n _kfor integer;

$\mathrm{\&omega;}=\frac{2\mathrm{\&pi;}{n}_{\mathrm{\&omega;}}{f}_s}{N_t},$ 0≤n _ω≤ N _t-1, and n _ωfor integer;

Wherein, Δ x is the step-length that transducer moves in the x-direction, f _sfor sample frequency (its value is systemic presupposition), N _x=1+L _x/ Δ x is that transducer detection is in the x-direction counted, L _x/ Δ x is that transducer moves step number, L _xfor the length (x direction) of testee, N _tfor transducer is at each sampling number of sensing point place to echoed signal.Because its working model is considered as explosive reflection model by Phase shift ultrasonic imaging technique, only consider that sound wave is from the reverberation part that upwards (along z axle negative direction) propagated, so at Phase shift amount α (Δ z, k, ω, v _z) calculating in only get π f _s≤ ω <2 π f _sthe value of part, and make 0≤ω < π f _sthe Phase shift amount of part is 0, if while b<0 gets

In step (5.2), become velocity of wave Phase shift amount α (Δ z, k, ω, v _z(x)) be the function of x, velocity of wave v _znon-constant, so step (5.2.1) can not directly be calculated again, and need to merge in the Fourier inversion of step (5.2.2),

$s(z+\mathrm{\&Delta;z},x,t_0)=\frac{1}{2\mathrm{\&pi;}}\&Integral;\{\frac{1}{2\mathrm{\&pi;}}\&Integral;[S(z,k,\mathrm{\&omega;})\mathrm{\&alpha;}(\mathrm{\&Delta;z},k,\mathrm{\&omega;},v_z\left(x\right))\left]e^{\mathrm{ikx}}\mathrm{dk}\right\}{e}^{\mathrm{i\&omega;t}}\mathrm{d\&omega;}---\left(8\right)$

Outer integration in formula (8) can utilize numerical integration method (as Gauss Gaussian integration) to try to achieve, internal layer integral operation (suc as formula (9)) needs to try to achieve (suc as formula (10)) by amendment discrete fourier inverse transformation computing formula

$S(z+\mathrm{\&Delta;z},x,\mathrm{\&omega;})=\frac{1}{2\mathrm{\&pi;}}\&Integral;\left[S(z,k,\mathrm{\&omega;})\mathrm{\&alpha;}(\mathrm{\&Delta;z},k,\mathrm{\&omega;},v_z\left(x\right))\right]e^{\mathrm{ikx}}\mathrm{dk}---\left(9\right)$

$S(z+\mathrm{\&Delta;z},x=n_x\&CenterDot;\mathrm{\&Delta;x},\mathrm{\&omega;})=\frac{1}{N_x}\underset{n_k=0}{\overset{N_x-1}{\mathrm{\&Sigma;}}}\left[S(z,k=\frac{2\mathrm{\&pi;}{n}_k}{N_x\mathrm{\&Delta;x}},\mathrm{\&omega;})\mathrm{\&alpha;}(\mathrm{\&Delta;z},k=\frac{2\mathrm{\&pi;}{n}_k}{N_x\mathrm{\&Delta;x}},\mathrm{\&omega;},v_z\left(n_x\right))\right]e^{i\frac{2\mathrm{\&pi;}}{N_x}{n}_k{n}_x}---\left(10\right)$

In formula (10), n _xfor the sensing point at the transducer place sequence number in x direction, span is same as n _k, and α (Δ z, k, ω, v _z(x) be) that independent variable is n _kand n _xfunction, can not directly use discrete fourier inverse transformation (IDFT) or Fast Fourier Transform Inverse (FFTI) (IFFT) method.Therefore, in the present invention, use N single output Fast Fourier Transform Inverse (FFTI) method (NsIFFTPruning).The thought of the method based on butterfly computation in Fast Fourier Transform (FFT), but a value in transform domain is only exported in conversion each time.NsIFFTPruning specific algorithm is as follows:

Forn _x=0toN _x-1

Forn _k=0toN _x-1

$G(n_x,n_k)=S(z,k=\frac{2\mathrm{\&pi;}{n}_k}{N_x\mathrm{\&Delta;x}},\mathrm{\&omega;})\mathrm{\&alpha;}(\mathrm{\&Delta;z},k=\frac{2\mathrm{\&pi;}{n}_k}{N_x\mathrm{\&Delta;x}},\mathrm{\&omega;},v_z\left(n_x\right))$

end

S(z+Δz,x＝n _x·Δx,ω)=IFFTPruning(G,N _x,n _x)

end

Wherein, work as n _xwhen fixing, function IFFTPruning (G, N _x, n _x) calculating G (n _x, n _k), n _k=0,1 ..., N _xthe discrete fourier inverse transformation of-1 sequence, but only get n _xthe value of this time domain point,

$S(z+\mathrm{\&Delta;z},x=n_x\&CenterDot;\mathrm{\&Delta;x},\mathrm{\&omega;})=\frac{1}{N}\underset{n_k=0}{\overset{N_x-1}{\mathrm{\&Sigma;}}}G(n_x,n_k)e^{i\frac{2\mathrm{\&pi;}}{N_x}{n}_k{n}_x}---\left(11\right)$

Formula (11) can be calculated with IDFT or IFFT algorithm, but different from general IDFT problem: formula (11) only needs to calculate G (n _x, n _k) at n _xthe value at this time domain point place of=h, if directly use IDFT or IFFT algorithm, certainly exists a large amount of redundant computation, for example, and in Figure 10, if only need calculate n _x=2 point, in figure, the butterfly computation of dotted portion is redundancy.In order to reduce calculated amount, the present invention uses the fft algorithm (IFFTPruning) of cutting output.This algorithm, still based on fft algorithm, adopts butterfly computation method, and its principle is as follows: according to butterfly group technology in fft algorithm, calculating is divided into H step (Stage), wherein H meets 2 ^h=N _x.Step h (h=1 ..., H) in, by N _xpoint list entries G (n _x, n _k), n _k=0,1 ..., N _x-1 is equally divided into 2 ^h-hindividual group, increase progressively 1 numbering b since 0 respectively in each group _h(b _h=0,1 ..., 2 ^h-1), in every group of butterfly computation except being numbered b _h=b _h+1mod2 ^h-1point outside, remaining point is all redundancy.In each calculation procedure, only calculate in each group of butterfly computation and be numbered b _hpoint.As shown in figure 10, N _x=8, n _x=2 o'clock, in step (1), calculate in every group of butterfly computation and be numbered b ₁=0 point, calculates in every group of butterfly computation and is numbered b in step (2) ₂=2 point, calculates in every group of butterfly computation and is numbered b in step (3) ₃=2 point.

The picture traverse calculating according to aforementioned image-forming step is N _x, but not the wide L of object under test _x, therefore, need to carry out x to interpolation arithmetic in the final step of imaging calculation stages.Detailed process is: insert z-1 point of inst=Δ x/ Δ at x to linearity between every adjacent two points.

Image shows that the two-dimensional image data obtaining by imaging calculation stages is presented on display, display gray scale image or coloured image as required.

For to be measured, to liking the situation of regular layering object, the separatrix c (x, z) that in image-forming step (5.1), use Canny arithmetic operators extracts is straight line, and imaging calculation procedure is constant.

For the situation (be illustrated as contact ultrasonic measurement mode, otherwise should be considered as regular layering object) that only contains Single Medium in object to be measured, z _min==z _max==z _depth, imaging stops after calculating execution of step (3).

Claims (1)

1. the multi-layer body Non-Destructive Testing ultrasonic imaging method based on becoming velocity of wave Phase shift, is characterized in that, contains successively following steps:

Step (1): build of being formed by a computing machine, ultrasonic transducer, a set of register control and analog to digital converter based on become velocity of wave Phase shift for to the system of doing not damaged ultrasonic imaging on the vertical section forming at the degree of depth and horizontal both direction at horizontal and vertical three layers of heterogeneous object that all have a foreign medium, wherein:

Described ultrasonic transducer is provided with: the pulse signal input terminal being connected with the output terminal of described register control; the input end of described register control is connected with the corresponding positioning control signal output part of described computing machine; described ultrasonic transducer is also provided with: the echoed signal output terminal being connected with the input end of described analog to digital converter, and the output terminal of described analog to digital converter is connected with the echo samples signal input part of described computing machine; Described transducer, by described register control control, moves with 1 step-length/ms fixed rate at three layers of heterogeneous body surface, and described register control is to control the gearing of described ultrasonic transducer shift position, and its parameter is by described computer input,

Three layers of heterogeneous object are L along the axial horizontal length of x, are divided into x of L/ Δ interval, and Δ x is burst length, be also described ultrasonic transducer from the point 0 x axle to terminal N _x-1 step-length only at every turn moving, N _x=1+L/ Δ x, the each mobile point reaching of described ultrasonic transducer is called sensing point, total N _xindividual, sequence number n _x=0,1 ..., N _x-1, described register control produces a TTL pulse at each sensing point place, trigger depth direction z perpendicular to the x axle transmitting driving pulse of described ultrasonic transducer to three layers of heterogeneous object, transducer transfers receiving mode to and starts timing subsequently, receive the echoed signal from testee reflection, described analog to digital converter to described transducer at sensing point n _xthe echoed signal that place receives is carried out N _tinferior sampling is also stored in computing machine, sampling sequence number n _t=0,1 ..., N _t-1, sample frequency is f _s, its value is default for analog to digital converter, and note s (z, x, t) is is x=n at horizontal ordinate _xthe echoed signal that Δ x place receives is at t=n _t/ f _sthe sampled value in moment;

Step (2): described computing machine is from n _t=0 starts sequentially to read sensing point n _xthe sampled value at=0 place, and record sampled value occur for the first time fluctuation sampling sequence number subscript 0 represents sensing point sequence number, then, repeats above-mentioned steps and reads successively n _x=1 ..., N _xthe sampled value at-1 each sensing point place, and record each point place sampled value occur for the first time fluctuation sequence number obtain N _xn altogether on individual sensing point _xthere is for the first time the sequence number of fluctuation in individual sampled value from in choose minimum sequence number and calculate this sequence number corresponding depth coordinate value in the z of described three layers of heterogeneous object direction with from in choose maximum sequence number and calculate this sequence number corresponding depth coordinate value in the z of described three layers of heterogeneous object direction wherein v ₁for the velocity of wave in ground floor medium, make two horizontal lines that are parallel to x axle by these two depth coordinate values respectively, separatrix between ground floor medium and second layer medium is surrounded, forms one and be positioned at first bounding box of simultaneously crossing over first, second two layer medium on described three layers of heterogeneous object x-z vertical section;

Step (3): along time shaft t one by one the sampled value s to each sensing point place (z, x, t) carry out one dimension Fast Fourier Transform (FFT) 1D-FFT _t(s (z, x, t)), gets the initial value z of z ₀=0, obtain x-ω territory 2D signal:

S(z ₀,x,ω)＝1D-FFT _t(s(z ₀,x,t))

Wherein, x is that sensing point sequence number is n _xthe level at place is to coordinate figure, and ω is sampling angular frequency, then, and to above-mentioned x-ω territory 2D signal S (z ₀, x, ω) and carry out successively one dimension Fast Fourier Transform (FFT) 1D-FFT along transverse axis x _x(S (z ₀, x, ω)), obtain z ₀the k-ω territory 2-d spectrum S (z at place ₀, k, ω):

S(z ₀,k,ω)＝1D-FFT _x(1D-FFT _t(s(z ₀,x,t)))

Wherein, the x that k is wave-number vector is to component, and sequence number is n _k, described wave-number vector represents the spatial gradient of wave phase,

0≤n _k≤ N _x-1, and n _kfor integer;

0≤n _ω≤ N _t-1, and n _ωfor integer;

Wherein, n _ωfor the sequence number of sampling angular frequency, value is at 0～(N _t-1), between, get t=t ₀=0, one dimension sampled data s (z ₀, x, t ₀) be z=z on three layers of heterogeneous object vertical section ₀the pixel value of the each pixel of=0 row;

Step (4): generate according to the following steps successively z direction homogeneity district z ₀<z≤z _minthe interior each row image being formed by each row pixel respectively:

Step (4.1): utilize following Phase shift formula to calculate described vertical section at z _min2-d spectrum S (the z at place _min, k, ω):

Wherein,

Wherein, i is imaginary unit, z ₀=0, ω is sampling angular frequency, v _zfor testee acoustic velocity in the capable medium of z on depth direction, at z≤z _mintime, v _z=v ₁, v ₁for the velocity of wave of sound wave in ground floor medium, Δ z is the spacing distance between upper adjacent two pixels of depth direction z, α (Δ z, k, ω, v _z) be the Phase shift amount under constant velocity of wave,

Step (4.2): carry out successively following steps:

Step (4.2.1): get t=t ₀=0, the S (z that described step (4.1) is obtained _min, k, ω) and make the Gauss integration of independent variable ω,

Step (4.2.2): the Gauss integration result to described step (4.2.1) is got conjugation, then k is done to one dimension Fast Fourier Transform (FFT), is multiplied by 1/N after result of calculation is got to conjugation again _x, obtaining coordinate figure on the interior depth direction of ground floor medium is z _minthe image that forms of one-row pixels point:

At z direction ground floor homogeneity district z ₀<z≤z _minin, each row image is all identical:

s(z,x,t ₀)＝s(z _min,x,t ₀),z ₀<z≤z _min；

Step (5): get the v in step (4) _z=v ₁, z value is circulated taking Δ z as step-length and carries out following step (5.1), step (5.2), until z>=z _maxonly, generate depth direction z in the first bounding box _min<z≤z _maxthe interval each row image being formed by each row pixel respectively, its step is as follows:

Step (5.1): utilize the 2-d spectrum at the Phase shift formula calculating z+ Δ z place under described constant velocity of wave, z _min<z+ Δ z≤z _max:

S(z+Δz,k,ω)＝S(z,k,ω)·α(Δz,k,ω,v _z)，

Step (5.2): carry out successively following steps:

Step (5.2.1): get t=t ₀=0, the S (z+ Δ z, k, ω) that described step (5.1) is obtained makes the Gauss integration of independent variable ω,

Step (5.2.2): the Gauss integration result to described step (5.2.1) is got conjugation, then k is done to one dimension Fast Fourier Transform (FFT), is multiplied by 1/N after result of calculation is got to conjugation again _x, obtain the image that on depth direction, coordinate figure forms for the one-row pixels point of z+ Δ z, z _min<z+ Δ z≤z _max:

Step (6): revise according to the following steps depth direction z in described the first bounding box _minto z _maxinterval z _min<z≤z _maxin each row image, the error causing to eliminate the not homogeneity between the second layer and ground floor medium:

Step (6.1): the z obtaining with step (5) _minto z _maxinterval image block, as input quantity, uses Canny operator Boundary extracting algorithm to extract described z _min, z _maxbetween ground floor medium and the separatrix c of second layer medium ₁(x, z),

Step (6.2): from z=z _minstart, until z>=z _maxonly, the z obtaining with described step (4.1) _min2-d spectrum S (the z at place _min, k, ω) and make initial value, circulation execution step (6.2.1) and the step (6.2.2) taking Δ z as step-length:

Step (6.2.1): use single output Fast Fourier Transform Inverse (FFTI) method N time, be calculated as follows the x-ω territory 2D signal that obtains z+ Δ z place in described first bounding box:

Wherein, the value of N equals N _x, n _x=0,1 ..., N _x-1, n _k=0,1 ..., N _x-1, α (Δ z, k, ω, v _z(x)) for becoming the Phase shift amount of velocity of wave:

Wherein, the same step of the value of k, ω (3), v _z(x) be the velocity of wave of sound wave in the medium located of coordinate (x, z):

V ₁for being positioned at separatrix c ₁sound velocity of wave propagation in ground floor medium on (x, z), v ₂for being positioned at separatrix c ₁sound velocity of wave propagation in second layer medium under (x, z),

Calculation procedure is as follows:

Step (a): respectively to different n _x, n _k(n _x=0,1 ..., N _x-1; n _k=0,1 ..., N _x-1) calculate

Step (b): to different n _x(n _x=0,1 ..., N _x-1), with G (n _x, n _k), n _k=0,1 ..., N _x-1 sequence is as the input quantity of cutting output fft algorithm, to n _kmake G (n _x, n _k) discrete fourier inverse transformation, but only get n _xthe value of this time domain point is the x-ω territory 2D signal at z+ Δ z place as coordinate figure on depth direction:

Step (6.2.2): get t=t ₀=0, the S (z+ Δ z, x, ω) that described step (6.2.1) is obtained makes the Gauss integration of independent variable ω, obtains the image that on depth direction, coordinate figure forms for the one-row pixels point of z+ Δ z, z _min<z+ Δ z<z _max:

Meanwhile, preserve S (z+ Δ z, x, ω) and obtain S (z+ Δ z, k, ω) do one-dimensional Fourier transform, as the initial value of next round iteration;

Step (7): by method step (5) described, get v _z=v ₂, z value is circulated taking Δ z as step-length and carries out described step (5.1), step (5.2), until z>=z _depthonly, generate depth direction z _max<z≤z _deptheach row image that interval each row pixel forms, v ₂for sound velocity of wave propagation in second layer medium, z _depthrepresenting the maximum coordinates value of three layers of heterogeneous object on depth direction z, is systemic presupposition value;

Step (8): the z obtaining with step (7) _maxto z _depthinterval image block, as input quantity, uses Canny operator Boundary extracting algorithm to extract described z _max, z _depthbetween second layer medium and the separatrix c of the 3rd layer of medium ₂(x, z), asks separatrix c ₂(x, z) min coordinates value z' on depth direction z _minwith maximum coordinates value z' _max, respectively by described z' _minand z' _maxtwo depth coordinate values are made two horizontal lines that are parallel to x axle, separatrix between second layer medium and the 3rd layer of medium is surrounded, forms one and be positioned at second bounding box simultaneously crossing over second, third two layer medium on described three layers of heterogeneous object x-z vertical section;

Step (9): according to depth direction z' in method correction the second bounding box step (6.2) Suo Shu _minto z' _maxinterval z' _min<z≤z' _maxin each row image, with the error of eliminating the not homogeneity between the 3rd layer and second layer medium and cause;

Step (10): by method step (5) described, get v _z=v ₃, to z value with z' _maxfor initial value, taking Δ z as step-length, described step (5.1), step (5.2) are carried out in circulation, until z>=z _depthonly, generate each row image that in the 3rd layer of medium, the each row pixel of part forms except described the second bounding box, v ₃it is sound velocity of wave propagation in the 3rd layer of medium.