全 涌， 張秉超， 顧 明

(同濟大學 土木工程防災國家重點實驗室，上海 200092)

Davenport[1-5]指出，極值風速的方向性對建筑風荷載有著重要的影響.Cook[6]是最早在處理風氣候數(shù)據(jù)時考慮極值風速的風向問題的.他基于每年可能發(fā)生上百次的獨立風暴記錄開發(fā)出一種計算各風向極值風速的穩(wěn)定的估算方法，該方法使得可用的樣本不僅限于年極值，使得樣本數(shù)量大大增加，每個風向下的風速極值規(guī)律可以進行獨立統(tǒng)計.之后Cook[7-8]又對這種方法進行了不斷更新改進，使之得到逐步完善.近年來，全涌等[9]對Cook方法的極值模型部分作出改進，使其由多元隨機變量更準確地轉化為一維隨機變量，從而得出更精確的結果.但是，此類方法不便于討論各方向極值風速的相關性.然而在此問題上，目前研究成果還很少.為了分析極值風速的風向相關性對總體風荷載評估的影響，Haraguchi和Kanda[10]以及Kanda和Itoi[11]分別提出了各風向極值風速的互相關Gumbel概率模型，但并沒有將它運用到極值估計中.Haraguchi和Kanda[10]比較了這兩個模型，指出上述兩種方法還需要進一步優(yōu)化.最近，Zhang和Chen[12]又提出一種應用高斯Copula函數(shù)的聯(lián)合隨機分布模型來考慮極值風速的風向相關性.這種方法用Copula函數(shù)將邊緣分布函數(shù)組合成聯(lián)合概率分布，得出多元極值風速模型.其后Zhang和Chen[13]考慮了使用不同的邊緣概率分布模型對預測結果的影響，并討論了風壓的風向因子特性.

另一方面，在臺風多發(fā)區(qū)，由于臺風與良態(tài)風的極端風速概率特性完全不同，仍然對混合氣候風速取樣并擬合，通常會得到不真實的結果.Gomes和Vickery[14]給出了一種基于混合氣候風速年極值估計方法，Cook，Harris和Whiting[15]對該方法進行了改進，Cook[16]進一步分析了根據(jù)該方法得到的混合氣候風速極值的置信區(qū)間.但由于臺風極值風速觀測數(shù)據(jù)不易獲得，目前還沒有考慮風向的混合氣候極值估計方法.現(xiàn)階段，臺風模擬的方法已經(jīng)比較成熟.比如，Meng[17-18]建立的臺風風場模型，并演化出了完整的解析解采用臺風模擬的方法進行極值估計已成為可能.

本文首先根據(jù)風速風向分布模型，給出一種新的考慮風向的極值風速標準值.然后用Monte Carlo方法分別獲取臺風和良態(tài)風氣候的風速極值.其中，臺風氣候極值風速采用Yan Meng模型模擬，良態(tài)風氣候則利用氣象站數(shù)據(jù)擬合極值分布后再模擬一定數(shù)量的極值風速數(shù)據(jù)，并對兩者進行組合.最后以上海地區(qū)為例，分析了混合氣候極值風速設計值，并對該方法中涉及到的極值風速的風向性和相關性等問題進行討論.

1 風速風向理論與極值估計

1.1 極值風速標準值定義[19]

結合上述定義及不同風向極值風速的相關性可得

(1)

在參數(shù)γ的值及Vd,max的分布的已知情況下，根據(jù)式(1)就可以確定各個風向的極值風速vR,d.若每個風向都能獲得一組風速極值數(shù)據(jù)且數(shù)據(jù)量足夠大，Vd,max的經(jīng)驗分布很容易給出，本文將采用Monte Carlo方法來獲得大量風速極值樣本.關鍵在于如何估計參數(shù)γ的值.

1.2 相關性參數(shù)估計方法

參數(shù)γ是用于描述所有風向同時發(fā)生指定重現(xiàn)期極值風速的可能性，顯然不同風向的指定重現(xiàn)期極值風速不一定相同，因此這個數(shù)其實與指定重現(xiàn)期極值風速的大小沒有直接關系，而只與極值風速值在各自風向內的概率分位數(shù)有關.因此，考慮可以用概率變換的方法從概率的角度計算參數(shù)γ值更加方便.

這里采用一種簡單的概率變換，由于年極值風速值與不超越概率值可以一一對應，因此可將不超越概率值Pd視為一個隨機變量，并替換各風向的年極值風速Vd,max.因而事件Vd,max

Pr{Vd,max

Pr{Pd

(2)

由于每年全風向失效概率為1/R，可知，

(3)

觀察上式，可以發(fā)現(xiàn)該式中只有一個未知數(shù)p，非常容易求解.實際應用中，可以采用觀測比例代替概率值.例如，若每個風向都能獲得一組風速極值數(shù)據(jù)且數(shù)據(jù)量足夠大，可以采用如下方法處理.

記vd,i為編號為d的風向內第i年的極值風速樣本，其中d=1,2,…,M，i=1,2,…,N，N為年極值風速樣本數(shù)量.對每個風向的極值風速樣本vd,i做概率變換成pd,i.對于某個指定的風向d，pd,i的取值可表達如下：

(4)

式中：md,i為風向d內極值樣本的排序編號.若vd,i為風向d中的最小值，則md,i為1；若vd,i為風向d中的第2小的值，則md,i為2；以此類推.然后，對第i年的M個不超越概率值取最大值，即

(5)

再將序列{pi,i=1,2,…,N}進行從小到大排序，取第round((N+1)(1-1/R))個值就是p=(1-1/R)1/Mγ的估計值，其中round(·)表示四舍五入取整.經(jīng)過簡單的代數(shù)運算即可求出γ的值.

2 不同氣候下極值風速的Monte Carlo模擬

由于臺風與良態(tài)風的極值風速概率特性完全不同，若仍然對混合氣候風速取樣并擬合，通常會得到不真實的結果.另外，相比良態(tài)風而言，臺風的數(shù)據(jù)記錄十分稀少，不同風向之間的相關性信息極難獲取，因此本文試圖采用Monte Carlo的方法，分別模擬臺風和良態(tài)風的極值風速后將兩者混合，再根據(jù)本文風速風向模型估算混合氣候極值風速.本節(jié)將以上海地區(qū)的混合氣候極值風速估計為例，介紹該方法.

2.1 臺風風速極值模擬

由氣象學的知識知道，雖然臺風的真實結構非常復雜，但是所有臺風的基本結構是非常相似的，臺風移動過程中可以用一些關鍵參數(shù)及移動路徑來描述，可以基于這些隨機參數(shù)進行Monte Carlo抽樣，得到給定位置的臺風風速序列.

肖玉鳳[20]總結了上海地區(qū)的臺風關鍵參數(shù)的統(tǒng)計模型.其數(shù)據(jù)來源為中國氣象局(CMA)提供的《臺風年鑒》(1989年以后更名為《熱帶氣旋年鑒》)(1949—2008)以及國家海洋局(SOA)第一海洋研究所.這些數(shù)據(jù)包括臺風間隔6 h的位置和強度信息.這里涉及的具體關鍵參數(shù)包括：臺風年發(fā)生率、臺風移動方向、臺風移動速度、中心壓差、最小距離、最大風速半徑、衰減模型、風壓徑向分布參數(shù).

在確定了所有臺風關鍵參數(shù)及其統(tǒng)計模型以后，就可基于Monte Carlo思想，采用Yan Meng[17-18]臺風模型模擬目標地區(qū)年臺風發(fā)生個數(shù)及每個臺風在目標地區(qū)的風速大小時程.依據(jù)《建筑結構荷載規(guī)范》GB50009—2012的相關規(guī)定，本文以B類地貌10m高度處的平均風速為研究對象.模擬一年內臺風極值風速，具體過程如下：

(1)根據(jù)臺風年發(fā)生率的參數(shù)分布抽取該年將會發(fā)生的臺風個數(shù)，并為每個臺風抽取其關鍵參數(shù).

(2)將每個臺風按直線路徑經(jīng)過目標地點周圍500 km的模擬圓，每隔1 h計算目標地點的風速，記錄下風速大小和風向.

(3)取一年內每個風向的最大值作為臺風在該年的風速極值樣本，記錄下這些風速值.例如本文中將風向劃分為16個，則此處需要記下16個風速值，若某風向無風，則將風速記為零.

2.2 良態(tài)風風速極值模擬

為了與臺風模擬進行數(shù)據(jù)組合，本文對良態(tài)風的年極值風速也將以Monte Carlo方法進行模擬.

本文所用的上海地區(qū)的風速風向樣本來自上海龍華氣象站1956-1-1—1990-12-31全部12 784 d的日最大風速值(10 min平均時距)和相應風向記錄，由于當?shù)亟?jīng)濟發(fā)展，1990年以后氣象站周圍環(huán)境有所變化，我們沒有采用了該站1990年以后的風速數(shù)據(jù).

本文采用文獻[19]中的考慮風向的獨立風暴法進行統(tǒng)計處理.該方法與Cook[6]的獨立風暴法略有不同，如圖1所示，需要將一個獨立風暴內每個風向的最大值作為各個風向的極值樣本.由于不會在一個獨立風暴內取到同一風向的其他風速值，所以對于同一風向而言，取到的極值樣本是獨立的.另外，這樣做的好處是，既提高了每個風向的樣本數(shù)量，又能在樣本中包含不同風向間極值風速的相關性信息，同時還避免了無法取到獨立風暴內其他風向的風速較大值的問題，圖1中NE、ESE、NNE等均表示風向.

圖1 本文方法與Cook獨立風暴法比較Fig.1 Comparison between the proposed MISand Cook’s MIS

本文取閾值為6.7 m·s-1，采用XIMIS方法[22]擬合各個風向年極值的分布函數(shù).這里需要注意的是，在得到獨立風暴樣本后，須根據(jù)臺風歷史數(shù)據(jù)，去除上海地區(qū)受臺風影響時的獨立風暴樣本，然后再進行統(tǒng)計.

根據(jù)各個風向良態(tài)風年極值的分布信息，可以對良態(tài)風做Monte Carlo模擬，即每個風向按對應的年極值分布隨機抽取一個樣本，記下16個風向的極值風速年極值.

2.3 混合氣候極值樣本獲取

對于每個風向，取上述Monte Carlo模擬得到良態(tài)風年極值和臺風年極值中較大者為混合氣候年極值，并記錄下這些值.重復上述過程1×106次，我們就得到1×106年各個風向混合氣候年極值樣本，可保存成一個1×106行16列的矩陣.

在Monte Carlo模擬中需要說明以下兩點：

(1)良態(tài)風Monte Carlo模擬中不考慮風向之間的相關性，即默認良態(tài)風各個風向極值風速不相關，相關度參數(shù)γ值為1.文獻[19]證實了良態(tài)風各個風向的風速是漸進獨立的，這意味著在討論極值問題時，不同風向的極值風速可以假設為完全獨立.同時還證明了這一假設是偏于保守的，并且極值分析結果對微弱的相關性不敏感.

(2)良態(tài)風Monte Carlo模擬中去除了受臺風影響的日最大風速數(shù)據(jù)，這些臺風風速數(shù)據(jù)在考慮風向相關性的極值估計中是沒有用的.這是因為一個強臺風在1 d內就可以使多個風向的風速達到極值的水平，用日最大風速來研究極值風向性是不合適的.但是，這些日最大風速可以用來檢驗臺風模擬結果是否與實際相符.

2.4 混合氣候極值樣本合理性檢驗

可以將混合氣候極值樣本與實測數(shù)據(jù)進行對比.根據(jù)上海龍華氣象站所有日最大風速，取8 d最大風速或獨立風暴極值風速[6](不考慮風向)，將這組數(shù)據(jù)與混合氣候年極值風速(同一年內取所有風向中風速最大值)一同繪制在Gumbel概率圖中，如圖2和圖3所示，圖中橫坐標中的P為超越概率值.這里橫坐標值采用Cook[16]描述的方法來計算.從圖中看出，兩者的經(jīng)驗分布相差不大.

除了從圖中直接觀察外，還可以用假設檢驗的方法來檢驗樣本合理性.預先假設混合氣候模擬結果為真實的極值風速分布，并設定的檢驗水準為α=0.05.令v1

圖2 模擬結果與實測數(shù)據(jù)對比(8 d極值)

Fig.2Comparisonbetweentheempiricaldistributionsofthesimulatedresultsandthemeasureddata(eight-daymaxima)

圖3 模擬結果與實測數(shù)據(jù)對比(獨立風暴極值)

Fig.3Comparisonbetweentheempiricaldistributionsofthesimulatedresultsandthemeasureddata(stormmaxima)

(6)

其中Γ(·)為Gamma函數(shù).pm的95%置信區(qū)間上下限值Pm,u和Pm,l分別由下列兩式確定

(7)

(8)

式中，Ix(z,w)稱為不完整的Beta函數(shù)(incomplete beta function)，其表達式為

(9)

在Matlab軟件中定義了Ix(z,w)的反函數(shù)(函數(shù)名為“betaincinv”)，所以Pm,u和Pm,l很容易用Matlab編程求出，其對應的一年內的不超越概率為(Pm,u)r和(Pm,l)r，其中r為實測風速樣本年平均個數(shù).再結合混合氣候年極值的經(jīng)驗分布(由模擬結果得到)，可求出vm的95%置信區(qū)間上下限值vm,u和vm,l，畫在圖2和圖3中.從圖中看出，實測數(shù)據(jù)均不在拒絕域內，因此混合氣候模擬結果是可信的.

3 考慮風向的混合氣候極值風速估計方法

通過Monte Carlo模擬獲得混合氣候各個風向年極值樣本，該樣本數(shù)量巨大，因此可以很容易地給出每個風向年極值風速的經(jīng)驗分布.另一方面，采用上文方法估算相關度參數(shù)γ值.確定上述兩者后，就可以應用式(1)確定各個風向的重現(xiàn)期為R年的極值風速vR,d.實際計算時，可簡單地將每個風向的年極值風速樣本從小到大排序，取第round((N+1)p)個值就是此風向極值風速vR,d值，其中p=(1-1/R)1/Mγ為單風向不超越概率，N為模擬年份數(shù)量，本文中N=1×106.

經(jīng)計算，上海地區(qū)重現(xiàn)期為50年時，全風向極值風速標準值為26.74 m·s-1，極值風速標準值的最大值發(fā)生在東偏南(ESE)方向，大小為29.57 m·s-1，最小值發(fā)生在西偏北(WNW)方向，大小為22.87 m·s-1.圖4為上海地區(qū)重現(xiàn)期為50年的極值風速標準值玫瑰圖，圖中N、NNE、NE等均表示風向.該結果十分符合上海地區(qū)東南沿海的地理特性.

圖4 上海地區(qū)重現(xiàn)期為50年的極值風速玫瑰圖

Fig.4Characteristicvaluesofwindspeedunderthereturnperiodof50yearsinShanghai

4 討論

4.1 臺風氣候與良態(tài)風氣候分開統(tǒng)計的優(yōu)勢

對于受臺風影響較為顯著的區(qū)域而言，不區(qū)分臺風氣候和良態(tài)風氣候進行混合取樣將導致不合理的結果，這是因為臺風與良態(tài)風風速極值分布規(guī)律通常是有很大區(qū)別的.圖5和圖6以階段極值法為例，給出了上海地區(qū)全風向本文方法(臺風氣候與良態(tài)風氣候分開統(tǒng)計)與傳統(tǒng)階段極值擬合方法(直接統(tǒng)計混合氣候數(shù)據(jù)，不區(qū)分臺風氣候與良態(tài)風氣候)結果比較.從圖中可以看出，基于8 d極值風速樣本的結果在上尾段明顯上揚，傳統(tǒng)的階段極值法很難處理這一現(xiàn)象，并且造成了觀察時距越大極值風速估計值越大的錯誤結論.這是因為上海極值風速上尾部分由臺風控制，而下尾部分由良態(tài)風控制，且良態(tài)風占大多數(shù)，延長觀察時距實際上就是減小下尾部分的樣本數(shù)量，也就是減小良態(tài)風數(shù)據(jù)樣本在極值風速數(shù)據(jù)總樣本中的權重，因而出現(xiàn)了觀察時距越大極值風速估計值越大.本文方法由于臺風和良態(tài)風分別處理，其結果較好地展現(xiàn)了兩種氣候的極值分布形態(tài)，且有合理的過渡，這種方法更加符合實際情況.

圖5 本文方法與不同觀察時距下的階段極法結果比較(年極值分布對比)

Fig.5Comparisonbetweentheproposedmethodandthetraditionalmethod(distributionsontheGumbelprobabilitypaper)

現(xiàn)有的各類極值方法往往以一個單一的分布來描述極值，因而在不同氣候混合而成的階段性風速極值分布時可能出現(xiàn)困難.而Monte Carlo方法具有靈活的優(yōu)點，其不依賴特定的極值分布.此外，極值方法的誤差來源于函數(shù)擬合誤差或分析誤差以及樣本誤差，不同的方法性質不同因而各有特色及適用范圍.Monte Carlo方法特色在于，只要模擬數(shù)量足夠多，則可將分析誤差無限縮小，在條件允許的情況下可以做到結果誤差僅來源于樣本誤差.也就是說，本文中只要模擬樣本數(shù)量足夠多，則極值估計結果精度僅取決于輸入數(shù)據(jù)的質量，即根據(jù)現(xiàn)有數(shù)據(jù)所得到的良態(tài)風和臺風模擬樣本的質量.

圖6 本文方法與不同觀察時距下的階段極法結果比較(50年重現(xiàn)期基本風速對比)

Fig.6Comparisonbetweentheproposedmethodandthetraditionalmethod(characteristicvaluesofwindspeedunderthereturnperiodof50years)

4.2 各個風向之間風速極值相關度參數(shù)

從求解γ值的過程中可以看出，該值是會隨著重現(xiàn)期R的大小而變化的.圖7展示了當單個風向的樣本容量N=1×106時，計算得到的各種風氣候下γ值隨重現(xiàn)期R的變化情況.從圖中可看出，對于較高的極值風速，混合氣候的γ值主要由臺風部分控制，并且顯示出較強的相關性.

圖7 各種氣候下γ值隨重現(xiàn)期R的變化Fig.7 Change in the value of γ with the returnperiod R

需要指出的是，本次模擬中良態(tài)風各個風向的年極值風速始終默認為互相獨立，即良態(tài)風的γ應始終為1.圖中當R>1 000時，良態(tài)風的γ產(chǎn)生波動而不再保持為1.即當N/R<1 000時，由于極值數(shù)量有限γ值將不再準確，進而使得極值風速產(chǎn)生誤差.因此若要計算重現(xiàn)期為R年的各個風向的極值風速，建議將單個風向的樣本容量N設定為1 000R以上.

5 結語

對于目前尚沒有考慮風向的混合氣候極值估計方法的現(xiàn)狀，本文試圖采用Monte Carlo方法分別獲取臺風和良態(tài)風氣候的風速極值樣本，并結合風速風向分布模型，給出一種新的考慮風向的極值風速估計方法.該方法能較好地展現(xiàn)臺風和良態(tài)風兩種氣候的極值分布形態(tài)，同時有效地結合了兩者各自的極值風速特性，給出的結果較傳統(tǒng)方法更為合理可信.

致謝

感謝中國氣象科學數(shù)據(jù)共享服務網(wǎng)(http://cdc.cma.gov.cn)提供氣象數(shù)據(jù)支持.

參考文獻:

[1] DAVENPORT A G. The application of statistical concepts to the wind loading of structures [J]. Proceedings of the Institution of Civil Engineers, 1961, 19(4): 449.

[2] DAVENPORT A G. Note on the distribution of the largest value of a random function with application to gust loading [J]. Proceedings of the Institution of Civil Engineers, 1964, 28(2): 187.

硝鹽是國內外肉制品中廣為使用的食品添加劑，可同時實現(xiàn)發(fā)色、抑菌、防腐、增香等作用，但過量或不當添加使用，將存在中毒甚至致癌隱患，特別是對于小作坊式加工，因缺乏應有的專業(yè)知識和必要的監(jiān)控手段，安全性難以確保。

[3] DAVENPORT A G. The buffeting of structures by gusts. Wind Effects on Buildings and Structures [C]//Proceedings of the Conference Held at the National Physical Laboratory. London, HMSO, 1965: 357-391.

[4] DAVENPORT A G. The relationship of reliability to wind load [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1983, 13(1-3): 3.

[5] DAVENPORT A G，SURRY D，STATHOPOULOS T. Wind loads on low rise buildings: final report of phases I and II BLWT-SS8-1977[R]. London: University of Western Ontario, 1977.

[6] COOK N J. Towards better estimation of extreme winds [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1982, 9(3): 295.

[7] COOK N J. Note on directional and seasonal assessment of extreme winds for design [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1983, 12(3): 365.

[8] COOK N J，MILLER C A. Further note on directional assessment of extreme winds for design [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1999, 79(3): 201.

[9] 全涌，劉磊，顧明. 考慮風向的極值風速估計的Cook方法改進 [J]. 同濟大學學報(自然科學版), 2015, 43(2): 189.

QUAN Yong, LIU Lei, GU Ming. Improvement of Cook method considering directional extreme wind speed [J]. Journal of Tongji University(Natural Science), 2015, 43(2): 189.

[10] HARAGUCHI K，KANDA J. Probability model for annual maximum wind speeds in multi-direction [C]//Wind Engineering into 21st Century(ICWE). Copenhagen: [s.n.], 1999: 205-212.

[11] KANDA J， ITOI T. Correlated Gumbel probability model for directional wind speeds[C]// Structural Safety and Reliability: Proceedings of the Eighth International Conference, ICOSSAR'01. Newport Beach: CRC Press, 2001:8-12.

[12] ZHANG X，CHEN X. Assessing probabilistic wind load effects via a multivariate extreme wind speed model: A unified framework to consider directionality and uncertainty [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2015, 147: 30.

[13] ZHANG X，CHEN X. Influence of dependence of directional extreme wind speeds on wind load effects with various mean recurrence intervals [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2016, 148: 45.

[14] GOMES L，VICKERY B J. Extreme wind speeds in mixed wind climates [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1978, 2(4): 331.

[15] COOK N J，HARRIS R I，WHITING R. Extreme wind speeds in mixed climates revisited [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2003, 91(3): 403.

[16] COOK N J. Confidence limits for extreme wind speeds in mixed climates [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2004, 92(1): 41.

[17] MENG Y，MATSUI M，HIBI K. An analytical model for simulation of the wind field in a typhoon boundary layer [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1995, 56(2-3): 291.

[18] MENG Y，MATSUI M，HIBI K. A numerical study of the wind field in a typhoon boundary layer [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1997, 67-68: 437.

[19] 張秉超. 建筑圍護結構極值風壓的全概率分析方法研究 [D]. 上海: 同濟大學, 2014.

ZHANG Bingchao. Research on the fully probabilistic method of the extreme wind pressure on buildings [D]. Shanghai: Tongji University, 2014.

[20] 肖玉鳳. 基于數(shù)值模擬的東南沿海臺風危險性分析及輕鋼結構風災易損性研究 [D]. 哈爾濱: 哈爾濱工業(yè)大學, 2011.

XIAO Yufeng. Typhoon wind hazard analysis based on numerical simulation and fragility of light-gauge steel structure in southeast china coastal regions [D]. Harbin: Harbin Institute of Technology, 2011.

[21] 中華人民共和國國家標準.建筑結構荷載規(guī)范：GB50009-2012 [S]. 北京: 中國建筑工業(yè)出版社, 2012.

Professional Standard of the People’s Republic of China. Load code for the design of building structures: GB50009-2012 [S]. Beijing: China Architecture and Building Press, 2012.

[22] HARRIS R I. XIMIS, a penultimate extreme value method suitable for all types of wind climate [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2009, 97(5/6): 271.

[23] HARRIS R I. Gumbel re-visited-a new look at extreme value statistics applied to wind speeds [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1996, 59(1): 1.