基于先验估计和相关相位的实时超声弹性成像

崔少国 彭彩碧 胡 明 刘向阳

1(湖北医药学院计算机中心，十堰 442300)

2(湖北医药学院附属太和医院，十堰 442300)

引言

超声弹性成像是生物组织成像的新模式，是当前超声医学工程学领域的前沿和热点研究。该成像技术可以探测出人体组织与病理变化相关的硬度改变，对肿瘤(癌症)的早期检测和诊断提供重要依据。弹性成像研究自1991年被Ophir等提出以来的20年中取得了快速发展，在临床中已逐渐显示出重要诊断价值，如在乳房癌、前列腺癌、甲状腺结节、肝纤维化、动脉粥样硬化等的检测中取得了重要应用［1－4］。

基于时间域的互相关成像技术主要分为两步:位移估计和应变计算。位移估计主要是通过互相关技术从组织压缩后的回波信号里找到与压缩前匹配的信号，匹配窗之间的时延即代表着组织的位移;应变计算主要是通过位移的梯度操作或差分滤波器来实现。这两步中位移估计是关键［5］。为了达到亚采样点级精度位移估计，Céspedes等提出了抛物线插值、余弦插值和重构插值等方法［6］。这些方法在一定程度上提高了位移估计精度，但插值操作计算代价比较大，降低了弹性成像的实时性。另外，位移估计的前提是从压缩后的回波信号里去搜索压缩前信号窗运动到的新位置，即运动追踪。使用穷举式搜索或启发式穷举搜索进行运动追踪搜索范围都比较大，是较耗时的追踪策略［7］。针对基于互相关成像算法中存在的这些问题，本研究提出了使用先验估计在较小的范围内进行运动追踪，使用复互相关的相位进行亚采样点级位移精度估计的新算法。

1 算法

提出的新算法主要分为两部分:基于先验估计的运动追踪和基于相关相位的亚采样点级形变估计。基于先验估计运动追踪的基本思想是根据连续性原则，物理位置上相邻点其运动的位移也是近似的，使用相邻位置上已估计出来的位移做为它的初始种子位移，然后在其附近较小的范围内追踪从而找到其真正的位移。基于相关相位的亚采样点级精度位移估计是先使用互相关技术在基带信号上进行位移初步估计(只精确到整数个采样点)，然后在该点使用复互相关函数的相位角进行亚采样点级精度位移估计。

1.1 运动追踪

1.1.1 种子位移选择

先验估计引导追踪可以使用轴向上相邻的上一已估计点的位移做为该点的种子位移(轴向引导)，也可以使用侧向上左边或右边相邻已估计点的位移作为该点的种子位移(侧向引导)。轴向引导算法示意图如图1(a)。用seed(i，j)表示第(i，j)个估计位置的初始种子位移，u(i，j)表示第(i，j)个估计位置的位移，则种子选择方法为:

位移提炼则使用:

式中，2≤i≤M，1≤j≤N，M 和 N 是位移估计图像的行和列维数。－S≤r≤S，S是相邻两个估计位置可能的最大位移差，［－S，S］是局部搜索范围，搜索步长1个采样点。Coefficient＜＞是计算互相关函数的系数，argmax()是求当互相关系数值最大时局部位移值r，round()是四舍五入函数。

侧向引导算法如图1(b)所示，假设最中间一列位移已被估计，则种子选择可按如下公式进行

位移提炼同样使用公式(2)进行。

1.1.2 局部搜索范围确定

种子位移确定后，需要确定局部搜索范围。局部搜索范围可以依据下面公式进行:

图1 轴向引导和侧向引导运动追踪示意(～u表示种子位移，û代表以～u为种子估计出来的位移，空心箭头代表种子选择来源，实心箭头代表种子传递方向)。(a)轴向引导;(b)侧向引导。Fig.1 Illustrations for various guidance tracking strategies，where～u denotes the seed displacement and û denotes the estimated displacement using the～u as seed.(a)Axial Guidance;(b)Lateral Guidance.

式中，Rup、Rdown、Rleft、Rright分别是局部搜索范围的向上、向下、向左、向右的维数。εmax是系统检测到的最大可能应变(属于系统设定参数，一般为6%)，Axial_shift和Lateral_shift分别是估计窗轴向和侧向移动的维数。一维轴向运动追踪时只采用式(4)，二维运动追踪时采用式(4)和式(5)。

1.1.3 局部启发式搜索

如果组织的全局运动方向已知，局部搜索范围可进一步缩小。全局运动方向的获取可以先随机地选取几个估计点，这些点的位移可通过启发式穷举搜索获得。然后根据这些位移的正负号确定组织全局运动方向。如果向上(压缩)，则局部搜索范围

如果向下(拉伸)，则

向左和向右搜索维数Rleft、Rright不变，使用式(5)。

1.2 位移估计

设压缩前的回波信号为x(t)，压缩后的回波信号为y(t)。则x(t)和y(t)的复互相关函数为

2-D离散化的归一化复互相关函数为

根据式(9)，当复互相关函数 f(i，j，d)的值最大时所求得的d即为整数采样点级位移估计。那么亚采样点级精度的位移估计u为

式中，angle()是求复互相关函数 f(i，j，d)的相位角，ω0是中心角频率。

1.3 应变计算

应变计算使用罗建文等提出的数字低通差分滤波器SG-I DD［8］，轴向应变计算方法为

式中，s表示应变，u表示位移，M是滤波器参数(滤波长度为2M+1)，Δk=u(n+k)－u(n－k)。

2 实验

在MatLab 7.0实现上述提出的成像算法，并在仿真数据集和弹性体模上获取的真实回波数据集上运行，比较提出的算法和经典的互相关插值法的成像速度和成像质量。插值法选用抛物线插值，比较的具体项目包括成像一帧应变图像所需要的时间、成像效果、应变图像的性噪比SNRe。成像所使用的计算机是一台 Pentium(R)4，CPU 2.4 GHz，RAM 1 GB，操作系统是 Microsoft Windows XP。

2.1 仿真

本实验仿真中使用的中心频率是5 MHz，采样频率是20 MHz，声速设定为常数1 540 m/s，两种模型被仿真:均匀弹性组织模型和含有两个硬包容物的组织模型。组织大小均为40 mm×40 mm，包容物的直径是10 mm。散射子浓度设为45(个散射子/脉冲宽度)，硬包容物的弹性模量是背景的10倍。1D应变模型被采用(设 dy=0)。仿真的每帧数据是1039(采样点)×128(线)，然后加高斯白噪声使其性噪比为50 dB。对每种模型产生不同整体应变的数据，应变分别是:0.01%，0.1%，0.5%，1%，2%，3%，4%，5%。

2.2 体模实验

为了验证提出的算法的真实成像效果，使用弹性体模上采集的超声回波数据进行实验。体模采用CIRS Inc.研发的专用于弹性成像研究的Model 049弹性体模。体模里在背景介质中放有4个10 mm和4个20 mm直径的球形包容物，球体和背景是等回声的，在B模式图像中不能区分。背景的弹性模量是30 kPa，，每种直径的4个包容物各有不同的弹性模量，它们分别是 8 kPa(TypeⅠ，软)，18 kPa(TypeⅡ，软)，44 kPa(TypeⅢ，硬)和 63 kPa(TypeⅣ，硬)。实际体模中不同直径的包容物在轴向(上下方向，即体模压缩方向)上是错开的。本实验中成像对象选择直径是10 mm、弹性模量是63 kPa的较硬包容物。

实验采用的探头是型号为SA5L38B的128阵元的线阵探头，中心频率是5 MHz，75%分数阶带宽。实验使用的超声系统是声泰特(成都)公司研发的iMago C21超声机，系统采样频率设为 40 MHz。使用手持探头均匀轴向压缩体模，在压缩过程中共采集序列回波300帧。RF信号使用积分解调器解调到基带I/Q信号，然后下采样到10 MHz。帧间应变大约0.7%，帧大小是 512(采样点)×188(线)。

3 结果及分析

3.1 时间比较

新算法与经典的互相关插值算法的时间性能比较分别如表1和表2。从表中可看出，不论是仿真还是体模实验，新算法比经典算法成像速度都大大提高，这主要是因为使用了先验位移进行位移估计，使运动追踪的范围大大减小。经典互相关插值算法实时性较差，每秒只能达到约2～3帧，而新算法实时性较好，每秒可以达30～50帧，可满足临床超声成像的实时性要求。

3.2 成像效果

仿真成像效果如图2所示。使用含两个硬包容物仿真模型，成像中估计窗大小为40个采样点，窗轴向重叠率0.75。从图2可以看出，新算法所得到的应变图像(图2(b))比较细腻平滑，噪声显著减少。这主要是因为采样相关相位进行亚采样点级精度位移估计比抛物线插值法更准确。

表1 仿真实验中经典算法和新算法的时间性能比较Tab.1 Comparison of computation time for conventional and presented algorithms in simulations

表2 体模实验中经典算法和新算法的时间性能比较Tab.2 Comparison of computation time for conventional and presented algorithms in phantom

图2 仿真弹性成像。(a)抛物线插值法;(b)复相关相位法Fig.2 Simulation strain images. (a)Parabola interpolation;(b)Cross－correlation phase

体模成像效果如图3所示，成像中估计窗大小20个采样点，窗轴向重叠率0.75。从图3可以看出，体模实验结果与仿真实验结果一致，新算法得到的应变图像(图3(b))有更好的性能，弹性图像中的噪声明显减少。

图3 体模弹性成像。(a)抛物线插值法;(b)复相关相位法Fig.3 Phantom strain images. (a)Parabola interpolation;(b)Cross－correlation phase

3.3 应变图像的性噪比

图4显示了在仿真的均匀弹性组织模型数据集上应用提出的新算法所得应变图像的性噪比与抛物线插值法的比较。性噪比计算公式为

式中μs和 σs分别是均匀弹性组织应变图像的应变平均值和应变标准方差。图中曲线上每一个点是5次独立仿真实验数据的平均值。图中最上面的一条曲线是使用相关相位(新算法)法在各种应变数据上成像的性噪比，中间一条是抛物线插值法对应的曲线，最下面一条是没有进行亚采样点级位移估计(估计精度是整数个采样点)的产生的曲线。从图中可以看出当应变s＜1%时，性噪比随应变减小而减小，这主要是因为应变较小时，回波信号的噪声是应变估计的主要噪声源，回波噪声使较小应变的数据窗匹配产生较大误差;当应变 s＞3%时，性噪比随应变增大而下降，这主要是因为应变较大时，解相关是应变估计的主要噪声源，高的解相关使位移(时延)估计产生较大的方差。另外，从图中可以明显看出，在任何应变下，新算法的性噪比SNRe最高，抛物线插值法次之，整数个采样点级算法所得的性噪比最低。尤其在1%和2%的应变时，新算法的性噪比较抛物插值法显著提高。

图4 抛物线插值法、相关相位法在仿真的各应变数据上应变成像的性噪比Fig.4 SNRe of strain images from the parabola interpolation and correlation phase algorithms in simulations

4 结论

本研究提出了一种新的实时高精度弹性成像算法，该算法根据连续性原则使用先验估计进行运动追踪，大大缩小运动追踪范围，提高成像速度;同时使用复相关函数相位进行亚采样点级精度位移估计，使应变估计精度大大提高，图像品质明显改善。仿真实验和体模实验显示，新算法与传统的基于互相关的抛物线插值法相比在成像速度、应变图像的性噪比等方面均有较大提高。实验证明了提出算法的正确性和高效性。该算法的提出将对组织弹性实时成像系统研发和临床推广应用具有重要意义。

［1］Ophir J，Céspedes I，Ponnekanti H，et al.Elastography:a quantitative method for imaging the elasticity of biological tissue［J］.Ultrasonic Imaging，1991，13(3):111 －134.

［2］Céspedes I，Ophir J，Ponnekanti H，et al.Elasticity imaging using ultrasound with application to muscle and breast in vivo［J］.Ultrasonic Imaging，1993，15(2):73－88.

［3］RighettiR.Elastographic characterization ofHIFU-induced lesions in canine livers［J］. Ultrasound in Medicine and Biology，1999，25(7):1099－1113.

［4］Konofagou E，Ophir J.Myocardial elastography—A feasibility study［J］.Ultrasound in Medicine and Biology，2002，28(4):1113－1124.

［5］Lindop JE，Treece GE，Gee A，et al.Phase-based ultrasonic deformation estimation［J］.IEEE Transactions on Ultrasonics，Ferroelectrics，and Frequency Control，2008，55(1):94 －111.

［6］Céspedes I，Huang Y， Ophir J.Methods for estimation of subsample time delays of digitized echo signals［J］.Ultrasonic Imaging，1995，17:142 －171.

［7］Chen Lujie，Treece GM， Lindop JE.A quality-guided displacement tracking algorithm for ultrasonic elasticity imaging［J］.Medical Image Analysis，2009，13:286 －296.

［8］Luo Jianwen，Bai Jing，He Ping，et al.Axial Strain Calculation Using a Low-Pass Digital Differentiator in Ultrasound Elastography ［J］. IEEE Transactions on Ultrasonics，Ferroelectrics，and Frequency Control，2004，51(9):1119 －1127.