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

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

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

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

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

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

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

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

1.1 反投影算法

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

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

1.2 牛頓一步重建法

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

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

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

1.3 靈敏度矩陣重建算法

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

Geselowitz的阻抗電場(chǎng)理論為：

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

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

2 仿真研究方法

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

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

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

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

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

3 結(jié)果及討論

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

3.1 圖像分辨率的比較

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

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

在圖3中，NOSER算法(“NS”)的重建圖像在噪聲強(qiáng)度達(dá)到0.05%時(shí)，其圖像分辨率已顯著下降；噪聲強(qiáng)度達(dá)到0.5%時(shí)，從重建圖像上已無法分辨目標(biāo)，重建失敗。靈敏度矩陣算法(“SE”)受噪聲影響更為嚴(yán)重，在0.1%的噪聲強(qiáng)度下就無法對(duì)擾動(dòng)目標(biāo)成像。反投影算法(“BP”)在噪聲強(qiáng)度小于0.1%時(shí)，擾動(dòng)目標(biāo)均可被有效重建成像，即使噪聲強(qiáng)度達(dá)到0.5%，EIT成像才略受影響，但從重建圖像上仍可分辨出擾動(dòng)目標(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個(gè)有限單元)下，進(jìn)行一次圖像重建，反投影算法耗時(shí)不足1秒，而NOSER算法和靈敏度矩陣算法至少需要十幾分鐘、甚至幾十分鐘，主要原因在于它們均須進(jìn)行一個(gè)大矩陣的正則化和求逆運(yùn)算。

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

4 結(jié)論

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