本发明涉及洪水频率分析领域，具体是一种考虑历史洪水的洪水频率估计方法。
背景技术：
在各种水利工程的规划设计中，都需要洪水频率分析计算工作，为工程规模和建筑物尺寸的确定提供依据。常采用频率分析方法进行计算，该方法对样本数据具有很高的敏感性。然而，在我国多数水文站始建于20世纪50年代后，实测洪水序列较短，即使通过插补延长，资料长度一般仅有40-60年。采用仅有的几十年长度的极值样本估计百年一遇、千年一遇设计洪水时产生的误差较大，不利于水利工程规模和建筑物尺寸的确定，可能会给工程带来风险。因此，从有限的实测资料中获取更多洪水信息是提高洪水频率计算精度的关键。
对实测资料的利用上，一般采用历史洪水、古洪水和区域洪水等信息插补或延长样本序列，提高分析结果的代表性。众多研究表明，在实测洪水资料中加入历史洪水构成不连续样本，能有效提高洪水频率分析精度、稳定性和合理性。然而，历史洪水信息记录并不详尽，仅有部分历史洪水信息记录相对完善，仍有部分历史洪水信息相关记录相对较少，导致其洪峰流量难以定量估计，也有部分历史洪水通过调查、分析和测量后仅知其量级或彼此之间的位置关系，无法估计洪峰流量具体数值，此类历史洪水信息研究人员称之为不定量历史洪水。当前相关研究中对不定量历史洪水对频率分析的影响研究较少，根据当前署名发明人的工作发现，现有技术中仅将定量或不定量历史洪水信息用于洪水频率分析，难以同时将两类历史洪水信息纳入样本序列用于频率分析计算，给参数估计带来很大困难，洪水频率分析结果的代表性和准确性仍然存在提高的空间。
技术实现要素：
发明目的：针对现有技术的不足，本发明提出了一种考虑历史洪水的洪水频率估计方法，本方法即可考虑定量历史洪水信息，也可考虑不定量的历史洪水信息，可提高洪水频率分析的样本代表性，以及参数估计的准确性。
技术方案：一种考虑历史洪水的洪水频率估计方法，包括如下步骤：
s1、根据流域特性或规范要求，选取合适的频率线型，确定待估参数；
s2、收集水文站监测的流量资料，构建实测洪峰样本序列，并计算确定待估参数的上、下边界，所述实测洪峰样本序列为由s个年最大洪峰流量数值构成的序列；
s3、收集历史洪水信息，构建历史洪水序列，包括定量和不定量历史洪水信息，所述定量历史洪水信息由k1个确定的洪峰流量数值构成的序列；所述不定量历史洪水信息由k2个不定量历史洪水的流量区间构成的序列；
s4、将历史洪水序列与实测洪水序列纳入计算样本，构建考虑定量及不定量历史洪水信息的参数估计似然函数；
s4、将似然函数作为目标函数，待估参数边界作为约束条件，建立优化问题，并采用遗传算法进行求解，进而得到频率线型的待估参数；
s5、根据选定的线型及估计参数，计算不同重现期的设计洪水。
本发明以《水利水电工程设计洪水计算规范》(sl44-2006)推荐的皮尔逊ⅲ型为例进行说明的。所述步骤s1中选用皮尔逊ⅲ型作为频率线型，确定待估参数θ，包括均值
有益效果：
1、本方法能够将定量及不定量历史洪水信息均纳入样本序列，可增加样本的容量和数量，有助于改善样本序列的代表性，以提高参数估计和设计洪水的准确性。
2、本方法能够有效利用定量及不定量历史洪水信息，在洪水频率分析中具有更为广泛的适用性。
3、采用遗传算法求解似然函数，可快速对参数进行寻优，具有更高的计算效率。
附图说明
图1是根据本发明实施例的洪水频率估计方法流程图。
具体实施方式
下面结合附图对本发明的技术方案作更进一步的说明。
本发明公开了一种考虑历史洪水的洪水频率估计方法，能够将不定量及定量历史洪水信息均纳入样本序列，增加样本序列容量或数量，改善样本序列的代表性，并用遗传算法进行参数快速寻优，可有效提高参数估计和设计洪水的准确性。参照图1，方法包括以下步骤：
(s1)根据流域特性或规范要求，选取合适的线型。
本发明选择《水利水电工程设计洪水计算规范》(sl44-2006)推荐的皮尔逊ⅲ型作为频率线型，确定待估参数θ，包括均值
式中，f(·)为皮尔逊ⅲ型分布概率密度函数，α为形状参数，β为尺度参数，a0为位置参数，γ(·)为伽马函数。
(s2)建立实测样本序列。
收集流域的水文站实测流量资料，采用年最大值法选取样本，组成实测洪峰样本序列，对其进行可靠性、一致性和代表性“三性”审查。本发明仅以洪峰流量序列为例进行说明，应当理解也可适用于洪量频率计算。采用矩法计算实测洪峰样本序列待估参数的抽样误差，以此来确定待估参数的上、下边界。
式中，s为实测样本序列个数；xi为实测洪峰样本序列，即洪峰流量构建的序列，每年一个值，监测多少年，就有多少个数；
式中，
用正态分布描述参数估计的抽样误差，以实测样本资料计算的均值和方差确定待估参数的上下边界：
式中，
(s3)获取历史洪水信息，包括定量和不定量历史洪水信息。
搜集本流域的建站以前关于历史洪水相关资料，包括文献档案资料、史志、碑记、传说、洪痕等，采用水位流量关系法、比降法、水面曲线法等已有方法估计历史洪水洪峰流量的大小，并检查其合理性。若已有历史洪水资料调查整编成果，可直接采用。定量历史洪水信息包含历史上已发生过的较大洪水，可提高样本序列的代表性。当根据现有方法难以直接估计出某一个具体数值时，我们称之为不定量历史洪水，可以根据现有信息、位置关系、大于或小于某一洪水(如：某年洪水超过已知的某一年，而低于另外某一年的)等资料，确定其上下边界(yl，yu)，构成不定量历史洪水信息。也就是，用这个上下边界表示该历史洪水信息，就是洪峰流量的具体范围，也纳入计算样本中用于分析计算，但与实测样本序列不一样，在似然函数中单独计算。
(s4)根据已收集的调查资料，分析洪水发生的次数和量级，合理确定重现期h。
根据历史文献记录中有关洪水淹没造成的破坏程度、情节等，与已有的文字描述及定量的调查洪水对比，可以分析各次洪水的量级范围和大小序位，合理确定计算序列中历史洪水的重现期。
(s5)将历史洪水与实测洪水信息纳入计算样本，构建考虑定量及不定量历史洪水信息的参数估计似然函数。
式中，l(d|θ)为参数估计似然函数，ls(d|θ)ls为实测洪水估计参数似然函数，lh(d|θ)lh为历史洪水似然函数，二者乘积作为最终参数估计似然函数，f(·)是皮尔逊ⅲ型概率密度函数；k为历史洪水个数，其中k1为可定量洪水个数，k2不定量历史洪水个数；x0为某一洪水标准(或称阈值)；f(·)为皮尔逊ⅲ型分布累计密度函数；h(＝h-s)为调查洪水历史时期长度。
(s6)建立优化问题并求解。
将似然函数作为目标函数，s2中通过实测样本计算得到的待估参数边界作为约束条件，采用遗传算法进行寻优，得到频率线型的参数，结合洪水的特性进行合理性分析。
具体地，将前文建立的似然函数作为遗传算法的目标函数；将前文的实测洪峰流量、定量历史洪水信息和不定量历史洪水信息都代入到似然函数中进行计算。根据待估参数的上、下边界，随机生成初始群体，计算群体中各个个体的适应度，选择群体中适应度高的个体，直接遗传至下一代，通过配对交叉产生新的个体再遗传至下一代，并对群体中的个体进行变异产生新一代，进而得到下一代群体，以进化过程中所得到的具有最大适应度个体作为最优解输出，终止计算，得到待估参数值。合理性分析是洪水频率计算的通用步骤，不管用何种方法计算得到的结果，均进行合理性分析，一般根据统计参数和设计洪水值的变化规律，以及上下游、干支流和邻近流域各水文站的成果进行合理性检查。
(s7)根据选定的线型及估计参数，计算不同重现期(t)的设计洪水。
根据上述方法估计的均值
xp＝β·g-1(1-p,α,1) a0
式中，xp为重现期t的设计洪水，频率
下面通过两个实例验证本发明所提方法的效果。
(1)模拟序列
采用蒙特卡洛方法生成序列长度为50年年最大洪峰序列。模拟序列服从p-ⅲ型分布，其统计参数分别为：均值
采用概率权重矩法(宋德敦，丁晶.概率权重矩法及其在p-ⅲ分布中的应用[j].水利学报，1988，19(3)：1-11)、stedinger和cohn提出的方法(stedingerjr,cohnta.floodfrequency-analysiswithhistoricalandpaleofloodinformation[j].waterresourcesresearch,1986,22(5):785-793)，以及本发明方法对模拟序列进行参数估计，并计算百年一遇(t-100)和千年一遇(t-1000)设计洪水，具体见表1。其中概率权重矩法、stedinger和cohn提出的方法均为确定性方法，仅定量历史洪水信息能够被利用，并以此进行参数估计与设计洪水计算。
表1模拟序列几种方法参数估计结果和设计洪水对比表
据表1，相对于概率权重矩法、stedinger和cohn提出的方法，本发明方法能够利用更多的洪水信息资料，待估参数与设计洪水估计值均与真值之间最为接近，相对误差也最小，其中待估参数(
概率权重矩法估计的参数和设计洪水误差最大，其中计算的千年一遇设计洪水小于百年一遇洪水的真值，stedinger和cohn提出的方法计算的设计洪水与真值相比，其相对误差均超过10％，并随着重现期的增加，相对误差的增速也加快，进而影响到参数估计与设计洪水的代表性和准确性，不利于工程规模的确定。
考虑到现阶段，我国实测洪水序列较短，资料长度一般仅有40-60年，相对于设计洪水重现期而言，样本数量相对较少，其代表性也存在一定的不确定性，定量与不定量历史洪水信息可有效改善样本的代表性，利用该资料可提高参数估计与设计洪水的准确性和代表性。
(2)实测序列
采用宜昌站水文资料对发明方法适用情况进行分析。宜昌站自1877年设立海关水尺以来，至2018年已有142年实测水文资料。在该河段，1153-1877年的725年间，共调查发生了23次较大洪水，仅有8次洪水可定量估算，其流量大于80000m3/s，发生的年份为1870年、1227年、1560年、1153年、1860年、1788年、1796年和1613年，对应的洪峰流量分别为105000m3/s、96300m3/s、93600m3/s、92800m3/s、92500m3/s、86000m3/s、82200m3/s、81000m3/s。洪峰流量超过71000m3/s(大于宜昌实测系列的首项洪水1896年)，但小于80000m3/s，发生的年份有1520年、1613年、1700年、1761年及1840年等5年，其他10个年份(1310年、1463年、1478年、1513年、1520年、1550年、1574年、1658年、1672年、1681年、1859年)的历史洪水仅能定性为超过1896年洪水(洪峰流量为71100m3/s)。
三峡水利枢纽设计依据的宜昌站1877-1990年114实测序列、8场定量历史洪水资料。为此，本发明也采用相同资料进行计算，并增加15场不定量历史洪水资料，但彼此之间位置关系难以明确，按照前述介绍，此15场不定量历史洪水其上边界均定位80000m3/s，下边界定位71100m3/s。调查期1153-1877年的725年内，23场次定量及不定量历史洪水信息均大于宜昌实测系列的首项洪水1896年，故将阈值x0确定为71100m3/s。
采用概率权重矩法、stedinger和cohn提出的方法，以及本发明方法对宜昌站洪峰序列资料进行参数估计，并计算百年一遇(t-100)、两百年一遇(t-200)、五百年一遇(t-500)、千年一遇(t-1000)、五千年一遇(t-5000)和万年一遇(t-10000)设计洪水，具体见表2。其中率权重矩法、stedinger和cohn提出的方法仅能利用确定性资料，共包括8场定量历史洪水资料及1877-1990年114年实测序列，并以此进行参数估计与设计洪水计算。本发明方法在此基础上，将15场难以定量历史洪水纳入计算样本，用于参数估计与设计洪水计算。
表2宜昌站几种方法参数估计结果与设计洪水对比表
相对于模拟序列而言，宜昌站资料相对较多，调查历史洪水资料相对丰富，均有效地增加的样本容量，改善了样本的代表性，提高参数估计与设计洪水的代表性和准确性。
鉴于实测资料难以知晓序列统计参数和不同重现期设计洪水的真值，在对比分析中，增加三峡水库设计洪水的相关成果。在这之后，随着水文资料积累增多，曹广昌与王俊对三峡设计洪水进行复核(曹广昌，王俊.长江三峡工程水文泥沙观测与研究[m].2015,244-245.)，增加了1998年和1999年大洪水资料，并延长至2012年，经复核表明，三峡工程设计洪水成果较为稳定，p-ⅲ型分布参数(
鉴于此，将三种计算方法所得参数估计、设计洪水计算成果与三峡水库设计成果进行对比。由表2可见，本发明的方法计算与三峡水库设计值之间的相对误差最小，其中待估参数(
此外，将本发明方法计算的设计洪水值与定量历史洪水相比，1870年长江洪水时历史上所调查到的毁灭性很大的洪水，经过14c测定，2500年以来，从未发现有大于1870年历史调查洪水更大的洪水(曹广昌，王俊.长江三峡工程水文泥沙观测与研究[m].2015,242.)，其洪峰流量为105000m3/s，本发明方法计算万年一遇设计洪水设计值为101520m3/s，与之较为接近，而其他两种方法相差较大，可见，本发明方法计算设计洪水成果较为合理。
最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求保护范围之内。
1.一种考虑历史洪水的洪水频率估计方法，其特征在于，包括以下步骤：
s1、根据流域特性或规范要求，选取合适的频率线型，确定待估参数；
s2、收集水文站监测的流量资料，构建实测洪峰样本序列，并计算确定待估参数的上、下边界，所述实测样本序列为由s个年最大洪峰流量数值构成的序列；
s3、收集历史洪水信息，构建历史洪水序列，包括定量和不定量历史洪水信息，所述定量历史洪水信息由k1个确定的洪峰流量数值构成的序列，不定量历史洪水信息由k2个不定量历史洪水的流量区间构成的序列；
s4、将历史洪水序列与实测洪峰序列纳入计算样本，构建考虑定量及不定量历史洪水信息的参数估计似然函数；
s5、将似然函数作为目标函数，待估参数边界作为约束条件，建立优化问题，并采用遗传算法进行求解，进而得到频率线型的待估参数；
s6、根据选定的线型及估计参数，计算不同重现期的设计洪水。
2.根据权利要求1所述的考虑历史洪水的洪水频率估计方法，其特征在于，所述步骤s1采用《水利水电工程设计洪水计算规范》(sl44-2006)推荐使用的皮尔逊ⅲ型作为频率线型，待估参数θ包括均值
式中，f(·)为皮尔逊ⅲ型分布概率密度函数，α为形状参数，β为尺度参数，a0为位置参数，γ(·)为伽马函数。
3.根据权利要求2所述的考虑历史洪水的洪水频率估计方法，其特征在于，所述步骤s2中采用矩法计算实测洪峰样本序列待估参数的抽样误差，以此来确定待估参数的上、下边界。
4.根据权利要求3所述的考虑历史洪水的洪水频率估计方法，其特征在于，所述步骤s4将定量及不定量历史洪水均纳入计算样本，增加样本数量，似然函数如下：
式中，ls(d|θ)为实测洪水估计参数似然函数，lh(d|θ)为历史洪水似然函数，二者乘积作为最终参数估计似然函数l(d|θ)，f(·)是皮尔逊ⅲ型概率密度函数；k为历史洪水个数，其中k1为可定量洪水个数，k2不定量历史洪水个数；x0为洪峰阈值；f(·)为皮尔逊ⅲ型分布累计密度函数；h(＝h-s)为调查洪水历史时期长度，h为洪水重现期。
5.根据权利要求4所述的考虑历史洪水的洪水频率估计方法，其特征在于，所述步骤s5建立优化问题如下：将步骤s4构建的似然函数作为目标函数，将步骤s2计算得到的待估参数的边界作为约束条件，使用遗传算法进行寻优求解，通过选择、交叉、变异操作，输出寻优过程最优解，进而得到频率线型的待估参数。
6.根据权利要求5所述的考虑历史洪水的洪水频率估计方法，其特征在于，所述遗传算法进行寻优求解的过程如下：根据待估参数的上、下边界，随机生成初始群体，计算群体中各个个体的适应度，选择群体中适应度高的个体，直接遗传至下一代，通过配对交叉产生新的个体再遗传至下一代，并对群体中的个体进行变异产生新一代，进而得到下一代群体，以进化过程中所得到的具有最大适应度个体作为最优解输出，终止计算，得到待估参数值。
7.根据权利要求2所述的考虑历史洪水的洪水频率估计方法，其特征在于，所述步骤s6计算公式如下：
xp＝β·g-1(1-p,α,1) a0
式中，xp为重现期t的设计洪水，频率
