CN112067698B - 一种时频结合快速全聚焦超声成像方法 - Google Patents

Abstract

本发明公开了一种时频结合快速全聚焦超声成像方法，首先，将采集到的三维全矩阵数据拆分为N个二维矩阵切片。其次，借助二维正/逆傅里叶变换，逐级实现某一切片上每个深度层上的接收波场外推。然后，通过计算每个像素点到发射阵元的延时矩阵，对外推后的时域信号进行插值，实现该切片上的聚焦。最后，对N个二维切片重复上述相移加插值操作并累加，获得时频结合的全聚焦图像。本发明的有益效果：与传统时域全聚焦相比，该方法复杂度低、成像速度快，有望解决全聚焦难以实时成像问题，可有效实现材料内部缺陷的高精度定量、定位、定性，具有良好的推广及应用前景。

Description

一种时频结合快速全聚焦超声成像方法

技术领域

本发明涉及超声无损检测领域，特别涉及一种时频结合快速全聚焦超声成像方法。

背景技术

全聚焦(total focusing method，TFM)是近年来兴起的一种相控阵超声后处理成像技术。与传统相控阵成像技术相比，该技术是一种依赖于全矩阵数据采集(full matrixcapture，FMC)的离线成像技术，克服了传统技术中发射声束数量受限的缺点，并通过离线处理实现整个被检区域的聚焦。研究表明，与传统相控阵所提供的B型、C型、D型、S型视图相比，TFM图像具有更高的检测信噪比和分辨率，其缺陷检测能力和定量精度更高。正因如此，相关学者将全聚焦称之为相控阵超声检测中的“黄金标准”成像技术。

全矩阵数据采集以一发全收的模式多方位的采集检测信号，形成了由发射阵元，接收阵元和时间采样点所组成的三维矩阵。以128个阵元为例，探头按照1、2、…、128号的顺序逐次对单个阵元进行激发。每次激发后，声波经过检测目标并被1～128号全部阵元接收。这样，经128次激发后，共计形成A×128×128个超声脉冲反射信号以三维矩阵形式保存，其中A为由上千个点组成的采样点数。与一般超声成像所用数据相比，用于全聚焦成像的全矩阵数据所携带的检测信号量十分庞大。

时域全聚焦成像时需要重复、繁重的迭代运算，难以在短时间内处理庞大的数据量，导致成像时间十分漫长。因此，当下主流高端相控阵超声检测设备中，全聚焦通常被用作离线图像处理技术，无法满足实时成像的需求。基于此，如何提高全聚焦技术的效率是当下业界亟待解决的问题之一。

对此，相关研究者从算法和硬件两方面入手，试图有效提高全聚焦成像的运算效率。硬件方面，研究者们提出了基于多核心CPU和GPU并行运算方法，即利用多核心硬件进行多路并行迭代运算缩短运算时间。虽然多核心硬件加速方法能够显著提升运算效率，但现阶段这些高性能硬件仅配置在工作站电脑上，短时期内难以普及到便携式超声相控阵探伤仪中。算法方面，研究者们根据全矩阵数据具有对称性的特点，将其简化为半矩阵或三角阵等稀疏矩阵，用以减少数据量进而缩短运算时间。然而，这种减少数据量的方法依旧未能避免逐像素点的迭代运算，仅能将运算时间缩短几分之一，运算效率提升不够明显。

综上，成像速度慢是现阶段时域全聚焦技术发展中亟待解决的问题，更有效的快速超声全聚焦成像方法期待被提出。

发明内容

本发明的目的在于提供一种时频结合快速全聚焦超声成像方法，利用相移加插值的时频结合方法代替全时域繁重的延时叠加运算，显著缩短全聚焦图像的计算时间，为材料内部缺陷的高精度定量、定位、定性提供快速、有效的方法。

本发明的目的是这样实现的，一种时频结合快速全聚焦超声成像方法，包括如下步骤：

1)将已采集到的全矩阵数据N_t×N_R×N_T拆分为N_T个N_t×N_R的二维子矩阵切片p(t,x)，每个切片均表示单阵元发射，分别被1～N号阵元接收的N个信号，其中，N_t为采样点个数，N_R为接收阵元个数，N_T为发射阵元个数，在所述子矩阵数据p(t,x)中，每一行表示水平坐标x上信号幅值的变化，即信号随接收阵元位置的变化,每一列表示深度z上信号幅值的变化，即不同接收采样点上信号幅值的变化；

2)利用相移加插值的方法，分别对各子矩阵数据进行虚拟聚焦；

3)将每一个发射阵元下接收到的回波数据p_i(t,x)循环操作步骤3)，共可获得1～32号子矩阵二维聚焦图像，最后将1～32号所有子矩阵二维聚焦图像进行图像融合，获得基于相移加插值的时频结合超声全聚焦图像。

进一步的,所述步骤二利用相移加插值的方法，分别对各子矩阵数据进行虚拟聚焦，具体实施方法如下：

2.1将上述所得到的子矩阵p_i(t,x)进行二维傅里叶变换，将时域空间p_i(t,x)转换为频域空间p_i(ω,k_x)，其中i为发射阵元序号，时域空间p_i(t,x)中，t为信号沿垂直坐标z上的传播时间，x为阵元水平方向位置上的变化，频域空间矩阵p_i(ω,k_x)中，ω为相位沿垂直方向z上的变化率，即角频率变化，k_x为频域空间中水平方向x上的波数矢量；

2.2利用已知量ω、k_x以及声速c计算接收信号不给发射延时下垂直方向z上的波数矢量k_z，k_z＝k(1+cos(θ)),其中k＝ω/c，kcos(θ)＝sqrt((ω/c)²-k_x ²),θ为像素点与接收阵元之间的夹角，则此时的迁移因子F＝exp(jk_z*dz)＝exp[j(ω/c+sqrt((ω/c)²-k_x ²)]*dz，式中，j为虚数单位，dz为垂直坐标z上的相邻信号采样点之间的间隔，即为每次外推的深度间隔，可由采样频率f_s、纵波声速c算出，满足如下关系式：dz＝c/f_s；

2.3将迁移因子F与p_i(ω,k_x)相乘并进行二维傅里叶逆变换，得到p_i(ω,k_x)外推深度dz的计算结果p_i(t,x)，将每一层相乘的频域波场作为下一层外推的起始波场逐层外推，直至最终的外推深度层Z_end，每一层迁移后的时域波场保存于成像矩阵IP_i中；

2.4对得到的IP_i矩阵进行发射延时，根据阵元位置、像素点位置以及距离公式，可得发射延时t_delay，所求延时矩阵的大小与像素点个数一致，利用t_delay对IP_i矩阵进行时间维度t上的插值，即得到经过发射延时之后的频域波场外推的结果。

本发明的有益效果是：基于相移和插值技术的时频结合超声全聚焦成像技术，与传统时域上的全聚焦成像相比，算法复杂度低、成像速度快，应用于实践有望解决现有相控阵超声波探伤仪难以实时全聚焦成像的问题，可有效实现材料内部缺陷的高精度定量、定位、定性，具有良好的应用前景。

附图说明

图1是本发明中超声全矩阵数据采集系统示意图；

图2是本发明中时频结合的全聚焦算法流程图；

图3是本发明基于相移加插值的时频结合超声全聚焦图像；

图4是常规全聚焦图像。

具体实施方式

以下结合附图和实施例对本发明作进一步说明。参见图1至图4，一种时频结合快速全聚焦成像技术，以32阵元相控阵探头，对60mm厚铝制试块中30mm深度的Φ2边钻孔进行检测为例，其中操作步骤如下：

1)如图1所示，通过相控阵探头1、多通道选择器2、超声信号发射接收器3、数字示波器4和计算机终端5搭建全矩阵超声信号采集系统。其中，相控阵探头1与多通道选择器2相连接，多通道选择器2与超声信号发射接收器3之间通过发射接口TX、接收接口RX连接，多通道选择器2与计算机终端5之间通过网线LAN连接，超声信号发射接收器3与数字示波器4之间通过OUT和SYMC接口连接，相控阵探头1放置于检测工件6上，超声波在检测工件6中的声速c为6300m/s。本例中相控阵探头1阵元数为32，阵元宽度为0.9mm，阵中心间距为1mm，中心频率为5MHz。在计算机终端5的系统控制界面中设置采样频率20MHz进行全矩阵数据采集。按顺序逐次激发1、2…n…32号阵元，每次激发后，声波经过检测目标后分别被1～32号阵元接收。32次激发后共计形成32×32个超声脉冲A扫反射信号，每个A扫信号共768个采样点，以三维全矩阵数据768×32×32形式保存于计算机终端5内，该数据命名为fmc_data。矩阵中第一维度表示采样点时间序列，第二维度表示接收阵元序列，第三维度表示发射阵元序列；

2)利用MATLAB软件读取全矩阵数据fmc_data，其数据量为768×32×32。利用for循环i＝1:32将三维全矩阵数据进行拆分，定义data＝fmc_data(:,:,i)，其中n代表发射阵元序号，取值范围0～32。通过上述操作，将三维全矩阵数据fmc_data拆分为32个二维子矩阵切片data，每个子矩阵切片中含有32个超声脉冲反射信号，用以表示单阵元n发出，被1～32阵元号全部阵元接收的回波数据。子矩阵数据data中，每一行表示水平坐标x上信号幅值的变化，即信号随接收阵元位置的变化；每一列表示深度z上信号幅值的变化，即信号随接收采样点的变化；

3)步骤2)中的data是幅值关于时间t和水平位置x的函数，表示为p(t,x),在fori＝1:32循环内，通过时频结合的全聚焦算法对data(i)进行聚焦，计算流程如图2所示，具体操作步骤实施如4)～7)所描述；

4)利用MATLAB中的fft2函数对p_i(t,x)进行二维傅里叶变换，将时域空间p_i(t,x)转换为频域空间p_i(ω,k_x)其中i为发射阵元序号。时域空间p_i(t,x)中，t为信号沿垂直坐标z上的传播时间，x为阵元水平方向位置上的变化，频域空间矩阵p_i(ω,k_x)中，ω为相位沿垂直方向z上的变化率，即角频率变化，k_x为频域空间中水平方向x上的波数矢量。其中，ω与t之间满足关系式ω＝2πf_s/N_t＝2π/t/N_t，f_s为采样频率，N_t为采样点的个数的2的整数次幂，本例的采样频率f_s为20MHz，N_t＝2^nextpow2(768)为1024。k_x与x之间满足关系式k_x＝2πx_s/N_x，x_s为相邻两阵元之间的阵中心间距，N_x为接收阵元的个数的2的整数次幂，本例中相邻两阵元之间的中心距离为1mm，N_x＝2^nextpow2(32)为32，此时频域空间内p_i(ω,k_x)矩阵大小为1024×32；

5)利用已知量ω、k_x以及声速c计算接收信号不给发射延时下垂直方向z上的波数矢量k_z，k_z＝k(1+cos(θ)),其中k＝ω/c，kcos(θ)＝sqrt((ω/c)²-k_x ²),θ为像素点与接收阵元之间的夹角，则此时的迁移因子F＝exp(jω(1+cos(θ))/c)*dz＝exp[j(ω/c+sqrt((ω/c)²-k_x ²)]*dz，式中，j为虚数单位，dz为垂直坐标z上的相邻信号采样点之间的间隔，即为每次外推的深度间隔，可由采样频率f_s、纵波声速c算出，满足如下关系式：dz＝c/f_s。将迁移因子F与p_i(ω,k_x)相乘并进行二维傅里叶逆变换，得到p_i(ω,k_x)外推深度dz的计算结果p_i(t,x)，将每一层相乘的频域波场作为下一层外推的起始波场逐层外推，直至最终的外推深度层Z_end，本例中外推的次数n_z＝Z_end/dz为768，每一层迁移后的时域波场保存于成像矩阵IP_i中，IP_i矩阵的大小为768×32；

6)对步骤5)得到的IP_i矩阵进行发射延时，根据相邻两阵元之间的阵中心间距可计算得到阵元的水平位置x_p＝(0:(32-1))*xs，像素点之间的水平间隔Δx＝xs和深度间隔Δz＝c/f_s，则像素点的水平位置x_m＝(0:Δx:l)和深度位置z_m＝(0:Δz:h),l和h分别为成像区域的宽度和深度。利用meshgrid函数将水平位置和深度位置网格化，得到矩阵化的像素点位置。根据得到的像素点位置以及距离公式，可得发射延时t_delay＝[sqrt((x_m-x_p)²+(z_p)²)+(z_p)]/c，所求延时矩阵的大小与像素点个数一致。利用t_delay对IP_i矩阵利用interp(t,IP_i,tdalay)函数进行时间维度t上的插值，即得到经过发射延时之后的频域波场外推的结果；

7)将每一个发射阵元下接收到的回波数据循环操作步骤4)-6)，共可获得1～32号子矩阵二维聚焦图像，最后将1～32号所有子矩阵二维聚焦图像进行图像融合，获得如图3所示的图像，即为本发明基于相移加插值的时频结合超声全聚焦图像。与图4所示的常规全聚焦图像相比，时频结合得到的全聚焦图像具有更高的横向分辨率，成像质量得到了改善。以I5 4310U型CPU，4G计算内存为测试平台对比两种方法的运行时间，MATLAB 2015a软件单核CPU计算环境下，本发明提出的算法仅需0.47s，而常规全聚焦算法需要362.15s。若考虑后续算法优化和硬件加速，本发明能够满足高质量实时成像条件，有望解决现有相控阵超声波探伤仪难以实时全聚焦成像的问题。

Claims (1)

1.一种时频结合快速全聚焦超声成像方法，其特征在于包括如下步骤：

第一步，将已采集到的全矩阵数据N_t×N_R×N_T拆分为N_T个N_t×N_R的二维子矩阵切片p(t,x)，每个切片均表示单阵元发射，分别被1～N号阵元接收的N个信号，其中，N_t为采样点个数，N_R为接收阵元个数，N_T为发射阵元个数，在所述子矩阵数据p(t,x)中，每一行表示水平坐标x上信号幅值的变化，即信号随接收阵元位置的变化,每一列表示深度z上信号幅值的变化，即不同接收采样点上信号幅值的变化；

第二步，利用相移加插值的方法，分别对各子矩阵数据进行虚拟聚焦；

第三步，将每一个发射阵元下接收到的回波数据p_i(t,x)循环操作步骤二，共可获得1～32号子矩阵二维聚焦图像，最后将1～32号所有子矩阵二维聚焦图像进行图像融合，获得基于相移加插值的时频结合超声全聚焦图像；

步骤二利用相移加插值的方法，分别对各子矩阵数据进行虚拟聚焦，具体实施方法如下：

2.1将步骤一所得到的子矩阵p_i(t,x)进行二维傅里叶变换，将时域空间p_i(t,x)转换为频域空间p_i(ω,k_x)，其中i为发射阵元序号，时域空间p_i(t,x)中，t为信号沿垂直坐标z上的传播时间，x为阵元水平方向位置上的变化，频域空间矩阵p_i(ω,k_x)中，ω为相位沿垂直方向z上的变化率，即角频率变化，k_x为频域空间中水平方向x上的波数矢量；

2.2利用已知量ω、k_x以及声速c计算接收信号不给发射延时下垂直方向z上的波数矢量k_z，k_z＝k(1+cos(θ)),其中k＝ω/c，kcos(θ)＝sqrt((ω/c)²-k_x ²),θ为像素点与接收阵元之间的夹角，则此时的迁移因子F＝exp(jk_z*dz)＝exp[j(ω/c+sqrt((ω/c)²-k_x ²)]*dz，式中，j为虚数单位，dz为垂直坐标z上的相邻信号采样点之间的间隔，即为每次外推的深度间隔，可由采样频率f_s、纵波声速c算出，满足如下关系式：dz＝c/f_s；

2.3将迁移因子F与p_i(ω,k_x)相乘并进行二维傅里叶逆变换，得到p_i(ω,k_x)外推深度dz的计算结果p_i(t,x)，将每一层相乘的频域波场作为下一层外推的起始波场逐层外推，直至最终的外推深度层Z_end，每一层迁移后的时域波场保存于成像矩阵IP_i中；

2.4对得到的IP_i矩阵进行发射延时，根据阵元位置、像素点位置以及距离公式，可得发射延时t_delay，所求延时矩阵的大小与像素点个数一致，利用t_delay对IP_i矩阵进行时间维度t上的插值，即得到经过发射延时之后的频域波场外推的结果。