基于椭圆弧扫描转换的相控阵超声回波成像方法

摘要

基于椭圆弧扫描转换的相控阵超声波回波成像方法属于快速超声波回波成像技术领域，其特征在于，利用椭圆弧的扫描转换技术，对现有时域相控阵中的每个发送/接收阵元组获得的回波数据进行椭圆弧绘制操作，目标反射物为多个椭圆弧的交点，使其像素值有效叠加，高于其周围像素，借此得到目标反射物的位置。依次由构造一个由计算机、相控阵探头和模/数转换器串接而成的成像系统，各次探测过程中数据的采集，利用椭圆弧扫描转换技术计算以发送/接收阵元为焦点的纵断面上的椭圆弧构成，最终实现对被探测物体的超声波反射信号的快速准确的成像。

说明书

技术领域

本发明涉及超声波无损探测技术、超声成像技术、相控阵技术和椭圆的扫描转换技术，实现对物体表面形状和内部结构的快速准确成像。

背景技术

超声波相控阵技术(Phased Array)是超声波无损探测领域一项新的超声成像技术，相对于单个换能器，相控阵更加灵活，并且对物体表面形状和内部缺陷更加敏感，因此备受研究者亲睐，成为近年来研究的热点，在重工业、油气管道、能源和航空等领域均有广泛应用。

相控阵技术的基本思想是利用脉冲-回波(pulse-echo)测量机制，使用一组超声相控阵阵元(超声换能器)阵列，每次一个相控阵阵元向探测物体发送脉冲扫描信号，并采用延时叠加(DAS)方法(时间延迟或相位延迟)对得到的脉冲回波信号进行聚焦成像。全矩阵捕获(Full Matrix Capture)是面向相控阵的一种新型数据获取方式，能够获取所有发送/接收阵元组的时域数据，能够获得更高的图像分辨率，并能够对物体内部缺陷特征进行更准确的判断。全矩阵捕获的工作模型如图1(a)所示，每次选取一个相控阵阵元向物体以与X轴的夹角为θ的方向发射超声波信号，对物体内部进行探测，同时整个相控阵阵列接收从物体的分界面或者内部反射物反射回的回波信号并采样保存，重复这种“一个发送，所有接收”的过程直至相控阵所有阵元均发送过声波信号，最后对全矩阵捕获所获得的全矩阵数据进行后处理并显示图像。多种后处理算法根据处理方式的不同可分为时域和频域两大类：

时域相控阵超声成像算法中应用最为广泛且成像效果最佳的为全聚焦算法，该算法流程如图1(b)所示，为了在目标图像任一像素点(x,z)聚焦，时域相控阵技术将相控阵阵列中所有(发送,接收)阵元组合接收的回波信号进行延时叠加处理：设为第i步采样过程中，相控阵阵元E_j接收到的采样回波信号(即对应发送阵元/接收阵元组(E_i,E_j)的回波信号)，为采样序号，相控阵阵元E_j关于像素点(x,z)的采样延迟为

其中，v为超声波在介质中的传播速度，r_i,j(x,z)为像素点(x,z)距相控阵阵元E_i和E_j的距离之和。对应的采样序号为f_s为采样频率。相控阵有效长度L内所有的延时构成一条延时曲线。L的计算公式为

L＝0.84λz/d(2)λ为超声在介质中的波长，d为相控阵阵元的直径，则在(x,z)处的像素值为

其中，N为相控阵中阵元数目，ω_i,j(x,z)为变迹函数，由下式计算

ω_i,j(x,z)＝ω_i(x,z)·ω_j(x,z)

其中，

对于阵元E_i(x_i,0)，在深度z时，其有效孔径的范围为而对于位于有效范围之外的点，变迹函数值为0，其像素值不会受E_i影响。如图1(d)所示，红色区域表示阵元E₀的声场范围，绿色区域表示阵元E₁的声场范围，黄色区域表示两者公共声场范围，对于发送接收阵元组(E₀,E₁)所对应的采样信号，仅能够影响黄色区域内的像素点。

频域相控阵算法基于波动方程反演理论。近期有学者将反射地震学(Reflection Seismology)中的迁移技术(Migration Technique)与频域相控阵技术相结合，得到了相位迁移超声成像方法(Phase Shift Migration)。该方法将在相控阵阵元获得的回波信号作为波动方程的边界条件，然后使用频域中的相移符号对位于不同深度的声场进行假设检测，从而外推得到整个声场。算法主要包括两个步骤：第一步对相控阵阵元获得的时域回波信号数据进行二维傅里叶变换，得到二维频谱；第二步是获得在图像边界的频谱，然后对上次计算得到的二维频谱作相位平移，然后进行傅里叶反变换，依次得到不同深度值的图像。

从时域相控阵技术的原理可知，在成像过程中需要计算图像中的所有像素点与相控阵阵元扫描位置点各点对之间的距离，算法复杂度高，并且涉及到均方根运算，使得算法的运行时间很长。而在频域相位迁移技术中，计算过程中需要对数据进行正向和逆向的傅里叶变换，计算效率也不太高。

因此，与相控阵相关的成像计算都比较耗时且成像精度误差较大，不能满足实际工程应用对成像速度和准确度的要求。所以，相控阵技术还需要进一步提高成像效率和成像精度。

分析时域相控阵技术的原理，可以发现该成像方法其实是一种逆向的计算过程：需要先确定图像上的像素点，然后，对于所有的发送/接收相控阵阵元对，再计算像素点到两个阵元的距离之和r_i,j及传输延时t_i,j，将延时对应的采样信号值S_i,j(n_t)累加得到像素值。计算过程中，被测物体中反射物的位置是未知的，通过对图像中所有的像素点进行动态聚焦和延时累加，获得整幅图像。位于目标反射物的像素点的延时采样信号得到一致的叠加，实现信号的聚焦。而对于周围像素点，信号的累加是无序的，因此在图像中，反射物所对应处累加得到的像素值要明显大于其他像素点。

全聚焦算法与合成孔径聚焦算法类似，主要基于延时叠加原理。不同之处在合成孔径聚焦方法中，发送阵元与接收阵元为同一阵元，则对于含有N个阵元的相控阵列，合成孔径聚焦将获得N组A扫数据。为了提升合成孔径聚焦算法的运行效率，有研究人员提出了基于圆弧扫描转换的合成孔径聚焦算法(DAB-SAFT)。从正向来对该过程重新解释，即计算所有获得的A扫数据对目标图像的影响结果，将整个成像过程转换为多次圆弧扫描转换操作，使合成孔径聚焦获得了极大的性能提升。但是全矩阵捕获方法获得了相控阵中所有发送接收阵元组的所对应的A扫信号，即对于含有N个阵元的相控阵列，共有N²组A扫信号，称为全矩阵数据。对于发送接收阵元不同的情况，该算法无法处理，因此无法应用到相控阵超声回波成像中。

基于椭圆弧扫描转换的相控阵成像方法，有如下特点:(1)该算法可以应用到发送-接收阵元不同的情况，也能对发送接收阵元相同的情况进行处理，DAB-SAFT仅是本算法的一个特例。(2)本算法通过正向解释全聚焦过程，将成像过程转换为多次椭圆弧扫描转换操作，去除了耗时的均方根运算，同时利用回波数据的稀疏性降低了算法复杂度，极大提升成像速度。

对于每个发送/接收阵元对采样得到的每个数据，不仅作用于反射物对应的像素点，也作用于无反射物的像素点。如图1(c)中，对应发送阵元为E_i、接收阵元为E_j且采样延时为t_i,j的采样数据不仅参与了p(x,z)点的成像计算，也参与了曲线段arc_i,j(r_i,j)上其他点如p’(x’,z’)的计算，r_i,j＝v·t_i,j/accuracy。其中，arc_i,j(r_i,j)是位于阵元E_i和E_j的声场范围内的一段曲>i,j(n_ti,j)的影响。曲线arc_i,j(r_i,j)上的点满足以下条件：曲线上的点距阵元E_i和E_j的距离之和r_i,j均相同、延时均为t_i,j。所以曲线段arc_i,j(r_i,j)是一段以E_i和E_j为焦点，长轴长度为r_i,j的椭圆弧。

另一方面，若从正向去综合理解整体计算过程，采样数据S_i,j(n_t)在整幅图像的成像计算中的功效就相当于以其数据值ω_i,js_i,j(n_t)/r_i,j在图像中画了一段椭圆弧arc_i,j(r_i,j),即椭圆弧arc_i,j(r_i,j)上各像素点的像素值累加ω_i,js_i,j(n_t)/r_i,j。因此，S_i,j(n_t)所在的延时曲线上的所有数据(S_0,0(n_t)至S_N-1,N-1(n_t))所对应的椭圆弧(arc_0,0(r_0,0)至arc_N-1,N-1(r_N-1,N-1))在图像中的交点即为原逆向计算过程中的目标像素点(x,z)。对于目标反射物，会有多个椭圆弧经过该反射物，使其像素值高于周围像素，以突出反射物所在位置，如图1(d)所示，当发送阵元为E₀时，目标反射物(x,z)有arc_0,0(r_0,0)，arc_0,1(r_0,1)，arc_0,2(r_0,2)等椭圆弧经过该点，该点像素值得到有效叠加。

若在整幅图像中仅考虑第i步探测过程中，阵元E_j的所有采样数据S_i,j的功效，则对应于一幅子图像I_i,j(如图1(e))，图像里包括多个以E_i和E_j为焦点，但长轴不同的椭圆弧。原DAS计算公式(3)则可以重新理解为各扫描位置处所对应的扇形图的叠加，图像中像素点值为：

结果与(3)式相同，其中，N为相控阵中阵元的总数。

在原算法过程中，图像上的任一像素点p(x,z)，在第i步探测过程中阵元E_j接收的采样数据S_i,j对该点像素值的贡献为为使新算法与原算法有相同的结果，在子图像I_i,j中的所有像素点p(x,z)，有且仅有一个椭圆弧经过该点。如图2所示，子图像I_t,r中的椭圆弧均关于直线AM：x＝x_c＝(x_i+x_j)/2对称，其中A(x_c,0)，M(x_c,Z_depth-1)，线段AM上所有像素点P_k(x_c,z_k),z_k＝0,1,…,Z_depth-1，均有一个椭圆弧经过。所以子图像I_t,r中共包括Z_depth个椭圆弧，每个椭圆弧以(x_i,0)和(x_j,0)为焦点，且经过点P_k(x_c,z_k),z_k＝0,1,…,Z_depth-1。

因此，原时域相控阵技术可以通过按照各发送/接收阵元对的位置在图像上依长轴r_i,j的大小顺序逐次画椭圆弧的方法来实现。而在计算机图形学中，为了在光栅显示设备上显示几何椭圆弧，已开发了多种成熟的椭圆弧扫描转换的技术。其中应用最为广泛的是中点画椭圆算法，利用该算法，可以快速得到椭圆弧上各个像素点的位置，无需计算成像点与各阵元位置点之间的距离，避免了均方根运算，能节省大量的计算开销。

算法中每个发射和接收换能器的循环叠加计算完成之后，将每个像素的值显示出来；整个图像在内部缺陷和分界面的像素处会出现明显的高亮度值，并会形成连续的白色边界线；再通过图像的边缘检测和提取(算法)便可得到被测物体的内部缺陷或分界面的形状曲线。

发明内容

本发明的目的在于提出一种基于椭圆弧扫描转换的时域相控阵超声成像技术的实现方法，提高成像速度和成像精度。

本发明的特征在于，将椭圆的扫描转换技术扩展为超声波的椭圆弧扫描转换技术，并使用椭圆弧段的光栅化方法实现时域相控阵超声快速成像技术，无需计算成像点与各阵元位置点之间的距离，避免了均方根运算，节省了大量的计算操作。

本发明的基于椭圆弧扫描转换的时域相控阵超声波回波成像方法，其特征在于，依次含有以下步骤：

步骤(1):构建一个基于椭圆弧扫描转换的能在被测物体的深度方向形成的纵断面上作无 损伤检测的时域相控阵超声波回波成像系统，以下简称系统，它包括一台计算机、一个相控阵探头、一个相控阵控制器和一个模数转换器组成，其中：

所述的相控阵控制器的输入端与计算机控制信号的输出端相连，相控阵探头的输入端与相控阵控制器的输出端相连，相控阵的回波信号输出端与所述模数转换器的输入端相连，所述模数转换器的输出端与所述计算机的回波采样信号输入端相连，

在被测物体的上表面设置X-Z坐标系，被测物体沿X轴方向的水平长度为M_x，向右为正，纵轴Z向下为正，作为目标反射物被测物体纵断面上深度方向的坐标为(x,z)

所述相控阵探头共有N个相控阵阵元(即：超声换能器)，以下简称阵元，各阵元沿X轴正方向等间距Δx排列组成一个相控阵阵列，各阵元用E_i表示，i＝0,1,2,…,N-1，i既是阵元编号，也是检测过程中各个探测步骤的序号，图像中相邻两个像素点的间距称为图像精度accuracy，阵列中首个阵元的坐标为(x₀,0)，第E_i个阵元的坐标为(x_i＝x₀+i*Δx/accuracy,0)，目标图像水平方向的像素点数目X_length＝M_x/accuracy，

计算超声换能器的半功率波束角β_0.5＝0.84λ/d，λ为超声在被测物体中传播时的波长，d为超声换能器的直径，对于所有的阵元均相同，在深度为z时，阵元有效孔径长度为：L_z＝z*β_0.5；

步骤(2)：所述系统在基于椭圆弧扫描转换的时域相控阵超声波成像方法控制下对被测物体沿深度方向形成的纵断面依次按如下步骤实现：

步骤(2.1)：第i＝0步，按如下步骤探测:

步骤(2.1.1)：所述计算机控制所述相控阵控制器产生一个晶体管-晶体管逻辑电平TTL脉冲，触发首个阵元E_i＝E₀，使其向被测物体的垂直于X轴的深度方向Z发射一个激励脉冲；

步骤(2.1.2)：包括E_i在内的所有N个阵元都转为接收模式并开始计时，接收从被测物体反射的回波信号，所述模/数转换器对本次探测过程中各阵元E_j接收到的回波信号进行采样，共计N_t次，采样序号为n_t，n_t＝0,1,2,….,N_t-1，采样频率为f_s，f_s的值为模数转换器预设，记S_i,j(n_t)为相控阵在第i＝0步探测过程中，第j个阵元E_j第n_t次采样得到的采样值；每个阵元的总采样次数N_t为：

N_t＝2Z_depth·accuracy·f_s/v

v为超声波在被测物体中的传播速度，

Z_depth为生成目标图像的垂直方向的像素数目，对于不同的聚焦深度

本算法同样适用，通过设置Z_depth来实现；

步骤(2.2)：所述计算机重复步骤(2.1.1)-步骤(2.1.2)，依次读取第i＝1,2,…,N-1步所述模/数转换器输入的全部N个相控阵阵元的回波信号采样值；

步骤(3)：对于第i＝0步探测过程中,第j＝0个阵元接收的回波采样信号S_i,j(n_t)，n_t＝0,1,…,N_t-1，求解其对应的子图像I_t,r，即采样信号对目标图像上所有像素点的贡献值。

步骤(3.1)：依次按以下步骤计算所述纵断面上的以坐标点(x_i,0)和(x_j,0)为焦点的一系列椭圆弧。

步骤(3.1.1)：如图2所示，取深度方向上的初始坐标值z_α＝0，椭圆中心点坐标为(x_c＝(x_i+x_j)/2,0)，计算路径T_i→p(x_c,z_α)→T_j的传输距离：

计算超声在被测物体中传播距离为r_i,j(x_c,z_α)时所需要的时间

t_i,j(x_c,z_α)＝r_i,j(x_c,z_α)·accuracy/v

及该时间所对应的采样序号n_t＝t_i,i(x_c,z_α)·f_s；

步骤(3.1.2)：判断采样信号S_i,j(n_t)是否为非0值，若是，执行步骤(3.1.3)，否则，跳过步骤(3.1.3)直接执行步骤(3.1.4)；

步骤(3.1.3)：采样信号S_i,j(n_t)所对应的椭圆弧为arc_i,j(r_i,j(x_c,z_α))，该椭圆弧以(x_i,0)和(x_j,0)为焦点且经过点(x_c,z_α)。其中椭圆长半轴为a＝r_i,j(x_c,zα)/2，短半轴为b＝0，椭圆弧方程为

f(x,z)＝b²(x-x_c)²+a²z²-a²b²＝0

调用中点画椭圆算法得到该椭圆弧上所有像素点的坐标，对于所有满足z≥0的像素点p(x_k,z_k)，k＝0,1,2,…,n_k-1，n_k为满足条件的像素点数目，其像素点的值累加其中，变迹函数ω_i,j(x,z)为：

ω_i,j(x,z)＝ω_i(x,z)·ω_j(x,z)

步骤(3.1.4)：依次取深度方向上的坐标值z_α＝1,2,…,Z_depth-1，重复执行步骤(3.1.1)至步骤(3.1.3)；

步骤(4)：依次取接收阵元j＝1,2,…,N-1，执行步骤(3)；

步骤(5)：依次对于探测步骤i＝1,2,…,N-1,执行步骤(3)-步骤(4),生成X_length×Z_depth的纵断面图像。

步骤(6)：运用边缘检测算法提取出纵断面图像中像素值显著变化的点，这些像素点构成了被测物体中缺陷(分界面)的形状和位置。

算法的计算复杂性分析：

假设探测物体的水平长度包括像素数目为X_length，垂直方向包括像素数目为Z_depth，相控阵探头中阵元的数目为N。以一个像素点的计算为例，原时域相控阵算法对于相控阵列中所有发送/接收阵元组合，都要计算“发送阵元→像素点→接收阵元”的传输距离，然后运算得到信号传输时间，根据时间查找回波信号的强度值。每一次距离计算都需要作均方根运算，设每次均方根运算的计算时间为T_RMS，每次查找回波信号值的时间为T_L，则原算法的计算时间为O(X_length*Z_depth*N²*(T_L+T_RMS))。而运用本算法，每次仅需对阵元接收的非零信号进行运算，并>s，每次光栅化椭圆弧所涉及的均方根运算只有一次(计算椭圆弧初始点的位置)。若椭圆弧上像素点的平均数目为n_e，椭圆弧扫描转换时，光栅化每个像素的加法运行时间为T_A，因此每次绘制椭圆弧操作的时间为T_RMS+n_e*T_A，则本算法的运行时间为O(N²*n_s*(T_RMS+n_e*T_A))。假设X_length，Z_depth属于同一数量级(用O(X_length)表示)，N为常数(通常N<<X_length),信号最远传输距离为最近传输距离为0，为简化计算，假设平均传输距离为则因此n_e与X_length，Z_depth属于同一数量级(用O(X_length)表示),n_s为常数(用O(1)表示)。原算法的复杂度为O(N²*X_length²)，本算法复杂度为O(N²*X_length)。实验表明，本方法的计算速度相对于传统的时域相控阵成像方法有1个数量级的提升。

实验结果分析：

本发明与现有的技术相比，优点在于成像速度快且成像更准确。对于图3(a)中的被测物体，物体中有五个等间距排列的直径为1mm的孔洞。原算法和新算法在Core^TMi5-4210CPU@1.7GHz和4GB内存的机器上运行。图3(b)为原算法的运行结果，图3(c)为新算法的运行结果。两图中的高像素值区域(即白色区域)表明此处存在缺陷，对比3(b)和3(c)中的低像素区域位置，发现两图结果完全相同。图3(d)为上述结果图像的灰度对比图，可以看出，两者灰度图也几乎完全相同，从而说明新算法与原算法有相同的运行结果。而原算法运行时间为10.423s，基于椭圆弧扫描转换的相控阵成像方法的运行时间为0.224s，相对于原算法有46倍的性能提升。

附图说明

图1是时域相控阵超声成像技术的工作模型及原理图：1(a)是全矩阵捕获流程；1(b)是原时域相控阵技术的逆向计算过程原理图；1(c)是对时域相控阵技术的正向解释说明图；1(d)是第0步探测过程中，经过目标反射物的椭圆弧：其中红色区域表示阵元E₀的声场范围，绿色区域表示阵元E1的声场范围，蓝色区域表示阵元E2的声场范围，黄色区域表示E₀和E₁公共声场范围，紫色区域表示E₀和E₂公共声场范围；1(e)是第i步探测过程中，阵元E_j采样得到的所有数据在原时域相控阵技术中的功效图。

图2是计算子图像中包含的椭圆弧的数目示意图。

图3是实验示意图：3(a)为被测物体结构示意图；3(b)为原算法结果图像；3(c)为基于椭圆弧扫描的成像方法结果图像；3(d)为图像(b)和(c)的灰度直方对比图。

图4是本超声成像系统流程示意图。

图5是本超声成像硬件系统结构图。

图6是超声阵元工作示意图。

图7是基于椭圆弧扫描转换的时域相控阵超声成像方法总体流程图。

图8是基于椭圆弧扫描转换的时域相控阵超声成像方法伪代码。

具体实施方式

本发明的具体实施过程包含三部分(如图4)：超声数据获取、成像计算和图像显示。硬件平台系统结构图如图5所示，超声成像系统由一台计算机、一个相控阵探头、一个相控阵控 制器和一个模数转换器组成，所述的相控阵控制器的输入端与计算机控制信号的输出端相连，相控阵探头的输入端与相控阵控制器的输出端相连，相控阵的回波信号输出端与所述模数转换器的输入端相连，所述模数转换器的输出端与所述计算机的回波采样信号输入端相连。

被测物体沿X轴方向的水平长度为M_x，相控阵探头中一维相控阵列阵元的数目为N，阵元为E_i(i＝0,1,2,,…,N-1)，相邻阵元的间距为Δx，相控阵第0个阵元的坐标为E₀(x₀＝0,0)，则第i个阵元E_i的坐标为(x_i＝x₀+i*Δx/accuracy,0)，其中accuracy为图像精度，物体水平方向的像素值为X_length＝M_x/accuracy。探测过程共分为N步。第i(i＝0,1,2,…,N-1)步探测过程中，相控阵控制器产生一个TTL(晶体管-晶体管逻辑电平)脉冲，触发相控阵中阵元E_i向被测物体的垂直于X轴的深度方向Z发射一个激励脉冲，随后所有阵元(E₀,E₁,…,E_N-1)转为接收模式并开始计时，接收从被测物体反射的回波信号，所述模数转换器对每个阵元(E₀,E₁,…,E_N-1)在第i步探测过程中接收到的回波信号进行N_t次采样并存储到计算机中,采样序号n_t＝0,1,…,N_t-1，采样频率为f_s，其值为模数转换器预设，记s_i,j(n_t)为第i步探测过程中，阵元E_j接收到的回波信号在t＝n_t/f_s时刻的采样值。其中N_t＝2Z_depth·accuracy·f_s/v，v为超声在被测物体中的传播速度，Z_depth为被测物体垂直方向的像素值，其值为系统预设。

成像计算就是将全矩阵捕获方式获得的采样数据作为计算机输入，然后按前述的成像步骤计算被测物体的纵断面图像，总体流程图参见图7。

成像步骤(3.1)中，计算传输延时所对应的采样序号n_t＝f_s·t＝f_s·2a/v时，得到的值n_t可能为小数，此时获得与之相邻的两个整数对应的采样信号进行插值得到采样信号。例如，假设得到的采样序号为n_t＝30.4，与之相邻的整数对应的采样信号为s_i,j(30)和s_i,j(31)，则n_t对应的采样信号为s_i,j(30)+(s_i,j(31)-s_i,j(30))*(30.4-30)。

图3所述的实验中，相控阵列中阵元数目为32，第0个阵元坐标为(0,0)，阵元间隔为1.2mm，图像精度为accuracy＝0.05mm，阵元E_i的坐标为(i*24,0)。阵元直径为d＝1mm，阵元频率为f_t＝3.00MHz，半功率波束角为β_0.5＝0.84λ/d＝0.431，声波在水中传播速度为v＝1540m/s，采样频率为f_s＝100.00MHz。首先使用图1(a)的流程获得所有(发送,接收)阵元组合的回波采样数据，然后利用图8的成像算法对采样数据进行处理，获得最终X_length×Z_depth的结果图像。

图像显示即把计算得到的二维图像显示在相应的显示设备上。用边缘检测算法对结果图像进行处理，得到缺陷所处的位置。