基于Copula函數(shù)的葉爾羌河流域洪水要素聯(lián)合分布研究

何 英 ，彭 亮，鄭淑文，馬林瀟

(新疆農(nóng)業(yè)大學水利與土木工程學院，烏魯木齊 830052)

近年來，用Copula函數(shù)進行洪水過程分析成為水文計算領域的一個研究熱點[1]，楊衛(wèi)等選用Gumbel、Clayton和Frank 3種Copula函數(shù)建立降雨-洪量極值聯(lián)合分布模型[2]；范嘉煒等基于Copula函數(shù)分析了潖江河大廟峽流域洪峰流量與洪水歷時的聯(lián)合頻率分布特征[3]；侯蕓蕓以陜北地區(qū)洪水資料為研究對象，運用Copula函數(shù)建立了洪水特征變量的聯(lián)合概率分布和條件概率分布模型[4]；張冬冬等應用Archimedean Copula函數(shù)探討洪水多要素的聯(lián)合概率分布和條件概率分布[5]；高玉琴等應用Gumbel-Hougaard Copula函數(shù)進行秦淮河流域洪水風險分析[6]；李天元等采用金沙江屏山站和岷江高場站的日徑流資料，探討了雙參數(shù)Archimedean Copula函數(shù)在洪水聯(lián)合分布中的應用[6]；林嫻等基于Copula函數(shù)原理，利用武江流域實測水文資料，構建了流域組合變量Copula概率分布模型，分析了洪峰與洪量、洪量與洪水歷時、洪峰與洪水歷時的聯(lián)合概率分布[7]；李大鳴等采用阿基米德Gumbel-Hougaard Copula函數(shù)，構建了入庫洪峰與洪量的兩變量聯(lián)合分布模型，對桃林口水庫進行防洪風險分析[8]，姚瑞虎等基于Copula函數(shù)論述并推求了四種洪水重現(xiàn)期的理念及其計算，并以岷江支流馬邊河馬邊水文站為例，對比分析各重現(xiàn)期之間的關系[9]。

本文以葉爾羌河流域為例，利用卡群水文站實測水文資料，基于洪峰和洪量均服從P-Ⅲ型分布，應用Archimedean Copula函數(shù)，探討葉爾羌河流域洪峰、洪量聯(lián)合概率分布模型，以期為水利工程規(guī)劃設計和風險評估提供科學依據(jù)。

1 數(shù)據(jù)與方法

1.1 研究區(qū)域與數(shù)據(jù)來源

1.1.1 研究區(qū)域

葉爾羌，維吾爾語中意為“土地寬廣的地方”。葉爾羌河是塔里木河四源之一，位于中國最西部的新疆維吾爾自治區(qū)內(nèi)。洶涌的山區(qū)急流出了昆侖峽谷后向北流，形成許多分支，散布在沖積扇上，灌溉著葉爾羌綠洲；葉爾羌綠洲是新疆最大的綠洲之一，流出綠洲后的葉爾羌河繞過塔克拉瑪干沙漠西緣，流向東北，在阿克蘇綠洲南部匯集喀什噶爾河、阿克蘇河及和田河，形成塔里木河，是喀什地區(qū)的第一大河流。河流全長1 179 km，水源一是來自喬戈里峰的冰雪融水；二是河床西岸巖層中涌出的泉水；三是雨水，多年平均年徑流量75.71 億m3，年均向塔里木河輸水1.7 億m3，灌溉塔什庫爾干、葉城、澤普、莎車、麥蓋提、巴楚6個縣和農(nóng)三師10個團場共28.889 萬hm2耕地。葉爾羌河流域面積為10.8 萬km2，屬干旱性內(nèi)陸氣候區(qū)，降水少，蒸發(fā)強烈，多年平均降水量僅為55 mm左右，多年平均蒸發(fā)量為2 400 mm左右。區(qū)內(nèi)干旱、沙暴、洪水等災害發(fā)生較頻繁，自然生態(tài)環(huán)境脆弱，流域地理位置示意圖見圖1。

圖1 葉爾羌河流域地理位置示意圖Fig.1 Geographical location schematic diagram of the Yeerqiang River Basin

1.1.2 數(shù)據(jù)來源及資料代表性分析

本文研究主要采用葉爾羌河流域卡群水文站1954-2011年共58 a 日流量資料，采用年最大值選樣法(AM)，選取該水文站的年最大洪水所對應的洪峰流量、1日洪水總量為研究變量，運用Coupla函數(shù)對葉爾羌河流域卡群水文站站址處洪水要素進行聯(lián)合分布研究。

繪制卡群水文站洪峰流量模比系數(shù)差積曲線見圖2。由圖2可以看出，卡群水文站實測洪水系列中包含了較完整的漲落水過程，所以，卡群水文站58 a的最大洪峰流量樣本系列包括有豐、平、枯各種年份，具有較好的代表性。

圖2 卡群水文站洪峰流量模比系數(shù)差積曲線Fig.2 The difference product curve of flood peak flow modulus coefficient in Kaqun hydro-logical station

1.2 邊緣分布函數(shù)

P-Ⅲ型分布常用來模擬降水量、徑流量或洪峰流量等，其概率密度函數(shù)見公式(1)，P-Ⅲ分布參數(shù)可選用線性矩法進行估計，通過對樣本參數(shù)的無偏估計，便可得出P-Ⅲ曲線的參數(shù)：

(1)

(2)

(3)

式中：α、β、a0分別為P-Ⅲ分布的形狀、尺度和位置參數(shù)；Γ(·)為伽馬函數(shù)。

1.3 Copula函數(shù)與參數(shù)估計

H(x,y)=C[F(x),G(y)]

(4)

如果F和G是連續(xù)的，則C是唯一的。否則，C在RanF×RanG上(Ran表示值域)唯一確定。多種Coupla函數(shù)都可應用于建立水文變量的多維聯(lián)合分布。其中Archimedean Copula函數(shù)結構簡單，計算簡便，可以構造出形式多樣、適應性強的多變量聯(lián)合分布函數(shù)，能夠滿足大多數(shù)領域的應用要求，在實際應用中占有很重要的地位。常用的二維Archimedean Copula函數(shù)主要包括：Gumbel-Hougaard(G-H) Copula、Clayton Copula、Ali-Mikhail-Haq(AMH)Copula、Frank Copula，Nelson[9]對Archimedean Copula函數(shù)及其性質進行了詳細的介紹，見表1。

表1 函數(shù)類型與參數(shù)Tab.1 Function types and parameters

估計Copula函數(shù)參數(shù)的方法有多種，本文采用Kendall相關系數(shù)法[10]對Copula函數(shù)中的參數(shù)α進行估計。建立Kendall秩相關系數(shù)τ與α的關系，見表1，其中Kendall秩相關系數(shù)表示為：

(5)

式中：τ為秩相關系數(shù)；(xi，yi)為測點據(jù)；sgn(·)為符號函數(shù)，n為系列長度。

由實測資料計算得到Kendall秩相關系數(shù)后，可以很容易地求得Copula函數(shù)的參數(shù)α。

1.4 擬合優(yōu)度評價

在運用Archimedean Copula函數(shù)研究實際問題的過程中，為選擇合理的Copula函數(shù)來正確描述變量間相關結構，需要對各函數(shù)進行擬合度評價和檢驗。常見檢驗方法有均方根誤差RMSE法、AIC信息準則法、Genest-Rivest圖形分析法。

(1)均方根誤差RMSE法。Yue等最先采用基于理論值與實測值的擬合曲線來直觀地檢驗多維分布的擬合情況[11,12]。而后Beersma和Buishand開始將該方法用于Copula函數(shù)的擬合檢驗[13]。Zhang和Singh基于該方法，通過計算理論值和實測值的均方根誤差RMSE，定量地評估了擬合誤差大小。RMSE的計算公式為[14]：

(6)

式中：E為數(shù)學期望；N為樣本容量；pc為實測概率值；p0為Copula函數(shù)計算得到的理論概率值。

(2)AIC信息準則法。Zhang和Singh引入AIC指標(Akaike’s information criterion)，用于選擇合適的Copula函數(shù)[15]，其計算公式為：

AIC=Nln(MSE)+2k

(7)

式中：k是Copula函數(shù)參數(shù)的個數(shù)；AIC值越小證明擬合效果越好。

(3)Genest-Rivest圖形分析法。Genest-Rivest圖形分析法是Genest和Rivest提出的一種評價Copula函數(shù)擬合效果的方法，分別計算理論估計值Kc(t)和Ke(t)，然后點繪Kc(t)和Ke(t)關系圖，如果圖上的點都落在45°對角線上，那么表明Kc(t)和Ke(t)完全相等，即Copula函數(shù)擬合得很好[16]。

1.5 聯(lián)合分布重現(xiàn)期

在洪水頻率分析計算時，我們更關心的是洪峰、洪水總量等洪水變量超過某一特定值或者洪水多要素共同超過某一特定值時的頻率，即重現(xiàn)期。設FX(x)和FY(y)分別為洪水特征變量X和Y的分布函數(shù)，兩變量同時超過某一特定值的重現(xiàn)期如下式所示。

(8)

兩變量中任一變量超過某一特定值時的重現(xiàn)期如式(9)所示。

(9)

其中用Txy記T(X>x,Y>y)；用Tx/y記T(X>xorY>y)。

2 結果與討論

2.1 洪水特征變量的邊緣分布

運用矩法公式(3)初步估計洪峰、最大1日洪量服從的P-Ⅲ分布函數(shù)的參數(shù)值，可得洪峰和最大1日洪量的統(tǒng)計參數(shù)，以AIC準則為依據(jù)，通過計算機適線法得到卡群水文站洪峰流量、最大1日洪量頻率曲線。洪峰流量采用Q均值=2 052.59 m3/s，Cv=0.62，Cs/Cv=4.5；最大1日洪量采用W均值=1.228 億m3，Cv=0.32，Cs/Cv=5。繪制變量邊際分布經(jīng)驗頻率與理論頻率關系圖，如圖3所示，由圖3可以看出，點據(jù)均分布在1∶1斜線附近，洪峰流量相關方程為y=0.839 2x+0.043 4，相關系數(shù)R2為0.971 2；洪水總量相關方程為y=0.999 4x+0.016，相關系數(shù)R2為0.977 8；證明擬合程度較好。

圖3 邊緣分布擬合圖Fig.3 Marginal distribution fitting map

2.2 Copula函數(shù)參數(shù)估計及擬合優(yōu)度檢驗

本文選擇Copula函數(shù)中常見的4種二維Archimedean Copula函數(shù)進行參數(shù)估計，并從中選擇最優(yōu)函數(shù)擬合變量。采用Pearson′γ，Spearman′ρ和Kendall′τ進行洪峰和洪量間的相依性度量，經(jīng)計算γ=0.83，ρ=0.68，τ=0.57，說明洪峰與洪量之間存在較強相關性，因而可進行洪水變量間的聯(lián)合概率特性分析。

各Copula函數(shù)的參數(shù)估計值及擬合優(yōu)度檢驗指標列于表2，圖4 為幾種Copula函數(shù)的Kc～Ke關系圖。表2中的REMS和AIC值顯示，G-H Copula函數(shù)擬合最優(yōu)，圖4顯示G-H Copula函數(shù)的Kc～Ke關系圖上的點基本落于45°對角線附近，擬合效果優(yōu)于其他3種常用的Archimedean Copula函數(shù)。

表2 Archimedean Copula函數(shù)參數(shù)估計及擬合優(yōu)度檢驗指標Tab.2 Parameter estimation goodness of fit test indicator of Archimedean Copula function

圖4 Archimedean Copula函數(shù)Kc～Ke關系圖Fig.4 Archimedean Copula function Kc～Ke diagram

2.3 組合變量分布函數(shù)與單變量分布函數(shù)對比分析

(1)同場洪水重現(xiàn)期結果對比。由前述擬合度優(yōu)檢驗，可知葉爾羌河流域洪峰與洪量聯(lián)合分布擬合最優(yōu)的copula概率分布函數(shù)為Gumbel Copula函數(shù)，由Gumbel Copula聯(lián)合概率分布函數(shù)即可推算出不同頻率下的洪水特征值，利用Matlab繪圖工具繪制參數(shù)α為2.32的Gumbel Copula的聯(lián)合概率密度函數(shù)圖和聯(lián)合概率分布函數(shù)圖(見圖5)，由公式(8)和公式(9)可求得洪峰、洪量兩變量組合下的聯(lián)合重現(xiàn)期[T(x/y)]和同現(xiàn)重現(xiàn)期[T(xy)]，繪出兩變量等值線圖(見圖6)，根據(jù)兩變量聯(lián)合分布重現(xiàn)期等值線圖可得到兩變量組合下的重現(xiàn)期。

由于組合洪水變量采用Gumbel Copula函數(shù)擬合，而單一洪水變量均采用P-Ⅲ曲線擬合，同場次洪水的重現(xiàn)期結果不一樣。以1999年洪水為例，年最大洪峰流量為6 070 m3/s，一日最大洪量為2.29 億m3，由圖4和圖5可得該場次洪水各要素聯(lián)合分布概率及重現(xiàn)期如表3所示。如按單變量計算，發(fā)生洪峰流量為6 070 m3/s的大洪水重現(xiàn)期為51年，發(fā)生一日最大洪量為2.29 億m3的重現(xiàn)期為44年，平均每44年發(fā)生如此大的洪水，顯然是不符合實際的。按洪峰流量與一日最大洪量兩變量聯(lián)合分布計算，同時發(fā)生最大洪峰流量為6 070 m3/s，一日最大洪量為2.29 億m3大洪水的重現(xiàn)期為182年，即平均每182年發(fā)生一次，較為合理。

圖5 兩變量聯(lián)合概率密度函數(shù)和分布函數(shù)圖Fig.5 Probability density function and distribution function graph of two joint variable

圖6 兩維聯(lián)合分布重現(xiàn)期等值線平面圖Fig.6 contour map of joint distribution recurrence period

表3 1999年場次洪水特征值聯(lián)合分布重現(xiàn)期Tab.3 The recurrence period of flood eigenvalues joint distribution in 1999

由表3聯(lián)合分布重現(xiàn)期計算結果可以看出，采用單一變量擬合分布函數(shù)時，該場次洪水洪峰流量的重現(xiàn)期為51年，洪量的重現(xiàn)期為44年；采用組合變量分布函數(shù)時，該場次洪水洪峰和洪量的聯(lián)合重現(xiàn)期[T(x/y)](即相同大小的洪峰或洪量任一出現(xiàn))為27年，同現(xiàn)重現(xiàn)期[T(xy)](即相同大小的洪峰和洪量同時出現(xiàn))達182年。

(2)同重現(xiàn)期洪水參數(shù)對比。比較同一重現(xiàn)期下，單變量設計值與組合變量聯(lián)合分布條件下設計值的不同。以葉爾羌河流域洪峰-洪量聯(lián)合概率分布為例，統(tǒng)計分析了同一重現(xiàn)期下，單變量設計值與聯(lián)合變量分布設計值的區(qū)別。由重現(xiàn)期計算公式反推出某一重現(xiàn)期下，葉爾羌河流域洪峰與洪量單變量設計值與聯(lián)合變量分布設計值見表4。

表4 單變量設計值與峰量聯(lián)合分布下設計值Tab.4 Design value of single variable and joint distribution

由表4單變量設計值與峰量聯(lián)合分布下設計值計算結果可知，對于百年一遇洪水，單變量分布函數(shù)洪峰和洪量設計值分別為7 103.5 m3/s和2.56 億m3，而組合變量分布函數(shù)的設計值為7 478.83 m3/s和2.64 億m3，由此可知，在同一重現(xiàn)期下，兩變量聯(lián)合分布法推求的洪峰流量和洪量設計值較單變量設計值偏大。

3 結 語

(1)本文對葉爾羌河卡群水文站歷年水文資料進行代表性分析，認為卡群水文站具有較好的代表性。采用年最大法進行選樣，提取洪峰流量與最大1日洪量相關信息，以皮爾遜Ⅲ型曲線為基礎建立洪峰、洪量的單變量分布函數(shù)，繪制經(jīng)驗頻率與理論頻率關系圖，得到洪峰流量相關方程為y=0.839 2x+0.043 4，相關系數(shù)R2為0.971 2；洪水總量相關方程為y=0.999 4x+0.016，相關系數(shù)R2為0.977 8；證明擬合程度較好。

(2)采用RMSE法、AIC信息準則法、Genest-Rivest圖形分析法進行擬合優(yōu)度評價，3種評價方法均得出一致的結果，認為Gumbel Copula函數(shù)擬合度最優(yōu)，以此為基礎建立基于Gumbel Copula函數(shù)的洪峰與洪量聯(lián)合分布函數(shù)。

(3)采用MATLAB軟件，利用Gumbel Copula聯(lián)合概率分布函數(shù)推求了洪峰、洪量兩變量組合下的聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期，繪制了洪峰、洪量組合下的聯(lián)合分布圖及兩種重現(xiàn)期的等值線圖，比較分析了同一重現(xiàn)期下，單變量設計值與組合變量聯(lián)合分布條件下設計值的不同。結果顯示聯(lián)合分布時的設計值大于單變量設計值，即利用聯(lián)合分布設計值對水利工程偏于安全。

由此可見，基于Copula函數(shù)的兩變量聯(lián)合分布的洪水頻率分析方法能更好地描述洪水的特征，提供了分析變量間相互聯(lián)系的大量信息，是單變量擬合所不具備的，且Copula函數(shù)對于邊緣分布類型沒有限制，大大增加了這種方法的實用性和推廣性。本文研究結果可對未來葉爾羌河流域水利工程規(guī)劃設計和風險評估工作提供參考依據(jù)。