王伊凡， 朱海瀅， 熊嘯岳， 顏 華， 周英鋼

(沈陽(yáng)工業(yè)大學(xué) 信息科學(xué)與工程學(xué)院， 沈陽(yáng) 110870)

電容層析成像(electrical capacitance tomography，簡(jiǎn)稱ECT)技術(shù)是過(guò)程層析成像技術(shù)中的一種，其研究始于上世紀(jì)80年代末期.該技術(shù)基于電容敏感原理，通過(guò)測(cè)量放置在封閉管道或容器外不同極板間的電容值作為投影數(shù)據(jù)，采用適當(dāng)?shù)膱D像重建算法重建出封閉物場(chǎng)內(nèi)部介質(zhì)(介電常數(shù))分布情況的二維或三維圖像[1-2].ECT技術(shù)在諸多方面得到廣泛應(yīng)用，例如鼓泡硫化床氣/固兩相流成像、流動(dòng)的油/水兩項(xiàng)介質(zhì)的動(dòng)態(tài)分布成像、循環(huán)流化床氣/固兩相流動(dòng)三維成像以及熱虹吸管冷凝實(shí)驗(yàn)成像等[3-6].該技術(shù)具有結(jié)構(gòu)簡(jiǎn)單、費(fèi)用低、響應(yīng)速度快、安全無(wú)輻射等優(yōu)點(diǎn)，被認(rèn)為是一種極具發(fā)展前景的可視化測(cè)量技術(shù).

傳統(tǒng)的ECT重建圖像是管道的二維截面圖.近幾年興起的三維ECT重建圖像以三維圖像形式給出.由于三維ECT可以獲得介質(zhì)分布的幾何形狀、相對(duì)位置和體積等信息，因此受到廣泛關(guān)注.但是與傳統(tǒng)的二維ECT相比，三維ECT系統(tǒng)測(cè)量信號(hào)的動(dòng)態(tài)范圍和敏感場(chǎng)強(qiáng)度差異較大，導(dǎo)致靈敏度矩陣的條件數(shù)變大，病態(tài)問(wèn)題也更嚴(yán)重[7]，因此對(duì)于重建算法的要求就更高.

工業(yè)中對(duì)重建速度與精度要求同等重要，因此，對(duì)ECT圖像重建算法的速度和精度要求更高[8].目前，電容層析成像技術(shù)的圖像重建算法總體上可分為非迭代算法和迭代算法兩大類.迭代算法相比非迭代算法在大多數(shù)情況下可獲得較好的重建效果，但是其也具有兩個(gè)比較顯著的缺點(diǎn)：1)運(yùn)算時(shí)算法的收斂速度比較慢，通常需要很多次的迭代運(yùn)算才能獲得比較高質(zhì)量的重建圖像；2)在靈敏度比較高的成像區(qū)域，重建圖像會(huì)產(chǎn)生偽跡現(xiàn)象.特別是在三維ECT圖像重建時(shí)，這兩個(gè)問(wèn)題就變得更為嚴(yán)重.目前，國(guó)內(nèi)外學(xué)者對(duì)于算法的研究主要集中在算法迭代公式的改進(jìn)或靈敏度矩陣的優(yōu)化上，如Gunes等[2]采用離散傅里葉方法進(jìn)行偽三維成像；吳新杰等[9]將卷積神經(jīng)網(wǎng)絡(luò)應(yīng)用在ECT圖像重建上；張立峰[10]將壓縮感知理論和快速傅里葉變換用于圖像重建；Moura等[11]提出了一種稀疏重建方案，采用冗余靈敏度矩陣來(lái)改善重建效果等.而閾值濾波是減少偽跡的一種簡(jiǎn)單而有效的方法，目前，迭代濾波方面的研究工作相對(duì)較少，并且將各改進(jìn)算法與傳統(tǒng)Landweber迭代算法重建結(jié)果進(jìn)行比較，而傳統(tǒng)Landweber迭代算法均采用簡(jiǎn)單的0-1方式進(jìn)行閾值濾波.

本文將Landweber迭代算法和閾值濾波技術(shù)相融合，提出了一種基于自適應(yīng)閾值濾波的三維ECT重建算法(three-dimensional reconstruction algorithm based on adaptive threshold filtering，簡(jiǎn)稱ATF算法).在迭代重建過(guò)程中使用閾值濾波對(duì)重建圖像進(jìn)行強(qiáng)化偽跡抑制操作，并定期自動(dòng)更新濾波閾值.為了獲得最佳濾波閾值，首先定義基于香農(nóng)熵函數(shù)的圖像模糊度量，然后找到使圖像模糊度量最小的閾值作為最佳濾波閾值.

1 三維ECT圖像重建系統(tǒng)和原理

電容層析成像系統(tǒng)由傳感器系統(tǒng)、數(shù)據(jù)采集系統(tǒng)和圖像重建系統(tǒng)三部分組成，如圖1所示.

圖1 ECT系統(tǒng)結(jié)構(gòu)示意圖Fig.1 Schematic ECT system structure

假設(shè)屏蔽罩內(nèi)沒(méi)有自由電荷，則ECT傳感器的靜電場(chǎng)由拉普拉斯方程[12]描述，即

·(ε0ε(x，y，z)φ(x，y，z))=0

(1)

式中：ε0為真空介電常數(shù)；ε(x，y，z)為相對(duì)介電常數(shù)分布；φ(x，y，z)為電勢(shì)分布.電容與介電常數(shù)之間的關(guān)系式為

(2)

式中：C0為電容值；Q為檢測(cè)極板上的電荷；U為電勢(shì)差；E為電場(chǎng)強(qiáng)度；Γ為包圍極板的封閉曲面.測(cè)量電容值和介質(zhì)分布之間的關(guān)系是非線性的，為了計(jì)算方便，將二者的非線性關(guān)系按照線性化近似可得

C=SG

(3)

式中：C∈RN×1為歸一化的電容矢量；G∈RM×1為歸一化的介電常數(shù)矢量；S∈RN×M為歸一化的靈敏度矩陣.對(duì)于J極板構(gòu)成的ECT系統(tǒng)，可獲得獨(dú)立測(cè)量電容數(shù)目為N=J(J-1)/2.

2 帶自適應(yīng)閾值濾波的三維ECT重建算法

本文算法的重建公式可描述為

G0=P(STC)

(4)

Gk=P(Gk-1+αST(C-SGk-1))

(5)

(6)

式中：η為濾波閾值；ST為M×N維的歸一化靈敏度矩陣的轉(zhuǎn)置陣；α為迭代步長(zhǎng).本文取α=2/λmax，λmax為STS的最大特征值，當(dāng)α<1時(shí)，則取α=1.5.式(4)用來(lái)獲取重建圖像的初始值，式(5)為迭代重建公式，式(6)為閾值濾波算子，算法具體實(shí)施步驟如下：

1) 令η=0，用式(4)、(5)計(jì)算圖像初始值G0=[G0(1)，G0(2)，…，G0(M)]，并令k=1，其中，M為歸一化的介電常數(shù)矢量個(gè)數(shù).

2) 計(jì)算Gk=P(Gk-1+αST(C-Gk-1)).

3) 如果((k))q≠1，((k))q表示k對(duì)q取余數(shù)，q為預(yù)先設(shè)定好的一個(gè)正整數(shù)，則轉(zhuǎn)到步驟4)；否則按下述方法求η.

將0～1的灰度值均勻地劃分為0～L-1，即L(L為預(yù)先設(shè)定好的一個(gè)正整數(shù)，且為2的冪指數(shù))個(gè)灰度級(jí)；令gmax、gmin分別為圖像的最大、最小灰度級(jí)，h(g)為灰度級(jí)g在圖像中出現(xiàn)的次數(shù)；對(duì)于給定閾值t，基于香農(nóng)熵定義閾值t所對(duì)應(yīng)的圖像模糊度量E(t)為

(7)

F(g)=-μ(g)ln[μ(g)]-[1-μ(g)]·

ln[1-μ(g)]

(8)

(9)

(10)

(11)

式中，int[x]為對(duì)實(shí)數(shù)x四舍五入取整數(shù)值.令E=1、t=gmin、topt=gmin，取值采用如下步驟：

① 利用式(7)～(11)計(jì)算閾值t所對(duì)應(yīng)的圖像模糊度量E(t)，如果E(t)

② 如果t=gmax-1，則轉(zhuǎn)到步驟③；如果t

4) 利用式(12)對(duì)圖像Gk進(jìn)行強(qiáng)化偽跡抑制，其表達(dá)式為

Gk=P(Gk)

(12)

通過(guò)仿真實(shí)驗(yàn)發(fā)現(xiàn)，η值計(jì)算時(shí)間的長(zhǎng)短與灰度級(jí)L的大小有關(guān)，即灰度級(jí)L的取值越大，η值的計(jì)算時(shí)間就越長(zhǎng).但是算法的重建精度并不總是隨著灰度級(jí)L的取值增大而減少.因此，通過(guò)仿真實(shí)驗(yàn)確定灰度級(jí)L=64.

本文中，q值越小，對(duì)偽跡抑制效果越好，但是所消耗的重建時(shí)間也越長(zhǎng).通過(guò)仿真實(shí)驗(yàn)，本文取q=10.

3 算法驗(yàn)證

對(duì)本文提出的三維重建方法分別進(jìn)行了仿真和實(shí)驗(yàn)驗(yàn)證.采用12極板三維ECT傳感器，實(shí)物圖如圖2所示.傳感器的極板軸向上分為三層，每一層上布置4個(gè)極板，相鄰兩層極板水平方向上旋轉(zhuǎn)45°.每一個(gè)極板均由寬為32.725 mm，高為9 mm的銅箔構(gòu)成，水平方向上兩個(gè)相鄰極板間的距離為6.545 mm，軸向上相鄰兩層極板間的距離為9 mm；實(shí)驗(yàn)使用的絕緣有機(jī)玻璃管道內(nèi)半徑為23 mm，外半徑為25 mm，管壁厚度為2 mm，管壁相對(duì)介電常數(shù)約為4；外部屏蔽罩的內(nèi)半徑為33 mm.

圖2 三層ECT傳感器實(shí)物剖面圖Fig.2 Profile of three-layer ECT sensor

3.1 仿真重建

為了驗(yàn)證算法的有效性，采用本文算法與傳統(tǒng)的Landweber迭代算法對(duì)U、H、L三種字母型復(fù)雜分布進(jìn)行了仿真重建.表1為原型分布圖，以及傳統(tǒng)Landweber迭代算法和本文算法的重建圖像.其中，傳統(tǒng)Landweber迭代算法的迭代次數(shù)為3 000次，本文算法迭代次數(shù)為800次.兩種算法的迭代步長(zhǎng)皆為α=1.334 2.表1列出了兩種算法的空間圖像誤差(SIE)、體積誤差(VE)和介電常數(shù)誤差(PE).

表1 仿真成像及誤差Tab.1 Simulated images and errors

極板覆蓋區(qū)域稱為成像區(qū)域，設(shè)e為成像區(qū)域中的第e個(gè)像素單元，T為成像區(qū)域內(nèi)總的像素單元數(shù)，即e=1，2，…，T，則空間圖像誤差(SIE)的計(jì)算公式為

(13)

式中：

空間圖像誤差(SIE)包括重建物體在形狀、位置和體積上的空間誤差信息，因此，SIE的值越小越好.

體積誤差(VE)反映了實(shí)際物體和三維重建物體在體積上的差異.VE的絕對(duì)值越小越好，其計(jì)算公式為

(14)

式中，V(e)為成像區(qū)域內(nèi)第e個(gè)像素單元的體積.

介電常數(shù)誤差(PE)反映了實(shí)際物體的歸一化介電常數(shù)平均值與三維重建物體歸一化介電常數(shù)之間的差異.介電常數(shù)誤差的絕對(duì)值越小，重建圖像的質(zhì)量就越好，其計(jì)算公式為

(15)

本文提出的三維重建算法迭代800次的重建偽跡明顯少于Landweber算法迭代3 000次時(shí)的重建偽跡，重建圖像更接近于原型分布.從重建誤差來(lái)看，本文提出的三維重建算法的SIE、VE和PE值均優(yōu)于Landweber算法值.通過(guò)仿真分析，本文提出的三維重建算法迭代800次和Landweber算法迭代3 000次時(shí)，對(duì)于同一重建對(duì)象的程序運(yùn)行時(shí)間比值約為1∶3.因此，本文提出的三維重建算法對(duì)于偽跡抑制是有效的，既提高了重建圖像的質(zhì)量，又加快了重建速度.

3.2 實(shí)驗(yàn)重建

實(shí)驗(yàn)所用ECT系統(tǒng)如圖3所示，實(shí)驗(yàn)所使用的測(cè)試物是由U、H、L型字母分布的介電常數(shù)約為3.7的石英砂構(gòu)成.圖4為三種字母型測(cè)試物的形狀以及在管道中與極板的相對(duì)位置.

圖3 基于WK65120B的ECT系統(tǒng)Fig.3 ECT system based on WK65120B

表2為傳統(tǒng)Landweber迭代算法和本文所提三維重建算法的重建圖像，可以看出，兩種算法均可以重建出復(fù)雜分布的基本形狀.相比之下，本文所提算法的重建偽跡相對(duì)較少，進(jìn)一步驗(yàn)證了本文所提算法是有效的、可行的.

4 結(jié) 論

本文將Landweber迭代算法和模糊閾值濾波技術(shù)相結(jié)合，提出了基于自適應(yīng)閾值濾波的三維ECT重建算法，在迭代過(guò)程中采用閾值濾波方式

圖4 重建實(shí)物圖Fig.4 Reconstruction of real objects

表2 實(shí)驗(yàn)重建圖像Tab.2 Reconstructed images in experiments

對(duì)所獲得的重建圖像進(jìn)行強(qiáng)化偽跡抑制操作，并定期調(diào)整濾波閾值.仿真和實(shí)驗(yàn)結(jié)果表明，與傳統(tǒng)Landweber迭代算法相比，具有偽跡抑制作用的重建算法能夠顯著減少重建圖像中的偽跡，縮短重建時(shí)間，該方法具有良好的應(yīng)用潛力.