CN103018333A - Synthetic aperture focused ultrasonic imaging method of layered object - Google Patents

Abstract

The invention relates to a synthetic aperture focused ultrasonic imaging method of a layered object, which belongs to the technical field of fast imaging of an inner structure of a multilayered heterogenic object. The method is characterized by comprising the steps: based on the condition that a ultrasonic transducer moves at equal distance to probe in the horizontal direction of the multilayered heterogenic object, changing the implementation process of original SAFT to a process of drawing a track curve on an image by sampling data based on synthetic aperture focusing technology (SAFT) and refraction law; combining the refraction law to calculate the refraction vector positively; obtaining pixel points on a refraction path by a linear scanning and converting technology, so as to improve the imaging speed. The method is not only suitable for ultrasonic imaging of a regular layered object with heterogenic media in the depth direction, but also suitable for ultrasonic imaging of irregular and regular layered objects with heterogenic media in the depth and horizontal directions. The imaging speed is fast, and the imaging precision is high.

Description

The ultrasonic imaging by synthetic aperture focusing method of layering object

Technical field

The present invention relates to Ultrasonic Nondestructive technology, synthetic aperture focusing imaging technique, refraction law and line scan conversion technology, realize surface configuration and inner structure fast imaging to comprising regular interphase or irregular interfacial layering object.

Background technology

At present, in the Ultrasonic Nondestructive field, (the Ultrasonic Detection mode comprises contact and two kinds of patterns of liquid immersion type for the ultrasonic imaging that contains interfacial layering object, if the couplant layer is considered as the part of testee, then the liquid immersion type Ultrasonic Detection also can be considered the contact Ultrasonic Detection to the layering object, therefore among the present invention both unifications are called Ultrasonic Detection and the imaging of layering object), mainly contain two schemes: the one, the method that adopts synthetic aperture focusing (SAFT) ultrasonic imaging technique 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.

Synthetic aperture focusing technology (SAFT) stems from Synthetic Aperture Radar Technique (SAR), be introduced in the ultrasonic imaging field in early 1970s, its have be not subjected to that near field region restriction, azimuthal resolution are high, resolution only relevant with the ultrasonic transducer size and with the advantage such as range-independence.The ultimate principle of SAFT ultrasonic imaging technique is to utilize pulse-echo (pulse-echo) measurement mechanism, use a ultrasonic transducer along fixation locus testee to be carried out orderly scanning, and adopt time-delay stack (DAS) method that the pulse echo signal that scanning obtains is carried out focal imaging, reach the purpose that the ultrasonic transducer that utilizes single smaller aperture due is simulated large array of apertures.SAFT ultrasonic imaging working model is shown in Fig. 1 (a), ultrasonic transducer along the direction of scanning (X-direction) do equally spaced movement at body surface, at depth direction (Z-direction) the emission ultrasonic signal of each scanning position to object, interior of articles is surveyed, ultrasonic transducer receives echoed signal and the Sampling hold that object internal reflection thing is reflected back simultaneously, and the sampled data that at last all scanning position places is obtained carries out DAS and the multiple spot dynamic focusing is processed and the demonstration image.

DAS need to calculate different time lag curve to the focus point on the different depth of target imaging zone from the multiple spot dynamic focus technology.Shown in Fig. 1 (a), in order to focus at target reverberation place (x, z), the echoed signal that the SAFT technology obtains ultrasonic transducer each scanning position place in its synthetic aperture effective length L overlap-add procedure of delaying time: establish s _i(t) be that ultrasonic transducer is at u _iThe echoed signal that the place receives, t is sampling instant, u _iThe place about the time-delay of target reverberation (x, z) is

$t_i=\frac{2r_i}{v}=\frac{2}{\mathrm{<figure-callout\; id="2"\; label="v\; z"\; filenames="HDA00002540433900031.png"\; state="\{\{state\}\}">v</figure-callout>}}\sqrt{z^{}}\sqrt{{}^2+{(x-u_i)}^2},i=\mathrm{0,1},...,L-1.---\left(1\right)$

Wherein, v is ultrasonic velocity of propagation in medium, r _iFor (x, z) puts apart from u _iAir line distance.Time-delays all in the synthetic aperture effective length L consist of a duration curve, and this curve is one section hyperbolic curve.The computing formula of L is

L=0.84λz/D (2)

λ is ultrasonic wavelength in medium, and D is the diameter of ultrasonic transducer, then being imaged as that (x, z) locates

$I(x,z)=\underset{i=0}{\overset{L-1}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&omega;}}_i{s}_i\left(t_i\right)}{r_i}---\left(3\right)$

Wherein, ω _iBe apodizing function.

From the principle of DAS as can be known, the SAFT algorithm need to be determined time-delay to the relative distance between the reverberation and velocity of wave according to ultrasonic transducer.And when the testee layering, shown in Fig. 1 (b), the ultrasonic signal of ultrasonic transducer emission can reflect at the interphase place, and the travel path of sound wave can change, target reverberation (x, z) and u _iApart from r _iNo longer be the length of point-to-point transmission straight-line segment.In addition, the medium of each layer of layering object is also different, and sound wave is usually also unequal in the velocity of propagation of each layer, therefore, the calculating of time delay can not adopt relative distance divided by the simple computation method of velocity of wave again among the DAS, and duration curve also no longer is the hyperbolic shape of rule.In order to try to achieve time-delay, need to obtain first sound wave at u _iWith the travel path between the target reverberation (x, z), then obtain piecemeal the length of each route segment and travel-time (such as r among Fig. 1 (b) _I1, r _I2And t _i).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 purpose just.So what Tracing Technology was very natural is incorporated in the SAFT technology, to realize the ultrasonic imaging to the layering object.

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.So, such as Fig. 1 (b), in the method for SAFT in conjunction with ray trace, when calculating u _iWith the Acoustic Wave Propagation time-delay t between the target reverberation (x, z) _iThe time, need to carry out iteration along all points on pair line of separatrix c (x, z)=0, calculate each and put corresponding two sections travel path r _I1, r _I2Length and travel-time thereof, find out so that the shortest point of travel-time, this point is refraction point.

To be Tracing Technology do not have special requirement to the interphase shape of layering object to the advantage of the method, can not only be applicable to have foreign medium and horizontal direction is medium of the same race at depth direction, namely layer and interphase between layer medium are the ultrasonic imaging of the regular layering object on level or plane parallel to each other and so on, also can be applicable on depth direction and horizontal direction, all have foreign medium, namely comprise the ultrasonic imaging of the irregular layering object of irregular interphase (such as curved surface).But its serious shortcoming is that ray tracing algorithm comprises interative computation, time complexity is high, in addition, comprise a root mean square computing in the formula (1) of calculation delay in the former SAFT technology, computing cost is very large, and after introducing Tracing Technology, the calculating of time delay then comprises repeatedly root mean square computing, in Fig. 1 (b), when calculating u _iWith the Acoustic Wave Propagation time-delay t between the target reverberation (x, z) _iThe time, need to each point on the separatrix, calculate respectively rms twice to obtain two sections path r _I1, r _I2Length, its computing cost is larger.So the method is very consuming time, shown in Fig. 3 (b), with the method testee (4.4cm*4.9cm) imaging shown in Fig. 3 (a) is just needed 30 minutes.In addition, because hyperacoustic wavelength is relevant with medium, behind the dielectric stratifying of testee, the computing formula (2) of synthetic aperture effective length L is also no longer applicable, yet the degree of accuracy of SAFT technology imaging and the accuracy of L are closely related, if L is excessive, the imaging noise is more, and signal to noise ratio (S/N ratio) is low; If L is too small, imaging can be lost the details of testee, and precision is not high.But for the irregular layering object of complexity, the formula that lacks at present accurate calculating L, in the method for SAFT in conjunction with ray trace, be difficult to find out accurately L value, cause the imaging effect of the method also bad, for example among Fig. 3 (b), the image that the method becomes the testee longitudinal profile fails to reconstruct preferably the curvature portion in the second layer interphase.

Phase shift (Phase Shift Migration) ultrasonic imaging method is that the migrating technology (Migration Technique) in the reflection seismology (Reflection Seismology) is incorporated into the ultrasonic imaging field, and the ultrasonic imaging method under a kind of frequency domain that obtains.The method is considered as the explosive reflection model with the supersonic sounding model, supposes that the reverberation in the object under test explodes constantly simultaneously at t=0, and the bursting strength of each reverberation is proportional to its reflectivity, and whole field intensity is measured with one group of ultrasonic transducer.Its main thought is according to extrapolating to calculate the sound field of other positions of depth direction from horizontal level (the being the depth direction the first row) sound field that observes.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 _i,t))

Second step is the circulation at depth direction, and the 2-d spectrum that first the last time circulation is obtained is done phase shift, then makes two-dimensional inverse Fourier transform and gets t0, obtains delegation's image,

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

Wherein,

Be the corresponding Phase shift amount of Δ z step-length on the depth direction, c is ultrasonic velocity of propagation in medium, and its value is constant.

For the imaging of multi-layer body, the Phase shift ultrasonic imaging method is without the demand raypath, and image taking speed can significantly improve, and only needs 62s as generating the image shown in Fig. 3 (c).But, because the velocity of wave of hypothesis sound wave is constant in the Phase shift technology, horizontal direction must be the ultrasonic imaging of the regular layering object of medium of the same race so that the method can only be applicable to have foreign medium on the depth direction, and for comprising irregular interfacial object, all there is foreign medium at depth direction and horizontal direction, and cause testing result accurate not, shown in Fig. 3 (c), grave error has appearred 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.

In conjunction with finding out the principle of Tracing Technology, the method part the most consuming time is to seek the interative computation step of optimum refraction point from SAFT.From Tracing Technology and DAS principle as can be known, the method is a kind of reverse computation process, need to determine first a pixel on the image, then find out all scanning position points corresponding to this point, again each scanning position point is found out optimum refraction point to obtain distance and the time-delay between picture point and this scanning position point in iteration on the separatrix, the data point on the duration curve is added up obtain the pixel value of this picture point at last.Because in advance and do not know the particular location of testee internal object reverberation, use dynamic focus technology that all pixels in the image are all carried out this reverse computation process one time with regard to needs, to generate entire image, thereby so that at pixel place corresponding to reverberation, ultrasonic signal obtains consistent stack, reach the superposed strength maximization, realize focusing on, and at other pixel places, the stack of ultrasonic signal is chaotic, and superposed strength is difficult to maximization.Therefore, in image, the cumulative pixel value that obtains in reverberation institute corresponding position will be obviously greater than other pixels.

If the computation process of all pixels on the image is done whole the investigation, can find that each data that each scanning position place sampling obtains not only act on pixel corresponding to reverberation, also act on the pixel of no reflection events thing.In Fig. 1 (c), u _iThe sampled data s at place _i(t _i) not only participated in the imaging calculating that (x, z) puts, also participated in segment of curve The calculating of upper other points.

U _iOne section curve in the sound field scope of place's ultrasonic transducer emission, the point on this curve is apart from u _iTravel-time t _iAll identical, but apart from r _iNot necessarily identical.In the ground floor medium, ultrasonic transducer is at u _iThe sound field at place is positioned at the half-power beam angle β of this transducer _0.5In (namely in figure from u _iBetween two imaginary point lines that set out in the place) because synthetic aperture effective length L also is defined as L=z * β simultaneously _0.5, in conjunction with formula (2), so β _0.5=0.84 λ/D.And in second layer medium, because medium changes, sound wave reflects at the interphase place, former half-power beam angle β _0.5Two boundary lines also can reflect (such as the imaginary point line among Fig. 1 (c)), thereby u _iThe sound field scope of place's ultrasonic transducer also can change, but sound field is still between the represented boundary line of two imaginary point lines.So, segment of curve Can be regarded as when refraction point from left to right moves in the sound field scope along the separatrix between ground floor and the second layer medium, ultrasonic transducer is at u _iThe sampled data s at place _i(t _i) track of on image, passing by.When second layer medium and ground floor identical, then without layering, this moment this track be one section with u _iFor the center of circle, angle are β _0.5Circular arc; And work as two layer medium not simultaneously, the shape of this track is then relevant with marginal shape.

Therefore, remove Integrated Understanding overall calculation process, sampled data s from forward _i(t _i) effect in the imaging of entire image is calculated just is equivalent to data value ω _is _i(t _i)/r _iIn image, drawn one section curve

It is the pixel value of each pixel on the curve

In former L scope, s _i(t _i) all data (s on the duration curve at place ₀(t ₀) to s _L-1(t _L-1)) corresponding segment of curve (

Extremely

) intersection point in image is the target reflection object point (x, z) in the former reverse computation process.

If in entire image, only consider u _iThere is sampled data s in the place _i(t) effect, then corresponding to a width of cloth with u _iCentered by irregular sector diagram I _Ui(such as Fig. 1 (d)), namely

$I_{u_i}=\underset{t_i=0}{\overset{N}{\mathrm{\&Sigma;}}}I\left({\mathrm{tra}}_{s_i\left(t_i\right)}\right)---\left(4\right)$

N is that ultrasonic transducer is at the sampling number at each scanning position place.Former SAFT then can be interpreted as the stack of each corresponding sector diagram in scanning position place again in conjunction with the DAS computing formula (3) in the ray tracing method, that is:

$I(x,z)=\underset{i=0}{\overset{M-1}{\mathrm{\&Sigma;}}}{I}_{u_i}(x,z)---\left(5\right)$

Wherein, M is the sum of ultrasonic transducer scanning position point.

So former SAFT can realize by draw geometric locus with each sampled data at each scanning position place on image in conjunction with the method for ray trace.In order accurately to calculate geometric locus

Need to try to achieve the function expression of this curve.Such as Fig. 1 (c), postulated point P (x _P, z _P) be geometric locus On any point, ultrasound wave is from ultrasonic transducer position U (u _i, 0) locate to launch the some R (x on the separatrix _R, z _R) to locate to reflect the rear arrival P used time of point be t _i/ 2, namely

$\left|\stackrel{\&RightArrow;}{\mathrm{RU}}\right|/v_{\mathrm{water}}+\left|\stackrel{\&RightArrow;}{\mathrm{RP}}\right|/v_{\mathrm{object}}=t_i/2---\left(6\right)$

Wherein, v _WaterBe the velocity of propagation of sound wave in water, v _ObjectBe the velocity of propagation of sound wave in testee.Simultaneously, according to refraction law:

Eta=sin α _i/ sin α _r=v _Water/ v _Object(7) wherein, α _iBe incident angle, α _rBe the refraction angle, eta is refractive index.Coordinate by formula (6) and formula (7) can be ordered in the hope of P namely gets geometric locus Function, but easier in order to calculate, can calculate the refraction vector with the formula in the computer graphics

$\stackrel{\&RightArrow;}{\mathrm{RP}}=\mathrm{eta}\&times;\stackrel{\&RightArrow;}{\mathrm{RU}}-[\mathrm{eta}\&times;(\stackrel{\&RightArrow;}{\mathrm{RU}}\&CenterDot;\stackrel{\&RightArrow;}{e})+\sqrt{a}]\&times;\stackrel{\&RightArrow;}{e}---\left(8\right)$

Wherein, To put R (x on the separatrix _R, z _R) normal vector located,

The computing of " " expression dot product.

Because R is known, is obtained by formula (8)

After, can try to achieve the coordinate that P is ordered.Along with the R point is mobile in the separatrix, geometric locus

Upper all points all can accurately be obtained.But in actual applications, if at a scanning position place, paint one section curve

After draw again next section curve

The method of this continuous picture geometric locus can comprise repeatedly the double counting of the represented refraction vector of formula (8), and for example among Fig. 1 (d), R is a refraction point of the sound field scope that is positioned on the separatrix, and P and P' are respectively curves

With

On point, at trace

In time, need to calculate once

And trace

In time, need to calculate once equally

Yet

With

Normalized vector identical, so

Calculating be repetition.In fact, for fear of this double counting, we can be at refraction point R place along the refraction vector

Handle all pixels that are positioned at image on the refracted ray, then refraction point is moved to next refraction point R _gCalculate new refraction vector and process pixel on the refracted ray at the place.Be equivalent to like this process simultaneously all geometric locus that the scanning position place is corresponding, but at each refraction point, only drawn a point on each bar geometric locus, after refraction points all in the sound field scope had traveled through, points all on all geometric locuses all can paint.So handle all scanning position points, just obtained entire image.

At refraction point R place along the refraction vector Processing all pixel essence that is positioned at image on the refracted ray is exactly from the R point, with ω corresponding to each pixel at image _is _i(t _i)/r _iValue is along the refraction vector

Draw the straight line section.And in field of Computer Graphics, in order to show straight-line segment at raster display device, developed multiple line scan conversion technology, wherein, most popular is the Bresenham algorithm.The principle of this algorithm as shown in Figure 2, known pixels point (x _p, z _p) be on the straight line a bit, but next pixel has two kinds of selected element P ₁(x _p+ 1, z _p) or P ₂(x _p+ 1, z _p+ 1), the symbol of algorithm use error item ε determines that next pixel gets the right point or lower-right most point:

If straight-line equation is z=kx+b, z is arranged then _P+1=z _p+ k (x _P+1-x _p)=z _p+ k, wherein, k is straight slope.If | k|≤1, can find out that from the coordinate of two candidate points of next pixel the horizontal ordinate of next pixel is x _p+ 1, and ordinate or still be z _p, or increasing progressively 1 is z _p+ 1, whether increase 1 value that depends on error term ε.The initial value of ε is 0, and when the every increase by 1 of horizontal ordinate, the corresponding incremental change of ε value is the slope k of straight line, i.e. ε=ε+k.

If | k|〉1, then should consider x _pThe no incremental change that increases 1, ε value is 1/k.When ε 〉=0.5, the intersection point of straight line and grid is from P ₂Nearer, ordinate increases 1, simultaneously with P ₂As the new reference point that calculates next time, the ε value is corresponding to deduct 1; When ε＜0.5, the intersection point of straight line and grid is from P ₁Nearer, ordinate is constant, and the reference point that next time calculates is also constant.

Calculate for convenient, make δ=ε-0.5, the initial value of δ is-0.5, and increment is k, and when δ 〉=0, next pixel is got (x _p+ 1, z _p+ 1), the δ value reduces 1; When δ＜0, next pixel is got (x _p+ 1, z _p), the δ value is constant.Then for the practical application scene shown in Fig. 1 among the present invention (d), the initial point (x on the Bresenham algorithm cathetus _p, z _p)=(x _R, z _R), k is the refraction vector Slope, the pixel of scanning on the conversion linear is until reach the border of image.So, utilize the Bresenham algorithm can easily scan conversion refraction vector

Be positioned at the straight-line segment of image.

This formula of utilizing is calculated the refraction vector, and utilizes the line scan conversion technology to ask pixel on the refraction path, is a kind of process of finding the solution refraction path of forward with the method for indirect realization picture geometric locus.Thereby the method need not iteration with time-delay as known conditions with refraction point and asks refraction point, also need not to use the root mean square computing to ask acoustic wave propagation path length, has saved a large amount of computing costs.In addition, the method will be asked synthetic aperture effective length L to replace with to ask half-power beam angle β in the former SAFT technology _0.5, and only relevant with the ground floor medium, all the other levels all need not to calculate this value again, so the imaging results of the method can be more accurate in conjunction with the method for ray trace than former SAFT.

Summary of the invention

The object of the invention is to propose a kind of ultrasonic imaging method, realize the fast and accurately imaging to the layering object.

The invention is characterized in, can be to containing horizontal interphase or contain interfacial regular layering object imaging parallel to each other, simultaneously also can be to containing non-level and non-interfacial irregular layering object imaging parallel to each other, and image taking speed is fast.

The invention is characterized in, based on synthetic aperture focusing (SAFT) ultrasonic imaging technique, in conjunction with refraction law, use refraction vector calculation formula forward to calculate refraction path, and utilize the Bresenham algorithm of linear scan conversion to scan each pixel of changing on the refraction path, former SAFT is embodied as the process of using sampled data to draw geometric locus at image in conjunction with the method for ray trace, avoided the root mean square computing asking the interative computation of refraction point and ask the Acoustic Wave Propagation distance, saved a large amount of calculating operations, and substituted the synthetic aperture effective length with the half-power beam angle, the imaging results error that the inaccuracy of having avoided the synthetic aperture effective length to calculate is brought in stratified medium.

The invention is characterized in, contain successively following steps:

Step (1) makes up of being comprised of a computing machine, ultrasonic transducer, a cover register control and analog to digital converter and makes the system of not damaged ultrasonic imaging based on the vertical section for the layering object is formed at the degree of depth and horizontal both direction of synthetic aperture focusing technology and refraction law, wherein:

Described ultrasonic transducer is provided with: the pulse signal input terminal that links to each other with the output terminal of described register control, the input end of described register control links to each other with the corresponding positioning control signal output part of described computing machine, described ultrasonic transducer also is provided with: the echoed signal output terminal that links to each other with the input end of described analog to digital converter, the output terminal of described analog to digital converter links to each other with the echo samples signal input part of described computing machine, described ultrasonic transducer is controlled by described register control, fixed rate on the testee surface with 1 step-length/ms moves, described register control is the gearing of the described ultrasonic transducer of control shift position, its parameter is by described computer input

Testee is X along the horizontal length of X-direction _Length, be divided into X _Length/ Δ x is interval, and Δ x is burst length, also be described ultrasonic transducer along X-axis from coordinate points (0,0) rise to terminal point (X _Length, 0) and end each mobile step-length, the each mobile point that reaches of described ultrasonic transducer is called sensing point, and total M, M=1+X _Length/ Δ x, sequence number m=0,1, M-1, described register control produces a TTL transistor-transistor logic level pulse at each sensing point place, trigger described ultrasonic transducer to driving pulse of depth direction Z emission perpendicular to X-axis of testee, ultrasonic transducer transfers receiving mode to and begins timing subsequently, receives from the echoed signal of testee reflection, the echoed signal that described analog to digital converter receives at sensing point m place described ultrasonic transducer is carried out N sampling and is stored in the computing machine, sampling sequence number n=0,1 ... N-1, sample frequency is f _s, f _sValue be that analog to digital converter is default, note s _m(n) sampled value that obtains for the n time sampling of ultrasonic transducer at m sensing point place, s _m(n) sampling instant is t=n/f _s

Step (2): described computing machine begins sequentially to read the sampled value of sensing point m=0 from n=0, and then, repeat this process and read successively m=1 ..., the sampled value at each sensing point place of M-1;

Step (3): get v=v ₁, v ₁Be ultrasonic velocity of propagation in the ground floor medium of testee, use the ultrasonic imaging by synthetic aperture focusing software package to generate z on the depth direction ₀=0 to Z _DepthThe skiagraph picture in-1 interval, Z _DepthBe the length of default synthetic image, the depth value that the image that namely generates represents with pixel count in the vertical;

Step (4): the z that obtains with step (3) ₀To Z _Depth-1 interval image block uses Canny operator edge extracting software package to extract the separatrix c of ground floor medium and second layer medium as input quantity ₁(x, z);

Step (5): revise according to the following steps on the described vertical section at separatrix c ₁Under (x, z) to Z _DepthThe image in-1 interval, the error that causes to eliminate the not homogeneity between ground floor and other each layer media:

Step (5.1): get m=0, remember that m sensing point is U (x _u, 0), wherein, x _u=m Δ x/accuracy, accuracy is the precision of images, the spacing of adjacent two pixels on the image that namely generates is calculated m sensing point U (x successively according to the following steps _u, 0) and locate corresponding sector image:

Step (5.1.1): the half-power beam angle β that calculates ultrasonic transducer _0.5=0.84 λ/D, λ are ultrasonic wavelength when propagating in testee, and D is the diameter of ultrasonic transducer, and the left and right boundary line of calculating described half-power beam angle respectively with separatrix c ₁Two intersection points B about (x, z) _l(x _l, z _l), B _r(x _r, z _r), wherein, x _l=x _u-z _l* tg (0.5 β _0.5) and satisfy c ₁(x _l, z _l)=0, x _r=x _u+ z _r* tg (0.5 β _0.5) and satisfy c ₁(x _r, z _r)=0;

Step (5.1.2): at separatrix c ₁Get refraction point R (x on (x, z) _R, z _R)=B _l(x _l, z _l), calculate normalized refraction vector

Wherein, eta be in the described testee with separatrix c ₁Relative index of refraction between (x, z) adjacent two layer medium,

To put R (x on the separatrix _R, z _R) the unit normal vector located, $a=1.0-\mathrm{eta}\&times;\mathrm{eta}\&times;[1.0-(\stackrel{\&RightArrow;}{\mathrm{RU}}\&CenterDot;\stackrel{\&RightArrow;}{e})\&times;(\stackrel{\&RightArrow;}{\mathrm{RU}}\&CenterDot;\stackrel{\&RightArrow;}{e})],$ " " is the dot product computing,

For from described sensing point U (x _u, 0) and to refraction point R (x _R, z _R) normalized incident vector;

Step (5.1.3): calculate the refraction vector

Slope k, the initial point on the cut-off line is R (x _R, z _R), utilize Bresenham line scan conversion software package, from refraction point R, calculate all pixels on the refraction path until arrive left margin x=0 or right margin x=(M-1) Δ x/accuracy or the lower boundary z=Z of image _DepthTill-1, remember that current refraction point R is R _f, i.e. R _f(x _f, z _f)=R (x _R, z _R), and the note refraction path is R _fE _f, E _fBe R _fE _fTerminal point in image, i.e. refraction path R _fE _fIntersection point with image boundary;

Step (5.1.4): get x _g=x _f+ 1, search separatrix c ₁(x, z) upper current refraction point R _fThe next pixel R of position _g(x _g, z _g), some R _gCoordinate satisfy c ₁(x _g, z _gR is got in)=0 _gBe new refraction point, i.e. R (x _R, z _R)=R _g(x _g, z _g), execution in step (5.1.2) and step (5.1.3), the refraction path R that the scanning conversion is new _gE _g, the note terminal point is E _g

Step (5.1.5): calculate current refraction point R _gRefraction path R _gE _gTerminal point E _gWith a upper refraction point R _fRefraction path R _fE _fTerminal point E _fBetween distance, delta d=|E _fE _g|, with separatrix c ₁(x, z) upper current refraction point R _gWith a upper refraction point R _fBetween segment of curve be divided into Δ d part, namely insert Δ d-1 point, the insertion point sequence number is designated as τ=1,2 ..., Δ d-1, to the τ value take 1 as step-length take τ=1 carry out following step (5.1.5.1) until τ=Δ d ends as the initial value circulation:

Step (5.1.5.1): τ insertion point R of interpolation calculation _τCoordinate, x _τ=x _f+ τ/Δ d, z _τ=z _f+ (z _g-z _f) τ/Δ d, get refraction point R (x _R, z _R)=R _τ(x _τ, z _τ), execution in step (5.1.2) and step (5.1.3), the scanning conversion is from a R _τThe refraction path that sets out;

Step (5.1.6): get R _f=R _g, E _f=E _g, execution in step (5.1.4) and step (5.1.5);

Step (5.1.7): repeated execution of steps (5.1.6) is until x _g=x _r+ 1 ends, and namely handles separatrix c ₁(x, z) is upper between B _l(x _l, z _l) and B _r(x _r, z _r) between all refraction point and refraction paths, obtain m sensing point U (x _u, 0) and locate corresponding sector image;

Step (5.2): get successively m=1 ..., M-1, repeated execution of steps (5.1) generates on the described vertical section at separatrix c ₁Under (x, z) to Z _DepthThe image in-1 interval;

Step (6): the separatrix c that obtains with step (5) ₁Under (x, z) to Z _DepthThe image block in-1 interval uses Canny operator edge extracting software package to extract described separatrix c as input quantity ₁(x, z) is to Z _DepthThe separatrix c of the second layer medium between-1 and the 3rd layer of medium ₂(x, z), according on the described vertical section of the described method correction of step (5) at separatrix c ₂Under (x, z) to Z _DepthThe image in-1 interval, the error that causes to eliminate the not homogeneity between the second layer and its lower each layer medium;

Step (7): by the described method of step (6), process each bar separatrix of residue, until only handle all separatrix in the described vertical section, generating width is that x/accuracy+1 pixel of (M-1) Δ, length are Z _DepthThe image of the described vertical section of individual pixel.

The present invention compared with prior art, advantage is and can carries out imaging to containing regular interphase or containing irregular interfacial layering object, image taking speed is fast and imaging is more accurate.For example: for the testee shown in Fig. 3 (a), if the ultrasonic transducer diameter is 0.5mm, the ultrasonic transducer moving step length is 0.7mm, the ultrasound wave centre frequency of ultrasonic transducer emission is 5MHz, sample frequency 100MHz, imaging precision accuracy is taken as 0.05mm, utilize method of the present invention to generate Fig. 3 (d) at the experiment machine of Intel Core Duo2.66GHz CPU, 2.0GB RAM and only need 16s, be about former SAFT under the same case in conjunction with the used time of ray tracing method 1/112, the Phase shift technology 1/4.The maximum error of upper surface is 0.5mm (such as Fig. 4 (a)) among Fig. 3 (d), and the maximum error of lower surface is 0.9mm (such as Fig. 4 (b)), and imaging precision is higher.And can find that with Fig. 3 (b) and Fig. 3 (c) contrast method of the present invention is better than front two kinds of methods to the successful of lower surface imaging.

Description of drawings

Fig. 1 is working model and the schematic diagram of SAFT ultrasonic imaging technique: 1 (a) is that testee is in the situation of Single Medium, the reverse computation process schematic diagram of SAFT technology; 1 (b) is that testee is in the situation of stratified medium, the reverse computation process schematic diagram of SAFT technology; 1 (c) is that testee is in the situation of stratified medium, and the forward of SAFT technology is explained key diagram; 1 (d) is that testee is in the situation of stratified medium, and ultrasonic transducer is at scanning position u _iThe effect figure of all data in former SAFT technology that place's sampling obtains.

Fig. 2 is pixel and the error term synoptic diagram that the per step iteration of Bresenham algorithm of linear scan conversion relates to.

Fig. 3 is the comparison diagram of each ultrasonic imaging method under the identical experiment environment: 3 (a) are the longitudinal diagram of the testee of irregular layering; 3 (b) are the image that SAFT generates in conjunction with Tracing Technology, and imaging time 30 minutes is not good enough to the effect of lower surface imaging; The image that 3 (c) generate for the Phase shift ultrasonic imaging technique, imaging time 62s gross error occurs to the lower surface imaging; The image that 3 (d) generate for this formation method, imaging time 16s, best to the effect of lower surface imaging, among the figure the white curve in the upper and lower surperficial imaging is respectively the upper and lower separatrix that experiment obtains.

Fig. 4 is the analysis of experimental data figure to Fig. 3 (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, the marginal graph of errors of the 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, the second layer separatrix that the imaginary point curve obtains for this formation method, 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 synoptic diagram.

Fig. 8 is the algorithm flow chart that the ultrasonic imaging by synthetic aperture focusing of layering object calculates.

Embodiment

Specific implementation process of the present invention comprises three parts (such as Fig. 5): ultrasound data obtains, imaging calculating and image show.The hardware platform system structural drawing as shown in Figure 6, ultrasonic image-forming system is comprised 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. surface of the direct contact measured object of ultrasonic transducer; Perhaps liquid immersion type, namely object under test places liquid, and ultrasonic transducer is measured at liquid surface.According to the difference of detection mode, ultrasonic transducer can use contact probe or liquid immersion type probe.System uses single transmitting/receiving ultrasonic transducer, and ultrasonic transducer moves with about 1 step/ms fixed rate along directions X (such as Fig. 7) with uniform step delta x in testee surface or liquid by register control.TTL(transistor-transistor logic level of moment generation that controller is stable in each target location) pulse, this pulse is used for triggering ultrasonic transducer to the perpendicular driving pulse of depth direction Z emission of testee and directions X, ultrasonic transducer transfers receiving mode to and begins timing subsequently, receives from the echo of testee reflection.Each position of ultrasonic transducer transponder pulse and reception echo is sensing point.The echoed signal that ultrasonic transducer receives is by the analog to digital converter collection and be stored in the storer.Ultrasonic transducer moving step length Δ x needs to require to determine that comprehensively its value is less that the image of generation is more accurate according to the actual size of object under test and imaging precision, but computing time is also longer.

Imaging calculate be exactly sampled data with testee each sensing point place on a vertical section as computer input, then by the skiagraph picture of aforesaid image-forming step calculating testee, overview flow chart is referring to Fig. 8.

In the specific implementation program, the ultrasonic transducer moving step length Δ x in the aforementioned image-forming step (1) can be converted into pixel count, only needs to carry out Δ x ← Δ x/accuracy and gets final product, and then can directly calculate abscissa value x in the step (5.1) _u=m Δ x, and avoid need to divided by the precision of images to be converted into the pixel coordinate on the image, saving the division arithmetic expense to each m.

Because in advance and do not know marginal accurate information in the layering object,, step (3) extracts the separatrix between ground floor and second layer medium so using the ultrasonic imaging algorithm of Single Medium to generate piece image.Because what use in the algorithm is the velocity of sound in the ground floor medium, so the image that the ground floor medium is become is accurately, and the separatrix between ground floor and second layer medium also belongs to the part of ground floor medium, so this separatrix of extracting in the step (4) also is accurately.In fact, the ultrasonic imaging algorithm that uses in the step (3) both can be former SAFT technology, also can be that Phase shift technology etc. is any can be to the technology of Single Medium imaging.

Need to use the separatrix between the Canny operator Boundary extracting algorithm extraction medium in the image-forming step (4).The computation process of Canny arithmetic operators was divided into for four steps:

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

Step (b): the gradient G of each pixel of computed image I (x, z) and direction F.Adopt 2 * 2 templates as to directions X and

The first approximation of the partial differential of Z direction, namely

$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 G and direction F are

$G=\sqrt{p\&times;p+q\&times;q},$ $F=\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 G (x, z) illustrates that less than the gradient magnitude along two consecutive point of its gradient direction F (x, z) this point is not marginal point, then the gray scale with I (x, z) is made as 0;

Step (d): the dual threashold value is processed.Set the dual threshold method and detect and the low threshold value Low that is connected the edge needs and high threshold High, gradient image is carried out the dual threashold value process.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 the edge pixel that whether exists in eight neighbors of this pixel greater than high threshold High, then is edge pixel if exist, otherwise is not.

In whole imaging is calculated, the half-power beam angle β in the image-forming step (5.1.1) _0.5Need to use M time, however should value and the intrinsic relating to parameters of ultrasonic transducer, and the present invention only uses a ultrasonic transducer to survey, and is fixed value so this is worth, and for fear of double counting, can calculate in advance β in program _0.5← tg (0.5 * 0.84 λ/D), and save as global variable, then in step (5.1.1), directly call this variable and calculate x _l=x _u-z _l* β _0.5And x _r=x _u+ z _r* β _0.5

In addition, the left and right boundary line that needs to calculate the half-power beam angle in the step (5.1.1) respectively with separatrix c ₁The intersection points B of (x, z) _l(x _l, z _l), B _r(x _r, z _r).In reality realizes, the separatrix c of Canny operator extraction ₁(x, z) saves as the pixel point set with x coordinate figure growth order, calculating intersection points B _lAnd B _rThe time, only need these pixels of sequential search to get final product, if separatrix c ₁The coordinate of certain pixel (x, z) on (x, z) satisfies x=x _u-z * β _0.5Then be B _lIf satisfy x=x _u+ z * β _0.5Then be B _rSince in whole imaging is calculated, to each m, two intersection points B about all need calculating _lAnd B _r, but the corresponding half-power beam of each m angle is (such as Fig. 1 (c)) parallel to each other, the then corresponding left intersection points B of current m _lMust be in the right of the left intersection point of a upper m, the corresponding right intersection points B of current m _rAlso certainty, only needs from last round of B so each is taken turns when calculating left intersection point in the step (5.1.1) on the right of the right intersection point of a upper m _lThe place begins to search for separatrix c ₁Pixel on (x, z), corresponding, each is taken turns when calculating right intersection point, only needs from last round of B _rThe place begins search.

The specific algorithm of imaging computation process is divided into the ultrasonic imaging by synthetic aperture focusing algorithm master routine SAFT-Refraction () of Bresenham line scan conversion subroutine BresenhamLine () and layering object:

Image show to be about to the two-dimensional image data that the imaging calculation stages obtains and to be presented on the display, as required display gray scale image or coloured image.

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

Claims (1)

1. the ultrasonic imaging by synthetic aperture focusing method of layering object is characterized in that, contains successively following steps:

Step (1): make up of being formed by a computing machine, ultrasonic transducer, a cover register control and analog to digital converter and make the system of not damaged ultrasonic imaging based on the vertical section for the layering object is formed at the degree of depth and horizontal both direction of synthetic aperture focusing technology and refraction law, wherein:

Described ultrasonic transducer is provided with: the pulse signal input terminal that links to each other with the output terminal of described register control, the input end of described register control links to each other with the corresponding positioning control signal output part of described computing machine, described ultrasonic transducer also is provided with: the echoed signal output terminal that links to each other with the input end of described analog to digital converter, the output terminal of described analog to digital converter links to each other with the echo samples signal input part of described computing machine, described ultrasonic transducer is controlled by described register control, fixed rate on the testee surface with 1 step-length/ms moves, described register control is the gearing of the described ultrasonic transducer of control shift position, its parameter is by described computer input

Testee is X along the horizontal length of X-direction _Length, be divided into X _Length/ Δ x is interval, and Δ x is burst length, also be described ultrasonic transducer along X-axis from coordinate points (0,0) rise to terminal point (X _Length, 0) and end each mobile step-length, the each mobile point that reaches of described ultrasonic transducer is called sensing point, and total M, M=1+X _Length/ Δ x, sequence number m=0,1, M-1, described register control produces a TTL transistor-transistor logic level pulse at each sensing point place, trigger described ultrasonic transducer to driving pulse of depth direction Z emission perpendicular to X-axis of testee, ultrasonic transducer transfers receiving mode to and begins timing subsequently, receives from the echoed signal of testee reflection, the echoed signal that described analog to digital converter receives at sensing point m place described ultrasonic transducer is carried out N sampling and is stored in the computing machine, sampling sequence number n=0,1 ... N-1, sample frequency is f _s, f _sValue be that analog to digital converter is default, note s _m(n) sampled value that obtains for the n time sampling of ultrasonic transducer at m sensing point place, s _m(n) sampling instant is t=n/f _s

Step (2): described computing machine begins sequentially to read the sampled value of sensing point m=0 from n=0, and then, repeat this process and read successively m=1 ..., the sampled value at each sensing point place of M-1;

Step (3): get v=v ₁, v ₁Be ultrasonic velocity of propagation in the ground floor medium of testee, use the ultrasonic imaging by synthetic aperture focusing software package to generate z on the depth direction ₀=0 to Z _DepthThe skiagraph picture in-1 interval, Z _DepthBe the length of default synthetic image, the depth value that the image that namely generates represents with pixel count in the vertical;

Step (4): the z that obtains with step (3) ₀To Z _Depth-1 interval image block uses Canny operator edge extracting software package to extract the separatrix c of ground floor medium and second layer medium as input quantity ₁(x, z);

Step (5): revise according to the following steps on the described vertical section at separatrix c ₁Under (x, z) to Z _DepthThe image in-1 interval, the error that causes to eliminate the not homogeneity between ground floor and other each layer media:

Step (5.1): get m=0, remember that m sensing point is U (x _u, 0), wherein, x _u=m Δ x/accuracy, accuracy is the precision of images, the spacing of adjacent two pixels on the image that namely generates is calculated m sensing point U (x successively according to the following steps _u, 0) and locate corresponding sector image:

Step (5.1.1): the half-power beam angle β that calculates ultrasonic transducer _0.5=0.84 λ/D, λ are ultrasonic wavelength when propagating in testee, and D is the diameter of ultrasonic transducer, and the left and right boundary line of calculating described half-power beam angle respectively with separatrix c ₁Two intersection points B about (x, z) _l(x _l, z _l), B _r(x _r, z _r), wherein, x _l=x _u-z _l* tg (0.5 β _0.5) and satisfy c ₁(x _l, z _l)=0, x _r=x _u+ z _r* tg (0.5 β _0.5) and satisfy c ₁(x _r, z _r)=0;

Step (5.1.2): at separatrix c ₁Get refraction point R (x on (x, z) _R, z _R)=B _l(x _l, z _l), calculate normalized refraction vector Wherein, eta be in the described testee with separatrix c ₁Relative index of refraction between (x, z) adjacent two layer medium,

To put R (x on the separatrix _R, z _R) the unit normal vector located, $a=1.0-\mathrm{eta}\&times;\mathrm{eta}\&times;[1.0-(\stackrel{\&RightArrow;}{\mathrm{RU}}\&CenterDot;\stackrel{\&RightArrow;}{e})\&times;(\stackrel{\&RightArrow;}{\mathrm{RU}}\&CenterDot;\stackrel{\&RightArrow;}{e})],$ " " is the dot product computing,

For from described sensing point U (x _u, 0) and to refraction point R (x _R, z _R) normalized incident vector;

Step (5.1.3): calculate the refraction vector

Slope k, the initial point on the cut-off line is R (x _R, z _R), utilize Bresenham line scan conversion software package, from refraction point R, calculate all pixels on the refraction path until arrive left margin x=0 or right margin x=(M-1) Δ x/accuracy or the lower boundary z=Z of image _DepthTill-1, remember that current refraction point R is R _f, i.e. R _f(x _f, z _f)=R (x _R, z _R), and the note refraction path is R _fE _f, E _fBe R _fE _fTerminal point in image, i.e. refraction path R _fE _fIntersection point with image boundary;

Step (5.1.4): get x _g=x _f+ 1, search separatrix c ₁(x, z) upper current refraction point R _fThe next pixel R of position _g(x _g, z _g), some R _gCoordinate satisfy c ₁(x _g, z _gR is got in)=0 _gBe new refraction point, i.e. R (x _R, z _R)=R _g(x _g, z _g), execution in step (5.1.2) and step (5.1.3), the refraction path R that the scanning conversion is new _gE _g, the note terminal point is E _g

Step (5.1.5): calculate current refraction point R _gRefraction path R _gE _gTerminal point E _gWith a upper refraction point R _fRefraction path R _fE _fTerminal point E _fBetween distance, delta d=|E _fE _g|, with separatrix c ₁(x, z) upper current refraction point R _gWith a upper refraction point R _fBetween segment of curve be divided into Δ d part, namely insert Δ d-1 point, the insertion point sequence number is designated as τ=1,2 ..., Δ d-1, to the τ value take 1 as step-length take τ=1 carry out following step (5.1.5.1) until τ=Δ d ends as the initial value circulation:

Step (5.1.5.1): τ insertion point R of interpolation calculation _τCoordinate, x _τ=x _f+ τ/Δ d, z _τ=z _f+ (z _g-z _f) τ/Δ d, get refraction point R (x _R, z _R)=R _τ(x _τ, z _τ), execution in step (5.1.2) and step (5.1.3), the scanning conversion is from a R _τThe refraction path that sets out;

Step (5.1.6): get R _f=R _g, E _f=E _g, execution in step (5.1.4) and step (5.1.5);

Step (5.1.7): repeated execution of steps (5.1.6) is until x _g=x _r+ 1 ends, and namely handles separatrix c ₁(x, z) is upper between B _l(x _l, z _l) and B _r(x _r, z _r) between all refraction point and refraction paths, obtain m sensing point U (x _u, 0) and locate corresponding sector image;

Step (5.2): get successively m=1 ..., M-1, repeated execution of steps (5.1) generates on the described vertical section at separatrix c ₁Under (x, z) to Z _DepthThe image in-1 interval;

Step (6): the separatrix c that obtains with step (5) ₁Under (x, z) to Z _DepthThe image block in-1 interval uses Canny operator edge extracting software package to extract described separatrix c as input quantity ₁(x, z) is to Z _DepthThe separatrix c of the second layer medium between-1 and the 3rd layer of medium ₂(x, z), according on the described vertical section of the described method correction of step (5) at separatrix c ₂Under (x, z) to Z _DepthThe image in-1 interval, the error that causes to eliminate the not homogeneity between the second layer and its lower each layer medium;

Step (7): by the described method of step (6), process each bar separatrix of residue, until only handle all separatrix in the described vertical section, generating width is that x/accuracy+1 pixel of (M-1) Δ, length are Z _DepthThe image of the described vertical section of individual pixel.

Images