洪 敏,艾 萍,2,岳 兆 新
(1.河海大學(xué) 水文水資源學(xué)院,江蘇 南京 210098; 2.河海大學(xué) 計(jì)算機(jī)與信息學(xué)院,江蘇 南京 211100; 3.南京工業(yè)職業(yè)技術(shù)大學(xué) 計(jì)算機(jī)與軟件學(xué)院,江蘇 南京 210023)
流域水文系統(tǒng)是一個(gè)復(fù)雜的動(dòng)態(tài)系統(tǒng),影響中長(zhǎng)期徑流過程變化的水文要素具有復(fù)雜的時(shí)空變異性。目前還沒有一個(gè)通用、完善的預(yù)測(cè)模型可適用于所有情況及地區(qū)的中長(zhǎng)期徑流預(yù)測(cè),且預(yù)測(cè)模型中的智能算法選擇也在很大程度上影響中長(zhǎng)期徑流預(yù)測(cè)效果。一直以來,以BP神經(jīng)網(wǎng)絡(luò)(Backpropagation Neural Networks,BPNN)為代表的智能算法廣泛應(yīng)用于徑流預(yù)測(cè)領(lǐng)域[1-3],但在訓(xùn)練過程中,模型參數(shù)難以準(zhǔn)確確定,因而整體預(yù)測(cè)效果欠佳。相比上述模型,極限學(xué)習(xí)機(jī)[4-5](Extreme Learning Machine,ELM)具有模型結(jié)構(gòu)簡(jiǎn)單、通用性好、計(jì)算速度快等優(yōu)點(diǎn),被廣泛應(yīng)用于干旱、水文預(yù)報(bào)等諸多領(lǐng)域。但該模型的參數(shù)選取具有隨機(jī)性,導(dǎo)致模型的魯棒性不強(qiáng),尤其是部分隱含層節(jié)點(diǎn)在實(shí)際應(yīng)用中可能無效[6]。Yang[7]受自然界中花朵授粉行為啟發(fā),提出了花授粉算法(Flower Pollination Algorithm,F(xiàn)PA)。該算法具有結(jié)構(gòu)簡(jiǎn)單、控制參數(shù)少、魯棒性強(qiáng)及計(jì)算效率高等優(yōu)點(diǎn)[8]。相比于遺傳算法(Genetic Algorithm,GA)和粒子群(Particle Swarm Optimization,PSO)等生物啟發(fā)式算法,F(xiàn)PA不僅具有較強(qiáng)的收斂速度和尋優(yōu)能力,而且在收斂精度和搜索能力等方面更有優(yōu)勢(shì)[9]。
因此,本文提出了一種花授粉算法優(yōu)化極限學(xué)習(xí)機(jī)模型的中長(zhǎng)期徑流預(yù)測(cè)方法(FPA-ELM),并應(yīng)用于雅礱江流域中長(zhǎng)期徑流預(yù)測(cè)中。通過與BPNN、支持向量機(jī)(Support Vector Machine,SVM)、ELM和GA-ELM等數(shù)據(jù)驅(qū)動(dòng)模型進(jìn)行比較分析,驗(yàn)證了本文所提算法具有更好的預(yù)測(cè)效果,可為基于智能算法的中長(zhǎng)期徑流變化趨勢(shì)預(yù)測(cè)提供借鑒。
本文以中長(zhǎng)期徑流預(yù)測(cè)分析為目標(biāo),主要分為數(shù)據(jù)組織、因子篩選、預(yù)測(cè)模型構(gòu)建和結(jié)果評(píng)價(jià)4個(gè)方面。數(shù)據(jù)組織主要包括流域徑流整體趨勢(shì)變化因子構(gòu)造和氣候及降雨數(shù)據(jù)處理;因子篩選是通過信息熵方法計(jì)算得到影響徑流的顯著因子集合;預(yù)測(cè)模型為花授粉算法優(yōu)化極限學(xué)習(xí)機(jī)模型,構(gòu)建FPA-ELM模型,完成中長(zhǎng)期徑流預(yù)測(cè);模型評(píng)估是選用水文預(yù)報(bào)領(lǐng)域常用的評(píng)價(jià)指標(biāo),綜合評(píng)價(jià)預(yù)測(cè)模型的性能。
模型的輸入數(shù)據(jù)主要包括兩個(gè)來源:① 基于流域地形與站點(diǎn)空間分布的位置特征,構(gòu)建能夠反映流域徑流整體趨勢(shì)變化的綜合因子;② 影響流域徑流產(chǎn)生的降雨數(shù)據(jù)與氣候相關(guān)數(shù)據(jù)。
依據(jù)粒計(jì)算理論,以水文站點(diǎn)的流域控制面積占比為權(quán)重乘以每個(gè)站點(diǎn)的月平均徑流量,構(gòu)造流域徑流整體趨勢(shì)變化因子(Comprehensive Runoff Index,COM),用以描述流域徑流的情勢(shì)變化。流域徑流整體趨勢(shì)變化因子構(gòu)造方法如下[10-11]:
(1)
(2)
式中:Wi為第i個(gè)水文站點(diǎn)的權(quán)重;Qi為第i個(gè)水文站點(diǎn)的控制面積百分比;m為流域內(nèi)月均徑流一致性較好的水文站點(diǎn)個(gè)數(shù);Cj為第j個(gè)月的徑流趨勢(shì)變化因子;Cij為第i個(gè)水文站點(diǎn)第j個(gè)月的月均徑流量。
1.2.1偏互信息原理
當(dāng)前,影響徑流過程變化的關(guān)鍵因子的篩選方法主要包括相關(guān)系數(shù)法、先驗(yàn)知識(shí)法、信息熵法及主成分分析法[12-14],且每種方法都有各自適用的領(lǐng)域范圍。其中,偏互信息法是基于信息熵的因子篩選方法,適用于備選因子間的線性和非線性相關(guān)關(guān)系。相較于上述幾個(gè)常用的因子篩選方法,偏互信息法具有減少變量的冗余度、提高因子的篩選速度等優(yōu)點(diǎn),因此更適用于中長(zhǎng)期徑流過程變化的因子選擇。
偏互信息計(jì)算方法如下:
(3)
x′=x-E[x|z]
(4)
y′=y-E[y|z]
(5)
式中:PMI為偏互信息;fX′,Y′(x′,y′)為變量X′與Y′的聯(lián)合概率密度函數(shù);fX′(x′)為X′的邊緣概率密度函數(shù);fY′(y′)為Y′的邊緣概率密度函數(shù);E表示數(shù)學(xué)期望值;x表示輸入變量;y表示預(yù)測(cè)對(duì)象;z表示已入選的輸入變量集合。
偏互信息離散計(jì)算方法如下:
(6)
式中:N為離散樣本個(gè)數(shù);i為觀測(cè)樣本編號(hào);fX′,Y′(xi′,yi′)為(xi′,yi′)處的聯(lián)合概率密度估計(jì);fX′(xi′)為xi′處的邊緣概率密度估計(jì);fY′(yi′)為yi′處的邊緣概率密度估計(jì)。
1.2.2基于PMI的變量選擇步驟
當(dāng)輸入變量有多個(gè)時(shí),由于多個(gè)變量之間可能存在某種相關(guān)關(guān)系,例如X,Z為輸入變量,Y為輸出變量,則輸入變量X與Z之間可能存在相關(guān)關(guān)系,則互信息I(Y,Z)的值可能會(huì)大于實(shí)際值。因此,本文基于條件期望方法將變量Y和Z中包含的有關(guān)X信息剔除,變量間的相關(guān)度通過偏互信息PMI來度量。剔除X后,Y記為u,Z記為v,具體定義為
(7)
u=Y-mY(x)
(8)
v=Z-mZ(x)
(9)
Y和Z之間的偏互信息轉(zhuǎn)化為
PMI(Z,Y)=I(v,u)
(10)
假設(shè)C為輸入變量集,Y為輸出變量,S為預(yù)測(cè)模型的關(guān)鍵輸入變量集合,Cs為候選變量,對(duì)應(yīng)PMI值最大時(shí)的輸入變量集合,那么基于PMI的變量篩選步驟如下:
(1) 初始化S,且S為空集;
(2) 如果C是空集,返回步驟(1);
(3) 計(jì)算u=Y-mY(S);
(4) 對(duì)于Cj∈C的每個(gè)元素,計(jì)算v=Cj-mCj(S);
(5) 如果I(v,u)最大時(shí),選取候選變量Cs;
(6) 計(jì)算赤池信息準(zhǔn)則AIC值,如果AIC值降低,則將Cs移到最優(yōu)輸入變量集為S,返回步驟(2);若AIC值增大則終止篩選。
AIC計(jì)算公式如下:
(11)
式中:ui為根據(jù)已選變量計(jì)算的Y回歸殘差;p為已選變量個(gè)數(shù);n為采樣個(gè)數(shù)。
AIC值隨著自變量的篩選不斷減小,當(dāng)AIC為最小時(shí),最優(yōu)自變量集合篩選完畢。
1.3.1極限學(xué)習(xí)機(jī)
假設(shè)任意給定N個(gè)不同樣本(Xi,ti)。其中,Xi=[xi1,xi2,…,xin]T∈Rn,ti=[ti1,ti2,…,tim]T∈Rm,目標(biāo)函數(shù)定義如下[4]:
(12)
式中:N為樣本總量;g(x)為激活函數(shù);Wi為輸入層與隱含層的權(quán)重矩陣,且Wi=[wi1,wi2,…,win]T;βi為隱含層與輸出層之間的權(quán)重矩陣,且βi=[βi1,βi2,…,βim]T;bi為第i個(gè)隱含層神經(jīng)元的偏置;oj為第j個(gè)樣本的網(wǎng)絡(luò)輸出值;Wi·Xj為Wi和Xj的內(nèi)積;C為隱含層神經(jīng)元個(gè)數(shù)。
預(yù)測(cè)值與真實(shí)值誤差最小,可表示為
(13)
即存在βi,bi,Wi使得:
(14)
用矩陣表示為
Hβ=T
(15)
(16)
1.3.2結(jié)合K折交叉驗(yàn)證與花授粉算法的ELM參數(shù)優(yōu)化
花授粉算法(FPA)是基于花粉傳播的自然過程,依靠其他生物體的攜帶,彼此之間形成一種合作共生關(guān)系[15]。FPA由概率常數(shù)調(diào)節(jié)全局搜索和局部搜索之間的轉(zhuǎn)換,且概率常數(shù)取值范圍為0~1。FPA算法流程如圖1所示,具體算法如下:
(1) 全局搜索的數(shù)學(xué)定義為
(17)
(2) 局部搜索的數(shù)學(xué)定義為
(18)
本文以上述花授粉算法為基礎(chǔ),提出結(jié)合K折交叉驗(yàn)證(K-fold Cross Validation,K-CV)與花授粉算法優(yōu)化極限學(xué)習(xí)機(jī)模型參數(shù)的方法(算法流程如圖2所示),包括以下幾個(gè)部分:
(1) 參數(shù)初始化。假設(shè)訓(xùn)練樣本為[xi,yi](xi∈Rn,n為極限學(xué)習(xí)機(jī)模型的輸入神經(jīng)元數(shù)量,i=1,2,…,N,N為總樣本數(shù)量),構(gòu)造極限學(xué)習(xí)機(jī)的激勵(lì)函數(shù)以及設(shè)置隱含層節(jié)點(diǎn)數(shù)量,其中C為隱含層節(jié)點(diǎn)個(gè)數(shù),g為極限學(xué)習(xí)機(jī)模型的迭代次數(shù)。
(2) 構(gòu)造極限學(xué)習(xí)機(jī)模型的適應(yīng)度函數(shù)。以K-CV的均方根誤差(RMSE)作為極限學(xué)習(xí)機(jī)模型的適應(yīng)度,尋找平均RMSE最小的個(gè)體。
(3) 迭代更新。計(jì)算極限學(xué)習(xí)機(jī)模型的適應(yīng)度大小,并據(jù)此更新個(gè)體。
(4) 極限學(xué)習(xí)機(jī)模型最優(yōu)參數(shù)生成。判斷是否達(dá)到預(yù)先設(shè)定的算法終止要求,如果達(dá)到,則獲得極限學(xué)習(xí)機(jī)模型的最優(yōu)參數(shù)組合;否則,回到步驟(2)。
圖2 FPA-ELM模型流程Fig.2 Flowchart of the FPA-ELM model
模型預(yù)測(cè)性能評(píng)價(jià)通過選用水文預(yù)報(bào)中常用的5個(gè)指標(biāo):平均絕對(duì)百分比誤差(Emape)、均方根誤差(Ermse)、確定性系數(shù)(Edc)、合格率(Eqr)及運(yùn)算時(shí)間(T)等,具體計(jì)算見文獻(xiàn)[10]。
雅礱江流域位于青藏高原東部,流域地形落差大、水量豐沛,水資源豐富,因此開展中長(zhǎng)期徑流預(yù)測(cè)具有十分重要的現(xiàn)實(shí)意義和應(yīng)用價(jià)值。雅礱江流域及其站點(diǎn)分布如圖3所示。
圖3 雅礱江流域及其站點(diǎn)分布Fig.3 Yalong River Basin and its stations
本文的試驗(yàn)數(shù)據(jù)包括前期月徑流整體趨勢(shì)變化因子、氣候、降水等624組數(shù)據(jù),時(shí)間為1960年1月至2011年12月。具體包括兩河口、錦屏、官地和二灘等4個(gè)水文站的徑流數(shù)據(jù)、22個(gè)氣象站的降雨數(shù)據(jù)以及相關(guān)氣候數(shù)據(jù)。其中具體的氣候數(shù)據(jù)如表1所列。
選取兩河口、錦屏、官地和二灘4個(gè)水文站點(diǎn)的流域控制面積占比,根據(jù)公式(1)和(2)計(jì)算每個(gè)水文站點(diǎn)的權(quán)重值(見表2),再以權(quán)重值分別乘以各自站點(diǎn)的月平均徑流量(細(xì)粒度),最終求和得到流域徑流整體趨勢(shì)變化因子(粗粒度)構(gòu)造。
表1 相關(guān)氣候數(shù)據(jù)Tab.1 Related climate data
表2 基于測(cè)站控制面積的流域徑流整體趨勢(shì)變化因子 構(gòu)建權(quán)重Tab.2 Construction of weight of Comprehensive Runoff Index of watershed runoff based on the area controlled by stations
雅礱江流域中長(zhǎng)期徑流預(yù)測(cè)的主要特征因子包括:徑流整體趨勢(shì)變化因子fcom(fcom(t-1),…,fcom(t-11),fcom(t-12));降雨因子fr(fr(t-1),…,fr(t-11),fr(t-12));氣候因子fc1(fc1(t-1),…,fc1(t-11),fc1(t-12)),fc2(fc2(t-1),…,fc2(t-11),fc2(t-12)),…,fc21(fc21(t-1),…,fc21(t-11),fc21(t-12)),共計(jì)23×12=276個(gè)。
本文在偏互信息方法基礎(chǔ)上,結(jié)合了人工篩選,獲得預(yù)測(cè)模型的關(guān)鍵特征因子集,具體流程如下:
(1) 基于PMI方法計(jì)算相關(guān)性排名前20的備選因子(見表3)。
表3 基于PMI方法的備選因子相關(guān)性大小(排名前20)Tab.3 Correlation of candidate factors based on PMI method(top 20)
(2) 以上述結(jié)果為基礎(chǔ),通過人工挑選方式分別從徑流整體趨勢(shì)變化因子、降雨因子和氣候因子三類對(duì)象中選取排序前5的關(guān)鍵特征因子集合,形成預(yù)測(cè)模型的最終特征因子輸入,共計(jì)13個(gè):徑流整體趨勢(shì)變化因子fcom(t-12),fcom(t-1),fcom(t-11),fcom(t-2),fcom(t-3);降雨因子fr(t-1),fr(t-7),fr(t-12);氣候因子fc1(t-1),fc15(t-6),fc3(t-7),fc13(t-8),fc16(t-1)。
2.4.1數(shù)據(jù)集劃分
基于FPA-ELM模型的月徑流預(yù)測(cè)數(shù)據(jù)集劃分為兩個(gè)部分(7折交叉驗(yàn)證和模型測(cè)試)。其中,用于7-CV的數(shù)據(jù)為1960年1月至2001年12月共504組樣本數(shù)據(jù)(隨機(jī)選取6組用于訓(xùn)練,余下1組用于驗(yàn)證模型),測(cè)試期2002年1月至2011年12月共120組樣本數(shù)據(jù)。
2.4.2參數(shù)設(shè)置
不同算法的參數(shù)初始化設(shè)置如下:
(1) FPA:種群大小為90,最大迭代次數(shù)為600,常數(shù)P為0.8,γ為1,λ為1.5,適應(yīng)度函數(shù)采用K-CV的平均RMSE,極限學(xué)習(xí)機(jī)模型的激活函數(shù)選擇“sigmoid”函數(shù)。
(2) GA:種群大小、適應(yīng)度函數(shù)及ELM的激活函數(shù)類似FPA,最大遺傳代數(shù)設(shè)置為600,交叉概率為0.7,變異概率為0.01。
(3) BPNN采用7-CV方法,且結(jié)構(gòu)與極限學(xué)習(xí)機(jī)模型一致,BPNN模型的訓(xùn)練函數(shù)采用“tansig”函數(shù),學(xué)習(xí)函數(shù)采用“l(fā)ogsig”函數(shù),最大訓(xùn)練次數(shù)設(shè)置為1 000,學(xué)習(xí)速率設(shè)置為0.1,訓(xùn)練算法采用LM(Levenberg-Marquardt)算法,動(dòng)量因子大小設(shè)置為0.9,模型的期望誤差設(shè)置為0.001。
(4) SVM模型采用RBF(Radial Basis Function)核函數(shù),其中σ=0.5,懲罰參數(shù)C=1,ε=0.001。
2.4.3預(yù)測(cè)結(jié)果及對(duì)比分析
FPA-ELM與BPNN,SVM,ELM和GA-ELM等4種對(duì)比模型的交叉驗(yàn)證期和測(cè)試期的性能對(duì)比如表4所列。不同模型的預(yù)測(cè)結(jié)果對(duì)比如圖4所示,預(yù)測(cè)值和觀測(cè)值及確定性系數(shù)Edc對(duì)比如圖5所示。
表4 不同模型在交叉驗(yàn)證期和測(cè)試期的性能比較Tab.4 Performance comparison of different models in cross validation period and test period
圖4 不同模型的預(yù)測(cè)結(jié)果對(duì)比Fig.4 Comparison of prediction results of different models
圖5 不同模型的預(yù)測(cè)值和觀測(cè)值對(duì)比Fig.5 Comparison on observed and predicted values of different models
上述試驗(yàn)結(jié)果表明:5種數(shù)據(jù)驅(qū)動(dòng)預(yù)測(cè)模型均具有較好的預(yù)測(cè)效果,其中生物啟發(fā)式算法優(yōu)化極限學(xué)習(xí)機(jī)模型在Emape,Ermse和Edc這幾個(gè)指標(biāo)中整體上優(yōu)于另外3種模型,特別是FPA-ELM模型的性能最優(yōu),驗(yàn)證了本文所提算法的優(yōu)越性。其中:在Emape評(píng)價(jià)指標(biāo)中,F(xiàn)PA-ELM模型效果最好,BPNN和SVM效果較差;在Ermse評(píng)價(jià)指標(biāo)中,GA-ELM和FPA-ELM模型較小,SVM模型較大;在Edc評(píng)價(jià)指標(biāo)中,5種數(shù)據(jù)驅(qū)動(dòng)模型的確定系數(shù)都超過了0.85,顯示了基于信息熵因子篩選方法的優(yōu)勢(shì),其中FPA-ELM模型在該項(xiàng)指標(biāo)方面表現(xiàn)最佳;在Eqr評(píng)價(jià)指標(biāo)中,盡管FPA-ELM模型可以應(yīng)用于水文作業(yè),但5種預(yù)測(cè)模型的整體合格率較低,其主要原因在于月徑流預(yù)測(cè)的預(yù)見期較長(zhǎng)且影響中長(zhǎng)期徑流預(yù)測(cè)的主要對(duì)象及其特征因子較多,導(dǎo)致預(yù)測(cè)的不確定性提高,因而影響整體的中長(zhǎng)期徑流預(yù)報(bào)合格率;在運(yùn)算速度評(píng)價(jià)中,F(xiàn)PA-ELM和GA-ELM兩種生物啟發(fā)式模型明顯優(yōu)于其他3種模型,其中FPA-ELM模型表現(xiàn)最佳。
綜上所述,相比其他4種常用的數(shù)據(jù)驅(qū)動(dòng)模型,本文所提出的花授粉算法優(yōu)化極限學(xué)習(xí)機(jī)模型(FPA-ELM)在中長(zhǎng)期徑流預(yù)測(cè)方面更有優(yōu)勢(shì)。主要原因在于相較于傳統(tǒng)的非線性預(yù)測(cè)模型(BPNN,SVM),ELM具有模型結(jié)構(gòu)簡(jiǎn)單、通用性好、計(jì)算速度快等優(yōu)點(diǎn)。此外,本文結(jié)合K折交叉驗(yàn)證與花授粉算法,克服了傳統(tǒng)ELM的不足,且相比于GA算法,F(xiàn)PA具有較強(qiáng)的收斂速度和尋優(yōu)能力,能夠快速搜索ELM的最優(yōu)參數(shù),因而整體預(yù)測(cè)效果較好。
(1) 構(gòu)造了反映流域水情豐枯變化的流域徑流整體趨勢(shì)變化因子(COM),并采用偏互信息法獲得了影響中長(zhǎng)期徑流過程變化的關(guān)鍵因子集,形成中長(zhǎng)期徑流預(yù)測(cè)模型的關(guān)鍵特征因子輸入。
(2) 相較于BPNN,SVM,ELM和GA-ELM等數(shù)據(jù)驅(qū)動(dòng)模型,結(jié)合K折交叉驗(yàn)證與花授粉算法優(yōu)化ELM參數(shù)構(gòu)建的FPA-ELM模型,在Emape,Ermse,Edc,Eqr和運(yùn)算時(shí)間等性能評(píng)價(jià)指標(biāo)方面整體上優(yōu)于上述4種預(yù)測(cè)模型,具有更好的預(yù)測(cè)效果,可為基于智能算法的中長(zhǎng)期徑流預(yù)測(cè)提供借鑒。