基于L1/2正則項的磁共振圖像稀疏重構(gòu)

馬杰，葛嶺嶺，苑煥朝，張婷婷

（河北工業(yè)大學電子信息工程學院，天津 300401）

基于L1/2正則項的磁共振圖像稀疏重構(gòu)

馬杰，葛嶺嶺，苑煥朝，張婷婷

（河北工業(yè)大學電子信息工程學院，天津 300401）

磁共振圖像可以利用壓縮感知從數(shù)量非常有限的觀測數(shù)據(jù)集合中重構(gòu)出，然而為了能夠做到這一點，必須要解決定義在大量數(shù)據(jù)集合上的非光滑函數(shù)的最小化這一困難問題．通常L1范數(shù)能夠產(chǎn)生稀疏解,但它往往與真實稀疏解（L0的解）差距甚大．針對該問題，研究一種基于變量分裂的圖像重構(gòu)模型，引入待重構(gòu)圖像的L1/2范數(shù)作為新正則項，采用交替增廣拉格朗日乘子法進行求解．為考察方法的穩(wěn)定性和重構(gòu)效果，結(jié)合不同參數(shù)等評價標準與現(xiàn)有的圖像重構(gòu)模型進行比較．實驗結(jié)果表明，L1/2范數(shù)作為正則子的圖像重構(gòu)模型相對于原有模型，圖像重構(gòu)結(jié)果穩(wěn)定性好，可以獲得更高的信噪比．

磁共振；壓縮感知；L1范數(shù)；L1/2范數(shù)；變量分裂

壓縮感知[1]（CS）是Donoho和Candes等提出的一種新興的信號獲取與處理理論．如果圖像在某個變換域可以稀疏表示，那么可以通過求解相關(guān)的優(yōu)化問題，就可由隨機采樣的稀疏系數(shù)來進行重構(gòu)，并在一定程度上保持原有圖像的重構(gòu)質(zhì)量．核磁共振成像（簡稱NMRI），也稱磁共振成像（MagneticResonanceImaging，簡稱MRI），是利用核磁共振原理，繪制物體內(nèi)部的結(jié)構(gòu)圖像．傅里葉變換法是最基本的MRI重建方法之一，但由于該方法運算量大，重建速度慢，實際中很少采用這種方法．基于壓縮感知的MRI重構(gòu)算法利用MRI稀疏表示或局部光滑的先驗知識，通過求解相應(yīng)的優(yōu)化問題來實現(xiàn)圖像重構(gòu)．目前已有多種算法解決此類優(yōu)化問題，例如Trzasko等提出SL0算法[2]（smooth L0norm），它利用圖像小波變換域稀疏性，用光滑L0范數(shù)近似L0范數(shù)重構(gòu)圖像，有效解決了范數(shù)優(yōu)化中的NP難問題，Bhaskar提出當0＜p＜1時基于Lp范數(shù)的迭代加權(quán)最小二乘（iterativelyreweightedleastsquare，IRLS）算法[3]，實驗證明了基于Lp范數(shù)優(yōu)化算法在重構(gòu)信號效果及可靠性方面優(yōu)于L1及L0范數(shù)優(yōu)化算法，Bioucas-Dias等提出TVMM（Total variation based majorization minimization）算法[4]，該算法利用圖像局部光滑特性，采用全變分正則化重構(gòu)圖像．文獻[5]采用非線性共軛梯度下降算法求解加權(quán)范數(shù)最小化問題，在迭代過程中，根據(jù)所求得的圖像稀疏表示來更新權(quán)值矩陣，增強MR圖像的稀疏性，提高了圖像的重建質(zhì)量．四川大學的李青[6]提出了一種基于壓縮感知的自適應(yīng)正則化MRI重建算法，該方法結(jié)合了圖像稀疏性和圖像局部光滑性的先驗知識，在非線性共軛梯度下降算法的基礎(chǔ)上，同時自適應(yīng)地改變局部正則化參數(shù)，實驗證明該算法可以較好地恢復(fù)圖像邊緣．研究中發(fā)現(xiàn)，當采樣數(shù)逐漸減少時，上述算法重構(gòu)效果都不是很理想，重構(gòu)的圖像也不穩(wěn)定．針對這一問題，徐宗本院士[7]等提出嘗試引用L1/2范數(shù)作為新的正則項求解，通過仿真實驗表明，所提出的基于變量分裂的圖像重構(gòu)新模型提高了待重構(gòu)圖像的穩(wěn)定性，在采樣數(shù)減少的情況下也能很好的重構(gòu)出圖像．

1 MRI圖像重構(gòu)模型及解法

1.1 圖像重構(gòu)的基本模型

磁共振成像技術(shù)是一種對人體無損傷的體外影像檢查手段，可以獲得組織和器官的解剖結(jié)構(gòu)、生理功能和病變信息．磁共振掃描的速度是磁共振成像在許多應(yīng)用中的一大瓶頸，緩慢的掃描速度導(dǎo)致運動偽影的產(chǎn)生，也給動態(tài)成像帶來了困難．縮短原始數(shù)據(jù)采集時間，減少采樣數(shù)據(jù)，不僅可以提高效率、減少病患的不適，還是實現(xiàn)心血管檢查、功能信息獲取等動態(tài)成像的關(guān)鍵所在．由于MRI信號可以通過TV模型、頻域模型和字典模型等方法進行稀疏變換將信號轉(zhuǎn)換為稀疏信號，故MRI圖像本身具有稀疏性，因此MRI圖像重構(gòu)構(gòu)成壓縮傳感理論應(yīng)用的重要分支．利用壓縮傳感理論，可以極大減少對MRI圖像Fourier域的采樣數(shù)目．研究發(fā)現(xiàn)，目前的磁共振圖像重建算法大多是基于的凸優(yōu)化問題的求解，雖然基于的凸優(yōu)化算法是極其有效的，但它們進行重建所需的觀測值與理論最小值的差距甚大，特別是在采樣數(shù)少的情況下，不能高效穩(wěn)定的重構(gòu)出原始圖像，針對這一問題需要新的改進方法．

圖像重構(gòu)屬于圖像處理范疇中的逆問題，其模型表示如下

但是L0問題為NP組合難問題，對較大規(guī)模數(shù)據(jù)無法直接求解．針對這一問題，研究者們提出了一系列尋找次優(yōu)解的貪婪算法：匹配追蹤算法，正交匹配追蹤算法[9]等．但是貪婪算法時間代價過高，無法保證收斂到全局最優(yōu)．為使求解稀疏性問題變得可行，需要尋找合適的范數(shù)近似求解上述優(yōu)化問題，當前一般性的作法是將L0范數(shù)最小化問題松弛到L1范數(shù)最小化問題[10]，從而將一個組合優(yōu)化問題松弛到一個凸優(yōu)化問題來求解．即

1.2 基于交替增廣拉格朗日乘子法的算法

對于式(5)求解已經(jīng)存在很多算法，其中基追蹤BP[11]算法是最早提出的一種算方法，該算法在每次迭代時尋求最匹配的原子，文獻[12]提出的凸集投影POCS算法在兩個凸平面上進行交替投影來實現(xiàn)圖像的重構(gòu)．文獻[13]提出一種簡單的迭代收縮算法IST，只要確定步長值和閾值就可求解，但該算法收斂速度慢．文獻[14，15]提出一種兩步迭代收縮算法TwIST，該算法在每次迭代時，不像IST算法利用前一次的估計值，而是尋找新的估計值，新的估計值取決于前兩次的估計值，但收斂速度也不夠快．快速閾值收縮算法[16]FIST是對基于梯度的光滑凸優(yōu)化問題的Nesterov優(yōu)化算法的非光滑變形，雖然在時間上有所提高，但是還不能到達理想效果，文獻[17]提出一種針對稀疏結(jié)構(gòu)的稀疏重構(gòu)算法SpaRSA，實驗表明該算法的重構(gòu)效果明顯高于IST算法，但運行時間有待提高．

增廣拉格朗日乘子法[18]（ALM）相對與其他算法，重構(gòu)時間

快而且可以達到更高的精度，ALM的優(yōu)化模型如下

表1 ALM算法Tab.1ALM algorithm

ALM算法的流程圖如表1所示．

迭代過程中，當目標函數(shù)的相對變化低于停止閾值時迭代停止，默認的停止閾值為0.01．

與ALM算法相比，在交替最小化方法[19]的基礎(chǔ)上提出的交替增廣拉格朗日乘子法（ADMM）[20]，其最大的優(yōu)越性在于將原問題分解為若干個更容易得到全局解的交替的極小化子問題，充分利用了目標函數(shù)的可分離性．利用ADMM方法，Afonso[21]等提出了一種基于變量分裂法的圖像復(fù)原算法（SALSA算法）

2 改進的MRI稀疏磁共振圖像重構(gòu)方法

2.1 L1/2正則化

研究中發(fā)現(xiàn)：L1求解框架不能保證獲得滿意的稀疏解，它往往與真實稀疏解（L0的解）差距甚大[7]．并且當采樣數(shù)逐漸減少時，L1求解框架重構(gòu)效果不理想，對于含有重尾分布的誤差數(shù)據(jù)往往不能取得較好的效果，重構(gòu)的圖像也不穩(wěn)定．因此一個自然改進方法是使用Lq框架（0＜q＜1），但是Lq是非凸問題，理論分析與求解都非常困難．對于這一問題，徐宗本[6]等人根據(jù)空間理論得到啟示．

1/p+1/q=1在對偶框架下，LP1p形成一個以L2為中心的局部凸空間體系．在p+q=1的對偶框架下，Lq0q1會形成一個以為中心的非局部凸空間體系，L1/2應(yīng)該有特殊性．從圖1單位球的形狀可以看出，L1/2正則子所產(chǎn)生的解會比L1正則子產(chǎn)生的解更具有稀疏性．

圖1 不同正則子產(chǎn)生稀疏解得示意圖Fig.1Different from L1/2regularizationand L1regularizations

根據(jù)徐宗本院士的基本思想，研究一種改進的圖像重構(gòu)模型，懲罰項中將L1換成L1/2，即

2.2 MRI稀疏磁共振圖像重構(gòu)的改進模型

根據(jù)推導(dǎo)過程，求解流程如表2所示．

本文的重構(gòu)算法中將停止閾值設(shè)定為5×106．

如果式(16)和式(17)有準確解，根據(jù)Eckstein-Bertsekas定理[22]可保證此算法的收斂性．

定理1Eckstein-Bertsekas定理

式(8)的非約束優(yōu)化形式為

表2 基于正則項的磁共振圖像重構(gòu)算法Tab.2Magnetic resonance image reconstruction algorithm based on1/2regularization

表2 基于正則項的磁共振圖像重構(gòu)算法Tab.2Magnetic resonance image reconstruction algorithm based on1/2regularization

No算法1 set k=0，choose＞0,v0,d02 repeat 3xk+1=argminx1/2‖Ax y‖22+/2‖x vkdk‖22 4 dk+1=dkHzk+1b 5vk+1=argminv‖v‖1/21/2+/2‖xk+1v dk‖22 6 kk+1 7 until stopping criterion is satisfied

3 實驗結(jié)果與分析

3.1 評價標準

為驗證本文方法的有效性，采用Shepp-Logan腦部模型的磁共振圖像進行驗證，在實驗中采樣模式有好多種，例如非笛卡爾采樣、徑向采樣、螺旋采樣、變密度采樣等，但非笛卡爾采樣具有采樣速度快、對流動不敏感等優(yōu)點，它在腦功能成像、心臟冠狀動脈成像等方面得到了人們的關(guān)注，故本文采用非笛卡爾輻射狀采樣，在腦部模型的Fourier頻譜表示圖上，均勻取L條射線，然后在每條射線上高斯采樣．并且根據(jù)CSMRI理論要求可知，對k空間進行的欠采樣要滿足是隨機的，即采樣模式要滿足變換點分布函數(shù)不相關(guān)性，除了本文采用的非笛卡爾采樣外，其他幾種采樣模式通常都滿足不相關(guān)欠采樣的性質(zhì)，因此，本方法在其他采樣模式下依然適用．現(xiàn)利用CPU為2.0GHz，內(nèi)存為2G的計算機，通過MATLAB進行編碼．本文分別通過以下幾個參數(shù)對改進的模型與L1范數(shù)的進行分析比較．

式中：g為原始圖像；f為重構(gòu)圖像；d為含噪圖像；m和n為圖像的維數(shù)；MSE代表差錯率．從上式可以看出重構(gòu)出的圖像與原始圖像越接近，則‖fg‖F(xiàn)就越接近零，那么差錯率就越?。甋NR與PSNR分別代表信噪比和峰值信噪比，SNR與PSNR越大說明重構(gòu)出來的圖像越逼近原始圖像．

3.2 不同模型的比較

將正則項分別是范數(shù)L1和L1/2范數(shù)的實驗結(jié)果進行比較，L1范數(shù)的求解方法采用軟閾值方法，實驗結(jié)果如表3所示，運用兩種方法，在不用分辨率與采樣數(shù)的情況下，計算出原始圖像和重構(gòu)圖像的相對誤差，以及算法的迭代數(shù)和運行時間．從仿真結(jié)果表3中可以看出，當分辨率較低時，雖然改進模型在重構(gòu)時間上較L1模型慢，但是準確度較高，隨著分辨率的增加改進模型的重構(gòu)速度明顯高于原來的L1模型，當分辨率為256×256，采樣數(shù)為52時，改進模型的重構(gòu)時間是55.14s，L1模型的重構(gòu)時間是184.97 s，改進模型的重構(gòu)速度比原來L1模型的速度高出3倍．在采樣數(shù)減少的情況下，改進的模型的準確度明顯高于L1的模型，其相對誤差至少高出一個百分點，改進的模型在采樣射線為22時重構(gòu)出的圖像的差錯率是9.02×10-7，L1的模型在采樣射線為52時的差錯率是2.38×10-6，兩種模型的差錯率相近．

圖2給出在分辨率為32×32情況下，當采樣數(shù)不同時改進模型和L1模型差錯率的對比圖，從圖中可以看出采樣數(shù)越少，改進的模型重構(gòu)圖像的效果較于L1模型重構(gòu)的效果越好，當采樣射線減少到22時，改進的模型重構(gòu)圖像的準確度比L1模型高出大約4個百分點，與采樣數(shù)為52時重構(gòu)圖像的準確度相近．圖3給出在分辨率為128×128情況下兩種方法差錯率的對比圖，與圖2進行比較可以看出，當分辨率增加時，L1模型重構(gòu)的準確度有了一定程度的下降，改進模型基本沒有太大的波動，說明改進模型重構(gòu)圖像的穩(wěn)定性要高于原來的L1模型．

表3 不同正則化方法隨不同采樣數(shù)的圖像恢復(fù)結(jié)果Tab.3Image restoration results from different regularizations with different sampling number

圖3 兩種方法在恢復(fù)差錯率上的對比（128×128）Fig.3 Comparisonoftwomethodsintherecoveryerrorrate(128×128)

圖4是兩種算法恢復(fù)效果的PSNR和SNR，從圖中可以看出改進的模型有很大的優(yōu)勢，在采樣數(shù)為22時，改進模型的SNR是48 bB，L1模型是23 dB，比L1模型高出25 dB，改進模型的PSNR是119 dB，L1模型是100 dB，比L1模型高出19 dB．圖5給出重構(gòu)時間的對比圖，可以看出改進模型在重構(gòu)時間上有了很大的提高，尤其是采樣數(shù)少的情況下，改進模型的重構(gòu)時間是16.14 s，L1模型重構(gòu)時間是164.09s．圖5和圖6分別是MRIscan圖像的原始圖像和分辨率為64×64，采樣數(shù)為52的Fourier頻譜圖．圖7和圖8分別是L1模型和改進模型在采樣數(shù)為52的情況下的重構(gòu)結(jié)果，可以看出L1模型只是重構(gòu)出圖像的大致輪廓，差錯率為3.48×10-5，而改進模型重構(gòu)的差錯率是2.72×10-7，重構(gòu)效果高于原來的L1模型2個百分點．

圖4 兩種方法在PSNR與SNR上的對比（128×128）Fig.4Comparison of two methods in PSNR and SNR（128×128）

圖5 兩種方法在運行時間上的比較（128×128）Fig.5Comparison of two methods in PSNR and SNR(128×128)

圖6 原始圖像Fig.6The original image

圖7 采樣數(shù)為52的Fourier頻譜圖Fig.7The Fourier spectrum with the sampling number 52

圖8 模型重構(gòu)結(jié)果Fig.8Model reconstruction result

圖9 改進模型重構(gòu)結(jié)果Fig.9Improved model reconstruction result

4 結(jié)論

本文研究了一種磁共振圖像重構(gòu)的新方法，可以穩(wěn)定有效地從稀疏采樣中重構(gòu)原圖像，該模型選取L1/2范數(shù)作為正則項，解決了L1范數(shù)作為正則項時重構(gòu)圖像不穩(wěn)定和稀疏度較差的問題，為磁共振圖像重構(gòu)方法提供了新的思路．由實驗結(jié)果可以看出，該模型可以有效對稀疏磁共振圖像進行重構(gòu)，重構(gòu)效果明顯高于原來的模型，恢復(fù)時間以及信噪比都有了很大的提高，并且差錯率有了明顯的下降．數(shù)據(jù)實驗表明在采樣數(shù)減少的情況下，本文的方法可以有效重構(gòu)出圖像，增強了圖像重構(gòu)的穩(wěn)定性，提高了重構(gòu)圖像的信噪比并且降低了重構(gòu)圖像的差錯率，達到高信噪比低差錯率的理想效果．

[1]BLUMENSATHT.Compressedsensingwithnonlinearobservationsandrelatednonlinear optimizationproblem[J]．IEEEtransactionsonInformation theory[J]．2013，59（6）：3466．

[2]PENGX，ZHANGM，ZHANGJ，etal．AnalternativerecoveryalgorithmbasedonSL0formultibandsignal[C]//Instrumenta-tionandMeasurement Technology Conference（I2MTC），2013 IEEE International．IEEE，2013：114-117.

[3]BHASKAR D RAO，KENNETH KREUTZ-DELGADO．An affine scaling methodology for best basis selection[J]．IEEE Transactions on signal processing，1999，47（1）：1.

[4]BIOUCAS-DIAS J M，F(xiàn)IGUEIREDO M A T，Oliveira J P．Total variation based image deconvolution:a majorization minimization approach C// Toulouse：IEEE Conference，2006，2：278-281.

[5]李紅，楊曉梅，李青.基于加權(quán)壓縮感知的MR圖像重建方法[J].中國組織工程研究與臨床康復(fù)，2012，16（26）：4817-4821．

[6]李青，楊曉梅，李紅．基于壓縮感知的自適應(yīng)正則化磁共振圖像重建[J].計算機應(yīng)用，2012，32（2）：541-544.

[7]XU ZONGBE，CHANG XIANGYU，XU FENGMINe.L-1/2 Regularization：A thresholding representation theory and a fast solver[J]．IEEE Transactions on neural net-works and learning systems，2012，23（7）：1013-1027.

[8]CANDES EJ，ROMBERR J．Robust UncertaintyPrinciples：Exact signal reconstructionfrom highly iIncomplete frequencyinformation[J]．IEEE Transaction on Information Theory，2006，52（2）：489-509.

[9]TROPP J A，GILBERT A C．Signal recovery from random measurements via orthogonal matching pursuit[J]．IEEE Transactions on Information Theory，2007，53（12）：465-466.

[10]CANDES E J，WAKIN M.An Introduction to compressive sampling[J]．IEEE Signal Processing Magazine，2008，25（2）：21-30．

[11]HE W，QU T．Audio lossless coding/decoding method using basis pursuit algorithm[C]//Acoustics，Speech and Signal Processing（ICASSP），2013 IEEE International Conference on．IEEE，2013：552-555.

[12]PIRES R G，PEREIRA L A M，MANSANO A F，et al．A hybrid image restoration algorithm based on projections onto convex sets and harmony search[C]//Circuits and Systems（ISCAS），2013 IEEE International Symposium on．IEEE，2013：2824-2827.

[13]MA Z．Sparse principal component analysis and iterative thresholding[J]．The Annals of Statistics，2013，41（2）：772-801.

[14]BIOUCAS-DIASJM，MARIOATFIUEIREDO．ANewTwIST：Two-stepiterativeshrinkage/thresholdingalgorithmsforimagerestoration[J]．IEEE Transactions on Image processi processing，2007，16（12）：2992-3994.

[15]BIOUCAS-DIAS J M，MARIO A T FIUEIREDO．Two-step algorithms for linear inverse problem with non-quadratic regularization proceeding of ICIP[J]．IEEE International Conference on Image Processing，2007（1）：I-105-I-108.

[16]MENG X，DUAN S F，MA S X．An adaptive regularized fast iterative shrinkage-thresholding algorithm for image reconstruction in compressed sensing[J]．Advanced Materials Research，2013，710：593-597.

[17]WRIGHTS，NOWAKR，F(xiàn)IUEIREDOM．Sparsereconstructionbyseparableapproximation[J]．IEEETransactiononSignalProcessing，2009，57：2479-2493．

[18]LIN Z C，CHEN M M，MA Y．The augmented lagrange multiplier method for exact recovery of a corrupted low-rank matrix[EB/OL]．http：// arxiv.org/abs/1009.5055．2013．

[19]JAIN P，NETRAPALLI P，SANGHAVI S．Low-rank matrix completionusing alternating minimization[C]//Proceedings of theforty-fifth annual ACM symposium on Theory of computing．ACM，2013：665-674．

[20]STEPHEN BOYD，NEAL PARIKH．Distributed optimization and statistical learning via the alternating direction method of multipliers[R]．Foundation and Trends in Machine Learning，2010，3（1）：1-122．

[21]AFONSO M，BIOUCAS-DIAS J，F(xiàn)IUEIREDO M．An augmented lagrangian approach to the constrained optimization formulation of imaging inverse problems[J]．IEEE Transactions on Image Processing，2011，20（3）：681-695．

[22]AFONSOM，BIOUCAS-DIAS J，F(xiàn)IUEIREDOM．Fastimagerecoveryusingvariablesplittingandconstrainedoptimization[J]．IEEE Transactions on Image Processing，2010，19（9）：2345-2356．

[責任編輯 代俊秋]

Sparse MRI reconstruction based on L1/2regularization

MA Jie，GE Lingling，YUAN Huanchao，ZHANG Tingting

(School of Electronic Information Engineering,Hebei University of Technology,Tianjin 300401,China)

Magnetic resonance images can be reconstructed close to the original in which we use the compressed sensing from a very limited number of observation data set.However,to achieve this,we must solve a problem:the nonsmooth functions that definition in the data set to minimize.Usually the norm can produce a sparse solution,but it is different from the true sparse solution(solution).To solve the problem,a method is proposed based on variable splitting image reconstruction model.The method is introduced the norm as a new regular of the reconstructed images,which was solved by alternating theaugmentedLagrange multipliermethod.Finally,to verifythe stability andreconstructioneffect,selecting different parameters and evaluation standard are compared with existing models ofimage reconstruction.The experiments indicate that the improved model has excellent stability and gets a higher signal-to-noise ratio.

magnetic resonance imaging;compressed sensing;L1norm;L1/2norm;variable splitting

TP391.9

A

1007-2373(2015)04-0001-08

10.14081/j.cnki.hgdxb.2015.04.001

2014-12-22

國家自然科學基金（61203245）；河北省自然科學基金（F2012202027）

馬杰（1978-），男（回族），副教授．

數(shù)字出版日期：2015-06-30

數(shù)字出版網(wǎng)址：http://www.cnki.net/kcms/detail/13.1208.T.20150630.1615.001.html