一种设计条件下双峰型洪水过程线形状推求方法与流程

本发明属于水库防洪规划设计领域，特别涉及一种设计条件下双峰型洪水过程线形状推求方法。

背景技术：

我国大中型水库的防洪安全设计一般以设计洪水过程线为依据。现行拟定设计洪水过程线方法的程序是先通过洪水频率分析得到某一指定重现期的设计洪峰流量和设计时段洪量，然后从实测资料中选择一种典型洪水过程线形状，采用同倍比或同频率方法放大推求设计洪水过程线^[1]。这就使得典型洪水过程线形状的选择结果将对水利工程建设规模及运行管理产生重要影响。

在工程实践中，通常基于“峰高量大、主峰靠后、对水库防洪偏于不利”等定性原则选取典型洪水过程线形状，理论依据不足，且选择的典型往往因人而异，任意性较大^[2-3]。事实上，洪水过程线的形状千变万化，在稀遇频率(如千年一遇)的设计条件下由于设计暴雨往往超过实测的常遇暴雨很多，其时空分布和产汇流规律与实测的典型暴雨不同，从而导致产生的设计洪水形状同典型洪水之间可能存在较大的差异^[4]。因此，推求符合设计条件下产汇流特性的洪水过程线形状并作为工程设计依据更加科学合理。

考虑到随机变量概率密度函数的形状与洪水过程线具有很大的相似性，有水文学者提出采用概率密度函数来表示洪水过程线的形状^[5]。该法采用概率密度函数对归一化的单位洪水过程线进行拟合，并可以通过调整概率分布的参数组来控制概率密度函数的形状，进而模拟多种可能的洪水过程形状类型。此外，概率密度函数法的另一个优点是将洪水过程线形状转化为易于处理的形状参数随机变量来表征，通过对形状参数变量进行频率分析就可以计算得到设计条件下的形状参数值，在此基础上推求设计洪水过程线形状，统计理论基础较强^[6]。

然而，现有技术只能提供设计条件下单峰型洪水过程线形状，存在较大的局限性，无法满足大中型水库工程设计中对双峰型洪水过程的实际需求。目前，没有文献研究了设计条件下双峰型洪水过程线形状的推求方法。

本发明涉及的参考文献如下：

[1]张波,张敏,范瑞琪,等.推求设计洪水过程线的方法,CN104182650A.

[2]李敏,薛志春,张艳平.水库典型洪水选择的可变模糊模式识别方法[J].水利科技与经济,2013,19(6):3-5.

[3]钟平安,张阳,张宇,等.一种保持典型洪水形态的设计洪水过程解析推求方法及系统,CN104598676A.

[4]葛慧,黄振平,王银堂,等.基于模糊识别理论的典型洪水选择[J].水电能源科学,2011,29(3):54-56.

[5]Yue S,Ouarda T B M J,Bobée B,et al.Approach for describing statistical properties of flood hydrograph[J].Journal of Hydrologic Engineering,2002,7(2):147-153.

[6]肖义,郭生练,方彬,等.设计洪水过程线方法研究进展与评价[J].水力发电,2006,32(7):61-63.

技术实现要素：

针对现有技术存在的不足，本发明提供了一种设计条件下双峰型洪水过程线形状推求方法。

为解决上述技术问题，本发明采用如下的技术方案：

一种设计条件下双峰型洪水过程线形状推求方法，其特征在于包括以下步骤：

步骤1，收集水库历年最大洪水过程线资料；

步骤2，根据步骤1中控制时段水库流量数据，选定历年最大洪水过程线，进而计算出历年最大洪水过程的时段洪量，对历年最大洪水过程线进行归一化得到对应的单位洪水过程线；

步骤3，根据步骤2中的单位洪水过程线，选取混合Beta分布作为混合概率密度函数线型，并估计混合概率密度函数的参数，获得洪水形状参数系列；

步骤4，对步骤3中的洪水形状参数系列进行频率分析，选取Gumbel分布作为洪水形状参数的概率分布函数线型，并估计概率分布函数的参数，计算设计条件下的洪水形状参数值；

步骤5，将步骤4中设计条件下的洪水形状参数值代入步骤3中的混合概率密度函数线型，推求设计条件下双峰型洪水过程线形状。

作为优选，所述步骤1中收集水库历年最大洪水过程线资料的时间尺度为△t，所述△t范围为0-24小时。

作为优选，所述步骤3中，估计混合概率密度函数的参数采用非线性优化算法、模拟退火算法、极大似然法、以及极大似然的EM算法和最大熵准则法的任意一种或几种组合。

作为优选，所述步骤4中，估计概率分布函数的参数采用的方法为线性矩法、权函数法或概率权重矩法。

本发明在采用混合概率密度函数对归一化的单位洪水过程线进行优化拟合得到洪水形状参数系列的基础上，对洪水形状参数系列进行频率分析，计算获得设计条件下的洪水形状参数值，进而推求设计条件下双峰型洪水过程线形状。

与现有技术相比，本发明的有益效果在于：提供了基于混合概率密度函数的设计条件下双峰型洪水过程线形状的推求方法，统计理论基础较强，更加符合设计条件下的暴雨洪水产汇流特性，更好地满足大中型水库工程设计中对双峰型洪水过程的实际需求。

附图说明

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

图2是模拟退火算法的流程图。

图3是设计条件下双峰型洪水过程线形状示意图。

具体实施方式

下面通过实施例，并结合附图对本发明作进一步说明。

如图1-图3所示，一种设计条件下双峰型洪水过程线形状推求方法，收集水库历年最大洪水过程线资料，统计得到相应的时段洪量，并归一化得到对应的单位洪水过程线，在采用混合概率密度函数对归一化的单位洪水过程线进行优化拟合得到洪水形状参数系列的基础上，对洪水形状参数系列进行频率分析，计算获得设计条件下的洪水形状参数值，进而推求设计条件下双峰型洪水过程线形状。图1是本实施例的计算流程图，按照以下步骤进行：

1.收集水库历年最大洪水过程线资料。

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

2.归一化得到对应的单位洪水过程线。

根据步骤1中控制时段水库流量数据，选定历年最大洪水过程线，进而计算出历年最大洪水过程线的时段洪量，本具体实施中控制时段为D＝m·△t，有n年最大洪水过程线系列，采用复化梯形公式计算第j年的时段洪量V_j(j＝1,2,…,n)：

其中，q_ij为第j年，第i时刻(i＝0,1,…,m)的洪水流量，i、j、n和m均为自然数。

由于混合Beta概率密度函数的定义域为[0，1]，且纵坐标概率密度函数值与横轴包围的面积为1个单位，为了使历年最大洪水过程线满足混合Beta概率密度函数这些统计特性，需要将历年最大洪水过程线进行归一化得到对应的单位洪水过程线。

本具体实施中将横坐标除以控制时段总历时，纵坐标乘以控制时段总历时与时段洪量之比，计算公式分别如下：

其中，x_i为归一化的横坐标，f_ij为归一化的第j年，第i时刻(i＝0,1,…,m)的纵坐标洪水流量。

3.确定混合概率密度函数线型和参数，得到洪水形状参数系列。

本步骤包括两个子步骤：

3.1选择混合概率密度函数线型

由于单位洪水过程线的概率密度线型是未知的，通常选用能较好拟合实测单位洪水过程线的线型。经过水文学者的比较与研究，发现Beta概率密度函数适应性强，拟合使用灵活，能够给出多种可能的形状类型。

本具体实施中采用由2个Beta概率密度函数加权组成的混合Beta概率密度函数作为混合概率密度函数线型，表达式为

f(x)＝αf₁(x)+(1-α)f₂(x)(0≤x≤1) (4)

其中，f(x)为混合Beta概率密度函数；α为权重系数；f₁(x)、f₂(x)分别为Beta概率密度函数，表达式分别为

其中，B(a,b)为Beta函数，计算公式如下

其中，a₁和b₁、a₂和b₂分别为Beta函数的参数。

3.2估计混合概率密度函数线型的参数

在子步骤3.1中选定混合Beta概率密度函数线型后，接着需要进行估计混合Beta概率密度函数线型的参数。每年的单位洪水过程线采用一个混合Beta概率密度函数进行表示和模拟，每个混合Beta概率密度函数共有α、a₁、b₁、a₂和b₂等5个参数需要进行估计。

混合分布的参数估计方法主要有非线性优化算法、极大似然法、极大似然的EM算法和最大熵准则法等。模拟退火算法是一种有效的全局优化算法，在组合优化问题的求解中，与以往的近似算法相比，其具有描述简单、使用灵活、运用广泛、运行效率高及受初始条件限制较少等特点。模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率l收敛于全局最优解的全局优化算法。

本具体实施中采用模拟退火算法估计混合Beta概率密度函数的参数。寻求最优化参数的目标函数为

其中，E_j表示第j年的目标函数值；f_ij^、f_ij分别为第j年，第i时刻模拟和实测的单位洪水过程线流量值。

如图2所示，给出了模拟退火算法的流程图。基本步骤如下：

(1)初始化：初始温度T(充分大)，初始解状态S(算法迭代的起点)，每个T值的迭代次数L；

(2)对k＝1,…,L做第3至第6步：

(3)产生新解S′；

(4)计算增量ΔT＝C(S′)-C(S)，其中C(S)为目标函数；

(5)若ΔT<0则接受S′作为新的当前解，否则按Metropolis准则接受新解，即以概率exp(-ΔT/T)接受S′作为新的当前解；

(6)如果满足终止条件则输出当前解作为最优解，结束程序；终止条件通常取为连续若干个新解都没有被接受时终止算法；

(7)T逐渐减少，且T＝T-1且T>0，然后转第2步。

对历年的单位洪水过程线，分别估计其混合Beta概率密度函数的参数，就可以获得样本容量为n的洪水形状参数系列：{α_j}，{a_1j}，{b_1j}，{a_2j}和{b_2j}(j＝1,2,…,n)。

4.对洪水形状参数系列进行频率分析，计算设计条件下的洪水形状参数值。

根据对步骤3中的洪水形状参数系列进行频率分析，选取适当的概率分布函数线型，并估计概率分布函数的参数，计算设计条件下的洪水形状参数值，本步骤包括三个子步骤：

4.1选择概率分布函数线型

由于洪水形状参数的总体分布频率线型是未知的，通常选用能较好拟合洪水形状参数资料系列的线型。经过研究发现，Gumbel分布对于洪水形状参数资料拟合较好。

本具体实施中采用Gumbel分布作为洪水形状参数的概率分布函数线型，其概率分布函数表达式为

其中μ和β为Gumbel分布函数的参数。

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

在子步骤4.1中采用Gumbel分布作为洪水形状参数线型后，接着需要进行估计Gumbel分布的参数。线性矩法是目前国内外公认的有效参数估计方法，最大特点是对序列的极大值和极小值没有常规矩那么敏感，求得的参数估计值比较稳健。

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

4.3计算设计条件下的洪水形状参数值

经过子步骤4.1和4.2后就可以得到对应的洪水形状参数的概率分布函数曲线。对于指定的设计重现期T，设计条件下的洪水形状参数值由下式计算：

y_T＝μ-βln(lnT) (10)

5.推求设计条件下双峰型洪水过程线形状。

将步骤4中设计条件下的洪水形状参数值α_T，a_1T，b_1T，a_2T和b_2T代入步骤3中的混合Beta概率密度函数式(4)、(5)和(6)，推求设计条件下双峰型洪水过程线形状。

设计条件下混合Beta概率密度函数表达式为

f(x)＝α_Tf₁(x)+(1-α_T)f₂(x)(0≤x≤1) (11)

如图3所示，给出了设计条件下双峰型洪水过程线形状示意图。

综上，本发明通过收集水库历年最大洪水过程线资料，统计得到相应的时段洪量，并归一化得到对应的单位洪水过程线，在采用混合概率密度函数对归一化的单位洪水过程线进行优化拟合得到洪水形状参数系列的基础上，对洪水形状参数系列进行频率分析，计算获得设计条件下的洪水形状参数值，进而推求设计条件下双峰型洪水过程线形状。本发明提供了基于混合概率密度函数的设计条件下双峰型洪水过程线形状的推求方法，统计理论基础较强，更加符合设计条件下的暴雨洪水产汇流特性，更好地满足大中型水库工程设计中对双峰型洪水过程的实际需求。