CN109118553A - 基于几何形状约束的电阻抗层析成像内含物边界重建方法 - Google Patents

Abstract

本发明涉及一种基于几何形状约束的电阻抗层析成像内含物边界重建方法，包括：1)采用局部弧长参数x(s)表征目标内含物的边界；2)基于残差函数和几何形状约束构建形状反演问题的能量函数：3)根据变分法原理，使能量函数ε(x)最小化的最优边界估计x满足拉格朗日方程；对能量函数中的变量进行离散化处理，由一系列离散点[x(s1),x(s2),…,x(sN)]表征内含物边界；4)采用半隐式方法迭代求解拉格朗日方程；5)经过多次迭代后边界估计值上的采样点可逐步逼近目标内含物的真实边界，实现内含物边界重建。

Description

基于几何形状约束的电阻抗层析成像内含物边界重建方法

技术领域

本发明属于电阻抗层析成像技术领域，涉及一种基于几何形状约束的内含物边界重建方法。

背景技术

电阻抗层析成像(Electrical Impedance Tomography，简称EIT)是一种具有非侵入或非扰动特点的过程可视化在线监测技术。它通过安置在待测敏感场的阵列式传感器，向目标物场施加电学激励信号，并能够获得反映敏感场内电导率分布信息的电学响应信号，进而实现场内介质分布的二维/三维可视化。该技术具有便携、低成本和高时间分辨率等优点，在工业和生物医学方面具有广阔的应用价值。然而，EIT的图像重建问题具有严重的非线性和病态性，这导致EIT的空间分辨率较低，且易受噪声干扰影响。对于非线性的解决方法通常是采用迭代线性化或直接非线性的方法。而病态性则可以通过对先验信息的利用得到改善,常用的约束方式有平滑约束，稀疏约束等。从其他测量模式，如超声测量，计算流体力学，动态数据序列等，获得的先验信息也有助于提升EIT的空间分辨率。

EIT可以应用于如乳腺肿瘤诊断、长期呼吸监测以及流体中的气泡测量等领域，在这些应用中观测域内的电导率分布近似为分段常数，因此EIT的目标是重建有限个嵌入于均匀背景电导率中的单连通子区域，这些单连通子区域就是所需重建的内含物。这种分段常数约束具有保留内含物形状以及增强内含物边界分辨率的优势，因此吸引了EIT领域的诸多关注。

2007年，S.Babaeizadeh等人发表于《IEEE Transactions OnMedical Imaging(IEEE医学影像处理)》第26卷，第637-647页，题为《Electrical impedance tomographyfor piecewise constant domains using boundary element shape-based inversesolutions(基于边界元形状的EIT分段常数域逆问题解决方法)》，提出了一种利用分段常数约束的内含物边界重建法，该方法通过一组形状系数参数化目标内含物的边界，之后通过最小化形状能量函数计算得到最优的内含物边界估计。由于只需重建内含物的边界，该方法也被称为基于形状的重建方法，它相较于基于像素的重建方法减少了一维自由度，因此内含物边界重建法在改善EIT的病态性，提升EIT的空间分辨率方面有广阔的前景。典型的边界重建方法有2014年H.Haddar等人发表于《Complex Variables and EllipticEquations(复变函数与椭圆方程)》第59卷，第863-882页，题为《A conformal mappingmethod in inverse obstacle scattering(障碍反散射中的保形变换方法)》提出的保形变换法，2012年F.Cakoni等人发表于《Inverse Problems(逆问题)》第29卷，第015005-015027页，题为《Integral equation methods for the inverse obstacle problem withgeneralized impedance boundary condition(广义阻抗边界条件下的逆问题的积分方程方法)》提出的积分法，以及基于边界元法的一些方法。

在内含物边界重建方法的设计中，关键因素是对边界参数化模型的改进。通用的边界参数化模型分为两种，全局边界模型以及局部边界模型。全局边界模型，如：2007年S.Babaeizadeh等人发表于《IEEE Transactions On Medical Imaging(IEEE医学影像处理)》第26卷，第637-647页，题为《Electrical impedance tomography for piecewiseconstant domains using boundary element shape-based inverse solutions(基于边界元形状的EIT分段常数域逆问题解决方法)》采用的球面谐波模型，以及1999年D.K.Han等人发表于《Journal ofComputational Physics(计算物理学杂志)》第155卷，第75-95页，题为《A shape decomposition technique in electrical impedance tomography(电阻抗层析成像中的形状分解技术)》采用的傅里叶模型，可以在重建过程中体现目标边界的全局几何特性，如平滑性等。而一些局部几何特性如曲率等则通常定义于局部边界模型，比如，2007年S.Babaeizadeh等人发表于《IEEE Transactions On Medical Imaging(IEEE医学影像处理)》第26卷，第637-647页，题为《Electrical impedance tomography for piecewiseconstant domains using boundary element shape-based inverse solutions(基于边界元形状的EIT分段常数域逆问题解决方法)》采用的B-样条模型，以及2017年M.M.Zhang等人发表于《IEEE Sensors Journal(IEEE传感器杂志)》第17卷，第8263-8270页，题为《Quantitative reconstruction ofthe exteriorboundary shape of metallicinclusions using electrical capacitance tomography(采用电容层析成像定量重建金属夹杂物的外边界形状)》采用的建立于离散边界点上的模型。两种边界模型具有各自的优势以及局限性。一方面，由于全局边界模型包含了正则化的过程，基于全局边界模型的形状重建方法要比基于局部边界模型的方法更加稳定。另一方面，由于局部边界模型具有更高的局部形变自由度，因此局部边界模型比全局边界模型更加的灵活。对于一些复杂的形状如凹形状等，采用局部边界模型可以更有效地进行表征。

发明内容

本发明针对内含物边界重建方法中，全局和局部参数化边界模型不能同时满足高稳定性和高精确度的边界重建的问题，提出一种基于几何形状约束的内含物边界重建方法，该方法采用几何形状约束对局部边界参数化模型进行正则化约束，并设计了新的能量最小化方程用来约束边界重建问题，以达到改善电阻抗层析成像问题的病态性、同时提高边界重建稳定性和精度的目的。技术方案如下：

一种基于几何形状约束的电阻抗层析成像内含物边界重建方法，包括下列步骤：

1)采用局部弧长参数x(s)表征目标内含物的边界，其中x∈Γ表示目标边界上的点集，s∈[0,1]为局部弧长参数，构建灵敏度矩阵J并获取边界测量值向量U，采用一组已知电导率分布的测量值作为参考电压U_ref，用以对所测边界测量值U进行归一化，归一化后的边界测量值向量为

2)基于残差函数和几何形状约束构建形状反演问题的能量函数：

其中，R(x)为残差项，P_Γ(x)为几何形状约束项，U(x)是由估计的内含物边界计算而出的边界电压的估计值，表示2-范数的平方，符号’和”表示x(s)对s的一阶微分和二阶微分，表示沿目标边界Γ的曲线积分，超参数α和β用来调整几何形状约束的程度，张力项约束边界上两相邻点间的距离变化，用来控制边界的拉伸程度，刚度项约束边界上的扰动，用来控制边界的曲率；

3)根据变分法原理，使能量函数ε(x)最小化的最优边界估计x满足以下拉格朗日方程：

其中，符号””和”表示x(s)对s的四阶微分和二阶微分；M为边界测量值的总数；n为内含物边界的外法线方向向量；Jⁱ表示测量值Uⁱ关于单位边界在点x∈Γ处沿外法线方向位移的灵敏度，其计算公式为：

其中，κ＝σ_k/σ_b是内含物电导率σ_k与背景介质电导率σ_b的比值；φ是由正问题计算出的边界电势；是沿边界外法线方向的微分算子；是沿边界切线方向的微分算子；Iⁱ＝[I₁,I₂,…,I_L]^T是给定的激励电流向量，L为电极个数，Mⁱ是由激励测量策略决定的测量向量；

对能量函数中的变量进行离散化处理，由一系列离散点[x(s₁),x(s₂),…,x(s_N)]表征内含物边界，其中N为边界点的总数，拉格朗日方程的矩阵形式表示如下：

其中，N是N×N维的对角矩阵，对角元素为[n₁n₂,…,n_N]；J是M×N维的灵敏度矩阵；A为循环五对角带状矩阵：

其中，a＝α/δs²，b＝β/δs²，c＝-a-4b，d＝2a+6b；

4)采用半隐式方法迭代求解拉格朗日方程；

5)经过多次迭代后边界估计值上的采样点[x(s₁),x(s₂),…,x(s_N)]可逐步逼近目标内含物的真实边界，之后经过离散点拟合得到重建内含物边界的局部弧长参数x(s)，进而实现内含物边界重建。

附图说明

图1为具体实施例中采用的传感器结构和待重建的目标内含物；

图2为基于几何形状约束的内含物边界重建方法的原理图；

图3为仿真得到的在不同电导率对比度和不同信噪比下的边界重建结果；

图4为实验得到的基于像素方法和形状方法下对单内含物重建的对比结果；

图5为实验得到的基于像素方法和形状方法下对双内含物重建的对比结果。

具体实施方式

以下详细描述本发明所涉及方法的实现步骤，旨在作为本发明的实施例描述，并非是本发明实现的唯一形式，对其他可实现相同结构和功能的实施例也应包括在本发明的范围内。

在具体实施例中，所涉及的EIT系统传感器如图1(a)所示，系统包含16个电极，均匀分布于被测场域外。为简化描述，本实施例中针对单内含物边界重建进行描述。对于多内含物边界重建可由本实施例中的方法直接扩展。EIT系统采用电流激励电压测量且激励电极不测量的模式，采集循环激励循环测量下各个电极上的边界电压构成边界测量值向量U。

该实施例的具体实现流程主要包含以下几个步骤：

(1)参数化表征目标内含物边界

待重建目标内含物如图1(b)中所示，其边界Γ可由一系列局部弧长参数化表征为x(s)，其中x∈Γ表示目标边界上的点集，s∈[0,1]为局部弧长参数。

本发明所提方法是一种绝对成像方法，为了降低该方法对模型误差和噪声干扰的敏感性，本发明提出一种虚拟标定的方法。采用一组已知电导率分布的测量值作为参考电压U_ref，用以对所测边界测量值U进行归一化。归一化后的边界测量值向量可表示为：

边界重建方法的目的就是通过边界电压测量值来估计待重建目标内含物的边界分布x(s)。

(2)基于残差函数和几何形状约束构建形状反演问题的能量函数

根据变分原理，求解边界重建问题等同于最小化残差函数。不同于常用的能量函数构成，本发明结合逼近误差，并额外采用几何形状约束构建了新的能量函数，如下式所示：

ε(x)＝R(x)+P_Γ(x)

其中，R(x)为残差项，P_Γ(x)为几何形状约束项。

残差项测量了估计电压值与实测电压值之间的逼近误差：

其中为边界电压测量值，U(x)是由估计的内含物边界计算而出的边界电压的估计值。表示2-范数的平方。残差项通过最小化估计电压值与实测电压值之间的差值来改进对目标内含物边界分布的估计。但是由于EIT问题的病态性，该过程极易受到测量噪声和模型误差的影响，从而导致不准确的估计结果。

为了改善以上问题，本发明采用几何形状约束作为正则化项，用来对能量函数加以约束。几何形状约束通过张力和刚度来进行阐述：

其中，符号’和”表示x(s)对s的一阶微分和二阶微分，表示沿目标边界Γ的曲线积分，超参数α和β用来调整几何形状约束的程度。张力项约束边界上两相邻点间的距离变化，它用来控制边界的拉伸程度。刚度项约束边界上的扰动，它用来控制边界的曲率。

(3)对能量函数中的变量进行离散化处理

采用本发明所示约束的情况下，最小化能量函数ε(x)就是寻找完全拟合测量值与保持稳定边界之间的最佳折衷。根据变分法原理，使能量函数ε(x)最小化的最优边界估计x满足以下拉格朗日方程：

其中，符号””和”表示x(s)对s的四阶微分和二阶微分；M为边界测量值的总数；n为内含物边界的外法线方向向量；Jⁱ表示测量值Uⁱ关于单位边界在点x∈Γ处沿外法线方向位移的灵敏度，其计算公式为：

其中，κ＝σ_k/σ_b是内含物电导率σ_k与背景介质电导率σ_b的比值；φ是由正问题计算出的边界电势；是沿边界外法线方向的微分算子；是沿边界切线方向的微分算子；Iⁱ＝[I₁,I₂,…,I_L]^T是给定的激励电流向量，L为电极个数，Mⁱ是由激励测量策略决定的测量向量。

为求解以上拉格朗日方程，需对变量进行离散化处理。采用一个恒定的间隔δs对局部弧长参数s进行采样，则待重建内含物边界可由一系列离散点[x(s₁),x(s₂),…,x(s_N)]表示，其中N为边界点的总数。考虑到上述拉格朗日方程中各向量是独立且可分离的，可将各向量分离为x和y部分，用u_j指示x(s_j)或y(s_j)，用n_j指示n(s_j)的x分量或y分量。拉格朗日方程的矩阵形式表示如下：

其中，N是N×N维的对角矩阵，对角元素为[n₁n₂,…,n_N]；J是M×N维的灵敏度矩阵；A为循环五对角带状矩阵：

其中，a＝α/δs²，b＝β/δs²，c＝-a-4b，d＝2a+6b。上述矩阵方程中，第一项是由残差项R导出，第二项是由几何形状约束项P_Γ导出。灵敏度矩阵J表征了EIT测量值与内含物边界变化的局部线性关系。为防止计算过程中产生的截断误差对残差项R和梯度r的影响，在本实施例中对灵敏度矩阵进行了缩放处理，缩放后的灵敏度矩阵等同于J/max(J)。

(4)迭代求解拉格朗日方程

本发明采用半隐式方法迭代求解以上拉格朗日方程。首先，假定方程等式右边的项时不变而等式左边的项随迭代步数改变。其次，假定在每一次迭代过程中r为常量且等同于r^t-1。最后，假定矩阵A在时间t是已知的。

根据上述假设可得如下边界演化方程：

其中，δt为时间步长。由此推导出本发明所提边界重建方法的迭代形式：

其中S＝δtA+E是由几何约束项中计算得出的滤波器，E为单位矩阵。

图2展示了本发明所提方法的计算原理。在每一次迭代过程中，本发明所提边界重建求解器第一步通过梯度下降法最小化残差函数R计算得到最优位移梯度r^t-1，之后第二步立即使用平滑滤波器S调整边界估计x^t-1+δtr^t-1。在迭代过程中，超参数α和β以及采样步长δs都是恒定不变的。时间步长δt的选择需综合考虑稳定性和收敛速度，根据经验，为实现光滑且稳定的边界逼近过程，时间步长的选择需满足δs²＜δtmax(Nr)＜δs。迭代过程在残差足够小或近似为常数时停止，因此，停止条件设定为：

R^t≤τ₁or|R^t-R^t-1|/R^t-1≤τ₂

其中，τ₁是残差的容忍限；τ₂残差变化量的容忍限。容忍限根据试错法选择。

(5)内含物边界轮廓拟合

根据上述方法，经过多次迭代后边界估计值上的采样点[x(s₁),x(s₂),…,x(s_N)]可逐步逼近目标内含物的真实边界，之后经过离散点拟合即可实现具有高稳定性和高精度的内含物边界重建。

实施结果：对上述实施方案进行仿真和实验测试。图3为仿真得到的在不同电导率对比度和不同信噪比下的边界重建结果，结果对比了范围为1.25-1000的电导率对比度以及加入了60dB，40dB和20dB信噪比噪声和无噪声下的边界重建结果，可以看出虽然重建结果随着噪声的增加而变差，但是本发明所提方法仍具有良好的抗噪性，在20dB噪声下，所提方法仍可以较好的重建出目标内含物的位置，尺寸和形状。同时，所提方法的抗噪性可以通过提高电导率对比度来改善。图4为实验得到的基于像素方法和形状方法下对单内含物重建的对比结果，结果对比了NOSER算法成像结果，基于L1范数的稀疏成像结果，TV正则化算法成像结果，以及所提出的基于几何形状约束的边界重建方法(简称GCBR)的成像结果。可以看出，基于像素的方法只能粗略的重建内含物的形状，其中，NOSER算法几乎不能识别内含物的不同形状，而L1和TV正则化只能大致区分椭圆和三角形内含物，不能区分方形和心形的内含物。而所提方法则可以较准确的重建出目标内含物的几何形状。图5为实验得到的对双内含物重建的对比结果，可以看出，相较于基于像素的图像重建方法，所提出的GCBR方法对非圆内含物有更好的重建结果。

Claims (1)

1.一种基于几何形状约束的电阻抗层析成像内含物边界重建方法，包括下列步骤：

1)采用局部弧长参数x(s)表征目标内含物的边界，其中x∈Γ表示目标边界上的点集，s∈[0,1]为局部弧长参数，构建灵敏度矩阵J并获取边界测量值向量U，采用一组已知电导率分布的测量值作为参考电压U_ref，用以对所测边界测量值U进行归一化，归一化后的边界测量值向量为

2)基于残差函数和几何形状约束构建形状反演问题的能量函数：

其中，R(x)为残差项，P_Γ(x)为几何形状约束项，U(x)是由估计的内含物边界计算而出的边界电压的估计值，表示2-范数的平方，符号’和”表示x(s)对s的一阶微分和二阶微分，∮_Γ·ds表示沿目标边界Γ的曲线积分，超参数α和β用来调整几何形状约束的程度，张力项约束边界上两相邻点间的距离变化，用来控制边界的拉伸程度，刚度项约束边界上的扰动，用来控制边界的曲率；

3)根据变分法原理，使能量函数ε(x)最小化的最优边界估计x满足以下拉格朗日方程：

其中，符号””和”表示x(s)对s的四阶微分和二阶微分；M为边界测量值的总数；n为内含物边界的外法线方向向量；Jⁱ表示测量值Uⁱ关于单位边界在点x∈Γ处沿外法线方向位移的灵敏度，其计算公式为：

其中，κ＝σ_k/σ_b是内含物电导率σ_k与背景介质电导率σ_b的比值；φ是由正问题计算出的边界电势；▽_n是沿边界外法线方向的微分算子；▽_t是沿边界切线方向的微分算子；Iⁱ＝[I₁,I₂,…,I_L]^T是给定的激励电流向量，L为电极个数，Mⁱ是由激励测量策略决定的测量向量；

对能量函数中的变量进行离散化处理，由一系列离散点[x(s₁),x(s₂),…,x(s_N)]表征内含物边界，其中N为边界点的总数，拉格朗日方程的矩阵形式表示如下：

其中，N是N×N维的对角矩阵，对角元素为[n₁n₂,…,n_N]；J是M×N维的灵敏度矩阵；A为循环五对角带状矩阵：

其中，a＝α/δs²，b＝β/δs²，c＝-a-4b，d＝2a+6b；

4)采用半隐式方法迭代求解拉格朗日方程；

5)经过多次迭代后边界估计值上的采样点[x(s₁),x(s₂),…,x(s_N)]可逐步逼近目标内含物的真实边界，之后经过离散点拟合得到重建内含物边界的局部弧长参数x(s)，进而实现内含物边界重建。