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

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

引言

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

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

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

邊界條件為

式中，σ為目標(biāo)內(nèi)部電導(dǎo)率的分布，φ0為給定邊界上的電位，Jn為給定邊界上外加激勵的電流密度（無注入電流時為零）。

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

又因?yàn)?/p>

可以得到

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

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

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

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

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

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

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

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

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

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

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

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

圖2 邊界電位擬合曲面

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

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

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

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

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

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

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

2 算法驗(yàn)證實(shí)驗(yàn)

2.1 仿真實(shí)驗(yàn)

建立了半徑10 cm、高10 cm的圓柱體仿真實(shí)驗(yàn)?zāi)Ｐ?，分別在距底面2、4、6、8 cm的截面處每層都放置16個電極，從而形成16×4的電極陣列，如圖3所示。

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

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

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

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

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

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

進(jìn)行曲線擬合的推導(dǎo)，可得

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

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

為了對算法重建結(jié)果進(jìn)行更好的說明，引入總體誤差［10］，即

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

2.2 物理模型實(shí)驗(yàn)

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

圖6 實(shí)驗(yàn)?zāi)Ｐ虵ig.6 Experimental phantom with target

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

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

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

3 結(jié)果

3.1 仿真實(shí)驗(yàn)結(jié)果

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

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

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

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

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

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

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

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

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

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

3.2 物理模型實(shí)驗(yàn)結(jié)果

用節(jié)點(diǎn)反投影方法對圖6所示的實(shí)驗(yàn)?zāi)Ｐ瓦M(jìn)行重建，結(jié)果如圖9所示。

圖9 實(shí)驗(yàn)?zāi)Ｐ椭亟▓D像。（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é)果很好地將水槽內(nèi)的目標(biāo)物體棒顯示出來，說明節(jié)點(diǎn)反投影方法具有良好的空間分辨能力。

4 討論與結(jié)論

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

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

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

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

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

［2］羅辭勇.基于快速牛頓一步誤差重構(gòu)的電阻抗成像算法和實(shí)驗(yàn)研究［D］.重慶：重慶大學(xué)，2005.

［3］徐管鑫.電阻抗成像技術(shù)理論及應(yīng)用研究［D］.重慶：重慶大學(xué)，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é)點(diǎn)反投影法電阻抗成像研究［D］.天津：河北工業(yè)大學(xué)，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ù)學(xué)辭海（第五卷）［M］.北京：中國科學(xué)技術(shù)出版社，2002：86.

［10］吳煥麗.三維電阻抗成像問題研究［D］.天津：河北工業(yè)大學(xué)，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.