0 引言

作用在建筑物表面的风压极值计算一直是风工程中倍受关注的问题，研究人员和工程师们将得到的风荷载极值用以验算如幕墙等围护结构的局部安全性。由于风荷载对于高层建筑的影响尤为突出，且高层建筑表面存在较多的围护结构或装饰结构，这些围护结构的存在对于表面风场分布有较大的影响，因此对高层建筑结构表面风压极值的计算就显得尤为重要。

在风洞试验中，普遍应用的风荷载极值计算基本方法是Davenport^[1]提出的基于高斯分布假定的峰值因子法。该方法基于脉动风压服从高斯分布为理论基础，通过将平均风压加上峰值因子与脉动风压均方根的乘积得到风压的极值。对于建筑物表面的风压分布，有些区域风压服从高斯分布，但有些区域表面(如气流分离区)风压分布表现出明显的非高斯特性，峰值因子法的假定就不成立，此时峰值因子法的应用就受到很大的限制。

对于偏离高斯分布的风压极值计算方法, 研究人员采用很多方法做了改进。其中具有代表性的有：Kwon等^[2]和Yang等^[3]分别先后提出了基于Hermite多项式^[4]模型的改进峰值因子法；Simiu等^[5]利用三参数Gumma分布和高斯分布拟合非高斯时程尾部的方法；Cook^[6]提出了计算极值风速的独立风暴法，简记为MIS法，Harris^[7-8]在此基础上对MIS法进行改进。改进峰值因子法需要转换公式满足单调性的条件，但峰度和偏度有时不满足该条件，限制了该方法的应用范围。Simiu等的映射方法在计算极值的过程中，采用了样本的零值穿越率代替高斯过程的穿越率，并没有从根本上摆脱高斯假定带来的误差。全涌等^[9]采用了经典极值理论，并根据风压时程样本的相关系数发展了Gumbel法，采用全部观察极值样本计算极值分布参数，充分利用了样本的内部信息，且保证了风压极值样本之间的独立性；Quan等^[10]采用广义极值理论扩展了该方法的实用范围。

本文基于上述研究，根据广义Pareto分布和广义极值分布之间的关系，得到新的极值估计方法。通过高层建筑风洞试验对高层建筑表面的风压进行多次采样，得到足够多的样本数据，对现有几种经典的极值计算方法进行比较，验证新方法的有效性。

1 基于Pareto分布风压极值计算方法 1.1 基于Gumbel分布的极值计算方法

已有的研究表明，Gumbel分布由于拟合所需要的数据较少，且符合工程的要求，因此经常采用Gumbel分布拟合极值分布。假定一个时间序列X，其概率密度函数如式(1)所示。

$ F\left( {x;t} \right) = \exp \left\{ { - \exp \left[ { - \left( {\frac{{x - {b_t}}}{{{\sigma _t}}}} \right)} \right]} \right\} $

(1)

其中，σ_t和b_t分别为观察时距t对应下的尺度参数和位置参数，x和t分别为风压系数和观察时距。

为了拟合上述参数，传统方法需要大量的样本，先从N个观察时距为t的样本中获取N个观察极值，进而获取极值的概率分布函数。为了在有限的试验样本中获取足够数量的观察极值，Cook等^[11]提出将不同观察时距下参数转化的方法，对于不同观察时距t和T下的位置参数和尺度参数，其关系如式(3)、式(4)所示。全涌等^[9]根据时间序列的自相关系数与时间延迟的关系，确定了最小观察时距，这样就可以充分利用有限的数据拟合得到可靠的参数，从而得到极值的概率密度函数，进而得到不同保证率下的风压极值系数。

$ {F_T}\left( x \right) = F_t^n\left( x \right) $

(2)

$ {b_T} = {b_t} $

(3)

$ {\sigma _T} = {\sigma _t} + {b_t}\ln \left( {T/t} \right) $

(4)

1.2 极值概率分布模型的建立

通过式(2)我们可以知道，在短时距拓展到长时距的过程中，长时距下的极值分位数对应于短时距下的极值高分位数，时距拓展的越长，极值对短时距下的极值分布尾部越敏感。图 1为高层建筑2个不同测点的风压极值分布，当采用Gumbel分布拟合时，由于兼顾了全部的极值数据，导致对极值分布高尾部分估计存在误差，导致对长时距下极值估计的偏差放大，因此，可以采用更为合理的方法对极值分布尾部做出更为准确的估计。

根据经典极值理论，对于式(1)，即Gumbel分布，设有足够大的临界点u，对于X>u的情况，有(X-u)服从广义Pareto分布Ⅰ型分布，两者存在如下关系：

$ G = 1 + \lg \left( H \right) $

(5)

式中，H代表广义极值分布Ⅰ型分布，即Gumbel分布，G代表广义Pareto分布Ⅰ型分布。

基于广义极值Ⅰ型分布和广义Pareto分布中的Ⅰ型分布之间的关系，本文采用广义Pareto分布中的Ⅰ型分布来拟合极值分布的尾部数据，从而根据有限的数据得到风压极值的合理估计，在此简记为S01法。在Gumbel分布的基础上，对于式(1)，设定足够大的临界点u，对于X>u的情况，有(X-u)服从Pareto分布：

$ \begin{array}{l} G\left( {X \le x\left| {X \ge u} \right.} \right) = \frac{{P\left( {u \le X \le x} \right)}}{{P\left( {X \ge u} \right)}}\\ \;\;\;\;\; = \frac{{{F_X}\left( x \right) - {F_X}\left( u \right)}}{{1 - {F_X}\left( u \right)}} = 1 - \exp \left( { - \frac{{x - \mu }}{\sigma }} \right) \end{array} $

(6)

其中，u为阈值，σ、μ分别为尺度参数和位置参数。

根据式(6)，确定极值的概率分布函数为：

$ {F_X}\left( x \right) = {F_X}\left( u \right) + \left[ {1 - {F_X}\left( u \right)} \right]\left[ {1 - \exp \left( { - \frac{{x - \mu }}{\sigma }} \right)} \right] $

(7)

然后根据式(2)建立短时距与标准时距之间的关系，从而得到指定保证率p对应的极值x_p：

$ {x_p} = - \ln \left( {1 - \frac{{{p^{\frac{1}{k}}} + {F_X}\left( u \right)}}{{1 - {F_X}\left( u \right)}}} \right)\sigma + u $

(8)

上述公式用于风压极大值的计算。对于风压极小值，对风压样本取相反数，按求取极大值的方法得到风压的极小值。

1.3 极值分布参数的确定

基于Pareto分布的极值计算方法，需要确定式(8)中阈值u、尺度参数σ的值。对于阈值u的选取，本文采用Hasofer^[12]建议取值：

$ {k_{{\rm{Tail}}}} = 1.5\sqrt k $

(9)

其中：k为全部的观察极值的个数，k_Tail为尾部数据的个数，对应于阈值u的位置。阈值u确定以后，可以得到：

$ {F_X}\left( u \right) = \frac{{{k_u}}}{{k + 1}},其中\;{k_u} = k - {k_{{\rm{Tail}}}} $

(10)

然后对式(7)进行变换得到关于分布参数的线性函数：

$ - \ln \left( {\frac{{1 - {F_X}\left( x \right)}}{{1 - {F_X}\left( u \right)}}} \right) = ax - b $

(11)

其中：$a = \frac{1}{\sigma },b = \frac{\mu }{\sigma }$。

对选取的高于阈值的k_Tail个独立峰值按由小到大进行排序，记为x_v，对应的可以得到F_X(x)的估计量：

$ {p_v} = v/\left( {k + 1} \right),\left( {v = {k_u} + 1,{k_u} + 2, \cdots ,k} \right) $

(12)

式(12)采用经验概率分布函数代替真实概率分布函数，这一过程会产生系统误差。在求取未知参数前，为消除系统误差，进行$ - \ln \left( {\frac{{1 - {F_X}\left( x \right)}}{{1 - {F_X}\left( u \right)}}} \right)$的无偏估计。根据文献[15]，令y=F_X(x), 则对v阶$ - \ln \left( {\frac{{1 - {F_X}\left( x \right)}}{{1 - {F_X}\left( u \right)}}} \right)$的无偏估计函数$\overline {{\xi _v}} $表达式为：

$ \overline {{\zeta _v}} = \int_0^1 {\left[ {\left( {1 - y} \right)/\left( {1 - {F_X}\left( u \right)} \right.} \right]{f_v}\left( y \right){\rm{d}}y} $

(13)

式(13)中f_v(y)的函数表达式为：

$ {f_v}\left( y \right) = \frac{{k!}}{{\left( {v - 1} \right)!\left( {k - v} \right)!}}{y^{v - 1}}{\left( {1 - y} \right)^{k - v}} $

(14)

通过对式(11)进行最小二乘拟合即可得到分布参数的估计量使差平方和s²最小。对于平方和中每一项的作用大小，采用权系数w_v来调整各项在平方差中的权重。差平方和的表达式如下：

$ {s^2} = \sum\limits_{v = {k_u} + 1}^k {{w_v}{{\left( {\overline {{\zeta _v}} - a{x_v} + b} \right)}^2}} $

(15)

权系数${w_v} = \frac{1}{{\sigma _v^2}}/\sum\limits_{v = {k_u} + 1}^k {\frac{1}{{\sigma _v^2}}} $，其中σ_v为ζ_v的标准差，表达式为${\sigma _v} = {\left[ {\overline {\xi _v^2} - {{\left( {\overline {{\xi _v}} } \right)}^2}} \right]^{\frac{1}{2}}}$。

为了使式(15)中对应于自变量a、b的函数取最小值，分别令自变量的一阶偏导数为0，得到关于自变量a、b线性方程组，然后根据线性方程组即可得到估计量：

$ a = \frac{{\sum\limits_{v = {k_u} + 1}^k {{w_v}\overline {{\zeta _v}} {x_v}} - \left( {\sum\limits_{v = {k_u} + 1}^k {{w_v}\overline {{\zeta _v}} } } \right)\left( {\sum\limits_{v = {k_u} + 1}^k {{w_v}{x_v}} } \right)}}{{\sum\limits_{v = {k_u} + 1}^k {{w_v}x_v^2} - {{\left( {\sum\limits_{v = {k_u} + 1}^k {{w_v}{x_v}} } \right)}^2}}} $

(16)

$ b = a\sum\limits_{v = {k_u} + 1}^k {{w_v}{x_v}} - \sum\limits_{v = {k_u} + 1}^k {{w_v}\overline {{\zeta _v}} } $

(17)

2 风洞试验概况

试验在湖南大学风洞实验室的HD-3大气边界层风洞中进行。试验段长10m，截面宽2.5m，高3m，转盘直径1.8m，试验段风速0~20m/s连续可调。试验模型为厦门沿海某高层建筑，测点沿模型高度分层布置，除3层布置22个测点外，其余层布置15个测点。

试验采用A类风场, 试验的几何比尺、风速比尺和时间比尺分别为1/200、1/8和1/25。风压数据采样频率为312.5Hz。试验中共进行80次重复独立采样，每个样本的数据为7500个，对应于实际采样时间600s, 即10min。图 2为侧风面测点111的风压系数样本。

各测点的风压系数由下式计算：

$ {C_{pi}}\left( t \right) = \frac{{{P_i}\left( t \right)}}{{0.5\rho V_\infty ^2}} $

(18)

式中：C_pi(t)为试验模型上第i测压孔所在位置的风压系数时间序列；P_i(t)为该位置处测得的表面风压的时间序列；V_∞为参考点高度处的来流平均风速。

3 基于高层建筑风洞试验数据的极值计算结果比较

已有研究表明，常用的峰值因子法对于非高斯风压的估计明显存在误差，因此本文只对改进峰值因子法、Sadek-Simiu法、改进Gumbel法(Quan法)和本文方法进行比较，这四种方法依次简记为KWO、SAD、H01和G01。

1) 对于H01法，通过对试验数据的分析，确定最小观察时距为4s，因此从每个风压样本中得到150个独立峰值风压进行拟合分析。

2) 采用KWO法计算时，要求偏度和峰度满足要求：3－r₄+(1.25r₃)²≤0，其中r₃、r₄分别为风压样本的偏度和峰度。对于不符合要求的点，通过调整偏度使3－r₄+(1.25r₃)²=0，满足该方法要求。

3) 对于SAD法，计算三参数Gamma分布有极大似然估计法和矩估计法，由于两者差别不大^[5]，因此本文采用矩估计法。

将试验得到的80个样本的观察极值x_i进行排序，并得到经验极值概率分布函数F_emp：

$ {F_{{\rm{emp}}}} = i/\left( {k + 1} \right),其中\;k = 80 $

(19)

根据经验极值概率分布函数式(19)得到指定保证率下的极值，并将此极值作为标准极值，同时将观察极值的平均值作为标准的极值期望值。

首先，验证上述四种方法对于风压极值计算的准确性。图 3中分析的测点为高层建筑角部边缘的侧风面测点111，具有明显的非高斯特性，分别采用上述四种风压极值计算方法对该点的80个采样样本进行计算。由于Gumbel分布的均值对应于其57%分位数值，因此对保证率57%^[15]下的计算极值与标准值进行对比，其结果如图 3所示。图 3中的红色横线为通过式(19)得到的57%保证率极值的标准值，分别采用箱形图表示通过上述四种方法计算得到的极值，从图中可以明显看到上述四种方法对于极值计算的准确性。对于极大值，采用SAD和G01计算得到的极值更接近于标准值，H01和KWO则明显偏离标准值；对于极小值，采用G01、H01和KWO计算得到的极值接近于标准值，其中G01最为接近，而SAD则偏离较为严重。

其次，通过全部测点来验证计算方法的稳定性。对计算极值和标准值之间的误差进行量化处理，本文定义平均标准差ASD^[16]作为衡量标准：

$ ASD = \frac{1}{k}\sum\limits_{i = 1}^k {\sqrt {{{\left( {C_i^{{\rm{esm}}} - C_i^{{\rm{em}}}} \right)}^2}} } $

(20)

式中，C_i^esm代表第i个测点计算极值的平均值，C_i^em代表第i个测点的标准值，计算结果如表 1所示。

从表 1中可以发现，采用本文方法得到的极值和标准极值最为接近，平均标准差ASD最小，H01法次之，SAD法误差最大。

4 结论

通过对高层建筑风洞试验多次采样得到的数据，提出了用于计算高层建筑表面风压极值的计算方法。主要结论有：

1) 对于高层建筑整体的风压，由于选用的测点比较多，且测点包含高斯与非高斯风压，采用本文方法可以获得更为准确的风压极值，研究方法具有一般性的意义，可靠性较高。

2) 本文只是对高层建筑表面的风压进行研究，对于其他建筑类型表面的风压，尚需要更多试验的验证。