乳腺腫瘤同倫LM算法EIT圖像重建

秦偉剛，王 超

（1. 天津大學(xué)電氣與自動化工程學(xué)院，天津 300072；2. 天津工業(yè)大學(xué)電氣工程與自動化學(xué)院，天津 300387）

乳腺腫瘤同倫LM算法EIT圖像重建

秦偉剛1,2，王 超1

（1. 天津大學(xué)電氣與自動化工程學(xué)院，天津 300072；2. 天津工業(yè)大學(xué)電氣工程與自動化學(xué)院，天津 300387）

將同倫延拓方法與Levenberg-Marquardt（LM）迭代算法相結(jié)合，對不同的乳腺腫瘤模型進(jìn)行圖像重建，并與LM算法進(jìn)行比較.在不同迭代初始值下，同倫LM算法與LM算法的平均電導(dǎo)率相對誤差相比，最優(yōu)值能降低接近7%,，運(yùn)行時(shí)間最大能減少46%,的計(jì)算開銷，電壓相對誤差要優(yōu)于LM算法7%,以上.實(shí)驗(yàn)結(jié)果表明，同倫LM算法比LM算法能重建更好的圖像，提高了重建圖像的精度和收斂速度.

電阻抗層析成像；同倫延拓法；Levenberg-Marquardt算法

與其他腫瘤相比，女性乳腺腫瘤的發(fā)病率最高［1］.乳腺腫瘤檢測方法主要分為兩大類：侵入式和非侵入式.在非侵入式方法中，電阻抗層析成像（electrical impedance tomography，EIT）是一種新的無損傷功能成像技術(shù)［2］，在生物醫(yī)學(xué)工程領(lǐng)域被應(yīng)用于如乳腺腫瘤檢測［3］、肺通氣功能檢測［4-5］和腦功能檢測［6］等的臨床研究中.由此可見，EIT可以用來研究病變組織的變化情況，對其成像研究仍然是有價(jià)值的.

通常，EIT逆問題具有非線性和欠定性，因此迭代算法被廣泛采用，如正則化牛頓法、Levenberg-Marquardt（LM）法和Landweber法［7-9］.這些圖像重構(gòu)算法在一定約束條件下都可以使求解問題收斂，但這些算法都是局部收斂的 .

LM算法是結(jié)合二次收斂牛頓法與梯度下降法的一種非線性優(yōu)化迭代方法，在EIT逆問題求解中獲得了很多應(yīng)用［10］.對LM算法來說，如果初始值不在收斂域內(nèi)，全局最小值就很難獲得.許秋平［11］使用LM算法求解非線性方程，雖然可以解決系數(shù)矩陣奇異的問題，但當(dāng)選取的初始值與真實(shí)值相距較遠(yuǎn)時(shí)，出現(xiàn)不收斂的情況.Kleefeld等［12］將LM算法用于參數(shù)估計(jì)問題，由于實(shí)際的電導(dǎo)率未知，初始值的選取較為困難.

同倫延拓法（homotopy continuation method）是一種解決非線性逆問題的強(qiáng)有力工具，對初值要求十分寬松，且具有全局收斂特性［13-14］，可以被用于EIT逆問題的求解［15］，并且與其他方法結(jié)合能進(jìn)一步提高計(jì)算效率.吳煥麗等［16］將牛頓法和同倫方法進(jìn)行結(jié)合，F(xiàn)u等［17］提出了一種Tikhonov正則化-同倫方法，Wu［18］設(shè)計(jì)了一種同倫Newton-Raphson方法，提高了收斂速度和計(jì)算精度.

由上可知，用LM算法來解非線性方程的解，要求初始值應(yīng)接近真實(shí)解，才能更好地收斂到真實(shí)解，實(shí)際上真實(shí)解并不知道.本文基于同倫方法和LM算法，給出了一種同倫LM算法（稱為HLM算法），對不同乳腺腫瘤模型進(jìn)行圖像重建.首先利用同倫延拓方法為EIT逆問題求解的LM算法產(chǎn)生一個較好的迭代初始值，然后再通過LM迭代算法來收斂到全局最優(yōu)解，并討論了平均電導(dǎo)率相對誤差、運(yùn)行時(shí)間和電壓相對誤差隨初始值的變化情況.

1　EIT數(shù)學(xué)描述

假設(shè)所研究的敏感場為似穩(wěn)場，內(nèi)部無電流源，根據(jù)Maxwell方程，在敏感場內(nèi)滿足Poisson方程

式中：σ為電導(dǎo)率分布；φ為敏感域內(nèi)的電勢分布.

在施加電流的電極上滿足

式中：jm為由電極注入的電流密度；n為外方向單位矢量；N為電極數(shù).

在測量電極上滿足

當(dāng)媒質(zhì)各向同性、電導(dǎo)率σ均勻分布時(shí)，方程（1）變?yōu)長aplace方程，即

2　EIT逆問題求解算法

2.1LM算法

EIT逆問題可表示為非線性方程

或者最小二乘問題，即將計(jì)算電壓和測量電壓的誤差2范數(shù)作為代價(jià)函數(shù)，記R（σ）=f（σ）-U ，則

式中：f（σ）表示注入電流后，對應(yīng)電導(dǎo)率分布σ的邊界計(jì)算電壓；U表示邊界測量電壓；表示歐幾里德范數(shù).

令代價(jià)函數(shù)

要滿足式（7），可通過標(biāo)準(zhǔn)Gauss-Newton法進(jìn)行迭代，即

式中：k為迭代數(shù)，k=1，2，…；J為Jacobian矩陣；JTJ稱為Hessian矩陣，且為對稱陣；M為有限元剖分單元數(shù).

為克服Hessian矩陣的嚴(yán)重病態(tài)性，以及避免步長過大，重新將代價(jià)函數(shù)寫為

求式（10）對σk+1的極小點(diǎn)，可得

式（11）就是所要求解的LM算法的迭代公式，其中λ＞0稱為阻尼因子，I為單位矩陣.當(dāng)λ較大時(shí)LM算法接近于梯度下降法，當(dāng)λ較小時(shí)接近于Gauss-Newton法.

2.2同倫延拓法

由于LM算法具有局部收斂性，同倫延拓法具有全局收斂性，且對初始值不敏感.構(gòu)造定點(diǎn)同倫方程［19］為

式中t為同倫參數(shù)，t∈［0,1］.式（12）滿足2個邊界條件

式中g(shù)（σ）為一簡單函數(shù)，與G（σ）為同倫映射.

2.3HLM算法

定義泛函

為了尋求式（14）的解，讓G（σ）表示為T（σ）的一階導(dǎo)數(shù)，即

同倫是拓?fù)鋵W(xué)中的基本概念，定點(diǎn)同倫的表達(dá)式可寫為

式中σ0表示初始電導(dǎo)率，它也是式（17）的解.

同倫方程（16）的解為σ隨t由0到1變化的一條連續(xù)曲線.根據(jù)薩德（Sard）定理，同倫方程對初始值的選擇幾乎都有解，不可解的測度為零［20］.

為簡化計(jì)算，方程（16）中t∈［0,1］取等間隔分布tk=k/K，同倫LM優(yōu)化算法的迭代公式為

當(dāng)?shù)螖?shù)小于K時(shí)，按照式（18）進(jìn)行迭代；大于K時(shí)，按照式（19）進(jìn)行迭代.目的是經(jīng)過式（18）的計(jì)算為式（19）尋求一個較好的初值，最終經(jīng)式（19）收斂到真實(shí)解.

由于逆問題中Jacobian矩陣的條件數(shù)很大，為提高成像質(zhì)量和收斂速度，對參數(shù)λ采用可調(diào)模式，初始值定義為

式中τ根據(jù)拇指規(guī)則［20］確定，取τ=0.1.

為解決同倫算法計(jì)算效率不高的問題，此處對式（18）中λ的取值為

式（19）中λ的取值為

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

本文采用16電極結(jié)構(gòu)，等間隔排列在敏感場的外表面，采取“相鄰電流激勵-相鄰電壓測量”模式.根據(jù)有限元法對其進(jìn)行三角元網(wǎng)格剖分，進(jìn)行圖像重建，剖分結(jié)果包含1,024個三角元和609個節(jié)點(diǎn).

實(shí)驗(yàn)運(yùn)行在Inter Core2 CPU2.4,GHz、RAM 2,GB的計(jì)算機(jī)上.以計(jì)算電壓和測量電壓差值的2范數(shù)作為停止迭代條件，允許誤差大小設(shè)為0.000,01，最大迭代次數(shù)設(shè)為30.根據(jù)惡性腫瘤和正常組織的差異性，病變組織的電導(dǎo)率要高于正常組織，文中設(shè)背景電導(dǎo)率為1,S/m，目標(biāo)電導(dǎo)率為5,S/m，給出了4種不同乳腺腫瘤模型分布，如圖1～圖4所示.上述算法通過MATLAB實(shí)現(xiàn)，圖1～圖4中第2列為LM算法重建的灰度圖像，第3列為HLM算法重建的灰度圖像.

圖1　模型分布1及其圖像重建Fig.1 Model distribution 1 and its image reconstruction

圖2　模型分布2及其圖像重建Fig.2 Model distribution 2 and its image reconstruction

圖3　模型分布3及其圖像重建Fig.3 Model distribution 3 and its image reconstruction

圖4　模型分布4及其圖像重建Fig.4 Model distribution 4 and its image reconstruction

4　分析和討論

為定量評價(jià)重建圖像的質(zhì)量，定義平均電導(dǎo)率相對誤差

式中：σc為計(jì)算電導(dǎo)率分布；σr為實(shí)際電導(dǎo)率分布.定義參數(shù)evr表示電壓相對誤差，即

由于同倫方法計(jì)算效率不高，針對不同的乳腺腫瘤模型，圖5給出了同倫法計(jì)算的不同初始電導(dǎo)率時(shí)的電壓相對誤差.

接下來，通過LM算法和HLM算法對平均電導(dǎo)率相對誤差進(jìn)行了對比，圖6描述了兩種算法在不同初始電導(dǎo)率時(shí)的計(jì)算結(jié)果.結(jié)果顯示，HLM算法計(jì)算的電導(dǎo)率相對誤差要小于LM算法的電導(dǎo)率相對誤差，并且兩種算法得到結(jié)果的變化趨勢是一致的，也說明HLM算法提高了重建圖像的精度.

雖然同倫方法可以擴(kuò)大算法的收斂域，但是并沒有消除對初值的依賴性，從圖5的結(jié)果來看，不同的初值下，前期迭代的電壓相對誤差變化較大，隨著迭代次數(shù)增加，電壓相對誤差變化一直處于下降的趨勢.對同一種模型分布來說，不同的初值情況，隨著迭代次數(shù)的增加，電壓相對誤差數(shù)值的下降也是比較明顯的.從圖6的結(jié)果來看，不同的模型分布取不同的初值迭代時(shí)，HLM算法與LM算法的平均電導(dǎo)率相對誤差相比，最大能降低接近7%,.實(shí)驗(yàn)發(fā)現(xiàn)，對于HLM算法來說，不同的初值，從同倫法切換到LM算法的迭代時(shí)刻不同，目的是為LM算法提供一個較好的初始值，這也說明引入同倫法后，迭代算法仍然是穩(wěn)定的.

為進(jìn)一步說明HLM算法的有效性，分別取不同電導(dǎo)率迭代初始值，利用LM算法和HLM算法計(jì)算不同模型分布的運(yùn)行時(shí)間，如圖7所示.可以發(fā)現(xiàn)，使用HLM算法計(jì)算的運(yùn)行時(shí)間比LM算法計(jì)算的運(yùn)行時(shí)間有所減少，最大能減少46%,的計(jì)算開銷，這也表明HLM算法進(jìn)一步提高了重建圖像的收斂速度.

最后，對兩種算法的電壓相對誤差進(jìn)行了比較，如圖8所示.同樣，HLM算法計(jì)算的電壓相對誤差明顯要優(yōu)于LM算法7%,以上，也體現(xiàn)了HLM算法優(yōu)于LM算法.

圖5　同倫法計(jì)算的不同初始電導(dǎo)率時(shí)的電壓相對誤差Fig.5 Relative errors of the voltage calculated by homotopy method for different initial conductivities

圖6　LM算法和HLM算法計(jì)算的平均電導(dǎo)率相對誤差Fig.6 Relative errors of the average conductivity calcu-lated by LM and HLM algorithms

從上述結(jié)果可以看出，與LM算法相比，HLM算法能提高圖像重建的精度和收斂速度，能獲得更好的圖像重建結(jié)果.

另一方面，由于EIT逆問題的軟場特性，即靈敏場分布隨著電導(dǎo)率的變化而改變，這種非線性增加了圖像重建的難度.文中HLM算法克服了LM算法的缺點(diǎn)，使得重建圖像更加穩(wěn)定.當(dāng)然，為了獲得更精確的解，像電極數(shù)目、電極形狀以及網(wǎng)格剖分的數(shù)量都需要進(jìn)一步的綜合考慮，以提高敏感場的靈敏度，為EIT逆問題求解提供可靠、有效的測量信息.

5　結(jié)　語

文中給出了一種優(yōu)化的圖像重建迭代算法——HLM算法，對不同的乳腺腫瘤模型進(jìn)行研究，利用具有大范圍收斂域的同倫延拓法進(jìn)行電導(dǎo)率修正，然后作為初始值代入LM算法進(jìn)行迭代，直至最后收斂，說明了HLM算法的可行性和有效性.從平均電導(dǎo)率相對誤差、運(yùn)行時(shí)間和電壓相對誤差來看，HLM算法提高了重建圖像的精度和收斂速度.另外，提高敏感場的靈敏度和基于阻抗成像算法仍然是研究的重點(diǎn)，并且值得進(jìn)一步探索.

［1］ Siegel R，Ma J M，Zou Z H，et al. Cancer statistics CA：A cancer［J］. Journal for Clinicians，2014，64：9-29.

［2］ Mamatjan Y，Grychtol B，Gaggero P，et al. Evaluation and real-time monitoring of data quality in electrical impedance tomography［J］. IEEE Transactions on Medical Imaging，2013，32（11）：1997-2005.

［3］ Kantartzis P，Abdi M，Liatsis P. Stimulation and measurement patterns versus prior information for fast 3D EIT：A breast screening case study［J］. Signal Processing，2013，93：2838-2850.

［4］ Camporota L，Smith J，Barrett N，et al. Assessment of regional lung mechanics with electrical impedance tomography can determine the requirement for ECMO in patients with severe ARDS［J］. Intensive Care Med，2012，38：2086-2087.

［5］ 王 超，白瑞峰，魏藝明，等. 基于小波變換肺部氣血阻抗分離及能量分析［J］. 天津大學(xué)學(xué)報(bào)：自然科學(xué)與工程技術(shù)版，2014，47（3）：195-199. Wang Chao，Bai Ruifeng，Wei Yiming，et al. Separation of pulmonary gas-blood bio-impedance and energy analysis based on wavelet transformation［J］. Journal of Tianjin University：Science and Technology，2014，47（3）：195-199（in Chinese）.

［6］ Bera T K，Biswas S K，Rajan K，et al. Projection error propagation-based regularization（PEPR）method for resistivity reconstruction in electrical impedance tomography（EIT）［J］. Measurement，2014，49：329-350.

［7］ Wang C，He X R，Bai R F. Trackability evaluation of reconstruction algorithms to the change of measured objects in electrical tomography［J］. Physiological Measurement，2014，35：583-596.

［8］ Kanzow C，Yamashita N，F(xiàn)ukushima M. Levenberg-Marquardt methods with strong local convergence properties for solving nonlinear equations with convex constraints［J］. Journal of Computational and Applied Mathematics，2005，173：321-343.

［9］ Wang H X，Wang C，Yin W L. A pre-iteration method for the inverse problem in electrical impedance tomography［J］. IEEE Trans on Instrum Meas，2004，53：1093-1096.

［10］ Barra L P S，Telles J C F. A geometric inverse problem identification procedure for detection of cavities［J］. Engineering Analysis with Boundary Elements，2013，37：1401-1407.

［11］ 許秋平. 電力調(diào)度系統(tǒng)中拓?fù)浞治雠c潮流計(jì)算方法的研究［D］. 濟(jì)南：山東大學(xué)數(shù)學(xué)學(xué)院，2011. Xu Qiuping. A Study of Topological Analysis and Flow Calculation in Power System Dispatching［D］. Jinan：School of Mathematics，Shandong University，2011（in Chinese）.

［12］ Kleefeld A，Rei?el M. The Levenberg-Marquardt method applied to a parameter estimation problem arising from electrical resistivity tomography［J］. Applied Mathematics and Computation，2011，217：4490-4501.

［13］ 黃象鼎，曾鐘鋼，馬亞南. 非線性數(shù)值分析的理論與方法［M］. 武漢：武漢大學(xué)出版社，2004. Huang Xiangding，Zeng Zhonggang，Ma Yanan. The Theory and Method for Nonlinear Numerical Analysis［M］. Wuhan：Wuhan University Press，2004（in Chinese）.

［14］ Huang Q Q，Zhu Z B，Wang X L. A predictor-corrector algorithm combined conjugate gradient with homotopy interior point for general nonlinear programming ［J］. Applied Mathematics and Computation，2013，219：4379-4386.

［15］ Liu S，Lei J，Li Z H，et al. Generalized flow pattern image reconstruction algorithm for electrical capacitance tomography［J］. Nuclear Engineering and Design，2011，241：1970-1980.

［16］ 吳煥麗，徐桂芝，張 帥，等. 基于圓柱模型的三維電阻抗成像問題研究［J］. 山東大學(xué)學(xué)報(bào)：理學(xué)版，2009，44（5）：45-48. Wu Huanli，Xu Guizhi，Zhang Shuai，et al. Research of the three-dimensional electrical impedance tomography based on the cylinder model ［J］. Journal of Shandong University：Natural Science，2009，44（5）：45-48（in Chinese）.

［17］ Fu H S，Han B. Tikhonov regulation-homotopy method for electrical impedance tomography［J］. Journal of Natural Science of Heilongjiang University，2011，28（3）：319-323.

［18］ Wu T M. Solving the nonlinear equations by the Newtonhomotopy continuation method with adjustable auxiliary homotopy function［J］. Applied Mathematics and Computation，2006，173：383-388.

［19］ 徐桂芝，李 穎，楊 碩，等. 生物醫(yī)學(xué)電阻抗成像技術(shù)［M］. 北京：機(jī)械工業(yè)出版社，2010. Xu Guizhi，Li Ying，Yang Shuo，et al. Electrical Impedance Tomography in Biomedical Engineering［M］. Beijing：China Machine Press，2010（in Chinese）.

［20］ Madson K，Nielsen H B，Tingleff O. Methods for Non-Linear Least Squares Problems（Informatics and Mathematical Modelling）［M］. 2nd ed.Copenhagen，Denmark：Technical University of Denmark，2004.

（責(zé)任編輯：孫立華）

Image Reconstruction of Breast Tumors Using Homotopy LM Algorithm in Electrical Impedance Tomography

Qin Weigang1,2，Wang Chao1
（1.School of Electrical Engineering and Automation，Tianjin University，Tianjin 300072，China；2.School of Electrical Engineering and Automation，Tianjin Polytechnic University，Tianjin 300387，China）

A homotopy Levenberg-Marquardt（HLM）algorithm，which is the combination of homotopy continuation method and Levenberg-Marquardt（LM）algorithm，was presented.The reconstructed images of different breast tumor models were acquired by HLM algorithm and compared with those by LM algorithm.The relative errors of the average conductivity were obtained by means of HLM algorithm and LM algorithm.The optimal value of the relative errors of the average conductivity by HLM algorithm was reduced by nearly 7 percent for different initial values.The best running time and the relative errors of the voltage are decreased by 46 percent and 7 percent respectively.The results show that the better images are reconstructed by HLM algorithm，and the accuracy and convergence rate are improved.

electrical impedance tomography（EIT）；homotopy continuation method；Levenberg-Marquardt（LM）algorithm

R318

A

0493-2137（2016）05-0506-07

10.11784/tdxbz201504037

2015-04-10；

2015-08-26.

國家自然科學(xué)基金重點(diǎn)資助項(xiàng)目（50937005）.

秦偉剛（1978—），男，博士研究生，qinweigang@tjpu.edu.cn.

王 超，wangchao@tju.edu.cn.

網(wǎng)絡(luò)出版時(shí)間：2015-11-11. 網(wǎng)絡(luò)出版地址：http://www.cnki.net/kcms/detail/12.1127.N.20151111.1802.010.html.