基于空间位移复合的超声弹性成像抑噪

邵党国 钟 明 易三莉 相 艳 马 磊 贺建峰

（1.昆明理工大学信息工程与自动化学院，昆明，650500；2.昆明理工大学国际学院，昆明，650500）

引 言

触诊是一种主观地检测人体组织中硬块的方法，超声弹性成像通过计算和显示组织弹性应变，弥补了触诊的这一缺点。超声弹性成像是一种生物组织成像的新模式，该成像技术可以探测出人体组织与病理变化相关的硬度改变，对肿瘤（癌症）的早期检测和诊断具有重要意义。弹性成像的研究自1991年被Ophir等人［1］提出以来，20年中取得了快速发展，在临床中已逐渐显示出重要诊断价值，比如在乳腺癌、前列腺癌、甲状腺结节、肝纤维化、动脉粥样硬化等的检测中取得了重要的应用。其成像的基本原理是：对组织施加一个动态或者准静态的压力，组织将产生一个响应，在相同的应力下硬的组织形变较小，而软的组织形变较大；使用信号处理和图像处理相关技术计算出组织的应变大小，进而重构出组织的弹性模量；对组织的应变或弹性模量进行成像就是弹性成像［1］。

弹性成像的生成过程包括两个主要步骤：位移估计和应变估计。一幅弹性成像图像可以通过显示的位移场空间而衍生。一般来说，只需要考虑轴向应变信息，而忽略横向应变信息。但衍生过程趋于噪声放大，产生伪影噪音。弹性成像过程中产生的伪影从物理特性上不同于B模式图像的斑点［2］，弹性成像的伪影主要由两个原因造成：解相关误差和幅度调制误差。解相关误差主要是由于组织的压缩、形变、滑动使运动估计不能遵守时延估计模型；幅度调制误差主要是由于信号幅度的随机波动使估计的位移值不是估计窗口中心点的位移。一些研究小组已经提出了一些抑制这种伪影的方法，如：自适应拉伸［3］、小波去噪［4］、幅度调制修正［5］等。

在传统的B模式超声成像中，复合成像技术由于能有效抑制斑点噪声提高图像质量而被广泛应用［6-8］。该技术主要是利用不同子图像的斑点模式解相关性，将多个图像进行加权复合得到噪声抑制。该技术主要有两类：频率复合成像和空间复合成像。频率复合是使用不同频率的超声波对同一感兴趣区域（Region of intevest，ROI）进行成像，然后使用某一加权策略对个子图像进行复合；空间复合是从不同空间角度对同一ROI进行成像，然后使用某一加权策略对各子图像进行复合。以Varghese和Techavipoo为代表的许多研究小组将空间复合［9-13］技术应用到弹性成像中，其基本原理是将不同扫描角度产生的弹性图像转换后加权平均产生了一个复合后的弹性图像。实验结果表明，该方法能有效抑制弹性成像的伪影噪声。本文从伪影产生的原理出发，提出位移空间复合技术，对运动估计的位移进行复合校正，从而达到抑制伪影噪声的目的。本文研究的位移空间复合方法是通过线性探头交替斜角扫描帧组成的，这些帧数据从声泰特C21系列超声系统中获得，此方法使用商业弹性体模进行了测试，通过计算信噪比（SNRe）和对比噪声比（CNRe）的分析，验证了此方法可以抑制伪影噪声，从而改进超声弹性图像质量。

1 空间复合原理

超声图像的质量通常通过信噪比（SNRe）来评价［13］

式中：μ表示均匀组织区域的平均值；σ表示均匀组织区域的标准差。

先设有两个子图像P和Q，它们的大小均为m×n，它们之间的相关性系数ρ（0≤ρ≤1）表示为

式中：Pij和Qij分别为图像P和Q中位置（i，j）上的灰度值；和分别为图像P和Q的灰度平均值。

定义N个随机变量Xi（i＝1，2，…，N，N≥2）表示N个子图（来自于同一个感兴趣区域ROI）。将这N个子图使用平均加权策略复合成一个新的图像X，则图像X的方差为

式中Cov（Xi，Xj）表示的协方差。假设Var（Xi）近似等于σ2，ρ（Xi，Xj）近似等于ρ，则式（3）可简化为

如果把图像的平均值看作一个常数μ，则N个子图像平均复合后所得图像的数学期望也是这个值，即

由式（1，4，5）得复合后图像X的信噪比为

从式（6）可以得出如下结论：如果各子图像间完全解相关（即ρ＝0，随机变量Xi完全独立），则复合后图像的信噪比是原来的倍；如果各子图像间完全相关（即ρ＝1），则复合后图像的信噪比没有变化；如果各子图像间部分相关（即0＜ρ＜1），则复合后图像的信噪比是原来的倍。

通过以上的推导可以得出，只要各子图像之间具有解相关性，复合后的图像信噪比就能得到一定的提高，且解相关性越强，信噪比提高越显著。

2 算法描述

2.1 运动位移估计

在超声弹性成像中，位移估计是第一个关键的步骤。迭代零-相位算法［14］具有计算效率高和估计准确的特点。本文采用迭代的零-相位算法来求解位移。在沿着扫描线方向，假设前一时刻的扫描线信号为y1，后一个时刻同一位置的扫描线信号为y2。前后两个时刻的同一位置的扫描线信号仅存在刚性平移，则x2（t＋τ）＝x1（t）。两个信号的互相关函数定义为

根据互相关函数的性质，当延时时间为0时，也即两个信号完全重合时，互相关函数取得最大值。

要求的运动位移即为C×τ，这里C为超声波声速，是个常量。利用互相关函数的定义和性质，可以在对超声斑点跟踪的过程中，通过寻找同一超声斑点的相邻运动帧中的两个信号的相关函数的最大值，从而求得相对应的时间延迟τ。设两个解析信号的相位函数为

由自相关函数零延时的值是正实数的性质，可以知道此时它的相位为0，由此可以推知解析信号x1a和x2a的互相关函数的相位函数在t＝-τ时为0。

由解析信号互相关函数的性质可知，在其最大值的附近，互相关函数的相位函数近似于一个线性函数，其斜率为探头载波的中心角频率ω0，利用这一性质，用牛顿迭代法来求解互相关最大时的延时τ

式中：n表示迭代的次数；τ0一般取0。

2.2 位移空间复合

一幅图像的解相关子图像，可以通过不同频率扫描获得，也可以通过不同空间角度扫描获得，在算法的第二步中，本文采用后者。本文先后改变探头线性扫描序列的转向角，以3个角为例，按顺序发送角度A，B和C，连续序列为“…ABCABCABC…”。在位移估计中，应用相同的连续帧。这里使用轴向扫描转换1/（cos2θ-ν＊sin2θ），ν是泊松比，对于不可压缩的组织ν＝0.495［13］。扫描转换之前，运动估计代表着沿探头扫描线方向的运动位移；通过扫描转换，希望得到统一的轴向方向的运动位移。在统一方向的前提下，对不同空间角度的运动估计标量进行平均复合，达到校正位移从而抑制伪影噪声的目的。复合算法如下

2.3 弹性图像

弹性成像的最后一步，是从位移估计计算得到弹性应变的分布图，本文使用2.2节中复合后的位移估计，这有别于在应变图像上进行复合［13］。应变和位移之间的关系，可以用梯度来描述，梯度的处理过程可以用数字差分滤波器。定义SG-I DD［15］如下

式中：x（n）是在点n处的位移估计；y（n）是在点n处的应变估计；M是最大差分距离。

3 算法验证

评价超声图像质量的另外一个指标是对比度信噪比（CNRe）［13］

式中：ut和ub分别为目标囊肿和背景组织图像的平均值；σt和σb分别为它们的标准差。

试验所使用的超声成像系统是某公司研发的超声机，采用的探头是信号为SA5L38B的128阵元的线阵探头，阵元之间的间隙为0.3mm。探头的中心频率是5MHz，75%分数阶带宽。试验过程没有安装额外的装置，使用徒手准静态压缩组织的方式进行数据采集和成像。在压缩和释放过程中，尽量保持探头匀速轴向运动，用力均匀，确保扫描在轴向同一平面。

实验中，体模（Phantom）采用美国某公司研发的专用弹性成像研究的弹性体模。体模中，在背景介质中放有4个直径为10mm和4个直径为20 mm的球形包容物，球形包容物和背景是等回声的，在B模式图像中不能区分。背景介质的弹性模量是30kPa，由仿人体组织的均匀弹性各向异性介质组成。两种直径的4个球形包容物均有不同的弹性模量，它们分别是8kPa（TypeⅠ，软），18kPa（Type Ⅱ，软），44kPa（Type Ⅲ，硬）和63kPa（TypeⅣ，硬），如图1所示。本文成像选择的是直径为10mm且弹性模量是63kPa的硬包容物。

在试验中，使用5MHz的线性探头交替发射-5，0，＋5角度的扫描线，手动选择连续6帧来计算3个角度的位移，从而进行复合和求出复合后的应变图像。在试验中，窗口大小为10，窗口重叠率为0.75，牛顿迭代次数为2，低通滤波器参数M为7。

图2显示了弹性成像中位移估计结果图。图2（a）是0角度位移图；图2（b）是＋5角度位移图；图2（c）是-5角度位移图；图2（d）是位移复合后位移图。图3显示了弹性成像中应变结果图。图3（a）是0角度应变图；图3（b）是＋5角度应变图；图3（c）是-5角度应变图；图3（d）是应变复合后的图像［13］；图3（e）是位移复合后应变图。可以看出，复合后的弹性应变图比复合前任何一个角度的应变图都要平滑，伪影噪声得到一定抑制。

图1 Model 049弹性体模内部结构示意图Fig.1 Structure figure of Model 049phantom

图2 轴向弹性成像位移图Fig.2 Displacement of axial strain imaging

在表1中，给出了超声弹性成像质量量化表，其中独立的方框（图3（d））是100×250个采样点的窗口，用来计算信噪比（SNRe）；带椭圆的方框用来计算对比度信噪比（CNRe），这里方框里的椭圆区域是目标囊肿，方框里的其余区域为组织背景。从表1中可以看出，复合后的信噪比和对比度信噪比均有提高，且本文位移复合校正后应变图像的信噪比和对比度信噪比，均比在应变图像上进行复合［13］要高，这一实验结果也证明了论文提出的方法的有效性。

图3 轴向弹性成像应变图Fig.3 Strain image of axial strain imaging

表1 各角度子弹性图像与复合后弹性图像的信噪比和对比度信噪比Table 1 SNRe and CNRe of sub-elastogram and compounded images

另外，应变复合算法［13］和本文的位移复合算法在步骤上进行比较：应变复合算法先得到各个角度的位移估计，然后对各个角度的位移图像进行数字差分滤波来求得各个角度的应变图像，最后对各个角度的应变图像进行复合得到抑噪后的应变图像；而位移复合算法是得到各个角度的位移估计后，对多个位移图像进行复合，然后只需进行一次数字差分滤波即可得到最后的抑噪应变图像。很显然，本文的位移复合算法只需进行一次数字差分滤波，即可得到抑噪的应变图像，而应变复合算法需要对每个角度的位移估计进行数字差分滤波。在实时系统中，本文的位移复合算法比应变复合算法占有优势。

4 结束语

本文描述了位移空间角度复合技术在超声弹性成像中的应用，这种方法在不减小分辨率的情况下减少伪影噪声的影响，提高了弹性成像图像质量。本文在原理上说明了这种方法的正确性。在某公司研发的超声系统上，采用美国某公司研发弹性体模进行试验，实验结果也表明了该方法的有效性，信噪比和对比度信噪比的提高，也显示了这一点。另外，在实时系统中，本文的位移复合算法比应变复合算法占有优势。

［1］Konofagou E，Ophir J.Myocardial elastography—A feasibility study［J］.Ultrasound in Medicine and Biology，2002，28（4）：1113-1124.

［2］许慰玲，沈民奋，杨金耀.基于DTCWT方向信息的超声图像斑点噪声消除［J］.数据采集与处理，2011，26（6）：671-675.

Xu Weiling，Shen Minfen，Yang Jinyao.Speckle reduction for ultrasound image using edge directions in DTCWT domain［J］.Journal of Data Acquisition and Processing，2011，26（6）：671-675.

［3］Alam S K，Ophir J，Konofagou E E.An adaptive strain estimator for elastography［J］.IEEE Trans Ultrason Ferroel Freq Cont，1998，45（2）：461-472.

［4］Techavipoo U，Varghese T.Wavelet de-noising of displacement estimates in elastography［J］.Ultrason Med Biol，2004，30（4）：477-491.

［5］Lindop J E，Treece G M，Gee A H，et al.Estimation of displacement location for enhanced strain imaging［J］.IEEE Trans Ultrason Ferroel Freq Cont，2007，54（9）：1751-1771.

［6］Dantas R G，Costa E T.Ultrasound speckle reduction using modified Gabor fillters［J］.IEEE Trans Ultrason Ferroel Freq Cont，2007，54（3）：530-538.

［7］Liu P，Liu D C.Oriented demodulation and frequency splitting for directive filtering based compounding［C］//2008IEEE Ultrasonics Symposium.Beijing：［s.n.］，2008.

［8］Liu P，Liu D C.Directive filtering schemes for frequency compounding in ultrasound speckle reduction［C］//Proceedings of 2008International Pre-olympic Congress on Computer Science.Nanjing：［s.n.］，2008：227-231.

［9］Tanter M，Bercoff J，Sandrin L，et al.Ultrafast compound imaging for 2-D motion vector estimation：Application to transient elastography［J］.IEEE Trans Ultrason Ferroel Freq Cont，2002，49（3）：1363-1374.

［10］Techavipoo U，Chen Q，Varghese T，et al.Noise reduction using spatial-angular compounding for elastography［J］.IEEE Trans Ultrason Ferroel Freq Cont，2004，51（5）：510-520.

［11］Chen Q，Gerig A L，Techavipoo U，et al.Correlation of RF signals during angular compounding［J］.IEEE Trans Ultrason Ferroel Freq Cont，2005，52（6）：961-970.

［12］Rao M，Varghese T.Correlation analysis for angular compounding in strain imaging［J］.IEEE Trans Ultrason Ferroel Freq Cont，2007，54（9）：1903-1907.

［13］Zhang Kun，Liu P，Liu D C.Strain image with spatial angular compounding［C］//International Conference on Bioinformatics and Biomedical Enginerring.Chengdu：［s.n.］，2010.

［14］Andreas Pesavento，Christian Perrey，Martin Krueger，et al.A time-efficient and accurate strain estimation concept for ultrasonic elastography using iterative phase zero estimation［J］.IEEE Transactions on Ultrasonics，Ferroelectrics and Frequency Control，1999，46（5）：1057-1067.

［15］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）：1429-1434.