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

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

(1.廣東工業大學環境生態工程研究院，廣東 廣州 510006；2.廣東省流域水環境治理與水生態修復重點實驗室，廣東 廣州 510006；3.中山大學軟件工程學院，廣東 珠海 519082；4.中山大學水資源與環境研究中心，廣東 廣州 510275)

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

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

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

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

1 方法與數據

1.1 研究區概況與基礎數據

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

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

圖1 研究區域與水系

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

1.2 方法原理

1.2.1Copula函數

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

Gassian Copula函數的公式為：

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

(1)

式中R——相關系數矩陣；ΦR——n元標準正態分布的分布函數；Φ-1——標準正態分布函數的逆函數。

多元t-Copula函數的公式為：

(2)

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

Clayton Copula函數的公式為：

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

un-θ-1)-1/θ

(3)

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

Frank Copula函數的公式為：

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

(4)

式中θ——函數擬合參數，θ∈R。

Gumbel Copula函數的公式為：

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

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

(5)

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

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

(6)

AIC=nln(MSE)+2l

(7)

(8)

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

RMSE與AIC的數值越小，說明Copula函數計算的聯合分布概率值越靠近經驗概率值，擬合效果越好。

1.2.2邊緣分布

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

(9)

(10)

(11)

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

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

1.2.3豐枯遭遇分析

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

在構建的Copula函數為二維結構時，設徑流量為X，P頻率對應的水量為Xp，Pk對應豐枯分級的數值偏小端，Pf對應豐枯分級的數值偏大端，當二維結構下兩河流流量分別是X和Y，邊緣分布分別是u和v。可以推導出二維聯合分布概率的公式為：

P(Xpk

(12)

同上，在三維Copula結構下，3條河流流量分別為X、Y、Z，邊緣分布分別為u、v、w，可以推導出三維結構下的聯合概率公式為：

P(Xpk

(13)

2 結果與分析

2.1 徑流量邊緣分布的確定

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

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

表1 邊緣分布函數擬合顯著性

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

a)年徑流

d)年徑流

g)年徑流

a)年徑流

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

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

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

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

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

2.2 徑流豐枯組合分析

2.2.1二維Copula徑流聯合分布

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

依據式(12)計算出二維來水豐枯組合的概率值(圖6)。可以發現，同豐枯組合1、7、13、19、25的聯合概率較其他組合大，說明當一條河流出現某一種豐枯類型，另一條河流出現相同的豐枯類型的概率總是最大的。 西江和北江、西江和東江、北江和東江的年徑流同步概率分別為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。說明西江和北江年徑流及汛期徑流同步性比西江和東江、北江和東江大。而在非汛期，西江和東江、北江和東江的豐枯同步性更加強烈。二維豐枯同步的5種組合概率和都接近或者超過0.3，相比二維豐枯異步其他組合的概率，呈現了較強的豐枯同步性，證明了西江、北江和東江之間較強的豐枯同步關系。

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

a)F1-F2

a)年徑流

d)年徑流

g)年徑流

2.2.2三維Copula徑流聯合分布

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

由表4可知，Frank Copula函數的RMSE和AIC最低，因此選用Frank Copula函數作為三維Copula函數進行來水豐枯組合分析。采用式(13)計算125種三維來水豐枯組合的概率值，結果見圖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可知，同偏枯、同平水和同偏豐的組合概率較其他組合概率都偏大，但由于豐枯異步組合數量多，因此同豐枯組合總體概率較低，即三江徑流同步豐枯相對而言屬于小概率事件。且較高概率的豐枯組合大部分集中于32—94的組合區間，說明任意組合中，有2條河流出現偏枯或平水或偏豐是較大概率事件。在年、汛期和非汛期徑流尺度上，同特枯和同特豐都是概率值極小的組合，相對而言，出現“特枯-特枯-偏枯”和“特枯-偏枯-偏枯”的概率值更大，是枯水年更有可能的豐枯組合。

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

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

a)F1

b)F2

c)F3

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

b)F2

c)F3

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

3 結論

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

a)在二維與三維Copula函數聯合分布擬合中，Frank Copula函數的RMSE和AIC檢驗值最小，且對三江所有徑流要素都呈現較好的擬合度。

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

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

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