周愛珍 梅穎潔 王 聰 莫紀(jì)江 陳武凡 馮衍秋

（南方醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院，廣州 510515）

引言

并行磁共振成像（PMRI）概念的提出，突破了傳統(tǒng)MR成像時(shí)間受到磁場梯度以及射頻場硬件性能約束的限制。多線圈部分采集方式是：使用多個(gè)射頻接收線圈同時(shí)接收感應(yīng)信號(hào)，利用各個(gè)線圈的空間敏感度差異來編碼空間信息，降低成像所需的梯度編碼步數(shù)，提高數(shù)據(jù)獲取速度?，F(xiàn)有的并行磁共振成像方法主要有3類［1，3］：第一類是基于圖像域的成像方法，如敏感度編碼方法（SENSE）和基于局部化敏感度分布的并行成像方法（PILS）；第二類是基于k空間域的成像方法，運(yùn)用比較成熟的是GRAPPA方法，另外還有空間諧波并行采集技術(shù)（SMASH），以及在該算法上改進(jìn)出的AUTO-SMASH和VD-AUTO-SMASH方法；第三類是同時(shí)基于k空間域和圖像域的重建方法，如SPACE-RIP。在已知線圈敏感度的情況下，SENSE重建方法可以得到相對較好的圖像。在SENSE中，重建問題是解一系列的線性方程［2］，即

式中，E為編碼矩陣，；sγ為第γ個(gè)線圈的線圈敏感度；kκ為k空間第κ個(gè)采樣位置；rρ為第ρ個(gè)體素的位置；f為待求的全視野圖像；d為由各個(gè)線圈采集的部分k空間數(shù)據(jù)形成的矢量。

通常采用最小化殘差的L2范數(shù)的平方方法求最優(yōu)解。但是，方程的解非常依賴于線圈形狀和采樣軌跡，特別是當(dāng)欠采樣因子比較大時(shí)，方程往往是病態(tài)的，重建結(jié)果的信噪比（SNR）會(huì)降低。為了提高圖像信噪比，通常在重建方程中添加懲罰項(xiàng)。TV約束和Tikhonov約束是兩種常用的約束方式［4］，PM模型也是較好的約束方式。

SENSE重建TV約束模型為freg=

SENSE重建Tikhonov約束模型［2］為freg=，其中，ξ為正則化參數(shù)。和TV約束方式不同，L2范數(shù)懲罰項(xiàng)值的大小和作用域的整個(gè)趨勢都有關(guān)。當(dāng)圖像存在“跳變”時(shí)，L2范數(shù)值會(huì)很大，約束最小化會(huì)使“跳變”被約束，這樣能很好地去除平滑區(qū)域內(nèi)的“跳變”噪聲點(diǎn)，但也同時(shí)會(huì)導(dǎo)致圖像邊緣信息被模糊。

綜上可知，TV模型能保護(hù)邊緣信息，但會(huì)產(chǎn)生“階梯效應(yīng)”；Tikhonov約束可以在較大程度上平滑噪聲及避免“階梯效應(yīng)”，但是會(huì)導(dǎo)致圖像邊緣模糊［16］。為了獲得更好的約束方式，考慮將這兩種平滑性度量結(jié)合起來，形成某種“混合”測度。筆者提出，由先驗(yàn)圖像的梯度特征決定懲罰項(xiàng)，對重建方程約束方式進(jìn)行自適應(yīng)調(diào)整，并借鑒PM模型的思想，進(jìn)一步調(diào)整約束方式，以便在梯度值較大的邊緣區(qū)懲罰項(xiàng)使用TV約束方式，在梯度值較小的平滑區(qū)使用Tikhonov約束方式。這樣，可以在保留圖像細(xì)節(jié)的同時(shí)抑制噪聲，增加圖像信噪比。

1 自適應(yīng)約束SENSE重建理論

對采樣數(shù)據(jù)進(jìn)行網(wǎng)格化處理，然后利用平方和重建方法［18］進(jìn)行圖像重建，將得到的結(jié)果作為先驗(yàn)圖像。根據(jù)先驗(yàn)圖像的梯度特征決定懲罰項(xiàng)［14，17］，從而對多通道非笛卡爾軌跡采樣數(shù)據(jù)進(jìn)行SENSE重建，重建模型為：

式中，p（x′）=）是關(guān)于先驗(yàn)圖像m0的梯度幅度下降函數(shù)，且，。也就是說，在梯度較大的區(qū)域p取1，在梯度較小的區(qū)域p取2。

這里，定義p=2-M，且M定義為［5］

式中，m0為先驗(yàn)圖像，G為高斯卷積核，λ為正則化參數(shù)，x′為像素位置，為高斯平滑后的初始圖像的梯度模值的直方圖中元素個(gè)數(shù)最多的值，*代表卷積運(yùn)算。

從式（3）可以推出，在梯度值較小區(qū)域，表達(dá)式分母中的分子項(xiàng)值較小，M接近于0；在梯度值較大區(qū)域，表達(dá)式分母中分子項(xiàng)值較大，M接近于1。根據(jù)前述基于L1范數(shù)和L2范數(shù)正則化重建方法的優(yōu)缺點(diǎn)，為了得到較好的結(jié)果，應(yīng)該在梯度較大的邊緣區(qū)域采用TV約束方式（2），在平滑區(qū)域采用Tikhonov約束方式。模型在梯度較大區(qū)域p=1，在平滑區(qū)域p=2，符合要求。

另外，由于梯度模值較小的部分多是較平滑的區(qū)域，梯度模值較大的部分往往是邊緣區(qū)，如果在這兩種區(qū)域里也采用自適應(yīng)約束，則擴(kuò)散速度較慢。筆者借鑒PM模型的思想，進(jìn)一步調(diào)整約束方式，根據(jù)圖像的不同設(shè)定兩個(gè)邊界值，對于梯度模值小于較小邊界值的部分直接采用Tikhonov約束方式，梯度模值大于較大邊界值的部分直接采用TV約束方式，其余部分采用自適應(yīng)約束方式，即

式中：α1、α2為設(shè)定的兩個(gè)邊界值，可以根據(jù)圖像效果選擇；p的值可以在迭代過程中根據(jù)新得到的中間結(jié)果更新，也可以由初始圖像給出后不再更新，比較簡單的方法是中間過程不再進(jìn)行p的更新。

2 方法

對于多通道笛卡爾坐標(biāo)系采樣數(shù)據(jù)，SENSE可以對混疊圖像逐點(diǎn)處理，求其對應(yīng)系統(tǒng)矩陣的偽逆，得出原始圖像中發(fā)生混疊的像素值。但是，對于非笛卡爾坐標(biāo)系采樣數(shù)據(jù)，由于無法直接估計(jì)出混疊像素點(diǎn)之間的相互作用，所以不能使用上述方法求解。筆者采用非線性共軛梯度下降方法迭代解方程［6，13］，流程如圖1所示。

圖1 重建過程流程Fig.1 Schematic diagram of reconstruction process

此過程即將敏感度加權(quán)后的圖像首先進(jìn)行傅里葉變換，然后按非笛卡爾采樣軌跡進(jìn)行k空間數(shù)據(jù)重采樣，見圖2（a）。另外，有

此過程先將數(shù)據(jù)網(wǎng)格化至笛卡爾坐標(biāo)系，然后進(jìn)行傅里葉逆變換至圖像域得到混疊圖像，再對混疊圖像進(jìn)行敏感度共軛加權(quán)求和，見圖2（b）。所以，E*并不是簡單的E的逆［7，17］。

在上述兩個(gè)過程中，涉及網(wǎng)格化和逆網(wǎng)格化。網(wǎng)格化即利用卷積核將非笛卡爾軌跡數(shù)據(jù)插值到笛卡爾坐標(biāo)系，逆網(wǎng)格化過程即將笛卡爾軌跡數(shù)據(jù)插值到非笛卡爾采樣軌跡。卷積核通常采用Kaiser-Bessel窗［8-12］，有

圖2 圖像域與k空間域數(shù)據(jù)變換過程。（a）圖像域數(shù)據(jù)至k空間域非笛卡爾軌跡數(shù)據(jù)的變換過程；（b）k空間域非笛卡爾軌跡數(shù)據(jù)至圖像域數(shù)據(jù)的變換過程Fig.2 Transformation of data between image domain and k-space.（a）transformation of data from image domain to non-Cartesian trajectory k-space；（b）transformation of data from non-Cartesian trajectoryk-space to image domain

另外，由于受k空間采樣模式的影響，得到的圖像往往會(huì)變模糊，如對螺旋采樣軌跡數(shù)據(jù)來說會(huì)產(chǎn)生螺旋環(huán)狀邊葉。所以，要將得到的圖像進(jìn)行去極化，即除以卷積核的傅里葉變換。這樣做盡管會(huì)使圖像產(chǎn)生一些陰影，但是它能很好地抑制由于數(shù)據(jù)和采樣函數(shù)卷積產(chǎn)生的邊葉。去極化核近似為

卷積核的寬度和采樣密度相關(guān)，比較寬的核在成像中往往有平滑信息的作用，這樣也可以平滑掉一部分噪聲，比較窄的核有銳化圖像的作用，可以更好地保留小的細(xì)節(jié)信息，但是不利于消除噪聲。在本實(shí)驗(yàn)中，取L=4，α=2。

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

實(shí)驗(yàn)中采用的數(shù)據(jù)來自2010 ISMRM會(huì)議，為人體動(dòng)靜脈畸形瘤動(dòng)脈注射的X線照射結(jié)果。數(shù)據(jù)采樣軌跡為可變密度螺旋軌跡，數(shù)據(jù)采集方式為8通道并行采集，每個(gè)通道采集200條螺旋，每條螺旋上采集2 000個(gè)數(shù)據(jù)點(diǎn)，線性欠采樣螺旋按golden angle角度旋轉(zhuǎn)，在13個(gè)TR后，在k空間邊緣的欠采樣因子近似為15。各個(gè)線圈的敏感度值由軸向?qū)佑?jì)算得來，是已經(jīng)提供的數(shù)據(jù)。在每個(gè)通道采集的數(shù)據(jù)中，都加入了獨(dú)立噪聲。圖像大小設(shè)定為512像素×512像素，本實(shí)驗(yàn)采用了其中的50個(gè)（序號(hào)61～110）個(gè)螺旋作為采樣數(shù)據(jù)，欠采樣倍數(shù)2.6。分別用平方和（sum-of-squares，SOS）重建方法、無約束SENSE重建方法、TV約束SENSE重建方法，以及所提出的自適應(yīng)約束方法進(jìn)行了重建，結(jié)果如圖3～圖5所示。

從實(shí)驗(yàn)結(jié)果中可以得出，同平方和重建方法、無約束SENSE重建方法相比，約束重建方法得到的結(jié)果邊緣更加清晰，與傳統(tǒng)的約束重建方法相比，本研究提出的算法同時(shí)保持了去噪及保護(hù)細(xì)節(jié)的功能。圖4和圖5顯示，無約束重建結(jié)果中比較模糊的小細(xì)節(jié)在本重建算法中得到保護(hù)，平滑區(qū)域的噪聲也得到抑制。但是，也可以從圖中看到，部分小細(xì)節(jié)模糊，這主要是由于懲罰項(xiàng)約束方式p值的大小受先驗(yàn)圖像（見圖3（a））梯度的影響。如果先驗(yàn)圖像中某些部位圖像質(zhì)量過差，那么相應(yīng)部位的p值會(huì)受到影響，從而可能導(dǎo)致相應(yīng)部位重建圖像質(zhì)量降低。在本實(shí)驗(yàn)中，先驗(yàn)圖像在圖4所示部位含有的噪聲較多，而在圖5所示部位的含噪聲較少，p值受到先驗(yàn)圖像中噪聲影響，從而圖4（d）中的部分小細(xì)節(jié)恢復(fù)得不是很好，但是圖5（d）中的小細(xì)節(jié)普遍恢復(fù)得很好。

圖3 4種方法重建結(jié)果。（a）平方和重建結(jié)果；（b）無約束SENSE重建結(jié)果；（c）TV約束SENSE重建結(jié)果；（d）自適應(yīng)約束SENSE重建結(jié)果Fig.3 Reconstruction results of four methods.（a）the result of SOS method；（b）the result of SENSE；（c）the result of SENSE with TV constraint；（d）the result of SENSE with adaptive constraint

圖4 圖3（a）中實(shí)線矩形框內(nèi)對應(yīng)部分的放大圖。（a）平方和重建結(jié)果部分放大；（b）無約束SENSE重建結(jié)果部分放大；（c）TV約束SENSE重建結(jié)果部分放大；（d）自適應(yīng)約束SENSE重建結(jié)果部分放大Fig.4 Enlargement of solid rectangular ROI in Fig.3（a）.（a）the enlargement of SOS method；（b）the enlargement of SENSE；（c）the enlargement of SENSE with TV constraint；（d）the enlargement of SENSE with adaptive constraint

圖5 圖3（a）中虛線矩形框內(nèi)對應(yīng)部分的放大圖。（a）平方和重建結(jié)果部分放大；（b）無約束SENSE重建結(jié)果部分放大；（c）TV約束SENSE重建結(jié)果部分放大；（d）自適應(yīng)約束SENSE重建結(jié)果部分放大Fig.5 Enlargement of dotted rectangular ROI in Fig.3（a）.（a）the enlargement of SOS method；（b）the enlargement of SENSE；（c）the enlargement of SENSE with TV constraint；（d）the enlargement of SENSE with adaptive constraint

4 結(jié)論

在多通道磁共振并行成像中，非笛卡爾軌跡部分采樣k空間數(shù)據(jù)重建仍是一個(gè)難點(diǎn)，在重建中需要考慮噪聲和偽影的去除以及計(jì)算量問題。傅里葉變換和部分采樣多通道非笛卡爾軌跡并行成像技術(shù)能大幅度地減少梯度編碼步數(shù)，提高成像掃描速度，相比單線圈成像技術(shù)能更加有效地縮短數(shù)據(jù)獲取時(shí)間，從而達(dá)到優(yōu)化成像的目的。但是，由于部分采樣數(shù)據(jù)的不完整性及敏感度信息獲取的不準(zhǔn)確性，導(dǎo)致單純的并行重建出現(xiàn)偽影和噪聲。筆者在本文中提出利用先驗(yàn)圖像的梯度特征決定懲罰方式的自適應(yīng)約束方法，第一部分介紹了常用的多通道非笛卡爾采樣部分?jǐn)?shù)據(jù)的各種SENSE約束重建方法，第二部分綜合各種約束重建方式的特點(diǎn)提出了自適應(yīng)算法。

實(shí)驗(yàn)結(jié)果顯示，相比原來的方法，所提出的方法可以有效地抑制噪聲、保留細(xì)節(jié)。但是，由于重建方法需要多次迭代完成，所以重建時(shí)間上有一定的延長。如何進(jìn)一步的優(yōu)化算法，是今后研究的一個(gè)重點(diǎn)。

［1］Liu B，King K.Regularized sensitivity encoding（SENSE）reconstruction using bregman iterations［J］.Magnetic Resonance in Medicine，2009，61（1）：145-152.

［2］Pruessmann KP，Weiger.M.SENSE：Sensitivity Encoding for fast MRI［J］.Magnetic Resonance in Medicine，1999，42（5）：952-962.

［3］陳武凡.并行磁共振成像的回顧、現(xiàn)狀與發(fā)展前景［J］.中國生物醫(yī)學(xué)工程學(xué)報(bào)，2005，24（6）：649-654.

［4］Vijayakumar S，Duensing R，Huang F.G-factor and gradient weighted denoising with edge restoration（g-denoiser）for SENSE reconstruction MR images［C］//Christophe J，Marin O，eds.2008 5th IEEE International Symposium on Biomedical Imaging：from Nano to Macro.Paris：IEEE，2008：472-475.

［5］Perona P，Malik J.Scale space and edge detection using anisotropic diffusion［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1990，12（7）：629-639.

［6］Lustig M，Donoho D，Pauly JM.Sparse MRI：the application of compressed sensing for rapid MR imaging［J］.Magnetic Resonance in Medicine，2007，58（6）：1182-1195.

［7］Pruessmann KP，Weiger M.Advances in sensitivity encoding with arbitrary k-space trajectories［J］.Magnetic Resonance in Medicine，2001，46（4）：638-651.

［8］Rasche V，Proksa R，Sinkus R，et al.Resampling of data between arbitrary grids using convolution interpolation［J］.IEEE Trans Med Imaging，1999，18（5）：385-392.

［9］Jackson JI，Meyer CH，Nishimura DG.Selection of a convolution function for fourier inversion using gridding［J］.IEEE Trans Med Imaging，1991，10（3）：473-478.

［10］Osulllivan J.A fast sinc function gridding algorithm for fourier inversion in computed tomography［J］.IEEE Trans Med Imaging，1985，4（4）：200-207.

［11］?ukur T，Santos JM，Nishimura DG，et al.Varying kernel-extent gridding reconstruction for under sampled variable-density spirals［J］.Magnetic Resonance in Medicine，2008，59（1）：196-201.

［12］Hoge RD，Kwan RKS，Pike GB.Density compensation functions for spiral MRI［J］.MRM，1997，38（1）：117-128.

［13］Beatty PJ，Nishimura DG.Rapid gridding reconstruction with aminimal oversampling ratio［J］.IEEE Trans Med Imaging，2005，24（6）：799-808.

［14］Osher S，Burger M，Goldfarb D.An iterative regularization method for total variation-based image restoration［J］.Multiscale Model Simul，2005，4（2）：460-489.

［15］Huang F，Chen YM，Duensing GR.Application of partial differential equation-based in painting on sensitivity maps［J］.Magnetic Resonance in Medicine，2005，53（2）：388-397.

［16］Rudin LI，Osher S，F(xiàn)atemi E.Nonlinear total variation based noise removal algorithms［J］.Physica，1992，60（1-4）：259-268.

［17］Block KT，Uecker M，F(xiàn)rahm J.Under sampled radial MRI with multiple coils.Iterative image reconstruction using a total variation constraint［J］.Magnetic Resonance in Medicine，2007，57（6）：1086-1098.

［18］Larsson EG，Erdogmus D，Yan R，et al.SNR-optimality of sum of-squares reconstruction for phased-array magnetic resonance imaging［J］.Journal of Magnetic Resonance，2003，163（1）：121-123.