热疗消融具有副作用小、无创或微创等优点，已经成为继手术、放疗、化疗和免疫疗法之后的第5大肿瘤疗法^[1]。为保障热疗消融的安全性和有效性，必须对消融过程进行实时准确的温度监测。磁共振温度成像(magnetic resonance temperature imaging，MRTI)由于具有非介入式、时间和空间分辨率较高、能够任意方向成像等优点，是当前热疗消融温度监控的最佳技术手段^[2-3]。

磁共振温度成像可以基于多种参数进行温度测量，主要包括：纵向弛豫时间T₁、横向弛豫时间T₂、水质子共振相移(proton resonance frequency shift, PRFS)等^[2]。在上述参数中，PRFS具有最好的温度线性关系，因而被最为广泛地使用。基于PRFS参数，一系列的测温方法被研发出来^[4-6]。然而，磁共振温度成像普遍面临的挑战是，测温过程面临着时间与空间分辨率、视场(field of view, FOV)大小和测温精度之间的权衡。当选用更高时间和空间分辨率的测温序列、采集更大FOV的图像时，测温的精准度和温度图像的信噪比往往难以得到保证^[7]。

针对这一问题，有学者提出了使用0范数约束进行降噪的Hybrid方法。Hybrid方法假定消融升温区域仅占FOV的极少部分，而其他大部分区域在消融过程中温度不变，从而利用升温区域的稀疏性，有效降低磁共振温度成像过程中的噪声^[6]。然而，该方法不仅对约束项的权重选择具有较高的依赖性，权重选择不当会导致重建效果变差，还会将温度升高较小的像素点误判为噪声点，从而抑制为0。在Hybrid方法基础上，Zhang等通过对热疗消融过程进行生物传热(bio-heat transfer, BHT)建模，使用了模拟温度场对Hybrid正则项的权重进行自动选择，以提高迭代重建收敛速度和磁共振温度图的信噪比^[8]。De Senneville等将生物传热过程视作非线性过程，引入了扩展Kalman滤波方法将模拟温度场和MRTI测量温度场进行结合，对磁共振信号极低区域的温度进行估测^[9]。上述2种将生物传热模型模拟温度场和磁共振测量温度场结合的方法对于提高温度图像信噪比和测温鲁棒性均有一定的促进作用。然而，上述2种基于生物传热模型的方法中，Zhang等所提出的方法基于单个像素域进行生物传热模拟^[8]，而De Senneville的方法则将生物传热过程进行了非线性近似^[9]，并使用扩展Kalman滤波方法进行结合，这些处理在一定程度上引入了生物传热模型温度场模拟的误差，从而在最终估测温度中引入了新的偏差。

针对目前存在的问题，本文首先对生物传热模型的Pennes方程进行差分形式的改写，将Pennes差分方程改写为Kalman滤波方法的状态转移矩阵，最终严格在Kalman滤波方法的框架下完成与磁共振温度图像的结合。为了探究所提出方法的准确性与鲁棒性，本文分别设计了临床数据仿真实验和琼脂仿体加热实验，并且分别使用所提出的方法和PRFS方法对实验数据进行重建。

1 原理 1.1 基于PRFS的磁共振测温

在磁共振成像过程中，多种参数会随温度发生改变，从而衍生出一系列的测温方法。水质子共振相移与温度具有最好的线性关系。该方法的主要思想是在使用梯度回波(gradient recall echo, GRE)序列采集数据的过程中，温度变化会导致采集图像的相位相较于基准图像发生线性变化。产生的相位变化与温度变化的关系可以描述为

$ \Delta \varphi=\Delta T \cdot \alpha \cdot \gamma \cdot B_{0} \cdot T_{\mathrm{E}}. $

(1)

其中：ΔT表示温度变化，α=10^-7℃^-1是比例系数，γ为旋磁比，B₀是主磁场的磁感应强度，T_E是回波时间。因此，在磁共振测温过程中，在已知初始温度T₀时，实时测量的温度为

$ T=\frac{\Delta \varphi}{\alpha \cdot \gamma \cdot B_{0} \cdot T_{\mathrm{E}}}+T_{0}. $

(2)

1.2 生物传热模型

生物传热模型的基础是Pennes于1948年提出的Pennes方程^[10]，

$ \rho c \frac{\partial T}{\partial t}=\lambda \nabla^{2} T-W_{b} c_{b}\left(T-T_{b}\right)+Q. $

(3)

其中：ρ为组织密度，kg/m³；c为组织比热容，c_b为血液比热容，J/(kg·℃)；T为温度，℃；t为时间，s；λ为导热系数，W/(m·℃)；W_b为毛细血管系统的灌注系数，kg/(m³·s)；T_b为血液温度，在人体中认为温度恒为37 ℃；Q为热源项，W/m³。在实际仿真模拟中，为了提高计算效率，会将微血管组织灌注项并入导热系数项，视作等效导热系数^[11]，因而式(3)可以改写为

$ \rho c \frac{\partial T}{\partial t}=\lambda^{\prime} \nabla^{2} T+Q. $

(4)

其中λ′是等效导热系数。考虑二维平面，将式(4)进行离散化处理，可以改写成差分形式^[12]，

$ \begin{aligned} T_{i, j}^{n+1}=& \frac{\Delta t \lambda_{i, j}^{\prime}}{\rho_{i, j} c_{i, j} \Delta x^{2}}\left(T_{i-1, j}^{n}+T_{i+1, j}^{n}+T_{i, j-1}^{n}+T_{i, j+1}^{n}\right)+\\ &\left(1-\frac{4 \Delta t \lambda_{i, j}^{\prime}}{\rho_{i, j} c_{i, j} \Delta x^{2}}\right) T_{i, j}^{n}+\frac{\Delta t \phi_{i, j}}{\rho_{i, j} c_{i, j}}. \end{aligned} $

(5)

式(5)中：下标i, j表示节点的坐标，上标则表示离散时间点。

1.3 Kalman滤波方法

Kalman滤波方法能够基于状态空间模型，从含噪声的观测数据中提取信号。在本文中，Kalman滤波方法将利用生物传热模型模拟的温度和MRTI测量的温度，递推地估计出待测区域的实际温度。

Kalman滤波方法的核心思想为：离散时间动态系统可由描述状态向量的过程方程和描述观测向量的观测方程共同表示，通过对状态向量和观测向量进行线性组合，并求解使均方误差最小的加权矩阵，可以得到最终的估测向量^[13]。

将Kalman滤波方法运用到本问题中，生物传热模型对组织温度场的模拟过程可以视作状态转移过程，而采用MRTI方法对组织温度场测量的过程则可视作观测过程。下面分别写出观测方程和过程方程的公式。

观测方程：

$ \boldsymbol{T}_{\mathrm{m}}(n)=\boldsymbol{C}(n) \boldsymbol{T}(n)+\boldsymbol{v}_{1}(n). $

(6)

M×1向量T(n)表示系统在离散时间的状态值向量，实际上无法得知其真实值；M×M矩阵C(n)是观测矩阵，能够将状态向量T(n)转化成M×1的观测向量T_m(n)；M×1向量v₁(n)是观测过程的噪声向量，一般认为是Gauss白噪声。

在使用MRTI方法测温的过程中，C(n)是M×M单位矩阵，因此观测向量T_m(n)为N×1向量，即测量温度矩阵展开得到的向量。

过程方程：

$ \boldsymbol{T}(n+1)=\boldsymbol{F}(n+1, n) \boldsymbol{T}(n)+\boldsymbol{u}(n)+\boldsymbol{v}_{2}(n). $

(7)

式(7)中：M×M矩阵F(n+1, n)为状态转移矩阵，描述了动态系统从n时刻的T(n)状态转移到T(n+1)状态；N×1向量u(n)为系统输入；N×1向量v₂(n)为对应的过程噪声向量。

现确定生物传热模型中的状态转移矩阵F。观察式(5)可知，节点(i, j)在n+1时刻的温度是n时刻该节点(i, j)本身温度以及其相邻4个节点(i-1, j)、(i+1, j)、(i, j-1)、(i, j+1)温度的线性组合加上一热源相关项，因此可通过将二维温度矩阵变换成为温度列向量的方法，将式(5)的差分方程改写成状态方程的形式。二维温度矩阵变换成温度列向量的方法如图 1所示，对应的状态方程矩阵形式为

$ \boldsymbol{T}(n+1)= \boldsymbol{F T}(n)+ \boldsymbol{u}(n). $

(8)

原始的温度图为m×m的二维矩阵，将二维的温度图按照图 1所示的方式进行变换，得到M×1的温度列向量T(n)，其中M=m×m。F为M×M的状态转移矩阵；u(n)为M×1系统输入列向量，表示外界热源输入。

为了最终确定状态转移矩阵F各位置的元素值，取矩阵F第a₀行进行分析。具体计算步骤为(见图 2)：

步骤1   首先通过式(9)确认列向量T中第a₀行各元素所在原二维矩阵的行列坐标值(i, j)；

步骤2   求解二维矩阵中与节点(i, j)相邻的4个节点的坐标，并根据式(10)计算出4个节点转换到列向量T后的行坐标值，分别记为a₁、a₂、a₃、a₄；

步骤3   根据矩阵和向量的乘法原理，T(n+1)第a₀行的值为F第a₀行向量的各元素与T(n)各元素乘积之和，因而可以通过式(5)反推出F矩阵第a₀行第a₁、a₂、a₃、a₄列的值。

使用列向量T的行坐标值计算原二维矩阵节点坐标值，

$ \begin{array}{c} x=\bmod (a, m)+1 , \\ y=\operatorname{ceil}\left(\frac{a}{m}\right). \end{array} $

(9)

式(9)中：a表示该元素所在列向量的行坐标值，x和y分别表示该元素在二维矩阵中的行和列坐标，mod为求余符号，ceil为向上取整符号，m为二维方阵的阶数。反之，通过二维矩阵节点值计算对应到温度列向量行坐标值的公式为

$ a=(x-1)·m+y. $

(10)

最后，根据式(5)可知，状态转移矩阵F中第a₀行共有5个非零元素，分别为第a₀、a₁、a₂、a₃、a₄列，其中a₁、a₂、a₃、a₄列的值相同，为了简便，对应的坐标记为(a₀, a_{1(2, 3, 4)})，因此有：

$ \left\{\begin{array}{l} \boldsymbol{F}\left(a_{0}, a_{0}\right)=1-\frac{4 \Delta t \lambda_{i, j}^{\prime}}{\rho_{i, j} c_{i, j} \Delta x^{2}} , \\ \boldsymbol{F}\left(a_{0}, a_{1(2, 3, 4)}\right)=\frac{\Delta t \lambda_{i, j}^{\prime}}{\rho_{i, j} c_{i, j} \Delta x^{2}}. \end{array}\right. $

(11)

特别地，对处于二维矩阵边界的节点(i, j)，考虑它为第1边界节点，温度恒定不变，因此对应的F(a₀, a₀)值为1。至此，矩阵F的所有元素值均计算完毕。

对于系统输入向量uⁿ，可通过式(12)求得，

$ \boldsymbol{u}^{n}=\frac{\Delta t \boldsymbol{\phi}}{\boldsymbol{\rho} \cdot \boldsymbol{c}}. $

(12)

其中：ф、ρ、c均为二维矩阵通过图 1方式转换后得到的列向量。式(12)中的列向量间的乘除法计算均为对应元素之间的乘除法。

图 3给出了使用Kalman滤波和生物传热模型模拟温度场对磁共振温度成像进行滤波并计算最佳估测温度图的完整算法。算法的完整推导参考文[12]。

图 3中，系统噪声协方差矩阵Q_w和观测噪声协方差矩阵Q_v均为M×M的主对角方阵(M为温度列向量的行数)，对角元素值分别为对应列向量中分量的方差值，在实际计算过程中，可通过预测量计算得到。

图 3所示方法得到的T(n)即为每一轮计算后的温度列向量估计值，之后将温度向量按照图 1方法进行一次反变换，即可得到估计的温度分布图。

2 实验设计 2.1 临床数据仿真实验

仿真实验数据通过对事先采集的人体前列腺相位图进行加热仿真获得。数据在GE 3T磁共振平台使用EPI序列对前列腺位置采集得到，共采集12帧图像作为基础图像，图像同时包含幅值信息和相位信息。采集时具体参数为：FOV=256 mm×256 mm×10.8 mm，Slices=3，Slice Thickness=3.6 mm，TR=160 ms，TE=11.2 ms，Flip Angle=35°，Scans=12，Time Resolution=6 s。采集到的幅值图和相位图分别如图 4a和4b所示。

之后，如图 4c所示，在前列腺给定区域使用二维Gauss热源进行模拟加热，热源总功率为800 W/m；为了获得作为真实值的对比图像，设置人体组织区域的热力学参数分别为等效导热系数为1.28 W/(m·K)、比热容为4.2×10³ J/(kg·℃)、密度为10³ kg/m³进行生物传热模型模拟，冷却水恒定为35 ℃，人体远离加热区域的温度恒定为37 ℃。模拟得到的不同时间帧的温度场通过式(1)转换成相位场, 并叠加到基础相位图上，得到仿真相位图。

图 5a是通过生物传热模型在给定热源和组织热力学参数下，模拟出的时间帧1、4、7、10时刻的温度场。图 5b为根据对应时刻温度场生成的相位图; 观察黑色箭头区域可以发现，随着温度上升，对应的相位值也在发生变化。

最后，分别使用基于PRFS的单基准相位相减法(single baseline subtraction, SBS)和本文提出的方法对仿真相位图进行重建得到温度场，并对比分析重建结果。

2.2 仿体实验

为了探究真实情况下本文方法的可适用性，使用琼脂仿体设计实验。图 6a是仿体实验设计图，其中立方体为琼脂仿体，在琼脂仿体中插入管道和加热囊(见图 6b)，用来通过热水进行加热。仿体使用Agar粉和纯净水按照1:100的比例进行配置，管道和加热囊的材质均为聚氯乙烯(PVC)塑料。完成安装后，在加热囊周围插入精度为0.1 ℃的光纤探针测温，探针插入位置如图 6c所示。

插入光纤探针后，把整个仿体放入磁共振扫描仪(Philips Achieva 3.0T TX)进行数据采集。数据采集过程中，首先使用Spin Echo序列对仿体进行结构扫描，扫描主要参数为：FOV=150 mm×150 mm×80 mm，Recon Voxel Size=0.446 4 mm，TR=600~700 ms，TE=10 ms，Flip Angle=70°。图 6b为Spin Echo序列扫描得到的测温层结构图；结构图采集完成后，在加热囊中通过60 ℃的恒温热水，同时打开光纤测温仪器，并使用磁共振快速场回波(fast field echo, FFE)测温序列进行动态数据采集，磁共振参数设置为：FOV=100 mm×100 mm×12 mm，Slices=3，Slice Thickness=4 mm，TR=50 ms，TE=10 ms，Flip Angle =30°，Scans=300，Time Resolution=6.2 s。

临床仿真实验和仿体实验的数据处理使用的运算环境为Anaconda 3.5+Pytorch 1.0.0，硬件平台为Dell Precision T7910工作站，显卡为Nvidia GTX 1080Ti。使用本文所提的BHT-Kalman方法时，单张图片的计算时间为4.5 s。

3 实验结果 3.1 临床数据仿真实验结果

图 7是临床数据仿真实验的结果图。观察图 7可知，相较于传统的SBS方法，本文方法(BHT-Kalman)在引入模拟温度场、使用Kalman滤波方法对图像进行滤波后，重建的温度图(图 7c)整体信噪比更高。将两种方法的重建结果均减去金标准温度图，得到的误差图分别如图 7d和7e所示。

结果表明，在感兴趣区域(field of interest, FOI)内，本文方法重建误差在2 ℃以内，而SBS方法重建图像信噪比较低，FOI内部分低信号强度区域重建温度误差甚至会超过10 ℃。

图 8是2种方法重建温度图与金标准图之间在不同时间帧的均方根误差(root mean square error，RMSE)值的曲线图。观察图 8可知，传统SBS方法重建结果与金标准之间存在6~8 ℃的RMSE值；而本文所提出的方法重建结果与金标准之间的RMSE值大幅下降，并且在不同时间帧均较为稳定，在2 ℃以下。

3.2 仿体实验结果

图 9是仿体实验结果图。对比图 9a和9b可以发现，使用本文的BHT-Kalman方法重建的温度图更平滑。从图 9c和9d升温区域放大图中可以发现，SBS方法重建结果中测温噪点较多，等温线边缘粗糙；而使用BHT-Kalman方法重建结果由于生物传热模型的引入，较好地消除了测温噪点，并且等温线边缘较为均匀，更符合真实情况。

图 10是光纤测温点温度随时间变化图。红色曲线是测温光纤测量得到的温度，可以视为真实温度；黑点为使用SBS方法重建温度；蓝色曲线为使用BHT-Kalman方法重建温度。从图 10可以发现，SBS方法重建结果虽然能够测量出温度变化趋势，但是测温偏差相对较大；而采用BHT-Kalman方法进行重建后，测温精度有了显著提升。计算300幅温度图的测温点的重建温度与真实温度的RMSE值，可得SBS方法和BHT-Kalman方法的RMSE值分别为1.927 ℃和0.735 ℃。

4 结论

本文基于Kalman滤波理论，提出了使用生物传热模型模拟温度场来对磁共振温度成像进行滤波的方法，并设计了临床数据仿真实验和仿体实验对所提方法的性能进行了验证。

通过临床数据仿真实验和仿体实验可以发现，本文提出的BHT-Kalman方法具有抗噪性能好、精度较高的优点，尤其是在磁共振温度图像的信号强度较低的情况下，传统的SBS方法重建结果会存在最高10 ℃和平均7 ℃的均方根误差，而BHT-Kalman方法能够将均方根误差控制在2 ℃以内。然而，BHT-Kalman方法也具有一定的局限性，其精度在一定程度上依赖于对组织热力学参数的选择和对温度变化过程的精确模拟。例如，该方法用于真实的人体组织，需要根据对流换热情况对大血管组织进行处理，对人体中骨骼等疏松孔状结构的热力学建模也需要进一步研究。