利用场量旋转变换的电容层析成像灵敏度系数简易计算方法

唐凯豪,胡红利,李林,王小鑫

(1.西安交通大学电气工程学院,710049,西安;2.西安石油大学电子工程学院,710065,西安)

电容层析成像(Electrical Capacitance Tomography,ECT)技术是一种过程层析成像技术,通过测量安装在管道外壁的电容传感器阵列的电容值来反演管道内的介质分布。该技术具有传感器非侵入、结构简单、造价低廉等优点,因而被广泛应用于需要对多相流参数监测与可视化的工业过程中,如电力工业、石油工业、化工等行业。

ECT技术包括正问题与逆问题:正问题是根据已知的介质分布求解各电极对之间的互电容;逆问题是根据测得的互电容值反演管道内介质分布,即图像重构问题。基于灵敏度理论的图像重构是ECT逆问题的一种较广应用,因此,传感器灵敏度系数的计算成为影响该逆问题求解质量的关键因素之一。

目前,通过求解场域内电势分布来求解灵敏度系数是一种较常见的方法[1-2],即场量提取法。该方法的思想可追溯到Geselowitz于1970年提出的阻抗灵敏度的计算方法[3],实质上属于微扰法[4]的思想。另外,灵敏度系数也可通过灵敏度系数的定义直接计算,即扰动法。与场量提取法相比,扰动法的计算过程更复杂、需要进行的有限元(或其他数值方法)求解次数也更多。虽然基于场量提取法的灵敏度系数求解仅需要进行电极个次数的有限元(或其他数值方法)求解,但在计算机性能不佳、管道结构较复杂、需要精确建模等情况时,算法的时间复杂度仍然较高,进行灵敏度系数矩阵计算也有一定的烦琐性。尤其是对于电极数目较多的传感器的研究(相比12电极传感器),例如:Mohamad等人研究的16电极ECT棕榈油监测设备[5],Peng等人研究的4、8、12、16、20、24和32电极的传感器图像重构质量[6],Yang等人研究的24电极ECT传感器及其电极组合策略[7],Liu等人研究的旋转ECT传感器[8],范优飞等提出的电极组合可配置的48电极ECT系统[9]等。这些研究的灵敏度系数矩阵的求解复杂度较高,工作量较大。

考虑到ECT传感器电极长度通常不小于管道直径的一倍[10],因此在敏感区域中部的电势函数可近似认为与轴向无关,从而三维场问题可以简化为二维场问题。所以,本文提出了基于电势函数旋转变换的二维灵敏度系数求解方法(以下称为旋转变换法)。该方法是基于电势的灵敏度系数求解法的扩展,只需进行一次数值方法计算即可得到全部灵敏度系数,较大程度简化了灵敏度系数的求解过程。

首先,分析了ECT传感器的对称结构,得到了ECT不同激励模式下电势函数之间的旋转对称关系,基于该关系,运用矢量分析理论得到了电场强度之间的旋转对称性,并根据该对称关系给出了只运用一种激励模式下的电势函数计算灵敏度系数的方法。其次,对比了运用该方法得到的灵敏度系数矩阵和运用传统方法即场量提取法得到的灵敏度系数矩阵之间的均方误差,并且用两种灵敏度系数矩阵进行了图像重构实验(仿真实验和实测实验)。最后,从重构图像与原物场图像的相关系数和相对均方误差两方面对比了两种灵敏度系数矩阵的重构图像,验证了提出的方法的可靠性。

1 电容层析成像灵敏度系数与静电场的关系

1.1 电容层析成像的物理模型

12电极ECT传感器的剖面图如图1所示。

图1 ECT传感器剖面图

根据ECT的测量原理,当i号电极施加激励、j号电极检测时(称为i激励模式),其余电极均应接地或与地等电势[11]。因此,容易得到当管道内介质不带电时,电极数为N的ECT传感器(不包括屏蔽罩)的物理模型,可描述为电势函数的狄里克莱问题,其公式为

(1)

式中:φi是电极i作为激励电极时,传感器敏感区域内的电势函数;Ω是传感器的敏感区域;∂Ωi是第i个电极的表面,∂Ωj是第j个电极的表面;VE是传感器的激励电压,该激励模式下的电容测量值为

(2)

其中,Ci,j是电极i和电极j之间的互电容,qj是电极j上的感应电荷,ds是第二类面积分的积分元。

通过数值方法求解式(1)和式(2),即可得到ECT正问题所需的全部量值。

1.2 电容层析成像传感器敏感场的旋转对称性

如图1所示,注意到ECT传感器的电极等大小、等间距排布,电压激励相同,因而当管道内介质线性、均匀、各向同性时,i、j激励模式的电势函数φi(x,y)和φj(x,y)之间只是绕z轴旋转了一个角度的关系,即不同激励模式情况下的正问题不需要重复求解。设电势函数φi(x,y)旋转θ后的函数表达式为

φj(x′,y′)=φi(x,y)

(3)

(4)

其中

(5)

同理,也可得到

φj(x,y)=φj(x)=φi(A-1x)=φi(ATx)

(6)

由于A是正交阵,A-1=AT,因此式(6)中的逆阵均用转置表示。式(5)和式(6)揭示了不同激励模式下电势函数的旋转对称性。

根据电磁场理论,电场E和电势函数φ的关系为

E=-φ

(7)

将式(4)带入式(7),并根据矩阵函数微分理论[12]可得

Ej(Ax)=-xφj(Ax)=-AAxφj(Ax)=

-Axφi(x)=AEi(x)

(8)

式中:Ei表示当电极i为激励电极时,传感器敏感区域内的电场;x表示作用于x的哈密顿算子。上述过程为“正推”过程,即已知φi(x)求其对应的φj(Ax)的表达式。同理,可得“逆推”表达式为

Ej(x)=-xφj(x)=-ATATxφi(ATx)=

ATEi(ATx)

(9)

式(8)和式(9)即为场强旋转对称性的表达式。在线性、均匀、各向同性介质中,电通密度有同样的形式。

2 基于单一电势函数的灵敏度系数求解

2.1 理论推导

ECT正问题求解[13]的灵敏度系数为si,j(k),表示当第i个电极施加激励、第j个电极作为检测电极时,场域内第k个单元的介质对互电容Ci,j的贡献,其求解公式为

(10)

(11)

当成像区域均匀剖分,且单元数足够多时,式(11)中的积分可表示为坐标点函数值与区域面积的乘积。由于逆问题求解之前灵敏度系数需要归一化处理,因此在均匀剖分时,区域面积与激励电压构成的系数可用1代替,从而式(11)可被进一步简化为

si,j(xk,yk)=-φi(xk,yk)φj(xk,yk)

(12)

式中(xk,yk)为第k个单元的中心坐标。将式(9)代入式(12)得到

si,j(xk,yk)=-φi(xk)φj(xk)=

-Ei(xk)[ATEi(ATxk)]

(13)

至此,得到了新的二维ECT灵敏度系数表达式,即本文提出的旋转变换法表达式。

利用商业数值模拟软件很容易得到式(13)中的电势与场强。以COMSOL Multiphysics软件为例,在使用该软件的静电接口进行有限元计算时,电势分布与电场强度各方向分量为默认求解变量,无需用户再进行数据后处理。

2.2 逆问题剖分算法

图像重构时,根据式(5),若被测截面采用直角坐标进行剖分,旋转后的坐标有可能会落在计算点之外,即式(13)中的φi(ATxk)可能不存在于φi的函数值列表中,此时需要进行二元函数差值运算,这将带来不便。

为此,本文提出一种均匀极坐标剖分的逆问题剖分规则。极坐标剖分时,在圆周r=rk上取格点定为单元中心,每圈的点数随圆周半径增加而变化。该规则能保证按照式(13)计算灵敏度系数时,旋转后的函数值仍位于原函数值列表中。规则描述如下。

位于第k圈的任一单元和该单元中的唯一格点(rk,θ)满足:

条件1,所有单元面积都为s0=Ap/np,式中np为逆问题剖分预定单元数,每圈格点等角分布;

条件5,每圈直径Rk和单元弧长Lk=2πRk/nk与第一圈单元的上述指标的误差在阈值tg内,即(Rk-R1)/R1≤tg,(Lk-L1)/L1≤tg;

条件6,最终单元剖分数n与np的差在阈值tu内,即(np-n)/np≤tu,且最后一圈半径Rn≤Rp,Rp为被测截面的半径。

满足上述条件1～6的逆问题剖分算法的流程图如图2所示。

图2 逆问题均匀极坐标剖分算法

本文的研究对象是直径50 mm、壁厚5 mm的石英玻璃管道。预设剖分单元数为3 204,并设置tg=0.5,tu=0.05,利用图2所示算法可将成像区域剖分为面积相等的3 180个曲边四边形单元和1个位于圆心的圆形单元,如图3所示。当设置该阈值时,能保证最终剖分单元数与预期剖分单元数的偏差在5%以内(本文中为0.75%),每个单元弧长与L1最大差距在50%以内(本文中为7.99%)。可见,该算法能保证逆问题单元剖分的均匀性。

图3 成像单元极坐标格点

将图3所示的格点坐标按列排布,输入COMSOL Multiphysics软件即可导出这些格点上的电势、电场强度值。

3 仿真与实验验证

3.1 仿真模型建立

采用COMSOL Multiphysics软件作为有限元仿真工具,对12电极ECT传感器进行正问题求解。建立管道直径传感器仿真模型,如图4所示。

图4 ECT传感器仿真模型截面图

管道材料为石英玻璃,传感器模型的仿真参数如表1所示。

将脆弱性mvul的变换周期记为interval，根据定义3，interval=min{Δti}.将该脆弱性的变换空间记为W，根据笛卡尔积定义，W=P1×P2×…×Pn，将W大小记为|W|，则|W|=|P1|·|P2|·…·|Pn|，其中|Pi|(1≤i≤n)为第i个属性的值域空间大小.

表1 传感器模型的仿真参数

仿真时,激励电极的边界条件为电势3.3 V,其余电极的边界条件均为接地(即零电势)。

需要说明的是,传统灵敏度系数计算方法需要依次设置12个电极为激励电极,即进行12次有限元计算,而本文方法只需设置1号电极作为激励电极,即进行1次有限元计算。

有限元计算时采用三角形自适应网格剖分,网格剖分情况如图5所示。

图5 有限元计算网格剖分

3.2 灵敏度系数计算结果对比

限于篇幅,选取S1,2、S1,7、S1,10三个灵敏度作为对象进行对比。采用场量提取法和式(13)所示的旋转变换法计算灵敏度系数,上述三个灵敏度的计算结果分别如图6和图7所示。由于灵敏度是二维曲面,可读性较差,因此以灰度图的方式呈现,图中越白的地方表示灵敏度系数越大,越黑的地方表示灵敏度系数越小。从图中可见,两种方法的灵敏度计算结果相似。

(a)S1,2 (b)S1,7 (c)S1,10图6 场量提取法灵敏度计算结果

(a)S1,2 (b)S1,7 (c)S1,10图7 旋转变换法灵敏度计算结果

灵敏度相对均方误差Es的计算公式为

Es=‖SF-SR‖2/‖SF‖2

(14)

式中:SF是场量提取法计算得到的灵敏度;SR是旋转变换法计算得到的灵敏度。上述两种方法计算结果的相对均方误差分别为0.013 2、0.020 6、0.005 3。可见,本文所述方法与传统方法的计算结果高度吻合。

3.3 图像重构实测实验

采用12电极ECT传感器进行重构图像实测实验。电容测量电路为文献[15]中所述的交流法微小电容测量电路。实验时,将不同直径的有机玻璃棒摆放在管道内不同位置或在管道内注水来模拟不同流形时管道内物场分布。设置三种物场分布代表三种流型进行实验,分别为:一根粗介质棒(核心流,直径20 mm)、一根偏心介质棒(偏心流,直径20 mm)、三根介质棒(每根直径16 mm)。这三种物场分布的二值灰度图像如图8所示。

(a)中心介质棒 (b)偏心介质棒 (c)对称三介质棒图8 物场分布设置

采取广泛使用的OIOR-Landweber(Offline Iteration Online Reconstruction Landweber)[16]法对测得的电容值进行图像重构,用两种方法计算的灵敏度系数矩阵重构的图像如图9所示,迭代次数为3 000次。

(a) 场量提取法重构(b) 场量提取法重构(c) 场量提取法重构的中心介质棒 的偏心介质棒 的对称三介质棒

(d) 旋转变换法重构(e) 旋转变换法重构(f) 旋转变换法重构的中心介质棒 的偏心介质棒 的对称三介质棒图9 重构图像

规定在图8所示的物场分布图中,介质存在区域灰度为1,其余区域灰度用于对比两种灵敏度系数矩阵重构的图像。为对比重构图像中非零最小元素,采用较常用的相关系数与相对均方误差[17]两项指标作为评价标准,结果如表2所示。

表2 图像重构情况对比

可见,基于上述两种灵敏度系数矩阵的图像重构结果吻合度高,证明了前述理论分析的正确性。

4 结 论

(1)基于矩阵函数理论,提出了一种基于矩阵函数旋转变换的ECT二维灵敏度系数求解方法,即旋转变换法,并到了该方法的数学表达式。运用该方法,仅进行一次数值计算即可获得ECT传感器的全部灵敏度系数。该方法与传统的场量提取法的灵敏度求解结果具有高度一致性,二者计算结果的相对偏差在2.06%以下,证明旋转变换法可以替代传统的场量提取法,以减少ECT正问题求解的工作量。

(2)针对所提出的旋转变换法,设计了一种均匀的极坐标成像区域剖分算法。该算法可保证对电势(或电场)函数做旋转变换的过程中不产生新的插值节点,从而进一步减少了灵敏度矩阵的计算量。

(3)在三种物场分布情况下进行了实验,分别运用传统的场量提取法和旋转变化法计算得到的灵敏度系数矩阵,对相同的ECT测量数据进行了图像重构。实验结果显示,两种方法的图像相关系数和相对误差高度吻合,表明在实际工程应用中,可以使用旋转变化法进行ECT灵敏度系数求解。

(4)ECT/ERT双模态系统是电学层析成像技术的研究热点之一,考虑到ERT的灵敏度系数表达与ECT灵敏度系数表达有相同的形式[10],因此旋转变化法对于双模态问题的灵敏度求解同样有效。