无需心电门控与呼吸控制的心脏磁共振实时成像

郑翊宸，陈硕，孙爱琦，赵波，李睿*

1.清华大学医学院生物医学工程系生物医学影像研究中心，北京 100084；2.上海东软医疗科技有限公司，上海 200241；3.麻省总医院，马萨诸塞州波士顿 02114；4.哈佛大学医学院，马萨诸塞州波士顿 02115；*通讯作者 李睿 leerui@tsinghua.edu.cn

心脏病是一种严重危害人类健康的常见疾病，具有高患病率、高病死率、高治疗费用等特点[1]。心脏磁共振电影成像（cine cardiovascular magnetic resonance，CMR）已初步成为心功能评估的“金标准”，并在心脏疾病的精准诊断中发挥重要作用[2]。然而由于MRI扫描速度较慢，且受心脏自身节律性运动及呼吸运动的影响，其位置与形态在成像期间不断变化，CMR检查时往往需要患者屏气配合并使用心电门控，通过采集多个心动周期的数据，最终得到心脏在一个心动周期内随时间变化的运动情况。然而，该方法对于心律不齐患者存在扫描效率较低，甚至无法正常扫描等问题；同时该方法依赖于患者屏气状况，故不完全屏气会造成严重的运动伪影；此外，电影成像结果为经过平均的心动周期，无法反映不同心动周期间，尤其是心率异常周期内的心脏形态及功能。

为克服传统电影方法的问题，研究者们提出了一系列新的心脏成像方法，如并行成像、径向采集、压缩感知等，以实现心脏磁共振快速成像[3]。然而，这些方法存在时间分辨率不足、适用范围受限或依旧需要门控配合等问题，难以获得多心动周期、高时空分辨率的心脏图像，也难以用于心律不齐患者。

本研究拟通过基于低秩模型[4-6]和稀疏约束[7]联合并行成像的方法实现实时心脏成像，实现在不使用心电门控及屏气配合条件下的心脏磁共振采集，从而获得高时间分辨率多个心动周期内心脏动态图像，并且能够捕捉异常心动周期的心脏搏动情况。

1 资料与方法

1.1 理论分析 为实现心脏实时成像，本研究使用低秩模型与稀疏约束对心脏信号进行重建恢复。通常，所采集到关于心脏的磁共振信号可表达为公式（1）。

其中，di(k,t)代表第i个接收线圈在t时刻采集的k空间数据（i的取值范围是1-Nc，Nc代表接收线圈的个数），ρ(r,t)表示t时刻r位置的心脏图像数据，Si表示第i个接收线圈的敏感度图。ηi(k,t)表示第i个接收线圈在t时刻采集的k空间的测量噪声。图像重建是从采集获得的数据di(k,t)求解心脏动态图像ρ(r,t)的逆问题。

考虑离散的图像模型，每个图像序列可通过基于时间空间分布的Casorati矩阵表示，如公式（2）。

其中，M代表时间帧的数目，N代表一个时间帧内对应获取的空间域中像素个数。

由于心脏动态图像中时间与空间高度相关，故Casorati矩阵C具有低秩特性，即公式（3）。

在求解此低秩问题时，对矩阵C进行分解有：

其中，

在这种低秩表示中，U的列和V的行分别代表C的空间子空间和时间子空间。

根据理论分析，得到动态图像需要空间子空间与 时间子空间两部分，在数据采集中数据也分为训练数据与成像数据两种。此处使用了特殊的交错采集轨迹，数据高度降采（图1）。

图1 实时成像方法的采集轨迹。紫色“●”代表训练数据，位于k空间中心低频区域，其时间分辨率较高，满足奈奎斯特抽样定律；红色“○”代表成像数据，位于k空间外周高频区域，通过均匀随机采样获得，空间分辨率较高；虚线将不同时间帧进行分隔，根据每条训练数据与成像数据组成的数据对即可对应重建出一帧图像

将公式（1）写成矩阵向量符号表示的形式，见公式（6）。

其中，di表示从第i个接收线圈中采集到的测量数据，Ω是稀疏采样算子，Fs是空间傅里叶变换矩阵；Si是包括第i个线圈敏感度图的对角矩阵；C是表示心脏动态图像的Casorati矩阵；ni是第i个线圈的测量噪声。

使用联合低秩约束与稀疏约束，结合并行成像方法，图像重建问题可以表示为公式（7）。

其中Rs(UV)表示稀疏约束。

根据设计的采集轨迹，使用仅包含k空间中心但时间分辨率较高的训练数据，通过主成分分析估计时间子空间V，然后根据高空间分辨率的成像数据与时间子空间一起估计空间子空间U。

在求解得到U后，根据C=UV即可求出动态图像序列的Casorati矩阵。

1.2 研究对象 前瞻性纳入2019年8—12月6名无任何心血管疾病的健康志愿者，其中男5例，女1例；年龄21～26岁，平均（23.5±1.6）岁；同时纳入1名26岁女性心律不齐患者进行心脏MRI扫描。其中2名健康志愿者用于对成像方法的优化；4名健康志愿者与1名心律不齐患者用于对实时成像方法进行验证。本研究获得清华大学伦理审查委员会批准，所有受试者均签署知情同意书。

1.3 仪器与方法 使用Philips Ingenia 3.0T MR扫描仪，配合32通道心脏线圈进行扫描。扫描序列包括常规心脏电影方法及本研究所使用的实时扫描方法，对心脏短轴位进行成像，扫描序列：快速场回波（fast field echo，FFE），视野300 mm×208 mm，空间分辨率2 mm×2.4 mm，层厚8 mm，层数7，TR 5.4 ms，TE 3.2 ms，翻转角15°。实时成像方法无需心电门控，在自由呼吸状态下采集数据，每层采集时间51.3 s，时间分辨率10.8 ms。传统电影方法使用心电触发，配合患者屏气，采集26个心脏相位，每层采集时间约18.5 s，时间分辨率约30 ms。在实际扫描过程中，采集时间与分辨率随患者心律存在细微差异。

1.4 研究方法

1.4.1 仿真实验 由于实时成像方法结果受模型参数影响，故需要选择合适的参数。本研究招募2名健康志愿者的6层短轴位数据进行重建参数优化。利用低秩模型及稀疏约束，使用不同的低秩模型阶数L与稀疏约束系数λ，在具有24 GB随机存取存储器、3.40 GHz CPU的工作站上进行重建。

使用MATLAB R2018a计算不同重建结果的信噪比（SNR），进行图像质量评估。使用Image J 1.6.0_24软件分别针对传统电影方法与不同参数下实时成像仿真重建结果在收缩末期与舒张末期的左心室内膜进行勾画，获得左心室面积大小后计算实时方法与电影方法的面积比。在保证左心室面积测量一致的前提下，优化SNR，确定最优参数。

1.4.2 健康志愿者试验 优化成像序列后对健康志愿者进行扫描。由1名具有3年心脏磁共振研究经验的研究者对实时方法与电影方法的图像质量进行随机盲法判读，综合考虑图像的伪影以及结构的可 见度。图像评分包括4个等级：4=优秀，3=良好，2=差，1=无可见度。使用MATLAB R2018a计算实时成像方法与传统电影方法在收缩末期及舒张末期每一层的左心室心内膜与心外膜的面积，进而计算射血分数（ejection fraction，EF）、心输出量（cardiac output，CO）、每搏输出量指数（stroke volume index，SVI）及左心室质量（left ventricle mass，LVM）。由于实时方法可以捕捉到多个心动周期，因此对随机选择的3个心动周期进行计算，将其结果均值与电影方法进行对比，验证两者在心功能评估方面的一致性。

1.4.3 心律不齐患者实验 通过对纳入的心律不齐患者进行实时心脏MRI，验证该方法对捕捉心动异常周期发生时心脏逐节拍变化的可行性，并使用MATLAB R2018a计算心律不齐发生时左心室面积的变化情况，与正常心动周期进行对比。

1.5 统计学方法 使用MedCalc 18.2.1软件，健康志愿者收缩末期及舒张末期的左心室面积、射血分数、心输出量、每搏输出量指数和左心室质量等计量资料以表示。使用Bland-Altman分析及Pearson相关分析实时成像方法和传统电影方法结果的一致性。r>0.9表示一致性良好。P<0.05表示差异有统计学意义。

2 结果

2.1 仿真实验 不同重建参数条件下SNR、左心室面积比的变化见图2。随着低秩模型阶数L的增大，SNR快速下降后逐渐平缓；左心室面积比在舒张末期保持基本一致，在收缩末期下降后保持不变，随着稀疏约束系数λ的增大，SNR快速增大后缓慢上升；左心室面积比在舒张末期保持基本一致，在收缩末期缓慢增加后急速上升。通过比较上述参数的变化，选择L=25和λ=0.5时左心室面积比接近1，同时拥有较高的SNR。

图2 SNR、左心室面积比随模型参数的变化。A.SNR随模型阶数L的变化；B.SNR随稀疏约束系数λ的变化；C.实时方法与电影方法左心室面积比随模型阶数L的变化；D.实时方法与电影方法左心室面积比随稀疏约束系数λ的变化。C与D各点标注了±5%的误差线

2.2 健康志愿者实验 实时成像方法与传统电影成像方法在收缩末期与舒张末期的典型图像见图3。2种方法具有相似的形态。实时成像方法与电影成像方法的图像质量评分分别为（3.75±0.43）分、（3.50±0.50）分，两者图像质量接近，差异无统计学意义（P>0.05）。

实时成像方法与传统电影成像方法每一层的左心室心内膜与心外膜面积在收缩末期与舒张末期的Bland-Altman结果见图4。两者心内膜面积在收缩末期的均数差为-0.1 cm2（95%CI-2.2～2.0 cm2），在舒张末期的均数差为0.2 cm2（95%CI-3.4～3.8 cm2）；心外膜面积在收缩末期的均数差为-0.5 cm2（95%CI-5.4～4.4 cm2），在舒张末期的均数差为-0.1 cm2（95%CI-5.2～4.9 cm2）。两种方法的均数差接近于0，偏倚较小，故两者在左心室面积测定方面一致性良好。

图3 电影方法与实时方法在收缩末期与舒张末期的示例图像。A～D分别为电影方法收缩末期、电影方法舒张末期、实时方法收缩末期、实时方法舒张末期的图像

图4 实时方法与电影方法左心室面积在收缩末期与舒张末期测量心外膜面积一致性的Bland-Altman分析结果

针对分别由实时成像方法与传统电影成像方法获得的心功能参数EF、CO、SVI和LVM的Bland-Altman一致性分析结果见图5。两者EF的均数差为0%（95%CI-8%～8%），CO的均数差为160.3 ml（95%CI-653.1～973.7ml），SVI的均数差为0.9 ml/m2（95%CI-4.9～6.7 ml/m2），LVM的均数差为-2.9 g（95%CI-14.2～8.4g）。4种参数的所有点均位于95%一致性界限内。由实时成像方法与传统电影成像方法获得的心功能参数统计信息及相关系数见表1。4种心功能参数的相关系数均>0.99。由此判断，实时成像方法计算所得心功能参数与传统电影成像方法所得结果具有极高的相关性与一致性。

图5 实时方法与电影方法心功能参数EF（A）、CO（B）、SVI（C）与LVM（D）的Bland-Altman分析结果

表1 实时方法与电影方法心功能参数的相关性分析（±s）

表1 实时方法与电影方法心功能参数的相关性分析（±s）

注：EF为射血分数，CO为心输出量，SVI为每搏输出量指数，LVM为左心室质量

项目 EF（%） CO（ml） SVI（ml/m2） LVM（g） 实时方法 63.10±5.56 5212±1908 35.98±10.81 104.30±32.56 102.10±28.35 r值 0.9934 0.9997 0.9936 0.9916 电影方法 62.56±9.53 5312±2320 35.44±13.45 P值 ＜0.05 ＜0.05 ＜0.05 ＜0.05

2.3 心律不齐患者试验 心律不齐患者不同心动周期间存在差异。心跳时长在正常心动周期后有2次较短的心动周期，紧接1个较长的心动周期，然后恢复正常心跳。发生心律不齐时，舒张末期面积小于正常心动周期，心脏未完全舒张即开始下一次心跳（图6）。

图6 心律不齐患者结果。A.穿过左心室的线在17 s内的变化，箭示发生心律不齐的心动周期；B.发生心律不齐时左心室面积的变化

3 讨论

本研究基于低秩模型与稀疏约束的方法，实现在无需心电门控与屏气的条件下进行心脏实时成像，能够在10.8 ms内（2倍TR）采集一帧CMR图像，时 间分辨率优于传统电影成像方法（约30 ms）。同时，本研究初步验证了实时成像方法对心律不齐患者异常心动周期的结构及功能观测的可行性。

由于实时心脏成像重建结果与低秩和稀疏重建参数有关，本研究首先应用不同重建参数进行实时成像重建，通过优化SNR和左心室面积得到最优实时成像重建参数。低秩约束能够提高图像SNR，但约束过于严格时会抑制心脏运动信息：当阶数参数L较小时，图像的SNR较高，但由于较小的L不能代表心脏的搏动情况，心脏运动受到影响，因而在收缩末期左心室面积被高估；L较大时，可正确反映心脏运动，左心室面积与传统电影方法一致。然而，较大的L会导致迭代时的病态问题，因而SNR较低，图像质量受到影响。稀疏模型中，λ较小时，约束被抑制，无法得到相对平滑的结果，因而SNR相对较低。λ较大时，过强的约束会导致结果的过度平滑，抑制了心脏的运动，收缩期与舒张期的差距变小，因而在收缩末期时左心室的面积被高估。本研究中，低秩约束阶数L=25和稀疏约束系数λ=0.5保证了SNR和心脏运动（左心室面积）的平衡。

本研究通过分析健康志愿者心功能参数以验证实时成像方法与传统电影成像方法的一致性。左心室面积的Bland-Altman一致性分析中，2种方法的偏倚较小，充分证明实时方法左心室心内膜与心外膜面积与电影方法高度一致。EF、CO、SVI和LVM心功能参数的Bland-Altman一致性分析中，所有点均位于95%一致性界限内，实时方法与电影方法差值的均值接近0，归一化后的百分比分别是0.2%、1.1%、1.1%和-2%，两者均值差百分比均≤2%，较小的差距进一步体现其一致性。结合r>0.99，表明实时方法在心功能参数上可代替传统的电影成像方法。

本研究通过对心律不齐患者的扫描，验证了实时方法捕捉心脏逐节拍变化，特别是心律不齐发生时的可行性。本研究所纳入的心律不齐患者心律不齐发生频率较低。传统的心脏电影成像方法无法捕捉并评估心律不齐发生时心脏结构和功能的变化。本研究所提出的实时心脏成像能够捕捉并反映心律不齐发生时患者左心室的结构和面积变化，且发现在发生心律不齐时患者心脏先经历2次较短的心动周期，紧接着经历1个较长的心动周期，这对解释和观察心律不齐患者疾病的发生及发展规律具有潜在的辅助作用。

多种快速成像方法也应用于心脏成像。衡阳等[8]应用图像序列差分正则项实现了CMR图像重构；陈思吉等[9]基于稀疏和低秩先验分离的方法进行快速动态CMR重建；但这2种方法的欠采样率有限，难以满足高时间分辨率的要求；另外均为回顾性仿真，未进行前瞻性验证，评价指标仅反映图像质量，而无心功能相关研究。陈辉等[10]应用优化实时并行采集技术的成像序列对心律失常患者进行评估，但图像的时间分辨率超过60 ms，难以观察到发生心律不齐时心脏的细微运动情况。本研究通过联合低秩约束与稀疏约束，结合并行成像方法实现了无需心电门控和屏气配合的实时电影成像，解决了临床心脏成像前准备复杂、患者配合欠佳、心律不齐无法评估等多种问题。

本研究从仿真、健康志愿者实验、心律不齐患者实验方面反映实时心脏成像方法的可行性与准确性，但也存在以下局限性：①纳入人数相对较少，未来需要纳入更多的志愿者进一步佐证方法的稳定性和可重复性；同时需要纳入更多不同病情心律不齐的患者，验证该方法对不同种类的心脏异常运动捕捉的可行性，并通过进一步计算发生心律不齐时心动周期的定量参数，为临床提供更直观的依据。②计算的临床参数有限，未来可以计算更多种类的临床参数（如与右心室相关参数等）；另外采集方位也不仅限于短轴位，为临床提供更多信息。③采集时间相对较长，由于无需门控配合，使用较多的数据量以迭代恢复动态图像，未来需要优化数据量与图像质量的关系。④重建速度相对较慢，由于数据量较大，重建单层50 s的连续图像需要10～20 min，未来通过优化重建方法及使用高性能计算机能够提升重建速度，便于临床使用。

总之，基于低秩模型与稀疏约束的心脏实时成像方法能够在无需心电门控与呼吸控制的情况下获得高时间、空间分辨率的心脏动态图像，且与传统电影方法具有良好的一致性，并能够捕捉到心脏的异常搏动。实时方法具有扫描容易、适用性广等特点，值得进一步研究与应用。