李寿科，李寿英，陈政清，孙洪鑫

（1.湖南科技大学土木工程学院，湖南 湘潭 411201；2.北京交通大学结构风工程与城市风环境北京市重点实验室，北京：100044；3.湖南大学风工程试验研究中心，长沙 410082）

基于非高斯仿真的风压系数极值计算方法

李寿科1，2，李寿英3，陈政清3，孙洪鑫1

（1.湖南科技大学土木工程学院，湖南 湘潭 411201；2.北京交通大学结构风工程与城市风环境北京市重点实验室，北京：100044；3.湖南大学风工程试验研究中心，长沙 410082）

以多变量相关非高斯过程仿真方法为基础，发展了一种基于单次采样的多变量非高斯仿真极值计算方法。首先介绍开孔屋盖的风洞试验概况和多变量相关非高斯过程仿真的基本理论，对屋盖上一组测点风压进行了非高斯仿真，结果表明基于谱修正的多变量相关非高斯过程仿真方法得到的时程在功率谱密度，相干函数，高阶矩三方面与目标值接近，仿真效果较好，然后采用经典极值理论对多次仿真的非高斯时程进行极值计算，将该方法得到的峰值因子与以往常用方法的结果进行比较，结果表明：Davenport峰值因子法高估气流分离区左偏风压的正峰值因子60%，低估负峰值因子43%；Sadek－Simiu峰值因子法低估了高峰度风压的峰值因子50%；而基于单次样本进行仿真的非高斯仿真峰值因子法，其估计的开孔屋盖的峰值因子最为准确，与观察峰值因子总体上最为接近。

非高斯仿真；风压极值；峰值因子；风洞试验

根据风灾调查，建筑围护结构的破坏在风灾破坏损失中占很大的比例。建筑围护结构的设计由其表面的极值风荷载决定，通常得到极值风荷载的方法主要有两类。第一类为基于多次独立采样的经典极值理论方法，该方法对每次采样的极值进行分布拟合，从而确定出具有一定保证率的极值风压。该方法对非高斯过程和高斯过程具有同样的适应性，被认为是较准确的一种极值统计方法，但此种方法需对实际信号进行多次重复独立的采样，将耗费较大的人力和物力，在实际过程中却往往较少采用。第二类为基于单次采样的零值穿越理论的方法。对于服从高斯分布的随机信号，Davenport［1］在零值穿越理论的基础上，获得了服从窄带高斯分布随机信号的峰值因子，此种方法使用简单方便，在早期工程实际中被广泛采用。但很多时候风压并不服从通常的高斯分布，风压的非高斯特征对风压极值的估计具有较大的影响，Holmes［2］及Gioffre［3－4］对此做了相关研究，结果表明来流分离区风压的非高斯特性对风压极值的影响很大，按照高斯假定估计的峰值因子明显偏小。非高斯风压具有比高斯风压更大的破坏性，Holmes［2］发现非高斯分布风压导致的结构破坏要比高斯分布的风压大15%～30%，所以对于屋盖表面风压极值的估计需考虑其概率分布特性，兼顾其非高斯特性。对于非高斯过程的极值风压，Sadek［5］以Gamma分布和高斯分布作为母体分布，提出了一种非高斯转换过程的极值风压计算方法，而Kareem等［6］在Davenport方法的基础上，通过将非高斯过程展开为标准高斯变量的Hermite级数，将仅适用于窄带高斯过程的基于零值穿越理论的峰值因子法扩展到了非高斯过程。全涌等［7］提出了一种对单次标准长度的非高斯风压时程数据进行分段，然后通过子段的极值分布规律估算出母段的期望极值的改进经典极值法，在此方法中选择合适的采样长度较为重要。

Gurley［8］在进行非高斯风压场仿真研究时采用了Hermite级数来描述非高斯风压的概率分布，结果表明Hermite级数可以较好的体现随机过程的高阶矩特性从而达到可描述非高斯过程概率密度函数的目的。本文针对多次独立采样的经典极值方法需耗费大量人力物力的缺点，研究多变量的相关非高斯随机过程仿真方法，提出一种基于单次采样结合多变量相关非高斯随机过程仿真的峰值因子估计方法。最后，基于一个屋盖开孔建筑的刚性模型测压试验数据，通过与以往几种方法的比较，显示本方法的进步之处。

1 刚性模型测压试验概况

试验在湖南大学HD－2风洞的高速试验段进行。试验模型采用有机玻璃制作，实际结构在主要尺寸上与TTU建筑（13.72×9.14×3.96 m）保持几何相似，在屋盖的中心进行15%的开孔，模型的几何缩尺比为1∶50，试验照片如图1所示。模型立墙表面布置56个测点，屋盖上下表面布置120个测点，上下表面测点位置对应，详细测点布置见图2。试验的风向角定义见图2，风向角间隔5°，共72个测试风向角。采样时长33 s，采样频率330 Hz，共采集10 000个数据点，所有采样后数据采用管道频响函数进行畸变修正，且对典型风向角0°、45°、90°、270°、315°进行10次独立重复采样。试验风速11.0 m/s，参考高度为8 cm，相当于实际高度4 m。试验模拟了B类地貌风场［9］，风场比例为1∶50，平均风剖面指数为0.15，湍流度剖面也与实际大气中的情况基本一致，10 m高度处为0.20左右。

图1 试验模型照片Fig.1 Photo of testmodel

图2 试验模型测点布置图Fig.2 Tap location of testmodel

2 数据处理方法

风压系数是结构风压的无量纲表现形式，测点i的风压系数CPi（t）定义如下：

其中：Pi（t）为风洞试验中压力扫描阀测得的风压时程；P0为风洞试验段处的静压，采用皮托管测得；ρ为空气密度，取ρ＝1.225 kg/m3；uh为屋盖最高点处的平均风速。CPi（t）的平均值CPi＿mean为平均风压系数，CPi（t）的脉动值CPi＿rms为脉动风压系数，极大值风压系数为)CPi和极小值风压系数为(CPi。风压系数的正负峰值因子定义为

上表面减下表面风压系数为净风压系数，根据式（2）可得相应的峰值因子。

3 基于多变量相关非高斯过程仿真的极值计算理论

3.1 基于谱修正的多变量相关非高斯过程的仿真

仿真多变量非高斯相关过程，需要保证仿真结果的前四阶矩、功率谱密度和相干函数与测量目标值一致，首先需要得到多变量的高斯相关过程，继而进行非高斯过程转换，其具体步骤如下。

（1）基于谱分解方法的多变量高斯仿真

已知n个相关过程的自谱和它们之间的互谱，其互谱密度矩阵可以用切比雪夫分解表示为：

式（3）中Gii（ω）为第i个变量频率点ω处的自谱，Gij（ω）为第i个变量和第j个变量频率点ω处的互谱，其中：

式（4）中Rij（ω）为互相干函数，参考Gurley的方法［8］，产生两个独立的0均值和Δω根方差白噪声序列η和ζ组成复数序列ξ＝η＋kζ，则这n个变量的傅里叶序列可以表示为：

对式（5）进行傅里叶逆变换则可得到n个相关的高斯时间序列Yi（t）。

（2）单变量基于谱修正方法的非高斯仿真

在得到n个相关的高斯过程后，需将n个高斯过程利用Hermite级数转换方法转换到非高斯过程，转换后的变量的前四阶矩和功率谱密度需要保证与目标变量值一致，如不满足则需要进一步迭代修正，具体步骤如下。

（ⅰ）向前非高斯转换

应用向前Hermite转换方法将高斯过程Y（t）产生非高斯过程Xng（t）：

转换后的非高斯变量满足目标值的前四阶矩，但其功率谱密度将被扭曲。

（ⅱ）功率谱修正

保持非高斯过程的相位不变，利用目标谱值进行修正形成新的非高斯序列：

式（7）中GT为目标功率谱值。修正后序列的偏度和峰度将可能不再满足目标矩，设定偏度和峰度的误差判别标准为：

设定一误差限，假如误差超过其误差限，则进行下一步骤（ⅲ），如果在误差限内则退出，则得到了多条相关性可能不满足的非高斯过程，继续步骤（ⅲ）。

（ⅲ）向后高斯转换

由于没有达到目标非高斯过程，需再次进行仿真迭代，则需将步骤（ⅱ）得到的非高斯过程采用式（6）的反函数进行Hermite高斯过程转换，然后再返回到步骤（ⅰ），进行迭代求解，直至达到误差限。将式（6）求反函数，可得u（x）为：

（3）相干函数修正

采用单变量非高斯仿真方法可以得到多条非高斯过程，但各个变量之间的相关性将得不到满足，此时可对其相干函数也进行迭代修正，其设定误差判别标准为：

式（11）中：RXng（f）为迭代前非高斯过程的相干函数；RT

ijij（f）为目标相干函数；tol为相干函数误差限。如果结果在误差限内则退出，如果结果不在误差范围内，更新相干函数，返回到开始步骤（i）继续仿真，其更新相干函数表达如下：

3.2 非高斯仿真峰值因子计算

为利用多变量相关过程非高斯仿真方法获得峰值因子，首先可利用多变量非高斯仿真方法对标准化测点时程进行多次仿真，可得到多条时程，对每条时程求得观察极值，对所有观察极值进行经典极值分布（Gumbel分布）拟合，则可得到一定保证率下的峰值因子，在后文中称为非高斯仿真峰值因子，进行平均即可得到与Davenport峰值因子法［1］、Sadek-Simiu方法［5］对应的峰值因子。如果仿真时程不进行标准化，则可直接得到相应的风压系数极大和极小值。

4 非高斯仿真极值算例

4.1 多变量相关非高斯过程仿真实例

为验证基于谱修正多变量相关过程的非高斯仿真的有效性，基于本文试验结果，对开孔屋盖0°风向角时的76、77、78测点的风压时程进行仿真，仿真测点布置如图2所示，设定误差上限为10%，当然，测点数越多，误差上限设定越小，需要迭代计算的时间越长。图3给出了76、77、78测点的风压系数功率谱和相干函数的仿真和试验结果比较。由图3可以看出，测点风压之间的相干函数仿真结果和试验结果在低频处非常接近，在高频处稍有偏差，测点风压的功率谱密度函数的仿真结果和试验结果比较一致。表1给出了仿真时程的高阶矩与目标矩之间的对比，从表1可以看出，仿真风压的偏度和峰值均与目标试验值较为接近，其误差在设定的误差限10%以内。由此可以看出，多变量的相关非高斯过程仿真方法可以较准确的仿真出多条相关的非高斯风压时程，其仿真结果的功率谱密度、相干函数以及高阶矩均在目标偏差范围内。

图3 三变量相关过程非高斯仿真结果和测量结果频域统计比较Fig.3 Comparison between simulation and target in frequency domain

表1 目标时程和仿真时程的2、3、4阶矩比较Tab.1 Com parsion of deviation，skewness，ku rtosis between simulation and target

4.2 非高斯仿真峰值因子法与以往峰值因子法的比较

图4给出了开孔屋盖在0°、45°、90°风向角时的横向中轴线测点净风压系数的峰度和偏度分布规律。由图4可以看出，对于0°风向角，迎风区测点风压左偏，尾流区测点风压右偏，边缘分离处测点风压出现大负偏度（－1.23）和高峰度（6.21），表现出明显非高斯特性，背风尾流区测点表现出中等非高斯（偏度约为0.5，峰度值大于4.2）；45°风向角时，测点风压偏度较0°风向角小，其相应的峰度值也减小；90°风向角时，屋盖横向中轴线测点风压偏度均为负，在三个典型风向角中达到最小，此时测点处风压分离程度最小，但峰度值明显偏离3，在三个典型风向角中（0°、45°、90°）为最大。由此可以看出，屋盖中轴线测点明显偏离标准正态高斯分布（偏度为0峰度为3），风压的偏度主要由气流的分离程度决定，气流分离越明显，其风压偏度越大。

图4 屋盖横向中轴线测点风压系数偏度和峰度分布规律Fig.4 Skewness and kurtosis ofwind pressures on opening roof

图5给出了开孔屋盖在0°、45°、90°风向角时的横向中轴线测点净风压在基于多变量相关过程非高斯仿真峰值因子法、Sadek-Simiu的转换过程法、Davenport峰值因子法和观察极值法下的正负峰值因子分布规律，各种方法对应的结果为非高斯仿真峰值因子、Sadek-Simiu峰值因子、Davenport峰值因子和观察峰值因子，其中非高斯仿真峰值因子为通过多变量相关非高斯过程仿真方法进行非高斯仿真16次取极值的平均值获得。在0°风向角时，迎风区分离测点风压左偏，即风压概率密度函数左边表现为长尾部，由图5（a）可以看出，Davenport峰值因子法低估了负峰值因子达43%，高估了正峰值因子达60%，非高斯仿真峰值因子稍大于观察峰值因子，而Sadek-Simu峰值因子则与观察峰值因子较为接近；估计负峰值因子时，Davenport峰值因子法明显低估了迎风区屋盖的负峰值因子，Sadek-Simu峰值因子小于观察峰值因子，非高斯仿真峰值因子与观察峰值因子较为接近。在45°风向角时，视风压为高斯分布的Davenport峰值因子与观察峰值因子的误差较大，Sadek-Simu峰值因子则小于观察峰值因子，非高斯仿真峰值因子与观察峰值因子较为接近；在90°风向角时，Davenport峰值因子明显低估负偏度高峰度的屋盖风压正负峰值因子，其偏差高达50%，而Sadek-Simu峰值因子则明显小于观察峰值因子，非高斯仿真峰值因子与观察峰值因子较为接近。

综合不同方法对开孔屋盖的峰值因子计算结果比较可以得出以下结论：①Sadek-Simiu峰值因子法对0°风向角的正峰值因子具有较好的估计，但在一定程度上低估了开孔屋盖多数风向角下的峰值因子，且对于高峰度的风压峰值有明显的低估；②Davenport峰值因子法由于忽略风压的非高斯特性，对开孔屋盖的左偏风压，高估正峰值因子，低估负峰值因子；③基于观察样本进行仿真的非高斯仿真峰值因子法，其估计的开孔屋盖的峰值因子最为准确，与观察峰值因子总体上最为接近。

图5 不同计算方法的峰值因子比较Fig.5 Comparison of peak factors by differentmethods

5 结 论

（1）基于谱修正的多变量相关非高斯过程仿真方法仿真出的多条非高斯时程的功率谱密度，相干函数，根方差，偏度和峰度与目标值吻合得很好，误差在目标范围内；

（2）开孔屋盖测点风压明显偏离标准正态高斯分布，风压的偏度主要由气流的分离程度决定，气流分离越剧烈，其风压偏度越大；

（3）Davenport峰值因子法由于忽略风压的非高斯特性，对开孔屋盖气流分离区的左偏风压，高估正峰值因子60%，低估负峰值因子43%；

（4）Sadek-Simiu峰值因子法对于高峰度的风压峰值有明显的低估，偏差高达50%；基于观察样本进行仿真的非高斯仿真峰值因子法，其估计的开孔屋盖的峰值因子最为准确，与观察峰值因子总体上最为接近。

［1］Davenport A G.Note on the distribution of the largest value of a random function with application to gust loading［J］.In：Proc.ICE，1964（28）：187－195.

［2］Holmes JD.Wind action on glass and Brown's integral［J］.Engineering Structures，1985，7（4）：226－230.

［3］GioffrèM，Gusella V，Grigoriu M.Non-Gaussian wind pressure on prismatic buildings.I：Stochastic field［J］.Journal of Structural Engineering，2001，127（9）：981－989.

［4］GioffrèM，Gusella V，Grigoriu M.Non-Gaussian wind pressure on prismatic buildings.II：Numerical simulation［J］.Journal of Structural Engineering，2001，127（9）：990－995.

［5］Sadek F E S.Peak non-gaussian wind effects for databaseassisted low-rise building design［J］.Journal of Engineering Mechanics-Asce，2002，5（128）：530－539.

［6］Kareem A，Zhao J，TognarelliMA.Surge response statistics of tension leg platforms under wind and wave loads：a statistical quadratization approach［J］.Probabilistic Engineering Mechanics，1995，10（4）：225－240.

［7］全涌，顾明，陈斌.非高斯风压的极值计算方法［J］.力学学报，2010（3）：560－566.

QUAN Yong，GUMing，CHEN Bin.Study on the extreme value estimationmethod of non-gaussian wind pressure［J］.Chinese Journal of Theoretical and App lied Mechanics，2010，42（3）：560－566.

［8］Gurley K R.Modelling and simulation of non-Gaussian processes［M］.University of Notre Dame，Notre Dame.1997：166－203.

［9］中华人民共和国建设部.建筑结构荷载规范50009－2012［S］.北京：建筑结构出版社，2012：36－50.

Prediction of wind pressure peak factor with non-gaussian simulation

LIShou-ke1，2，LIShou-ying3，CHEN Zheng-qing3，SUN Hong-xing1
（1.School of Civil Engineering，Hunan University of Science and Technology，Xiangtan 411201，China；
2.Beijing's Key Laboratory of StructuralWind Engineering and Urban Wind Environment，Beijing Jiaotong University，Beijing100044，China；
3.Wind Engineering Research Center，Hunan University，Changsha 410082，China）

Based on the Non-Gaussian Simulation of multivariate stochastic processes method，one prediction method for wind pressure extreme value was proposed.The wind pressure time histories of several opening roofs were simulated with multivariate Non-Gaussian simulation method based on wind tunnel test data.It was shown that power spectral density，coherence，deviation，skewness and kurtosis of simulated Non-Gaussian time histories are very close to the destination values.Then the peak factors ofwind pressures on the opening roofswere predicted from themultivariate non-gaussian simulation time histories for several times with the typical extreme value theory，and the results were compared with those several generalmethods.Itwas shown that Davenportmethod overestimates the positive peak factor by 60%，its skewness is negative，and it underestimates the negative peak factor by 43%；Sadek-Simiu method underestimates the peak factor by 50%，it has a higher kurtosis；the proposed method can predict the peak factor effectively，and the overall error is smallest.

non-Gaussian simulation；extreme value；peak factor；wind tunnel test

TU119＋.21

A

10.13465/j.cnki.jvs.2014.24.020

国家自然科学基金（51248001）资助；湖南省教育厅科学研究一般项目（14C0431）资助；湖南省高校创新平台开放基金（湘教通（2012）595号）资助

2013－09－23 修改稿收到日期：2014－01－02

李寿科男，博士，讲师，1981年生