朱巖 翟丹華 吳志鵬 張焱

(重慶市氣象臺(tái),重慶 401147)

引言

短時(shí)強(qiáng)降水是重慶地區(qū)最嚴(yán)重的災(zāi)害性天氣之一[1]，由于降水效率高、雨量大，局地性、突發(fā)性強(qiáng)，且極易引發(fā)山洪、泥石流和內(nèi)澇等次生災(zāi)害，是汛期短臨預(yù)報(bào)預(yù)警的重點(diǎn)和難點(diǎn)，因此，當(dāng)前亟需更多預(yù)報(bào)產(chǎn)品為短時(shí)強(qiáng)降水的預(yù)報(bào)提供技術(shù)支撐。前人針對短時(shí)強(qiáng)降水的氣候統(tǒng)計(jì)、個(gè)例分析和復(fù)雜地形影響[2]等方面的特征條件已經(jīng)做了很多有益的探索。李強(qiáng)等[3]分析了川渝地區(qū)汛期短時(shí)強(qiáng)降水事件日變化特征，指出其多發(fā)并加強(qiáng)于夜間，且歷時(shí)較長，極值降水多持續(xù)2 h等特征。袁晨[4]、周北平[5]等人對川渝地區(qū)短時(shí)強(qiáng)降水時(shí)空分布特征和監(jiān)測預(yù)警技術(shù)也有研究。

短時(shí)強(qiáng)降水臨近預(yù)報(bào)方法主要包括基于衛(wèi)星云參數(shù)預(yù)報(bào)[6]、雷達(dá)定量估測降水、基于配料法的統(tǒng)計(jì)回歸預(yù)報(bào)[7]等傳統(tǒng)方法以及新興的人工智能預(yù)報(bào)?；谂淞戏ê徒y(tǒng)計(jì)回歸方法的短時(shí)強(qiáng)降水預(yù)報(bào)技術(shù)已得到長足發(fā)展，并取得了一系列較好成果[8-14]。但是隨著氣象大數(shù)據(jù)的累積以及預(yù)報(bào)業(yè)務(wù)對準(zhǔn)確率和精細(xì)化程度不斷提高的要求，這些傳統(tǒng)方法在應(yīng)對高時(shí)空分辨率、多變量和高度非線性等挑戰(zhàn)上逐漸顯露出不足。

隨著人工智能在各行各業(yè)突飛猛進(jìn)的發(fā)展，氣象預(yù)報(bào)領(lǐng)域也逐漸引入機(jī)器學(xué)習(xí)算法，并在諸如污染物、氣溫和能見度預(yù)報(bào)、強(qiáng)對流天氣識別、回波預(yù)報(bào)等方面取得了一些成果。在強(qiáng)降水預(yù)報(bào)方面，Gagne等[15]建立了基于邏輯回歸和隨機(jī)森林的定量降水校正概率預(yù)報(bào)。Ma Liang等[16]基于葵花衛(wèi)星資料和數(shù)字高程資料，利用梯度提升樹(GBDT)算法建立了降水落區(qū)預(yù)報(bào)模型。劉媛媛等[17]采用動(dòng)態(tài)聚類算法研究了北京城區(qū)近12年來短歷時(shí)暴雨的時(shí)空分布規(guī)律。楊通曉等[18]采用粒子群優(yōu)化算法，結(jié)合支持向量機(jī)開發(fā)了雙偏振雷達(dá)對流降水類型識別方法?；跈C(jī)器學(xué)習(xí)算法、針對短時(shí)強(qiáng)降水的預(yù)報(bào)技術(shù)研究也陸續(xù)見諸學(xué)界。路志英等[19]建立了一個(gè)深度信念網(wǎng)絡(luò)(DBNs)，利用地面加密觀測數(shù)據(jù)對天津市短時(shí)強(qiáng)降水進(jìn)行預(yù)報(bào)，檢驗(yàn)結(jié)果表明DBNs預(yù)報(bào)明顯優(yōu)于支持向量機(jī)(SVM)和邏輯回歸(LR)。白曉平等[20]分別運(yùn)用改進(jìn)的二元Logistic回歸法和綜合多指標(biāo)疊加法建立了西北地區(qū)東部的短時(shí)強(qiáng)降水預(yù)報(bào)模型，二元Logistic回歸法對獨(dú)立樣本預(yù)報(bào)的TS評分高達(dá)46.6%。鐘海燕[21]引入支持向量機(jī)(SVM)、梯度提升樹(GBDT)、極限提升樹(eXtreme Gradient Boosting, XGBoost)3種機(jī)器學(xué)習(xí)方法應(yīng)用于雷達(dá)降雨產(chǎn)品的臨近降雨預(yù)報(bào)，實(shí)驗(yàn)表明，基于SVM和GBDT的方法都具有較好效果，而XGBoost方法與PPLK方法相結(jié)合的預(yù)報(bào)性能在所有實(shí)驗(yàn)成員中效果最好。孫俊奎等[22]基于概率神經(jīng)網(wǎng)絡(luò)(PNN)、支持向量機(jī)(SVM)和邏輯回歸3種機(jī)器學(xué)習(xí)算法對石林地區(qū)的逐3 h間隔降水量的8個(gè)等級進(jìn)行回歸建模，業(yè)務(wù)測試表明，PNN和SVM模型優(yōu)于邏輯回歸模型，TS在45%左右，中雨以上量級預(yù)報(bào)，3種模型的TS評分達(dá)28%。

目前，機(jī)器學(xué)習(xí)算法在短時(shí)強(qiáng)降水預(yù)報(bào)中的應(yīng)用研究為數(shù)尚少，而重慶也缺乏本地化的短時(shí)強(qiáng)降水客觀預(yù)報(bào)方法及產(chǎn)品，因此本文將針對重慶地區(qū)短時(shí)強(qiáng)降水，基于二源或三源融合的逐小時(shí)格點(diǎn)降水資料和EC細(xì)網(wǎng)格再分析資料，通過箱線圖差異指數(shù)尋找預(yù)報(bào)變量閾值并建立消空規(guī)則，結(jié)合K均值聚類算法和Relief算法重建訓(xùn)練集并優(yōu)選建模變量，最后建立基于Xgboost算法的短時(shí)強(qiáng)降水客觀預(yù)報(bào)模型，為重慶和周邊地區(qū)短時(shí)強(qiáng)降水預(yù)報(bào)提供技術(shù)支撐和客觀預(yù)報(bào)產(chǎn)品。

本文基于配料法選取了兩方面預(yù)報(bào)變量，一是短時(shí)強(qiáng)降水短臨預(yù)報(bào)業(yè)務(wù)中常用的基本診斷量，如可降水量、K指數(shù)和水汽通量散度等，二是由這些基本量衍生得到的具有綜合指示意義的復(fù)合因子，如濕位渦、濕熱力平流等，它們在高守亭等[23-24]對暴雨預(yù)報(bào)的理論和預(yù)報(bào)系統(tǒng)研發(fā)中得到了檢驗(yàn)和推廣。

1 資料和方法

1.1 觀測和預(yù)報(bào)資料

實(shí)況降水?dāng)?shù)據(jù)取自中國氣象局國家氣象信息中心的高分辨率融合降水產(chǎn)品[25]?？紤]到重慶地區(qū)短時(shí)強(qiáng)降水主要集中在每年的5—9月，因此選取的數(shù)據(jù)時(shí)段為2011年至2015年的5—9月。其中2011年至2014年數(shù)據(jù)為地面-衛(wèi)星二源融合產(chǎn)品，水平分辨率0.1°×0.1°，2015年數(shù)據(jù)為地面-衛(wèi)星-雷達(dá)三源融合產(chǎn)品，水平分辨率0.05°×0.05°。為統(tǒng)一分辨率，從0.05°×0.05°降水場中抽稀取出與0.1°×0.1°格點(diǎn)場相同的格點(diǎn)，并截取重慶和周邊地區(qū)(28°～32.5°N，105°～110.5°E)作為研究區(qū)域，共2475個(gè)格點(diǎn)。根據(jù)中央氣象臺(tái)標(biāo)準(zhǔn)，短時(shí)強(qiáng)降水閾值為1 h降水量≥20 mm，為方便數(shù)據(jù)處理和計(jì)算，這里取逐小時(shí)的整點(diǎn)值。由于本文預(yù)報(bào)對象為短時(shí)強(qiáng)降水是否發(fā)生，因此樣本的標(biāo)記采用二分類，達(dá)到閾值的樣本記為1，未達(dá)到的記為0，文中分別稱為正例和負(fù)例。實(shí)際業(yè)務(wù)中，在不利于短時(shí)強(qiáng)降水的大尺度天氣形勢下，中小尺度局地環(huán)流或者較好的環(huán)境場配置也會(huì)引發(fā)極少數(shù)孤立的短時(shí)強(qiáng)降水。以達(dá)標(biāo)站數(shù)是否滿足總格點(diǎn)數(shù)的1%為判斷標(biāo)準(zhǔn)，剔除了低于1%的時(shí)次，減少了短時(shí)強(qiáng)降水比例過低時(shí)次對建模過程的干擾。

預(yù)報(bào)場采用EC細(xì)網(wǎng)格(0.25°×0.25°)2011—2015年的5—9月逐日02:00、08:00、14:00、20:00(北京時(shí)，下同)起報(bào)的再分析資料，包括位勢高度、溫度、相對濕度、風(fēng)速、風(fēng)向等基本氣象要素用于計(jì)算物理量場，以及相同網(wǎng)格分辨率的逐小時(shí)總降水量再分析場作為基準(zhǔn)，用于對比評估本方法相對模式預(yù)報(bào)的效果。另外，EC資料中不包含地形數(shù)據(jù)，而重慶地區(qū)山脈縱橫，地形復(fù)雜，地形對強(qiáng)降水的觸發(fā)作用應(yīng)當(dāng)考慮在內(nèi)。因此，地形數(shù)據(jù)采用重慶市氣象局與美國Oklahoma大學(xué)共同研發(fā)的中尺度數(shù)值模式輸出的水平分辨率3 km地形高度值。空間上，采用雙線性插值將EC細(xì)網(wǎng)格上的要素插值到0.1°格點(diǎn)場上；時(shí)間上，以距離當(dāng)前降水時(shí)次最近的前一再分析場時(shí)次作為相匹配的預(yù)報(bào)場，經(jīng)過整合后形成一套完整統(tǒng)一的樣本集。其中2011—2014年樣本作為訓(xùn)練集，2015年樣本作為測試集。除地形高度以外所有變量的計(jì)算均基于EC細(xì)網(wǎng)格高空形勢場的基本要素，中英文變量名及其單位見表1。

表1 預(yù)報(bào)變量Ibd值

1.2 樣本中變量和預(yù)報(bào)對象的三維分布特征

為初步觀察預(yù)報(bào)變量對短時(shí)強(qiáng)降水的區(qū)分程度，分別從訓(xùn)練集的基本和衍生診斷量中選取K指數(shù)、整層可降水量、850 hPa水汽通量散度和對流有效位能、濕熱力平流參數(shù)、風(fēng)暴強(qiáng)度指數(shù)兩組變量，繪制正、負(fù)例在變量構(gòu)成的特征空間中的分布。由于各變量量級差異以及自身變化幅度均較大，為方便對比數(shù)據(jù)和可視化的美觀，將以上變量標(biāo)準(zhǔn)化為無量綱數(shù)值。

如圖1所示，選取的變量是具有一定識別作用的，部分正例處于負(fù)例之外，但是也存在很多重疊、交叉樣本，且負(fù)例總量遠(yuǎn)遠(yuǎn)大于正例。因此，需要根據(jù)變量對正負(fù)例的區(qū)分度初步篩選樣本消除負(fù)例，并選擇預(yù)測能力較強(qiáng)的變量進(jìn)行建模，以增強(qiáng)分類效果。另外，嚴(yán)重的類別不平衡性也要求在建模和預(yù)報(bào)時(shí)必須根據(jù)樣本權(quán)重在重采樣、誤分代價(jià)、判定閾值等至少某一環(huán)節(jié)中采取措施，以免模型的響應(yīng)能力被多數(shù)類樣本“淹沒”。

圖1 2011—2014年樣本集變量和預(yù)報(bào)對象的三維分布:(a)K指數(shù)(K),整層可降水量(PW)和850 hPa水汽通量散度(QF850),(b)對流有效位能(CAPE),濕熱力平流參數(shù)(GMTP)和風(fēng)暴強(qiáng)度指數(shù)(SSI)(紅點(diǎn)為正例，綠點(diǎn)為負(fù)例)

1.3 方法

1.3.1 閾值法和消空規(guī)則

由于樣本集存在的嚴(yán)重類別不平衡性會(huì)削弱分類效果，因此需要在建模前尋找特征變量進(jìn)行消空，通過剔除負(fù)樣本在一定程度上削弱不平衡性。為考察所有預(yù)報(bào)變量對短時(shí)強(qiáng)降水的區(qū)分度，引入Fu等[26]定義的箱線圖差異指數(shù)(box difference index,Ibd)，對于變量的選取和預(yù)報(bào)前的消空具有很好作用：

(1)

式(1)中，M1表示短時(shí)強(qiáng)降水發(fā)生時(shí)某變量的均值，M0為無短時(shí)強(qiáng)降水時(shí)的相應(yīng)值，σ1和σ0分別為該變量在兩種情況下相應(yīng)的標(biāo)準(zhǔn)差。如果事件發(fā)生和不發(fā)生時(shí)變量的平均值差異大而總方差小，也即Ibd較大，則該變量對事件是否發(fā)生的區(qū)分度高，適于在預(yù)報(bào)前初步消空潛勢低的樣本。

各預(yù)報(bào)變量的Ibd如表1所示，左右兩列已按照Ibd絕對值降序排列。在與短時(shí)強(qiáng)降水呈反相關(guān)的變量中，地形高度、濕熱力平流參數(shù)、700 hPa散度和水汽通量散度的Ibd相對較高，絕對值超過0.1；在與短時(shí)強(qiáng)降水呈正相關(guān)變量中，850 hPa渦度、修正K指數(shù)和K指數(shù)等變量的Ibd相對較高。這些具有相對大值Ibd的變量將被賦予更高權(quán)重加入消空過程，以體現(xiàn)其對短時(shí)強(qiáng)降水更強(qiáng)的區(qū)分能力。另外，衍生量雖然綜合反映了動(dòng)、熱力和水汽條件，具有更全面的指示意義，但I(xiàn)bd都相對較低。

如圖2所示，以K指數(shù)和700 hPa散度為例，訓(xùn)練集負(fù)例對應(yīng)的兩個(gè)變量值域分別為12～48 ℃、-10～32 s-1，且該值域完全包含了正例的值域。負(fù)例對應(yīng)的變量奇異值可能為模式的錯(cuò)誤預(yù)報(bào)，在實(shí)際情況中極少出現(xiàn)。因此如果分別選用訓(xùn)練集中的最小、最大變量值(通常為負(fù)例的最值)作為短時(shí)強(qiáng)降水發(fā)生的臨界閾值，而無視利于降水發(fā)生的最小、最大變量值(正例的最值)，就會(huì)引入大量負(fù)例從而造成空報(bào)過多。為此引入消空所依據(jù)的閾值法：選取短時(shí)強(qiáng)降水發(fā)生時(shí)某變量剔除異常值后的最小值和最大值作為閾值。異常值的定義為：令Q3，Q1分別對應(yīng)某診斷量箱的上下邊界值，Q3-Q1則為上下四分位距，當(dāng)某值小于最小閾值Q1-1.5(Q3-Q1)或大于最大閾值Q3+1.5(Q3-Q1)時(shí)，則為異常值。根據(jù)以上規(guī)則，如圖2所示，K指數(shù)最小、最大閾值分別為32 ℃、44 ℃，700 hPa散度最小、最大閾值分別為-4 s-1、4 s-1。閾值法消空會(huì)以剔除一小部分正例為代價(jià)，大幅降低空報(bào)率，使整體TS評分得到提高。最后，即使K指數(shù)的Ibd在所有正相關(guān)變量中名列第3，已具有較好的區(qū)分能力，正、負(fù)例的K指數(shù)箱線仍非常接近，可見滿足單一或者少數(shù)幾個(gè)變量的閾值仍然難以區(qū)分兩類樣本并有效消空，因此應(yīng)建立多變量閾值規(guī)則判斷降水潛勢，剔除潛勢較弱的格點(diǎn)或滿足潛勢條件的格點(diǎn)數(shù)占總體格點(diǎn)數(shù)比例過低的時(shí)次。應(yīng)該注意到，由于不同地形高度處均有發(fā)生短時(shí)強(qiáng)降水可能性，如果只預(yù)報(bào)滿足某高度閾值的格點(diǎn)，則會(huì)造成低于此高度閾值格點(diǎn)從來不進(jìn)入預(yù)報(bào)的不合理現(xiàn)象，故地形高度不進(jìn)入消空環(huán)節(jié)。

圖2 訓(xùn)練集短時(shí)強(qiáng)降水發(fā)生(1)和不發(fā)生(0)時(shí)K指數(shù)(a)和700 hPa散度(b)分布箱線

具體消空規(guī)則為：對于某一時(shí)次預(yù)報(bào)場，計(jì)算每個(gè)格點(diǎn)的多變量閾值條件加權(quán)和。對每個(gè)格點(diǎn)有：

(2)

式(2)中，|IBD|為各變量的IBD值絕對值，Bool為1或0的布爾值，表示是否滿足該變量的閾值，n為變量數(shù)，S為所有變量的IBD絕對值在該格點(diǎn)上的加權(quán)和。判斷總體潛勢是否達(dá)到消空標(biāo)準(zhǔn)，即S值非零的格點(diǎn)數(shù)是否達(dá)到總格點(diǎn)的1%，若未達(dá)到，則認(rèn)為該時(shí)次發(fā)生短時(shí)強(qiáng)降水潛勢很低，不進(jìn)入模型預(yù)報(bào)，所有格點(diǎn)預(yù)報(bào)為無短時(shí)強(qiáng)降水發(fā)生。表2為訓(xùn)練集和測試集在經(jīng)過消空前后的正、負(fù)例數(shù)量和比例。消空后，正例占比都有一定升高，但相對負(fù)例仍然很低，類別不平衡現(xiàn)象依然明顯。

表2 訓(xùn)練集和測試集消空前后正、負(fù)例數(shù)量和比例變化

1.3.2 類別平衡和特征選擇

初選的預(yù)報(bào)變量為半經(jīng)驗(yàn)性，且各因子之間存在一定程度共線性，可能會(huì)對機(jī)器學(xué)習(xí)過程帶來負(fù)擔(dān)。由Ibd值可見，若干變量對短時(shí)強(qiáng)降水的區(qū)分度很低，其預(yù)測能力可能也較低，需要進(jìn)一步篩選預(yù)報(bào)變量。而且，經(jīng)過初步消空以后正例所占比例仍然很低，兩類樣本依然很不平衡。處理類別不平衡數(shù)據(jù)分類問題主要有兩大途徑：在算法層面，采用集成學(xué)習(xí)穩(wěn)定整個(gè)分類器的預(yù)測性能并提高泛化能力，或者引入懲罰機(jī)制，根據(jù)不同類別的風(fēng)險(xiǎn)偏好對錯(cuò)誤分類設(shè)置相應(yīng)的代價(jià)函數(shù)加深學(xué)習(xí)“記憶”；在數(shù)據(jù)層面，采用過采樣、欠采樣、人工合成新樣本等方式平衡類別或者設(shè)置非常規(guī)的判定閾值。

(3)

其中，diff(xi(j),yi(j))為樣本與其臨近樣本在第j維特征上的距離函數(shù)。M為隨機(jī)抽樣次數(shù)。

當(dāng)?shù)趈維特征為非數(shù)值型變量時(shí)，距離函數(shù)為：

(4)

當(dāng)?shù)趈維特征為數(shù)值型變量時(shí)，距離函數(shù)為：

(5)

其中，max(j)、min(j)分別為第j維特征數(shù)值的最大、最小值。

由式(3)可見，對于第j維特征，當(dāng)樣本到不同類最近樣本的距離大于同類最近樣本的距離時(shí)，該樣本在第j維特征上有利于分類器分類，因此w(j)累加正數(shù)，否則該特征不利于分類，w(j)累加負(fù)數(shù)。

然而在不平衡數(shù)據(jù)集上，多數(shù)類樣本被選中的概率遠(yuǎn)遠(yuǎn)高于少數(shù)類樣本，等比例更新權(quán)重可能使特征權(quán)重偽偏大從而影響分類效果。針對此問題，改進(jìn)了Relief算法以面向不平衡數(shù)據(jù)。首先使用K均值算法將多數(shù)類樣本聚類為q類，然后分別在q個(gè)子類和少數(shù)類，總共q+1類樣本中按比例釆樣，形成兩大類別大致平衡的訓(xùn)練集。同時(shí)引入樣本權(quán)重到Relief算法中，判斷隨機(jī)選擇的樣本，若為少數(shù)類樣本，更新特征權(quán)重時(shí)乘以大于1的因子，反之則乘以小于1的因子。

當(dāng)隨機(jī)選擇的樣本屬于少數(shù)類樣本時(shí)，特征的權(quán)重更新公式為：

(6)

屬于多數(shù)類樣本時(shí)：

(7)

其中，|S|、|L|分別表示少數(shù)類和多數(shù)類樣本數(shù)。決定特征去留的閾值由下式計(jì)算：

(8)

式中,M是隨機(jī)抽樣次數(shù)：3678，α是第1類錯(cuò)誤的概率，取0.05，得τ≈0.073。將經(jīng)過以上樣本平衡和特征篩選步驟后形成的樣本作為最終訓(xùn)練集。

各個(gè)變量的平均權(quán)重如圖3所示，平均權(quán)重大于閾值的變量在圖中以黑色柱體顯示，有風(fēng)暴強(qiáng)度指數(shù)、對流有效位能、總指數(shù)、修正K指數(shù)、整層可降水量、抬升指數(shù)、500 hPa渦度平流和700 hPa渦度。采用這些入選變量建模。和Ibd值比較可見，部分高Ibd值的變量平均權(quán)重較低，少數(shù)高Ibd值變量仍具有較高權(quán)重，如總指數(shù)、修正K指數(shù)和抬升指數(shù)，因此Ibd不能完全代表變量的預(yù)測能力。

圖3 各預(yù)報(bào)變量平均權(quán)重(黑色為大于閾值的變量)

1.3.3 評價(jià)指標(biāo)

本文結(jié)合機(jī)器學(xué)習(xí)領(lǐng)域常用的AUC值以及氣象預(yù)報(bào)領(lǐng)域的TS評分、空報(bào)率、漏報(bào)率作為評價(jià)指標(biāo)。

ROC曲線(受試者工作特征曲線)是指在特定刺激條件下，以被試對象在不同判斷標(biāo)準(zhǔn)下所得的空報(bào)概率為橫坐標(biāo)，以命中概率為縱坐標(biāo)，連接各點(diǎn)而成的連線。AUC(Area Under Curve)是衡量二分類模型優(yōu)劣的評價(jià)指標(biāo)，為ROC曲線下方與坐標(biāo)軸圍成的面積，取值范圍在[0.5,1]，越接近1，分類器越完美，越接近0.5，分類器越接近隨機(jī)猜測。AUC同時(shí)考慮了分類器對于正例和負(fù)例的分類能力，在樣本不平衡的情況下，依然能夠?qū)Ψ诸惼髯鞒龊侠淼脑u價(jià)，但AUC反應(yīng)了太過籠統(tǒng)的信息,無法反應(yīng)召回率、精確率等在實(shí)際業(yè)務(wù)中經(jīng)常關(guān)心的指標(biāo)。

有鑒于此，根據(jù)TS評分、命中率POD和空報(bào)率FAR的定義，結(jié)合機(jī)器學(xué)習(xí)領(lǐng)域常用評價(jià)指標(biāo)，得到三者的計(jì)算公式為：

TS=NA/(NA+NB+NC)

(9)

POD=NA/(NA+NB)

(10)

FAR=NC/(NA+NC)

(11)

其中，NA，ND，NC，NB為模型對二分類問題正確和錯(cuò)誤判斷的樣本數(shù)量，其意義見表3。

表3 機(jī)器學(xué)習(xí)性能指標(biāo)的意義

值得注意的是，AUC表示模型總體上對所有類別的預(yù)測性能。當(dāng)類別不均衡時(shí)，容易因?yàn)楦叩腁UC而忽略實(shí)際上對少數(shù)類預(yù)報(bào)效果并不理想的情況，因此命中率和空報(bào)率對于類別不平衡問題是更清晰的衡量指標(biāo)。

1.4 Xgboost建模方法

集成學(xué)習(xí)通過構(gòu)建并結(jié)合多個(gè)學(xué)習(xí)器來完成學(xué)習(xí)任務(wù)，可獲得比單一學(xué)習(xí)器更好的泛化性能。作為Boosting集成學(xué)習(xí)算法家族中的一員，Xgboost是一個(gè)樹集成模型，將K個(gè)CART回歸樹的結(jié)果進(jìn)行求和，作為最終的預(yù)測值：

(12)

式中，Xi為第i個(gè)樣本，f(x)為單個(gè)樹的結(jié)構(gòu)和葉節(jié)點(diǎn)權(quán)重，Φ為所有k個(gè)樹的集成。不同于傳統(tǒng)集成決策樹算法，Xgboost能夠在節(jié)點(diǎn)內(nèi)選擇最佳分裂點(diǎn)，候選分裂點(diǎn)計(jì)算增益用多線程并行，訓(xùn)練速度快。其代價(jià)函數(shù)為：

(13)

(14)

Xgboost有包括正則化項(xiàng)、學(xué)習(xí)率和決策樹數(shù)量和樹結(jié)構(gòu)屬性等眾多超參數(shù)。超參數(shù)定義了模型的復(fù)雜度或?qū)W習(xí)能力等特定基本屬性，是在開始學(xué)習(xí)過程之前需要確定的參數(shù)，而不是在學(xué)習(xí)過程中由算法習(xí)得的參數(shù)，如權(quán)重和偏置。調(diào)節(jié)超參數(shù)的意義在于最小化期望風(fēng)險(xiǎn)，使模型優(yōu)化度和模型復(fù)雜度達(dá)到平衡，盡可能同時(shí)避免欠擬合和過擬合。網(wǎng)格搜索是應(yīng)用最廣泛的建立在交叉驗(yàn)證基礎(chǔ)上的超參數(shù)搜索算法，這種窮舉式調(diào)參算法通過循環(huán)遍歷嘗試每種參數(shù)組合的可能性，找出表現(xiàn)最好的組合，能夠找到全局最大或最小值。過程中采用了5折交叉驗(yàn)證，也即將訓(xùn)練集5等分，取其中一份為驗(yàn)證集，其余4份為新訓(xùn)練集，經(jīng)過5次在不同驗(yàn)證集上的測試，取最優(yōu)結(jié)果所對應(yīng)的超參數(shù)組合。最佳關(guān)鍵超參數(shù)組合如表4所示：

表4 Xgboost最優(yōu)超參數(shù)配置

2 試驗(yàn)結(jié)果

2.1 總體檢驗(yàn)結(jié)果

將2015年獨(dú)立樣本經(jīng)過以上閾值消空和變量選擇步驟后投入模型預(yù)報(bào)。模型可以輸出概率預(yù)報(bào)也可以輸出確定性預(yù)報(bào)。經(jīng)過聚類后的新訓(xùn)練集類別大致平衡，Xgboost的缺省樣本不平衡度為1，也即類別平衡，而測試集中仍為類別不平衡，短時(shí)強(qiáng)降水出現(xiàn)概率小于0.5，如果采用模型在類別平衡時(shí)默認(rèn)的0.5概率值，確定性預(yù)報(bào)會(huì)產(chǎn)生大量漏報(bào)，因此采用了不同的概率閾值生成相應(yīng)的確定性預(yù)報(bào)以觀察效果。

如圖4所示，AUC由概率預(yù)報(bào)產(chǎn)生的概率值本身決定而不隨概率閾值變化，因此保持0.92，表明模型的分類性能總體上較好。但對于類別不平衡的預(yù)測問題，還需要考察模型對各個(gè)類別的預(yù)測準(zhǔn)確率。隨著概率值從0.01開始增加，TS評分上升，POD和FAR下降，在閾值為0.1左右時(shí)三者達(dá)到穩(wěn)定，分別為0.30、0.88和0.69。在閾值從0.1增至0.35的過程中，由于除個(gè)別樣本的概率值接近0.35以外，其余均小于0.1，各項(xiàng)指標(biāo)因此不隨概率閾值的增加而變化。當(dāng)閾值超過0.35后，由于模型未能識別出任何正例(NA、NC均為0)而導(dǎo)致TS和POD陡降至0，F(xiàn)AR為無意義的除零數(shù)(圖中置為0)，因此最佳的概率閾值為0.1，對應(yīng)TS為0.3。實(shí)際應(yīng)用中，可根據(jù)對空報(bào)和漏報(bào)率的敏感程度調(diào)節(jié)閾值得到用戶自定義的確定性預(yù)報(bào)結(jié)果。在相同測試集上，EC未能體現(xiàn)出實(shí)況的任何短時(shí)強(qiáng)降水事件，且全部偏離實(shí)況，因此其TS和POD評分為0，空報(bào)、漏報(bào)率為1(圖略)。模型相對EC對于2015年測試集上的短時(shí)強(qiáng)降水預(yù)報(bào)具有一定優(yōu)勢。

圖4 各項(xiàng)檢驗(yàn)指標(biāo)在不同概率閾值上的分布

2.2 獨(dú)立樣本檢驗(yàn)期間兩次短時(shí)強(qiáng)降水預(yù)報(bào)分析

2.2.1 渝西南短時(shí)強(qiáng)降水過程分析

2015年6月28日夜間至30日，受高空波動(dòng)槽、中低層低渦切變和低層暖平流影響，重慶長江沿線以北地區(qū)大雨到暴雨，西部區(qū)縣的部分鄉(xiāng)鎮(zhèn)達(dá)大暴雨，并有短時(shí)強(qiáng)降水等強(qiáng)對流天氣，24 h累積雨量如圖5a所示。

圖5 重慶市2015年6月28日20:00至29日20:00(a) 和7月21日08:00至22日08:00(b)實(shí)況雨量

如圖6a、c所示，雖然EC再分析場在6月29日02:00—05:00的累積降水量主雨帶形態(tài)與實(shí)況較吻合，重慶西北部預(yù)報(bào)的3～10 mm降水落區(qū)(紅色實(shí)線圈所示)對實(shí)況的相應(yīng)強(qiáng)降水區(qū)域有所反映，但預(yù)報(bào)較實(shí)況相比明顯偏弱，主雨帶中未預(yù)報(bào)10 mm以上落區(qū)，對于川北地區(qū)孤立的20 mm以上強(qiáng)降水中心(紅色虛線圈所示)的預(yù)報(bào)位置也有明顯偏西，預(yù)報(bào)效果不佳。Xgboost模型于29日02:00起報(bào)的短時(shí)強(qiáng)降水概率高值區(qū)(圖6b)與實(shí)況的吻合度大大提升，除少數(shù)地區(qū)的空、漏報(bào)外，Xgboost預(yù)報(bào)不僅抓住了從重慶西北部到重慶東北部的西南—東北向主雨帶形態(tài)特征，還對川北地區(qū)的短時(shí)強(qiáng)降水落區(qū)有所反應(yīng)(如圖中紅色虛線圈所示)，而全球模式往往很難預(yù)報(bào)出此類相對較弱的降雨中心或次雨帶。02:00—03:00(圖6d)，遂寧、潼南已產(chǎn)生短時(shí)強(qiáng)降水，且中心強(qiáng)度在30 mm以上，從南充到城口為一條斷裂為南北兩段的強(qiáng)降水帶，隨后遂寧—潼南強(qiáng)降水中心向東、向南發(fā)展進(jìn)入合川、銅梁，南充—城口強(qiáng)降水帶在保持基本形態(tài)的前提下向東發(fā)展，有若干小中心生消演變(圖6e)。川北地區(qū)在04:00—05:00新生強(qiáng)降水中心(圖6f)。Xgboost預(yù)報(bào)的概率高值區(qū)基本包含了這些時(shí)段的強(qiáng)降水落區(qū)，且其西南-東北向的大值區(qū)中具有兩條主線(如圖中紅色實(shí)曲線所示)，分別與03:00—05:00強(qiáng)降水發(fā)展演變所形成的兩條主雨帶(如圖中黑色實(shí)曲線所示)形成對應(yīng)。到29日下午，如圖7所示，Xgboost預(yù)報(bào)仍好于EC再分析場，對雨帶強(qiáng)度和形態(tài)的刻畫均更準(zhǔn)確。Xgboost的高概率區(qū)不僅分布在重慶西部的重慶主城、銅梁、璧山等已經(jīng)發(fā)生了短時(shí)強(qiáng)降水的地區(qū)，在廣安和江津也有分布(分別為紅色實(shí)線、虛線圈所示)。隨著降水系統(tǒng)的移動(dòng)和演變，到20:00(圖7f)，18:00初生于廣安的降水增強(qiáng)到20 mm以上，而江津也產(chǎn)生短時(shí)強(qiáng)降水，與Xgboost在14:00預(yù)報(bào)的概率高值區(qū)吻合。此過程的兩個(gè)時(shí)段中，如圖10所示，Xgboost預(yù)報(bào)的TS、POD、FAR分別在0.2～0.4、0.6～0.9和0.6～0.8之間，EC的TS為0，未在圖中顯示。因此Xgboost對于此次過程兩個(gè)時(shí)段的短時(shí)強(qiáng)降水預(yù)報(bào)好于EC。

圖6 EC細(xì)網(wǎng)格再分析場的6月29日02:00—05:00降水量(a)、Xgboost模型29日02:00起報(bào)的短時(shí)強(qiáng)降水客觀概率預(yù)報(bào)(b)和29日02:00—05:00實(shí)況降水量(c)以及02:00—03:00(d)、03:00—04:00(e)、04:00—05:00(f)小時(shí)降水量

圖7 EC細(xì)網(wǎng)格再分析場的6月29日17:00—20:00降水量(a)、Xgboost模型29日14:00起報(bào)的短時(shí)強(qiáng)降水客觀概率預(yù)報(bào)(b)和29日17:00—20:00實(shí)況降水量(c)以及17:00—18:00(d)、18:00—19:00(e)、19:00—20:00(f)小時(shí)降水量

2.2.2 渝西短時(shí)強(qiáng)降水過程分析

2015年7月21日至22日，受高空槽冷平流和中低層低渦切變影響，中西部和東南部地區(qū)及東北部偏南地區(qū)普降大雨到暴雨，局部大暴雨，并伴有短時(shí)強(qiáng)降水等強(qiáng)對流天氣。24 h累積雨量如圖5b所示。

如圖8a、c所示，EC預(yù)報(bào)較上一次過程更好，雖然總體上仍然偏弱，但在重慶西部與四川交界地區(qū)預(yù)報(bào)了20 mm以上強(qiáng)降水，降水落區(qū)的形態(tài)與實(shí)況在一定程度上吻合。14:00起報(bào)的Xgboost則不僅預(yù)報(bào)了重慶偏西地區(qū)的概率高值區(qū)，對應(yīng)著實(shí)況3 h累積雨量在50 mm以上的強(qiáng)降水中心，在合川和永川分別也有概率高值區(qū)(分別為紅色實(shí)線、虛線圈所示)，從圖8d～f的逐時(shí)降水量演變可見，雨團(tuán)東移發(fā)展并逐漸體現(xiàn)較清晰的人字形切變形態(tài)，19:00—20:00在合川和永川也產(chǎn)生了短時(shí)強(qiáng)降水。到21日夜間，如圖9所示，EC在重慶西南部預(yù)報(bào)的強(qiáng)降水落區(qū)明顯落后于實(shí)況，重慶中部地區(qū)(紅色實(shí)線圈所示)的落區(qū)預(yù)報(bào)較準(zhǔn)確。Xgboost仍然抓住了主雨帶動(dòng)態(tài)，較好預(yù)報(bào)了黔江、彭水地區(qū)的強(qiáng)降水(如圖中曲線所示)。對于切變線上的雨帶東移和湖南西北部新生的強(qiáng)降水區(qū)(紅色虛線圈)，Xgboost都有所體現(xiàn)，即彭水—黔江和務(wù)川—酉陽—咸豐(紅色曲線所示)分別有概率高值區(qū)與未來3 h內(nèi)相應(yīng)地區(qū)的實(shí)況短時(shí)強(qiáng)降水(黑色曲線所示)相對應(yīng)。

圖8 EC細(xì)網(wǎng)格再分析場在7月21日17:00—20:00降水量(a)、Xgboost模型21日14:00起報(bào)的短時(shí)強(qiáng)降水客觀概率預(yù)報(bào)(b)和21日17:00—20:00實(shí)況降水量(c)以及17:00—18:00(d)、18:00—19:00(e)、19:00—20:00(f)小時(shí)降水量

圖9 EC細(xì)網(wǎng)格再分析場在7月22日02:00—05:00降水量(a)、Xgboost模型22日02:00起報(bào)的短時(shí)強(qiáng)降水客觀概率預(yù)報(bào)(b)和22日02:00—05:00實(shí)況降水量(c)以及02:00—03:00(d)、03:00—04:00(e)、04:00—05:00(f)小時(shí)降水量

如圖10所示，在此次過程的兩個(gè)時(shí)段中，Xgboost預(yù)報(bào)的TS、POD、FAR分別在0.2～0.4、0.6～1和0.6～0.8之間。TS和FAR評分與前一過程總體持平，POD略高于前一過程，但兩次過程的Xgboost預(yù)報(bào)無論從定量還是定性結(jié)果來看都優(yōu)于EC。從以上個(gè)例分析可見，該方法可以較好預(yù)報(bào)短時(shí)強(qiáng)降水發(fā)生的概率及落區(qū)，對短臨預(yù)報(bào)預(yù)警具有一定參考價(jià)值。

圖10 兩次短時(shí)強(qiáng)降水個(gè)例的Xgboost預(yù)報(bào)的各項(xiàng)評分隨時(shí)次變化(起報(bào)時(shí)刻分別為6月29日02:00、14:00以及7月21日14:00和22日02:00)

2.3 近年短時(shí)強(qiáng)降水過程回報(bào)檢驗(yàn)

收集了2016、2017和2019年的主要短時(shí)強(qiáng)降水過程，并使用本文創(chuàng)建的預(yù)報(bào)模型進(jìn)行回報(bào)檢驗(yàn)。同時(shí)，檢驗(yàn)了EC再分析場在相應(yīng)時(shí)次的預(yù)報(bào)效果，結(jié)果見表5。

表5 Xgboost模型預(yù)報(bào)和EC預(yù)報(bào)在2016—2019年幾次短時(shí)強(qiáng)降水過程中的檢驗(yàn)

由表5檢驗(yàn)結(jié)果可見，Xgboost預(yù)報(bào)的TS評分在0.1以上，POD在0.2左右，相對2015年獨(dú)立樣本測試的結(jié)果較低，但在各個(gè)樣本集和指標(biāo)上均優(yōu)于EC細(xì)網(wǎng)格。已有研究表明[28]，常規(guī)業(yè)務(wù)中，短時(shí)強(qiáng)降水預(yù)報(bào)在第1小時(shí)的TS在0.1以下，并隨時(shí)效遞減,就此次回報(bào)檢驗(yàn)而言，其TS評分略高于文獻(xiàn)指出的常規(guī)業(yè)務(wù)水平。FAR總體較高，達(dá)到0.7左右，模型的空報(bào)較多。值得說明的是，本文采用了嚴(yán)格的時(shí)空點(diǎn)對點(diǎn)二分類檢驗(yàn)，即預(yù)報(bào)了短時(shí)強(qiáng)降水的小時(shí)時(shí)段以及格點(diǎn)與實(shí)況完全一致時(shí)才判斷為命中，且EC細(xì)網(wǎng)格再分析場具有全球模式對中小尺度對流性強(qiáng)降水預(yù)報(bào)偏弱的固有局限性，另外在插值模式降水統(tǒng)一分辨率的過程中也會(huì)削弱、平滑強(qiáng)降水，以上因素可能導(dǎo)致了檢驗(yàn)結(jié)果中EC表現(xiàn)差。綜上，Xgboost模型相對EC細(xì)網(wǎng)格在短時(shí)強(qiáng)降水預(yù)報(bào)上具有明顯優(yōu)勢，業(yè)務(wù)應(yīng)用中也具有一定參考意義。

3 結(jié)論和討論

盡管高分辨率數(shù)值模式不斷發(fā)展，但其對短時(shí)強(qiáng)降水等強(qiáng)對流天氣的預(yù)報(bào)能力仍然嚴(yán)重不足。隨著大數(shù)據(jù)挖掘和機(jī)器學(xué)習(xí)在各個(gè)領(lǐng)域的大放異彩，基于高時(shí)空分辨率模式、結(jié)合了機(jī)器學(xué)習(xí)算法的短時(shí)強(qiáng)降水客觀預(yù)報(bào)技術(shù)成為一種非常值得嘗試的思路。本文選取了重慶地區(qū)2011—2015年 5—9月間逐小時(shí)實(shí)況格點(diǎn)降水場，在利用EC細(xì)網(wǎng)格模式的再分析資料計(jì)算預(yù)報(bào)變量的基礎(chǔ)上，通過觀察變量的箱線圖差異指數(shù)確定了閾值進(jìn)行消空，初步剔除了噪音樣本并建立訓(xùn)練集，然后結(jié)合K均值聚類和改進(jìn)的Relief算法，構(gòu)造了類別平衡的訓(xùn)練集并挑選出預(yù)測能力較強(qiáng)的預(yù)報(bào)變量最終進(jìn)入模型，最后依托Xgboost算法建立起短時(shí)強(qiáng)降水預(yù)報(bào)模型，可以輸出概率預(yù)報(bào)或用戶自定義概率閾值生成的確定性預(yù)報(bào)，可對目前業(yè)務(wù)中的雷達(dá)降水估測和模式預(yù)報(bào)形成補(bǔ)充。2015年獨(dú)立樣本測試和近年來短時(shí)強(qiáng)降水過程的回報(bào)檢驗(yàn)表明，該方法對重慶地區(qū)的短時(shí)強(qiáng)降水預(yù)報(bào)較EC細(xì)網(wǎng)格更好，其產(chǎn)品在實(shí)際業(yè)務(wù)中也具有一定參考意義。本文主要結(jié)論歸納如下：

(1)基于EC細(xì)網(wǎng)格再分析資料計(jì)算了短時(shí)強(qiáng)降水預(yù)報(bào)變量，并根據(jù)各預(yù)報(bào)變量的箱線圖差異指數(shù)IBD制定了閾值法消空規(guī)則，通過剔除短時(shí)強(qiáng)降水潛勢過低的時(shí)次來提高短時(shí)強(qiáng)降水樣本在總體樣本中的占比，達(dá)到消空目的并做出初步預(yù)報(bào)。其中850 hPa渦度、K指數(shù)、修正K指數(shù)、700 hPa散度、700 hPa水汽通量散度以及地形高度的IBD絕對值大于0.2，相對其他變量較高。

(2)結(jié)合K均值聚類算法和改進(jìn)的Relief算法，建立了類別平衡的新訓(xùn)練集，并考察了變量對短時(shí)強(qiáng)降水的預(yù)測能力。變量平均權(quán)重表明，抬升指數(shù)、整層總降水量、修正K指數(shù)、總指數(shù)、對流有效位能和風(fēng)暴強(qiáng)度指數(shù)等變量預(yù)測能力較突出，因此將其納入建模過程。

(3)在樣本初步消空、預(yù)報(bào)因子篩選和重建類別平衡訓(xùn)練集的基礎(chǔ)上，初始化了Xgboost算法，并通過網(wǎng)格搜索調(diào)參確立了最佳超參數(shù)，建立起Xgboost短時(shí)強(qiáng)降水客觀概率預(yù)報(bào)模型。

(4)2015年獨(dú)立樣本測試表明，當(dāng)概率閾值取0.1時(shí)，模型的AUC為0.92，總體分類效果較好，全體樣本的短時(shí)強(qiáng)降水TS評分可達(dá)0.3，好于EC再分析場。對其中兩次個(gè)例分析表明，Xgboost短時(shí)強(qiáng)降水客觀概率預(yù)報(bào)能更好刻畫強(qiáng)降水發(fā)生的概率和落區(qū)，逐時(shí)次的預(yù)報(bào)效果仍優(yōu)于EC，TS評分在0.2～0.4之間。近年來短時(shí)強(qiáng)降水過程的回報(bào)檢驗(yàn)中模型預(yù)報(bào)的TS雖有降低，但仍高于EC再分析場，與業(yè)務(wù)水平持平，具有一定參考意義。

同時(shí)，本研究也存在以下幾點(diǎn)不足需要注意：受資料所限，本文采用6h間隔的再分析資料作為起報(bào)場計(jì)算預(yù)報(bào)變量，起報(bào)后的2～6 h內(nèi)預(yù)報(bào)效果會(huì)逐漸變差，但隨著該模型與業(yè)務(wù)EC細(xì)網(wǎng)格模式的對接，上述問題有望得到緩解；其次，氣候背景異常，如厄爾尼諾年的異常降水，以及模式在分辨率、參數(shù)化方案和同化方案等方面的變換更新都可能對模型穩(wěn)定性和預(yù)報(bào)結(jié)果產(chǎn)生不利影響，因此預(yù)報(bào)概率閾值需要根據(jù)氣候背景和模式升級重新確定；最后，本方法只能對短時(shí)強(qiáng)降水有無做出預(yù)報(bào)，無法指示其具體量級或強(qiáng)度。這些都是下一步工作中需要改進(jìn)和注意的。