基于逐步多元線性回歸和隨機(jī)森林模型預(yù)測(cè)黃河流域極端氣溫事件

陳俊清,李 毅,王 斌,楊雪寧,劉峰貴

(1. 西北農(nóng)林科技大學(xué) 水利與建筑工程學(xué)院,教育部旱區(qū)農(nóng)業(yè)水土工程重點(diǎn)實(shí)驗(yàn)室,陜西 楊凌 712100;2. 水利部黃河流域水治理與水安全重點(diǎn)實(shí)驗(yàn)室(籌),河南 鄭州 450003; 3. 澳大利亞新南威爾士州初級(jí)產(chǎn)業(yè)部,澳大利亞 新南威爾士州 NSW 2650; 4. 青海師范大學(xué) 地理科學(xué)學(xué)院,青海 西寧 810016)

0 引言

聯(lián)合國(guó)政府間氣候變化專門委員會(huì)(Intergovernmental Panel on Climate Change, IPCC)第六次報(bào)告指出,2001—2020年全球平均地表溫度較1850—1900年升高約0.99 ℃[1]。全球變暖背景下,干旱、暴雨和熱浪等極端事件呈增加態(tài)勢(shì),不僅加劇了旱澇災(zāi)害、造成農(nóng)作物減產(chǎn),而且嚴(yán)重危害了人體健康[2]。極端氣候事件是指天氣變量觀測(cè)值超過(guò)某一閾值時(shí)的小概率事件[3]。氣候變化檢測(cè)和指數(shù)專家組(Expert Team on Climate Change Detection and Indices, ETCCDI)定義了27個(gè)具有代表性的極端氣候指數(shù)[4],包括極端氣溫和極端降水兩大類指數(shù),已被廣泛應(yīng)用于全球或區(qū)域尺度極端氣候研究[5]。其中ETCCDI定義了16個(gè)極端氣溫指數(shù)(extreme temperature index, ETI),可從不同角度評(píng)估極端氣溫事件,例如強(qiáng)度指數(shù)有最高氣溫(TXx)和最低氣溫(TNn),相對(duì)閾值指數(shù)有暖晝天數(shù)(TX90p)和冷夜天數(shù)(TN10p),以及絕對(duì)閾值指數(shù)霜凍天數(shù)(FD0)等[6-7],而酷熱天數(shù)(TD30)作為表征極端氣溫事件的絕對(duì)閾值暖指數(shù),在極端氣溫事件評(píng)估中具有很好的適用性[8]。

大量研究人員對(duì)極端氣溫指數(shù)進(jìn)行分析預(yù)測(cè),為緩解氣候變化可能造成的不利影響提供了參考。研究者們采用多種方法對(duì)ETI與其關(guān)鍵影響因子進(jìn)行建模,并預(yù)報(bào)極端氣溫事件。如逐步多元線性回歸模型(stepwise multiple linear regression, SMLR)可較好地描述極端氣溫指數(shù)與多個(gè)影響因子之間的定量關(guān)系,如CHEN等[9]利用SMLR模型量化分析了湖北省植被對(duì)極端氣候指數(shù)的敏感性。隨著大數(shù)據(jù)時(shí)代的到來(lái),各種機(jī)器學(xué)習(xí)模型被廣泛應(yīng)用,其中隨機(jī)森林(random forest, RF)模型是BREIMAN[10]建立的一種組合式機(jī)器學(xué)習(xí)方法,具有準(zhǔn)確性高、不易過(guò)擬合的優(yōu)點(diǎn),能有效處理大數(shù)據(jù)集并判定變量的重要性。LI等[11]利用機(jī)器學(xué)習(xí)方法對(duì)北京地區(qū)網(wǎng)格溫度進(jìn)行了預(yù)報(bào)。ZHOU等[12]對(duì)比了多元線性回歸(multiple linear regression, MLR)、梯度提升決策樹、k-最近鄰法、隨機(jī)森林、極端梯度提升(extreme gradient boosting, XGB)和深度神經(jīng)網(wǎng)絡(luò)(deep neural network, DNN)等6種機(jī)器學(xué)習(xí)方法對(duì)近地面氣溫的預(yù)報(bào)。

此外,極端氣溫事件與大氣環(huán)流變化的密切聯(lián)系為極端氣溫事件的預(yù)報(bào)提供了新的研究方向。大氣環(huán)流會(huì)促進(jìn)高低緯度之間及海陸之間的水熱交換,是影響氣候的重要因素。大氣環(huán)流異常往往會(huì)導(dǎo)致極端氣溫事件的發(fā)生,因此,常用環(huán)流指數(shù)來(lái)定量表征大氣環(huán)流的變化,而大氣環(huán)流通常分為大氣類、海溫類和其他三類。目前,已有諸多學(xué)者開展了大量針對(duì)極端氣溫事件影響因素的研究,發(fā)現(xiàn)大氣環(huán)流對(duì)極端氣候事件的發(fā)生具有重要影響[13-14]。梁蘇潔等[15]分析并采用了北極濤動(dòng)、北大西洋濤動(dòng)、太平洋—北美型和平流層極渦4種大氣遙相關(guān)型指數(shù)構(gòu)建了多元回歸統(tǒng)計(jì)方程對(duì)華北地區(qū)冬季極端低溫事件進(jìn)行了預(yù)測(cè)。戴聲佩等[13]對(duì)華南地區(qū)極端氣溫事件進(jìn)行大尺度環(huán)流影響分析,指出厄爾尼諾—南方濤動(dòng)(El Nio-Southern Oscillation, ENSO)異常與華南地區(qū)極端氣溫事件關(guān)系密切,南海副高強(qiáng)度指數(shù)和西太平洋副高強(qiáng)度指數(shù)也對(duì)極端氣溫事件有明顯貢獻(xiàn)。

黃河發(fā)源于青藏高原,是我國(guó)第二大河,流經(jīng)內(nèi)蒙古高原、黃土高原和黃淮海平原,最終匯入太平洋。黃河流域位于32°N～42°N,96°E～119°E,面積約7.95×105km2,主要是干旱半干旱地區(qū),其年平均降水量為123～1021 mm,年平均氣溫為-4～14 ℃[16]。全球變暖背景下,黃河流域極端氣溫事件頻發(fā),苗書玲等[17]發(fā)現(xiàn)黃河流域極端氣溫指數(shù)在流域西北部有明顯的下降趨勢(shì),東部則上升趨勢(shì)明顯。張克新等[18]也發(fā)現(xiàn)黃河流域極端氣溫事件和厄爾尼諾—南方濤動(dòng)(ENSO)指數(shù)存在多尺度的顯著相關(guān)性。盡管目前國(guó)內(nèi)對(duì)極端氣溫事件的時(shí)空變化及其成因分析已有相關(guān)研究,但研究大多采用主成分分析[19]、基于SPSS軟件的因子分析[13]等統(tǒng)計(jì)分析方法,對(duì)于RF等機(jī)器學(xué)習(xí)模型的應(yīng)用還較少。此外,諸多學(xué)者在單個(gè)大氣環(huán)流因子如ENSO對(duì)極端氣溫事件的影響方面取得了成果[20],但極端氣溫事件的成因復(fù)雜,受到諸多因素的共同作用,綜合考慮多種大氣環(huán)流因子對(duì)黃河流域極端氣溫事件的預(yù)測(cè)研究還不夠深入。黃河流域是氣候變化的敏感區(qū),長(zhǎng)期以來(lái),頻發(fā)的極端氣溫事件對(duì)該區(qū)域生態(tài)系統(tǒng)和人類生活造成嚴(yán)重影響。探究黃河流域極端氣溫事件的氣候驅(qū)動(dòng)要素及其預(yù)測(cè),能夠?yàn)檫m應(yīng)氣候變化、流域生態(tài)環(huán)境保護(hù)和防震減災(zāi)提供重要的科學(xué)依據(jù)。

本文在收集1961年1月—2022年10月環(huán)流指數(shù)數(shù)據(jù)和1961年1月1日—2020年12月31日黃河流域80個(gè)站點(diǎn)氣溫?cái)?shù)據(jù)的基礎(chǔ)上,分析黃河流域極端氣溫事件的時(shí)空分布特征及其對(duì)環(huán)流的響應(yīng)規(guī)律,運(yùn)用SMLR和RF模型模擬極端氣溫指數(shù),并利用環(huán)流指數(shù)對(duì)未來(lái)極端氣溫指數(shù)進(jìn)行預(yù)測(cè),以期揭示黃河流域極端氣候事件的變化規(guī)律及成因,為極端氣溫的預(yù)報(bào)提供科學(xué)依據(jù)。

1 材料和方法

1.1 數(shù)據(jù)收集

黃河流域80個(gè)氣象站點(diǎn)1961年1月1日—2020年12月31日的逐日氣溫?cái)?shù)據(jù)來(lái)源于中國(guó)氣象科學(xué)數(shù)據(jù)中心(http://data.cma.cn),經(jīng)檢查該數(shù)據(jù)具有一致性和完整性,數(shù)據(jù)缺失率<1%,如圖1所示。

圖1 黃河流域氣象站點(diǎn)空間分布Fig.1 Spatial distribution of weather stations in the Yellow River Basin

研究選取1961年1月—2022年10月130個(gè)環(huán)流指數(shù)的逐月數(shù)據(jù),該數(shù)據(jù)來(lái)源于國(guó)家氣候中心(http://cmdp.ncc-cma.net/cn/download.htm),為保證數(shù)據(jù)的連續(xù)性,將連續(xù)缺失3個(gè)月及以上的環(huán)流指數(shù)剔除后保留了100個(gè)環(huán)流指數(shù)。

1.2 極端氣溫事件關(guān)鍵驅(qū)動(dòng)環(huán)流指數(shù)篩選

1.2.1 極端氣溫指數(shù)的選取和計(jì)算

選用描述頻率和強(qiáng)度的6個(gè)極端氣溫指數(shù),其中5個(gè)是ETCCDI定義的(http://etccdi.pacificclimate.org/list_27_indices.shtml),包括TXx、TX90p、TNn、TN10p和FD0。另一個(gè)選用的極端氣溫指數(shù)TD30具有生物學(xué)重要性且被廣泛應(yīng)用于極端氣溫事件研究[8]。TXx、TX90p和TD30是極端氣溫暖指數(shù),而TNn、TN10p和FD0屬于冷指數(shù)。6個(gè)ETI的詳細(xì)說(shuō)明如表1所示。

表1 極端氣溫指數(shù)說(shuō)明Table 1 Description of extreme temperature indexes

根據(jù)逐日氣象數(shù)據(jù)提取了各站點(diǎn)的逐月TXx、TX90p、TNn、TN10p、FD0和TD30。

1.2.2 環(huán)流指數(shù)的多重共線性分析

選留的100個(gè)環(huán)流指數(shù)之間可能有較高的相關(guān)性,因此在建模前需排除相互不獨(dú)立的環(huán)流指數(shù)。利用方差膨脹系數(shù)(variance inflation factor, VIF)對(duì)環(huán)流指數(shù)進(jìn)行多重共線性診斷[21]。VIF的計(jì)算式為

(1)

重復(fù)多次共線性診斷,逐步剔除VIF>10的具有顯著共線性[22]的環(huán)流指數(shù)。經(jīng)分析,最終保留了57個(gè)環(huán)流指數(shù)。

1.2.3 關(guān)鍵環(huán)流指數(shù)的篩選

單站點(diǎn)6個(gè)ETI對(duì)應(yīng)的關(guān)鍵環(huán)流指數(shù)具有差異,但相鄰站點(diǎn)關(guān)鍵環(huán)流指數(shù)較相似,故對(duì)各站點(diǎn)進(jìn)行關(guān)鍵環(huán)流指數(shù)的篩選,以便用于建模。用Pearson相關(guān)系數(shù)r表征每個(gè)ETI與各環(huán)流指數(shù)之間的相關(guān)程度[23]。初步分析兩者關(guān)系發(fā)現(xiàn)ETI和環(huán)流指數(shù)之間的相關(guān)性基本具有12個(gè)月的周期,因此,分別計(jì)算滯后0～12個(gè)月情況下、各ETI與經(jīng)過(guò)多重共線性分析篩選的57個(gè)環(huán)流指數(shù)之間的r值。r的取值范圍是-1～1,r>0(<0)表示正(負(fù))相關(guān),|r|越大相關(guān)性越高。根據(jù)r的臨界值表[24],當(dāng)顯著性水平P<0.01、樣本數(shù)為720時(shí),r的臨界值為0.097,即具有99%顯著性水平的關(guān)鍵環(huán)流指數(shù)的篩選條件為r>0.097。

1.3 ETI模型的建立及模擬

1.3.1 逐步多元線性回歸模型(SMLR)

由于ETI與各關(guān)鍵環(huán)流指數(shù)的相關(guān)系數(shù)r隨滯后時(shí)間(月份)變化,因此,取r最大(rmax)時(shí)的滯后時(shí)間對(duì)應(yīng)的多個(gè)關(guān)鍵環(huán)流指數(shù)進(jìn)行建模。TX90p、TN10p、FD0和TD30表征極端氣溫事件發(fā)生的頻率,單位為d、且只能為正值。SMLR模型采用原始數(shù)據(jù)得到的模擬值存在負(fù)值,因此將各站點(diǎn)TX90p、TN10p、FD0和TD30的序列進(jìn)行對(duì)數(shù)轉(zhuǎn)換。以1961—2010年的數(shù)據(jù)作為率定集,2011—2020年的數(shù)據(jù)為驗(yàn)證集,建立各對(duì)數(shù)轉(zhuǎn)換后的ETI序列和特定滯后時(shí)間下關(guān)鍵環(huán)流指數(shù)之間的SMLR模型[25],該模型表示為

y=a1x1,j+a2x2,j+aixi,j+…+anxn,j+b

(2)

式中:y為ETI,其中TXx和TNn采用原始數(shù)據(jù),而TX90p、TN10p、FD0和TD30為對(duì)數(shù)轉(zhuǎn)換后的序列;xi,j為對(duì)應(yīng)rmax的滯后時(shí)間為j個(gè)月的第i個(gè)環(huán)流指數(shù);ai為偏回歸系數(shù);b為截距。

1.3.2 隨機(jī)森林RF模型

隨機(jī)森林模型是通過(guò)Bootstrap重抽樣構(gòu)建決策樹的組合分類器模型[10],具有良好的分類篩選與預(yù)測(cè)性能。利用與SMLR模型數(shù)量及滯后時(shí)間一致的環(huán)流指數(shù),對(duì)各ETI采用RF建模,以1961—2010年的數(shù)據(jù)作為率定集,2011—2020年的數(shù)據(jù)為驗(yàn)證集,并與SMLR模型的建模效果進(jìn)行比較。

在RF模型訓(xùn)練中需確定2個(gè)參數(shù):決策樹數(shù)量(ntree)和每個(gè)樹節(jié)點(diǎn)上隨機(jī)變量的數(shù)量(mtry)。對(duì)于回歸問(wèn)題,ntree一般取500;mtry設(shè)置為回歸模型所有預(yù)測(cè)變量數(shù)的1/3[26]。RF模型通過(guò)均方誤差增量(increase in mean squared error, IncMSE)評(píng)估變量的重要性,對(duì)每棵決策樹的輸入變量xi進(jìn)行隨機(jī)擾動(dòng)后,檢驗(yàn)袋外數(shù)據(jù)(out of band, OOB)的均方誤差,IncMSE值越大表示變量越重要[27],計(jì)算公式為

(3)

式中:IncMSEi為第i個(gè)變量的IncMSE(%);MSEi1,OOB為未受到擾動(dòng)的袋外數(shù)據(jù)的均方誤差;MSEi2,OOB為隨機(jī)擾動(dòng)后袋外數(shù)據(jù)的均方誤差。

由于模型實(shí)際包含的環(huán)流指數(shù)數(shù)量及其重要性在黃河流域80個(gè)站點(diǎn)存在差異,因此,提出新的重要性評(píng)估指數(shù)RankYR(%),以表征各環(huán)流指數(shù)在整個(gè)黃河流域上對(duì)極端氣溫指數(shù)的影響程度,表達(dá)式為

(4)

式中:Ranki,YR為第i個(gè)環(huán)流指數(shù)的黃河流域總體重要性(%);SIii為站點(diǎn)ii模型實(shí)際包含的環(huán)流指數(shù)的數(shù)量;Ranki,ii為第i個(gè)環(huán)流指數(shù)在站點(diǎn)ii中的變量重要性名次;Impi,ii為第i個(gè)環(huán)流指數(shù)在站點(diǎn)ii的變量重要性IncMSE(%);SImax為黃河流域?qū)嶋H包含環(huán)流指數(shù)最多的站點(diǎn)的環(huán)流指數(shù)數(shù)量。

RankYR值越大,說(shuō)明該環(huán)流指數(shù)對(duì)整個(gè)黃河流域ETI的影響越大,反之亦然。

1.4 模型效果評(píng)價(jià)

用決定系數(shù)(R-square,R2)、林氏一致性相關(guān)系數(shù)[28](Lin’s concordance correlation coefficient, LCCC)、平均絕對(duì)誤差(mean absolute error, MAE)以及均方根誤差(root mean square error, RMSE)來(lái)評(píng)估SMLR和RF模型率定和驗(yàn)證的效果。評(píng)價(jià)指標(biāo)的計(jì)算式為

(5)

(6)

(7)

(8)

2 結(jié)果與分析

2.1 極端氣溫指數(shù)ETI的時(shí)空變化

黃河流域1961—2020年6個(gè)ETI的年際變化曲線及年極端氣溫事件的多年平均情況如圖2所示。從時(shí)間變化來(lái)看,極端氣溫暖指數(shù)TXx,TX90p和TD30呈現(xiàn)出波動(dòng)上升的趨勢(shì),說(shuō)明黃河流域極端高溫事件的強(qiáng)度和頻率均呈增加趨勢(shì),對(duì)全球變暖呈正響應(yīng)。極端氣溫冷指數(shù)FD0、TN10p呈下降趨勢(shì),TNn則呈上升趨勢(shì),表征黃河流域極端低溫事件的發(fā)生頻率在減小,強(qiáng)度也有所降低。

在空間分布上,由圖2可知,TXx、TX90p和TD30的值在黃河流域下游低海拔地區(qū)均較高;同時(shí)該區(qū)域的TNn較大,TN10p和FD0的值則較低,表明該地區(qū)極端高溫事件的強(qiáng)度和頻率均較大,而極端低溫事件的強(qiáng)度較弱且頻率較低,這可能是受到緯度和地形的影響。總的來(lái)說(shuō),黃河流域極端高溫事件強(qiáng)度和頻率的空間分布情況與極端低溫事件大致相反。

2.2 關(guān)鍵環(huán)流指數(shù)及其滯后時(shí)間

分析在滯后0～48個(gè)月的情況下,ETI與經(jīng)過(guò)多重共線性分析的57個(gè)環(huán)流指數(shù)的Pearson相關(guān)系數(shù)r的變化,發(fā)現(xiàn)其具有明顯的周期性,周期為12個(gè)月。因此對(duì)滯后時(shí)間為0～12個(gè)月的情況進(jìn)行具體分析。57個(gè)環(huán)流指數(shù)中大氣環(huán)流指數(shù)共35個(gè),分別是北非副高面積指數(shù)(North African Subtropical High Area Index, NAHAI)、西太平洋副高強(qiáng)度指數(shù)(Western Pacific Subtropical High Intensity Index, WPSH)、東太平洋副高強(qiáng)度指數(shù)(Eastern Pacific Subtropical High Intensity Index, EPSH)、太平洋區(qū)極渦強(qiáng)度指數(shù)(Pacific Polar Vortex Intensity Index, PPVI)、北大西洋—?dú)W洲區(qū)極渦強(qiáng)度指數(shù)(Atlantic-European Polar Vortex Intensity Index, AEPVI)、北半球極渦中心經(jīng)向位置指數(shù)(Northern Hemisphere Polar Vortex Central Longitude Index, NVCL)、北半球極渦中心緯向位置指數(shù)(Northern Hemisphere Polar Vortex Central Latitude Index, NVCLI)、歐亞經(jīng)向環(huán)流指數(shù)(Eurasian Meridional Circulation Index, EMC)、亞洲緯向環(huán)流指數(shù)(Asian Zonal Circulation Index, AZC)、亞洲經(jīng)向環(huán)流指數(shù)(Asian Meridional Circulation Index, AMC)、東亞槽位置指數(shù)(East Asian Trough Position Index, EATP)、印緬槽強(qiáng)度指數(shù)(India-Burma Trough Intensity Index, IBTI)、北極濤動(dòng)(Arctic Oscillation, AO)、南極濤動(dòng)指數(shù)(Antarctic Oscillation, AAO)、北大西洋濤動(dòng)(North Atlantic Oscillation, NAO)、北非—北大西洋—北美副高脊線位置指數(shù)(North African-North Atlantic-North American Subtropical High Ridge Position Index, NANRP)、東太平洋副高脊線位置指數(shù)(Eastern Pacific Subtropical High Ridge Position Index, EPRP)、南海副高脊線位置指數(shù)(South China Sea Subtropical High Ridge Position Index, SSRP)、西太平洋副高西伸脊點(diǎn)指數(shù)(Western Pacific Sub Tropical High Western Ridge Point Index, WHWRP)、亞洲區(qū)極渦面積指數(shù)(Asia Polar Vortex Area Index, APVA)、太平洋區(qū)極渦面積指數(shù)(Pacific Polar Vortex Area Index, PPVA)、北美區(qū)極渦面積指數(shù)(North American Polar Vortex Area Index, NAPVA)、大西洋歐洲區(qū)極渦面積指數(shù)(Atlantic-European Polar Vortex Area Index, AEPVA)、太平洋—北美遙相關(guān)型指數(shù)(Pacific—North American Pattern, PNA)、東大西洋遙相關(guān)型指數(shù)(East Atlantic Pattern, EA)、西太平洋遙相關(guān)型指數(shù)(West Pacific Pattern, WP)、東大西洋—西俄羅斯遙相關(guān)型指數(shù)(East Atlantic-West Russia Pattern, EAWR)、極地—?dú)W亞遙相關(guān)型指數(shù)(Polar-Eurasia Pattern, POL)、斯堪的納維亞遙相關(guān)型指數(shù)(Scandinavia Pattern, SCA)、50 hPa緯向風(fēng)指數(shù)(50 hPa zonal wind Index, X50ZW)、赤道中東太平洋200 hPa緯向風(fēng)指數(shù)(Mid-Eastern Pacific 200 hPa Zonal Wind Index, MPZW)、850 hPa西太平洋信風(fēng)指數(shù)(West Pacific 850 hPa Trade Wind Index, WPTW)、850 hPa中太平洋信風(fēng)指數(shù)(Central Pacific 850 hPa Trade Wind Index, CPTW)、850 hPa東太平洋信風(fēng)指數(shù)(East Pacific 850 hPa Trade Wind Index, EPTA)和北大西洋—?dú)W洲環(huán)流型C型指數(shù)(Atlantic-European Circulation C Pattern Index, ACCP)。15個(gè)海溫指數(shù)分別為NINO 1+2區(qū)海表溫度距平指數(shù)(NINO 1+2 SSTA Index, NINO1+2)、NINO W區(qū)海表溫度距平指數(shù)(NINO W SSTA Index, NINOW)、NINO A區(qū)海表溫度距平指數(shù)(NINO A SSTA Index, NINOA)、西太平洋暖池面積指數(shù)(Western Pacific Warm Pool Area Index, WPWPA)、大西洋多年代際振蕩指數(shù)(Atlantic Multi-decadal Oscillation Index, AMO)、親潮區(qū)海溫指數(shù)(Oyashio Current SST Index, OC)、西風(fēng)漂流區(qū)海溫指數(shù)(West Wind Drift Current SST Index, WWDC)、NINO B區(qū)海表溫度距平指數(shù)(NINO B SSTA Index, NINOB)、熱帶南大西洋海溫指數(shù)(Tropical Southern Atlantic SST Index, TSA)、西半球暖池指數(shù)(Western Hemisphere Warm Pool Index, WHWP)、印度洋暖池面積指數(shù)(Indian Ocean Warm Pool Area Index, IOWPA)、類ENSO指數(shù)(ENSO Modoki Index, EM)、東部型ENSO指數(shù)(Nino Eastern Pacific index, NE)、熱帶印度洋海溫偶極子指數(shù)(Tropic Indian Ocean Dipole Index, TIOD)和副熱帶南印度洋偶極子指數(shù)(South Indian Ocean Dipole Index, SIOD)。其他指數(shù)共7個(gè),分別為西太平洋編號(hào)臺(tái)風(fēng)數(shù)(Western North Pacific Typhoon Number, WNPTN)、登陸中國(guó)臺(tái)風(fēng)數(shù)(Number of Landing Typhoon on China, NLTC)、太陽(yáng)黑子指數(shù)(Total Sunspot Number Index, TSN)、南方濤動(dòng)指數(shù)(Southern Oscillation Index, SOI)、大西洋經(jīng)向模海溫指數(shù)(Atlantic Meridional Mode SST Index, AMM)、準(zhǔn)兩年振蕩指數(shù)(Quasi-Biennial Oscillation Index, QBO)、大西洋海溫三極子指數(shù)(North Atlantic Triple index, NAT)。

圖3以靖遠(yuǎn)站TXx為例,對(duì)比了在滯后0～12個(gè)月的情況下,TXx與經(jīng)過(guò)多重共線性分析的57個(gè)環(huán)流指數(shù)的Pearson相關(guān)系數(shù)r的變化,其中圖3(a)、(b)區(qū)分了對(duì)TXx具有較大(|rmax|>0.097)和較小(|rmax|<0.097)影響的環(huán)流指數(shù)組。由圖3可知:TXx受到多種環(huán)流指數(shù)的綜合影響,不同環(huán)流指數(shù)rmax對(duì)應(yīng)的滯后時(shí)間不同,大致都在0～12個(gè)月的范圍內(nèi)。說(shuō)明環(huán)流指數(shù)對(duì)TXx的影響存在滯后性,且該滯后效應(yīng)存在明顯的周期性。通過(guò)這種方式,在建模過(guò)程中可直接剔除圖3(b)中的15個(gè)環(huán)流指數(shù)。其他站點(diǎn)的r值各有不同,不再贅述。

進(jìn)一步對(duì)環(huán)流指數(shù)與TXx在不同滯后時(shí)間的Pearson相關(guān)關(guān)系進(jìn)行顯著性檢驗(yàn),并選擇rmax對(duì)應(yīng)的滯后時(shí)間,篩選出|r|>0.097且P<0.01的關(guān)鍵環(huán)流指數(shù),如圖4所示。關(guān)鍵環(huán)流指數(shù)與TXx在不同滯后時(shí)間的|r|的峰值介于0.098～0.893之間,平均為0.441,說(shuō)明大部分關(guān)鍵環(huán)流指數(shù)對(duì)TXx具有較為顯著的影響,rmax對(duì)應(yīng)的滯后時(shí)間分別為5、6、11、12個(gè)月。其他站點(diǎn)的環(huán)流指數(shù)滯后時(shí)間各有不同,不再贅述。

圖4 靖遠(yuǎn)站TXx與關(guān)鍵環(huán)流指數(shù)在不同滯后時(shí)間下的相關(guān)系數(shù)熱圖Fig.4 Heat map of correlations between TXx and key circulation indexes at Jingyuan station under different lag times

2.3 模型性能和變量重要性

2.3.1 ETI的對(duì)數(shù)轉(zhuǎn)換及描述性統(tǒng)計(jì)特征

因站點(diǎn)較多,此處以靖遠(yuǎn)站為例說(shuō)明4個(gè)以d為單位且只能為正值的ETI的原始值及其自然對(duì)數(shù)轉(zhuǎn)換值的描述性統(tǒng)計(jì)特征,如表2所示。TX90p、TD30、TN10p和FD0原始值的偏度系數(shù)分別為1.87、1.75、2.25和0.47,峰度系數(shù)分別為2.51、1.91、4.09、-1.60,呈現(xiàn)高度正偏態(tài)分布且在峰度上偏離正態(tài)分布。對(duì)原始值進(jìn)行對(duì)數(shù)轉(zhuǎn)換后,偏度系數(shù)分別為1.09、1.01、1.48、0.09,而峰度系數(shù)分別為-0.45、-0.65、0.51、-1.81,說(shuō)明對(duì)數(shù)轉(zhuǎn)換值總體上更接近正態(tài)分布。其他站點(diǎn)的ETI對(duì)數(shù)轉(zhuǎn)換序列統(tǒng)計(jì)特征有各自的特點(diǎn),此處不再贅述。

表2 靖遠(yuǎn)站TX90p、TD30、TN10p和FD0的原始值及其自然對(duì)數(shù)轉(zhuǎn)換值的描述性統(tǒng)計(jì)Table 2 Descriptive statistics of original values and natural logarithmic values for TX90p, TD30, TN10p and FD0 at Jingyuan station

2.3.2 建立的SMLR模型參數(shù)特征

基于篩選的關(guān)鍵環(huán)流指數(shù)進(jìn)行SMLR模型建模,并對(duì)黃河流域各站點(diǎn)的ETI進(jìn)行了模擬,得到具體方程,以靖遠(yuǎn)站為例展示了SMLR模型中對(duì)應(yīng)關(guān)鍵環(huán)流指數(shù)的偏回歸系數(shù)以及ETI與環(huán)流指數(shù)之間的相關(guān)性為rmax時(shí)對(duì)應(yīng)的滯后月數(shù),如表3所示。TX90p和TD30與各環(huán)流指數(shù)的rmax對(duì)應(yīng)的滯后時(shí)間完全相同。TXx和TNn及FD0和TN10p與環(huán)流指數(shù)的rmax對(duì)應(yīng)的滯后時(shí)間也較為接近,與TX90p和TD30的rmax對(duì)應(yīng)的滯后時(shí)間接近或間隔6個(gè)月。其他站SMLR模型的參數(shù)特征也已經(jīng)獲得,不再贅述。

表3 靖遠(yuǎn)站逐步多元線性回歸模型方程的偏回歸系數(shù)、截距及rmax對(duì)應(yīng)的滯后時(shí)間Table 3 Partial regression coefficient, intercept and the lag time for obtaining the maximum value of r in SMLR model equation at Jingyuan station

2.3.3 基于RF的環(huán)流指數(shù)重要性排序

以靖遠(yuǎn)站為例,基于RF模型獲得的環(huán)流指數(shù)的重要性IncMSE排序如圖5所示。環(huán)流指數(shù)對(duì)6個(gè)ETI的影響程度均不同,其中太平洋區(qū)極渦強(qiáng)度指數(shù)(PPVI)是影響TNn和FD0的最重要指數(shù),印度洋暖池面積指數(shù)(IOWPA)則是影響TX90p和TD30的最重要指數(shù),北大西洋—?dú)W洲區(qū)極渦強(qiáng)度指數(shù)(AEPVI)和北非—北大西洋—北美副高脊線位置指數(shù)(NANRP)分別是影響TXx和TN10p的最重要指數(shù)。太平洋區(qū)極渦強(qiáng)度指數(shù)(PPVI)、北非—北大西洋—北美副高脊線位置指數(shù)(NANRP)、南海副高脊線位置指數(shù)(SSRP)、東太平洋副高脊線位置指數(shù)(EPRP)和北非副高面積指數(shù)(NAHAI)對(duì)ETI的重要性排序最靠前。其他站的重要性排序結(jié)果也已獲得(數(shù)據(jù)略)。

圖5 靖遠(yuǎn)站影響極端氣溫指數(shù)的環(huán)流指數(shù)的重要性排序Fig.5 Rank of importance of circulation indexes influencing extreme temperature indexes at Jingyuan station

IncMSE僅體現(xiàn)環(huán)流指數(shù)在單站點(diǎn)的重要性,RankYR考慮了黃河流域全部80個(gè)站點(diǎn)上選取的關(guān)鍵環(huán)流指數(shù)對(duì)ETI的重要性。黃河流域RankYR>5%的環(huán)流指數(shù)如表4所示。PPVI是影響黃河流域TXx、TNn、TX90p和FD0的最重要的環(huán)流因子;NANRP對(duì)黃河流域TN10p和TD30的影響最大。

表4 影響黃河流域各極端氣溫指數(shù)的環(huán)流指數(shù)重要性RankYRTable 4 RankYR of circulation indexes affecting extreme temperature indexes in the Yellow River Basin %

2.3.4 模型性能評(píng)價(jià)

SMLR模型和RF模型模擬的黃河流域80個(gè)站點(diǎn)的ETI值的率定和驗(yàn)證結(jié)果如圖6所示。由圖可知:①在率定期,SMLR和RF模型R2的中位數(shù)變化范圍是74.7%～95.4%和74.4%～95.9%;LCCC中位數(shù)的變化范圍是0.83～0.98和0.78～0.98。SMLR模型MAE中位數(shù)的變化范圍是1.6～1.8 ℃和1.3～2.3 d,而RF的是1.6～1.8 ℃和1.3～2.2 d。SMLR模型RMSE的中位數(shù)變化范圍是2.1～2.3 ℃和2.7～4.0 d,RF模型的是2.0～2.3 ℃和2.8～3.5 d。②在驗(yàn)證期,SMLR和RF模型分別解釋了ETI變化的53.4%～94.6%和64.3%～95.4%;而LCCC中位數(shù)的變化范圍是0.71～0.97和0.69～0.97。SMLR模型的MAE中位數(shù)的變化范圍是1.7～1.8 ℃和1.5～2.9 d,而RF模型的是1.7～1.9 ℃和1.5～3 d。SMLR模型RMSE的中位數(shù)變化范圍是2.2～2.3 ℃和3.3～4.7 d,而RF模型的是2.1～2.3 ℃和2.8～4.4 d。2種模型對(duì)ETI的預(yù)測(cè)能力都比較好,其中RF模型模擬TNn的R2、LCCC、MAE和RMSE均優(yōu)于SMLR模型。RF模型模擬TX90p、TD30和TN10p時(shí)僅LCCC稍差,其他情況下RF模型均優(yōu)于SMLR模型。而RF模型模擬FD0的R2明顯優(yōu)于SMLR模型,LCCC、MAE、RMSE則略微差于SMLR模型。SMLR模型僅模擬6個(gè)ETI中TXx的效果優(yōu)于RF模型。③總的來(lái)說(shuō),RF模型的模擬效果優(yōu)于SMLR模型。驗(yàn)證期2個(gè)模型模擬ETI的R2,LCCC,MAE和RMSE的表現(xiàn)整體上略低于率定期,但未出現(xiàn)明顯的模擬性能下降,且RF模型模擬TN10p的MAE和RMSE在驗(yàn)證期的表現(xiàn)略優(yōu)于率定期。

圖6 逐步多元線性回歸和隨機(jī)森林模型在黃河流域80個(gè)站點(diǎn)模擬極端氣溫指數(shù)的率定和驗(yàn)證結(jié)果Fig.6 Calibrated and validated results of each extreme temperature index simulated by the SMLR and RF models at 80 sites in the Yellow River Basin

2.4 極端氣溫指數(shù)的預(yù)測(cè)

由于RF模型模擬ETI的性能總體優(yōu)于SMLR模型,且各環(huán)流指數(shù)對(duì)黃河流域極端氣溫指數(shù)的影響具有0～12個(gè)月的時(shí)滯特征,故基于RF模型利用2021年10月—2022年10月的環(huán)流指數(shù)數(shù)據(jù)預(yù)測(cè)了2022年11月黃河流域80個(gè)站點(diǎn)的ETI,其空間分布如圖7所示。由圖可知:2022年11月,極端高溫事件的強(qiáng)度(TXx)和發(fā)生頻率(TX90p和TD30)較高的站點(diǎn)主要集中在流域下游的低海拔地區(qū)和局部河谷地帶。極端低溫事件的頻率指數(shù)(TN10p、FD0)與極端高溫事件的頻率指數(shù)(TX90p、TD30)空間分布規(guī)律基本相反;TNn和TXx的空間分布特征基本相同,最低氣溫(TNn)越小,表征極端低溫事件的強(qiáng)度越大,即極端高溫事件與極端低溫事件的強(qiáng)度具有相反的空間分布規(guī)律。總體上,預(yù)測(cè)的2022年11月黃河流域各ETI的空間分布情況與多年平均情況相似。待新的環(huán)流指數(shù)更新后,可采用類似方法對(duì)黃河流域多站點(diǎn)的ETI進(jìn)行實(shí)時(shí)動(dòng)態(tài)預(yù)測(cè),從而對(duì)流域未來(lái)極端氣溫事件進(jìn)行預(yù)警。

圖7 基于隨機(jī)森林模型預(yù)測(cè)的2022年11月黃河流域極端氣溫指數(shù)Fig.7 Forecasted extreme temperature indexes in November 2022 based on the RF models established in the Yellow River Basin

3 討論

極端氣溫事件受多種大氣環(huán)流因子的綜合影響。丁之勇等[14]發(fā)現(xiàn)我國(guó)北疆地區(qū)極端氣溫指數(shù)受AO的影響最大,其次是NAO,而ENSO對(duì)其影響較小。龔道溢等[29]認(rèn)為NAO主要影響北美、西歐以及北非個(gè)別地區(qū)的氣溫,北太平洋濤動(dòng)(north Pacific oscillation, NPO)主要影響北太平洋和北美個(gè)別地區(qū)的氣溫,南方濤動(dòng)指數(shù)(SOI)則主要對(duì)太平洋熱帶地區(qū)氣溫產(chǎn)生影響。但過(guò)去研究多集中于單個(gè)或幾個(gè)環(huán)流指數(shù)對(duì)研究區(qū)域內(nèi)極端氣溫事件的影響,而本研究篩選了對(duì)極端氣溫事件有顯著影響的多個(gè)環(huán)流指數(shù),充分探究極端氣溫事件對(duì)多種大氣環(huán)流因子的綜合響應(yīng),本文從130個(gè)環(huán)流指數(shù)中,通過(guò)數(shù)據(jù)預(yù)處理、VIF多重共線性分析、考慮滯后性的Pearson相關(guān)分析及顯著性檢驗(yàn)等方法,嚴(yán)格地篩選出6個(gè)ETI的關(guān)鍵環(huán)流指數(shù),最終構(gòu)建SMLR和RF模型并進(jìn)行環(huán)流指數(shù)重要性排序,克服了大氣環(huán)流因子選取不客觀、不全面的局限性。

大氣環(huán)流對(duì)極端氣溫的影響具有滯后性。汪寧等[30]研究發(fā)現(xiàn)我國(guó)北部、東北部、南部地區(qū)的氣溫和降水對(duì)歐亞遙相關(guān)型的響應(yīng)約有2～4 d的滯后時(shí)間。黃浩等[31]發(fā)現(xiàn)AO和ENSO對(duì)河?xùn)|地區(qū)冬、夏兩季的極端氣溫指數(shù)有顯著影響,且呈現(xiàn)出0～1 a的時(shí)滯特征。馬乃孚等[32]提出北太平洋海溫與華中夏季氣溫存在較好的滯后相關(guān)關(guān)系,通過(guò)赤道東太平洋海溫異常預(yù)測(cè)華中地區(qū)夏季氣溫的可預(yù)報(bào)尺度為1～3個(gè)月。本文發(fā)現(xiàn)環(huán)流指數(shù)對(duì)黃河流域極端氣溫指數(shù)的影響具有0～12個(gè)月的時(shí)滯特征,且不同環(huán)流指數(shù)的最佳滯后時(shí)間也有差異,這與黃浩等[31]的結(jié)果部分一致。同時(shí)本研究以時(shí)滯特征為切入點(diǎn)計(jì)算了不同環(huán)流指數(shù)的最佳滯后時(shí)間,并通過(guò)最佳滯后時(shí)間序列數(shù)據(jù)建立模型預(yù)測(cè)了黃河流域ETI,考慮了環(huán)流指數(shù)對(duì)ETI影響的時(shí)滯性,從而提高了模擬的效率和精度。

黃河流域ETI影響因子的重要性在不同站點(diǎn)上具有差異。RF模型具有良好的模擬性能,已被應(yīng)用于很多領(lǐng)域。本研究采用2種模型對(duì)黃河流域ETI進(jìn)行實(shí)時(shí)預(yù)測(cè),對(duì)比分析發(fā)現(xiàn)在相同輸入下,RF模型除對(duì)TXx的模擬效果稍弱外,對(duì)其他5個(gè)ETI的模擬效果均優(yōu)于SMLR模型。由于RF模型能夠給出單個(gè)站點(diǎn)的變量重要性排序,且環(huán)流對(duì)ETI的影響具有區(qū)域差異[33],本文基于RF模型給出的變量重要性指數(shù)IncMSE,創(chuàng)新性地提出了影響整個(gè)黃河流域80個(gè)站點(diǎn)的各ETI的環(huán)流指數(shù)重要性排序方法RankYR,并進(jìn)行了具體分析,從而擴(kuò)展了重要性排序的應(yīng)用領(lǐng)域,得出了對(duì)黃河流域極端氣溫事件提供準(zhǔn)確可靠預(yù)報(bào)的新方法,為將來(lái)研究分析黃河流域氣候變化中極端氣溫變化可能造成的影響提供了參考。

4 結(jié)論

本文基于關(guān)鍵環(huán)流指數(shù)的最佳滯后時(shí)間序列數(shù)據(jù)構(gòu)建了6個(gè)極端氣溫指數(shù)的SMLR和RF模型,探究了極端氣溫事件對(duì)多種大氣環(huán)流因子的綜合響應(yīng),得出以下結(jié)論:

1)黃河流域1961—2020年的ETI中,TXx、TX90p、TD30和TNn呈波動(dòng)上升的趨勢(shì),而TN10p和FD0呈現(xiàn)下降趨勢(shì);極端高溫事件的強(qiáng)度和發(fā)生頻率的空間分布特征與極端低溫事件基本相反。

2)以靖遠(yuǎn)站TXx為例,極端氣溫事件受多種環(huán)流指數(shù)的綜合影響。環(huán)流指數(shù)對(duì)TXx的影響存在滯后性和明顯的周期性。關(guān)鍵環(huán)流指數(shù)對(duì)極端氣溫指標(biāo)具有顯著影響,其rmax對(duì)應(yīng)的滯后時(shí)間主要為5、6、11和12個(gè)月。同時(shí),環(huán)流指數(shù)對(duì)ETI的影響程度存在差異,PPVI對(duì)TNn和FD0的影響最大,IOWPA對(duì)TX90p和TD30影響最大,AEPVI和NANRP對(duì)TXx和TN10p影響最大。對(duì)整個(gè)黃河流域來(lái)說(shuō),PPVI對(duì)TXx、TNn、TX90p和FD0影響最大,而NANRP對(duì)TN10p和TD30影響最大。

3)SMLR和RF模型對(duì)ETI的預(yù)測(cè)能力較好,驗(yàn)證期R2分別為0.53～0.95和0.64～0.95。RF模型除對(duì)TXx的模擬效果稍差外,對(duì)其他5個(gè)ETI的模擬效果均優(yōu)于SMLR模型,因此被選擇用于預(yù)測(cè)黃河流域多站點(diǎn)的實(shí)時(shí)ETI。預(yù)測(cè)的2022年11月,黃河流域ETI的空間分布特征與多年平均情況基本相似。極端高溫事件的強(qiáng)度(TXx)和發(fā)生頻率(TX90p和TD30)較高的站點(diǎn)主要集中在流域下游的低海拔地區(qū)和局部河谷地帶,極端低溫事件強(qiáng)度和發(fā)生頻率的空間分布規(guī)律與基本極端高溫事件相反。