基于壓縮感知的核磁共振成像重構(gòu)算法

袁 靜

（1.宿遷學院計算機科學系，江蘇 宿遷223800；2.南京郵電大學信號處理與傳輸研究院，南京210003）

基于壓縮感知的核磁共振成像重構(gòu)算法

袁 靜1，2

（1.宿遷學院計算機科學系，江蘇 宿遷223800；2.南京郵電大學信號處理與傳輸研究院，南京210003）

在核磁共振成像的應用中，一般采用聯(lián)合方式求解L1范數(shù)算子和全變差分算子，而聯(lián)合正則算子的求解模型比較復雜，為此，利用算子分裂技術(shù)求解聯(lián)合正則算子，以降低求解模型的復雜度。在此基礎(chǔ)上，提出一種迭代加權(quán)的壓縮感知核磁共振重構(gòu)算法，根據(jù)圖像在離散傅里葉變換下系數(shù)的先驗統(tǒng)計特性優(yōu)化觀測矩陣。仿真結(jié)果表明，該重構(gòu)算法不僅提高了算法的重構(gòu)精度而且減少了重構(gòu)時間。

壓縮感知；迭代加權(quán)；核磁共振成像；全變差分L1壓縮；算子分裂

DO I：10.3969/j.issn.1000-3428.2015.10.051

1 概述

近年來，Donoho［1］和Candes［2］共同提出了一種新的信號獲取理論——壓縮感知（Com pressed Sensing，CS）。壓縮感知理論的創(chuàng)新之處在于對信號進行采樣的同時對數(shù)據(jù)進行壓縮，其采樣速率遠低于奈奎斯特頻率，使得對于海量數(shù)據(jù)的信號采集變?yōu)榭赡埽?-4］。因此，國內(nèi)外學者從各個方面對這一理論提出了改進和應用［5-7］。

隨著人類醫(yī)療水平的不斷提高，核磁共振成像（Magnetic Resonance Imaging，MRI）已經(jīng)逐漸成為一種非常重要的醫(yī)療輔助技術(shù)，而如何加快核磁共振的成像速度一直是該方向的研究熱點問題。MRI自身含有壓縮感知理論成功應用的2個最為關(guān)鍵的要求［8］：（1）醫(yī)學圖像在某個適當?shù)淖儞Q域（比如小波變換等）上是可壓縮的；（2）MRI的數(shù)據(jù)直接在傅里葉域上采樣，而不是直接在像素域上［8］。將壓縮傳感理論應用到核磁共振成像中，可以極大地減少MRI成像所需要的采樣數(shù)目，從而縮短掃描時間，提高成像效率。文獻［9］提出了迭代加權(quán)L1范數(shù)法。該算法用一個加權(quán)的L1范數(shù)取代了原來的L1范數(shù)問題。由于加權(quán)系數(shù)的緣故，幅值很小的系數(shù)在迭代中會更加趨近于0，這樣迭代加權(quán)L1

范數(shù)法可以獲得比L1范數(shù)法更加稀疏的解，求解結(jié)果更加精確。文獻［9］中同時提出了迭代加權(quán)TV模型，基本思想和迭代加權(quán)L1范數(shù)法基本一致。文中的結(jié)果表明：迭代加權(quán)TV模型的重構(gòu)效果明顯優(yōu)于一般的TV模型。文獻［10］提出了將全變差分正則算子加入到壓縮感知重構(gòu)核磁共振圖像的求解模型中，取得了更好的重構(gòu)效果。文獻［11］利用算子分裂思想［12］求解聯(lián)合正則算子的問題，提出了一種全變差分L1壓縮核磁共振（Total Variation L1 Compressed MR Imaging，TVCMRI）算法。

本文引入迭代加權(quán)思想對TVCMRI思想進行改進，提出一種重構(gòu)效果更佳的算法。

2 基于迭代加權(quán)的改進求解模型

2.1 求解模型

對于一幅大小為n維的核磁共振圖像u，可以通過小波變換Φ獲得它的一個稀疏表示：

根據(jù)壓縮感知的理論，首先通過一個 m×n（m＜n）維的觀測矩陣R對其進行觀測可得：

再考慮式（1）的關(guān)系，可以得到：

其中，A=RΦ-1。

已知信號的觀測值 b以及觀測矩陣 R，就可以利用壓縮感知的重構(gòu)算法求解式（4）的方程來獲得圖像在小波變換下的稀疏表示χ：

然后，通過對 χ進行小波逆變換就可以重構(gòu)出原圖像，然而在實際應用中，觀測值b中總是會含有觀測噪聲，所以，要對式（4）進行微小的改動，轉(zhuǎn)為求解：

或者式（5）的拉格朗日形式：

其中，σ和μ為參數(shù)。

由于核磁共振圖像具有分段平滑的性質(zhì)，也就是說不同的器官其本身應該有一致的結(jié)構(gòu)，因此核磁共振圖像應該有很小的全變差分常數(shù)。圖像的全變差分定義的離散化形式為：

其中，▽1和▽2分別為第一維和第二維上的前向差分算子。

文獻［8］提出了將全變差分正則算子加入到壓縮感知重構(gòu)核磁共振圖像的求解模型中，取得了更好的重構(gòu)效果。聯(lián)合 L1范數(shù)和全變差分正則算子的求解模型為：

其中，α和β為2個正的參數(shù)。

考慮到基于迭代加權(quán)思想的改進 L1范數(shù)法和TV模型法相比原方法都有很大的提高，本文考慮利用迭代加權(quán)思想對式（8）的求解模型進行改進：其中，W2為一個對角矩陣，其對角線元素2n為加權(quán)系數(shù)。

由于全變差分正則算子和 L1范數(shù)正則算子對于χ都是非平滑的，因此式（9）所示的求解模型比單獨含有L1范數(shù)正則算子或者單獨含有全變差分正則算子的求解模型要復雜得多?？紤]利用算子分裂技術(shù)來求解聯(lián)合正則算子的問題。

2.2 算法推導

在算法推導之前，首先集中介紹推導中用到的幾個重要變量，其他變量將在推導中第一次出現(xiàn)時給予解釋。變量u∈Rn1×n2表示一幅二維的像素大小為n1×n2的核磁共振圖像。算子 L=（▽1，▽2）：Rn1×n2→Rn1×n2×Rn1×n2表示沿著第一維坐標和第二維坐標的有限差分算子，它的子算子并且令Ψ=Φ-1，對于任意的正交變換Φ，Ψ=Φ*為Φ的伴隨算子。其中，A=

RΨ。根據(jù)以上定義的符號，式（9）可以寫為：

另外，在推導中還要反復用到泛函的一個性質(zhì)，對于一個凸函數(shù)f（χ）存在以下性質(zhì)：

下面給出利用算子分裂理論求解式（10）的詳細推導過程：

如果算法存在最優(yōu)解χ，那么必然滿足：

其中，?χ表示對元素χ求導，▽χh（χ）=A*（Aχ-b）= ΦR*（RΨχ-b）。

令Lij（Ψχ）=z，則式（12）變?yōu)椋?/p>

由式（13）可以得出下面2個關(guān)系：

根據(jù)式（11）的關(guān)系，式（15）可以化為：

利用算子分裂理論，可以分別對式（14）和式（16）進行處理：

其中，τ1和τ2為2個大于0的參數(shù)。

這樣，由式（17）～式（20）已經(jīng)可以比較清楚地看出算法的基本結(jié)構(gòu)了。給定χi和yij，可以分別根據(jù)式（18）和式（20）求出si和tij；相反，若已知si和同樣可以分別根據(jù)式（17）和式（19）算出 χi和 yij。下面推導式（17）和式（19）顯示的表達式。

由式（17），有：

對于絕對值的求導，令：

對于式（19），有：

根據(jù)式（11）的關(guān)系，有：

移項后，有：

對于2-范數(shù)的求導，令：

對式（25）進行推導可得：

本文所提出的算法詳細步驟為：

輸入 核磁共振成像的部分采樣值b，算子L和它的伴隨算子 L*，部分傅里葉觀測矩陣 R，正交小波變換算子Ψ及其逆變換Φ，參數(shù)α和 β的衰減速率ηα和 ηβ，以及它們的終值α-和β-，最大迭代次數(shù)Maχ，算法終止標準ftol，加權(quán)系數(shù)中的調(diào)整參數(shù)ε1和 ε2，步長 τ1和 τ2

輸出 重構(gòu)圖像u=Ψχ

2.3 觀測矩陣的優(yōu)化

在基于壓縮感知的核磁共振成像中，觀測矩陣R為部分離散傅里葉變換矩陣，可以通過隨機選擇n×n的離散傅里葉變換矩陣的 m（m＜n）行作為觀測矩陣R。文獻［10］指出，由于核磁共振成像的傅里葉變換的低頻能量較高而高頻能量較低，因此在

低頻選擇較多的采樣點顯得更為有效，并且在文中提出了變密度的采樣矩陣。文獻［11］中由于離散傅里葉變換的對稱性提出了將變換矩陣的上半部分全部（除第一行右半部分）不采樣，在下半部分采用變密度采樣：越靠近右下和左下的部分采樣越多，越靠近中心的部分采樣越少。經(jīng)過比較，發(fā)現(xiàn)文獻［11］提出的這種觀測矩陣的效果要優(yōu)于文獻［10］提出的觀測矩陣。

文獻［13］提出基于圖像在不同變換下的變換系數(shù)能量的統(tǒng)計特征的先驗知識來優(yōu)化觀測矩陣，取得了更好的重構(gòu)效果。本文受文獻［13］的啟發(fā)，結(jié)合文獻［11］中的采樣矩陣的結(jié)構(gòu)，提出一種改進的觀測矩陣。

對于一幅大小為M×N的圖像，首先根據(jù)圖像在離散傅里葉變換下系數(shù)能量的統(tǒng)計特征的先驗知識，確定圖像中每個點被采樣的概率：

圖1 采樣率為30%時的觀測矩陣

本文所提出的觀測矩陣確定后，可以保存連續(xù)使用。由于利用了圖片在離散傅里葉變換下系數(shù)的先驗統(tǒng)計特性，因此更能夠采樣到對重構(gòu)貢獻大的系數(shù)，獲得較好的重構(gòu)效果。

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

本節(jié)通過仿真實驗證明算法的有效性，所有的程序均在M atlab環(huán)境下運行。實驗中，設(shè)置程序中參數(shù)α和β的終值為α-=1×10-3，β-=3.5×10-2，它們的衰減速度分別為ηα=0.25和ηβ=0.25。程序最大的迭代次數(shù)Maχ=20，程序終止條件ftol=1× 10-3，程序步長τ1=τ2=0.8，加權(quán)系數(shù)中的調(diào)整系數(shù)分別為ε1=0.2，ε2=20。觀測中混入的高斯白噪聲方差為σ2=1×10-4。

利用仿真程序測試了4種不同人體器官的核磁共振圖像：210×210的腦部圖像，220×220的胸部圖像，220×220的腎臟圖像以及208×924的人體全身圖像。

圖2顯示采樣率為 20%時胸部的重構(gòu)圖像，圖2（a）為原始圖像，圖2（b）為TVCMRI算法的重構(gòu)圖像，圖2（c）為利用本文提出的改進觀測矩陣TVCMRI算法的重構(gòu)圖像，圖2（d）為本文提出算法的重構(gòu)圖像。可以看出，TVCMRI算法作為一種求解聯(lián)合L1范數(shù)和TV模型的有效算法，在采樣率為20%的情況下可以很好地對原圖像進行重構(gòu)。利用本文提出的優(yōu)化的觀測矩陣，TVCMRI算法（在后面的圖中用TCM 1表示該算法）的重構(gòu)圖像更加明亮，表明優(yōu)化的觀測矩陣的采樣能量更為豐富。而本文算法的重構(gòu)效果顯然更好，人體器官的一些復雜的細節(jié)也得到了較好的重構(gòu)。

圖2 采樣率為20%時的重構(gòu)圖像

表1記錄了仿真實驗中3種算法在不同采樣率下不同部位圖像重構(gòu)的峰值性噪比和重構(gòu)時間。圖3顯示了不同采樣率下3種算法重構(gòu)圖像的峰值信噪比的平均值。從圖中可以看出，利用本文提出的優(yōu)化的采樣矩陣，TVCMRI重構(gòu)圖像的峰值信噪比可以獲得2 dB左右的提高。由于引入了迭代加權(quán)思想，本文提出的改進算法較之原來的TVCMRI算法在采樣率相同并且都使用優(yōu)化的觀測矩陣的情況下，峰值信噪比可以獲得3 dB～5 dB的提高。從圖中可以發(fā)現(xiàn)，利用本文提出的改進算法和優(yōu)化的

觀測矩陣，相比于原來的TVCMRI算法和觀測矩陣，在相同的重構(gòu)效果要求下可以節(jié)省10%的觀測數(shù)據(jù)量。

表1 圖像的峰值信噪比和重構(gòu)時間比較

圖3 不同采樣率下算法重構(gòu)圖像的峰值信噪比

從表1中也可以看到，本文所提出的改進算法的計算時間大概為原TVCMRI算法的2倍。但是由于該算法在相同的重構(gòu)要求下可以節(jié)省10%的采樣數(shù)據(jù)，也就是可以大幅減少患者的核磁共振掃描時間，相比后端重構(gòu)所多出的這點時間可以忽略不計。

4 結(jié)束語

壓縮感知技術(shù)在核磁共振成像中有著廣闊的應用前景。本文利用算子分裂技術(shù)和迭代加權(quán)思想改進了求解模型，提出了迭代加權(quán)的聯(lián)合L1和全差變分正則算子模型，并且推導出了基于算子分裂技術(shù)的重構(gòu)算法，同時優(yōu)化了觀測矩陣。仿真實驗證明了改進算法的有效性，但從仿真實驗中也可以看出改進算法的計算時間有所延長，今后將進一步利用平滑技術(shù)對重構(gòu)算法進行優(yōu)化，以降低算法重構(gòu)的運行時間。

［1］ Donoho D L.Compressed Sensing［J］.IEEE Transactions on Information Theory，2006，52（4）：1289-1306.

［2］ Candes E.Compressive Sampling［C］//Proceedings of International Congress of Mathematicians，Madrid，Spain：［s.n.］，2006：1433-1452.

［3］ 李少東，楊 軍，馬曉巖.基于壓縮感知的ISAR高分辨成像算法［J］.通信學報，2013，34（9）：150-157.

［4］ 寧方立，何碧靜，韋 娟.基于lp范數(shù)的壓縮感知圖像重建算法研究［J］.物理學報，2013，62（17）：174-212.

［5］ 馮 振，郭 禾，王宇新，等.CS-MRI中稀疏信號支撐集混合檢測方法［J］.計算機工程，2014，40（5）：164-167.

［6］ 史久根，吳文婷，劉 勝.基于壓縮感知的圖像重構(gòu)算法［J］.計算機工程，2014，40（2）：229-232.

［7］ 閆 鵬，王阿川.基于壓縮感知的CoSaMP算法自適應性改進［J］.計算機工程，2013，39（6）：28-33.

［8］ Lustig M，Donoho D L，Santos J M，et al.Com pressed Sensing MRI［J］.IEEE Signal Processing Magazine，2008，25（3）：72-82.

［9］ Candes E J，Wakin M B，Body S P.Enhancing Sparsity by Reweighted L1 Minimization［J］.Journal of Fourier Analysis and Applications，2008，14（5）：877-905.

［10］ Lustig M，Donoho D L，Pauly J M.Sparse MRI：The Application of Com pressed Sensing for Rapid MR Imaging［J］.Magnetic Resonance Medicine，2007，58（6）：1182-1195.

［11］ Ma Shiqian，Yin Wotao，Zhang Yin，et al.An Efficient Algorithm for Com pressed MR Imaging Using Tatal Variation and Wavelets［C］//Proceedings of IEEE Conference on Computer Vision and Pattern Recognition.Washington D.C.，USA：IEEE Press，2008：1-8.

［12］ Lions P.L，Mercier B.Splitting Algorithms for the Sun of Two Nonlinear Operators［J］.SIAM Journal on Numerical Analysis，1979，16（3）：964-979.

［13］ Wang Zhongmin，Arce G R.Variable Density Compressed Sampling［J］.IEEE Transactions on Image Processing，2010，19（1）：264-270.

編輯 顧逸斐

Magnetic Resonance Imaging Reconstruction Algorithm Based on Compressed Sensing

YUAN Jing1，2
（1.Department of Computer Science，Suqian College，Suqian 223800，China；2.Institute of Signal Processing and Transmission，Nanjing University of Posts and Telecommunications，Nanjing 210003，China）

In the application of Magnetic Resonance Imaging（MRI），it is common to solve the problem by combining L1 norm with total variation operator.Because the model of solving compound regularizer is more complicated，the operator splitting technique is used to solve the problem of compound regularizer，which is in order to lower the complexity of the solving model，and puts forward a reconstruction method which is iterative weighted.The observation matrix is optimized，according to the priori statistical properties of imaging，which is under different transformations. Simulation results show that this image reconstruction algorithm not only enhances the reconstruction accuracy，but also decreases the time for the reconstruction.

Compressed Sensing（CS）；iterative weighted；Magnetic Resonance Imaging（MRI）；total variation L1 compression；operator splitting

袁 靜.基于壓縮感知的核磁共振成像重構(gòu)算法［J］.計算機工程，2015，41（10）：270-274.

英文引用格式：Yuan Jing.Magnetic Resonance Imaging Reconstruction Algorithm Based on Compressed Sensing［J］. Computer Engineering，2015，41（10）：270-274.

1000-3428（2015）10-0270-05

A

TP301.6

國家自然科學基金資助項目（61372122）；宿遷市科技創(chuàng)新專項基金資助項目（Z201209）。

袁 靜（1979-），女，講師、碩士，主研方向：信號處理。

2014-09-17

2014-10-31E-mail：yuanxiaojing1979@163.com