李 佳， 岳士弘， 王亞茹

(天津大學(xué) 電氣自動化與信息工程學(xué)院 自動化系，天津 300072)

0 引 言

醫(yī)學(xué)電阻抗層析成像(electrical impedance tomography,EIT)是當(dāng)前生物醫(yī)學(xué)工程的研究熱點[1]。目前，針對醫(yī)學(xué)EIT技術(shù)的研究包括肺部呼吸成像、乳腺癌診斷、癲癇病灶診斷等[2,3]。在肺癌早期，癌變組織尚未形成結(jié)節(jié)，傳統(tǒng)的影像技術(shù)無法完成準(zhǔn)確診斷，且頻繁使用會對人體造成傷害，但電導(dǎo)率(介電常數(shù))已隨組織的病變發(fā)生變化，可以利用EIT技術(shù)通過對電導(dǎo)率(介電常數(shù))的檢測，完成對早期肺癌的診斷和預(yù)警。

EIT技術(shù)具有“軟場”效應(yīng)和不適定性，導(dǎo)致重建圖像分辨率較低，為提高圖像質(zhì)量，本文根據(jù)人體胸腔的計算機(jī)斷層掃描(computed tomography,CT)掃描圖片，利用計算機(jī)輔助設(shè)計(computer-aided design,CAD)軟件構(gòu)建胸腔器官的結(jié)構(gòu)模型，并結(jié)合器官、組織的電導(dǎo)率等先驗信息構(gòu)建正則化矩陣，引入到擴(kuò)展Kalman濾波(extended Kalman filtering,EKF)的狀態(tài)空間表達(dá)式中，完成肺癌組織的圖像重建。

1 算 法

1.1 電阻抗層析成像原理

電阻抗層析運算過程主要分為正問題和逆問題。圖像重建即屬于逆問題的求解。EIT的線性模型為[4]

ΔU=U′(σ0)Δσ=JΔσ

(1)

式中 ΔU為邊界測量電壓變化量,Δσ為電導(dǎo)率變化量,J為Jacobian矩陣。

1.2 建立含有先驗信息的正則化矩陣

對肺部進(jìn)行圖像重建時，真實的胸腔器官結(jié)構(gòu)信息有利于提高重建圖像精度。因此，根據(jù)人體的CT胸腔掃描圖片，利用CAD軟件，對人體胸腔器官組織進(jìn)行建模，如圖1所示。主要包括心臟、肺、脊髓、脂肪4部分，及其電導(dǎo)率σ(Ω/m)變化范圍[5]。

圖1 人體胸腔模型及器官電導(dǎo)率變化范圍

在模型數(shù)值變化范圍內(nèi)隨機(jī)取值，構(gòu)建500個可能的電導(dǎo)率分布向量σn(n=1,2,3,…,500)作為學(xué)習(xí)集。利用主成分分析(principal component analysis,PCA)法對學(xué)習(xí)集進(jìn)行降維處理。將特征值按從大到小順序排列，前5個最大特征值中λ1=3.7570，與其他特征值相差很大，λ1所對應(yīng)的特征向量w1中已包含了絕大部分胸腔電導(dǎo)率分布的先驗信息。

構(gòu)建正則化矩陣L，使L的零空間與向量w1張成的子空間Sw相等[6，7]，即

L=I-WWT

(2)

式中I為單位向量,W的列向量由特征向量w1構(gòu)成。

1.3 含先驗信息的EKF算法

EKF方法是一種狀態(tài)估計次優(yōu)方法[8]，利用EKF對肺部癌變組織進(jìn)行圖像重建時，需建立狀態(tài)空間模型

σk=Fk/k-1σk+ωk，Uk=Vk(σk)+υk

(3)

式中σk為tk時刻的電導(dǎo)率；Fk/k-1為tk-1～tk時刻的狀態(tài)轉(zhuǎn)移矩陣，一般取單位矩陣；Uk為tk時刻的電壓測量值；Vk為電壓值與電導(dǎo)率的函數(shù)關(guān)系；ωk與υk分別為tk時刻的系統(tǒng)噪聲和測量噪聲，為0均值、互不相關(guān)的白噪聲序列，相應(yīng)的協(xié)方差矩陣分別為Qk，Rk。

對Uk關(guān)于初始電導(dǎo)率σ0進(jìn)行泰勒展開，得線性形式

(4)

EKF基于狀態(tài)誤差協(xié)方差最小的原則對狀態(tài)值進(jìn)行估計，加入含先驗信息的正則化項后，構(gòu)建新的代價函數(shù)[9]

(5)

式中Pk/k-1為誤差協(xié)方差矩陣；σ*為先驗電導(dǎo)率分布；α為正則化參數(shù)，通常由經(jīng)驗獲得；L為正則化矩陣，由式(2)計算得到。

式(5)變形并最小化，得到含先驗信息EKF濾波方程

(6)

2 仿真實驗及結(jié)果分析

仿真條件：酷睿i7 3.4 GHz CPU，12 GB內(nèi)存臺式機(jī)。將人體胸腔模型導(dǎo)入Comsol 3.5a軟件中，進(jìn)行物理條件仿真和有限元剖分，并利用MATLAB R2009a軟件實現(xiàn)正問題的求解和圖像重建。EIT系統(tǒng)共有16個電極，采用相鄰電流激勵、相鄰電壓測量的激勵—測量模式，共獲得208個測量值，胸腔被剖分為1 283個單元。

利用Comsol 3.5a建立3種肺癌模型，如圖2(a)，圖中圓表示癌變組織。仿真時，模型中電導(dǎo)率設(shè)置心臟0.7，肺0.12，脊髓0.006，脂肪0.037，肺癌組織0.2，同時向電壓測量值中添加均值為0，協(xié)方差矩陣R為10-2I208×208的白噪聲，分別采用共軛梯度法(conjugate gradient,CG)、EKF方法和含先驗信息的EKF方法進(jìn)行圖像重建。EKF和含先驗信息的EKF方法的初值按經(jīng)驗選取[10]，分別為電導(dǎo)率初值取0 Ω/m，誤差協(xié)方差矩陣Po初值取5I208×208，另外系統(tǒng)噪聲協(xié)方差矩陣Q取10-2I1 283×1 283。CG算法取迭代500次時重建結(jié)果，EKF算法和含先驗信息的EKF算法均取迭代兩次時的重建結(jié)果，各方法重建圖像如圖2。

圖2 各方法仿真實驗重建圖像

可以看出，對模型1，當(dāng)癌變組織較少時，三種方法均可較為準(zhǔn)確地重建出癌變組織的位置和形狀，其中EKF和含先驗信息的EKF算法的形狀與癌變形狀更為接近。對模型2，CG算法的重建圖像位置開始發(fā)生偏移，且形狀也與真實形狀相差較大，而EKF算法和含先驗信息的EKF算法仍可較為準(zhǔn)確地重建出癌變組織位置，其中含先驗信息的EKF算法重建圖像偽影更小。對模型3，CG算法已發(fā)生嚴(yán)重形變，無法分開左肺中的2個癌變組織，但EKF算法和含先驗信息的EKF算法仍可較為準(zhǔn)確地重建出癌變組織位置，可以分隔開左肺的2個癌變組織，但右肺的成像也與真實形狀相差較大，比較2種方法，含先驗信息的EKF算法的偽影更小。

為定量比較重建圖像的質(zhì)量，采用圖像相對誤差e作為評價參數(shù)

(7)

經(jīng)計算3種模型中含先驗信息的EKF方法圖像的相對誤差最小，說明其重建圖像質(zhì)量最好。其中CG算法與其他2種方法的相對誤差值相差較大，說明CG算法重建圖像質(zhì)量最差。對模型3，癌變組織較多時，含先驗信息的EKF方法較EKF方法的誤差更小，重建圖像質(zhì)量更好。

3 結(jié) 論

利用人體胸腔器官分布結(jié)構(gòu)和器官組織電導(dǎo)率等先驗信息建立正則化矩陣，并將正則化項引入EKF中，進(jìn)行肺癌組織的圖像重建。仿真實驗表明：與共軛梯度算法和傳統(tǒng)的EKF算法相比，含先驗信息的EKF算法得到的重建圖像相對誤差更小，圖像質(zhì)量更好。