A Parametric Level Set based Approach to Difference Imaging inA Parametric Level Set based Approach

Liu D, Smyl D, Du J. A parametric level set-based approach to difference imaging in electrical impedance tomography[J]. IEEE transactions on medical imaging, 2018, 38(1): 145-155.

A Parametric Level Set based Approach to Difference Imaging in Electrical Impedance Tomography

本文通过引入形状先验信息，在一定程度上缓解了差分成像重建问题的病态性。受作品和的启发，我们在中扩展了之前的工作，在中开发了一种基于参数水平集(PLS)的绝对EIT方法。在这项工作中，我们考虑使用一种基于形状的方法来差分EIT基于PLS公式。在基于PLS的差分成像方法中，我们假设要重构的电导率变化可以被描述为一个或多个未知的分段常数电导率变化值的嵌入对象，并将异常的几何形状用高斯径向基函数(GRBF)表示。通过使用GRBF来表示偏最小二乘(PLS)函数，在用较少的未知数描述大量形状时提供了灵活性。这至少带来了以下好处:1大大减少了未知的总数量;2改进反问题的条件数;3提高了该技术的计算效率。为了解决基于pls的反问题，采用了一种旨在提高可靠重构电导率变化概率的迭代机制。

1、介绍

EIT的主要挑战之一是图像重建是一个严重非线性不适定逆问题;因此，需要有鲁棒的方法来处理噪声和建模误差，以及一个稳定促进的正则化策略。解决EIT问题的方法一般可分为绝对成像和差分成像。

1)、绝对成像

在绝对成像中，目标是根据单个状态对应的测量值来估计绝对电导率分布。重构绝对电导率分布需要精确的辅助模型参数信息，即电极位置、大小和形状，以及总体的体形或域边界。因此，大多数算法都假设电极信息和车身形状的边界是已知的，唯一的未知参数是电导率分布。在实际应用中，这些辅助模型参数的认识是不完整和不确定的，特别是在医学EIT应用中。例如，在胸部应用中，胸部形状会随着呼吸和测量过程中患者体位的变化而变化。结果表明，这些辅助模型参数的误差，会导致绝对重建中的出现严重的误差。

2)、动态（差分）成像

另一方面，在差分成像中，根据导电性变化前后采集的数据集，估算参考点导电性分布的变化。例如，在肺成像中，可以采集呼吸吸气时变化前的数据，呼气呼气时变化后的数据，这样差像显示了吸气和呼气时电导率的变化。

3)、成像方法

目前，解决差分EIT图像重建问题的算法有很多，如下文所述。

①基于逆散射问题[12]的分解方法[9]-[11]。因子分解方法对检测电导率变化是有用的;然而，在使用小数据集的情况下，它们的有效性可能值得怀疑，目前还不清楚如何纳入系统的先验信息。

②直接(非迭代)算法，如Block方法[13]和D-bar方法[14]-[19]，解决了完整的非线性EIT问题——不需要从正演模型中对电导率进行任何中间估计。Block方法假设待成像对象为二维矩形，由相同大小的固定块组成，在实际应用中存在一定的困难。另一方面，D-bar方法是根据电阻抗成像数据计算电导率的非线性傅里叶变换，然后反变换。该方法及其实现的详细描述，如[20]。对于讨论D-bar方法实用程序的不同应用程序，我们请读者参考[2]，[21]-[23]。为此，最近在[24]中也演示了Dbar方法的实时功能。需要注意的是，在D-bar方法中，正则化是通过非线性频域的低通滤波来实现的。然而，使用傅里叶数据的低通滤波通常会导致重建图像的清晰度下降。

③基于线性化的方法，如NOSER[25]、back projection method[26]和GREIT[27]，是应用最广泛的一些方法。在基于线性化的方法中，电导率与测量变化之间的关系通常是通过对观测模型进行全局线性化来建立线性化模型。线性差分重构在一定程度上可以容忍由辅助模型参数的不确定性引起的建模误差。当未知辅助参数与时间无关时，会出现该特征，从而在计算[28]测量值的差值时部分抵消建模误差。

2、EIT正向问题

设被研究的成像体占据一个二维或三维区域Ω⊂Rq, q = 2,3，其边界用∂Ω表示。EIT正向问题可以这样表述:给定感兴趣区域Ω内的电导率分布σ(x)和电流注入模式Iℓ，计算测量电极上的电势u(x)。在本研究中，所谓的完整电极模型(CEM)采用[43]。在数学上，CEM的正问题表述如下:

其中x∈Ω为空间坐标，zℓ为电极eℓ接触阻抗;Uℓ为电极eℓ测得的电势;L和n分别表示电极数和垂直于边界向量的一个向外单位。

要用CEM边界条件(2-4)数值解拉普拉斯方程(1)，需要使用数值技术。有限元方法[44]是常用的选择方法，本文也采用了这种方法。假设测量噪声是加性的、高斯的，EIT的观测模型可以写成：

其中V是测量电压的矢量，U(σ)是基于fem的正解，e是平均值e∗和协方差Γe的高斯分布噪声。

3、传统线性差分成像

在传统的线性差分成像方法中，初始状态σ1和最终状态σ2通过电导率变化∆σ线性相关;即σ2 = σ1+∆σ。最后的观测模型被写了出来：

其中ei ~ N(e∗i， Γei)， i = 1,2。请注意，噪声通常是固定的，即e∗i = e∗和Γei = Γe。本文也采用了这一模型。

差分图像重建是通过差分数据∆V = V2−V1重建状态间电导率变化∆σ = σ2−σ1。通常，∆σ的重构是基于(6&7)中观测模型的一阶泰勒近似：

其中J =∂U/∂σ(σ0)是根据背景导电性σ0∈R的初始猜测计算的正向映射的雅可比矩阵(灵敏度):

利用线性化和V2-V1得到观测模型:

电导率变化∆σ的Tikhonov正则化解为:

这里，L∆e定义为:

其中，噪声项∆e的协方差Γ∆e为：Γ∆e = Γe1 + Γe2 = 2Γe。术语p∆σ(∆σ)是一个正则函数，例如，最常用的ℓ2范数稳定的反演。

4、基于参数水平集的差分成像

本节介绍了一种基于参数水平集的差分成像方法来估计电导率的变化。为简单表示，设域名Ω存在一条边界Γ⊂Ω，将其划分为Ω+和Ω−两个部分，

即Ω=Ω+ U Ω-。我们考虑到电导率变化∆σ(x)是一个未知的电导率变化值 ∆σ(x) =∆σ1 ，X∈Ω+和∆σ(x) =∆σ0，X∈Ω-的分段常数函数，如图1所示：

图1 在二维空间中自由边界的对象和水平集表示的图示

对PLS技术在绝对电阻抗成像中的应用进行了综合研究。但是，为了方便读者，我们将在下面简要介绍它的基本理论和对差分EIT的扩展。

水平集方法将区域之间的边界Γ = {x: f(x) = 0}表示为高维函数f(x)的零轮廓，称为水平集函数(LSF)，需要满足：

在LSF f(x)中，∆σ表示为：

在式中，后一项表示异常，H(s)为Heaviside函数，当s < 0时H(s) = 0，否则H(s) = 1。

在实践中，为了能够数值解决反问题并使LSF的更新成为可能，人们通常使用Heaviside函数的平滑近似。其中一个选择是：

在这里，参数ε定义了一个宽度为2ε的频带，在这个频带内，Heaviside函数被平滑了[45]。

在大多数当代基于形状的方法中，LSF f(x)通常被选择为带符号的距离函数[38]，它与x-空间的离散化有关。现在考虑基集P = {p1, p2，···pN}中参数化表示的LSF f(x):

其中，N表示径向基函数(rbf) pi(x)的个数，µ=[µ1，µ2···，µN]∈RN是确定rbf权值的未知PLS参数向量。P基集的可能选择包括高斯、多二次、多谐样条和薄板样条多项式。和[42]一样，我们使用高斯RBF。也就是说，

其中γ为高斯宽度，xi为RBF中心，详情见V-D节，||·||为欧几里德范数。

根据LSF的参数化，可以将式(12)修改为：

这里，c是一个较小的正数。

根据式(17)，式(13)中电导率变化的分布可表示为：

实际上，该模型将未知区域Ω+映射到未知PLS参数µ的空间中，与传统的基于像素/体素的非线性迭代方法相比，极大地减少了已知问题的总未知数，显著提高了算法的效率。采用Tikhonov风格正则化[46]和平滑先验[8]的最大后验估计的改进Gauss-Newton算法。

(10)中的观测模型可以表示为：

然后，对基于PLS的差分成像中的分段常数∆σ0和∆σ1进行形状重建和估计，即可解决最小化问题：

这里I是单位矩阵。u∗是预定的常量值(参见章节V-D)。在不同成像中，∆σ∗j通常设置为0，因为∆σ的减小和增大是等可能的。(20)中增加的惩罚项可以提高算法的收敛速度，同时保持算法的稳定性。我们注意到，在最小化问题(20)中，未知的∆σ0和∆σ1被附加到未知的PLS参数µ和µ同时被估计。

A.水平集灵敏度分析

由于最小化问题(20)的非线性，需要迭代计算解，迭代过程中需要得到雅可比矩阵J(µ，∆σ0，∆σ1)。根据链式法则，PLS参数µµ的导数∆U(µ，∆σ0，∆σ1) w.r.t可以分解为三个偏导数的乘积：

式中，δ(·)为狄拉克函数。同样，我们得到：

为了解决式(20)中的极小化问题，我们采用了带有直线搜索的高斯-牛顿方法。

5、方法

本节通过数值模拟和实验测试了基于PLS的差分成像方法的性能。阐述了计算方法中的测试用例、估计、实现问题、参数选择和实验评估。有关结果和讨论，见第六节。

A、测试用例

为了研究提出的方法的性能，进行了以下研究。在每个测试用例中，模拟或收集两个测量集:V1是初始电导率σ1的测量参考集，V2是变化后的电导率σ2的测量参考集。注意，在本研究中不要求使用同质物体进行测量。

表1 电导率值基于从模拟中指定的工作[47]中获得的典型值

1)模拟胸部成像:首先，我们考虑三个对应于胸部成像不同条件的模拟测试案例1-3来演示我们基于PLS的差分成像方法的一般行为。如图4所示，对于情形1，初始状态σ1对应终止期相位，σ2变化后的电导率对应终止期相位。在案例2中，我们模拟了一个更真实的肺监测案例，在这个案例中，部分左肺萎陷。例3中，初始状态σ1模拟心脏在收缩末期的状态，变化后的传导率σ2模拟心脏在舒张末期的状态。值得注意的是，在模拟中，我们没有使用PLS函数来表示分配电导率分布的区域边界。相反，肺、心脏和主动脉的区域使用结构网来定义。

图4 利用传统的线性差分成像(E1)和基于pls的差分成像(E2)方法对模拟数据进行重建

目标形状是通过人体胸部计算机断层扫描得到的。在胸腔边缘放置16个长度为2cm的电极。在模拟中使用的组织的电导率列在表I中，测量数据是在振幅为1mA的相邻电流刺激和相邻测量下模拟的。为了模拟真实情况，我们在模拟数据中加入高斯噪声，其标准差为无噪声测量数据的最大值与最小值之差的0.1%。所选噪声电平对应信噪比SNR=42.63 dB。

2)水箱病例:在幻影实验研究中(病例4-7)，实验采用人胸形水箱，见图7。L = 16个宽度为2 cm的相同的矩形金属电极贴附在油箱的内侧面。水箱的水平和垂直半径分别为Rh = 17.5 cm和Rv = 14 cm。水箱里装满了自来水，容器里放置了不同材料(琼脂、塑料或铜)制成的物体，以模拟不同的电导率分布。为了使夹杂物易于阅读，将每个夹杂物的材料类型用大写字母标记，例如，用字母“a”表示“由琼脂制成”的夹杂物，详情见图7。所有的包裹体在垂直方向上都是均匀的，并通过水面延伸。使用KIT4系统[48]进行测量，采用相邻电流刺激(频率10 kHz，幅值1 mA)和相邻测量。

图7 利用传统的线性差分成像(E1)和基于pls的差分成像(E2)方法对水箱数据进行重建。幻影照片中标记的字母“A”、“C”和“P”分别表示由琼脂、铜和塑料材料制成的夹杂物