CN102764139B

CN102764139B - 基于特征空间分析和区域判别的医学超声波束形成方法

Info

Publication number: CN102764139B
Application number: CN201210240918.4A
Authority: CN
Inventors: 汪源源; 陈呈; 余锦华; 曾星
Original assignee: Fudan University
Current assignee: Fudan University
Priority date: 2012-07-12
Filing date: 2012-07-12
Publication date: 2014-10-29
Anticipated expiration: 2032-07-12
Also published as: CN102764139A

Abstract

本发明属于医学超声工程技术领域，具体为基于特征空间分析与区域判别的医学超声波束形成。本发明一方面根据超声回波数据特征矢量与导向矢量的投影关系，预估计期望信号子空间以及期望信号强度；另一方面利用延时叠加法获得一幅预扫描超声图像，根据其峰值和有效动态范围，计算出一个信号强度阈值；最后将期望信号强度估计值与该阈值比较，对超声低回声区与一般散射区进行判别，低回声区采用主分量最小方差算法，一般散射区采用传统最小方差算法。本发明适用于高帧率的平面波和合成孔径超声成像方式，能够显著改善成像对比度，同时较好地保持图像的散斑特性，降低人工干扰和失真。

Description

基于特征空间分析和区域判别的医学超声波束形成方法

技术领域

本发明属于医学超声工程技术领域，具体为基于特征空间分析与区域判别的医学超声波束形成方法。

背景技术

超声成像具有无创、实时、可重复性强、成本低廉等优点，在医学临床诊断中得到了广泛的应用。波束形成器是整个系统中影响成像质量的关键部分，传统的延时叠加(delay-and-sum, DAS)法[1]虽然快速简便，但获得的图像质量十分不理想，主要体现在横向分辨率和对比度差，这也是目前医学超声诊断面临的主要难题。

为改善图像质量，近来基于最小方差(minimum variance, MV)[2]的自适应波束形成方法被引入超声成像领域。其中空间平滑最小方差[3]、前后向平均最小方差[4]虽然显著改善了分辨率，但是对比度性能的提升并不明显；主分量最小方差[5]和相干系数处理[6]虽然能够有效提高超声低回声区与散斑背景的对比度，但是高强度目标周围易产生人工干扰和失真，原因来自于方法内在的缺陷。

针对上述问题，本发明通过进一步的特征空间分析，对成像目标的低回声区和散斑背景以及高回声区加以判别，并对不同区域采用不同的处理方法，以结合传统最小方差算法和主分量最小方差算法的优势。本发明为医学超声波束形成提供了一种的新的处理框架，能够有效地实现高分辨率和高对比度，同时保持良好的散斑特性，降低人工干扰和失真。

发明内容

本发明的目的是提出一种可以有效改善超声成像对比度并保持散斑特性和降低失真的医学超声波束形成方法。

本发明提出的医学超声波束形成方法，其具体步骤为：

1、对超声回波数据协方差矩阵进行特征分解，根据特征矢量与导向矢量的投影关系，预估计期望信号子空间以及期望信号强度；

2、利用延时叠加法获得的预扫描超声图像，根据图像峰值和有效动态范围计算出一个信号强度阈值；

3、将期望信号强度估计值与该阈值比较，对超声低回声区与一般散射区进行判别，两类区域采用不同的波束形成方法：低回声区采用主分量最小方差算法，一般散射区采用传统最小方差算法。

医学超声波束形成器的工作原理如图1所示，首先对探头的不同阵元施加适当的延时，实现对成像点的聚焦，再对聚焦后的阵列信号以适当的权重进行加权平均，输出信号即聚焦点散射信号的估计，取其幅值作为该点的成像值，反映的是该点的散射强度。波束形成器的输出可表示为：

(1)

其中o代表成像点，x _d(o)表示经过聚焦延时的输入信号，w为权矢量，H为共轭转置。(1)式表示权矢量可以随着成像点变化，即自适应的。对于平面波(plane wave, PW)和合成孔径(synthetic aperture, SA)成像，采用的是逐点合成聚焦，对所有成像点重复(1)式后得到最终图像，因此不失一般性，后文将省略o。

为对超声信号去相关，采用空间平滑和对角加载技术估计得到的数据协方差矩阵为[3]：

(2)

其中M为探头阵元数，L为空间平滑子阵长度，x _d ^l为子阵l的数据矢量，I为单位阵，δ为对角加载因子。

本发明方法的流程可总结如图2。下面就本发明方法各个步骤涉及的相关技术细节作进一步的描述。

1、数据协方差矩阵R可按下式特征分解：

(3)

其中U和Λ分别是L×L维的复数矩阵，λ _i为矩阵R的降序排列的第i个特征值，u _i为其对应的L×1维的复数特征矢量。

考虑u _i与L×1维的复数导向矢量a的归一化投影值：

(4)

在医学超声成像中，经过延时聚焦后，a的元素将全为1。

选取最大的若干p _j，使得：

(5)

其中T_project为阈值，一般可以取0.8，后文也讨论了该参数对方法性能的影响。(5)式表示导向矢量与取到的特征矢量方向性足够接近，这些u _j构成的子空间可近似为期望信号子空间，对应的特征值均值作为期望信号强度的估计值：

(6)

2、现有的主分量最小方差算法虽然能够较好地抑制低回声区的干扰信号，但是在散斑背景中，尤其是高强度目标周围，容易把相对较弱的有用信号一并滤除，这会造成严重的失真，原因在于其统一的处理方式不能兼顾这两类区域的特性。针对该问题，本发明试图从区域判别的思路入手，对两类区域采用不同的处理方法，以实现更高的综合性能。

超声低回声区和一般散射区的重要差异在于期望信号的绝对强度，在低回声区，期望信号强度明显低于一般散射区。这里利用事先获得的DAS预扫描超声图像确定一个阈值，当期望信号强度低于该阈值时，认为成像点位于低回声区；反之则认为成像点位于一般散射区。

首先根据预扫描超声图像峰值推算整个成像区域的全局最大特征值(最大散射强度)。波束形成器的输出能量可表示为：

(7)

其中w是L×1维的复数波束形成加权矢量，cos<·>表示两矢量夹角的余弦值。

对于传统延时叠加法，w=(1/L)[1 1 … 1]^T。图像峰值处是整个成像区域中散射强度最大的位置，全局最大特征值也应当在该位置。对该点成像时，可近似认为最大特征值对应的特征矢量u ^max与导向矢量方向是一致的，即也与w一致，从而cos<u ^max, w>≈1，且其他特征矢量满足cos<u _i, w> ≈ 0，因此有：

(8)

其中λ ^max即得到的全局最大特征值，y _max为预扫描超声图像峰值。

接着估计预扫描图像相对于峰值的有效动态范围。利用图像直方图，寻找到一个灰度值，使得像素值在此之上的像素总能量与图像总能量的比值足够大，设大于T_intensity (0< T_intensity<1)。该灰度值对应了一个相对于峰值的动态范围，该范围内包含了图像的大部分能量。一般可采用T_intensity = 0.8，下文也讨论了该参数对性能的影响。

最后将全局最大特征值与有效动态范围相除，即得到所需的强度阈值，在此阈值之下的期望信号被认为来自超声低回声区，反之为一般散射区。

3、将步骤1得到的期望信号强度估计值与步骤2得到的强度阈值比较：若估计值小于阈值，则认为当前成像点为超声低回声区；若估计值大于阈值，则认为当前成像点位于一般散射区。

对于一般散射区，采用传统最小方差算法，能够很好地保持散斑特性和降低失真。求解下面的最优化问题：

(9)

通过拉格朗日乘子法解得MV权矢量为[2]：

(10)

对于超声低回声区，采用主分量最小方差算法，以达到更好的干扰抑制效果，提高对比度。具体为将R的特征值降序排列，设最大特征值为λ ₀，取大于λ ₀/10的特征值对应的特征矢量构成子空间U _s，将传统MV权矢量向U _s投影即得到主分量MV权矢量[5]：

(11)

得到波束形成权矢量后，该成像点的最终输出为：

(12)。

附图说明

图1 医学超声波束形成器示意图。

图2 本发明方法流程图。

图3 散斑背景中点散射子的仿真结果：(a)DAS矩形窗，(b)DAS海明窗，(c)MV，(d)前后向平均MV，(e)主分量MV，(f)本发明方法。

图4 散斑背景中的点散射子仿真：不同方法轴距37 mm处的横向幅值比较。

图5 散斑背景中吸声囊肿的仿真结果：(a)DAS矩形窗，(b)DAS海明窗，(c)MV，(d)前后向平均MV，(e)主分量MV，(f)本发明方法。

图6 散斑背景中的吸声囊肿仿真：不同方法对比度指标(CR、CNR)比较。

图7 散斑背景中的吸声囊肿仿真：本发明方法随参数变化的性能比较，(a)-(b) CR、CNR随T_intensity的变化，(c)-(d) CR、CNR随T_project的变化。

图8 一种更复杂情况的仿真结果：(a)DAS矩形窗，(b)DAS海明窗，(c)MV，(d)前后向平均MV，(e)主分量MV，(f)本发明方法。

图9 铜丝实验结果：(a) DAS矩形窗，(b) DAS海明窗，(c)MV，(d)前后向平均MV，(e)主分量MV，(f)本发明方法。

具体实施方式

对本发明提出的波束形成方法进行Field II [7]仿真和实际实验测试，将本发明方法与传统DAS(矩形窗、海明窗)、MV、前后向平均MV以及主分量MV等方法进行了性能比较。未特别指出时，本发明方法中参数T_project和T_intensity均取0.8。

仿真采用96阵元、0.308 mm间距的线阵和5 MHz、60%带宽的超声脉冲，平面波成像。

图3为不同方法对散斑背景中点散射子的仿真结果，图4为轴距37 mm处的横向幅值比较。主分量MV产生了散斑背景的失真，而本发明方法保持了散斑特性和较窄的散射子主瓣，体现了比DAS更好的分辨率。在该仿真中，本发明方法与传统MV以及前后向平滑MV具有相似的性能。

图5为不同方法对散斑背景中吸声囊肿的仿真结果，图6为这些结果的对比度指标比较。这里采用了对比度比值(contrast ratio, CR)[8]和对比度噪声比值(contrast-to-noise ratio, CNR)[9]两个指标来衡量图像的对比度性能。其中CR定义为吸声囊肿内外平均幅值之差：

(13)

CNR定义为CR与散斑背景幅值标准差之比：

(14)

DAS方法得到的图像对比度和分辨率最差；传统MV和前后向平均MV的分辨率有明显提升，对比度也有所增强；主分量MV的分辨率并没有进一步的改善，而囊肿的对比度有显著提高，但是在强散射子周围却产生了严重失真，破坏了原有的散斑特性；本发明方法消除了这种失真，而且基本维持了主分量MV的高对比度，CNR更是所有方法中最高的。

图7显示了本发明方法中参数T_project和T_intensity对性能的影响，考察了CR和CNR随两个参数从0.3增加到1.0的变化情况。首先可以看到，两个参数在0.8以下时对CR的影响都不大，其中T_project对CNR也几乎没有影响，可见在取值并不极端的情况下本发明方法对参数T_project的敏感性较低，取值不高于0.8即可。而参数T_intensity与图像动态范围的计算有关，取得太小不能包含图像的主要能量，取得太大则会包含过多的噪声成分，对最终结果的CNR有一定影响，从图7中可以看到0.7 ~ 0.8的取值是比较合适的。

图8为不同方法对一种更复杂情况的仿真结果，仿真对象包含了超声低回声区、高回声区和离散强散射子。结果进一步验证了本发明方法对复杂情况的鲁棒性，能够保持良好的散斑特性、降低失真，并维持较高的对比度。

此外，还将本发明方法应用于仪器实验采集的数据。实验中使用128阵元、0.308 mm间距的线阵，5 MHz超声脉冲，采样率为40 MHz，平面波成像。目标为均匀介质中的若干根细铜丝。

图9为不同方法得到的超声图像。可以看到本发明方法能够较好地消除主分量MV产生的失真，保护了原有的散斑特性。此外还发现，对于实际数据，本发明方法与MV以及前后向平均MV一样，点目标的分辨力相比于DAS方法，并没有仿真中那样明显的改善，这是可能是由实际应用中探头的误差、噪声以及超声传播的衰减所致。

综上所述，本发明作为一种新的超声波束形成方法，能够在维持较高分辨率、显著增强对比度的同时，保持超声图像原有的散斑特性，降低甚至消除人工干扰和失真。

参考文献

[1] Thomenius K E. Evolution of ultrasound beamformers [C]. Proceedings of the 1996 IEEE Ultrasonics Symposium, 1996: 1615-1622.

[2] Capon J. High resolution frequency-wavenumber spectrum analysis [C]. Proceedings of IEEE, 1969: 1408-1418.

[3] Synnevag J F, Austeng A, Holm S. Adaptive beamforming applied to medical ultrasound imaging [J]. IEEE Transaction on Ultrasonics, Ferroelectrics, and Frequency Control, 2007, (8): 1606-1613.

[4] Asl B M, Mahloojifar A. Contrast enhancement and robustness improvement of adaptive ultrasound imaging using forward-backward minimum variance beamforming [J]. IEEE Transaction on Ultrasonics, Ferroelectrics, and Frequency Control, 2011, (4): 858-867.

[5] Asl B M, Mahloojifar A. Eigenspace-based Minimum Variance beamforming applied to medical ultrasound imaging [J]. IEEE Transaction on Ultrasonics, Ferroelectrics, and Frequency Control, 2010, (11): 2381-2390.

[6] Nilsen C C, Holm S. Wiener beamforming and the coherence factor in ultrasound imaging [J]. IEEE Transaction on Ultrasonics, Ferroelectrics, and Frequency Control, 2010, (6): 1329-1346.

[7] Jensen J A. Field: a program for simulating ultrasound systems [J]. Medical and Biological Engineering and Computing, 1996, (1): 351-353.

[8] O’Donnell M, Flax S W. Phase-aberration correction using signals from point reflectors and diffuse scatterers: measurements [J]. IEEE Transaction on Ultrasonics, Ferroelectrics, and Frequency Control, 1988, (6): 768-774.

[9] Krishnan S, Rigby K W, O’Donnell M. Improved estimation of phase aberration profiles [J]. IEEE Transaction on Ultrasonics, Ferroelectrics, and Frequency Control, 1997, (3): 701-713.。

Claims

1.基于特征空间分析和区域判别的医学超声波束形成方法，其特征在于具体步骤为：

(1) 对超声回波数据协方差矩阵进行特征分解，根据特征矢量与导向矢量的投影关系，预估计期望信号子空间以及期望信号强度；

(2) 利用延时叠加法获得的预扫描超声图像，根据图像峰值和有效动态范围计算出一个信号强度阈值；

(3) 将期望信号强度估计值与该阈值比较，对超声低回声区与一般散射区进行判别，两类区域采用不同的波束形成方法：低回声区采用主分量最小方差算法，一般散射区采用传统最小方差算法；

在步骤(1)中，所述对超声回波数据协方差矩阵进行特征分解，根据特征矢量与导向矢量的投影关系，选取投影值较大的若干特征矢量构成期望信号子空间，对应的特征值均值作为期望信号强度估计值，具体步骤为：

数据协方差矩阵R按下式特征分解：

(1)

其中U和Λ分别是L×L维的复数矩阵，λ _i为矩阵R的降序排列的第i个特征值，u _i为其对应的L×1维的复数特征矢量，H为共轭转置；

考虑u _i与L×1维的复数导向矢量a的归一化投影值p _i：

(2)

在医学超声成像中，经过延时聚焦后，a的元素将全为1；

选取最大的若干p _j，使得：

(3)

其中T_project为阈值， (3)式表示导向矢量与取到的特征矢量方向性足够接近，这些u _j构成的子空间近似为期望信号子空间，对应的特征值均值作为期望信号强度的估计值：

(4)；

在步骤(2)中，所述利用延时叠加法获得的预扫描超声图像，根据图像峰值和有效动态范围计算出一个信号强度阈值，具体步骤为：

首先，根据预扫描超声图像峰值推算整个成像区域的全局最大特征值，波束形成器的输出能量表示为：

(5)

其中w是L×1维的复数波束形成加权矢量，cos<·>表示两矢量夹角的余弦值；

对该点成像时，近似认为最大特征值对应的特征矢量u ^max与导向矢量方向一致，即也与w一致，从而cos<u ^max, w> ≈ 1，且其他特征矢量满足cos<u _i, w> ≈ 0，因此有：

(6)

其中λ ^max即得到的全局最大特征值，y _max为预扫描超声图像峰值；

接着，估计预扫描超声图像相对于峰值的有效动态范围，利用图像直方图，寻找到一个灰度值，使得像素值在此之上的像素总能量与图像总能量的比值足够大，设该比值大于阈值T_intensity ，0 < T_intensity< 1；寻找到的灰度值对应一个相对于峰值的动态范围，该范围内包含了图像的大部分能量；

最后，将全局最大特征值与有效动态范围相除，即得到所需的强度阈值，在此阈值之下的期望信号被认为来自超声低回声区，反之为一般散射区。

2.根据权利要求1所述的基于特征空间分析和区域判别的医学超声波束形成方法，其特征在于步骤(3)中，所述将步骤(1)得到的期望信号强度估计值与步骤(2)得到的强度阈值比较：若估计值小于阈值，则认为当前成像点为超声低回声区；若估计值大于阈值，则认为当前成像点位于一般散射区；

对于一般散射区，采用传统最小方差算法，即求解下面的最优化问题：

(7)

其中w为波束形成加权矢量，R和a分别为数据协方差矩阵和阵列聚焦后的导向矢量；

最小方差权矢量可通过拉格朗日乘子法解得：

(8)

对于超声低回声区，采用主分量最小方差算法，具体是将R的特征值降序排列，最大特征值为λ ₀，取大于λ ₀/10的特征值对应的特征矢量构成子空间U _s，传统最小方差权矢量在U _s上的投影即主分量最小方差权矢量：

(9)。