基于鯨魚優化算法的長短期記憶模型水庫洪水預報*

丁藝鼎,蔣名亮,徐力剛,范宏翔,呂海深

(1:河海大學水文與水資源學院,南京 210098) (2:中國科學院南京地理與湖泊研究所,南京 210018) (3:中國科學院大學南京學院,南京 211135) (4:江西省鄱陽湖流域生態水利技術創新中心,南昌 330029)

洪水預報的目的是預測短、中、長期河道洪水的發生與變化趨勢。它為防汛搶險和防洪系統調度等決策提供依據,為水資源的合理利用和保護、水利工程的建設和管理以及工農業的安全生產服務。我國降水時空分布不均勻,加之全球氣候變暖導致的極端天氣事件,使得我國近年來洪水災害發生頻繁,這種洪水災害導致的巨額經濟損失不容小覷,總而言之精準及時的洪水預報方案顯得尤為重要。由于城市化的影響,流域中的下墊面情況較自然狀態發生了巨大變化,這種變化使得以經驗公式為基礎的概念模型不得不重新考慮下墊面變化對經驗公式的影響。其次,各種水利工程例如水壩等對自然的徑流過程也有著較大影響,如Nilsson等[1]指出全球大部分的大型河流水文情勢受到大壩建設的影響;而在我國已建庫容在10萬m3以上的水庫工程98112座,其中大型水庫744座,中型水庫3938座[2],所以當前形勢下考慮水利工程影響的洪水預報有著不言而喻的意義,精確的洪水預報可以為水庫群之間防洪減災和水資源調控提供實際支撐。

目前應用于洪水預報的水文模型主要有3種:(1)數學物理模型;(2)概念模型;(3)數據驅動模型(黑箱)[3-4]。隨著計算機技術的發展,人工智能神經網絡一家獨大,而深度學習又逐漸成為人工智能機器學習領域的主流技術,將深度學習理論引入水文預報已經成為當下洪水預報新的發展方向[5-8]。長短期記憶模型(LSTM)相較于傳統水文模型有著更強的非線性擬合和處理長期記憶的能力,在處理具有時序關系的問題上表現更加優秀[9-12],同時LSTM網絡可以自動學習洪水和降雨的特征,具有更高的自適應性和普適性。但是LSTM模型的參數較多,率定參數工作量大且傳統的率定方法主要依賴于經驗或者人工的方法,容易導致模型的參數落入局部最優解以及泛化能力不強等問題[13],如何提高數據驅動模型在洪水預報中的準確性、可解釋性,并增強對未來洪水預測的能力,是目前洪水預報領域的研究熱點。針對這類問題本文采用優化算法代替人機交互率定。鯨魚優化算法(whale optimization algorithms,WOA)有操作簡單、參數少以及跳出局部最優解能力強等優點,在求解優化問題中展現出優良的性能和廣泛的應用前景[14]。

因此,本文以浙江省橫錦水庫為研究對象并且選用新安江模型[15]為對比模型,基于1986-1997年洪水資料對該流域進行次洪模擬,引入優化算法和長短期記憶網絡構建浙江橫錦水庫洪水預報模型,通過置換特征重要性(permutation feature importance)、SHAP(Shapley additive explanations)方法確定特征權重,再結合傳統水文方法分析模擬結果的物理機制,以期提高神經網絡模型的預報可解釋性,為提高模型的外延性[16]提供一種新的思路。

1 數據與方法

1.1 研究區域及數據

1.1.1 研究區數據資料 橫錦水庫(29°0′～29°20′N,120°00′～121°00′E)位于我國東南沿海,屬于錢塘江水系,流域多山,土層較薄,并且水系發育良好。流域屬典型的亞熱帶季風氣候,受東亞季風影響,盛行東南風,降水有明顯的季節變化,雨量豐沛但時空分布不均,空氣濕潤,雨熱季節變化同步。橫錦流域控制斷面以上有兩條支流匯入,控制站以上集水面積為378 km2,河長81 km[15]。選取該流域中7個雨量站(史姆、西塢、銅錢、窈村、龍壇、八達、橫錦)1986-1997年的日蒸發量以及日降雨量和1 h時段降雨量。1個斷面流量監測站(橫錦站)的日流量和時段流量資料,時段長1 h,流域雨量站、流量站分布見圖1。本文以1986-1993作為模型的率定期,1994-1997年作為模型的檢驗期,具體率定期和檢驗期次洪資料的選擇見圖2。

圖2 次洪資料率定期和檢驗期的劃分Fig.2 Division of regular and test period in sub-flood data rate

1.1.2 數據輸入與輸出 研究中將水文站的上3個時刻時段流量值、6個雨量站和一個水文站的前3個時段和當前時段的雨量資料作為7個特征值輸入模型作為輸入數據,后第N時段的流量作為輸出值,N為預見期,本文選預見期為1 h。

1.2 長短期記憶模型計算原理

LSTM是一種對于傳統的循環神經網絡(recurrent neural network,RNN)的改進,在傳統的RNN網絡中信息由于梯度消失[17],無法學習到數據中的依賴關系。在LSTM模型中使用記憶細胞替換原先RNN中的細胞單元,從而能有效地提取出時間周期的序列數據特征值[18]。LSTM隱藏層中記憶細胞單元結構如圖3所示。

圖3 LSTM隱藏層中記憶細胞單元結構Fig.3 The structure of memory cell unit in the hidden layer of LSTM

LSTM區別于其他神經網絡之處在于其擁有獨特的單元狀態(細胞狀態),Ct表示神經網絡對于t+1時刻前所有輸入信息的總結。輸入信息經過記憶細胞單元具體過程如下:

輸入信息通過遺忘門,遺忘門決定了記憶細胞單元從上一狀態中遺忘多少信息,其計算式為:

ft=σ[Wf(ht-1,xt)+bf]

(1)

式中,Wf為遺忘門的權重矩陣;bf為遺忘門的偏移向量;σ為sigmoid激活函數,ft經過sigmoid函數運算之后就被約束在(0,1),1表示“全部記憶”,0表示“全部遺忘”。

經過輸入門,輸入門決定哪些信息被記憶細胞單元添加到記憶中,其計算式為:

it=σ[Wf(ht-1,xt)+bi]

(2)

(3)

決定出遺忘和記憶的信息后,對記憶單元狀態變量進行更新計算,其計算式為:

(4)

式中,Ct為更新后的記憶單元狀態變量,⊙是點乘運算。

處理后的數據最后通過輸出門Ot,輸出門的計算決定了隱藏層狀態變量ht,其計算式為:

Ot=σ[W0(ht-1,xt)+b0]

(5)

ht=Ot⊙tanh(Ct)

(6)

式中,W0為輸出門權重矩陣,b0為輸出門偏移量。

1.3 鯨魚優化算法及其改進

WOA算法是受到自然界座頭鯨捕食行為啟發的一種智能優化算法,該算法較傳統優化算法有精準度高且收斂速度快的特點。WOA算法的基本工作流程為:(Ⅰ)生成[0,1]的隨機數p,當p>0.5時采用泡泡網攻擊方式,否則當|A|<1時進入獵物包圍階段,|A|≥1時進入搜索獵物階段;(Ⅱ)隨著迭代次數不斷增加,a逐漸由2減小到0,|A|也逐漸收斂到0,算法從搜索獵物階段過渡到包圍獵物階段[19]。但是由于WOA算法易陷入局部最優并且模型內參數調整方法簡單從而使得全局探索和局部開發能力難以平衡[20]。本文使用基于混沌的正余弦策略來改進鯨魚優化算法,具體改進步驟如下:(1)采用在WOA算法中尋優能力和收斂速度較好的Tent映射來生成初始種群[21];(2)使用正余弦算法(Sine cosine algorithm,SCA)作為全局優化算法篩選首領鯨魚位置;(3)添加基于混沌理論的自適應權重,使算法的全局搜索和局部開發得到較好的平衡,一定程度上也增加了算法跳出局部最優區域的概率。該策略增加模型跳出局部最優的能力,減小了原有模型缺陷[21]。其改進的自適應權重計算式為:

(7)

X(t+1)=ω·X*(t)-A·D

(8)

X(t+1)=ω·Xrand-A·Drand

(9)

X(t+1)=D′·ebl·cos(2πl)+ωX*(t)

(10)

式中,Tmax表示最大迭代次數;ωs表示慣性權重初始下限;ωe是最大迭代次數的慣性權重上限。研究表明,當ω在[0.4,0.9]之間變化時,算法具有較好的動態尋優能力。通常較大的慣性權重能夠讓算法具有較好的全局探索能力,而較小的慣性權重會使算法具有較好的局部開發能力,隨著迭代的變化,將慣性權重從0.9非線性遞減至0.4,表明模型實現以全局探索為主向局部探索能力為主的動態轉變[22],本文所提改進的WOA算法的流程如圖4所示。

圖4 改進算法流程Fig.4 The flow of the improved algorithm

1.4 WOA-LSTM預報模型的建立

1.4.1 模型結構 本文選擇使用TensorFlow-GPU 2x深度學習框架來搭建網絡結構,通過高階應用程序接口Keras來進行深度學習模型的設計、調試、評估、應用。構建網絡主要步驟為:

1)數據預處理:本研究使用的是流量和降雨資料,由于數據的類型不同,數據差異較大,在輸入模型之前需要將輸入數據資料進行標準化。由于降雨和流量數據會時而出現異常值(極端事件),而異常值對極差標準化方法影響巨大,故使用標準差標準化公式(11)。對訓練數據和測試數據分別進行標準化。

(11)

2)LSTM優化器選擇:Bi-LSTM的激活函數選擇ReLU,全連接層的激活函數為linear。并且優化器選擇Adam算法。

3)損失函數的選擇:深度學習中可供選擇的損失函數有多種,本次研究最終選擇了較為常用的均方誤差損失(mean squared error, MSE)[23],計算公式為:

(12)

式中,y為實測值;yper為預測值;n為資料序列長度。

4)數據校準計算:在輸出結果通過馬氏檢驗后,使用馬爾科夫鏈對最后模型輸出的數據進行殘差校正。

1.4.2 模型參數選取 本研究中模型的神經網絡層數、學習率(learning-rate)、迭代輪次(epochs)、單次訓練數據窗口大小(batch-size)、各層的節點數(units)由經驗確定范圍后由WOA優化算法迭代計算出。分別取Bi-LSTM為1、2、3、4層,在預見期為1 h時,分別獨立模擬10次計算的納什系數(NSE),模擬結果見圖5。本文選用模擬效果最好的3層Bi-LSTM構建LSTM模型。

圖5 不同Bi-LSTM層數模擬結果對比Fig.5 Comparison of simulation results with different numbers of Bi-LSTM layers

1.4.3 模型評價指標 預報計算過程與實測過程之間的吻合程度可用確定性系數NSE來評定。NSE是一個標準化統計指標,其值介于(-∞,1)之間,等于1表示模型結果完美擬合實測值。計算公式為:

(13)

另外,還選擇洪峰相對誤差、徑流深相對誤差和峰現時差3個指標來評定預報精度。根據《水文情報預報規范》(GB/T 22482-2008)來最后評價預報方案精度。

1.5 神經網絡解釋方法

本文采用在Python3.x中使用置換特征方法(permutation feature importance method,PFIM)和SHAP法來計算特征值的重要性,置換特征方法計算步驟為:(1)選擇一個特征值,并且在數據集上對該特征值的所有數據進行隨機排列。(2)計算新的模型預測結果。(3)對比新舊模擬結果的差異,差異越大說明該特征值越重要。

SHAP是一種用于解釋機器學習模型預測結果的方法,它基于Shapley值的概念,給出了一種針對單個特征或特征組合的影響分析方法。在SHAP中,將每個特征視為博弈中的參與者,并將模型的預測視為博弈中的收益。通過使用多個樣本來計算每個特征的Shapley值,SHAP提供了每個特征對模型預測的影響的全面解釋。SHAP方法可以解釋機器學習模型中每個特征的相對重要性,即在模型預測中每個特征對結果的貢獻大小。

2 結果與討論

2.1 預報精度總體效果評價

2.1.1 模型率定期參數設置及模擬結果 文中的LSTM模型和WOA-LSTM模型的參數率定結果見表1。LSTM的參數使用人機交互率定,而WOA-LSTM參數使用自動率定,選用NSE作為率定的適配度。為對比研究WOA-LSTM的模擬計算效果,對比模型選用在該地區運用較好的洪水預報模型——新安江模型[24],通過人機交互率定法得到新安江模型參數結果見表2。

表1 神經網絡參數率定結果Tab.1 The results of parameter tuning for neural network

表2 新安江最終參數Tab.2 The final parameters of the Xin’an River

圖6展示了率定期18場洪水的模擬結果,并采用NSE、均方根誤差(RMSE)以及洪峰相對誤差(絕對值)等指標評估了模型在率定期的性能表現。計算結果表明,在3個指標下,WOA-LSTM模型的性能表現均顯著優于其他兩個模型。具體而言,在率定期內,WOA-LSTM模型的RMSE值為7.23 m3/s,NSE值高達0.997,這表明該模型不僅能夠準確地模擬洪峰流量,而且對洪水形狀的刻畫也非常精準。

圖6 率定期洪水模擬對比Fig.6 Comparison of rate-constant flood simulation

2.1.2 模型檢驗期模擬結果 使用模型對檢驗期3年(1994-1997年)資料的7場洪水進行模擬計算。不同模型模擬計算結果見圖7,統計得到在檢驗期中LSTM、WOA-LSTM和XAJ 3個模型的NSE值分別為0.946、0.926 和0.858,從NSE數值來看,神經網絡模型能更好地表達出實測值和模擬流量之間的關系。圖7b顯示了實測數據和模擬值的頻率累積分布,4條線幾乎一致表明3個模型都可以較好地模擬出該流域的洪水過程,但是由于模型本身和資料誤差會導致模型在模擬過程中出現異常值。圖7a得出LSTM和WOA-LSTM在模擬過程中都出現了異常值,盡管LSTM模型能夠更好地擬合小流量數據,在大流量模擬時表現欠佳。初步分析表明,這可能是由于模型結構復雜度和參數未達到最優配置,從而導致洪水預測偏大。為了進一步驗證這一結論,表3統計了檢驗期7次洪水事件的洪峰流量相對誤差、NSE、徑流深相對誤差以及峰現時差。

表3 目標函數計算Tab.3 Objective function

圖7 檢驗期7場洪水模擬結果對比Fig.7 The comparison of simulation results for seven flood events during the verification period

根據《水文情報規范》(GB/T 22482-2008),3個模型模擬的7場洪水,WOA-LSTM和XAJ模擬結果洪峰相對誤差、徑流深相對誤差和峰現時差全部小于誤差許可,除LSTM模擬的31950426號洪水峰現時差大于誤差許可3 h外,其余預測洪水全部達標。洪峰相對誤差的絕對值全部小于誤差許可20%,但是LSTM模型預測的7場洪水中有6場洪水的預報洪峰大于實測值,存在系統誤差,相對誤差最大值為17.58,接近誤差許可,這與圖7a得到的結論相吻合,說明LSTM在預測大流量時確實存在過擬合問題。而本文所使用的WOA-LSTM模型中由于優化的參數選取和MC矯正流量殘差,雖然小幅度增大了洪水過程線的波動而減小了NSE值,但是卻極大地減小了LSTM自身模擬洪峰偏大的誤差。從模擬效果來看,WOA-LSTM的NSE平均值為0.915、LSTM的NSE平均值為0.928,而洪峰流量誤差絕對均值卻由LSTM模擬的8%下降到4.9%,降幅約為40%,可見鯨魚優化算法、長短期記憶算法和馬爾科夫鏈結合使用在保證流量預測精度的同時,避免了預測大流量時存在的過擬合問題,提升了洪峰流量預測精度。綜上所述,WOA-LSTM擁有更好的模擬效果。

2.2 典型洪水模擬分析

當相同的模型應用于不同的洪水過程時,其表現往往取決于降雨情況的復雜性。多峰洪水相比于單峰洪水,其洪峰流量的大小、時間、洪水過程等方面都更為復雜,而大流量洪峰由于短時間內迅速形成,數據采集和傳輸過程中可能存在延遲或遺漏,兩種情況都會導致預報的準確性下降,本節對比分析了3個模型對檢驗期這兩種典型洪水的模擬結果,分別為多次降雨造成典型多峰洪水(31940609、31950426)和單次強降雨帶來的典型單峰洪水(31970818)。

從3次典型洪水的流量過程線(圖8)可以看出,XAJ模型在模擬復雜降雨過程帶來的多峰洪水時,往往誤差較大,尤其體現在洪峰流量和退水過程。而LSTM和WOA-LSTM模型卻能較好地模擬多峰洪水的復雜過程。而在大洪水預報中,由于歷史大洪水資料少,長短期模型能學習到的特征也會變少,從而影響到預測精度,在本文中加入了馬爾科夫鏈殘差矯正步驟可以在一定程度上減少由于學習資料少而帶來的精度下降問題。圖8顯示,31970818號洪水在檢驗期中洪峰流量最大,LSTM的模擬結果無論是洪峰流量還是峰現時間都有較明顯的偏差,而這一偏差在WOA-LSTM中則變得微乎其微,說明這種改進策略在橫錦水庫的洪水預報中是可行的。

圖8 典型洪水流量過程線Fig.8 Typical flood flow process line

2.3 不同預見期下模擬結果與比較

本節分析了3種模型在不考慮未來降雨情況下的預報能力,具體方法是從預報時刻開始,不考慮未來降雨資料,僅利用現有的實測降雨資料和流量資料作為模型輸入,每1 h驅動一次模型,對洪水過程進行逐時段預報分析。

表4表明,3種模型在預見期增加后預報能力都有所下降,當預見期為4 h時,WOA-LSTM、LSTM和XAJ對率定期洪水的綜合NSE分別為0.88、0.84和 0.74,WOA-LSTM模型有5場NSE優于其他兩個模型,平均洪峰相對誤差XAJ模型較低,但是總體相差不大。當預見期為6 h時,WOA-LSTM和LSTM的率定期綜合NSE分別為0.72和0.58,WOA-LSTM模型要遠好于LSTM,對于其他3個指標二者相差不大。當預見期為8 h時,WOA-LSTM和LSTM的率定期綜合NSE分別為0.55和0.34,峰現時差分別為3和5 h,平均洪峰相對誤差分別為22.9%和27.0%,平均徑流深相對誤差為21%和22%。在這4個指標中,WOA-LSTM模型的表現都要優于LSTM模型。雖然隨著預見期增加WOA-LSTM模型預報精度有所下降,但是下降幅度要小于LSTM?？偠灾?WOA-LSTM模型在預見期2～8 h內模擬結果都要好于其他兩個模型。

表4 不同預見期預報結果比較Tab.4 Comparison of forecast results at different lead times

從不同預見期預測結果分析中可以看出神經網絡模型的優勢在于,即使沒有實時降雨資料輸入,它仍能夠有效地考慮未來降雨對系統的影響。通過利用過去的降雨數據,神經網絡能夠推斷未來降雨的潛在模式和趨勢,從而提供對系統響應和預測的有價值信息。因此,該模型不僅在實時數據可獲得時表現出色,而且在沒有實時降雨資料時仍能提供有意義的預測和決策支持。例如在預見期6 h時,XAJ模型模擬的31940820號洪水洪峰相對誤差高達48.4%,而其他兩個模型模擬的洪峰相對誤差在6%以內,仍然可以較好地刻畫出洪峰流量,說明較傳統水文模型在長預見期且無未來降水資料輸入時表現不佳,神經網絡模型依舊可以預測到一定的洪水特征。但是當預見期增加時,數據關聯性減小,可能會使得模型學到無用特征而產生過擬合現象,導致模型預報精度下降[25]。在預見期為8 h時,WOA-LSTM和LSTM產生的loss曲線見圖9。在LSTM模型中隨著迭代輪次增加,測試集的損失逐漸增大,證明在關聯度差的數據集中LSTM模型確實會存在過擬合現象,而這種現象在WOA-LSTM模型訓練時卻并沒有這么明顯。證明由于WOA優化算法,進一步優化了神經網絡結構,提高了模型的穩定性,降低了過擬合的風險。盡管LSTM在納什效率方面表現更好,但在其他方面卻不如WOA-LSTM,WOA-LSTM在綜合效果上更出色。

圖9 模型擬合損失曲線Fig.9 The model fitting loss curve

綜上所述,WOA-LSTM在無需輸入降雨資料的較長預見期內,可以一定程度上避免過擬合,從而保持良好的預測精度和魯棒性。

2.4 WOA-LSTM可解釋性的討論

神經網絡模型是典型的黑箱子模型,由于本身不是構建于物理規律基礎上的,所以模型可解釋性較差,探究其參數和模型構造背后的物理意義逐漸成為神經網絡模型新的發展方向[26]。

本文先使用PFIM法和SHAP法分別計算各特征值的重要性,再結合水文模型和氣象學中使用較多的泰森多邊形法(thiessen polygon method,TPM)和實際流域特性驗證得到的雨量站面積權重的合理性。

由表5泰森多邊形計算結果表明橫錦流域各雨量站重要性相似,結合圖1分析,說明研究流域中各雨量站分布均勻,面積權重應該相似。圖10對比了不同方法得到的特征值重要性值,除了SHAP法給橫錦站賦予了較低的權重外,SHAP、泰森多邊形、PFMI所計算的權重值分布大致相同,表明WOA-LSTM模型可以通過數據資料學習到流域中一定的雨量站面積權重特征,而這種面積權重是符合現實地理分布的。

表5 不同方法的特征值權重Tab.5 Eigenvalue weights of different methods

圖10 站點面積權重與特征值的重要性評估Fig.10 Importance evaluation of the station area weight and characteristic value

2.5 結語

本文以橫錦流域作為研究區域,結合改進的鯨魚優化算法,馬爾科夫鏈和長短期記憶神經網絡共同構建洪水預報模型,討論了在不同預見期下模型的預測性能、結合置換特征,SHAP方法增加了模型的可解釋性,得到結果如下:

1)在對場次洪水預報中,當預見期同為1 h時,檢驗期1994-1997年7場洪水的NSE值為0.926,徑流深相對誤差、洪峰相對誤差和峰現時差全部在誤差許可以內,且優于對比模型。

2)WOA-LSTM模型中WOA優化算法通過改變神經網絡結構,從而降低模型在復雜數據集上的過擬合風險。

3)通過泰森多邊形,SHAP和置換特征值方法相互驗證,表明WOA-LSTM模型可以在訓練模型時學習到一定的物理機制,增加模型的可解釋性。

綜上所述,通過將優化算法和神經網絡相結合,在一定程度上避免了人工率定超參數陷入局部解、過擬合等問題,并且在水庫洪水預報過程中取得較好的成果,同時通過計算機方法和水文方法驗證了神經網絡在訓練模型時可以學到一定的物理特征作為數據預測的基礎,表明在以后的研究中可以結合置換特征值等一些新方法進一步增加神經網絡的可解釋性,增加模型計算結果的說服力。但是本文對于神經網絡物理機制的研究還不夠深入,今后可以通過更多的計算機方法,從神經細胞層面剖析模型的物理機制,增加模型的泛化能力,也可以結合水文統計等方法進一步提高模型預報結果精度。