基于Archimedean Copula 模型的风浪联合分布第二重现期

徐龙军，陈祉宏，周道成，谢礼立

(1. 哈尔滨工业大学(威海)土木工程系，威海 264209；2. 大连理工大学深海工程研究中心，大连 116024）

基于概率的结构设计方法中，环境荷载的计算与重现期的概念密切相关[1]．工程上，现有的单变量重现期计算方法不具备考虑荷载之间相关性的能力．在结构的可靠度研究中，相关双变量重现期的计算因不能真实地反映某一概率对应的荷载组合，使得相关双变量荷载的计算结果与实际情况差别较大，已不能满足结构基于可靠度设计的要求．因此，较准确地估计相关荷载的重现期对结构可靠度相关理论的发展和应用均具有重要作用．

Copula 函数最初被应用于金融、保险、生物建模等领域，但由于它能够灵活地构造多维分布，且当随机变量拟合成不同分布时，表示相关性的系数恒定不变，这种特性较好地解决了相关荷载构造联合概率分布问题．周道成等[2]利用Copula Gumbel 分布对涠洲岛海洋站的极值风速和有效波高实测数据进行了分析．刘德辅等[3]提出了二维泊松混合Copula Gumbel分布，并计算了嵊泗海区的台风海浪荷载．冯平等[4]利用Copula Gumbel 分布计算了在不同种类重现期下不同程度降雨和入境水遭遇组合的频率．

本文首先讨论了一般重现期和基于Archimedean Copula 中C-测度(C-measure)[5]的第二重现期(the secondary return period)[6]的基本特征．在对重现期安全域进行定义的基础上，将第二重现期应用到风浪联合概率分布的数据分析中，并与一般重现期及荷载效应重现期进行了比较，认为第二重现期较好地解决了相关双变量重现期的计算问题．

1 Archimedean Copula函数

二元Copula 函数C(u，v)的定义及其性质如下：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)为二维递增函数．

Archimedean Copula 函数定义为：令函数φ 是一个连续的、严格递减函数，定义域为[0，1]，值域为[0，∞]，φ(1)＝0，φ−1为函数φ 的反函数．称函数φ 为Archimedean Copula 函数生成元，定义φ 的伪逆函数为φ[−1]，即

φ[−1](t)在[0，∞]为连续非增，在[0，φ(0)]为严格递减，φ(φ[−1](t) ） ＝ min{t ，φ(0)} ，而 且 在[0 ，1] 上φ[−1](φ(u)） ＝u ，将 具 有 C(u ，v) ＝φ[−1](φ(u) ＋φ(v)）形式的Copula 函数称为Archimedean Copula函数．

2 第二重现期

单变量荷载效应对应的重现期为

式中F(s)为结构荷载效应s 的概率分布函数．令ts＝F(s)，如果用等效荷载模型表示[7]，概率分布函数为

式中：fXY(x，y)可以看作风速和波高的联合概率密度函数；Ω 为积分区间，Ω＝ {(x，y)|Ax2＋By2≤s}．

将任意随机变量超越其对应的给定值的概率算入超越概率中，对应的为联合重现期T1，即

将仅当所有随机变量同时超越各自给定值时的概率算入超越概率中，对应的为同现重现期T2，即

令 t1＝P(X≤x ，Y≤y) ，t2＝1−P(X≥x ，Y≥y)．如果随机变量X 和Y 的联合概率分布为Copula 函数，那么

图1给出了T1和T2的等值线图．图2 为T1和T2的计算范围，同心圆为联合概率密度函数等值线，图2(a)阴影部分的体积为1 − t1，图2(b)阴影部分的体积为1－t2．从图2 中可以看出，对于同一个样本(x，y)来说，1−t1＞1−t2，即t1＜t2，这时T1＜T2．当(x，y)在原点时，1−t1＝1−t2，即t1＝t2，这时T1＝T2，由于把(x，y)取在原点是没有讨论意义的，所以，可以认为T1＜T2．对于同一样本，t1与t2的关系为u＋v−t1＝t2．

图1 重现期T1和T2Fig.1 Return period for T1 and T2

图2 重现期T1和T2的计算范围Fig.2 Calculation ranges of return period for T1 and T2

某一样本(x，y)的重现期T 的安全域Ω 定义为：在某种重现期T 下，(x，y)对应的重现期值为将t 的积分区间称作样本(x，y)的重现期T 的安全域．

对于Ts，Ωs＝{(x，y)|Ax2＋By2≤s}，即第一象限椭圆与双坐标轴构成的区域；对于T1，Ω1＝{(X，Y)|X≤x，Y≤y}，即图2(a)中白色区域；对于T2，Ω2＝{(X，Y)|X≤x，Y≤y}∪{(X，Y)|X≥x，Y≤y}∪{(X，Y)|X≤x，Y≥y}，即图2(b)中白色区域．

取T1等值线上另一个样本(x*，y*)，其中x*＜x，y*＞y，如图2(a)所示，(x*，y*)的安全域Ω1*为[0，x*]×[0，y*]，与(x，y)对应的Ω1并不相同，Ω1里的部分样本不会落在Ω1*里，Ω1*里的部分样本也不会落在Ω1里，即当(x，y)≠(x*，y*)时，Ω1≠ Ω1*．所以T1的缺点为：①没能考虑计算样本联合重现期安全域之外的安全样本，将导致计算的重现期值小于实际荷载效应重现期；②相同重现期值的不同样本所包含的联合重现期安全域不同，该定义不够严格．

取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里，即当(x ，y)≠(x*，y*) 时，Ω2≠Ω2*．所以T2的缺点为：①计算的同现重现期大于实际荷载效应重现期，由同现重现期安全域中的[x，∞]×[0，y] [0∪ ，x]×[y，∞]部分可知，其定义范围过大，当一个随机变量取无穷大时，结构必定破坏；②相同重现期值的不同样本所包含的同现重现期安全域不同，该定义不明确．

通过以上分析可知，对于同一个样本，T1、T2和荷载效应重现期Ts的关系为T1＜Ts＜T2．

指定Archimedean Copula 函数生成元，单侧导数φ′(t-)和φ′(t＋)的定义域为(0，1]和[0，1)，令 t ∈ (0,1]，w=ϕ(t)，n 为一不变的正整数，把[0,w ]等分为n 份{0,w / n,…,kw / n, …,w}，通过 w =ϕ(t)把w 的等分点映射到[ t ,1]区间，结果为{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)表示联合概率密度函数 Rk矩形内的体积．再令

当n →∞时

那么，图5 中阴影部分的面积之和就可以表示成

由于KC(t)为集合{(u，v) [0∈ ，1]2|C(u，v)≤t}的C-测度，则第二重现期定义为[8]

于是，某一样本的第二重现期安全域Ως就如图5 中阴影部分所示．

图4 C-测度Fig.4 C-measure

图5 第二重现期ς 的安全域Fig.5 Security domain of the secondary return period ς

对于图2(a)中的样本(,)x y ，它的重现期分别为ς、T1与T2，Ως的物理意义可表述为所有重现期等于T1所对应的无数个Ω1的并集．与T1和T2比较可以发现，Ως的边界为一条曲线，比Ω1和Ω2的线性边界更明确，而且在ς 中相同重现期值的不同样本对应的Ως相同．

3 风浪资料分析

采用数据来自涠洲岛海洋站实测记录的风速和有效波高同步观测资料[2]．应用极值Ⅰ型分布拟合实际风速、浪高概率分布，并选用Copula Gumbel 函数构造联合概率分布，其生成元为φ(t)＝[ −ln(t)]α，随机变量X 和Y 的联合分布函数为

随机变量X 和Y 的联合概率密度函数为

式中α 为相关参数，它与Kendall 秩相关系数τ 的关系[8-11]为

随机变量(X，Y)的n 组观测值组成的样本为[(x1，y1)，(x2，y2)，…，(xn，yn)]，那么其Kendall 秩相关系数[8-9]为

根据式(13)和式(14)可得

根据实测的风浪资料，计算了子样1 和子样2 不同重现期(T1、T2、ς、Ts)值分别为50、100,a 的等值线，见图6．可以看出，ς 与Ts的等值线均介于T1和T2之间．从安全域的角度看，在同一重现期值的条件下，随着等值线上样本的不同，Ωs与Ως的范围不变．Ω1与Ω2则相反，不同的样本对应的范围完全不同．T1等值线虽然完全包含了Ωs，但等值线上某一样本(x，y)所对应的Ω1仅仅包含矩形范围[0，x] ×[0，y]，比Ωs小很多．而且无论怎么选取，在此样本附近的样本一定为Ωs之外的危险样本；T2等值线虽然被Ωs完全包含在内，但是等值线上某一样本(x，y)所对应的Ω2包含所对应的几个矩形叠加的范围[0，x]×[0，y] [∪x，∞]×[0，y] [0∪ ，x]×[y，∞]，比Ωs大很多，而且无论怎么选取，在此样本附近的样本总比Ωs所限制的最大值小很多；与T1和T2不同，ς 等值线所包含的范围就是Ως，虽然其两侧的范围超出了Ωs，但是在对角线附近，略小于Ωs，由风浪相关性可知，其样本主要分布在主对角线附近，因此与T1和T2相比，ς 更适合描述风浪重现期．

图6 4种重现期等值线对比Fig.6 Comparison of isolines for four kinds of return period

进一步比较子样1 和子样2 的4 种重现期，结果见图7．可以看出，ς 总是位于T1和T2之间，最接近Ts．计算了各重现期的累积误差ERR和均方误差SS，列入表1 中．容易发现ς 的ERR 和SS 均为最小值．

图7 子样样本重现期比较Fig.7 Comparison between return periods of subexamples

表1 各子样重现期误差Tab.1 Errors for different return periods and subexamples a

将100,a、50,a、30,a、10,a 的Ts等值线和风浪联合概率密度fG(x，y)等值线绘制到一起，如图8 所示．Ts等值线与fG(x，y)等值线相切的点为在此重现期值的条件下最接近fG(x，y)众值的点，即最可能出现的风浪组合值，表2 列出了不同重现期条件下子样1 与子样2 最可能出现的风浪组合值．

表2 最可能出现的风浪组合值Tab.2 The most possible wind and wave combination values

计算了表2 中各组合值在其他重现期(T1、T2、ς)下的重现期值示于表3，可以看出ς 最接近Ts．

图8 Ts等值线和风浪联合概率密度fG(x，y)等值线Fig.8 Isolines for Ts and for wind-wave joint probability density fG(x，y)

表3 最可能出现风浪组合值在T1、T2和ς 中的重现期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 的最可能出现风浪组合值求出来并列入表4 中．将表4 与表2 的结果进行了差值分析，结果见表5．可以看出，对于不同重现期对应的最可能出现风速值来说，ς 的误差在2.5%以内，比T1、T2都小很多．对于不同重现期值的最可能出现波高值来说，ς的最大误差不超过7.7%，同样比T1、T2小．

表4 T1、T2和 ς 最可能出现的风浪组合值Tab.4 The most possible wind and wave combination values of T1，T2 and ς

表5 最可能出现风浪组合值的比较Tab.5 Comparison of the most possible wind and wave combination values

4 结 语

为了较准确计算风浪联合分布的荷载效应，引入了基于Archimedean Copula 函数C-测度的第二重现期概念，第二重现期的曲线边界具有比一般重现期线性边界更确切的特点，根据定义的安全域概念，采用实测风浪资料对四种重现期进行了对比分析．认为：一般重现期中相同重现期的不同样本所包含的安全域不同，而第二重现期和荷载效应重现期的相同重现期对应的不同样本对应的安全域相同，说明第二重现期比一般重现期定义严谨．实际风浪资料分析也证明了这一点．根据实际风浪组合值计算了4 种重现期并与荷载效应重现期对比，第二重现期的累积误差和均方误差均为最小值．在最可能出现的风浪组合值比较中，相同重现期条件下，第二重现期最可能出现风速的误差(2.5%) 和最可能出现波高的误差(7.7%)都是最小的．说明第二重现期较一般重现期更接近荷载效应重现期，其应用简单，而且能较准确获得海洋结构在设计环境荷载重现期下承受的风浪荷载．

［1］董 胜，侯 鹏，冯春健. 考虑设计寿命的风浪联合重现期标准初探[J]. 中国造船，2009，50(增)：215-221.Dong Sheng，Hou Peng，Feng 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］周道成，段忠东. 耿贝尔逻辑模型在极值风速和有效波高联合概率分布中的应用[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］刘德辅，温书勤，王利萍. 泊松-混合冈贝尔复合极值分布及其应用[J]. 科学通报，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］冯 平，毛慧慧，王 勇. 多变量情况下的水文频率分析方法及其应用[J]. 水利学报，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］欧进萍，段忠东，肖仪清. 海洋平台结构安全评定：理论、方法与应用[M]. 北京：科学出版社，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.