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

(1.寧波市水利水電規(guī)劃設(shè)計(jì)研究院，浙江 寧波 315192；2.河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210098)

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

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

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

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

1.1.1 邊緣分布函數(shù)

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

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

(1)

s

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

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

K-S檢驗(yàn)：

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

(3)

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

W檢驗(yàn)：

(4)

當(dāng)W>Wα=0.05時(shí)，認(rèn)為在α=0.05的水平上不拒絕正態(tài)性假設(shè)。

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

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

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

(5)

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

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

1.1.3 擬合優(yōu)度評(píng)價(jià)

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

表1 常見(jiàn)的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)最小二乘法評(píng)價(jià)(OLS)。

(6)

(2)最小信息準(zhǔn)則(AIC)。

(7)

AIC=nlnMSE+2k

(8)

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

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

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

(9)

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

1.2 設(shè)計(jì)組合風(fēng)險(xiǎn)分析

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

(10)

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

1.3 條件最可能組合

設(shè)計(jì)潮位過(guò)程的確定需要確定高潮位與潮差的組合方式，即需在給定高潮位Z的條件下，確定高潮位對(duì)應(yīng)的最大潮差R。當(dāng)高潮位一定時(shí)，每一個(gè)可能的同期潮差對(duì)應(yīng)的概率分布為：

(11)

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

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

(12)

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

當(dāng)條件概率密度f(wàn)R|Z(r)取到最大值時(shí)，對(duì)應(yīng)的高潮位與潮差組合(z,r0)即稱為條件最可能組合，r0稱為條件最可能潮差。此外，可以給定某個(gè)顯著性水平α，計(jì)算得到條件最可能同期潮差r0的區(qū)間估計(jì)[r1,r2]，通過(guò)區(qū)間估計(jì)可以定量評(píng)價(jià)條件最可能同期潮差的不確定性。

2 鎮(zhèn)海站設(shè)計(jì)潮位過(guò)程計(jì)算

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

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

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

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

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

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

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

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

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

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

為選擇最優(yōu)Copula 函數(shù)，根據(jù)式(6)～(8)對(duì)計(jì)算出4種Copula函數(shù)的優(yōu)選度評(píng)價(jià)指標(biāo)見(jiàn)表4。

表4 4種Copula函數(shù)的優(yōu)選度評(píng)價(jià)指標(biāo)Tab.4 Optimization index of four kinds of Copula functions

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

F(Z,R)=-

(13)

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

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

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

2.2 設(shè)計(jì)組合風(fēng)險(xiǎn)分析

由式(12)分別計(jì)算出設(shè)計(jì)頻率50%、20%、10%、5%、2%、1%的高潮位和同期潮差的聯(lián)合分布概率，由式(10)計(jì)算出高潮位和同期潮差的組合風(fēng)險(xiǎn)率見(jiàn)表5。

表5 設(shè)計(jì)高潮位與同期潮差組合風(fēng)險(xiǎn)率Tab.5 The portfolio risk rate of design high tide andtidal range of corresponding period

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

2.3 設(shè)計(jì)潮位過(guò)程線推求

因?yàn)楦叱蔽慌c潮差的發(fā)生不具有同頻率的特征，且對(duì)防洪標(biāo)準(zhǔn)影響更大的是高潮位，所以潮位過(guò)程線的設(shè)計(jì)應(yīng)以高潮位設(shè)計(jì)為主，同時(shí)考慮在給定高潮位的條件下，出現(xiàn)幾率最大的潮差。由式(12)計(jì)算得到年最高潮位相應(yīng)的條件最可能潮差見(jiàn)表6。

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

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

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

圖4 設(shè)計(jì)潮位過(guò)程Fig.4 Design tidal process

3 結(jié) 論

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

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

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

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

□

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

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

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

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

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