CN103018333A - 分层物体的合成孔径聚焦超声成像方法 - Google Patents

Abstract

分层物体的合成孔径聚焦超声成像方法属于多层异质物体的内部结构快速成像技术领域，其特征在于，基于使用一个超声换能器在多层异质物体水平方向作等间距移动探测的条件下，基于合成孔径聚焦技术(SAFT)和折射定律，将原SAFT的实现过程改为以采样数据在图像上画轨迹曲线的过程，并结合折射定律正向计算折射向量，利用直线扫描转换技术得到折射路径上的各像素点，从而提高成像速度。本发明既适用于仅在深度方向存在异种介质的规则分层物体的超声成像，也适用于在深度方向和水平方向上均存在异种介质的非规则分层物体的超声成像，且成像速度快，成像精度高。

Description

分层物体的合成孔径聚焦超声成像方法

技术领域

本发明涉及超声波无损检测技术、合成孔径聚焦成像技术、折射定律和直线扫描转换技术，实现对包含规则分界面或非规则分界面的分层物体的表面形状及内部结构快速成像。

背景技术

目前，在超声波无损检测领域，对于含有分界面的分层物体的超声成像（超声检测方式包含接触式和液浸式两种模式，若将耦合剂层视为被测物体的一部分，则液浸式超声检测也可视为对分层物体的接触式超声检测，因此本发明中将两者统一称为分层物体的超声检测及成像），主要有两种方案：一是采用合成孔径聚焦(SAFT)超声成像技术与射线跟踪(Ray Tracing)技术相结合的方法，二是采用基于相位迁移(Phase Shift Migration)技术的相位迁移超声成像方法。

合成孔径聚焦技术(SAFT)源自于合成孔径雷达技术(SAR)，于20世纪70年代初被引入到超声成像领域，其具有不受近场区限制、方位分辨率高、分辨率只与超声换能器尺寸有关而与距离无关等优点。SAFT超声成像技术的基本原理是利用脉冲—回波(pulse-echo)测量机制，使用一个超声换能器沿着固定轨迹对被测物体进行有序的扫描，并采用延时叠加(DAS)方法对扫描得到的脉冲回波信号进行聚焦成像，达到利用单一的较小孔径的超声换能器来模拟大的孔径阵列的目的。SAFT超声成像工作模型如图1(a)所示，超声换能器沿着扫描方向（X向）在物体表面作等间距的移动，在每一个扫描位置向物体的深度方向（Z向）发射超声信号，对物体内部进行探测，同时超声换能器接收物体内部反射物反射回的回波信号并采样保存，最后对所有扫描位置处得到的采样数据进行DAS与多点动态聚焦处理并显示图像。

DAS与多点动态聚焦技术需要对目标成像区域不同深度上的聚焦点计算不同的时延曲线。如图1(a)所示，为了在目标反射物处(x，z)聚焦，SAFT技术将超声换能器在其合成孔径有效长度L内的每一个扫描位置处获得的回波信号进行延时叠加处理：设s_i(t)为超声换能器在u_i处接收到的回波信号，t为采样时刻，u_i处关于目标反射物(x,z)的延时为

$t_i=\frac{2r_i}{v}=\frac{2}{v}\sqrt{z^2+{(x-u_i)}^2},i=\mathrm{0,1},...,L-1.---\left(1\right)$

其中，v为超声在介质中的传播速度，r_i为(x，z)点距u_i的直线距离。合成孔径有效长度L内所有的延时构成一条延时曲线，该曲线为一段双曲线。L的计算公式为

L=0.84λz/D    (2)

λ为超声在介质中的波长，D为超声换能器的直径，则在(x,z)处的成像为

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

其中，ω_i为变迹函数。

从DAS的原理可知，SAFT算法需要根据超声换能器到反射物之间的相对距离以及波速来确定延时。而当被测物体分层时，如图1(b)所示，超声换能器发射的超声信号会在分界面处发生折射，声波的传播路径会发生改变，目标反射物(x,z)与u_i的距离r_i不再是两点间直线段的长度。此外，分层物体各层的介质也不一样，声波在各层的传播速度通常也不相等，因此，DAS中时延的计算不能再采用相对距离除以波速的简单计算方法，延时曲线也不再是规则的双曲线形状。为了求得延时，需要先获得声波在u_i与目标反射物(x，z)之间的传播路径，然后逐段求出各路径段的长度及传播时间(如图1(b)中r_i1、r_i2和t_i)。其重点和难点在于快速而准确地找出波束的传播路径，而射线跟踪技术正是实现此目的之最佳方法。所以，射线跟踪技术很自然的被结合到SAFT技术中，以实现对分层物体的超声成像。

射线跟踪的原理主要基于Snell定理或Fermat原理，通过迭代计算找出费时最短的声波传播路径。所以，如图1(b)，在SAFT结合射线跟踪的方法中，当计算u_i与目标反射物(x,z)之间的声波传播延时t_i时，需要沿着分界线c(x，z)=0对线上所有的点进行迭代，计算每个点所对应的两段传播路径r_i1、r_i2的长度及其传播时间，找出使得传播时间最短的点，该点即为折射点。

此方法的优点是射线跟踪技术对分层物体的分界面形状并没有特殊的要求，不仅能适用于在深度方向上存在异种介质而水平方向为同种介质，即层与层介质之间的分界面为水平或者互相平行的平面之类的规则分层物体的超声成像，也能适用于在深度方向和水平方向上均存在异种介质，即包含非规则分界面(如曲面)的非规则分层物体的超声成像。但其严重的缺点是射线跟踪算法包含迭代运算，时间复杂度高，此外，原SAFT技术中计算时延的公式(1)中包含一次均方根运算，计算开销已经很大，而当引入射线跟踪技术后，时延的计算则包含多次均方根运算，如图1(b)中，当计算u_i与目标反射物(x,z)之间的声波传播延时t_i时，需要对分界线上的每一个点，分别计算两次均方根以得到两段路径r_i1、r_i2的长度，其计算开销更大。所以该方法非常耗时，如图3(b)所示，用该方法对图3(a)所示的被测物体（4.4cm*4.9cm）成像就需要30分钟。另外，由于超声波的波长与介质有关，当被测物体的介质分层后，合成孔径有效长度L的计算公式(2)也不再适用，然而SAFT技术所成像的精确度与L的准确度紧密相关，如果L过大，成像噪声较多，信噪比低；如果L过小，所成像会丢失被测物体的细节，精度不高。但是对于复杂非规则分层物体，目前缺少准确计算L的公式，在SAFT结合射线跟踪的方法中很难找出准确的L值，导致该方法的成像效果并不好，例如图3(b)中，该方法对被测物体纵剖面所成的图像未能较好的重构出第二层分界面中的曲面部分。

相位迁移(Phase Shift Migration)超声成像方法是将反射地震学(Reflection Seismology)中的迁移技术(Migration Technique)引入到超声成像领域，而得到的一种频域下的超声成像方法。该方法将超声探测模型视为***物反射模型，假设待测物体中的反射物在t=0时刻同时***，每一个反射物的***强度正比于其反射率，整个场强用一组超声换能器来测量。其主要思想是根据从水平位置（即深度方向第一行）处观测到的声场外推以计算出深度方向其他位置处的声场。具体算法包含两个主要的步骤：第一步对时域数据进行二维傅里叶变换，得到二维频谱，

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

第二步是在深度方向的循环，先对上一次循环得到的二维频谱作相移，然后作二维傅里叶反变换并取t0，得到一行图像，

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

其中，

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

对于多层物体的成像，相位迁移超声成像方法无需求射线路径，成像速度可大幅提高，如生成图3(c)所示的图像仅需62s。但是，由于相位迁移技术中假设声波的波速恒定，使得该方法只能适用于深度方向上存在异种介质而水平方向须为同种介质的规则分层物体的超声成像，而对于包含非规则分界面的物体，在深度方向和水平方向均存在异种介质，而导致检测结果不够精确，如图3(c)所示，第二层分界面的成像结果已经出现严重的错误。

因此，在对含有非规则分界面的复杂分层物体的超声成像方面，目前还缺乏快速且准确有效的方法。

从SAFT结合射线跟踪技术的原理中可以看出，该方法最耗时之处在于寻找最优折射点的迭代运算步骤。从射线跟踪技术和DAS原理可知，该方法是一种逆向的计算过程，需要先确定图像上的一个像素点，然后找出该点对应的所有扫描位置点，再对每一个扫描位置点在分界线上迭代找出最优的折射点以得到图像点与该扫描位置点间的距离及延时，最后对延时曲线上的数据点进行累加得到该图像点的像素值。由于事先并不知道被测物体内目标反射物的具***置，就需要使用动态聚焦技术对图像内所有像素点都进行一遍该逆向计算过程，以生成整幅图像，从而使得在反射物对应的像素点处，超声信号得到一致的叠加，达到累加强度最大化，实现聚焦，而在其他像素点处，超声信号的叠加是混乱的，累加强度难以最大化。因此，在图像中，反射物所对应处累加得到的像素值要明显大于其他像素点。

如果将图像上的所有像素点的计算过程作整体考察，可以发现每个扫描位置处采样得到的每个数据不仅作用于反射物对应的像素点，也作用于无反射物的像素点。如图1(c)中，u_i处的采样数据s_i(t_i)不仅参与了(x,z)点的成像计算，也参与了曲线段上其他点的计算。

是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的圆弧；而当两层介质不同时，该轨迹的形状则与分界线的形状有关。

因此，从正向去综合理解整体计算过程，采样数据s_i(t_i)在整幅图像的成像计算中的功效就相当于以数据值ω_is_i(t_i)/r_i在图像中画了一段曲线

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

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

至

）在图像中的交点即为原逆向计算过程中的目标反射物点(x,z)。

若在整幅图像中仅考虑u_i处所有采样数据s_i(t)的功效，则对应于一幅以u_i为中心的不规则的扇形图I_ui（如图1(d)），即

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

N为超声换能器在每个扫描位置处的采样次数。原SAFT结合射线跟踪方法中的DAS计算公式(3)则可以重新理解为各扫描位置处所对应的扇形图的叠加，即：

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

其中，M为超声换能器扫描位置点的总数。

所以，原SAFT结合射线跟踪的方法可以通过在图像上以每一个扫描位置处的每一个采样数据画轨迹曲线来实现。为了准确计算轨迹曲线

需要求得该曲线的函数表达式。如图1(c)，假设点P(x_P,z_P)是轨迹曲线上的任意一点，超声波从超声换能器所在位置U(u_i,0)处发射，在分界线上的点R(x_R,z_R)处折射后到达P点所用的时间为t_i/2，即

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

其中，v_water为声波在水中的传播速度，v_object为声波在被测物体中的传播速度。同时，根据折射定律：

eta＝sinα_i/sinα_r＝v_water/v_object    (7)其中，α_i为入射角，α_r为折射角，eta为折射率。由式(6)和式(7)可以求得P点的坐标，即得轨迹曲线的函数，但为了计算更简便，可以使用计算机图形学中的折射公式来计算折射向量

$\stackrel{\→}{\mathrm{RP}}=\mathrm{eta}\×\stackrel{\→}{\mathrm{RU}}-[\mathrm{eta}\×(\stackrel{\→}{\mathrm{RU}}\·\stackrel{\→}{e})+\sqrt{a}]\×\stackrel{\→}{e}---\left(8\right)$

其中，是分界线上点R(x_R,z_R)处的法向矢量，

“·”表示向量点积运算。

由于R已知，由式(8)得到

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

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

后再画下一段曲线

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

和

上的点，在画曲线

时需要计算一次

而画曲线

时同样需要计算一次

然而

和

的归一化矢量相同，所以

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

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

在折射点R处沿着折射矢量处理折射线上位于图像内的所有的像素点实质就是在图像内从R点出发，以各像素点对应的ω_is_i(t_i)/r_i值，沿着折射矢量

画一条直线段。而在计算机图形学领域，为了在光栅显示设备上显示直线段，已开发了多种直线扫描转换技术，其中，使用最广泛的是Bresenham算法。该算法的原理如图2所示，已知像素点(x_p,z_p)为直线上一点，下一个像素点有两种可选择点P₁(x_p+1,z_p)或P₂(x_p+1,z_p+1)，算法使用误差项ε的符号决定下一个像素取右边点还是右下点：

设直线方程为z=kx+b，则有z_p+1=z_p+k(x_p+1-x_p)=z_p+k，其中，k为直线斜率。如果|k|≤1，从下一个像素的两个候选点的坐标可以看出，下一个像素的横坐标为x_p+1，而纵坐标要么仍为z_p，要么递增1为z_p+1，是否增1取决于误差项ε的值。ε的初值为0，当横坐标每增加1，ε值的相应递增量为直线的斜率k，即ε=ε+k。

如果|k|>1，则应考虑x_p否增1，ε值的递增量为1/k。当ε≥0.5时，直线与网格的交点离P₂较近，纵坐标增1，同时以P₂作为下一次计算的新基准点，ε值相应减去1；当ε<0.5时，直线与网格的交点离P₁较近，纵坐标不变，下一次计算的基准点也不变。

为方便计算，令δ＝ε-0.5，δ的初值为-0.5，增量为k，当δ≥0时，下一个像素取(x_p+1,z_p+1)，δ值减小1；当δ<0时，下一个像素取(x_p+1,z_p)，δ值不变。则对于本发明中图1(d)所示的实际应用场景，Bresenham算法中直线上的初始点(x_p,z_p)=(x_R,z_R)，k为折射矢量的斜率，扫描转换直线上的像素点直到达到图像的边界为止。所以，利用Bresenham算法能方便的扫描转换折射矢量

位于图像内的直线段。

这种利用折射公式计算折射向量，并利用直线扫描转换技术求折射路径上的像素点，以间接实现画轨迹曲线的方法是一种正向的求解折射路径的过程。此方法将折射点和延时作为已知条件从而无需迭代求折射点，也无需使用均方根运算求声波传播路径长度，节省了大量的计算开销。此外，此方法将原SAFT技术中的求合成孔径有效长度L替换为求半功率波束角β_0.5，并且只与第一层介质有关，其余层次都无需再计算该值，所以此方法的成像结果会比原SAFT结合射线跟踪的方法更准确。

发明内容

本发明的目的在于提出一种超声成像方法，实现对分层物体的快速、准确的成像。

本发明的特征在于，能够对含有水平分界面或者含有互相平行的分界面的规则分层物体成像，同时也能对含非水平且非互相平行的分界面的非规则分层物体成像，且成像速度快。

本发明的特征在于，基于合成孔径聚焦(SAFT)超声成像技术，结合折射定律，使用折射向量计算公式正向计算折射路径，并利用Bresenham直线扫描转换算法扫描转换折射路径上的各像素点，将原SAFT结合射线跟踪的方法实现为使用采样数据在图像上画轨迹曲线的过程，避免了求折射点的迭代运算以及求声波传播距离的均方根运算，节省了大量的计算操作，并且用半功率波束角替代了合成孔径有效长度，避免了合成孔径有效长度在分层介质中计算的不准确性所带来的成像结果误差。

本发明的特征在于，依次含有以下步骤：

步骤(1)构建一个由一台计算机、一个超声换能器、一套定位控制器和一个模数转换器组成的一个基于合成孔径聚焦技术和折射定律的用于对分层物体在深度和水平两个方向形成的纵断面上作无损伤超声成像的***，其中：

所述超声换能器设有：与所述定位控制器的输出端相连的脉冲信号输入端，所述定位控制器的输入端与所述计算机相应的定位控制信号输出端相连，所述超声换能器还设有：与所述模数转换器的输入端相连的回波信号输出端，所述模数转换器的输出端与所述计算机的回波采样信号输入端相连，所述超声换能器由所述定位控制器控制，在被测物体表面以1步长/ms的固定速率移动，所述定位控制器是控制所述超声换能器移动位置的传动装置，其参数由所述计算机输入，

被测物体沿X轴方向的水平长度为X_length，均分为X_length/Δx个区间，Δx为区间长度，也是所述超声换能器沿X轴从坐标点(0,0)起向终点(X_length,0)止每次移动的步长，所述超声换能器每次移动所达到的点称为探测点，共有M个，M＝1+X_length/Δx，序号m=0,1,…，M-1，所述定位控制器在每一个探测点处产生一个TTL晶体管-晶体管逻辑电平脉冲，触发所述超声换能器向被测物体的垂直于X轴的深度方向Z发射一个激励脉冲，随后超声换能器转为接收模式并开始计时，接收从被测物体反射的回波信号，所述模数转换器对所述超声换能器在探测点m处接收到的回波信号进行N次采样并存储到计算机中，采样序号n=0,1，…，N-1，采样频率为f_s，f_s的值为模数转换器预设，记s_m(n)为超声换能器在第m个探测点处的第n次采样得到的采样值，s_m(n)的采样时刻为t＝n/f_s；

步骤(2):所述计算机从n=0开始依序读取探测点m=0处的采样值，然后，重复该过程依次读取m=1，…,M-1各探测点处的采样值；

步骤(3):取v=v₁，v₁为超声在被测物体的第一层介质中的传播速度，使用合成孔径聚焦超声成像软件包生成深度方向上z₀=0至Z_depth-1区间的纵断面图像，Z_depth为预设的生成图像的长度，即生成的图像在纵向上用像素数表示的深度值；

步骤(4):以步骤(3)得到的z₀至Z_depth-1区间图像块作为输入量，使用Canny算子边缘提取软件包提取第一层介质与第二层介质的分界线c₁(x，z)；

步骤(5):按以下步骤修正所述纵断面上在分界线c₁(x,z)之下至Z_depth-1区间的图像，以消除第一层与其他各层介质间的不同质而导致的误差：

步骤(5.1):取m=0，记第m个探测点为U(x_u,0)，其中，x_u=mΔx/accuracy，accuracy为图像精度，即所生成的图像上相邻两个像素点的间距，依次按以下步骤计算第m个探测点U(x_u,0)处所对应的扇形图像：

步骤(5.1.1):计算超声换能器的半功率波束角β_0.5=0.84λ/D，λ为超声在被测物体中传播时的波长，D为超声换能器的直径，并计算所述半功率波束角的左、右边界线分别与分界线c₁(x，z)的左右两个交点B_l(x_l,z_l)、B_r(x_r，z_r)，其中，x_l=x_u-z_l×tg(0.5β_0.5)且满足c₁(x_l,z_l)=0，x_r=x_u+z_r×tg(0.5β_0.5)且满足c₁(x_r,z_r)=0；

步骤(5.1.2):在分界线c₁(x，z)上取折射点R(x_R,z_R)＝B_l(x_l,z_l)，计算归一化的折射向量

其中，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)的归一化的入射向量；

步骤(5.1.3):计算折射矢量

的斜率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与图像边界的交点；

步骤(5.1.4):取x_g=x_f+1，查找分界线c₁(x，z)上当前折射点R_f所在位置的下一个像素点R_g(x_g,z_g)，点R_g的坐标满足c₁(x_g,z_g)=0，取R_g为新的折射点，即R(x_R,z_R)＝R_g(x_g,z_g)，执行步骤(5.1.2)和步骤(5.1.3)，扫描转换新的折射路径R_gE_g，记终点为E_g；

步骤(5.1.5):计算当前折射点R_g的折射路径R_gE_g的终点E_g与上一个折射点R_f的折射路径R_fE_f的终点E_f之间的距离Δd＝|E_fE_g|，将分界线c₁(x，z)上当前折射点R_g与上一个折射点R_f之间的曲线段均分为Δd份，即***Δd-1个点，***点序号记为τ=1,2,…,Δd-1，对τ值以1为步长以τ=1为初始值循环执行下述步骤(5.1.5.1)直到τ=Δd止：

步骤(5.1.5.1):插值计算第τ个***点R_τ的坐标，x_τ=x_f+τ/Δd，z_τ＝z_f+(z_g-z_f)τ/Δd，取折射点R(x_R,z_R)＝R_τ(x_τ,z_τ)，执行步骤(5.1.2)和步骤(5.1.3)，扫描转换从点R_τ出发的折射路径;

步骤(5.1.6):取R_f＝R_g，E_f=E_g，执行步骤(5.1.4)和步骤(5.1.5)；

步骤(5.1.7):重复执行步骤(5.1.6)直到x_g=x_r+1止，即处理完分界线c₁(x，z)上介于B_l(x_l,z_l)和B_r(x_r,z_r)之间所有的折射点和折射路径，得到第m个探测点U(x_u,0)处所对应的扇形图像；

步骤(5.2):依次取m=1,…,M-1，重复执行步骤(5.1)，生成所述纵断面上在分界线c₁(x,z)之下至Z_depth-1区间的图像；

步骤(6):以步骤(5)得到的分界线c₁(x,z)之下至Z_depth-1区间的图像块作为输入量，使用Canny算子边缘提取软件包提取所述分界线c₁(x,z)至Z_depth-1之间的第二层介质与第三层介质的分界线c₂(x,z)，按照步骤(5)所述的方法修正所述纵断面上在分界线c₂(x,z)之下至Z_depth-1区间的图像，以消除第二层与其下各层介质间的不同质而导致的误差；

步骤(7):按步骤(6)所述的方法，处理剩余各条分界线，直到所述纵断面中所有分界线处理完止，生成宽度为(M-1)Δx/accuracy+1个像素、长度为Z_depth个像素的所述纵断面的图像。

本发明与现有的技术相比，优点在于能对含有规则分界面或含有非规则分界面的分层物体进行成像，成像速度快且成像更准确。例如：对于图3(a)所示的被测物体，若超声换能器直径为0.5mm，超声换能器移动步长为0.7mm，超声换能器发射的超声波中心频率为5MHz，采样频率100MHz，成像精度accuracy取为0.05mm，在Intel Core Duo2.66GHz CPU、2.0GBRAM的实验机器上利用本发明的方法生成图3(d)只需16s，约为相同情况下原SAFT结合射线跟踪方法所用时间的1/112、相位迁移技术的1/4。图3(d)中上表面的最大误差为0.5mm(如图4(a))，下表面的最大误差为0.9mm(如图4(b))，成像精度较高。而且与图3(b)和图3(c)对比可以发现，本发明的方法对下表面所成像的效果明显优于前两种方法。

附图说明

图1是SAFT超声成像技术的工作模型及原理图：1(a)是被测物体为单一介质的情况下，SAFT技术的逆向计算过程原理图；1(b)是被测物体为分层介质的情况下，SAFT技术的逆向计算过程原理图；1(c)是被测物体为分层介质的情况下，SAFT技术的正向解释说明图；1(d)是被测物体为分层介质的情况下，超声换能器在扫描位置u_i处采样得到的所有数据在原SAFT技术中的功效图。

图2是Bresenham直线扫描转换算法每步迭代涉及的像素和误差项示意图。

图3是各超声成像方法在相同实验环境下的对比图:3(a)为非规则分层的被测物体的纵剖面图；3(b)为SAFT结合射线跟踪技术生成的图像，成像时间30分钟，对下表面所成像的效果欠佳；3(c)为相位迁移超声成像技术生成的图像，成像时间62s，对下表面所成像出现严重错误；3(d)为本成像方法生成的图像，成像时间16s，对下表面所成像的效果最好，图中对上、下表面所成像中的白色曲线分别为实验获得的上、下分界线。

图4是对图3(d)的实验数据分析图：4(a)为第一层分界线误差图，实曲线为被测物体第一层分界线的标准值，虚点曲线为本成像方法获得的第一层分界线，短划线曲线为本成像方法所获得的第一层分界线的误差曲线；4(b)为第二层分界线误差图，实曲线为被测物体第二层分界线的标准值，虚点曲线为本成像方法获得的第二层分界线，短划线曲线为本成像方法所获得的第二层分界线的误差曲线。

图5是本超声成像***流程示意图。

图6是本超声成像硬件***结构图。

图7是超声换能器工作示意图。

图8是分层物体的合成孔径聚焦超声成像计算的算法流程图。

具体实施方式

本发明的具体实施过程包含三部分（如图5）：超声数据获取、成像计算和图像显示。硬件平台***结构图如图6所示，超声成像***由一台计算机、一个超声换能器、一套定位控制器和一个模数转换器组成，超声换能器的脉冲信号输入端与定位控制器的输出端相连，定位控制器的输入端与计算机的定位控制信号输出端相连。超声换能器的回波信号输出端与模数转换器的输入端相连，模数转换器的输出端与计算机的回波采样信号输入端相连。

超声数据获取的方式可选择：接触式，即超声换能器直接接触待测物体的表面；或者液浸式，即待测物体置于液体中，超声换能器在液体表面进行测量。根据检测方式的不同，超声换能器可以使用接触式探头或者液浸式探头。***使用单个发射/接收超声换能器，超声换能器通过定位控制器在被测物体表面或者液体中以均匀的步长Δx沿X方向（如图7）以约1步/ms固定速率移动。控制器在每一个目标位置稳定的瞬间产生一个TTL（晶体管-晶体管逻辑电平）脉冲，该脉冲用来触发超声换能器向被测物体与X方向相垂直的深度方向Z发射一个激励脉冲，随后超声换能器转为接收模式并开始计时，接收从被测物体反射的回波。超声换能器发射脉冲及接收回波的各个位置处为探测点。超声换能器接收到的回波信号由模数转换器采集并存储在存储器中。超声换能器移动步长Δx需根据待测物体的实际大小和成像精度要求来综合确定，其值越小，生成的图像越精确，但计算时间也越长。

成像计算就是以被测物体在一个纵断面上各探测点处的采样数据作为计算机输入，然后按前述的成像步骤计算被测物体的纵断面图像，总体流程图参见图8。

在具体实现程序中，前述成像步骤(1)中的超声换能器移动步长Δx可以转化为像素数，只需执行Δx←Δx/accuracy即可，则步骤(5.1)中可以直接计算横坐标值x_u=mΔx，而避免对每一个m都需要除以图像精度以转化为图像上的像素坐标，节省除法运算开销。

由于事先并不知道分层物体中的分界线的准确信息，所以步骤(3)使用单一介质的超声成像算法生成一幅图像以便提取第一层与第二层介质间的分界线。由于算法中使用的是第一层介质中的声速，所以对第一层介质所成的图像是准确的，而第一层与第二层介质间的分界线也属于第一层介质的一部分，所以，步骤(4)中提取的该分界线也是准确的。事实上，步骤(3)中使用的超声成像算法既可以是原SAFT技术，也可以是相位迁移技术等任何能对单一介质成像的技术。

成像步骤(4)中需要使用Canny算子边缘提取算法提取介质之间的分界线。Canny边缘提取算子的计算过程分为四步：

步骤(a):图像平滑。对原图像用二维高斯函数平滑滤波得到图像I(x，z)；

步骤(b):计算图像I(x,z)每个像素的梯度G和方向F。采用2×2模板作为对X方向和

Z方向的偏微分的一阶近似，即

$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)$

则梯度大小G和方向F为

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

步骤(c):梯度图像的非极大值抑制。对每个像素点I(x，z)，若其梯度值G(x，z)小于沿着其梯度方向F(x，z)的两个相邻点的梯度幅值，说明该点不是边缘点，则将I(x，z)的灰度设为0；

步骤(d):双阈值化处理。设定好双阈值方法检测和连接边缘需要的低阈值Low和高阈值High，对梯度图像进行双阈值化处理。梯度幅值大于高阈值High的是边缘；梯度幅值小于低阈值Low的不是边缘；梯度幅值介于两者之间的，则判断该像素的八个相邻像素中是否存在大于高阈值High的边缘像素，若存在则是边缘像素，否则不是。

在整个成像计算中，成像步骤(5.1.1)中的半功率波束角β_0.5需要使用M次，然而该值只与超声换能器固有参数有关，本发明只使用一个超声换能器进行探测，所以该值为固定值，为了避免重复计算，可以在程序中预先计算β_0.5←tg(0.5×0.84λ/D)，并保存为全局变量，然后在步骤(5.1.1)中直接调用该变量来计算x_l=x_u-z_l×β_0.5和x_r=x_u+z_r×β_0.5。

此外，步骤(5.1.1)中需要计算半功率波束角的左、右边界线分别与分界线c₁(x，z)的交点B_l(x_l,z_l)、B_r(x_r,z_r)。在实际实现中，Canny算子提取的分界线c₁(x,z)保存为以x坐标值增长顺序的像素点集，在计算交点B_l和B_r时，只需顺序搜索这些像素点即可，如果分界线c₁(x,z)上的某个像素点(x,z)的坐标满足x=x_u-z×β_0.5则为B_l；如果满足x=x_u+z×β_0.5则为B_r。由于在整个成像计算中，对每一个m，均需计算左右两个交点B_l和B_r，但每个m所对应的半功率波束角是互相平行的（如图1(c)），则当前m所对应的左交点B_l必在上一个m的左交点的右边，当前m所对应的右交点B_r也必然在上一个m的右交点的右边，所以步骤(5.1.1)中每一轮计算左交点时，只需从上一轮的B_l处开始搜索分界线c₁(x，z)上的像素点，相应的，每一轮计算右交点时，只需从上一轮的B_r处开始搜索。

成像计算过程的具体算法分为Bresenham直线扫描转换子程序BresenhamLine()和分层物体的合成孔径聚焦超声成像算法主程序SAFT-Refraction()：

图像显示即将成像计算阶段得到的二维图像数据显示在显示器上，可根据需要显示灰度图像或彩色图像。

对于待测对象是规则分层物体的情况，成像步骤(4)和步骤(6)中使用Canny算子边缘提取软件包提取的所有分界线c(x,z)为直线，成像计算步骤不变。

Claims (1)

1.分层物体的合成孔径聚焦超声成像方法，其特征在于，依次含有以下步骤：

步骤(1):构建一个由一台计算机、一个超声换能器、一套定位控制器和一个模数转换器组成的一个基于合成孔径聚焦技术和折射定律的用于对分层物体在深度和水平两个方向形成的纵断面上作无损伤超声成像的***，其中：

所述超声换能器设有：与所述定位控制器的输出端相连的脉冲信号输入端，所述定位控制器的输入端与所述计算机相应的定位控制信号输出端相连，所述超声换能器还设有：与所述模数转换器的输入端相连的回波信号输出端，所述模数转换器的输出端与所述计算机的回波采样信号输入端相连，所述超声换能器由所述定位控制器控制，在被测物体表面以1步长/ms的固定速率移动，所述定位控制器是控制所述超声换能器移动位置的传动装置，其参数由所述计算机输入，

被测物体沿X轴方向的水平长度为X_length，均分为X_length/Δx个区间，Δx为区间长度，也是所述超声换能器沿X轴从坐标点(0,0)起向终点(X_length,0)止每次移动的步长，所述超声换能器每次移动所达到的点称为探测点，共有M个，M＝1+X_length/Δx，序号m=0,1,…，M-1，所述定位控制器在每一个探测点处产生一个TTL晶体管-晶体管逻辑电平脉冲，触发所述超声换能器向被测物体的垂直于X轴的深度方向Z发射一个激励脉冲，随后超声换能器转为接收模式并开始计时，接收从被测物体反射的回波信号，所述模数转换器对所述超声换能器在探测点m处接收到的回波信号进行N次采样并存储到计算机中，采样序号n=0,1，…,N-1，采样频率为f_s，f_s的值为模数转换器预设，记s_m(n)为超声换能器在第m个探测点处的第n次采样得到的采样值，s_m(n)的采样时刻为t＝n/f_s；

步骤(2):所述计算机从n=0开始依序读取探测点m=0处的采样值，然后，重复该过程依次读取m=1，…,M-1各探测点处的采样值；

步骤(3):取v＝v₁，v₁为超声在被测物体的第一层介质中的传播速度，使用合成孔径聚焦超声成像软件包生成深度方向上z₀=0至Z_depth-1区间的纵断面图像，Z_depth为预设的生成图像的长度，即生成的图像在纵向上用像素数表示的深度值；

步骤(4):以步骤(3)得到的z₀至Z_depth-1区间图像块作为输入量，使用Canny算子边缘提取软件包提取第一层介质与第二层介质的分界线c₁(x，z)；

步骤(5):按以下步骤修正所述纵断面上在分界线c₁(x,z)之下至Z_depth-1区间的图像，以消除第一层与其他各层介质间的不同质而导致的误差：

步骤(5.1):取m=0，记第m个探测点为U(x_u,0)，其中，x_u=mΔx/accuracy，accuracy为图像精度，即所生成的图像上相邻两个像素点的间距，依次按以下步骤计算第m个探测点U(x_u,0)处所对应的扇形图像：

步骤(5.1.1):计算超声换能器的半功率波束角β_0.5=0.84λ/D，λ为超声在被测物体中传播时的波长，D为超声换能器的直径，并计算所述半功率波束角的左、右边界线分别与分界线c₁(x，z)的左右两个交点B_l(x_l,z_l)、B_r(x_r，z_r)，其中，x_l=x_u-z_l×tg(0.5β_0.5)且满足c₁(x_l,z_l)=0，x_r=x_u+z_r×tg(0.5β_0.5)且满足c₁(x_r,z_r)=0；

步骤(5.1.2):在分界线c₁(x，z)上取折射点R(x_R,z_R)＝B_l(x_l,z_l)，计算归一化的折射向量其中，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)的归一化的入射向量；

步骤(5.1.3):计算折射矢量

的斜率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与图像边界的交点；

步骤(5.1.4):取x_g=x_f+1，查找分界线c₁(x，z)上当前折射点R_f所在位置的下一个像素点R_g(x_g,z_g)，点R_g的坐标满足c₁(x_g,z_g)=0，取R_g为新的折射点，即R(x_R,z_R)＝R_g(x_g,z_g)，执行步骤(5.1.2)和步骤(5.1.3)，扫描转换新的折射路径R_gE_g，记终点为E_g；

步骤(5.1.5):计算当前折射点R_g的折射路径R_gE_g的终点E_g与上一个折射点R_f的折射路径R_fE_f的终点E_f之间的距离Δd＝|E_fE_g|，将分界线c₁(x，z)上当前折射点R_g与上一个折射点R_f之间的曲线段均分为Δd份，即***Δd-1个点，***点序号记为τ=1,2,…,Δd-1，对τ值以1为步长以τ=1为初始值循环执行下述步骤(5.1.5.1)直到τ=Δd止：

步骤(5.1.5.1):插值计算第τ个***点R_τ的坐标，x_τ=x_f+τ/Δd，z_τ＝z_f+(z_g-z_f)τ/Δd，取折射点R(x_R,z_R)＝R_τ(x_τ,z_τ)，执行步骤(5.1.2)和步骤(5.1.3)，扫描转换从点R_τ出发的折射路径;

步骤(5.1.6):取R_f=R_g，E_f=E_g，执行步骤(5.1.4)和步骤(5.1.5)；

步骤(5.1.7):重复执行步骤(5.1.6)直到x_g=x_r+1止，即处理完分界线c₁(x，z)上介于B_l(x_l,z_l)和B_r(x_r，z_r)之间所有的折射点和折射路径，得到第m个探测点U(x_u,0)处所对应的扇形图像；

步骤(5.2):依次取m=1,…,M-1，重复执行步骤(5.1)，生成所述纵断面上在分界线c₁(x,z)之下至Z_depth-1区间的图像；

步骤(6):以步骤(5)得到的分界线c₁(x,z)之下至Z_depth-1区间的图像块作为输入量，使用Canny算子边缘提取软件包提取所述分界线c₁(x,z)至Z_depth-1之间的第二层介质与第三层介质的分界线c₂(x,z)，按照步骤(5)所述的方法修正所述纵断面上在分界线c₂(x,z)之下至Z_depth-1区间的图像，以消除第二层与其下各层介质间的不同质而导致的误差；

步骤(7):按步骤(6)所述的方法，处理剩余各条分界线，直到所述纵断面中所有分界线处理完止，生成宽度为(M-1)Δx/accuracy+1个像素、长度为Z_depth个像素的所述纵断面的图像。

Images