珠江流域下游近60年輸沙率年際與年內變化特征

2021-02-05 09:06:12楊遠東王永紅蔡斯龍楊清書

水土保持研究 2021年2期

楊遠東, 王永紅, 蔡斯龍, 劉鋒, 楊清書

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

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

1 研究區概況

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

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

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

2 數據來源與研究方法

2.1 數據來源

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

2.2 研究方法

2.2.1 年際變化特征研究方法

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

(1)

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

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

Sw=Wmax/Wmin

(2)

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

(3)

其中：

(4)

在時間序列隨機獨立的假設下，Sk的均值與方差分別為：

(5)

將Sk標準化為：

(6)

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

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

(7)

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

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

(8)

2.2.2 年內變化特征研究方法

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

(9)

RCPyear=arctan(Rx/Ry)

(10)

(11)

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

表1 RCPyear計算各月份包含及代表角度值

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

3 結果與分析

3.1 數據分析

(1) 輸沙率年際變差系數(Sv)與輸沙率年際極值比(Sw)。珠江流域西、北、東江1960—2017年徑流年際變化特征值計算結果見表2，西、北、東江輸沙率變差系數(Sv)值在0.5左右，總體數值較大，說明珠江流域輸沙率年際變率較大，輸沙率穩定性較差，穩定性大小順序為北江大于西江大于東江，西、北、東江輸沙率變差系數(Sv)值與輸沙率年際極值比(Sw)值變化趨勢有所不同，輸沙率年際極值比(Sw)值介于7.15～16.15，北江與西江Sw值較大且較為接近，說明北江與東江更易出現輸沙率的極大值或極小值點，綜合輸沙率變差系數(Sv)值與輸沙率年際極值比(Sw)值可以發現，東江為三江中輸沙率變異性最大的支流。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 輸沙率年際與年內變化原因分析

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

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

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

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

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

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

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