基于壓縮感知的核磁共振成像重構算法

袁 靜

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

基于壓縮感知的核磁共振成像重構算法

袁 靜1，2

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

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

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

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

1 概述

近年來，Donoho［1］和Candes［2］共同提出了一種新的信號獲取理論——壓縮感知（Com pressed Sensing，CS）。壓縮感知理論的創新之處在于對信號進行采樣的同時對數據進行壓縮，其采樣速率遠低于奈奎斯特頻率，使得對于海量數據的信號采集變為可能［3-4］。因此，國內外學者從各個方面對這一理論提出了改進和應用［5-7］。

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

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

本文引入迭代加權思想對TVCMRI思想進行改進，提出一種重構效果更佳的算法。

2 基于迭代加權的改進求解模型

2.1 求解模型

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

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

再考慮式（1）的關系，可以得到：

其中，A=RΦ-1。

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

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

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

其中，σ和μ為參數。

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

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

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

其中，α和β為2個正的參數。

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

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

2.2 算法推導

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

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

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

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

如果算法存在最優解χ，那么必然滿足：

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

令Lij（Ψχ）=z，則式（12）變為：

由式（13）可以得出下面2個關系：

根據式（11）的關系，式（15）可以化為：

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

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

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

由式（17），有：

對于絕對值的求導，令：

對于式（19），有：

根據式（11）的關系，有：

移項后，有：

對于2-范數的求導，令：

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

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

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

輸出 重構圖像u=Ψχ

2.3 觀測矩陣的優化

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

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

文獻［13］提出基于圖像在不同變換下的變換系數能量的統計特征的先驗知識來優化觀測矩陣，取得了更好的重構效果。本文受文獻［13］的啟發，結合文獻［11］中的采樣矩陣的結構，提出一種改進的觀測矩陣。

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

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

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

3 實驗結果與分析

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

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

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

圖2 采樣率為20%時的重構圖像

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

觀測矩陣，相比于原來的TVCMRI算法和觀測矩陣，在相同的重構效果要求下可以節省10%的觀測數據量。

表1 圖像的峰值信噪比和重構時間比較

圖3 不同采樣率下算法重構圖像的峰值信噪比

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

4 結束語

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

［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范數的壓縮感知圖像重建算法研究［J］.物理學報，2013，62（17）：174-212.

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

［6］ 史久根，吳文婷，劉 勝.基于壓縮感知的圖像重構算法［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

袁 靜.基于壓縮感知的核磁共振成像重構算法［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）；宿遷市科技創新專項基金資助項目（Z201209）。

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

2014-09-17

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