張玲玲 ，王化祥，范文茹

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

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

1 ERT圖像重建

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

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

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

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

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

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

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

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

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

1.3 ERT圖像重建

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

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

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

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

2 算法分析

2.1 最小二乘法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

從而式(14)變換為

對式(15)中的限制條件 ? xi≤ ui≤ xi定義對數(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；迭代的步長采用回追步長搜索方法．

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

3 仿真實驗結(jié)果

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

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

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

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

4 結(jié) 語

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

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

［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]. 天津大學(xué)學(xué)報，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］ 畢德顯. 電磁場理論[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ù)測量 ERT圖像重建算法的研究[J]. 浙江大學(xué)學(xué)報：工學(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.