基于Copula方法的干旱歷時和烈度的聯(lián)合概率分析

許月萍，張慶慶，樓章華，劉德地

(浙江大學(xué)建筑工程學(xué)院水文與水資源研究所，杭州 310058)

水文干旱極限分析一直是水文學(xué)家及水資源決策者共同關(guān)心的問題之一，近幾年來，隨著全球氣候變化，干旱事件發(fā)生的頻率和強度增加，研究極限干旱事件越發(fā)重要．干旱事件通常有 3個重要指標，即干旱歷時、烈度和間隔時間．然而，分析干旱事件往往只關(guān)注干旱的歷時或干旱烈度，如馮國章[1]應(yīng)用獨立同分布隨機變量序列最大負游程長的概率分布函數(shù)進行極限干旱歷時最大干旱持續(xù)時間的概率分析；丁晶等[2]分析了作為水文干旱定量指標負輪長的統(tǒng)計變化特性，并以哈爾濱和陜縣站為例，應(yīng)用隨機模擬方法探討了嚴重干旱出現(xiàn)可能性的定量估計方法；朱廷舉等[3]探討了嚴重水文干旱事件的發(fā)生規(guī)律和確定黃河上中游地區(qū)連續(xù)枯水段的發(fā)生概率和重現(xiàn)期；袁超等[4]研究了 P-Ⅲ分布和馬爾可夫過程年徑流序列極限水文干旱歷時概率分布隨水文參數(shù)的變化規(guī)律．對于干旱事件，干旱歷時、干旱烈度甚至間隔時間都非常重要，目前針對這三者的多元變量分析在國內(nèi)外已經(jīng)引起了充分重視，如González等[5]應(yīng)用經(jīng)典多元變量分析方法建立了邊際分布分別為幾何和伽馬分布的干旱歷時和干旱烈度兩元概率分布；Shiau等[6]則應(yīng)用 Copula方法建立了邊際分布分別為混合指數(shù)和伽馬分布的歷時和強度兩元概率分布．鑒于干旱數(shù)據(jù)短缺，筆者利用自回歸馬爾可夫模型對水文序列進行延長，在此基礎(chǔ)上，應(yīng)用Copula方法對干旱發(fā)生的歷時和強度的概率分布進行定量研究，并以浙江省錢塘江流域諸暨水文站的年降雨序列為例，分析了此方法的實際可行性．

1 研究方法

因干旱時間序列短，直接對短序列進行分析往往不確定性比較大．采用兩步驟來對干旱進行分析．第1步是通過自回歸馬爾可夫模型對水文序列進行延長，以解決短序列不確定性大的問題．第 2步是利用Copula方法對多變量水文數(shù)據(jù)進行分析，考慮多變量之間的相依性．

1.1 干旱的定義

對水文干旱進行分析，一般采用 Yevjevich[7]的游程理論．圖 1為干旱定義的示意圖，水文極限干旱事件一般定義為所關(guān)注的水文氣象變量處于某一截取水平R0以下的時期．此截取水平可為常數(shù)，也可隨時間變化．可以看出干旱事件 3個重要指標即干旱歷時、烈度和間隔時間的定義．

圖1 干旱事件的定義Fig.1 Definition of drought event

1.2 自回歸馬爾可夫模型

研究極限事件，最主要的困難是資料不足，水文干旱事件同樣如此．為揭示水文干旱的基本規(guī)律，利用現(xiàn)有觀測序列所包含的統(tǒng)計信息進行隨機模擬．假定水文序列是一個馬爾可夫過程，并服從正態(tài)分布．利用下面的模型就可以生成均值為μ、方差為2σ、逐年相關(guān)系數(shù)為ρ的正態(tài)分布的徑流量序列．

式中：y1v+為獨立標準正態(tài)隨機變量；ρ為相關(guān)系數(shù)；σ為標準差；μ為均值；1yQ+為 1y+時刻的徑流量．

然而水文序列有時候并不適宜用正態(tài)分布來描述，調(diào)整式(1)中隨機成分的偏差以使模型生成的水文序列具有相應(yīng)的均值、方差和偏差系數(shù)，也即保持序列的部分統(tǒng)計特性[8]．這里利用標準化的χ2分布隨機數(shù)ξ來實現(xiàn)．當χ2分布的自由度大于 30時，ξ與標準正態(tài)的隨機數(shù)之間的關(guān)系為

式中ξγ為ξ的偏差系數(shù)，與水文序列偏差系數(shù)Qγ的關(guān)系為

因此，偏態(tài)水文序列的自回歸馬爾科夫模型為

1.3 Copula方法

1.3.1 Copula方法簡述

Copula作為一種全新的相關(guān)性度量工具，由Sklar在1959第一次提出．Sklar定理是Copula領(lǐng)域中最重要的一個定理，它證明了 Copula存在的唯一性．假設(shè)X和 Y是兩個隨機變量，它們的聯(lián)合分布函數(shù)為G(x，y)，邊緣分布函數(shù)分別為Fx(x)和Fy(y)，則一定存在著一個聯(lián)結(jié)函數(shù) Copula C，滿足 G(x，y)＝C(Fx(x)，F(xiàn)y(y))．如果 Fx(x)和 Fy(y)是連續(xù)的，那么C 是唯一的．相反，如果 C 是一個 Copula，F(xiàn)x(x)和Fy(y)是分布函數(shù)，則由 G(x，y)＝C(Fx(x)，F(xiàn)y(y))所定義的 G(x，y)是一個聯(lián)合分布函數(shù)，而其邊緣分布函數(shù)分別是Fx(x)和Fy(y)．

可用來描述極限情況下的多變量相依關(guān)系的Copula 有 Archimedean Copula、Farlie-Gumbel-Morgenstern Copula、Joe BB1 Copula和 Joe BB5 Copula等．Copula方法因數(shù)學(xué)表達上比較簡潔，各個變量之間的相依結(jié)構(gòu)和邊緣分布可以獨立考慮，所以其在國內(nèi)外水文學(xué)中的應(yīng)用日漸增加[9-12]．作為水文干旱極限分析探索性的研究，筆者應(yīng)用 Archimedean聯(lián)結(jié)函數(shù)中的 Clayton Copula．它可以用來描述變量之間的正相依性，而且適合用來模擬尾部，也即極限事件．以雙變量為例，Clayton Copula的概率分布為

而經(jīng)驗Copula定義為

式中：n為序列長度；l為指示函數(shù)．

1.3.2 Copula擬合檢驗

判斷Clayton Copula方法能否很好地描述相依變量間的概率分布，通常是通過比較所用 Copula跟經(jīng)驗 Copula的距離來完成的，用 Kolmogorov-Smirnov(KS)統(tǒng)計量來做擬合檢驗．KS統(tǒng)計量計算公式為

利用自助抽樣 Bootstrap法來計算 KS統(tǒng)計量，具體步驟為：

根據(jù)觀測樣本計算Clayton Copula的參數(shù)θ；

根據(jù)式(5)計算經(jīng)驗CopulaCn的值；

根據(jù)式(6)計算KS統(tǒng)計量 Sn的值；

選用一個較大的整數(shù) N，對每一個 k ∈ { 1,… ,N }重復(fù)以下步驟，即

(1)從Cθ里面產(chǎn)生隨機數(shù)，并計算相應(yīng)的，這里為生成序列的秩；

如果通過自助抽樣法計算出來的P值太大，則假設(shè)不成立，即所用的 Copula模型不能很好地模擬多變量之間的相關(guān)關(guān)系．P值是進行檢驗決策的依據(jù)，即概率反映某一事件發(fā)生的可能性大?。?/p>

1.3.3 聯(lián)合概率分布計算

利用 Copula函數(shù)，計算X和Y聯(lián)合概率分布為

2 應(yīng) 用

采用浙江省錢塘江流域上諸暨站的年降雨量數(shù)據(jù)，從 1951—2003年一共 53年．多年來浙江省資源型、水質(zhì)型和工程型缺水并存，據(jù)預(yù)測，浙江省在現(xiàn)狀情況下，如遭遇一般干旱年，將缺水 6億 m3，如遭遇大旱年份，將缺水23億m3．到 2010年，隨著工業(yè)的發(fā)展、城市化率和人民生活水平的提高，如不增加水源建設(shè)，遭遇一般干旱年，將缺水71億m3，如遭遇大旱年份，將缺水90億m3．因此，對水文干旱情況進行分析至關(guān)重要，可為浙江省水資源的有效管理提供技術(shù)支持．

2.1 干旱歷時和烈度數(shù)據(jù)獲取

圖2所示為錢塘江諸暨站的年降水量序列，其多年平均降雨量約為1,400,mm．從圖2中可以看出，諸暨站枯水和豐水常成組出現(xiàn)，這說明年降雨量時間序列存在著正相依性．假如以多年平均降雨量為截取水平，可以看到諸暨站很多年份都處于缺水狀態(tài)．

圖2 諸暨站的年降雨量時間序列Fig.2 Annual precipitation time series at Zhuji station

圖3所示為諸暨站利用自回歸馬爾可夫模型生成的1,000年年降雨數(shù)據(jù)(截取水平為1,400,mm)，表1為模擬數(shù)據(jù)跟實際數(shù)據(jù)主要統(tǒng)計特征值．可以看出，自回歸馬爾可夫模型較合理地保持了實測年降雨序列的基本統(tǒng)計信息值．

圖3 自回歸馬爾可夫模型模擬結(jié)果Fig.3 Simulation results of autoregressive Markov model

表1 模擬結(jié)果與實測數(shù)據(jù)基本統(tǒng)計特性比較Tab.1 Comparison of simulation and actual data statistics

選取多年平均降雨量作為截取水平，由圖1可獲得干旱歷時、干旱烈度和間隔時間數(shù)據(jù)，系列長度為243年．表 2為干旱歷時、烈度和間隔時間的一些基本統(tǒng)計信息，可看出這些變量并不遵循正態(tài)分布．

表2 干旱歷時、烈度和間隔時間的一些基本統(tǒng)計特性Tab.2 Basic statistics of drought duration，severity and Tab.2 interval time

選取干旱歷時和烈度兩個變量來做干旱分析，圖4所示為兩變量之間的相關(guān)關(guān)系，可以看出兩者呈現(xiàn)較大的正相關(guān)關(guān)系．經(jīng)計算，其Kendall秩相關(guān)系數(shù)τ為 0.70．

圖4 干旱歷時和烈度相關(guān)關(guān)系Fig.4 Relationship of drought duration and severity

2.2 擬合檢驗

采用Joe[13]提出的IFM方法，先估計邊際函數(shù)及其參數(shù)，然后估計Copula的參數(shù)．應(yīng)用線性距比例圖得出，干旱歷時符合皮爾遜Ⅲ型曲線，干旱烈度符合伽馬分布函數(shù)．皮爾遜Ⅲ型和伽馬分布函數(shù)的公式分別為和()f x=其中Γ為伽瑪函數(shù)，α、β、a0均為皮爾遜Ⅲ型曲線的3個參數(shù)，α1、β1分別為伽馬分布函數(shù)的兩個參數(shù)．用 Kolmogorov-Smirnov檢驗兩者的適用性，計算得出P分別為0.01和0.04．取置信水平 5%，則兩者都小于臨界值，因此假設(shè)成立．利用線性距方法估計得出：α =1.52,β =0.77,a0=0.13(皮爾遜Ⅲ型曲線的3個參數(shù))，α1= 3 54.04，β1= 1 .15(伽馬分布函數(shù)的兩個參數(shù))．

利用自助抽樣法計算Kolmogorov-Smirnov統(tǒng)計量來檢驗Clayton Copula的適用性，采用N＝1 000，得出P值為0.001，由此可見，Clayton Copula可以很好地描述干旱歷時和烈度的相依性．因 Kendall秩相關(guān)系數(shù)τ約為0.70，Clayton Copula的參數(shù)為

2.3 結(jié)果分析

以皮爾遜曲線和伽馬函數(shù)分別為干旱歷時和干旱烈度的邊際函數(shù)、Clayton Copula為聯(lián)結(jié)函數(shù)，對錢塘江諸暨站的干旱情況做出模擬．圖 5和圖 6分別為 Copula函數(shù)及干旱歷時和干旱烈度的兩元聯(lián)合分布等值線圖，可以看出兩變量存在強相依關(guān)系．從圖5和圖 6即可查出干旱歷時和強度發(fā)生各種組合的聯(lián)合分布，例如，可以查出干旱歷時為3年時，干旱烈度為600,mm的概率為 69%，從而可以定量得出干旱烈度和干旱歷時發(fā)生各種組合的概率．

圖5 Copula函數(shù)及相應(yīng)的等值線Fig.5 Copula and corresponding level curves

圖6 干旱歷時和干旱烈度的兩元聯(lián)合分布及其等值線Fig.6 Joint distribution of drought duration and severity and corresponding level curves

表 3為一些典型干旱歷時與烈度的聯(lián)合概率值，即 PX,Y(X ≤ x , Y ≤ y )的值．此聯(lián)合概率值可為流域/地區(qū)的水資源的有效管理提供強有力的技術(shù)支持．

表3 不同干旱歷時與烈度的聯(lián)合分布Tab.3 Joint distribution for various combinations of drought durations and severity

3 結(jié) 論

(1)利用自回歸馬爾可夫模型來延長水文干旱數(shù)據(jù)，可以解決干旱數(shù)據(jù)稀少的問題，在本文中，自回歸馬爾可夫模型較好地模擬了偏態(tài)的年降雨序列．

(2)Clayton Copula可以很好地模擬正相依的干旱數(shù)據(jù)，是分析極限干旱事件的有用工具．

(3)選用的截取水平也可根據(jù)具體的水資源情況來確定，如實際氣象、農(nóng)業(yè)或水文干旱發(fā)生時的降雨值，不一定是個定值．

［1］ 馮國章. 極限水文干旱歷時概率分析[J]. 水利學(xué)報，1995，6：37-41，50.Feng Guozhang. An analysis of frequency of critical hydrologic drought duration[J]. Journal of Hydraulic Engineering，1995，6：37-41，50(in Chinese).

［2］ 丁 晶，袁 鵬，楊榮富，等. 中國主要河流干旱特性的統(tǒng)計分析[J]. 地理學(xué)報，1997，52(4)：374-381.Ding Jing，Yuan Peng，Yang Rongfu，et al. Stochastic analysis of the drought properties of thema in rivers in China[J]. Acta Geographica Sinica，1997，52(4)：374-381(in Chinese).

［3］ 朱廷舉，胡和平. 基于隨機模擬和模糊聚類的水文干旱特性分析[J]. 清華大學(xué)學(xué)報：自然科學(xué)版，2001，41(8)：101-106.Zhu Tingju，Hu Heping. Drought analysis based on stochastic simulation and fuzzy classification[J]. Journal of Tsinghua University：Natural Science Edition，2001，41(8)：101-106(in Chinese).

［4］ 袁 超，宋松柏，荊 萍. 極限水文干旱歷時概率分布解析法研究[J]. 西北農(nóng)林科技大學(xué)學(xué)報：自然科學(xué)版，2008，36(7)：212-218.Yuan Chao，Song Songbai，Jing Ping. Analytical study on probability distribution of critical hydrological drought duration[J]. Journal of Northwest A and F University：Natural Science Edition，2008，36(7)：212-218(in Chinese).

［5］ González J，Valdés J B. Bivariate drought recurrence analysis using tree ring reconstructions[J]. Journal of Hydrological Engineering，2003，8(5)：247-258.

［6］ Shiau J T，F(xiàn)eng S，Nadarajah S. Assessment of hydrological droughts for the Yellow River，China，using Copulas[J]. Hydrological Processes，2007，21(16)：2157-2163.

［7］ Yevjevich V. Stochastic Processes in Hydrology[M]. Fort Collins：Water Resources Publications，1972.

［8］ 翁文斌，王忠靜，趙建世. 現(xiàn)代水資源規(guī)劃：理論、方法和技術(shù)[M]．北京：清華大學(xué)出版社，2007.Weng Wenbin，Wang Zhongjing，Zhao Jianshi. Advanced Water Resources Planning：Theory，Method and Technology [M]. Beijing：Tsinghua University Press，2007(in Chinese).

［9］ 熊立華，郭生練，肖 義，等. Copula聯(lián)結(jié)函數(shù)在多變量水文頻率分析中的應(yīng)用[J]. 武漢大學(xué)學(xué)報：工學(xué)版，2005，38(6)：16-19.Xiong Lihua，Guo Shenglian，Xiao Yi，et al. Application of Copula to multivariate hydrological frequency analysis[J]. Engineering Journal of Wuhan University：Engineering Science，2005，38(6)：16-19(in Chinese).

［10］ Zhang L，Singh V P. Bivariate flood frequency analysis using the Copula method[J]. Journal of Hydrological Engineering，2006，11(2)：150-164.

［11］ 肖 義，郭生練，熊立華，等. 一種新的洪水過程隨機模擬方法研究[J]. 四川大學(xué)學(xué)報：工程科學(xué)版，2007，39(2)：55-60.Xiao Yi，Guo Shenglian，Xiong Lihua，et al. A new random simulation method for constructing synthetic flood hydrographs[J]. Journal of Sichuan University：Engineering Science Edition，2007，39(2)：55-60(in Chinese).

［12］ 許月萍，李 佳，曹飛鳳，等. Copula在水文極限事件分析中的應(yīng)用[J]. 浙江大學(xué)學(xué)報：工學(xué)版，2008，42(7)：1119-1122.Xu Yueping，Li Jia，Cao Feifeng，et al. Applications of Copula in hydrological extreme analysis[J]. Journal of Zhejiang University：Engineering Science，2008，42(7)：1119-1122(in Chinese).

［13］ Joe H. Asymptotic efficiency of the two-stage estimation method for Copula-based models [J]. Journal of Multivariate Analysis，2005，94(2)：401-419.