方晟 吳文川 應葵 郭華?

1 引言

縮短數(shù)據(jù)采集時間一直是磁共振成像發(fā)展的核心目標之一.傳統(tǒng)快速成像方法通過提高磁場梯度性能加快k空間數(shù)據(jù)采集速度.然而,受制于生理因素和硬件制約,這類方法加速效果已經(jīng)接近極限.而近年誕生的壓縮感知成像理論,提供了通過欠采樣縮短數(shù)據(jù)采集時間的新途徑.壓縮感知成像理論指出：如果欠采樣導致混迭偽影是空間非連貫的,通過采用合適的非線性重建方法,可以在抑制混迭偽影的同時,完整重建圖像信息[1-3].

非均勻螺旋線序列便是一種符合壓縮感知成像理論關(guān)于混迭偽影的要求的數(shù)據(jù)采集方法[4].非均勻螺旋線序列在k空間中心附近進行過采樣,在k空間周圍進行正常采樣[5].這使得欠采樣情況下,非均勻螺旋線序列的混迭偽影在空間分布呈現(xiàn)出近似不連貫的分布,不會嚴重破壞圖像結(jié)構(gòu)信息[6].同時非均勻螺旋線序列本身便是一種數(shù)據(jù)采集時間短的超快速數(shù)據(jù)采集方法[5],在心臟成像、灌注成像等方面已獲得應用[6-8],并有望在磁共振彈性成像[9,10]等新領域獲得進一步的應用.所以,非均勻螺旋線序列非常適合作為包括壓縮感知成像在內(nèi)的快速磁共振成像的數(shù)據(jù)采集手段.

實現(xiàn)壓縮感知成像的另一個必備條件是能夠去除欠采樣所引起混迭偽影的非線性圖像重建算法.目前在磁共振成像領域上應用較多的非線性重建方法主要有小波域稀疏性約束重建方法[3]和全變分重建方法[11].其中全變分方法在去除噪聲和混迭偽影方面獲得了一定的成功.然而,全變分方法在抑制欠采樣混迭偽影的同時,會造成圖像細節(jié)結(jié)構(gòu)的丟失[12,13].

為了解決這一問題,本文結(jié)合非均勻螺旋線磁共振成像序列和布雷格曼迭代[14]重建,提出一種新的快速磁共振成像方法,并通過數(shù)值模擬實驗和在體磁共振成像進行了實驗驗證.實驗結(jié)果表明,布雷格曼迭代能有效消除混迭偽影和吉布斯振鈴偽影,并保持準確的圖像結(jié)構(gòu)信息,這對于建立基于非均勻螺旋線的快速磁共振成像方法具有重要推動作用.

2 非均勻螺旋線磁共振成像

2.1 非均勻螺旋線數(shù)據(jù)采集序列設計

非均勻螺旋線磁共振數(shù)據(jù)采集序列可以靈活控制k空間的數(shù)據(jù)采集密度,改變欠采樣引起混迭偽影的形態(tài),并針對性地獲取實時掃描信息.其k空間軌跡方程為[5]

其中,k(τ)是螺旋線軌跡在k空間的位置;τ是關(guān)于時間的方程,決定了螺旋線軌跡的形狀;n是螺旋線的匝數(shù);λ是由視野范圍FOV和成像矩陣大小決定的常數(shù),λ=N/(2×FOV).

與大多數(shù)核磁共振序列一樣,非均勻螺旋線序列在磁場梯度爬升過程中會受到最大磁場梯度爬升率的限制;隨后當梯度達到硬件允許的最大幅度后,限制條件變?yōu)橛布试S的最大磁場梯度幅度.針對這兩種不同的邊界條件,將(1)式中的τ(t)描述如下：

其中,γ是磁化率;sm和gm分別是最大磁場梯度爬升率和磁場幅度;Ts2α是邊界條件由最大磁場梯度爬升率轉(zhuǎn)變?yōu)樽畲蟠艌鎏荻确鹊臅r刻;Tes是邊界條件持續(xù)為最大磁場梯度爬升率的時間;Tea是非均勻螺旋線序列的結(jié)束時刻(起始時刻為t=0).

α是控制采樣密度變化的常數(shù)(α≥1).α=1對應于均勻螺旋線采樣軌跡,整個k空間采樣率等于奈奎斯特采樣率.欠采樣時,混迭偽影形態(tài)是連貫的.α＞1時,k空間中心采樣率高于奈奎斯特采樣率.α越大,k空間中心采樣密度越大,欠采樣情況下,混迭偽影空間分布連貫性越低,越接近噪聲,越有利于運用非線性重建算法去除混迭偽影.

然而,增大α會造成采集時間的增加,所以,非均勻螺旋線需要在縮短采集時間和降低混迭偽影的空間連貫性之間進行均衡取舍.通過實驗比對,本文中設計和編寫了α=3的非均勻螺旋線數(shù)據(jù)采集序列,用于在體成像實驗數(shù)據(jù)的采集.

2.2 非均勻螺旋線磁共振成像的重建

非均勻螺旋線磁共振成像可以通過以下的線性方程表示：

其中,s(k)是所采集的k空間信號,ρ(r)是待重建的磁共振圖像;k是k空間的位置向量,二維情況下,k=[kx,ky],kx=Re(k(τ)),ky=Im(k(τ));r 是空間向量r=[rx,ry].離散化后,方程(3)可以表示成

其中,s是所采集的k空間信號向量,ρ是待重建的磁共振圖像向量;E是由非均勻螺旋線采樣軌跡決定的圖像編碼采樣矩陣,具體形式為

其中N為圖像矩陣一行的像素數(shù),M為采樣數(shù)據(jù)總數(shù),km是離散化后非均勻螺旋線軌跡上的離散采樣點位置向量,rn是離散化后的空間位置向量.圖像重建過程便是對(4)式進行矩陣求逆的過程,經(jīng)典的線性方法包括非直角坐標傅里葉變換[15]和網(wǎng)格插值法[16].但在欠采樣條件下,這些線性重建方法無法消除混迭偽影.

3 壓縮感知成像框架下的非線性重建

壓縮感知成像理論提供了欠采樣條件下,利用非線性方法完整重建圖像信號的理論框架.根據(jù)壓縮感知成像理論[2],對于長度為N2的磁共振向量信號ρ,假設經(jīng)過線性變換x=Ψρ后,可以得到非零值個數(shù)為K稀疏信號x.此時,如果采樣矩陣E同時滿足約束等距性準則：

則總共只需采集M=O(Klog(N2/K))·N2個數(shù)據(jù),即可以完整重建出信號.在當前實用磁共振采樣軌跡中,非均勻螺旋線數(shù)據(jù)采集軌跡正是與約束等距性條件符合最好的數(shù)據(jù)采集方式之一[3].

壓縮感知成像的重建可以通過求解下述優(yōu)化問題來實現(xiàn)：

其中J(·)是提升解的稀疏性的非線性泛函;λ是正則化因子.對于目前應用較多的全變分重建方法而言,Ψ 為差分算子矩陣 ?,J(·)= ‖·‖1為 1-范數(shù),其表達式為

4 基于布雷格曼迭代的重建方法

布雷格曼迭代是在全變分的基礎上,通過多次迭代求解全變分泛函J(·)的布雷格曼距離來減小圖像結(jié)構(gòu)損失的問題,其表達式為

(9)式有兩個重要的性質(zhì)：首先,在(9)式的迭代過程中,‖s-Eρ‖2單調(diào)遞減,直至為 0.其次,在 (9)式的迭代求解中,只要滿足 ‖s-Eρk‖2＞‖s-Eρtrue‖2,ρ 單調(diào)地趨近真實解 ρtrue.因此,布雷格曼方法具有穩(wěn)定的收斂性;并且通過選擇合適的迭代次數(shù),即可以獲得混迭偽影得到充分抑制,并且圖像結(jié)構(gòu)信息保存較完好的求解結(jié)果.

然而,由于J(·)為非線性泛函,qk-1的計算較為復雜.為了簡化計算,我們結(jié)合(9)式的特點對(10)式進行變換.

當(9)式收斂到穩(wěn)定最優(yōu)解時,有

展開可得

同時,收斂到穩(wěn)定最優(yōu)解時,有ρk=ρk-1,代入(13)式可得：

即

根據(jù)(16)式便可以只通過矩陣-向量運算,通過迭代的方法計算出qk-1,避免了計算帶來復雜度.因此,(9)式可以通過如下迭代過程求解.

初始值：u0=0,?J(u0)=0.

對于第k次迭代,有

5 實驗系統(tǒng)

為了驗證所建立的快速成像方法,分別進行了水模成像實驗和在體成像實驗,并比較了線性重建、全變分重建和布雷格曼迭代重建三種方法的重建結(jié)果.

水模成像實驗和在體成像實驗均在Philips公司的3 T磁共振掃描儀上進行,采用所編寫的自旋回波非均勻螺旋線序列采集數(shù)據(jù),具體參數(shù)如下：α=3,翻轉(zhuǎn)角度為90°,回波時間TE=20 ms,重復時間TR=500 ms,成像視野220 mm×220 mm,圖像采集矩陣大小為256×256.信號采集前有實行高階勻場.所編寫序列的梯度時序波形如圖1(a)所示,其對應k空間軌跡如圖1(b)所示.

水模成像實驗和在體成像實驗均首先進行全采樣,驗證所設計序列的可行性,并同時以全采樣數(shù)據(jù)的重建結(jié)果作為比較評價算法的金標準;然后再進行欠采樣,僅采集1/3的數(shù)據(jù),比較不同方法的重建效果.

圖1 所編寫非均勻螺旋線序列的波形與軌跡 (a)梯度波形圖;(b)k空間軌跡

6 實驗結(jié)果分析

圖2顯示水模成像實驗結(jié)果.如圖2(a)所示,全采樣情況下,所設計的非均勻螺旋線序列的成像結(jié)果具有較高的信噪比與分辨率,通過高階勻場,偏共振引起的模糊效應也限制在最低.在欠采樣情況下(加速倍數(shù)=3),標準線性重建結(jié)果出現(xiàn)了混迭偽影.然而這些混迭偽影在空間中分布均勻,并且強度較低,對圖像的真實結(jié)構(gòu)信息破壞較少(圖2(b)),達到了通過非均勻螺旋線軌跡降低混迭偽影空間連貫性的設計目標.

圖2 水模成像實驗結(jié)果 (a)全采樣水模圖像;(b)欠采樣數(shù)據(jù)的線性方法重建結(jié)果;(c)欠采樣數(shù)據(jù)的全變分方法重建結(jié)果;(d)欠采樣數(shù)據(jù)的布雷格曼迭代重建結(jié)果

全變分重建結(jié)果雖然整體上偽影水平要低于線性重建結(jié)果,但是在水模下方的殘余混迭偽影依然存在(圖2(c)中黑色箭頭所指).同時,全變分重建方法還導致了水模中細節(jié)結(jié)構(gòu)的模糊(圖2(c)中白色箭頭所指).而布雷格曼迭代則有效減弱殘余混迭偽影的空間強度(圖2(d)中黑色箭頭所指).并且,布雷格曼重建結(jié)果的細節(jié)結(jié)構(gòu)清晰(圖2(d)中黑色箭頭所指).與其他兩種重建方法相比,布雷格曼迭代重建結(jié)果與全采樣水模圖像(圖2(a))最為符合.

圖3顯示在體成像實驗結(jié)果.如圖3(a)所示,全采樣情況下,所采集的頭部圖像結(jié)構(gòu)清晰,信噪比較高,并且沒有明顯的偏共振引起的模糊問題,具有較好的成像質(zhì)量.在欠采樣情況下(加速倍數(shù)=3),標準線性重建結(jié)果出現(xiàn)了呈振鈴狀的混迭偽影(圖3(b)).然而,由于強度較低,混迭偽影并沒有掩蓋圖像的真實結(jié)構(gòu)信息,對圖像的診斷信息破壞較少.

圖3 在體成像實驗結(jié)果 (a)全采樣圖像;(b)欠采樣數(shù)據(jù)的線性方法重建結(jié)果;(c)欠采樣數(shù)據(jù)的全變分方法重建結(jié)果;(d)欠采樣數(shù)據(jù)的布雷格曼迭代重建結(jié)果

全變分重建方法基本上去除了混迭偽影,然而,全變分重建方法同時也帶來了細節(jié)結(jié)構(gòu)處圖像分辨率的下降問題,出現(xiàn)了一定的模糊效應(圖3(c)中箭頭所指).布雷格曼迭代則在去除混迭偽影的同時,保持了圖像的分辨率,具有清晰銳利的細節(jié)結(jié)構(gòu)(圖3(d)中箭頭所指),和參考圖像(圖3(a))非常接近.

表1給出了三種重建算法的均方根誤差(以全采樣圖像為金標準).無論是水模成像實驗還是在體成像實驗,線性重建方法由于混迭偽影的影響,均方根誤差誤差最大.由于能夠有效去除偽影,全變分重建方法的均方根誤差小于線性重建方法.而布雷格曼迭代同時解決了線性重建方法的混迭偽影問題和全變分方法的細節(jié)結(jié)構(gòu)模糊問題,其均方根誤差最小.

表1 三種重建算法的均方根誤差比較

7 結(jié)論

本文基于壓縮感知成像理論,將非均勻螺旋線成像序列和布雷格曼迭代重建方法結(jié)合起來,提出了一種新的快速磁共振成像方法,并通過磁共振水模成像實驗和在體成像實驗對該方法進行了驗證.實驗結(jié)果表明：所提出快速成像方法能夠在只采集1/3數(shù)據(jù)的情況下,有效去除混迭偽影并保持精細的圖像細節(jié)信息,獲得與全采樣圖像符合良好的成像結(jié)果,從而保證了圖像質(zhì)量.該方法顯著縮短了磁共振的成像時間,可以解決磁共振成像某些時間分辨率要求高的醫(yī)學成像應用中的瓶頸問題,有望在心血管成像、腦功能成像和灌注成像等方面中獲得進一步的應用.

[1]Cand`es E J,Romberg J,Tao T 2006 IEEE Tran.Inform.Theory 52 489

[2]Donoho D L 2006 IEEE Tran.Inform.Theory 52 1289

[3]Lustig M,Donoho D,Pauly J M 2007 Magn Reson.Med.58 1182

[4]Lustig M,Donoho D L,Santos J M,Pauly J M 2008 IEEE Signal Proc.Mag.25 72

[5]Kim D,Adalsteinsson E,Spielman D M 2003 Magnet.Reson.Med.50 214

[6]Tsai C M,Nishimura D G 2000 Magnet.Reson.Med.43 452

[7]LeeJH,HargreavesBA,HuBS,NishimuraDG2003Magnet.Reson.Med.50 1276

[8]Lu G,Liu M L,Li L Y,Ye C H 2002 Chin.Phys.Lett.19 1385

[9]Bensamoun S F,Glaser K J,Ringleb S I,Chen Q,Ehman R L,An K N 2008 J.Magn.Reson.Imaging.27 1083

[10]Wang H Z,Xu L F,Yu J,Huang Q M,Wang X Y,Lu L 2010 Acta Phys.Sin.59 7463(in Chinese)[汪紅志,許凌峰,俞捷,黃清明,王曉琰,陸倫2010物理學報59 7463]

[11]Rudin L,Osher S,Fatemi E 1992 Physica D 60 259

[12]Chan T,Esedoglu S,Park F,Yip A 2006 Mathematical Models of Computer Vision：the Handbook(Boston,MA：Springer)p176

[13]Block K T,Uecker M,Frahm J 2007 Magnet.Reson.Med.57 1086

[14]Osher S,Burger M,Goldfarb D,Xu J,Yin W 2005 Multiscale Model.Sim.4 460

[15]Sha L,Guo H,Song A W 2003 J.Magnet.Reson.162 250

[16]Jackson J I,Meyer C H,Nishimura D G,Macovski A 1991 IEEE Trans.Med.Imaging.10 473