基于自监督深度学习的心肌磁共振T1量化成像运动校正算法研究

吴春燕，李雨泽，丁海艳，陈慧军

清华大学 医学院 生物医学工程系，北京 100084

引言

心脏磁共振T1量化成像可以指示心肌组织的一系列病理改变，特别是心脏的炎症、急（慢）性心肌瘢痕以及纤维化引起的细胞外间质重组[1-2]。临床上，改良Look-Locker 反转恢复（Modified Look-Locker Inversion Recovery，MOLLI）序列是标准的心肌T1参数量化序列，但是扫描过程中心脏及呼吸运动将导致图像质量降低，需对采集的不同帧图像进行配准以提高图像质量及T1量化精度。

国际上已有许多研究者对心肌的T1量化提出了运动校正方法，Roujol等[3]提出一个基于光流法的局部非刚性配准架构，同时估计运动场和信号强度变化；El-Rewaidy等[4]利用一种活动形状模型分割得到的心肌轮廓实现配准，提高了各心肌节段T1参数量化的准确度和精度；Tilborghs等[5]提出一种结合数据驱动初始化和基于模型的精细改良的配准方法，实现了更加鲁棒的T1参数量化。但是，这些传统配准方法需要多次迭代优化，耗时长，降低效率[6-7]。Fahmy等[8]提出一种基于深度学习的自动心肌T1参数量化方法，然而，这一深度学习算法为监督式的，需要人工标注数据，且在自动分割效果不佳时需要额外的后处理操作。

为了克服前述算法的缺点，本研究提出一种基于深度学习的心肌T1参数量化成像运动校正算法，无需传统方法的多次迭代优化，可大幅度降低运算时间；同时，采用自监督的训练方法，无需人工标注的金标准图像，大大降低了医生的工作量，为大规模临床应用提供了可能。

1 材料与方法

1.1 数据集

实验数据来自47例健康志愿者（30例男性，17例女性，年龄30～64岁），纳入标准：① 既往无心源性或系统性疾病；② 既往CMR扫描无心脏形态和功能异常；③ 无磁共振检查禁忌证。所有受试者均签署了知情同意书。

实验使用3.0 T磁共振扫描仪（Achieva TX，Philips，荷兰），32通道的心脏线圈，所有受试者在屏气情况下进行成像。采用MOLLI序列进行扫描，每名志愿者采集1～3层图像，扫描方位为左心室短轴位，每层图像包含8帧T1加权图像。序列的具体扫描参数如下：重复时间2.59 ms；回波时间1.31 ms；触发延迟（Trigger Delay，TD）时间为舒张末期；完成一次扫描的时间为8～34 s；翻转角35°；扫描视野30 cm×15 cm；分辨率1.25 mm×1.25 mm；层厚8 mm。

1.2 基于深度学习的心肌T1参数量化成像运动校正算法

本方法的总体流程主要包括两个部分：① 利用基于深度学习的自监督配准网络将MOLLI图像中除固定图像以外的所有图像配准到固定图像上；② 对配准后的T1加权图像进行参数拟合得到T1量化图像。

1.2.1 基于自监督深度学习的非刚性配准算法

配准流程图如图1所示。图像配准通常是将一幅移动（或待配准）图像配准到固定（或参考）图像上。令m、f分别表示定义在2D空间域Ω∩R2上的移动、固定图像，利用卷积神经网络（Convolutional Neural Network，CNN）对配准运动场u建模，即gθ(f,m)=u，其中g为神经网络，θ是可学习参数。网络可学习每一组输入对(m,f)之间的配准变换；接着，运动场利用空间变换函数(º)作用到移动图像m上，可得到运动校正后的图像m'。具体地，对于成像空间中每个像素p，经运动场u的空间交换(°)得到的m´(p)是其周围8个像素点的线性插值，见式(1)。

图1 基于自监督深度学习的非刚性配准算法

其中，Ζ8(p')表示p'的8像素邻域，p'=p+u(p)。

本方法中的自监督体现在，由于固定图像来源于直接采集到的图像，无需人为标注或额外的处理，经过运动场作用后的图像m'与固定图像进行损失函数计算并反向传播优化网络，达到自监督学习的目的。

1.2.2 神经网络结构

图2为本方法使用的神经网络的结构。该网络采用Unet结构[9]，主要由编码器和解码器两个部分组成：编码器利用卷积层增加CNN的感受野，降低CNN的空间尺寸，同时学习不同尺度下的图像信息；解码器采用上采样恢复图像尺寸，并利用跳跃连接融合编码器和解码器部分的多尺度信息。编码和解码阶段都采用尺寸为3×3、步长为2的卷积和Leaky ReLU激活函数（小于0部分的斜率为0.2）。

图2 运动场估计的网络结构图

1.2.3 损失函数

本方法的损失函数由两部分构成，一部分为自监督损失函数，即图1中的sim，sim衡量的是校正后图像与固定图像的相似度。考虑到T1加权图像不同帧之间对比度相差大，直接采用基于信号强度的相似性测度效果不佳[10]，im选择对于对比度变化不敏感的互相关（Cross-Correlation，CC）和互信息（Mutual Information，MI）作为相似性测度，见式(2)～(3)。

P(a)是f中信号强度等于a的像素比例，P(b)是m'中信号强度等于b的像素比例，P(a,b)是f和m'信号强度的联合分布。为了使损失函数可微分以便于神经网络的反向传播，使用Parzen窗法近似公式中的概率分布[11]，见式(4)～(6)。

其中，窗函数W选择高斯窗。

另外，为了对比采用基于信号强度的相似性测度的运动校正效果，实验也验证了sim取均方误差（Mean Square Error，MSE）的情况，见式(7)。

综上，网络的总损失函数如式(9)所示。

本研究测试不同的相似性测度（MSE、CC、MI）以寻找最优的损失函数。

1.2.4 训练与测试

网络使用TensorFlow 2.1.0框架实现，优化器为Adam，学习率设置为1e-4，beta1=0.9，beta2=0.999；采用的计算服务器配置为I7-6950K CPU, NVIDIA TITAN Xp 12 GB GPU及32 GB内存。在本研究中，固定图像为MOLLI序列采集到的T1加权图像中的第8帧，因为此时心肌-血池对比度较高，移动图像则为其他帧的图像。47例受试者数据中，36例（2035张图像）用于网络训练，11例（650张图像）用于测试。

此外，本文对比了一种传统配准方法：基于互相关的B样条自由形变模型配准算法[12]，由Matlab医学图像配准工具箱（Medical Image Registration Toolbox，MIRT）实现。

1.3 T1参数拟合

图像经过配准后，还需经过T1拟合才可得到T1量化图像。如图3所示，本实验采用的MOLLI序列由两次ECG触发的Look-Locker（LL）实验（LL1、LL2）组成，分别由三、五次单次激发读出。每个LL实验中，反转脉冲和TD之后，在递增的反转恢复时间（Inversion Time，TI）心脏舒张末期采集信号，LL1和LL2之间至少间隔4 s保证完全的信号恢复。

图3 MOLLI序列信号采集示意图

本文采用三参数模型进行T1拟合[13]，见式(10)。

对于MOLLI序列，拟合得到的T1*经过如下的修正后得到T1：T1=T1*(-B/A-1)。本文采用叶猛[14]提出的结合RD（Reduced Dimension Nonlinear Least Square with Rolarity Restorations）算法和LM（Levenberg-Marquardt）算法的拟合流程，即对运动校正后的T1加权图像加mask去除背景后，先进行RD拟合，再对拟合残差过大的点用LM拟合。

1.4 评价指标及统计学分析

本文评价了T1加权图像的配准效果及T1量化误差。采用Dice相似系数（Dice Similarity Coefficient，DSC）和平均边界误差（Mean Boundary Error，MBE）作为量化评估配准误差的指标。这两个指标的获得需要在原始加权图像上手动勾画出心肌内、外轮廓。指标具体计算可参考Tilborghs等[5]的方法。采用T1量化图像的标准差（Standard Deviation，SD）[15]评估T1参数量化的准确度，可利用T1拟合残差和测量噪声的SD近似得到待估计参数的协方差矩阵，从而得到T1参数量化的SD。除了SD 量化图，本文还计算了心肌部分的平均SD。

本文使用Matlab R2020a软件进行统计学分析。DSC、MBE和心肌SD表示为±s；对不同方法测量得到的DSC、MBE和心肌SD进行Wilcoxon符号秩检验，P<0.05表示具有统计学意义。

2 结果

2.1 T1加权图像配准结果

经过本文提出的深度学习算法（以CC和MI作为相似性测度）配准后，T1加权图像的心肌部分得到了较好的运动校正。T1加权图像配准效果如图4所示。未配准时（图4a、 e、 i中蓝色箭头所示），心肌位置相对于固定图像发生了明显的运动偏移，采用深度学习算法配准后，相似性测度选择CC和MI时（图4c、g、k和图4d、h、l），心肌位置与固定图像的基本一致。当心肌和血池对比度很低时（图4i、k、l），算法也能较好地配准；相似性测度选择MSE时，配准后T1加权图像出现明显的变形、失真（图4b、f、j中红色箭头所示）；当心肌和血池对比度很低时（图4j中红色箭头所示），失真情况更加严重。

图4 三例未配准和深度学习算法（MSE、CC、 MI）配准后的T1加权图像

定量评估算法的配准准确度结果如表1所示。深度学习（CC）配准后DSC较原始未配准显著提高（P<0.05），而采用传统算法和深度学习（MI）配准后DSC指标没有显著改善（P=0.13和P=0.98）。传统算法、深度学习算法（CC、MI）配准后MBE均显著降低（P<0.05），且深度学习算法（CC、MI）的MBE显著低于传统算法（P<0.05）。而深度学习（MSE）配准后，DSC和MBE均比未配准时差。综合DSC和MBE两个指标，本文提出的深度学习算法（CC）配准准确度最高。

表1 未配准和传统算法、深度学习算法（MSE、CC、MI）配准后的DSC和MBE

2.2 T1参数量化结果

图5展示了不同配准算法运动校正结果的实例。从图5中放大部分可以观察到，原始T1参数量化图心肌内膜附近（图5中箭头所示）存在运动伪影，降低了T1量化的准确性。传统算法和深度学习算法（CC、MI）在不同程度上减少了运动伪影，其中深度学习（CC）的运动校正效果最佳。观察对应的SD量化图，除深度学习（MSE）之外，应用不同算法配准后，心肌部分的SD在整体上均减小了，深度学习（CC）的心肌SD整体上最小。

图5 未配准和传统算法、深度学习算法（MSE、CC、MI）配准后的T1参数量化图和对应的SD量化图

心肌部分的T1参数量化误差如图6所示。传统算法、深度学习（CC）、深度学习（MI）配准后心肌部分的平均SD分别为（68.74±31.27）、（59.22±29.26）、（60.35±28.47）ms，均小于原始图像心肌部分的平均SD（74.37±33.29） ms，差异有统计学意义（P<0.05），且深度学习（CC）和深度学习（MI）的平均SD明显小于传统算法（P<0.05），深度学习（CC）和深度学习（MI）的平均SD无统计学差异（P=0.056）。深度学习（MSE）配准后心肌部分的平均SD增加。

3 讨论

本研究提出并验证了一种基于自监督深度学习的心肌T1参数量化成像运动校正算法。该方法突出优点为基于深度学习的配准算法，提高效率及准确度；另一个优势为采用了自监督学习框架，无需人工标注数据，大大节省了医生的工作量，为心脏T1量化在临床上大面积使用提供可能。

图6 未配准图像和传统算法、深度学习算法（MSE、CC、MI）配准后心肌T1参数量化的平均SD的均值和标准差

首先，深度学习是本方法的优势之一。深度学习适合处理大数据问题，在图像识别、分类、分割、配准等计算机视觉任务上取得了许多重大研究成果。区别于传统配准算法依赖于迭代优化，将深度学习应用于T1加权图像配准，可达到一步配准，大幅提高计算效率。经过测算，本研究提出的方法完成一例配准平均用时0.4 s，比传统算法快20倍。在配准精度及T1量化准确度方面，本研究在MOLLI序列采集到的47例受试者的左心室短轴切片数据上验证了提出的运动校正算法，与传统的基于互相关的B样条自由形变模型配准算法相比，深度学习（CC）算法的MBE降低了12%，心肌T1量化SD降低了14%，心肌T1加权图像的配准准确度和T1参数量化准确度都有所提高。另外一个值得讨论的地方是，本研究测试了不同的损失函数，以寻求最适宜本问题的方案。本研究测试的三种损失函数中，CC和MI测度可以有效校正心肌T1参数量化图中的运动伪影，CC测度的方法也得到了最佳的效果；然而，选择MSE作为相似性测度时，本文提出的算法无法提供有效的运动校正，可能的原因为误差无法解决信号强度和对比度变化的图像配准问题，而CC和MI在衡量图像相似程度时摒除了图像对比度变化带来的影响。

其次，采用自监督学习框架是本方法的另一个优势。深度学习算法的一个挑战是训练真值的获得，与Fahmy等[8]提出的方法需要大量预先勾画的心肌轮廓作为训练集不同，本文提出的方法是自监督的，网络训练过程无需真值，实际应用中也不需要预处理和用户交互。该方法提出的自监督配准网络引入了空间变换层，使得在训练阶段不需要运动变换真值就可进行损失函数的计算和网络优化，大大降低了医生的工作量，可以将本方法自由应用在任何相似的图像处理应用当中，如造影剂增强的磁共振心脏、肝脏成像等需要运动配准的成像方式中。

在后续的研究中，为了提升算法在运动偏移幅度过大情况下的鲁棒性和准确度，将考虑以下改进策略：① 将心肌结构相似性损失加入到目前的损失函数中[16]；② 在目前非刚性配准网络上增加一个仿射配准子网络[17]以克服运动幅度过大的问题。

4 结论

本文提出了一种基于自监督深度学习的非刚性配准算法，为心肌T1参数量化成像提供了快速、有效的运动校正，增加了T1参数量化成像大规模临床应用的可能性。