方晟 吳文川 應(yīng)葵 郭華?

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

展開可得

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

即

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

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

對(duì)于第k次迭代,有

5 實(shí)驗(yàn)系統(tǒng)

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

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

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

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

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

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

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

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

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

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

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

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

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

7 結(jié)論

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

[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物理學(xué)報(bào)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