1. 清华大学 工程物理系, 粒子技术与辐射成像教育部重点实验室, 北京 100084;
2. 清华大学 医学院, 生物医学工程系, 北京 100084;
3. 航天中心医院 骨科, 北京 100048
收稿日期：2019-04-10
基金项目：国家自然科学基金项目（61571257）
作者简介：吴锦超(1994—), 男, 硕士研究生
通讯作者：应葵, 副教授, E-mail:yingkui@tsinghua.edu.cn
摘要：磁共振温度成像能够为肿瘤热疗消融提供实时、全局的温度场监控，是保障消融安全、有效进行的重要技术手段。然而，临床中磁共振温度成像信噪比较低，该现象在选用快速成像序列时尤为严重。针对这一问题，该文提出了基于Kalman滤波和生物传热模型对磁共振温度成像进行滤波的方法。该方法将生物传热模型改写成Kalman滤波方法的状态转移矩阵形式，并将模拟温度和磁共振测温值相结合，以获得具有高精度和高信噪比的估计温度。临床仿真实验表明，该方法能将测温均方根误差从6℃降至2℃；仿体实验表明，测量值与真实值的均方根误差从1.927℃下降到0.735℃，且信噪比显著提升。
关键词：核磁共振成像磁共振温度成像生物传热模型Kalman滤波
Kalman filtering and bio-heat transfer model based real-time MR temperature imaging for increased accuracy
WU Jinchao1, QIU Shihan2, LI Muheng1, WEI Xing3, CHEN Bingyao3, YING Kui1
1. Key Laboratory of Particle and Radiation Imaging of Ministry of Education, Department of Engineering Physics, Tsinghua University, Beijing 100084, China;
2. Department of Biomedical Engineering, School of Medicine, Tsinghua University, Beijing 100084, China;
3. Department of Orthopedics, Aero Space Center Hospital, Beijing 100048, China
Abstract: Magnetic resonance temperature imaging is an important technique to ensure safe and effective use of tumor hyperthermia ablation by providing real-time, global temperature field monitoring. In clinical trials, however, the signal-to-noise ratio (SNR) of magnetic resonance temperature imaging is relatively low, with the signal quality degrading more with fast imaging sequences. To solve this problem, a bio-heat transfer based Kalman filtering model is developed for magnetic resonance temperature imaging where the bio-heat transfer equation is transformed into the form of a Kalman state transition matrix. Then the simulated temperature is combined with the measured temperature to create an accurate, high SNR estimated temperature. Clinical simulations show that this method reduces the temperature measurement root mean square error from 6℃ to 2℃ and the physical phantom experiment shows that this method reduces the root mean square error of the measured temperature and the true temperature from 1.927℃ to 0.735℃ while significantly improving the SNR.
Key words: magnetic resonance imagingmagnetic resonance temperature imagingbio-heat transfer modelKalman filter
热疗消融具有副作用小、无创或微创等优点，已经成为继手术、放疗、化疗和免疫疗法之后的第5大肿瘤疗法[1]。为保障热疗消融的安全性和有效性，必须对消融过程进行实时准确的温度监测。磁共振温度成像(magnetic resonance temperature imaging，MRTI)由于具有非介入式、时间和空间分辨率较高、能够任意方向成像等优点，是当前热疗消融温度监控的最佳技术手段[2-3]。
磁共振温度成像可以基于多种参数进行温度测量，主要包括：纵向弛豫时间T1、横向弛豫时间T2、水质子共振相移(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=\frac{\Delta \varphi}{\alpha \cdot \gamma \cdot B_{0} \cdot T_{\mathrm{E}}}+T_{0}.$ | (2) |
| $\rho c \frac{\partial T}{\partial t}=\lambda \nabla^{2} T-W_{b} c_{b}\left(T-T_{b}\right)+Q.$ | (3) |
| $\rho c \frac{\partial T}{\partial t}=\lambda^{\prime} \nabla^{2} T+Q.$ | (4) |
| $\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) |
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) |
在使用MRTI方法测温的过程中，C(n)是M×M单位矩阵，因此观测向量Tm(n)为N×1向量，即测量温度矩阵展开得到的向量。
过程方程：
| $\boldsymbol{T}(n+1)=\boldsymbol{F}(n+1, n) \boldsymbol{T}(n)+\boldsymbol{u}(n)+\boldsymbol{v}_{2}(n).$ | (7) |
现确定生物传热模型中的状态转移矩阵F。观察式(5)可知，节点(i, j)在n+1时刻的温度是n时刻该节点(i, j)本身温度以及其相邻4个节点(i-1, j)、(i+1, j)、(i, j-1)、(i, j+1)温度的线性组合加上一热源相关项，因此可通过将二维温度矩阵变换成为温度列向量的方法，将式(5)的差分方程改写成状态方程的形式。二维温度矩阵变换成温度列向量的方法如图 1所示，对应的状态方程矩阵形式为
 |
| 图 1 (网络版彩图)温度矩阵变换成温度列向量示意图 |
| 图选项 |
| $\boldsymbol{T}(n+1)= \boldsymbol{F T}(n)+ \boldsymbol{u}(n).$ | (8) |
为了最终确定状态转移矩阵F各位置的元素值，取矩阵F第a0行进行分析。具体计算步骤为(见图 2)：
 |
| 图 2 (网络版彩图)矩阵计算过程示意图 |
| 图选项 |
步骤1 ??首先通过式(9)确认列向量T中第a0行各元素所在原二维矩阵的行列坐标值(i, j)；
步骤2 ??求解二维矩阵中与节点(i, j)相邻的4个节点的坐标，并根据式(10)计算出4个节点转换到列向量T后的行坐标值，分别记为a1、a2、a3、a4；
步骤3 ??根据矩阵和向量的乘法原理，T(n+1)第a0行的值为F第a0行向量的各元素与T(n)各元素乘积之和，因而可以通过式(5)反推出F矩阵第a0行第a1、a2、a3、a4列的值。
使用列向量T的行坐标值计算原二维矩阵节点坐标值，
| $\begin{array}{c}x=\bmod (a, m)+1 , \\y=\operatorname{ceil}\left(\frac{a}{m}\right).\end{array}$ | (9) |
| $a=(x-1)·m+y.$ | (10) |
| $\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) |
对于系统输入向量un，可通过式(12)求得，
| $\boldsymbol{u}^{n}=\frac{\Delta t \boldsymbol{\phi}}{\boldsymbol{\rho} \cdot \boldsymbol{c}}.$ | (12) |
图 3给出了使用Kalman滤波和生物传热模型模拟温度场对磁共振温度成像进行滤波并计算最佳估测温度图的完整算法。算法的完整推导参考文[12]。
 |
| 图 3 生物传热模型与磁共振测温Kalman滤波算法 |
| 图选项 |
图 3中，系统噪声协方差矩阵Qw和观测噪声协方差矩阵Qv均为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所示。
 |
| 图 4 (网络版彩图)临床数据仿真实验示意图 |
| 图选项 |
之后，如图 4c所示，在前列腺给定区域使用二维Gauss热源进行模拟加热，热源总功率为800 W/m；为了获得作为真实值的对比图像，设置人体组织区域的热力学参数分别为等效导热系数为1.28 W/(m·K)、比热容为4.2×103 J/(kg·℃)、密度为103 kg/m3进行生物传热模型模拟，冷却水恒定为35 ℃，人体远离加热区域的温度恒定为37 ℃。模拟得到的不同时间帧的温度场通过式(1)转换成相位场, 并叠加到基础相位图上，得到仿真相位图。
图 5a是通过生物传热模型在给定热源和组织热力学参数下，模拟出的时间帧1、4、7、10时刻的温度场。图 5b为根据对应时刻温度场生成的相位图; 观察黑色箭头区域可以发现，随着温度上升，对应的相位值也在发生变化。
 |
| 图 5 (网络版彩图)相位模拟示意图 |
| 图选项 |
最后，分别使用基于PRFS的单基准相位相减法(single baseline subtraction, SBS)和本文提出的方法对仿真相位图进行重建得到温度场，并对比分析重建结果。
2.2 仿体实验为了探究真实情况下本文方法的可适用性，使用琼脂仿体设计实验。图 6a是仿体实验设计图，其中立方体为琼脂仿体，在琼脂仿体中插入管道和加热囊(见图 6b)，用来通过热水进行加热。仿体使用Agar粉和纯净水按照1:100的比例进行配置，管道和加热囊的材质均为聚氯乙烯(PVC)塑料。完成安装后，在加热囊周围插入精度为0.1 ℃的光纤探针测温，探针插入位置如图 6c所示。
 |
| 图 6 仿体实验设计图 |
| 图选项 |
插入光纤探针后，把整个仿体放入磁共振扫描仪(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所示。
 |
| 图 7 (网络版彩图)临床数据仿真实验结果对比 |
| 图选项 |
结果表明，在感兴趣区域(field of interest, FOI)内，本文方法重建误差在2 ℃以内，而SBS方法重建图像信噪比较低，FOI内部分低信号强度区域重建温度误差甚至会超过10 ℃。
图 8是2种方法重建温度图与金标准图之间在不同时间帧的均方根误差(root mean square error，RMSE)值的曲线图。观察图 8可知，传统SBS方法重建结果与金标准之间存在6~8 ℃的RMSE值；而本文所提出的方法重建结果与金标准之间的RMSE值大幅下降，并且在不同时间帧均较为稳定，在2 ℃以下。
 |
| 图 8 SBS方法和BHT-Kalman方法重建结果与金标准图之间在不同时间帧的RMSE值 |
| 图选项 |
3.2 仿体实验结果图 9是仿体实验结果图。对比图 9a和9b可以发现，使用本文的BHT-Kalman方法重建的温度图更平滑。从图 9c和9d升温区域放大图中可以发现，SBS方法重建结果中测温噪点较多，等温线边缘粗糙；而使用BHT-Kalman方法重建结果由于生物传热模型的引入，较好地消除了测温噪点，并且等温线边缘较为均匀，更符合真实情况。
 |
| 图 9 (网络版彩图)仿体实验重建结果对比 |
| 图选项 |
图 10是光纤测温点温度随时间变化图。红色曲线是测温光纤测量得到的温度，可以视为真实温度；黑点为使用SBS方法重建温度；蓝色曲线为使用BHT-Kalman方法重建温度。从图 10可以发现，SBS方法重建结果虽然能够测量出温度变化趋势，但是测温偏差相对较大；而采用BHT-Kalman方法进行重建后，测温精度有了显著提升。计算300幅温度图的测温点的重建温度与真实温度的RMSE值，可得SBS方法和BHT-Kalman方法的RMSE值分别为1.927 ℃和0.735 ℃。
 |
| 图 10 (网络版彩图)仿体实验温度曲线 |
| 图选项 |
4 结论本文基于Kalman滤波理论，提出了使用生物传热模型模拟温度场来对磁共振温度成像进行滤波的方法，并设计了临床数据仿真实验和仿体实验对所提方法的性能进行了验证。
通过临床数据仿真实验和仿体实验可以发现，本文提出的BHT-Kalman方法具有抗噪性能好、精度较高的优点，尤其是在磁共振温度图像的信号强度较低的情况下，传统的SBS方法重建结果会存在最高10 ℃和平均7 ℃的均方根误差，而BHT-Kalman方法能够将均方根误差控制在2 ℃以内。然而，BHT-Kalman方法也具有一定的局限性，其精度在一定程度上依赖于对组织热力学参数的选择和对温度变化过程的精确模拟。例如，该方法用于真实的人体组织，需要根据对流换热情况对大血管组织进行处理，对人体中骨骼等疏松孔状结构的热力学建模也需要进一步研究。
参考文献
| [1] | IRVING J, CEPEDA M F J, VALDéS F, et al. Microwave ablation:State-of-the-art review[J]. OncoTargets and Therapy, 2015, 8: 1627-1632. |
| [2] | RIEKE V, PAULY K B. MR thermometry[J]. Journal of Magnetic Resonance Imaging, 2008, 27(2): 376-390. |
| [3] | 刘静, 邓中山. 肿瘤热疗物理学[M]. 北京: 科学出版社, 2008. LIU J, DENG Z S. Physics of tumor hyperthermia[M]. Beijing: Science Press, 2008. (in Chinese) |
| [4] | GRISSOM W A, LUSTIG M, HOLBROOK A B, et al. Reweighted ?1 referenceless PRF shift thermometry[J]. Magnetic Resonance in Medicine, 2010, 64(4): 1068-1077. DOI:10.1002/mrm.22502 |
| [5] | DE SENNEVILLE B D, MOUGENOT C, MOONEN C T W. Real-time adaptive methods for treatment of mobile organs by MRI-controlled high-intensity focused ultrasound[J]. Magnetic Resonance in Medicine, 2007, 57(2): 319-330. DOI:10.1002/mrm.21124 |
| [6] | GRISSOM W A, RIEKE V, HOLBROOK A B, et al. Hybrid referenceless and multibaselinesubtraction MR thermometry for monitoring thermal therapies in moving organs[J]. Medical Physics, 2010, 37(9): 5014-5026. DOI:10.1118/1.3475943 |
| [7] | TODD N, PAYNE A, PARKER D L. Model predictive filtering for improved temporal resolution in MRI temperature imaging[J]. Magnetic Resonance in Medicine, 2010, 63(5): 1269-1279. DOI:10.1002/mrm.22321 |
| [8] | ZHANG Y X, CHEN S, DENG K X, et al. Kalman filtered bio heat transfer model based self-adaptive hybrid magnetic resonance thermometry[J]. IEEE Transactions on Medical Imaging, 2017, 36(1): 194-202. DOI:10.1109/TMI.2016.2601440 |
| [9] | DE SENNEVILLE B D, ROUJOL S, HEY S, et al. Extended Kalman filtering for continuous volumetric MR-temperature imaging[J]. IEEE Transactions on Medical Imaging, 2013, 32(4): 711-718. DOI:10.1109/TMI.2012.2234760 |
| [10] | XUAN Y M, ROETZEL W. Bioheat equation of the human thermal system[J]. Chemical Engineering & Technology, 1997, 20(4): 268-276. |
| [11] | YUAN P. Numerical analysis of an equivalent heat transfer coefficient in a porous model for simulating a biological tissue in a hyperthermia therapy[J]. International Journal of Heat and Mass Transfer, 2009, 52(7-8): 1734-1740. DOI:10.1016/j.ijheatmasstransfer.2008.09.033 |
| [12] | 杨世铭, 陶文铨. 传热学[M]. 4版. 北京: 高等教育出版社, 2006. YANG S M, TAO W Q. Heat transfer[M]. 4th ed. Beijing: Higher Education Press, 2006. (in Chinese) |
| [13] | 张贤达. 现代信号处理[M]. 北京: 清华大学出版社, 1995. ZHANG X D. Modern signal processing[M]. Beijing: Tsinghua University Press, 1995. (in Chinese) |
