ABSTRACT

Magnetic induction tomography (MIT) is an imaging technique based on the measurement of the magnetic field perturbation due to eddy currents induced in conducting objects exposed to an external magnetic excitation field. When the driving frequency is significantly high relative to the frequency range that MIT normally operates, a highly conducting and permeable metallic target between the coils can be treated as perfect electric conductors (PEC). In this scenario, the penetration depth of the magnetic field into the target is extremely small and Boundary Element Method (BEM) based on integral formulations becomes an effective way to analyze this kind of scattering problems. However, BEM is most efficient when the scatter object is small rather than distributed. This is not suitable for the case when a surrounding shield of significant size and surface area is incorporated in the BEM model, which negates any benefits due to the BEM formulation. For this reason, we proposed a hybrid method in this paper, which combines the BEM and the Finite Element Method (FEM). FEM is used to calculate the shield effect while BEM is used to calculate the perturbation due to a small PEC inside the sensing space. In this way, we are able to compute sensitivity maps for MIT due to the effect of different shield configurations. An appropriate index related to the sensitivity maps was proposed to obtain the optimal parameters of the shield.

Keywords:Magnetic Induction Tomography, Perfect Electric Conductor, BEM-FEM, Sensitivity Map, Shield

利用BEM-FEM混合法对高频电磁层析成像
系统中金属屏蔽层的优化分析

尹武良^1,2，赵 倩¹

¹天津大学电气与自动化工程学院，天津

²曼彻斯特大学电气与电子工程学院，曼彻斯特，英国

Email:shmilyshenzhen@163.com

收稿日期：2014年3月1日；修回日期：2014年4月2日；录用日期：2014年4月12日

摘 要

电磁层析成像技术是一种基于电磁感应原理的成像技术。外加激励电源在导电物体中产生涡流，进而生成二次磁场，该磁场反过来会影响到原有磁场的分布。当激励频率非常高时，导电物体可以看作是完美导体。此时，磁场对物体的渗透深度非常浅，从而基于边界积分方程的边界元法便成为一种有效的分析方法。但是如果电磁层析成像系统中存在一个无法忽略尺寸的屏蔽层时，边界元法就失去了它的优势。针对这种情况，我们提出了一种边界元–有限元混合法。其中，有限元法用于计算屏蔽层对主磁场的影响，边界元法用于计算完美导体在被测区域内移动时接收线圈上的感应电压。通过这种方法，可以计算出在不同参数的屏蔽层的影响下，电磁层析成像系统的灵敏度分布，进而通过合适的灵敏度评估参数对屏蔽层进行优化设置。

关键词

电磁层析成像，完美导体，边界元–有限元法，灵敏度，屏蔽层

1. 引言

电磁层析成像技术(Magnetic induction tomography，简称MIT)是一种非侵入非接触的成像技术，在工业和医学上都有很广泛的应用前景[1] -[8] 。MIT通过构建被测物体的电学参数来得到其分布。MIT的正问题对于重建图像有非常重要的作用，因此引起了相关研究学者的关注[9] 。对于一般的工业MIT系统，激励频率通常在1 kHz~100 kHz之间。然而在医学中，由于被测溶液的电导率往往比较低，所以为了得到更为精确的图像，设计出了高频MIT系统，即激励频率 > MHz [1] [8] 。

对于金属物体而言，当激励频率很高时，它可以看作是完美导体(Perfect electric conductor，简称PEC)。此时，入射磁场基本上全部被反射，磁场的渗透深度非常浅。有限元法(Finite elements method，简称FEM)就不再适合求解这类物体的电磁问题，因为需要大量的网格以取得高精度的结果，造成很重的计算负担。相较之下，基于边界积分方程的边界元法(Boundary element method，简称BEM)只需要在物体表面进行网格划分，所以求解起来更快速简单[10] -[12] 。迄今为止，已有不少学者对利用积分方程求解MIT的正问题进行了研究[13] -[15] 。

因为高频MIT系统中激励线圈会与被测物体间产生电学耦合现象，给求解过程造成一定的困难。对此，一个最常用的解决方法就是在传感器阵列的外面放置一个导电却不导磁的金属屏蔽层，通常选用铝或铜[16] 。由于屏蔽层的加入导致BEM的计算优势丧失，为了解决此问题，我们提出了一种边界元-有限元混合方法。

这种方法的原创之处在于，作者并非自行推导有限元法理论，而是使用了一款有限元商业软件MAXWELL ANSOFT作为有限元求解的工具，与自行推导的边界元软件进行了接口，使得既利用了商业软件处理一些复杂形状的能力，又综合了边界元的一些特性。是一种快捷简单的混合方法。这种方法有别于传统意义上的边界元–有限元混合法，但仍然可以称作边界元–有限元混合法。

首先，我们通过MAXWELL ANSOFT求得主磁场的分布及其相关的标量磁势。进而在被测区域中利用物体表面离散点处的标量磁势构建积分方程，根据一定的边界条件，可以求出物体表面的标量磁势，在此基础上，被测物体外面的磁场分布以及接收线圈上的感应电压就可以求出。当PEC在被测区域内移动时，可以计算出不同参数的屏蔽层对系统灵敏度分布的影响。通过对评估参数的分析，可以对屏蔽层进行优化处理。

2. PEC问题的MIT模型

2.1. 边界积分方程

总磁场包括两部分：由激励线圈产生的主磁场以及由涡流产生的二次磁场

(1)

麦克斯韦–安培方程

(2)

这里H和E分别代表磁场强度和电场强度。上标e、pr和sc分别表示与外部空间、主磁场和二次磁场相关的量。

针对一般的MIT系统模型，我们有以下假设：

1) 外部磁场是似稳的，故而可以忽略位移电流。

2) 物体外部的空间具有绝缘性，并且被测区域不包含附加电源。

基于以上假设，由公式(2)可得，外部磁场的旋度为0，可用标量磁势来表示[12] 。同时由散度定理

(3)

可知满足Laplace方程，解式(3)可得边界积分方程为：

(4)

其中是Laplace方程的格林函数，是主磁场的标量磁势。矢径r和r'分别表示场点

和源点的位置。表示物体的表面，其单位法向量记为n。参数c的值决定于场点的位置。当场点位于光滑表面上时，c = 1/2；位于内时，c = 1；在其他位置时，c = 0[13] 。

由于PEC内部不存在磁场，故根据物体表面处的边界条件可将式(4)简化为关于的积分方程：

(5)

在数值求解过程中，物体表面被划分为NE个三角形平面网格，第e个网格记为。这时，式(5)便可以转化为离散点处的积分方程组：

(6)

其中NP为离散点的数目。对进行线性插值得：

(7)

其中N为插值函数，为上第k个点处的标量磁势。

为了将局部编号k转化为整体编号j，我们引入了一个转换矩阵：

(8)

其中是整体坐标系中第j个点处的标量磁势。故可得最终的积分方程为：

(9)

求解式(9)涉及到若干个复杂的积分，当网格数目很多时，会给计算机带来很大的负担。为了简化计算，我们在三角形网格中建立局部坐标系，将原有积分转换为新坐标系中的积分[17] [18] 。

2.2. ANSOFT仿真

ANSOFT是一个工程中用途较为广泛的电磁场仿真软件，通过使用精确的有限元法来求解电磁场分布。

2.2.1. 系统模型

图1描述了MIT的系统结构及金属屏蔽层的形状参数。激励线圈的半径为0.1 m，中心位于点(0, 0, 0.6 m)处。系统中有4个接收线圈，半径均为0.1 m，位置如图1(c)所示。激励线圈和接收线圈组成传感器阵列。MIT系统中激励电流的频率为10 MHz，幅值为1 A。保持金属屏蔽层和传感器的位置不变，在被测区域(,)内移动PEC，可以得到物体表面离散点处的磁场分布。

屏蔽层的主要参数有：

t：厚度，变化范围是0.05:0.05:0.25，其中中间的0.05表示变化间距，下同。

h：高度，变化范围是0.4:0.2:1.4。

d：内半径，变化范围是0.7:0.1:1.2。

对于每一种屏蔽层的参数配置，我们均会通过仿真计算出其对MIT系统灵敏度分布的影响。首先得到PEC表面的离散点集，当PEC在被测区域内移动时，点集也会发生改变。将包含三个坐标值的点集输入到ANSOFT软件中，利用场计算器可以求解出主磁场的分布。

2.2.2. 主磁场的标量磁势

一旦求出主磁场，则我们可以根据下面的步骤计算出标量磁势。

图2中的图形用来求解与主磁场相关的标量磁势，当激励线圈平行于XY平面时，选取球体表面的一条曲线为例。在球面上共有38个离散点，除了上下顶点，还有9 ´ 4个点。其中9代表不同的z值，4代表每个z值对应的四个点。图2中的曲线上共有11个点，从上到下依次排列。令和分别代表第i点处的主磁场强度以及标量磁势。

图1. MIT系统结构及屏蔽层的形状。(a) MIT系统的俯视图；(b) 屏蔽层的侧视图；(c) ANSOFT模型

图2. 求解标量磁势的物体模型

在这里，我们假设点6处的标量磁势为0，即。由于，故

(10)

其中l是从点a到点b的积分路径，假设点a到点b间的磁场是恒定的，并定义为a，b两点处磁场强度的平均值，则可以得到

(11)

由式(11)，可以求出所有离散点处的。进而通过式(9)求得。

2.3. 计算灵敏度分布

分别对式(9)的两端求梯度，可以得到物体外部的磁场分布：

(12)

其中为主磁场，为物体表面第k点处的标量磁势。接收线圈上的感应电压为：

(13)

其中U为接收线圈上的感应电压。为接收线圈围成的平面，为其单位法向量。将划分为N个小网格，并将第i个小网格记为。

灵敏度分布可以通过下式来求[15] ：

(14)

其中S_en为被测区域的灵敏度，，和分别为总电压、由二次磁场产生的电压以及由主磁场产生的电压。在被测区域内移动PEC可以得到不同的接收电压，进而通过式(14)求得MIT系统的灵敏度分布。

3. 仿真及结果分析

3.1. 仿真方法及优化参数

在这里，我们采用正交阵列法[19] 来推导MIT系统中金属屏蔽层的最优设计。三个评估参数分别为：

平均灵敏度：

灵敏度的标准差：

灵敏度的相对标准差：

其中为第e个网格处的灵敏度值。下标i和j分别为激励线圈和接收线圈的编号。一般情况下，K越小则表示灵敏度分布越均匀。

如果将四类典型灵敏度的相对标准差分别记为、、和，它们的和记为K：

(15)

其中下标1表示激励线圈，2~5表示四个不同位置的接收线圈。由于当、、和达到最小值时，t往往不同，故为了得到一个更全面的结果，我们使用参数K进行分析。

3.2. 仿真步骤

仿真通过以下四步完成：

1) 厚度t的影响

首先将参数h和d固定，在这里，选取h和d为其变化范围的中间值，即h = 0.9，d = 0.9。t的变化范围是0.05:0.05:0.25，通过对比不同t时的K可以得到最优的t，并记为t = t_opt。同时选取K值较小时对应的三个t值，分别记为t₁、t₂和t₃。仿真结果如图3所示。

通过图3可以看出，当、、和达到最小值时对应不同的t值，综合考虑下，它们之和K达到最小值10.9614时，t = 0.2，故可知t_opt= 0.2。

2) 高度h的影响

首先将参数t和d固定，t = t_opt，d = 0.9。h的变化范围是0.4:0.2:1.4。同(1)中的仿真原理类似，可得h_opt= 0.8，同时，选取三个较优的h记为h₁、h₂和h₃。仿真结果如图4所示。

3) 内半径d的影响

首先将参数t和h固定，t = t_opt，h = h_opt。d的变化范围是0.7:0.1:1.2。由图5可以得到最优的d_opt= 0.8，

图3. 不同t值对应的K

图4. 不同h值对应的K

同时选取较优的三个d记为d₁、d₂和d₃。

4) 对三个较优t、h和d的正交测试

由图3、图4和图5可以看出，参数t、h和d对灵敏度分布的作用是非线性的。为了得到最优的组合，我们采用正交测试法，具体的仿真步骤如表1所示。

图5. 不同d值对应的K

通过步骤1~3，我们可以得到较优的三个t、h和d分别记为：

正交测试的结果如表2所示。

参数Rp表示参数对结果的影响程度，由表2可以看出，屏蔽层高度h对K的影响最大，其次是厚度t，内半径d的影响最小。

由图6可以看出，最优的参数组合是(t₃, h₂, d₃)：t₃= 0.2；h₂= 0.6；d₃= 0.9。在9种不同参数的屏蔽层

图6. 屏蔽层参数与优化指标间的关系

存在的情况下，MIT的灵敏度图如图7所示。

图7. 正交测试中得到的四种典型灵敏度图

4. 结论

本文中提出了一种BEM-FEM混合法，该方法可以对高频MIT系统中金属屏蔽层参数进行有效优化。为了设计出最优的屏蔽层，需要考虑很多因素，例如屏蔽层的形状和位置等。在FEM软件MAXWELLANSOFT的帮助下，我们首先求得了在不同参数屏蔽层存在时的主磁场分布，接着利用BEM和扰动法计算出MIT系统的灵敏度分布。通过对灵敏度图的分析，可以得到不同屏蔽层参数对灵敏度分布的影响。选取合适的优化指标，就可以得到最优的屏蔽层参数，同时由仿真结果也可以看出，屏蔽层的设计对系统的灵敏度分布有较大影响。

参考文献 (References)