劉海，閆丹丹，荊會(huì)敏

中國計(jì)量大學(xué)信息工程學(xué)院，浙江杭州 310018；

生物電阻抗斷層成像（electrical impedance tomography，EIT）是繼形態(tài)、結(jié)構(gòu)成像后的新一代無創(chuàng)功能性成像技術(shù)。它通過向體表注入較小的安全交變電流，由體表電極陣列測(cè)量得到體外邊界電壓（電勢(shì)），并采用重建算法重建組織內(nèi)部電阻率分布[1]。與正常的生物阻抗相比，生物組織在發(fā)生病變時(shí)，其阻抗值變化十分顯著[2]。EIT技術(shù)的實(shí)質(zhì)是利用反映組織生理及功能狀態(tài)的電阻抗信息進(jìn)行阻抗重建。因此，EIT技術(shù)可用于腫瘤、癲癇等多種疾病的早期發(fā)現(xiàn)與診斷；并且EIT技術(shù)不產(chǎn)生電離輻射、對(duì)人體無害、相對(duì)成本較低、可多次重復(fù)使用、便攜，能對(duì)患者進(jìn)行長(zhǎng)期、實(shí)時(shí)動(dòng)態(tài)監(jiān)護(hù)?；谝陨咸攸c(diǎn)，盡管目前應(yīng)用EIT技術(shù)與主流醫(yī)學(xué)成像技術(shù)在圖像分辨率上還有一定差距，但其在醫(yī)學(xué)工程及相關(guān)領(lǐng)域仍引起極大的關(guān)注[3-4]。

目前主要的重建算法有修正的牛頓-拉夫遜法、牛頓單步誤差重構(gòu)法、層析法、雙限定法等，上述方法收斂速度較快，但計(jì)算量大且過程復(fù)雜。近年來基于粒子群算法、神經(jīng)網(wǎng)絡(luò)算法和遺傳算法等相關(guān)的智能尋優(yōu)算法逐漸應(yīng)用于EIT圖像重建，這些智能算法與傳統(tǒng)算法相比，計(jì)算量較小、較易修正局部最優(yōu)值，具有較快的收斂速度和全局搜索能力等優(yōu)勢(shì)[5]。

本文采用改進(jìn)的差分進(jìn)化（differential evolution，DE）算法對(duì)頭部電阻抗分布進(jìn)行重建，在三維頭球模型上進(jìn)行仿真實(shí)驗(yàn)，并與其他進(jìn)化類算法進(jìn)行對(duì)比。

1 EIT的基本原理

在EIT研究中，電阻抗分布與電流之間的關(guān)系由麥克斯韋方程給出。EIT測(cè)量中采用較低（10～100 kHz）的激勵(lì)源頻率，可忽略介電常數(shù)的影響，此時(shí)電流場(chǎng)可當(dāng)作穩(wěn)態(tài)電流場(chǎng)。為簡(jiǎn)化計(jì)算，忽略電極與皮膚間的接觸阻抗，電壓與電流間的非線性關(guān)系可由滿足相應(yīng)邊界條件的拉普拉斯方程確定，其求解場(chǎng)域的數(shù)學(xué)模型[6-7]見公式（1），應(yīng)滿足的第一類或第二類邊界條件（、）見公式（2）。

其中，表示散度，表示梯度，ρ為目標(biāo)內(nèi)部電阻率分布，φ為場(chǎng)域中電位分布函數(shù)，φ0為邊界電位，n為方向向量，為邊界電位分布的法向?qū)?shù)，Jn為電流密度，無電流注入時(shí)為0。

EIT正問題是指由人體內(nèi)部阻抗分布及邊界條件，求解體表或內(nèi)部電壓分布，這是求解逆問題的基礎(chǔ)。EIT逆問題是指由體表電壓分布及邊界條件，求解人體內(nèi)部阻抗分布。經(jīng)過有限元法離散后，EIT逆問題可轉(zhuǎn)化為公式（3）中求電位分布?xì)埐畹南鄬?duì)均方誤差最小值[RMS(ρ)]問題。

其中，上標(biāo)的T表示矩陣的轉(zhuǎn)置，ρ為目標(biāo)內(nèi)部電阻率分布，φ(ρ)是邊界電位的計(jì)算值，φ0是邊界電位“測(cè)量值”。

因此，EIT逆問題是找到合適的電阻抗分布ρ使目標(biāo)函數(shù)RMS(ρ)達(dá)到最小，即搜尋一個(gè)m維參數(shù)矢量使RMS(ρ)最小，其中m為有限元剖分單元數(shù)，即電阻率分布自由度。從而，EIT逆問題即轉(zhuǎn)化為多參數(shù)函數(shù)優(yōu)化問題。

2 基于改進(jìn)DE的EIT重建算法

2.1 基本 DE的算法 DE算法是一種基于種群進(jìn)化的全局啟發(fā)式算法[8]，它在許多優(yōu)化問題上的效果均優(yōu)于遺傳算法[9]、粒子群算法[10]、模擬退火算法[11]等。該算法一般包括變異、交叉、選擇3個(gè)基本操作。首先在整個(gè)搜索空間隨機(jī)生成一組初始種群，然后對(duì)種群中每個(gè)個(gè)體進(jìn)行變異、交叉操作得到臨時(shí)種群，最后用臨時(shí)種群中適應(yīng)性更好的個(gè)體替代原始種群中相應(yīng)的個(gè)體。整個(gè)算法的步驟見圖1。

圖1 DE基本算法流程

在整個(gè)算法迭代過程中，變異是其中的核心過程，它采用隨機(jī)偏差擾動(dòng)生成變異向量的方式能夠增加種群多樣性，增強(qiáng)算法全局搜索能力。根據(jù)生成變異向量方法的不同，能夠組合成不同的變異策略，其中基本變異策略DE/rand/1/bin方程見公式（4）。

此外，其他典型的變異策略 DE/best/1/bin和DE/rand-to-best/1/bin的方程分別見公式（5）、（6）。

其中，νi,G+1表示生成的變量向量，i表示種群個(gè)體序列號(hào)，G表示種群進(jìn)化代數(shù)，r1≠r2≠r3≠i，且i、r1、r2、r3∈[1,NP]，NP表示種群規(guī)模，xbest,G為第G代種群中的最優(yōu)個(gè)體，F(xiàn)、λ是縮放因子，為常數(shù)。

2.2 改進(jìn)DE算法 本文在基本DE算法的基礎(chǔ)上，引入隨機(jī)最佳變異和局部增強(qiáng)算子，提出一種能在保證種群多樣性前提下加速收斂、提高局部搜索能力的改進(jìn)DE算法。具體實(shí)現(xiàn)如下：

2.2.1 隨機(jī)最佳變異 在基本 DE算法中，為克服DE/rand/*/*變異策略收斂速度較慢的缺點(diǎn)，普遍采用DE/best/*/*變異策略，在變異過程中，種群中個(gè)體均采用相同的最佳個(gè)體，對(duì)加快算法收斂性有良好的引導(dǎo)作用，但最佳個(gè)體周圍集中大量的尋優(yōu)個(gè)體，種群多樣性遭到破壞，容易導(dǎo)致算法局部收斂。

為使變異向量具有較好的多樣性和引導(dǎo)性，本文采用一種隨機(jī)最佳變異策略，見公式（7）。

其中，r1≠r2≠r3≠i，且i、r1、r2、b∈[1,NP]，RMS(xb,G)≤min{RMS(xr1,G), RMS(xr2,G)}。

隨機(jī)最佳變異策略將3個(gè)隨機(jī)選擇的個(gè)體中的局部最優(yōu)個(gè)體當(dāng)作變異算子的基向量，而另外2個(gè)向量作為差向量對(duì)，這樣既保證種群多樣性，又為種群優(yōu)化提供了一定的引導(dǎo)作用，避免種群大量聚集于最佳個(gè)體周圍，從而避免了算法的“早熟”。

2.2.2 局部增強(qiáng)算子 引入隨機(jī)最佳變異策略，保證了算法收斂速度和種群多樣性之間的相對(duì)平衡，但算法收斂速度相對(duì)較慢，為進(jìn)一步加快算法收斂，引入局部增強(qiáng)算子，對(duì)除最佳個(gè)體外的其他種群個(gè)體以比例因數(shù)CP（0＜CP＜1）的概率進(jìn)行重新賦值，使部分個(gè)體向最佳個(gè)體靠攏，以增強(qiáng)這些個(gè)體的貪婪性，加快算法收斂速度，局部增強(qiáng)算子見公式（8）。

其中，r1≠r2≠i，且i、r1、r2∈[1,NP]，?i,G+1是增強(qiáng)后的新個(gè)體，用來代替?zhèn)€體xi,G+1，xbest,G+1是新種群中的最佳個(gè)體，g為算法迭代次數(shù)，改進(jìn)算法引入新的比例因數(shù)CP，當(dāng)CP=1時(shí)，除最佳個(gè)體外，其余所有個(gè)體向量都需重新賦值；當(dāng)CP=0時(shí)，未選中個(gè)體進(jìn)行局部增強(qiáng)。

局部增強(qiáng)算子的目的是在增加部分個(gè)體貪婪性的同時(shí)保證種群多樣性，以使算法達(dá)到既快又好地逼近到全局最優(yōu)解。算法在增加局部增強(qiáng)算子后，每次迭代僅按一定比例使部分個(gè)體變得更加貪婪，從而在一定程度上限制了算法的貪婪性，避免影響算法的收斂性。需注意的是，隨著迭代過程的進(jìn)行，種群中這部分個(gè)體隨機(jī)擾動(dòng)范圍也在動(dòng)態(tài)縮小，這種方式對(duì)算法局部搜索能力有較大的提高，尤其是能減少收斂到全局最優(yōu)解所需的迭代次數(shù)，但算法也有陷入局部最優(yōu)的可能。

3 仿真實(shí)驗(yàn)與結(jié)果

3.1 創(chuàng)建仿真模型 在ANSYS 10.0軟件上采用有限元法對(duì)頭部進(jìn)行建模，并對(duì)EIT正問題進(jìn)行求解。三層有限元頭球模型如圖 2A，由 14 332個(gè)單元和20 193個(gè)節(jié)點(diǎn)組成。頭部模型分為頭皮、顱骨、大腦及病變組織4個(gè)部分[12]。為簡(jiǎn)化計(jì)算，假定頭部模型為各向同性均質(zhì)導(dǎo)體，且電阻率值已知。根據(jù)先驗(yàn)知識(shí)，將模型 4個(gè)部分的電阻率按比例設(shè)定為：頭皮∶顱骨∶大腦∶病變組織=1∶15∶1∶2[13]，其目標(biāo)電阻率分布的剖面圖如圖2B，其中淺色區(qū)域?yàn)椴∽兘M織。

圖 2 采用有限元法對(duì)頭部進(jìn)行建模。A.頭球模型；B.目標(biāo)電阻率分布

3.2 仿真參數(shù)設(shè)定 本文對(duì)三維頭球模型進(jìn)行仿真實(shí)驗(yàn)，并采用相對(duì)電極注入方式向模型左右注入±5 mA的電流。在逆問題重建實(shí)驗(yàn)中，將頭球模型的頭皮、顱骨、大腦及病變組織4個(gè)部分的上限設(shè)為[1.20，17.40，1.20，2.25]，下限設(shè)為[0.81，13.60，0.81，1.75]，初始種群需從此約束邊界范圍內(nèi)隨機(jī)生成。由于大腦和頭皮電阻率的設(shè)定值相等，為簡(jiǎn)化計(jì)算，假定變量個(gè)數(shù)D=3，種群規(guī)模NP=60，目標(biāo)函數(shù)為公式（3）。算法終止條件為迭代次數(shù)達(dá)到100次。根據(jù)文獻(xiàn)將比例因數(shù)F、交叉因數(shù)CR分別設(shè)定為0.8、0.9[13]。局部增強(qiáng)算子的比例因數(shù)CP等到下文分析此參數(shù)對(duì)重建結(jié)果的影響后再另行確定。

3.3 仿真實(shí)驗(yàn)與結(jié)果 為了對(duì)算法圖像重建質(zhì)量進(jìn)行定量分析，分別引入相對(duì)誤差RE和誤差總和TE準(zhǔn)則，見公式（9）、（10）。

其中，ρ*是指模型各個(gè)部分的重建電阻率值，ρ0是指相應(yīng)的模型各個(gè)部分的目標(biāo)電阻率值。

其中，m為電阻率分布自由度，ρitar和ρirec分別代表單元i的目標(biāo)電阻率值和重建電阻率值。

圖3為改進(jìn)DE算法的重建圖像誤差總和TE隨局部增強(qiáng)因子的比例因數(shù)CP變化的曲線，當(dāng)CP=5/60時(shí)，重建圖像的整體誤差最小，算法效果最佳；當(dāng)CP為1/60～5/60時(shí)，在相同的迭代次數(shù)下，隨著CP取值的增加，重建圖像的誤差逐漸減小，說明當(dāng)聚集在最優(yōu)個(gè)體周圍的個(gè)體在一定范圍內(nèi)時(shí)，會(huì)在一定程度上加快算法收斂的速度，提高重建圖像的精度；當(dāng)CP為5/60～15/60時(shí)，隨著CP取值的增加，重建圖像的誤差反而逐漸增加，說明當(dāng)聚集在最優(yōu)個(gè)體周圍的個(gè)體超出一定范圍后，算法的貪婪性過大，反而使算法陷入“早熟”。因此建議采用此改進(jìn)方法進(jìn)行三維腦部阻抗圖像重建時(shí)，將CP的取值范圍確定在 3/60～7/60，即可能獲得更好的圖像重建效果。在本文的仿真對(duì)比實(shí)驗(yàn)中選擇局部增強(qiáng)因子CP=4/60。

圖3 誤差總和隨局部增強(qiáng)因子的比例因數(shù)變化的關(guān)系

采用上述實(shí)驗(yàn)參數(shù)設(shè)置，在相同的初始條件下，將本文提出的改進(jìn)算法與其他算法進(jìn)行仿真對(duì)比實(shí)驗(yàn)。不同算法的電阻抗重建圖像結(jié)果見圖4。

圖4 不同算法電阻率重建分布

其中算法DE-rand、DE-best分別為采用文獻(xiàn)[13]中的DE/rand/1/bin變異策略和DE/best/1/bin變異策略，算法PSO采用文獻(xiàn)[14]中的粒子群優(yōu)化算法。比較圖4及圖2B發(fā)現(xiàn)，這些算法均能在電阻抗重建過程中得到較為清晰的圖像，且能對(duì)病變組織進(jìn)行準(zhǔn)確定位，但不同算法的電阻率重建圖像與目標(biāo)電阻率的分布圖像僅在顏色深淺上有細(xì)微差別，用肉眼無法看出明顯差異，因此需要對(duì)阻抗重建結(jié)果進(jìn)行定量分析。

表 1是經(jīng)過 100輪迭代后不同算法的圖像重建結(jié)果，其中t表示各個(gè)算法單次循環(huán)所需平均時(shí)間，其單位為分鐘。由表1可見，改進(jìn)算法相對(duì)于算法DE-rand、PSO、DE-best而言，其模型各個(gè)部分整體的重建誤差總和TE分別降低了66.996%、32.723%、49.834%；其單次循環(huán)所需平均時(shí)間分別降低了 8.331%、2.017%、4.793%，表明在相同條件下改進(jìn)算法圖像重建精度更高、迭代效率更高，算法整體性能更好。

表1 不同算法的圖像重建實(shí)驗(yàn)結(jié)果

由表2可見，與其他部分相比，各個(gè)算法對(duì)病變組織部分的重建誤差更大，表明在求解EIT正問題時(shí)，大腦顱骨會(huì)阻礙一部分電流的流入，導(dǎo)致病變部分的重建誤差變大，大體上符合理論情況。從總體上看，改進(jìn)算法相比算法 DE-rand、PSO、DE-best，其模型的大腦、顱骨、頭皮、病變部分電阻率重建誤差分別降低86.794%、68.690%、86.720%、86.874%，-16.084%、33.431%、-18.571%、81.164%，80.378%、53.955%、80.238%、70.146%，表明改進(jìn)算法有效降低了圖像重建的誤差，提高了圖像重建的精度，算法性能得到一定的提升。

表2 頭球模型不同部分的圖像重建相對(duì)誤差（%）

由圖 5可見，所有算法在迭代初期收斂速度較快，隨著迭代次數(shù)的增加，收斂速度慢慢趨于平緩。當(dāng)?shù)螖?shù)＜18時(shí)，改進(jìn)DE算法較其他算法無明顯優(yōu)勢(shì)，說明改進(jìn)算法在種群多樣性和算法收斂速度上保持了較好的水平；當(dāng)?shù)螖?shù)為 18～27時(shí)，改進(jìn)算法的局部增強(qiáng)算子逐漸發(fā)揮作用，使改進(jìn)算法的重建效果逐漸優(yōu)于算法1和算法2；當(dāng)算法迭代次數(shù)為18～100時(shí)，改進(jìn)算法明顯優(yōu)于其他算法的效果和精度。

圖5 種群平均目標(biāo)函數(shù)值隨迭代次數(shù)變化關(guān)系

4 結(jié)束語

電阻抗圖像重建算法是 EIT研究的關(guān)鍵技術(shù)之一。本文提出了一種新的用于電阻抗重建的改進(jìn) DE算法。該算法通過引入隨機(jī)最佳變異策略和局部增強(qiáng)算子，在算法收斂速度和種群多樣性之間取得良好的平衡，并在一定程度上加快了算法收斂速度，增強(qiáng)了算法局部搜索能力。結(jié)果表明，改進(jìn)算法能夠獲得較為清晰的電阻抗重建圖像，且能對(duì)病變組織進(jìn)行準(zhǔn)確定位。同時(shí)，與其他算法相比，改進(jìn)算法的圖像重建誤差減小30%以上，有效提高了算法的整體性能，為EIT臨床研究提供了一定的參考。

由于DE算法是一種相對(duì)較新的全局類智能尋優(yōu)算法，仍有較多的相關(guān)理論問題亟需完善。改進(jìn)算法由于每次均需對(duì)所有種群個(gè)體進(jìn)行目標(biāo)函數(shù)評(píng)價(jià)，仍存在測(cè)量數(shù)據(jù)、計(jì)算時(shí)間較長(zhǎng)的缺點(diǎn)。因此，如何在一定的實(shí)驗(yàn)條件下增加測(cè)量數(shù)據(jù)并減少仿真計(jì)算時(shí)間是下一步研究的重點(diǎn)。

[1] Silvera-Tawil D, Rye D, Soleimani MA. Electrical impedance tomography for artificial sensitive robotic skin: a review. IEEE Sens J, 2015, 15(4): 2001-2016.

[2] 陳會(huì), 閆丹丹. 基于磁共振技術(shù)的生物電阻抗成像研究進(jìn)展. 中國醫(yī)學(xué)影像學(xué)雜志, 2014, 22(9): 714-717.

[3] Wagenaar J, Adler A. Electrical impedance tomography in 3D using two electrode planes: characterization and evaluation. Physiol Meas, 2016, 37(6): 922-937.

[4] Schullcke B, Gong B, Moeller K. Steps towards 3D electrical impedance tomography. Conf Proc IEEE Eng Med Biol Soc, 2015: 5323-5326.

[5] 萇飛霸, 張和華, 顏樂先, 等. 電阻抗斷層成像技術(shù)研究.中國醫(yī)療器械雜志, 2016, 40(1): 52-54.

[6] Trepte CJ, Phillips CR, Solà J, et al. Electrical impedance tomography (EIT) for quantification of pulmonary edema in acute lung injury. Crit Care, 2016, 20(1): 18.

[7] 張輝, 李穎, 王西明, 等. 電阻抗斷層成像的MPSOMNR算法研究. 計(jì)算機(jī)工程與應(yīng)用, 2013, 49(9): 29-32.

[8] Ribeiro RR, Feitosa AS, Souza RD, et al. A modified differential evolution algorithm for the Reconstruction of electrical impedance tomography images//Biosignals and Biorobotics Conference. Salvador: IEEE, 2014: 1-6.

[9] Rashid A, Kim S, Liu D, et al. A dynamic oppositional biogeography-based optimization approach for time-varying electrical impedance tomography. Physiol Meas, 2016,37(6): 820-842.

[10] Martin S, Choi CT. Nonlinear electrical impedance tomography reconstruction using artificial neural networks and particle swarm optimization. IEEE Trans Magn, 2016,52(3): 1-4.

[11] Zhang XL, Wang W, Sze G, et al. An image reconstruction algorithm for 3-D electrical impedance mammography.IEEE Trans Med Imaging, 2014, 33(12): 2223-2241.

[12] Hamilton SJ, Lassas M, Siltanen S. A direct reconstruction method for anisotropic electrical impedance tomography.Inverse Probl, 2014, 30(7): 1375-1378.

[13] 閆丹丹, 陳會(huì). 基于電阻抗成像技術(shù)的腦病變組織檢測(cè)仿真研究. 中國醫(yī)學(xué)影像技術(shù), 2014, 30(7): 1113-1116.

[14] 陳民鈾, 楊艷利, 何為, 等. 基于粒子群優(yōu)化算法的電阻抗圖像重建. 重慶大學(xué)學(xué)報(bào), 2011, 34(1): 82-87.