基于Copula函數的葉爾羌河流域洪水要素聯合分布研究

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

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

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

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

1 數據與方法

1.1 研究區域與數據來源

1.1.1 研究區域

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

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

1.1.2 數據來源及資料代表性分析

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

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

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

1.2 邊緣分布函數

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

(1)

(2)

(3)

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

1.3 Copula函數與參數估計

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

(4)

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

表1 函數類型與參數Tab.1 Function types and parameters

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

(5)

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

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

1.4 擬合優度評價

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

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

(6)

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

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

AIC=Nln(MSE)+2k

(7)

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

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

1.5 聯合分布重現期

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

(8)

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

(9)

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

2 結果與討論

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

運用矩法公式(3)初步估計洪峰、最大1日洪量服從的P-Ⅲ分布函數的參數值，可得洪峰和最大1日洪量的統計參數，以AIC準則為依據，通過計算機適線法得到卡群水文站洪峰流量、最大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。繪制變量邊際分布經驗頻率與理論頻率關系圖，如圖3所示，由圖3可以看出，點據均分布在1∶1斜線附近，洪峰流量相關方程為y=0.839 2x+0.043 4，相關系數R2為0.971 2；洪水總量相關方程為y=0.999 4x+0.016，相關系數R2為0.977 8；證明擬合程度較好。

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

2.2 Copula函數參數估計及擬合優度檢驗

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

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

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

圖4 Archimedean Copula函數Kc～Ke關系圖Fig.4 Archimedean Copula function Kc～Ke diagram

2.3 組合變量分布函數與單變量分布函數對比分析

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

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

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

圖6 兩維聯合分布重現期等值線平面圖Fig.6 contour map of joint distribution recurrence period

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

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

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

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

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

3 結 語

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

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

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

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