袁 靜

（1.宿遷學(xué)院計(jì)算機(jī)科學(xué)系，江蘇 宿遷223800；2.南京郵電大學(xué)信號(hào)處理與傳輸研究院，南京210003）

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

袁 靜1，2

（1.宿遷學(xué)院計(jì)算機(jī)科學(xué)系，江蘇 宿遷223800；2.南京郵電大學(xué)信號(hào)處理與傳輸研究院，南京210003）

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

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

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

1 概述

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

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

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

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

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

2.1 求解模型

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

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

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

其中，A=RΦ-1。

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

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

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

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

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

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

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

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

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

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

2.2 算法推導(dǎo)

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

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

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

下面給出利用算子分裂理論求解式（10）的詳細(xì)推導(dǎo)過(guò)程：

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

其中，?χ表示對(duì)元素χ求導(dǎo)，▽?duì)謍（χ）=A*（Aχ-b）= ΦR*（RΨχ-b）。

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

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

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

利用算子分裂理論，可以分別對(duì)式（14）和式（16）進(jìn)行處理：

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

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

由式（17），有：

對(duì)于絕對(duì)值的求導(dǎo)，令：

對(duì)于式（19），有：

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

移項(xiàng)后，有：

對(duì)于2-范數(shù)的求導(dǎo)，令：

對(duì)式（25）進(jìn)行推導(dǎo)可得：

本文所提出的算法詳細(xì)步驟為：

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

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

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

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

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

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

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

圖1 采樣率為30%時(shí)的觀測(cè)矩陣

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

3 實(shí)驗(yàn)結(jié)果與分析

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

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

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

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

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

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

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

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

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

4 結(jié)束語(yǔ)

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

［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］.通信學(xué)報(bào)，2013，34（9）：150-157.

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

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

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

［7］ 閆 鵬，王阿川.基于壓縮感知的CoSaMP算法自適應(yīng)性改進(jìn)［J］.計(jì)算機(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］.計(jì)算機(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

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

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

2014-09-17

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