基于DREAM_ZS 算法的EIT 電阻率反演方法研究

李穎 ，馬重蕾 ，趙營鴿 ，王冠雄 ，郝虎鵬

（1.河北工業大學 生命科學與健康工程學院，天津 300130；2.河北工業大學 河北省生物電磁與神經工程重點實驗室，天津 300130；3.河北工業大學 天津市生物電工與智能健康重點實驗室，天津 300130；4.新鄉醫學院三全學院 智能醫學工程學院，河南 新鄉 453003）

電阻抗成像（EIT）的基本原理是在人體表面施加安全的激勵電流，測量由此電流產生的目標表面的電壓信號，重構出被測目標內部的阻抗圖像分布［1］.由于它具有操作簡單、成本低廉、無損傷等特點，在醫學成像方面獲得了廣泛的關注.在對EIT 的研究中，為了簡化計算，通常假設生物組織的阻抗分布是不變的，但實際情況中，阻抗分布是變化的，這就導致圖像重建結果是不確定的［2］.因此研究電阻率的變化對EIT逆問題的影響具有重要的意義.

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

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

另外，在貝葉斯方法的實現中，對后驗分布的計算需要調用大量與正問題模型相關的似然函數計算，正問題原始模型計算耗時過長會大大影響整體效率，所以采用計算效率高的模型來替代原始模型是很有必要的.常用的替代模型有神經網絡類模型［17-18］、高斯過程回歸模型［19］、克里金（Kriging）模型［20］等.反向傳播（Back Propagation，BP）神經網絡由于其非線性映射能力強、計算效率高，被廣泛應用于復雜非線性模型的替代模型.

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

1 理論與方法

1.1 EIT逆問題的求解模型

EIT研究的是一個特殊的電流場問題，可以將待測目標看作成一個導體，電流的流動受場域內電導率分布的影響［21］.EIT 數學物理模型可以用一個二階橢圓型偏微分方程來表示［22］.設物體所在區域為Ω，邊界為Γ，假設物體內部無源，僅考慮電導率σ分布時，有

EIT 逆問題是根據測得的邊界電壓分布反演內部阻抗分布.由于電流密度是阻抗分布的函數，所以EIT 逆問題是個非線性問題.在實際測量中，電極數量和電流注入方式是有限的，這就導致了從邊界電極獲得的信息不足以確定目標內部的阻抗分布，而且，電極位置對阻抗分布的敏感性不同.這些使得EIT逆問題成為一個高度病態的非線性問題.

由于EIT 逆問題的不適定性和病態性十分嚴重，在求解時需要引入先驗信息.結合了先驗信息和觀測數據的貝葉斯方法可以量化逆問題解的不確定性，有些情況下可以降低不適定問題的不適定程度［10］.

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

貝葉斯理論描述的是由果及因，即利用事情發生后的結果推出導致事情發生的原因.它的數學表達式為

式中，p（m|d）為后驗分布，p（d|m）為似然函數，p（m）為先驗分布，p（d）為測量數據的概率密度.在參數反演問題中m為待反演參數，d為測量數據，由于p（d）與反演參數無關，所以可看成一個常量.那么上式也可以表示為

可見參數的后驗分布完全是由先驗分布以及似然函數決定的.先驗參數分布是由人的日常生活經驗或者科學家及權威科研資料所給出的.似然函數表征的是理論數據與實測數據之間的差異程度.最常用的似然函數是高斯似然函數，即假設誤差向量δi服從獨立標準正態分布，其似然函數計算表達式為

式中，L（x）為似然函數，xi為反演參數，（fxi）為正問題計算結果，即模擬計算數據，udi為實測數據，nd為參數數目.

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

當該馬爾科夫鏈的轉移概率滿足此特性時，說明該鏈已經達到可靠的收斂程度，此時該鏈的分布即為參數的穩定后驗分布.

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

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

2）對第k代每一個個體進行變異操作.

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

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

式中，U是一個取值在0～1 的隨機數，交叉概率CR∈［0，1］，根據抽樣自適應性發生變化.

4）對新樣本計算接受概率.若接受新樣本，則取變異后的樣本點；若不接受，保持原位置狀態，使平行序列繼續進化.

式中：q為用于評價的馬爾科夫鏈的條數；B為鏈內方差；為q條馬爾科夫鏈均值的方差；為uj的均值；W為q條馬爾科夫鏈方差Sj的均值.若值達標，表明馬爾科夫鏈已可靠收斂；若不達標，則繼續進行序列進化.

1.3 基于貝葉斯框架的EIT參數反演

對于EIT 參數反演，電阻率ρ和測量的電壓U被假設為具有某種聯合概率的多元隨機變量.貝葉斯方法的目標是在給定電位測量值的情況下求解電阻率分布的后驗概率.

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

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

1）根據電阻率ρ的先驗分布產生N組初始樣本；

2）將樣本值輸入EIT 正問題計算模型，計算出理論電位值U（ρ）；

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

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

5）采用DREAM_ZS抽樣算法對后驗分布進行抽樣，得出各個電阻率的邊緣后驗密度；

6）對電阻率的反演結果進行不確定性分析.

1.4 EIT正問題計算的BP神經網絡替代模型

在EIT 正問題的研究中，往往采用有限元法進行求解.如果在EIT 模型中反復調用有限元來進行計算，將會導致計算效率很低.BP 神經網絡結構簡單、非線性映射能力強，可以明顯提高計算速度，故本文采用BP神經網絡來替代有限元模型.

BP 神經網絡主要由三部分構成：輸入層、隱含層和輸出層［25］.在神經網絡的訓練過程中進行了兩種計算，一種是前向計算，另一種是逆向反饋計算.

具體地，前向計算過程為

式中，ωij表示節點i與節點j之間的權值，bj為節點j的閾值，xi為節點輸入值，f為激活函數，yj為節點輸出值.

逆向反饋過程為

式中，dj為假定的輸出結果，e(ω，b)表示誤差函數，表示修正后的權值，表示修正后的閾值.當達到目標誤差期望函數時，BP神經網絡訓練完成.

2 仿真模型的構建與正問題計算

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

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

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

圖2 高維參數模型Fig.2 High dimensional parameter model.

為了提高有限元模型的計算效率，我們采用BP神經網絡模型作為有限元的替代模型.首先將有限元模型的電阻率值作為網絡的輸入，計算出的電極電位結果作為網絡的輸出，隨機生成1 000組輸入輸出數據，選擇數據的70%作為網絡的訓練數據、15%作為修正數據、15%作為測試數據，訓練BP 人工神經網絡.

利用均方差來度量BP 神經網絡的性能，均方誤差MSE的計算公式為

式中，N為樣本個數，Yi為驗證集輸出結果，yi為測試集輸出結果.圖3 顯示了BP 神經網絡訓練過程中的均方差變化.

圖3 BP神經網絡均方差Fig.3 BP neural network variance diagram.

由圖3 可以看出各個測試集最終都趨于收斂，且在迭代完成時四層同心圓模型的均方差低于10-8，高維參數模型的均方差低于10-5，可見BP 神經網絡擬合得很好.

進一步對訓練好的網絡進行驗證，按均值±30%隨機生成20 組四層電阻率值，分別代入有限元模型以及BP 神經網絡模型.計算16 個電極20 組數據的平均相對誤差值，如圖4 所示，由于9 號電極為參考電極，圖中忽略不計.

圖4 平均相對誤差Fig.4 Average relative error.

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

3 參數反演結果

3.1 基于同心圓模型的低維參數反演

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

表2 各激勵方式下參數后驗分布統計信息表Tab.2 Statistical information table of posterior distribution of parameters under various excitation modes

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

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

三種電流激勵方式下的電阻率參數后驗分布的邊緣概率密度曲線如圖5 所示，圖中藍色叉點為MAP值.

圖5 各激勵方式下參數后驗分布邊緣概率密度曲線Fig.5 Marginal density curve of posterior distribution of parameters under various excitation modes

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

為了定量描述參數的后驗分布信息，計算各激勵方式下標準差.其中標準差計算為

圖6 不同激勵方式下參數的標準差Fig.6 Standard deviation of parameters under different excitation modes

從圖6 可以看出，對于三種激勵模式，頭皮的標準差均最小，其次是顱骨，而大腦和腦脊液的標準差大一些.說明四個參數中頭皮的后驗分布最集中，不確定性最小，敏感性最強.其次是顱骨，大腦和腦脊液的不確定性較大.

3.2 基于二維圓模型的高維參數反演

通過前面的結論可以得知在相對激勵方式下對逆問題電阻率反演效果最好，所以在高維參數反演中電流采用相對激勵方式注入.根據貝葉斯理論可知，先驗分布的不同在一定程度上會影響后驗結果的分布，所以接下來研究電阻率在不同先驗分布下（即均勻分布與正態分布下）參數的反演結果.

3.2.1 先驗為均勻分布時參數的反演結果

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

圖7 先驗為均勻分布時的反演結果示意圖Fig.7 Schematic diagram of inversion results when prior distribution is uniform distribution

由圖7 可知，各單元先驗分布和后驗分布的均值相差不大，經算法識別后的單元均值接近于真值，說明DREAM_ZS算法可以識別出各個單元參數.

圖8 顯示了單元9、單元10、單元25 的后驗分布邊緣密度曲線.對于非病灶單元，由于單元數目過多，我們選取幾個典型單元的后驗分布邊緣密度曲線展示于圖9中.

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

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

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

在圖9 中，單元49 和單元50 的峰值最接近于MAP，其次是單元26，單元2 的峰值與MAP 相差最大.說明DREAM_ZS算法對單元49和單元50的識別效果最強，單元26也可以有效地識別出來，單元2的識別效果最差.單元49 和單元50 的不確定性最小，單元2的不確定性最大.

3.2.2 先驗為正態分布時參數的反演結果

采用與以上均勻分布相同的高維參數模型.假設病灶單元參數服從均值為2 Ω·m，標準差為0.05的正態分布，其余參數均服從均值為1 Ω·m，標準差為0.05 的正態分布.選取9 號、10 號、25 號參數區域作為病灶單元進行參數反演.各單元先驗分布和后驗分布的均值結果如圖10所示.

圖10 先驗分布為正態分布時的反演結果示意圖Fig.10 Schematic diagram of inversion results when prior distribution is normal distribution

由圖10 可知，各單元先驗分布和后驗分布的均值相差不大，經算法識別后的單元均值接近于真值，說明DREAM_ZS算法可以識別出各個單元參數.

部分反演結果如圖11、圖12 所示.由圖8、圖9、圖11、圖12 的曲線趨勢可以看出：參數的變化波動在先驗為正態分布時比為均勻分布時要小，曲線更為規律，更接近于正態分布，而且各參數的峰值離MAP 值更近.說明參數的先驗分布為正態時要比為均勻時識別效果更好.

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

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

為了更直觀地說明參數的后驗標準差分布，作出各單元標準差的灰度圖，如圖13所示.

圖13 不同先驗分布下后驗分布標準差灰度圖Fig.13 Gray scale of standard deviation of posterior distribution under different prior distributions

由圖13 可以直觀地看出：不論先驗是均勻分布還是正態分布，內層16 個單元的標準差均大于外層標準差.這說明圓模型的內層參數的不確定性要大于外層.分析可知，這是由于外層參數單元距離電極更近，對電極的變化也就更敏感.前16 個單元后驗分布的標準差都要大于先驗的標準差，其余單元后驗分布的標準差都小于先驗的標準差，這說明這些單元的后驗分布相對于先驗分布變得集中，參數經DREAM_ZS算法識別降低了不確定性.

4 結論

1）利用BP 神經網絡模型替代有限元模型完成EIT 正問題的計算，取得了計算精度高的結果，且大大提高計算效率，為EIT參數反演打下基礎.

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

3）以二維圓模型為仿真算例研究EIT 高維反演問題，基于貝葉斯框架的DREAM_ZS 抽樣算法能夠準確識別病灶單元.對比不同的先驗分布，發現正態分布比均勻分布的反演結果更優.另外，單元與電極的相對位置會影響其反演結果和不確定性.當先驗為正態分布時，單元距離電極越近，其參數不確定性越小，算法識別程度越高.

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