基于G-H Copula的水庫防洪安全計算方法研究

范嘉煒，黃錦林，袁明道，馬妍博

(1.廣東省水利水電科學(xué)研究院，廣州 510635；2. 廣東省水利水電科學(xué)研究院 河口水利技術(shù)國家地方聯(lián)合工程實驗室，廣州 510635)

0 引 言

水庫確立防洪調(diào)度標(biāo)準(zhǔn)與設(shè)計洪水緊密相關(guān)，且設(shè)計洪水過程線的大小會直接影響水庫的防洪安全設(shè)計。而推求設(shè)計洪水過程線，目前國內(nèi)通常采用的是典型洪水放大方法：同倍比放大法和同頻率放大法。推求設(shè)計洪水通常采用典型洪水放大方法，而水庫確立防洪調(diào)度標(biāo)準(zhǔn)是與設(shè)計洪水緊密相關(guān)的，即計算所得的設(shè)計洪水過程線的大小會直接影響水庫的防洪安全設(shè)計。目前國內(nèi)在推求設(shè)計洪水過程線時通常采用同倍比放大法和同頻率放大法。由于洪水事件包含洪峰流量、洪量、洪水歷時等多個特征量，因此同倍比法在進(jìn)行洪水放大時需要先選定某個洪水特征量再進(jìn)行計算，直接忽略了其他特征量對洪水過程的影響；同頻率放大法與同倍比法相比，雖然能夠同時考慮兩個及多個特征變量，但在推求特征量的設(shè)計值時仍采用單變量的分布函數(shù)，分別計算在設(shè)定重現(xiàn)期條件下的變量設(shè)計值[1,2]。然而在一場洪水中，特征量均為同一重現(xiàn)期的情況僅為一種假設(shè)條件下的理想狀況，忽略了特征量間的關(guān)聯(lián)性，將導(dǎo)致推求設(shè)計洪水過程線出現(xiàn)偏差，進(jìn)而影響水庫的防洪安全設(shè)計[3]。本文通過G-H Copula函數(shù)，以洪峰和72 h洪量為研究對象建立兩變量聯(lián)合模型，推求在考慮兩變量相關(guān)性情況下的設(shè)計洪水過程線[4,5]，進(jìn)而對比分析峰量聯(lián)合法與傳統(tǒng)方法對水庫調(diào)洪的計算結(jié)果。分析結(jié)論可為水庫防洪安全設(shè)計提供有效的科學(xué)依據(jù)。

1 G-H Copula函數(shù)與聯(lián)合重現(xiàn)期

Gumbel-Hougaard (G-H) Copula函數(shù)屬Archimedean-Copula函數(shù)中的一種[6]，它能夠?qū)⒍鄠€變量的邊緣分布聯(lián)結(jié)起來，反映變量之間的相關(guān)關(guān)系。目前在水文領(lǐng)域被廣泛應(yīng)用。Copula函數(shù)由于表達(dá)形式簡潔易懂，且由于在其構(gòu)造基于多變量的聯(lián)合分布函數(shù)時可以變換多種形式，使得在實際的學(xué)術(shù)研究中具有較高的科研價值，在水利工程中也被廣泛應(yīng)用于分析實際問題[7,8]。

Gumbel-Hougaard(G-H)Copula的表達(dá)形式為：

C(u1,u2)=exp{-[(-lnu1)θ+(-lnu2)θ]1/θ},θ∈[1,∞)

(1)

式中：θ為函數(shù)參數(shù)；u1和u2為邊緣分布函數(shù)。

參數(shù)θ的估計通常采用相關(guān)性指標(biāo)法，通過建立kendall秩相關(guān)系數(shù)τ與θ的關(guān)系推求θ值。τ與θ的關(guān)系可由式(2)表示為

(2)

而τ可由式(3)推求：

(3)

式中：τ為kendall秩相關(guān)系數(shù)；sign為符號函數(shù)；(xi,yi)為測點據(jù)；n為實測序列長度。

聯(lián)合重現(xiàn)期目前被水文界普遍認(rèn)可[9,10]。變量的重現(xiàn)期與它們的聯(lián)合分布函數(shù)C相關(guān)聯(lián)，C值越大，說明發(fā)生該水文事件的可能性越小，相應(yīng)的重現(xiàn)期就越大。

若以洪峰流量與72 h洪量為研究對象，即若取組合序列Q和W72，Q和W72為研究對象，將邊緣分布函數(shù)表示為：

(4)

重現(xiàn)期是描述設(shè)計洪水的重要概念之一，眾所周知，在傳統(tǒng)的洪水頻率分析當(dāng)中，洪水特征量重現(xiàn)期定義為其頻率的倒數(shù)，公式表示為：

(5)

式中：TQ、TW分別表示洪峰與72 h洪量的重現(xiàn)期。

傳統(tǒng)的頻率分析忽略了洪峰和洪量之間緊密的正相關(guān)性，為了能夠更加準(zhǔn)確合理地描述洪水事件，考慮洪水峰量的內(nèi)在聯(lián)系，建立聯(lián)合分布是十分必要的。

對于Q和W72，它們的聯(lián)合重現(xiàn)期可以用式(6)表示，即Q或W72中至少有一個超過設(shè)定閾值。聯(lián)合重現(xiàn)期可以用式(6)表示兩個變量之一超過設(shè)定閾值時的工況。

(6)

式中：x和y為設(shè)定的閾值，u、v、C(u,v)意義同前。

如果設(shè)洪峰xT與洪量的yT分別為T年一遇在單變量方法下求得的設(shè)計值，易知聯(lián)合重現(xiàn)期與單變量重現(xiàn)期存在以下關(guān)系：

T=min[(TQ(xT),TW(yT)]≥TY(xT,yT)

(7)

(8)

2 實例研究

2.1 邊緣分布函數(shù)確定

鶴地水庫位于九洲江中游，九洲江為南海北部灣水系，發(fā)源于廣西陸川縣，流經(jīng)陸川縣、博白縣、廣東廉江市，注入北部灣。全長162 km，集雨面積2 137 km2，為廉江市最長和支流最多的河流。在1959-2013年中每年選取鶴地水庫1～3場大洪水作為研究對象，每場洪水歷時一般在72 h左右，可以包括一次洪水起漲落平的全過程。由于在推求設(shè)計洪水過程線以及水庫的防洪安全設(shè)計時，一般選取歷時較長的洪水過程，因此，本文重點分析洪峰與3 d洪量的聯(lián)合分布。

本次以鶴地水庫1958-2013年洪水監(jiān)測資料為依據(jù)，選取年最大洪峰流量Q和W72為研究序列，建立Q-W72峰量組合樣本。根據(jù)《水利水電工程設(shè)計洪水規(guī)范》(SL44-2006)，考慮歷史洪水，并對實測系列進(jìn)行特大洪水處理后樣本數(shù)據(jù)采用P-Ⅲ型曲線頻率適線，分析結(jié)果如表1、圖1所示。

表1鶴地站設(shè)計洪水參數(shù)統(tǒng)計

特征量參數(shù)估計ExCvCs設(shè)計值50 a100 a200 a500 aQ/(m3·s-1)1 5540.621.584 225.224 807.305 383.056 136.69W72/億m31.350.611.623.644.154.655.30

圖1 鶴地水庫年最大洪峰與72 h洪量P-Ⅲ型曲線

2.2 G-H Copula峰量聯(lián)合分布模型

通過式(3)計算得到Q-W72組合的秩相關(guān)系數(shù)τ為0.630，因此兩組變量間存在一定的正相關(guān)性，故建立峰量組合的聯(lián)合分布函數(shù)，其中函數(shù)參數(shù)θ由τ-θ關(guān)系推求：

(9)

為了確認(rèn)G-H Copula函數(shù)是否能準(zhǔn)確擬合變量間的相關(guān)關(guān)系，本文采用較為直觀的圖形分析法進(jìn)行擬合優(yōu)度評價[11]。該方法將經(jīng)驗概率和理論概率分別標(biāo)繪于同一圖中，若當(dāng)散點均分布在45°線周圍時，則變量組合具有較好的相關(guān)性。變量組合的經(jīng)驗概率值可由Femp(xi1,xi2)表示為：

Femp(xi,yj)=P(Q≤xi,W72≤yj)=nij/(N+1)

(10)

式中：N為實測序列的總對數(shù)；nij為序列中小于等于(xi,yj)的對數(shù)；(xi,yj)為Q和W72的聯(lián)合觀測值。

將以上兩組變量為研究對象計算得出的理論頻率值和式(10)得到的經(jīng)驗頻率值點繪于圖2中，觀察可得，Gumbel-Hougaard Copula構(gòu)造的兩組降雨序列聯(lián)合分布的經(jīng)驗點值與理論計算值均勻分布在45°線附近，因此G-H Copula擬合效果良好。

圖2 Q-W72組合理論頻率與經(jīng)驗頻率擬合

由此建立洪水峰量聯(lián)合分布概率模型，并繪于圖3，根據(jù)式(9)建立兩變量聯(lián)合重現(xiàn)期模型以及相應(yīng)等值線圖如圖4(圖中T取2～500 a)。

圖3 洪峰與72 h洪量聯(lián)合分布模型

圖4 洪峰與72 h洪量聯(lián)合重現(xiàn)期分布模型

推求設(shè)計洪水過程線前需要分別確定洪峰流量和洪量的設(shè)計值(表2)，從圖4中可知在選定的重現(xiàn)期下存在無數(shù)多種滿足條件的峰量組合。這里采用同頻率方法，即規(guī)定式(9)中u=v，如圖4(b)中所示。可確定所求設(shè)計值組合同頻率值，如式(11)，進(jìn)而通過邊緣分布函數(shù)求得對應(yīng)的洪峰和洪量設(shè)計值，如式(12)。

(11)

表2 不同重現(xiàn)期下的Q-W72設(shè)計值

(12)

2.3 調(diào)洪計算分析

調(diào)洪計算選取對水庫工程安全最為不利的“94·6”典型洪水過程。 “94·6”洪水是由強臺風(fēng)和冷空氣的共同作用下造成的，而臺風(fēng)雨是影響本流域主要的天氣系統(tǒng)。因此，選定“94·6”為典型洪水過程線是比較理想的。

圖5列出了“94·6”典型洪水過程以及基于兩種方法的洪水放大過程。由于峰量聯(lián)合放大法考慮了兩變量間的關(guān)聯(lián)性，基于聯(lián)合重現(xiàn)期推求設(shè)計洪水過程線，因此設(shè)計結(jié)果較傳統(tǒng)的同頻率法是偏安全的，從圖5中也可以明顯看出，在不同的設(shè)計重現(xiàn)期下，峰量聯(lián)合法的洪水過程線略大于傳統(tǒng)方法。對于整個洪水過程來說，防洪標(biāo)準(zhǔn)的制定不應(yīng)只針對某一個特征量，而在聯(lián)合重現(xiàn)期條件下的設(shè)計值組合是基于峰量相關(guān)性的一種組合，因此，推求的設(shè)計洪水過程線更符合水文現(xiàn)象的內(nèi)在規(guī)律，設(shè)計值組合更加合理。對比之下，目前以單變量同頻放大法推求的設(shè)計洪水過程線達(dá)不到防洪標(biāo)準(zhǔn)，即考慮多變量相關(guān)性下的設(shè)計洪水可以更準(zhǔn)確地反映洪水的整體過程。

表3列出了基于傳統(tǒng)同頻率法與峰量聯(lián)合同頻率法的水庫調(diào)洪計算結(jié)果，分別統(tǒng)計了在50、100、200、500 a一遇洪水下的水庫調(diào)洪最高庫水位和最大泄流量。從表3中可以看出，基于峰量同頻法的水庫調(diào)洪最高庫水位均高于傳統(tǒng)同頻率法的計算值，在水庫防洪安全設(shè)計中，采用兩變量同頻設(shè)計法不僅可以更全面地考慮整個洪水過程，并且從防洪設(shè)計安全性的角度考慮，采用新方法更為科學(xué)合理。

圖5 重現(xiàn)期下設(shè)計洪水過程線推求結(jié)果

表3 不同重現(xiàn)期下的水庫調(diào)洪計算成果(39.3 m起調(diào))

3 結(jié) 語

(1)基于單變量的設(shè)計方法只能分別控制各個特征變量的重現(xiàn)期等于防洪標(biāo)準(zhǔn)，忽略了洪水過程的整體性。本文采用G-H Copula函數(shù)構(gòu)造了P-Ⅲ型峰量分布的聯(lián)合分布模型，將設(shè)計洪水簡化成洪峰、洪量兩個特征量，考慮了洪水特征量之間的關(guān)聯(lián)性。

(2)在水庫防洪設(shè)計中，以兩變量方法推求的設(shè)計洪水過程線是偏安全的，計算所得調(diào)洪最高庫水位和最大泄流量也略大于單變量方法下的值，其中最高庫水位升高0.13～0.15 m；最大泄量增加30～60 m3/s。

(3)新方法在保持典型洪水形態(tài)的基礎(chǔ)上略優(yōu)于傳統(tǒng)計算方法，從安全性的角度考慮，采用新方法更為科學(xué)合理，為水庫防洪標(biāo)準(zhǔn)設(shè)定提供了新的研究思路。

□