利用場量旋轉變換的電容層析成像靈敏度系數簡易計算方法

唐凱豪,胡紅利,李林,王小鑫

(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],因此旋轉變化法對于雙模態問題的靈敏度求解同樣有效。