基于模糊算子的电学层析成像算法

岳士弘，荣西拉，马海涛，卢 剑

(天津大学电气自动化与信息工程学院，天津 300072)

电学层析(electrical tomography，ET)是一种新型、先进、非侵入式的可视化测量技术，它具有成本低、响应快、安全性高等优点，在工业和医学领域具有广阔的应用前景[1-2]．一个典型的ET 系统主要通过被测场域边界电极阵列的轮流激励和测量；测量值通过信号数据采集单元；最后上传到上位机进行图像重建与显示可视化．一个ET 系统结构简单易于实现．然而，由于ET 过程固有的“欠定性”和“软场效应”[3-4]，其成像的空间分辨率很低，难以满足很多应用中的基本需求．在大多数情况下，ET 成像算法依赖于一个基于目标函数的优化问题．从最先提出的核心优化算法——高斯-牛顿(Gaussian-Newton，GN)法[5-7]、吉洪诺夫正则化(Tikhonov regularization，TR)[8]法到最近的基于Lp范数各种变形形式[9-11]以及Bregman 散度[12]，它们都有各自可应用的范围以及明显的局限性．如何面向一个具体的应用对象，采用有效的方法提高ET 的空间分辨率和稳定性，几十年来一直是该领域一个重点和难点问题．

针对上述问题，本文基于模糊算子[13]提出一个一般的ET 成像算法．它具有两个优势：①具有近百个功能各异的模糊算子可供使用者选择，特别是目前ET 中目标函数常用的加权和形式只不过是模糊算子的一个特例；②一个易于实现的非传统优化的求解方法；比较而言，很多已有的ET 优化算法往往要求可导或至少是连续的，面临着如何优化和实现的问题，而本文方法只要求目标函数在可行域上每点有唯一的函数值就可以实现优化；这种要求更加符合工程实际．最后，通过重建图像质量和评价参数对方法的有效性进行了验证．

1 研究背景与相关工作

1.1 电学层析原理

ET 技术主要分为电阻层析成像[14]、电容层析成像[15]和电磁层析成像[16]．这3 类技术虽存在一定差异，其重构原理是相似的．本文仅以二维电阻层析成像为例进行说明．

将调查的场域采用有限元刨分和离散化，并将边界测量与内部场域关系进行线性化后得

式中：m 为被测场域内的独立测量数；n 为像素点个数；Jm×n为雅可比矩阵，表示测量电极对每个像素点单位电导率增量的响应程度或灵敏度；σ 为被测场域的电导率分布．记U为边界测量构成的向量，则归一化式(1)后得到

该方法求解速度快并且易于实现．但是方程(2)对应的逆问题求解本质上是非线性、欠定和病态的，因此LBP 算法很难获得准确的解．大多数ET 算法把问题转化为以下优化问题，即

式中||·||为任意范数．过去大多数研究集中在使用不同形式的范数对式(4)进行求解；但另一种改变是式(4)左边增加限制项，其中最典型的是TR 算法，它具有下列形式，即

式中λ 为正则化系数．

这些方法对于具体应用往往有特有的优势．但是，这些方法的使用往往带有一定的假设且泛化能力不强难以移植[18]．此外，对于任意一个目标函数如何求解也是关键问题[19]，克服这些问题是当前ET 研究领域的重点和难点．

1.2 模糊算子

本文模糊算子是指基于模糊关系的合成运算[20].模糊算子种类丰富而多样，在模糊推理过程中，算子选取的不同，会导致推理方法的可适用范围和推理结果的不同[21]．对于任意a，b∈[0，1]，模糊算子T(a，b)(常称为T模)必须满足3 个条件：

常用的模糊算子有以下几种：

上述“max”算子就是在两个变量之间取较大的一个，这个算子应用最为简单和方便，但容易丢失信息．相比于max 算子，式(7)中其他算子能更多地保留运算前的信息．如T2和T4实际上是求两个向量的夹角正弦和正切值等；它们对于模糊推理和模糊运算能够起到不同的作用．

2 基于模糊三角算子的图像重建算法

2.1 基本动机和模型

本文提出算法的核心思想是根据测量值找到一个合适的模糊算子，通过这个模糊算子代替式(2)中的求和运算，实现下述一般化转换，即

式中：角标k和j对应第k个方程(测量)以及第j个变量，[skj]∈S；“·”表示对于任何两项使用模糊三角算子做运算．使用式(8)的主要动机包括以下几个方面．

(1) 由于软场效应和非线性，S 中每个元素skj所反映的任意第k个测量值uk对于第j个像素σj具有高度不确定性，从而看成模糊量是必要和可行的；将S 中的值离散和模糊化之后，一定程度上能反映灵敏度系数本身不准确性以及测量值中噪声的影响，具有较好的可解释性．

(2) 大量具有不同可适用范围和功能的模糊算子为式(8)的表达和拓展提供了丰富的选择．实际上，式(2)相当于模糊算子中取T3形式的一个特殊情况；相对于一般的工程需求和复杂工况，有潜力实现更好和更有效的计算．

(3) 式(8)中取模糊算子构造的任何形式目标函数都必须通过优化才能实现电学层析过程．在第2.2节将使用一个一般的优化方法使式(8)的求解和实现成为可能．

(4) 工程中很多先验信息是用语言规则表示的.模糊系统正是基于经验和规则进行推理和分析的，能够较好地将先验信息和工程实践经验导入模糊系统，从而抑制噪声并提高系统的鲁棒性．特别是对于ET系统面临的具有不一致、不确定和不完整特点的工程问题，精确模型及其优化往往耗费极大，而模糊模型能够提供更好的选择．

2.2 实现算法

假设ET 的边界测量值已知且灵敏度矩阵S 已经预先确定．基于一般的模糊算子对应的ET 任务是如何求解最优电导率的分布，使得式(8)对应的目标函数误差最小化．

假设过程的m个测量值uk的被归一化为一组离散值a1，a2，…，am，满足0＜a1＜a2＜…＜am-1＜am＜1，那么求解式(8)右边的目标函数等价于近似求解一组σ1，σ2，…，σn的值，使得

式中b1，b2，…，bm是在0，a1，a2，…，am之间的任意插值，满足0＜b1＜a1＜b2＜…＜am＜bm． 对于满足此条件的ai和bi和任意的模糊算子 “·”，必然存在一个严格增函数f[22]，满足f(0)＝0，且

式(10)的意义在于：对于式(9)中任意模糊算子的运算结果，都可以通过一个一元增函数f及其简单的加法运算代替来计算．

如果f(br)的值均为已知，那么为了求出σj的值，必须根据式(10)求解未知量f(skjσj)，j＝1，2，…，n．这类不等式有很多求解方法，本文使用一个易于实现和操作的纠错算法[23]来实现．

式中[s/2]、[(s＋1)/2]中的符号[·]表示向下取整，即不超过s/2 以及(s＋1)/2 的最大整数．

如果设Qk=(q1, q2, …,q2m)，Rk=(r1, r2, …, r2m)，W =(u1, u2, u3, …,u2m)，Es表示第s列为1 其余列均为0 的1×2m单位矩阵，式(12)、式(16)和式(17)的求解可以被写为

求式(18)的一个一般解W*的迭代步骤为

式中：μ 是一个权重系数，若无先验信息可取为1；Ai是将Qk、Rk、Es按照迭代循环排列后的向量；||Ai||是取定的矩阵标准范数．通常序列Ai的迭代可以是无限的，但该算法必然会在有限的次数内收敛．事实上，式(19)的迭代次数不超过上限[23]

式中f-1为单调函数f的逆算子．为了最终求出各个像素的σ分布值，必须将单调函数f具体表达出来．

考虑到ET 过程的近似性，可以与已知的测量值uk构造数据点列(uk，f(uk))进行样条插值[24]，i＝1，2，…，m，设获得的插值函数为F(x)，而作为一个单调单元函数，它的逆函数F-1(x)值即为F(x)．它也是

f--1(x)的近似表达式，这里取为F(x)的倒数．因此，F-1(x)的表达式为

最后，各个像素的σ的分布值为

由于灵敏度系数skj的值是完全确定的，式(25)作为选择ET 中关键的并算子一般通式．

当一组模糊算子可以适用和选择时，由式(23)可以定义不同的结果，为了检测哪个算子是最适合的结果，可以通过下面的式子求解，即

使得式(26)获得最小值的模糊三角算子为最优解．

最后，当任意一个模糊三角算子给定后，基于上述各个公式，可以得出以下关于电学层析的新优化算法基本实施步骤．

步骤1获得边界电压测量数据并归一化后升序排列为u1，u2，…，um；

步骤2确定空场的灵敏度系数skj，k＝1，2，…，m；j＝1，2，…，n；

步骤3根据式(16)和式(17)计算值；

步骤4根据式(19)求解不等式(18)得到一个最优解W*；

步骤 5根据式(22)求解f(skjσ)和f(bj)；

步骤 6根据式(25)计算电导率分布．

3 仿真实验和结果

利用有限元仿真软件COMSOL 对场域建模，圆形场域半径为10 cm，场域表面分别安放16 个电极，电极长为1 cm，背景电导率为1 S·m，目标电导率为0.01 S·m．原始模型和重建图像如表1 所示．采用常用的相邻激励和相邻测量方法，可获得的独立样本数为208 个，为此场域刨分为812 个单元．

表1 中第1 列对应原型，其中红色像素构成目标，蓝色区域像素对应背景．这些原型又两两分为3组：1 组为两个目标离散的模型；2 组为人体胸部截面模型；3 组为两个目标连续的模型；这些原型具有不同的代表性．第2 到第6 列为式(10)中采用模糊算子T2、T4、T5以及使用TR 得到ET 图像；其中T3相当于使用了LBP 算法；而TR 算法的结果受正则化参数λ 的影响．为此，实验中遍历了参数λ在一定范围内的所有可能值，并取了对应值下最佳空间分辨率的重构结果进行比较．

本实验实施过程中做以下设定．

(1) 归一化处理．对灵敏度系数矩阵S 与测量值U 相乘作为整体进行归一化，即先把灵敏度系数矩阵的每一行都乘以该行对应的最大值，然后再进行归一化．由于作为求解依据的测量值不是完全准确的而是存在误差，对于模拟气泡分布比较特殊的样本，求解得到的W 是有可能出现个别元素有异常值，因此需要足够数量和不同分布的样本来避免这些偶然因素．

(2) 离散化的尺度．为使求解不等式不过于复杂同时照顾实际问题中的精度要求，aj和bj各使用30 个不同的但等间隔的取值，即aj＝1/30，1/15，1/10，2/15，…，1；bj＝1/60，1/20，1/12，7/60，…，59/60．考虑算法的边值的范围，要求aj和bj的每个值不能等于 0 或等于 1，最后一个 a 值可改为0.999 9．这相当于输出的所有灰度值可能为30 个不同的离散值，即把连续的灰度值转化为30 个不同的灰度等级，这对于仅有812 个像素的ET 应用于实际问题是足够的．

(3)W的初值为全1 向量．虽然理论上求解W的纠错算法必然会收敛，但实际收敛的速度因不同样本而异可能会十分缓慢．为保证效率，规定序列Aj经过1 000 个完整周期之后，必须终止．同时，给每一次W的修正量取不同的修正权重μ值，随着序列Aj的循环次数不断增加，每次检测收敛条件时，W的修正量不断减少．第i 次循环的权重μi为100/(i＋99)．例如序列Aj已经在最后一个周期时，W的修正能力只有W1的修正能力的1/10．

(4) 阈值设定．为了避免异常值，先计算前50个样本求解得到的W向量的每一个元素的平均值和标准差E，此后对于每一个新样本得到的新的W向量，计算每一个元素和该元素平均值的偏差，如果超过5E 则视为无效值．

表1 显示了仿真实验所测试的6 个原型在不同的模糊算子下的成像结果，按照3 组不同的类型具体说明如下．

(1) 对于一组中目标离散分布的模型，可以直观地观察到T2和T5目标成像更加接近于原型，T4次之而T3(LBP)成像结果总体是最差的，这是由于LBP不包含任何迭代过程，成像结果最差，因此，T2和T5更具有分辨离散分布目标的能力．

表1 使用模糊算子T2、T3、T4、T5 的ET成像对比Tab.1 Comparison of ET images by using fuzzy operators T2，T3，T4，and T5

对于TR 算法而言，其重构图像中的目标伪迹很小，不仅目标位置正确而且边缘清楚．

(2) 对于2 组中肺部界面模型而言，由于设置了3 种电导率分布的目标，它们分别对应肺部组织(包含2 个和4 个结节)、心脏和外部骨骼，导致ET 成像结果都不太理想；各个组织的边界都不够清楚，几乎无法发现心脏；但是，模拟医学诊断过程的结节能够一定程度上呈现．

(3) 对于最后两个目标联系的界面模型而言，分辨率是有边界的清晰度决定；其中，第2 个模型对应有两层边界；从ET 结果图中可以看出，基于T3的算法实际输出与原型比较更加接近，所有边界基本上可以辨识出来．这说明T3(LBP)在重构真实的边界时是有优势的．比较而言．基于模糊算子T4的结果空间分辨率比T3有所提高，但仍然存在着大量伪迹；TR 重构的界面是最不清楚的．

上述4 个模糊算子对应的空间分辨率可以用相关系数λR(correlation coefficient)和RE(relative error，相对误差er)指标来进一步定量说明和评价，分别定义为

式中：G和G′分别表示重建图像和实际图像的灰度值向量；σi和σi*分别为σ 和σ*的第i 个像素；σ0和σ0*分别为σ 和σ*的平均值；n 为表示像素总数．er越小表明ET 算法重构的图像越接近真实图像，从而成像分辨率越高；er值在极端情况下不小于0．λ指标越大表示重建图像的灰度与实际灰度相关程度越大，重建图像的质量越高．

表2 显示了3 个算法的λ值与er值．显然，模糊算子的λR值与er值显示结果与它们的可视化图像是一致的．依据不同算子和模型，er值最大可以减小0.146 2，且平均减小0.093 1．总体而言，当目标离散时，T3对应结果是最劣的．而T5目标函数比模糊算子的目标函数值大近一个数量级．当目标连续时，T3对应结果是有优势的．对于像模拟人肺部断面这样的复杂分布的情况，几种算子各有特点；这说明不同的模糊三角3 种有各自可适用的范围．TR 算法由于带有约束项，对于离散的小目标辨识能力很高；但是对于连续的目标边缘(如两种介质边界面)则几乎失效，而T3(LBP)算法则正好相反．使用Matlab 中的Tic 和Toc 函数对6 个模型统计得：T2～T5的平均运行时间基本是一致的，分别为0.201 s、0.223 s、0.230 s和0.219 s，TR 由于要求逆矩阵时间为0.332 s，略长于其他方法．

表2 使用模糊算子T2、T3、T4、T5 和TR的精度比较Tab.2 Comparison of accuracy by using fuzzy operators T2，T3，T4，T5，and TR

4 结 语

本文将模糊理论中的模糊算子运用在ET 核心方程组的求解上．首先从ET 系统的欠定性和病态性出发，分析了ET 图像重建中采用模糊算子代替矩阵乘法中的求和运算的理论依据与合理性，然后针对ET 系统的特点，提出了一种可代替传统求和运算的模糊算子的算法，并给出了采用模糊算子的ET 成像方法，最后通过仿真实验，证明了不同的模糊三角算子对于不同的可视化目标有各自适用范围，从而拓宽了ET 技术的可应用范围并提高了针对性．特别情况下，已有的ET 算法所常用的加和形式只不过是本文算法取特殊模糊算子的一个特例，证明本文算法具有一定泛化能力．

考虑到模糊算子的一般性和普遍性，本文仅为提高ET 成像的空间分辨率提供一条新的途径或者说打开一个新的窗口；所使用的实验验证也是有限的.本文将来的研究聚焦在真实和复杂环境中，探讨本文提出方法的可应用范围．