嚴春華,張 玥,姜中清,梅 杰,高志發
(江蘇省水利勘測設計研究院有限公司,江蘇 揚州 225000)
徑流量預測是洪水預報的關鍵環節,對于水庫大壩、堤防、水閘等水利工程的安全穩定及正常運行有著至關重要的作用,關系到流域人民生命財產安全[1]。隨著人工神經網絡和深度學習的發展,基于數學模型的徑流量預測的研究成果也越來越多,如王軍等基于TCN-Attention模型對黃河徑流量進行預測[2],李代華構建了基于PCA-SHO-SVM和PCA-SHO-BP模型進行徑流預測[3]。這些組合模型可以較好地實現對徑流量的預測,但忽略了徑流量的時間序列的滯后性和誤差累積,從而導致誤差的增加和收斂速度的下降。因此,需要綜合考慮組合模型的特性,進一步深度研究算法的優缺點,進而構建高精度和穩定性的徑流量預測數學模型。
為了解決徑流量時間序列的滯后性,本文引入偏自相關函數進行研究[4]。偏自相關函數在處理時間序列數據時可以研究監測序列和本身的k階滯后序列兩者的相關性,從而得出徑流量序列對某時刻的影響作用程度。在智能學習算法中,目前以神經網絡模型(如CNN、RNN和LSTM等網絡)的預測能耐較好,其缺點為運算時間長,如果訓練和預測集過大的話,迭代速度會大大減慢、收斂時間長。為了解決和克服常用神經網絡算法的缺點,2014年Cho等[5]優化了LSTM算法模型,構建了GRU算法模型,借助設置隱藏層提高訓練學習的速度。GRU相比于LSTM網絡優化了網絡參數,計算收斂速度更快,并且兩個算法預測精度相差不大,所以本文選擇GRU網絡構建徑流量預測模型。為了進一步地提高徑流量預測模型的精度,需要解決預測過程中的誤差修正問題。由于徑流量的非穩定性,任何模型模擬出的徑流量均具有一定偏差,并且誤差的積累是一種非線性問題,而1998年Huang等[6]提出了一種可以較為高效地解決非線性與非平穩信號時頻分解難題的方法——經驗模態分解(EMD),因而采用EMD方法對殘差部分進行預處理。該處理時間序列的方法被廣泛運用到很多研究領域并有不錯的結果精度,如測繪空間定位[7]、風力預測[8]、軍事雷達[9]、工程爆破[10]研究等領域。因此,本文提出了以GRU為主要建模手段的徑流量預測方法,采用偏自相關方法確定時滯步長的前處理和EMD進行誤差修正的后處理,分別解決徑流量預測中存在的收斂速度慢、時序滯后性和誤差積累的問題。通過工程案例以及和其他預測模型進行對比,驗證了該徑流量預測模型的精度和可靠性。
為了尋取徑流量時間序列的滯后特性,采用偏自相關分析方法,并通過分析結果分析輸入量與預測值之間的關系[11]。偏自相關函數主要研究yt和yt+k之間的相關性,而排除了它們之間元素yt+1,yt+2,…,yt+k-1對yt的影響,其表達式如下[12]
(1)

GRU模型保留了LSTM模型的基本網絡構造,同時增加隱藏層,模型的基本架構由細胞狀態、更新門和重置門組成,該算法模型架構如圖1所示。

圖1 GRU網絡模型的架構
GRU算法的隱藏層的輸入由上一個時間節點的隱藏層狀態輸出h(t-1)和當前時間序列實際輸入x(t)組成。模型內部的更新門z(t)和重置門r(t)是優化時序預測模型的內部結果和減少參數數量的關鍵組成部分[13]。
更新門的作用是對當前時間節點輸入值的梯度判斷,與LSTM算法輸入門的功能基本一致,能夠極大程度克服梯度消失的難題,其門結構為
z(t)=σ(Wz·h(t-1)+Uzx(t)+bz)
(2)
式中,σ為激活函數,一般是sigmoid函數;Wz、Uz分別為更新門結構中的權重;bz為更新門結構中的偏置。
重置門的作用是判斷前一時間節點隱藏層狀態的h(t-1)對當前時間節點的x(t)的影響程度。如果程度不大,從而打開重置門,其門結構為
r(t)=σ(Wr·h(t-1)+Urx(t)+br)
(3)
式中,σ為激活函數,一般是sigmoid函數;Wr、Ur分別為重置門結構中的權重;br為重置門結構中的偏置。
確定更新隱藏層中的細胞狀態。細胞狀態的更新由更新門與重置門的輸出確定,其中隱藏層狀態中最右側的h~(t)為
h~(t)=tanh(W(r(t)°h(t-1))+Ux(t)+b)
(4)
式中,tanh為該部分的激活函數;W、U分別為更新h~(t)結構中的權重;b為更新h~(t)結構中的偏置。
當前時間節點隱藏層的輸出h(t),結合更新門的輸出更新細胞狀態得到當前時間節點隱藏層的輸出h(t)
h(t)=(1-z(t))°h(t-1)+z(t)°h~(t)
(5)
計算預測值的輸出是對當前時間節點的預測值輸出,則輸出門的輸出為
y(t)=σ(Vh(t)+c)
(6)
式中,V、c分別為權重參數及偏置。
EMD是一種自適應分解法,把具有時序特性的信號分解成不同頻率特征的本征模態函數(Intrinsic Mode Function,IMF)和一個長期趨勢項。EMD法根據原始信號的上包絡線和下包絡線的均值求解瞬時均值,進而實現對原始信號的特征分解[8]。給定原始數據信號xt,EMD法的實現基本分為5個步驟。
第1步,求解原始信號xt的上包絡線u(t)與下包絡線l(t)。計算xt中所有的局部極值點,再采用三次樣條插值法對極值點插值處理,從而得到原始信號的上下包絡線。
第2步,計算瞬時均值。把第1步得到的上包絡線u(t)和下包絡線l(t)取均值,即為瞬時均值m(t),
(7)
第3步,把原始信號xt減掉瞬時均值m(t),可得新序列f1(t),
f1(t)=xt-m(t)
(8)
第4步,得到IMF分量一定要有兩個條件:其一為序列f1(t)的局部極值點與過原點的個數的差值必須是0或1;其二為任何時刻,瞬時均值都為0。基于構成IMF的條件對序列f1(t)進行判別,如果達到IMF特征條件,序列f1(t)為IMF分量;如果序列f1(t)達不到條件,把f1(t)作為原始信號序列重復步驟1到3,一直分解得到的序列滿足IMF的特征條件。而且設置c1(t)=f1(t),得到信號xt和c1(t)的差值,從而一階殘差量為r1(t),再把殘差量作為原始信號重復步驟1到3。
第5步,當有第1個IMF分量f1(t),使c1(t)=f1(t),把原始信號xt減掉c1(t)可以得到1個殘差量r1(t),為
r1(t)=xt-c1(t)
(9)
上述5步完成了EMD法的首次信號分解。進一步地,把殘差量r1(t)重新重復上述步驟從而可以分解得到若干個IMF分量。當第n個殘差量rn(t)只有單個極值點或者為單調序列時分解停止,EMD分解流程結束。
原始信號xt經過分解可以得到若干IMF分量為
(10)
式中,fi(t)為分解的第i個IMF分量;rn(t)為第n次分解的殘差量;n為原始信號的分解次數,同時是分解得到的IMF分量數。
綜上理論,構建基于GRU的徑流量預測模型,該模型的基本實現流程如下:
第1步,確定徑流量序列的訓練學習集和預測集。
第2步,通過偏自相關性方法確定徑流量序列的時滯步長。
第3步,基于GRU算法進行學習訓練,輸出徑流量預測值。
第4步,提取殘差并基于EMD分解后利用GRU對各分量進行殘差修正。
第5步,并與其他常用算法對比,驗證本文模型預測結果的精度和模型的可靠性。
將某河流的徑流量采集數據作為本章方法的建模對象,研究時段為2015年1月1日至2019年12月31日,共計5 a,該時段的徑流讀取頻次為每日一次,其徑流量原始數據如圖2所示。在進行建模分析時,首先需要對徑流量數據進行數據集的劃分,本文以2015年1月1日~2018年12月31日4年的日徑流量數據作為訓練樣本集,2019年1月1日~2019年12月31日1年的日徑流量數據作為測試樣本,以分析提出的徑流量預測模型的有效性。

圖2 徑流量原始變化數據
由理論部分可知,對河流徑流量建立GRU模型進行徑流時序特性分析之前,首先應當確定時間序列的時滯步長,這對于正確建立歷史徑流量與當前徑流量之間的關系具有重要意義。若是確定的時滯步長過小,則不足以有效建立徑流量的時序關系;若是時滯步長過大,則將增加過多無用信息,GRU模型在學習徑流量時序關系時的效果也將變差。因此,本文按照理論部分提出的偏自相關分析方法進行時滯步長確定,若將徑流量的當前時刻確定為t,本案例中分析了t-1,t-2,t-3,…,t-12共計12個時間點處徑流量與當前時刻徑流量之間的偏自相關性,如圖3所示。

圖3 相關性分析
由圖3可以看出,當徑流量時間序列的時滯步長小于8時,歷史徑流數據與當前徑流量之間的相關性是逐漸下降的;當時滯步長的取值大于8時,歷史徑流數據與當前徑流量之間的相關性趨于平穩且處于較低相關程度。因此,在后面的時序建模中,將時滯步長的取值確定為8。
通過2.2節設定的徑流量時間序列時滯步長,可以對原始徑流量時間序列進行GRU建模。在此過程中,以2015年1月1日~2018年12月31日4年的日徑流量數據作為訓練樣本集,2019年1月1日~2019年12月31日1年的日徑流量數據作為測試樣本,即訓練集樣本數和測試集樣本數的比例為4∶1。首先通過訓練集樣本訓練GRU模型,然后將訓練完成的模型應用在測試集上,則徑流量GRU模型在測試集上的分析結果如圖4所示。

圖4 基于GRU模型的徑流量預測結果
通過圖4可以看出,該方法能夠大體模擬徑流量的變化趨勢,然而由于徑流量波動性較大,導致分析結果不能夠準確估計其波動性,因此,本文提出了殘差修正的方法對訓練集部分擬合殘差進行進一步分析,訓練集部分擬合結果如圖5所示,由此得到的徑流量殘差如圖6所示。

圖5 基于GRU模型的訓練集徑流量擬合結果
由理論部分可知,本小節主要針對GRU分析結果殘差進行修正,目的是得到更為準確的徑流量預測結果。首先采用EMD方法對上圖中的徑流量殘差進行分解,從而降低徑流量殘差序列的復雜性以提升建模分析結果的精度,EMD分析結果如圖7所示。

圖7 基于EMD的徑流量殘差分解結果
由圖7可以看出,徑流量殘差經過EMD處理后,原本具有很強不規則特性的殘差變成了具有不同變化頻率的殘差分量,這對于殘差建模分析提供了基礎。與直接對徑流量包含殘差的建模分析相比可知,通過對具有較低頻且規律較為明顯的分量進行建模,更容易學習其包含的徑流量信息,殘差修正的效果更明顯。通過對徑流量各殘差(IMF1~IMF11)分別建立GRU模型,并對各GRU模型的預測結果進行綜合得到最終的徑流量預測結果,如圖8所示。

圖8 基于殘差修正的徑流量預測結果
對比圖4和圖8徑流量測試階段的預測結果,很顯然在經過殘差修正后,徑流量的預測結果更加接近徑流量真實值,分析結果的波動性較之殘差修正之前也更接近真實值,驗證了本文提出的徑流量預測方法的有效性。
上文已分析得到基于偏自相關分析和殘差修正的徑流量預測模型的預測結果,為說明該方法在徑流預測領域的優越性和可靠性,選取ARIMA和LSTM這2種常用的時間關系分析模型進行對比說明。其中在建立上述2個模型時,徑流數據集的劃分同上述案例。3種模型的預測結果對比如圖9所示。

圖9 3種方法徑流量預測結果
由圖9可以看出,ARIMA的預測結果較為平滑,且與徑流量的真實值相差甚遠;而具有時序學習能力的LSTM方法相對來說具有一定模擬原始徑流量波動程度的能力,相比之下文中方法預測結果與真實的徑流量變化過程一致,徑流預測精度高。上述3種方法分析得到的徑流量預測結果殘差對比圖,如圖10所示。由圖10可以看出,ARIMA的預測誤差明顯是最大的,甚至在某些時間內其誤差高達3 000 m3/s,說明該模型并不適合應用于徑流量預測問題;LSTM相對誤差較小,而文中提出的方法徑流量預測結果誤差最小。3種方法徑流量預測結果評價指標如表1所示。從表1可以看出,文中方法分析得到的徑流量預測結果對應的RMSE值和MAE值是3種方法中最小的,其R2值也明顯優于另外2種方法。以上關于3種方法的對比分析說明了本文所提出方法的優越性,這為洪水預報提供了高可靠度的技術支撐。

表1 3種方法徑流量預測結果評價指標

圖10 徑流量預測結果殘差對比
本文通過偏自相關方法的前處理和EMD誤差的后處理,構建了基于GRU模型的徑流量預測模型,其結論如下:
(1)基于GRU模型的徑流量預測結果更接近于徑流量真實值,計算結果的波動性更接近真實值,預測誤差更小,說明該模型的工程實用性高。
(2)通過與ARIMA和LSTM預測模型的對比表明該模型的優越性,驗證了該徑流量預測模型的有效性和可靠度,對洪水預報提供了有力的技術支撐。
(3)由于本文采用的是固定時滯步長,而河流在不同時間下的步長是變化的,可以進一步研究動態滯后步長對預測值的影響。