馬 敏，孫美娟

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

1 引 言

電學(xué)層析成像(electrical tomography，ET)是20世紀80年代被提出并逐漸發(fā)展起來的一種多相流檢測技術(shù)。它根據(jù)研究不同的電學(xué)特性，劃分為電容層析成像(electrical capacitance tomography，ECT)、電阻層析成像(electrical resistance tomography，ERT)和電磁層析成像(electromagnetic tomography，EMT)。其中，電容層析成像因其在檢測過程中具有非侵入性、無害性、成像速度快等特點，應(yīng)用前景十分廣泛[1]。目前主要應(yīng)用于實時監(jiān)測多相流的流型、工業(yè)管道內(nèi)流體情況及航空器中的滑油等[2]。

電容層析成像主要包括傳感器陣列、數(shù)據(jù)采集單元和圖像重建三部分[3]。圖像重建作為ECT系統(tǒng)的主要組成部分，其重建速度和質(zhì)量對整個系統(tǒng)有著重要影響。目前ECT圖像重建算法可分為非迭代類算法和迭代類算法。非迭代類算法主要有線性反投影，Tikhonov正則化、奇異值分解法(singular value decomposition，SVD)等，其共同點是一步成像，圖像重建速度快，但重建精度不高。迭代算法主要有Landweber、共軛梯度(conjugate gradient，CG)等，其重建圖像分辨率較好，但收斂速度較慢，實時性不佳。

Donoho和其他人引入了Bregman迭代并應(yīng)用于解決基追蹤問題[4,5]。由Osher S，Burger M，Goldfarb D，Xu J，Yin W等給出的線性Bregman迭代算法(linearized Bregman algorithm，LBA)主要用來求解稀疏問題，應(yīng)用于圖像去噪領(lǐng)域，之后又被用于基于壓縮感知的圖像重建問題，是一種快速、有效的優(yōu)化算法。

本文將線性Bregman迭代應(yīng)用到ECT圖像重建當中，為進一步提高成像質(zhì)量，減少成像時間，提出兩種改進線性Bregman算法(modifified linearized Bregman algorithm，MLBA)，即通過奇異值分解和牛頓二階迭代逼近廣義逆，與線性Bregman迭代結(jié)合形成的兩種混合迭代算法。然后利用圖像重建速度、相對誤差、相關(guān)系數(shù)等評價指標對常用的成像算法成像結(jié)果與改進后的算法成像結(jié)果進行對比分析。

2 ECT系統(tǒng)模型

電容層析成像基本工作原理是通過安裝在被測管道外側(cè)的敏感電容傳感器，對管道內(nèi)部介質(zhì)進行測試并獲取相應(yīng)的介質(zhì)信息，再借助于計算機圖像重建算法，獲取管道內(nèi)介質(zhì)分布圖像。在實際操作過程中，由于不同流體的介電常數(shù)不同，當流體通過傳感器陣列時會引起極板間介電常數(shù)的改變，從而引起電容值的變化，采集單元進行電容值數(shù)據(jù)采集，然后將數(shù)據(jù)處理后傳送到成像計算機，成像計算機再選擇合適的成像算法完成圖像重建[6]。其系統(tǒng)模型如圖1所示。

圖1 ECT系統(tǒng)模型Fig.1 System model of ECT

ECT系統(tǒng)圖像重建的線性模型如下：

SG=C

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

影響圖像重建的主要因素有：(1)軟場特性：敏感場分布易受被測介質(zhì)分布的影響；(2)欠定性：反問題求解過程中，獨立測量的電容值個數(shù)遠小于未知量的個數(shù)，導(dǎo)致解不唯一；(3)不適定性：邊界測量的變化對場域內(nèi)介質(zhì)變化不敏感，邊界電位值的微小變化會引起解的較大變化，導(dǎo)致求解過程不穩(wěn)定[7]。

在實際工程中，ECT技術(shù)中的非線性問題一般采用局部線性化的方法解決，因此，它的圖像重建問題轉(zhuǎn)變?yōu)橐环N優(yōu)化問題。由于數(shù)據(jù)采集條件的限制，往往導(dǎo)致所采集的數(shù)據(jù)量不能滿足恢復(fù)原始圖像的要求，此時相應(yīng)的重建方程SG=C就變成了一組欠定方程[8]。

3 線性Bregman迭代算法

在所有解決欠定線性系統(tǒng)Au=f的方案中，通過最小化l1范數(shù)‖u‖1來恢復(fù)u，這個問題被稱為最小化問題：

(1)

線性Bregman算法在解決了Bregman每一步需要最小化的問題，實質(zhì)上是一階導(dǎo)數(shù)的線性逼近[9]。因此，非常適用于求解l1范數(shù)正則化問題，也符合ECT圖像重建的求解問題。

3.1 線性Bregman迭代

定義3.1(Bregman距離)：凸函數(shù)J在u，v兩點的Bregman距離D(u,v)定義如下：

式中：p∈?J(v)是凸函數(shù)J在V點的次梯度。

利用Bregman距離的概念，線性Bregman迭代旨在解決式(1)中的優(yōu)化問題:

(2)

式中：J(u)=‖u‖1；δ是常量；p0=u0=0。

由式(2)產(chǎn)生的序列{uk}k∈N的極限是式(3)的唯一解：

(3)

盡管式(3)與式(1)不同，但是當μ→∞時，式(3)趨近于式(1)，這個迭代方程可以寫為如下形式：

(4)

這里u0=v0=0，軟閾值算子

Tλ(w)=[tλ(w1),tλ(w2),…,tλ(wn)]T

(5)

式中：

(6)

在此基礎(chǔ)上，把迭代規(guī)則式(4)推廣到式(7)

(7)

式中：A∈Cm×n,m≤n，是任意矩陣，0<δ<2/‖AAT‖[10]。

且文獻[10]中證明了當μ→∞，0<δ<1時，式(7)中序列極限趨向于式(1)的解，并且在所有解中最接近Au=f的最小l2范數(shù)解。

4 基于MLBA的ECT圖像重建算法

盡管線性Bregman可以解決式(1)問題，但是由于線性Bregman的收斂速度與矩陣A的條件數(shù)有關(guān)，條件數(shù)越小，其收斂速度越快；式(7)中涉及求解A的廣義逆矩陣，目前MATLAB中存在很多通過迭代方式逼近廣義逆的方法[11]，但大部分都是一階迭代，其收斂速度不快且誤差較大。對此，本文通過減小A的條件數(shù)、提高求A+的迭代階數(shù)來改進線性Bregman算法。

4.1 基于奇異值分解的線性Bregman重建算法

定義4.1(奇異值分解)：對于m×n的矩陣A(m

A=UDVT

式中：U是m×m階酋矩陣；V是n×n階酋矩陣；D是m×n的半正定對角矩陣，它的前r列的對角線元素包含了S的r個奇異值，即為矩陣S的全部非零奇異值，分別記作σ1,σ2,…,σr(σ1≥σ2≥…≥σr>0)，其中r是矩陣S的秩，其對于矩陣的影響也隨著數(shù)值遞減相應(yīng)地減弱，其余對角元素為0。另外，U的列向量是AAT的特征向量，V的列向量是ATA的列向量。所以，A+=VD-1UT。

雖然SVD分解可以減小A的條件數(shù)，加速線性Bregman算法收斂速度[12]，但SVD分解耗費時長較長，且ECT系統(tǒng)中靈敏度矩陣的奇異值逐漸趨向于零。對此，有學(xué)者提出在D-1的前r個主對角元素加入修正因子的方法來改善A+的穩(wěn)定性[13]，本文用此SVD思想來改進線性Bregman迭代，從而達到既提高收斂速度，也增加穩(wěn)定性的目的。形成的改進算法1(MLBA1)如下：

(8)

式中：u0=0；f0=0；0<σ<1；kmax=20；k=0,1,2,…。

4.2 基于二階迭代的線性Bregman重建算法

引理4.1：設(shè)給定矩陣A∈Am×n≠0，初始矩陣V0∈An×m滿足[14]

(1)V0∈μ(A*,A*)

(2)ρ(I-AV0)<1

式中：I為m×m單位矩陣；ρ(A)是矩陣A的譜半徑。則(9)產(chǎn)生的序列{Vq}q∈N收斂于A+。

Vq+1=Vq+V0(I-AVq),q=0,1,…

(9)

由于逼近引理4.1中的廣義逆是一階迭代，令：

Wq+1=Wq+Wq(I-AWq),q=0,1,…

(10)

式中：W0可以是滿足ρ(AW0-AA+)<1的任意矩陣，這里取W0=V0=αAT；迭代序列Vq和Wq之間的關(guān)系是Wq=V2q-1，從而把一階迭代變?yōu)槎A迭代。

由文獻[4]可知，‖A+-V2q-1‖≤‖A+-V0‖‖I-αAA*‖2q-1，又因為‖I-αAA*‖<1，所以‖A+-V2q-1‖<‖A+-V0‖。

因此，V2q-1比V0更接近A+，二階迭代將比一階迭代提供更多關(guān)于A+的信息，且迭代收斂速度增加，從而不僅能提高圖像重建速度，也提升了重建質(zhì)量。用二階迭代改進線性Bregman迭代，形成改進算法2(MLBA2)如下：

(11)

式中：u0=0；ξ0=Wqf0；f0=0，0<δ<1；0<α<1/‖A‖2；W0=αAT。基于MLBA2的求解步驟如下：

1)初始化：

k=0，u0=0，ξ0=Wqf0，f0=0，α=1/2‖A‖2，δ=0.3，kmax=100，W0=αAT，ε=0.18，μ=0.5

2)先計算Wq：

Wq+1=Wq(2I-AWq),q=0,1,2

3)執(zhí)行步驟4～7，直到滿足停止標準。

4)計算：fk+1=fk+(f-Auk)。

5)計算：ξk+1=ξk+Wq(fk+1-Aξk)。

6)計算：uk+1=δTμ(ξk+1)

7)如果滿足任何一個停止迭代標準，則迭代停止。

相對誤差：‖uk+1-uk‖/‖uk‖<ε

絕對誤差：‖uk+1-uk‖<ε

最大迭代次數(shù)：k=kmax

8)返回輸出uk+1。

5 仿真及結(jié)果分析

為驗證本文提出的兩種MLBA算法圖像重構(gòu)效果，實驗選取12電機的仿真模型，管道外徑50 mm、內(nèi)徑46 mm，對氣固兩相流進行建模分析，各類介質(zhì)介電常數(shù)分別?。嚎諝?,銅2.2,塑料5.8,玻璃4.2。

依次對核心流、三泡流、四泡流、環(huán)流和層流5種流型進行仿真，仿真原型見表1中第1行所示。使用COMSOL Multiphysics 5.3軟件建立模型，采用網(wǎng)格自動剖分法，對正問題進行求解；然后通過Matlab 2016a對反問題求解，進行圖像重建仿真實驗。實驗選取Landweber、CG、SVD、線性Bregman迭代算法與MLBA算法進行仿真對比，各個參數(shù)取值為：q=2[15]，μ=5[16]，δ=0.3，kmax=100，ε=18，μ=0.5圖像評價標準采用圖像重建速度、相關(guān)系數(shù)和相關(guān)誤差[17]。

各種算法圖像重建效果如表1所示，圖像重建速度如表2所示，相關(guān)系數(shù)和相對誤差如表3、表4所示。

由表1圖像可知，LBA可以用于ECT圖像重建，但僅在核心流中成像效果好，對其他復(fù)雜流型成像效果很差。MLBA1和MLBA2算法在成像質(zhì)量上相對于Landweber、CG、SVD、LBA算法有明顯提升，其中，核心流的偽影得到改善；三泡流和四泡流減少了圖像粘連情況，目標位置、輪廓較為清晰；環(huán)流不僅形狀得到改善，而且偽影減少。且MLBA1相比MLBA2算法圖像偽影更小，成像更清晰。雖然MLBA1成像效果最好，但由表2可知，MLBA1在時間上相比Landweber、CG、SVD并無很大優(yōu)勢，且圖像越復(fù)雜比其他算法耗費時間越多。而MLBA2算法，在成像速度上有很大提升，只需Landweber算法成像時間的近十分之一，MLBA1成像時間的約20分之一。因此從系統(tǒng)實時性和重建質(zhì)量兩方面綜合考慮，MLBA2算法的優(yōu)勢更明顯。

表1 不同算法重建效果Tab.1 Reconstruction effects of different algorithm

表2 圖像重建速度比較Tab.2 Comparison of image reconstruction speed s

從表3和表4的誤差指標數(shù)據(jù)可知，MLBA1與MLBA2的圖像誤差得到了很好的控制，圖像相對系數(shù)也相對得到了提高。故改進后的算法提高了圖像成像質(zhì)量。

表3 圖像相關(guān)系數(shù)Tab.3 Comparison of image correlation coefficient

表4 圖像相對誤差Tab.4 Comparison of image relative error

6 結(jié) 論

本文在分析ECT欠定性問題的基礎(chǔ)上，將線性Bregman迭代算法應(yīng)用到ECT圖像重建過程，并對該算法進行了改進，使之不僅具有線性Bregman算法的優(yōu)點，而且可以加快計算速度，提高準確性。仿真實驗表明，基于線性Bregman迭代算法能夠?qū)崿F(xiàn)ECT圖像的重構(gòu)，但對復(fù)雜流型成像質(zhì)量不高。本文提出的兩種改進算法，相比Landwebr、CG、SVD、線性Bregman算法，具有更好的圖像分辨率和更快的成像速度。但是，由于仿真的流型有限且線性Bregman迭代中涉及廣義逆矩陣的求解，后續(xù)研究可以繼續(xù)優(yōu)化廣義逆的求解方法，使該算法在電容層析成像技術(shù)中更有優(yōu)勢、效果更佳。