基于Copula函数的桂平航运枢纽水闸设计洪水位
2016-03-01甘富万黄永俊孙晋东
甘富万，黄永俊，熊 健，2，孙晋东，胡 泊
(1.广西大学土木建筑工程学院， 广西南宁530004；2.广西水利电力勘测设计研究院， 广西南宁530023)
基于Copula函数的桂平航运枢纽水闸设计洪水位
甘富万1，黄永俊1，熊健1，2，孙晋东1，胡泊1
(1.广西大学土木建筑工程学院， 广西南宁530004；2.广西水利电力勘测设计研究院， 广西南宁530023)
摘要：针对桂平航运枢纽水闸原设计洪水位由经验方法确定、缺乏理论依据的情况，利用Copula函数理论建立二变量设计洪水组合及洪水位的确定方法，并应用于桂平航运枢纽水闸设计洪水位的确定。结果表明：建立的4种常用Copula函数中，Clayton Copula函数拟合效果最好；设计洪水组合应该包括Copula同现重现期等于设计重现期的所有洪水组合；按照建立的方法计算得到桂平航运枢纽水闸设计洪水位为43.26 m，通过经验方法确定的原设计洪水位为43.48 m，可达到防洪要求，且略偏安全。
关键词：Copula函数；桂平航运枢纽水闸；重现期；设计洪水组合；设计洪水位
桂平航运枢纽位于广西桂平市郁江河段，是一个集航运、发电、灌溉、交通于一体的综合利用性航运枢纽工程，处于郁江与西江干流(包括黔江、浔江)汇合口附近的郁江段，主要由航运船闸、溢洪坝及水闸组成[1]，见图1。非洪水时期，溢流坝水闸关闭，枢纽坝上水位抬高以保证航运水位以及发电；洪水时期，溢洪坝水闸打开，洪水通过溢洪坝水闸下泄。桂平航运枢纽同时受到郁江和西江干流洪水的作用，汇合口的郁江上游附近有贵港水文站、干流下游附近有大湟江口站水文站[2]。桂平航运枢纽水闸设计重现期为100 a一遇，设计时采用大湟江口站重现期为100 a的洪峰流量与0.2倍该洪峰值作为溢洪坝坝址设计洪水组合，并经过调洪计算得到溢洪坝水闸设计洪水位为43.48 m。该设计洪水位的确定方法为经验方法，缺乏理论依据，所得设计洪水位是否达到设计重现期的要求，有待商榷。
图1 桂平航运枢纽位置及组成
桂平航运枢纽水闸设计洪水问题是典型二变量洪水组合问题。Copula函数可描述变量间的相关性，是一类将联合分布函数与边缘分布函数连接在一起的函数，在多变量水文研究中已经被证明有着广泛的推广应用价值，相关理论参见文献[3]～文献[11]。本文尝试利用Copula函数对桂平航运枢纽的设计洪水组合进行计算，并以其为基础进一步确定桂平航运枢纽溢洪坝水闸的设计洪水位。
1Copula函数的建立
1.1函数的选择
首先选择已被证明在多变量水文问题中可以有效应用的四种Copula函数，然后根据各种Copula函数的要求，利用大湟江口站和贵港站历年实测洪峰系列资料建立Copula函数并进行对比评价，以选择最为合适的Copula函数形式。初步选择的4种Copula函数形式见表1。
表1 常用copula函数形式
表1中各式的C为Copula连接函数，u1和u2为变量X1和X2的边缘分布函数，u1=F1(x1)，u2=F2(x2)，θ为Copula函数的参数。
1.2函数的建立
①边缘分布
边缘分布统一采用为P-III型分布。为充分考虑郁江、西江干流洪水的遭遇特点，本次分两种情况分别建立Copula函数：以西江干流洪峰系列与相应时间的郁江洪峰系列进行组合(方案1)；以郁江洪峰系列与相应时间的西江干流洪峰系列进行组合(方案2)。
方案1具体做法：西江干流选择大湟江口站1951年～2010年实测洪水系列，郁江以大湟江口站发生最大洪峰流量的时间为控制，统计贵港站相应场次洪水的相应洪峰流量并进行频率计算和P-Ⅲ型曲线适线。
方案2具体做法：郁江选择贵港站1936年～2010年实测洪水系列资料，西江干流以贵港站发生最大洪峰流量的时间为控制，统计大湟江口站同场洪水的相应洪峰流量，并以此统计P-Ⅲ型曲线参数。
P-III型分布曲线也称伽马分布，它的概率密度函数[12]是：
(1)
式中，Γ(α)是α的伽玛函数；α、β、a0分别是P-Ⅲ型分布函数的形状、尺度和位置参数，α>0，β>0。
(2)
经过频率计算和P-Ⅲ型曲线适线分析，得到边缘分布参数统计结果见表2和表3。
表2 边缘分布函数的统计参数值
表3 边缘分布函数的α、β、a0参数值
②Copula函数参数估计
采用相关性指标法[12]计算Copula函数的计算参数θ，即首先计算各变量序列的Kendall秩τ，然后再利用τ计算θ。各Copula函数参数θ与Kendall秩相关系数τ之间的关系见表4。
Kendall秩相关系数τ的计算方法可参考文献[13]，分别利用实测数据可求得方案1和方案2的τ为0.147 6和0.219 6，并进一步计算得到各Copula函数的参数θ值见表4。
表4 Copula函数参数θ与 Kendall 秩相关系数τ转换关系
1.3评价分析
本文利用Copula函数常见的AIC信息准则法以及离差平方和最小准则(OLS)法进行评价[13]。
AIC信息准则形式如下：
(3)
AIC=nln(MSE)+2k，
(4)
式中：Femp(xi1,xi2,…,xim)表示经验频率值，C(ui1,ui2,…,uim)表示理论频率值，m为函数的维数，k为模型参数的个数。
离差平方和最小准则形式如下：
(5)
式中各符号意义同上。
AIC值、OLS值越小，Copula函数拟合越好。运用AIC信息准则以及离差平方和最小准则进行计算，结果见表5。由计算成果可见，各copula函数的AIC值、OLS值非常接近，可见拟合效果都比较好。其中AIC值、OLS值最小的函数为Clayton Copula函数，下文主要利用Clayton Copula函数成果进行分析。
表5 AIC值和OLS值计算成果
2基于Copula函数的设计洪水组合
2.1设计洪水重现期
桂平航运枢纽水闸设计洪水重现期为100 a一遇。Copula函数的重现期分为联合重现期、同现重现期和条件重现期三种，不同类型重现期对应的洪水组合不同，必须明确哪种重现期可以作为水闸设计洪水的重现期。
假定洪水事件中的洪峰系列分别为X1和X2。u1=F1(x1)=P(X1>x1)、u2=F2(x2)=P(X2>x2)为对应洪水边缘分布函数，则对于系列X1的单变量洪水重现期可表示为：
T(x1)=1/P(X1>x1)。
(6)
当以此重现期作为设计重现期T设，设计洪峰为(x1)，可以保证洪水事件(X1
比较单变量设计洪水重现期的定义，在二变量洪水组合问题中，如果设计的洪水组合为(x1,x2)，则发生洪水事件(X1
T设=T(x1,x2)=1/P(X1>x1,X2>x2)。
(7)
Copula函数同现重现期是指两个变量同时超过特定值时事件发生的重现期，计算公式如下：
T0(x1,x2)=1/P(X1>x1,X2>x2)=1/(1-u1-u2+C(u1,u2))。
(8)
对比可知，应该以Copula函数的同现重现期作为本次二变量设计洪水重现期T设。
2.2设计洪水组合
对于单变量洪水事件，洪峰—重现期关系为一条曲线，如图2所示。设计洪水重现期T设(如100 a一遇)对应的洪峰为曲线上该重现期对应的一个值，以该洪峰值作为设计洪峰，可保证重现期小于100 a一遇的洪峰出现时的行洪安全。
在本次研究的情况中，水闸同时受到郁江、西江干流洪水作用，是二变量洪水组合问题。T设(如100 a一遇)对应的洪水不是一个洪峰，而是郁江、西江干流的洪峰组合，即设计洪水组合。郁江、西江干流洪峰与重现期的关系为一个曲面，如图3所示。T设(如100 a一遇)对应的洪峰组合则构成了一条曲线，而重现期小于100 a一遇对应的洪水组合构成一个曲面。要使得设计洪水重现期达到100 a一遇，则必须保证所有重现期小于100 a一遇的洪水组合的行洪安全，此时设计洪水组合应该包含设计重现期(如100 a一遇)对应的曲线上的所有洪水组合。
图2单变量洪水事件设计洪峰
Fig.2Design flood of single variable flood
图3二变量洪水事件设计洪水组合
Fig.3Design flood combination of two variable flood
3基于Copula函数的水闸设计洪水位
3.1计算方法
利用建立的Clayton Copula函数计算得到桂平航运枢纽同现重现期为100 a一遇的洪水组合，并进行调洪计算可得到水闸的设计洪水位。由于同现重现期100 a一遇对应的洪水组合有无数个，本次采用抽取大量代表洪水组合进行调洪计算，并取调洪水位最大值作为设计值的方法。
设计洪水组合下的调洪计算方法直接参考枢纽原设计的计算方法[1]。
水闸泄洪能力的计算公式如下：
(9)
图4 坝址水位流量关系Fig.4 Water level and flux relationship of dam site
坝址自然水位同时受到郁江和西江干流洪水的作用，根据不同的设计洪水组合，利用“贵港站流量—大湟江口站流量—坝址自然水位”关系插值计算得到不同的设计洪水组合下的闸下水位，再根据公式(9)计算得到闸上水位，即为不同设计洪水组合对应的的水闸设计洪水位。
利用实测资料整理溢洪坝坝址的自然水位流量关系[1]，见图4。
3.2计算成果
针对方案1和方案2两种方案建立的Clayton Copula函数，计算大量重现期为100 a一遇的洪水组合并进行调洪计算，得到各设计洪水组合的闸上洪水位，典型设计洪水组合计算结果见表6。
表6 典型设计洪水组合调洪计算成果(T设=100 a)
由计算成果可见，桂平航运枢纽水闸闸上水位主要受到西江干流洪水的控制，在各设计洪水组合中，西江干流洪水越大，闸上洪水位越高。洪水位最大的洪水组合为西江干流洪水接近100 a一遇(对应西江干流洪峰48 200 m3/s，郁江洪峰966 m3/s)的洪水组合，对应洪水位为43.26 m，此值即为桂平航运枢纽水闸的设计洪水位。
桂平航运枢纽水闸原设计采用的设计洪水组合为西江干流48 200 m3/s、郁江9 640 m3/s，对应原设计洪水位43.48 m，高于本次利用Copula函数计算得到的设计洪水位0.22 m。可见，原设计洪水位虽然按照经验方法确定，但达到了100 a一遇设计重现期的要求，而且略偏安全。
4结语
本文建立了郁江和西江干流洪水的Copula函数，并探讨了基于Copula函数的二变量设计洪水组合及设计洪水位的计算方法，最后对桂平航运枢纽水闸的设计洪水位进行了计算。结论如下：
①建立的4种常用Copula函数中，Clayton Copula函数拟合效果最好。
②设计洪水组合应该包括Copula同现重现期等于设计重现期的所有洪水组合。
③按照本文建立的方法计算得到桂平航运枢纽水闸设计洪水位为43.26 m。通过经验方法确定的原设计洪水位43.48 m可达到防洪要求，且略偏安全
参考文献：
[1]徐永乐.西江航运建设第一期工程广西段初步设计报告(第一册 桂平航运梯级)[R]. 南宁：广西壮族自治区交通勘测设计院，1983.
[2]甘富万，胡秀英，刘欣，等.广西境内西江流域洪水特性分析[J]. 广西大学学报(自然科学版)，2015，40(1)：244-250.
[3]GRIMALDI S, SERINALDI F.Design hyetographs analysis with 3 -Copula function [J]. Hydrological Sciences, 2006, 51(2): 223-238.
[4]SHIAU J T, SONG F, NADARAJAH S.Assessment of hydrological droughts for the Yellow River, China, using Copulas [J]. Hydrological Processes, 2007, 21: 2157-2163.
[5]DEMICHELE C, SALVADORI G, CANOSSI M, et al.Bivariate statistical approach to check adequacy of dam spillway[J]. Journal of Hydrologic Engineering, 2005, 10(1): 50-57.
[6]郭生练，闫宝伟，肖义，等.Copula函数在多变量水文分析计算中的应用及研究进展[J]. 水文，2008，28(3)：1-7.
[7]张娜，郭生练，闫宝伟，等.Copula函数在分期设计洪水中的应用研究[J]. 水文，2008，28(5)：28-32.
[8]张娜，郭生练，刘攀，等.基于Copula函数法推求分期设计洪水和汛限水位[J]. 武汉大学学报(工学版)，2008，41(6)：33-40.
[9]闫宝伟，郭生练，郭靖，等.基于Copula 函数的设计洪水地区组成研究[J]. 水力发电学报，2010，29(6)：60-65.
[10]陆桂华，张亚洲，肖恒，等.气候变化背景下蚌埠市暴雨与淮河上游洪水遭遇概率分析[J]. 气候变化研究进展，2015，11 (1)：31-37.
[11]刘章君，郭生练，李天元，等.梯级水库设计洪水最可能地区组成法计算通式[J]. 水科学进展，2014，25(4)：575-584.
[12]侯芸芸，宋松柏，赵丽娜，等.基于Copula 函数的3 变量洪水频率研究[J]. 西北农林科技大学学报(自然科学版)，2010，38(2)：219-228.
[13]于延胜，陈兴伟.基于Mann_Kendall法的径流丰枯变化过程划分[J]. 水资源与水工程学，2013，24(1)：60-63.
(责任编辑唐汉民梁碧芬)
Design flod level of sluice of Guiping Shipping Hub based on Copula function
GAN Fu-wan1, HUANG Yong-jun1, XIONG Jian2, SUN Jin-dong， HU Bo1
(1. College of Civil Engineering and Architecture, Guangxi University, Nanning 530004, China;
2. Guangxi Water and Power Design Institute, Nanning 530023, China)
Abstract:In view of that the design flood level of the Guiping Shipping Hub was calculated by an empirical method which lacked of theoretical basis, a new method based on Copula function to calculate the design flood combination of two variables and the flood level was established and applied to the sluice of the Guiping Shipping Hub. The result shows that Clayton Copula has the best fitting effect among the four common Copula functions; that the design flood combination should include all combinations, in which Copula function occurrence recurrence period is equal to the design recurrence period of the Guiping Shipping Hub sluice; that the design flood level of the sluice of Guiping Shipping Hub is 43.26 m according to the new method; that the original design flood level (43.48 m) is slightly higher and meet the design standard.
Key words:Copula function; Sluice of Guiping Shipping Hub; recurrence period；design flood combination; design flood level
中图分类号：TV122+.3
文献标识码：A
文章编号：1001-7445(2016)01-0284-07
doi:10.13624/j.cnki.issn.1001-7445.2016.0284
通讯作者：甘富万(1981—)，男，广西北流人，广西大学副教授，博士；E-mail:gxugfw@126.com。
基金项目：国家自然科学基金资助项目(51178124)；广西自然科学基金资助项目(2014GXNSFBA118263)；广西大学科研基金资助项目(土科基-2012-02)
收稿日期：2015-11-20；
修订日期：2015-12-21
引文格式：甘富万，黄永俊，熊健,等.基于Copula函数的桂平航运枢纽水闸设计洪水位[J].广西大学学报(自然科学版)，2016，41(1)：284-290.
