基于联合代数重建的电容层析成像图像重建

宋亚杰， 张立峰， 朱炎峰

(华北电力大学 控制与计算机工程学院，河北 保定 071003)

0 引言

电容层析成像(Electrical Capacitance Tomography，ECT)技术是一种新型的多相流检测技术[1]。其原理是通过安装在管道外侧的电极阵列采集数据，对管道内部介电常数分布进行实时可视化测量。ECT系统由于其非侵入、响应速度快、结构简单、无辐射的特点，已在多相流领域逐渐得到应用[2，3]。

ECT图像重建算法是一个逆问题求解的过程，但是逆问题具有非线性和病态性[4]，而非线性和病态性都将直接影响图像精度。目前常用的ECT重建算法可分为非迭代算法和迭代算法，非迭代算法如线性反投影算法(Linear Back Projection,LBP)虽然成像速度快，但很难取得较高精度的重建图像；迭代算法成像精度明显优于非迭代算法，如Landweber算法、Tikhonov正则化算法等，但一般收敛速度较慢[5]。

代数重建技术(Algebraic Reconstruction Technique，ART)是投影重建的典型方式之一，文献[6]将ART算法应用在ECT图像重建上，取得了较好的成像效果，但ART算法计算量大，重建速度较慢。联合代数重建技术(Simultaneous Algebraic Reconstruction Technique，SART)是在ART算法基础上改进而来，可有效克服逆问题的病态性和不适定性，不仅加快了迭代速度，重建图像的质量也得到进一步提高，且具有一定的抗噪能力[7]。本文将SART算法应用于ECT图像重建，采用16电极ECT系统，针对4种常见流型，分别使用ART算法和SART算法进行仿真并重建图像，同时用成像效果、迭代次数、迭代时间和相对误差等评价指标进行对比分析。最后，使用天津大学研制的数字化TJUET测量系统，模拟油/气两相流进行了静态实验，基于两种算法得到重建图像，并对比了分析测试结果。

1 理论基础

1.1 ECT数学模型

ECT系统测量过程中，依次选择其中一个电极作为激励电极，其它电极作为测量电极用以获取所有不同电极对之间的电容数据。由此可获得独立测量电容数据个数为：

M=K(K-1)/2

(1)

式中：K为电极数。本文使用16电极传感器，因此将会得到120个独立测量电容数据。

ECT系统中，被测圆形管道内部电场分布方程为：

·[ε(x,y)]φ(x,y)=0

(2)

式中：ε(x,y)为相对介电常数；φ(x,y)为电势；φ(x,y)为电场强度。对式(1)进行离散化，即可得到ECT系统近似线性模型：

C=SG+e

(3)

式中：C为N×1维归一化后的电容测量值；S为M×N维归一化后的灵敏度系数矩阵；G为需要N×1维归一化后的原始图像灰度值向量；e是噪声向量；N为剖分的网格数。独立测量数据个数往往小于网格数，即M

1.2 代数重建算法(ART)

代数重建算法[8]广泛应用于图像重建领域。该算法要求待重建的截面图像f(x,y)分成N个离散的像素方格，而每个像素方格内部的f(x,y)都是常数，记作fi。然后，根据成像的物理过程和相应的数学模型，得到重建图像和投影数据之间关系的线性方程组[9]。

(4)

式中：pi为第i号射线的测量数据；wij为反映像元fi对投影pi贡献的权重因子；M为投影射线总数；N是像素方格个数；pi可视为归一化的测量电容值；wij可认为是灵敏度矩阵单元；fi便是待重建的像素单元灰度值。观察可知，式(4)可近似表示为ECT系统的线性模型方程(3)，因此可用于ECT系统的图像重建。图像重建的任务即是由w、p求取f。

ART算法的具体步骤为：

步骤1：设定图像初始值。初始值的设定将会直接影响到图像重建迭代次数，一般情况下可设初始值如式(5)；为提高重建速率，本文则以LBP算法得到的解作为初始值。

(5)

步骤2：根据下式计算图像的修正值

(6)

式中：k为迭代次数；1≤i≤M；1≤j≤N；λ为松弛因子(0<λ<2)。为更好地抑制伪影，对于像元fi< 0，则使fi= 0。当所有的投影数据都被使用完，便完成了一次完整的迭代。

步骤3：反复进行第2步中的运算，对图像进行修正。

步骤4：检验是否满足收敛准则。常用的收敛标准是：

(7)

当(U(k)N)1/2<ε(0<ε<1)时，停止迭代，也可设置迭代次数作为收敛停止条件。

1.3 联合代数重建算法(SART)

作为ART算法的改进技术，SART算法仍然以式(4)为基础进行改进。由上文可知，式(4)可近似表示为ECT系统的线性模型，因此SART算法克服了逆问题求解的非线性问题；在SART算法中，每次迭代对像素的校正值并非只与一个投影数据有关，而是所有投影数据校正值的平均值，一些随机误差被平均掉了，这也就使得SART算法能够更好地抑制带状伪影，得到比ART算法更加平滑的重建图像[10～12]。

SART算法的具体步骤为：

步骤1：设置图像初始值(见ART算法第1步)。

步骤2：SART算法子迭代过程可分解为以下两步：

(8)

(9)

步骤3：反复进行第2步中的运算，对图像进行修正。

步骤4：检验是否满足收敛准则(见ART算法第4步)。

2 实验结果及分析

2.1 仿真实验

仿真实验平台基于Windows8 Intel Core i7 2.60 GHz操作系统，内存为4 GB RAM；仿真软件采用MATLAB R2009 a和COMSOL3.5 a有限元软件。网格剖分采用自适应三角形网格剖分，图像重建采用正方形网格剖分，将成像区域分为812个单元。采用16电极ECT系统，电极均匀排列在圆形管道外侧，屏蔽罩和径向屏蔽电极接地，用于屏蔽外界电磁干扰对测量信号的影响。其中电极张角θ=12°，管道内径为80 mm，如图1所示。

图1 16电极ECT传感器结构图

仿真实验以油/气相流为研究对象，设置了如图2所示4种常见流型。其中油和气的相对介电常数分别设为3和1。使用COMSOL3.5 a有限元软件获得仿真测量电容值数据，并在此基础上加入±2%的高斯白噪声；利用MATLAB软件编写ART算法和SART算法完成ECT图像重建。

图像重建算法参数设置：为增加ART算法和SART算法重建效果的对比有效性，选取线性反投影算法(LBP)的重建结果作为重建算法的图像初始值；迭代次数设置为150次；固定松弛因子λ=1。重建图像如图2所示。

图2 仿真重建图像

为定量评价重建图像质量，采用相对误差(IE)和相关系数(CC)作为评价标准。

(10)

(11)

表2 重建图像相关系数

由表1和表2可知，对所选4种常见流型，ART算法和SART算法的相对误差和相关系数较为接近；结合图2的重建图像，不难发现SART算法能对不同物体有效地区分，且边缘信息保真度高，具有较少的伪影。

流型1和流型2中，ART算法成像的物体略微变形，而SART算法成像的物体更接近原流型；此外可以看到，在ART算法重建图像中，管壁部分存在伪影，而SART算法图像中几乎不存在；对于流型3的层流，两者成像效果都不理想，尤其是靠近管壁的部分失真严重，相比较而言，SART算法的成像效果略好一些；而对于流型4的环流，如果不考虑伪影的存在，ART算法优于SART算法，成像更接近原流型。

表3为采用ART和SART的图像重建耗时对比。

表3 图像重建耗时 s

另外，由表3可以看出，对于4种常见流型，SART算法的成像时间最短。尤其对流型4，成像速度比ART算法快了接近一倍。

2.2 静态实验

采用天津大学研制的数字化TJUET测量系统进行了静态实验，数字化TJUET系统如图3所示。

图3 数字化TJUET系统

两相流模拟实验中采用传感器模型为内径 60 mm的圆形有机玻璃管道，16电极电容敏感阵列。流型1’、流型2’和流型3’使用的是直径分别为5 mm、10 mm、20 mm的塑料棒，而流型4’使用的是内径30 mm、厚度2 mm的空心PC塑料棒，空场及满场标定介质分别为空气及塑料颗粒，管内剖分812个单元。实验结果如图4所示。

图4 静态实验重建图像

由图4可知，对于流型1’和流型2’，ART算法成像的物体变形严重，且与管壁粘连，而SART算法成像边缘信息保留较好，图像更接近真实分布；对于流型3’的层流，两者成像效果仍不理想，高介电物质呈离散分布；而对于流型4’的环流，两种算法都取得了较好成像效果，相比较下，ART算法成像物体大小更接近原流型。

3 结论

本文针对ECT图像重建中逆问题的非线性和病态性，将基于ART算法的SART算法应用于ECT图像重建。通过仿真和实验测试显示，ART算法和SART算法皆可有效地对物体进行重建，而SART算法的成像物体边缘信息保真度高，且能抑制大部分伪影，具有一定的抗噪性能；另外，在成像效果相差不大的情况下，SART算法的成像速度要优于ART算法。该算法为解决ECT成像中的逆问题提供了一条新的途径。