基于广义帕累托分布模型的短期风速资料抗风设计风速计算研究

王晓惠，潘晓春，张 洋

( 中国能源建设集团江苏省电力设计院有限公司，江苏 南京 211102)

0 引言

设计风速是涉及结构抗风安全的重要参数。实际生产工作中往往因为气象观测站点观测资料年限不够等问题，导致工程设计风速取值困难，为此许多学者及工程设计人员进行了有益的研究。

许多学者[1-6]将临时测站与气象站的同步观测资料进行相关性分析，并以此相关关系在气象站的统计风速基础上进行修正，获得工程区域的设计风速。由于许多情况下，同步资料的相关关系并不能达到修正的预期等原因，不少学者试着仅采用短期实测资料进行分析研究。董安正[7]等针对短期风速样本不足的缺陷，采用由Monte Carlo 模拟方法得到的“理想”数据样本对人工神经网络进行训练，然后由“非理想”的短期风速原始资料估计设计风速。然而搭建人工神经网络模型相对复杂，且不具有普适性。随着数值计算的发展，庞加斌[8]等、罗颖[9]等采用虚拟气象站法，根据周围气象站基本风速与海拔高度之间的关系来推算虚拟气象站的设计风速，但在复杂地形区域数值模型难以得到可靠的风速与海拔的关系。

基于极值理论的数理模型是具有长期观测资料的设计风速计算的基本方法。许多学者尝试在短期资料中改进极值理论以求获得期望的结果。张秀芝[10]等采用二项—对数正态复合极值分布拟合短期风速资料计算海面设计风速基准值，并以沿岸站长年大风经验公式计算订正值，基准值与订正值叠加作为设计风速。丘海珊[11]等利用短期风速及经验频率的对数关系进行线性拟合，计算风电场的设计风速。广义帕累托分布模型在极值模型的基础上，通过阈值设置增加样本容量，充分利用短期观测资料，模拟效果优于极值模型，越来越受到青睐[12]。国外学者[13-16]用GPD 模型，采用小时、天、周及月最大风速求解重现期为小时、天、周、月的设计风速。罗乃东[17]利用月最大风速及年最大风速构建新的样本，用GPD 模型尝试计算短期风速资料的设计风速。GPD 模型中阈值的选择直接影响设计风速取值，而常用的平均剩余函数图法[18-19]、Hill 图法[20]在阈值选择时存在较大的主观性，使得设计风速的取值存在不稳定的现象[21]。在已有的研究中，难见普适的、可靠的短期资料计算设计风速的数理统计方法。

本文基于短期风速观测资料，拟采用GPD模型推导重现期单位为年的设计风速计算方法，并对现有的阈值方法进行改进。采用多座气象站实测资料进行短期资料的年设计风速计算成果与长期资料的Gumbel 分布模型分析成果对比，验证方法的可靠性与适用性。

1 短期风速资料的GPD模型

1.1 GPD函数

利用长期测风资料计算设计风速时常采用年最大风速序列进行分析。年最大风速序列一般符合广义极值分布(generalized extreme value，GEV)。当形状参数为0 时，GEV 便退化为Gumbel(亦称极值I 型)分布，常被用于工程结构的抗风设计。

假设年最大风速序列{Xi}，阈值为u0，则超阈值系列{X｜X＞u0}的频率分布模型近似符合GPD，见式(1)。

式中：u0为阈值；σ为尺度参数，且恒为正值。

由于风速的随机性，可以认为超越阈值的风速发生也是随机的。因此超过阈值的大风年发生次数m为一随机变量，且服从Poisson 分布[22]，见式(2)。

式中：λ为超过阈值年平均发生次数，λ＞0；k为任意整数。

按超阈值模型的经典假设，每年超过阈值的大风发生次数服从以λ为参数的Poisson 分布，则在R年中共发生λR个超过阈值的大风。对应R年一遇大风的重现期见式(3)。

式中：Pe为GPD 模型设计风速对应的频率。

对应R年一遇设计风速VR可由式(4)计算。

GPD 的参数估计方法有多种，如矩估计、最小二乘估计、基于分位数估计、近似广义最小二乘估计等。考虑到采用的是尺度参数为0时GPD 模型，本文采用矩估计，得到尺度参数的计算公式如式(5)。

1.2 短期风速序列推算设计风速

Grigoriu 认为，对于短期数据，可采用月最大风速来推算年最大风速分布[23]。假定各月的月最大风速满足独立同分布的条件，此时风速的年最大分布和月最大分布满足式(6)[24]。

式中：Fa(x)和Fm(x)分别为风速年最大和月最大分布。

定义λ为超过阈值u0的月平均发生次数，FG(Y)为短期风速资料的GPD 分布，则有式(7)。

按照式(6)，以短期风速资料求解年设计风速VR，可得到式(8)。

式中，σ为短期风速资料GPD 模型的尺度参数。

经方程式变换得到式(9)：

考虑到R＞10 时，1-(1-1/R)1/12≈1/(12R)，则得到不同年重现期的设计风速近似计算公式，如式(10)。

1.3 阈值选取的改进

在GPD 模型中阈值u0的选取是准确估计分布参数的前提，并且对模型的有效性至关重要。常用的阈值选择方法有Hill 图法、平均超额函数(Mean Excess Function，MEF)图法[25]。

Hill 图法、MEF 图法在分别寻找稳定区域、近似线性区域的判断均属定性判断，具有主观性，阈值选择的误差较大，而Moriarty 提出的经验阈值法容易造成阈值选择过大，剔除较多的有效样本。本文在Moriarty 提出的方法基础上进行改进，推荐阈值按“样本均值加Δ倍的标准差”取值，Δ值的选取兼顾样本长度足够长及同时通过“K-S拟合优度检验”“以λ为参数的Poisson 分布检验”的最大值，并以江苏省的主要气象站为例进行验证分析。

2 阈值选取方法对比分析

2.1 资料的选择

为避免观测环境变化对风速观测成果的影响，确保短期资料具有较好的代表性，本文以江苏某海岛气象站点2009 年7 月1 日—2011 年6 月30 日整两年的日最大10 min 平均风速观测资料为基础，如图1 所示。考虑到连续的风速样本间可能存在较强的相关性，同时结合江苏地区大风过程不超过4 d 的特征，将日最大10 min 平均风速数据4 d 作为一个周期进行取样。选取每个周期的最大值，且相邻的最大值的间距不能小于半个周期。

图1 某海岛气象站日最大10 min平均风速变化图

2.2 现有阈值选取方法的可靠性分析

按照Hill 图法的绘制方式，将选取的日最大10 min 平均风速按降序进行排列，计算极值指 数Hk，点 绘 集 合{(k，Hk-1)；1 ≤k≤n-1}得 到Hill 图， 如 图2 所 示。 由 于Loretan、Phillips 验证了Hill 图法选取阈值的范围不会超过样本容量的10%[29]，因此Hk-1趋于稳定区间为[32，45]，相对应的阈值的选取范围为[14.8，15.3]，阈值取该范围的起始值14.8。按照MEF图法的绘制原理，得到图形如图3 所示，按照平均超额函数e(x) =E(X-x|X＞x)是线性的且斜率为正的特点，在[16.8， 17.4]范围内平均超额函数近似为斜率为正的直线，阈值即为16.8。另外可发现MEF 图法的成果与Hill 图推荐的范围[14.8， 15.3]并无交集。按照Moriarty 提出的经验阈值法，该海岛气象站日最大10 min 平均风速序列的均值、标准差分别为9.38、3.11，按照经验阈值法阈值取值为13.7。

图2 日最大10 min平均风速的Hill图

图3 日最大10 min平均风速的MEF图

根据现有不同的阈值选取方法得到对应的阈值及λ，本文估算了GPD 模型的尺度参数，发现尺度参数的估算值、样本容量随阈值减小而增大，如表1 所示。三种方法中，Hill 图法、Moriarty 经验法的λ值均大于1，且超越阈值的样本数大于月最大10 min 平均风速序列的样本数，较月最大风速样本的容量有所扩大；MEF图法的λ值小于1，且样本容量仅为11 个，不足月最大风速样本的50%，这一特点与使用GPD 模型增大样本容量的初衷不一致。

表1 不同阈值取值方法的GPD模型参数

本文分别采用K-S拟合优度检验法、χ2检验法对GPD 模型的拟合程度、λ值是否符合Poisson 分布进行了检验。如表2 所示的检验结果可知，现有的三种阈值选取方法均通过了置信度95%的检验。

表2 不同阈值取值的GPD模型参数的K-S检验、χ2检验

为说明不同阈值选取方法得到的设计风速取值准确性，本文列举了短期、长期资料采用Gumbel 分布模型的分析计算成果。短期资料采用同期的月最大风速序列进行分析，长期资料采用1974—2014 年年最大风速序列计算，其中以长期资料的计算成果为标准进行各个方法的误差分析，如表3 所示。

表3 不同方法的设计风速成果对比m/s

现有阈值选取方法的GPD 模型的设计风速计算成果表明，三种方法计算成果均明显偏小，误差自大到小依次为Hill 图法、Moriarty 经验法、MEF 图法，计算误差均在9.3%以上，Moriarty 经验法、Hill 图法的计算误差达到13.5%～17.1%。同为基于Gumbel 分布模型，采用短期月最大风速序列采用得到的设计风速明显小于长期年最大风速序列的成果，且小了约20%。

2.3 本文推荐方法的阈值选取分析

本文推荐阈值按“样本均值加Δ倍的标准差”取值，Δ值的选取兼顾样本长度足够长及同时通过“K-S拟合优度检验”、“以λ为参数的Poisson 分布检验”的最大值。具体的选取方法既u0=X-+Δ·s，其中Δ∈[-1.4，1.4]，变化步长为0.1，由大至小试算至GPD 模型不能同时通过“K-S拟合优度检验”及“以λ为参数的Poisson 分布检验”的第一个拐点的前一个Δ值，具体的计算成果如表4 所示。尺度参数的估算成果随着Δ值减小而增大，对应的不同重现期风速也增大，以100 a 一遇设计风速为例，Δ从1.4降低至-1.4 时，尺度参数由1.856 增大至4.578，设计风速由29.0 m/s 增加至52.9 m/s。以置信度95%，对GPD 模型进行K-S拟合优度检验，没有呈现出明显的规律，在Δ值由大变小时，检验依次呈现出通过、不通过、通过、不通过的交替规律。概率以λ为参数的Poisson 分布检验也呈现出同样的规律，但两者的通过检验的Δ值区间不同，具体如图4 所示，按照前文所述的取值方法，Δ按0.6 取值。

图4 两者检验置信度95%的通过情况及Δ值的寻找

表4 本文推荐的不同阈值取值对应的设计风速成果m/s

续表

通过阈值选取、设计风速误差对比分析，可见基于本文推荐阈值选取方法的GPD 模型得到的设计风速成果准确度更高，且通过了K-S拟合优度检验，λ值亦符合Poisson 分布。

通过如图5 所示P-P 图(概率图)、如图6所示Q-Q 图(分位数图)，可见基于本文推荐阈值选取方法的GPD 模型对样本的拟合效果较好，能可靠地反应样本的真实分布。

图5 本文推荐方法的P-P图

图6 本文推荐方法的Q-Q图

3 GPD模型应用分析

选取江苏10 座主要国家基本气象站，分别进行基于短期资料、长期资料的设计风速计算，并进行对比分析，用以验证GPD 模型进行短期资料设计风速计算成果的可靠性，GPD 模型的阈值选取采用本文推荐的方法。短期资料采用各气象站2009 年7 月1 日—2011 年6 月30 日两年的日最大10 min 平均风速资料。考虑到短期资料采用自计仪器获得，为确保可比性，长期资料亦选取自记资料，并使用Gumbel 分布模型计算设计风速。

10 座气象站均按照推荐的阈值选取方法找到了满足条件的阈值，具体如表5 所示。满足两种假设检验的情况下，以长期资料的设计风速计算成果为依据，分析短期资料的设计风速取值误差，如表6 所示。误差在±5%以内的气象站有5 座，不同重现期设计风速的取值误差在-3.7%～3.0%，其中4 座气象站的误差在±3%以内；误差介于±5%～±10%之间的气象站5 座，误差在-6.8%～6.7%。

表5 江苏13座气象站长、短期资料的设计风速成果m/s

表6 短期资料设计最大风速成果的误差分析

江苏省气象站实测年最大风速呈现逐年下降的趋势，主要因周边的探测环境受城市化影响而发生了变化。从而导致短期资料脱离了长期资料的样本特征，使得短期资料的分析成果误差较大。因此，本文推荐的GPD 模型在采用此类气象站点的实测资料计算设计风速时准确度不佳，而周围相对空旷、周边探测环境变化较小的气象站点的短期实测资料在采用本文推荐的方法时，计算误差基本满足实际工程的使用要求，计算成果基本可靠。

4 结论

通过短期资料进行设计风速计算的GPD 模型推导及阈值选取方法的比较，采用江苏10 座气象站的长期、短期风速资料的设计风速成果对比，得出如下的结论。

1)在使用GPD 模型计算短期资料的设计风速时，采用合理的阈值，可得到精度可靠的设计风速取值。

2)只有2 年的短期资料时，本文推荐采用均值与标准差的组合来确定阈值的取值，即u0=X-+Δ·s，其中Δ∈[-1.4，1.4]，变化步长为0.1，由大至小试算至GPD 模型不能同时通过“K-S拟合优度检验”及“以λ为参数的Poisson 分布检验”的第一个拐点的前一个Δ值。推荐的方法较常规的Hill 图法、MEF 图法，Moriarty 经验法精度更高。

3)采用GPD 模型进行短期资料的设计风速计算时，周围相对空旷、周边探测环境变化较小的地区相比城市化影响较大的区域，成果更加可靠。

本文研究GPD 模型在短期资料设计风速计算时采用了2 年的日最大10 min 平均风速序列，得出满意的结果。然而，对1 年、3 年等不同长度的短期资料时本文的阈值选取方法及其可靠性尚需进一步研究。