李祿德，崔東文

(1.云南省水文水資源局文山分局，云南 文山 663000；2.云南省文山州水務(wù)局，云南 文山 663000)

水文時間序列預(yù)報是揭示水文自身演變規(guī)律性的一種有效而可靠的方法，同時是水文預(yù)報研究領(lǐng)域的重要內(nèi)容之一。提高水文流時間序列預(yù)報精度對區(qū)域水資源開發(fā)利用、防洪抗旱規(guī)劃、水資源管理保護、水庫優(yōu)化調(diào)度等具有重要意義。為提高水文預(yù)測精度，多元回歸分析[1]、馬爾科夫鏈模型[2]、ARIMA模型[3]、人工神經(jīng)網(wǎng)絡(luò)[4]、支持向量回歸機(SVM)[5]、隨機森林(RF)[6]、相關(guān)向量機(RVM)[7]、長短期記憶網(wǎng)絡(luò)(LSTM)[8]、自適應(yīng)神經(jīng)模糊推理系統(tǒng)[9]、組合模型[10]等模型及方法都被廣泛應(yīng)用于徑流與降水預(yù)報。

水文時間序列往往表現(xiàn)出多尺度、非線性等特征，傳統(tǒng)單一預(yù)測方法難以獲得理想的預(yù)報效果。當(dāng)前，基于“分解-預(yù)測-重構(gòu)”思想的多種方法組合預(yù)報模型被廣泛應(yīng)用于水文時間序列預(yù)報，取得較好的預(yù)報效果。桑宇婷等[11]利用互補集合經(jīng)驗?zāi)B(tài)分解(CEEMD)方法和BP神經(jīng)網(wǎng)絡(luò)建立組合預(yù)測模型，將其應(yīng)用于汾河上游月徑流預(yù)測；劉祖發(fā)等[12]利用小波分解(WD)方法和秩次集對模型建立組合預(yù)測模型，將其應(yīng)用于馬口站年總徑流量以及深圳市年總降雨量預(yù)測；王麗麗等[13]融合奇異譜分析(SSA)方法、灰狼優(yōu)化算法、回歸支持向量機模型，提出SSA-GWO-SVR月徑流組合預(yù)測模型；張以晨等[14]建立奇異譜分解-支持向量回歸機耦合模型，將其應(yīng)用于吉林省西部某氣象站月降水量預(yù)測；呂晗芳等[15]建立變分模態(tài)分解(VMD)-最小二乘支持向量回歸機(LSSVM)耦合模型，將其應(yīng)用于上靜游站等多個水文站月徑流預(yù)測；徐冬梅等[16]融合小波包分解(WPD)方法、LSSVM、ARIMA模型，提出WPD-LSSVM-ARIMA降水量組合預(yù)測模型；李繼清等[17]將極點對稱模態(tài)分解(ESMD)方法與Elman神經(jīng)網(wǎng)絡(luò)相結(jié)合，建立ESMD-Elman模型對長江干、支流8站的年、月徑流進行預(yù)報。

然而，上述組合模型[11-17]徑流預(yù)測主要存在以下問題和不足。①分解不充分或分解子序列分量過多。如CEEMD分解雖然能在一定程度上解決模型混疊問題，但仍然存在計算量大、復(fù)雜度高等問題；WD僅對低頻信號再次分解，難以有效將原序列分解為更具規(guī)律的子序列分量；SSA存在子序列分解數(shù)量多、模型計算量大等不足；VMD缺點在于模態(tài)個數(shù)和懲罰因子難以確定，很大程度上影響了VMD的分解精度。②常用于徑流預(yù)測的模型中，BP、Elman神經(jīng)網(wǎng)絡(luò)存在調(diào)節(jié)參數(shù)多、易陷入局部最優(yōu)等缺點；SVM、LSSVM能較好地解決非線性和局部最小值問題，但處理大規(guī)模數(shù)據(jù)時預(yù)測精度較低。

為進一步提高水文時間序列預(yù)報精度，拓展組合模型應(yīng)用范疇，基于小波包分解(Wavelet Packet Decomposition，WPD)和相空間重構(gòu)原理，提出松鼠搜索算法(Squirrel Search Algorithm，SSA)-極限學(xué)習(xí)機(Extreme Learning Machine，ELM)預(yù)報模型，并構(gòu)建鯨魚優(yōu)化算法(WOA)-ELM、灰狼優(yōu)化(GWO)算法-ELM、粒子群優(yōu)化(PSO)算法-ELM模型及未經(jīng)優(yōu)化的ELM模型作對比模型，并通過云南省上果水文站月徑流、月降水量時間序列預(yù)報實例對5種模型進行驗證。

1 研究方法

1.1 小波包分解(WPD)

小波包分解(WPD)是在小波分解(WD)的基礎(chǔ)上改進和發(fā)展而來的一種有效分解方法。WPD能同時對高頻、低頻部分實施信號分解，并能根據(jù)信號特性和分析要求自適應(yīng)地選擇相應(yīng)頻帶與信號頻譜相匹配。對于波動信號，采用WPD能夠凸出信號的細(xì)節(jié)特征，理論上3層小波包分解能夠提取絕大部分信號中的有效信息，逼近任意非線性函數(shù)，從而達到解決實際問題的目的[18-20]。

小波包分解算法公式[18-20]為：

(1)

重構(gòu)算法為：

(2)

1.2 相空間重構(gòu)

相空間重構(gòu)的目的在于恢復(fù)高維空間中時間序列的規(guī)律。根據(jù)塔肯斯定理(Takens theorem)[21]，只要選擇合適的時間延遲和嵌入維，就能重構(gòu)出一個與原系統(tǒng)具有相同拓?fù)湫再|(zhì)的動力學(xué)系統(tǒng)。重構(gòu)后時間序列的長度為M=S-(m-1)τ，S、m、τ分別為時間序列長度、嵌入維數(shù)和延遲時間[22]。重構(gòu)后的相點計算如下：

X=[xi,xi+τ,xi+2τ，…，xi+(m-1)τ]

(3)

式中X——維數(shù)為M×m的矩陣；i=1,2,…，M；x——徑流或降水時間序列。

相空間重構(gòu)的對象是經(jīng)WPD分解的水文序列各分量，為便于將各分量預(yù)測結(jié)果疊加重構(gòu)，確定延遲時間τ=1，即在延遲時間τ=1條件下，利用Cao方法[23]確定各分量嵌入維數(shù)m。Cao是虛假鄰近法(FNN)的改進方法，其優(yōu)點在于：計算時只需要延遲時間τ一個參數(shù)，能夠有效區(qū)分隨機系統(tǒng)和確定性系統(tǒng)，使用較少的數(shù)據(jù)就可以求得嵌入維數(shù)，目前已在各行業(yè)領(lǐng)域得到廣泛應(yīng)用。Cao方法主要通過式(4)計算相空間中的點在不同嵌入維數(shù)下的最鄰近點的距離變化α(i,m)；再利用式(5)計算α(i,m)的均值E(m)，并通過式(6)計算E1(m)的變化情況；最后通過繪制E1(m)-m曲線，當(dāng)E1(m)的變化逐漸達到穩(wěn)定時，由穩(wěn)定處的m即為所求[23]。

(4)

(5)

(6)

式中E(m)——α(i,m)的均值；E1(m)——E(m)的變化情況；τ——延遲時間，取1。

1.3 松鼠搜索算法(SSA)

2018年，Jain M等[24]受飛行松鼠靈活多變的覓食策略啟發(fā)，提出松鼠搜索算法(SSA)，該算法具有尋優(yōu)能力強、收斂速度快等優(yōu)點，已在盲源分離[25]、故障診斷[26]等方面得到應(yīng)用。SSA基本以下假設(shè)：①森林中有n只飛行松鼠，假設(shè)一只松鼠在一棵樹上，飛行松鼠并采取靈活多變的覓食策略從普通樹、(一般食物源)、橡樹(較好食物源)、核桃樹(最佳食物源)上獲得食物資源；②每只飛行松鼠都通過自己的動態(tài)覓食策略來尋找食物；③只有3種類型的樹——普通樹、橡樹和核桃樹。

參考文獻[24-25]，SSA數(shù)學(xué)描述如下。

a)初始化。隨機初始化松鼠在搜索空間中的位置，并進行適應(yīng)度計算和排序。SSA采用式(7)初始化松鼠的位置：

FSi,j=FSL+rand()×(FSU-FSL)

(7)

式中 FSi,j——第i只松鼠在第j維空間的位置；FSU、FSL——搜索空間上、下限值；rand()——[0,1]范圍均勻分布的隨機數(shù)。

b)通過滑行更新位置。SSA中，松鼠通過3種不同的覓食策略進行位置更新。

策略一：處于橡樹上的松鼠通過滑行向核桃樹上移動，該策略下松鼠位置更新描述如下：

(8)

策略二：處于普通樹上的松鼠通過滑行向橡樹上移動，該策略下松鼠位置更新描述如下：

(9)

策略三：處于普通樹上的松鼠通過滑行直接向核桃樹上移動，該策略下松鼠位置更新描述如下：

(10)

式中R3——[0,1]范圍均勻分布的隨機數(shù)；其他參數(shù)意義同上。

c)檢查季節(jié)監(jiān)控條件。松鼠覓食策略受季節(jié)變化的影響而改變，因此引入季節(jié)監(jiān)控條件，防止算法陷入局部最優(yōu)。當(dāng)滿足季節(jié)常數(shù)Sc小于季節(jié)常數(shù)最小值Smin時，說明冬季結(jié)束，將無法在森林中尋找到冬季食物源的普通松鼠個體按式(11)進行隨機遷移：

(11)

1.4 極限學(xué)習(xí)機(ELM)

極限學(xué)習(xí)機(ELM)是一種廣義的單隱層前饋神經(jīng)網(wǎng)絡(luò)，具有較快的學(xué)習(xí)速度和良好的泛化能力。給定M個樣本Xk={xk,yk}，k=1，2，…，M，其中xk為輸入數(shù)據(jù)，yk為真實值，f(·)為激活函數(shù)，隱層節(jié)點為m個，ELM輸出可表示為[27-29]：

(12)

式中oj——輸出值；Wi={ωi1，ωi2，…，ωim}′——輸入層節(jié)點與第i個隱含層節(jié)點的連接權(quán)值；bi——第i個輸入節(jié)點和隱含層節(jié)點的偏值；λi——第i個隱含層節(jié)點與輸出節(jié)點的連接權(quán)值。

1.5 建模流程

基于WPD分解與相空間重構(gòu)的SSA-ELM模型預(yù)測步驟如下。

步驟一小波包分解層數(shù)的確定是小波包分解的關(guān)鍵，分解層數(shù)過大則導(dǎo)致預(yù)報復(fù)雜度增加，分解層數(shù)過小則不利于將水文時序數(shù)據(jù)中的不同頻率特征分量分離。經(jīng)分析，對于實例1960年1月至2005年12月月徑流時序數(shù)據(jù)，利用2層WPD進行分解，得到4個子序列分量[2,0]—[2,3]，見圖1、2；對于實例1959年1月至2005年12月月降水量時序數(shù)據(jù)，在兼顧模型計算規(guī)模和預(yù)報精度條件下，先采用2層WPD進行分解，得到4個子序列分量[2,0]—[2,3]。由于子序列分量[2,1]的預(yù)測效果不理想，再利用WPD繼續(xù)對[2,1]分量進行分解，共得到5個子序列分量[2,0]、[2,1,1]、[2,1,2]、[2,2]、[2,3]，見圖3、4。

圖1 月徑流時序分解

圖2 月徑流WPD分解3D效果

圖3 月降水量時序分解

圖4 月降水量WPD分解3D效果

步驟二為便于各分量預(yù)報結(jié)果重構(gòu)，在延遲時間為1的條件下，采用Cao方法確定各子序列分量的輸入向量(嵌入維)，見表1。并利用月徑流、月降水量前417～432組數(shù)據(jù)作為訓(xùn)練樣本，后120組作為預(yù)報樣本。

表1 WPD分解各子序列分量嵌入維數(shù)及序列長度

步驟三利用訓(xùn)練樣本的均方誤差作為SSA、WOA、GWO、PSO算法優(yōu)化ELM輸入層權(quán)值和隱含層偏值(超參數(shù))的目標(biāo)函數(shù)：

(13)

ELM輸入層權(quán)值、隱含層偏值待優(yōu)化維度表示為：

Dim=input×hid+hid+output×hid+output

(14)

式中 Dim——待優(yōu)化維度；input——ELM網(wǎng)絡(luò)輸入維數(shù)；hid——隱含層神經(jīng)元數(shù)；output——ELM網(wǎng)絡(luò)輸出維數(shù)。

步驟四設(shè)置松鼠種群規(guī)模n，最大迭代次數(shù)tm，捕食者出現(xiàn)概率Pdp，滑向常數(shù)Gc；利用式(7)隨機初始化松鼠位置(WOA、GWO、PSO算法具體優(yōu)化步驟可參考SSA實現(xiàn))。

步驟五根據(jù)式(13)計算每只松鼠位置的適應(yīng)度值并按升序進行排序，具有最小適應(yīng)度值的松鼠處在核桃樹上FSht，適應(yīng)度值較小的3只松鼠處在橡樹上FSat，其余松鼠處在普通樹上FSnt。令迭代次數(shù)t=1。

步驟六基于式(8)—(10)，松鼠通過3種策略進行位置更新。

步驟七檢查季節(jié)監(jiān)控條件，當(dāng)Sc

步驟八計算每只松鼠位置的適應(yīng)度值，比較并保存最佳松鼠位置。

步驟九令t=t+1。判斷t是否等于tm，若是，輸出最佳松鼠位置，算法結(jié)束；否則重復(fù)步驟六至九。

步驟十輸出最佳松鼠位置，即ELM輸入層權(quán)值和隱含層偏值矩陣。利用SSA-ELM模型對各分量進行預(yù)報，預(yù)報結(jié)果疊加重構(gòu)后即得到實例月徑流、月降水量最終預(yù)報結(jié)果。

步驟十一模型評估。利用平均相對誤差(MAPE)、平均絕對誤差(MAE)、均方根誤差(RMSE)、預(yù)報合格率(QR)對各模型進行評估，見式(15)。

(15)

2 實例應(yīng)用

上果水文站設(shè)立于1958年11月，位于云南省西疇縣西灑鎮(zhèn)上果村董金河干流，觀測項目有水位、流量、降雨。董金河屬紅河流域瀘江水系，發(fā)源于硯山縣龍所街水頭寨，流經(jīng)硯山、西疇、廣南、麻栗坡、富寧，注入越南。董金河河長51.9 km，流域集水面積704 km2，氣候以亞熱帶低緯高原季風(fēng)山地氣候為主，降水豐沛，但時空分布不均，導(dǎo)致水旱災(zāi)害頻發(fā)和水資源供需矛盾突出。因此，開展上果水文站徑流和降水時間序列預(yù)報研究，對于緩解流域水資源供需矛盾、提高防洪減災(zāi)應(yīng)對能力、保護水環(huán)境水生態(tài)安全具有重要意義。

2.1 參數(shù)設(shè)置

a)SSA-ELM、WOA-ELM、GWO-ELM、PSO-ELM模型。設(shè)置SSA、WOA、GWO、PSO算法最大迭代次數(shù)tm=100，種群規(guī)模n=50。其中SSA捕食者出現(xiàn)概率Pdp=0.1，滑向常數(shù)Gc=1.9；WOA對數(shù)螺旋形狀常數(shù)b=2；PSO算法慣性權(quán)重wmax、wmin分別取值0.9、0.1，自我學(xué)習(xí)因子、社會學(xué)習(xí)因子c1、c2均取值1.494 45；其他參數(shù)設(shè)置采用各算法默認(rèn)值。ELM網(wǎng)絡(luò)激活函數(shù)選擇sigmoid函數(shù)，隱層數(shù)為2×input-1(input為輸入維數(shù))，輸入層權(quán)值和隱含層偏值搜索空間設(shè)置為[-1,1]，數(shù)據(jù)采用[-1,1]進行歸一化處理。

b)ELM模型。為驗證群體智能算法優(yōu)化ELM網(wǎng)絡(luò)輸入層權(quán)值和隱含層偏值對提高網(wǎng)絡(luò)預(yù)報精度的影響，除算法部分外，ELM參數(shù)設(shè)置同上述SSA-ELM等4種模型。

2.2 目標(biāo)函數(shù)優(yōu)化

基于步驟二、三，利用SSA、WOA、GWO、PSO算法對目標(biāo)函數(shù)minf(w,b)進行尋優(yōu)，尋優(yōu)結(jié)果見表2。

表2 目標(biāo)函數(shù)minf(w,b)尋優(yōu)結(jié)果

從表2可以看出，除月徑流分解分量[2,1]外，SSA對各分解分量目標(biāo)函數(shù)的尋優(yōu)精度均優(yōu)于WOA、GWO、PSO算法，具有較好的尋優(yōu)效果。整體而言，4種算法尋優(yōu)性能由優(yōu)至劣依次是SSA>WOA>GWO>PSO算法。

2.3 預(yù)報結(jié)果及分析

利用SSA、WOA、GWO、PSO算法對目標(biāo)函數(shù)minf(w,b)尋優(yōu)獲得的輸入層權(quán)值和隱含層偏值矩陣分別構(gòu)建SSA-ELM、WOA-ELM、GWO-ELM、PSO-ELM預(yù)報模型，并構(gòu)建未經(jīng)優(yōu)化的ELM模型對各分量進行預(yù)報，將各分量預(yù)報結(jié)果加和重構(gòu)后即得到最終預(yù)報結(jié)果。各模型采用上述MAPE(%)、MAE(m3/s或mm)、RMSE(m3/s或mm)、QP(%)進行評估，結(jié)果見表3；預(yù)報相對誤差、絕對誤差效果見圖5—8。

表3 實例月徑流、月降水量時間序列預(yù)報結(jié)果對比

圖5 月徑流時間序列預(yù)報相對誤差3D

圖6 月徑流時間序列預(yù)報絕對誤差3D

圖7 月降水量時間序列預(yù)報相對誤差3D

圖8 月降水量時間序列預(yù)報絕對誤差3D

依據(jù)表3、圖5—8可以得出以下結(jié)論。

a)SSA-ELM模型對實例月徑流預(yù)報的MAPE為5.32%，分別較WOA-ELM、GWO-ELM、PSO-ELM、ELM對比模型提高22.0%、23.1%、66.8%、88.1%；預(yù)報的MAE為0.078 m3/s，分別較對比模型提高23.5%、24.3%、65.5%、89.8%；預(yù)報的RMSE為0.103 m3/s，分別較對比模型提高21.4%、25.9%、63.1%、89.9%；預(yù)報的QP達97.5%，分別較對比模型提高5.0%、2.5%、15.8%、42.5%。對實例月降水量預(yù)報的MAPE為3.84%，分別較對比模型提高81.3%、83.1%、93.2%、98.8%；預(yù)報的MAE為0.169 mm，分別較對比模型提高68.8%、68.9%、90.8%、97.3%；預(yù)報的RMSE為0.209 mm，分別較對比模型提高68.9%、70.7%、91.4%、97.2%；預(yù)報的QP達95.8%，分別較對比模型提高3.3%、4.1%、10.0%、25.8%。5種模型預(yù)報精度由高至低依次是：SSA-ELM>WOA-ELM>GWO-ELM>PSO-ELM>ELM，與各算法對目標(biāo)函數(shù)優(yōu)化性能排序一致?？梢姡ㄟ^智能算法優(yōu)化ELM輸入層權(quán)值和隱層偏值，可有效提高ELM網(wǎng)絡(luò)性能；智能算法尋優(yōu)效果越好，尋優(yōu)獲得的ELM超參數(shù)越佳，所構(gòu)建的模型預(yù)報精度越高。

b)SSA-ELM、WOA-ELM、GWO-ELM模型對實例月徑流、月降水量均具有較好的預(yù)報效果，表明通過WPD對原序列進行分解，可有效弱化了復(fù)雜環(huán)境對徑流及降水時間序列的影響，降低預(yù)報復(fù)雜度；通過Cao方法確定各分量嵌入維數(shù)，可以科學(xué)構(gòu)建輸入向量；利用SSA、WOA、GWO優(yōu)化ELM超參數(shù)，可以得到更好的ELM輸入層權(quán)值和隱層偏值，將WPD、智能算法和ELM三者結(jié)合，可大大提高水文時間序列預(yù)報精度。其中，SSA-ELM模型預(yù)報效果最佳，尤其對于月降水量預(yù)報。

c)從5種模型綜合預(yù)報效果來看，對于月徑流預(yù)報，WOA-ELM、GWO-ELM、PSO-ELM模型預(yù)報效果優(yōu)于PSO-ELM模型，遠(yuǎn)優(yōu)于ELM模型；對于月降水量預(yù)報，WOA-ELM模型預(yù)報效果優(yōu)于GWO-ELM、PSO-ELM模型，遠(yuǎn)優(yōu)于PSO-ELM模型，遠(yuǎn)遠(yuǎn)優(yōu)于ELM模型。

3 結(jié)論

基于WPD分解和相空間重構(gòu)原理，提出SSA-ELM預(yù)報模型，并建立WOA-ELM、GWO-ELM、PSO-ELM、ELM作對比分析模型，通過云南省上果水文站月徑流和月降水時間序列預(yù)報實例對各模型進行驗證，得到以下結(jié)論。

a)對于實例高維優(yōu)化問題，除月徑流分解分量[2,1]外，SSA對各分解分量目標(biāo)函數(shù)的尋優(yōu)精度均優(yōu)于WOA、GWO、PSO算法，具有較好的尋優(yōu)精度。智能算法尋優(yōu)效果越好，尋優(yōu)獲得的ELM超參數(shù)越佳，所構(gòu)建的模型預(yù)報精度越高。

b)SSA-ELM模型對實例月徑流和月降水量的預(yù)報精度優(yōu)于WOA-ELM、GWO-ELM模型，遠(yuǎn)優(yōu)于PSO-ELM、ELM模型。WPD分解可有效弱化復(fù)雜環(huán)境對水文時間序列的影響，降低預(yù)報復(fù)雜度；Cao方法可科學(xué)確定嵌入維數(shù)和構(gòu)建輸入向量；SSA優(yōu)化ELM超參數(shù)，可以獲得更佳的ELM輸入層權(quán)值和隱層偏值，將三者結(jié)合，可大大提高水文時間序列預(yù)報精度。

c)通過對預(yù)報精度不理想的子序列分量進行再分解，不但有效提高模型預(yù)報精度，而且相對于3層WPD分解，大大降低了模型的復(fù)雜程度，使預(yù)報模型簡潔、高效，具有較好的實際應(yīng)用價值。