丁倩倩， 劉友存， 焦克勤， 邊曉輝， 劉 燕

(1.江西理工大學 建筑與測繪工程學院, 江西 贛州 341000; 2.江西理工大學 資源與環(huán)境工程學院,江西 贛州 341000； 3.中國科學院 西北生態(tài)環(huán)境資源研究院， 甘肅 蘭州 730000)

河流的水沙變化不僅受流域內(nèi)降水、下墊面和河床演變等自然因素的影響，亦是水利工程和河道改造等人為因素的反映[1]。如何定量區(qū)分自然因素和人為因素對河流水沙變化的影響，已成為國內(nèi)外河流水沙研究的熱點。同時，水沙作為天然水體中所攜帶的固體顆粒，不僅對河道形態(tài)、洪水特征、河道生態(tài)和河口沉積地貌等有著很大的影響[2-3]，而且是水利和防洪工程、水資源合理利用和保護以及水土保持等方面必須考慮的重要因素[4-5]。自1753年布朗姆斯(Brahms)提出水沙起動流速與重量的1/6次方成正比關系以來[6]，國際上很多學者[7-8]對河流水沙變化及其影響因素以及河流輸沙變化所引起的河口地區(qū)的沖淤演變和環(huán)境變化等進行了觀測研究。1950s以來，隨著我國對河流水沙治理的日益重視，使得水沙學科體系逐步得到完善，并在非均勻懸移質(zhì)輸沙[9]、泥沙運動統(tǒng)計理論[10]、水流挾沙能力[11-12]、水庫淤積和河流模擬[13-14]等方面取得了一批重要成果。近年來，圍繞黃河小浪底[15]和長江三峽[16]工程的水沙問題以及黃河水沙調(diào)控[17]等重大工程問題的研究也取得了引人注目的進展。然而，針對江南典型流域輸沙過程的時間序列特征及成因分析，以及進一步探討河流輸沙與氣候變化相互關系的研究較少，尤其是定量地分析流域的徑流和輸沙時空變化特征，并探討水沙過程對自然因素與人類活動的響應研究更少。為此，本研究以贛江重要的支流桃江流域為研究對象，運用數(shù)理統(tǒng)計分析(雙累積曲線和Mann-Kendall法)和3S技術等方法，依據(jù)流域內(nèi)過去36 a逐日實測降水、徑流和輸沙資料，分析了桃江流域輸沙過程的時間序列特征及其與氣候變化的相互關系，進一步探討桃江河流輸沙與氣候和人類活動的關系，以期為桃江乃至贛江流域的社會、經(jīng)濟與生態(tài)環(huán)境的可持續(xù)發(fā)展，特別是流域內(nèi)土壤侵蝕治理、土地利用和功能調(diào)整等提供科學的依據(jù)。

1 研究區(qū)概況

桃江流域位于24°28′—25°53′N，114°11′—115°19′E，流域面積約8 440 km2。桃江流域主要由中低山和丘陵組成，海拔從南至北逐漸降低[17]。流域上游和東側(cè)邊緣山地海拔較高，呈弧形，最高海拔為1 383 m；中下游地勢雖有起伏，但較為和緩，平均海拔約300 m，而與貢江交匯處海拔低于200 m。桃江流域?qū)儆趤啛釒駶櫦撅L氣候區(qū)，氣候溫潤，光照充足，雨量豐沛。年均降水量為1 580.6 mm，降水年內(nèi)分配不均勻，3—8月(汛期)約占全年降水量的74%左右。降水量年際差別較大，豐水年約為枯水年的2倍。桃江的豐水期、平水期和枯水期特征很明顯，而徑流的年內(nèi)和年際變化受降水控制。

2 資料與研究方法

2.1 資料來源

選取了桃江流域控制站—居龍灘水文站36 a的逐日徑流量和輸沙量實測資料；以及流域內(nèi)27個站點的逐日降水量實測資料,運用泰森多邊形的插值方法由各站點降水量推算出整個流域的降水量，從而得到1980—2015年桃江流域的逐日降水量。

2.2 研究方法

2.2.1 水文法 水文法是用流域水沙的觀測資料分析水土保持措施減水減沙作用的一種方法[19]。它以治理前(基準期)的實測水沙資料，建立降水產(chǎn)水產(chǎn)沙量數(shù)學模型，然后將治理后的降水因子代入所建模型，計算出相當于治理前的產(chǎn)水產(chǎn)沙量，再與治理后的實測水沙量進行對比，其差值即為經(jīng)過治理后減少的水沙量。如果將治理前的實測水沙量視為天然產(chǎn)水產(chǎn)沙量，那么，根據(jù)治理后的降水因子用產(chǎn)流產(chǎn)沙模型計算出的產(chǎn)流產(chǎn)沙量，就相當于治理后降水條件下所產(chǎn)生的天然產(chǎn)流產(chǎn)沙量，兩個時段天然產(chǎn)流產(chǎn)沙量的差值即為降水變化對產(chǎn)流產(chǎn)沙的影響量。相應地，如果將模型計算的天然產(chǎn)流產(chǎn)沙量與同一時段實測的水沙量相減，即可視為人類活動對產(chǎn)流產(chǎn)沙的影響量。

2.2.2 雙累積曲線法 雙累積曲線法是檢驗兩個參數(shù)間關系一致性及其變化的常用方法。Searcy等[20]認為，如果雙積累曲線的斜率發(fā)生突變(Break)則意味著2個變量之間的比例常數(shù)發(fā)生了改變或者其對應的累積值的比可能根本就不是常數(shù)。若接受2個變量累積值之間直線斜率已發(fā)生改變，那么斜率發(fā)生突變點所對應的年份就是兩個變量累積關系出現(xiàn)突變的時間。

設有2個變量X(參考變量或基準變量)和Y(被檢驗變量)，在N年的觀測期內(nèi)，有觀測值Xi和Yi，其中i=1，2，3,…,N。首先對變量X和變量Y按年序計算各自的累積值，得到新的逐年累積序列Xi′和Yi′，其中i=1，2，3,…,N，即：

(1)

(2)

然后，在直角坐標系中繪制兩個變量所對應點累積值的關系曲線。繪制的曲線圖一般以被檢驗的變量為縱坐標(即Y軸)、參考變量或基準變量為橫坐標(即X軸)。

2.2.3 Mann-Kendall法(M-K法) M-K法[33]是用來評估水文氣候要素的時間序列趨勢的檢驗方法。

設一時間序列為：X1，X2，X3, …,Xn，構造一序列mi，mi表示Xi>Xj(1≤j≤i)的樣本累積數(shù)。定義dk：

(3)

其中：dk的均值、方差定義為：

(4)

(5)

在時間序列隨機獨立的假定下，定義統(tǒng)計量：

(6)

式中：UFk——標準正態(tài)分布，給定顯著性水平a，查正態(tài)分布表得到臨界值t，當UFk>t，表明序列存在一個顯著的增長或減少趨勢，繪制時間序列UFk的曲線L1，通過顯著性檢驗可得知其是否具有顯著性趨勢。將時間序列X按逆序排列，再重復上述的計算過程，并將計算值乘以-1，得到UBk,UBk的時間序列曲線用L2表示，當曲線L2超出置信區(qū)間，即表示存在明顯的變化趨勢。若L1和L2在置信區(qū)間內(nèi)存在交點，則此點為突變點的開始。

2.2.4 連續(xù)小波變換(CWT)法 如果函數(shù)Ψ(t)滿足以下容許條件：

(7)

則稱Ψ(t)為一個基本小波[21]，并定義如下的積分變換：

(8)

式中：f(t)——Ψ(t)為基本的積分連續(xù)小波變換；a——尺度因子，表示與頻率相關的伸縮；b——時間平移因子。

3 分析與討論

3.1 水沙變化

3.1.1 水沙年際變化 依據(jù)桃江流域1980—2015年的年均徑流量和輸沙率的實測數(shù)據(jù)，在1980—2015年的時間序列上，徑流量和輸沙率總體上呈降低趨勢，但徑流量的降低趨勢不是很明顯，只是略有下降。輸沙率的下降趨勢則較為明顯，降低幅度比較大。在輸沙變化過程中，1991和2004年達到極小值， 同期的徑流量也為極小值。從2007—2015年輸沙量呈顯著下降趨勢，而同期的徑流量并沒有明顯的減小。另外，年平均輸沙率通過了0.01顯著性水平的相關性檢驗，而年平均徑流量則未通過顯著性相關檢驗。

3.1.2 水沙年內(nèi)變化 依據(jù)桃江流域1980—2015年的年內(nèi)月平均徑流量和輸沙率的實測數(shù)據(jù)，桃江流域的徑流量在上半年不斷增加，6月達到最大值，下半年則不斷減小?？傮w上，徑流量在3—6月較大，為豐水期。輸沙率在3—6月遠遠大于其他月份，其中4月的輸沙率達到最大值為94.50 kg/s，5—6月稍有遞減，分別為92.83和87.58 kg/s。

3.1.3 水沙變化的突變分析 采用M-K法對桃江流域逐年徑流量和輸沙率做M-K突變檢驗，結(jié)果如圖1所示。由圖1可以看出，徑流量的UFk統(tǒng)計量在1997—1998年,2001—2003年為正值，其余年份均為負值，表明徑流量呈下降趨勢。但是UFk統(tǒng)計量值除了1989和1991年超過了顯著性水平為 的臨界值，其余年份的統(tǒng)計量值均未超過該顯著性水平的臨界值，故徑流量的下降趨勢并不顯著。從圖2總體來看，徑流量的UFk統(tǒng)計量值波動變化較多。在1991年之前，雖然其值上下波動過3次，但總體為下降趨勢，其變化值域沒有超過0。在1991—1998年UFk呈上升趨勢，而在1997—1998年超過了0值，表明徑流量為增加，緊接著在1998—2002年UFk呈先下降后上升。從2002年開始UFk開始呈明顯的下降，但其值并未超過顯著水平 的臨界值，說明徑流量雖有下降，但并不顯著。正逆序列的M-K統(tǒng)計量值UFk和UBk在2002年后出現(xiàn)了7個交叉點(分別為2003—2004年，2005年，2007年，2009—2010年，2010—2011年，2011—2012年和2012—2013年)。結(jié)合桃江流域水利工程情況(表2)可知，五渡港水庫于2004—2005年加固維修，其庫容量為3.33×107m3，桃江電站于2004年建成投入使用，其庫容量為3.71×107m3，而居龍灘電站于2007年建成投入運行，其庫容量為7.76×107m3，3個水電站都屬于大型水庫，其攔水攔沙作用比較明顯，故徑流量在上述幾個年份出現(xiàn)突變是合理的。采用M-K法對桃江流域逐年輸沙率進行的突變檢驗，由圖1可以看出，輸沙率的UFk統(tǒng)計量值除了1996年為正值，其余年份均為負值，表明輸沙率在整個時間序列中以下降為主。在2007年輸沙率的UFk統(tǒng)計量值超過了顯著性水平 的臨界值，其下降趨勢更為顯著。正逆序列的M-K統(tǒng)計量值UFk和UBk在2005—2006年間出現(xiàn)一個交叉點，且交叉點位于顯著性水平臨界值的范圍之內(nèi)，因此，可以認為2005年即為輸沙率的突變年份。

圖1 桃江流域徑流量和輸沙率的M-K突變檢驗

總體看來，徑流量和輸沙率的M-K統(tǒng)計量值絕大部分為負值，故徑流量和輸沙率均呈現(xiàn)減小趨勢。在1980—1982年徑流量和輸沙率的M-K統(tǒng)計量值UFk不僅為負值，而且呈現(xiàn)遞減趨勢，到了1983年徑流量和輸沙率的M-K統(tǒng)計量值UFk又上升到1980年的位置。1983—1991年徑流量的M-K統(tǒng)計量值UFk呈現(xiàn)波動減小，同期的輸沙率M-K統(tǒng)計量值UFk同樣呈現(xiàn)波動減小，但在其變化過程中波動起伏很小，波動起伏的時間節(jié)點基本與徑流量M-K統(tǒng)計量值UFk的時間節(jié)點相吻合。1991—1998年徑流量的M-K統(tǒng)計量值UFk呈遞增趨勢，同期的輸沙率M-K統(tǒng)計量值UFk同樣呈現(xiàn)上升趨勢，但其上升的速度較緩慢，且上升幅度沒有前者大。從1999年開始，徑流量和輸沙率的M-K統(tǒng)計量值UFk均呈現(xiàn)下降趨勢，前者的下降幅度相對較小，且起伏波動較多。后者則呈現(xiàn)遞減趨勢，1999—2006年雖有較小幅度的起伏，但總體下降幅度較大。從整個時間序列來看，徑流量和輸沙率的M-K統(tǒng)計量值UFk的變化趨勢基本上是同步的，只是在1999年前，前者的變化幅度大，起伏波動較明顯；而在1999年后，后者的減小幅度要比前者大很多。

3.1.4 水沙變化周期分析 圖2中細實線表示影響錐(the cone of influence， COI)，在該曲線以外的能量譜由于受到邊界效應的影響而不予考慮，粗實線表示通過95%置信水平顯著性檢驗的臨界值。桃江流域1980—2015年的逐月降水量、徑流量和輸沙量在時頻域中都存在不同尺度且通過0.05置信水平檢驗的顯著周期(圖2)。降水量在1980—1990年和1991—2015年均在時頻域低能區(qū)存在9—15個月的周期。徑流量在1980—1990年間、1991—1997年和2004—2013年在時頻域低能區(qū)分別存在6—16個月，10—16個月和10—15個月的周期。輸沙量在1981—1990年和1991—1997年在時頻域低能區(qū)存在10—16個月的周期。

3.2 驅(qū)動因素分析

3.2.1 降水 桃江流域在1980—2015年時段內(nèi)逐年降水量如圖3所示。從圖3可以看出，在1983年達到最大值為2 007.9 mm。從1983—1989年為遞減波動趨勢，在1989年達到最小值為1 277.2 mm，1990年又上升到1 704.4 mm，而到1991年再次降到整個時間序列的最小值為1 106.9 mm。1992年上升到1 912.4 mm，在1992—2002年時段內(nèi)除1999年為1 404.6 mm外，其余年份均達到1 500 mm以上。到2003和2004年降水量再次降低，且均低于1 200 mm。2005—2015年時段內(nèi)呈現(xiàn)波動變化，但在1 500 mm上下波動，所有值均大于1 200 mm。從1980—2015年的降水量變化過程來看，呈現(xiàn)微弱的遞減趨勢，其變化趨勢可用線性擬合(圖3)，但年平均降水量沒有通過顯著性相關檢驗。

圖2 桃江流域連續(xù)小波變換能量譜

桃江流域的年平均降水量為1 580.63 mm，由流域降水量的年內(nèi)分布情況可知(圖3)，3—8月為豐水期，占全年降水量的73.7%；10—1月為枯水期，占全年降水量的12.8%;2和9月為平水期，占全年降水量的13.5%。

由1980—2015年桃江流域月降水量的時間序列單累積曲線可以看出(圖4a)，桃江流域的降水量隨著時間其累積值與時間呈線性相關關系，其線性擬合方程的擬合優(yōu)度R2達到0.999 7，說明桃江流域的降水量隨著時間的變化并沒有顯著性的增加或減少。由1980—2015年桃江流域月徑流量的時間序列單累積曲線可知(圖4b)，桃江流域徑流量總體減小趨勢不明顯，但存在小幅度的下降，結(jié)合圖1，可將2003年確定為徑流量的一個分界點，并以此求出徑流量的減少量。第一階段(1980—2002年)的線性擬合的擬合優(yōu)度R2為0.998 3，第二階段(2003—2015年)的R2為0.996 9，第二階段相對第一階段的擬合曲線其斜率有所減小，說明徑流量在第二階段相對第一階段有所減少。從2004年開始，年均徑流量減少7.26×1010m3。由流域1980—2015年的逐月輸沙率的時間序列單累積曲線可知(圖4c)，輸沙率在2005年發(fā)生了突變，其輸沙量明顯存在兩個階段的變化。第一階段(1980—2004年)的線性擬合的擬合優(yōu)度R2為0.995 7，第二階段(2005—2015年)的R2為0.896 2。第二階段的線性擬合曲線斜率明顯小于第一階段，也就是說，從2005年開始每年輸沙率的累積值減小，平均每年的輸沙量減少量為8.82×107t。由桃江流域1980—2015年降水量和徑流量的時間序列雙累積曲線可知(圖4d)，雙累積曲線在2003年發(fā)生突變，第一階段(1980—2002年)擬合曲線的擬合優(yōu)度R2為0.999 6，第二階段(2003—2015年)的R2為0.999，第二階段擬合曲線的斜率要小于第一階段的。根據(jù)“水文法”計算得到(表1)，2003—2015年桃江流域人類活動造成的年均徑流量較1980—2002年的年均徑流量減少了5.77×1010m3，占年均總徑流量減少量的76.77%；其次因降水量減少引起的徑流量減少了1.49×1010m3，占年均徑流量減少量的23.23%。人類活動和降水量變化對徑流量的影響之比約為0.8∶0.2。

表1 桃江流域減水減沙量

圖3 桃江流域降水量年際和年內(nèi)變化

圖4 降水量、徑流量和輸沙率單累積和雙累積曲線

由桃江流域1980—2015年降水量和輸沙率的時間序列雙累積曲線可知(圖4e)，雙累積曲線在2005年發(fā)生突變，可將曲線分為兩個階段。第一階段(1980—2004年)的線性擬合曲線的擬合優(yōu)度R2為0.998 3，第二階段(2005—2015年)的R2為0.993 1。由表1可以看出，人類活動造成桃江流域第二階段的輸沙量比第一階段年均減少量為8.44×107t，占總體輸沙量減少量的95.7%；其次為降水量減少造成的輸沙量減少量為3.82×106t，占總體年均輸沙量減少量的4.3%。人類活動和降水量變化對輸沙量的影響之比約為0.95∶0.05。由流域1980—2015年徑流量和輸沙量的時間序列雙累積曲線曲線可知(圖4f)，整個趨勢與降水量和輸沙率的時間序列雙累積曲線圖相一致，在2005年發(fā)生突變。1980—2004年的線性擬合曲線的擬合優(yōu)度R2為0.997 3，2005—2015年的R2為0.9174，且第二階段的擬合曲線斜率小于第一階段的斜率。由徑流量和輸沙量間的線性關系可知，2005—2015年年均徑流產(chǎn)生的輸沙量減少量為6.82×107t。

由圖4可知，輸沙量的時間序列單累積值從2005年開始呈現(xiàn)階段性的躍升變化。而降水量和輸沙率、徑流量和輸沙率的雙累積圖也呈現(xiàn)同樣的變化，故這些線性曲線的擬合優(yōu)度均不是高。結(jié)合桃江流域輸沙率的實測資料可知，輸沙率在2005年6月,2006年7月,2007年6月,2008年6月,2010年6月,2011年5月,2012年3月,2013年5月,2014年5月和2015年5月均為極大值，分別為167,140,163,66.6,79.9,57.2,48.6,50.8,52.5和52.6 kg/s, 最小輸沙率為0.37 kg/s，除去這些極大值之后的月平均輸沙率為4.66 kg/s，輸沙率的極大值要比平常值高出幾十倍甚至數(shù)百倍。故而，輸沙率的累積值會出現(xiàn)階段性的躍升變化，進而導致降水量和輸沙率、徑流量和輸沙率的雙累積曲線呈現(xiàn)同樣的變化。

對輸沙率做歸一化處理，剔除極值的影響，重新做時間序列單累積圖以及與降水量和徑流量的雙累積曲線，得到這些線性曲線的擬合優(yōu)度均在0.99以上(圖5)。

圖5 桃江流域輸沙率歸一化處理的相關性分析

3.2.2 水利工程 在人類活動對河道水沙變化的影響因素中，水土保持措施和水利工程建設是影響河道徑流和輸沙變化的主要因素。本文主要討論水利工程措施對河道來水來沙量的影響。

桃江流域的水利工程至2015年共有9項，其中水電站和水庫分別為3座(居龍灘、龍頭灘和桃江電站)和6座(龍井、上逕、白蘭、五渡港、走馬垅和中村水庫)，其主要特征參數(shù)表2。而對1980—2015年間徑流量和輸沙量有影響的水利工程，主要是1980年后建成運行的龍頭灘、居龍灘和桃江電站。由前面分析得到的徑流量和輸沙量分別在2003和2005年開始發(fā)生突變。而3座電站(龍頭灘、桃江和居龍灘電站)分別于2000,2004和2007年建成。由于龍頭灘電站的庫容較小(1.38×107m3)，故對桃江流域的徑流量和輸沙量影響也較小；居龍灘和桃江電站的庫容較大，尤其是居龍灘電站的庫容達到7.76×107m3，控制面積約為7 739 km2，占桃江流域總面積的91.7%，故對桃江流域的徑流量和輸沙量的影響也較大。綜合來看，電站建成時間與徑流量和輸沙量的突變時間基本一致。

從表1—2可知，水庫單位庫容的減水量和減沙量分別為4.49×107m3/a和6.57 t/a。運用趨勢前推法可知，1980年之前建的走馬垅、白蘭、五渡港、上逕、龍井和中村水庫對桃江流域的減水減沙量分別為4.73×1010m3/a和6.91×107t/a。由此可見，水庫庫容的大小與其攔水攔沙的能力成正相關關系。

表2 桃江流域水利工程情況

4 結(jié) 論

(1) 桃江流域1980—2015年的年均輸沙率呈現(xiàn)明顯的下降趨勢，而年均徑流量下降則不明顯趨勢；年內(nèi)月均徑流量在3—6月較大為豐水期，年內(nèi)月均輸沙率在4月達到全年最大值為94.50 kg/s，5—6月稍有遞減。

(2) 降水量在1980—1990和1991—2015年的時頻域低能區(qū)存在9—15個月的周期；徑流量在1980—1990，1991—1997和2004—2013年的時頻域低能區(qū)分別存在6—16個月、10—16個月和10—15個月的周期；輸沙量在1981—1990和1991—1997年的時頻域低能區(qū)均存在10—16個月的周期。

(3) 桃江流域1980—2015年的徑流量和輸沙率的M-K突變檢驗中得知，輸沙率在2005年發(fā)生突變，并呈現(xiàn)明顯的下降趨勢；徑流量則呈現(xiàn)不明顯的下降趨勢，2003年為突變開始點。桃江流域的徑流量從2004年開始，年均徑流量減少量為7.26×1010m3。輸沙量從2005年開始，年均減少量為8.82×107t。

(4) 2003—2015年桃江流域人類活動造成的年均徑流量較1980—2002年減少了5.77×1010m3，占年均徑流量的76.8%；而因降水量減少引起的徑流量減少量為1.49×1010m3，占年均徑流量的23.2%。 水利工程造成2005—2015年輸沙量比1980—2004年年均減少量為8.44×107t，占總輸沙量減少量的95.7%；而降水量減少造成輸沙量減少量為3.82×106t，占年均輸沙量減少量的4.3%。

(5) 在桃江流域人類活動對水沙變化的影響因素中，只考慮水利工程建設的水土保持治理措施。桃江流域的水庫單位庫容的減水量和減沙量分別為4.49×106m3/a和6 569 t/a。1980年前所建的走馬垅、白蘭、五渡港、上逕、龍井和中村水庫總體對桃江流域的減水量和減沙量分別為4.73×1010m3/a和6.91×107t/a。