不同重现期标准双变量设计洪水计算方法

晋 恬闻 昕方国华黄显峰袁 玉

(1.河海大学水利水电学院,江苏南京 210098;2.中国水利水电科学研究院流域水循环模拟与调控国家重点实验室,北京 100038)

设计洪水是我国水利工程在规划、设计、运行、管理等阶段的重要防洪依据,常通过同倍比、同频率放大典型洪水的方法得到。同倍比法仅考虑单个洪水特征变量因素,同频率法虽然考虑了洪峰流量和洪量两个因素,但未考虑二者之间的联系,对于感潮河流等洪水成因复杂的流域存在较大的防洪隐患和风险。双变量设计洪水则考虑了洪水成因的复杂性以及洪水事件本身的多属性特征,可为工程和区域防洪安全提供新的视角和科学依据。目前,对于双变量设计洪水的计算方法已有较多探索和研究,但这些研究多针对单一重现期标准或单一计算方法,不同重现期标准和不同计算方法,特别是不同双变量重现期标准适用的设计洪水值计算方法尚不明晰。

Copula函数[1]是研究变量相关性测度的工具,可以刻画变量之间的非线性关系,具有边缘分布灵活以及其形式不受边缘分布限制的特点,可用于洪水、暴雨、干旱等多变量水文事件的联合概率分布分析[2鄄4]。 Archimedean Copula 函数由于具有构造简单、仅含有一个参数的优点,水文事件联合概率分布多选用该类函数进行拟合分析。目前常用于描述双变量洪水频率的重现期标准有联合重现期和同现重现期两种,Salvadori等[5鄄7]根据描述危险区域的不同提出了第二重现期,并将其与联合重现期进行了比较分析,确定两者有较强的指数相关关系;史黎翔等[8]进一步研究了上述3种重现期标准的差异,论证了第二重现期在危险区域描述方面的合理性。目前确定设计洪水值的方法[9鄄19]大致分为两类:第一类为基于重现期等值线确定设计洪水值的方法,常用的有最可能超越法、极大似然法、双变量同频率法;另一类为先确定满足条件的洪水峰量组合线,再求出该曲线与洪水重现期等值线交点(即为对应重现期的设计洪水值)的方法,常用的有条件最可能组合法、条件期望组合法、双变量同频率组合法。

本文基于Archimedean Copula函数建立峰量组合二维分布模型,选择C测度、联合、同现重现期,采用极大似然法、条件最可能组合法和双变量同频率法,研究不同重现期标准下双变量设计洪水的计算方法。

1 基于Copula函数的洪水峰量联合分布模型

洪水具有多特征变量的属性,如洪峰流量、洪量、历时、洪水过程线等,本文选择洪峰流量、洪量两个变量构建洪水联合分布模型。

定义q、w分别为一场洪水过程的洪峰流量和洪量,设FQ(q)、FW(w)为洪峰流量、洪量的边缘分布函数,其联合分布函数为F(q,w),按照Sklar定理[1],洪峰流量、洪量的联合分布函数可以用对称Archimedean Copula函数来表示:

其中

式中:兹为 Copula函数的参数;u、v分别为洪峰流量、洪量的边缘分布函数。

2 双变量重现期标准计算方法

目前描述双变量的常见重现期标准为联合重现期和同现重现期。不同重现期标准的差异本质为定义危险事件方式的不同。

若以洪峰流量和洪量两个特征量来表述洪水事件,对于联合重现期T胰,其危险事件E胰表示洪峰流量或洪量其中一个变量超过设定的阈值(q0或w0),对应的危险域AE胰={(q,w)沂R:q>q0胰w>w0}。联合重现期可用Copula函数表示为

对于同现重现期T疑,其危险事件E疑表示洪峰流量和洪量同时超过设定的阈值,对应的危险域为同现重现期可用Copula函数表示为 假设一个洪水事件的双变量重现期为T,联合重现期和同现重现期的危险区域如图1所示,图中的实线即为重现期等值线。图1(a)中玉区域为在联合重现期等值线上任选一点A(u1,v1)所对应的危险区域,图1(b)中域区域为在同现重现期等值线上任选一点B(u2,v2)所对应的危险区域,通过公式u=FQ(q)、v=FW(w)可以推求出等值线上不同点所对应的洪峰流量q以及洪量w。结合图1(a)和(b)进行分析,在重现期等值线上选择不同的点会得到相对应的危险区域,选定的点不同,对应的危险区域也不同,对于上述两种重现期标准,其重现期水平与危险区域并非一一对应关系。另外,联合、同现重现期对危险区域的定义也存在不合理性,一般而言,重现期越大相应的危险区域越小,重现期大的事件的危险区域应包含于重现期小的事件中,但在联合、同现重现期标准下存在不满足上述规律的情况,结合图2可看出,B点处的重现期大于A点处的重现期(TB>TA),B点处的危险区域包括事件E,即事件E在B点处被识别为危险事件,但是在A点,其危险区域并不包括危险事件E。采用上述两种重现期标准,均出现了重现期较大的洪水组合识别的危险事件在重现期较小的洪水组合中被识别为非危险事件的矛盾的情景。

图1 重现期的危险区域

图2 重现期危险区域识别分析

基于上述两种重现期标准的局限性,将危险事件定义为洪水峰量组合在重现期等值包络线以上的区域,即图1(c)中的芋区域。假设重现期等值线满足C兹(u,v)=t,其对应的危险区域为AEC(t)={(u,v)沂玉2:C兹(u,v)>t},其中0t)可以用C测度函数KC(t)[13]来表示:

式中(t)为Copula函数的生成元。

在C测度函数下对应的重现期为C测度重现期TC,表达式为 C测度重现期等值线上任意一点对应的危险区域均相同,且重现期水平越大,危险区域越小,重现期水平小的危险区域必定会覆盖重现期水平大的危险区域。可见C测度重现期对危险区域的定义更为合理。

C测度重现期是以联合重现期等值线为基础通过修正危险区域得到的新的重现期标准,因此其等值线线型与联合重现期类似。C测度函数是以Archimedean Copula函数为前提条件推导得到的,在此基础上定义的C测度重现期仅适用于以Archimedean Copula函数构建的联合分布模型,对于其他类型的Copula函数的适用性还需论证。

3 设计洪水值计算

对于选定的洪水双变量重现期水平,存在一系列峰量组合满足防洪标准要求,需根据一定的准则从中选择设计洪水值组合,本文选择极大似然法、条件最可能组合法以及双变量同频率法作为推求设计洪水值的方法,并以Archimedean Copula函数中的G鄄H Copula函数为例进行分析。

3.1 极大似然法

在相应重现期等值线上,要求设计洪水值组合(qd,wd)满足权重值GML最大,其目标函数为

式中:c兹(u,v)为 Copula函数的联合密度函数;fQ(q)、fW(w)分别为洪峰流量、洪量对应的边缘密度函数。

式(6)应满足重现期等值线约束,即保证所求设计洪水值组合(qd,wd)均在对应重现期等值线上,约束条件为

3.2 条件最可能组合法

当洪峰流量q取qi时,要求洪量w对应的密度函数值fW|Q(w)最大,其目标函数为

求出洪峰流量qi对应的洪量wi,其对应的洪水峰量组合为(qi,wi),得出一系列的峰量组合点连成条件最可能洪水峰量组合线,其与对应重现期等值线的交点即为设计洪水值组合(qd,wd)。

式(8)以洪峰流量为条件进行分析计算。

3.3 双变量同频率法

采用双变量同频率法求解设计洪水值组合,即要在相应重现期T上求u=v的点为设计洪水值组合(qd,wd),应满足以下条件:

也可先求出双变量同频率洪水降量组合线,该线洪峰流量与洪量的边缘分布函数相同,其目标函数为

双变量同频率洪水峰量组合线与对应重现期等值线的交点即为设计洪水值组合(qd,wd)。

4 算例分析

新安江水库建在钱塘江支流新安江铜官峡,位于建德市新安江镇以北约4郾5 km,坝址以上河长为323 km,坝址控制流域面积10 442 km2,主要承担发电、防洪、供水等任务[20鄄21]。罗桐埠站是新安江水库的设计依据站,本文选取罗桐埠站1961—2008年的洪水资料,根据新安江流域的洪水特征,提取历时为48 a的年最大洪峰流量q和年最大10 d洪量w10的峰量数据,进行峰量联合设计洪水研究。

4.1 边缘分布

首先分别建立年最大洪峰流量q和年最大10 d洪量w10的边缘分布函数。GEV分布被广泛用于洪峰流量、洪量的概率统计分析拟合中并取得了较好的拟合效果,本文选择GEV概率分布函数作为边缘分布函数,采用极大似然法估计 FQ(q)和FW(w)的统计参数,运用K鄄S拟合优度检验方法对边缘分布进行检验,计算q和w10的样本检验统计值和显著水平琢=0郾05的检验临界值,结果见表1。对于q和w10,样本检验统计值均小于相应的检验临界值,表明FQ(q)和FW(w)均服从GEV分布。

表1 洪水特征量边缘分布统计参数

4.2 联合分布

选择Archimedean Copula函数中的G鄄H Copula函数、Clayton Copula函数和 Frank Copula函数对q鄄w10联合分布进行模拟,采用极大似然估计法计算参数兹,并采用赤池信息量准则(Akaike information criterion,AIC)、贝叶斯信息量准则(Bayesian information criterion,BIC)和均方根误差(RMSE)3种拟合优度检验方法选择最优Copula函数,AIC、BIC、RMSE值越小代表Copula函数的拟合度越好。计算结果见表2,其中G鄄H Copula的AIC、BIC值均最小,RMSE值略大于Frank Copula函数,基于拟合结果选择G鄄H Copula构建q鄄w10联合分布。图3为经验联合分布频率与理论联合分布(G鄄H Copula函数)频率拟合结果,样本数据均在45毅线附近,表明G鄄H Copula函数拟合度较好。

表2 q鄄w10联合分布Copula函数相关参数

图3 经验联合分布频率与理论联合分布频率拟合结果

4.3 不同重现期标准比较分析

根据同现重现期、联合重现期、C测度重现期计算公式,绘制不同重现期标准下新安江流域10年、20年、50年、100年、200年一遇洪水的重现期等值线如图4所示。由图4分析可知,对相同重现期标准,不同重现期水平的等值线走势一致;对不同重现期标准,联合重现期等值线与C测度重现期等值线走势相同,联合、C测度重现期与同现重现期等值线走势相反。从危险事件定义的角度分析,同现重现期与联合重现期等值线的走势差异性体现在前者考虑的为并集危险事件,后者考虑的为交集危险事件;而C测度重现期是基于联合重现期等值线提出的修正,因此二者走势相似。

图4 不同重现期标准下重现期等值线及设计洪水值组合

图5 不同重现期水平下重现期等值线及设计洪水值组合

在给定重现期水平下,不同重现期标准之间的关系见图5。同现、联合和C测度重现期在同一重现期水平下其重现期等值线无交点;对于联合重现期与同现重现期,其重现期等值线的渐近线均为单变量(q或w10)在相同重现期水平下的设计洪水值;对于C测度重现期,其重现期渐近线均小于单变量(q或w10)在相同重现期水平下的设计洪水值。对于同一洪峰流量值(洪量值),联合重现期等值线对应的洪量值(洪峰流量值)最大,C测度重现期等值线对应的次之,同现重现期等值线对应的最小。这是由于同现重现期标准下的峰量组合值仅考虑了并集危险事件发生的概率,使得其峰量组合值均偏小;而联合重现期标准下的峰量组合值却相对偏大,因为其考虑的是交集危险事件发生的概率,即放大了实际的危险区域;C测度重现期危险区域的范围在二者之间,因此其对应的峰量组合值也在二者范围之内。峰量组合值的分布区域与危险区域的划分方式呈对应关系,选择不同的重现期计算标准即确定了不同的峰量组合值的分布区域,是确定设计洪水值的前提条件。

4.4 设计洪水值计算

以 10 a、20 a、50 a、100 a、200 a 为重现期,采用极大似然法、条件最可能组合法和双变量同频率法,分别计算同现重现期、联合重现期、C测度重现期标准下的设计洪水值组合,结果见图4。

4.4.1 不同重现期标准设计洪水

对不同重现期标准,比较相同方法推求设计洪水值组合的差异性。图5中A、B、C点分别为在不同重现期标准下利用极大似然法所推求的设计洪水值,以百年一遇洪水为例,同现、联合、C测度重现期的设计洪水峰量值组合分别为(15 606郾06 m3/s,511郾76 亿 m3)、 (18741郾60 m3/s,580郾81 亿 m3)和(16436郾38 m3/s,524郾24 亿 m3),联合重现期的设计洪水值组合中的峰量值最大,C测度重现期次之,同现重现期的最小。由图5可看出,对于任意给定的重现期水平,3种重现期标准的设计洪水值组合间的大小关系均满足上述规律。

重现期的选择对设计洪水值组合起着决定性的作用。在同现重现期标准下得到的设计洪水值组合最小,危险范围考虑得最小,相应的设计标准、工程造价最低,安全性最低。在联合重现期标准下得到的设计洪水值组合最大,危险范围考虑最大,相应的设计标准和工程造价最高,因此安全性也最高。水利工程设计过程中采用同现重现期标准,即以洪峰流量、洪量值同时大于某一峰量设计值作为危险事件,忽略了短时强降雨以及连续降雨引发的洪水的危险因素,低估了洪水的危险范围,确定的防洪设计标准偏低。若以洪峰流量或洪量值中的其中一个变量大于某一峰量设计值作为危险事件,在此基础上确定设计标准,以其对应的重现期为联合重现期,在危险事件的确定上,则考虑了洪峰流量较大而洪量较小、洪峰流量较小而洪量较大这类不构成相应洪水威胁的条件,过分强调安全性,提高了设计标准。上述两种水利工程设计标准都过于极端,仅适用于有特殊设计要求的工程,未能体现出双变量联合分布在水利工程设计中的优势。从定义角度分析,C测度重现期是同现重现期和联合重现期危险事件的折衷考虑,既考虑了洪水特征变量之间的联系,又达到了经济性与安全性的平衡。在双变量频率分析中,推荐采用C测度重现期作为双变量联合分布的重现期计算标准。

4.4.2 不同计算方法设计洪水

以双变量同频率法为基准,比较采用极大似然法、条件最可能组合法推求设计洪水值组合差异,结果见表3~5以及图4。从整体上看,各组合的差异均维持在2%左右,最低为0郾37%,最高为2郾51%,3种设计洪水计算方法具有较强的稳定性,差异性较小。在同现重现期标准下采用3种不同的方法,由极大似然法推求的设计洪水值组合位于双变量同频率法和条件最可能组合法之间,且随着重现期水平的提高,极大似然法推求的设计洪水值组合越接近双变量同频率法推求的设计洪水值组合。在联合重现期以及C测度重现期标准下采用3种不同的方法,由双变量同频率法推求的设计洪水值组合位于极大似然法和条件最可能组合法之间,且极大似然法与双变量同频率法的接近程度更高。可通过选择中间值的方法来确定设计洪水值的计算方法,即在同现重现期标准下,可采用极大似然法,在联合重现期及C测度重现期标准下,可采用双变量同频率法。

表3 同现重现期标准下设计洪水值组合差异

表4 联合重现期标准下设计洪水值组合差异

表5 C测度重现期标准下设计洪水值组合差异

图6为采用双变量同频率法与条件最可能组合法的峰量组合曲线,条件最可能组合法的设计洪水值组合均大于双变量同频率法的对应值,且差距越来越大,但在重现期等值线上的差异维持在2%左右。

图6 峰量组合曲线

5 结 论

a.采用C测度重现期作为双变量频率分析的标准,可确定合理的危险事件的范围,有效规避以同现重现期为标准带来的安全隐患,解决以联合重现期为标准造成的性价比低的问题。

b.在不同的重现期标准下确定不同重现期水平下的设计洪水值组合时,采用的3种方法计算结果较为接近,可通过选择中间值的方法来确定设计洪水值的计算方法。在同现重现期标准下,可选择极大似然法作为双变量设计洪水的计算方法,在联合重现期及C测度重现期标准下,可选择双变量同频率法作为双变量设计洪水的计算方法。

[1]SKLAR M.Functions de repartition a n dimensions et leurs marges[J].Publ Inst Statist Univ Paris,1959(8):229鄄231.

[2]王鹏新,冯明悦,孙辉涛,等.基于主成分分析和Copula函数的干旱影响评估研究[J].农业机械学报,2016,47(9):334鄄340.(WANG Pengxin,FENG Mingyue,SUN Huitao,etal.Droughtimpactassessmentbased on principal component analysis and copula function[J].Transactions ofthe Chinese Society forAgricultural Machinery,2016,47(9):334鄄340.(in Chinese))

[3]关帅,林颖妍,查悉妮,等.基于Copula函数的韩江流域干支流洪水遭遇分析[J].中山大学学报(自然科学版),2015,54(5):130鄄137.(GUAN Shuai,LIN Yingyan,ZHA Xini,et al.Copula function鄄based flood coincidence probability analysis for mainstream and tributary of the Hanjiang River Basin[J].Acta Scientiarum Naturalium Universitatis Sunyatseni,2015,54(5):130鄄137.(in Chinese))

[4]高玉琴,陈钇西,赵立梅,等.秦淮河流域不同频率降雨联合概率分析[J].水电能源科学,2016(3):1鄄5.(GAO Yuqin,CHEN Yixi,ZHAO Limei,et al.Joint probability analysis of different frequencies rainfall of Qinhuaihe River Basin[J].Water Resources and Power,2016(3):1鄄5.(in Chinese))

[5]SALVADORI G,MICHELEC D,DURANTE F.On the return period and design in a multivariate framework[J].Hydrology and Earth System Sciences,2011,15(11):3293鄄3305.

[6]SALVADORI G,DE MICHELE C.Multivariate multiparameter extreme value models and return periods:a copula approach[J].Water Resources Research,2010,46(10):W10501.

[7]SALVADORI G.Bivariate return periods via 2鄄Copulas[J].Statistical Methodology,2004,1(1):129鄄144.

[8]史黎翔,宋松柏.基于Copula函数的两变量洪水重现期与设计值计算研究[J].水力发电学报,2015,34(10):27鄄34.(SHI Lixiang,SONG Songbai.Calculation of return periods and design values of bivariate floods based on Copula function[J].Journal of Hydroelectric Engineering,2015,34(10):27鄄34.(in Chinese))

[9]刘和昌,梁忠民,姚轶,等.基于Copula函数的水文变量条件组合分析[J].水力发电,2014,40(5):13鄄16.(LIU Hechang,LIANG Zhongmin,YAO Yi,et al.Analysis on conditional compositions of hydrological variants based on Copula function[J].Water Power,2014,40(5):13鄄16.(in Chinese))

[10]李天元,郭生练,刘章君,等.基于峰量联合分布推求设计洪水[J].水利学报,2014,45(3):269鄄276.(LI Tianyuan,GUO Shenglian,LIU Zhangjun,et al.Design flood estimation based on bivariate joint distribution of flood peak and volume[J].JournalofHydraulic Engineering,2014,45(3):269鄄276.(in Chinese))

[11]郭生练,刘章君,熊立华.设计洪水计算方法研究进展与评价[J].水利学报,2016,47(3):302鄄314.(GUO Shenglian,LIU Zhangjun,XIONG Lihua.The research progress and the estimate about design flood calculation method[J].Journal of Hydraulic Engineering,2016,47(3):302鄄314.(in Chinese))

[12]芮孝芳,张超.论设计洪水计算[J].水利水电科技进展,2014,34(1):20鄄26.(RUI Xiaofang,ZHANG Chao.Research on design flood calculation[J].Advances in Science and Technology of Water Resources,2014,34(1):20鄄26.(in Chinese))

[13]陈甜,董增川,贾本有,等.考虑上游梯级水库群调度的小南海水电站坝址设计洪水分析[J].河海大学学报(自然科学版),2014,42(6):476鄄480.(CHEN Tian,DONG Zengchuan,JIA Benyou,et al.Analysis of design flood of dam site of Xiaonanhai Hydropower Station considering flood control operation of upstream cascade reservoir group[J].Journal of Hohai University(Natural Sciences),2014,42(6):476鄄480.(in Chinese))

[14]黄振平,林小丽,侯云青,等.P鄄芋型分布参数的矩估计与设计洪水的计算频率[J].河海大学学报(自然科学版),2005,33(1):49鄄51.(HUANG Zhenping,LIN Xiaoli,HOU Yunqing,etal.Momentestimation for parameters of Pearson Type鄄芋 distribution and calculated frequency of design flood[J].Journal of Hohai University(Natural Sciences),2005,33(1):49鄄51.(in Chinese))

[15]李杰友,陆波,穆方驰.棉花滩水库建成后对青溪水库设计洪水影响研究[J].水利水电科技进展,2005,25(1):25鄄27.(LI Youjie,LU Bo,MU Fangchi.Effect of Mianhuatan Reservoir on the design flood of Qingxi Reservoir[J].Advances in Science and Technology of Water Resources,2005,25(1):25鄄27.(in Chinese))

[16]徐长江,郭生练,熊立华,等.英国设计洪水研究进展[J].水利水电科技进展,2004,24(2):66鄄69.(XU Changjiang,GUO Shenglian,XIONG Lihua,et al.Research advances in England design flood[J].Advances in Science and Technology of Water Resources,2004,24(2):66鄄69.(in Chinese))

[17]黄振平,王春霞,马军建.历史洪水重现期的误差对设计洪水的影响[J].河海大学学报(自然科学版),2002,30(1):79鄄82.(HUANG Zhenping,WANG Chunxia,MA Junjian.Effectofestimation errorof historical flood recurrence interval on design flood[J].Journal of Hohai University(Natural Sciences),2002,30(1):79鄄82.(in Chinese))

[18]华家鹏.小流域设计洪水的计算方法[J].河海大学学报 (自 然 科 学 版),1998,26(2):100鄄104.(HUA Jiapeng.Calculation method of design flood for small watershed[J].Journal of Hohai University(Natural Sciences),1998,26(2):100鄄104.(in Chinese))

[19]NELSENR B.An introduction to Copulas[M].New York:Springer,1999.

[20]吴书悦,赵建世,雷晓辉,等.气候变化对新安江水库调度影响与适应性对策[J].水力发电学报,2017,36(1):50鄄58.(WU Shuyue,ZHAO Jianshi,LEI Xiaohui,et al.Impacts of climate change on operation of Xin爷an River Rreservoirand adaption strategies[J].Journalof Hydroelectric Engineering,2017,36(1):50鄄58.(in Chinese))

[21]郑艳妮,闻昕,方国华.新安江流域气候变化及径流响应研究[J].水资源与水工程学报,2015,26(1):106鄄110.(ZHENG Yanni,WEN Xin,FANG Guohua.Research on climate change and runoff response in Xin爷an River Basin[J].JournalofWaterResources & Water Engineering,2015,26(1):106鄄110.(in Chinese))