基于图像序列差分正则项的实时心脏磁共振重构

衡阳，徐剑峰，陈峰，汤敏*

1.南通大学电子信息学院，江苏南通 226007；2.通科微电子学院，江苏南通 226007；3.南通大学-南通智能信息技术联合研究中心，江苏南通 226007；4.南通大学附属医院医学影像科，江苏南通 226007；5.南通大学电气工程学院，江苏南通 226007；

心脏磁共振（cardiac magnetic resonance，CMR）具有多平面成像、多序列成像、无电离辐射、图像质量好、可重复性强等优点，能够全面提供心脏及其邻近结构的信息，可对心脏形态结构、心室功能、心肌活性等生理方面做出准确评价[1-4]。然而，CMR成像时间有时长达1 h以上，制约了其临床应用。随着场强、梯度硬件、脉冲序列等性能上的改进，成像速度有所提高，但是高速切换的梯度场制造成本昂贵，还会刺激受检者的神经与肌肉，引起涡流扰动。因此，从硬件设备角度提高成像速度已基本达到极限，从算法角度入手是目前解决这一瓶颈问题的突破口[5-6]。

压缩感知（compressed sensing，CS）理论是指在某一变换域内具有稀疏表达的信号，以远低于奈奎斯特采样定理标准的方式采集数据，通过与变换基非相干的随机投影采样，运用合适的优化算法高概率精确重构原始信号[7]。由于 CMR成像属于动态 MRI范畴，时空相关度较高，即序列图像上保持不变的背景部分较多，变化的动态部分较少，相邻像素的相似程度较高。因此，可将目前常见的CMR图像重构算法大致分为串行和并行两大类（图1）。串行重构是指在重构算法运行之前必须搜集所有帧的K空间数据，利用整个数据的冗余性重构图像序列，这类算法通常重构精度较高，但是需要调整的参数较多，因此速度较慢，难以达到实时性的要求；并行重构则是假设相邻两帧的图像之间差异很小，通常在密度域或小波域进行，通过稀疏正则化重构并更新图像，这类算法的特点是计算速度较快，但是由于存在误差累积，易导致重构结果不够精确，因此不适用于相对较长的序列图像的重构。

图1 CMR图像重构算法分类及流程。红色实线表示并行重构，蓝色虚线表示串行重构

鉴于上述两类算法的差异，结合串行和并行算法的各自优点，本研究提出基于图像序列差分正则项的实时CMR图像重构算法，力争在重构结果精确性和实时性之间取得较好的平衡。该算法的基本思路是：以第1帧图像作为参考图像，从K空间数据重构出第1帧图像后，利用图像序列差分正则化使得后续帧的重构均以第1帧图像为参考，从而同时利用时间域和空间域的稀疏性。采用改进的NESTA算法对模型进行求解，提高算法的实时性。以峰值信噪比（peak signal to noise ratio，PSNR）、相对误差（relativel2norm error，RLNE）和均方根误差（root mean square error，RMSE）作为定量评价指标，结合人眼视觉感受以及局部区域放大，对比分析本算法与 kt FOCUSS[8]、基于运动补偿的 kt FOCUSS[9]、kt SLR[10]、DTV[11]等主流算法的性能优劣。

1 算法关键技术

1.1 CS-MRI一般模型 CS重构算法的一般形式见公式（1），亦可写作拉格朗日形式，见公式（2）。

一般而言，CS重构算法主要分为贪婪算法、约束优化算法、迫近算法和同伦算法4大类。其中，贪婪算法采用l0范数最小化实现，因其非多项式难解问题，一般只能寻找局部最优解，算法简单快捷，但对于具体应用条件有严格规定。约束优化算法采用l1范数最小化来替代原始的无约束问题，尽管重构效果良好，但是计算量剧增。迫近算法将原始问题转换为相应的迫近算子，计算复杂度低，易于实现，被视为解决非光滑、有约束的大规模数据优化问题的有力工具，但重构精度有待提高。同伦算法通过调整同伦参数，用于l1范数最小化不断迭代寻找最优解。对于算法特点和数据量而言，压缩感知重构算法可以分成一阶和二阶算法两大类，其中一阶算法最多只求一阶导数，计算速度较快，但是不够精确；二阶算法在求解时较为准确，但是计算复杂度很高，难以适应和满足大规模计算的应用场合。

1.2 图像序列差分正则项 TV模型利用相邻像素值差别小的特征对图像进行稀疏表示，在去除噪声和不需要的小尺度细节信息外，保留图像不连续的边缘信息，提高二维图像重构精度，一般模型见公式（3）[12-13]。

其中，下标i和j表示像素x的空间位置。

尽管传统TV在动态MRI重构算法中已取得一定的积极进展[14-16]，但是传统TV是基于图像的梯度模稀疏这一假设前提的，对于细节和纹理信息较少的图像具有较好的重构结果；但对于细节和纹理信息丰富的CMR图像而言，传统TV时常导致过度平滑、边缘模糊、细节失真、阶梯效应等缺点。本研究提出图像序列差分（total variation of image sequence，TVIS）正则项的思路，以图像序列的第1帧作为参考图像，其余后续帧的重构均以此作为比较，其定义见公式（4）。

其中，r表示参考图像，即整个图像序列的第 1帧图像；∇x和∇y分别表示沿x和y方向的梯度。比较公式（3）和公式（4），与传统TV方法不同的是，本研究提出的TVIS可以同时利用时间域和空间域的稀疏性，在CMR图像的实时重构中应该能发挥更好的作用。因此，基于TVIS模型的压缩感知重构算法的一般形式可写作公式（5）。

1.3 改进的 NESTA算法 NESTA算法是在 Yuri Nesterov思想的基础上拓展而来的，是属于迫近算法的一阶重构算法，可应用于求解有约束或无约束的优化问题，理论收敛速度快，重构结果精确，鲁棒性好，而且算法所需调节参数少[17-18]。NESTA算法用于求解如下问题：

以下重点研究yk和zk的计算。因为只要计算出它们的值，xk便是两者的加权平均。首先，yk的求解因其目标函数是凸函数，可行域是凸集，因此可采用标准的拉格朗日乘法实现。构造拉格朗日函数是拉格朗日问题的对偶解，则满足 Karush-Kuhn-Tucker（KKT）条件：

根据公式（8）的第4个条件和Sherman-Mprrison-Woodbury求逆公式，可得yk解如公式（9）所示。同理，zk解如公式（10）所示。将公式（9）、（10）代入改进的NESTA算法的计算步骤，便可完成NESTA算法的迭代过程。

综上所述，改进后的NESTA算法迭代步骤主要分成3步：求出的基础上，加权平均得到其核心思路是采用平滑技术，把不光滑的函数平滑化，然后再用高速便捷的迭代算法来求解。改进后的NESTA算法收敛速度可达级，基本逼近一阶算法的理论极限。

1.4 性能评价指标 采用峰值信噪比（peak signal to noise ratio，PSNR）、相对误差（RLNE）、均方根误差（root mean square error，RMSE）3个指标定量分析和评价算法的性能特点，其中PSNR是一种评价图像品质的客观标准，其值越大表示图像失真越小；RLNE用于比较完整采样图像和重构图像的相对l2范数误差，其值越小越好；RMSE说明完整采样图像和重构图像之间样本的离散程度，能够较好地反映重构图像的精密程度，其值越小越好。上述3个指标的计算分别见公式（11）～（13）。

2 结果

本实验在HP工作站Z600上采用MATLAB 2016编程实现，硬件配置为：Intel Xeon CPU E5606@2.13GHz，24 GB内存，64位Windows 7操作系统。使用的完整采样数据来自 kt FOCUSS程序中的full_rad数据，大小为384 384 32× ×。该数据集采用径向采样模式，沿通过K空间中心的径向直线采集数据，具有以下优点：①径向采样数据的每条线含有等量的低频到高频信息，有利于图像的欠采样重构；②径向采样模式对决定图像主要信息的中心数据过采样，因此径向釆样对物体运动没有笛卡尔采样那样敏感，更适用于CMR图像重构。

2.1 对比算法及参数设置 对比算法包括：直接填0的DFFT算法、kt FOCUSS算法[8]、基于运动补偿的kt FOCUSS 算法[9]、kt SLR 算法[10]、DTV 算法[11]。各个算法的相关参数设置见表1。

表1 本实验算法与5种对比算法的相关参数设置

2.2 算法精确性分析 本研究算法中，由于第1帧图像将作为参考图像，后续帧图像的压缩感知重构均以此作为参照，因此第1帧图像的欠采样率不能过低，可选定为40%左右，后续各帧的欠采样率设定为 12%。完整采样图像序列的第 16帧图像见图2A，其中红色正方形方框放大见图2B。本研究算法和上述对比算法同时选取第 16帧及其局部区域放大分别见图3、4，其中A～F依次为直接填0的DFFT算法、kt FOCUSS算法、基于运动补偿的kt FOCUSS算法、kt SLR算法、DTV算法和本研究算法结果图。

将本研究算法和上述对比算法所得结果，沿图2A中黄色虚线展开可得图5。

图2 完整采样图像第 16帧及其局部区域放大显示。A为完整采样图像第 16帧，其中红色虚线框放大为B图

图3 本研究算法和其他对比算法第16帧重构结果。A～F依次为直接填0的DFFT算法、kt FOCUSS算法、基于运动补偿的kt FOCUSS算法、kt SLR算法、DTV算法、本研究算法

图4 本研究算法和对比算法第16帧重构结果局部区域放大图。A～F依次为直接填0的DFFT算法、kt FOCUSS算法、基于运动补偿的kt FOCUSS算法、kt SLR算法、DTV算法、本研究算法

图5 本研究算法和对比算法沿图2黄色虚线展开图。A～G依次为完整采样图像、直接填0的DFFT算法、kt FOCUSS算法、基于运动补偿的kt FOCUSS算法、kt SLR算法、DTV算法、本研究算法。与完整采样图像相比，直接填0的DFFT方法展开图中存在不少伪影；箭所指位置，TV方法的展开图中缺失了上方的灰色条带，kt FOCUSS和基于运动补偿的kt FOCUSS算法中上方的灰色条带不够连续，kt SLR算法中下方箭所指位置的灰色条带明显偏暗

由图2～4可见，本研究算法得到的重构结果人眼视觉效果最佳，尤其在图4中标出的圆形和矩形框中，黑色像素点的位置及大小与图2B中所示的完整采样图像放大结果最接近。相比之下，直接填 0的DFFT方法、kt SLR算法和TV算法放大区域明显缺少黑色像素点聚集，而kt FOCUSS和基于运动补偿的 kt FOCUSS算法圆形框中的明亮像素明显缺失。本研究算法所得重构图像在人眼视觉系统的定性评价上性能最优，最接近完整采样图像。

从定量指标上进行分析，上述算法32帧图像的性能指标对比见图6，性能指标的平均值见表2。由于直接填0的DFFT方法效果最差，PSNR最低，因此在RLNE和RMSE的对比图上未参与显示，重点比较其他几种算法的RLNE和RMSE指标。本研究算法几乎在各帧图像重构中均能得到最高的PSNR值和最小的RMSE值，32帧的性能指标平均值也是如此；仅RLNE指标，是基于运动补偿的kt FOCUSS算法性能最优，可能因为在运动补偿环节对重构图像进行了必要的修正所致。因此，从性能指标定量分析来看，本研究算法总体性能仍为最佳（图6、表2）。

图6 本研究算法和对比算法的性能指标对比。A为PSNR指标，B为RLNE指标，C为RMSE指标

表2 本研究算法和对比算法的性能指标平均值

2.3 首帧图像欠采样率对算法性能的影响 如前所述，第 1帧图像作为整个图像序列的参考图像，其重构精度将显著影响后续帧的重构精度，因此有必要考查首帧图像欠采样率对整个算法性能的影响。首先，固定第1帧图像的欠采样率为40%，后续帧的欠采样率分别设置为10%、15%、20%、25%、40%，算法性能指标见图7。然后，将第1帧图像的欠采样率分别设置为100%、75%、50%、25%、15%，而将后续帧的欠采样率固定为 25%，算法的性能指标见图8。

图7 首帧图像欠采样率对算法性能的影响。固定第1帧图像欠采样率为40%，后续帧欠采样率可调。A为PSNR指标，B为RLNE指标，C为RMSE指标

图8 首帧图像欠采样率对算法性能的影响。第1帧图像欠采样率可调，后续帧欠采样率固定为25%。A为PSNR指标，B为RLNE指标，C为RMSE指标

由图7可见，首帧图像欠采样率固定为40%，即使后续帧欠采样率仅10%时，本研究算法依然能取得较好的性能指标，其PSNR均值为44.6268 dB，RLNE均值0.0374，RMSE均值0.0059，完全能够满足临床辅助诊疗的需要。随着后续帧欠采样率上升，各项性能指标逐步趋优，当后续帧欠采样率也达40%时，性能指标最高可达PSNR均值为49.3388 dB，RLNE均值0.0185，RMSE均值0.0034。由图8可见，当固定后续帧欠采样率为25%时，首帧图像的欠采样率对图像重构质量的影响不是特别明显。若首帧图像完整采样，则此时的性能指标最高可达 PSNR均值为 49.1890 dB，RLNE均值0.0149，RMSE均值0.0035；若首帧图像欠采样率仅为15%，此时的性能指标也能达到PSNR均值为45.9503 dB，RLNE均值0.0303，RMSE均值0.0051。综上所述，为使本研究算法的精确度能够达到临床辅助诊疗的作用，可以将首帧图像欠采样率设置为25%～40%，后续帧的欠采样率设置为15%～20%，这样即能确保RLNE和RMSE分别小于0.05和0.005。2.4 算法实时性分析 本研究算法和对比算法的运行时间结果见表2。相比而言，尽管直接填0的DFFT算法运行时间短到可以忽略不计，但由于其性能太差，无法应用于医学图像临床辅助诊疗，因此排除在外。其余算法中，kt SLR的运行时间达2977 s，约需50 min；kt FOCUSS算法运行222 s，基本满足实时性的要求；基于运动补偿的kt FOCUSS算法因运动补偿环节，其运行时间达882 s，即15 min左右；本研究算法运行时间仅36 s，平均每帧图像重构仅需1.125 s，完全满足医学图像重构实时性的要求。

3 讨论

本研究将压缩感知理论引入 CMR，提出基于图像序列差分正则项的实时CMR图像重构算法，在取得精确重构结果的同时，提高了程序运行的实时性。

在精确性方面，首帧图像欠采样率40%，后续各帧欠采样率12%时，相比其他5种对比算法，本研究算法获得最佳的性能评价指标，同时局部区域放大后黑色像素点的位置和大小最接近完整采样图像，即人眼视觉系统的主观感受最佳。

首帧图像欠采样率对算法性能影响不大。一般来说，首帧图像欠采样率为 25%～40%，后续帧欠采样率为 15%～20%，即能取得令人满意的重构结果，完全能够满足临床辅助诊断的需求。

实时性方面，对于实验使用的32帧CMR图像而言，本研究算法重构仅需36 s，比DTV算法运行时间略短，与kt FOCUSS算法、基于运动补偿的kt FOCUSS算法以及kt SLR算法相比，本研究算法的运行时间仅为其1/10甚至1/100。因此，本研究算法可以在获得精确重构结果的同时满足实时性的要求，在临床上具有广阔的应用前景。

当然，本研究算法尚有进一步提升的空间。近年来，基于低秩稀疏分解的CS-MRI图像重构研究方兴未艾。有学者首次将低秩稀疏分解模型应用于MRI图像重构，初步应用于心脏灌注、心脏电影、血管成像、胸腹部 MRI等[19]。有学者利用鲁棒主成分分析提取动态 MRI中的动态组织部分，提出基于稀疏和低秩先验分离的重构算法[20]。有人提出基于部分可分离模型和压缩感知相结合的Stepped-Sparse PS算法[21]。还有研究者初步提出多尺度和多维低秩矩阵分解，并应用于人脸图像、监控视频、动态 MRI等的多尺度建模和协同滤波[22-23]。

今后的研究方向包括：①将本研究算法与基于低秩稀疏分解算法进行性能对比和分析；②将本研究算法应用于更多CMR序列的快速重构，检验算法的普适性；③将本研究算法应用于儿童心肌病、复杂先天性心脏病、心脏二尖瓣病变图像重构中，并与临床医师的主观判断相对照，从而检验算法在临床疾病辅助诊疗中的实际作用。

综上所述，压缩感知理论引入 CMR，可以加快成像速度，提高成像质量，提高检查舒适度和检查质量，促进循环系统疾病诊疗新技术的发展。本研究结合串行和并行算法的各自优点，提出基于图像序列差分正则项的实时CMR图像重构算法，基本思路是以第1帧图像作为参考图像，从K空间数据重构出第1帧图像，此后利用图像序列差分正则化使得后续帧的重构均以第1帧图像为参考，从而同时利用时间域和空间域的稀疏性。采用改进的NESTA算法对模型进行求解，提高算法的实时性。以PSNR、RLNE和RMSE作为定量评价指标，结合人眼视觉感受、局部区域放大和沿时间轴展开等定性分析方法，充分说明本研究算法较kt FOCUSS、基于运动补偿的kt FOCUSS、kt SLR、DTV等主流算法的性能优势，能在获取较为精确的重构结果的同时，很好地满足医学图像处理实时性的要求。