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

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

0 引 言

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

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

1 材料與方法

1.1 數據來源

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

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

1.2 干旱事件識別

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

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

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

(1)

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

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

1.3 Mann-Kendall 趨勢檢驗法

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

(2)

假設各變量獨立同分布，則統計量S近似服從正態分布[13]。其期望和方差分別如下：

(5)

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

構造M-K統計量Uc：

(5)

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

1.4 Copula函數及其參數估計與選擇

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

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

(6)

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

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

表2 Kendall-τ與copula函數的參數θ對應關系表Tab.2 The corresponding relation between the Kendall-τ and the parameter of Copula Function

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

(7)

式中：n為序列樣本數；pc為Copula函數在原始點的理論概率分布值；pi為聯合經驗概率分布值。

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

(9)

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

2 結果與分析

2.1 單站逐月SPI趨勢分析

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

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

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

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

2.2 干旱歷時和嚴重程度趨勢分析

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

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

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

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

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

2.3 典型干旱情景分析

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

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

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

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

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

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

3 結 語

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

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

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

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

□

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

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

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

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

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

[6] 韓海濤, 胡文超, 陳學君, 等. 三種氣象干旱指標的應用比較研究[J]. 干旱地區農業研究, 2009,(1):237-241.

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

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

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

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

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

[12] 肖名忠, 張 強, 陳曉宏. 基于多變量概率分析的珠江流域干旱特征研究[J]. 地理學報, 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函數的鄱陽湖流域極值流量遭遇頻率及災害風險[J]. 湖泊科學, 2011,(2):183-190.