馬 敏，孫美娟，李 明

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

1 引 言

電容層析成像(electrical capacitance tomography， ECT)是上個世紀80年代發(fā)展起來的一種用于工業(yè)管道檢測的過程層析成像技術。由于其快速、安全、非侵入、無損害的特點已被應用于氣液兩相流空隙率測量、流型識別、氣固兩相流濃度分布測量、航空發(fā)動機滑油檢測等多個領域[1]。

ECT傳感器系統(tǒng)主要有傳感器陣列、數(shù)據(jù)采集單元和成像計算機3部分。在實際操作過程中，由于不同流體的介電常數(shù)不同，當流體通過傳感器陣列時會引起極板間介電常數(shù)的改變，從而測量引起電容值的變化，采集單元進行電容值數(shù)據(jù)采集，然后將數(shù)據(jù)處理后傳送到成像計算機。成像計算機再選擇合適的成像算法完成圖像重建[2]。

目前ECT圖像重建算法可分為非迭代算法和迭代算法。非迭代算法主要有線性反投影算法、Tikhonov正則化算法等，其共同點是一步成像，圖像重建速度快，但是重建精度不高。迭代算法主要有Landweber算法、共軛梯度算法等，其重建圖像分辨率較好，但收斂速度較慢，實時性不佳。

2006年，Candes及其所在團隊提出壓縮感知的概念。該理論指出如果信號本身是稀疏的或者可壓縮時，那么可以通過線性測量，對信號進行低維表示，因此可以獲得比奈奎思特采樣更少的采樣數(shù)，并且能實現(xiàn)對原始信號的精確重構。該理論極大地降低了數(shù)據(jù)采集、存儲及傳輸過程的壓力，同時也極大地提高了信號處理效率[3]。

2 壓縮感知基本原理

2.1 信號的稀疏表示

壓縮感知理論認為只要信號是稀疏或者可壓縮的，就可以經(jīng)過某種線性觀測，再利用重構算法從少量觀測值中恢復出原始信號，不同于傳統(tǒng)的信號采樣，基于壓縮感知理論的的信號處理方式將壓縮與采樣同時進行，采用這種方式會丟棄大量冗余或無用的信息，不僅提升了采集速度，更減少了信號存儲和傳輸過程的壓力。壓縮感知信號處理過程如圖1所示。

圖1 壓縮感知理論框架Fig.1 Theoretical framework of compressed sensing

根據(jù)壓縮感知理論，壓縮測量實際上是線性測量，通過觀測矩陣對信號進行觀測，因此觀測矩陣的設計需要保證原始信息盡可能多的包含在觀測值當中，這樣才能保證原始信號由觀測值精確地重構出來[4]。

假設長度為N的一維信號x∈RN,由信號稀疏理論可知，該信號可用一組正交向量基ψ=[ψ1,ψ2,ψ3,…,ψN]線性表示為[5]：

(1)

式中:ψi(i=1,2,3,…,N)為N維列向量；α為長度為N的一維向量；α中大部分系數(shù)為0或者接近0，且僅有K個非零元素，滿足K?N。此時信號x在正交基ψ下是稀疏的，稀疏系數(shù)為α，稀疏度為K。在壓縮感知理論中，只要信號具有稀疏性或者在經(jīng)過某種變換后變得稀疏，就可以把信號進行壓縮。

常用的稀疏基有:離散Fourier變換基(DFT)、離散余弦變換基(DCT)、離散正弦變換基(DST)、離散小波變換基(DWT)等。

2.2 觀測矩陣的設計

基于壓縮感知理論的信號處理是通過觀測矩陣進行信號的壓縮與采樣，因此需要設計一個隨機觀測矩陣φ∈RM×N(M?N)，觀測矩陣和稀疏基不相關，信號的投影過程如下：

y=φx

(2)

式中:y∈RM×1,可以看出投影得到的向量維數(shù)遠小于原始信號。將式(1)代入式(2)，整個壓縮感知過程可表示為：

y=φx=φψα=Acsα

(3)

式中：y為M×1維觀測值向量；φ為M×N維觀測矩陣;令Acs=φφ，稱Acs為壓縮感知矩陣。

對于給定的y從式(3)中求出x是一個線性規(guī)劃問題，但由于M?N，使得方程個數(shù)少于未知數(shù)的個數(shù)，故此方程為不定方程，通常無確定解。但如果原始信號x在某組基下的系數(shù)向量是稀疏的，則可以運用壓縮感知理論求其確定解。文獻[3]指出，當壓縮感知矩陣Acs滿足限制等距條件(restricted isometry property，簡稱RIP)時，系數(shù)向量α的K個非零系數(shù)便可以從M個測量值中準確恢復，從而保證重建算法收斂。然而判斷給定矩陣Acs是否滿足RIP性質是一個較為復雜的過程[6]。文獻[7]指出，當觀測矩陣和稀疏基不相關時，Acs可以在很大概率上滿足RIP條件。

2.3 重構算法

最后就是信號的重構過程，對于式(3),矩陣Acs中列數(shù)大于行數(shù)，因此該線性方程有無數(shù)多個解[8]，為了滿足壓縮感知理論信號稀疏性的前提，因此系數(shù)向量α應盡可能稀疏，該問題就變成一個系數(shù)向量盡可能稀疏的最優(yōu)問題。理論證明，在測量矩陣滿足限制等距條件時，原始信號可以通過下式優(yōu)化問題進行重構[9]：

(4)

由于l0范數(shù)具有非凸性，使得直接求解式(4)的數(shù)值計算極不穩(wěn)定而且是NP困難問題。Donoho和Chen等人指出，當觀測矩陣和稀疏基不相關時，求解一個更加簡單的最小l1范數(shù)問題可以得到同等的解：

(5)

基于l1-范數(shù)的算法雖然重構效果好但計算復雜度高而且傳輸數(shù)據(jù)有較大冗余。文獻[10]中進一步改進重構算法的性能，使壓縮感知重構算法能更加具有一般性，提出基于lp(0

(6)

3 基于自適應閾值迭代的ECT圖像重建

3.1 稀疏基的選取

CS理論要求觀測系統(tǒng)對原始信號的觀測過程為線性過程，ECT系統(tǒng)線性模型如下：

C=SG

(7)

式中:C為M×1維歸一化電容測量值矩陣；G為N×1維的歸一化介電常數(shù)分布矩陣，在圖像重建過程中代表圖像灰度值；S為M×N維矩陣，反映了電容C受物質分布變化G的影響，稱為敏感場。

圖2 離散余弦變換基稀疏效果Fig.2 Sparse effect of discrete cosine transform basis

在ECT系統(tǒng)中，原始信號為圖像灰度信號，不具有稀疏性，根據(jù)壓縮感知理論，必須對圖像灰度信號進行稀疏化處理，使其滿足壓縮感知理論的前提，本文以典型的四泡流為例，分別使用常用的離散傅里葉變換基、余弦變換基、小波變換基、正弦變換基進行處理,結果表明離散余弦變換基稀疏效果較好，如圖2所示。因此將離散余弦變換基作為ECT系統(tǒng)的稀疏基。

x=ψDCTα

(8)

式中:x為N×1維圖像灰度向量；ψDCT為N×N維離散余弦變換基；α為N×1的稀疏系數(shù)向量。

3.2 基于SVD的觀測矩陣優(yōu)化

在ECT系統(tǒng)中，重建圖像信號可以經(jīng)過稀疏處理變得稀疏，但靈敏度矩陣S是否滿足RIP難以判定。Donoho指出，如果觀測矩陣的列向量獨立性越強，則觀測矩陣與稀疏基之間的非相關性就越強，從而壓縮感知重構的效果越好[10];Szarek等證明了矩陣的列相關性與矩陣的最小奇異值相關，矩陣的最小奇異值越大，列向量間的獨立性越強[11]。因此可以通過增大矩陣的最小奇異值來增加觀測矩陣列向量之間的獨立性，從而對觀測矩陣進行優(yōu)化。

SVD分解是將矩陣A分解為左奇異矩陣U、對角矩陣、右奇異矩陣V，表達式如下：

A=U[Δ0]VT

式中:A為M×N維矩陣；U為M×M維正交矩陣；[Δ0]為M×N維對角矩陣；V為N×N維正交矩陣；Δ=diag(δ1,δ2,…,δM),δ1,δ2,…,δM為矩陣A的奇異值。

通過矩陣的奇異值分解將靈敏度矩陣進行分解，然后計算對角矩陣的奇異值的平均值，并將奇異值都賦值為平均值，處理方式如下，令：

式中:δmean為δ1,δ2,δ3…,δM的均值。

Snew=QUTS=[Δ*0]VT

(10)

式中:S為未優(yōu)化前的靈敏度矩陣；

相應的電容測量值也變換為：

Cnew=QUTC

(11)

變換后可得到新的ECT系統(tǒng)數(shù)學模型：

Cnew=Snewg=SnewψDCTα=Acsα

(12)

式中:Snew為優(yōu)化后的M×N維系統(tǒng)觀測矩陣；Cnew為M×1維的電容值向量；Acs為壓縮感知矩陣。

3.3 自適應閾值迭代算法

近年來，在稀疏信號重建方面，隨著非凸松弛方法的發(fā)展，許多學者證明了lp(0

(13)

式中:ε=(ε1,ε2,ε3,…εn)，通過調(diào)整參數(shù)εi>0，可以得到

(14)

則

(15)

原來的最優(yōu)化模型可以更改為：

(16)

根據(jù)對新模型建立閾值表示，并提出改進的迭代閾值算法進行求解。

對于任意的參數(shù)λ>0，μ>0，α>0，令

易知，對于任意的μ>0，H2(x,x)=H1(x)

證明過程如下：

≥μH1(x*)

=H2(x*,x*)

(17)

式中:Bμ(x*)=x*+μAT(b-Ax*)。

證明過程如下：

由引理可得，令x*∈Rn為式(17)的解，則：

(18)

式中:[Bμ(x*)]i為向量Bμ(x*)的第i項。

結合閾值表示理論，用于求解式(17)的迭代閾值算法可以被定義為：

=sign([Bμ(xk)]i)·

(19)

通常對于正則化模型，圖像的重建質量取決于正則化參數(shù)的選取，然而人工選擇合適的參數(shù)非常困難，因此本文選用自適應的方法進行參數(shù)的選取，令x*為模型的最優(yōu)解，稀疏度為r,假設

|[Bμ(x*)]1|≥|[Bμ(x*)]2|≥

…≥|[Bμ(x*)]n|，

根據(jù)式(19)有以下不等式：

可得：

則令：

(20)

在本文中，通過自適應的方法來選擇最優(yōu)參數(shù)，使得p值能夠根據(jù)輸入?yún)?shù)自動確定其大小，為了使p值收斂敏感定義相對誤差:

則:

rRMSEk=rRMSEk-rRMSEk-1

最終p值的更新為：

(21)

式中：β為步長，通過這種方法可以有效求解式(16)。

4 仿真實驗

為了驗證范數(shù)值p取不同時對重構率的影響，固定的觀測矩陣A∈R120×3 228,隨機生成稀疏向量x0，通過Ax0=b得到向量b，通過本文算法求解原始信號x0，實驗取Tol=10-5，β=1。由于該算法對參數(shù)ε的依賴性較強，這里設置εi=max{0.7|[μAT(b-Axk)]i|,10-3}。定義重構率為原始信號和重構信號之間的相關系數(shù)。同時為驗證該算法在ECT圖像重建過程中p值的選取對重建結果的影響，選取三泡流、四泡流兩種流型進行不同p值時的圖像重建。

圖3 不同p值及稀疏度時的重構率Fig.2 Reconstruction rate with different p values and sparsity

由圖3所示，顯示了不同p值時，本文算法對不同稀疏度信號的重構能力，比較不同p值時算法的重構能力，可知當稀疏度為[0,40]時，p=0.7的重構效果較好，當稀疏度為[40，80]時，p=0.1的重構效果較好，當稀疏度為[80,95]時，p=0.9的重構效果較好，由表1可知，四泡流中，p=0.7時的重建效果偽影較少圖像較清晰，而三泡流中，p=0.3時圖像質量較好，目標位置輪廓準確。綜上可知p值的選取影響著最終的重構效果。

表1不同p值重建結果
Tab.1Reconstruction results with different value

為驗證本文提出的自適應與之迭代算法的重構效果，本組實驗選取三泡流、十字流、環(huán)流、層流4個相對復雜的流型進行仿真實驗，選取Landweber算法、IRLS算法進行仿真實驗對比，重建結果如表2所示，圖像重建速度、圖像相對誤差、圖像相關系數(shù)如表3、表4、表5所示。

表2 不同算法重建結果Tab.2 Reconstruction results with different algorithm

表3 圖像重建速度比較Tab.3 Image reconstruction speed comparison s

表4 圖像相對誤差Tab.4 Image relative error

表5 圖像相關系數(shù)Tab.5 Image correlation coefficient

由仿真實驗結果可以看出，本文提出的自適應閾值迭代算法在成像質量上，相對于Landweber迭代算法有明顯提升，三泡流粘連情況得到改善，目標位置、輪廓較為清晰;十字流成像結果雖然未完全反映原始流型信息，但是改進算法可以明顯辨別出目標的形狀;改善最為明顯的是環(huán)流，邊緣效應的影響減少，且偽影較少，目標輪廓明顯,從數(shù)據(jù)反應看，圖像誤差得到了很好的控制。在成像速度方面，由于算法引入了閾值表示理論，在計算量方面得到了很好的改善，成像速度也得到了明顯提高。

5 結 論

(1) 本文在分析ECT不定問題的基礎上，提出將壓縮感知理論應用到ECT的圖像重建中去。仿真實驗表明，基于壓縮感知的ECT圖像重建算法能夠利用少量數(shù)據(jù)實現(xiàn)ECT圖像的精確重構。

(2) 針對取不同值時，對隨機生成數(shù)據(jù)及三泡流、四泡流的重構結果的成像質量進行比較。實驗表明，基于lp(0

(3) 本文提出的基于lp(0