基于混合全变差正则化算法的电阻抗成像技术研究

左从磊， 李 静

(中国科学技术大学 工程科学学院，安徽 合肥 230027)

0 引 言

电阻抗成像(electrical impedance tomography,EIT)技术作为一种新型的成像技术，具有无辐射、非侵入性、成本低等优势，在医学成像、地质勘察、流体成像、材料检测等领域都有广阔的应用前景[1,2]。EIT的目标是重构检测对象的电导率，一般在边界放置若干个测量电极，并向成像物体中注入一定幅值的安全激励电流，通过测量并处理边界电压来重构电导率分布。由于从边界电压获得的数据数目远小于未知的电导率分布，EIT重构结果具有不唯一性和不适定性；另一方面，测量容易受到噪声的影响，边界电压的微小扰动会导致重构过程发生很大变化，使得EIT成为一个严重病态的逆问题。常用的L2范数正则化算法使逆问题的病态性得到改善，但由于L2范数的处处连续性质，使重构图像的边界不理想，同时重构图像往往具有比较严重的伪迹和重影，这些都制约了L2范数正则化的运用[3～5]。近年来，L1范数的使用被越来越多的研究者所考虑[6,7]。

本文基于原始对偶内点方法[8]，通过组合L1和L2范数，实现混合全变差(hybrid total variation,HTV)正则化重构算法，有效提高了图像重构质量，并在此基础上进一步研究了算法的抗噪声性能。

1 相关理论

在各项研究中，共轭梯度(conjugate gradient，CG)法[9]、范数正则化等方法被广泛应用于EIT逆问题求解。其中,常用的L2范数正则化方法通过添加罚函数项对解起到一定的阻尼作用，使逆问题解保持稳定，可以表示为

式中lj为相邻边长度，gm(j)和gn(j)分别为相邻单元元素值。

将L1和L2正则化组成混合罚函数项，以发挥L1和L2范数的优势。定义混合HTV正则化算法的目标函数为

α2‖L2(g-g0)‖2}

(3)

式中Lj为变差矩阵，为式(2)的矩阵形式，L2为NOSER正则化矩阵，α1和α2分别为L1范数和L2范数的正则化参数。

2 算法推导

式(3)是一个凸优化问题，由于L1范数的不可微分性质，求解需要引入平滑参数。而传统的算法对于平滑参数的选择十分苛刻，其值过大或者过小对图像的质量都有很大影响[12]。Anderson K D等人[8]已经证明原始对偶内点法(primal-dual interior-point method,PDIPM)在求解L1范数上具有巨大优势。PDIPM算法通过引入对偶变量将不等式约束转化为等式约束[7]，对平滑参数的选择具有高的鲁棒性。本文基于PDIPM推导了HTV算法的求解框架如下。把式(3)标记为原始(P)问题

α2‖L2(g-g0)‖2}

(4)

根据PDIPM理论，对偶(D)问题可表示为

α2‖L2(g-g0)‖2]},‖x‖≤1

(5)

式中x为引入的对偶变量，L1为二维EIT有限元问题的L1范数矩阵。为得到内部最小化条件，式(5)对g求一阶偏导数并令其值为0

原始问题是一个最小化问题，而对偶问题是一个最大化问题。通过消除其互补间隙GPD来求取可行域内的最优点

当且仅当式(8)成立，互补间隙GPD为0

Ljg-xi|Ljg|=0,∀i,j

(8)

如果相邻单元的元素值相同，Ljg为0，上述条件不可微分。引入平滑参数β，克服不可微的问题

使用Newton方法进行计算。式(6)和式(9)分别对原始变量g和对偶变量x求一阶偏导数。为方便描述，式(6)记为Eq6，式(9)记为Eq9

最终得到HTV算法的PDIPM求解框架如式(11)

基于式(11)可以对原始变量和对偶变量进行迭代更新求解。g的更新步长使用优化的线搜索方法，其中,k为迭代次数，λk为线搜索得到的步长因子

gk+1=gk+λkδgk

(12)

为保证x的更新方向是对偶问题的下降方向，使用步长更新规则[8]

xk+1=xk+min(1,φ*)δxk

式中φ*为第k次迭代时xi与可行域边界+1和-1之间的最小距离，N为x的维度。

3 性能参数

如式(3)所示，g的更新取决于三个部分：基于二次函数的残差项、L1罚函数项以及L2罚函数项。对L1和L2罚项的系数做一定调整，可以改变罚项对逆问题解的惩罚作用。通过引入权重系数对罚项系数做如下的变形操作

α(η|L1g|+(1-η)‖L2(g-g0)‖2) }

(14)

算法的思路是平衡L1与L2罚项的阻尼作用，结合L2罚项的光滑和L1罚项的变分优势，重构出分辨率高的图像。算法在成像初期，L2罚项占据重构的优势地位，而在成像后期L1罚项占据重构的优势地位。在成像初期，考虑使逆问题的解快速收敛到真值的附近范围;在成像后期，则考虑利用L1的变差特性增加边界锐利度。一个重要的问题是正则化参数α的选择，α过大或者过小都会对图像的重构结果产生影响。在实际成像中，α值具有较高的鲁棒性[12]。基于先验知识，α值在一定范围内的成像结果都有好的表现。本文HTV算法中α的较好范围是10-4～10-6，仿真算例的重构是在这个基础上进行研究的。对于其他的算例，α的取值范围会有所不同。权重参数η设定为[0,1]之间的等间距数，本文设置了20等分。选用三个指标以定量评判图像重构的质量，一是图像重构误差IRE，二是图像相关系数ICC，三是结构相似度(structural similarity,SSIM)。

1)图像重构误差(image reconstruction error,IRE)

IRE=‖G-‖/‖G‖

(15)

式中G为设定电导率，为重构电导率。IRE越小，表明重构误差越小。

2)图像相关系数(image correlation coefficient,ICC)

3)SSIM

式中uM和uN分别为模型M和重构结果N的均值，ρM和ρN分别为M和N的标准差，ρMN为协方差，表征N相对于M的非线性变化。SSIM(M,N)∈[0,1]，越接近1表明重构质量越好，其对于边缘跳变有较好的区分能力。

4 仿真与讨论

4.1 仿真平台搭建

基于MATLAB平台开发二维圆域仿真算例，将背景区域电导率设定为1 S/m，目标区域设定为2 S/m。为测试算法对不同目标的重构情况，建立了单目标模型、双目标模型、多目标模型，不同模型的目标形状和位置分布也不尽相同。为保证测量条件相同，均采用16电极模型，激励电流幅值为1 mA，激励模式和测量模式均为相邻模式。通过二维仿真数值模拟，测试HTV算法在电导率分布重构中的性能表现，并与TK法、CG法、NOSER法进行了比较。

4.2 无噪声情况下重构成像

图1对比了无噪声时HTV算法以及几种经典的重构算法的性能表现。HTV算法在单目标、双目标、多目标模型中都展现出更好的重构效果。1)HTV算法的成像边界最清晰；Tikhonov法、CG法和NOSER法的成像边界模糊，证明L1范数在重构空间阶跃性变化时具有独特的优势。2)HTV算法的形状重构最准确;三种模型的目标区域及形状被都被较好地定位和重构，相比之下其余方法形变较为严重。3)HTV算法的成像伪迹最少;TK法和CG法的锯齿较严重，NOSER法的重构目标范围过大且伪迹多。

图1 无噪声时不同算法的重构图像

表1给出了不同算法重构性能的定量对比。从表1和图2结果来看，HTV算法的重构结果在三个参数上都表现最好，其次是TK算法。HTV算法三种模型的重构图像IRE值分别是TK法的约78 %，58 %，95 %；ICC值分别是TK法的约105 %，121 %，107 %；SSIM值分别是TK法的约110 %，128 %，113 %。而CG法和NOSER法重构在ICC上表现尚好，但IRE值均在0.9以上，显著高于HTV和TK法，SSIM值相比HTV法也非常小，究其原因是这两种方法虽然可以一定程度上重构分布范围，但重构电导率与真实值相差大。从总体上看，随着目标区域的增大和形状的复杂程度增加，重构图像相关系数都出现了一定的下降。

表1 无噪声时不同算法的重构性能参数

图2 无噪声时不同算法的重构性能参数柱状图

4.3 含噪声情况下重构成像

在实际应用中，噪声是不可避免的。为了测试算法的抗噪声性能，在正问题生成的测量电压中加入了SNR为40 dB的随机高斯噪声。

图3和表2给出了当SNR为40 dB时算法的重构结果对比。NOSER算法的抗噪性能最好，与理想无噪声情况下相比重构图像变化最小，三种模型的重构结果中ICC高于其余算法。而HTV算法因为结合了NOSER正则化罚项，抗噪性能较好。相较于无噪声情况，当目标靠近边界时出现了一定的形变，但是HTV算法的IRE值在几种算法中最小，抑制了伪迹且有清晰的边界。在医学成像领域，不同医学组织之间的边界是具有价值的图像信息[6]。从这一点上而言，HTV算法以抗噪性能略微下降的代价获得更清晰的图像是值得的。Tikhonov法的重构结果表现不佳，模型2和模型3的重构目标区域已经很难分辨，图像的伪迹和锯齿比较严重。CG法表现不佳，其重构电导率值与真实值差异大，相关系数和结构相似度均较小。

图3 含噪声时不同算法的重构图像

表2 含噪声时不同算法的重构性能参数

4.4 真实肺部图像重构

人体在进行肺部呼吸时，电导率会发生一定的变化。Andler A等人[13]采集了人体呼吸时肺部EIT数据。16个测量电极被等间距地放置在一个成年人的胸部周围，并以相邻模式进行电流激励和电压测量。幅值为10 mA的激励电流通过相邻的电极对注入，并在其余相邻的电极之间测量出电压差。采用Andler A等人采集的数据，应用HTV算法和其余算法对肺部图像进行重构。

图4是不同算法的重构图像结果对比，结果表明HTV算法相较于其他方法成像质量最高，重构图像边界清晰、伪迹较少，形状最接近真实肺部。因此HTV算法提高了肺部重构图像的重构质量，在实际应用中有显著优势和良好的应用前景。

图4 真实肺部图像重构图像结果对比

5 结 论

本文创新使用L2范数的NOSER罚函数项和L1罚函数项的HTV正则化重构算法，并基于原始对偶内点法推导出算法的求解框架。通过建立不同大小和形状的仿真测试模型，在EIT正问题解的基础上，使用HTV算法实现了图像重建，并将其与经典Tikhonov、共轭梯度、牛顿一步误差重构算法进行了比较。结果表明：HTV算法在理想无噪声情况下的重构表现优于其他算法，重构图像的边界更为清晰且形状最准确。图像的重构误差IRE最小，ICC和SSIM最高。在含噪声情况下(SNR为40 dB)时HTV算法也有良好的抗噪性能，可以重构出高质量的图像。基于人体肺部呼吸EIT采集数据进行重构，证明HTV算法重构成像清晰，有较好的实用价值和应用前景。