譚永杰,王現(xiàn)勛*,段茗續(xù),劉亞茹,姚華明,2,3

(1.長江大學(xué)資源與環(huán)境學(xué)院 油氣地球化學(xué)與環(huán)境湖北省重點實驗室,湖北 武漢 430100;2.中國長江電力股份有限公司,湖北 宜昌 443000;3.智慧長江與水電科學(xué)湖北省重點實驗室,湖北 宜昌 443000)

精準可靠的徑流預(yù)報是防洪減災(zāi)、水資源的配置與管理、電力生產(chǎn)等流域管理和調(diào)度決策的關(guān)鍵依據(jù),也是全球水文水資源領(lǐng)域研究的熱點[1]。徑流過程通常受降水、蒸發(fā)、地形、土地利用及大尺度氣候異常因子等諸多因素的影響,本質(zhì)上是高度動態(tài)的,甚至是混沌的,同時也是非線性和非平穩(wěn)的[2]。近年來,受氣候變化和人類活動的影響,許多流域的水資源狀況發(fā)生了巨大的改變,徑流的一致性因此改變的愈加明顯,精準捕捉日徑流時間序列的動態(tài)過程變得更加困難[3]。為了準確地預(yù)測未來時段的徑流,國內(nèi)外水文預(yù)報研究者提出了多種預(yù)報模型以及預(yù)報方法。在這些水文預(yù)報模型方法中大致可以分為3個類型。第一種類型是基于統(tǒng)計學(xué)原理的水文預(yù)報模型,該模型主要是通過對歷史水文數(shù)據(jù)進行分析和擬合,來預(yù)測未來水文過程的變化趨勢。常用的統(tǒng)計學(xué)方法包括回歸分析[4]、時間序列分析[5]、聚類分析[6]等。這種模型具有簡單易用的優(yōu)點,適用于數(shù)據(jù)量較少的情況,但存在著較大的誤差和不確定性。第二種類型是基于物理學(xué)原理的水文預(yù)報模型,該模型主要是基于水文過程的物理規(guī)律,通過建立數(shù)學(xué)公式來描述水文的變化過程。常用的模型包括新安江模型[7-8]、水箱模型[9]、HBV模型[10]等。這種模型具有較高的準確性和可靠性,但需要大量的數(shù)據(jù)和對水文過程的深入理解。第三種類型是基于人工神經(jīng)網(wǎng)絡(luò)的水文預(yù)報模型,該模型主要是通過模擬人腦神經(jīng)元的工作原理,來建立水文過程的非線性模型。常用的神經(jīng)網(wǎng)絡(luò)模型包括(Back Propagation,BP)反向傳播神經(jīng)網(wǎng)絡(luò)[11]、(Extreme Learning Machines,ELM) 極限學(xué)習(xí)機[12]、(Support Vector Networks,SVM)支持向量機[13]、(Long Short Term Memory,LSTM)長短期記憶神經(jīng)網(wǎng)絡(luò)[14]等。由于深度學(xué)習(xí)算法能更好的模擬非線性、非平穩(wěn)的時間序列,可以更加準確地進行徑流預(yù)報,逐漸成為近年來新的研究熱點。

隨著近幾年學(xué)者的深入研究,不少學(xué)者發(fā)現(xiàn)徑流時間序列是含有多種頻率和趨勢成分的非線性非平穩(wěn)時間序列,其非平穩(wěn)性和非線性增加了預(yù)報的預(yù)測難度[15-17]。不少研究表明將時域分析的數(shù)據(jù)預(yù)處理方法與神經(jīng)網(wǎng)絡(luò)模型組合可以提高徑流預(yù)報的精度。李福興等[18]采用經(jīng)驗?zāi)B(tài)分解算法(Empirical Mode Decomposition,EMD)提取徑流序列中具有物理含義的信號并得到不同時間尺度的多個固有模態(tài)函數(shù)(Intrinsic Mode Functions,IMF),在利用(Autoregressive Integrated Moving Average,ARIMA)自回歸移動平均模型與(Generalized Regression Neural Network,GRNN)廣義回歸神經(jīng)網(wǎng)絡(luò)分別對不同時間尺度的IMF進行模擬,分析徑流未來變化趨勢,預(yù)測結(jié)果表明,其提出結(jié)合EMD模態(tài)分解算法的組合預(yù)測模型能取得較好的預(yù)測效果;王秀杰等[19]利用時變?yōu)V波器的經(jīng)驗?zāi)B(tài)分解(Time-Varying-Filter,Empirical Mode Decomposition,TVF-EMD)分解原始徑流序列,再利用(Grey Relation Analysis,Light Gradient Boosting Machine,GRA-LightGBM)組合模型進行預(yù)測,結(jié)果表明結(jié)合模態(tài)分解算法的組合預(yù)測模型較其他模型擁有更高的精度。劉揚等[20]提出一種改進的集合經(jīng)驗?zāi)B(tài)分解算法(EEMD)將徑流信號分解為若干模態(tài)分量,在將各分量分別輸入門控循環(huán)單元(Gate Recurrent Unit,GRU)神經(jīng)網(wǎng)絡(luò)進行預(yù)測,預(yù)測結(jié)果表明,該組合模型擁有更高的預(yù)測精度,可進一步減少防洪區(qū)的經(jīng)濟損失。但上述所描述的模態(tài)分解算法仍存在一定缺陷,EMD受獨特的分解原理限制,基于EMD的分解理論存在的模態(tài)混疊現(xiàn)象和端點效應(yīng)等問題,且只能被削弱而難以根除;EEMD模型在添加噪音時使用了隨機數(shù),將導(dǎo)致每次分解結(jié)果都存在微弱偏差,結(jié)果難以完全復(fù)現(xiàn)。而基于自適應(yīng)噪聲的完整集合經(jīng)驗?zāi)B(tài)分解(Complete Ensemble Empirical Mode Decomposition with Adaptive Noise,CEEMDAN)[21]方法是一種建立在EMD分解和EEMD分解基礎(chǔ)上的自適應(yīng)白噪聲分解方法,在分解過程中的每一階段都添加有限次的自適應(yīng)高斯白噪聲,不僅克服了EMD分解模態(tài)混疊問題,還解決了EEMD方法在較少的實驗平均次數(shù)下導(dǎo)致重構(gòu)誤差較大以及計算效率低的缺點。胡斯曼[22]將CEEMDAN模態(tài)分解算法運用于月徑流、年徑流的預(yù)測中,且取得了不錯的預(yù)測效果,但在更具非平穩(wěn)性、非線性的日徑流時間序列的運用卻少有涉及,在應(yīng)對更加復(fù)雜的日徑流序列中,還有待進一步探究。

針對上述問題,結(jié)合不同的時域分解算法和LSTM模型,構(gòu)建基于經(jīng)驗?zāi)B(tài)分解(EMD)算法、集合經(jīng)驗?zāi)B(tài)分解(EEMD)算法、完整經(jīng)驗?zāi)B(tài)分解(CEEMDAN)算法并利用(Grid Search,GS)網(wǎng)格搜索參數(shù)尋優(yōu)算法[23]對LSTM模型的超參數(shù)進行尋優(yōu)的組合預(yù)測模型,對日徑流進行預(yù)報。以宜昌水文站為研究區(qū)域,其日徑流時間序列為研究對象,比較不同模態(tài)分解算法構(gòu)成的組合預(yù)測模型的預(yù)測精度,結(jié)果表明結(jié)合CEEMDAN的GS-LSTM組合預(yù)測模型擁有更高的預(yù)測精度。

1 研究方法

1.1 經(jīng)驗?zāi)B(tài)分解(EMD)和集合經(jīng)驗?zāi)B(tài)分解(EEMD)

在處理非平穩(wěn)和非線性時間序列數(shù)據(jù)方面,HUANG等[24]于1998年提出的EMD在處理信號時,將會產(chǎn)生一組固有模態(tài)函數(shù),消除了寄生諧波來表示非線性和非平穩(wěn)信號的需要。眾多文獻表明其在處理如徑流[25]、風(fēng)速[26]非線性時序數(shù)據(jù)時具有明顯的優(yōu)勢。EMD的技術(shù)步驟如下。

步驟一找到原始徑流時間序列S(t)的所有局部最大值和局部最小值,然后用三次樣條曲線擬合S(t)的上包絡(luò)線U(t)和下包絡(luò)線L(t),并求出上下包絡(luò)線的均值系列M(t),見式(1):

(1)

步驟二使用原始時間序列S(t)減去均值序列M(t)得到一個新的時間序列H(t),見式(2):

H(t)=S(t)-M(t)

(2)

步驟三如果M(t)和H(t)滿足終止標(biāo)準,可以得到第一個IMF為c1(t)=M(t)和第一個殘差R為r1(t)=H(t),其終止標(biāo)準為①M(t)接近于0;②H(t)中局部最小極值點和零交叉點的數(shù)量必須相同或最多相差一個。

步驟四如果不滿足終止標(biāo)準,對H(t)重復(fù)步驟一—三,直到得到c1(t)和r1(t)。

步驟五對r1(t)重復(fù)上述步驟一—四,直到得到所有的IMF和殘差。因此,原始時間序列S(t)將按照式(3)被分解。

(3)

式中ci(t)——第i個IMF;rk(t)——第k個殘差項。

盡管EMD技術(shù)在實際應(yīng)用中取得了不錯的效果,但其分解得到的IMF存在模態(tài)混疊現(xiàn)象,給IMF的重構(gòu)帶來一定的困擾。Wu等[27]于2009年提出了一種名為EEMD的噪聲輔助數(shù)據(jù)分析方法,其最大優(yōu)點在于其能夠以自適應(yīng)方式提取信號的各分量及變化趨勢,同時顯著地減少了經(jīng)驗?zāi)B(tài)分解(EMD)方法中存在的模態(tài)混疊現(xiàn)象,克服了小波變換容易造成很多虛假的沒有原有物理意義的諧波的缺點[28]。其通過在分解的過程中多次引入均勻分布的白噪聲,將信號本身的噪聲通過多次人為添加的噪聲掩蓋過去,從而得到更加精準的上下包絡(luò)線。同時對分解結(jié)果進行平均處理,平均處理次數(shù)越多,噪聲給分解帶來的影響就越小。EEMD的技術(shù)步驟如下:①在原時間序列中加入高斯白噪聲,然后得到一個新的時間序列;②分解新的時間序列以獲得每個IMF成分;③連續(xù)重復(fù)①和②,但每次都使用不同的高斯白噪聲;④以分解的相應(yīng)IMF的合計平均值作為最終結(jié)果。

1.2 完整集合經(jīng)驗?zāi)B(tài)分解(CEEMDAN)

基于自適應(yīng)噪聲的完整集合經(jīng)驗?zāi)B(tài)分解(CEEMDAN)方法是一種建立在EMD分解和EEMD分解基礎(chǔ)上的自適應(yīng)白噪聲分解方法,在分解過程中的每一階段都添加有限次的自適應(yīng)高斯白噪聲,不僅克服了EMD分解的模態(tài)混疊問題,還克服了EEMD方法在較少的實驗平均次數(shù)下導(dǎo)致重構(gòu)誤差較大以及計算效率低的缺點。

步驟一通過添加一個白噪聲wi(t)(i=1,2,……,n),信噪比為ε0到原始日徑流時間序列S(t),得到的S(t)被用來進行一次分解,見式(4)。其中t代表不同的時間點,i代表第i次加入的白噪聲,n代表加入白噪聲的總數(shù)。

Si(t)=S(t)+ε0ωi(t)

(4)

(5)

(6)

步驟三自適應(yīng)噪聲項是通過EMD獲得的第一個IMF,加上白噪聲ωi(t)且信噪比為ε1然后加入自適應(yīng)噪聲項得到第一個殘差R1(t),并得到一個新的時間序列。隨即利用式(7)分解得到CEEMDAN的第二個IMF,第二個殘差R2(t)根據(jù)式(8)計算得到。

(7)

(8)

步驟四重復(fù)步驟三,加入新的自適應(yīng)噪聲得到新的時間序列與殘差項,然后分解得到CEEMDAN的第K個IMF,其中ωi(t),(i=1,2,…n),εk,(k=2,3,…K),計算步驟見式(9)、(10)。

EMDk-1(ωi(t)))

(9)

(10)

步驟五最后,CEEMDAN算法的殘差項不能繼續(xù)分解,因為不會超過2個極值點,因此最終殘差項R(t)是一個獨特的趨勢項。獲得的完整IMF和R(t)與原始日徑流時間序列見式(11)。

(11)

1.3 LSTM模型

LSTM[29]網(wǎng)絡(luò)是循環(huán)神經(jīng)網(wǎng)絡(luò)(Recurrent Neural Network,RNN)的一種改進網(wǎng)絡(luò),在時間反饋機制基本不變的情況下,引入了記憶單元(memory cell)和門(gate)的機制,來實現(xiàn)對信息流的存儲和控制。該機制保障了信息在其特殊的神經(jīng)元結(jié)構(gòu)中實現(xiàn)長期性的篩選與保留,從而有效避免橫向深度帶來的信息丟失問題。此外,記憶單元的門控機制也使得梯度下降時,梯度的衰減得到緩解,有效緩解了模型訓(xùn)練時的梯度消失與梯度爆炸的問題。

圖1所示,LSTM單元中有3種類型的門控,分別為輸入門、遺忘門和輸出門。圖中分別使用i、f和o來表示輸入、遺忘和輸出門,⊙代表對應(yīng)元素相乘,W和b表示網(wǎng)絡(luò)的權(quán)重矩陣和偏置向量。

LSTM的計算過程大致可以用式(12)—(16)表示,在時間步為t時,LSTM的隱藏層的輸入和輸出向量分別用xt和ht表示,記憶單元用ct表示。輸入門用于控制網(wǎng)絡(luò)當(dāng)前輸入數(shù)據(jù)xt流入記憶單元的多少,即ct可以保存多少有效信息,見式(12):

it=σ(Wxixt+Whiht-1+bi)

(12)

遺忘門是LSTM單元的核心組成部分,其可以保留部分信息并遺忘部分信息。遺忘門控制自連接單元,可以決定歷史信息中的哪些部分會被丟棄,即上一時刻記憶單元ct-1中的信息對當(dāng)前記憶單元ct的影響見式(13)、(14)。

ft=σ(Wxfxt+Whfht-1+bf)

(13)

ct=ft⊙ct-1+it⊙tanh(Wxcxt+Whcht-1+bc)

(14)

輸出門控制記憶單元ct對當(dāng)前輸出值ht的影響,即記憶單元中的部分信息會在時間步t輸出。輸出門的表達見式(15),LSTM單元在t時刻的輸出ht可以通過式(16)得到。

ot=σ(Wxoxt+Whoht-1+bo)

(15)

ht=ot⊙tanh(ct)

(16)

2 實例分析

2.1 研究區(qū)域及數(shù)據(jù)

本文以長江流域中的宜昌水文站為研究對象,宜昌水文站位于三峽工程下游距三峽水利工程約44 km處,地處長江上游與中游的交界處,被視為干流上中游的咽喉。已知宜昌水文站控制流域面積占到整個流域面積的55%,該站所測徑流多年平均徑流量為4 131億m3,其中汛期平均流量為3 412億m3,占全年徑流量的79%左右[30],研究區(qū)域見圖2。

圖2 宜昌水文站地理位置示意

本文采用長江流域宜昌水文站1952年1月1日至2020年12月31日實測(日)徑流資料進行分析,其徑流數(shù)據(jù)見圖3。

圖3 徑流數(shù)據(jù)

2.2 評估指標(biāo)

在評估模型的性能時,均方根誤差(Root Mean Square Error,RMSE)是預(yù)測值與實際值之間差異的平方的平均值的平方根;平均絕對誤差(Mean Absolute Error,MAE)是預(yù)測值與實際值之間差異的絕對值的平均值;平均絕對百分比誤差(Mean Absolute Percentage Error,MAPE)是預(yù)測值與實際值之間差異的絕對值與實際值之間的比率的平均值。本文采用RMSE、MAE和MAPE統(tǒng)計指標(biāo)來進行評價[31],見式(17)—(19):

(17)

(18)

(19)

2.3 模態(tài)分解

利用上述EMD、EEMD、CEEMDAN 3種模態(tài)分解算法對宜昌水文站1952年1月1日至2020年12月31日的日徑流時間序列進行分解,使用EMD將原始序列自適應(yīng)的分解為若干個子模態(tài),結(jié)果見圖4;使用EEMD、CEEMDAN對徑流時間序列進行分解,分解內(nèi)容見圖5、6。

圖4 EMD模態(tài)分解

圖5 EEMD模態(tài)分解

圖6 CEEMDAN模態(tài)分解

2.4 預(yù)報結(jié)果與分析

利用上述模態(tài)分解算法對宜昌水文站1952年1月1日至2020年12月31日的日徑流數(shù)據(jù)進行分解,按照1952年1月1日至1999年10月22日為訓(xùn)練數(shù)據(jù),1999年10月23日至2020年12月31日為測試集,并將分解后的模態(tài)分別輸入LSTM模型進行預(yù)測,并將各模態(tài)的預(yù)測結(jié)果進行加和得到最終的預(yù)測值,本算例直接選取10 d的歷史徑流數(shù)據(jù)預(yù)測下1 d的徑流,考慮到預(yù)見期過短,無法提供太多參考價值,因此采用滾動預(yù)報的方式去增長日徑流預(yù)見期,并利用網(wǎng)格搜索參數(shù)尋優(yōu)算法選擇的超參數(shù)設(shè)置見表1,其預(yù)測流程見圖7。

表1 LSTM超參數(shù)及取值

圖7 預(yù)測流程

為驗證模型預(yù)測性能,對以下4種模型的結(jié)果分別進行討論,以下4種模型均按照上述網(wǎng)格搜索算法尋優(yōu)的超參數(shù)進行設(shè)置,4種模型的評價指標(biāo)結(jié)果見表2、圖8。

表2 4種模型評價指標(biāo)比較

圖8 預(yù)測結(jié)果對比

模型一(GS-LSTM):將原始日徑流時間序列直接輸入LSTM模型進行預(yù)測。

模型二(EMD-GS-LSTM):運用EMD模態(tài)分解算法先將原始日徑流時間序列進行分解,然后分別輸入進LSTM模型進行預(yù)測,最終將各個預(yù)測值進行加和得到徑流預(yù)測值。

模型三(EEMD-GS-LSTM):使用EEMD模態(tài)分解算法先將原始日徑流時間序列進行分解,然后分別輸入進LSTM模型進行預(yù)測,最終將各個預(yù)測值進行加和得到徑流預(yù)測值。

模型四(CEEMDAN-GS-LSTM):利用CEEMDAN模態(tài)分解算法先將原始日徑流時間序列進行分解,然后分別輸入進LSTM模型進行預(yù)測,最終將各個預(yù)測值進行加和得到徑流預(yù)測值。

從表2、圖8可以看出,RMSE在463.46～1 325.23,MAE在327.58～786.47,MAPE在2.90%～5.78%,均處于較低水平,LSTM模型能較好地預(yù)測日徑流的總體趨勢,此外可以看出模型四(CEEMDAN-LSTM)的預(yù)測精度優(yōu)于其他4個模型的精度,說明其在日徑流預(yù)報中擁有更好的表現(xiàn)。

a)將模型一(GS-LSTM)的預(yù)測結(jié)果與其他模型預(yù)測結(jié)果對比,可以看出將模態(tài)分解算法加入到日徑流的預(yù)測,可以提高模型的預(yù)測精度,與模型四相比,其RMSE、MAE和MAPE分別降低了65.02%、58.35%與2.88%,有效降低了徑流時間序列的非平穩(wěn)性與波動性,使得徑流序列更易被LSTM模型捕捉,從而取得更好的預(yù)報效果。

b)將模型二(EMD-GS-LSTM)與模型三(EEMD-GS-LSTM)的預(yù)測結(jié)果進行比較,模型二的預(yù)測結(jié)果略微優(yōu)于模型三,其RMSE、MAE和MAPE分別降低了20.91%、39.23%與1.57%。在使用EEMD算法對徑流時間序列進行分解時,由于隨機噪聲的加入,導(dǎo)致在模態(tài)分解時存在微弱的偏差,進而影響到預(yù)測結(jié)果,這表明傳統(tǒng)的模態(tài)分解算法在進行時域分解時仍存在一定的缺陷,若對模態(tài)分解算法進行改進,可以提高模型的預(yù)測精度。

c)對比模型四(CEEMDNA-GS-LSTM)與其他3個模型的預(yù)測結(jié)果,克服模態(tài)混疊和重構(gòu)誤差問題的CEEMDAN-GS-LSTM預(yù)測模型預(yù)測效果優(yōu)于EMD-GS-LSTM、EEMD-GS-LSTM預(yù)測模型。與模型二相比,其RMSE、MAE和MAPE分別降低了33.65%、27.57%與0.76%,對比模型三,其RMSE、MAE和MAPE分別降低了47.53%、55.99%與2.33%。采用CEEMDAN模態(tài)分解算法的預(yù)測模型均優(yōu)于EMD、EEMD,合適的模態(tài)分解算法可以更有效地追蹤徑流時間序列,進而提高模型預(yù)測精度。

綜上所述,在面對更具波動性的日徑流預(yù)報中,運用模態(tài)分解算法去預(yù)測日徑流可以提高模型預(yù)報精度,且在本文4個模型中,模型四(CEEMDAN-GS-LSTM)模型預(yù)報效果優(yōu)于其他模態(tài)分解算法結(jié)合LSTM的組合預(yù)測模型,CEEMDAN模態(tài)分解算法在日徑流預(yù)報將時間序列自適應(yīng)地將其分解為一系列不同尺度的模態(tài)分量,與EEMD、EMD方法相比,其能實現(xiàn)對幅值較大數(shù)據(jù)分解的完整性,有利于進一步提高預(yù)測精確度。

3 結(jié)論

本文針對日徑流預(yù)報構(gòu)建了基于CEEMDAN與GS-LSTM的組合預(yù)報模型,并對比分析了不同模態(tài)分解算法在日徑流預(yù)報中的影響,通過算例分析得到以下結(jié)論。

a)本文首先利用傳統(tǒng)模態(tài)分解算法EMD和EEMD對徑流數(shù)據(jù)進行分解,再分別對各個模態(tài)進行預(yù)測,最后對預(yù)測值進行重構(gòu)得到最終的預(yù)測值,其預(yù)報結(jié)果明顯優(yōu)于單一預(yù)報模型,EMD-GS-LSTM組合預(yù)測模型其RMSE、MAE和MAPE分別降低了625.32%、334.21%與2.12%,EEMD-GS-LSTM組合預(yù)測模型其RMSE、MAE和MAPE分別降低了440.15%、42.22%與0.55%,表明將模態(tài)分解算法運用日徑流預(yù)報,可有效提高預(yù)報精度。

b)本文所使用的克服模態(tài)混疊和重構(gòu)誤差CEEMDAN模態(tài)分解算法對日徑流進行分解,并使用網(wǎng)格搜索參數(shù)尋優(yōu)算法對LSTM模型的超參數(shù)進行優(yōu)化,其預(yù)報結(jié)果優(yōu)于EMD、EEMD模態(tài)分解算法,其預(yù)報效果較EMD-GS-LSTM組合預(yù)測模型,其RMSE、MAE和MAPE分別降低了236.45%、124.68%與0.76%,表明該模型可以更準確地對日徑流進行預(yù)報。

以“分解-預(yù)測-重構(gòu)”這個徑流預(yù)報方法而言,模態(tài)分解算法的選取將會明顯影響到日徑流的預(yù)報精度,在未來的研究中,改進模態(tài)分解算法,對徑流時間序列進行更加有效且精準的分解,這將獲得更加理想的預(yù)測效果。