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

Synthetic Aperture Focused Ultrasound Imaging Method for Layered Objects

technical field

The invention relates to ultrasonic non-destructive testing technology, synthetic aperture focusing imaging technology, refraction law and linear scanning conversion technology, and realizes rapid imaging of the surface shape and internal structure of layered objects containing regular or irregular interfaces.

Background technique

At present, in the field of ultrasonic nondestructive testing, for ultrasonic imaging of layered objects containing interfaces (ultrasonic testing methods include contact and liquid immersion modes, if the coupling agent layer is regarded as a part of the measured object, liquid immersion Ultrasonic detection can also be regarded as contact ultrasonic detection of layered objects, so in the present invention, the two are collectively referred to as ultrasonic detection and imaging of layered objects), there are mainly two schemes: one is to use synthetic aperture focusing ( SAFT) ultrasonic imaging technology combined with Ray Tracing (Ray Tracing) technology, the second is the phase shift ultrasonic imaging method based on phase shift (Phase Shift Migration) technology.

Synthetic Aperture Focusing Technology (SAFT) originated from Synthetic Aperture Radar Technology (SAR) and was introduced into the field of ultrasound imaging in the early 1970s. The transducer size has nothing to do with distance and other advantages. The basic principle of SAFT ultrasonic imaging technology is to use the pulse-echo (pulse-echo) measurement mechanism, use an ultrasonic transducer to scan the measured object in an orderly manner along a fixed trajectory, and use the delay superposition (DAS) method The pulse-echo signal obtained by scanning is focused and imaged to achieve the purpose of simulating a large aperture array by using a single ultrasonic transducer with a smaller aperture. The working model of SAFT ultrasonic imaging is shown in Figure 1(a). The ultrasonic transducer moves equidistantly along the scanning direction (X direction) on the surface of the object, and emits light in the depth direction (Z direction) of the object at each scanning position. Ultrasonic signal, to detect the inside of the object, while the ultrasonic transducer receives the echo signal reflected by the internal reflector of the object and samples and saves it, and finally performs DAS and multi-point dynamic focusing processing on the sampled data obtained at all scanning positions and displays the image .

DAS and multi-point dynamic focusing technology need to calculate different delay curves for the focus points at different depths of the target imaging area. As shown in Figure 1(a), in order to focus on the target reflector (x, z), the SAFT technique delays the echo signal obtained by the ultrasonic transducer at each scanning position within the effective length L of the synthetic aperture. Time superposition processing: Let s _i (t) be the echo signal received by the ultrasonic transducer at u _i , t is the sampling time, and the delay at u _i about the target reflector (x, z) is

${\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}{\mathrm{<span\; class="notranslate">\; v</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; <figure-callout\; id="2"\; label="v\; z"\; filenames="HDA00002540433900031.png"\; state="\{\{state\}\}">v</figure-callout></span><figure-callout\; id="2"\; label="v\; z"\; filenames="HDA00002540433900031.png"\; state="\{\{state\}\}">\; </figure-callout>}}\sqrt{{\mathrm{<span\; class="notranslate">\; </span>}}^{}}\sqrt{{\mathrm{<span\; class="notranslate">z</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0,1</span>}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; .</span><span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; .</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$

Among them, v is the propagation speed of ultrasound in the medium, and _ri is the straight-line distance between point (x, z) and u _i . All the time delays within the effective length L of the synthetic aperture constitute a time delay curve, which is a hyperbola. The formula for calculating L is

L=0.84λz/D (2)

λ is the wavelength of ultrasound in the medium, D is the diameter of the ultrasonic transducer, then the imaging at (x,z) is

$\mathrm{<span\; class="notranslate">\; I</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; z</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; L</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\frac{{\mathrm{<span\; class="notranslate">\; \&omega;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span>}{{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$

Among them, ω _i is the apodization function.

It can be seen from the principle of DAS that the SAFT algorithm needs to determine the delay according to the relative distance between the ultrasonic transducer and the reflector and the wave speed. When the measured object is layered, as shown in Figure 1(b), the ultrasonic signal emitted by the ultrasonic transducer will be refracted at the interface, the propagation path of the sound wave will change, and the target reflector (x, z) The distance r _i from u _i is no longer the length of the straight line segment between two points. In addition, the medium of each layer of a layered object is different, and the propagation speed of sound waves in each layer is usually not equal. Therefore, the calculation of time delay in DAS can no longer use the simple calculation method of dividing the relative distance by the wave speed, and the time delay curve is also different. No longer a regular hyperbolic shape. In order to obtain the delay time, it is necessary to first obtain the propagation path of the sound wave between u _i and the target reflector (x, z), and then calculate the length and propagation time of each path segment segment by segment (as shown in Figure 1(b) r _i1 , r _i2 and t _i ). The important and difficult point is to find out the propagation path of the beam quickly and accurately, and the ray tracing technology is the best way to realize this purpose. Therefore, ray tracing technology is naturally combined into SAFT technology to realize ultrasonic imaging of layered objects.

The principle of ray tracing is mainly based on Snell's theorem or Fermat's principle, through iterative calculation to find out the shortest sound wave propagation path. Therefore, as shown in Figure 1(b), in the method of SAFT combined with ray tracing, when calculating the acoustic wave propagation delay t _i between u _i and the target reflector (x, z), it is necessary to follow the dividing line c(x , z)=0 to iterate all the points on the line, calculate the length and propagation time of the two propagation paths r _i1 and r _i2 corresponding to each point, and find the point that makes the propagation time the shortest, which is refraction point.

The advantage of this method is that ray tracing technology has no special requirements on the shape of the interface of layered objects, and it is not only applicable to the presence of heterogeneous media in the depth direction but the same type of media in the horizontal direction, that is, the separation between layer and layer media. Ultrasonic imaging of regularly layered objects with interfaces that are horizontal or parallel to each other can also be applied to heterogeneous media in both the depth direction and the horizontal direction, that is, irregular layers that contain irregular interfaces (such as curved surfaces) Ultrasound imaging of objects. But its serious disadvantage is that the ray tracing algorithm contains iterative operations, which has high time complexity. In addition, the formula (1) for calculating the time delay in the original SAFT technology contains a root mean square operation, and the calculation cost is already very large. When the ray tracing algorithm is introduced After the tracking technology, the calculation of time delay includes multiple root mean square calculations, as shown in Figure 1(b), when calculating the acoustic wave propagation delay t _i between u _i and the target reflector (x, z), it is necessary For each point on the boundary line, the root mean square is calculated twice to obtain the lengths of the two paths r _i1 and r _i2 , and the calculation overhead is larger. Therefore, this method is very time-consuming. As shown in Figure 3(b), it takes 30 minutes to image the measured object (4.4cm*4.9cm) shown in Figure 3(a) with this method. In addition, because the wavelength of the ultrasonic wave is related to the medium, when the medium of the measured object is layered, the calculation formula (2) of the effective length L of the synthetic aperture is no longer applicable. However, the accuracy of the imaging of the SAFT technology is closely related to the accuracy of L Related, if L is too large, there will be more imaging noise and low signal-to-noise ratio; if L is too small, the details of the measured object will be lost in the imaging, and the accuracy will not be high. However, for complex and irregularly layered objects, there is currently a lack of an accurate formula for calculating L, and it is difficult to find an accurate L value in the method of SAFT combined with ray tracing, resulting in poor imaging effects of this method, as shown in Figure 3(b) In this method, the image formed by the longitudinal section of the measured object cannot reconstruct the curved surface part in the second layer interface well.

The Phase Shift Migration ultrasonic imaging method is an ultrasonic imaging method in the frequency domain obtained by introducing the Migration Technique in Reflection Seismology into the field of ultrasonic imaging. This method regards the ultrasonic detection model as an explosive reflection model, assuming that the reflectors in the object to be tested explode at the same time at t=0, and the explosion intensity of each reflector is proportional to its reflectivity. energy meter to measure. The main idea is to extrapolate the sound field observed from the horizontal position (ie, the first row in the depth direction) to calculate the sound field at other positions in the depth direction. The specific algorithm includes two main steps: the first step is to perform two-dimensional Fourier transform on the time domain data to obtain the two-dimensional spectrum,

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

The second step is the cycle in the depth direction. First, phase shift the two-dimensional spectrum obtained in the previous cycle, and then perform two-dimensional inverse Fourier transform and take t0 to obtain a row of images.

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

为深度方向上Δz步长所对应的相位迁移量，c为超声在介质中的传播速度，其取值恒定。in,

is the phase shift corresponding to the Δz step in the depth direction, and c is the propagation speed of ultrasound in the medium, and its value is constant.

For the imaging of multi-layer objects, the phase shift ultrasonic imaging method does not require a ray path, and the imaging speed can be greatly improved. For example, it only takes 62s to generate the image shown in Figure 3(c). However, since the wave velocity of the sound wave is assumed to be constant in the phase shift technique, this method can only be applied to ultrasonic imaging of regularly layered objects with heterogeneous media in the depth direction and the same kind of media in the horizontal direction. There are heterogeneous media in both the depth direction and the horizontal direction of the object, which leads to inaccurate detection results. As shown in Figure 3(c), serious errors have occurred in the imaging results of the second layer interface.

Therefore, there is currently a lack of rapid, accurate and effective methods for ultrasound imaging of complex layered objects with irregular interfaces.

From the principle of SAFT combined with ray tracing technology, it can be seen that the most time-consuming part of this method is the iterative operation step of finding the optimal refraction point. From the ray tracing technology and the principle of DAS, it can be known that this method is a reverse calculation process. It is necessary to first determine a pixel point on the image, and then find out all the scanning position points corresponding to this point, and then calculate each scanning position point in the Iteratively find the optimal refraction point on the boundary to obtain the distance and delay between the image point and the scanning position point, and finally accumulate the data points on the delay curve to obtain the pixel value of the image point. Since the specific position of the target reflector in the measured object is not known in advance, it is necessary to use dynamic focusing technology to perform the reverse calculation process on all pixels in the image to generate the entire image, so that the pixels corresponding to the reflector At the pixel point, the ultrasonic signal is superimposed consistently to maximize the cumulative intensity and achieve focus, while at other pixel points, the superposition of the ultrasonic signal is chaotic, and it is difficult to maximize the cumulative intensity. Therefore, in the image, the accumulated pixel value corresponding to the reflective object is significantly larger than other pixel points.

是u_i处超声换能器发射的声场范围内的一段曲线，该曲线上的点距u_i的传播时间t_i均相同，但距离r_i不一定相同。在第一层介质中，超声换能器在u_i处的声场位于该换能器的半功率波束角β_0.5内（即介于图中从u_i处出发的两条虚点线之间），由于合成孔径有效长度L同时也定义为L=z×β_0.5，结合公式(2)，因此β_0.5=0.84λ/D。而在第二层介质中，由于介质发生改变，声波在分界面处发生折射，原半功率波束角β_0.5的两条边界线也会发生折射（如图1(c)中的虚点线），从而u_i处超声换能器的声场范围也会发生变化，但声场仍位于两条虚点线所表示的边界线之间。所以，曲线段可以看作是当折射点沿着第一层与第二层介质之间的分界线在声场范围内从左到右移动时，超声换能器在u_i处的采样数据s_i(t_i)在图像上所走过的轨迹。当第二层介质与第一层的相同时，则无分层，此时该轨迹为一段以u_i为圆心、角度为β_0.5的圆弧；而当两层介质不同时，该轨迹的形状则与分界线的形状有关。If the calculation process of all pixels on the image is considered as a whole, it can be found that each data sampled at each scanning position not only acts on the pixel corresponding to the reflector, but also acts on the pixel without the reflector. As shown in Figure 1(c), the sampling data s _i (t _i ) at u _i not only participates in the imaging calculation of point (x, z), but also participates in the curve segment calculations on other points.

is a curve within the range of the sound field emitted by the ultrasonic transducer at u _i , and the propagation time t _i of the points on the curve from u _i are all the same, but the distance r _i is not necessarily the same. In the first layer of medium, the sound field of the ultrasonic transducer at u _i is located within the half-power beam angle β _0.5 of the transducer (that is, between the two dotted lines starting from u _i in the figure) , since the effective length L of the synthetic aperture is also defined as L=z×β _0.5 , combined with formula (2), therefore β _0.5 =0.84λ/D. In the second layer of medium, due to the change of the medium, the sound wave is refracted at the interface, and the two boundary lines with the original half-power beam angle β _0.5 will also be refracted (as shown in the dotted line in Figure 1(c)) , so that the sound field range of the ultrasonic transducer at u _i will also change, but the sound field is still located between the boundary lines indicated by the two dotted lines. So, the curve segment It can be seen as when the refraction point moves from left to right in the sound field along the boundary line between the first layer and the second layer, the sampling data s _i (t _i ) of the ultrasonic transducer at u _i The trajectory traveled on the image. When the second layer of medium is the same as that of the first layer, there is no delamination. At this time, the trajectory is a circular arc with u _i as the center and an angle of β _0.5 ; and when the two layers of medium are different, the shape of the trajectory It is related to the shape of the dividing line.

即曲线上各像素点的像素值

在原L范围内，s_i(t_i)所在的延时曲线上的所有数据（s₀(t₀)至s_L-1(t_L-1)）所对应的曲线段（

至

）在图像中的交点即为原逆向计算过程中的目标反射物点(x,z)。Therefore, from a positive comprehensive understanding of the overall calculation process, the effect of sampling data s _i (t _i ) in the imaging calculation of the entire image is equivalent to using the data value ω _i s _i (t _i )/r _i in the image draw a curve

That is, the pixel value of each pixel point on the curve

Within the original L range, the curve segment corresponding to all the data (s ₀ ( _{t 0 )} to s _L-1 (t _L-1 )) on the delay curve where s _i (t i ) is located (

to

) in the image is the target reflector point (x,z) in the original reverse calculation process.

If only the effect of all sampling data s _i (t) at u _i is considered in the entire image, it corresponds to an irregular fan-shaped graph I _ui centered on u _i (as shown in Figure 1(d)), namely

${\mathrm{<span\; class="notranslate">\; I</span>}}_{{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}<span\; class="notranslate">\; =</span>\underset{{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}\mathrm{<span\; class="notranslate">\; I</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; tra</span>}}_{{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 4</span>}<span\; class="notranslate">\; )</span>$

N is the number of samples of the ultrasonic transducer at each scanning position. The DAS calculation formula (3) in the original SAFT combined with ray tracing method can be re-understood as the superposition of the fan diagrams corresponding to each scanning position, namely:

$\mathrm{<span\; class="notranslate">\; I</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; z</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; I</span>}}_{{\mathrm{<span\; class="notranslate">\; u</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; z</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 5</span>}<span\; class="notranslate">\; )</span>$

Among them, M is the total number of ultrasonic transducer scanning position points.

需要求得该曲线的函数表达式。如图1(c)，假设点P(x_P,z_P)是轨迹曲线上的任意一点，超声波从超声换能器所在位置U(u_i,0)处发射，在分界线上的点R(x_R,z_R)处折射后到达P点所用的时间为t_i/2，即Therefore, the original SAFT method combined with ray tracing can be realized by drawing a trajectory curve on the image with each sampling data at each scanning position. In order to accurately calculate the trajectory curve

The function expression of the curve needs to be obtained. As shown in Figure 1(c), it is assumed that the point P(x _P , z _P ) is a trajectory curve At any point on , the ultrasonic wave is emitted from the position U(u _i ,0) where the ultrasonic transducer is located, and the time it takes to reach point P after being refracted at point R(x _R ,z _R ) on the dividing line is t _i / 2, namely

$<span\; class="notranslate">\; |</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; RU</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; /</span>{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; the\; water</span>}}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; |</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; RP</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; /</span>{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; object</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 6</span>}<span\; class="notranslate">\; )</span>$

Among them, v _water is the propagation speed of sound wave in water, and v _object is the propagation speed of sound wave in the measured object. Meanwhile, according to the law of refraction:

eta=sinα _i /sinα _r =v _water /v _object (7) where α _i is the incident angle, α _r is the refraction angle, and eta is the refractive index. The coordinates of point P can be obtained from formula (6) and formula (7), that is, the trajectory curve function, but for easier calculation, you can use the refraction formula in computer graphics to calculate the refraction vector

$\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; RP</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; eta</span>}<span\; class="notranslate">\; \&times;</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; RU</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; [</span>\mathrm{<span\; class="notranslate">\; eta</span>}<span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; (</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; RU</span>}}<span\; class="notranslate">\; \&Center\; Dot;</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\sqrt{\mathrm{<span\; class="notranslate">\; a</span>}}<span\; class="notranslate">\; ]</span><span\; class="notranslate">\; \&times;</span>\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 8</span>}<span\; class="notranslate">\; )</span>$

“·”表示向量点积运算。in, is the normal vector at point R(x _R , z _R ) on the dividing line,

"·" means vector dot product operation.

后，即可求得P点的坐标。随着R点在分界线上移动，轨迹曲线

上所有的点都可准确求出。但是在实际应用中，若在一个扫描位置处，画完一段曲线

后再画下一段曲线

这种连续画轨迹曲线的方法会包含多次式(8)所表示的折射向量的重复计算，例如图1(d)中，R是分界线上的位于声场范围内的一个折射点，P和P'分别是曲线

和

上的点，在画曲线

时需要计算一次

而画曲线

时同样需要计算一次

然而

和

的归一化矢量相同，所以

的计算是重复的。事实上，为了避免这种重复计算，我们可以在折射点R处沿着折射矢量

处理完折射线上位于图像内的所有的像素点，然后将折射点移至下一个折射点R_g处，计算新的折射矢量并处理折射线上像素点。这样相当于同时处理一个扫描位置处对应的所有的轨迹曲线，但在每一个折射点，只画了各条轨迹曲线上的一个点，当声场范围内所有的折射点遍历完后，所有的轨迹曲线上所有的点都会画完。如此处理完所有的扫描位置点，就获得了整幅图像。Since R is known, it can be obtained from formula (8)

After that, the coordinates of point P can be obtained. As point R moves on the dividing line, the trajectory curve

All points above can be obtained exactly. However, in practical applications, if a section of the curve is drawn at a scanning position

Then draw the next curve

This method of continuously drawing trajectory curves will include repeated calculations of the refraction vector represented by the multiple equation (8). For example, in Figure 1(d), R is a refraction point on the boundary line within the range of the sound field, and P and P' are the curves

and

points on the curve

needs to be calculated once

while drawing a curve

also needs to be calculated once

However

and

The normalized vectors of are the same, so

calculations are repeated. In fact, to avoid this double counting, we can follow the refraction vector at the refraction point R

After processing all the pixels on the refraction line in the image, then move the refraction point to the next refraction point _Rg , calculate a new refraction vector and process the pixels on the refraction line. This is equivalent to processing all the trajectory curves corresponding to a scanning position at the same time, but at each refraction point, only one point on each trajectory curve is drawn. After all the refraction points in the sound field have been traversed, all the trajectory All points on the curve will be drawn. After processing all the scanning position points in this way, the whole image is obtained.

画一条直线段。而在计算机图形学领域，为了在光栅显示设备上显示直线段，已开发了多种直线扫描转换技术，其中，使用最广泛的是Bresenham算法。该算法的原理如图2所示，已知像素点(x_p,z_p)为直线上一点，下一个像素点有两种可选择点P₁(x_p+1,z_p)或P₂(x_p+1,z_p+1)，算法使用误差项ε的符号决定下一个像素取右边点还是右下点：At the refraction point R along the refraction vector The essence of processing all the pixels located in the image on the refraction line is to start from point R in the image, and use the value of ω _i s _i (t _i )/r _i corresponding to each pixel, along the refraction vector

Draw a straight line segment. In the field of computer graphics, in order to display straight line segments on raster display devices, a variety of line scan conversion techniques have been developed, among which the Bresenham algorithm is the most widely used. The principle of the algorithm is shown in Figure 2. The known pixel point (x _p , z _p ) is a point on the line, and the next pixel point has two options: P ₁ (x _p +1, z _p ) or P ₂ (x _p +1, z _p +1), the algorithm uses the sign of the error term ε to determine whether the next pixel is the right point or the bottom right point:

If the equation of the straight line is z=kx+b, then z _p+1 =z _p +k(x _p+1 -x _p )=z _p +k, where k is the slope of the straight line. If |k|≤1, it can be seen from the coordinates of the two candidate points of the next pixel that the abscissa of the next pixel is x _p +1, and the ordinate is either still z _p or incremented by 1 to z _p + 1, whether to increase by 1 depends on the value of the error term ε. The initial value of ε is 0, when the abscissa increases by 1, the corresponding incremental value of ε is the slope k of the straight line, that is, ε=ε+k.

If |k|>1, it should be considered whether x _p increases by 1, and the incremental value of ε is 1/k. When ε≥0.5, the intersection point of the straight line and the grid is closer to P ₂ , and the ordinate is increased by 1. At the same time, P ₂ is used as the new reference point for the next calculation, and the ε value is correspondingly subtracted by 1; when ε<0.5, the straight line The intersection point with the grid is closer to P ₁ , the ordinate remains unchanged, and the reference point for the next calculation also remains unchanged.

位于图像内的直线段。For the convenience of calculation, set δ=ε-0.5, the initial value of δ is -0.5, and the increment is k. When δ≥0, the next pixel is (x _p +1, z _p +1), and the value of δ decreases 1; when δ<0, the next pixel takes (x _p +1, z _p ), and the value of δ remains unchanged. Then for the actual application scenario shown in Fig. 1(d) in the present invention, the initial point (x _p , z _p )=(x _R , z _R ) on the straight line in the Bresenham algorithm, k is the refraction vector The slope of , scan converts pixels on the straight line until reaching the boundary of the image. Therefore, the Bresenham algorithm can be used to conveniently scan-convert the refraction vector

Line segments that lie within the image.

This method of using the refraction formula to calculate the refraction vector, and using the linear scan conversion technology to obtain the pixel points on the refraction path to indirectly realize the drawing of the trajectory curve is a positive process of solving the refraction path. This method takes the refraction point and delay as known conditions, so it does not need to iteratively calculate the refraction point, and does not need to use the root mean square calculation to find the length of the acoustic wave propagation path, which saves a lot of computational overhead. In addition, this method replaces the calculation of the effective length L of the synthetic aperture in the original SAFT technology with the calculation of the half-power beam angle β _0.5 , and it is only related to the first layer of media, and the other layers do not need to calculate this value, so the imaging results of this method It will be more accurate than the original SAFT combined with ray tracing method.

Contents of the invention

The purpose of the present invention is to propose an ultrasonic imaging method to realize rapid and accurate imaging of layered objects.

The present invention is characterized in that it can image a regularly layered object containing a horizontal interface or an interface parallel to each other, and can also image an irregularly layered object containing a non-horizontal and non-parallel interface, and the imaging high speed.

The present invention is characterized in that, based on the synthetic aperture focusing (SAFT) ultrasonic imaging technology, combined with the law of refraction, the refraction vector calculation formula is used to calculate the refraction path in a forward direction, and the Bresenham linear scanning conversion algorithm is used to scan and convert each pixel point on the refraction path. The original SAFT method combined with ray tracing is implemented as a process of drawing trajectory curves on the image using sampled data, which avoids the iterative calculation of refraction points and the root mean square calculation of sound wave propagation distance, saves a lot of calculation operations, and uses half The power beam angle replaces the effective length of the synthetic aperture, which avoids the error of the imaging result caused by the inaccuracy of the calculation of the effective length of the synthetic aperture in the layered medium.

The present invention is characterized in that it contains the following steps in sequence:

Step (1) Construct a computer, an ultrasonic transducer, a set of positioning controller and an analog-to-digital converter based on synthetic aperture focusing technology and the law of refraction for the depth and level of layered objects A system for noninvasive ultrasound imaging on a longitudinal section formed in two directions, wherein:

The ultrasonic transducer is provided with: a pulse signal input end connected to the output end of the positioning controller, the input end of the positioning controller is connected to the corresponding positioning control signal output end of the computer, and the ultrasonic transducer The energy device is also provided with: an echo signal output end connected to the input end of the analog-to-digital converter, the output end of the analog-to-digital converter is connected to the echo sampling signal input end of the computer, and the ultrasonic transducer The transducer is controlled by the positioning controller, and moves at a fixed rate of 1 step/ms on the surface of the measured object. The positioning controller is a transmission device that controls the moving position of the ultrasonic transducer, and its parameters are determined by the computer input,

The horizontal length of the measured object along the X-axis is X _length , which is equally divided into X _length /Δx intervals, and Δx is the interval length, which is also the distance from the coordinate point (0,0) to the end point of the ultrasonic transducer along the X-axis. (X _length , 0) is the step size of each movement, the point reached by the ultrasonic transducer each time is called the detection point, there are M, M=1+X _length /Δx, serial number m=0, 1,..., M-1, the positioning controller generates a TTL transistor-transistor logic level pulse at each detection point, triggering the ultrasonic transducer to move to the depth direction Z perpendicular to the X-axis of the measured object An excitation pulse is transmitted, and then the ultrasonic transducer is switched to a receiving mode and starts timing to receive the echo signal reflected from the measured object, and the analog-to-digital converter performs an The echo signal is sampled for N times and stored in the computer, the sampling number n=0,1,...,N-1, the sampling frequency is f _s , the value of f _s is preset by the analog-to-digital converter, denote s _m (n ) is the sampling value obtained by the nth sampling of the ultrasonic transducer at the m detection point, and the sampling moment of s _m (n) is t=n/f _s ;

Step (2): The computer reads the sampling value at the detection point m=0 in sequence from n=0, and then repeats the process to read m=1 successively,..., the sampling at each detection point of M-1 value;

Step (3): Take v=v ₁ , where v ₁ is the propagation velocity of ultrasound in the first medium of the object to be measured, and use a synthetic aperture focused ultrasound imaging software package to generate z ₀ =0 to Z _depth -1 in the depth direction The longitudinal section image of the interval, Z _depth is the preset length of the generated image, that is, the depth value of the generated image represented by the number of pixels in the vertical direction;

Step (4): Using the z ₀ to Z _depth -1 interval image blocks obtained in step (3) as input, use the Canny operator edge extraction software package to extract the boundary line c ₁ between the first layer medium and the second layer medium ( x,z);

Step (5): according to the following steps, correct the image on the vertical section below the dividing line c ₁ (x, z) to the Z _depth -1 interval, to eliminate the heterogeneity between the first layer and other layers of media. Resulting error:

Step (5.1): Take m=0, record the mth detection point as U(x _u ,0), where x _u =mΔx/accuracy, accuracy is the image accuracy, that is, two adjacent pixels on the generated image Point spacing, follow the steps below to calculate the fan-shaped image corresponding to the mth detection point U(x _u ,0):

Step (5.1.1): calculate the half-power beam angle β _0.5 =0.84λ/D of the ultrasonic transducer, λ is the wavelength when the ultrasonic wave propagates in the measured object, and D is the diameter of the ultrasonic transducer, and calculate the The left and right boundary lines of the half-power beam angle and the left and right intersection points B _l (x _l , z l ) and B _r (x _r , _{z r} ) of the dividing line c ₁ (x, z) respectively, where 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;

其中，eta为所述被测物体中与分界线c₁(x，z)相邻的两层介质间的相对折射率，

是分界线上点R(x_R,z_R)处的单位法向矢量， $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})],$ “·”为向量点积运算，

为从所述探测点U(x_u,0)到折射点R(x_R,z_R)的归一化的入射向量；Step (5.1.2): Take the refraction point R(x _R , z _R )=B _l (x _l , z _l ) on the dividing line c ₁ (x, z), and calculate the normalized refraction vector

Wherein, eta is the relative refractive index between the two layers of media adjacent to the boundary line c ₁ (x, z) in the measured object,

is the unit normal vector at point R(x _R , z _R ) on the dividing line, $a=1.0-\mathrm{eta}\&times;\mathrm{eta}\&times;[1.0-(\stackrel{\&Right\; Arrow;}{\mathrm{RU}}\&CenterDot;\stackrel{\&Right\; Arrow;}{e})\&times;(\stackrel{\&Right\; Arrow;}{\mathrm{RU}}\&Center\; Dot;\stackrel{\&Right\; Arrow;}{e})],$ "·" is the vector dot product operation,

is the normalized incident vector from the detection point U(x _u ,0) to the refraction point R(x _R ,z _R );

的斜率k，取直线上的初始点为R(x_R,z_R)，利用Bresenham直线扫描转换软件包，从折射点R开始，计算折射路径上的所有像素点直到到达图像的左边界x=0或者右边界x=(M-1)Δx/accuracy或者下边界z=Z_depth-1为止，记当前折射点R为R_f，即R_f(x_f,z_f)=R(x_R,z_R)，并记折射路径为R_fE_f，E_f为R_fE_f在图像内的终点，即折射路径R_fE_f与图像边界的交点；Step (5.1.3): Calculate the refraction vector

The slope k of , take the initial point on the straight line as R(x _R , z _R ), use the Bresenham linear scan conversion software package, start from the refraction point R, and calculate all the pixels on the refraction path until reaching the left boundary of the image x= 0 or the right boundary x=(M-1)Δx/accuracy or the lower boundary z=Z _depth -1, record the current refraction point R as R _f , that is, R _f (x _f ,z _f )=R(x _R , z _R ), and record the refraction path as R _f E _f , E _f is the end point of R _f E _f in the image, that is, the intersection point of the refraction path R _f E _f and the image boundary;

Step (5.1.4): Take x _g =x _f +1, find the next pixel point R _g (x g , z g ) at the position of the current refraction point R _f on the dividing line c ₁ ( _x , z ₎ , point The coordinates of R _g satisfy c ₁ (x _g , z _g )=0, take R _g as the new refraction point, that is, R(x _R ,z _R )=R _g (x _g ,z _g ), execute step (5.1 .2) and step (5.1.3), scan and convert the new refraction path R _g E _g , record the end point as E _g ;

Step (5.1.5): Calculate the distance Δd=|E between the end point E _g of the refraction path R _g E _g of the current refraction point R _g and the end point E _f of the refraction path R _f E _f of the previous refraction point R _{f f} E _g |, divide the curve segment between the current refraction point R _g and the previous refraction point R _f on the dividing line c ₁ (x, z) into Δd parts, that is, insert Δd-1 points, and insert the point number Denoted as τ=1,2,...,Δd-1, the following steps (5.1.5.1) are cyclically executed for the value of τ with 1 as the step size and τ=1 as the initial value until τ=Δd:

Step (5.1.5.1): Calculate the coordinates of the τth insertion point R _τ by interpolation, x _τ = x _f +τ/Δd, z _τ = z _f +(z _g -z _f )τ/Δd, take the refraction point R (x _R , z _R )=R _τ (x _τ , z _τ ), execute step (5.1.2) and step (5.1.3), scan and convert the refraction path starting from point R _τ ;

Step (5.1.6): take R _f =R _g , E _f =E _g , execute step (5.1.4) and step (5.1.5);

Step (5.1.7): Repeat step (5.1.6) until x _g = x _r +1, that is, after processing the boundary line c ₁ (x, z) between B _l (x _l , z _l ) and All refraction points and refraction paths between B _r (x _r , z _r ), get the fan-shaped image corresponding to the mth detection point U(x _u ,0);

Step (5.2): Take m=1,...,M-1 in sequence, repeat step (5.1), and generate an image from the section below the dividing line c ₁ (x,z) to Z _depth -1 on the longitudinal section ;

Step (6): With the boundary line c ₁ (x, z) obtained in step (5) below to the image block of the Z _depth -1 interval as input, use the Canny operator edge extraction software package to extract the boundary line c The boundary c ₂ (x, z) between the second layer medium and the third layer medium between ₁ (x, z) and Z _depth -1, according to the method described in step (5), correct the The image of the interval between the boundary line c ₂ (x, z) and Z _depth -1, to eliminate the error caused by the heterogeneity between the second layer and the layers below it;

Step (7): According to the method described in step (6), process the remaining boundary lines until all the boundary lines in the longitudinal section are processed, and the generated width is (M-1)Δx/accuracy+1 pixels , an image of the longitudinal section whose length is Z _depth pixels.

Compared with the prior art, the invention has the advantages that it can image the layered object with regular interface or irregular interface, and the imaging speed is fast and the imaging is more accurate. For example: for the measured object shown in Figure 3(a), if the diameter of the ultrasonic transducer is 0.5mm, the moving step of the ultrasonic transducer is 0.7mm, the center frequency of the ultrasonic waves emitted by the ultrasonic transducer is 5MHz, and the sampling frequency 100MHz, the imaging accuracy is taken as 0.05mm, on the experimental machine of Intel Core Duo2.66GHz CPU, 2.0GBRAM, it only takes 16s to generate Figure 3(d) by using the method of the present invention, which is about the same as the original SAFT combined with the ray tracing method 1/112 of the time, 1/4 of the phase shift technology. In Figure 3(d), the maximum error of the upper surface is 0.5mm (as shown in Figure 4(a)), and the maximum error of the lower surface is 0.9mm (as shown in Figure 4(b)), and the imaging accuracy is high. Moreover, comparing with FIG. 3( b ) and FIG. 3( c ), it can be found that the method of the present invention has a significantly better imaging effect on the lower surface than the former two methods.

Description of drawings

Figure 1 is the working model and schematic diagram of SAFT ultrasonic imaging technology: 1(a) is the schematic diagram of the reverse calculation process of SAFT technology when the measured object is a single medium; 1(b) is the measured object is a layered medium In the case of , the schematic diagram of the reverse calculation process of SAFT technology; 1(c) is the forward explanation diagram of SAFT technology when the measured object is a layered medium; 1(d) is the measured object is a layered medium In the case of , the efficacy diagram of all the data sampled by the ultrasonic transducer at the scanning position u _i in the original SAFT technique.

Fig. 2 is a schematic diagram of pixels and error terms involved in each iteration of the Bresenham linear scan conversion algorithm.

Fig. 3 is a comparison diagram of each ultrasonic imaging method under the same experimental environment: 3 (a) is a longitudinal section view of an irregularly layered object; 3 (b) is an image generated by SAFT combined with ray tracing technology, and the imaging time 30 minutes, the imaging effect on the lower surface is not good; 3(c) is the image generated by phase shift ultrasonic imaging technology, the imaging time is 62s, and there are serious errors in the imaging of the lower surface; 3(d) is the image generated by this imaging method For the image, the imaging time is 16s, and the imaging effect on the lower surface is the best. The white curves in the imaging of the upper and lower surfaces in the figure are the upper and lower dividing lines obtained by the experiment, respectively.

Fig. 4 is an analysis diagram of the experimental data of Fig. 3(d): 4(a) is the error map of the first layer boundary line, the solid curve is the standard value of the first layer boundary line of the measured object, and the dotted point curve is the imaging method The obtained first layer boundary line, the dashed line curve is the error curve of the first layer boundary line obtained by this imaging method; 4(b) is the second layer boundary line error map, and the solid curve is the second layer of the measured object The standard value of the boundary line, the dashed point curve is the second layer boundary line obtained by this imaging method, and the dashed line curve is the error curve of the second layer boundary line obtained by this imaging method.

FIG. 5 is a schematic flow chart of the ultrasound imaging system.

FIG. 6 is a structural diagram of the ultrasonic imaging hardware system.

Fig. 7 is a schematic diagram of the operation of the ultrasonic transducer.

Fig. 8 is a flowchart of an algorithm for calculating synthetic aperture focused ultrasound imaging of a layered object.

Detailed ways

The specific implementation process of the present invention includes three parts (as shown in Figure 5): ultrasonic data acquisition, imaging calculation and image display. The hardware platform system structure diagram is shown in Figure 6. The ultrasonic imaging system consists of a computer, an ultrasonic transducer, a positioning controller and an analog-to-digital converter. The pulse signal input terminal of the ultrasonic transducer is connected to the positioning control The output end of the controller is connected, and the input end of the positioning controller is connected with the output end of the positioning control signal of the computer. The echo signal output end of the ultrasonic transducer is connected with the input end of the analog-digital converter, and the output end of the analog-digital converter is connected with the echo sampling signal input end of the computer.

The method of ultrasonic data acquisition can be selected: contact type, that is, the ultrasonic transducer directly contacts the surface of the object to be measured; or liquid immersion type, that is, the object to be measured is placed in the liquid, and the ultrasonic transducer measures on the liquid surface. Depending on the detection method, the ultrasonic transducer can use a contact probe or a liquid immersion probe. The system uses a single transmitting/receiving ultrasonic transducer, and the ultrasonic transducer moves with a uniform step size Δx along the X direction (as shown in Figure 7) at a fixed rate of about 1 step/ms through the positioning controller on the surface of the measured object or in the liquid . The controller generates a TTL (transistor-transistor logic level) pulse at every moment when the target position is stable, which is used to trigger the ultrasonic transducer to emit an excitation pulse to the depth direction Z perpendicular to the X direction of the measured object, The ultrasonic transducer then switches to receive mode and starts timing, receiving echoes reflected from the object under test. The locations where the ultrasonic transducer emits pulses and receives echoes are detection points. The echo signal received by the ultrasonic transducer is collected by an analog-to-digital converter and stored in a memory. The ultrasonic transducer moving step Δx needs to be comprehensively determined according to the actual size of the object to be measured and the imaging accuracy requirements. The smaller the value, the more accurate the generated image, but the longer the calculation time.

Imaging calculation is to use the sampling data of the measured object at each detection point on a longitudinal section as computer input, and then calculate the longitudinal section image of the measured object according to the aforementioned imaging steps. The overall flow chart is shown in Figure 8.

In the specific implementation program, the ultrasonic transducer movement step Δx in the aforementioned imaging step (1) can be converted into the number of pixels, and only need to execute Δx←Δx/accuracy, then the abscissa can be directly calculated in step (5.1) Value x _u =mΔx, avoiding the need to divide by the image precision for each m to be converted into pixel coordinates on the image, saving division operation overhead.

Since the exact information of the boundary line in the layered object is not known in advance, step (3) generates an image using a single medium ultrasonic imaging algorithm to extract the boundary line between the first layer and the second layer medium. Since the speed of sound in the first layer medium is used in the algorithm, the image formed on the first layer medium is accurate, and the boundary line between the first layer and the second layer medium is also a part of the first layer medium. Therefore, the boundary line extracted in step (4) is also accurate. In fact, the ultrasonic imaging algorithm used in step (3) can be either the original SAFT technique, or any technique capable of imaging a single medium, such as the phase shift technique.

In the imaging step (4), it is necessary to use the Canny operator edge extraction algorithm to extract the boundary between the media. The calculation process of the Canny edge extraction operator is divided into four steps:

Step (a): Image smoothing. The original image is smoothed and filtered with a two-dimensional Gaussian function to obtain an image I(x, z);

Step (b): Calculate the gradient G and direction F of each pixel of the image I(x,z). A 2×2 template is used for the X direction and

A first-order approximation of the partial differential in the Z direction, namely

$\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}\left(\begin{array}{cc}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}\\ <span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}\end{array}\right)<span\; class="notranslate">\; ,</span>$ $\mathrm{<span\; class="notranslate">\; q</span>}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}\left(\begin{array}{cc}\mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}\\ <span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}& <span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}\end{array}\right)$

Then the gradient size G and direction F are

$\mathrm{<span\; class="notranslate">\; G</span>}<span\; class="notranslate">\; =</span>\sqrt{\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; \&times;</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; q</span>}<span\; class="notranslate">\; \&times;</span>\mathrm{<span\; class="notranslate">\; q</span>}}<span\; class="notranslate">\; ,</span>$ $\mathrm{<span\; class="notranslate">\; f</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; arctan</span>}<span\; class="notranslate">\; (</span>\frac{\mathrm{<span\; class="notranslate">\; q</span>}}{\mathrm{<span\; class="notranslate">\; p</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ;</span>$

Step (c): Non-maximum suppression of the gradient image. For each pixel point I(x, z), if its gradient value G(x, z) is smaller than the gradient magnitude of two adjacent points along its gradient direction F(x, z), it means that the point is not an edge point, set the gray level of I(x, z) to 0;

Step (d): double thresholding. Set the low threshold Low and high threshold High required by the dual threshold method to detect and connect edges, and perform dual thresholding processing on the gradient image. If the gradient magnitude is greater than the high threshold High, it is an edge; if the gradient magnitude is smaller than the low threshold Low, it is not an edge; if the gradient magnitude is between the two, it is judged whether there is an edge greater than the high threshold High among the eight adjacent pixels of the pixel. If there is an edge pixel, it is an edge pixel, otherwise it is not.

In the entire imaging calculation, the half-power beam angle β _0.5 in the imaging step (5.1.1) needs to be used M times, but this value is only related to the intrinsic parameters of the ultrasonic transducer, and the present invention only uses one ultrasonic transducer for detection , so this value is a fixed value, in order to avoid repeated calculations, you can pre-calculate β _0.5 ←tg(0.5×0.84λ/D) in the program, save it as a global variable, and then call this directly in step (5.1.1) variables to calculate x _l = x _u - z _l × β _0.5 and x _r = x _u + z _r × β _0.5 .

In addition, in step (5.1.1), it is necessary to calculate the intersection points _{B l (} x _l ,z _l ) and B _r (x _r ,z _r ). In actual implementation, the boundary line c ₁ (x, z) extracted by the Canny operator is saved as a set of pixel points in increasing order of the x coordinate value. When calculating the intersection points B _l and B _r , it is only necessary to sequentially search for these pixel points. Yes, if the coordinates of a certain pixel point (x, z) on the boundary line c ₁ (x, z) satisfy x=x _u -z×β _0.5 , then it is B _l ; if it satisfies x=x _u +z×β _0.5 is B _r . Because in the whole imaging calculation, for each m, two left and right intersection points B _l and B _r need to be calculated, but the half-power beam angles corresponding to each m are parallel to each other (as shown in Figure 1(c)), then The left intersection point B _l corresponding to the current m must be on the right side of the left intersection point of the previous m, and the right intersection point B _r corresponding to the current m must also be on the right side of the right intersection point of the previous m, so each step (5.1.1) When calculating the left intersection in one round, you only need to start searching for the pixel points on the boundary line c ₁ (x, z) from B _l in the previous round. Correspondingly, when calculating the right intersection in each round, you only need to start from the previous round Start the search at B _r .

The specific algorithm of the imaging calculation process is divided into the Bresenham line scan conversion subroutine BresenhamLine() and the main program SAFT-Refraction() of the synthetic aperture focused ultrasound imaging algorithm for layered objects:

Image display is to display the two-dimensional image data obtained in the imaging calculation stage on the monitor, and grayscale images or color images can be displayed as required.

For the case where the object to be measured is a regular layered object, all the dividing lines c(x,z) extracted by the Canny operator edge extraction software package in the imaging steps (4) and (6) are straight lines, and the imaging calculation steps remain unchanged .

Claims (1)

1. The synthetic aperture focused ultrasonic imaging method of layered object, it is characterized in that, contains following steps successively: Step (1): Construct a computer, an ultrasonic transducer, a positioning controller and an analog-to-digital converter based on the synthetic aperture focusing technology and the law of refraction for the depth and depth of the layered object A system for non-invasive ultrasonic imaging on longitudinal sections formed in two horizontal directions, in which: The ultrasonic transducer is provided with: a pulse signal input end connected to the output end of the positioning controller, the input end of the positioning controller is connected to the corresponding positioning control signal output end of the computer, and the ultrasonic transducer The energy device is also provided with: an echo signal output end connected to the input end of the analog-to-digital converter, the output end of the analog-to-digital converter is connected to the echo sampling signal input end of the computer, and the ultrasonic transducer The transducer is controlled by the positioning controller, and moves at a fixed rate of 1 step/ms on the surface of the measured object. The positioning controller is a transmission device that controls the moving position of the ultrasonic transducer, and its parameters are determined by the computer input, The horizontal length of the measured object along the X-axis is X _length , which is equally divided into X _length /Δx intervals, and Δx is the interval length, which is also the distance from the coordinate point (0,0) to the end point of the ultrasonic transducer along the X-axis. (X _length , 0) is the step size of each movement, the point reached by the ultrasonic transducer each time is called the detection point, there are M, M=1+X _length /Δx, serial number m=0, 1,..., M-1, the positioning controller generates a TTL transistor-transistor logic level pulse at each detection point, triggering the ultrasonic transducer to move to the depth direction Z perpendicular to the X-axis of the measured object An excitation pulse is transmitted, and then the ultrasonic transducer is switched to a receiving mode and starts timing to receive the echo signal reflected from the measured object, and the analog-to-digital converter performs an The echo signal is sampled for N times and stored in the computer, the sampling number n=0,1,…,N-1, the sampling frequency is f _s , the value of f _s is preset by the analog-to-digital converter, denote s _m (n ) is the sampling value obtained by the nth sampling of the ultrasonic transducer at the m detection point, and the sampling moment of s _m (n) is t=n/f _s ; Step (2): The computer reads the sampling value at the detection point m=0 in sequence from n=0, and then repeats the process to read m=1 successively,..., the sampling at each detection point of M-1 value; Step (3): Take v=v ₁ , where v ₁ is the propagation velocity of ultrasound in the first layer medium of the object to be measured, and use a synthetic aperture focused ultrasound imaging software package to generate z ₀ =0 to Z _depth -1 in the depth direction The longitudinal section image of the interval, Z _depth is the preset length of the generated image, that is, the depth value of the generated image represented by the number of pixels in the vertical direction; Step (4): Using the z ₀ to Z _depth -1 interval image blocks obtained in step (3) as input, use the Canny operator edge extraction software package to extract the boundary line c ₁ between the first layer medium and the second layer medium ( x,z); Step (5): according to the following steps, correct the image on the vertical section below the dividing line c ₁ (x, z) to the Z _depth -1 interval, to eliminate the heterogeneity between the first layer and other layers of media. Resulting error: Step (5.1): Take m=0, record the mth detection point as U(x _u ,0), where x _u =mΔx/accuracy, accuracy is the image accuracy, that is, two adjacent pixels on the generated image Point spacing, follow the steps below to calculate the fan-shaped image corresponding to the mth detection point U(x _u ,0): Step (5.1.1): calculate the half-power beam angle β _0.5 =0.84λ/D of the ultrasonic transducer, λ is the wavelength when the ultrasonic wave propagates in the measured object, and D is the diameter of the ultrasonic transducer, and calculate the The left and right boundary lines of the half-power beam angle and the left and right intersection points B _l (x _l , z l ) and B _r (x _r , _{z r} ) of the dividing line c ₁ (x, z) respectively, where 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): Take the refraction point R(x _R , z _R )=B _l (x _l , z _l ) on the dividing line c ₁ (x, z), and calculate the normalized refraction vector Wherein, eta is the relative refractive index between the two layers of media adjacent to the boundary line c ₁ (x, z) in the measured object,

is the unit normal vector at point R(x _R , z _R ) on the dividing line, $a=1.0-\mathrm{eta}\&times;\mathrm{eta}\&times;[1.0-(\stackrel{\&Right\; Arrow;}{\mathrm{RU}}\&Center\; Dot;\stackrel{\&Right\; Arrow;}{e})\&times;(\stackrel{\&Right\; Arrow;}{\mathrm{RU}}\&Center\; Dot;\stackrel{\&Right\; Arrow;}{e})],$ "·" is the vector dot product operation,

is the normalized incident vector from the detection point U(x _u ,0) to the refraction point R(x _R ,z _R ); Step (5.1.3): Calculate the refraction vector

The slope k of , take the initial point on the straight line as R(x _R , z _R ), use the Bresenham linear scan conversion software package, start from the refraction point R, and calculate all the pixels on the refraction path until reaching the left boundary of the image x= 0 or the right boundary x=(M-1)Δx/accuracy or the lower boundary z=Z _depth -1, record the current refraction point R as R _f , that is, R _f (x _f ,z _f )=R(x _R , z _R ), and record the refraction path as R _f E _f , E _f is the end point of R _f E _f in the image, that is, the intersection point of the refraction path R _f E _f and the image boundary; Step (5.1.4): Take x _g =x _f +1, find the next pixel point R _g (x g , z g ) at the position of the current refraction point R _f on the dividing line c ₁ ( _x , z ₎ , point The coordinates of R _g satisfy c ₁ (x _g , z _g )=0, take R _g as the new refraction point, that is, R(x _R ,z _R )=R _g (x _g ,z _g ), execute step (5.1 .2) and step (5.1.3), scan and convert the new refraction path R _g E _g , record the end point as E _g ; Step (5.1.5): Calculate the distance Δd=|E between the end point E _g of the refraction path R _g E _g of the current refraction point R _g and the end point E _f of the refraction path R _f E _f of the previous refraction point R _{f f} E _g |, divide the curve segment between the current refraction point R _g and the previous refraction point R _f on the dividing line c ₁ (x, z) into Δd parts, that is, insert Δd-1 points, and insert the point number Denoted as τ=1,2,...,Δd-1, the following steps (5.1.5.1) are cyclically executed for the value of τ with 1 as the step size and τ=1 as the initial value until τ=Δd: Step (5.1.5.1): Calculate the coordinates of the τth insertion point R _τ by interpolation, x _τ = x _f +τ/Δd, z _τ = z _f +(z _g -z _f )τ/Δd, take the refraction point R (x _R , z _R )=R _τ (x _τ , z _τ ), execute step (5.1.2) and step (5.1.3), scan and convert the refraction path starting from point R _τ ; Step (5.1.6): take R _f =R _g , E _f =E _g , execute step (5.1.4) and step (5.1.5); Step (5.1.7): Repeat step (5.1.6) until x _g = x _r +1, that is, after processing the boundary line c ₁ (x, z) between B _l (x _l , z _l ) and All refraction points and refraction paths between B _r (x _r , z _r ), get the fan-shaped image corresponding to the mth detection point U(x _u ,0); Step (5.2): Take m=1,...,M-1 in sequence, repeat step (5.1), and generate an image from the section below the dividing line c ₁ (x,z) to Z _depth -1 on the longitudinal section ; Step (6): With the boundary line c ₁ (x, z) obtained in step (5) below to the image block of the Z _depth -1 interval as input, use the Canny operator edge extraction software package to extract the boundary line c The boundary c ₂ (x, z) between the second layer medium and the third layer medium between ₁ (x, z) and Z _depth -1, according to the method described in step (5), correct the The image of the interval between the boundary line c ₂ (x, z) and Z _depth -1, to eliminate the error caused by the heterogeneity between the second layer and the layers below it; Step (7): According to the method described in step (6), process the remaining boundary lines until all the boundary lines in the longitudinal section are processed, and the generated width is (M-1)Δx/accuracy+1 pixels , an image of the longitudinal section whose length is Z _depth pixels.

Images