楊遠(yuǎn)東, 王永紅, 蔡斯龍, 劉 鋒, 楊清書

(1.中國海洋大學(xué) 海洋地球科學(xué)學(xué)院 海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,山東 青島 266100; 2.廣東省水文局, 廣州 510150; 3.中山大學(xué) 海洋科學(xué)學(xué)院, 廣州 510275)

河流入海水沙是三角洲的物質(zhì)基礎(chǔ)[1]，河流通過向海洋輸送陸源物質(zhì)從而影響河口、海岸及邊緣海，河流的輸送量占全球陸—海沉積物通量的95%[2]，近年來，隨著人類活動的加劇，全球各主要河流輸沙率均發(fā)生巨大的變化，Milliman[3]指出在過去50 a來歐洲眾多河流入海泥沙通量發(fā)生了明顯的下降。Meade等[4]對Mississippi河研究發(fā)現(xiàn)，其輸沙量在幾十年間下降了75%；Loisel等[5]對湄公河進(jìn)行研究，發(fā)現(xiàn)目前湄公河懸浮泥沙濃度(SSC)現(xiàn)在正在以每年5%的趨勢減少。與國外河流相似，近幾十年河流水沙異變在我國同樣顯著[6]，彭濤等[7]對長江干流寸灘、宜昌、漢口和大通4個(gè)控制站水沙進(jìn)行研究，發(fā)現(xiàn)近60 a來長江干流4個(gè)站點(diǎn)徑流量年際變化相對平穩(wěn)，而輸沙量年際變化較為劇烈，大多數(shù)站點(diǎn)水沙在1968年前后和21世紀(jì)初期發(fā)生了突變：府仁壽等[8]提出近幾十年受氣候變化和人類活動共同影響，長江流域輸沙量大幅減少，由此導(dǎo)致了長江流域中下游嚴(yán)重的生態(tài)環(huán)境問題：張金萍等[9]對黃河三門峽水庫入庫潼關(guān)水文站1919—2015年共計(jì)97 a徑流與輸沙量數(shù)據(jù)進(jìn)行研究，表明潼關(guān)水文站徑流量與輸沙量具有顯著下降趨勢，且在未來一段時(shí)間內(nèi)這種變化趨勢將持續(xù)。珠江作為我國南方地區(qū)的重要河流，其河口處的珠江三角洲由于獨(dú)特的地理位置在中國經(jīng)濟(jì)發(fā)展中起著舉足輕重的作用，近年來，珠江流域及三角洲地區(qū)受人類活動影響顯著，各種研究也較多，但是對這一地區(qū)輸沙率的歷史變化研究較少，本文基于1960—2017年珠江流域下游主要水文站(西江高要、北江石角、東江博羅)年際及年內(nèi)月度輸沙率數(shù)據(jù)，運(yùn)用多種方法分析珠江流域下游地區(qū)輸沙率年際與年內(nèi)變化特征，并分析人類活動對其造成的影響，以期為流域的綜合治理提供參考。

1 研究區(qū)概況

珠江發(fā)源于云南省馬雄山，由西江、北江與東江三大支流構(gòu)成，地理坐標(biāo)位為北緯21°31′—26°49′，東經(jīng)102°14′—115°53′，全長2 400 km，流域面積45萬km2，流經(jīng)云南、貴州、廣西、廣東、湖南、江西6個(gè)省(區(qū))和越南的北部，是中國第二大河，境內(nèi)第三長河。珠江流域地處亞熱帶，北回歸線橫貫流域的中部，氣候溫和多雨，多年平均溫度為14～22℃，多年平均降雨量1 200～2 200 mm，多年平均徑流量3 100多億m3(1954—2016年)，多年平均輸沙率2 800多萬t(1954—2016年)，徑流水沙季節(jié)分配很不均勻，常年4—9月份洪水期多年平均徑流量約占全年的80%，輸沙率占91%～95%。

珠江流域西北兩江分別在廣東省三水市思賢滘東西兩側(cè)、東江在廣東省東莞市石龍鎮(zhèn)匯入珠江三角洲。高要(北緯23.02°, 東經(jīng)112.27°)、石角(北緯23.11°，東經(jīng)112.52°)、博羅(北緯23.11°，東經(jīng)114.17°)分別為進(jìn)入三角洲之前西、北、東江最主要的水文控制站(圖1)。

圖1 珠江流域及高要、石角、博羅水文站位置

2 數(shù)據(jù)來源與研究方法

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

本文收集了珠江流域進(jìn)入珠江三角洲之前西、北、東江的主要水文控制站——高要站、石角站、博羅站1960—2017年年度及月度輸沙率數(shù)據(jù)，長時(shí)間序列輸沙率數(shù)據(jù)來源于廣東省水文局(歷史測量數(shù)據(jù))，中山大學(xué)(歷史保存數(shù)據(jù))與《中國河流泥沙公報(bào)》(2002—2017年)。珠江流域主要大壩建設(shè)數(shù)據(jù)來源于中國大壩工程學(xué)會網(wǎng)站。

2.2 研究方法

2.2.1 年際變化特征研究方法

(1) 輸沙率年際變差系數(shù)(Sv)與輸沙率年際極值比(Sw)。為研究輸沙率的年際變化特征，與徑流年際變化類似[徑流年際變異性常用變差系數(shù)(Cv)與年際極值比(w)來表示]，引入輸沙率變差系數(shù)(Sv)與輸沙率年際極值比(Sw)，輸沙率變差系數(shù)(Sv)用來說明水文變量長期變化的穩(wěn)定程度，其值為輸沙率樣本序列的標(biāo)準(zhǔn)差與算術(shù)平均值的比值，Sv值越大，年輸沙率的年際變化越大，可以作為反映輸沙率年際變化不均勻性的重要指標(biāo)。輸沙率變差系數(shù)(Sv)計(jì)算公式如下:

(1)

式中:ki為模比系數(shù),即第i年的年輸沙率與平均年輸沙率的比值;n為系列長度。

年際極值比(Sw)是指徑流量變化的絕對值比例，其大小用多年最大(Wmax)、最小(Wmin)年輸沙率比值表示，年際極值比計(jì)算公式如下：

Sw=Wmax/Wmin

(2)

(2) M-K趨勢分析法與滑動平均法。M-K趨勢分析法(Mann-Kendall法)是一種非參數(shù)統(tǒng)計(jì)檢驗(yàn)方法，其優(yōu)點(diǎn)是樣本不需要遵循一定的分布，也不受少數(shù)異常值的干擾，更適用于類型變量與順序變量，計(jì)算也比較簡便，最初由曼和肯德爾提出原理并發(fā)展了這一方法?，F(xiàn)已成為世界氣象組織推薦并已廣泛使用的非參數(shù)檢驗(yàn)方法，適用于水文、氣象等非正態(tài)分布的數(shù)據(jù)。其計(jì)算公式見公式(3)—(6)，設(shè)徑流序列為x1，x2，…，xn；Sk表示第i個(gè)樣本xi>xj(1≤j≤i)的累計(jì)數(shù)。

(3)

其中：

(4)

在時(shí)間序列隨機(jī)獨(dú)立的假設(shè)下，Sk的均值與方差分別為：

(5)

將Sk標(biāo)準(zhǔn)化為：

(6)

式中:UF1=0,給定顯著性水平α,若|UF|>Uα,表明序列存在明顯的趨勢變化。將此方法引用到反序列,序列順序變?yōu)閤n,xn-1,…,x1,反序列標(biāo)準(zhǔn)化由UBk表示;UBk,UFk可組成兩條曲線,在這里設(shè)定α值為0.05,臨界值U0.05為±1.96。

滑動平均法因其簡單直觀性，在水文學(xué)研究中應(yīng)用較為廣泛。其主要思路是將時(shí)間序列x1，x2，…，xn的幾個(gè)前期值和后期值取平均，求出新的序列yt。其計(jì)算公式如公式(7)[10]：

(7)

式中:yt為新序列;k為滑動平均尺度;xt+i為參與此次生成新序列的舊序列值;k=1時(shí)為3 a滑動平均值;k=2時(shí)為5 a滑動平均值;若Xn有趨勢性,選擇合適的k,即可通過將其趨勢清晰地表示出來。

(3) 累積距平法。為了更直觀地觀察輸沙率的年度變化特征，引入輸沙率累計(jì)距平方法[11-12]，輸沙率累積距平計(jì)算方法如下[13]：

(8)

2.2.2 年內(nèi)變化特征研究方法

(1) 集中度與集中期法。目前對于長序列輸沙率年內(nèi)變化的分析方法較少，本文借鑒長序列徑流年內(nèi)分配的分析方法，對珠江三角洲西江高要站、北江石角站、東江博羅站近幾十年輸沙率的年內(nèi)變化特征進(jìn)行分析，具體原理與計(jì)算過程見文獻(xiàn)[14]。計(jì)算公式如公式(9)—(11)：

(9)

RCPyear=arctan(Rx/Ry)

(10)

(11)

式中:RCDyear為年輸沙率集中度;RCPyear為年輸沙率集中期;Ryear為年輸沙率;Rx,Ry分別為12個(gè)月的分量之和所構(gòu)成的水平、垂直分量;ri為第i月的輸沙率;θi為第i月輸沙率的矢量角度。RCDyear越大表示輸沙率越分散,越小輸沙率越集中;RCPyear表示1年中最大輸沙率出現(xiàn)的時(shí)間,表1為各個(gè)月份包含及代表角度值。

表1 RCPyear計(jì)算各月份包含及代表角度值

(2) 等值線圖法。等值線圖又稱等量線圖，是以相等數(shù)值點(diǎn)的連線表示連續(xù)分布且逐漸變化的數(shù)量特征的一種圖型，具有更強(qiáng)的直觀性。

3 結(jié)果與分析

3.1 數(shù)據(jù)分析

(1) 輸沙率年際變差系數(shù)(Sv)與輸沙率年際極值比(Sw)。珠江流域西、北、東江1960—2017年徑流年際變化特征值計(jì)算結(jié)果見表2，西、北、東江輸沙率變差系數(shù)(Sv)值在0.5左右，總體數(shù)值較大，說明珠江流域輸沙率年際變率較大，輸沙率穩(wěn)定性較差，穩(wěn)定性大小順序?yàn)楸苯笥谖鹘笥跂|江，西、北、東江輸沙率變差系數(shù)(Sv)值與輸沙率年際極值比(Sw)值變化趨勢有所不同，輸沙率年際極值比(Sw)值介于7.15～16.15，北江與西江Sw值較大且較為接近，說明北江與東江更易出現(xiàn)輸沙率的極大值或極小值點(diǎn)，綜合輸沙率變差系數(shù)(Sv)值與輸沙率年際極值比(Sw)值可以發(fā)現(xiàn)，東江為三江中輸沙率變異性最大的支流。

表2 1960-2017年珠江流域西、北、東江輸沙率年際變化特征

(2) M-K趨勢分析法與滑動平均法。用于診斷連續(xù)水文序列趨勢性方法較多，其中應(yīng)用最廣范、最為簡便直觀的方法為M-K趨勢分析法與滑動平均法，本文分別用這兩種方法對珠江流域西、北、東江1960—2017年輸沙率趨勢性進(jìn)行分析。

由圖2可知，高要、石角、博羅站歷史輸沙率變化幅度較大，西江、東江輸沙率明顯下降，通過M-K趨勢檢驗(yàn)法檢測發(fā)現(xiàn)，西江高要站輸沙率變化曲線UF與UB在2001年前后相交，交點(diǎn)位于0.05置信區(qū)間內(nèi)，證明高要站輸沙率在2001年開始出現(xiàn)下降趨勢，UF在2004年之后始終位于0.05置信區(qū)間之外，輸沙率持續(xù)下降。高要站2001—2016年平均輸沙率僅為1970—2000年的30%，在流量基本保持穩(wěn)定的情況下，反映了西江含沙量的快速下降；石角站2000年之后輸沙率雖有下降，但下降趨勢不明顯，UF始終位于0.05置信區(qū)間之內(nèi)，2001—2016年平均輸沙率為1960—2000年的80%。東江博羅站自1987年前后開始下降，1990年之后UF始終位于0.05置信區(qū)間之外，下降趨勢較明顯，2001—2016年平均輸沙率為1960—2000年的65%。

圖2 1960-2017年高要、石角、博羅輸沙率M-K檢驗(yàn)

珠江流域高要、石角、博羅1960—2017年輸沙率3 a滑動平均值與5 a滑動平均值見圖3?？梢姡咭⑹恰⒉┝_水文站1960—2017年輸沙率3 a滑動平均值與5 a滑動平均值較為相似，近幾十年均呈現(xiàn)“上升—下降—上升”的波動過程，但3站輸沙率總體呈下降趨勢，尤其以西江高要站下降最為迅速，其次為東江博羅站，下降最緩慢的為北江石角站。對高要、石角、博羅1960—2017年輸沙率數(shù)據(jù)作線性回歸分析，研究其58 a間總體變化特征，線性回歸方程見圖3，高要、石角、博羅線性分析相關(guān)系數(shù)r值分別為0.33，0.04，0.34，由線性回歸統(tǒng)計(jì)學(xué)可知，當(dāng)α=0.05，n=58時(shí)，查表得相關(guān)系數(shù)r=0.259，3個(gè)流域中西江、東江|r|值大于r0.05，證明這兩個(gè)支流輸沙率下降趨勢明顯，北江|r|值小于r0.05，輸沙率下降趨勢不明顯。綜合圖3逐年輸沙率與3 a，5 a輸沙率滑動平均值可知，西江高要站輸沙率下降始于1997年前后，北江石角站下降趨勢不明顯，東江博羅站輸沙率下降始于1984年前后。

綜合比較滑動平均值(3 a，5 a)法與M-K趨勢檢驗(yàn)法，可以發(fā)現(xiàn)其反映的總體規(guī)律基本相似，兩種方法都顯示西江高要站與東江博羅站輸沙率存在下降趨勢，而北江石角站輸沙率總體變化趨勢不大。但其變化的時(shí)間存在差異，M-K趨勢檢驗(yàn)法檢測西江高要站輸沙率在2001年開始出現(xiàn)下降趨勢，東江博羅站自1987年前后開始下降，這分別比滑動平均值法遲3～4 a，其主要原因是由于滑動平均法只是將3 a或5 a進(jìn)行滑動平均，使曲線變得平緩，但并不能避免下降過程中會出現(xiàn)小的波動，從這一方面來說，M-K趨勢在統(tǒng)計(jì)學(xué)上講意義更大，同時(shí)滑動平均值法可以起到很好的檢驗(yàn)作用，綜合分析兩種結(jié)果，認(rèn)為1960—2017年珠江三角洲頂點(diǎn)處西江、東江輸沙率下降趨勢顯著，起開始下降年份分別為2001年、1987年，北江下降趨勢不明顯。

圖3 1960-2017高要、石角、博羅年輸沙率及3 a，5 a滑動平均值

(3) 累積距平法。高要、石角、博羅3站1960—2017年輸沙率累計(jì)距平見圖4，高要、博羅累積距平形狀較為相似，都經(jīng)歷一個(gè)上升—下降過程，其上升、下降的轉(zhuǎn)折點(diǎn)分別在1998年、1987年前后，而石角站累積距平曲線形狀較為曲折，存在多個(gè)轉(zhuǎn)折點(diǎn)。

(4) 集中度與集中期法。珠江流域高要、石角、博羅水文站1960—2017年10 a尺度的輸沙率集中度與集中期計(jì)算結(jié)果見表3，由表3可知，輸沙率的年內(nèi)集中程度為西江高于北江高于東江，西江高要站年內(nèi)輸沙率集中度在1960—2009年較為穩(wěn)定，基本維持在0.8左右，2010—2017年迅速減小為0.7，北江石角站年內(nèi)輸沙率集中度自1960年以來基本呈增大趨勢，東江博羅站年內(nèi)輸沙率集中度自1960年以來基本呈減小趨勢，其RCDyear值較西江高要站與北江石角站略小。

圖4 1960-2017年高要、石角、博羅輸沙率累積距平

表3 高要、石角、博羅水文站年輸沙率集中度與集中期

根據(jù)年輸沙率集中期對應(yīng)的角度，可以發(fā)現(xiàn)西江高要站年內(nèi)輸沙率主要集中在7月份，但近幾十年其輸沙率集中度對應(yīng)數(shù)值在緩慢變小，說明其最大輸沙率對應(yīng)的月份不變，但在最大輸沙月內(nèi)最大輸沙日期向前移動；北江石角站輸沙率集中期在5月、6月份交替出現(xiàn)，1970—1989年期間最大輸沙發(fā)生在5月份，其余時(shí)間最大輸沙發(fā)生在6月份，近幾十年在季節(jié)尺度無明顯向前或向后移動的趨勢；東江博羅站輸沙率集中期在6月、7月份交替出現(xiàn)，近幾十年其輸沙率集中度對應(yīng)數(shù)值在緩慢變大，說明其最大輸沙率對應(yīng)的時(shí)間在向后移動，1960—1999年期間最大輸沙基本發(fā)生在6月份，此后最大輸沙多發(fā)生在7月份。

(5) 等值線圖法。對1960—2017年高要、石角、博羅水文站月度輸沙率作等值線圖(圖5)，從等值線圖中可以清楚地看出輸沙率年內(nèi)的分散程度與最大、最小徑流出現(xiàn)的時(shí)間。輸沙率等值線圖與輸沙率集中度與集中期表格演變階段具有一致性，等值線圖中年內(nèi)輸沙率較分散年份，輸沙率集中度往往也較低，反之亦然，而輸沙率集中期角度的波動與等值線圖最大值的上下移動也有很好的對應(yīng)性。

圖5 1960-2017年高要、石角、博羅輸沙率等值線

3.2 輸沙率年際與年內(nèi)變化原因分析

3.2.1 輸沙率年際變化原因分析 影響河流年際輸沙率變化的原因眾多，而珠江三角洲輸沙率年際變化受人類活動影響程度較大[15-16]，對珠江流域西、北、東江徑流量近幾十年分析，發(fā)現(xiàn)其徑流量并未發(fā)生較大變化[17]，因此造成輸沙率劇烈變化的主要原因是河流含沙量的變化。

圖6為高要、石角、博羅1960—2017年河流含沙量的變化圖。由圖6可知，近幾十年西江高要站、東江博羅站河流含沙量呈下降趨勢，且西江下降率高于東江，北江石角站含沙量變化不大，含沙量的變化趨勢與輸沙率變化趨勢相一致，且變化時(shí)間也基本吻合，西江高要站含沙量在1997年開始下降，東江博羅站含沙量在1987年前后開始下降。

圖6 1960-2017年高要、石角、博羅站含沙量變化

人類活動主要通過影響河流含沙量的變化從而影響流域的輸沙率，通過對國內(nèi)外人類活動對流域輸沙率的干預(yù)狀況研究發(fā)現(xiàn)，人類活動主要通過兩種方式改變河流的含沙量，一是流域植被覆蓋率的變化，植被覆蓋率影響土壤侵蝕，從而對河流含沙量產(chǎn)生影響，而植被覆蓋率與國家或地區(qū)有關(guān)的環(huán)境政策與流域的綜合管理或治理情況有關(guān)；二是流域大壩、水庫的修建，大壩、水庫在完成其蓄水的同時(shí)，也使大量的泥沙在大壩前或水庫中積聚，從而改變下游地區(qū)的含沙量。

據(jù)珠江水利委員會水農(nóng)處通過遙感資料統(tǒng)計(jì)，1984—1988年相比1954年，水土流失土地的面積增大了38.9%，從而使大量泥沙下泄，增大了河流的含沙量，這與圖5對應(yīng)較好，由于西江來沙量占珠江流域總的來沙量的絕大部分，在這兩個(gè)時(shí)間點(diǎn)其變化特征也最大，北江偏小，東江雖呈現(xiàn)下降趨勢，但由于其輸沙量在珠江流域所占比例較小，不影響總體的變化。1982年全國第4次水土保持工作會議以后，水土保持治理工作逐漸開展，這是珠江流域含沙量下降的一個(gè)原因。

但隨著新建大壩、水庫數(shù)量及庫容量的增大，這些大工程也越來越成為影響流域河流含沙量的一個(gè)重要因素[18]，本文統(tǒng)計(jì)了珠江流域較大規(guī)模大壩、水庫的修建情況，包括大壩、水庫的庫容及建成時(shí)間，具體數(shù)據(jù)見表4。西江流域較大的水利工程一般在1997年之后，1997年建成的天勝橋水電站其庫容量高于前期建成的較大規(guī)模的水庫庫容量之和，而2006年建成的龍灘水電站，是當(dāng)今珠江流域最大的水利工程，同時(shí)也是繼三峽水電站之后國內(nèi)第二大水電站，其庫容量高達(dá)273億m3，對于泥沙的積聚作用巨大[18]。由圖5可以看出，含沙量在這兩個(gè)時(shí)間節(jié)點(diǎn)都有明顯的下降過程，可見大壩、水庫的修建對于含沙量的減小具有明顯的作用。北江水電站規(guī)模較小，且修建時(shí)間較早，大壩修建引起的含沙量的降低與水土流失造成的含沙量的增大互相抵消一部分，導(dǎo)致含沙量變化不大，而東江1969年修建了新豐江水電站，庫容量高達(dá)139億m3，相比于東江較小的輸沙率占比(圖7)，其庫容量相對較大，在1969年之后博羅站含沙量有一個(gè)下降過程，此后隨著楓樹壩與白盆珠水電站的建成，其含沙量又經(jīng)歷了一個(gè)下降過程。

3.2.2 輸沙率年內(nèi)變化原因分析 珠江流域高要、石角、博羅輸沙率年內(nèi)變化主要與兩個(gè)因素有關(guān)，一是徑流的年內(nèi)分配，二是年內(nèi)河流含沙量的變化。通過對比高要、石角、博羅1960—2017年徑流年內(nèi)集中度與集中期[17]，發(fā)現(xiàn)高要站輸沙率集中度與徑流集中度變化特征相似，兩者都表現(xiàn)為1960—2009年較為穩(wěn)定，2010—2017年迅速減小，徑流集中期與輸沙率集中期基本相似，但輸沙率集中度遠(yuǎn)遠(yuǎn)大于徑流集中度，表現(xiàn)為輸沙率年內(nèi)更加集中，北江石角、東江博羅站輸沙率集中度、集中期表現(xiàn)為與徑流集中度、集中期密切相關(guān)。而分析河流年內(nèi)含沙量分析，發(fā)現(xiàn)雖然整體含沙量呈下降趨勢，但每個(gè)月相對全年所占比例變化不大，說明輸沙率的年內(nèi)變化主要與徑流的年內(nèi)變化有關(guān)，而徑流的年內(nèi)變化又與上文提到的大壩、水庫的修建有關(guān)。

表4 珠江流域較大規(guī)模水庫建設(shè)情況

圖7 1960-2017年高要、石角、博羅分沙比

4 結(jié) 論

(1) 1960—2017年珠江流域主要水文站高要、石角、博羅站輸沙率變差系數(shù)(Sv)值在0.5左右，總體數(shù)值較大，說明珠江流域輸沙率年際變率較大，輸沙率穩(wěn)定性較差，而輸沙率年際極值比顯示北江與東江更易出現(xiàn)輸沙率的極大值或極小值點(diǎn)，綜合分析東江為三江中輸沙率變異性最大的支流。

(2) M-K趨勢分析、滑動平均法、累積距平法均顯示1960—2017年珠江流域西江、東江輸沙率下降明顯，北江輸沙率變化不大，高要站輸沙率在2001年開始出現(xiàn)下降趨勢，東江博羅站自1987年前后開始下降。

(3) 珠江流域西、北、東江輸沙率年內(nèi)分配極不均勻，輸沙率的年內(nèi)集中程度為西江高于北江高于東江，西江高要站年內(nèi)輸沙率主要集中在7月份，北江石角站輸沙率集中期在5月、6月份交替出現(xiàn)，東江博羅站輸沙率集中期在6月、7月份交替出現(xiàn)。

(4) 近幾十年人類活動在珠江流域輸沙率變化方面起著很大的作用，主要表現(xiàn)在流域大壩、水庫的修建與由于政策與經(jīng)濟(jì)原因?qū)е碌乃亮魇У某潭炔町悺?/p>