基于Copula函數的韓江流域干支流洪水遭遇分析*

關 帥，林穎妍，查悉妮，丁 波，陶 楠，2，劉祖發

(1.中山大學水資源與環境研究中心//華南地區水循環與水安全廣東省普通高校重點實驗室，廣東廣州510275;2.廣東省水利廳農村機電局，廣東廣州510635)

研究洪水的遭遇問題，對流域內防洪減災和洪水規律的認識有著重要的意義。根據研究區域歷年同步洪水資料進行統計仍是目前經常使用的方法，但這只能對已發生洪水進行簡單的分析計算［1－2］，無法定量估算百年一遇或千年一遇等設計洪水的遭遇風險。洪水遭遇是一個多變量的頻率組合事件，因此可以采用多變量的分析方法對其進行研究。目前，多變量水文分析方法主要有經驗頻率法、非參數法、特定邊緣分布構成的聯合分布法、多元正態分布法、將多維轉換成一維的方法等，但各法均存在各自的局限和不足［3－5］，無法準確描述水文事件的內在規律。Copula函數克服了這些方法的不足，它通過任意邊緣分布和相關性結構來構造多維聯合概率分布，具有很大的靈活性和適應性［6－8］。近年來，Copula 函數廣泛應用于暴雨［9－10］、洪水［11－12］和干旱的多特征屬性頻率分析［13－14］、雨洪遭遇［15］以及洪水遭遇［16］等問題。

韓江流域是廣東省的第二大流域，韓江下游是人口稠密經濟發達的潮汕平原。流域內水系主要由梅江、汀江和韓江三大水系組成，其中梅江是韓江的主流上游，汀江是韓江的一級支流，韓江的主要洪水就來自梅江和汀江，若兩江相遇或梅江、汀江連續洪水疊加，都會造成韓江下游發生大洪水。根據歷史文獻記載統計，韓江平均每8年就發一次大水，對流域內人民的生命和財產安全造成了極大的威脅。目前對于韓江流域洪水的研究，主要集中在單次的暴雨洪水分析［17］、暴雨洪水的移植計算［18］以及年徑流演變過程［19］等方面，而對于韓江流域內三大水系聯合分布，尤其是干流 (上游)和支流(下游)洪水量級之間的相互作用尚無人研究過。本文利用二維和三維Copula函數，構建梅江、汀江和韓江之間的聯合分布，分析韓江流域干支流洪水遭遇的規律和特征，以期為韓江流域洪水規律的認識和防洪減災提供理論參考和決策支持。

1 研究區域概況及數據來源

1.1 研究區域概況

韓江流域是廣東境內僅次于珠江的第二大流域，在廣東的經濟社會發展中占有十分重要的地位。流域范圍包括廣東、福建、江西3省共22個縣市，流域面積約3萬km2;北面接江西省贛江流域，東面接閩江黃崗河流域，南面接榕江流域，西面接東江流域，地理位置相當優越。韓江發源于廣東省紫金縣的七星崠，上游稱琴江，流至五華水寨后稱梅江，由西南向東北至大埔縣的三河壩與由福建省流來的汀江匯合后稱韓江，見圖1。在潮州市分為東、西、北溪經汕頭市的五大出?？诹魅肽虾?，干流長470 km。

圖1 研究區域和水文站位置Fig.1 Research area and hydrological station locations

韓江流域屬亞熱帶氣候，水量豐富，多年平均降雨量高達1620 mm。由于流域內以熱帶氣旋雨影響為主，暴雨大且集中，所以洪水峰高量大，且多發生在6月和8月。通過對實測資料及歷史洪水分析，韓江流域洪水具有較明顯的時空性，干支流具有各自不同特點的洪水過程，從總體上看，梅江洪水量大，洪水過程線較肥胖，汀江峰高呈尖瘦型［18］。

1.2 數據來源

本文收集了1954－2003年梅江橫山站的實測逐日平均流量資料，1959－2010年汀江溪口站的實測逐日平均流量資料，1947－2011年韓江潮安站的實測逐日平均流量資料。由于位于汀江上游的棉花灘水庫于2001年竣工，為了避免棉花灘水庫影響選取資料的一致性，同時保持三站數據資料的同步，本文選取1959－2000年逐日平均流量，采用年最大值取樣法進行獨立取樣，得到1959－2000年橫山、溪口和潮安站最大日平均流量。

2 研究方法

2.1 Copula函數

Copula［20］函數是定義域為 ［0，1］均勻分布的多維聯合分布函數，它將聯合分布分為變量的邊緣分布和變量間的相關性結構分別處理，而且不要求變量同分布，所以可將多個任意形式的邊緣分布連接起來，生成一個多變量聯合概率分布模型。Copula函數有多種類型，在水文及相關領域常見的也是本文用來構造梅江、汀江和韓江站兩兩間聯合分布的4種二維Archimedean Copula函數分別是:Gumbel-Hougaard(GH)Copula、 FrankCopula、Clayton Copula和 Ali-Mikhail-Haq(AMH)Copula。三維Archimedean Copula函數可以由二維Archimedean Copula函數通過一重嵌套得到，本文選擇三維Gumbel-Hougaard(GH)Copula來構造梅江、汀江和韓江3者之間的聯合分布，它又分為對稱型和不對稱型兩種。

2.2 邊緣分布的確定與檢驗

確定不同變量的邊緣分布函數是構建Copula函數的第一步。我國常采用Pearson-III分布作為水文變量的分布線型，但是為了計算的準確性，本文選取國內外應用較為廣泛的四種概率分布函數，分別為Pearson-III分布 (P-III)，指數分布 (EXP)，廣義極值分布 (GEV)和對數正態分布 (LOGN)，并利用比較穩健的線性距法［21］進行參數估計。采用 Kolmogorov-Smirnov(K-S)［22］方法檢驗樣本理論分布與經驗分布的擬合程度，并用概率點據相關系數法 (PPCC)［23］、均方根誤差 (RMSE)準則和AIC最小準則法［24］評價確定出與各站數據擬合效果最好的邊緣分布。利用以下公式計算邊緣分布的經驗頻率。

式中，P為X≤xi的經驗頻率;i為樣本從小到大排列后的序號;n為樣本容量

2.3 參數估計與擬合優度檢驗

對于二維Copula函數的參數，本文選擇用相關性指標法來計算，該方法是根據兩變量間的Kendall相關系數與二維Copula函數的參數θ之間的關系來計算θ;對于三維及以上Copula函數，相關性指標法已不再適用，一般用極大似然法估計參數。其中Kendall相關系數的計算公式為:

式中，{(x1，y1)，(x2，y2)，…(xn，yn)}為隨機樣本;sign為符號函數。

然后采用AIC最小準則法和離差平方和最小(OLS)準則評價確定出最優的Copula函數。

2.4 聯合重現期和同現重現期

在實際工程應用中，常引入“重現期”概念。所謂重現期是指隨機變量的取值在長時期內平均多少年出現一次，又稱多少年一遇。重現期T與累計頻率P(水文統計上一般采用超值累計頻率)有一定的對應關系，對于暴雨洪水事件，T=1/P。對于多變量來說，有聯合重現期和同現重現期兩種定義:聯合重現期是指多個變量中至少有一個超過某一特定值時，事件發生的重現期;同現重現期是指多個變量同時都超過特定值時，事件發生的重現期。

根據定義，兩變量X、Y的聯合重現期To(x，y)計算公式為:

兩變量X、Y的同現重現期Ta(x，y)計算公式為:

式中，u、v為邊緣分布函數，C(u，v)為聯合分布函數。

2.5 條件概率分布及重現期

當采用多個變量來描述水文事件時，得到多變量聯合概率分布后，給定某一變量的范圍，另一變量在某一區間發生的概率大小問題就是條件概率的問題。比如用潮安站、橫山站和溪口站的流量作為特征變量，來描述韓江流域上下游地區的洪水組成，我們想知道當干流韓江潮安站發生大于某一重現期的洪水時，梅江橫山站和汀江溪口站發生某一量級的洪水的條件概率多大;或者當梅江橫山站和汀江溪口站發生大于某一重現期的洪水時，干流韓江潮安站發生某一量級的洪水的條件概率是多大。

用兩變量X、Y來描述一事件時，當X＞x時，Y＞y的條件概率為:

相應的條件重現期為:

用三變量X、Y、Z來描述一事件時，當X＞x時，Y＞y和Z＞z的條件概率:

相應的條件重現期為:

同理可以定義當Y＞y時，X＞x和Z＞z的條件概率及重現期，當Z＞z時，Y＞y和X＞x的條件概率及重現期。

3 結果與討論

3.1 最大日平均流量邊緣分布的確定

Copula函數不限定變量的邊緣分布，所以本文選取P-III，EXP，GEV和LOGN等4種應用較為廣泛的概率分布函數對梅江、汀江和韓江的最大日平均流量進行擬合，并進行檢驗與優選，確定出擬合效果最優的邊緣分布，擬合度檢驗結果見表1。其中4種分布的K-S檢驗值均小于臨界值0.2098，即都通過了K-S檢驗。最終選擇PPCC值最大、均方根誤差 (RMSE)最小并且AIC值最小的分布作為最優的邊緣分布，在表1中用加粗表示。從表1可以看出，潮安站和橫山站服從P-III型分布，溪口站服從指數分布。

表1 韓江流域三大水文站邊緣分布擬合檢驗結果1)Table 1 The fittest results of the marginal distributions of the three stations

3.2 最大日平均流量聯合分布模型的確定

對于二維Copula函數構造的3大水文站兩兩間的聯合分布，首先根據公式 (2)估算樣本的Kendall相關系數τ，然后再根據τ計算3大水文站最大日平均流量二維Copula聯合分布的參數，見表2。從中可以看出，潮安站與橫山站的Kendall相關系數較大，表明兩者關系密切，相關性較好。然后采用AIC最小準則和離差平方和最小準則(OLS)確定擬合效果最好的Copula函數。結果發現，三大水文站兩兩間最優Copula函數均為GH Copula。

表2 韓江流域三大水文站二維及三維聯合分布參數估計結果Table 2 Parameter estimations of 2-dimensional and 3-dimensional joint distributions of the three stations

本文采用三維Gumbel-Hougaard(GH)Copula構造3大水文站的聯合分布，首先利用極大似然法進行參數估計，然后同樣采用AIC最小準則和OLS最小準則從對稱型和非對稱型兩種函數優選出擬合效果最好的Copula函數。從表2可以看出最優的為三維對稱型GH Copula。

3.3 聯合概率分布及重現期

圖2 韓江流域三大水文站聯合概率分布Fig.2 The joint distribution of the three stations

圖3 韓江流域三水文站兩兩間聯合重現期等值線圖 (a、b、c)和同現重現期等值線圖 (d、e、f)Fig.3 The joint return period(a，b，c)and current return period(d，e，f)of Hanjiang Basin

根據三大水文站最大日平均流量兩兩間的Copula函數，可以得到相應的二維聯合分布，并據此計算特定條件下的聯合概率，見圖2(省略了橫山站和溪口站的聯合分布圖)。根據公式 (3)和(4)，計算其聯合重現期和同現重現期，并繪制其等值線圖，見圖3。根據圖2和圖3，可以得到三站之間任何一個兩兩組合的聯合分布值、聯合重現期和同現重現期。如以1960年6月的大洪水為例，潮安于11日出現年最大流量，橫山和溪口站于10日出現年最大流量，三站最大日洪峰流量分別為13300，6810，3900 m3/s，其邊緣分布重現期為53.2、50.0、5.6年，該場洪水主要是由梅江地區6001號臺風帶來的暴雨引起的。潮安和橫山站的聯合分布函數值為0.9754，聯合重現期為40.7年，同現重現期為70.4 a，意味著韓江和梅江其中之一發生超過1960年最大洪水的重現期為40.7 a，兩站均發生超過1960年最大洪水的重現期為70.4 a;同理可求得潮安和溪口站聯合分布函數值為0.8201，聯合重現期為5.6 a，同現重現期為60.6 a;橫山和溪口站的聯合分布函數值為0.8162，聯合重現期為5.4 a，同現重現期為72.5 a;潮安、橫山和溪口三站的聯合分布函數值為0.8205，聯合重現期為5.6 a，同現重現期為131 a。由以上分析可知，潮安和溪口站，橫山和溪口站，以及潮安、橫山和溪口三站的聯合重現期比較小，且都和溪口站的單變量的重現期相差不多，說明這場洪水中汀江并未發生大洪水，主要洪水是由梅江引起的，這與實際情況相吻合。

3.4 三江洪水遭遇分析

韓江、梅江與汀江兩兩之間以及三江之間的遭遇組合概率可定義為兩江 (三江)同時發生大于某一重現期的洪水，即P(Q1＞qT1，Q2＞qT2)、P(Q1＞qT1，Q2＞qT2，Q3＞ qT3)，其中 Q1、Q2、Q3分別代表洪水發生的量級，qT1、qT2、qT3分別代表T年一遇的設計值洪水。表4列出了潮安、橫山與溪口重現期分別為1000、200、100、50、20、10 a情況下的設計洪水值和兩兩之間的遭遇概率。潮安與橫山、潮安與溪口、橫山與溪口遭遇千年一遇洪水的概率分別為0.0731%、0.0465%、0.0284%，遭遇10年一遇的概率分別為7.4843%、5.0595、3.4574%，低重現期的洪水遭遇比高重現期的洪水遭遇概率大;相同重現期洪水遭遇組合，潮安與橫山遭遇的概率大于潮安與溪口，潮安與溪口遭遇的概率大于橫山與溪口。

表3 韓江流域三大水文站兩兩之間洪水發生量級組合遭遇概率Table 3 Coincidence risk analysis of flood magnitudes of two stations %

表4列出了潮安、橫山和溪口重現期分別為1000、200、100、50、20、10年情況下三站之間遭遇概率。潮安、橫山和溪口同時出現1000、200、100、50、20、10年一遇洪水的遭遇概率分別為0.0243%、0.1226%、0.2479%、0.5063%、1.3453%、2.9610%，比兩站同時發生相應重現期的洪水概率要小。對比分析表3和表4可以發現，潮安、橫山和溪口同時發生1000、200、100、50、20、10年一遇洪水的遭遇概率與橫山和溪口同時發生相應重現期洪水的概率很接近，即當梅江與汀江同時發生洪水時，韓江潮安站出現大洪水幾乎是必然事件，這與劉樹鋒［18］等基于韓江流域50年的降雨洪水資料分析的結果是一致的。

表4 韓江流域各站洪水同時發生量級組合遭遇概率Table 4 Coincidence risk analysis of flood magnitudes of the three stations %

3.5 條件概率分布及重現期

當韓江潮安站發生大于某一重現期的洪水時，干流梅江橫山站和汀江溪口站發生某一量級的洪水的條件概率有多大，對研究韓江流域洪水組成規律的認識和防洪安全等都具有十分重要的價值。表5列出了韓江潮安站發生1000、200、100、50、20、10年一遇洪水時，梅江橫山站、汀江溪口站發生各種重現期洪水的條件概率。從表5可以看出，當潮安站發生某一重現期的洪水時，橫山、溪口分別發生低重現期洪水的可能性與發生高重現期洪水的可能性非常接近，低重現期洪水的可能性略大;橫山站發生各種重現期洪水的條件概率比溪口站要大，如潮安發生1000、200、100、50、20、10 a一遇洪水時，橫山站發生同頻率的洪水的條件概率分別為 73.11%、73.18%、73.26%、73.43%、73.96%、74.84%，溪口站發生同頻率的洪水的條件概率分別為 46.46%、46.62%、46.83%、47.24%、48.49%、50.60%，這說明，潮安站與橫山站洪水遭遇組合的概率比與溪口站遭遇組合的概率要高，說明梅江的洪量占潮安站的比重比汀江大，這與橫山站的集雨面積較大，且梅江洪量較大有關。同理可以得到，潮安站發生大于某一重現期的洪水時，橫山與溪口站各種重現期洪水同時發生的條件概率以及橫山與溪口站各種重現期洪水聯合發生 (至少有一個發生)的條件概率，由于篇幅限制此結果未列出。

表5 潮安站和橫山、溪口站T年一遇洪水遭遇的條件概率Table 5 Condition probabilities of T-year flood coincidence of Chaoan and Hengshan station，Chaoan and Xikou station %

4 結論

利用Copula函數，構建了韓江流域的韓江潮安站、梅江橫山站和汀江溪口站年最大日平均流量的二維和三維聯合分布，對其洪水遭遇和條件概率進行分析，主要結論如下:

1)通過構建三大水文站的聯合分布，可以獲得水文站之間不同流量條件下的概率，以及特定聯合重現期和同現重現期下，不同水文站最大日流量的可能組合，對韓江流域的洪水組成規律的認識和防洪安全等具有十分重要的理論與實踐價值。

2)三大水文站之間洪水遭遇組合分析表明:相同重現期洪水遭遇組合，潮安與橫山遭遇的概率大于潮安與溪口，潮安與溪口遭遇的概率大于橫山與溪口;當梅江與汀江同時發生洪水時，韓江潮安站出現大洪水幾乎是必然事件。

3)三大水文站之間的條件概率表明:潮安發生1000、200、100、50、20、10 a一遇的條件下，橫山站發生同頻率的洪水的可能性均超過73.11%，溪口站發生同頻率的洪水的可能性在46.46%和50.60%之間，說明梅江的洪量占潮安站的比重比汀江大。

［1］王政祥，張明波.南水北調中線水源與受水區降水豐枯遭遇分析［J］.人民長江，2008，39(17):103－105.

［2］肖天國.金沙江、岷江洪水遭遇分析［J］.人民長江，2001，32(1):30－32.

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

［4］閆寶偉，郭生練，郭靖，等.多變量水文分析計算方法的比較［J］.武漢大學學報:工學版，2009，42(1):10－15.

［5］陳子燊，劉曾美，路劍飛.廣義極值分布參數估計方法的對比分析［J］.中山大學學報:自然科學版，2010，49(6):105－109.

［6］閆寶偉，郭生練，陳璐，等.長江和清江洪水遭遇風險分析［J］.水利學報，2010，41(4):67－73.

［7］馮平，李新.基于Copula函數的非一致性洪水峰量聯合分析［J］.水利學報，2013，44(10):1137－1147.

［8］GENEST C，FAVRE A C.Everything you wanted to know about copula modeling but were afraid to ask［J］.Journal of Hydrologic Engineering，2007，12(4):347 ～368.

［9］VANDENBERGHE S，VERHOEST N E C，ONOF C，et al.A comparative copula-based bivariate frequency analysis of observed and simulated storm events:A case study on Bartlett-Lewis modeled rainfall［J］.Water Resource Research，2011，47:W07529.

［10］SUBIMAL G.Modelling bivariate rainfall distribution and generating bivariate correlated rainfall date in neighbouring meteorological subdivisions using copula［J］.Hydrological Process，2010，24:3558－3567.

［11］CHOWDHARY H，ESCOBAR L A，SINGH V P.Identification of suitable copulas for bivariate frequency analysis of flood peak and flood volume data［J］.Hydrology Reaearch，2011，42(2/3):193－215.

［12］ZHANG L，SINGH V P.Bivariate flood frequency analysis using the copula method［J］.Journal Hydrologic Engineering ASCE，2006，11(2):150－164.

［13］KAO S C ，GOVINDARAJU R S.A copula-based joint deficit index for droughts［J］.Journal of Hydrology，2010，380:121－134.

［14］WONG G，LAMBERT M F，LEONARD M，et al.Drought analysis using trivariate copulas conditional on climatic states［J］.Journal of Hydrologic Engineering，2010，15(2):129－141.

［15］劉曾美，陳子燊.區間暴雨和外江洪水位遭遇組合的風險［J］.水科學進展，2009，20(5):619－625.

［16］康玲，何小聰.南水北調中線降水豐枯遭遇風險分析［J］.水科學進展，2011，22(1):44－50.

［17］李俊偉.韓江流域“03.5”與“03.6”暴雨洪水分析［J］.廣東水利水電，2005，(1):50－52.

［18］劉樹鋒，邱靜，黃本勝，等.韓江流域可能最大暴雨洪水移置計算［J］.水文，2009，29(1):42－46.

［19］李俊偉.基于小波變換的韓江年徑流多時間尺度演變特征分析［J］.廣東水利水電，2014，(4):5－10.

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

［21］HOSKING J R M，WALLIS J R.Regional frequency analysis:An approach based on L-moments［M］.UK:Cambridge University，1997.

［22］MASSEY J R F J.The Kolmogorov-Smirnov test for goodness of fit［J］.Journal of the American statistical Association，1951，46(253):68－78.

［23］FILLIBEN J J.The probability plot correlation coefficient test for normality ［J］.Technometrics，1975，17(1):111－117.

［24］BOZDOGAN H.Model selection and Akaike's information criterion(AIC):The general theory and its analytical extensions［J］.Psychometrika，1987，52(3):345－370.