基于Copula的汾河上游水文干旱频率分析

任　璐，赵雪花

基于Copula的汾河上游水文干旱频率分析

任璐，赵雪花

（太原理工大学水利科学与工程学院，山西太原 030024）

运用游程理论对汾河上游4个水文站的月径流资料进行干旱识别；选择4种分布函数拟合干旱特征变量，通过Kolmogorov-Smirnov检验法优选单变量边缘分布；利用Copula函数建立干旱特征变量的二维联合分布，计算重现期。结果表明，汾河上游干旱特征变量的最佳边缘分布是对数正态分布；汾河上游大多数干旱重现期小于5年，且干旱历时与烈度的单变量重现期与二维重现期较为接近；干旱历时与烈度峰值、烈度与烈度峰值的单变量重现期与二维重现期差距较大。

水文干旱；频率分析；Copula函数；月径流；汾河上游

0　引言

干旱是一种复杂、多因素并具有全球影响力的自然灾害［1］。它经常发生却甚少为人了解。我国属于东亚季风气候，年季间季风的不稳定性是造成我国干旱频发的一个主要原因［2］。气候变化和人类活动的共同作用，导致干旱的发生频率不断上升、造成损失不断加大，严重制约了我国社会经济的发展。缓解现今严峻的干旱形势迫在眉睫，因此对于干旱事件的频率分析是极其必要的。

最初的多变量干旱频率的研究建立在一系列假定的基础上。即，各变量的边缘分布必须服从同种分布［3-6］或者联合分布必须采用正态分布和经过数学转换得到的正态分布［6-8］。然而，实际情况往往不能满足以上假定，难以应用到干旱频率分析中。随后出现的非参数研究方法［9-10］虽然可避免以上假定，但是其数学推导及参数计算都太过复杂，在实际应用中难以推广。相比以上几种方法，Copula函数法既可以构建服从不同边缘分布的特征变量的联合分布函数［11］，又便于计算，更符合干旱频率分析的实际情况；因此可用于水文干旱频率分析中。由于水文现象具有很明显的区域性，不同地区水文情况不同，干旱特征变量的边缘分布类型也不尽相同。部分研究者直接拟定单变量边缘分布类型，可能会给后续研究造成不必要的误差。前人对于水文干旱频率的研究多集中于广东、新疆等省份或者黑河、鄱阳湖等流域，鲜有对汾河上游的水文干旱频率的研究。因此，本文先选择4种水文分析中常见的分布对汾河上游的干旱特征变量进行拟合，再从中优选出适合各特征变量的边缘分布，在此基础上利用Copula函数构建二维干旱变量的联合分布函数，并计算从干旱序列中提取出的代表数对应的单变量重现期、二维联合重现期和同现重现期，以此来分析汾河上游的水文干旱现象。

本文选取汾河上游的4个典型水文站：上静游（1955年 ～2012年）、汾河水库（1958年 ～2000年）、寨上（1956年 ～2000年）、兰村（1951年 ～2012年）的月径流资料，对汾河上游的水文干旱情况进行分析，以为该区域今后的农业发展、工业生产以及防旱抗旱工作提供科学依据。

1　研究方法

1.1干旱事件的定义与识别

利用游程理论［12］提取干旱特征变量：选取月径流距平百分比为-10%作为截取水平R1，当径流序列Ri（t=1，2，…，N）在一个或多个时段内连续小于R1时，出现负游程，即认为发生干旱；负游程的长度作为干旱历时D（月）；负游程的阴影部分面积作为干旱烈度S（万m3）；负游程的极值作为干旱烈度峰值M（万m3）。

实际分析中常出现这样一种情况：一些历时很短的非干旱事件夹杂在一场严重干旱事件之间，导致一场严重干旱被分隔为多场一般干旱，不能反映干旱的实际情况，因此需要按照一定的准则和方法对这些干旱进行合并［13］。取R0为径流距平百分比为0时所对应的径流值，R2为径流距平百分比为 -50%时所对应的径流值。如图1所示，R0、R1、R2将径流量划分为4部分，则按以下方法判断是否发生干旱：①当径流序列中出现单月径流量R＜R2时，认为这个月为一场干旱事件（如图中AB段）；②当径流序列中出现单月径流量R2＜R＜R1时，并不认为这个月是一次干旱事件（如图中IJ段）；③两次干旱事件之间只间隔一个月且该月的月径流量R＜R0时，认为这两个干旱事件有关联，应合并为一场干旱事件（如图中CD段和EF段应合并为一场干旱事件，干旱历时D2=DCD+DDE+DEF，干旱烈度S= S2+S2’）；④当两次干旱事件之间只间隔一个月且该月的月径流量R＞R0时，认为其仍为两场干旱事件（如图中EF段和GH段不需要合并为一场干旱）。

图1　游程理论识别干旱变量示意

1.2干旱特征变量的边缘分布

本文选取对数正态分布（两参数，即2P，下同）、指数分布（2P）、伽马分布（2P）和广义帕累托分布（3P）4种分布对汾河上游的4个水文站的干旱历时、烈度和烈度峰值进行拟合，具体表达式如下：

式中，α为尺度参数；ε为位置参数；k为形状参数；erf为误差函数，即为伽马函数

采用矩法进行参数估计，利用 Kolmogorovmirnov检验法（K-S法）［14］进行拟合优度检验，得到最优的单变量边缘分布函数。

1.3二维干旱特征变量模型

1.3.1Copula函数定义

Copula可以连接一维的边缘分布，形成在［0，1］上的多元分布的函数［14］。设X，Y为连续的随机变量，其边缘分布函数为FX和FY，F（x，y）为变量X和Y的联合分布函数，则存在唯一的Copula函数C，使得

式中，θ为待定参数。

1.3.2二维干旱特征变量的函数选择、参数估计及拟合度检验

本文选用Clayton Copula、Gumbel-Hougaard（GH）Copula、Ali-Mikhail-Haq（AMH）Copula和 Nelson NO.2 Copula 4种Copula函数拟合二维干旱特征变量的联合分布。利用 Kendall秩相关系数 τ和Archimedean Copula函数的参数 θ的解析关系［14，15］，对参数θ进行估计。利用均方根误差（Root Mean Square Error，RMSE）准则和赤池信息量准则（Akaike Information Criterion，AIC）［14］进行拟合优度检验，依据是RMSE和AIC值最小，表达式：

式中，n为样本容量；m为Copula参数估算数；Pc为Copula多元联合分布计算值；P0为多元联合分布经验值。

1.4干旱重现期计算

1.4.1单变量重现期

若干旱事件系列（N年）发生n次，则发生一次干旱事件的年数（重现期）［10］可记为

式中，F（x）为干旱特征变量的边缘分布，N为研究站点的径流资料长度，n为干旱事件次数。

可推导出干旱历时、烈度和烈度峰值对应的单变量重现期，见文献［16］。

1.4.2二维干旱变量重现期

二维干旱变量重现期包括联合重现期与同现重现期。联合重现期Tα与同现重现期T0的计算见文献［14］。

2　实例分析

2.1单变量的边缘分布

根据游程理论对汾河上游的上静游、汾河水库、寨上、兰村4个水文站的月径流资料进行干旱识别，结果见表1。

由表1可看出，在调查年限内汾河上游的4个水文站都曾发生过多场干旱，并且都发生过历时超过24个月或烈度超过1亿m3或峰值超过300万m3的严重干旱，越靠近中下游，发生过的最严重的干旱烈度越大，峰值越高。

采用对数正态分布、指数分布、伽马分布和广义帕累托分布分别对汾河上游的4个水文站的干旱历时、烈度和烈度峰值进行拟合，计算公式见式（1）～式（4）。采用矩法进行参数估计，利用K-S法进行拟合优度检验，依据是K-S统计量的值小于临界值，结果见表2。

表1　4个水文站的基本信息及干旱识别结果

表2　干旱特征变量的4种概率分布的K-S检验结果

由表2可看出，①对于汾河上游4个站的干旱特征变量的拟合，对数正态分布全部通过K-S检验，拟合效果最佳；指数分布有50%通过K-S检验，广义帕累托分布有33.3%通过K-S检验，拟合效果次之，伽马分布全部未通过K-S检验，拟合效果最差。②上静游的干旱历时、烈度和烈度峰值只能由对数正态分布拟合。汾河水库的干旱历时和烈度可由对数正态分布、指数分布和广义帕累托分布拟合，其中对数正态分布对二者的拟合效果最好；烈度峰值只可由对数正态分布拟合。寨上的干旱历时和烈度可由对数正态分布、指数分布和广义帕累托分布拟合，其中历时由对数正态分布拟合效果最好；烈度由广义帕累托分布拟合效果最好；烈度峰值只可由对数正态分布拟合。兰村的干旱历时和烈度可由对数正态分布和指数分布拟合，指数分布对二者的拟合效果最好；烈度峰值只可由对数正态分布拟合。

利用矩法对各站的干旱特征变量边缘分布进行参数估计，结果见表3。

表3　干旱特征变量边缘分布及其参数

2.2干旱特征变量的二维联合分布

选取Clayton Copula、GH Copula、AMH Copula 和Nelson No.2 Copula 4种Copula函数建立干旱特征变量的二维联合分布。根据Kendall秩相关系数τ 和Archimedean Copula函数的参数θ的解析关系确定参数θ，再根据式（6）进行拟合优度评价，依据是RMSE与AIC值越小，拟合效果越好（见表4）。

由表4可看出，4个站的干旱历时、烈度与烈度峰值两两之间的Kendall秩相关系数τ分别介于0.815～0.889、0.583～0.663、0.687～0.729之间，表明干旱特征变量间具有明显的相关关系，可用Copula函数构建其联合分布。上静游站，Clayton Copula对于干旱历时、烈度和烈度峰值两两之间的拟合的效果都最好；汾河水库站，GH Copula对于干旱历时和烈度之间拟合效果最好，Clayton Copula对于干旱历时和烈度峰值、烈度和烈度峰值之间的拟合的效果最好；寨上站，Nelson NO.2 Copula对于干旱历时和烈度之间拟合效果最好，Clayton Copula对于干旱历时和烈度峰值之间拟合效果最好，GH Copula对于烈度和烈度峰值之间的拟合的效果最好；兰村站，GH Copula对于干旱历时、烈度和烈度峰值两两之间的拟合的效果都最好。由于篇幅所限，以寨上站为例，将优选的Copula函数的经验频率和理论频率对比如图2。

表4　二维干旱特征变量的参数及拟合优度评价

由图2可看出，横坐标表示经验频率，纵坐标表示理论频率的点据均匀分布在横纵坐标轴的对角线上及两侧，表明经验频率与理论频率近似相等，说明Copula函数对于寨上的干旱特征变量的二维分布拟合效果较好。

图2　寨上站的干旱历时与干旱烈度、干旱历时与烈度峰值、干旱烈度与烈度峰值拟合

2.3重现期计算

一组数据中最具代表性且最能体现数据特征的5个数是中位数Me，下四分位数Q1，上四分位数Q3，最小值Min和最大值max。因此，选取各站提取的各干旱特征变量序列中的这5个代表性数据，分别计算其单变量重现期，再计算二维Copula函数值，求得对应的联合重现期（Tα）与同现重现期（T0），结果见表5。

表5　基于最优Copula函数的水文干旱特征变量的重现期

从表5中的单变量重现期可看出：①4个站的干旱特征变量序列的最小值min、下四分位数Q1、中位数me和上四分位数Q3对应的重现期分别为0.74～1.12年、0.9～1.49年、1.37～2.86年和2.18～5.10年。说明这4个站在调查年份中所发生的干旱大多数重现期小于5年。②4个站的干旱特征变量序列的最大值max对应的重现期变化范围很大，其中烈度峰值最大值对应的重现期较小，从4.71～17.15 a，干旱历时与烈度最大值对应的重现期较大，分别为 73.73～1045.4 a和 72.95～990.09 a。说明所发生过的干旱事件中，烈度峰值对应重现期比干旱历时与烈度对应的重现期小很多。③从干旱历时与烈度看，兰村站在调查年分内干旱最严重，已发生过超越或接近千年一遇的大旱。

从表5的三种重现期对比可看出：①4个站的二维联合重现期小于等于其对应的单变量重现期；单变量重现期小于或等于与其对应的二维同现重现期。②由于4个站干旱历时与烈度的相关性很强，二者的三种重现期都较为接近；干旱历时与烈度峰值、烈度与烈度峰值的相关性较弱，它们的三种重现期差距较大。这给今后汾河上游的水文干旱的分析预测提供了一定依据。

3　结论

（1）对数正态分布对于汾河上游的干旱特征变量的边缘分布拟合效果最佳，指数分布和广义帕累托分布次之，伽马分布不适合汾河上游的干旱特征变量的边缘分布拟合。

（2）汾河上游的4个水文站的单变量重现期介于对应的二维联合重现期与同现重现期之间。

（3）汾河上游所发生的干旱大多数重现期较短，且所发生的极少数严重干旱的特点是历时长、烈度大、峰值并不太高。

（4）汾河上游干旱历时与烈度的三种重现期都较为接近；干旱历时与烈度峰值、烈度与烈度峰值的三种重现期差距较大。

［1］SHIAU J T，MODARRES R.Copula-based drought severity-durationfrequency analysis in Iran［J］.Meteorological applications，2009，16 （4）：481-489.

［2］张强，潘学标，马柱国，等.干旱［M］.北京：气象出版社，2009：30.

［3］BACCHI B，BECCIU G，KOTTEGODA N T.Bivariate exponential model applied to intensities and durations of extreme rainfall［J］. Journal of Hydrology，1994，55（1/2）：225-236.

［4］Yue S，Ouarda T B M J，Bobee B.A review of bivariate gamma distributions for hydrological application［J］.Journal of Hydrology，2001，246（1）：1-18.

［5］SHIAU J T.Return period of bivariate distributed hydrological events ［J］.Stochastic Environmental Research and Risk Assessment，2003，17（1/2）：42-57.

［6］YUE S.Applying bivariate normal distribution to flood frequency analysis［J］.Water International，1999，124（3）：248-254.

［7］BOX G，COX D.An analysis of transformations［J］.Journal of the Royal Statistical Society，1964，26（2）：211-252.

［8］戴昌军，梁忠民.多维联合分布计算方法及其在水文中的应用［J］.水利学报，2006，37（2）：160-165.

［9］王文圣，丁晶.基于核估计的多变量非参数随机模型初步研究［J］.水利学报，2003，34（2）：9-14.

［10］KIM T W，VALDES J B，YOO C.Nonparametric approach for estimating return periods of droughts in arid regions［J］.Journal of Hydrologic Engineering，2003，8（5）：237-246.

［11］刘和昌，梁忠民，姚轶，等.基于Copula函数的水文变量条件组合分析［J］.水力发电，2014，40（5）：13-16.

［12］YEVJEVICH V M.An objective approach to definitions and investigations of continental hydrologic droughts［M］.Colorado State University，1967.

［13］TALLAKSEN L M，MADSEN H，CLAUSEN B.On the definition and modelling of streamflow drought duration and deficit volume［J］. Hydrological Sciences Journal，1997，42（1）：15-33.

［14］宋松柏，蔡焕杰，金菊良，等.Copulas函数及其在水文中的应用［M］.北京：科学出版社，2012.

［15］NELSON R B.An Introduction to Copulas［M］.New York:Springer，2006.

［16］SHIAU J T，SHEN H W.Recurrence analysis of hydrologic droughts of differing severity［J］.Journal of Water Resources Planning and Management，2001，127（1）：30-40.

（责任编辑陈萍）

锦屏一级水电站左岸基础处理工程通过验收

锦屏一级水电站左岸基础处理工程顺利通过完工验收。电站左岸基础处理工程的成功建设，较好地解决了复杂地质条件制约特高拱坝建设的关键技术难题，标志着我院复杂地基处理技术水平处于世界领先水平。

锦屏一级水电站拱坝高305m，是世界最高拱坝。基础工程地质条件极为复杂，左岸抗力体发育有断层、煌斑岩脉、深部裂缝及低波速岩带、层间挤压错动带等，这些软弱结构面规模大、性状差，对坝体受力状态和变形稳定产生较大影响。如此复杂的地质条件，极大地增加了拱坝设计和基础处理设计难度，其设计技术与施工超出了现行设计规范和已有工程经验。

左岸基础处理工程采取了防渗帷幕灌浆、坝基排水及抗力体排水、左岸混凝土垫座、左岸抗力体固结灌浆、断层及煌斑岩脉混凝土网格置换、基础洞室回填、软弱岩带水泥-化学复合灌浆等综合加固处理措施，处理措施的规模与难度在水电界少有。其经验和成果对促进复杂地基条件下高拱坝的建设，具有较高的借鉴价值和指导意义。

（中国电建集团成都勘测设计研究院有限公司）

Copula-based Analysis of Hydrologic Drought Frequency in the Upper Reaches of Fenhe River

REN Lu，ZHAO Xuehua
（College of Water Resources Science and Engineering，Taiyuan University of Technology，Taiyuan 030024，Shanxi，China）

Monthly runoff data taking from four hydrometric stations in the upper reaches of Fenhe River are applied to identify the droughts by using Run Theory.The drought characteristic variables are fitted by four distribution functions，and the Kolmogorov-Smirnov method is employed to evaluate the goodness-of-fit of univariate marginal distributions.The Copula functions are employed to construct two-dimensional joint distribution models of drought characteristic variables，and then their return periods are determined.The results show that:（a）the lognormal distribution is the candidate marginal distribution function with highest goodness-of-fit in fitting drought characteristic variables in the upper reaches of Fenhe River；（b）the return periods of most of the droughts are less than 5 years，and the univariate return periods and the two-dimensional return periods of drought duration and severity are close to each other；and（c）there are large differences between the univariate return periods and the two-dimensional return periods of drought duration and peak intensity，as well as the drought severity and peak intensity.

hydrologic drought；frequency analysis；Copula function；monthly runoff；upper reaches of Fenhe River

P333.3

A

0559-9342（2016）02-0011-06

2015-03-09

国家自然科学基金资助项目（40901018）；山西省科技攻关项目（20140313023-4）；太原理工大学校青年团队启动项目（2013T039）

任璐（1991—），女，山西运城人，硕士研究生，研究方向为水资源系统工程；赵雪花（通讯作者）.