趙 倩，郝建娜，尹武良

(天津大學(xué)電氣與自動(dòng)化工程學(xué)院，天津 300072)

電磁層析成像系統(tǒng)中標(biāo)量磁勢(shì)的數(shù)值解法

趙 倩，郝建娜，尹武良

(天津大學(xué)電氣與自動(dòng)化工程學(xué)院，天津 300072)

電磁層析成像技術(shù)是一種基于電磁感應(yīng)定律的工業(yè)過程成像技術(shù)，激勵(lì)線圈產(chǎn)生的交變磁場(chǎng)在目標(biāo)物體中產(chǎn)生渦流，進(jìn)而產(chǎn)生二次磁場(chǎng)．接收線圈檢測(cè)到感應(yīng)電壓后利用重建算法可以得到物場(chǎng)的分布信息．邊界元法以積分方程為數(shù)學(xué)基礎(chǔ)，同時(shí)采用了與有限元法相似的劃分單元離散技術(shù)，將邊界積分方程離散為代數(shù)方程組后用數(shù)值方法求解．邊界元法在電磁層析成像技術(shù)中已有一定的應(yīng)用，而基于標(biāo)量磁勢(shì)的邊界元法使得求解過程的速度和效率顯著提高．本文針對(duì)一個(gè)簡(jiǎn)化的電磁層析成像系統(tǒng)模型，利用3種數(shù)值方法求解目標(biāo)物體表面的標(biāo)量磁勢(shì)，并利用Matlab編程得到仿真結(jié)果．對(duì)結(jié)果進(jìn)行對(duì)比分析后，選出最優(yōu)解法．

電磁層析成像；邊界元法；標(biāo)量磁勢(shì)；數(shù)值仿真

電磁層析成像(electromagnetic tomography，EMT)技術(shù)是20世紀(jì)90年代發(fā)展起來(lái)的一種新型電學(xué)層析成像(electrical tomography，ET)技術(shù)，EMT技術(shù)基于電磁感應(yīng)原理，載流線圈在目標(biāo)物體中產(chǎn)生渦流，進(jìn)而激勵(lì)出二次磁場(chǎng)，這樣接收線圈可以測(cè)量到感應(yīng)電壓．最后采用一定的圖像重建算法得到被測(cè)場(chǎng)域內(nèi)的分布信息．它具有非接觸、非侵入、響應(yīng)快速、信息量豐富等特點(diǎn)，因而在工業(yè)生產(chǎn)過程和生物醫(yī)學(xué)領(lǐng)域都有著廣泛的應(yīng)用前景[1-7]．

邊界元法(boundary element method，BEM)具有邊界積分方程的深厚的數(shù)學(xué)根基．它以邊界積分方程為數(shù)學(xué)基礎(chǔ)，通過單元離散技術(shù)將邊界離散為邊界元，將邊界積分方程離散為代數(shù)方程組．再用數(shù)值方法求解代數(shù)方程組，從而得到原問題邊界積分方程的解．邊界元法的優(yōu)點(diǎn)為可以降低求解問題的維數(shù)，提高了求解效率，并且由于求解過程中采用了解析基本解，所以具有更高的精度[8-14]．

BEM在EMT中的應(yīng)用已經(jīng)有了許多的研究[15].在邊界元法求解電磁問題的過程中，標(biāo)量磁勢(shì)作為基本變量具有重要意義．在求解外部磁場(chǎng)分布之前，需要先求出目標(biāo)物體表面的標(biāo)量磁勢(shì)，以此為基礎(chǔ)再進(jìn)行后續(xù)計(jì)算．目前對(duì)于標(biāo)量磁勢(shì)的求解方法一般表示為對(duì)要求點(diǎn)處求立體角[16]，該方法同樣涉及如何求解積分的問題．對(duì)于復(fù)雜不規(guī)則物體，快速多極算法是一種快速有效的計(jì)算方法[17]．但是該算法的步驟較為繁瑣，對(duì)于一些簡(jiǎn)單問題反倒不適用．

筆者通過建立一個(gè)簡(jiǎn)單的EMT系統(tǒng)模型，推導(dǎo)出標(biāo)量磁勢(shì)的最終表達(dá)式．通過3種不同的數(shù)值方法——曲面積分定義法、高斯-勒讓德積分公式法和積分坐標(biāo)變換法，分別仿真得到不同的結(jié)果．對(duì)這些數(shù)據(jù)進(jìn)行分析可以得出，在保證其他條件不變和相同求解精度的前提下，高斯-勒讓德積分公式法是最簡(jiǎn)單快速的方法．

1 建立模型并推導(dǎo)公式

假設(shè)一個(gè)簡(jiǎn)單的EMT系統(tǒng)，只含有1個(gè)激勵(lì)線圈和1個(gè)接收線圈．激勵(lì)線圈和接收線圈的中心分別位于(0，0，4,m)和(0，0，-4,m)處，半徑均為1,m．目標(biāo)球體為一般金屬，中心位于原點(diǎn)，半徑為1,m．

對(duì)于恒定磁場(chǎng)，由磁場(chǎng)環(huán)路定理可得

式中L為曲面S的邊界．這種情況下如果想引入標(biāo)量磁勢(shì)，其前提條件是，對(duì)于求解區(qū)域內(nèi)的任何閉合回路，都有

根據(jù)上面推導(dǎo)出的公式，這里采用了3種方法來(lái)求解球體表面離散點(diǎn)處的標(biāo)量磁勢(shì)．之后利用Matlab仿真，對(duì)這3種方法的求解速度和所得數(shù)據(jù)精度進(jìn)行了分析對(duì)比．

圖1 一個(gè)簡(jiǎn)單的EMT系統(tǒng)Fig.1 A simple EMT system

在實(shí)際應(yīng)用中，對(duì)于單連通，不存在自由電流的區(qū)域，全區(qū)域皆可以引入標(biāo)量磁勢(shì)描述磁場(chǎng)．如果區(qū)域中存在自由電流，則挖去電流及電流線圍著的一個(gè)曲面S，剩余空間里磁場(chǎng)是無(wú)旋的[18-19]，可以假設(shè)為標(biāo)量磁勢(shì)，A．

圖1所示的EMT系統(tǒng)中，載流線圈在物體表面產(chǎn)生的磁場(chǎng)可由畢奧-薩伐兒定理求出[19]，即

式中：上標(biāo)pr表示外部激勵(lì)產(chǎn)生的場(chǎng)量；B為磁感應(yīng)強(qiáng)度；r=r1-r2，r1為球面上的點(diǎn)到原點(diǎn)的矢徑，r2為激勵(lì)線圈所圍平面上的點(diǎn)到原點(diǎn)的矢徑．

式中：S為由激勵(lì)線圈L圍成的平面；n0為S的單位法向量．

對(duì)場(chǎng)域中的球體表面進(jìn)行三角形網(wǎng)格劃分，如圖2所示，可以得到38個(gè)離散點(diǎn)和72個(gè)網(wǎng)格．劃分時(shí)，縱軸z的范圍為1,m到-1,m，每隔0.2,m取一組，對(duì)應(yīng)11個(gè)不同的z值．在每個(gè)z確定的圓上，角度范圍為0到2π，間隔為π/2．

圖2 EMT系統(tǒng)中球體表面的網(wǎng)格劃分Fig.2 Mesh of the target surface in the EMT system

若激勵(lì)線圈上載有單位電流，即I=| I |=1,A，則球體表面由外加磁場(chǎng)產(chǎn)生的標(biāo)量磁勢(shì)為

2 Matlab數(shù)值仿真

2.1 由曲面積分的定義求解(方法1)

由曲面積分的定義可得[20]，通過將積分區(qū)域劃分成正方形網(wǎng)格，可將積分離散化，有

式中：iSΔ為第i個(gè)正方形網(wǎng)格的面積；2ir為激勵(lì)線圈所圍平面上的點(diǎn)到原點(diǎn)的矢徑．在式(6)中，積分平面為激勵(lì)線圈所包圍的圓面，將該平面分成N個(gè)小正方形，r2i為矢徑取第i個(gè)正方形的中心到原點(diǎn)的矢徑．首先對(duì)r2i的位置進(jìn)行判斷，若是在圓內(nèi)則將其帶入積分公式；若在圓外，則其對(duì)積分的作用為0；再將所有的積分進(jìn)行求和，最終得到在球體表面每個(gè)離散點(diǎn)處的標(biāo)量磁勢(shì)．

由仿真結(jié)果可得，隨著劃分網(wǎng)格數(shù)目的增多，積分結(jié)果成收斂趨勢(shì)．當(dāng)網(wǎng)格數(shù)為6,400時(shí)，積分基本不再變，計(jì)算時(shí)間為t=0.556,122,s．38個(gè)離散點(diǎn)對(duì)應(yīng)11個(gè)z值，通過仿真可知，每個(gè)z值相同的點(diǎn)具有相等的標(biāo)量磁勢(shì)，故可用11個(gè)不同z值上的標(biāo)量磁勢(shì)來(lái)表示所有離散點(diǎn)處的標(biāo)量磁勢(shì)．

表1 方法1得到的標(biāo)量磁勢(shì)Tab.1 Magnetic scalar potential obtained by the first method

以表1第1行為例，可以畫出隨著網(wǎng)格數(shù)增加標(biāo)量磁勢(shì)逐漸收斂的仿真圖形，如圖3所示．

由于每個(gè)確定的z值對(duì)應(yīng)的圓上的4個(gè)點(diǎn)處具有相同的標(biāo)量磁勢(shì)，所以下面的表格中可以只列出不同z值處的標(biāo)量磁勢(shì)．圖4為x=0處的球體切面上的等勢(shì)線，可以看出，球體從上往下標(biāo)量磁勢(shì)逐漸減小，并且在每個(gè)確定z值處具有相同的標(biāo)量磁勢(shì)．

圖3 隨著網(wǎng)格數(shù)標(biāo)量磁勢(shì)數(shù)值逐漸收斂Fig.3Convergence of the magnetic scalar potential as the increase of mesh

圖4 x=0處的切面上等勢(shì)線分布Fig.4 Isopotential lines on the section where x=0

2.2 利用N段高斯-勒讓德積分公式來(lái)求解(方法2)

首先將面積分轉(zhuǎn)換成二重積分[21]為

則原積分可以簡(jiǎn)寫為

積分變量θ和r′分別表示角度和半徑，積分函數(shù)為

對(duì)積分區(qū)間進(jìn)行標(biāo)準(zhǔn)化得

變換方程為

則

式中：n為節(jié)點(diǎn)數(shù)；k為插值節(jié)點(diǎn)數(shù)；Ai和Aj(i，j=1，2，…，R)均為系數(shù)矩陣A中第i個(gè)和第j個(gè)值．具體的節(jié)點(diǎn)和對(duì)應(yīng)的系數(shù)矩陣如表2所示．為了使精度更高，在這里先將積分區(qū)域劃分為N段小區(qū)域，然后在每段小區(qū)域上進(jìn)行插值積分運(yùn)算．運(yùn)行程序得到的結(jié)果表明，當(dāng)N=2，即分成2段，然后進(jìn)行5點(diǎn)插值可得高精度的解，用時(shí)為t=0.109,799,s．

表2 高斯-勒讓德積分公式的節(jié)點(diǎn)及系數(shù)矩陣Tab.2 Nodes and coefficient matrices of Gaussian-Legendre integral formula

由表3可知，對(duì)于N(N=1，2，4)段積分公式而言，所求得的標(biāo)量磁勢(shì)基本與N無(wú)關(guān)．只有在z=0時(shí)，不同N對(duì)應(yīng)的求解結(jié)果才有細(xì)微差別．故可知，利用N段高斯-勒讓德積分公式來(lái)求解標(biāo)量磁勢(shì)時(shí)，段數(shù)N對(duì)結(jié)果精度的影響甚微，從而為了節(jié)約計(jì)算時(shí)間，可以直接使用1段5點(diǎn)積分公式來(lái)求解．

表3 方法2得到的標(biāo)量磁勢(shì)Tab.3 Magnetic scalar potential obtained by the second method

2.3 利用積分坐標(biāo)變換法來(lái)求解(方法3)

這種方法首先將積分平面劃分為三角形網(wǎng)格，然后利用坐標(biāo)系轉(zhuǎn)換得到的積分公式來(lái)求解[22-24]．

劃分網(wǎng)格時(shí)，先將圓面上的點(diǎn)進(jìn)行離散，然后按照逆時(shí)針順序?qū)Ⅻc(diǎn)依次連接為三角形．根據(jù)其中心點(diǎn)的位置，可以判斷出三角形是否在圓內(nèi)．為了方便求解，在這里將全局坐標(biāo)系轉(zhuǎn)換為局部坐標(biāo)系．對(duì)每個(gè)網(wǎng)格而言，新坐標(biāo)系平行于該網(wǎng)格所在的平面，并且以其中的一邊為其橫坐標(biāo)．坐標(biāo)系的各參數(shù)項(xiàng)如圖5所示．

圖5 局部坐標(biāo)系與全局坐標(biāo)系的關(guān)系及局部坐標(biāo)系中的參數(shù)Fig.5 Relationship between the local coordinate system and the global coordinate system and the parameters of the local coordinate system

采用與文獻(xiàn)[22]中類似的參數(shù)和函數(shù)，這里以一個(gè)三角形網(wǎng)格為例．

V1、V2和V3為第e個(gè)三角形網(wǎng)格eΓ的3個(gè)頂點(diǎn)．iS?為eΓ的第i條邊．li、si和mi(i=1，2，3)分別是iS?的長(zhǎng)度、單位切向量和單位法向量．則局部坐標(biāo)系可以定義為

根據(jù)以上所列參數(shù)，則原積分可化為

其中n0=(0,0,1)，ξ′=0．故有

式中Ne為三角形網(wǎng)格的個(gè)數(shù)．

由文獻(xiàn)[22]可得eΓ上的積分為

數(shù)值計(jì)算表明，當(dāng)網(wǎng)格數(shù)為12,800時(shí)，標(biāo)量磁勢(shì)可以收斂，但用時(shí)較長(zhǎng)，t=157.126,357,s．

表4 方法3得到的標(biāo)量磁勢(shì)Tab.4 Magnetic scalar potential obtained by the third method

2.4 利用Matlab進(jìn)行數(shù)據(jù)對(duì)比分析

利用Matlab軟件對(duì)3種方法得到的數(shù)據(jù)進(jìn)行分析，如圖6所示．

從圖6可以看出，3種方法得到的數(shù)值解具有很好的一致性，即3種方法都可以得到精確的解，從程序運(yùn)行時(shí)間上看，高斯-勒讓德積分公式用時(shí)最短，t=0.109,799,s，積分坐標(biāo)變換法用時(shí)最長(zhǎng)，t= 157.126,357,s．

圖6 3種方法得到的數(shù)據(jù)曲線Fig.6 Results obtained by the three methods

3 結(jié) 語(yǔ)

快速求解標(biāo)量磁勢(shì)為后續(xù)計(jì)算磁場(chǎng)分布提供了良好基礎(chǔ)．本文采用了3種數(shù)值方法來(lái)求解物體表面的標(biāo)量磁勢(shì)，通過以上計(jì)算可以看出，以上3種方法都能有效求解目標(biāo)物體表面的標(biāo)量磁勢(shì)．若綜合考慮的話，高斯-勒讓德積分公式是最好的計(jì)算方法，不僅精度高，而且速度快．

［1］ Griffiths H. Magnetic induction tomography[J]. Meas Sci Technol，2001，12：1126-1131.

［2］ Watson S，Wee H C，Griffiths H，et al. A highly phase-stable differential detector amplifier for magnetic induction tomography[J]. Physiol Meas，2011，32：917-926.

［3］ Watson S，Williams R J，Gough W A，et al. A magnetic induction tomography system for samples with conductivities less than 10 S·m-1[J]. Meas Sci Technol，2008，19：1-11.

［4］ Yin W L，Peyton A J. A planar EMT system for the detection of faults on thin metallic plates[J]. Meas Sci Technol，2006，17：2130-2135.

［5］ Soleimani M，Lionheart W R B，Riedel C H，et al. Forward problem in 3,D magnetic induction tomography(MIT)[C]//Proc 3rd World Congr Industrial Process Tomography. Banff，Canada，2003：275-280.

［6］ Yin W L，Peyton A J. Sensitivity formulation including velocity effects for electromagnetic induction systems[J]. IEEE Transactions on Magnetics，2010，46(5)：1172-1176.

［7］ Yin W L，Peyton A J，Zysko G，et al. Simultaneous noncontact measurement of water-level and conductivity [J]. IEEE Transactions on Instrumentation and Measurement，2008，57：2665-2669.

［8］ Brebbia C A，Telles J C F，Wrobel L C. Boundary Element Techniques：Theory and Applications in Engineering[M]. New York：Springer -Verlag，1984.

［9］ 姚振漢，王海濤. 邊界元法[M]. 北京：高等教育出版社，2010.

Yao Zhenhan，Wang Haitao. The Boundary Element Method[M]. Beijing：Higher Education Press，2010(in Chinese).

［10］ 劉永健，姚振漢. 三維接觸邊界元法的一種誤差直接估計(jì)[J]. 清華大學(xué)學(xué)報(bào)：自然科學(xué)版，2003，43(11)：1499-1502.

Liu Yongjian，Yao Zhenhan. Direct error estimate for the boundary element method for a 3-D contact problem[J]. Journal of Tsinghua Univ：Sci & Tech，2003，43(11)：1499-1502(in Chinese).

［11］ 單 磊，董天臨. 應(yīng)用邊界元法對(duì)電磁場(chǎng)計(jì)算的分析及優(yōu)化[J]. 信息通信，2007，33(1)：31-33.

Shan Lei，Dong Tianlin. Analysis and optimization of boundary element method to irregular electromagnetic field[J]. Information & Communications，2007，33(1)：31-33(in Chinese).

［12］ Harrington R F. Field Computation by Moment Methods[M]. New York：Acmillan Company，1968.

［13］ Banerjee P K. The Boundary Element Methods in Engineering[M]. New York：McGraw-Hill Book Company，1994.

［14］ Chen W，F(xiàn)u Z. A novel numerical method for infinite domain potential problems[J]. Chinese Science Bulletin，2010，55(16)：1598-1603.

［15］ Pham M H，Peyton A J. A model for the forward problem in magnetic induction tomography using boundary integral equations[J]. IEEE Transactions on Magnetics，2008，44(10)：2262 -2267.

［16］ 張澤明，鄧亦仁. 運(yùn)用畢-沙定律導(dǎo)出磁標(biāo)量位表達(dá)式的新方法[J]. 重慶建筑工程學(xué)院學(xué)報(bào)，1985(4)：99-103.

Zhang Zeming，Deng Yiren. A new method using Biot-Savart Law to derive magnetic scalar potential notation[J]. Journal of Chongqing Institute of Civil Engineering and Architecture，1985(4)：99-103(in Chinese).

［17］ Hafla W，Buchau A，Rucker W M. Efficient computation of source magnetic scalar potential [J]. Adv Radio Sci，2006，4：59-63.

［18］ 阮放鳴. 有限區(qū)域里的磁標(biāo)勢(shì)和磁多極展開[J]. 貴州師范大學(xué)學(xué)報(bào)：自然科學(xué)版，1997，15(3)：102-107.

Ruan Fangming. Magnetic scalar potentials in finite regious and the multipole expansion in material media[J]. Journal of Guizhou Normal University：Natural Science，1997，15(3)：102-107(in Chinese).

［19］ 謝處方，饒克謹(jǐn). 電磁場(chǎng)與電磁波[M]. 北京：高等教育出版社，1987.

Xie Chufang，Rao Kejin. The Electromagnetic Field and Waves[M]. Beijing：Higher Education Press，1987(in Chinese).

［20］ 同濟(jì)大學(xué)數(shù)學(xué)教研室. 高等數(shù)學(xué) [M]. 4版. 北京：高等教育出版社，1996.

The Mathematical Teaching Team of Tongji University. Advanced Mathematics[M]. 4th ed. Beijing：Higher Education Press，1996(in Chinese).

［21］ 張 軍. 數(shù)值計(jì)算[M]. 北京：清華大學(xué)出版社，2008.

Zhang Jun. Numerical Computation[M]. Beijing：Tsinghua University Press，2008(in Chinese).

［22］ Graglia R D. On the numerical integration of the linear shape functions times the 3-D Green's function or its gradient on a plane triangle [J]. IEEE Transactions on Antennas and Propagation，1993，41(10)：1448-1455.

［23］ Morrison J A. Integral equations for electromagnetic scattering by perfect conductors with two-dimensional geometry [J]. Bell Syst Tech J，1979，58(2)：409-425.

［24］ Graglia R D. Static and dynamic potential integrals for linearly varying source distributions in two- and threedimensional problems [J]. IEEE Transactions on Antennas and Propagation，1987，AP-35(6)：662-669.

(責(zé)任編輯：孫立華)

Numerical Approach for the Magnetic Scalar Potential of Electromagnetic Tomography System

Zhao Qian，Hao Jianna，Yin Wuliang
(School of Electrical Engineering and Automation，Tianjin University，Tianjin 300072，China)

Electromagnetic tomography (EMT) is an industrial process tomography based on the Faraday law of electromagnetic induction. In EMT, current-carrying coils are used to induce eddy currents in the object which in turn radiates a scattered field. The induced voltages are sensed with receiving coils and image reconstruction algorithms are used to solve for the spatial distribution of the complex conductivity in the target. Boundary element method (BEM) based on integral formulations is an effective way to analyze the problems of EMT system. By adopting the similar mesh technique as FEM, the boundary integral equations can be transformed into linear equations. Then numerical method is used to solve the linear equations and the solution of the original integral equations can be obtained. The magnetic scalar potential was the important part of the application of BEM into EMT which improves efficiency and speed. Three numerical approaches were used to calculate the magnetic scalar potential on the surface of the target for a simple EMT model. With the aid of Matlab, the optimal method was chosen by comparing the results obtained by the three methods.

electromagnetic tomography；boundary element method；magnetic scalar potential；numerical simulation

O241.8

A

0493-2137(2014)07-0613-06

10.11784/tdxbz201210003

2012-10-06；

2013-01-23.

國(guó)家自然科學(xué)基金國(guó)際重大合作資助項(xiàng)目(60910001)；教育部博士點(diǎn)基金資助項(xiàng)目(20090032110062).

趙 倩（1987— ），女，博士研究生，shmilyshenzhen@163.com.

尹武良，wuliang.yin@gmail.com.