基于SOA优化的电磁层析成像图像重建研究

王亚东，刘 泽，霍继伟，赵鹏飞，王成飞，袁 伟

(1.北京交通大学 电子信息工程学院，北京 100044;2.中国铁道科学研究院集团有限公司 标准计量研究所，北京 100081)

0 引言

电磁层析成像(EMT,electromagnetic tomography)是一种近二、三十年发展起来的基于电磁感应原理的过程层析成像技术[1]，其目的是研究具有电磁特性的物质在空间的分布。典型的EMT系统主要由4个部分组成：传感器阵列(即检测激励线圈)、激励信号产生和分配单元、数据采集和调制单元以及图像重建和信息获取单元。通过对激励线圈通入交变电流，在被测物场中产生激励磁场；当被测物场内有具有导磁或导电的物体放入时，原有磁场分布将会改变，线圈之间的感性耦合也将随之改变[2-3]。通过切换不同的激励线圈，可获得多个投影方向的检测信息，再运用合适的图像重建算法就可以重建出放入物场中物体的电导率或者磁导率的分布图，从而获得放入物体的信息。由此可见，EMT具有非侵入、非接触和快速性的特点，目前在生物医学、无损检测、钢轨探伤、异物检测和工业炼钢等领域均已展开相应研究，具有很好的应用以及发展前景[4-6]。

在EMT系统中，图像重建算法是EMT技术研究中的重点和难点，图像重建算法的好坏将直接影响最终成像的准确性。目前的图像重建算法主要有线性反投影(LBP)[7]、Tikhonov正则化[8]、Landweber迭代[9]等算法。在传统的EMT系统中，边界测量值的微小变化就会导致求解出的物场内的介质分布有很大变化，导致求解过程不稳定，因此EMT技术本身具有病态性，严重影响了重建图像的质量。目前常采用截断奇异值算法和Tikhonov正则化算法来改善灵敏度矩阵的病态性[10-11]，但都存在参数选取问题。Tikhonov正则化算法主要依赖于正则化因子，但是正则化因子的选取较为复杂，在实际应用中往往依据经验进行选取，具有较大的局限性。如果正则化参数选择较小，原始的问题就不能得到很好地近似，逆问题的病态性的改进效果不显著；如果正则化参数选择过大，虽然对误差的敏感性会减少，但最终解会与正常值相差较大，达不到效果。

电磁层析成像的逆问题的病态性的根源是灵敏度矩阵的条件数较大，本文通过利用人群搜索算法(SOA，seeker optimization algorithm)对物场为空时的灵敏度矩阵进行优化，能够明显降低其条件数。将优化过后的灵敏度矩阵应用于Landweber迭代算法进行图像重建，仿真实验结果表明，该方法能够明显提高图像重建的效果以及精度。

1 理论基础

1.1 电磁层析成像数学模型

电磁层析成像的研究包括正问题以及逆问题。正问题是具有电磁特性的物质在敏感场中的分布位置已知，求解其对激励场的调制的结果；逆问题则相反，是利用在正问题中已经得出的调制结果，来求解具有电磁特性的物质的在物场的分布。在EMT中，其调制结果往往表现为边界测量值，即线圈的感应电压。二者的关系可由图1所示。

图1 正问题与逆问题关系图

EMT正问题中的模型中线圈感应电压可用式(1)描述：

Udec=f(Be(x,y),Bo(x,y),μ(x,y),σ(x,y),ε(x,y))

(1)

式中,Udec为测量电压，Be为激励主磁场，Bo为由涡流产生的次级磁场，σ为电导率，μ为磁导率，ε为介电常数。

EMT的逆问题就是利用已知的灵敏度矩阵S以及边界测量值U来解出灰度矩阵G，再根据所求出的灰度矩阵G即可对灵敏场中被测物体进行图像重建。由于检测线圈获得的边界测量值与物场空间之间的关系是非线性的，因此EMT的逆问题可近似归结为线性方程组式(2)的求解[12]。

U=S·G

(2)

式中,U∈RM为归一化之后的边界测量值矩阵，S∈RM×N为归一化后的灵敏度矩阵，G∈RN为归一化之后的灰度矩阵。M为独立测量值的个数，N是逆问题剖分单元数。

1.2 EMT灵敏度矩阵分析

在EMT逆问题的成像计算中，灵敏度矩阵作为先验信息，反应了特定的激励检测对于不同位置单元介质变化的敏感程度，是图像重建的基础。在实际应用中，获取灵敏度矩阵主要有3种方法：模型扰动法[13]、测量扰动法[14]和场量提取法[15]。因模型扰动法及测量扰动法计算量大且较为繁琐，因此在本文中，将采用简便易行，计算速度快的场量提取法。

采用场量提取法进行灵敏度矩阵的求取并画出典型激励检测组合的灵敏度分布如图2所示。

图2 典型灵敏度矩阵分布图

1.3 人群搜索算法应用分析

电磁层析成像中的图像重建属于第一类Fredholm积分方程问题，其解具有不适定性，这主要表现在图像重建的过程中灵敏度矩阵的病态性。条件数可以用来衡量灵敏度矩阵的病态程度，因此，可以采取一定的算法令灵敏度矩阵的条件数达到最小，从而改善其病态性，即转化为一个优化问题。

人群搜索算法[16]能够直接模拟人的智能搜索行为，立足传统的直接搜索算法，以搜寻队伍为种群，搜寻者位置为候选解，通过模拟人类搜寻“经验梯度”和不确定推理进行搜索方向和搜索步长的更新，并不断更新位置，完成对问题的最优求解。

人群搜索算法主要包含搜索步长、搜索方向以及个体位置的更新[17]，为了简化模型，提高运算效率，在本文中，对目标函数值进行排序时采用线性隶属函数，定义如式(3)和式(4)所示：

(3)

Uij=Ui+(1-Ui)·rand,j=1,2,,D

(4)

式中,Ui是目标函数i的隶属度；Uij是j维搜索空间目标函数值i的隶属度；P为种群大小；D为搜索空间维数；Ki为种群函数值按照降序排列之后的序列编号；rand()能够生成0～1之间的随机数，用来模仿人的随机行为。

1.4 Landweber迭代算法

由于边界测量值与敏感场中的电导率分布是非线性的关系，因此通过任何非迭代的算法都几乎不可能找到精确解。目前Landweber迭代算法具有良好的成像质量以及精度，是目前应用最广泛的图像重建算法之一。

从原理上来讲，Landweber算法是最速下降法的一种变形，其本质是通过迭代来计算灵敏度矩阵的广义逆矩阵[18]。一般地，Landweber迭代法更多的用来解决Fredholm积分方程的典型病态问题。

本文将采用LBP重建的图像作为迭代的初值，Landweber迭代算法如式(5)所示：

(5)

式中,β是增益因子，用来控制收敛速度，gk+1为第k步迭代的图像灰度值，g0为利用线性反投影算法计算得到的灰度初始值，S为灵敏度矩阵，U为测量电压值。

2 基于人群搜索算法优化的EMT图像重建

基于人群搜索算法优化的电磁层析图像重建的流程主要分为两个部分：第一部分是利用人群搜索算法对利用场量提取法得到的灵敏度矩阵进行预处理，降低其条件数；第二部分则是根据优化之后灵敏度矩阵利用图像重建算法进行成像。针对第二部分，本文所采用图像重建算法为Landweber迭代算法，为了便于表示，采用此方式进行图像重建结果的算法命名为O_Landweber迭代算法，流程如下：

1)参数的初始化，在人群搜索算法中，设定种群规模为sizepop，迭代次数为maxgen，空间维数为D，最大隶属度值Umax，最小隶属度值Umin，权重最大值Wmax，权重最小值Wmin，种群个体最大值为popmax，种群个体最小值为popmin，EMT逆问题剖分网格数为N，Landweber迭代算法的迭代次数为V。

2)按照1.2节的方法计算灵敏度矩阵S，维度为D*N。

3)按照参数初始化中定义的种群中个体的最大值和最小值，完全随机化生成初始种群T，维度为1*D。

4)根据人群搜索算法以适应度值最小为最优进行搜索的特点[19]，对文献[19]中提出的适应度值的求取公式进行取负操作，并且以矩阵T中的元素为对角线元素生成对角阵T1，然后计算种群中的个体的适应度值。定义如式(6)和式(7)：

T1=diag(T)

(6)

(7)

其中:P是惩罚函数，定义如式(8)所示：

(8)

式中,cond(T1)表示T1的条件数。

5)按照1.3节的方法对搜索步长、搜索方向进行更新。

6)如果达到迭代次数，进行下一步，否则返回第(4)步。

7)我们依据所改善的灵敏度矩阵的条件数占初始灵敏度条件数的比例来最终衡量灵敏度矩阵的优化效果，定义如式(9)所示，delta越大表明优化效果越好。在人群搜索算法中，本文采取双停止策略，即若连续循环3次delta的变化量都不超过0.1%，或者已经循环优化了20次，则进入(8)，否则进入(4)继续更新T。

(9)

式中,T1best表示该次循环后种群中的最佳个体。

8)在得到最终优化的T1best之后，根据式(10)更新灵敏度矩阵以及测量电压。本文将按照线性反投影算法得到的初值进行Landweber迭代算法的计算，将更新后的灵敏度矩阵以及测量电压应用到式(5)所定义的Landweber迭代算法中进行迭代，直到达到迭代次数V时 停止，进行图像重建。

(10)

式中,U1和S1分别表示优化之后的测量电压以及灵敏度矩阵。

基于SOA优化的EMT图像重建的流程图如图3所示。

图3 基于SOA优化的EMT图像重建流程图

3 仿真实验与分析

为了验证本文所提出的利用人群搜索算法改善灵敏度矩阵的病态程度，然后进行图像重建的有效性，主要选取了6种物场分布，采用传统的8线圈EMT系统结构进行仿真[20-21]。为了比较不同算法的仿真效果，本文采用LBP，Tikhonov正则化，Landweber迭代算法以及第2节提出的O_Landweber迭代算法进行成像。

3.1 仿真实验环境设定

与二维仿真相比，三维仿真更能够模拟实际空间的复杂情况，准确地表达出各物体之间相互关系，得到的仿真结果也就更为准确，故本文采用电磁场有限元仿真软件对仿真模型进行三维建模。利用Matlab来控制在仿真过程中激励线圈的切换，激励电流的设置以及敏感场中物理量的获取。

在建模时设定敏感场半径为50 mm，线圈半径为16 mm，激励频率为200 kHz，激励电流为10 A，线圈材质为铜。在进行实验时，所放入的圆柱形物质为半径为8 mm的铜棒，其余为空气，本文所建立的模型如图4所示。

图4 三维仿真模型

设定电磁层析成像的剖分单元格数N为322，如图5所示；在人群搜索算法中，设定种群规模为sizepop为100，迭代次数为maxgen为50，空间维数D为28，最大隶属度值Umax为0.950 0，最小隶属度值Umin为0.011 1，权重最大值Wmax为0.9，权重最小值Wmin为0.1，种群个体最大值为popmax为1，种群个体最小值popmin为0；设定Landweber迭代算法的迭代次数V为1 000次[22]。

图5 逆问题网格剖分结果

3.2 仿真实验结果与分析

本文的仿真是在Matlab以及电磁场有限元仿真软件中进行的，利用Matlab来控制有限元仿真软件的仿真进程、仿真数据获取以及图像重建。仿真过程主要分为3个部分：第一部分为灵敏度矩阵的计算，第二部分为边界测量值获取，第三部分为图像重建。首先获取灵敏度矩阵，计算每个线圈在单位电流激励下在被测物场中产生的磁场，然后确定激励和检测线圈的位置，从场量中提取灵敏度[15]，从而得到整个物场的灵敏度矩阵；获取边界测量值则采用单线圈激励的策略，在电磁场有限元仿真软件中，对3.1节中建立的模型中的8个铜线圈依次激励，共可得到8×8个测量值；最后依据灵敏度矩阵和边界测量值采用不同的图像重建算法在MATLAB中进行图像重建。

本文通过采用循环利用人群搜索算法来改善灵敏度矩阵的病态程度，适应度值和灵敏度矩阵条件数的收敛曲线如图5所示，本次实验共循环了11次，每次循环次数maxgen为50，为了将所有的循环过程都呈现在一副图片上，设定横坐标为总迭代次数，共迭代11*50=550代，纵坐标为每次迭代所对应的适应度值以及灵敏度矩阵条件数。

图6 适应度值和灵敏度矩阵条件数收敛曲线

由图6可知，根据适应度值收敛曲线可以看出，迭代次数小于300时，随着迭代次数的增加适应度值下降明显，当迭代次数大于300时，适应度值下降趋于平缓，逐渐收敛，其表现为灵敏度矩阵条件数也趋于收敛。根据灵敏度矩阵条件数收敛曲线可以看出，用人群搜索算法对灵敏度矩阵进行优化，能够有效降低其条件数。在迭代次数小于100时灵敏度矩阵条件数下降明显，当迭代次数大于300时，适应度值趋于稳定，与适应度值收敛曲线表现基本一致。

灵敏度矩阵条件数的改善量收敛曲线如图7所示。

图7 灵敏度矩阵条件数改善量收敛曲线

由图7可知，在人群搜索算法循环迭代的过程中，灵敏度矩阵条件数改善量随循环次数的增加不断增加；在增加的过程中，曲线的斜率逐渐减小，表明灵敏度矩阵条件数的改善速率不断减小，数值差异也不断缩小，即可以停止迭代优化。结合图6来看灵敏度矩阵的最初的条件数为632.621 8，优化之后降至77.009 9，共降低了87.83%，有效地改善了灵敏度矩阵的病态程度仿真得到的图像重建结果如表1所示。

表1 仿真图像重建结果

如表1所示，第一列为仿真的物体的真实分布图像，第二列为采用LBP算法重建所得到的图像，第三列为采用Tikhonov正则化算法重建所得到的图像，第四列为采用Landweber迭代算法重建所得到的图像，第五列为本文提出O_Landweber迭代算法重建所得到的图像。相对于LBP算法以及Tikhonov正则化算法，通过Landweber迭代算法以及本文所提出O_Landweber算法重建所得到的图像能够较为准确的重建出物体所在的位置，并且能够减少伪影。由表中第四列Landweber算法重建的图像可知，在物体周围存在较多的伪影，物体的位置还不够精确，本文提出的O_Landweber迭代算法能够进一步减少物体周围的伪影，还原物体的实际位置。

为了以定量的方式来对以上4种算法的成像质量以及成像精度进行比较，本文采用图像相对误差(IE,image error)以及相关系数(CC,correlation coefficient)作为评价指标进行比较[23]。图像相关系数是衡量变量之间线性相关程度的指标；图像误差是重建图像与原始图像之间的差异程度。若将采用两种不同成像算法同一种物场分布的仿真结果进行比较，如果其中一种仿真图像的图像相对误差更小，而相关系数更大，那么就说明该仿真结果的图像重建的质量最好。定义如式(11)和(12)所示：

(11)

(12)

表2和表3给出了6中不同分布下4种不同的算法所得到的图像重建结果的图像误差以及图像相关系数。

表2 不同图像重建算法对各类型分布的图像相关系数

表3 不同图像重建算法对各典型分布的图像误差模型

由表2、表3可知，当被测物场中含有一个物体时，本文提出的O_Landweber迭代算法相比Landweber迭代算法的图像相关系数能够提升22%以上，图像误差能够降低56%以上；当被测物场中含有两个物体时，不论是图像相关系数还是图像误差，改善幅度均没有一个物体时明显，但仍有提升，图像相关系数能够提升8%以上，图像误差能够降低28%以上。本文提出的O_Landweber迭代算法相对于LBP、Tikhonov正则化以及Landweber迭代算法的图像相关系数更大，图像误差更小，这表明其所成图像与原始分布的图像误差较小且线性相关程度较高，与仿真图像重建结果表现的基本一致。

4 结束语

由于EMT自身特性，导致灵敏度矩阵往往都具有病态性。为了进一步提高电磁层析成像的成像质量，本文提出先通过人群算法降低灵敏度矩阵病态程度，然后再进行图像重建来提高成像质量。通过循环利用人群搜索算法对灵敏度矩阵进行不断优化的方法，使灵敏度矩阵的条件数能够在迭代300代之后趋于收敛，相对于初始灵敏度矩阵条件数能够降低87.83%，灵敏度矩阵病态程度得到明显改善。仿真实验表明：与传统的直接用于成像的电磁层析成像算法相比，本文所提出的O_Landweber迭代算法能够将图像相关系数最高提升45%，降低图像误差最高75%；与LBP算法、Tikhonov正则化算法，Landweber迭代算法相比，本文所提出的算法能够明显消除伪影，提高成像质量，进一步精确物体的位置。