CN110334406B - 一种考虑风速特大值的极值风速重现期确定方法和装置 - Google Patents

Abstract

本发明公开了一种考虑风速特大值的极值风速重现期确定方法和装置。本发明方法在历年年最大10min平均风速序列中，检验并挑选序列中的特大值；通过构建含有特大值情况下的期望与均方差计算方法，在经验频率计算的基础上，估算耿贝尔极值分布模型的位置参数、尺度参数，并迭代试算特大值出现的重现期，直到符合误差要求，得出最终的重现期及分布模型参数。本发明方法可有效解决采用频率分布模型分析历年年最大10min平均风速序列时，特大值重现期估计不准确导致的设计风速取值不可靠的问题。本发明方法提供了可靠的特大值重现期估算方法，并在此基础上获得的设计风速可确保工程抗风可靠。

Description

一种考虑风速特大值的极值风速重现期确定方法和装置

技术领域

本发明涉及结构抗风技术领域，具体涉及一种考虑风速特大值的极值风速重现期确定方法和装置。

背景技术

在工程抗风领域体系中，需计算设计风速确保抗风结构设计的可靠性。设计风速的取值一般采用工程所在区域的气象站历年年最大10min平均风速序列进行频率分析计算，《建筑结构荷载规范》GB50009-2012推荐耿贝尔(Gumbel)极值分布模型。

随着极端天气的出现，往往在历年年最大10min平均风速序列中会出现偏离序列点群中心的特大值。如果不对特大值重现期做特殊的处理或准确的估计，耿贝尔(Gumbel)极值分布模型拟合风速值与经验频率时无法全面顾及到所有的序列值，所得到的设计风速成果明显难以符合人们的常规认识，而特大值重现期估计的稍有偏差便会导致设计风速成果的较大或大范围的偏离。

在水利工程设计中，面对洪水特大值的出现，可以通过考证文献、历史洪水的调查进一步确定特大值的重现期范围；或者假定雨洪同频率，通过更长观测年限的暴雨资料，寻找到洪水特大值出现时的暴雨参数，由暴雨的重现期来替代洪水特大值的重现期。在设计风速计算时，无论是直接进行风灾调查还是间接的考证，均难以找到定量的标准来直接或间接的确定风速特大值的重现期。工程实践中面对特大值的出现时，人们常常在序列中将其剔除后再进行频率分布模型的参数估计，计算设计风速取值，人为的选择忽视它们的存在，虽然得到了理想的频率分布模型，但忽视了特大值出现对样本分布的影响，会使得设计风速偏小、结构抗风安全的可靠性不足；或者人们直接根据历史记录时长来判断特大值的重现期，这样处理虽然考虑了特大值的作用，但往往因为记录时长较短，使得特大值的重现期及频率脱离整体的样本，使得设计风速成果偏大、抗风结构投资大量浪费。

发明内容

本发明的目的在于克服目前采用含有特大值的历年最大10min平均风速序列进行频率分析计算设计风速时，特大值的重现期难以考证或直接剔除特大值忽视其重现期对样本频率分布的影响，导致频率分布模型采用年最大风速极值序列分析设计风速不准确的问题，提出一种考虑风速特大值的极值风速重现期确定方法和装置。

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

在一方面，本发明提供了一种考虑风速特大值的极值风速重现期确定方法，包含以下步骤：

步骤1，采集气象站历年实测离地10m高年最大10min平均风速构成年最大风速极值序列；

步骤2：将年最大风速极值序列中的风速从大到小排列；将年最大风速极值序列中序列值满足预设的特大值判定条件的风速作为特大值；将所有特大值从大到小排列形成特大值序列，按照样本容量假设得到特大值序列对应的重现期序列，并作为第一次特大值重现期序列；

步骤3，求取年最大风速极值序列每个风速值对应的经验频率值，在此基础上求解耿贝尔极值分布模型的位置参数α和尺度参数β；根据求得的位置参数α和尺度参数β确定特大值序列对应的重现期序列，将其作为第二次特大值重现期序列；

步骤4，将第一次特大值重现期序列与第二次特大值重现期序列内的对应序列值之差的绝对值求和，当求和结果大于设定阈值δ时，则将第二次特大值重现期序列代替第一次特大值重现期序列，重复步骤3～步骤4；直到上述求和结果小于等于设定阈值，则根据第二次特大值重现期序列确定极值风速重现期及其对应的位置参数α和尺度参数β。

在以上技术方案中，所述特大值判定条件如下式：

上式中，x_i为年最大风速极值序列的序列值，XT_i为特大值序列中的特大值，n为年最大风速极值序列的样本容量，τ_θ＝0.01为显著性水平为0.01时t分布的临界值。

进一步地，特大值从大到小排列形成的特大值序列

(XT₁，XT₂，……，XT_l)，假定特大值序列对应的重现期序列分别为(T₁₁，T₁₂，……，T_1l)，其中

T_1k＝n-k+2 k＝1,2,......,l 式(1-4)

上式中，k为特大值的排序顺序，l为特大值序列的样本容量，即特大值个数。

进一步地，年最大风速极值序列中不含特大值的其它每个风速值对应的经验频率值的计算公式如下：

其中l为极大值个数，n为年最大风速极值序列的样本容量，T_1l为特大值序列中最后一个序列值对应的重现值，i为年最大风速极值序列中每个风速值序号。

进一步地，步骤3中求解耿贝尔极值分布模型的位置参数α和尺度参数β时，构建新的序列y_i，则x、y的期望、均方差的计算方法如下：

y_i＝-ln[-ln(1-p_i)] 式(1-6)

在所得到的期望、均方差的基础上通过下式求解位置参数α和尺度参数β的估计值：

α＝S(y)/S(x) 式(1-11)

β＝E(x)-E(y)S(x)/S(y) 式(1-12)。

其中l为特大值的个数，p_i为年最大风速极值序列第i个风速值对应的经验频率值，x_i为年最大风速极值序列中序列值，XT_k为特大值序列中的序列值，n为年最大风速极值序列的样本容量。

进一步地，根据位置参数α和尺度参数β确定特大值序列对应的重现期序列即第二次特大值重现期序列，其表达式如下：：

其中XT_k为特大值序列中的序列值，l为特大值个数。

进一步地，所述方法还包括：根据确定的极值风速重现期确定耿贝尔极值分布模型的位置参数、尺度参数；利用耿贝尔极值分布模型计算工程设计人员所需重现期T的设计风速，采用下式：

其中α为位置参数、β为尺度参数。

在另一方面，本发明提供一种考虑风速特大值的极值风速重现期确定的装置，包括：

数据采集模块，用于采集气象站历史实测离地10m高年最大10min平均风速构成年最大风速极值序列；

特大值确定模块，用于将年最大风速极值序列中的风速从大到小排列；将年最大风速极值序列中风速的重现期满足预设的特大值判定条件的风速作为特大值；

第一次特大值重现期确定模块，用于将所有特大值从大到小排列形成特大值序列，按照样本容量假设得到特大值序列对应的重现期序列，并作为第一次特大值重现期序列；

第二次特大值重现期确定模块，用于求取年最大风速极值年最大风速极值序列每个风速值对应的经验频率值，在此基础上求解耿贝尔极值分布模型的位置参数α和尺度参数β；根据求得的位置参数α和尺度参数β确定特大值序列对应的重现期序列，将其作为第二次特大值重现期序列；

参数确定模块，用于将第一次特大值重现期序列与第二次特大值重现期序列内的序列值之差的绝对值求和，当求和结果大于设定阈值时，则将第二次特大值重现期序列代替第一次特大值重现期序列，重复步骤3～步骤4；直到上述求和结果小于等于设定阈值，则根据第二次特大值重现期序列输出极值风速重现期及其对应的位置参数α和尺度参数β。

进一步地，所述特大值确定模块包括特大值判定模块，其用于执行以下判定条件：

上式中，x_i为年最大风速极值序列的序列值，XT_i为特大值序列中的序列值，n为年最大风速极值序列的样本容量，τ_θ＝0.01为显著性水平为0.01时t分布的临界值。

与现有技术相比，本发明所达到的有益效果是：

本发明方法将特大值、非特大值形成单独的两个序列进行分析，但又综合考虑频率分布模型同时拟合两个序列的精度，提供了一种新的特大值处理方法，既通过样本整体的特征确定了特大值的重现期，又考虑了特大值对样本整体的影响。可有效解决含有特大值的历年最大10min平均风速序列进行频率分析计算设计风速时，风速特大值重现期估算不准确的问题，使本发明方法适用度更广更加灵活；

本发明对特大值处理后的拟合效果明显优于未对特大值处理的效果，无论是特大值序列、年最大值序列还是点群中心均做到了兼顾，是可靠的特大值重现期估算方法，并在此基础上获得的设计风速可确保工程抗风可靠。

附图说明

图1所示是本发明具体实施例历史最大10min平均风速变化曲线；

图2所示是本发明具体实施例特大值未处理时的耿贝尔极值模型拟合效果；

图3所示是本发明具体实施例特大值重现期处理后的耿贝尔极值模型拟合效果；

图4所示是本发明具体实施例方法流程图。

具体实施方式

下面结合附图、实例对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

本发明选取了气象站作为案例，该气象站始建于1959年，1960～2014年历年实测离地10m高最大10min平均风速如图1所示，离地10m高历年最大10min平均风速如表1所示。

表1某气象站离地10m高历年最大10min平均风速序列

(1)挑选特大值序列

对历年最大10min平均风速序列按照从大到小的顺序进行排列，具体如表2所示。

表2某气象站离地10m高历年最大10min平均风速排序后的序列

按照式(1-1)、(1-2)计算得到

σ(x)＝3.11，样本容量n＝55。根据t分布表显著性水平为0.01时，τ_θ＝0.01＝2.40，根据式(1-3)得到判别特大值的标准为

对序列值进行逐一判别，仅排序为1的37.06为特大值，则l＝1，XT₁＝37.06。按照式(1-4)，假设特大值XT₁的重现期T₁₁为56年，对应经验频率为0.0178。

(2)估算耿贝尔极值分布的参数及重现期试算

进行第一次试算，按照式(1-5)计算经验频率，见表3第3列；构建y_i，并按照式(1-6)计算其值，见表3第4列。

按照式(1-7)、(1-8)计算

按照式(1-9)、(1-10)计算

按照式(1-11)、(1-12)计算α＝S(y)/S(x)＝1.236/3.13＝0.395，β＝E(x)-E(y)S(x)/S(y)＝24.08-0.648×3.13/1.236＝22.44。

按照式(1-13)计算各个序列值对应的重现期，见表3第5列。

表3第一次计算成果汇总表

由上表可知，特大值XT₁的计算重现期T₂₁＝322.7935年，按照判断标准δ＝10l＝10×1＝10，

因此重复上述计算过程进行第二次试算，并T₁₁＝T₂₁＝322.7935，对应特大值XT₁的经验频率为1/322.7935＝0.0031。

按照式(1-7)、(1-8)计算

按照式(1-9)、(1-10)计算

按照式(1-11)、(1-12)计算α＝S(y)/S(x)＝1.199/2.687＝0.446，β＝E(x)-E(y)S(x)/S(y)＝23.884-0.573×2.687/1.199＝22.60。

按照式(1-13)计算各个序列值对应的重现期，见表4第5列。

表4第二次计算成果汇总表

由上表可知，特大值XT₁的计算重现期T₂₁＝634.5193年，按照判断标准

因此重复上述计算过程进行第三次试算，并T₁₁＝T₂₁＝634.5193，对应特大值XT₁的经验频率为1/634.5193＝0.00158。

按照式(1-7)、(1-8)计算

按照式(1-9)、(1-10)计算

按照式(1-11)、(1-12)计算α＝S(y)/S(x)＝1.188/2.637＝0.451，β＝E(x)-E(y)S(x)/S(y)＝23.864-0.563×2.637/1.188＝22.62。

按照式(1-13)计算各个序列值对应的重现期，见表5第5列。

表5第三次计算成果汇总表

由上表可知，特大值XT₁的计算重现期T₂₁＝671.7177年，按照判断标准

因此重复上述计算过程进行第四次试算，并T₁₁＝T₂₁＝671.7177，对应特大值XT₁的经验频率为1/671.7177＝0.00149。

按照式(1-7)、(1-8)计算

按照式(1-9)、(1-10)计算

按照式(1-11)、(1-12)计算α＝S(y)/S(x)＝1.187/2.634＝0.4508，β＝E(x)-E(y)S(x)/S(y)＝23.863-0.562×2.634/1.187＝22.616。

按照式(1-13)计算各个序列值对应的重现期，见表6第5列。

表6第四次计算成果汇总表

由上表可知，特大值XT₁的计算重现期T₂₁＝673.5062年，按照判断标准

满足条件，考虑到阈值误差，拾位数取整，则特大值XT₁的重现期为T₂₁＝670年，对应的耿贝尔极值分布模型的位置参数α、尺度参数β取值分别为0.4508、22.616。

通过经验频率、序列值点绘后，用求解的耿贝尔分布模型拟合对特大值处理前后的拟合效果对比分析(如图2、图3)，可以发现通过本发明对特大值处理后的拟合效果明显优于未对特大值处理的效果，无论是特大值、大值还是点群中心均做到了兼顾。

(3)设计风速的计算

假设重现期T如下表第(1)列(根据规程规范规定，不同级别工程所对应的设计风速重现期)，采用上一步骤最终的位置参数α、尺度参数β取值分别为0.4508、22.616，利用耿贝尔极值分布模型按照式(1-14)计算各个重现期的设计风速，如下表第(2)列所示。

表7不同重现期的设计风速计算成果

列(1)	列(2)
		重现期(年)	设计风速取值(m/s)
1000	37.94
		500	36.40
300	35.26
		200	34.36
100	32.82
		50	31.27
30	30.12
		20	29.20
10	27.61
		5	25.94
2	23.43

本发明提供了可靠的特大值重现期估算方法，并在此基础上获得的设计风速可确保工程抗风可靠。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (10)

1.一种考虑风速特大值的极值风速重现期确定方法，其特征在于，包含以下步骤：

步骤1，采集气象站历年实测离地10m高年最大10min平均风速构成年最大风速极值序列；

步骤2：将年最大风速极值序列中的风速从大到小排列；将年最大风速极值序列中序列值满足预设的特大值判定条件的风速作为特大值；将所有特大值从大到小排列形成特大值序列，按照样本容量假设得到特大值序列对应的重现期序列，并作为第一次特大值重现期序列；

步骤3，求取年最大风速极值序列每个风速值对应的经验频率值，在此基础上求解耿贝尔极值分布模型的位置参数α和尺度参数β；根据求得的位置参数α和尺度参数β确定特大值序列对应的重现期序列，将其作为第二次特大值重现期序列；

步骤4，将第一次特大值重现期序列与第二次特大值重现期序列内的对应序列值之差的绝对值求和，当求和结果大于设定阈值δ时，则将第二次特大值重现期序列代替第一次特大值重现期序列，重复步骤3～步骤4；直到上述求和结果小于等于设定阈值δ，则根据第二次特大值重现期序列确定极值风速重现期及其对应的位置参数α和尺度参数β。

2.根据权利要求 1所述的一种考虑风速特大值的极值风速重现期确定方法，其特征在于，所述特大值判定条件如下式：

上式中，x_i为年最大风速极值序列的序列值，XT_i为特大值序列中的特大值，n为年最大风速极值序列的样本容量，τ_θ＝0.01为显著性水平为0.01时t分布的临界值。

3.根据权利要求 1所述的一种考虑风速特大值的极值风速重现期确定方法，其特征在于，特大值从大到小排列形成的特大值序列(XT₁，XT₂，……，XT_l)，假定特大值序列对应的重现期序列分别为(T₁₁，T₁₂，……，T_1l)，即第一次特大值重现期序列，其中

T_1k＝n-k+2 k＝1,2,......,l 式(1-4)

上式中，k为特大值的排序顺序，l为特大值序列的样本容量，即特大值个数。

4.根据权利要求 1所述的一种考虑风速特大值的极值风速重现期确定方法，其特征在于，年最大风速极值序列中不含特大值的其它每个风速值对应的经验频率值的计算公式如下：

其中n为年最大风速极值序列的样本容量，T_1l为特大值序列中最后一个序列值对应的重现值，i为年最大风速极值序列中每个风速值序号。

5.根据权利要求 1所述的一种考虑风速特大值的极值风速重现期确定方法，其特征在于，步骤3中求解耿贝尔极值分布模型的位置参数α和尺度参数β时，包括：构建新的序列y_i，则x、y的期望、均方差的计算方法如下：

y_i＝-ln[-ln(1-p_i)] 式(1-6)

在所得到的期望、均方差的基础上通过下式求解位置参数α和尺度参数β的估计值：

α＝S(y)/S(x) 式(1-11)

β＝E(x)-E(y)S(x)/S(y) 式(1-12)

其中l为特大值的个数，p_i为年最大风速极值序列第i个风速值对应的经验频率值，x_i为年最大风速极值序列中序列值，XT_k为特大值序列中的序列值，n为年最大风速极值序列的样本容量。

6.根据权利要求 1所述的一种考虑风速特大值的极值风速重现期确定方法，其特征在于，根据位置参数α和尺度参数β确定特大值序列对应的重现期序列即第二次特大值重现期序列，其表达式如下：

其中l为特大值的个数，XT_k为特大值序列中的序列值。

7.根据权利要求 1所述的一种考虑风速特大值的极值风速重现期确定方法，其特征在于，步骤4中所述阈值δ的表达式如下：δ＝10l，其中l为特大值的个数。

8.根据权利要求 1所述的一种考虑风速特大值的极值风速重现期确定方法，其特征在于，所述方法还包括：根据确定的极值风速重现期确定耿贝尔极值分布模型的位置参数、尺度参数；利用耿贝尔极值分布模型计算工程设计人员所需重现期T的设计风速，采用下式：

其中α为位置参数、β为尺度参数。

9.一种考虑风速特大值的极值风速重现期确定的装置，其特征在于，包括：

数据采集模块，用于采集气象站历史实测离地10m高年最大10min平均风速构成年最大风速极值序列；

特大值确定模块，用于将年最大风速极值序列中的风速从大到小排列；将年最大风速极值序列中风速的重现期满足预设的特大值判定条件的风速作为特大值；

第一次特大值重现期确定模块，用于将所有特大值从大到小排列形成特大值序列，按照样本容量假设得到特大值序列对应的重现期序列，并作为第一次特大值重现期序列；

第二次特大值重现期确定模块，用于求取年最大风速极值序列每个风速值对应的经验频率值，在此基础上求解耿贝尔极值分布模型的位置参数α和尺度参数β；根据求得的位置参数α和尺度参数β确定特大值序列对应的重现期序列，将其作为第二次特大值重现期序列；

参数确定模块，用于将第一次特大值重现期序列与第二次特大值重现期序列内的序列值之差的绝对值求和，当求和结果大于设定阈值δ时，则将第二次特大值重现期序列代替第一次特大值重现期序列，重复以下步骤:求取年最大风速极值序列每个风速值对应的经验频率值，在此基础上求解耿贝尔极值分布模型的位置参数α和尺度参数β；根据求得的位置参数α和尺度参数β确定特大值序列对应的重现期序列，将其作为第二次特大值重现期序列；将第一次特大值重现期序列与第二次特大值重现期序列内的对应序列值之差的绝对值求和，当求和结果大于设定阈值δ时，则将第二次特大值重现期序列代替第一次特大值重现期序列；直到上述求和结果小于等于设定阈值δ，则根据第二次特大值重现期序列输出极值风速重现期及其对应的位置参数α和尺度参数β。

10.根据权利要求9所述的一种考虑风速特大值的极值风速重现期确定装置，其特征在于，包括：所述特大值确定模块包括特大值判定模块，其用于执行以下判定条件：

上式中，x_i为年最大风速极值年最大风速极值序列的序列值，XT_i为特大值序列中的序列值，n为年最大风速极值年最大风速极值序列的样本容量，τ_θ＝0.01为显著性水平为0.01时t分布的临界值。

Images