极值风速风向的联合概率密度函数

发布日期：2025-01-04 15:21    点击次数：125

在建筑设计中, 风荷载是计算风振响应、确定抗风设计的基础.对于高层建筑、大跨度结构等柔性结构, 风荷载是主要的控制性荷载, 合理的风荷载值关乎工程建设的安全性和经济性.作为典型的随机动力荷载, 风荷载通常是以极值风速为变量通过建立风速的概率密度函数而确定的, 如我国建筑结构荷载规范(GB 50009-2012)[1]规定, 以极值Ⅰ型分布作为极值风速的概率分布模型确定建筑设计风荷载.国内外学者开展了大量关于极值风速概率密度分布数学模型和计算理论的研究[2-4], 提出了有效的样本筛选、模型建立和参数优化方法.在样本筛选方面, 提出了年最大值法、跨阈法和独立风暴法等抽样方法；在数学模型方面, 建立了Gumbel分布、Frechet分布、Weibull分布、广义Pareto分布[5]等极值分布理论；基于数学模型和参数特征, 发展了矩法、极大似然法、粒子群算法等参数优化方法.尽管这些方法理论不一, 但在以极值风速为单值变量的基本风速预测中, 均具有较好的计算优度.然而, 风向角也是极值风速的重要特征参数, 仅考虑风速大小而忽略风向角是不合理的.Simiu等[6]研究了飓风区的极值风荷载, 指出若不考虑风向的影响, 50 a重现期的极值风荷载明显偏保守.Goyal等[7]研究了风向对混合住宅群的结构可靠性影响, 结果表明忽略风向将高估计算值.王钦华等[8]研究了风向对某超高层建筑等效静力风荷载的影响, 结果表明不考虑风向影响的建筑物等效风荷载偏保守, 应考虑风向对基本风压的影响.日本风荷载规范考虑风向对建筑物效应的影响, 给出了建筑结构设计的风向折减因子[9].因此, 以极值风速和风向角为双变量, 建立极值风速风向的联合概率密度(joint probabilistic density function, JPDF)模型, 对更精确、合理地计算结构风荷载具有重要的现实意义. 有关极值风速风向的JPDF模型研究较少, 主要原因是风向角变量是周期性的, 且观测记录多为非连续的方位角.Johnson等[10]基于谐波函数建立了离散角变量的连续概率密度结构, 能有效地拟合方向角变量的概率密度直方图, 但该模型采用多个三角函数拟合, 形式繁琐, 且无法给出固定的概率分布.目前, 比较常见的角变量分布有均匀分布、心形分布、包柯西分布、缠绕正态分布和von Mises分布等；其中, von Mises分布和混合von Mises分布被认为是最有效的描述角变量统计特性的分布[11], 在图像分析[12]、大气污染防治[13]、风能评估[14]等领域得到了广泛的应用. 在风速风向的JPDF建模方面, 陈隽等[15]采用谐波函数模拟风向角分布, 并基于不同风向间风速分布相互独立的假定建立了风速风向的JPDF分析方法；范文亮等[16]基于乘法定理导出了离散-连续混合联合分布模型, 并建立了风向风速的二维连续联合分布模型.这些模型需要对每个方位角的极值风速样本一一给出概率分布, 并对不同风向下的极值风速分布相关性做出假定.然而, 由于极值风速样本数量稀少, 分布在某些方位角下的极值风速样本量更少, 一般较难获得各方位角下的极值风速分布.事实上, 风速和风向的联合分布属于角度-线性分布.Johnson等[10]从理论上推导得出, 当给定角度变量和线性变量的边缘分布, 可以根据最大熵原理导出2个变量的JPDF结构.目前国外已有学者将该理论应用于风能预测[17]领域. 本文根据我国某地的月极值风速和对应风向记录, 分别建立极值风速的Gumbel分布和风向的二阶混合von Mises分布模型；在此基础上, 使用最大熵原理构建极值风速风向的联合概率密度结构, 并与Copula函数建立相互关联；采用该联合分布模型计算不同重现期下的基本风速, 与建筑结构荷载规范(GB 50009-2012) 中的规范值进行对比；最后, 分别采用Spearman秩相关系数和线性-角度变量相关系数对模型的相关性予以验证, 探究模型的有效性. 1 极值风速风向联合分布理论 1.1 极值风速分布 在不同的抽样方法下, 学者分别发展了相适应的数学模型和参数估计理论用于重现期下的基本风速计算[4].目前, 运用最多的2种抽样方法分别是年最大值法和跨阈法[5].一般认为, 以年最大值形成的风速样本服从极值Ⅰ型(即Gumbel)分布, 我国建筑结构荷载规范(GB 50009-2012) 即采用该方法计算基本风速.跨阈法通过设置特定的阈值, 建立由超越阈值风速形成的极值风速样本, 并假定该极值风速样本服从广义帕累托(generalized pareto distribution, GPD)分布.为准确计算极值风速的边缘分布, 采用上述2种模型进行对比. Gumbel分布和GPD分布的数学模型如下： ${F_v}\left( v \right) = {\rm{exp}}[ - {\rm{exp}}\left( { - y} \right)].$ (1) ${G_v}\left( v \right) = 1 - {\left( {1 + by} \right)^{ - 1/b}},y = \left( {v - u} \right)/a.$ (2) 式中：a、u和b分别为尺度参数、位置参数和形状参数, v为极值风速度量. 对Gumbel分布, 采用极大似然法确定a、u的参数估计；对GPD分布, 采用核拟合优度统计量法确定u, 再根据极大似然函数法计算a、u的最佳近似值. 1.2 风向圆周分布 气象站一般以有限方位角形式记录风向, 如国家基本站和一般站采用16个方位角整编风向.由于风向角的非连续性, 通常基于风向概率密度直方图对风向分布进行建模. 假设有一组离散风向角数据, 分散于单位圆(0≤θ＜2π)的不同角度区间.将单位圆等分为m份, 集中在第k(1≤k≤m)区间的观测点数为nk, 则各区间的风向角频度Pk和概率密度f(θk)分别为 ${P_k} = {n_k}/\sum\limits_{k = 1}^m {{n_k}} .$ (3) ${f_\theta }({\theta _k}) = \frac{m}{{2{\rm{\pi }}{P_k}}}.$ (4) 式中：x=Fv(v)、y=Fθ(θ). Johnson等[10]建议采用谐波函数对风向角直方图进行拟合, 为防止概率密度函数出现负值或旁瓣过大, 建议同时引入Bartlett窗函数作为权重函数.这种方法需要将与方位角数相等的谐波进行叠加, 因此, 其数学形式复杂, 且不具有对风向分布的普适性.实际上, 对风向的圆周分布, 目前应用更为广泛的是混合von Mises分布[14].von Mises分布常被称为圆周正态分布, 是最主要的用于描述方向数据的一种模型，其概率密度函数为 $f\left( \theta \right) = \frac{{{\rm{exp}}[\kappa {\rm{cos}}\left( {\theta - \mu } \right)]}}{{2{\rm{\pi }}{I_0}\left( \kappa \right)}}.$ (5) 式中：-π＜μ≤π, κ>0；I0(κ)为零阶修正贝塞尔函数, 计算式为 ${I_0}\left( \kappa \right) = {\left( {2{\rm{\pi }}} \right)^{ - 1}}\smallint _0^{2{\rm{\pi }}}{\rm{exp}}\left[ {\kappa {\rm{cos}}\left( {\theta - \mu } \right)} \right]{\rm{d}}\theta .$ (6) 研究表明, 风向角变量的概率密度分布一般具有多峰值, 单一von Mises分布不能完整描述风向圆周分布特性, 通常采用多阶混合von Mises分布表示风向角变量的分布[18], 其具体计算式为： $f\left( \theta \right) = \sum\limits_{i = 1}^c {({\omega _i}{f_i}\left( \theta \right)).} $ (7) 式中：c为混合von Mises分布的阶数, 根据风向峰值分布特征确定；ωi为各阶von Mises分布的权重系数；φi={ωi, μi, κi}为各阶参数. 目前, 关于von Mises分布的参数估计方法比较多, 如MSBC算法[19]、SU算法[18]等.然而, 这些迭代算法计算过程均较为繁琐, 且涉及手动查表, 实用性较差.因此, 本文采用高效简捷的LM算法(Levenberg-Marquardt algorithm)进行参数估计.为保证良好的收敛性, 需要对各参数进行初始赋值, 具体赋值参见文献[14]. 1.3 基于最大熵原理的联合分布 极值风速风向的联合分布属于角度-线性分布(angular-linear distribution, AL), Johnson等[10]基于最大熵原理导出了AL分布的概率密度结构, 可有效用于描述风速风向的联合分布.假设fv(v)和fθ(θ)分别为极值风速和相应风向的概率密度函数, 对应的分布函数为Fv(v)和Fθ(θ), 则根据最大熵原理可导出极值风速风向的JPDF形式如下式所示： ${f_{v,\theta }}\left( {v,\theta } \right) = 2{\rm{\pi }}g\left( \xi \right){f_v}\left( v \right){f_\theta }\left( \theta \right),$ (8) $\xi = 2{\rm{\pi }}\left[ {{F_v}\left( x \right) - {F_\theta }\left( \theta \right)} \right]; - \infty ＜ x ＜ \infty .$ (9) 式中：g(·)是角变量ξ的函数, 本文采用二阶混合von Mises函数形式予以表示. 从式(8) 可以看出, 最大熵原理实际上是将风速风向的联合分布函数与其各自的边缘分布函数连接在一起.早在1959年, Sklar指出：可以将一个联合分布分解为多个边缘分布和一个Copula函数, 这个Copula函数描述了变量间的相关性[21].根据该理论, 极值风速变量与风向变量间的Copula函数存在如下关系： $\frac{{\partial {C^2}}}{{\partial v\partial \theta }} = 2{\rm{\pi }}g\left( \xi \right).$ (10) 式中：C= C[Fv(v), Fθ(θ)], 为极值风速风向变量间的Copula函数. 1.4 拟合优度检验 为检验文中JPDF的拟合优度, 采用确定系数计算理论值与实际值的差异： ${\eta ^2} = 1 - \frac{{\sum\limits_{k = 1}^m {{{({P_k} - {T_k})}^2}} }}{{\sum\limits_{k = 1}^m {{{({P_k} - T)}^2}} }}.$ (11) 式中：Tk为理论累积频度值；T为Tk的平均值.η2介于0~1, 其值越大, 拟合效果越好. 2 极值风速风向实测数据 本研究的极值风速风向样本源自我国某地1971年1月1日——2000年12月31日全部360 m的月最大风速和相应风向记录(以正北方向为0°, 顺时针为正).月最大风速指每个月内10 min风速样本中的最大值, 风向记录包含16个方位角.极值风速频度分布和风向玫瑰分别如图 1、2所示.从图中可以看出, 月极值风速主要集中在0~20 m/s以内, 对应风向记录主要在NNE和SSW方向, 具有典型的双峰值分布特征.因此, 本文采用二阶混合von Mises拟合风向分布函数, 即c=2. 图 1 月极值风速频度分布直方图 Fig. 1 Frequency distribution histogram of monthlyextreme wind speeds 图 2 月极值风速对应风向玫瑰图 Fig. 2 Wind direction rose of monthly extreme wind 根据实测数据, 计算得到实测风速风向的联合概率分布, 结果如图 3所示(图中θ表示风向).从图中可以看出, 联合概率分布函数主要集中在10~15 m/s风速区间, 风向峰值在25°和200°附近. 图 3 实测极值风速风向联合概率分布图 Fig. 3 Joint probabilistic distribution or measuredextreme wind speed and direction 3 极值风速风向的联合概率分布模型 对极值风速样本分别使用Gumbel分布和GPD分布建立模型, 拟合参数如表 1所示, 分布曲线如图 4所示, 图中Fv为极值风速的累积分布值.可以看出, Gumbel分布与月极值风速样本的实测分布拟合得很好, 计算拟合结果的确定系数接近于1, 说明由极大似然估计法确定的Gumbel分布与实际累积分布的差异很小；GPD分布风速分位值略大于Gumbel分布, 且与实测累积分布差异较大.因此, 采用Gumbel模型建立极值风速风向联合分布. 表 1 极值风速模型的拟合参数 Table 1 Fitted parameters in extreme wind speed models 图 4 极值风速分布模型与实测值的分布曲线对比 Fig. 4 Distribution curves' comparison of probabilistic models and measured data for extreme speed 分别采用谐波函数和二阶混合von Mises分布建立风向圆周分布模型, 各参数估计值如表 2所示, 拟合曲线如图 5所示.可以看出, 在峰值风向处, 混合von Mises分布的计算结果略大于谐波函数值.由确定系数计算得到的两类函数拟合优度分别为RF2=0.998和R2von=0.988, 说明这2种模型都能较好地表征风向圆周分布.下文将采用二阶混合von Mises分布模型建立极值风速风向联合分布. 表 2 风向圆周分布和JPDF下二阶混合von Mises分布的拟合参数 Table 2 Fitted parameters of second order mixture von Mises distributions for wind circumferential direction and JPDF 图 5 风向圆周分布模型分布曲线与实测数据对比 Fig. 5 Distribution comparison of circumferential distribution models and measured data for winddirection 根据极值风速和风向的边缘分布模型, 基于最大熵理论建立联合概率分布JPDF, 其中g(ξ)仍采用二阶von Mises函数, 拟合参数如表 2所示, 联合分布如图 6所示.对比图 3和图 6可以看出, 联合分布模型呈现2个明显的峰值, 与实测分布吻合良好, 但模型的连续、光滑处理导致实测数据局部“毛刺”现象消失. 图 6 风速和风向的联合分布模型拟合结果 Fig. 6 Fitted result of proposed joint probabilistic model for extreme wind speed and direction 4 结果与讨论 我国建筑结构荷载规范(GB 50009-2012) 规定采用极值Ⅰ型分布确定不同重现期下的基本风速.设重现期为N, 则基本风速UNG的计算式为 $U_N^{\rm{G}} = u - a{\rm{ln}}\left( {{\rm{ln}}\frac{N}{{N - 1}}} \right).$ (12) 式中:a、u的具体取值参见表 1. 由式(12) 计算得到的是不考虑风向角分布的基本风速.采用联合分布模型同样可计算得到全风向角下的基本风速UNJ, 其计算公式为 $1 - \frac{1}{N} = \smallint _0^{U_N^{\rm{J}}}\smallint _0^{2{\rm{\pi }}}{f_{v,\theta }}\left( {v,\theta } \right){\rm{d}}\theta {\rm{d}}v.$ (13) 理论上, 由式(13) 与式(12) 计算获得的解析解应满足UNJ=UNG；在特定风向角下, 联合分布模型计算的极值风速边缘分布应与由相应风向角的极值风速Gumbel分布模型一致, 即Fv, θ(v|θ)=G(v|θ).采用上述原理可检验联合分布模型建立过程是否正确.如图 7、8所示分别为全风向角(即不考虑风向角的概率分布)下不同重现期的基本风速曲线, 及风向角为NNE时的概率分布模型对比.从图中可以看出：1) 由联合分布模型JPDF计算得到的基本风速-重现期曲线与由Gumbel分布模型计算得到的基本风速-重现期曲线完全吻合, 说明参数优化过程正确有效；2) 当风向角为NNE时, 联合分布模型JPDF确定的极值风速边缘分布与实测结果及Gumbel分布模型几乎一致.上述检验结果表明：JPDF模型能够有效表征实际风速风向的概率分布特征. 图 7 全风向角下, 联合分布模型和Gumbel模型不同重现期的基本风速对比 Fig. 7 Comparison of reference wind speeds in all directions vary with time obtained by joint distribution model and Gumbel model 图 8 风向角为NNE时, 不同模型的极值风速分布与实测结果的对比 Fig. 8 Comparison of extreme wind speed distributions in models and measured data when wind direction is NNE 现采用极值风速风向的JPDF分布模型计算不同风向角的基本风速.由定义确定不同风向角下基本风速U的计算式如下： $1 - \frac{1}{N} = \frac{{\smallint _0^v{F_{v,\theta }}\left( {v|\theta } \right)}}{{\smallint _0^\infty {F_{v,\theta }}\left( {v|\theta } \right)}}.$ (14) 式中：fv, θ(v|θ)表示风向角为θ时, 极值风速的条件分布函数.式(14) 较复杂, 难以获得解析解.本文采用数值迭代法求解基本风速. 如图 9所示为当重现期N=100 m时基本风速随风向角的变化规律.可以看出, 在不同风向角下, 基本风速具有一定波动, 其中风向角θ=200 °附近时, 基本风速值最大, 该结果与图 3所呈现处的现象一致.如图 10所示为当重现期分别为N=5, 10, 50, 100, 600 m时的基本风速玫瑰图, 图中同时列出由式计算得到的基本风速以示对比.从图中可知, 在30°~210°区间(顺时针), 由JPDF分布模型(式(14))计算得到基本风速值大于按规范计算得到的结果(式(12))；而在区间210°~30°(顺时针)内, 现象与之相反. 图 9 当重现期为100 m时, 联合分布模型的基本风速随风向的变化 Fig. 9 Reference wind speeds vary with wind direction in joint distribution model when return period is 100 months 图 10 当重现期为5, 10, 50, 100, 600时, 基本风速的玫瑰图 Fig. 10 Reference wind speed rose when return period is 5, 10, 50, 100, 600 months, respectively 值得注意的是, 虽然由联合分布模型JPDF计算得到的基本风速具有波动性, 但波动的幅度比较小.计算表明对50 a重现期即N=600 m时, 波动幅度小于5%.这说明, 计算样本中极值风速与风向的相关性较差.为验证上述结论, 下文对极值风速风向的相关性进行计算. 根据Sklar提出的Copula理论, Copula函数(式(10))表征2个变量间的相关性.基于Copula函数计算2个变量间的相关性测度, 采用Spearman秩相关系数ρ进行衡量，计算式[21]为 $\rho = 12\smallint _0^1\smallint _0^1xy{\rm{d}}C\left( {x,y} \right) - 3.$ (15) 式中：x=Fv(v)、y=Fθ(θ ). 根据Copula函数的拟合结果, 计算得到Spearman秩相关系数|ρ|=0.066 8, 说明极值风速和风向正向变化的一致性较差. Spearman秩相关系数仅对变量间的变化方向一致性进行度量, 并不能衡量变量间相关性程度.Mardia[22]提出采用线性-角度变量相关系数计算极值风速风向的相关性： ${r^2} = \frac{{r_{{\rm{vc}}}^2 + r_{{\rm{vs}}}^2 - 2{r_{{\rm{vc}}}}{r_{{\rm{vs}}}}}}{{1 - r_{{\rm{cs}}}^2}}.$ (16) 式中：rvc=fcorr (v, cos θ), rvs=fcorr (v, sin θ), rvc=fcorr (cos θ, sin θ), fcorr(·)表示序列间的相关系数. 计算结果表明, 文中采用的极值风速与对应风向间的相关系数为r2=0.001 3, 说明极值风速与对应风向不仅正向变化一致性较差, 其相关性程度也较弱.这一结论与由联合分布模型计算得到的基本风速风向间的相关性表现一致.因此, 极值风速风向联合分布模型能够较好地表征实际样本间的相关性.总体而言, 该模型能比较有效地预测极值风速风向的概率分布. 5 结论 (1) 二阶混合von Mises分布能较好地表征风向角的圆周分布, 基于最大熵原理建立的极值风速风向联合概率密度函数与实测值吻合良好. (2) 采用本文所建立的极值风速风向联合概率密度函数JPDF可以较有效地计算不同重现期下的基本风速, 且在特定风向角下的极值风速分布与建筑结构荷载规范(GB 50009-2012) 的分布规律相一致. (3) 由Spearman秩相关系数和线性-角度变量相关系数计算结果表明, 本文采用的极值风速与风向序列相关性较差, 符合极值风速风向联合概率密度函数确定的变量相关性, 验证了模型的有效性. 需要说明的是, 本文在建立极值风速风向的联合概率分布模型时, 由于非主风向角和高风速值的数据量较少, 对高保证率(如50年一遇、100年一遇)的基本风速值精度尚有待考证, 后续工作应对数据量更丰富的风速风向样本进行研究.