基于空間貝葉斯層次模型的淮河流域氣候極值特征分析

王懷軍 ，潘瑩萍 ，馮 如，肖明賢

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

0 引 言

【研究意義】氣候變化背景下，極端氣候事件發(fā)生的頻率和強(qiáng)度增加，對(duì)生態(tài)系統(tǒng)、人類健康、糧食生產(chǎn)和基礎(chǔ)設(shè)施安全均會(huì)產(chǎn)生重大影響[1-2]。全球范圍內(nèi)，觀測(cè)到的氣候變化記錄表明，冷日天數(shù)和冷夜天數(shù)正在逐漸減少，暖日天數(shù)和暖夜天數(shù)正在逐漸增加[3]。在歐洲、亞洲和澳大利亞的大部分地區(qū)，熱浪的頻率呈增加趨勢(shì)[4]。對(duì)于極端強(qiáng)降水，北半球中高緯度地區(qū)表現(xiàn)為增加趨勢(shì)，但這種變化不像極端溫度一樣，具有全球一致性[5]。一般來(lái)說(shuō)，與世界上降水量減少的地區(qū)相比，大多數(shù)地區(qū)的強(qiáng)降水事件仍在增加[6]。極端氣候事件的時(shí)空規(guī)律對(duì)于深刻認(rèn)識(shí)氣候?yàn)?zāi)害演變過(guò)程和防災(zāi)減災(zāi)具有重要的意義?！狙芯窟M(jìn)展】站點(diǎn)尺度上的極值頻率分析往往局限于有限的觀測(cè)數(shù)據(jù)[7]。在傳統(tǒng)地統(tǒng)計(jì)學(xué)中，常使用半變異函數(shù)來(lái)確定空間相關(guān)性并使用線性預(yù)測(cè)進(jìn)行變量空間插值[8]；然而，傳統(tǒng)地統(tǒng)計(jì)學(xué)不能很好地用于分析非高斯分布的極端氣候數(shù)據(jù)[7]。同時(shí)考慮多個(gè)站點(diǎn)的極端氣候事件分布，則需要將極值理論擴(kuò)展到多維變量[9]。最常見(jiàn)的空間極值模型有Max-stable 模型，它是多元極值分布的無(wú)限維泛化[10-11]。盡管Max-stable 模型能很好地模擬極端事件的空間結(jié)構(gòu)，但其模擬的各重現(xiàn)水平空間分布過(guò)于平滑，模擬精度下降[12]。此外，Max-stable 模型一般適用于描述年最大值序列（AM），而對(duì)于超門(mén)限峰值序列（POT），則需要選擇其他空間極值模型[13-14]?？臻g貝葉斯層次模型（Spatial Bayesian hierarchical model）是量化上述不確定性的較好方法[15-17]。例如，Gaetan[18]提出了一個(gè)用于分析意大利極端降雨的空間貝葉斯層次模型。Schliep 等[19]利用空間貝葉斯模型分析了區(qū)域氣候模式輸出所產(chǎn)生的極端降水量。Hrafnkelsson 等[16]利用貝葉斯層次模型分析了年最低氣溫和最高氣溫?？臻g貝葉斯層次模型在極端事件中的運(yùn)用還包括極端海浪高度[20-21]，高溫?zé)崂薣14]，地震[22]，和風(fēng)速[23]等。上述研究均表明，空間貝葉斯模型可以很好地捕捉極端事件的時(shí)空變化。

【切入點(diǎn)】淮河流域地處我國(guó)南北氣候過(guò)渡帶，是我國(guó)重要的地理生態(tài)分界線和生態(tài)環(huán)境脆弱區(qū)。極端氣候事件引發(fā)的氣象災(zāi)害已經(jīng)成為制約流域經(jīng)濟(jì)、生態(tài)、社會(huì)可持續(xù)發(fā)展的因素之一。據(jù)統(tǒng)計(jì)，1949年至今，淮河流域發(fā)生了數(shù)10 次全流域性洪澇和干旱災(zāi)害。其中，1991、2003、2007 年的洪澇災(zāi)害造成的直接經(jīng)濟(jì)損失分別達(dá)339.6 億、286.0 億和155.2億元；1949—2010 年，淮河流域累計(jì)干旱受災(zāi)面積達(dá)1.67 億hm2，損失糧食13.96 億kg，平均每年有269.8 萬(wàn)hm2農(nóng)作物受旱災(zāi)影響[24]。目前，淮河流域的氣候極值時(shí)空研究大多是基于站點(diǎn)尺度的趨勢(shì)分析，通過(guò)空間插值得到其極端氣候的空間分布[25-26]。然而，經(jīng)典地統(tǒng)計(jì)學(xué)理論基于多元正態(tài)分布，這并不能完全適用于具有尾部特征的極端數(shù)據(jù)建模[27]。且1 次極端氣候事件會(huì)影響多個(gè)站點(diǎn)，因此需要考慮站點(diǎn)之間的空間相關(guān)性[11]。單變量極值函數(shù)可用于描述極端數(shù)據(jù)，但實(shí)現(xiàn)高維數(shù)據(jù)的空間極值建模是一項(xiàng)挑戰(zhàn)?！緮M解決的關(guān)鍵問(wèn)題】本文將極值理論與地統(tǒng)計(jì)學(xué)理論相結(jié)合，評(píng)估了空間貝葉斯層次模型對(duì)淮河流域極端氣候事件的模擬能力。系統(tǒng)分析了極端氣候事件不同重現(xiàn)期下重現(xiàn)水平的空間分布，以期為極端氣候事件時(shí)空分析提供一種新的思路。

1 研究方法與數(shù)據(jù)來(lái)源

1.1 研究區(qū)概況

淮河流域地處中國(guó)東部，位于長(zhǎng)江與黃河2 大流域之間，經(jīng)緯度為111°55'—121°25'E，30°55'—36°36'N（圖1）。流域東西長(zhǎng)約700 km，南北寬約400 km，面積為2.7×105km2?；春恿饔蛭鞑?，南部和東北部是山地和丘陵地區(qū)，約占流域總面積的1/3，其余地區(qū)為廣闊的平原地帶?；春恿饔蚴侵袊?guó)重要的南北氣候過(guò)渡地帶，淮河以北為暖溫帶氣候區(qū)，淮河以南為北亞熱帶氣候區(qū)。流域南部年降水量達(dá)1 300～1 400 mm，流域北部減少至600～700 mm。1 年中6—9 月是雨季，占全年降水量的60%～70%。地表徑流年平均深度約為 240 mm，年徑流量為 622×108m3, 屬于我國(guó)水資源短缺地區(qū)。流域種植面積約為0.14×108km2，人口1.8 億（2010 年），流域內(nèi)主要種植作物有小麥、水稻、玉米、薯類、大豆、棉花和油菜。

1.2 數(shù)據(jù)來(lái)源

所需的溫度、降水?dāng)?shù)據(jù)來(lái)自中國(guó)氣象數(shù)據(jù)服務(wù)中心（http://data.cma.cn/en），時(shí)間跨度為1960—2015年。采用標(biāo)準(zhǔn)正態(tài)均一性檢驗(yàn)（SNHT）對(duì)日數(shù)據(jù)進(jìn)行一致性檢驗(yàn)，并修正。選用RclimDex 軟件計(jì)算模型所需的日最大降水量（RX1day）、日最高氣溫（TXx）分別作為降水極值和溫度極值。最終選擇淮河流域及其周邊的43 個(gè)氣象臺(tái)站作為研究對(duì)象（圖1）。數(shù)字高程模型數(shù)據(jù)（DEM）來(lái)自地理空間數(shù)據(jù)云（http://www.gscloud.cn/），分辨率為30 m×30 m，為了提高計(jì)算速度，將其重新采樣到0.15°×0.15°（圖1）。

圖1 淮河流域氣象站點(diǎn)分布 Fig.1 Study area and meteorological station distributed in Huaihe river basin

1.3 計(jì)算方法

空間貝葉斯層次模型的推斷過(guò)程如下[28-29]。Y(s)是空間區(qū)域的組最大值（Block maxima）過(guò)程（本文為1日為最大降水量（RX1day）和日最高氣溫（TXx）），假設(shè)這個(gè)過(guò)程是Max-stable過(guò)程，那么所有的邊際分布必然服從廣義極值分布（GEV）即：

式中μ(s)、ξ(s)、σ(s)分別為位置，形狀和尺度參數(shù)。

位置和尺度參數(shù)在空間上變化，故GEV位置參數(shù)和尺度參數(shù)可由（2）式表示：

式中：c1(s)、c2(s)、c3(s)分別表示經(jīng)度、緯度和高程協(xié)變量，β為線性方程系數(shù)。由于模型是通過(guò)多個(gè)層次結(jié)構(gòu)化，貝葉斯推斷是首選，模型后驗(yàn)參數(shù)由馬爾可夫鏈蒙特卡羅（MCMC）算法得到。

由于GEV 參數(shù)是隱高斯過(guò)程建模（ Latent Gaussian processes），目標(biāo)位置Sn的GEV邊緣分布參數(shù)用普通克里金法插值得到，每次 MCMC 迭代將通過(guò)克里金法插值得到一個(gè)目標(biāo)位置 GEV 參數(shù)，取所有迭代次數(shù)的均值作為目標(biāo)位置的最終參數(shù)。 假定目標(biāo)位置Sn的GEV最終參數(shù)為 μ、σ、ξ，則T 年重現(xiàn)水平可由式（3）計(jì)算：

2 結(jié)果與分析

2.1 空間貝葉斯層次模型的建立與驗(yàn)證

淮河流域氣候極值在站點(diǎn)尺度上的廣義極值分布（GEV）參數(shù)與站點(diǎn)環(huán)境變量（高程、經(jīng)度和緯度）之間的相關(guān)系數(shù)如表1所示。其中TXx的GEV位置參數(shù)（Loc）與經(jīng)度顯著相關(guān)，這表明經(jīng)度對(duì)TXx均值的空間分布具有顯著影響；尺度參數(shù)（Scale）與緯度、經(jīng)度均顯著相關(guān)，表明TXx的方差分布同時(shí)受經(jīng)緯度的影響。而對(duì)于RX1day，其尺度參數(shù)只與海拔顯著相關(guān)；同時(shí)，RX1day的位置參數(shù)與高程、經(jīng)度和緯度均存在顯著相關(guān)性。對(duì)于RX1day和TXx，其形狀參數(shù)（Shape）與高程、緯度和經(jīng)度相關(guān)性均不強(qiáng)，這意味著該參數(shù)對(duì)氣候極值空間分布及變化的影響較小。因此，我們?cè)诳臻g極值建模時(shí)，在邊際函數(shù)中將其設(shè)定為常數(shù)。

由上述結(jié)果確定海拔、經(jīng)度和緯度均被納入RX1day和TXx空間貝葉斯層次模型中。在空間極值建模過(guò)程中，MCMC迭代次數(shù)為5 000次，通過(guò)GEV各參數(shù)（Loc、Scale、Shape）的跡線圖來(lái)評(píng)估模型的收斂性（圖2）。當(dāng)?shù)螖?shù)為400時(shí)，TXx的各GEV參數(shù)基本穩(wěn)定下來(lái)（圖2（a）—圖2（d））。與TXx相比，RX1day的迭代次數(shù)達(dá)到2 400才能保持GEV各參數(shù)跡線圖平穩(wěn)（圖2（e）—圖2（h））。TXx和RX1day的對(duì)數(shù)似然比也表明空間貝葉斯層次模型可以很好地模擬空間極值的變化，400和2 400的迭代可以使TXx和RX1day參數(shù)收斂到真實(shí)值，因此取400和2 400以后的參數(shù)均值作為模型參數(shù)。最終計(jì)算的模型參數(shù)如表2所示。

表1 觀測(cè)數(shù)據(jù)GEV參數(shù)與高程、緯度、經(jīng)度的相關(guān)系數(shù) Table1 The correlation coefficients between GEV parameters of observed data and elevation, latitude and longitude

表2 TXx 和RX1day 建立的空間貝葉斯極值模型參數(shù) Table2 Estimated parameters of the spatial Bayesian hierarchical model fitted by the TXx and RX1day data set

圖2 空間貝葉斯層次模型模擬的GEV參數(shù)與對(duì)數(shù)似然比隨迭代次數(shù)變化跡線圖 Fig. 2 Variation of the GEV parameter and the log likelihood ratio with the number of iterations simulated by the spatial Bayesian hierarchical model

為評(píng)估模型的準(zhǔn)確性，將空間極值模型參數(shù)按站點(diǎn)提取，并與直接基于GEV函數(shù)分析的站點(diǎn)參數(shù)進(jìn)行對(duì)比（圖3）。結(jié)果顯示，TXx參數(shù)之間具有很強(qiáng)的相關(guān)性，位置參數(shù)和尺度參數(shù)的決定系數(shù)分別達(dá)到了0.98和0.80（圖3（a）—圖3（b））；對(duì)于RX1day，位置和尺度參數(shù)的決定系數(shù)分別為0.79和0.24（圖3（c）—圖3（d）），這表明空間極值模型可以很好地模擬各站點(diǎn)的降水/溫度極值變化。對(duì)空間極值模型的結(jié)果按站點(diǎn)進(jìn)行提取，與站點(diǎn)觀測(cè)數(shù)據(jù)亦進(jìn)行了比較，每個(gè)臺(tái)站模擬的次數(shù)為5 000次（圖4，圖5）。對(duì)于TXx, 所有觀測(cè)數(shù)據(jù)都包含在模擬數(shù)據(jù)中，表明空間貝葉斯層次模型結(jié)果可以很好地模擬觀測(cè)數(shù)據(jù)（圖4）。同樣，空間貝葉斯層次模型亦可以較好的模擬RX1day，但模擬效果總體上不如TXx（圖5）。

2.2 氣候極值不同重現(xiàn)水平的空間分布特征

由建立的空間貝葉斯層次模型計(jì)算淮河流域氣溫與降水極值不同重現(xiàn)期重現(xiàn)水平的空間分布（圖6，圖7）。對(duì)于氣溫極值TXx，其重現(xiàn)水平顯示出2種分布模式，即流域東部重現(xiàn)水平相對(duì)較低，而流域西部重現(xiàn)水平高（圖6），這表明TXx具有典型的經(jīng)向地帶性，這可能與距離海洋的遠(yuǎn)近有關(guān)。流域西部高海拔地區(qū)也有較低的重現(xiàn)水平，這可能是高程相依性的影響。RX1day重現(xiàn)水平的空間格局與TXx有所不同（圖7），流域西北部為低值區(qū)，流域東南為高值區(qū)。值得注意的是，西部高海拔地區(qū)具有相對(duì)較高的重現(xiàn)水平，表明海拔高度同樣對(duì)降水極值分布有重要影響。

圖3 GEV參數(shù)與空間貝葉斯層次模型GEV參數(shù)散點(diǎn)對(duì)比 Fig. 3 The scatter plot of GEV parameter and that derived from the spatial Bayesian hierarchical model at individual station

圖5 淮河流域站點(diǎn)RX1day模擬數(shù)據(jù)（藍(lán)色散點(diǎn)）和觀測(cè)數(shù)據(jù)（紅色散點(diǎn)）對(duì)比散點(diǎn)圖 Fig. 5 The simulated RX1day (Blue scatter) and observed RX1day (Red scatter) in certain stations in Huaihe River basi

圖6 淮河流域TXx不同重現(xiàn)期重現(xiàn)水平空間分布 Fig. 6 Spatial pattern of return level for TXx in Huaihe River basin

圖7 淮河流域RX1day不同重現(xiàn)期重現(xiàn)水平空間分布 Fig. 7 Spatial pattern of return level for RX1day in Huaihe River basin

3 討 論

溫度極值TXx流域具有典型的經(jīng)向地帶性（圖6），這可能與距離海洋的遠(yuǎn)近有關(guān)。海洋和陸地的物理性質(zhì)不同，對(duì)太陽(yáng)輻射的吸收和發(fā)射以及熱容量大小都有顯著影響。陸地比海洋的熱容量小，夏季升溫快，溫度比海洋高，且距海洋越近，受海洋的影響越大。流域東部TXx總體偏低，自東往西逐漸升高，反映了海洋對(duì)流域東部地區(qū)氣溫的調(diào)節(jié)作用[30]。流域西部高海拔地區(qū)也有較低的重現(xiàn)水平，這可能是氣溫直減率的影響，氣溫梯度成為這一響應(yīng)的關(guān)鍵參數(shù)，在對(duì)流層，地形對(duì)溫度的影響主要為氣溫隨海拔升高而降低，且海拔越高下降率越大。RX1day的重現(xiàn)水平在流域西北部為低值區(qū)，流域東南為高值區(qū)（圖7），這與淮河流域主要受東亞季風(fēng)影響有關(guān)。由于淮河地區(qū)處于強(qiáng)季風(fēng)區(qū)，降水主要發(fā)生在夏季。在從流域東南進(jìn)入西北過(guò)程中，隨著距離的增加，受夏季風(fēng)的影響越來(lái)越小，所以降水量由東南向西北遞減。值得注意的是，西部高海拔地區(qū)具有相對(duì)較高的重現(xiàn)水平，表明海拔高度同樣對(duì)降水極值分布有重要影響。不同海拔和地形因素，對(duì)降水影響有很大差別，這種差別主要體現(xiàn)在降水隨海拔的變化率上。在山區(qū)，降水量受海拔影響明顯，在一定范圍海拔高度上，降水隨高度增加而增加。大量觀測(cè)研究表明，不同海拔、不同地形因素（坡向和坡度）的山區(qū)，每上升100 m，降水量增加24.9～144.9 mm，該變化率涵蓋了大氣環(huán)流、地理經(jīng)、緯度、海拔和地形因素對(duì)降水的影響[31]。

空間極值模型的不同重現(xiàn)水平能很好地?cái)M合單站點(diǎn)GEV分布擬合的重現(xiàn)水平，特別是TXx，其所有的重現(xiàn)水平?jīng)Q定系數(shù)都達(dá)到了0.9以上（圖8）。RX1day的決定系數(shù)相對(duì)較小，但亦達(dá)到了顯著水平（圖9）。另外，所有回歸方程的斜率都不等于1，這表明與實(shí)際重現(xiàn)水平相比，貝葉斯模擬的重現(xiàn)水平被低估或高估。推斷其他因素，如人類活動(dòng)，距海洋的遠(yuǎn)近可能影響極端氣候的空間格局，然而一本研究并沒(méi)有考慮這些因素，這可能導(dǎo)致重現(xiàn)水平的低估或高估。此外，隨著全球氣候變化，極端氣候事件的時(shí)間序列顯示為非平穩(wěn)性，在建模過(guò)程中亦沒(méi)有考慮非平穩(wěn)性對(duì)時(shí)空格局的影響。今后的研究應(yīng)進(jìn)一步考慮更多的協(xié)變量和將非平穩(wěn)性納入貝葉斯空間模型框架中。

圖8 淮河流域單站點(diǎn)TXx GEV不同重現(xiàn)期重現(xiàn)水平和相同位置空間貝葉斯模型重現(xiàn)水平散點(diǎn)圖 Fig. 8 The scatter plot of return level of GEV and simulated return level derived from spatial Bayesian hierarchical model for TXx at individual station in Huaihe River basin

圖9 淮河流域單站點(diǎn)RX1day GEV不同重現(xiàn)期重現(xiàn)水平和相同位置空間貝葉斯模型重現(xiàn)水平散點(diǎn)圖 Fig. 9 The scatter plot of return level of GEV and simulated return level derived from the spatial Bayesian hierarchical model for RX1day at individual station in Huaihe River basin

4 結(jié) 論

1）氣候極值在站點(diǎn)尺度上的廣義極值分布（GEV）位置參數(shù)和尺度參數(shù)與環(huán)境變量（高程，經(jīng)度和緯度）之間具有顯著相關(guān)關(guān)系。在一定迭代次數(shù)情況下，GEV各參數(shù)跡線圖均能達(dá)到收斂?？臻g極值模型可以很好地?cái)M合單站點(diǎn)GEV分布參數(shù)、重現(xiàn)水平以及站點(diǎn)觀測(cè)數(shù)據(jù)。

2）淮河流域東部氣溫極值重現(xiàn)水平較低，西部重現(xiàn)水平較高；對(duì)于降水極值RX1day，流域西北具有較低重現(xiàn)水平，流域東南具有較高重現(xiàn)水平。