基于EEMD分解與BOA算法優(yōu)化神經(jīng)網(wǎng)絡(luò)的密云水庫大閣水文站徑流預(yù)測

陳 芳，張志強，李 扉*，孫愷琦

(1.北京林業(yè)大學(xué) 理學(xué)院，北京 100083；2.北京林業(yè)大學(xué) 水土保持學(xué)院，北京 100083)

密云水庫作為北京飲用水源地，肩負著向全市供給生活用水和工農(nóng)業(yè)生產(chǎn)用水的重任[1]。配置水資源，以最大程度地滿足城市用水需求，預(yù)測潮河流域水資源時空分布對流域經(jīng)濟持續(xù)發(fā)展、生態(tài)系統(tǒng)良性循環(huán)和水資源合理利用具有指導(dǎo)意義。

目前，對于徑流預(yù)測方法，研究提出了2大類模型：一是過程驅(qū)動模型方法，二是數(shù)據(jù)驅(qū)動模型方法[2]。過程驅(qū)動模型以水文學(xué)概念為基礎(chǔ)，通過模擬河流演進過程和徑流的產(chǎn)生過程來預(yù)測徑流量，可以做到精確預(yù)測，極具代表性的有SWAT模型和新安江模型[3]。應(yīng)用較為廣泛的數(shù)據(jù)驅(qū)動模型主要包括以序列自相關(guān)分析為基礎(chǔ)的時間序列模型(AR，MA，ARMA，ARIMA等)[4]、以人工神經(jīng)網(wǎng)絡(luò)(ANN)[5]、支持向量機(SVM)[6]等為代表的機器學(xué)習模型以及多學(xué)科方法交叉的非線性綜合預(yù)報模型[7]等。人工神經(jīng)網(wǎng)絡(luò)是通過數(shù)學(xué)方法抽象和模擬生物神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)和功能而得到的一種信息處理系統(tǒng)，應(yīng)用于解決復(fù)雜的非線性問題。Muhammad Tayyab等[8]將人工智能模型與分解方法相結(jié)合，在離散小波變換(DWT)和集成經(jīng)驗?zāi)Ｊ椒纸?EEMD)的基礎(chǔ)上，建立了包含前饋-反向傳播(FFBP)和徑向基函數(shù)(RBF)的混合模型，結(jié)果發(fā)現(xiàn)結(jié)合EEMD的RBFNN具有更好的預(yù)測能力，EEMD-RBF能更準確地捕捉汛期徑流時間序列的非線性特征；D.K.Ghose[9]針對婆羅門河的Govindpur流域建立了BP神經(jīng)網(wǎng)絡(luò)(BPNN)徑流預(yù)測模型，得到了良好的預(yù)測結(jié)果；A.B.Dariane等[10]將小波變換和遺傳算法應(yīng)用于輸入選擇，發(fā)現(xiàn)反向傳播神經(jīng)網(wǎng)絡(luò)(ANN-BP)和ELM算法(ANN-ELM)的Nash-Sutcliffe效率(NSE)均得到了提高。

本研究在針對參數(shù)優(yōu)化神經(jīng)網(wǎng)絡(luò)模型方面，摒棄傳統(tǒng)的粒子群算法(PSO)、模擬退火算法(SA)等優(yōu)化算法，首次將新型優(yōu)化算法——蝴蝶優(yōu)化算法(BOA)用于優(yōu)化BP神經(jīng)網(wǎng)絡(luò)建立預(yù)測模型，研究其是否能更好地提高預(yù)測精度，以期為密云水庫水資源的最優(yōu)利用提供理論基礎(chǔ)和科學(xué)依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

密云水庫位于北京市密云區(qū)，是亞洲最大的人工湖，面積約為180 km2。水庫主要由白河和潮河兩大入庫河流組成，潮白河流域是北京市密云水庫的集水地區(qū)，其中潮河是潮白河的重要支流之一[11]。潮河流域位于116°17′-117°12′E,40°27′-41°32′N,發(fā)源于河北省豐寧縣上黃旗哈拉海溝分水嶺,經(jīng)豐寧縣城、灤平縣部分鄉(xiāng)鎮(zhèn)與潮白河水系的另一條支流安達木河匯合后流入密云水庫，干流全長253 km，在承德市境內(nèi)185.8 km，流域面積4 787.39 km2，是北京重要的水源地和生態(tài)屏障[12](圖1)。

圖1 潮河流域概況

1.2 數(shù)據(jù)來源

大閣水文站作為潮河背山區(qū)代表站，控制集水面積約1 850 km2。選取大閣水文站為研究對象，研究數(shù)據(jù)為大閣站年/月均徑流量，選取的時間段為1969-2013年，圖2為大閣站年均徑流量趨勢，圖3為45 a間共計540個月的月均徑流量走勢。

圖2 大閣站年徑流量走勢

圖3 大閣站月度流量走勢

1.3 研究方法

1.3.1 R/S分析法 重標極差R/S分析方法最早是由英國水文學(xué)家H.E.Hurst[13]提出，主要用于對時間序列的持續(xù)性和長程記憶性進行分析，并對其未來發(fā)展趨勢進行預(yù)估[14]。基本步驟為：1)將序列均分為a個長度均為n的相鄰子序列，即a×n=N。記任一子序列為Iα，α=1,2,3，…，a，Iα中的元素為Xm,α，m=1,2,3,...,n。

4)n從3開始取值，重復(fù)步驟(1)～(3)，直到增加至N，從而得到{(R/S)n}序列。

5)取lg(R/S)為因變量，lg(n)為自變量進行OLS線性回歸，即得到lg(R/S)=lg(C)+H×lg(n)，其中H即為估計的Hurst指數(shù)(0

6)當0

1.3.2 蝴蝶算法優(yōu)化BP神經(jīng)網(wǎng)絡(luò)

1.3.2.1 BP神經(jīng)網(wǎng)絡(luò) BP神經(jīng)網(wǎng)絡(luò)是一種具有3層或3層以上結(jié)構(gòu)的層內(nèi)無互連的網(wǎng)絡(luò)[15]，由信息正向傳播和誤差逆向傳播2個過程組成，具有較強的自學(xué)習和自適應(yīng)能力[16]，結(jié)構(gòu)如圖4所示。輸入層負責接收外部信息，傳達給隱含層；信息在隱含層被轉(zhuǎn)換加工后再傳到輸出層；最后信息經(jīng)過輸出層向外界輸出。當期望輸出和實際輸出間存在誤差時，模型開始誤差逆向傳播。誤差經(jīng)輸出層，按誤差梯度下降方法來修正各層閾值和權(quán)值，向隱含層、輸入層逐層反傳。循環(huán)反復(fù)的信息正向傳播和誤差逆向傳播，是各層閾值和權(quán)值不斷調(diào)整的過程，也是網(wǎng)絡(luò)學(xué)習的過程(圖5)，直到誤差減小到可以接受的程度，或者迭代至最大學(xué)習次數(shù)。

圖4 BP神經(jīng)網(wǎng)絡(luò)拓撲結(jié)構(gòu)

圖5 BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練流程

1.3.2.2 蝴蝶算法 蝴蝶優(yōu)化算法(butterfly optimization algorithm，簡稱BOA)是Sankalap Arora等[17]于2018年提出的一種模擬自然界中蝴蝶食物(花蜜)搜尋和交配行為而衍生出的新型群體智能優(yōu)化算法。從生物學(xué)上講，蝴蝶用感官感受器感知空氣中的氣味，來確定花蜜或交配伙伴的位置，并通過覓食過程中的全局搜索策略和局部搜索策略不斷迭代，以在搜索空間中獲得食物源或伴侶的最佳位置。

蝴蝶優(yōu)化算法分為3個階段。

1)初始化階段：定義目標函數(shù)及其解空間，并為BOA中參數(shù)賦值。賦值后算法創(chuàng)建蝴蝶的初始種群以進行優(yōu)化。如生成n個初始解：

xi=lb+r1(ub-lb)

(1)

式中，xi為種群中第i只蝴蝶(i=1,2,...,n)，lb為搜索空間下界，ub為搜索空間上界，r1為[0,1]的隨機數(shù)。

2)迭代階段：首先定義目標函數(shù)，即：

f=cIa

(2)

式中，f表示香味被其他蝴蝶感知到的強度，c是感覺形態(tài)，I是刺激強度，a是依賴于形態(tài)的冪指數(shù)，在[0,1]取值。每次迭代過程中，解空間中的所有蝴蝶都會移動到新的位置，然后評估它們的適應(yīng)度值。首先，計算解空間中不同位置上所有蝴蝶的適應(yīng)度值；然后利用公式計算蝴蝶產(chǎn)生的香味，在全局搜索階段，公式為：

(3)

在局部搜索階段，公式為：

(4)

BOA中使用切換概率p實現(xiàn)全局搜索和局部搜索之間的切換。通過全局搜索策略和局部搜索策略不斷進行迭代，更新蝴蝶位置并獲得最優(yōu)解。

3)最終階段：當?shù)A段滿足終止條件結(jié)束時，算法輸出具有最佳適應(yīng)度的最優(yōu)解。

用蝴蝶優(yōu)化算法優(yōu)化BP神經(jīng)網(wǎng)絡(luò)的初始權(quán)值和閾值，操作步驟為：1)設(shè)置種群規(guī)模、感覺模態(tài)、冪指數(shù)和切換概率等基本參數(shù)；2)將BP神經(jīng)網(wǎng)絡(luò)的全部權(quán)值和閾值組合成向量，映射到種群的不同個體中去，初始化初始解的空間位置；3)建立BP神經(jīng)網(wǎng)絡(luò)，用BP神經(jīng)網(wǎng)絡(luò)對訓(xùn)練集進行訓(xùn)練；4)計算每一個個體的適應(yīng)度值；5)考察這個適應(yīng)度值是否已經(jīng)在可接受的范圍內(nèi)或達到最大更新次數(shù)，若是，則結(jié)束優(yōu)化，輸出權(quán)值和閾值；6)若否，則按蝴蝶優(yōu)化算法繼續(xù)生成新一代個體重復(fù)步驟(4)～(5)，進行循環(huán)。

1.3.3 評價指標 (1)均方根誤差RMSE：

(5)

(2)平均絕對誤差MAE：

(6)

(3)平均絕對百分比誤差MAPE：

(7)

(4)預(yù)報合格率與預(yù)報等級：

(8)

表1 預(yù)報等級判定標準

2 結(jié)果與分析

2.1 R/S分析法研究徑流變化趨勢的持續(xù)性

根據(jù)R/S分析法對大閣站的年徑流變化趨勢持續(xù)性進行研究，作出lg(n)-lg(R/S)線性擬合圖像，并得到Hurst指數(shù)(圖6)。

由圖6可知，決定系數(shù)R2為0.9672，所以大閣站年徑流序列擬合效果較為理想，年徑流Hurst指數(shù)為0.803 0，介于0.65～1，表明大閣站年徑流序列具有強持續(xù)性，年徑流序列未來總體趨勢會與過去相同，即未來還將繼續(xù)呈現(xiàn)與過去一樣的下降趨勢。

根據(jù)分析法對大閣站的月徑流變化趨勢持續(xù)性進行研究，得到各月的lg(n)-lg(R/S)線性擬合結(jié)果從而得到各月Hurst指數(shù)和決定系數(shù)R2(表2)。

表2 大閣站月徑流時間序列Hurst指數(shù)和決定系數(shù)

對表2進行分析，絕大部分R2都>0.93，大閣站月徑流序列擬合效果很理想，除7月外，其余月份的Hurst指數(shù)均>0.65。結(jié)果表明，除7月外，大閣站各月徑流量序列具有強持續(xù)性，7月份的月徑流量序列具有弱持續(xù)性。綜合來看，各月徑流序列未來總體趨勢均會與過去相同，即未來還將繼續(xù)呈現(xiàn)與過去一樣的下降趨勢。最后，由大閣站月徑流和年徑流的R/S分析結(jié)果，可以看出大閣站徑流總體呈下降趨勢。

2.2 大閣站月徑流量預(yù)測結(jié)果分析

2.2.1 BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型 本研究數(shù)據(jù)處理工具為MATLAB R2016a軟件。利用G-P算法計算得到此序列的嵌入維數(shù)m為7，時間延遲τ為5，由此利用相空間重構(gòu)原理對540個月的大閣站月徑流量(m3/s)數(shù)據(jù)進行重構(gòu)，得到新的輸入輸出變量X、Y，帶入BP網(wǎng)絡(luò)進行訓(xùn)練。其中X每一行為1組輸入數(shù)據(jù)，每組輸入數(shù)據(jù)由7個數(shù)組成;Y每一行為X對應(yīng)行的輸出數(shù)據(jù)，每組輸出數(shù)據(jù)由1個數(shù)組成。

取前449個變量當作訓(xùn)練集，后60個當作測試集，建立BP神經(jīng)網(wǎng)絡(luò)模型，結(jié)果如圖7所示。觀察圖7，模型預(yù)測結(jié)果較差，準確率一般，經(jīng)計算，預(yù)報等級為丙等，僅可作為參考使用。

圖7 BPNN模型預(yù)測輸出

2.2.2 EEMD-BOA-BP預(yù)測模型 對1969-2013年540個月的大閣站月徑流量(m3/s)數(shù)據(jù)進行EEMD分解，得到8個IMF分量和1個殘差R(圖8)；分別對IMF分量和殘差R建立BP神經(jīng)網(wǎng)絡(luò)模型，模型輸入、輸出層節(jié)點數(shù)同3.2.1，隱含層節(jié)點數(shù)分別通過試錯法確定，IMF1～IMF8和R的模型結(jié)構(gòu)見表3。

圖8 大閣站月流量數(shù)據(jù)EEMD分解

表3 IMF1～IMF8和R的模型結(jié)構(gòu)

分別取前449個數(shù)據(jù)作為訓(xùn)練數(shù)據(jù),后60個數(shù)據(jù)作為測試數(shù)據(jù)，通過蝴蝶優(yōu)化算法來優(yōu)化分量和殘差的所構(gòu)建的BP神經(jīng)網(wǎng)絡(luò)的初始權(quán)值和閾值。設(shè)置種群規(guī)模為30，進化代數(shù)為200，切換概率為0.8，感覺模態(tài)為0.01，冪指數(shù)為0.1。分量和殘差預(yù)測結(jié)果見圖9。

由圖9可知，IMF3～R這幾個分量擬和效果較好，IMF1與IMF2擬合較差，但相對未分解模型RMSE、MAE和MAPE均有所減小(表5)；將預(yù)測分量和殘差結(jié)果求和得到EEMD-BOA-BPNN模型的最終預(yù)測結(jié)果(圖10)。

圖9 IMF1～R的EEMD-BOA-BP預(yù)測結(jié)果

圖10 EEMD-BOA-BP模型預(yù)測輸出

為了體現(xiàn)BOA算法優(yōu)化BP模型的高效性，圖11為未經(jīng)過優(yōu)化的EEMD-BP模型最終預(yù)測結(jié)果。

圖11 EEMD-BPNN模型預(yù)測輸出

為了更直觀地觀察3種模型預(yù)測性能的優(yōu)劣，表4為3種模型的評價指標的對比。

從表4可以看出，在模型預(yù)測效果上，EEMD-BOA-BP>EEMD-BP>BPNN。其中，BPNN模型預(yù)報等級為丙等，雖不能用于預(yù)報作業(yè)，但可以作為參考使用，同時EEMD-BPNN模型的RMSE、MAE、MAPE等指標均小于BPNN模型，說明對月徑流進行EEMD分解，再進行重構(gòu)-預(yù)測，能夠提高一定的預(yù)測精度。經(jīng)過BOA算法優(yōu)化后的模型預(yù)報等級為乙等，能夠用于預(yù)報作業(yè)，并且EEMD-BOA-BPNN模型的RMSE、MAE、MAPE等指標值最小。

表4 3種月徑流模型預(yù)測評價指標對比

3 結(jié)論與討論

水文時間序列徑流預(yù)測作為如今的研究熱點，耦合信號分解方法的機器學(xué)習組合模型進行時序預(yù)測成為了一種新的趨勢[18]。由于徑流序列非平穩(wěn)的特性，利用EEMD分解法可以自適應(yīng)地將非平穩(wěn)序列分解成不同時間尺度的分量，使徑流數(shù)據(jù)趨于平穩(wěn)化[19]。BP神經(jīng)網(wǎng)絡(luò)基于經(jīng)驗風險最小化原則，以誤差反向傳播算法為基礎(chǔ)，具有很好的非線性映射能力。但由于BP神經(jīng)網(wǎng)絡(luò)參數(shù)選取對預(yù)測結(jié)果影響較大，并且傳統(tǒng)的參數(shù)選取耗時耗力，過分依賴于經(jīng)驗選取，具有很大的隨機性，造成預(yù)報結(jié)果不準確，引入BOA算法對BP模型進行參數(shù)尋優(yōu)，可以很大程度地避免上述缺點。

本研究以潮河流域大閣水文站數(shù)據(jù)為研究對象，利用非線性分析方法研究了徑流序列的趨勢和持續(xù)性，運用G-P算法與關(guān)聯(lián)維數(shù)法求得水文站月徑流時間序列相空間重構(gòu)的延遲時間和嵌入維數(shù)，反應(yīng)出徑流動力系統(tǒng)內(nèi)部的混沌特性[20]。建立3層BP神經(jīng)網(wǎng)絡(luò)徑流預(yù)測模型，以及結(jié)構(gòu)為分解-重構(gòu)-預(yù)測的組合模型，并利用BOA算法優(yōu)化BP神經(jīng)網(wǎng)絡(luò)模型，結(jié)果表明：1)單一的BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型預(yù)測結(jié)果較差，僅供參考；2)組合模型EEMD-BPNN模型相較于單一BP模型預(yù)測精度提高；3)經(jīng)過優(yōu)化后的新型組合預(yù)測模型：EEMD-BOA-BP，預(yù)測精度最高，且模型復(fù)雜度適中，易于實現(xiàn)，提高了徑流預(yù)測模型的效率和精度，可為今后開展潮河及其他流域的徑流預(yù)測提供科學(xué)依據(jù)，具有廣闊的應(yīng)用前景。

本研究還存在一定的改進空間，如何進一步改進蝴蝶優(yōu)化算法，減少模型運行時間，提高模型精度是接下來需要進一步研究的內(nèi)容。