基于Copula函數(shù)的風(fēng)浪多方向極限狀態(tài)曲線

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

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

ECM(environmental contour method)是多維隨機(jī)環(huán)境變量聯(lián)合作用下結(jié)構(gòu)極限荷載效應(yīng)估計(jì)的常用方法，其核心思想是在具有指定超越概率的環(huán)境變量極限狀態(tài)曲線上搜尋荷載效應(yīng)極值作為結(jié)構(gòu)的極限荷載效應(yīng)[1]，在海洋工程[2-3]、地震工程[4]和風(fēng)工程[5]中均有廣泛的應(yīng)用。其中構(gòu)造具有指定超越概率的多維隨機(jī)環(huán)境變量極限狀態(tài)曲線是ECM的基礎(chǔ)。

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

風(fēng)向和波向是海洋環(huán)境中重要的隨機(jī)變量，對海洋結(jié)構(gòu)設(shè)計(jì)有重要影響[19-20]。Zhang等[21]研究發(fā)現(xiàn)各方向的平均風(fēng)速和有效波高可能服從不同的邊緣分布模型，二者的相關(guān)關(guān)系也不盡相同。季新然等[22]認(rèn)為波向分布對波浪爬高具有較大的影響，在實(shí)際設(shè)計(jì)中若不考慮波向的影響可能低估或高估群墩周圍的波浪爬高。吳家鳴等[23]在研究波向?qū)?dǎo)管架平臺(tái)主樁應(yīng)力響應(yīng)的影響時(shí)發(fā)現(xiàn)，存在一個(gè)最不利波向使主樁應(yīng)力最大。同時(shí)大量風(fēng)向與建筑物抗風(fēng)可靠性的研究表明若不考慮風(fēng)向效應(yīng)，極值風(fēng)速和結(jié)構(gòu)設(shè)計(jì)風(fēng)荷載均有可能被高估[24-27]。然而目前與風(fēng)向和波向?qū)︼L(fēng)浪極限狀態(tài)曲線及海洋結(jié)構(gòu)極限荷載效應(yīng)影響相關(guān)的報(bào)道還不多見。

本文在Rosenblatt變換的基礎(chǔ)上引入Copula函數(shù)，提出了二維隨機(jī)變量極限狀態(tài)曲線構(gòu)造的改進(jìn)方法；結(jié)合風(fēng)浪同步觀測數(shù)據(jù)研究了風(fēng)向與波向?qū)︼L(fēng)浪極限狀態(tài)曲線的影響；估計(jì)了某導(dǎo)管架平臺(tái)各方向的基底剪力極限荷載效應(yīng)及對應(yīng)的平均風(fēng)速和有效波高取值，并得到了幾點(diǎn)有益的結(jié)論可供海洋結(jié)構(gòu)設(shè)計(jì)參考。

1 基于Copula函數(shù)的極限狀態(tài)曲線

1.1 Rosenblatt變換

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

(1)

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

(2)

在標(biāo)準(zhǔn)正態(tài)空間中極限狀態(tài)曲線是半徑等于可靠指標(biāo)β的圓，即

(3)

1.2 Copula函數(shù)的基本理論

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

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

(4)

式中：C(·)為Copula函數(shù)；F2(x2)為變量X2的邊緣分布函數(shù)。若F1(x1),F2(x2)均連續(xù)，則C(·)唯一確定。由式(4)可知，當(dāng)隨機(jī)變量的邊緣分布已知時(shí)，Copula函數(shù)的作用在于確定變量之間的相關(guān)結(jié)構(gòu)。

表1 常用二維Copula函數(shù)族及其典型函數(shù)

采用Copula函數(shù)構(gòu)造隨機(jī)變量的聯(lián)合分布模型時(shí)，同時(shí)滿足以下兩個(gè)條件的Copula函數(shù)稱為最優(yōu)Copula函數(shù)：①Ui=Fi(xi)能準(zhǔn)確地描述Xi的邊緣分布特性；②由Copula函數(shù)確定的相關(guān)結(jié)構(gòu)能準(zhǔn)確描述各隨機(jī)變量間的相關(guān)關(guān)系。目前常用的最優(yōu)Copula函數(shù)評價(jià)準(zhǔn)則為AIC(Akaike information criterion)，其表達(dá)式為[28-29]

(5)

1.3 改進(jìn)方法

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

(6)

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

e2=C2 |1(u2|u1)

(7)

那么，

(8)

2 風(fēng)浪實(shí)測數(shù)據(jù)與樣本選取

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

圖1 風(fēng)向與波向之差角直方圖

圖2 日極值風(fēng)速對應(yīng)風(fēng)向玫瑰圖

圖3 日極值風(fēng)速及對應(yīng)有效波高散點(diǎn)圖

3 風(fēng)浪多方向極限狀態(tài)曲線

3.1 風(fēng)浪邊緣分布

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

Weibull

(9)

Gamma

(10)

Lognormal

(11)

表2 SE方向上日極值風(fēng)速和有效波高邊緣分布估計(jì)結(jié)果

圖4 各風(fēng)向下日極值風(fēng)速的邊緣分布擬合

圖5 各風(fēng)向下同步有效波高的邊緣分布擬合

3.2 風(fēng)浪聯(lián)合分布

從表1中選擇最優(yōu)Copula函數(shù)來構(gòu)造日極值風(fēng)速和同步有效波高的聯(lián)合分布模型。采用偽極大似然估計(jì)法估計(jì)各Copula函數(shù)的參數(shù)[31]，并根據(jù)AIC準(zhǔn)則選擇最優(yōu)Copula函數(shù)，結(jié)果如表3所示，其中括號外的數(shù)據(jù)為參數(shù)估計(jì)值，括號內(nèi)的數(shù)據(jù)為AIC值。由表3可知，除NW方向的最優(yōu)Copula函數(shù)是Plackett Copula外，其他方向的最優(yōu)Copula函數(shù)均為Gumbel Copula。

表3 日極值風(fēng)速和同步有效波高聯(lián)合分布模型參數(shù)估計(jì)值及AIC值

SE，W和All 3個(gè)方向的日極值風(fēng)速和同步有效波高聯(lián)合分布，如圖6所示。由圖6可知，3個(gè)方向的峰形均呈現(xiàn)出非對稱性但并不相同，具體表現(xiàn)為：①峰值的取值坐標(biāo)不同，SE，W和All方向峰值對應(yīng)的日極值風(fēng)速和同步有效波高分別為(6.5 m/s，1 m)，(11.5 m/s，1.6 m)，(8.7 m/s，1.2 m)；②與SE和All方向相比，W方向的峰形較陡且狹長。以上不同是因?yàn)楦鞣较虻娜諛O值風(fēng)速和同步有效波高的邊緣分布及相關(guān)性均不相同，需要采用不同的模型來描述其聯(lián)合分布。

圖6 日極值風(fēng)速和同步有效波高的聯(lián)合分布

3.3 改進(jìn)方法的驗(yàn)證

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

圖7 不同步長下T=50 a時(shí)SE方向日極值風(fēng)速和同步有效波高的極限狀態(tài)曲線

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

(12)

(13)

各重現(xiàn)期下SE方向上基于Gumbel Copula函數(shù)的日極值風(fēng)速和同步有效波高極限狀態(tài)曲線，見圖8?？傮w而言，各重現(xiàn)期下基于Copula函數(shù)的極限狀態(tài)曲線與基于MCS的極限狀態(tài)曲線差異較小。在其他方向上作相同的對比可觀察到類似的現(xiàn)象。這表明本文提出的改進(jìn)方法具有較高的可靠性。同時(shí)從整個(gè)計(jì)算過程來看，改進(jìn)方法的數(shù)值算法具有較高的可操作性，克服了MCS的局限性。

圖8 SE方向日極值風(fēng)速和同步有效波高的極限狀態(tài)曲線

3.4 風(fēng)浪極限狀態(tài)曲線

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

圖9 基于Copula函數(shù)的日極值風(fēng)速和同步有效波高的多方向極限狀態(tài)曲線

采用ECM估計(jì)某導(dǎo)管架平臺(tái)基底剪力的極限荷載效應(yīng)。根據(jù)Liu等[32]的研究，某導(dǎo)管架平臺(tái)的基底剪力Q、平均風(fēng)速v和有效波高Hs的關(guān)系為

(14)

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

表4 某導(dǎo)管架平臺(tái)基底剪力極限荷載效應(yīng)及對應(yīng)的平均風(fēng)速和有效波高取值

表5 某導(dǎo)管架平臺(tái)各方向基底剪力極限荷載效應(yīng)及對應(yīng)平均風(fēng)速和有效波高取值的誤差

4 結(jié) 論

本文在Rosenblatt變換的基礎(chǔ)上引入Copula函數(shù)，提出了二維隨機(jī)環(huán)境變量極限狀態(tài)曲線構(gòu)造的改進(jìn)方法及其數(shù)值算法；結(jié)合風(fēng)浪同步觀測數(shù)據(jù)研究了風(fēng)浪多方向極限狀態(tài)曲線，并在此基礎(chǔ)上估計(jì)了某導(dǎo)管架平臺(tái)的基底剪力極限荷載效應(yīng)及對應(yīng)的平均風(fēng)速和有效波高取值，得到了以下結(jié)論：

(1) 基于Copula函數(shù)的極限狀態(tài)曲線構(gòu)造改進(jìn)方法及其數(shù)值算法將隨機(jī)環(huán)境變量的條件分布函數(shù)及其逆函數(shù)的求解轉(zhuǎn)化為條件Copula函數(shù)及其逆函數(shù)的求解，具有較高的可靠性和可操作性。

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

(3) 各方向?qū)Ч芗芷脚_(tái)基底剪力極限荷載效應(yīng)及對應(yīng)的平均風(fēng)速和有效波高取值差異較大。與全方向相比，極限狀態(tài)曲線較大方向的平均風(fēng)速、有效波高和基底剪力極限荷載效應(yīng)取值也較大。