电磁层析成像(Electromagnetic Tomography，EMT)技术已经有几十年的发展历史，随着研究的深入，在多个领域有所建树，其主要应用于现场状况不好条件下需要长期进行工作的稳定监测系统，在地质勘查、异物检测与定位、化工分离过程的监测等领域有了长足进展。
针对EMT逆问题^[1-5]，现有的算法有线性反投影(Linear Back-Projection，LBP)法、Tikhonov正则化法、Landweber迭代算法、共轭梯度法等，并在原有的算法基础上，提出了参数优化、迭代算法的初始值优化、迭代速度优化等多种优化方法^[6-10]。本文通过研究灵敏度矩阵中信息冗余特性，提出先降维再计算广义逆的新思路，在该算法中又提出图像相关系数最大化算法(与图像相对误差最小化算法具有一致性)，确定保留协方差矩阵特征值个数。首先，介绍了灵敏度矩阵降维原理和必要性，以及本文成像算法和参数选取算法原理；然后，通过仿真，从成像质量和成像速度2方面对本文算法和传统算法进行了对比，阐明本文算法的优势和劣势；最后，应用电动平移台和EMT硬件系统对比各算法的实际成像效果。
1 灵敏度矩阵降维原理 对于灵敏度矩阵，电导率灵敏度场的计算可由式(1)得到，感应电压V对电导率σ的灵敏度为

(1)

式中：ω为系统工作频率；Ω为扰动体积；A_i和A_j分别为激励线圈和检测线圈在电流I激励时的矢量磁位分布，最终可得电导率灵敏度图。
灵敏度矩阵的行相关性是造成其病态性、解不稳定的重要原因，因此分析其协方差矩阵，将各列的样本在新坐标基下做空间变换，保留大方差方向的信息，去除行间相关性以提高解稳定性。
通过保留低阶主成分，忽略高阶主成分，达到降低数据维数，同时保持数据集中对方差贡献最大的特征的目的。通过正交线性变换，把原始数据变换到一个新坐标中，使得这一数据的任何投影的第一大方差在第一个坐标上，第二大方差在第二个坐标上。设去均值后的矩阵为X，存在变换矩阵P，X的主成分w₁被定义为

(2)

为了得到第k个主成分，必须先从X中减去前面的k-1个主成分：

(3)

然后把求得的第k个主成分代入数据集，得到新的数据集，继续寻找主成分：

(4)

具体分解过程如下：存在每行元素去均值后的矩阵X_m×n，其协方差矩阵为，通过分析其协方差矩阵，将原矩阵的主要信息保留，去掉噪声信息，并达到降维的效果，使降维后矩阵Y的各行间协方差为零，每行的方差尽可能大，使数据的分散程度较高。式(5)和式(6)分别为零均值化后的灵敏度矩阵及其协方差矩阵。

(5)

(6)

利用矩阵P为X降维，得到新的灵敏度矩阵Y=PX，设新的灵敏度矩阵的协方差矩阵为

将上式转化为

(7)

使cov(Y)满足：

由于cov(X)为实对称矩阵，则可对角化为

(8)

式中：U=[q₁??q₂??…??q_m]，满足关系：

由式(7)、式(8)可以看出：P=U^T，其中q_i为cov(X)的特征向量。分析仿真灵敏度矩阵q_i的协方差的特征值，由图 1可以看出，较大数值的特征值个数较少，表明大部分原灵敏度矩阵对应行的方差较小，该行数据对灵敏场的特性表征作用不明显，通过特征值的截取，使灵敏度矩阵从行与行之间相关性较大、存在单行方差很小的特性，转化为保留方差较大信息，去除无表征作用的信息，小特征值使方程的解g有较高自由度，这也是解不稳定的主要原因。解的不稳定性也是病态矩阵亟须解决的问题。通过保留满足特定要求的特征值，完成新空间的映射，从而达到降维和去噪的目的。

图 1 仿真数据和实验数据的协方差矩阵及其特征值分布 Fig. 1 Covariance matrix and eigenvalue distribution of simulation data and experimental data

图选项

将特征值从大到小排列后，通过保留前k个特征值，得到新的U，即U′=[q₁??q₂??…??q_k]，最终完成Y=(U′)^TX的降维。对于公式z=Sg，其中z为物场添加激励信号后，检测线圈检测到的测量向量，S为灵敏度矩阵，g为设定物场空间电导率和磁导率的分布。对公式两边做同样处理，保证方程成立。
对于去均值的灵敏度矩阵，利用其协方差矩阵降维，被称为主成分分析(Principal Component Analysis，PCA)^[11]，此方法在数据存储、图像压缩、分类问题的特征降维等问题中都有所应用，对均值化后的灵敏度矩阵做协方差分析，其本质和奇异值分解(Singular Value Decomposition，SVD)是具有一致性的。
2 特征值个数选取算法原理 特征值保留个数主要考虑噪声、保留成像特征和解稳定性3方面的影响。对于此问题，有很多研究算法，如L曲线法、广义交叉验证法等。如果保留个数过少，很可能导致虽然噪声减少，但有物体部位的图像成像信息丢失严重，即使是很小的特征值部分，仍未有理论表明其为噪声特性还是正常成像特性；如果保留个数过多，又会导致噪声信息过多，该方向自由度过大，致使解不稳定。因此，三者的权衡十分重要。
本文提出了一种基于多样本特性的图像相关系数最大化算法。设已知m个检测空间，每个检测空间分别放入铜棒并得到仿真数据z和g。然后将协方差矩阵特征值从大到小排列，在保留不同个数特征值条件下分别计算其图像重建质量评价参数：图像相对误差(R_IE)和图像相关系数(R_CC)，其具体计算公式分别为

(9)

(10)

式中：为仿真结果分布；g和分别为g和的平均值; i为保留特征值个数。利用m个z求得m个，分别用g和求出R_IE和R_CC，求出保留不同个数协方差矩阵特征值条件下m个样本累计R_IE和R_CC的平均值，分别取使其均值最小和最大的k值：

(11)

(12)

式中：r为原矩阵的秩。
画出保留不同个数特征值R_IE或R_CC的均值分布，以R_CC为例，如图 2所示。

图 2 多样本的图像相关系数均值分布 Fig. 2 Image correlation coefficient mean distribution with multi-sample

图选项

图 2综合了灵敏度矩阵包含的318个样本，以及下文中测试的5个样本，则仿真中m值为323，将协方差矩阵的特征值从大到小排列，得到保留不同特征值个数条件下所有样本的图像相关系数均值。由图 2可以看出，最大值在30~36之间取得，由于此区段内数据存在波动，考虑尽可能减小成像信息丢失，最终选取了36个最大的特征值对应的特征向量组成的矩阵作为灵敏度矩阵的变换矩阵。由于仿真情况和实际情况噪声来源等多种条件不同，此方法的优点是：直接从优化目标出发，确定参变量大小，在不同外界数据获取条件下更具有适应性，利用灵敏度矩阵中包含的多样本特性，避免偶然性，从而得到特征值最优截取个数。
3 奇异值分解原理 对于实矩阵A，可将其分解为A=FΣV^T，F为m×m维单位正交矩阵，其列向量为AA^T的特征向量，Σ为m×n维半正定对角矩阵，V^T为n×n维单位正交矩阵，其行向量为A^TA的特征向量。


式中：σ₁>σ₂>…>σ_r>0，对于未降维处理的病态灵敏度矩阵，σ₁/σ_r较大，仿真表明，可达到10¹⁰的数量级，矩阵病态程度较高，对于方程Ax=b，当A或b有较小改变，容易引起x的较大变动，从而造成解的不稳定。例如，σ_i很小，解x在v_i方向上有较高自由度，自由度越大，解越不稳定，但对灵敏度矩阵降维后，奇异值分解中不会出现这种情况，对于仿真数据，降维后，σ_max/σ_min < 10⁴，降低了矩阵病态程度，提高了解稳定性。
然后利用SVD求广义逆，S⁺=VΣ⁺F^T，最终得=S⁺z，是对实际灰度值向量的估计，也是最终求得的最优解。
基于SVD广义逆原理^[12-14]和前期灵敏度矩阵降维两步处理，对于仿真过程，降维前后成像效果对比如表 1所示。可以看出，降维去除了中心附近的误成像噪声(仿真原始图见表 2)。
表 1 灵敏度矩阵降维前后多物体成像对比 Table 1 Multi-object imaging contrast before and after dimension reduction of sensitivity matrix


表选项






表 2 各算法多物体仿真成像对比 Table 2 Comparison of multi-algorithm and multi-object simulation imaging


表选项






4 仿真和实验 4.1 仿真 综合考虑仿真验证速度，不同传感器参数对成像质量的影响和硬件实验条件等因素，本文做了318个剖分单元格和8线圈^[15]检测的仿真实验，分别用LBP法、Tikhonov正则化法^[16-17]、Landweber迭代算法^[18-19](Tikhonov做初值，迭代500次，并用自适应时刻估计算法最小化目标函数，以提高收敛速度，其参数分别为?=0.7，β₁=0.9，β₂=0.995)和降维SVD做对比实验，仿真结果如表 2所示。
由表 2~表 4可知，对于成像质量：LBP法存在原理上的缺陷；Tikhonov正则化法相比LBP法有所提高，通过参数优化，仍有上升空间；Landweber迭代算法具有较好的成像质量，在迭代方式和优化函数2个方向也有可研究价值；本文的降维SVD与LBP法、Tikhonov正则化法相比在图像相对误差和图像相关系数2个成像指标上都具有优势，与Landweber迭代算法的成像效果相似。对于耗时(其中时间计算未计算数据加载过程和成像过程，为连续计算灰度值向量1 000次的平均时间，计算资源：Intel(R)Core(TM)i3-3220CPU，3.30 GHz, 计算平台：PyCharm, 计算语言：Python)，由表 5可知：LBP法、Tikhonov正则化法和降维SVD耗时相近，速度较快；降维SVD相比传统算法在时间上仍具有优势，但Landweber迭代算法相对而言耗时较长，这也与步长、学习率、计算停止判断条件等优化指标有关，非迭代算法的先天优势为快速高质量成像提供了可能。
表 3 仿真图像相关系数数据 Table 3 Image correlation coefficient data in simulation

样本类型	LBP法	Tikhonov 正则化法	Landweber 迭代算法	降维 SVD
1物体	-0.010 0	0.506 4	0.787 1	0.801 5
2物体	-0.007 0	0.595 4	0.761 9	0.688 5
3物体	-0.027 4	0.538 6	0.707 9	0.652 5
4物体	-0.040 6	0.424 5	0.644 9	0.606 1
5长物体	-0.026 5	0.279 6	0.319 4	0.326 4

表选项






表 4 仿真图像相对误差数据 Table 4 Image relative error data in simulation

样本类型	LBP法	Tikhonov 正则化法	Landweber 迭代算法	降维 SV
1物体	1.188 7	1.046 7	0.753 2	0.751 1
2物体	1.102 8	0.726 5	0.675 2	0.795 5
3物体	1.042 9	0.729 7	0.637 6	0.668 2
4物体	1.069 5	0.788 3	0.662 2	0.684 0
5长物体	3.422 0	0.886 5	0.868 8	0.875 2

表选项






表 5 各算法计算时间 Table 5 Calculation time of each algorithm

ms
样本类型	LBP法	Tikhonov 正则化法	Landweber 迭代算法	降维 SVD
1物体	1.1	0.4	200.1	0.2
2物体	1.2	0.4	306.4	0.2
3物体	1.1	0.4	440.4	0.2
4物体	1.2	0.3	280.8	0.2
5长物体	1.2	0.4	453.2	0.2

表选项






4.2 实验 本文应用8线圈电动平移台系统(见图 3)和8线圈传感器激励检测系统(见图 4)采集原始数据，分别利用4种算法得到实验成像图。

图 3 电动平移台及其控制器 Fig. 3 Electric moving stage and its controller

图选项

图 4 8线圈传感器激励检测系统 Fig. 4 8-coil sensor excitation and detection system

图选项

实验中，利用如下公式计算实验数据信噪比(Signal-Noise Ratio, SNR)：

(13)

式中：D_noise为某路信号多次采集后求得的平均噪声值；D_signal为某路信号多次采集后求得的平均有用信号值。

连续采集n帧数据x_i，并计算其均值μ和标准差γ，求出某路数据的信噪比SNR，然后利用如下公式, 在8线圈分别激励条件下，求出采集64路数据的平均信噪比：

(14)

实验得：=50.7 dB。
对于本文的特征值选取问题，由于剖分了408个有限元，得到的灵敏度矩阵即为408个训练样本，从而利用图像相关系数最大化算法确定选取特征值个数，避免了实测物体不规则、无法得到实际灰度值向量的问题。由表 6~表 8可以看出，LBP法相对其他算法成像质量较差，改进空间较小；Tikhonov正则化法噪点较多，干扰较强；Landweber迭代算法和本文的降维SVD成像效果相近，仍有改进空间。
表 6 各算法实验成像对比 Table 6 Comparison of multi-algorithm experimental imaging


表选项






表 7 实验图像相关系数数据 Table 7 Image correlation coefficient data in experiment

样本类型	LBP法	Tikhonov 正则化法	Landweber 迭代算法	降维 SVD
1物体	0.531 6	0.612 2	0.652 3	0.647 5
2物体	0.453 4	0.543 3	0.608 4	0.599 0
3物体	0.467 1	0.519 1	0.587 9	0.590 4

表选项






表 8 实验图像相对误差数据 Table 8 Image relative error data in experiment

样本类型	LBP法	Tikhonov 正则化法	Landweber 迭代算法	降维 SVD
1物体	0.960 3	0.787 0	0.674 1	0.685 0
2物体	0.959 7	0.696 6	0.667 2	0.671 1
3物体	0.953 7	0.778 3	0.693 0	0.687 6

表选项






5 结论 本文针对EMT逆问题中，病态矩阵方程求解问题，提出了降维灵敏度矩阵算法，在灵敏度矩阵的协方差矩阵特征值个数选取中，提出图像相关系数最大化算法，并对其原理进行了阐述。仿真和硬件实验结果表明：
1) 相较于传统的LBP法、Tikhonov正则化法，本文算法在成像质量上具有优势。
2) 相较于传统的Landweber迭代算法、LBP法、Tikhonov正则化法，本文算法在成像时间上具有优势。
3) 灵敏度矩阵中多样本的利用，对电磁层析成像各算法的参数选取问题具有一定借鉴意义，可一定程度上避免偶然性。

参考文献

[1]	PEYTON A J, BACK M S, BORGES A R, et al. Development of electromagnetic tomography(EMT)for industrial applications. Part1: Sensor design and instrumentation[C]//1st World Congress on Industrial Process Tomography, 1999: 306-317.
[2]	RAMLI S, PEYTON A J. Feasibility study of planar-array electromagnetic inductance tomography(EMT)[C]//1st World Congress on Industrial Process Tomography, 1999: 502-510.
[3]	王化祥, 尹武良. 用于多相界面检测系统的阵列式容栅电容传感器的优化设计[J].仪器仪表学报, 1996, 17(1): 8–13. WANG H X, YIN W L. Optimum design of an array of segmented capacitance sensor for multi-phase interface measuring system[J].Chinese Journal of Scientific Instrument, 1996, 17(1): 8–13.DOI:10.3321/j.issn:0254-3087.1996.01.002(in Chinese)
[4]	LIU C, XU L J, CHEN J L, et al. Development of a fan-beam TDLAS-based tomographic sensor for rapid imaging of temperature and gas concentration[J].Optics Express, 2015, 23(17): 494–511.
[5]	刘泽, 薛芳其, 杨国银, 等. 电磁层析成像灵敏度矩阵实验测试方法[J].仪器仪表学报, 2013, 34(9): 1982–1988. LIU Z, XUE F Q, YANG G Y, et al. Experimental measurement method of sensitivity matrix for electromagnetic tomography[J].Chinese Journal of Scientific Instrument, 2013, 34(9): 1982–1988.(in Chinese)
[6]	刘泽, 何敏, 徐苓安, 等. 多激励模式的电磁层析成像系统[J].仪器仪表学报, 2001, 22(6): 614–617. LIU Z, HE M, XU L A, et al. Multi-mode excitation electromagnetic tomography (EMT) system[J].Chinese Journal of Scientific Instrument, 2001, 22(6): 614–617.DOI:10.3321/j.issn:0254-3087.2001.06.019(in Chinese)
[7]	彭黎辉, 陆耿, 杨五强. 电容成像图像重建算法原理及评价[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)
[8]	吴杰, 李明峰, 余腾. 测量数据处理中病态矩阵和正则化方法[J].大地测量和地球动力学, 2010, 30(4): 102–105. WU J, LI M F, YU T. Ill matrix and regularization method in surveying data processing[J].Journal of Geodesy and Geodynamics, 2010, 30(4): 102–105.(in Chinese)
[9]	徐驰, 韩磊, 张书第, 等. 最小二乘矩阵滤波器设计与性能分析[J].舰船科学技术, 2011, 33(4): 72–76. XU C, HAN L, ZHANG S D, et al. Design and performance analysis of least-square matrix filter[J].Ship Science and Technology, 2011, 33(4): 72–76.DOI:10.3404/j.issn.1672-7649.2011.04.014(in Chinese)
[10]	王友华. 病态矩阵的本质及其解决方法[J].矿山测量, 2005(2): 62–63. WANG Y H. Essence of ill-conditioned matrix and corresponding solutions[J].Mine Surveying, 2005(2): 62–63.DOI:10.3969/j.issn.1001-358X.2005.02.024(in Chinese)
[11]	阮越, 陈汉武, 刘志昊, 等. 量子主成分分析算法[J].计算机学报, 2014, 37(3): 666–676. RUAN Y, CHEN H W, LIU Z H, et al. Quantum principal component analysis algorithm[J].Chinese Journal of Computers, 2014, 37(3): 666–676.(in Chinese)
[12]	周笋, 吉国力. 基于改进BP算法的广义逆矩阵求解方法[J].厦门大学学报, 2003, 42(3): 386–389. ZHOU S, JI G L. Generalized inverse matrix method based on BP algorithm[J].Journal of Xiamen University, 2003, 42(3): 386–389.DOI:10.3321/j.issn:0438-0479.2003.03.027(in Chinese)
[13]	李高明, 李海鹏. 求矩阵奇异值分解的一种新方法[J].数学实践与认识, 2011, 41(7): 212–215. LI G M, LI H P. A new methods of matrix singular values decomposition[J].Mathematics in Practice and Theory, 2011, 41(7): 212–215.(in Chinese)
[14]	朱晓临, 李雪艳, 邢燕, 等. 基于小波和奇异值分解的图像边缘检测[J].图学学报, 2014, 35(4): 563–570. ZHU X L, LI X Y, XING Y, et al. Image edge detection based on wavelet and singular value decomposition[J].Journal of Graphics, 2014, 35(4): 563–570.DOI:10.3969/j.issn.2095-302X.2014.04.012(in Chinese)
[15]	PENG L, YE J, LU G, et al. Evaluation of effect of number of electrodes in ECT sensors on image quality[J].IEEE Sensors Journal, 2012, 12(5): 1554–1564.
[16]	冯洁婷, 颜刚, 王涛, 等. Tikhonov正则化参数选择对高速率刺激听觉诱发电位重建的影响[J].航天医学与医学工程, 2012, 25(1): 54–60. FENG J T, YAN G, WANG T, et al. Effects of parameter selection of Tikhonov regularization on reconstruction of high-rate auditory evoked potentials[J].Space Medicine & Medical Engineering, 2012, 25(1): 54–60.(in Chinese)
[17]	余瑞艳. 基于混沌粒子群算法的Tikhonov正则化参数选取[J].数学研究, 2011, 44(1): 101–106. YU R Y. Optimal choosing of Tikhonov regularization parameter based on chaos particle swarm optimization algorithm[J].Journal of Mathematical Study, 2011, 44(1): 101–106.(in Chinese)
[18]	唐利民, 朱建军. 基于Landweber迭代的秩亏非线性最小二乘问题算法研究[J].大地测量与地球动力学, 2010, 30(1): 95–98. TANG L M, ZHU J J. Algorithm based on Landweber iteration for solving rank deficiency nonlinear least squares proble[J].Journal of Geodesy and Geodynamics, 2010, 30(1): 95–98.(in Chinese)
[19]	王旭, 王静文, 王柯元. 阈值Landweber在MIT图像重建中的应用[J].东北大学学报(自然科学版), 2016, 37(4): 477–480. WANG X, WANG J W, WANG K Y. Application of threshold-ing Landweber algorithm in MIT image reconstruction[J].Journal of Northeastern University (Natural Science), 2016, 37(4): 477–480.DOI:10.3969/j.issn.1005-3026.2016.04.005(in Chinese)