近70年珠江水沙变化特征及人类活动影响因素分析

杜以超，罗孝文*,3,4，王 峻，崔家馨

(1.成都理工大学 沉积地质研究院，四川 成都 610059; 2.自然资源部第二海洋研究所，自然资源部海底科学重点实验室，浙江 杭州 310012；

3.自然资源部海洋空间资源管理技术重点实验室，浙江 杭州 310012；

4.浙江省海洋科学院，浙江 杭州 310012)

河流是陆源物质(如淡水、泥沙、化学物质)向海洋输送的主要途径，是全球物质循环的重要环节。据统计，每年有3.6万km3淡水和13～20 Gt沉积物经河流输运到海洋[1-3]。近几十年来，人类活动对河流系统的影响日益加剧，导致许多大河的入海泥沙通量显著减少[4-5]。河流入海泥沙是河口三角洲形成的重要物质基础，泥沙通量的变化不仅影响流域地表形态，也对河口三角洲演变及其附近海岸带的环境产生影响[6]。由于入海泥沙通量的减少，有些河流的河口三角洲已发生侵蚀，例如长江[7]、黄河[8]、尼罗河[9]等。因此，研究人类活动对河流入海泥沙通量的影响对流域河道整治与海岸带环境监测至关重要。

珠江是中国华南地区最大的河系，其下游的珠江三角洲是我国经济最发达的地区之一。近几十年，随着珠江流域经济的迅速发展，受土地开垦、水库建设、水土保持、退耕还林等人类活动的影响，珠江入海泥沙通量已经发生了显著变化。一方面流域内大坝、水库等水利工程的建设导致珠江入海泥沙通量减少，如DAI et al[10]研究表明水土流失面积扩大是20世纪80年代之前珠江入海泥沙通量增加的主要原因，而90年代水库建设后，珠江输沙量明显下降。另一方面土地利用方式和植被覆盖的变化也会影响土壤侵蚀强度和河流输沙量[11-13]，如高海东 等[14]研究发现2000—2011年黄河流域的河龙区间输沙减少量的54%是由植被恢复引起的。目前，对珠江流域植被覆盖变化与河流水沙关系的研究相对缺乏。戴仕宝 等[15]研究发现珠江流域的水土流失面积在1995年与2004年基本未发生变化，而输沙量呈下降趋势，认为1995—2004年珠江入海泥沙的减少与流域水土保持关系不大。而TAN et al[16]通过21世纪初珠江流域植被覆盖面积较20世纪90年代有所减少的现象，认为森林砍伐增多导致西江主要支流柳江输沙量增加。

本文基于珠江流域1954—2019年的年径流量与年输沙量数据以及2000—2019年珠江流域的年降水量与归一化植被指数(Normalized Difference Vegetation Index，NDVI)对珠江流域近70年来的输沙量变化特征以及2000—2019年流域植被覆盖变化与输沙量变化的关系进行研究，讨论了不同时期人类活动对珠江输沙量的影响。

珠江发源于云南马雄山，全长2 400 km，总流域面积为45.37万km2(图1)，流经我国的云南、贵州、广西、广东、湖南和江西6个省(区)及越南北部，经八大口门注入南海，是中国径流量第二大的河流[10]。珠江流域由西江、北江、东江及三角洲诸河四个水系组成。西江是珠江的主干流，全长2 214 km，流域面积35.15万km2，由南盘江、红水河、黔江、浔江及西江等河段组成，其多年平均径流量和平均输沙量分别占珠江流域入海总径流量与总输沙量的77%和89%[10]。北江与东江的河流长度分别为468 km与520 km，流域面积分别是3.84万km2与2.53万km2[11]。高要站、石角站、博罗站分别为西江、北江、东江最主要的水文控制站，3个站实测径流量与输沙量之和分别占珠江入海总径流量与总输沙量的80%和95%左右[12]。珠江流域建有多个不同规模的水库，主要水库信息统计见表1。

图1 珠江流域主要水系、水文站、气象站和水库分布图[17]Fig.1 Distribution data of major river systems, hydrological stations, meteorological stations and reservoirs in the Pearl River Basin[17]

表1 珠江流域的主要水库Tab.1 Large dams and reservoirs in the Pearl River Basin

2.1 数据来源与处理

本文收集了1954—2019年高要站、石角站和博罗站3个水文站的年径流量和年输沙量资料，数据来源于中国河流泥沙公报[19]。气象数据来自中国气象局科学数据共享服务网(http:// www.cma.gov.cn/)，收集了1954—2019年珠江流域内38个气象站的年降水量数据。NDVI数据来源于美国国家航空航天局的中分辨率成像光谱仪(Moderate-resolution Imaging Spectroradiometer, MODIS)的植被指数产品数据MOD13Q1-NDVI，空间分辨率为500 m，时间分辨率为16 d，时间范围为2000年3月—2019年12月，与珠江防护林第二期和第三期工程基本同步。

NDVI是表征植被生长状况的指标，代表了地表植被覆盖变化。采用最大值合成法(Maximum Value Composites, MVC)获取2000—2019年逐年最大化NDVI。MVC法对云层、双向反射噪声移除效果明显，运算效率高，是目前使用最广泛的植被指数合成算法[20]。年最大化NDVI可作为流域植被覆盖指标，但可能存在受极端天气的影响，个别地区NDVI值偏高的问题[21-22]。

2.2 研究方法

2.2.1 Mann-Kendall 非参数检验与Theil-Sen Median趋势分析

Mann-Kendall非参数检验[23](简称M-K检验)可用于径流量和输沙量时间序列的突变分析，同时可用于检验NDVI时间序列数据变化趋势的显著性，计算简便且不受少数异常值干扰，广泛应用于水文、气象等领域。M-K检验用于径流量和输沙量时间序列的突变分析时，可将样本数量为n的时间序列(x1,x2, …,xn)组成一个秩序列Sk(表示第i个样本大于第j个样本的累计数)：

(1)

(2)

式中：xi和xj分别代表第i年和第j年的数值。Sk的均值E(Sk)和方差Var(Sk)分别为：

(3)

(4)

在时间序列随机独立的假设下，定义时间序列统计量UF:

(5)

当UF>0时，表明序列呈上升趋势；
UF<0时，表明序列呈下降趋势。取UF与置信水平α(本文取0.05)对应的临界值(本文为±1.96)进行对比，当 |UF|>1.96，表明趋势变化显著。用公式(5)计算时间序列的逆序列统计量UB，当UB=-UF且曲线交点在95%置信区间内，表明该时间节点可能发生突变。若时间区间内有多个交点，且UF曲线未发生显著变化则不能完全判断为突变点。具体方法参考文献[24-25]。

M-K检验用于检验NDVI时间序列数据变化趋势的显著性时，在0.05置信水平上划分显著性变化，判定趋势变化的显著性。设定2000—2019年的流域植被指数为NDVIt(t=2000,2001,…,2019)，用下式计算统计检验值Z：

(6)

(7)

(8)

(9)

式中：
NDVIt和NDVIb分别表示像元t年与b年的NDVI值，n为时间序列长度，sgn为符号函数。Z值为正值，表明为增长趋势；
Z值为负值，表明为下降趋势。

Theil-Sen Median趋势分析是一种稳健的非参数计算方法[26]，该方法已被广泛应用于NDVI时间序列的分析中[27-28]。用Theil-Sen Median趋势分析法计算n(n-1)/ 2个数据组合的斜率中位数(SNDVI)，其计算公式为：

(10)

若SNDVI>0，则NDVI时间序列呈上升趋势，反之呈下降趋势。由于基本不存在SNDVI严格等于0的区域，因此取Theil-Sen Median分析结果等于±0.000 5作为分界值：SNDVI≥0.000 5 为植被改善区;SNDVI≤-0.000 5 为植被退化区;-0.000 5

将M-K检验在0.05置信水平上的显著性检验结果划分为变化显著 (Z>1.96或Z<-1.96)和变化不显著(-1.96≤Z≤1.96)。将Theil-Sen Median趋势分析的分级结果与M-K检验的Z值结合可以有效反映珠江流域NDVI变化的空间分布特征：当SNDVI≥0.000 5且Z>1.96时，植被覆盖明显改善；
当SNDVI≥0.000 5且-1.96≤Z≤1.96时，植被覆盖轻微改善；
当-0.000 5当SNDVI≤-0.000 5且-1.96≤Z≤1.96时，植被覆盖轻微退化；
当SNDVI≤-0.000 5且Z<-1.96时，植被覆盖严重退化。

2.2.2 累积距平法

绘制输沙量与径流量的累积距平曲线，根据曲线变化趋势可大致区分径流量与输沙量的阶段性变化[30]。年累积距平值的计算公式如下：

(11)

(12)

2.2.3 双累积曲线分析

双累积曲线是一种检验两种参数一致性及其变化趋势最直观的分析方法，如今广泛应用于水文与气象因素的趋势性变化及其强度分析[31-32]。该方法基于水文要素之间的统计关系评价人类活动对水文过程的影响。在直角坐标系中绘制同时期径流量与输沙量的连续累积关系曲线，当曲线呈线性关系时表明该时期流域产汇流条件相对一致，如果斜率发生改变则表示两变量之间的关系受人类活动的影响。

图2为1954—2019年西江(高要站)年径流量与年输沙量的M-K检验统计结果与累积距平曲线。由 图2a 可知，西江(高要站)的年径流量变化大致分为5个阶段：其中1954—1967年、1989—1994年、2005—2019年的年径流量UF曲线为负值，年径流量呈下降趋势；
1968—1988年和1995—2004年的年径流量UF曲线为正值，年径流量呈缓慢上升趋势。西江(高要站)年径流量的UF曲线与UB曲线有多个交点，但UF曲线未超过95%的置信区间，因此不存在显著突变点。累积距平曲线显示1967年、1983年、1992年、1999年和2014年为年径流量的丰枯阶段性趋势变化的转折年，但变化趋势并不显著。由图2b可知，西江(高要站)的年输沙量UF曲线在1954—1968年和2000—2019年为负值，表明该时段年输沙量呈减少趋势；
UF曲线在1969—1999年为正值，表明该时段的年输沙量呈缓慢增加趋势。UF曲线与UB曲线相交于2003年，且在2007年UF曲线超过临界值(-1.96)，表明年输沙量显著下降，由此可以判断年输沙量变化趋势于2003年发生了突变。西江(高要站)1954—2002年的年平均输沙量为69.3 Mt，2003—2019年的年平均输沙量为19.9 Mt，突变后年输沙量下降了71.30%。累积距平曲线显示1998年为年输沙量丰枯阶段性趋势变化的转折年。

图2 1954—2019年西江(高要站)年径流量(a)与年输沙量(b)的M-K检验统计结果与累积距平曲线Fig.2 Mann-Kendall test statistic result and accumulative anomaly curves of the annual runoff (a) and sediment discharge (b) in the Xijiang River (Gaoyao Station) from 1954 to 2019

图3为1954—2019年北江(石角站)年径流量与年输沙量的M-K检验统计结果与累积距平曲线。由图3a可知，北江(石角站)年径流量UF曲线变化复杂，年径流量变化有4个缓慢上升期(1959—1962年、1975—1987年、1994—2008年和2016—2019年)和4个缓慢下降期(1954—1958年、1963—1974年、1988—1993年和2009—2015年)，变化趋势并不显著。北江(石角站)的年径流量UF曲线与UB曲线有多个交点，但均没有超过临界值，因此没有显著突变点。结合累积距平曲线判断1962年、1972年、1984年、1991年和2002年是年径流量丰枯阶段性趋势变化的转折点。由图3b可知，1954—2000年的年输沙量UF曲线基本呈波动上升趋势，从2000年至2019年呈下降趋势，但均未超过95%的置信区间，虽然UF曲线与UB曲线于1998年、2005年、2016年相交，但不能判断这些年份发生突变。累积距平曲线显示1972年、1984年、1998年和2005年为年输沙量丰枯阶段性趋势变化的转折年。

图3 1954—2019年北江(石角站)年径流量(a)与年输沙量(b)M-K的检验统计结果与累积距平曲线Fig.3 Mann-Kendall test statistic result and accumulative anomaly curves of the annual runoff (a) and sediment discharge (b) in the Beijiang River (Shijiao Station) from 1954 to 2019

东江(博罗站)年径流量的变化趋势与突变分析结果(图4a)显示1954—2010年东江年径流量主要呈波动上升趋势，但变化趋势并不显著，2010年以后UF曲线在0值上下波动，年径流量不存在显著突变点。年径流量累积距平曲线显示年径流量变化趋势转折出现在1962年、1972年、1985年、2004年和2008年。从图4b可知，东江(博罗站)年输沙量从1969年左右开始减少，UF曲线与UB曲线相交于1991年，且于1994年UF曲线突破临界值(-1.96)，说明输沙量呈显著下降趋势。因此可以判断年输沙量于1991年发生突变，1991—2019年东江(博罗站)的年均输沙量相较1954—1990年下降了48.88%。

图4 1954—2019年东江(博罗站)年径流量(a)与年输沙量(b)M-K的检验统计结果与累积距平曲线Fig.4 Mann-Kendall test statistic result and accumulative anomaly curves of the annual runoff (a) and sediment discharge (b) in the Dongjiang River (Boluo Station) from 1954 to 2019

4.1 植被覆盖的时空变化

近20年来，西江、北江和东江流域的NDVI整体呈增长趋势，表明珠江流域植被覆盖率增加。2019年西江流域的NDVI值比2000年增加了0.065，其中在2000—2007年和2012—2015年增长较快，在2008—2011年和2016—2019年变化较平稳(图5a)。北江流域NDVI值在2000—2008年呈波动变化，2008年后明显增加，并于2015年后趋于稳定(图5b)。东江流域的NDVI值在2008年以后呈波动式上升趋势(图5c)。

图5 2000—2019年珠江流域NDVI年际变化Fig.5 Inter-annual variation of NDVI in the Pearl River Basin from 2000 to 2019

将Theil-Sen Median趋势分析的分级结果与M-K检验两种分析结果叠加得到珠江流域像元尺度上的NDVI变化趋势图(图6)，并将植被变化分为5种类型(表2)。据统计，2000—2019年珠江流域植被覆盖状况改善的区域占珠江流域总面积的89.27%。由于不同流域的地表环境与植被恢复情况具有差异，将珠江三大流域植被覆盖变化分开统计，如表2所示，西江、北江和东江流域植被覆盖状况改善的区域分别为90.97%，90.77%和89.67%，整个珠江流域植被覆盖以改善为主。

表2 西江、北江和东江流域的植被覆盖变化统计Tab.2 Statistics of vegetation cover in Xijiang River Basin, Beijiang River Basin and Dongjiang River Basin

从图6也可以看出，过去20 a，珠江流域地表植被改善区域面积远大于植被退化区域面积。虽然在西江上游的昆明、曲靖、玉溪、红河哈尼族彝族自治州等地区出现轻微退化甚至严重退化区域，但大部分地区植被覆盖明显改善;西江中游的柳州、桂林、南宁、贵港等地出现了小面积植被退化和稳定不变的区域;严重退化的区域主要分布在下游珠江三角洲地区的广州、佛山、江门、中山等地。西江上游是珠江泥沙供给的主要区域，下游地区不是河流泥沙的主要供给区，因此下游的植被退化现象对输沙量变化影响有限。

图6 2000—2019年珠江流域植被覆盖变化分布Fig.6 Distribution of vegetation cover of the Pearl River Basin from 2000 to 2019

4.2 植被覆盖与水沙的相关性

通过双累积曲线法分析2000—2019年植被覆盖变化对水沙的影响。图7中3个水文站的年降水量-径流量累积曲线几乎均为一条直线，说明径流量主要受流域降水影响。

图7 2000—2019年珠江流域年降水量-径流量和年降水量-输沙量双累积曲线Fig.7 The double cumulative curve of annual precipitation-runoff and annual precipitation-sediment discharge in the Pearl River Basin from 2000 to 2019

各流域年降水量-输沙量累积曲线表现不同。西江流域(高要站)，曲线从2007年开始向下偏转，2007—2019年流域的年均输沙量相较2000—2006年下降了48.39%。北江流域(石角站)，曲线未发生明显偏转，表明输沙量主要受降水的影响。东江流域(博罗站)，从2009年开始输沙量下降，曲线向下发生偏转。

在龙滩水库建设前(2000—2006年)，西江流域(高要站)年输沙量与年降水量呈显著正相关(R=0.93，p<0.01)，而NDVI与年降水量和年输沙量均不相关(图8a)，表明输沙量变化主要受到降水量影响。龙滩水库建成后的2 a内(2007—2008年),输沙量变化剧烈，2009年输沙量变化开始趋于平稳，表明水库建设导致的输沙量变化逐渐稳定。2009—2019年西江流域NDVI与年降水量呈正相关(R=0.79，p<0.01)，与年输沙量不相关(R=0.55，p<0.1)，表明植被改善受降水量影响但未起到减沙的作用。

2000—2019年北江流域(石角站)的年降水量与年输沙量呈正相关(R=0.61，p<0.01)，与NDVI不相关(R=0.37，p>0.1)。2008年开始，北江流域NDVI呈明显上升趋势。2000—2008年该流域输沙量与NDVI不相关，与年降水量呈显著正相关(R=0.90，p<0.01)，表明该时期北江流域的输沙量主要受降水影响。2009—2019年NDVI与年输沙量呈不显著负相关(R=-0.52，p<0.1)，与河流含沙量(输沙量与径流量的比值)呈显著负相关(R=-0.84，p<0.01)，说明植被改善是该时期北江含沙量减少的主要原因，但其对输沙量的影响并不显著(图8b)。

2000—2019年东江流域(博罗站)年降水量与年输沙量呈正相关(R=0.56，p<0.05)，与NDVI不相关(R=0.26，p>0.1)。其中，2000—2008年东江流域的年输沙量与NDVI不相关，与年降水量呈显著正相关(R=0.84，p<0.01)，说明该时期东江流域的输沙量受降水影响显著。2009—2019年东江流域的NDVI与降水量呈正相关(R=0.60，p<0.05)，与输沙量呈不显著正相关(R=0.55，p<0.1)，表明降水有利于植被恢复，而植被恢复并未影响输沙量(图8c)。

图8 2000—2019年珠江流域NDVI与降水量和输沙量的关系Fig.8 Relationship between annual NDVI and annual precipitation and sediment discharge in the Pearl River Basin from 2000 to 2019

根据以上分析，2000—2019年间珠江流域年降水量、径流量、NDVI以及输沙量之间存在着相互关系，径流量主要受降水影响，NDVI在某些时期与降水量具有正相关性。在不同时期，西江、北江和东江流域NDVI与输沙量的相关性有差异。西江与东江流域植被改善未起到减沙作用，北江流域植被改善对输沙量的影响也不显著。因此2000—2019年整个珠江流域NDVI变化对输沙量影响不明显。

4.3 1954—2019年人类活动对输沙量的影响

分别利用3个水文站1954—2019年的年径流量和年输沙量的累积值绘制双累积曲线(图9)。

1954—1982年西江(高要站)的径流量与输沙量双累积曲线基本为直线，表明输沙量主要受径流量影响，人类活动对输沙量影响较小(图9a)。1983—1992年曲线向上偏转，表明水沙关系未发生明显变化的情况下，人类活动使输沙量增加。1993—1997年、1998—2006年、2007—2019年3个时期水沙累积曲线均为向下偏转，表明西江输沙量受人类活动影响逐渐减少。20世纪80年代开始，我国实行家庭联产承包责任制，大规模的毁林开荒活动破坏了地表植被，加剧了流域的水土流失。据统计，1984—1988年珠江的水土流失面积达到5.71×104km2，相比1954年，水土流失面积增加了38.9%[33]。根据图9a中方程式计算，1982—1992年西江河流输沙量累积增加205 Mt。1993—1997年与1998—2006年斜率有所减小，与岩滩水库(1992年)和天生桥水库(1997年)开始蓄水拦沙的时间吻合，龙滩水库建成后(2006年)输沙量进一步减少，说明水库拦沙作用显著。总体上，1954—2019年，西江流域河流输沙量随着流域水库建设而显著减少。

北江(石角站)的径流量与输沙量双累积曲线(图9b)整体变化趋势不明显，1954—1980年曲线基本为直线；
1981—1987年曲线向上偏转，这主要归因于北江流域森林砍伐导致的水土流失；
1999—2005年输沙量向下偏转，这与北江飞来峡水电站的建成时间(1999年)相吻合，说明该时期输沙量下降主要受到水库拦截蓄水影响；
虽然斜率在2006—2019年较1999—2005年变大，但总体上输沙量仍呈缓慢下降趋势。北江流域的水库建设较少，且库容不大，对水沙的影响总体不大。

东江(博罗站)的径流量与输沙量双累积曲线(图9c)分为1954—1973年、1974—1991年、1992—2010年3个阶段，斜率逐渐变小。东江流域位于东南沿海，森林和耕地面积相对较少，与西江和北江流域相比，东江流域的森林砍伐对流域输沙影响较小，东江输沙量主要受到流域大型水库建设蓄水拦沙影响。东江流域1954—2019年多年平均输沙量仅为 2.25 Mt，远小于西江多年平均输沙量，仅为北江多年平均输沙量的41%。水库建设对输沙量的影响不尽相同：新丰江水库于1962年建成，但年径流量-输沙量双累积曲线在1962年没有明显转折；
1973年枫树坝水电站建成投入使用后曲线发生偏转。在各级水库建设、航道整治及挖砂等人类活动影响下，东江输沙量呈现逐级递减趋势。

图9 1954—2019年珠江流域的年径流量-输沙量双累积曲线Fig.9 The double cumulative curve of annual runoff and sediment discharge in the Pearl River Basin from 1954 to 2019

通过对1954—2019年珠江流域水沙数据的分析可知，西江(高要站)、北江(石角站)、东江(博罗站)的径流量变化趋势不明显，输沙量整体呈下降趋势，其中西江(高要站)和东江(博罗站)的年输沙量序列分别于2003年和1991年发生突变，突变后多年平均输沙量较突变前分别减少了71.30%和48.88%，北江(石角站)的年输沙量变化趋势未发生显著突变。

2000—2019年珠江流域植被覆盖状况改善的区域占流域总面积的89.27%，其中西江、北江、东江流域植被覆盖状况改善的区域分别达90.97%、90.77%和89.67%。2000—2019年西江和东江流域植被阻沙作用不明显性。北江流域仅在2009—2019年，植被覆盖状况改善对输沙量有一定的减少作用，但对整个珠江流域的输沙量影响有限。因此植被覆盖状况改善对整个珠江流域输沙量的影响不明显。

通过年径流量-输沙量双累积曲线分析可知，1980年以前输沙量与径流量高度相关，人类活动对输沙量影响不明显。1980年以后剧烈的人类活动导致输沙量发生变化，其中森林砍伐等人类活动加剧了流域水土流失，使输沙量增大；
水库建设、植被恢复等人类活动使输沙量减小。2000年以后，龙滩水库建设是西江输沙量下降的主要原因。