CN108955652A - 一种水文非一致性序列重构方法 - Google Patents

Abstract

本发明涉及水文频率分析技术领域，更具体地，涉及一种水文非一致性序列重构方法。包括以下步骤：S1.收集流域水文站点实测长序列流量或水位数据，采用年最大方法或超门限方法选取样本构成时间序列；S2.采用非一致性识别方法对样本的趋势性及突变性进行检验；S3.基于水文时序的非一致性特征，构建并筛选时间变化模型；S4.采用L矩法进行模型参数估计与设计值计算。本发明提供的一种水文非一致性序列重构方法，充分考虑了水文极值序列在环境变化影响下存在非一致性的特点，可提高设计结果可靠性与设计计算精度，可广泛应用于水文统计计算中。

Description

一种水文非一致性序列重构方法

技术领域

本发明涉及水文频率分析技术领域，更具体地，涉及一种水文非一致性序列重构方法。

背景技术

极端水文事件频率分析是基于实测样本数据提取极值序列，采用频率统计分析方法计算得出一定量级水文事件出现重现期。正确的重现期预估是地区防洪排涝政策制定及工程措施建设等的理论基础，对保障地区经济、生产安全具有重要作用。目前对于极端水文事件的设计值计算方法主要有四个步骤：一、提取极值样本序列；二、优选分布类型；三、估计模型参数；四、配线及设计值计算。总的来说，现有设计值计算方法没有考虑气候及人类活动造成样本统计特性改变，甚至不属于同一总体，即非一致性产生的问题。采用现有设计计算方法所得计算结果与变化的现实情况存在偏差，导致无法准确得出设计结果。

发明内容

本发明的目的在于克服现有技术的不足，提供一种水文非一致性序列重构方法，考虑了水文极值序列在环境变化影响下存在非一致性的特点，可提高设计结果可靠性与设计计算精度。

为解决上述问题，本发明提供的技术方案为：一种水文非一致性序列重构方法，考虑了水文序列的非一致性，包括以下步骤：

S1.收集流域水文站点实测长序列流量或水位数据，采用年最大方法或超门限方法选取样本构成时间序列；

S2.采用非一致性识别方法对样本的趋势性及突变性进行检验；

S3.基于水文时序的非一致性特征，构建并筛选时间变化模型；

S4.采用L矩法进行模型参数估计与设计值计算。

进一步地，所述的S1步骤包括：

S11.确定取样方法，若采用年最大法，则选取年最大值构成时序；若采用超门限方法则根S12中的方法确定超门限取样阀值；

S12.首先，判断样本的独立性，其判断方式为：θ＞5+ln(A)且X_min＜0.75min[Q₁,Q₂]，式中，θ为两个极值样本的间隔时间(单位：d)；A为流域面积(单位：km²)；Qi为第i次样本的对应量级；X_min为两个极值样本间的最小过程值；

然后，门限设定需满足以下三种要求：a)超门限样本的均值超过门限的部分是门限值S的线性函数；b)门限值的选择应使样本的分散指数在一个置信区间内，以确保超定量发生次数服从泊松Poisson分布；c)年均超定量发生次数μ>1.65，年均超定量发生次数μ在(1.65,3]的区间范围内，采用指数分布拟合超门限样本模型，μ过大易影响样本独立性；

最后，超门限样本的年发生次数序列服从泊松分布，同时超门限样本序列服从广义帕累托分布；

S13.对样本长度进行检验，样本个数少于30的对数据进行查补延长。

进一步地，所述的S3步骤包括：

S31.根据样本存在非一致性趋势特性构建时间变化模型，用简单的连续函数表述样本前二阶矩存在的趋势，考虑存在线性趋势和抛物线趋势的情况，共分析了六种时间-趋势种类：

1)均值具有线性趋势——(AL)；

2)标准差具有线性趋势——(BL)；

3)均值和标准差都具有线性趋势，并以一固定值作为比例(变差系数Cv)相关——(CL)；

4)均值和标准差都具有线性趋势，并且不相关——(DL)；

5)均值具有抛物线趋势——(AP)；

6)均值、标准差都有抛物线趋势，并以一固定值为比例(变差系数Cv)相关——(CP)。

还有一种基本情况是前两阶矩都不具有趋势，即分布参数不随时间变化，就是稳定状态情况(S)。

POT模型需要考虑的有两个序列，一个是用Poisson分布拟合的超定量年发生次数序列，一个是用GP分布拟合的超定量洪水序列，两个序列中任何一个序列存在趋势，都将影响AMS模型GEV分布的参数计算结果。由于Poisson分布参数λ为超定量年发生次数样本的数学期望，故超定量年发生次数序列只需考虑一阶矩的趋势；POT序列考虑前二阶矩的趋势。两序列结合考虑，可衍生出多种趋势模型，本发明共考虑15种，模型名称“ALS”表示超定量年发生次数样本均值具有线性趋势，超定量洪水序列前二阶矩均不具有趋势，“SAL”表示超定量年发生次数样本均值不具有趋势，超定量洪水序列均值具有线性趋势，其余模型同理。

S32.根据AIC准则计算得到最适合的时变模型，AIC值最小的趋势模型为最优模型，计算公式：

AIC＝-2lnML+2k

式中，ML为似然函数的最大值，为极大似然参数估计结果对应的似然函数值；k为模型参数个数；

基于最大熵原理AIC准则为最优趋势模型的选择标准。该准则考虑两部分内容，一是模型对样本的拟合效果，用似然函数值反应；二是模型稳定性，通过对模型的参数个数进行惩罚来实现。最终选取对数据拟合较好而参数数量尽可能少的模型为最优模型。增加模型参数可提高对样本的拟合效果，但可能因为过分强调对样本的拟合效果而降低曲线外延性，由于洪水频率分析更关心曲线的外延部分，所以在模型选择上，应尽量简化参数。

可检验出不同模型间差异显著性，且能综合权衡模型适用性和参数个数之间的关系，计算简单、客观。

进一步地，所述的S4步骤包括：

S41.采用L矩法根据样本直接得到样本L-矩(l_i，i＝1,2,3,4)、L-变差系数t₂、L-偏态系数t₃和L-峰度系数t₄；

研究样本序列的L-变差系数t₂，L-偏态系数t₃，L-峰度系数t₄分别定义如下：

其中，前四阶样本L-矩(l_i，i＝1,2,3,4)由bj(j＝0，1，2，3，4)计算得到：

S42.在L矩法计算样本矩及各类系数的基础上，分布参数的计算方法参照宋德敦提出的转化计算方法，根据L-矩和分布函数参数的关系估计参数bj(j＝0，1，2，3，4)，采用如下方式计算：

其中，n表示样本容量，x_k耀n表示样本序列中第k小的值，即x_1耀n≤x_2耀n≤…≤x_n耀n。

与现有技术相比，有益效果是：本发明提供的一种水文非一致性序列重构方法，充分考虑了水文极值序列在环境变化影响下存在非一致性的特点，可提高设计结果可靠性与设计计算精度，可广泛应用于水文统计计算中。

附图说明

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

图2为本发明实施例东江三站AMS洪峰流量序列CSDMC检验图。

图3为本发明实施例东江三站POT洪峰流量序列CSDMC检验图。

图4为本发明实施例龙川站基于时间基准点的洪水频率曲线变化情况。

图5为本发明实施例河源站基于时间基准点的洪水频率曲线变化情况。

图6为本发明实施例博罗站基于时间基准点的洪水频率曲线变化情况。

具体实施方式

附图仅用于示例性说明，不能理解为对本专利的限制；为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。附图中描述位置关系仅用于示例性说明，不能理解为对本专利的限制。

如图1所示，一种水文非一致性序列重构方法，包括以下步骤：

步骤1.收集流域水文站点实测长序列流量或水位数据，采用年最大方法或超门限方法选取样本构成时间序列；

S11.确定取样方法，如采用年最大方法，则选取年最大值构成时序，如采用超门限方法则采用S12确定超门限取样阈值；

S12.根据三条原则确定超门限取样阈值：

①样本独立；其判断方式如下所示：

θ＞5+ln(A)且X_min＜0.75min[Q₁,Q₂]

式中，θ为两个极值样本的间隔时间(单位：d)；A为流域面积(单位：km²)；Qi为第i次样本的对应量级；X_min为两个极值样本间的最小过程值。

②门限设定同时符合以下三种要求：a)超门限样本的均值超过门限的部分(X_S-S)是门限值S的线性函数；b)门限值的选择应使样本的分散指数在一个置信区间内，以确保超定量发生次数服从Poisson分布；c)年均超定量发生次数μ>1.65，此时采用指数分布拟合超门限样本模型能取得较好的效果，μ过大易影响样本独立性。

③超门限样本的年发生次数序列服从泊松(Poisson)分布，同时超门限样本序列服从广义帕累托分布(GP)。

S13.对样本长度进行检验，样本个数少于30的对数据进行查补延长。

步骤2.采用Mann-Kendall等非一致性识别方法对样本的趋势性及突变性进行检验；

步骤3.基于水文时序的非一致性特征，构建并筛选时间变化模型；

根据样本存在非一致性趋势特性构建时间变化模型，具体如下：

用简单的连续函数表述样本前二阶矩存在的趋势，考虑存在线性趋势和抛物线趋势的情况，共分析了六种时间-趋势种类：

1、均值具有线性趋势——(AL)；

2、标准差具有线性趋势——(BL)；

3、均值和标准差都具有线性趋势，并以一固定值作为比例(变差系数Cv)相关——(CL)；

4、均值和标准差都具有线性趋势，并且不相关——(DL)；

5、均值具有抛物线趋势——(AP)；

6、均值、标准差都有抛物线趋势，并以一固定值为比例(变差系数Cv)相关——(CP)。

还有一种基本情况是前两阶矩都不具有趋势，即分布参数不随时间变化，就是稳定状态情况(S)。

POT模型需要考虑的有两个序列，一个是用Poisson分布拟合的超定量年发生次数序列，一个是用GP分布拟合的超定量洪水序列，两个序列中任何一个序列存在趋势，都将影响AMS模型GEV分布的参数计算结果。由于Poisson分布参数λ为超定量年发生次数样本的数学期望，故超定量年发生次数序列只需考虑一阶矩的趋势；POT序列考虑前二阶矩的趋势。两序列结合考虑，可衍生出多种趋势模型，本发明共考虑15种，模型名称“ALS”表示超定量年发生次数样本均值具有线性趋势，超定量洪水序列前二阶矩均不具有趋势，“SAL”表示超定量年发生次数样本均值不具有趋势，超定量洪水序列均值具有线性趋势，其余模型同理。

S32、根据AIC准则计算得到最适合的时变模型，具体如下：

基于最大熵原理AIC准则为最优趋势模型的选择标准。该准则考虑两部分内容，一是模型对样本的拟合效果，用似然函数值反应；二是模型稳定性，通过对模型的参数个数进行惩罚来实现。最终选取对数据拟合较好而参数数量尽可能少的模型为最优模型。增加模型参数可提高对样本的拟合效果，但可能因为过分强调对样本的拟合效果而降低曲线外延性，由于洪水频率分析更关心曲线的外延部分，所以在模型选择上，应尽量简化参数。

可检验出不同模型间差异显著性，且能综合权衡模型适用性和参数个数之间的关系，计算简单、客观。计算公式：

AIC＝-2lnML+2k式中，ML为似然函数的最大值，为极大似然参数估计结果对应的似然函数值；k为模型参数个数。AIC值最小的趋势模型为最优模型。

步骤4.采用L矩法进行模型参数估计与设计值计算；

S41.采用L矩法根据样本直接得到样本L-矩(l_i，i＝1,2,3,4)、L-变差系数t₂、L-偏态系数t₃和L-峰度系数t₄(Hosking，1990)。具体如下：

研究样本序列的L-变差系数t₂，L-偏态系数t₃，L-峰度系数t₄分别定义如下：

其中，前四阶样本L-矩(l_i，i＝1,2,3,4)由bj(j＝0，1，2，3，4)计算得到：

S42.在L矩法计算样本矩及各类系数的基础上，分布参数的计算方法参照宋德敦(1989)提出的转化计算方法根据L-矩和分布函数参数的关系估计参数bj(j＝0，1，2，3，4)采用如下方式计算：

其中，n表示样本容量，x_k:n表示样本序列中第k小的值(x₁:n≤x₂:n≤…≤x_n:n)。

实时例1

以东江流域龙川、河源、博罗站点为例，采用本发明的计算方法，计算其对应10、20、30、50、100年重现期的洪峰流量设计值。

步骤1、收集流域水文站点实测长序列流量或水位数据，采用年最大方法或超门限方法选取样本构成时间序列。

本例采用东江流域龙川、河源、博罗站点的1954.1.1-2009.12.31逐日流量数据，从中提取年最大(AMS)与超门限样本(POT)。

步骤2、采用Mann-Kendall等非一致性识别方法对样本的趋势性及突变性进行检验。

本例采用Mann-Kendall和CSDMC两种方法识别样本存在趋势性非一致性。

M-K趋势检验结果见表1，结果表明东江三站年最大洪峰流量序列、超定量洪水年发生次数序列、超定量洪峰流量序列均存在下降趋势，除了博罗站的POT年发生次数序列和POT序列，其他所有样本的变化趋势均达到5％的显著性水平；龙川站三个样本和河源站AMS序列、POT年发生次数序列更是达到1％的显著性水平。

表1东江三站各洪水样本M-K趋势检验结果

AMS序列的CSDMC变化过程，如图2所示。龙川年最大洪峰流量序列的CSDMC曲线在1954-1970年上升，从1974年开始下降；河源站CSDMC曲线在1976年前上升，1976年后开始下降；博罗站1962年前CSDMC曲线斜率为正，1966-1987年较平稳，1987后开始下降；这些与枫树坝、新丰江、白盆珠水库的建造时间基本吻合；三个水文站AMS序列的CSDMC曲线均在下降至2004年后出现轻微回升。龙川站AMS序列可分为两个丰水段，即1954-1971年、2005-2009年；一个枯水段，即1972-2004。河源站洪峰流量序列可分为一个丰水段，1954-1976年；一个枯水段，1977-2004年；一个平水段，2004-2009年；博罗站洪峰流量可分为一个丰水段，1954-1962年；一个平水段，1962-1987年；一个枯水段，1988-2009年。

POT序列的CSDMC曲线，如图3所示，龙川站POT洪峰流量在1970年前大于均值，1970-2004年洪峰流量小于均值，2005-2009年洪峰流量有所增大。1966年以前，河源站POT洪峰流量大于均值；1966-1973年CSDMC曲线趋势不明显，POT洪峰流量约等于均值；1974-2009年曲线斜率总体上为负，洪峰流量小于均值，其中2004年-2009年波动较大。博罗站1962年前CSDMC曲线上升，1962-1987年间曲线趋势不明显，1987-2004年曲线总体存在下降趋势，2004年后有所回升。三个水文站POT洪峰序列的CSDMC曲线对比表明，龙川站POT洪峰流量的变化趋势比河源站明显，河源站波动较大，博罗站变化趋势在三站中最不明显，中间有一段时间POT洪峰流量较稳定。POT洪峰流量的变化过程及CSDMC分析结果与AMS序列较一致，龙川站洪峰流量变化较明显的时间点约为1970、2004年，河源站约为1974和2004年，博罗站约为1962、1987、2004年。受到水利工程调蓄影响下，东江洪水序列呈现出非一致性。

步骤3、基于水文时序的非一致性特征，构建并筛选时间变化模型。

结合AIC准则检验结果和各站超定量年发生次数序列均值和POT序列均值、标准差的变化过程，选择最优趋势模型。采用AIC准则检验GP-Poisson分布稳定模型及各趋势模型对超定量年发生次数序列、POT序列的拟合程度，如表2所示，选取AIC值最小的为各站最优趋势模型，龙川站拟合最优的为ALCP趋势模型，河源、博罗站为ALCL趋势模型。从趋势模型角度分析模型选择的合理性。根据各个站点最优趋势模型下的趋势模型参数，得到超定量年发生次数序列均值、POT序列均值和标准差的变化方程，绘制出各站超定量年发生次数序列均值和POT序列均值、标准差的变化过程。

龙川站POT年发生次数序列均值具有线性下降趋势，故该序列适合用AL模型拟合。POT序列均值具有“下降-平稳－上升”的抛物线趋势，而根据序列的波动情况，标准差也存在“下降-平稳”的抛物线趋势，且均值和标准差趋于平稳的时间较同步，故均值、标准差以C_v为比例相关具有抛物线趋势的CP模型能较好地描述龙川站POT序列变化趋势。两序列相结合表明，ALCP趋势模型作为龙川站最优趋势模型是合理的。

河源站POT年发生次数序列均值也具有线性下降趋势，采用AL模型描述。POT序列的均值和标准差均具有下降的线性趋势，采用CL模型描述POT序列变化趋势较合理。两序列趋势结合的ALCL模型与河源站洪峰变化情况较吻合。

博罗站POT年发生次数序列均值具有线性下降趋势，与龙川、河源对比趋势特征不显著，采用AL模型描述，斜率较小。POT序列的均值和标准差均具有下降的线性趋势，采用CL模型描述POT序列变化趋势较合理。ALCL模型反应的趋势特征能较好地描述博罗站洪峰变化情况。

表2龙川、河源、博罗站GP-Poisson各趋势模型AIC拟合检验值

步骤4、采用L矩法进行模型参数估计与设计值计算。

(1)龙川站

龙川站的3个时间基准点(1962、1989、2009)分别代表枫树坝水库建库前、建库后以及现状情况。洪峰流量未受枫树坝水库调蓄影响之前的阶段1954-1970，以1962年为代表时间基准点：与平稳条件下的AMS及GEV洪水频率曲线对比，以该年为时间基准点重构的AMS样本点据均在实测AMS样本点据上方，超过概率较高及较低部分的样本点，重构的值与实测值量级上差异较小，而中间部分重构值与实测值差异较大。因此，根据重构AMS序列绘制的GEV-ALCP曲线位于根据实测AMS绘制的GEV-S曲线的上方，GEV-ALCP曲线与GEV-S对比高水部分尾部较平缓，两曲线的中部位置差异最大。

枫树坝水库建成后，以1989和2009年为代表时间基准点，这两个时间基点对应的GEV-ALCP曲线位置均在GEV-S曲线下方，可见同量级洪水发生概率降低，说明龙川站年最大洪峰流量受水库调蓄削减洪峰作用影响明显。GEV-ALCP线型方面，1989年的频率曲线尾部与2009年比较陡峭，1989年的GEV-ALCP曲线低水端在2009年频率曲线上方，两曲线中高水端基本重合，在频率小到一定值后1989年频率曲线高于2009年。从频率曲线参数估计角度上看，1989年与2009年的频率曲线位置及中高水端较相似，与POT样本均值和标准差变化趋势在1975-2009年内较平稳有关，如图4所示，而频率曲线差异主要是年均超过次数序列的线性趋势导致的。

表3龙川、河源、博罗站最优趋势模型参数估计结果

(2)河源站

如图5所示，河源站的4个时间基准点为1956、1966、1989、2009年。在新丰江水库建库前，河源站流量不受水库调蓄影响，1956年时间基准点重构AMS序列及GEV-ALCL曲线位于各时间基点频率曲线的最上方，大量级洪水出现概率较大。新丰江水库蓄水、枫树坝水库建库前，受新丰江水库调蓄影响，1966年时间基点的GEV-ALCL曲线有所下降。枫树坝水库建成后，受两大水库蓄水影响，1989和2009的频率曲线均在GEV-S曲线下方，大洪水出现概率显著下降。

(3)博罗站

博罗站划分的时间区间较多，五个时间基准点分别为1956、1966、1978、1994、2009年。由于博罗站线性趋势的斜率小于河源站，各时间基准点的GEV-ALCL曲线线型较相似，随时间推移，GEV-ALCL曲线高水部分尾部逐渐下降，大量级洪水出现概率降低。1978年的GEV-ALCL曲线与GEV-S曲线基本重合，这与博罗站AMS序列的CSDMC曲线1966-1987年较平稳，如图6所示，即洪峰流量约为序列均值的非一致性诊断结果相吻合。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (4)

1.一种水文非一致性序列重构方法，其特征在于，考虑了水文序列的非一致性，包括以下步骤：

S1.收集流域水文站点实测长序列流量或水位数据，采用年最大方法或超门限方法选取样本构成时间序列；

S2.采用非一致性识别方法对样本的趋势性及突变性进行检验；

S3.基于水文时序的非一致性特征，构建并筛选时间变化模型；

S4.采用L矩法进行模型参数估计与设计值计算。

2.根据权利要求1所述的一种水文非一致性序列重构方法，其特征在于，所述的S1步骤包括：

S11.确定取样方法，若采用年最大法，则选取年最大值构成时序；若采用超门限方法则根S12中的方法确定超门限取样阀值；

S12.首先，判断样本的独立性，其判断方式为：θ＞5+ln(A)且X_min＜0.75min[Q₁,Q₂]，式中，θ为两个极值样本的间隔时间(单位：d)；A为流域面积(单位：km²)；Qi为第i次样本的对应量级；X_min为两个极值样本间的最小过程值；

然后，门限设定需满足以下三种要求：a)超门限样本的均值超过门限的部分是门限值S的线性函数；b)门限值的选择应使样本的分散指数在一个置信区间内，以确保超定量发生次数服从Poisson分布；c)年年均超定量发生次数μ在(1.65,3]的区间范围内，采用指数分布拟合超门限样本模型；

最后，超门限样本的年发生次数序列服从泊松分布，同时超门限样本序列服从广义帕累托分布；

S13.对样本长度进行检验，样本个数少于30的对数据进行查补延长。

3.根据权利要求1所述的一种水文非一致性序列重构方法，其特征在于，所述的S3步骤包括：

S31.根据样本存在非一致性趋势特性构建时间变化模型，

S32.根据AIC准则计算得到最适合的时变模型，AIC值最小的趋势模型为最优模型，计算公式：

AIC＝-2ln ML+2k

式中，ML为似然函数的最大值，为极大似然参数估计结果对应的似然函数值；k为模型参数个数。

4.根据权利要求1所述的一种水文非一致性序列重构方法，其特征在于，所述的S4步骤包括：

S41.采用L矩法根据样本直接得到样本L-矩(l_i，i＝1,2,3,4)、L-变差系数t₂、L-偏态系数t₃和L-峰度系数t₄；

研究样本序列的L-变差系数t₂，L-偏态系数t₃，L-峰度系数t₄分别定义如下：

其中，前四阶样本L-矩(l_i，i＝1,2,3,4)由bj(j＝0，1，2，3，4)计算得到：

S42.在L矩法计算样本矩及各类系数的基础上，分布参数的计算方法参照宋德敦提出的转化计算方法，根据L-矩和分布函数参数的关系估计参数bj(j＝0，1，2，3，4)，采用如下方式计算：

其中，n表示样本容量，x_k：n表示样本序列中第k小的值，即x_1：n≤x_2：n≤…≤x_n：n。