李穎 ，馬重蕾 ，趙營(yíng)鴿 ，王冠雄 ，郝虎鵬

（1.河北工業(yè)大學(xué) 生命科學(xué)與健康工程學(xué)院，天津 300130；2.河北工業(yè)大學(xué) 河北省生物電磁與神經(jīng)工程重點(diǎn)實(shí)驗(yàn)室，天津 300130；3.河北工業(yè)大學(xué) 天津市生物電工與智能健康重點(diǎn)實(shí)驗(yàn)室，天津 300130；4.新鄉(xiāng)醫(yī)學(xué)院三全學(xué)院 智能醫(yī)學(xué)工程學(xué)院，河南 新鄉(xiāng) 453003）

電阻抗成像（EIT）的基本原理是在人體表面施加安全的激勵(lì)電流，測(cè)量由此電流產(chǎn)生的目標(biāo)表面的電壓信號(hào)，重構(gòu)出被測(cè)目標(biāo)內(nèi)部的阻抗圖像分布［1］.由于它具有操作簡(jiǎn)單、成本低廉、無(wú)損傷等特點(diǎn)，在醫(yī)學(xué)成像方面獲得了廣泛的關(guān)注.在對(duì)EIT 的研究中，為了簡(jiǎn)化計(jì)算，通常假設(shè)生物組織的阻抗分布是不變的，但實(shí)際情況中，阻抗分布是變化的，這就導(dǎo)致圖像重建結(jié)果是不確定的［2］.因此研究電阻率的變化對(duì)EIT逆問(wèn)題的影響具有重要的意義.

目前，EIT逆問(wèn)題求解的很多算法都是在確定性框架下的，如牛頓迭代類(lèi)方法、神經(jīng)網(wǎng)絡(luò)類(lèi)算法、進(jìn)化類(lèi)全局優(yōu)化算法等，以及各類(lèi)算法的改進(jìn)形式，這些算法可對(duì)目標(biāo)內(nèi)的阻抗分布進(jìn)行重構(gòu)，并試圖提高重構(gòu)圖像的分辨率［3-5］，但它們對(duì)于電阻率分布的不確定性欠考慮.為了使逆問(wèn)題計(jì)算結(jié)果更加全面、精確，考慮不確定性因素是必要的.對(duì)逆問(wèn)題進(jìn)行不確定性量化可以歸納為統(tǒng)計(jì)推斷問(wèn)題，貝葉斯方法是概率統(tǒng)計(jì)范疇的方法，它不僅可以對(duì)參數(shù)進(jìn)行反演，還能對(duì)輸出結(jié)果進(jìn)行不確定性量化，因此在逆問(wèn)題求解中的應(yīng)用很廣泛［6-10］.

貝葉斯方法是一種結(jié)合先驗(yàn)信息和似然函數(shù)來(lái)獲得后驗(yàn)分布的方法.在一些非線性的復(fù)雜模型中，計(jì)算后驗(yàn)分布時(shí)往往難以形成解析式，通常采用有效的抽樣算法來(lái)獲得穩(wěn)定的后驗(yàn)分布.馬爾科夫鏈蒙特卡洛（Markov Chain Monte Carlo，MCMC）算法憑借其良好的收斂?jī)?yōu)勢(shì)成為目前應(yīng)用最為廣泛的抽樣算法［11］.MCMC算法是一種通過(guò)構(gòu)建馬爾科夫鏈來(lái)獲得隨機(jī)參數(shù)樣本的間接抽樣方法，在復(fù)雜的非線性模型中，需要大量反復(fù)地進(jìn)行抽樣計(jì)算，雖然可以取得較好的效果，但其計(jì)算效率很低，尤其是面對(duì)高維度參數(shù)時(shí)［12］.為了解決這個(gè)問(wèn)題，可以在MCMC 算法中引入自適應(yīng)的思想，例如自適應(yīng)Metropolis 算法和延遲拒絕自適應(yīng)Metropolis 算法等［13-14］，但這些算法都是單鏈運(yùn)行.有研究者提出一種差分進(jìn)化自適應(yīng)Metropolis 算法，被命名為DREAM 算法，它是一種利用多條隱馬爾科夫鏈并列運(yùn)行的改進(jìn)MCMC 算法，利用差分進(jìn)化使其收斂和計(jì)算速度得到了很大的提升，而進(jìn)一步發(fā)展的自適應(yīng)差分進(jìn)化Metropolis 抽樣（Differential Evolution Adaptive Metropolis，DREAM_ZS）算法在DREAM 算法基礎(chǔ)上進(jìn)行了改進(jìn)，在參數(shù)反演及其不確定性量化方面得到了廣泛的應(yīng)用［15-16］.

另外，在貝葉斯方法的實(shí)現(xiàn)中，對(duì)后驗(yàn)分布的計(jì)算需要調(diào)用大量與正問(wèn)題模型相關(guān)的似然函數(shù)計(jì)算，正問(wèn)題原始模型計(jì)算耗時(shí)過(guò)長(zhǎng)會(huì)大大影響整體效率，所以采用計(jì)算效率高的模型來(lái)替代原始模型是很有必要的.常用的替代模型有神經(jīng)網(wǎng)絡(luò)類(lèi)模型［17-18］、高斯過(guò)程回歸模型［19］、克里金（Kriging）模型［20］等.反向傳播（Back Propagation，BP）神經(jīng)網(wǎng)絡(luò)由于其非線性映射能力強(qiáng)、計(jì)算效率高，被廣泛應(yīng)用于復(fù)雜非線性模型的替代模型.

本文針對(duì)EIT 逆問(wèn)題中電阻率分布不確定的情況，采用BP 神經(jīng)網(wǎng)絡(luò)模型替代求解EIT 正問(wèn)題的有限元模型，基于貝葉斯理論，采用DREAM_ZS算法求得電阻率的后驗(yàn)概率密度，從而對(duì)EIT 逆問(wèn)題的不確定性進(jìn)行量化.

1 理論與方法

1.1 EIT逆問(wèn)題的求解模型

EIT研究的是一個(gè)特殊的電流場(chǎng)問(wèn)題，可以將待測(cè)目標(biāo)看作成一個(gè)導(dǎo)體，電流的流動(dòng)受場(chǎng)域內(nèi)電導(dǎo)率分布的影響［21］.EIT 數(shù)學(xué)物理模型可以用一個(gè)二階橢圓型偏微分方程來(lái)表示［22］.設(shè)物體所在區(qū)域?yàn)棣?，邊界為Γ，假設(shè)物體內(nèi)部無(wú)源，僅考慮電導(dǎo)率σ分布時(shí)，有

EIT 逆問(wèn)題是根據(jù)測(cè)得的邊界電壓分布反演內(nèi)部阻抗分布.由于電流密度是阻抗分布的函數(shù)，所以EIT 逆問(wèn)題是個(gè)非線性問(wèn)題.在實(shí)際測(cè)量中，電極數(shù)量和電流注入方式是有限的，這就導(dǎo)致了從邊界電極獲得的信息不足以確定目標(biāo)內(nèi)部的阻抗分布，而且，電極位置對(duì)阻抗分布的敏感性不同.這些使得EIT逆問(wèn)題成為一個(gè)高度病態(tài)的非線性問(wèn)題.

由于EIT 逆問(wèn)題的不適定性和病態(tài)性十分嚴(yán)重，在求解時(shí)需要引入先驗(yàn)信息.結(jié)合了先驗(yàn)信息和觀測(cè)數(shù)據(jù)的貝葉斯方法可以量化逆問(wèn)題解的不確定性，有些情況下可以降低不適定問(wèn)題的不適定程度［10］.

1.2 貝葉斯反演的基本原理及方法

貝葉斯理論描述的是由果及因，即利用事情發(fā)生后的結(jié)果推出導(dǎo)致事情發(fā)生的原因.它的數(shù)學(xué)表達(dá)式為

式中，p（m|d）為后驗(yàn)分布，p（d|m）為似然函數(shù)，p（m）為先驗(yàn)分布，p（d）為測(cè)量數(shù)據(jù)的概率密度.在參數(shù)反演問(wèn)題中m為待反演參數(shù)，d為測(cè)量數(shù)據(jù)，由于p（d）與反演參數(shù)無(wú)關(guān)，所以可看成一個(gè)常量.那么上式也可以表示為

可見(jiàn)參數(shù)的后驗(yàn)分布完全是由先驗(yàn)分布以及似然函數(shù)決定的.先驗(yàn)參數(shù)分布是由人的日常生活經(jīng)驗(yàn)或者科學(xué)家及權(quán)威科研資料所給出的.似然函數(shù)表征的是理論數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)之間的差異程度.最常用的似然函數(shù)是高斯似然函數(shù)，即假設(shè)誤差向量δi服從獨(dú)立標(biāo)準(zhǔn)正態(tài)分布，其似然函數(shù)計(jì)算表達(dá)式為

式中，L（x）為似然函數(shù)，xi為反演參數(shù)，（fxi）為正問(wèn)題計(jì)算結(jié)果，即模擬計(jì)算數(shù)據(jù)，udi為實(shí)測(cè)數(shù)據(jù)，nd為參數(shù)數(shù)目.

MCMC 算法通常用于非線性模型的后驗(yàn)分布計(jì)算.MCMC 算法通過(guò)構(gòu)造一條具有轉(zhuǎn)移特性的隱馬爾科夫鏈來(lái)探尋未知參數(shù)的后驗(yàn)分布空間。假設(shè)該鏈為x={xt；t＞0}，在抽樣計(jì)算過(guò)程中，該鏈在t+1時(shí)刻的概率值p(xt+1)僅僅只與t時(shí)刻的概率值p(xt)有關(guān).此轉(zhuǎn)移特性表示為

當(dāng)該馬爾科夫鏈的轉(zhuǎn)移概率滿足此特性時(shí)，說(shuō)明該鏈已經(jīng)達(dá)到可靠的收斂程度，此時(shí)該鏈的分布即為參數(shù)的穩(wěn)定后驗(yàn)分布.

但MCMC 算法在面對(duì)高維參數(shù)模型時(shí)，其收斂性和收斂速度都受到了很大的限制.DREAM_ZS 是一種多鏈MCMC 算法，它在DREAM 算法的基礎(chǔ)上增加了簡(jiǎn)化離散和組合搜索參數(shù)空間的擴(kuò)展功能，可以通過(guò)差分進(jìn)化、子空間探索、從過(guò)去狀態(tài)采樣、斯諾克更新和多重嘗試Metropolis 采樣的組合來(lái)有效地探索高維和多模態(tài)后驗(yàn)概率密度，具有較強(qiáng)的收斂能力和魯棒性［23］.DREAM_ZS算法基本步驟為：

1）根據(jù)n個(gè)待反演參數(shù)的先驗(yàn)分布抽樣出N個(gè)初始樣本，作為初始種群（k=1），代入正演模型，利用式（3）、式（4）計(jì)算樣本的后驗(yàn)概率分布.

2）對(duì)第k代每一個(gè)個(gè)體進(jìn)行變異操作.

式 中，vk為變異 后個(gè)體，γ為縮放因子，一 般γ=為模型參數(shù)，ε為方差.

3）按照交叉概率CR進(jìn)行交叉操作，之后判斷是否取代新樣本.若U≤1 -CR，則令未變異樣本點(diǎn)取代變異后的樣本點(diǎn)，反之不取代.

式中，U是一個(gè)取值在0～1 的隨機(jī)數(shù)，交叉概率CR∈［0，1］，根據(jù)抽樣自適應(yīng)性發(fā)生變化.

4）對(duì)新樣本計(jì)算接受概率.若接受新樣本，則取變異后的樣本點(diǎn)；若不接受，保持原位置狀態(tài)，使平行序列繼續(xù)進(jìn)化.

式中：q為用于評(píng)價(jià)的馬爾科夫鏈的條數(shù)；B為鏈內(nèi)方差；為q條馬爾科夫鏈均值的方差；為uj的均值；W為q條馬爾科夫鏈方差Sj的均值.若值達(dá)標(biāo)，表明馬爾科夫鏈已可靠收斂；若不達(dá)標(biāo)，則繼續(xù)進(jìn)行序列進(jìn)化.

1.3 基于貝葉斯框架的EIT參數(shù)反演

對(duì)于EIT 參數(shù)反演，電阻率ρ和測(cè)量的電壓U被假設(shè)為具有某種聯(lián)合概率的多元隨機(jī)變量.貝葉斯方法的目標(biāo)是在給定電位測(cè)量值的情況下求解電阻率分布的后驗(yàn)概率.

使用貝葉斯定理，電阻率后驗(yàn)概率表示為p(ρ|U)，表征測(cè)量電位數(shù)據(jù)的似然函數(shù)表示為p(U|ρ)，p(ρ)表示為電阻率的先驗(yàn)分布，U(ρ)為EIT正問(wèn)題計(jì)算得到的模擬電位值，udi為理論電位值.引入高斯誤差項(xiàng)ε～N（μδi，σδi2）來(lái)表征模擬電位值和實(shí)測(cè)電位值之間的關(guān)系.

綜上所述，基于貝葉斯框架的EIT 參數(shù)反演有以下步驟：

1）根據(jù)電阻率ρ的先驗(yàn)分布產(chǎn)生N組初始樣本；

2）將樣本值輸入EIT 正問(wèn)題計(jì)算模型，計(jì)算出理論電位值U（ρ）；

3）引入高斯誤差項(xiàng)，根據(jù)式（13）計(jì)算似然函數(shù)p(U|ρ)；

4）求解式（14）得出電阻率的后驗(yàn)分布函數(shù)p(ρ|U)；

5）采用DREAM_ZS抽樣算法對(duì)后驗(yàn)分布進(jìn)行抽樣，得出各個(gè)電阻率的邊緣后驗(yàn)密度；

6）對(duì)電阻率的反演結(jié)果進(jìn)行不確定性分析.

1.4 EIT正問(wèn)題計(jì)算的BP神經(jīng)網(wǎng)絡(luò)替代模型

在EIT 正問(wèn)題的研究中，往往采用有限元法進(jìn)行求解.如果在EIT 模型中反復(fù)調(diào)用有限元來(lái)進(jìn)行計(jì)算，將會(huì)導(dǎo)致計(jì)算效率很低.BP 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)簡(jiǎn)單、非線性映射能力強(qiáng)，可以明顯提高計(jì)算速度，故本文采用BP神經(jīng)網(wǎng)絡(luò)來(lái)替代有限元模型.

BP 神經(jīng)網(wǎng)絡(luò)主要由三部分構(gòu)成：輸入層、隱含層和輸出層［25］.在神經(jīng)網(wǎng)絡(luò)的訓(xùn)練過(guò)程中進(jìn)行了兩種計(jì)算，一種是前向計(jì)算，另一種是逆向反饋計(jì)算.

具體地，前向計(jì)算過(guò)程為

式中，ωij表示節(jié)點(diǎn)i與節(jié)點(diǎn)j之間的權(quán)值，bj為節(jié)點(diǎn)j的閾值，xi為節(jié)點(diǎn)輸入值，f為激活函數(shù)，yj為節(jié)點(diǎn)輸出值.

逆向反饋過(guò)程為

式中，dj為假定的輸出結(jié)果，e(ω，b)表示誤差函數(shù)，表示修正后的權(quán)值，表示修正后的閾值.當(dāng)達(dá)到目標(biāo)誤差期望函數(shù)時(shí)，BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練完成.

2 仿真模型的構(gòu)建與正問(wèn)題計(jì)算

在EIT 正問(wèn)題的計(jì)算過(guò)程中，采用有限元網(wǎng)絡(luò)對(duì)模型進(jìn)行剖分，仿真模型如圖1 和圖2 所示.圖1為四層同心圓模型，由內(nèi)到外分別代表大腦、腦脊液、顱骨和頭皮，每層的電阻率一致，其參考值如表1所示［26］.圖2 為高維參數(shù)模型，它是將圓域剖分成64個(gè)小三角單元，每一個(gè)單元的電阻率都不同，總計(jì)64個(gè)電阻率值.兩個(gè)仿真模型均采用16 個(gè)電極，均勻放置于目標(biāo)表面，目標(biāo)外界半徑為10 cm.在仿真實(shí)驗(yàn)中，設(shè)置激勵(lì)電流為1 mA，50 kHz.

表1 四層同心圓模型參數(shù)取值表Tab.1 Parameter values of the four-layer concentric circle model

圖1 四層同心圓模型Fig.1 Four-layer concentric circle model.

圖2 高維參數(shù)模型Fig.2 High dimensional parameter model.

為了提高有限元模型的計(jì)算效率，我們采用BP神經(jīng)網(wǎng)絡(luò)模型作為有限元的替代模型.首先將有限元模型的電阻率值作為網(wǎng)絡(luò)的輸入，計(jì)算出的電極電位結(jié)果作為網(wǎng)絡(luò)的輸出，隨機(jī)生成1 000組輸入輸出數(shù)據(jù)，選擇數(shù)據(jù)的70%作為網(wǎng)絡(luò)的訓(xùn)練數(shù)據(jù)、15%作為修正數(shù)據(jù)、15%作為測(cè)試數(shù)據(jù)，訓(xùn)練BP 人工神經(jīng)網(wǎng)絡(luò).

利用均方差來(lái)度量BP 神經(jīng)網(wǎng)絡(luò)的性能，均方誤差MSE的計(jì)算公式為

式中，N為樣本個(gè)數(shù)，Yi為驗(yàn)證集輸出結(jié)果，yi為測(cè)試集輸出結(jié)果.圖3 顯示了BP 神經(jīng)網(wǎng)絡(luò)訓(xùn)練過(guò)程中的均方差變化.

圖3 BP神經(jīng)網(wǎng)絡(luò)均方差Fig.3 BP neural network variance diagram.

由圖3 可以看出各個(gè)測(cè)試集最終都趨于收斂，且在迭代完成時(shí)四層同心圓模型的均方差低于10-8，高維參數(shù)模型的均方差低于10-5，可見(jiàn)BP 神經(jīng)網(wǎng)絡(luò)擬合得很好.

進(jìn)一步對(duì)訓(xùn)練好的網(wǎng)絡(luò)進(jìn)行驗(yàn)證，按均值±30%隨機(jī)生成20 組四層電阻率值，分別代入有限元模型以及BP 神經(jīng)網(wǎng)絡(luò)模型.計(jì)算16 個(gè)電極20 組數(shù)據(jù)的平均相對(duì)誤差值，如圖4 所示，由于9 號(hào)電極為參考電極，圖中忽略不計(jì).

圖4 平均相對(duì)誤差Fig.4 Average relative error.

由圖4 可以看出，在四層同心圓模型中BP 神經(jīng)網(wǎng)絡(luò)模型與有限元模型的誤差在10-5左右，二維圓模型中誤差為10-3，說(shuō)明兩種計(jì)算模型擬合得很好.在程序運(yùn)行中，生成1 000 組輸入數(shù)據(jù)，四層同心圓模型用時(shí)116.446 s，BP 神經(jīng)網(wǎng)絡(luò)用時(shí)3.281 s；高維參數(shù)模型用時(shí)56 875.324 s，BP 神經(jīng)網(wǎng)絡(luò)用時(shí)389.732 s，說(shuō)明使用BP 神經(jīng)網(wǎng)絡(luò)替代模型極大地提升了計(jì)算速度.故BP 神經(jīng)網(wǎng)絡(luò)可以作為有限元的替代模型進(jìn)行后續(xù)的參數(shù)反演工作.

3 參數(shù)反演結(jié)果

3.1 基于同心圓模型的低維參數(shù)反演

針對(duì)同心圓模型進(jìn)行參數(shù)反演，取四層電阻率的先驗(yàn)分布為均勻分布，變化范圍取均值的±30%，如表1所示，并假設(shè)噪聲服從獨(dú)立的正態(tài)分布N（0，22）.分別選取相鄰激勵(lì)、相對(duì)激勵(lì)以及間隔六電極激勵(lì)三種電流激勵(lì)方式，來(lái)探討不同的激勵(lì)方式對(duì)電阻率反演結(jié)果的影響.其中，在相鄰激勵(lì)方式下電流從1號(hào)電極注入，2 號(hào)電極流出；在相對(duì)激勵(lì)方式下電流從1 號(hào)電極注入，9 號(hào)電極流出；在間隔六電極方式下電流從1 號(hào)電極注入，8 號(hào)電極流出.利用DREAM_ZS 抽樣算法進(jìn)行抽樣，其中DREAM_ZS 算法的參數(shù)設(shè)置為：馬爾科夫鏈并行運(yùn)行的數(shù)量為3條，生成候選點(diǎn)的鏈對(duì)數(shù)量為1 對(duì)，使用3%的交叉值.三種激勵(lì)方式下電阻率的后驗(yàn)均值以及與真值之間的誤差如表2所示.

表2 各激勵(lì)方式下參數(shù)后驗(yàn)分布統(tǒng)計(jì)信息表Tab.2 Statistical information table of posterior distribution of parameters under various excitation modes

由表2 可知，每層電阻率的后驗(yàn)均值與真值之間的誤差很小，說(shuō)明DREAM_ZS 算法可以識(shí)別各個(gè)參數(shù).對(duì)比不同的激勵(lì)方式，當(dāng)電流激勵(lì)方式為相對(duì)激勵(lì)時(shí)誤差最小，其次是相鄰激勵(lì)，而間隔六電極激勵(lì)方式的誤差較其它兩種方式大很多，這說(shuō)明在相對(duì)激勵(lì)方式下參數(shù)反演效果是最優(yōu)的，而相隔六電極激勵(lì)不太理想.在三種激勵(lì)方式下，頭皮的誤差均為最小，說(shuō)明頭皮的不確定性最小.在相對(duì)激勵(lì)和相鄰激勵(lì)方式下，顱骨和大腦的誤差次之，誤差最大的是腦脊液；在間隔六電極激勵(lì)方式下，顱骨的誤差次之，大腦和腦脊液的誤差較大.

為了進(jìn)一步對(duì)后驗(yàn)分布的概率值進(jìn)行衡量，我們引入最大后驗(yàn)概率（Maximum a Posteriori，MAP）值.MAP 是理論上后驗(yàn)分布概率取最大時(shí)的參數(shù)值，可以判定算法對(duì)參數(shù)的識(shí)別程度.參數(shù)的后驗(yàn)分布峰值越接近于MAP，說(shuō)明對(duì)算法的識(shí)別效果越好.

三種電流激勵(lì)方式下的電阻率參數(shù)后驗(yàn)分布的邊緣概率密度曲線如圖5 所示，圖中藍(lán)色叉點(diǎn)為MAP值.

圖5 各激勵(lì)方式下參數(shù)后驗(yàn)分布邊緣概率密度曲線Fig.5 Marginal density curve of posterior distribution of parameters under various excitation modes

通過(guò)圖5 可以直觀地看出：間隔六電極激勵(lì)方式下，四層電阻率參數(shù)后驗(yàn)分布均未出現(xiàn)單個(gè)峰值的正態(tài)分布形狀；相鄰激勵(lì)方式下，頭皮和顱骨有近似的單個(gè)峰值的正態(tài)分布形狀；相對(duì)激勵(lì)方式下，頭皮、顱骨和腦脊液均有近似的單個(gè)峰值的正態(tài)分布形狀。在反演效果最佳的相對(duì)激勵(lì)方式下，頭皮和顱骨的峰值明顯，且接近MAP，腦脊液次之，而大腦的整體趨勢(shì)較為平緩，且峰值離MAP 較遠(yuǎn)，說(shuō)明DREAM_ZS 算法對(duì)頭皮、顱骨的識(shí)別能力強(qiáng)，其不確定性小，腦脊液和大腦的不確定性要大.

為了定量描述參數(shù)的后驗(yàn)分布信息，計(jì)算各激勵(lì)方式下標(biāo)準(zhǔn)差.其中標(biāo)準(zhǔn)差計(jì)算為

圖6 不同激勵(lì)方式下參數(shù)的標(biāo)準(zhǔn)差Fig.6 Standard deviation of parameters under different excitation modes

從圖6 可以看出，對(duì)于三種激勵(lì)模式，頭皮的標(biāo)準(zhǔn)差均最小，其次是顱骨，而大腦和腦脊液的標(biāo)準(zhǔn)差大一些.說(shuō)明四個(gè)參數(shù)中頭皮的后驗(yàn)分布最集中，不確定性最小，敏感性最強(qiáng).其次是顱骨，大腦和腦脊液的不確定性較大.

3.2 基于二維圓模型的高維參數(shù)反演

通過(guò)前面的結(jié)論可以得知在相對(duì)激勵(lì)方式下對(duì)逆問(wèn)題電阻率反演效果最好，所以在高維參數(shù)反演中電流采用相對(duì)激勵(lì)方式注入.根據(jù)貝葉斯理論可知，先驗(yàn)分布的不同在一定程度上會(huì)影響后驗(yàn)結(jié)果的分布，所以接下來(lái)研究電阻率在不同先驗(yàn)分布下（即均勻分布與正態(tài)分布下）參數(shù)的反演結(jié)果.

3.2.1 先驗(yàn)為均勻分布時(shí)參數(shù)的反演結(jié)果

在對(duì)高維參數(shù)模型進(jìn)行的仿真實(shí)驗(yàn)中，假設(shè)非病灶單元均值為1 Ω·m，而病灶單元電阻率為非病灶單元的2 倍.假設(shè)病灶單元和非病灶單元均服從變化范圍為均值±20%的均勻分布，即非病灶單元參數(shù)服從（0.8，1.2）上的均勻分布，病灶單元參數(shù)均服從（1.6，2.4）上的均勻分布.如圖7（a）所示，選取9號(hào)、10 號(hào)、25 號(hào)參數(shù)區(qū)域作為病灶單元進(jìn)行參數(shù)反演.采用DREAM_ZS算法進(jìn)行抽樣，可以得到所有單元的反演結(jié)果。各單元先驗(yàn)分布和后驗(yàn)分布的均值結(jié)果如圖7（b）、圖7（c）所示.

圖7 先驗(yàn)為均勻分布時(shí)的反演結(jié)果示意圖Fig.7 Schematic diagram of inversion results when prior distribution is uniform distribution

由圖7 可知，各單元先驗(yàn)分布和后驗(yàn)分布的均值相差不大，經(jīng)算法識(shí)別后的單元均值接近于真值，說(shuō)明DREAM_ZS算法可以識(shí)別出各個(gè)單元參數(shù).

圖8 顯示了單元9、單元10、單元25 的后驗(yàn)分布邊緣密度曲線.對(duì)于非病灶單元，由于單元數(shù)目過(guò)多，我們選取幾個(gè)典型單元的后驗(yàn)分布邊緣密度曲線展示于圖9中.

圖8 病灶單元后驗(yàn)分布邊緣密度曲線Fig.8 Marginal density curve of posterior distribution of parameters when prior is uniform distribution.

圖9 典型單元后驗(yàn)分布邊緣密度曲線Fig.9 Marginal density curve of posterior distribution of typical unit

由圖8 的曲線趨勢(shì)可以看出，單元10 的后驗(yàn)分布峰值最接近于MAP，其次是單元9，單元25的峰值與MAP相差最大.說(shuō)明DREAM_ZS算法對(duì)單元10的識(shí)別效果更強(qiáng)，對(duì)單元9 也可有效識(shí)別對(duì)單元25 的識(shí)別效果較差.單元10 的不確定性要小于單元9 和單元25，敏感性最強(qiáng).

在圖9 中，單元49 和單元50 的峰值最接近于MAP，其次是單元26，單元2 的峰值與MAP 相差最大.說(shuō)明DREAM_ZS算法對(duì)單元49和單元50的識(shí)別效果最強(qiáng)，單元26也可以有效地識(shí)別出來(lái)，單元2的識(shí)別效果最差.單元49 和單元50 的不確定性最小，單元2的不確定性最大.

3.2.2 先驗(yàn)為正態(tài)分布時(shí)參數(shù)的反演結(jié)果

采用與以上均勻分布相同的高維參數(shù)模型.假設(shè)病灶單元參數(shù)服從均值為2 Ω·m，標(biāo)準(zhǔn)差為0.05的正態(tài)分布，其余參數(shù)均服從均值為1 Ω·m，標(biāo)準(zhǔn)差為0.05 的正態(tài)分布.選取9 號(hào)、10 號(hào)、25 號(hào)參數(shù)區(qū)域作為病灶單元進(jìn)行參數(shù)反演.各單元先驗(yàn)分布和后驗(yàn)分布的均值結(jié)果如圖10所示.

圖10 先驗(yàn)分布為正態(tài)分布時(shí)的反演結(jié)果示意圖Fig.10 Schematic diagram of inversion results when prior distribution is normal distribution

由圖10 可知，各單元先驗(yàn)分布和后驗(yàn)分布的均值相差不大，經(jīng)算法識(shí)別后的單元均值接近于真值，說(shuō)明DREAM_ZS算法可以識(shí)別出各個(gè)單元參數(shù).

部分反演結(jié)果如圖11、圖12 所示.由圖8、圖9、圖11、圖12 的曲線趨勢(shì)可以看出：參數(shù)的變化波動(dòng)在先驗(yàn)為正態(tài)分布時(shí)比為均勻分布時(shí)要小，曲線更為規(guī)律，更接近于正態(tài)分布，而且各參數(shù)的峰值離MAP 值更近.說(shuō)明參數(shù)的先驗(yàn)分布為正態(tài)時(shí)要比為均勻時(shí)識(shí)別效果更好.

圖12 典型單元后驗(yàn)分布邊緣密度曲線Fig.12 Marginal density curve of posterior distribution of typical unit

在圖11 中，單元25 的后驗(yàn)分布峰值最接近于MAP，其次是單元10 和單元9.說(shuō)明DREAM_ZS 算法對(duì)單元25 的識(shí)別效果更強(qiáng)，對(duì)單元9 和單元10 的識(shí)別效果次之.單元9 和單元10 的不確定性要大于單元25，敏感程度低.在圖12 中，單元49 和單元50 的峰值最接近于MAP，其次是單元26，單元2 的峰值與MAP相差最大.說(shuō)明DREAM_ZS 算法對(duì)單元49和單元50的識(shí)別效果最強(qiáng)，單元26也可以有效地識(shí)別出來(lái)，單元2 的識(shí)別效果最差.單元49 和單元50 的不確定性最小，單元2 的不確定性最大.

為了更直觀地說(shuō)明參數(shù)的后驗(yàn)標(biāo)準(zhǔn)差分布，作出各單元標(biāo)準(zhǔn)差的灰度圖，如圖13所示.

圖13 不同先驗(yàn)分布下后驗(yàn)分布標(biāo)準(zhǔn)差灰度圖Fig.13 Gray scale of standard deviation of posterior distribution under different prior distributions

由圖13 可以直觀地看出：不論先驗(yàn)是均勻分布還是正態(tài)分布，內(nèi)層16 個(gè)單元的標(biāo)準(zhǔn)差均大于外層標(biāo)準(zhǔn)差.這說(shuō)明圓模型的內(nèi)層參數(shù)的不確定性要大于外層.分析可知，這是由于外層參數(shù)單元距離電極更近，對(duì)電極的變化也就更敏感.前16 個(gè)單元后驗(yàn)分布的標(biāo)準(zhǔn)差都要大于先驗(yàn)的標(biāo)準(zhǔn)差，其余單元后驗(yàn)分布的標(biāo)準(zhǔn)差都小于先驗(yàn)的標(biāo)準(zhǔn)差，這說(shuō)明這些單元的后驗(yàn)分布相對(duì)于先驗(yàn)分布變得集中，參數(shù)經(jīng)DREAM_ZS算法識(shí)別降低了不確定性.

4 結(jié)論

1）利用BP 神經(jīng)網(wǎng)絡(luò)模型替代有限元模型完成EIT 正問(wèn)題的計(jì)算，取得了計(jì)算精度高的結(jié)果，且大大提高計(jì)算效率，為EIT參數(shù)反演打下基礎(chǔ).

2）對(duì)于EIT 低維反演問(wèn)題，基于四層同心圓模型，在貝葉斯框架下，采用DREAM_ZS 抽樣算法，對(duì)四個(gè)參數(shù)進(jìn)行了有效地反演，并分析了各個(gè)參數(shù)的不確定性，發(fā)現(xiàn)頭皮的不確定性最小，其次是顱骨，大腦和腦脊液的不確定性較大.對(duì)比不同的電流激勵(lì)模式發(fā)現(xiàn)，相對(duì)激勵(lì)的參數(shù)反演效果最優(yōu).

3）以二維圓模型為仿真算例研究EIT 高維反演問(wèn)題，基于貝葉斯框架的DREAM_ZS 抽樣算法能夠準(zhǔn)確識(shí)別病灶單元.對(duì)比不同的先驗(yàn)分布，發(fā)現(xiàn)正態(tài)分布比均勻分布的反演結(jié)果更優(yōu).另外，單元與電極的相對(duì)位置會(huì)影響其反演結(jié)果和不確定性.當(dāng)先驗(yàn)為正態(tài)分布時(shí)，單元距離電極越近，其參數(shù)不確定性越小，算法識(shí)別程度越高.

4）基于貝葉斯理論的參數(shù)反演方法在EIT 不確定性量化研究中具有可靠性高、適用范圍廣的優(yōu)點(diǎn).在已知參數(shù)不確定性程度的基礎(chǔ)上，如何合理地給出改進(jìn)不確定性較大參數(shù)的方法，這需要在今后的研究中做進(jìn)一步探索.