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

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

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

3 海軍總醫(yī)院全軍航海航空醫(yī)學專科中心，北京，100048

0 前言

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

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

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

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

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

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

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

1.1 反投影算法

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

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

1.2 牛頓一步重建法

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

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

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

1.3 靈敏度矩陣重建算法

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

Geselowitz的阻抗電場理論為：

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

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

2 仿真研究方法

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

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

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

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

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

3 結果及討論

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

3.1 圖像分辨率的比較

從圖2可見，三種算法對于徑向不同位置的目標均能進行有效重建(前三列)；對于位置相距較遠的多目標，其重建圖像均能反映出目標的位置信息(后兩列)；但對于相距較近的目標，重建圖像僅能顯示一個目標，無法分辨出實際目標個數(shù)(第四列)。相比較而言，NOSER算法的重建圖像分辨率高，目標區(qū)域突出，靈敏度矩陣算法次之，但在二者圖像的目標區(qū)域附近存在明顯亮斑偽影。反投影算法的目標區(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ù)是通過計算機模擬計算得到的，它們不含有任何的誤差和噪聲，是一種理想的測量數(shù)據(jù)。在實際應用中，測量數(shù)據(jù)不可能不耦合一定的噪聲。因此，在一定噪聲背景下研究算法的重建效果很有必要。以下按公式(6)為邊界測量數(shù)據(jù)添加噪聲,取噪聲強度系數(shù)分別為0、0.01%、0.05%、0.1%和0.5%，其中0表示不添加噪聲，如圖3各列所示。

在圖3中，NOSER算法(“NS”)的重建圖像在噪聲強度達到0.05%時，其圖像分辨率已顯著下降；噪聲強度達到0.5%時，從重建圖像上已無法分辨目標，重建失敗。靈敏度矩陣算法(“SE”)受噪聲影響更為嚴重，在0.1%的噪聲強度下就無法對擾動目標成像。反投影算法(“BP”)在噪聲強度小于0.1%時，擾動目標均可被有效重建成像，即使噪聲強度達到0.5%，EIT成像才略受影響，但從重建圖像上仍可分辨出擾動目標。因此，就抗噪性能而言，反投影算法最佳，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內存)和剖分規(guī)模(512個有限單元)下，進行一次圖像重建，反投影算法耗時不足1秒，而NOSER算法和靈敏度矩陣算法至少需要十幾分鐘、甚至幾十分鐘，主要原因在于它們均須進行一個大矩陣的正則化和求逆運算。

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

4 結論

以上仿真研究結果表明：反投影算法的重建圖像雖然圖像分辨率不高，但其圖像平滑，無明顯偽影，且該算法具有較強的抗噪性能和快速的運算能力，在臨床檢測的強噪聲背景下，應優(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]劉國強.醫(yī)學電磁成像[M].北京: 科學出版社,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]羅辭勇.基于快速牛頓一步誤差重構的電阻抗成像算法和實驗研究[D].重慶大學,博士學位論文,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.