基于改進半閾值迭代算法的ECT圖像重建

馬 敏，劉一斐，劉亞楠

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

1 引 言

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

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

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

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

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

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

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

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

2 壓縮感知理論

2.1 信號的稀疏表示

信號的稀疏性可以簡單理解為信號向量具有較少的非零元素，一般的真實信號并非絕對稀疏，而是在某個變換域下近似稀疏，這樣的信號稱為可壓縮信號。

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

(1)

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

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

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

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

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

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

觀測器的采樣目的是在獲得M個觀測值的情況下，同時保證其可以輸出一個長度為N的的信號x對于信號測量的過程可以記為：

y=Φx

(2)

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

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

y=Φx=ΦΨα

(3)

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

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

y=S1Ψα=Aα

(4)

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

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

s.t.y=Aα

(5)

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

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

針對y=Aα此類具有欠定性特點的優(yōu)化求解問題，利用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ù)閾值t*>0使得下式成立：

(9)

(10)

g(Bμ(α))=

(11)

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

(12)

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

(13)

4 l1/2正則化的改進半閾值迭代算法

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

(14)

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

(15)

(16)

于是得到原點處的迭代公式：

(17)

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

(18)

(19)

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

(20)

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

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

(21)

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

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

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

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

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

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

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

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

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

基于改進半閾值迭代算法流程圖如圖2所示。

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

5 算法成像效果

在實驗仿真的過程中，將LBP算法獲得的灰度矩陣作為初始信號矩陣，選擇DCT基作為初始信號的稀疏基對初始信號進行稀疏；將靈敏度矩陣的各行按照高斯隨機順序重新排列，并引入截斷奇異值的思想，增大奇異值中較小的數(shù)值，從而得到新的奇異值矩陣和觀測值，重新建立了基于壓縮感知理論的ECT圖像重建問題。

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

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

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

選取直接算法中的Tikhonov正則化、迭代類算法中的Landweber和基于壓縮感知理論的半閾值迭代類算法進行對比實驗。

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

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

(22)

(23)

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

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

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

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

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

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

圖5 各算法圖像成像時間對比Fig.5 Comparison of imaging time of each algorithm

6 結(jié) 論

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