一种电磁层析图像快速重建算法

刘泽， 肖君， 刘向龙， 赵鹏飞， 李勇， 霍继伟

(北京交通大学 电子信息工程学院, 北京 100044)

电磁层析成像(Electromagnetic Tomography，EMT)技术已经有几十年的发展历史，随着研究的深入，在多个领域有所建树，其主要应用于现场状况不好条件下需要长期进行工作的稳定监测系统，在地质勘查、异物检测与定位、化工分离过程的监测等领域有了长足进展。

针对EMT逆问题[1-5]，现有的算法有线性反投影(Linear Back-Projection，LBP)法、Tikhonov正则化法、Landweber迭代算法、共轭梯度法等，并在原有的算法基础上，提出了参数优化、迭代算法的初始值优化、迭代速度优化等多种优化方法[6-10]。本文通过研究灵敏度矩阵中信息冗余特性，提出先降维再计算广义逆的新思路，在该算法中又提出图像相关系数最大化算法(与图像相对误差最小化算法具有一致性)，确定保留协方差矩阵特征值个数。首先，介绍了灵敏度矩阵降维原理和必要性，以及本文成像算法和参数选取算法原理；然后，通过仿真，从成像质量和成像速度2方面对本文算法和传统算法进行了对比，阐明本文算法的优势和劣势；最后，应用电动平移台和EMT硬件系统对比各算法的实际成像效果。

1 灵敏度矩阵降维原理

对于灵敏度矩阵，电导率灵敏度场的计算可由式(1)得到，感应电压V对电导率σ的灵敏度为

(1)

式中：ω为系统工作频率；Ω为扰动体积；Ai和Aj分别为激励线圈和检测线圈在电流I激励时的矢量磁位分布，最终可得电导率灵敏度图。

灵敏度矩阵的行相关性是造成其病态性、解不稳定的重要原因，因此分析其协方差矩阵，将各列的样本在新坐标基下做空间变换，保留大方差方向的信息，去除行间相关性以提高解稳定性。

通过保留低阶主成分，忽略高阶主成分，达到降低数据维数，同时保持数据集中对方差贡献最大的特征的目的。通过正交线性变换，把原始数据变换到一个新坐标中，使得这一数据的任何投影的第一大方差在第一个坐标上，第二大方差在第二个坐标上。设去均值后的矩阵为X，存在变换矩阵P，X的主成分w1被定义为

(2)

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

(3)

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

(4)

(5)

(6)

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

将上式转化为

cov(X)=PTcov(Y)P

(7)

使cov(Y)满足：

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

(8)

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

对于去均值的灵敏度矩阵，利用其协方差矩阵降维，被称为主成分分析(Principal Component Analysis，PCA)[11]，此方法在数据存储、图像压缩、分类问题的特征降维等问题中都有所应用，对均值化后的灵敏度矩阵做协方差分析，其本质和奇异值分解(Singular Value Decomposition，SVD)是具有一致性的。

2 特征值个数选取算法原理

特征值保留个数主要考虑噪声、保留成像特征和解稳定性3方面的影响。对于此问题，有很多研究算法，如L曲线法、广义交叉验证法等。如果保留个数过少，很可能导致虽然噪声减少，但有物体部位的图像成像信息丢失严重，即使是很小的特征值部分，仍未有理论表明其为噪声特性还是正常成像特性；如果保留个数过多，又会导致噪声信息过多，该方向自由度过大，致使解不稳定。因此，三者的权衡十分重要。

本文提出了一种基于多样本特性的图像相关系数最大化算法。设已知m个检测空间，每个检测空间分别放入铜棒并得到仿真数据z和g。然后将协方差矩阵特征值从大到小排列，在保留不同个数特征值条件下分别计算其图像重建质量评价参数：图像相对误差(RIE)和图像相关系数(RCC)，其具体计算公式分别为

(9)

(10)

(11)

(12)

式中：r为原矩阵的秩。

画出保留不同个数特征值RIE或RCC的均值分布，以RCC为例，如图2所示。

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

3 奇异值分解原理

对于实矩阵A，可将其分解为A=FΣVT，F为m×m维单位正交矩阵，其列向量为AAT的特征向量，Σ为m×n维半正定对角矩阵，VT为n×n维单位正交矩阵，其行向量为ATA的特征向量。

式中：σ1>σ2>…>σr>0，对于未降维处理的病态灵敏度矩阵，σ1/σr较大，仿真表明，可达到1010的数量级，矩阵病态程度较高，对于方程Ax=b，当A或b有较小改变，容易引起x的较大变动，从而造成解的不稳定。例如，σi很小，解x在vi方向上有较高自由度，自由度越大，解越不稳定，但对灵敏度矩阵降维后，奇异值分解中不会出现这种情况，对于仿真数据，降维后，σmax/σmin<104，降低了矩阵病态程度，提高了解稳定性。

基于SVD广义逆原理[12-14]和前期灵敏度矩阵降维两步处理，对于仿真过程，降维前后成像效果对比如表1所示。可以看出，降维去除了中心附近的误成像噪声(仿真原始图见表2)。

表1 灵敏度矩阵降维前后多物体成像对比

表2 各算法多物体仿真成像对比

4 仿真和实验

4.1 仿 真

综合考虑仿真验证速度，不同传感器参数对成像质量的影响和硬件实验条件等因素，本文做了318个剖分单元格和8线圈[15]检测的仿真实验，分别用LBP法、Tikhonov正则化法[16-17]、Landweber迭代算法[18-19](Tikhonov做初值，迭代500次，并用自适应时刻估计算法最小化目标函数，以提高收敛速度，其参数分别为∂=0.7，β1=0.9，β2=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 仿真图像相关系数数据

表4 仿真图像相对误差数据

表5 各算法计算时间

4.2 实 验

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

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

(13)

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

(14)

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

表6 各算法实验成像对比

表7 实验图像相关系数数据

表8 实验图像相对误差数据

5 结 论

本文针对EMT逆问题中，病态矩阵方程求解问题，提出了降维灵敏度矩阵算法，在灵敏度矩阵的协方差矩阵特征值个数选取中，提出图像相关系数最大化算法，并对其原理进行了阐述。仿真和硬件实验结果表明：

1) 相较于传统的LBP法、Tikhonov正则化法，本文算法在成像质量上具有优势。

2) 相较于传统的Landweber迭代算法、LBP法、Tikhonov正则化法，本文算法在成像时间上具有优势。

3) 灵敏度矩阵中多样本的利用，对电磁层析成像各算法的参数选取问题具有一定借鉴意义，可一定程度上避免偶然性。