張立峰, 盧棟臣, 劉衛(wèi)亮

(華北電力大學(xué) 自動(dòng)化系，河北 保定 071003)

1 引 言

兩相及多相分布廣泛存在于石油、化工、生物制藥等工業(yè)生產(chǎn)中。電容層析成像(electrical capacitance tomography, ECT)技術(shù)作為一種無損檢測技術(shù),能有效對(duì)各相分布進(jìn)行可視化重構(gòu)。該系統(tǒng)具有非侵入、結(jié)構(gòu)簡單、快速成像等特點(diǎn)[1]。ECT“軟場”特性導(dǎo)致邊界電容值與介電常數(shù)分布為非線性關(guān)系,微小擾動(dòng)將導(dǎo)致逆問題求解的病態(tài)性。如何精確、快速地獲得重建圖像是ECT系統(tǒng)亟待解決的問題。傳統(tǒng)重建算法可分為兩類:非迭代與迭代類算法[2]。前者實(shí)時(shí)性好,但其所獲得的圖像精度低、形狀保真度差,代表算法為LBP算法;后者則成像精度明顯提升,但對(duì)于中心-邊緣分布及諸多復(fù)雜分布所獲得的成像質(zhì)量仍難以滿足要求,成像速度較慢,代表算法有Landweber算法[3]。近年來,基于壓縮感知原理的圖像重建得到了快速發(fā)展[4]。在實(shí)際工業(yè)生產(chǎn)中許多兩相流分布具有稀疏性,王琦等根據(jù)這一原理對(duì)ECT/CT雙模態(tài)圖像進(jìn)行重構(gòu),使重建質(zhì)量得到提升[5];Ye J M等提出了構(gòu)建介電常數(shù)稀疏基的方式,將ECT待重建圖像表示為基的線性組合,以提升相鄰多分布的重建相關(guān)度[6];文獻(xiàn)[7]采用徑向基函數(shù)的方法,將相應(yīng)的優(yōu)化問題轉(zhuǎn)化為尋找最佳相位邊界和分布的實(shí)際介電常數(shù),該方法既可保證重建質(zhì)量又能增強(qiáng)重建的魯棒性。

壓縮感知在解決稀疏性不適定性問題中通常施加懲罰項(xiàng),可確保解的穩(wěn)定性,由此衍生出許多基于正則化的ECT重建方法,如LASSO方法,全變分(total variation,TV)正則化,均取得了較好的重建效果[8]。而松弛因子法是一種通過構(gòu)建目標(biāo)函數(shù)來量化解的稀疏度的方法,近年來得到眾多學(xué)者的青睞。Huang S等將熵函數(shù)作為松弛因子以降低稀疏信號(hào)恢復(fù)的誤差[9];Zhang L F等利用貝葉斯重加權(quán)的方法,構(gòu)造新的松弛目標(biāo)函數(shù)并結(jié)合ECT反問題模型進(jìn)行求解,為ECT圖像重建提出了一種新方法[10]。由于所得稀疏解并不是絕對(duì)稀疏的,通過設(shè)置閾值的方法對(duì)重建得到的廣義稀疏解進(jìn)行處理,能進(jìn)一步提高重建圖像精度。常見的優(yōu)化方法有二值處理、設(shè)置濾波器及聚類尋優(yōu)等[11]。在電學(xué)層析成像領(lǐng)域,He等[12]通過密度聚類的方法對(duì)灰度值進(jìn)行閾值處理,提高了重建圖像質(zhì)量。

高斯混合模型是通過計(jì)算樣本點(diǎn)在任一模型中的后驗(yàn)概率來判斷樣本點(diǎn)的類別[13],而且,它在處理樣本點(diǎn)時(shí)可保證各維度的不相關(guān)性,因此對(duì)于ECT稀疏解,直接在一維空間聚類尋優(yōu),可減少計(jì)算復(fù)雜度,提升聚類效果[14]。本文基于松弛因子恢復(fù)稀疏信號(hào)算法和聚類尋優(yōu)降噪的思想,提出一種將非凸熵函數(shù)作為松弛因子并結(jié)合高斯混合模型進(jìn)行后處理的方法。先將原始問題轉(zhuǎn)換為一系列重加權(quán)子問題進(jìn)行求解,而后利用高斯混合模型聚類尋優(yōu)確定相應(yīng)閾值點(diǎn),從而提升重建圖像的效果。仿真及實(shí)驗(yàn)均驗(yàn)證了本文方法的有效性。

2 基本理論

2.1 壓縮感知ECT模型

12電極ECT工作原理如圖1所示。對(duì)于n電極ECT傳感器,其獨(dú)立測量值數(shù)目為N=n(n-1)/2,12電極ECT傳感器,可得66個(gè)不同的測量電容值。經(jīng)離散化、線性化后得到ECT數(shù)學(xué)模型:

圖1 12電極ECT工作原理示意圖Fig.1 Schematic diagram of working principle 12-electrode ECT

C=Sg+e

(1)

式中:C為N×1維的歸一化電容測量值矩陣;S為N×M維的歸一化靈敏度矩陣;g為M×1維的歸一化灰度值;e為N×1維的歸一化測量噪聲。M、N分別表示有限元剖分網(wǎng)格數(shù)和測量電容數(shù)。

由于大多數(shù)分布中歸一化灰度值是稀疏的,因此可通過壓縮感知原理重構(gòu)原始信號(hào)。常見壓縮感知的數(shù)學(xué)模型如下所示:

(2)

式中:y為測量信號(hào);A為測量矩陣;x為待恢復(fù)信號(hào);λ為正則化參數(shù),用來衡量正則化項(xiàng)在數(shù)據(jù)保真度中的權(quán)重。

l0范數(shù)可以很好地保證解的稀疏性,但l0求解為非確定多項(xiàng)式(NP,nondeterministic polynomially)問題,通常采用l1范數(shù)近似代替l0范數(shù)來度量稀疏向量的稀疏度,將原問題轉(zhuǎn)換為凸優(yōu)化問題進(jìn)行求解,由此可得基于壓縮感知的ECT求解模型為:

(3)

從式(3)可以看出,通過壓縮感知求解模型可得到滿足條件的稀疏解,緩解ECT病態(tài)性逆問題。但通過該方式所得重建圖像仍不能滿足要求,因此研究人員通過不斷改進(jìn)模型或提出新的求解算法以提升重建質(zhì)量[15]。本文介紹了一種通過引入松弛變量作為懲罰項(xiàng)的ECT重建模型,并采用FISTA算法對(duì)模型進(jìn)行求解。

2.2 非凸熵ECT模型

20世紀(jì)40年代,香農(nóng)借鑒了熱力學(xué)的概念,把信息中排除了冗余后的平均信息量稱為“信息熵”,并給出了計(jì)算信息熵的數(shù)學(xué)表達(dá)式[16]:

(4)

式中:Q為隨機(jī)變量,由{q1,q2,…,qn}組成;p(i) 為Q=qi時(shí)的概率;log用來保證信息的可加性。由于熵函數(shù)是1個(gè)凸函數(shù),Donoho等試圖通過最小化香農(nóng)熵函數(shù)來解決下面的問題[17]:

(5)

通過lp范數(shù)歸一化X構(gòu)造廣義的香農(nóng)熵函數(shù)可以有效解決稀疏恢復(fù)問題。

(6)

目標(biāo)函數(shù)hp(X)是非凸的。文獻(xiàn)[9]中證明了最大化hp(X)可保證解的稀疏性,將目標(biāo)函數(shù)與ECT壓縮感知模型相結(jié)合得到基于非凸熵(nonconvex entropy, NE)ECT求解模型:

(7)

2.3 高斯混合模型

高斯混合模型(Gaussian mixture model, GMM)使用高斯分布作為參數(shù)模型[18],通過計(jì)算每個(gè)樣本后驗(yàn)概率可以對(duì)數(shù)據(jù)進(jìn)行無監(jiān)督分類,相應(yīng)的概率模型為:

(8)

(9)

式中:σ為數(shù)據(jù)標(biāo)準(zhǔn)差;μ為期望值。分別計(jì)算每個(gè)樣本點(diǎn)在R個(gè)高斯分布中的后驗(yàn)概率,對(duì)應(yīng)概率值大的樣本點(diǎn)即為相應(yīng)的高斯類別,通過期望最大化算法可不斷更新模型參數(shù),以得到準(zhǔn)確的分類結(jié)果。

基于2.2節(jié)及2.3節(jié)的2種模型,構(gòu)建了應(yīng)用于ECT圖像重建的非凸熵與高斯混合模型(NE-GMM)算法。

3 NE-GMM算法求解

3.1 快速迭代閾值收縮算法

f(gt+1)≈f(gt)+〈gt+1-gt,f(gt)〉

(10)

2) 直接最小化非凸函數(shù)hp(g)較為困難,在t+1次迭代后,可通過|gt|一階展開來近似表達(dá):

hp(gt+1)≈hp(|gt|)

+〈|gt+1|-|gt|,h(|gt|)〉

(11)

對(duì)|gi|求梯度可得:

在近似轉(zhuǎn)換完成后,結(jié)合式(10)、式(11)可以進(jìn)一步得到便于求解的代價(jià)函數(shù):

(12)

式(12)可轉(zhuǎn)換為一系列獨(dú)立的一維問題,求解得出:

(13)

3) 引入Nesterov加速可以加快ISTA算法的收斂速度,得到具有快速收斂特性的FISTA算法,在式(13)求得解的基礎(chǔ)上可優(yōu)化得到:

(14)

(15)

3.2 期望最大化算法

高斯混合模型求解是對(duì)R個(gè)高斯分布,求解最佳均值、方差及權(quán)重的過程。期望最大化算法(expectation maximization algorithm, E-M)可有效解決該問題。

在模型中每個(gè)觀測數(shù)據(jù)gn都對(duì)應(yīng)1個(gè)隱變量Z∈RN,Z即為類別變量,且Zk∈{0,1},疊加起來概率為1,p(Zk=1)=αk。進(jìn)而可得其似然函數(shù)完全數(shù)據(jù)概率分布為:

lnp(g,Z|μ,σ,α)

計(jì)算每個(gè)樣本gn對(duì)應(yīng)隱變量zn的后驗(yàn)分布:

γ(znk)=p(zn=k|gn;μt,σt)

(16)

令參數(shù)θt=(αt,μt,σt),可得出關(guān)于參數(shù)θt的函數(shù)表達(dá)式:

(17)

(18)

(19)

(20)

圖2 NE-GMM算法流程圖Fig.2 Flowchart of NE-GMM algorithm

4 實(shí)驗(yàn)結(jié)果與分析

4.1 仿真實(shí)驗(yàn)

為驗(yàn)證方法的有效性,本文依托Intel(R) Core(TM) i5-6300HQ CPU,4G仿真平臺(tái),對(duì)油、氣兩相進(jìn)行模擬仿真,設(shè)置油和氣的相對(duì)介電常數(shù)分別為3和1。采用MATLAB設(shè)計(jì)仿真程序,共模擬8種分布來考察本文算法對(duì)不同分布的保真性。通過有限元法對(duì)每種分布進(jìn)行正問題求解,可得每種分布的66個(gè)獨(dú)立測量電容值,采用不同算法在MATLAB2021a上進(jìn)行重建測試,其中待重建的像素點(diǎn)個(gè)數(shù)為768。

非凸熵模型中需要對(duì)算法參數(shù)初始化,分別設(shè)置:λ=0.000 01;l為梯度的利普希茨常數(shù);p的選取可根據(jù)文獻(xiàn)[19]的方法。本文選取p=1;r1=1,將LBP算法求得的解作為迭代初始解(g0)。由于FISTA算法求得的近似解中含有部分廣義稀疏數(shù)據(jù),影響重建效果,采用高斯混合模型進(jìn)行聚類,找到相應(yīng)的閾值點(diǎn),對(duì)數(shù)據(jù)進(jìn)行閾值降噪處理,可以提升圖像的重建質(zhì)量。

圖3左側(cè)表示預(yù)先設(shè)定的灰度分布,包含768個(gè)像素點(diǎn)(其值為0或1),右側(cè)表示通過FISTA算法所得重建圖像像素對(duì)應(yīng)的灰度值,非凸熵模型重建所得的解可分為:0或1代表邊界上的解(即兩相的歸一化稀疏解);介于0～1之間的解(即影響重建質(zhì)量的廣義稀疏解)。因此,將高斯聚類類別設(shè)置成三類,即R=3,設(shè)置初始權(quán)重αk=1/3;μk、σk初始值分別在FISTA重建得到的灰度值數(shù)據(jù)中隨機(jī)獲得。同時(shí)為保證所得閾值的準(zhǔn)確性,采用兩次聚類中最大均值(max(μk))點(diǎn)的加權(quán)平均作為閾值點(diǎn)。

圖3 FISTA所得解的分布圖Fig.3 Distribution diagram of FISTA solution

將NE-GMM算法與LBP、Landweber、IHT、ADMM-L1及NE算法進(jìn)行重建對(duì)比,結(jié)果如圖4所示。

圖4 仿真重建結(jié)果Fig.4 Simulation reconstruction results

圖4中,Landweber、IHT、ADMM-L1及NE算法最大迭代次數(shù)分別為500、200、200及200。為定量評(píng)價(jià)重建圖像質(zhì)量,引入圖像相對(duì)誤差Er及相關(guān)系數(shù)Cc作為評(píng)價(jià)指標(biāo)。相對(duì)誤差越小,相關(guān)系數(shù)越大,則重建圖像質(zhì)量越佳。

2種評(píng)價(jià)指標(biāo)具體計(jì)算公式如下所示[20]:

(21)

(22)

表1 圖像相對(duì)誤差Tab.1 Relative error of image

表2 圖像相關(guān)系數(shù)Tab.2 Image correlation coefficient

由圖4可知,LBP算法成像精度普遍較低,圖像容易產(chǎn)生粘連;Landweber算法成像精度優(yōu)于LBP算法,對(duì)簡單分布成像精度較好,成像偽影較少,但對(duì)相鄰較近分布及復(fù)雜分布成像精度低,重建圖像偽影較大;IHT算法作為1種貪婪算法,對(duì)邊緣多物體分布的成像精度優(yōu)于Landweber算法,但對(duì)于相鄰較近的中心分布以及稀疏度較大的復(fù)雜分布,其重建圖像同樣存在形狀保真度低、重建偽影大等問題;ADMM-L1正則化算法是常用的稀疏求解算法,該算法得到的重建圖像在形狀保真度上均優(yōu)于上述3種算法,且重建精度也得到進(jìn)一步提升,但隨著迭代次數(shù)的增加,重建圖像會(huì)產(chǎn)生較大的偽影,且對(duì)稀疏度較大的分布(分布8)往往難以精確重構(gòu)。

NE算法在簡單中心分布(如分布1)成像精度低于Landweber算法,但對(duì)于中心多物體分布和復(fù)雜分布,其重建圖像精度高、形狀保真度好、油氣兩相分割清晰優(yōu)于其他四種算法;同時(shí),與ADMM-L1算法類似,隨著迭代次數(shù)增加,NE算法重建圖像也將產(chǎn)生較大偽影。而本文所提出NE-GMM算法成像質(zhì)量明顯提升。重建圖像偽影明顯減少,且各相分布分割更為清晰,同時(shí)圖像的形狀保真度也進(jìn)一步得到提升。

從表1和表2的評(píng)價(jià)指標(biāo)定量分析得出,對(duì)于每種設(shè)定分布,NE-GMM算法計(jì)算得到的圖像相對(duì)誤差最小,相關(guān)系數(shù)也更接近于1,表明了采用非凸熵模型與高斯混合模型結(jié)合得到的重建圖像質(zhì)量更好。

4.2 靜態(tài)實(shí)驗(yàn)

采用華北電力大學(xué)先進(jìn)測量實(shí)驗(yàn)室的12電極ECT系統(tǒng),通過單激勵(lì)測量模式進(jìn)行靜態(tài)實(shí)驗(yàn)。在靜態(tài)實(shí)驗(yàn)中,采用有機(jī)玻璃棒和塑料顆粒模擬油相,通過12電極電容傳感陣列進(jìn)行數(shù)據(jù)采集,設(shè)置單棒、多棒、環(huán)狀及層狀模型,通過隨機(jī)分布方式來模擬6種分布,分別采用6種算法進(jìn)行成像,重建結(jié)果如圖5所示。

圖5 靜態(tài)實(shí)驗(yàn)結(jié)果Fig.5 Static experimental results

由圖5可見,采用NE-GMM算法重建得到的圖像偽影更小,重建精度更優(yōu),各分布形狀保真度更好,優(yōu)于其他5種算法成像結(jié)果。尤其是對(duì)于包含中心物體的多物體分布,由于敏感場的軟場特性,中心分布一般難以成像,而本文采用的算法,可實(shí)現(xiàn)對(duì)中心物體的重建。

5 結(jié) 論

1) 單一的非凸熵算法不能滿足重建精度的要求,在此基礎(chǔ)上引入高斯混合模型,通過閾值尋優(yōu)的方式來進(jìn)一步提升重建質(zhì)量,重建結(jié)果表明:該算法不僅提升對(duì)中心及多物體分布的重建精度,同時(shí)對(duì)于每種分布的形狀保真度也得到極大的提升。

2) 相較于LBP算法、Landweber算法、IHT算法、ADMM-L1算法以及非凸熵算法,GMM-非凸熵算法在重建精度方面均優(yōu)于其他五種算法,但在實(shí)時(shí)性方向仍需進(jìn)一步改進(jìn)。在后續(xù)研究中將重點(diǎn)研究改進(jìn)該算法,以滿足實(shí)時(shí)性的需求,并探索其在三維ECT重建中的應(yīng)用。