基于低秩稀疏分解快速算法的動態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-），男，博士，副教授.研究領域：圖像處理與模式識別.