基于非精确Newton迭代正则化的EIT图像重构算法

赵一帆,王 静

(黑龙江大学 数学科学学院,哈尔滨 150080)

0 引 言

电阻抗成像(Electrical impedance tomography,EIT)是一种新兴的无损功能成像技术,根据目标体表面的边界测量数据来估计目标体内部的电导率分布情况,具有非侵入性、设备轻便、成本低廉以及无损检测等特点,在医疗诊断、地理探测和工业检测等领域有着广泛的应用前景。

图像重构是EIT的核心技术,直接影响着成像的空间分辨率和实时性。EIT图像重构问题本身存在严重的不适定性,即观测数据的微小变化会导致重构参数发生较大的变化,而有限的测量数据以及噪声的干扰都加剧了这种不适定性。为克服重构问题的不适定性,往往需要借助不同的正则化技术来提高重构图像的质量。国际上关于EIT图像重构方法的研究非常多,大致可分为3类：Tikhonov正则化为代表的直接方法、Landweber迭代为代表的迭代正则化方法,以及神经网络为代表的智能方法。Tikhonov正则化[1-2]是克服EIT重构问题不适定性和病态性的有效方法,但正则化参数的选取较为复杂,直接影响成像效果,实际一般采用经验值,具有一定的局限性。神经网络[3-5]是智能方法中研究较多的,也是当前研究的热点,可直接建立边界测量电压数据与电导率分布之间的非线性映射关系,但其需要大量的样本集来训练网络模型,对训练样本的完整性要求较高,且对未知模型无法预测,在一定程度上也限制了其应用。

迭代正则化方法中,最为经典的Landweber迭代(Landweber iteration,LDI)[6-8]因其在实现图像重构的过程中只需用到数据拟合的梯度信息,易于实现且计算成本低,稳定性好,也是目前图像重构中应用最为成功的一种迭代方法。然而,灵敏度矩阵的病态性使得LDI迭代存在收敛速度较慢的不足,需要大量的迭代步来实现高质量的成像,影响了其应用,因此需要寻找更加快速的方法。近来,备受关注的非精确Newton迭代正则化方法[9-13]具有较快的收敛速度,是一种先线性化后正则化的方法。但需注意的是,线性化过程并不能消除问题本身的不适定性,仍需借助正则化策略处理局部线性化方程。该方法包括内外两层迭代：内层迭代应用迭代正则化方法求解局部线性化问题产生近似迭代序列,外层迭代为非精确Newton法用于更新迭代点列。稳定性较好且格式简单的LDI迭代法[14]常常被用作为内层迭代正则化策略,但鉴于LDI迭代法收敛速度较慢的缺陷,若将其作为内层迭代策略,外层每一步迭代往往需要大量的内层迭代次数来满足停止准则,造成内层迭代计算量较大。2009年,韩波等首次提出将同伦摄动迭代(Homotopy perturbation iteration,简称HPI)用于求解非线性不适定反问题[15]。随后,HPI迭代法被应用于多个领域,如电阻抗成像[16]、地震勘探测井约束反演[17]等,数值结果表明：当达到相似重构精度时,HPI迭代法仅需LDI迭代法大约一半的计算时间。近来,HPI迭代法作为一类加速策略得到了推广应用[18-19]。基于此,本文考虑采用HPI迭代作为内层迭代并结合适当的步长选取策略来加速内层迭代速度,从而获得一种更高效的非精确Newton迭代正则化方法。

1 EIT数学模型

目前广泛采用的EIT数学模型为比较贴近实际的全电极数学模型[20]：

(1)

EIT问题的研究主要包括正问题和反问题。EIT正问题可归结为求解边值问题(1),即给定内部电导率分布、注入电流和接触阻抗,计算内部电位u和电极电位U。正问题的求解可借助于有限元方法[21]对边值问题(1)进行离散化,可得到如下离散形式：

A(σ)b=f

(2)

其中,A(σ)为有限元离散矩阵,b为待定系数 (用于计算u和U),f为包含电流激励源的向量。

然而,在EIT实际应用中,电导率分布σ是未知的,仅已知的是电极处的注入电流数据和测量电压数据,而测量电压数据在数据采集过程中又不可避免地含有噪声,这里记为Uδ。EIT观测模型响应可表示为如下的非线性算子方程

Uδ=F(σ)+e

(3)

这里,F表示模型空间到测量数据空间的投影,即求解正问题的过程,也称之为正演算子,e表示测量误差向量。因此,根据边界电极处测量电压数据Uδ反推内部电导率分布σ的过程就称之为EIT反问题,由于需要通过图像的形式呈现出来,故也称之为图像重构。

2 非精确Newton迭代正则化方法

2.1 LDI迭代法和HPI迭代法

EIT图像重构问题本身存在严重的不适定性,最小二乘法是处理该问题最为流行的方法,即极小化数据拟合项

(4)

最为经典的LDI迭代格式为

σn+1=σn+μnF′(σn)*(Uδ-F(σn))

(5)

可视为(4)的最速下降法;HPI迭代法则是通过由同伦方法求解(4)的欧拉方程F′(σ)*(F(σ)-Uδ)=0得到的,其迭代格式为

σn+1=σn+μnF′(σn)*(2I-F′(σn)F′(σn)*)(Uδ-F(σn))

(6)

式中μn为适当选取的步长因子,这里采用比较适宜在EIT领域使用的变步长因子选取准则,满足

(7)

2.2 INLDI迭代法和INHPI迭代法

鉴于非精确Newton迭代正则化方法具有较快的收敛速度,考虑将其用于处理EIT图像重构问题。该方法包括内外两层迭代：外层迭代为非精确Newton方法且以偏差原则作为停止准则,内层迭代这里考虑采用LDI迭代或HPI迭代并结合选取适当的步长来加速。

首先构造内层迭代。假设已知当前迭代步为σn,通过极小化问题

(8)

构造相应的迭代格式

σn,k+1=σn,k+μn,kdn,k

(9)

(10)

若选取

(11)

则为LDI迭代格式;若选取

(12)

‖sn,k‖≤γ‖Uδ-F(σn)‖

(13)

‖Uδ-F(σn*)‖≤τδ≤‖Uδ-F(σn)‖,0≤n≤n*

(14)

作为停止准则,其中n*=n*(δ,Uδ)为由(14)所确定的停止指标。

算法1 EIT图像重构的非精确Newton迭代正则化方法Algorithm 1 Inexact Newton iterative regularization method for EIT image reconstruction

最后,以LDI迭代作为内层迭代策略,引入基于LDI迭代的非精确Newton方法(Inexact Newton-LDI,INLDI),其格式为

(15)

以HPI迭代作为内层迭代策略,引入基于HPI迭代的非精确Newton方法(Inexact Newton-HPI,INHPI),其格式为

(16)

以上求解EIT图像重构问题的非精确Newton迭代正则化方法的具体过程如算法1所示。

3 数值模拟

仿真模拟基于Matlab环境下Vauhkonen开发的EIDORS 2D软件[22],建立16个电极的二维圆形EIT电极模型,采用传统的相邻电极电流激励和相邻电极电压测量模式。设目标区域Ω是坐标系上规则的一个圆域,半径为14,半径无需标注单位,仅作背景和目标大小对比。为了避免反问题中的“Inverse Crime”,考虑了两种不同规模的有限元网格剖分用于正反问题的求解,即用于正问题的网格由1 968个三角单元和1 049个结点构成,如图1(左一)所示,用于反问题的网格由492个三角单元和279个结点构成,如图1(左二)所示,其中网格周围的绿色条纹表示理想电极的位置。重构过程中,电极处的接触阻抗假设已知,设置接触阻抗大小为0.05。

图1 有限元网格：细网格(左一)和粗网格(左二);三种测试模型：(a) 模型1;(b) 模型2;(c) 模型3Fig.1 FEM meshes：fine mesh (left one) and coarse mesh (left two),three test models：(a) model 1;(b) model 2;(3) model 3

为验证本文提出的INLDI迭代法和 INHPI迭代法在非线性EIT图像重构问题中的重构性能,构造了含有不同异物的三种测试模型作为研究对象,如图1(a)～(c)所示,模型1包含1个圆形异物,模型2 包含2个方形异物,模型3包含5个不同异物,其中的标尺显示了各区域电导率σ的倒数(电阻率),蓝色异物电阻率取值为1Ω·m,红色背景电阻率取值为6Ω·m。为了便于对比分析,选取传统的LDI迭代法、HPI迭代法与本文算法进行比较。

真实测量数据由合成数据Usyn=F(σ*)(σ*为真实的电导率分布)计算得到,噪声测量数据通过在合成数据上添加随机噪声产生,即Uδ=Usyn+δ·n,其中δ为噪声水平,随机变量n为服从[0,1]上的Gaussian分布。此外,为了比较EIT图像重构质量,选择较常用的相对误差(RE)及相关系数(CC)作为衡量算法的客观评价指标,其计算公式如下：

数值模拟的相关参数选取为τ=2.01,γ=0.8,初值σ0选取为背景值。由于过长的计算时间并不符合EIT图像重构的需求,因此设置迭代的最大步数为nmax=10 000,kmax=50 000,即意味着即使未达到停止准则,当迭代进行到最大迭代步时停止迭代。

针对三种测试模型,为了测试四种方法的重构效果,图2描述了噪声数据下相对误差随迭代时间的变化曲线,图3描述了噪声数据下相关系数随迭代时间的变化曲线,用以分析算法的稳定性和有效性。观察图2和图3可以看到：四种方法可达到同等效果的RE值和CC值,说明具有相似的重构精度和重构质量,但相比LDI和HPI迭代法,INLDI和INHPI这两种迭代法较早地达到了稳定状态,即RE曲线下降的更快、CC曲线上升的更快,这正体现了Newton迭代法具有收敛特性上的优势,具有更快的收敛速度,而且也注意到,相比INLDI迭代法,INHPI迭代法具有更快的收敛速度,随着噪声水平的减少,INHPI迭代法在计算效率上的优势更加明显,主要是内层HPI迭代步引起的(具体见表1中的内迭代步比较)。另一方面,从图2和图3中也可注意到：随着噪声水平的减少,四种方法获得的RE值在变小而CC值在增大,说明重构图像的重构精度和重构质量变得更好,但需要的计算时间却大幅度增加了,这说明收敛速度与噪声水平有关。

表1 LDI、HPI、INLDI以及INHPI迭代法重构模型2的数值比较Table 1 Comparison by LDI,HPI,INLDI and INHPI for model 2

图2 相对误差曲线：(a) 模型1;(b)模型2;(c)模型3Fig.2 Relative error curves：(a) model 1;(b) model 2;(3) model 3

图3 相关系数变化曲线：(a) 模型1;(b)模型;(c)模型3Fig.3 Correlation coefficient curves：(a) model 1;(b) model 2;(3) model 3

针对图2和图3中的结果,这里仅以模型2为例(见图2(b)和图3(b)),表1记录了在不同噪声水平下,四种方法的详细数值结果比较,其中n*表示LDI迭代法与HPI迭代法的迭代终止步,nouter表示INHPI方法与INLDI方法的外层Newton迭代步,ninner表示INHPI方法与INLDI方法的内层总迭代步,Time(s)表示计算所需时间。从表1可以看出,四种方法具有相似的重构精度,但重构效率明显不同。相比LDI和HPI迭代法,INHPI和INLDI迭代法的重构时间显著减少,大大提高了重构效率。另一方面,比较LDI和HPI迭代法,验证了HPI迭代法只需LDI大约一半的迭代步和计算时间;但比较INLDI和INHPI迭代法,INHPI迭代法总的内层迭代次数比INLDI迭代法减少了约一半,高噪声水平下的重构效率虽有所提高,但并未像预期的那样提高约一倍,主要原因在于需要调用正问题的外层迭代为收敛速度较快的Newton步(外层每步迭代耗时多,但需要的外层迭代步数少),内层每步迭代耗时较少,迭代步数虽节约了一半左右,但在内层迭代迭代步数基数较小的情况下,计算效率无法凸显出来,但在低噪声水平下,随着需要的内层迭代步数的剧增,相比INLDI迭代法,INHPI迭代法在计算效率上的优势就呈现出来了,改进了成像效率。

最后,为了可视化重构结果,鉴于四种方法具有相似的重构精度和重构质量,针对三种测试模型,这里仅给出不同噪声水平下INHPI迭代法的重构图像,如图4所示。从图4可以看出,在低噪声水平下,重构图像均能够较准确地反映目标体的大小、位置和性质,随着噪声水平的增加,重构结果比较稳定,说明对噪声较为鲁棒,不过随着噪声水平的增加,重构图像的分辨率有所降低。

图4 基于不同噪声数据的INHPI迭代法所得重构图像,从左到右依次为：真实模型,以及δ=5%,δ=2%,δ=0.5%,δ=0.05%噪声数据下对应的重构图像Fig.4 Reconstructed images for three test models by INHPI method under different noise levels.Reconstructions from left to right：true model,reconstructed images under δ=5%,δ=2%,δ=0.5%,δ=0.05%,respectively

4 结 论

基于LDI迭代法和HPI迭代法,提出了两种非精确Newton迭代算法：INLDI迭代法和 INHPI迭代法,用于处理非线性EIT图像重构问题。其基本思想是应用LDI迭代法或HPI迭代法逐步求解局部线性化的EIT观测响应模型来构造内层迭代序列,并结合适当的停止准则终止内层迭代,进而将得到的近似重构解作为新的外层迭代点并以偏差原则作为外层停止准则。对LDI迭代法和HPI迭代法进行了比较,并以重构图像的相对误差和相关系数作为评价标准。结果表明,本文算法具有节约迭代步以及提高计算效率方面的优越性,具有更快的收敛速度,可有效改进成像效率。