基于Xgboost算法的短時強降水預(yù)報方法

朱巖 翟丹華 吳志鵬 張焱

(重慶市氣象臺,重慶 401147)

引言

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

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

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

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

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

1 資料和方法

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

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

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

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

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

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

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

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

1.3 方法

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

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

(1)

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

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

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

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

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

(2)

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

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

1.3.2 類別平衡和特征選擇

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

(3)

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

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

(4)

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

(5)

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

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

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

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

(6)

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

(7)

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

(8)

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

各個變量的平均權(quán)重如圖3所示，平均權(quán)重大于閾值的變量在圖中以黑色柱體顯示，有風(fē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ù)報變量平均權(quán)重(黑色為大于閾值的變量)

1.3.3 評價指標(biāo)

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

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

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

TS=NA/(NA+NB+NC)

(9)

POD=NA/(NA+NB)

(10)

FAR=NC/(NA+NC)

(11)

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

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

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

1.4 Xgboost建模方法

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

(12)

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

(13)

(14)

Xgboost有包括正則化項、學(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)險，使模型優(yōu)化度和模型復(fù)雜度達到平衡，盡可能同時避免欠擬合和過擬合。網(wǎng)格搜索是應(yīng)用最廣泛的建立在交叉驗證基礎(chǔ)上的超參數(shù)搜索算法，這種窮舉式調(diào)參算法通過循環(huán)遍歷嘗試每種參數(shù)組合的可能性，找出表現(xiàn)最好的組合，能夠找到全局最大或最小值。過程中采用了5折交叉驗證，也即將訓(xùn)練集5等分，取其中一份為驗證集，其余4份為新訓(xùn)練集，經(jīng)過5次在不同驗證集上的測試，取最優(yōu)結(jié)果所對應(yīng)的超參數(shù)組合。最佳關(guān)鍵超參數(shù)組合如表4所示：

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

2 試驗結(jié)果

2.1 總體檢驗結(jié)果

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

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

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

2.2 獨立樣本檢驗期間兩次短時強降水預(yù)報分析

2.2.1 渝西南短時強降水過程分析

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

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

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

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

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

2.2.2 渝西短時強降水過程分析

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

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

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

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

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

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

2.3 近年短時強降水過程回報檢驗

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

表5 Xgboost模型預(yù)報和EC預(yù)報在2016—2019年幾次短時強降水過程中的檢驗

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

3 結(jié)論和討論

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

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

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

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

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

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