基于梯度電感約束和圖形處理單元加速算法的磁共振雙平面自屏蔽梯度線圈設(shè)計(jì)方法

黃清明,鄭 剛

(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)系