劉 海，閆丹丹，荊會敏

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

基于微分進化的EIT圖像重建控制參數(shù)研究

劉 海，閆丹丹，荊會敏

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

微分進化算法主要有三個隨機參數(shù)：種群大小(NP),縮放因子(F),交叉因數(shù)(CR). 這些參數(shù)的取值對EIT圖像重建效果的好壞起著重要的作用. 但當(dāng)前微分進化算法參數(shù)選擇具有隨機性，大多數(shù)的參數(shù)研究是通過標(biāo)準(zhǔn)測試函數(shù)進行，沒有具體到特定的領(lǐng)域. 針對這些問題，文章以頭部EIT圖像重建為例，在給定目標(biāo)函數(shù)和終止條件的基礎(chǔ)上，通過大量的仿真實驗，分析了各個參數(shù)對圖像重構(gòu)結(jié)果的影響，并給出了這些參數(shù)的合理選取區(qū)間，從而為微分進化算法在EIT圖像重建中的應(yīng)用提供了有效的依據(jù).

電阻抗成像；微分進化算法；有限元模型；參數(shù)設(shè)置

ControlparametersstudyofEITimagereconstruction

Abstract: Population size(NP), scaling factor (F) and crossover factor (CR) are the main random parameters of a differential evolution algorithm. The value of the parameters plays an important role in the EIT image reconstruction effect. However, the current parameter selection is random and nonspecific. To solve these problems, based on given objective functions and termination conditions, and a large number of simulation experiments of the head EIT image reconstruction, the influence of each parameter on image reconstruction was analyzed, and the reasonable selection interval of these parameters was given. It provides an effective basis for the application of differential evolution algorithms in EIT image reconstruction.

Keywords: electrical impedance tomography(EIT); differential evolution(DE) algorithm; finite element model; parameter setting

電阻抗斷層成像(electrical impedance tomography, EIT)是現(xiàn)今國內(nèi)外生物醫(yī)學(xué)領(lǐng)域的重要研究課題之一，是一種非侵入性的無損功能性成像技術(shù). 首先將對人體安全的交變電流注入到附著在被測組織表面的電極上，然后測量得到相應(yīng)的邊界電位，最后通過重構(gòu)算法求解邊界值問題，重構(gòu)出被測組織內(nèi)部的電阻率分布[1]. 與傳統(tǒng)的醫(yī)療成像技術(shù)相比，EIT技術(shù)能在疾病發(fā)生初期檢測到人體組織電阻率分布的變化，這有助于疾病的早期發(fā)現(xiàn)與診斷，并且可多次反復(fù)使用、無電離輻射[2]、成本相對較低、能長期對疾病進行監(jiān)測和觀察，因而具有較大潛力，應(yīng)用領(lǐng)域十分廣泛，是一種具有較大研究價值的醫(yī)學(xué)成像技術(shù)[3-4].

EIT圖像重構(gòu)問題是一個高度病態(tài)的非線性逆問題. 國內(nèi)外研究者提出了許多經(jīng)典的圖像重構(gòu)算法，如修正牛頓-拉夫遜法、雙限定法、剝層算法等. 目前修正的牛頓-拉夫遜法是最成熟的重構(gòu)算法之一，應(yīng)用十分廣泛. 但此方法對初值的選擇較為苛刻，一旦初值選擇不當(dāng)，則較易陷入局部最優(yōu)，并且需進行恰當(dāng)?shù)恼齽t化處理并計算雅克比矩陣，計算復(fù)雜，繁瑣. 近來研究者逐漸將智能優(yōu)化類算法(如神經(jīng)網(wǎng)絡(luò)算法、粒子群算法、遺傳算法等)應(yīng)用到EIT圖像重構(gòu)中，由于這些算法能直接對目標(biāo)函數(shù)進行評價，無需求微分，對函數(shù)是否連續(xù)沒有要求等優(yōu)勢而獲得了廣泛的應(yīng)用[5].

微分進化(differential evolution, DE)是一種基于種群進化的全局啟發(fā)式算法，易于實現(xiàn)[6]，控制參數(shù)較少，且在很多優(yōu)化類的問題上都有優(yōu)于其他算法的效果. DE算法的幾個主要控制參數(shù)為：種群大小(NP)、交叉因數(shù)(CR)、縮放因子(F). DE性能對控制參數(shù)的選取非常敏感，但在DE中很難獲得這些控制參數(shù)值的最佳組合，并且在不同的應(yīng)用問題中需要不同的參數(shù)組合，所以需要大量反復(fù)的實驗[7]. 研究表明，DE中動態(tài)參數(shù)控制有優(yōu)于恒定參數(shù)控制的效果. 這些參數(shù)控制中的大多數(shù)是通過標(biāo)準(zhǔn)測試函數(shù)進行分析的，很少具體到特定的領(lǐng)域[6,8]. 因而EIT中DE算法的參數(shù)選取還沒有嚴(yán)格的理論依據(jù). 在本研究中，采用DE算法進行EIT圖像重建，對影響圖像重建質(zhì)量和重建效率的DE算法的三個控制參數(shù)進行了詳細(xì)的分析，同時對參數(shù)初始合理范圍的選取給出了一些建議.

1 EIT原理

EIT圖像重構(gòu)理論上是一個低頻電流場計算的逆問題. 通常忽略皮膚與電極間的接觸阻抗，在EIT測量中激勵源的頻率較低，可將電流場看作穩(wěn)態(tài)電流場. 場域內(nèi)的電阻率分布函數(shù)ρ與電位分布φ由滿足相應(yīng)邊界條件的拉普拉斯方程確定[9]：

▽·ρ-1(x,y)▽φ(x,y)=0 (x,y)∈Ω.

(1)

其邊界條件為：

φ(x,y)=φ0(x,y) (x,y)∈?Ω.

(2)

(3)

其中?Ω為區(qū)域Ω的邊界，ρ為目標(biāo)內(nèi)部電阻率分布，φ為場域內(nèi)電位分布，φ0為已知的邊界電位，n是區(qū)域邊界外法線方向，j為流入?yún)^(qū)域的電流密度.

EIT正問題是在已知電阻率分布ρ以及激勵電流模式j(luò)的情況下求解邊界電位φ，而逆問題則是由已知邊界處的φ0和j求解目標(biāo)內(nèi)部電阻率分布.

在正問題求解中，本實驗通過有限元分析軟件ANSYS進行建模、剖分，采用相對注入方式對模型注入安全電流，然后通過ANSYS“測量”出邊界電位值；在逆問題求解中，DE算法通過不斷地迭代求解使目標(biāo)函數(shù)達(dá)到最小的電阻率分布，目標(biāo)函數(shù)通常選為電位計算值與電位“測量”值之間殘差的函數(shù)，其公式如下：

(4)

式中(·)T表示矩陣的轉(zhuǎn)置，φ(ρ)是邊界電位計算值，φ0是邊界電位“測量值”.

2 DE算法

DE算法是一種基于群體的隨機啟發(fā)式搜索算法. 其基本思想是：首先由隨機產(chǎn)生的初始種群開始，然后由父代個體間的變異操作產(chǎn)生變異個體，接著父代個體和子代個體之間按一定的概率進行交叉操作，生成實驗個體，最后試驗個體與父代個體之間依據(jù)目標(biāo)函數(shù)值的大小進行貪婪選擇操作，保留較優(yōu)的個體，從而實現(xiàn)種群的進化. 其具體過程如下：

Step 1 種群初始化. DE算法一般采用隨機函數(shù)產(chǎn)生NP個維數(shù)為D的初始種群xi,G，i= 1，2，…，NP，其中i表示種群個體序列號，G表示種群進化代數(shù)，NP表示種群規(guī)模.DE算法的初始種群是在給定邊界約束內(nèi)隨機產(chǎn)生的，然后通過變異、交叉、選擇等操作產(chǎn)生下一代種群.

Step 2 變異. 對第G代種群中的每個目標(biāo)向量xi,G，根據(jù)生成變異向量方法的不同，能夠組合成不同的變異策略，其中基本的變異策略DE/rand/1/bin方程為

vi,G+1=xr1,G+F·(xr2,G-xr3,G).

(5)

其中r1 ≠r2 ≠r3 ≠i，且r1,r2,r3∈[1,NP]，F(xiàn)是縮放因子.

Step 3 交叉. 交叉操作是將變異操作得到的變異向量vi,G+1和當(dāng)前目標(biāo)向量xi,G按一定規(guī)則進行雜交得到實驗向量ui,G+1，交叉操作過程如式(6)：

uij,G+1=

(6)

其中i=1, 2, …,NP，j=1, 2,…,D，randb(j)是(0, 1)之間均勻分布的隨機數(shù)，rnbr(i)是[1,D] 之間隨機選擇的整數(shù)序列，CR是交叉因數(shù).

Step 4 選擇. DE算法根據(jù)貪婪準(zhǔn)則決定實驗向量ui,G+1是否成為下一代個體成員. 對實驗向量和相應(yīng)目標(biāo)向量進行目標(biāo)函數(shù)評價，然后根據(jù)式(7)決定是否將實驗向量確定為下一代種群個體成員.

(7)

其中f(·)表示其目標(biāo)函數(shù)值.

Step 5 邊界條件的處理. 在大多數(shù)的問題中，必須確保新產(chǎn)生個體位于邊界約束內(nèi).為保證這一點，用在可行域中隨機生成的個體向量替代超出邊界約束的新個體.

3 仿真實驗與結(jié)果

3.1 創(chuàng)建仿真模型

由于人體頭部近似球體，對二維EIT研究，可采用同心圓模型進行仿真實驗研究 本文在ANSYS10.0軟件上采用有限元法對頭部進行二維建模，其三層同心圓模型如圖1(a)，該模型由430個三角形單元和901個節(jié)點組成.模型由頭皮、顱骨、大腦、病變組織四個部分[10]組成. 為簡化計算，將該二維模型看成各向同性均質(zhì)導(dǎo)體，且電阻率值已知. 在仿真試驗中，將四個部分的電阻率值設(shè)置為：頭皮∶顱骨∶大腦∶病變組織=1∶15∶1∶2[11]. 其目標(biāo)電阻率分布如圖1(b)，其中淺色區(qū)域為病變組織.

圖1 仿真模型Figure 1 Simulation model

3.2 仿真實驗條件

本文采用圖1(a)所示的三層同心圓模型進行仿真實驗，采用相對電極注入方式在模型的左右兩側(cè)注入±5 mA的電流. 在圖像逆問題重構(gòu)實驗中，將此模型各個部分電阻率比值的上下限分別設(shè)為[1.20, 17.40, 1.20, 2.25]，[0.81, 13.60, 0.81, 1.75]. 由于大腦和頭皮部分的電阻率比值相等，為簡化計算，將DE算法的變量個數(shù)設(shè)定為D= 3,目標(biāo)函數(shù)為式(4)，算法終止條件為：最大迭代次數(shù)為130或個體目標(biāo)函數(shù)值f(ρ)小于10-6，如果在算法終止前沒有搜索到近似最優(yōu)解則認(rèn)為搜索失敗.

為了對算法圖像重構(gòu)質(zhì)量進行定量分析，引入誤差總和TE準(zhǔn)則，定義如下：

(8)

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

為確定DE基本策略DE/rand/1/bin中各個參數(shù)對圖像重構(gòu)的影響，對三個參數(shù)值的選取進行了仿真試驗. 種群大小NP設(shè)置為10到60，步長為5；縮放因子范圍F為0.1到1.5，步長0.1；交叉因數(shù)CR范圍為0.1到1.0，步長為0.1. 與其他進化算法類似，DE是通過隨機搜索來獲得最優(yōu)解的，每次DE的輸出都是不確定的. 因此需要進行多次重復(fù)試驗以獲得平均輸出結(jié)果.

為確定參數(shù)組合的重復(fù)試驗次數(shù)，隨機選擇6組參數(shù)組合進行300次重復(fù)試驗，計算前n(n為1 ～ 300之間的整數(shù))次誤差總和和求解目標(biāo)函數(shù)次數(shù)的平均值. 其結(jié)果如圖2、3，其中，用 (x1,x2,x3) 代表參數(shù)NP、F、CR的值分別為x1、x2、x3時的一組參數(shù)組合，求解目標(biāo)函數(shù)次數(shù)nf=NP× (g+1 )，g為算法收斂代數(shù).

圖2 平均誤差總和隨運行次數(shù)變化曲線Figure 2 Change of average total error with run times

圖3 平均求解目標(biāo)函數(shù)次數(shù)隨運行次數(shù)變化曲線Figure 3 Change of average number of objective function solved with run times

從圖2、3中可看出，當(dāng)每組參數(shù)組合的運行次數(shù)在120以下時，平均誤差總和和平均求解目標(biāo)函數(shù)次數(shù)曲線變動較為劇烈，而運行120次以上時，其曲線基本上為一條直線，說明當(dāng)實驗次數(shù)達(dá)到120以后，增加實驗次數(shù)，基本上不能提高圖像重建算法的性能. 因此，在以下的試驗中，將每組參數(shù)組合的重復(fù)實驗次數(shù)設(shè)置為120，取其搜索成功率，并計算在搜索成功情況下求解目標(biāo)函數(shù)次數(shù)以及誤差總和的平均值.

3.3 種群規(guī)模對重構(gòu)結(jié)果的影響

從算法的復(fù)雜度來看，種群規(guī)模NP越大，算法搜索范圍越大，得到全局最優(yōu)解的可能性越大，但所需的計算時間也越長. 因此種群規(guī)模的合理選取對提高算法的搜索效率具有重要的作用. 為確定NP對圖像重構(gòu)的影響，在三組F與CR的值分別設(shè)定為0.3、0.7，0.5、0.5，0.8、0.3的情況下，將NP由5增加到60，步長為5. 其重構(gòu)結(jié)果如圖4～6.

圖4 搜索成功率隨種群規(guī)模變化曲線Figure 4 Change of search success rate with population size

圖5 平均誤差總和隨種群規(guī)模變化曲線Figure 5 Change of average total error with population size

圖6 平均求解目標(biāo)函數(shù)次數(shù)隨種群規(guī)模變化曲線Figure 6 Change of average number of objective function solved with population size

由圖4可知，隨著NP的增大，算法搜索成功率越來越高.NP小于10時，算法的搜索成功率都小于80%，而當(dāng)NP在25以上時，算法的搜索成功率都達(dá)到了100%. 由圖5可知，NP增大時，誤差總和呈無規(guī)律波動，但所有平均誤差總和都處于2.5×10-3～ 2.8×10-3之間，解的精度基本不變. 但從圖6可看出，隨著NP的增大，平均求解目標(biāo)函數(shù)的次數(shù)不斷增大，收斂速度越來越慢. 因此，在給定最大迭代次數(shù)的情況下，NP在10 ～ 40之間時，在保證搜索成功率和圖像重構(gòu)精度的情況下，算法能夠獲得較高的搜索效率，性能較好.

3.4 縮放因子對重構(gòu)結(jié)果的影響

由式(4)可知，在一定程度上縮放因子F決定了種群進化的方向，對算法的局部搜索及全局搜索能力起到一定的調(diào)節(jié)作用.F較小時，種群多樣性較低，但能夠提高算法局部搜索能力同時加快算法的收斂；F較大時，有利于提高種群的多樣性，算法類似于隨機搜索，降低了算法的收斂速度.

為研究F對圖像重構(gòu)結(jié)果的影響，在三組NP與CR的值分別設(shè)定為20、0.4，30、0.5，40、0.9的情況下，將F由0.1增加到1.5，步長為0.1. 其重構(gòu)結(jié)果如圖7～9.

由圖7可知，當(dāng)F< 0.4時，在一定的迭代次數(shù)下，算法找到最優(yōu)解概率低于100%，但隨著F的增大，其搜索成功率逐漸增大. 當(dāng)F≥ 0.4時，其搜索成功率達(dá)到100%. 從圖8可知，隨著F的改變，平均誤差總和呈無規(guī)則變化，但都處于2.2×10-3～ 3.0×10-3之間，解的精度變化不大. 但從圖9可看出，當(dāng)F≥ 0.4時，隨著F的增大，雖然解的精度不一定降低，其求解目標(biāo)函數(shù)次數(shù)不斷增大，算法收斂速度逐漸減慢. 因此，在一定的迭代次數(shù)下，為確保算法在種群多樣性和收斂速度上取得良好的平衡，將F的初始合適區(qū)間取為0.4 ～ 1.0.

圖7 搜索成功率隨縮放因子變化曲線Figure 7 Change of search success rate with scaling factor

圖8 平均誤差總和隨縮放因子變化曲線Figure 8 Change of average total error with scaling factor

圖9 平均求解目標(biāo)函數(shù)次數(shù)隨縮放因子變化曲線Figure 9 Change of average number of objective function solved with scaling factor

3.5 交叉因數(shù)對重構(gòu)結(jié)果的影響

對交叉因數(shù)CR對重構(gòu)結(jié)果的試驗中，按照類似的方法，先固定NP和F，然后將CR的值從0.1增加到1.0，步長為0.1. 其部分實驗結(jié)果如圖10、11，其中所有參數(shù)組合設(shè)置的搜索成功率均為100%.

圖10 平均誤差總和隨交叉因數(shù)變化曲線Figure 10 Change of average total error with crossover factor

圖11 平均求解目標(biāo)函數(shù)次數(shù)隨交叉因數(shù)變化曲線Figure 11 Change of average number of objective function solved with crossover factor

從圖10、11可看出，隨著CR的增大，整體的誤差總和呈現(xiàn)下降的趨勢，且平均求解目標(biāo)函數(shù)次數(shù)不斷下降. 說明CR較小時，種群進化能力降低，收斂速度慢、所需計算時間較長；而CR較大時，種群多樣性增強，收斂速度快、所需計算時間短. 因此，為提高圖像重構(gòu)效率和精度，將CR取值設(shè)定在0.4 ～ 1.0區(qū)間上.

4 結(jié) 論

DE算法是一種易于種群進化的啟發(fā)式算法，逐漸應(yīng)用于EIT圖像重建中. 本實驗在二維同心圓模型上采用DE算法的基本策略DE/rand/1/bin對頭部電阻率進行重構(gòu)實驗，在給定目標(biāo)函數(shù)和終止條件的基礎(chǔ)上，對影響重建質(zhì)量和效率的三個主要控制參數(shù)進行了詳細(xì)的分析與研究，并給出了電阻抗重構(gòu)中DE算法控制參數(shù)合理選取范圍. 結(jié)果表明，基于DE算法的EIT重構(gòu)算法的性能對參數(shù)的選擇非常敏感. 而NP、F、CR分別設(shè)定在10～40、0.4～1.0、0.4～1.0區(qū)間上時，在保證搜索成功率的前提下，圖像重構(gòu)精度和重構(gòu)效率較高.

本文工作為EIT圖像重構(gòu)中DE算法控制參數(shù)的選取提供了可靠的定量分析數(shù)據(jù)，對EIT圖像重構(gòu)中DE算法的改進與參數(shù)優(yōu)化具有一定的指導(dǎo)意義. 下一步將對參數(shù)間的相互影響以及三維EIT圖像重構(gòu)中的參數(shù)選取問題展開研究. 此外在本文的基礎(chǔ)上，對DE算法的改進還有許多問題亟需解決.

[1] SILVERA-TAWIL D, RYE D, SOLEIMANI M, et al. Electrical impedance tomography for artificial sensitive robotic skin: a review[J].IEEESensorsJournal, 2015, 15(4): 2001-2016.

[2] FEITOSA A R S, RIBEIRO R R, BARBOSA V A F, et al. Reconstruction of electrical impedance tomography images using particle swarm optimization, genetic algorithms and non-blind search[C]//BiosignalsandBioroboticsConference. Salvador: IEEE, 2014:1-6.

[3] WAGENAAR J, ADLER A. Electrical impedance tomography in 3D using two electrode planes: Characterization and evaluation[J].PhysiologicalMeasurement, 2016, 37(6): 922-937.

[4] SCHULLCKE B, GONG B, MOELLER K. Steps towards 3D Electrical Impedance Tomography[C]//InternationalConferenceoftheIEEEEngineeringinMedicineandBiologySociety. Milan: IEEE, 2015: 5323-5326.

[5] MARTIN S, CHOI C T M. Nonlinear electrical impedance tomography reconstruction using artificial neural networks and particle swarm optimization[J].IEEETransactionsonMagnetics, 2016, 52(3): 1-4.

[6] FAN Q, YAN X. Self-adaptive differential evolution algorithm with discrete mutation control parameters[J].ExpertSystemswithApplications, 2015, 42 (3): 1551-1572.

[7] YIT K K, PARVATHY R. Differential-evolution control parameter optimization for unmanned aerial vehicle path planning[J].PlosOne, 2016, 11(3): 1-12.

[8] LEON M, XIONG N. Adapting differential evolution algorithms for continuous optimization via greedy adjustment of control parameters[J].JournalofArtificialIntelligence&SoftComputingResearch, 2016, 6(2): 103-118.

[10] HAMILTON S J, LASSAS M, SILTA-NEN S. A direct reconstruction method for anisotropic electrical impedance tomography[J].InverseProblems, 2014, 30(7): 1375-1378.

[11] 閆丹丹, 陳會. 基于電阻抗成像技術(shù)的腦病變組織檢測仿真研究[J]. 中國醫(yī)學(xué)影像技術(shù), 2014, 30(7): 1113-1116. YAN D D, CHEN H. Detection of brain anomaly based on electrical impedance imaging simulation research [J].ChineseJournalofMedicalImagingTechnology, 2014, 30(7): 1113-1116.

algorithmbasedondifferentialevolution

LIU Hai, YAN Dandan, JING Huimin

(College of Information Engineering, China Jiliang University, Hangzhou 310018, China)

2096-2835(2017)03-0345-07

10.3969/j.issn.2096-2835.2017.03.013

2017-06-15 《中國計量大學(xué)學(xué)報》網(wǎng)址zgjl.cbpt.cnki.net

國家自然科學(xué)基金資助項目(No.51107130).

TP39；R3

A