針對(duì)腦部電阻抗成像的四種正則參數(shù)選取方法的比較研究*

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

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

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

在EIT領(lǐng)域正則化參數(shù)選取方法研究上，Graham等[11]用仿真的方法，在二維圓域且內(nèi)部阻抗呈均勻分布的情況下仿真比較了L型曲線法(L-curve)、廣義交叉校驗(yàn)法(generalized cross-validation，GCV)、固定噪聲系數(shù)法和偏差原理法等4種方法，得出固定噪聲系數(shù)法(fixed noise figure，F(xiàn)NF)是較好的正則化參數(shù)選取方法；Abascal等[12]對(duì)真實(shí)邊界的三維模型進(jìn)行了剖分和仿真，比較了L型曲線法等4種方法，得出L型曲線法和廣義交叉校驗(yàn)法效果好，并將其運(yùn)用到嬰兒腦部電阻抗成像中，得到了較好的成像結(jié)果。研究顯示，Graham等[11]的研究是針對(duì)圓域均勻模型，Abascal等[12]的研究是針對(duì)三維真實(shí)邊界阻抗均勻分布的模型，二者均未考慮高阻抗顱骨對(duì)正則化參數(shù)選取方法的影響，且得到的結(jié)論并不相同。因此，本研究擬通過含高阻抗顱骨的真實(shí)顱腦模型，比較L型曲線法、廣義交叉校驗(yàn)法、固定噪聲系數(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)的先驗(yàn)信息有關(guān)，一般可取單位對(duì)角矩陣、高通高斯濾波矩陣、拉普拉斯濾波矩陣、對(duì)角矩陣[R]對(duì)角矩陣等[13]。

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

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

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

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

3.1 L型曲線法

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

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

3.2 廣義交叉校驗(yàn)法

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

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

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

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

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

3.4 偏差原理

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

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

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

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

4 不同算法的量化比較

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

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

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

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

5 結(jié)果

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

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

圖2 全腦剖分模型

5.2 一例重建結(jié)果

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

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

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

6 討論

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

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

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

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

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

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

7 結(jié)論

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

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

[2]史學(xué)濤,霍旭陽,尤富生,等.顱內(nèi)出血電阻抗成像系統(tǒng)及初步動(dòng)物實(shí)驗(yàn)[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.