基于C藤Copula的多維隨機環境變量極限狀態曲面

涂志斌， 姚劍鋒， 黃銘楓， 樓文娟

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

EC(environmental contour,EC)法是多維隨機環境變量聯合作用下結構極限荷載效應估計的有效方法，其中心思想是將結構響應與環境變量分開，僅在具有指定超越概率的極限狀態曲面上搜索荷載效應極值作為結構極限荷載效應[1]，在海洋工程[2]、地震工程[3]和風工程[4]中有廣泛應用。

極限狀態曲面是具有相同超越概率的多維隨機環境變量組合點的集合，其計算基礎為多維隨機環境變量的聯合分布。常用的聯合分布模型有條件聯合分布模型、傳統聯合分布模型和Copula函數。在條件聯合分布模型中聯合分布表達為主變量的邊緣分布和次變量的條件分布，并將條件分布的參數視為主變量的函數，如Huseby等[5-6]、Jonathan等[7]和Saranyasoontorn等[8-9]的研究。由于條件分布的參數與主變量的函數關系不具有普遍性，條件聯合分布模型的適用范圍受到了限制。傳統聯合分布模型包含皮爾遜III型分布、Nataf分布等多種分布模型，其特點是將一維分布拓展為多維分布，采用線性相關系數表達變量之間的相關結構，如董勝等[10]、Yue[11]和Silva-González等[12]的研究。然而多維隨機變量可能服從不同種類的邊緣分布，不滿足線性相關關系，傳統聯合分布模型的應用也具有一定的局限性。Copula函數是一種建立聯合分布的工具，可建立具有不同邊緣分布和任意相關關系的多維隨機變量聯合分布函數，并將Gumbel分布、Nataf分布等多種分布納入其中[13]。涂志斌等[14]采用二維Copula函數建立高層建筑二維風荷載分量的聯合分布。Zhang等[15-16]和黃銘楓等[17]采用多種高維Copula函數建立了多風向極值風速的聯合分布。為彌補多維Copula函數采用同種函數結構表達各變量間的相關關系這一缺陷，Bedford等[18-19]提出了藤Copula。藤Copula將多維Copula密度函數表達為一系列二維pair Copula密度函數的乘積，且每個pair Copula可選用不同的二維Copula函數。常用的藤Copula有C藤Copula和D藤Copula[20]，其中C藤Copula在海洋工程、土木工程和金融工程等領域的應用更為廣泛。Montes-Iturrizaga等[21]以有效波高為主變量，平均風速和譜峰周期為次變量，采用3種二維Copula鏈接成C藤Copula建立聯合分布；樊學平等[22]采用二維Gaussian Copula鏈接成C藤Copula 考察了天津富民橋主梁5個截面極值應力的相關性；陳好杰等[23]采用5種二維Copula鏈接成多種C藤Copula討論了3個大型風電場群之間的出力相關性；李磊等[24]利用C藤Copula估計了連續幾個交易日收益率之間的自相依結構。以上研究表明當多維隨機變量由一個主變量和多個次變量組成時，基于C藤Copula的多維隨機變量聯合分布模型能較好地表達各變量之間的相關結構。

極限狀態曲面構造的一般方法為Rosenblatt變換。該方法先根據可靠指標在標準正態空間中構造極限狀態曲面，再將其映射到物理空間中。根據該變換Silva-González等、Valamanesh 等[25]和Li等[26]構造了不同海域平均風速、有效波高和譜峰周期的極限狀態曲面。然而采用Rosenblatt變換構造極限狀態曲面涉及條件分布函數的求解和求逆，當多維隨機變量的聯合分布較為復雜時，以上計算難以實施。

本文提出了基于C藤Copula函數的多維隨機變量極限狀態曲面計算公式及數值算法，結合風浪同步觀測數據構造了平均風速、有效波高和譜峰周期的極限狀態曲面，提出了根據極限狀態曲面等值線設置風浪參數組合值，討論了不同Copula函數對極限狀態曲面和組合值的影響，并得到幾點有益的結論可供海洋結構設計參考。

1 Copula函數的基本理論

設X=(X1,X2,…,Xn)為n維隨機環境變量，其邊緣分布為Fi(xi)(i=1,…,n)；n維隨機變量U=(U1,U2,…,Un)在[0,1]n范圍內獨立均勻分布，且ui=Fi(xi)；若Fi(xi)單調連續，那么根據Sklar定理，X的聯合分布函數F(x)可表達為[20]：

F(x)=C[F1(x1),F2(x2),…,Fn(xn)]=

C(u1,u2,…,un)

(1)

表1 常用Copula函數

設c(u)為Copula密度函數，其表達式為

(2)

X的聯合概率密度函數f(x)與c(u)的關系為

f(x)=f(x1,x2,…,xn)=

(3)

(4)

2 基于C藤Copula的聯合分布模型

2.1 聯合分布模型

藤是藤Copula的鏈接方式，由樹、節點和邊構成，不同藤的邏輯結構有所不同。對于5維(n=5)隨機變量，C藤的邏輯結構如圖1所示。

圖1 n=5時C藤的邏輯結構圖

(5)

式中，cj,j+i|1,…,j-1為變量Xj和Xj+i基于變量X1,…,Xj-1的二維Copula密度函數，表達了變量Xj|X1,…,Xj-1與變量Xj+i|X1,…,Xj-1的相關結構。將式(5)代入式(3)，n維隨機變量的聯合概率密度函數可表達為

(6)

(7)

特殊地，對于二維隨機變量，式(7)可寫作

(8)

式中，α21為F1(x1)和F2(x2)的相關參數。對于三維隨機變量，式(7)可寫作

(9)

2.2 參數估計

?

(10)

?

(11)

2.3 最優Copula

在眾多Copula函數中，為C藤中的每一條邊指定最優Copula作為二維pair Copula。在采用相同邊緣分布的前提下，最優Copula能準確表達變量之間的相關結構。目前常用的最優Copula評價準則為AIC準則，其計算公式為

(12)

3 基于C藤Copula的極限狀態曲面

設Z=(Z1,Z2,…,Zn)為標準正態空間中相互獨立的n維隨機變量。當超越概率為Pe時，在標準正態空間中極限狀態曲面可表達為

|z|=β

(13)

式中，β為可靠指標，β=-Φ-1(Pe)。對于三維變量，在球坐標系中式(13)可寫作

z1=βsinθcosφ

z2=βsinθsinφ

z3=βcosθ

(14)

式中，θ和φ分別仰角和方位角。根據Rosenblatt變換，物理空間中的極限狀態曲面與標準正態空間中的極限狀態曲面一一對應，且具有相同的超越概率，可表達為

(15)

其中ei=Φ(zi)。對式(15)進行改寫，可得：

(16)

將式(7)代入式(16)

e1=F1(x1)

?

(17)

在式(17)中，多維變量的聯合分布模型由C藤Copula函數建立，Rosenblatt變換表達為二維Copula的偏導數并采用數值算法求解，避免了式(15)中高維條件分布函數的求解和求逆，具有良好的可行性。

4 工程應用

4.1 風浪實測數據及樣本選取

本文以紐約東海岸附近63115號海洋觀測站(北緯40.50°，西經-72.92°)的風參數和波浪參數觀測記錄為統計樣本。各參數同步記錄，每小時記錄一次。提取有效波高日極值Hs及平均風速v、譜峰周期Tp、風向φv及波向φH同步值作為分析樣本。圖2為風向與波向之差(φv-φH)的角直方圖。圖2表明風向和波向的差異較小。根據IEC規范[29]和Zhang等[30]的建議，當風向和波向差異較小時，在風浪聯合作用分析中可認為風向和波向保持一致。將波向歸類到8個角度區間，其玫瑰圖如圖3所示。由圖3可知，波向記錄在SE、S、SW和W方向上的頻率較高，分別為14.28%、21.48%、16.98%和16.71%，而在NW、N、NE和E方向上的頻率較低，分別為8.75%、8.48%、7.81%和5.51%。由于波向與風向的偏差較小，圖3近似表達了風向在各方向上的分布頻率。由于S方向風浪參數的分布頻率較大，選擇該方向的統計數據作為分析樣本。

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

圖3 波向玫瑰圖

圖4為S方向有效波高日極值與平均風速、譜峰周期同步值散點圖，表2為3變量的Pearson線性相關系數。由表2可知，有效波高與平均風速有較強的正相關性，與譜峰周期有較弱的正相關性，平均風速與譜峰周期有較弱的負相關性。因此采用C藤Copula構造風浪參數的聯合分布模型時取(X1,X2,X3)=(Hs,v,Tp)。

圖4 S方向有效波高日極值及平均風速、譜峰周期同步值散點圖

表2 S方向風浪參數的相關系數

4.2 邊緣分布擬合

采用Lognormal分布擬合有效波高的邊緣分布，采用Weibull分布擬合平均風速和譜峰周期的邊緣分布，表達式如式(18)和(19)，擬合結果如圖5所示。

Lognormal：

(18)

Weibull：

(19)

(a) 有效波高

4.3 風浪參數聯合分布模型

表3 風浪參數聯合分布模型參數估計值及AIC值

4.4 風浪極限狀態曲面

由于有效波高日極值的統計時長為24 h，重現期T與超越概率Pe的換算關系為

Pe=1/365.25/T

(20)

因此T=20 a時Pe=1.37×10-4，可靠指標β=-Φ-1(Pe)=3.64。在[0,π)和[0,2π)范圍內對仰角θ和方位角φ進行均勻離散，離散點數分別為d=60、r=120。分別根據數值算法和參考文獻[12]計算基于C藤Copula和三維Gaussian Copula的三維風浪參數極限狀態曲面，結果如圖6和7所示。曲面上各點為T=20 a時可能發生的平均風速、有效波高和譜峰周期組合。由于相關結構不同，兩種Copula構造的極限狀態曲面不同，特別是在中部右側處，前者凸出，后者凹陷。考察曲面形狀對風浪參數極值的影響，結果如表4。其中Hs,20、v20和Tp,20分別為T=20 a時有效波高、平均風速和譜峰周期的極值。由表4可知，基于C藤Copula的平均風速極值和譜峰周期極值略微偏小，原因是在C藤Copula中平均風速和譜峰周期的相關性度量樣本為基于有效波高的條件經驗分布。但總體而言，基于C藤Copula和基于三維Gaussian Copula的風浪參數極值基本與邊緣分布一致，受曲面形狀的影響較小。

圖6 T=20 a時基于C藤Copula的風浪參數極限狀態曲面

圖7 T=20 a時基于三維Gaussian Copula的風浪

表4 T=20 a時的風浪參數極值

圖8為重現期T=20 a時基于C藤Copula的風浪參數極限狀態曲面等值線，最大輪廓線為二維風浪參數極限狀態曲線，分別由對應的二維最優Copula構造。在圖8(a)、(b)中，由于在二維和三維風浪參數聯合分布模型中(Hs,v)、(Hs,Tp)的相關結構均采用Galambos Copula表達，相關性度量樣本均為各變量的經驗分布，等值線分布在二維曲線形成的外包絡內；在圖8(c)中盡管在二維和三維聯合分布模型中(v,Tp)的相關結構均采用Frank Copula表達，但前者的相關性度量樣本為經驗分布，后者為條件經驗分布，導致等值線不再分布在二維曲線內，而是略有超出。

(a) 固定Tp

圖9為重現期T=20 a時基于三維Gaussian Copula的風浪參數極限狀態曲面等值線，最大輪廓線與圖8相同。由于三維Gaussian Copula由二維Gaussian Copula拓展而來，任意兩變量的相關結構均由二維Gaussian Copula表達，與二維最優Copula的相關結構差異較大，在圖9(a)～(c)中等值線和二維曲線明顯不同。造成圖8、9不同的原因是，C藤Copula由各變量的最優Copula鏈接而來，在局部和整體均能最優表達變量間的相關性，由此形成的極限狀態曲面有效納入了二維變量的相關結構，其等值線的外包絡是二維極限狀態曲線；而三維Gaussian Copula僅在整體上最優表達了各變量的相關性，不能顧及局部，由此形成的極限狀態曲面也不能準確納入二維變量的相關結構。因此對于三維及高維變量，更適合采用C藤Copula建立聯合分布模型。

(a) 固定Tp

在EC法中，結構荷載效應計算僅在極限狀態曲面上可能產生極限荷載效應的隨機變量組合點處進行。為了保證風、浪參數各自及聯合重現期，本文提出在使平均風速取極值的等值線上設置風浪參數組合值。根據基于C藤Copula和三維Gaussian Copula的風浪參數極限狀態曲面，提取各自平均風速極值所在的等值線，重現期T=20 a時有效波高、平均風速和譜峰周期的組合值Hs,c、vc和Tp,c如表5所示，其中err為各組合值與基于邊緣分布的極值(見表4)之間的誤差。由表5可知：①vc即為基于C藤Copula或三維Gaussian Copula的平均風速極值，與基于邊緣分布的極值間的誤差較小；② 基于C藤Copula的Hs,c有所降低，幅度不超過9%；③ 基于C藤Copula的Tp,c顯著降低，幅度接近60%；④ 與風浪參數極值相比，采用極限狀態曲面等值線確定的組合值在保證多維變量各自及聯合重現期的前提下有所降低，可使結構設計更加經濟合理；⑤ 基于三維Gaussian Copula的Hs,c和Tp,c均顯著降低，表明極限狀態曲面的形狀對組合值有較大影響，準確地構造極限狀態曲面十分必要。

表5 T=20 a時風浪參數組合值

5 結 論

本文提出了基于C藤Copula函數的多維隨機變量極限狀態曲面計算公式及其數值算法。在該方法中聯合分布模型由C藤Copula函數建立，Rosenblatt變換表達為二維Copula的偏導數并采用數值算法求解。結合風浪同步觀測數據構造了平均風速、有效波高和譜峰周期的極限狀態曲面，提出了根據極限狀態曲面確定風浪參數組合值，討論了不同Copula的影響，并得到了以下結論：

(1) 與多維Copula函數相比，C藤Copula函數在局部和整體上均能最優表達變量間的相關結構，更適用于建立多維隨機變量的聯合分布模型。

(2) 基于C藤Copula函數的極限狀態曲面有效納入了二維變量的相關結構，其等值線的外包絡是二維極限狀態曲線。

(3) 極限狀態曲面的形狀對風浪參數組合值有較大的影響。與風浪參數極值相比，根據等值線確定的組合值在保證多維變量各自及聯合重現期的前提下有所降低，可使結構設計更加經濟合理。