电学层析成像(Electrical Tomography，ET)技术作为一种高新检测技术，包括电阻层析成像(Electrical Resistance Tomography，ERT)^[1-2]、电磁层析成像(Electromagnetic Tomography，EMT)^[3-4]、电容层析成像(Electromagnetic Capacitance Tomography，ECT)^[5-6]、电阻抗层析成像(Electromagnetic Impedance Tomography，EIT)^[7-8]。其中，电磁层析成像技术是一种基于电磁感应原理的层析成像技术，具有非接触、非介入等特点，在多相流检测、无损检测、生物医学等领域得到了广泛研究^{[3-4, 9-10]}。
对于一个电磁层析成像系统，影响重建图像质量的因素包括激励信号频率、传感器结构、图像重建算法等。在众多因素中，图像重建算法是研究的热点之一^[11-12]。目前，常用的电磁层析成像图像重建算法可分为2类：一类是快速成像算法，包括线性反投影算法(Linear Back-Projection, LBP)、Tikhonov正则化算法等^[13]，该类算法图像重建速度较快，但重建的图像质量较差；另一类是迭代成像算法，包括Landweber迭代算法、代数迭代算法(Algebraic Reconstruction Technique, ART)、Tikhonov迭代算法等^[14]，该类算法图像重建质量较高，但耗时较长。
为了进一步提高Landweber迭代算法的成像质量与成像速度，本文提出了一种优化Landweber迭代快速图像重建算法。首先, 简要讨论了电磁层析成像的逆问题，给出了逆问题的求解公式。其次，讨论了本文所用灵敏度矩阵的计算方法，并给出了典型的灵敏度分布图。然后，给出了本文算法的实现方法。最后，通过仿真实验，从成像质量和成像速度上分别对本文算法与传统算法进行了对比。
1 电磁层析成像逆问题描述 电磁层析成像的逆问题是通过测得的感应线圈电压求解被测区域电磁特性物质的分布。若将敏感场看作为二维场，则感应线圈电压V与物场分布f间的关系可表示为^[15]

(1)

式中：μ(x, y)和σ(x, y)分别为被测区域的磁导率和电导率分布；B为与物场分布有关的系统特性函数。
当被测物场离散化后，式(1)可近似表示为线性矩阵形式：

(2)

式中：U为测量电压向量；S为灵敏度矩阵；g为待求物场电导率(磁导率)分布的灰度值矩阵。
2 灵敏度矩阵计算 由第1节讨论可知，灵敏度矩阵是图像重建的基础。灵敏度矩阵被定义为检测线圈的电压变化和引起检测线圈电压变化的单一像素单元的电磁参数变化之比，可表示为^[16]

(3)

式中：ΔV为检测线圈电压的变化量；I_exc为激励线圈电流；Δκ为复电导率的变化量。
本文采用场量提取法，通过电磁场仿真软件对灵敏度矩阵进行计算，可表示为^[17]

(4)

式中：ω为激励信号频率；H_A为线圈A作为激励线圈时，被测区域的磁场强度分布；H_B为线圈B作为激励线圈时，被测区域的磁场强度分布。
采用该方法所计算的典型的灵敏度矩阵如图 1所示。

图 1 典型灵敏度矩阵分布 Fig. 1 Typical sensitivity matrix distribution

图选项

3 优化Landweber迭代快速图像重建算法 3.1 灵敏度矩阵降维 灵敏度矩阵S具有很强的行相关性^[18]。通过矩阵S的协方差矩阵，将各列数据变换到新坐标下，保留大方差方向信息，可去除行相关性，降低矩阵维度，从而减少成像计算量，进而提高成像速度。
设灵敏度矩阵S表示为

(5)

式中：s_ij为灵敏度矩阵S中元素。
灵敏度矩阵均值矩阵为

(6)

式中：。
由此，可得矩阵X表示为

(7)

式中：。
矩阵X协方差矩阵C_X可表示为

(8)

设矩阵X可通过变换矩阵P降维，则降维后得到的新灵敏度矩阵Y可表示为

(9)

则Y的协方差矩阵C_Y可表示为

(10)

从而C_X可表示为

(11)

使得C_Y满足:

(12)

由于C_X为实对称矩阵，可进行对角化：

(13)

式中：V=[v₁, v₂, …, v_m]，满足条件:

(14)

其中：v_i为C_X的特征向量；λ_i为v_i对应的特征值。将需要保留的k个特征值从大到小按递减顺序排列，将相对应的特征向量组合成矩阵: V′=[v₁, v₂, …, v_k]，通过式(15)求得降维后的灵敏度矩阵：

(15)

对等式(2)两边做同样运算以保证等式成立，即

(16)

进而采用已有的成像算法估计g。
图 2给出了原始灵敏度矩阵的协方差矩阵的特征值。可知，排序后的后8个特征值较小，表明该行数据对灵敏度矩阵的表征不明显，故本文保留特征值较大的20个特征向量构成变换矩阵V′，经过变换矩阵V′降维后，灵敏度矩阵维度由28×322降至20×322，降低了激励检测的组合数，灵敏度矩阵的条件数由610降至158。

图 2 原始协方差矩阵特征值 Fig. 2 Eigenvalues of original covariance matrix

图选项

因降维后的灵敏度矩阵存在部分数据缺失，若直接采用降维后的数据进行图像重建，成像质量难以提升，在某些分布情况下，成像质量会降低。故需对降维后的灵敏度矩阵采用人群搜索算法(Seeker Optimization Algorithm, SOA)进行优化，从而提高成像质量。
3.2 人群搜索算法降低灵敏度矩阵条件数 为了进一步改善灵敏度矩阵的不适定性，本文采用人群搜索算法在降维灵敏度矩阵基础上进一步降低灵敏度矩阵的条件数，从而提高成像质量。人群搜索算法是通过建模人的利己、利他、预动等行为特征，计算优化问题的搜索步长和方向，从而得到最优解的智能优化方法^[19]，具体流程如下：
步骤1 ??设定初始化参数，包括种群规模、迭代次数、种群个体最小值、种群个体最大值、最小隶属度值、最大隶属度值、权重最大值及权重最小值。
步骤2 ??根据初始化的种群个体最大值、种群个体最小值，随机生成初始种群F。
步骤3?? 以矩阵F中的元素作为对角元素构成对角阵M，由式(17)确定适应度值：

(17)

式中：M为以矩阵F中的元素作为对角元素构成对角阵；q为适应度值；cond为求取条件数；P为惩罚函数，表示为

(18)

步骤4?? 采用线性隶属函数，确定种群i在j维搜索空间的步长。线性隶属函数可表示为

(19)

式中：U_i为种群i的隶属度；U_ij为j维搜索空间种群i的隶属度；D为搜索空间维数；N为种群大小；rand为随机数；K_i为种群函数值按降序排列后的序列编号。
种群i在j维搜索空间的步长可表示为

(20)

式中：δ_ij为高斯隶属函数参数，可由式(21)确定：

(21)

其中：x_min和x_max分别为同一种群中最小和最大适应度函数值的位置；b为惯性权值，可由式(22)确定：

(22)

式中：gen为当前迭代次数；gen_max为最大迭代次数。
步骤5 ??通过式(23)确定种群中每个个体搜索方向：

(23)

式中：φ₁、φ₂为区间[0, 1]内任意的实数；sign为符号函数；d_i(t)为第i个搜索因子的搜索方向；d_{i, ego}(t)为第i个搜索因子的利己方向；d_{i, alt}(t)为第i个搜索因子的利他方向；d_{i, pro}(t)为第i个搜索因子的预动方向。

(24)

式中：p_{i, best}为第i个搜索因子的历史最优位置；g_{i, best}为第i个搜索因子所在邻域内的历史最优位置；x_i(t₁)、x_i(t₂)分别为{x_i(t-2)，x_i(t-1)，x_i(t)}中的最优点。
步骤6?? 通过式(25)获得新的位置：

(25)

步骤7?? 若不满足迭代次数，返回步骤3；若满足，通过式(26)计算e，若e连续3次变化不超过0.1%或优化循环20次，进入步骤8，否则返回步骤3。

(26)

式中：M_best为该次循环后种群的最佳个体。
步骤8?? 对式(16)两边做相同的矩阵变换，保证等式成立：

(27)

图 3给出了灵敏度矩阵的条件数的变化曲线。可知，随着迭代次数的增加，条件数逐渐下降。优化前的灵敏度条件数为158，优化后降为38，有效改善了灵敏度矩阵的病态程度。

图 3 灵敏度矩阵条件数收敛曲线 Fig. 3 Convergence curve of condition number of sensitivity matrix

图选项

3.3 Landweber迭代成像 Landweber迭代算法是最速下降算法的一种变形^[14]，其迭代形式为

(28)

式中：g₀为迭代初值；β为增益因子；g_k为第k次迭代的图像灰度值。
Landweber迭代算法的增益因子可由其收敛判据确定：

(29)

本文的增益因子β取为2/(λ_max+λ_min)，λ_max和λ_min分别为矩阵S^TS的最大特征值和最小特征值。
4 仿真实验结果与分析 4.1 仿真实验环境 本文仿真环境为Intel Core i3-3260 3.3 GHz CPU，8 GB RAM。采用电磁场仿真软件对模型进行三维建模，线圈数目为8，线圈直径设为16 mm，激励频率设为200 kHz，圆形检测区域半径为50 mm，激励电流强度为1 A，检测物体为铜棒，检测区域剖分为322个单元。
4.2 无噪声情况下成像分析 为验证本文给出的成像算法的成像效果，设置如表 1第2列所示的10种典型分布用于仿真实验。采用不同成像算法(LBP算法、Tikhonov正则化算法、Landweber迭代算法及本文提出的优化Landweber迭代快速图像重建算法)的成像结果如表 1所示。本文采用后验策略确定Tikhonov算法的正则化因子和Landweber迭代算法的迭代次数。在计算正则化因子时，通过已知模型计算不同正则化因子情况下成像质量，选取最优。在确定Landweber迭代算法的迭代次数时，通过已知模型计算不同迭代次数情况下成像质量，选取最优。在表 1中，分布1~8中圆柱的半径为8 mm；分布9中位于左侧的圆柱直径为10 mm，位于右侧的圆柱直径为15 mm；分布10中为三角形柱体。
表 1 无噪声情况下图像重建结果 Table 1 Results of image reconstruction without noise


表选项






为了定量描述各成像算法的成像质量，本文采用图像误差(Image Error，IE)和相关系数(Correlation Coefficient，CC)作为评价指标。图像误差定义为

(30)

相关系数定义为

(31)

式中：g为设定的电导率分布；为计算得到的电导率分布；为中元素；为g的均值；为的均值；L为剖分单元个数。
表 2和表 3分别给出了10种物场分布情况下，不考虑噪声时，采用不同成像算法重建图像的相关系数和图像误差。
表 2 无噪声情况下不同成像算法相关系数比较 Table 2 Comparison of correlation coefficients among different imaging algorithms without noise

模型序号	相关系数
	LBP算法	Tikhonov算法	Landweber迭代算法	本文算法
1	0.166 5	0.394 5	0.719 1	0.810 9
2	0.066 1	0.354 5	0.732 8	0.868 9
3	0.112 6	0.385 8	0.719 2	0.754 9
4	0.114 7	0.354 4	0.690 2	0.720 1
5	0.153 1	0.433 1	0.693 2	0.750 7
6	0.075 7	0.146 3	0.435 6	0.794 8
7	0.112 1	0.289 8	0.575 3	0.788 8
8	0.148 2	0.178 1	0.568 1	0.701 8
9	0.239 3	0.279 3	0.539 6	0.642 8
10	0.378 2	0.569 7	0.687 1	0.730 7

表选项






表 3 无噪声情况下不同成像算法图像误差比较 Table 3 Comparison of image error among different imaging algorithms without noise

模型序号	图像误差
	LBP算法	Tikhonov算法	Landweber迭代算法	本文算法
1	1.200 2	1.025 6	0.718 6	0.652 2
2	1.172 5	0.961 8	0.676 1	0.535 1
3	1.195 8	1.034 2	0.792 3	0.756 4
4	1.205 9	1.001 6	0.782 9	0.765 7
5	1.340 6	1.157 3	0.949 7	0.792 5
6	1.071 8	0.974 7	0.894 1	0.608 7
7	1.399 6	0.901 1	0.872 4	0.858 9
8	1.279 4	0.956 3	0.810 4	0.725 2
9	1.589 4	0.907 4	0.850 9	0.862 4
10	0.964 2	0.911 1	0.936 6	0.788 0

表选项






由表 1~表 3可见，对于不同物场分布，本文算法的图像误差均小于其他算法，相关系数均大于其他算法，即用本文给出的图像重建算法所重建的图像与原图像的吻合最好，能够较为准确地反映物体位置。
4.3 噪声情况下成像分析 为了研究本文算法的抗干扰性能，本文对测量信号加入干扰，研究算法在不同分布下的成像效果。对检测的物场信号和空场信号分别加入干扰：

(32)

(33)

式中；U_f为测量的物场向量；U_e为测量的空场向量；min()用于求取向量的最小值，控制干扰量级；rand()用于获取随机向量；α为设定值，用于控制干扰大小。当参数α=0.03时，重建图像如表 4所示。表 5和表 6分别给出了存在噪声时不同物场分布情况下各成像算法重建图像的相关系数和图像误差。
表 4 噪声情况下图像重建结果 Table 4 Results of image reconstruction with noise


表选项






表 5 噪声情况下不同成像算法相关系数比较 Table 5 Comparison of correlation coefficients among different imaging algorithms with noise

模型序号	相关系数
	LBP算法	Tikhonov算法	Landweber迭代算法	本文算法
1	0.166 8	0.394 2	0.714 3	0.803 1
2	0.066 4	0.351 5	0.728 1	0.862 5
3	0.112 2	0.387 1	0.716 7	0.757 1
4	0.114 1	0.354 7	0.687 2	0.757 1
5	0.152 7	0.433 6	0.695 4	0.756 9
6	0.017 8	0.154 91	0.451 3	0.689 2
7	0.125 9	0.289 2	0.574 3	0.783 3
8	0.141 9	0.176 1	0.565 4	0.691 3
9	0.240 4	0.277 1	0.536 8	0.637 7
10	0.378 0	0.571 3	0.685 1	0.729 4

表选项






表 6 噪声情况下不同成像算法图像误差比较 Table 6 Comparison of image error among different imaging algorithms with noise

模型序号	图像误差
	LBP算法	Tikhonov算法	Landweber迭代算法	本文算法
1	1.200 7	1.021 9	0.725 6	0.684 5
2	1.180 3	0.961 9	0.680 8	0.543 2
3	1.191 0	1.033 7	0.763 0	0.804 4
4	1.208 4	0.999 2	0.784 8	0.749 5
5	1.337 8	1.149 5	0.937 7	0.810 4
6	1.006 8	0.974 1	0.886 6	0.712 0
7	1.381 5	0.901 3	0.874 5	0.883 8
8	1.265 1	0.956 8	0.816 5	0.739 9
9	1.597 5	0.908 1	0.850 8	0.874 9
10	0.963 7	0.909 6	0.941 2	0.782 7

表选项






对比表 1~表 6可见，对于不同物场分布，当检测信号中含有一定噪声时，各算法的成像质量均有所下降，但本文算法的成像质量总体优于其他算法，体现出一定的抗噪声性能。
4.4 算法计算量分析 由第3节讨论可知，本文算法中的灵敏度矩阵降维与人群搜索算法降低灵敏度矩阵条件数均为离线计算，并不影响在线成像的计算时间。本文算法的灵敏度矩阵维度为20×322，Landweber迭代算法的灵敏度矩阵维度为28×322，通过式(28)可确定每次迭代运算的计算量。
若灵敏度矩阵为m×n维，则根据式(28)可知，每次迭代中，计算项U-Sg_k需要m×(n-1)次加法、m×n次乘法、n次减法；计算项g_k+1额外需要n次减法、2(m×n)次乘法、n×(m-1)次加法；降维后灵敏度矩阵变化的是m，故减法数量不变。
表 7给出了本文算法及Landweber迭代算法每次迭代所需的运算量。由表 7可知，本文算法相比于Landweber迭代算法，在每次迭代运算中，加法运算量减少了约29%，乘法运算量减少了约29%。
表 7 两种成像算法运算量比较 Table 7 Comparison of calculation load between two imaging algorithms

算法	加法运算量	减法运算量	乘法运算量
本文算法	12 538	644	19 320
Landweber迭代算法	17 682	644	27 048

表选项






5 结论 本文通过对灵敏度矩阵预处理，给出了一种Landweber迭代快速图像重建算法，提高了电磁层析成像重建图像的重建质量与重建速度。仿真实验结果表明：
1) 本文算法可有效降低灵敏度矩阵条件数，可有效改善灵敏度矩阵的病态性。初始灵敏度矩阵条件数为610，改善后灵敏度矩阵条件数为38，共降低了约94%。
2) 本文算法在一定噪声情况下，成像质量优于传统的LBP算法、Tikhonov正则化算法、Landweber迭代算法。
3) 本文算法可有效减少图像重建运算量，提高成像速度。在每次迭代中，加法的运算量减少约29%，乘法的运算量减少约29%。

参考文献

[1]	李峰, 谭超, 董峰. 全连接深度网络的电学层析成像算法[J]. 工程热物理学报, 2019, 40(7): 1526-1531. LI F, TAN C, DONG F. Fully connected deep network algorithm for electrical tomography[J]. Journal of Engineering Thermophysics, 2019, 40(7): 1526-1531. (in Chinese)
[2]	肖理庆, 王化祥. 基于聚类电阻层析成像静态图像重建算法[J]. 仪器仪表学报, 2016, 37(6): 1258-1266. XIAO L Q, WANG H X. Absolute image reconstruction algorithm based on clustering for ERT[J]. Chinese Journal of Scientific Instrument, 2016, 37(6): 1258-1266. DOI:10.3969/j.issn.0254-3087.2016.06.008 (in Chinese)
[3]	MA X, PEYTON A J, HIGSON S R, et al. Hardware and software design for an electromagnetic induction tomography (EMT) system for high contrast metal process applications[J]. Measurement Science and Technology, 2006, 17(1): 111-118. DOI:10.1088/0957-0233/17/1/018
[4]	WATSON S, WILLIAMS R J, GOUGH W, et al. A magnetic induction tomography system for samples with conductivities below 10 S·m-1[J]. Measurement Science and Technology, 2008, 19(4): 045501. DOI:10.1088/0957-0233/19/4/045501
[5]	HUANG A, CAO Z, SUN S, et al. An agile electrical capacitance tomography system with improved frame rates[J]. IEEE Sensors Journal, 2019, 19(4): 1416-1425. DOI:10.1109/JSEN.2018.2880999
[6]	SUN S, CAO Z, HUANG A, et al. A high-speed digital electrical capacitance tomography system combining digital recursive demodulation and parallel capacitance measurement[J]. IEEE Sensors Journal, 2017, 17(20): 6690-6698. DOI:10.1109/JSEN.2017.2750741
[7]	BROWN B H. Medical impedance tomography and process impedance tomography: A brief review[J]. Measurement Science and Technology, 2008, 19(9): 1-9.
[8]	范文茹, 王勃, 李靓瑶, 等. 基于电阻抗层析成像的CFRP结构损伤检测[J]. 北京航空航天大学学报, 2019, 45(11): 2177-2183. FAN W R, WANG B, LI J Y, et al. Damage detection of CFRP structure based on electrical impedance tomography[J]. Journal of Beijing University of Aeronautics and Astronautics, 2019, 45(11): 2177-2183. (in Chinese)
[9]	曾星星, 何敏, 张健, 等. EMT用于金属结构裂纹图像重建的仿真研究[J]. 电子测量与仪器学报, 2020, 34(1): 186-192. ZENG X X, HE M, ZHANG J, et al. Simulation study on EMT image reconstruction of metal structure flaw[J]. Journal of Electronic Measurement and Instrumentation, 2020, 34(1): 186-192. (in Chinese)
[10]	LIU Z, LI W, XUE F, et al. Electromagnetic tomography rail defect inspection[J]. IEEE Transactions on Magnetics, 2015, 51(10): 1-7.
[11]	何敏, 柴孟阳. 三种电磁无损检测方法综述[J]. 测控技术, 2012, 31(3): 1-4. HE M, CHAI M Y. Summary of three electromagnetic nondestructive testing methods[J]. Measurement & Control Technology, 2012, 31(3): 1-4. DOI:10.3969/j.issn.1000-8829.2012.03.001 (in Chinese)
[12]	何敏, 刘泽, 董峰, 等. 电磁层析成像技术图像重建的仿真研究[J]. 天津大学学报(自然科学与工程技术版), 2001, 34(4): 435-438. HE M, LIU Z, DONG F, et al. Simulation study on image reconstruction of electromagnetic tomography[J]. Journal of Tianjin University (Science and Technology), 2001, 34(4): 435-438. DOI:10.3969/j.issn.0493-2137.2001.04.005 (in Chinese)
[13]	YANG W Q, SPINK D M, YORK T A, et al. An image-reconstruction algorithm based on Landweber's iteration method for electrical-capacitance tomography[J]. Measurement Science and Technology, 1999, 10(11): 1065-1069. DOI:10.1088/0957-0233/10/11/315
[14]	彭黎辉, 陆耿, 杨五强. 电容成像图像重建算法原理及评价[J]. 清华大学学报(自然科学版), 2004, 44(4): 478-484. PENG L H, LU G, YANG W Q. Image reconstruction algorithms for electrical capacitance tomography: State of the art[J]. Journal of Tsinghua University (Science and Technology), 2004, 44(4): 478-484. DOI:10.3321/j.issn:1000-0054.2004.04.013 (in Chinese)
[15]	何元胜, 何敏, 李杰仁. 电磁层析成像图像重建算法[J]. 电气技术, 2007(4): 43-48. HE Y S, HE M, LI J R. Review about image reconstruction arithmetic for electromagnetic tomography[J]. Electrical Engineering, 2007(4): 43-48. DOI:10.3969/j.issn.1673-3800.2007.04.011 (in Chinese)
[16]	ROSELL J, CASA?AS R, SCHARFETTER H. Sensitivity maps and system requirements for magnetic induction tomography using a planar gradiometer[J]. Physiological Measurement, 2001, 22(1): 121-130. DOI:10.1088/0967-3334/22/1/316
[17]	徐凯, 陈广, 尹武良, 等. 基于场量提取法的电磁层析成像系统的灵敏度推算[J]. 传感技术学报, 2011, 24(4): 543-547. XU K, CHEN G, YIN W L, et al. Sensitivity derivation and calculation of electromagnetic tomography (EMT) sensor based on field value extraction[J]. Chinese Journal of Sensors and Actuators, 2011, 24(4): 543-547. DOI:10.3969/j.issn.1004-1699.2011.04.015 (in Chinese)
[18]	刘泽, 肖君, 刘向龙, 等. 一种电磁层析图像快速重建算法[J]. 北京航空航天大学学报, 2018, 44(8): 1569-1576. LIU Z, XIAO J, LIU X L, et al. An algorithm for fast reconstruction of electromagnetic tomography images[J]. Journal of Beijing University of Aeronautics and Astronautics, 2018, 44(8): 1569-1576. (in Chinese)
[19]	张连强, 王东风. 基于改进人群搜索算法的PID参数优化[J]. 计算机工程与设计, 2016, 37(12): 3389-3393. ZHANG L Q, WANG D F. Optimization of PID parameters based on improved seeker optimization[J]. Computer Engineering and Design, 2016, 37(12): 3389-3393. (in Chinese)