李　帆,鄭　騫,張　磊

(1. 揚州大學水利與能源動力工程學院,江蘇 揚州　225000;2. 衢州市水利局,浙江 衢州　324000; 3. 徐州市水利建筑設計研究院,江蘇 徐州　221116)

?

Copula熵方法及其在三變量洪水頻率計算中的應用

李帆1,鄭騫2,張磊3

(1. 揚州大學水利與能源動力工程學院,江蘇 揚州225000;2. 衢州市水利局,浙江 衢州324000; 3. 徐州市水利建筑設計研究院,江蘇 徐州221116)

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

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

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

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ù)熵理論,對于一個連續(xù)的隨機變量X∈[a1,a2],若概率密度函數(shù)為f(x),則其熵可定義為

(2)

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

(3)

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

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

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

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

(4)

(5)

(6)

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

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

(7)

其中

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

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

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

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

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

(8)

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

(9)

(10)

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

2　結果與分析

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

表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ù)分布對觀測值進行擬合,并通過計算歸一化均方差(Normalized Mean Square Error,NMSE)與經(jīng)驗分布進行比較。結果表明P-Ⅲ分布的擬合效果最好(表1)。圖1為P-Ⅲ理論分布曲線與經(jīng)驗點據(jù)的擬合情況。

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

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

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

(11)

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

(12)

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

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

2.3三變量洪水事件模擬

由于使用了超閾值方法進行洪水事件的選取,所以每年發(fā)生洪水事件的次數(shù)也需要進行估計。筆者使用經(jīng)驗頻率來估計每一年洪水的發(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。結合1.3節(jié)中的Gibbs采樣方法,得到1 000 a的模擬洪水事件(圖2)。

圖2　1 000 a的洪水模擬事件與觀測值以及兩變量聯(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)度評估

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

(13)

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

表3　卡方檢驗與K-S檢驗結果Table 3　Results of Chi squared test and K-S test

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

3　結論與展望

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

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

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

[1] 郭生練,閆寶偉,肖義,等. Copula 函數(shù)在多變量水文分析計算中的應用及研究進展[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ù)及其在水文中的應用[M]. 北京: 科學出版社, 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]. 四川大學學報(工程科學版), 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)合分布推求設計洪水[J]. 水利學報, 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] 李計,李毅,宋松柏,等. 基于Copulas 函數(shù)的多維干旱變量聯(lián)合分布[J]. 自然資源學報, 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] 王棟,朱元甡. 最大熵原理在水文水資源科學中的應用[J]. 水科學進展,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

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

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

P333.2

A

1000-1980(2016)05-0443-06