李昊庭，徐燦華，劉本源，楊琳，董秀珍，付峰

第四軍醫(yī)大學(xué) 生物醫(yī)學(xué)工程系，陜西西安 710032

稀疏加權(quán)算法與GREIT算法在顱腦電阻抗成像中的對比研究

李昊庭，徐燦華，劉本源，楊琳，董秀珍，付峰

第四軍醫(yī)大學(xué) 生物醫(yī)學(xué)工程系，陜西西安 710032

稀疏加權(quán)算法（L1算法）與GREIT算法是近年來兩種較為熱門的顱腦電阻抗成像（Electrical Impedance Tomography，EIT）圖像重建優(yōu)化算法，基于不同的數(shù)學(xué)模型，這兩種算法性能不同。為了改善顱腦EIT重建圖像質(zhì)量并為其算法優(yōu)選提供依據(jù)，本文基于三維顱腦仿真模型開展了仿真研究，對比了傳統(tǒng)二階范數(shù)加權(quán)算法（L2算法）、稀疏加權(quán)算法以及GREIT算法在顱腦電阻抗圖像重建中的性能。仿真結(jié)果表明，相比傳統(tǒng)二階范數(shù)加權(quán)算法，兩種優(yōu)化的算法對顱腦EIT均有改善，而稀疏加權(quán)算法在圖像噪聲、形變誤差、位置誤差3項指標(biāo)的評價中性能最優(yōu)。在閾值函數(shù)的作用下，稀疏加權(quán)算法經(jīng)多步迭代抑制了圖像噪聲，突出了重建目標(biāo)。稀疏加權(quán)算法可以大大改善顱腦EIT效果，適用于顱腦EIT且對未來顱腦電阻抗重建算法擴展研究有重要參考價值。

顱腦電阻抗斷層成像；動態(tài)成像；圖像重建算法；算法比較；三維顱腦模型

引言

電阻抗成像（Electrical Impedance Tomography，EIT）技術(shù)的原理為通過安裝在物體表面的電極對物體有規(guī)律地施加激勵，若物體內(nèi)部存在阻抗變化，則會引起表面測量電極電位的變化，基于測量電極電位變化，結(jié)合相應(yīng)重建算法可以獲得物體內(nèi)部阻抗變化圖像。EIT技術(shù)應(yīng)用領(lǐng)域廣泛，涉及肺功能成像，腦部功能成像，乳腺癌檢測，腹部臟器功能成像等多個領(lǐng)域[1]。其中，顱腦電阻抗技術(shù)可以檢測出腦功能性活動引起的阻抗變化，也可以檢測出顱腦疾病如腦出血等引起的阻抗變化，具有誘人的應(yīng)用前景[2]。

第四軍醫(yī)大學(xué)課題組長期致力于顱腦EIT研究，在國際上率先開展了顱腦EIT的臨床實驗研究，并取得了一系列的研究成果[3-4]。課題組研究發(fā)現(xiàn)相比實驗室環(huán)境，臨床實驗中顱腦EIT面臨更為嚴(yán)重的噪聲與干擾。隨機噪聲、電極位移、接觸阻抗的變化等多種因素均會影響到EIT的重建圖像質(zhì)量，故在臨床實驗背景下對重建算法的性能提出了更高的要求。在過去的幾年中，圍繞著提高重建圖像質(zhì)量和穩(wěn)定性，國內(nèi)外多個研究小組對重建算法開展了深入研究，在傳統(tǒng)二階范數(shù)加權(quán)算法（L2算法）的基礎(chǔ)上提出了一系列優(yōu)化算法。隨著稀疏重建理論的發(fā)展，稀疏加權(quán)算法（L1算法）被應(yīng)用到EIT，相比二階范數(shù)加權(quán)算法，稀疏加權(quán)算法的抗噪聲能力強，可以避免對重建圖像邊界造成過多模糊[5-7]。2009年Adler小組提出一種名為GREIT的優(yōu)化算法并將其應(yīng)用到肺部阻抗成像[8]。

基于簡單仿真和部分實驗研究，這些優(yōu)化算法相對于傳統(tǒng)算法，對EIT圖像質(zhì)量均有改善。然而，在顱腦EIT成像背景下針對這兩種優(yōu)化算法開展的算法比較工作尚未有文獻(xiàn)報道。對此，本文基于三維顱腦仿真模型，模擬了EIT監(jiān)測腦出血過程。我們在顱腦仿真模型不同位置設(shè)置了與血液電導(dǎo)率相同的球形目標(biāo)，在添加一定的高斯噪聲后，使用傳統(tǒng)二階范數(shù)加權(quán)算法、稀疏加權(quán)算法與GREIT算法進(jìn)行圖像重建，并基于圖像噪聲、位置誤差、形變誤差3項指標(biāo)對重建結(jié)果進(jìn)行了量化評價。根據(jù)重建圖像與評價指標(biāo)，我們總結(jié)了3種算法的不同特點并針對哪種算法更適合顱腦EIT進(jìn)行了分析，為以后顱腦電阻抗算法的優(yōu)選與優(yōu)化打下了基礎(chǔ)。

1 材料與方法

1.1 EIT算法原理

EIT算法可以分為正問題和逆問題求解。正問題指的是已知場域內(nèi)部阻抗分布和邊界條件，求解場域電勢分布，其常用的求解方法為有限元方法[9]。逆問題又稱作圖像重建，則是給定邊界測量數(shù)據(jù)和邊界條件，求解場域內(nèi)部阻抗分布。在EIT中，當(dāng)場域內(nèi)部阻抗變化較小時，邊界測量數(shù)據(jù)與場域阻抗變化分布如下線性方程：

其中ρ表示場域電導(dǎo)率變化分布，v為邊界測量數(shù)據(jù)，S是敏感系數(shù)矩陣。

由于EIT技術(shù)中測量數(shù)遠(yuǎn)遠(yuǎn)超過未知數(shù)，因此逆問題求解存在病態(tài)性，不能直接對（1）式求逆，采取的解決方法是引入正則項來構(gòu)造如下的目標(biāo)方程，然后求滿足下式的最優(yōu)解：

ρ^為阻抗變化分布的最優(yōu)估計，λ為正則化參數(shù)，Φ(ρ)是正則化項。

1.2 常用的電阻抗重建算法

隨著EIT算法的發(fā)展，形成了多種重建算法，不同算法具有不同的成像特點。本文選取了傳統(tǒng)二階范數(shù)加權(quán)算法、稀疏加權(quán)算法和GREIT算法3種常用的重建算法用于算法比較。

1.2.1 傳統(tǒng)二階范數(shù)加權(quán)算法

按照二階范數(shù)正則化加權(quán)方法[10]，電阻抗目標(biāo)函數(shù)有如下形式：

對于這種求二階范數(shù)加權(quán)最優(yōu)化問題，常用的阻尼最小二乘解[11]的形式為：

W為正則化矩陣，可以控制解的一些特性。

1.2.2 稀疏加權(quán)算法

按照稀疏重建理論，其目標(biāo)函數(shù)如下：

由于一階范數(shù)的存在，不能直接對上式進(jìn)行求導(dǎo)，針對處理不可導(dǎo)項產(chǎn)生了一系列的稀疏加權(quán)算法，本文選取的是成像較好、收斂速度較快的SplitBregman算法（SB算法）[12]。

根據(jù)SplitBregman算法[13]，電導(dǎo)率分布可以按照以下3個步驟迭代求解：

其中(6)為二階范數(shù)優(yōu)化問題通過求導(dǎo)可得：

(7)為L1正則問題，可通過軟閾值算子策略進(jìn)行求解

其中Shrink為軟閾值算子

1.2.3 GREIT算法

EIT線性算法一般滿足如下形式：

其中x為所求電導(dǎo)率分布，y為邊界電位變化，R為重構(gòu)矩陣。為了尋取最優(yōu)的重構(gòu)矩陣R，GREIT算法首先定義一系列的訓(xùn)練目標(biāo)t(i)，然后尋求使下式誤差最小的R作為優(yōu)化的重構(gòu)矩陣：

對于每一個訓(xùn)練目標(biāo)，均可以生成與之對應(yīng)的邊界電位y和理想的重構(gòu)值分布x。w是加權(quán)矩陣，代表訓(xùn)練目標(biāo)集合t中每一個目標(biāo)的權(quán)重。

對上式進(jìn)行求導(dǎo)：

由上式，重構(gòu)矩陣R為：

已知訓(xùn)練目標(biāo)集合t與理想重構(gòu)分布x、邊界電位y滿足如下關(guān)系：

D為訓(xùn)練目標(biāo)集合與理想重構(gòu)分布之間的映射矩陣，J為敏感系數(shù)矩陣，將上述兩個 式子帶入

由(19)可知，相比傳統(tǒng)算法，GREIT算法的核心在于映射矩陣D。D的作用為使訓(xùn)練目標(biāo)集合t可以映射到符合EIT特點且均勻性較好的理想重構(gòu)分布x。

1.3 三維真實顱腦仿真模型構(gòu)建

基于Mimics軟件對顱腦CT進(jìn)行分割和三維表面重建，然后使用SolidWorks對三維表面模型實體轉(zhuǎn)換，將實體化的模型導(dǎo)入Comsol4.4生成可用于顱腦EIT分析的有限元模型[14-15]，見圖1。構(gòu)建出的模型結(jié)構(gòu)完整，包括頭皮層、顱骨、腦脊液、腦實質(zhì)以及腦室，見圖2。

圖1 病人顱腦CT圖像

圖2 導(dǎo)入Comsol后生成的三維有限元模型

分別對不同組織賦予不同的電導(dǎo)率[16]，基于該模型，可以較為真實的模擬顱腦疾病及顱腦電阻抗技術(shù)應(yīng)用過程。

1.4 仿真實驗設(shè)計

1.4.1 真實顱腦邊界重構(gòu)模板

選取電極所在平面對三維模型進(jìn)行切割，導(dǎo)出切割面。手動選取中心點與邊界點，然后對顱腦切面圖進(jìn)行自適應(yīng)剖分，生成重構(gòu)模板，見圖3。

圖3 真實顱腦邊界重構(gòu)模板

1.4.2 圖像評價指標(biāo)

我們參考了Adler等在2009年提出的算法評價指標(biāo)體系并使用圖像噪聲（Image Noise，IN）、位置誤差（Location Error，LE）、形變誤差（Shape Error，SE）3項指標(biāo)對重建圖像進(jìn)行評價。首先，定義重建圖像中單元像素值大于最大像素值四分之一的集合為感興趣區(qū)域集合，記作ΩROI，感興趣區(qū)域以外的區(qū)域為Ωother。

dPre為參考圖像中預(yù)設(shè)目標(biāo)距離中心點的距離，dRe為重建圖像重建目標(biāo)距離中心點的距離，l為顱腦長軸長度，LE越小位置誤差越小。

2 結(jié)果

首先使用簡單圓域模型，在理想條件下進(jìn)行仿真。在Comsol中構(gòu)建圓域模型，圓域半徑為0.2，背景電導(dǎo)率為1，見圖4(a)。在圓域1/2處設(shè)置電導(dǎo)率為1.2，半徑為0.01的圓形目標(biāo)，生成邊界電位數(shù)據(jù)后分別使用L2加權(quán)算法、L1-SB算法、GREIT算法進(jìn)行圖像重建，見圖4。

圖4 基于簡單圓域的仿真實驗結(jié)果

由圖4簡單圓域仿真結(jié)果可知，理想情況下，3種重建算法均能獲得質(zhì)量較好的重建圖像，可以對目標(biāo)準(zhǔn)確定位。但由重建圖像對應(yīng)的擴散函數(shù)來看，3種算法重建特征存在較大差異。稀疏重建算法擴散函數(shù)在背景處幾乎為0，而在目標(biāo)邊界處變化極為劇烈，對應(yīng)的重建圖像背景偽影少，重建目標(biāo)邊界明確；GREIT算法在目標(biāo)內(nèi)部擴散函數(shù)變化平緩，對應(yīng)的重建目標(biāo)較為均勻。

然后基于三維真實顱腦模型，進(jìn)行仿真。分別在腦實質(zhì)不同位置設(shè)置體積為5 cm3，電導(dǎo)率為0.7 Ω/m的球形目標(biāo)模擬腦出血，并在生成的重建數(shù)據(jù)中添加0.1%的高斯噪聲，使用L2算法、稀疏加權(quán)算法、GREIT算法進(jìn)行圖像重建，見圖5。

由圖5重建結(jié)果可知，相比理想重構(gòu)， 基于顱腦三維仿真模型并添加高斯噪聲生成的數(shù)據(jù)經(jīng)重建后圖像出現(xiàn)了擾動與偽影，尤其是在中心處的目標(biāo)發(fā)生了較大的形變。對于3種算法的重建結(jié)果，傳統(tǒng)L2加權(quán)算法獲取的重建圖像質(zhì)量低于另外兩種優(yōu)化的算法；L1-SB算法對應(yīng)的圖像偽影最少，目標(biāo)形變也較??；GRETI重建目標(biāo)均勻性好。

圖5 基于三維真實顱腦模型的仿真結(jié)果

與重建圖像結(jié)果對應(yīng)，在圖像噪聲、位置誤差和形狀誤差三項指標(biāo)上，稀疏加權(quán)算法的性能要優(yōu)于其他兩種算法，見圖6。

3 討論

本文在顱腦EIT應(yīng)用背景下，基于三維顱腦仿真模型比較了傳統(tǒng)二階加權(quán)算法、稀疏加權(quán)算法與GREIT算法的性能。根據(jù)仿真實驗結(jié)果可知，相比于傳統(tǒng)的二階范數(shù)加權(quán)算法，兩種優(yōu)化的算法對顱腦電阻抗圖像重建均有一定的改善作用，但是基于不同的數(shù)學(xué)模型，兩者的重建特性不同。稀疏重建算法的突出特點為抗噪性能好，可以保留邊界，抑制由重建帶來的形狀誤差，但是相比其他兩種算法稀疏重建算法的需要調(diào)節(jié)的參數(shù)較多，在實際使用中增加了調(diào)參的負(fù)擔(dān)[17]；而GREIT算法在加權(quán)矩陣D的作用下其解分布比較均勻，噪聲相對較少。對于顱腦EIT成像，相比均勻性，抑制偽影與減少目標(biāo)形變對于臨床應(yīng)用更為重要，故我們認(rèn)為稀疏加權(quán)重建算法更適用于顱腦EIT。

另外，我們的對比研究中存在幾點需要進(jìn)一步完善的地方。首先，本文所有的參數(shù)選擇均為手動選擇，基于該方法選取的參數(shù)可能不是最優(yōu)。但只要保證3種算法進(jìn)行重建時采用的參數(shù)相同，算法對比工作就不會受到本質(zhì)影響。其次，本文在仿真模型上驗證了算法的性能，下一步我們將開展物理模型實驗、活體實驗來驗證算法性能。

圖6 圖像重建誤差堆狀圖

[1] 徐燦華,董秀珍.生物電阻抗斷層成像技術(shù)及其臨床研究進(jìn)展[J].高電壓技術(shù),2014,40(12):3738-3745.

[2] Brown BH,Barber DC,Seagar AD.Applied potential tomography: possible clinical applications[J].Clin Phys Physiol Meas,1985,6 (2):109-121.

[3] Fu F,Li B,Dai M,et al.Use of electrical impedance tomography to monitor regional cerebral edema during clinical dehydration treatment[J].PloS One,2014,9(12):0113202.

[4] Dai M,Li B,Hu S,et al.In vivo imaging of twist drill drainage for subdura l hematoma: a clinical feasibility study on electrical impedance tomography for measuring intracranial bleeding in humans[J].PLoS One,2013,8(1):055020.

[5] Borsic A,Graham BM,Adler A,et al.Total Variation Regularization in Electrical Impedance Tomography[J].Inverse Probl, 2007,99(A12):A12.

[6] 范文茹,王化祥,郝魁紅,等.基于兩步迭代TV正則化的電阻抗圖像重建算法[J].儀器儀表學(xué)報,2012,33(3):625-630.

[7] 常甜甜,魏雯婷,叢偉杰,等.電阻抗成像的稀疏重建算法[J].西安郵電學(xué)院學(xué)報,2013,18(2):92-96,110.

[8] Adler A,Arnold J H,Bayford R,et al.GREIT: a unif ed approach to 2D linear EIT reconstruction of lung images[J].Physiol Meas, 2009,30(6):S35-S55.

[9] Murai T,Kagawa Y.Electrical impedance computed tomography based on a finite element model[J].IEEE Trans Biomed Eng, 1985,(3):177-184.

[10] Lukaschewitsch M,Maass P,Pidcock M.Tikhonov regularization for electrical impedance tomography on unbounded domains[J]. Inverse Probl,2003,19(3):585.

[11] 徐燦華.床旁電阻抗動態(tài)圖像監(jiān)測的重構(gòu)算法及實驗研究[D].西安:第四軍醫(yī)大學(xué),2010.

[12] Zhou Z,Dos Santos GS,Dowrick T,et al.Comparison of total variation algorithms for electrical impedance tomography[J]. Physiol Meas,2015,36(6):1193.

[13] Wang J,Ma J,Han B,et al.Split Bregman iterative algorithm for sparse reconstruction of electrical impedance tomography[J]. Signal Process,2012,92(12):2952-2961.

[14] 梁樂,徐燦華,沈志鵬,等.用于生物電磁場分析的腦部三層有限元模型構(gòu)建[J].醫(yī)療衛(wèi)生裝備,2014,35(9):12-14.

[15] Markus J,James A,Emma M,et al.Correcting electrode modelling errors in EIT on realistic 3D head models[J].Physiol Meas,2015,36(12):2423.

[16] Tang C,You F,Cheng G,et al.Correlation between structure and resistivity variations of the live human skull[J].IEEE Trans Biomed Eng,2008,55(9):2286-2292.

[17] Tehrani JN,McEwan A,Jin C,et al.L1 regularization method in electrical im pedance tomography by using the L1-curve (Pareto frontier curve)[J].Appl Math Model,2012,36(3):1095-1105.

Comparative Study on the Application of Sparse Constrained Algorithm and GREIT in Brain Electrical Impedance Tomography

LI Hao-ting, XU Can-hua, LIU Ben-yuan, YANG Lin, DONG Xiu-zhen, FU Feng
Faculty of Biomedical Engineering, the Forth Military Medical University, Xi'an Shaanxi 710032, China

Sparse constrained and GREIT algorithms were two hot optimized reconstruction algorithms for brain Electrical Impedance Tomography (EIT). With different mathematical models, these two algorithms had different performances. To give recommendations to brain EIT algorithm selection, based on the three dimensional head model, a comparison of the performances of conventional quadratic norm regularization algorithm (L2 algorithm), sparse constrained algorithm and GREIT algorithm was made in this paper. On the basis of evaluation results, it could be concluded that sparse constrained algorithm and GREIT demonstrated better performances than L2 algorithm in image noise, location error and shape error. With threshold function, the sparse constrained algorithm suppressed the image noise by multi-step iteration, and highlighted the reconstruction object. Sparse constrained algorithm could greatly improve the performances of brain EIT and was suitable for brain EIT, which was of great reference for the researches on the brain EIT reconstruction algorithm EIT in the future.

brain electrical impedance tomography; dynamic imaging; reconstruction algorithm; algorithm comparison; three dimensional head model

R318；TM934.7

A

10.3969/j.issn.1674-1633.2016.11.005

1674-1633(2016)11-0023-05

2016-09-28

2016-10-10

國家自然科學(xué)基金（51477176）；軍隊重大課題（AWS14C006）。

付峰，教授，主要研究方向生物醫(yī)學(xué)電阻抗成像。

通訊作者郵箱：fengfu@fmmu.edu.cn