CN101822550B - 超声彩色血流成像中基于动态区域划分的杂波抑制方法 - Google Patents

Abstract

本发明属于超声彩色血流成像技术领域，具体涉及一种基于动态区域划分的杂波抑制方法。本方法先根据回波信号的幅度对血流与组织区域做出初步划分，再根据杂波运动的非平稳性对组织区域作进一步调整，最后对各区域用特征向量滤波器进行杂波抑制。本方法在抑制非平稳杂波的同时，能较好地保持血流流速剖面的完整性，解决了传统均匀分段特征向量滤波器选取区间长度时的矛盾，是彩色血流成像中一种有效的杂波抑制方法。

Description

超声彩色血流成像中基于动态区域划分的杂波抑制方法

技术领域

本发明属于超声彩色血流成像技术领域，具体为一种基于动态区域划分的杂波抑制方法。

背景技术

超声彩色血流成像技术(CFI)能够显示待测剖面上的二维血流速度分布，是诊断血管疾病的重要依据。CFI首先利用超声换能器沿某个方向重复发射M次脉冲(间隔为T_prf)。接收到的M个回波即携带了该扫描线方向上各深度处的运动信息。对回波信号进行处理后可得到该方向上的血流速度分布。最后将各条扫描线上的速度估计排列成矩阵形式，以彩色编码显示，即得到了整个待测剖面上的流速分布。

在接收的回波信号中往往包含了来自管壁搏动和慢速组织运动的干扰(称为杂波)，其功率通常比血流信号高出40到100dB。这将给最终血流速度的正确估计带来很大的困难，所以必须采用高性能的滤波器来抑制杂波的影响。

传统多普勒系统中常采用高通滤波器(HPF)来抑制杂波。现代CFI系统为保证图像帧率等原因，重复脉冲次数(即数据的有效长度)将受到限制。所以不得不采用相对低阶的HPF，因此滤波性能也就难以得到保障。

特征向量滤波器，也称为奇异值分解法，是一种适用于CFI的自适应杂波抑制方法。该方法要求杂波运动具有平稳性。实际应用中往往通过对回波数据均匀分段来近似满足这一条件。

发明内容

本发明的目的是针对特征向量滤波器，提供一种动态区域划分的方法，以取代传统的均匀划分。尤其涉及超声彩色血流成像中基于动态区域划分的杂波抑制方法。

本发明的目的通过下述方法和步骤实现：

先根据初始段长2K₀+1对回波信号均匀划分，计算各区间内的平均能量并与总平均能量比较，得到血流与组织区域的初步划分。对位于组织区域的信号，计算得到指标矢量δ来表征各区间内杂波运动的非平稳程度。利用δ的最值和K₀得到δ^(l)与对应区间长度的函数关系，并对段长进行调整，得到最终动态区域划分后的结果。

下面对各步骤作进一步具体描述。

设N和M分别为某一扫描线上的纵向采样点数和重复发射次数。若第i次回波信号正交解调后的数据为[X_i，1 X_i，2...X_i，l]，则由同一扫描线上M个回波构成的复数回波矩阵可表示为：

$X=\left[\begin{array}{cccc}{x}_{\mathrm{1,1}}& x_{\mathrm{1,2}}& \mathrm{\Λ}& x_{1,N}\\ x_{\mathrm{2,1}}& x_{\mathrm{2,2}}& \mathrm{\Λ}& x_{2,N}\\ M& M& O& M\\ x_{M,1}& x_{M,2}& \mathrm{\Λ}& x_{M,N}\end{array}\right]---\left(1\right)$

将X沿深度方向均匀分段，区间长度初始为2K₀+1，共可分为

段：

对于 $x_{i,j}^{\left(l\right)}\∈X^{\left(l\right)}(i=1,...,M;j=1,...,{2K}_0+1),$ 存在关系 $x_{i,j}^{\left(l\right)}=x_{i,(l-1)({2K}_0+1)+j}.$

计算整个回波信号的平均能量E_mean：

$E_{\mathrm{mean}}=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\left|x_{i,j}\right|}^2---\left(3\right)$

对于每个区间l＝1，2，...，L，计算区间内信号平均能量E^(l)：

$E^{\left(l\right)}=\frac{1}{M({2K}_0+1)}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{{2K}_0+1}{\mathrm{\Σ}}}{\left|x_{i,j}^{\left(l\right)}\right|}^2---\left(4\right)$

由于血流区域的回波幅度远小于组织区域的回波幅度，所以可以利用阈值来区分血流与组织区域。平均能量大于该阈值的区间归入组织区域，否则归入血流区域：

其中β为常数，满足0＜β＜1，这里取0.1。

接着计算组织区域的非平稳性指标矢量δ。对于X^(l)，若之前被归入血流区域，则不做处理，l←l+1；若被归入组织区域，则进一步计算δ，过程如下：

令X₁ ^(l)等于X^(l)去掉首行，X₂ ^(l)等于X^(l)去掉尾行，即：

计算互相关矩阵R^(l)：

$R^{\left(l\right)}=X_1^{\left(l\right)H}{X}_2^{\left(l\right)}---\left(7\right)$

其中R^(l)为(2K₀+1)×(2K₀+1)的方阵。对每个 $r_{i,j}^{\left(l\right)}\∈R^{\left(l\right)},$ 计算其幅角^[7]：

${\mathrm{\φ}}_{i,j}^{\left(l\right)}=\arctan \left(\frac{\mathrm{Im}\left(r_{i,j}^{\left(l\right)}\right)}{\mathrm{Re}\left(r_{i,j}^{\left(l\right)}\right)}\right),i,j=1,...,{2K}_0+1---\left(8\right)$

计算非平稳性指标δ^(l)：

${\mathrm{\δ}}^{\left(l\right)}=\sqrt{\mathrm{var}\left(\mathrm{diag}\left({\mathrm{\Φ}}^{\left(l\right)}\right)\right)}---\left(9\right)$

其中Ф^(l)是由φ_i，j ^(l)构成的方阵。设共有L_c(L_c＜L)个组织区域，即可构成指标矢量：

$\mathrm{\δ}=\left[\begin{array}{ccccc}{\mathrm{\δ}}^{\left(1\right)}& \mathrm{\Λ}& {\mathrm{\δ}}^{\left(l\right)}& \mathrm{\Λ}& {\mathrm{\δ}}^{\left(L_c\right)}\end{array}\right].---\left(10\right)$

式(8)与自相关法速度估计公式仅相差一个系数，方阵Ф^(l)对角线上的元素就对应了组织区域l内各点的杂波运动速度。式(9)用标准差的形式来衡量某组织区域内各点速度的变化剧烈程度，即对应了杂波运动的非平稳性。

矢量δ的元素δ^(l)表征了第l个组织区间内的非平稳程度：δ^(l)越小，杂波运动越平稳，区间长度可取大些；反之，δ^(l)越大，杂波运动非平稳性越明显，区间长度应相对小。根据上述原则即可建立非平稳指标δ^(l)与区间长度参数K^(l)的对应关系：

设δ^(l)与区间长度参数K^(l)近似为一线性关系(如图1所示)。令δ_max＝max(δ)处K^(l)＝K₀；δ_min＝min(δ)处K^(l)＝2K₀，这样就保证了K^(l)与δ^(l)成反比关系，且数据重叠部分不会超过2K₀，即相邻区间的中心不会重叠。

由两点坐标P₁(δ_max，K₀)和P₂(δ_min，2K₀)得直线斜率：

$k=\frac{K_0}{{\mathrm{\&delta;}}_{\min }-{\mathrm{\&delta;}}_{\max }}<0---\left(11\right)$

所以非平稳指标δ^(l)与区间长度参数K^(l)的函数对应关系为：

$K^{\left(l\right)}=F\left({\mathrm{\&delta;}}^{\left(l\right)}\right)=K_0+k({\mathrm{\&delta;}}^{\left(l\right)}-{\mathrm{\&delta;}}_{\max })$

$=K_0(1+\frac{{\mathrm{\&delta;}}^{\left(l\right)}-{\mathrm{\&delta;}}_{\max }}{{\mathrm{\&delta;}}_{\min }-{\mathrm{\&delta;}}_{\max }})---\left(12\right)$

由式(12)可见，F(δ^(l))为减函数，非平稳性越明显，δ^(l)越大，K^(l)就越小；反之平稳性越强，K^(l)就越大。这样就得到了最终动态划分后的各区间长度，它将随着杂波的运动情况而改变，且首尾可以部分重叠。

综上所述，本发明提出的动态区域划分算法的基本流程可以概括为：首先根据初始段长2K₀+1对回波信号均匀划分，由式(3)和(4)计算各区间平均能量和总平均能量，得到血流与组织区域的初步划分；再由式(9)(10)计算指标矢量δ，进一步利用式(12)得到δ与对应区间长度的函数关系，并以此为依据对段长进行调整。最后对动态划分后的各区间数据用特征向量滤波器进行杂波抑制，重叠部分取平均值即可。

在各种非平稳杂波情况下，本方法仍能较好地抑制杂波并保证血流速度剖面的完整性，是一种快速有效的杂波抑制方法。

附图说明

图1、非平稳指标δ^(l)与区间长度参数K^(l)的函数对应关系。

图2、仿真散射子模型。

图3、五种杂波运动模式：静止、匀速、线性增加、非线性变化、不连续变化。

图4、杂波运动模式为静止时的杂波抑制结果。

图5、杂波运动模式为匀速时的杂波抑制结果。

图6、杂波运动模式为线性增加时的杂波抑制结果。

图7、杂波运动模式为非线性变化时的杂波抑制结果。

图8、杂波运动模式为不连续变化的杂波抑制结果。

图9、杂波抑制能力最强的三种方法的血流速度剖面完整度比较。

图10、血流速度剖面完整度最好的三种方法的杂波抑制能力比较。：

图11、人体颈动脉彩色血流成像结果：(a)HPF(b)均匀分段特征向量滤波器(c)HankelSVD法(d)动态区域划分法。

具体实施方案

以下结合具体的实施例，对本发明做进一步的阐述。实施例仅用于对本发明做说明而不是对本发明的限制。

实施例1

在计算机上用MATLAB进行仿真实验，运行环境为Pentium Core Dual1.8GHz。利用Jensen等开发的FieldII软件包仿真不同环境下的回波信号。图2为采用的散射子模型，虚线为声束方向。表1列出了仿真参数设置。

表1 FieldII仿真参数设置

同时，仿真采用了5种杂波运动模式，如图3所示。参与比较的方法有：高通滤波器(HPF)、特征向量滤波器(均匀分段)、HankelSVD法和动态区域划分法。对于均匀分段特征向量滤波器中段长参数K的不同取值如表2所示：

表2 K的不同取值

图4～8列出了5种杂波运动模式下的仿真结果，后三种模式有着较明显的非平稳性，下面以图7为例进行分析。

图7(a)中比较了四种杂波抑制方法：HPF、特征向量滤波器(K₀＝K_max即不分段)、HankelSVD法和本发明提出的方法，同时还给出了标准流速剖面。由于HPF的阶数受到限制(3阶)，阻带衰减有限导致了杂波抑制效果不佳。不分段的特征向量滤波器能使血流区域邻近的杂波得到抑制，而远离中心处却效果不佳，这是由于组织的非平稳运动引起的，不同区域的统计特性差异使抑制效果受到了影响。HankelSVD有着很好的杂波抑制效果，但在对流速剖面影响较大，其原因是HankelSVD法是基于单数据集(single-ensemble)的方法，参与协方差矩阵估计的样本远远少于特征向量滤波器法。而本发明提出的动态区域划分法能够较好地抑制非平稳杂波运动，同时又保证了流速剖面的完整性。

图7(b)中进一步比较了不同段长下(K_small＜K_median＜K_large)，特征向量滤波器的杂波抑制结果。可以看到段长越小，杂波运动能得到更好的抑制，但流速剖面将受到影响；段长越大，流速剖面相对完整，但杂波抑制效果变差。动态区域划分法能够针对不同的杂波运动情况改变段长，从而解决了均匀分段时的矛盾。

以上针对图7的分析可以推广到其余四种杂波运动模式下的结果。

进一步定量分析各种杂波抑制方法：用非血流区域内估计值与理论值的均方误差(MSE)来衡量杂波抑制能力；同时用理论流速剖面与估计流速剖面的相关系数来衡量滤波后速度剖面的完整度。

图9比较了杂波抑制能力最强的三种方法(HankelSVD法、动态区域划分法和K_small均匀分段)的血流速度剖面完整度。横坐标为五种不同的杂波运动模式，纵坐标为剖面完整度(Corr.)。可见动态区域划分法结果最优。图10比较了速度剖面完整度最好的三种方法(K_large、K_median均匀分段、动态区域划分法)的杂波抑制能力(MSE)。综合MSE和Corr.两方面的表现，动态区域划分法具有较明显的优势。

表3比较了四种杂波抑制方法运行速度。动态区域划分法的耗时与HPF和均匀分段特征向量滤波器相近，大大小于HankelSVD法。

表3  四种杂波抑制方法运行速度比较

图11为实际人体颈动脉彩色血流成像结果，图中白线为人为给定的血流区域参考线。为定量比较，计算图11中各血流区域内外信号能量及其比值(BV/MT)，结果如表4所示。可见动态区域划分法的BV/MT结果最大，效果最好。

表4组织区域和血流区域的能量及能量比

由仿真和实验结果可见，动态区域划分法解决了传统均匀分段特征向量滤波器在选取区间长度时的矛盾，在杂波抑制能力和保证流速完整度两方面都表现突出。因此，本发明提出的动态区域划分法是彩色血流成像中一种快速有效的杂波抑制算法。

Claims (1)

1.一种超声彩色血流成像中基于动态区域划分的杂波抑制方法，其特征在于先以初始段长2K₀+1对回波信号均匀划分，获得回波区间；计算各区间内的平均能量并与总平均能量比较，得到血流与组织区域的初步划分；对位于组织区域的信号，通过计算互相关矩阵R^(l)得到指标矢量δ，以该指标矢量δ表征各区间内杂波运动的非平稳程度；利用所述的指标矢量δ的最值和初始值K₀计算得到δ与对应区间长度的函数关系假设为线性关系，并进行区间长度的调整；最终得到动态区域划分后的结果；

所述的血流与组织区域的初步划分通过下述计算公式：

$E_{\mathrm{mean}}=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left|x_{i,j}\right|}^2---\left(1\right)$

其中E_mean为总平均能量，M为同一扫描线上的重复发射次数，N为深度采样点数，x_i，j为第i个脉冲回波中，深度j处的回波信号值；

$E^{\left(l\right)}=\frac{1}{M({2K}_0+1)}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{{2K}_0+1}{\mathrm{\&Sigma;}}}{\left|x_{i,j}^{\left(l\right)}\right|}^2---\left(2\right)$

其中E^(l)为区间l内的回波平均能量，2K₀+1为该区间的长度，根据E^(l)与E_mean的关系对血流和组织区域作初步划分，其中，平均能量大于该阈值的区间归入组织区域，否则归入血流区域：

其中β为常数，满足0＜β＜1；

所述的互相关矩阵R^(l)如下式所示：

$R^{\left(l\right)}=X_1^{\left(l\right)H}{X}_2^{\left(l\right)}---\left(4\right)$

当X^(l)为杂波区域l内回波信号构成的矩阵，则为X^(l)去掉首行，

为X^(l)去掉尾行；

${\mathrm{\&phi;}}_{i,j}^{\left(l\right)}=\arctan \left(\frac{\mathrm{Im}\left(r_{i,j}^{\left(l\right)}\right)}{\mathrm{Re}\left(r_{i,j}^{\left(l\right)}\right)}\right),i,j=1,...,{2K}_0+1---\left(5\right)$

其中 $r_{i,j}^{\left(l\right)}\&Element;R^{\left(l\right)};$

${\mathrm{\&delta;}}^{\left(l\right)}=\sqrt{\mathrm{var}\left(\mathrm{diag}\left({\mathrm{\&Phi;}}^{\left(l\right)}\right)\right)}---\left(6\right)$

其中Ф^(l)是由

构成的矩阵；

当共有L_c。个组织区域，则构成指标矢量δ：

$\mathrm{\&delta;}=[{\mathrm{\&delta;}}^{\left(1\right)}\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&delta;}}^{\left(l\right)}\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&delta;}}^{\left(L_c\right)}];---\left(7\right)$

指标矢量6与区间长度参数K^(l)的关系为：

K^(l)＝F(δ^(l)＝K₀+K₀(δ^(l)-δ_max)/(δ_min-δ_max)    (8)

其中K₀为初始段长参数；δ_min＝min(δ)，δ_max＝max(δ)。

Images