基于Copula函數(shù)的甬江流域設計潮位過程研究

薛曉鵬，張衛(wèi)國，朱從飛，趙思遠，郝振純

(1.寧波市水利水電規(guī)劃設計研究院，浙江 寧波 315192；2.河海大學水文水資源與水利工程科學國家重點實驗室，南京 210098)

甬江流域地處寧波市，位于我國東海岸，為獨流入海的小流域，流域內(nèi)最大的排洪河道——甬江在鎮(zhèn)?？谥苯优c外海相連，城市防洪易受潮汐影響。鎮(zhèn)海站設計潮位過程是影響寧波市城市河道堤防設計的主要因素，因此研究鎮(zhèn)海站設計潮位對寧波市城市防洪具有重要意義。目前工程設計中多采用設計高潮位與設計潮差同頻率放大的方法，選取典型潮型，然后將高潮位和潮差按照設計值縮放典型潮型。但是這種方法沒有考慮設計高潮位與潮差的遭遇可能性大小。

研究不同水文系列遭遇頻率常用多維聯(lián)合分布的方法。目前，多維聯(lián)合分布模型的研究和應用已經(jīng)比較成熟，Copula函數(shù)因形式多樣且易于求解，廣泛地應用在水文系列的多維聯(lián)合分析中。比如，楊志勇等[1]采用Copula函數(shù)擬合降水距平百分率序列，計算出灤河流域各站點旱澇組合事件發(fā)生的概率；張冬冬等[2]采用Copula函數(shù)探討洪水多要素的聯(lián)合概率分布和條件概率分布；李天元等[3]構(gòu)造了分別金沙江和岷江2江洪峰流量的單參數(shù)和雙參數(shù)Copula聯(lián)合分布，表明雙參數(shù)Copula函數(shù)擬合效果優(yōu)于單參數(shù)Copula函數(shù)。黃鋒華等[4]研究分析了東江干支流年最大洪水遭遇風險特征，結(jié)果表明：Copula函數(shù)能較好地擬合東江干支流年最大洪水(出現(xiàn)時間和量級)聯(lián)合分布。但是，目前的研究多停留在探究水文系列聯(lián)合分布上，對應用高潮位與同期潮差聯(lián)合分布推求潮位過程線的研究尚顯不足。

1 高潮位與同期潮差組合分析方法

1.1 二維Copula聯(lián)合分布模型L

1.1.1 邊緣分布函數(shù)

設年最高潮位為Z，與之相對應的統(tǒng)一潮位過程的較大潮差為R，根據(jù)統(tǒng)計特性，假定年最高潮位和同期潮差分別服從Gumbel分布和正態(tài)分布，其分布函數(shù)分別為：

FZ(z)=exp {-exp [-(z-μ)/σ]}

(1)

s

式中：μ、σ為統(tǒng)計參數(shù)，可用極大似然法對其進行估計。

對于假設的Gumbel分布和正態(tài)分布是否符合實際要求，需要對其進行假設檢驗。對于年最高潮位的Gumbel分布檢驗常使用K-S檢驗法，指標的表達式見式(3)；國際標準化組織組織統(tǒng)計學家對幾十種正態(tài)性檢驗方法進行比較后，認為Wilk-Shapiro的W檢驗效果最好[5]。因此選用W檢驗來判斷同期潮差的分布是否符合正態(tài)分布，指標表達式見式(4)。

K-S檢驗：

Dn=max {|F(x)-Fn(x)|}

(3)

式中：Fn(x)=i/(n+1)是等于或小于x的所有觀察結(jié)果數(shù)目；F(x)表示理論分布的累積概率分布函數(shù)。

W檢驗：

(4)

當W>Wα=0.05時，認為在α=0.05的水平上不拒絕正態(tài)性假設。

1.1.2 二維Copula聯(lián)合分布模型

Copula函數(shù)的優(yōu)點在于分別考慮了變量的邊緣分布和變量間的相關(guān)性結(jié)構(gòu)，且能描述變量間的非線性相關(guān)關(guān)系。目前在水文領域常用的Copula函數(shù)主要有：Gumbel-Hougaard Copula函數(shù)、Clayton Copula函數(shù)、Frank Copula函數(shù)和AMH Copula函數(shù)。這4種Copula函數(shù)各有特點，可以描述各種相關(guān)關(guān)系的水文系列。

Copula 函數(shù)的參數(shù)θ可由Kendall 秩相關(guān)系數(shù)來推算。Kendall 秩相關(guān)系數(shù)用來描述年最高潮位和同期潮差之間的相依關(guān)系，計算公式為：

(5)

式中：τ為 Kendall秩相關(guān)系數(shù)，-1≤τ≤1，τ的絕對值越大表示變量間的相關(guān)關(guān)系越強；xi-xj為實測點據(jù)；n為系列長度。

4種常用的Copula函數(shù)表達式及Kendall秩相關(guān)系數(shù)τ與Copula 函數(shù)的參數(shù)θ間的關(guān)系見表1。

1.1.3 擬合優(yōu)度評價

上述4種Copula函數(shù)對變量間相關(guān)性的描述各有特點，可通過擬合優(yōu)度評價選擇最優(yōu)Copula 函數(shù)。常用的擬合優(yōu)度評價方法有以下3種：

表1 常見的Copula函數(shù)及Kendall秩相關(guān)系數(shù)τ與 θ間的關(guān)系Tab.1 Common Copula functions and the relationship betweenrank correlation coefficients τ and θ

注：u和v均為邊緣分布函數(shù)，u=FZ(z)，v=FR(r)；θ為Copula函數(shù)的參數(shù)；τ為 Kendall 秩相關(guān)系數(shù)。

(1)最小二乘法評價(OLS)。

(6)

(2)最小信息準則(AIC)。

(7)

AIC=nlnMSE+2k

(8)

式中：Femp(xi,yi)為經(jīng)驗頻率；C(ui,vi)為理論頻率；k為模型參數(shù)的個數(shù)。

計算出的OLS或MSE越小，說明擬合優(yōu)度越高。

(3) 理論聯(lián)合概率與經(jīng)驗聯(lián)合頻率比較。繪制由理論聯(lián)合分布和實測點據(jù)計算所得的經(jīng)驗聯(lián)合分布的散點圖進行比較。如果圖上的點都落在45°線附近，表明所選的Copula函數(shù)擬合效果良好。經(jīng)驗聯(lián)合分布的計算公式為：

(9)

式中：xi,yi為實測點據(jù)；nm,k為同時滿足X≤xi,Y≤yi時聯(lián)合觀測值的個數(shù)；n為聯(lián)合觀測值的總個數(shù)。

1.2 設計組合風險分析

在防洪設計中，影響安全的主要因素是高潮位。當設計標準的高潮位發(fā)生時，相應出現(xiàn)的潮差有多種可能，且不同的潮差出現(xiàn)的概率也不相同。因此需研究選定設計高潮位時，不同頻率潮差與之組合的概率風險情況。在滿足設計高潮位z條件下，將潮差超過r發(fā)生的概率作為其組合風險率：

(10)

式中：Z和R分別為高潮位和潮差變量；FZ(z)為高潮位的邊緣分布函數(shù)。

1.3 條件最可能組合

設計潮位過程的確定需要確定高潮位與潮差的組合方式，即需在給定高潮位Z的條件下，確定高潮位對應的最大潮差R。當高潮位一定時，每一個可能的同期潮差對應的概率分布為：

(11)

相應的條件密度函數(shù)為：

fR|Z(r)=cθ(u,v)fR(r)

(12)

式中：cθ(u,v)為Cθ(u,v)的密度函數(shù)；fR(r)為FR(r)的密度函數(shù)。

當條件概率密度fR|Z(r)取到最大值時，對應的高潮位與潮差組合(z,r0)即稱為條件最可能組合，r0稱為條件最可能潮差。此外，可以給定某個顯著性水平α，計算得到條件最可能同期潮差r0的區(qū)間估計[r1,r2]，通過區(qū)間估計可以定量評價條件最可能同期潮差的不確定性。

2 鎮(zhèn)海站設計潮位過程計算

2.1 年最高潮位與同期潮差二維聯(lián)合Copula聯(lián)合分布模型

選取鎮(zhèn)海站1971-2015年實測潮位系列，假定年最高潮位服從Gumbel分布，同期潮差服從和Normal分布，對邊緣分布函數(shù)的參數(shù)進行最大似然估計，結(jié)果見表2。

表2 邊緣分布函數(shù)參數(shù)最大似然估計值Tab.2 Maximum likelihood estimation of marginaldistribution function parameters

對Gumbel分布函數(shù)進行K-S檢驗，計算得到Dn=0.023 78，查表得Dα=0.05n=42=0.209 85，易知Dn

對Normal分布函數(shù)進行W檢驗，計算得到W=0.967 998，查表Wα=0.05=0.942(n=42)，易知W>Wα=0.05，因此在α=0.05水平上潮差分布不拒絕正態(tài)性假設。

因此，認為在 0.05 的顯著水平時，年最高潮位服從 Gumbel 分布，同期潮差服從正態(tài)分布。圖1、圖2分別為年最高潮位和同期潮差的概率分布圖，可以看出，各單變量的理論分布與經(jīng)驗點據(jù)均擬合良好，證明采用的分布函數(shù)合理。

圖1 年最高潮位理論分布與經(jīng)驗分布對比Fig.1 Comparison of theoretical distribution and empirical distribution of annual maximum tidal levels

圖2 同期潮差理論分布與經(jīng)驗分布對比Fig.2 Comparison of tidal theoretical distribution and the empirical distribution over the same period

為了計算Copula函數(shù)的參數(shù)θ，首先根據(jù)式(5)計算得到秩相關(guān)系數(shù)τ=0.277 962，可見2者為弱正相關(guān)關(guān)系。由 Kendall 秩相關(guān)系數(shù)分別求出4種Copula函數(shù)的參數(shù)值，見表3。

表3 4種Archimedean Copula的參數(shù)θ值Tab.3 The parameter θ value of four kinds ofArchimedean Copula functions

為選擇最優(yōu)Copula 函數(shù)，根據(jù)式(6)～(8)對計算出4種Copula函數(shù)的優(yōu)選度評價指標見表4。

表4 4種Copula函數(shù)的優(yōu)選度評價指標Tab.4 Optimization index of four kinds of Copula functions

由表2可看出Frank Copula的OLS與AIC值均最小，因此選用Frank Copula擬合程度最好。由此可得年最高潮位和同期潮差的聯(lián)合分布為：

F(Z,R)=-

(13)

為進一步檢查所選的Frank Copula函數(shù)的合理性，對所選函數(shù)進行理論聯(lián)合分布與經(jīng)驗聯(lián)合分布相關(guān)性分析，并對點距進行線性擬合，見圖3。

圖3 理論聯(lián)合分布與經(jīng)驗聯(lián)合分布相關(guān)性分析Fig.3 Correlation analysis of theoretical joint distribution and empirical joint distribution

由圖3可看出理論聯(lián)合分布與經(jīng)驗聯(lián)合分布的相關(guān)系數(shù)高達0.988 6，且均落在45°線兩側(cè)，說明擬合質(zhì)量較高，所選用的Frank Copula函數(shù)是合適的。

2.2 設計組合風險分析

由式(12)分別計算出設計頻率50%、20%、10%、5%、2%、1%的高潮位和同期潮差的聯(lián)合分布概率，由式(10)計算出高潮位和同期潮差的組合風險率見表5。

表5 設計高潮位與同期潮差組合風險率Tab.5 The portfolio risk rate of design high tide andtidal range of corresponding period

由表5可知，鎮(zhèn)海潮位站低頻率潮位與低頻率潮差遭遇的概率極小。其原因在于高潮位的出現(xiàn)是由風暴潮增水引起，增水過程引起高潮位升高的同時也會一定程度上引起低潮位提升，因而潮差不一定會增大。

2.3 設計潮位過程線推求

因為高潮位與潮差的發(fā)生不具有同頻率的特征，且對防洪標準影響更大的是高潮位，所以潮位過程線的設計應以高潮位設計為主，同時考慮在給定高潮位的條件下，出現(xiàn)幾率最大的潮差。由式(12)計算得到年最高潮位相應的條件最可能潮差見表6。

表6 設計高潮位及條件最可能潮差Tab.6 The design high tide level and conditionalmost likely tidal range

由表6可知，鎮(zhèn)海站高潮位與條件最可能潮差不具有同頻率性，各重現(xiàn)期的高潮位對應的條件最可能潮差均為2～5 a一遇。這種低頻高潮位、高頻低潮差的組合更不利于河道行洪，容易造成城市洪澇災害。比如，近5 a的2012年“?？薄?013年“菲特”、2015年“燦鴻”、“杜鵑”等臺風期間的潮位過程均是低頻高潮位與高頻潮差的組合，對甬江洪水外排造成外海潮位持續(xù)頂托作用，導致寧波城市澇水難以外排，造成嚴重的洪水災害。

對防洪排澇影響較大的是高潮位，推求設計潮位線首先選定設計高潮位，然后根據(jù)表6，按照設計年最高潮位和條件最可能潮差雙重控制，等比例縮放典型潮型。以50 a一遇設計潮位過程為例，選取2000年“桑美”臺風期間的實測潮位過程作為典型潮型，縮放得到的設計潮位過程見圖4。這種方法推求的潮位過程線充分考慮了高潮位和潮差組合的因素，為水利計算提供更合理的下邊界，使水利工程設計也更安全。

圖4 設計潮位過程Fig.4 Design tidal process

3 結(jié) 論

(1)根據(jù)優(yōu)選度評價，選用Frank Copula函數(shù)構(gòu)建鎮(zhèn)海站年最高潮位與同期潮差的二維聯(lián)合分布合理可行。

(2)由設計組合分析可知，鎮(zhèn)海潮位站高潮位與潮差遭遇不具有同頻率性，低頻率高潮位與低頻率潮差遭遇的可能性極小。

(3)通過條件最可能組合分析，鎮(zhèn)海站各重現(xiàn)期的高潮位對應的條件最可能潮差均為2～5 a一遇。

(4)提出了一種設計高潮位和條件最可能潮差雙重控制，等比例縮放典型的方法來推求設計潮位過程線，這樣充分考慮的潮差的影響，使工程設計更合理。

□

[1] 楊志勇,袁 喆,方宏陽,等. 基于Copula函數(shù)的灤河流域旱澇組合事件概率特征分析[J]. 水利學報,2013,(5):556-561,569.

[2] 張冬冬,魯 帆,嚴登華,等. 基于Archimedean Copula函數(shù)的洪水多要素聯(lián)合概率分布研究[J]. 中國農(nóng)村水利水電,2015,(1):68-74,79.

[3] 李天元,郭生練,劉章君,等. 基于峰量聯(lián)合分布推求設計洪水[J]. 水利學報,2014,(3):269-276.

[4] 黃鋒華,黃本勝,郭 磊,等. 東江干流與支流河涌洪水遭遇風險研究[J]. 中國農(nóng)村水利水電,2016,(3):144-148.

[5] 茆詩松,濮曉龍,程依明.概率論與數(shù)理統(tǒng)計[M].高等教育出版社,2000:331.