基于非對稱Archimedean Copula的三變量風浪重現水平分析

陳子燊，路劍飛，于吉濤

（1.中山大學水資源與環境系，廣東廣州510275；2.國土資源部廣州海洋地質調查局，廣東廣州510075；3.河南理工大學測繪與國土信息工程學院，河南焦作454000）

基于非對稱Archimedean Copula的三變量風浪重現水平分析

陳子燊1，路劍飛2，于吉濤3

（1.中山大學水資源與環境系，廣東廣州510275；2.國土資源部廣州海洋地質調查局，廣東廣州510075；3.河南理工大學測繪與國土信息工程學院，河南焦作454000）

采用非對稱Archimedean Copula函數與Kendall分布函數分析極端波況下的波高、周期和風速三變量聯合概率分布與風險率及其設計分位數，為海岸海洋工程設計和風險評估提供參考依據。以粵東汕尾海域的實測風浪數據為例，使用非對稱Gumbel-Hougaard Copula函數計算三變量風浪聯合分布的“或”重現期、“且”重現期和二次重現期及其最可能的風浪設計值。主要結論如下：對比不同設計風浪重現期顯示，“或”重現期的風險率偏高，“且”重現期的風險率偏低，二次重現期更準確地反映了特定設計頻率情況下三變量風浪的風險率；按目前有關規范設計要求的單變量風浪要素設計值已經達到安全標準，按三變量“或”重現期和三變量同頻率設計值推算的風浪設計值偏高，以最大可能概率推算的三變量風浪要素的二次重現期設計值可為相關工程安全與風險管理提供新的選擇。

非對稱Archimedean Copula；風浪風險評估；Kendall分布函數；二次重現期

極端風浪影響下的海岸和海洋工程存在著高損毀風險。極端海浪屬于相互關聯的多維隨機變量，如波高、波周期和相應風速。風速和波高相關性對設計荷載的確定有顯著影響（周道成等，2003），波高和周期的聯合分布可用來計算海洋結構物對波浪的響應（方鐘圣等，1989），實際海岸和海洋工程的損毀通常是其聯合作用的結果。鑒于海岸和海洋工程失事代價很高，探索多變量海浪的重現期，推算其聯合設計水平值，對于相關工程的設計與風險管理具有重要的應用參考價值。

描述多變量相關結構的Copula函數的提出大大地促進了多變量聯合分布及其風險概率在海岸海洋工程領域的應用。目前，大多采用可有效地描述兩變量聯合分布的阿基米德Copulas函數構建風速與波高或波高與周期聯合分布（董勝等，2011；秦振江等，2007；徐龍軍等，2013；陳子燊等，2011，2012）。但是，對于高維（3變量及以上）分布，由于各變量間的相依性是不對等的，一個參數的高維copulas不足以完整地描述變量間的相依性，這種情況下非對稱阿基米德Copulas（Asymmetric Archimedean Copulas）顯得更適用。由于求解相對復雜，目前非對稱阿基米德Copulas除在洪水頻率分析有少數研究外（Grimaldi et al，2006；Ganguli et al，2013；李天元等，2013），其它領域與方向研究很少。

至今，“或”和“且”重現期還是最常用的兩種多變量重現期定義方法，Salvadori等（2004）指出，在安全與危險事件的判定上兩者都存在著較大局限性。Salvadori等（2011）針對“或”重現期的不足提出了劃分安全與危險臨界域的新的多變量重現期——二次重現期（Secondary Return Periods，或稱Kendall重現期，Kendall Return Periods），其相關原理已在海岸工程設計研究中得到初步應用（Corbellaetal，2012；Salvadorietal，2013）。

目前，國內尚未見到基于Copula函數的三變量風浪聯合分布的研究。本文擬應用三維非對稱Archimedean極值copula，通過實例研究，推算極值波高及其相應波周期和風速聯合分布的重現期及其設計值，對比分析二次重現期（Secondary return periods）和“或”（“OR”）、“且”（“AND”）重現期及其危險率。此研究有助于深入認識極端風浪對海岸海洋工程的致災風險，有助于相關工程的風險管理，并為工程設計提供參考。

1 理論與方法

1.1 Copula函數

根據Sklar定理，若F（·）是一個n維隨機變量（X1，X2，…，Xn）的累積分布函數，其邊緣分布函數是連續函數F1，…，Fn，則有唯一的Copula函數C使得：

1.2 Archimedean copula極值分布函數

Archimedean copula是分析極端水文氣象事件最常用的copula函數族，其基本形式為：

式中，uj∈[0，1]（j＞1）為邊緣分布，φθ為Archimedean copula生成元，θ為參數。

Nelson（2006）和Salvadori等（2010）證明當且僅當邊緣分布和Copula函數均為極值分布時，構造的聯合分布才是極值分布，而Gumbel－Hougaard Copula是Archimedean copula函數族中的唯一多變量極值Copula函數，適用于極值事件的頻率分析。三變量對稱Archimedean Copula結構形式為：

三變量非對稱Archimedean Copula具有以下形式（Grimaldi et al，2006）：

式中，符號“o”表示函數組合。三種常用的非對稱Archimedean族Copula函數如下：

式中，θ1，θ2代表相關參數；u1，u2，u3表示邊緣累積分布。

1.3 三變量聯合分布重現期

假定E（u1，u2，u3）為三變量極端事件，對于預設閾值

“或”重現期和“且”重現期又統稱為首次重現期（Primary return periods）（Requena et al，2013）。然而，“或”重現期和“且”重現期存在對安全事件與危險事件識別的局限性，以二維分布定義以下兩類事件為危險事件：1）三變量中任一個超過預設閾值（記為事件E∨）；2）三變量同時超過預設閾值（記為事件E∧）。若將E∨作為危險事件，則危險率為P（U1＞若將E∧作為危險事件，則危險率為P

稱E∨的重現期為“或”重現期TO，三變量“或”重現期為：

稱E∧的重現期則為“且”重現期TA，三變量“且”重現期為：

（黃強等，2015）圖解說明如下：如圖1a所示，對于“或”重現期，A點處的重現期大于B點（p1＞p2），而A點識別的危險域卻有一部分在B點識別的危險域之外（深灰色部分），即當事件E*發生時，若以B點的組合值作為災害防御的標準則可以認為E*是安全的，造成的效應在防御的標準下；若以重現期更大的A點組合值作為災害防御的標準則E*顯然是危險事件。高風險的防御標準識別的安全事件卻被更低風險的防御標準識別為危險事件，這顯然是矛盾的。對于“且”重現期（圖1b），A點處的重現期小于B點（p1＜p2），但A點識別的危險域并沒有完全覆蓋B點的危險域，即當事件E*發生時，若以B點的組合值作為災害防御的標準，E*造成的效應超過了防御的標準，可認為E*是危險事件；若以A點的組合值作為災害防御的標準則E*造成的效應在防御的標準下，E*是安全事件。低風險的防御標準識別的危險事件在更高風險的防御標準里卻是安全事件，這顯然也是矛盾的。因此，無論是“或”重現期還是“且”重現期，都存在對安全事件與危險事件錯誤的識別，將造成工程設計與風險管理上的問題。

圖1 （a）“或”重現期和（b）“且”重現期對安全危險域識別的局限性

1.4 Kendall分布函數與二次重現期

鑒于“或”重現期和“且”重現期在安全事件與危險事件的判定上，兩者都存在著較大的局限性，Salvadori等（2010）利用Kendall分布函數劃分出亞臨界、臨界和超臨界三種情景。將三變量極端事件點構成的曲面分成三個部分：（1）次臨界域；（2）臨界面；（3）超臨界域。以等曲面作為極端事件風險的臨界面，任何落在次臨界域內的極端事件可視為特定臨界條件下的安全事件，落在超臨界域內的極端事件可視為特定臨界條件下的危險事件，落在臨界面上的極端事件則可視為警戒事件。在此基礎上提出了更精確的估計風險的多變量重現期——二次重現期，以區別于“或”重現期和“且”重現期這兩類首次重現期。

對給定的累積概率，必定存在多個極端事件與之對應，通過求累積概率小于或等于的概率可將多維的極端事件投射為一維分布。定義三維Archimedean copula的臨界概率條件下的重現期為：

稱TK為二次重現期，KC為Kendall分布函數。KC將多維的信息轉化為單維信息，因此利用KC可以方便地對多維的變量進行單維分析。對于三維copula函數，式（10）中的KC沒有顯式，可由蒙特卡洛模擬構建經驗Kendall函數作為替代進行計算：Kn（p）=m/n。式中，m為聯合累積概率小于或等于q的樣本個數，n為樣本總數.

1.5 多變量設計值

在單變量頻率分析中，對于給定的重現期可通過概率分布的反函數推算對應的設計分位值。對于多變量的情形，同一個重現期可以有無數的分位值組合與之對應，并且這些分位值組合不能通過概率分布的反函數直接估算。具有相同重現期Tp（C（u1，u2，u3）=p）的波高H、周期T與風速W分位值組合構成了一個三維點集（等曲面）。在這些不同的分位值組合中，必然存在一個組合（hml，tml，wml）使聯合概率密度f（h，t，w）達到最大值，即該組合出現的可能性最大。因此在設定的重現期下，出現可能性最大的H－T－W組合可作為海岸工程設計中的優選參考標準。式中，c為三維Archimedean Copula的概率密度函數，f（h）、f（t）和f（w）分別為波高、周期與風速的概率密度函數。

2 示例研究

2.1 研究海域與基本數據

研究數據采用粵東汕尾海域遮浪海洋站1972－1992年實測的歷年熱帶氣旋影響期間S向浪的最大波高H與相應的波周期T和風速W。期間最大波高H為7.3 m，相應波周期和風速分別為8 s和18.5 m/s。最大風速45 m/s，相應的波高和周期分別為6.2 m和6.8 s。

2.2 邊緣分布與聯合分布

3個樣本的邊緣分布均采用廣義極值（GEV）分布函數：

式中，ξ，β，μ分別為形態參數、尺度參數和位置參數。參數估計使用線性矩法。經驗頻率分布使用Gringorten公式計算。擬合結果采用均方根誤差（RMSE）、經驗頻率和理論頻率擬合誤差平方和（Q）和概率點據相關系數（PPCC）檢驗其擬合優度，結果見表1、圖2。顯然，年最大H和相應平均周期T與極值風速W均服從GEV分布。

表1 樣本的GEV分布參數與優度檢驗值

計算波高與相應周期與風速兩兩間的Kendall相關系數τ，波高與周期的相關系數為0.659；波高與風速的相關系數為0.304；風速與周期的相關系數為0.126。可見三變量之間存在著明顯的非對稱相關性。

構建非對稱Archimedean Gumbel－Hougaard copula、Frank copula和Clayton copula。首先計算波高與周期的Kendall相關系數τ，繼而求波高分布函數u1與周期分布函數u2的copula參數θ2和內copula函數C1，再利用C1和風速分布函數u3求得參數θ1和外copula函數C2。作為對比，采用離差平方和最小準則（OLS）的適線法估計對稱Copula聯合分布的參數θ（張雨等，2011）。構造的三變量風浪非對稱與對稱Copula聯合分布的擬合優度檢驗結果見表2、圖3。可見，此3種copula擬合結果都很好，同類copula中非對稱Copula擬合效果略優于對稱Copula。考慮到三維非對稱形式的極值copula更適用于極端風浪分析，因此選用非對稱Archimedean Gumbel-Hougaard建立H-T-W之間的三維聯合分布：

圖2 三變量邊緣分布圖

表2 Copula參數估計及擬合優度評價

圖33 種非對稱Archimedean copula的概率分布擬合

2.3 聯合分布重現期和危險率

設定重現期標準下H-T-W組合的3種重現期及其超值累積概率（P）計算結果見表3。由式（8）、（9）和（10）以及C的非遞減性可知，單變量重現期（T）、”或”重現期（TO）、”且”重現期（TA）和二次重現期（TK）四者的關系為：TO≤T≤TK≤TA，即，對于設定的重現期，二次重現期介于“或”重現期和“且”重現期之間，TK大于設定的重現期。重現期大則危險率小，反之則危險率大。由表3可見，對于設定的重現期，“或”重現期最小且小于設定的重現期標準，以任一風浪要素超標可能致災的“或”重現期為標準的危險率最大，“且”重現期最大且大于設定的重現期標準，按三變量風浪要素同時超標的“且”重現期設計海岸海洋工程時出現的危險率最小。然而由于“或”重現期和“且”重現期都存在確定危險域不準確的問題，如按此兩種首次重現期標準設計海岸海洋工程將存在對風險高估或低估的問題。從安全角度考慮，對三變量風浪任一要素超標致災的重現期標準宜采用大于“或”重現期和小于“且”重現期的的二次重現期更合理。以此極端風浪序列為例，預定100年的重現期標準，三變量風浪任一要素超標致災宜采用345年和危險率0.002 9的二次重現期標準，而非49年和危險率0.020 4的“或”重現期標準或風浪三要素同時超標可能致災的395年的“且”重現期和危險率0.002 5的標準。

表3 不同風浪變量不對稱H-T-W組合聯合分布重現期

2.4 風浪設計值

按設定重現期推算波高、周期和風速單變量設計值，以出現概率最大的原理推算三變量“或”重現期和二次重現期設計值分別列于表4。結果顯示，對于設定的5～100年重現期，按二次重現期推算的設計波高和設計周期分別小于“或”重現期設計值和邊緣分布設計值。按二次重現期推算的設計波高、設計周期和設計風速分別小于對應的邊緣分布設計值，設計波高、設計周期和設計風速的相對誤差為－14.9%～－19.1%、－3.4%～－5.3%和－10%～－14.5%。按“或”重現期推算的設計波高、設計周期和設計風速分別大于對應的邊緣分布設計值，二者的相對誤差可達8.5%～13.1%、3.4%～8.6%和13.8%～21.2%。這一結果也表明，按目前有關規范設計要求的單變量風浪設計值已可達到安全標準，按三變量“或”重現期推算的風浪設計值則過于保守將導致高估。盡管按三變量同時超閾值的“且”重現期標準推算的風浪設計值最低，但由于存在低估風險問題，如何更合理地估計同時超閾值條件下的風險與推算設計值有待作進一步分析。

表4 不同重現期H-T-W設計值

為了比較，采用非對稱Gumbel－Hougaard copula推算波高－周期－風速三變量同頻率分布聯合設計值：

3 結論

（1）非對稱Archimedean Copula能更好地反映汕尾海域三變量風浪之間的不對稱相關結構，其中作為非對稱極值copula更適用于分析此海域極端風浪重現水平。

（2）對比不同設計風浪重現期顯示，“或”重現期的風險概率偏高，“且”聯合重現期的風險概率偏低，根據臨界條件C（u1，u2，u3）=P和Kendall函數定義的二次重現期更準確地代表了特定設計頻率情況下風浪的風險率。

（3）按目前有關規范設計要求的單變量風浪設計值已可達到設計標準，依據出現概率最大的原則推算的三變量“或”重現期和三變量同頻率設計值推算的風浪設計值偏高，采用二次重現期設計值可為海岸海洋工程安全與風險設計提供新的選擇。

Corbella S,Stretch D D,2012.Multivariate return periods of sea storms for coastal erosion risk assessment.Natural Hazards and Earth System Sciences,12:2699-2708.

Ganguli P,Reddy M,2013.Probabilistic assessment of flood risks using trivariate copulas.Theor Appl Climatol 111:341-360.

Grimaldi S,Serincesco F,2006.Asymmetric copula in multivariate flood frequency analysis.Advances in Water Resources,29(8):150-164.

Nelson R B,2006.An introduction to copulas(Second edition).New York:Springer.

Requena A I,Mediero L,Garrote L,2013.Bivariate return period based on copulas for hydrologic dam design:comparison of theoretical and empirical approach.Hydrol.Earth Syst.Sci.Discuss.,10,557-596.

Salvadori G,De Michele C,2004.Frequency analysis via copulas:Theoretical aspects and applications to hydrological events.Water Resources Research,40,W12511,doi:10.1029/2004WR003133.

Salvadori G,De Michele C,2010.Multivariate multiparameter extreme value models and return periods:A copula approach.Water Resource Research,46,W10501,doi:10.1029/2009WR009040.

Salvadori G,Michele D C,Durante F,2011.On the return period and design in a multivariate framework.Hydrology and Earth System Sciences,15:3293-3305.

Salvadori G,Tomasicchio G R,Alessandro F D,2013.Multivariate approach to design coastal and off-shore structures.Journal of Coastal Research,65:386-391.

陳子燊,曹深西,2012.基于Copula函數的波高與周期長期聯合分布.海洋通報,31(6):8-13.

陳子燊,2011.波高與風速聯合概率分布研究.海洋通報,30(2):158-163.

董勝,周沖,陶山山,等,2011.基于Clayton copula函數的二維Gumbel模型及其在海洋平臺設計中的應用.中國海洋大學學報,41(10):117-120.

方鐘圣,戴順孫,金承儀,1989.海洋特征波高和周期的長期聯合分布及其應用.海洋學報,11(5):535-543.

黃強,陳子燊,2015.基于二次重現期的多變量洪水風險評估.湖泊科學,27(2):352-360.

李天元,郭生練,閆寶偉,等,2013.基于多變量聯合分布推求設計洪水過程線的新方法.水力發電學報,32(3):10-15.

秦振江,孫廣華,閆同新,等,2007.基于Copula函數的聯合概率法在海洋工程中的應用.海洋預報,24(2):83-90.

徐龍軍,陳祉宏,周道成,等,2013.基于Archimedean Copula模型的風浪聯合分布第二重現期.天津大學學報(自然科學與工程技術版),46(2):114-120.

張雨,宋松柏,2011.基于Archimedean Copula的三維干旱特征變量聯合分布研究.中國農村水利水電,1:65-68.

周道成,段忠東,2003.耿貝爾邏輯模型在極值風速和有效波高聯合概率分布中的應用.海洋工程,21(2):45-51.

Analysis on return levels of trivariate stormy waves based on asymmetric Archimedean copula function

CHEN Zi-shen1,LU Jian-fei2,YU Ji-tao3

（1.Department of Water Resource and Environment,Sun Yat-sen University,Guangzhou 510275,China;2.Guangzhou Marine Geological Survey,Guangzhou 510075,China;3.School of Surveying and Land Information Engineering,Henan Polytechnic University,Jiaozuo 454000,China）

The joint probability distribution for triple variables（i.e.,wave height,wave period and wind speed under extreme wave condition）was analyzed by using 3-dimensional asymmetric Archimedean copula functions.And the trivariate stormy wave risk probabilities were further explored based on the Kendall distribution function in order to provide the basis for the coastal marine engineering design and risk assessment.A case is studied by using annual maximum wave heights and the relevant wave periods and wind speeds measured in Shanwei sea waters of eastern Guangdong,the primary return periods and secondary return periods of trivariate stormy wave joint distribution and the most likely design values were computed by using the asymmetric Gumbel-Hougaard copula.The main conclusions of this study can be summarized as follows：Comparing the risk probabilities of trivariate stormy waves among the different design return periods,the'OR'joint return period showed higher risk probabilities,and the risk probabilities of'AND'joint return period were lower,while the secondary return periods were more accurately depicted the risk probabilities under these specific design frequencies.According to the relevant specifications of the current design requirements,the univariate wave design values have reached safety standards.The estimated storm wave design values with the trivariate'OR'return period and triple variables with the same frequencies were obviously higher than expected.The most-likely design realization of the secondary return period can serve as the new selec-tion for coastal marine engineering and risk management.

asymmetric Archimedean copulas;risk assessment of stormy waves;Kendall distribution function;secondary return periods

P429

A

1001－6932（2017）06－0631－07

10.11840/j.issn.1001-6392.2017.06.004

2017-03-27；

2017-06-13

國家自然科學基金（41371498）。

陳子燊(1952-)，教授，從事工程水文與極端水文氣象事件與風險研究。電子郵箱:eesczs@mail.sysu.edu.cn。

袁澤軼）