基于低秩稀疏分解快速算法的动态MRI重建

杨青海 杨敏

摘  要：通过低秩加稀疏矩阵分解模型重建欠采样动态磁共振图像时，常采用变量分裂算法来求解。针对共轭梯度法在二次项更新中迭代计算较为复杂的问题，为了加快重建速度，提出一种考虑数据采集算子形式的高效变量分裂方案，将数据采集算子根据欠采样掩码矩阵、傅里叶变换算子和线圈灵敏度矩阵进行拆分，简化算法子问题中二次项更新所涉及的矩阵逆运算，达到加快算法收敛速度的目的。仿真实验结果表明：与迭代软阈值法和共轭梯度法相比，所提算法在心电影数据集中收敛速度分别提高了57.9%和83.0%，结构相似性分别提升了3.3%和1.4%;在心脏灌注数据集中收敛速度分别提高了55.5%和79.6%，结构相似性分别提升了1.5%和0.4%。

关键词：动态磁共振成像;压缩感知;低秩稀疏分解;变量分裂

中图分类号：TP391     文献标识码：A

Dynamic MRI Reconstruction based on Low-rank

Sparse Decomposition Fast Algorithm

YANG Qinghai， YANG Min

（School of Automation， Nanjing University of Posts and Telecommunications， Nanjing 210023， China）

1789321617@qq.com; yangm@njupt.edu.cn

Abstract： Variable splitting algorithm is often used to solve the problem of under-sampled dynamic MRI （Magnetic Resonance Imaging） reconstruction by low rank and sparse matrix decomposition model. Aiming at the complex iterative calculation of conjugate gradient method in quadratic term updating， in order to speed up the reconstruction， this paper proposes an efficient variable splitting scheme considering the form of the data acquisition operator. The data acquisition operator is divided according to the under-sampled mask code matrix， the Fourier transform operator and the coil sensitivity matrix， which simplifies the matrix inverse operation involved in the update of the quadratic term in the algorithm sub-problem， and achieves the purpose of accelerating the convergence speed the algorithm. Simulation results show that compared with the iterative soft threshold method and the conjugate gradient method， the convergence speed of the proposed algorithm in the cardiac cine dataset has increased by 57.9% and 83.0% respectively， and the structural similarity has increased by 3.3% and 1.4% respectively. In the cardiac perfusion dataset， the convergence speed has increased by 55.5% and 79.6% respectively， and the structural similarity has increased by 1.5% and 0.4% respectively.

Keywords： dynamic magnetic resonance imaging; compressed sensing; low rank sparse decomposition; variable

splitting

1   引  言（Introduction）

對于动态磁共振图像，可以把同一层空间不同时间帧的图像看作一个列向量，作为时空矩阵的一列，由多个时间帧构建成一个低秩矩阵，从而把动态磁共振成像（MRI）重建问题转化为低秩矩阵恢复问题。OTAZO等[1]和GAO等[2]基于L+S模型，通过IST算法对临床医学图像进行重建，并取得了可观的重建效果。此前的大量研究表明，基于L+S模型对动态MRI进行重建整体上可以取得较好的重建精度，但重建速度仍有提升空间。

本文考虑到在常规的变量分解优化方案中，二次项的更新一般需要共轭梯度法来迭代求解，提出基于L+S框架下结合数据采集算子的快速变量分裂法。该方法采用与欠采样、傅里叶编码和灵敏度映射相关联的矩阵结构，从而加快重建速度。实验表明，此方法能有效实现动态MRI，重建速度较快。

2   研究现状（Research status）

压缩感知（CS）已在MRI中得到广泛应用，提高了数据采集效率[3-4]。由于动态磁共振图像具有高度的时空相关性，CS-MRI模型[5]同样可以用于动态MRI重建。此外，压缩感知与灵敏度编码（SENSitivity Encoding， SENSE）[6]等并行MRI技术相结合，通过多个线圈收集更多的数据，达到改善重建图像的时空分辨率平衡的效果。JUNG等根据FOCUSS算法在时间变换域进行稀疏约束，提出包含运动估计和补偿的k-t FOCUSS[7]方法，成功应用于心电影MRI重建中。但在非周期运动情况下，稀疏化残差信号阻碍了预测方案的发展。近年来，研究人员不再简单地利用向量的稀疏性，同时也在矩阵的低秩性方面做了大量工作。LINGALA等提出k-t SLR算法[8]，利用卡-洛变换（Karhunen-Louve Transform， KLT）的低秩先验和全局稀疏性进行动态MRI重建。但该算法并没有考虑磁共振图像的结构稀疏性，限制了算法的发展。同时，有研究提出基于patch的字典学习方法用于动态MRI重建[9-10]。然而由于动态MRI序列一般都比较大，字典学习对大数据集的效率很低，重建时间过长，也有通过局部低秩再加稀疏约束（LLRS）来进行动态MRI重建[11]，但对局部块的大小要求较高，重建精度有待进一步提高。压缩感知和低秩矩阵结合的思想为动态MRI重建提供了新方向。文献[12]提出将原始数据矩阵拆分成一个低秩矩阵和稀疏矩阵的模型用于解决鲁棒主成分分析（RPCA）问题。

但是，目前基于L+S框架的重建方法在求解过程中的二次项求解涉及多次迭代计算，虽然整体上能取得较好的重建质量，但在重建速度上仍有较大提升空间。针对此问题，本文提出一种结合数据采集算子进行变量分裂的方法优化分裂模型，该方法能有效重建动态磁共振图像。

3  L+S矩阵分解模型及其解法（L+S matrix factorization model and its solution）

3.1   L+S分解模型

低秩加稀疏分解模型的目的是将原始输入数据矩阵分解成一个低秩矩阵和一个稀疏矩阵相叠加的形式。可以用经典的鲁棒主成分分析法求解此类问题。

（1）

式（1）中的优化问题是NP（非确定多项式）难问题，需要利用凸松弛将非凸函数变换为凸函数，矩阵的秩用核范数代替，范数用范数代替[13]。由式（1）可以转化为如下凸优化问题：

（2）

式（2）中，代表核范数，代表矩阵的范数。

在动态MRI重建中，可把图像序列分解成静态的背景和动态的前景两个部分，利用相应的低秩及稀疏行先验知识，通过鲁棒主成分分析法分别进行重建，再叠加得到重建图像。将L+S分解模型应用于动态MRI重建的前提条件是能稳定地区分图像序列的背景和动态部分，需要满足“不相关”准则，即代表背景的低秩部分是不稀疏的或者说稀疏部分是非低秩的。在实际的医学图像运用中，可能无法严格地满足不相关性，但由于的秩通常远小于矩阵的秩，且的奇异值远大于，因此大部分的背景信息仍保留在低秩矩阵中。

L+S模型對于动态MRI，目标是恢复一个未知的图像。文献[1]采用以下正则化凸优化方案：

（3）

式（3）中，是考虑了线圈灵敏度和欠采样的傅里叶变换的数据采集算子，相当于对图像的空间编码;是基于图像稀疏域先验假设的稀疏变换算子，这种稀疏变换已广泛应用于动态MRI重建的稀疏化[14-15];是欠采样的空间数据;和是数据平衡参数，用以平衡范数、核范数和范数。

3.2   变量分裂方案

经验表明，基于增广拉格朗日（Augmented Lagrangian， AL）框架的变量分裂法与近端梯度法（PGM）相比，可以在更少的迭代中达到更高的精度。文献[12]中通过变量分裂方案来求解L+S分解问题，用辅助变量U、W作为约束条件，将式（3）重新表述为：

（4）

其中，。再用乘子法对式（4）进行无约束转化，得到相应的增广拉格朗日函数为：

（5）

其中，V1、V2是拉格朗日乘子数组，、是相应的AL惩罚参数。利用交替方向法（Alternating Direction Method of Multipliers， ADMM）可将式（5）转化为四个变量子问题的求解。和的求解涉及范数，计算过程中涉及项，由于在并行MRI中数据采集算子通常包含线圈敏感性的额外信息，因此的求解需要用到共轭梯度法。考虑到此迭代方法计算复杂，本文提出基于L+S模型的计算效率更高的变量分裂方案，提高并行MRI的计算速率。

4  基于变量分裂的加速方案（Acceleration scheme based on variable split）

针对L+S模型下动态并行MRI重建问题，本文提出基于数据采集算子的分裂方案。数据采集算子，其中为所有帧的欠采样模式，为傅里叶编码矩阵，为接收线圈的灵敏度。

式（3）中，和的更新是二次的，需要计算。由于在并行MRI中是不循环的，因此涉及的矩阵逆的求解使用的是共轭梯度法，计算复杂，重建速度较慢。本文方法在常规的变量分裂方案框架下，将数据采集算子根据欠采样掩码矩阵、傅里叶变换算子和线圈灵敏度矩阵进行拆分，虽然同样也会涉及对四个子问题的求解，但在变量分裂框架中利用数据采集算子的等价公式，将采样矩阵、傅里叶编码矩阵以及灵敏度矩阵带入分裂方案中，可以把更新中涉及的替换成，其中欠采样掩码矩阵用克罗内克积表示。因为这里的是对角矩阵，所以的计算更为简单，可以大大加快二次项的求解速度。同时，因为模型考虑了时间傅里叶变换的稀疏性，使得算法可以规范线圈灵敏度矩阵，让归一化，在不影响低秩分量的秩的情况进一步加快算法的重建速度。被规范化为：。在此条件下，式（4）可以被约束表示为：

（6）

其中，，。结合拉格朗日乘子，利用乘子法对式（6）进行无约束转化，得到修正的AL函数：

（7）

的更新涉及核范数，其近端映射通过奇异值阈值求解：

（8）

的更新涉及范数，通过软阈值得出。在这里使用的是一个酉算子，设置变量，则的更新表示为：

（9）

的更新表示为：

（10）

的更新表示为：

（11）

的更新都涉及二次项，其中是傅里叶编码矩阵，且。

5  仿真实验与分析（Simulation experiment and analysis）

5.1   评价指标

为了验证本文提出的方法，分别用迭代收缩阈值算法（Iterative Shrinkage Thresholding Algorithm，ISTA）、涉及共轭梯度法的AL-CG和本文提出算法对心电影成像和心脏灌注数据集进行重建。对于各算法的重建结果，除主观判断外，本文通过计算每一次迭代收敛图像的标准化均方根差（NRMSD）来验证其收敛速度，定义为：

（12）

其中，。

用图像的结构相似性（SSIM）来度量每个时间帧的重建图像与全采样图像间的差值，定义为：

（13）

SSIM值的范围是[0，1]，当值为1时，代表两张图完全一致。

用图像的峰值信噪比（PSNR）衡量重建图像的精度，定义为：

（14）

其中，为原始图像，为重建图像，为图像大小。

5.2   数据集及参数设置

实验在i5-10210U CPU、Windows 10操作系统的笔记本下采用MATLAB R2020a进行仿真。实验中采用笛卡尔采样模式，为了保证能快速收敛，对于两个数据集ISTA的步长都设置为0.99。

心电影数据集大小为256×256、24帧、12线圈，采样加速因子为8。实验过程中设置=0.01，=0.025;将采用共轭梯度方案的分裂方案记为AL-CG，设置=0.1，=0.05;将本文所提方案记为AL-E，设置=0.1，=0.01。

心脏灌注数据集大小为128×128、40帧、12线圈，采样加速因子为10。参数设置为=0.01;对于AL-CG算法，设置=0.2;对于AL-E算法，设置=0.1，=0.02。

5.3   实验结果

图1（a）为心电影数据集三种算法分别运行3，500 s的收敛性分析，图1（b）为心脏灌注数据集三种算法分别运行2，000 s的收敛性分析。图1中用每100 次迭代标记点来体现算法的相对速度。从图1中可以看出，AL-CG是三种算法中收敛最慢的，由于步长选择的原因，ISTA在开始阶段比其他两种方法收敛更快，但AL-E在总体上是最快的。

图2为心电影成像分别抽取第2、8、14、20 帧的重建结果;图3是心脏灌注成像分别抽取第5、15、25、35 帧的重建结果。

从图3可以看出，本文方法重建的图像在低秩部分比AL-CG和ISTA更加清晰。如圖3（a）箭头区域所示，本文方法器官边界信息的重建更加明显，细节部分的重建精度更高。

图4和图5为三种方法在两个数据集下重建结果的展示。结合两组数据集的结果，从总体上看三种算法的SSIM在两个数据集中差异不大，都在一定区间内波动。可以看到AL-E算法的SSIM值总体上比ISTA和AL-CG高，仅在个别帧下效果有浮动，其重建还原度比另两种算法更高。虽然PSNR相较常规的共轭梯度法略低，但在结构相似性上表现效果更好。

表1是三种方法在两个数据集下分别迭代50 次所耗时间及所有时间帧的峰值信噪比和结构相似性的平均值。由于步长选择的原因，ISTA方法在较少的迭代次数下会有较快的收敛速度，总体上本文所提方法相较共轭梯度法在重建速度上有较大提升。

6   结论（Conclusion）

本文在低秩稀疏矩阵分解模型的基础上，提出考虑数据采集算子形式的变量分裂方案，将算法二次项更新时所涉及的矩阵逆运算简化，提高算法的收敛速度。通过与AL-CG算法以及ISTA算法的仿真结果对比，本文提出的算法重建速度更快，且边缘信息的重建效果更好。但是本文方法需要额外优化两个拉格朗日参数，在收敛速度上仍有较大的提升空间，有待进一步研究。

参考文献（References）

[1] OTAZO R， CANDES E， SODICKSON D K. Low-rank plus sparse matrix decomposition for accelerated dynamic MRI with separation of background and dynamic components[J]. Magnetic Resonance in Medicine， 2015， 73（3）：1125-1136.

[2] GAO H， LI L， HU X. Compressive diffusion MRI[C]// Iidaka T， Miyakoshi M， Harada T， et al. Proceedings of the 21th Annual Meeting of International Society for Magnetic Resonance in Medicine. Salt Lake， USA： ISMRM， 2013：610.

[3] LUSTIG M， DONOHO D， PAULY J M. Sparse MRI： The application of compressed sensing for rapid MR imaging[J]. Magnetic Resonance in Medicine， 2007， 58（6）：1182-1195.

[4] RAVISHANKAR S， BRESLER Y. MR image reconstruction from highly undersampled k-space data by dictionary learning[J]. IEEE Transactions Med Imaging， 2011， 30（5）：

1028-1041.

[5] TREMOULHEAC B， DIKAIOS N， ATKINSON D， et al.

Dynamic MR image reconstruction-separation from undersampled k-t-space via low rank plus sparse prior[J]. IEEE Transactions Med Imaging， 2014， 33（8）：1689-1701.

[6] PRUESSMANN K P， WEIGER M， SCHEIDEGGER M B，

et al. Sense： Sensitivity encoding for fast MRI[J]. Magnetic Resonance in Medicine， 1999， 42（5）：952-962.

[7] JUNG H， SUNG K， NAYAK K S， et al. K-T Focuss： A general compressed sensing framework for high resolution dynamic MRI[J]. Magnetic Resonance in Medicine， 2009， 61（1）：

103-116.

[8] LINGALA S G， HU Y， DIBELLA E， et al. Accelerated dynamic MRI exploiting sparsity and low-rank structure： K-t SLR[J]. IEEE Transactions on Medical Imaging， 2011， 30（5）：

1042-1054.

[9] AWATE S P， DIBELLA E V R. Spatiotemporal dictionary learning for undersampled dynamic MRI reconstruction via joint frame-based anddictionary-based sparsity[C]// Alejandro F，

Andrés S， Raimund O， et al. 2012 9th IEEE International Symposium on Biomedical Imaging （ISBI）. Barcelona， Spain：

IEEE， 2012：318-321.

[10] CABALLERO J， PRICE A N， RUECKERT D， et al. Dictionary learning and time sparsity for dynamic MR data reconstruction[J]. IEEE Transactions on Medical Imaging， 2014， 33（4）：979-994.

[11] CHANDRASEKARAN V， SANGHAVI S， PARRILO P A， et al. Rank-sparsity incoherence for matrix decomposition[J]. SIAM Journal on Optimization， 2011， 21（2）：572-596.

[12] KAFALI S G， SHIH S F， RUAN D， et al. Adaptive locally low rank and sparsity constrained reconstruction for accelerated dynamic MRI[C]// MATHEWS J， JONG C Y. 2020 IEEE 17th International Symposium on Biomedical Imaging （ISBI）. Iowa， USA： IEEE， 2020：930-934.

[13] 馬杰，王晓云，张志伟，等.一种基于全变分正则化低秩稀疏分解的动态MRI重建方法[J].光电子·激光，2016，27（1）：87-96.

[14] LUSTIG M， PAULY J M. Spirit： Iterative self-consistent parallel imaging reconstruction from arbitrary k-space[J]. Magnetic Resonance in Medicine， 2010， 64（2）：457-471.

[15] 史加荣，郑秀云，周水生.矩阵补全算法研究进展[J].计算机科学，2014，41（4）：13-20.

作者简介：

杨青海（1995-），男，硕士生.研究领域：图像处理.

杨  敏（1969-），男，博士，副教授.研究领域：图像处理与模式识别.