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

基于三維模型的改進(jìn)正則化ERT成像算法

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

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

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

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

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

ERT成像算法包括非迭代算法和迭代算法[4]．非迭代算法包括線性反投影、Tikhonov正則化算法、基于截?cái)嗥娈愔捣纸獾闹苯铀惴?；迭代算法包括landweber迭代算法、牛頓-拉夫遜迭代算法以及共軛梯度迭代算法．線性反投影算法雖然簡(jiǎn)單快速，但該算法成像精度較低．landweber迭代算法與最速梯度下降法類似，但其收斂性比較差[5]．牛頓-拉夫遜迭代算法雖然收斂快但計(jì)算量較大且具有局部收斂的缺點(diǎn)．Tikhonov正則化算法的成像質(zhì)量主要依靠正則化參數(shù)的選擇．最優(yōu)的正則化參數(shù)將得到質(zhì)量高的圖像．但該種算法中正則參數(shù)一般不易確定，通常需要通過(guò)實(shí)驗(yàn)觀測(cè)確定，而且確定的參數(shù)并非最優(yōu)值.

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

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

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

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

2 算 法

2.1Tikhonov正則化方法

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

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

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

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

2.2同倫映射

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

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

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

式(7)滿足

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

應(yīng)用正則化方法進(jìn)行求解

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

2.3改進(jìn)正則化算法

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

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

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

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

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

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

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

為驗(yàn)證改進(jìn)正則化方法的有效性，通過(guò)建立三維模型，分別運(yùn)用Tikhonov正則化方法(TIK)、改進(jìn)的Tikhonov正則化方法(HTIK)和截?cái)嗥娈愔捣纸獾闹苯臃椒?TSVD)對(duì)3種電導(dǎo)率分布模型進(jìn)行仿真[11]與實(shí)驗(yàn)．

3.1ERT仿真實(shí)驗(yàn)

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

反問(wèn)題求解采用方形網(wǎng)格，網(wǎng)格數(shù)為812個(gè)，分別運(yùn)用TIK算法、TSVD算法和HTIK算法進(jìn)行成像．根據(jù)經(jīng)驗(yàn)選擇能夠得出理想圖像的μ值．在仿真中對(duì)TIK算法選擇不同μ值，模型(Ⅰ)的μ值選為10-6，模型(Ⅱ)和模型(Ⅲ)的μ值選為0.02．迭代次數(shù)都選為20次，進(jìn)行成像．表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

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

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

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

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

3.2ERT實(shí)驗(yàn)結(jié)果

圖4 ERT成像實(shí)驗(yàn)Fig.4 ERT imaging experiment

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

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

表3 模型和實(shí)驗(yàn)成像效果Tab.3 Models and reconstructed images for experimental results

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

4 結(jié) 語(yǔ)

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

［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ī)學(xué)工程學(xué)雜志，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］ 彭黎輝，陸 耿. 電容成像圖像重建算法原理及評(píng)價(jià)[J]. 清華大學(xué)學(xué)報(bào)：自然科學(xué)版，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]. 儀器儀表學(xué)報(bào)，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］ 崔 凱，楊國(guó)偉，李興斯. 用同倫方法反演非飽和土中溶質(zhì)遷移參數(shù)[J]. 力學(xué)學(xué)報(bào)，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］ 王化祥，朱學(xué)明，張立峰. 用于電容層析成像技術(shù)的共軛梯度算法[J]. 天津大學(xué)學(xué)報(bào)，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.

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

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

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