基于改進(jìn)線性Bregman算法的ECT圖像重建算法

馬 敏，孫美娟

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

1 引 言

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

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

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

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

2 ECT系統(tǒng)模型

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

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

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

SG=C

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

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

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

3 線性Bregman迭代算法

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

(1)

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

3.1 線性Bregman迭代

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

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

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

(2)

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

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

(3)

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

(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]。

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

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

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

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

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

A=UDVT

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

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

由文獻(xiàn)[4]可知，‖A+-V2q-1‖≤‖A+-V0‖‖I-αAA*‖2q-1，又因?yàn)椤琁-αAA*‖<1，所以‖A+-V2q-1‖<‖A+-V0‖。

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

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

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

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

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

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

7)如果滿足任何一個(gè)停止迭代標(biāo)準(zhǔn)，則迭代停止。

相對(duì)誤差：‖uk+1-uk‖/‖uk‖<ε

絕對(duì)誤差：‖uk+1-uk‖<ε

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

8)返回輸出uk+1。

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

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

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

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

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

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

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

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

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

表4 圖像相對(duì)誤差Tab.4 Comparison of image relative error

6 結(jié) 論

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