一种超声弹性成像试验数据的快速仿真算法

崔少国 LIU Pauls 刘东权(四川大学计算机学院，四川成都 610065)

2(郧阳医学院计算机中心，湖北十堰 442000)

3(声泰特成都科技有限公司，四川成都 610041)

引言

超声弹性成像是一种生物组织成像的新模式，是当前超声医学工程学领域的前沿和热点。该成像技术可以探测出人体组织与病理变化相关的硬度改变，对肿瘤(癌症)的早期检测和诊断具有重要意义。

弹性成像的研究自1991 年被Ophir 等人提出以来，19 年中取得了快速发展，在临床中已逐渐显示出重要诊断价值，比如在乳腺癌、前列腺癌、甲状腺结节、肝纤维化、动脉粥样硬化等的检测中取得了重要的应用［1-4］。其成像的基本原理是:对组织施加一个动态或者准静态的压力，组织将产生一个响应。在相同的应力下硬的组织形变较小，而软的组织形变较大。使用信号处理和图像处理相关技术计算出组织的应变大小，进而重构出组织的弹性模量。对组织的应变或弹性模量进行成像就是弹性成像［2］。

1 弹性成像试验数据仿真

在弹性成像算法设计、性能评估和成像参数选择的研究中，由于从实际超声系统中获取的数据具有事先不可预知的噪声或其它干扰因素，所以一种新算法的验证通常都需要在仿真的数据上进行预实验，然后再进行体模实验和在体实验。弹性成像试验数据的仿真在超声弹性成像研究中具有重要的意义。

当前，超声弹性成像试验数据的仿真一般采用有限元分析(FEA)和Field Ⅱ声场模拟程序来实现。这种方法把组织仿真成在某种基质中分布着与基质声学特性有显著差异的散射子，各散射子的位置服从均匀随机分布，而其反射强度服从高斯正态分布，从而建立计算机人体模型［5］(以下简称计算机体模)。对于压缩前回波信号，依据线性声学理论，使用Field Ⅱ程序在仿真的计算机体模上生成。对于压缩后回波信号，先使用有限元分析出各散射子在一定应力作用下运动到达的新位置，使用它对计算机体模中散射子位置进行重新调整，然后再使用Field Ⅱ程序生成压缩后回波信号。该仿真方法与实际成像过程很相似，结果也较准确，但仿真数据生成速度较慢。一般来说，仿真一帧射频回波信号需要6 小时左右。在弹性成像研究中，通常需要仿真不同组织模型、不同应变量和不同信噪比的序列回波数据帧，数据量较大，使用该方法来仿真需要较长时间。

2 一种快速仿真算法

笔者在弹性成像研究过程中探讨出了一种快速仿真弹性成像试验数据的算法，该算法主要包括以下三步:参考帧数据生成，压缩后回波数据生成和数据加噪处理。

2.1 参考帧数据生成

参考帧就是压缩过程中采集的序列帧中的第一帧。对参考帧回波数据的仿真可以使用两种方法:基于计算机体模和Field Ⅱ的方法和使用超声线性系统理论的方法。实际应用中可以选择其中任一种方法。

2.1.1 使用计算机体模和Field Ⅱ声场模拟程序

计算机体模使用文献［4］中介绍的方法生成，采用组织分立散射模型，各散射子的位置服从均匀随机分布，而其反射强度服从高斯正态分布。声场计算使用Tupholme 和Stepanishen 提出的算法。

探头在声场中位置为→r 处产生的声压被描述为空间冲击响应［6］，表示如下:

这里，δ 是冲击函数，c 是超声波在组织中传播速度，s 是探头的表面积。那么探头发射波的声压场为［7］

这里Vn(t)是探头的表面速度，ρ 是组织密度，“* ”代表卷积。

因此，从探头接收到的超声回波信号可模拟为［8］

式中，Vpe是回波冲击响应，fm项是由于组织的密度变化和声速扰动引起的散射信号，hpe是脉冲回波的冲击响应，分别为

2.1.2 使用超声线性系统理论

使用超声简单线性系统理论，参考帧回波信号可以用如下的公式产生:

这里f(x，y)是射频回波信号，p(x，y)是系统点扩散函数(PSF)，h(x，y)是散射函数，散射子的位置服从均匀随机分布，散射子的散射强度服从高斯分布，“* ”代表卷积。二维点扩散函数p(x，y)在轴向(x方向)上可以模拟成高斯调制的余弦函数，在侧向(y 方向)上可以模似成一个高斯包络，表示如下

其中，ω0是探头的中心角频率，σx和σy分别是轴向和侧向的标准方差。散射函数h(x，y)可以用服从均匀分布的随机函数生成。

2.2 压缩后回波数据生成

如果组织内各区域弹性模量一致，当在相同应力作用下，组织各部分发生相同的应变，这种称为均匀应变模型。此模型主要用来仿真均匀的无肿瘤人体组织，用来研究弹性成像算法的信噪比SNRe。如果组织中某些区域弹性模量与周围区域不同，当在均匀应力作用下，不同弹性模量的组织将产生不同的应变，这种模型称为非均匀应变模型。此模型主要用来仿真均匀组织中存在一硬癌块或软囊肿，用来研究弹性成像算法的对比度(contrast)或对比度噪声比CNRc。

算法1 是仿真组织压缩后回波信号的快速生成算法。设均匀组织各点的弹性模量m 是常数C，非均匀组织中各点的弹性模量m 是e(i，j)(i，j 是组织区域二维空间坐标)，模量大小随空间位置坐标(i，j)而变化。设组织为平面应变状态，组织的各部分均匀受力，各点的应力为F。根据弹性力学及胡克定律，各点的应变

根据线弹性体中位移连续性原则，则各点的位移d(i，j)为

设RF_pre 和RF_post 是组织压缩前和压缩后的回波RF 信号，IQ_pre 和IQ_post 是其对应的基带I/Q信号，则算法可以简单描述如下:

算法1:组织压缩后回波信号快速生成算法

1)设计预仿真组织的弹性模量分布图。对于均匀组织各区域设成相同的模量，即m(i，j)=C(C是常数);对于非均匀组织根据预仿真的模型，不同的空间位置设成不同的弹性模量，即

2)根据式(7)由弹性模量图计算组织应变分布图。

3)根据式(8)由应变图计算各点在应力作用下的位移。

4)根据位移图求出插值位置图，使用式p(i，j)= d(i，j)+i 来完成。

5)使用式(9)对压缩前信号在新的位置上进行插值，生成压缩后RF 回波信号RF_post。

6)使用希尔伯特变换生成复信号，Analytic_post=Hilbert(RF_post)。

7)生成基带I/Q 信号，IQ_post = Analytic_post×exp(j ×w0 ×p)。

8)如果还需生成压缩过程中的下一帧数据，则RF_pre =RF_post，转4);否则结束。

这里插值方法可以选用计算代价较小的线性插值，也可以选用计算代价较大但较精确的三次样条插值。这里采用线性插值技术，插值函数如下:

2.3 数据加噪处理

为了仿真真实超声系统中获取的不同信噪比的回波数据，需要对以上产生的信号进行加噪处理。

根据信噪比的定义

式中，Ps是信号的功率，Pu是噪声的功率，则

因此，所需要加的高斯白噪声为

这里，Rand(m，n)是产生m 行n 列均值为0 方差为1 的高斯随机信号值函数。加噪算法可简单描述如下:

算法2:数据加噪算法

1)设定数据的信噪比SNR，计算原信号的功率Ps。

2)根据式(10)计算需加的噪声功率Pu。

3)根据式(11)计算所需的高斯白噪声u(m，n)。

4)对仿真的回波数据f(i，j)加噪处理，fnoise(i，j)= f(i，j)+ u(i，j);

5)结束。

3 算法验证

为了验证以上提出的快速仿真算法的正确性和高效性，笔者编程实现了上述仿真算法，并生成了一系列数据集。然后选用一标准超声弹性成像程序，在仿真的数据集上运行，以验证能否产生预期的应变图及成像效果。

3.1 参考帧生成

实验对两种生成参考帧的方法都进行了验证。计算机体模是仿真40 mm ×40 mm ×10 mm 并均匀分布有200 000 个散射子的人体组织，探头中心频率5 MHz，系统采样频率是20 MHz，使用Field Ⅱ和式(1)～式(4)生成128 条A-lines，帧大小128 像素×1039 像素。使用超声线性系统理论产生参考帧时，使用式(5)和式(6)。

3.2 应变模型设计及压缩后回波数据生成

本实验中预仿真两种组织应变模型:均匀组织应变模型和软背景中包含一球形硬包容物非均匀组织应变模型，模型的弹性模量分布如图1 所示(弹性模量用256 级灰度值表示，平均模量映射成128)。图1(a)是均匀组织弹性模量分布图，组织内各点的弹性模均是30 kPa。图1(b)是非均匀组织应变模型的弹性模量分布图，图中硬物的直径为20 mm，其弹性模是300 kPa，背景的弹性模量是30 kPa，硬物是背景硬度的10 倍。使用算法1 对每种模型产生20 帧应变序列。

实验针对均匀组织应变模型仿真不同整体应变的回波序列，包括0.01%，0.1%，0.5%，1%，2%，3%，4%，5%，用来验证算法产生的数据是否满足弹性成像应变滤波器。对非均匀模型仿真了背景应变为1%、硬物应变为0.1%的数据。

图1 预仿真的组织模型弹性模量分布图。(a)均匀组织;(b)含一硬包容物的非均匀组织Fig.1 Simulated tissue modulus distributing images.(a)uniform model;(b)one-lesion model

3.3 数据加噪

使用算法2 对产生的各帧仿真数据添加均值为零的不同能量高斯白噪声，使其信噪比分别为5、10、15、20、25、30、35、40、45、50 dB。

3.4 成像算法的选择

实验中选用英国剑桥大学Joel. E Lindop 于2008 年提出的基于相位的WPS(weighted phase separation)实时弹性成像算法［9］来验证仿真的数据，该算法已被验证是正确可行的，能在实际超声系统上进行弹性成像。位移估计使用式(12)迭代进行

其中，WA(t)= Ax(t)Ay(t + dk)，

应变计算使用低通数字差分滤波器SG-I DD［10］，方法为

式中，s 为应变，Δk= d(n + k)- d(n - k)

3.5 时间效率和数据仿真效果检测方法

为了对新算法产生弹性成像试验数据的时间效率和仿真效果进行检测，本研究中使用了三种方法产生弹性成像试验数据并进行比较。方法1 是:参考帧使用计算机体模和Field Ⅱ生成，压缩后的回波数据先使用有限元分析，进而对体模中的散射子的位置进行调整，然后使用Field Ⅱ产生回波，数据加噪使用算法2;方法2 是:参考帧使用计算机体模和Field Ⅱ生成，压缩后的回波数据使用提出的算法1、算法2 生成;方法3 是:参考帧使用简单线性系统理论卷积生成，压缩后的回波数据使用提出的算法1、算法2 生成。运行算法程序所使用的计算 机 配 置 是:Pentium (R) 4，CPU 2.4 GHz，RAM 1GB。

弹性成像试验数据效果的比较，主要以在该数据集上产生的对应弹性图像的信噪比来衡量。信噪比越高说明仿真与实际越接近，仿真越准确。该项实验使用均匀应变模型的数据，数据没有进行加噪处理。

4 结果

4.1 应变图像及性能

图2 是在仿真的均匀组织应变模型的数据集上使用WPS 算法产生的应变图。成像中窗长为60 个采样点，窗口重叠率为0.75，迭代3 次。运动追踪在1D 轴向进行，轴向应变计算采用M =4 的低通数字差分滤波器SG-I DD 生成。图2(a)是由仿真的1%应变的无噪声数据集产生，图2(b)是由仿真的1%应变且信噪比为10 dB 的数据集产生。应变映射成256 级灰度值，平均应变映射成128(下同)。

图2 仿真的均匀组织应变模型数据集上产生的应变图像。(a)无噪声;(b)SNR =10 dBFig. 2 Elastography obtained from the simulated echoes of uniform strain profile. (a)no noise;(b)SNR =10 dB

图3 是在仿真的非均匀组织应变模型(含一硬包容物)的数据集上使用WPS 算法产生的应变图。成像参数与方法同图2。图3(a)是由仿真的1%整体应变的无噪声数据集产生，图3(b)是由仿真的1%整体应变且信噪比为10 dB 的数据集产生。

图3 仿真的非均匀应变模型数据集上产生的应变图像。(a)无噪声;(b)SNR =10 dBFig. 3 Elastography obtained from the simulated datum of one-lesion strain profile. (a)no noise;(b)SNR =10 dB

从图2 和图3 可以看出，在用提出的快速算法所仿真的数据集上可以生成预期的应变图像。图2产生的应变图像与预仿真的均匀组织应变模型相符，图3 也与预仿真的的非均匀组织应变模型相符(较硬的圆形区域里应变约为背景的1/10)。在无或有噪声的仿真数据上采用相同的成像算法，和成像参数所生成的应变图有明显差别，在有噪声的仿真数据上产生的应变图噪声较大，这与弹性成像的实际情况相符。

图4 显示了应变图像中最中间的一条线(轴向)上的应变数据。从图4 中也可以看出，在仿真的无噪声数据集和加噪声数据集上产生的应变数据是有效的，与弹性成像理论相符。

图5 显示了在仿真的各种应变数据集上生成的应变图像的信噪比。从图中可以看出，信噪比曲线满足应变滤波器特性，类似于一带通滤波器。而且在低信噪比(10 dB)的数据集上生成的应变图像信噪比明显低于无噪声数据集。

4.2 时间效率与仿真效果比较

各种产生弹性成像试验数据方法的时间性能如表1。从表1 中可以看出，采用基于计算机体模和Field II 方法，产生1 帧形变数据大约需要6.3 h，而采用新方法(算法1)，产生1 帧相同的数据只需要约35 s，速度大大提高。

表1 3 种方法产生弹性成像试验数据所需时间(h/帧)Tab.1 Comparison time of three methods for producing strain echoes(h/frame)

各种方法产生的弹性成像试验数据的效果见表2。表2 显示，在新方法(方法3)产生的试验数据集上所获得的应变图像的信噪比略低于传统的基于有限元和Field II 方法(方法1)，尤其在大应变时较为明显。方法2 和方法3 的信噪比几乎没有差别。

表2 3 种方法数据集上产生的弹性图像信噪比比较Tab.2 Comparison of SNR obtained from data by using three methods

图4 一条应变数据线。(a)均匀应变;(b)10 dB 信噪比数据的均匀应变;(c)含一硬包容物的应变;(d)10 dB 信噪比数据且含一硬包容物的应变Fig.4 A strain line. (a)uniform strain;(b)uniform strain with a 10 dB SNR;(c)one-lesion strain;(d)onelesion strain with a 10 dB SNR

5 讨论和结论

图5 各种不同应变的仿真数据集上产生的应变图像信噪比Fig.5 SNRe of strain images obtained from various strain simulation data

本研究提出的一种超声弹性成像试验数据的仿真算法是可行的、高效的，能灵活快速地生成各种应变、各种信噪比、各种组织模型的仿真超声回波数据。从成像的可行性来看，使用新算法产生的试验数据能够生成预期的均匀组织的弹性图像和含有硬包容物组织的弹性图像;从时间效率来看，新算法产生弹性成像试验数据比Field Ⅱ方法快得多;从产生的试验数据效果来看，与使用Field Ⅱ方法产生的数据相比除了在较大应变时略低外，其它几乎无差别。

新方法的提出对弹性成像仿真研究具有重要意义，为弹性成像算法验证、成像各参数的选择等研究提供了重要手段，为超声弹性成像试验数据仿真提供了一种新的高效快捷、方便可行的方法。该方法也能对任何复杂组织进行仿真，此时只需要设计好复杂组织的弹性量分布图即可。该方法的局限性是对大应变的试验数据仿真较不准确，这主要是因为当应变较大时，使用有限元方法对散射子位置进行估计要比使用插值法精确。但在实际弹性成像中一般应用的应变不会太大(一般小于0.01)，因为太大会引起严重的解相关，使当前主流时延估计算法失效。在低应变时新算法所获得的数据效果略微偏低，但仍能实行有效的仿真任务，然而时间效率的大幅度提高是值得肯定的。另外一个局限性是此算法仅适用于1D 应变模型。对于2D 应变模型，必须考虑组织散射子在应力作用下的三维运动。2D 应变模型的快速仿真算法的研究是一个有价值的课题。

本研究提出了一种弹性成像试验数据仿真的快速算法。通过成像实验、应变数据分析、信噪比实验、时耗和数据效果比较等证实，用该算法生成的仿真数据与预期目标相符，并且生成数据的速度大大提高。该算法的提出为弹性成像算法的验证、性能的比较、参数优化、成像效果的评价提供了有力工具，也为缺少真实超声实验环境的研究人员进行弹性成像前期研究提供了重要手段。

［1 ］ Ophir J，é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 ］ Righetti R. Elastographic characterization of HIFU-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 ］ Jensen JA. Computer phantoms for simulating ultrasound B-mode and CFM images［J］. Acoustical Imaging Symposium，1997，23:75 -80.

［6 ］ Jensen JA. Calculation of pressure fields from arbitrarily shaped，apodized，and excited ultrasound transducers ［J ］. IEEE Transactions on Ultrasonics， Ferroelectrics， and Frequency Control，1992，39(5):262 –267.

［7 ］ Jensen JA. Fast Simulation of Ultrasound Images［J］. IEEE Ultrasonics Symposium，2000，16:1721 -1724.

［8 ］ Jensen JA. Field:A program for simulating ultrasound systems［J］. Med. Biol. Eng. Comp.，1996，4:351 -353.

［9 ］ 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.

［10］ 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.