基于Hankel-SVD的非平稳超声血流成像杂波抑制技术研究

王录涛王 微金 钢

①(电子科技大学自动化工程学院 成都 611731)

②(国家空气动力研究与发展中心 绵阳 621000)

基于Hankel-SVD的非平稳超声血流成像杂波抑制技术研究

王录涛*①王 微②金 钢②

①(电子科技大学自动化工程学院 成都 611731)

②(国家空气动力研究与发展中心 绵阳 621000)

有效抑制由血管或血管周围组织时变运动引起的非平稳杂波对于提高诊断超声彩色血流成像中血流动力学参数描述的准确性有着极其重要的意义。该文基于奇异值滤波技术提出一种改进的非平稳杂波自适应抑制方法。该方法逐次利用单个慢时多普勒回波采样矢量构建Hankel矩阵，然后根据奇异值分解后得到的正交Hankel主成份所代表的频域内容，动态选取高阶Hankel主成份重构多普勒血流信号，实现非平稳杂波的有效抑制。为验证算法的有效性，分别对多普勒回波仿真模型合成数据与利用彩色超声设备(Sonix RP)采集的颈动脉血流基带回波信号进行滤波处理，然后采用滞一自相关估计法计算血流平均速度与功率并进行成像。处理结果表明，相对于传统IIR滤波方法与多项式回归滤波技术，利用该文所提算法可对高强度、非平稳杂波进行充分抑制，提高血流估计精度，此外，该算法具有空间自适应性，无需人为设定阈值参数以估计杂波空间维数，与现有基于特征分解的自适应滤波方法相比，可以有效提高组织空间高强度时变运动时血流与组织的区分能力。

彩色血流成像；奇异值分解；自适应杂波抑制；非平稳杂波

1 引言

在超声彩色血流成像技术中，杂波指由血管或血管周围慢动组织反射引起的多普勒回波分量，该分量的强度通常比血流信号强度要高40～100 dB，即使当采样容积位于血管中心位置处，由于采样容积的3维特性、合成声束旁瓣泄露、声波多重反射以及距离模糊等因素的影响，血流回波信号中也包含大量的杂波信号[1,2]。杂波信号的存在导致血流动力学参数估计偏差，严重降低了彩色血流成像系统的血流鉴别能力，因此充分抑制杂波成份对于提高彩色超声成像质量有着至关重要的意义[3,4]。

传统的杂波抑制滤波器是一种高通滤波器，包括有限冲激响应(FIR)滤波器、无限冲激响应(IIR)滤波器与多项式回归滤波器等。在彩色超声血流成像系统中，有限的多普勒回波矢量采样数据导致FIR滤波器的幅频特性无法满足高强度杂波抑制需求。与FIR滤波器相比，相同阶数的IIR滤波器具有更加陡峭的过渡带滚降特性与阻带抑制特性，但当样本数据有限时，IIR滤波器的幅频响应由暂态响应主导[5]。多项式回归型滤波器采用正交多项式对杂波分量进行最小二乘拟合以消除杂波成份的影响，为杂波滤波器的设计提供了新的途径[6]。当杂波为平稳信号时，传统静态杂波滤波器可获得良好的杂波抑制效果，但当组织存在加速运动时，不充分的杂波抑制将导致血流参数估计产生较大偏差[7]。

针对杂波信号的非平稳特性，文献[8]提出了基于特征值分解的杂波滤波器(Eigen filter)设计方法。该方法采用对回波多普勒信号自相关矩阵进行特征分解获得的特征矢量构造杂波滤波器，将特征矢量作为回归滤波器的正交基矢量可使得期望信号的拟合具有最小均方误差性质，但特征值分解滤波器假设整个观察范围内的多普勒矢量信号是平稳的，不具有空间自适应性，此外，杂波空间维数需要根据经验值选择，不恰当的选择将导致血流参数估计严重偏离真实值[9−11]。

作为一种数据特征提取的有效方法，奇异值分解技术可以将数据主导特征用少数分解系数集中体现，具有良好的稳定性和不变性，近年来在特征提取、信号去噪等领域的应用引起了广泛重视[12,13]。本文介绍了一种基于Hankel奇异值分解(Hankel-SVD)的非平稳杂波抑制技术。该技术采用单个慢时回波多普勒矢量构建Hankel矩阵，具有空间自适应性；根据Hankel主成份代表的频率内容自适应选择杂波主成份阶数，通过摒弃杂波Hankel主成分，保留血流Hankel主成份并重构血流多普勒信号，实现非平稳杂波的有效抑制。最后采用基于多普勒回波仿真模型合成数据的血流参数估计结果与人体颈动脉血流成像实例，对本文所提算法针对高强度、非平稳杂波抑制的有效性进行了验证。

2 Hankel-SVD滤波

超声血流成像技术利用超声声束对待成像区域进行扫描，在每一个扫描方向形成有限个发射声束，通过接收血红细胞散射的多普勒信号估计血流动力学参数。Hankel-SVD滤波器采用KL(Karhunen-Loeve)变换将解调后的回波多普勒矢量信号表示为一组正交基矢量的线性叠加，通过保留与血流信号相对应的正交成份并进行信号重构，可有效降低杂波对血流参数估计的影响[11]。对于包含N个采样值的慢时回波多普勒矢量x，其Hankel-SVD逼近可表示为

其中P(P＜N)为用于逼近x的正交基矢量的阶数，T表示转置，γk,ϕk分别为k阶展开系数与正交基矢量，满足下述正交条件：

为计算式(1)中的正交成份，将多普勒矢量x分割成部分重叠的数据段并重新组织以构造Hankel矩阵：

式(3)的奇异值分解可表示为

其中，Ak表示秩为1的第k阶Hankel成份，σk,uk, vk分别为奇异值、左奇异矢量(N维)、右奇异矢量(N−P+1维)。式(4)中的P个奇异值满足σ1＞σ2＞…＞σP，因此Ak又可称为第k阶Hankel主成份。

超声血流多普勒回波信号由杂波、血流信号与噪声3种成份构成，其中杂波成份是由心壁、瓣膜、血管壁及静止的或慢速运动的肌肉引起的，平均频率较低但强度较高，在回波信号中占主导地位，而血流信号是由血红细胞散射的多普勒信号构成，强度较低，但频率成份高于杂波信号[1]。文献[6]对超声多普勒回波自相关矩阵进行特征值分解，得到了超声多普勒信号特征值谱，并结合超声回波多普勒信号特性对特征值谱分布与回波信号不同成份的能量、频率特性构成之间的关系进行了研究，认为与包含低频成份的特征矢量相对应的较大特征值反映了杂波成份构成，与包含高频成份的特征矢量相对

其中Kc为杂波Hankel主成份阶数。

分析Hankel矩阵的构造过程可知，在构造Hankel矩阵A时，输入信号矢量中的每个元素被沿反对角线方向重复配置，虽然经奇异值分解，反对角线方向的元素不再相等，但仍可沿Ak反对角线方向对代表同一输入数据的所有元素求平均值来重建血流多普勒回波信号所需的正交基矢量γkϕk。应的数值较小的高阶奇异值则集中体现了血流多普勒成份特征，而噪声信号能量则均匀分布在所有奇异值上。因此，当杂波分量包含的Hankel主成份阶数已知时，可由Ak计算得到的正交基矢量γkϕk，进而按照式(1)实现血流多普勒信号的准确重构：

3 杂波Hankel主成份阶数估计

为充分抑制杂波成份，必须准确判断Hankel主成份Ak是否为杂波成份。基于特征值分解的杂波滤波器设计方法根据杂波成份与血流成份对奇异值谱的贡献不同，通过将分解获得奇异值与预先设置的经验阈值比较，估计杂波Hankel主成份阶数[11]。当杂波成份在回波信号中处于主导时，可实现杂波子空间维数的准确估计，但对于高频超声血流成像，由于血流信号的Rayleigh散射特性，杂波与血流信号的功率比值(CBR)接近0 dB，采用阈值法将导致杂波子空间维数估计产生较大偏差[15]。此外，组织空间生理特性差异与组织运动的空间、时间不确定性，导致回波信号随空间、时间变化存在巨大差异，因此很难通过一个特定的阈值准确实现杂波信号与血流多普勒信号的分离。

由超声多普勒回波成份构成与奇异值分解原理可知，Hankel-SVD的正交基矢量包含了多普勒回波不同成份的频率信息，杂波成份通常为低频分量，而血流信号则由高频多普勒成份构成[6,16]。因此，可根据ϕk所反映的频率信息判定Hankel主成份Ak是否包含于杂波子空间，进而估计杂波Hankel主成份的阶数。首先利用ϕk计算第k阶平均多普勒频率：

其中，TPRI为发射脉冲重复间隔。fD,1为一阶Hankel主成份，可用于区分当前多普勒矢量是否由杂波成份所主导，如为杂波成份主导，则可以根据当前多普勒矢量的瞬时杂波频率估计杂波带宽：其中，fi,ct为瞬时杂波频率，可由当前采样的复多普勒矢量计算得出：

其中I与Q分别为复多普勒信号采样值的实部与虚部。根据fD,1与fth,ct可确定杂波频率范围，进而对Hankel主成份Ak的性质进行判决：

当fD,k(k＞1)满足式(9)时，说明Ak仍为杂波成份，此时，需增加k值进行下一步搜索，即将fD,k+1代入式(9)，验证fD,k+1是否位于杂波频率范围内，杂波Hankel主成份阶数为使式(9)不成立的最小k值，即Kc=k。增加k值使得杂波成份逐次从回波多普勒信号中抽取出来，当k＞Kc时，则杂波成份被有效剥离，Ak以血流成份为主导，根据式(5)可准确重构血流多普勒信号，而累加大于Kc的奇异值并进行归一化，则可得到滤波后的血流平均功率：

4 血流参数估计仿真

为了对本文设计的滤波器的杂波抑制性能进行量化评估，采用文献[16,17]提出的回波合成模型产生仿真多普勒回波信号，仿真参数如表1所示，滤波处理后的数据经滞一自相关估计处理以估计血流平均多普勒频率与功率。为分析不同流速时的滤波器杂波抑制性能，设计了20组血流信号，其平均多普勒频率在0～0.5发射脉冲重复频率(PRF)范围内均匀分布。

表1 壁滤波器性能比较仿真参数取值表

图1给出了200次独立仿真实验得到的血流平均多普勒频率与功率估计结果。在图1中，IIR-Prj为六阶投影初始化Chebyshev IIR滤波器，通带归一化截止频率为0.15，为保证相近的阻带抑制性能，Reg-Pol多项式回归滤波器采用五阶Legendre 多项式重构血流信号，Eigen滤波器杂波空间维数设为5。

图1(a)给出了估计平均多普勒频率比较结果，虚线表示理想估计平均多普勒频率。从图1(a)中可以看出采用Reg-Pol与Eigen得到的平均多普勒频率估计精度优于IIR-Prj，这是因为采用正交基函数重构血流信号可有效提升滤波器的杂波抑制性能，有效降低因杂波残留导致的频率估计偏差。而对比Reg-Pol与Eigen则可以发现，Eigen可以获得更小的逼近误差，这是由于Eigen将回波信号作为随机过程处理，根据杂波的统计特性构造杂波子空间的基矢量，当杂波子空间维数小于信号空间维数时，可使得杂波逼近最小均方误差最小[6]。相对于其它3种被计较滤波算法，当输入血流平均频率从100 Hz增加到200 Hz时， Hankel-SVD可以显著提高血流速度估计精度。

图1(b)示出了平均功率估计比较结果，其中IIR-Prj滤波处理后的数据仍保留较多杂波成份，导致功率估计高于其它3种杂波滤波器。Reg-Pol滤波器具有最强的阻带抑制性能，但过渡带较宽，而在更宽的频率范围内，Hankel-SVD功率估计值更接近于理想值。

5 颈动脉血流成像比较

为验证本文算法在活体诊断应用中的有效性，利用Sonix RP (Ultrasonix Inc., Vancouver, Cannada)超声成像系统采集人体颈动脉基带多普勒回波数据，探头中心频率为6.6 MHz，脉冲重复频率为2.5 kHz，发射脉冲重复次数为10，系统采样频率为40 MHz，轴向采样750点，侧向扫描50线。对杂波抑制后的采样数据，采用滞一自相关估计法计算每个空间采样点的血流速度与功率并进行成像。

图2 给出了一个完整心动周期50帧数据中心脏收缩期的一帧数据经4种杂波滤波器处理后的，未设置速率门限与功率门限的彩色多普勒血流成像结果。 其中，IIR-Prj为二阶、归一化截止频率为0.05的投影初始化IIR滤波器，Reg-Pol为二阶Legendre多项式回归滤波器，构建Eigen滤波器自相关矩阵的慢时多普勒矢量数为750，杂波空间维数固定为2。

在图2中，不同的颜色代表不同的速度估计结果。在远离血管区域杂波成份被充分抑制，噪声成份使得速度估计趋向于系统抗混叠速度，在图像中表现为红、蓝交错分布。绿色区域表示未能被充分抑制的杂波成份，在血管壁、血管附近组织以及靠近内血管壁的血流区域均有所分布。对比图2可知，图2(d) 中杂波残留区域最小，图2(c)中最大，图2 (b)优于图2(a)。这是因为Reg-Pol具有更加陡峭的过渡带，因而杂波性能优于IIR-Prj，但在血管壁上缘区域仍残留一定数量的杂波； Eigen滤波器假设杂波沿声束轴向是平稳的，当过多的多普勒矢量用于构建自相关矩阵且滤波器阶数选择不当时，杂波抑制性能严重恶化，大量残留的杂波导致血流估计速度明显低于其它3种方法；Hankel-SVD根据慢时方向单个多普勒回波矢量的杂波成份所反映的频率特性，动态调整Hankel杂波主成份阶数，因而能有效抑制由于组织空间高强度时变运动导致的非平稳杂波，提高血流与组织区分能力。

图3为滤波后功率多普勒成像比较，显示为归一化多普勒功率，滤波器参数设置同图2。从图3中可以看出，经4种滤波器进行滤波处理后，杂波成份均得到了不同程度的抑制，远离血管区域的杂波抑制最为充分，而在血管壁及其附近区域，血管的非匀速搏动使得在血管壁边缘仍残留一定部分的较强功率杂波。图3(d)为Hankel-SVD滤波处理结果，可以看出杂波功率抑制最为充分，因而血流区域平均功率值明显高于其它3种滤波算法。

图1 滤波器杂波抑制性能比较

图2 彩色多普勒血流成像

图3 功率多普勒成像

6 结论

针对超声彩色血流成像中的非平稳杂波，本文提出了一种改进的基于奇异值分解的杂波抑制算法。该方法逐次采用慢时方向单个回波多普勒矢量构建Hankel矩阵，利用Hankel-SVD分解获得的Hankel主成分作为回归滤波器的正交基函数重构血流多普勒信号，可使杂波拟合最小均方误差最小，因而杂波抑制性能优于传统采用多项式函数作为正交基函数的回归型杂波抑制方法。根据Hankel主成份包含的频率内容区分杂波与血流信号，当杂波成份为主导时，通过估计杂波瞬时带宽与一阶Hankel主成分的平均多普勒频率确定杂波多普勒频率范围，并与高阶Hankel主成分平均多普勒频率进行比较，动态选择杂波Hankel主成份阶数，实现非平稳杂波的有效抑制。基于多普勒回波模型合成数据的仿真实验与颈动脉血流与功率成像结果表明该算法在有效抑制由血管加速运动导致的非平稳、强杂波的同时充分保留血流成份，改善低速血流参数的估计精度。由于仅采用单个慢时方向的多普勒回波矢量构建Hankel矩阵，本文算法具有空间自适应性，与Eigen杂波滤波器相比，无需人为设定阈值参数以估计杂波空间维数，因而可有效提升组织空间高强度时变运动时的血流成像质量。

[1] Oeltze S, Lehmann D J, and Kuhn A. Blood flow clustering and applications in virtual stinting of intracranial aneurysms[J]. IEEE Transactions on Visualization and Computer Graphics, 2014, 20(5): 686-701.

[2] Park G Y, Yeo S M, and Yoon C H. New adaptive clutter rejection based on spectral decomposition and tissue acceleration for ultrasound color Doppler imaging[C]. IEEE International Ultrasonics Symposium, Prague, 2013: 1484-1487.

[3] Masuda K, Bossard A, and Sugano Y. Reconstruction and error detection of blood vessel network from ultrasound volume data[C]. IEEE 26th International Symposium on Computer-Based Medical Systems, Porto, 2013: 497-501.

[4] Yap C H, Thiele K, and Wei Q F. Novel method of measuring vascular regurgitation using three-dimensional nonlinear curve fitting of Doppler signals within the flow convergence zone[J]. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2013, 60(7): 1295-1311.

[5] Kadi A and Loupas T. On the performance of regression and step-initialized iir clutter filters for color Doppler systems in diagnostic medical ultrasound[J]. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 1995, 42(3): 927-937.

[6] Bjaerum S, Torp H, and Kristoffersen K. Clutter filters design for ultrasound color flow imaging[J]. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2002, 49(2): 693-704.

[7] Yoo Y and Kim Y. New adaptive clutter rejection for ultrasound color Doppler imaging: in vivo study[J]. Ultrasound in Medicine & Biology, 2010, 36(3): 480-487.

[8] Kruse D E and Ferrara K W. A new high resolution color flow system using an eigen decomposition-based adaptive filter for clutter rejection[J]. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2002, 49(12): 1384-1399.

[9] Yu A C H. Eigen-based signal processing methods for ultrasound color flow imaging[D]. [Ph.D. dissertation], University of Toronto, 2007.

[10] Song F, Zhang D, and Gong X. Performance evaluation of eigendecomposition-based adaptive clutter filter for color flow imaging[J]. Ultrasonics, 2006, 44(12): 67-71.

[11] Yu A C H and Lovstaken L. Eigen-based clutter filter design for ultrasound color flow imaging: a review[J]. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2010, 57(5): 1096-1111.

[12] 刘松, 张水莲, 李世凯, 等. 一种低信噪比损失SVD 滤波部分频带干扰抑制算法[J]. 电路与系统学报, 2013, 18(1): 70-75.

Liu Song, Zhang Shui-lian, Li Shi-kai, et al.. A low distortion algorithm based singular value decomposition for partialband interference rejection[J]. Journal of Circuits and Systems, 2013, 18(1): 70-75.

[13] 徐彦凯, 双凯. 自适应奇异值分解瞬变信号检测研究[J]. 电子与信息学报, 2014, 36(3): 583-588.

Xu Yan-kai and Shuang Kai. Detection of transient signal based on adaptive singular value decomposition[J]. Journal of Electronics & Information Technology, 2014, 36(3): 583-588.

[14] Yu A C H and Lovstakken L. Eigen-based clutter filters for clolor flow imaging: single-ensemble vs. multi-ensemble approaches[C]. IEEE International Ultrasonics Symposium, New York, 2007: 1101-1104.

[15] Park G Y, Yeo S M, and Yoon C H. New adaptive clutter rejection based on spectral analysis for ultrasound color Doppler imaging: phantom and in vivo abdominal study[J]. IEEE Transactions on Biomedical Engineering, 2014, 61(1): 55-63.

[16] Yu A C H, Johnston K W, and Cobbold R S C. Frequency based signal processing in ultrasound color flow imaging[J]. Canadian Acoustics, 2007, 35(2): 11-23,

[17] Yu A C H and Cobbold R S C. Single ensemble based eigenprocessing methods for color flow imaging part 1 the Hankel-SVD filter[J]. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2008, 55(3): 559-573.

王录涛： 男，1979年生，博士后，研究方向为阵列信号处理与医学成像处理.

王 微： 男，1982年生，博士生，研究方向为雷达信号处理.

金 钢： 男，1958年生，研究员，博士生导师，从事雷达信号处

理、图像处理与流体力学方向研究.

Non-stationary Clutter Rejection Based on Hankel-SVD for Ultrasound Color Flow Imaging

Wang Lu-tao①Wang Wei②Jin Gang②

①(School of Automation Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China)

②(China Aerodynamics Research & Development Center, Mianyang 621000, China)

Effective rejection of the time-varying clutter originating from slowly moving vessels and surrounding tissues is very important for depicting hemodynamics in ultrasound color Doppler imaging. In this paper, a new adaptive clutter rejection method based on Hankel Singular Value Decomposition (Hankel-SVD) is presented for suppressing non-stationary clutter. In the proposed method, a Hankel data matrix is created for each slow-time ensemble. Then the orthogonal principle Hankel components can be obtained through the SVD of the Hankel data matrix. It achieves non-stationary clutter suppression by reconstructing the flow signal with only the high order principle Hankel components, which are estimated from the frequency content carried by the principle Hankel components. To assess its efficiency, the proposed Hankel-SVD based method is applied to synthetic slow-time data obtained from a Doppler flow model and carotid arterial complex baseband data acquired by a commercial ultrasound system (Sonix RP). The resulting flow and power images show that the proposed method outperforms the traditional IIR and polynomials regression filter in attenuation of high intense non-stationary clutter signal. It is also adaptive to highly spatially-varying tissue motion and can automatically select the order of the filter, which leads to improved distinguishing between blood and tissue regions compared to other eigen-based filters.

Color flow imaging; Singular Value Decomposition (SVD); Adaptive clutter rejection; Non-stationary clutter

TN911.7; R445.1

： A

：1009-5896(2015)04-0830-06

10.11999/JEIT140893

2014-07-07收到，2014-11-24改回

*通信作者：王录涛 wltuestc@163.com