基于Copula 的干旱频率分析方法在唐乃亥流域中的应用

杨 涛，奚圆圆，刘 伟

（1.云南省水利水电勘测设计研究院，云南 昆明 650021；2.云南水利水电职业学院，云南 昆明 650021）

1 引言

人类生存的环境经常遭受各类自然灾害的袭击，如干旱、泥石流、台风、洪涝、地震、火山喷发等等。据世界气象组织统计，在世界范围内，气象灾害导致的损失占各类自然灾害总损失的85%之多，其中，约50%的气象灾害损失是由干旱引起的，全球每年因干旱造成的经济损失高达60 亿美元～80 亿美元，居各类气象灾害之首[1]。就我国而言，由于降水在时间和空间上分布的高度不均匀性，进而引起水文极端干旱事件频繁发生，干旱在中国已成为影响区域最广、发生最频繁、造成经济损失最大的气象灾害之一。历史旱灾资料统计表明，从公元前206 年至公元1949 年的2155 年间，共发生1056 次干旱事件，平均每两年就发生一次。在新中国成立以后的1950 年～2010 年共61 年间，全国发生严重、特大旱灾24 次，发生频次为40%，平均每2.5 年发生一次，平均受灾面积达到2.16×107hm2，平均每年成灾面积达9.61×107hm2，平均因灾损失粮食1.61×107t[2]。

黄河上游唐乃亥站以上区域位于青藏高原东北部，流域面积12.20 万km2，占黄河流域面积的16.2%，多年平均来水量200 亿m3，占黄河流域多年平均来水量的34.5%。流域地势总的趋势是西高东低，境内分布着大小冰川40 多条，冰川总面积120.57 km2，是黄河流域的天然固体水库和主要水源。该区多年平均降水量400 mm～700 mm，雨量分布不均，干旱寒冷，年平均气温为-4℃，具有显著的寒缺氧、气温低、光辐射强、昼夜温差大等典型的高原大陆性气候特点。据统计，从1990 至今，黄河已出现多次断流：1996 年和1998 年，黄河在“黄河源头第一县”的玛多县出现两次断流；1999 年5 月2 日至6 月3 日，扎陵湖和鄂陵湖之间的河道出现断流，河道干涸达8 km；1998 年从鄂陵湖以下黄河出现60 km 的断流。2001 年玛多县遭遇了百年不遇的特大旱情，7 月份以后降水比上年同期偏少近71%，大多牧草未结草籽便提前枯死，之前绿草如茵的地方，现已变成了荒漠和沙地，此次干旱造成玛多县1599 户7640 名牧民和38 万头牲畜受灾，2095 万亩夏秋草场可利用率不到60%，30 万牲畜无法安全越冬。

干旱的发生是自然变化和区域人类活动共同作用的产物，它已成为阻碍社会发展、扰乱人们正常生活的焦点问题，迫切需要加强对抗旱减灾和旱灾风险管理的研究，而目前研究者们对该区域水文干旱和频率分析的研究较少。干旱是包括干旱历时和干旱烈度等多个相关变量的极值水文事件，在干旱频率分析研究中，干旱历时与干旱烈度之间的关系是重点，必须加以考虑，而copula 函数是一种描述变量之间相关性的有效方法[3]。

2 数据和方法

2.1 数据

本文数据主要是黄河源区唐乃亥水文站45 a（1956 年～2000 年）的逐月天然径流资料，数据来源于水文统计年鉴。图1 是唐乃亥水文站1956 年～2000 年年径流过程。

图1 唐乃亥水文站1956 年～2000 年年径流过程

2.2 水文干旱识别

游程理论是目前水文干旱识别的主要理论，它因为无需假定变量所服从的概率分布，均可直接从简单的统计游程现象的方法入手，揭示游程现象发生的概率[4～5]。在实际应用中，依据径流量距平百分率确定月径流量的三个阈值R0，R1，R2，当单月径流量小于R1初步判定为干旱，有a，b，c，d 四个干旱过程；对历时等于1 的干旱过程，若指标小于R2则视为一次干旱，若小于R1但大于R2则不视为一次干旱过程；对于两次干旱过程，若间隔等于1 且该指标值小于R0，则合并为一次干旱过程，其干旱烈度S=Sb+Sc，反之为两次干旱过程，根据《水文情报预报规范》中的标准，径流量距平百分率的三个阈值分别为R0=0，R1=-20%，R2=-40%。

2.3 基于Copula 函数的干旱频率分析

（1）Copula 函数

Copula 是连接一维边际分布形成在[0,1]上的多元分布的联合分布函数。Copula 函数能够将二元联合分布划分为变量的边缘分布和变量之间的相关性结构来分别进行处理，其中用Copula 函数来描述变量间的相关性结构。任意形式的边缘分布均可通过Copula 函数来构造其联合分布函数，形式灵活多样，计算求解较易。由于边缘分布包含变量的所有信息，在转换过程中信息不会失真，被广泛地应用在多变量水文频率分析计算中[6]。设x、y 为连续的随机变量，其边缘分布函数分别为Fx和Fy，F（x,y）为变量x 和变量y 的联合分布函数，那么存在唯一的Copula 函数C，使得F（x,y）=Cθ（Fx（x）,Fy（y）），式中：（Fx（x）,Fy（y））为Copula 函数，θ 为待定参数，可以通过与Kendall 秩相关系数τ 的关系确定。

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

一般情况下，假设某一特征变量服从某种分布，通过参数估计得到该分布概率密度函数，再由假设检验判断该假设的正确性。假定干旱历时服从指数分布，干旱烈度服从2 参数的Gamma 分布[7]。经过计算得出干旱历时和干旱烈度的参数，然后构造统计量，采用χ2检验验证拟合度。

（3）Copula 联合分布函数的选择

采用Kolmogorov-Smimov（K-S）检验来评价联合分布计算与联合观测值的拟合程度[8]，表征描述干旱特征变量之间的相关性，采用均方根误差RMSE 作为拟合优度的评价指标。

①Copula 函数的拟合检验，采用Kolmogorov-Smimov（K-S）检验，其统计量D 计算如下：

式中：F（xi，yi）为（x，y）的联合分布；m（i）为观测样本中满足x≤xi，y≤yi条件的联合观测值个数。

②采用均方根误差（RMSE）来评定Copula 函数拟合结果优劣，其表达式为：

式中：Pc（i）是Copula 函数计算得的第i 个联合观测值（di，si）的联合分布概率值；n 为联合观测样本数；P0（i）为多元联合分布经验值。

式中：mi表示联合观测样本中满足条件D≤di且S≤si的联合观测值的个数。

（4）干旱联合重现期

干旱重现期作为描述随机事件在时间上稀遇程度的统计量，是水文极值事件的重点研究对象，也是干旱频率分析和水文频率分析的重要内容[9]。一般定义为：具有某种属性的随机事件出现一次的平均间隔时间[10]。然而，干旱事件不同于其它水文极值事件，持续时间长，发展缓慢且不易被发现，一次干旱事件可能持续多年或者一年中发生多次干旱过程。根据现有的研究成果，常采用超定量法选样来描述其重现期。据此，Shiau 和Shen[11]推导出了干旱重现期的计算公式。干旱历时重现期：TD=E（L）/[1-FD（d）]；干旱烈度TS重现期：TS=E（L）/[1-FS（s）]，式中，E（L）为干旱间隔的期望值，它等于干旱历时和非干旱历时的期望值之和。

当同时考虑干旱历时与干旱烈度时，干旱事件的重现期分两种情况计算：

①在干旱历时D＞d，或者干旱烈度S＞s 时的干旱重现期为：

②当干旱历时D＞d，且干旱烈度S＞s 时，干旱重现期为：

当识别出干旱特征系列后，将对应的干旱历时和干旱烈度值代入上式即可求得对应的干旱重现期。

3 结果分析

3.1 水文干旱识别结果

基于游程理论的水文干旱识别方法，本文对唐乃亥站1956 年～2000 年逐月天然径流数据进行干旱序列提取，其结果见表1 和表2。由表1 和表2 可知：（1）唐乃亥站以上黄河源区1956 年～2000 年45 a 内共发生29 次干旱，从1956 年7 月～1957 年12 月发生了持续18 个月的严重干旱，干旱烈度是45年内最大的；干旱烈度值最小值为30093 万m3，发生在1999年4 月和5 月，历时为两个月。干旱的历时均值为4.7 个月，干旱间隔期望值为18.62 个月；（2）仅1958 年、1961 年、1964 年、1967 年、1968 年、1973 年～1976 年、1981 年～1984 年、1989 年、1992 年、1993 年共计16 年未发生干旱，其中1967 年～1968 年年连续2 年无干旱事件，1973 年～1976 年连续4 年无旱，1981 年～1984 年连续4 年无旱；（3）1994 年～2000 年每年均发生干旱，平均每年干旱持续时间为5 个月，平均干旱烈度为183397 万m3，较45 a 年的多年平均干旱烈度高25098 万m3。

表1 唐乃亥站1956 年～2000 年月尺度干旱统计结果

表2 月尺度干旱特征值统计结果

3.2 二维联合Copula 函数结果分析

（1）边缘分布函数的确定

在假定干旱历时服从指数分布，干旱烈度服从二参数Gamma 分布的前提下，分别将干旱历时系列和干旱烈度系列样本采用矩法进行参数估计，得出参数α、β。并选择统计量对干旱特征系列进行假设检验。经过对干旱特征变量边际分布假设检验，其结果表明干旱历时服从参数为α=4.66 的指数分布，干旱烈度服从参数α=1.21、参数β=130913 二参数Gamma 分布。

表3 唐乃亥干旱特征变量的参数估计

表4 干旱特征变量边际分布假设检验

（2）Copula 联合分布函数的确定

Copula 函数具有构造方便、容易使用等优点，本文选用Gumbel-Copul、Clayton-Copula 和Frank-Copula 这三 类Copula函数构建干旱历时和干旱烈度的联合分布函数，采用K-S 检验来评价联合分布计算值与联合观测值的拟合程度，表征干旱特征变量之间的相关性，采用均方根误差RMSE 作为拟合优度的评价指标。K-S 的值越低，其理论分布与样本序列的经验分布拟合越好，即无显著差异；均方根误差RMSE 值越小，表明该Copula 函数最优。从表5 数据对比可知，构建的三类Copula 函数的kendall 相关系数均为0.537，再从均方根误差RMSE 值来判断，Gumbel-Copula 的均方根误差RMSE 值为0.044，在所有Copula 函数中最小，且根据图3 干旱特征序列点据与45°斜线拟合程度最好，故Gumbel-Copula 函数拟合结果相对较优，更适合唐乃亥流域干旱特征变量的联合分布概率函数，其二维联合分布概率函数见图4。

表5 基于K-S 和RMSE 的Copula 函数最优选择

图3 三类Copula 函数干旱特征变量的理论与经验联合概率分布对比图

图4 基于Gumbel-Copula 函数的干旱历时与干旱烈度的联合概率分布

表6 干旱特征系列的重现期结果统计 单位干旱历时：月；干旱烈度：万m3；联合重现期：月

（3）干旱联合重现期

针对提取出来的干旱特征值系列，将近45 年来发生的29次干旱事件的各次干旱历时与干旱烈度系列代入式（4）和式（5），经计算，干旱间隔的期望值E（L）=18.62 月。由表6 中数据可知：①1956 年7 月～1957 年12 月发生的这次干旱是唐乃亥站45 a 来最严重的一次干旱，干旱历时持续18 个月，其联合重现期为115 个月；1959 年9 月～1960 年6 月的这次干旱联合重现期为31 个月，在这45 年排在第二位。②29 次干旱事件的重现期平均值为11 个月。

4 结语

（1）唐乃亥站以上黄河源区1956 年～2000 年45 a 内共发生29 次干旱，干旱的历时均值为4.7 个月，干旱间隔期望值为18.62 个月。从1956 年7 月～1957 年12 月发生了持续18 个月的严重干旱，干旱烈度是45 年内最大的。干旱烈度值最小值为30093 万m3，发生在1999 年4 月和5 月，历时为2 个月。另外根据统计结果发现1994 年～2000 年每年均发生干旱，平均每年干旱持续时间为5 个月。

（2）通过对Gumbel-Copula、Clayton-Copula、Frank-Copula三种连接函数用K-S 和RMSE 对比分析发现Gumbel-Copula函数拟合结果相对较优，故Gumbel-Copula 函数是用于研究唐乃亥站以上流域水文干旱特征的最佳分布概率函数。

（3）通过对干旱重现期的计算，结果表明29 次干旱事件的重现期平均值为11 个月，即发生一次干旱事件的平均间隔时间为11 个月。1956.7～1957.12 发生的这次干旱是唐乃亥站以上流域45 a 来最严重的一次干旱，干旱历时持续18 个月，其重现期最长，为115 个月。