基于Copula函數的風浪多方向極限狀態曲線

涂志斌， 黃銘楓， 樓文娟， 李 蓓

(1. 浙江水利水電學院 建筑工程學院，杭州 310018； 2. 浙江大學 建筑工程學院結構工程研究所，杭州 310058)

ECM(environmental contour method)是多維隨機環境變量聯合作用下結構極限荷載效應估計的常用方法，其核心思想是在具有指定超越概率的環境變量極限狀態曲線上搜尋荷載效應極值作為結構的極限荷載效應[1]，在海洋工程[2-3]、地震工程[4]和風工程[5]中均有廣泛的應用。其中構造具有指定超越概率的多維隨機環境變量極限狀態曲線是ECM的基礎。

Rosenblatt變換是極限狀態曲線構造的一般方法。該方法首先在標準正態空間中構造一個與可靠指標對應的圓，再將該圓映射到物理空間中得到具有指定超越概率的隨機環境變量極限狀態曲線[6]。根據該變換，Agarwal等[6]和Valamanesh等[7]分別構造了某海域平均風速和有效波高的極限狀態曲線；Saranyasoontorn等和Haver等[8]分別構造了某海域有效波高和譜峰周期的極限狀態曲線。在以上研究中隨機變量的聯合分布由主變量的邊緣分布和其他變量的條件分布表達，采用Rosenblatt變換構造極限狀態曲線較為簡便。若隨機變量的聯合分布由傳統聯合分布模型(如皮爾遜Ⅲ型分布[9]、二維Gumbel邏輯分布[10]、Nataf分布[11-12]等)或Copula函數[13-15]表達，在Rosenblatt變換之前則需先根據聯合分布求解條件分布。然而由于以上聯合分布模型采用了較為復雜的相關結構，條件分布的求解較為困難。Huseby等[16-17]基于蒙特卡洛模擬(Monte Carlo simulation，MCS)提出通過求解隨機樣本具有指定超越概率的凸集的超平面來構造極限狀態曲線，直接避免了Rosenblatt變換需要求解隨機變量條件分布這一難題，提高了極限狀態曲線構造的可行性。然而針對該方法有三點需要說明：①隨機樣本的生成依賴于隨機變量的聯合分布模型，若聯合分布模型復雜，則隨機樣本的生成算法也復雜；②求解精度和效率與樣本容量密切相關，若樣本容量較小，極限狀態曲線在局部將出現封閉環；③該方法將凸集邊界近似為所有超平面構成的封閉曲線，使所求極限狀態曲線略大于真實極限狀態曲線。Vanem等[18]詳細比較了兩種構造極限狀態曲線的方法，提出了各自的適用范圍。

風向和波向是海洋環境中重要的隨機變量，對海洋結構設計有重要影響[19-20]。Zhang等[21]研究發現各方向的平均風速和有效波高可能服從不同的邊緣分布模型，二者的相關關系也不盡相同。季新然等[22]認為波向分布對波浪爬高具有較大的影響，在實際設計中若不考慮波向的影響可能低估或高估群墩周圍的波浪爬高。吳家鳴等[23]在研究波向對導管架平臺主樁應力響應的影響時發現，存在一個最不利波向使主樁應力最大。同時大量風向與建筑物抗風可靠性的研究表明若不考慮風向效應，極值風速和結構設計風荷載均有可能被高估[24-27]。然而目前與風向和波向對風浪極限狀態曲線及海洋結構極限荷載效應影響相關的報道還不多見。

本文在Rosenblatt變換的基礎上引入Copula函數，提出了二維隨機變量極限狀態曲線構造的改進方法；結合風浪同步觀測數據研究了風向與波向對風浪極限狀態曲線的影響；估計了某導管架平臺各方向的基底剪力極限荷載效應及對應的平均風速和有效波高取值，并得到了幾點有益的結論可供海洋結構設計參考。

1 基于Copula函數的極限狀態曲線

1.1 Rosenblatt變換

設二維隨機變量E=(E1,E2)服從[0,1]2內的獨立均勻分布；Z=(Z1,Z2)為標準正態空間中相互獨立的二維隨機變量；ei=Φ(zi)(i=1,2)，Φ(·)為一維標準正態分布函數；X=(X1,X2)為物理空間中具有任意相關關系的二維隨機變量。根據Rosenblatt變換有

(1)

式中：F1(x1)為變量X1的邊緣分布函數；F2|1(x2|x1)為變量X2的條件分布函數，由X2的邊緣分布和X1與X2的相關結構決定。對式(1)求逆可得

(2)

在標準正態空間中極限狀態曲線是半徑等于可靠指標β的圓，即

(3)

1.2 Copula函數的基本理論

Copula函數是一種構造聯合分布模型的工具系統，可構造具有任意相關關系和邊緣分布的多維隨機變量的聯合分布函數，并將現有的多種聯合分布模型納入其中，如Gumbel分布、Nataf分布等。設二維隨機變量U=(U1,U2)服從[0,1]2內的獨立均勻分布，且ui=Fi(xi)(i=1,2)，那么根據Copula函數，二維隨機變量X的聯合分布函數F(x1,x2)可表達為

F(x1,x2)=C[F1(x1),F1(x2)]=C(u1,u2)

(4)

式中：C(·)為Copula函數；F2(x2)為變量X2的邊緣分布函數。若F1(x1),F2(x2)均連續，則C(·)唯一確定。由式(4)可知，當隨機變量的邊緣分布已知時，Copula函數的作用在于確定變量之間的相關結構。

表1 常用二維Copula函數族及其典型函數

采用Copula函數構造隨機變量的聯合分布模型時，同時滿足以下兩個條件的Copula函數稱為最優Copula函數：①Ui=Fi(xi)能準確地描述Xi的邊緣分布特性；②由Copula函數確定的相關結構能準確描述各隨機變量間的相關關系。目前常用的最優Copula函數評價準則為AIC(Akaike information criterion)，其表達式為[28-29]

(5)

1.3 改進方法

根據Copula函數理論的相關研究，條件分布函數F2|1(x2|x1)與Copula函數的關系為[30]

(6)

式中，C2|1(u2|u1)為條件Copula函數。將式(6)代入式(1)中的第2式可得

e2=C2 |1(u2|u1)

(7)

那么，

(8)

2 風浪實測數據與樣本選取

本文以紐約東海岸附近63115號海洋觀測站(北緯40.50°，西經-72.92°)的風參數和波浪參數觀測記錄為統計樣本，數據來源于美國陸軍工程兵團工程研究發展中心官網。該統計樣本的記錄時間為1980年1月1日—2014年12月31日，但部分時間段的統計樣本有缺失。各參數均同步記錄，每小時記錄一次。同步提取平均風速日極值v和對應的有效波高Hs、風向φv及波向φH作為風浪多方向極限狀態曲線分析的樣本。風向與波向之差(φv-φH)的角直方圖，如圖1所示。其中風向和波向均為記錄時間處的瞬時值。圖1表明風向和波向的差異較小。根據規范IEC 61400-3和Zhang等的建議，當風向和波向差異較小時，在風浪聯合作用分析中可認為風向和波向保持一致。將風向歸類到8個角度區間，其玫瑰圖如圖2所示。由圖2可知：日極值風速的風向記錄在SE，S，SW和W方向上的頻率較高，分別為16.57%，17.42%，13.37%和18.17%；而在NW，N，NE和E方向上的頻率較低，分別為10.94%，9.52%，7.59%和6.40%。由于波向與風向的偏差較小，圖2也表達了波向在各方向上的分布頻率。SE(深色)、W(深色)及全方向(淺色)3個方向上的日極值風速及對應有效波高散點圖，如圖3所示。由圖3可知：①不同方向上散點的尾部特性有所不同，在W方向上尾部較為厚重，分布范圍為[15 m/s, 20 m/s]，在SE方向上尾部輕且較為離散，分布范圍為[10 m/s, 20 m/s]；②在3個方向上平均風速和有效波高的線性相關系數分別為ρSE=0.770 1，ρW=0.826 3和ρAll=0.744 5，均具有較強的相關性，但相關程度又有所不同。這說明不同方向上平均風速和有效波高的邊緣分布特性和相關關系均有所不同。

圖1 風向與波向之差角直方圖

圖2 日極值風速對應風向玫瑰圖

圖3 日極值風速及對應有效波高散點圖

3 風浪多方向極限狀態曲線

3.1 風浪邊緣分布

目前平均風速和有效波高極值的邊緣分布多采用Weibull分布、Gamma分布和Lognormal分布來描述，其表達式分別為：

Weibull

(9)

Gamma

(10)

Lognormal

(11)

表2 SE方向上日極值風速和有效波高邊緣分布估計結果

圖4 各風向下日極值風速的邊緣分布擬合

圖5 各風向下同步有效波高的邊緣分布擬合

3.2 風浪聯合分布

從表1中選擇最優Copula函數來構造日極值風速和同步有效波高的聯合分布模型。采用偽極大似然估計法估計各Copula函數的參數[31]，并根據AIC準則選擇最優Copula函數，結果如表3所示，其中括號外的數據為參數估計值，括號內的數據為AIC值。由表3可知，除NW方向的最優Copula函數是Plackett Copula外，其他方向的最優Copula函數均為Gumbel Copula。

表3 日極值風速和同步有效波高聯合分布模型參數估計值及AIC值

SE，W和All 3個方向的日極值風速和同步有效波高聯合分布，如圖6所示。由圖6可知，3個方向的峰形均呈現出非對稱性但并不相同，具體表現為：①峰值的取值坐標不同，SE，W和All方向峰值對應的日極值風速和同步有效波高分別為(6.5 m/s，1 m)，(11.5 m/s，1.6 m)，(8.7 m/s，1.2 m)；②與SE和All方向相比，W方向的峰形較陡且狹長。以上不同是因為各方向的日極值風速和同步有效波高的邊緣分布及相關性均不相同，需要采用不同的模型來描述其聯合分布。

圖6 日極值風速和同步有效波高的聯合分布

3.3 改進方法的驗證

考察數值算法中角度離散點數n對基于Copula函數的二維隨機變量極限狀態曲線的影響。當步長Δθ=3°，Δθ=6°，Δθ=10°和Δθ=15°(對應n=120，n=60，n=36，n=24)時，SE方向日極值風速和同步有效波高的極限狀態曲線(重現期T= 50 a)，如圖7所示。由圖7可知，當n=60時極限狀態曲線趨于光滑。對其他各方向、各重現期下的極限狀態曲線進行考察可得到相同的結論。因此在本文中采用數值算法計算基于Copula函數的日極值風速和同步有效波高的極限狀態曲線時角度離散點數取n=60。

圖7 不同步長下T=50 a時SE方向日極值風速和同步有效波高的極限狀態曲線

本文將對比基于Copula函數的極限狀態曲線和基于MCS的極限狀態曲線，以驗證改進方法的準確性。基于MCS的極限狀態曲線構造方法由Huseby等提出，其基本步驟為：①根據二維隨機變量的聯合分布模型生成大量統計樣本；②在統計樣本的基礎上構造一個超越概率為Pe的凸集，并在[0,2π)內尋找其超平面；③求解相鄰兩超平面的交點并依次連接即可得到極限狀態曲線。按照以上步驟構造重現期T=1 a，T=20 a，T=50 a時各方向上日極值風速和同步有效波高的極限狀態曲線。由于日極值風速對應的統計時長為24 h，重現期T與超越概率Pe的換算公式為Pe=1/365.25/T，T=1 a，T=20 a，T=50 a對應的超越概率分別為Pe=2.7×10-3，Pe=1.37×10-4，Pe=5.48×10-5。根據日極值風速和同步有效波高的最優邊緣分布及最優Copula函數生成統計樣本，超平面及其交點的求解公式為

(12)

(13)

各重現期下SE方向上基于Gumbel Copula函數的日極值風速和同步有效波高極限狀態曲線，見圖8。總體而言，各重現期下基于Copula函數的極限狀態曲線與基于MCS的極限狀態曲線差異較小。在其他方向上作相同的對比可觀察到類似的現象。這表明本文提出的改進方法具有較高的可靠性。同時從整個計算過程來看，改進方法的數值算法具有較高的可操作性，克服了MCS的局限性。

圖8 SE方向日極值風速和同步有效波高的極限狀態曲線

3.4 風浪極限狀態曲線

當重現期T=1 a，T=20 a和T=50 a時，各方向基于Copula函數的日極值風速和同步有效波高的極限狀態曲線，如圖9所示。由圖9可知：①隨著重現期的增加，極限狀態曲線的形狀基本不變，但軌跡向外擴大；②除NW外，其他方向的極限狀態曲線因采用同一類型的邊緣分布和Copula函數構造聯合分布模型而形狀相似，但因邊緣分布和Copula函數的參數不同而軌跡相差較大，特別是在日極值風速和同步有效波高聯合取值較大的區域；③從圖示箭頭方向來看，在各重現期下極限狀態曲線由外向內依次為NE，E，NW，SE，N，All，W，S，SW，其中All方向的極限狀態曲線處于各曲線的中間位置。

圖9 基于Copula函數的日極值風速和同步有效波高的多方向極限狀態曲線

采用ECM估計某導管架平臺基底剪力的極限荷載效應。根據Liu等[32]的研究，某導管架平臺的基底剪力Q、平均風速v和有效波高Hs的關系為

(14)

根據式(14)和圖9所示的風浪極限狀態曲線估計各方向不同重現期下該導管架平臺的v，Hs和Q取值，結果如表4所示。與All方向相比，其他方向v,Hs和Q的取值誤差r，如表5所示。由表4和表5可知：①隨著重現期的增加，各方向v，Hs和Q的取值均有所增加；②在同一重現期下各方向v，Hs和Q的取值各不相同，其中NE方向的取值最大(最不利方向)，All方向的取值居中，NW方向的取值最小，這與各方向風浪極限狀態曲線的軌跡分布相符；③與All方向相比，NE，E，SE和W 4個方向v，Hs和Q的取值偏大，其中在NE方向，當T=50 a時三者的誤差分別為8.98%，60.00%和141.18%；④與All方向相比，N，S，SW和NW 4個方向v，Hs和Q的取值偏小，其中在NW方向，當T=50 a時三者的誤差分別為-8.59%，-78.82%和-82.35%。盡管根據不同的風浪觀測數據，最不利方向和與All方向相比各方向上v，Hs和Q的取值偏差會有所不同，但若某方向上的極限狀態曲線大于All方向，那么在該方向上的v，Hs和Q的取值將偏大，反之亦然。由此可知，對導管架平臺進行結構設計時，不對風、浪的方向加以區分而直接采用全方向的統計數據構造極限狀態曲線并估計結構極限荷載效應，會造成不同程度的高估或低估。出現高估或低估的原因是在主風(波)向上風浪較大而在其他方向上風浪較小，直接構造All方向的極限狀態曲線忽略了這一特性，得到了一個較為平均的極限狀態曲線。而分方向進行極限狀態曲線的構造和結構極限荷載效應的估計則考慮了風浪在各方向上的統計規律，可得到更加精確的結果。

表4 某導管架平臺基底剪力極限荷載效應及對應的平均風速和有效波高取值

表5 某導管架平臺各方向基底剪力極限荷載效應及對應平均風速和有效波高取值的誤差

4 結 論

本文在Rosenblatt變換的基礎上引入Copula函數，提出了二維隨機環境變量極限狀態曲線構造的改進方法及其數值算法；結合風浪同步觀測數據研究了風浪多方向極限狀態曲線，并在此基礎上估計了某導管架平臺的基底剪力極限荷載效應及對應的平均風速和有效波高取值，得到了以下結論：

(1) 基于Copula函數的極限狀態曲線構造改進方法及其數值算法將隨機環境變量的條件分布函數及其逆函數的求解轉化為條件Copula函數及其逆函數的求解，具有較高的可靠性和可操作性。

(2) 各方向的日極值風速和同步有效波高極限狀態曲線軌跡呈現出較大的差異，而全方向的極限狀態曲線處于各條曲線的中間位置。

(3) 各方向導管架平臺基底剪力極限荷載效應及對應的平均風速和有效波高取值差異較大。與全方向相比，極限狀態曲線較大方向的平均風速、有效波高和基底剪力極限荷載效應取值也較大。