基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法

王 琦 陳曉靜 汪劍鳴* 李秀艷 段曉杰 王化祥

1(天津工業(yè)大學電子與信息工程學院，天津 300387)2(天津市光電檢測技術(shù)與系統(tǒng)重點實驗室，天津 300387)3(天津大學電氣與自動化工程學院，天津 300072)

引言

電阻抗成像(electrical impedance tomography,EIT)是臨床使用的一種新的成像工具，它的基本原理是：通過配置于人體體表的一組電極陣列，施加一定頻率的低伏交變安全電流；再通過掃描電極測量電壓數(shù)據(jù)，提取與人體生理、病理狀態(tài)相關(guān)組織或器官的電阻抗信息；經(jīng)數(shù)據(jù)采集單元并行處理后，送至重構(gòu)計算機；運用求解被測物場電磁問題的圖像重建算法，重建被測組織和器官的二維/三維圖像[1-2]。由于EIT具有非侵入性、實時性、設(shè)備簡單、功能性成像的優(yōu)勢，因此在醫(yī)學、地質(zhì)勘探和工業(yè)過程成像領(lǐng)域獲得廣泛應用[3]。

在當代醫(yī)學中，人體肺部 EIT 檢測，對于發(fā)現(xiàn)某些肺部疾病具有十分重要的科學研究意義和臨床應用價值[4-5]。EIT 不僅可以檢測肺通氣能力，而且可提供與人體胸腔生理和病理變化直接相關(guān)的功能圖像。該技術(shù)比應用肺活量計和流量計更方便，更適合于長期的臨床監(jiān)護，成為 ICU 病房中輔助呼吸機進行呼吸監(jiān)測的有效手段。此外，由于EIT 具有無輻射、實時性強及便攜性等優(yōu)點，因此在肺損傷病人早期診斷和床旁實時監(jiān)測方面具有不可比擬的優(yōu)勢[6]。

EIT問題包括正問題和逆問題。正問題是由生物模型內(nèi)部的阻抗分布及邊界激勵信號，求物體內(nèi)部或表面的電壓和電流分布；逆問題是通過表面的電壓、電流分布及邊界激勵信號，求生物體模型內(nèi)部的阻抗分布。圖像重建過程屬于逆問題，它是EIT的核心技術(shù)，圖像重建模型中的靈敏度矩陣構(gòu)建依賴于正問題求解。通常采用有限元剖分方法求解電阻抗成像正、逆問題[7]，其中正問題通過網(wǎng)格剖分獲得場域節(jié)點電勢，逆問題通過網(wǎng)格剖分獲得生物體模型內(nèi)部電導率分布[8]。網(wǎng)格剖分需要已知重建對象的邊界形狀，才能得到更加逼近生物體真實電阻抗分布的剖分模型，進而提高重建圖像的質(zhì)量[9-10]。建立不同邊界形狀的胸部模型會影響成像的精度，因此構(gòu)造符合人體真實邊界形狀的胸部模型是求解EIT問題的前提。目前國內(nèi)外對于胸部電阻抗的研究多用圓形代替胸部結(jié)構(gòu)，Saka提出，在實際EIT應用中，根據(jù)生物體邊界形狀建立橢圓幾何模型，比圓形幾何模型更能準確分析EIT問題[11]。嚴佩敏等用橢圓形模型代替胸部結(jié)構(gòu)進行電阻抗成像[12]，但是由于胸部輪廓具有特異性，用近似人體胸部形狀形成的統(tǒng)一模型[13-14]會引入測量誤差以及成像誤差。因此，獲得正確反映人體胸部真實結(jié)構(gòu)的電阻抗成像方法仍是一個需要探究的問題。

本研究提出一種基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法。通過評價指標——肺部區(qū)域比例(lung regional ratio，LRR)，評價該方法、基于橢圓形模型成像方法、基于圓形模型成像方法的成像效果，并用仿真實驗進行了驗證。結(jié)果表明：基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法能獲得與人體結(jié)構(gòu)更為接近的肺部區(qū)域比例，有效降低傳統(tǒng)模型成像與人體真實胸部結(jié)構(gòu)成像的誤差，提高成像效果，更加真實地反映人體胸部結(jié)構(gòu)。

1 方法

1.1 正問題建模方法

基于先驗信息的正問題建模包括兩部分。第一部分是對屏息時獲得的胸部CT圖像進行圖像處理，得到人體胸、肺部形狀和輪廓相對準確的信息，尤其是邊界信息，為構(gòu)建較準確的胸、肺部模型奠定基礎(chǔ)，圖1即為電極所在中心平面對應的胸部CT圖像。第二部分是將圖像處理后提取的輪廓邊界信息轉(zhuǎn)化為正問題模型。

1.1.1圖像處理

為了準確地從CT圖像中提取有用信息，需要對原始CT圖像進行處理，以突出有用的圖像信息。首先對圖像進行預處理，包括去噪、濾波等；其次通過圖像分割提取出感興趣區(qū)域，該效果會直接影響到重建后模型的準確性；最后對圖像分割后的圖形進行輪廓提取，作為仿真軟件模型的場域輪廓線。具體的圖像處理流程如圖2所示。

圖1 胸部CT圖像Fig.1 Chest CT tomography

圖2 圖像處理流程
Fig.2 Image processing flow chart

1)圖像預處理。為了有效消除邊緣噪聲，提高輪廓線提取精度，本研究采用基于偏微分方程(partial differential equation′s，PDE′s)的圖像方法。目前，以Perona&Malik模型[15]為代表的各向異性平滑方法已在邊緣檢測、圖像增強、圖像分割以及目標識別等領(lǐng)域得到廣泛應用。本研究采用改進的Perona&Malik模型[15]，設(shè)u0(i,j)為原始CT灰度圖像，引入時間變量t∈[0,T]，在求解偏微分方程的迭代過程中，將得到一組逐漸平滑的圖像u(i,j,t)，即

(1)

式中：Gσ是高斯平滑模板；|u|是圖像梯度模；擴散系數(shù)c(·)是關(guān)于|u|的非負函數(shù)，用于控制擴散速度。

2)閾值分割。為了得到肺部及胸部區(qū)域的輪廓線，需要對圖像預處理后的CT圖像進行閾值分割。本研究采用一種基于圖像形態(tài)學中腐蝕、膨脹的閾值分割算法[16]。CT圖像的分割步驟是：首先對預處理后的CT圖像進行二值化處理，然后去除檢查床提取胸部輪廓，再利用去除檢查床的圖像減去胸廓圖像，即可得到肺部輪廓圖像。最后通過對肺部輪廓和胸部輪廓的融合，得到肺部及胸部邊緣。具體的閾值分割流程如圖3所示。

圖3 閾值分割流程Fig.3 Threshold segmentation flow chart

圖4 胸部及肺部輪廓曲線Fig.4 Chest and lung contours

3)提取輪廓線。為了提高輪廓提取的精準度，本研究采用基于二值圖像的輪廓提取算法[17]。設(shè)e(i,j)為二值化后CT圖像的像素點，E(i,j)為應用規(guī)則處理后的CT圖像像素點，判斷當前像素是否為邊界像素的規(guī)則如下：

(2)

E(i,j)=0時，當前像素不是邊界像素，刪除；E(i,j)=1時，當前像素是邊界像素，予以保留。應用此規(guī)則對二值化CT圖像的每個像素進行處理，保留在圖像中的像素為圖像邊界，即可完成輪廓提取，最終得到較為平滑的胸部及肺部輪廓曲線如圖4所示。

1.1.2構(gòu)建模型

圖5 胸部仿真模型Fig.5 Chest simulation model

基于本文第1.1.1節(jié)獲得的胸部及肺部輪廓圖像，構(gòu)建胸部仿真正問題模型。將胸部及肺部輪廓曲線轉(zhuǎn)化為CAD矢量圖，運用有限元仿真軟件Comsol，分別建立基于人體結(jié)構(gòu)先驗信息的胸部仿真模型、基于橢圓形模型的胸部仿真模型和基于圓形模型的胸部仿真模型，進行EIT正問題計算，從單元劃分到計算都可以在短時間內(nèi)得到計算結(jié)果數(shù)據(jù)。需要指出的是，對于圓形仿真模型，其直徑是CT圖片中最高點與最低點坐標之差；對于橢圓形仿真模型，其長軸是CT圖片中X軸最大值與最小值坐標之差，短軸是CT圖片中最高點與最低點坐標之差。3種重構(gòu)模型的電導率是一致的，正問題初始電阻抗分布：胸腔是均勻分布的，不考慮內(nèi)部器官組織。設(shè)置胸腔組織電導率為0.024 S/m，設(shè)置肺部的電導率為0.1 S/m[18]。胸腔組織電導率為0.024 S/m，脂肪電導率為0.024 S/m[18]仿真實驗采用的是相鄰驅(qū)動模式，對于有N個電極的EIT系統(tǒng)，相鄰模式可以采集到N×(N-3)邊界電壓。本研究采用16個電極，胸部仿真模型如圖5所示。

1.2 逆問題剖分方法

逆問題剖分要解決兩個問題：首先，要有一個有效的圖像邊界，該邊界可以通過正問題模型獲??；其次，在該邊界的內(nèi)部進行逆問題剖分。為了避免逆問題風暴(inverse crime)[19]，正逆問題采用不同的網(wǎng)格剖分，本研究的正問題采用三角剖分，逆問題采用矩形剖分，即每個剖分單元為矩形。

1.2.1構(gòu)建n×n的矩形網(wǎng)格

(3)

(4)

圖6 n×n矩形網(wǎng)格Fig.6 n×nrectangular grid

圖7 胸部模型上下兩部分邊界坐標點Fig.7 Chest model of the upper and lower parts of the boundary coordinates

1.2.2邊界線性插值坐標

由本文第1.2.1節(jié)可知，剖分網(wǎng)格可以把模型沿X軸方向做n等分，獲得剖分點橫坐標集合A={x1,x2,x3…xn}，則模型上下兩部分邊界剖分點坐標集合分別為B={(x1,yU1),(x2,yUM2),…，(xn,yUn)}和C={(x1,yD1),(x2,yD2),…，(xn,yDn)}，其中yU1,yU2,…，yUn和yD1,yD2,…，yDn分別是上下邊界剖分點的縱坐標。邊界剖分點縱坐標可根據(jù)已知的邊界剖分點橫坐標A及邊界點坐標(xu,yu)∈Rk1×2和(xd,yd)∈Rk2×2通過插值方法獲得。本研究采用線性插值的方法，上下兩部分邊界剖分點縱坐標均需要通過插值方法得到。以獲得上邊界剖分點縱坐標為例：設(shè)(xui,yui)和(xu(i+1),yu(i+1))是兩個相鄰的邊界點，若兩邊界點間隔s個邊界剖分點，則每個邊界剖分點的縱坐標如下：

yU(m+l)=yui+k(x(m+l)-xui)

(5)

最終得到上下兩部分邊界剖分點坐標集合B和C，線性插值后的坐標如圖8紅點所示。模型內(nèi)部的坐標可通過模型邊界坐標獲得，具體方法見本文第1.2.3節(jié)。

圖8 線性插值后的坐標點Fig.8 Coordinates after linear interpolation

線性插值后的邊界坐標點圍成的輪廓與原圖像輪廓的對比如圖9所示。其中第一行是三個CT樣本，第二行中黑色線代表CT圖像邊界輪廓，紅色點表示線性插值后的坐標點。由圖9可知，線性插值后的邊界坐標點圍成的輪廓與原輪廓基本一致。

圖9 線性插值后的邊界坐標點圍成的輪廓與原輪廓對比Fig.9 Comparison between the original contour and the contour enclosed by the boundary points after linear interpolation

圖10 基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法的有效網(wǎng)格篩選方法Fig.10 An effective mesh selection for chest electrical impedance tomography based on prior information of human body

1.2.3篩選有效網(wǎng)格

將包含在模型邊界點坐標范圍內(nèi)的網(wǎng)格定義為有效網(wǎng)格,定義變量H(Xi,Yi)為網(wǎng)格的有效值，有

(6)

若H(Xi,Yi)的值為1，說明該矩形網(wǎng)格在模型邊界內(nèi)部記為有效網(wǎng)格；若H(Xi,Yi)的值為0，則其為無效網(wǎng)格。

通過對圖9中的樣本1進行邊界獲取及剖分，然后篩選得到有效網(wǎng)格，如圖10所示。其中，灰色網(wǎng)格為有效網(wǎng)格，白色網(wǎng)格為無效網(wǎng)格。將有效網(wǎng)格的坐標數(shù)據(jù)、編號以及數(shù)量保存下來，用于后期的圖像重建。為了更加直觀地顯示網(wǎng)格剖分過程，繪制流程如圖11所示。

圖11 基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法的網(wǎng)格剖分流程Fig.11 Mesh generation flow chart of chest electrical impedance tomography based on prior information of human body

網(wǎng)格精細程度將決定圖像分辨率，即n值越大，圖像分辨率越高。本研究采用n=32的分辨率坐標值，共1 024個矩形剖分單元。圖12(a)為剖分后繪制的矩形網(wǎng)格線，其中每個小矩形就是一個剖分單元，它比較直觀地顯示剖分效果，且能較直觀地描繪胸部模型的形狀。為了對比剖分效果，64×64的矩形網(wǎng)格線如圖12(b)所示，128×128的矩形網(wǎng)格線如圖12(c)所示。由上述剖分網(wǎng)格線可知，在有限元剖分時，隨著網(wǎng)格精細程度的提高，圖像輪廓更加接近真實場域形狀。

圖12 不同精細程度的剖分網(wǎng)格。(a)32×32；(b)64×64；(c)128×128矩形剖分網(wǎng)格Fig.12 Different square mesh.(a)32×32；(b)64×64；(c)128×128

1.3 評價方法

為評估不同方法構(gòu)成的胸部模型與人體胸部真實結(jié)構(gòu)的差異，引入一種基于重建圖像結(jié)果的評價指標——肺部區(qū)域比例(lung regional ratio，LRR)，其定義如下：

(7)

式中，TV表示肺部通氣面積，SV表示胸部截面除肺部以外區(qū)域的面積[20]。

若用F∈Rl×1、E∈Rl×1分別表示肺部區(qū)域像素點的集合、胸部截面除肺部以外區(qū)域像素點的集合(l表示有效的剖分點數(shù))，則集合F、E可分別表示為

(8)

(9)

式中：t(x,y)是像素點(x,y)的像素值；T為肺部區(qū)域阻抗閾值，表示肺部區(qū)域最低阻抗值。

那么，TV和SV可分別表示為

(10)

(11)

本研究采用最大類間方差算法[21]確定區(qū)分肺部區(qū)域和非肺部區(qū)域的閾值。設(shè)f(x,y)為重建圖像的位置(x,y)處的灰度值，具有L個灰度級(本研究中L選為20，即將歸一化的電導率范圍平均分為20份)，灰度級用變量b表示。設(shè)灰度級為bi的所有像素個數(shù)為fi，其中i=0,1,…,L-1，且b0

(12)

選擇第t個灰度級作為圖像像素的閾值，用于區(qū)分胸腔C0和肺部組織C1，此時胸腔C0的灰度級包括b0,b1，…，bt-1，肺部組織C1的灰度級包括bt,bt+1，…，bL-1，則胸腔與肺部組織的分割閾值bt可表示為

bt=arg.max{?2(t)} (0≤t

(13)

式中，?2(t)是類間方差，可表示為

?2(t)=ω0(λ0-λ)2+ω1(λ1-λ)2

(14)

分別計算不同成像方法重建圖像的LRR，與基于CT圖像計算的人體真實LRR進行比較，并計算相對誤差，用來評價本研究提出方法的有效性。相對誤差的計算公式為

(15)

式中,l為成像方法計算得出的LRR，L為基于CT圖像計算的人體真實LRR。

筆者共研究30個肺部健康的人體CT圖像樣本，樣本取自5個年齡段的人群，分別為：20～30，30～40，40～50，50～60，60～70歲；每個年齡段共6個樣本，其中男性、女性各3個。每個樣本分別采用基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法、基于橢圓形模型成像方法、基于圓形模型成像方法進行成像。在重建圖像中，重構(gòu)斷面為電極所在的中心平面，重建的灰度值根據(jù)電導率的真實值歸一化至[1,3]范圍內(nèi)。由于動態(tài)成像方法具有相對簡單、測量誤差小、成像質(zhì)量更高的優(yōu)點[22]，所以本研究采用基于共軛梯度的動態(tài)成像方法來求解電阻抗成像的逆問題。逆問題的剖分及成像采用Matlab2012a軟件，所有的結(jié)果在處理器為Intel(R)Core(TM)i5-3 230 M CPU、內(nèi)存為4GB的計算機上獲得。

1.4 顯著性檢驗

對研究的30個樣本進行顯著性檢驗，根據(jù)實驗的特點使用配對樣本t檢驗的方法[23]，分別對3種成像方法計算的LRR值以及3種成像方法與人體真實LRR值的相對誤差做顯著性檢驗。以3種成像方法算出的LRR值為例，顯著性檢驗具體步驟如下：

1) 建立假設(shè)檢驗，確定檢驗水準。H0表示基于成像方法算出的LRR與真實值不存在顯著性差異，H1表示基于成像方法算出的LRR與真實值存在顯著性差異。雙側(cè)檢驗，檢驗水準α=0.05。

2) 計算統(tǒng)計量為

(16)

3) 確定P值，做出統(tǒng)計推斷。運用統(tǒng)計軟件確定P值：若P>α，則H0假設(shè)成立；若P<α，則H1假設(shè)成立。

2 結(jié)果

用基于不同模型的成像方法，對樣本的圖像進行重建比較。以圖9中的3個樣本為例，3種成像方法重建圖像結(jié)果如圖13所示，其中分別用黑色實線、紅色實線標出了肺部區(qū)域輪廓、胸部區(qū)域輪廓。不同成像方法的統(tǒng)計學分析如表1所示。

圖13 不同成像方法重建圖像結(jié)果Fig.13 Reconstruction of image results with different imaging methods(From the left to the right: sample 1～sample 3)

從圖13的結(jié)果可看出，邊界形狀對成像效果有一定影響。對比3種方法重建圖像的結(jié)果，發(fā)現(xiàn)存在明顯的差異，其中基于人體先驗信息的胸部電阻抗成像方法可獲得更加接近人體胸部的邊界形狀。

由表1可見，對于3種基于不同模型方法成像結(jié)果的LRR均與基于CT圖像計算的人體真實LRR存在差異，但是基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法得到的LRR更加接近真實值。經(jīng)計算可知：基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法與人體真實LRR的相對誤差最小。通過顯著性分析發(fā)現(xiàn)：基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法得到的LRR與基于CT圖像計算的人體真實LRR相比，P>0.05，大于顯著性水平；基于橢圓形模型成像方法、基于圓形模型成像方法與基于CT圖像計算LRR相比，P<0.05；基于橢圓形模型成像方法、基于圓形模型成像方法與基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法的相對誤差相比，P<0.05，均小于顯著性水平。由此說明，基于人體結(jié)構(gòu)先驗信息的成像方法能有效降低成像誤差，提高成像效果。

表1 不同成像方法計算LRR統(tǒng)計學分析(均值±標準差)Tab.1 Statistical analysis of LRR by different imaging method (mean±SD)

注：與LRR真實值相比,a:P>0.05, b:P<0.05; 與基于人體結(jié)構(gòu)先驗信息成像方法的相對誤差相比，c:P<0.05

Note: Compared with real LRR, a:P>0.05, b:P<0.05; Compared with the relative error of imaging method based on priori information of human body structure, c:P<0.05

3 討論

EIT在人體胸部應用中，邊界形狀的不確定性是一個難點，由于胸部輪廓具有特異性，建立不同邊界形狀會影響性能和重建效果[24-25]，所以獲取真實的人體胸部邊界形狀是求解EIT問題的前提。

在胸部電阻抗成像問題中，建立符合人體真實胸部及肺部邊界形狀的模型并進行有效剖分，對成像效果有著積極作用。傳統(tǒng)圓形模型、橢圓形模型成像方法均不能真實反映人體胸部結(jié)構(gòu)，成像誤差較大。所以，如何正確獲取人體胸部及肺部邊界形狀來提高成像效果是一個難點。本研究提出的基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法，通過對CT圖片進行圖像處理提取胸部及肺部輪廓，為正問題和逆問題提供圖像邊界的先驗信息，基于邊界先驗信息進行正問題胸部建模，避免了統(tǒng)一模型引入的測量誤差。同時，基于邊界先驗信息提出一種有效的圖像逆問題剖分方法，在減小誤差的同時提高成像效果，從而使重建圖像形狀更接近真實情況。

研究結(jié)果反映出所提出的基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法的有效性，與傳統(tǒng)基于橢圓形模型成像方法、基于圓形模型成像方法相比，該方法在成像效果上更加接近人體胸部邊界形狀；不同成像方法的評價指標——肺部區(qū)域比例(LRR)的計算結(jié)果表明，該方法可獲得與人體真實結(jié)構(gòu)更接近的肺部區(qū)域比例。所提出方法的LRR與基于CT圖像計算的真實LRR之間的相對誤差(3.71%±1.77%)顯著小于基于橢圓形模型(10.29%±3.30%)和基于圓形模型(12.74%±2.87%)這兩種成像方法。該方法顯著降低了傳統(tǒng)模型成像方法與人體真實胸部結(jié)構(gòu)的成像誤差，成像效果得到明顯提高。由于本研究提出方法的標準差與真實LRR相對誤差的標準差最小，所以該成像方法更穩(wěn)定。分別就這三種成像結(jié)果的LRR值及其與真實值的相對誤差進行統(tǒng)計學分析，結(jié)果表明：基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法與真實的LRR值不存在顯著性差異，基于橢圓形模型成像方法、基于圓形模型成像方法與真實的LRR值存在顯著性差異，基于橢圓形模型成像方法、基于圓形模型成像方法與基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法的相對誤差存在顯著性差異。這是由于胸部輪廓具有特異性，用圓形模型或橢圓形模型代替人體胸部模型會引入測量誤差以及成像誤差?；谌梭w結(jié)構(gòu)先驗信息的胸部電阻抗成像方法可有效解決胸部EIT圖像重建的邊界輪廓特異性問題，提高成像質(zhì)量。

在本研究中，利用屏息狀態(tài)下的CT圖像獲取胸部及肺部的邊界信息。由于呼吸過程中肺部區(qū)域、胸部區(qū)域會隨著呼吸而變化，所以無法獲取確定的輪廓邊界信息。在下一步工作中，將考慮呼吸過程對胸廓形狀的影響，融入形變先驗信息，并且進一步提高剖分精度，改進成像質(zhì)量。

4 結(jié)論

本研究提出一種基于人體結(jié)構(gòu)先驗信息的胸部電阻抗成像方法，通過現(xiàn)有圖像處理方法獲得胸部、肺部輪廓，為正問題和逆問題提供圖像邊界的先驗信息，同時基于邊界先驗信息提出一種有效的圖像逆問題剖分方法，從而使重建圖像形狀更接近人體胸部真實情況，有效降低用圓形、橢圓形代替胸部結(jié)構(gòu)或用近似人體胸部形狀形成統(tǒng)一模型引入的誤差，改善成像效果。