馬 敏，劉一斐，劉亞楠

(中國(guó)民航大學(xué) 電子信息與自動(dòng)化學(xué)院，天津 300300)

1 引 言

電容層析成像(electrical capacitance tomography，ECT)技術(shù)問(wèn)世于20世紀(jì)80年代中期。由于其具有非侵入性、安全性、響應(yīng)速度快、安裝便捷以及成本低廉等優(yōu)點(diǎn)，并且與傳統(tǒng)的儀器儀表相比較，電學(xué)層析技術(shù)具有更高的檢測(cè)水平，所以ECT技術(shù)得到了國(guó)內(nèi)外相關(guān)研究人員的廣泛關(guān)注，并應(yīng)用于國(guó)民經(jīng)濟(jì)多個(gè)領(lǐng)域[1～3]。

ECT技術(shù)在成像的過(guò)程中主要分為正問(wèn)題和反問(wèn)題[2]，本文著重于反問(wèn)題的討論。一般在求解數(shù)學(xué)物理反問(wèn)題的過(guò)程中，有兩個(gè)實(shí)際的困難：一是原始數(shù)據(jù)可能與所論問(wèn)題無(wú)關(guān)，則經(jīng)典意義下的近似解可能不存在；二是原始數(shù)據(jù)中小的觀測(cè)誤差會(huì)導(dǎo)致近似解與真實(shí)解嚴(yán)重偏離，因而反問(wèn)題常常是不適定的，很難得到合理的解答與答案。在ECT技術(shù)中成像是通過(guò)將測(cè)量得出的所有電極對(duì)之間的電容值作為投影數(shù)據(jù)，并根據(jù)被測(cè)區(qū)域電磁場(chǎng)分布形成的靈敏度矩陣，利用兩者反演出場(chǎng)域內(nèi)各相的介電常數(shù)分布[4]。由于ECT傳感器提供的數(shù)據(jù)(受電極數(shù)目限制)遠(yuǎn)少于成像像素?cái)?shù)目，導(dǎo)致ECT成像質(zhì)量偏低，分辨率較差，具有嚴(yán)重的欠定性問(wèn)題[5]。

因此在重建圖像的過(guò)程中運(yùn)用合適的圖像重建算法尤為重要。目前已知的典型重建算法大部分是基于l2-范數(shù)的的凸優(yōu)化算法，其特點(diǎn)是通過(guò)平滑圖像使求解穩(wěn)定，并沒(méi)有產(chǎn)生更稀疏的解，導(dǎo)致圖像分辨能力較差[7]。為了得到更加稀疏的解，本文將壓縮感知理論[8]引入ECT圖像重建當(dāng)中，進(jìn)一步提高了圖像分辨能力。

由于傳統(tǒng)信號(hào)的采集與傳輸需要滿(mǎn)足奈奎斯特(Nyquist)采樣定理，即采樣頻率必須大于原始信號(hào)最高頻率的2倍，這樣能夠有效恢復(fù)原始信號(hào)，保證信息的高效傳輸。2006年，加利福尼亞理工學(xué)院的Candes E等人證明了：在圖像的傅里葉表示中，若含有K個(gè)非零傅里葉系數(shù)，滿(mǎn)足系數(shù)個(gè)數(shù)M≥2K，就可以將原始圖像準(zhǔn)確重構(gòu)出來(lái)[9～11]。同年，斯坦福大學(xué)Donoho D L，以及加州大學(xué)的Tao T等在稀疏表示的基礎(chǔ)上，將壓縮和采樣過(guò)程進(jìn)行了合并，壓縮感知理論應(yīng)運(yùn)而生。壓縮感知理論包括3個(gè)步驟：1) 信號(hào)的稀疏表示；2) 設(shè)計(jì)測(cè)量矩陣，要在降低維數(shù)的同時(shí)保證原始信號(hào)x的信息損失最??；3) 設(shè)計(jì)信號(hào)恢復(fù)算法，利用M個(gè)觀測(cè)值無(wú)失真地恢復(fù)出長(zhǎng)度為N的原始信號(hào)。

近年，壓縮感知理論被應(yīng)用到圖像處理、頻譜分析、人臉識(shí)別、無(wú)線(xiàn)通信以及生物醫(yī)學(xué)等眾多領(lǐng)域，在電容層析成像領(lǐng)域的研究成果主要有：

2013年，吳新杰提出一種基于CS理論的ECT流型辨識(shí)方法，根據(jù)待測(cè)樣本和標(biāo)準(zhǔn)樣本稀疏解之間的線(xiàn)性相關(guān)程度來(lái)確定歸屬流型，為ECT流型辨識(shí)技術(shù)的研究提供了一種新的手段[12]。

2014年，王琦在ECT/CT雙模態(tài)成像中引入了壓縮感知理論，建立了一種隨機(jī)采樣模型，并結(jié)合相應(yīng)的l1正則化算法進(jìn)行動(dòng)態(tài)成像，通過(guò)具體實(shí)驗(yàn)得到了多相流介質(zhì)分布的真實(shí)情況[13]。

2017年，張立峰利用FFT基進(jìn)行圖像灰度的稀疏化；并基于高斯隨機(jī)方法設(shè)計(jì)了ECT系統(tǒng)的觀測(cè)矩陣；最后計(jì)算相應(yīng)的雅可比矩陣；研究了基于內(nèi)點(diǎn)法、GPSR算法以等CS重建算法，進(jìn)行ECT圖像重建[14,21]。該文首先選取了合適的稀疏基將原始灰度信號(hào)轉(zhuǎn)化為相應(yīng)的稀疏信號(hào)；其次針對(duì)ECT系統(tǒng)中靈敏度矩陣與稀疏基不相關(guān)性較弱的問(wèn)題，基于高斯隨機(jī)陣對(duì)靈敏度矩陣的各行重新排列，隨后進(jìn)行SVD分解，并加入修正因子得到了列獨(dú)立性更高的觀測(cè)矩陣；最后，采用基于l1/2-范數(shù)的半閾值迭代算法對(duì)ECT逆問(wèn)題進(jìn)行求解，考慮到懲罰項(xiàng)在原點(diǎn)處并不光滑，于是在罰函數(shù)中又加入了l2約束項(xiàng)，改善了原有算法在原點(diǎn)處的光滑性。通過(guò)改進(jìn)的半閾值迭代算法求解，既縮小了圖像誤差又兼顧了成像速度，驗(yàn)證了本文所提算法在ECT成像過(guò)程中的可行性與優(yōu)越性。

2 壓縮感知理論

2.1 信號(hào)的稀疏表示

信號(hào)的稀疏性可以簡(jiǎn)單理解為信號(hào)向量具有較少的非零元素，一般的真實(shí)信號(hào)并非絕對(duì)稀疏，而是在某個(gè)變換域下近似稀疏，這樣的信號(hào)稱(chēng)為可壓縮信號(hào)。

假設(shè)在RN空間中存在一組向量基φi，構(gòu)建一個(gè)基矩陣Ψ=[φ1，φ2，…，φN]，在其上存在一個(gè)具有稀疏度K(即有K個(gè)非零值)的采樣信號(hào)x，x∈RN，即

(1)

式中：α為X在Ψ域上的稀疏表示，α的稀疏度為k，若存在k≤K?N，則說(shuō)明x是稀疏的(即可壓縮的)，Ψ則稱(chēng)為信號(hào)x的稀疏基。要使信號(hào)的稀疏度盡可能減少,需要選取恰當(dāng)?shù)南∈杌?/p>

本文的研究對(duì)象為ECT系統(tǒng)，需要對(duì)管道中的氣固兩相流進(jìn)行圖像仿真重建，初始信號(hào)選擇被測(cè)管道中的圖像灰度信號(hào)，該信號(hào)包含眾場(chǎng)域介質(zhì)分布信息并非稀疏信號(hào)，根據(jù)壓縮感知理論，需要找到合適的稀疏基將初始灰度信號(hào)變?yōu)橄∈栊盘?hào)。

常見(jiàn)的稀疏基有離散余弦變換基(DCT)、離散正弦變換基(DST)、離散傅里葉變換基(DFT)將LBP算法獲得的灰度矩陣作為初始信號(hào)矩陣，經(jīng)不同的稀疏基處理之后效果如圖1所示。

圖1 不同稀疏基的稀疏效果圖Fig.1 Sparse effects of different sparse matrix

橫坐標(biāo)表示為原始信號(hào)經(jīng)稀疏基稀疏后的有用信號(hào)數(shù)量，縱坐標(biāo)表示為有用信號(hào)的數(shù)值。以矩陣零范數(shù)度量信號(hào)的稀疏度，稀疏度在55～78，稀疏效果大致相同。離散余弦變換相當(dāng)于一個(gè)長(zhǎng)度大概是它兩倍的離散傅里葉變換，由于實(shí)偶序列的DFT基是實(shí)對(duì)稱(chēng)的，因此數(shù)據(jù)存在一半的冗余[15]。注意到本文的初始灰度信號(hào)并非實(shí)偶矩陣，經(jīng)過(guò)FFT基變換之后會(huì)產(chǎn)生虛數(shù)，從而對(duì)重構(gòu)算法的要求更高。同DST基相比，DCT基的稀疏效果更好，因此本文選擇DCT基作為初始信號(hào)的稀疏基。

2.2 觀測(cè)矩陣與重構(gòu)算法

觀測(cè)器的采樣目的是在獲得M個(gè)觀測(cè)值的情況下，同時(shí)保證其可以輸出一個(gè)長(zhǎng)度為N的的信號(hào)x對(duì)于信號(hào)測(cè)量的過(guò)程可以記為：

y=Φx

(2)

式中：Φ∈RM×N表示為觀測(cè)矩陣，在ECT系統(tǒng)中則為靈敏度矩陣；y為觀測(cè)值(即電容值)。

如果x是稀疏的，則按照式(1)將x變換到Ψ域，式(2)可寫(xiě)成：

y=Φx=ΦΨα

(3)

選取非相關(guān)性作為ECT系統(tǒng)中靈敏度矩陣的優(yōu)化準(zhǔn)則。因此，需要在靈敏度矩陣的基礎(chǔ)上進(jìn)一步增大與稀疏基的不相關(guān)性，同時(shí)提高靈敏度矩陣的列獨(dú)立性。

考慮到高斯矩陣與大多數(shù)的稀疏基并不相關(guān)，如果將靈敏度矩陣的各行按照高斯隨機(jī)順序重新排列[14]，便可以增大靈敏度矩陣與稀疏基的不相關(guān)性，對(duì)應(yīng)的電容矩陣也做相同變換，因此選取非相關(guān)性作為ECT系統(tǒng)中靈敏度矩陣的優(yōu)化準(zhǔn)則，將靈敏度矩陣的各行按照高斯隨機(jī)順序重新排列，并引入截?cái)嗥娈愔档乃枷隱16]，增大奇異值中較小的數(shù)值，對(duì)變換后的靈敏度矩陣進(jìn)行奇異值分解重新得到：

y=S1Ψα=Aα

(4)

式中：y表示變換后的電容矩陣;S1表示變換后的靈敏度矩陣即新的測(cè)量矩陣;Ψ代表稀疏基;A=S1Ψ代表觀測(cè)矩陣;α是初始灰度信號(hào)G0變換后的稀疏向量。

對(duì)于式(4)的求解過(guò)程稱(chēng)為算法的重構(gòu)，轉(zhuǎn)換為如下形式：

s.t.y=Aα

(5)

由于l0-范數(shù)的求解過(guò)程是把重構(gòu)信號(hào)α的非零項(xiàng)位置逐一列出，利用線(xiàn)性組合求解，是一個(gè)NP-hard問(wèn)題。通常將其轉(zhuǎn)化為其他替代模型利用不同的重構(gòu)算法進(jìn)行求解。

3 l1/2正則化的半閾值迭代算法

針對(duì)y=Aα此類(lèi)具有欠定性特點(diǎn)的優(yōu)化求解問(wèn)題，利用lp代替l0構(gòu)建新的泛函模型，并且p取0.5，則基于l1/2正則化[17]的泛函模型如下：

(6)

(7)

αn+1=H1/2(α+μAΤ(y-Aα))

(8)

H1/2(·)代表閾值函數(shù)，設(shè)其自變量Bμ(α)=α+μAΤ(y-Aα)，存在正實(shí)數(shù)閾值t*>0使得下式成立：

(9)

(10)

g(Bμ(α))=

(11)

文獻(xiàn)[18]中提出若用Bμ(α)k代表Bμ(α)進(jìn)行降序排列的第k個(gè)值，據(jù)此求得λ的取值范圍：

(12)

根據(jù)文獻(xiàn)[19]中的經(jīng)驗(yàn)值，參數(shù)λ被確定為：

(13)

4 l1/2正則化的改進(jìn)半閾值迭代算法

考慮到l1/2正則子在α=0處并不是光滑的，需要對(duì)罰函數(shù)項(xiàng)進(jìn)一步約束，在原點(diǎn)處加入了更加平滑的l2約束項(xiàng)，對(duì)式(11)進(jìn)行修改，重新構(gòu)建基于非凸正則子的ECT模型：

(14)

根據(jù)據(jù)第3節(jié)中的步驟重新計(jì)算一階最優(yōu)性條件，并加入稀疏信號(hào)變換可得：

(15)

(16)

于是得到原點(diǎn)處的迭代公式：

(17)

根據(jù)文獻(xiàn)[20]，原點(diǎn)處的優(yōu)化計(jì)0算結(jié)果如下：

(18)

(19)

λ1為l2-范數(shù)優(yōu)化約束下的正則化參數(shù)，在ECT系統(tǒng)中，通常情況下按照經(jīng)驗(yàn)值取λ1=0.000 3?；趌1/2正則子的參數(shù)λ可利用交叉驗(yàn)證法進(jìn)行選取，將Bμ(α)代入式(10)中，可知存在半閾值函數(shù)需滿(mǎn)足：

(20)

若用Bμ(α)k代表Bμ(α)進(jìn)行降序排列的第k個(gè)值，并結(jié)合ECT系統(tǒng)中的數(shù)據(jù)計(jì)算，進(jìn)一步可得λ=0.001 5。

如果ε代表一個(gè)較小的常數(shù)值，改進(jìn)的半閾值迭代算法最終可寫(xiě)成如下形式：

(21)

ECT系統(tǒng)中基于l1/2正則化的半閾值迭代算法求解步驟如下：

1) LBP算法得到的灰度矩陣作為初始矩陣，并采用DCT基進(jìn)行稀疏化；

2) 將靈敏度矩陣S的各行向量按高斯隨機(jī)順序重排，隨后進(jìn)行奇異值分解，加入修正因子增大觀測(cè)矩陣的奇異值，從而得到優(yōu)化后的觀測(cè)矩陣；

3) 改變罰函數(shù)，構(gòu)建基于l1/2正則化的ECT圖像優(yōu)化模型；

4) 初始化，設(shè)置初始迭代矩陣為零矩陣；

5) 采用半閾值迭代算法對(duì)優(yōu)化模型進(jìn)行求解；

6) 結(jié)合匹配追蹤策略，加速迭代收斂過(guò)程；

7) 設(shè)置迭代終止條件α(k+1)-α(k)

8) 將精確解α*代入G=ΨDCTα*得到重建后的圖像灰度信號(hào)。

基于改進(jìn)半閾值迭代算法流程圖如圖2所示。

圖2 基于l1/2正則化的改進(jìn)半閾值迭代算法的流程圖Fig.2 Diagram of improved half threshold iterative algorithm based on l1/2 regularization

5 算法成像效果

在實(shí)驗(yàn)仿真的過(guò)程中，將LBP算法獲得的灰度矩陣作為初始信號(hào)矩陣，選擇DCT基作為初始信號(hào)的稀疏基對(duì)初始信號(hào)進(jìn)行稀疏；將靈敏度矩陣的各行按照高斯隨機(jī)順序重新排列，并引入截?cái)嗥娈愔档乃枷?，增大奇異值中較小的數(shù)值，從而得到新的奇異值矩陣和觀測(cè)值，重新建立了基于壓縮感知理論的ECT圖像重建問(wèn)題。

本次仿真使用COMOSL 5.5軟件對(duì)16電極管道的ECT傳感系統(tǒng)建模，空?qǐng)鲋胁牧?空氣)的相對(duì)介電常數(shù)設(shè)置為1.0；電極和屏蔽罩(銅)的相對(duì)介電常數(shù)設(shè)置為2.2；管道材料(塑料)相對(duì)介電常數(shù)設(shè)置為5.8；管道內(nèi)被測(cè)物體的相對(duì)介電常數(shù)設(shè)置為3.0。分割網(wǎng)絡(luò)設(shè)置為64*64，共3 228個(gè)有效單元。利用MATLAB 2014a進(jìn)行圖像重建和評(píng)估，仿真流型選擇核心流、二泡流、四泡流、層流、扇形流、環(huán)流和雙U型。

不同算法的仿真成像效果如表1所示。

表1 改進(jìn)半閾值算法與其他算法成像效果對(duì)比Tab.1 Comparison of image effect of improved half threshold algorithm and other algorithm

選取直接算法中的Tikhonov正則化、迭代類(lèi)算法中的Landweber和基于壓縮感知理論的半閾值迭代類(lèi)算法進(jìn)行對(duì)比實(shí)驗(yàn)。

比較不同算法成像效果可得：ECT典型成像算法中的Tikhonov正則化和Landweber迭代算法在泡狀流和層流中的成像效果較好，但是受到周?chē)姌O影響，圖像會(huì)產(chǎn)生明顯的粘連以及偽影現(xiàn)象；在環(huán)形、雙U型中成像效果較差，無(wú)法正確分辨被測(cè)物體的實(shí)際形狀和大小。半閾值迭代算法在各個(gè)流型中均不會(huì)產(chǎn)生偽影和粘連現(xiàn)象，幾乎可以完整復(fù)現(xiàn)被測(cè)物體的實(shí)際位置和形狀，但是在多泡流中所呈現(xiàn)圖形的大小稍有欠缺。而改進(jìn)的半閾值迭代算法同上述傳統(tǒng)的算法相比可以有效改善多泡流中的偽影粘連現(xiàn)象，與半閾值迭代算法相比可以更準(zhǔn)確復(fù)現(xiàn)被測(cè)物體的位置和大小。

為了更加客觀地評(píng)價(jià)不同算法的成像效果，定量反映成像質(zhì)量的優(yōu)劣，選取圖像誤差(image error，IME)和圖像相關(guān)系數(shù)(correlation coefficient，CORR)作為評(píng)價(jià)指標(biāo)，并計(jì)算了算法運(yùn)行時(shí)間。

(22)

(23)

表2 不同算法的圖像誤差Tab.2 Image error of different algorithms

表3 不同算法的圖像相關(guān)系數(shù)Tab.3 Image correlation coefficients of different algorithms

表4 不同算法的成像時(shí)間Tab.4 Imaging time of different algorithms s

通過(guò)分析表2～表4可得：大多數(shù)流型中，Landweber迭代算法成像時(shí)間最長(zhǎng)，Tikhonov正則化算法次之；所有流型中，半閾值迭代算法最短，改進(jìn)半閾值算法與半閾值時(shí)間相近。同時(shí)Landweber迭代算法和Tikhonov正則化法圖像誤差相差不大，但Tikhonov正則化的實(shí)時(shí)性更好。半閾值迭代算法相較于傳統(tǒng)的算法，圖像誤差較小，成像的精度更高，而且成像速度更快，即使與非迭代算法中的Tikhonov正則化相比，成像時(shí)間也僅為后者的六分之一，改進(jìn)后的半閾值迭代算法，雖然成像時(shí)間稍微有所增加，但是在一定程度上減小了圖像誤差，且成像時(shí)間仍舊比其他算法的量級(jí)小很多，從圖3～圖5可以更直觀看出改進(jìn)半閾值迭代算法的優(yōu)勢(shì)。綜合考慮系統(tǒng)實(shí)時(shí)性和重建質(zhì)量的要求，改進(jìn)半閾值迭代算法的優(yōu)勢(shì)更為明顯，改進(jìn)的半閾值迭代算法與所有對(duì)照組的對(duì)比結(jié)果最優(yōu)，證明了該算法在求取稀疏解時(shí)強(qiáng)大的尋優(yōu)能力。

圖3 各算法圖像誤差圖對(duì)比Fig.3 Comparison of image error of each algorithm

圖4 各算法圖像相關(guān)系數(shù)對(duì)比Fig.4 Comparison of image correlation coefficient of each algorithm

圖5 各算法圖像成像時(shí)間對(duì)比Fig.5 Comparison of imaging time of each algorithm

6 結(jié) 論

電容層析成像逆問(wèn)題的求解存在欠定性和病態(tài)性，傳統(tǒng)的直接算法和迭代類(lèi)算法往往無(wú)法兼顧求解速度和成像質(zhì)量。本文將壓縮感知理論引入ECT成像逆問(wèn)題的研究中，通過(guò)稀疏化初始信號(hào)、重構(gòu)觀測(cè)矩陣以及重構(gòu)算法等過(guò)程完成了ECT圖像重建的任務(wù)，有效改善了系統(tǒng)的欠定性問(wèn)題。通過(guò)對(duì)ECT系統(tǒng)的泛函模型進(jìn)行修正，采用改進(jìn)的半閾值迭代算法求解，精確復(fù)現(xiàn)了初始信號(hào)，且系統(tǒng)具有較好的實(shí)時(shí)性。同傳統(tǒng)算法和部分壓縮感知中的算法相比，成像精度較高，同時(shí)能夠兼顧成像速度，顯示出本文所提算法在ECT系統(tǒng)中的良好應(yīng)用。