李　帆,鄭　騫,張　磊

(1. 揚(yáng)州大學(xué)水利與能源動(dòng)力工程學(xué)院,江蘇 揚(yáng)州　225000;2. 衢州市水利局,浙江 衢州　324000; 3. 徐州市水利建筑設(shè)計(jì)研究院,江蘇 徐州　221116)

?

Copula熵方法及其在三變量洪水頻率計(jì)算中的應(yīng)用

李帆1,鄭騫2,張磊3

(1. 揚(yáng)州大學(xué)水利與能源動(dòng)力工程學(xué)院,江蘇 揚(yáng)州225000;2. 衢州市水利局,浙江 衢州324000; 3. 徐州市水利建筑設(shè)計(jì)研究院,江蘇 徐州221116)

為解決傳統(tǒng)Copula方法在進(jìn)行聯(lián)合概率分布擬合過程中要先進(jìn)行函數(shù)類型選擇的問題,將Copula函數(shù)和最大熵原理進(jìn)行耦合,通過求解具有最大熵的Copula方程,求得二維聯(lián)合分布函數(shù),即Copula熵方法。用求得的Copula函數(shù)對(duì)洪水事件的3個(gè)相關(guān)變量(洪峰流量、洪水總量和洪水歷時(shí))進(jìn)行兩兩配對(duì)的二維聯(lián)合分布擬合,并利用Gibbs采樣方法和Copula函數(shù)實(shí)現(xiàn)三變量洪水事件的隨機(jī)模擬。以淮河魯臺(tái)子水文站的實(shí)測(cè)洪水資料為研究對(duì)象,進(jìn)行實(shí)例分析,并通過擬合優(yōu)度的計(jì)算,證明Copula熵方法對(duì)多維相關(guān)變量概率擬合的有效性以及Gibbs采樣方法在三變量洪水事件模擬過程中的有效性。

洪水頻率;聯(lián)合分布函數(shù);Copula函數(shù);最大熵原理;Copula熵方法;淮河;魯臺(tái)子水文站

洪水事件的概率估計(jì)是防洪規(guī)劃、水利工程建設(shè)和水資源管理的重要依據(jù)。一個(gè)洪水事件一般可以由3個(gè)具有相關(guān)關(guān)系的隨機(jī)變量來描述,即洪峰流量(P)、洪水總量(V)和洪水歷時(shí)(D)。Copula函數(shù)方法是一種實(shí)現(xiàn)多變量洪水事件聯(lián)合分布擬合的有效途徑,并在水文領(lǐng)域中被廣泛應(yīng)用。郭生練等[1]對(duì)Copula函數(shù)在水文中的應(yīng)用進(jìn)行了綜述;宋松柏等[2]系統(tǒng)地探討了各種Copula函數(shù)結(jié)構(gòu)在西北典型流域中的應(yīng)用,分析了Copula函數(shù)在多變量水文計(jì)算領(lǐng)域的適用性及優(yōu)越性。Copula函數(shù)存在多種形式,其中二維Archimedean Copula函數(shù)應(yīng)用最廣泛。常用的Archimedean Copula函數(shù)包括Frank Copula、Gumbel Copula和Clayton Copula等類型,此外還有屬于橢圓Copula函數(shù)的Gaussian Copula等。針對(duì)具體問題選取合適的Copula函數(shù)進(jìn)行多維聯(lián)合概率分布擬合。Salvadori等[3]用Archimedean Copula函數(shù)建立了最高洪水位與洪水總量、洪峰流量與洪水總量、降水強(qiáng)度與降水歷時(shí)的3組二維聯(lián)合分布,分析了Gumbel Copula和Frank Copula函數(shù)在水文二維極值聯(lián)合分布構(gòu)建上的有效性。Zhang 等[4]用4種Archimedean Copula 函數(shù)構(gòu)造了洪峰與洪水總量、洪水總量與洪水歷時(shí)的2組二維聯(lián)合分布,并比較了他們的擬合優(yōu)度,認(rèn)為Gumbel Copula對(duì)所選河流洪水的擬合最好。肖義等[5]、李天元等[6]用Gumbel Copula函數(shù)構(gòu)造邊緣分布為P-Ⅲ型分布的兩變量聯(lián)合分布,用以描述洪峰和洪水總量。Renard等[7]探討了Gaussian Copula在洪水頻率分析上的應(yīng)用。李計(jì)等[8]比較了4種對(duì)稱的Archimedean Copula和5種非對(duì)稱的Archimedean Copula對(duì)干旱歷時(shí)、干旱烈度和烈度峰值3個(gè)干旱特征變量聯(lián)合分布的擬合效果。雖然Copula函數(shù)的應(yīng)用已十分廣泛,但在進(jìn)行Copula函數(shù)參數(shù)的估算之前,必須先要對(duì)Copula函數(shù)的類型進(jìn)行選擇,然后再進(jìn)行擬合分析。這種選型過程往往并不完全客觀,從而造成計(jì)算結(jié)果不準(zhǔn)確。為了避免這種由于選型而造成的偏差,本文將最大熵原理和Copula函數(shù)理論進(jìn)行耦合,即Copula熵方法,用以求得聯(lián)合概率分布密度函數(shù)。另外,傳統(tǒng)的Copula方程都是兩變量的聯(lián)合分布方程,在向三維變量拓展時(shí)往往比較困難,計(jì)算也比較復(fù)雜。筆者基于蒙特卡洛算法,利用Gibbs采樣方法和Copula函數(shù),生成大量洪水事件,從而實(shí)現(xiàn)對(duì)三變量的洪水事件的隨機(jī)模擬,進(jìn)而分析洪水事件的發(fā)生概率，并以淮河魯臺(tái)子水文站實(shí)測(cè)資料為對(duì)象,進(jìn)行了示例研究。

1　理論與方法

1.1Copula函數(shù)

由Sklar定理可得二維Copula函數(shù):

(1)

式中:F(x1,x2)——聯(lián)合概率分布函數(shù);F1(x1)=u、F2(x2)=v——邊緣分布函數(shù)。

根據(jù)Sklar定理,若邊緣分布函數(shù)是連續(xù)的,則Copula函數(shù)C(u,v)唯一。

1.2最大熵原理與Copula熵

根據(jù)熵理論,對(duì)于一個(gè)連續(xù)的隨機(jī)變量X∈[a1,a2],若概率密度函數(shù)為f(x),則其熵可定義為

(2)

根據(jù)最大熵原理(POME),當(dāng)H取得最大值時(shí),人為假定最小,從而可以得到最合乎自然、最為超然、偏差最小的f(x)[9]。對(duì)于Copula函數(shù),其對(duì)應(yīng)的Copula熵為

(3)

式中:c(u,v)——二維Copula概率密度函數(shù)。

根據(jù)POME,式(3)中使得W最大的c(u,v)即為最適合的聯(lián)合概率密度函數(shù)。

1.3求解Copula概率密度函數(shù)

為計(jì)算最優(yōu)的c(u,v),可以利用拉格朗日乘子法對(duì)式(3)進(jìn)行求解。根據(jù)Copula概率密度函數(shù)的性質(zhì),可列出以下約束條件:

(4)

(5)

(6)

式中:r——r階矩;ρ——Spearman秩相關(guān)系數(shù)。

Chu[10]證明,使得式(3)取得最大值的c(u,v)的表達(dá)式為

(7)

其中

式中:α、β、λ——待求系數(shù);n——矩的階數(shù)。

式(7)沒有解析解,可以使用牛頓迭代法計(jì)算待求系數(shù)[11],進(jìn)而求得Copula概率密度函數(shù)的表達(dá)式。

1.4Gibbs采樣方法生成三變量洪水事件

Copula熵方法在理論上可以直接從二維擴(kuò)展到三維,得到三變量Copula概率密度函數(shù),但計(jì)算量較大,比較耗時(shí)。筆者采用Gibbs采樣方法[12]和二維Copula函數(shù)進(jìn)行三維隨機(jī)變量模擬,這種方法原理簡(jiǎn)單,計(jì)算速度快,并已在海浪波況的隨機(jī)模擬中成功運(yùn)用[13-14]。

令X、Y、Z為洪峰流量、洪水總量和洪水歷時(shí)的隨機(jī)變量,x、y、z分別表示洪峰流量、洪水總量、洪水歷時(shí)的累計(jì)頻率,cXY和cXZ分別表示洪峰流量與洪水總量的Copula密度函數(shù)和洪峰流量與洪水歷時(shí)的Copula密度函數(shù)。隨機(jī)生成4個(gè)0～1之間的隨機(jī)數(shù)(a,b,c,d),令y0=a,并使得

(8)

由式(8)可求得x1,把x1代入

(9)

(10)

可分別求得y1和z1。重復(fù)以上步驟可以得到任意多組三維模擬值(x,y,z),然后利用邊緣概率分布函數(shù)的反函數(shù)將(x,y,z)轉(zhuǎn)換成其所對(duì)應(yīng)的物理值,得到任意多個(gè)三變量洪水事件,從而利用模擬的洪水事件進(jìn)行概率的估計(jì)和分析。

2　結(jié)果與分析

以淮河干流魯臺(tái)子水文站多年(1954—1998年)洪水觀測(cè)資料為研究對(duì)象,驗(yàn)證Copula熵方法的有效性。魯臺(tái)子水文站位于淮河中游,控制面積88 630 km2,多年平均流量693 m3/s,資料長(zhǎng)度45 a。為了充分利用觀測(cè)數(shù)據(jù),采用超閾值法選取洪水事件。流量超過2 000 m3/s的時(shí)刻,定義為洪水事件的起始時(shí)間,流量低于2 000 m3/s的時(shí)刻,定義為洪水事件的結(jié)束時(shí)間。為了保證洪水事件的獨(dú)立性,規(guī)定2個(gè)洪水事件之間的間隔不小于3 d,否則,兩場(chǎng)洪水將合并為一場(chǎng)洪水。根據(jù)該定義方法,共從原始觀測(cè)資料中提取出56場(chǎng)洪水事件。

表1　各分布函數(shù)的擬合優(yōu)度Table 1　Goodness-of-fit for different distribution functions

2.1邊緣分布擬合

分別用P-Ⅲ分布、廣義帕累托分布(Generalized Pareto Distribution,GPD)、廣義極值分布(Generalized Extreme Value distribution,GEV)、Gamma分布和指數(shù)分布對(duì)觀測(cè)值進(jìn)行擬合,并通過計(jì)算歸一化均方差(Normalized Mean Square Error,NMSE)與經(jīng)驗(yàn)分布進(jìn)行比較。結(jié)果表明P-Ⅲ分布的擬合效果最好(表1)。圖1為P-Ⅲ理論分布曲線與經(jīng)驗(yàn)點(diǎn)據(jù)的擬合情況。

圖1　洪水事件三要素的理論分布曲線與經(jīng)驗(yàn)點(diǎn)據(jù)Fig. 1　Theoretical distribution curves and empirical cumulative distribution functions of flood variables

2.2Copula密度函數(shù)計(jì)算

按照1.3中方法進(jìn)行系數(shù)估計(jì),當(dāng)式(7)中的n分別取2和3時(shí),計(jì)算得到的各系數(shù)見表2。當(dāng)n=3時(shí),系數(shù)α3和β3非常小,表明ε(u,v)只需展開到n=2即可。于是可得到3個(gè)二元Copula密度函數(shù):

(11)

則Copula累積概率分布函數(shù)可表示為

(12)

表2　兩變量Copula函數(shù)系數(shù)Table 2　Parameters of bivariate copula functions

利用得到的累計(jì)概率分布函數(shù)及聯(lián)合重現(xiàn)期公式得到圖2中的聯(lián)合重現(xiàn)期等值線(重現(xiàn)期單位為a)。聯(lián)合重現(xiàn)期表示其中一個(gè)變量超過某一定值的概率,其公式見文獻(xiàn)[3],本文平均每年發(fā)生的洪水次數(shù)取為56/45。

2.3三變量洪水事件模擬

由于使用了超閾值方法進(jìn)行洪水事件的選取,所以每年發(fā)生洪水事件的次數(shù)也需要進(jìn)行估計(jì)。筆者使用經(jīng)驗(yàn)頻率來估計(jì)每一年洪水的發(fā)生次數(shù)。在45 a(1954—1998)間,沒有發(fā)生超閾值洪水的年數(shù)是10,發(fā)生1次的有14 a,發(fā)生2次的有17 a,發(fā)生3次的有6 a,即P(n=0)=10/45,P(n=1)=17/45,P(n=2)=12/45,P(n=3)=6/45。結(jié)合1.3節(jié)中的Gibbs采樣方法,得到1 000 a的模擬洪水事件(圖2)。

圖2　1 000 a的洪水模擬事件與觀測(cè)值以及兩變量聯(lián)合重現(xiàn)期等值線Fig. 2　Observations and simulated flood events with a 1000-year return period and bivariate joint return period contours

2.4洪水事件的擬合優(yōu)度評(píng)估

為定量估計(jì)Copula熵方法與Gibbs采樣方法的有效性,分別計(jì)算觀測(cè)值與模擬值的經(jīng)驗(yàn)頻率,并進(jìn)行比較。其中,三變量經(jīng)驗(yàn)頻率定義為[15-16]

(13)

式中:K——滿足的樣本個(gè)數(shù),m——總的樣本個(gè)數(shù)。

表3　卡方檢驗(yàn)與K-S檢驗(yàn)結(jié)果Table 3　Results of Chi squared test and K-S test

兩變量經(jīng)驗(yàn)頻率的定義與式(13)相似。利用卡方檢驗(yàn)與K-S(Kolmogorov-Sminov)檢驗(yàn)來進(jìn)行模擬值與實(shí)測(cè)值的比較。表3給出了2種檢驗(yàn)方法的計(jì)算結(jié)果,如果卡方統(tǒng)計(jì)值小于臨界值,K-S檢驗(yàn)的P值大于0.05,則表明在5%的顯著水平下通過假設(shè)檢驗(yàn),說明擬合情況較好,反之則不通過假設(shè)檢驗(yàn),說明擬合情況無法達(dá)到要求。計(jì)算結(jié)果表明,模擬值在5%的顯著性水平下與樣本的擬合較好,能夠反映實(shí)測(cè)洪水的統(tǒng)計(jì)特征。

3　結(jié)論與展望

a. Copula熵方法與Copula方法類似,可以通過任意邊緣分布和變量間的相關(guān)性構(gòu)建二維的聯(lián)合概率分布函數(shù),與單純的Copula方法比較,Copula熵方法省去了Copula方程類型選擇這一過程,計(jì)算效率更高。

b. 結(jié)合Gibbs采樣方法和二維Copula函數(shù),可以利用二維聯(lián)合概率分布函數(shù)和條件概率函數(shù)進(jìn)行三維隨機(jī)變量的模擬,從而實(shí)現(xiàn)對(duì)三變量洪水事件的概率模擬。以淮河干流魯臺(tái)子水文站為研究對(duì)象,應(yīng)用上述方法結(jié)合1954—1998年間的洪水實(shí)測(cè)資料,對(duì)洪水事件進(jìn)行了概率模擬,模擬結(jié)果和檢驗(yàn)結(jié)果顯示,利用Copula熵方法對(duì)二維聯(lián)合分布進(jìn)行擬合適用且有效。此外,Gibbs采樣方法可以實(shí)現(xiàn)對(duì)三變量洪水事件的隨機(jī)模擬,模擬結(jié)果通過了擬合優(yōu)度檢驗(yàn)。該方法可以為該地區(qū)的水利工程規(guī)劃、洪水風(fēng)險(xiǎn)評(píng)估以及水資源合理利用提供參考依據(jù)。

c. Copula熵方法在水文領(lǐng)域中的應(yīng)用還處于起步階段,相關(guān)的研究還比較少,有必要對(duì)該方法進(jìn)行更為深入的研究。本文僅對(duì)Copula熵方法進(jìn)行了簡(jiǎn)單的介紹,并用一個(gè)實(shí)例進(jìn)行了驗(yàn)證。Copula熵方法在其他地區(qū)及領(lǐng)域的有效性以及與傳統(tǒng)Copula函數(shù)擬合效果的比較還需要進(jìn)一步研究與分析。此外,Copula在求解聯(lián)合分布方程的過程中,不同約束條件對(duì)擬合結(jié)果的影響等問題也需要進(jìn)一步討論。

[1] 郭生練,閆寶偉,肖義,等. Copula 函數(shù)在多變量水文分析計(jì)算中的應(yīng)用及研究進(jìn)展[J]. 水文, 2008,28(3): 1-7. (GUO Shenglian, YAN Baowei, XIAO Yi, et al. Multivariate hydrological analysis and estimation[J]. Journal of China Hydrology, 2008, 28(3): 1-7. (in Chinese))

[2] 宋松柏, 蔡煥杰, 金菊良,等. Copulas函數(shù)及其在水文中的應(yīng)用[M]. 北京: 科學(xué)出版社, 2012.

[3] SALVADORI G, DE MICHELE C. Frequency analysis via Copulas: theoretical aspects and applications to hydrological events [J]. Water Resources Research, 2004, 40(12):1-17.

[4] ZHANG L, SINGH V P. Bivariate flood frequency analysis using the Copula method [J]. Journal of Hydrologic Engineering, 2006, 11(2): 150-164.

[5] 肖義, 郭生練, 熊立華, 等. 一種新的洪水過程隨機(jī)模擬方法研究[J]. 四川大學(xué)學(xué)報(bào)(工程科學(xué)版), 2007, 39(2): 55-60. (XIAO Yi, GUO Shenglian, XIONG Lihua, et al. A new random simulation method for constructing synthetic flood hydrographs[J]. Journal of Sichuan University(Engineering Science Edition), 2007, 39(2): 55-60. (in Chinese))

[6] 李天元, 郭生練, 劉章君, 等. 基于峰量聯(lián)合分布推求設(shè)計(jì)洪水[J]. 水利學(xué)報(bào), 2014, 45(3):269-276. (LI Tianyuan, GUO Shenglian, LIU Zhangjun et al. Design flood estimation based on bivariate joint distribution of flood peak and volume [J]. Journal of Hydraulic Engineering, 2014, 45(3): 269-276. (in Chinese))

[7] RENARD B, LANG M. Use of a Gaussian Copula for multivariate extreme value analysis: some case studies in hydrology [J]. Advances in Water Resources, 2007,30: 897-912.

[8] 李計(jì),李毅,宋松柏,等. 基于Copulas 函數(shù)的多維干旱變量聯(lián)合分布[J]. 自然資源學(xué)報(bào), 2013,28(2): 312-320. (LI Ji, LI Yi, SONG Songbai, et al. Multivariate joint distributions of drought variables based on Copulas function[J].Journal of Natural Resources, 2013, 28(2): 312-320. (in Chinese))

[9] 王棟,朱元甡. 最大熵原理在水文水資源科學(xué)中的應(yīng)用[J]. 水科學(xué)進(jìn)展,2001,12(3):424-430. (WANG Dong, ZHU Yuansheng. Principle of maximum entropy and its application in hydrology and water resources[J]. Advances in Water Science, 2001, 12(3): 424-430. (in Chinese))

[10] CHU B. Recovering Copulas from limited information and an application to asset allocation[J]. Journal of Banking and Finance, 2011, 35: 1824-42.

[11] AGHAKOUCHAK A. Entropy-Copula in hydrology and climatology [J]. Journal of Hydrometeorology, 2014,15: 1-13.

[12] GEMAN S, GEMAN D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1984,20: 25-62.

[13] LI F, VAN GELDER P H A J M, RANASINGHE R, CALLAGHAN D P, JONGEJAN R B. Probabilistic modelling of extreme storms along the Dutch coast[J]. Coastal Engineering, 2014,86(4): 1-13.

[14] CALLAGHAN D P, NIELSEN P, SHORT A, et al. Statistical simulation of wave climate and extreme beach erosion[J]. Coastal Engineering, 2008, 55(5): 375-390.

[15] SALVADORI G, DE MICHELE C. Frequency analysis via Copulas: theoretical aspects and applications to hydrological events[J]. Water Resources Research, 2004, 40: 1-17.

[16] ZHANG L, SINGH V P. Trivariate flood frequency analysis using the Gumbel-Hougaard Copula[J]. Journal of Hydrological Engineering, 2007,12: 431-439.

Copula entropy method and its application to trivariate flood frequency calculation

LI Fan1, ZHENG Qian2, ZHANG Lei3

(1.SchoolofHydraulic,Energy,andPowerEngineering,YangzhouUniversity,Yangzhou225000,China;2.QuzhouWaterConservancyBureau,Quzhou324000,China;3.XuzhouHydraulicConstructionandDesignResearchInstitute,Xuzhou221116,China)

In order to avoid selection of the type of functions when using traditional copula methods to fit the joint probability distribution, the copula function and the maximum entropy principle were jointly used (in a method also called the copula entropy method), and the two-dimensional joint distribution functions of flood variates were obtained by solving the copula function with the maximum entropy. With the solved copula function, the two-dimensional joint distribution functions of flood variates (peak discharge, flood volume, and flood duration) were mutually constructed. The Gibbs sampling method and copula functions were used to stochastically simulate trivariate flood events. A case study was conducted using the observed flood data from the Lutaizi Hydrological Station on the Huaihe River. A goodness-of-fit analysis shows that the copula entropy method is effective in fitting multivariate probability distributions and the Gibbs sampling method is effective in simulating trivariate flood events.

flood frequency; joint distribution function; copula function; maximum entropy principle; copula entropy method; Huaihe River; Lutaizi Hydrological Station

10.3876/j.issn.1000-1980.2016.05.011

2015-12-07

江蘇省高校自然科學(xué)研究項(xiàng)目(15KJD170003);揚(yáng)州大學(xué)科技創(chuàng)新培育基金(2015CXJ027)

李帆(1982—),男,遼寧營(yíng)口人,講師,博士,主要從事水文統(tǒng)計(jì)研究。E-mail:lifan@yzu.edu.cn

P333.2

A

1000-1980(2016)05-0443-06