基于集成學習的巢湖面雨量計算方法

王 杰，徐 祥，羅曉丹，張 萌，黃 澈，洪冠中，汪 翔

（1.安徽省農(nóng)村綜合經(jīng)濟信息中心，安徽 合肥 230031；2.廣東省清遠市氣象局，廣東 清遠 511515；3.安徽省公共氣象服務(wù)中心，安徽 合肥 230031）

0 引 言

巢湖是中國五大淡水湖之一，地處中緯度地帶，由于其位于南北氣候過渡帶，冷暖氣團常在此交匯而造成暴雨等惡劣天氣頻發(fā)［1］。建國以來，巢湖發(fā)生較大洪澇災(zāi)害的年份共計18 年［2］，尤其2020 年最為嚴重，與巢湖相關(guān)的主要河湖洪水位均突破歷史極值，嚴重威脅人民群眾的生命財產(chǎn)安全。巢湖各條支流、水庫、堤防的抗洪能力都與流域和湖面面雨量息息相關(guān)，因此開展面雨量計算和監(jiān)測可為應(yīng)急管理部門提供數(shù)據(jù)信息支撐，為防汛抗洪的科學決策和合理部署提供重要依據(jù)。

目前，面雨量的計算方法主要有等值線法、數(shù)值法、算術(shù)平均法和網(wǎng)格點平均法等［3］。其中，等值線法計算較為復雜，難以滿足實際應(yīng)用要求。數(shù)值法是指泰森多邊形法或三角形法，該方法通過確定區(qū)域內(nèi)各氣象觀測站的權(quán)重系數(shù)計算面雨量，適用于氣象觀測站或降雨量分布不均勻的區(qū)域。算術(shù)平均法是對區(qū)域內(nèi)各站測得的同期雨量求算術(shù)平均從而計算面雨量，僅適用于區(qū)域面積小、地形起伏不大且觀測站多而分布又較均勻的場景［4］。網(wǎng)格點平均法［5］是利用GIS 技術(shù)將區(qū)域按設(shè)定的經(jīng)緯度分辨率進行網(wǎng)格劃分，通過計算網(wǎng)格點上的降雨量實現(xiàn)面雨量計算，其中網(wǎng)格點降雨量多基于天氣雷達、氣象衛(wèi)星和地面觀測等資料的某種或多種的融合進行估算，從時空分辨率精度及空間覆蓋范圍方面考慮，采用雷達數(shù)據(jù)估測網(wǎng)格點降水具有明顯優(yōu)勢。

傳統(tǒng)的雷達估測降水常采用Z-R 關(guān)系法構(gòu)建雷達降水估測模型［6］，通過雷達反射率因子Z 直接反演降水R。Z-R 關(guān)系隨地區(qū)、季節(jié)、降水類型而變化，甚至同一次降水過程中亦可能發(fā)生變化，這種變化直接影響雷達降水估測精度［7］。前人研究多基于計算復雜、專業(yè)性較強的數(shù)值模式擬合Z-R 關(guān)系，而以集成學習為代表的人工智能技術(shù)在雷達估測降水中的應(yīng)用鮮有報道。近年來，隨著人工智能技術(shù)的發(fā)展，集成學習因其優(yōu)越的泛化性能，被學者們應(yīng)用于氣象方面研究。張鈞民等［8］利用XGBoost 算法開展多源降水數(shù)據(jù)融合方法研究；于霞等［9］以雷雨大風的46 個相關(guān)特征屬性數(shù)據(jù)為輸入，降雨和大風的數(shù)值為輸出，利用XGBoost 算法建立雷雨大風預(yù)測模型；陳志月等［10］選用GPR、CatBoost、XGBoost 這3 種集成學習算法預(yù)測江西地區(qū)水面蒸發(fā)量；陳雙等［11］基于云頂溫度、中層融化參數(shù)等構(gòu)建的決策樹判別模型，可較好地提升臨界氣溫下雨、雪的判別準確率。集成學習雖在氣象領(lǐng)域不同場景開展應(yīng)用研究，但此類算法在雷達降水估測的適用性還需進一步探究。

本文以巢湖為研究對象，提出一種基于集成學習算法和網(wǎng)格點平均法相結(jié)合的面雨量計算方法，通過對巢湖流域邊界內(nèi)的雷達三維拼圖數(shù)據(jù)和氣象觀測站降水數(shù)據(jù)進行預(yù)處理，構(gòu)建以雷達多層拼圖數(shù)據(jù)為特征、氣象觀測站小時降水量為標簽的數(shù)據(jù)集，利用所構(gòu)建數(shù)據(jù)集對模型進行訓練迭代，結(jié)合GIS 技術(shù)，將巢湖區(qū)域按經(jīng)緯度網(wǎng)格劃分并與雷達拼圖在空間上進行疊加，根據(jù)QPE 模型計算各網(wǎng)格點降雨量，從而實現(xiàn)巢湖面雨量的自動估算，以期為提升巢湖防災(zāi)減災(zāi)能力提供重要信息支撐。

1 資料與方法

1.1 研究區(qū)域概況

巢湖流域地處安徽省中部，位于東經(jīng)116°24'～118°0'，北緯30°58'～32°06'，面積為13486 km2［12-13］。巢湖作為巢湖流域的中心，從南、西、北三面匯入，西面、北面均為山區(qū)，是巢湖重要支流的發(fā)源地，東面經(jīng)巢湖閘由裕溪河及牛屯河注入長江（見圖1）。巢湖及流域地勢輪廓都是東西長南北窄，流域內(nèi)地形復雜多樣，有山地、平原、丘陵等多種地形類型，且西高東低，中間地勢低洼平坦，區(qū)域地形地貌極易形成山洪沖擊、江洪頂托等不利因素，從而使得巢湖防汛抗洪工作面臨嚴峻挑戰(zhàn)［13］。

圖1 巢湖流域及氣象觀測站點分布圖

1.2 數(shù)據(jù)來源

本 文以2019 年6 月5 日01 時—6 日01 時 的一次降水過程為例，開展巢湖面雨量計算方法研究。數(shù)據(jù)包括雷達三維拼圖數(shù)據(jù)、氣象觀測站降水數(shù)據(jù)、巢湖流域DEM 高程數(shù)據(jù)和基于氣象陸面數(shù)據(jù)同化系統(tǒng)（CMA Land Data Assimilation System，CLDAS）的降水數(shù)據(jù)［14］。其中，雷達三維拼圖數(shù)據(jù)來自災(zāi)害天氣短時臨近預(yù)報系統(tǒng)（Severe Weather Automatic Nowcast System，SWAN）［15］，主要包含經(jīng)緯度、時間、層數(shù)、站點信息等29 個字段信息，拼圖按高度0.5 km～19 km 排列共21 層，每層包含800×900 個網(wǎng)格，空間分辨率0.01°×0.01°，時間分辨率為6 min。圖2為3 km 高度巢湖流域的雷達拼圖，6 月5 日14 時雷達回波強度普遍小于21時，這表明6月5日21時降水概率和強度普遍大于14時。降水數(shù)據(jù)和CLDAS數(shù)據(jù)均來自氣象大數(shù)據(jù)云平臺“天擎”，降水數(shù)據(jù)包括巢湖流域范圍內(nèi)375 個氣象觀測站，其中國家站10 個，區(qū)域站365 個（見圖1）；CLDAS 數(shù)據(jù)為空間分辨率為0.0625°×0.0625°，時間分辨率為1 h 的降水資料。本文雷達三維拼圖數(shù)據(jù)和氣象觀測站降水數(shù)據(jù)作為構(gòu)建雷達降水估測模型的數(shù)據(jù)集，基于CLDAS 得到的面雨量數(shù)據(jù)用于檢驗本文模型的精度。

圖2 3 km雷達拼圖回波變化

1.3 面雨量計算方法

1.3.1 基于集成學習算法構(gòu)建雷達QPE模型

鑒于集成學習算法能較好地反應(yīng)地面降水與雷達拼圖各層降水數(shù)據(jù)之間的非線性關(guān)系，為提高雷達估測降水的精度，本文選用3 種經(jīng)典集成學習算法RF、XGBoost、LightGBM 開展雷達估測降水方法研究，以期獲得性能較好的雷達QPE模型。

RF 屬于Bagging 類算法，通過自助法（Bootstrap）對訓練集樣本和特征有放回地隨機采樣，構(gòu)建許多相互不關(guān)聯(lián)的回歸決策樹［16］。RF 并行訓練所有決策樹，每棵決策樹得到一個預(yù)測結(jié)果，將所有樹的預(yù)測結(jié)果結(jié)合取平均值，從而得到整個模型的預(yù)測結(jié)果。XGBoost、LightGBM 屬于Boosting 類算法，兩者是梯度提 升 決 策 樹（Gradient Boosting Decision Tree，GBDT）的進化版本［17］，其中，XGBoost 是在決策樹構(gòu)建過程中加入了正則項，還使用一階、二階梯度函數(shù)來優(yōu)化損失函數(shù)，這樣做是為了避免過擬合和降低模型復雜度［17-18］；而LightGBM 是通過采用直方圖算法（Histogram）、基于梯度的單邊采樣算法（GOSS）、互斥特征捆綁算法（EFB）分別對分裂點數(shù)、樣本量和特征數(shù)量3個變量進行優(yōu)化［19］，并在葉子生長策略上采用帶深度限制的按葉子生長（leaf-wise）算法和對緩存命中率進行優(yōu)化，在提高預(yù)測效率和準確率的同時防止模型過擬合。引入3 種集成學習算法建立的雷達QPE模型如式（1）：

其中，j表示RF、XGBoost、LightGBM表示采用集成學習算法獲得的降水預(yù)測數(shù)據(jù)；fEM表示基于集成學習算法構(gòu)建的非線性回歸預(yù)測模型；Radarlayeri表示雷達拼圖的第i層降水數(shù)據(jù)。本文選擇雷達拼圖的5 個低仰角數(shù)據(jù)作為訓練集的特征輸入，這樣做一方面是由于雷達降水信息主要反映在雷達體掃5 個低仰角上，另一方面可以減少模型訓練數(shù)據(jù)量。

1.3.2 基于網(wǎng)格點降雨量計算面雨量

巢湖面雨量是整個湖面區(qū)域單位面積上的平均降水量［4］?；诰W(wǎng)格點降雨量的計算方法可以利用區(qū)域內(nèi)一切可獲得的降水信息，能較精確地計算面雨量，故本文采用基于GIS 技術(shù)的網(wǎng)格點平均法計算巢湖面雨量。其具體實現(xiàn)步驟如下：

1）通過GIS 技術(shù)，將巢湖區(qū)域按0.01°×0.01°經(jīng)緯度網(wǎng)格劃分并與雷達拼圖在空間上進行疊加，建立雷達拼圖經(jīng)緯度網(wǎng)格與巢湖區(qū)域的映射關(guān)系，并對湖面邊界內(nèi)網(wǎng)格點進行標記（見圖3）。

圖3 雷達拼圖與巢湖邊界的映射關(guān)系

2）根據(jù)雷達QPE模型計算各網(wǎng)格點的降雨量。

3）基于步驟2 的計算結(jié)果，對湖面區(qū)域所有網(wǎng)格點降雨量進行算術(shù)平均，以實現(xiàn)湖面面雨量估算。具體計算公式如式（2）所示，其中n為網(wǎng)格數(shù)為雷達QPE模型預(yù)測結(jié)果。

1.4 模型超參數(shù)調(diào)優(yōu)

不同超參數(shù)會對模型性能造成較大差異［20］，鑒于超參數(shù)調(diào)優(yōu)需要花費大量的計算時間和資源，本文選擇網(wǎng)格搜索法提高超參數(shù)搜索效率，通過在Python語言環(huán)境下加載sklearn（Version: 0.23.2）開源庫的GridSearchCV函數(shù)實現(xiàn)［21-22］。具體步驟如下：

1）根據(jù)RF、XGBoost、LightGBM 這3 種算法模型特點和超參數(shù)的原理選擇需要調(diào)整的超參數(shù)組合，并將每一個超參數(shù)值設(shè)置一個比較寬松的調(diào)參區(qū)間［23］。

2）采用網(wǎng)格搜索法進行訓練迭代，每次訓練結(jié)束后根據(jù)模型評價指標和基于CLDAS 得到的面雨量數(shù)據(jù)評估調(diào)參效果，逐漸縮小調(diào)參范圍，直至得到最優(yōu)調(diào)優(yōu)結(jié)果。

3）根據(jù)步驟2 中調(diào)優(yōu)結(jié)果，在超參數(shù)組合中選擇某一超參數(shù)為變量，其他超參數(shù)值按步驟2 中結(jié)果保持不變。改變此超參數(shù)的值開展訓練迭代驗證對模型性能的影響，若影響程度較小則可將此超參數(shù)剔出最佳參數(shù)組合，后續(xù)訓練中使用超參數(shù)默認值。反之，此超參數(shù)則是影響模型性能關(guān)鍵超參數(shù)。其他超參數(shù)按此方法依次進行驗證。

1.5 模型評價指標

本文以決定系數(shù)（R2）、平均絕對誤差（MAE）和均方根誤差（RMSE）作為雷達QPE 模型的評價指標［24］，計算公式為：

2 實驗與結(jié)果分析

2.1 實驗環(huán)境

實驗硬件環(huán)境為16 GB NVIDIA Tesla V100 GPU，操作系統(tǒng)為CentOS Linux7.5，選用的機器學習開源架構(gòu)為Tensorflow（Version: 2.3.0），通過Tensorflow調(diào)用GPU訓練［25］雷達QPE回歸預(yù)測模型。

2.2 數(shù)據(jù)集制作

表1 給出了氣象觀測站降水數(shù)據(jù)預(yù)處理部分數(shù)據(jù)實例。

表1 氣象觀測站降水數(shù)據(jù)預(yù)處理部分數(shù)據(jù)實例

1）氣象觀測站降水數(shù)據(jù)預(yù)處理：流域內(nèi)氣象觀測站降水原始數(shù)據(jù)有9377 條，為了使訓練的結(jié)果更加真實有效，剔除降水缺測記錄、異常記錄、沒有通過質(zhì)量控制的數(shù)據(jù)和不完整觀測記錄。同時，由于0.1 mm 的降水數(shù)據(jù)為非有效數(shù)據(jù)，會影響訓練模型的評價指標［26］，因此刪除降水量為0.1 mm 的觀測記錄。經(jīng)過上述操作后最終篩選出7376條有效數(shù)據(jù)。

2）雷達三維拼圖數(shù)據(jù)預(yù)處理：提取流域內(nèi)氣象觀測站空間經(jīng)緯度對應(yīng)的雷達21 個仰角反射率，將6 min/次的拼圖在特定時次內(nèi)求算術(shù)平均得到逐小時平均反射率。表2 中列舉了處理后的逐小時平均反射率（layer0～layer4）部分數(shù)據(jù)實例。其中，dBZ 表示雷達反射率因子的單位。

表2 模型數(shù)據(jù)集部分實例

3）制作以小時降水量為標簽的數(shù)據(jù)集：提取氣象觀測站各時次降水量，按站點和時間，標注雷達數(shù)據(jù)，形成以雷達5 個低仰角逐小時平均反射率（layer 0～layer4）為特征，以小時降水量（pre）為標簽的數(shù)據(jù)集，建立氣象觀測站小時降水量與雷達逐小時平均反射率映射。

2.3 模型訓練

數(shù)據(jù)集分為訓練集和驗證集，訓練集和驗證集按8:2 劃分。模型訓練流程如圖4 所示，采用RF、XGBoost、LightGBM 這3 種集成學習算法模型分別開展訓練迭代，并基于模型評價指標進行不同集成算法下模型性能對比。設(shè)定不同超參數(shù)開展訓練對比實驗，訓練完成記錄每次訓練時長、GPU 內(nèi)存消耗、評價指標值、輸出2019 年6 月5 日01 時—6 月6 日00 時巢湖實況面雨量并保存雷達QPE模型。

圖4 基于集成學習的巢湖面雨量訓練流程圖

2.4 超參數(shù)調(diào)優(yōu)結(jié)果

將每次訓練結(jié)束后超參數(shù)的調(diào)優(yōu)結(jié)果和效果進行對比，得到RF、XGBoost、LightGBM 模型最佳超參數(shù)組合和調(diào)優(yōu)結(jié)果。表3以數(shù)據(jù)形式直觀給出了3種算法模型最佳超參數(shù)組合、獲得調(diào)優(yōu)結(jié)果前確定的調(diào)參范圍和得到調(diào)優(yōu)結(jié)果的訓練時間。其中，n_estimators 表示回歸樹數(shù)量；max_depth 表示樹的深度；min_samples_split 表示拆分內(nèi)部節(jié)點所需最少樣本數(shù)；min_samples_leaf 表示在葉節(jié)點處需要最小樣本數(shù)；min_child_weight 表示子節(jié)點所需的實例權(quán)重的最小總和；learning_rate 表示學習率；num_leaves 表示每棵回歸樹上的葉子數(shù)。調(diào)參前超參數(shù)數(shù)值代表含義，以n_estimators 為例，（100，500，10）表示調(diào)參范圍為100～500，步長為10。

表3 模型超參數(shù)組合和訓練時間

本文以RF 模型的n_estimators 為例，驗證表3 調(diào)優(yōu)結(jié)果的可信度和超參數(shù)是否為模型關(guān)鍵超參數(shù)。圖5 為n_estimators 值改變時，MAE、RMSE 和R2這3個評價指標的變化情況。由圖4可見，當n_estimators變化時，對模型評價指標影響明顯，結(jié)合n_estimators=480 時RF 模型面雨量預(yù)測結(jié)果（見表4），可以判定n_estimators 是RF 模型的關(guān)鍵超參數(shù)，表3 中RF 模型的n_estimators調(diào)優(yōu)結(jié)果可信。

表4 集成學習算法+網(wǎng)格點平均法與CLDAS面雨量計算結(jié)果比較 單位：mm

圖5 不同超參數(shù)值對模型性能的影響

2.5 集成學習算法模型評估

經(jīng)過多次訓練迭代，采用RF、XGBoost 和Light-GBM 這3 種集成學習算法模型預(yù)測降水的評價指標結(jié)果如表5所示。

表5 不同算法模型評價指標結(jié)果

由表5 可知，3 種算法模型訓練時長和內(nèi)存開銷差別很小。在訓練期最優(yōu)模型是RF 預(yù)測模型（MAE=0.492 mm/h，RMSE=1.254 mm/h，R2=0.933）；在驗證期XGBoost 的R2值（R2=0.697）雖然略高，但MAE 和RMSE 表現(xiàn)不如其他2 種模型，綜合分析各項評價指標，RF 預(yù)測模型（MAE=1.342 mm/h，RMSE=3.435 mm/h，R2=0.689）優(yōu) 于XGBoost 和LightGBM 模 型。綜上，基于RF 算法的雷達QPE 估測模型在本次降水過程中表現(xiàn)最佳。

2.6 實驗結(jié)果驗證

輸出RF、XGBoost、LightGBM 這3 種模型2019 年6 月5 日01—6 月6 日00 時巢湖面雨量的計算結(jié)果，并與基于CLDAS 得到的面雨量結(jié)果對比，以此驗證基于集成學習算法和網(wǎng)格點平均法相結(jié)合的方法在實際場景中應(yīng)用效果。

由表4 可知，在2019 年6 月5 日01 時—12 時，基于RF 的雷達QPE 模型和網(wǎng)格點平均法相結(jié)合的面雨量計算結(jié)果與CLDAS 中面雨量計算結(jié)果均為0 mm；在13 時—21 時和23 時，兩者結(jié)果相差在1 mm以內(nèi)，但在22 時和2019 年6 月6 日00 時兩者計算結(jié)果相差均大于1 mm，這點在圖6中體現(xiàn)得更為直觀。圖6是4種模型面雨量計算結(jié)果趨勢圖，由圖可見，基于RF 的雷達QPE 模型和網(wǎng)格點平均法相結(jié)合的面雨量計算方法與基于CLDAS 的計算結(jié)果更加接近，但在后面3 個時刻雖然4 種計算方法結(jié)果趨勢一致，但離散程度增大，XGBoost偏離程度最大。

圖6 4種巢湖面雨量計算方法結(jié)果趨勢圖

3 結(jié)束語

本文提出了一種基于集成學習算法與網(wǎng)格點平均法相結(jié)合的巢湖面雨量計算方法，通過對巢湖流域內(nèi)雷達三維拼圖數(shù)據(jù)和氣象觀測站降水數(shù)據(jù)進行預(yù)處理，建立了模型訓練數(shù)據(jù)集，構(gòu)建了3 種基于集成學習算法的雷達QPE 模型，使用GridSearchCV 函數(shù)對模型超參數(shù)進行了調(diào)優(yōu)，評價了3 種集成學習算法在雷達降水預(yù)測中的表現(xiàn)，最后引入基于CLDAS 的面雨量計算方法對本文面雨量計算方法進行了驗證。結(jié)果表明：

1）RF、XGBoost 和LightGBM 這3 種網(wǎng)絡(luò)模型中，RF模型在驗證期表現(xiàn)最佳，評價指標MAE、RMSE、R2分別為1.342 mm/h、3.435 mm/h、0.689，為適用于本次降水過程的預(yù)測模型。

2）對比分析了本文與CLDAS 這2 種方法得到的湖面面雨量計算結(jié)果，計算結(jié)果數(shù)值雖有差異，但總體趨勢一致，說明使用本文提出的方法計算湖面面雨量可行。

本文為進一步開展巢湖及其流域防汛抗洪方法研究、建立科學合理的湖面面雨量閾值、增強防汛抗洪監(jiān)測預(yù)警能力［27］提供了新的方法，當前可作為數(shù)值模式面雨量計算的有力補充。但本文方法還存在一些不足，對比本文方法得到的面雨量結(jié)果與數(shù)值模式計算結(jié)果在后面3 個時刻發(fā)現(xiàn)，計算結(jié)果與實際值偏離程度增大，說明有些因素沒有考慮到。下一步將研究地形、降水相態(tài)、氣象觀測站與雷達拼圖在空間不同對應(yīng)關(guān)系等因素對模型的影響，以期進一步提高巢湖面雨量計算精度。