图1 水文干旱的识别过程
马建琴1,2，和鹏飞1，彭高辉1，崔弼峰1，刘 蕾1
（1.华北水利水电大学，郑州 450045；2.水资源高效利用与保障工程河南省协同创新中心，郑州450045）
摘 要：为进一步研究水文干旱特征及其演变规律，选取沙颍河流域2个水文站1951—2008年逐月径流资料，采用游程理论提取水文干旱特征变量，选取4种对称Archimedean Copula函数和对应的非对称Archimedean Copula函数拟合干旱特征变量的联合分布并计算联合重现期和同现重现期。经拟合检验，Frank Copula对干旱特征变量二维联合分布的拟合效果最优，M3 Copula函数对三维干旱变量的联合分布拟合效果最优。干旱特征的单变量重现期值大于二、三维变量的联合重现期，小于二、三维变量的同现重现期。因此，Copula函数能够较好地拟合沙颍河流域水文干旱特征变量间的联合分布。
关键词：沙颍河流域；标准化径流指数；游程理论；Copula函数；重现期
干旱是由于持续时间段内降水缺少导致水分收支或供求失衡形成的水资源亏缺[1]，且造成的灾害对区域经济和社会发展危害较大，干旱的形成机制和演变规律研究显得十分重要[2-3]。Copula函数能在不考虑单变量边缘分布类型情况下构建多变量间的联合分布，相对传统的多变量分析方法更有优势[4-5]，在水文、气象领域中得到广泛应用[6]。国内外学者对二维联合分布的研究较多，诸如提取了极端降水事件的特征变量，运用二维Copula函数构建了雨量、历时与平均强度的联合分布模型[7]；以汉江上游干旱事件为例，用Copula函数研究了干旱历时与干旱烈度的联合概率分布[8]；运用椭圆型Copulas函数模型分析了西安站干旱特征[9]；利用非对称Archimedean Copula函数研究洪水频率，分析了变量间相关性较低时的联合分布模型[10]。随着理论的发展，高维变量间的联合分布研究也随之增多，如应用5种非对称Archimedean Copula函数对径流序列进行了多变量的水文干旱分析，并计算了3维变量的概率与重现期[11]；基于Pair-copula函数对南盘江流域的水文干旱频率分析并计算了不同站点的联合重现期[12]。表明Copula函数对二、三维干旱特征变量间联合分布的拟合有较好的适用性。然而，现有研究运用特定水文干旱指数进行干旱频率分析的不多。标准化径流指数能够准确地分析流域水文干旱情况，具有较强适用性[13]，故选择标准化径流指数作为识别水文干旱的指标。
沙颍河发源于河南省伏牛山区，至安徽正阳汇入淮河，全长624 km，流域面积39 880 km2。沙颍河流域受南北过渡气候的影响，流域降水时空差异较大。由于干旱事件的影响，制约了流域农业生产和社会经济的发展。选取位于沙颍河流域河南段上游和下游的漯河与周口水文站，基于二站1951—2008年的逐月径流资料计算标准化径流指数，依据游程理论提取水文干旱特征变量，用Copula函数建立干旱变量联合分布模型，计算对应的各个变量联合重现期和同现重现期，以此来推求沙颍河流域的水文干旱演变特征，为流域农业生产、水资源规划、内河航运和防旱抗旱提供一定的理论支撑。
SRI[13]不仅可以描述多时间尺度的干旱分析，在资料缺乏和地形复杂的流域适用性也较好，计算方法比较简单。计算原理和干旱等级划分标准见文献[14-15]。
根据游程理论[16]提取水文干旱特征变量，干旱历时D为一次干旱事件持续的时间；干旱烈度S是干旱历时时段内SRI的累加，即阴影部分的面积；干旱烈度峰值P为一次干旱事件中最小SRI，即负游程的极值。
识别的过程：选取SRI阈值R1=-1、R0=-0.5、R2=0（图1），若SRI小于R0，初步判断此月为干旱，即时段a、b、d、e；当干旱历时为1个时段(a、e)，若SRI小于R1，则为1次干旱事件(a)，反之不计为干旱事件(e)；若相邻2次干旱事件(b、d)间隔期的SRI2小于R2，即Sc＜0,则合并为1次干旱事件，干旱历时D=Db+ Dd+1，干旱烈度（阴影面积）S=Sb+Sd，否则为2次干旱事件。
图1 水文干旱的识别过程
选择二参数函数—伽马(Gamma)函数、指数(Exponential)函数、对数正态(log-normal)函数、威布尔(Weibull)分布函数分别拟合沙颍河流域2个水文站的干旱历时、烈度及烈度峰值的边缘分布，利用极大似然估计参数，用Kolmogorov-Sm irnov(K-S)检验法评价拟合度，计算原理详见文献[17]。K-S检验统计量越低，表示变量的理论分布和经验分布拟合度越高，以此选出最佳分布函数。
Copula是连接单变量边缘的分布，形成在[0，1]区间上服从均匀分布的多元函数[18]。设F1，F2，…,Fn为连续随机变量，则有唯一Copula函数C对任意x∈Rn满足[19]：
式中：x1,x2,…,xn为观测样本；F(x)为边缘分布函数。
选取Frank、Clayton、Gumbel-Hougaard(Gumbel)、Joe Copula 4种对称的Archimedean Copula函数及4种非对称Archimedean Copula函数M3、M 4、M 5、M 12 Copula，拟合三维干旱变量的联合分布。采用半参数极大估计法[17]对参数θ进行估计。采用均方根误差(RMSE)准则和赤池信息量准则(AIC)进行拟合优度检验，RMSE和AIC最小则拟合结果最优，其计算公式为：
式中：n为样本容量；k为Copula参数个数；i为样本序号；Pi为第i个Copula函数联合分布计算值；Pei为第i个经验值；MSE为均方误差。
根据重现期理论[20]，可得出单变量重现期的计算公式，即：
式中：F(x)为干旱特征变量的分布函数；N为站点资料的时间长度；m为干旱事件发生的次数。
二维及三维干旱特征变量重现期分为联合重现期(Ta)和同现重现期(To)，其计算原理见文献[21]。
以沙颍河流域漯河水文站、周口水文站1951—2008年的月径流资料为例，用游程理论进行干旱识别，提取干旱特征变量，并用伽马(Gamma)、指数(Exponential)、对数正态(log-normal)、威布尔(Weibull)分布分别拟合2个水文站的干旱历时、烈度及烈度峰值，选用极大似然法进行参数估计，结果见表1。利用K-S检验法对各变量边缘分布进行拟合优度检验，并用Kendall秩相关系数表示变量间的相依程度，结果见表2。
表1 单变量边缘分布及其参数估计
表2 边缘分布函数的拟合检验统计量及相关系数
从表2可看出，在显著性水平α=0.05下，2个站点干旱变量参数的检验统计量均比临界值小，p均大于0.05，干旱特征变量分布函数拟合较好。干旱特征变量间的秩相关系数τ均大于0.5，说明干旱特征变量间存在较高的相关性，可以用Copula函数构建其联合分布。
选取4种Copula函数，参数θ用半参数估计法计算，再用式(2)评价拟合优度，结果见表3。
表3 二维干旱变量参数值及拟合度评价
由表3可知，根据RMSE和AIC评判准则，漯河站Frank Copula对干旱历时、烈度间的拟合效果较好，Gumbel Copula对干旱历时与烈度峰值的拟合效果较好，ClaytonCopula拟合烈度、烈度峰值联合分布的效果最优；周口站Frank Copula对二维联合分布拟合效果较好。以周口站为例，将最优Copula函数（Frank Copula）的经验频率及理论频率进行对比，结果见图2。
由图2可知，Frank Copula函数经验频率与理论频率点均匀分布在45°对角线附近，拟合效果较好，可以作为模拟二维干旱特征变量联合分布的连接函数。
图2 周口站干旱特征变量二维联合分布拟合结果
选取4种对称Archimedean Copula函数和对应构造的非对称Archimedean Copula函数拟合干旱变量的三维联合分布，Copula函数参数θ用半参数估计法计算，再用式（2）评价拟合度，结果见表4。
表4 三维干旱变量参数值及拟合度评价
由表4可知，根据RMSE和AIC评判准则，漯河站Frank Copula对干旱历时、烈度与烈度峰值三维联合分布拟合的效果最优；周口站M3Copula函数对干旱历时、烈度与烈度峰值三维联合分布拟合的效果最优。在干旱频率分析中，通常要考虑干旱特征变量的尾部相关性，以确定某个干旱特征变量的取值大小是否会对其他变量的取值产生影响，从而提高干旱风险评估的精确性[22]。
将漯河站与周口站的最优Copula函数的经验频率及理论频率进行对比可知（图3），Frank Copula、M3 Copula函数点均匀在45°对角线附近，M3 Copula函数的上尾相关性较好，即Frank Copula、M3 Copula函数可以作为模拟三维干旱特征变量联合分布的连接函数。
图3 沙颍河流域干旱特征变量三维联合分布拟合结果
取单变量重现期（T）水平为2、5、10、20和50 a，分别求出各干旱特征的单变量值，根据文献[21]计算对应二、三维Copula函数的值以及给定单变量重现期水平下对应的二、三维联合（Ta）重现期和同现（To）重现期的值，结果见表5。
表5 周口站干旱特征变量的重现期
由表5可知，同一单变量重现期（T）水平下，二维变量联合分布中干旱历时D与烈度S的联合重现期（Ta）最大，对应同现重现期（To）最小，说明二者之间相关程度最高。干旱历时D与烈度峰值P的联合重现期（Ta）最小且对应同现重现期（To）最大，说明二者之间相关程度最低。干旱特征变量两两间相关程度越高，干旱事件中某个特征变量发生时，另一个特征变量同时发生的概率越大，干旱风险性就越高。根据变量间相关系数τ可知，沙颍河周口站的干旱风险高于漯河站。
单变量重现期介于对应的二、三维干旱变量的联合重现期（Ta）及同现重现期（To），说明联合重现期（Ta）与同现重现期（To）可以作为单变量重现期（T）的上下阈值，以此可推求出干旱特征联合概率，有助于干旱预测和旱情分析。从各重现期水平来看，周口站干旱历时小于漯河站，而干旱烈度、干旱烈度峰值大于漯河站，说明干旱事件发生时下游的旱情比上游更严重。
研究选取的干旱指标与传统水文干旱指标降水距平百分率和径流距平百分率不同，选取标准化径流指数作为干旱研究和识别的指标，计算的稳定性和适应性比前2种指标相对较好，能准确地反映整个区域内的干旱情况[13]。Copula函数在干旱研究中的应用多集中在建立干旱历时和干旱烈度的二维联合分布模型，但干旱特征变量还包括干旱烈度峰值及其他多个相关变量，二维Copula函数模型已不能全面反映干旱的特征，研究构建了水文干旱特征变量干旱历时、干旱烈度、干旱烈度峰值三维联合分布，结合多个特征变量综合分析水文干旱特征，保证了研究结果的全面性；在三维Copula函数模型参数选择和计算中，单参数的对称Archimedean Copula模型对干旱特征多变量间的相依性描述有一定的局限性[11]，对多变量间非对称相依结构描述不足，因此，选择有更多参数的非对称Archimedean Copula函数构建干旱特征变量的三维模型十分有必要，本研究的结果也表明非对称Archimedean Copula函数构建干旱特征变量的三维模型中取得了更好地拟合效果。
目前已有的Copula函数有上百种，只有少数在研究应用中比较成熟，而应用于干旱研究的更少，其他类型的Copula函数将是以后研究的重点，随着干旱研究的变量增多，高维度Copula函数将成为未来新的研究方向，同时也要注意到随着Copula函数维度的增高，参数估计以及计算难度都会相应地增加。
1）二维干旱变量联合分布拟合函数中，Frank Copula函数的RMSE、AIC值最小，函数经验频率与理论频率点在45°对角线附近分布均匀，且上尾相关性好，适合模拟沙颍河流域水文干旱特征变量的二维联合分。
2）M3Copula函数对干旱历时、烈度与烈度峰值三维联合分布拟合的效果最优。将漯河站与周口站的最优Copula函数的经验频率及理论频率进行对比可知（图3），Frank Copula函数中部和上尾的相关性较好；M3Copula函数显示出下尾和上尾较好的相关性
3）沙颍河2个水文站的单变量重现期均介于对应的联合重现期(Ta)和同现重现期(To)之间，联合重现期(Ta)和同现重现期(To)可以作为单变量重现期(T)的上下阈值，用于干旱预测和旱情分析。
4）上游漯河站干旱特征变量间的相关系数小于下游周口站，相同重现期水平内上游的干旱历时比下游长，干旱烈度及烈度峰值小于下游，下游的干旱风险较高，且旱情比上游更严重。
参考文献：
[1]董前进,谢平.水文干旱研究进展[J].水文,2014,34(4):1-7.
[2]李天水,王顺,庄文化,等.游程理论和Copula函数在二维干旱变量联合分布中的应用[J].干旱区资源与环境,2016,30(6):77-82.
[3]耿鸿江,沈必成.水文干旱的定义及其意义[J].干旱地区农业研究,1992(4):95-98.
[4]李计,李毅,贺缠生.基于Copula函数的黑河流域干旱频率分析[J].西北农林科技大学学报（自然科学版），2013,41(1):213-220.
[5]谢华,黄介生.两变量水文频率分布模型研究述评[J].水科学进展,2008,19(3):4-4.
[6]郭生练,闫宝伟,肖义,等.Copula函数在多变量水文分析计算中的应用及研究进展[J].水文,2008,28(3):1-7.
[7]ZHANG L,SINGH V P.Bivariate rainfall frequency distributionsusing Archimedean copulas[J].Journalof Hydrology,2007,332(1/2):93-109.
[8]闫宝伟,郭生练,陈璐,等.长江和清江洪水遭遇风险分析[J].水利学报,2010,39(5):553-559.
[9]马明卫,宋松柏.椭圆型Copulas函数在西安站干旱特征分析中的应用[J].水文,2010,30(4):36-42.
[10]GRIMALDIS,SERINALDIF.Asymmetric copula inmultivariate flood frequency analysis[J].Advances inWater Resources,2006,29(8):1155-1167.
[11]宋松柏,聂荣.基于非对称阿基米德Copula的多变量水文干旱联合概率研究[J].水力发电学报,2011,30(4):20-29.
[12]杨茂灵,王龙,余航,等.基于Pair-copula函数和标准化径流指数的水文干旱频率分析——以南盘江流域为例[J].长江流域资源与环境,2014,23 (9):1 315-1 321.
[13] SHUKLA S,WOOD AW.Use of a standardized runoff index for characterizing hydrologic drought[J].Geophysical Research Letters,2008,35(2): 226-236.
[14]黄晚华,杨晓光,李茂松,等.基于标准化降水指数的中国南方季节性干旱近58a演变特征[J].农业工程学报,2010,26(7):50-59.
[15]GB/T20481-2006，气象干旱等级[S].北京:中国标准出版社，2006.
[16]YEVJEVICH V.An objective approach to definitions and investigations of continental hydrologic droughts[M].Fort Collins,Colorado:Colorado State University,1967.
[17]宋松柏.Copulas函数及其在水文中的应用[M].北京：科学出版社,2012.
[18]张翔,冉啟香,夏军,等.基于Copula函数的水量水质联合分布函数[J].水利学报,2011,42(4):483-489.
[19]NELSEN RB.An Introduction to Copulas[M].New York:Springer,1998.
[20]SHIAU JT,SHEN HW.Recurrence Analysis of Hydrologic Droughts of Differing Severity[J].Journal of Water Resources Planning&Management, 2001,127(1):30-40.
[21]李计,李毅,宋松柏,等.基于Copulas函数的二维干旱变量联合分布[J].自然资源学报,2013,32(2):312-320.
[22]刘占明,陈子遷,黄强,等.基于Copula函数和SPEI的广东北江流域干旱频率分析[J].中山大学学报（自然科学版）,2014,53(5):118-125.
Analyzing the Drought Frequency in Shaying River BasinUsing the Three-dimensional Copula Function
MA Jianqin1,2,HEPengfei1,PENGGaohui1,CUIBifeng1,LIU Lei1
(1.North China University ofWater Resourcesand Electric Power，Zhengzhou 450045,China；2.Collaborative Innovation CenterofWater Resources EfficientUtilization and Guarantee Engineering,Henan Province，Zhengzhou 450045,China）
Abstract:This paper analyzed the drought characteristics based on datameasured in 1951—2008 from two hydrologicalstations in the Shaying River Basin.We took standardized runoff index as an indicator of the hydrologicaldrought,calculating the hydrological drought characteristics including droughtduration,droughtseverity,and peak of the droughtseverity,using the theory of run test.The jointdistributionsof the drought characteristic variableswere analyzed using four types of symmetric and asymmetric Archimedean Copulas functions;the goodnessof the fitting of the functionswasevaluated using the Kolmogorov-Sm irnov testmethod.The results showed that the Frank Copula function gave the best fitting for two-dimensional joint distributions of two characteristic droughtvariables,while the M3 copula function worked better for fitting three-dimensional joint distributions of three characteristic droughtvariables.Itwasalso found that the return period calculated from single variablewas greater than that calculated from the two-and three-dimensional distribution for the joint return periods,but less than that from the two-and three-dimensional distributions for the co-occurrence return periods.Thus,the Copula functions are able to describe themultivariate jointdistributions of hydrological droughtvariables in the Shaying River Basin.
Key words:Shaying River Basin;Standardized Runoff Index;theory of run;Copulas function;return period
中图分类号：S276
文献标志码：A
doi：10.13522/j.cnki.ggps.2017.09.018
责任编辑：刘春成
马建琴，和鹏飞，彭高辉，等.基于三维Copula函数的沙颍河流域水文干旱频率分析[J].灌溉排水学报，2017，36（9）：102-107.
文章编号：1672-3317（2017）09-0102-06
收稿日期：2016-11-29
基金项目：河南省高校科技创新人才支持计划项目（15HASTIT046）；河南省科技攻关项目（152102110095）；河南省高等学校重点科研项目（15A570008）；国家自然科学基金项目（41601019）
作者简介：马建琴（1973-），女。教授，研究方向为干旱地区水资源管理、农业水资源可持续利用、区域水资源优化配置、水环境。E-mail：393204148@qq.com
通信作者：和鹏飞（1987-），男。硕士研究生，研究方向为水文水资源。E-mail：hexu_19870831@163.com