1. 引言
电磁学自产生起就不断给社会带来巨大的变革,影响到生活的方方面面,小到指南针,发电机、电动机,大到通信系统、航天技术等。麦克斯韦方程组描述电磁场的基本规律,是电磁学的精华,也是电磁学广泛应用的理论基础[1] 。从数学和物理的角度看,电磁场相关的应用研究中涉及正反两类问题。正问题指已知介质分布,求电场分布,该问题已形成完整的理论体系,有各种解析方法和数值方法可供选择。逆问题指已知场分布,求介质分布,在大多数情况下都是非适定的,求解非常困难且没有普适性的方法。电阻抗成像(Electrical Impedance Tomography,简称EIT)是电磁学的应用之一,它是通过给导电物体注入电流后测量表面的电势分布来确定导电物体内部的阻抗分布,本质是电磁学的逆问题,在地球物理,医学成像以及无损检测等领域有重要应用[2] -[4] 。
电阻抗成像方法包括动态成像[5] 、静态成像[6] [7] 以及准静态成像[8] [9] (多频成像);动态成像是利用“两个不同时刻测量数据的差值”来重构电导率的差分图像,该成像方法能消除噪声误差,对系统精度要求较低,计算量小,但在两个时刻电导率不变时没法成像,应用范围较窄。静态成像直接对电导率绝对值成像,应用范围广,但计算量大,对噪声敏感,图像的对比度和分辨率差。准静态成像是在同一时刻向被测生物组织施加多种频率的激励,利用不同频率下的测量数据进行差分成像。
本文围绕静态成像问题开展研究。静态电阻抗成像的过程是先假设“阻抗分布”计算出边界电势分布(求解非线性方程组),用计算值和测量得到的边界值比较,反复进行“假设-计算-比较”直到误差小到一个预设的范围为止,然后对最后假设的阻抗分布成像。在求解过程中,对初始阻抗分布的选取很重要,选的初始值合适,计算量会大幅减小,且能提高成像质量,初值不合适,不光计算量大,有时还可能算不出结果。如果能够先锁定目标区域的大致范围,则选好初值的概率将大幅提升。本文利用对称性分析确定单目标成像区域的大致范围,可为EIT初值选取提供良好的参考,为定性的无损检测提供依据。
此外,现有的静态成像对中心区域的成像效果很差,结合对称性分析则能很好地克服这个问题。
对称性分析是物理学的一种重要研究手段,已形成完整的理论体系——群论[10] 。对称性通常都对应着一些重要规律。本文用Comsol multiphysics 4.3a设计实验系统,并计算出均匀区域内分布单目标时边界电极的电势,通过分析对称点之间的电势差和目标位置、大小的关系,找到边界电势的对称性与目标区域定位信息的关联,最终能通过边界电势对称性的分析确定目标区域的位置、大小,或者减小确定目标区域的位置和大小的工作量。
2. 系统描述
2.1. 实验模型
本文采用的模拟实验系统包括一个电导率固定的大薄圆盘(大圆),在大圆上挖一个小圆并补上不同电导率的介质,大圆边界上每隔π/8一个电极,分别标上1,2,...,16,1-9,2-10,3-11,4-12,5-13,6-14,7-15,8-16分别对应每隔π/8弧度直径的端点。大圆半径R = 10 cm,厚度0.1 mm,电导率σ = 1 s/m,输入电流为10 mA,目标小圆电导率σ0 =0.001 s/m,半径为a,圆心的极坐标为(r, θ)。具体如图1所示。
2.2. 对称性
本文的对称性指关于x轴对称的电极点之间的电势差,对16电极系统共有7组,分别是V2-V16,V3-V15,V4-V14,V5-V13,V6-V12,V7-V11,V8-V10,它们的值为零时对称性最高,绝对值越大对称性越差。当目标圆心在x轴上时,这7组电势差的理论值为零,由于计算时的离散化会产生一定的误差,该误差精度可用θ=0时测到的电势差最大值来确定,本文算得该精度值为1.32‰。
3. 用对称性分析确定目标小圆圆心的区域
确定目标小圆圆心区域需要三步,第一步确定在哪个π/2区域,第二步确定在哪个π/4区域,第三步确定在哪个π/8区域,确定的过程见图2。红色箭头数字小的一个表示电流注入方向,下面叙述具体过程。
3.1. 判断目标圆心在哪个π/2区域
电流加在电极1→9之间,见图2(a),由于假定目标圆心的电导率小于大圆的电导率,根据欧姆定律和圆的对称性可知,当目标圆心在I或IV区域时V5-V9的差值大于V1-V5的值,若在II或III区域时则情况相反,所以根据(V5-V9) − (V1-V5)的差值就能确定目标圆心在哪个半圆面。等势线分析可知V5-V13 > 0时,目标圆心在I或III区域;V5-V13 < 0时,目标圆心在II或IV区域。两者结合可确定目标圆心在哪个π/2区域,判断逻辑见表1。
假定测得小圆圆心在第I个π/2区域。
3.2. 判断目标圆心在哪个π/4区域
电流加在电极3→11之间,见图2(b),根据(V7-V11) − (V3-V7)的值确定小圆圆心在哪个π/4区域,由于有3.1节的结论,故不需计算V7-V15的结果,判断逻辑见表2。
假定测得小圆圆心在(0, π/4)区域。
3.3. 判断目标圆心在哪个π/8区域
电流加在电极4→12之间,见图2(c),根据(V8-V12) − (V4-V8)的值确定小圆圆心在哪个π/8区域,判断逻辑见表3。

Figure 1. Figure of electrical impedance tomography with single target
图1. 单目标电阻抗成像示意图


(a) (b) (c)
Figure 2. Sketch map of determining the target’s center area, (a)-(c) are used to determine the target’s center is in which π/2, π/4, π/8 area
图2. 确定目标圆心区域示意图,(a)~(c)分别用于确定目标圆心在那个π/2,π/4,π/8区域

Table 1. The logic table of determining target center is in which π/2 area
表1. 确定目标圆心在哪个π/2区域的逻辑表

Table 2. The logic table of determining target center is in which π/4 area
表2. 确定目标圆心在哪个π/4区域的逻辑表

Table 3. The logic table of determining target center is in which π/8 area
表3. 确定目标圆心在哪个π/8区域的逻辑表
假定测得小圆圆心在(0, π/8)。
经过本节的分析,可以确定目标小圆在哪个π/8区域,该结果可极大提升赋初值的准确性,而且获得这个结论仅运用了减法运算和逻辑判断,只要采集了数据,瞬时就能给出结果。
4. 目标小圆半径和圆心坐标对对称性的影响
为了获得小圆半径和圆心坐标对对称性的影响,测量了V2-V16,V3-V15,V4-V14,V5-V13,V6-V12,V7-V11,V8-V10,7个电势差在a = 1,2,3,4;r = 1,2,3,4,5和θ = kπ/72 (k = 0, 1, 2, …, 9)的数值,然后对不同的变量做出了相应的二维图,按简单和复杂两种情况分析了对称性,下面详细讨论。
4.1. 对称性简单变化的情形
4.1.1. 电势差大于零
V2-V16,V3-V15,V4-V14,V5-V13四个电势差值变化趋势类似且均大于零,只是数值大小有差异,表4列出了r = 5时,不同a时四个电势差的最大值。本节以V5-V13为代表说明。
V5-V13随着r,θ,a的增加而增加,即对称性变差,具体情况见图3和图4。a = 1,2时,可用2次多项式拟合,a = 3,4时可用3次多项式拟合,拟合度可达到或超过0.9998,拟合结果见表5,其中θ = 0时,V5-V13的理论值为零,拟合多项式没有意义,但其测量值可用来标记系统的计算精度。
4.1.2. 电势差小于零
V8-V10的值小于零,当a = 1,2,3时,其值近线性单调变化,随θ,r增加对称性变差(V8-V10的绝对值增加);a = 4,r = 1,2,3,4时,V8-V10随θ增加对称性近线性的单调变差,r = 5较特殊,V8-V10的对称性随θ增加近抛物线式的变差,具体见图5和图6。
4.2. 对称性复杂变化情形
4.2.1. V6-V12
对相同的r,不同半径a的小圆,V6-V12对称性变化的趋势相似。r = 1时,V6-V12<0,对称性随a和θ的增大而减小;r = 2时的对称性非常高,此时V6-V12随θ增加先增加后减小;r = 3,4,5时V6-V12 > 0,对称性随r和θ增加而变小,见图7和图8。对不同的θ而言,固定a,V6-V12随r变化在r = 2附近存在一个临界点,r在临界点左侧,V6-V12小于零;r在临界点右侧,V6-V12大于零,见图9和图10。可根据V6-V12大于或小于零,大概判断r的大小,这一信息对成像而言很有价值。

Table 4. The maximal electrical potential difference for different a (r = 5)
表4. 不同a (r = 5)时电势差的最大值

Figure 3. V5-V13 changes with the polar angle variety of target center
图3. V5-V13随小圆圆心极角的变化

Figure 4. V5-V13 changes with the radius vector variety of target center
图4. V5-V13随小圆圆心矢径的变化

Table 5. The polynomial fitting of V5-V13
表5. V5-V13的多项式拟合结果

Figure 5. V8-V10 changes with the polar angle variety of target center
图5. V8-V10随小圆圆心极角的变化

Figure 6. V8-V10 changes with the radius vector variety of target center
图6. V8-V10随小圆圆心矢径的变化

Figure 7. V6-V12 changes with the polar angle variety of target center
图7. V6-V12随小圆圆心极角的变化
4.2.2. V7-V11
V7-V11的变化较复杂,见图11和图12,现分述如下
1) a=1
r = 1,2,3时,对称性随θ的增加变差,r越小对称性越高;r = 4时,V7-V11 < 0,其值随θ的增加而减小,在θ < 6 π/72时的对称性好于r = 1,2,3,在6 π/72 < θ时,对称性好于r = 2,3;不同的a时,r = 5时V7-V11随θ的变化趋势类似,均随着θ的增加先增后减,极值随a的增加向θ增加的方向漂移。
2) a = 2
r = 1,2,3时,对称性随θ的增加变差,r越小对称性越高;r = 4,V7-V11随θ的增加先增加后减

Figure 8. V6-V12 changes with the polar angle variety of target center
图8. V6-V12随小圆圆心极角的变化

Figure 9. V6-V12 changes with the radius vector variety of target center
图9. V6-V12随小圆圆心矢径的变化

Figure 10. V6-V12 changes with the polar angle variety of target center
图10. V6-V12随小圆圆心极角的变化
小,其对称性在θ < 8 π/72时好于r = 1,2,3,在θ = 9 π/72时的对称性好于r = 2和3,θ = π/72,2 π/72,3 π/72时V7-V11 > 0,其它θ时V7-V11 < 0;
3) a = 3
r = 1时的对称性最高;r = 2和r = 3的对称性相近,在θ < 7 π/72时,r = 3的对称性好于r = 2,7 π/72 ≤ θ,r = 3的对称性差于r = 2;r = 4,5时,V7-V11随θ的增加呈抛物型变化,极值点随a的增加向θ增加的方向移动。r = 4时,θ = 8 π/72,9 π/72时V7-V11<0,其它θ时V7-V11 > 0。

Figure 11. V7-V11 changes with the polar angle variety of target center
图11. V7-V11随小圆圆心极角的变化

Figure 12. V7-V11 changes with the radius vector variety of target center
图12. V7-V11随小圆圆心矢径的变化
4) a = 4
r = 1的对称性最高;相同的θ时,r = 3的对称性好于r = 2;r = 4时,V7-V11随θ的增加先增加后减小,其值均大于零。
r = 1,2,3时,V7-V11 < 0;r = 5时,V7-V11 > 0;r = 4时,V7-V11的值有正有负,见图12。
由图12看到r = 4,5时,V7-V11随θ的变化比较复杂。图13是重排数据后画出的不同a时相应的电势差图,表6给出了相应的多项式拟合系数,拟合度均在0.9998及以上。r = 4时,V7-V11 < 0的情况有a = 1;a = 2且3 π/72 < θ;a = 3且7 π/72 < θ。r = 4,a = 2,3,4时,V7-V11随θ的增加先增加后减小,极值随a的增加向θ增加的方向偏移。r = 4时,对a = 2,3,V7-V11均存在一个由正值转负值的临

Figure 13. V7-V11 changes with the polar angle variety of target center
图13. V7-V11随小圆圆心极角的变化

Table 6. The polynomial fitting of V7-V11 at r = 4, 5
表6. V7-V11在r = 4,5时的多项式拟合结果
界θ,这条信息也很有价值。
5. 结论
利用对称性分析,能够较快确定单目标区域中心的范围区间,对16电极可确定目标中心区域的范围为π/8(若是32电极,目标中心区域范围可缩小至π/16)。结合具体情况,还可在较大范围内消除部分不确定性。在a = 1,2,3,4;r = 1,2,3,4,5和θ = kπ/72 (k = 0,1,2,…,9) 条件下有:1) V2-16,V3-V15,V4-V14,V5-V13均大于零,它们对相同的参量变化趋势类似;2) V8-V10小于零,其值在a = 4和r = 5时随θ的变化以抛物线的方式减小,在其它a,r值时随θ的变化趋势变化类似,均以近线性的方式减小,(θ, r, a)越大,V8-V10越小;3) V6-V12,V7-V11比较特殊,对不同的a,r有较大差异,存在电势差值正负转变的临界点,是消除不确定性重点考虑的因素;4)对称点电势差V(θ, r, a)的单变量情形可用2次或3次多项式拟合。
本文的计算是在均匀背景、目标区域规则的条件下得到,r和a的取值不大,忽略边界效应的影响,该方法经过适当修改可处理非均匀背景的情况。对于定性判定“有或没有”有绝对的参考价值。另外对称性分析对于背景中心区域的目标成像有重要意义,因为目标接近背景中心时,对称点间的电势差会趋近于零,结合对称性分析能够很好地改进那些不能对中心区域成像的算法。
下一步的工作:1) 寻找V(θ, r, a)的拟合通式,由于单变量均可拟合,在理论上应该存在一个多变量拟合的结果,这样会极大促进定量描述目标圆心位矢和半径;2) 寻找对称点电势差正负转换的极值点;3) 改进算法以适合非均匀背景。
基金项目
国家自然科学基金(Grant Nos. 11275099,11475135)。