基于POT法确定风压系数极值的自动阈值选取与参数估计

陈元坤, 毛 丹, 李寿科, 刘 敏, 陈晓强, 陈 俊, 孙洪鑫

(1. 中南建筑设计院有限公司,武汉 430071; 2. 武汉大学 土木建筑工程学院,武汉 430072; 3. 湖南科技大学 土木工程学院,湖南 湘潭 411201)

风压系数极值是确定建筑围护结构设计风荷载的重要变量。区组最大极值方法(block maximum,BM)是风压系数极值估计常用方法,其理论基于三种渐进型极值分布(极值Ⅰ型、Ⅱ型、Ⅲ型),对多个区组的样本极值进行概率分布拟合估计极值,在数据样本较大的情况下被认为是精确方法。BM方法使用过程中忽略了区组中的第二大、第三大、等峰值数据,数据使用率较低,在数据样本量较少的情况下,BM方法难以给出准确的极值结果。因此,学者们提出了一些基于短时程样本的极值计算方法[1-4]。同时,为充分利用区组的峰值数据,Simiu等[5]提出可采用超越阈值模型(peak over threshold,POT)来解决数据样本较少的问题,用于极值风速估计。POT模型极值计算方法基于广义Pareto分布(generalized pareto distribution,GPD)理论。数据独立同分布是POT模型使用的前提,Simiu等提出相邻独立峰值间距必须要满足一定的时间间隔(4天或者8天),Duthinh等[6]使用均值穿越获得独立峰值数据,全涌等[7]则建议使用自相关分析方法确定独立峰值数据。而在概率分布参数估计方法方面, Haan[8]提出的GPD模型参数估计的无偏估计方法。Naess等[9]通过C0估计器改进Haan估计与矩估计方法。Hosking等[10]验证极大似然估计,矩估计,概率加权矩估计方法对GPD模型两参数估计的可靠性,用于水文领域。而对于风压系数极值估计领域,当前缺乏GPD参数估计适用性研究。阈值选取可决定POT极值估计结果的稳定性和准确性[11-12]。Davison等[13]提出用平均剩余寿命图衡量超出量函数随阈值变化的曲线,超出量函数趋于线性时,线性范围内阈值都为最佳阈值,该方法具有较大的主观性。Coles等[14-15]利用似然比检验证明形状参数在最佳阈值之后具有不变规律,此方法为图形诊断方法,具有一定的主观性,选择的阈值不唯一,存在人工观测过程。Duthinh基于W统计量方法研究自动选择阈值,并用Bootstrap方法验证阈值选取方法的准确性。Liang等[16]提出基于峰值稳定性的阈值自动法,实现自动选取最佳阈值,用于水文处理。李正农等[17-18]基于POT方法研究了风压系数极值求取问题,但缺乏对其阈值自动选取方面的研究。

本文主要针对POT在高层建筑风压系数极值计算中的参数估计和阈值选取方面展开研究,研究不同GPD分布参数估计方法的适用性,确定合适的参数估计方法用于风压系数峰值的GPD分布拟合;将变点理论引入阈值选取,使用局部比较法确定最佳阈值,建立新的阈值自动选取方法,对比不同方法的性能,完善当前风压系数极值估计的POT方法。

1 POT极值计算方法

1.1 GPD概率分布

对于独立同分布的样本峰值数据x,设定样本阈值为μ,超过阈值的数据样本(超阈值样本),其概率分布将向GPD分布收敛[19]。

(1)

式中:μ为位置参数;ξ为形状参数;σ为尺度参数。GPD分布尾部分位数xp-GPD及重现期由式(2)和式(3)给出,即

(2)

(3)

式中:xp-GPD为GPD分布概率为p的尾部分位数;n为独立数据总数;Nu为超阈值数目总数;λ为单位时间独立同分布数据数目;R为重现期。

1.2 独立峰值样本获取

极值事件服从泊松点过程,峰值样本独立同分布。风洞试验采集的风压系数时程样本,采样频率较高,样本之间具有明显的相关性。采取一定的方法提取独立峰值样本,才可依据POT方法进行风压系数极值估计。常用的独立峰值提取方法为均值超越峰值提取方法,指定一个特定值—均值,将平均值以下的风压时程舍去,平均值以上的风压时程将与上交叉点、下交叉点作为起点与终点进行分组,单个组提取峰值,组成独立数据,该方法已被证明是一性能优越方法[20]。

1.3 参数估计

1.3.1 矩估计

矩估计方法可用于GPD分布参数估计,根据超出量的前两阶矩即可计算GPD参数

(4)

(5)

式中:E(x-μ)为超出量的平均数;s(x-μ)为超出量的标准差;μ为阈值,同时也是位置参数。

1.3.2 概率加权矩估计

GPD分布参数可用概率加权矩估计方法计算,式(6)～式(10)给出概率加权矩估计计算方法,位置参数μ为阈值

(6)

(7)

(8)

式中,Wr为总体矩的样本估计量,对于有限的样本大小,一致的矩估计(Wr)可以计算为

(9)

总体矩估计可以用Landwehr提出样本估计量Wr来替代

(10)

式中:xi,n为超阈值的数据样本;n为样本数量。

1.3.3 De Haan估计

Haan提出了GPD分布参数估计方法基于以下两个估计量,即

(11)

(12)

式中,k为数据高于阈值的观测次数,阈值代表第(k+1)个最高数据点。最高,第二高,…,第k高,(k+1)高数据分别由xn,xn-1,…,xn-(k+1),xn-k表示。GPD参数由式(13)、式(14)给出

σ=ρxn-kHk,n

(13)

(14)

这里当ξ≥0时,ρ=1,当ξ<0时,ρ=1-ξ。

1.3.4 极大似然估计

将GPD分布分布函数作为似然函数,得出对数似然函数,即

(15)

1.4 阈值选取

本文给出三种阈值选取方法,基于W统计量的阈值选取方法和基于形状参数或极值结果稳定性结合变点理论的阈值选取方法。

1.4.1 基于W统计量的阈值自动选取方法

基于W统计量的阈值选取方法是根据阈值所选择数据对GPD分布的拟合优度来确定阈值,依据GPD分布拟合最佳效果确定最佳阈值。也可通过将数据转换为一种统计量,检验统计量与其对应分布拟合优度,W统计量就属于后一种转换统计量,与W统计量所对应的分布为标准指数分布,在式(16)中给出W统计量的计算方法

(16)

式中:ui为i时刻的阈值;Wi为对应的W统计量;ξi为此时的形状参数;σi为尺寸参数;μi为位置参数;Yi超出量。以下给出使用W统计量阈值选取方法的具体步骤:

步骤1计算W统计量;

步骤2给出W统计量与标准极值分布图,以超阈值计算标准极值分布的值为横坐标,对应的W统计量为纵坐标,绘制W统计量与标准极值分布图,并做横纵坐标均为标准极值分布的直线;

步骤3计算图上点到标准极值分布直线的距离之和,最小距离和时对应的阈值为最佳阈值。

1.4.2 基于形状参数或极值分位数稳定性-变点阈值自动选取方法

本文结合变点理论以及形状参数/极值分位数的稳定性提出阈值自动选取新方法。由变点统计理论,指出“模型中的某个或某些量起突然变化的点”为变点,变点能反映事物某种质的变化[21]。具体使用流程如下:首要任务就是确定阈值范围,确定一个有高概率存在最佳阈值的阈值范围。确定候选阈值范围需要考虑以下三点。 ①总体阈值范围。总体阈值范围起点为所有独立峰值风压系数的最小值,终点为保证超阈值样本数据量大于5。②候选阈值间隔。在确定总体阈值范围后,基于等间距方法确立候选阈值。候选阈值间隔过大会忽略掉隐藏在大间隔中的合适阈值,阈值间隔过小则使计算时间大幅增加。在本文中,基于独立峰值样本的大小,将候选阈值间隔确定为峰值样本最小间距,且保证每次计算得到的形状参数有变化。③总体阈值范围分段。依据峰值风压系数分布范围,将总体阈值范围均分为4个区间,从这些小区间中,选择其中的区间作为候选阈值区间。用形状参数差值比来衡量候选阈值区间是否存在最优阈值。在候选阈值区间内,形状参数差值比值保持为正值或负值,则表示其形状参数在递增或递减,并未趋于稳定,存在最优阈值的概率较小;如形状参数差值比围绕0值波动,则表示形状参数在波动,此候选区间存在合适阈值。

在确定候选阈值范围后,在候选阈值区间计算对应的形状参数或极值分位数。进而计算相邻3个阈值点对应的局部比较量的方差变化值,在方差变化最大值所对应的位置,即为变点位置,其对应的阈值,即为候选最优阈值。计算重现期为4.5组的POT极值分位数,基于概率相等原则与Cook &Maney的78%极值分位数对应进行比较。

2 极值分析数据来源

本文主要研究POT极值估计方法在风压系数极值估计中的应用,风压系数数据来源于CAARC高层建筑刚性模型测压风洞试验。试验在湖南科技大学风工程试验中心进行,模拟GB 50009—2012《建筑结构荷载规范》[22]规定的D类地貌。模型缩尺比1∶400,试验时风洞中参考高度与模型顶部高度一致,为45.72 cm,参考高度处风速为10.2 m/s。CAARC高层建筑足尺尺寸长×宽×高为45.72 m×30.48 m×182.88 m,模型缩尺比为1∶400,试验风向角定义以及模型测点布置如图1所示,总共布置有308个测压点。选取0°典型风向角,独立重复采样表面风压时程200次,单次采样时间为30 s,采样频率为330 Hz,每个测点单次采样数为10 000个。风速比为1/3.3,可得时间缩尺比为1/120,可得实验室采样约为足尺时长60 min。

图1 CAARC建筑模型示意图(cm)Fig.1 CAARC building model (cm)

测点风压系数时程Cp,ij

(17)

3 风压系数极值POT模型估计

3.1 独立峰值提取分析

3.1.1 提取独立峰值量

选择0°风向时峰度最大测点29,将其单次风压系数时程,使用均值超越法提取独立峰值。图2给出测点29原始风压系数时程,以及使用独立峰值提取方法后的风压时程图。

图2 风压系数时程图Fig.2 Wind pressure coefficient time chart

图2 中的原始数据个数为10 000个,使用独立峰值提取方法后独立峰值个数分别为2 251个。后续进一步验证提取峰值后样本的相关性和同分布性。

3.1.2 独立峰值自相关性分析

图3中给出Tap 29的原始数据,以及使用独立峰值提取方法提取的独立峰值自相关函数图。图3中可以明显看到Tap 29原始数据的自相关性较强,在时间间隔为大于3 s时,自相关系数降至0.2以下;使用均值超越法提取的独立峰值在非常小的时间间隔内,自相关系函数幅值下降至0.05以下,在-0.05～0.05波动,方法提取得到独立峰值相关性低,独立性强。

图3 自相关系数图Fig.3 Autocorrelation coefficient

3.1.3 独立峰值齐次点泊松过程检验

适合拟合为GPD分布的独立峰值,其相邻时间间隔t需要满足齐次点泊松过程。推测相邻独立峰值时间间隔是否满足齐次点泊松过程,通过对点泊松分布式(18)求反函数,将独立峰值时间间隔的经验分位数代入,图4中给出独立峰值的齐次点泊松过程检验。图4中独立峰值由圆点表示,超越均值法所提取的峰值与拟合直线的拟合优度良好。

图4 齐次点泊松过程检验图Fig.4 Test of Poisson point process

在提取独立峰值数量,以及自相关系数检验,齐次点泊松过程检验3个方面,均值超越法在独立峰值提取方面表现良好。

P(t)=1-e-λt

(18)

3.2 GPD分布参数估计方法分析

3.2.1 参数估计方法性能研究及评价标准

GPD参数估计方法多为二参数估计方法,需关注参数估计方法对形状参数、尺寸参数的估计不确定性,同时衡量对样本量的依赖性,本文采取蒙特卡洛方法进行参数估计方法性能评估研究,具体步骤如下所示:

步骤1对于检验参数,设定一个取值范围,形状参数设定为[-1.5,0.5],尺寸参数设定为(0,0.6],样本量设定为[10,500],固定除检验参数以外的参数,具体参数设置如表1所示,参数的变化幅值为0.1;

表1 参数估计方法参数设置表Tab.1 Parameter estimation method parameter setting

步骤2根据设定参数值,使用蒙特卡洛模拟GPD分布,生成1 000组数据,单次数据长度为1 000,采用相应的参数估计方法(概率加权矩、矩估计、De-Haan估计、极大似然估计)计算对应的GPD参数;

步骤3将设定参数GPD分布的99.9%保证率极值分位数,与蒙特卡洛模拟的1 000组数据所得分位数平均值相减,得到偏差,用以衡量参数估计方法准确性;同时计算均方根误差(均方根误差=偏差2+方差),判断参数估计方法的稳定性。

3.2.2 参数估计方法对形状参数估计的不确定性

图5中给出了基于形状参数变化下的不同参数估计方法得到的极值分位数偏差、均方根误差图。图5(a)中,除概率加权矩估计方法与De Haan估计方法,其余方法得到的POT极值分位数偏差均随着形状参数增大偏差增大,其中极大似然估计对于形状参数的变化最为敏感,在形状参数为[-0.3～0.5]内偏差开始呈指数形式上升,最大偏差为-15.09;其次影响较大的是矩估计,最大偏差为-9.30;概率加权矩估计与De Haan估计对于形状参数的变化不敏感。图5(b)中显示在形状参数为正值时,分位数均方根误差随形状参数增大而增大,均方根误差概率加权矩估计>De Haan估计>极大似然估计>矩估计。参数估计方法准确性易受到形状参数的影响,当形状参数为正值时,概率加权矩估计以及De Haan估计的准确性较好,而矩估计和极大似然估计的准确性则较差;对于估计结果稳定性,矩估计的稳定性较好,其他参数估计方法的稳定性比较接近。因此,建议使用概率加权矩估计方法、De Haan估计方法计算GPD分布的形状参数。

图5 基于形状参数变化的分位数比较图Fig.5 Quantile comparison based on the change of shape parameters

3.2.3 参数估计方法对尺寸参数估计的不确定性

图6 给出了基于尺度参数变化下的不同参数估计方法得到的极值分位数偏差、均方根误差图。由图6(a)可以看出,极值分位数偏差总体随着尺寸参数的增大而增大,其中De Haan估计偏差绝对值最大,最大为-0.26,其次为极大似然估计,最大偏差为-0.12,概率加权矩估计与矩估计参数估计方法的偏差总体较小,接近于0值。由图6(b)可以看出,针对于估计方法结果的稳定性,极值分位数的均方根误差随形状参数增大而增大,De Haan估计受到影响最大,最大值为 0.11,其次概率加权矩估计>矩估计>极大似然估计。

图6 基于尺寸参数变化的分位数比较图Fig.6 Quantile comparison based on the change of size parameters

综上所述,De Haan参数估计方法受到尺寸参数的影响较大,随着尺寸参数的增大(正值范围内),偏差与均方根误差均增大,即准确性与稳定性降低;其次为极大似然估计方法,偏差较大,准确性较低;综合看来,概率加权矩估计方法对于GPD概率模型参数估计准确性与稳定性均较优。

3.2.4 参数估计方法对样本大小的依赖性

图7(a)中给出样本量变化时采用不同参数估计方法得到的GPD分布99.9%极值分位数偏差图。概率加权矩估计、矩估计、极大似然估计随样本量增大而分位数偏差减小,偏差范围为[0,0.5];De Haan估计的极值分位数偏差几乎不变,误差稳定在[-0.13,-0.10]。图7(b)给出样本量变化时不同参数估计方法得到的极值分位数均方根误差图,与偏差图一致,概率加权矩估计方法在样本量较小(<40)的情况下效果不理想,在样本量≥50的情况下,参数估计方法均方根误差十分接近,并且较小。然而,在极值风压系数采用POT超阈值方法进行估计,其独立峰值样本量均在100以上,因此对于样本量敏感性方面,本文推荐的这几个参数估计方法均可适用。

图7 基于样本量变化的分位数比较图Fig.7 Quantile comparison based on the change of sample size

综合本节不同参数估计方法对形状参数、尺寸参数、样本量的不确定性研究,可以看出,在进行GPD参数估计时,最佳参数估计方法宜选用概率加权矩估计方法。

3.3 阈值自动选取方法分析

3.3.1 阈值选取方法性能评估算例设置

本节将研究阈值选取方法对样本大小、样本非高斯特性的敏感程度。对于样本大小,将衡量单组,以及10,20,…,200组数据对阈值选取方法的依赖程度,比较POT估计极值与BM估计极值之间的偏差大小,从而衡量阈值选取方法的准确性及性能。对于样本非高斯特性方面,依据本文风洞试验结果,将风压系数时程的峰度与偏度划分成不同范围,研究不同非高斯程度样本基于POT方法和BM方法估计极值之间的误差,由此给出阈值选取方法对样本非高斯特性的依赖程度,具体算例参数设定如表2所示。

表2 阈值选取方法性能研究参数设置表Tab.2 Threshold selection method specific parameter setting

3.3.2 阈值选取方法对样本大小的依赖性

单组采样长度30 s,样本量在[1,200]组变化时,采用不同阈值选取方法选取阈值,依据POT极值估计方法得到极值分位数,与标准极值(基于200组的风压系数时程采用BM极值估计方法的78%保证率分位数)进行比较,图8中给出比较结果。

图8 样本量影响下不同阈值选取方法分位数比较图Fig.8 Quantiles comparison of different threshold selection methods under the influence of sample size

图8中给出典型测点Tap 29依据不同阈值选取方法进行POT极值估计结果,与BM方法得到的标准极值(黑色水平线)比较结果。由图8(a)可以看出,三种阈值选取方法在样本量为1～10组时,POT极值分位数均快速逼近标准极值,在样本量为10～100组时,三个方法所得结果非常接近,而样本量为100～200组时,基于W统计量的阈值选取方法反而会偏离标准极值,而基于变点自动阈值选取方法的POT极值依然非常接近标准极值。图8(b)给出所有测点在样本容量为1～10组时依据不同阈值选取方法得到的POT极值分位数与标准极值的误差率平均值,可以看出,三种阈值选取方法的POT极值误差率平均值均约5%,并随样本量的增加而减小,在样本量达到5组后趋于平稳,基于极值稳定性-变点阈值自动选取方法的POT极值误差平均值略高于其他阈值选取方法,基于形状参数稳定性-变点阈值选取方法的POT极值与基于W统计量阈值选取方法的POT极值总体误差接近。

综上所述,基于W统计量的阈值选取方法计算的POT极值会在样本量较大的情况下,有所偏差。基于极值稳定性-变点阈值自动选取方法计算的POT极值,在小样本容量下误差率略大于其他方法。基于形状参数稳定性-变点阈值自动选取方法计算的POT极值,在样本量变化的情况下有较高准确性,且稳定性较好。因此,从样本大小适用性可以看出,基于形状参数稳定性-变点阈值自动选取方法可较好的实现POT极值估计方法的阈值选择。

3.3.3 阈值选取方法对样本非高斯特征的适用性

为研究阈值选取方法对样本非高斯特征的敏感性,本节将选取0°风向所有测点的风压系数数据,样本长度约5个足尺1小时,依据其峰度偏度分布范围划分为5段,使用不同阈值选取方法,依据POT极值估计方法计算风压系数极值与标准极值(BM-GEV分布78%保证率极值)进行比较,表3～表5给出三种阈值选取方法计算的POT极值误差率。同时,将表中所示的偏度与峰度范围分为三种情况,低峰度-低偏度(偏度范围为[-0.1 , 0.3]且峰度范围为[2 , 3.5]),中等峰度-中等偏度(偏度范围为[0.3 , 0.8]且峰度范围为[3.5 , 5]),高峰度-高偏度(偏度>0.8或峰度>5)。

表3 基于形参稳定性的变点阈值选取方法极值误差率Tab.3 Based on shape parameter stability change-point threshold selection method extreme value error rate

表4 基于极值稳定性的变点阈值选取方法极值误差率Tab.4 Based on extreme value stability change-point threshold selection method extreme value error rate

表5 基于W统计量的阈值选取方法极值误差率Tab.5 Based on W statistic threshold selection method extreme value error rate

图9给出所有测点基于形状参数稳定性-变点阈值选取方法(简称“形参—变点阈值选取方法”)、基于极值稳定性-变点阈值自动选取方法(简称“极值—变点阈值选取方法”)、基于W统计量的阈值选取方法计算POT极值与标准极值的误差率平均值比较结果。

图9 不同非高斯特征的测点POT极值分位数误差Fig.9 Quantile error under the change of sample non-Gaussian characteristics

由图9可以看出,随偏度与峰度增大而POT极值误差率增大。在低偏度—低峰度、中偏度—中峰度区域的测点,形参—变点阈值选取方法与W统计量阈值选取方法计算的POT极值的误差率在3%～4%左右,低于极值—变点阈值选取方法的5%误差率。在高偏度—高峰度区域的测点,极值—变点阈值选取方法的误差率较低,为10.5%。综上所述,基于形状参数稳定性的变点-阈值选取方法适用于低峰度-低偏度、中峰度-中偏度的测点;基于极值稳定性的变点-阈值选取方法适用于高峰度-高偏度的强非高斯测点。

4 结 论

本文针对超越阈值模型极值估计中的关键步骤——参数估计方法和阈值选取方法开展了参数化研究,确定了不同方法的适用性。

(1) 均值超越独立峰值提取方法提取的独立峰值样本的自相关性低,泊松点过程拟合效果效果好,较好的满足了独立同分布。

(2) 通过蒙特卡洛方法衡量了形状参数、尺寸参数、样本量对参数估计方法的影响,依据GPD参数估计方法准确性、稳定性、方法效果方面的研究结果,建议使用概率加权矩估计方法估计广义Pareto分布参数。

(3) 基于变点理论提出了阈值选取方法,该阈值选取方法可实现阈值自动选取,受到样本量以及样本非高斯特性的影响较小,由此构建的改进超越阈值模型计算风压系数极值与标准极值的总体偏差小于5%。

(4) 对比研究了不同阈值选取方法的适用性,结果表明基于形状参数稳定性的变点-阈值选取方法适用于低峰度-低偏度、中峰度-中偏度的测点;基于极值稳定性的变点-阈值选取方法适用于高峰度-高偏度测点。