壓縮感知在ECT分相含率檢測中的應用

唐晨暉,胡紅利,王 格，張 肖

(1.西安交通大學 電力設備電氣絕緣國家重點實驗室，陜西 西安 710049；2.西安市產品質量監督檢驗院，陜西 西安 710065)

多相流作為一種混合流動體系廣泛存在于自然界和工業生產中,尤其在電力、石油、化工、冶金、航空航天等關系到國計民生的專業領域。而多相流中又以兩相流最為常見,如氣/液兩相流,氣/固兩相流,油/水兩相流等,其中，油水兩相流廣泛存在于石油開采、管道輸送等工業過程中[1-2],屬液液兩相流范疇，在國內其研究進展遠遠落后于氣液兩相流、氣固兩相流及液固兩相流[3]。

電容層析成像(electrical capacitance tomography,ECT)技術通過測量隨流過極板間流體組分含率變化而變化的極板間電容值來重構流型圖。當流過極板間流體組分含率變化時,混合流體的等價相對介電常數會隨之發生變化,從而導致極板間的電容值發生相應的改變,便可將兩相流中分相含率的信息以電容值的變化表現出來[4]。

圖像重建算法是ECT的關鍵,直接關系到成像精度和成像速度。由ECT靈敏度理論可知,介質分布的微小變化對敏感場分布的影響可忽略,于是，介質分布與測得電容值之間可用簡易的線性映射關系表示[5]。根據電磁場理論,由微小介質擾動方法可推導出靈敏度系數矩陣的場域表達式,從而計算靈敏度系數矩陣,即ECT正問題的求解。而ECT逆問題的求解是根據靈敏度系數矩陣計算介質分布向量,從而進行圖像重構。但由于靈敏度系數矩陣往往非可逆矩陣,因此，ECT圖像重構算法需要逼近靈敏度系數矩陣的廣義逆，或利用迭代算法計算出介質分布向量的近似解。

壓縮感知理論(compressed sensing, CS)由Donoho于2006年基于稀疏分解及逼近論[6]提出,其主要內容可表述為:若一個信號是可壓縮的或在某個變換域內稀疏,則可用一個與變換矩陣非相干的測量矩陣將高維信號線性投影為低維觀測向量,然后，通過求解一個稀疏最優化問題，便能夠從低維觀測向量中精確地重構出原始信號[7]。壓縮感知理論的提出,一定條件下突破了采樣定律,在遠小于Nyquist采樣頻率的條件下用隨機采樣獲取的離散樣本,通過非線性重建算法便可恢復出原信號。目前,CS理論的研究已涉及軍事、醫療、工業、社會生活等諸多領域[8]。

目前，多相流含水率在線測量研究較多的是密度法、射線法、電導率法、電容法以及電磁波法[9],幾種測量方法都有一定的適用范圍和使用條件。相比于其他方法,基于介電常數原理的電容法由于其結構簡單,非侵入,成本低,測量精度高等優點在含水率檢測中被廣泛采用。

本文將CS理論運用到ECT圖像重構問題中,將低維觀測向量轉換為低維稀疏向量,通過求解稀疏最優化問題，較為精確地重構出原信號,從而解決ECT系統的欠定性問題。利用最優閾值算法對重構圖像進行灰度處理,由灰度圖像對分相含率進行計算。

1 壓縮感知基本理論

1.1 信號的稀疏化處理

如果N維信號x∈R最多只有k個分量不為零,即‖x‖0≤k,則該信號是k稀疏的;或者存在一個稀疏域Ψ,x可經其轉換,得到一個稀疏信號,即x=Ψz,其中，‖z‖0≤k[10]。設{Ψ1,Ψ2,…,ΨN}為RN的一組正交基向量,則RN空間中的任何一個N×1的離散信號x都可以線性地表示為

(1)

其中，Ψ={Ψ1,Ψ2,…,ΨN}為N×N的稀疏基矩陣,z為x在Ψ域中的N×1維稀疏化表示向量;若z中只有少部分元素取值很大,大部分取值很小,那么可以丟棄取值較小的元素,只需要少量取值大的元素就可以較好地逼近信號x,此時，稱x在Ψ域中是可壓縮的。

1.2 信號的壓縮采樣

假設Φ為Μ×N維(Μ?N)的觀測矩陣,則使用Φ對Μ×1的高維信號x進行采樣,可以得

y=Φx=ΦΨz。

(2)

其中，y為M×1的低維觀測向量。令A=ΦΨ作為傳感矩陣,則式(2)可寫為

y=Αz。

(3)

高維信號的求解過程就是由M×1維觀測向量y先從式(3)中求出N×1維的稀疏向量z,再根據式(1)線性求解N×1維原始信號x的過程。當Μ

Candes指出,當觀測矩陣Φ與稀疏基Ψ不相關,傳感矩陣Α滿足等距約束條件[11](restricted isometry property, RIP)，從而可求解欠定方程式(3)。

傳感矩陣Α的等距約束條件可表示為

(4)

其中，δk為使得上式成立的最小值,稱為傳感矩陣的等距約束常數。

文獻[12]證明,選擇高斯隨機矩陣作為觀測矩陣Φ可以滿足其與稀疏基Ψ的不相關性,從而使傳感矩陣Α滿足RIP。

1.3 最優化問題求解

最優化問題求解是壓縮感知理論的最后一個階段。當x在Ψ域中k稀疏,且當傳感矩陣Α的2k階等距約束常數δ2k<1時,式y=Αz可轉化為L0范數約束最優化求解問題:

(5)

由式(5)求解出稀疏向量z的唯一精確解是最直接的方法,但由于L0范數具有高度非凸性,為NP-hard問題[13]。這樣的問題需要組合搜索,無法求出數值解。因而,需要尋求其他的替代模型及其相應的重建算法來求解出唯一精確的稀疏系數向量z。

松弛方法是基于L1范數最小化的思想。當傳感矩陣A滿足RIP條件,求解凸松弛的L1范數最小化與求解非凸的L0范數最小化問題是等價的[14],從而把式(5)轉化為L1范數約束最優化求解問題,即

(6)

對于凸優化問題的求解算法目前常用的有內點法、梯度投影算法(gradient projectionfor sparse reconstruction, GPSR)以及不動點連續算法(fixed point continuation,FPC)。FPC作為一種新的基于算子分裂和連續的迭代算法,計算復雜度低、收斂速度快,能夠解決大規模問題[15]。因此，本文選用FPC算法進行凸優化問題的求解。

不動點連續算法主要解決L1范數問題。通過引入正則化參數μ將式(6)約束化,得到去約束凸優化問題:

(7)

H(z)=g(z)+μf(z)=

(8)

則根據凸優化理論[16],H(z)的最小值等價于求T(z)=?H(z)/?z=0。根據算子分裂理論,若H(z)可被分裂成兩個凸函數相加的形式,即H=H1+H2,則有T=T1+T2。根據不動點定理，

0∈T?0∈(z+τT1z)-(z-τT2z)，

(9)

推導得

z=(I+τT1(z))-1(I-τT2(z))z，

(10)

由此得到稀疏向量的迭代公式為

zn+1=(I+τT1(zn))-1(I-τT2(zn))zn。

(11)

其中，T1(z)=g(z),T2(z)=μf(z),I為恒等映射。對式(11)中的微分進行求解得

(12)

以及

(13)

根據凸優化理論,當且僅當

z=sgn(z-τf(z))⊙

(14)

為式(7)的最優解。因此,由FPC算法求解L1范數最優化問題的迭代公式可簡化為

zn+1=sgn(zn-τf(zn))⊙

(15)

其中，τ，μ均為常數。

2 基于CS的ECT圖像重構

將壓縮感知理論運用于電容層析成像系統時,圖像的灰度向量g便是需要重構的原始信號。為了能滿足CS的使用條件,需要尋找一組合適的正交基向量對需要重建的灰度向量進行稀疏化處理。目前,常用的稀疏基主要包括離散傅里葉變換基(DFT)、離散正弦變換基(DST)、離散余弦變換基(DCT)、離散小波變換基(DWT)等。文獻[17]指出,DFT基可以較好地稀疏化ECT中多種典型流型的灰度向量。而經比較分析,時域基能更好地對觀測矩陣進行稀疏化。因此，本文選取時域基作為稀疏基矩陣,可得

g=ΨTz，

(16)

ECT系統的線性模型為

c=Sg，

(17)

將式(16)代入式(17)得到

c=Sg=SΨTz。

(18)

其中：c為ECT系統中的M×1維電容測量值,作為觀測投影向量;S為M×N維靈敏度矩陣,作為壓縮采樣中的觀測矩陣。

(19)

基于CS的ECT圖像重構問題是典型的最小L0范數問題,可由1.3節所述的最優化問題求解算法FPC進行求解。

重構的圖像向量gopt的元素大小范圍為[0,1],對其進行二值化處理,可使重構圖像更加清晰。本文采用最優閾值算法,搜尋閾值后,利用重構圖像計算出新的電容值,并與測量電容進行比較,不斷逼近以至電容誤差最小。

3 仿真及結果分析

為比較基于CS的ECT圖像重構效果,本文分別選用線性反投影算法(linear back projection, LBP)、經典迭代算法Landweber以及壓縮感知背景下的不動點連續算法(fixed point continuation,FPC)進行圖像重構。利用COMSOL有限元計算仿真軟件建立ECT傳感器模型，如圖1所示。最外層為接地屏蔽層;中間是PVC塑料管道,12塊電極片位于管道外周;最內層是待求的成像區域。ECT傳感器的幾何參數如表1所示。

圖1 ECT傳感器仿真模型Fig.1 ECT sensor simulation model

參數 數量 描述 Rp25mm管道內徑Dp3mm管道厚度Di10mm絕緣層厚度L75mm電極長度α26°電極張角β4°電極間距角εp3.2管壁介電常數VE3.3V激勵電壓

分別利用線性反投影算法(LBP)、經典迭代算法Landweber以及壓縮感知背景下的不動點連續算法(FPC)對分層流、環狀流、偏心流3種典型流型進行圖像重構,結果如表2所示。

為了評價重建圖像的效果,采用相對電容殘差、圖像相關系數、重構時間作為評價指標[19]。分別計算3種算法的重建圖像評價指標，結果如表3～5所示。

表2 3種不同算法重構圖像結果Tab.2 Results of reconstruct image using three different algorithms

表3 3種重構圖像的相對電容殘差

Tab.3 Relative capacitance residuals of reconstructed images under three algorithms

流型電容殘差LBPLandweberFPC 層流2.665 11.502 20.941 7環流1.790 60.132 10.224 6偏心流3.887 21.044 20.937 3

表4 3種重構圖像的圖像相關系數

Tab.4 Image correlation coefficient of reconstructed images under three algorithms

流型圖像相關系數LBPLandweberFPC 層流0.165 10.737 20.786 5環流0.090 60.816 10.811 2偏心流0.107 20.803 20.832 2

表5 3種重構圖像的重構時間

Tab.5 Reconstructing time of reconstructed images under three algorithms

流型重構時間/sLBPLandweberFPC 層流0.015.390.81 環流0.014.720.79偏心流0.015.730.81

由重構圖像及上述評價指標可看出LBP算法的重構圖像基本無法用肉眼識別,辨識度很差;采用Landweber算法重構的圖像邊緣較為模糊,辨識度較好;而基于FPC算法的重構圖像具有更清晰的邊緣信息,重建圖像偽跡較少,辨識度更高,更接近真實分布。此外,FPC算法的重構時間遠遠低于Landweber算法,滿足工業成像對實時性的要求。

為進一步提升圖像的重構質量,采用最優閾值算法對Landweber及FPC算法重構的圖像進行圖像灰度變換。處理后的效果及圖像相關系數如表6所示。

表6 兩種算法的灰度處理圖像與相關系數

Tab.6 Grayscale processing image and correlation coefficient of two algorithms

流型Landweber FPC 灰度處理相關系數灰度處理相關系數層流0.79210.8354環流0.83320.8043偏心流0.76540.9133

從灰度圖像可以看出,經最優閾值算法灰度處理后的圖像效果有所改善,且FPC算法重構圖像優于Landweber算法。

4 截面含水率測量

本文中的分相含率測量以二相流截面含水率為例,用管道截面處待定相的面積與管道截面積之比表示。油/水兩相流中截面含水率的計算公式為

(20)

其中，Aw表示水的截面區域面積,Ao表示油的截面區域面積。則通過重構圖像計算截面含水率的公式可表示[20-21]為

(21)

fi是經最優閾值算法灰度處理后管道截面的像素值,其中，待定相的像素值為1,其余區域像素值為0;M是管道截面總的像素個數。

為了分析基于CS的ECT重構圖像對于不同分相含率流體的計算能力,分別選用含有8mm，5mm，2mm三種不同內徑的待定相流型進行仿真,仿真結果如表7所示。

表7 3種不同分相含率的圖像重構

Tab.7 Image reconstruction of three different phase separation ratios

預設物場FPC灰度處理

預設物場的實際分相含率分別設定為10.24%，4%，0.64%;通過灰度處理后的重構圖像計算得到3種分相含率的大小分別是:9.91%，3.61%，1.9%;相對誤差分別是:3.22%，9.75%，70.31%。由表7中重構圖像和分相含率計算結果可以看出,隨著分相含率的減小,重構圖像邊界開始模糊,物場變形較為嚴重,這極大影響灰度處理中最優閾值算法求解的二值化閾值,從而導致分相含率的相對誤差隨著真實含率的減小而增大。

5 結 論

1)研究了CS理論在ECT問題中的應用，提出了稀疏基的選擇方案，設計了觀測矩陣；并將不動點迭代算法FPC應用于凸優化問題的求解，得出基于FPC的迭代算法。

2)建立了ECT傳感器模型，并分別對3種典型兩相流流型進行仿真重構，仿真實驗說明，基于CS的ECT重構算法在只有少量測量數據的情況下也能精確地重建圖像，解決了ECT的欠定性問題，其成像速度大幅提升，且成像質量遠高于LBP，并優于Landweber迭代算法。

3)為了分析基于CS的ECT重構圖像對于不同分相含率流體的計算能力，以二相流截面含水率為例，分別選用了3種不同內徑的待定相流進行仿真計算，利用最優閾值算法對兩種算法下的重構圖像進行灰度處理。最終得出，分相含率的相對誤差隨著真實含率的減小而增大。因此，最優閾值算法的優化是下一步工作重點。