黄河中游河南段设计洪水及其地区组成研究

程 旭，马细霞，肖 遥，2，王倩丽

（1.郑州大学水利科学与工程学院，郑州450001；2.中水东北勘测设计研究有限责任公司，长春130021）

0 引言

设计洪水是水利工程规划设计、防洪策略和水资源管理的重要水文数据，由流量资料推求设计洪水是设计洪水计算的主要方法［1，2］。传统上，通过拟合各历史洪水要素观测值概率分布函数以估计其不同重现期下的单变量设计洪水值，以其峰值、洪量和持续时间这3个内部存在物理联系而并不相互独立的随机变量分析洪水过程，对于一个完整的洪水事件，仅依靠单变量概率分布分析存在一定的局限性，不能完整描述洪水过程［3，4］。因此要进行更完整的设计洪水频率分析，就需要对各洪水要素进行多变量联合统计分析［5，6］。设计洪水地区组成是流域开发方案设计以及上下游、干支流工程联合调洪运行的重要依据［7，8］，目前，常用的地区组成方法有同频率法和典型年法，该类方法主观地考虑“对防洪不利”［9-11］。虽然便于使用和规划，但更多是主观臆断而缺少客观的判断，不能充分反映各分区之间的空间相关性，也无法定量地描述各分区与设计断面之间的相关程度［12］。Copula 函数是将变量联合累积分布函数同变量边缘累积分布函数连接起来的函数，能够构建各洪水要素之间的联合分布，也能够构建各分区洪水之间的联合分布［13-15］。

为求得同时考虑洪水要素相关性和空间相关性的设计洪水结果，本文以黄河中游河南段洪水地区组成复杂的花园口水文站上游为研究对象，选取12天洪量、洪峰流量两个洪水要素，基于Copula 函数建立两个洪水要素以及不同断面之间洪水要素的联合分布，并与遗传算法相结合构建Copula-GA（Copula-Genetic Algorithms）算法，推求最有可能的设计洪水地区组成，以期为防洪预案制定提供参考依据。

1 研究区概况及数据来源

花园口水文站位于黄河中游，控制流域面积73 万km2，占黄河流域总面积的92%［16］。花园口断面上游干流河南段建有三门峡和小浪底两座大型水库，流域水系图如图1所示，研究区水系纵横，呈不对称分布，水患灾害频发，各量级洪水以上大洪水和上下同时大洪水为主占95.2%，即以全流域性洪水为主［17］。各防洪断面洪水组成成分的确定对于水库间防洪优化联调工作的开展起着至关重要的作用［18］。花园口断面的洪水概率分布直接受到干支流水库调节作用的影响［19］。所以，应与黄河中游各水系实际情况相结合，确定花园口洪水地区组成方案。本文将以花园口站过境洪量为控制标准，探求花园口以上，小浪底和相应区间洪水的相关关系及其地区组成。

图1 研究区水系图Fig.1 Water system diagram of the study area

本文基本资料为1950-2015年黄河干流河南段的花园口、小浪底以及小-花区间各站点的实测流量。1958年7月份洪水是花园口水文站记载以来的最大洪水，花园口站洪峰流量高达22 300 m3/s，该场洪水预见期短，对下游威胁较大，具有一定的代表性，对该场实测洪水过程进行插值，对洪峰所在时段取实测洪峰值，得到分辨率为6 h的洪水典型洪水过程。

2 研究方法

本文基于Copula 函数特性以及研究区洪水情况，建立花园口站洪量和洪峰流量的联合分布，计算相应重现期的时段洪量和洪峰流量，并推求设计洪水过程线；采用Copula 函数构造各干支流洪水要素之间的联合分布，并与遗传算法结合构建Copula-GA算法，推算各干支流设计洪水的最可能地区组成。

2.1 两变量重现期

选取时段洪量、洪峰流量两个变量作为洪水特征要素，应用一组随机变量（Qi，Wi）描述相互独立的洪水事件系列｛I1，I2，...，In｝中的每一场洪水事件Ii，假设其联合分布为F（q，w）。在设计洪水要素推求过程中，洪峰和洪量的表示方式为：

将以上单变量事件整合为联合变量事件，表达方式为：

式中：Q和W分别表示洪水事件的洪峰流量和洪量；洪水事件Iq，w表示Q，W至少一个超过设定阈值。

若Q和W为年最大值，其单变量重现期为Tq、Tw［20］：

此时，Eq，w事件的重现期为Tq，w：

FQ（q）和FW（w）均大于或等于F（q，w），由（4）和（5）可得：

假定采用单变量方法得到的T年一遇设计洪峰qT和洪量wT，即变量Q和W分别取T年一遇的设计值，根据式（6）可得：

从上式可知，若要求Tq，w（qT，wT），则必定q和w分别大于或等于qT和wT。

2.2 建立洪峰流量和洪量的联合分布

选取洪峰流量以及洪量两个洪水要素，使用Copula 函数建立两者的联合分布，本文以GH-Copula函数为例［21，22］：

式中：FQ（q）和FW（w）分别为洪峰流量Q和洪量W的边缘分布；θ为函数参数，变量之间的相关性越强，其值越大，θ的参数估计为：

式中：τ为Kendall 秩相关系数；（xi-xj）为观测点据；sign（*）为符号函数，当（xi-xj）（yi-yj）＞0 时，sign（*）=1；当（xi-xj）（yi-yj）＜0 时，sign（*）=-1；当（xi-xj）（yi-yj）=0时，sign（*）=0。

对于特定重现期Tq，w，将式（8）带入式（5）得［23］：

式中：Tu，v=Tq，w。可以满足式（13）的q，w的无数种组合，在平面上绘制q，w关系图，会出现相应重现期下的等值线，借助于式（8）～（10），式（13）可以表示为：

式中：v=L（u）为相应重现期下等值线的表达式。

根据式（14），可以得到重现期为Tq，w洪峰流量q与洪量w的等值线。该等值线上所有点的重现期均为Tq，w，当Tq，w=T时，该等值线上点（q，w）分别大于相应的qT和wT值。在设计洪水的洪峰流量频率u与洪量频率v相等时，有唯一解，即等值线v=L（u）与直线v=u的唯一交点，其解析表达式为：

式中：qTO和wTO为两变量重现期Tq，w（Tq，w=Tu，v）对应的Q和W的设计值。

2.3 推求设计洪水过程线

根据所得到的设计洪峰流量qTO、设计洪量wTO以及典型洪水过程，采用变倍比放大法获得设计洪水过程［24］：

式中：DF（t）、TF（t）分别为设计洪水和典型洪水在t时刻的流量；QTD和WTD分别为典型洪水过程的洪峰流量和历时为TD的洪量。

2.4 设计洪水地区组成及遗传算法的优化

2.4.1 洪水地区组成的数学描述

设计洪水地区组成的核心内容是在汇流洪量一定的情况下，确定洪水来源的各分区洪水所占比例［25，26］。根据流域水系分布情况，本文重点研究花园口水文站、小浪底水库及小-花区间的设计洪水地区组成。如图2所示，A为小浪底水库，X为小浪底入库洪水；B为小浪底-花园口区间，Y小浪底-花园口区间洪水；Z为花园口水文站断面，其中Z=X+Y。

图2 研究区内洪水地区组成示意图Fig.2 Schematic diagram of the flood region composition of study areas

本文各分区独立选样并采用P-III型分布作为频率曲线，其概率密度函数为［27］：

式中：α、β、a0分别为P-III 分布的形状、尺度以及位置参数，α＞0，β＞0；Γ(α)为α的伽马函数。

α、β、a0与总体的统计参数Xˉ、Cv、Cs的关系为：

2.4.2 同频率地区组成法

根据防洪要求，上游分区的洪量与下游设计断面洪量重现期（设计频率）相同。对于本研究，同频率法可以有以下的两种组合形式，分别为：当花园口断面发生重现期为T（频率为P）的洪水zp时，小浪底水库也发生重现期为T的洪水xp，小-花区间发生相应洪水y，有y=zp-xp；当花园口断面发生重现期为T（频率为P）的洪水zp时，小-花区间也发生重现期为T的洪水yp，小浪底水库发生相应洪水x，即x=zp-yp。

2.4.3 最可能洪水地区组成

针对传统洪水地区组成算法缺乏对空间相关性的考虑，本文利用GH-Copula 函数构建研究区内小浪底洪水X与小-花区间洪水Y之间的联合分布，其中，各个断面独立选样并计算边缘分布，以此推求在花园口站发生一定量级的洪水时，小浪底和小-花区间最可能的设计洪水过程。

利用GH-Copula 函数建立小浪底水库洪水X和小-花区间洪水Y的联合分布如下式：

式中：FX（x）和FY（y）分别为小浪底水库洪水与小-花区间洪水的P-III型分布函数。

由式（19）可得联合密度函数如下：

当下游花园口断面发生设计洪水zp时，寻求上游小浪底和小-花区间洪水过程的最可能组合，即在满足x+y=zp的条件下，求解概率密度函数f（x，y）的最大值点。将y=zp-xp代入f（x，y）中，可式（19）转化为单变量x的函数，如下式：

上式对x求导可得：

令df/dx，整理得：

式中：αx、βx、a0.x和αy、βy、a0，y分别为小浪底与小花区间洪水要素的P-III分布FX（x）和FY（y）对应参数。

对上述的方程求解，可以得到在给定的花园口站设计洪水值zp时，小-花区间洪水与小浪底的最可能的洪水地区组合（x，y）。

2.4.4 构建Copula-GA算法

对于已经构建的基于Copula 洪水地区组成模型，求得最有可能的组合情况，即在整个变量序列里求得Copula 密度函数的极值所对应的设计洪量值。但是由于Copula 函数特性，该模型求解过程繁琐，不能有效和准确的推求其全局最有可能的极值组合。因此，本文将遗传算法（Genetic Algorithms）对与构建的基于Copula 洪水地区组成模型相耦合，以探寻既客观又高效的地区组成结果［28-30］。

此处所构建的Copula-GA（Copula-Genetic Algorithms）模型，以式（23）作为遗传算法的目标函数，即适应度函数，计算个体适应度：

Copula-GA模型的运行过程可概述如下：

（1）建立Copula洪水地区组成模型；

（2）确定目标函数f（x），即率密度函数式（23）；

（3）初始化种群；

（4）随机生成初始种群；

（5）利用式（23）计算个体适应度，并保留适应度值较高的染色体；

（6）对种群交叉、变异得到下一代群体；

（7）重复（5）～（6），满足终止条件后输出适应度最大的染色体。

3 结果分析

3.1 花园口站设计洪水

本文的研究区属于半湿润区，暴雨强度大，一次暴雨历时2～3 d，最长达5 d，洪水过程历时长，可达10～12 d，因此本文构建年最大12 d 洪量W12与年最大洪峰流量Q的联合分布。Q与W12的边缘分布函数分别为FQ（q）和FW（w12）的P-III 型分布。使用线性矩法计算FQ（q）、FW（w12）的参数，成果如表1所示。

表1 花园口站洪水统计特征值成果Tab.1 The results of flood statistics characteristic value of Huayuankou Hydrological Station

由花园口站实测流量资料，计算得花园口站洪峰流量和12天洪量的Kendall 秩相关系数τ= 0.776，对应θ≈4.464。使用联合变量同频率法以及单变量同频率法对花园口站设计洪峰流量以及12天洪量进行推求，成果见图3、4以及表2所示。

图3 花园口站洪峰流量和洪量联合频率分布图Fig.3 Joint probability diagram of flood peak and flood volume at Huayuankou Station

由表2可知，采用联合变量同频率法，计算得到的各个重现期的设计洪峰流量Q和设计洪量W12大于通过单变量同频率法计算获得的设计洪峰流量Q和设计洪量W12。即使传统单变量法能够单独控制变量的重现期，但是缺乏对变量之间相关性的考虑，然而联合变量重现期所得的设计值不会低于设计标准。鉴于对防洪尤为不利的情况，取联合变量重现期作为本次推求花园口站的设计洪水特征值，将会加大花园口洪水威胁程度，而由此设计结果优化得到的调度方案会提高花园口站整体抵御洪水的能力，从而降低花园口站遭受洪灾的风险，按变倍比放大获得花园口站设计洪水过程线，结果如图5所示。

图5 花园口站设计洪水过程线Fig.5 Design flood process of Huayuankou Hydrological Station

表2 花园口水文站设计洪峰和洪量设计成果Tab.2 Design results of flood peak and flood volume of Huayuankou Hydrological Station

3.2 设计洪水地区组成

小浪底坝址、小-花区间采用1955-2015年连续还原天然状态下洪水系列资料，通过矩法估计获得小浪底水库以及小-花区间洪水统计参数，见表3。

表3 小浪底和小-花区间洪水统计参数Tab.3 Results of flood statistics in Xiaolangdi and Xiao-hua interval

小浪底坝址和小-花区间及花园口的12天洪量边缘分布分别采用P-III型分布，参数估计结果如表4所示。

表4 花园口、小浪底以及小-花区间P-III型分布参数估计结果Tab.4 Estimation results of P-III distribution parameters in Huayuankou Hydrological Station，Xiaolangdi and Xiao-hua interval

根据现有资料系列，分别采用同频率地区组成法和基于GH-Copula函数和遗传算法构建的Copula-GA 算法对小浪底与小-花区间洪水的地区组成进行计算。

（1）同频率地区组成法通过假设小浪底洪水重现期与花园口相同推求小浪底设计洪水特征值，再利用干支流关系计算小-花区间洪水特征值，计算结果如表5所示。

表5 小浪底与花园口同频率地区组成计算结果Tab.5 Flood region composition results of Xiaolangdi and Huayuankou with the same frequency

（2）基于Copula 函数建立小浪底和小-花区间联合分布，并推求其Copula 密度函数，利用干支流洪量组合关系，将在特定重现期下的Copula 密度函数转变为单变量密度函数。然后采用GA 算法寻优求解，以推求其密度函数最大极值的可能组合，并推求小浪底与小-花区间12 天洪量的对应值，得到花园口站设计洪水最可能地区组成，成果见图6以及表6。

图4 花园口站洪峰和洪量联合重现期概率图Fig.4 Probability Diagram of Joint Return Period of Flood Peak and Flood Volume at Huayuankou Hydrological Station

图6 小浪底和小-花区间联合分布函数图及小浪底洪量概率密度图（P=2%）Fig.6 Xiaolangdi and Xiao-Hua interval joint distribution function diagram and Xiaolangdi flood probability density diagram（P=2%）

通过表6 可以看出，随着花园口洪量的增大，小浪底与小-花区间洪量量级都相应的随之增大，并且小浪底洪量所占的比例明显加大。由此可以得出：小浪底的洪水对花园口洪水的发生起到了决定性作用，如果仅发生小-花区间为主的洪水，花园口一般不会出现灾害性的洪水；只有发生全流域型洪水时，即小-花区间与小浪底同时发生洪水时花园口才会出现灾害性大洪水，这与黄河中下游洪水成因及来源相印证［17，31］。

表6 花园口站设计洪水最可能地区组成成果Tab.6 The most likely flood region composition results of the design flood at Huayuankou Station

由花园口站设计洪水最可能地区组成成果可知，将小浪底在各频率地区组成所分配的洪量值，由小浪底洪量频率曲线推求出其相对应的频率值；然后再根据洪峰、洪量同频率的这一假设，分别推求出与小浪底各洪量频率对应相同的洪峰值，成果见表7。

表7 小浪底坝址设计洪水洪峰、洪量成果Tab.7 Design flood peak and flood volume results of Xiaolangdi dam site

对同频率地区组成法和Copula-GA 法计算的小浪底的洪峰流量和洪量结果分析可知：当花园口发生的洪水重现期为200年（P=1%）、100年（P=1%）以及50年（P=2%）时，小浪底最可能发生的洪水重现期约为52 a（P=1.92%）、29 a（P=3.44%）和18 a（P=5.76%）；通过假设小浪底洪水重现期与花园口相同计算得到的小浪底设计洪水，相较Copula-GA 法计算得到的最可能的设计洪水地区组成，前者所推求的小浪底洪水特征值较大，进而导致小-花区间洪水洪水特征值较小，不能达到相应的设计频率。

3.3 小浪底及小-花区间设计洪水过程

根据求得的小浪底对应花园口各设计频率的洪峰、洪量设计值，放大所选取的典型洪水，推求小浪底52年一遇、29年一遇、18年一遇设计洪水过程，成果见图7。

图7 小浪底设计洪水过程线Fig.7 Xiaolangdi design flood process

根据小浪底设计洪水过程，由马斯京根方法将小浪底洪水演进至花园口断面。通过花园口各重现期设计洪水过程线逐时刻减去小浪底演进至花园口过程线，即可得到小-花区间设计洪水过程线。

4 结论

本文使用Copula 函数构建了花园口站洪峰流量与洪量的联合分布，并计算得到了特定重现期的洪峰流量以及洪量，充分地考虑了洪峰流量与洪量之间的相关性，通过对比分析现有设计成果发现，兼顾洪峰和洪量相关性计算得到的各重现期下的洪水峰值流量和时段洪量设计值略高于现有的设计成果值［13］，更全面的考虑洪水要素之间的相互作用，也将更适用于受洪峰和洪量共同控制的区域。

在计算花园口站各重现期下的设计洪水地区组成时，针对传统的设计洪水地区组成计算方法中存在问题，本文采用Copula 函数构建了小浪底与小-花区间12 天洪量的联合分布，并与遗传算法相耦合构建了Copula-GA 模型，以快速准确的推求Copula 密度函数的最大值组合，即设计洪水最可能地区组成，得到的结果更具备客观性；当下游防洪断面将流经特定重现期洪水时，通过本研究最可能的地区组成方法所推得上游各断面洪水重现期，相比传统同频率方法，也更符合研究区洪水的实际情况。

洪水历时作为另一个洪水特征量，并与洪峰流量、洪量存在一定的相关性［21］。本文仅考虑了洪峰流量、洪量的相关性，建立了洪峰流量、洪量的联合分布，在后续研究中还应建立洪峰流量、洪量以及洪水历时间的联合分布，推求不同重现期下的洪峰流量、洪量、洪水历时组合。□