楊青海 楊敏

摘 ?要：通過(guò)低秩加稀疏矩陣分解模型重建欠采樣動(dòng)態(tài)磁共振圖像時(shí)，常采用變量分裂算法來(lái)求解。針對(duì)共軛梯度法在二次項(xiàng)更新中迭代計(jì)算較為復(fù)雜的問(wèn)題，為了加快重建速度，提出一種考慮數(shù)據(jù)采集算子形式的高效變量分裂方案，將數(shù)據(jù)采集算子根據(jù)欠采樣掩碼矩陣、傅里葉變換算子和線圈靈敏度矩陣進(jìn)行拆分，簡(jiǎn)化算法子問(wèn)題中二次項(xiàng)更新所涉及的矩陣逆運(yùn)算，達(dá)到加快算法收斂速度的目的。仿真實(shí)驗(yàn)結(jié)果表明：與迭代軟閾值法和共軛梯度法相比，所提算法在心電影數(shù)據(jù)集中收斂速度分別提高了57.9%和83.0%，結(jié)構(gòu)相似性分別提升了3.3%和1.4%;在心臟灌注數(shù)據(jù)集中收斂速度分別提高了55.5%和79.6%，結(jié)構(gòu)相似性分別提升了1.5%和0.4%。

關(guān)鍵詞：動(dòng)態(tài)磁共振成像;壓縮感知;低秩稀疏分解;變量分裂

中圖分類號(hào)：TP391 ? ? 文獻(xiàn)標(biāo)識(shí)碼：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）

對(duì)于動(dòng)態(tài)磁共振圖像，可以把同一層空間不同時(shí)間幀的圖像看作一個(gè)列向量，作為時(shí)空矩陣的一列，由多個(gè)時(shí)間幀構(gòu)建成一個(gè)低秩矩陣，從而把動(dòng)態(tài)磁共振成像（MRI）重建問(wèn)題轉(zhuǎn)化為低秩矩陣恢復(fù)問(wèn)題。OTAZO等[1]和GAO等[2]基于L+S模型，通過(guò)IST算法對(duì)臨床醫(yī)學(xué)圖像進(jìn)行重建，并取得了可觀的重建效果。此前的大量研究表明，基于L+S模型對(duì)動(dòng)態(tài)MRI進(jìn)行重建整體上可以取得較好的重建精度，但重建速度仍有提升空間。

本文考慮到在常規(guī)的變量分解優(yōu)化方案中，二次項(xiàng)的更新一般需要共軛梯度法來(lái)迭代求解，提出基于L+S框架下結(jié)合數(shù)據(jù)采集算子的快速變量分裂法。該方法采用與欠采樣、傅里葉編碼和靈敏度映射相關(guān)聯(lián)的矩陣結(jié)構(gòu)，從而加快重建速度。實(shí)驗(yàn)表明，此方法能有效實(shí)現(xiàn)動(dòng)態(tài)MRI，重建速度較快。

2 ? 研究現(xiàn)狀（Research status）

壓縮感知（CS）已在MRI中得到廣泛應(yīng)用，提高了數(shù)據(jù)采集效率[3-4]。由于動(dòng)態(tài)磁共振圖像具有高度的時(shí)空相關(guān)性，CS-MRI模型[5]同樣可以用于動(dòng)態(tài)MRI重建。此外，壓縮感知與靈敏度編碼（SENSitivity Encoding， SENSE）[6]等并行MRI技術(shù)相結(jié)合，通過(guò)多個(gè)線圈收集更多的數(shù)據(jù)，達(dá)到改善重建圖像的時(shí)空分辨率平衡的效果。JUNG等根據(jù)FOCUSS算法在時(shí)間變換域進(jìn)行稀疏約束，提出包含運(yùn)動(dòng)估計(jì)和補(bǔ)償?shù)膋-t FOCUSS[7]方法，成功應(yīng)用于心電影MRI重建中。但在非周期運(yùn)動(dòng)情況下，稀疏化殘差信號(hào)阻礙了預(yù)測(cè)方案的發(fā)展。近年來(lái)，研究人員不再簡(jiǎn)單地利用向量的稀疏性，同時(shí)也在矩陣的低秩性方面做了大量工作。LINGALA等提出k-t SLR算法[8]，利用卡-洛變換（Karhunen-Louve Transform， KLT）的低秩先驗(yàn)和全局稀疏性進(jìn)行動(dòng)態(tài)MRI重建。但該算法并沒有考慮磁共振圖像的結(jié)構(gòu)稀疏性，限制了算法的發(fā)展。同時(shí)，有研究提出基于patch的字典學(xué)習(xí)方法用于動(dòng)態(tài)MRI重建[9-10]。然而由于動(dòng)態(tài)MRI序列一般都比較大，字典學(xué)習(xí)對(duì)大數(shù)據(jù)集的效率很低，重建時(shí)間過(guò)長(zhǎng)，也有通過(guò)局部低秩再加稀疏約束（LLRS）來(lái)進(jìn)行動(dòng)態(tài)MRI重建[11]，但對(duì)局部塊的大小要求較高，重建精度有待進(jìn)一步提高。壓縮感知和低秩矩陣結(jié)合的思想為動(dòng)態(tài)MRI重建提供了新方向。文獻(xiàn)[12]提出將原始數(shù)據(jù)矩陣拆分成一個(gè)低秩矩陣和稀疏矩陣的模型用于解決魯棒主成分分析（RPCA）問(wèn)題。

但是，目前基于L+S框架的重建方法在求解過(guò)程中的二次項(xiàng)求解涉及多次迭代計(jì)算，雖然整體上能取得較好的重建質(zhì)量，但在重建速度上仍有較大提升空間。針對(duì)此問(wèn)題，本文提出一種結(jié)合數(shù)據(jù)采集算子進(jìn)行變量分裂的方法優(yōu)化分裂模型，該方法能有效重建動(dòng)態(tài)磁共振圖像。

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

3.1 ? L+S分解模型

低秩加稀疏分解模型的目的是將原始輸入數(shù)據(jù)矩陣分解成一個(gè)低秩矩陣和一個(gè)稀疏矩陣相疊加的形式?？梢杂媒?jīng)典的魯棒主成分分析法求解此類問(wèn)題。

（1）

式（1）中的優(yōu)化問(wèn)題是NP（非確定多項(xiàng)式）難問(wèn)題，需要利用凸松弛將非凸函數(shù)變換為凸函數(shù)，矩陣的秩用核范數(shù)代替，范數(shù)用范數(shù)代替[13]。由式（1）可以轉(zhuǎn)化為如下凸優(yōu)化問(wèn)題：

（2）

式（2）中，代表核范數(shù)，代表矩陣的范數(shù)。

在動(dòng)態(tài)MRI重建中，可把圖像序列分解成靜態(tài)的背景和動(dòng)態(tài)的前景兩個(gè)部分，利用相應(yīng)的低秩及稀疏行先驗(yàn)知識(shí)，通過(guò)魯棒主成分分析法分別進(jìn)行重建，再疊加得到重建圖像。將L+S分解模型應(yīng)用于動(dòng)態(tài)MRI重建的前提條件是能穩(wěn)定地區(qū)分圖像序列的背景和動(dòng)態(tài)部分，需要滿足“不相關(guān)”準(zhǔn)則，即代表背景的低秩部分是不稀疏的或者說(shuō)稀疏部分是非低秩的。在實(shí)際的醫(yī)學(xué)圖像運(yùn)用中，可能無(wú)法嚴(yán)格地滿足不相關(guān)性，但由于的秩通常遠(yuǎn)小于矩陣的秩，且的奇異值遠(yuǎn)大于，因此大部分的背景信息仍保留在低秩矩陣中。

L+S模型對(duì)于動(dòng)態(tài)MRI，目標(biāo)是恢復(fù)一個(gè)未知的圖像。文獻(xiàn)[1]采用以下正則化凸優(yōu)化方案：

（3）

式（3）中，是考慮了線圈靈敏度和欠采樣的傅里葉變換的數(shù)據(jù)采集算子，相當(dāng)于對(duì)圖像的空間編碼;是基于圖像稀疏域先驗(yàn)假設(shè)的稀疏變換算子，這種稀疏變換已廣泛應(yīng)用于動(dòng)態(tài)MRI重建的稀疏化[14-15];是欠采樣的空間數(shù)據(jù);和是數(shù)據(jù)平衡參數(shù)，用以平衡范數(shù)、核范數(shù)和范數(shù)。

3.2 ? 變量分裂方案

經(jīng)驗(yàn)表明，基于增廣拉格朗日（Augmented Lagrangian， AL）框架的變量分裂法與近端梯度法（PGM）相比，可以在更少的迭代中達(dá)到更高的精度。文獻(xiàn)[12]中通過(guò)變量分裂方案來(lái)求解L+S分解問(wèn)題，用輔助變量U、W作為約束條件，將式（3）重新表述為：

（4）

其中，。再用乘子法對(duì)式（4）進(jìn)行無(wú)約束轉(zhuǎn)化，得到相應(yīng)的增廣拉格朗日函數(shù)為：

（5）

其中，V1、V2是拉格朗日乘子數(shù)組，、是相應(yīng)的AL懲罰參數(shù)。利用交替方向法（Alternating Direction Method of Multipliers， ADMM）可將式（5）轉(zhuǎn)化為四個(gè)變量子問(wèn)題的求解。和的求解涉及范數(shù)，計(jì)算過(guò)程中涉及項(xiàng)，由于在并行MRI中數(shù)據(jù)采集算子通常包含線圈敏感性的額外信息，因此的求解需要用到共軛梯度法。考慮到此迭代方法計(jì)算復(fù)雜，本文提出基于L+S模型的計(jì)算效率更高的變量分裂方案，提高并行MRI的計(jì)算速率。

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

針對(duì)L+S模型下動(dòng)態(tài)并行MRI重建問(wèn)題，本文提出基于數(shù)據(jù)采集算子的分裂方案。數(shù)據(jù)采集算子，其中為所有幀的欠采樣模式，為傅里葉編碼矩陣，為接收線圈的靈敏度。

式（3）中，和的更新是二次的，需要計(jì)算。由于在并行MRI中是不循環(huán)的，因此涉及的矩陣逆的求解使用的是共軛梯度法，計(jì)算復(fù)雜，重建速度較慢。本文方法在常規(guī)的變量分裂方案框架下，將數(shù)據(jù)采集算子根據(jù)欠采樣掩碼矩陣、傅里葉變換算子和線圈靈敏度矩陣進(jìn)行拆分，雖然同樣也會(huì)涉及對(duì)四個(gè)子問(wèn)題的求解，但在變量分裂框架中利用數(shù)據(jù)采集算子的等價(jià)公式，將采樣矩陣、傅里葉編碼矩陣以及靈敏度矩陣帶入分裂方案中，可以把更新中涉及的替換成，其中欠采樣掩碼矩陣用克羅內(nèi)克積表示。因?yàn)檫@里的是對(duì)角矩陣，所以的計(jì)算更為簡(jiǎn)單，可以大大加快二次項(xiàng)的求解速度。同時(shí)，因?yàn)槟Ｐ涂紤]了時(shí)間傅里葉變換的稀疏性，使得算法可以規(guī)范線圈靈敏度矩陣，讓歸一化，在不影響低秩分量的秩的情況進(jìn)一步加快算法的重建速度。被規(guī)范化為：。在此條件下，式（4）可以被約束表示為：

（6）

其中，，。結(jié)合拉格朗日乘子，利用乘子法對(duì)式（6）進(jìn)行無(wú)約束轉(zhuǎn)化，得到修正的AL函數(shù)：

（7）

的更新涉及核范數(shù)，其近端映射通過(guò)奇異值閾值求解：

（8）

的更新涉及范數(shù)，通過(guò)軟閾值得出。在這里使用的是一個(gè)酉算子，設(shè)置變量，則的更新表示為：

（9）

的更新表示為：

（10）

的更新表示為：

（11）

的更新都涉及二次項(xiàng)，其中是傅里葉編碼矩陣，且。

5 ?仿真實(shí)驗(yàn)與分析（Simulation experiment and analysis）

5.1 ? 評(píng)價(jià)指標(biāo)

為了驗(yàn)證本文提出的方法，分別用迭代收縮閾值算法（Iterative Shrinkage Thresholding Algorithm，ISTA）、涉及共軛梯度法的AL-CG和本文提出算法對(duì)心電影成像和心臟灌注數(shù)據(jù)集進(jìn)行重建。對(duì)于各算法的重建結(jié)果，除主觀判斷外，本文通過(guò)計(jì)算每一次迭代收斂圖像的標(biāo)準(zhǔn)化均方根差（NRMSD）來(lái)驗(yàn)證其收斂速度，定義為：

（12）

其中，。

用圖像的結(jié)構(gòu)相似性（SSIM）來(lái)度量每個(gè)時(shí)間幀的重建圖像與全采樣圖像間的差值，定義為：

（13）

SSIM值的范圍是[0，1]，當(dāng)值為1時(shí)，代表兩張圖完全一致。

用圖像的峰值信噪比（PSNR）衡量重建圖像的精度，定義為：

（14）

其中，為原始圖像，為重建圖像，為圖像大小。

5.2 ? 數(shù)據(jù)集及參數(shù)設(shè)置

實(shí)驗(yàn)在i5-10210U CPU、Windows 10操作系統(tǒng)的筆記本下采用MATLAB R2020a進(jìn)行仿真。實(shí)驗(yàn)中采用笛卡爾采樣模式，為了保證能快速收斂，對(duì)于兩個(gè)數(shù)據(jù)集ISTA的步長(zhǎng)都設(shè)置為0.99。

心電影數(shù)據(jù)集大小為256×256、24幀、12線圈，采樣加速因子為8。實(shí)驗(yàn)過(guò)程中設(shè)置=0.01，=0.025;將采用共軛梯度方案的分裂方案記為AL-CG，設(shè)置=0.1，=0.05;將本文所提方案記為AL-E，設(shè)置=0.1，=0.01。

心臟灌注數(shù)據(jù)集大小為128×128、40幀、12線圈，采樣加速因子為10。參數(shù)設(shè)置為=0.01;對(duì)于AL-CG算法，設(shè)置=0.2;對(duì)于AL-E算法，設(shè)置=0.1，=0.02。

5.3 ? 實(shí)驗(yàn)結(jié)果

圖1（a）為心電影數(shù)據(jù)集三種算法分別運(yùn)行3，500 s的收斂性分析，圖1（b）為心臟灌注數(shù)據(jù)集三種算法分別運(yùn)行2，000 s的收斂性分析。圖1中用每100 次迭代標(biāo)記點(diǎn)來(lái)體現(xiàn)算法的相對(duì)速度。從圖1中可以看出，AL-CG是三種算法中收斂最慢的，由于步長(zhǎng)選擇的原因，ISTA在開始階段比其他兩種方法收斂更快，但AL-E在總體上是最快的。

圖2為心電影成像分別抽取第2、8、14、20 幀的重建結(jié)果;圖3是心臟灌注成像分別抽取第5、15、25、35 幀的重建結(jié)果。

從圖3可以看出，本文方法重建的圖像在低秩部分比AL-CG和ISTA更加清晰。如圖3（a）箭頭區(qū)域所示，本文方法器官邊界信息的重建更加明顯，細(xì)節(jié)部分的重建精度更高。

圖4和圖5為三種方法在兩個(gè)數(shù)據(jù)集下重建結(jié)果的展示。結(jié)合兩組數(shù)據(jù)集的結(jié)果，從總體上看三種算法的SSIM在兩個(gè)數(shù)據(jù)集中差異不大，都在一定區(qū)間內(nèi)波動(dòng)。可以看到AL-E算法的SSIM值總體上比ISTA和AL-CG高，僅在個(gè)別幀下效果有浮動(dòng)，其重建還原度比另兩種算法更高。雖然PSNR相較常規(guī)的共軛梯度法略低，但在結(jié)構(gòu)相似性上表現(xiàn)效果更好。

表1是三種方法在兩個(gè)數(shù)據(jù)集下分別迭代50 次所耗時(shí)間及所有時(shí)間幀的峰值信噪比和結(jié)構(gòu)相似性的平均值。由于步長(zhǎng)選擇的原因，ISTA方法在較少的迭代次數(shù)下會(huì)有較快的收斂速度，總體上本文所提方法相較共軛梯度法在重建速度上有較大提升。

6 ? 結(jié)論（Conclusion）

本文在低秩稀疏矩陣分解模型的基礎(chǔ)上，提出考慮數(shù)據(jù)采集算子形式的變量分裂方案，將算法二次項(xiàng)更新時(shí)所涉及的矩陣逆運(yùn)算簡(jiǎn)化，提高算法的收斂速度。通過(guò)與AL-CG算法以及ISTA算法的仿真結(jié)果對(duì)比，本文提出的算法重建速度更快，且邊緣信息的重建效果更好。但是本文方法需要額外優(yōu)化兩個(gè)拉格朗日參數(shù)，在收斂速度上仍有較大的提升空間，有待進(jìn)一步研究。

參考文獻(xiàn)（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] 馬杰，王曉云，張志偉，等.一種基于全變分正則化低秩稀疏分解的動(dòng)態(tài)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] 史加榮，鄭秀云，周水生.矩陣補(bǔ)全算法研究進(jìn)展[J].計(jì)算機(jī)科學(xué)，2014，41（4）：13-20.

作者簡(jiǎn)介：

楊青海（1995-），男，碩士生.研究領(lǐng)域：圖像處理.

楊 ?敏（1969-），男，博士，副教授.研究領(lǐng)域：圖像處理與模式識(shí)別.