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

(1.山西農(nóng)業(yè)大學(xué) 城鄉(xiāng)建設(shè)學(xué)院，山西 太谷 030801；2.武漢大學(xué) 水利水電學(xué)院，湖北 武漢430072)

1 研究背景

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

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

2 研究區(qū)域概況及數(shù)據(jù)來源

2.1 研究區(qū)域概況

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

圖1 三峽水庫庫區(qū)主要水系圖

2.2 水文測點概況

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

2.3 數(shù)據(jù)來源

收集了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年年降雨量和年徑流量按照降序排列后得到二者的關(guān)系如圖2所示。

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

由圖2可知，在三峽水庫蓄水前后，庫區(qū)內(nèi)的降雨量-徑流量關(guān)系總體呈現(xiàn)出一定的線性相關(guān)趨勢，說明庫區(qū)內(nèi)降雨與徑流之間仍具有水文參數(shù)的一般規(guī)律。

2.4 徑流序列平穩(wěn)性檢測

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

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

表1 1977-2017年不同時段研究區(qū)徑流序列統(tǒng)計參數(shù)對比

3 降雨-徑流Copula函數(shù)的建立

3.1 邊緣分布函數(shù)

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

皮爾遜Ⅲ型分布函數(shù)：

(1)

(2)

(3)

(4)

對數(shù)正態(tài)分布函數(shù)：

(5)

(x≥0)

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

Gumbel分布函數(shù)：

(6)

(7)

μ=K-rη

(8)

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

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

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

表2 1977-2017年不同階段各邊緣分布函數(shù)降雨量和徑流量對應(yīng)參數(shù)的估計結(jié)果

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

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

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

3.2 評價指標

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

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

3.3 Copula聯(lián)合分布的建立

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

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

表4 不同階段兩種Copula聯(lián)合分布函數(shù)參數(shù)θ值

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

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

(9)

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

(10)

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

(11)

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

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

4.1 邊緣分布的改變

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

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

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

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

4.2 聯(lián)合分布的改變

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

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

圖4 1977-2017年不同階段研究區(qū)降雨量-徑流量Copula聯(lián)合分布模擬結(jié)果

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

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

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

4.3 邊緣—聯(lián)合分布趨勢分析

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

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

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

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

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

4.4 研究區(qū)未來降雨量-徑流量變化預(yù)測

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

此外，根據(jù)上述關(guān)系還可以推斷出不同保證率情況下，對應(yīng)的降雨量和徑流量。

5 討 論

5.1 降雨量變化的討論

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

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

5.2 徑流量減小的原因

三峽水庫蓄水后，萬州段水文關(guān)系出現(xiàn)變化還可能與以下因素有關(guān)。

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

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

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

6 結(jié) 論

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

(2)三峽水庫萬州段在天然河道階段，降雨量和徑流量均采用皮爾遜III型分布最為合理，其Copula聯(lián)合分布采用θ=3.24的Clayton Copula分布較好；施工期和初步蓄水階段，降雨量分布變?yōu)镚umbel分布，徑流變?yōu)閷?shù)正態(tài)分布，其聯(lián)合分布變?yōu)棣?2.73的Gumbel-Hougaard Copula分布；試驗性蓄水階段，降雨量分布變?yōu)閷?shù)正態(tài)分布，徑流量恢復(fù)至皮爾遜III型分布，其聯(lián)合分布變?yōu)棣?2.02的Gumbel-Hougaard Copula分布。

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

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