基于Copula函数的洪水峰量联合设计方法研究

苗正伟，丁志宏，路 梅，杨学智

(1.河北水利电力学院 水利工程学院，河北 沧州 061001； 2.中水北方勘测设计研究有限责任公司，天津 300222)

1 研究背景

洪水频率是水文统计的重点内容，传统的单变量频率分析或者研究洪峰频率，或者研究洪量频率，或者假定二者同频率，这既违背了洪水是多变量随机事件的科学事实，也忽略了洪水特征变量之间普遍存在的较强相关性。有不少学者基于多变量频率分析方法来描述洪水、干旱等水文事件，但传统多变量联合分布对边缘分布的苛刻要求与实际水文现象的高度复杂性格格不入，极大地限制了传统多变量联合分布的应用[1]。近年来，基于Copula函数的联合分布在水文领域越来越受关注，Sklar[2]在1959年将Copula概念引入统计领域，继而由Genest等[3]、Harry等[4]进一步发展、完善了Copula相关理论，随后，Nelsen[5]在其专著中对Copula函数理论与成果进行了总结，Marius Hofert等[6]进一步总结了Copula函数的性质、应用及其R语言的实现。众多研究结果表明，Copula函数很适于描述水文特征变量间的相关性，包括非线性相关关系，且对边缘分布没有限制，已被广泛应用于暴雨、洪水、干旱等水文事件的多变量联合分布研究[7-11]。

重现期是水文设计的核心概念之一。在统计意义上，重现期是发生指定事件的时间间隔的平均值。对单变量事件而言，重现期和设计值是一一对应的。但对于多变量事件，因变量间存在无穷多种组合，其重现期的定义以及设计值的计算都要复杂得多，这已成为近年来水文领域的热点问题[12-14]。定义双变量洪水重现期的核心问题是危险域的界定，所谓危险域[15]就是所有危险事件的集合，危险域之外的二维实数域即为安全域，但一般将安全域与危险域的分界线称为临界域。给定洪水峰量(峰量即洪峰流量和洪量)组合，若其落在危险域内，则将其判定为超标洪水，属于危险事件，反之则是安全的或处于临界状态。根据双变量框架下危险域界定方式的不同，目前对双变量洪水事件重现期的定义也有很多种，其中，“或”(“OR”)和 “且”(“AND”)重现期应用最广泛， OR重现期将洪水峰量中至少有一个超过指定阈值的事件集合定义为危险域；而AND重现期将危险域定义为两变量同时超过指定阈值的事件集合。根据不同的危险域进而计算相应的重现期。但是，这2种重现期对于危险事件的判断均不尽合理[16-18]。为此，Salvadori 等[17]根据Kendall函数得到Copula等值线，并将峰量联合分布概率超过给定阈值的事件集合定义为危险域，进而提出了Kendall 重现期的概念，这使得危险事件的判断不再有逻辑上的矛盾[19-21]。但该方法给出的安全域却是无限的，这与事实明显不符。因此，Salvadori 等[17]于2013年又提出了生存 Kendall重现期的概念，将危险域定义为洪水峰量的联合生存概率小于给定阈值的事件集合。生存Kendall重现期不仅对危险域的判定科学合理，而且它给出的安全域也是有界的，与现实相符，是当前双变量水文事件相对严谨的重现期定义方式。

与单变量情况不同，给定的双变量重现期会对应无穷多种变量组合，由此构成重现期等值线，这意味着，双变量框架下联合洪水设计值并非唯一，而不同的洪水峰量组合将得到不同的洪水过程线，从而对调洪演算的结果产生直接影响，并进一步影响水库特征水位及水库规模，因此，多变量框架下如何科学合理地计算联合洪水设计值是至关重要的[22]。目前，计算联合洪水设计值的常用方法主要包括同频率法[23]、极大似然法[18,24]、条件最可能组合法[25-27]等。

当前关于洪水的双变量设计主要基于OR、AND以及Kendall重现期，较少涉及生存Kendall重现期，而且缺乏洪水设计值计算方法的对比分析。为此，本文以岗南水库为例，基于Copula函数，采用OR、AND、Kendall及生存Kendall 4种重现期标准以及同频率、极大似然、条件最可能组合3种洪水设计值计算方法，对洪水的峰量联合分布展开研究，以期为水文联合设计提供科学参考，为岗南水库的防洪安全提供理论依据。

2 数据与方法

2.1 研究区与数据

位于海河流域子牙河系滹沱河中游的岗南水库控制流域面积15 900 km2(如图1所示)，总库容15.71亿m3，其中9.71亿m3为防洪库容，主要担负着下游京九铁路、华北油田及石家庄等的防洪保护任务，同时承担向灌区及石家庄等地区供水的任务，该水库对支撑河北省水资源可持续利用及社会经济发展至关重要。本研究所用数据主要是1949—2012年岗南水库流域年最大洪峰流量和六日洪量，此外还包括一场重现期为3 000 a的古洪水及11场历史特大洪水，数据来自水文年鉴及海河水利委员会。

图1 岗南水库流域Fig.1 Drainage map of Gangnan reservoir

2.2 研究方法

3 结果与分析

3.1 洪峰和六日洪量的边缘分布

对岗南水库流域洪峰流量及六日洪量系列，由非连续样本的线性矩法初估参数进而进行适线，经验频率点据、拟合的统计参数及P-III分布曲线如图2所示。由此可见，适线效果良好。

图2 P-III分布拟合结果Fig.2 Fitting results of P-III distribution

3.2 洪峰和六日洪量的联合分布

表1 Copula函数参数估计和统计检验结果Table 1 Estimated parameters of Copula functions and statistical test results

图3 联合观测值的经验分布与理论分布比较Fig.3 Comparison between empirical and theoretical distributions of joint observations

3.3 洪水峰量联合分布重现期比较分析

给定重现期水平20、50、100、200 a，在同一坐标系中绘制4种重现期标准的等值线，如图4所示。由此可知，对于不同标准的重现期，OR与Kendall重现期等值线走势一致；AND与生存Kendall重现期走势一致。尽管重现期等值线趋势相似，但它们确定的危险域截然不同，OR、AND重现期等值线上每一点都对应不同的危险域，正因如此，这二者对危险域的判断存在逻辑上的矛盾，而Kendall、生存Kendall重现期等值线上各点都有一样的危险域，这使得两者对危险域的判定更为合理。图4还显示： Kendall重现期给出的安全域是无界的，这显然与实际相悖，与之相反，生存Kendall重现期给出了有界的安全域，更契合工程安全的内在要求，综上所述，4种重现期标准中，生存Kendall重现期最为科学合理。由图4还可知，4种给定重现期水平下，不同重现期标准等值线间的关系是一致的，均表现为：OR重现期等值线在最上、AND重现期等值线在最下，Kendall与生存Kendall重现期等值线居中；AND、OR、Kendall 3种重现期等值线无交点；生存Kendall重现期等值线与Kendall重现期等值线有2个对称的交点，与OR、AND重现期等值线无交点。

图4 重现期等值线Fig.4 Isolines of return periods

3.4 设计值比较分析

3.4.1 不同计算方法的洪水设计值比较

采用同频率、极大似然、条件最可能组合3种方法计算各重现期标准的设计值，结果如表2、表3所示。由此可知,对所有给定的重现期水平和重现期标准，无论是洪峰流量还是六日洪量，都有如下规律成立：同频率法与极大似然法所计算的设计值差别非常小，条件最可能组合法计算的洪峰流量略大于极大似然法和同频率法的洪峰流量，而六日洪量略小于极大似然法和同频率法的六日洪量。不失一般性，以100 a生存Kendall重现期为例，同频率法与极大似然法求得的洪峰流量设计值仅相差0.334 m3/s，条件最可能组合法比同频率法的结果大差值为65.62 m3/s，相对差异仅为0.69%。六日洪量设计值的差异与洪峰设计值的差异类似。由此可见，对于洪水设计而言，极大似然法、条件最可能组合法、同频率法3种不同计算方法对洪水设计值并无明显影响，其中尤以同频率法计算最为简单，因此建议使用同频率法进行双变量洪水设计值计算。

3.4.2 不同重现期标准的洪水设计值比较

由表2、表3可知，给定重现期水平，无论何种洪水设计值计算方法，计算所得的洪峰、洪量设计值都遵循以下规律：OR重现期标准的洪水设计值>生存Kendall重现期洪水设计值>单变量洪水设计值>Kendall重现期洪水设计值>AND重现期洪水设计值；各重现期标准的洪水设计值与单变量洪水设计值的相对差异基本呈现随重现期水平增大而减小的趋势，而且AND重现期、OR重现期洪水设计值的相对差异明显大于Kendall和生存Kendall重现期洪水设计值的相对差异。因此，以OR重现期标准设计工程将偏于保守，工程规模偏大，投资高；而以AND重现期标准进行设计，则工程规模偏小，投资少，但安全性降低。由此可见，AND和OR是工程设计的两种极端情形，它们在工程规模与投资、安全性方面是失衡的，相较而言，根据生存Kendall或Kendall重现期进行工程设计更为科学合理。又因为生存Kendall重现期的定义更加严谨，因此，生存Kendall重现期是兼顾设计洪水计算的合理性以及工程的经济性、安全性的一个较好选择，更利于洪水的风险管理。

4 讨 论

由表2、表3可知：本文计算的不同重现期标准的洪水设计值与单变量洪水设计值相对差异的绝对值绝大部分在9%以内，差别较小，尤其是当重现期水平超过1 000 a时，各重现期标准、各计算方法所得的联合洪水设计值与单变量洪水设计值的相对差异都在3%以内。当将Copula参数减小为3.000 0，联合洪水设计值与单变量洪水设计值相对差异见表4和表5(联合洪水设计值计算方法为同频率法)。由此可见，当Copula参数减小时，相对单变量洪水设计值而言，Kendall、AND、OR 3种重现期标准联合洪水设计值的差异都明显增大，增幅在100%左右，而生存Kendall重现期的联合洪水设计值差异也有小幅增加。根据二维对称Gumbel Copula参数与Kendall相关系数的关系进一步分析可知，当Copula参数减小时，Kendall相关系数也是减小的，即洪水的峰量相关性减弱，而本研究采用的Copula参数为5.286 2，相应的Kendall相关系数高达0.810 8，当Copula参数为3.000 0时，Kendall相关系数减小为0.666 7，这说明：洪水峰、量相关性越弱，基于Copula的双变量联合洪水设计值与单变量洪水设计值的差异越明显。

表2 洪峰流量设计值比较Table 2 Comparison of design peak discharge

表3 六日洪量设计值比较Table 3 Comparison of design 6d flood volume

表4 Copula参数对洪峰流量设计值的影响Table 4 Influence of parameter of Copula on design peak discharge

表5 Copula参数对六日洪量设计值的影响Table 5 Influence of parameter of Copula on design 6 d flood volume

图4显示，生存Kendall重现期等值线与Kendall重现期等值线围成的闭合区域在Kendall重现期概念下，被判定为危险域，而在生存Kendall概念下被判定为安全域，相较于洪量很大而洪峰很小、或者洪峰很大而洪量很小的情形，该交集区域内的峰量关系更契合现实洪水峰量的相关性，即这类洪水事件更可能在现实中发生。因此，该区域究竟被判定为危险域还是安全域，对防洪工程的规模具有直接影响。因为生存Kendall重现期将该闭合区域判定为安全域，所以生存Kendall重现期确定的洪水设计值比Kendall重现期的洪水设计值更大(图4中的同频率设计值及表2、表3的结果都证明了该种情况)。因此，基于生存Kendall重现期确定的工程规模也比Kendall重现期的更大，即以生存Kendall重现期进行工程设计比Kendall重现期进行工程设计更安全。

5 结 论

本文以岗南水库为例，采用4种重现期标准、3种联合洪水设计值计算方法对双变量洪水设计值进行了研究，主要结论如下。

(1)AND和OR重现期在危险域和安全域的识别上存在局限性，相对而言，Kendall重现期更合理，但其安全域是无界的，这与实际不符。生存Kendall重现期则界定了有界的安全域，使得重现期的概念在逻辑上更科学合理。

(2)同频率法、极大似然法、条件最可能组合法3种洪水设计值计算方法所得结果相差不大，尤其是同频率法与极大似然法，二者结果几乎相等，从简单实用的角度出发，推荐采用同频率法计算联合洪水设计值。

(3)基于Copula计算给定重现期水平的洪水设计值，既受重现期标准的影响，也受设计值计算方法的影响。但不同设计值计算方法对结果的影响较小，因此，设计结果主要取决于重现期标准，且遵循以下规律：OR重现期标准的洪水设计值>生存Kendall重现期洪水设计值>单变量洪水设计值>Kendall重现期洪水设计值>AND重现期洪水设计值。从计算过程科学合理、结果安全可靠且兼顾经济性的角度来说，推荐采用生存Kendall重现期进行双变量洪水设计。

(4)联合洪水设计值与单变量洪水设计值的差异受变量间相关性的影响较大，且变量相关性越弱，差异越大。