甘富萬，黃永俊，熊　健，2，孫晉東，胡　泊

(1.廣西大學(xué)土木建筑工程學(xué)院， 廣西南寧530004；2.廣西水利電力勘測設(shè)計研究院， 廣西南寧530023)

?

基于Copula函數(shù)的桂平航運樞紐水閘設(shè)計洪水位

甘富萬1，黃永俊1，熊健1，2，孫晉東1，胡泊1

(1.廣西大學(xué)土木建筑工程學(xué)院， 廣西南寧530004；2.廣西水利電力勘測設(shè)計研究院， 廣西南寧530023)

摘要：針對桂平航運樞紐水閘原設(shè)計洪水位由經(jīng)驗方法確定、缺乏理論依據(jù)的情況，利用Copula函數(shù)理論建立二變量設(shè)計洪水組合及洪水位的確定方法，并應(yīng)用于桂平航運樞紐水閘設(shè)計洪水位的確定。結(jié)果表明：建立的4種常用Copula函數(shù)中，Clayton Copula函數(shù)擬合效果最好；設(shè)計洪水組合應(yīng)該包括Copula同現(xiàn)重現(xiàn)期等于設(shè)計重現(xiàn)期的所有洪水組合；按照建立的方法計算得到桂平航運樞紐水閘設(shè)計洪水位為43.26 m，通過經(jīng)驗方法確定的原設(shè)計洪水位為43.48 m，可達到防洪要求，且略偏安全。

關(guān)鍵詞：Copula函數(shù)；桂平航運樞紐水閘；重現(xiàn)期；設(shè)計洪水組合；設(shè)計洪水位

桂平航運樞紐位于廣西桂平市郁江河段，是一個集航運、發(fā)電、灌溉、交通于一體的綜合利用性航運樞紐工程，處于郁江與西江干流(包括黔江、潯江)匯合口附近的郁江段，主要由航運船閘、溢洪壩及水閘組成[1]，見圖1。非洪水時期，溢流壩水閘關(guān)閉，樞紐壩上水位抬高以保證航運水位以及發(fā)電；洪水時期，溢洪壩水閘打開，洪水通過溢洪壩水閘下泄。桂平航運樞紐同時受到郁江和西江干流洪水的作用，匯合口的郁江上游附近有貴港水文站、干流下游附近有大湟江口站水文站[2]。桂平航運樞紐水閘設(shè)計重現(xiàn)期為100 a一遇，設(shè)計時采用大湟江口站重現(xiàn)期為100 a的洪峰流量與0.2倍該洪峰值作為溢洪壩壩址設(shè)計洪水組合，并經(jīng)過調(diào)洪計算得到溢洪壩水閘設(shè)計洪水位為43.48 m。該設(shè)計洪水位的確定方法為經(jīng)驗方法，缺乏理論依據(jù)，所得設(shè)計洪水位是否達到設(shè)計重現(xiàn)期的要求，有待商榷。

圖1　桂平航運樞紐位置及組成

桂平航運樞紐水閘設(shè)計洪水問題是典型二變量洪水組合問題。Copula函數(shù)可描述變量間的相關(guān)性，是一類將聯(lián)合分布函數(shù)與邊緣分布函數(shù)連接在一起的函數(shù)，在多變量水文研究中已經(jīng)被證明有著廣泛的推廣應(yīng)用價值，相關(guān)理論參見文獻[3]～文獻[11]。本文嘗試?yán)肅opula函數(shù)對桂平航運樞紐的設(shè)計洪水組合進行計算，并以其為基礎(chǔ)進一步確定桂平航運樞紐溢洪壩水閘的設(shè)計洪水位。

1Copula函數(shù)的建立

1.1函數(shù)的選擇

首先選擇已被證明在多變量水文問題中可以有效應(yīng)用的四種Copula函數(shù)，然后根據(jù)各種Copula函數(shù)的要求，利用大湟江口站和貴港站歷年實測洪峰系列資料建立Copula函數(shù)并進行對比評價，以選擇最為合適的Copula函數(shù)形式。初步選擇的4種Copula函數(shù)形式見表1。

表1　常用copula函數(shù)形式

表1中各式的C為Copula連接函數(shù)，u1和u2為變量X1和X2的邊緣分布函數(shù)，u1=F1(x1)，u2=F2(x2)，θ為Copula函數(shù)的參數(shù)。

1.2函數(shù)的建立

①邊緣分布

邊緣分布統(tǒng)一采用為P-III型分布。為充分考慮郁江、西江干流洪水的遭遇特點，本次分兩種情況分別建立Copula函數(shù)：以西江干流洪峰系列與相應(yīng)時間的郁江洪峰系列進行組合(方案1)；以郁江洪峰系列與相應(yīng)時間的西江干流洪峰系列進行組合(方案2)。

方案1具體做法：西江干流選擇大湟江口站1951年～2010年實測洪水系列，郁江以大湟江口站發(fā)生最大洪峰流量的時間為控制，統(tǒng)計貴港站相應(yīng)場次洪水的相應(yīng)洪峰流量并進行頻率計算和P-Ⅲ型曲線適線。

方案2具體做法：郁江選擇貴港站1936年～2010年實測洪水系列資料，西江干流以貴港站發(fā)生最大洪峰流量的時間為控制，統(tǒng)計大湟江口站同場洪水的相應(yīng)洪峰流量，并以此統(tǒng)計P-Ⅲ型曲線參數(shù)。

P-III型分布曲線也稱伽馬分布，它的概率密度函數(shù)[12]是：

(1)

式中，Γ(α)是α的伽瑪函數(shù)；α、β、a0分別是P-Ⅲ型分布函數(shù)的形狀、尺度和位置參數(shù)，α>0，β>0。

(2)

經(jīng)過頻率計算和P-Ⅲ型曲線適線分析，得到邊緣分布參數(shù)統(tǒng)計結(jié)果見表2和表3。

表2　邊緣分布函數(shù)的統(tǒng)計參數(shù)值

表3　邊緣分布函數(shù)的α、β、a0參數(shù)值

②Copula函數(shù)參數(shù)估計

采用相關(guān)性指標(biāo)法[12]計算Copula函數(shù)的計算參數(shù)θ，即首先計算各變量序列的Kendall秩τ，然后再利用τ計算θ。各Copula函數(shù)參數(shù)θ與Kendall秩相關(guān)系數(shù)τ之間的關(guān)系見表4。

Kendall秩相關(guān)系數(shù)τ的計算方法可參考文獻[13]，分別利用實測數(shù)據(jù)可求得方案1和方案2的τ為0.147 6和0.219 6，并進一步計算得到各Copula函數(shù)的參數(shù)θ值見表4。

表4　Copula函數(shù)參數(shù)θ與 Kendall 秩相關(guān)系數(shù)τ轉(zhuǎn)換關(guān)系

1.3評價分析

本文利用Copula函數(shù)常見的AIC信息準(zhǔn)則法以及離差平方和最小準(zhǔn)則(OLS)法進行評價[13]。

AIC信息準(zhǔn)則形式如下：

(3)

AIC=nln(MSE)+2k，

(4)

式中：Femp(xi1,xi2,…,xim)表示經(jīng)驗頻率值，C(ui1,ui2,…,uim)表示理論頻率值，m為函數(shù)的維數(shù)，k為模型參數(shù)的個數(shù)。

離差平方和最小準(zhǔn)則形式如下：

(5)

式中各符號意義同上。

AIC值、OLS值越小，Copula函數(shù)擬合越好。運用AIC信息準(zhǔn)則以及離差平方和最小準(zhǔn)則進行計算，結(jié)果見表5。由計算成果可見，各copula函數(shù)的AIC值、OLS值非常接近，可見擬合效果都比較好。其中AIC值、OLS值最小的函數(shù)為Clayton Copula函數(shù)，下文主要利用Clayton Copula函數(shù)成果進行分析。

表5　AIC值和OLS值計算成果

2基于Copula函數(shù)的設(shè)計洪水組合

2.1設(shè)計洪水重現(xiàn)期

桂平航運樞紐水閘設(shè)計洪水重現(xiàn)期為100 a一遇。Copula函數(shù)的重現(xiàn)期分為聯(lián)合重現(xiàn)期、同現(xiàn)重現(xiàn)期和條件重現(xiàn)期三種，不同類型重現(xiàn)期對應(yīng)的洪水組合不同，必須明確哪種重現(xiàn)期可以作為水閘設(shè)計洪水的重現(xiàn)期。

假定洪水事件中的洪峰系列分別為X1和X2。u1=F1(x1)=P(X1>x1)、u2=F2(x2)=P(X2>x2)為對應(yīng)洪水邊緣分布函數(shù)，則對于系列X1的單變量洪水重現(xiàn)期可表示為：

T(x1)=1/P(X1>x1)。

(6)

當(dāng)以此重現(xiàn)期作為設(shè)計重現(xiàn)期T設(shè)，設(shè)計洪峰為(x1)，可以保證洪水事件(X1x1)，行洪安全無法得到保證。

比較單變量設(shè)計洪水重現(xiàn)期的定義，在二變量洪水組合問題中，如果設(shè)計的洪水組合為(x1,x2)，則發(fā)生洪水事件(X1x1,X2>x2)時，行洪不能保證按設(shè)計工況安全通過。此時設(shè)計洪水重現(xiàn)期T設(shè)為：

T設(shè)=T(x1,x2)=1/P(X1>x1,X2>x2)。

(7)

Copula函數(shù)同現(xiàn)重現(xiàn)期是指兩個變量同時超過特定值時事件發(fā)生的重現(xiàn)期，計算公式如下：

T0(x1,x2)=1/P(X1>x1,X2>x2)=1/(1-u1-u2+C(u1,u2))。

(8)

對比可知，應(yīng)該以Copula函數(shù)的同現(xiàn)重現(xiàn)期作為本次二變量設(shè)計洪水重現(xiàn)期T設(shè)。

2.2設(shè)計洪水組合

對于單變量洪水事件，洪峰—重現(xiàn)期關(guān)系為一條曲線，如圖2所示。設(shè)計洪水重現(xiàn)期T設(shè)(如100 a一遇)對應(yīng)的洪峰為曲線上該重現(xiàn)期對應(yīng)的一個值，以該洪峰值作為設(shè)計洪峰，可保證重現(xiàn)期小于100 a一遇的洪峰出現(xiàn)時的行洪安全。

在本次研究的情況中，水閘同時受到郁江、西江干流洪水作用，是二變量洪水組合問題。T設(shè)(如100 a一遇)對應(yīng)的洪水不是一個洪峰，而是郁江、西江干流的洪峰組合，即設(shè)計洪水組合。郁江、西江干流洪峰與重現(xiàn)期的關(guān)系為一個曲面，如圖3所示。T設(shè)(如100 a一遇)對應(yīng)的洪峰組合則構(gòu)成了一條曲線，而重現(xiàn)期小于100 a一遇對應(yīng)的洪水組合構(gòu)成一個曲面。要使得設(shè)計洪水重現(xiàn)期達到100 a一遇，則必須保證所有重現(xiàn)期小于100 a一遇的洪水組合的行洪安全，此時設(shè)計洪水組合應(yīng)該包含設(shè)計重現(xiàn)期(如100 a一遇)對應(yīng)的曲線上的所有洪水組合。

圖2單變量洪水事件設(shè)計洪峰

Fig.2Design flood of single variable flood

圖3二變量洪水事件設(shè)計洪水組合

Fig.3Design flood combination of two variable flood

3基于Copula函數(shù)的水閘設(shè)計洪水位

3.1計算方法

利用建立的Clayton Copula函數(shù)計算得到桂平航運樞紐同現(xiàn)重現(xiàn)期為100 a一遇的洪水組合，并進行調(diào)洪計算可得到水閘的設(shè)計洪水位。由于同現(xiàn)重現(xiàn)期100 a一遇對應(yīng)的洪水組合有無數(shù)個，本次采用抽取大量代表洪水組合進行調(diào)洪計算，并取調(diào)洪水位最大值作為設(shè)計值的方法。

設(shè)計洪水組合下的調(diào)洪計算方法直接參考樞紐原設(shè)計的計算方法[1]。

水閘泄洪能力的計算公式如下：

(9)

圖4　壩址水位流量關(guān)系Fig.4　Water level and flux relationship of dam site

壩址自然水位同時受到郁江和西江干流洪水的作用，根據(jù)不同的設(shè)計洪水組合，利用“貴港站流量—大湟江口站流量—壩址自然水位”關(guān)系插值計算得到不同的設(shè)計洪水組合下的閘下水位，再根據(jù)公式(9)計算得到閘上水位，即為不同設(shè)計洪水組合對應(yīng)的的水閘設(shè)計洪水位。

利用實測資料整理溢洪壩壩址的自然水位流量關(guān)系[1]，見圖4。

3.2計算成果

針對方案1和方案2兩種方案建立的Clayton Copula函數(shù)，計算大量重現(xiàn)期為100 a一遇的洪水組合并進行調(diào)洪計算，得到各設(shè)計洪水組合的閘上洪水位，典型設(shè)計洪水組合計算結(jié)果見表6。

表6　典型設(shè)計洪水組合調(diào)洪計算成果(T設(shè)=100 a)

由計算成果可見，桂平航運樞紐水閘閘上水位主要受到西江干流洪水的控制，在各設(shè)計洪水組合中，西江干流洪水越大，閘上洪水位越高。洪水位最大的洪水組合為西江干流洪水接近100 a一遇(對應(yīng)西江干流洪峰48 200 m3/s，郁江洪峰966 m3/s)的洪水組合，對應(yīng)洪水位為43.26 m，此值即為桂平航運樞紐水閘的設(shè)計洪水位。

桂平航運樞紐水閘原設(shè)計采用的設(shè)計洪水組合為西江干流48 200 m3/s、郁江9 640 m3/s，對應(yīng)原設(shè)計洪水位43.48 m，高于本次利用Copula函數(shù)計算得到的設(shè)計洪水位0.22 m?？梢?，原設(shè)計洪水位雖然按照經(jīng)驗方法確定，但達到了100 a一遇設(shè)計重現(xiàn)期的要求，而且略偏安全。

4結(jié)語

本文建立了郁江和西江干流洪水的Copula函數(shù)，并探討了基于Copula函數(shù)的二變量設(shè)計洪水組合及設(shè)計洪水位的計算方法，最后對桂平航運樞紐水閘的設(shè)計洪水位進行了計算。結(jié)論如下：

①建立的4種常用Copula函數(shù)中，Clayton Copula函數(shù)擬合效果最好。

②設(shè)計洪水組合應(yīng)該包括Copula同現(xiàn)重現(xiàn)期等于設(shè)計重現(xiàn)期的所有洪水組合。

③按照本文建立的方法計算得到桂平航運樞紐水閘設(shè)計洪水位為43.26 m。通過經(jīng)驗方法確定的原設(shè)計洪水位43.48 m可達到防洪要求，且略偏安全

參考文獻：

[1]徐永樂.西江航運建設(shè)第一期工程廣西段初步設(shè)計報告(第一冊 桂平航運梯級)[R]. 南寧：廣西壯族自治區(qū)交通勘測設(shè)計院，1983.

[2]甘富萬，胡秀英，劉欣，等.廣西境內(nèi)西江流域洪水特性分析[J]. 廣西大學(xué)學(xué)報(自然科學(xué)版)，2015，40(1)：244-250.

[3]GRIMALDI S, SERINALDI F.Design hyetographs analysis with 3 -Copula function [J]. Hydrological Sciences, 2006, 51(2): 223-238.

[4]SHIAU J T, SONG F, NADARAJAH S.Assessment of hydrological droughts for the Yellow River, China, using Copulas [J]. Hydrological Processes, 2007, 21: 2157-2163.

[5]DEMICHELE C, SALVADORI G, CANOSSI M, et al.Bivariate statistical approach to check adequacy of dam spillway[J]. Journal of Hydrologic Engineering, 2005, 10(1): 50-57.

[6]郭生練，閆寶偉，肖義，等.Copula函數(shù)在多變量水文分析計算中的應(yīng)用及研究進展[J]. 水文，2008，28(3)：1-7.

[7]張娜，郭生練，閆寶偉，等.Copula函數(shù)在分期設(shè)計洪水中的應(yīng)用研究[J]. 水文，2008，28(5)：28-32.

[8]張娜，郭生練，劉攀，等.基于Copula函數(shù)法推求分期設(shè)計洪水和汛限水位[J]. 武漢大學(xué)學(xué)報(工學(xué)版)，2008，41(6)：33-40.

[9]閆寶偉，郭生練，郭靖，等.基于Copula 函數(shù)的設(shè)計洪水地區(qū)組成研究[J]. 水力發(fā)電學(xué)報，2010，29(6)：60-65.

[10]陸桂華，張亞洲，肖恒，等.氣候變化背景下蚌埠市暴雨與淮河上游洪水遭遇概率分析[J]. 氣候變化研究進展，2015，11 (1)：31-37.

[11]劉章君，郭生練，李天元，等.梯級水庫設(shè)計洪水最可能地區(qū)組成法計算通式[J]. 水科學(xué)進展，2014，25(4)：575-584.

[12]侯蕓蕓，宋松柏，趙麗娜，等.基于Copula 函數(shù)的3 變量洪水頻率研究[J]. 西北農(nóng)林科技大學(xué)學(xué)報(自然科學(xué)版)，2010，38(2)：219-228.

[13]于延勝，陳興偉.基于Mann_Kendall法的徑流豐枯變化過程劃分[J]. 水資源與水工程學(xué)，2013，24(1)：60-63.

(責(zé)任編輯唐漢民梁碧芬)

Design flod level of sluice of Guiping Shipping Hub based on Copula function

GAN Fu-wan1, HUANG Yong-jun1, XIONG Jian2, SUN Jin-dong， HU Bo1

(1. College of Civil Engineering and Architecture, Guangxi University, Nanning 530004, China;

2. Guangxi Water and Power Design Institute, Nanning 530023, China)

Abstract:In view of that the design flood level of the Guiping Shipping Hub was calculated by an empirical method which lacked of theoretical basis, a new method based on Copula function to calculate the design flood combination of two variables and the flood level was established and applied to the sluice of the Guiping Shipping Hub. The result shows that Clayton Copula has the best fitting effect among the four common Copula functions; that the design flood combination should include all combinations, in which Copula function occurrence recurrence period is equal to the design recurrence period of the Guiping Shipping Hub sluice; that the design flood level of the sluice of Guiping Shipping Hub is 43.26 m according to the new method; that the original design flood level (43.48 m) is slightly higher and meet the design standard.

Key words:Copula function; Sluice of Guiping Shipping Hub; recurrence period；design flood combination; design flood level

中圖分類號：TV122+.3

文獻標(biāo)識碼：A

文章編號：1001-7445(2016)01-0284-07

doi:10.13624/j.cnki.issn.1001-7445.2016.0284

通訊作者：甘富萬(1981—)，男，廣西北流人，廣西大學(xué)副教授，博士；E-mail:gxugfw@126.com。

基金項目：國家自然科學(xué)基金資助項目(51178124)；廣西自然科學(xué)基金資助項目(2014GXNSFBA118263)；廣西大學(xué)科研基金資助項目(土科基-2012-02)

收稿日期：2015-11-20；

修訂日期：2015-12-21

引文格式：甘富萬，黃永俊，熊健,等.基于Copula函數(shù)的桂平航運樞紐水閘設(shè)計洪水位[J].廣西大學(xué)學(xué)報(自然科學(xué)版)，2016，41(1)：284-290.