基于Copula函數(shù)的鐵山港雨潮遭遇組合分析

甘富萬，邵帥杰，黃宇明，胡秀英，劉 旻，陳 帥

（1.廣西大學土木建筑工程學院，南寧 530004；2.廣西防災減災與工程安全重點實驗室，南寧 530004；3.南京水利科學研究院港口航道泥沙工程交通行業(yè)重點實驗室，南京 210029；4.南寧學院，南寧 530004）

0 引 言

洪、潮組合問題一直是沿海地區(qū)研究的熱點問題。近年來，Copula 函數(shù)廣泛應用于水文氣象領(lǐng)域，對于具有任意邊緣分布、任意相關(guān)性結(jié)構(gòu)的多個變量，Copula 函數(shù)可以方便地構(gòu)造出其聯(lián)合概率分布，并用于分析多變量水文事件。涂新軍等［1］基于Copula方法模擬24 h暴雨遭遇日高潮位的聯(lián)合分布特征，對比雨潮遭遇傳統(tǒng)重現(xiàn)期和二次重現(xiàn)期差異，根據(jù)同頻法和權(quán)函數(shù)法反推計算雨潮設(shè)計組合值，為濱海地區(qū)防洪潮設(shè)計提供參考。楊星等［2］采用Clayton Copula 構(gòu)建了深圳市洪潮遭遇組合的風險分析模型，得到了深圳市洪潮組合的風險率，可用于深圳地區(qū)河道防洪治理。劉曾美等［3］采用Copula函數(shù)構(gòu)建感潮河段的年最大洪水流量和相應潮位的聯(lián)合分布以及年最高潮位和相應洪水流量的聯(lián)合分布，可為感潮河段洪潮遭遇組合的合理選取提供科學依據(jù)。

鑒于此，本文基于Copula 函數(shù)，構(gòu)建鐵山港海灣暴雨與石頭埠站潮位的聯(lián)合概率分布，引入二次重新期概念，在雨潮組合風險分析的基礎(chǔ)上，為鐵山港海灣防洪治澇規(guī)劃提供科學依據(jù)。

1 雨潮遭遇組合分析方法

1.1 Copula函數(shù)基本理論

Copula 函數(shù)是定義域為［0，1］均勻分布的多維聯(lián)合分布函數(shù)，它可以通過連接多個隨機變量的邊緣分布來構(gòu)造聯(lián)合分布。令H為聯(lián)合分布函數(shù)，F(xiàn)和G為其邊緣累積分布函數(shù)，那么存在唯一的Copula函數(shù)C使得對?x，y∈有：

如果F和G是連續(xù)的，則C是唯一的。否則，C在RanF×RanG（Ran 表示值域）唯一確定。相反地，如果C是一個Copula函數(shù)，F(xiàn)和G是邊緣累積分布函數(shù)，則上式定義H 是具有邊緣累積分布F和G的一個二維聯(lián)合累積分布函數(shù)。

現(xiàn)階段水文氣象領(lǐng)域的Copula 函數(shù)主要有橢圓Copula 函數(shù)［4］、Plackett Copula 函數(shù)［5］以及阿基米德Copula 函數(shù)［6，7］，相較于其他兩種，阿基米德Copula 函數(shù)具有構(gòu)造簡單且相關(guān)參數(shù)θ計算較為簡便等特點，使得阿基米德Copula 函數(shù)應用最為廣泛［8，9］。本文選取阿基米德Copula 函數(shù)中應用最為廣泛的Frank、Clayton、GH、AMH Copula、No.16 Copula 以及No.17 Copula 函數(shù)為候選函數(shù)，并對其函數(shù)模型進行檢驗，選擇最合適的Copula函數(shù)模型。

1.2 雨潮遭遇設(shè)計組合

在流域防洪規(guī)劃中，通常需要確定雨潮組合方式，即在給定流域暴雨的條件下，確定河口潮位取值情況。但在實際工程中，當流域暴雨一定時，河口潮位取值并非唯一，并且在實際計算中不同的設(shè)計組合值所形成的防雨、潮能力也不盡相同，在現(xiàn)階段針對多變量水文事件無法確定設(shè)計組合的情況，一般是取特殊的設(shè)計組合如同頻組合、最可能組合等進行分析計算。

雨、潮組合的同頻組合是指在計算中令雨、潮的遭遇事件同頻發(fā)生，即令其邊緣分布u=v；另外一種組合就是最可能組合，其表達式如下，

式中：c（u，v）為雨潮聯(lián)合分布的概率密度函數(shù)；f（x）和f（y）則是其邊緣分布的概率密度函數(shù)。

1.3 雨潮組合風險

為了定量分析不同的雨潮組合發(fā)生的概率，本文構(gòu)建兩種不同的風險概率模型，用于評價洪潮遭遇風險：①同現(xiàn)風險率。即降雨和潮位同時超過某一設(shè)計標準的概率，該情況是最不利，危險最大的情況。②條件風險率。即當流域發(fā)生超過某一標準暴雨時，潮位也發(fā)生超過某一標準潮位的概率。設(shè)P表示流域暴雨，Z表示潮位，借助Copula 函數(shù)，同現(xiàn)風險率和條件風險率的計算公式分別為：

式中：u為P的邊緣分布函數(shù)；v為Z的邊緣分布函數(shù)。

2 實例分析

鐵山港地區(qū)地勢較為低洼，容易受到暴雨洪水以及風暴潮的影響。以鐵山港工業(yè)區(qū)為對象的研究區(qū)域的匯水面積為19.41 km2，屬于小流域匯水區(qū)域?，F(xiàn)階段，針對雨、潮組合的水文事件研究選樣一般有兩種方法：①以雨為主的24 h 暴雨遭遇日高潮位；②以潮為主的日高潮位遭遇日降水。參考年最大值法［1］，選擇北海站1967-2008年的年最大24 h 降雨序列作為本次研究的雨量資料，再根據(jù)北海站年最大24 h 降雨事件，選擇對應時段石頭埠站石頭埠站1967-2008年的42年日最高潮位序列，作為雨、潮事件進行分析計算。

2.1 雨潮聯(lián)合分布構(gòu)建

根據(jù)北海站連續(xù)42 a 最大24 h 雨資料序列，以及石頭埠潮位站對應的日最高潮位的Cs和Cv值，計算雨量資料以及潮位資料各自的四種邊緣分布函數(shù)：皮爾遜III 型（PE3）、廣義極值分布（GEV）、廣義Pareto 分布（GPD）和廣義Logistic 分布（GLO）的相關(guān)參數(shù)，并采用RMSE 準則和AIC 準則對雨量、潮位的邊緣分布函數(shù)進行優(yōu)選。計算結(jié)果表明雨量資料序列的最優(yōu)邊緣分布模型為皮爾遜III 型（PE3），而潮位資料序列的最優(yōu)邊緣分布模型則為廣義Pareto 分布（GPD）函數(shù)模型，因此分別選取皮爾遜III 型和廣義Pareto 分布雨量以及相應潮位的邊緣分布模型，并進行下一步計算。

采用Kendall 秩相關(guān)系數(shù)方法來計算不同Copula 函數(shù)模型的相關(guān)參數(shù)得到北海站歷年最大24 h 降雨量與石頭埠站相應的潮位資料之間的Kendall秩相關(guān)系數(shù)τ的值為0.217 2，證明雨量資料與潮位資料之間存在一定的相關(guān)性。

計算Frank Copula、Clayton Copula、GH Copula、AMH Copula、No.16 Copula 以及No.17 Copula 的相關(guān)參數(shù)，計算結(jié)果見表2。

表2 不同Copula函數(shù)的參數(shù)計算結(jié)果Tab.2 The result of the parameter calculation of different Copula functions

通過對比不同Copula 函數(shù)參數(shù)的取值范圍，得到上述6 種Copula 函數(shù)的參數(shù)計算結(jié)果均在合理的范圍，可以進行下一步的計算。

選擇Frank、Clayton、GH、AMH Copula、No.16 Copula 以及No.17 Copula 函數(shù)作為候選函數(shù)，采用圖形評價法，AIC 準則［10，11］，BIC準則［12］和OLS準則［13］進行函數(shù)模型評價。

表1 雨、潮邊緣分布計算成果Tab.1 Calculated results of rain and tide margin distribution

圖形評價法的計算結(jié)果如圖1。

圖1 雨量與潮位聯(lián)合分布理論頻率與經(jīng)驗頻率擬合結(jié)果Fig.1 Fitting results of theoretical frequency and empirical frequency of joint distribution of rainfall and tide level

由圖1可以看出不同的聯(lián)合分布模型均比較分散地分布在45°線的兩側(cè)，同時根據(jù)理論頻率與經(jīng)驗頻率的散點分布情況，并不能直觀明顯地判斷出哪個聯(lián)合分布模型是最優(yōu)的雨、潮聯(lián)合模型，因此需要通過分析其他優(yōu)度評價結(jié)果才可以選擇出最合適Copula函數(shù)模型。

表3計算結(jié)果顯示，6 種不同的雨、潮聯(lián)合分布模型的分別對應的三種函數(shù)優(yōu)度評價計算結(jié)果差異不大，均比較接近。通過對比，可以發(fā)現(xiàn)GH Copula函數(shù)的AIC值，BIC值以及RMSE值均為6 個Copula 函數(shù)中的最小值，根據(jù)函數(shù)優(yōu)度評價方法——評價準則計算結(jié)果越小，函數(shù)擬合優(yōu)度越好的依據(jù)，證明了GH Copula 函數(shù)是本研究中雨、潮聯(lián)合分布的最優(yōu)模型。根據(jù)上述結(jié)果，雨、潮的聯(lián)合分布模型為：

表3 Copula函數(shù)聯(lián)合分布模型優(yōu)度評價結(jié)果Tab.3 Results of the evaluation of the superiority of the joint distribution model of Copula function

根據(jù)雨、潮聯(lián)合分布計算結(jié)果，采用最優(yōu)的GH Copula 函數(shù)作為雨、潮聯(lián)合分布模型，并分析雨、潮聯(lián)合重現(xiàn)期以及同現(xiàn)重現(xiàn)期下的同頻組合以及最可能組合等的設(shè)計組合值。

分別采用變量X和Y來表示北海站的雨量資料和石頭埠站潮位資料序列，以GH Copula 函數(shù)構(gòu)建雨、潮組合的聯(lián)合分布模型如圖2。

圖2 雨量與潮位資料聯(lián)合概率分布模型Fig.2 Joint probability distribution model of rainfall and tide level data

2.2 雨潮組合設(shè)計重現(xiàn)期

現(xiàn)階段對于雨、潮組合的重現(xiàn)期的研究關(guān)注較多的是同現(xiàn)重現(xiàn)期以及二次重現(xiàn)期。在單變量條件下，雨量以及潮位超過某一個閾值的計算可以表示如下：

雨、潮同現(xiàn)重現(xiàn)期（“And”重現(xiàn)期）表示為降雨量以及潮位同時超過某一個特定設(shè)計閾值的情況。聯(lián)合Copula 函數(shù)，則根據(jù)水文重現(xiàn)期概率計算方法，雨、潮組合的同現(xiàn)重現(xiàn)期可以表示如下：

雨、潮聯(lián)合重現(xiàn)期（“OR”重現(xiàn)期）表示為降雨量或潮位其中一個超過某一個特定設(shè)計閾值的情況。聯(lián)合Copula 函數(shù)，則根據(jù)水文重現(xiàn)期概率計算方法，雨、潮組合的聯(lián)合重現(xiàn)期可以表示如下：

近年來一些專家學者認為傳統(tǒng)重現(xiàn)期并不是很好地描述多變量危險事件，于是提出了二次重現(xiàn)期［14-17］。生存Kendall重現(xiàn)期［14］是在同現(xiàn)重現(xiàn)期的基礎(chǔ)上，以雨、潮聯(lián)合分布概率相等作為的一系列事件作為臨界事件，并以此來劃分雨、潮聯(lián)合分布的危險域和安全域。一般情況下，采用基于C（u，v）的Kendall分布函數(shù)KC作為生存Kendall 重現(xiàn)期對應的超閾值的聯(lián)合分布概率函數(shù)［15］，因此可以表示為：

式中，q∈(0，1)。KC的分布函數(shù)如下：

則對于雨、潮組合來說，其生存Kendall重現(xiàn)期可以表示為：

式中：NE是暴雨事件年平均發(fā)生次數(shù)，在本文研究中，取歷年最大24 h 降雨量及其相對時刻的潮位作為雨、潮序列組合，因此NE取值為1。

同時，針對聯(lián)合（“OR”）重現(xiàn)期在描述多變量水文事件時存在局限，Salvadori［16，17］等在聯(lián)合重現(xiàn)期的基礎(chǔ)上提出采用Kendall重現(xiàn)期來計算多變量水文研究。即：

式中：φ(p)是阿基米德Copula 函數(shù)的生成元；φ′(p)是生成元函數(shù)的導數(shù)。

Kendall重現(xiàn)期可以定義如下，

2.3 雨潮遭遇設(shè)計組合分析

對傳統(tǒng)重現(xiàn)期和二次重現(xiàn)期的雨、潮同頻組合以及最可能組合進行計算，從表4的計算結(jié)果可以看出，同頻組合條件下或者最可能組合的二次重現(xiàn)期對應的100、50、20、10、5以及2 a重現(xiàn)期水平下的雨、潮組合均明顯大于傳統(tǒng)重現(xiàn)期所對應的各重現(xiàn)期下的組合。同現(xiàn)重現(xiàn)期條件下，最可能組合的各設(shè)計重現(xiàn)期水平年對應的潮位的重現(xiàn)期范圍為91.97～1.60 a，而同頻組合下的相應的重現(xiàn)期范圍是30.03～1.47 a，且最可能組合對應的各重現(xiàn)期條件下的潮位值均大于同頻組合的設(shè)計潮位值；同時，可以看到最可能組合條件的各重現(xiàn)期水平年對應的雨量設(shè)計的重現(xiàn)期范圍為1.85～1.34 a，均明顯小于同頻組合下雨量設(shè)計重現(xiàn)期。生存Kendall重現(xiàn)期條件下，最可能組合的各設(shè)計重現(xiàn)期水平年對應的潮位的重現(xiàn)期范圍為316.26～3.06 a，而同頻組合下的相應的重現(xiàn)期范圍是95.24～2.27 a，生存Kendall重現(xiàn)期條件下不同設(shè)計水平年的最可能組合的潮位設(shè)計值均明顯大于同頻組合條件的潮位各水平年設(shè)計值。由此可以推斷在相同的雨量、潮位組合對應下的同現(xiàn)重現(xiàn)期大于生存Kendall重現(xiàn)期，也就是說在相同組合條件下，生存Kendall 重現(xiàn)期要比同現(xiàn)重現(xiàn)期相對來說更加安全。

表4 雨、潮遭遇組合不同重現(xiàn)期設(shè)計值Tab.4 Different reproducibility design values for rain and tide encounter combinations

聯(lián)合重現(xiàn)期同頻組合的計算重現(xiàn)期為172.41～3.02 a，相同條件下Kendall 重現(xiàn)期的變化范圍為75.19～1.67 a。相同的設(shè)計重現(xiàn)期下，聯(lián)合重現(xiàn)期的各設(shè)計值均比Kendall 重現(xiàn)期的雨、潮設(shè)計組合要大。即相同的雨、潮設(shè)計組合，采用Kendall 重現(xiàn)期的設(shè)計標準要比聯(lián)合重現(xiàn)期高，說明Kendall重現(xiàn)期比聯(lián)合重現(xiàn)期更加安全。

進一步計算給定設(shè)計雨量條件下的各設(shè)計水平年條件下的條件風險概率以及同現(xiàn)風險概率，如表5所示。

表5 雨、潮組合同現(xiàn)風險以及條件風險概率Tab.5 Rain，tide combination co occurrence risk and conditional risk probability

表5表明，雨、潮不同設(shè)計重現(xiàn)期水平年組合的同現(xiàn)風險概率均比較小，且雨量（或潮位）遭遇重現(xiàn)期越大，遭遇潮位（或雨量）越大的同現(xiàn)風險概率越??；而當雨、潮遭遇均較小的時，雨潮同現(xiàn)風險概率則相對較大，這說明雨、潮遭遇組合在重現(xiàn)期較大的水平下，其相關(guān)性不高，但當雨、潮遭遇水平較小的情況下，其相關(guān)性相對較大。從條件風險概率來分析，也可以看出當雨量設(shè)計重現(xiàn)期較大時，遭遇較大重現(xiàn)期水平的潮位也較大時的條件概率較小，其中當雨量設(shè)計重現(xiàn)期為100 a 時，潮位遭遇也為100 a的條件風險概率為0.286；在相同條件下，潮位遭遇為2 a時的條件風險概率則為0.879。綜上，可以得到雨、潮組合各自遭遇重現(xiàn)期越小，其同現(xiàn)風險概率越高；而當雨量遭遇重現(xiàn)期越大的條件下，其遭遇潮位重現(xiàn)期越小的事件的條件概率風險越大。

3 結(jié) 論

（1）通過計算雨量資料以及潮位資料各自的四種邊緣分布函數(shù)的相關(guān)參數(shù)，并采用RMSE 準則和AIC 準則對雨量、潮位的邊緣分布函數(shù)進行優(yōu)選，最終選定雨量和潮位資料序列最優(yōu)邊緣分布函數(shù)分別為皮爾遜III 型（PE3）和廣義Pareto 分布（GPD）函數(shù)模型。

（2）在圖形評價法的基礎(chǔ)上，通過AIC值、BIC值、RMSE值，結(jié)合函數(shù)優(yōu)度評價方法，綜合考慮，最終選定GH Copula 函數(shù)構(gòu)建鐵山港區(qū)域雨潮聯(lián)合分布。

（3）通過分析雨、潮組合的傳統(tǒng)重現(xiàn)期、二次重現(xiàn)期以及同頻組合，最可能組合等可以看出，二次重現(xiàn)期相對于傳統(tǒng)重現(xiàn)期可以更好地描述雨、潮組合的設(shè)計情況，采用二次重現(xiàn)期能夠進行防洪減災設(shè)計會更加安全、合理。

（4）分別計算了不同暴雨潮位組合下的同現(xiàn)風險率和條件風險率，分析了各風險率的變化規(guī)律，為鐵山港區(qū)域防洪治澇規(guī)劃提供科學的依據(jù)。