萬家全
(桂林理工大學 環境科學與工程學院,廣西 桂林 541000)
珠江三角洲地區是中國經濟最發達地區之一,西江是珠江流域的干流,是珠江流域最大的水系,其水沙變化直接影響著珠江的水沙變化[1]。開展西江流域干流的水文特征研究對珠三角地區的經濟發展、航道水運、生態環境治理及相關科研具有指導性意義。
前人對西江流域水沙變化的研究多基于單個干流水文控制站或者干流的中下游水文站的資料,運用常規方法分析得到控制站附近河流航道水沙變化[2],本文基于西江流域干流上游紅河水的遷江站、中游潯江的大湟江口站以及中下游梧州站和高要站的同步水文資料和采用模擬退火算法進行非線性回歸分析,最后創新性地將4個站虛擬成“西江站”總站,沿程上全面分析其徑流量和輸沙量的變化規律,進而初步探尋影響西江流域干流水沙產生突變的因素。
西江流域干流遷江站、大湟江口站和高要站的徑流泥沙數據來自中華人民共和國水利部刊發的2002—2016年《中國河流泥沙公報》,梧州水文站輸沙量和徑流量數據來自于梧州市水文水資源局。數據為1957—2016年的年時間序列,時間序列無缺失且各站數據時間一致。
采用Mann-Kendall趨勢檢驗法[3-4]、五年滑動平均法和累積距平法[5-7]定量與定性分析4個代表站年徑流量和年輸沙量的演變趨勢;結合T檢驗法[8-9]、累積距平法和有序聚類分析方法[10],定性和定量分析年徑流量與年輸沙量的突變點;采用模擬退火法進行非線性擬合,獲得年輸沙量和年徑流量間的定性關系,進而分析其年徑流和年輸沙量的變化規律。
統計分析西江流域干流上中下游主要控制站1957—2016年平均徑流量實測資料,見表1。1957—2016年4個水文站多年平均徑流量的偏差系數較低,均在0.177~0.216之間。其徑流極值比介于2.628~2.760,表明西江干流主要水文站年徑流量年際變化不大。

表1 4個水文站多年平均徑流量和輸沙量及其偏差系數
參照SL 250—2000《水文情報預報規范》,將年徑流的距平值百分率Pi劃分為5個級別:枯水年Pi< -20%,偏枯年Pi在-20%~-10%,平水年Pi在-10%~10%,偏豐年Pi在10%~20%,豐水年Pi> 20%。4個水文站年徑流量距平值百分率見圖1。

圖1 4個水文站年徑流量距平值百分率
由圖1可知,干流各水文站年徑流量豐、枯年波動規律基本一致。1955—1960年4個站平水年和枯水年次數較多,可視為平水期和枯水期;同理,1960—1980年為平水期和豐水期;20世紀80年代末—90年代初為枯水期;1990—2000年為豐水期和平水期;2000—2010年回到平水期和枯水期。2015年后除了大湟江口站,其他3站出現平水年和豐水年的趨勢。可概括出西江干流豐枯年份出現規律可分3個階段:平枯期→平豐期→枯水期。遷江站的豐枯變化最為明顯,2013年Pi≤-40%的高概率,1968和1979年Pi≥50%,變化趨勢最為顯著,因此遷江站是4個站中極易發生干旱和洪澇災害的。
從站點年徑流量序列過程線可直接觀察到西江流域干流徑流量年際變化不明顯,大部分圍繞于趨勢線上下,較多年份的徑流量均高于多年平均值;其中遷江站多年徑流量總體上有下降趨勢,且年徑流量趨勢線斜率最大,減少趨勢較其他水文控制站更明顯。
根據趨勢線判斷,大湟江口站、梧州站和高要站年徑流量變化不明顯,結合Mann-Kendall檢驗結果(表2),大湟江口的標準統計量值Z為負值,|Z|< |Za/2|=1.96,表明該站年徑流量增加趨勢不顯著;同理分析出梧州站和高要站年徑流量在近59年均呈下降趨勢,但下降趨勢并不顯著。

表2 1955—2016年徑流變化趨勢的Mann-Kendall檢驗結果
注:Z為標準統計量值;α為顯著性參數。
4個水文站的年徑流量序列過程線見圖2。根據5年滑動平均分析,遷江站年徑流量在1965年前呈下降幅度稍緩趨勢,1965—1987年為平緩上升階段,1987—2000年遷江站的徑流量先平緩下降后平緩上升,2005年后急劇下降,遷江站在研究時間系列上總體呈下降趨勢。大湟江口站、梧州站以及高要站的年徑流量在20世紀70年代之前呈緩慢下降變化;第2階段為1968—1985年,此階段年徑流量的變幅較小,呈上升趨勢;第3階段大湟江口站、梧州站和高要站的下降幅度均稍大;第4階段為1995—2004年,3個水文站流量的變化趨勢均為較快上升;第5階段為2004年之后,年徑流量呈加速下降趨勢。結合一元趨勢線和Mann-Kendall趨勢檢驗法,可知大湟江口站、梧州站以及高要站的年徑流量在研究時間系列上的總體變化趨勢下降且不顯著。




圖2 年徑流量序列過程線
滑動T檢驗法(表3)和有序聚類法分析(圖3)得到2002年為年徑流量的顯著突變點,而累積距平法計算的結果則是2003年,再結合圖2確定遷江站年徑流量變化的顯著突變點為2002年,且突變前均值為677.83億m3,突變后均值543.53億m3,突變量為134.3億m3。由圖4也可以得到5年滑動平均值判斷結果,均呈緩慢下降→平緩上升→下降→加速上升→急劇下降。

表3 4個水文站年徑流量滑動T檢驗結果
注:t0為徑流量滑動T檢驗初始值。




圖3 年徑流量有序聚類曲線




圖4 年徑流量累積距平曲線
由圖3可知,1963年出現大湟江口站年徑流量最小值,同時在圖中存在2003年和1964年2個極值,利用滑動T檢驗時大湟江口站年徑流量的統計量值小于置信度的臨界值(α=0.05),表明其沒有顯著的突變點,變化趨勢緩慢(表3)。同理,梧州站與高要站均無突變點。
西江流域干流輸沙量年際變化較大,各站的輸沙量數據與控制站的面積有一定的關系,輸沙量的空間分布不均。
干流年輸沙量的極值比為15.757~551.877,大湟江口站、梧州站和高要站的極值比分別為15.757、26.310、32.510,遷江站的極值比551.877,變幅最大;4個水文站的最大值出現時段基本相同(圖5),遷江站、大湟江口站、高要站出現于1983年,梧州站出現于1982年,其主要原因是20世紀80年代西江流域干流不合理的土地利用、亂砍濫伐以及造成水土流失的工程建設;而年輸沙量的最小值都在2011—2013年間出現,主要原因是各站上游興修水利工程與實行水土保持方案以及大規模采砂現象,直接導致各站輸沙量發生了巨大變化,總體上為豐水多沙、少水少沙。與年徑流量變化趨勢相比,輸沙量過程曲線均在20世紀90年代后下降幅度變大,輸沙量年際變化顯著。




圖5 年輸沙量過程線
由圖5可知,遷江站年輸沙量在1964年前呈下降趨勢,1964—1992年平緩上升,1992年后急劇下降;大湟江口站、梧州站和高要站的年輸沙量均呈平緩上升趨勢后急劇下降。總體上,西江流域干流遷江站、大湟江口站、梧州站和高要站的年輸沙量總體上均呈下降趨勢。
西江流域干流年輸沙量變化趨勢的Mann-Kendall檢驗結果見表4。4個站年輸沙量統計值|Z|均大于置信度的臨界值(α=0.05),4個水文站年輸沙量變化趨勢均為顯著;由于4個站的統計量值Z均為負數,因此年輸沙量呈下降趨勢。

表4 西江流域干流年輸沙量變化趨勢的Mann-Kendall檢驗結果
3.3.1遷江站年輸沙量的突變點
遷江站年輸沙量有序聚類曲線和累積距平曲線見圖6。遷江站有序聚類值在1988—1997年處于較低水平,在1991年時有序聚類值最低;利用累積距平曲線分析遷江站年輸沙量,得1989—1995年累積距平值較高,1991年出現極大值點;兩種方法結果皆說明遷江站年輸沙量的跳躍時間點為1991年。通過滑動T檢驗法計算出的輸沙量跳躍時間點也是1991年(表5),且突變前輸沙量平均值為5 035.59萬t,突變后輸沙量平均值為1 014.59萬t,該點的突變量為4 021.01萬t。由此可知遷江站年輸沙量序列在1991年時發生顯著突變。


圖6 遷江站輸沙量有序聚類曲線和累積距平曲線

表5 西江流域干流年輸沙量突變性的滑動T檢驗結果
注:T為年輸沙量滑動T檢驗初始值。
3.3.2大湟江口站年輸沙量的突變點
大湟江口站年輸沙量有序聚類曲線和累積距平曲線見圖7。有序聚類分析得到的最小值出現時間點為1997年,累積距平值在較寬的時間域處于高水平階段,出現最高值5.539 78億t時間點為1998年,且最高值與1997和1999年對應的累計值相差較小;綜合分析年輸沙量的突變時間點為1997年。利用滑動T檢驗法再次確定1997年為突變點,且跳躍前后輸沙量的均值分別為6 164.63萬和1 960.63萬t,突變量為4 204萬t。


圖7 大湟江口站輸沙量有序聚類曲線和累積距平曲線
3.3.3梧州站年輸沙量的突變點
梧州站年輸沙量有序聚類曲線和累積距平曲線見圖8。梧州站在1996年時出現最低值,最低值與左右聚類值之差較小;距平值結果顯示1988和1994年的累積極大值相近,1994年值為最大;滑動T檢驗的結果與累積距平值方法結果一致,年輸沙量的突變點均為1994年。結合實際資料以及有序聚類分析結果中1996年的最小值與兩邊(1988年和1994年)之差小,確認1994年為梧州站年輸沙量序列的顯著突變點,突變前后輸沙量多年平均值為7 010.72萬和2 598.04萬t,突變量為4 412.67萬t。


圖8 梧州站輸沙量有序聚類曲線和累積距平曲線
3.3.4高要站年輸沙量的突變點
高要站年輸沙量有序聚類曲線和累積距平曲線見圖9。年輸沙量序列的統計量值|T|大于置信度的臨界值(α=0.05),由此可知高要站年輸沙量序列存在突變點,時間為1998年,跳躍前后輸沙量均值分別為7 249.33萬和2 532.17萬t,突變量為4 717.17萬t;有序聚類分析結果與滑動T檢驗結果相同,高要站年輸沙量有序聚類曲線在1998年出現最小值;累積距平值曲線的最高點所對應的時間為1999年,其值為5.943 63億t,1998年的年輸沙量累積距平值為5.900 55億t,而2000年的值為5.747 91億t,與1999前后累積距平值之差分別為0.043 08億與0.195 72億t。綜合分析3種方法,確認高要站年輸沙量的顯著跳躍點為1999年。


圖9 高要站輸沙量有序聚類曲線和累積距平曲線
假設存在一個水文站能夠控制整個西江流域干流,將此虛擬站命名西江站。通過模擬退火算法進行年徑流量非線性擬合,年徑流量變化擬合結果見圖10a),擬合得到的曲線相關系數R2為0.500 83。
由圖10a)可知,西江流域干流年徑流量呈波動變化,變化趨勢較平緩,趨勢變化具體可細分為:1957—1974年緩慢上升趨勢,1974—1986年緩慢下降趨勢,1986—1998年呈上升趨勢,1998—2011年以較快速度下降,而后至2016年呈急劇上升趨勢。
西江流域干流西江站年輸沙量擬合結果見圖10b)。經過多次調節參數,最終求得年輸沙量模擬曲線的相關系數R2為0.87321,擬合度較高。西江流域干流西江站年輸沙量變化趨勢不復雜,總體上呈下降趨勢:1970年以前呈波動變化,1970年后呈顯著下降。1960—1970年年輸沙量急劇增大。


注:年份1~59代表1957—2016年。
非線性回歸曲線如圖11所示,其擬合相關系數R2為0.507 99,說明1957—2015年西江總站年徑流量與輸沙量呈顯著的相關關系。

圖11 西江站年徑流量與年輸沙量非線性擬合曲線
模擬退火法計算的西江站年徑流量與年輸沙量非線性擬合所得的方程為:
y=p1+p2x0.5+p3x+p4x1.5+p5x2+p6x2.5+p7x3
(1)
式中:y為年輸沙量;p1~p7均為參數;x為年徑流量。
多次調整參數,最終確定出圖11擬合曲線參數值,利用式(1)計算出1957—2015年輸沙量的回歸值。統計實際年輸沙量和計算值的偏差系數Cv,其結果分別為0.529和0.269,說明實際輸沙量的變化幅度比回歸值的大。
西江站年輸沙量與計算年輸沙量過程線見圖12。年輸沙量、計算年輸沙量均與年徑流呈波動狀態,變化過程基本同步;然而實際年輸沙量變化量均比年徑流量的變化幅度大,1990年之后實際輸沙量下降速度顯著,另一方面回歸推算的年輸沙量與徑流量的變化幅度相對較小,而且是豐水多沙、少水少沙的規律。

圖12 西江站年輸沙量與計算年輸沙量過程線
有研究表明在氣候條件和下墊面綜合作用下流域徑流和侵蝕產沙[11]。氣候變化可能多方面影響徑流過程和侵蝕產沙過程,其中,最直接的影響因素是降水變化[12]。人類活動如水利工程建設和水土保持措施等均改變流域下墊面,也影響流域產匯流和侵蝕產沙過程[13]。若流域的輸沙量變化只與氣候變化相關,則徑流量與輸沙量的變化幅度不大。而由圖11可知,西江流域干流實際年輸沙量變化幅度遠大于年徑流量的變幅,其非參數回歸模擬后計算得到的輸沙量變化與年徑流變化幅度相對較小,因此,該時間段內人類活動作用應當是西江流域干流西江站輸沙量減少的主要因素。
1992年庫容為34.3億m3的巖灘水庫開始發揮攔水蓄沙作用,龍灘水電站也于2006年建成并投入使用,1990—2010年以來水利工程建成的總庫容為493.03億m3,占1960—2010年總水利工程建造庫容的94%,可見這些大型水利工程對西江流域干流輸沙量下降的影響力不容忽視[14]。此外1991年國家通過《中華人民共和國水土保持法》,運用法律規范實施水土保持措施,廣西壯族自治區人民政府出臺各種控制水土流失的政策,1979—1998年期間,廣西壯族自治區總造林面積達3 591 km2;2003—2009年,廣西壯族自治區內每年采砂量為120萬t。
根據上述資料,西江流域干流輸沙量自1990年后持續下降原因主要是受人類活動作用。
1)年輸沙量與年徑流量的關系為多水多沙、少水少沙。4個站豐、枯變化規律為:平枯期→平豐期→枯水期;遷江站極易發生洪水與干旱災害。
2)4個站年輸沙量變化均為顯著下降趨勢,均存在顯著突變點,遷江、大湟江口、梧州和高要站的年輸沙量突變點為1991、1997、1994和1999年。輸沙量最大值在20世紀80年代初出現,最小值在2011—2013年間出現。
3)遷江站年徑流量變化呈顯著下降趨勢;其他3個站的年徑流量呈平緩下降趨勢;遷江站年徑流量顯著突變點為2002年,突變前后的多年徑流均值之差為134.3億m3,其他站無顯著突變點。
4)徑流量、輸沙量的非線性擬合曲線的相關系數R2分別為0.500 83、0.873 21。年徑流量呈微波動,年際變動幅度較小;年輸沙量總體上呈下降趨勢,先急劇上升后持續下降,突變點為1990年。模擬退火法計算徑流量與輸沙量的擬合相關系數R2=0.507 99。人類活動作用導致年輸沙量發生巨大變化。