黃清明,鄭 剛

(1.上海健康醫(yī)學(xué)院醫(yī)學(xué)影像學(xué)院，上海 201318；2.上海理工大學(xué)光電信息與計(jì)算機(jī)工程學(xué)院，4.醫(yī)學(xué)影像工程研究所，上海 200093；3.上海市分子影像重點(diǎn)實(shí)驗(yàn)室，上海 201318)

MR成像過程中，梯度線圈在靈敏區(qū)內(nèi)產(chǎn)生期望的梯度磁場(chǎng)強(qiáng)度，但渦流使靈敏區(qū)實(shí)際梯度磁場(chǎng)為激發(fā)梯度磁場(chǎng)和渦流產(chǎn)生磁場(chǎng)的合成，導(dǎo)致脈沖上升和下降沿拖長(zhǎng)，抑制梯度磁場(chǎng)快速變化，影響梯度切換率[1]；空間定位編碼會(huì)使樣品空間定位線性變差，導(dǎo)致圖像畸變和偽影等失真現(xiàn)象，延長(zhǎng)系統(tǒng)回波時(shí)間，減緩MR成像時(shí)間。為加快MR成像速度、獲得高質(zhì)量圖像，必須減弱甚至消除渦流的影響[2-4]。

梯度線圈設(shè)計(jì)的實(shí)質(zhì)是根據(jù)已知磁場(chǎng)分布計(jì)算梯度線圈的繞線模式，使其通電后產(chǎn)生期望磁場(chǎng)強(qiáng)度。1981年，BANGERT等[5]提出梯度線圈最簡(jiǎn)單模型和制作梯度線圈的基本方法，使MRI得以實(shí)現(xiàn)；1993年，TURNER等[6]總結(jié)梯度線圈設(shè)計(jì)的矩陣求逆法[7]、流函數(shù)法[8]和目標(biāo)場(chǎng)點(diǎn)法[9]等，利用目標(biāo)場(chǎng)點(diǎn)法預(yù)設(shè)圓柱表面磁場(chǎng)大小，通過計(jì)算機(jī)輔助軟件計(jì)算圓柱形電流密度分布函數(shù)，經(jīng)離散處理得到梯度線圈的繞線模型。之后，LIU等[10]提出雙平面梯度線圈設(shè)計(jì)方法，計(jì)算出梯度電感能量分布和占電感總值的比例；GREEN等[11]提出單平面梯度線圈設(shè)計(jì)理論，以傅里葉變換電流密度函數(shù)得到梯度線圈功耗系數(shù)最小值。本研究采用目標(biāo)場(chǎng)點(diǎn)法，增加梯度電感最小化和梯度磁場(chǎng)自屏蔽等約束條件建立數(shù)學(xué)模型，利用數(shù)值計(jì)算軟件計(jì)算梯度線圈的理論參數(shù)性能和閉環(huán)繞線模式，以圖形處理單元(graphic processing unit， GPU)加速算法有效解決了上述復(fù)雜數(shù)學(xué)模型中電流密度系數(shù)矩陣運(yùn)算時(shí)間長(zhǎng)的問題，并以電磁場(chǎng)仿真軟件對(duì)設(shè)計(jì)方案進(jìn)行正演驗(yàn)證，可在一定程度上減小梯度線圈電感值和線圈能耗，使梯度線圈具有良好梯度切換率[12]。

1 梯度線圈性能指標(biāo)

MR通過梯度線圈產(chǎn)生的梯度磁場(chǎng)實(shí)現(xiàn)信號(hào)空間定位，主要性能指標(biāo)包括梯度強(qiáng)度、梯度線性度、梯度切換率和梯度電感等[13]。

1.1 梯度線圈性能指標(biāo)

1.1.1 梯度強(qiáng)度 指梯度磁場(chǎng)強(qiáng)度大小，即每單位長(zhǎng)度磁場(chǎng)強(qiáng)度的變化?？焖俪上窈瞳@取高分辨率MRI需要較高梯度強(qiáng)度[14]，梯度強(qiáng)度越高，選擇脈沖序列越靈活。

(1)

式中γ為原子核的磁旋比，N為某個(gè)方向采集的信號(hào)數(shù)，BW為射頻帶寬，Gmax為最大梯度磁場(chǎng)強(qiáng)度，χ為分辨像素的最小限度。選定射頻帶寬BW后，梯度強(qiáng)度Gmax越大，則切面層厚越薄，圖像空間分辨率越高，掃描視野越小。

1.1.2 梯度線性度 指梯度強(qiáng)度變化率的大小，用梯度磁場(chǎng)強(qiáng)度變化的斜率表示。

(2)

式中γmax、γmin和γave分別為靈敏區(qū)內(nèi)各點(diǎn)梯度強(qiáng)度變化斜率最大值、最小值和平均值。梯度線性度決定圖像空間定位的準(zhǔn)確性，λ值越小表示梯度磁場(chǎng)越精確，空間定位越準(zhǔn)確圖像質(zhì)量越好，反之圖像幾何失真度越大。

1.1.3 梯度切換率 指梯度強(qiáng)度從零增加到某一預(yù)定梯度強(qiáng)度值的速度。梯度磁場(chǎng)變化時(shí)，掃描序列處于等待狀態(tài)，梯度切換率增加，則回波時(shí)間和重復(fù)時(shí)間減小，梯度變化速度越快，成像時(shí)間越短[15-16]。

1.1.4 梯度電感 指梯度線圈電感量的大小，主要取決于線圈匝數(shù)、繞制方式和材料等。

(3)

式中μ0為真空磁導(dǎo)率，I為電流，J表示電流密度，v和v′表示電流密度分布的線圈空間。梯度線圈通電時(shí)，線圈導(dǎo)線周圍產(chǎn)生電磁場(chǎng)，處于電磁場(chǎng)范圍內(nèi)的線圈導(dǎo)線發(fā)生“自感”和“互感”作用，在電導(dǎo)體中感應(yīng)出渦流，抑制梯度磁場(chǎng)快速變化。

1.2 梯度線圈的品質(zhì)因素 梯度線圈品質(zhì)因數(shù)Q與梯度線圈效率η、梯度電感L和梯度均勻性δ之間關(guān)系如下[17]：

(4)

2 梯度線圈的設(shè)計(jì)理論

梯度線圈設(shè)計(jì)方法有分離導(dǎo)線法和分布電流密度法。本研究采用分布電流密度法(即目標(biāo)場(chǎng)點(diǎn)法)，加入梯度電感最小化和梯度磁場(chǎng)自屏蔽等約束條件進(jìn)行優(yōu)化。

2.1 梯度線圈目標(biāo)場(chǎng)點(diǎn)法電流密度 設(shè)雙平面梯度線圈位于平面z=±a，梯度線圈半徑ρ滿足ρmin≤ρ≤ρmax，ρmin為線圈最小半徑，ρmax為線圈最大半徑。極坐標(biāo)系上，通電梯度線圈的電流密度J(ρ,φ)分解成徑向分量Jρ(ρ，φ)和切向分量Jφ(ρ，φ)，對(duì)電流密度J(ρ,φ)進(jìn)行傅里葉變換[18]：

(5)

其中，c=π/(ρmax-ρmin)，Uq為電流密度系數(shù)，Q為展開的級(jí)數(shù)，q和k是整數(shù)。k=0時(shí)，為縱向梯度線圈z，產(chǎn)生縱向梯度磁場(chǎng)，電流密度表達(dá)式：

(6)

k=1時(shí)，為橫向梯度線圈x或y，產(chǎn)生橫向梯度磁場(chǎng)，電流密度表達(dá)式：

(7)

k≥2為對(duì)應(yīng)階次的勻場(chǎng)線圈。

2.2 自屏蔽約束雙平面梯度線圈 梯度線圈開合使電導(dǎo)體產(chǎn)生渦流，致梯度磁場(chǎng)不穩(wěn)定，引起圖像畸變[19]。為抑制渦流產(chǎn)生，在梯度線圈外增加一組屏蔽線圈，梯度線圈產(chǎn)生MR梯度磁場(chǎng)，與屏蔽線圈通反向電流產(chǎn)生的磁場(chǎng)疊加，抵消梯度系統(tǒng)渦流產(chǎn)生的磁場(chǎng)，這種結(jié)構(gòu)稱為自屏蔽線圈。設(shè)屏蔽線圈位置為z=±b(b>a)，屏蔽線圈電流密度函數(shù)表達(dá)式類似于梯度線圈的電流密度表達(dá)式，根據(jù)畢奧-薩伐爾定律可知：

(8)

Dq和Eq代表P(x,y,z)的功能函數(shù)，預(yù)設(shè)的目標(biāo)場(chǎng)內(nèi)的點(diǎn)Bz在目標(biāo)場(chǎng)點(diǎn)的靈敏區(qū)和屏蔽區(qū)內(nèi)求解(8)式，得Q和P2個(gè)未知變量。給出目標(biāo)點(diǎn)和對(duì)應(yīng)磁場(chǎng)強(qiáng)度值Bi(i=1， 2， …，Q+P)，將通過(8)式求解出的DiQ和EiQ寫成矩陣形式：

(9)

求得梯度線圈和屏蔽線圈的電流密度函數(shù)，得到兩組線圈繞線模式。

2.3 構(gòu)建梯度電感最小約束函數(shù) 設(shè)計(jì)梯度線圈時(shí)，既要找到最佳電流密度分布，滿足靈敏區(qū)磁場(chǎng)強(qiáng)度要求，更應(yīng)減輕梯度電感對(duì)渦流的影響[20]。梯度線圈為電感L和電阻R組成的電磁場(chǎng)回路，梯度電感L影響電流從零值上升到最大值63%的速度，即時(shí)間常數(shù)τ=L/R。τ越小越好，故應(yīng)增大電阻R或減小電感L，而增大電阻R會(huì)導(dǎo)致線圈功率消耗增加，因此選擇梯度電感L盡可能小。增加梯度電感L最小化約束條件的磁能表示為[21]：

(10)

(11)

其中，U為電流密度表達(dá)式的系數(shù)。用拉格朗日乘子法構(gòu)造梯度電感最小約束函數(shù)：

(12)

對(duì)(12)式電流密度函數(shù)求偏導(dǎo)，推導(dǎo)出電流密度函數(shù)為：

U=W-1D(DTW-1D)-1B

(13)

根據(jù)電流密度散度為零的特點(diǎn)，雙平面線圈的流函數(shù)I(ρ,φ)滿足如下關(guān)系：

(14)

圖1 電流密度系數(shù)矩陣GPU加速運(yùn)算并行運(yùn)算流程

對(duì)于橫向梯度線圈y或x的流函數(shù)表示為：

I=-∑Usin[qc(ρ-ρmin)]cosφ

(15)

對(duì)于縱向梯度線圈z的流函數(shù)表示為：

(16)

電流密度的等量線圖為：

I(ρ,φ)=Imin+(i+1/2)I0， (i=0， 1， 2， …，N-1)

(17)

式中I0=(Imax-Imin)/N，Imax為線圈平面內(nèi)流函數(shù)最大值，Imin為線圈平面內(nèi)流函數(shù)最小值，N為離散的線圈匝數(shù)，繪制(17)式等高線的曲線分布，保證每匝線圈中電流為I0，繞線位置與等高線重合，即得到滿足場(chǎng)點(diǎn)要求的梯度線圈繞線[22]。

3 電流密度系數(shù)矩陣GPU加速運(yùn)算

GPU具有高顯存帶寬、良好的浮點(diǎn)計(jì)算和并行運(yùn)算能力。對(duì)雙平面梯度線圈目標(biāo)場(chǎng)點(diǎn)法加入梯度電感最小和自屏蔽線圈等約束條件，反演運(yùn)算矩陣的求解屬于大數(shù)據(jù)量并行處理計(jì)算過程，執(zhí)行相同流程的并行化浮點(diǎn)運(yùn)算。利用MATLAB平臺(tái)在多核CPU基礎(chǔ)上搭載多塊GPU計(jì)算異構(gòu)，即以分布式并行計(jì)算方式實(shí)現(xiàn)計(jì)算梯度線圈電流密度系數(shù)矩陣[23]。

電流密度系數(shù)矩陣GPU加速運(yùn)算并行運(yùn)算流程見圖1。CPU按照預(yù)設(shè)梯度線圈參數(shù)選擇目標(biāo)場(chǎng)點(diǎn)，對(duì)GPU進(jìn)行顯存容量分配并復(fù)制原始數(shù)據(jù)，調(diào)用Matrix函數(shù)執(zhí)行流函數(shù)和繞線模式運(yùn)算處理，循環(huán)進(jìn)行上述運(yùn)算流程直至結(jié)束。將GPU運(yùn)算所得結(jié)果數(shù)據(jù)復(fù)制到計(jì)算機(jī)內(nèi)存并保存，釋放顯存中的數(shù)據(jù)[24]。

4 實(shí)現(xiàn)梯度線圈設(shè)計(jì)算法

基于梯度電感約束和GPU加速算法的梯度線圈的算法實(shí)現(xiàn)見圖2。算法流程：①設(shè)定真空磁導(dǎo)率、梯度線圈平面的最小半徑等參數(shù)值；②設(shè)定梯度線圈平面的最大直徑、線圈兩平面間距、靈敏區(qū)直徑、目標(biāo)場(chǎng)點(diǎn)個(gè)數(shù)、流函數(shù)展開級(jí)數(shù)、線圈匝數(shù)和預(yù)期梯度強(qiáng)度等參數(shù)；③選取目標(biāo)場(chǎng)點(diǎn)；④給定極坐標(biāo)下電流密度函數(shù)和磁場(chǎng)強(qiáng)度的函數(shù)表達(dá)式；⑤GPU調(diào)用矩陣函數(shù)加速運(yùn)算后，將計(jì)算結(jié)果傳回CPU，得到系數(shù)矩陣的求解結(jié)果；⑥在梯度電感最小化約束條件下計(jì)算電流密度函數(shù)的系數(shù)；⑦比較實(shí)際磁場(chǎng)強(qiáng)度與預(yù)期梯度強(qiáng)度，計(jì)算目標(biāo)場(chǎng)點(diǎn)的梯度均勻性，若在允許范圍5%以內(nèi)跳轉(zhuǎn)到流程⑧，否則跳轉(zhuǎn)至流程③，重新選取目標(biāo)場(chǎng)點(diǎn)；⑧離散化處理流函數(shù)，得到線圈繞線模型；⑨結(jié)束計(jì)算[25]。

由于梯度線圈的空間對(duì)稱性，流程③選取目標(biāo)場(chǎng)點(diǎn)時(shí)，應(yīng)將目標(biāo)場(chǎng)點(diǎn)設(shè)置在第1象限，且所取任意3個(gè)點(diǎn)不可在同一平面上，否則運(yùn)算電流密度系數(shù)矩陣時(shí)將出現(xiàn)無解。

在梯度線圈目標(biāo)場(chǎng)點(diǎn)法電流密度數(shù)學(xué)模型基礎(chǔ)上，加入梯度電感最小化和自屏蔽線圈約束條件，采用GPU加速算法實(shí)現(xiàn)梯度線圈的流函數(shù)和計(jì)算繞線分布。利用數(shù)值計(jì)算工具M(jìn)ATLAB軟件，實(shí)現(xiàn)編譯算法梯度線圈參數(shù)化設(shè)計(jì)的圖形用戶界面(graphical user interface， GUI)如圖3。

運(yùn)用上述基于GUI的梯度線圈參數(shù)化設(shè)計(jì)軟件得到梯度線圈繞線模型，通過電磁仿真軟件ANSOFT正演計(jì)算梯度線圈各項(xiàng)性能指標(biāo)，并對(duì)線圈邊界進(jìn)行修正，得到實(shí)際應(yīng)用意義上的線圈繞線分布(圖4)。

圖2 基于梯度電感約束和GPU加速算法的梯度線圈算法實(shí)現(xiàn)流程圖

5 實(shí)驗(yàn)與討論

梯度電感是渦流產(chǎn)生的施予方。梯度線圈工作電流為脈沖電流，為使其上升前沿陡峭、加快MR成像速度，在滿足其他性能指標(biāo)的前提下，梯度電感值要盡可能小。根據(jù)梯度線圈的品質(zhì)因素公式(4)，梯度線圈品質(zhì)因數(shù)與梯度線圈效率、梯度電感和梯度均勻性之間存在一定函數(shù)關(guān)系，故以梯度線圈品質(zhì)因素作為評(píng)價(jià)其性能的依據(jù)。理想梯度線圈要求梯度線圈效率盡可能高，梯度電感盡可能小，梯度磁場(chǎng)均勻性較好。事實(shí)上，梯度線圈繞線模式的關(guān)鍵在于選擇梯度線圈品質(zhì)因素Q值，Q值不同，得到的梯度繞線模式不同。以不同梯度線圈品質(zhì)因素Q模擬仿真，對(duì)測(cè)得的梯度非線性度、梯度切換率、梯度電感和梯度電路電阻相關(guān)數(shù)據(jù)進(jìn)行分析，發(fā)現(xiàn)梯度線圈品質(zhì)因素Q與線圈性能指標(biāo)的關(guān)系如圖5。

隨著梯度線圈品質(zhì)因素Q值增大，梯度線性度指標(biāo)改善，梯度切換率變小，但線圈能量損耗增大、系統(tǒng)效率降低，使梯度線圈整體性能變差。因此，適當(dāng)選擇Q值對(duì)梯度線圈性能指標(biāo)而言非常重要。

6 結(jié)論

梯度電感引起的渦流導(dǎo)致梯度磁場(chǎng)性能指標(biāo)變差，進(jìn)而影響MR成像速度和圖像質(zhì)量。為緩解渦流對(duì)MRI系統(tǒng)的不利因素，本研究基于梯度線圈目標(biāo)場(chǎng)點(diǎn)法，建立梯度電感最小和梯度磁場(chǎng)自屏蔽等約束條件下的數(shù)學(xué)模型，將電磁場(chǎng)計(jì)算的物理反問題轉(zhuǎn)換為約束條件下數(shù)學(xué)模型的優(yōu)化數(shù)值計(jì)算問題，采用GPU并行運(yùn)算加速算法，縮短了系數(shù)矩陣的計(jì)算時(shí)間。梯度線圈性能分析結(jié)果表明，采用最小化梯度電感能有效降低渦流影響，實(shí)現(xiàn)磁場(chǎng)自屏蔽，具有良好的梯度切換率，改善了梯度線圈整體性能。

圖3 梯度線圈參數(shù)化設(shè)計(jì)的GUI A.縱向梯度線圈設(shè)計(jì)的GUI； B.橫向梯度線圈設(shè)計(jì)的GUI

圖5 梯度線圈品質(zhì)因素Q與線圈性能指標(biāo)的關(guān)系 A.線圈品質(zhì)因素Q與梯度非線性度的關(guān)系； B.線圈品質(zhì)因素Q與梯度切換率的關(guān)系； C.線圈品質(zhì)因素Q與梯度電感的關(guān)系； D.線圈品質(zhì)因素Q與回路電阻的關(guān)系