萬家全

(桂林理工大學(xué) 環(huán)境科學(xué)與工程學(xué)院，廣西 桂林 541000)

珠江三角洲地區(qū)是中國經(jīng)濟(jì)最發(fā)達(dá)地區(qū)之一，西江是珠江流域的干流，是珠江流域最大的水系，其水沙變化直接影響著珠江的水沙變化[1]。開展西江流域干流的水文特征研究對珠三角地區(qū)的經(jīng)濟(jì)發(fā)展、航道水運(yùn)、生態(tài)環(huán)境治理及相關(guān)科研具有指導(dǎo)性意義。

前人對西江流域水沙變化的研究多基于單個干流水文控制站或者干流的中下游水文站的資料，運(yùn)用常規(guī)方法分析得到控制站附近河流航道水沙變化[2]，本文基于西江流域干流上游紅河水的遷江站、中游潯江的大湟江口站以及中下游梧州站和高要站的同步水文資料和采用模擬退火算法進(jìn)行非線性回歸分析，最后創(chuàng)新性地將4個站虛擬成“西江站”總站，沿程上全面分析其徑流量和輸沙量的變化規(guī)律，進(jìn)而初步探尋影響西江流域干流水沙產(chǎn)生突變的因素。

1 數(shù)據(jù)來源及分析方法

西江流域干流遷江站、大湟江口站和高要站的徑流泥沙數(shù)據(jù)來自中華人民共和國水利部刊發(fā)的2002—2016年《中國河流泥沙公報(bào)》，梧州水文站輸沙量和徑流量數(shù)據(jù)來自于梧州市水文水資源局。數(shù)據(jù)為1957—2016年的年時間序列，時間序列無缺失且各站數(shù)據(jù)時間一致。

采用Mann-Kendall趨勢檢驗(yàn)法[3-4]、五年滑動平均法和累積距平法[5-7]定量與定性分析4個代表站年徑流量和年輸沙量的演變趨勢；結(jié)合T檢驗(yàn)法[8-9]、累積距平法和有序聚類分析方法[10]，定性和定量分析年徑流量與年輸沙量的突變點(diǎn)；采用模擬退火法進(jìn)行非線性擬合，獲得年輸沙量和年徑流量間的定性關(guān)系，進(jìn)而分析其年徑流和年輸沙量的變化規(guī)律。

2 西江流域干流徑流年際變化

2.1 徑流年際變化的特征值分析

統(tǒng)計(jì)分析西江流域干流上中下游主要控制站1957—2016年平均徑流量實(shí)測資料，見表1。1957—2016年4個水文站多年平均徑流量的偏差系數(shù)較低，均在0.177～0.216之間。其徑流極值比介于2.628～2.760，表明西江干流主要水文站年徑流量年際變化不大。

表1 4個水文站多年平均徑流量和輸沙量及其偏差系數(shù)

參照SL 250—2000《水文情報(bào)預(yù)報(bào)規(guī)范》，將年徑流的距平值百分率Pi劃分為5個級別：枯水年P(guān)i< -20%，偏枯年P(guān)i在-20%～-10%，平水年P(guān)i在-10%～10%，偏豐年P(guān)i在10%～20%，豐水年P(guān)i> 20%。4個水文站年徑流量距平值百分率見圖1。

圖1 4個水文站年徑流量距平值百分率

由圖1可知，干流各水文站年徑流量豐、枯年波動規(guī)律基本一致。1955—1960年4個站平水年和枯水年次數(shù)較多，可視為平水期和枯水期；同理，1960—1980年為平水期和豐水期；20世紀(jì)80年代末—90年代初為枯水期；1990—2000年為豐水期和平水期；2000—2010年回到平水期和枯水期。2015年后除了大湟江口站，其他3站出現(xiàn)平水年和豐水年的趨勢。可概括出西江干流豐枯年份出現(xiàn)規(guī)律可分3個階段：平枯期→平豐期→枯水期。遷江站的豐枯變化最為明顯，2013年P(guān)i≤-40%的高概率，1968和1979年P(guān)i≥50%，變化趨勢最為顯著，因此遷江站是4個站中極易發(fā)生干旱和洪澇災(zāi)害的。

2.2 徑流序列變化的演變趨勢

從站點(diǎn)年徑流量序列過程線可直接觀察到西江流域干流徑流量年際變化不明顯，大部分圍繞于趨勢線上下，較多年份的徑流量均高于多年平均值；其中遷江站多年徑流量總體上有下降趨勢，且年徑流量趨勢線斜率最大，減少趨勢較其他水文控制站更明顯。

根據(jù)趨勢線判斷，大湟江口站、梧州站和高要站年徑流量變化不明顯，結(jié)合Mann-Kendall檢驗(yàn)結(jié)果(表2)，大湟江口的標(biāo)準(zhǔn)統(tǒng)計(jì)量值Z為負(fù)值，|Z|< |Za/2|=1.96，表明該站年徑流量增加趨勢不顯著；同理分析出梧州站和高要站年徑流量在近59年均呈下降趨勢，但下降趨勢并不顯著。

表2 1955—2016年徑流變化趨勢的Mann-Kendall檢驗(yàn)結(jié)果

注：Z為標(biāo)準(zhǔn)統(tǒng)計(jì)量值；α為顯著性參數(shù)。

4個水文站的年徑流量序列過程線見圖2。根據(jù)5年滑動平均分析，遷江站年徑流量在1965年前呈下降幅度稍緩趨勢，1965—1987年為平緩上升階段，1987—2000年遷江站的徑流量先平緩下降后平緩上升，2005年后急劇下降，遷江站在研究時間系列上總體呈下降趨勢。大湟江口站、梧州站以及高要站的年徑流量在20世紀(jì)70年代之前呈緩慢下降變化；第2階段為1968—1985年，此階段年徑流量的變幅較小，呈上升趨勢；第3階段大湟江口站、梧州站和高要站的下降幅度均稍大；第4階段為1995—2004年，3個水文站流量的變化趨勢均為較快上升；第5階段為2004年之后，年徑流量呈加速下降趨勢。結(jié)合一元趨勢線和Mann-Kendall趨勢檢驗(yàn)法，可知大湟江口站、梧州站以及高要站的年徑流量在研究時間系列上的總體變化趨勢下降且不顯著。

圖2 年徑流量序列過程線

2.3 徑流系列突變性分析

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

表3 4個水文站年徑流量滑動T檢驗(yàn)結(jié)果

注：t0為徑流量滑動T檢驗(yàn)初始值。

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

圖4 年徑流量累積距平曲線

由圖3可知，1963年出現(xiàn)大湟江口站年徑流量最小值，同時在圖中存在2003年和1964年2個極值，利用滑動T檢驗(yàn)時大湟江口站年徑流量的統(tǒng)計(jì)量值小于置信度的臨界值(α=0.05)，表明其沒有顯著的突變點(diǎn)，變化趨勢緩慢(表3)。同理，梧州站與高要站均無突變點(diǎn)。

3 西江流域干流輸沙量年際變化分析

3.1 輸沙量時間變化特征分析

西江流域干流輸沙量年際變化較大，各站的輸沙量數(shù)據(jù)與控制站的面積有一定的關(guān)系，輸沙量的空間分布不均。

干流年輸沙量的極值比為15.757～551.877，大湟江口站、梧州站和高要站的極值比分別為15.757、26.310、32.510，遷江站的極值比551.877，變幅最大；4個水文站的最大值出現(xiàn)時段基本相同(圖5)，遷江站、大湟江口站、高要站出現(xiàn)于1983年，梧州站出現(xiàn)于1982年，其主要原因是20世紀(jì)80年代西江流域干流不合理的土地利用、亂砍濫伐以及造成水土流失的工程建設(shè)；而年輸沙量的最小值都在2011—2013年間出現(xiàn)，主要原因是各站上游興修水利工程與實(shí)行水土保持方案以及大規(guī)模采砂現(xiàn)象，直接導(dǎo)致各站輸沙量發(fā)生了巨大變化，總體上為豐水多沙、少水少沙。與年徑流量變化趨勢相比，輸沙量過程曲線均在20世紀(jì)90年代后下降幅度變大，輸沙量年際變化顯著。

圖5 年輸沙量過程線

3.2 輸沙量時間演變趨勢

由圖5可知，遷江站年輸沙量在1964年前呈下降趨勢，1964—1992年平緩上升，1992年后急劇下降；大湟江口站、梧州站和高要站的年輸沙量均呈平緩上升趨勢后急劇下降?？傮w上，西江流域干流遷江站、大湟江口站、梧州站和高要站的年輸沙量總體上均呈下降趨勢。

西江流域干流年輸沙量變化趨勢的Mann-Kendall檢驗(yàn)結(jié)果見表4。4個站年輸沙量統(tǒng)計(jì)值|Z|均大于置信度的臨界值(α=0.05)，4個水文站年輸沙量變化趨勢均為顯著；由于4個站的統(tǒng)計(jì)量值Z均為負(fù)數(shù)，因此年輸沙量呈下降趨勢。

表4 西江流域干流年輸沙量變化趨勢的Mann-Kendall檢驗(yàn)結(jié)果

3.3 輸沙量突變類型分析

3.3.1遷江站年輸沙量的突變點(diǎn)

遷江站年輸沙量有序聚類曲線和累積距平曲線見圖6。遷江站有序聚類值在1988—1997年處于較低水平，在1991年時有序聚類值最低；利用累積距平曲線分析遷江站年輸沙量，得1989—1995年累積距平值較高，1991年出現(xiàn)極大值點(diǎn)；兩種方法結(jié)果皆說明遷江站年輸沙量的跳躍時間點(diǎn)為1991年。通過滑動T檢驗(yàn)法計(jì)算出的輸沙量跳躍時間點(diǎn)也是1991年(表5)，且突變前輸沙量平均值為5 035.59萬t，突變后輸沙量平均值為1 014.59萬t，該點(diǎn)的突變量為4 021.01萬t。由此可知遷江站年輸沙量序列在1991年時發(fā)生顯著突變。

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

表5 西江流域干流年輸沙量突變性的滑動T檢驗(yàn)結(jié)果

注：T為年輸沙量滑動T檢驗(yàn)初始值。

3.3.2大湟江口站年輸沙量的突變點(diǎn)

大湟江口站年輸沙量有序聚類曲線和累積距平曲線見圖7。有序聚類分析得到的最小值出現(xiàn)時間點(diǎn)為1997年，累積距平值在較寬的時間域處于高水平階段，出現(xiàn)最高值5.539 78億t時間點(diǎn)為1998年，且最高值與1997和1999年對應(yīng)的累計(jì)值相差較??；綜合分析年輸沙量的突變時間點(diǎn)為1997年。利用滑動T檢驗(yàn)法再次確定1997年為突變點(diǎn)，且跳躍前后輸沙量的均值分別為6 164.63萬和1 960.63萬t，突變量為4 204萬t。

圖7 大湟江口站輸沙量有序聚類曲線和累積距平曲線

3.3.3梧州站年輸沙量的突變點(diǎn)

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

圖8 梧州站輸沙量有序聚類曲線和累積距平曲線

3.3.4高要站年輸沙量的突變點(diǎn)

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

圖9 高要站輸沙量有序聚類曲線和累積距平曲線

4 水沙變化趨勢分析

4.1 西江站整體水沙變化趨勢

假設(shè)存在一個水文站能夠控制整個西江流域干流，將此虛擬站命名西江站。通過模擬退火算法進(jìn)行年徑流量非線性擬合，年徑流量變化擬合結(jié)果見圖10a)，擬合得到的曲線相關(guān)系數(shù)R2為0.500 83。

由圖10a)可知，西江流域干流年徑流量呈波動變化，變化趨勢較平緩，趨勢變化具體可細(xì)分為：1957—1974年緩慢上升趨勢，1974—1986年緩慢下降趨勢，1986—1998年呈上升趨勢，1998—2011年以較快速度下降，而后至2016年呈急劇上升趨勢。

西江流域干流西江站年輸沙量擬合結(jié)果見圖10b)。經(jīng)過多次調(diào)節(jié)參數(shù)，最終求得年輸沙量模擬曲線的相關(guān)系數(shù)R2為0.87321，擬合度較高。西江流域干流西江站年輸沙量變化趨勢不復(fù)雜，總體上呈下降趨勢：1970年以前呈波動變化，1970年后呈顯著下降。1960—1970年年輸沙量急劇增大。

注：年份1～59代表1957—2016年。

4.2 西江站年輸沙量與年徑流量的關(guān)系

非線性回歸曲線如圖11所示，其擬合相關(guān)系數(shù)R2為0.507 99，說明1957—2015年西江總站年徑流量與輸沙量呈顯著的相關(guān)關(guān)系。

圖11 西江站年徑流量與年輸沙量非線性擬合曲線

模擬退火法計(jì)算的西江站年徑流量與年輸沙量非線性擬合所得的方程為：

y=p1+p2x0.5+p3x+p4x1.5+p5x2+p6x2.5+p7x3

(1)

式中：y為年輸沙量；p1～p7均為參數(shù)；x為年徑流量。

多次調(diào)整參數(shù)，最終確定出圖11擬合曲線參數(shù)值，利用式(1)計(jì)算出1957—2015年輸沙量的回歸值。統(tǒng)計(jì)實(shí)際年輸沙量和計(jì)算值的偏差系數(shù)Cv，其結(jié)果分別為0.529和0.269，說明實(shí)際輸沙量的變化幅度比回歸值的大。

西江站年輸沙量與計(jì)算年輸沙量過程線見圖12。年輸沙量、計(jì)算年輸沙量均與年徑流呈波動狀態(tài)，變化過程基本同步；然而實(shí)際年輸沙量變化量均比年徑流量的變化幅度大，1990年之后實(shí)際輸沙量下降速度顯著，另一方面回歸推算的年輸沙量與徑流量的變化幅度相對較小，而且是豐水多沙、少水少沙的規(guī)律。

圖12 西江站年輸沙量與計(jì)算年輸沙量過程線

有研究表明在氣候條件和下墊面綜合作用下流域徑流和侵蝕產(chǎn)沙[11]。氣候變化可能多方面影響徑流過程和侵蝕產(chǎn)沙過程，其中，最直接的影響因素是降水變化[12]。人類活動如水利工程建設(shè)和水土保持措施等均改變流域下墊面，也影響流域產(chǎn)匯流和侵蝕產(chǎn)沙過程[13]。若流域的輸沙量變化只與氣候變化相關(guān)，則徑流量與輸沙量的變化幅度不大。而由圖11可知，西江流域干流實(shí)際年輸沙量變化幅度遠(yuǎn)大于年徑流量的變幅，其非參數(shù)回歸模擬后計(jì)算得到的輸沙量變化與年徑流變化幅度相對較小，因此，該時間段內(nèi)人類活動作用應(yīng)當(dāng)是西江流域干流西江站輸沙量減少的主要因素。

1992年庫容為34.3億m3的巖灘水庫開始發(fā)揮攔水蓄沙作用，龍灘水電站也于2006年建成并投入使用，1990—2010年以來水利工程建成的總庫容為493.03億m3，占1960—2010年總水利工程建造庫容的94%，可見這些大型水利工程對西江流域干流輸沙量下降的影響力不容忽視[14]。此外1991年國家通過《中華人民共和國水土保持法》，運(yùn)用法律規(guī)范實(shí)施水土保持措施，廣西壯族自治區(qū)人民政府出臺各種控制水土流失的政策，1979—1998年期間，廣西壯族自治區(qū)總造林面積達(dá)3 591 km2；2003—2009年，廣西壯族自治區(qū)內(nèi)每年采砂量為120萬t。

根據(jù)上述資料，西江流域干流輸沙量自1990年后持續(xù)下降原因主要是受人類活動作用。

5 結(jié)論

1)年輸沙量與年徑流量的關(guān)系為多水多沙、少水少沙。4個站豐、枯變化規(guī)律為：平枯期→平豐期→枯水期；遷江站極易發(fā)生洪水與干旱災(zāi)害。

2)4個站年輸沙量變化均為顯著下降趨勢，均存在顯著突變點(diǎn)，遷江、大湟江口、梧州和高要站的年輸沙量突變點(diǎn)為1991、1997、1994和1999年。輸沙量最大值在20世紀(jì)80年代初出現(xiàn)，最小值在2011—2013年間出現(xiàn)。

3)遷江站年徑流量變化呈顯著下降趨勢；其他3個站的年徑流量呈平緩下降趨勢；遷江站年徑流量顯著突變點(diǎn)為2002年，突變前后的多年徑流均值之差為134.3億m3，其他站無顯著突變點(diǎn)。

4)徑流量、輸沙量的非線性擬合曲線的相關(guān)系數(shù)R2分別為0.500 83、0.873 21。年徑流量呈微波動，年際變動幅度較??；年輸沙量總體上呈下降趨勢，先急劇上升后持續(xù)下降，突變點(diǎn)為1990年。模擬退火法計(jì)算徑流量與輸沙量的擬合相關(guān)系數(shù)R2=0.507 99。人類活動作用導(dǎo)致年輸沙量發(fā)生巨大變化。