王宏斌 徐桂芝 張 帥 張劍軍 李 穎 郭 磊 顏威利

（河北工業(yè)大學電磁場與電器可靠性省部共建重點實驗室，天津 300130）

引言

電阻抗成像技術通過簡單易行的測量手段，由邊界電位信息重構出目標區(qū)域內的阻抗分布情況，即由較少的先驗信息達到獲得豐富內部功能信息的目的。這使得EIT技術具有低廉、簡易等優(yōu)勢，同時在兼具無創(chuàng)、無放射性的基礎上，可以進行其他成像手段所無法實現(xiàn)的功能成像。不過，其逆問題的欠定、非線性以及嚴重病態(tài)增添了算法研究的困難，各國研究人員提出了許多獨特、巧妙的方法克服障礙，以實現(xiàn)阻抗成像。美國 Rensselaer Polytechnic Institute將這些方法歸納為 3種［1］：一是求解線性反問題，適用于求解場域內電阻率（電導率）變化微小的情況，典型算法如牛頓一步誤差重構算法（NOSER）［2］和反投影方法［3］；二是通過迭代的方法求解非線性逆問題，如牛頓迭代法［4］；三是直接進行逆問題求解，如D-BAR算法［5］。本研究所描述的方法屬于第一類，在傳統(tǒng)反投影方法的投影疊加思想上舍棄了對剖分單元的計算，選擇節(jié)點為計算對象實現(xiàn)阻抗成像。這樣做減少了運算量，提高了成像速度，在三維圓柱體模型的仿真實驗和水槽物理模型的實驗中取得了較好的效果。

1 節(jié)點反投影方法

對EIT問題求解場域的數(shù)學描述用Laplace方程，有

邊界條件為

式中，σ為目標內部電導率的分布，φ0為給定邊界上的電位，Jn為給定邊界上外加激勵的電流密度（無注入電流時為零）。

區(qū)域內的電流場任意位置的電流方向總是與等位線垂直，因此可以將EIT系統(tǒng)看作是一個硬場。由于電流場滿足似穩(wěn)條件，電場E可由標量電位函數(shù)φ來表示，即

又因為

可以得到

從式（6）可以看出，目標內電阻率 ρ（ρ=σ－1）與電位梯度成正比。當電阻率從ρ0變?yōu)棣?時，其變化量為

即阻抗的變化與電位梯度的變化成正比。

在二維情況下，邊界上各個電極測得的電位數(shù)據(jù)是一維的。由于通過等位線的電流總和是一定的，那么可以等效地假定存在一個電流源，電阻率ρ分布即為通過各等位線 l上各節(jié)點的電阻率總和，即

于是，可以將圓域的目標模型映射到一維情況中［6］，如圖1所示，其中x為電極所在邊界位置。

圖1 EIT二維圓域模型的一維電流等效模型Fig.1 1D equivalent circuit model for 2D EIT circle model

假設等位線的等效電阻率與邊界位置存在投影關系，其函數(shù)表示為

由此得到圖1所示等效電路的函數(shù)，即

從而由式（7）得到一維映射等效模型簡化形式為

這樣，就可以通過邊界電位的變化情況，得到區(qū)域內部電阻率的變化情況。將這一變化沿電阻率均勻分布時確定的投影路徑，反投影至所對應的區(qū)域內部等位線上的各節(jié)點，通過多次激勵測量和疊加運算，得到區(qū)域內部的阻抗變化情況，從而實現(xiàn)了動態(tài)阻抗成像。因此，該方法得名為節(jié)點反投影方法（node back projection algorithm，NBPA）［7］。

在三維情況下，圓柱體邊界展開后為一個二維平面，由式（7）展開后其二維映射等效形式為電位變化量對二維坐標的導數(shù)（即梯度形式），故仍由式（7）形式表示為宜。

由于節(jié)點反投影方法需通過相鄰激勵模式在場域內形成投影域，因而其邊界電壓分布為一單調曲線。針對節(jié)點反投影方法對邊界坐標的函數(shù)關系的單調特點，多項式曲線擬合方法［8］通過有限的測量數(shù)據(jù)，擬合出邊界電位曲線。

在三維圓柱體模型中，采用4層同時激勵的相鄰激勵模式，即每次激勵4層中編號相同的電極對，依次進行1-2，2-3，…，16-1的4層同時激勵。這樣，使圓柱體邊界上電壓分布在應用最小二乘法擬合時情況較為特殊。將圓柱體側面展開為二維xoy平面，取電極處電壓為z軸上的坐標值，建立直角坐標系。4層電極上電壓數(shù)據(jù)構成的平面近乎平行于xoy平面。由最小二乘法計算擬合曲面，多項式包含z項時無法通過改變階數(shù)將邊界電壓數(shù)據(jù)擬合為一適合平面。因此，考慮了4層同時激勵模式的特殊性，在擬合多項式中舍去了 z項，得到擬合曲面，如圖2所示。

圖2 邊界電位擬合曲面

式中，n為擬合多項式的預設階數(shù)，x為邊界位置坐標，pi（i=1，2，…，n+1）為各項系數(shù)。

由式（12），可以計算出內部各節(jié)點電位對應邊界電位在擬合曲線上的坐標。當區(qū)域內部阻抗發(fā)生變化時，邊界電位亦發(fā)生變化。節(jié)點反投影方法即由此來實現(xiàn)逆問題求解。

在區(qū)域內阻抗發(fā)生變化的t1時刻測得邊界電位，經曲線擬合后得到多項式為

式中，n為擬合多項式預設階數(shù)，x為邊界位置坐標，pi（i=1，2，…，n+1）為各項系數(shù)。

將式（12）和式（13）代入式（11），可得

由此，可以得到區(qū)域內部各個節(jié)點處的電阻率變化情況，從而實現(xiàn)圖像重構。

由文獻［9］中對先驗信息的定義，在由正問題計算投影位置時，利用邊界電壓在4層同時相鄰激勵時的單調性，對處于已知電壓的電極之間的邊界點進行了擬合近似，因此可以說通過邊界電壓擬合，豐富了先驗信息。

2 算法驗證實驗

2.1 仿真實驗

建立了半徑10 cm、高10 cm的圓柱體仿真實驗模型，分別在距底面2、4、6、8 cm的截面處每層都放置16個電極，從而形成16×4的電極陣列，如圖3所示。

圖3 三維圓柱體仿真模型。（a）側視圖；（b）截面圖Fig.3 3D cylindrical simulation model.（a）lateral view；（b）section

應用有限元方法（finite element method，F(xiàn)EM）進行三棱柱剖分，自下而上均勻分為10層，每個截面剖分為9層。在圓柱體模型中，剖分單元5 040個，節(jié)點3 311個。

在對圓柱體模型進行仿真時，將測得的16組、4層電極上的電壓數(shù)據(jù)，用節(jié)點反投影方法進行處理，得到4層截面內部阻抗信息的重建圖像。

在仿真實驗中，對式（14）進行處理后應用于編程可簡化計算。由于外部激勵的電流是恒定的，即一維等效電路中的電流是恒定的，因此可以假定從而使得式（11）簡化為

通過式（15）的假定得到式（17），在一定程度上簡化了編程。

在模型中放置目標Ⅰ，設定其電阻率為3 Ω·m，圓柱體區(qū)域內其他部分的介質電阻率為2 Ω·m，如圖4所示。

進行曲線擬合的推導，可得

圖4 目標Ⅰ的仿真模型。（a）側視圖；（b）俯視圖；（c）平視圖Fig.4 Simulation model of target.（a）lateral view；（b）top view；（c）front view

為了考察節(jié)點反投影方法對三維空間的適應性，在圓柱體模型中放置目標Ⅱ，如圖5所示。設定目標Ⅱ的兩物體電阻率為3 Ω·m，圓柱體區(qū)域內其他部分的介質電阻率為2 Ω·m。目標Ⅱ與目標Ⅰ的不同之處在于其兩物體的水平高度不同，在豎直方向上形成了一部分交錯。左側目標柱體底部距離圓柱體模型底部1 cm，高4 cm，即其跨越自下而上的第一、二層電極所在平面；右側目標柱體底部距離自下而上第二層電極所在平面1 cm，高4 cm，即其跨越自下而上第三、四層電極所在平面。如此，除自下而上第二層和第三層電極平面之間外，各層電極所在平面上僅存在單一的目標物體。

為了對算法重建結果進行更好的說明，引入總體誤差［10］，即

式中，ρc（i）為計算獲得重建電阻率分布，ρt（i）為設定電阻率分布，m為電阻率分布矢量維數(shù)。

2.2 物理模型實驗

實驗室設計了128通道（4層電極，每層16個復合電極）的 EIT實驗系統(tǒng)［11］，利用節(jié)點反投影方法進行物理模型實驗。實驗水槽直徑30 cm，高30 cm，如圖6所示。在其邊界上的16×4電極陣列中，采用的是4 cm×2 cm的矩形銅片電極，其激勵與測量部分之間由絕緣子隔離。在豎直方向上，相鄰兩個電極相距2 cm。水槽底部距離第1層電極下沿所在平面4 cm，水槽頂部距離第4層電極上沿所在平面4 cm。

圖6 實驗模型Fig.6 Experimental phantom with target

在實驗模型中盛滿 NaCl溶液（電導率約為0.08S/m），然后將不同電導率的材料（環(huán)氧樹脂棒等）放入模型中進行實驗。采用4層同步激勵模式且同時測量，各層均采用相鄰激勵-相鄰測量模式。每次激勵施加于各層之上的電流幅值小于0.5 mA，頻率為100 kHz。由于信號頻率較低，故虛部影響可以忽略，僅考慮阻抗實部進行計算。

在實際操作中，以未放入任何目標時模型內NaCl溶液的電導率分布作為參考數(shù)據(jù)。在放入具有不同電導率的材料之后，系統(tǒng)進行數(shù)據(jù)采集，將這些數(shù)據(jù)與參考數(shù)據(jù)輸入計算機，由所編寫的節(jié)點反投影成像程序進行斷層圖像重建。

將長度為18 cm、直徑為2.5 cm的透明樹脂棒垂直浸入模型，直達底部，樹脂棒上部距離第4層（與仿真模型一致，電極依舊自下而上分布在1～4層）電極所在平面6 cm，

3 結果

3.1 仿真實驗結果

由式（17）對圖4所示目標Ⅰ進行節(jié)點反投影圖像重建，自下而上的各截面結果如圖7所示。

在圓柱體模型中，目標Ⅰ兩物體均為與模型等高的柱體，且其位置中心對稱。圖7（a）是節(jié)點反投影方法的成像結果，相對于圖7（b）傳統(tǒng)反投影方法的重建圖像，不僅很好地顯示出了目標Ⅰ兩物體的形狀和位置，而且偽影也較少。圖7（b）中目標物體附近的偽影擴散較為嚴重，邊界較模糊，而圖7（a）中目標物體附近沒有擴散現(xiàn)象，邊界圓潤，與設定中的目標物體更為相似。同時，模型內距目標物體較遠的區(qū)域中圖像幾乎沒有雜亂的偽影，這與圖7（b）中的偽影分布有很大不同，說明成像效果有所提高。實驗結果表明，節(jié)點反投影方法可以通過電極激勵截面測得的電位信息，對區(qū)域內部的多個三維體進行重建，而且效果較傳統(tǒng)反投影方法有所提高。

圖7 目標Ⅰ的仿真重建圖像（每列自下而上依次為第1～4層）。（a）節(jié)點反投影方法的成像結果；（b）反投影方法成像結果Fig.7 Reconstruction of target Ⅰ（bottom-up section 1～4）.（a）reconstruction images by NBPA；（b）reconstruction images by BPA

用節(jié)點反投影方法對圖5所示目標Ⅱ進行圖像重建，各截面結果如圖8自下而上所示。

在圖8中，兩種方法均清晰地顯示出了目標Ⅱ中的兩物體，但豎直方向上的差異體現(xiàn)得略有不同。圖8（b）的傳統(tǒng)反投影方法重建圖像依舊存在目標物體附近偽影擴散、形狀與設定差距稍大的問題，同時在分辨豎直方向上目標物體的差異時敏感性稍遜。在圖8（a）中，第2層上右側目標物體和第3層上左側目標物體的圖像更為清晰。節(jié)點反投影方法可以更好地重建目標Ⅱ中兩個物體的空間結構，說明其具有較好的三維空間識別能力。

在圖8中，兩種方法重建圖像均存在對第1、2層右側目標物體和對第3、4層左側目標物體的顯示痕跡，這是因為在模型內部放入電阻率不同的目標物體后，改變了電流場的分布情況。由于電流場的軟場特性，發(fā)生變化的范圍會影響到其他層上，從而在跨層截面上出現(xiàn)對目標的反應，是正?，F(xiàn)象。式（7）正是這種現(xiàn)象的數(shù)學表達。而節(jié)點反投影方法分辨目標物體的敏感性，使其在重建圖像中此部分偽影稍濃。如將兩種方法相結合，則可取長補短地解決分辨敏感性與偽影的矛盾。

圖8 目標Ⅱ的仿真重建圖像（每列自下而上依次為第1～4層）。（a）節(jié)點反投影方法的成像結果；（b）反投影方法的成像結果Fig.8 Reconstruction of target Ⅰ（bottom-up section 1～4）.（a）reconstruction images by NBPA；（b）reconstruction images by BPA

用式（18）所述的總體誤差對兩種方法的重建圖像進行比較，結果如表1所示。

表1 節(jié)點反投影方法與傳統(tǒng)反投影方法總體誤差比較Tab.1 Comparison of NBPA and BPA in total error

從表1中可以看出，節(jié)點反投影方法的重建結果比傳統(tǒng)反投影方法重建結果分別降低總體誤差8.87％和6.85％。

3.2 物理模型實驗結果

用節(jié)點反投影方法對圖6所示的實驗模型進行重建，結果如圖9所示。

圖9 實驗模型重建圖像。（a）第1層；（b）第2層；（c）第3層；（d）第4層Fig.9 Reconstruction images of experimental model.（a）section 1；（b）section 2；（c）section 3；（d）section 4

圖9中的重建結果很好地將水槽內的目標物體棒顯示出來，說明節(jié)點反投影方法具有良好的空間分辨能力。

4 討論與結論

節(jié)點反投影方法選擇了節(jié)點作為計算對象，這與大多數(shù)方法針對剖分單元進行計算的方式不同。一般情況下，有限元方法剖分所得到的節(jié)點數(shù)遠小于單元數(shù)。選擇節(jié)點為計算對象，降低了計算中矩陣的階數(shù)，減少了運算量，大大降低了對計算機硬件和數(shù)學軟件的要求，在節(jié)省運算時間的基礎上相對提高了計算精度。

在EIT問題中，內部阻抗的變化會引起邊界電位的變化。相對區(qū)域內豐富的阻抗信息，有限的電極數(shù)所測得的電位明顯單薄。邊界電位的擬合，增加了可用的信息。通過推導，得到邊界電位與內部節(jié)點電位的關系，使得邊界電位擬合不單起到增加先驗信息的作用。將離散的電極（文中采用的模型為每周16個）信息擴展為一條曲線，為數(shù)學計算架起了橋梁。

仿真實驗顯示，節(jié)點反投影方法在三維模型中具有一定的適應性。模型中所放置的多個物體的重建圖像，基本反映了其形狀和位置。與傳統(tǒng)反投影方法相比可以看出，節(jié)點反投影方法具有更好的成像效果、更為準確的目標物體形狀重建能力。水槽模型實驗印證了節(jié)點反投影方法在三維物理模型中的適應能力，重建圖像中目標物體處的邊界柔和、偽影不多，為將來進行去除偽影的工作提供了良好的條件。重建圖像尤其體現(xiàn)了節(jié)點反投影方法具有一定的三維空間識別能力，較為準確地重構了目標物體在豎直方向的形狀和位置。其存在的誤差分析推測為：目標物體的存在，使電流走向發(fā)生了偏轉。針對由此造成的空間分辨問題也是下一步的研究重點之一，與偽影消除有著緊密的聯(lián)系。

本算法是一種動態(tài)反投影算法，需要首先采集背景邊界電壓的分布情況，方可進行差分成像。因此，與其他動態(tài)算法一樣，適合于對人體肺部成像。對于人體內部電阻率不發(fā)生周期性變化的部位，如頭部、手臂等，該方法則無法通過多次測量求平均來得到與實時變化有差異的背景電壓，因而不能進行差分成像。此外，節(jié)點反投影方法在推導中的一些假設可以更加貼近實際情況，有些地方可以通過更加完善的數(shù)學理論來加以雕琢。研究表明，不同層間電極信息與節(jié)點的關系為下一步的研究重點。綜上所述，節(jié)點反投影方法具有深入挖掘的潛力。

［1］Siltanen S.Electricalimpedance tomography and Faddeev Green’s functions［D］.Helsinki：Helsinki University of Technology，1999.

［2］羅辭勇.基于快速牛頓一步誤差重構的電阻抗成像算法和實驗研究［D］.重慶：重慶大學，2005.

［3］徐管鑫.電阻抗成像技術理論及應用研究［D］.重慶：重慶大學，2004.

［4］Cheney M，IsaacsonD，NewellJC.Electricalimpedance tomography［J］.SIAM Review，Society for Industrial and Applied Mathematics，1999，41（1）：85-101.

［5］Siltanen S，Jennifer M，Isaacson D.An implementation of the reconstruction algorithm of A Nachman for the 2D inverse conductivity problem［J］.Inverse Problems，2000，16（3）：681－699

［6］張劍軍.節(jié)點反投影法電阻抗成像研究［D］.天津：河北工業(yè)大學，2008.

［7］Zhang Jianjun，Yan Weili，Xu Guizhi，et al.A new algorithm to reconstructEIT images：node-back-projection algorithm［C］//Proceedings-29th Annual International Conference of the IEEE-EMBS.Piscataway：IEEE，2007：4390－4393.

［8］Zhang Jianjun， Xu Guizhi， Zhao Quanming， et al. Using polynomial curve fitting method to improve image quality in EIT［C］//Proceedings-28th Annual International Conference of the IEEE-EMBS.Piscataway：IEEE，2006：6769-6772.

［9］程民德，何思謙，張景中，等.數(shù)學辭海（第五卷）［M］.北京：中國科學技術出版社，2002：86.

［10］吳煥麗.三維電阻抗成像問題研究［D］.天津：河北工業(yè)大學，2009.

［11］Xu Guizh，Wang Renping，Zhang Shuai，et al.A 128-electrode three dimensionalelectrical impedance tomography system［C］//Proceedings-29th Annual International Conference of the IEEE-EMBS.Piscataway：IEEE，2007：4386－4389.