CN102608205A - 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 the 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 inner structure fast imaging the non-regular layering object of complicacy.

Background technology

At present; In the Ultrasonic Nondestructive field; (comprise contact ultrasonic imaging and immersion method ultrasonic imaging for the ultrasonic imaging that contains interfacial layering object; If the couplant layer is regarded as the part of testee; Then the immersion method ultrasonic imaging also can be considered the contact ultrasonic imaging of layering object; Therefore among the present invention with both unified ultrasonic imagings that is called the layering object), currently mainly contain two kinds of schemes: 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 based on Phase shift (Phase Shift Migration) technology and the technological Phase shift ultrasonic imaging method of frequency domain ultrasonic imaging by synthetic aperture focusing.

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 orderly scanning, and adopt time-delay stack (DAS) method (time delay or phase delay) that the pulse echo signal that scanning obtains is carried out focal imaging.The SAFT working model is as shown in Figure 1, and wherein, Z is to being depth direction, and X is to being the transducer direction of scanning, when transducer is positioned at u _nDuring the place, the distance of target reflection thing (x ', z ') to transducer is r _nThe SAFT technology has the branch of time domain and frequency domain:

Time domain SAFT calculates different time delays based on DAS principle and multiple spot dynamic focus technology to the focus point on the different depth of target imaging zone.Like Fig. 1, for (x ', z ') focuses at target reflection thing place, the echoed signal that time domain SAFT technology obtains transducer each scanning position place in the length of synthetic aperture L overlap-add procedure of delaying time: establish s (u _n, be that transducer is at u t) _nThe echoed signal that the place receives, t is sampling instant, u _nThe place about the time-delay of target reflection thing (x ', z ') does

${\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, then is imaged as what (x ', z ') located

$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, ω _nBe apodizing function, the δ function is the unit sampling sequence.

Frequency domain SAFT is called as wave number algorithm (wavenumber algorithm), mainly is based on the method that phase delay generates wave beam, and it stems from the wave equation inversion theory.Ultimate principle is that Green function (Green function) is carried out Fourier decomposition.Algorithmic procedure is for carrying out two-dimensional Fourier transform to ultrasound data earlier; Obtain 2-d spectrum; Utilize Nonlinear Mapping (Stolt mapping) and interpolation to realize coordinate transform then; At last the frequency spectrum after the conversion is done 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 the reflection seismology (Reflection Seismology) is combined with frequency domain SAFT, and a kind of ultrasonic imaging method that obtains.The Phase shift technology is proposed in 1978 by Gazdag the earliest, is used for earthquake/radar imagery always, and the Phase shift ultrasonic imaging method that proposes has recently just applied to the Ultrasonic NDT field with this technology.

The Phase shift ultrasonic imaging method is regarded as the explosive reflection model with the supersonic sounding model, supposes that the reverberation in the object under test explodes at t=0 constantly 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 to calculate the sound field of other positions of depth direction according to the sound field extrapolation that observes from horizontal level (being depth direction first row).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 was the circulation at depth direction: the 2-d spectrum that earlier the last time circulation is obtained is done phase shift, makes two-dimensional inverse Fourier transform then and gets t=0, obtains delegation's image,

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

Wherein,

is the pairing Phase shift amount of Δ z step-length on the depth direction; C is ultrasonic velocity of propagation in medium; Its value is constant, and therefore, what this ultrasonic imaging method used is the Phase shift technology of constant velocity of wave.

Can know that from the principle of DAS time domain SAFT algorithm need be confirmed time-delay to relative distance between the reverberation and velocity of wave according to transducer.And acoustic characteristics such as for multi-layer body, sound wave can reflect at the interphase place of layer and layer, reflection, thereby travel path can change.In addition, the medium of each layer of object is different, and sound wave is also unequal usually in the velocity of propagation of each layer; Therefore; The calculating of time delay can not be adopted the simple computation method of relative distance divided by velocity of wave again, but needs to obtain earlier the travel path of wave beam, obtains the travel-time then piecemeal.Its emphasis and difficult point are to find out quickly and accurately the travel path of wave beam, and the ray trace technology realizes the best approach of this purpose 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 through iterative computation.To be this technology do not have special requirement to the interphase shape of multi-layer body to the advantage of ray trace technology, can be applicable to the ultrasonic imaging of arbitrary shaped body.But its important disadvantages is a ray tracing algorithm comprises interative computation, and time complexity is high, and simultaneously, the DAS principle is actually a kind of convolution algorithm, and its calculated amount is also bigger, and imaging time is long.Shown in Fig. 2 (b), this scheme generates this image (5cm*5cm) just needs 2 hours.This shortcoming seriously limited should technology promote the use.

For the imaging of multi-layer body, the Phase shift ultrasonic imaging method does not have demand raypath and convolution algorithm, and image taking speed can significantly improve, and only needs 4 minutes as generating the image shown in Fig. 2 (c).And the Phase shift technology is a kind of image reconstructing method based on frequency domain, than the SAFT technology under the time domain have higher distance to the orientation of resolution and Geng Gao to resolution and lower secondary lobe, its image quality is higher.But; Because the velocity of wave of hypothesis sound wave is constant in the Phase shift technology, make this method to be applicable to have foreign medium on the depth direction and horizontal direction must be medium of the same race that promptly layer is the regular situation on level straight line/plane perhaps parallel to each other and so on interphase between the layer medium; And for the object that all has foreign medium on depth direction and the horizontal direction; The object that promptly comprises non-regular interphase (curve/curved surface) since sound wave the velocity of wave in the medium is different throughout, and cause testing result accurate inadequately; Shown in Fig. 2 (c), grave error has appearred in the interfacial imaging results of the second layer.

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

Summary of the invention

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

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

The invention is characterized in; Adopt the Phase shift ultrasonic imaging technique; And the Phase shift technology improved, on the depth direction of testee and horizontal direction, considered to exist the situation of foreign medium respectively, in the 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 non-regular layering object.One-piece construction is as shown in Figure 3; In imaging process to the separatrix individual processing; Outside the depth range at place, separatrix; Use the Phase shift formula and the Fast Fourier Transform (FFT) technology of constant velocity of wave, and in this scope, use Phase shift technology and N the single output Fast Fourier Transform (FFT) technology that becomes velocity of wave, to reduce calculated amount.

The invention is characterized in, contain following steps successively:

Step (1): make up one forms by a computing machine, ultrasonic transducer, a cover register control and analog to digital converter one being used for based on change velocity of wave Phase shift on the vertical section that the degree of depth and horizontal both direction form, make the system of not damaged ultrasonic imaging at horizontal and vertical three layers of heterogeneous object that all have a foreign medium, wherein:

Said ultrasonic transducer is provided with: the pulse signal input terminal that links to each other with the output terminal of said register control; The input end of said register control links to each other with the corresponding positioning control signal output part of said computing machine; Said ultrasonic transducer also is provided with: the echoed signal output terminal that links to each other with the input end of said analog to digital converter, the output terminal of said analog to digital converter links to each other with the echo samples signal input part of said computing machine.Said transducer is controlled by said register control, moves with 1 step-length/ms fixed rate at three layers of heterogeneous body surface, and said register control is the gearing of the said ultrasonic transducer of control shift position, and its parameter is imported by said computing machine,

Three layers of heterogeneous object are L along the axial horizontal length of x, and it is interval to be divided into x of L/ Δ, and Δ x is a burst length, and the point 0 that also is said ultrasonic transducer from the x axle is to terminal point N _x-1 ends each step-length that moves, N _x=1+L/ Δ x, said ultrasonic transducer move the point that is reached at every turn and are called sensing point, total N _xIndividual, sequence number n _x=0,1 ..., N _x-1; Said register control produces a TTL pulse at each sensing point place; Trigger said ultrasonic transducer driving pulse of depth direction z emission perpendicular to the x axle to three layers of heterogeneous object; Transducer transfers receiving mode to and picks up counting subsequently, receive from the echoed signal of testee reflection, said analog to digital converter to said transducer at sensing point n _xThe echoed signal that the place receives is carried out N _tInferior sampling is also stored in the computing machine sampling sequence number n into _t=0,1 ..., N _t-1, SF is f _s, its value is preset for analog to digital converter, and (z, x are to be x=n at horizontal ordinate t) to note s _xThe echoed signal that Δ x place receives is at t=n _t/ f _sSampled value constantly;

Step (2): said computing machine is from n _tSensing point n is read in=0 beginning in regular turn _xThe sampled value at=0 place, and write down the sampling sequence number that sampled value fluctuates for the first time

Subscript 0 expression sensing point sequence number then, repeats above-mentioned steps and reads n successively _x=1 ..., N _xThe sampled value at-1 each sensing point place, and write down the sequence number that each point place sampled value fluctuates for the first time

Obtain N _xN altogether on the individual sensing point _xThe sequence number that individual sampled value fluctuates for the first time

From

In choose minimum sequence number

And calculate this sequence number pairing depth coordinate value on the z of said three layers of heterogeneous object direction

With from

In choose maximum sequence number

And calculate this sequence number pairing depth coordinate value on the z of said three layers of heterogeneous object direction

V wherein ₁Be the velocity of wave in the ground floor medium; Make two horizontal lines that are parallel to the x axle through these two depth coordinate values respectively; Be surrounded the separatrix between ground floor medium and second layer medium, form one and be positioned at first bounding box of crossing over first, second two layer medium on said three layers of heterogeneous object x-z vertical section simultaneously;

Step (3): (z, x t) carry out one dimension Fast Fourier Transform (FFT) 1D-FFT to the sampled value s at each sensing point place one by one along time shaft t _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 n for the sensing point sequence number _xThe level at place is to coordinate figure, and ω is the sampling angular frequency, then, and to above-mentioned x-ω territory 2D signal S (z ₀, x ω) carries out one dimension Fast Fourier Transform (FFT) 1D-FFT along transverse axis x successively _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, k be the x of wave-number vector to component, sequence number is n _k, said wave-number vector is represented the spatial gradient of wave phase,

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

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

Wherein, n _ωBe the sequence number of sampling angular frequency, value is at 0～(N _t-1) between, gets 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=0 each pixel of row;

Step (4): generate z direction homogeneity district z successively according to the following steps ₀＜z≤z _MinInterior each row image that forms by each row pixel respectively:

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

$S(z_{\min },k,\mathrm{\ω})=S(z_0,k,\mathrm{\ω})\underset{z=1}{\overset{z_{\min }}{\mathrm{\Π}}}\mathrm{\α}(\mathrm{\Δz},k,\mathrm{\ω},v_z)$

Wherein,

Wherein, i is an 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 _MinThe time, v _z=v ₁, v ₁Be the velocity of wave of sound wave in the ground floor medium, Δ z is the spacing distance between last adjacent two pixels of depth direction z, α (Δ z, k, ω, v _z) be the Phase shift amount under the constant velocity of wave,

Step (4.2): carry out following steps successively:

Step (4.2.1): get t=t ₀=0, the S (z that said step (4.1) is obtained _Min, k ω) makes the Gauss integration of independent variable ω,

Step (4.2.2): the Gauss integration result to said step (4.2.1) gets conjugation, again k is done the one dimension Fast Fourier Transform (FFT), multiply by 1/N after result of calculation is got conjugation again _x, obtain in the ground floor medium that coordinate figure is z on the depth direction _MinOne-row pixels put formed image:

$s(z_{\min },x,t_0)=\frac{1}{N_x}\{1D-{\mathrm{FFT}}_k\{\left[\frac{1}{2\mathrm{\π}}{\∫}_{-\mathrm{\π}}^{\mathrm{\π}}S(z_{\min },k,\mathrm{\ω})\mathrm{d\ω}{]}^{*}\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 the step (4) _z=v ₁, be that following step (5.1), step (5.2) are carried out in the step-length circulation to the z value with Δ z, up to z>=z _MaxEnd, generate depth direction z in first bounding box _Min＜z≤z _MaxBy formed each the row image of each row pixel, its step is following respectively in the interval:

Step (5.1): utilize the 2-d spectrum at the Phase shift formula calculating z+ Δ z place under the said 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 following steps successively:

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

Step (5.2.2): the Gauss integration result to said step (5.2.1) gets conjugation, again k is done the one dimension Fast Fourier Transform (FFT), multiply by 1/N after result of calculation is got conjugation again _x, the one-row pixels that obtains coordinate figure on the depth direction and be z+ Δ z is put formed image, z _Min＜z+ Δ z≤z _Max:

$s(z+\mathrm{\Δz},x,t_0)=\frac{1}{N_x}\{1D-{\mathrm{FFT}}_k\{\left[\frac{1}{2\mathrm{\π}}{\∫}_{-\mathrm{\π}}^{\mathrm{\π}}S(z+\mathrm{\Δz},k,\mathrm{\ω})\mathrm{d\ω}{]}^{*}\right){\}}^{*};$

Step (6): revise depth direction z in said first bounding box according to the following steps _MinTo z _MaxInterval z _Min＜z≤z _MaxIn each row image, the error that causes to eliminate the not homogeneity between the second layer and ground floor medium:

Step (6.1): the z that obtains with step (5) _MinTo z _MaxInterval image block uses Canny operator edge extracting algorithm to extract said z as input quantity _Min, z _MaxBetween ground floor medium and the separatrix c of second layer medium ₁(x, z),

Step (6.2): from z=z _MinBeginning is up to z>=z _MaxOnly, the z that obtains with said step (4.1) _Min2-d spectrum S (the z at place _Min, k ω) makes initial value, is step-length circulation execution in step (6.2.1) and step (6.2.2) with Δ z:

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 said 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]{\mathrm{<figure-callout\; id="2"\; label="e\; i"\; filenames="HDA0000138237050000012.png,HDA0000138237050000022.png"\; state="\{\{state\}\}">e</figure-callout>}}^i$

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 coordinate (x, the velocity of wave of sound wave in the medium of z) locating:

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

Calculation procedure is following:

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

$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) the 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 the 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, (z+ Δ z, x ω) make the Gauss integration of independent variable ω to the S that said step (6.2.1) is obtained, and the one-row pixels that obtains coordinate figure on the depth direction and be z+ Δ z is put formed image, z _Min＜z+ Δ z＜z _Max:

$s(z+\mathrm{\&Delta;z},x,t_0)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\mathrm{\&pi;}}^{\mathrm{\&pi;}}S(z+\mathrm{\&Delta;z},x,\mathrm{\&omega;})\mathrm{d\&omega;}$

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

Step (7): (5) described method set by step, get v _z=v ₂, be that said step (5.1), step (5.2) are carried out in the step-length circulation to the z value with Δ z, up to z>=z _DepthEnd, generate depth direction z _Max＜z≤z _DepthFormed each the row image of interval each row pixel, v ₂Be the velocity of propagation of sound wave in the second layer medium, z _DepthRepresent the maximum coordinates value of three layers of heterogeneous object on depth direction z, be the systemic presupposition value;

Step (8): the z that obtains with step (7) _MaxTo z _DepthInterval image block uses Canny operator edge extracting algorithm to extract said z as input quantity _Max, z _DepthBetween second layer medium and the separatrix c of the 3rd layer of medium ₂(x z), asks separatrix c ₂(x, z) the min coordinates value z ' on depth direction z _MinWith maximum coordinates value z ' _Max, respectively through said z ' _MinAnd z ' _MaxTwo depth coordinate values are made two horizontal lines that are parallel to the x axle; Be surrounded the separatrix between second layer medium and the 3rd layer of medium, form one and be positioned at second bounding box crossing over second, third two layer medium on said three layers of heterogeneous object x-z vertical section simultaneously;

Step (9): according to depth direction z ' in described method correction second bounding box of step (6.2) _MinTo z ' _MaxInterval z ' _Min＜z≤z ' _MaxIn each row image, the error that causes to eliminate the not homogeneity between the 3rd layer and second layer medium;

Step (10): (5) described method set by step, get v _z=v ₃, to the z value with z ' _MaxBeing initial value, is that said step (5.1), step (5.2) are carried out in the step-length circulation with Δ z, up to z>=z _DepthEnd, generate in the 3rd layer of medium except said second bounding box with formed each the row image of each row pixel of exterior portions, v ₃It is the velocity of propagation of sound wave 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, the transducer moving step length is 0.5mm, and transducer ultrasonic waves transmitted centre frequency is 5MHz; SF 100MHz; The step delta z of depth direction is taken as 0.05mm, and then ground floor interphase error can be controlled in the 0.4mm, and second layer interphase can be controlled at 1.0mm with interior (like Fig. 4).

The present invention has following advantage compared with prior art:

1. can carry out fast ultrasonic imaging to containing the interfacial object of non-rule;

2. the resolution and the precision of imaging are higher.

Description of drawings

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

The comparison diagram of Fig. 2 right and wrong rule layering object ultrasonic imaging method under the identical experiment environment: 2 (a) are the sectional view of testee; 2 (b) are that time domain SAFT and ray trace technology combines image that 2 (a) are generated, and imaging time is more than 2 hours, and to the unusual sensitivity of the parameter of algorithm, robustness is not fine; 2 (c) are the image that the Phase shift ultrasonic imaging technique of constant velocity of wave generates, imaging time 4 minutes, and its error is bigger; The image that 2 (d) generate for this formation method, imaging time 6 minutes.

Fig. 3 is this formation method overall process synoptic 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; The ground floor separatrix that the imaginary point curve obtains for this formation method, dashed curve are the marginal graph of errors of ground floor that this formation method obtained; 4 (b) are second layer separatrix Error Graph, and solid-line curve is the marginal standard value of the testee second layer, the second layer separatrix that the imaginary point curve obtains for this formation method, and dashed curve is the marginal graph of errors of the second layer that this formation method obtained.

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

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

Fig. 7 is a ultrasonic transducer work synoptic diagram.

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

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

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

Embodiment

Practical implementation process of the present invention comprises three parts (like Fig. 5): ultrasound data obtains, forms images and calculates and the image demonstration.The hardware platform system structural drawing is as shown in Figure 6; Ultrasonic image-forming system is made up of a computing machine, ultrasonic transducer, a cover register control and an analog to digital converter; The pulse signal input terminal of ultrasonic transducer links to each other with the output terminal of register control, and the input end of register control links to each other with the positioning control signal output part of computing machine.The echoed signal output terminal of ultrasonic transducer links to each other with the input end of analog to digital converter, and the output terminal of analog to digital converter links to each other with the echo samples signal input part of computing machine.

The mode that ultrasound data obtains can be selected: contact, the i.e. direct contact measured object surfaces of ultrasonic transducer; Perhaps liquid immersion type, promptly object under test places 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 uses single emission/receiving transducer, and transducer moves with about 1 step/ms fixed rate along x direction (like Fig. 7) with uniform step delta x in testee surface or liquid through register control.A moment generation TTL (transistor-transistor logic level) pulse that controller is stable in each target location; This pulse is used for triggering ultrasonic transducer to testee and the perpendicular driving pulse of depth direction emission of x direction; Transducer transfers receiving mode to and picks up counting subsequently, receives from the echo of testee reflection.Each position of transducer transponder pulse and reception echo is a sensing point.The echoed signal that transducer receives is by the analog to digital converter collection and be stored in the storer.Δ x needs to require to confirm that its value is more little that the image of generation is accurate more based on the actual size of object under test and imaging precision, but computing time is also long more.

Imaging is calculated and is imported as computing machine with the sampled data at testee each sensing point place on a vertical section exactly; Calculate the skiagraph picture of testee then by following image-forming step (process flow diagram is referring to Fig. 8); S (z wherein; X t) is illustrated in the signals sampling value that the transducer at horizontal direction x place receives at t constantly; Z representes the coordinate figure of depth direction; Δ z representes the step-length (relevant with precision) of depth direction, and the image that is promptly generated is in the spacing of vertical (that is depth direction) last adjacent two pixels; z _DepthBe preset image maximum depth value in the vertical:

Step (1): the ground floor medium is the isomorphism medium, and its ultrasonic reflection signal is identical, and at ground floor and second layer boundary, ultrasonic reflection signal can produce fluctuation.Then, search for the moment of each sensing point place generation signals fluctuation for the first time, and from these fluctuations constantly, find out the minimum t that constantly is worth according to the sampled data of ultrasound echo signal _MinBe worth t constantly with maximum _Max, then according to the velocity of wave v in the 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(like Fig. 9);

Step (2): get depth direction initial value z=z ₀=0, parameter x and t are carried out 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, ω is that the x of wave-number vector (wavenumber vector) is to component (wave-number vector refers to the spatial gradient of wave phase) for sampling angular frequency, k;

Step (3): generate depth direction z ₀To z _MinInterval (z ₀≤z≤z _Min) image.Because the medium of testee in this section interval is identical, so this section image is also just the same:

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

$S(z_{\min },k,\mathrm{\&omega;})=S({\mathrm{<figure-callout\; id="0"\; label="z"\; filenames="HDA0000138237050000022.png"\; state="\{\{state\}\}">z</figure-callout>}}_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 an imaginary unit, v _zFor testee acoustic velocity in the capable medium of z on depth direction, as z≤z _MinThe time, v _z=v ₁, v ₁Velocity of wave for sound wave in the ground floor medium;

Step (3.2): make two-dimensional inverse Fourier transform, make t=t ₀=0, obtaining the along slope coordinate value is delegation's 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;},({\mathrm{<figure-callout\; id="0"\; label="z"\; filenames="HDA0000138237050000022.png"\; state="\{\{state\}\}">z</figure-callout>}}_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 the z value with Δ z step-length circulation execution in step (4.1) and step (4.2) up to z>=z _Max:

Step (4.1): utilize the 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 the along slope coordinate value is delegation's 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 velocity of wave of ground floor has only been used in the calculating in the step (4), soon should still be regarded as medium of the same race in the interval, the result of calculation that obtains certainly exists error, so this step needs to revise.Specifically comprised for two steps:

Step (5.1): the z that step (4) is obtained _MinTo z _MaxInterval image uses Canny edge extracting operator extraction z _MinWith z _MaxBetween separatrix c (x, z), z _Min≤z≤z _MaxThe computation process of Canny edge extracting operator was divided into for four steps:

Step (a): image smoothing.To original image with the two-dimensional Gaussian function smothing filtering obtain image I (x, z);

Step (b): computed image I (x, z) the gradient M of each pixel and direction Q.Adopt the first approximation of 2 * 2 templates conduct, promptly to the partial differential of 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)$

Then gradient magnitude M and direction Q do

$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 (x, z) less than along its gradient direction Q (x, the gradient magnitude of two consecutive point z) explain that this point is not a marginal point, then with I (x, gray scale z) is made as 0;

Step (d): the dual threashold value is handled.Configure the dual threshold method detects and adjoining edge needs low threshold value Low and high threshold High, gradient image is carried out the dual threashold value handle.Gradient magnitude is the edge greater than high threshold High's; Gradient magnitude is not the edge less than low threshold value Low's; Gradient magnitude is marginal, then judges in eight neighbors of this pixel whether have the edge pixel greater than high threshold High, then is edge pixel if exist, otherwise is not;

Step (5.2): to the z value from z=z _MinBeginning, with Δ z step-length circulation execution in step (5.2.1) and step (5.2.2) up to z>=z _Max:

Step (5.2.1): utilize the Phase shift formula that becomes velocity of wave to calculate the k-ω territory signal at z+ Δ z place:

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:

The value of vz (x) is following:

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

Step (5.2.2): make z ← z+ Δ z, t=t ₀=0, make two-dimensional inverse Fourier transform, obtaining the along slope coordinate value is delegation's 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 ₂, the method for (4) is calculated z set by step _MaxTo z _DepthEach interval row image;

Step (7): handle follow-up separatrix.Image section (the z that step (6) is generated _Max＜z≤z _Depth) (x z), makes v to use the new separatrix c of Canny edge extracting operator extraction ₁, v ₂Be respectively separatrix c (x, the z) velocity of propagation of sound wave in the two layer medium up and down at place, and make z _MinAnd z _Max(x, the z) up-and-down boundary on depth direction form new bounding box, if z to be respectively new separatrix c _Min==z _Max==z _Depth, then finish; Otherwise execution in step (5.2) is revised z _MinTo z _MaxInterval image;

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

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

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

Need make two-dimensional inverse Fourier transform and get t=t in step (3.2) and the step (4.2) ₀=0 value then can shilling t=t as view data ₀=0, and ω made integration, and then k is made 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 (like Gauss Gaussian integration) to try to achieve fast.And the Fourier inversion of k is realized through one dimension Fast Fourier Transform (FFT) (1D-FFT) still the data that just earlier numerical integration obtained are got conjugation ([Integration _ω(S (z, k, ω))] *), directly utilize the 1D-FFT subroutine to carry out computing then, again operation result is got conjugation one time, and multiply by 1/N _x, can obtain.N _xFor transducer moves step number, i.e. N _x=1+L _x/ Δ x, wherein L _xBe the length of object under test in the x direction.Therefore, the concrete implementation procedure of formula (3) (5) is expressed as:

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

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

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

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

Wherein, Δ x is the step-length that transducer moves along the x direction, f _sBe SF (its value is systemic presupposition), N _x=1+L _x/ Δ x is that transducer is counted L along the detection of x direction _x/ Δ x is that transducer moves step number, L _xBe the length (x direction) of testee, N _tFor transducer at each sensing point place the sampling number to echoed signal.Because the Phase shift ultrasonic imaging technique is regarded as the explosive reflection model with its working model, only consider that sound wave is from reverberation (along z axle negative direction) part of propagating upwards, so at Phase shift amount α (Δ z, k, ω, v _z) calculating in only get the value of ω＜0 part, be 0 and make the Phase shift amount of ω>=0 part, simultaneously if b＜0 then get

Know by the periodicity of ω value, when requiring ω＜0, desirable N in the calculating _t/ 2≤n _ω＜N _tPart.

In the 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 merge in the Fourier inversion of step (5.2.2), promptly

$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 the formula (8) can utilize numerical integration method (like Gauss Gaussian integration) to try to achieve, and internal layer integral operation (suc as formula (9)) then needs to try to achieve (suc as formula (10)) through revising 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]{\mathrm{<figure-callout\; id="2"\; label="e\; i"\; filenames="HDA0000138237050000012.png,HDA0000138237050000022.png"\; state="\{\{state\}\}">e</figure-callout>}}^i---\left(10\right)$

In the formula (10), n _xBe the sequence number of sensing point on the x direction at transducer place, span is same as n _k, and α (Δ z, k, ω, v ( _zX)) 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, use single output Fast Fourier Transform Inverse (FFTI) method (NsIFFTPruning) among the present invention N time.This method is based on the thought of butterfly computation in the Fast Fourier Transform (FFT), but a value in the transform domain is only exported in conversion each time.The NsIFFTPruning specific algorithm is following:

For?n _x＝0?to?N _x-1

For?n _k＝0?to?N _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 _xIn the time of fixedly, 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, promptly

$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}---\left(11\right)$

Formula (11) can use IDFT or IFFT algorithm to calculate, but different with general IDFT problem is: formula (11) only need be calculated G (n _x, n _k) at n _xThe value at this time domain point place of=h if directly use IDFT or IFFT algorithm, then certainly exists a large amount of redundant computation, for example, and in Figure 10, if only need to calculate n _x=2 point, then the butterfly computation of dotted portion is redundant among the figure.In order to reduce calculated amount, the present invention uses the fft algorithm (IFFTPruning) of cutting output.This algorithm adopts the butterfly computation method still based on fft algorithm, and its principle is following: according to butterfly group technology in the fft algorithm, calculating is divided into H step (Stage), wherein H satisfies 2 ^H=N _xStep h (h=1 ..., H) in, with 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+1Mod 2 ^H-1Point outside, remaining the point all be the redundancy.In each calculation procedure, only calculating is respectively organized and is numbered b in the butterfly computation _hThe point.Shown in figure 10, N _x=8, n _x=2 o'clock, calculate in every group of butterfly computation in the step (1) and be numbered b ₁Calculate in every group of butterfly computation in=0 the point, step (2) and be numbered b ₂Calculate in every group of butterfly computation in=2 the point, step (3) and be numbered b ₃=2 point.

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

Image shows that the two-dimensional image data that the calculation stages that is about to form images obtains is presented on the display, as required display gray scale image or coloured image.

To liking the situation of regular layering object, (x z) is straight line to the separatrix c of use Canny edge extracting operator extraction, and the imaging calculation procedure is constant in the image-forming step (5.1) for to be measured.

For the situation that only contains single medium in the object to be measured (be illustrated as contact ultrasonic measurement mode, otherwise should be regarded as regular layering object), then z _Min==z _Max==z _Depth, imaging promptly stops after calculating execution of step (3).

Claims (1)

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

Step (1): make up one forms by a computing machine, ultrasonic transducer, a cover register control and analog to digital converter one being used for based on change velocity of wave Phase shift on the vertical section that the degree of depth and horizontal both direction form, make the system of not damaged ultrasonic imaging at horizontal and vertical three layers of heterogeneous object that all have a foreign medium, wherein:

Said ultrasonic transducer is provided with: the pulse signal input terminal that links to each other with the output terminal of said register control; The input end of said register control links to each other with the corresponding positioning control signal output part of said computing machine; Said ultrasonic transducer also is provided with: the echoed signal output terminal that links to each other with the input end of said analog to digital converter, the output terminal of said analog to digital converter links to each other with the echo samples signal input part of said computing machine.Said transducer is controlled by said register control, moves with 1 step-length/ms fixed rate at three layers of heterogeneous body surface, and said register control is the gearing of the said ultrasonic transducer of control shift position, and its parameter is imported by said computing machine,

Three layers of heterogeneous object are L along the axial horizontal length of x, and it is interval to be divided into x of L/ Δ, and Δ x is a burst length, and the point 0 that also is said ultrasonic transducer from the x axle is to terminal point N _x-1 ends each step-length that moves, N _x=1+L/ Δ x, said ultrasonic transducer move the point that is reached at every turn and are called sensing point, total N _xIndividual, sequence number n _x=0,1 ..., N _x-1; Said register control produces a TTL pulse at each sensing point place; Trigger said ultrasonic transducer driving pulse of depth direction z emission perpendicular to the x axle to three layers of heterogeneous object; Transducer transfers receiving mode to and picks up counting subsequently, receive from the echoed signal of testee reflection, said analog to digital converter to said transducer at sensing point n _xThe echoed signal that the place receives is carried out N _tInferior sampling is also stored in the computing machine sampling sequence number n into _t=0,1 ..., N _t-1, SF is f _s, its value is preset for analog to digital converter, and (z, x are to be x=n at horizontal ordinate t) to note s _xThe echoed signal that Δ x place receives is at t=n _t/ f _sSampled value constantly;

Step (2): said computing machine is from n _tSensing point n is read in=0 beginning in regular turn _xThe sampled value at=0 place, and write down the sampling sequence number that sampled value fluctuates for the first time

Subscript 0 expression sensing point sequence number then, repeats above-mentioned steps and reads n successively _x=1 ..., N _xThe sampled value at-1 each sensing point place, and write down the sequence number that each point place sampled value fluctuates for the first time

Obtain N _xN altogether on the individual sensing point _xThe sequence number that individual sampled value fluctuates for the first time

From

In choose minimum sequence number

And calculate this sequence number pairing depth coordinate value on the z of said three layers of heterogeneous object direction

With from

In choose maximum sequence number

And calculate this sequence number pairing depth coordinate value on the z of said three layers of heterogeneous object direction

V wherein ₁Be the velocity of wave in the ground floor medium; Make two horizontal lines that are parallel to the x axle through these two depth coordinate values respectively; Be surrounded the separatrix between ground floor medium and second layer medium, form one and be positioned at first bounding box of crossing over first, second two layer medium on said three layers of heterogeneous object x-z vertical section simultaneously;

Step (3): (z, x t) carry out one dimension Fast Fourier Transform (FFT) 1D-FFT to the sampled value s at each sensing point place one by one along time shaft t _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 n for the sensing point sequence number _xThe level at place is to coordinate figure, and ω is the sampling angular frequency, then, and to above-mentioned x-ω territory 2D signal S (z ₀, x ω) carries out one dimension Fast Fourier Transform (FFT) 1D-FFT along transverse axis x successively _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, k be the x of wave-number vector to component, sequence number is n _k, said wave-number vector is represented the spatial gradient of wave phase,

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

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

Wherein, n _ωBe the sequence number of sampling angular frequency, value is at 0～(N _t-1) between, gets 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=0 each pixel of row;

Step (4): generate z direction homogeneity district z successively according to the following steps ₀＜z≤z _MinInterior each row image that forms by each row pixel respectively:

Step (4.1): utilize following Phase shift formula to calculate said 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 an 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 _MinThe time, v _z=v ₁, v ₁Be the velocity of wave of sound wave in the ground floor medium, Δ z is the spacing distance between last adjacent two pixels of depth direction z, α (Δ z, k, ω, v _z) be the Phase shift amount under the constant velocity of wave,

Step (4.2): carry out following steps successively:

Step (4.2.1): get t=t ₀=0, the S (z that said step (4.1) is obtained _Min, k ω) makes the Gauss integration of independent variable ω,

Step (4.2.2): the Gauss integration result to said step (4.2.1) gets conjugation, again k is done the one dimension Fast Fourier Transform (FFT), multiply by 1/N after result of calculation is got conjugation again _x, obtain in the ground floor medium that coordinate figure is z on the depth direction _MinOne-row pixels put formed image:

$s(z_{\min },x,t_0)=\frac{1}{N_x}\{1D-{\mathrm{FFT}}_k\{\left[\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\mathrm{\&pi;}}^{\mathrm{\&pi;}}S(z_{\min },k,\mathrm{\&omega;})\mathrm{d\&omega;}{]}^{*}\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 the step (4) _z=v ₁, be that following step (5.1), step (5.2) are carried out in the step-length circulation to the z value with Δ z, up to z>=z _MaxEnd, generate depth direction z in first bounding box _Min＜z≤z _MaxBy formed each the row image of each row pixel, its step is following respectively in the interval:

Step (5.1): utilize the 2-d spectrum at the Phase shift formula calculating z+ Δ z place under the said 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 following steps successively:

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

Step (5.2.2): the Gauss integration result to said step (5.2.1) gets conjugation, again k is done the one dimension Fast Fourier Transform (FFT), multiply by 1/N after result of calculation is got conjugation again _x, the one-row pixels that obtains coordinate figure on the depth direction and be z+ Δ z is put formed image, z _Min＜z+ Δ z≤z _Max:

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

Step (6): revise depth direction z in said first bounding box according to the following steps _MinTo z _MaxInterval z _Min＜z≤z _MaxIn each row image, the error that causes to eliminate the not homogeneity between the second layer and ground floor medium:

Step (6.1): the z that obtains with step (5) _MinTo z _MaxInterval image block uses Canny operator edge extracting algorithm to extract said z as input quantity _Min, z _MaxBetween ground floor medium and the separatrix c of second layer medium ₁(x, z),

Step (6.2): from z=z _MinBeginning is up to z>=z _MaxOnly, the z that obtains with said step (4.1) _Min2-d spectrum S (the z at place _Min, k ω) makes initial value, is step-length circulation execution in step (6.2.1) and step (6.2.2) with Δ z:

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 said 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 coordinate (x, the velocity of wave of sound wave in the medium of z) locating:

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

Calculation procedure is following:

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

$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) the 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 the 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, (z+ Δ z, x ω) make the Gauss integration of independent variable ω to the S that said step (6.2.1) is obtained, and the one-row pixels that obtains coordinate figure on the depth direction and be z+ Δ z is put formed image, z _Min＜z+ Δ z＜z _Max:

$s(z+\mathrm{\&Delta;z},x,t_0)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\mathrm{\&pi;}}^{\mathrm{\&pi;}}S(z+\mathrm{\&Delta;z},x,\mathrm{\&omega;})\mathrm{d\&omega;}$

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

Step (7): (5) described method set by step, get v _z=v ₂, be that said step (5.1), step (5.2) are carried out in the step-length circulation to the z value with Δ z, up to z>=z _DepthEnd, generate depth direction z _Max＜z≤z _DepthFormed each the row image of interval each row pixel, v ₂Be the velocity of propagation of sound wave in the second layer medium, z _DepthRepresent the maximum coordinates value of three layers of heterogeneous object on depth direction z, be the systemic presupposition value;

Step (8): the z that obtains with step (7) _MaxTo z _DepthInterval image block uses Canny operator edge extracting algorithm to extract said z as input quantity _Max, z _DepthBetween second layer medium and the separatrix c of the 3rd layer of medium ₂(x z), asks separatrix c ₂(x, z) the min coordinates value z ' on depth direction z _MinWith maximum coordinates value z ' _Max, respectively through said z ' _MinAnd z ' _MaxTwo depth coordinate values are made two horizontal lines that are parallel to the x axle; Be surrounded the separatrix between second layer medium and the 3rd layer of medium, form one and be positioned at second bounding box crossing over second, third two layer medium on said three layers of heterogeneous object x-z vertical section simultaneously;

Step (9): according to depth direction z ' in described method correction second bounding box of step (6.2) _MinTo z ' _MaxInterval z ' _Min＜z≤z ' _MaxIn each row image, the error that causes to eliminate the not homogeneity between the 3rd layer and second layer medium;

Step (10): (5) described method set by step, get v _z=v ₃, to the z value with z ' _MaxBeing initial value, is that said step (5.1), step (5.2) are carried out in the step-length circulation with Δ z, up to z>=z _DepthEnd, generate in the 3rd layer of medium except said second bounding box with formed each the row image of each row pixel of exterior portions, v ₃It is the velocity of propagation of sound wave in the 3rd layer of medium.

Images