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

薛曉鵬，張衛國，朱從飛，趙思遠，郝振純

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

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

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

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

1.1 二維Copula聯合分布模型L

1.1.1 邊緣分布函數

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

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

(1)

s

式中：μ、σ為統計參數，可用極大似然法對其進行估計。

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

K-S檢驗：

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

(3)

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

W檢驗：

(4)

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

1.1.2 二維Copula聯合分布模型

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

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

(5)

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

4種常用的Copula函數表達式及Kendall秩相關系數τ與Copula 函數的參數θ間的關系見表1。

1.1.3 擬合優度評價

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

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

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

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

(6)

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

(7)

AIC=nlnMSE+2k

(8)

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

計算出的OLS或MSE越小，說明擬合優度越高。

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

(9)

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

1.2 設計組合風險分析

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

(10)

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

1.3 條件最可能組合

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

(11)

相應的條件密度函數為：

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

(12)

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

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

2 鎮海站設計潮位過程計算

2.1 年最高潮位與同期潮差二維聯合Copula聯合分布模型

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

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

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

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

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

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

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

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

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

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

表4 4種Copula函數的優選度評價指標Tab.4 Optimization index of four kinds of Copula functions

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

F(Z,R)=-

(13)

為進一步檢查所選的Frank Copula函數的合理性，對所選函數進行理論聯合分布與經驗聯合分布相關性分析，并對點距進行線性擬合，見圖3。

圖3 理論聯合分布與經驗聯合分布相關性分析Fig.3 Correlation analysis of theoretical joint distribution and empirical joint distribution

由圖3可看出理論聯合分布與經驗聯合分布的相關系數高達0.988 6，且均落在45°線兩側，說明擬合質量較高，所選用的Frank Copula函數是合適的。

2.2 設計組合風險分析

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

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

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

2.3 設計潮位過程線推求

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

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

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

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

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

3 結 論

(1)根據優選度評價，選用Frank Copula函數構建鎮海站年最高潮位與同期潮差的二維聯合分布合理可行。

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

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

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

□

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

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

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

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

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