基于Copula函數(shù)的韓-榕-練三江徑流豐枯遭遇分析

吳宇彤，鐘平安，萬新宇，王 涵

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

0 引 言

韓江、榕江和練江屬于粵東沿海諸河，由北向南依次分布(見圖1)，各流域面積分別為30 112、4 628 km2和1 353 km2，均獨流入海。三江流域人均水資源分布不均，人均水資源量分別為2 750、1 138 m3/人和365 m3/人，水資源開發(fā)利用率分別為20.3%、22.4%和52%。三江水系連通工程是優(yōu)化韓江、榕江、練江等水資源配置的戰(zhàn)略性工程。三江來水的豐枯互補程度是決定工程效益的重要因素，豐枯遭遇分析是三江水系連通工程可行性研究的重要任務(wù)。

圖1 韓江-榕江-練江位置分布

徑流(或洪水)遭遇分析方法很多，鄭笑天等[1]利用頻率分析法分析水文序列的豐枯遭遇，得出王快和西大洋兩水庫的豐枯遭遇組合頻率；薛樹文等[2]基于FGM法建立了黃河與小清河年徑流量的聯(lián)合分布函數(shù)，得到了二者年徑流量豐枯遭遇的概率。在諸多方法中常見的是利用多變量聯(lián)合概率和條件概率表述遭遇特征[3]。Copula函數(shù)又稱連接函數(shù)，它可通過連接函數(shù)表達任意邊緣分布的聯(lián)合分布，并且可以描述變量間的復(fù)雜相關(guān)結(jié)構(gòu)。為此被廣泛應(yīng)用于暴雨、洪水、干旱等水文時間分析，以及水資源調(diào)度、徑流豐枯遭遇分析中。尤綺英等[4]采用不同參數(shù)估計法擬合了西安供水河流年徑流的Copula函數(shù)，得到Kendall系數(shù)估計法和極大似然估計法均可行的結(jié)果；石衛(wèi)等[5]利用Copula函數(shù)構(gòu)建了南水北調(diào)中線水源區(qū)與受水區(qū)降水的聯(lián)合分布，對干旱遭遇規(guī)律開展了探討；劉祖發(fā)等[6]利用Copula函數(shù)構(gòu)建廣東東江流域三大水庫入庫流量的聯(lián)合分布，得到三大水庫可豐枯互補利用的結(jié)論。

本文基于韓江、榕江、練江的降水和徑流資料，采用Copula函數(shù)分析三江來水豐枯遭遇規(guī)律，為三江水系連通工程的建設(shè)和運行提供參考。

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

1.1 基礎(chǔ)數(shù)據(jù)

本文收集了韓江流域的潮安站、榕江北河的赤坎站、榕江南河的東橋園站1968年～2010年的逐月徑流資料，榕江南河8個和練江6個氣象站的同期降水資料。

練江流域缺乏徑流資料，榕江南河與練江比鄰，下墊面條件接近，本文首先利用榕江南河的降雨與徑流資料，建立降雨徑流關(guān)系[7]，再根據(jù)練江流域的降雨資料得到練江的徑流過程系列。

1.2 方法原理

1.2.1Copula函數(shù)及其參數(shù)估計

利用Copula函數(shù)構(gòu)建多維聯(lián)合分布模型，分析多條河流間徑流的豐枯遭遇關(guān)系是目前認為較好的方法。Copula函數(shù)又稱連接函數(shù)。Copula函數(shù)基于Sklar定理[8]：設(shè)X、Y為連續(xù)的隨機變量，邊緣分布函數(shù)分別為FX和FY，F(xiàn)(x，y)為變量X和Y的聯(lián)合分布函數(shù)[9]，如果FX和FY連續(xù)，則存在唯一的函數(shù)Cθ(u，v)使得

F(x，y)=Cθ(Fθ(x)，F(xiàn)θ(y))，?x，y

(1)

式中，Cθ(u，v)為Copula函數(shù)；θ為待定參數(shù)。

Gumbel Copula函數(shù)[2]、Clayton Copula函數(shù)[9]和Frank Copula函數(shù)[10]是水文領(lǐng)域常用的3種Copula 函數(shù)，詳見表1。

表1 水文領(lǐng)域常用的3種Copula函數(shù)

Copula函數(shù)的參數(shù)估計一般采用Kendall系數(shù)估計法或極大似然法，后者包括全參數(shù)極大似然估計法、半?yún)?shù)極大似然估計法和逐步極大似然估計法[9-11]。

本文將采用半?yún)?shù)極大似然法，其原理為用邊緣分布的經(jīng)驗函數(shù)值代替理論函數(shù)值，利用τ與θ的關(guān)系計算得到θ，不同Copula函數(shù)的τ與θ參數(shù)均不相同。

τ為描述二維變量間非線性相關(guān)系數(shù)的Kendall相關(guān)系數(shù)，即

(2)

(3)

當(dāng)τ可知時，可由表1中的公式根據(jù)τ和θ的關(guān)系計算得到θ。

1.2.2邊緣分布優(yōu)選

本文采用GAMMA分布、LOGISTIC分布、對數(shù)正態(tài)分布(Log-Normal Distribution)對韓江、榕江、練江的包括年、汛期、非汛期的3個時間尺度，共9組徑流變量進行擬合，利用穩(wěn)健性較好、計算誤差較小的線性矩法[12](MME法)進行參數(shù)估計。

通過K-S方法檢驗邊緣分布、Copula函數(shù)擬合結(jié)果，并用均方根誤差法(RMSE)和AIC信息準(zhǔn)則法計算評價指標(biāo)，根據(jù)AIC值大小選取最優(yōu)的邊緣分布。計算公式如下

(4)

(5)

式中，pei和pi分別為經(jīng)驗頻率和理論頻率；N為徑流序列的樣品大??；M為時變邊緣分布和Copula模型的評估參數(shù)個數(shù)。

1.2.3豐枯遭遇分析

單一河流的徑流豐枯劃分呈現(xiàn)的是徑流的概率區(qū)間。水文中的徑流豐枯等級劃分，一般基于距平百分率或者頻率[13-14]。本文將年徑流豐枯等級劃分為豐水、平水、枯水3級，設(shè)枯水、平水的分界值流量為Xpk，平水、豐水的分界值流量為Xpf。若流量Xi

豐枯遭遇是指多條河流的徑流在不同概率區(qū)間的組合事件[15-17]。多條河流流量的豐枯遭遇本質(zhì)是聯(lián)合概率問題。

設(shè)2組數(shù)據(jù)分別為X和Y，Xpk、Xpf、Ypk、Ypf分別為X和Y枯平、平豐的分界值，其對應(yīng)的邊緣分布結(jié)果為upk、upf、vpk和vpf，則這兩組數(shù)據(jù)相應(yīng)不同豐枯組合概率的計算公式如下：

(1)X豐Y豐，P(X≥Xpf，Y≥Ypf)=1-upf-vpf-C(upf，vpf)；

(2)X豐Y枯，P(X≥Xpf，Y

(3)X平Y(jié)豐，P(Xpk≤X

(4)X平Y(jié)枯，P(Xpk≤X

(5)X平Y(jié)平，P(Xpk≤X

(6)X枯Y枯，P(X

C(u，v)指u和v的Copula函數(shù)，各豐枯組合的計算公式可類推。

2 結(jié)果分析與結(jié)論

2.1 流量邊緣分布的確定

選用GAMMA、Logistic和對數(shù)正態(tài)分布對韓江、榕江、練江的年、汛期、非汛期的共9組徑流變量進行擬合，采用MME法擬合參數(shù)，并利用K-S方法進行檢驗評價，適線結(jié)果的RMSE值和AIC值、K-S檢驗的P值和D值如表2、3所示。選取原則為AIC值越小越好，計算得到的P值均通過0.05的顯著性水平，D值的臨界值為DC=0.204。9組數(shù)據(jù)的邊緣分布函數(shù)均選取GAMMA分布函數(shù)。以練江年徑流為例，其頻率分布及年徑流關(guān)系見圖2、圖3。

表2 邊緣分布適線結(jié)果

圖2 練江年徑流頻率分布

圖3 練江實際與理論年徑流關(guān)系

不同尺度數(shù)據(jù)的GAMMA分布函數(shù)參數(shù)見表4。

表3 K-S檢驗

表4 GAMMA分布參數(shù)

2.2 Copula函數(shù)的確定

構(gòu)建Clayton Copula，F(xiàn)rank Copula和Gumbel Copula 3種二元Copula聯(lián)合分布模型，用半?yún)?shù)極大似然法估計Copula函數(shù)的參數(shù)，9組數(shù)據(jù)的參數(shù)估計結(jié)果和評價結(jié)果見表5。從韓江、榕江和練江年、汛期和非汛期3個時期的徑流量兩兩Copula聯(lián)合分布的擬合評價的結(jié)果來看，不同組合的不同函數(shù)的RMSE和AIC值均不同；選擇評價指標(biāo)較小的Copula函數(shù)類型，可以得到較高的模型的擬合精度。表5中帶下劃線的數(shù)值即對應(yīng)其所取的Copula函數(shù)。由表5可知，有7個組合均Gumbel Copula函數(shù)評價指標(biāo)占優(yōu)，只有汛期的韓江-練江和韓江-榕江組合Clayton Copula函數(shù)略優(yōu)?？紤]到水文領(lǐng)域相關(guān)研究中Gumbel Copula函數(shù)應(yīng)用最為廣泛，所以本文也采用Gumbel Copula函數(shù)。

表5 不同Copula函數(shù)參數(shù)估計及評價結(jié)果

2.3 三江兩兩豐枯遭遇分析

以年徑流頻率37.5%和62.5%將徑流劃分為豐水、平水、枯水3級，根據(jù)年、汛期、非汛期的韓江、榕江、練江流量的同一時間尺度兩兩間的Copula函數(shù)，可以計算得出相應(yīng)的二維聯(lián)合分布，得到不同豐枯遭遇組合情景的概率，即不同聯(lián)合分布在特定邊際分布下的聯(lián)合概率(見表6)。各時期不同組合徑流豐枯遭遇的聯(lián)合分布等值線圖見圖4。從圖4中可以查出水源區(qū)與受水區(qū)徑流遭遇發(fā)生各種徑流組合的頻率。例如，由圖4韓江-練江非汛期徑流遭遇等值線圖可以查出當(dāng)水源區(qū)韓江非汛期徑流為500 m3/s時，練江出現(xiàn)非汛期徑流為24 m3/s的概率為0.109 6。從而可以定量得出韓江、榕江、練江之間發(fā)生各種徑流組合遭遇的概率，其豐枯遭遇表現(xiàn)出以下特征。

圖4 三江流域徑流豐枯遭遇的聯(lián)合分布等值線

表6 三江豐枯遭遇頻率

(1)從空間看，榕江-練江的豐枯同步的概率大于韓江-榕江、韓江-練江，榕江-練江的年尺度同步概率為0.465 4，大于韓江-榕江、韓江-練江的0.444 2、0.415 8，榕江-練江的汛期尺度同步概率為0.547 4，大于韓江-榕江、韓江-練江的0.458 2、0.496 6。練江和榕江的同步性較好是因為練江和榕江主要水量來源的榕江南河地理距離較為接近，下墊面條件和氣候條件相似。而韓江與練江距離較遠，且本次研究的韓江潮安站位于韓江下游，與榕江、練江的下游較為接近，而韓江的主要水量來自上游，流域面積較大，故來水量較大時，距離較近的榕江-練江的同步概率較大。

(2)分析韓江、榕江、練江的豐枯遭遇，得到“豐-平”、“豐-枯”、“平-豐”等利于調(diào)水的組合概率，討論水量較大的水源區(qū)韓江、榕江對受水區(qū)練江補水有利的豐枯組合。從時間看，對于年時間尺度：①練江作為受水區(qū)，與韓江、榕江的豐枯互補概率較大，存在韓江、榕江作為水源對練江進行補水的可能性；②韓江作為第一水源，對榕江和練江利于調(diào)水的概率分別為0.466 1、0.470 1，榕江作為第二水源及過水通道，練江作為受水區(qū)，利于調(diào)水的概率為0.511 6。三者間豐枯異步且利于調(diào)水的組合概率較大，存在豐枯補償?shù)奶卣鳎邆淙?lián)合調(diào)水的可能性。韓江來水較豐的情況下，韓江可在滿足自身下游河道內(nèi)外需水的前提下，通過過水通道榕江為練江調(diào)水，同時可利用余水供給榕江，榕江也可在來水較豐的情況下，補償練江，韓江、榕江兩水源具有相機補水的特征。

(3)對于調(diào)水的主要時段非汛期而言，韓江-練江、榕江-練江的利于調(diào)水概率分別為0.443 3、0.464 1，概率越大，來水充沛的上游越可以補償來水偏少的下游，越有利于調(diào)水，這一情況對持續(xù)缺水的受水區(qū)練江而言是非常有利的。

3 結(jié) 語

本文利用Copula函數(shù)構(gòu)建了廣東粵東中小河流域韓江、榕江、練江三江徑流的二維聯(lián)合分布模型，分析其豐枯遭遇，主要結(jié)論如下：

(1)練江和榕江地理位置相近，下墊面條件和氣候條件相似，所以表現(xiàn)出較強的豐枯同步性。韓江潮安站徑流主要來自韓江上游，所以韓江與榕江、練江同步性稍差。

(2)三江兩兩豐枯異步的概率明顯大于豐枯同步的概率，韓江與練江的豐枯異步概率大于榕江與練江的概率，但均大于相應(yīng)的豐枯同步概率。

(3)非汛期韓江-練江、榕江-練江的利于調(diào)水概率分別為0.443 3、0.464 1，有利于跨流域補水。