李彥東 楊 濱 徐燦華 代 萌 付 峰 董秀珍

電阻抗成像(electrical impedance tomography，EIT)是一種無創(chuàng)的醫(yī)學(xué)成像技術(shù)，其基本原理是根據(jù)人體內(nèi)不同組織在不同的生理、病理狀態(tài)下具有不同的電阻(導(dǎo))率，通過各種方法給人體施加小的安全驅(qū)動電流(電壓)，在體外測量響應(yīng)電壓(電流)信號以重建人體內(nèi)部的電阻率分布或其變化的圖像[1]。與傳統(tǒng)醫(yī)學(xué)成像技術(shù)CT、MRI等相比，EIT具有成本低、無創(chuàng)無損傷以及能夠?qū)颊哌M行長時間實時監(jiān)護等優(yōu)點。因此，EIT在腦部持續(xù)性、遲發(fā)性出血類疾病中可展現(xiàn)出良好的臨床應(yīng)用前景[2-6]。

由于EIT的逆問題具有嚴重的不適定性，圖像重構(gòu)矩陣具有很大的條件數(shù)，同時由于系統(tǒng)噪聲和建模誤差的存在，測量邊界電壓信號的微小變化也會造成重構(gòu)阻抗分布的劇烈變化。高阻抗顱骨對激勵電流的阻擋作用加大了逆問題的不適定性，使得內(nèi)部的電阻抗變化更難體現(xiàn)在邊界電壓的變化上[7-9]。目前，廣泛使用的方法是通過各種正則化方法(如先驗正則化方法、圖像平滑等)降低重構(gòu)矩陣的條件數(shù)，盡可能求取接近精確解的近似解，從而得到穩(wěn)定而精確的圖像[10]。該類方法重建圖像的穩(wěn)定程度和精確程度依賴于正則化參數(shù)，因此選取適當?shù)恼齽t化參數(shù)則成為獲取最優(yōu)近似解的重要途徑。

在EIT領(lǐng)域正則化參數(shù)選取方法研究上，Graham等[11]用仿真的方法，在二維圓域且內(nèi)部阻抗呈均勻分布的情況下仿真比較了L型曲線法(L-curve)、廣義交叉校驗法(generalized cross-validation，GCV)、固定噪聲系數(shù)法和偏差原理法等4種方法，得出固定噪聲系數(shù)法(fixed noise figure，F(xiàn)NF)是較好的正則化參數(shù)選取方法；Abascal等[12]對真實邊界的三維模型進行了剖分和仿真，比較了L型曲線法等4種方法，得出L型曲線法和廣義交叉校驗法效果好，并將其運用到嬰兒腦部電阻抗成像中，得到了較好的成像結(jié)果。研究顯示，Graham等[11]的研究是針對圓域均勻模型，Abascal等[12]的研究是針對三維真實邊界阻抗均勻分布的模型，二者均未考慮高阻抗顱骨對正則化參數(shù)選取方法的影響，且得到的結(jié)論并不相同。因此，本研究擬通過含高阻抗顱骨的真實顱腦模型，比較L型曲線法、廣義交叉校驗法、固定噪聲系數(shù)法、偏差原理(discrepancy principle，DP)4種方法的效果，希望得出一種適用于腦部EIT的正則參數(shù)選取方法。

1 EIT圖像重建

EIT圖像重建是通過邊界電壓求解內(nèi)部阻抗分布或阻抗變化，即通過方程Ax=y求取其中的x。由于逆問題的不適定性，通常采用正則化技術(shù)求解。目前，廣泛采用的L2正則化技術(shù)，是通過最小化泛函(公式1)：

來求解反問題。式中，L為特定連續(xù)二乘函數(shù)形成的正則化矩陣，x0表示模型參數(shù)x的初值，α為正則化參數(shù)。其相應(yīng)的正則化解為(公式2)：

式中W為加權(quán)矩陣，R為正則化矩陣，該矩陣與相應(yīng)的先驗信息有關(guān)，一般可取單位對角矩陣、高通高斯濾波矩陣、拉普拉斯濾波矩陣、對角矩陣[R]對角矩陣等[13]。

2 正則參數(shù)對圖像重建的影響

區(qū)域內(nèi)重構(gòu)阻抗大于半高寬的單元可以量化反映重構(gòu)圖像的效果，因此，僅將重構(gòu)值大于半高寬的單元顯示出來(黑色)，當正則參數(shù)過小時，大于半高寬的單元分布比較分散；隨著正則參數(shù)的增大，大于半高寬的單元逐漸聚攏，當正則參數(shù)取得合適的值時，大于半高寬的單元面積取得最??；而當正則參數(shù)繼續(xù)增大，大于半高寬的單元不再聚攏而是擴散。因此，選取合適的正則參數(shù)才能重建出理想的圖像(如圖1所示)。

圖1 不同正則參數(shù)下的重構(gòu)圖像

3 四種正則化參數(shù)選取方法

3.1 L型曲線法

該方法是最為廣知的選取正則參數(shù)的方法，其以解范數(shù)的對數(shù)lg‖xα‖和殘余范數(shù)的對數(shù)lg‖Axα-y‖的關(guān)系曲線對比來確定正則參數(shù)。該曲線往往呈一個明顯的L型，故名L型曲線法。運用L型曲線準則的關(guān)鍵是求取曲線的隅角，該隅角反映了極小化殘差‖Axα-y‖與極小化解模‖xα‖之間的平衡。令ρ=lg‖Axα-y‖，θ=lg‖xα‖，則曲率關(guān)于正則參數(shù)α函數(shù)可定義為(公式3)：

其中“'”表示關(guān)于α的微分。曲率等于最大值對應(yīng)的α即為最優(yōu)正則化參數(shù)。

3.2 廣義交叉校驗法

廣義交叉驗證法的基本思想是，如果將任意一個觀測值Li從觀測序列L中刪除，則此時由剩余觀測值求得的正則化解應(yīng)能夠較好地預(yù)測L中被去掉的這一觀測值Li。廣義交叉驗證可以等效為求解最小GCV問題(公式4)：

式中，“Tr(·)”表示求矩陣的跡，B使得xα=By得到滿足。當GCV函數(shù)達到極小值時，對應(yīng)的α即為最優(yōu)正則化參數(shù)。

3.3 固定噪聲系數(shù)法

固定噪聲系數(shù)法是基于1996年Adler和Guardo提出的噪聲系數(shù)(Noise Figure，NF)[14]，該參數(shù)被定義為測量信噪比和圖像信噪比的比值(公式5)：

上式中信號定義為yc=Axc，xc為一個小的在區(qū)域中心的阻抗變化。一般情況下，取NF=0.5～2所對應(yīng)的正則化參數(shù)。

3.4 偏差原理

很多情況下，原始資料y中的誤差水平的參數(shù)δ往往是可以獲取或近似獲取的，在這種情況下，我們可以通過偏差方程(公式6)：

記p(α)=‖Axα-y‖2-δ2，對p(α)在α=αn作泰勒展開并在三階處截斷，則有(公式7)：

結(jié)合式(6)、(7)得出偏差原理的迭代公式，且該式是三階收斂的[15](公式8)：

設(shè)置合適的初始正則參數(shù)，運用(8)式迭代數(shù)次即得最優(yōu)化正則參數(shù)。

4 不同算法的量化比較

不同的方法是通過模糊半徑和重建位置誤差來進行比較的(公式9)：

式中Aq為重建阻抗值大于重建阻抗最大值一半的面積，A為整個腦部的面積，分辨率反映了圖像重建的精度。

位置誤差(position error，PE)反映了重構(gòu)目標與真實目標的位置偏移程度。

為使結(jié)果有統(tǒng)計學(xué)意義，對邊界電壓施加50次不同的0.1%的高斯白噪聲，對L型曲線法、GCV、NF法和偏差原理法4種方法進行比較，結(jié)果以(均值±均方差)的形式給出。

5 結(jié)果

5.1 數(shù)據(jù)來源

研究所采用的仿真模型是用有限元方法生成的，模型如圖2(a)所示。為了模擬高阻抗顱骨的真實顱腦模型，模型中的阻抗分布是非均勻的，電阻率值設(shè)置為∶皮膚∶顱骨∶腦實質(zhì)=2.27 Ωm∶25.64 Ωm∶3.33 Ωm[4]，在顱腦的邊界等間距均勻設(shè)置16個電極，采用對向驅(qū)動模式，激勵電流為1mA。根據(jù)重建需要，在真實顱腦模型中2/3半徑處設(shè)置了低于腦實質(zhì)阻抗的電阻抗變化，電阻率大小設(shè)置為1.48 Ωm(模擬顱內(nèi)出血)[4]，如圖2(b)所示。由于正逆問題采用同樣的模型會產(chǎn)生逆問題補償，因此本研究采用圖2(c)中的逆問題模型重建圖像。

圖2 全腦剖分模型

5.2 一例重建結(jié)果

為進一步闡述不同方法的結(jié)果，取出一例重建結(jié)果進行展示。如圖3所示，展示了正則參數(shù)連續(xù)變化的影響。對于小的α，殘差‖Axα-y‖較小(圖3a)但是解?！瑇α‖較大(圖3b)；與之相應(yīng)，對于大的α，殘差‖Axα-y‖較大(圖3-a)但是解?！瑇α‖較小(圖3-b)；圖3-c為L型曲線的曲率與正則參數(shù)的關(guān)系，圖中曲率最大值對應(yīng)的正則參數(shù)即為最優(yōu)正則參數(shù)；圖3-d為GCV函數(shù)與正則參數(shù)的關(guān)系，圖中函數(shù)極小值對應(yīng)的正則參數(shù)為最優(yōu)正則參數(shù)。

4種方法求得的最優(yōu)正則參數(shù)在L型曲線上的位置如圖4所示。

圖3 向邊界電壓添加0.1%高斯白噪聲的一例重建結(jié)果

圖4 4種方法選取的正則參數(shù)在L型曲線上的位置

4種方法選取的正則參數(shù)下重建得到的圖像如圖5所示。

圖5 4種方法重構(gòu)內(nèi)部阻抗的歸一化顯示

5.3 不同方法的統(tǒng)計結(jié)果

本研究添加50次不同的0.1%高斯白噪聲的不同方法的統(tǒng)計結(jié)果見表1。

表1 添加50次不同的0.1%高斯白噪聲的不同方法的統(tǒng)計結(jié)果

6 討論

腦部電阻抗圖像重建曾經(jīng)使用的是固定正則參數(shù)的方法，該方法已不適應(yīng)實時最優(yōu)化圖像的需要。圖1結(jié)合圖3說明，在正則參數(shù)α很小時，‖Axα-y‖很小，說明xα逼近了真實解，但是‖xα‖很大，說明解是不穩(wěn)定的；當正則參數(shù)α較大時，‖xα‖很小，表明解是穩(wěn)定的，但是‖Axα-y‖很大時則表明xα與真實解相去甚遠。最優(yōu)的正則化參數(shù)必須同時兼顧解的“精確性”和“穩(wěn)定性”。為此，本研究在含高阻抗顱骨的顱腦模型中，量化比較了4種正則化參數(shù)選取方法：L型曲線法、廣義交叉校驗法、固定噪聲系數(shù)法和偏差原理法，希望找到一種適用于腦部EIT的選取正則參數(shù)的算法。

(1)L型曲線法：統(tǒng)計結(jié)果表明，雖然L型曲線法取得了最小的模糊半徑，但L型曲線法重建所得的目標位置誤差比較大，達到0.0354±0.0154，表明L型曲線法重建目標位置不準確；此外，當邊界電壓噪聲較大時，曲線并不呈現(xiàn)明顯的L型，而是變得扁平，在這種情況下曲線法求出的拐點并不準確甚至無法求得拐點，即無法求得最優(yōu)的正則參數(shù)。

(2)廣義交叉校驗法：廣義交叉校驗法重建所得的目標位置誤差為0.0187±0.0175，位置誤差的均方差較大，表明通過該方法重建目標的位置不穩(wěn)定。Graham等[11]的研究亦表明，GCV方法在電阻抗斷層成像中不可靠。

(3)固定噪聲系數(shù)法：本研究中在使用固定噪聲系數(shù)法取NF=0.5時取得了較好的成像結(jié)果。該方法提供了一種新的選取正則參數(shù)的方法，可重復(fù)性好，并且比L型曲線法、GCV方法更穩(wěn)定，但是不同的NF取值會對重建圖像結(jié)果有較大影響，而NF取值存在很大的主觀因素。因此，固定噪聲系數(shù)法不夠客觀。

(4)偏差原理法：偏差原理法結(jié)果較好，模糊半徑處于合理水平，且多次成像位置誤差小，位置誤差的方差較小，成像結(jié)果穩(wěn)定，是較好的正則參數(shù)選取方法。

研究顯示，偏差原理法成像效果比較依賴對于誤差的估計，目前一般采用先采集一段數(shù)據(jù)，然后根據(jù)預(yù)采集數(shù)據(jù)來估計誤差，但這種方法有一定的局限性；此外，重建模型的剖分規(guī)模的增大，會導(dǎo)致偏差原理法耗時指數(shù)式增長。因此，提高誤差估計精確性和利用更快的迭代算法加快偏差原理法的運行速度有待進一步研究。

7 結(jié)論

本研究通過統(tǒng)計比較4種選取正則參數(shù)的方法得出結(jié)論：偏差原理法取得了最好的結(jié)果，重建圖像的結(jié)果最客觀、最穩(wěn)定和最接近真實解，可采用偏差原理法作為腦部EIT正則參數(shù)的選取算法。

[1]董秀珍.生物電阻抗成像研究的現(xiàn)狀與挑戰(zhàn)[J].中國生物醫(yī)學(xué)工程年報,2008,27(5):641-643.

[2]史學(xué)濤,霍旭陽,尤富生,等.顱內(nèi)出血電阻抗成像系統(tǒng)及初步動物實驗[J].航天醫(yī)學(xué)與醫(yī)學(xué)工程,2007,20(1):24-27.

[3]劉麗旭,董為偉,賈建平,等.無創(chuàng)性腦電阻抗測定在腦出血患者腦水腫監(jiān)測中的應(yīng)用[J].中華神經(jīng)科雜志,2007,40(6):383-386.

[4]Dai M,Wang L,Xu C,et al.Real-time imaging of subarachnoid hemorrhage in piglets with electrical impedance tomography[J].Physiol Meas,2010,31(9):1229-1239.

[5]Xu C,Dai M,You F,et al.An optimized strategy for real-time hemorrhage monitoring with electrical impedance tomography[J].Physiol Meas,2011,32(5):585-598.

[6]Dai M,Li B,Hu S,et al.In vivo imaging of twist drill drainage for subdural hematoma:a clinical feasibility study on electrical impedance tomography for measuring intracranial bleeding in humans[J].PloS One,2013,8(1):e55020.

[7]倪安勝,楊國勝,付峰,等.基于非均勻顱骨頭模型的電阻抗斷層成像圖像重構(gòu)[J].航天醫(yī)學(xué)與醫(yī)學(xué)工程,2008,21(1):56-60.

[8]Sadleir RJ,Argibay A.Modeling skull electrical properties[J].Ann Biomed Eng,2007,35(10):1699-1712.

[9]Tang C,You F,Cheng G,et al.Modeling the frequency dependence of the electrical properties of the live human skull[J].Physiol Meas,2009,30(12):1293-1301.

[10]Vauhkonen M,Vadász D,Karjalainen PA,et al.Tikhonov regularization and prior information in electrical impedance tomography[J].IEEE Trans Med Imaging,1998,17(2):285-293.

[11]Graham BM,Adler A.Objective selection of hyperparameter for EIT[J].Physiol Meas,2006,27(5):S65-S79.

[12]Abascal JF,Arridge SR,Bayford RH,et al.Comparison of methods for optimal choice of the regularization parameter for linear electrical impedance tomography of brain function[J].Physiol Meas,2008,29(11):1319-1334.

[13]Adler A,Dai T,Lionheart WR.Temporal image reconstruction in electrical impedance tomography[J].Physiol Meas,2007,28(7):S1-S11.

[14]Adler A,Guardo R.Electrical impedance tomography:regularized imaging and contrast detection[J].IEEE Trans Med Imaging,1996,15(2):170-179.

[15]Wang Y,Xiao T.Fast realization algorithms for determining regularization parameters in linear inverse problems[J].Inverse problems,2001,17(2):281-291.