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

(1. 濟南大學(xué)水利與環(huán)境學(xué)院， 山東濟南 250022； 2. 中國水利水電科學(xué)研究院， 北京 100038； 3. 水發(fā)規(guī)劃設(shè)計有限公司， 山東濟南 250100)

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

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

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

1 基本方法

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

圖1 山洪災(zāi)害臨界雨量預(yù)估模型計算流程

1.1 水文水力學(xué)方法

水文水力學(xué)法主要包括水位流量反推法、降雨驅(qū)動指標(biāo)法、降雨徑流關(guān)系曲線插值法等。本文中在計算臨界雨量時采用的是水位流量反推法。

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

1.1.1 控制斷面水位流量關(guān)系計算

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

(1)

式中：Q為流量， m3/s；n為河道糙率，無量綱；A為過水?dāng)嗝婷娣e，m2；R為濕周，m；J為河道比降，無量綱。

1.1.2 設(shè)計暴雨計算與產(chǎn)匯流計算

根據(jù)當(dāng)?shù)氐乃氖謨耘c暴雨圖集，查得多年得平均雨量值，經(jīng)產(chǎn)流計算(降雨徑流關(guān)系曲線法)、匯流計算(推理公式法)，獲得不同頻率洪水過程，并繪制設(shè)計暴雨頻率p-洪峰流量Qmp(m為匯流參數(shù))曲線圖。

Htp=HtKP

，

(2)

(3)

(4)

m=0.061 4θ0.75

，

(5)

(6)

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

1.1.3 臨界雨量反推計算

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

1.2 隨機森林回歸算法

1.2.1 分類回歸樹算法

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

1.2.2 隨機森林回歸算法原理

多棵回歸決策樹構(gòu)成了隨機森林回歸算法?；诩蓪W(xué)習(xí)的思想，取各回歸決策樹的均值作為預(yù)測結(jié)果，即

(7)

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

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

(8)

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

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

1.2.3 算法流程

1)利用Bagging方法，對原始訓(xùn)練集隨機抽取樣本，構(gòu)造出k個樣本子集。

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

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

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

2 臨界雨量預(yù)估模型構(gòu)建

2.1 研究區(qū)概況及數(shù)據(jù)準(zhǔn)備

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

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

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

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

2.2 臨界雨量影響參數(shù)選擇與獲取

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

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

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

各臨界雨量影響參數(shù)的范圍、平均值及標(biāo)準(zhǔn)差見表1所示。

2.3 模型設(shè)置

本文中采用Python語言自帶的Scikit-learn機器學(xué)習(xí)庫，以選擇的臨界雨量影響參數(shù)為自變量，臨界雨量為因變量，利用隨機森林回歸算法，構(gòu)建臨界雨量預(yù)估模型，實現(xiàn)臨界雨量預(yù)估。

2.3.1 訓(xùn)練集、測試集

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

表1 臨界雨量影響參數(shù)統(tǒng)計表

2.3.2 模型參數(shù)優(yōu)化

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

(9)

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

圖4 決策樹數(shù)量與解釋方差回歸得分關(guān)系曲線

3 結(jié)果與分析

3.1 模型評價指標(biāo)

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

(10)

(11)

(12)

表2 臨界雨量預(yù)估模型參數(shù)表

3.2 臨界雨量預(yù)估結(jié)果分析

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

為進一步比較隨機森林回歸算法在臨界雨量預(yù)估問題上的優(yōu)劣，分別采用隨機森林回歸算法與誤差逆?zhèn)鞑?BP)神經(jīng)網(wǎng)絡(luò)算法預(yù)估臨界雨量。BP神經(jīng)網(wǎng)絡(luò)算法使用3層結(jié)構(gòu)，經(jīng)交叉驗證優(yōu)化后，選擇tansig函數(shù)作為輸入層到隱含層傳遞函數(shù)，purelin函數(shù)作為隱含層到輸出層傳遞函數(shù)，采用levenberg-marquardt算法進行訓(xùn)練，迭代次數(shù)為3 000。

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

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

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

表3 不同臨界雨量預(yù)估模型的測試集預(yù)估結(jié)果

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

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

3.3 重要度分析

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

表4 采用不同算法的模型評價指標(biāo)

圖8 臨界雨量影響參數(shù)重要度

4 結(jié)論

本文中以山東省臨朐縣山丘區(qū)237個沿河村落為例，基于隨機森林回歸算法構(gòu)建山洪災(zāi)害臨界雨量預(yù)估模型，并得到以下結(jié)論：

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

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

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

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