三種電阻抗斷層動態(tài)成像算法的仿真比較研究

【作 者】帥萬鈞，董秀珍，呂曉寧，晁勇

1 解放軍總醫(yī)院第一附屬醫(yī)院醫(yī)學(xué)工程科，北京，100048

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

3 海軍總醫(yī)院全軍航海航空醫(yī)學(xué)?？浦行模本?，100048

0 前言

電阻抗斷層成像技術(shù)（Electrical Impedance Tomography，EIT）是新興的功能成像技術(shù)，它根據(jù)生物體內(nèi)不同組織在不同生理或病理條件下具有不同的電阻抗特性這一原理，通過給目標(biāo)體內(nèi)注入安全微弱的交流電信號，測量體表電極間電位差來重建體內(nèi)的電阻抗分布或其變化的圖像。與現(xiàn)有醫(yī)學(xué)影像技術(shù)相比，該技術(shù)不使用核素或放射線，對人體無電離或輻射損傷，可多次測量、重復(fù)使用。而且，其設(shè)備結(jié)構(gòu)簡單、測量簡便、成本低廉，是一種具有誘人應(yīng)用前景的新型無損傷功能成像技術(shù)[1,2]。

EIT按成像目標(biāo)劃分，可分為靜態(tài)成像和動態(tài)成像。靜態(tài)成像是以電阻抗分布的絕對值為成像目標(biāo)，而動態(tài)成像是以電阻抗分布的相對值（差異）為成像目標(biāo)。一般而言，靜態(tài)成像能夠獲得電阻抗的絕對分布信息，但對硬件采集系統(tǒng)要求苛刻，重建算法復(fù)雜，計算量大，目前僅限于實驗室研究。動態(tài)成像以變化量為重建目標(biāo)，使用測量數(shù)據(jù)的差值進(jìn)行成像，減少了測量的系統(tǒng)誤差，使該技術(shù)有望由實驗室研究進(jìn)入臨床應(yīng)用。

EIT成像算法一直以來就是研究的熱點之一[3,4]。其動態(tài)成像算法主要有反投影算法(Back-Projection Algorithm)、牛頓一步重建法(Newton’s One-Step Error Reconstructor,NOSER)和靈敏度矩陣重建算法(Sensitivity Matrix Algorithm)。本文基于仿真研究，對此三種EIT動態(tài)成像算法的性能進(jìn)行了比較，以探討各自的優(yōu)劣和適應(yīng)性。

1 三種EIT動態(tài)成像算法

EIT所研究的是一個具有特殊邊界條件的電場，其場域的數(shù)學(xué)描述可用一組偏微分方程來表達(dá)，EIT的逆問題求解可分為迭代數(shù)值計算和非迭代近似計算兩大類[4]。非迭代線性近似方法常用于動態(tài)EIT成像，可表達(dá)為：

其中，△U為邊界電壓變化向量，△C為電導(dǎo)率變化向量，S為敏感矩陣，它表示△C與△U之間的線性敏感關(guān)系，由式（1）進(jìn)一步可得：

式(2)表明可通過邊界測量電壓的變化量計算區(qū)域內(nèi)擾動電導(dǎo)率的分布，此即為動態(tài)EIT重建的原理。其中敏感矩陣S是動態(tài)的，電極的位置和激勵測量模式以及場域的邊界外形都直接影響矩陣S求逆的穩(wěn)定性。

1.1 反投影算法

等位線反投影算法思想源于X射線計算機(jī)斷層成像（X-CT）的反投影技術(shù)[5]。其思路是：在電導(dǎo)率變化較小時，等電位線間的電導(dǎo)率的平均變化與相應(yīng)等位線間的電位差變化成線性關(guān)系，并認(rèn)為等位線分布在電導(dǎo)率變化前后不變。因而以相鄰兩條等位線為一個投影區(qū)域，將等位線間的電位差變化在此區(qū)域內(nèi)投影，再將不同反投影區(qū)域內(nèi)的反投影結(jié)果疊加起來，就可以得到電導(dǎo)率分布變化的結(jié)果。

該處理使得EIT重建不必對動態(tài)敏感矩陣S進(jìn)行求逆運算，可通過近似線性處理直接得到反投影矩陣B，從而有△C=B·△U，使得計算量小，重建迅速。反投影矩陣B的計算方法如下：在將場域離散后得到的有限元模型上，假定各個有限單元的電導(dǎo)率均勻，利用有限元法計算EIT正問題，得到場域內(nèi)所有節(jié)點的電位分布；然后根據(jù)每次激勵下各個電極的電位和各個單元節(jié)點的電位大小，將各個單元劃分到相應(yīng)測量電極對的投影區(qū)域中。設(shè)有限元模型的單元數(shù)為M，EIT的電極對為N，則B為(M×N)矩陣，有△Ci=Bij·△Uj,(i=1,2,Λ,M; j=1,2,Λ,N)，Bij表示第j個邊界電壓變化投影到第i個有限單元時的系數(shù)。當(dāng)該單元電位位于第j對測量電極電位之間時，令Bij=1，否則Bij=0。

1.2 牛頓一步重建法

NOSER算法是由美國Cheney教授于1990年首次提出的[6]。該算法實質(zhì)是Newton法一步。Cheney認(rèn)為，如果Newton法中設(shè)定的初始電導(dǎo)率足夠接近于真實電導(dǎo)率分布，則Newton法中的迭代次數(shù)會大大減小，甚至只需要在初始電導(dǎo)率的基礎(chǔ)上迭代一次就可以得到場域電導(dǎo)率分布的近似解。而且，Newton法之所以運算緩慢的原因在于每次迭代時都需要重新計算一次雅可比矩陣，如果只計算一步，運算量將大大降低。

NOSER算法確定初始電導(dǎo)率的策略是：在假設(shè)被測場域的電導(dǎo)率均勻分布的前提下，利用最小二乘法求得當(dāng)前邊界電壓下所對應(yīng)的場域電導(dǎo)率的分布。然后再計算雅可比矩陣，進(jìn)行圖像重建。

以此為基礎(chǔ)，重慶大學(xué)羅辭勇博士進(jìn)一步提出了FNOSER算法[7]，即快速(fast)NOSER算法。由于不需迭代，所以可事先設(shè)定初始電導(dǎo)率的分布，計算出雅可比矩陣。在實際檢測中，直接讀取該矩陣數(shù)據(jù)進(jìn)行圖像重建計算，從而節(jié)省大量運算時間。

1.3 靈敏度矩陣重建算法

靈敏度矩陣重建算法是由Murai和Kagawa于1985年提出來的[8]。其主要思想是基于Geselowitz的阻抗電場理論，利用有限元方法將場域等效為一電阻網(wǎng)絡(luò)模型，根據(jù)這一網(wǎng)絡(luò)模型建立總體傳導(dǎo)阻抗變化與單個單元中電阻抗變化之間的對應(yīng)關(guān)系，從而利用這一對應(yīng)關(guān)系來修正阻抗分布。

Geselowitz的阻抗電場理論為：

其中σ為電導(dǎo)率初始分布，其相應(yīng)的電位分布為φ(σ)，同時在電導(dǎo)率分布變化后的邊界電位為ψφ(σ+△σ)，且△σ為電導(dǎo)率分布變化量。EIT成像的主要問題是通過測量邊界電流驅(qū)動所導(dǎo)致的邊界電壓或邊界節(jié)點間的傳導(dǎo)阻抗Z來估計電導(dǎo)率分布σ+△σ，電位梯度 和 是電導(dǎo)率分布的函數(shù)。若對于已知的σ，則只需求出△σ即可得到新的電導(dǎo)率分布。進(jìn)一步近似處理，有：

上式中的Sij表示第i對電極對于第j個單元的電導(dǎo)率靈敏度系數(shù)。S被稱為靈敏度矩陣，對其求逆，即可得電導(dǎo)率變化值△σ。

2 仿真研究方法

所建立的場域有限元仿真模型如圖1所示，在圓模型邊界等間隔放置16個電極。為盡量模擬人體參數(shù),模型半徑設(shè)為14 cm，采用三角單元剖分，共有512個剖分單元。首先設(shè)置各個剖分單元電導(dǎo)率均勻分布，大小為0.2 S/m。從邊界注入1 mA激勵電流，通過數(shù)值計算得到擾動前的邊界電極上的電壓數(shù)據(jù)，以模擬實際檢測中邊界上所測量得到的電壓數(shù)據(jù)Uo，然后設(shè)置擾動目標(biāo)，其阻抗擾動量的大小按背景電導(dǎo)率的百分比計算，再通過數(shù)值計算獲得存在擾動后的邊界電極上的電壓數(shù)據(jù)，記為U1。在實際應(yīng)用中，U0和U1即為兩個不同時刻檢測的邊界電壓。對U0和U1按EIT算法進(jìn)行圖像重建，可獲得到擾動前后所發(fā)生的阻抗分布變化圖像。

圖1 有限元仿真模型及電極分布Fig.1 Finite element Model and electrodes distribution

另外,為進(jìn)一步討論各算法的抗噪性能，在實際測量中用數(shù)據(jù)噪聲進(jìn)行了模擬。邊界測量數(shù)據(jù)的加噪方法如下：

其中，Un為加噪前邊界測量數(shù)據(jù)，U'n加噪后邊界數(shù)據(jù)，Gn為高斯分布序列，其均值為0，方差為1，Kn為噪聲幅度系數(shù)，a為噪聲強(qiáng)度系數(shù)。利用計算機(jī)可以快速有效生成Gn，Kn計算方法如下：

其中，為邊界數(shù)據(jù)的獨立測量數(shù)。

3 結(jié)果及討論

有限元仿真模型的參數(shù)如上所述，目標(biāo)的擾動量設(shè)為10%，其位置和尺寸如圖2第一行所示。圖2和圖3中“BP”、“NS”和“SE”分別表示該行圖像為“反投影算法”、“NOSER”算法和“靈敏度矩陣算法”的成像。

3.1 圖像分辨率的比較

從圖2可見，三種算法對于徑向不同位置的目標(biāo)均能進(jìn)行有效重建(前三列)；對于位置相距較遠(yuǎn)的多目標(biāo)，其重建圖像均能反映出目標(biāo)的位置信息(后兩列)；但對于相距較近的目標(biāo)，重建圖像僅能顯示一個目標(biāo)，無法分辨出實際目標(biāo)個數(shù)(第四列)。相比較而言，NOSER算法的重建圖像分辨率高，目標(biāo)區(qū)域突出，靈敏度矩陣算法次之，但在二者圖像的目標(biāo)區(qū)域附近存在明顯亮斑偽影。反投影算法的目標(biāo)區(qū)域模糊，圖像分辨率較低，但無亮斑偽影?？傮w上講，在理想無噪聲條件下，NOSER算法和靈敏度矩陣算法的圖像分辨率較高；反投影算法的圖像平滑，偽影相對較少。

3.2 抗噪性能的比較

圖2 仿真模型上的反投影算法、NOSER算法和靈敏度矩陣算法的重建圖像Fig.2 Images reconstructed by back-projection algorithm,Newton’s one-step error reconstructor algorithm and sensitivity matrix algorithm based on finite element model

前面的仿真研究中，邊界電壓數(shù)據(jù)是通過計算機(jī)模擬計算得到的，它們不含有任何的誤差和噪聲，是一種理想的測量數(shù)據(jù)。在實際應(yīng)用中，測量數(shù)據(jù)不可能不耦合一定的噪聲。因此，在一定噪聲背景下研究算法的重建效果很有必要。以下按公式(6)為邊界測量數(shù)據(jù)添加噪聲,取噪聲強(qiáng)度系數(shù)分別為0、0.01%、0.05%、0.1%和0.5%，其中0表示不添加噪聲，如圖3各列所示。

在圖3中，NOSER算法(“NS”)的重建圖像在噪聲強(qiáng)度達(dá)到0.05%時，其圖像分辨率已顯著下降；噪聲強(qiáng)度達(dá)到0.5%時，從重建圖像上已無法分辨目標(biāo)，重建失敗。靈敏度矩陣算法(“SE”)受噪聲影響更為嚴(yán)重，在0.1%的噪聲強(qiáng)度下就無法對擾動目標(biāo)成像。反投影算法(“BP”)在噪聲強(qiáng)度小于0.1%時，擾動目標(biāo)均可被有效重建成像，即使噪聲強(qiáng)度達(dá)到0.5%，EIT成像才略受影響，但從重建圖像上仍可分辨出擾動目標(biāo)。因此，就抗噪性能而言，反投影算法最佳，NOSER算法次之，靈敏度矩陣算法最差。

圖3 不同噪聲條件下反投影、NOSER和靈敏度矩陣算法的重建成像Fig.3 Images reconstructed by back-projection algorithm,newton’s one-step error reconstructor algorithm and sensitivity matrix algorithm with different noise conditions

3.3 重建速度的比較

在現(xiàn)有條件(Intel Pentium 4 CPU 2.8 G，256 M內(nèi)存)和剖分規(guī)模(512個有限單元)下，進(jìn)行一次圖像重建，反投影算法耗時不足1秒，而NOSER算法和靈敏度矩陣算法至少需要十幾分鐘、甚至幾十分鐘，主要原因在于它們均須進(jìn)行一個大矩陣的正則化和求逆運算。

當(dāng)然，可以按照快速(fast)NOSER算法思想[7]，事先以文件方式保存雅可比矩陣和靈敏度矩陣，重建時直接讀取矩陣數(shù)據(jù)而不再計算，以提高運算效率。但NOSER算法和靈敏度矩陣算法對電導(dǎo)率分布初值要求較高，事先計算好的相應(yīng)矩陣在實際應(yīng)用中難以使用。而反投影算法對電導(dǎo)率分布初值要求較低，卻可以事先按均勻分布計算好反投影矩陣，從而進(jìn)一步減小應(yīng)用中圖像重建時間。

4 結(jié)論

以上仿真研究結(jié)果表明：反投影算法的重建圖像雖然圖像分辨率不高，但其圖像平滑，無明顯偽影，且該算法具有較強(qiáng)的抗噪性能和快速的運算能力，在臨床檢測的強(qiáng)噪聲背景下，應(yīng)優(yōu)選此算法。NOSER算法的圖像分辨率最好，且抗噪性能適中，在對運算速度要求不高時可作為首選。靈敏度矩陣算法是一種整體性能適中的算法，可作為EIT成像的備選算法。

[1]RH Bayford.Bioimpedance tomography (electrical impedance tomography)[J].Annu Rev Biomed Eng,2006,8: 63-91.

[2]T Muders,H Luepschen,C Putensen.Impedance tomography as a new monitoring technique[J].Curr Opin Crit Care,2010,16(3):269-275.

[3]M Bodenstein,M David,K Markstaller.Principles of electrical impedance tomography and its clinical application[J].Crti Care Med,2009,37(2): 713-724.

[4]劉國強(qiáng).醫(yī)學(xué)電磁成像[M].北京: 科學(xué)出版社,2006

[5]DC Barrber,AD Seagar.Fast reconstruction of resistance images[J].Clin Phys Physiol Meas,1987,8(Suppl A): 47-54.

[6]M Cheney,D Isaacson,JC Newell,et al.NOSER: An algorithm for solving the inverse conductivity problem[J].Int J Imag Syst Technol,1990,2: 66-75.

[7]羅辭勇.基于快速牛頓一步誤差重構(gòu)的電阻抗成像算法和實驗研究[D].重慶大學(xué),博士學(xué)位論文,2005

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