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

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

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

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

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

本研究在針對參數(shù)優(yōu)化神經(jīng)網(wǎng)絡模型方面，摒棄傳統(tǒng)的粒子群算法(PSO)、模擬退火算法(SA)等優(yōu)化算法，首次將新型優(yōu)化算法——蝴蝶優(yōu)化算法(BOA)用于優(yōu)化BP神經(jīng)網(wǎng)絡建立預測模型，研究其是否能更好地提高預測精度，以期為密云水庫水資源的最優(yōu)利用提供理論基礎和科學依據(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分析方法最早是由英國水文學家H.E.Hurst[13]提出，主要用于對時間序列的持續(xù)性和長程記憶性進行分析，并對其未來發(fā)展趨勢進行預估[14]。基本步驟為：1)將序列均分為a個長度均為n的相鄰子序列，即a×n=N。記任一子序列為Iα，α=1,2,3，…，a，Iα中的元素為Xm,α，m=1,2,3,...,n。

4)n從3開始取值，重復步驟(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)絡

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

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

圖5 BP神經(jīng)網(wǎng)絡訓練流程

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

蝴蝶優(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]取值。每次迭代過程中，解空間中的所有蝴蝶都會移動到新的位置，然后評估它們的適應度值。首先，計算解空間中不同位置上所有蝴蝶的適應度值；然后利用公式計算蝴蝶產(chǎn)生的香味，在全局搜索階段，公式為：

(3)

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

(4)

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

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

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

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

(5)

(2)平均絕對誤差MAE：

(6)

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

(7)

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

(8)

表1 預報等級判定標準

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 大閣站月徑流量預測結(jié)果分析

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

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

圖7 BPNN模型預測輸出

2.2.2 EEMD-BOA-BP預測模型 對1969-2013年540個月的大閣站月徑流量(m3/s)數(shù)據(jù)進行EEMD分解，得到8個IMF分量和1個殘差R(圖8)；分別對IMF分量和殘差R建立BP神經(jīng)網(wǎng)絡模型，模型輸入、輸出層節(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ù)作為訓練數(shù)據(jù),后60個數(shù)據(jù)作為測試數(shù)據(jù)，通過蝴蝶優(yōu)化算法來優(yōu)化分量和殘差的所構(gòu)建的BP神經(jīng)網(wǎng)絡的初始權(quán)值和閾值。設置種群規(guī)模為30，進化代數(shù)為200，切換概率為0.8，感覺模態(tài)為0.01，冪指數(shù)為0.1。分量和殘差預測結(jié)果見圖9。

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

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

圖10 EEMD-BOA-BP模型預測輸出

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

圖11 EEMD-BPNN模型預測輸出

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

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

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

3 結(jié)論與討論

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

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

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