基于<bold>Vine Copula</bold>函数的风浪要素联合概率分布模型

王望，朱金，康锐，李永乐; WANG Wang; ZHU Jin; KANG Rui; LI Yongle

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

摘要

随着全球气候变暖的加剧，极端气候现象发生的频率和强度均可能加大，这对海岸和近海结构的安全不利。基于中国东海的连云港海洋观测站实测风浪数据和Vine Copula理论，建立风浪要素中风速、波高、波浪周期、风向和波向五维随机变量之间的联合概率分布模型。采用极大似然法确定各风浪要素边缘分布模型参数，通过AIC信息准则和均方根误差RMSE进行拟合优度评价，由此建立风浪要素的边缘分布。采用带有基于残差的高斯似然函数的贝叶斯框架估计二维Copula函数的参数，结合AIC信息准则进行拟合优度评价并确定最优Copula函数。绘制最优联合分布概率密度图，与二维频率直方图进行对比以评价模型效果。采用Vine Copula函数建立多维联合概率模型并结合AIC值评价其拟合优度。研究结果表明：建立的Vine Copula模型可以较好地刻画风速、波高、波浪周期、风向和波向五维随机变量之间的联合概率分布。

关键词

风浪联合概率分布模型; 风浪荷载; 参数估计; 拟合优度检验

海洋环境恶劣，往往存在大风和巨浪的组合作用，同时，跨海大桥具有基础结构尺寸大、主跨轻柔、阻尼小、刚度小等特点，导致跨海大桥对风浪等海洋环境荷载非常敏感，极端风浪荷载成为威胁跨海大桥安全的主要因素之一。目前，跨海大桥多灾害作用的研究中通常研究单一灾害（如风、波浪、地震、冲刷等）对桥梁结构的作用，未能很好地考虑各灾害之间的相关性。特别是风浪要素，其包含了众多影响跨海大桥动力特性的参数，如风速、风向、波高、波向、波浪周期等，各参数是随时间变化的随机变量，且相互之间具有不同的尾部相关结构（即非线性）。因此，构造风浪要素的联合分布函数是研究海洋环境中风浪耦合效应及进一步探究风浪各要素对跨海大桥动力响应特征影响规律的基本前提。

Copula函数将多维联合分布分解为相应的边缘分布和Copula函数之积，可以灵活地表示两两随机变量之间的相关结构。近年来，Copula函数在土木工程领域逐渐得到关注。Li等^[

1]基于C-Vine Copula理论，研究了阿拉斯加南部海岸风浪相关性，并由建立的模型推算了特定重现期的长期荷载。Zhang等^{[参考文献 2

百度学术}2]采用椭圆Copula和三维Placket Copula建立了风速、风向和温度的三维联合分布。Bai等^{[参考文献 3

百度学术}3]采用混合二维Copula来描述风浪环境变量的相关性并与采用单一Copula的情况进行对比，研究表明，混合Copula能更好地描述多维变量之间复杂的相关关系。Wang等^{[参考文献 4

百度学术}4]基于Copula理论，提出了一种利用气象资料研究大跨桥梁风温联合作用的方法并应用于常泰长江大桥。Zhang等^{[参考文献 5

百度学术}5]基于非对称Copula函数，模拟了海洋环境要素的相关关系，侧重于捕捉环境要素的非对称相关关系。Li等^{[参考文献 6

百度学术}6]在浮式风力发电机组的概率疲劳评估研究中，基于Copula函数建立了多维概率模型，用来定义风浪环境参数相关性。Vanem^{[参考文献 7

百度学术}7]对特征波高和波浪周期的同时进行了多维分布研究，结果表明，使用非对称Copula函数能较好地模拟非对称的相关关系。Yang等^{[参考文献 8

百度学术}8]基于三维Copula函数建立了风速、风暴潮和暴雨的联合概率分布，提出了一种有效的PSO算法来估计边缘分布和联合分布的参数，并与极大似然法和对数似然法进行对比验证。Wang等^{[参考文献 9

百度学术}9]通过建立挪威Sulafjord的短期近海数据和桥址区数据之间的定量关系，推算了桥址区长期风浪数据，从而建立了风浪联合分布。Zhang等^{[参考文献 10

百度学术}10]通过基于Copula的多维概率模型，研究了考虑多种环境因素的近海结构长期荷载，由此建立了多种常见海况参数之间的联合概率模型。Xu等^{[参考文献 11

百度学术}11]将Copula函数应用于海洋结构可靠度的分析，描述了波高、波峰和平均风速3个近海环境参数之间的相关性，研究表明，该方法可用于考虑长期疲劳荷载和极端响应的海洋结构可靠度分析。陈子燊等^{[参考文献 12

百度学术}12]采用非对称Archimedean Copula函数分析极端波况下波高、风速和周期的三变量联合概率分布，计算了“或”重现期、“且”重现期和二次重现期，探讨了规范设计值的安全性并给出了建议。Copula理论不仅在海洋环境要素相关性研究方面有广泛的应用，在结构可靠度等其他领域也有重要的应用。陈建兵等^{[参考文献 13

百度学术}13]基于Copula理论，建立了混凝土参数之间的相关性模型，由生成的样本计算了混凝土本构全曲线。樊学平等^{[参考文献 14

百度学术}14]基于R-Vine Copula理论对大跨桥梁主梁检测点失效概率的相关性进行了研究，建立了检测点的相关性模型。刘月飞等^{[参考文献 15

百度学术}15]基于桥梁检测点的极值应力数据，建立了描述检测变量非线性相关的Vine Copula模型。宋帅等^{[参考文献 16

百度学术}16]将混合Copula方法应用于桥梁系统地震易损性分析中，准确描述了构件地震需求的非对称相关关系，简化了联合分布模型的建立过程。

综上所述，以上基于Copula函数的海洋环境参数联合分布的研究极大促进了跨海桥梁的建设，但是目前相关研究多针对二维及三维的海洋环境参数，这对准确模拟复杂多变的海洋环境来说是不够的。如前所述，海洋风浪要素中的风速、风向、波高、波向、波浪周期对于跨海大桥的动力响应均有重要影响，然而，目前针对风浪要素多维联合分布的研究还鲜有报道。笔者在单一Copula函数的基础上，基于Vine Copula理论建立了海洋风浪要素中风速、风向、波高、波向、波浪周期五维变量之间的联合分布模型，从而准确刻画了风浪要素之间的相关关系。首先，建立风浪各个要素的边缘概率模型，采用均方根误差（RMSE）进行拟合优度评价；在得到风浪各要素边缘分布的基础上，基于Copula理论，建立风浪要素两两之间的二维联合概率分布模型，通过AIC信息准则和均方根误差RMSE进行拟合优度评价，并考察风浪要素两两之间的相关性；基于Vine Copula理论，采用C-Vine结构构建了风浪要素中风速、风向、波高、波向、波浪周期五维变量之间的联合概率分布模型。通过AIC准则对模型进行拟合优度评价。

1 边缘分布

1.1 数据说明

采用位于中国东海的连云港海洋观测站2016—2020年波浪和风场观测数据，数据由中国国家科技资源共享服务平台——国家海洋科学数据中心(http://mds.nmdis.org.cn/)提供数据支撑。选用的风浪要素包括10 m高度处最小平均风速、特征波高、波浪周期、风向和波向，测量频率为每小时测量一次，站点的经纬度为34°47′0″N 119°26′0″E，最大风速为22 m，达到强风等级，最大波高为2.4 m。风浪数据信息如表1所示。需要说明的是，笔者在进行数据处理时发现，海洋站的观测数据中有很少一部分数据存在缺失的情况，即有个别或多个要素的数据没有观测到。针对这种情况，将缺失的样本数据予以剔除，尽可能多地保留其余数据样本，最后得到的样本总量为29 363个。

表1 风浪要素数据信息

Table 1 Information of wind and wave data

风浪要素	说明
风速	10 min平均风速，单位m/s
波高	特征波高，单位m
波周期	单位s
风向	单位度(°)，正北为0°，度数沿顺时针增加
波向	单位度(°)，正北为0°，度数沿顺时针增加

1.2 边缘分布

首先，需要建立风浪要素的边缘分布模型。研究中发现，风速、波高、波周期样本具有单峰分布的特征（图1（a）~（c）），采用常见的单峰分布模型进行拟合，包括Weibull、广义极值分布（GEV）、含有尺度参数和位置参数的t分布等。风向和波向具有多峰分布的特征（图1（d）、（e）），采用混合模型进行拟合，包括混合Gaussian分布、混合Gamma分布、混合Weibull分布。

1）Weibull分布

f (x, λ, k) = \frac{k}{λ} {(\frac{x}{λ})}^{k - 1} e^{- {(x / λ)}^{k}}

（1）

式中：λ为尺度参数；k为形状参数。

2）广义极值分布（GEV）

\begin{array}{l} F_{X} (x) = P (X < x) = \\ e x p - \{- {[1 - ξ (\frac{x - μ}{σ})]}^{1 / ξ}\}, ξ \neq 0 \end{array}

（2）

式中：ξ、β、μ为分布函数参数；ξ为形态参数；σ为尺度参数；μ为位置参数。

3）含有尺度参数和位置参数的t分布

f (x) = \frac{Γ (\frac{ν + 1}{2})}{σ \sqrt[]{ν π} Γ (\frac{ν}{2})} {[\frac{ν + {(\frac{x - μ}{σ})}^{2}}{ν}]}^{- (\frac{ν + 1}{2})}

（3）

式中：ν、σ、μ均为分布函数参数，其中：ν为形态参数；σ为尺度参数；μ为位置参数。

4）混合Gaussian模型

f (θ) = \sum_{l = 1}^{3} w_{l} \frac{1}{σ \sqrt[]{2 π}} e x p \{- \frac{1}{2} \frac{(θ - μ_{l})_{}^{2}}{σ_{l}^{2}}\}

（4）

式中： $\sum_{l = 1}^{3} w_{l} = 1$ ；μ_l、σ_l分别为相应部分Gaussian分布的均值和方差。

5）混合Gamma模型

f (θ) = \sum_{l = 1}^{c} w_{l} \frac{{β_{l}}^{α_{l}}}{Γ (α_{l})} θ^{α_{l} - 1} e^{- β_{l} θ}, θ > 0

（5）

式中： $\sum_{l = 1}^{3} w_{l} = 1$ ；α_l、β_l分别为相应部分Gamma分布的形状参数和尺度参数。

当采用上述概率分布模型对风浪要素进行拟合时，概率分布模型的参数估计采用极大似然法。另外，为了评价不同概率分布模型的拟合效果，采用AIC、BIC和RMSE对概率分布模型进行拟合优度评价，并据此选取最优的概率分布模型。AIC、BIC和RMSE的计算式为

A I C = - 2 \sum_{i = 1}^{n} l n f (x_{i}) + 2 k

（6）

式中：x_i为样本值；n为样本数量；f (x_i)为备选边缘分布函数的密度函数；k为备选边缘分布函数中分布参数的数量。

B I C = - 2 \sum_{i = 1}^{n} l n f (x_{i}) + k l n n

（7）

式中：x_i为样本值；n为样本数量；f (x_i)为备选边缘分布函数的密度函数；k为备选边缘分布函数中分布参数的数量。

R M S E = \sqrt[]{\frac{1}{n} {\sum_{i = 1}^{n} [P_{c} (i) - P_{0} (i)]}^{2}}

（8）

式中：n为样本数量；P_c为多维Copula联合分布理论频率值。RMSE值越小，拟合的效果越好。

通过观察风向概率直方图（图1（d））可以看出，风向的概率分布有3个较为明显的峰值，因此，采用的混合概率模型中应包含3个单峰分布，故采用混合维度为3的混合模型（3个单峰分布函数混合）来拟合。与风向类似，波向频率直方图（图1（e））也有多个峰值，分别采取二维和三维混合模型对波向的概率分布进行拟合，通过对比RMSE发现，三维混合模型的概率拟合效果较好。表2和表3给出了风速、波高、波浪周期、风向、波向5个风浪要素的最优边缘概率分布类型和相应参数。此外，图1还给出了这5个风浪要素样本的直方图及最优边缘概率分布曲线。由图1可以看出，选取的最优边缘概率分布曲线对5个风浪要素样本的拟合效果较好。拟合优度评价结果表明：风速、波高、波周期的最优拟合分布分别为Weibull分布、广义极值分布（GEV）、含有尺度参数和位置参数的t分布；而风向和波向的最优拟合分布均为混合Gaussian分布。

表2 风速、波高、波周期最优边缘概率分布

Table 2 The optimal marginal probability distribution of wind speed, wave height and wave period

环境要素	分布	BIC	AIC	尺度参数	形态参数	位置参数	RMSE
风速V_w	Weibull	140 950	140 930	5.496 2	1.728 8		0.006 3
波高H_s	Generalized Extreme Value	-7 914	-7 939	0.150 4	0.313 2	0.205 8	0.052 7
波周期T_p	七分布	101 130	101 110	0.959 7	3.090 1	4.863 0	0.015 3

表3 风向、波向边缘分布

Table 3 The marginal distribution of wind direction and wave direction

环境要素	函数类型	拟合参数										RMSE
风向θ_w	混合Gaussian	参数	w₁	w₂	w₃	μ₁	μ₂	μ₃	σ₁	σ₂	σ₃	0.006 3
	混合Gaussian	取值	0.149 7	0.552 3	0.298 0	0.461 0	2.080 3	4.923 6	0.225 5	0.793 8	0.672 0	0.006 3
	混合Gamma	参数	w₁	w₂	w₃	λ₁	λ₂	λ₃	k₁	k₂	k₃	0.007 3
	混合Gamma	取值	0.381 1	0.346 2	0.272 6	0.212	0.827	0.075 9	11.428 3	1.388 9	65.938 1	0.007 3
波向θ_s	混合Gaussian	参数	w₁	w₂	w₃	μ₁	μ₂	μ₃	σ₁	σ₂	σ₃	0.046 5
	混合Gaussian	取值	0.646 6	0.212 8	0.140 7	1.173 2	1.590 6	5.202 6	0.677 6	0.283 1	0.564 1	0.046 5
	混合Weibull	参数	w₁	w₂	w₃	λ₁	λ₂	λ₃	k₁	k₂	k₃	0.048 5
	混合Weibull	取值	0.799 2	0.142 4	0.058 4	1.533 5	5.433 6	0.000 1	2.748 1	10.909 7	63.643 4	0.048 5

（a）风速

（b）波高

（c）波周期

（d）风向

（e）波向

图1 风浪要素概率分布直方图及最优边缘概率分布曲线

Fig. 1 Histogram of wind and wave data and the optimal marginal probability distribution curve

2 二维联合分布

2.1 二维Copula理论

根据Sklar定理^[

17]，多维联合分布和其边缘分布可以写为

F（x₁ … x_i）=C（F₁（x₁）… F_n（x_n） ; θ）

（9）

式中：F_i(x_i)为随机变量x_i的边缘概率分布函数；C(·)为Copula函数；θ为Copula函数的参数。

将海洋环境变量定义为一个n维随机变量X=（x₁, x₂… x_i… x_n），基于Copula理论，联合概率密度可以表示为

f （x₁, x₂… x_i… x_n） =c（F₁（x₁）, F₂（x₂） … F_i（x_i） … F_n（x_n））

\prod_{i = 1}^{n} f (x_{i})

（10）

式中：F_i（x_i）、f （x_i）分别为随机变量x_i的边缘概率分布的分布函数和概率密度函数；c为Copula密度函数。

由式（10）可得二维联合分布概率密度函数公式

f （x₁,x₂）=c（F₁（x₁）, F₂（x₂））f （x₁）f （x₂）

（11）

式中：c为二维Copula函数的密度函数。

由式（11）可知，为了研究风浪要素之间的二维联合概率分布，首先应建立风浪要素两两变量之间的Copula函数。选取风速—波高、波高—波周期、风速—风向、波高—波向4个随机变量对来研究并建立二维Copula函数。目前，研究风浪要素联合分布的大多数文献均将这些变量对作为研究对象。其次，这些变量对中的两个变量之间往往存在较大的相关性，而这些相关性在进行结构（如海上风机、跨海桥梁等）的动力响应分析时又至关重要。样本数据经过式（12）可以变换为范围0~1的标准分布，样本的标准分布忽略了样本数值大小的差异，只保留数据的相对大小关系，以便于更清晰地研究数据之间的相关性。图2是风速—波高经验Copula频率直方图，可以看出，风速和波高具有较强的尾部相关性。

U = \frac{n - R + 0.5}{n}

（12）

式中：U为将样本变换为范围为0~1的标准分布后的随机变量；n为样本个数；R为样本点在所有样本中的排序。

图2 风速—波高标准分布二维频率直方图

Fig. 2 Binary frequency histogram of standard distribution of wind and wave

常用的二维Copula函数类型有：Gaussian、T、Clayton、Gumbel、Frank，如表4所示。

表4 5种典型的Copula函数

Table 4 Five typical Copula functions

Copula函数	函数公式	ρ范围
Gaussian	Φ(Φ^-1(u₁),Φ^-1(u₂)\| ρ)	[-1,1]
T	t_ρ,n(t_n^-1(u₁),t_n^-1(u₂)\|ρ)	[-1,1]
Clayton	(u₁^-ρ+u₂^-ρ-1)^-1^/ρ	(0, ∞)
Gumbel	exp(-((-ln u₁) ^ρ+(-ln u₂)^ρ)¹^/ρ)	[1, ∞)
Frank	$- \frac{1}{ρ} l n (1 + \frac{(e^{- ρ u_{1}} - 1) (e^{- ρ u_{2}} - 1)}{e^{- ρ} - 1})$	(-∞,∞)/｛0｝

2.2 拟合优度评价

要进行拟合优度评价，首先要计算经验Copula，经验Copula可以通过式（13）计算^[

18]。

C_{E} (u_{1}, u_{2}) = \frac{1}{n} \sum_{i = 1}^{n} I (u_{1 i} \leq u_{1}, u_{2 i} \leq u_{2})

（13）

式中：n为样本的大小，对每一个1≤i≤n，满足 $u_{1 i} \leq u_{1}, u_{2 i} \leq u_{2}$ 时，I ( $u_{1 i} \leq u_{1}, u_{2 i} \leq u_{2}$ )=1。

对应三维的情况，经验Copula可通过式（14）计算。

C_{E} (u_{1}, u_{2}, u_{3}) = \frac{1}{n} \sum_{i = 1}^{n} I (u_{1 i} \leq u_{1}, u_{2 i} \leq u_{2}, u_{3 i} \leq u_{3})

（14）

式中：n为研究样本的大小，对每一个1≤i≤n，满足 $u_{1 i} \leq u_{1}, u_{2 i} \leq u_{2}, u_{3 i} \leq u_{3}$ 时，I( $u_{1 i} \leq u_{1},$ $u_{2 i} \leq u_{2}$ ， $u_{3 i} \leq u_{3}$ )=1。

为了对选取的Copula函数进行拟合优度评价，选用均方根误差（RMSE）作为评价标准来评价模型的优劣。

R M S E = \sqrt[]{\frac{1}{n} {\sum_{i = 1}^{n} |\overset{}{C_{E}} [(u_{1})_{i}, (u_{2})_{i}] - C [(u_{1})_{i}, (u_{2})_{i}]|}^{2}}

（15）

式中：C_E为Copula经验值；C为建立的模型计算出的理论值。

S_n越小,说明拟合的Copula模型与经验Copula越接近，拟合效果越好。

建立二维联合分布可以对多种Copula函数分别进行参数估计，采用贝叶斯框架和基于残差的高斯似然函数进行参数估计^[

19]。

贝叶斯分析已在多个领域应用于参数估计。当获得新信息时，贝叶斯理论更新假设的先验概率，将所有建模的不确定因素归因于参数，通过式（16）估计模型参数的后验分布。

p (θ | \tilde{Y}) = \frac{p (θ) p (\tilde{Y} | θ)}{p (\tilde{Y})}

（16）

式中： $p (θ)$ 和 $p (θ | \tilde{Y})$ 分别表示参数的先验和后验分布； $p (\tilde{Y} | θ)$ 可代表似然函数； $p (\tilde{Y}) = \int_{θ} p (\tilde{Y} | θ) d θ$ 为证据，证据在每次建模中都是常数。

在缺乏参数先验分布的有效信息时，可以采用均匀先验，假设残差不相关，同方差、均值为零的高斯分布，那么似然函数就可以通过式（17）表示。

L (θ | \tilde{Y}) = \prod_{i = 1}^{n} \frac{1}{\sqrt[]{2 π {\tilde{σ}}^{2}}} e x p \{- \frac{1}{2} {\tilde{σ}}^{- 2} {[\tilde{y_{i}} - y_{i} (θ)]}^{2}\}

（17）

式中： $\tilde{σ}$ 为测量误差的标准偏差估计值； $\tilde{y_{i}}$ 为观测变量的联合概率； $y_{i} (θ)$ 为Copula预测概率。为了简洁和数值稳定，式（17）可以采用对数转换成为

l (θ | \tilde{Y}) = - \frac{n}{2} l n (2 π) - \frac{n}{2} l n {\tilde{σ}}^{2} - \frac{1}{2} {\tilde{σ}}^{- 2} {\sum_{i = 1}^{n} [\tilde{y_{i}} - y_{i} (θ)]}^{2}

（18）

通过AIC、BIC及RMSE对不同类型的Copula模型进行拟合优度评价，拟合优度评价结果见表5，参数估计值见表6。风速—波高的最优联合分布为θ=0.846 1的Gaussian Copula；波高—波浪周期的最优联合分布为θ=0.465 3的Gaussian Copula；风速—风向的最优联合分布为θ=-1.145 7的Frank Copula；波高—波向的最优联合分布为θ=0.541 0的Gaussian Copula。图3为θ=0.846 1的Gaussian Copula函数图像。

表5 二维Copula拟合优度评价

Table 5 The performance evaluation of 2-dimensional Copula distribution

随机变量	Copula类型	AIC	BIC	ML	RMSE
风速-波高	Gaussian	-191 776	-191 768	95 889	6.54
	Clayton	-187 416	-187 408	93 709	7.05
	Frank	-190 322	-190 314	95 162	6.71
	Gumbel	-190 648	-190 640	95 325	6.67
	Joe	-187 587	-187 579	93 795	7.02
	Plackett	-189 264	-189 256	94 633	6.83
波高-周期	Gaussian	-197 802	-197 794	98 902	5.90
	Clayton	-193 411	-193 402	96 706	6.36
	Frank	-197 103	-197 095	98 552	5.97
	Gumbel	-196 943	-196 934	98 472	5.99
	Joe	-195 444	-195 435	97 722	6.15
	Plackett	-196 352	-196 343	98 176	6.05
风速-风向	Gaussian	-271 710	-271 701	135 855	1.68
	Clayton	-225 937	-225 929	112 969	3.66
	Frank	-274 369	-274 361	137 185	1.60
	Gumbel	-225 919	-225 911	112 960	3.66
	Joe	-225 954	-225 946	112 978	3.66
	Plackett	-274 571	-274 562	137 286	1.60
波高-波向	Gaussian	-180 147	-180 139	90 074	7.97
	Clayton	-176 902	-176 894	88 452	8.43
	Frank	-178 757	-178 749	89 379	8.16
	Gumbel	-179 820	-179 812	89 910	8.02
	Joe	-179 358	-179 350	89 680	8.08
	Plackett	-178 097	-178 089	89 049	8.26

表6 最优二维Copula分布及拟合参数

Table 6 The optimal 2-dimensional Copula distribution and fitting parameters

联合分布	Copula类型	θ
风速-波高	Gaussian	0.846 1
波高-周期	Gaussian	0.465 3
风速-风向	Frank	-1.145 7
波高-波向	Gaussian	0.541 0

图3 Gaussian Copula（θ=0.846 1）概率密度图

Fig. 3 Probability density of Gaussian Copula (θ=0.846 1)

由图2和图3可以看出，Gaussian Copula能较好地模拟风速—波高的相关性。此外，样本数据的二维频率直方图可以表示联合分布的经验图像，与图4不同的是，图3只表示两随机变量之间的相关关系，略去了两随机变量之间数值大小的影响，通过图3能更清晰地了解二者的相关关系。图4则表示出两随机变量之间联合分布的图像，保留了样本数值大小的信息。通过二维频率直方图和联合分布的概率密度图可以直观比较本文选取的二维Copula模型的拟合效果，如图4所示。由图4可以看出，对于风速—波高、波高—周期、风速—波向，其二维联合概率密度图与经验图像的峰值和形状都比较接近，说明相关性模拟较为合理，较好地考虑了风浪要素之间的相关性。同时，从图4也可以看到波，高—波向的联合分布与经验直方图的峰值大小存在一定差异，这是由于该样本波高—波向之间的相关性比较复杂，采用单一的Copula函数模拟它们之间相关性的效果有限，在后续研究中，可以采用混合Copula函数来进行模拟，并优化拟合效果。由于不同海域的风浪要素之间的相关关系可能会有较大的差异，本研究侧重于探究建立风浪联合分布的方法和思路，因此未对波高—波向的二维联合概率模型作进一步优化。

（a）风速-波高频率直方图

（b）风速-波高（Gaussian，θ=0.846 1）

（c）波高-波周期频率直方图

（d）波高-波周期（Gaussian，θ=0.465 3）

（e）风速-风向频率直方图

（f）风速-风向（Frank，θ=-1.145 7）

（g）波高-波向频率直方图

（h）波高-波向（Gaussian，θ=0.541）

图4 风浪要素两两变量之间的频率直方图及最优二维Copula函数

Fig. 4 Binary frequency histogram of wind and wave data and the corresponding optimal 2-dimensional Copula distribution

3 多维联合分布

可以采用多维Copula函数建立多维随机变量的联合分布，但可供选用的多维Copula函数类型有限，且灵活性较弱。Joe^[

20]提出了采用Vine Copula来构建多维随机变量的联合分布的方法，Vine Copula具有结构多样灵活的优点。

3.1 Pair-Copula理论

Bedford等^[

21]提出了通过Pair-Copula构建多维随机变量联合分布概率模型的方法，基于条件概率可以将多维随机变量分解成一系列的Pair-Copula结构，分解时采用的逻辑结构不同，可以构造出不同的模型。Pair-Copula结构推动了Copula理论在多维随机变量应用中的发展。

根据条件概率，多维联合分布概率密度函数可以表示为

\begin{array}{l} f (x_{1} \dots x_{n}) = f_{n} (x_{n}) \cdot f (x_{n - 1} | x_{n}) \cdot \\ f (x_{n - 2} | x_{n - 1}, x_{n}) \dots f (x_{1} | x_{2} \dots x_{n}) \end{array}

（19）

由式（19）可以得到二维随机变量的联合概率密度函数为

f (x_a, x_j)=c_{a j} (F_a(x_a) , F_j(x_j) ) f (x_a) f (x_j)

（20）

式中：a, j=1, 2 … n，且a≠j；c_{a j}(F_a(x_a), F_j(x_j))为x_a和x_j的二维Copula密度函数。

由式（20）可以推导x_j已知的情况下x_a的条件概率密度

f (x_a|x_j)= c_{a j}(F_a(x_a), F_j(x_j)) f (x_a)

（21）

由式（21）可得在n维随机变量u已知的条件下，任意随机变量x的条件密度函数

f (x | u) = c_{x u_{a} | u_{- a}} (F (x | u_{- a}), F (u_{a} | u_{- a})) \cdot f (x | u_{- a})

（22）

式中：u_a是n维随机变量u中的一个分量；u_-a是n维随机变量u中去掉u_a之后的n-1维分量。

以三维联合分布为例说明分解过程，根据式（19），三维的联合分布函数可写为

f (x_{1}, x_{2}, x_{3}) = f_{1} (x_{3}) \cdot f (x_{2} | x_{3}) \cdot f (x_{1} | x_{2}, x_{3})

（23）

由式（22）可以得到

f (x_{1} | x_{2}, x_{3}) = c_{13 | 2} {F (x_{1} | x_{2}), F (x_{3} | x_{2})} \cdot f (x_{1} | x_{2})

（24）

将f（x₁|x₂）= c₁₂｛F₁ (x₁), F₂(x₂)｝f₁(x₁)带入式（24）可得

f (x_{1} | x_{2}, x_{3}) = c_{13 | 2} {F (x_{1} | x_{2}), F (x_{3} | x_{2})} \cdot c_{12} {F (x_{1}), F (x_{1})} \cdot f (x_{1})

（25）

将f（x₂ | x₃）= c₂₃｛F₂ (x₂), F₃(x₃)｝f₂(x₂)和式（25）代入式（23）即可将多维联合分布分解为Pair-Copula结构和边缘分布的概率密度函数的乘积

\begin{array}{l} f (x_{1}, x_{2}, x_{3}) = f_{3} (x_{3}) \cdot c_{23} (F_{2} (x_{2}), F_{3} (x_{3})) \cdot f_{2} (x_{2}) \\ c_{13 | 2} {F (x_{1} | x_{2}), F (x_{3} | x_{2})} \cdot c_{12} {F_{1} (x_{1}), F_{2} (x_{1})} \cdot f (x_{1}) \end{array}

（26）

3.2 Vine结构

高维Copula函数可以通过R-Vine结构来建立。基于两两随机变量之间的相依组合，结合条件概率可以建立多维Copula函数。R-Vine中有两类特殊的Vine：C-Vine和D-Vine，这两类Vine有各自的逻辑结构，用于建立高维变量之间的联合分布。

C-Vine结构的特点是每层树都有一个主节点，主节点连接其他所有节点。为方便阐述构造原理，图5给出了四维的C-Vine结构图，该C-Vine的结构有3棵树，每棵树有一个主节点，主节点和其他节点（node）相连形成边（edge），边在下一棵树中作为主节点，以此类推，直到连接所有树的节点。对任意维n维C-Vine，有n-1棵树，有n(n-1)/2条边，每条边对应一个Pair-Copula结构。

图5 四维C-Vine分解结构

Fig. 5 Decomposition structure for four-dimensionalC-Vine model

D-Vine结构的特点是每层树的节点依次相连，呈直线状。图6给出了四维D-Vine的结构，该D-Vine结构有3棵树，每棵树中节点一次连接，形成边，边在下一棵树中成为节点，用同样的方式类推，直至连接所有节点。对任意维n维D-Vine，同样有n-1棵树，有n(n-1)/2条边，每条边对应一个Pair-Copula结构。

图6 四维D-Vine分解结构

Fig. 6 Decomposition structure for four-dimensionalD-Vine model

C-Vine的多变量分解结构的确定需要先确定根节点和其他节点的顺序，C-Vine的根节点通常选取与其他变量相关性最强的节点，可以采用Aas等^[

22]与Dißmann等^{[参考文献 23

百度学术}23]提出的序惯估计法来确定。

1）从第1层树开始，计算所有随机变量两两组合的经验Kendall相关系数τ，如表7所示，选取和其他随机变量的相关系数τ之和最大的随机变量作为根节点，这样就确定了第1层树的结构。

表7 经验Kendall’s τ矩阵及Kendall’s τ之和

Table 7 The empirical Kendall’s τ and the sum of it

参数	风速	波高	波周期	风向	波向	τ_sum
风速	1	0.415 3	0.081 4	0.147 6	0.064 2	1.708 5
波高	0.415 3	1	0.037 5	0.274 5	0.162 5	1.889 8
波周期	0.081 4	0.037 5	1	0.008 6	0.001 6	1.129 1
风向	0.147 6	0.274 5	0.008 6	1	0.166 9	1.597 6
波向	0.642 0	0.162 5	0.001 6	0.166 9	1	1.973 0

2）选择第1层树中二维随机变量的Copula函数种类并使用极大似然法估计参数θ。

3）结合第2）步确定的Copula函数及参数θ，通过式（25）将原随机变量转换为条件随机变量，即树2的随机变量。

F (x | v) = \frac{\partial C_{x v_{j} | v_{- j}} (F (x | v_{- j}), F (v_{j} | v_{- j}))}{\partial F (v_{j} | v_{- j})}

（27）

式中：v_-j为v中除去v_j的n-1维向量，F(x|v)为条件分布函数， $C_{x v_{j} | v_{- j}}$ 为连接 $F (x | v_{- j})$ 与 $F (v_{j} | v_{- j})$ 的二维Copula函数。

4）使用确定树1结构的方法确定剩余的所有树。

采用C-Vine Copula构建多维联合分布，图7给出了最优C-Vine的构造，1~5分别表示风速、波高、周期、风向、波向，五维C-Vine Copula有4层树，在第1层树中，随机变量2（波高）与其余4个变量分别通过二维Copula连接，形成4条边；在第2层树中，这4条边作为节点，选取一个主节点，再一次通过二维Copula进行连接，以此类推，直至连接所有节点。

图 7 维C-Vine Copula结构

Fig. 7 The optimal five-dimensional C-Vine model

由式（10）和图7可得，五维随机变量(x₁,x₂,x₃,x₄,x₅)的C-Vine联合概率密度函数如式（28）所示。

f(x₁,x₂,x₃,x₄,x₅)=f₁(x₁)f₂(x₂)f₃(x₃)f₄(x₄)f₅(x₅)c₁₂(F₁(x₁),F₂(x₂))×c₂₃(F₂(x₂),F₃(x₃))×c₂₄(F₂(x₂),F₄(x₄))c₂₅(F₂(x₂),F₅(x₅)) × c_13|2(F_1|2(x₁|x₂),F_3|2(x₃|x₂)) × c_14|2(F_1|2(x₁|x₂),F_4|2(x₄|x₂))×c_15|2(F_1|2(x₁|x₂),F_5|2(x₅|x₂))·c_34|12(F_3|12(x₃|x₁,x₂),F_4|12(x₄|x₁,x₂))×c_45|12(F_4|12(x₄|x₁,x₂),F_5|12(x₅|x₁,x₂))×c_35|124(F_3|124(x₃|x₁,x₂,x₄),F_5|124(x₅|x₁,x₂,x₄))

（28）

式中：F_i(·)为每个随机变量的累积分布函数；F_·|·(·|·)为条件分布函数；c_·，·(·,·)为Copula密度函数；c_·|·(·|·)为Copula密度函数。

可以通过比较AIC值进行C-Vine Copula模型的拟合优度评价，AIC的值越小说明拟合效果越好。一般来说，C-Vine和D-Vine均可用于构建数据的联合分布，具体采用哪一种模型应该取决于数据本身，通常应通过试算来最终确定较优的模型。分别采用C-Vine和D-Vine来构建风浪联合分布模型，并对比两种模型的AIC值，发现C-Vine能更好地模拟该海洋站点的风浪联合分布模型。因此，最终采用C-Vine来构建风浪的联合分布模型。此外，采用序惯估计法得到风浪联合概率分布的C-Vine模型。通过序惯估计法得出的最优根节点顺序为2、1、4、5、3（1~5分别表示风速、波高、周期、风向、波向），AIC值为-27 508，如表8所示。采用遍历所有根节点顺序的方法（即考虑了所有可能的根节点顺序），对由序惯估计法得出的C-Vine风浪联合概率分布模型进行拟合优度评价。五维随机变量的根节点顺序共有 $A_{5}^{3} = 60$ 种，针对每一种根节点顺序，首先采用AIC准则建立相应的C-Vine模型，并通过AIC值的大小对60个C-Vine模型进行排序，结果如表9所示。限于篇幅，仅给出了前20个较优的C-Vine模型及相应的根节点顺序。通过比较发现，通过序惯估计法建立的C-Vine模型的AIC值在60种情况中排第3，并且其AIC值与前面2个模型的AIC值非常接近（相差0.31%以内）。因此，序惯估计法能高效且准确的找出较优的C-Vine联合概率模型，该模型能很好地描述风浪要素间的相关关系。随着变量维度的提升，根节点的组合会迅速增加，采用遍历所有根节点顺序的方法会极大地增加计算负担，由此也体现出序惯估计法的优越性。

表8 最优C-Vine参数及AIC值

Table 8 The optimal C-Vine parameters and AIC values

树	连接边	Copula	参数1	参数2	AIC	总AIC
T1	2，1	BB8	2.07	0.99	-13 473	-27 508
	2，3	t	0.05	3.54	-2 143
	2，4	rotated BB8 Copula (270 degrees)	-2.21	-0.83	-5 116
	2，5	t	-0.23	3.6	-2 570
T2	1，3\|2	rotated Tawn type 2 Copula (270 degrees)	-1.45	0.37	-2 024
	1，4\|2	t	0.04	9.06	-348
	1，5\|2	BB8	1.18	0.92	-376
T3	4，5\|2，1	rotated BB8 Copula (180 degrees; “survival BB8”)	1.48	0.85	-996
T3	4，3\|2，1	BB8	1.25	0.88	-436
T4	3，5\|2，1，4	Gaussian	0.03		-22

表9 不同C-Vine模型的根节点顺序和AIC值

Table 9 The root node order and AIC value of different C-Vine models

根节点顺序	AIC	根节点顺序	AIC
2、1、5、4、3	-27 593.8	2、3、4、5、1	-27 338.8
2、1、3、5、4	-27 555.6	2、5、4、3、1	-27 244.4
2、1、4、5、3	-27 508.1	3、2、1、5、4	-26 133.6
2、5、1、4、3	-27 454.9	3、2、5、4、1	-26 128.8
2、4、1、3、5	-27 424.2	3、2、4、5、1	-26 080.5
2、4、3、5、1	-27 409.7	5、2、1、3、4	-26 048.4
2、3、1、5、4	-27 385.6	5、2、3、4、1	-26 029.7
2、5、3、4、1	-27 384.4	1、2、4、3、5	-26 007.7
2、3、5、4、1	-27 375.2	1、2、5、4、3	-25 977.3
2、4、5、3、1	-27 368.1	1、2、3、5、4	-25 945.7

通过建立的C-Vine模型仿真了风速、波高、波浪周期、风向和波向五维随机变量之间的累积概率密度（CDF）。为了更直观地观察拟合效果，以风速、波高、风向为例，结合式（14），图8给出了其三维联合分布的累积概率密度，黑色网格是原始数据的经验Copula图像，彩色图像是由C-Vine模拟产生的图像。由图8可知，建立的C-Vine模型能较好地模拟风速、波高、风向的相关关系。

图8 三维累积概率密度图

Fig. 8 Diagram of three-dimensional cumulative probability density

4 结论

基于Vine Copula函数研究了中国东海连云港海洋观测站的风浪要素之间的联合概率分布，得出以下结论：

1）风速、波高、波周期的概率分布为单峰分布，最优拟合分布分别为Weibull分布、广义极值分布、t分布；而风向和波向的概率分布为多峰分布，最优拟合分布均为混合Gaussian分布。

2）风浪要素中两两随机变量之间的联合概率分布研究表明，风速-波高、波高—波周期、风速—风向、波高—波向4个二维变量对的最优二维联合概率分布分别为Gaussian Copula、Gaussian Copula、Frank Copula和Gaussian Copula。

3）建立的C-Vine模型可以较好地刻画风速、波高、波浪周期、风向和波向五维随机变量之间的联合概率分布。

4）采用Vine Copula函数建立了东海连云港海洋观测站风浪要素之间的联合概率分布模型，对于中国其他海域海洋观测站风浪要素之间的联合分布规律还有待进一步研究。但该研究方法和思路可为中国其他海域海洋观测站风浪要素之间的联合分布研究提供借鉴。

致谢

感谢中国国家科技资源共享服务平台——国家海洋科学数据中心(http://mds.nmdis.org.cn/)提供数据支撑。参考文献

参考文献

LI X, ZHANG W. Long-term assessment of a floating offshore wind turbine under environmental conditions with multivariate dependence structures [J]. Renewable Energy, 2020, 147: 764-775. [百度学术]

ZHANG W M, WANG Z W, LIU Z. Joint distribution of wind speed, wind direction, and air temperature actions on long-span bridges derived via trivariate metaelliptical and plackett copulas [J]. Journal of Bridge Engineering, 2020, 25(9): 04020069. [百度学术]

BAI X Y, JIANG H, LI C, et al. Joint probability distribution of coastal winds and waves using a logtransformed kernel density estimation and mixed copula approach [J]. Ocean Engineering, 2020, 216: 107937. [百度学术]

WANG Z W, ZHANG W M, TIAN G M, et al. Joint values determination of wind and temperature actions on long-span bridges: Copula-based analysis using longterm meteorological data [J]. Engineering Structures, 2020, 219: 110866. [百度学术]

ZHANG Y, KIM C W, BEER M, et al. Modeling multivariate ocean data using asymmetric copulas [J]. Coastal Engineering, 2018, 135: 91-111. [百度学术]

LI X, ZHANG W. Probabilistic fatigue evaluation of floating wind turbine using combination of surrogate model and copula model [C]//AIAA Scitech 2019 Forum. San Diego, California. Reston, Virginia: AIAA, 2019. [百度学术]

VANEM E. Joint statistical models for significant wave height and wave period in a changing climate [J]. Marine Structures, 2016, 49: 180-205. [百度学术]

YANG X, QIAN J. Joint occurrence probability analysis of typhoon-induced storm surges and rainstorms using trivariate Archimedean copulas [J]. Ocean Engineering, 2019, 171: 533-539. [百度学术]

WANG J G, LI L, JAKOBSEN J B, et al. Metocean conditions in a Norwegian fjord in relation to the floating bridge design [J]. Journal of Offshore Mechanics and Arctic Engineering, 2019, 141(2): 021604. [百度学术]

ZHANG Y, BEER M, QUEK S T. Long-term performance assessment and design of offshore structures [J]. Computers & Structures, 2015, 154: 101-115. [百度学术]

XU X, LAI F, LI G J, et al. A novel vibration suppression device for floating offshore wind generator [C]//Proceedings of ASME 2018 12th International Conference on Energy Sustainability Collocated With the ASME 2018 Power Conference and the ASME 2018 Nuclear Forum. Lake Buena Vista, Florida, USA: 2018. [百度学术]

陈子燊, 路剑飞, 于吉涛. 基于非对称Archimedean Copula的三变量风浪重现水平分析[J]. 海洋通报, 2017, 36(6): 631-637. [百度学术]

CHEN Z S, LU J F, YU J T. Analysis on return levels of trivariate stormy waves based on asymmetric Archimedean Copula function [J]. Marine Science Bulletin, 2017, 36(6): 631-637. (in Chinese) [百度学术]

陈建兵, 陶金聚, 任晓丹, 等. 基于Copula理论的混凝土受压本构全曲线参数相关性研究[J]. 土木工程学报, 2020, 53(7): 52-63. [百度学术]

CHEN J B, TAO J J, REN X D, et al. Research on the dependence between parameters in complete stress-strain constitutive curve of concrete based on copula theory [J]. China Civil Engineering Journal, 2020, 53(7): 52-63. (in Chinese) [百度学术]

樊学平, 杨光红, 肖青凯, 等. 大跨桥梁主梁失效概率分析的最优R-Vine Copula [J]. 吉林大学学报(工学版), 2021, 51(4): 1296-1305. [百度学术]

FAN X P, YANG G H, XIAO Q K, et al. Optimal R-Vine Copula information fusion for failure probability analysis of long-span bridge girder [J]. Journal of Jilin University (Engineering and Technology Edition), 2021, 51(4): 1296-1305. (in Chinese) [百度学术]

刘月飞, 樊学平. 失效非线性相关的桥梁截面可靠性Vine-Copula数据融合[J]. 同济大学学报(自然科学版), 2019, 47(3): 315-321. [百度学术]

LIU Y F, FAN X P. Data fusion about Vine-Copula for bridge section reliability considering nonlinear correlation of failure modes [J]. Journal of Tongji University (Natural Science), 2019, 47(3): 315-321. (in Chinese) [百度学术]

宋帅, 钱永久, 吴刚. 基于多元Copula函数的桥梁体系地震易损性分析方法研究[J]. 振动与冲击, 2017, 36(9): 122-129, 208. [百度学术]

SONG S, QIAN Y J, WU G. Seismic fragility analysis of a bridge system based on multivariate Copula function [J]. Journal of Vibration and Shock, 2017, 36(9): 122-129, 208. (in Chinese) [百度学术]

SKLAR A. Fonctions de repartition a n dimensions et leurs marges [J]. Publications de l,Institut de statistique de l'Université de Paris, 1959, 8: 229-231. [百度学术]

CARNICERO J A, AUSÍN M C, WIPER M P. Non-parametric copulas for circular-linear and circular-circular data: An application to wind directions [J]. Stochastic Environmental Research and Risk Assessment, 2013, 27(8): 1991-2002. [百度学术]

SADEGH M, RAGNO E, AGHAKOUCHAK A. Multivariate Copula Analysis Toolbox (MvCAT): Describing dependence and underlying uncertainty using a Bayesian framework [J]. Water Resources Research, 2017, 53(6): 5166-5183. [百度学术]

JOE H. Families of m-variate distributions with given margins and m(m-1)/2 bivariate dependence parameters [M]//Institute of Mathematical Statistics Lecture Notes-Monograph Series. Hayward, CA: Institute of Mathematical Statistics, 1996: 120-141. [百度学术]

BEDFORD T, COOKE R M. Probability density decomposition for conditionally dependent random variables modeled by vines [J]. Annals of Mathematics and Artificial Intelligence, 2001, 32(1/2/3/4): 245-268. [百度学术]

AAS K, CZADO C, FRIGESSI A, et al. Pair-copula constructions of multiple dependence [J]. Insurance: Mathematics and Economics, 2009, 44(2): 182-198. [百度学术]

DIßMANN J, BRECHMANN E C, CZADO C, et al. Selecting and estimating regular vine copulae and application to financial returns [J]. Computational Statistics & Data Analysis, 2013, 59: 52-69. [百度学术]

基于Vine Copula函数的风浪要素联合概率分布模型 PDF

摘要

关键词