CN111651427A - 一种基于gamlss模型的非一致性水文频率计算方法 - Google Patents

Abstract

本发明提供一种基于GAMLSS模型的非一致性水文频率计算方法，包括构建分布函数库；构建拟合模式库；运用GAMLSS模型在每种拟合模式下对非一致性水文序列进行拟合；比较不同模式下拟合结果的AIC值，选取AIC值最小的模式作为备用最优模式；采用似然比检验法，检验备用最优模式的拟合效果是否明显大于其他更简单的模式，确定最优模式；通过残差分析，评价最优模式的拟合效果；判断是否需要重新构建分布函数库和拟合模式库；根据最优模式及其拟合结果，分析得到水文序列的非一致性频率。本发明构建了能适用于绝大多数水文序列的分布函数库以及能精确描述分布参数与时间的定量关系的拟合模式库，提出的非一致性水文频率计算流程和方法科学合理、操作性强。

Description

一种基于GAMLSS模型的非一致性水文频率计算方法

技术领域

本发明涉及水文频率分析领域，具体的说是一种基于GAMLSS模型的非一致性水文频率计算方法。

背景技术

水文序列是一定时期内反映水文规律现象的资料，受气候条件、自然地理条件和人类活动等多个因素的影响。近年来，由于气候变化和人类活动导致的下垫面变化，造成了水文序列的物理成因发生了变化，表现为水文序列的概率分布型式或分布参数在研究时间范围内发生了显著变化，传统水文频率分析中的一致性假设不在成立。准确分析非一致性条件下水文序列的频率，对灾害预测、防洪抗旱决策、水利工程设计施工、资源开发利用以及用水安全保障等具有重要意义。

相比于还原还现法、分解合成法、水文模型法以及全概率公式法等方法，时变参数法可以定量分析水文序列分布参数随时间的线性或非线性变化特征，是分析非一致性水文频率的有效手段。其中，GAMLSS(Generalized Additive Model for Location，Scale andShape)是一种可以描述随机变量概率分布参数与解释变量间的线性或非线性关系的(半)参数回归模型，包含了适用于绝大部分随机变量的分布函数，可以拟合过离散、高偏度以及高峰度的随机变量，具有适应性强、模拟精度高等优点，被逐渐用于非一致性水文频率计算。目前，采用GAMLSS模型分析非一致性水文频率还没有一套标准的流程和成熟的方法，需要进一步研究。

发明内容

针对现有技术的不足，本发明提出了提供一种基于GAMLSS模型的非一致性水文频率计算方法，给出了最优模式不满足要求时的解决方案，降低了技术风险。

本发明一种基于GAMLSS模型的非一致性水文频率计算方法，包括以下具体步骤：

步骤1，选择八种常用水文序列拟合分布函数，构建分布函数库；

步骤2，针对每种分布函数，设定分布参数的拟合函数类型，组合形成48种拟合模式，构建拟合模式库；

步骤3，基于R软件编程，以时间为解释变量，运用GAMLSS模型在每种拟合模式下对非一致性水文序列进行拟合，通过模型参数估计，得到不同模式下的拟合结果；

步骤4，比较不同模式下拟合结果的AIC值，选取AIC值最小的模式作为备用最优模式；

步骤5，采用似然比检验法，检验备用最优模式的拟合效果是否明显大于其他更简单的模式，确定最优模式；

步骤6，通过残差分析，评价最优模式的拟合效果；

步骤7，判断最优模式的模拟效果是否满足要求，如果不满足要求，需重新选取其他的分布函数构建分布函数库和拟合模式库，重复上述步骤3-6；

步骤8，根据上述步骤选择的最优模式及其拟合结果，得到非一致性水文序列的拟合分布函数以及其分布参数与时间的函数关系，即水文序列的非一致性频率。

所述步骤1中构建了由Gamma、Gumbel、Logistic、Weibull四个两参数分布函数和Box-Cox Cole and Green、PeasonⅢ、Generalized Gamma、Generalized InverseGaussian四个三参数分布函数组成的分布函数库。

所述步骤2中拟合函数类型设定为，

两参数分布函数设定以下5种拟合函数类型：1)μ、σ都不随时间变化；2)μ随时间线性变化，σ不随时间变化；3)μ、σ都随时间线性变化；4)μ随时间非线性变化，σ不随时间变化；5)μ、σ都随时间非线性变化；

三参数分布设定以下7种拟合函数类型：1)μ、σ、ν都不随时间变化；2)μ随时间线性变化，σ、ν不随时间变化；3)μ、σ随时间线性变化，ν不随时间变化；4)μ、σ、ν都随时间线性变化；5)μ随时间非线性变化，σ、ν不随时间变化；6)μ、σ随时间非线性变化，ν不随时间变化；7)μ、σ、ν都随时间非线性变化。

所述步骤3中模型参数估计方法如下：

对于一个三参数分布函数，首先构造关于模型参数的似然函数：

式中，f(y_i|μ,σ,ν)表示水文序列第i个观测值y_i(i＝1,2,…,N)服从的概率密度函数，然后对两边取对数得：

最后以lnL(μ,σ,ν)取值最大为目标，估计模型参数的最优值。

所述步骤4中采用AIC值法选取备用最优模式，其方法为，

AIC＝GD+2df (3)

式中，

为全局拟合偏差，

为最优模型参数所对应的对数似然函数；df为模型的整体自由度，在数值上等于模型参数的个数；选取AIC值最小的模式作为备用最优模式。

所述步骤5中采用似然比检验法确定最优模式，其方法为，

首先设定零假设(H₀)和备择假设(H₁)，H₀：简单模式拟合效果更好，H₁：备用最优模式拟合效果更好；

然后构造检验统计量

其中

和

分别为简单模式和备用最优模式在各自已拟合的模型参数下的似然函数；

最后计算p_value值，如果p_value小于设定的显著性水平，则放弃零假设，认为备用最优模式的拟合效果更好，并将备用最优模式作为最优模式；否则，则接受零假设，认为简单模式的拟合效果更好，并选择简单模式作为最优模式。

所述步骤6中采用残差分析评价最优模式的拟合效果，操作方法为，

首先计算原序列各观测值所对应的顺序统计值，得到经验残差r_i'；然后将经验残差和模拟得到的理论残差(r_i')组成的数据对(r_i,r_i')点绘在平面直角坐标系中；最后绘出1:1的辅助线作为参考直线；数据点与参考直线的偏差越小，表明理论残差和经验残差越接近，最优模式的拟合效果越好。

本发明的技术效果：(1)本发明所构建的分布函数库包括了Gamma分布在内的两参数分布函数和PeasonⅢ分布在内的三参数分布函数，能适用于绝大多数水文序列。

(2)本发明所构建的拟合模式库既考虑了分布函数不随时间变化，也考虑分布函数随时间的线性或非线性变化，能更加精确描述分布参数与时间的定量关系。

(3)本发明在对比分析AIC值的基础上，采用似然比法选择最优模式，避免选择的最优模式过度复杂。

(4)本发明采用残差分析进一步评价最优模式拟合效果，并给出了最优模式不满足要求时的解决方案，降低了技术风险。

(5)本发明提出的基于GAMLSS模型的非一致性水文频率计算流程和方法科学合理、操作性强。

附图说明

图1为基于GAMLSS模型的非一致性水文频率计算方法流程图；

图2为构建的分布函数库；

图3为构建的拟合模式库；

图4为残差分析示意图；

图5为入流量序列最优模式拟合残差分析图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1至图5所示，一种基于GAMLSS模型的非一致性水文频率计算方法，包括以下步骤：

步骤1，选择八种常用水文序列拟合分布函数，构建分布函数库；

步骤2，针对每种分布函数，设定分布参数的拟合函数类型，组合形成48种拟合模式，构建拟合模式库；

步骤3，基于R软件编程，以时间为解释变量，运用GAMLSS模型在每种拟合模式下对非一致性水文序列进行拟合，通过模型参数估计，得到不同模式下的拟合结果；

步骤4，比较不同模式下拟合结果的AIC值，选取AIC值最小的模式作为备用最优模式；

步骤5，采用似然比检验法，检验备用最优模式的拟合效果是否明显大于其他更简单的模式，确定最优模式；

步骤6，通过残差分析，评价最优模式的拟合效果；

步骤7，判断最优模式的模拟效果是否满足要求，如果不满足要求，需重新选取其他的分布函数构建分布函数库和拟合模式库，重复上述步骤3-6；

步骤8，根据上述步骤选择的最优模式及其拟合结果，得到非一致性水文序列的拟合分布函数以及其分布参数与时间的函数关系，即水文序列的非一致性频率。

具体实施例：

运用本发明计算湖北省漳河水库1963-2016年各季度入库径流量序列的非一致性频率。

步骤1，构建图2所示的分布函数库。

步骤2，构建图3所示的拟合模式库。

步骤3，拟合水文序列

对漳河水库冬季、春季、夏季、秋季入流量序列，基于R软件编程，运用GAMLSS模型在每种拟合模式下进行拟合，得到不同模式下的拟合结果。

步骤4，选择备用最优模式

由步骤三得到了每种模式下各入流量序列的拟合结果(见表1)，对冬季、春季、夏季、秋季入流量序列，分别选择AIC值最小的作为备用最优模式，即模式47、模式26、模式38以及模式36。

表1入流量序列拟合结果(AIC值)

注：-表示水文序列在该模式下拟合出错；下划线表示水文序列在该模式下拟合的AIC值最小。

步骤5，确定最优拟合模式

对每个入流量序列，采用似然比检验法检验备用最优模式的拟合效果是否明显优于其他更简单的模式，计算得到的p_value(见表2)，并取显著性水平为0.05，进行分析。

对冬季入流量、春季入流量和夏季入流量序列，所有的p_value均小于0.05，因此采用备用最优模式作为这三个序列的最优拟合模式。对秋季入流量序列，备用最优模式与模式21、模式23和模式24相比的p_value为0.083、0.083和0.099，都大于0.05，说明备用最优模式没有明显优于这三种简单模式。此外，在这三种模式中，模式24的AIC值最小，因此采用模式24作为秋季入流量序列的最优拟合模式。因此，对于冬季、春季、夏季和秋季入流量序列的最优拟合模式分别为模式47、模式26、模式38以及模式24，对应的分布函数及参数组合类型分别为：(1)μ、σ、ν都随时间非线性变化的GG分布；(2)μ随时间线性变化，σ、ν不随时间变化的P-Ⅲ分布；(3)μ随时间非线性变化，σ、ν不随时间变化的P-Ⅲ分布；(4)μ、σ、ν都不随时间变化的GIG分布。

表2入流量序列备用最优模式检验结果(p_value)

注：-表示水文序列在该模式下拟合出错；**表示该模式为备用最优模式；*表示该模式比备用最优模式更复杂。

步骤6，评价最优拟合模式

采用残差分析方法对每个入流量序列的最优拟合模式的拟合效果进行评价，得到的结果如图5所示。数据点与1:1辅助线的偏差越大，表明最优模式对序列的拟合效果越差。如果偏差过大，表明最优模式对序列的拟合效果不能满足要求，须放弃该最优模式。由图5可看出，各入流量序列的最优模式都具有很好的较好的拟合效果，数据点基本上在辅助线附近。

步骤7，计算入流量概率分布

根据最优模式的拟合结果，得到了冬季、春季、夏季、秋季入流量序列的概率分布参数与时间的函数关系，如表3所示。

表3入流量序列最优拟合模式的分布参数

注：t代表的是公元纪年；入流量的单位是万m³。

根据表3，则可计算每个入流量序列在过去或未来某年的概率分布。例如，计算2017年漳河水库冬季入流量的概率分布时，可令t＝2017代入表3第一行的相关公式，可得到分布参数为μ＝4.37×10⁷,σ＝0.0704,ν＝-3.1014，并代入分布函数库中GG的概率密度函数公式，即可得到2017年漳河水库冬季入流量的概率密分布。

此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (7)

1.一种基于GAMLSS模型的非一致性水文频率计算方法，其特征在于，包括以下具体步骤：

步骤1，选择八种常用水文序列拟合分布函数，构建分布函数库；

步骤2，针对每种分布函数，设定分布参数的拟合函数类型，组合形成48种拟合模式，构建拟合模式库；

步骤3，基于R软件编程，以时间为解释变量，运用GAMLSS模型在每种拟合模式下对非一致性水文序列进行拟合，通过模型参数估计，得到不同模式下的拟合结果；

步骤4，比较不同模式下拟合结果的AIC值，选取AIC值最小的模式作为备用最优模式；

步骤5，采用似然比检验法，检验备用最优模式的拟合效果是否明显大于其他更简单的模式，确定最优模式；

步骤6，通过残差分析，评价最优模式的拟合效果；

步骤7，判断最优模式的模拟效果是否满足要求，如果不满足要求，需重新选取其他的分布函数构建分布函数库和拟合模式库，重复上述步骤3-6；

步骤8，根据上述步骤选择的最优模式及其拟合结果，得到非一致性水文序列的拟合分布函数以及其分布参数与时间的函数关系，即水文序列的非一致性频率。

2.根据权利要求1所述的一种基于GAMLSS模型的非一致性水文频率计算方法，其特征在于：所述步骤1中构建了由Gamma、Gumbel、Logistic、Weibull四个两参数分布函数和Box-Cox Cole and Green、PeasonⅢ、Generalized Gamma、Generalized Inverse Gaussian四个三参数分布函数组成的分布函数库。

3.根据权利要求1所述的一种基于GAMLSS模型的非一致性水文频率计算方法，其特征在于：所述步骤2中拟合函数类型设定为，

两参数分布函数设定以下5种拟合函数类型：1)μ、σ都不随时间变化；2)μ随时间线性变化，σ不随时间变化；3)μ、σ都随时间线性变化；4)μ随时间非线性变化，σ不随时间变化；5)μ、σ都随时间非线性变化；

三参数分布设定以下7种拟合函数类型：1)μ、σ、ν都不随时间变化；2)μ随时间线性变化，σ、ν不随时间变化；3)μ、σ随时间线性变化，ν不随时间变化；4)μ、σ、ν都随时间线性变化；5)μ随时间非线性变化，σ、ν不随时间变化；6)μ、σ随时间非线性变化，ν不随时间变化；7)μ、σ、ν都随时间非线性变化。

4.根据权利要求1所述的一种基于GAMLSS模型的非一致性水文频率计算方法，其特征在于：所述步骤3中模型参数估计方法如下：

对于一个三参数分布函数，首先构造关于模型参数的似然函数：

式中，f(y_i|μ,σ,ν)表示水文序列第i个观测值y_i(i＝1,2,…,N)服从的概率密度函数，然后对两边取对数得：

最后以lnL(μ,σ,ν)取值最大为目标，估计模型参数的最优值。

5.根据权利要求1所述的一种基于GAMLSS模型的非一致性水文频率计算方法，其特征在于：所述步骤4中采用AIC值法选取备用最优模式，其方法为，

AIC＝GD+2df (3)

式中，

为全局拟合偏差，

为最优模型参数所对应的对数似然函数；df为模型的整体自由度，在数值上等于模型参数的个数；选取AIC值最小的模式作为备用最优模式。

6.根据权利要求1所述的一种基于GAMLSS模型的非一致性水文频率计算方法，其特征在于：所述步骤5中采用似然比检验法确定最优模式，其方法为，

首先设定零假设(H₀)和备择假设(H₁)，H₀：简单模式拟合效果更好，H₁：备用最优模式拟合效果更好；

然后构造检验统计量

其中

和

分别为简单模式和备用最优模式在各自已拟合的模型参数下的似然函数；

最后计算p_value值，如果p_value小于设定的显著性水平，则放弃零假设，认为备用最优模式的拟合效果更好，并将备用最优模式作为最优模式；否则，则接受零假设，认为简单模式的拟合效果更好，并选择简单模式作为最优模式。

7.根据权利要求1所述的一种基于GAMLSS模型的非一致性水文频率计算方法，其特征在于：所述步骤6中采用残差分析评价最优模式的拟合效果，操作方法为，

首先计算原序列各观测值所对应的顺序统计值，得到经验残差r′_i；然后将经验残差和模拟得到的理论残差(r′_i)组成的数据对(r_i,r′_i)点绘在平面直角坐标系中；最后绘出1:1的辅助线作为参考直线；数据点与参考直线的偏差越小，表明理论残差和经验残差越接近，最优模式的拟合效果越好。

Images