基于三維模型的改進正則化ERT成像算法

李守曉，王化祥，范文茹，張玲玲

基于三維模型的改進正則化ERT成像算法

李守曉1，王化祥1，范文茹2，張玲玲3

(1. 天津大學電氣與自動化工程學院，天津 300072；2. 中國民航大學航空自動化學院，天津 300300；3. 天津大學理學院，天津 300072)

電阻層析成像(ERT)通過對被測場邊界注入電流，測量被測場電壓變化，重建物場內(nèi)電導率．針對ERT成像分辨率低，提出一種基于三維模型的改進Tikhonov迭代電阻成像算法．針對Tikhonov正則化參數(shù)選擇問題，提出基于同倫映射的方法，并利用非線性函數(shù)Sigmoid調(diào)節(jié)正則化參數(shù)，以獲得的圖像灰度值作為Tikhonov迭代法的初始值進行迭代，重建敏感場圖像．仿真及實驗結(jié)果表明，該方法有效地改進了ERT圖像質(zhì)量．

電阻層析成像；圖像重建；正則化；同倫映射

電阻層析成像[1](electrical resistance tomography，ERT)技術(shù)是基于邊界電壓測量值反映物場內(nèi)部電導率分布的變化[2]．由于該技術(shù)具有非侵入、無輻射、響應快、結(jié)構(gòu)簡單以及成本低廉等優(yōu)點，在石油、化工等領(lǐng)域具有廣闊的應用前景[3]．

ERT成像算法包括非迭代算法和迭代算法[4]．非迭代算法包括線性反投影、Tikhonov正則化算法、基于截斷奇異值分解的直接算法；迭代算法包括landweber迭代算法、牛頓-拉夫遜迭代算法以及共軛梯度迭代算法．線性反投影算法雖然簡單快速，但該算法成像精度較低．landweber迭代算法與最速梯度下降法類似，但其收斂性比較差[5]．牛頓-拉夫遜迭代算法雖然收斂快但計算量較大且具有局部收斂的缺點．Tikhonov正則化算法的成像質(zhì)量主要依靠正則化參數(shù)的選擇．最優(yōu)的正則化參數(shù)將得到質(zhì)量高的圖像．但該種算法中正則參數(shù)一般不易確定，通常需要通過實驗觀測確定，而且確定的參數(shù)并非最優(yōu)值.

為了提高正則化方法的效率，減少由于試探正則化參數(shù)所帶來的不便，筆者基于同倫映射，構(gòu)造了同倫映射函數(shù)，將Tikhonov正則化方法引入了一種連續(xù)可調(diào)的正則化參數(shù)修正方法，節(jié)省了時間．隨著迭代的進行，采用固定點正則化參數(shù)修正，用所得圖像灰度值作為Tikhonov迭代法的初始值進行成像，提高成像質(zhì)量．

1 ERT系統(tǒng)基本原理

為電阻層析成像的反問題求解，即通過已經(jīng)測量的電壓值反求被測管道截面的電導率分布情況．電阻層析成像線性模型為[6]

式中：g為電導率分布矢量，在圖像重建中代表圖像灰度值；z為電壓矢量，即電導率變化引起邊界的電壓變化；S為靈敏度矩陣．ERT反問題的求解即由已知的測量電壓z求解未知g的過程．電阻層析成像的不適定性直接影響成像質(zhì)量．運用正則化的方法可以使求解更穩(wěn)定、成像更精確．

2 算 法

2.1Tikhonov正則化方法

Tikhonov方法是ERT中求解反問題最常用的方法之一．Tikhonov正則化的基本思想是通過約束解的范數(shù)以保證解不發(fā)散．其最小化目標函數(shù)為

式中：μ為正則化參數(shù)；L(g)稱為正則化函數(shù)．應用牛頓-拉夫遜迭代法，標準的迭代公式為

式中I為nn×單位陣．將最小化目標函數(shù)變?yōu)?/p>

采用牛頓-拉夫遜迭代可得到改變的Tikhonov正則化迭代形式

2.2同倫映射

設(shè)X和Y是Rn的非空子集，f，p：X→Y是光滑映射．如果對任意α∈[0,1]，(α,x)∈[0,1]×X都成立，H(α,x)=(1-α)f(x)+αp(x)∈Y 則稱光滑映射H:[0,1]×X→Y是f和p之間的一個線性同倫[7]，α為同倫參數(shù)．

根據(jù)定義，若需要求解的原方程為

基于同倫映射，構(gòu)造同倫函數(shù)

式(7)滿足

由式(7)可知，α為同倫參數(shù)，α∈[0,1]．當α= 0時式(7)等價于式(6)．若從式(7)出發(fā)，在求解的過程中通過連續(xù)減少α的值，最終將求得式(6)的解．由于H(α,x)≥0，所以可按最優(yōu)擬合的原則將式(7)的求解轉(zhuǎn)化為無約束最優(yōu)化問題，即

應用正則化方法進行求解

將式(11)應用于牛頓-拉夫遜迭代將得到與式(5)相似的迭代形式，其中α與μ有著相似的關(guān)系．由于同倫映射可跟蹤α的求解，因此可解決正則化參數(shù)μ難以確定的問題．

2.3改進正則化算法

對于正則化參數(shù)的確定，目前已發(fā)展了Arcangeli和誤差極小化等，上述方法的共同點是需要預先估計原始數(shù)據(jù)的誤差水平．但在過程參數(shù)檢測中，被測參數(shù)分布較復雜，先驗知識較難獲取．傳統(tǒng)試探的方法雖然理論上能獲得最優(yōu)的μ值，但計算量大．針對這種情形，本文中提出改進的正則化方法．由同倫參數(shù)與正則化參數(shù)相似性，可以采用同倫路徑追蹤的方法確定正則化參數(shù)μ．同倫路徑追蹤有很多方法，如高階插值預估-Newton校正算法和歐拉預估-Newton校正算法等．因正則化參數(shù)介于[0，1]之間，考慮到同倫映射，同倫參數(shù)需要從1或接近l的數(shù)逐漸逼近0．參考人工神經(jīng)元網(wǎng)絡的非線性作用函數(shù)，即Sigmoid函數(shù)，引入連續(xù)同倫參數(shù)調(diào)整公式，即

式中：N為迭代次數(shù)，當a=1，b=0時為典型的Sigmoid函數(shù)[8]．式(12)中的同倫參數(shù)即是正則化參數(shù)．

對于本文的問題經(jīng)過多次實驗，a取0.5，b為[-10，10]之間的實數(shù)，本文b=-0.5．由于同倫參數(shù)的跟蹤變化是連續(xù)的，可以保證迭代穩(wěn)定地進行．

引入原始圖像與重建圖像的相關(guān)系數(shù)[9]

式中：g*是被測區(qū)域內(nèi)的真實電導率分布；g是電導率計算值；g*和g分別是g*和g的平均值．由于式(12)隨著N值的越來越大，同倫參數(shù)的變化很緩慢，算法效率低．故在圖像最大相關(guān)系數(shù)下降5%時，μ值取0.001，將由式(13)的算法所成圖像作為式(5)迭代圖像灰度初始值，提高算法的實時性[10]和圖像質(zhì)量．

改進正則化方法采用式(5)的迭代形式，其中μ值選取采用分段函數(shù)形式為

3 實驗與仿真結(jié)果

為驗證改進正則化方法的有效性，通過建立三維模型，分別運用Tikhonov正則化方法(TIK)、改進的Tikhonov正則化方法(HTIK)和截斷奇異值分解的直接方法(TSVD)對3種電導率分布模型進行仿真[11]與實驗．

3.1ERT仿真實驗

采用Matlab軟件和有限元軟件COMSOL求解正問題．首先采用有限元軟件COMSOL建立三維模型，并用其剖分網(wǎng)格．采用3種不同電導率分布的三維模型，有限元單元采用四面體單元如圖1所示．模型(Ⅰ)共有13,885個節(jié)點，7,068個四面體單元；模型(Ⅱ)共有14,699個節(jié)點，7,711個四面體單元；模型(Ⅲ)共有15,703個節(jié)點，8,434個四面體單元．ERT的測量方式選擇相鄰驅(qū)動模式，電導率對比度1∶2．

反問題求解采用方形網(wǎng)格，網(wǎng)格數(shù)為812個，分別運用TIK算法、TSVD算法和HTIK算法進行成像．根據(jù)經(jīng)驗選擇能夠得出理想圖像的μ值．在仿真中對TIK算法選擇不同μ值，模型(Ⅰ)的μ值選為10-6，模型(Ⅱ)和模型(Ⅲ)的μ值選為0.02．迭代次數(shù)都選為20次，進行成像．表1給出了3種算法的成像結(jié)果，成像圖選取3維模型的中間截面如圖2所示．可以看出HTIK法成像效果明顯好于TIK法和TSVD法．

圖1 ERT三維模型Fig.1 3D model of ERT

圖2 三維模型截面示意Fig.2 Section of 3D model

表1 仿真模型和成像效果Tab.1 Simulation models and reconstructed images

由于仿真中實際電導率分布已知，所以采用圖像的相關(guān)系數(shù)對重建圖像質(zhì)量進行評判．對3種不同電導率分布模型分別采用TIK和HTIK進行成像比較．圖3為取50次迭代所成圖像的相關(guān)系數(shù)．從圖3中可以看出HTIK的圖像相關(guān)系數(shù)更高．TIK方法要耗費大量的時間進行μ值的選擇，而且圖像最大相關(guān)系數(shù)低于HTIK．

圖3 模型的圖像相關(guān)系數(shù)曲線Fig.3 Curves of correlation coefficient of models

對反問題采用不同數(shù)量的網(wǎng)格剖分[12]，分別對兩種成像算法的計算時間進行比較(計算機配置：Intel 酷睿2雙核P,7350，主頻2.0,GHz，內(nèi)存2.0,GB，操作系統(tǒng)為Windows7)，如表2所示．可以看出隨著逆問題網(wǎng)格剖分數(shù)量的增加，HTIK所花時間比TIK略少．由于傳統(tǒng)TIK算法對μ值的選擇要進行多次實驗觀測確定，所以實際上HTIK方法節(jié)省了大量時間．

表2 不同算法圖像成像實時性比較(50次迭代)Tab.2 Comparison of different algorithms in terms of computation time(iteration 50)

3.2ERT實驗結(jié)果

圖4 ERT成像實驗Fig.4 ERT imaging experiment

電阻成像實驗系統(tǒng)如圖4所示，激勵電流小于5,mA，激勵頻率范圍10,kHz～1,MHz(實驗激勵頻率為100,kHz)，激勵模式為相鄰激勵，采用相鄰測量的測量模式，解調(diào)方式為數(shù)字相敏解調(diào)，測量信噪比60,dB，數(shù)據(jù)采集速度為120幀/s，每幅含208個數(shù)據(jù)，USB2.0數(shù)據(jù)接口，上位機軟件采用VC6.0的界面，運用OpenGL進行圖像顯示．水槽直徑為200,mm；16個電極為直徑10,mm的圓形不銹鋼電極；溶液為NaCl溶液；成像物體為直徑25,mm有機玻璃棒；1,m長單層屏蔽電纜．

將NaCl加入水中測得溶液電導率為0.2,s/m，運用實驗裝置進行實驗．實驗模型參照仿真模型，有機玻璃棒放在鹽水里的分布見表3．分別采用HTIK(25次迭代)、TIK(100次迭代)、TSVD 3種方法進行成像．其中TIK中的正則化參數(shù)μ根據(jù)經(jīng)驗選擇，模型(Ⅰ)選為10-6、模型(Ⅱ)(Ⅲ)選為0.02．成像結(jié)果如表3所示．

表3 模型和實驗成像效果Tab.3 Models and reconstructed images for experimental results

在實際的測量中，由于測量誤差以及儀器本身的誤差導致實驗數(shù)據(jù)與仿真有差異[13]．由表3可以看出HTIK成像場域更加均勻，離散相更接近真實分布．HTIK成像效果優(yōu)于TIK和TSVD．

4 結(jié) 語

提出了一種改進的Tikhonov正則化方法用于ERT成像．基于同倫映射方法對正則化參數(shù)進行連續(xù)調(diào)節(jié)，然后用所得圖像灰度值作為Tikhonov迭代法的初始值進行成像，提高成像質(zhì)量．同傳統(tǒng)的Tikhonov正則化方法比較，HTIK具有成像質(zhì)量高，花費時間少的優(yōu)點，通過仿真和實驗證明了圖像重建的可行性．

［1］ Dickin F，Wang Mi. Electrical resistance tomography for process applications[J]. Meas Sci Technol，1996，7(3)：247-260.

［2］ Zlochiver S，F(xiàn)reimark D，Arad M，et al. ParametricEIT for monitoring cardiac stroke volume[J]. Physiol Meas，2006，27(5)：139-146.

［3］ Hu Li，Wang Huaxiang，Zhao Bo，et al. A hybrid reconstruction algorithm for electrical impedance tomography[J]. Meas Sci Technol，2007，18(3)：813-818.

［4］ 王化祥，范文茹，胡 理. 基于GMRES和Tikhonov正則化的生物電阻抗圖像重建算法[J]. 生物醫(yī)學工程學雜志，2009，26(4)：701-705.

Wang Huaxiang，F(xiàn)an Wenru，Hu Li. A hybrid reconstruction method in electrical impedance tomography based on GMRES and Tikhonov regularization[J]. Journal of Biomedical Engineering，2009，26(4)：701-705(in Chinese).

［5］ 彭黎輝，陸 耿. 電容成像圖像重建算法原理及評價[J]. 清華大學學報：自然科學版，2004，44(4)：478-484.

Peng Lihui，Lu Geng. Image reconstruction algorithms for electrical capacitance tomography：State of the art[J]. Journal of Tsinghua University：Science and Technology，2004，44(4)：478-484(in Chinese).

［6］ 王化祥，唐 磊，閆 勇. 電容層析成像圖像重建的總變差正則化算法[J]. 儀器儀表學報，2007，28(11)：2014-2018.

Wang Huaxiang，Tang Lei，Yan Yong. Total variation regularization algorithm for electrical capacitance tomography[J]. Chinese Journal of Scientific Instrument，2007，28(11)：2014-2018(in Chinese).

［7］ 王則柯. 同倫方法引論[M]. 重慶：重慶出版社，1990.

Wang Zeke. An Introduction to Homotopy Methods[M]. Chongqing：Chongqing Publishing House，1990(in Chinese).

［8］ 崔 凱，楊國偉，李興斯. 用同倫方法反演非飽和土中溶質(zhì)遷移參數(shù)[J]. 力學學報，2005，37(3)：307-312.

Cui Kai，Yang Guowei，Li Xingsi. A homotopy method for parameter inversion of solute transport through unsaturated soils[J]. Acta Mechanica Sinica，2005，37(3)：307-312(in Chinese).

［9］ Li Yi，Yang Wuqiang. Virtual electrical capacitance tomography sensor[J]. Journal of Physics：Conference Series，2005，15(1)：183-188.

［10］ Sun Meng，Liu Shi，Lei Jing，et al. Mass flow measurement of pneumatically conveyed solids using electrical capacitance tomography[J]. Measurement Science and Technology，2008，19(4)：1-6.

［11］ Nick P，William R B L. A Matlab toolkit for threedimensional electrical impedance tomography：A contribution to the electrical impedance and diffuse optical reconstruction software project[J]. Meas Sci Technol，2002，13(12)：1871-1883.

［12］ 王化祥，朱學明，張立峰. 用于電容層析成像技術(shù)的共軛梯度算法[J]. 天津大學學報，2005，38(1)：1-4.

Wang Huaxiang，Zhu Xueming，Zhang Lifeng. Conjugate gradient algorithm for electrical capacitance tomography[J]. Journal of Tianjin University，2005，38(1)：1-4(in Chinese).

［13］ Fagerberg A，Stenqvist O，Aneman A. Electrical impedance tomography applied to assess matching of pulmonary ventilation and perfusion in a porcine experimental model[J]. Crit Care，2009，13(2)：852-857.

Improved Regularization Reconstruction Algorithm Based on 3D Model for ERT System

LI Shou-xiao1，WANG Hua-xiang1，F(xiàn)AN Wen-ru2，ZHANG Ling-ling3
(1. School of Electrical Engineering and Automation，Tianjin University，Tianjin 300072，China；2. College of Aeronautical Automation，Civil Aviation University of China，Tianjin 300300，China；3. School of Sciences，Tianjin University，Tianjin 300072，China）

Electrical resistance tomography (ERT) is a technique for reconstructing the conductivity distribution inside an inhomogeneous distribution by injecting currents at the boundary ofa subject and measuring the resulting changes in voltage. To improve the poor imaging quality for ERT, an improved regularization reconstruction algorithm(HTIK)based on 3D modelling was presented. In general, the regularization coefficient of Tikhonov iterative algorithm was difficult to choose, homotopic mapping is adopted. A quasi-Sigmoid method is used to adjust the regularization parameter. Thereafter sensitive field image can be reconstructed by deploying the resultant image vector as the initial value of Tikhonov iteration algorithm. Both simulation and experimental results demonstrate that the proposed method can improve imaging quality obviously.

electrical resistance tomography；image reconstruction；regularization；homotopic mapping

TP212

A

0493-2137(2012)03-0215-06

2011-01-22；

2011-04-24.

國家自然科學基金重大國際合作資助項目(60820106002)；國家自然科學基金資助項目(60532020，50937005)；國家高技術(shù)研究發(fā)展計劃(863計劃)資助項目(2006BAI O3A00)；天津市自然科學基金資助項目(11JCYBJC06900).

李守曉（1985— ），男，博士研究生，lsxfly@tju.edu.cn.

王化祥，hxwang@tju.edu.cn.