程雨春,楊傳國(guó),郝振純,楊海燕
(1. 河海大學(xué)水文水資源學(xué)院, 南京 210098;2. 河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 南京 210098)
干旱是一種水分缺乏導(dǎo)致的常見(jiàn)自然現(xiàn)象,有些干旱危害不大,甚至難以察覺(jué),有些干旱則會(huì)嚴(yán)重影響自然環(huán)境和社會(huì)活動(dòng),特別是對(duì)于經(jīng)濟(jì)發(fā)達(dá)、人口密集的地區(qū),干旱更易造成危害。北京和上海是我國(guó)的政治和經(jīng)濟(jì)中心,在全國(guó)大發(fā)展的格局中發(fā)揮著不可替代的作用,研究這兩個(gè)地區(qū)近百年干旱事件發(fā)生規(guī)律,有助于這些地區(qū)預(yù)防和減輕干旱所帶來(lái)的危害。
國(guó)內(nèi)外研究干旱方法很多,如時(shí)間序列分析法,王壬等用線性趨勢(shì)、M-K檢驗(yàn)等分析了雷州半島季節(jié)性氣象干旱時(shí)空特征[1];又如統(tǒng)計(jì)方法,張強(qiáng)等借助Copula函數(shù)分析了東江流域的干旱風(fēng)險(xiǎn)并對(duì)該地區(qū)的水資源管理提出了建議[2],Masud等使用單變量和雙變量概率統(tǒng)計(jì)方法研究了加拿大Saskatchewan河流域的干旱狀況[3]。另外也有許多研究將這兩種方法結(jié)合起來(lái)使用。本文使用Copula函數(shù)等統(tǒng)計(jì)學(xué)方法,以重現(xiàn)期為評(píng)估干旱狀況的主要指標(biāo)。
Copula函數(shù)是研究多元概率分布的重要工具。由于研究干旱事件時(shí)常同時(shí)考慮多個(gè)特征量,如干旱的歷時(shí)、強(qiáng)度、面積等,故Copula函數(shù)在干旱研究中得到了廣泛應(yīng)用。干旱事件各特征量中以干旱歷時(shí)最為特殊,它只能以旬、月等計(jì),是離散變量,而其他均可表示成連續(xù)變量。而正因?yàn)楦珊禋v時(shí)為離散變量,其經(jīng)驗(yàn)點(diǎn)距圖呈現(xiàn)“階梯狀”,這就產(chǎn)生了新的問(wèn)題,即擬合干旱歷時(shí)時(shí)以“階梯中心”還是以“階梯頂部”為標(biāo)準(zhǔn)來(lái)評(píng)估擬合效果好壞的問(wèn)題。因此相對(duì)而言干旱歷時(shí)處理辦法更多,且不同方法會(huì)采用不同擬合優(yōu)度評(píng)價(jià)標(biāo)準(zhǔn)。本文收集了北京和上海近百年降水資料,使用了兩種干旱歷時(shí)擬合方法——參數(shù)和非參數(shù)方法,并比較了不同的擬合優(yōu)度評(píng)價(jià)標(biāo)準(zhǔn),進(jìn)而評(píng)估這些因素對(duì)干旱多變量分析的影響。
本文降水?dāng)?shù)據(jù)來(lái)自《中國(guó)地面氣候資料日值數(shù)據(jù)集》和《中國(guó)降水資料》[4]。將兩者資料整合,得上海徐家匯站1873-2013年和北京站1940-2013年逐月降水序列。通過(guò)1960-2013年54年降雨序列分析,上海站與周邊蘇州站、南通站、平湖站(嘉興市)多年平均降雨量分別為97.6、95.3、91.3、101.9 mm,月降雨序列相關(guān)系數(shù)分別達(dá)到0.86、0.71、0.82。北京站與周邊天津站、滄州站多年平均降雨量分別為49.8、45.1、51.2 mm,月降雨序列相關(guān)系數(shù)分別達(dá)到0.81、0.89。上述代表性分析表明所選站點(diǎn)具有典型的區(qū)域代表性,能夠代表該區(qū)域的降雨年代際變化特征。
本文僅使用降水?dāng)?shù)據(jù),所研究的干旱均指氣象干旱[5]。為研究干旱狀況需定義干旱事件,為此又需干濕程度評(píng)價(jià)指標(biāo)。本文選取SPI指數(shù)評(píng)估干濕狀況,它由McKee提出[6],一直以來(lái)都被廣泛應(yīng)用[7,8],并成為我國(guó)氣象干旱等級(jí)國(guó)家標(biāo)準(zhǔn)推薦的指標(biāo)之一[9]。SPI指數(shù)由降水量的概率分布通過(guò)正態(tài)標(biāo)準(zhǔn)化算出,反映研究區(qū)相對(duì)于平均水平的干濕狀況。圖1給出了北京和上海SPI指數(shù)的霍弗莫勒?qǐng)D,其中黃色表示偏干,藍(lán)色表示偏濕。由圖可見(jiàn)兩地干旱事件發(fā)生頻繁,并存在顯著的年際和年代際變化。北京站在20世紀(jì)60年代夏秋旱頻繁,2000年前后發(fā)生了嚴(yán)重的秋旱災(zāi)害;上海站在1930年前后發(fā)生了嚴(yán)重的春旱和夏旱,1960年代的秋旱也很頻繁。近年來(lái)兩地處于偏濕狀態(tài),干旱災(zāi)害不顯著。
圖1 北京和上海兩站SPI指數(shù)變化分布圖Fig.1 Variation of SPI of Beijing and Shanghai
本文采用McKee的定義[10]:一場(chǎng)干旱事件從SPI的值小于0開(kāi)始,到此后SPI第一次大于0結(jié)束,且其中至少一個(gè)SPI值小于等于-1。定義了干旱事件后,即可定義干旱歷時(shí)D、干旱強(qiáng)度S和干旱間隔時(shí)間L。干旱歷時(shí)即事件持續(xù)的時(shí)間;干旱強(qiáng)度即事件中各時(shí)刻SPI指數(shù)之和的絕對(duì)值;間隔時(shí)間即為兩次干旱事件開(kāi)始時(shí)間之差。
(1)
(2)
(3)
(4)
式中:E(L)為間隔時(shí)間的數(shù)學(xué)期望;FS(x)和FD(x)分別為干旱強(qiáng)度和歷時(shí)的分布函數(shù);FDS(x,y)是歷時(shí)和強(qiáng)度的聯(lián)合分布函數(shù);Ts和Td分別為強(qiáng)度S≥s的干旱和歷時(shí)D≥d的干旱的單變量重現(xiàn)期,可如下計(jì)算:
(5)
(6)
給出一組歷時(shí)和強(qiáng)度(d,s),即可計(jì)算該組干旱特征值對(duì)應(yīng)的雙變量重現(xiàn)期,但不同干旱特征值組可能對(duì)應(yīng)相同雙變量重現(xiàn)期,只有作雙變量重現(xiàn)期等值線圖,才能全面掌握干旱狀況。因此,研究邊際分布擬合方法對(duì)雙變量等值線圖的影響有重要的意義。
不論單變量還是多變量重現(xiàn)期,都需先求出干旱歷時(shí)和強(qiáng)度的分布FD(x)和FS(x)。干旱歷時(shí)可用參數(shù)和非參數(shù)方法處理。參數(shù)方法指?jìng)鹘y(tǒng)統(tǒng)計(jì)學(xué)中的線型擬合,本文取指數(shù)分布為代表。非參數(shù)方法包括核密度估計(jì)[13]、正交級(jí)數(shù)估計(jì)[14]、馬秀峰游程理論等[15]。本文取游程理論為代表,此方法中歷時(shí)的累積概率分布可如下計(jì)算:
(7)
其中:
(8)
式中:g(d)是干旱歷時(shí)不小于d的頻次;n為實(shí)測(cè)的樣本長(zhǎng)度。
計(jì)算雙變量重現(xiàn)期需聯(lián)合分布函數(shù)FDS(x,y),它可用Copula函數(shù)計(jì)算。Copula函數(shù)C:[0,1]2→[0,1]定義如下:若兩隨機(jī)變量的分布函數(shù)為F(x)和G(y),則它們的聯(lián)合分布函數(shù)是:
H(x,y)=C[F(x),G(y)],x,y∈R
(9)
諸Copula函數(shù)中,以對(duì)稱Archimedean Copula函數(shù)類在水文分析中最為常用[16],對(duì)干旱洪澇等極端事件的擬合效果較好。對(duì)稱Archimedean Copula定義為:
C(u,v)=γ[-1][γ(u)+γ(v) ]u,v∈[0,1]
(10)
其中函數(shù)γ(·)稱為生成函數(shù),γ[-1]表示γ(·)的廣義逆函數(shù)。生成函數(shù)γ(·)有一個(gè)參數(shù)θ,同時(shí)也是Copula函數(shù)的參數(shù)。參數(shù)估計(jì)時(shí)常用Kendall相關(guān)系數(shù)τ,因?yàn)槔碚撋峡赏茖?dǎo)出τ~θ關(guān)系。本文先用樣本Kendall相關(guān)系數(shù)τn估計(jì)τ:
(11)
再用τ~θ關(guān)系估計(jì)θ??梢宰C明τn是τ的漸近無(wú)偏估計(jì)量[17]。Gumbel、Clayton和Frank Copula的基本信息列于表1。Copula函數(shù)的擬合優(yōu)度檢驗(yàn)基于經(jīng)驗(yàn)Copula和Cramer-von Mises統(tǒng)計(jì)量, 若多個(gè)Copula函數(shù)都能通過(guò),則要用留一法交叉檢驗(yàn)(leave-one-out cross validation)選出最優(yōu)者。
表1 對(duì)稱Archimedean Copula函數(shù)表Tab.1 Symmetric Archimedean Copulas used in this paper
3.1.1 歷時(shí)擬合結(jié)果
在上海和北京兩站分別用兩種方法擬合干旱歷時(shí),得如圖2~3的結(jié)果??梢钥闯?,指數(shù)分布的擬合結(jié)果穿過(guò)經(jīng)驗(yàn)點(diǎn)距群的中心,而游程理論的擬合結(jié)果接近經(jīng)驗(yàn)點(diǎn)距圖像中各階梯的頂部,梁忠民等稱之為經(jīng)驗(yàn)累積概率點(diǎn)距[18]。干旱分析中,以經(jīng)驗(yàn)點(diǎn)距群的中心還是經(jīng)驗(yàn)累積概率點(diǎn)距為標(biāo)準(zhǔn)擬合干旱歷時(shí)仍有爭(zhēng)論,本文將比較這兩者在干旱風(fēng)險(xiǎn)評(píng)估中的表現(xiàn)。
圖2 上海干旱歷時(shí)擬合結(jié)果Fig.2 Results of fitting drought duration of Shanghai
圖3 北京干旱歷時(shí)擬合結(jié)果Fig.3 Results of fitting drought duration of Beijing
用參數(shù)方法可以計(jì)算任意大的干旱歷時(shí)對(duì)應(yīng)的概率,而游程方法則不能隨意地超出實(shí)測(cè)數(shù)據(jù),例如,如果干旱歷時(shí)大于樣本長(zhǎng)度,則用式(7)~(8)算出的概率將可能大于1。因此,本文中凡使用游程方法的地方都進(jìn)行了插值和外延,并嚴(yán)格控制了外延的范圍。其他非參數(shù)擬合方法也會(huì)遇到類似的問(wèn)題,特別是對(duì)特大的干旱歷時(shí)它們都需進(jìn)行外延,所以非參數(shù)方法在處理極端值時(shí)不如參數(shù)方法,即傳統(tǒng)線型擬合方法有優(yōu)勢(shì)。
3.1.2 強(qiáng)度擬合結(jié)果
強(qiáng)度為連續(xù)變量,一般用參數(shù)方法即線型擬合即可。試驗(yàn)表明,用P-III分布擬合強(qiáng)度就能得到很好的效果,見(jiàn)圖4-5。
圖4 上海干旱強(qiáng)度擬合結(jié)果Fig.4 Results of fitting drought severity of Shanghai
圖5 北京干旱強(qiáng)度擬合結(jié)果Fig.5 Results of fitting drought severity of Beijing
3.1.3Copula函數(shù)擬合結(jié)果
Copula函數(shù)的擬合優(yōu)度檢驗(yàn)和交叉檢驗(yàn)法結(jié)果如表2。在95%置信度水平下,P值達(dá)到5%即可認(rèn)為通過(guò),而交叉檢驗(yàn)值越高則表示擬合效果越好??梢钥闯?,北京和上海兩站均以GumbelCopula為最優(yōu)。
表2 Copula函數(shù)擬合效果比較表Tab.2 Comparison of goodness-of-fit of the three Copulas
圖6 上海和北京雙變量重現(xiàn)期等值線圖Fig.6 Isolines of bivariate return periods of Shanghai and Beijing
由上海、北京的“且”重現(xiàn)期和“或”重現(xiàn)期等值線圖可以看出,游程方法得到的等值線圖與指數(shù)分布得到的等值線圖中的每根等值線形狀上相差一個(gè)沿橫軸方向的拉伸??梢?jiàn)干旱歷時(shí)擬合方法對(duì)Copula函數(shù)等值線圖的影響作用在歷時(shí)坐標(biāo)軸方向上,對(duì)強(qiáng)度方向沒(méi)有影響。由對(duì)稱性[見(jiàn)式(1)、(2)]可以想到,若干旱強(qiáng)度擬合方法不同,那么對(duì)應(yīng)等值線的差異體現(xiàn)在強(qiáng)度方向上,對(duì)歷時(shí)方向沒(méi)有影響。因此,對(duì)于“且”重現(xiàn)期和“或”重現(xiàn)期等值線圖而言,歷時(shí)擬合對(duì)其造成的影響與強(qiáng)度擬合造成的影響是相互獨(dú)立的。
由圖6(a)~(b)看出,對(duì)所有的實(shí)測(cè)干旱事件,游程方法計(jì)算出的“且”重現(xiàn)期都不小于指數(shù)分布計(jì)算出的“且”重現(xiàn)期;同樣,由圖6(c)~(d),游程方法算出的“或”重現(xiàn)期都也不小于指數(shù)分布算出的“或”重現(xiàn)期,由此可見(jiàn),不論是認(rèn)為干旱的兩個(gè)方面同樣重要,任一方面都不能超過(guò)危險(xiǎn)閾值,還是只關(guān)心歷時(shí)和強(qiáng)度都超過(guò)危險(xiǎn)值的大旱事件,以累積經(jīng)驗(yàn)概率點(diǎn)距為標(biāo)準(zhǔn)擬合歷時(shí)評(píng)估得到的干旱風(fēng)險(xiǎn)總是低于以經(jīng)驗(yàn)點(diǎn)距群的中心為標(biāo)準(zhǔn)擬合干旱歷時(shí)所得的評(píng)估結(jié)果。
但條件重現(xiàn)期等值線圖6(e)~(f)的另一特征與“且”和“或”重現(xiàn)期等值線圖相同,即游程方法擬合(以累積經(jīng)驗(yàn)概率點(diǎn)距為標(biāo)準(zhǔn)擬合)歷時(shí)得到的等值線總是被指數(shù)分布擬合(以經(jīng)驗(yàn)點(diǎn)距群的中心為標(biāo)準(zhǔn)擬合)歷時(shí)得到的相應(yīng)等值線包圍,故對(duì)于相同的干旱狀況,前者總是會(huì)有更大的重現(xiàn)期,從而評(píng)估得到更低的風(fēng)險(xiǎn)。
本文利用北京和上海近百年逐月降水序列,求出了兩者SPI指數(shù)序列,識(shí)別了干旱事件及它們的歷時(shí)和強(qiáng)度,計(jì)算了雙變量重現(xiàn)期,比較了參數(shù)方法和非參數(shù)方法進(jìn)行干旱歷時(shí)擬合對(duì)雙變量重現(xiàn)期計(jì)算的影響。主要結(jié)論有:
(1)指數(shù)分布方法擬合歷時(shí)以經(jīng)驗(yàn)點(diǎn)距群中心為擬合優(yōu)度評(píng)估標(biāo)準(zhǔn),而馬秀峰游程理論以經(jīng)驗(yàn)累積概率點(diǎn)距為擬合優(yōu)度評(píng)估標(biāo)準(zhǔn),兩者計(jì)算結(jié)果均有其合理性。但非參數(shù)擬合方法在處理特大干旱歷時(shí)時(shí)都需要進(jìn)行外延,而參數(shù)方法不需要,因而在這方面優(yōu)于非參數(shù)方法。
(2)不論以“且”重現(xiàn)期、“或”重現(xiàn)期和條件重現(xiàn)期中的哪一種作為設(shè)計(jì)標(biāo)準(zhǔn),在干旱歷時(shí)和強(qiáng)度都不是極端大的情況下,游程理論擬合歷時(shí)評(píng)估出的風(fēng)險(xiǎn)總小于指數(shù)分布擬合歷時(shí)評(píng)估出的風(fēng)險(xiǎn)。
□