基于非高斯仿真的風壓系數極值計算方法

李壽科，李壽英，陳政清，孫洪鑫

（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年生