基于隨機森林回歸算法的山洪災害臨界雨量預估模型

趙 龍， 桑國慶， 武 瑋， 劉昌軍， 王君諾

(1. 濟南大學水利與環境學院， 山東濟南 250022； 2. 中國水利水電科學研究院， 北京 100038； 3. 水發規劃設計有限公司， 山東濟南 250100)

臨界雨量是山洪災害重要的預警預報指標[1]。 當某個區域在一段時間內降雨量達到某一量級時， 該區域會發生山洪災害， 這時的降雨量即稱為該區域的臨界雨量[2]。 隨著山洪災害防御精細化發展， 臨界雨量的研究對象逐漸從較大區域發展為沿河村落等居民聚居區[3]。 按照計算原理的不同， 可以將臨界雨量的計算方法分為兩大類， 一類是基于災變物理機制的水文水力學法， 另一類是基于數據驅動的統計歸納法[4]。 許多學者對臨界雨量采用水文水力學法進行了研究。 葉勇等[5]對臨界雨量進行研究， 提出了水位反推法， 此方法簡單實用且便于推廣。 原文林等[6]利用流域數字水系與美國陸軍工程兵團水文中心水文模型系統(HEC-HMS)建立模型， 獲得了精確的降雨徑流關系并計算了臨界雨量， 結果表明該模型預警精度較高。 彭萬兵等[7]采用降雨驅動指標法建立了山洪災害預警模型， 模型適用性較好。 另外， 許多學者對臨界雨量采用統計歸納法進行了研究。 張玉龍等[8]以云南省山洪災害典型區為例， 采用徑向基函數插值法、 反距離加權插值法、 克里金空間插值法研究了區域臨界雨量的變化規律， 結果表明克里金空間插值法表現較好。 樊建勇等[9]提出與臨界雨量有較大相關性的參數有主溝道比降、 流域面積和主溝道長度等小流域參數， 同時構建了統計模型推算山洪災害臨界雨量， 但該統計模型考慮參數較少。

綜上，在理論上水文水力學法較為嚴謹，但對降雨、下墊面條件、災害頻次、河道特性等資料有一定要求，且計算過程較為復雜，應用水文水力學法計算資料匱乏地區的臨界雨量存在一定困難。相對而言，統計歸納法理論較為簡單，不涉及山洪物理過程機制，多在于數理統計與規律總結，但較少對已有的臨界雨量成果數據深入挖掘。機器學習在數據的挖掘與處理方面具有較大優勢，但目前在臨界雨量預估方面應用較少。隨著全國山洪災害防治項目的推進，有關部門積累了大量采用水文水力學法計算的臨界雨量數據，在此基礎上利用機器學習算法對臨界雨量預估是可行的。

本文中以山東省臨朐縣山丘區237個受山洪災害威脅的沿河村落作為研究對象，將水文水力學法與統計歸納法緊密結合，基于隨機森林回歸算法構建臨界雨量預估模型。以190個沿河村落的臨界雨量成果數據作為訓練集訓練模型，其余47個沿河村落的臨界雨量成果數據作為測試集，對模型效果進行驗證，并分析所選參數對模型的重要度，以期為臨界雨量計算提供參考。

1 基本方法

本文中綜合水文水力學法與統計歸納法，將沿河村落作為研究對象，首先選取資料充足，計算條件好的沿河村落，采用水文水力學法計算臨界雨量，成果數據經驗證合格后作為原始數據集；其次，綜合考慮降雨特征、流域特征、沿河村落及河道特征等參數，利用隨機森林回歸算法建立山洪災害臨界雨量預估模型，計算無實測資料沿河村落的臨界雨量，計算流程如圖1所示。

圖1 山洪災害臨界雨量預估模型計算流程

1.1 水文水力學方法

水文水力學法主要包括水位流量反推法、降雨驅動指標法、降雨徑流關系曲線插值法等。本文中在計算臨界雨量時采用的是水位流量反推法。

水位流量反推法是假定降雨與徑流同頻的情況下計算控制斷面水位流量關系，確定成災水位，獲得成災流量，再經過產流、匯流計算獲得不同頻率的洪水過程，用成災流量的頻率確定降雨頻率，最后反推臨界雨量。

1.1.1 控制斷面水位流量關系計算

基于王維林等[10]對確定山丘區沿河村落成災水位方法的研究，定位沿河村落防洪能力最薄弱的居民戶的位置，確定成災水位，采用曼寧公式計算控制斷面水位流量關系，即

(1)

式中：Q為流量， m3/s；n為河道糙率，無量綱；A為過水斷面面積，m2；R為濕周，m；J為河道比降，無量綱。

1.1.2 設計暴雨計算與產匯流計算

根據當地的水文手冊與暴雨圖集，查得多年得平均雨量值，經產流計算(降雨徑流關系曲線法)、匯流計算(推理公式法)，獲得不同頻率洪水過程，并繪制設計暴雨頻率p-洪峰流量Qmp(m為匯流參數)曲線圖。

Htp=HtKP

，

(2)

(3)

(4)

m=0.061 4θ0.75

，

(5)

(6)

式中：Htp為t時間、 設計暴雨頻率為p的點雨量值，mm；Ht為t時間的平均降雨量，mm；KP由皮爾遜Ⅲ型線表查得；τ為匯流時間，h；L為河道長度，m；Qmp為洪峰流量，m3/s；F為流域面積，km2；t為時間，h；H為凈雨量，mm；θ為匯流參數。

1.1.3 臨界雨量反推計算

成災流量是成災水位對應的流量，在p-Qmp曲線圖上查出成災流量的成災頻率。假定降雨與徑流同頻，利用成災頻率設計出該頻率下典型時段的降雨量，此降雨量即為臨界雨量。

1.2 隨機森林回歸算法

1.2.1 分類回歸樹算法

隨機森林回歸(RFR)算法作為一種機器學習算法， 已被應用于多個領域[11-13]， 它是由分類回歸樹(CART)算法組合形成的。 CART是一種形成二叉樹模型的技術。 根據輸出變量類型的不同， 可將CART分為分類樹與回歸樹， 本文中使用的是回歸樹。

1.2.2 隨機森林回歸算法原理

多棵回歸決策樹構成了隨機森林回歸算法。基于集成學習的思想，取各回歸決策樹的均值作為預測結果，即

(7)

作為一種基于統計學理論的機器學習算法，隨機森林回歸算法引入了Bagging方法和隨機子空間方法[14]，避免了單棵決策樹模型容易過擬合且精度不高的問題。

1)Bagging方法[15]又被稱為自助聚集(bootstrap aggregating)，是基于自舉的統計方法。該方法以可重復的隨機采樣為基礎，通過Bootstrap重采樣方法形成多個預測器。假定原始樣本中共有N個樣本，反復抽取N次，組成新訓練樣本。當N趨向無窮大時，每個樣本不被抽中概率為

(8)

在同一棵樹的訓練樣本中將近有36.8%的原始樣本不會出現，未被抽中的樣本稱為袋外數據(OOB)[16]。通過Bagging方法可避免回歸決策樹局部最優解的產生。

2)隨機子空間方法。構建回歸決策樹時需要選取隨機特征。選取隨機特征是指從全體屬性集合中隨機選擇部分特征屬性，根據最小均方差原則選取節點分裂最優特征，讓每棵樹不剪枝以實現最大限度的生長。訓練樣本的隨機取樣根據特征屬性隨機選取，能最大程度地保證回歸決策樹的多樣性。

1.2.3 算法流程

1)利用Bagging方法，對原始訓練集隨機抽取樣本，構造出k個樣本子集。

2)利用隨機子空間方法，隨機在所有X個特征屬性中抽取部分特征屬性，進行節點分裂，構建單棵回歸決策樹。

3)重復步驟1)、2)，建立多棵回歸決策樹，并使每棵樹最大程度地生長，形成森林。

4)最終預測結果是對所有回歸決策樹的預測結果取平均值。隨機森林回歸算法計算流程如圖2所示。

2 臨界雨量預估模型構建

2.1 研究區概況及數據準備

臨朐縣地處山東省濰坊市西南部，位于彌河上游、沂山北麓，是山洪災害易發地區，尤其是隨著近年來極端天氣事件增加，暴雨天氣增多，更易形成洪水，歷史上曾在1963、 1984、 1986、 1998、 2000、 2001、 2010、 2012、 2019年發生過較大洪水。該縣30 m分辨率數字高程數據(DEM)從地理空間數據云平臺(http://www.gscloud.cn/)獲取。全縣237個研究對象分布情況，如圖3所示。

圖2 隨機森林回歸算法計算流程

從國家標準地圖網站下載，地圖審批號為GS(2019)3266號(http://bzdt.ch.mnr.gov.cn/browse.html?picId=%224o28b0625501ad13015501ad2bfc0211%22)，結合資源環境科學與數據中心2015年中國縣級行政邊界數據(http://www.resdc.cn/data.aspx?DATAID=202)， 經過ArcGIS 10.2軟件數字化處理后得到。圖3 山東省臨朐縣研究對象分布圖

根據《山東省水文圖集(1975)》查得臨朐縣土壤最大蓄水量Wmax為60 mm， 本文中將土壤含水量分為較干(0.2Wmax)、 一般(0.5Wmax)、 較濕(0.8Wmax)3種情況， 選取1、 3、 6 h作為典型時間，對全縣237個研究對象采用水位流量反推法計算3種土壤含水量下1、 3、 6 h的臨界雨量。通過2017、 2018、 2019年實測降雨資料與山洪災害調查資料對不同時間段臨界雨量計算成果檢驗復核，臨界雨量有效預警率為82.3%，數據精度較高，可作為原始數據集輸入模型。

2.2 臨界雨量影響參數選擇與獲取

臨界雨量受降雨、 土壤含水量、 下墊面情況等因素影響[17]， 并借鑒山洪災害風險評價中的指標[18-19]， 考慮數據的易獲取性選取臨界雨量影響參數。 本文中將臨界雨量影響參數分為降雨特征參數、 流域特征參數、 沿河村落及河道特征參數。

2)流域特征參數包括流域面積F、 流域平均坡度Javg、 流域最長匯流路徑Lmax、 最長匯流路徑比降Jmax、 流域形狀系數C。以上流域特征參數可以通過ArcGIS中的統計工具對DEM分析計算獲得。

3)沿河村落及河道特征參數。河道糙率n、 控制斷面寬度B、 村落河道比降Jr、 河道到居民戶水平距離S。河道糙率n為過水斷面的綜合糙率，可采用實測資料推算或參考《水工建筑物與堰槽測流規范》(SL 537—2011)表K.0.4中的內容；村落河道比降Jr可根據DEM分析計算獲得。控制斷面河道寬度B、 河岸至居民戶水平距離S，以天地圖在線影像為底圖，利用ArcGIS測量工具提取。

各臨界雨量影響參數的范圍、平均值及標準差見表1所示。

2.3 模型設置

本文中采用Python語言自帶的Scikit-learn機器學習庫，以選擇的臨界雨量影響參數為自變量，臨界雨量為因變量，利用隨機森林回歸算法，構建臨界雨量預估模型，實現臨界雨量預估。

2.3.1 訓練集、測試集

因篇幅所限，本文中以土壤含水量0.5Wmax為例，對1、 3、 6 h典型時間段分別建立模型， 將原始數據集按4∶1的比例劃分為訓練集和外部獨立測試集，使用訓練集數據訓練模型，外部獨立測試集數據僅用于模型效果評估，不用于模型建立。本文中訓練集為190組，外部獨立測試集為47組。

表1 臨界雨量影響參數統計表

2.3.2 模型參數優化

采用網格搜索方法(grid search)對模型參數尋優， 并通過K折交叉驗證確定最佳參數。 當K值較大時， 交叉驗證結果傾向于更好， 但計算時間也會更長， 綜合考慮計算時間和方差， 將K設置為5。交叉驗證中的打分器(scorer)可以通過使用評分(scoring)參數指定一個分數衡量指標對訓練結果進行評分。 分數衡量指標遵循的原則是較大的返回值比較小的返回值更好。 采用的評分指標為解釋方差回歸得分。 解釋方差回歸得分滿分為1分， 分值越小說明模型預測能力越差。 計算表達式為

(9)

隨機森林回歸算法在訓練集上的表現是通過5輪訓練和驗證的表現進行平均得到的。 模型的主要參數有決策樹數量(n_estimator)、 決策樹最大深度(max_depth)、 最大特征數(max_features)。 其中模型決策樹的數量越多， 訓練效果越好， 同時訓練所需時間也越長。 圖4所示為訓練集中決策樹數量與解釋方差回歸得分的關系曲線。 由圖可知， 決策樹數量達到300后， 解釋回歸方差得分趨向穩定， 因此在不影響計算效率的情況下， 將決策樹的數量定為300。 模型其他參數通過網格搜索方法調整， 優化后的模型參數見表2。 由表可知，1、 3、 6 h臨界雨量預估模型的得分分別為0.964、 0.960、 0.951， 表明模型訓練效果較好。

圖4 決策樹數量與解釋方差回歸得分關系曲線

3 結果與分析

3.1 模型評價指標

為了對模型進行客觀評估，采用均方根誤差RRMSE、 平均絕對誤差RMAE以及決定系數r2作為模型模擬結果與實際值的擬合程度的衡量標準，其中，當RRMSE、RMAE越小，r2越大時，模型擬合效果越好，精度越高。

(10)

(11)

(12)

表2 臨界雨量預估模型參數表

3.2 臨界雨量預估結果分析

將測試集的臨界雨量影響參數輸入到訓練好的臨界雨量預估模型中，對測試集中沿河村落的臨界雨量進行預估。以水文水力學法計算的臨界雨量作為實際值，模型計算的臨界雨量作為預估值，分別對1、 3、 6 h臨界雨量預估模型訓練集、測試集的實際值與預估值繪制散點圖，如圖5—7所示。對1、 3、 6 h臨界雨量預估模型的訓練集擬合效果分析，預估值與實際值接近，模型擬合效果較好。同時，3組模型的訓練結果顯示，當臨界雨量實際值偏大時，訓練樣本數量較少，模型的擬合效果相對較差。1、 3、 6 h臨界雨量預估模型的測試集均顯示臨界雨量預估值與實際值分布在對角線附近，說明模型預估效果較好。表3列出了測試集中部分沿河村落的臨界雨量預估結果。

為進一步比較隨機森林回歸算法在臨界雨量預估問題上的優劣，分別采用隨機森林回歸算法與誤差逆傳播(BP)神經網絡算法預估臨界雨量。BP神經網絡算法使用3層結構，經交叉驗證優化后，選擇tansig函數作為輸入層到隱含層傳遞函數，purelin函數作為隱含層到輸出層傳遞函數，采用levenberg-marquardt算法進行訓練，迭代次數為3 000。

(a)訓練集(b)測試集圖5 1 h臨界雨量預估模型預估值與實際值的對比

(a)訓練集(b)測試集圖6 3 h臨界雨量預估模型預估值與實際值的對比

(a)訓練集(b)測試集圖7 6 h臨界雨量預估模型預估值與實際值的對比

表3 不同臨界雨量預估模型的測試集預估結果

表4給出了隨機森林回歸算法與BP神經網絡算法在預估臨界雨量時訓練集﹑測試集的模型評價指標。由表可以看出，采用隨機森林回歸算法的1、 3、 6 h臨界雨量預估模型中， 訓練集的r2均大于0.952， 測試集的r2均大于0.946，說明所建立的隨機森林回歸算法能較好地擬合建模數據；對于訓練集與測試集的誤差，以實際值的10%為許可誤差，訓練集與測試集的預估值合格率均大于80%，且RRMSE與RMAE均小于10，誤差可接受。此外，測試集與訓練集相比，模型評價指標相差不大，說明模型未出現過擬合，泛化性能較好。

通過比較隨機森林回歸算法與BP神經網絡算法的模型評價指標可知， 各時段采用隨機森林回歸算法的訓練集與測試集的r2比BP神經網絡算法的至少增大了6%， 擬合效果更好； 各時段隨機森林回歸算法的RRMAE與RMAE均小于BP神經網絡算法的， 且預估值合格率更高。 綜上， 基于隨機森林回歸算法的臨界雨量預估模型預估結果與BP神經網絡算法相比， 預估精度更高， 能更好地進行臨界雨量預估。

3.3 重要度分析

在臨界雨量預估模型中，對每個臨界雨量影響參數在每棵樹上的貢獻求平均值，經比較后可度量每個臨界雨量影響參數的重要度。在Python語言中使用Scikit-learn機器學習庫的feature_importances函數可以直接得到。圖8所示為各臨界雨量影響參數的重要度。由圖可知，臨界雨量影響參數中最重要的是流域面積，重要度值為0.320，流域最長匯流路徑比降、 24 h降雨均值、流域平均坡度的重要度也相對較高，相比之下，河道糙率、村落河道比降的重要度相對較低。

表4 采用不同算法的模型評價指標

圖8 臨界雨量影響參數重要度

4 結論

本文中以山東省臨朐縣山丘區237個沿河村落為例，基于隨機森林回歸算法構建山洪災害臨界雨量預估模型，并得到以下結論：

1)利用水文水力學法計算沿河村落臨界雨量數據與相應的臨界雨量影響參數，采用隨機森林回歸算法構建臨界雨量預估模型，將無實測資料的沿河村落臨界雨量影響參數輸入到模型中，實現了對各村落臨界雨量的預估，為預估臨界雨量提供一種有效的方法。

2)通過對模型進行訓練與測試表明，該模型具有很好的預估能力，在訓練集和測試集上的r2均大于0.9，說明該模型擬合情況較好；測試集中臨界雨量預估值與水位流量反推法計算的臨界雨量實際值接近，且RMAE和RRMSE均小于10，模型預估結果精度較高，能滿足實際工作的需要。同時隨機森林回歸算法結構簡單，需要調整的參數較少，能夠針對無實測資料的沿河村落進行快速、批量預估臨界雨量。

3)模型影響參數中流域面積、流域最長匯流路徑比降、降雨均值、流域平均坡度的重要度相對較高，因此應確保這些數據的精度，提高模型的準確性。

4)本文中僅對山東省臨朐縣沿河村落的臨界雨量進行預估，樣本數據較少，隨著樣本數據的增加，模型的性能還需進一步研究。