基于Copula函數(shù)的珠江三角洲上游來水豐枯遭遇分析

張趙毅，何艷虎，2*，林柱良，陳曉宏

(1.廣東工業(yè)大學(xué)環(huán)境生態(tài)工程研究院，廣東 廣州 510006；2.廣東省流域水環(huán)境治理與水生態(tài)修復(fù)重點(diǎn)實(shí)驗室，廣東 廣州 510006；3.中山大學(xué)軟件工程學(xué)院，廣東 珠海 519082；4.中山大學(xué)水資源與環(huán)境研究中心，廣東 廣州 510275)

河流徑流量的豐枯變化對水資源的配置與管理有著重要的影響[1]。一方面，受自然地理特性和人類活動等影響，流域內(nèi)徑流的豐枯變化存在著客觀的差異性和不確定性；另一方面，由于全球氣候變化影響的深入，極端氣象事件發(fā)生的強(qiáng)度和頻率增加，降水量空間異質(zhì)性增強(qiáng)，從而使得流域徑流量的時空分布發(fā)生變化[2]，徑流豐枯遭遇特性復(fù)雜多變。因此，研究流域徑流豐枯組合遭遇對變化環(huán)境下的水資源配置與管理有著重要的理論與現(xiàn)實(shí)意義。

常用的河川徑流豐枯組合分析方法有Moran法、EFM法、FEI法、貝葉斯網(wǎng)絡(luò)理論方法、集對分析法、Copula函數(shù)方法等。如陳長清等[3]基于貝葉斯網(wǎng)絡(luò)對塔里木河流域三源流和塔里木干流的年徑流量等進(jìn)行研究；李繼清等[4]采用集對分析法和常規(guī)法對長江中上游各地區(qū)的徑流進(jìn)行了豐枯分類；馬秀峰等[5]使用游程理論分析了丹江口水庫和鄭州地區(qū)的4種徑流豐枯遭遇情況。不難發(fā)現(xiàn)，這些方法均能計算出研究區(qū)域的豐枯遭遇情況，但對于水文事件中多變量的處理能力略有不足。由于地理位置與氣候條件相近，同一流域下的河流之間，其徑流過程多具有一定的相依性。因此，構(gòu)建多維聯(lián)合分布模型使對同一流域下多條河流進(jìn)行徑流豐枯組合分析成為可能。

多維Copula函數(shù)作為一種多變量連接函數(shù)，相較于聯(lián)合分布法、多元正態(tài)分布法、非參數(shù)方法等[6]多變量水文分析方法，其對分布特性不同的要素具有同等的適應(yīng)性，打破了其他多變量分析方法的局限，已經(jīng)在多維聯(lián)合分布中得到了廣泛應(yīng)用。涂新軍等[7]利用阿基米德Copula函數(shù)構(gòu)建了枯水徑流聯(lián)合分布模型，對珠三角西水東調(diào)的西江缺水風(fēng)險進(jìn)行了分析；吳海鷗等[8]基于Copula函數(shù)對鄱陽湖水系徑流豐枯遭遇進(jìn)行了多維分析，探討了鄱陽湖水系多維豐枯遭遇同步聯(lián)合概率的變化特征；張倩等[9]使用Copula函數(shù)對鄭州市南水北調(diào)工程及黃河來水進(jìn)行了補(bǔ)償特性計算分析，認(rèn)為兩水同枯的風(fēng)險最大，易造成供水短缺問題；Thilakarathne等[10]利用Meta-elliptical Copula函數(shù)對湄公河下游流域進(jìn)行了干旱分析，認(rèn)為未來湄公河下游流域干旱事件將更加嚴(yán)重。

珠江三角洲城市群是粵港澳大灣區(qū)重要組成部分，水系發(fā)達(dá)、產(chǎn)業(yè)集聚、人口稠密，枯水期供水受上游西江、北江、東江來水變化和下游咸潮上溯影響，水資源供給安全保障面臨嚴(yán)峻挑戰(zhàn)。研究三江來水的豐枯遭遇組合，有利于在枯水期進(jìn)行水資源分配和調(diào)度，保障三角洲城市群供水安全。因此，本文選用Copula函數(shù)對三江年、汛期及非汛期徑流的二維及三維豐枯組合進(jìn)行計算分析，研究三江的徑流豐枯同步性及枯水期可能的豐枯組合，并求出對應(yīng)豐枯條件下的各月徑流量，以期為珠三角城市群水資源配置提供科學(xué)依據(jù)。

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

1.1 研究區(qū)概況與基礎(chǔ)數(shù)據(jù)

珠江三角洲(下稱珠三角)位于廣東省中南部，是由復(fù)雜河網(wǎng)形成的沖擊平原，區(qū)域內(nèi)包括廣州、深圳、珠海、佛山、江門、肇慶、東莞、惠州、中山9個城市，總面積達(dá)5.5萬km2，是中國經(jīng)濟(jì)和社會開放度最高、活躍度最強(qiáng)的地區(qū)之一，在國家發(fā)展戰(zhàn)略中具有重要地位[11]。珠三角地區(qū)的主要水系為西江、北江、東江和珠三角諸河水系(圖1)，珠江流域面積達(dá)45萬km2，西江、北江、東江為珠江的三大支流[12]，珠江三角洲位于三江的河流下游。珠江三角洲水資源豐富，其總量達(dá)3 402億m3[13]，其中來自西江、北江、東江的過境水量達(dá)2 941億m3。

研究分別選取馬口、石角、博羅水文站作為西江、北江和東江的控制站。馬口站集水面積為353 120 km2，地處西、北江三角洲的頂端，是西江進(jìn)入珠三角的首個控制站；石角站集水面積為38 363 km2，占北江流域總面積的82.1%，是北江中下游總控制站[14]；博羅站集水面積為25 325 km2，占東江流域總面積的71.7%[15]，是東江下游控制水文站。所選用水文站點(diǎn)均具有代表性。

圖1 研究區(qū)域與水系

根據(jù)實(shí)際情況，以4—9月為西江、北江和東江的汛期，以10月至次年3月為非汛期。以1959—2008年的月平均流量作為基礎(chǔ)數(shù)據(jù)，并在此基礎(chǔ)上計算出三江河流的年、汛期和非汛期徑流的徑流模數(shù)、徑流模比系數(shù)參數(shù)值。徑流模數(shù)常用于對不同流域的徑流進(jìn)行比較，經(jīng)計算得出，西江、北江和東江的多年平均徑流模數(shù)分別為0.020 4、0.025 5和0.021 4 m3/(s·km2)，且三江汛期徑流模數(shù)均約為非汛期徑流模數(shù)的3倍，因此三江的豐枯特征具有一定的關(guān)聯(lián)性。

1.2 方法原理

1.2.1Copula函數(shù)

Copula函數(shù)基于Sklar[16]定理，可看作邊緣分布在[0，1]上均勻分布的隨機(jī)向量的聯(lián)合分布函數(shù)。常用的多元Copula函數(shù)有多元正態(tài)(Gassian Copula)函數(shù)[17]、多元t-Copula函數(shù)[18]和阿基米德Copula函數(shù)3類，其中阿基米德Copula函數(shù)又包括Clayton Copula[19]、Frank Copula[20]和Gumbel Copula[21]3種主要函數(shù)。5種Copula函數(shù)公式如下。

Gassian Copula函數(shù)的公式為：

C(u1,u2,…,un;R)=ΦR(Φ-1(u1),Φ-1(u2),…,Φ-1(un))

(1)

式中R——相關(guān)系數(shù)矩陣；ΦR——n元標(biāo)準(zhǔn)正態(tài)分布的分布函數(shù)；Φ-1——標(biāo)準(zhǔn)正態(tài)分布函數(shù)的逆函數(shù)。

多元t-Copula函數(shù)的公式為：

(2)

式中R——相關(guān)系數(shù)矩陣；tR,k——n元t分布的分布函數(shù)；k——自由度。

Clayton Copula函數(shù)的公式為：

C(u1,u2,…,un;θ)=(u1-θ+u2-θ+…+

un-θ-1)-1/θ

(3)

式中θ——函數(shù)擬合參數(shù)，θ∈[0，+∞)。

Frank Copula函數(shù)的公式為：

C(u1,u2,…,un;θ)=

(4)

式中θ——函數(shù)擬合參數(shù)，θ∈R。

Gumbel Copula函數(shù)的公式為：

C(u1,u2,…,un;θ)=exp{-[(-lnu1)θ+

(-lnu2)θ+…+(-lnun)θ]1/θ}

(5)

式中θ——函數(shù)擬合參數(shù)，θ∈[1，+∞)。

采用均方根誤差RMSE、赤池信息量準(zhǔn)則(Akaike Information Criterion，AIC)進(jìn)行Copula函數(shù)優(yōu)度檢驗，以確定最優(yōu)的Copula函數(shù)。檢驗方式的函數(shù)公式為：

(6)

AIC=nln(MSE)+2l

(7)

(8)

式中Pei——經(jīng)驗頻率；Pi——理論頻率；n——徑流序列的樣品大小；l——邊緣分布函數(shù)參數(shù)個數(shù)。

RMSE與AIC的數(shù)值越小，說明Copula函數(shù)計算的聯(lián)合分布概率值越靠近經(jīng)驗概率值，擬合效果越好。

1.2.2邊緣分布

為應(yīng)用Copula函數(shù)，首先用適宜的邊緣分布描述各單個變量。常用的邊緣分布函數(shù)及分布擬合曲線有正態(tài)分布 (NORM)、泊松分布(POISS)、極值分布(EV)、廣義極值分布(GEV)、伽馬分布(GAMMA)、韋布爾分布(WBL)、廣義帕累托分布(GP)等函數(shù)。采用極大似然公式進(jìn)行三江徑流分布函數(shù)的參數(shù)估計，公式為：

(9)

(10)

(11)

式中L(θ)——似然函數(shù)；F(xi;θ)——邊緣分布密度函數(shù)，θ為待估算參數(shù)。

應(yīng)用Kolmogorov-Smirnov(K-S)作各邊緣分布擬合優(yōu)度檢驗，確定邊緣分布函數(shù)。根據(jù)顯著性檢驗原則，取0.05顯著性水平，P值大于0.05，則肯定原假設(shè)，認(rèn)為徑流分布與選定擬合的分布函數(shù)相同，選定分布函數(shù)可用作徑流分布擬合。

1.2.3豐枯遭遇分析

常用的豐枯劃分方法一般基于距平百分率劃分枯、平、豐來水，或者基于徑流累積頻率(P)，如以頻率37.5%、62.5%或頻率25%、75% 作為分界線劃分枯、平、豐3種水平年。為了細(xì)分豐枯來水，本研究采用水利部頒布的《地表水資源調(diào)查和統(tǒng)計分析細(xì)則》中的相關(guān)規(guī)定，以徑流累積頻率(P)12.5%、37.5%、62.5%、87.5%為分界線將徑流劃分為特枯來水(0%

在構(gòu)建的Copula函數(shù)為二維結(jié)構(gòu)時，設(shè)徑流量為X，P頻率對應(yīng)的水量為Xp，Pk對應(yīng)豐枯分級的數(shù)值偏小端，Pf對應(yīng)豐枯分級的數(shù)值偏大端，當(dāng)二維結(jié)構(gòu)下兩河流流量分別是X和Y，邊緣分布分別是u和v。可以推導(dǎo)出二維聯(lián)合分布概率的公式為：

P(Xpk

(12)

同上，在三維Copula結(jié)構(gòu)下，3條河流流量分別為X、Y、Z，邊緣分布分別為u、v、w，可以推導(dǎo)出三維結(jié)構(gòu)下的聯(lián)合概率公式為：

P(Xpk

(13)

2 結(jié)果與分析

2.1 徑流量邊緣分布的確定

2.1.1年、汛期和非汛期徑流量邊緣分布的確定

運(yùn)用式(9)—(11)對三江的年、汛期和非汛期徑流的7種邊緣分布函數(shù)進(jìn)行擬合，并利用K-S方法進(jìn)行顯著性檢驗，得到的P值見表1，P值越大，顯著性越高，擬合效果越好。

表1 邊緣分布函數(shù)擬合顯著性

由表1可知，泊松分布和廣義帕累托分布均沒有通過顯著性檢驗，其他分布函數(shù)對變量要素分別有不同的擬合效果。對西江年徑流擬合較好的有正態(tài)分布、廣義極值分布和韋布爾分布；對西江汛期徑流擬合效果較好的有正態(tài)分布、廣義極值分布和伽馬分布；對西江非汛期徑流擬合效果較好的只有廣義極值分布。對北江年、汛期和非汛期徑流擬合效果較好的有廣義極值和伽馬分布函數(shù)。對于東江年徑流和汛期徑流擬合較好的有正態(tài)分布和廣義極值分布函數(shù)；而東江非汛期徑流擬合效果較好的有廣義極值分布和伽馬分布。分布函數(shù)擬合見圖2，根據(jù)擬合情況與P值大小選取每個徑流要素的最佳邊緣分布函數(shù)，并對最終選用的函數(shù)作密度概率函數(shù)圖(圖3)。

a)年徑流

d)年徑流

g)年徑流

a)年徑流

徑流邊緣分布函數(shù)的選用及參數(shù)見表2。表2中分別代表三江徑流年、汛期和非汛期徑流的模比系數(shù)值。

表2 西江(F1)、北江(F2)、東江(F3)徑流邊緣分布函數(shù)及參數(shù)

2.1.2月徑流量邊緣分布的確定

對三江月徑流進(jìn)行邊緣分布函數(shù)擬合，用于確定不同豐枯組合下的西江、北江和東江月徑流量。根據(jù)邊緣分布函數(shù)的顯著性檢驗結(jié)果(表1和圖2)，廣義極值分布對三江徑流要素都具有較好的擬合度，因此使用廣義極值分布對三江月徑流進(jìn)行擬合，月徑流邊緣函數(shù)擬合曲線見圖4。圖中藍(lán)色曲線為月平均模比系數(shù)的經(jīng)驗累積概率曲線，紅色曲線則是對應(yīng)的邊緣函數(shù)擬合曲線。

圖4 西江(F1)、北江(F2)和東江(F3)月徑流分布函數(shù)擬合

2.2 徑流豐枯組合分析

2.2.1二維Copula徑流聯(lián)合分布

將西江、北江和東江的年、汛期和非汛期徑流分別兩兩組合，依次采用五種Copula函數(shù)對聯(lián)合分布函數(shù)進(jìn)行參數(shù)估計和擬合檢驗(表3)。在二維Copula擬合中，F(xiàn)rank Copula函數(shù)的檢驗值比其他Copula函數(shù)都小，且對所有的徑流要素都呈現(xiàn)較好的擬合度。因此選用Frank Copula函數(shù)作為二維徑流豐枯擬合函數(shù)。通過圖5可知，三江年徑流豐枯具有較高的同步性，當(dāng)一條河流出現(xiàn)特枯(或特豐)時，其他河流出現(xiàn)特豐(或特枯)的概率幾乎為零。

依據(jù)式(12)計算出二維來水豐枯組合的概率值(圖6)。可以發(fā)現(xiàn)，同豐枯組合1、7、13、19、25的聯(lián)合概率較其他組合大，說明當(dāng)一條河流出現(xiàn)某一種豐枯類型，另一條河流出現(xiàn)相同的豐枯類型的概率總是最大的。 西江和北江、西江和東江、北江和東江的年徑流同步概率分別為0.448 2、0.340 4、0.350 5；汛期徑流同步概率分別為0.475 0、0.369 9、0.246 7；非汛期徑流同步概率分別為0.264 4、0.423 6、0.273 9。說明西江和北江年徑流及汛期徑流同步性比西江和東江、北江和東江大。而在非汛期，西江和東江、北江和東江的豐枯同步性更加強(qiáng)烈。二維豐枯同步的5種組合概率和都接近或者超過0.3，相比二維豐枯異步其他組合的概率，呈現(xiàn)了較強(qiáng)的豐枯同步性，證明了西江、北江和東江之間較強(qiáng)的豐枯同步關(guān)系。

表3 西江(F1)、北江(F2)和東江(F3)Copula函數(shù)二維聯(lián)合分布參數(shù)估計及擬合檢驗結(jié)果

a)F1-F2

a)年徑流

d)年徑流

g)年徑流

2.2.2三維Copula徑流聯(lián)合分布

選取了3種阿基米德Copula函數(shù)進(jìn)行三維來水豐枯聯(lián)合分布參數(shù)估計及擬合檢驗(表4)。

由表4可知，F(xiàn)rank Copula函數(shù)的RMSE和AIC最低，因此選用Frank Copula函數(shù)作為三維Copula函數(shù)進(jìn)行來水豐枯組合分析。采用式(13)計算125種三維來水豐枯組合的概率值，結(jié)果見圖7。在年、汛期和非汛期徑流尺度上，三江徑流同豐枯的概率分別為0.085 7、0.086 0、0.068 5，比二維組合下的同豐枯概率值小，說明三江來水同豐枯的概率較二江來水同豐枯概率降低。三維組合中，1(同特枯)、32(同偏枯)、63(同平水)、94(同偏豐)、125(同特豐)是同豐枯組合，年徑流的概率值分別為0.006 5、0.024 7、0.018 1、0.027 3、0.009 2，并且從圖7可知，同偏枯、同平水和同偏豐的組合概率較其他組合概率都偏大，但由于豐枯異步組合數(shù)量多，因此同豐枯組合總體概率較低，即三江徑流同步豐枯相對而言屬于小概率事件。且較高概率的豐枯組合大部分集中于32—94的組合區(qū)間，說明任意組合中，有2條河流出現(xiàn)偏枯或平水或偏豐是較大概率事件。在年、汛期和非汛期徑流尺度上，同特枯和同特豐都是概率值極小的組合，相對而言，出現(xiàn)“特枯-特枯-偏枯”和“特枯-偏枯-偏枯”的概率值更大，是枯水年更有可能的豐枯組合。

表4 西江(F1)、北江(F2)、東江(F3)Copula函數(shù)三維聯(lián)合分布參數(shù)估計及擬合檢驗結(jié)果

根據(jù)邊緣分布函數(shù)和三維豐枯聯(lián)合概率計算結(jié)果，得到三江同步豐枯組合下的聯(lián)合概率、年平均流量(表5)和月平均流量(圖8)。

a)F1

b)F2

c)F3

表5 西江、北江和東江年徑流同步豐枯組合

b)F2

c)F3

由表5與圖8可知，在同特枯或同偏枯來水條件下，三江年平均流量能達(dá)到同平水來水條件下的55%或77%，但三江非汛期(1—3月、10—12月)流量處于較低水平，說明三江在汛期仍能提供較為充足的來水，但在非汛期三江同枯遭遇，上游來水減少，珠三角城市群面臨水資源短缺風(fēng)險。考慮到三江具有一定的豐枯同步性，且三維豐枯組合中“特枯-特枯-偏枯”和“特枯-偏枯-偏枯”的概率值較高，因此如何在上述兩種豐枯組合下的非汛期，滿足珠三角城市群的用水需求，將是一個值得注意的研究方向。

3 結(jié)論

本文利用Copula函數(shù)構(gòu)建了珠江流域三大支流西江、北江、東江的二維與三維聯(lián)合分布模型，分析了三江年、汛期和非汛期徑流的豐枯遭遇，計算了同豐枯組合下的三江年和月徑流量。主要結(jié)論如下。

a)在二維與三維Copula函數(shù)聯(lián)合分布擬合中，F(xiàn)rank Copula函數(shù)的RMSE和AIC檢驗值最小，且對三江所有徑流要素都呈現(xiàn)較好的擬合度。

b)三江兩兩間徑流豐枯同步的概率大于豐枯異步的概率。其中在全年和汛期，西江和北江的徑流豐枯同步性較好；而在非汛期，西江和北江、北江和東江之間的徑流豐枯同步性更加強(qiáng)烈。

c)三江三維Copula函數(shù)聯(lián)合分布結(jié)果表明，隨著維數(shù)的增加，三江來水同豐枯組合(同特枯、同偏枯、同平水、同偏豐、同特豐)的概率降低，豐枯異步的概率升高。同時，在年、汛期和非汛期徑流角度，同特枯和同特豐都是概率值極小的組合，“特枯-特枯-偏枯”和“特枯-偏枯-偏枯”的概率值更大。

d)在“特枯-特枯-偏枯”和“特枯-偏枯-偏枯”來水條件下，三江在非汛期來水量較少，因此在非汛期，珠三角城市群將面臨著嚴(yán)峻的供水安全問題。