基于Copula函數的洪峰流量與洪水歷時聯合分布研究

范嘉煒，黃錦林

(1.天津大學建筑工程學院，天津 300072；2. 廣東省水利水電科學研究院，廣州 510635)

0 引 言

洪水是一種常見的自然現象，洪水頻率分析是防御洪水的關鍵技術內容。我國水文工作者在借鑒國外有關經驗的同時，結合國內水文資料和實際情況進行了大量的水文頻率分析工作。在實踐中，洪峰流量作為一個重要的特征量，常用來表征水文事件的整體過程，通過點繪頻率曲線、確定線形、估計參數等推求設計值，進而估算設計洪水。洪水歷時作為另一個洪水特征量，同樣在一定程度上反映了洪水的演進過程，并與洪峰流量具有較強的相關性，洪水歷時較短時，洪峰傳遞過程削減較快，洪峰衰減，反之，洪水歷時較長時，洪峰流量傳遞過程延長，峰值增大。要想全面地了解其具有的統計規律，就必須通過多個方面的特征屬性來對它進行定義和描述。但是，在實踐中人們往往只進行單個變量的頻率分析，即用一個特征量來代表洪水事件的整體特征: 我國水文界通過大量分析研究，得出P-Ⅲ分布函數對洪峰流量擬合效果較好[1]，作為我國較為普遍的線形；對于洪水歷時的邊際分布擬合的主要應用有Gamma分布、Gumbel分布[2]、指數分布[3]以及對數正態分布[4]等。隨著研究問題的復雜性，單變量分析將難以達到設計的要求[5]。

近年來，用Copula函數[6,7]進行洪水過程分析成為水文計算領域的一個研究熱點，這種分析方法能更好地描述洪水隨機變量間的內在規律和關系。目前Copula在水文頻率分析中主要應用于以下方面：暴雨、洪水、干旱等極值水文事件的多變量聯合分布研究[8,9]，分析防洪系統中設計洪水地區組成[10]，神經網絡徑流預報的預報因子選擇方法優化[11]，不同季節旱澇組合概率特征分析[12]，考慮氣候變化的區域性干旱預測[13]等問題，而對于洪峰流量和洪水歷時的聯合頻率分布的研究以及變量間的組合風險分析相對較少。本文以潖江流域的大廟峽為例，選取洪峰流量與洪水歷時2變量，通過數據分析后選擇P-Ⅲ分布和對數正態分布分別構建洪峰和歷時的邊緣分布函數，利用G-H Copula建立2者的聯合與同現概率分布，并探討了2變量在不同重現期情況下的組合概率分布情況，以期為防洪減災和水利工程規劃提供參考依據。

1 研究方法

1.1 邊緣分布函數

(1)P-Ⅲ型分布函數。對于水文變量的頻率曲線選擇問題，我國學者們通過大量研究發現P -Ⅲ分布函數對于洪水要素的擬合效果較好，目前應用也較為普遍。P -Ⅲ分布參數可選用概率權重法進行估計，其中H、R是與Cs有關的參數，通過對樣本的前3階距M0,M1,M2的估計，便可得出P -Ⅲ曲線的參數：

(2)

(2)對數正態分布函數。對數正態分布是一種連續型分布。它可用于描述某些呈偏態分布的資料。如果隨機變量經對數變換后服從正態分布，就說明此變量服從對數正態分布。對數正態分布擬合洪水歷時等變量具有較好的效果：

(3)

式中：μ、σ分別為變量對數值lnx系列的均值和標準差。

對數正態分布參數可選用極大似然估計法進行估計，首先建立極大似然函數，根據定義，解極大似然方程得到的參數值則為極大似然估計值，由于篇幅所限，此處只給出推導結果如下：

lL(μ,σ|x1,x2,…,xn)=

(4)

1.2 Copula函數與參數估計

Copula是定義域為[0,1]均勻分布的多維聯合分布函數，它可以將多個隨機變量的邊緣分布連接起來得到它們的聯合分布，早在1959年即被提出，但直到20世紀90年代該方法才得以迅速發展，成為統計學中的一個新的課題。

Sklar[14]定理:令H為一個n維分布函數，其邊緣分布為F1,F2,…,Fn，則存在一個n-Copula函數C，使得對任意x∈R有：

H(x1,x2,…,xn)=C[F1(x1),F2(x2),…,Fn(xn)]

(5)

Copula函數的優點在于不必要求具有相同的邊緣分布，任意邊緣分布經過Copula函數連接都可構造成聯合分布，由于變量的所有信息都包含在邊緣分布里，在轉換過程中不會產生信息失真。

Copula函數總體上可以劃分為3類:橢圓型、二次型和Archimedean型，以第3類Archimedean型函數應用最為廣泛。多維Archimedean Copula函數的構造通常是基于二維的，其結構簡單，可以構造出多種形式多樣、適應性強的多變量聯合分布函數，具有廣泛的應用價值。常見的二維Archimedean Copula函數包括： Clayton Copula、Gumbel-Hougaard(G-H)Copula 、Ali-Mikhail-Haq (AMH)Copula、Frank Copula，Nelson[15]對Archimedean Copula函數及其性質進行了詳細的介紹，見表1。

表1 函數類型與參數Tab.1 Function type and parameter

本文采用相關性指標法[16]對Copula函數中的參數θ進行估計。建立kendall秩相關系數τ與θ的關系， 其中kendall秩相關系數[17]表示為：

(6)

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

估計各種Copula函數的參數θ后，可以建立洪峰和歷時的Gumbel-Hougaard(G-H) Copula函數、Clayton Copula函數、Ali-Mikhail-Haq (AMH)Copula函數和Frank Copula函數的聯合頻率分布。

1.3 擬合優度評價

要知道哪一種Copula函數反映洪水事件的真實特性，能更準確地描述變量間的相關結構，就需要對Copula函數進行擬合優度評價，以選擇最合適的Copula函數來描述變量之間的相關結構，最終確定最優的Copula函數。常用的擬合優度評價的方法有離差平方和準則法(OLS)[18]、AIC信息準則法[19,20]等。本文采用AIC和OLS進行Copula函數的擬合優度評價以確定最優的Copula函數。

(1)AIC準則。即：

(7)

AIC=nln (MSE)+2k

(2)OLS準則。即：

(8)

式中：Femp(xi1,xi2,…,xim)為經驗頻率值；C(ui1,ui2,…,uim為理論頻率值；k為參數個數；m為維數。

Copula函數擬合越好，則AIC值與OLS值越小。

2 實例分析

2.1 背景資料

大廟峽位于廣東佛岡縣石角鎮南面的北江河段，屬山區性河流，全長約6 km，兩岸高山對峙，山勢險要，上游段河道比降0.38%～0.98%，洪水匯流快，洪峰尖瘦，10 a一遇洪水洪峰模數6～10 m3/(s·km2)，下游為丘陵平原區，主河道洪峰模數較少，約3～5 m3/(s·km2)。由于境內河流都屬山區型，集雨區山地陡峭，河床比降大。該區域全年降雨天數多達166 d，次數多，來勢猛，是廣東多暴雨的地區之一。汛期暴雨后，河道渲泄不暢，易造成洪澇災害。該地區山洪特點是暴漲暴落、歷時短，一般山洪以單峰為主，單峰型山洪歷時一般為1～3 d。山洪特性可歸納為“四大、兩快、一短”，即山洪流速大、沖刷力大、含沙量大及破壞力大；“兩快一短”是指山洪漲得快、落的快、歷時短。1960年1月1日設立省級水文站——大廟峽水文站。大廟峽水文站設立以來，在其轄下的10多個雨量站的協調配合下，至今已積累了50 a以上的降雨、水位、流量等水文測驗資料，代表性和穩定性較好。

本文采用潖江流域大廟峽水文站1960-2014年的洪水資料為例，選取該水文站的年最大洪水所對應的洪峰流量和洪水歷時作為所研究的特征變量，運用Copula函數中的G-H Copula對潖江河流域大廟峽水文站址處洪水的洪峰和歷時進行聯合分布研究，推求2變量聯合分布函數及聯合重現期，并與單變量洪水頻率計算進行比較分析，分析結果對潖江河流域的防洪減災、調度管理、洪水資源合理規劃等有重要的意義。

2.2 邊緣分布函數確定

運用概率權重法公式(2)估計洪峰服從的P-Ⅲ分布函數的參數值，運用極大似然法公式(4)估計歷時服從的對數正態分布函數的參數值，可得洪峰和歷時的統計參數，可得大廟峽水文站洪峰流量P-Ⅲ頻率曲線。采用Q均值=785 m3/s，Cv=0.65，Cs=1.22，倍比Cs/Cv=3.5。頻率曲線適線考慮了整體系列經驗頻率，不單以配合高水點據為原則。洪水歷時服從對數正態分布，即給定的變量取對數后服從正態分布。

采用文獻[21]中的極大似然法估計對數正態分布函數的參數，得到洪水歷時分布的參數為μ=3.2,σ=0.29。圖1、圖2為擬合的2變量頻率累積曲線。

圖1 洪峰流量頻率累積曲線Fig.1 The peak flow frequency curve

圖2 洪水歷時頻率累積曲線Fig.2 Flood duration frequency curve

2.3 最優Copula函數選擇

在本例中選擇Copula函數中常見的4種二維Archimedean Copula函數進行參數估計，并從中選擇最優函數擬合變量。首先，運用公式(6)計算洪峰與歷時的Kendall秩相關系數得到τ=0.713 6，運用相關性指標法計算公式計算洪峰和歷時2變量聯合頻率分布模型的參數，得到4種二維Archimedean Copula函數的參數θ，見表2。

表2 Archimedean Copula參數θTab.2 Archimedean Copula parameter θ

采用AIC信息準則法和OLS離差平方和準則法進行Copula函數的擬合優度評價以確定最優的Copula函數，計算結果見表3。

表3 AIC與OLS擬合優度評價Tab.3 Evaluation of AIC and OLS

根據表3數據顯示，Gumbel-Hougaard(G-H)Copula對于洪峰歷時的聯合分布擬合效果最好，因此，將計算得到的參數代入，建立潖江河大廟峽流域洪峰和歷時2變量G-H Copula函數的聯合頻率分布函數。

2.4 洪峰歷時聯合分析

運用前文分析計算得出的最優Copula函數——G-H Copula，建立洪峰和歷時兩特征變量的聯合分布、同現分布以及2變量的組合頻率分布[22,23]。同現頻率用公式(9)表示，組合頻率分布用公式(10)表示，聯合重現期和同現重現期分別用公式(11)、(12)表示：

H′(x1,y1)=P(X1>x1,Y1>y1)=

1-u1-u2+C(u1,u2)

(9)

(10)

(11)

(12)

式中，X1和Y1為假定的具有相關關系的特征變量序列；u1,u2分別為邊緣分布函數；C(u1,u2)為聯合頻率分布函數。

由G-H Copula函數擬合后得到2變量的聯合頻率分布，根據公式(9)可以得到2變量同現頻率分布，見圖3、圖4。根據其等值線圖(圖5、圖6)可分別查到給定2變量條件下洪水變量發生的頻率。以2變量聯合頻率分布為例，洪峰小于705 m3/s情況下，歷時在30 h以下的頻率為0.5；洪峰小于1 095 m3/s情況下，歷時在30 h以下的頻率為0.7。同理，以2變量同現頻率分布為例，洪峰大于1 155 m3/s且洪水歷時在25 h以上的頻率為0.2；洪峰大于465 m3/s且洪水歷時在25 h以上的頻率為0.5。

圖3 洪峰、歷時聯合頻率分布圖Fig.3 Distribution of joint frequency

圖4 洪峰、歷時同現頻率分布圖 Fig.4 Distribution of simultaneous frequency

圖5 洪峰、歷時聯合頻率分布等值線圖Fig.5 Contour of joint frequency

圖6 洪峰、歷時同現頻率分布等值線圖Fig.6 Contour of simultaneous frequency

根據公式(11)、(12)可以得到2變量聯合重現期和同現重現期分布，見圖7、圖8。同時可以得到2變量在不同重現期下的各種組合，與單變量的結果相比能更全面地反映洪水要素之間的相關關系。基于前面得到的洪峰流量與洪水歷時的邊緣分布，計算得到不同設計重現期下的單變量設計值。利用聯合重現期等值線圖(圖9)和同現重現期等值線圖(圖10)可以分別得出在此設計值下的重現期。計算結果見表4。例如，給定邊緣分布設計重現期為200 a，由邊緣分布函數求得單變量情況下的洪峰設計值為2 850 m3/s，歷時設計值為52.3 h。由圖9查得此單變量設計值下的聯合重現期為187.5 a，由圖10查得同現重現期為308.9 a。又如設計重現期為50 a，求得單變量情況下的洪峰設計值為2 150 m3/s，歷時設計值為45.0 h。此單變量設計值下聯合重現期為40.9 a，同現重現期為63.7 a。由此可以看出，2變量聯合重現期均低于單變量設計重現期，而同現重現期均高于單變量設計重現期。在實際的工程應用中，聯合重現期和同現重現期分別適用于2種遭遇情況，若洪峰

圖7 洪峰、歷時聯合重現期分布圖Fig.7 Distribution of joint return period

圖8 洪峰、歷時同現重現期分布圖Fig.8 Distribution of simultaneous return period

圖9 洪峰、歷時聯合重現期等值線圖Fig.9 Contour of joint return period

和歷時的設計值2者中有一個被超過，則被認為是破壞，此時應該用聯合重現期來描述實際重現期；若洪峰和歷時的設計值都被超過時才被認為是破壞，則應采用同現重現期來描述實際重現期。

圖10 洪峰、歷時同現重現期等值線圖Fig.10 Contour of simultaneous return period

邊緣分布設計重現期T0/a單變量設計值洪峰/(m3·s-1)歷時/h聯合重現期T1/a同現重現期T2/a1000354060．6851．51224．8500312057．1409．1639．1200285052．3187．5308．9100244549．086．8135．350215045．040．963．720177240．116．425．310148036．28．212．45116631．94．36．2

采用洪峰流量與洪水歷時同頻率的假定，在給定一定重現期下，基于G-H Copula函數反求得2變量的同頻率值，再通過邊緣分布函數分別得到相應的設計洪峰和設計歷時值，結果見表5。對比表4和表5可以得出，基于2變量聯合分布的設計值均高于單變量情況下計算得到的設計值。以1 000 a設計重現期為例，單變量情況下洪峰設計值為3 540 m3/s，歷時設計值為60.6 h；2變量聯合分布情況下洪峰設計值為3 675 m3/s，歷時設計值為62 h。可以看出，基于單變量推算的設計值實際上達不到所要求的標準，基于2變量聯合分布的設計結果與單變量設計結果相比更為安全、可靠。

表5 2變量聯合分布設計值Tab.5 Design value of joint distribution

根據式(10)定義，建立洪水歷時與洪峰流量的組合頻率分布(見圖11、圖12)，即當歷時小于某一設定值T1時，對應洪峰流量可能超過設定值Q1的概率。若以(T1，Q1)作為工程設計值，可得出雖然歷時未超過設定值，但仍有受災的風險。

圖11 洪峰、歷時組合風險分布圖Fig.11 Distribution of combined risk probability

圖12 洪峰、歷時組合風險等值線圖Fig.12 Contour of combined risk probability

由表6可以查出小于任一重現期的洪水歷時下，遭遇不同洪峰流量的概率。可以看出，同頻率下的洪峰與歷時的組合風險率都遠遠低于單變量情況下的頻率值，如兩變量設計重現期為200 a和100 a一遇時，組合風險率分別為0.11%和0.21%，表明2變量在同頻率設計情況下發生組合受災風險的概率很低。同一重現期下的洪水歷時與不同重現期下的洪峰流量組合頻率隨重現期的減小而增大，如洪水歷時為200 a一遇時，洪峰重現期為100 a和10 a的組合風險率分別為0.5%和9.73%；而在洪峰流量重現期一定時，較短的洪水歷時則對應較小的組合風險率，例如，當洪峰重現期為5 a一遇時，不同重現期的洪水歷時所對應的組合風險率分別為19.41%、18.96%、18.20%、15.60%、12.59%、3.6%。由于從15.6%至3.6%有顯著降低，可見5 a一遇洪峰對應的歷時在31.7～40.2 h的可能性較大，且歷時小于31.7 h時遭遇5 a一遇以上的洪峰概率較低。同理對于10 a一遇的洪峰，由于歷時在45.0～52.3 h的組合風險率為8.37%～9.73%，變化幅度很小，可知歷時在45.0 h以上的概率較小。另外，當歷時重現期低于200 a一遇時，流域遭遇200 a一遇以上洪峰的概率極低。根據以上的分析結果，可以為洪峰歷時遭遇組合概率的選取提供較為合理的科學依據。

表6 洪峰、歷時組合風險率計算成果Tab.6 Rusults of combined risk probability

3 結論和展望

選取潖江河大廟峽水文站歷年水文資料，從資料中提取洪峰流量與洪水歷時的相關信息，分別以皮爾遜Ⅲ型和對數正態分布曲線建立了洪峰、歷時的單變量分布函數，根據數據資料檢驗得到理論頻率與經驗頻率基本相符。之后以AIC、OLS等進行擬合優度評價，建立了基于G-H Copula的洪峰與歷時聯合分布。通過比較在不同洪峰、歷時組合情況下的單變量與2變量分布函數發現，后者可以更加全面地考慮洪水要素間的相關性，能更全面地反映洪水事件的整個過程，計算結果更為可靠，采用此方法進行工程規劃時也更加安全，并且Copula函數對于邊緣分布函數的無限制也大大增強了這種方法的實用性和拓展空間，對未來的水利工程設計有一定的參考價值，也為今后的多變量洪水頻率計算研究提供了新的思路。

□

[1] 馮 平,崔廣濤,胡明罡.暴雨洪水共同作用下的多變量防洪計算問題[J].水利學報, 2000,31(2):0 049-0 054.

[2] 謝 華，康英軍，丁夏平.多特征變量洪水頻率分析[J].灌溉排水學報,2012,31(6):11-14.

[3] 張 娜，郭生練，方 彬，等.暴雨事件中兩變量聯合分布研究[J].人民長江，2008,39(13):33-38.

[4] 侯蕓蕓. 基于Copula函數的多變量洪水頻率計算研究[D]. 陜西楊凌:西北農林科技大學,2010:23-26.

[5] 謝 華，黃介生.兩變量水文頻率分布模型研究述評[J].水科學進展，2008,(3):443-452.

[6] D J Dupuis. Using Copulas in hydrology: benefits cautions and issues[J].Journal of Hydrologic Engineering, 2007,(4):381-393.

[7] G Salvadori C D Michele. On the use of Copulas in hydrology: theory and practice[J].Journal of Hdrologic Engineering, 2007,(4):369-380.

[8] 郭生練，閆寶偉，肖 義，等.Copula函數在多變量水文分析計算中的應用及研究進展[J].水文,2008,28(3):1-7.

[9] 劉和昌，梁忠民，姚 軼，等.基于Copula函數的水文變量條件組合分析[J].水力發電,2014,40(5):13-16.

[10] 閆寶偉，郭生練，郭 靖，等.基于Copula函數的設計洪水地區組成研究[J].水力發電學報，2010,29(6):60-65.

[11] 陳 璐，葉 磊，盧韋偉，等.基于Copula熵的神經網絡徑流預報模型預報因子選擇[J].水力發電學報，2014,33(6):25-29.

[12] 楊志勇，袁 喆，方宏陽，等.基于Copula函數的灤河流域旱澇組合事件概率特征分析[J].水利學報，2013,44(5):556-561.

[13] Ganguli P, Reddy M J. Ensemble prediction of regional droughts using climate inputs and the SVM-copula approach [J]. Hydrological Processes, 2014,28(19):4 989-5 009.

[14] Sklar M. Fonctions de répartitionn dimensions et leurs marges[M] . Université Paris,1959- 08:229-231.

[15] Nelson R B. An introduction to Copulas[M]. New York:Springer, 1999．

[16] 杜 江，陳希鎮，于 波.Achimedean Copula函數的參數估計[J].科學技術與工程，2009,(3):637-640.

[17] Christian Genest Anne一Catherine Favre. Everything you always wanted to know about Copula modeling but were afraid to ask[J].Journal of Hydrologic Engineering, 2007,(4):347-368.

[18] 宋松柏，康 艷，荊 萍.水文頻率曲線參數優化估計研究[J].西北農林科技大學學報(自然科學版)，2008,36(4):193-198.

[19] 單國莉，陳東峰.一種確定最優Copula的方法及應用[J].山東大學學報(理學版)，2005,(4): 66-69.

[20] Cheng Wang, Ni-Bin Chang, Gour-Tsyhyeh. Copula-based flood frequency analysis at the confluences of confluences of river system[J].Hydrological Processes, 2009,(2):7 273-7 288.

[21] 張冬冬，魯 帆，嚴登華，等.云南省干旱的時空演變規律及季節連旱的概率特征分析[J].應用基礎與工程科學學報，2014,22(3):705-717.

[22] 史黎翔，宋松柏.基于Copula函數的兩變量洪水重現期與設計值計算研究[J].水力發電學報，2015,(10):27-33.

[23] 李天元，郭生練，劉章君，等.基于峰量聯合分布推求設計洪水[J].水利學報，2014,32(3):269-276.