電阻抗斷層成像技術重建方法對比研究

孔令睿，賓光宇，吳水才

北京工業(yè)大學 環(huán)境與生命學部，北京 100124

引言

電阻抗成像（Electrical Impedance Tomography，EIT）是一種根據(jù)場域內阻抗分布變化重建場域內部圖像的成像技術。在測量場周圍設置一定數(shù)量的電極，通過對其中一組電極施加安全范圍內的激勵電流或激勵電壓信號，測量其余電極的電位分布，將測量得到的電位分布重建為電導率分布的圖像。

由于人體在呼吸時肺部的通氣量會產(chǎn)生明顯改變，吸入的空氣不導電，導致空氣的電導率與肺部組織的電導率有明顯差異，肺部區(qū)域的阻抗分布也會發(fā)生顯著變化[1-2]，所以EIT 技術能夠有效描繪肺部通氣區(qū)域的動態(tài)變化[3]，可以應用于肺部呼吸功能的檢查，這也是目前EIT 技術主要應用的領域[4-6]。EIT 技術還具有對人體無害、無創(chuàng)、成像速度快、成本低、可多次重復測量等優(yōu)點，在臨床監(jiān)測領域具有廣闊的應用前景[7-8]。

EIT 成像技術中通常會遇到傳統(tǒng)算法計算復雜、噪聲環(huán)境中重建圖像質量差的問題[9]，臨床應用中EIT 成像技術會面臨更多干擾，如被測者動作導致較大的噪聲、電極接觸面的阻抗變化、計算量大導致成像速度慢等都會對重建圖像質量產(chǎn)生影響[10]，所以在臨床應用中需要對傳統(tǒng)算法進行優(yōu)化。2009 年，Adler 等[11]提出一種依靠數(shù)據(jù)驅動的線性重建方法—GREIT，并在肺部阻抗成像中得到應用。2012 年，Borsic 等[12]設計了一個使用L1 范數(shù)的原始對偶內點法框架用于EIT 圖像重建，證明可以有效處理數(shù)據(jù)異常值，對帶有噪聲的數(shù)據(jù)有更好的重建效果。2014 年，Grychtol 等[13]采用8 只麻醉豬進行實驗，通過給豬肺插管連接至呼吸機控制通氣量，在不同通氣效果下進行EIT的數(shù)據(jù)測量，使用反投影（Back Projection，BP）算法、GREIT 算法、截斷奇異值分解法、一步高斯牛頓法及其幾種變體等十二種2D 重建方法對測量數(shù)據(jù)進行重建，比較各算法重建效果。比較結果表明，十二種方法中的先進優(yōu)化算法比傳統(tǒng)算法效果更好，但傳統(tǒng)算法中的BP 算法在實時成像效果中表現(xiàn)較好。

交替方向乘子法（Alternating Direction Method of Multipliers，ADMM）是一種求解具有可分離的凸優(yōu)化問題的計算框架，其處理速度快、收斂性能好[14]，可以快速并準確地求解EIT 逆問題。本研究設計了三組實驗，選取了EIT 傳統(tǒng)算法牛頓法中的一步高斯牛頓法[15-16]，適用于實時成像的BP 算法[17-18]以及新式優(yōu)化算法GREIT 算法[11]與ADMM 算法做圖像重建。實驗建立仿真模型，在模型內部的不同位置設置圓形目標，添加一定的高斯噪聲，將4 種算法重建出的圖像進行直觀對比，并采用5 種評價參數(shù)[11]對重建結果的質量進行量化評價。根據(jù)重建結果圖像和評價參數(shù)對比驗證ADMM 算法的有效性，并分析幾種算法的優(yōu)劣性。

1 電阻抗成像數(shù)學模型

EIT 問題的求解分為求解正問題及逆問題兩部分。正問題是根據(jù)麥克斯韋方程組建立數(shù)學模型，采用有限元法在已知場域內電導率分布和邊界激勵信號的情況下，求解電勢分布近似解[19]。逆問題則是根據(jù)測得的電壓以及已知的邊界激勵信號求解場域內阻抗分布。對于EIT 問題，邊界處的電壓測量值與阻抗分布變化的對應數(shù)學模型如式（1）所示。

其中，U表示邊界處電極上測量的電壓數(shù)據(jù)；S表示正問題求得的敏感系數(shù)矩陣；σ表示阻抗分布。

逆問題求解時方程的解可表現(xiàn)為σ=S-1U，由于S不是一個方陣，方程無法直接求逆，由電磁場的軟場特性也可得知測量電壓與電導率分布存在非線性關系，且求解方程存在病態(tài)性，因為在實際測量中邊界電壓是在離散點上采樣，電極數(shù)目有限，由邊界電極獲得的信息不足以決定目標內部的阻抗分布，方程為欠定方程，存在無窮多組解，求得的解不唯一且存在不確定性。所以EIT 逆問題是一個高度非線性病態(tài)問題，求解具有一定的難度，逆問題的求解也就是圖像重建方法也成了多數(shù)研究的重點。

由于逆問題是一種非線性問題，傳統(tǒng)的EIT 成像算法通常是將逆問題近似為一種線性問題，構造目標函數(shù)，見式（2）。逆問題求解就是求解該目標方程的最優(yōu)解。

2 電阻抗成像重建方法

2.1 BP算法

BP 算法[17]是一種基于等位線間電位差變化投影的圖像重建方法。算法的基本原理是在電導率變化較小時，在兩條相鄰等位線間的投影區(qū)域內將等位線間的阻抗分布和電勢分布的變化等價為線性關系，利用兩個時刻之間的測量電壓數(shù)據(jù)的差值，進行歸一化然后再沿等位線作反投影，再將投影區(qū)域內的結果疊加起來，得到阻抗值變化的結果，就可以作出差分重建圖像。BP 算法的計算公式如式（3）所示。

其中，A為阻抗分布與電勢分布線性關系系數(shù)；ΔU為電勢分布變化；Δσ為阻抗分布變化；j為測量電極編號。

因為在差分計算的過程中測量數(shù)據(jù)里的部分噪聲和誤差可以消減一部分，所以BP 算法的抗干擾性較好。BP 算法還有重建算法易實現(xiàn)、成像速度快、便于實時成像的優(yōu)點。但BP 算法的成像精度不夠理想，達不到臨床應用的精度要求。

2.2 一步高斯牛頓法

一步高斯牛頓法[15]是一種基于牛頓迭代法衍生的一步快速算法。一步高斯牛頓法的基本思想原理是在假定電導率分布發(fā)生一個較小的擾動的條件下，認為作差分求解后求得的電導率分布與真實的電導率分布之間的誤差極小，可以忽略，因為動態(tài)成像過程中主要是獲得兩個時刻間數(shù)據(jù)的差值，所以放棄牛頓迭代法中的迭代過程，只取一步計算。其目標函數(shù)如式（4）所示。

其中，W為測量電壓值的協(xié)方差矩陣。

為了克服病態(tài)性問題，需要在目標函數(shù)中引入正則化項，使目標函數(shù)的求解有約束后能夠變得收斂穩(wěn)定。加入正則化項的目標函數(shù)表現(xiàn)形式如式（5）所示。

其中，R為正則化矩陣；λ 為正則化參數(shù)，λ 通常根據(jù)經(jīng)驗定義。

對這種形式的目標函數(shù)求解用最小二乘法，目標函數(shù)最小二乘解的形式如式（6）所示。

其中，Δσ為兩個時刻前后電導率變化；J為Jacobian矩陣。

2.3 GREIT算法

GREIT 算法是一種依靠數(shù)據(jù)驅動的重建算法[11]。其目標函數(shù)如式（7）所示。

其中，σ為阻抗分布；U為電勢分布；R為重建矩陣。

GREIT 算法的基本思想是將一系列的仿真數(shù)據(jù)作為訓練目標進行訓練后，得到使重建圖像誤差最小的重建矩陣R，再對其他數(shù)據(jù)進行重建時，直接利用重建矩陣即可得到阻抗分布。求解最優(yōu)重建矩陣的計算公式如式（8）所示。

其中，ε為圖像重建誤差；k為訓練樣本數(shù)；w為加權矩陣，是每一個訓練目標在計算重建矩陣中的權重。

對公式（8）進行求導可得公式如式（9）所示。

由公式（9）可得重建矩陣R如式（10）所示。

2.4 ADMM算法

ADMM 算法是一種將對偶上升法的可分解性與增廣拉格朗日乘子法的收斂性結合起來的算法，在求解分布式凸優(yōu)化問題中得到廣泛應用。在EIT 應用中可以用來求解式（11）：

其中，J為Jacobian 矩陣；F為正則化矩陣；λ 為正則化參數(shù)。

分3 步做迭代求解，見式（12）～（14）。

其中，公式（13）為2 范數(shù)求解問題，求解可得式（15）。

公式（14）為1 范數(shù)正則問題，可以通過軟閾值公式進行求解得式（16）。

其中，Shrink稱為軟閾值算子，其含義如式（17）所示。

軟閾值公式也稱軟閾值收縮公式，因求解的過程可以看作收縮的過程。這種求解思路在很多其他類似的問題上（如硬閾值、TV 正則化項等）都得到應用。

3 實驗設計

為了全面對比不同算法的重建性能，共設計了3 組實驗，實驗基于Matlab R2018a 版本和EIDORS 軟件包完成。第一組實驗是在普通圓場域中的不同位置放置電導率相同、半徑相同的單目標，分別在理想環(huán)境和噪聲環(huán)境下用不同算法進行重建，對比不同算法在場域不同位置處的目標的重建能力。第二組實驗是在普通圓場域中的不同位置放置電導率不同、半徑不同的雙目標，同樣分別在理想環(huán)境和噪聲環(huán)境下進行重建，對比不同算法對雙目標物體的重建能力。第三組實驗是選取胸腔輪廓模型，將成人在呼吸過程中的測量數(shù)據(jù)進行重建，對比不同算法在實際應用中的性能。

第一組、第二組實驗使用的都是半徑為1 的二維圓形有限元模型（圖1），在邊界周圍設有16 個電極，采取相鄰激勵-相鄰測量模式，施加幅值為1 mA 的電流，模型被劃分為46194 個單元，背景電導率設置為1。然后使用軟件包中的公式來計算正問題，選擇一對電極對作為電流的注入電極，隨后測量剩余電極對之間的電壓。最終獲得208（16×13）個邊界電壓值。

圖1 圓形二維有限元模型

第一組實驗在沿圓域水平半徑上0～0.9 區(qū)域內以平均間隔取50 個點作為圓心，依次放置半徑為0.05、電導率為2的目標物體。分別使用GREIT 算法、BP 算法、一步高斯牛頓法、ADMM 算法進行圖像重建。第二組實驗設計了雙目標仿真數(shù)據(jù)重建，來比較不同數(shù)目的目標物體的重建效果。在簡單圓形場域內隨機生成兩個大小不同、電導率不同的目標。目標的半徑大小和位置都是隨機的，半徑的取值范圍為0.05～0.1，其中一個目標的電導率設置為10，另一個目標電導率設置為0.1。兩個目標的位置設置不會使兩個目標完全重合，但可能會有一部分的重疊。第三組實驗設計使用模擬肺部輪廓的有限元模型對成人肺部呼吸實際測量數(shù)據(jù)進行重建，來驗證幾種重建方法在實際應用中的重建效果。數(shù)據(jù)來自EIDORS 軟件包，受試者采取仰臥位，0 號電極放置在受試者背部中心位置，采取相鄰激勵-相鄰測量模式，在受試者正常呼吸時以7 幀/s 的速度采集數(shù)據(jù)，這些數(shù)據(jù)體現(xiàn)了一次完整的呼吸過程，共有34 組。

3.1 評價參數(shù)

研究選取了振幅響應（Amplitude Response，AR）、位置誤差（Position Error，PE）、振鈴（Ringing，RNG）、分辨率（Resolution，RES）、形狀變形（Shape Deformation，SD）5 項評價參數(shù)對重建算法進行量化的評價。

3.1.1 AR

AR 反映了重建圖像成像的穩(wěn)定性，其公式如式（18）所示。

其中，k為一幅圖像像素的個數(shù)，x為圖像像素，Vt為目標的體積，σt為目標的阻抗值，σr為背景區(qū)域的阻抗值，Δσ=σt-σr。

對于AR 的評價標準，在整個場域內任何位置的目標的AR 越恒定，說明重建效果越好。恒定的AR 是評價參數(shù)中最重要的一項，因為當AR 不恒定時，說明目標處于場域內不同位置時將對圖像產(chǎn)生不同的貢獻，將會難以評價重建圖像。

3.1.2 PE

PE 反映了重建圖像目標位置與實際目標位置的誤差。其公式如式（19）所示。

其中，rt為重建圖像中目標到場域重心的距離；rq為xq到場域重心的距離；xq為所有大于1/4 最大振幅的圖像像素，其計算公式如式（20）所示。

PE 的值為正時，說明重建圖像目標位置對比實際目標位置更靠近中心位置。

對于PE 的評價標準，PE 越小就是誤差越小，說明重建圖像越準確；并且在位于不同徑向位置目標的PE 越穩(wěn)定說明重建圖像效果越好。

3.1.3 RNG

RNG 是測量重建圖像中圓C外相反符號的圖像振幅與圓C內圖像振幅之比。其公式如式（21）所示。

其中，圓C為以xq的重心為圓心的圓。對于RNG 的評價標準，RNG 應該盡可能小且均勻，過沖會導致圖像產(chǎn)生不正確的解釋。

3.1.4 RES

RES 反映了重建圖像中目標的大小，與點擴散函數(shù)、的大小相同，其公式如式（22）所示。

其中，Aq=Σk[xq]k，為圓C內像素值的和；AO為重建圖像的所有像素疊加。

對于RES 的評價標準，RES 應該盡可能小且均勻，才能更準確地表示目標阻抗分布的形狀。如果不同位置的目標RES 不均勻，會導致面積較大的目標重建位置不準確。

3.1.5 SD

SD 是測量重建圖像目標的偽影，即重建后的圖像中的目標像素xq不在圓C內的部分占比，其公式如式（23）所示。

對于SD 的評價標準，SD 應該盡可能小且均勻，SD越大說明重建目標位置與原始目標位置誤差越大，說明重建目標偽影大，結果不準確，解釋重建圖像時可能會導致解釋錯誤。

3.2 實驗結果

程序運行的計算機為榮耀MagicBook Pro，運行內存16.0 GB，CPU 為Intel（R）Core（TM）i5-8265U@1.60 GHz，顯卡為Nvidia GeForce MX250。由于GREIT 方法需要先使用仿真數(shù)據(jù)訓練模型，再用訓練好的模型進行圖像重建，但每個模型只需訓練一次，所以對模型的訓練時間單獨計算，模型訓練大約需要21.02 s。每種方法的重建時間如表1 所示。

表1 圖像重建方法運行時間比較

從50 個數(shù)據(jù)中取選10 個進行重建圖像結果展示，結果如圖2 所示。其中ADMM 算法采用的是L1 正則化，一步高斯牛頓法采用的是L2 正則化。可以由圖2 看出4 種算法在無噪聲環(huán)境中都是在場域中心位置的偽影較大，在邊緣位置的偽影較小。

圖2 單目標無噪聲數(shù)據(jù)圖像重建結果

為了更貼近實際情況環(huán)境，測試本文重建方法的抗干擾性，實驗在仿真數(shù)據(jù)中加入12 dB 高斯噪聲，使用帶有噪聲干擾的模擬電位進行圖像重建。重建結果如圖3 所示。由圖3 可知，相比其他算法，ADMM 算法在加入噪聲后重建結果總體看來受到的擾動最小，中心區(qū)域的目標圖像擾動較大，目標的形變比較大，而邊緣區(qū)域的重建效果較好，幾乎與無噪聲環(huán)境中的重建效果相同。而一步高斯牛頓法、BP 算法和GREIT 算法在噪聲環(huán)境中的目標圖像受到噪聲的影響都比較大，抗干擾能力較差。

圖3 單目標12 dB噪聲數(shù)據(jù)圖像重建結果

各重建方法對50 個數(shù)據(jù)的重建評價參數(shù)的對比結果如表2～3 所示。由表2～3 的評價參數(shù)對比可看出，在無噪聲環(huán)境中GREIT 算法的各評價參數(shù)表現(xiàn)相對較優(yōu)，AR 最接近于1，所有參數(shù)都相對較小而且較均勻穩(wěn)定。而對噪聲數(shù)據(jù)的重建，ADMM 算法的PE、RNG 和SD 的值都相對較小，說明ADMM 算法相比其他算法對加入噪聲的數(shù)據(jù)重建結果更準確，具有明顯的抗噪聲干擾優(yōu)勢。

根據(jù)表2 的評價指標，GREIT 算法在無噪聲時RNG、RES 和SD 的均值為0.168、0.281、0.067，均值都是最小的，并且方差為0.005、0.001、0.001，方差較小，變化幅度最穩(wěn)定，AR 為1.095±0.025 也較為接近1，所以在無噪聲情況下GREIT 算法的表現(xiàn)較好。根據(jù)表3 的評價指標，ADMM 算法在加入12 dB 噪聲后PE、RNG 和SD 的均值為0.018、0.405、0.202，相比其他算法最小，方差為0.013、0.066、0.009，相比其他算法是最小的，變化幅度最穩(wěn)定，所以ADMM 算法相比其他算法對加入噪聲的數(shù)據(jù)重建結果更準確。對比表2 和表3 也可看出，在加入12 dB 噪聲后，ADMM 算法的各評價指標相比無噪聲時的評價指標變化最小，變化最大的參數(shù)是AR，在加入噪聲后均值降低了0.242，其他算法的參數(shù)均值變化最大值約為0.8，所以ADMM 算法具有明顯的抗噪聲干擾優(yōu)勢。

表2 無噪聲數(shù)據(jù)重建圖像評價參數(shù)（±s）

表2 無噪聲數(shù)據(jù)重建圖像評價參數(shù)（±s）

方法ARPERNGRESSD一步高斯牛頓法0.732±0.0690.012±0.0010.285±0.0370.356±0.0080.089±0.002 BP算法0.801±0.012-0.009±0.0010.224±0.0690.283±0.0010.082±0.003 GREIT算法1.095±0.0250.019±0.00020.168±0.0050.281±0.0010.067±0.001 ADMM算法0.923±0.0340.010±0.0040.250±0.0030.414±0.0010.111±0.015

表3 加入12 dB噪聲數(shù)據(jù)重建圖像評價參數(shù)（±s）

表3 加入12 dB噪聲數(shù)據(jù)重建圖像評價參數(shù)（±s）

方法ARPERNGRESSD一步高斯牛頓法0.710±0.1200.052±0.0070.627±0.0870.412±0.0200.258±0.015 BP算法1.035±0.0390.168±0.0171.031±0.1240.287±0.0020.536±0.026 GREIT算法1.020±0.0330.203±0.0260.936±0.1620.346±0.0030.512±0.039 ADMM算法0.681±0.0870.018±0.0130.405±0.0660.467±0.0130.202±0.009

除了測試重建算法的抗噪能力，還需要對重建算法的泛化能力進行評估，因此第二組實驗設置了兩個位置隨機的重建目標，兩個目標的電導率分別設置為0.1 S/m 和10 S/m，第二組實驗重建結果如圖4 所示。由圖4 可以看出，雙目標重建結果中BP 算法和GREIT 算法可以將選展示結果中的各個位置的目標全部明確重建出來，一步高斯牛頓法和ADMM 法某些部位的重建結果比較模糊。

圖4 雙目標無噪聲數(shù)據(jù)圖像重建結果

在雙目標數(shù)據(jù)中同樣加入12 dB 高斯噪聲，測試雙目標條件下各重建算法的抗干擾性，重建結果如圖5 所示。由圖5 可以看出，雙目標重建在加入噪聲后一步高斯牛頓法、BP 算法、GREIT 算法都受到了噪聲的干擾，重建目標的振鈴現(xiàn)象增多，目標形變增多，而ADMM 算法受到噪聲干擾最小。

圖5 雙目標12 dB噪聲數(shù)據(jù)圖像重建結果

共生成雙目標數(shù)據(jù)50 個，50 個數(shù)據(jù)重建結果的PE、RNG、RES 和SD 的對比結果如表4～5 所示（AR 僅可對單目標數(shù)據(jù)重建結果進行計算）。由表4～5 的評價參數(shù)對比可看出， ADMM 算法重建結果的位移誤差都相對較小，無噪聲與加入噪聲的重建結果評價參數(shù)變化較小，說明ADMM算法有一定抗噪聲干擾能力。

根據(jù)表4 的評價指標，GREIT 算法在無噪聲時RNG、RES 和SD 的均值與方差的值都相對較小，重建圖像的偽影較少，GREIT 算法對無噪聲數(shù)據(jù)的重建效果較好。根據(jù)表5 的評價指標，ADMM 算法在加入12dB 噪聲后PE 的均值和方差相比其他算法最小，加入噪聲后重建目標的位置是最準確的，而振鈴和形變的表現(xiàn)較差，說明ADMM算法的泛化能力不足。對比表4 和表5 參數(shù)同樣可看出，在加入12 dB 噪聲后，ADMM 算法的各評價指標相比無噪聲時的評價指標同樣變化最小，所以可以再次驗證ADMM算法較強的抗干擾性能。

表4 無噪聲雙目標數(shù)據(jù)重建圖像評價參數(shù)（±s）

表4 無噪聲雙目標數(shù)據(jù)重建圖像評價參數(shù)（±s）

方法PERNGRESSD一步高斯牛頓法0.103±0.0291.358±0.8360.341±0.0100.193±0.017 BP算法0.102±0.0231.473±1.3910.275±0.0020.148±0.020 GREIT0.150±0.0161.258±1.0010.28600.111±0.003 ADMM0.103±0.0252.120±3.4300.365±0.0030.189±0.029

表5 12 dB噪聲雙目標數(shù)據(jù)重建圖像評價參數(shù)（±s）

表5 12 dB噪聲雙目標數(shù)據(jù)重建圖像評價參數(shù)（±s）

方法PERNGRESSD一步高斯牛頓法0.105±0.0291.365±0.8220.343±0.0100.195±0.018 BP算法0.112±0.0241.552±1.4690.278±0.0020.171±0.021 GREIT0.198±0.0231.269±0.6950.304±0.0010.200±0.020 ADMM0.107±0.0242.090±3.0380.365±0.0040.198±0.033

為了驗證幾種重建算法在對實測數(shù)據(jù)的重建中是否有效，采用了在成人呼吸過程中測量的肺部電壓數(shù)據(jù)進行重建，選取了一段完整的呼吸過程，第三組實驗重建結果如圖6 所示。由圖6 可知， GREIT 與ADMM 算法的重建效果較好，成像目標的位置和大小的誤差都相對較小，重建圖像內容解釋準確，重建目標邊界更清晰。

4 討論

在近幾年的相關研究中[20-22]，在比較算法性能時大多數(shù)研究僅采用了仿真數(shù)據(jù)進行實驗，一些采用實測數(shù)據(jù)的實驗也是以模擬環(huán)境進行實驗[20]，為了能夠更準確驗證重建算法的實用性，本文采用了成人呼吸過程中對肺部的測量電壓進行重建。由圖6 可得，對于實測數(shù)據(jù)的重建效果，GREIT 算法的重建結果中肺部輪廓較清晰，偽影少，ADMM 算法對肺部位置和大小的重建都比較準確，兩種算法都可以觀察到呼吸過程中肺部通氣的變化。

圖6 實測數(shù)據(jù)圖像重建結果

但這兩種方法還存在一些需要進一步研究和解決的問題。由于ADMM 算法需要進行迭代運算，相比其他算法運行速度慢了很多。達到最好的重建效果時，需要迭代的次數(shù)更多，用時會更長；并且目前實驗都是基于二維有限元模型進行重建，如果應用于三維模型中，有限元劃分更多導致計算量會更大。而GREIT 算法無法很好地抑制噪聲的干擾。在保持算法精度的前提下，進一步提升ADMM 算法的運行速度和GREIT 算法的抗噪性能是需要繼續(xù)研究改進的問題。

5 結論

本文基于圓形模型和胸腔輪廓模型比較了ADMM 算法、GREIT 算法、一步高斯牛頓法和BP 算法的重建性能，設計了三組實驗。根據(jù)實驗結果可知，ADMM 算法的重建性能表現(xiàn)較為優(yōu)秀，對單目標和多目標物體的面積與位置的重建都比較準確，重建圖像精度高，最突出的特點是抗干擾性能較好，但是重建目標的邊界不夠清晰；而GREIT算法在無噪聲的情況下重建效果比較好，重建目標的均勻性好，并且可以保留邊界。對于肺部功能成像應用來說，這兩種方法都可以做到有效輔助醫(yī)護人員對肺部通氣情況進行判斷，具有一定的臨床應用價值。