干旱歷時(shí)擬合方法對(duì)基于Copula函數(shù)的干旱分析影響研究

程雨春，楊傳國(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ì)干旱多變量分析的影響。

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

1.1 降水資料

本文降水?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ū)域的降雨年代際變化特征。

1.2 SPI指數(shù)

本文僅使用降水?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

1.3 干旱事件定義

本文采用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í)間之差。

2 研究方法

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

(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ì)雙變量等值線圖的影響有重要的意義。

2.2 干旱歷時(shí)和強(qiáng)度的擬合

不論單變量還是多變量重現(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)度。

2.3 Copula函數(shù)

計(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 干旱特征分析

3.1 擬合結(jié)果分析

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

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

圖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)。

4 結(jié) 論

本文利用北京和上海近百年逐月降水序列，求出了兩者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)。

□