基于Vine Copula函數的風浪要素聯合概率分布模型

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

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

海洋環境惡劣，往往存在大風和巨浪的組合作用，同時，跨海大橋具有基礎結構尺寸大、主跨輕柔、阻尼小、剛度小等特點，導致跨海大橋對風浪等海洋環境荷載非常敏感，極端風浪荷載成為威脅跨海大橋安全的主要因素之一。目前，跨海大橋多災害作用的研究中通常研究單一災害（如風、波浪、地震、沖刷等）對橋梁結構的作用，未能很好地考慮各災害之間的相關性。特別是風浪要素，其包含了眾多影響跨海大橋動力特性的參數，如風速、風向、波高、波向、波浪周期等，各參數是隨時間變化的隨機變量，且相互之間具有不同的尾部相關結構（即非線性）。因此，構造風浪要素的聯合分布函數是研究海洋環境中風浪耦合效應及進一步探究風浪各要素對跨海大橋動力響應特征影響規律的基本前提。

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

1.2 邊緣分布

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

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

1）Weibull 分布

式中：λ為尺度參數；k為形狀參數。

2）廣義極值分布（GEV）

式中：ξ、β、μ為分布函數參數；ξ為形態參數；σ為尺度參數；μ為位置參數。

3）含有尺度參數和位置參數的t分布

式中：ν、σ、μ 均為分布函數參 數，其 中：ν為形態參數；σ為尺度參數；μ為位置參數。

4）混合Gaussian 模型

5）混合Gamma 模型

當采用上述概率分布模型對風浪要素進行擬合時，概率分布模型的參數估計采用極大似然法。另外，為了評價不同概率分布模型的擬合效果，采用AIC、BIC 和RMSE 對概率分布模型進行擬合優度評價，并據此選取最優的概率分布模型。AIC、BIC 和RMSE 的計算 式為

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

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

式中：n為樣本數量；Pc為多維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

表3 風向、波向邊緣分布Table 3 The marginal distribution of wind direction and wave direction

2 二維聯合分布

2.1 二維Copula 理論

根據Sklar 定理[17]，多維聯合分布和其邊緣分布可以寫為

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

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

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

由式（10）可得二維聯合分布概率密度函數公式

式中：c為二維Copula 函數的密度函數。

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

圖2 風速—波高標準分布二維頻率直方圖Fig.2 Binary frequency histogram of standard distribution of wind and wave

式中：U為將樣本變換為范圍為0～1的標準分布后的隨機變量；n為樣本個數；R為樣本點在所有樣本中的排序。

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

表4 5 種典型的Copula 函數Table 4 Five typical Copula functions

2.2 擬合優度評價

要進行擬合優度評價，首先要計算經驗Copula，經驗Copula 可以通過式（13）計算[18]。

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

對應三維的情況，經驗Copula 可通過式（14）計算。

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

為了對選取的Copula 函數進行擬合優度評價，選用均方根誤差（RMSE）作為評價標準來評價模型的優劣。

Sn越小,說明擬合的Copula 模型與經驗Copula越接近，擬合效果越好。

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

貝葉斯分析已在多個領域應用于參數估計。當獲得新信息時，貝葉斯理論更新假設的先驗概率，將所有建模的不確定因素歸因于參數，通過式（16）估計模型參數的后驗分布。

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

通過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 函數圖像。

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

表5 二維Copula 擬合優度評價Table 5 The performance evaluation of 2-dimensional Copula distribution

表6 最優二維Copula 分布及擬合參數Table 6 The optimal 2-dimensional Copula distribution and fitting parameters

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

圖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 理論在多維隨機變量應用中的發展。

根據條件概率，多維聯合分布概率密度函數可以表示為

由式（19）可以得到二維隨機變量的聯合概率密度函數為

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

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

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

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

以三維聯合分布為例說明分解過程，根據式（19），三維的聯合分布函數可寫為

由式（22）可以得到

代入式（23）即可將多維聯合分布分解為Pair-Copula 結構和邊緣分布的概率密度函數的乘積

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-dimensional C-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-dimensional D-Vine model

C-Vine 的多變量分解結構的確定需要先確定根節點和其他節點的順序，C-Vine 的根節點通常選取與其他變量相關性最強的節點，可以采用Aas等[22]與Di?mann 等[23]提出的序慣估計法來確定。

1）從第1 層樹開始，計算所有隨機變量兩兩組合的經驗Kendall 相關系數τ，如表7 所示，選取和其他隨機變量的相關系數τ之和最大的隨機變量作為根節點，這樣就確定了第1 層樹的結構。

表7 經驗Kendall’s τ 矩陣及Kendall’s τ 之和Table 7 The empirical Kendall’s τ and the sum of it

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

3）結合第2）步確定的Copula 函數及參數θ，通過式（25）將原隨機變量轉換為條件隨機變量，即樹2 的隨機變量。

式中：v-j為v中除去vj的n-1 維向量，F(x|v)為條件分布函數，Cxvj|v-j為連接F(x|v-j)與F(vj|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 可得，五維隨機變量(x1,x2,x3,x4,x5)的C-Vine 聯合概率密度函數如式（28）所示。

式中：Fi(·)為每個隨機變量的累積分布函數；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 風浪聯合概率分布模型進行擬合優度評價。五維隨機變量的根節點順序共有=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

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

通過建立的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/)提供數據支撐。