靈敏性編碼重建磁共振成像與非笛卡爾k空間運動軌跡的研究

張連俊， 劉　剛

(山東理工大學 計算機科學與技術學院， 山東淄博 255049)

?

靈敏性編碼重建磁共振成像與非笛卡爾k空間運動軌跡的研究

張連俊， 劉剛

(山東理工大學 計算機科學與技術學院， 山東淄博 255049)

摘要：對敏感性編碼的磁共振成像并行重建方案進行了研究，實現了非笛卡爾采樣k空間軌跡分析.在敏感圖像信息中，由接收器線圈陣列采集的螺旋k空間數據被應用到欠采樣數據集中以重構建無混疊重建圖像.此外，由于網格算法具有前向估值和非正交采樣數據伴隨運算能力，被用來實現敏感性編碼重建圖像.研究發現，敏感性編碼可減少用于磁共振成像實驗的掃描時間，同時保持與下采樣數據集的圖像分辨率.最后，利用真實數據組計算驗證了敏感性編碼的性能，并分析了進一步提高計算速度的方法.

關鍵詞：并行成像; 敏感性編碼; 網格; 磁共振成像

磁共振成像(MRI)是一項新穎的醫療成像診斷技術.但由于制造快速切換的磁場梯度的技術困難，當前MRI掃描儀達到了數據采集速度的物理極限.未來的發展方向是通過減小掃描和重建的時間來提高運行速度.采取平行相控陣列線圈和平行重建方法可達到更快MRI成像的成像技術.1987年相控陣列線圈首次被用于減少相位編碼[1]，可成比例的減少掃描時間量.在隨后二十年MRI研究中，其它并行成像方案也有了長足的進步.在這些方案中最熟知的是敏感性編碼法[2](SENSE)和空間諧波同時采集法[3](SMASH).SENSE與SMASH相比的優點是提供并行成像任意線圈的構造，但在信噪比(SNR)方面有額外損耗.盡管這兩種技術都可以顯著壓縮掃描時間，但是對于欠采樣數據變換圖像，特別是在非正交的數據集的情況下，計算數據的密集的逆過程是必需的.針對此問題有人提出了一種正交重建方案[4]，該方案結合了快速傅立葉變換(FFT)以正向和反向網格運算，計算復雜度由N4降低到N2logN.這種方法使靈敏編碼的影像實現了任意k空間采集模式[5].此外，為了優化SENSE的計算效率，計算中采用迭代方法[6-7].這種方法不會對SENSE成像方案的整體計算性能有很大的影響.本文利用網格化算法的敏感性編碼來實現與非笛卡爾k空間軌跡的重建，對比其與“雙瞳”的數據設置的性能.這些方法將完善目前SENSE的實施，并為實現三維核磁共振成像數據集的重建打下堅實的基礎.

1　基于SENSE的MRI并行成像

磁共振并行成像常常是通過收器線圈陣列收物體數據成像的、收器線圈陣列主要用于提高圖像的信噪比，減少數據采集時間，其中每個獨立空間中的不同朝向的接收器線圈接收不同的數據集.公式(1)表示一個典型的核磁共振實驗中的基帶圖像數據.

ρ(t)=∫d(r)S(r)e-iω(r)te-i2πk(t)tdr+ε(t)

(1)

式中：d(r) 是目標的坐標r處的磁化軌跡連續函數；沿著r的運動軌跡直接從磁共振成像掃描儀在k空間得到，S(r) 是接收器線圈的空間靈敏度分布； ω(r) 是位置r處B0磁場不均勻性；k(t) 是時間t的輸入數據軌跡; ε(t) 是噪聲項.公式(1) 表示k空間數據，線圈靈敏度分布和磁場的不均勻性的信息可用于并行圖像形成三個基本組成部分，從多個獨立的線圈所獲得的數據可以集成到一個系統中.每個線圈接收的數據由Sl(r)加權, l = 1,..., L, 等等. 平行成像模型中的綜合矩陣形式因此可表示如下(取l=3為例)：

(2)

式中：ρ1表示信號矢量從線圈1接收到的數據和圖像矢量ρ通過疊加ρ1形成； SL是對角矩陣保持從第1個線圈的空間靈敏度曲線；F表示成像系統矩陣.

從上面的數學模型可知，并行成像是利用信號處理方法，減少掃描時間并獲得k空間中的數據量，同時仍保持相同的空間分辨率，即k空間的相同區域仍有效的覆蓋奈奎斯特敏感點與從相控陣列線圈所獲取的數據.根據采樣理論通用信號處理知識，任何并行成像方案將使得相位編碼線的數目減少，之間的間隔是1/視場(FOV).在疊加計算重構圖像時，采用了傳統的FFT變換.圖1中完全采樣數據集和欠采樣數據集重建的兩個圖像，相互比較證明了相位編碼線數目的減少.

圖1　完全采樣和欠采樣k空間分析

SENSE并行成像方案[5]結合上述疊加效應，可使多個線圈混疊信號具有更快的成像速度.以笛卡爾采樣k空間為例，相位編碼線之間的位置和距離在MRI實驗前通常已知.從位置和距離能夠準確的計算疊加在每個線圈的重構圖像的量.假定掃描時間減小，加速因子為R，在FOV中的重建的圖像將被縮小或折疊.對于相同因子R，FOV/R的位移像素值可得到一個重構象素的信號強度，疊加在目標對象疊加位置的強度為R-1倍.這使得SENSE能夠求出圖像的疊加強度并篩選出所需的值，如圖2所示.

圖2 SENSE并行成像方案數據重構

SENSE成像方案：為了獲得全視場圖像數據所需的圖像，像素的圖像值Ik(x，y)在相應的位置加權，線圈靈敏度分布為(x,yl)，其中l= 1，...，R，在位置(x，y)的與第k個線圈的信息的準確重構象素值如下：

Ik(x,y)=Sk(x,y1)ρ(x,y1)+Sk(x,y2)ρ(x,y2)+

…+Sk(x,yR)ρ(x,yR).

(3)

該方程可以轉化成矩陣表示法，并形成重構圖像，顯示如圖3所示.

(4)

為了簡化計算公式，在公式(4)中暫不考慮噪聲相關性問題.根據矩陣理論，得出視場圖像的靈敏度矩陣S.

(5)

從方程(5)可以看到，SENSE過程的“展開”功能由從廣義矩陣求逆來實現.另外，值得注意的是SENSE重建算法需要FOV的圖像的每一個像素位置，最終重構一個完整的視場圖像.在SENSE重建的實際實施中，由于存儲器的限制、計算速度瓶頸的影響，不采用傳統 的矩陣S和SH廣義矩陣求逆的解決方案，而采用以共軛梯度迭代的方法求解.

2　網格化算法

該網格算法通常采用不屬于在k空間定期笛卡爾網格數據集進行處理.網格化方法允許收集和處理非笛卡爾數據集用于笛卡爾數據重建方法.如基于FFT的方案的運算方式，可以在沒有多大修改時仍能應用.網格算法中所獲得的螺旋形軌跡的數據是由第一次重新采樣到笛卡爾采樣網格的FFTW庫來完成圖像重建的.在重新采樣的過程中，網格算法提供近似的伴隨算子和向前計算算子.伴隨算子和向前計算算子的概念表示如下：

(6)

(7)

式中：d(m) 為降界的k空間采樣的數據；a(t) 表示插值窗函數的時間分割值；f(xn) 表示的是B0場的不均勻性參數；km表示第k-空間中的軌跡值.在這兩個計算途徑中，共軛算子計算是以圖像數據的值和k-空間的軌跡作為輸入，以笛卡爾網格點的k空間數據值作為輸出.用向前算子計算時與此相反，取k空間數據的值和k-空間的軌跡作為輸入，并輸出每個對應的像素位置的圖像值.雖然這兩種運算不同，但圖像重建的評值時都采用了共軛梯度法.

伴隨算子的計算可以有以下三個步驟：第一步，在k空間中的每個數據點進行采樣，以任意的采樣函數來采樣數據點，然后卷積用網格化的內核.在本文中利用凱塞貝塞爾函數被用作網格內核.第一步可以表示為下面的等式：

(8)

式中：d(kx,ky) 表示原始k空間數據點；S(kx,ky) 表示第k-空間采樣函數；K(kx,ky) 表示凱澤貝塞爾網格內核. 理想的圖像d(x,y) 實際上是通過與傅里葉卷積網格化模糊化的凱撒貝塞爾網格化內核來實現重建.第二步，對網格卷積的結果進行傅里葉變換，然后在笛卡爾網格點重新采樣.這個步驟是將原始的螺旋k空間數據點和重新采樣值變換到到相鄰的笛卡爾網格點，在圖像域中的相應的計算如下：

(9)

式中：d(x,y) 表示理想的圖像數據；s(x,y) 表示圖像域取樣函數；k(x,y) 表示在圖像領域的凱撒貝塞爾網格內核；ψ(x,y) 表示圖像域重新采樣函數.在第二步會出現不理想的因素，即螺旋形旁瓣通常引起采樣圖案的點擴散.為了抑制旁瓣的引入，在第三步驟修正圖像的傅立葉變換的網格內核.

(10)

第三步引入陰影中的圖像從而抑制從第二步中產生旁瓣.最后計算結果即是以共軛算子網格算法計算出的重建圖像.

向前算子網格算法即是實現逆網格化，反轉網格化計算中的每一個步驟.先從共軛梯度求解圖像變換值，再估計在k空間相應的變換數據.從上面的網格化的過程中，很明顯看出逆網格化的第一步驟應該是變跡，即取消在網格算法的最后一步的修正圖像的效果.

(11)

逆網格的其余步驟是簡單地在k空間中卷積和重新采樣.與在共軛算子的那些對應步驟相對，這些步驟的唯一區別是從笛卡爾網格點估算在k空間中的螺旋軌跡的值.在k空間中數據的表達式如下.

K(kx,ky)]S(kx,ky)

(12)

有些文獻[8-9]說明，密度補償是必要的.而本文中使用的特定數據集沒有采用密度補償，而特定數據集采用密度補償重建圖像的分辨率，信噪比沒有顯著改善.在混疊失真較大時應加補償，在計算中加上補償因子即可實現.

3　數據分析與實驗結果

圖像重建的數據設定的是來自“ISMRM重建”的“雙瞳”數據[10].12軸位圖像的“雙瞳”的數據源在腹部用軀干相控陣線圈采集.原始矩陣尺寸為320×320，圖像采集采用梯度回波圖像TE=3方式采集，用96×96采集矩陣調整為四種不同的尺寸，分別為64×64矩陣、128×128矩陣、256×256矩陣和512×512矩陣，以測試重建的有效性.同時還提供數據集與原來的數據集合成8螺旋行軌跡，具有20000像素點.k空間軌跡重建前后的圖像見圖4和圖5，圖中的灰色區域表示交錯覆蓋序列中的區域.四種不同尺寸的重構圖像顯示見圖6至圖9.

圖4 圖像重建前掃描軌跡圖5 圖像重建后掃描軌跡

圖6　64× 64圖像重建圖 7 　128× 128圖像重建

圖8　256× 256圖像重建圖9　512× 512圖像重建

上述實驗分析表明，采用網格算法的SENSE并行成像方案能夠有效地處理非笛卡爾k空間軌跡.在圖7和圖8中，對應于數據集的矩陣大小的圖像分辨率比圖6有所增加.然而在圖8和圖9之間圖像的分辨率提高不顯著，這是由于原始采集的數據集僅有320×320矩陣的大小.通過內插在k空間縮放矩陣大小不會顯著改善圖像分辨率.表1給出了SENSE方法實現的計算效率.估值計算時間大約是實際計算時間的90%以上，計算時間的差異主要是由于存儲器大小影響.

表1SENSE實現的計算效率

圖像尺寸/像素伴隨算子/s前向算子/s估值時間/s實際用時/s估值時間實際用時/%512×51247.5342.3410992.012056.691.17256×25615.9914.153690.43827.496.41128×1285.133.241080.01170.992.2464×641.760.76358.4372.396.24

4 結束語

利用網格化算法的SENSE并行成像，可使從任意k空間軌跡重建圖像的效果提高，并提高重建圖像的速度.通過來自多個接收器線圈的采樣數據，可以有效地補償失真.采用共軛梯度算法、正向和反向網格和FFT變換等組合方法可以降低計算復雜度.然而從圖6到圖9中可以看出，重構圖像中包含少量的來源于給定數據集中的高斯白噪聲，這可以通過將正則化與SENSE相結合來消除.這將不斷完善目前SENSE處理方法，并為實現SENSE三維核磁共振成像數據重建打下堅實的基礎.

參考文獻:

[1]CarlsonJ.AnalgorithmforNMRimagingreconstructionbasedonmultipleRFreceivercoils[J].MagnReson, 1987,74: 376-380.

[2]PruessmannK.SensitivityencodingforfastMRI[J].MagnResonMed, 1999,42: 952-962.

[3]SodicksonD.Simultaneousacquisitionofspatialharmonicsfastimagingwithradiofrequencycoilarrays[J].MagnResonMed, 1997, 38: 591-603.

[4]PruessmannK.Advancesinsensitivityencodingwitharbitaryk-spacetrajectories[J].MagnResonMed, 2001,46: 638-651.

[5]BlaimerM,SmashS,GrappaP.Howtochoosetheoptimalmethod[J].TopMagnResonImaging,2004,15:223-236.

[6]LinF.Parallelimagingreconstructionusingautomaticregularization[J].MagnResonMed, 2004,51:559-567.

[7]FesslerJ.Toeplitz-basediterativeimagereconstructionforMRIwithcorrectionformagneticfieldinhomogeneity[J].IEEETransSigProc, 2005,53: 3393-3402.

[8]PaulyJ.Non-CartesianReconstruction[J].Lecturenote, 2007, 86:759-767.

[9]FangS,YingK,ZhaoL,etal.CohereregularizationforSENSEreconstructionwithanonlocaloperator[J].MagneticResonanceReview, 2010, 64(5):1414-1426.

[10]WuB,MillaneR,WattsR, et al.Priorestimate-basedcompressedsensinginparallelMRI[J].MagneticResonanceReview,2011,2011,65(1):83-95.

(編輯：姚佳良)

Researchofsensitivityencodingreconstructionfor

MRIwithnon-cartesiank-spacetrajectories

ZHANGLian-jun,LIUGang

(SchoolofComputerScienceandTechnology,ShandongUniversityofTechnology，Zibo255049,China)

Abstract:The sensitivity encoding (SENSE) parallel reconstruction scheme for magnetic resonance imaging (MRI) was studied and implemented with non-Cartesian sampled k-space trajectories in this report. The sensitivity map profile, field map information and the spiral k-space data collected from an array of receiver coils were used to reconstruct un-aliased images from under-sampled data. The gridding algorithm was also implemented with SENSE due to its ability in evaluating forward and adjoins operator with non-Cartesian sampled data. SENSE has the special capability to reduce the scanning time for MRI experiments while maintaining the image resolution with under-sampling data sets.We also analyzed the performance of SENSE with real data set and identified the computational issues that need to be improved for further research.

Key words:parallel imaging; SENSE; gridding; MRI

中圖分類號：391.4

文獻標志碼：A

文章編號：1672-6197(2015)04-0021-04

作者簡介：張連俊，男，zhlj@sdut.edu.cn

收稿日期：2014-09-29