佘緯,安智,田莎莎,張相來(lái)
(1 中南民族大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)學(xué)院,武漢430074;2 中南民族大學(xué) 計(jì)算機(jī)科學(xué)學(xué)院,武漢430074)
徑流是陸地上水量最直觀的反映,為了能夠更好地利用和保護(hù)水資源,我們不僅要了解徑流的長(zhǎng)期變化規(guī)律,還要做好徑流的預(yù)測(cè)工作,來(lái)指導(dǎo)人們的生產(chǎn)和生活[1,2].由于地理位置偏僻或者自然災(zāi)害,有時(shí)會(huì)出現(xiàn)水文站水文數(shù)據(jù)缺失的現(xiàn)象,這就需要我們能夠根據(jù)鄰近站點(diǎn)的水文數(shù)據(jù)來(lái)預(yù)測(cè)出缺失水文站的水文數(shù)據(jù).由于水文環(huán)境和水文現(xiàn)象的復(fù)雜性,單變量的水文分析已經(jīng)不能滿足工程設(shè)計(jì)和人們生產(chǎn)生活的需要,在進(jìn)行水文分析時(shí),往往會(huì)涉及到多個(gè)參數(shù),而Copula函數(shù)的優(yōu)點(diǎn)正是可以描述變量間的相關(guān)性,而且可以通過(guò)連接每個(gè)變量的邊緣分布,得到這些變量的聯(lián)合分布函數(shù).因此,Copula函數(shù)得到了水文學(xué)家的青睞,越來(lái)越多的用于水文學(xué)的相關(guān)分析和預(yù)測(cè)問(wèn)題中.近年來(lái)很多學(xué)者將Copula函數(shù)應(yīng)用于多變量水文頻率分析中[3,4].2007年,《Journal of Hydrologic Engineering》雜志對(duì)Copula函數(shù)在水文學(xué)中的應(yīng)用進(jìn)行了??榻B[5-7].SONG利用三變量Plackett copula函數(shù)模擬了干旱持續(xù)時(shí)間、嚴(yán)重程度和到達(dá)時(shí)間的聯(lián)合概率分布,結(jié)果表明Plackett copula能夠產(chǎn)生相關(guān)干旱變量的雙變量和三變量概率分布[8].CHEN采用多變量Copula方法對(duì)干旱期進(jìn)行量化,以中國(guó)漢江上游流域?yàn)槔?,研究了各干旱狀態(tài)下的依賴結(jié)構(gòu),計(jì)算并分析了各干旱狀態(tài)下的干旱概率和恢復(fù)周期[9].REDDY研究了兩種元啟發(fā)式算法在Copula模型參數(shù)估計(jì)和干旱嚴(yán)重持續(xù)頻率曲線推導(dǎo)中的應(yīng)用[10].MASUD通過(guò)Copula函數(shù)建立二維分布,對(duì)加拿大Saskatchewan流域的干旱風(fēng)險(xiǎn)進(jìn)行了分析[11].YUE利用混合Gumbel模型建立了相關(guān)洪峰和洪量的聯(lián)合概率分布,以及相關(guān)洪峰和洪峰持續(xù)時(shí)間的聯(lián)合概率分布[12].REQUENA利用二元Copula模型得到洪峰與洪量的二元聯(lián)合分布,模擬了大壩溢流的回歸周期,為大壩設(shè)計(jì)風(fēng)險(xiǎn)提供了參考[13].JEONG研究了氣候變化對(duì)加拿大東北部21個(gè)流域春季(3 -6月)洪峰、洪量和洪峰持續(xù)時(shí)間的影響,對(duì)各洪水特征進(jìn)行常規(guī)的單變量頻率分析,對(duì)相互關(guān)聯(lián)的洪水特征對(duì)(峰量、峰時(shí)、量時(shí))進(jìn)行基于Copula的雙變量頻率分析[14].JEONG將Copula模型與PMC模型相結(jié)合,提出了一種可行的間歇月流時(shí)間序列模擬方法[15].
單一的Copula函數(shù)只能描述變量間某一方面的相關(guān)性,為了全面的描述變量之間的相關(guān)性,本文構(gòu)建了混合Copula函數(shù)來(lái)預(yù)測(cè)鄰近水文站點(diǎn)的年徑流量.構(gòu)建混合Copula函數(shù)的關(guān)鍵是參數(shù)的尋優(yōu)問(wèn)題,為此深入研究了通用性強(qiáng)、模型簡(jiǎn)單,搜索速度快的布谷鳥搜索算法.布谷鳥搜索算法是YANG等通過(guò)模擬布谷鳥的產(chǎn)卵習(xí)性而提出的一種與萊維飛行結(jié)合的智能搜索算法[16].該算法在算法運(yùn)行前期搜索速度較快,但是由于其完全依賴于隨機(jī)游走機(jī)制,所以算法的局部搜索精度不高,容易陷入早熟收斂,收斂速度較慢.為了改進(jìn)布谷鳥搜索算法的性能,很多學(xué)者對(duì)該算法進(jìn)行了研究.為了進(jìn)行混沌系統(tǒng)的估計(jì)和連續(xù)函數(shù)優(yōu)化問(wèn)題,LI等人利用正交學(xué)習(xí)策略提高了布谷鳥搜索算法的挖掘能力[17].OUAARAB等人提出了一種改進(jìn)的離散CS算法來(lái)求解著名的旅行商問(wèn)題[18].本文提出了自適應(yīng)搜索平衡布谷鳥搜索算法(Adaptive Search Balanced Cuckoo Search,ASBCS)來(lái)進(jìn)行混合Copula函數(shù)的參數(shù)尋優(yōu),并使用該混合Copula函數(shù),以漢口水文站的年徑流量為因變量,對(duì)宜昌水文站的年徑流量進(jìn)行了預(yù)測(cè).
本文要研究的是兩變量的水文頻率分析和預(yù)測(cè)問(wèn)題,主要使用的是二元Copula函數(shù),所以這里只敘述二元Copula函數(shù)的Sklar定理.
二元Copula函數(shù)Sklar定理:令H(X,Y)為具有邊緣分布F(X)和G(X)的二元聯(lián)合分布函數(shù),則存在一個(gè)Copula函數(shù)C(u,v),滿足
H(X,Y)=C[F(X),G(Y)].
其中u,v均服從[0,1]上的均勻分布;H(X,Y)是具有邊緣分布F(X)和G(X)的二元聯(lián)合分布函數(shù);C(u,v)對(duì)每一個(gè)變量而言,都是單調(diào)非降的,若一個(gè)變量保持不變,則C(u,v)將隨另一個(gè)變量的增大(減小、不變)而增大(減小、不變).
在實(shí)際應(yīng)用中,常見的二維Copula函數(shù)有三種:二維Clayton Copula函數(shù)、二維Gumbel Copula函數(shù)和二維Frank Copula函數(shù).它們的表達(dá)式分別如式(1)、(2)和(3)所示,式中的θ是與Copula函數(shù)的Kendall秩相關(guān)系數(shù)和尾部相關(guān)性相關(guān)的一個(gè)重要參數(shù).
(1)
(2)
(3)
在圖1的(a)、(b)和(c)三個(gè)子圖中,分別展示了三個(gè)函數(shù)的概率密度圖,從圖1中可以看出Clayton Copula函數(shù)的尾部特征不對(duì)稱,其密度函數(shù)很像一個(gè)“L”,上尾低下尾高, 所以它善于敏銳地捕捉尾部不對(duì)稱的數(shù)據(jù)在下尾處的變化,得到下尾處的相關(guān)性特征.Gumbel Copula函數(shù)的尾部特征也不對(duì)稱,其密度函數(shù)很像一個(gè)“丁”,下尾低上尾高, 所以它善于敏銳地捕捉尾部不對(duì)稱的數(shù)據(jù)在上尾處的變化,得到上尾處的相關(guān)性特征.Frank Copula函數(shù)的尾部特征是對(duì)稱的,所以無(wú)法捕捉到數(shù)據(jù)尾部不對(duì)稱的相關(guān)關(guān)系,反而對(duì)于對(duì)稱的尾部關(guān)系很敏感.
(a) Clayton (b) Gumbel (c) Frank Copula圖1 三個(gè)函數(shù)的概率密度圖Fig.1 Probability density graph of three functions
通過(guò)上述分析可知,三種函數(shù)分別對(duì)對(duì)稱的、上尾高和下尾高的數(shù)據(jù)之間尾部相關(guān)關(guān)系比較敏感,所以它們分別有不同的應(yīng)用場(chǎng)景.然而,現(xiàn)實(shí)生活中的數(shù)據(jù)都是復(fù)雜的,單一的Copula函數(shù)往往無(wú)法全面的描述數(shù)據(jù)之間的相關(guān)關(guān)系,基于此,本文將根據(jù)三個(gè)函數(shù)的特性,設(shè)計(jì)混合Copula函數(shù),來(lái)應(yīng)對(duì)現(xiàn)實(shí)世界中復(fù)雜多樣的數(shù)據(jù).
布谷鳥搜索算法必須滿足如下3個(gè)假設(shè):①布谷鳥每次產(chǎn)一枚卵,并隨機(jī)將其放入一個(gè)鳥巢;②最好的巢才能被保留到下一代;③布谷鳥可以搜索的巢的數(shù)量是固定的, 宿主鳥巢主人發(fā)現(xiàn)外來(lái)卵的發(fā)現(xiàn)概率是P, 鳥蛋如果被發(fā)現(xiàn),宿主將移除該蛋或直接遺棄該鳥巢.基于上述假設(shè),布谷鳥尋找鳥巢位置的更新公式為:
xt+1,i=xt,i+α?levy(β),i=1,2,…,n.
(4)
其中xt,i表示經(jīng)過(guò)t代更新之后第i個(gè)鳥巢的位置;?表示點(diǎn)乘;levy(β)表示Levy飛行隨機(jī)游走的路線,可以用一個(gè)簡(jiǎn)單的冪律公式表示為levy(β)~μ=t-1-β;α表示步長(zhǎng)因子,設(shè)α0是一個(gè)常數(shù),xbest是全局最優(yōu)位置,步長(zhǎng)因子可以用式(5)來(lái)表示:
α=α0?(xt,i-xbest),
(5)
位置更新之后,對(duì)每一個(gè)xt+1,i,隨機(jī)產(chǎn)生一個(gè)0到1之間的數(shù)r,與P進(jìn)行比較,如果P大于r,令xi=xt+1,i,將xt+1,i丟棄,然后根據(jù)式(6)的隨機(jī)游走策略求得一個(gè)新的位置對(duì)它進(jìn)行替換:
xt+1,i=xi+q(xt+1,k-xt+1,j),
(6)
其中,q是0到1之間的隨機(jī)數(shù),xt+1,k和xt+1,j是位置更新之后的任意兩個(gè)位置.
標(biāo)準(zhǔn)布谷鳥搜索算法的偽代碼如下所示.
Algorithm 1 Cuckoo search algorithm via Levy flightsInput: Objective function f(x)Output: the optimal solution Randomly initialize population of n host nests xiwhile (t 標(biāo)準(zhǔn)的CS算法是一種模型簡(jiǎn)單,算法前期搜索速度較快的群智能算法,它具有如下的幾個(gè)優(yōu)點(diǎn): (1)CS算法中的位置更新是基于Levy飛行策略的.在Levy飛行策略中,短距離的探索和偶爾較長(zhǎng)距離的開發(fā)是相間的.這種位置更新策略能有效的擴(kuò)大搜索范圍,易于跳出局部最優(yōu)點(diǎn),增加種群的多樣性; (2)相比于其他的智能算法(例如遺傳算法、粒子群算法),CS算法參數(shù)少,除了種群規(guī)模,只有1個(gè)參數(shù)P,且算法的收斂速度對(duì)P不敏感,因此其操作簡(jiǎn)單,通用性更強(qiáng). 同時(shí),經(jīng)過(guò)深入研究,標(biāo)準(zhǔn)的CS算法也具有如下一些缺點(diǎn): (1)Levy飛行策略執(zhí)行全局搜索,隨機(jī)游走策略執(zhí)行局部搜索,發(fā)現(xiàn)概率P用于調(diào)節(jié)二者之間的平衡,但是在算法執(zhí)行的后期,由于個(gè)體差異的縮小,算法將比較偏重于局部搜索,從而很容易令算法陷入局部最優(yōu); (2)發(fā)現(xiàn)概率P和步長(zhǎng)α都不能根據(jù)算法收斂的程度而自適應(yīng)變化,因此標(biāo)準(zhǔn)算法的收斂速度慢,且個(gè)體的多樣性差. 基于上述分析,本文提出了ASBCS算法,就如下幾個(gè)方面對(duì)布谷鳥搜索算法進(jìn)行了改進(jìn). (1)對(duì)步長(zhǎng)因子進(jìn)行了改進(jìn),令步長(zhǎng)因子隨著算法的迭代次數(shù)自適應(yīng)的變化. 標(biāo)準(zhǔn)布谷鳥搜索算法的步長(zhǎng)公式如式(5)所示, 其中α0是一個(gè)常數(shù),步長(zhǎng)因子只能隨著位置的更新而發(fā)生不規(guī)律的變化,不能隨著算法迭代次數(shù)的增加而自適應(yīng)調(diào)整,因此,本文在步長(zhǎng)公式里增加了迭代次數(shù)因子,具體如式(7)所示,其中α0max和α0min是根據(jù)實(shí)際情況給出的步長(zhǎng)的最大值和最小值,T是總的迭代次數(shù),t是本次的迭代次數(shù). (7) (2)在算法后期,對(duì)全局搜索和局部搜索進(jìn)行了平衡. 在標(biāo)準(zhǔn)布谷鳥搜索算法的迭代后期,因?yàn)閭€(gè)體之間的差距變小,造成搜索偏向于局部搜索,容易令算法陷入局部最優(yōu),為了改善標(biāo)準(zhǔn)布谷鳥搜索算法的這一缺點(diǎn),對(duì)公式(7)進(jìn)行了改進(jìn),改進(jìn)之后的步長(zhǎng)公式如式(8)所示. (8) 相對(duì)于CS算法,ASBCS算法在兩個(gè)方面進(jìn)行了改進(jìn),提高了算法的性能,主要表現(xiàn)在如下幾個(gè)方面. 首先,步長(zhǎng)因子的調(diào)整,令算法中解的跳躍,不止依賴于當(dāng)前解和全局最優(yōu)解,還與迭代次數(shù)成反比,隨著迭代次數(shù)的增加,令算法更加趨向于局部搜索,對(duì)局部的解區(qū)間進(jìn)行深度挖掘,提高了算法的收斂速度. 其次,在提高收斂速度的同時(shí),為了提高解的多樣性,避免算法陷入局部最優(yōu),在算法迭代的后期采用了全局搜索和局部搜索的平衡策略,以一定的概率P1來(lái)平衡兩種搜索方式,一方面加快了收斂速度,另一方面又令算法不斷的跳出局部最優(yōu),增加解的多樣性,令算法在全局和局部搜索之間靈活的轉(zhuǎn)換,從而達(dá)到兩種搜索的平衡. 為了比較標(biāo)準(zhǔn)CS算法和本文所提出的ASBCS算法的性能,分別采用兩種算法對(duì)一個(gè)經(jīng)典函數(shù)進(jìn)行尋優(yōu),該函數(shù)如式(9)所示,它在(0,0)處取最小值. (9) 用兩種算法對(duì)該函數(shù)進(jìn)行尋優(yōu)的實(shí)驗(yàn)結(jié)果如圖2所示.從圖中可以看出,ASBCS算法比CS算法具有更快的收斂速度和更好的精確度. 圖2 標(biāo)準(zhǔn)CS算法和ASBCS算法性能比較Fig.2 Comparison of performance between standard CS algorithm and ASBCS algorithm 通過(guò)深入分析二維Clayton Copula函數(shù)、二維Gumbel Copula函數(shù)和二維Frank Copula函數(shù)的概率密度圖,發(fā)現(xiàn)這三種函數(shù)分別適合的是下尾相關(guān)、上尾相關(guān)和尾部對(duì)象的兩變量相關(guān)關(guān)系,利用它們各自的優(yōu)勢(shì),可以將它們?nèi)诤显谝黄?,設(shè)計(jì)成混合Copula函數(shù),從而令Copula函數(shù)適合描述各種復(fù)雜的兩變量關(guān)系.混合之后Copula函數(shù)可以描述為式(10): C=ω1CFrank(u,v,α)+ω2CGumbel(u,v,β)+ ω3CClayton(u,v,γ), (10) 將式(1)、(2)和(3)式代入,可以得到具體的混合Copula函數(shù)表達(dá)式如式(11)所示. (11) 式中α、β、γ、ω1、ω2、ω3由改進(jìn)的布谷鳥搜索算法尋優(yōu)得到. 尋優(yōu)時(shí)改進(jìn)的布谷鳥搜索算法的適應(yīng)度函數(shù)為式(12)所示: (12) 本文以漢口站的年徑流量作為自變量,宜昌站的年徑流數(shù)據(jù)作為因變量,以它們1952-2000年的數(shù)據(jù)作為歷史數(shù)據(jù)進(jìn)行訓(xùn)練,以漢口站2001-2010年的數(shù)據(jù)作為自變量,來(lái)預(yù)測(cè)宜昌站2001-2010年的年徑流量.分別采用Clayton Copula函數(shù)、Gumbel Copula函數(shù)、Frank Copula函數(shù)以及本文所提出的混合Copula函數(shù)進(jìn)行預(yù)測(cè),實(shí)測(cè)值與預(yù)測(cè)值的比較如圖3、圖4、圖5和圖6所示. 圖3 Clayton Copula函數(shù)求取的宜昌站年徑流預(yù)測(cè)值與實(shí)測(cè)值對(duì)比Fig.3 Comparison between the predicted annual runoff calculated by Clayton Copula function and the measured annual runoffin Yichang station 圖4 Gumbel Copula函數(shù)求取的宜昌站年徑流預(yù)測(cè)值與實(shí)測(cè)值對(duì)比Fig.4 Comparison between the predicted annual runoff calculated by Gumbel Copula function and the measured annual runoffin Yichang station 從圖3、圖4和圖5中可以看出單一的Copula函數(shù)在某些年份預(yù)測(cè)的精度非常高,實(shí)測(cè)值與預(yù)測(cè)值基本重合,但是在另一些年份預(yù)測(cè)精度又很低,極大一部分原因是由于氣候、人為引流或者建造工程等因素造成的,但是也有一部分原因是由于Copula函數(shù)的性能造成. 圖5 Frank Copula函數(shù)求取的宜昌站年徑流預(yù)測(cè)值與實(shí)測(cè)值對(duì)比Fig.5 Comparison between the predicted annual runoff calculated by Frank Copula function and the measured annual runoffin Yichang station 比較圖6和圖3、圖4、圖5可以發(fā)現(xiàn),相比Clayton、Gumbel和Frank三個(gè)單一的Copula函數(shù),基于ASBCS算法的混合Copula函數(shù)具有更高的預(yù)測(cè)精度. 為了更好地比較四個(gè)Copula函數(shù),將四個(gè)Copula函數(shù)預(yù)測(cè)宜昌站2001-2010年的年徑流量的相對(duì)誤差全部列于表1中,從表1可以看出,混合Copula函數(shù)可以彌補(bǔ)三個(gè)單一Copula函數(shù)的缺陷,比如在2001年Clayton的預(yù)測(cè)精度遠(yuǎn)遠(yuǎn)低于Gumbel和Frank函數(shù),但是在2002年Gumbel和Frank函數(shù)的相對(duì)誤差又遠(yuǎn)遠(yuǎn)高于Clayton函數(shù),而混合Copula函數(shù)因?yàn)閷⑷齻€(gè)函數(shù)結(jié)合起來(lái),取其優(yōu)點(diǎn),所以很好的融合了三種函數(shù),在2001年和2002都取得了更好的預(yù)測(cè)效果.比較四個(gè)預(yù)測(cè)方法的平均相對(duì)誤差,發(fā)現(xiàn)基于ASBCS算法的混合Copula函數(shù)具有更好的預(yù)測(cè)精度. 圖6 基于ASBCS算法的混合Copula函數(shù)求取的宜昌站年徑流預(yù)測(cè)值與實(shí)測(cè)值對(duì)比Fig.6 Comparison between the predicted annual runoff calculated by hybrid Copula function based on ASBCS algorithm and themeasured annual runoff in Yichang station 表1 四個(gè)函數(shù)預(yù)測(cè)宜昌站2001-2010年的年徑流量的相對(duì)誤差比較Tab.1 The relative error comparison of the annual runoff in Yichang station from 2001 to 2010 predicted by four functions 本文提出了基于自適應(yīng)搜索平衡的布谷鳥搜索算法,對(duì)混合Copula函數(shù)進(jìn)行參數(shù)尋優(yōu),以漢口水文站和宜昌水文站的歷史數(shù)據(jù)作為訓(xùn)練集,以10年間的漢口水文站的年徑流量為自變量,對(duì)宜昌水文站10年間的年徑流量進(jìn)行了預(yù)測(cè).本文所構(gòu)建的混合Copula函數(shù)是在深入研究了三種二元Copula函數(shù),詳細(xì)分析了三種函數(shù)的參數(shù)設(shè)置以及適合應(yīng)對(duì)的變量關(guān)系之后所構(gòu)建的.實(shí)驗(yàn)結(jié)果證明了ASBCS算法收斂速度快,尋優(yōu)精度高,用于混合Copula函數(shù)尋優(yōu)之后,對(duì)宜昌水文站10年間的年徑流量的預(yù)測(cè)精度明顯高于三個(gè)單一的Copula函數(shù).本文所提出ASBCS算法對(duì)其他應(yīng)用中的參數(shù)尋優(yōu)具有重要的參考意義,該鄰近水文站的年徑流量預(yù)測(cè)方法可以推廣至任何兩個(gè)地域上較鄰近水文站的年徑流預(yù)測(cè).在今后的研究中,將進(jìn)一步混合Copula函數(shù),研究更優(yōu)的群智能算法來(lái)進(jìn)行混合Copula函數(shù)的參數(shù)尋優(yōu).3 ASBCS算法
3.1 標(biāo)準(zhǔn)CS算法性能分析
1.2 ASBCS算法
3.3 ASBCS算法的性能分析
3.4 ASBCS算法與CS算法的性能比較
4 基于ASBCS算法的混合Copula函數(shù)
5 實(shí)驗(yàn)數(shù)據(jù)與分析
6 結(jié)語(yǔ)