摘" 要：利用南方地區(qū)15年地面觀測(cè)凍雨數(shù)據(jù)，結(jié)合有明顯凍雨形成物理機(jī)制的16種氣象因子，構(gòu)建具有可解釋性的XGBoost凍雨模型，其能夠很好地模擬非訓(xùn)練期間的凍雨個(gè)例和年凍雨日數(shù)的時(shí)間變化。利用1954—2023年的ERA5再分析資料驅(qū)動(dòng)該模型，獲得南方地區(qū)凍雨日數(shù)70年的長(zhǎng)時(shí)間序列。在長(zhǎng)達(dá)70年中，凍雨日數(shù)存在2～3年的顯著周期，且隨著全球溫度的升高，凍雨日數(shù)呈明顯的下降趨勢(shì)，其有3個(gè)突變點(diǎn)，分別為1958年、1967年和1990年，其中最長(zhǎng)的下降時(shí)段為1990—2023年。該研究為凍雨的氣候研究提供堅(jiān)實(shí)的數(shù)據(jù)支撐。

關(guān)鍵詞：凍雨日數(shù)；預(yù)報(bào)模型；可解釋性；機(jī)器學(xué)習(xí)；南方地區(qū)

中圖分類號(hào)：P426.6" " " 文獻(xiàn)標(biāo)志碼：A" " " " " 文章編號(hào)：2095-2945（2024）32-0087-04

Abstract： An interpretable XGBoost model of freezing rain was constructed by using surface freezing rain data observed in southern China over the past 15 years， and combining 16 meteorological factors with obvious physical mechanisms of freezing rain formation. The model was able to simulate well the temporal variation of individual freezing rain events and the annual number of freezing rain days during non-training periods. A 70-year long time series of freezing rain days in the southern region was obtained by XGBoost freezing rain model which driven using ERA5 reanalysis data from 1954 to 2023. The freezing rain days showed an obvious decreasing trend with the increase of global temperature， and existed a significant cycle of 2-3 years. There were three mutation points in the long 70-year period， which are 1958， 1967 and 1990， with the longest period of decline from 1990 to 2023. This study provides solid data support for climate change studies of freezing rain.

Keywords： freezing rain days; forecast model; interpretability; machine learning; southern region

凍雨是我國(guó)南方地區(qū)冬季發(fā)生的一種極端天氣現(xiàn)象，其量雖不及雨和雪，但嚴(yán)重的凍雨一旦發(fā)生，對(duì)自然環(huán)境系統(tǒng)和社會(huì)經(jīng)濟(jì)將造成巨大的危害[1-2]。如2008年冬季南方地區(qū)發(fā)生的凍雨事件，覆蓋范圍廣泛、持續(xù)時(shí)間長(zhǎng)，導(dǎo)致湖南、貴州、江西等地大面積斷電，交通中斷，農(nóng)作物顆粒無(wú)收，造成超千億元的經(jīng)濟(jì)損失[3]。

目前，凍雨預(yù)報(bào)主要依靠耦合了凍雨參數(shù)化方案的數(shù)值模式，如耦合在WRF模式中的凍雨參數(shù)化方案如Ramer方案、AFWA方案、Thompson方案和RUC方案等，它們根據(jù)凍雨發(fā)生時(shí)氣象要素的統(tǒng)計(jì)閾值和各降水相態(tài)形成的物理機(jī)制所得，但由于凍雨具有高度復(fù)雜性和多變性，目前的凍雨物理及經(jīng)驗(yàn)?zāi)Ｐ驮谔幚砩婕皟鲇甑闹T多氣象要素的復(fù)雜性和不確定性時(shí)仍面臨諸多挑戰(zhàn)，憑借經(jīng)驗(yàn)統(tǒng)計(jì)得到的氣象要素閾值范圍存在偏差，導(dǎo)致這些凍雨參數(shù)化方案對(duì)凍雨的預(yù)報(bào)評(píng)分相對(duì)較低。機(jī)器學(xué)習(xí)模型能夠更好地適應(yīng)非線性、多維度、非穩(wěn)態(tài)的氣象數(shù)據(jù)，通過(guò)對(duì)大量氣象數(shù)據(jù)的學(xué)習(xí)，提取氣象要素與凍雨之間的潛在非線性關(guān)系，這為凍雨預(yù)報(bào)提供了新的途徑。目前，已有學(xué)者利用機(jī)器學(xué)習(xí)模型構(gòu)建了我國(guó)凍雨預(yù)報(bào)模型，如王珊珊等[4]使用探空、地面觀測(cè)和高分辨率的模式預(yù)報(bào)產(chǎn)品作為特征參量，并使用決策樹的集成學(xué)習(xí)算法，構(gòu)建了長(zhǎng)江中下游地區(qū)的凍雨預(yù)報(bào)模型，比耦合了凍雨參數(shù)化方案的數(shù)值模式預(yù)報(bào)的準(zhǔn)確率有所提高，但這種機(jī)器學(xué)習(xí)模型仍缺乏凍雨物理機(jī)制，具有不可解釋性。另外，我國(guó)現(xiàn)有的凍雨觀測(cè)資料數(shù)據(jù)質(zhì)量差、時(shí)空分辨率低，凍雨序列不足25年（1995—2019年），僅有學(xué)者[5]根據(jù)風(fēng)速、相對(duì)濕度和氣溫等簡(jiǎn)單氣象要素構(gòu)建了1954—2008年凍雨的長(zhǎng)時(shí)間的序列，但因凍雨的物理機(jī)制非常復(fù)雜，其反演的凍雨序列具有很大的不確定性。加上地面測(cè)站分布稀疏且部分測(cè)站存在缺測(cè)、漏測(cè)的現(xiàn)象，導(dǎo)致凍雨數(shù)據(jù)質(zhì)量差、時(shí)空分辨率低，急需找到一種能夠較好模擬南方地區(qū)凍雨的預(yù)報(bào)模型，并利用該模型重新構(gòu)建長(zhǎng)期的、高時(shí)空分辨率的凍雨數(shù)據(jù)集，為我國(guó)南方地區(qū)凍雨長(zhǎng)期預(yù)報(bào)提供科學(xué)支撐。

針對(duì)上述存在的問(wèn)題，本文利用大量的凍雨個(gè)例和同期的ERA5再分析資料產(chǎn)品，構(gòu)建了基于凍雨物理機(jī)制的XGBoost凍雨預(yù)報(bào)模型，使其具有可解釋性。并通過(guò)TS、命中率、空?qǐng)?bào)率和漏報(bào)率評(píng)價(jià)該模型對(duì)凍雨預(yù)報(bào)的能力，并對(duì)1954—2023年逐小時(shí)南方凍雨日數(shù)進(jìn)行回溯，獲得高時(shí)空分辨率的凍雨數(shù)據(jù)集，構(gòu)建70年南方地區(qū)凍雨日數(shù)的時(shí)間序列，進(jìn)而分析其變化趨勢(shì)、突變性和周期性。

1" XGBoost凍雨模型的構(gòu)建

1.1" 數(shù)據(jù)的來(lái)源和預(yù)處理

資料來(lái)源：地面降水相態(tài)觀測(cè)資料來(lái)源于中國(guó)氣象局國(guó)家信息中心整編發(fā)布的Micaps資料，研究區(qū)域范圍為20～31°N，97～117°E，資料包含西藏、云南、四川、貴州、重慶、湖南、廣東、廣西、江西、福建、安徽和湖北省共729個(gè)氣象站點(diǎn)每日2時(shí)、8時(shí)、14時(shí)、20時(shí)（北京時(shí)）的小時(shí)變溫與變壓、風(fēng)向、風(fēng)速、過(guò)去天氣、現(xiàn)在天氣、露點(diǎn)溫度、總云量和本站氣壓等。探空觀測(cè)資料來(lái)源于美國(guó)懷俄明大學(xué)的天氣數(shù)據(jù)網(wǎng)站（http：//weather.uwyo.edu/wyoming/），研究區(qū)域所包含的騰沖、昆明、威寧和西昌等22個(gè)探空站點(diǎn)每日8時(shí)和20時(shí)的探空觀測(cè)結(jié)果。此外，對(duì)探空資料進(jìn)行插值處理，除每個(gè)探空資料處特性層外，間隔50 hPa進(jìn)行插值，并采用Poore等[6]的溫度露點(diǎn)差判斷方法，計(jì)算云頂高度及對(duì)應(yīng)的云頂溫度。ERA5再分析資料來(lái)源于歐洲中期氣象預(yù)報(bào)中心（ECMWF）基于數(shù)值模式和資料同化技術(shù)構(gòu)建的高時(shí)空分辨率氣象再分析數(shù)據(jù)集，該資料逐小時(shí)進(jìn)行更新，水平網(wǎng)格分辨率為0.25°×0.25°。

數(shù)據(jù)預(yù)處理：提取南方地區(qū)2000—2015年地面觀測(cè)降水相態(tài)數(shù)據(jù)，整理得到凍雨數(shù)據(jù)4 954條，按雨、雪、凍雨和雨夾雪的順序標(biāo)號(hào)為1—4。ERA5再分析產(chǎn)品中的氣象要素采用離差標(biāo)準(zhǔn)化（即min-max標(biāo)準(zhǔn)化），使原始數(shù)據(jù)在0～1區(qū)間波動(dòng)；同時(shí)，根據(jù)對(duì)凍雨物理機(jī)制的研究，暖層、融化層和再凍結(jié)層是否存在等這些因子對(duì)凍雨的形成非常關(guān)鍵，因此，對(duì)這些定性的氣象數(shù)據(jù)如是否存在再凍結(jié)層、融化層和暖層等采用“是”為1，“否”為0進(jìn)行處理。

為了解決因子之間多重共線性的問(wèn)題，采用方差膨脹因子（VIF）來(lái)衡量獨(dú)立變量中的多重共線性。對(duì)47種氣象要素進(jìn)行多重共線性檢查，當(dāng)VIF顯著高于其他變量時(shí)，直接丟棄該變量，當(dāng)VIFlt;5，則認(rèn)為不存在顯著的共線性關(guān)系[7]，最終選定16種氣象因子作為輸入量，分別為降水量、氣溫、露點(diǎn)溫度、濕球溫度、云頂溫度、能見度、液水含量、海拔高度、相對(duì)濕度、最小相對(duì)濕度、云頂高度、次凍結(jié)層頂部高度、暖層高度、風(fēng)向、含冰率和探空廓線類型。

這些變量對(duì)凍雨的形成均有明顯的物理機(jī)制，如相對(duì)濕度能準(zhǔn)確劃分云頂和云底的邊界[8]，氣溫、濕球溫度、露點(diǎn)溫度、海拔高度、含冰率、暖層高度、降水量、液水含量和風(fēng)向等要素均是凍雨參數(shù)化方案中的關(guān)鍵要素[9-13]，將其作為輸入量構(gòu)建的XGBoost凍雨模型具有了很好的解釋性。

1.2" 預(yù)報(bào)效果評(píng)估指標(biāo)

預(yù)報(bào)效果評(píng)估采用小概率成功指標(biāo)即風(fēng)險(xiǎn)評(píng)分（Threat Score），簡(jiǎn)稱為TS評(píng)分。因?yàn)閮鲇晔且环N小概率事件，采用TS評(píng)分必須報(bào)“有”且實(shí)況為“有”時(shí)才有記分，其有效避免了報(bào)“無(wú)”而導(dǎo)致預(yù)報(bào)準(zhǔn)確率很高的弊病。其算法如下

小概率成功指標(biāo)TS

命中率POD

空?qǐng)?bào)率FAR

漏報(bào)率MAR

式中：NA為實(shí)際發(fā)生且預(yù)報(bào)發(fā)生次數(shù)；NB為實(shí)際不發(fā)生但預(yù)報(bào)發(fā)生次數(shù)；NC為實(shí)際發(fā)生但預(yù)報(bào)不發(fā)生次數(shù)。

1.3" XGBoost凍雨模型的構(gòu)建

在XGBoost對(duì)凍雨模型的訓(xùn)練中，以7∶3的比例將凍雨數(shù)據(jù)隨機(jī)分配成訓(xùn)練集和測(cè)試集，即3 467條凍雨數(shù)據(jù)作為訓(xùn)練集，1 487條數(shù)據(jù)數(shù)據(jù)作為測(cè)試集。同時(shí)，匹配同時(shí)刻ERA5再分析資料產(chǎn)品中的16種氣象因子作為輸入量。

XGBoost凍雨方案中各參數(shù)設(shè)置如下：學(xué)習(xí)率learning_rate（0.05），控制每一步迭代中對(duì)模型參數(shù)的更新幅度；gamma（0.1）為減少損失最小值，alpha（0.01）為L(zhǎng)1正則化權(quán)重項(xiàng)，lamb-da（0.6）為L(zhǎng)2正則化權(quán)重項(xiàng)；max_depth（5）為每棵樹的最大深度，控制樹的復(fù)雜度；min_child_weight（2）為決定子樹停止的條件閾值，subsample（0.8）為每棵樹所用的子樣本的比例，colsample_bytree（0.8）為弱學(xué)習(xí)器所用特征的比例， num_class（5）為分類的類別數(shù)量，objective（multi：softmax）為優(yōu)化目標(biāo)此處設(shè)置為多分類問(wèn)題，n_estimators（2000）為弱學(xué)習(xí)器（即決策樹）的數(shù)量，num_boost_ro und（20）為迭代次數(shù)，即訓(xùn)練樹的數(shù)量。

2" XGBoost凍雨模型預(yù)報(bào)效果檢驗(yàn)

選擇2016—2019年南方22個(gè)探空測(cè)站的大氣溫濕廓線與同一地面站點(diǎn)及附近的凍雨觀測(cè)結(jié)果進(jìn)行匹配，共獲取凍雨發(fā)生時(shí)的探空廓線樣本156個(gè)，首先對(duì) XGBoost凍雨模型預(yù)報(bào)的凍雨個(gè)例與觀測(cè)進(jìn)行了TS評(píng)分（圖1），從圖1中可以看到TS評(píng)分為0.6，命中率0.92，漏報(bào)率0.05，顯示了非常好的凍雨個(gè)例預(yù)報(bào)能力。然后，使用XGBoost凍雨方案，對(duì)南方地區(qū)年凍雨日數(shù)進(jìn)行回溯。由于將2000—2015年的凍雨數(shù)據(jù)隨機(jī)取樣分為訓(xùn)練集和測(cè)試集，因此僅對(duì)1995—1999年和2016—2019年南方地區(qū)凍雨日數(shù)時(shí)間序列進(jìn)行回溯。模擬結(jié)果（圖2）顯示，南方地區(qū)凍雨日數(shù)與觀測(cè)凍雨日數(shù)二者的相關(guān)系數(shù)均達(dá)到了0.99，通過(guò)了0.001的信度檢驗(yàn)，且均方根誤差平均在0.1左右，表明XGBoost凍雨方案模擬的年凍雨日數(shù)和與實(shí)際非常一致。

3" 南方地區(qū)70年凍雨日數(shù)時(shí)間變化特征

由于XGBoost凍雨模型能夠很好地模擬南方地區(qū)非訓(xùn)練期間的凍雨個(gè)例和年凍雨日數(shù)的時(shí)間變化，故而利用ERA5再分析資料驅(qū)動(dòng)XGBoost凍雨模型，構(gòu)建了1954—2023年逐小時(shí)中國(guó)南方地區(qū)3 km×3 km凍雨災(zāi)害數(shù)據(jù)集，進(jìn)而獲得了南方區(qū)域凍雨日數(shù)70年的長(zhǎng)時(shí)間序列。圖3（a）為1954—2023年南方地區(qū)凍雨日數(shù)的年際變化，70年平均凍雨發(fā)生天數(shù)為9.3 d，存在以2～3年為周期的顯著變化。70年間南方地區(qū)年凍雨日數(shù)呈下降趨勢(shì)，線性變化趨勢(shì)系數(shù)為-0.6/10 a，這種趨勢(shì)同全球氣候變暖密切相關(guān)，全球平均溫度與70年南方地區(qū)區(qū)域凍雨的時(shí)間序列在1～2月呈現(xiàn)顯著的負(fù)相關(guān)，負(fù)相關(guān)系數(shù)達(dá)0.37，通過(guò)了99.99%的相關(guān)信度檢驗(yàn)。凍雨日數(shù)距平值的年際變化（圖3（b））顯示，1990年以前凍雨發(fā)生天數(shù)以正距平為主，其中1977年正距平值最大為19.7 d，次大值出現(xiàn)在1984年為16.7 d；凍雨日數(shù)的負(fù)距平主要出現(xiàn)的時(shí)段在1958—1966年和1990—2023年，其中1991、2004年和2020年無(wú)凍雨出現(xiàn)，其負(fù)距平值最大為-9.3 d，次大值分別出現(xiàn)在1997年、2002年和2017年。但在以負(fù)距平為主的1990—2023年，仍出現(xiàn)了較高的正距平如2008年的8.7 d和2005年的6.7 d。為了分析其是否發(fā)生突變，使用累計(jì)距平法，發(fā)現(xiàn)凍雨日數(shù)在1958年、1967年、1990年存在突變點(diǎn)，并得到了Mann-Kendall檢驗(yàn)的證實(shí)，最終確定1958年、1967年、1990年為真實(shí)突變點(diǎn)，其中1990年為顯著突變點(diǎn)。最終確定有2個(gè)顯著凍雨日數(shù)下降時(shí)段為1958—1967年、1990—2023年；3個(gè)顯著凍雨日數(shù)上升時(shí)段為1954—1958、1967—1973年、1980—1990年。

4" 結(jié)論

本研究提取了南方地區(qū)2000—2015年地面觀測(cè)凍雨數(shù)據(jù)4 954條，通過(guò)數(shù)據(jù)的預(yù)處理，獲得了對(duì)凍雨的形成有明顯的物理機(jī)制的16種氣象因子，然后以7∶3的比例將凍雨數(shù)據(jù)隨機(jī)分配成訓(xùn)練集和測(cè)試集，匹配同時(shí)刻16種氣象因子作為輸入量構(gòu)建了XGBoost凍雨模型。該模型不但能夠很好地模擬南方地區(qū)非訓(xùn)練期間的凍雨個(gè)例和年凍雨日數(shù)的時(shí)間變化，而且構(gòu)建的XGBoost凍雨模型具有了很好的解釋性。利用ERA5再分析資料驅(qū)動(dòng)該模型，獲得了南方地區(qū)凍雨日數(shù)70年的長(zhǎng)時(shí)間序列。在70年間，凍雨日數(shù)存在2～3年的顯著周期，且隨著全球溫度的升高，凍雨日數(shù)呈明顯的下降趨勢(shì)，其中1958年、1967年和1990年為突變點(diǎn)，最長(zhǎng)的下降時(shí)段為1990—2023年。

參考文獻(xiàn)：

[1] 陶詩(shī)言，衛(wèi)捷.2008年1月我國(guó)南方嚴(yán)重冰雪災(zāi)害過(guò)程分析[J].氣候與環(huán)境研究，2008（4）：337-350.

[2] 黃強(qiáng)，王家紅，歐名勇.2005年湖南電網(wǎng)冰災(zāi)事故分析及其應(yīng)對(duì)措施[J].電網(wǎng)技術(shù)，2005（24）：16-19.

[3] 曾明劍，陸維松，梁信忠，等.2008年初中國(guó)南方持續(xù)性冰凍雨雪災(zāi)害形成的溫度場(chǎng)結(jié)構(gòu)分析[J].氣象學(xué)報(bào)，2008，66（6）：1043-1052.

[4] 王珊珊，雷彥森，方鴻斌，等.基于CatBoost算法的長(zhǎng)江中游冬季降水相態(tài)預(yù)報(bào)方法研究[J].氣象，2022，48（9）：1153-1161.

[5] 王遵婭.中國(guó)冰凍日數(shù)的氣候及變化特征分析[J].大氣科學(xué)，2011，35（3）：411-421.

[6] POORE K D， WANG J， ROSSOW W B. Cloud layer thicknesses from a combination of surface and upper-air observations [J]. Journal of Climate， 1995，8（3）：550-568.

[7] BELSLEY D A， KUH E， WELSCH R E. Regression diagnostics： Identifying influential data and sources of collinearity[M]. John Wiley amp; Sons， 2005.

[8] 周毓荃，歐建軍.利用探空數(shù)據(jù)分析云垂直結(jié)構(gòu)的方法及其應(yīng)用研究[J].氣象，2010，36（11）：50-58.

[9] RAUBER R M， OLTHOFF L S， RAMAMURTHY M K，et al. The relative importance of warm rain and melting processes in freezing precipitation events[J]. Journal of Applied Meteorology， 2000，39（7）：1185-1195.

[10] RAMER J. An empirical technique for diagnosing precipitation type from model output[C]// International Conference on Aviation Weather Systems， 5 th， Vienna， VA. 1993： 227-230.

[11] 漆梁波.我國(guó)冬季凍雨和冰粒天氣的形成機(jī)制及預(yù)報(bào)著眼點(diǎn)[J].氣象，2012，38（7）：769-778.

[12] 趙培娟，邵宇翔，李周，等.凍雨形成的天氣條件分析[J].氣象與環(huán)境科學(xué)，2008，31（4）：36-39.

[13] 牛生杰，周悅，賈然，等.電線積冰微物理機(jī)制初步研究：觀測(cè)和模擬[J].中國(guó)科學(xué)：地球科學(xué)，2011，41（12）：1812-1821.