何為，賀中華

重慶大學(xué) 電氣工程學(xué)院，重慶 400044

基于2D電阻抗成像的位置誤差和相對數(shù)量指數(shù)的后處理方法

何為，賀中華

重慶大學(xué) 電氣工程學(xué)院，重慶 400044

本文提出了一種基于2D電阻抗成像的位置誤差和相對數(shù)量指數(shù)的后處理方法，用于多種電導(dǎo)率分布的確定、組織成分的識別以及位置信息的確定。此外，還介紹了一種基于現(xiàn)有兩種重建算法的混合正則化算法，并將三種重建算法重建的結(jié)果進行了比較。結(jié)果證明，這種后處理方法可以提高2D電阻抗成像的位置精度和空間分辨率，使得電阻抗成像適用于臨床檢查和過程監(jiān)測。

電阻抗成像；一步牛頓法；正則化算法；數(shù)量指數(shù)；位置誤差

0 前言

電阻抗成像技術(shù)（Eletrical Impedamce Tomography，EIT）通過在體表放置電極注入激勵電流并檢測電壓以獲得電導(dǎo)率的分布圖像[1]。很多研究表明EIT在醫(yī)學(xué)的應(yīng)用前景廣闊[2-4]。然而，EIT成像的固有缺點包括空間分辨率低、檢測精度低和可重復(fù)性差，妨礙了EIT的臨床應(yīng)用。此外，EIT測量時的微弱干擾可導(dǎo)致完全不同的電導(dǎo)率分布圖像[5]。這種病態(tài)的逆問題通常采用正則化算法來求解，以獲得盡可能接近真實值的電導(dǎo)率分布。

本文提出了一種基于二維EIT的位置誤差和相對數(shù)量指數(shù)的后處理方法，以提高EIT的空間分辨率和位置精度。另外，為求解病態(tài)的逆問題，在現(xiàn)有正規(guī)化算法的基礎(chǔ)上，提出了一種混合正則化算法。由于骨組織的電學(xué)特性相對復(fù)雜，故選擇了一種三種電導(dǎo)率分布的簡化的小腿模型來進行分析，以驗證所提出的后處理方法和混合正則化算法。由于該模型有助于進一步了解骨組織的電性能和邊界電壓的關(guān)系，使得采用EIT診斷骨密度成為可能[6]。

1 方法

1.1 建立模型

本文采用一種簡化的二維人體小腿模型。人體小腿實際上有五種組織，包括皮膚、脂肪、肌肉、骨密質(zhì)和骨髓。在仿真模型中，只考慮到了脛骨和三種組織，即肌肉、骨密質(zhì)和骨髓。該簡化模型及其三種組織電導(dǎo)率分布，見圖1。各組織在20 kHz下的電導(dǎo)率值為皮膚：0.0002 S/m，脂肪：0.0240 S/m，肌肉：0.3448 S/m，骨密度：0.0205 S/m，骨髓：0.0028 S/m[6]。肌肉和骨髓的電導(dǎo)率設(shè)定為常數(shù)。Dinc等[7]表明骨密度從0.3332 g/cm3降至0.1666 g/cm3對應(yīng)骨礦物質(zhì)流失，是典型的骨質(zhì)疏松，因此假定骨密質(zhì)電導(dǎo)率從0.0205 S/m降至0.01025 S/m，以此代表骨質(zhì)疏松。但是這種骨密質(zhì)的電導(dǎo)率變化與典型骨質(zhì)疏松的骨密度變化之間對應(yīng)關(guān)系還不確定。

圖1 含三種電導(dǎo)率分布的簡化的小腿模型。（a）用COMSOL仿真的初始電導(dǎo)率分布；（b）用Matlab仿真的電導(dǎo)率分布，其中圓環(huán)部分（代表骨密質(zhì)部分）電導(dǎo)率發(fā)生改變。

采用全電極激勵方式，每個電極環(huán)繞2D模型的邊界均勻放置。該模型檢測電極有16個，分布在半徑為2.8cm的圓上，則相鄰兩電極之間的角度為22.5°。電導(dǎo)率異常區(qū)域（目標(biāo)區(qū)域）為直徑為2.4cm的外環(huán)和直徑為0.8cm的內(nèi)環(huán)圍繞而成的圓環(huán)形區(qū)域。

1.2 測量策略

在相對驅(qū)動模式下目標(biāo)區(qū)域內(nèi)部有較高的靈敏度，因此采用相對驅(qū)動模式[8]。由徑向相對的電極對注入電流，由除注入電流電極對外的相鄰電極測量電壓。當(dāng)電極數(shù)為16個時，共有192（16×12）個電壓測量值。

1.3 重建算法

為了提高重建圖像的穩(wěn)定性，需要克服隨時間變化的電導(dǎo)率擾動引起的微小的測量電壓變化。邊界測量電壓是電導(dǎo)率分布的函數(shù)，表示為。線性近似表示為，可用低階泰勒級數(shù)表示，其中代表離散的電導(dǎo)率矢量，表示靈敏度矩陣。計算在時間間隔上的電導(dǎo)率差異[5]。

1.3.1 一步牛頓誤差重建算法和正則化算法

在一步牛頓誤差重建算法（NOSER）中，正則化矩陣是一個簡單的對角矩陣[9]，該NOSER算法表示為：

正則化算法的通用表達式為：

1.3.2 混合正則化算法

NOSER重建2D圖像效果比較好，但它對病態(tài)逆問題的信號噪聲很敏感。正則化算法對微弱噪聲不敏感，然而它會使邊界平滑化模糊化。因此，將這兩種重建算法結(jié)合起來，條件數(shù)減小，推測混合正則化算法（Combined Regularizoction，CR）能改善重建圖像質(zhì)量[11]。CR算法可以表示為：

1.4 后處理方法

該后處理方法用于提高2D EIT重建圖像的質(zhì)量。在重建過程中所需對兩個參數(shù)進行分析：位置誤差和相對數(shù)量指數(shù)。

1.4.1 位置誤差

異物的位置定義類似于質(zhì)心，即位于重建圖像中的半振幅集合的“質(zhì)量中心”，在此物質(zhì)的重建異常電導(dǎo)率代替物質(zhì)的密度，半幅度集合是指重建圖像中那些異常電導(dǎo)率大于最大異常電導(dǎo)率的一半的所有元素。重構(gòu)圖像的PE越小，與目標(biāo)對象的中心越接近。異物在成像區(qū)域內(nèi)的位置被定義為如下關(guān)系式[14]

1.4.2 數(shù)量指數(shù)

數(shù)量指數(shù)（Quantity Index, QI），通常也可看作圖像指數(shù)，用于評估范圍內(nèi)由于異物引起的電導(dǎo)率相對變化。

1.4.3 后處理過程

后處理過程的步驟如下：

第一步，基于上述3種不同的重建算法重建2D圖像。

第二步，根據(jù)式（7）～（8）算出2D品質(zhì)因數(shù)和位置誤差，并以2D圖像顯示。

第三步，若重建出電導(dǎo)率變化值小于某個閾值，或者超出某個范圍，式（10）中的值將作為新的重構(gòu)電導(dǎo)率值。的值由重建電導(dǎo)率（）除以相對品質(zhì)因數(shù)值（）所得。

2 數(shù)值仿真

數(shù)值仿真包括模型的建立，采用軟件COMSOL Multiphysics 4.2進行正問題仿真，采用Matlab（版本：7.11.0）進行電壓分析、參數(shù)評估和重建圖像。

2.1 正問題仿真

如圖1所示，最大圓、中間圓和最小圓的半徑分別等于2.8cm、2.2cm和1.8cm，由三個圓分隔開的三塊區(qū)域的電導(dǎo)率分別為0.3448 S/m、0.0205 S/m和0.0028 S/m。電流源注入和電壓測量采用16電極的全電極模型。激勵源電流為5mA，頻率為20 kHz。

2.2 數(shù)據(jù)預(yù)處理

對該模型進行有限元網(wǎng)格剖分，得到3136個三角元網(wǎng)格和1625個網(wǎng)格節(jié)點。仿真模型的正則化參數(shù)λ為9.4277×10-11，NOSER的正則化參數(shù)設(shè)定為0.3。

2.3 仿真結(jié)果

采用三種算法重建2D圖像來進行比較，以此檢驗后處理方法和混合正則化算法是否能獲得更好的重建圖像。2.3.1 條件數(shù)

一個矩陣（H）的條件數(shù)表示該矩陣對噪聲的靈敏度，條件數(shù)越大說明逆問題的病態(tài)性越嚴(yán)重。矩陣（H）的奇異值分解（SVD）中，U和V分別是左、右邊的向量，是由H的奇異值由大到小排列而成的奇異值矩陣。H的條件數(shù)表示為。

CR算法的條件數(shù)比其它兩種重建算法都要小（表1）。因此，推測CR重建的效果更好。

2.3.2 定性評價

進行后處理方法的第1步，分別采用3種重建算法重建出圖像，如圖2所示。雖然所有的重建算法都能顯示相對正確的電導(dǎo)率異常位置，但其邊緣都不太清晰。正則化算法和混合正則化算法重建的電導(dǎo)率異常部分的圖像近似成圓形，如圖2（b）和圖2（c）所示，對應(yīng)電導(dǎo)率變化的骨密質(zhì)部分和背景（肌肉）之間的對比較明顯。三種重建算法的2D重建圖像都比骨密質(zhì)部分的實際尺寸偏小，而在形狀上相似，均是環(huán)形。

在圖2中，三種算法都不能確定骨髓對應(yīng)的那部分圓形區(qū)域電導(dǎo)率是否變化，用深藍色表示，可以理解為電導(dǎo)率異常變化很大，也可以理解為這個圓形深藍色區(qū)域的初始電導(dǎo)率本身就很低。圖2重建的圖像與在正文中描述的只有兩種電導(dǎo)率分布的模型極其相似，故圓形深藍色區(qū)域不能判斷為初始電導(dǎo)率本身就很低的骨髓部分。因此，這種較復(fù)雜的3種電導(dǎo)率分布模型，不能區(qū)分和識別骨密質(zhì)和骨髓的確切位置和邊緣，降低了重建圖像的分辨率。

圖2 三種重建算法的重建結(jié)果。（a）一步牛頓誤差重建算法；（b）正則化算法；（c）混合正則化算。

2.3.3 位置誤差

接著進行后處理方法的第二步，重建出2D位置誤差和相對數(shù)量指數(shù)的圖像，見圖3～4。從圖3可以看出，當(dāng)電導(dǎo)率異常區(qū)域（即目標(biāo)區(qū)域）和檢測電極所在的邊界之間的距離越遠時，位置誤差也就越大，原因是電導(dǎo)率異常區(qū)域離檢測電極所在的邊界越遠，電流密度分布越低，因而算法靈敏度變差，使得位置誤差增大。盡管結(jié)果顯示三種重建算法重建出的不同位置的位置誤差呈現(xiàn)出非線性，但其百分比均低于5％。NOSER算法的最大位置誤差值（用深紅色表示）比其它兩種重建算法的都要小，且NOSER算法產(chǎn)生位置誤差的區(qū)域范圍也要小。

圖3 三種重建算法的位置誤差。（a）一步牛頓誤差重建算法；（b）正則化算法；（c）混合正則化算。

圖4 三種重建算法的相對數(shù)量指數(shù)。（a）一步牛頓誤差重建算法；（b）正則化算法；（c）混合正則化算。

2.3.4 相對數(shù)量指數(shù)

由圖4可看出，三種重建方法的相對數(shù)量指數(shù)在骨髓對應(yīng)的區(qū)域變化顯著。而圖3中盡管顯示出電導(dǎo)率異常區(qū)域出現(xiàn)了較小的位置誤差，而邊界比較模糊。將圖4與圖2和圖3比較，從圖4不僅看出環(huán)狀骨密質(zhì)的重建圖像，更能清晰地看到骨髓的真實輪廓。骨髓的電導(dǎo)率比較低，在這里相對數(shù)量指數(shù)數(shù)值也很小，用深藍色表示，其最小值在-800左右。經(jīng)分析主要有兩個原因影響電導(dǎo)率的分布：第一，骨髓的導(dǎo)電性低，僅為肌肉的0.812％，比骨密質(zhì)也要低得多，因此，注入骨髓的電流非常小，檢測電極很難反映出骨髓的實際電導(dǎo)率。骨髓的電導(dǎo)率實際上不變，但所有的重構(gòu)圖像顯示出變化非常大（用深綠色表示）；第二，從目標(biāo)區(qū)域到檢測電極所在的邊界的距離也影響電導(dǎo)率分布,骨髓區(qū)域內(nèi)部相對數(shù)量指數(shù)不恒定，由中心到邊界逐漸變化。從圖4看出，CR重建的相對數(shù)量指數(shù)深藍色部分區(qū)域最小，說明骨髓對應(yīng)區(qū)域內(nèi)由于異物引起的電導(dǎo)率變化的范圍也最小。

表1 雅可比矩陣的條件數(shù)（=0.3,=9.4277×10-11）

表1 雅可比矩陣的條件數(shù)（=0.3,=9.4277×10-11）

矩陣H雅可比矩陣一步牛頓誤差重建算法正則化算法混合正則化算法計算式條件數(shù)1.568×10202.789×10135.077×10132.532×1012

在進行后處理方法第三步之前，先建立另一個仿真模型，只有兩種電導(dǎo)率分布，如圖5（a）所示，這有助于理解為什么在式（9）中要選擇更新的電導(dǎo)率。在這個模型中，是一個圓內(nèi)全部區(qū)域電導(dǎo)率異常，而不是圖1所示的一個圓環(huán)面內(nèi)電導(dǎo)率異常。重建圖像、位置誤差圖像和相對數(shù)量指數(shù)圖像分別如圖5（b）～（d）所示。將圖5和圖2、圖3比較，可以看出這兩種模型的重建圖像與位置誤差圖像均相似，無法區(qū)分這兩種模型。幸運的是與圖4比較，可發(fā)現(xiàn)相對數(shù)量指數(shù)圖像差異十分顯著。從圖4不僅可以清晰地看到骨髓所對應(yīng)的第三種電導(dǎo)率分布，而且相對數(shù)量指數(shù)的最小值約為-800，而圖5看不到骨髓所對應(yīng)的輪廓線和圖形，并且相對數(shù)量指數(shù)的最小值僅約為-40。因此，很容易區(qū)分圖1與圖5（a）兩種不同電導(dǎo)率分布的模型。如果選擇一個適當(dāng)?shù)拈撝?，例?40，采用式（10）處理并更新電導(dǎo)率，重建電導(dǎo)率分布，就能顯示出真實的電導(dǎo)率分布。

圖5 采用后處理法更新電導(dǎo)率前建模驗證。（a）只有兩種電導(dǎo)率分布的模型，圓形區(qū)域內(nèi)電導(dǎo)率異常；（b）重建圖像；（c）位置誤差圖像；（d）相對數(shù)量指數(shù)圖像，采用混合正則化算法。

后處理方法的最后一步是采用更新后的電導(dǎo)率值再次重建2D圖像。采用三種重建算法的重建結(jié)果如圖6，可以看出，骨髓部分（用深紅色顯示）的電導(dǎo)率沒有改變且形狀保持不變，骨密質(zhì)的重建形狀為環(huán)形，但出現(xiàn)了位置偏差，若以骨髓成像部分作為參照物，位置偏差的大小顯而易見。與圖5（b）中的只有兩種電導(dǎo)率分布的重建結(jié)果比較，差異很顯著，從而很容易分辨出是兩種電導(dǎo)率分布還是三種電導(dǎo)率分布。

圖6 采用后處理方法更新電導(dǎo)率后的重建結(jié)果。（a）一步牛頓誤差重建算法；（b）正則化算法；（c）混合正則化算。

3 討論

與現(xiàn)有的2D 電阻抗成像相比，經(jīng)本文提出的后處理方法重構(gòu)出的2D 電阻抗成像更加接近于真實目標(biāo)，并且在復(fù)雜的電導(dǎo)率分布模型中也能檢測出隱藏的電導(dǎo)率分布。這種方法簡單，計算量小，很容易應(yīng)用到復(fù)雜的電導(dǎo)率分布模型中，比如人的小腿模型。結(jié)果表明，這種后處理方法還能對多種電導(dǎo)率成分進行準(zhǔn)確識別。因此，該方法有較大的應(yīng)用前景。不足的是，經(jīng)這種后處理方法重建的2D電阻抗圖像仍然顯示出明顯的位置偏差，而且重構(gòu)目標(biāo)的邊緣不清晰。在今后的研究中，要根據(jù)位置偏差的有利信息進行位置校正，以獲得更接近真實的EIT圖像。

4 結(jié)論

盡管本文提出的CR的條件數(shù)比其它算法更小，但三種重建算法的重建圖像看來，混合正則化算法的優(yōu)勢不明顯?？傮w來說，EIT圖像重建是一個非線性和病態(tài)的逆問題。逆問題的屬性決定了解的不穩(wěn)定性。下一步的研究將嘗試采用其它重建算法，如總變差正則化算法[15]，以改善成像質(zhì)量?？傊?，以后的研究方向應(yīng)不斷優(yōu)化和提高EIT重建圖像的空間分辨率、穩(wěn)定性。

致謝

本文工作得到了中國國家自然科學(xué)基金（編號：51377186）和973計劃（編號：2014CB541600）的資助。

[1]Holder D.Electrical tomography for industrial applications electrical impedance tomography：Methods, history and applications[J].Med Phys Biomed Eng,2004,295-347.

[2]Lee E,Ts ME,Seo JK,et al.Breast EIT using a new projected imagereconstruction method with multi-frequencymeasurements[J].Physiol Meas,2012,33（5）：751-765.

[3]Stojadinovic A,Nissan A,Gallimidi Z,et al.Electrical impedance scanning for the early detection of breast cancer in young women：preliminary results of a multicenter prospective clinical trial[J].J Clin Oncol,2005,23（12）：2703-2715.

[4]Gao J,Yue S,Chen J,et al.Classification of normal and cancerous lung tissues by electrical impendence tomography[J].Biomed Mater Eng,24（6）：2229-2241.

[5]Polydorides N.Image Reconstruction Algorithms for Soft-Field Tomography[D].University of Manchester Institute of Science and Technology,2002.

[6]Katz S,Zlochiver S,Abboud S.Induced current bio-impedance technique for monitoring bone mineral density--a simulation model[J].Ann Biomed Eng,2006,34（8）：1332-1342.

[7]Din? H,Sadiko?lu Y,Savci G,et al.Bone mineral density measurement by quantitative computed tomography in a normal Turkish population[J].Eur J Radiol,1995,21（2）：79-83.

[8]Polydorides N,Aghasi A,Miller EL.High-order Regularized Regression in Electrical Impedance Tomography[J].SIAM J Imaging Sci,2012,5（3）：912-943.

[9]Cheney M,Isaacson D,Newell JC,et al.NOSER：An algorithm for solving the inverse conductivity problem[J].J Imaging Syst Technol,1990,2（2）：66-75.

[10]OH S.Compensation of shape change artifacts and spatiallyvariant image reconstruction problems in electrical impedance tomography[D].University of Florida,2009.

[11]He W,Li B,Xu Z,et al.A combined regularization algorithm for electrical impedance tomography systemusing rectangular electrodes array[J].Biomed Eng-App Bas C,2012,24（4）：313-322.

[12]Hansen PC,Jensen TK.An adaptive pruning algorithm for the discrete L-curve criterion[J].J Comput Appl Math,2007,198（2）：483-492.

[13]Nasehi-Tehrani J,McEwan A,Jin C,et al.L1 regularization method in electrical impedance tomography by using the L1-curve[J].Appl Math Model,2012,36（3）：1095-1105.

[14]He W,Ran P,Xu Z,et al.3D electrical impedance tomography represented by reconstructed planes in a semispherical electrode array model[J].Int J Appl Electrom, 2013,41：443-446.

[15]Hao L,Li G,Xu L.Magnetic detection electrical impedance tomography with total variation regularization[J].Biomed Mater Eng,2014,24（6）：2857-2864.

A Post-Process Method for 2D-EIT-Based Position Error and Relative Quantity Index

HE Wei, HE Zhong-hua
School of Electric Engineering, Chongqing University, Chongqing 400044, China

This paper proposed a post-process method for 2D-EIT-based（Two Dimensional Electrical Impedance Tomography）PE（Position Error）and relative QI（Quantity Index）, which was intended to identify the distribution of multiple electrical conductivity, organ compositions as well as position information.Additionally, an integrated TR（Tikhonov Regularization）based on two present construction algorithms was also introduced and compared.According to the results, this post-process method could improve the placement accuracy and spatial resolution of 2D EIT, which made EIT applicable for clinical examinations and process monitoring.

electrical impedance tomography；newton’s one-step error reconstruction；tikhonov regularization；quantity index；position errors

TM934.7

A

10.3969/j.issn.1674-1633.2015.07.003

1674-1633（2015）07-0008-05

2015-06-28

國家自然科學(xué)基金（51377186）；國家973計劃（2014CB541600）。

作者郵箱：hzh186@cqu.edu.cn