基于1范數(shù)的電阻層析成像圖像重建算法

張玲玲 ，王化祥，范文茹

(1. 天津大學(xué)理學(xué)院，天津 300072；2. 天津大學(xué)電氣與自動(dòng)化工程學(xué)院，天津 300072)

電阻層析成像(electrical resistance tomography，ERT)可對(duì)封閉的管道或過程容器設(shè)備內(nèi)部多相介質(zhì)進(jìn)行可視化測(cè)量，非常適用于以液相為連續(xù)相的多相流檢測(cè)[1-3]．圖像重建是ERT系統(tǒng)的關(guān)鍵問題之一．傳統(tǒng)的圖像重構(gòu)方法，如共軛梯度法、Tikhonov正則化方法等均是基于2范數(shù)的優(yōu)化方法．2范數(shù)處理適用于具有一定光滑性的信號(hào)重建，但不能有效地表征信號(hào)的稀疏性．ERT圖像重建的信號(hào)通常具有很強(qiáng)的稀疏性，基于 2范數(shù)重建圖像邊界模糊，圖像質(zhì)量不高．為能更好地反映ERT重建信號(hào)的稀疏性，獲得相對(duì)理想的重建效果，筆者采用基于 1范數(shù)的1l正則化最小二乘法(1l-regularized least-squares programs，LSPs)進(jìn)行ERT圖像重建，引入?yún)?shù)λ的自適應(yīng)選擇標(biāo)準(zhǔn)，將其與Tikhonov正則化方法進(jìn)行比較，并通過仿真實(shí)驗(yàn)驗(yàn)證其有效性．

1 ERT圖像重建

1.1 ERT系統(tǒng)工作模式

筆者的研究組開發(fā)的 ERT系統(tǒng)， 設(shè)計(jì) 16敏感電極陣列，采用相鄰激勵(lì)模式，如圖 1所示．敏感陣列獲取被測(cè)物場(chǎng)分布的投影信息，測(cè)量和數(shù)據(jù)收集所獲得的物場(chǎng)信息將傳給主控計(jì)算機(jī)，進(jìn)行圖像重建與顯示，從而實(shí)現(xiàn)對(duì)被測(cè)物場(chǎng)的監(jiān)測(cè)和控制[3-4]．

圖1 ERT系統(tǒng)結(jié)構(gòu)Fig.1 ERT system structure

1.2 ERT問題的數(shù)學(xué)模型

ERT技術(shù)是根據(jù)敏感場(chǎng)的電導(dǎo)率分布獲得的物場(chǎng)媒質(zhì)分布信息．在敏感長(zhǎng)邊界施加激勵(lì)電流，當(dāng)場(chǎng)內(nèi)電導(dǎo)率分布變換時(shí)，導(dǎo)致場(chǎng)內(nèi)電勢(shì)分布變換，從而使場(chǎng)域邊界上的測(cè)量電壓發(fā)生變化．通過一定的圖像重建算法可重建出場(chǎng)內(nèi)的電導(dǎo)率分布．

根據(jù)似穩(wěn)場(chǎng)理論的麥克斯韋方程，對(duì)于 ERT敏感場(chǎng)內(nèi)任意一點(diǎn)滿足[5]

式中：S為電流密度；σ為電導(dǎo)率；E為電場(chǎng)強(qiáng)度；φ為場(chǎng)內(nèi)電勢(shì)分布.由式(1)～式(3)，φ滿足

ERT滿足第2類邊界條件，其形式可表示為

式中l(wèi)I為電極l的注入電流．

1.3 ERT圖像重建

ERT實(shí)現(xiàn)一般基于數(shù)值方法，分正問題和逆問題．正問題是已知敏感場(chǎng)內(nèi)被測(cè)介質(zhì)的電導(dǎo)率分布和敏感場(chǎng)的邊界條件(外部激勵(lì)電流)，求解電磁場(chǎng)的電勢(shì)分布．通過有限元仿真軟件(COMSOL)獲得正問題的解．對(duì)正問題的求解可獲取場(chǎng)內(nèi)所有節(jié)點(diǎn)的電位信息，根據(jù)敏感電極邊界的測(cè)量信息，求解敏感場(chǎng)內(nèi)介質(zhì)的電導(dǎo)率分布并重建出物場(chǎng)分布圖像，即為ERT逆問題．逆問題由于測(cè)量數(shù)據(jù)少，存在嚴(yán)重的病態(tài)性，測(cè)量數(shù)據(jù)的微小擾動(dòng)可導(dǎo)致解的劇烈變化甚至發(fā)散．因此，ERT問題的解決主要依賴于逆問題求解，即圖像重建算法的計(jì)算精度和速度．

通常，電學(xué)系統(tǒng)方程進(jìn)行線性化處理，ERT系統(tǒng)對(duì)應(yīng)的方程為

式中：J是Jacobian矩陣；δU為隨σ變化的邊界電極間電壓差：δσ為內(nèi)部場(chǎng)域的電導(dǎo)率的變化值．為方便起見，式(7)轉(zhuǎn)化為線性系統(tǒng)方程

由于矩陣 A 通常為大型稀疏的欠定矩陣，從而式(8)的解不存在或者不唯一．

2 算法分析

2.1 最小二乘法

求解病態(tài)方程組 =Ax y比較有效的算法是最小二乘法．對(duì)于 ERT系統(tǒng)來講，靈敏度矩陣 A是欠定矩陣，方程組的解不存在或者不唯一，通過求最小范數(shù)最小二乘解作為其最優(yōu)解，即

由于矩陣A不對(duì)稱，求廣義逆過程較復(fù)雜． 為簡(jiǎn)化求解，將欠定矩陣對(duì)應(yīng)的線性方程組求解問題轉(zhuǎn)化為對(duì)稱矩陣對(duì)應(yīng)的線性方程組求解，即將問題(9)轉(zhuǎn)化為

求解的解析方法由下面的定理給出．

定理 1[6]對(duì)任意非零矩陣A ∈，方程組(10)解存在且 x = A+b為方程組的解．這里 A+是矩陣 A的廣義逆．

問題轉(zhuǎn)化成凸規(guī)劃問題，通過相應(yīng)的優(yōu)化算法實(shí)現(xiàn)解此問題常用的算法有共軛梯度法和Landweber迭代法等[7-8]．

由于ERT系統(tǒng)的靈敏度矩陣A具有很強(qiáng)的稀疏性和病態(tài)性，對(duì)稱矩陣TAA的譜趨近于零，導(dǎo)致問題的嚴(yán)重不適定性．為了克服病態(tài)性對(duì)問題求解精度的影響，一般采用 Tikhonov正則化方法[9]，其優(yōu)化算法的形式為

式中λ為正則化參數(shù)．

式(12)的解析解表達(dá)式為

當(dāng)λ是 ATA的正則點(diǎn)時(shí)，ATA?λI可逆．λ選擇適當(dāng)時(shí)，此方法可有效求解方程(9)的不適定性，但參數(shù)λ的選取需要憑經(jīng)驗(yàn)選定．式(12)是用二次規(guī)劃的優(yōu)化算法得到式(13)的數(shù)值解，計(jì)算速度快，實(shí)時(shí)性強(qiáng)，對(duì)于ERT圖像重建問題具有一定的效果．

以上算法是基于2l范數(shù)誤差估計(jì)的計(jì)算方法，不可避免地會(huì)將原問題過度光滑，從而使不光滑的信息丟失，空間的分辨率不高．對(duì)于 ERT應(yīng)用的大多數(shù)領(lǐng)域，被測(cè)區(qū)域的介質(zhì)往往存在不連續(xù)性，靈敏度矩陣 A具有很強(qiáng)的稀疏性．基于 2范數(shù)的最小二乘法無法反映 ERT系統(tǒng)的稀疏性和不連續(xù)性，為了改善圖像重建質(zhì)量，借鑒Tikhonov正則化思想，改變誤差估計(jì)方式，采用1l正則化的最小二乘法進(jìn)行圖像重建．

2.2 LSPs方法及其實(shí)現(xiàn)

對(duì)于方程(8)，當(dāng)矩陣 A是稀疏的或其解具有可壓縮性時(shí)，將方程(8)的求解轉(zhuǎn)化[10]為

式(14)與式(12)不同之處在于罰函數(shù)用 1范數(shù)代替 2范數(shù)．LSPs方法是基于壓縮傳感理論發(fā)展起來的，特別針對(duì)大型稀疏線性方程組求解的新方法．既能克服問題的不適定性，同時(shí)能反映矩陣的稀疏性，避免2范數(shù)將問題過度平滑的缺點(diǎn)[11-12]．

內(nèi)點(diǎn)法是目前求解 1范數(shù)優(yōu)化問題最精確的算法，該方法由加州理工的壓縮傳感研究組提出[13]．整個(gè)過程是把 1 范數(shù)轉(zhuǎn)化為不等式限制，再把不等式轉(zhuǎn)化為門限函數(shù)進(jìn)而去掉不等式．然后用迭代法進(jìn)行求解．

為了去掉式(14)中的 1范數(shù)，引進(jìn)新的變量iu，使其滿足

從而式(14)變換為

對(duì)式(15)中的限制條件 ? xi≤ ui≤ xi定義對(duì)數(shù)門限函數(shù)(logarithmic barrier function)為

設(shè)置門限函數(shù)是模擬不等式的限制，從而有約束條件的二次規(guī)劃問題式(15)變成無約束二次規(guī)劃問題，即

參數(shù)t的變化范圍為 0～∞．采用牛頓迭代法為基本原理，每次迭代令 t增大，增大的方式是按照序列，…, 其中μ＞1．

搜索方向可通過所表示的Newton系統(tǒng)的精確求解得到，即

式中：H為Hession矩陣；H = ?2φt(x,u ) ∈Rn×2n；迭代的步長(zhǎng)采用回追步長(zhǎng)搜索方法．

該算法中，參數(shù)λ的選擇是問題的關(guān)鍵，決定計(jì)算速度和成像質(zhì)量．從前面分析可知，λ在之間變化，當(dāng)0λ→時(shí)接近問題的解，當(dāng)時(shí)， 0x→ ．如果選擇過小，必然導(dǎo)致計(jì)算時(shí)間過長(zhǎng)，實(shí)時(shí)性較差；如果λ選擇過大，將會(huì)嚴(yán)重影響計(jì)算結(jié)果的準(zhǔn)確度．事實(shí)上，λ取決于問題的稀疏程度．對(duì)于確定的方程(8)，的大小反映了問題的稀疏度，選擇即可．不失一般性，令．如果矩陣的稀疏性較差，為得到更好的效果，令

3 仿真實(shí)驗(yàn)結(jié)果

采用有限元方法模擬不同流型分布驗(yàn)證本文提出的方法．ERT系統(tǒng)采用 16個(gè)電極，圖像采用 812像素剖分．計(jì)算機(jī)配置為 Pentium(R)，2.93 GHz CPU、1G內(nèi)存，以 MATLAB7.6為測(cè)試平臺(tái)．分別用Landweber迭代算法、共軛梯度法、Tikhonov正則化方法和LSPs進(jìn)行圖像重建．表1是4種不同算法的重建效果，表2是4種算法的迭代次數(shù)與耗時(shí)．

由實(shí)驗(yàn)結(jié)果不難看出，Tikhonov正則化方法所獲得的重建圖像形狀接近真實(shí)，但是邊界模糊，且圖像中存在干擾，需要進(jìn)一步借助先驗(yàn)知識(shí)增強(qiáng)重構(gòu)圖像的魯棒性．而LSPs算法中，可清楚分離不同物質(zhì)，沒有干擾信息，控制參數(shù)λ的選擇具有自適應(yīng)性．模型(1)、(2)，參數(shù)取值，模型(3)、(4)參數(shù)取值．LSPs算法迭代次數(shù)相對(duì)Tikhonov正則化方法迭代次數(shù)增多，實(shí)時(shí)性較差．有待進(jìn)一步改進(jìn)算法提高其實(shí)時(shí)性．

表1 圖像重建算法的仿真實(shí)驗(yàn)結(jié)果Tab.1 Simulation experimental results of image reconstruction algorithm

表2 迭代次數(shù)與耗時(shí)Tab.2 Iterations and elapsed time

4 結(jié) 語(yǔ)

針對(duì) ERT逆問題的嚴(yán)重不適定性和信號(hào)的稀疏性，采用基于1l正則化的最小二乘法，將罰函數(shù)項(xiàng)由傳統(tǒng)的2范數(shù)改成1范數(shù)，改善了Tikhonov算法將問題過度平滑的缺點(diǎn)．采用該方法得到的重建圖像，在不同介質(zhì)分割的邊緣處沒有偽影，重建效果理想．

此外，基于1l正則化的最小二乘法中，參數(shù)λ的選擇具有一定的自適應(yīng)性，相對(duì)于較小即可．一般選擇

［1］ York T. Status of electrical tomography in industrial applications[J]. Journal of Electronic Imaging，2001，10(3)：608-619.

［2］ Daily W，Ramirez A. Electrical resistance tomography[J]. Society of Exploration Geophysicists，2004，23(5)：438-442.

［3］ Dong Feng，Xu Yanbin. Application of electrical resistance tomography in two-phase flow measurement [J].Journal of Engineering Thermophysics，2006，27(5)：791-794.

［4］ 王化祥，汪 婧，胡 理，等. ERT／ECT雙模態(tài)敏感陣列電極優(yōu)化設(shè)計(jì)[J]. 天津大學(xué)學(xué)報(bào)，2008，41(8)：911-918.

Wang Huaxiang，Wang Jing，Hu Li，et al. Optimal design of ERT／ECT dual-modality sensing electrode array[J]. Journal of Tianjin University，2008，41(8)：911-918(in Chinese).

［5］ 畢德顯. 電磁場(chǎng)理論[M]. 北京：電子工業(yè)出版社，1985.

Bi Dexian. Theory of Electromagnetics[M]. Beijing：Publishing House of Electronics Industry，1985(in Chinese).

［6］ 熊洪允. 應(yīng)用數(shù)學(xué)基礎(chǔ)[M]. 天津：天津大學(xué)出版社，2005.

Xiong Hongyun. Basis of Applied Mathematics[M].Tianjin：Tianjin University Press，2005(in Chinese).

［7］ 李 英，黃志光，冀海峰，等. 兩相流參數(shù)測(cè)量 ERT圖像重建算法的研究[J]. 浙江大學(xué)學(xué)報(bào)：工學(xué)版，2003，37(4)：382-385.

Li Ying，Huang Zhiguang，Ji Haifeng，et al. Study on image reconstruction algorithm of electrical resistant tomography for two-phase flow measurement[J]. Journal of Zhejiang University：Engineering Science，2003，37(4)：382-385(in Chinese).

［8］ 袁亞湘，孫文瑜. 最優(yōu)化理論與方法[M]. 北京：科學(xué)出版社，2007.

Yuan Yaxiang，Sun Wenyu. Optimization Theory and Method[M]. Beijing：Science Press，2007(in Chinese).

［9］ 劉繼軍. 不適定問題的正則化方法及應(yīng)用[M]. 北京：科學(xué)出版社，2005.

Liu Jijun. Regularization Methods and Application of Illposed Problems[M]. Beijing：Science Press，2005(in Chinese).

［10］ Kim S J，Koh K，Lustig M，et al. An interior-point method for large-scale 1-regularized least squares[J].IEEE Journal of Selected Topics in Signal Processing，2007，1(4)：606-617.

［11］ Candès E，Romberg J，Tao T. Stable signal recovery from incomplete and inaccurate information[J]. Communications on Pure Applied Mathematics，2005，59(8)：1207-1233.

［12］ Donoho D L. For most large underdetermined systems of linear equations，the minimal l1-norm solution is also the sparsest solution[J]. Comunications on Pure and Applied Mathematics，2006，59(6)：907-934.

［13］ Donoho D L. Compressed Sensing[J]. IEEE Transactions on Information Theory，2006，52(4)：1289-1306.