CN105887751A - 一种考虑水库调洪演算的多变量洪水重现期计算方法 - Google Patents

Abstract

本发明公开了一种考虑水库调洪演算的多变量洪水重现期计算方法，在利用遗传规划进行符号回归推求水库调洪函数的显式表达式，并采用Copula函数求解联合概率密度函数的基础上，对联合概率密度函数在安全区域内进行积分运算获取指定年最高水位的累积概率，从而计算指定年最高水位的单变量重现期，作为考虑水库调洪演算的多变量洪水重现期。本发明考虑了水库调洪演算，以年最高水位等值线定义危险区域，能够正确表征水库大坝的水文失事机理，适用于具有任意复杂程度的洪水调度规则，为水库防洪安全设计提供更加科学合理的依据。

Description

一种考虑水库调洪演算的多变量洪水重现期计算方法

技术领域

本发明涉及防洪安全技术领域，具体地说是一种考虑水库调洪演算的多变量洪水重现期计算方法。

背景技术

洪水灾害是世界上最严重的自然灾害之一，如何进行科学有效的洪水风险管理以减少损失具有重大意义。洪水重现期是指某一特定量级的洪水事件连续两次发生的平均时间间隔(以年计)，常用来衡量洪水大小和等级，是确定水利水电工程规模的主要依据(参见郭生练等，一种梯级水库下游设计洪水的推求方法，专利号ZL2013101982328)。传统的洪水重现期往往只基于某一特征变量(洪峰或者某一时段洪量)进行单变量频率计算得到。然而，洪水事件本质是一种多变量随机水文事件，需要同时考虑多个特征变量才能完整描述，导致单变量洪水重现期无法真实反映洪水事件的量级(参见张波等，推求设计洪水过程线的方法，公开号CN104182650A)。近年来，多变量洪水重现期越来越多地被国内外水文工作者应用于洪水事件的描述和建模中。

多变量洪水重现期的计算方法由于涉及到多个洪水变量的组合，因而比单变量情形更加复杂。目前使用最为广泛的多变量洪水重现期主要有“或”和“且”重现期，“或”重现期认为只要洪峰或洪量其中一个超过设定的阈值就构成危险事件，“且”重现期则需要洪峰和洪量同时都超过设定的阈值(参见黄强等，基于二次重现期的多变量洪水风险评估[J].湖泊科学,2015,27(2):352-360)。一般来说，防洪标准越高，相应的危险率就越小，指定防御标准下水工建筑物遭遇破坏的风险就越低。这两种重现期均根据联合概率等值线上的某一特定联合设计值定义危险事件，由此导致的多变量重现期与危险区域范围大小的矛盾会造成对安全事件与危险事件的错误识别，局限性较为明显(参见史黎翔等，基于Copula函数的两变量洪水重现期与设计值计算研究[J].水力发电学报,2015,34(10):27-34)。为了解决上述问题，Salvadori等根据联合概率等值线定义危险区域提出了Kendall重现期计算方法，认为只有当洪峰洪量联合概率值超过了设定的阈值时才是危险事件(参见Salvadori G,Michele C D,DuranteF.On the return period and design in a multivariate framework[J].Hydrology and Earth System Sciences,2011,15(11):3293-3305)。Kendall重现期使得洪峰、洪量联合概率等值线上的任意一点都具有相同的危险区域，对安全事件与危险事件的判定更加合理。

事实上，就水库防洪安全而言，最重要的因素是坝前最高水位(参见丁晶等.水库防洪安全设计时设计洪水过程线法适用性的探讨[J].水科学进展,1992,3(1):45-52)。前面三种多变量重现期所判定的危险事件不符合实际情况，真正的危险事件应该是坝前最高水位超过某一指定高程。因此，现有的多变量重现期计算方法没有考虑水库调洪演算，无法正确表征水库大坝的水文失事机理。

发明内容

本发明的目的是克服现有技术存在的不足，提供一种考虑水库调洪演算的多变量洪水重现期计算方法。

本发明一种考虑水库调洪演算的多变量洪水重现期计算方法，包括如下步骤：

步骤1，收集水库历年最大洪水过程资料和洪峰、洪量系列；

步骤2，将步骤1中的历年最大洪水过程输入水库按照调度规则进行调洪演算，从输出的库水位过程线统计得到年最高水位系列；

步骤3，根据步骤1中的洪峰、洪量系列和步骤2中得到的年最高水位系列，利用遗传规划进行符号回归，推求水库调洪函数的显式表达式，并确定指定年最高水位值对应的安全区域；

步骤4，根据步骤1中的洪峰、洪量系列，选取边缘概率分布函数线型，并估计边缘概率分布函数的参数；

步骤5，采用Copula函数构造洪峰、洪量的联合概率分布函数，并估计Copula函数的参数，求解联合概率密度函数；

步骤6，对步骤5得到的联合概率密度函数在步骤3中确定的安全区域内进行积分运算，获取指定年最高水位的累积概率；

步骤7，根据步骤6中得到的累积概率计算指定年最高水位的单变量重现期，作为考虑水库调洪演算的多变量洪水重现期。

所述步骤2中，水库调洪演算采用四阶龙格-库塔算法进行数值求解。

所述步骤4中，将P-III型分布作为洪峰和洪量的边缘概率分布函数线型，采用线性矩法估计边缘概率分布函数的参数。

所述步骤5中，采用Gumbel-Hougaard Copula函数构造洪峰和洪量的联合概率分布函数，采用Kendall秩相关性系数法估计Gumbel-HougaardCopula函数的参数。

所述步骤6中，利用蒙特卡洛模拟方法计算定积分的值。

本发明在利用遗传规划进行符号回归推求水库调洪函数的显式表达式，并采用Copula函数求解联合概率密度函数的基础上，对联合概率密度函数在安全区域内进行积分运算获取指定年最高水位的累积概率，从而计算指定年最高水位的单变量重现期，作为考虑水库调洪演算的多变量洪水重现期。

与现有技术相比，本发明方法的有益效果在于：

1、本发明考虑了水库调洪演算，以年最高水位等值线定义危险区域，能够正确表征水库大坝的水文失事机理；

2、本发明利用遗传规划进行符号回归，可以自动得到回归曲线的函数形式及其参数，适用于具有任意复杂程度的水库洪水调度规则，为挖掘水库调洪函数的显式表达式提供了有力工具。

附图说明

图1为本发明方法的流程图。

图2为利用遗传规划进行符号回归的流程图。

图3为指定年最高水位值对应的安全区域示意图。

具体实施方式

下面结合附图对本发明进行进一步说明。

如图1-图3所示，一种考虑水库调洪演算的多变量洪水重现期计算方法，收集水库历年最大洪水过程资料和洪峰、洪量系列，输入水库按照调度规则进行调洪演算统计得到年最高水位系列，在利用遗传规划进行符号回归推求水库调洪函数的显式表达式，并采用Copula函数求解联合概率密度函数的基础上，对联合概率密度函数在安全区域内进行积分运算获取指定年最高水位的累积概率，从而计算指定年最高水位的单变量重现期，作为考虑水库调洪演算的多变量洪水重现期。图1是本实施例的计算流程图，按照以下步骤进行：

1.收集水库历年最大洪水过程资料和洪峰、洪量系列。

收集的水库流量数据的时间尺度为Δt(一般宜小于日，如1小时，3小时等)。历年最大洪水过程的控制时段根据流域暴雨洪水特性和水库调洪能力确定。本具体实施中洪峰系列由历年的最大洪水流量组成，洪量系列由历年的最大控制时段洪水总量组成。

2.将历年最大洪水过程输入水库调洪演算，得到年最高水位系列。

水库调洪演算的基本依据是水量平衡方程和槽蓄动力方程：

$I_t-O_t=\frac{{\mathrm{dS}}_t}{dt}---\left(1\right)$

O_t＝f(S_t) (2)

式中：I_t、O_t分别为t时刻的入库流量和出库流量(m³/s)；S_t为t时刻的水库蓄水量(m³)；函数f(S_t)依据水库库容曲线和泄流曲线综合确定。

目前调洪演算求解方法主要有列表试算法、图解法和数值解法等。其中，数值解法既不必试算，也不需绘制辅助曲线，不仅可方便地变动计算步长，且能简捷地进行变化闸门开度的计算，尤其对复杂情况下水库的调洪计算，较其他方法有一定的优势。

本具体实施中水库调洪演算采用四阶龙格-库塔算法进行数值求解。

将步骤1中的历年最大洪水过程输入水库按照调度规则进行调洪演算，从汛期防洪限制水位起调，采用四阶龙格-库塔算法进行数值求解，从输出的库水位过程线统计得到年最高水位系列。

3.利用遗传规划进行符号回归，推求水库调洪函数的显示表达式，确定安全区域。

本步骤包括两个子步骤：

3.1推求水库调洪函数的显示表达式

遗传规划是在遗传算法的基础上提出的一种新的进化计算方法，具有强大的启发式自动搜索寻优能力，其优点是能自动找到因变量与自变量间的变化规律，无需事先确定各变量之间的函数关系，简单易行，因此在数据挖掘、自动控制、发现工程经验公式、微分方程求解、符号回归等方面得到了广泛的应用。

遗传规划的基本思想是根据达尔文“适者生存”的原则对问题进行结构化处理，随机产生一个适用于问题环境的初始种群，计算群体中每个个体的适应度，采用复制、杂交和变异等遗传算子对种群不断进行迭代优化，直到问题的最优解在某一代出现。

令H、Q和W分别表示年最高水位、洪峰和洪量，相应的取值为h、q和w.本具体实施根据步骤1中的洪峰、洪量系列和步骤2中得到的年最高水位系列，利用遗传规划进行符号回归，推求水库调洪函数的显式表达式h＝g(q,w)。如图2所示，给出了利用遗传规划进行符号回归的流程图。基本步骤如下：

(1)确定个体表达式结构，包括F(函数集)和T(终止符集)。本具体实施中F＝{+,-,×,/,√,log,exp}，T＝{q,w}。

(2)产生初始群体。本具体实施中利用混合法随机生成，生长法与完全法各50％，从函数集F及终止符集T中选择字符组成不同的公式作为初始个体。本具体实施中设定种群数量为M。

(3)计算个体适应度。本具体实施中采用均方根误差作为适应度函数对个体的好坏进行评判，适应度的值越小，该个体越好。

(4)生成新一代群体。执行遗传操作，生成新的个体，主要遗传操作包括：①复制，将已有的优良个体复制，加入新群体，并相应删除劣质个体；②交换，将选出的两个个体的部分节点进行交换，所产生的两个新个体加入新群体；③突变，随机改变个体的某一部分，将新个体插入新群体中。本具体实施中复制概率P_r，交叉概率为P_c，变异概率P_m。

(5)重复(3)、(4)，直至满足终止条件(达到最大迭代代数或最佳个体适应度达到预定值)，选择最好结果作为最终解。本具体实施中G代表迭代次数，初始种群为第0代，终止准则为达到最大迭代代数G_max。

3.2确定安全区域

在子步骤3.1中确定了水库调洪函数后，给定指定年最高水位值h，以年最高水位等值线h＝g(q,w)判定安全事件和危险事件，满足g(q,w)<h的洪峰、洪量组合(q,w)为安全事件。

满足g(q,w)<h所有点(q,w)组成的集合即为指定年最高水位值h相应的安全区域为D_h＝{(q,w):g(q,w)<h}。如图3所示，给出了指定年最高水位值对应的安全区域示意图。

4.确定洪峰和洪量的边缘概率分布函数。

根据步骤1中的洪峰、洪量数据资料，选取适当的边缘概率分布函数线型，并估计边缘概率分布函数的参数，本步骤包括两个子步骤：

4.1选择边缘概率分布函数线型

由于洪峰和洪量的总体分布频率线型是未知的，通常选用能较好拟合多数洪水样本资料系列的线型。经过多年分析比较与研究，发现P-III型分布对于我国大部分河流的洪水资料拟合较好，推荐在工程实践中采用。

本具体实施中采用P-III型分布作为洪峰和洪量的边缘概率分布函数线型。

4.2估计边缘概率分布函数线型的参数

在子步骤4.1中采用P-III型分布作为洪峰和洪量边缘概率分布函数线型后，接着需要进行估计P-III型分布的参数。线性矩法(L-矩法)是目前国内外公认的有效参数估计方法，最大特点是对序列的极大值和极小值没有常规矩那么敏感，求得的参数估计值比较稳健。

本具体实施中采用L-矩法估计边缘概率分布函数线型的参数。

5.采用Copula函数构造洪峰、洪量的联合概率分布函数，求解联合概率密度函数。

根据步骤1中的洪峰、洪量数据资料以及步骤4中估计的边缘概率分布函数，选取适当的Copula函数作为连接函数构造洪峰和洪量的联合概率分布函数，并估计其参数，在此基础上求解联合概率密度函数，本步骤包括三个子步骤：

5.1选择Copula函数

假设F_Q(q)、F_W(w)分别表示Q和W的边缘概率分布函数，相应的概率密度函数为f_Q(q)、f_W(w)。由Sklar定理可知，Q,W的联合概率分布函数可以用一个二维Copula函数C表示：

F_Q,W(q,w)＝C_θ(F_Q(q),F_W(w))＝C_θ(u,v) (3)

其中，θ为Copula函数的参数；u＝F_Q(q),v＝F_W(w)为边缘概率分布函数。

本具体实施中，采用Gumbel-Hougaard Copula函数构造洪峰和洪量的联合概率分布函数，其表达式如下：

$C_{\mathrm{\&theta;}}(u,v)=\exp \{-{\&lsqb;{(-lnu)}^{\mathrm{\&theta;}}+{(-lnv)}^{\mathrm{\&theta;}}\&rsqb;}^{\frac{1}{\mathrm{\&theta;}}}\},\mathrm{\&theta;}\&GreaterEqual;1---\left(4\right)$

5.2估计Copula函数的参数

本具体实施中，采用Kendall秩相关性系数法估计Gumbel-HougaardCopula函数的参数。Kendall相关系数τ与参数θ的关系为：

$\mathrm{\&theta;}=\frac{1}{1-\mathrm{\&tau;}}---\left(5\right)$

令{(x₁,y₁),…,(x_n,y_n)}表示从连续随机变量(X,Y)中抽取的n个观测值的随机样本，则在样本中有种不同的观测值组合(x_i,y_i)和(x_j,y_j)。样本Kendall秩相关系数τ通过下式计算

$\mathrm{\&tau;}=\frac{2}{n(n-1)}\underset{i=1}{\overset{n-1}{\&Sigma;}}\underset{j=i+1}{\overset{n}{\&Sigma;}}sign\&lsqb;(x_i-x_j)(y_i-y_j)\&rsqb;---\left(6\right)$

其中，sign(·)是符号函数。

5.3求解联合概率密度函数

在子步骤5.1和5.2中得到洪峰和洪量的联合概率分布函数后，就可以通过式(7)求解相应的联合概率密度函数：

f_Q,W(q,w)＝c_θ(u,v)f_Q(q)f_W(w) (7)

其中，c_θ(u,v)为Copula函数的密度函数，解析表达式为：

$c_{\mathrm{\&theta;}}(u,v)=C_{\mathrm{\&theta;}}(u,v)\frac{{\left(\ln u\ln v\right)}^{\mathrm{\&theta;}-1}}{uv}\{(\mathrm{\&theta;}-1){\&lsqb;-{\mathrm{lnC}}_{\mathrm{\&theta;}}(u,v)\&rsqb;}^{1-2\mathrm{\&theta;}}+{\&lsqb;-{\mathrm{lnC}}_{\mathrm{\&theta;}}(u,v)\&rsqb;}^{2-2\mathrm{\&theta;}}\}---\left(8\right)$

6.对联合概率密度函数在安全区域内进行积分运算，获取指定年最高水位的累积概率。

对步骤5得到的联合概率密度函数f_Q,W(q,w)在步骤3中确定的安全区域D_h内进行积分运算，通过式(9)获取指定年最高水位h的累积概率：

$F_H\left(h\right)=\underset{D_h:g(q,w)<h}{\&Integral;\&Integral;}{f}_{Q,W}(q,w)dqdw---\left(9\right)$

实际计算时，考虑到联合概率密度函数f_Q,W(q,w)和安全区域D_h都比较复杂，因而难以得到式(9)的解析解。本具体实施中利用蒙特卡洛模拟方法计算定积分的值。

7.计算考虑水库调洪演算的多变量洪水重现期。

根据步骤6中得到的累积概率F_H(h)，通过式(10)计算指定年最高水位h的单变量重现期：

$T_R\left(h\right)=\frac{1}{1-F_H\left(h\right)}---\left(10\right)$

T_R(h)即为所有满足h＝g(q,w)的洪峰、洪量组合(q,w)所对应的考虑水库调洪演算的多变量洪水重现期。

综上，本发明通过收集水库历年最大洪水过程资料和洪峰、洪量系列，输入水库按照调度规则进行调洪演算统计得到年最高水位系列，在利用遗传规划进行符号回归推求水库调洪函数的显式表达式，并采用Copula函数求解联合概率密度函数的基础上，对联合概率密度函数在安全区域内进行积分运算获取指定年最高水位的累积概率，从而计算指定年最高水位的单变量重现期，作为考虑水库调洪演算的多变量洪水重现期。本发明考虑了水库调洪演算，以年最高水位等值线定义危险区域，能够正确表征水库大坝的水文失事机理，适用于具有任意复杂程度的洪水调度规则，为水库防洪安全设计提供更加科学合理的依据。

Claims (5)

1.一种考虑水库调洪演算的多变量洪水重现期计算方法，其特征在于包括如下步骤：

步骤1，收集水库历年最大洪水过程资料和洪峰、洪量系列；

步骤2，将步骤1中的历年最大洪水过程输入水库按照调度规则进行调洪演算，从输出的库水位过程线统计得到年最高水位系列；

步骤3，根据步骤1中的洪峰、洪量系列和步骤2中得到的年最高水位系列，利用遗传规划进行符号回归，推求水库调洪函数的显式表达式，并确定指定年最高水位值对应的安全区域；

步骤4，根据步骤1中的洪峰、洪量系列，选取边缘概率分布函数线型，并估计边缘概率分布函数的参数；

步骤5，采用Copula函数构造洪峰、洪量的联合概率分布函数，并估计Copula函数的参数，求解联合概率密度函数；

步骤6，对步骤5得到的联合概率密度函数在步骤3中确定的安全区域内进行积分运算，获取指定年最高水位的累积概率；

步骤7，根据步骤6中得到的累积概率计算指定年最高水位的单变量重现期，作为考虑水库调洪演算的多变量洪水重现期。

2.如权利要求1所述的方法，其特征在于：所述步骤2中，水库调洪演算采用四阶龙格-库塔算法进行数值求解。

3.如权利要求1所述的方法，其特征在于：所述步骤4中，将P-III型分布作为洪峰和洪量的边缘概率分布函数线型，采用线性矩法估计边缘概率分布函数的参数。

4.如权利要求1所述的方法，其特征在于：所述步骤5中，采用Gumbel-Hougaard Copula函数构造洪峰和洪量的联合概率分布函数，采用Kendall秩相关性系数法估计Gumbel-Hougaard Copula函数的参数。

5.如权利要求1所述的方法，其特征在于：所述步骤6中，利用蒙特卡洛模拟方法计算定积分的值。