基于Vine Copula函數(shù)的風(fēng)浪要素聯(lián)合概率分布模型

王望,朱金,康銳,李永樂

（西南交通大學(xué) 橋梁工程系，成都 610031）

海洋環(huán)境惡劣，往往存在大風(fēng)和巨浪的組合作用，同時(shí)，跨海大橋具有基礎(chǔ)結(jié)構(gòu)尺寸大、主跨輕柔、阻尼小、剛度小等特點(diǎn)，導(dǎo)致跨海大橋?qū)︼L(fēng)浪等海洋環(huán)境荷載非常敏感，極端風(fēng)浪荷載成為威脅跨海大橋安全的主要因素之一。目前，跨海大橋多災(zāi)害作用的研究中通常研究單一災(zāi)害（如風(fēng)、波浪、地震、沖刷等）對(duì)橋梁結(jié)構(gòu)的作用，未能很好地考慮各災(zāi)害之間的相關(guān)性。特別是風(fēng)浪要素，其包含了眾多影響跨海大橋動(dòng)力特性的參數(shù)，如風(fēng)速、風(fēng)向、波高、波向、波浪周期等，各參數(shù)是隨時(shí)間變化的隨機(jī)變量，且相互之間具有不同的尾部相關(guān)結(jié)構(gòu)（即非線性）。因此，構(gòu)造風(fēng)浪要素的聯(lián)合分布函數(shù)是研究海洋環(huán)境中風(fēng)浪耦合效應(yīng)及進(jìn)一步探究風(fēng)浪各要素對(duì)跨海大橋動(dòng)力響應(yīng)特征影響規(guī)律的基本前提。

Copula 函數(shù)將多維聯(lián)合分布分解為相應(yīng)的邊緣分布和Copula 函數(shù)之積，可以靈活地表示兩兩隨機(jī)變量之間的相關(guān)結(jié)構(gòu)。近年來，Copula 函數(shù)在土木工程領(lǐng)域逐漸得到關(guān)注。Li 等[1]基于C-Vine Copula理論，研究了阿拉斯加南部海岸風(fēng)浪相關(guān)性，并由建立的模型推算了特定重現(xiàn)期的長期荷載。Zhang 等[2]采用橢圓Copula 和三維Placket Copula 建立了風(fēng)速、風(fēng)向和溫度的三維聯(lián)合分布。Bai 等[3]采用混合二維Copula 來描述風(fēng)浪環(huán)境變量的相關(guān)性并與采用單一Copula 的情況進(jìn)行對(duì)比，研究表明，混合Copula 能更好地描述多維變量之間復(fù)雜的相關(guān)關(guān)系。Wang 等[4]基于Copula 理論，提出了一種利用氣象資料研究大跨橋梁風(fēng)溫聯(lián)合作用的方法并應(yīng)用于常泰長江大橋。Zhang 等[5]基于非對(duì)稱Copula 函數(shù)，模擬了海洋環(huán)境要素的相關(guān)關(guān)系，側(cè)重于捕捉環(huán)境要素的非對(duì)稱相關(guān)關(guān)系。Li 等[6]在浮式風(fēng)力發(fā)電機(jī)組的概率疲勞評(píng)估研究中，基于Copula 函數(shù)建立了多維概率模型，用來定義風(fēng)浪環(huán)境參數(shù)相關(guān)性。Vanem[7]對(duì)特征波高和波浪周期的同時(shí)進(jìn)行了多維分布研究，結(jié)果表明，使用非對(duì)稱Copula 函數(shù)能較好地模擬非對(duì)稱的相關(guān)關(guān)系。Yang 等[8]基于三維Copula 函數(shù)建立了風(fēng)速、風(fēng)暴潮和暴雨的聯(lián)合概率分布，提出了一種有效的PSO 算法來估計(jì)邊緣分布和聯(lián)合分布的參數(shù)，并與極大似然法和對(duì)數(shù)似然法進(jìn)行對(duì)比驗(yàn)證。Wang 等[9]通過建立挪威Sulafjord 的短期近海數(shù)據(jù)和橋址區(qū)數(shù)據(jù)之間的定量關(guān)系，推算了橋址區(qū)長期風(fēng)浪數(shù)據(jù)，從而建立了風(fēng)浪聯(lián)合分布。Zhang 等[10]通過基于Copula 的多維概率模型，研究了考慮多種環(huán)境因素的近海結(jié)構(gòu)長期荷載，由此建立了多種常見海況參數(shù)之間的聯(lián)合概率模型。Xu 等[11]將Copula 函數(shù)應(yīng)用于海洋結(jié)構(gòu)可靠度的分析，描述了波高、波峰和平均風(fēng)速3 個(gè)近海環(huán)境參數(shù)之間的相關(guān)性，研究表明，該方法可用于考慮長期疲勞荷載和極端響應(yīng)的海洋結(jié)構(gòu)可靠度分析。陳子燊等[12]采用非對(duì)稱Archimedean Copula 函數(shù)分析極端波況下波高、風(fēng)速和周期的三變量聯(lián)合概率分布，計(jì)算了“或”重現(xiàn)期、“且”重現(xiàn)期和二次重現(xiàn)期，探討了規(guī)范設(shè)計(jì)值的安全性并給出了建議。Copula 理論不僅在海洋環(huán)境要素相關(guān)性研究方面有廣泛的應(yīng)用，在結(jié)構(gòu)可靠度等其他領(lǐng)域也有重要的應(yīng)用。陳建兵等[13]基于Copula理論，建立了混凝土參數(shù)之間的相關(guān)性模型，由生成的樣本計(jì)算了混凝土本構(gòu)全曲線。樊學(xué)平等[14]基于R-Vine Copula 理論對(duì)大跨橋梁主梁檢測點(diǎn)失效概率的相關(guān)性進(jìn)行了研究，建立了檢測點(diǎn)的相關(guān)性模型。劉月飛等[15]基于橋梁檢測點(diǎn)的極值應(yīng)力數(shù)據(jù)，建立了描述檢測變量非線性相關(guān)的Vine Copula 模型。宋帥等[16]將混合Copula 方法應(yīng)用于橋梁系統(tǒng)地震易損性分析中，準(zhǔn)確描述了構(gòu)件地震需求的非對(duì)稱相關(guān)關(guān)系，簡化了聯(lián)合分布模型的建立過程。

綜上所述，以上基于Copula函數(shù)的海洋環(huán)境參數(shù)聯(lián)合分布的研究極大促進(jìn)了跨海橋梁的建設(shè)，但是目前相關(guān)研究多針對(duì)二維及三維的海洋環(huán)境參數(shù)，這對(duì)準(zhǔn)確模擬復(fù)雜多變的海洋環(huán)境來說是不夠的。如前所述，海洋風(fēng)浪要素中的風(fēng)速、風(fēng)向、波高、波向、波浪周期對(duì)于跨海大橋的動(dòng)力響應(yīng)均有重要影響，然而，目前針對(duì)風(fēng)浪要素多維聯(lián)合分布的研究還鮮有報(bào)道。筆者在單一Copula 函數(shù)的基礎(chǔ)上，基于Vine Copula理論建立了海洋風(fēng)浪要素中風(fēng)速、風(fēng)向、波高、波向、波浪周期五維變量之間的聯(lián)合分布模型，從而準(zhǔn)確刻畫了風(fēng)浪要素之間的相關(guān)關(guān)系。首先，建立風(fēng)浪各個(gè)要素的邊緣概率模型，采用均方根誤差（RMSE）進(jìn)行擬合優(yōu)度評(píng)價(jià)；在得到風(fēng)浪各要素邊緣分布的基礎(chǔ)上，基于Copula 理論，建立風(fēng)浪要素兩兩之間的二維聯(lián)合概率分布模型，通過AIC 信息準(zhǔn)則和均方根誤差RMSE 進(jìn)行擬合優(yōu)度評(píng)價(jià)，并考察風(fēng)浪要素兩兩之間的相關(guān)性；基于Vine Copula理論，采用C-Vine結(jié)構(gòu)構(gòu)建了風(fēng)浪要素中風(fēng)速、風(fēng)向、波高、波向、波浪周期五維變量之間的聯(lián)合概率分布模型。通過AIC 準(zhǔn)則對(duì)模型進(jìn)行擬合優(yōu)度評(píng)價(jià)。

1 邊緣分布

1.1 數(shù)據(jù)說明

采用位于中國東海的連云港海洋觀測站2016—2020 年波浪和風(fēng)場觀測數(shù)據(jù)，數(shù)據(jù)由中國國家科技資源共享服務(wù)平臺(tái)——國家海洋科學(xué)數(shù)據(jù)中心(http://mds.nmdis.org.cn/)提供數(shù)據(jù)支撐。選用的風(fēng)浪要素包括10 m 高度處最小平均風(fēng)速、特征波高、波浪周期、風(fēng)向和波向，測量頻率為每小時(shí)測量一次，站點(diǎn)的經(jīng)緯度為34°47′0″N 119°26′0″E，最大風(fēng)速為22 m，達(dá)到強(qiáng)風(fēng)等級(jí)，最大波高為2.4 m。風(fēng)浪數(shù)據(jù)信息如表1 所示。需要說明的是，筆者在進(jìn)行數(shù)據(jù)處理時(shí)發(fā)現(xiàn)，海洋站的觀測數(shù)據(jù)中有很少一部分?jǐn)?shù)據(jù)存在缺失的情況，即有個(gè)別或多個(gè)要素的數(shù)據(jù)沒有觀測到。針對(duì)這種情況，將缺失的樣本數(shù)據(jù)予以剔除，盡可能多地保留其余數(shù)據(jù)樣本，最后得到的樣本總量為29 363 個(gè)。

表1 風(fēng)浪要素?cái)?shù)據(jù)信息Table 1 Information of wind and wave data

1.2 邊緣分布

首先，需要建立風(fēng)浪要素的邊緣分布模型。研究中發(fā)現(xiàn)，風(fēng)速、波高、波周期樣本具有單峰分布的特征（圖1（a）～（c）），采用常見的單峰分布模型進(jìn)行擬合，包括Weibull、廣義極值分布（GEV）、含有尺度參數(shù)和位置參數(shù)的t分布等。風(fēng)向和波向具有多峰分布的特征（圖1（d）、（e）），采用混合模型進(jìn)行擬合，包括混合Gaussian 分布、混合Gamma 分布、混合Weibull 分布。

圖1 風(fēng)浪要素概率分布直方圖及最優(yōu)邊緣概率分布曲線Fig.1 Histogram of wind and wave data and the optimal marginal probability distribution curve

1）Weibull 分布

式中：λ為尺度參數(shù)；k為形狀參數(shù)。

2）廣義極值分布（GEV）

式中：ξ、β、μ為分布函數(shù)參數(shù)；ξ為形態(tài)參數(shù)；σ為尺度參數(shù)；μ為位置參數(shù)。

3）含有尺度參數(shù)和位置參數(shù)的t分布

式中：ν、σ、μ 均為分布函數(shù)參 數(shù)，其 中：ν為形態(tài)參數(shù)；σ為尺度參數(shù)；μ為位置參數(shù)。

4）混合Gaussian 模型

5）混合Gamma 模型

當(dāng)采用上述概率分布模型對(duì)風(fēng)浪要素進(jìn)行擬合時(shí)，概率分布模型的參數(shù)估計(jì)采用極大似然法。另外，為了評(píng)價(jià)不同概率分布模型的擬合效果，采用AIC、BIC 和RMSE 對(duì)概率分布模型進(jìn)行擬合優(yōu)度評(píng)價(jià)，并據(jù)此選取最優(yōu)的概率分布模型。AIC、BIC 和RMSE 的計(jì)算 式為

式中：xi為樣本值；n為樣本數(shù)量；f(xi)為備選邊緣分布函數(shù)的密度函數(shù)；k為備選邊緣分布函數(shù)中分布參數(shù)的數(shù)量。

式中：xi為樣本值；n為樣本數(shù)量；f(xi)為備選邊緣分布函數(shù)的密度函數(shù)；k為備選邊緣分布函數(shù)中分布參數(shù)的數(shù)量。

式中：n為樣本數(shù)量；Pc為多維Copula 聯(lián)合分布理論頻率值。RMSE 值越小，擬合的效果越好。

通過觀察風(fēng)向概率直方圖（圖1（d））可以看出，風(fēng)向的概率分布有3 個(gè)較為明顯的峰值，因此，采用的混合概率模型中應(yīng)包含3 個(gè)單峰分布，故采用混合維度為3 的混合模型（3 個(gè)單峰分布函數(shù)混合）來擬合。與風(fēng)向類似，波向頻率直方圖（圖1（e））也有多個(gè)峰值，分別采取二維和三維混合模型對(duì)波向的概率分布進(jìn)行擬合，通過對(duì)比RMSE 發(fā)現(xiàn)，三維混合模型的概率擬合效果較好。表2 和表3 給出了風(fēng)速、波高、波浪周期、風(fēng)向、波向5 個(gè)風(fēng)浪要素的最優(yōu)邊緣概率分布類型和相應(yīng)參數(shù)。此外，圖1 還給出了這5 個(gè)風(fēng)浪要素樣本的直方圖及最優(yōu)邊緣概率分布曲線。由圖1 可以看出，選取的最優(yōu)邊緣概率分布曲線對(duì)5 個(gè)風(fēng)浪要素樣本的擬合效果較好。擬合優(yōu)度評(píng)價(jià)結(jié)果表明：風(fēng)速、波高、波周期的最優(yōu)擬合分布分別為Weibull 分布、廣義極值分布（GEV）、含有尺度參數(shù)和位置參數(shù)的t分布；而風(fēng)向和波向的最優(yōu)擬合分布均為混合Gaussian 分布。

表2 風(fēng)速、波高、波周期最優(yōu)邊緣概率分布Table 2 The optimal marginal probability distribution of wind speed,wave height and wave period

表3 風(fēng)向、波向邊緣分布Table 3 The marginal distribution of wind direction and wave direction

2 二維聯(lián)合分布

2.1 二維Copula 理論

根據(jù)Sklar 定理[17]，多維聯(lián)合分布和其邊緣分布可以寫為

式中：Fi(xi)為隨機(jī)變量xi的邊緣概率分布函數(shù)；C(·)為Copula 函數(shù)；θ為Copula 函數(shù)的參數(shù)。

將海洋環(huán)境變量定義為一個(gè)n維隨機(jī)變量X=（x1,x2…xi…xn），基于Copula 理論，聯(lián)合概率密度可以表示為

式中：Fi（xi）、f（xi）分別為隨機(jī)變量xi的邊緣概率分布的分布函數(shù)和概率密度函數(shù)；c為Copula 密度函數(shù)。

由式（10）可得二維聯(lián)合分布概率密度函數(shù)公式

式中：c為二維Copula 函數(shù)的密度函數(shù)。

由式（11）可知，為了研究風(fēng)浪要素之間的二維聯(lián)合概率分布，首先應(yīng)建立風(fēng)浪要素兩兩變量之間的Copula 函數(shù)。選取風(fēng)速—波高、波高—波周期、風(fēng)速—風(fēng)向、波高—波向4 個(gè)隨機(jī)變量對(duì)來研究并建立二維Copula 函數(shù)。目前，研究風(fēng)浪要素聯(lián)合分布的大多數(shù)文獻(xiàn)均將這些變量對(duì)作為研究對(duì)象。其次，這些變量對(duì)中的兩個(gè)變量之間往往存在較大的相關(guān)性，而這些相關(guān)性在進(jìn)行結(jié)構(gòu)（如海上風(fēng)機(jī)、跨海橋梁等）的動(dòng)力響應(yīng)分析時(shí)又至關(guān)重要。樣本數(shù)據(jù)經(jīng)過式（12）可以變換為范圍0～1 的標(biāo)準(zhǔn)分布，樣本的標(biāo)準(zhǔn)分布忽略了樣本數(shù)值大小的差異，只保留數(shù)據(jù)的相對(duì)大小關(guān)系，以便于更清晰地研究數(shù)據(jù)之間的相關(guān)性。圖2 是風(fēng)速—波高經(jīng)驗(yàn)Copula 頻率直方圖，可以看出，風(fēng)速和波高具有較強(qiáng)的尾部相關(guān)性。

圖2 風(fēng)速—波高標(biāo)準(zhǔn)分布二維頻率直方圖Fig.2 Binary frequency histogram of standard distribution of wind and wave

式中：U為將樣本變換為范圍為0～1的標(biāo)準(zhǔn)分布后的隨機(jī)變量；n為樣本個(gè)數(shù)；R為樣本點(diǎn)在所有樣本中的排序。

常用的二維Copula 函數(shù)類型有：Gaussian、T、Clayton、Gumbel、Frank，如表4 所示。

表4 5 種典型的Copula 函數(shù)Table 4 Five typical Copula functions

2.2 擬合優(yōu)度評(píng)價(jià)

要進(jìn)行擬合優(yōu)度評(píng)價(jià)，首先要計(jì)算經(jīng)驗(yàn)Copula，經(jīng)驗(yàn)Copula 可以通過式（13）計(jì)算[18]。

式中：n為樣本的大小，對(duì)每一個(gè)1≤i≤n，滿足u1i≤u1,u2i≤u2時(shí)，I(u1i≤u1,u2i≤u2)=1。

對(duì)應(yīng)三維的情況，經(jīng)驗(yàn)Copula 可通過式（14）計(jì)算。

式中：n為研究樣本的大小，對(duì)每一個(gè)1≤i≤n，滿足u1i≤u1,u2i≤u2,u3i≤u3時(shí)，I(u1i≤u1,u2i≤u2，u3i≤u3)=1。

為了對(duì)選取的Copula 函數(shù)進(jìn)行擬合優(yōu)度評(píng)價(jià)，選用均方根誤差（RMSE）作為評(píng)價(jià)標(biāo)準(zhǔn)來評(píng)價(jià)模型的優(yōu)劣。

Sn越小,說明擬合的Copula 模型與經(jīng)驗(yàn)Copula越接近，擬合效果越好。

建立二維聯(lián)合分布可以對(duì)多種Copula 函數(shù)分別進(jìn)行參數(shù)估計(jì)，采用貝葉斯框架和基于殘差的高斯似然函數(shù)進(jìn)行參數(shù)估計(jì)[19]。

貝葉斯分析已在多個(gè)領(lǐng)域應(yīng)用于參數(shù)估計(jì)。當(dāng)獲得新信息時(shí)，貝葉斯理論更新假設(shè)的先驗(yàn)概率，將所有建模的不確定因素歸因于參數(shù)，通過式（16）估計(jì)模型參數(shù)的后驗(yàn)分布。

在缺乏參數(shù)先驗(yàn)分布的有效信息時(shí)，可以采用均勻先驗(yàn)，假設(shè)殘差不相關(guān)，同方差、均值為零的高斯分布，那么似然函數(shù)就可以通過式（17）表示。

通過AIC、BIC 及RMSE 對(duì)不同類型的Copula模型進(jìn)行擬合優(yōu)度評(píng)價(jià)，擬合優(yōu)度評(píng)價(jià)結(jié)果見表5，參數(shù)估計(jì)值見表6。風(fēng)速—波高的最優(yōu)聯(lián)合分布為θ=0.846 1 的Gaussian Copula；波高—波浪周期的最優(yōu)聯(lián)合分布為θ=0.465 3 的Gaussian Copula；風(fēng)速—風(fēng)向的最優(yōu)聯(lián)合分布為θ=-1.145 7 的Frank Copula；波高—波向的最優(yōu)聯(lián)合分布為θ=0.541 0 的Gaussian Copula。圖3 為θ=0.846 1 的Gaussian Copula 函數(shù)圖像。

圖3 Gaussian Copula（θ=0.846 1）概率密度圖Fig.3 Probability density of Gaussian Copula (θ=0.846 1)

表5 二維Copula 擬合優(yōu)度評(píng)價(jià)Table 5 The performance evaluation of 2-dimensional Copula distribution

表6 最優(yōu)二維Copula 分布及擬合參數(shù)Table 6 The optimal 2-dimensional Copula distribution and fitting parameters

由圖2 和圖3 可以看出，Gaussian Copula 能較好地模擬風(fēng)速—波高的相關(guān)性。此外，樣本數(shù)據(jù)的二維頻率直方圖可以表示聯(lián)合分布的經(jīng)驗(yàn)圖像，與圖4 不同的是，圖3 只表示兩隨機(jī)變量之間的相關(guān)關(guān)系，略去了兩隨機(jī)變量之間數(shù)值大小的影響，通過圖3能更清晰地了解二者的相關(guān)關(guān)系。圖4 則表示出兩隨機(jī)變量之間聯(lián)合分布的圖像，保留了樣本數(shù)值大小的信息。通過二維頻率直方圖和聯(lián)合分布的概率密度圖可以直觀比較本文選取的二維Copula 模型的擬合效果，如圖4 所示。由圖4 可以看出，對(duì)于風(fēng)速—波高、波高—周期、風(fēng)速—波向，其二維聯(lián)合概率密度圖與經(jīng)驗(yàn)圖像的峰值和形狀都比較接近，說明相關(guān)性模擬較為合理，較好地考慮了風(fēng)浪要素之間的相關(guān)性。同時(shí)，從圖4 也可以看到波，高—波向的聯(lián)合分布與經(jīng)驗(yàn)直方圖的峰值大小存在一定差異，這是由于該樣本波高—波向之間的相關(guān)性比較復(fù)雜，采用單一的Copula函數(shù)模擬它們之間相關(guān)性的效果有限，在后續(xù)研究中，可以采用混合Copula 函數(shù)來進(jìn)行模擬，并優(yōu)化擬合效果。由于不同海域的風(fēng)浪要素之間的相關(guān)關(guān)系可能會(huì)有較大的差異，本研究側(cè)重于探究建立風(fēng)浪聯(lián)合分布的方法和思路，因此未對(duì)波高—波向的二維聯(lián)合概率模型作進(jìn)一步優(yōu)化。

圖4 風(fēng)浪要素兩兩變量之間的頻率直方圖及最優(yōu)二維Copula 函數(shù)Fig.4 Binary frequency histogram of wind and wave data and the corresponding optimal 2-dimensional Copula distribution

3 多維聯(lián)合分布

可以采用多維Copula 函數(shù)建立多維隨機(jī)變量的聯(lián)合分布，但可供選用的多維Copula 函數(shù)類型有限，且靈活性較弱。Joe[20]提出了采用Vine Copula來構(gòu)建多維隨機(jī)變量的聯(lián)合分布的方法，Vine Copula 具有結(jié)構(gòu)多樣靈活的優(yōu)點(diǎn)。

3.1 Pair-Copula 理論

Bedford 等[21]提出了通過Pair-Copula 構(gòu)建多維隨機(jī)變量聯(lián)合分布概率模型的方法，基于條件概率可以將多維隨機(jī)變量分解成一系列的Pair-Copula結(jié)構(gòu)，分解時(shí)采用的邏輯結(jié)構(gòu)不同，可以構(gòu)造出不同的模型。Pair-Copula 結(jié)構(gòu)推動(dòng)了Copula 理論在多維隨機(jī)變量應(yīng)用中的發(fā)展。

根據(jù)條件概率，多維聯(lián)合分布概率密度函數(shù)可以表示為

由式（19）可以得到二維隨機(jī)變量的聯(lián)合概率密度函數(shù)為

式中：a,j=1,2 …n，且a≠j；ca j(Fa(xa),Fj(xj))為xa和xj的二維Copula 密度函數(shù)。

由式（20）可以推導(dǎo)xj已知的情況下xa的條件概率密度

由式（21）可得在n維隨機(jī)變量u已知的條件下，任意隨機(jī)變量x的條件密度函數(shù)

式中：ua是n維隨機(jī)變量u中的一個(gè)分量；u-a是n維隨機(jī)變量u中去掉ua之后的n-1 維分量。

以三維聯(lián)合分布為例說明分解過程，根據(jù)式（19），三維的聯(lián)合分布函數(shù)可寫為

由式（22）可以得到

代入式（23）即可將多維聯(lián)合分布分解為Pair-Copula 結(jié)構(gòu)和邊緣分布的概率密度函數(shù)的乘積

3.2 Vine 結(jié)構(gòu)

高維Copula 函數(shù)可以通過R-Vine 結(jié)構(gòu)來建立?；趦蓛呻S機(jī)變量之間的相依組合，結(jié)合條件概率可以建立多維Copula 函數(shù)。R-Vine 中有兩類特殊的Vine：C-Vine 和D-Vine，這兩類Vine 有各自的邏輯結(jié)構(gòu)，用于建立高維變量之間的聯(lián)合分布。

C-Vine 結(jié)構(gòu)的特點(diǎn)是每層樹都有一個(gè)主節(jié)點(diǎn)，主節(jié)點(diǎn)連接其他所有節(jié)點(diǎn)。為方便闡述構(gòu)造原理，圖5 給出了四維的C-Vine 結(jié)構(gòu)圖，該C-Vine 的結(jié)構(gòu)有3 棵樹，每棵樹有一個(gè)主節(jié)點(diǎn)，主節(jié)點(diǎn)和其他節(jié)點(diǎn)（node）相連形成邊（edge），邊在下一棵樹中作為主節(jié)點(diǎn)，以此類推，直到連接所有樹的節(jié)點(diǎn)。對(duì)任意維n維C-Vine，有n-1 棵樹，有n(n-1)/2 條邊，每條邊對(duì)應(yīng)一個(gè)Pair-Copula 結(jié)構(gòu)。

圖5 四維C-Vine 分解結(jié)構(gòu)Fig.5 Decomposition structure for four-dimensional C-Vine model

D-Vine結(jié)構(gòu)的特點(diǎn)是每層樹的節(jié)點(diǎn)依次相連，呈直線狀。圖6給出了四維D-Vine的結(jié)構(gòu)，該D-Vine結(jié)構(gòu)有3棵樹，每棵樹中節(jié)點(diǎn)一次連接，形成邊，邊在下一棵樹中成為節(jié)點(diǎn)，用同樣的方式類推，直至連接所有節(jié)點(diǎn)。對(duì)任意維n維D-Vine，同樣有n-1棵樹，有n(n-1)/2條邊，每條邊對(duì)應(yīng)一個(gè)Pair-Copula結(jié)構(gòu)。

圖6 四維D-Vine 分解結(jié)構(gòu)Fig.6 Decomposition structure for four-dimensional D-Vine model

C-Vine 的多變量分解結(jié)構(gòu)的確定需要先確定根節(jié)點(diǎn)和其他節(jié)點(diǎn)的順序，C-Vine 的根節(jié)點(diǎn)通常選取與其他變量相關(guān)性最強(qiáng)的節(jié)點(diǎn)，可以采用Aas等[22]與Di?mann 等[23]提出的序慣估計(jì)法來確定。

1）從第1 層樹開始，計(jì)算所有隨機(jī)變量兩兩組合的經(jīng)驗(yàn)Kendall 相關(guān)系數(shù)τ，如表7 所示，選取和其他隨機(jī)變量的相關(guān)系數(shù)τ之和最大的隨機(jī)變量作為根節(jié)點(diǎn)，這樣就確定了第1 層樹的結(jié)構(gòu)。

表7 經(jīng)驗(yàn)Kendall’s τ 矩陣及Kendall’s τ 之和Table 7 The empirical Kendall’s τ and the sum of it

2）選擇第1 層樹中二維隨機(jī)變量的Copula 函數(shù)種類并使用極大似然法估計(jì)參數(shù)θ。

3）結(jié)合第2）步確定的Copula 函數(shù)及參數(shù)θ，通過式（25）將原隨機(jī)變量轉(zhuǎn)換為條件隨機(jī)變量，即樹2 的隨機(jī)變量。

式中：v-j為v中除去vj的n-1 維向量，F(xiàn)(x|v)為條件分布函數(shù)，Cxvj|v-j為連接F(x|v-j)與F(vj|v-j)的二維Copula 函數(shù)。

4）使用確定樹1 結(jié)構(gòu)的方法確定剩余的所有樹。

采用C-Vine Copula 構(gòu)建多維聯(lián)合分布，圖7 給出了最優(yōu)C-Vine 的構(gòu)造，1～5 分別表示風(fēng)速、波高、周期、風(fēng)向、波向，五維C-Vine Copula 有4 層樹，在第1 層樹中，隨機(jī)變量2（波高）與其余4 個(gè)變量分別通過二維Copula 連接，形成4 條邊；在第2 層樹中，這4 條邊作為節(jié)點(diǎn)，選取一個(gè)主節(jié)點(diǎn)，再一次通過二維Copula進(jìn)行連接，以此類推，直至連接所有節(jié)點(diǎn)。

圖7 維C-Vine Copula 結(jié)構(gòu)Fig.7 The optimal five-dimensional C-Vine model

由式（10）和圖7 可得，五維隨機(jī)變量(x1,x2,x3,x4,x5)的C-Vine 聯(lián)合概率密度函數(shù)如式（28）所示。

式中：Fi(·)為每個(gè)隨機(jī)變量的累積分布函數(shù)；F·|·(·|·)為條件分布函數(shù)；c·，·(·,·)為Copula 密度函數(shù)；c·|·(·|·)為Copula 密度函數(shù)。

可以通過比較AIC 值進(jìn)行C-Vine Copula 模型的擬合優(yōu)度評(píng)價(jià)，AIC 的值越小說明擬合效果越好。一般來說，C-Vine 和D-Vine 均可用于構(gòu)建數(shù)據(jù)的聯(lián)合分布，具體采用哪一種模型應(yīng)該取決于數(shù)據(jù)本身，通常應(yīng)通過試算來最終確定較優(yōu)的模型。分別采用C-Vine 和D-Vine 來構(gòu)建風(fēng)浪聯(lián)合分布模型，并對(duì)比兩種模型的AIC 值，發(fā)現(xiàn)C-Vine 能更好地模擬該海洋站點(diǎn)的風(fēng)浪聯(lián)合分布模型。因此，最終采用C-Vine 來構(gòu)建風(fēng)浪的聯(lián)合分布模型。此外，采用序慣估計(jì)法得到風(fēng)浪聯(lián)合概率分布的C-Vine 模型。通過序慣估計(jì)法得出的最優(yōu)根節(jié)點(diǎn)順序?yàn)?、1、4、5、3（1～5 分別表示風(fēng)速、波高、周期、風(fēng)向、波向），AIC 值為-27 508，如表8 所示。采用遍歷所有根節(jié)點(diǎn)順序的方法（即考慮了所有可能的根節(jié)點(diǎn)順序），對(duì)由序慣估計(jì)法得出的C-Vine 風(fēng)浪聯(lián)合概率分布模型進(jìn)行擬合優(yōu)度評(píng)價(jià)。五維隨機(jī)變量的根節(jié)點(diǎn)順序共有=60 種，針對(duì)每一種根節(jié)點(diǎn)順序，首先采用AIC 準(zhǔn)則建立相應(yīng)的C-Vine 模型，并通過AIC值的大小對(duì)60 個(gè)C-Vine 模型進(jìn)行排序，結(jié)果如表9所示。限于篇幅，僅給出了前20 個(gè)較優(yōu)的C-Vine模型及相應(yīng)的根節(jié)點(diǎn)順序。通過比較發(fā)現(xiàn)，通過序慣估計(jì)法建立的C-Vine 模型的AIC 值在60 種情況中排第3，并且其AIC 值與前面2 個(gè)模型的AIC 值非常接近（相差0.31%以內(nèi)）。因此，序慣估計(jì)法能高效且準(zhǔn)確的找出較優(yōu)的C-Vine 聯(lián)合概率模型，該模型能很好地描述風(fēng)浪要素間的相關(guān)關(guān)系。隨著變量維度的提升，根節(jié)點(diǎn)的組合會(huì)迅速增加，采用遍歷所有根節(jié)點(diǎn)順序的方法會(huì)極大地增加計(jì)算負(fù)擔(dān)，由此也體現(xiàn)出序慣估計(jì)法的優(yōu)越性。

表8 最優(yōu)C-Vine 參數(shù)及AIC 值Table 8 The optimal C-Vine parameters and AIC values

表9 不同C-Vine 模型的根節(jié)點(diǎn)順序和AIC 值Table 9 The root node order and AIC value of different C-Vine models

通過建立的C-Vine 模型仿真了風(fēng)速、波高、波浪周期、風(fēng)向和波向五維隨機(jī)變量之間的累積概率密度（CDF）。為了更直觀地觀察擬合效果，以風(fēng)速、波高、風(fēng)向?yàn)槔?，結(jié)合式（14），圖8 給出了其三維聯(lián)合分布的累積概率密度，黑色網(wǎng)格是原始數(shù)據(jù)的經(jīng)驗(yàn)Copula 圖像，彩色圖像是由C-Vine 模擬產(chǎn)生的圖像。由圖8 可知，建立的C-Vine 模型能較好地模擬風(fēng)速、波高、風(fēng)向的相關(guān)關(guān)系。

圖8 三維累積概率密度圖Fig.8 Diagram of three-dimensional cumulative probability density

4 結(jié)論

基于Vine Copula 函數(shù)研究了中國東海連云港海洋觀測站的風(fēng)浪要素之間的聯(lián)合概率分布，得出以下結(jié)論：

1）風(fēng)速、波高、波周期的概率分布為單峰分布，最優(yōu)擬合分布分別為Weibull 分布、廣義極值分布、t分布；而風(fēng)向和波向的概率分布為多峰分布，最優(yōu)擬合分布均為混合Gaussian 分布。

2）風(fēng)浪要素中兩兩隨機(jī)變量之間的聯(lián)合概率分布研究表明，風(fēng)速-波高、波高—波周期、風(fēng)速—風(fēng)向、波高—波向4 個(gè)二維變量對(duì)的最優(yōu)二維聯(lián)合概率分布分別為Gaussian Copula、Gaussian Copula、Frank Copula 和Gaussian Copula。

3）建立的C-Vine 模型可以較好地刻畫風(fēng)速、波高、波浪周期、風(fēng)向和波向五維隨機(jī)變量之間的聯(lián)合概率分布。

4）采用Vine Copula 函數(shù)建立了東海連云港海洋觀測站風(fēng)浪要素之間的聯(lián)合概率分布模型，對(duì)于中國其他海域海洋觀測站風(fēng)浪要素之間的聯(lián)合分布規(guī)律還有待進(jìn)一步研究。但該研究方法和思路可為中國其他海域海洋觀測站風(fēng)浪要素之間的聯(lián)合分布研究提供借鑒。

致謝：感謝中國國家科技資源共享服務(wù)平臺(tái)——國家海洋科學(xué)數(shù)據(jù)中心(http://mds.nmdis.org.cn/)提供數(shù)據(jù)支撐。