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

(1.西安交通大學(xué)電氣工程學(xué)院,710049,西安;2.西安石油大學(xué)電子工程學(xué)院,710065,西安)

電容層析成像(Electrical Capacitance Tomography,ECT)技術(shù)是一種過程層析成像技術(shù),通過測量安裝在管道外壁的電容傳感器陣列的電容值來反演管道內(nèi)的介質(zhì)分布。該技術(shù)具有傳感器非侵入、結(jié)構(gòu)簡單、造價低廉等優(yōu)點,因而被廣泛應(yīng)用于需要對多相流參數(shù)監(jiān)測與可視化的工業(yè)過程中,如電力工業(yè)、石油工業(yè)、化工等行業(yè)。

ECT技術(shù)包括正問題與逆問題:正問題是根據(jù)已知的介質(zhì)分布求解各電極對之間的互電容;逆問題是根據(jù)測得的互電容值反演管道內(nèi)介質(zhì)分布,即圖像重構(gòu)問題?；陟`敏度理論的圖像重構(gòu)是ECT逆問題的一種較廣應(yīng)用,因此,傳感器靈敏度系數(shù)的計算成為影響該逆問題求解質(zhì)量的關(guān)鍵因素之一。

目前,通過求解場域內(nèi)電勢分布來求解靈敏度系數(shù)是一種較常見的方法[1-2],即場量提取法。該方法的思想可追溯到Geselowitz于1970年提出的阻抗靈敏度的計算方法[3],實質(zhì)上屬于微擾法[4]的思想。另外,靈敏度系數(shù)也可通過靈敏度系數(shù)的定義直接計算,即擾動法。與場量提取法相比,擾動法的計算過程更復(fù)雜、需要進行的有限元(或其他數(shù)值方法)求解次數(shù)也更多。雖然基于場量提取法的靈敏度系數(shù)求解僅需要進行電極個次數(shù)的有限元(或其他數(shù)值方法)求解,但在計算機性能不佳、管道結(jié)構(gòu)較復(fù)雜、需要精確建模等情況時,算法的時間復(fù)雜度仍然較高,進行靈敏度系數(shù)矩陣計算也有一定的煩瑣性。尤其是對于電極數(shù)目較多的傳感器的研究(相比12電極傳感器),例如:Mohamad等人研究的16電極ECT棕櫚油監(jiān)測設(shè)備[5],Peng等人研究的4、8、12、16、20、24和32電極的傳感器圖像重構(gòu)質(zhì)量[6],Yang等人研究的24電極ECT傳感器及其電極組合策略[7],Liu等人研究的旋轉(zhuǎn)ECT傳感器[8],范優(yōu)飛等提出的電極組合可配置的48電極ECT系統(tǒng)[9]等。這些研究的靈敏度系數(shù)矩陣的求解復(fù)雜度較高,工作量較大。

考慮到ECT傳感器電極長度通常不小于管道直徑的一倍[10],因此在敏感區(qū)域中部的電勢函數(shù)可近似認為與軸向無關(guān),從而三維場問題可以簡化為二維場問題。所以,本文提出了基于電勢函數(shù)旋轉(zhuǎn)變換的二維靈敏度系數(shù)求解方法(以下稱為旋轉(zhuǎn)變換法)。該方法是基于電勢的靈敏度系數(shù)求解法的擴展,只需進行一次數(shù)值方法計算即可得到全部靈敏度系數(shù),較大程度簡化了靈敏度系數(shù)的求解過程。

首先,分析了ECT傳感器的對稱結(jié)構(gòu),得到了ECT不同激勵模式下電勢函數(shù)之間的旋轉(zhuǎn)對稱關(guān)系,基于該關(guān)系,運用矢量分析理論得到了電場強度之間的旋轉(zhuǎn)對稱性,并根據(jù)該對稱關(guān)系給出了只運用一種激勵模式下的電勢函數(shù)計算靈敏度系數(shù)的方法。其次,對比了運用該方法得到的靈敏度系數(shù)矩陣和運用傳統(tǒng)方法即場量提取法得到的靈敏度系數(shù)矩陣之間的均方誤差,并且用兩種靈敏度系數(shù)矩陣進行了圖像重構(gòu)實驗(仿真實驗和實測實驗)。最后,從重構(gòu)圖像與原物場圖像的相關(guān)系數(shù)和相對均方誤差兩方面對比了兩種靈敏度系數(shù)矩陣的重構(gòu)圖像,驗證了提出的方法的可靠性。

1 電容層析成像靈敏度系數(shù)與靜電場的關(guān)系

1.1 電容層析成像的物理模型

12電極ECT傳感器的剖面圖如圖1所示。

圖1 ECT傳感器剖面圖

根據(jù)ECT的測量原理,當i號電極施加激勵、j號電極檢測時(稱為i激勵模式),其余電極均應(yīng)接地或與地等電勢[11]。因此,容易得到當管道內(nèi)介質(zhì)不帶電時,電極數(shù)為N的ECT傳感器(不包括屏蔽罩)的物理模型,可描述為電勢函數(shù)的狄里克萊問題,其公式為

(1)

式中:φi是電極i作為激勵電極時,傳感器敏感區(qū)域內(nèi)的電勢函數(shù);Ω是傳感器的敏感區(qū)域;?Ωi是第i個電極的表面,?Ωj是第j個電極的表面;VE是傳感器的激勵電壓,該激勵模式下的電容測量值為

(2)

其中,Ci,j是電極i和電極j之間的互電容,qj是電極j上的感應(yīng)電荷,ds是第二類面積分的積分元。

通過數(shù)值方法求解式(1)和式(2),即可得到ECT正問題所需的全部量值。

1.2 電容層析成像傳感器敏感場的旋轉(zhuǎn)對稱性

如圖1所示,注意到ECT傳感器的電極等大小、等間距排布,電壓激勵相同,因而當管道內(nèi)介質(zhì)線性、均勻、各向同性時,i、j激勵模式的電勢函數(shù)φi(x,y)和φj(x,y)之間只是繞z軸旋轉(zhuǎn)了一個角度的關(guān)系,即不同激勵模式情況下的正問題不需要重復(fù)求解。設(shè)電勢函數(shù)φi(x,y)旋轉(zhuǎn)θ后的函數(shù)表達式為

φ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)中的逆陣均用轉(zhuǎn)置表示。式(5)和式(6)揭示了不同激勵模式下電勢函數(shù)的旋轉(zhuǎn)對稱性。

根據(jù)電磁場理論,電場E和電勢函數(shù)φ的關(guān)系為

E=-φ

(7)

將式(4)帶入式(7),并根據(jù)矩陣函數(shù)微分理論[12]可得

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

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

(8)

式中:Ei表示當電極i為激勵電極時,傳感器敏感區(qū)域內(nèi)的電場;x表示作用于x的哈密頓算子。上述過程為“正推”過程,即已知φi(x)求其對應(yīng)的φj(Ax)的表達式。同理,可得“逆推”表達式為

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

ATEi(ATx)

(9)

式(8)和式(9)即為場強旋轉(zhuǎn)對稱性的表達式。在線性、均勻、各向同性介質(zhì)中,電通密度有同樣的形式。

2 基于單一電勢函數(shù)的靈敏度系數(shù)求解

2.1 理論推導(dǎo)

ECT正問題求解[13]的靈敏度系數(shù)為si,j(k),表示當?shù)趇個電極施加激勵、第j個電極作為檢測電極時,場域內(nèi)第k個單元的介質(zhì)對互電容Ci,j的貢獻,其求解公式為

(10)

(11)

當成像區(qū)域均勻剖分,且單元數(shù)足夠多時,式(11)中的積分可表示為坐標點函數(shù)值與區(qū)域面積的乘積。由于逆問題求解之前靈敏度系數(shù)需要歸一化處理,因此在均勻剖分時,區(qū)域面積與激勵電壓構(gòu)成的系數(shù)可用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靈敏度系數(shù)表達式,即本文提出的旋轉(zhuǎn)變換法表達式。

利用商業(yè)數(shù)值模擬軟件很容易得到式(13)中的電勢與場強。以COMSOL Multiphysics軟件為例,在使用該軟件的靜電接口進行有限元計算時,電勢分布與電場強度各方向分量為默認求解變量,無需用戶再進行數(shù)據(jù)后處理。

2.2 逆問題剖分算法

圖像重構(gòu)時,根據(jù)式(5),若被測截面采用直角坐標進行剖分,旋轉(zhuǎn)后的坐標有可能會落在計算點之外,即式(13)中的φi(ATxk)可能不存在于φi的函數(shù)值列表中,此時需要進行二元函數(shù)差值運算,這將帶來不便。

為此,本文提出一種均勻極坐標剖分的逆問題剖分規(guī)則。極坐標剖分時,在圓周r=rk上取格點定為單元中心,每圈的點數(shù)隨圓周半徑增加而變化。該規(guī)則能保證按照式(13)計算靈敏度系數(shù)時,旋轉(zhuǎn)后的函數(shù)值仍位于原函數(shù)值列表中。規(guī)則描述如下。

位于第k圈的任一單元和該單元中的唯一格點(rk,θ)滿足:

條件1,所有單元面積都為s0=Ap/np,式中np為逆問題剖分預(yù)定單元數(shù),每圈格點等角分布;

條件5,每圈直徑Rk和單元弧長Lk=2πRk/nk與第一圈單元的上述指標的誤差在閾值tg內(nèi),即(Rk-R1)/R1≤tg,(Lk-L1)/L1≤tg;

條件6,最終單元剖分數(shù)n與np的差在閾值tu內(nèi),即(np-n)/np≤tu,且最后一圈半徑Rn≤Rp,Rp為被測截面的半徑。

滿足上述條件1～6的逆問題剖分算法的流程圖如圖2所示。

圖2 逆問題均勻極坐標剖分算法

本文的研究對象是直徑50 mm、壁厚5 mm的石英玻璃管道。預(yù)設(shè)剖分單元數(shù)為3 204,并設(shè)置tg=0.5,tu=0.05,利用圖2所示算法可將成像區(qū)域剖分為面積相等的3 180個曲邊四邊形單元和1個位于圓心的圓形單元,如圖3所示。當設(shè)置該閾值時,能保證最終剖分單元數(shù)與預(yù)期剖分單元數(shù)的偏差在5%以內(nèi)(本文中為0.75%),每個單元弧長與L1最大差距在50%以內(nèi)(本文中為7.99%)?？梢?該算法能保證逆問題單元剖分的均勻性。

圖3 成像單元極坐標格點

將圖3所示的格點坐標按列排布,輸入COMSOL Multiphysics軟件即可導(dǎo)出這些格點上的電勢、電場強度值。

3 仿真與實驗驗證

3.1 仿真模型建立

采用COMSOL Multiphysics軟件作為有限元仿真工具,對12電極ECT傳感器進行正問題求解。建立管道直徑傳感器仿真模型,如圖4所示。

圖4 ECT傳感器仿真模型截面圖

管道材料為石英玻璃,傳感器模型的仿真參數(shù)如表1所示。

將脆弱性mvul的變換周期記為interval，根據(jù)定義3，interval=min{Δti}.將該脆弱性的變換空間記為W，根據(jù)笛卡爾積定義，W=P1×P2×…×Pn，將W大小記為|W|，則|W|=|P1|·|P2|·…·|Pn|，其中|Pi|(1≤i≤n)為第i個屬性的值域空間大小.

表1 傳感器模型的仿真參數(shù)

仿真時,激勵電極的邊界條件為電勢3.3 V,其余電極的邊界條件均為接地(即零電勢)。

需要說明的是,傳統(tǒng)靈敏度系數(shù)計算方法需要依次設(shè)置12個電極為激勵電極,即進行12次有限元計算,而本文方法只需設(shè)置1號電極作為激勵電極,即進行1次有限元計算。

有限元計算時采用三角形自適應(yīng)網(wǎng)格剖分,網(wǎng)格剖分情況如圖5所示。

圖5 有限元計算網(wǎng)格剖分

3.2 靈敏度系數(shù)計算結(jié)果對比

限于篇幅,選取S1,2、S1,7、S1,10三個靈敏度作為對象進行對比。采用場量提取法和式(13)所示的旋轉(zhuǎn)變換法計算靈敏度系數(shù),上述三個靈敏度的計算結(jié)果分別如圖6和圖7所示。由于靈敏度是二維曲面,可讀性較差,因此以灰度圖的方式呈現(xiàn),圖中越白的地方表示靈敏度系數(shù)越大,越黑的地方表示靈敏度系數(shù)越小。從圖中可見,兩種方法的靈敏度計算結(jié)果相似。

(a)S1,2 (b)S1,7 (c)S1,10圖6 場量提取法靈敏度計算結(jié)果

(a)S1,2 (b)S1,7 (c)S1,10圖7 旋轉(zhuǎn)變換法靈敏度計算結(jié)果

靈敏度相對均方誤差Es的計算公式為

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

(14)

式中:SF是場量提取法計算得到的靈敏度;SR是旋轉(zhuǎn)變換法計算得到的靈敏度。上述兩種方法計算結(jié)果的相對均方誤差分別為0.013 2、0.020 6、0.005 3。可見,本文所述方法與傳統(tǒng)方法的計算結(jié)果高度吻合。

3.3 圖像重構(gòu)實測實驗

采用12電極ECT傳感器進行重構(gòu)圖像實測實驗。電容測量電路為文獻[15]中所述的交流法微小電容測量電路。實驗時,將不同直徑的有機玻璃棒擺放在管道內(nèi)不同位置或在管道內(nèi)注水來模擬不同流形時管道內(nèi)物場分布。設(shè)置三種物場分布代表三種流型進行實驗,分別為:一根粗介質(zhì)棒(核心流,直徑20 mm)、一根偏心介質(zhì)棒(偏心流,直徑20 mm)、三根介質(zhì)棒(每根直徑16 mm)。這三種物場分布的二值灰度圖像如圖8所示。

(a)中心介質(zhì)棒 (b)偏心介質(zhì)棒 (c)對稱三介質(zhì)棒圖8 物場分布設(shè)置

采取廣泛使用的OIOR-Landweber(Offline Iteration Online Reconstruction Landweber)[16]法對測得的電容值進行圖像重構(gòu),用兩種方法計算的靈敏度系數(shù)矩陣重構(gòu)的圖像如圖9所示,迭代次數(shù)為3 000次。

(a) 場量提取法重構(gòu)(b) 場量提取法重構(gòu)(c) 場量提取法重構(gòu)的中心介質(zhì)棒 的偏心介質(zhì)棒 的對稱三介質(zhì)棒

(d) 旋轉(zhuǎn)變換法重構(gòu)(e) 旋轉(zhuǎn)變換法重構(gòu)(f) 旋轉(zhuǎn)變換法重構(gòu)的中心介質(zhì)棒 的偏心介質(zhì)棒 的對稱三介質(zhì)棒圖9 重構(gòu)圖像

規(guī)定在圖8所示的物場分布圖中,介質(zhì)存在區(qū)域灰度為1,其余區(qū)域灰度用于對比兩種靈敏度系數(shù)矩陣重構(gòu)的圖像。為對比重構(gòu)圖像中非零最小元素,采用較常用的相關(guān)系數(shù)與相對均方誤差[17]兩項指標作為評價標準,結(jié)果如表2所示。

表2 圖像重構(gòu)情況對比

可見,基于上述兩種靈敏度系數(shù)矩陣的圖像重構(gòu)結(jié)果吻合度高,證明了前述理論分析的正確性。

4 結(jié) 論

(1)基于矩陣函數(shù)理論,提出了一種基于矩陣函數(shù)旋轉(zhuǎn)變換的ECT二維靈敏度系數(shù)求解方法,即旋轉(zhuǎn)變換法,并到了該方法的數(shù)學(xué)表達式。運用該方法,僅進行一次數(shù)值計算即可獲得ECT傳感器的全部靈敏度系數(shù)。該方法與傳統(tǒng)的場量提取法的靈敏度求解結(jié)果具有高度一致性,二者計算結(jié)果的相對偏差在2.06%以下,證明旋轉(zhuǎn)變換法可以替代傳統(tǒng)的場量提取法,以減少ECT正問題求解的工作量。

(2)針對所提出的旋轉(zhuǎn)變換法,設(shè)計了一種均勻的極坐標成像區(qū)域剖分算法。該算法可保證對電勢(或電場)函數(shù)做旋轉(zhuǎn)變換的過程中不產(chǎn)生新的插值節(jié)點,從而進一步減少了靈敏度矩陣的計算量。

(3)在三種物場分布情況下進行了實驗,分別運用傳統(tǒng)的場量提取法和旋轉(zhuǎn)變化法計算得到的靈敏度系數(shù)矩陣,對相同的ECT測量數(shù)據(jù)進行了圖像重構(gòu)。實驗結(jié)果顯示,兩種方法的圖像相關(guān)系數(shù)和相對誤差高度吻合,表明在實際工程應(yīng)用中,可以使用旋轉(zhuǎn)變化法進行ECT靈敏度系數(shù)求解。

(4)ECT/ERT雙模態(tài)系統(tǒng)是電學(xué)層析成像技術(shù)的研究熱點之一,考慮到ERT的靈敏度系數(shù)表達與ECT靈敏度系數(shù)表達有相同的形式[10],因此旋轉(zhuǎn)變化法對于雙模態(tài)問題的靈敏度求解同樣有效。