电阻抗成像(electrical impedance tomography，EIT)作为一种新型的功能成像技术^[1-2]，在体表通过电极施加弱的安全电流来测量体表引起的电压，然后将所获得的电压值送回计算机，依照一定的图像重建算法来重建关于被测对象内部的电导率(电阻抗)分布信息的断层图像。EIT具有无创无害、低成本、结构简单和功能成像及动态图像监护的优点^[3-4]，因而得到国内外众多研究者的广泛关注。

EIT图像重建的本质就是求解逆问题，即已知边界电压，求得目标区域内电导率(电阻抗)的分布，但是其技术困难在于逆问题的不适定性^[5]。因此在重建图像算法中通常采用正则化算法来改善其病态性，以得到适定的解^[6]。

EIT逆问题中最常用的正则化算法就是Tikhonov正则化算法。此算法的实质是构造特定形式的罚函数来实现对解的阻尼作用，从而达到使解稳定的目的。它的成功应用依赖于正则化参数的合适选取。目前选取正则化参数的方法有Arcangeli准则、偏差原理、岭迹法、广义交叉核实法(generalized cross-validating，GCV)和L曲线法等^[7-8]。使用前两种方法来确定正则化参数时，需要先对原始数据的误差水平做出估计。但是关于这些数据的先验性的信息仅便于理论分析，实际中参数的复杂分布导致其难于获取和验证。岭迹法具有直观的优点，但是缺乏严格的理论依据，并且具有主观性和随意性。GCV的缺点在于其函数的变化平缓，甚至难以收敛, 以至于求解的最优值不是全局最优值。L曲线法是目前广泛使用的确定正则化参数的方法，具有函数图像变化明显、定位准确的特点，在控制网评差、GPS快速定位和航天飞行器的精密轨道结算等方面已有应用^[9-10]，但是在有些情况下L曲线并没有拐点，因此在这种情况下L曲线法失效。

笔者提出了一种新的正则化参数的选取方法，通过极小化正则化算法中的残差范数(residual norm)和解范数(solution norm)的乘积(product of residual norm and solution norm，PRS)来求取正则化参数。理论上该方法较L曲线法更为简单。为验证该方法的有效性，结合二维场域模型，进行了相关的图像重构仿真实验和水槽实验，并与L曲线法进行了对比。结果表明，这种方法更具有优越性。

1 EIT数学模型

以二维EIT问题为例^[11]，EIT的图像重建问题实质上是一个低频电流场的逆问题求解，即已知边界电流和电压，求解内部的电导率分布。通常情况下，被研究区域的电流场被视为准静态场，该电场区域内的电位分布函数φ与电导率分布函数σ满足拉普拉斯方程：

$ \nabla \cdot \left[ {\sigma \left( {x,y} \right) \cdot \nabla \varphi \left( {x,y} \right)} \right] = 0,\left( {x,y} \right) \in \mathit{\Omega }; $

(1)

其边界条件为

$ \varphi \left( {x,y} \right) = f\left( {x,y} \right),\;\;\left( {x,y} \right) \in \partial \mathit{\Omega }; $

(2)

$ \sigma \left( {x,y} \right)\frac{{\partial \varphi \left( {x,y} \right)}}{{\partial \mathit{\boldsymbol{v}}}} = j\left( {x,y} \right),\;\;\left( {x,y} \right) \in \partial \mathit{\Omega }。$

(3)

式中：∂Ω为场域Ω的边界；f为已知边界电位；j为流入场域Ω的电流密度；v为场域边界∂Ω的外法向单位向量。电阻抗成像的正问题是在已知电导率σ分布和边界条件的情况下求取场域内的电位分布φ，而逆问题则是利用已知的边界条件(2)、(3)求解方程(1)中电导率σ的分布。

将场域离散成M个单元，电导率σ₀=(σ₁, σ₂, …, σ_M-1, σ_M)^T具有M个元素，其中每个元素代表每个单元的电导率。实际电流注入后，所测量的电极上的电压为V_ij(i, j=1, 2, …, N)，N为电极个数。使用有限元法计算每次电流注入时电极上的电压F(σ₀)，σ为求得的电导率分布。EIT的线性模型表示为

$ \mathit{\boldsymbol{F}}\left( {{\mathit{\boldsymbol{\sigma }}_0}} \right) - \mathit{\boldsymbol{V}} = \mathit{\boldsymbol{J}}\left( {{\mathit{\boldsymbol{\sigma }}_0} - \mathit{\boldsymbol{\sigma }}} \right), $

(4)

其中J为Jacobian矩阵。

应用最小二乘法的原理，将式(4)转换为

$ Z\left( \mathit{\boldsymbol{\sigma }} \right) = \min \left\| {J\left( {{\mathit{\boldsymbol{\sigma }}_0} - \mathit{\boldsymbol{\sigma }}} \right) - \left[ {\mathit{\boldsymbol{F}}\left( {{\mathit{\boldsymbol{\sigma }}_0}} \right) - \mathit{\boldsymbol{V}}} \right]} \right\|_2^2。$

(5)

电阻抗成像的逆问题就是寻找电导率σ的分布从而使得Ζ(σ)最小。

2 Tikhonov正则化及参数选取 2.1 Tikhonov正则化

EIT图像重构的难点是对严重病态的非线性逆问题进行求解，而该问题表现在对条件数很大的Jacobian矩阵求逆。因此对于这样的问题，在求解过程中往往需要利用正则化算法来降低Jacobian矩阵的条件数，改善病态性。

Tikhonov正则化算法是基于最小二乘法的原理，通过在目标函数中引入特定形式的罚函数来对Jacobian矩阵进行惩罚，以求得稳定的解。应用于EIT逆问题的Tikhonov正则化算法的形式^[12]为

$ \begin{array}{*{20}{c}} {E\left( {{\mathit{\boldsymbol{\sigma }}_0} - \mathit{\boldsymbol{\sigma }}} \right) = \min \left\{ {\left\| {J\left( {{\mathit{\boldsymbol{\sigma }}_0} - \mathit{\boldsymbol{\sigma }}} \right) - } \right.} \right.}\\ {\left. {\left. {\left[ {F\left( {{\mathit{\boldsymbol{\sigma }}_0}} \right) - \mathit{\boldsymbol{V}}} \right]} \right\|_2^2 + {\lambda ^2}\left\| {\mathit{\boldsymbol{L}}\left( {{\mathit{\boldsymbol{\sigma }}_0} - \mathit{\boldsymbol{\sigma }}} \right)} \right\|_2^2} \right\}。} \end{array} $

(6)

式中：λ＞0，为正则化参数，控制着右边残差范数‖J(σ₀-σ)-[F(σ₀)-V]‖和罚函数‖L(σ₀-σ)‖的相对比重；L为正则化矩阵，通常取L=I。

2.2 正则化参数的选取

令F(σ₀)-V=B，δσ=σ₀-σ，则式(6)可以表示为

$ \mathit{\boldsymbol{E}}\left( {\mathit{\boldsymbol{\delta \sigma }}} \right) = \min \left\{ {\left\| {\mathit{\boldsymbol{J}} \cdot \mathit{\boldsymbol{\delta \sigma }} - \mathit{\boldsymbol{B}}} \right\|_2^2 + {\lambda ^2}\left\| {\mathit{\boldsymbol{\delta \sigma }}} \right\|_2^2} \right\}。$

(7)

对式(7)求导可得到解为

$ \mathit{\boldsymbol{\delta \sigma }} = {\left( {{\mathit{\boldsymbol{J}}^{\rm{T}}}\mathit{\boldsymbol{J}} + {\lambda ^2}\mathit{\boldsymbol{I}}} \right)^{ - 1}}{\mathit{\boldsymbol{J}}^{\rm{T}}}\mathit{\boldsymbol{B}}。$

(8)

设J∈R^p×q，秩为r，对J进行奇异值分解得：

$ \mathit{\boldsymbol{J}} = \mathit{\boldsymbol{US}}{\mathit{\boldsymbol{V}}^{\rm{T}}} = \sum\limits_{i = 1}^r {{\mathit{\boldsymbol{u}}_i}{\alpha _i}\mathit{\boldsymbol{v}}_i^{\rm{T}}} , $

(9)

式中：U和V分别为S的左右奇异矩阵，且均为正交阵；u_i和v_i分别为α_i在U和V中对应的向量；S为由J的r个奇异值为对角元素组成的矩阵，且α₁≥α₂ ≥…≥α_r-1≥α_r＞0。则式(8)可写成：

$ \mathit{\boldsymbol{\delta \sigma }} = \sum\limits_{i = 1}^r {\frac{{\alpha _i^2}}{{\alpha _i^2 + {\lambda ^2}}}\frac{{\mathit{\boldsymbol{u}}_i^{\rm{T}}\mathit{\boldsymbol{B}}}}{{{\alpha _i}}}{\mathit{\boldsymbol{v}}_i}} 。$

(10)

称$\frac{{\alpha _i^2}}{{\alpha _i^2 + {\lambda ^2}}}$为滤波函数。显然，当λ→0时，Tikhonov正则化求得的解无限接近于基于最小二乘法原理建立的目标函数，即式(5)的解。

可见式(7)右边两项均是λ的函数，即

$ \begin{array}{l} \rho {\left( \lambda \right)^2} = \left\| {\mathit{\boldsymbol{J}} \cdot \mathit{\boldsymbol{\delta \sigma }} - \mathit{\boldsymbol{B}}} \right\|_2^2\\ \;\;\;\;\;\;\;\;\;\; = \sum\limits_{i = 1}^r {\frac{{{\lambda ^2}B_i^2}}{{{{\left( {\alpha _i^2 + {\lambda ^2}} \right)}^2}}}\left\| \gamma \right\|_2^2} , \end{array} $

(11)

$ \eta {\left( \lambda \right)^2} = \left\| {\mathit{\boldsymbol{\delta \sigma }}} \right\|_2^2 = \sum\limits_{i = 1}^r {\frac{{\alpha _i^2B_i^2}}{{{{\left( {\alpha _i^2 + {\lambda ^2}} \right)}^2}}}} 。$

(12)

其中$\gamma = \sum\limits_{i = r + 1}^p {{B_i}} \cdot \delta {\sigma _i},{B_i} = \mathit{\boldsymbol{u}}_i^{\rm{T}} \cdot \mathit{\boldsymbol{B}}$。而寻找λ的目的就是要使得式(7)中右边两项残差范数‖J·δσ-B‖₂和解函数‖δσ‖₂同时为小。

2.2.1 基于L曲线法的正则化参数选取

L曲线法是以正则化参数λ为参变量，通过横纵坐标的对比来确定正则化参数。它是基于残差范数‖J·δσ-B‖₂和罚函数‖δσ‖₂的比值的曲线图。曲线用log-log坐标，具有明显的“L”形状。

Hanke等认为L曲线的关键是定位曲线上曲率最大的那个点，选择其对应的λ为所求的正则化参数。设横坐标和纵坐标分别为

$ \hat \rho \left( \lambda \right) = \mathit{lg}{\left\| {\mathit{\boldsymbol{J}} \cdot \mathit{\boldsymbol{\delta \sigma }} - \mathit{\boldsymbol{B}}} \right\|_2},\hat \eta \left( \lambda \right) = \lg {\left\| {\mathit{\boldsymbol{\delta \sigma }}} \right\|_2}; $

(13)

根据曲率公式：

$ \kappa = \frac{{\hat \rho '\hat \eta '' - \hat \rho ''\hat \eta '}}{{{{\left[ {{{\left( {\hat \rho '} \right)}^2} + {{\left( {\hat \eta '} \right)}^2}} \right]}^{3/2}}}}, $

(14)

式中：$\overset{\hat{\ }\prime }{\mathop{\rho }}\,、\overset{\hat{\ }\prime }{\mathop{\eta }}\,、\overset{\hat{\ }\prime\prime }{\mathop{\rho }}\,、\overset{\hat{\ }\prime\prime }{\mathop{\eta }}\,$分别为$\hat{\rho }、\hat{\eta }$的一阶导数和二阶导数；计算出最大曲率κ_max，κ_max所对应的点处的λ即为L曲线所确定的正则化参数。

2.2.2 基于PRS法的正则化参数选取

正则化参数的选取应该综合考虑残差范数和解范数对目标函数的贡献，通过减小Jacobian矩阵较小的奇异值对噪声的影响来控制解的光滑程度。在大多数情况下，根据文献[13], 可以通过极小化泛函来实现。

$ \mathit{\Phi }\left( \lambda \right) = {\left\| {\mathit{\boldsymbol{J}} \cdot \mathit{\boldsymbol{\delta \sigma }} - \mathit{\boldsymbol{B}}} \right\|_2} \cdot {\left\| {\mathit{\boldsymbol{\delta \sigma }}} \right\|_2}。$

(15)

Engl^[13]等指出，在L曲线的拐点C($\hat \rho$(λ_*)，$\hat \eta$(λ_*))处，切线的斜率为-1。通过式(13)和式(15)，可以推导出Φ(λ)=10 ^{[$\hat \rho$(λ)+$\hat \eta$(λ)]} 。因此如果Φ(λ_*)在λ=λ_*处有局部极小值，则存在

$ \left[ {\hat \rho '\left( {{\lambda _ * }} \right) + \hat \eta '\left( {{\lambda _ * }} \right)} \right] \cdot \mathit{\Phi }\left( {{\lambda _ * }} \right) = 0; $

(16)

即$\overset{\hat{\ }\prime }{\mathop{\rho }}\,\left( {{\lambda }_{*}} \right)+\overset{\hat{\ }\prime }{\mathop{\eta }}\,\left( {{\lambda }_{*}} \right)=0$。该式表明，L曲线在C点处的切线与向量(1，1)^T正交。

因此通过对式(15)的求解，可以求得λ的非劣解。对于任何λ，该方法不需像L曲线计算导数和曲率，大大地减少了计算量，节省了运算时间，易于实现。

3 重建结果比较 3.1 模型及设定

建立如图 1所示的有限元二维模型，内有1 090个单元、578个节点，采用16电极相邻注入模式。

3.2 计算重建

仅仅观察重建图像，只能定性地感觉图像的好坏，不能定量评价图像的优劣，因此，依据图像的相关度(correlation degree，CD)^[14]，残差范数(residual norm，RN)和解的误差范数(error norm，EN)^[15]来综合评价不同λ所重建出的图像质量。其中，图像相关度的定义为

$ {\rm{CD}} = \frac{{\sum\limits_{i = 1}^M {\left( {{\sigma _{{\rm{true}} - i}} - \bar \sigma { - _{{\rm{true}}}}} \right)\left( {{\sigma _{{\rm{rec}} - i}} - {{\bar \sigma }_{{\rm{rec}}}}} \right)} }}{{\sqrt {\sum\limits_{i = 1}^M {{{\left( {{\sigma _{{\rm{true}} - i}} - {{\bar \sigma }_{{\rm{true}}}}} \right)}^2}} \cdot \sum\limits_{i = 1}^M {{{\left( {{\sigma _{{\rm{rec}} - i}} - {{\bar \sigma }_{{\rm{rec}}}}} \right)}^2}} } }}, $

(17)

式中：σ_true-i、σ_rec-i为第i个单元的真实电导率和重建计算出的电导率；σ_true、σ_rec分别为所有单元的真实电导率与重建计算出的电导率的平均值；M为总单元数。

相关度代表重建的图像与真实图像的相似程度^[16]。因此，相关度的值越大，表明重建图像与真实图像越接近，重建的图像质量越好。残差范数和解的误差范数均表示重建图像与真实图像的误差，理论上这两个值越小越好。

3.3 仿真结果

为了验证所选择的正则化参数方法的适用性，采用上述模型，利用笔者提出的优化方法和L曲线法进行对比。针对目标区域的电导率、位置、大小以及噪声等不同情况，将仿真实验分为5组。对用式(14)和式(15)所求的正则化参数(分别用λ_LC和λ_PRS表示)重建出的图像做比较，并根据上述3个评价指标对图像进行了评价。

3.3.1 改变目标区域位置

设置背景区域电导率为3 S/m，目标电导率为1 S/m，目标区域的中心位置分别位于(8, 0)，半径约为2，如图 2(a)所示。根据求得的λ_PRS和λ_LC重建的图像分别如图 2(b)和(c)所示。

需要说明的是，笔者采用的是动态成像，因此灰度条的数值为相对灰度值，并非为绝对电导率值。根据前面给出的3个评价指标来评价它们的质量。另外了还对比两种方法所用的时间，结果见表 1所示。

从上表可以看出，根据该优化方法选取的λ重建的图像与真实图像的相关度较大，同时残差范数和误差范数以及计算时间也较小。因此，利用这种方法能够快速找到相对最优的λ。

改变目标区域的位置至(0, 0)，其他参数不变，根据不同的λ重构出的图像和评价结果如图 3和表 2所示。

3.3.2 改变目标区域大小

与上文3.3.1中图 2对应，设置背景区域电导率为3 S/m，目标区域电导率为1 S/m，目标区域中心位于(6, 0)，半径约为4。根据不同的λ所重建出的图像分别如图 4所示，评价结果如表 3所示。

3.3.3 改变目标区域电导率

保持目标区域中心不变，仍位于(8, 0)，半径约为2，如图 2(a)所示。设置背景区域电导率为3 S/m，目标区域电导率分别为6 S/m。根据不同的λ所重建出的图像如图 5所示，评价指标的计算结果如表 4所示。

3.3.4 改变目标区域的数目

设置背景区域电导率为3 S/m，目标区域电导率为1 S/m，双目标区域的中心位置分别位于(-8, 0)和(8, 0)，半径分别约为2。根据不同的λ所重建出的图像分别如图 6所示，评价结果如表 5所示。

3.3.5 改变数据噪声强度

为了测试该优化方法的抗噪性，在原始数据的基础上添加不同强度的的高斯白噪声。根据不同的噪声强度，依据该方法计算出相应的λ，再利用所求的λ重建出图像。设置背景区域电导率为3 S/m，目标区域电导率为1 S/m，其中心位置位于(0，0)，半径约为2 (如图 3(a)所示)。重建的图像如图 7所示，同样，评价结果如表 6所示。

3.4 水槽实验

在实际物理测量中，不可避免会受到多方面因素影响，如硬件系统的测量噪声、工频干扰等，这直接影响最后的重构图像质量。为了验证该方法的有效性，笔者利用其所在团队开发出的二维水槽实时电阻抗成像系统所采集的电压数据，进行了相应的水槽实验，如图 8(a)所示。在水槽中装满水，放入一铝棒，针对图 8(b)两种不同的目标位置，与仿真实验类似，分别利用笔者所改进的方法和L曲线选取的λ进行图像重构，如图 8(c)-(d)所示。

需要说明的是，上图中，深色表示高电导率，浅色表示低电导率。由于铝棒是良导体，其电导率比水要高很多，因此在上图中以黑色表示。

3.5 分析和讨论

根据3.3节中前4组仿真实验的重建图像和评价结果，可以看出，当不存在噪声时，根据笔者提出的优化方法求取的λ重建出的图像质量普遍较好，用时较少, 能较好重现目标区域。从而验证了这种正则化参数选取方法的适用性。但需说明，第5组仿真实验中，随着噪声的增加，所重建的图像效果呈恶化趋势。当噪声强度为5%时，图像已经完全失真，背景区域扰动严重，出现有多个虚假的目标，重建图像失败。这是由于Tikhonov正则化算法采取的是2范数的函数为正则化函数，侧重于对适定性问题的解决，而忽略了对逆问题由此产生的有平滑效应，使得重建的图像呈平滑状态，空间分辨率不高，干扰噪声在平滑的过程中被传播和扩散。

在对铝棒的水槽成像实验中，从重构图像中能基本识别出铝棒的位置。横向对比图 8(c)-(d)，目标靠近边缘时的图像比目标接近中心区域的图像要好，重构出的目标大小更接近实际目标的大小。纵向对比可以发现采用笔者提出的方法PRS法所重构的图像质量比用L曲线法重构的图像高。但是图像边缘处不均匀，均存在伪影，造成以上现象有以下3点原因：1)硬件测量系统本身的噪声和外界环境中的干扰等因素。2)有限元剖分误差，为了保持全文的一致性，在水槽实验部分采用与仿真实验相同的有限元剖分模型。该有限元模型的剖分密度低，同时利用有限元进行正问题计算时存在模型误差。3)EIT的逆问题是欠定的，成像区域内部大的电导率变化在边界可能产生非常微小的变化，且边界测量值对边缘附近，尤其是电极附近的电导率变化很敏感，而对中心附近电导率变化不敏感，这就导致了目标在中心区域的重构图像质量没有目标在边缘处的图像质量好。

最后，笔者对这种方法作了不适定分析。参数设置如前文3.3.1中相同，λ_PRS=3.504 5e-9，使用Tikhonov正则化后的Jacobian矩阵的奇异值分解图如图 9所示。

从上图可以看出，对Jacobian矩阵进行正则化处理后，奇异值仍处于分散的状态，这说明正则化前后Jacobian矩阵的性质没有发生大的变动。经计算，当λ_PRS=3.504 5e9时，条件数由正则化前的7.672 9e+21减小到2.466 0e+4，说明原Jacobian矩阵的不适定性得到了明显的改善。

4 结语

笔者基于Tikhonov正则化算法，结合EIT的模型，利用极小化正则化算法中的残差范数和解范数的乘积的方法来求取正则化参数。并将这种方法引入到EIT图像重构算法中以解决正则化参数优化问题。通过相关的仿真实验和水槽实验，分析比较了根据该方法选取的λ和L曲线法选取的λ分别重建出的图像质量。实验结果表明，该方法具有简单实用和计算时间短的优点。采用所提出的优化方法选取的正则化参数重建的图像质量较好，能够较好地再现目标区域。这为EIT逆问题的病态方程的研究提供了一种实用而有效的参数优化方法，避免了以往选取正则化参数的随意性和主观性，提高了成像质量。