張連俊， 劉　剛

(山東理工大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院， 山東淄博 255049)

?

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

張連俊， 劉剛

(山東理工大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院， 山東淄博 255049)

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

關(guān)鍵詞：并行成像; 敏感性編碼; 網(wǎng)格; 磁共振成像

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

1　基于SENSE的MRI并行成像

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

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

(1)

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

(2)

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

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

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

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

圖2 SENSE并行成像方案數(shù)據(jù)重構(gòu)

SENSE成像方案：為了獲得全視場圖像數(shù)據(jù)所需的圖像，像素的圖像值Ik(x，y)在相應(yīng)的位置加權(quán)，線圈靈敏度分布為(x,yl)，其中l(wèi)= 1，...，R，在位置(x，y)的與第k個線圈的信息的準(zhǔn)確重構(gòu)象素值如下：

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

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

(3)

該方程可以轉(zhuǎn)化成矩陣表示法，并形成重構(gòu)圖像，顯示如圖3所示.

(4)

為了簡化計算公式，在公式(4)中暫不考慮噪聲相關(guān)性問題.根據(jù)矩陣?yán)碚摚贸鲆晥鰣D像的靈敏度矩陣S.

(5)

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

2　網(wǎng)格化算法

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

(6)

(7)

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

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

(8)

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

(9)

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

(10)

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

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

(11)

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

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

(12)

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

3　數(shù)據(jù)分析與實驗結(jié)果

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

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

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

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

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

表1SENSE實現(xiàn)的計算效率

圖像尺寸/像素伴隨算子/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 結(jié)束語

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

參考文獻(xiàn):

[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

文獻(xiàn)標(biāo)志碼：A

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

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

收稿日期：2014-09-29