基于统计推断的GEV分布水文频率计算充分样本长度研究

戈立婷，宋松柏，2

(1.西北农林科技大学 水利与建筑工程学院，陕西 杨凌 712100；2.西北农林科技大学 旱区农业水土工程教育部重点实验室，陕西 杨凌 712100)

1 研究背景

水文频率分析的主要目的是根据频率模型估计水文变量的设计值及其相关重现期，计算过程具有不确定性，进而影响分析结果的可靠性。认知不确定性是由于对研究系统的认识不完全而产生的，取决于数据长度和模型推理重现潜在水文过程的能力。模型不确定性源于线型选择和参数估计方法，其不确定性分析方法主要有贝叶斯理论和NUSAP法和抽样分布理论。而水文计算和分析模型需要足够长的数据，充分样本长度是可靠建模和准确参数估计的先决条件。Su等评估了洪水灾害中采样造成的不确定性，结果表明增加数据长度提高了Gumbel模型的可靠性，降低了采样误差和不确定性。Hao等建立多指标干旱监测模型，发现至少需要30年的数据可以减少计算偏差。设计方案的运行性能会受到样本长度的影响。Soukissian等采用不同参数估计方法，在不同样本长度下进行模拟研究，研究结果表明广义极值(GEV)分布参数估计方法的性能取决于参数估计方法和选用样本容量。Khastagir等研究表明由于序列长度的变化，GEV分布参数会有很大差异。

目前，世界许多水文学工作者采用概率论的经验法则，结合区域水文计算经验，获得了水文频率计算30年的标准。数理统计中30的样本长度普遍基于正态分布假设，可以利用正态分布的一些性质来推断总体参数。根据这一方法，李鸿雁等采用数值模拟方法对随机生成样本和代表性站点的均值、均方差进行分析，采用假设性检验证明30年样本长度的充分性。然而样本长度的充分性与研究问题有关，标准很难界定。水文概率分布的复杂程度、待估计参数数目和统计推断方法的选取都与样本长度密切相关。这一经验法则，对于估计均值来说是充分合理的，而论证相对复杂且具有多个参数的分布是否充分值得商榷。其次，不同统计推理方法、不同参数的收敛速度差别较大，因此，在=30的情况下，计算结果差别也很大。此外，30年的样本长度并不能为不同重现期的设计值计算提供有效参考。洪水频率分析中，计算较大重现期的设计值，要得到可靠的结果，需要的样本长度也相应更长。Katz等的分析结果表明，GEV分布的极端分位数估计在小样本(<25)情况下非常不稳定。Robson等认为采用GEV分布进行频率计算需要的样本量是重现期的4倍。Vanem等认为35年数据足以估计100年的设计值。Hu等结果显示估计重现期为100年的设计值，35年的序列长度比20年的误差小50%。Merz等的研究表明一定样本长度下，增加更多的数据并不意味着可以提高计算精度。Benson等指出不超过序列长度两倍的重现期，能得到可靠的结果。而Jeong等认为数据长度外推3至4倍的重现期，其频率计算是合理的。

以上研究表明：(1)样本长度不充分，难以支撑分析结果的可靠性，模型的运行性能也会因样本短缺而受到影响。(2)30年样本长度这一标准并不能为不同频率下的设计值计算提供有效参考。(3)模拟和实例分析大都是基于固定长度的样本，样本长度分析方法较为主观。目前，不同频率下样本长度的充分性没有共识，水文频率分析中缺少一种基于严谨的数理统计推理而非经验法则的分析方法进行给定频率计算所需的充分样本长度推求，且在这一过程中，缺乏充分考虑到线型、参数估计方法、重现期的影响。基于此，本文应用渐近正态性质、极限近似理论原理，选用GEV分布，研究不同参数下充分样本长度的确定方法，并验证其合理性。在此基础上，根据误差指标进行分析比较，采用数值模拟方法分析参数、频率、参数估计方法与样本长度的潜在关系及规律，可在实际工作中用于判断计算数据是否足够充分用于模型的建立和参数估计。文中方法以期为确定水文频率所需样本长度，提高水文设计值成果的可靠性等提供支撑。

2 GEV分布及其参数估计方法

广义极值分布(Generalized Extreme Value Distribution，GEV)的概率密度函数()及其分布函数()分别为：

(1)

(2)

极值Ⅰ型分布是GEV分布的一个特例，形状参数=0，极值Ⅰ型分布的概率密度函数()及其分布函数()分别为：

(3)

(4)

式中：为尺度参数；为位置参数；为形状参数，当<0时，为极值Ⅱ型分布(或Frechet分布)，当>0时，为极值Ⅲ型分布(或Weibull分布)，=0时，为极值Ⅰ型分布(或Gumbel分布)。

最大似然法(Maximun Likelihood Method，MLM)是根据对数似然函数最大化进行概率分布参数估计的方法，要求对数似然函数是可微的。MLM与大多数参数估计方法相比，其优势在于它直接为GEV分布参数估计提供近似的标准误差和误差边界，同时允许将协变量信息合并到参数估计中。此外，MLM具有良好的渐近性质。由于这些优点，MLM是目前最主要的参数估计方法之一，也是GEV分布参数估计中最为常用的估计方法。Smith等和Coles等指出>-0.5时，MLM具有一般的渐近性质；-1<<-0.5，MLM结果存在，但不具有渐近性质；<-1时，MLM结果通常不存在；>0.5时，不存在二阶矩和高阶矩，此时极值分布有一个非常短的有界右尾，这在极值分布的应用中极少出现。因此，极大似然理论的局限性在实际应用中影响不大。

矩法(Method of Moments，MOM)的基本假设为样本矩等于相应的总体矩，对未知参数进行求解。详细的求解过程参见文献[25-27]。

3 充分样本长度确定方法

以GEV分布为例，数值模拟可以假设总体规模足够大，足以进行充分样本长度的分析。然而，在实际的工程设计中，选用的数据有限，必须论证选用资料长度估计结果的准确性和所需的充分样本长度。基于bootstrap的数学分析方法采用Shapiro-Wilk方法检验设计值的正态性，计算充分样本长度，判断计算数据长度是否足以进行参数估计。原假设定义为：样本长度为时，MLM设计值服从正态分布。如果拒绝原假设，说明在给定显著性水平下，样本长度不足，需要采集更多的样本，并再次应用上述方法进行检验，直到增加样本至无法拒绝原假设。详细步骤如下。

步骤一：给定参数(，，)，生成=600服从GEV分布的随机样本=(，，…，)。

步骤二：根据bootstrap法，从中抽取-1组相同容量的样本，得到再生样本=(1，2，…，)，其中=2，…，。

(5)

(6)

曲线拟合法是基于极限近似理论，用解析表达式来逼近离散数据。选取95%的相对置信区间(Relative confidence interval，RCI)，平均相对误差(the mean relative error，RMAE)两种不确定性分析指标。在水文频率的不确定性分析中，通常根据置信区间来量化设计值的不确定性。置信区间的宽度随样本长度的增加而减小，但随重现期的增加而增大。相比较而言，使不同重现期下无法直接对比的置信区间变为可比。避免误差相互抵消，可以反映实际误差的大小，也具有可比性。由某一样本长度下的多组设计值，计算得到、，随样本长度不断增大，根据其变化趋势，拟合数据点绘制出曲线，确定出最佳拟合方程。具体方法描述如下。

步骤一：给定参数(，，)，生成=2000组服从GEV分布，样本长度为(=3，4，…，500)的随机样本。

(7)

用对不同样本长度下设计值误差进行定量评价，定义如下式所示：

(8)

式中：为样本长度、频率下的估计设计值；为GEV分布下的真实设计值。

步骤四：随样本长度不断增大，根据各指标的变化趋势，拟合数据点绘制出曲线，由评价其拟合效果，确定出最佳拟合方程。

步骤五：将3.1中计算出的代入拟合方程，得到MLM估计下拟合曲线对应的、值。给定值，根据离均系数计算公式建立～|Ф|的曲线关系，与、值的变化趋势来验证文中方法的正确性。通过分析MLM、MOM估计下、与样本长度的拟合曲线的交叉点，进一步分析样本长度与参数估计方法间的关系。

4 数值模拟结果

本节基于GEV分布，应用上述方法进行数值模拟实验。根据最大似然估计量的渐近正态性，采用正态检验方法确定频率计算所需的充分样本长度，并验证其合理性。在此基础上，选用、为误差指标，应用曲线拟合法建立拟合关系表达式，定量分析不同参数、频率、参数估计方法与样本长度间的关系，根据充分样本长度的规律性进一步论证方法的正确性。最后，分析、两种参数估计方法在同一误差标准下的优劣。

根据GEV分布设计值的计算公式，设计值的大小与位置参数、尺度参数成正比，但随形状参数的变化呈e指数增加，不同频率下设计值对的变化更为敏感，值对拟合结果有较大影响。其次，GEV分布根据形状参数值的不同，分为极值Ⅰ型、极值Ⅱ型、极值Ⅲ型，三种分布的曲线形状也大不相同。因此，本文根据不同值和频率，计算所需的充分样本长度。在水文应用中，的取值范围一般为[-0.3，0.3]。本文给定GEV分布参数=(，，)，取100，取10，的取值范围为[-0.3，0.3]，步长为0.1。

选取频率为0.5%，不同值下的模拟结果进行分析。随样本长度变化Shapiro-Wilk检验的平均α值及其95%的上下置信区间如图1所示。从图1可以看出，整体呈现上升趋势，水平线α=0.05，对应垂线为充分样本长度：样本长度小于时，下置信区间位于水平线之下，或上下波动，样本长度大于时，下置信区间均位于水平线之上，说明在>时，最大似然估计量服从正态分布，估计结果是渐近有效的，证明样本长度充分，那么的样本长度就足以估计该频率下的设计值。

图1 频率为0.5%，不同形状参数下，Shapiro-Wilk检验的平均α值及其95%置信区间

样本长度为模拟100次，计算频率为0.5%的设计值，给出设计值的正态概率图，如图2所示。图2表明，不同条件下，散点均匀分布在红色参照线上，且位于95%的上下置信区间内，表明该样本长度下，计算得到的0.5%设计值分布近似服从正态分布。进一步对样本长度下0.5%的设计值，进行卡方拟合优度检验，显著性水平取=0.05，检验结果同样表明，设计值的分布近似服从正态分布，说明是指定频率下设计值计算所需的充分样本长度。

图2 频率0.5%，不同形状参数设计值的正态概率图

以相同的方法计算得到不同值和频率下，计算所需的充分样本长度，如表1所示，频率为5%～99%范围内30的样本长度基本能够满足不同参数下的设计值计算，而0.1%～2%范围内，即重现期为50～1000年时，30年的样本长度用于设计值计算是远远不够的，且随重现期增大所需样本长度也增加。>0，为极值Ⅲ型分布，曲线的头部趋平，因此有上限值，在0.1%～2%范围内所需的样本长度较短，但不适合水文数据拟合。

以表1中=-0.3、-0.1、0.1、0.3为例，计算不同样本长度的、，采用曲线拟合法得到-、-的拟合曲线，确定出最佳拟合方程=1，=2。结果表明幂率函数拟合良好，均大于0.9。将代入拟合方程得到估计下拟合曲线对应的值、值，如图3所示。图中的、表示在频率下用充分样本长度进行频率计算，设计值结果对应的误差值。两种误差指标的变化趋势大致相同，0.1%～75%范围内，随值增大，、值变小；75%～99%范围内，随值增大，、值变大。研究洪水问题时，<50%，随值增大，重现期减小，估计所需的样本长度减小。75%～99%为研究枯水问题时的常用频率范围，随频率增大，重现期增大，估计设计值所需的样本长度增加。图中存在明显的规律性：75%前后趋势有所不同，在0.1%～75%范围内随值增大，误差曲线依次降低，而在75%～99%范围内，随值增大，误差曲线依次升高，显然这与的取值有关。

表1 不同形状参数下频率计算所需的充分样本长度na

图3 na对应拟合曲线中的RCI、RMAE值

进一步分析值对设计值的影响，根据谢欣荣GEV分布离均系数的推导，当≠0时，频率对应的离均系数Ф的计算公式为：

=±{Γ(1+)-[-ln(1-)]}/[Γ(1+2)-Γ(1+)]

(9)

式中：>0时，取正；<0时，取负。

当=0时，频率对应的离均系数Ф计算公式为：

(10)

给定值，由离均系数计算公式可建立～|Ф|的曲线关系，如图4所示，不同值下，离均系数|Ф|的变化与图3中、值的变化趋势相同，符合客观规律，进一步说明了文中方法的正确性。

图4 不同形状参数下不同频率的离均系数|ФP|

采用参数估计方法计算同一参数和频率下不同样本长度设计值的、，根据其变化趋势，拟合数据点得到-、-的拟合曲线，结果表明幂率函数拟合良好，均大于0.9。由图3中对应的、值，可以得到同一误差标准下MOM频率计算所需的样本长度。给出频率为0.5%、99%，两种参数估计方法的对比结果，如图5所示。

由图5可见，同一、和同一参数估计方法的、拟合曲线的参数、值较为接近，曲线变化趋势相同。随值增大，与的曲线越来越接近，两种参数估计方法的误差随之减小。随样本长度的增加，、下降较为明显，之后趋于稳定。起初样本长度较小，比的估计误差小，得到的结果误差很大但相对稳定，而需要一定量的样本才能进行参数估计，在重现期较大且样本长度极小的情况下无法得到合理的结果。这与Martins等和Hosking等的研究结果相一致，在样本量较小的情况下，GEV分布的最大似然估计值不稳定，有时会得到不合理的形状参数估计值，在极值上分位数估计中造成较大的偏差和均方根误差。0.5%、99%代表重现期较大时的洪水、枯水稀遇事件，=-0.3时，同一误差指标下需要较大样本量，基本上无观测站能满足这样的条件。

图5 MLM、MOM估计下RCI、RMAE与样本长度的拟合结果对比

达到一定样本量，两种参数估计方法的、拟合曲线出现交叉点，如图6所示，<时，采用进行参数估计得到的设计值误差较小，>采用误差较小。0.1%～5%范围内随频率增大减小，在=10%～75%时趋于平稳，=75%～99%，逐渐增大，且随值增大，误差曲线依次变高，与图3、图4表现一致。

图6 两种参数估计方法的RCI、RMAE拟合曲线交叉点nc

由图3中对应的、值，可以得到同一误差标准下频率计算所需的样本长度。取-0.3、-0.1、0.1、0.3时，给出不同频率下两种参数估计方法的对比结果，如表2所示。在频率为5%～99%时，较小在30上下波动，值对样本长度的影响不大，仅在=-0.3时估计需要40左右的样本量，取30基本能满足计算要求，与李鸿雁等认为水文频率计算需要30年的样本结论相同。而在为0.1%～2%范围内普遍需要更长的样本用于频率计算。且不同参数和频率下计算得到的不同，这与熊立华等基于P-Ⅲ型分布的研究结论相一致，因为P-Ⅲ分布与GEV分布拟合一些水文序列较接近。而在为75%～99%范围内逐渐增大，且同频率下>0比<0大。由极值Ⅱ型与极值Ⅲ型的曲线形状可以发现，Ⅲ型曲线尾部相较于Ⅱ型曲线坡度更大，得到相对稳定的结果需要的样本长度更长，因此也更大。

表2 同一RCI、RMAE标准下MOM、MLM频率计算充分样本长度对比

5 实例分析

将上述方法应用于径流、降水、洪峰数据序列，以民和、兰州年径流序列，西安、周至年降水序列，宜昌、长洲、缆子湾、沙坪年最大洪峰序列为例，分析数据长度是否足以进行相应频率下的水文设计值计算，序列长度及其统计参数、GEV分布参数如表3所示。

表3 实测序列长度及其统计参数、GEV分布参数

基于GEV分布，针对年径流、年降水序列，选取具有代表性的特丰水年、丰水年、平水年、枯水年、特枯水年作为设计代表年，设计频率分别选择为5%、25%、50%、75%、95%，年最大洪峰序列选取重现期分别为200、100、50、30、10、5年进行充分样本长度的计算。以长洲站年洪峰序列为例，不同样本长度变化Shapiro-Wilk检验的平均α值及其95%的上下置信区间计算结果如图7所示。图中水平线α=0.05，对应垂线为该重现期下频率计算所需的充分样本长度，重现期越大，频率计算所需样本长度越大。

图7 长洲站年最大洪峰序列不同重现期Shapiro-Wilk检验的平均α值及其95%置信区间

将上述方法应用于民和、兰州年径流序列，西安、周至年降水序列，宜昌、缆子湾、沙坪年最大洪峰序列，得到法进行频率计算所需的样本长度，并作假设验证。采用假设检验中的ANOVA来论证该样本下计算洪水频率的合理性。ANOVA的基本思想是设定了一种描述数据变异的指标，分析出样本长度和随机误差影响和贡献，从而确定可控因素(样本长度)对研究结果影响力的大小。方差分析采用的检验统计量是统计量，即值检验，是组间与组内的离差平方和与自由度的比值。计算统计量的观测值，给定显著性水平，可做出决策，>0.05表示差异性不显著。

将相应频率所需的样本长度，与长度为+10的样本计算得到的各20组设计值进行方差分析，如表4和表5。表5中列“”表示选取序列长度不足以进行该重现期下的频率计算。其余均通过方差齐性检验，显著性水平均大于0.05，在基础上增加样本长度计算得到的设计值，总体上与计算得到的结果并无显著差异，说明该样本长度下进行水文频率计算是充分且合理的。

表4 不同设计频率计算的充分样本长度及其假设检验结果

表5 不同重现期洪水频率计算的充分样本长度及其假设检验结果

6 结论与讨论

本文首先阐述了极大似然估计量的渐近正态性质，界定了充分样本长度的概念，将bootstrap、Shapiro-Wilk正态性检验方法用于确定给定水文序列频率计算所需的充分样本长度，判断水文数据长度是否足够充分。基于GEV分布，根据形状参数的取值变化，在不同参数和频率下进行数值模拟分析，推求频率计算的充分样本长度。基于极限近似理论，选取、为误差指标，根据曲线拟合法分析参数、频率、参数估计方法与样本长度间的潜在关系规律。并应用于径流、降水、洪峰数据，推求不同设计频率(重现期)下频率计算的充分样本长度，进行假设检验进一步验证方法的可行性。经理论推导、模拟分析与实际应用，获得以下几点结论：

(1)文中充分样本长度的推理方法具有严谨的数理统计基础，不依赖于模型及经验法则，可避免主观因素的影响。

(2)文中模型不同频率下设计值的充分样本长度表明，重现期越大，所需样本长度越大。表明在样本长度不够充分的情形下，对结果进行不确定性评估是至关重要的。

(3)数值模拟及实例计算结果均通过检验，论证了文中方法的合理性。充分样本长度对应的误差指标与的取值有关，且在不同指标下呈现的趋势相同，也与相应值、相应频率下设计值离均系数的曲线趋势一致，符合客观规律，进一步说明了方法的正确性。

(4)最大似然估计法需要一定长度的样本量才能得到合理的结果，当矩法与最大似然法估计误差都较大时，矩法表现更好，交叉点之后误差显著减小，最大似然估计优于矩法，表现出其良好的估计性能。

(5)在10%～75%的频率范围内，充分样本长度及交叉点基本小于30，或在30上下波动，说明在这一频率范围内30的样本量足以进行最大似然估计，表现出其良好的估计性能。而在该频率范围外，交叉点的样本长度大于30，矩法、最大似然估计法在30的样本长度下存在较大误差。

本文研究成果可为水文频率计算中充分样本长度的确定提供一定参考依据，定量论证了水文数据进行分析计算的充分长度，以期提高水文分析计算的严谨性。但其中存在几点值得进一步深入研究：(1)本文方法基于数理统计推断，没有考虑物理测量、水文成因在内的其他因素影响。(2)文中假定水文序列是平稳序列进行了充分样本长度的分析计算，而对于非平稳序列频率计算：可采用混合概率模型、分段平稳模型，序列分解法模型，结合本文中充分样本长度计算方法进行研究。而对于时变参数模型，需要进一步研究时变参数的抽样评定方法。(3)文中基于GEV分布充分样本长度的分析方法可应用于其他概率模型，以探究该方法的普适性。