粒子群優(yōu)化算法在確定越流含水層參數(shù)中的應(yīng)用

張鵠志,郭建青

(長(zhǎng)安大學(xué)環(huán)境科學(xué)與工程學(xué)院,陜西 西安 710054)

粒子群優(yōu)化(particle swarm optimization,PSO)算法[1-2]是由Eberhart R C和Kennedy J于1995年提出的一種基于迭代的智能隨機(jī)優(yōu)化算法,具有易于理解和實(shí)現(xiàn)、全局搜索能力強(qiáng)[3]和計(jì)算結(jié)果精度不易受人為因素影響等優(yōu)點(diǎn)。由于PSO算法具有這些優(yōu)點(diǎn),其在水利工程規(guī)劃設(shè)計(jì)領(lǐng)域中已經(jīng)得到了一定的應(yīng)用。例如,國(guó)內(nèi)學(xué)者將PSO算法應(yīng)用于水庫(kù)(群)優(yōu)化調(diào)度中,為求解水庫(kù)調(diào)度問題提供了一種新的有效解決途徑[4-5];還有的學(xué)者將PSO算法引入大壩安全監(jiān)控領(lǐng)域,并結(jié)合多元回歸統(tǒng)計(jì)模型建立基于粒子群算法的混凝土壩變形預(yù)報(bào)模型,并通過(guò)優(yōu)化迭代計(jì)算確定壩體變形統(tǒng)計(jì)模型中各回歸系數(shù)[6]。在開發(fā)利用地下水資源過(guò)程中,利用不同的方法分析抽水試驗(yàn)數(shù)據(jù),是確定含水層參數(shù)的主要途徑之一。標(biāo)準(zhǔn)曲線配線法是分析抽水試驗(yàn)數(shù)據(jù)的常用方法之一,但是,其在應(yīng)用過(guò)程中存在著人為隨意性,對(duì)計(jì)算結(jié)果的精度影響較大。為了彌補(bǔ)這個(gè)缺陷,人們已將智能優(yōu)化算法應(yīng)用于分析抽水試驗(yàn)數(shù)據(jù),確定滿足泰斯假設(shè)條件下的含水層參數(shù)[7]。但是,在利用PSO算法分析越流系統(tǒng)含水層中的抽水試驗(yàn)數(shù)據(jù)并確定含水層參數(shù)方面,以及PSO算法中的3個(gè)控制參數(shù)對(duì)算法收斂性和計(jì)算結(jié)果精度的影響方面開展的研究工作較少。鑒于此,筆者將PSO算法應(yīng)用于確定第1類越流系統(tǒng)的含水層參數(shù),并通過(guò)數(shù)值試驗(yàn),初步討論算法控制參數(shù)c1,c2和 ω對(duì)算法收斂性的影響。

1 PSO算法

PSO算法假想一群粒子的每個(gè)粒子都在搜索空間中以一定速度飛行,并根據(jù)本身以及同伴的飛行經(jīng)驗(yàn)動(dòng)態(tài)調(diào)整自身速度,從而形成群體尋優(yōu)的機(jī)制,最終搜索到問題的最優(yōu)解[1-2]。人們用粒子的位置表示待優(yōu)化問題的解,將每個(gè)粒子位置的優(yōu)劣程度稱為它的適應(yīng)度,由待優(yōu)化問題的目標(biāo)函數(shù)值來(lái)定量地表征。每個(gè)粒子由1個(gè)速度矢量決定其飛行方向和速率大小。設(shè)在1個(gè)d維的目標(biāo)搜索空間中,由m個(gè)粒子組成1個(gè)群體,其中,在第 t次迭代時(shí)粒子i的位置表示為Xi(t)=(xi1(t),xi2(t),…,xid(t)),相應(yīng)的飛行速度表示為 V i(t)=(vi1(t),vi2(t),…,vid(t))。開始運(yùn)算時(shí),首先隨機(jī)初始化這m個(gè)粒子的位置和速度,然后通過(guò)迭代尋找最優(yōu)解。在每次迭代過(guò)程中,各粒子通過(guò)跟蹤2個(gè)極值來(lái)更新自己的速度和位置:1個(gè)是該粒子本身迄今搜索到的最優(yōu)解,稱為個(gè)體極值 pbest,表示為P i(t)=(pi1(t),pi2(t),…,pid(t));另1個(gè)是整個(gè)粒子群迄今找到的最優(yōu)解,稱為全局極值gbest,表示為Pg(t)=(pg1(t),pg2(t),…,pgd(t))。將每個(gè)粒子更新個(gè)體極值的過(guò)程稱為歷史尋優(yōu),將各不同粒子通過(guò)相互比較來(lái)更新全局極值的過(guò)程稱為全局尋優(yōu)。在第t+1次迭代計(jì)算時(shí),粒子i根據(jù)下列規(guī)則來(lái)更新自己的速度和位置[3]:

式中:ω為慣性權(quán)重系數(shù);c1和c2為學(xué)習(xí)因子,文獻(xiàn)[3]中建議 c1=c2=2.0;rand1和 rand2分別為 2個(gè)在[0,1]區(qū)間內(nèi)變化的隨機(jī)數(shù);下標(biāo)k表示第k維。

式(1)右邊第1項(xiàng)為個(gè)體粒子飛行過(guò)程中的“慣性”部分,反映了粒子的運(yùn)動(dòng)“習(xí)慣”,代表粒子有維持自己先前速度的趨勢(shì);第2項(xiàng)為“認(rèn)知”部分,反映了粒子對(duì)自身歷史經(jīng)驗(yàn)的記憶,代表粒子有向自身個(gè)體極值逼近的趨勢(shì);第3項(xiàng)為“社會(huì)”部分,反映了粒子間協(xié)同合作與知識(shí)共享的群體歷史經(jīng)驗(yàn),代表粒子有向全局極值逼近的趨勢(shì)[3]。另外,粒子每一維的速度vik都應(yīng)被一個(gè)最大速度所限制。而為粒子位置設(shè)定閾值,從而保證各粒子在限定的搜索空間內(nèi)飛行,也是有待進(jìn)一步研究的問題。

2 PSO算法在確定第1類越流系統(tǒng)含水層參數(shù)中的應(yīng)用

2.1 基本公式

本文討論的第1類越流系統(tǒng)的模型與文獻(xiàn)[8]中相同,抽水過(guò)程中的水頭降深公式為

式中:s為水頭降深,[L];Q為抽水流量,[L3?T-1];T為導(dǎo)水系數(shù),[L2?T-1];r為距抽水主井的距離,[L];1/B為越流補(bǔ)給因子,[L-1]。

在應(yīng)用PSO算法確定含水層參數(shù)時(shí),要求式(4)表達(dá)的目標(biāo)函數(shù)值達(dá)到最小,即

式中:Xi(t)為第t次迭代時(shí)粒子i的位置,它以各含水層參數(shù)作為其各維上的坐標(biāo),在第1類越流系統(tǒng)中有3個(gè)待求的含水層參數(shù),分別是導(dǎo)水系數(shù)T、彈性釋水系數(shù)S和越流補(bǔ)給因子1/B,分別視為xi1(t),xi2(t),xi3(t),這樣粒子的位置坐標(biāo)就是X=(T,S,1/B);j為抽水試驗(yàn)過(guò)程中水頭降深觀測(cè)數(shù)據(jù)的序列號(hào),j=1,2,…,N;soj為實(shí)際觀測(cè)到的第j個(gè)降深值,[L];scj為利用式(3)計(jì)算的與 soj相同時(shí)刻的水頭降深值,[L]。

式(4)的意義為選取適當(dāng)?shù)牧Ｗ游恢?X=(T,S,1/B),使降深的計(jì)算值與其觀測(cè)值之間的殘差平方和的均值達(dá)到最小。此時(shí)的 T,S,1/B即為所求的最優(yōu)解。

2.2 算例

根據(jù)上述基本公式和PSO算法的流程,利用C語(yǔ)言編寫了相應(yīng)的計(jì)算程序,在運(yùn)行開始時(shí),輸入數(shù)據(jù)為相應(yīng)不同抽水時(shí)間的水頭降深值,輸出結(jié)果為含水層參數(shù)T,S,1/B的最優(yōu)值。

表1中給出了抽水流量Q=17.0m3/min時(shí),距抽水主井距離r=12.19 m的觀察孔中觀測(cè)到的降深隨時(shí)間的變化數(shù)據(jù)[10]。

表1 觀測(cè)數(shù)據(jù)

根據(jù)多次試算的經(jīng)驗(yàn),采用的粒子個(gè)數(shù)為50個(gè);收斂標(biāo)準(zhǔn)取目標(biāo)函數(shù)值小于0.00060842。權(quán)重系數(shù)按式(4)進(jìn)行線性衰減:式中:ωmax為最大慣性權(quán)重系數(shù),取 0.9;ωmin為最小慣性權(quán)重系數(shù),取0.1;t′為當(dāng)前迭代次數(shù);n為最大迭代次數(shù),取3000。

取學(xué)習(xí)因子c1=c2=2.0,經(jīng)過(guò)平均700次迭代得到計(jì)算結(jié)果,見表2。表2也給出了傳統(tǒng)配線法的計(jì)算結(jié)果。

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

2.3 可靠性檢驗(yàn)及與傳統(tǒng)配線法效果比較

分別利用由PSO算法和配線法求出的含水層參數(shù)來(lái)計(jì)算降深～時(shí)間曲線,并與實(shí)際的降深～時(shí)間曲線進(jìn)行比較,如圖1所示。由圖1可以看出,PSO算法所求結(jié)果對(duì)觀測(cè)數(shù)據(jù)的擬合程度與傳統(tǒng)配線法很接近,是可靠的。但配線法具有較大的主觀性,不同人員所得結(jié)果的精度可能不同,而PSO算法則可以穩(wěn)定地得到精度較高的結(jié)果。

圖1 PSO算法與配線法效果比較

3 PSO算法控制參數(shù)的作用和影響

由前文可以看出,有5個(gè)參數(shù)影響著PSO算法的運(yùn)算過(guò)程和結(jié)果,分別是:粒子個(gè)數(shù)、待求參數(shù)初值的取值范圍、收斂標(biāo)準(zhǔn)、c1和c2以及 ω。擬通過(guò)數(shù)值試驗(yàn)結(jié)果,分析 c1,c2和 ω對(duì)算法收斂性的影響。

3.1 學(xué)習(xí)因子

3.1.1 定性分析

以 c1為例,從式(1)和式(2)不難看出,c1rand 1[pik(t)-xik(t)]的作用是:在粒子更新位置的過(guò)程中,使粒子以“一定程度”向它的個(gè)體極值靠攏。而這個(gè)“一定程度”是由 c1rand1確定的。比如xik(t)+1.0rand1[pik(t)-xik(t)](即 c1=1.0)就表示粒子向其個(gè)體極值隨機(jī)地移動(dòng)某一距離;而xik(t)+2.0rand1[pik(t)-xik(t)](即 c1=2.0)表示粒子不但可以飛向個(gè)體極值,而且也可以飛過(guò)個(gè)體極值位置到達(dá)另一側(cè)(但不會(huì)超出原來(lái)相差的距離),如圖2所示。因此,當(dāng)c1=1.0時(shí),粒子直接以其個(gè)體極值為目標(biāo),靠攏的趨勢(shì)是最強(qiáng)的;當(dāng) c1=2.0時(shí),新位置在個(gè)體極值兩側(cè)的分布幾率是最均衡的。

圖2 粒子飛行示意圖

同理,c2rand2確定了粒子向全局極值移動(dòng)的程度:當(dāng)c2=1.0時(shí),粒子向全局極值靠攏的趨勢(shì)最強(qiáng);當(dāng)c2=2.0時(shí),粒子趨于在全局極值周圍幾率均衡地分布。

3.1.2 數(shù)值試驗(yàn)結(jié)果分析

文中主要對(duì)最具代表性的 ci=1.0和 ci=2.0(i=1,2)情況進(jìn)行試算。另由于不同目標(biāo)函數(shù)值對(duì)應(yīng)著不同的粒子位置,因此用目標(biāo)函數(shù)值的變化來(lái)反映粒子位置的變化情況。

筆者在其他參數(shù)取值相同(c2=1.0,ωmax=0.8,ωmin=0.6,粒子個(gè)數(shù)取 5個(gè),收斂標(biāo)準(zhǔn)取0.00060842)的條件下,含水層參數(shù)取值范圍為0.93 m2/min≤T≤3.72m2/min,0.002≤S≤0.005,0.0016 m-1≤1/B ≤0.0033 m-1時(shí),對(duì) c1=2.0和c1=1.0分別進(jìn)行3次試算,輸出每個(gè)粒子目標(biāo)函數(shù)值隨迭代次數(shù)增加的變化(篇幅所限未列出具體結(jié)果)。通過(guò)試算可以觀察到2種現(xiàn)象:①c1影響粒子的活躍程度,進(jìn)而影響粒子的搜索能力:當(dāng) c1=2.0時(shí),目標(biāo)函數(shù)值在迭代次數(shù)達(dá)到2000次以上時(shí)仍然有所變化,即粒子仍然移動(dòng),保持著搜索能力;當(dāng)c1=1.0時(shí),目標(biāo)函數(shù)值在迭代次數(shù)達(dá)到1000次以上時(shí)就已經(jīng)不再變化,即粒子基本失去搜索能力。②當(dāng)c1=2.0時(shí),3次試算目標(biāo)函數(shù)值分別穩(wěn)定于0.00061086,0.00061082和0.00061098左右;c1=1.0時(shí),3次試算目標(biāo)函數(shù)值分別穩(wěn)定于0.00086816,0.00071975和0.00101040,均大于c1=2.0時(shí)的結(jié)果。顯然 c1=2.0時(shí)的最終運(yùn)算結(jié)果具有更高的精度。

產(chǎn)生上述現(xiàn)象的原因是,當(dāng) c1=1.0時(shí),粒子在每次迭代中都強(qiáng)烈地趨向于它的某個(gè)個(gè)體極值,因此,在經(jīng)過(guò)較少次的迭代后粒子就會(huì)到達(dá)某個(gè)個(gè)體極值附近。這時(shí),負(fù)責(zé)歷史尋優(yōu)的c2rand1[pik(t)-xik(t)]項(xiàng)就會(huì)趨于零(因?yàn)檫@時(shí)pik(t)-xik(t)趨于零),從而使歷史尋優(yōu)機(jī)制形同虛設(shè),粒子便停滯在這個(gè)個(gè)體極值而不再隨迭代而移動(dòng)(實(shí)際上,試算中c2=1.0的取法加劇了這種作用,它使各粒子的個(gè)體極值迅速趨向于某一局部極值)。粒子這種過(guò)早的停滯使其不能再搜索到更優(yōu)的值,因此其最終結(jié)果的精度普遍較差。當(dāng)c1=2.0時(shí),粒子在每次迭代中都能幾率均衡地在個(gè)體極值兩側(cè)進(jìn)行歷史尋優(yōu),所以不易停滯,從而保持了更強(qiáng)的歷史尋優(yōu)能力,可以找到更優(yōu)的值。

對(duì)于c2的取值,在其他參數(shù)取值相同(c1=1.0,ωmax=0.9,ωmin=0.7,粒子個(gè)數(shù)取 5 個(gè),收斂標(biāo)準(zhǔn)取0.00060842)的條件下,含水層參數(shù)取值范圍為:0.93 m2/min≤T≤3.72m2/min,0.002≤S≤0.005,0.0016m-1≤1/B ≤0.0033m-1時(shí) ,對(duì) c2=2.0和c2=1.0分別進(jìn)行3次運(yùn)算,其觀察到的現(xiàn)象以及原因分析都與 c1時(shí)非常相似,不同的是c2直接對(duì)全局尋優(yōu)而不是歷史尋優(yōu)起作用。

實(shí)際上,c1和c2的作用是相互交織的,歷史尋優(yōu)的效果和全局尋優(yōu)的效果是相互影響的。

3.2 慣性權(quán)重系數(shù)

由式(1)可以看出,ω決定了對(duì)前一速度值繼承的程度。它影響粒子的搜索能力并影響搜索的精細(xì)程度。

筆者在其他參數(shù)取值相同(c1=c2=1.0,粒子個(gè)數(shù)取5個(gè),收斂標(biāo)準(zhǔn)取0.00060842)的條件下,含水層參數(shù)取值范圍為 0.93 m2/min≤T≤3.72m2/min,0.002≤S≤0.005,0.0016 m-1≤1/B≤0.0033 m-1時(shí),對(duì) 0.5≤ω≤1.0和0≤ω≤0.5分別進(jìn)行3次試算,輸出每個(gè)粒子目標(biāo)函數(shù)值隨迭代次數(shù)增加的變化。通過(guò)試算觀察到這樣的現(xiàn)象:當(dāng) ω較大(0.5≤ω≤1.0)、迭代次數(shù)達(dá)到1000次以上時(shí)各粒子的目標(biāo)函數(shù)值仍有較大變化,即粒子仍保持搜索能力,3次運(yùn)算的最終目標(biāo)函數(shù)值分別穩(wěn)定于0.00062041,0.00061140,0.00061102;當(dāng) ω較小(0≤ω≤0.5)、迭代次數(shù)不到1000次時(shí)各粒子已經(jīng)基本停滯,3次運(yùn)算的最終目標(biāo)函數(shù)值分別穩(wěn)定于0.00110930,0.00066110,0.00064549。產(chǎn)生上述現(xiàn)象的原因是,ω值越大,前一速度對(duì)當(dāng)前的速度貢獻(xiàn)越大,當(dāng)粒子聚集到個(gè)體極值和全局極值附近時(shí),式(1)中c1rand1[pik(t)-xik(t)]+c2rand2[pgk(t)-xik(t)]變得很小,這樣 ωvik(t)項(xiàng)就成為新速度的主要部分;如果 ω很小,那么新的速度就會(huì)迅速變小,導(dǎo)致粒子活動(dòng)范圍越來(lái)越小,全局搜索能力降低。但有研究認(rèn)為,ω較小時(shí)可增強(qiáng)局部搜索能力[3]。

4 結(jié) 語(yǔ)

文中將PSO算法應(yīng)用于分析第1類越流系統(tǒng)含水層中的抽水試驗(yàn)數(shù)據(jù),并確定含水層的參數(shù)。通過(guò)對(duì)數(shù)值試驗(yàn)結(jié)果的初步分析討論,可以得出:①PSO算法原理簡(jiǎn)單,易于編程運(yùn)算,計(jì)算結(jié)果精度不受人為影響;②PSO算法的計(jì)算結(jié)果精度高,可靠性強(qiáng);③學(xué)習(xí)因子c1,c2和慣性權(quán)重系數(shù) ω綜合影響著粒子的活躍程度,進(jìn)而影響粒子群的搜索能力;④建議算法控制參數(shù)取值為 c1=c2=2.0,ω由0.9隨迭代次數(shù)線性遞減至0.4。總之,PSO算法不失為一種分析抽水試驗(yàn)資料、確定第1類越流系統(tǒng)含水層參數(shù)的有效方法。

[1]KENNEDY J,EBERHART R C.Particle swarm optimization[C]//Proceedings of IEEE International Conference on Neural Networks.Piscataway:IEBE Service Center,1995(4):1942-1948.

[2]DOROGO M,MANIEZZO V,COLORI A.Ant system:optimization by a colony of cooperating agents[J].IEEE Trans on System,Man,and Cybernetics,1996,26(1):28-41.

[3]范娜,云慶夏.粒子群優(yōu)化算法及其應(yīng)用[J].信息技術(shù),2006(1):53-56.

[4]鄧顯羽,彭勇,葉碎高,等.粒子群算法在水庫(kù)(群)優(yōu)化調(diào)度研究中的應(yīng)用綜述[J].水利水電科技進(jìn)展,2010,30(5):90-94.

[5]向波,紀(jì)昌明,羅慶松.免疫粒子群算法及其在水庫(kù)優(yōu)化調(diào)度中的應(yīng)用[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,2008,36(2):198-202.

[6]王偉,王連慶.基于粒子群仿生算法的混凝土壩變形預(yù)報(bào)模型[J].水利水電科技進(jìn)展,2008,28(4):11-14.

[7]張娟娟.智能優(yōu)化算法在含水層參數(shù)反演中的應(yīng)用[D].西安:長(zhǎng)安大學(xué),2006.

[8]陳崇希,林敏.地下水動(dòng)力學(xué)[M].武漢:中國(guó)地質(zhì)大學(xué)出版社,1999:70-122.

[9]楊天行,傅澤周,劉金山,等.地下水向井的非穩(wěn)定流動(dòng)的原理與計(jì)算方法[M].北京:地質(zhì)出版社,1980:195-199.

[10]U.S.Department of the Interior(USDI).Ground watermAnual[M].Washington,D.C.:Bureau of Reclamation,1981:121-127.