基于SPI和Copula的湟水流域干旱趨勢研究

蘇夏羿，張 鑫，王 云，王文亞(西北農(nóng)林科技大學(xué)水利與建筑工程學(xué)院，陜西 楊凌 712100)

0 引 言

由于氣候變化導(dǎo)致的水資源在時間和空間上的重新分配，以及質(zhì)量和數(shù)量上的改變，區(qū)域旱澇災(zāi)害發(fā)生的頻率有所增加[1]。近年來全國旱澇災(zāi)害頻發(fā)，西北干旱地區(qū)更是首當(dāng)其沖。湟水河是黃河上游的一級支流，主要由湟水干流和大通河構(gòu)成，湟水干流全長為374 km，流域面積17 733 km2[2]。湟水流域位于黃土高原與青藏高原的交錯地帶，屬高原大陸性氣候，流域內(nèi)降水地區(qū)分布很不均勻，西北部高原地帶氣溫低、蒸發(fā)少，植被覆蓋率大，降水整體上從西北到東南部遞減，并且干、濕季節(jié)明顯，降水高度集中于5-9月[3]。盡管湟水流域面積僅占青海省土地面積的2.3%，卻養(yǎng)育了青海省56%的人口并貢獻了48%的GDP[4]。流域中部和東南部地勢相對較低，俗稱湟水谷地，土壤肥沃、農(nóng)業(yè)資源優(yōu)越，是青海省東部農(nóng)業(yè)區(qū)的主要生產(chǎn)基地，也是青海省政治、經(jīng)濟、文化、交通中心區(qū)域。楊發(fā)源[5]在對湟水流域研究時得出20世紀(jì)90年代后出現(xiàn)了氣溫升高、降水減少和蒸發(fā)量增大的干旱化趨勢的結(jié)論，可見，對湟水流域干旱趨勢的研究很有必要。

在眾多的干旱指標(biāo)當(dāng)中，標(biāo)準(zhǔn)降水指標(biāo)(SPI)計算簡單，能較好地反映干旱歷時和干旱嚴(yán)重程度，廣泛應(yīng)用于干旱監(jiān)測與評估中[6]。姚瑤[7]和劉蕾[8]已將SPI指數(shù)應(yīng)用于青海東部農(nóng)業(yè)區(qū)研究，通過與實際干旱事件的比對，得出SPI能反映該區(qū)氣象干旱的實際情況的結(jié)論。湟水流域中西部和南部均屬于青海農(nóng)業(yè)區(qū)，故本研究選取SPI作為氣象干旱檢測指標(biāo)，結(jié)合游程理論對降水序列進行干旱識別，獲得干旱歷時和相應(yīng)干旱嚴(yán)重程度序列。基于copula函數(shù)在極端水文事件中的廣泛應(yīng)用[9,10]，本文采用Copula函數(shù)將干旱歷時和干旱嚴(yán)重程度這兩個顯著相關(guān)但服從不同分布的極值事件聯(lián)合起來，計算干旱歷時和干旱嚴(yán)重程度中有一個大于某一定值的“或”聯(lián)合重現(xiàn)期和兩者同時大于等于某一定值的“且”聯(lián)合重現(xiàn)期，從而全面且系統(tǒng)地對湟水流域干旱趨勢分析，為當(dāng)?shù)氐母珊笛芯考八Y源管理提供科學(xué)有效的依據(jù)。

1 材料與方法

1.1 數(shù)據(jù)來源

本文選取湟水河流域及周邊17個氣象站的1960-2006年日降水序列作為研究對象，其中烏鞘嶺和皋蘭站數(shù)據(jù)來源于中國氣象數(shù)據(jù)共享網(wǎng)，其余15站數(shù)據(jù)均為氣象站實測資料。湟水河流域以及氣象站點所在位置如圖1所示。所有站的降水?dāng)?shù)據(jù)均通過了三性審查。

圖1 湟水流域以及氣象站點位置示意圖Fig.1 The schematic diagram of Huangshui River Basin and meteorological stations

1.2 干旱事件識別

不同時間尺度下的SPI值具有不同的物理意義，一般認(rèn)為6個月尺度的SPI值可以很好地反映氣象干旱。為了更好地研究湟水流域干旱的趨勢變化，本文采用SPI6對降水序列進行處理，具體計算過程參考《氣象干旱等級》[11]。SPI干旱等級劃分如表1所示。

表1 SPI干旱強度等級劃分Tab.1 Drought intensity grades of standardized precipitation index

本文選取干旱歷時與干旱嚴(yán)重程度來表征一場干旱事件，故需結(jié)合劃分的干旱等級和游程理論對降水序列進行干旱識別，即定義干旱歷時D是SPI6值持續(xù)低于-0.5的時間(包括極個別略大于-0.5的月份)，干旱嚴(yán)重程度S為干旱歷時內(nèi)SPI6超出-0.5絕對值的累計值[12]。

(1)

式中：Ii為時段i的SPI值。

計算可獲得干旱歷時序列和相應(yīng)的干旱嚴(yán)重程度序列，以及單變量的經(jīng)驗概率分布，此時為離散狀態(tài)，需通過函數(shù)擬合的辦法進行處理。本文采用指數(shù)分布對干旱歷時進行擬合，采用威布爾分布對干旱嚴(yán)重程度進行擬合。各個分布函數(shù)的參數(shù)均采用最大似然法進行參數(shù)估計，并使用Kolmogorov-Smirnov(K-S)法進行非參數(shù)檢驗。

1.3 Mann-Kendall 趨勢檢驗法

假設(shè)時間序列(x1,x2,…,xn)無趨勢，計算統(tǒng)計量S[13]，計算式如下：

(2)

假設(shè)各變量獨立同分布，則統(tǒng)計量S近似服從正態(tài)分布[13]。其期望和方差分別如下：

(5)

式中：m代表序列中元素相同的組個數(shù)；tj代表第j個元素相同組中所包含觀測值的個數(shù)。

構(gòu)造M-K統(tǒng)計量Uc：

(5)

假設(shè)序列無趨勢，當(dāng)樣本數(shù)量n比較大時，Uc近似服從標(biāo)準(zhǔn)正態(tài)分布，給定顯著性水平α，則可根據(jù)|Uc|與臨界值Uα/2的比較結(jié)果判定序列趨勢的統(tǒng)計顯著性[13]。本文取顯著性水平α=0.05。本文主要運用M-K法對17個氣象站12個月份的SPI6、干旱歷時以及干旱嚴(yán)重程度進行趨勢分析，經(jīng)計算表明這些數(shù)據(jù)均無明顯的自相關(guān)性，因此可以采用該方法對其進行趨勢分析。

1.4 Copula函數(shù)及其參數(shù)估計與選擇

Sklar定理[14]是Copula函數(shù)理論及其應(yīng)用的基礎(chǔ)，闡述了Copula函數(shù)構(gòu)造多元函數(shù)聯(lián)合分布的方法。由Sklar定理可知，F(xiàn)(x1,x2,…,xn)是隨機變量(X1,X2,…,Xn)的聯(lián)合分布函數(shù)，其邊緣分布函數(shù)分別F1(x1),F2(x2),…,Fn(xn)，則存在一個Copula函數(shù)C滿足：

F(x1,x2,…,xn)=C[F1(x1),F2(x2),…,Fn(xn)]

(6)

如果邊緣分布函數(shù)是連續(xù)的，那么C是唯一的。Copula函數(shù)可將不同邊緣分布的變量構(gòu)造成聯(lián)合分布函數(shù)，由邊際分布和copula函數(shù)來確定多變量的聯(lián)合分布函數(shù)。對于干旱歷時X與干旱嚴(yán)重程度Y，其邊緣分布分別為FX(x)與FY(y)，F(xiàn)(x，y)為其聯(lián)合分布。

Copula函數(shù)有很多種類，本文所使用的是只有一個參數(shù) 的二維阿基米德Copulas函數(shù)族，包括Gumbel Copulas，Clayton Copulas 和Frank Copulas，對湟水流域的干旱歷時和干旱嚴(yán)重程度的聯(lián)合分布進行擬合。這三種Copula函數(shù)的參數(shù)θ與干旱歷時和嚴(yán)重程度序列的Kendall相關(guān)系數(shù)系數(shù)τ之間有確切的計算表達式(表2)，可據(jù)此關(guān)系進行參數(shù)估算[15]。

表2 Kendall-τ與copula函數(shù)的參數(shù)θ對應(yīng)關(guān)系表Tab.2 The corresponding relation between the Kendall-τ and the parameter of Copula Function

本文評價上述3種Copula函數(shù)對干旱歷時和嚴(yán)重程度的擬合優(yōu)劣的指標(biāo)值采用均方根誤差(RMSE)：

(7)

式中：n為序列樣本數(shù)；pc為Copula函數(shù)在原始點的理論概率分布值；pi為聯(lián)合經(jīng)驗概率分布值。

優(yōu)選出擬合效果最好的Copula函數(shù)，即為能反映該流域干旱歷時和干旱嚴(yán)重程度特征的聯(lián)合分布函數(shù)。為更好反映湟水流域某次干旱事件的干旱歷時與干旱嚴(yán)重程度分別大于某一特定值發(fā)生的概率，本文計算對應(yīng)的重現(xiàn)期值，重現(xiàn)期計算公式如下:

(9)

式中:T(X>x, Y>y)為干旱歷時與干旱嚴(yán)重程度均超過特定值的重現(xiàn)期，下文稱為“且”重現(xiàn)期；T(X>x or Y>y)表示干旱歷時或者干旱嚴(yán)重程度超過某一特定值的重現(xiàn)期，稱為“或”重現(xiàn)期；P(X>x,Y>y)和P(X>x or Y>y)相應(yīng)的聯(lián)合概率密度。重現(xiàn)期大表示出現(xiàn)干旱風(fēng)險相對較低，重現(xiàn)期小表示出現(xiàn)干旱風(fēng)險相對較高。

2 結(jié)果與分析

2.1 單站逐月SPI趨勢分析

當(dāng)M-K趨勢檢驗法統(tǒng)計量的絕對值大于1.96 時，說明該序列的趨勢在95%的置信區(qū)間內(nèi)顯著。單點M-K統(tǒng)計量小于零表示該站該月的SPI6序列隨時間有減小的趨勢，反映在旱澇上是干旱的趨勢；大于零表示該站該月的SPI6序列隨時間有增大的趨勢，反映在旱澇上表示濕潤的趨勢。湟水河流域各個站點12個月份的SPI6用M-K趨勢分析結(jié)果如圖2所示。圖中17個氣象站的站點按照經(jīng)度的大小，從西到東依次排列。

圖2 湟水流域17站各月SPI序列趨勢分析圖Fig.2 Trend graph of SPI in each months at seventeen stations in Huangshui River Basin

由圖2可知，湟水流域西北部各站SPI6序列在各月M-K統(tǒng)計量均在-1.96～1.96之間，呈不明顯的變化趨勢；流域東南部除了華隆和民和站之外其他站點的SPI6序列，均出現(xiàn)了M-K統(tǒng)計量的絕對值大于1.96的情況，呈現(xiàn)顯著的變化趨勢，并且顯著變化的情況集中在樂都、互助、西寧三站，共占總體的55%。只考慮干旱情況時，流域南部的湟源、貴德以及東部的互助、華隆、樂都、民和站幾乎每個月份的SPI6值均呈現(xiàn)下降的干旱化趨勢。尤其是互助站，所有月份的SPI6序列M-K統(tǒng)計量都小于零，呈現(xiàn)干旱化趨勢；并且3月、11月和12月均小于-1.96，呈現(xiàn)顯著干旱化趨勢。

M-K統(tǒng)計量為大于零集中在6-9月，意味著隨時間有變濕潤趨勢；小于零集中在12月到次年3月，意味著隨時間有變干的趨勢，這與該流域降水高度集中在5-9月的背景相疊加，勢必會造成湟水流域出現(xiàn)干旱的季節(jié)更干旱，濕潤的季節(jié)更濕潤的情況，符合氣候變化對極端降雨和干旱事件的描述。尤其是樂都站，1月和2月 的M-K統(tǒng)計量均小于-2.3，呈顯著干旱化趨勢；5-7月SPI的M-K統(tǒng)計量均大于2.1，呈顯著濕潤化趨勢，將增加旱澇災(zāi)害的可能，不利于水資源管理。

2.2 干旱歷時和嚴(yán)重程度趨勢分析

結(jié)合SPI值以及游程理論進行干旱識別，獲得各站干旱歷時與干旱嚴(yán)重程度序列，干旱特征統(tǒng)計情況以及M-K統(tǒng)計量見表3。單站干旱歷時的M-K統(tǒng)計量小于零表示該站干旱歷時減小的趨勢，旱情減弱的的趨勢；大于零表示該站干旱歷時隨時間有增大的趨勢，反映出旱情加重趨勢。干旱嚴(yán)重程度的M-K統(tǒng)計量表征能力同干旱歷時。

表3 湟水流域各站干旱特征統(tǒng)計表Tab.3 Statistical table of drought character at each stations in Huangshui River Basin

由表3可知，除了剛察和烏鞘嶺兩站，整個湟水流域干旱歷時的M-K統(tǒng)計量均在-1.26～1.36之間，干旱嚴(yán)重程度的M-K統(tǒng)計量均在-1.71～1.25之間，均無顯著變化趨勢。其中流域中部和西部地區(qū)干旱歷時和嚴(yán)重程度的M-K統(tǒng)計量大部分落在-1～0之間，呈現(xiàn)不顯著下降趨勢，相比之下，流域東部地區(qū)的M-K統(tǒng)計量大于零，意味干旱歷時和干旱嚴(yán)重程度有上升的趨勢。雖然沒有顯著干旱等級加重的趨勢，但湟水流域這些區(qū)域大部分為青海農(nóng)業(yè)區(qū)和牧場，依然需要注意輕度等級干旱對農(nóng)牧業(yè)帶來的損害，做好防范應(yīng)對措施。

利用湟水流域各站點干旱歷時與干旱嚴(yán)重程度序列的M-K統(tǒng)計量，通過ArcGIS10.0軟件的克里格差值法得到整個湟水流域干旱歷時與干旱嚴(yán)重程度的增減趨勢的空間分布圖(圖3)。由圖可以看出，干旱歷時與干旱嚴(yán)重程度的變化趨勢在空間分布上存在較好的一致性，歷時越長，相對應(yīng)的干旱嚴(yán)重程度越大。

圖3 湟水流域干旱歷時與干旱嚴(yán)重程度趨勢分析圖Fig.3 Trend graph of the drought duration and drought severity in Huangshui River Basin

2.3 典型干旱情景分析

借助SPSS軟件單參數(shù)K-S檢驗方法，湟水流域各個站點干旱歷時的指數(shù)分布和嚴(yán)重程度的威布爾分布均通過了K-S檢驗。計算兩序列的Kendall系數(shù)，再依據(jù)表1的表達式，對3種Copula函數(shù)的參數(shù)進行估算，并通過式(5)進行copula函數(shù)優(yōu)選，經(jīng)計算Frank Copula函數(shù)對經(jīng)驗聯(lián)合分布的擬合效果最好，所以選取Frank Copula作為理論聯(lián)合分布函數(shù)，計算干旱情景下的重現(xiàn)期。

針對不同重現(xiàn)期的計算，結(jié)合前人研究，選擇兩個典型干旱情景：干旱歷時大于6個月，基于SPI指數(shù)的干旱嚴(yán)重程度大于4表示為聯(lián)合中度干旱；干旱歷時大于8個月，基于SPI指數(shù)的干旱嚴(yán)重程度大于6為聯(lián)合重度干旱。再結(jié)合克里格差值法，便可得到整體湟水流域聯(lián)合中度干旱和中度干旱的重現(xiàn)期空間分布(圖4)。

圖4 湟水流域典型聯(lián)合重現(xiàn)期空間分布圖Fig.4 Spatial distributions of joint return periods under typical scenarios in Huangshui River Basin

從圖4(a)和圖4(b)上看，在發(fā)生聯(lián)合中度干旱的情況下，流域東南部的干旱重現(xiàn)期要低于西北部，因此東南部發(fā)生中度干旱的風(fēng)險要高于西北部。湟水流域中度干旱“或”聯(lián)合重現(xiàn)期平均為13.2年一遇[圖4(a)]，“且”聯(lián)合重現(xiàn)期平均為25.6年一遇[圖4(b)]。而位于東南部的民和站，“或”聯(lián)合重現(xiàn)期和“且”聯(lián)合重現(xiàn)期分別低至8.9年一遇和17.9年一遇。相比較下，西北部的門源和湟源站“或”和“且”聯(lián)合重現(xiàn)期均大于14.9年一遇和31.6年一遇，處于干旱低發(fā)區(qū)。

從圖4(c)和圖4(d)上看，在發(fā)生聯(lián)合重度干旱的情況下，與中度干旱相似，流域東南部的干旱重現(xiàn)期要低于西北部，東南部發(fā)生聯(lián)合重度干旱的風(fēng)險要高于西北部。湟水流域重度干旱“或”聯(lián)合重現(xiàn)期平均為23.7年一遇[圖4(c)]，“且”聯(lián)合重現(xiàn)期平均為59.7年一遇[圖4(d)]。最低的民和站“或”和“且”聯(lián)合重現(xiàn)期分別低至15.5年一遇和34.5年一遇，該站點附近區(qū)域重度干旱風(fēng)險較大，相比之下流域北部和中西部地區(qū)重度干旱風(fēng)險較小。

整體上看，無論對于中度干旱還是重度干旱，湟水流域東南部四種干旱聯(lián)合重現(xiàn)期均小于西北部，流域東南部發(fā)生旱災(zāi)的概率均大于西北部。流域北部區(qū)域即大通河上游區(qū)域，因獨特的地理位置，高山環(huán)繞，溫度低，蒸發(fā)量小的原因，無論是何種情境的重現(xiàn)期都相對比較高，發(fā)生干旱的可能性較河流下游小。中西部的門源和湟源區(qū)域，“或”重現(xiàn)期處在中等水平，而“且”重現(xiàn)期處于大值范圍，所以該區(qū)域易發(fā)生低等級干旱。流域東南部兩種情景的“或”與“且”重現(xiàn)期都比較低，所以兩種程度干旱均極易發(fā)生，應(yīng)作為流域的干旱敏感區(qū)域?qū)Υ?/p>

3 結(jié) 語

(1)本研究采用的標(biāo)準(zhǔn)化降水指數(shù)(SPI)對湟水流域17站進行干旱趨勢分析，從時間上看，湟水流域出現(xiàn)了干旱的季節(jié)更干旱，濕潤的季節(jié)更濕潤的兩極化情況。在氣候變化的背景下，湟水流域夏旱的情況可能會加重，不利于水資源管理。

(2)從空間上看，流域的東部與南部地區(qū)干旱化趨勢嚴(yán)重，但通過對干旱歷時和干旱嚴(yán)重程度趨勢分析，流域各站干旱歷時和干旱嚴(yán)重程度的上升趨勢不明顯，可以得出低等程度的干旱在湟水流域東南部發(fā)生的可能性增大的結(jié)論。

(3)通過湟水流域17個站點構(gòu)建Copula函數(shù)描述干旱歷時和干旱嚴(yán)重程度超過閾值的組合事件重現(xiàn)期，F(xiàn)rank Copula函數(shù)對各站的干旱歷時和干旱嚴(yán)重程度擬合效果最好，可以用于該區(qū)域的干旱分析。從空間分布上看，不管是聯(lián)合中度干旱還是重度干旱，湟水流域東南部的干旱風(fēng)險要高于西北部。而流域東南部為需水量大、人口密集的農(nóng)牧業(yè)區(qū)，干旱造成的損失更為嚴(yán)重，因此相關(guān)部門應(yīng)加大該區(qū)域的干旱預(yù)警及監(jiān)測力度，提前做好應(yīng)對措施。

SPI指數(shù)只是基于降水的氣象干旱指標(biāo)，沒有考慮溫度、蒸發(fā)、地形等因素對干旱的影響，與實際情況存在一定的偏差；并且本文在Copula函數(shù)確定干旱評估指標(biāo)上的研究還不夠，有待進一步深入研究。

□

[1] 趙曉慎, 周 海, 王文川. 氣候變化對區(qū)域水循環(huán)系統(tǒng)影響的研究進展[J]. 華北水利水電學(xué)院學(xué)報, 2012,(2):46-49.

[2] 邱玉利, 周建中, 馬 林. 湟水流域地表水資源特征[J]. 水資源與水工程學(xué)報, 2007,(6):98-102.

[3] 郭 武. 青海省湟水流域河川徑流特征分析[J]. 干旱區(qū)研究, 1996,(2):25-30.

[4] 趙串串, 楊曉陽, 張鳳臣, 等. 基于生態(tài)足跡分析的青海湟水河流域可持續(xù)發(fā)展能力[J]. 干旱區(qū)研究, 2009,(3):326-332.

[5] 楊發(fā)源, 戴 升. 湟水河流域水資源變化特征及其對氣候的響應(yīng)[J]. 青海氣象, 2010,(3):12-18.

[6] 韓海濤, 胡文超, 陳學(xué)君, 等. 三種氣象干旱指標(biāo)的應(yīng)用比較研究[J]. 干旱地區(qū)農(nóng)業(yè)研究, 2009,(1):237-241.

[7] 姚 瑤, 張 鑫, 馬 全, 等. 基于SPI青海省東部農(nóng)業(yè)區(qū)季節(jié)干旱變化分析[J]. 灌溉排水學(xué)報, 2013,(6):96-100.

[8] 劉 蕾, 張 鑫, 徐 靜.SPI指數(shù)和Z指數(shù)在青海省東部農(nóng)業(yè)區(qū)應(yīng)用分析[J]. 中國農(nóng)村水利水電, 2015,(4):97-100.

[9] 陸桂華, 閆桂霞, 吳志勇, 等. 基于copula函數(shù)的區(qū)域干旱分析方法[J]. 水科學(xué)進展, 2010,(2):188-193.

[10] 陳永勤, 孫 鵬, 張 強, 等. 基于Copula的鄱陽湖流域水文干旱頻率分析[J]. 自然災(zāi)害學(xué)報, 2013,(1):75-84.

[11] GB/T 20481-2006，氣象干旱等級[S].

[12] 肖名忠, 張 強, 陳曉宏. 基于多變量概率分析的珠江流域干旱特征研究[J]. 地理學(xué)報, 2012,(1):83-92.

[13] Kendall M G. Rank correlation methods[M]. London: Griffin, 1975.

[14] Sklar A. Fonctions de repartition an dimensions et leurs marges[J]. Publ Inst Statist Univ Paris, 1959:229-231.

[15] 孫 鵬, 張 強, 陳曉宏. 基于Copula函數(shù)的鄱陽湖流域極值流量遭遇頻率及災(zāi)害風(fēng)險[J]. 湖泊科學(xué), 2011,(2):183-190.