駱國麗，肖亮

北京化工大學(xué)信息科學(xué)與技術(shù)學(xué)院，北京 100029

國內(nèi)低場磁共振成像系統(tǒng)大多采用永磁磁體，相較于高場超導(dǎo)系統(tǒng)，永磁系統(tǒng)場強(qiáng)低(場強(qiáng)一般為0.5 T及以下)，并且磁場均勻性和穩(wěn)定性都比較差，磁場強(qiáng)度受溫度的影響也比較大。主磁場非均勻性對成像的影響程度一直是研究人員、技術(shù)人員與使用者關(guān)心的問題，以往主要靠定性分析，缺乏直觀的圖像表達(dá)。

自1984年Bittoun等[1]建立第一個基于Bloch方程的MRI模擬器以來，迄今國內(nèi)外對其進(jìn)行了大量廣泛深入的研究，主要集中在以下幾個方面：模擬器從1D到2D、3D推廣及應(yīng)用[2-3]；組織解剖模型的建立[4-6]，根據(jù)需求可分為離散組織模型和模糊組織模型；為了加快計(jì)算速度，模擬器算法的優(yōu)化及并行化[7-9]。這些研究取得了良好的效果，使其可以應(yīng)用于低成本大規(guī)模磁共振成像教學(xué)、磁共振實(shí)驗(yàn)設(shè)計(jì)、MRI序列編寫、參數(shù)優(yōu)化、系統(tǒng)偽影等問題的解決[10-13]。不過，對于低場中主磁場不均勻性的模擬一直沒有系統(tǒng)的研究與模擬。

本文實(shí)現(xiàn)一種基于Bloch方程的MRI模擬器，通過對主磁場不均勻性分布進(jìn)行建模，提出了對主磁場不均性的模擬仿真方案。以自旋回波(spin echo，SE)和梯度回波(gradient echo，GRE)序列成像方法為例對主磁場的不均勻性進(jìn)行了模擬，并得到了預(yù)期的實(shí)驗(yàn)結(jié)果。

1 材料和方法

1.1 模擬器

模擬器構(gòu)成框圖如圖1所示，包括虛擬物體、靜磁場定義、MRI序列、磁化強(qiáng)度矢量計(jì)算核心，計(jì)算出射頻(radio frequency，RF)信號(k空間信號)，加入噪聲，最后用快速傅立葉變換(fast Fourier transform，F(xiàn)FT)重建圖像。

虛擬物體是對真實(shí)自旋系統(tǒng)的離散化描述。虛擬物體的每一個體素包含一組用于通過Bloch方程計(jì)算其自旋磁化矢量的物理量：質(zhì)子密度(ρ)、兩個弛豫常量T1和T2。

MRI序列主要包含3個事件：自由進(jìn)動、RF脈沖激發(fā)、信號采集。在MRI模擬器中通過排列組合以上3種事件、設(shè)置事件參數(shù)得到不同MRI序列。MRI序列實(shí)施過程中的磁化強(qiáng)度矢量計(jì)算是基于簡化的3D Bloch方程的解，Bloch方程描述了自旋磁化矢量隨時間的變化，如式(1)所示：

其中M0是自旋磁化矢量的平衡常數(shù)，由質(zhì)子密度決定；(T1,T2)是弛豫時間；γ是物體成分的旋磁比。磁場強(qiáng)度由式(2)表示：

圖1 MRI模擬器構(gòu)成框圖Fig.1 The overview frame of the MRI simulator

模擬計(jì)算核心通過旋轉(zhuǎn)矩陣和指數(shù)擴(kuò)展執(zhí)行Bloch方程的離散時間解，得到演變的磁化強(qiáng)度矢量。對于RF信號采集，一次激發(fā)可以得到一維離散復(fù)數(shù)信號，可以填充k空間的一條相位編碼線。RF信號的一點(diǎn)S[t]是對應(yīng)時間點(diǎn)虛擬物體所有體素的磁化強(qiáng)度矢量之和，如式(3)所示：

1.2 系統(tǒng)參數(shù)設(shè)置

創(chuàng)建一個2D虛擬水模，對模擬器的2D模擬成像效果進(jìn)行實(shí)驗(yàn)驗(yàn)證。2D虛擬水模模型：FOV 250 mm×250 mm，采樣矩陣和重建矩陣同為256×256；水模大小150 mm×150 mm，弛豫時間T1=100 ms，T2=100 ms(參考實(shí)驗(yàn)用水模：1L水中加入3.6 g NaCl和1.995 g CuSO4·5H2O)；主磁場強(qiáng)度B0為0.4 T；在k空間加入信噪比為40 db的高斯白噪聲。分別用SE序列和GRE序列進(jìn)行實(shí)驗(yàn)，SE序列參數(shù)為：TR 300 ms，TE 20 ms，采樣時間間隔為48 us；GRE序列參數(shù)為：TR 80 ms，TE 11.6 ms，射頻脈沖翻轉(zhuǎn)角為70°，采樣時間間隔為36 us。虛擬成像中不考慮渦流的影響，梯度的上升下降時間設(shè)置為零。在實(shí)際成像實(shí)驗(yàn)中除將梯度的上升下降時間均設(shè)為532 us外，其他參數(shù)設(shè)置與上述參數(shù)相同。

1.3 主磁場不均勻性建模

用0BΔ表示系統(tǒng)的勻場精度，(x,y)表示體素點(diǎn)的坐標(biāo)位置，圖像大小為N×N。永磁系統(tǒng)中主磁場的非均勻分布可以用基于球坐標(biāo)的解析法分解為一系列線性無關(guān)的空間諧波的線性組合[14]，在直角坐標(biāo)系中可以展開為一個關(guān)于空間坐標(biāo)的多項(xiàng)式。據(jù)此假定主磁場不均勻性分兩種模式：線性分布、拋物線形式分布。

模式一：線性分布函數(shù)如式(4)所示，當(dāng)場不均勻性沿x軸線性變化時，i代表體素的x坐標(biāo)；當(dāng)場不均勻性沿y軸線性變化時，i代表體素的y坐標(biāo)。

模式二：拋物線形式的分布函數(shù)如式(5)所示：

2 結(jié)果

2.1 主磁場均勻條件下的模擬成像

設(shè)置 ΔB0=0ppm，圖2給出了實(shí)際實(shí)驗(yàn)和虛擬MRI的SE和GRE序列成像結(jié)果。在實(shí)際成像中，因?yàn)閳A柱形水模中存在氣泡，實(shí)驗(yàn)結(jié)果中存在缺口。

實(shí)際實(shí)驗(yàn)中因?yàn)樯漕l場的不均勻、梯度的非線性、接收線圈的敏感度變化等因素得到的磁共振圖像信號強(qiáng)度并不是均勻的。相比而言，因?yàn)樘摂M成像中設(shè)定的都是理想化的實(shí)驗(yàn)條件，所以圖像各像素點(diǎn)的信號強(qiáng)度基本相同。

2.2 主磁場不均勻條件下的模擬成像

圖3分別給出了模式一下SE和GRE序列的虛擬圖像，其中A、C表示場不均勻性沿x軸(頻率編碼方向)線性變化時的圖像，B、D表示場不均勻性沿y軸(相位編碼方向)線性變化時的圖像。從圖3可以看出，隨0BΔ增大圖像幾何失真逐漸增強(qiáng)。當(dāng)0BΔ增加到50 ppm，圖像幾何失真會非常明顯。現(xiàn)在常用的永磁系統(tǒng)勻場精度能夠達(dá)到30 ppm，對圖像的影響一般視覺上是不可見的。脈沖序列成像對頻率編碼方向的場不均勻性較敏感，對相位編碼方向的場不均勻性不敏感。

實(shí)際實(shí)驗(yàn)中，分別將x方向的梯度值設(shè)定為正常梯度的99%、y方向梯度不變，y方向的梯度值設(shè)定為正常梯度的101%、x方向梯度不變，將SE序列和GRE序列的成像與相應(yīng)的虛擬成像比較，結(jié)果如圖4所示。它們的幾何失真幅度和模式幾乎相同，可見模擬器對線性主磁場不均勻具有良好的模擬效果。

圖2 SE和GRE序列的MRI圖像。A：虛擬成像結(jié)果；B：實(shí)驗(yàn)結(jié)果Fig.2 MRI images obtained with SE and GRE sequences.A:the simulated imaging results; B:the experimental imaging results.

圖3 模式一下的虛擬成像。A、B為SE序列圖像，C、D為GRE序列圖像；A、C為不均勻性沿x軸線性變化，B、D為不均勻性沿y軸線性變化；deltaB0表示勻場精度Fig.3 The virtual MRI with various shimming accuracy(deltaB0)of model one.A,B were obtained by SE sequence; C,D were obtained by GRE sequence; The main field inhomogeneity linearly related to the corresponding axis,A and C were x axis,B and D were y axis.

在場不均勻度相同的情況下，圖像的畸變會隨采樣帶寬增大而減小。圖5給出了在模式一x方向上的勻場精度為125 ppm情況下，采樣帶寬分別為20.83 kHz和41.67 kHz時SE序列的虛擬成像結(jié)果。在圖像中，頻率編碼方向水模大小從20.83 kHz采樣時的129像素變?yōu)?1.67 kHz采樣時的121像素，圖像畸變減小。

圖6給出了模式二下的SE和GRE序列虛擬成像。圖像的信號強(qiáng)度發(fā)生變化，并且GRE序列中信號變化比SE序列明顯。由于180°重聚相位脈沖存在，SE脈沖序列成像對場不均勻性引起的信號衰減不敏感。

圖4 虛擬實(shí)驗(yàn)(A)與實(shí)際成像實(shí)驗(yàn)(B)對比。(1)、(2)為SE序列成像，(3)、(4)為GRE序列成像；(1)、(3)為改變x方向的梯度，(2)、(4)為改變y方向的梯度Fig.4 Comparison of(A)experimental and(B)simulated images.(1)and(2)were SE sequence images,(3)and(4)were GRE sequence images;(1)and(3)were images with shimming of x_axis,(2)and(4)were images with shimming of y_axis.

圖5 采樣帶寬對MRI圖像的影響。采樣帶寬：A：20.83 kHz，B：40.67 kHz 圖6 模式二下的虛擬成像。A：SE序列圖像，B：GRE序列圖像，deltaB0表示勻場精度Fig.5 Impact of the bandwidth to MRI.the bandwidth of A is 20.83 kHz,the bandwidth of B is 40.67 kHz.Fig.6 The virtual MRI with various shimming accuracy(deltaB0)of model two.A was obtained by SE sequence,B was obtained by GRE sequence.

3 討論

通過與實(shí)際實(shí)驗(yàn)成像的對比可得本模擬器對線性主磁場不均勻性具有較好的模擬效果。不過在模擬中，只考慮了在空間上的場不均勻性，而沒有提到隨時間變化的場不均勻性?？紤]磁場隨時間的變化需要為磁場分布增加一個時間維度，增加了建模的復(fù)雜度，操作比較困難，但這可以是以后的一個研究方向，以模擬實(shí)際磁共振成像中主磁場隨溫度等周圍環(huán)境變化的情況。

實(shí)際MRI過程中，圖像還會受到RF的不均勻性、梯度的非線性及采集線圈的空間敏感度等因素的影響，因?yàn)閷λ鼈兊闹苯幽M不容易實(shí)現(xiàn)，本文沒有涉及。為了更真實(shí)準(zhǔn)確地模擬MRI成像過程，可以考慮將這些因素通過一些應(yīng)用的選擇實(shí)現(xiàn)允許用戶加入到模擬器中。

基于簡化的3D Bloch方程的離散時間解來計(jì)算k空間磁化矢量的方法相當(dāng)消耗時間。由式(3)可見，對于包含體素?cái)?shù)為N2、k空間柵格大小同樣為N2的情況，算法的復(fù)雜度為0(N4)。由此可見隨著虛擬物體包含的體素增加，算法的計(jì)算量增加會越來越快。對于上面的2D SE序列成像實(shí)驗(yàn)，256×256虛擬物體k空間磁化矢量的計(jì)算時間大概為50 min，所用的處理器為Intel Core i3，主頻為2.4 GHz。同時由式(3)還可以看出，磁化矢量計(jì)算核心具有很好的可并行性，利用GPU、Open-CL等方法實(shí)現(xiàn)程序并行化可以大幅度的減少程序運(yùn)行時間。這也是模擬器的主流研究方向之一。

本文提出的MRI模擬器是一種基于MRI事件的模擬器，它可以通過事件的排列仿真多種MRI序列；可并行化的磁化矢量計(jì)算核心算法，為程序的并行化發(fā)展提供了廣闊的可能性；通過對主磁場不均勻性進(jìn)行空間建模，并以此為帶入條件進(jìn)行模擬成像，提高了模擬器的對場不均勻性模擬的準(zhǔn)確性；可以通過設(shè)置與低場磁共振成像設(shè)備相同的參數(shù)進(jìn)行對實(shí)際情況的模擬仿真，有助于工程與科研中分析主磁場不均勻?qū)Τ上竦挠绊懀瑤椭芯咳藛T更好地進(jìn)行主磁場矯正。

[References]

[1]Bittoun J,Taquin J,Sauzade M.A computer algorithm for the simulation of any nuclear magnetic resonance(NMR)imaging method.Magn Reson Imaging,1984,2(2):113-120.

[2]Olsson MB,Wirestam R,Persson BR.A computer simulation program for MR imaging:application to RF and static magnetic fi eld imperfections.Magn Reson Med,1995,34(4):612-617.

[3]Benoit-Cattin H,Collewet G,Belaroussi B,et al.The SIMRI project:a versatile and interactive MRI simulator.J Magn Reson,2005,173(1):97-115.

[4]Kwan RK,Evans AC,Pike GB.MRI simulation-based evaluation of image-processing and classification methods.IEEE Trans Med Imaging,1999,18(11):1085-1097.

[5]Liu LR,Wu ZQ,Lin JR.The computer simulation of magnetic resonance imaging.J Biomed Eng,2002,19(4):611-615.劉麗榮,吳志嶠,林家瑞.磁共振成像的計(jì)算機(jī)仿真.生物醫(yī)學(xué)工程學(xué)雜志,2002,19(4):611-615.

[6]Herndon RC,Lancaster JL,Toga AW,et al.Quantification of white matter and gray matter volumes for T1 parametric images using fuzzy classifi ers.J Magn Reson Imaging,1996,6(3):425-435.

[7]Benoit-Cattin H,Bellet F,Montagnat J,et al.Magnetic resonance imaging(MRI)simulation on a grid computing architecture.Proceedings of the 3rd IEEE/ACM international symposium on cluster computing and the grid,2003.

[8]St?cker T,Vahedipour K,Pflugfelder D,et al.A high-performance computing MRI simulations.Magn Reson Med,2010,64(1):186-193.

[9]Xanthis CG,Venetis IE,Chalkias AV,et al.MRISIMUL:A GPU-based parallel approach to MRI simulations.Medical Imaging,2013,33(3):607-627.

[10]Kuang B,He CM.A perception on the development of the virtual scan technology for the magnetic resonance imaging.Chin J Magn Reson Imaging,2011,2(3):218-224.匡斌,何超明.磁共振成像虛擬掃描技術(shù)的發(fā)展及趨勢淺談.磁共振成像,2011,2(3):218-224.

[11]Chen SZ,Qi HF,Xie HB,et al.A universal MRI computer simulation software.Chinese Journal of Magnetic Resonance,1996,13(5):491-501.陳思中,戚海峰,謝海濱,等.一種通用的MRI計(jì)算機(jī)模擬軟件.波譜學(xué)雜志,1996,13(5):491-501.

[12]Shkarin P,Spencer R.Direct simulation of spin echoes by summation of isochromats.Concepts in Magnetic Resonance,1996,8(4):253-268.

[13]Hanson LG.A graphical simulator for teaching basic and advanced MR imaging techniques.Radio Graphics,2007,27(6):e27

[14]Vadovic R.Magnetic fi eld correction using magnetized shims.IEEE Trans on Magn,1989,25(4):3313-3319.