基于Copula函數(shù)的水文變量條件組合分析

劉和昌，梁忠民，姚 軼，唐甜甜

（河海大學(xué)水文與水資源學(xué)院，江蘇 南京 210098）

0 引 言

目前，在水文計(jì)算中采用的多是特定組合的分析方法，如最不利組合、同頻率組合等。理論上，進(jìn)行多變量的組合分析計(jì)算，必須給出變量間的聯(lián)合概率分布函數(shù)。目前，應(yīng)用最廣泛的是多元正態(tài)聯(lián)合分布；但其針對(duì)的是各變量邊際分布為正態(tài)分布的情形。非正態(tài)分布，原則上需要對(duì)變量進(jìn)行正態(tài)化變換，變換過(guò)程復(fù)雜且會(huì)使一些信息失真，限制了其使用范圍。

近年來(lái)興起的基于Copula函數(shù)的多維聯(lián)合分布理論[1]，將變量間的聯(lián)合分布分解為變量間的相關(guān)性結(jié)構(gòu)和變量的邊緣分布，為構(gòu)建多變量聯(lián)合分布提供了有力的數(shù)學(xué)分析工具。由于Copula函數(shù)對(duì)變量的邊際分布線型沒(méi)有任何限制，而且聯(lián)合分布函數(shù)可以解析表達(dá)，所以在水文領(lǐng)域中得到了較廣泛的研究和應(yīng)用。Favre等[2]首次利用Copula函數(shù)研究了洪水地區(qū)組成及峰量組成的計(jì)算問(wèn)題；Poulin等[3]針對(duì)多變量頻率分析中的尾部相關(guān)，研究了Copula函數(shù)的選型問(wèn)題；肖義等[4]應(yīng)用Copula計(jì)算洪水峰、量聯(lián)合分布，以并集聯(lián)合分布概率為設(shè)計(jì)標(biāo)準(zhǔn)作為控制，由此對(duì)典型洪水過(guò)程線進(jìn)行放大，獲得設(shè)計(jì)洪水過(guò)程線；郭生練等[5]對(duì)copula函數(shù)在水文中研究與應(yīng)用的進(jìn)展進(jìn)行了評(píng)述和總結(jié)；張娜等[6]通過(guò)Copula函數(shù)分別構(gòu)建最大日雨量與最大時(shí)段雨量的聯(lián)合分布，以推求設(shè)計(jì)暴雨過(guò)程；閆寶偉等[7]根據(jù)Copula聯(lián)合分布的條件期望組合概念，研究了設(shè)計(jì)洪水的地區(qū)組成計(jì)算方法；梁忠民等[8]采用三變量Copula函數(shù)，分析了三峽水庫(kù)預(yù)泄對(duì)鄱陽(yáng)湖防洪的影響。

本文基于Copula函數(shù)推導(dǎo)了條件最可能組合的計(jì)算公式，并對(duì)錢(qián)塘江水系的4個(gè)站點(diǎn)洪水頻率進(jìn)行了分析計(jì)算，提供了不同設(shè)計(jì)標(biāo)準(zhǔn)條件下最可能組合的峰、量設(shè)計(jì)值及其90%的置信區(qū)間，并與同頻率組合及條件期望組合的結(jié)果進(jìn)行了對(duì)比分析。

1 Copula函數(shù)和條件組合

1.1 Copula函數(shù)

Copula函數(shù)是定義域?yàn)閇0，1]均勻分布的多維聯(lián)合分布函數(shù)，二維Copula函數(shù)為[1]

式中，C為Copula函數(shù)；θ為Copula的參數(shù)；u=FX（x），v=FY（y）分別為隨機(jī)變量X和Y的邊緣分布。相應(yīng)的聯(lián)合概率密度函數(shù)為

式中，c（u，v）為Copula函數(shù)的密度函數(shù)；fX（x）和fY（y）分別為隨機(jī)變量X和Y的概率密度函數(shù)。

在水文領(lǐng)域，兩個(gè)變量單參數(shù)Archimedean Copula函數(shù)簇應(yīng)用最為廣泛，Gumbel-Hougaard Copula函數(shù)即為其中一種。Gumbel-Hougaard Copula函數(shù)適于變量間存在正相關(guān)的情形，并且能描述變量間的上尾相關(guān)性[3]，故本文選用Gumbel-Hougaard Copula作為聯(lián)合分布函數(shù)，其表達(dá)式為

其密度函數(shù)為

函數(shù)參數(shù)θ與Kendall秩相關(guān)τ系數(shù)間關(guān)系如下

采用離差平方和準(zhǔn)則 （OLS）來(lái)評(píng)價(jià)Copula方法的擬合度

式中，Pei和Pi分別為經(jīng)驗(yàn)頻率和理論頻率。

1.2 邊緣分布

隨機(jī)變量X和隨機(jī)變量Y均采用P-III分布，則其概率密度函數(shù)

式中，Γ（α）為伽瑪函數(shù)；α、β和a0分別為形狀、尺度和位置參數(shù)。

隨機(jī)變量X和隨機(jī)變量Y的分布函數(shù)

式中，x0和y0為大于0的任一值；fX（x）和fY（x）分別變量X和變量Y的概率密度函數(shù)。

1.3 條件組合

1.3.1 條件最可能組合

當(dāng)變量X取設(shè)計(jì)值xp時(shí)，變量Y所對(duì)應(yīng)的值并非唯一，可大可小，只是出現(xiàn)不同值的概率不一樣，即存在一個(gè)條件概率分布

此時(shí)，F(xiàn)Y|X（y）的密度函數(shù)fY|X（y）=c（u，v）fY（y），而當(dāng)fY|X（y）取最大值時(shí)對(duì)應(yīng)的yM值所形成組合（xp，yM），即為條件最可能組合。fY|X（y）為y的一元函數(shù)，將密度函數(shù)fY|X（y）對(duì)y求導(dǎo)，得

式中，c1=?c（u，v）/?v，c、u、v定義同上。 再令dfY|Xdy=0，整理得:

求解式（12），可得x取xp時(shí)，y的最可能取值yM。

1.3.2 條件期望組合

本節(jié)中，在給定變量X的值為xp時(shí)，y的期望值E（y|xp）對(duì)應(yīng)的組合（xp，E（y|xp））即為條件期望組合[7]。其中，E（y|xp）的計(jì)算如下:

式中，fY|X（y）=c（u，v）fY（y），為FY|X（y）的密度函數(shù)；（v）為v=FY（y）的反函數(shù)。

變量X取值為xp時(shí)，若令FY|X（y）分別等于α/2和1－α/2，求得的Y值分別為yl和yu；則已知xp時(shí)y的1－α的置信區(qū)間為[yl，yu]。

2 實(shí)例分析

以錢(qián)塘江湖南鎮(zhèn)、黃壇口、新安江和富春江4個(gè)站點(diǎn)40年的洪水資料為例，對(duì)本文所提方法進(jìn)行應(yīng)用分析，其中湖南鎮(zhèn)、黃壇口、新安江位于錢(qián)塘江支流上，富春江位于干流上。洪峰和時(shí)段洪量均按服從P-III型曲線分布，采用線性矩法估計(jì)分布參數(shù)，各站點(diǎn)年最大洪峰 （m3/s）和時(shí)段洪量 （億m3）統(tǒng)計(jì)參數(shù)估計(jì)結(jié)果見(jiàn)表1；洪峰和時(shí)段洪量間的Kendall秩相關(guān)系數(shù)τ由式（5）計(jì)算，其結(jié)果及對(duì)應(yīng)的OLS值見(jiàn)表2。

表1 各站點(diǎn)統(tǒng)計(jì)

表2 各站點(diǎn)τ值和OLS值

以富春江站為例，通過(guò)給定不同洪峰值，由式（12）和式 （10）可分別求得對(duì)應(yīng)的7天洪量條件最可能值和90%的置信區(qū)間。為了對(duì)比分析，也提供了條件期望值及同頻率值的估計(jì)結(jié)果 （見(jiàn)圖1，2及表3，4）。

圖1 富春江站以洪峰為條件的7天洪量估計(jì)結(jié)果

圖2 富春江站以7天洪量為條件的洪峰估計(jì)結(jié)果

表3 富春江站以洪峰為條件的7天洪量的各種估計(jì)值

表4 富春江站以7天洪量為條件的洪峰的各種估計(jì)值

考慮洪峰與7天洪量之間的相關(guān)性，條件最可能組合是給定洪峰值 （或7天洪量）時(shí)，7天洪量（或洪峰）取值的最可能情況，條件期望組合則代表的是其取值的平均情況，兩種組合方法均具有一定的統(tǒng)計(jì)基礎(chǔ)；一般認(rèn)為，目前常用的同頻率組合法是一種更偏安全或保守的方法。

由圖1～2可知，富春江站，在給定洪峰值情況下，以設(shè)計(jì)頻率等于17%為界 （對(duì)應(yīng)約6年一遇重現(xiàn)期），當(dāng)其小于17%時(shí)，7天洪量的條件最可能值介于同頻率值與條件期望值之間，同頻率值最大；在給定7天洪量情況下，當(dāng)設(shè)計(jì)頻率小于25%時(shí)（對(duì)應(yīng)4年一遇重現(xiàn)期），洪峰的條件最可能值位于同頻率值和條件期望值之間，同頻率值最大。因此，對(duì)富春江站而言，當(dāng)設(shè)計(jì)標(biāo)準(zhǔn)超過(guò)6年一遇時(shí)，按照最可能組合條件計(jì)算的峰、量設(shè)計(jì)值，其數(shù)值介于按同頻率組合與條件期望組合的估值之間，同頻率組合的估值最大。

湖南鎮(zhèn)、黃壇口和新安江3個(gè)站，也作了上述分析計(jì)算，結(jié)果類似。其中，湖南鎮(zhèn)的設(shè)計(jì)頻率臨界值為25%，黃壇口的為14%，新安江的為29%，即當(dāng)其設(shè)計(jì)頻率小于該值時(shí)，按照最可能組合條件計(jì)算的峰、量設(shè)計(jì)值，其數(shù)值介于按同頻率組合與條件期望組合的估值之間；另外，估計(jì)結(jié)果也表明，同頻率組合的估值均是最大的。

因此，可以初步認(rèn)為，對(duì)于錢(qián)塘江水系，當(dāng)設(shè)計(jì)頻率小于14%時(shí) （對(duì)應(yīng)于8年一遇以上的重現(xiàn)期），同頻率組合方式計(jì)算的設(shè)計(jì)洪水大于條件期望組合和條件最可能組合的結(jié)果。對(duì)于水庫(kù)工程而言，其設(shè)計(jì)標(biāo)準(zhǔn)一般都是較高的，如百年一遇、千年一遇等；所以按照同頻率組合的設(shè)計(jì)成果一般都是偏于安全的。本文的研究表明，最可能組合計(jì)算的設(shè)計(jì)值介于同頻率組合與期望組合的設(shè)計(jì)值之間。

3 結(jié) 論

（1）在聯(lián)合分布的基礎(chǔ)上，推導(dǎo)了條件最可能組合的計(jì)算公式，為計(jì)算峰、量條件最可能設(shè)計(jì)值提供了基礎(chǔ)。

（2）對(duì)錢(qián)塘江水系4個(gè)站點(diǎn)的洪水頻率分析結(jié)果表明，按照最可能組合條件計(jì)算的峰、量設(shè)計(jì)值，其數(shù)值介于按同頻率組合與按條件期望組合的估值之間。由此也證明了，目前的同頻率法設(shè)計(jì)洪水計(jì)算成果一般是偏于安全的。

（3）條件最可能組合代表的是變量間的一種最可能遭遇出現(xiàn)的方式，較適合于水文事件的組合情況，但本文的結(jié)果僅是初步的，值得作進(jìn)一步研究。

［1］NELSON R B.An introduction to Copulas［M］.New York:Springer,1999.

［2］FAVRE A-C,ADLOUNI S E,PERRAULT L,et al.Multivariate hydrological frequency analysis using Copulas［J］.Water Resource Research,2004,40（1）,W01101.

［3］POULIN A,HUARD D,FAVRE A C,et al.Importance of tail dependence in bivariate frequency analysis［J］.Journal of Hydrologic Engineering,2007,12（4）:394-403.

［4］肖義，郭生練，劉攀，等.基于Copula函數(shù)的設(shè)計(jì)洪水過(guò)程線方法［J］.武漢大學(xué)學(xué)報(bào): 工學(xué)版，2007，40（4）:13-17.

［5］郭生練，閆寶偉，肖義，等.Copula函數(shù)在多變量水文分析計(jì)算中的應(yīng)用及研究進(jìn)展［J］.水文，2008，28（3）:1-7.

［6］張娜，郭生練，肖義，等.基于聯(lián)合分布的設(shè)計(jì)暴雨方法［J］.水力發(fā)電，2008，34（1）:18-21.

［7］閆寶偉，郭生練，郭靖，等.基于Copula函數(shù)的設(shè)計(jì)洪水地區(qū)組成研究［J］.水力發(fā)電學(xué)報(bào)，2010，29（6）:60-65.

［8］梁忠民，郭彥，胡義明，等.基于Copula函數(shù)的三峽水庫(kù)預(yù)泄對(duì)鄱陽(yáng)湖防洪影響分析 ［J］.水科學(xué)進(jìn)展，2012，23（4）:485-492.