基于Max-stable模型的淮河流域氣候極值空間建模分析

王懷軍，趙卓怡，曹 蕾,潘瑩萍，，馮 如，楊雅雪

(1. 淮陰師范學(xué)院城市與環(huán)境學(xué)院，江蘇 淮安 223300；2. 南京水利科學(xué)研究院水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029；3. 南京水利科學(xué)研究院水利部應(yīng)對(duì)氣候變化研究中心，江蘇 南京 2100294；4. 北京師范大學(xué)地理學(xué)部，北京 100875)

全球氣候變化，特別是極端氣候事件的頻率加劇，可能對(duì)人類生命、財(cái)產(chǎn)和生態(tài)系統(tǒng)的可持續(xù)發(fā)展造成嚴(yán)重的影響[1- 3]。全球變暖導(dǎo)致極端氣候事件的頻率、強(qiáng)度、空間范圍和持續(xù)時(shí)間發(fā)生改變，也可能導(dǎo)致前所未有的極端氣候事件[4]。聯(lián)合國減災(zāi)辦公室(UNISDR)和災(zāi)害流行病學(xué)研究中心(CRED)報(bào)告顯示：1995—2015年，有33.94億人受到與天氣有關(guān)的災(zāi)害(如洪水、極端溫度和干旱)的影響，更嚴(yán)重的是造成34.3萬人死亡[5]。研究者試圖利用各種方法來解釋全球變暖對(duì)極端氣候的影響，且多數(shù)研究側(cè)重通過趨勢分析探測極端氣候的時(shí)空變化規(guī)律[6]。結(jié)果顯示，世界上大部分地區(qū)暖日和暖夜呈上升趨勢，冷日和冷夜呈下降趨勢[7]。同時(shí)，很多研究也指出極端氣候變化的不確定性遠(yuǎn)遠(yuǎn)超過平均氣候變化的不確定性[8]，這可能與空間數(shù)據(jù)建模以高斯過程實(shí)現(xiàn)有關(guān)[9]。此外，多元正態(tài)分布是為平均態(tài)而設(shè)計(jì)，其會(huì)低估極端氣候的概率，因此不適用于極值數(shù)據(jù)建模[10]。

近年來，Max-stable在極值空間建模中的應(yīng)用受到了相當(dāng)大的關(guān)注，因?yàn)樗嵌嘣獦O值分布的一種自然擴(kuò)展[11]。Max-stable模型在極端降水[12]和極端溫度[13]已有報(bào)道，結(jié)果均表明該模型可以很好地模擬氣候極值的空間變化特征。例如，Padoan等[14]和Davis等[15]認(rèn)為Max-stable可以描述極值數(shù)據(jù)的空間結(jié)構(gòu)和特征。Lee等[16]開發(fā)了一個(gè)通用的Max-stable模型，并演示了該模型在韓國年最高溫度變化的模擬能力。Zhang等[17]表明降水極值的空間分布受地形影響，增加海拔作為協(xié)變量可以提高M(jìn)ax-stable模型的估計(jì)精度。

淮河流域位于中國南北向氣候過渡帶，是中國南水北調(diào)工程實(shí)施的重要區(qū)域。近幾十年來，該地區(qū)極端降水的頻率和強(qiáng)度呈上升趨勢[18]。以往的研究大部分關(guān)注極端降水的概率分布特征。例如，She等[19]和Zheng等[20]觀察到，月降水量更適合于廣義帕累托分布。然而很少有研究考慮數(shù)據(jù)的空間相關(guān)性。由于極值空間分布和氣候變化風(fēng)險(xiǎn)是高度相關(guān)的，也是評(píng)估氣候變化背景下極值對(duì)環(huán)境危害影響的前提[21]。在這項(xiàng)研究中，我們研究了1日最大降水量(RX1day)和年最高氣溫(TXx)的Max-stable過程，以期通過選擇最優(yōu)空間極值模型來表征極端氣候的空間變化規(guī)律。

1 研究區(qū)域、數(shù)據(jù)和方法

1.1 研究區(qū)域和數(shù)據(jù)

淮河流域占地27萬km2，位于長江流域和黃河流域之間?；春恿饔蛭挥谖覈媳睔夂蜻^渡地帶，北起暖溫帶氣候，南至亞熱帶氣候；年平均氣溫11～16℃；月平均最高氣溫發(fā)生在7月，約為25℃，月平均最低氣溫發(fā)生在1月，約為0℃，最高氣溫達(dá)44.5℃。降水發(fā)生在6月至9月，年平均降水量約為920mm，一般由南向北、由山地向平原、由海岸向內(nèi)陸遞減。

淮河流域1960—2018年34個(gè)氣象站日值氣溫和降水?dāng)?shù)據(jù)來自中國氣象數(shù)據(jù)服務(wù)中心(http://data.cma.cn/en)。東亞夏季風(fēng)指數(shù)(EASM)來自http://ljp.gcess.cn/dctpage65577。氣候極值指數(shù)采用RClimDex軟件計(jì)算，其可以得到11個(gè)極端降水類指數(shù)和16個(gè)極端溫度類指數(shù)。本研究選擇RX1day和TXx分別為1日最大降水量和年極端最高氣溫。緯度、經(jīng)度、海拔和國內(nèi)生產(chǎn)總值(GDP)被視為Max-stable模型的空間協(xié)變量；緯度、經(jīng)度、海拔、GDP變化范圍分別為30.73°N-36.75°N、11.5°E-121.6°E、0～2160m、0～53178元。

1.2 方法

極端氣候事件趨勢的計(jì)算采用Mann-Kendall檢驗(yàn)[22- 23]，由于該檢驗(yàn)的結(jié)果容易受到時(shí)間序列自相關(guān)的影響，因此采用R語言包(“ZYP”)來消除lag- 1自相關(guān)的影響[24]。極端氣候事件的空間變化是通過Max-stable建模過程來計(jì)算，該模型充分考慮了極端氣候事件的空間相關(guān)性。Max-stable模型采用R語言包“SpatialExtremes”計(jì)算[25]，計(jì)算過程如圖1所示。

圖1 模型計(jì)算過程

下面簡要描述Max-stable建模過程。

隨機(jī)過程Z(·)是一個(gè)獨(dú)立同分布隨機(jī)場(Yi(x)，x∈Rd)上的最大值極限過程。如果存在常數(shù)an(x)>0、bn(x)∈R滿足

(1)

那么我們就能把Z(·)視為一個(gè)Max-stable過程，在區(qū)域極端氣候分析中，Z(·)可以看作是S個(gè)具有固定地理位置站點(diǎn)的氣候極值的一個(gè)隨機(jī)變量，其中用n表示時(shí)間段，一般為365天(1年)。同時(shí)，Z(·)的邊際分布遵循廣義極值分布(GEV)：

(2)

Pr[Z(x) ≤Z]=exp(-1/z),z>0

(3)

使用Brown-Resnick模型表示Max-stable過程的分布函數(shù)，其二元累積分布函數(shù)(CDF)表示如下：

Pr{Z(x1≤Z1,x1≤Z1)}=exp[-1/z1Φ(a/2+

1/alog(z2/z1))- 1/z2Φ(a/2+1/alog(z1/z2))]

(4)

式中，a=Var{Y(x1-x2)}，Y(·)是穩(wěn)定Gaussian過程；Φ—標(biāo)準(zhǔn)正態(tài)累計(jì)分布函數(shù)。

空間極值的相關(guān)性可以用極值系數(shù)(extremal coefficient)來表示，兩點(diǎn)之間的極值系數(shù)表示如下：

(5)

式中，γ—半變異方差；h—兩個(gè)地點(diǎn)之間的歐幾里得距離；極值系數(shù)的值在[1，2]中變化，1—完全相關(guān)，2—完全獨(dú)立。

由于傳統(tǒng)的似然方法需要對(duì)所有站點(diǎn)求取密度函數(shù)，因此不能直接用于Max-stable模型，研究采用成對(duì)似然法擬合Max-stable模型。

在考慮沒有觀測數(shù)據(jù)地區(qū)的空間極值時(shí)，可以將緯度、經(jīng)度和海拔的屬性數(shù)據(jù)視為GEV分布的協(xié)變量。本研究建立了16個(gè)以經(jīng)緯度為協(xié)變量的Max-stable模型，見表1。采用Takeuchi信息準(zhǔn)則(TIC)從以上模型中確定最佳擬合模型，TIC最小的模型確定為最佳擬合模型。然后，在最佳擬合模型中加入海拔、GDP和氣候指數(shù)，分析這些指數(shù)是否降低TIC值。如果TIC值降低，這些指標(biāo)將被認(rèn)為是Max-stable模型的協(xié)變量。重現(xiàn)水平計(jì)算如下：

(6)

式中，μ、σ、ξ—位置，尺度和形狀參數(shù)，yT—T年重現(xiàn)水平。

表1 廣義極值分布參數(shù)協(xié)變量模型

注：lon(x)和lat(x)分別是x站的經(jīng)度和緯度；alt(x)、GDP(x)分別是x站的高程和國內(nèi)生產(chǎn)總值；temp表示時(shí)間協(xié)變量；μ、σ分別表示位置和尺度參數(shù)。

圖2 中國淮河流域氣候極值平均值和趨勢強(qiáng)度

2 結(jié)果和分析

2.1 氣候極值空間變化

圖2展示了淮河流域氣候極值的空間均值和趨勢。TXx在35～39℃之間，且由東向西遞增，這是由于陸地升溫的速度比海洋快，離海岸越近，海洋對(duì)它的影響就越大。TXx趨勢為流域東南部呈上升趨勢，西北部呈下降趨勢。RX1day的均值范圍為80～115mm。由于受海洋和東亞季風(fēng)的影響，較西北地區(qū)而言，流域東南部的RX1day顯示出相對(duì)更高的值。RX1day的空間趨勢為-4～4mm/10年，但無顯著空間趨勢分布。TXx和RX1day均呈上升趨勢(TXx為0.02℃/10年，RX1day為0.55mm/10年)，但只有TXx在0.05顯著水平上顯著。

2.2 氣候極值空間建模

運(yùn)用單變量GEV模型對(duì)TXx和RX1day時(shí)間序列進(jìn)行擬合，并與環(huán)境參數(shù)進(jìn)行相關(guān)關(guān)系分析，見表2。對(duì)于TXx的位置參數(shù)，經(jīng)度相關(guān)系數(shù)為-0.67，說明TXx的空間分布與經(jīng)度密切相關(guān)，與圖2a2的空間分布趨勢具有一致性。此外，尺度參數(shù)與經(jīng)度之間也存在負(fù)相關(guān)關(guān)系，表明在經(jīng)度數(shù)值大的地方方差較小。與此相反，緯度與尺度參數(shù)呈現(xiàn)正相關(guān)系數(shù)，說明TXx在高緯度地區(qū)具有較高的方差。對(duì)于RX1day，位置參數(shù)與緯度/經(jīng)度之間的關(guān)系為負(fù)/正系數(shù)，表明高RX1day發(fā)生在低緯度、高經(jīng)度，這與圖2a1的結(jié)果近似。此外，RX1day的均值和方差隨著海拔的降低而增加。因此，可以推測經(jīng)度、緯度和海拔等變量可能影響淮河流域極端氣候事件的空間分布。

表2 站點(diǎn)GEV分布所得位置、尺度參數(shù)與海拔、緯度、經(jīng)度相關(guān)關(guān)系

為進(jìn)一步研究經(jīng)緯度對(duì)極端氣候事件的影響，以經(jīng)度和緯度作為Max-stable模型的協(xié)變量建立模型，選擇TIC值最低的模型為最佳模型，如3所示。M11是TXx和RX1day的最佳模型。將海拔加入到M11得到模型M17-M20中(表1)，結(jié)果顯示TIC值下降，說明海拔提高了模型模擬性能。除了經(jīng)度、緯度、海拔等因素外，人類活動(dòng)等其他因素也可能影響極端氣候事件的空間格局。此外，由于許多氣候時(shí)間序列是非平穩(wěn)的，加入非平穩(wěn)性的Max-stable模型可能比單一極值建模更為可靠[11]。因此，論文還關(guān)注了具有這些協(xié)變量(例如，反映人類活動(dòng)的GDP和氣候指數(shù))的Max-stable模型(表1)。從圖3可以看出，TXx加入GDP后(M21)TIC下降，而RX1day的TIC保持穩(wěn)定。夏季EASM作為GEV參數(shù)的時(shí)間協(xié)變量加入到Max-stable模型后，TIC值并沒有降低，表明添加時(shí)間協(xié)變量并不能提高模型的性能。因此，TXx選擇有TIC最低值9756的M21作為最佳模型，而RX1day選擇有TIC最低值26252的M20作為最佳模型，如圖3a、3b所示，見表3。

圖4展示了基于最優(yōu)Max-stable模型得到極值系數(shù)(TXx采用M21，RX1day采用M20)估計(jì)對(duì)(θ)的距離(h)散點(diǎn)圖。理論極值系數(shù)函數(shù)(用紅線表示)也表明所選擇的Max-stable模型適用于極端氣候事件的空間建模。除此之外，單個(gè)站點(diǎn)GEV分布和Max-stable模型估算的位置和尺度參數(shù)具有高度相關(guān)關(guān)系，特別是TXx，M20位置和尺度參數(shù)的決定系數(shù)分別達(dá)到0.66和0.32。TXxM21尺度參數(shù)的決定系數(shù)相比M20(0.32)增加到0.43，如圖5所示，說明人類活動(dòng)確實(shí)影響了TXx的空間格局。此外，我們還評(píng)估了單站點(diǎn)GEV重現(xiàn)水平和Max-stable模型重現(xiàn)水平在觀測站點(diǎn)的相關(guān)關(guān)系，見表4，所有決定系數(shù)均在0.05顯著水平上顯著。當(dāng)然，我們還注意到，當(dāng)TXx考慮GDP時(shí)，不同重現(xiàn)水平的決定系數(shù)均有所增加，而RX1day的決定系數(shù)未增加。這說明Max-stable模型能夠準(zhǔn)確模擬空間極值，而GDP這一協(xié)變量的添加提高了淮河流域TXx模擬的性能。

圖3 不同Max-stable模型下TXx和RX1day的TIC值散點(diǎn)圖

圖4 基于半變異函數(shù)的TXx和RX1day擬合極值系數(shù)隨距離變化圖

圖5 各站點(diǎn)TXx(a1-a4)和RX1day(b1-b4)Max-stable模型估計(jì)的GEV參數(shù)與最大似然估計(jì)估算的站點(diǎn)GEV參數(shù)關(guān)系散點(diǎn)圖

位置尺度TICM20(TXx)μ(x)=-1236.36-3.01 lat(x)-10.93 lon (x) +226.80lon 0.5(x) -0.0045 alt(x)+34.2 alt0.5 (x)σ(x)=66.62+2.01lat(x)+0.05lon0.5(x)+ 0.0003alt(x)-22.82alt0.5(x)9794M21(TXx)μ(x)=-1353.87-3.82 lat(x)-11.73 lon (x) +244.12lon 0.5(x) -0.0036 alt(x)+43.14 alt0.5 (x)+0.31GDPσ(x)=66.32+1.89lat(x)-0.49lon0.5(x)+ 0.0002alt(x)-20.66alt0.5(x)+0.12 GDP9756M20(RX1day)μ(x)=-22877.40-30.47 lat(x)-187.46lon (x) +4074.87lon 0.5(x) -0.0035 alt(x)+317.82 alt0.5 (x)σ(x)=891.95+25.55lat(x)+3.77lon0.5(x)-0.0116alt(x)-304.06alt0.5(x)26252M21(RX1day)μ(x)=-21585.65-22.10 lat(x)-178.67 lon (x) +3886.38lon 0.5(x) -0.0079 alt(x)+221.17 alt0.5 (x)-1.85GDPσ(x)=912.23+26.17lat(x)+3.59lon0.5(x)+ 0.0114alt(x)-310.87alt0.5(x)+0.084 GDP26252

表4 站點(diǎn)GEV重現(xiàn)水平與Max-stable模型對(duì)應(yīng)的站點(diǎn)模擬重現(xiàn)水平之間關(guān)系

2.3 氣候極值重現(xiàn)水平估計(jì)

利用選擇的最佳模型M20和M21計(jì)算Rx1day和TXx重現(xiàn)水平(2年、10年、50年和100年)的空間分布，如圖6—7所示。圖6b1—b4表明，TXx在考慮GDP后，大致空間分布沒有變化，但總體空間格局變得不平滑；圖6a1—b4表明以西北地區(qū)為代表的GDP低值地區(qū)重現(xiàn)水平出現(xiàn)下降。且觀測到的重現(xiàn)水平與模擬的重現(xiàn)水平平均絕對(duì)誤差空間差異更小，如圖8a1—b4所示，說明加入GDP后能夠更好地模擬TXx的空間格局。TXx 2年重現(xiàn)期的重現(xiàn)水平36 ～ 39℃，100年重現(xiàn)期的重現(xiàn)水平增長到39～43℃。不同重現(xiàn)期下的重現(xiàn)水平空間格局相似，東部為高溫地區(qū)，西部為低溫地區(qū)。這一結(jié)果表明，TXx主要表現(xiàn)為經(jīng)向分布，其空間格局主要由經(jīng)度決定。對(duì)于RX1day，當(dāng)考慮GDP時(shí)，空間平均絕對(duì)誤差MAE并沒有變小，如圖8c1—d4所示，這表明人類活動(dòng)對(duì)極端降水的影響較小。10年、50年和100年重現(xiàn)期的重現(xiàn)水平空間分布與2年重現(xiàn)期相似，都由西北向東南遞增，這個(gè)結(jié)果可以用經(jīng)緯度的綜合效應(yīng)來解釋。低緯度和高經(jīng)度地區(qū)易出現(xiàn)高RX1day值。由于淮河流域中氣象站數(shù)量較少，且氣象站之間的海拔差異較小，因此盡管將海拔納入了計(jì)算模型，研究中中海拔對(duì)極端氣候事件的空間格局的影響并不顯著。

圖6 TXx最佳Max-stable模型下的2、10、50和100年重現(xiàn)期重現(xiàn)水平空間分布(單位：℃)

圖7 RX1day最佳Max-stable模型下的2、10、50和100年重現(xiàn)期重現(xiàn)水平空間分布(單位：mm)

圖8 TXx(a1—b4)和RX1day(c1—d4)的GEV觀測數(shù)據(jù)與Max-stable模型模擬數(shù)據(jù)MAE重現(xiàn)水平空間差異

3 結(jié)論

本文基于Max-stable模型分析了淮河流域極端氣候事件的空間分布，主要結(jié)果如下：

(1)極值系數(shù)函數(shù)、Max-stable模型模型與GEV站點(diǎn)相關(guān)參數(shù)散點(diǎn)圖、平均絕對(duì)誤差(MAE)圖等標(biāo)準(zhǔn)表明，建立的Max-stable模型可以很好的擬合氣候極值的空間變化。緯度、經(jīng)度和海拔對(duì)極端氣候事件的空間模擬有重要影響。GDP可以提高溫度極值Max-stable模型的性能。

(2)TXx在不同重現(xiàn)期(2年、10年、50年和100年)下的重現(xiàn)水平呈經(jīng)向分布，且西向東逐漸升高。此外，溫度極值的空間分布也會(huì)受到以GDP為代表的人類活動(dòng)的影響。RX1day由西北向東南遞，這可以用經(jīng)度和緯度的綜合作用來解釋。

由于世界上許多地方的氣象數(shù)據(jù)沒有監(jiān)測，可以通過空間插值得到不同重現(xiàn)期的重現(xiàn)水平，其精度取決于插值方法和監(jiān)測站密度。Max-stable模型的一個(gè)優(yōu)點(diǎn)是，它可以在沒有測量數(shù)據(jù)的情況下，利用海拔、經(jīng)度、緯度等環(huán)境因素獲取重現(xiàn)水平。與以往的研究相似，本文使用的Max-stable模型只能得到重現(xiàn)水平的空間變化，而不能得到極端氣候事件的時(shí)空變化。此外，研究中使用的Max-stable模型僅適用于AM序列，而對(duì)于超門限序列(POT序列)還需要采用其它模型，如Bayesian層次模型進(jìn)行空間建模。目前，我們只增加了EASM和GDP作為協(xié)變量來代表氣候變化和人類活動(dòng)，更多的環(huán)境變量(如離海洋距離，ENSO)將在未來的工作中進(jìn)行深化。