CN109164453A - 一种融合高度相干滤波器的最小方差超声成像方法 - Google Patents

Abstract

本发明涉及一种融合高度相干滤波器的最小方差超声成像方法，属于超声成像技术领域。该方法包括：对超声阵元接收的回波信号进行延时补偿处理，得到超声回波信号；结合空间平滑和对角加载技术，获得样本协方差矩阵；根据最小方差原理，得到最小方差加权矢量；对协方差矩阵分解，得到信号子空间和噪声子空间；将最小方差加权矢量投影到信号子空间，得到特征空间最小方差加权矢量；根据相干因子得到子阵列回波相干总和与非相干总和，进而得到相干滤波系数；将特征空间最小方差加权矢量与相干滤波系数融合得到高度相干加权矢量，形成自适应波束信号。本发明提高了超声图像的分辨率、对比度以及对噪声的鲁棒性，从而在整体上提高超声成像的质量。

Description

一种融合高度相干滤波器的最小方差超声成像方法

技术领域

本发明属于超声成像技术领域，涉及一种融合高度相干滤波器的最小方差超声成像方法。

背景技术

超声成像中应用最为广泛的，也是最简单的波束形成技术即延时叠加算法(DelayAnd Sum，DAS)，它是根据阵元通道几何位置关系对所接收的回波信号进行延时量的计算，然后对延时后的数据对齐叠加。传统DAS算法复杂度低，成像速度快，但由于其采用固定窗函数加权导致主瓣宽度增加，分辨率较低。

近年来，为了提高波束形成算法的对比度和分辨率，自适应算法得到越来越多的研究。1969年Capon提出的最小方差(Minimum Variance,MV)波束形成算法是目前使用最为广泛的自适应算法。该方法依据保持期望方向增益不变，且使阵列输出能量达到最小的原则，通过动态地计算聚焦延时后的信号加权矢量，再将该矢量与输入信号相乘，提高了图像对比度和分辨率，但该算法的缺点是稳健性远不如传统的延时叠加算法，而且容易使有用信号相消，这在信噪比较低的情况下对图像质量有较大影响。因此，Asl等提出基于特征空间的最小方差(Eigenspace based Minimum Variance,ESBMV)算法，把估计协方差矩阵分解为信号子空间和噪声子空间，并将MV算法得到的权矢量投影到信号子空间得到新的加权矢量。但是，当超声回波信号信噪比较低时，回波中噪声含量高，相干系数低，这将会出现图像整体亮度降低，目标幅值减小等问题。

综上所述，急需发明一种能够在低信噪比条件下提高图像分辨率、对比度，并且保持算法稳健性的波束形成算法，以提高超声成像质量。

发明内容

有鉴于此，本发明的目的在于提供一种融合高度相干滤波器的最小方差超声成像方法，该方法能够在低信噪比条件下提高图像分辨率、对比度以及算法鲁棒性，有效克服了传统自适应波束形成算法在低信噪比条件下，不能显著提高图像对比度和分辨率等问题，从而全面提高了超声图像质量。

为达到上述目的，本发明提供如下技术方案：

一种融合高度相干滤波器的最小方差超声成像方法，包括以下步骤：

S1：对超声阵元接收的回波信号进行延时补偿处理，得到延时聚焦处理后的信号，即超声回波信号；

S2：对回波信号进行空间平滑和对角加载处理，获得样本协方差矩阵；

S3：根据最小方差原理，得到最小方差加权矢量；

S4：对步骤S2中获得的协方差矩阵进行特征分解，得到信号子空间和噪声子空间；

S5：将步骤S3中获得的最小方差加权矢量投影到信号子空间，得到特征空间最小方差加权矢量；

S6：根据相干因子得到子阵列回波相干总和与非相干总和，进而得到相干滤波系数；

S7：将步骤S5中得到的特征空间最小方差加权矢量与步骤S6中得到的相干滤波系数融合得到高度相干加权矢量，得到自适应波束信号。

进一步，所述步骤S1中，对超声阵元接收的回波信号进行延时补偿处理，得到延时聚焦处理后的信号X(t)：

X(t)＝[x₁(t),x₂(t),…,x_M(t)]^T

其中，M表示超声阵列的阵元个数，x_i(t)表示t时刻第i个阵元接收到的信号，t表示采样时间，X(t)里面元素表示在t时刻每个阵元接收的信号，[·]^T表示矩阵转置。

进一步，所述步骤S2中，对回波数据进行空间平滑和对角加载处理，获得样本协方差矩阵，具体包括以下步骤：

S21：将原始换能器的M个阵列划分成大小为L的重叠子阵列，计算各个子阵列的样本协方差矩阵R_l，根据空间平滑技术计算得到估计协方差矩阵为：

其中，L表示子阵列的阵元数，x_l(t)＝[x_l(t),x_l+1(t),…,x_l+L+1(t)]^T表示第l个子阵列接收的回波信号，t表示采样时间，为xl(t)的共轭转置；

S22：对估计协方差矩阵进行对角加载，得到对角加载后的协方差矩阵为：

其中，δ为空间噪声与信号功率之比，取值为0.1，为矩阵主对角元素之和，I为单位矩阵。

进一步，所述步骤S3中，根据最小方差原理，得到最小方差加权矢量，具体包括以下步骤：

S31：对回波信号X(t)添加一个大小为M×1的自适应加权向量ω，得到自适应波束信号y(t)为：

y(t)＝ω^HX(t)

其中，X(t)＝[x₁(t),x₂(t),…,x_l(t),…,x_M-L+1(t)]表示阵元接收的回波信号，M表示阵元总数，L表示子阵列的阵元数，ω^H为ω的共轭转置，t为采样时间；

S32：对于最小方差自适应波束合成，其最优化问题通过下式求解：

minE[(y(t))²]＝minω^HRω,s.t.ω^Hd＝1

其中，E[·]表示期望算子，d＝[1,1,…,1]^T表示长度为M的方向向量，R＝E[X(t)X^H(t)]为阵列数据的协方差矩阵，X^H(t)为X(t)的共轭转置；

S33：根据最小方差原理，得到最小方差加权矢量ω_MV为：

其中，R^-1为R的逆矩阵；

用协方差矩阵代替矩阵R，则最小方差加权矢量ω_MV变为：

其中，为的逆矩阵。

进一步，所述步骤S4中，对步骤S2中获得的协方差矩阵进行特征分解，得到信号子空间和噪声子空间为：

其中，L为子阵列的阵元数，λ_i为的特征值，且有λ₁＞λ₂＞…＞λ_L，v_i为λ_i对应的特征向量，为vi的共轭转置；Λ_s为较大特征值组成的对角矩阵，Λ_n为较小特征值组成的对角矩阵；E_s为较大特征值对应特征向量，E_n为较小特征值对应特征向量，E_s ^H和E_n ^H分别为E_s和E_n的共轭转置；R_s和R_n分别为信号子空间和噪声子空间。

进一步，所述步骤S5中，将步骤S3中获得的最小方差加权矢量投影到信号子空间，得到特征空间最小方差加权矢量为：

其中，ω_MV为最小方差自适应波束加权矢量，ω_ESBMV为特征空间最小方差加权矢量。

进一步，所述步骤S6中，根据相干因子得到子阵列回波相干总和与非相干总和，进而得到相干滤波系数，具体包括以下步骤：

S61：相干因子CF为回波信号相干总和CS与非相干总和ICS之比，表达式为：

S62：子阵列相干总和与非相干总和为所有子阵列上的平均值，通过相干因子CF表达式求得子阵列相干总和CS_sub与非相干总和ICS_sub分别为：

其中，为特征空间最小方差加权矢量ω_ESBMV的共轭转置，x_k(t)为t时刻第k个子阵接收的回波信号，为x_k(t)的共轭转置；

S63：根据子阵列相干总和与非相干总和求得相干滤波系数Q_Hcf为：

其中，γ为基于相干性的噪声加权系数，此处选取γ＝1/L，L为子阵列个数。

进一步，所述步骤S7中，将特征空间最小方差加权矢量与相干滤波系数融合得到高度相干加权矢量，得到自适应波束信号，具体包括以下步骤：

S71：将特征空间最小方差加权矢量ω_ESBMV与相干滤波系数Q_Hcf融合得到高度相干加权矢量ω_ESBMV-Hcf为：

ω_ESBMV-Hcf＝Q_Hcfω_ESBMV

S72：使用高度相干加权矢量ω_ESBMV-Hcf对采样信号进行加权求和，得到自适应波束信号：

其中，y(t)表示计算得到的自适应波束信号，表示ω_ESBMV-Hcf的共轭转置，x_l(t)表示第l个子阵的输出向量；然后根据自适应波束信号进行超声成像。

本发明的有益效果在于：本发明首先对协方差矩阵进行特征分解，得到信号子空间，再将最小方差算法得到的权矢量投影到信号子空间中提高成像的对比度，然后基于信号相干性设计滤波器，使得算法对噪声的鲁棒性进一步增加。因此，本发明所提算法在低信噪比条件下改善图像对比度、分辨率以及算法稳健性方面有较大提高，克服了传统自适应算法在信噪比较低条件下不能显著提高图像对比度和分辨率等问题。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为本发明所述方法的流程图；

图2为4种算法点目标成像结果；

图3为4种算法60mm焦点处横向分辨率曲线；

图4为加噪声后4种算法点目标成像结果；

图5为4种算法吸声斑成像结果；

图6为4种算法geabr_0数据成像结果；

图7为geabr_0实验15mm处散射点横向截面图。

具体实施方式

下面将结合附图，对本发明的优选实施例进行详细的描述。

图1为本发明的算法流程图，如图1所示，本发明提供一种融合高度相干滤波器的最小方差超声成像方法，包括以下步骤：

步骤S1：对超声阵元接收的回波信号进行延时补偿处理，以获得超声回波数据，得到延时聚焦处理之后的信号X(t)，X(t)表示为：

X(t)＝[x₁(t),x₂(t),…,x_M(t)]^T

其中，M表示超声阵列的阵元个数，x_i(t)表示t时刻第i个阵元接收到的信号，t表示采样时间，X(t)里面元素表示在t时刻每个阵元接收的信号，[·]^T表示矩阵转置。

步骤S2：对回波数据进行空间平滑和对角加载处理，获得样本协方差矩阵，具体包括以下步骤：

S21：将原始换能器的M个阵列划分成大小为L的重叠子阵列，计算各个子阵列的样本协方差矩阵R_l，根据空间平滑技术计算得到估计协方差矩阵为：

其中，L表示子阵列的阵元数，x_l(t)＝[x_l(t),x_l+1(t),…,x_l+L+1(t)]^T表示第l个子阵列接收的回波信号，t表示采样时间，为x_l(t)的共轭转置；

S22：对估计协方差矩阵进行对角加载，得到对角加载后的协方差矩阵为：

其中，δ为空间噪声与信号功率之比，取值为0.1，为矩阵主对角元素之和，I为单位矩阵。

步骤S3：根据最小方差原理，得到最小方差加权矢量，具体步骤如下：

S31：对回波信号X(t)添加一个大小为M×1的自适应加权向量ω，得到自适应波束信号y(t)为：

y(t)＝ω^HX(t)

其中，X(t)＝[x₁(t),x₂(t),…,x_l(t),…,x_M-L+1(t)]表示阵元接收的回波信号，M表示阵元总数，L表示子阵列的阵元数，ω^H为ω的共轭转置，t为采样时间；

S32：对于最小方差自适应波束合成，其最优化问题通过下式求解：

minE[(y(t))²]＝minω^HRω,s.t.ω^Hd＝1

其中，E[·]表示期望算子，d＝[1,1,…,1]^T表示长度为M的方向向量，R＝E[X(t)X^H(t)]为阵列数据的协方差矩阵，X^H(t)为X(t)的共轭转置；

S33：根据最小方差原理，得到最小方差加权矢量ω_MV为：

其中，R^-1为R的逆矩阵；

用协方差矩阵代替矩阵R，则最小方差加权矢量ω_MV变为：

其中，为的逆矩阵。

步骤S4：对协方差矩阵进行特征分解，得到信号子空间和噪声子空间：

对步骤S2中获得的协方差矩阵进行特征分解，得到信号子空间和噪声子空间为：

其中，L为子阵列的阵元数，λ_i为的特征值，且有λ₁＞λ₂＞…＞λ_L，v_i为λ_i对应的特征向量，为v_i的共轭转置；Λ_s为较大特征值组成的对角矩阵，Λ_n为较小特征值组成的对角矩阵；E_s为较大特征值对应特征向量，E_n为较小特征值对应特征向量，E_s ^H和E_n ^H分别为E_s和E_n的共轭转置；R_s和R_n分别为信号子空间和噪声子空间。

步骤S5：将获得的最小方差加权矢量投影到信号子空间，得到特征空间最小方差加权矢量：

其中，ω_MV为最小方差自适应波束加权矢量，ω_ESBMV为特征空间最小方差加权矢量。

步骤S6：根据相干因子得到子阵列回波相干总和与非相干总和，进而得到滤波系数，具体步骤如下：

S61：相干因子CF为回波信号相干总和CS与非相干总和ICS之比，表达式为：

S62：子阵列相干总和与非相干总和为所有子阵列上的平均值，通过相干因子CF表达式求得子阵列相干总和CS_sub与非相干总和ICS_sub分别为：

其中，为特征空间最小方差加权矢量ω_ESBMV的共轭转置，x_k(t)为t时刻第k个子阵接收的回波信号，为x_k(t)的共轭转置；

S63：根据子阵列相干总和与非相干总和求得相干滤波系数Q_Hcf为：

其中，γ为基于相干性的噪声加权系数，此处选取γ＝1/L，L为子阵列个数。

步骤S7：将特征空间最小方差加权矢量与相干滤波系数融合得到高度相干加权矢量，得到自适应波束信号，具体步骤如下：

S71：将特征空间最小方差加权矢量ω_ESBMV与相干滤波系数Q_Hcf融合得到高度相干加权矢量ω_ESBMV-Hcf为：

ω_ESBMV-Hcf＝Q_Hcfω_ESBMV

S72：使用高度相干加权矢量ω_ESBMV-Hcf对采样信号进行加权求和，得到自适应波束信号：

其中，y(t)表示计算得到的自适应波束信号，表示ω_ESBMV-Hcf的共轭转置，x_l(t)表示第l个子阵的输出向量；然后根据自适应波束信号进行超声成像。

Field II是丹麦理工大学基于声学原理开发的一款超声实验仿真平台，其在理论研究上获得了广泛的认可和使用。为验证所提算法的有效性，利用Field II对超声成像中常用的点散射目标和吸声斑目标进行成像并利用实际实验数据进行成像对比实验。在点目标仿真实验中，设置设定了15个目标点，轴向距离分布在30mm～80mm，每隔5mm设定目标点，其中40mm和60mm处分布3个目标点其余位置各有1个目标点，采用发射定点聚焦和接收动态聚焦方式，发射焦点固定在60mm处，设置图像的成像动态范围为60dB。同时，设一中心在25mm，半径为3mm的圆形区域吸声斑，外部随机分布着100000个散射点，并设定成像动态范围为80dB。实验所采用的阵元中心频率为3.33MHz，阵元数目为64个，间距为0.2413mm，采样频率为17.76MHz，声速为1500m/s，设成像动态范围为60dB。对上述三个实验目标采用延时叠加算法(DAS)，最小方差算法(MV)，特征空间最小方差算法(ESBMV)，本发明所述的融合高度相干滤波器的特征空间最小方差算法(ESBMV-Hcf)以进行对比成像实验。

图2给出了4种算法点目标成像结果。从图2中可以看出DAS算法成像质量最差，分辨率最低，相比于其他3种算法横向伪影最多，并排三个散射点已经相互干扰难以区分。MV算法较DAS算法旁瓣有所降低，在焦点处散射点已基本能够区分，但在其他深度处横向伪影仍较多，分辨率有待提高。ESBMV算法在整个深度范围内可以明显分辨出目标点。其中ESBMV-Hcf算法成像质量最优，点目标主瓣宽度最小。

图3给出了4种算法60mm焦点处横向分辨率曲线图。从图3中可以看出，DAS算法成像分辨率最差，主瓣宽度最宽且旁瓣等级最高。MV算法较DAS算法成像有所提高，其主瓣宽度和旁瓣等级都有所改善。ESBMV算法与DAS相比，主瓣宽度以及旁瓣等级已经改善明显，其中ESBMV-Hcf算法主瓣最窄，其-6dB主瓣宽度较MV算法减小65.6％，旁瓣等级最低，图像对比度最高。

图4给出了增加20dB高斯白噪声后4种算法点目标成像结果。从图4中可以看出，加入较强的噪声后，前三种算法成像背景区域明显观察出噪声白斑的存在，ESBMV-Hcf算法成像效果最好，背景区域白噪声明显减少，可以清楚地观测出目标点，对噪声鲁棒性最好。

图5给出4种算法吸声斑成像结果，表1给出4种算法对比度。从图5中可以看出，DAS算法相比于其他算法成像效果最差。ESBMV算法成像效果较DAS和MV算法有所提升，图像分辨率和对比度增加。ESBMV-Hcf算法效果最好，对比度最高，算法旁瓣抑制能力最强。由表1可见，DAS算法对比度最低，仅为21.43dB，背景方差最小，算法稳健性最好。MV算法提高了中心暗斑平均功率，但其外部平均功率也同时提高，对比度较DAS算法上升约2dB。ESBMV算法中心暗斑及背景功率分别在MV基础上有所提高。其中，ESBMV-Hcf算法对比度上升最多，对比度较DAS，MV，ESBMV，ESBMV-CF算法分别提高了9.09dB、8.01dB、2.08dB，且背景区域方差低于ESBMV。

表1 4种算法对比度

图6给出了5种算法geabr_0数据成像结果。图7给出geabr_0实验15mm处散射点横向截面图。从图6中可以看出，传统DAS算法成像效果最差，近场点目标受背景噪声干扰最为严重，采用自适应算法成像较DAS算法都好，其图像分辨率和对比度都有所改善，其中ESBMV-Hcf算法分辨率最高，对比度改善明显。从图7可以看出，ESBMV-Hcf算法分辨率和对比度最高，其主瓣宽度最窄，最大旁瓣幅值最低。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。

Claims (8)

1.一种融合高度相干滤波器的最小方差超声成像方法，其特征在于，该方法包括以下步骤：

S1：对超声阵元接收的回波信号进行延时补偿处理，得到延时聚焦处理后的信号，即超声回波信号；

S2：对回波信号进行空间平滑和对角加载处理，获得样本协方差矩阵；

S3：根据最小方差原理，得到最小方差加权矢量；

S4：对步骤S2中获得的协方差矩阵进行特征分解，得到信号子空间和噪声子空间；

S5：将步骤S3中获得的最小方差加权矢量投影到信号子空间，得到特征空间最小方差加权矢量；

S6：根据相干因子得到子阵列回波相干总和与非相干总和，进而得到相干滤波系数；

S7：将步骤S5中得到的特征空间最小方差加权矢量与步骤S6中得到的相干滤波系数融合得到高度相干加权矢量，进而得到自适应波束信号。

2.根据权利要求1所述的一种融合高度相干滤波器的最小方差超声成像方法，其特征在于，所述步骤S1中，对超声阵元接收的回波信号进行延时补偿处理，得到延时聚焦处理后的信号X(t)：

X(t)＝[x₁(t),x₂(t),…,x_M(t)]^T

其中，M表示超声阵列的阵元个数，x_i(t)表示t时刻第i个阵元接收到的信号，t表示采样时间，X(t)里面元素表示在t时刻每个阵元接收的信号，[·]^T表示矩阵转置。

3.根据权利要求2所述的一种融合高度相干滤波器的最小方差超声成像方法，其特征在于，所述步骤S2中，对回波数据进行空间平滑和对角加载处理，获得样本协方差矩阵，具体包括以下步骤：

S21：将原始换能器的M个阵列划分成大小为L的重叠子阵列，计算各个子阵列的样本协方差矩阵R_l，根据空间平滑技术计算得到估计协方差矩阵为：

其中，L表示子阵列的阵元数，x_l(t)＝[x_l(t),x_l+1(t),…,x_l+L+1(t)]^T表示第l个子阵列接收的回波信号，t表示采样时间，为x_l(t)的共轭转置；

S22：对估计协方差矩阵进行对角加载，得到对角加载后的协方差矩阵为：

其中，δ为空间噪声与信号功率之比，为矩阵主对角元素之和，I为单位矩阵。

4.根据权利要求3所述的一种融合高度相干滤波器的最小方差超声成像方法，其特征在于，所述步骤S3中，根据最小方差原理，得到最小方差加权矢量，具体包括以下步骤：

S31：对回波信号X(t)添加一个大小为M×1的自适应加权向量ω，得到自适应波束信号y(t)为：

y(t)＝ω^HX(t)

其中，X(t)＝[x₁(t),x₂(t),…,x_M(t)]^T表示阵元接收的回波信号，M表示阵元总数，ω^H为ω的共轭转置，t为采样时间；

S32：对于最小方差自适应波束合成，其最优化问题通过下式求解：

minE[(y(t))²]＝minω^HRω,s.t.ω^Hd＝1

其中，E[·]表示期望算子，d＝[1,1,…,1]^T表示长度为M的方向向量，R＝E[X(t)X^H(t)]为阵列数据的协方差矩阵，X^H(t)为X(t)的共轭转置；

S33：根据最小方差原理，得到最小方差加权矢量ω_MV为：

其中，R-¹为R的逆矩阵；

用协方差矩阵代替矩阵R，则最小方差加权矢量ω_MV变为：

其中，为的逆矩阵。

5.根据权利要求4所述的一种融合高度相干滤波器的最小方差超声成像方法，其特征在于，所述步骤S4中，对步骤S2中获得的协方差矩阵进行特征分解，得到信号子空间和噪声子空间为：

其中，L为子阵列的阵元数，λ_i为的特征值，且有λ₁＞λ₂＞…＞λ_L，v_i为λ_i对应的特征向量，为vi的共轭转置；Λ_s为较大特征值组成的对角矩阵，Λ_n为较小特征值组成的对角矩阵；E_s为较大特征值对应特征向量，E_n为较小特征值对应特征向量，E_s ^H和E_n ^H分别为E_s和E_n的共轭转置；R_s和R_n分别为信号子空间和噪声子空间。

6.根据权利要求5所述的一种融合高度相干滤波器的最小方差超声成像方法，其特征在于，所述步骤S5中，将步骤S3中获得的最小方差加权矢量投影到信号子空间，得到特征空间最小方差加权矢量为：

其中，ω_MV为最小方差自适应波束加权矢量，ω_ESBMV为特征空间最小方差加权矢量。

7.根据权利要求6所述的一种融合高度相干滤波器的最小方差超声成像方法，其特征在于，所述步骤S6中，根据相干因子得到子阵列回波相干总和与非相干总和，进而得到相干滤波系数，具体包括以下步骤：

S61：相干因子CF为回波信号相干总和CS与非相干总和ICS之比，表达式为：

S62：子阵列相干总和与非相干总和为所有子阵列上的平均值，通过相干因子CF表达式求得子阵列相干总和CS_sub与非相干总和ICS_sub分别为：

其中，为特征空间最小方差加权矢量ω_ESBMV的共轭转置，x_k(t)为t时刻第k个子阵接收的回波信号，为x_k(t)的共轭转置；

S63：根据子阵列相干总和与非相干总和求得相干滤波系数Q_Hcf为：

其中，γ为基于相干性的噪声加权系数。

8.根据权利要求7所述的一种融合高度相干滤波器的最小方差超声成像方法，其特征在于，所述步骤S7中，将特征空间最小方差加权矢量与相干滤波系数融合得到高度相干加权矢量，得到自适应波束信号，具体包括以下步骤：

S71：将特征空间最小方差加权矢量ω_ESBMV与相干滤波系数Q_Hcf融合得到高度相干加权矢量ω_ESBMV-Hcf为：

ω_ESBMV-Hcf＝Q_Hcfω_ESBMV

S72：使用高度相干加权矢量ω_ESBMV-Hcf对采样信号进行加权求和，得到自适应波束信号：

其中，y(t)表示计算得到的自适应波束信号，表示ω_ESBMV-Hcf的共轭转置，x_l(t)表示第l个子阵的输出向量；然后根据自适应波束信号进行超声成像。