基于隨機森林的水庫防洪調(diào)度研究

刁艷芳，王 蒙，2，王 昊，叢方杰，王 剛

（1.山東農(nóng)業(yè)大學水利土木工程學院，山東泰安 271018；2.華北水利水電大學水利學院，鄭州 450046；3.煙臺大學網(wǎng)絡與教育技術中心，山東煙臺 264005）

0 引 言

近幾年，極端氣候所導致的洪水災害頻繁發(fā)生，造成了巨大的社會、經(jīng)濟損失。水庫是保障防洪安全、提高洪水資源利用率和促進水生態(tài)文明建設的重要工程措施，科學合理的水庫防洪調(diào)度更能有效促進其作用的發(fā)揮。目前，水庫防洪調(diào)度的研究方法大體分為常規(guī)方法和系統(tǒng)分析方法兩類。常規(guī)調(diào)度方法是借助水庫防洪調(diào)度圖、抗洪能力圖等經(jīng)驗圖表實施操作的一種半經(jīng)驗、半理論方法［1，2］。隨著系統(tǒng)工程的發(fā)展和應用，許多學者嘗試采用系統(tǒng)分析方法解決水庫調(diào)度問題。早期主要是線性規(guī)劃［3，4］、動態(tài)規(guī)劃［5，6］、逐次優(yōu)化算法［7］和大系統(tǒng)分解協(xié)調(diào)原理［8］等；然而面對由水庫、湖泊、閘壩、河道及堤防、蓄滯洪區(qū)以及各種防洪工程組成的復雜防洪系統(tǒng)，基于生物學、物理學和人工智能的現(xiàn)代啟發(fā)式智能優(yōu)化算法迅速發(fā)展，其具有全局最優(yōu)化、魯棒性強、通用性強及高效性等優(yōu)點，主要包括遺傳算法［9］、蟻群算法［10］、粒子群算法［11］和人工神經(jīng)網(wǎng)絡［12］等。

然而無論是常規(guī)方法還是系統(tǒng)分析方法，均存在對歷史調(diào)度數(shù)據(jù)信息利用不夠的問題，而這些數(shù)據(jù)既包括了水庫流域水文特征與規(guī)律，又包括了決策者多年調(diào)度和配置經(jīng)驗。在當前水利部門已積累了大量流域水文、地貌以及洪水調(diào)度數(shù)據(jù)的基礎上，采用數(shù)據(jù)挖掘技術解決上述問題是主要手段。隨機森林算法是數(shù)據(jù)挖掘技術的一種，在降雨預測［13］、徑流預測［14，15］、水資源評價［16］以及洪災預測［17］等諸多領域得到了廣泛應用，但在水庫調(diào)度中應用極少。大量研究成果表明，隨機森林算法具有不易出現(xiàn)過擬合、泛化誤差低及預測精度高等優(yōu)點，被認為是當前最好的算法之一。因此，本文提出采用隨機森林算法擬定水庫防洪調(diào)度方案，并以山東省田莊水庫為例驗證了本方法的合理性。

1 研究區(qū)介紹

田莊水庫位于淮河流域沂河上游，流域面積424 km2，總庫容為1.305 7 億m3，興利庫容為0.684 億m3，是一座以防洪、灌溉為主，結(jié)合水力發(fā)電、水面養(yǎng)殖、工業(yè)供水等綜合利用的大（II）型水庫。水庫興建于1958年，1960年建成蓄水，按百年一遇設計，五千年一遇校核，其設計洪水位、校核洪水位、汛限水位、防洪高水位（P=2%）及死水位分別為312.38、315.07、310.64、312.33及293.64 m。水庫保護下游沂源縣城及南麻、悅莊、魯村等7 個鄉(xiāng)鎮(zhèn)30 萬人口及10 萬畝耕地。下游第一安全泄量為1 000 m3/s，田莊水庫控制下泄流量為600 m3/s，相應標準為20年一遇；下游第二安全泄量為2 000 m3/s，田莊水庫控制下泄流量為1 000 m3/s，相應標準為50年一遇。

田莊水庫流域地處中緯度，受東亞季風和歐亞大陸的影響，屬溫帶季風大陸性氣候，四季溫差大。流域多年平均降水量為730.6 mm，年際之間變化大，年內(nèi)分配也不均勻，主要集中在6-9月份，約占多年平均年降水量的74.9%。

2 資料與研究方法

2.1 降雨與洪水資料

根據(jù)田莊水庫1963-2019年的洪水調(diào)度數(shù)據(jù)，考慮到較大量級的洪水對防洪安全更具重要意義，選取了啟用溢洪閘且最大泄流量超過50 m3/s 的18 場洪水過程，整編出洪水過程逐時段的入庫流量、庫水位、蓄水量及出庫流量等數(shù)據(jù)資料，時段長為1小時，其洪號見表1所示。根據(jù)《水文情報預報規(guī)范》（GB/T 22482-2008），按照洪水總量W將18 場洪水劃分為小（W＜3 369 萬m3）、中（3 369 萬m3≤W＜6 470 萬m3）、大（6 470 萬m3≤W＜8 800 萬m3）、特大（W≥8 800 萬m3）4 個級別，列于表1 中。田莊水庫流域內(nèi)洪水主要由暴雨形成，由于天氣系統(tǒng)原因?qū)е卤┯瓿掷m(xù)時間一般較短；此外，田莊水庫流域為純山丘地區(qū)，山高坡陡，匯流歷時較短，因此，造成大部分洪水過程為形狀比較尖瘦的單峰過程，且洪水持續(xù)時間為3日左右。表1 所示的18場洪水過程均為符合上述特性的單峰洪水過程。田莊水庫流域內(nèi)共有6個雨量站，分別為包家莊站、草埠站、大張莊站、田莊水庫站、徐家莊站及朱家莊站，整理出與18 場洪水過程相對應的各雨量站逐時段的降雨量資料，鑒于雨量站分布較均勻，故采用算術平均法計算流域平均降水量。

表1 場次洪水及級別表Tab.1 Floods and their grades

2.2 研究方法

2.2.1 隨機森林算法

隨機森林（Random Forest，RF）算法是Breiman 在2001年提出的一種分類和預測的機器學習算法［18］，其以Bagging 集成學習算法和隨機空間算法為基礎。RF 算法分為RF 分類和RF 回歸兩類，本文采用后者進行水庫調(diào)度方案的擬定。與目前其他機器學習算法相比，RF 算法具有如下優(yōu)點：①較強非線性模擬能力，可以高效處理多變量和大數(shù)據(jù)量的問題，且無需進行變量的刪減和量綱統(tǒng)一處理，同時能夠評估所有變量的重要程度；②模型參數(shù)少、運算效率高、數(shù)據(jù)挖掘能力強及預測精度高等；③不易出現(xiàn)過擬合現(xiàn)象，對異常值、缺失值、干擾值的容忍度高，對數(shù)據(jù)集特征的挖掘具有很好的魯棒性。

（1）RF 回歸算法。RF 回歸算法的基本組成單元是回歸樹{t(x，θi)，i= 1，2，…，k}，其中x表示樣本數(shù)據(jù)集，k表示回歸樹的個數(shù)，θi表示第i棵回歸樹的參數(shù)向量，算法實現(xiàn)過程如下，具體流程見圖1所示。

圖1 RF回歸算法流程圖Fig.1 Flow chart of RF regression algorithm

①通過自助法（Bootstrap）重抽樣技術，從原始訓練樣本集N中有放回地重復隨機抽取n個樣本生成子訓練集k個，每個子訓練集構建1 棵回歸樹，每次未被抽到的樣本稱為袋外數(shù)據(jù)（Out-Of-Bag，OOB）。

②設原始訓練樣本集的變量個數(shù)為m，在每棵樹的葉節(jié)點處等概率隨機抽取mt個變量作為備選分枝變量，根據(jù)分枝優(yōu)度準則計算最佳分枝，且在建樹過程中mt保持不變。

③設葉節(jié)點的最小尺寸為nodesize，每棵樹最大限度地生長，不進行干預和剪裁。

④k棵回歸樹組成隨機森林回歸模型，將驗證樣本集代入該模型，以每棵回歸樹輸出結(jié)果的算數(shù)平均值作為預測結(jié)果。

（2）RF 回歸算法的參數(shù)優(yōu)化。由RF 回歸算法的計算過程可以看出，回歸樹個數(shù)k、節(jié)點備選分枝變量mt以及葉節(jié)點最小尺寸nodesize是3 個主要參數(shù)，其取值影響該算法的效率和精度。已有研究表明，k太小會使訓練不夠充分，降低RF 的隨機性和精度；太大會增加模型的計算量，導致計算效率低下、計算精度不高。mt太小會導致過擬合現(xiàn)象，預測分類的誤差變大及預測精度降低；mt太大會使建模效率降低。nodesize對RF 回歸算法性能的影響不顯著。由此可以看出，RF回歸算法的預測精度主要取決于k和mt，其取值采用OOB 均方誤差進行評價，計算公式為：

式中：MSEOOB為OOB 均方誤差；yi為OOB 中因變量的實際值；為OOB的預測值。

MSEOOB可以用來評估回歸樹的泛化誤差，Breiman 通過大量實驗數(shù)據(jù)證明OOB誤差是一種無偏估計［19］，不需要使用交叉檢驗進行估計。理論上，可以通過MSEOOB來選取最佳的k和mt數(shù)值組合；但是在具體實踐中，往往是設定k值或mt值，之后對另一個參數(shù)進行適應性調(diào)整。

2.2.2 隨機森林算法擬定調(diào)度方案的步驟

采用RF回歸算法推求場次洪水調(diào)度方案的步驟如下：

（1）洪水優(yōu)化調(diào)度數(shù)據(jù)集的生成。由于歷史洪水實測調(diào)度未必是最優(yōu)調(diào)度方案，故采用優(yōu)化算法對歷史洪水過程進行優(yōu)化調(diào)度，生成水庫洪水優(yōu)化調(diào)度數(shù)據(jù)集。

（2）水庫泄流量主要影響因素的選取及RF 回歸算法最優(yōu)參數(shù)的推求。通過查閱已有研究成果［20，21］可知，影響水庫泄流量的因素包括降雨量、凈雨量、庫水位、入庫流量及入庫水量等，整理出各場次洪水優(yōu)化調(diào)度過程逐時段的累計降雨量、累計凈雨量、庫水位、入庫流量、累計入庫水量及泄流量等數(shù)據(jù)，以泄流量作為決策屬性，以影響因素作為條件屬性，輸入RF 回歸算法中，由此計算出最優(yōu)的k和mt值。同時，RF 回歸算法可以計算出每個影響因素對每棵回歸樹的貢獻，取均值后即為影響因素對RF 回歸模型的貢獻程度，由此可以判斷出影響泄流量的主要因素。

（3）生成泄流量回歸樹。將所有場次洪水劃分為訓練樣本和驗證樣本，將訓練樣本洪水優(yōu)化調(diào)度過程逐時段的主要影響因素和泄流量輸入RF回歸算法生成泄流量回歸樹。

（4）由RF回歸算法擬定泄流量。

①設i=1；

②將第i時段的主要影響因素代入泄流量回歸樹中，計算泄流量qi，并由水量平衡方程求得zi；

③為保護下游防洪安全，判斷qi是否不大于該場洪水由水庫調(diào)度規(guī)則調(diào)洪后的最大泄流量q*。如果是，轉(zhuǎn)入第④步；否則，qi=q*；

④為確保水庫泄流能力約束，判斷qi是否不大于水位zi時的水庫泄流能力q（zi）。如果是，轉(zhuǎn)入第⑤步；否則，qi=q（zi）；

⑤令i=i+1，如果i≤T，T為洪水總時段數(shù)，轉(zhuǎn)入第②步；否則，調(diào)洪結(jié)束。

其流程見圖2所示。

圖2 RF回歸算法擬定調(diào)度方案流程圖Fig.2 Flow chart of operation scheme maked by RF regression algorithm

3 結(jié)果分析

3.1 RF回歸算法最優(yōu)參數(shù)的確定

田莊水庫為大（II）型水庫，調(diào)節(jié)性能較高且承擔下游防洪任務，故洪量在調(diào)洪過程中起主要控制作用。結(jié)合田莊水庫特性和已有研究成果，選取累計凈雨量、累計入庫水量、起調(diào)水位、實時水位以及入庫流量作為泄流量的影響因素。以田莊水庫下泄流量最小為目標函數(shù)，采用粒子群算法對18場洪水進行優(yōu)化調(diào)度，得到水庫洪水優(yōu)化調(diào)度數(shù)據(jù)集。以18場洪水優(yōu)化調(diào)度過程逐時段的累計凈雨量、累計入庫水量、起調(diào)水位、實時水位、入庫流量以及泄流量作為RF 回歸算法的輸入數(shù)據(jù)，回歸樹數(shù)量變化范圍為100～1 000，計算MSEOOB繪制于圖3 中。由圖3看出，當k=700 時MSEOOB最小，為32.6；當k=100 時MSEOOB最大，為37.3，兩者相差4.7，故確定最優(yōu)回歸樹的數(shù)量k=700。同時，由RF 回歸算法計算出的5 個影響因素對泄流量的影響程度分別為累計凈雨量21.8%、累計入庫水量16.1%、起調(diào)水位21.5%、實時水位16.0%、入庫流量24.6%，影響程度最大的為實時水位，最小的為累計入庫水量，兩者相差8.6%，差距較小，故將上述5個因素均作為RF回歸算法的輸入變量，即mt=5。

圖3 MSEOOB隨回歸樹數(shù)量的變化曲線Fig.3 The variation curve of MSEOOB with the number of regression trees

3.2 調(diào)洪結(jié)果分析

由表1 看出，小、中、特大洪水的數(shù)量分別為11 場、6 場和1場，故發(fā)生不同級別洪水的概率由大到小依次為小洪水、中洪水、特大洪水；洪水級別越大，對水庫及其上下游造成的威脅越大，因此，選擇洪號為19640912 與20010804 的兩場中洪水作為驗證樣本，其余16 場洪水作為訓練樣本，將訓練樣本逐時段的5個主要影響因素和泄流量作為輸入，由RF 回歸算法生成泄流量回歸樹，以此求解驗證樣本的泄流量。依據(jù)圖2的計算流程，求得19640912 與20010804 兩場洪水的泄流量及庫水位，與實測入庫流量、出庫流量及庫水位一并繪于圖4、5 中。本文比較了兩場洪水實測調(diào)度和RF 回歸算法調(diào)度結(jié)果的最大泄流量和最高水位兩個指標，計算了兩者的絕度偏差和相對偏差，見表2。由表2 可知：①兩場洪水最大泄流量的實測值均大于RF 回歸算法的計算值，19640912 洪水的絕度偏差和相對偏差的絕對值均小于20010804 洪水，由此可見，從保證下游安全而言，RF回歸算法的調(diào)度結(jié)果優(yōu)于實測調(diào)度；②兩場洪水最高水位的實測值均大于RF 回歸算法的計算值，19640912 洪水的絕度偏差和相對偏差的絕對值均大于20010804 洪水，由此可見，從保證水庫安全而言，RF回歸算法的調(diào)度結(jié)果同樣優(yōu)于實測調(diào)度。

表2 RF回歸算法計算值與實測值對比表Tab.2 Comparison between measured values and calculated values by RF regression algorithm

圖4 19640912洪水過程線Fig.4 19640912 flood hydrographs

圖5 20010804洪水過程線Fig.5 20010804 flood hydrographs

4 結(jié) 論

（1）本文首先對歷史洪水過程進行優(yōu)化調(diào)度生成優(yōu)化調(diào)度數(shù)據(jù)集，然后采用RF 回歸算法擬定水庫泄流量回歸樹，最后由泄流量回歸樹擬定水庫調(diào)度方案，指導實時防洪調(diào)度。該方法既能夠提取影響水庫泄流量的主要因素，又能夠充分挖掘歷史調(diào)度經(jīng)驗，開辟了水庫防洪調(diào)度的新方法。

（2）田莊水庫的實例可以看出，RF 回歸算法制定的水庫調(diào)度方案的最大泄流量、最高水位均小于實測值，表明RF 回歸算法擬定的調(diào)度方案較實測調(diào)度更能保證水庫本身及其下游防洪安全，故是一種合理的調(diào)度方案。鑒于田莊水庫暴雨洪水特性，本文僅對單峰洪水開展了研究，對于多峰洪水過程有待于進一步的驗證。 □