基于Copula 的干旱頻率分析方法在唐乃亥流域中的應(yīng)用

楊 濤，奚圓圓，劉 偉

（1.云南省水利水電勘測(cè)設(shè)計(jì)研究院，云南 昆明 650021；2.云南水利水電職業(yè)學(xué)院，云南 昆明 650021）

1 引言

人類生存的環(huán)境經(jīng)常遭受各類自然災(zāi)害的襲擊，如干旱、泥石流、臺(tái)風(fēng)、洪澇、地震、火山噴發(fā)等等。據(jù)世界氣象組織統(tǒng)計(jì)，在世界范圍內(nèi)，氣象災(zāi)害導(dǎo)致的損失占各類自然災(zāi)害總損失的85%之多，其中，約50%的氣象災(zāi)害損失是由干旱引起的，全球每年因干旱造成的經(jīng)濟(jì)損失高達(dá)60 億美元～80 億美元，居各類氣象災(zāi)害之首[1]。就我國(guó)而言，由于降水在時(shí)間和空間上分布的高度不均勻性，進(jìn)而引起水文極端干旱事件頻繁發(fā)生，干旱在中國(guó)已成為影響區(qū)域最廣、發(fā)生最頻繁、造成經(jīng)濟(jì)損失最大的氣象災(zāi)害之一。歷史旱災(zāi)資料統(tǒng)計(jì)表明，從公元前206 年至公元1949 年的2155 年間，共發(fā)生1056 次干旱事件，平均每?jī)赡昃桶l(fā)生一次。在新中國(guó)成立以后的1950 年～2010 年共61 年間，全國(guó)發(fā)生嚴(yán)重、特大旱災(zāi)24 次，發(fā)生頻次為40%，平均每2.5 年發(fā)生一次，平均受災(zāi)面積達(dá)到2.16×107hm2，平均每年成災(zāi)面積達(dá)9.61×107hm2，平均因?yàn)?zāi)損失糧食1.61×107t[2]。

黃河上游唐乃亥站以上區(qū)域位于青藏高原東北部，流域面積12.20 萬(wàn)km2，占黃河流域面積的16.2%，多年平均來(lái)水量200 億m3，占黃河流域多年平均來(lái)水量的34.5%。流域地勢(shì)總的趨勢(shì)是西高東低，境內(nèi)分布著大小冰川40 多條，冰川總面積120.57 km2，是黃河流域的天然固體水庫(kù)和主要水源。該區(qū)多年平均降水量400 mm～700 mm，雨量分布不均，干旱寒冷，年平均氣溫為-4℃，具有顯著的寒缺氧、氣溫低、光輻射強(qiáng)、晝夜溫差大等典型的高原大陸性氣候特點(diǎn)。據(jù)統(tǒng)計(jì)，從1990 至今，黃河已出現(xiàn)多次斷流：1996 年和1998 年，黃河在“黃河源頭第一縣”的瑪多縣出現(xiàn)兩次斷流；1999 年5 月2 日至6 月3 日，扎陵湖和鄂陵湖之間的河道出現(xiàn)斷流，河道干涸達(dá)8 km；1998 年從鄂陵湖以下黃河出現(xiàn)60 km 的斷流。2001 年瑪多縣遭遇了百年不遇的特大旱情，7 月份以后降水比上年同期偏少近71%，大多牧草未結(jié)草籽便提前枯死，之前綠草如茵的地方，現(xiàn)已變成了荒漠和沙地，此次干旱造成瑪多縣1599 戶7640 名牧民和38 萬(wàn)頭牲畜受災(zāi)，2095 萬(wàn)畝夏秋草場(chǎng)可利用率不到60%，30 萬(wàn)牲畜無(wú)法安全越冬。

干旱的發(fā)生是自然變化和區(qū)域人類活動(dòng)共同作用的產(chǎn)物，它已成為阻礙社會(huì)發(fā)展、擾亂人們正常生活的焦點(diǎn)問(wèn)題，迫切需要加強(qiáng)對(duì)抗旱減災(zāi)和旱災(zāi)風(fēng)險(xiǎn)管理的研究，而目前研究者們對(duì)該區(qū)域水文干旱和頻率分析的研究較少。干旱是包括干旱歷時(shí)和干旱烈度等多個(gè)相關(guān)變量的極值水文事件，在干旱頻率分析研究中，干旱歷時(shí)與干旱烈度之間的關(guān)系是重點(diǎn)，必須加以考慮，而copula 函數(shù)是一種描述變量之間相關(guān)性的有效方法[3]。

2 數(shù)據(jù)和方法

2.1 數(shù)據(jù)

本文數(shù)據(jù)主要是黃河源區(qū)唐乃亥水文站45 a（1956 年～2000 年）的逐月天然徑流資料，數(shù)據(jù)來(lái)源于水文統(tǒng)計(jì)年鑒。圖1 是唐乃亥水文站1956 年～2000 年年徑流過(guò)程。

圖1 唐乃亥水文站1956 年～2000 年年徑流過(guò)程

2.2 水文干旱識(shí)別

游程理論是目前水文干旱識(shí)別的主要理論，它因?yàn)闊o(wú)需假定變量所服從的概率分布，均可直接從簡(jiǎn)單的統(tǒng)計(jì)游程現(xiàn)象的方法入手，揭示游程現(xiàn)象發(fā)生的概率[4～5]。在實(shí)際應(yīng)用中，依據(jù)徑流量距平百分率確定月徑流量的三個(gè)閾值R0，R1，R2，當(dāng)單月徑流量小于R1初步判定為干旱，有a，b，c，d 四個(gè)干旱過(guò)程；對(duì)歷時(shí)等于1 的干旱過(guò)程，若指標(biāo)小于R2則視為一次干旱，若小于R1但大于R2則不視為一次干旱過(guò)程；對(duì)于兩次干旱過(guò)程，若間隔等于1 且該指標(biāo)值小于R0，則合并為一次干旱過(guò)程，其干旱烈度S=Sb+Sc，反之為兩次干旱過(guò)程，根據(jù)《水文情報(bào)預(yù)報(bào)規(guī)范》中的標(biāo)準(zhǔn)，徑流量距平百分率的三個(gè)閾值分別為R0=0，R1=-20%，R2=-40%。

2.3 基于Copula 函數(shù)的干旱頻率分析

（1）Copula 函數(shù)

Copula 是連接一維邊際分布形成在[0,1]上的多元分布的聯(lián)合分布函數(shù)。Copula 函數(shù)能夠?qū)⒍?lián)合分布劃分為變量的邊緣分布和變量之間的相關(guān)性結(jié)構(gòu)來(lái)分別進(jìn)行處理，其中用Copula 函數(shù)來(lái)描述變量間的相關(guān)性結(jié)構(gòu)。任意形式的邊緣分布均可通過(guò)Copula 函數(shù)來(lái)構(gòu)造其聯(lián)合分布函數(shù)，形式靈活多樣，計(jì)算求解較易。由于邊緣分布包含變量的所有信息，在轉(zhuǎn)換過(guò)程中信息不會(huì)失真，被廣泛地應(yīng)用在多變量水文頻率分析計(jì)算中[6]。設(shè)x、y 為連續(xù)的隨機(jī)變量，其邊緣分布函數(shù)分別為Fx和Fy，F(xiàn)（x,y）為變量x 和變量y 的聯(lián)合分布函數(shù)，那么存在唯一的Copula 函數(shù)C，使得F（x,y）=Cθ（Fx（x）,Fy（y）），式中：（Fx（x）,Fy（y））為Copula 函數(shù)，θ 為待定參數(shù)，可以通過(guò)與Kendall 秩相關(guān)系數(shù)τ 的關(guān)系確定。

（2）干旱特征變量概率分布

一般情況下，假設(shè)某一特征變量服從某種分布，通過(guò)參數(shù)估計(jì)得到該分布概率密度函數(shù)，再由假設(shè)檢驗(yàn)判斷該假設(shè)的正確性。假定干旱歷時(shí)服從指數(shù)分布，干旱烈度服從2 參數(shù)的Gamma 分布[7]。經(jīng)過(guò)計(jì)算得出干旱歷時(shí)和干旱烈度的參數(shù)，然后構(gòu)造統(tǒng)計(jì)量，采用χ2檢驗(yàn)驗(yàn)證擬合度。

（3）Copula 聯(lián)合分布函數(shù)的選擇

采用Kolmogorov-Smimov（K-S）檢驗(yàn)來(lái)評(píng)價(jià)聯(lián)合分布計(jì)算與聯(lián)合觀測(cè)值的擬合程度[8]，表征描述干旱特征變量之間的相關(guān)性，采用均方根誤差RMSE 作為擬合優(yōu)度的評(píng)價(jià)指標(biāo)。

①Copula 函數(shù)的擬合檢驗(yàn)，采用Kolmogorov-Smimov（K-S）檢驗(yàn)，其統(tǒng)計(jì)量D 計(jì)算如下：

式中：F（xi，yi）為（x，y）的聯(lián)合分布；m（i）為觀測(cè)樣本中滿足x≤xi，y≤yi條件的聯(lián)合觀測(cè)值個(gè)數(shù)。

②采用均方根誤差（RMSE）來(lái)評(píng)定Copula 函數(shù)擬合結(jié)果優(yōu)劣，其表達(dá)式為：

式中：Pc（i）是Copula 函數(shù)計(jì)算得的第i 個(gè)聯(lián)合觀測(cè)值（di，si）的聯(lián)合分布概率值；n 為聯(lián)合觀測(cè)樣本數(shù)；P0（i）為多元聯(lián)合分布經(jīng)驗(yàn)值。

式中：mi表示聯(lián)合觀測(cè)樣本中滿足條件D≤di且S≤si的聯(lián)合觀測(cè)值的個(gè)數(shù)。

（4）干旱聯(lián)合重現(xiàn)期

干旱重現(xiàn)期作為描述隨機(jī)事件在時(shí)間上稀遇程度的統(tǒng)計(jì)量，是水文極值事件的重點(diǎn)研究對(duì)象，也是干旱頻率分析和水文頻率分析的重要內(nèi)容[9]。一般定義為：具有某種屬性的隨機(jī)事件出現(xiàn)一次的平均間隔時(shí)間[10]。然而，干旱事件不同于其它水文極值事件，持續(xù)時(shí)間長(zhǎng)，發(fā)展緩慢且不易被發(fā)現(xiàn)，一次干旱事件可能持續(xù)多年或者一年中發(fā)生多次干旱過(guò)程。根據(jù)現(xiàn)有的研究成果，常采用超定量法選樣來(lái)描述其重現(xiàn)期。據(jù)此，Shiau 和Shen[11]推導(dǎo)出了干旱重現(xiàn)期的計(jì)算公式。干旱歷時(shí)重現(xiàn)期：TD=E（L）/[1-FD（d）]；干旱烈度TS重現(xiàn)期：TS=E（L）/[1-FS（s）]，式中，E（L）為干旱間隔的期望值，它等于干旱歷時(shí)和非干旱歷時(shí)的期望值之和。

當(dāng)同時(shí)考慮干旱歷時(shí)與干旱烈度時(shí)，干旱事件的重現(xiàn)期分兩種情況計(jì)算：

①在干旱歷時(shí)D＞d，或者干旱烈度S＞s 時(shí)的干旱重現(xiàn)期為：

②當(dāng)干旱歷時(shí)D＞d，且干旱烈度S＞s 時(shí)，干旱重現(xiàn)期為：

當(dāng)識(shí)別出干旱特征系列后，將對(duì)應(yīng)的干旱歷時(shí)和干旱烈度值代入上式即可求得對(duì)應(yīng)的干旱重現(xiàn)期。

3 結(jié)果分析

3.1 水文干旱識(shí)別結(jié)果

基于游程理論的水文干旱識(shí)別方法，本文對(duì)唐乃亥站1956 年～2000 年逐月天然徑流數(shù)據(jù)進(jìn)行干旱序列提取，其結(jié)果見(jiàn)表1 和表2。由表1 和表2 可知：（1）唐乃亥站以上黃河源區(qū)1956 年～2000 年45 a 內(nèi)共發(fā)生29 次干旱，從1956 年7 月～1957 年12 月發(fā)生了持續(xù)18 個(gè)月的嚴(yán)重干旱，干旱烈度是45年內(nèi)最大的；干旱烈度值最小值為30093 萬(wàn)m3，發(fā)生在1999年4 月和5 月，歷時(shí)為兩個(gè)月。干旱的歷時(shí)均值為4.7 個(gè)月，干旱間隔期望值為18.62 個(gè)月；（2）僅1958 年、1961 年、1964 年、1967 年、1968 年、1973 年～1976 年、1981 年～1984 年、1989 年、1992 年、1993 年共計(jì)16 年未發(fā)生干旱，其中1967 年～1968 年年連續(xù)2 年無(wú)干旱事件，1973 年～1976 年連續(xù)4 年無(wú)旱，1981 年～1984 年連續(xù)4 年無(wú)旱；（3）1994 年～2000 年每年均發(fā)生干旱，平均每年干旱持續(xù)時(shí)間為5 個(gè)月，平均干旱烈度為183397 萬(wàn)m3，較45 a 年的多年平均干旱烈度高25098 萬(wàn)m3。

表1 唐乃亥站1956 年～2000 年月尺度干旱統(tǒng)計(jì)結(jié)果

表2 月尺度干旱特征值統(tǒng)計(jì)結(jié)果

3.2 二維聯(lián)合Copula 函數(shù)結(jié)果分析

（1）邊緣分布函數(shù)的確定

在假定干旱歷時(shí)服從指數(shù)分布，干旱烈度服從二參數(shù)Gamma 分布的前提下，分別將干旱歷時(shí)系列和干旱烈度系列樣本采用矩法進(jìn)行參數(shù)估計(jì)，得出參數(shù)α、β。并選擇統(tǒng)計(jì)量對(duì)干旱特征系列進(jìn)行假設(shè)檢驗(yàn)。經(jīng)過(guò)對(duì)干旱特征變量邊際分布假設(shè)檢驗(yàn)，其結(jié)果表明干旱歷時(shí)服從參數(shù)為α=4.66 的指數(shù)分布，干旱烈度服從參數(shù)α=1.21、參數(shù)β=130913 二參數(shù)Gamma 分布。

表3 唐乃亥干旱特征變量的參數(shù)估計(jì)

表4 干旱特征變量邊際分布假設(shè)檢驗(yàn)

（2）Copula 聯(lián)合分布函數(shù)的確定

Copula 函數(shù)具有構(gòu)造方便、容易使用等優(yōu)點(diǎn)，本文選用Gumbel-Copul、Clayton-Copula 和Frank-Copula 這三 類Copula函數(shù)構(gòu)建干旱歷時(shí)和干旱烈度的聯(lián)合分布函數(shù)，采用K-S 檢驗(yàn)來(lái)評(píng)價(jià)聯(lián)合分布計(jì)算值與聯(lián)合觀測(cè)值的擬合程度，表征干旱特征變量之間的相關(guān)性，采用均方根誤差RMSE 作為擬合優(yōu)度的評(píng)價(jià)指標(biāo)。K-S 的值越低，其理論分布與樣本序列的經(jīng)驗(yàn)分布擬合越好，即無(wú)顯著差異；均方根誤差RMSE 值越小，表明該Copula 函數(shù)最優(yōu)。從表5 數(shù)據(jù)對(duì)比可知，構(gòu)建的三類Copula 函數(shù)的kendall 相關(guān)系數(shù)均為0.537，再?gòu)木礁`差RMSE 值來(lái)判斷，Gumbel-Copula 的均方根誤差RMSE 值為0.044，在所有Copula 函數(shù)中最小，且根據(jù)圖3 干旱特征序列點(diǎn)據(jù)與45°斜線擬合程度最好，故Gumbel-Copula 函數(shù)擬合結(jié)果相對(duì)較優(yōu)，更適合唐乃亥流域干旱特征變量的聯(lián)合分布概率函數(shù)，其二維聯(lián)合分布概率函數(shù)見(jiàn)圖4。

表5 基于K-S 和RMSE 的Copula 函數(shù)最優(yōu)選擇

圖3 三類Copula 函數(shù)干旱特征變量的理論與經(jīng)驗(yàn)聯(lián)合概率分布對(duì)比圖

圖4 基于Gumbel-Copula 函數(shù)的干旱歷時(shí)與干旱烈度的聯(lián)合概率分布

表6 干旱特征系列的重現(xiàn)期結(jié)果統(tǒng)計(jì) 單位干旱歷時(shí)：月；干旱烈度：萬(wàn)m3；聯(lián)合重現(xiàn)期：月

（3）干旱聯(lián)合重現(xiàn)期

針對(duì)提取出來(lái)的干旱特征值系列，將近45 年來(lái)發(fā)生的29次干旱事件的各次干旱歷時(shí)與干旱烈度系列代入式（4）和式（5），經(jīng)計(jì)算，干旱間隔的期望值E（L）=18.62 月。由表6 中數(shù)據(jù)可知：①1956 年7 月～1957 年12 月發(fā)生的這次干旱是唐乃亥站45 a 來(lái)最嚴(yán)重的一次干旱，干旱歷時(shí)持續(xù)18 個(gè)月，其聯(lián)合重現(xiàn)期為115 個(gè)月；1959 年9 月～1960 年6 月的這次干旱聯(lián)合重現(xiàn)期為31 個(gè)月，在這45 年排在第二位。②29 次干旱事件的重現(xiàn)期平均值為11 個(gè)月。

4 結(jié)語(yǔ)

（1）唐乃亥站以上黃河源區(qū)1956 年～2000 年45 a 內(nèi)共發(fā)生29 次干旱，干旱的歷時(shí)均值為4.7 個(gè)月，干旱間隔期望值為18.62 個(gè)月。從1956 年7 月～1957 年12 月發(fā)生了持續(xù)18 個(gè)月的嚴(yán)重干旱，干旱烈度是45 年內(nèi)最大的。干旱烈度值最小值為30093 萬(wàn)m3，發(fā)生在1999 年4 月和5 月，歷時(shí)為2 個(gè)月。另外根據(jù)統(tǒng)計(jì)結(jié)果發(fā)現(xiàn)1994 年～2000 年每年均發(fā)生干旱，平均每年干旱持續(xù)時(shí)間為5 個(gè)月。

（2）通過(guò)對(duì)Gumbel-Copula、Clayton-Copula、Frank-Copula三種連接函數(shù)用K-S 和RMSE 對(duì)比分析發(fā)現(xiàn)Gumbel-Copula函數(shù)擬合結(jié)果相對(duì)較優(yōu)，故Gumbel-Copula 函數(shù)是用于研究唐乃亥站以上流域水文干旱特征的最佳分布概率函數(shù)。

（3）通過(guò)對(duì)干旱重現(xiàn)期的計(jì)算，結(jié)果表明29 次干旱事件的重現(xiàn)期平均值為11 個(gè)月，即發(fā)生一次干旱事件的平均間隔時(shí)間為11 個(gè)月。1956.7～1957.12 發(fā)生的這次干旱是唐乃亥站以上流域45 a 來(lái)最嚴(yán)重的一次干旱，干旱歷時(shí)持續(xù)18 個(gè)月，其重現(xiàn)期最長(zhǎng)，為115 個(gè)月。