非笛卡爾并行磁共振成像數據的自適應約束重建新算法

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

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

引言

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

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

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

SENSE重建TV約束模型為freg=

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

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

1 自適應約束SENSE重建理論

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

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

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

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

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

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

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

2 方法

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

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

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

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

在上述兩個過程中，涉及網格化和逆網格化。網格化即利用卷積核將非笛卡爾軌跡數據插值到笛卡爾坐標系，逆網格化過程即將笛卡爾軌跡數據插值到非笛卡爾采樣軌跡。卷積核通常采用Kaiser-Bessel窗［8-12］，有

圖2 圖像域與k空間域數據變換過程。（a）圖像域數據至k空間域非笛卡爾軌跡數據的變換過程；（b）k空間域非笛卡爾軌跡數據至圖像域數據的變換過程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án)狀邊葉。所以，要將得到的圖像進行去極化，即除以卷積核的傅里葉變換。這樣做盡管會使圖像產生一些陰影，但是它能很好地抑制由于數據和采樣函數卷積產生的邊葉。去極化核近似為

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

3 實驗結果與分析

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

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

圖3 4種方法重建結果。（a）平方和重建結果；（b）無約束SENSE重建結果；（c）TV約束SENSE重建結果；（d）自適應約束SENSE重建結果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）中實線矩形框內對應部分的放大圖。（a）平方和重建結果部分放大；（b）無約束SENSE重建結果部分放大；（c）TV約束SENSE重建結果部分放大；（d）自適應約束SENSE重建結果部分放大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）中虛線矩形框內對應部分的放大圖。（a）平方和重建結果部分放大；（b）無約束SENSE重建結果部分放大；（c）TV約束SENSE重建結果部分放大；（d）自適應約束SENSE重建結果部分放大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 結論

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

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

［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］陳武凡.并行磁共振成像的回顧、現狀與發(fā)展前景［J］.中國生物醫(yī)學工程學報，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，Fatemi E.Nonlinear total variation based noise removal algorithms［J］.Physica，1992，60（1-4）：259-268.

［17］Block KT，Uecker M，Frahm 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.