基于伽遼金有限元法的磁感應(yīng)斷層成像正問(wèn)題仿真

柯 麗 龐 佩 佩 杜 強(qiáng)

（沈陽(yáng)工業(yè)大學(xué)生物醫(yī)學(xué)與電磁工程研究所，沈陽(yáng) 110870）

KE Li* PANG Pei-Pei DU Qiang

（Institute of Biomedical and Electromagnetic Engineering‘Shenyang University of Technology’Shenyang110870，China）

引言

磁感應(yīng)斷層成像（MIT）技術(shù)是一種利用電磁檢測(cè)原理測(cè)量生物組織電導(dǎo)率的技術(shù)［1－3］，該技術(shù)根據(jù)渦流感應(yīng)原理成像，在不接觸被測(cè)物的情況下，利用檢測(cè)線圈感應(yīng)成像區(qū)域內(nèi)的渦流磁場(chǎng)，并通過(guò)測(cè)量檢測(cè)線圈的相位差重構(gòu)出成像區(qū)域內(nèi)電導(dǎo)率的分布，進(jìn)而實(shí)現(xiàn)MIT的圖像重建。與現(xiàn)有的成像技術(shù)相比，MIT設(shè)備造價(jià)和檢測(cè)成本較為低廉，對(duì)人體安全、無(wú)創(chuàng)、無(wú)副作用，并可以功能成像，非常適用于長(zhǎng)時(shí)間的圖像監(jiān)護(hù)［4－5］。因此，MIT在人體生理功能、疾病診斷以及臨床監(jiān)護(hù)等方面具有重要的醫(yī)學(xué)價(jià)值。

磁感應(yīng)斷層成像正問(wèn)題是［6－7］在已知激勵(lì)源和場(chǎng)域內(nèi)電導(dǎo)率分布的情況下，研究場(chǎng)域內(nèi)渦流分布特性和相位延遲特性。正問(wèn)題的求解將為硬件系統(tǒng)中線圈的尺寸設(shè)計(jì)、頻率選擇以及測(cè)量數(shù)據(jù)正確性的判斷提供理論參考，并能為重建算法提供仿真數(shù)據(jù)。目前，國(guó)內(nèi)外有很多學(xué)者針對(duì)MIT正問(wèn)題進(jìn)行了研究，Scharfetter等使用軟件LC（基于時(shí)域有限差分法（FDTD））求解Maxwell方程［8］；Merwa等應(yīng)用軟件包HyperMesh采用有限元方法研究了磁感應(yīng)成像中的渦流問(wèn)題［9］；劉國(guó)強(qiáng)等在ANSYS中采用棱單元法求解磁感應(yīng)成像的正問(wèn)題［10］；董秀珍等采用里茲變分法推導(dǎo)了二維渦流問(wèn)題［11］。總之，目前對(duì)于正問(wèn)題的研究基本是利用商用軟件包進(jìn)行求解，而在MIT圖像重建算法的過(guò)程中，需要反復(fù)不斷的求解正問(wèn)題，占用了絕大部分的計(jì)算工作量。但MIT應(yīng)用于臨床醫(yī)學(xué)中的關(guān)鍵技術(shù)是人體二維層析成像的實(shí)時(shí)顯示，因此商用軟件包不適用于這種場(chǎng)合。

對(duì)于復(fù)雜、不規(guī)則的模型，正問(wèn)題不存在解析解，因此需要采用有限元方法求解。對(duì)于MIT成像來(lái)說(shuō)，由于體內(nèi)電導(dǎo)率分布的不均勻，在激勵(lì)場(chǎng)的作用下將會(huì)感應(yīng)渦流，即產(chǎn)生渦流算子，而渦流算子是非正定算子［12］，基于泛函理論的變分有限元要求微分算子是正定算子，因此無(wú)法采用該法精確求解。與變分有限元法相比，伽遼金有限元法可用來(lái)求解更為廣泛的電磁場(chǎng)問(wèn)題，本研究采用該方法求出了MIT正問(wèn)題的精確解并給出了相應(yīng)的仿真結(jié)果，可為MIT硬件系統(tǒng)的測(cè)量及重建算法的研究提供實(shí)驗(yàn)參考和理論依據(jù)。

1 MIT正問(wèn)題的數(shù)學(xué)模型

MIT正問(wèn)題的本質(zhì)是求解準(zhǔn)靜態(tài)時(shí)諧渦流場(chǎng)［6］的分布問(wèn)題，如圖1所示，二維正問(wèn)題可以描述如下：在激勵(lì)線圈中通入激勵(lì)頻率為ω、面密度為J的交變電流，與其有限距離內(nèi)，有一電導(dǎo)率為σ、磁導(dǎo)率為μ的目標(biāo)導(dǎo)體，求解此時(shí)導(dǎo)體與周圍空間中的磁場(chǎng)分布以及檢測(cè)線圈中的相位差。

圖1 MIT正問(wèn)題模型Fig.1 MIT forward model

根據(jù)Maxwell理論，MIT正問(wèn)題的微分方程：

邊界條件：n×H0=n×H1

式中，A為求解區(qū)域的矢量磁位，H0和H1為邊界分界處的磁場(chǎng)強(qiáng)度，這樣，MIT正問(wèn)題轉(zhuǎn)化成求解矢量磁位的非線性二階微分方程問(wèn)題。

2 基于伽遼金有限元法的正問(wèn)題求解

體內(nèi)組織各個(gè)斷層部位的電導(dǎo)率分布是不規(guī)則的，上述方程根本不存在解析解，再加上方程包含渦流算子jωσA，因此，本研究采用伽遼金有限元法對(duì)該方程進(jìn)行求解。

2.1 伽遼金有限元法的求解過(guò)程

伽遼金有限元法［13］是加權(quán)余量法的一種，基本思想是離散區(qū)域內(nèi)所有單元的殘數(shù)加權(quán)積分為零。首先，把成像區(qū)域進(jìn)行有限元?jiǎng)澐?，離散成M個(gè)三角單元，每個(gè)單元e的殘數(shù)加權(quán)

令

對(duì)于伽遼金有限元法，三角單元對(duì)應(yīng)的權(quán)函數(shù)等于其相應(yīng)的插值基函數(shù)，因此

式中，Δ為三角單元的面積，ai、bi、ci為與單元坐標(biāo)相關(guān)的系數(shù)。

將式（6）代入式（5）

式中，K包含成像區(qū)域的單元結(jié)構(gòu)信息及導(dǎo)體電導(dǎo)率的分布信息，F(xiàn)是右端向量，由電流源決定。

由上述過(guò)程得到M個(gè)線性方程組，應(yīng)用高斯消去法求解該方程組，得到成像區(qū)域內(nèi)矢量磁位A的分布。

2.2 磁感應(yīng)強(qiáng)度和渦流場(chǎng)的計(jì)算

根據(jù)場(chǎng)域內(nèi)矢量磁位A的分布，求解成像區(qū)域內(nèi)磁感應(yīng)強(qiáng)度和渦流密度的分布：每個(gè)三角單元的磁感應(yīng)強(qiáng)度：

式中，i，j=1、2、3

對(duì)所有單元組合形成總體剛度矩陣

對(duì)應(yīng)的渦流強(qiáng)度［14］：

2.3 相位差的計(jì)算

由MIT理論分析可知，成像區(qū)域內(nèi)電導(dǎo)率的分布與檢測(cè)線圈中的相位差相關(guān)［1］，通過(guò)以下推導(dǎo)過(guò)程，進(jìn)一步驗(yàn)證該關(guān)系。

成像區(qū)域內(nèi)電導(dǎo)率分布均勻時(shí)，激勵(lì)線圈在檢測(cè)線圈中產(chǎn)生的磁感應(yīng)強(qiáng)度B值大小

電導(dǎo)率分布不均勻時(shí)，在成像區(qū)域內(nèi)將感應(yīng)渦流，此時(shí)檢測(cè)線圈中磁場(chǎng)的變化量

式中，R、a、t、m是與模型相關(guān)的參數(shù)。

檢測(cè)線圈上感應(yīng)電流的相位差［15］

由于ωε0εrξ＜＜1

通過(guò)上式可以看出，檢測(cè)線圈的相位差和目標(biāo)導(dǎo)體的電導(dǎo)率與激勵(lì)源頻率的乘積成線性關(guān)系，對(duì)于確定的成像模型來(lái)說(shuō)，該比值是定值。

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

基于上述理論推導(dǎo)過(guò)程，求解MIT正問(wèn)題并進(jìn)行相應(yīng)的仿真實(shí)驗(yàn)，參數(shù)設(shè)置如下：成像模型直徑200 mm，區(qū)域剖分6層，最外層為遠(yuǎn)場(chǎng)區(qū)，求解區(qū)域4層，40個(gè)線圈均勻分布在圓周上，采用激勵(lì)—檢測(cè)線圈模式，激勵(lì)線圈位于成像區(qū)域最右側(cè)，激勵(lì)線圈中通入10 MHz、750 mA的交流電。

3.1 不同成像模型對(duì)磁場(chǎng)分布的影響

由于激勵(lì)源與目標(biāo)導(dǎo)體的相互作用，不同的成像模型將影響磁場(chǎng)的分布規(guī)律。本研究在保證其它參數(shù)不變的前提下，通過(guò)改變目標(biāo)導(dǎo)體的位置和數(shù)目分別觀察均勻場(chǎng)及擾動(dòng)場(chǎng)的參數(shù)分布圖。

由渦流與磁感應(yīng)強(qiáng)度的公式推導(dǎo)過(guò)程可以看出，兩者都是矢量磁位的線性組合，但是實(shí)部與虛部相反，因此本研究對(duì)磁感應(yīng)強(qiáng)度的虛部以及渦流強(qiáng)度的實(shí)部進(jìn)行了仿真。

圖2顯示了均勻場(chǎng)的成像模型以及磁感應(yīng)強(qiáng)度虛部的分布圖（由于均勻場(chǎng)時(shí)沒(méi)有目標(biāo)導(dǎo)體，因此不會(huì)產(chǎn)生渦流），均勻場(chǎng)的電導(dǎo)率設(shè)置為0.1×10－5S/m。圖3顯示了目標(biāo)導(dǎo)體位于左側(cè)、中心以及同時(shí)位于中心和左側(cè)兩個(gè)位置時(shí)的成像模型、磁感應(yīng)強(qiáng)度虛部、渦流實(shí)部的分布圖，其中目標(biāo)導(dǎo)體的電導(dǎo)率均設(shè)置為0.4S/m。從圖中可以看出：均勻場(chǎng)時(shí)，磁感應(yīng)強(qiáng)度的虛部隨著與激勵(lì)線圈位置的增大而減小。擾動(dòng)場(chǎng)時(shí)，磁感應(yīng)強(qiáng)度的虛部在目標(biāo)導(dǎo)體位置最大并能反映出目標(biāo)導(dǎo)體的位置，但同時(shí)在目標(biāo)導(dǎo)體的周圍也會(huì)產(chǎn)生偽影；由于螺旋管電流的注入方式，渦流的實(shí)部在目標(biāo)導(dǎo)體位置發(fā)生變化該參數(shù)能準(zhǔn)確反映目標(biāo)導(dǎo)體的位置。由仿真圖可以看出，對(duì)于單目標(biāo)以及多目標(biāo)的擾動(dòng)模型，磁感應(yīng)強(qiáng)度的虛部和渦流的實(shí)部均能夠反映出成像模型中目標(biāo)導(dǎo)體的位置信息；但是，對(duì)于多目標(biāo)擾動(dòng)來(lái)說(shuō)，與激勵(lì)線圈位置接近的目標(biāo)導(dǎo)體較其它位置對(duì)磁場(chǎng)以及渦流的分布影響更大。

圖2 成像模型及磁場(chǎng)分布圖。（a）均勻場(chǎng)模型；（b）磁感應(yīng)強(qiáng)度虛部值Fig.2 Imaging model and magnetic field distribution.（a）Uniform field model；（b）Imaginary value

3.2 電導(dǎo)率大小及位置對(duì)檢測(cè)線圈中相位差的影響

圖3 不同類型擾動(dòng)場(chǎng)成像模型及場(chǎng)分布圖。（a），（d），（g）擾動(dòng)模型；（b），（e），（h）相應(yīng)的磁感應(yīng)強(qiáng)度虛部值；（c），（f），（i）相應(yīng)的渦流強(qiáng)度實(shí)部Fig.3 Perturbation model and field distribution.（a），（d），（g）Perturbation models；（b），（e），（h）The imaginary values；（c），（f），（i）Eddy current in real part

為考察電導(dǎo)率大小對(duì)檢測(cè)線圈中相位差的影響，在其它仿真參數(shù)不變的情況下，使中心區(qū)域的電導(dǎo)率由1S/m依次增大到6S/m，如圖4所示，（a）圖顯示了電導(dǎo)率為1S/m、3S/m和5S/m時(shí)，邊界檢測(cè)線圈中的相位差隨電導(dǎo)率變化的曲線，從圖中可以看出：電導(dǎo)率越大，相位差的最大值越大；與激勵(lì)線圈180°的位置，相位差最大，并關(guān)于該點(diǎn)呈對(duì)稱分布；越接近激勵(lì)線圈，相位差越小。（b）圖顯示了與激勵(lì)線圈180°、252°、162°3個(gè)位置的檢測(cè)線圈，在中心區(qū)域的電導(dǎo)率由1～6S/m依次增加時(shí)相位差的變化，從圖中可以看出：相位差在三個(gè)位置的值依次減小，即檢測(cè)線圈與激勵(lì)線圈的相對(duì)位置越遠(yuǎn)，直線的斜率越大，即電導(dǎo)率與相位差的比值越大；對(duì)于確定位置的檢測(cè)線圈，電導(dǎo)率與相位差成線性關(guān)系。再次證實(shí)了上述理論推導(dǎo)的相位差與電導(dǎo)率呈線性關(guān)系的正確性。

圖4 電導(dǎo)率大小變化時(shí)的相位差。（a）檢測(cè)線圈的相位差隨電導(dǎo)率變化圖；（b）電導(dǎo)率依次增大時(shí)檢測(cè)線圈中的相位差Fig.4 Phase shift of conductivity change.（a）Every detecting coil phase shift with conductivity change；（b）Detecting coil phase shift change with conductivity increase

圖5顯示了目標(biāo)導(dǎo)體位于不同位置時(shí)邊界檢測(cè)線圈中相位差的分布圖，其中目標(biāo)導(dǎo)體的電導(dǎo)率均設(shè)置為1S/m。從圖中可以看出：當(dāng)目標(biāo)導(dǎo)體位于左側(cè)時(shí)，與目標(biāo)導(dǎo)體最近的180°位置檢測(cè)線圈相位差最大；位于中間位置時(shí)，與所有檢測(cè)線圈距離相同，此時(shí)，與激勵(lì)線圈相對(duì)位置最遠(yuǎn)的180°位置相位差最大；同理，右側(cè)位置的目標(biāo)導(dǎo)體最大值出現(xiàn)在180°檢測(cè)線圈的位置，由于激勵(lì)線圈的干擾，0°檢測(cè)線圈相位差值不穩(wěn)定。因此，當(dāng)目標(biāo)導(dǎo)體位置變化時(shí)，檢測(cè)線圈中的相位差與目標(biāo)導(dǎo)體的位置越近，激勵(lì)線圈的位置越遠(yuǎn)，值越大。

圖5 不同位置目標(biāo)導(dǎo)體相位差的分布。（a）目標(biāo)導(dǎo)體的位置；（b）邊界相位差分布圖Fig.5 Phase shift Distribution of different object conductor.（a）Object conductor location；（b）Phase shift in edge

4 結(jié)論

有效求解MIT正問(wèn)題是進(jìn)行MIT圖像重建的基礎(chǔ)，本文針對(duì)MIT正問(wèn)題中渦流場(chǎng)微分方程的非正定性，應(yīng)用伽遼金有限元法對(duì)該微分方程進(jìn)行求解，并在此基礎(chǔ)上推導(dǎo)了目標(biāo)導(dǎo)體的電導(dǎo)率與檢測(cè)線圈中相位差的關(guān)系。然后在6層有限元模型上對(duì)不同成像模型進(jìn)行了一系列仿真，結(jié)果表明：磁感應(yīng)強(qiáng)度虛部與渦流強(qiáng)度實(shí)部的分布能夠反映出目標(biāo)導(dǎo)體的位置信息，為MIT圖像重建算法的研究提供了理論依據(jù)；通過(guò)分析不同成像模型相位差的分布規(guī)律，總結(jié)了目標(biāo)導(dǎo)體電導(dǎo)率與檢測(cè)線圈位置的關(guān)系，證明了同一位置電導(dǎo)率與檢測(cè)線圈中相位差的線性關(guān)系，為MIT硬件系統(tǒng)的實(shí)驗(yàn)測(cè)量結(jié)果提供了參考。上述的仿真實(shí)驗(yàn)符合MIT的理論分析，進(jìn)而驗(yàn)證了本研究有限元法的有效性。本研究應(yīng)用的伽遼金有限元法同樣可以推廣到三維的情況。

［1］Griffiths H.Magnetic induction tomography［J］.Meas Sci Technol，2001，12（8）：1126－1131.

［2］Johannes N，Ewald F，SabineH.Contactlessimpedance measurementby magnetic induction-a possible method for investigation of brain impedance［J］.Physiol Meas，1993，14（4）：463－471.

［3］Morris A，Giffiths H，GoughW.A numerical model for induction tomographic measurements in biological tissues ［J］.Physiol Meas，2001，22（1）：113－119.

［4］王聰，董秀珍，秦明新.電磁感應(yīng)斷面成像研究的關(guān)鍵問(wèn)題［J］.CT理論與應(yīng)用研究，2003（3）：17－21.

［5］劉國(guó)強(qiáng)，霍小林.層狀生物組織磁感應(yīng)成像［J］.中國(guó)醫(yī)學(xué)物理學(xué)雜志，2003（1）：59－61.

［6］Merwa R，Hollaus K，Brandsttter B，et al..Numerical solution of the general 3D eddy current problem for magnetic induction tomography（spectroscopy）［J］.Physiol Meas，2003，24：545－554.

［7］Griffiths H，StewartWR，Gough W.Magneticinduction tomography—A measuring system for biological tissues［J］.Ann NY Acad Sci，1999，873：335－345.

［8］Hermann S，Pere R，Marcos P，et al.Sensitivity maps for lowcontrast perturbations within conducting background in magnetic induction tomography［J］.Physiol Meas，2002，23：195－202.

［9］Robert M，Karl H，Bernhard B，et al..Numerical solution of the general 3D eddy current problem for magnetic induction tomograph（spectroscopy）［J］.Physiol Meas，2003；24：545－554.

［10］劉國(guó)強(qiáng)，王濤.用棱單元方法求解磁感應(yīng)成像的正問(wèn)題［J］.中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào)，2006，25（2）：163－166.

［11］王聰，董秀珍.磁感應(yīng)斷層成像技術(shù)中渦流問(wèn)題的有限元法仿真研究［J］.航天醫(yī)學(xué)與醫(yī)學(xué)工程，2007，20（3）：219－222.

［12］魯滌強(qiáng)，黃學(xué)良，胡敏強(qiáng).低頻渦流電磁場(chǎng)非自伴變分問(wèn)題的研究［J］.東南大學(xué)學(xué)報(bào)，2004，34（1）：1－5.

［13］謝德馨，楊仕友.工程電磁場(chǎng)數(shù)值分析與綜合［M］.北京：機(jī)械工業(yè)出版社，2009：81－90.

［14］何為，李倩，徐征，等.頭部分層球模型磁感應(yīng)斷層成像正問(wèn)題的解析解［J］.計(jì)算物理，2010，27（6）：912－918.

［15］秦明新.檢測(cè)腦水腫的磁感應(yīng)成像測(cè)量方法研究［D］.西安：西安電子科技大學(xué)，2005.