基于COMSOL的磁共振成像雙平面梯度線圈的仿真研究

竇建輝，朱建明

(中國計(jì)量大學(xué)信息工程學(xué)院，杭州 310018)

0 引言

磁共振成像系統(tǒng) (magnetic resonance imaging，MRI)作為目前發(fā)展最快的成像技術(shù)，因其軟組織成像性能優(yōu)良，無輻射等優(yōu)點(diǎn)得到廣泛應(yīng)用，對其成像質(zhì)量和成像速度都提出了更高的要求。梯度線圈作為磁共振成像系統(tǒng)的核心部件之一，其主要作用是產(chǎn)生3個(gè)相互正交的梯度磁場對MRI信號的進(jìn)行空間定位編碼，即層面選擇、相位編碼和頻率編碼，因此才能實(shí)現(xiàn)磁共振圖像的重建［1］。隨著磁共振成像技術(shù)的飛速發(fā)展，圖像質(zhì)量要求越來越高，因此要求梯度線圈的梯度場場強(qiáng)更高、線性度更優(yōu)及切換速率更快。梯度線圈的性能在很大程度上決定了成像的速度及分辨率，因此設(shè)計(jì)性能更好的梯度線圈對提高磁共振成像系統(tǒng)的性能具有重要意義。

梯度線圈的設(shè)計(jì)屬于電磁場逆問題［2］，按求解方式可以分為兩類:1)直接求解，即分離導(dǎo)線法，首先選擇基本合適的線圈幾何形狀，將其產(chǎn)生的磁場按級數(shù)展開，然后根據(jù)可獲得最佳線性梯度原則來優(yōu)化線圈的位置，得到近似最優(yōu)解;2)逆推法，即目標(biāo)場方法，首先確定目標(biāo)區(qū)域內(nèi)的磁場強(qiáng)度與分布，根據(jù)畢奧－薩伐爾定律來計(jì)算產(chǎn)生此磁場的表面電流密度分布，然后用分離導(dǎo)線法近似電流密度分布。目標(biāo)場方法［3］是 R．Turner在1986年首先提出來的，并且實(shí)現(xiàn)了電感最小化［4］。隨著目標(biāo)場方法的進(jìn)一步發(fā)展，可以綜合考慮線圈儲能、線圈電感和自屏蔽線圈等多目標(biāo)約束條件，從而設(shè)計(jì)出結(jié)構(gòu)更加合理、性能更優(yōu)的梯度線圈。

目前梯度線圈的加工制作方法主要有矩形銅線繞制、澆筑金屬和切割金屬板3種［5］，之前由于加工技術(shù)水平的限制，分離導(dǎo)線的設(shè)計(jì)方法應(yīng)用更加廣泛。之前目標(biāo)場方法所設(shè)計(jì)的梯度線圈，在一定條件下，雖然線圈綜合性能得到優(yōu)化，但是其加工復(fù)雜度也隨之增大，限制了其實(shí)際的應(yīng)用范圍。但隨著加工技術(shù)水平的發(fā)展，如3D打印技術(shù)的應(yīng)用［6］，為復(fù)雜梯度線圈的加工制作提供了依據(jù)。

對梯度線圈進(jìn)行三維模型的有限元分析，可以更加接近實(shí)際情況。因此本研究采用當(dāng)前應(yīng)用非常廣泛的有限元仿真軟件COMSOLMultiphysics進(jìn)行分析研究，通過有限元計(jì)算得到電磁場分布，計(jì)算精度和速度可有效提高［7］。針對開放式磁共振系統(tǒng)的磁體結(jié)構(gòu)，進(jìn)行x線圈的設(shè)計(jì)及仿真，y線圈的設(shè)計(jì)與x線圈基本相同，x線圈旋轉(zhuǎn)90°即可得到y(tǒng)線圈。Z線圈的設(shè)計(jì)較為成熟，一般為一組規(guī)則的同心圓，設(shè)計(jì)方法與橫向線圈不同。本研究首先對改進(jìn)的目標(biāo)場設(shè)計(jì)方法進(jìn)行理論介紹，考慮線圈儲能及線圈能耗為約束進(jìn)行梯度線圈電流密度分布的計(jì)算，然后通過AUTOCAD建立三維模型，最后采用COMSOL有限元軟件對模型進(jìn)行仿真研究，分析了不同設(shè)計(jì)參數(shù)下梯度線圈的性能，為實(shí)際的梯度線圈設(shè)計(jì)提供借鑒指導(dǎo)。

1 目標(biāo)場方法

磁共振成像橫向雙平面梯度線圈模型如圖1所示。首先整個(gè)模型處于主磁場大小為B0的區(qū)域內(nèi)。橫向雙平面梯度線圈分別位于z=±L的兩平面上，由于線圈厚度與線圈半徑相比很小，因此可以將電流密度近似限制在圓形薄片區(qū)域內(nèi)，即電流密度分布在兩平面的直徑為p0≤p≤pm的圓形區(qū)域內(nèi)。以坐標(biāo)原點(diǎn)為中心的球?yàn)樘荻却艌龅母信d趣區(qū)域(region of interest，ROI)，最終求得球面區(qū)域內(nèi)為均勻梯度場。P為球面任意一點(diǎn)。以感興趣區(qū)域的圓心為坐標(biāo)原點(diǎn)建立三維坐標(biāo)系。假設(shè)面電流密度函數(shù)為J，在柱面坐標(biāo)系下，將其展開為三角函數(shù)形式為:

圖1 磁共振雙平面梯度線圈模型示意圖

式中，c=π/(pm－p0)。Un為第n階電流密度分量的權(quán)重系數(shù)，N為基函數(shù)的所有分量個(gè)數(shù)。理論上N的個(gè)數(shù)為無窮，但實(shí)際上限于計(jì)算精度及加工技術(shù)的限制，根據(jù)設(shè)計(jì)結(jié)果選擇一個(gè)足夠大且有限的數(shù)字來近似電流密度。

1.1 磁場強(qiáng)度

根據(jù)Biot－Savart定律，雙平面梯度線圈的電流密度函數(shù)珋J±在空間任意一點(diǎn)P(p1，φ1，z1)產(chǎn)生的沿Z軸磁感應(yīng)強(qiáng)度分量為［8］:

其中:μ0為真空磁導(dǎo)率，R±分別為上下極板中線圈源點(diǎn)到目標(biāo)場點(diǎn)的距離，

當(dāng)已知目標(biāo)區(qū)域場點(diǎn)坐標(biāo)時(shí)，可以將磁場z軸分量Bz(x，y，z)寫成矩陣形式:

其中:Bz為M個(gè)目標(biāo)區(qū)域的場點(diǎn)沿Z軸場強(qiáng)分量的列向量，MT可根據(jù)所設(shè)置的坐標(biāo)由公式 (3)直接計(jì)算出來，然后通過矩陣運(yùn)算求解出U中的N個(gè)基函數(shù)分量的未知權(quán)重系數(shù)Un，即求得近似的電流密度。

1.2 建立目標(biāo)函數(shù)

電感作為衡量梯度線圈性能好壞的一個(gè)重要指標(biāo)，對成像速度有著至關(guān)重要的影響。要得到更快的開關(guān)速度，必須使得線圈電感盡可能減小。根據(jù)基爾霍夫定律可知，為了減小電流的時(shí)間常數(shù)，需要低電感和高電阻。顯然電阻的增高受到梯度電源功率和散熱效率的限制，不可能無限制增大，因此需要盡可能的獲取最小電感。為了優(yōu)化提升梯度線圈的綜合性能，發(fā)揮目標(biāo)場設(shè)計(jì)方法的特點(diǎn)，以梯度磁場的非線性度、線圈電感儲能和線圈能耗為約束函數(shù)，建立目標(biāo)函數(shù)。

線圈儲能:

線圈功率:

式中，δ為梯度線圈的電阻系數(shù)，t為線圈厚度，S表示電流密度在平面上的分布區(qū)域。

求解電流密度分量的權(quán)重系數(shù)是一個(gè)經(jīng)典的問題。以所有設(shè)置的目標(biāo)場點(diǎn)的z軸磁場強(qiáng)度誤差平方和最小為目標(biāo)，同時(shí)考慮線圈儲能和線圈能耗的多因素約束，構(gòu)建目標(biāo)函數(shù):

式中，Bz，n表示仿真實(shí)驗(yàn)中得到的第n個(gè)場點(diǎn)的Z軸磁場強(qiáng)度分量，Bzdes，n表示理想設(shè)計(jì)中第n個(gè)場點(diǎn)的Z軸磁場強(qiáng)度分量。λ1、λ2分別表示線圈儲能和線圈能耗的權(quán)重系數(shù)。

1.3 電流密度離散繞線分布

電流密度函數(shù)滿足連續(xù)性方程，柱面坐標(biāo)下流函數(shù)ψ滿足如下條件:

將電流密度函數(shù)帶入得:

最后將流函數(shù)離散化得到梯度線圈的繞線分布。由于線圈是對稱的，設(shè)半圓平面上的線圈匝數(shù)為Nt，得到離散公式有:

其中，ψmin、ψmax分別為流函數(shù)在半平面區(qū)域內(nèi)的最小值和最大值。

2 三維仿真研究

COMSOL多物理場仿真軟件通過有限元方法模擬在科研和工程中能用偏微分方程 (PDE)描述的各種問題，用戶只需選擇或者自定義不同專業(yè)的偏微分方程進(jìn)行任意組合便可輕松實(shí)現(xiàn)多物理場的直接耦合分析，并提供了友好的圖形用戶界面，全面的第三方CAD導(dǎo)入功能，支持當(dāng)前主流CAD軟件格式文件的導(dǎo)入，同時(shí)擁有高效的計(jì)算性能和強(qiáng)大的多場雙向直接耦合能力，能夠?qū)崿F(xiàn)高精度的數(shù)值仿真［9］。

本文利用此有限元軟件分析梯度線圈施加電流后的磁場分布情況，線圈內(nèi)的電流密度分布情況。在COMSOL軟件中，通過調(diào)整和修改模型參數(shù)，材料屬性和施加電流的大小和方向，可以得到梯度線圈產(chǎn)生的不同磁場分布和磁場強(qiáng)度的大小。通過后處理，來分析和評價(jià)梯度線圈的綜合性能，然后通過模擬仿真的結(jié)果來調(diào)整梯度線圈的設(shè)計(jì)參數(shù)，以使得設(shè)計(jì)結(jié)果更加符合實(shí)際成像速度和質(zhì)量的要求。圖2為本文的整體設(shè)計(jì)分析流程圖。

結(jié)合COMSOL的強(qiáng)大計(jì)算分析能力，對所設(shè)計(jì)的梯度線圈模型進(jìn)行模擬仿真，來檢驗(yàn)設(shè)計(jì)結(jié)果的優(yōu)劣。

圖2 整體設(shè)計(jì)流程圖

2.1 導(dǎo)入模型

通過MATLAB計(jì)算線圈的繞線分布情況，得到梯度線圈離散的坐標(biāo)數(shù)據(jù)。采用AutoCAD建立梯度線圈的3D模型，最后將模型導(dǎo)入COMSOL軟件中進(jìn)行模型的檢測與調(diào)整，然后進(jìn)行仿真計(jì)算分析，圖3為雙平面x梯度線圈模型圖。

如圖3所示，雙平面x梯度線圈的材料為圓形銅線，也可選擇方形材料，線圈的上下為兩塊硅鋼極板，用來屏蔽線圈與主磁體，減小產(chǎn)生的渦流，中心為球形目標(biāo)區(qū)域。因?yàn)榍蠼鈺r(shí)還需添加求解區(qū)域，因此最后還要添加圓柱體區(qū)域包括整個(gè)模型。

圖3 橫向雙平面x梯度線圈三維模型圖

2.2 定義模型材料屬性

選擇COMSOL自帶材料庫中的材料，設(shè)置上下極板材料為硅鋼片，線圈材料為銅，中心球面目標(biāo)區(qū)域和整個(gè)圓柱求解域材料為真空。

2.3 添加物理場

添加所要研究的物理場，選擇磁場 (mf)。添加線圈電流激勵(lì)，設(shè)置導(dǎo)體模型為單匝線圈，線圈激勵(lì)類型為電流，最后在幾何分析中添加輸入，設(shè)置電流大小，并選擇電流方向，對稱雙平面梯度線圈的電流方向相反，因?yàn)榫€圈匝數(shù)較多，在添加電流激勵(lì)之前可以先進(jìn)行顯示定義，將線圈組定義為顯示線圈，這樣再添加電流激勵(lì)就可以了。

2.4 劃分網(wǎng)格

COMSOL軟件有物理場控制網(wǎng)格和用戶控制網(wǎng)格兩種模式。物理場控制網(wǎng)格是選擇自由四面體網(wǎng)格類型進(jìn)行剖分，網(wǎng)格單元質(zhì)量根據(jù)模型大小分為9個(gè)等級，用戶可自由選擇。用戶控制網(wǎng)格可以選擇網(wǎng)格的形狀，網(wǎng)格劃分方式，以及自定義網(wǎng)格單元的大小范圍。本文模型主要為圓柱線圈，形狀較為規(guī)則，因此采用物理場控制網(wǎng)格模式，選擇較細(xì)化的單元質(zhì)量進(jìn)行自由四面體類型的網(wǎng)格劃分，提高實(shí)驗(yàn)結(jié)果的精確度。

2.5 研究

添加研究進(jìn)行求解。因?yàn)檎w線圈為不規(guī)則圓形形狀，首先需要對線圈的幾何形狀進(jìn)行分析，以獲取電流在線圈中的電流分布狀況。所以研究分為兩步:首先要添加研究步驟1，添加線圈幾何分析研究，得到電流在線圈中的電流密度分布;然后選擇步驟2，添加穩(wěn)態(tài)研究，得到線圈的磁通密度分布。

2.6 后處理

研究計(jì)算完畢后得到線圈中的電流密度分布、磁場分布以及線圈電感等數(shù)據(jù)，COMSOL在計(jì)算過程中能夠?qū)λ袛?shù)據(jù)進(jìn)行保存，以便我們對其進(jìn)行后處理，得到表格或圖表等性能分析結(jié)果，如圖4。

圖4 磁通密度切面及磁感應(yīng)強(qiáng)度分布圖

3 仿真實(shí)驗(yàn)結(jié)果分析

本文所設(shè)計(jì)的橫向x梯度線圈，線圈為直徑為0.2 cm的圓形銅線，平面線圈半徑大小為42 cm，線圈中心留有直徑為1 cm的圓形孔洞，用來固定線圈，目標(biāo)區(qū)域半徑大小為15 cm，線圈距目標(biāo)區(qū)域圓心距離為23 cm，求解域?yàn)殚L60 cm，底面直徑80 cm的圓柱體，電流激勵(lì)大小為500 A，所設(shè)置目標(biāo)場點(diǎn)為在目標(biāo)球面區(qū)域上，均勻放置，放置個(gè)數(shù)為18，并計(jì)算其理想磁場強(qiáng)度z軸分量。主要仿真結(jié)果如圖4，觀察其磁通密度切面圖及磁感應(yīng)強(qiáng)度在xz平面上的z軸分量圖。

梯度線圈評價(jià)指標(biāo)主要采用梯度場的場強(qiáng)、非線性度、電感和電阻，綜合考慮4種指標(biāo)，選擇合適的設(shè)計(jì)方案。梯度場場強(qiáng)決定了空間分辨率，場強(qiáng)越高，圖像的分辨率越高。通常要求梯度場的非線性度小于5%，線性度越好，空間定位越準(zhǔn)確，圖像的質(zhì)量也就越好。對于快速成像和實(shí)時(shí)成像，線圈電感和電阻是影響其質(zhì)量的重要因素，成像速度越快，要求梯度線圈的開關(guān)速度也就越快，因此就要求線圈的電感越小，以便梯度場的快速產(chǎn)生和消失。由仿真結(jié)果計(jì)算所有線圈的自感與互感，得到線圈電感和線圈間的耦合關(guān)系，來評價(jià)設(shè)計(jì)結(jié)果是否符合要求。

在第一象限區(qū)域內(nèi)設(shè)置與x軸平行的16根直線，作為磁場的檢測線，對其進(jìn)行均勻采樣，由公式 (9)計(jì)算梯度場的非線性度。

其中，Bz(des，i)為理想設(shè)計(jì)情況下目標(biāo)區(qū)域點(diǎn)的Z軸磁場強(qiáng)度，Bz(cal，i)是根據(jù)畢奧薩伐爾定律計(jì)算得到的Z軸磁場強(qiáng)度。P為計(jì)算非線性度所取總的采樣點(diǎn)的個(gè)數(shù)。

圖5為磁感應(yīng)強(qiáng)度z軸分量在x方向的變化曲線，根據(jù)所設(shè)置的16根采樣線對目標(biāo)區(qū)域內(nèi)磁場強(qiáng)度z軸分量進(jìn)行采樣，繪制曲線圖。由結(jié)果可知，目標(biāo)區(qū)域中心位置梯度場較為均勻，線性度較好，符合設(shè)計(jì)預(yù)期，在邊緣處開始產(chǎn)生較大的震蕩。

圖5 磁感應(yīng)強(qiáng)度z軸分量

表1為在選擇不同線圈匝數(shù)時(shí)，梯度線圈的綜合性能評價(jià)?？梢钥闯?，隨著線圈匝數(shù)的變化，梯度場的非線性度、線圈儲能和線圈能耗的變化并非線性，需根據(jù)實(shí)際要求選擇合適的匝數(shù)。當(dāng)線圈匝數(shù)為13和14時(shí)，梯度場的非線性度才滿足小于5%的應(yīng)用要求，因此線圈匝數(shù)需從二者中進(jìn)行選擇。此時(shí)，線圈電阻和電感都能夠滿足需求，如果選擇較大的線圈匝數(shù)，梯度場場強(qiáng)較高，而梯度場的非線性度、電感和電阻均增加，然而其場強(qiáng)的增加程度是否優(yōu)于其非線性度、電感和電阻所帶來的不利影響。因此需要根據(jù)實(shí)際要求來綜合考慮設(shè)計(jì)參數(shù)。

表1 不同線圈匝數(shù)下的設(shè)計(jì)結(jié)果_

表2為在確定線圈匝數(shù)為13后，梯度線圈在不同基函數(shù)分量個(gè)數(shù)情況下的性能。由表2可以看出，當(dāng)基函數(shù)分量個(gè)數(shù)增加時(shí)，線圈的磁場梯度、電感及電阻變化沒有線圈匝數(shù)增加時(shí)變化的明顯，其非線性度減小，說明其設(shè)計(jì)結(jié)果對于圖像質(zhì)量有所改善。但是在使用MATLAB計(jì)算其線圈分布數(shù)據(jù)時(shí)，計(jì)算時(shí)間明顯延長。通過表1和表2可以看出，磁場梯度的增加主要與梯度線圈匝數(shù)變化有關(guān)，同時(shí)由于線圈匝數(shù)的增加，使得線圈電阻和電感都會增加，而梯度場非線性度的變化與線圈匝數(shù)及基函數(shù)分量個(gè)數(shù)都有關(guān)系。

表2 不同基函數(shù)分量個(gè)數(shù)下的設(shè)計(jì)結(jié)果

基于圖像質(zhì)量應(yīng)該作為設(shè)計(jì)結(jié)果的主要要求，因此首先應(yīng)該考慮線圈非線性度，作為篩選條件。實(shí)驗(yàn)表明，由目標(biāo)場方法設(shè)計(jì)，采用流函數(shù)方法對線圈分布進(jìn)行離散時(shí)，當(dāng)線圈匝數(shù)超過15匝時(shí)或者當(dāng)基函數(shù)分量個(gè)數(shù)超過5，如圖6，梯度線圈的繞線開始產(chǎn)生震蕩，之前限于加工技術(shù)的水平，震蕩小線圈的出現(xiàn)使得加工的難度大大增加，現(xiàn)在利用3D打印技術(shù)的優(yōu)勢［10］，可以進(jìn)行加工制作，這樣就能得到更高的場強(qiáng)。

圖6 梯度線圈繞線分布

根據(jù)實(shí)驗(yàn)結(jié)果，綜合考慮線圈儲能和線圈能耗的因素，當(dāng)線圈匝數(shù)為13匝，基函數(shù)分量個(gè)數(shù)為3時(shí)，通過仿真數(shù)據(jù)后處理計(jì)算得線圈電阻為0.159 4 Ω，電感為0.066 27 mH，磁場梯度為26 mT/m。

4 結(jié)論

對于本文采用目標(biāo)場方法設(shè)計(jì)的梯度線圈，由于將電流密度函數(shù)展開為三角函數(shù)，并且采用流函數(shù)方法進(jìn)行線圈離散化得到繞線分布，其電流密度基函數(shù)分量個(gè)數(shù)和線圈匝數(shù)以及線圈儲能和線圈功率，都會對梯度線圈的性能產(chǎn)生影響，需要通過前期仿真實(shí)驗(yàn)綜合考慮多種因素來確定設(shè)計(jì)參數(shù)的最優(yōu)值，滿足實(shí)際需求。

通過以上的仿真分析結(jié)果可知，所設(shè)計(jì)的梯度線圈產(chǎn)生的磁感應(yīng)強(qiáng)度在DSV內(nèi)具有良好的線性和均勻性，可滿足實(shí)際設(shè)計(jì)要求，對梯度線圈的設(shè)計(jì)優(yōu)化具有較好的理論和實(shí)踐指導(dǎo)意義。通過適當(dāng)修改相關(guān)參數(shù)，可以得到適應(yīng)不同應(yīng)用要求的梯度線圈模型，為梯度線圈的應(yīng)用提供依據(jù)。結(jié)合3D打印技術(shù)，進(jìn)行梯度線圈的制作。

渦流問題對成像質(zhì)量也有較大的影響，本文采用增加抗渦流板的無緣屏蔽方式，可在一定程度上減小渦流。但是抗渦流板不能完全將線圈和主磁體阻隔，仍然會對成像帶來影響。另一種解決方式是有緣屏蔽，采用自屏蔽梯度線圈，即在梯度線圈的外面添加一組電流方向與其相反的線圈，使得屏蔽區(qū)域以外的磁場強(qiáng)度為零。這種設(shè)計(jì)方式將在接下來的研究中進(jìn)行實(shí)驗(yàn)，從而設(shè)計(jì)出更好的梯度線圈系統(tǒng)。