用差異演化-粒子群混合算法確定含水層參數(shù)

段國(guó)榮，劉元會(huì)

(長(zhǎng)安大學(xué) 理學(xué)院，陜西 西安 710064)

0 引 言

在中國(guó)工業(yè)快速發(fā)展的進(jìn)程中，水資源問(wèn)題逐漸增多，地下水是水資源的重要組成部分，地下水資源的評(píng)價(jià)及開(kāi)發(fā)利用就變得特別重要。含水層參數(shù)是進(jìn)行地下水資源評(píng)價(jià)和開(kāi)發(fā)利用的基礎(chǔ)數(shù)據(jù)。開(kāi)發(fā)利用地下水資源時(shí)，往往會(huì)遇到含水層具有直線供水邊界的情況。標(biāo)準(zhǔn)曲線配線法[1]，時(shí)間定律法(Law of Times)和拐點(diǎn)法(Method of Inflection Point)[2]，非線性最小二乘法[3-4]，線性回歸法[5]，Su shil K.S.方法[6-7]，直線圖解法[8]和割離井法[9-10]等是分析直線供水邊界含水層抽水試驗(yàn)數(shù)據(jù)的傳統(tǒng)方法，但此類方法在實(shí)際應(yīng)用中受到一定的限制。近年來(lái)，人們廣泛使用智能優(yōu)化算法來(lái)分析抽水試驗(yàn)數(shù)據(jù)，求解含水層參數(shù)，例如粒子群優(yōu)化算法[11]，混沌人工魚(yú)群混合算法[12]，遺傳算法[13]，改進(jìn)粒子群算法[14]，單純形探索法[15]，單純形差分進(jìn)化混合算法[16]，混沌粒子群優(yōu)化算法[17]和單純形-粒子群混合算法[18]等方法。郭建青，李彥等將粒子群優(yōu)化算法[11]應(yīng)用于確定含水層參數(shù)，但該算法對(duì)粒子數(shù)目的大小依賴性很強(qiáng)。袁帆，劉元會(huì)等利用單純形-粒子群混合算法[18]計(jì)算含水層參數(shù)的問(wèn)題，但計(jì)算結(jié)果精度不高，且粒子數(shù)目的大小和待估參數(shù)取值范圍對(duì)算法的收斂性有較大的影響。粒子群優(yōu)化算法[11]操作簡(jiǎn)單且易實(shí)現(xiàn)，但種群多樣性的缺失使其易陷入局部最優(yōu)。因此，本研究對(duì)粒子群優(yōu)化算法(PSO算法)進(jìn)行改進(jìn)，提供一種差異演化-粒子群混合算法(DE-PSO算法)。將差異演化-粒子群混合算法用于分析直線供水邊界條件下含水層的抽水試驗(yàn)數(shù)據(jù)，求解含水層參數(shù)，計(jì)算觀測(cè)孔與虛擬映射井之間的距離。本研究重點(diǎn)討論了該算法的可靠性和有效性，并就待估參數(shù)取值范圍對(duì)算法收斂性的影響進(jìn)行數(shù)值實(shí)驗(yàn)分析。

1 優(yōu)化算法簡(jiǎn)介

1.1 DE-PSO算法的思想

PSO算法在迭代后期由于種群多樣性的下降存在收斂速度慢，精度差，且易陷入局部最優(yōu)等問(wèn)題，利用DE算法[19-21]中的3種遺傳操作，即變異操作，交叉操作和選擇操作，變異粒子中各個(gè)體的歷史最優(yōu)位置，不僅能使種群保持一定的多樣性，進(jìn)而避免出現(xiàn)早熟或停滯的現(xiàn)象，而且能保持PSO算法迭代前期收斂速度快的優(yōu)點(diǎn)。本研究引入聚集度因子[22]來(lái)確定粒子個(gè)體歷史最優(yōu)位置是否需要進(jìn)行變異

(1)

1.2 DE-PSO算法的步驟

步驟1：初始化算法的相關(guān)參數(shù)；確定粒子個(gè)數(shù)N，搜索空間維度D，最大迭代次數(shù)gen，最優(yōu)累計(jì)次數(shù)K[23]，各粒子的初始位置及初始速度、搜索空間的上下限、縮放因子F，交叉概率CR，初始化循環(huán)迭代次數(shù)t和最優(yōu)累計(jì)次數(shù)k，聚集度因子Fa，收斂精度e1以及計(jì)算相似度e2[23].

步驟2：計(jì)算初始種群中每個(gè)粒子的適應(yīng)度值，保存粒子局部和全局最優(yōu)位置，最優(yōu)值。

步驟3：若當(dāng)前迭代次數(shù)小于最大迭代次數(shù)gen且當(dāng)前的最優(yōu)累計(jì)次數(shù)k小于最優(yōu)累計(jì)次數(shù)K，則進(jìn)入循環(huán)，否則輸出最優(yōu)結(jié)果。

步驟4：更新每個(gè)粒子的速度和位置，且進(jìn)行粒子個(gè)體最優(yōu)更新及群體最優(yōu)更新。

步驟5：計(jì)算粒子的聚集度因子Fa，若Fa大于所給定的閾值，則進(jìn)入步驟6，否則返回步驟4.

步驟6：利用DE算法中的變異、交叉、選擇等遺傳操作，變異粒子中各個(gè)體的歷史最優(yōu)位置，并更新粒子歷史最優(yōu)位置。

步驟7：如果全局最優(yōu)解小于e1且全局最優(yōu)解的誤差小于相似度e2，則k=k+1，否則返回步驟3繼續(xù)執(zhí)行。

步驟8：輸出結(jié)果。

2 目標(biāo)函數(shù)與控制條件

2.1 降深表達(dá)式

由地下水動(dòng)力學(xué)[24-25]中的迭加原理可知，含水層中任一位置的降深可以由2部分迭加得到，即在直線供水邊界條件下，虛擬映射注水井注水和抽水實(shí)井抽水分別在該位置產(chǎn)生降深的疊加。虛擬映射注水井、抽水實(shí)井和直線供水邊界的位置關(guān)系[5]，如圖1所示。

圖1 直線供水邊界條件下抽水實(shí)井和虛擬映射注水井的位置關(guān)系Fig.1 Location relationship between the pumping well and the virtual mapping injection well in linear pervious boundary

該情況下，含水層中任一位置的降深表達(dá)式[26]如下

s=s1+s2=[w(u1)-w(u2)]Q/4πT

(2)

式中Q為抽水流量，L3/T1；T為含水層的導(dǎo)水系數(shù)，L2/T2；s1為直線供水邊界條件下實(shí)井抽水引起的降深，L；s2為直線供水邊界條件下虛擬映射注水井引起的降深，L；w(u1)和w(u2)分別為相應(yīng)的抽水實(shí)井和虛擬映射注水井的井函數(shù)，無(wú)量綱，其表達(dá)式分別為

(3)

其中無(wú)量綱時(shí)間參量u1和u2的表達(dá)式分別為

(4)

式中t為抽水持續(xù)時(shí)間,T；S為儲(chǔ)水系數(shù)，無(wú)量綱；r為抽水實(shí)井與觀測(cè)井的距離，L；r′為虛擬映射井與觀測(cè)井的距離，L.

③有法不依、執(zhí)法不嚴(yán)的問(wèn)題比較突出。由于管理體制不順，違反《塔里木河水資源管理?xiàng)l例》、自治區(qū)批準(zhǔn)的水量分配方案、水量調(diào)度管理辦法，不執(zhí)行水調(diào)指令，搶占、擠占生態(tài)水，不按塔里木河規(guī)劃確定的輸水目標(biāo)向塔里木河輸水的現(xiàn)象時(shí)有發(fā)生。塔里木河干流水權(quán)及生態(tài)用水水權(quán)得不到依法保護(hù)，塔管局由于現(xiàn)行的管理權(quán)限所限，難以依法進(jìn)行處罰。

由式(3)可知，在計(jì)算水位降深值時(shí)，需要對(duì)井函數(shù)進(jìn)行廣義積分求解。采用R.Srivastava給出的近似表達(dá)式[27]進(jìn)行計(jì)算

W(u)=-Inu+a0+a1u+a2u2+a3u3+a4u4+a5u5,u≤1

(5)

(6)

式中a0=-0.577 72，a1=0.999 99，a2=-0.249 1，a3=-0.055 19，a4=-0.009 76，a5=0.001 08，b0=0.267 77，b1=8.634 76，b2=18.059 02，b3=3.958 50，c0=3.958 50，c1=21.099 65，c2=25.632 96，c3=9.573 32.

2.2 目標(biāo)函數(shù)的構(gòu)成

采用DE-PSO算法求解含水層參數(shù)時(shí)，待估參數(shù)向量θ=(T,S,r′)應(yīng)滿足降深觀測(cè)值與計(jì)算值的方差達(dá)到最小，構(gòu)造如下目標(biāo)函數(shù)

(7)

2.3 算法參數(shù)控制

聚集度因子Fa∈[0,2]，根據(jù)多次試驗(yàn)結(jié)果，本研究選取Fa的閾值Fa=1.3，交叉概率CR∈[0,1]，縮放因子F∈[0,2]，粒子數(shù)目維度N=20，最大迭代次數(shù)gen=200，根據(jù)文獻(xiàn)[28],選取含水層導(dǎo)水系數(shù)T∈[2.5,3.5] m2/min，儲(chǔ)水系數(shù)S∈[0.050,0.070]，觀測(cè)井到虛擬映射井的距離r′∈[100,130] m，取e1=5×10-5，e2=10-9，K=10.

3 算例

3.1 不同方法計(jì)算結(jié)果的比較

表2中給出了DE-PSO算法和相關(guān)文獻(xiàn)中的其他方法對(duì)含水層參數(shù)的計(jì)算結(jié)果。

從表2可知，DE-PSO算法估計(jì)含水層參數(shù)值及觀測(cè)井與虛擬映射井之間的距離與其他方法接近，且由該算法求得的φ(θ)精度高于其他方法，因此，DE-PSO算法所求結(jié)果是可靠的。

3.2 算法的有效性

圖2表示的是DE-PSO算法的程序運(yùn)行100次目標(biāo)函數(shù)值的分布情況。從圖2能夠看出，在DE-PSO算法的程序運(yùn)行100次的過(guò)程中，有91次成功，其收斂率達(dá)到91%.由此可知，DE-PSO算法的收斂率較高，具有高效性。因此，利用DE-PSO算法能有效地估計(jì)含水層參數(shù)。

表1 原始抽水試驗(yàn)數(shù)據(jù)

表2 不同方法的計(jì)算結(jié)果

圖2 DE-PSO算法運(yùn)行100次目標(biāo)函數(shù)值分布Fig.2 Distribution of the objective function value after 100 runs of DE-PSO algorithm

3.3 算法的穩(wěn)定性

表3是在其他條件不變，待估參數(shù)取值范圍的上限擴(kuò)大到原來(lái)的2，4，8，10，12和14倍的情況下，利用DE-PSO算法求得含水層參數(shù)及觀測(cè)井與虛擬映射井之間距離的估計(jì)值。由表中信息可知，隨著待估參數(shù)范圍的增大，DE-PSO算法具有很好的收斂性，當(dāng)待估參數(shù)取值范圍的上限擴(kuò)大到原來(lái)的14倍時(shí)，收斂率高達(dá)97%，且得到目標(biāo)函數(shù)值基本一致，即DE-PSO算法對(duì)初值選取的敏感性低，尋優(yōu)能力強(qiáng)，穩(wěn)定性好。

表3 不同參數(shù)取值范圍下的計(jì)算結(jié)果

4 結(jié) 論

1)DE-PSO算法得到的目標(biāo)函數(shù)值φ(θ)精度高于其他方法。

2)估計(jì)含水層參數(shù)時(shí)具有較好的收斂性，具有較高的精度。

3)待估參數(shù)初始取值范圍的放寬對(duì)DE-PSO算法的收斂性影響較小，即DE-PSO算法對(duì)初始取值范圍選取的敏感性低，尋優(yōu)能力強(qiáng)，穩(wěn)定性好。因此，DE-PSO算法是分析抽水試驗(yàn)數(shù)據(jù)，求解含水層參數(shù)和計(jì)算觀測(cè)孔與虛擬映射井之間距離的有效方法。