基于Copula函數的三峽庫區萬州段蓄水前后降雨量-徑流量關系分析

王曉菊，毛海濤,2，黃慶豪，程龍飛

(1.山西農業大學 城鄉建設學院，山西 太谷 030801；2.武漢大學 水利水電學院，湖北 武漢430072)

1 研究背景

三峽水庫自2003年蓄水以來，改變了庫區內水文循環的時空分布，具有非線性、多尺度等突出特性[1-2]。近年來，多位學者對三峽庫區蓄水后流域內降雨-徑流關系進行了研究。如鄔昀等[3]在三峽水庫早期蓄水期對庫區內徑流量進行了回歸分析，推求了庫區各分流域未來徑流量；Bosshard等[4]對三峽水庫蓄水初期用半分布式水文模型PREVA對庫區內水文參數進行分析，表明三峽水庫在蓄水初期的度汛性能良好；Fang等[5]利用逐步多元回歸方法分析了三峽大壩建成后徑流變化情況，指出降雨是引起徑流變化的主要原因。上述學者均針對三峽庫區蓄水對水文情勢的影響進行了研究，但在時間或地域上有一定的局限性，未將三峽水庫各個蓄水期的水文參數進行對比研究，也未針對庫區內部的水文情勢進行探索。此外，三峽庫區降雨-徑流與水庫調蓄的內在關系還存在影響因素多且分布復雜等特點，需要探索新的方法進行研究。

三峽庫區萬州段位于長江上游下段，地處三峽庫區腹地，常年受三峽水庫調蓄影響，降雨-徑流關系與天然河道具有明顯的差別，其內在機理對于庫區內部水文情勢有著重要的價值。Copula函數具有較強的捕捉隨機變量間非線性關系等優點，可構建隨機變量的聯合分布，以診斷變量間的相關性結構變化[6]。因此，本文在構建三峽庫區萬州段降雨-徑流聯合分布的基礎上，采用Copula函數診斷流域降雨-徑流相關性結構的變異情況，以期為三峽庫區降雨-徑流等多水文要素關系變異診斷提供一條新的思路。

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

2.1 研究區域概況

三峽庫區萬州段位于長江中上游地區東經107°55′22″～108°53′25″、北緯30°24′00″～31°14′58″之間，屬亞熱帶季風濕潤帶，氣候四季分明，日照充足，雨量充沛，霜雪稀少，無霜期長。研究區年平均氣溫17.7 °C，年平均年日照時數1 484.4 h，年平均降水量1 243 mm，年平均蒸發總量達110×108m3。全區東西長97.25 km，南北長67.25 km，總面積為3 457 km2。境內河流縱橫，均屬長江水系，河流和溪澗切割深、落差大、高低懸殊，呈枝狀分布。三峽水庫庫區主要水系圖如圖1所示。

圖1 三峽水庫庫區主要水系圖

2.2 水文測點概況

萬縣水文站位于重慶市萬州區牌樓水廠處，東經108°24′，北緯30°47′，萬縣水文站為長江上游下段及三峽水庫庫區重要的控制站，是國家基本站(流量、泥沙一類精度站)，也是庫區內水文、泥沙的重要控制站和代表站，距三峽大壩約288 km，集水面積974 881 km2。萬縣水文站測船采用走航式ADCP常規測驗方法測定年徑流量，采用20cmJDZ05翻斗式自記雨量計觀測降雨量。其地理位置見圖1。

2.3 數據來源

收集了1977-2017年三峽水庫蓄水前后41 a的降雨和徑流資料。其中，最大年降雨量為1 678.5 mm，最小年降雨量為808.3 mm，多年平均降雨量為1 178 mm；最大年徑流量為4 518×108m3，最小年徑流量為2 573×108m3，平均年徑流量3 705×108m3。將1977-2017年年降雨量和年徑流量按照降序排列后得到二者的關系如圖2所示。

圖2 1977-2017年三峽水庫萬州段降雨量-徑流量關系

由圖2可知，在三峽水庫蓄水前后，庫區內的降雨量-徑流量關系總體呈現出一定的線性相關趨勢，說明庫區內降雨與徑流之間仍具有水文參數的一般規律。

2.4 徑流序列平穩性檢測

考慮到水文序列自相關性對檢測結果的影響，采取TFPW算法對各時段的實測降雨量和徑流量進行序列相關性檢驗，即代入Pettitt-BS耦合模型進行編譯檢驗[7-8]。經檢驗，三峽庫區徑流序列潛在變異點分別為2003、2007和2013年，與三峽水庫施工蓄水、初期蓄水和試驗性蓄水的時間相吻合。

為了驗證變異點的檢驗效果，對不同時段的徑流序列統計參數進行對比分析，如表1所示。由表1可見，總體而言，2007年前后的徑流量均值有一定差異，標準差、方差、偏態系數變化不大，一階偏自相關系數差異較大，尤其是水庫蓄水前后差異較為明顯。所以，庫區內部水文序列特性與天然河道相比存在變異性，三峽水庫蓄水過程對二者關系有一定的影響，需要進一步評估。

表1 1977-2017年不同時段研究區徑流序列統計參數對比

3 降雨-徑流Copula函數的建立

3.1 邊緣分布函數

由于Copula函數模型不受邊緣分布型式的限制，可以擇優采用邊緣分布[9]，對年降雨量、年徑流量在1977-2002(階段1)，2003-2012(階段2)，2013-2017(階段3)3個階段的邊緣分布分別采用皮爾遜Ⅲ型分布(P-Ⅲ)[10]、對數正態分布(logarithmic normal distribution, log n)[10]和Gumbel分布[12]進行經驗頻率擬合計算，其計算式分別如下。

皮爾遜Ⅲ型分布函數：

(1)

(2)

(3)

(4)

對數正態分布函數：

(5)

(x≥0)

式中：P為X≤x的擬合理論概率；u=Ex為樣本的對數均值；δ為樣本的對數均方差。

Gumbel分布函數：

(6)

(7)

μ=K-rη

(8)

式中：P為X≤x的擬合理論概率；μ為位置參數；η為尺度參數；S為樣本的均值；K為樣本的均方差；r為歐拉常數。

采用矩估計法對上述不同邊緣分布模型進行參數估計，分別計算出降雨量和徑流量對應公式(1)～(8)中的關鍵參數，結果見表2。

由表2確定出1977-2017年不同階段各邊緣分布函數，進而繪制萬縣水文站年降雨量和年徑流量的各邊緣分布曲線。同時將各實測的數據與邊緣分布曲線進行對比，如圖3所示。

表2 1977-2017年不同階段各邊緣分布函數降雨量和徑流量對應參數的估計結果

圖3 1977-2017年不同階段萬縣水文站年降雨量和年徑流量的邊緣分布曲線

由圖3可知，不同階段年降雨量和年徑流量的邊緣分布曲線呈現出不同的趨勢，但3種邊緣分布曲線與實際數據均能較好地擬合。

將上述3種邊緣分布函數擬合的經驗頻率與實際經驗頻率進行誤差分析，通過評價指標選擇擬合效果最好的(即誤差最小)的緣分布函數進行聯合分布計算。

3.2 評價指標

分別計算出各邊緣分布對應的均方根誤差(root mean squared error,RMSE)和赤池信息準則(Akaike information criterion，AIC)指標，模型的檢驗結果見表3。由表3可知，依據指標值最小準則，降雨分布模型擬合優度評價指標在階段1的最小值分別為78.2和80.7(RSME和AIC，下同)，階段2最小值分別為34.6和70.0，階段3最小值分別為112.3和134.4；徑流分布模型擬合優度評價指標在階段1的最小值分別為193.9和202.4，階段2最小值分別為132.3和193.2，階段3最小值分別為80.5和90.3。進而遴選出不同階段降雨量和徑流量的最佳的邊緣分布函數進行聯合分布計算。

表3 1977-2017年不同階段萬縣水文站降雨量和徑流量分布模型擬合優度檢驗結果

3.3 Copula聯合分布的建立

選用Clayton Copula和Gumbel-Hougaard Copula兩種構造函數來建立聯合分布函數[13]，結合Copula函數模型對其進行計算，以得到降雨-徑流的聯合分布函數。不同階段兩種Copula聯合分布函數參數θ的計算結果如表4所示。

由表4可見，庫區蓄水前后的降雨-徑流聯合分布有所改變，θ值也需采用RMSE與AIC指標進行檢驗，即選取較小的指標值來確定出不同階段最佳的聯合分布函數[14]。具體如下：

表4 不同階段兩種Copula聯合分布函數參數θ值

1977-2002年(階段1)，θ最佳值為3.24，降雨-徑流聯合分布服從Clayton Copula分布：

φ(x,y)=(u1-3.24+u2-3.24-1)1/-3.24

(9)

2003-2012年(階段2)，θ最佳值為2.73，降雨-徑流聯合分布服從Gumbel-Hougaard Copula分布：

(10)

2013-2017年(階段3)，θ最佳值為2.02，降雨-徑流聯合分布服從Gumbel-Hougaard Copula分布：

(11)

公式(9)、(10)和(11)即為三峽庫區萬州段在各蓄水階段的Copula聯合分布，聯合分布函數模型至此已經建立完畢，后文將基于上述公式和參數進行分析。

4 三峽水庫蓄水過程對萬州段降雨量和徑流量的影響

4.1 邊緣分布的改變

在三峽水庫的蓄水過程中，庫區內部萬州段的降雨量和徑流量的邊緣分布均發生了較大的改變。

階段1：該階段長江萬州段屬于天然河道，根據表3中降雨量和徑流量分布模型擬合優度的檢驗結果，未蓄水時長江中上游的降雨量和徑流量的邊緣分布均采用皮爾遜III型分布最為合理，其CV值分別為0.19和0.14，對應的CS/CV分別為3.1和2.6(表2)，符合我國水利行業規范所推薦的河道水文參數的頻率曲線，說明此時降雨量和徑流量邊緣分布滿足一般天然河道的特點。

階段2：三峽水庫處于相對低水位蓄水運行(蓄水位為135～156 m)，庫水位變動較大。根據表3，降雨量邊緣分布采用Gumbel分布最為合理，其尺度參數η=111.5，位置參數μ=1 069.5(表2)；徑流量邊緣分布最合理的是對數正態分布，其樣本對數均值Ex=8.2，對數均方差δ=0.1(表2)。這說明在施工期蓄水和初期蓄水階段，三峽水庫萬州段已經不同于天然河道，其降雨和徑流分布開始受到影響。

階段3：隨著庫水位175 m試驗性蓄水的進行，徑流量在水庫的調蓄下趨于穩定，降雨量呈現出新的分布情況。根據表3，此時降雨量的邊緣分布采用對數正態分布最合理，其樣本對數均值Ex=7.2，樣本均方差δ=0.2(表2)；而徑流量恢復至皮爾遜III型分布，但參數有所改變，此時CV=0.07，CS/CV=3.6。

4.2 聯合分布的改變

根據3.3節降雨量-徑流量聯合分布函數的模擬計算結果，繪制1977-2017年不同階段三峽庫區萬州段降雨量-徑流量聯合分布圖及等值線圖，如圖4所示。

由圖4并結合Copula的函數特性可知，聯合分布在單一平面上的投影是該平面對應的邊緣分布函數，等值線圖中的拐點為聯合分布的精確值。圖4中用虛線連接各聯合分布概率等值線的拐點值，并稱這一虛線為聯合概率分布趨勢線。此外，對比圖4各階段聯合分布圖和分布概率等值線圖不難發現，等值線的疏密反映了降雨量和徑流量年際波動的大小，等值線越稀疏，說明該參數年際波動越大；反之，則說明參數的波動和變化較小，聯合概率分布趨勢線會偏向波動較小的參數等值線[15]。為了定量分析庫區內蓄水前后水文參數的變化，將分布概率為0.2的等值線對應的值視為小概率值，其反映了某階段較難發生的水文情況；將分布概率為0.8的等值線對應的值視為大概率值；將分布概率為0.2～0.8的等值線范圍稱為分布區間，反映了某階段水文變化程度?；诖?，分析各階段變化規律如下：

圖4 1977-2017年不同階段研究區降雨量-徑流量Copula聯合分布模擬結果

階段1，降雨量和徑流量年際變化和波動均較大，聯合概率分布趨勢線基本位于分布概率等值線圖的對角線位置，符合天然河道水文參數變化規律。研究區域年降雨量均值為1 150 mm，Copula模擬年降雨量對應的小概率值和大概率值分別為971.1和1 358.2 mm。年徑流量的小概率值和大概率值分別為3 263.7×108和4 129.0×108m3。

階段2，徑流量受三峽水庫蓄水水位波動(蓄水位在135～156 m之間波動)的影響，年徑流量年際變化較大，而期間降雨并未受到較大影響。聯合概率分布趨勢線偏向降雨量等值線，斜率相對階段1變大。年降雨量均值為1 133.8 mm，年均降雨量較階段1減少1.4%，變化幅度較小。年降雨量對應的小概率值和大概率值分別為901.3和1 139.5 mm，較階段1的年降雨量小概率值和大概率值均有減少，其減幅分別為7.2%和16.1%，其降雨量聯合分布區間變幅減小了38.4%，減幅明顯。年徑流量對應的小概率值和大概率值分別為3 453.5×108和4 140.5×108m3，較階段1年徑流量的小概率值和大概率值分別增加5.8%和0.3%，年徑流量分布區間由865.3×108m3減小至687×108m3，減小了20.6%，減幅明顯。

階段3，三峽水庫試驗性蓄水后并出現常態化(蓄水位175 m)，受水庫調蓄作用的影響，年徑流量處于相對穩定狀態，年際變化較小，而年降雨量等值線間距大，年際變化也較大，降雨量處于整體增大的狀態。此時，聯合概率分布趨勢線明顯偏向徑流量等值線，斜率變小。年降雨量均值為1 317.4 mm，較階段1年降雨量均值增大14.6%，增幅明顯，對應的小概率值和大概率值分別為1 178.8和1 660.5 mm，較階段1均有增大，增幅大約20%，年降雨量分布區間變幅增加了24.5%。年徑流量對應的小概率值和大概率值分別為3 466.8×108和3 838.5×108m3，與其他兩個階段相比，其小概率值均增加而大概率值均減小，與階段1相比，小概率值增幅為6.2%；大概率值減幅為7.0%，變幅區間371.7×108m3，減少57%。

4.3 邊緣—聯合分布趨勢分析

為了更好地分析三峽水庫蓄水前后研究區降雨-徑流關系和預測未來水文變化趨勢，將邊緣分布和聯合分布進行整合分析，所得趨勢線分布圖如圖5所示。根據圖5，擬合出不同階段年徑流量和年降雨量邊緣—聯合分布趨勢線的線性表達式，如表5所示。

圖5 1977-2017年不同階段研究區降雨量和徑流量邊緣—聯合分布趨勢

表5 1977-2017年不同階段研究區年徑流量和年降雨量邊緣—聯合分布趨勢線的線性表達式

由表5中各擬合式的R2值可知，各擬合線的誤差均較小，可以作為預測的數據支撐，將各階段具體分析如下：

階段1，年降雨量和年徑流量的線性趨勢基本一致。階段2，兩線有一定的發散趨勢，其中年徑流量的擬合線斜率比階段1減小747，減幅為40.3%，初始值則增加284，增幅為9.3%。階段3的年徑流量趨勢線斜率較階段1、2明顯減小，變化趨緩，表明在三峽水庫的調控下，年徑流量已趨于穩定；階段3的年降雨量分布也同樣變化明顯，起始值比階段1增大27.8%，但斜率僅增大7.1%，表明研究區在三峽水庫蓄水完畢后降雨總量增大，但與自然階段相近。

4.4 研究區未來降雨量-徑流量變化預測

三峽庫區萬州段的年降雨量-徑流量受水庫調蓄等因素影響，將處于一個比較穩定的分布區間，2013年水庫試驗性蓄水以后，根據年徑流量線性擬合式y=571x+3 490可推測未來徑流分布，三峽庫區萬州段最小年徑流量大致為3 490×108m3(±5%)，最大年徑流量約為4 055×108m3(±5%)。同理根據年降雨量的線性擬合公式y=792.5x+1 048，可知該區域最小年降雨量為1 048 mm(±5%)，最大值不超過1 842 mm(±5%)。

此外，根據上述關系還可以推斷出不同保證率情況下，對應的降雨量和徑流量。

5 討 論

5.1 降雨量變化的討論

階段3降雨量較階段1、2均有所增加，且維持相對穩定的波動區間。降雨是水文循環的重要環節，降雨量主要受海陸位置、地形地貌等宏觀層面因素的影響，王文鵬等[1]研究發現三峽水庫蓄水后，上游流域短歷時降水集中出現的次數增加，降水強度增大，該類降水量占總降水量的比例較大，而長歷時降水出現頻次降低[16]。

長江流域與其他地區許多河流相比，整個流域的降雨量非常大，尤其在長江中上游地區，而年際變化很小。三峽水庫蓄水至175 m后，流域庫區水體表面積擴大了3～4倍[17]，從而產生豐富的水汽，在山區特殊的地形地貌下，促使水汽上升，加大了降水概率。當然，作為水文循環的一部分，區域降雨量的主要影響因素還是宏觀層面的，上述分析只是研究區在階段3降雨量增加的可能原因。

5.2 徑流量減小的原因

三峽水庫蓄水后，萬州段水文關系出現變化還可能與以下因素有關。

(1)上游水庫大壩的修建：除了宏觀水文循環因素，導致徑流量減小的主要原因還是人為作用。自2003年三峽水庫蓄水以來，水庫上游修建了大量的梯級電站和用于灌溉、調水等的水庫[18]。根據國際大壩委員會(International Commission on Large Dams, ICOLD)數據集[19]，長江流域有22 185座水壩，總庫容量為2 320×108m3，約占長江年徑流量的26.0%；其中，位于萬州上游的大壩數量約占總數的74.3%。在所有水壩中有840座為灌溉型水壩，另有762座水壩以灌溉為主要目的，兩者總庫容為55×108m3，占長江年徑流量的0.6%。

(2)水土保持：導致徑流減少的另一個因素為上游植被覆蓋率的變化。自三峽水庫蓄水后，近十多年來國家在水土保持和植被恢復方面投入力度越來越大[20]，長江上游植被覆蓋率明顯增大。植被覆蓋的增大使沿江地區的水土保持能力更好，降雨水量更多地被滯留在了陸地上。

(3)城鎮化建設和人口發展：萬州及其上游長江流域城鎮化建設發展迅猛，人口急劇增加，從而導致對水的需求量增大。但也有研究表明，在城鎮化較高地區，這些增加的生活用水量部分又以廢水形式最終返回河流，因而對流域徑流量的影響有限[21]。因此，關于該部分的影響評估仍需要開展大量的研究工作。此外，影響長江上游徑流量的因素還有很多，如跨流域調水等。

6 結 論

(1)Copula模型能較好診斷三峽庫區特定河段的水文關系，這為庫區河道多水文要素關系變異診斷提供了新思路。

(2)三峽水庫萬州段在天然河道階段，降雨量和徑流量均采用皮爾遜III型分布最為合理，其Copula聯合分布采用θ=3.24的Clayton Copula分布較好；施工期和初步蓄水階段，降雨量分布變為Gumbel分布，徑流變為對數正態分布，其聯合分布變為θ=2.73的Gumbel-Hougaard Copula分布；試驗性蓄水階段，降雨量分布變為對數正態分布，徑流量恢復至皮爾遜III型分布，其聯合分布變為θ=2.02的Gumbel-Hougaard Copula分布。

(3)施工期及初步蓄水階段較天然河道階段的降雨量聯合分布區間變幅減小38.4%；年徑流量變幅區間減少20.6%。試驗性蓄水階段的年降雨量較天然河道增大，變幅區間增加24.5%，年徑流量變幅區間減少57%，變化梯度趨緩。

(4)預測三峽庫區萬州段未來年徑流量應不小于3 490×108m3(±5%)，最大年徑流量不超過4 055×108m3(±5%)；年降雨量不小于1 048 mm(±5%)，最大年降雨量不超過1 842mm(±5%)。