基于改进的差分进化算法的电阻抗图像重建

刘海，闫丹丹，荆会敏

中国计量大学信息工程学院，浙江杭州 310018；

生物电阻抗断层成像（electrical impedance tomography，EIT）是继形态、结构成像后的新一代无创功能性成像技术。它通过向体表注入较小的安全交变电流，由体表电极阵列测量得到体外边界电压（电势），并采用重建算法重建组织内部电阻率分布[1]。与正常的生物阻抗相比，生物组织在发生病变时，其阻抗值变化十分显著[2]。EIT技术的实质是利用反映组织生理及功能状态的电阻抗信息进行阻抗重建。因此，EIT技术可用于肿瘤、癫痫等多种疾病的早期发现与诊断；并且EIT技术不产生电离辐射、对人体无害、相对成本较低、可多次重复使用、便携，能对患者进行长期、实时动态监护。基于以上特点，尽管目前应用EIT技术与主流医学成像技术在图像分辨率上还有一定差距，但其在医学工程及相关领域仍引起极大的关注[3-4]。

目前主要的重建算法有修正的牛顿-拉夫逊法、牛顿单步误差重构法、层析法、双限定法等，上述方法收敛速度较快，但计算量大且过程复杂。近年来基于粒子群算法、神经网络算法和遗传算法等相关的智能寻优算法逐渐应用于EIT图像重建，这些智能算法与传统算法相比，计算量较小、较易修正局部最优值，具有较快的收敛速度和全局搜索能力等优势[5]。

本文采用改进的差分进化（differential evolution，DE）算法对头部电阻抗分布进行重建，在三维头球模型上进行仿真实验，并与其他进化类算法进行对比。

1 EIT的基本原理

在EIT研究中，电阻抗分布与电流之间的关系由麦克斯韦方程给出。EIT测量中采用较低（10～100 kHz）的激励源频率，可忽略介电常数的影响，此时电流场可当作稳态电流场。为简化计算，忽略电极与皮肤间的接触阻抗，电压与电流间的非线性关系可由满足相应边界条件的拉普拉斯方程确定，其求解场域的数学模型[6-7]见公式（1），应满足的第一类或第二类边界条件（、）见公式（2）。

其中，表示散度，表示梯度，ρ为目标内部电阻率分布，φ为场域中电位分布函数，φ0为边界电位，n为方向向量，为边界电位分布的法向导数，Jn为电流密度，无电流注入时为0。

EIT正问题是指由人体内部阻抗分布及边界条件，求解体表或内部电压分布，这是求解逆问题的基础。EIT逆问题是指由体表电压分布及边界条件，求解人体内部阻抗分布。经过有限元法离散后，EIT逆问题可转化为公式（3）中求电位分布残差的相对均方误差最小值[RMS(ρ)]问题。

其中，上标的T表示矩阵的转置，ρ为目标内部电阻率分布，φ(ρ)是边界电位的计算值，φ0是边界电位“测量值”。

因此，EIT逆问题是找到合适的电阻抗分布ρ使目标函数RMS(ρ)达到最小，即搜寻一个m维参数矢量使RMS(ρ)最小，其中m为有限元剖分单元数，即电阻率分布自由度。从而，EIT逆问题即转化为多参数函数优化问题。

2 基于改进DE的EIT重建算法

2.1 基本 DE的算法 DE算法是一种基于种群进化的全局启发式算法[8]，它在许多优化问题上的效果均优于遗传算法[9]、粒子群算法[10]、模拟退火算法[11]等。该算法一般包括变异、交叉、选择3个基本操作。首先在整个搜索空间随机生成一组初始种群，然后对种群中每个个体进行变异、交叉操作得到临时种群，最后用临时种群中适应性更好的个体替代原始种群中相应的个体。整个算法的步骤见图1。

图1 DE基本算法流程

在整个算法迭代过程中，变异是其中的核心过程，它采用随机偏差扰动生成变异向量的方式能够增加种群多样性，增强算法全局搜索能力。根据生成变异向量方法的不同，能够组合成不同的变异策略，其中基本变异策略DE/rand/1/bin方程见公式（4）。

此外，其他典型的变异策略 DE/best/1/bin和DE/rand-to-best/1/bin的方程分别见公式（5）、（6）。

其中，νi,G+1表示生成的变量向量，i表示种群个体序列号，G表示种群进化代数，r1≠r2≠r3≠i，且i、r1、r2、r3∈[1,NP]，NP表示种群规模，xbest,G为第G代种群中的最优个体，F、λ是缩放因子，为常数。

2.2 改进DE算法 本文在基本DE算法的基础上，引入随机最佳变异和局部增强算子，提出一种能在保证种群多样性前提下加速收敛、提高局部搜索能力的改进DE算法。具体实现如下：

2.2.1 随机最佳变异 在基本 DE算法中，为克服DE/rand/*/*变异策略收敛速度较慢的缺点，普遍采用DE/best/*/*变异策略，在变异过程中，种群中个体均采用相同的最佳个体，对加快算法收敛性有良好的引导作用，但最佳个体周围集中大量的寻优个体，种群多样性遭到破坏，容易导致算法局部收敛。

为使变异向量具有较好的多样性和引导性，本文采用一种随机最佳变异策略，见公式（7）。

其中，r1≠r2≠r3≠i，且i、r1、r2、b∈[1,NP]，RMS(xb,G)≤min{RMS(xr1,G), RMS(xr2,G)}。

随机最佳变异策略将3个随机选择的个体中的局部最优个体当作变异算子的基向量，而另外2个向量作为差向量对，这样既保证种群多样性，又为种群优化提供了一定的引导作用，避免种群大量聚集于最佳个体周围，从而避免了算法的“早熟”。

2.2.2 局部增强算子 引入随机最佳变异策略，保证了算法收敛速度和种群多样性之间的相对平衡，但算法收敛速度相对较慢，为进一步加快算法收敛，引入局部增强算子，对除最佳个体外的其他种群个体以比例因数CP（0＜CP＜1）的概率进行重新赋值，使部分个体向最佳个体靠拢，以增强这些个体的贪婪性，加快算法收敛速度，局部增强算子见公式（8）。

其中，r1≠r2≠i，且i、r1、r2∈[1,NP]，ϲi,G+1是增强后的新个体，用来代替个体xi,G+1，xbest,G+1是新种群中的最佳个体，g为算法迭代次数，改进算法引入新的比例因数CP，当CP=1时，除最佳个体外，其余所有个体向量都需重新赋值；当CP=0时，未选中个体进行局部增强。

局部增强算子的目的是在增加部分个体贪婪性的同时保证种群多样性，以使算法达到既快又好地逼近到全局最优解。算法在增加局部增强算子后，每次迭代仅按一定比例使部分个体变得更加贪婪，从而在一定程度上限制了算法的贪婪性，避免影响算法的收敛性。需注意的是，随着迭代过程的进行，种群中这部分个体随机扰动范围也在动态缩小，这种方式对算法局部搜索能力有较大的提高，尤其是能减少收敛到全局最优解所需的迭代次数，但算法也有陷入局部最优的可能。

3 仿真实验与结果

3.1 创建仿真模型 在ANSYS 10.0软件上采用有限元法对头部进行建模，并对EIT正问题进行求解。三层有限元头球模型如图 2A，由 14 332个单元和20 193个节点组成。头部模型分为头皮、颅骨、大脑及病变组织4个部分[12]。为简化计算，假定头部模型为各向同性均质导体，且电阻率值已知。根据先验知识，将模型 4个部分的电阻率按比例设定为：头皮∶颅骨∶大脑∶病变组织=1∶15∶1∶2[13]，其目标电阻率分布的剖面图如图2B，其中浅色区域为病变组织。

图 2 采用有限元法对头部进行建模。A.头球模型；B.目标电阻率分布

3.2 仿真参数设定 本文对三维头球模型进行仿真实验，并采用相对电极注入方式向模型左右注入±5 mA的电流。在逆问题重建实验中，将头球模型的头皮、颅骨、大脑及病变组织4个部分的上限设为[1.20，17.40，1.20，2.25]，下限设为[0.81，13.60，0.81，1.75]，初始种群需从此约束边界范围内随机生成。由于大脑和头皮电阻率的设定值相等，为简化计算，假定变量个数D=3，种群规模NP=60，目标函数为公式（3）。算法终止条件为迭代次数达到100次。根据文献将比例因数F、交叉因数CR分别设定为0.8、0.9[13]。局部增强算子的比例因数CP等到下文分析此参数对重建结果的影响后再另行确定。

3.3 仿真实验与结果 为了对算法图像重建质量进行定量分析，分别引入相对误差RE和误差总和TE准则，见公式（9）、（10）。

其中，ρ*是指模型各个部分的重建电阻率值，ρ0是指相应的模型各个部分的目标电阻率值。

其中，m为电阻率分布自由度，ρitar和ρirec分别代表单元i的目标电阻率值和重建电阻率值。

图3为改进DE算法的重建图像误差总和TE随局部增强因子的比例因数CP变化的曲线，当CP=5/60时，重建图像的整体误差最小，算法效果最佳；当CP为1/60～5/60时，在相同的迭代次数下，随着CP取值的增加，重建图像的误差逐渐减小，说明当聚集在最优个体周围的个体在一定范围内时，会在一定程度上加快算法收敛的速度，提高重建图像的精度；当CP为5/60～15/60时，随着CP取值的增加，重建图像的误差反而逐渐增加，说明当聚集在最优个体周围的个体超出一定范围后，算法的贪婪性过大，反而使算法陷入“早熟”。因此建议采用此改进方法进行三维脑部阻抗图像重建时，将CP的取值范围确定在 3/60～7/60，即可能获得更好的图像重建效果。在本文的仿真对比实验中选择局部增强因子CP=4/60。

图3 误差总和随局部增强因子的比例因数变化的关系

采用上述实验参数设置，在相同的初始条件下，将本文提出的改进算法与其他算法进行仿真对比实验。不同算法的电阻抗重建图像结果见图4。

图4 不同算法电阻率重建分布

其中算法DE-rand、DE-best分别为采用文献[13]中的DE/rand/1/bin变异策略和DE/best/1/bin变异策略，算法PSO采用文献[14]中的粒子群优化算法。比较图4及图2B发现，这些算法均能在电阻抗重建过程中得到较为清晰的图像，且能对病变组织进行准确定位，但不同算法的电阻率重建图像与目标电阻率的分布图像仅在颜色深浅上有细微差别，用肉眼无法看出明显差异，因此需要对阻抗重建结果进行定量分析。

表 1是经过 100轮迭代后不同算法的图像重建结果，其中t表示各个算法单次循环所需平均时间，其单位为分钟。由表1可见，改进算法相对于算法DE-rand、PSO、DE-best而言，其模型各个部分整体的重建误差总和TE分别降低了66.996%、32.723%、49.834%；其单次循环所需平均时间分别降低了 8.331%、2.017%、4.793%，表明在相同条件下改进算法图像重建精度更高、迭代效率更高，算法整体性能更好。

表1 不同算法的图像重建实验结果

由表2可见，与其他部分相比，各个算法对病变组织部分的重建误差更大，表明在求解EIT正问题时，大脑颅骨会阻碍一部分电流的流入，导致病变部分的重建误差变大，大体上符合理论情况。从总体上看，改进算法相比算法 DE-rand、PSO、DE-best，其模型的大脑、颅骨、头皮、病变部分电阻率重建误差分别降低86.794%、68.690%、86.720%、86.874%，-16.084%、33.431%、-18.571%、81.164%，80.378%、53.955%、80.238%、70.146%，表明改进算法有效降低了图像重建的误差，提高了图像重建的精度，算法性能得到一定的提升。

表2 头球模型不同部分的图像重建相对误差（%）

由图 5可见，所有算法在迭代初期收敛速度较快，随着迭代次数的增加，收敛速度慢慢趋于平缓。当迭代次数＜18时，改进DE算法较其他算法无明显优势，说明改进算法在种群多样性和算法收敛速度上保持了较好的水平；当迭代次数为 18～27时，改进算法的局部增强算子逐渐发挥作用，使改进算法的重建效果逐渐优于算法1和算法2；当算法迭代次数为18～100时，改进算法明显优于其他算法的效果和精度。

图5 种群平均目标函数值随迭代次数变化关系

4 结束语

电阻抗图像重建算法是 EIT研究的关键技术之一。本文提出了一种新的用于电阻抗重建的改进 DE算法。该算法通过引入随机最佳变异策略和局部增强算子，在算法收敛速度和种群多样性之间取得良好的平衡，并在一定程度上加快了算法收敛速度，增强了算法局部搜索能力。结果表明，改进算法能够获得较为清晰的电阻抗重建图像，且能对病变组织进行准确定位。同时，与其他算法相比，改进算法的图像重建误差减小30%以上，有效提高了算法的整体性能，为EIT临床研究提供了一定的参考。

由于DE算法是一种相对较新的全局类智能寻优算法，仍有较多的相关理论问题亟需完善。改进算法由于每次均需对所有种群个体进行目标函数评价，仍存在测量数据、计算时间较长的缺点。因此，如何在一定的实验条件下增加测量数据并减少仿真计算时间是下一步研究的重点。

[1] Silvera-Tawil D, Rye D, Soleimani MA. Electrical impedance tomography for artificial sensitive robotic skin: a review. IEEE Sens J, 2015, 15(4): 2001-2016.

[2] 陈会, 闫丹丹. 基于磁共振技术的生物电阻抗成像研究进展. 中国医学影像学杂志, 2014, 22(9): 714-717.

[3] Wagenaar J, Adler A. Electrical impedance tomography in 3D using two electrode planes: characterization and evaluation. Physiol Meas, 2016, 37(6): 922-937.

[4] Schullcke B, Gong B, Moeller K. Steps towards 3D electrical impedance tomography. Conf Proc IEEE Eng Med Biol Soc, 2015: 5323-5326.

[5] 苌飞霸, 张和华, 颜乐先, 等. 电阻抗断层成像技术研究.中国医疗器械杂志, 2016, 40(1): 52-54.

[6] Trepte CJ, Phillips CR, Solà J, et al. Electrical impedance tomography (EIT) for quantification of pulmonary edema in acute lung injury. Crit Care, 2016, 20(1): 18.

[7] 张辉, 李颖, 王西明, 等. 电阻抗断层成像的MPSOMNR算法研究. 计算机工程与应用, 2013, 49(9): 29-32.

[8] Ribeiro RR, Feitosa AS, Souza RD, et al. A modified differential evolution algorithm for the Reconstruction of electrical impedance tomography images//Biosignals and Biorobotics Conference. Salvador: IEEE, 2014: 1-6.

[9] Rashid A, Kim S, Liu D, et al. A dynamic oppositional biogeography-based optimization approach for time-varying electrical impedance tomography. Physiol Meas, 2016,37(6): 820-842.

[10] Martin S, Choi CT. Nonlinear electrical impedance tomography reconstruction using artificial neural networks and particle swarm optimization. IEEE Trans Magn, 2016,52(3): 1-4.

[11] Zhang XL, Wang W, Sze G, et al. An image reconstruction algorithm for 3-D electrical impedance mammography.IEEE Trans Med Imaging, 2014, 33(12): 2223-2241.

[12] Hamilton SJ, Lassas M, Siltanen S. A direct reconstruction method for anisotropic electrical impedance tomography.Inverse Probl, 2014, 30(7): 1375-1378.

[13] 闫丹丹, 陈会. 基于电阻抗成像技术的脑病变组织检测仿真研究. 中国医学影像技术, 2014, 30(7): 1113-1116.

[14] 陈民铀, 杨艳利, 何为, 等. 基于粒子群优化算法的电阻抗图像重建. 重庆大学学报, 2011, 34(1): 82-87.