黄 鹏, 蓝鑫玥, 钟 奇

(同济大学 土木工程防灾国家重点实验室，上海 200092)

大量的风洞试验和现场实测数据表明[1-5]，低矮房屋处在大气边界层中湍流度较高的底部区域。屋面边、角、脊等位置受结构本身体型所引起特征湍流的影响，往往会产生幅值很大的风压“脉冲”。脉动风荷载产生的极值风压是低矮房屋破坏的重要原因，合理估计脉动风荷载产生的极值风压是计算低矮房屋风荷载的关键，是有效减少低矮房屋风致损坏及毁坏的重要途径。为此，研究人员基于零值穿越理论与经典极值理论提出了一系列的极值估计方法。

以Davenport[6]为代表的研究人员假定零均值的脉动风压服从标准高斯分布，根据零值穿越理论提出适用于高斯分布的脉动风压极值估计方法(下文称为Davenport法)。然而在高湍流风场的条件下或在风场的流动分离区，风荷载极值分布严重偏离高斯假设，采用Davenport法得到的极值估计值与实际结果差别较大。Sadek等[7]在Rice[8]零值穿越理论的基础上，假设时程的母体符合三参数 Gamma分布，应用Grigoriu[9]的“转换过程法(Translation Process Approach)”将非高斯风压时程的概率分布映射成标准高斯分布，由此提出了非高斯过程的极值计算方法(下文称为Sadek & Simiu法)。映射的标准高斯过程的概率分布函数中的零值穿越率必须由原非高斯时程获得，该过程对于非高斯较强的时程存在较大的误差。Quan等[10-12]以经典极值理论为基础，假定脉动极值风压时程满足极值Ⅰ型分布，利用样本的自相关系数确定母样本的分段长度，以子段的极值分布规律估算母样本的极值期望值(下文称为Quan et al法)。相对于零值穿越理论，该方法在理论上有很大的突破，对极值负压取得了较为满意的拟合效果。上述方法均以风洞试验的风压数据判断各种极值估计方法的合理性，而通过现场实测进行合理性验证的较少，此外，上述方法以零值穿越理论和经典极值理论下的极大值模型(Block Maxima Method, BMM)为基础。零值穿越理论的核心是假定脉动风压时程服从高斯分布，而极大值模型存在极值估计准确性受样本分布影响较大的问题。

因此，本文依托同济大学土木工程防灾国家重点实验室的风荷载实测基地，以15°屋面坡角低矮房屋在台风作用下的屋面风压实测数据为样本，基于经典极值理论下的超越阈值模型，提出一种以广义Pareto分布为拟合分布的脉动风压极值估计方法。首先确定阈值的上下限，将阈值等分为若干区间，并计算每一区间经验平均超额分布函数的斜率，以拟合残差平方和最小的方式选择合理的阈值区间，且将该区间第一个值作为最优阈值。本文分析了分段数的不同对拟合结果的影响，同时，探讨了本文方法在不同阈值选取方式中的优越性。此外，以10 min标准时距下观察风压极值的平均值作为标准极值，比较了本文方法所得极值估计值与标准极值的误差，并将本文方法与常用极值估计方法进行对比。

1 实测概况

实测基地依托同济大学土木工程防灾国家重点实验室，位于上海浦东国际机场附近沿海入海口区域，包括一栋屋面坡角可连续调节的足尺实测房，一座高约10 m和一座高约40 m的格构式测风塔[13]。其中，低矮房屋实测房长约10 m，宽约6 m，建筑檐口高度为8 m，代表了大部分村镇房屋的外形特征，屋面坡角可从0°～35°自由调节，满足对一般低矮房屋风荷载研究的需要。风压采集系统主要布置在屋面，其测点布置图如图1所示。

图1 浦东实测基地概况Fig.1 Pudong experimental base

本文使用实测风速样本为2016年第16号强台风“马勒卡”，其中心附近最大风力有14级(45 m/s)。下节将基于同济大学实测基地所采集到的在强台风“马勒卡”影响下低矮房屋15°坡角屋面风压结果，对典型位置关键测点的脉动风压时程进行分析。

2 脉动风压特性

结合“马勒卡”影响上海的实际时间，本文所选取的样本风速时程为2016年9月17日23时后48 h内实测基地10 m高度处的风速数据。

2.1 脉动风压分布特性

根据前人的研究成果，低矮房屋钝体形态明显，风荷载作用在低矮房屋的同时会形成明显的冲击、分离及旋涡脱落等特征湍流。因此，本文选取迎风面前端角部处测点74、背风面靠近屋脊位置测点36、背风面靠近山墙测点29及背风面后端中间位置测点6这4个具有代表性的测点进行脉动风压系数的研究。图2给出了典型测点在台风“马勒卡”影响下，平均风速为12.16～15.32 m/s湍流强度为0.35、湍流积分尺度为52 m、屋面坡度为15°时脉动风压系数随风向角的变化。

从图2 可以看出，当风向角在0°～50°内，背风面靠近屋脊位置测点36、背风面靠近山墙测点29的脉动风压系数随着风向角的增大而增大，且在风向角为30°～35°时数值增大速率较大。背风面后端中间位置测点6的脉动风压系数稳定在0.2～0.3，迎风面前端角部处测点74的脉动风压系数稳定在0.7～0.8。

可见，当风向角在0°～35°内，随着风向角度的增大，屋脊对来流的气动分离作用对背风面靠近屋脊位置和背风面靠近山墙区域屋面测点的脉动风压系数的影响逐渐加强，而整个屋面其他区域的脉动风压系数维持稳定。

图2 关键测点脉动风压系数随风向角的变化Fig.2 Curve of fluctuating wind pressure coefficient with wind direction angle at critical points

2.2 脉动风压概率特性

研究表明[14-17]，屋面风压具有非高斯性，特别是边、角、脊等典型位置的风压存在幅值较大的“脉冲”(见图3(a))，当然，屋面部分区域测点的非高斯性并不明显(见图3(b))。

(a)特征湍流作用区域测点

(b) 非特征湍流作用区域测点图3 实测风压系数时程Fig.3 Measured wind pressure coefficient time history

为直观分析低矮房屋屋面风压的非高斯性及概率分布特性，本文对10 min平均风速为13.56 m/s、平均风向角为39°的屋面关键测点实测脉动风压时程的概率分布采用Gauss分布、Lognormal分布、GEV分布与Gamma分布进行拟合，拟合结果见图4。

图4 关键测点风压系数时程概率分布特性及拟合效果Fig.4 Probability distribution characteristics and fitting effect of wind pressure coefficient time histories at critical points

从图4可以看到，测点74和测点34的概率分布呈现较为明显的负偏态，且均不服从Gauss分布、Gamma分布、Lognormal分布与GEV分布。在风压系数概率密度曲线的尾部即极端值区域的拟合效果并不理想。而测点67的风压实测数据与高斯分布拟合值相对误差不大，拟合效果优于测点74与测点34。其原因在于，测点74与测点34位于受特征湍流影响的边、角、脊位置，这些位置在旋涡脱落等湍流效应影响下，偏度系数与峰度系数绝对值较大(测点74的偏度系数为-1.82、峰度系数为8.16，测点34的偏度系数为-1.43、峰度系数为5.86)，风压具有较强的非高斯性，采用Gauss分布(标准高斯过程的偏度系数范围为-0.5～0.5，峰度系数范围为2.5～3.5)描述风压系数的概率特性是不合适的。而测点67位于迎风屋面中央，受特征湍流的影响小，风压基本围绕均值对称分布，非高斯性不显著，但Gauss分布对风压极值的拟合精度仍然不高。

由此可见，风压系数无法采用单一确定性函数较准确地描述其概率特性，风压极值的合理估计是比较困难的。以下采用基于经典极值理论下超越阈值模型的极值估计方法估计实测风压极值。

3 风压极值估计方法

3.1 超越阈值模型及广义Pareto分布

超越阈值模型以超过某一阈值的数据作为研究对象，有效地利用了有限的极端观察值，被广泛应用于极值分析、可靠性研究等领域。超越阈值模型分析超出阈值的数据(超阈值)的概率特性，即该模型主要针对母体样本分布的尾部区域。

选择适当的阈值，是准确估计广义Pareto分布形状参数与尺度参数的前提，也是正确分析数据的关键。研究人员发现，极值估计值对阈值非常敏感：阈值太

大，则样本数据容量有限，样本不具有代表性；阈值太小，则样本数据无法满足广义Pareto分布的要求，所得极值并不准确。目前，选取阈值的方法主要分为图解法和计算法两大类，计算法以峰度法及渐进均方误差法为代表。但上述两种方法计算效率较低，存在一定的局限性。

研究发现，超阈值服从广义Pareto分布。本节对广义Pareto分布进行介绍。

广义Pareto分布的概率分布表达式为

(1)

式中:μ∈R为位置参数;ξ∈R为形状参数;σ>0为尺度参数。

总体分位数函数为

(2)

极大似然估计是渐近正态分布的，且所得参数比较稳定，本文采用极大似然估计方法对形状参数ξ和尺度参数σ进行估计，具体方法如下：

超阈值的广义Pareto分布的概率密度函数为

(3)

式中:u为阈值。

那么对数似然函数为

(4)

分别对ξ和σ求偏导并令其等于0，可得

(5)

求解上述方程即可得到广义Pareto分布中参数ξ和σ的估计值。对于0

(6)

3.2 本文风压极值估计步骤

平均超额分布函数为超过阈值的样本均值，对于广义Pareto分布而言，其相应的表达式为

(7)

对于离散数据，可以定义经验平均超额分布函数为

(8)

本文方法基于两个原理。首先，若超阈值服从广义Pareto分布，则平均超额分布函数可描述为阈值的线性函数，如式(7)所示；其次，现有的经验平均超额分布函数图法(计算法)通过构造曲线肉眼选择使平均超额分布函数近似线性且斜率为正的阈值u0为合理阈值，属于定性判断，凭直觉和经验，存在主观性明显的缺陷，有时受样本数据结构影响甚至无法适用。本文提出以拟合残差平方和最小的方式对经验平均超额分布函数图法的合理改进，直接给出最优阈值u的大小，提高选取阈值的稳定性。

为了便于理解和采用，将此方法的计算步骤总结如下：

步骤1 根据浦东实测基地采集到的脉动风压时程数据，选择合适的阈值下限u1和上限un；

步骤2 将阈值区间[u1,un]进行等分，以k(k

步骤3 计算每一组数据中的第一个点和最后一个点连成直线的斜率，由定义有

(9)

步骤4 设q(li)为每一组数据对直线li的拟合残差平方和，由定义有

(10)

步骤5 选取拟合残差值最小的那组数据，并在其中选取第一个点的横坐标对应的脉动风压大小作为最优阈值；

3.3 分段数的确定

从拟合残差公式得到的参数转换关系式(10)可以看出，分段数(n-k+1)越大，每段中的数据越少，则越能反映经验平均超额分布与阈值之间线性关系，即所得到的拟合残差值越具有说服力。但是随着分段数的增大，则相邻段子样本之间的相关性较强，且明显地增大计算时间，降低计算效率。因此，有必要研究最优阈值与k值之间的相互关系。

(1)以现有方法的最优阈值为标准

图5给出了屋面15°坡角低矮房屋屋面典型测点脉动风压分别采用本文方法、峰度法及渐进均方误差法得到的最优阈值与k的关系曲线。本文对各测点10 min标准时距下的脉动风压时程样本采用本文方法将阈值区间分别按3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20个数据为一组等分得到最优阈值，并以峰度法及渐进均方误差法下的最优阈值作为参考。由图5可以看出，当k的值为5时，利用本方法得到的最优阈值与渐进均方误差法得到的最优阈值比较接近。更多测点的计算结果同样也证实了这一结论的正确性。

(a)测点34

(b)测点67图5 不同阈值选择方式最优阈值的比较Fig.5 Comparison of different threshold selection methods

(2)以K-S检验法作为拟合优度的评判标准

为进一步确定本文方法的合理性，本文选取Kolmogorov-Smirnov检验法作为判断超阈值是否服从广义Pareto分布的标准，并引入拟合优度判断指标来合理区分不同情况。其值为0时表明在显著型水平0.05下接受原假设，即认为超阈值服从广义Pareto分布；其值为1时表明在显著型水平0.05下拒绝原假设，即认为超阈值并不服从广义Pareto分布；其值为2时表明由于阈值较高，超阈值样本数量不足以与广义Pareto分布进行比较。

由Pickands[18]的理论研究成果可知，若最大值Mn近似服从广义极值分布，则超阈值X|X-u与超出量X-u近似服从广义Pareto分布，因此，可利用GPD(Generalized Pareto Distribution)拟合超阈值，从而利用超阈值分布求总体分布函数。

从图6中可以看出，当k=5时低矮房屋屋面各测点基于本文方法得到的超阈值与广义Pareto分布拟合程度较好，可以认为超阈值服从广义Pareto分布。更多的拟合结果同样也证实了这一结论的正确性。

图6 超阈值与广义Pareto分布的拟合效果Fig.6 Fitting effect of data beyond threshold and generalized Pareto distribution

3.4 不同阈值选取方式的比较

峰度法、渐进均方误差法是目前较为常用的超越阈值模型极值估计方法。本文采用标准时距(台风工况采用10 min)下脉动风压时程的观察极值作为标准极值，将上述两种方法及本文方法所得脉动风压极值估计值与标准极值进行比较，以标准误差率均值来衡量屋面关键测点在不同阈值选取方式下估计极值与标准极值之间的接近程度。

标准误差率平均值的表达式为

(11)

式中：M为标准时距样本数；N为所包含测点总数。

图7表示的是台风作用下超越阈值模型不同阈值选取方式所得极值估计值与低矮房屋屋面实测脉动风压的标准极值的比较。从误差率均值的角度而言，渐进均方误差法与本文方法所选取的阈值比较可信，而峰度法得到的阈值其估计结果与真实值误差率在-0.1～-0.5，并不稳定。而渐进均方误差法由于采用最小二乘法估计参数值，因此，在计算效率上比本文方法低得多，在工程及风洞试验的数据处理时其优越性并不突出。

3.5 不同极值估计方法的比较

基于零值穿越理论的Davenport法、Sadek & Simiu法以及基于经典极值理论的Quan et al法是目前较为常用的风压极值计算方法。本文用上述方法得到屋面测点脉动风压极值估计值，并与标准极值进行比较。

本文提出以拟合残差平方和最小的方法对经验平均超额分布函数图法进行改进，可直接得到最优阈值u，并将超阈值与广义Pareto分布进行拟合，以拟合广义Pareto分布的期望值作为估计极值，并与标准极值进行比较。

(a)估计极值与标准极值

(b)误差率均值图7 不同阈值选取方式的比较Fig.7 Comparison of different threshold selection methods

对于标准极值的选取，陶玲等[19-20]将多次采样获得的样本极大值的均值作为标准极值用于比较各种极值估计方法的准确性，并论证了选取均值作为比较标准是偏于保守的。本文根据上述结论截取浦东实测基地实测房屋面2 h脉动风压时程数据，以10 min为标准时距，将12个标准时距中的观察极值的期望值作为标准极值，以标准误差率均值来考察不同方法估计极值的准确性。

图8表示台风作用下常用极值估计方法所得极值估计值与低矮房屋屋面实测脉动风压的标准极值的比较。从图8中可知，Davenport法所给出的极值估计结果(绝对值)比标准极值小，误差均值基本集中在-0.4～-0.6，即采用Davenport法所拟合得到的脉动风压极值会极大地低估风压对结构的作用，特别对于极端值区域偏离程度更大；Sadek & Simiu法的极值估计效果优于Davenport法，误差均值基本保持在-0.4～0.1，离散程度较大，表明极值估计的效果不太好。Quan et al法给出的极值估计值结果较均匀地分布在标准极值上下，拟合结果较好，但离散程度仍较大。而本文方法离散程度最小，受来流扰动的影响较小，且极值估计值与标准极值误差率基本控制在5%以内。

表1给出了低矮房屋屋面特征湍流影响区域(迎风面前端角部、背风面靠近屋脊角部)关键测点标准极值与估计极值及其误差率。可见，Davenport法由于假定风压系数服从标准高斯分布，与实测风压系数的概率分布相差较大，因此，非高斯性较强区域风压系数极值估计值的误差也很大。Sadek & Simiu法以及Quan et al法极值估计结果相比Davenport法有较大改善，但对于非高斯性较强的风压系数极值估计仍不太理想。

(a)标准极值与估计极值

(b)误差率均值图8 常用极值估计方法的比较Fig.8 Comparison of common extreme value estimation methods

本文方法对风压较大区域的极值估计结果与实测值的误差率均值控制在-5%以内，因此，采用本文方法对低矮房屋屋面风压极值估计是可行的、有效的。

表1 标准极值与估计极值的误差Tab.1 Error of standard extreme values and estimated extreme values

4 结 论

通过对台风作用下位于上海浦东国际机场附近实测基地的实测房15°坡角屋面脉动风压的研究，采用极值理论下的超越阈值模型，并以广义Pareto分布作为风压极值的拟合分布合理估计屋面风压极值，可得到以下结论：

(1)Davenport法对脉动风压极值的估计结果误差较大，误差率均值达到40%～60%，因此，采用Davenport法作为极值估计方法会显著低估脉动风压极值对结构的影响；而Sadek & Simiu法对于非高斯性较强的时程所得极值估计值存在较大误差；Quan et al法虽然在一定程度上减小估计极值与标准极值的误差，但方法本身并不稳定，部分测点的误差率仍在10%以上。

(2)本文方法是对超越阈值模型的有益尝试，以广义Pareto分布作为风压极值的拟合分布得到的极值估计值与标准极值误差率控制在5%以内，特别对极值风压较大区域的估计结果较为理想，体现了该方法在现场实测条件下对屋面脉动风压极值估计准确性。

(3)超越阈值模型对阈值比较敏感，基于不同阈值选取方式得到的极值估计值相差极大。从误差率均值的角度分析，渐进均方误差法与本文方法所得到的最优阈值比较可信，而峰度法所得到的最优阈值无法作为超越阈值模型的合理阈值。