程雨春，楊傳國，郝振純，楊海燕

(1. 河海大學(xué)水文水資源學(xué)院, 南京 210098;2. 河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室, 南京 210098)

干旱是一種水分缺乏導(dǎo)致的常見自然現(xiàn)象，有些干旱危害不大，甚至難以察覺，有些干旱則會嚴(yán)重影響自然環(huán)境和社會活動，特別是對于經(jīng)濟發(fā)達、人口密集的地區(qū)，干旱更易造成危害。北京和上海是我國的政治和經(jīng)濟中心，在全國大發(fā)展的格局中發(fā)揮著不可替代的作用，研究這兩個地區(qū)近百年干旱事件發(fā)生規(guī)律，有助于這些地區(qū)預(yù)防和減輕干旱所帶來的危害。

國內(nèi)外研究干旱方法很多，如時間序列分析法，王壬等用線性趨勢、M-K檢驗等分析了雷州半島季節(jié)性氣象干旱時空特征[1]；又如統(tǒng)計方法，張強等借助Copula函數(shù)分析了東江流域的干旱風(fēng)險并對該地區(qū)的水資源管理提出了建議[2]，Masud等使用單變量和雙變量概率統(tǒng)計方法研究了加拿大Saskatchewan河流域的干旱狀況[3]。另外也有許多研究將這兩種方法結(jié)合起來使用。本文使用Copula函數(shù)等統(tǒng)計學(xué)方法，以重現(xiàn)期為評估干旱狀況的主要指標(biāo)。

Copula函數(shù)是研究多元概率分布的重要工具。由于研究干旱事件時常同時考慮多個特征量，如干旱的歷時、強度、面積等，故Copula函數(shù)在干旱研究中得到了廣泛應(yīng)用。干旱事件各特征量中以干旱歷時最為特殊，它只能以旬、月等計，是離散變量，而其他均可表示成連續(xù)變量。而正因為干旱歷時為離散變量，其經(jīng)驗點距圖呈現(xiàn)“階梯狀”，這就產(chǎn)生了新的問題，即擬合干旱歷時時以“階梯中心”還是以“階梯頂部”為標(biāo)準(zhǔn)來評估擬合效果好壞的問題。因此相對而言干旱歷時處理辦法更多，且不同方法會采用不同擬合優(yōu)度評價標(biāo)準(zhǔn)。本文收集了北京和上海近百年降水資料，使用了兩種干旱歷時擬合方法——參數(shù)和非參數(shù)方法，并比較了不同的擬合優(yōu)度評價標(biāo)準(zhǔn)，進而評估這些因素對干旱多變量分析的影響。

1 降水資料與干旱事件定義

1.1 降水資料

本文降水?dāng)?shù)據(jù)來自《中國地面氣候資料日值數(shù)據(jù)集》和《中國降水資料》[4]。將兩者資料整合，得上海徐家匯站1873-2013年和北京站1940-2013年逐月降水序列。通過1960-2013年54年降雨序列分析，上海站與周邊蘇州站、南通站、平湖站(嘉興市)多年平均降雨量分別為97.6、95.3、91.3、101.9 mm，月降雨序列相關(guān)系數(shù)分別達到0.86、0.71、0.82。北京站與周邊天津站、滄州站多年平均降雨量分別為49.8、45.1、51.2 mm，月降雨序列相關(guān)系數(shù)分別達到0.81、0.89。上述代表性分析表明所選站點具有典型的區(qū)域代表性，能夠代表該區(qū)域的降雨年代際變化特征。

1.2 SPI指數(shù)

本文僅使用降水?dāng)?shù)據(jù)，所研究的干旱均指氣象干旱[5]。為研究干旱狀況需定義干旱事件，為此又需干濕程度評價指標(biāo)。本文選取SPI指數(shù)評估干濕狀況，它由McKee提出[6]，一直以來都被廣泛應(yīng)用[7,8]，并成為我國氣象干旱等級國家標(biāo)準(zhǔn)推薦的指標(biāo)之一[9]。SPI指數(shù)由降水量的概率分布通過正態(tài)標(biāo)準(zhǔn)化算出，反映研究區(qū)相對于平均水平的干濕狀況。圖1給出了北京和上海SPI指數(shù)的霍弗莫勒圖，其中黃色表示偏干，藍色表示偏濕。由圖可見兩地干旱事件發(fā)生頻繁，并存在顯著的年際和年代際變化。北京站在20世紀(jì)60年代夏秋旱頻繁，2000年前后發(fā)生了嚴(yán)重的秋旱災(zāi)害；上海站在1930年前后發(fā)生了嚴(yán)重的春旱和夏旱，1960年代的秋旱也很頻繁。近年來兩地處于偏濕狀態(tài)，干旱災(zāi)害不顯著。

圖1 北京和上海兩站SPI指數(shù)變化分布圖Fig.1 Variation of SPI of Beijing and Shanghai

1.3 干旱事件定義

本文采用McKee的定義[10]：一場干旱事件從SPI的值小于0開始，到此后SPI第一次大于0結(jié)束，且其中至少一個SPI值小于等于-1。定義了干旱事件后，即可定義干旱歷時D、干旱強度S和干旱間隔時間L。干旱歷時即事件持續(xù)的時間；干旱強度即事件中各時刻SPI指數(shù)之和的絕對值；間隔時間即為兩次干旱事件開始時間之差。

2 研究方法

2.1 雙變量重現(xiàn)期

(1)

(2)

(3)

(4)

式中：E(L)為間隔時間的數(shù)學(xué)期望；FS(x)和FD(x)分別為干旱強度和歷時的分布函數(shù)；FDS(x,y)是歷時和強度的聯(lián)合分布函數(shù)；Ts和Td分別為強度S≥s的干旱和歷時D≥d的干旱的單變量重現(xiàn)期，可如下計算：

(5)

(6)

給出一組歷時和強度(d,s)，即可計算該組干旱特征值對應(yīng)的雙變量重現(xiàn)期，但不同干旱特征值組可能對應(yīng)相同雙變量重現(xiàn)期，只有作雙變量重現(xiàn)期等值線圖，才能全面掌握干旱狀況。因此，研究邊際分布擬合方法對雙變量等值線圖的影響有重要的意義。

2.2 干旱歷時和強度的擬合

不論單變量還是多變量重現(xiàn)期，都需先求出干旱歷時和強度的分布FD(x)和FS(x)。干旱歷時可用參數(shù)和非參數(shù)方法處理。參數(shù)方法指傳統(tǒng)統(tǒng)計學(xué)中的線型擬合，本文取指數(shù)分布為代表。非參數(shù)方法包括核密度估計[13]、正交級數(shù)估計[14]、馬秀峰游程理論等[15]。本文取游程理論為代表，此方法中歷時的累積概率分布可如下計算：

(7)

其中：

(8)

式中：g(d)是干旱歷時不小于d的頻次；n為實測的樣本長度。

2.3 Copula函數(shù)

計算雙變量重現(xiàn)期需聯(lián)合分布函數(shù)FDS(x,y)，它可用Copula函數(shù)計算。Copula函數(shù)C:[0,1]2→[0,1]定義如下：若兩隨機變量的分布函數(shù)為F(x)和G(y)，則它們的聯(lián)合分布函數(shù)是：

H(x,y)=C[F(x),G(y)],x,y∈R

(9)

諸Copula函數(shù)中，以對稱Archimedean Copula函數(shù)類在水文分析中最為常用[16]，對干旱洪澇等極端事件的擬合效果較好。對稱Archimedean Copula定義為：

C(u,v)=γ[-1][γ(u)+γ(v) ]u,v∈[0,1]

(10)

其中函數(shù)γ(·)稱為生成函數(shù)，γ[-1]表示γ(·)的廣義逆函數(shù)。生成函數(shù)γ(·)有一個參數(shù)θ，同時也是Copula函數(shù)的參數(shù)。參數(shù)估計時常用Kendall相關(guān)系數(shù)τ，因為理論上可推導(dǎo)出τ～θ關(guān)系。本文先用樣本Kendall相關(guān)系數(shù)τn估計τ：

(11)

再用τ～θ關(guān)系估計θ。可以證明τn是τ的漸近無偏估計量[17]。Gumbel、Clayton和Frank Copula的基本信息列于表1。Copula函數(shù)的擬合優(yōu)度檢驗基于經(jīng)驗Copula和Cramer-von Mises統(tǒng)計量， 若多個Copula函數(shù)都能通過，則要用留一法交叉檢驗(leave-one-out cross validation)選出最優(yōu)者。

表1 對稱Archimedean Copula函數(shù)表Tab.1 Symmetric Archimedean Copulas used in this paper

3 干旱特征分析

3.1 擬合結(jié)果分析

3.1.1 歷時擬合結(jié)果

在上海和北京兩站分別用兩種方法擬合干旱歷時，得如圖2～3的結(jié)果。可以看出，指數(shù)分布的擬合結(jié)果穿過經(jīng)驗點距群的中心，而游程理論的擬合結(jié)果接近經(jīng)驗點距圖像中各階梯的頂部，梁忠民等稱之為經(jīng)驗累積概率點距[18]。干旱分析中，以經(jīng)驗點距群的中心還是經(jīng)驗累積概率點距為標(biāo)準(zhǔn)擬合干旱歷時仍有爭論，本文將比較這兩者在干旱風(fēng)險評估中的表現(xiàn)。

圖2 上海干旱歷時擬合結(jié)果Fig.2 Results of fitting drought duration of Shanghai

圖3 北京干旱歷時擬合結(jié)果Fig.3 Results of fitting drought duration of Beijing

用參數(shù)方法可以計算任意大的干旱歷時對應(yīng)的概率，而游程方法則不能隨意地超出實測數(shù)據(jù)，例如，如果干旱歷時大于樣本長度，則用式(7)～(8)算出的概率將可能大于1。因此，本文中凡使用游程方法的地方都進行了插值和外延，并嚴(yán)格控制了外延的范圍。其他非參數(shù)擬合方法也會遇到類似的問題，特別是對特大的干旱歷時它們都需進行外延，所以非參數(shù)方法在處理極端值時不如參數(shù)方法，即傳統(tǒng)線型擬合方法有優(yōu)勢。

3.1.2 強度擬合結(jié)果

強度為連續(xù)變量，一般用參數(shù)方法即線型擬合即可。試驗表明，用P-III分布擬合強度就能得到很好的效果，見圖4-5。

圖4 上海干旱強度擬合結(jié)果Fig.4 Results of fitting drought severity of Shanghai

圖5 北京干旱強度擬合結(jié)果Fig.5 Results of fitting drought severity of Beijing

3.1.3Copula函數(shù)擬合結(jié)果

Copula函數(shù)的擬合優(yōu)度檢驗和交叉檢驗法結(jié)果如表2。在95%置信度水平下，P值達到5%即可認(rèn)為通過，而交叉檢驗值越高則表示擬合效果越好?？梢钥闯觯本┖蜕虾烧揪訥umbelCopula為最優(yōu)。

表2 Copula函數(shù)擬合效果比較表Tab.2 Comparison of goodness-of-fit of the three Copulas

3.2 雙變量重現(xiàn)期結(jié)果分析

圖6 上海和北京雙變量重現(xiàn)期等值線圖Fig.6 Isolines of bivariate return periods of Shanghai and Beijing

由上海、北京的“且”重現(xiàn)期和“或”重現(xiàn)期等值線圖可以看出，游程方法得到的等值線圖與指數(shù)分布得到的等值線圖中的每根等值線形狀上相差一個沿橫軸方向的拉伸?？梢姼珊禋v時擬合方法對Copula函數(shù)等值線圖的影響作用在歷時坐標(biāo)軸方向上，對強度方向沒有影響。由對稱性[見式(1)、(2)]可以想到，若干旱強度擬合方法不同，那么對應(yīng)等值線的差異體現(xiàn)在強度方向上，對歷時方向沒有影響。因此，對于“且”重現(xiàn)期和“或”重現(xiàn)期等值線圖而言，歷時擬合對其造成的影響與強度擬合造成的影響是相互獨立的。

由圖6(a)～(b)看出，對所有的實測干旱事件，游程方法計算出的“且”重現(xiàn)期都不小于指數(shù)分布計算出的“且”重現(xiàn)期；同樣，由圖6(c)～(d)，游程方法算出的“或”重現(xiàn)期都也不小于指數(shù)分布算出的“或”重現(xiàn)期，由此可見，不論是認(rèn)為干旱的兩個方面同樣重要，任一方面都不能超過危險閾值，還是只關(guān)心歷時和強度都超過危險值的大旱事件，以累積經(jīng)驗概率點距為標(biāo)準(zhǔn)擬合歷時評估得到的干旱風(fēng)險總是低于以經(jīng)驗點距群的中心為標(biāo)準(zhǔn)擬合干旱歷時所得的評估結(jié)果。

但條件重現(xiàn)期等值線圖6(e)～(f)的另一特征與“且”和“或”重現(xiàn)期等值線圖相同，即游程方法擬合(以累積經(jīng)驗概率點距為標(biāo)準(zhǔn)擬合)歷時得到的等值線總是被指數(shù)分布擬合(以經(jīng)驗點距群的中心為標(biāo)準(zhǔn)擬合)歷時得到的相應(yīng)等值線包圍，故對于相同的干旱狀況，前者總是會有更大的重現(xiàn)期，從而評估得到更低的風(fēng)險。

4 結(jié) 論

本文利用北京和上海近百年逐月降水序列，求出了兩者SPI指數(shù)序列，識別了干旱事件及它們的歷時和強度，計算了雙變量重現(xiàn)期，比較了參數(shù)方法和非參數(shù)方法進行干旱歷時擬合對雙變量重現(xiàn)期計算的影響。主要結(jié)論有：

(1)指數(shù)分布方法擬合歷時以經(jīng)驗點距群中心為擬合優(yōu)度評估標(biāo)準(zhǔn)，而馬秀峰游程理論以經(jīng)驗累積概率點距為擬合優(yōu)度評估標(biāo)準(zhǔn)，兩者計算結(jié)果均有其合理性。但非參數(shù)擬合方法在處理特大干旱歷時時都需要進行外延，而參數(shù)方法不需要，因而在這方面優(yōu)于非參數(shù)方法。

(2)不論以“且”重現(xiàn)期、“或”重現(xiàn)期和條件重現(xiàn)期中的哪一種作為設(shè)計標(biāo)準(zhǔn)，在干旱歷時和強度都不是極端大的情況下，游程理論擬合歷時評估出的風(fēng)險總小于指數(shù)分布擬合歷時評估出的風(fēng)險。

□