用于活体温度评估的非线性超声热应变模型∗

孙彪 郭霞生 屠娟,2 章东,2

(1 南京大学物理学院 南京 210093)

(2 中国科学院声学研究所 声场声信息国家重点实验室 北京 100190)

0 引言

在肿瘤和高血压等疾病的治疗研究中，热疗被认为是具有重要价值的微创或非侵入性技术[1−2]。考虑到热疗的安全性和有效性，在治疗过程中对生物组织的温度变化进行监控很有必要。其中，超声辐照肾周脂肪治疗高血压过程中的温度监控是本文作者团队近年来关注的重点。

与其他温度成像方法相比，超声测温技术具有无电离辐射、灵活、易操作、低成本、便携和兼容性好等优势，被认为是一种有前景的温度评估方法。其原理主要依赖于温度变化引起超声传播过程的变化，具体包括散射点信号强度变化[3−4]、回波信号的谱峰偏移[5]等。其中，基于超声回波偏移的热应变测温技术受到广泛关注[6−8]。

超声热应变测温基于组织热膨胀和声速随温度线性变化的假设。当组织温度发生变化时，由于热膨胀和声速变化引发回波时移，该时移的轴向梯度(即热应变)与温度变化量存在线性关系，可用于温度估计[9]。Simon等[10]对加热的凝胶进行实验，在4.22℃温度变化范围内实现的评估误差小于0.5℃。Liu等[11]对猪心脏组织加热，借助GPU运算实现了温度的实时监控。但当组织温度高于50℃(高于体温13℃)时，组织的性质通常发生较大变化，该线性温度评估模型不再适用。即对于生物组织，线性的声速变化和热膨胀仅在有限的温度范围适用[9]。因此，对于更大范围的温度评估，线性热应变测温模型需要改进。

对活体的温度变化进行评估时，热应变远小于组织运动导致的机械应变，容易被后者掩盖[9]。因此，有效抑制活体生理运动同样是临床应用之前要解决的关键问题。Daniels等[12]考虑活体呼吸的周期性，采取动态帧选取的方法来抑制呼吸运动，进而降低猪肾脏中温度评估的误差。Bayat等[13]用自适应运动补偿的方法来消除呼吸脉搏对活体超声评估的干扰。对于上述温度评估方法，报道的温度评估误差为4℃左右，并不能满足实际的要求。在之前的工作中，动态帧选取和自适应滤波结合的常系数温度估计(Constant coefficient model with dynamic frame selection and adaptive filtering，简称CDA)算法可以很好地抑制活体中生理运动对温度评估的干扰[14]。

本文进一步发展超声热应变测温理论，提出一种非线性热应变测温模型，结合之前发展的运动抑制算法对活体温度变化进行评估。以猪的肾周脂肪作为目标区域，进行温度评估实验。在3种不同温度变化范围的实验中，以最大相关系数动态帧选取和自适应滤波相结合的方法对活体猪的呼吸心跳的干扰进行抑制，并取得了很好的效果。本算法将推动超声测温在临床热疗中的应用，可显著提高热疗的安全性和有效性。

1 温度评估模型和算法

1.1 非线性热应变模型

超声成像中，轴向深度z处的散射子引起的时间延迟可以描述为声速和温度的函数：

其中，T0(ξ)是轴向位置ξ处的温度，c[ξ,T0(ξ)]为轴向位置ξ处温度为T0(ξ)时的声速。

加热时，目标区域组织发生热膨胀。若温度变化范围较大，热膨胀引起的体积变化与温度之间的关系是非线性的[11]，声路径上两个散射点之间的轴向距离l随温度变化的关系近似描述为

进一步，认为声速也是温度的非线性函数[11]，表示为

当组织温度从T0变化为T时，z处回波的时移(回波时延的变化量)为

其中，tf(z)和ti(z)分别为T和T0时刻的回波时延。将方程(1)和方程(4)代入方程(5)，得到

将式(6)对z微分，并忽略高阶小量，得到

进一步考虑到∣β2(z)δ2T(z)∣≪1+β1(z)δT(z)，得到

假设加热区的组织为单一组分，即区域内各处具有相同的基线温度，同时考虑到[β1(z)δT(z)]2≪1，由式(8)得到

其中，k1=α2(z)−β2(z)+β1(z)[α1(z)−β1(z)]和k2=α1(z)−β1(z)为非线性温度评估模型的两个测温系数。

通过求解方程(9)，得到时移与温度变化之间的关系：

此即为文中采用的非线性温度评估模型。如忽略其中的二阶小量，则方程(10)化简为

与文献[10,15]中描述的线性热应变测温模型一致。由此可见，对于小范围的温度变化的情况，式(10)中的二阶变化量很小。因此，非线性模型也适用于小范围的温度变化情况。

1.2 CDA算法

与线性模型的温度评估算法类似[10,15]，非线性温度评估模型也是基于组织回波偏移估计获得具体的温度估计。而对于活体组织的回波偏移估计，应考虑活体生理运动的干扰。对本文考虑的靶区组织肾周脂肪而言，其生理运动主要由呼吸和心跳运动共同引起[9]。本文采用CDA算法来对二维温度图像中的噪声信号进行抑制，该算法包含动态帧选取、热应变计算和自适应滤波3部分。动态帧选取算法通过计算超声射频(Radio frequency,RF)图像间的归一化互相关系数，并据此选择运动状态相似的图像，以抑制短周期噪声(心跳)引起的干扰。自适应滤波算法则基于呼吸运动的全局性，用自适应滤波器(Adaptive filter,AF)对长周期的呼吸干扰进行运动补偿。CDA算法流程图如图1所示，具体方法可以参考文献[14]。该算法大致按如下3个步骤进行：

图1 算法流程Fig.1 The f low diagram of the CDA algorithm

(1)对升温过程的超声RF图像，以心跳周期为选帧范围，使用动态帧选取算法挑选出最大相关性的RF图像序列I1。在升温前的RF图像中，对I1逐帧挑选相关性最好的图像作为非升温的参照图像序列I2。

(2)基于选取的RF图像I1、I2，通过散斑追踪计算RF图像中的散射点位移分布，再依次应用轴向数字差分滤波器和轴向、横向的低通滤波得到热应变分布S1及S2。

(3)设计归一化最小均方(Normalized least mean square,NLMS)自适应滤波器，通过S2所包含的机械应变信息，以步长0.1迭代以确定滤波器的系数。然后将滤波器应用于含噪声的热应变图像S1，得到降噪后的热应变分布，若已知温度系数k1和k2，则可得到温度估计结果。

1.3 实验设置

实验前将猪禁食12 h，完成清洗后，用异丙酚和七氟醚对猪进行麻醉。使用呼吸机对处于深度睡眠状态下猪进行辅助呼吸。实验装置示意图如图2所示。对猪备皮，使用B超成像探头对肾周脂肪定位，在平行于脊柱的位置使用高频电刀开出长约5 cm的创口。使用128阵元的L12-5线阵探头和便携式超声设备(定制，珠海医凯超声设备有限公司)对目标区域成像并采集RF数据。在B超成像画面的引导下，将通过模具固定的微波消融针(ECO-200G，南京亿高微波系统工程有限公司)和热电偶(HYP2,OMEGA)沿成像平面的法向插入肾周脂肪。微波消融针的发射单元和热电偶的感温区域均在尖端，并位于成像平面上。微波消融针和热电偶的间隔为1.3 mm。动物实验准备、麻醉和手术得到了江苏省人民医院的支持，且遵循美国国家卫生研究院《实验动物的护理和使用指南》(2011)，并经南京总医院伦理委员会批准(批准号20140510)。

图2 实验装置Fig.2 The experimental setup

实验中，在加热前采集至少两个呼吸周期长度的非升温RF信号。将微波消融针的发射功率设定在不同水平，在每个加热过程，同步获得热电偶测量的温度和B超图像的RF数据。热电偶采样率20 Hz，RF数据采样率40 MHz，帧率32 Hz。每帧RF图像尺寸为2608(轴向)×128(横向)像素。

2 结果

2.1 二维温度分布

根据多次微波加热活体猪肾周脂肪的实验，采集37℃～50℃升温过程的RF数据用以计算热应变，根据热电偶在RF图像中的位置提取热电偶处的热应变估计值，将其和热电偶测得的温度曲线进行最小二乘拟合确定非线性模型温度系数k1、k2。由此确定的非线性系数(二次项系数)为k1=0.0045，线性温度系数(一次项系数)k2=0.0930。一组超声CDA算法评估的结果如图3所示，其中给出了靶区在4个不同时刻的二维温度分布图像。由图3可知，随着时间的增加，组织温度逐渐升高，温升区域范围不断扩大。对于真实生物组织传热过程，无论是热传导方程计算的模拟结果，还是红外成像记录的真实温度扩散过程，加热点的温度变化都是由低到高，温升范围都是由中间向四周扩散。因此本文的超声温度评估算法的温度和真实的物理过程具有很好的一致性。

图3 二维温度变化序列Fig.3 Two-dimensional temperature change sequence diagram

2.2 两种不同模型的对比

图4分别给出了对应3种不同的加热功率(对应不同的温升范围)下，传统线性模型和本文非线性模型的温度评估结果，其中线性模型采用了式(11)描述的时移-温度关系。对于小范围的温度变化(12℃，满足线性模型温升的范围)，线性模型和非线性模型的最大温度偏差都在2℃左右，在这种情况下，两者的温度评估精度是相当的；当温度变化范围提高到17℃时，线性模型温度评估的最大评估误差为3.04℃，而非线性模型的最大误差为2.03℃。相较而言，非线性模型具有更高的温度评估精度。随着微波功率继续增加，组织最大温升达到25℃，线性模型的温度评估最大误差达到5.08℃，而非线性模型的温度评估误差最大为2.34℃，两种模型的性能差距进一步扩大。因此，对于较大温升范围过程中的温度评估，非线性模型更具有优势。

图4 线性模型与非线性模型升温曲线对比Fig.4 Temperature curves of linear model and nonlinear model

3 讨论

对于超声评估活体温度的算法，温度评估模型是影响温度评估准确性的重要因素。研究表明，无论是活体还是离体，只有在温度变化较小的情况，线性热应变模型的评估结果才会比较准确[8,16]。当生物体温度高于50℃，基于回波偏移的线性热应变模型无法对生物组织的升温过程进行有效的评估[10]。对于大范围的温度变化，生物组织的体积和声速随温度的变化存在非线性关系[17]。因此，本文提出将热膨胀和声速描述为温度的非线性函数，为较大范围温度变化下的温度评估问题提供了合适的解决方案。

需要指出的是，本文的模型仅考虑了轴向的回波时移，而实际组织受热时温度场是全向扩散的。事实上，当组织发生温度变化时，由于热膨胀和声速变化的共同作用产生轴向的回波偏移。随着加热过程的进行，RF图像在横向上仅存在由于热膨胀导致的图像变形，在轴向上存在由于热膨胀和声速变化导致的图像变化。但重要的是，热膨胀形成的回波偏移远小于声速变化导致的偏移[9]，因此文中所得热应变主要源于声速变化的贡献。尽管如此，如能同时从横向考虑组织受热导致的图像变形，是完全可能进一步提高精度的。但这通常很困难，因为B超图像的横向分辨率一般远低于轴向。

非线性热应变模型可在较大的温度变化范围内保证温度评估的准确性。本文的实验结果表明，对于25℃的温升，非线性模型的误差精度基本可以满足临床的测温需求。但对于活体脂肪组织，当热疗使其温度达到70℃时，组织可能发生变性，非线性模型也就不再适用。如图5所示，当温度变化达到40℃以上时(组织温度大于70℃)，虽然非线性模型的温度评估误差仍小于线性模型的评估误差，但是其最大估计偏差已大于4℃，难以满足临床的要求。这可能与脂肪组织的严重变性(如液化)有关[18]。因此，对于变性的组织，无论线性还是非线性模型温度评估模型都不再适用。鉴于此，Liu等[19]提出了自适应系数和红外成像相结合的算法对消融治疗中的温度(治疗温度超过90℃)进行评估，取得了一定的效果。但该方法须借助红外成像，难以应用在深部活体组织的测温中，因而其应用受到限制。

图5 温度变化40℃的线性和非线性模型温度估计Fig.5 Linear and non-linear model temperature estimation for temperature change of 40℃

对超声热应变测温在活体中的应用，生理运动的影响必须考虑在内。此前有很多关于运动补偿提高超声测温精度的研究。例如，基于在加热之前采集一系列参考帧的运动补偿算法[8,20]，利用靶区外的组织变形进行运动校正[13]等。这些算法在离体实验中显示出令人满意的结果，但在活体中的结果并不理想。本文基于心跳运动的准周期性，通过动态帧的选取来降低心跳周期性的运动干扰。此外，对于活体呼吸(血流等)带来的组织旋转和形变[8]，简单的平滑滤波和动态帧选取法并不足以解决运动步长问题，本文中采用自适应滤波得到了较好的结果。综合非线性模型和活体运动抑制算法，本文中对25℃温升的脂肪组织将测温误差控制在2.5℃以内，取得了令人满意的效果。

4 结论

本文利用生物组织升温过程中物理量变化的非线性，提出了用于活体温度评估的非线性热应变模型，对大范围的温度变化进行了有效的评估。采用动态选帧算法和自适应滤波相结合的方法来抑制活体的呼吸及心跳对温度评估算法的干扰，并取得了很好的效果。实验中以活体猪肾周脂肪为对象，对于小于25℃的温升过程，测温误差不高于2.5℃。但对于更大范围的温度变化(大于40℃)，非线性温度评估效果虽然比线性模型有明显的提升，但还无法满足临床的需求。本文提出的非线性模型提升了超声热应变测温的适用温度范围，拓展了超声温度评估技术的应用前景，有助于提高热疗相关的治疗过程的安全性。