基于Archimedean Copula 模型的風(fēng)浪聯(lián)合分布第二重現(xiàn)期

徐龍軍，陳祉宏，周道成，謝禮立

(1. 哈爾濱工業(yè)大學(xué)(威海)土木工程系，威海 264209；2. 大連理工大學(xué)深海工程研究中心，大連 116024）

基于概率的結(jié)構(gòu)設(shè)計方法中，環(huán)境荷載的計算與重現(xiàn)期的概念密切相關(guān)[1]．工程上，現(xiàn)有的單變量重現(xiàn)期計算方法不具備考慮荷載之間相關(guān)性的能力．在結(jié)構(gòu)的可靠度研究中，相關(guān)雙變量重現(xiàn)期的計算因不能真實地反映某一概率對應(yīng)的荷載組合，使得相關(guān)雙變量荷載的計算結(jié)果與實際情況差別較大，已不能滿足結(jié)構(gòu)基于可靠度設(shè)計的要求．因此，較準(zhǔn)確地估計相關(guān)荷載的重現(xiàn)期對結(jié)構(gòu)可靠度相關(guān)理論的發(fā)展和應(yīng)用均具有重要作用．

Copula 函數(shù)最初被應(yīng)用于金融、保險、生物建模等領(lǐng)域，但由于它能夠靈活地構(gòu)造多維分布，且當(dāng)隨機變量擬合成不同分布時，表示相關(guān)性的系數(shù)恒定不變，這種特性較好地解決了相關(guān)荷載構(gòu)造聯(lián)合概率分布問題．周道成等[2]利用Copula Gumbel 分布對潿洲島海洋站的極值風(fēng)速和有效波高實測數(shù)據(jù)進行了分析．劉德輔等[3]提出了二維泊松混合Copula Gumbel分布，并計算了嵊泗海區(qū)的臺風(fēng)海浪荷載．馮平等[4]利用Copula Gumbel 分布計算了在不同種類重現(xiàn)期下不同程度降雨和入境水遭遇組合的頻率．

本文首先討論了一般重現(xiàn)期和基于Archimedean Copula 中C-測度(C-measure)[5]的第二重現(xiàn)期(the secondary return period)[6]的基本特征．在對重現(xiàn)期安全域進行定義的基礎(chǔ)上，將第二重現(xiàn)期應(yīng)用到風(fēng)浪聯(lián)合概率分布的數(shù)據(jù)分析中，并與一般重現(xiàn)期及荷載效應(yīng)重現(xiàn)期進行了比較，認(rèn)為第二重現(xiàn)期較好地解決了相關(guān)雙變量重現(xiàn)期的計算問題．

1 Archimedean Copula函數(shù)

二元Copula 函數(shù)C(u，v)的定義及其性質(zhì)如下：C(u，v)的定義域為[0，1]2；對任意的u，v∈[0，1]，有C(u，0)＝0，C(u，1)＝u，C(0，v)＝0，C(1，v)＝v；C(u，v)為二維遞增函數(shù)．

Archimedean Copula 函數(shù)定義為：令函數(shù)φ 是一個連續(xù)的、嚴(yán)格遞減函數(shù)，定義域為[0，1]，值域為[0，∞]，φ(1)＝0，φ?1為函數(shù)φ 的反函數(shù)．稱函數(shù)φ 為Archimedean Copula 函數(shù)生成元，定義φ 的偽逆函數(shù)為φ[?1]，即

φ[?1](t)在[0，∞]為連續(xù)非增，在[0，φ(0)]為嚴(yán)格遞減，φ(φ[?1](t) ） ＝ min{t ，φ(0)} ，而 且 在[0 ，1] 上φ[?1](φ(u)） ＝u ，將 具 有 C(u ，v) ＝φ[?1](φ(u) ＋φ(v)）形式的Copula 函數(shù)稱為Archimedean Copula函數(shù)．

2 第二重現(xiàn)期

單變量荷載效應(yīng)對應(yīng)的重現(xiàn)期為

式中F(s)為結(jié)構(gòu)荷載效應(yīng)s 的概率分布函數(shù)．令ts＝F(s)，如果用等效荷載模型表示[7]，概率分布函數(shù)為

式中：fXY(x，y)可以看作風(fēng)速和波高的聯(lián)合概率密度函數(shù)；? 為積分區(qū)間，?＝ {(x，y)|Ax2＋By2≤s}．

將任意隨機變量超越其對應(yīng)的給定值的概率算入超越概率中，對應(yīng)的為聯(lián)合重現(xiàn)期T1，即

將僅當(dāng)所有隨機變量同時超越各自給定值時的概率算入超越概率中，對應(yīng)的為同現(xiàn)重現(xiàn)期T2，即

令 t1＝P(X≤x ，Y≤y) ，t2＝1?P(X≥x ，Y≥y)．如果隨機變量X 和Y 的聯(lián)合概率分布為Copula 函數(shù)，那么

圖1給出了T1和T2的等值線圖．圖2 為T1和T2的計算范圍，同心圓為聯(lián)合概率密度函數(shù)等值線，圖2(a)陰影部分的體積為1 ? t1，圖2(b)陰影部分的體積為1－t2．從圖2 中可以看出，對于同一個樣本(x，y)來說，1?t1＞1?t2，即t1＜t2，這時T1＜T2．當(dāng)(x，y)在原點時，1?t1＝1?t2，即t1＝t2，這時T1＝T2，由于把(x，y)取在原點是沒有討論意義的，所以，可以認(rèn)為T1＜T2．對于同一樣本，t1與t2的關(guān)系為u＋v?t1＝t2．

圖1 重現(xiàn)期T1和T2Fig.1 Return period for T1 and T2

圖2 重現(xiàn)期T1和T2的計算范圍Fig.2 Calculation ranges of return period for T1 and T2

某一樣本(x，y)的重現(xiàn)期T 的安全域? 定義為：在某種重現(xiàn)期T 下，(x，y)對應(yīng)的重現(xiàn)期值為將t 的積分區(qū)間稱作樣本(x，y)的重現(xiàn)期T 的安全域．

對于Ts，?s＝{(x，y)|Ax2＋By2≤s}，即第一象限橢圓與雙坐標(biāo)軸構(gòu)成的區(qū)域；對于T1，?1＝{(X，Y)|X≤x，Y≤y}，即圖2(a)中白色區(qū)域；對于T2，?2＝{(X，Y)|X≤x，Y≤y}∪{(X，Y)|X≥x，Y≤y}∪{(X，Y)|X≤x，Y≥y}，即圖2(b)中白色區(qū)域．

取T1等值線上另一個樣本(x*，y*)，其中x*＜x，y*＞y，如圖2(a)所示，(x*，y*)的安全域?1*為[0，x*]×[0，y*]，與(x，y)對應(yīng)的?1并不相同，?1里的部分樣本不會落在?1*里，?1*里的部分樣本也不會落在?1里，即當(dāng)(x，y)≠(x*，y*)時，?1≠ ?1*．所以T1的缺點為：①沒能考慮計算樣本聯(lián)合重現(xiàn)期安全域之外的安全樣本，將導(dǎo)致計算的重現(xiàn)期值小于實際荷載效應(yīng)重現(xiàn)期；②相同重現(xiàn)期值的不同樣本所包含的聯(lián)合重現(xiàn)期安全域不同，該定義不夠嚴(yán)格．

取T2等值線上另一個樣本(x*，y*)，其中x*＜x，y*＞y，如圖2(b)所示，(x*，y*)的安全域?2*為[0，x*]×[0，y*] [∪x*，∞] ×[0，y*] [0∪ ，x*]×[y*，∞]，?2之外的部分危險樣本落在?2*里，?2*之外的部分危險樣本也落在 ?2里，即當(dāng)(x ，y)≠(x*，y*) 時，?2≠?2*．所以T2的缺點為：①計算的同現(xiàn)重現(xiàn)期大于實際荷載效應(yīng)重現(xiàn)期，由同現(xiàn)重現(xiàn)期安全域中的[x，∞]×[0，y] [0∪ ，x]×[y，∞]部分可知，其定義范圍過大，當(dāng)一個隨機變量取無窮大時，結(jié)構(gòu)必定破壞；②相同重現(xiàn)期值的不同樣本所包含的同現(xiàn)重現(xiàn)期安全域不同，該定義不明確．

通過以上分析可知，對于同一個樣本，T1、T2和荷載效應(yīng)重現(xiàn)期Ts的關(guān)系為T1＜Ts＜T2．

指定Archimedean Copula 函數(shù)生成元，單側(cè)導(dǎo)數(shù)φ′(t-)和φ′(t＋)的定義域為(0，1]和[0，1)，令 t ∈ (0,1]，w=?(t)，n 為一不變的正整數(shù)，把[0,w ]等分為n 份{0,w / n,…,kw / n, …,w}，通過 w =?(t)把w 的等分點映射到[ t ,1]區(qū)間，結(jié)果為{t=t0,t1,…,tn?k,…,tn}，這里tn?k=?[?1](kw / n)，k＝0，1，…，n，如圖3 所示．

圖3 w 等分圖Fig.3 Divided w

因為 w＜?(0)，有

特別地

式中 VC(Rk)表示聯(lián)合概率密度函數(shù) Rk矩形內(nèi)的體積．再令

當(dāng)n →∞時

那么，圖5 中陰影部分的面積之和就可以表示成

由于KC(t)為集合{(u，v) [0∈ ，1]2|C(u，v)≤t}的C-測度，則第二重現(xiàn)期定義為[8]

于是，某一樣本的第二重現(xiàn)期安全域??就如圖5 中陰影部分所示．

圖4 C-測度Fig.4 C-measure

圖5 第二重現(xiàn)期? 的安全域Fig.5 Security domain of the secondary return period ?

對于圖2(a)中的樣本(,)x y ，它的重現(xiàn)期分別為?、T1與T2，??的物理意義可表述為所有重現(xiàn)期等于T1所對應(yīng)的無數(shù)個?1的并集．與T1和T2比較可以發(fā)現(xiàn)，??的邊界為一條曲線，比?1和?2的線性邊界更明確，而且在? 中相同重現(xiàn)期值的不同樣本對應(yīng)的??相同．

3 風(fēng)浪資料分析

采用數(shù)據(jù)來自潿洲島海洋站實測記錄的風(fēng)速和有效波高同步觀測資料[2]．應(yīng)用極值Ⅰ型分布擬合實際風(fēng)速、浪高概率分布，并選用Copula Gumbel 函數(shù)構(gòu)造聯(lián)合概率分布，其生成元為φ(t)＝[ ?ln(t)]α，隨機變量X 和Y 的聯(lián)合分布函數(shù)為

隨機變量X 和Y 的聯(lián)合概率密度函數(shù)為

式中α 為相關(guān)參數(shù)，它與Kendall 秩相關(guān)系數(shù)τ 的關(guān)系[8-11]為

隨機變量(X，Y)的n 組觀測值組成的樣本為[(x1，y1)，(x2，y2)，…，(xn，yn)]，那么其Kendall 秩相關(guān)系數(shù)[8-9]為

根據(jù)式(13)和式(14)可得

根據(jù)實測的風(fēng)浪資料，計算了子樣1 和子樣2 不同重現(xiàn)期(T1、T2、?、Ts)值分別為50、100,a 的等值線，見圖6．可以看出，? 與Ts的等值線均介于T1和T2之間．從安全域的角度看，在同一重現(xiàn)期值的條件下，隨著等值線上樣本的不同，?s與??的范圍不變．?1與?2則相反，不同的樣本對應(yīng)的范圍完全不同．T1等值線雖然完全包含了?s，但等值線上某一樣本(x，y)所對應(yīng)的?1僅僅包含矩形范圍[0，x] ×[0，y]，比?s小很多．而且無論怎么選取，在此樣本附近的樣本一定為?s之外的危險樣本；T2等值線雖然被?s完全包含在內(nèi)，但是等值線上某一樣本(x，y)所對應(yīng)的?2包含所對應(yīng)的幾個矩形疊加的范圍[0，x]×[0，y] [∪x，∞]×[0，y] [0∪ ，x]×[y，∞]，比?s大很多，而且無論怎么選取，在此樣本附近的樣本總比?s所限制的最大值小很多；與T1和T2不同，? 等值線所包含的范圍就是??，雖然其兩側(cè)的范圍超出了?s，但是在對角線附近，略小于?s，由風(fēng)浪相關(guān)性可知，其樣本主要分布在主對角線附近，因此與T1和T2相比，? 更適合描述風(fēng)浪重現(xiàn)期．

圖6 4種重現(xiàn)期等值線對比Fig.6 Comparison of isolines for four kinds of return period

進一步比較子樣1 和子樣2 的4 種重現(xiàn)期，結(jié)果見圖7．可以看出，? 總是位于T1和T2之間，最接近Ts．計算了各重現(xiàn)期的累積誤差ERR和均方誤差SS，列入表1 中．容易發(fā)現(xiàn)? 的ERR 和SS 均為最小值．

圖7 子樣樣本重現(xiàn)期比較Fig.7 Comparison between return periods of subexamples

表1 各子樣重現(xiàn)期誤差Tab.1 Errors for different return periods and subexamples a

將100,a、50,a、30,a、10,a 的Ts等值線和風(fēng)浪聯(lián)合概率密度fG(x，y)等值線繪制到一起，如圖8 所示．Ts等值線與fG(x，y)等值線相切的點為在此重現(xiàn)期值的條件下最接近fG(x，y)眾值的點，即最可能出現(xiàn)的風(fēng)浪組合值，表2 列出了不同重現(xiàn)期條件下子樣1 與子樣2 最可能出現(xiàn)的風(fēng)浪組合值．

表2 最可能出現(xiàn)的風(fēng)浪組合值Tab.2 The most possible wind and wave combination values

計算了表2 中各組合值在其他重現(xiàn)期(T1、T2、?)下的重現(xiàn)期值示于表3，可以看出? 最接近Ts．

圖8 Ts等值線和風(fēng)浪聯(lián)合概率密度fG(x，y)等值線Fig.8 Isolines for Ts and for wind-wave joint probability density fG(x，y)

表3 最可能出現(xiàn)風(fēng)浪組合值在T1、T2和? 中的重現(xiàn)期Tab.3 Return periods of T1，T2 and ? for the most possible wind and wave combination values a

采用同樣的方法，將T1、T2和? 在100,a、50,a、30,a、10,a 的最可能出現(xiàn)風(fēng)浪組合值求出來并列入表4 中．將表4 與表2 的結(jié)果進行了差值分析，結(jié)果見表5．可以看出，對于不同重現(xiàn)期對應(yīng)的最可能出現(xiàn)風(fēng)速值來說，? 的誤差在2.5%以內(nèi)，比T1、T2都小很多．對于不同重現(xiàn)期值的最可能出現(xiàn)波高值來說，?的最大誤差不超過7.7%，同樣比T1、T2?。?/p>

表4 T1、T2和 ? 最可能出現(xiàn)的風(fēng)浪組合值Tab.4 The most possible wind and wave combination values of T1，T2 and ?

表5 最可能出現(xiàn)風(fēng)浪組合值的比較Tab.5 Comparison of the most possible wind and wave combination values

4 結(jié) 語

為了較準(zhǔn)確計算風(fēng)浪聯(lián)合分布的荷載效應(yīng)，引入了基于Archimedean Copula 函數(shù)C-測度的第二重現(xiàn)期概念，第二重現(xiàn)期的曲線邊界具有比一般重現(xiàn)期線性邊界更確切的特點，根據(jù)定義的安全域概念，采用實測風(fēng)浪資料對四種重現(xiàn)期進行了對比分析．認(rèn)為：一般重現(xiàn)期中相同重現(xiàn)期的不同樣本所包含的安全域不同，而第二重現(xiàn)期和荷載效應(yīng)重現(xiàn)期的相同重現(xiàn)期對應(yīng)的不同樣本對應(yīng)的安全域相同，說明第二重現(xiàn)期比一般重現(xiàn)期定義嚴(yán)謹(jǐn)．實際風(fēng)浪資料分析也證明了這一點．根據(jù)實際風(fēng)浪組合值計算了4 種重現(xiàn)期并與荷載效應(yīng)重現(xiàn)期對比，第二重現(xiàn)期的累積誤差和均方誤差均為最小值．在最可能出現(xiàn)的風(fēng)浪組合值比較中，相同重現(xiàn)期條件下，第二重現(xiàn)期最可能出現(xiàn)風(fēng)速的誤差(2.5%) 和最可能出現(xiàn)波高的誤差(7.7%)都是最小的．說明第二重現(xiàn)期較一般重現(xiàn)期更接近荷載效應(yīng)重現(xiàn)期，其應(yīng)用簡單，而且能較準(zhǔn)確獲得海洋結(jié)構(gòu)在設(shè)計環(huán)境荷載重現(xiàn)期下承受的風(fēng)浪荷載．

［1］董 勝，侯 鵬，馮春健. 考慮設(shè)計壽命的風(fēng)浪聯(lián)合重現(xiàn)期標(biāo)準(zhǔn)初探[J]. 中國造船，2009，50(增)：215-221.Dong Sheng，Hou Peng，F(xiàn)eng Chunjian. Research on joint return period of wind speed and wave height considering project life[J]. Shipbuilding of China，2009，50(Suppl)：215-221(in Chinese).

［2］周道成，段忠東. 耿貝爾邏輯模型在極值風(fēng)速和有效波高聯(lián)合概率分布中的應(yīng)用[J]. 海洋工程，2003，21(2)：45-51.Zhou Daocheng，Duan Zhongdong. The Gumbel-Logistic model for joint probability distribution of extreme-value wind speeds and effective wave heights [J]. The Ocean Engineering，2003，21(2)：45-51(in Chinese).

［3］劉德輔，溫書勤，王利萍. 泊松-混合岡貝爾復(fù)合極值分布及其應(yīng)用[J]. 科學(xué)通報，2002，47(17)：1356-1360.Liu Defu，Wen Shuqin，Wang Liping. Poisson-Gumbel mixed compound extreme distrbution [J]. Chinese Science Bulletin and Its Application，2002，47(17)：1356-1360(in Chinese).

［4］馮 平，毛慧慧，王 勇. 多變量情況下的水文頻率分析方法及其應(yīng)用[J]. 水利學(xué)報，2009，40(1)：33-37.Feng Ping，Mao Huihui，Wang Yong. Method for hydrological reoccurrence frequency analysis under the condition of multivariate[J]. Journal of Hydraulic Engineering，2009，40(1)：33-37(in Chinese).

［5］Roger B N. An Introduction to Copulas [M]. 2nd ed.New York：Springer，2006.

［6］Salvadori G. Bivariate return periods via 2-copulas[J].Statistical Methodology，2004，l(1/2)：129-144.

［7］歐進萍，段忠東，肖儀清. 海洋平臺結(jié)構(gòu)安全評定：理論、方法與應(yīng)用[M]. 北京：科學(xué)出版社，2003.Ou Jinping，Duan Zhongdong，Xiao Yiqing. Safety Assessment of Offshore Platform Structure ： Theory ，Method and Application[M]. Beijing：Science Press，2003(in Chinese).

［8］Hollander M，Wolfe D A. Nonparametric Statistical Methods [M]. New York：John Wiley and Sons Ltd，1973.

［9］Lehmann E L. Nonparametrics ： Statistical Methods Based on Ranks [M]. San Francisco：HoldemDay，Inc，1975.

［10］Frees E W，Valdez E A. Understanding relationships using copulas[J]. North American Actuarial Journal，1998，2(1)：1-25.

［11］Patton A J. Skewness，Asymmetric Dependence，and Portfolios [Z]. London School of Economics and Political Science，2002.