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

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

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

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

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

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

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

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

引言

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

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

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

1 材料與方法

1.1 EIT算法原理

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

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

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

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

1.2 常用的電阻抗重建算法

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

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

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

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

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

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

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

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

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

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

(7)為L(zhǎng)1正則問題，可通過軟閾值算子策略進(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)矩陣：

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

對(duì)上式進(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ù)矩陣，將上述兩個(gè) 式子帶入

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

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

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

圖1 病人顱腦CT圖像

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

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

1.4 仿真實(shí)驗(yàn)設(shè)計(jì)

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

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

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

1.4.2 圖像評(píng)價(jià)指標(biāo)

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

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

2 結(jié)果

首先使用簡(jiǎn)單圓域模型，在理想條件下進(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ǎn)單圓域的仿真實(shí)驗(yàn)結(jié)果

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

然后基于三維真實(shí)顱腦模型，進(jìn)行仿真。分別在腦實(shí)質(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)了擾動(dòng)與偽影，尤其是在中心處的目標(biāo)發(fā)生了較大的形變。對(duì)于3種算法的重建結(jié)果，傳統(tǒng)L2加權(quán)算法獲取的重建圖像質(zhì)量低于另外兩種優(yōu)化的算法；L1-SB算法對(duì)應(yīng)的圖像偽影最少，目標(biāo)形變也較??；GRETI重建目標(biāo)均勻性好。

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

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

3 討論

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

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

圖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é)報(bào),2012,33(3):625-630.

[7] 常甜甜,魏雯婷,叢偉杰,等.電阻抗成像的稀疏重建算法[J].西安郵電學(xué)院學(xué)報(bào),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] 徐燦華.床旁電阻抗動(dòng)態(tài)圖像監(jiān)測(cè)的重構(gòu)算法及實(shí)驗(yàn)研究[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] 梁樂,徐燦華,沈志鵬,等.用于生物電磁場(chǎng)分析的腦部三層有限元模型構(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）；軍隊(duì)重大課題（AWS14C006）。

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

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