融合拉普拉斯变换下信号稀疏性的电容层析成像图像重建算法

张 政，吴 阳，车权齐，刘 凯，姚佳烽

(南京航空航天大学机电学院，江苏 南京 210016)

电容层析成像(electrical capacitance tomography，ECT)是一种通过测量边界电容的层析成像方式。由于ECT技术成本低、安全性高、无创性[1]好，使得其在气固两相流的流型识别、浓度分布测量[2]等工业应用中得到了广泛应用。虽然ECT图像重建算法发展迅速，其中Landweber算法[3]更是得到广泛应用，但是因为ECT图像重建的反问题具有病态性，所以很难获得高质量的重建图像。因此，在图像重构过程中，需要将先验知识作为正则项[4]加入到图像重构问题中。目前，典型的正则化重建算法有Tikhonov正则化算法、FOT正则化算法[5]和总变差算法[6]。

Tikhonov正则化算法的重建速度非常快，但重建质量不理想。FOT正则化算法则会使重建图像的边缘过度平滑[5]。总变差算法具有较好的图像重建效果，但计算复杂。当然，还有许多用其他方法设计的正则项，如结构感知先验[7]和机器学习[8]，但是它们通常局限于特定的框架。

本文将重构图像的边缘信号的稀疏性这一先验知识融入到逆问题求解中。稀疏信号的重建常常使用压缩感知理论(compressed sensing，CS)[9]，但是因为测量矩阵的限制等距性(RIP)难以验证，所以将边缘信号的稀疏性转变为1范数约束融入到ECT图像重构之中。因此，图像处理问题中的稀疏约束被期待有着更好的优化结果。

1 拉普拉斯变换下重构信号的稀疏性

1.1 三角网格下的拉普拉斯变换矩阵

在ECT图像重建过程中，本文采用如图1左侧所示的三角网格，该网格由4 128个三角形、2 129个点元素组成。图1右侧为网格中编号为i的三角形的示意图，图中标注了每个点上的横纵坐标x,y以及相对介电常数的变化值δε；li13,li32和li21代表第i个三角形三条边所表示的向量；θi1,θi2和θi3代表第i个三角形的3个角度。

图1 三角网格剖分图与第i个三角单元

使用余切权值表示的拉普拉斯转换矩阵[10]如下：

(1)

式中：δεi为第i个网格点上的相对介电常数变化量；L(δεi)为向量δε在拉普拉斯变换后在i点处的数值；N(i)为和第i个网格点相邻的点的标号集合；αij和βij分别为包含第i个和第j个网格点的线段在两个不同三角单元中的对角。

ECT逆问题可以简单地描述成式(2)：

δC=Sδε+n

(2)

式中：δC为测量电容值变化向量；δε为 相对介电常数在各个网格点上的变化量组成的向量；S为敏感矩阵，此处采用电容-相对介电常数分布关系函数在气相状态处的斜率来表示灵敏度[4];n为噪声向量。

FOT正则化[5]如式(3)所示：

(3)

式中：μ为正则项系数；K为转换矩阵。FOT正则项如式(4)所示：

(4)

式中：Ω为整个求解域，即所有三角形网格区域；x和y分别为积分域内点的笛卡尔坐标系下的横纵坐标。

在第i个三角形上有：

(5)

式中：Ki为第i个三角单元中的转换矩阵；Ωi为第i个三角形网格区域。

相对介电常数在第i个三角单圆上的线性插值表达式为式(6)，可以推得式(7)，其中a，b和c是由第i个三角形顶点的网格坐标决定的常数。

δεi(x,y)=ax+by+c

(6)

(7)

同样的Ki可以表示为：

Ki=

(8)

显然，将所有单元网格的转换矩阵拼接得到的K和依据式(1)将所有网格点上的权值整合得到的拉普拉斯转换矩阵是相同的。

1.2 拉普拉斯变换下重构信号的稀疏性

本文利用核状流与层流模型的仿真数据验证重构信号在K矩阵下的稀疏性。如图2所示，上方两幅图形为原始信号绘制的图像，灰色部分网格点数值为0，非0点数分别为909与1 085。下方两幅图形为经过拉普拉斯变换后信号的图像，浅灰色部分网格点数值为0，非0点数分别为220与112。

图2 拉普拉斯变换前后的核状流与层流图像

2 ECT图像重构问题求解

2.1 加入拉普拉斯1范数约束的ECT逆问题

本文在ECT图像重构问题中加入重构信号在拉普拉斯变化后的1范数约束得到式(9)，其中K为上述Laplace变化矩阵。

(9)

2.2 Split-Bregman算法求解

为了方便表达，根据式(10)、(11)和(12)定义如下函数关系：

φ(δε)=Kδε

(10)

(11)

J(δε,d)=‖d‖1+H(δε)

(12)

式中：d为过程变量，不代表其他含义。

最终，式(9)转化为式(13)：

(13)

式中：λ为调节参数。

Split-Bregman算法求解流程[11]如下。

1)初始化：k=0,δε0=0,b0=0。

2)执行下式：

bk+1=bk+(φ(δεk+1)-dk+1)

k=k+1

3)输出：δε=δεk。

其中：k为迭代次数；δε0为δε的迭代初值；b0为b的迭代初值；δεk+1为δε的(k+1)次迭代值；dk为d的k次迭代值；bk为b的k次迭代值；t为停止求解的阈值，大约为1×10-5；λ为10～20。

求解过程使用MATLAB和CVX工具包。

3 静态实验验证

3.1 实验设计

3.1.1实验设备

实验所用ECT采集系统由3部分组成：作为采集激励模块的火龙果(red pitaya)、为选通芯片提供I/O信号的STM32以及电容电压转换模块(如图3所示)。末端传感器由柔性电极制成。测量管为外径50 mm、厚度2 mm的亚克力管，铜电极片覆盖面积达85%，长度等于管径，如图4所示。

图3 实验设备原理图

图4 检测设备与末端传感器

3.1.2仿真实验设计

本文使用相对介电常数值为3～4的面粉和聚氨酯棒，模拟如图5所示的典型的4种流型：分层流、环形流、核状流、双核流。采用Landweber迭代算法、Tikhonov正则化算法、FOT正则化算法和拉普拉斯变换后的1范数约束算法进行图像重建。采用式(14)中定义的相关系数CC[12]和式(15)中定义的均方根误差RMSE[13]定量评价不同图像重构算法的重建质量。

图5 4种流型实验模型图

(14)

(15)

3.2 实验结果分析

图像重建时，做如下归一化处理：

(16)

(17)

由于ECT传感器只有8个电极，因此重建图像的CC较小。然而，这并不妨碍验证所提算法的有效性。从理论上讲，本文算法能够有效地减少图像中的伪影，突出图像边缘。

从图6中可以看出，用拉普拉斯变换后的1范数约束算法重建的图像CC得到提高，伪影更少，边缘更清晰。为了便于比较重构图像的质量，将评价指标CC和RMSE表示为如图7所示的柱状图，可以看出拉普拉斯变换后的1范数约束算法的重构图像质量得到了改善。

4 结束语

本文提出了用于ECT图像重建的拉普拉斯1范数约束算法，利用拉普拉斯变换下信号的稀疏表示作为1范数约束项，在ECT图像重构问题中有着较好的表现。该算法类似于总变分算法，并且在三角单元网格中具有更简便的数学表达，使用Spilt-Bregman迭代求解有着较好的结果，可进一步应用于其他图像处理问题中。