梯度提升決策樹(shù)在雷達(dá)定量降水估測(cè)中的應(yīng)用研究

陳訓(xùn)來(lái) 陳元昭 趙春陽(yáng) 張珂

（1 深圳市氣象局，深圳 518040；2 深圳南方強(qiáng)天氣研究重點(diǎn)實(shí)驗(yàn)室，深圳 518040；3 阿里云計(jì)算有限公司，杭州 310099）

0 引言

定量降水估測(cè)（QPE）是現(xiàn)代天氣業(yè)務(wù)體系的重要組成部分，一直以來(lái)是業(yè)務(wù)預(yù)報(bào)的難點(diǎn)和重點(diǎn)，也是定量降水預(yù)報(bào)、強(qiáng)降水短臨預(yù)警的基礎(chǔ)，對(duì)城市內(nèi)澇、地質(zhì)災(zāi)害精細(xì)化預(yù)報(bào)起重要支撐作用。隨著新一代多普勒天氣雷達(dá)觀測(cè)網(wǎng)的逐步建成和業(yè)務(wù)運(yùn)行，應(yīng)用多普勒天氣雷達(dá)聯(lián)合自動(dòng)氣象站雨量開(kāi)展定量降水估測(cè)已成為一種較為成熟可靠的實(shí)現(xiàn)手段[1]。

利用多普勒天氣雷達(dá)定量降水估測(cè)，相對(duì)于用自動(dòng)氣象站雨量觀測(cè)，具有測(cè)量范圍廣、時(shí)空分辨率高、可及時(shí)取得大面積定量降水資料的優(yōu)點(diǎn)。由于雷達(dá)回波（反射率因子）反映的是空中水滴、冰晶等對(duì)雷達(dá)電磁波的反射或散射，因此雷達(dá)回波強(qiáng)度（Z）與降水強(qiáng)度（I）之間存在密切的指數(shù)型關(guān)系，即Z=A×Ib[2]。通過(guò)Z與I之間的對(duì)應(yīng)關(guān)系就可以估測(cè)出雷達(dá)覆蓋范圍內(nèi)的降水強(qiáng)度。降水定量估測(cè)精準(zhǔn)度在很大程度上取決于Z-I關(guān)系式中系數(shù)A和b，降雨常具有時(shí)空特性強(qiáng)、變化快的特點(diǎn)，且不同季節(jié)、不同地點(diǎn)的降水性質(zhì)不同，其Z-I關(guān)系式也不同，常用來(lái)建立Z-I關(guān)系的方法有最優(yōu)化法[3]、概率密度法[4-5]、動(dòng)態(tài)關(guān)系法[6]、遺傳算法[7]和動(dòng)態(tài)分級(jí)法[8]等。為了得到更高精度的降水估計(jì)值，通常利用自動(dòng)氣象站雨量資料對(duì)雷達(dá)降水估測(cè)值做進(jìn)一步訂正，包括應(yīng)用平均校準(zhǔn)法[9]、變分同化法[10]、客觀插值分析法[11]和卡爾曼濾波校準(zhǔn)法[12]等。

隨著人工神經(jīng)網(wǎng)絡(luò)（ANN）等深度學(xué)習(xí)技術(shù)的發(fā)展，很多專家學(xué)者利用深度學(xué)習(xí)的算法開(kāi)展雷達(dá)定量降水估測(cè)的研究，由于深度學(xué)習(xí)是一種非參數(shù)運(yùn)算，且具有自適應(yīng)學(xué)習(xí)方法，能有效減少誤差，提高了降水估測(cè)的準(zhǔn)確性[13-16]。Liu等[13]研究發(fā)現(xiàn)基于人工神經(jīng)網(wǎng)絡(luò)的定量降水估測(cè)比用Z-I關(guān)系更精確。Chiang等[14]利用臺(tái)灣六次臺(tái)風(fēng)降雨過(guò)程研究表明，采用動(dòng)態(tài)人工神經(jīng)網(wǎng)絡(luò)比Z-I關(guān)系有更高的精度和穩(wěn)定性。邵月紅等[15]通過(guò)BP神經(jīng)網(wǎng)絡(luò)方法估測(cè)定量降雨，精度要優(yōu)于Z-I關(guān)系式，估測(cè)的降雨量與站點(diǎn)實(shí)測(cè)雨量吻合性較好。傅德勝等[16]基于徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)（RBF）建立的雷達(dá)定量估測(cè)降水模型，其對(duì)降水估測(cè)精度和穩(wěn)定性高于以變分法得到的Z-I關(guān)系。

目前針對(duì)華南地區(qū)開(kāi)展基于深度學(xué)習(xí)的定量降水估測(cè)算法的研究還比較少，本文主要應(yīng)用三維網(wǎng)格雷達(dá)回波資料和自動(dòng)氣象站數(shù)據(jù)開(kāi)展基于梯度提升決策樹(shù)（GBDT）算法在雷達(dá)定量估測(cè)降水中的研究和應(yīng)用，提高雷達(dá)定量降水估測(cè)的精度。

1 數(shù)據(jù)和預(yù)處理

1.1 數(shù)據(jù)

本研究使用的雷達(dá)資料為廣東省境內(nèi)廣州、深圳、韶關(guān)、清遠(yuǎn)、陽(yáng)江、河源、汕尾、梅州、湛江和肇慶等10部S波段多普勒雷達(dá)在VCP21探測(cè)模式下，每隔6 min對(duì)9個(gè)仰角進(jìn)行三維體掃描，為確保所有雷達(dá)能同步掃描，采用了北京敏視達(dá)公司的雷達(dá)同步控制和數(shù)據(jù)采集系統(tǒng)時(shí)間同步軟件進(jìn)行時(shí)間同步。為計(jì)算方便，在雷達(dá)三維體掃描原始數(shù)據(jù)基礎(chǔ)上，利用Cressman插值法將極坐標(biāo)格式的雷達(dá)探測(cè)資料插值到三維直角坐標(biāo)系中，并采用了蘭紅平等[17]提出的“膨脹-侵蝕”算法對(duì)回波缺陷進(jìn)行質(zhì)量控制。由于在雷達(dá)組網(wǎng)拼圖時(shí)，若低于1.5 km，任意兩部雷達(dá)之間幾乎沒(méi)有共同的覆蓋區(qū)域，觀測(cè)的范圍太小；若高于5 km，雷達(dá)覆蓋的范圍擴(kuò)大了，然而卻難以形成降水[18]。因此，本研究中用到經(jīng)過(guò)雷達(dá)回波質(zhì)量控制和拼圖處理后得到1.5、2.5、3.5和4.5 km共四層高度上CAPPI（constant altitude plan position indicator）數(shù)據(jù)，空間分辨率為0.01°×0.01°，水平網(wǎng)格點(diǎn)為約為500 km×300 km，數(shù)據(jù)覆蓋為2015—2017年的每年4—9月的降水過(guò)程。

另外，本研究中還用到該范圍內(nèi)區(qū)域自動(dòng)氣象站逐小時(shí)雨量觀測(cè)資料，自動(dòng)氣象站分布如圖1。

圖1 自動(dòng)氣象站分布圖（黑點(diǎn)） Fig. 1 Distribution of automatic weather stations(AWSs) (black dots)

1.2 預(yù)處理

由于雷達(dá)回波為空中水滴、冰晶等對(duì)雷達(dá)電磁波的反射或散射，自動(dòng)氣象站記錄的是落到地面的水滴，在選取雷達(dá)的回波強(qiáng)度數(shù)據(jù)時(shí)，需要綜合考慮多種因素的影響。梁建茵等[19]認(rèn)為，對(duì)低緯度地區(qū)，自動(dòng)氣象站點(diǎn)上降水強(qiáng)度與其上空約5 km×5 km區(qū)域內(nèi)的平均雷達(dá)回波強(qiáng)度關(guān)系較好，使用平均雷達(dá)回波強(qiáng)度的方法可以在一定程度上減少降水估測(cè)的誤差。因此，為了更有效的利用三維雷達(dá)回波強(qiáng)度數(shù)據(jù)進(jìn)行定量降水估測(cè)，本文建立以自動(dòng)氣象站為中心、1 h內(nèi)逐6 min的5 km×5 km格點(diǎn)范圍的雷達(dá)回波強(qiáng)度數(shù)據(jù)（圖2）。研究中用到的自動(dòng)氣象站降水資料為整點(diǎn)的1 h雨量，與雷達(dá)資料的匹配方式為：T時(shí)刻記錄的過(guò)去1 h雨量（T-1到T的1 h累計(jì)雨量）與T-1到T的逐6 min、10個(gè)時(shí)次5 km×5 km格點(diǎn)的雷達(dá)回波反射率相匹配，形成一條數(shù)據(jù)樣本，作為降雨定量估測(cè)模型的輸入數(shù)據(jù)。本文共建立100萬(wàn)組訓(xùn)練數(shù)據(jù)集和10萬(wàn)組測(cè)試數(shù)據(jù)集。

2 基于GBDT的雷達(dá)定量估測(cè)降水算法

梯度提升決策樹(shù)是一個(gè)應(yīng)用廣泛、效果較好的監(jiān)督式機(jī)器學(xué)習(xí)模型，它于2001年由Friedman[20]提出的組合決策樹(shù)算法，由決策樹(shù)（decision tree）[21]和梯度提升（gradient boosting）[22]組合而成，提高回歸預(yù)測(cè)的精度。

GBDT是一種迭代的決策樹(shù)算法，其迭代的目的是通過(guò)計(jì)算上一次模型的負(fù)梯度來(lái)改進(jìn)模型，每一次新模型的建立都減少了上一個(gè)模型的殘差，然后在殘差減少的梯度方向上建立新的決策樹(shù)；每個(gè)新模型的輸入都是上一個(gè)模型輸出，輸出是下一個(gè)模型的輸入，預(yù)測(cè)模型的最終結(jié)果為之前所有模型結(jié)果的相加。該算法由多個(gè)弱預(yù)測(cè)模型（通常是決策樹(shù)）通過(guò)集合的形式組合產(chǎn)生預(yù)測(cè)模型，所有樹(shù)的結(jié)論加起來(lái)作為最終答案。該方法的優(yōu)點(diǎn)在于可以盡量避免過(guò)擬合現(xiàn)象，每迭代一棵樹(shù)，殘差代表的都是全局最優(yōu)的絕對(duì)方向。由于GBDT算法考慮方差和偏差之間的權(quán)衡，具有較好的泛化能力，幾乎可以解決所有的數(shù)據(jù)回歸問(wèn)題并具有高檢測(cè)精度的特性，且在非線性數(shù)據(jù)處理方面具有良好的效果，被廣泛應(yīng)用消息檢測(cè)[23]、徑流預(yù)測(cè)[24]等研究領(lǐng)域。

本文將每一個(gè)數(shù)據(jù)樣本是具有一組特征屬性的雷達(dá)回波三維網(wǎng)格數(shù)據(jù)和自動(dòng)氣象站小時(shí)雨量數(shù)據(jù)組成，定義訓(xùn)練數(shù)據(jù)集：

圖2 以自動(dòng)氣象站為中心的雷達(dá)回波反射率數(shù)據(jù) Fig. 2 Three dimensional data of radar reflectivity for GBDT as input (The gauge location is at the center of the grid)

假設(shè)通過(guò)遞歸迭代分割將特征空間（即輸入變量）劃分為m個(gè)子集，則GBDT模型表示為m個(gè)模型函數(shù)相加來(lái)預(yù)測(cè)輸出：

基于GBDT的定量降水估測(cè)QPE模型算法的詳細(xì)步驟如下：

1）使用常量值初始化模型，估計(jì)使損失函數(shù)極小化的常數(shù)值，它是只有一個(gè)根節(jié)點(diǎn)的樹(shù)，即β是一個(gè)常數(shù)值：

迭代訓(xùn)練模型，在上一次模型損失函數(shù)的梯度下降方向上建立模型，從m=1到M（M為迭代次數(shù)）

計(jì)算損失函數(shù)的負(fù)梯度在當(dāng)前模型的值，將它作為殘差ri的估計(jì)值：

2）將式（5）中得到的估計(jì)殘差ri作為輸入，擬合到一顆回歸樹(shù)，求得回歸樹(shù)的葉節(jié)點(diǎn)區(qū)

3 個(gè)例分析

利用建立的訓(xùn)練數(shù)據(jù)集進(jìn)行GBDT算法的建模和訓(xùn)練，得到基于GBDT算法的雷達(dá)定量降水估測(cè)模型，并用測(cè)試數(shù)據(jù)集進(jìn)行檢驗(yàn)和誤差分析。作為比較，本文還利用測(cè)試數(shù)據(jù)集建立了固定Z-I關(guān)系法和動(dòng)態(tài)Z-I關(guān)系法進(jìn)行定量降水估測(cè)。其中固定Z-I關(guān)系法應(yīng)用汪瑛等[8]通過(guò)統(tǒng)計(jì)歷史上強(qiáng)降水日的雷達(dá)反射率因子和1 h降水資料，得到適合廣東地區(qū)的系數(shù)A和b的值，即Z=220.8I2.23，這個(gè)統(tǒng)計(jì)結(jié)果可直接用于雷達(dá)定量降水估測(cè)，且計(jì)算速度快。動(dòng)態(tài)Z-I關(guān)系法主要利用自動(dòng)氣象站雨量計(jì)降水資料，基于Z-I關(guān)系算法，實(shí)時(shí)計(jì)算當(dāng)前時(shí)刻最優(yōu)系數(shù)A和b的值，使雷達(dá)定量估測(cè)降水與區(qū)域自動(dòng)站觀測(cè)降水的誤差最小，作為下一時(shí)次的Z-I關(guān)系計(jì)算多普勒雷達(dá)估測(cè)降水，該方法的優(yōu)點(diǎn)是前期不需要收集大量的資料樣本做統(tǒng)計(jì)準(zhǔn)備，只需要統(tǒng)計(jì)上一個(gè)小時(shí)的最優(yōu)Z-I關(guān)系。

圖3給出了2017年4月21日07時(shí)的影響廣東區(qū)域一條颮線雷達(dá)定量估測(cè)降水和06—07時(shí)的1 h雨量實(shí)況。該小時(shí)降水的主要特點(diǎn)是有一條東北—西南向覆蓋廣東中部地區(qū)、小時(shí)雨量大于10 mm的降雨帶，在廣州和茂名存在兩個(gè)強(qiáng)降水中心（小時(shí)雨量在40 mm以上），其中最大小時(shí)雨量出現(xiàn)在廣州，為65.3 mm（圖3d）。圖3a為固定Z-I關(guān)系法的定量降水估測(cè)，與自動(dòng)氣象站實(shí)況雨量相比，降水分布和落區(qū)較為一致，但是降水量級(jí)估測(cè)偏小很多，對(duì)兩個(gè)強(qiáng)降水中心的降水估測(cè)較實(shí)況比偏小20～30 mm。動(dòng)態(tài)Z-I關(guān)系法對(duì)降水的估測(cè)比固定Z-I關(guān)系法有所改善（圖3b），對(duì)位于廣州的強(qiáng)降水中心估測(cè)在40 mm以上，雖然仍比實(shí)況偏小，但在降水量級(jí)上比固定Z-I關(guān)系法已有10 mm的提高。圖3c是基于GBDT算法反演的雷達(dá)定量估測(cè)降水，不同等級(jí)降水量的面積和強(qiáng)度相對(duì)固定Z-I關(guān)系法和動(dòng)態(tài)Z-I關(guān)系法都有所改進(jìn)，尤其對(duì)位于廣州的強(qiáng)降水中心的降水量在50～55 mm，與實(shí)況更接近，改善了降水量偏小的問(wèn)題。

4 雷達(dá)定量降水結(jié)果和誤差分析

為了描述定量估測(cè)降水的誤差和精度特征，在檢驗(yàn)和分析時(shí)，使用平均誤差（ME）、相對(duì)誤差（RE）和均方根誤差（RMSE）等進(jìn)行統(tǒng)計(jì)評(píng)價(jià)，具體計(jì)算公式如下：

圖3 2017年4月21日07時(shí)雷達(dá)1 h定量降水結(jié)果和自動(dòng)氣象站實(shí)測(cè)降水（a）固定Z-I關(guān)系法；（b）動(dòng)態(tài)Z-I關(guān)系法；（c）GBDT算法；（d）雨量實(shí)況 Fig. 3 Radar-derived 1h QPE and gauge spatial distributions on 06:00 BT 21 April 2017(a) 1 h QPE of Z-I relations; (b) 1 h QPE of dynamic Z-I relations; (c)1 h QPE of GBDT scheme; (d) observation

式中，n為總樣本數(shù)，i為QPE雨量計(jì)的序列，Pi和Gi分別為雷達(dá)估測(cè)雨量和自動(dòng)氣象站實(shí)測(cè)雨量。

用平均誤差評(píng)估數(shù)據(jù)的整體差異，其值越接近0則數(shù)據(jù)整體差異越小等評(píng)價(jià)方法；使用相對(duì)誤差評(píng)估數(shù)據(jù)的可信度，其值越小數(shù)據(jù)可信度越高；用均方根差評(píng)估雷達(dá)估測(cè)雨量與自動(dòng)氣象站觀測(cè)雨量的離散程度，其值越小分布越集中；從整體上評(píng)估不同方法對(duì)降水估測(cè)精度的高低，所有誤差都在站點(diǎn)上計(jì)算。為了使評(píng)估的結(jié)果具有客觀性，剔掉沒(méi)有通過(guò)質(zhì)量控制的自動(dòng)氣象站觀測(cè)雨量數(shù)據(jù)，將1 h雨量劃分為1～5、5～10、10～20、20～30、30～40、40～50 mm和50 mm以上的7個(gè)降水量級(jí)，分析這三種方法對(duì)降水的反演效果。

圖4 給出了這三種方法的雷達(dá)定量降水和自動(dòng)氣象站觀測(cè)雨量散點(diǎn)圖。越接近直線點(diǎn)表明雷達(dá)定量降水的效果越好，高于直線說(shuō)明高估，低于直線點(diǎn)為低估。可以看出，固定Z-I關(guān)系法對(duì)雷達(dá)定量降水的反演存在低估，隨著降水量級(jí)的增大，低估程度趨于明顯，與實(shí)測(cè)降水的相關(guān)系數(shù)為0.71；動(dòng)態(tài)Z-I關(guān)系法對(duì)雷達(dá)定量降水的估測(cè)準(zhǔn)確性有所改進(jìn)，雷達(dá)估測(cè)降水和實(shí)測(cè)降水的相關(guān)系數(shù)略有提高，為0.74；應(yīng)用GBDT算法定量估測(cè)的降水與實(shí)況觀測(cè)相關(guān)性較好，相關(guān)系數(shù)提高到0.80，散點(diǎn)圖顯示集中性較好，改善了降水估測(cè)的精度，但是仍然存在低量級(jí)降水高估、高量級(jí)降水低估的現(xiàn)象。

表1是對(duì)不同量級(jí)降水的雷達(dá)定量降水誤差分析，可以看出，固定Z-I關(guān)系法對(duì)10 mm/h以下的降水存在高估，而對(duì)10 mm/h以上的降水明顯低估，對(duì)50 mm/h以上的強(qiáng)降水低估平均誤差達(dá)到34.14 mm，相對(duì)誤差為56.95%。動(dòng)態(tài)Z-I關(guān)系法對(duì)雷達(dá)定量降水估測(cè)效果優(yōu)于固定Z-I關(guān)系法，各量級(jí)的平均誤差和相對(duì)誤差較固定Z-I關(guān)系法有所減小，但仍然存在對(duì)10 mm/h以下量級(jí)降水高估，而10 mm/h以上量級(jí)降水低估的現(xiàn)象。這可能是由于動(dòng)態(tài)Z-I關(guān)系法雖然根據(jù)自動(dòng)站觀測(cè)雨量和雷達(dá)回波關(guān)系不斷調(diào)整A和b值，提高雷達(dá)定量降水估測(cè)能力，但仍然存在對(duì)同一時(shí)刻的降水類型沒(méi)有進(jìn)行分類，所有降水類型都用同一個(gè)A、b值的Z-I關(guān)系的不足。

圖4 實(shí)測(cè)雨量和雷達(dá)定量估測(cè)降水散點(diǎn)圖：（a）固定Z-I關(guān)系法；（b）動(dòng)態(tài)Z-I關(guān)系法；（c）GBDT算法 Fig. 4 The scatter diagram of hourly rainfall estimation: (a)Z-I relations; (b) dynamic Z-I relations; (c) GBDT scheme

基于GBDT算法反演的降水平均誤差和相對(duì)誤差都小于固定Z-I關(guān)系法和動(dòng)態(tài)Z-I關(guān)系法的結(jié)果，尤其對(duì)30 mm/h以上量級(jí)降水的估測(cè)，相對(duì)誤差比動(dòng)態(tài)Z-I關(guān)系法都提高了5%以上，對(duì)降水反演的系統(tǒng)誤差更小，改善了降水估測(cè)的精度。RMSE的含義是雷達(dá)估測(cè)雨量與自動(dòng)氣象觀測(cè)的降雨量之間的離散程度，對(duì)5～20 mm/h量級(jí)的降水估測(cè)，動(dòng)態(tài)Z-I關(guān)系法RMSE最小，但對(duì)于20 mm/h以上的強(qiáng)降水，GBDT算法估測(cè)降雨量的RMSE都小于固定Z-I關(guān)系法和動(dòng)態(tài)Z-I關(guān)系法，表明GBDT算法對(duì)強(qiáng)降水的雨量估測(cè)擬合度較高。

表1 不同量級(jí)降水的雷達(dá)定量降水誤差分析 Table 1 Deviation statistics of different rainfall intensity categories

5 結(jié)論

利用廣東省新一代多普勒氣象雷達(dá)三維拼圖數(shù)據(jù)和地面自動(dòng)氣象站雨量資料，建立了一種基于梯度提升決策樹(shù)GBDT算法的雷達(dá)定量估測(cè)降雨模型，通過(guò)與固定Z-I關(guān)系法和動(dòng)態(tài)Z-I關(guān)系法進(jìn)行誤差分析和對(duì)比，結(jié)果表明：

1）固定Z-I關(guān)系法對(duì)小于10 mm/h的降水量級(jí)存在高估，但對(duì)10 mm以上量級(jí)的降水存在低估，隨著降水量級(jí)的增大，低估現(xiàn)象越明顯；動(dòng)態(tài)Z-I關(guān)系法優(yōu)于固定Z-I關(guān)系法，對(duì)雷達(dá)定量降水的估測(cè)準(zhǔn)確性有所改進(jìn)。

2）通過(guò)對(duì)平均誤差、均方根誤差和相對(duì)誤差等的統(tǒng)計(jì)結(jié)果表明，基于梯度提升決策樹(shù)GBDT算法估測(cè)的降雨量精度要高于固定Z-I關(guān)系法和動(dòng)態(tài)Z-I關(guān)系法，特別是對(duì)30 mm/h以上的強(qiáng)降水的估測(cè)，GBDT算法的強(qiáng)降水中心和分布與自動(dòng)氣象站觀測(cè)較為吻合，對(duì)短時(shí)強(qiáng)降水低估現(xiàn)象有一定的改進(jìn)作用。