长江源区1980—2020年水沙变化规律

姚春艳，刘洪鹄，刘 竞

（1.中国科学院教育部水土保持与生态环境研究中心，陕西 杨凌 712100；2.中国科学院水利部水土保持研究所，陕西 杨凌 712100；3.中国科学院大学，北京 100049；4.长江水利委员会长江科学院，湖北 武汉 430010）

近年来，随着气候变化和人类活动的加剧，极端气候事件频发，使得流域径流、泥沙等水文要素及其输移过程发生改变［1］。国内外许多学者对尼罗河［2］、密西西比河［3］、黄河［4］和长江［5］等重要河流开展了河流水沙情势变化规律及归因分析研究。其中，气候要素中降水对河流水沙的影响更为显著［5-6］，人类活动则是通过水土保持工程等的实施减少流域水土流失，进而影响水沙变化［5-7］。不可忽略的是，流域水土保持工程和水库建设是近年来造成河流输沙量下降的重要原因［1］，因此，开展流域水土流失与水沙规律的综合研究对流域水土保持和泥沙治理工作十分重要。

青藏高原是亚洲众多河流的发源地，被称为“亚洲水塔”［5］，对全球气候变化的响应极为敏感［8］。长江源区地处青藏高原腹地，是我国重要的水源涵养地与生态安全屏障区［8］，区内冻土、沼泽、湖泊和高寒湿地等广泛分布，生态环境极为脆弱［9］。近年来，随着全球气候变暖的持续作用，源区内水土流失量明显增加，径流量和泥沙等要素随之发生变化［10］，进而对下游河道和流域管理产生影响［5］。已有研究表明，源区内径流的季节性变化主要受降水的影响［11］，土地利用/覆盖和冰雪融水对径流变化也有一定调节作用［12］。水土流失和入河泥沙量增长的主要原因是降水的增加和频繁的人类活动［13］，而长时间序列内植被变化则与流域内泥沙变化相关性不大［14］。Teng等［15］和Wang等［16］基于模型模拟结果发现降水对源区水土流失影响较大，且未来水土流失量将会继续增加。由此可见，降水是引发源区土壤侵蚀和水沙变化的主要驱动力［15］，一方面直接影响坡面土壤的侵蚀，进而影响泥沙过程［16-17］，另一方面通过影响径流变化间接影响河流携沙能力［18］。尽管对长江源区径流和泥沙变化及影响因素进行了多项研究［11,14］，但是大/中流域尺度中侵蚀产沙与输移的综合研究仍较少，同时，气候变化是影响源区土壤侵蚀与河流水沙变化的重要因素。因此，有必要开展长江源区侵蚀与水沙变化及其与气候变化耦合的响应关系研究，以准确把握长江上游泥沙来源及流域侵蚀产沙对降雨在长时间尺度下的响应关系，为源区生态工程实施效果的研究与评价提供科学依据。

因此，本文基于长江源区1980—2020 年降水、径流和泥沙等资料，研究长江源区年土壤侵蚀量、径流量和泥沙变化特征，阐明长江源区土壤侵蚀变化规律，河道水沙变化对地表侵蚀的响应，以及长时间尺度下泥沙变化对侵蚀和降水的响应关系，以期为推动小流域科学治理和长江流域高质量发展提供理论依据。

1 数据与方法

1.1 研究区概况

长江源区为直门达水文站控制以上的流域（图1），地理位置位于90°～98°E、31°～36°N，总面积约为15.7×104km2。源区地势西高东低，平均海拔超过4000 m，太阳辐射较强，日温差大但年温差小［16］，多年降水量显著增加（P<0.05），年均约为370 mm，水系包括通天河水系、沱沱河水系，呈扇状分布。主要土壤类型为高山草原土和高山草甸土，其中高山草甸土分布面积最大，约占源区总面积63%，植被类型从东南到西北依次分布有高寒草甸、高寒草原、高寒荒漠，冻土面积约占青藏高原多年冻土面积的20%［18］。人口以藏族居民为主，人口密度不到1人·km-2，牧业是区域经济发展的主要生产方式。

图1 长江源区示意图Fig.1 Basic information of the headwaters of the Yangtze River

1.2 数据收集与来源

长江源区直门达水文站1980—2020 年逐年径流及泥沙数据由青海省水文局提供。从国家气象科学数据中心（http://data.cma.cn）下载长江源区及其周围共12 个气象站1980—2020 年逐日降水资料。长江源区高程数据（空间分辨率为1 km）、长江源区归一化植被指数（Normalized Difference Vegeta⁃tion Index，NDVI）数据（1998—2018年，空间分辨率为1 km）和1980年、1990年、1995年、2000年、2005年、2010年、2015年和2020年8期土地利用数据（图2）来源于资源环境科学数据中心（https://www.resdc.cn）。

图2 长江源区1980—2020年土地利用变化Fig.2 Landuse changes in the headwaters of the Yangtze River during 1980-2020

1.3 研究方法

1.3.1 土壤侵蚀量 利用中国土壤流失方程（CSLE）［19］估算长江源区1980—2020 年的土壤侵蚀量。在ArcGIS 10.2的支持下，对各个因子赋值叠加计算，得到土壤侵蚀量。模型如公式（1）所示：

式中：A是土壤侵蚀模数（t·hm-2·a-1）；R是降雨侵蚀力因子（MJ·mm·hm-2·h-1·a-1）；K是土壤可蚀性（t·hm2·h·hm-2·MJ-1·mm-1）；L是坡长因子，S是坡度因子，二者统称为地形因子，无量纲；B是生物措施因子；E是工程措施因子；T是耕作措施因子，B、E、T因子均无量纲（取值范围为0～1）。

（1）降雨侵蚀力因子（R）

基于章文波等［20］建立的降雨侵蚀力算法模型计算年降雨侵蚀力因子值（图3），

图3 长江源区1980—2020年年降水量(a)和降雨侵蚀力(b)变化Fig.3 Variations of annual average precipitation(a)and rainfall erosivity factor(b)in the headwaters of Yangtze River during 1980-2020

式中：R为半月周期的降雨侵蚀力（MJ·mm·hm-2·h-1·a-1）；k为一个半月周期（每个月的前15 d 为一个半月周期，该月的其余天数为另一个半月周期）；Pj是一个半月周期内超过12 mm 的日降水量（mm）；α和β分别为模型的待定参数，计算方程式如公式（3）和公式（4）：

式中：Pd12代表日雨量超过12 mm 的日降雨量（mm）；Py12代表一年内超过12 mm 的日雨量的和（mm）。

（2）土壤可蚀性因子（K）

土壤可蚀性因子（K）是评价土壤被降雨分离、冲蚀和搬运难易程度的一项指标［21］。源区土壤可蚀性因子栅格数据集是基于1∶1 000 000 中国土壤数据库、传统的“土壤发生分类”系统和Nomo 图法建立的K值与土壤性质之间的Nomo 模型得出。长江源区1000 m 分辨率土壤可蚀性K值来源于国家地球系统科学数据共享服务平台（http://www.geoda⁃ta.cn）。利用ArcGIS 10.2进行重采样和插值得到本研究土壤可蚀性因子分布（图4a）。

图4 长江源区K值(a)及LS值(b)空间分布Fig.4 Spatial distribution of K(a)and LS value(b)

（3）地形因子（LS）

基于DEM 数据计算地形因子［22］，得到LS值分布图（图4b）。方程如下：

式中：LS是地形因子；l是累计坡长（m）；β是坡度（°）；m是坡度变量，当坡度≤1°时，m=0.2；当坡度为1°～3°时，m=0.3；当坡度为3.5°～4.5°时，m=0.4；当坡度为≥5°时，m=0.5。

（4）生物措施因子（B）

生物措施因子（B）主要反映植被覆盖对土壤侵蚀的影响，其赋值范围为0～1，越接近0，表示防治越好，几乎不发生侵蚀；越接近1，表示未采取任何措施。具体包括退耕还林、还草等。本研究基于8 期土地利用数据，结合Liu 等［19］和Fu 等［23］提出的赋值原则（表1），利用ArcGIS 10.2 重分类得到B值的空间分布（图5）。

表1 不同土地利用类型的B值Tab.1 B value of different landuse types

图5 不同年份生物措施因子值空间分布Fig.5 Spatial distribution of biological practice factor at the different years

（5）工程措施因子（E）和耕作措施因子（T）

工程措施因子（E）反映了实施的工程措施对于流域内水土流失的影响。长江源区居民将旱地改为石坎水平梯田，通过查表可得石坎水平梯田E值为0.121。耕作措施（E）通过提高土壤抗蚀性，进而防治土壤侵蚀。长江源区属于北部中高原半干旱喜凉作物一熟区［24］，根据土地利用类型图将研究区耕地赋值为0.121 和0.488。其他区域分配为1，得到工程措施因子和耕作措施因子分布。

由于长江源区土地利用变化较小［25］，因此，各期土地利用可代表相邻5～10 a 变化，基于CSLE 模型，估算得到长江源区1980—2020 年土壤侵蚀模数，根据《土壤侵蚀分类分级标准》（SL190-2007）将长江源区土壤侵蚀模数划分为6个等级（0～5 t·hm-2·a-1为微度侵蚀，5～25 t·hm-2·a-1为轻度侵蚀，25～50 t·hm-2·a-1为中度侵蚀，50～80 t·hm-2·a-1为强烈侵蚀，80～150 t·hm-2·a-1为极强烈侵蚀，>150 t·hm-2·a-1为剧烈侵蚀），由于微度侵蚀被定义为小于土壤容许流失量，因此本文仅统计轻度及以上侵蚀强度等级的面积占比［19］。

1.3.2 偏最小二乘结构方程（PLS-SEM）偏最小二乘结构方程（PLS-SEM）是一种结合因子分析和回归的非参数方法［26］，通过迭代估计最大限度地解释独立潜变量的变化。本研究利用PLS-SEM分析，通过与多个潜在变量的相关性来解释降水、土壤侵蚀与水沙变化耦合模式。本文所提出的结构模型由降水量（P）、径流量（Q）、输沙量（S）和土壤侵蚀量（M）组成，使用SmartPLS 3.0软件绘制4个潜变量的影响路径，通过300 次迭代执行PLS 算法以获得路径系数（β），并使用5000次抽样进行bootstrapping以测试路径系数的显著性。当模型拟合中NFI>0.9，SRMR<0.12 时，表明PLS-SEM 模型具有预测有效性，能较强地解释因子间相互影响的变化。潜变量间路径是否通过检验的主要依据为显著性P值，当P<0.05时路径检验通过。

1.3.3 突变分析方法 利用累积距平法和滑动t检验法［27］判断长江源区1980—2020年径流量、含沙量和输沙量序列的突变情况。

2 结果与分析

2.1 土壤侵蚀的变化特征

图6 展示了长江源区1980—2020 年各侵蚀强度面积占比变化。如图所示，轻度侵蚀面积约占总侵蚀面积的70.27%，强烈及以上强度侵蚀面积约占总侵蚀面积的10.27%。长江源区1980—2020 年土壤侵蚀模数呈显著增加趋势（P<0.05）（图7），多年平均土壤侵蚀模数为4.71 t·hm-2·a-1。最大值为8.16 t·hm-2·a-1，出现在1985年；最小值为2.13 t·hm-2·a-1，出现在2015 年。土壤侵蚀模数变化规律与降雨侵蚀力变化规律一致（图3b）。

图6 长江源区1980—2020年土壤侵蚀强度面积占比变化Fig.6 Changes of soil loss intensity area proportion in the headwaters of Yangtze River during 1980-2020

图7 长江源区1980—2020年土壤侵蚀模数变化Fig.7 Change of annual soil loss modulus in the headwaters of Yangtze River during 1980-2020

图8展示了长江源区1980年、1990年、1995年、2000 年、2005 年、2010 年、2015 年和2020 年的土壤侵蚀强度的空间分布状况。从中发现源区1980—2020 年土壤侵蚀强度整体较低，以轻度侵蚀为主，侵蚀强度较大的地方主要分布在东南部，如玉树县。东南部降水量相对较多，且地形多以高山为主，坡面处极容易发生水力侵蚀，导致水土流失情况相对其他区域更为严重。如图7 和图8 所示，长江源区土壤侵蚀时空变化较为复杂。1980 年源区东南部土壤侵蚀面积较大且侵蚀强度较高，1990—1995 年侵蚀强度与侵蚀面积呈现减小趋势。土壤侵蚀模数降低至3.82 t·hm-2·a-1，强烈及以上侵蚀所占面积明显减少，土壤侵蚀得到一定遏制。然而，2000 年后侵蚀面积和侵蚀强度均较高且土壤侵蚀模数逐渐增加，轻度侵蚀面积下降，中度及以上侵蚀面积占比逐渐增加。2010 年后侵蚀量增加速度有所放缓，其中，2015 年土壤侵蚀面积最小且整体侵蚀强度最低，这与2015 年降雨侵蚀力较小有关，值仅为125.5 MJ·mm·hm-2·h-1·a-1（图3b）。

图8 长江源区1980—2020年土壤侵蚀模数空间分布Fig.8 Spatial distribution of soil loss modulus in the headwaters of Yangtze River during 1980-2020

2.2 年径流量、含沙量和输沙量的变化规律

直门达水文站1980—2020年年径流量、年含沙量和年输沙量的变化趋势如图9 所示，直门达水文站多年径流量显著增加（P<0.05），多年平均径流量为142.03×108m3；年径流量最大值在2009 年，为245.14×108m3，最小值为81.95×108m3，出现在1997年。相比于90 年代，21 世纪以来径流量明显增加，增幅约为79.7%。多年含沙量和输沙量均无显著变化趋势，多年平均含沙量为0.69 kg·m-3，最大值出现在1989年，为1.41 kg·m-3，最小值出现在1994年，为0.23 kg·m-3。多年平均输沙量为1063.41×104t；输沙量与含沙量变化较为一致，最大值均出现在1989年，为2980×104t，最小值均出现在1994年，为203×104t。

图9 长江源区1980—2020年径流量(a)、含沙量(b)和输沙量(c)的变化趋势Fig.9 Variations of runoff(a),suspend sediment concentrate(b),and sediment load(c)in the headwaters of Yangtze River during 1980-2020

2.3 降水、土壤侵蚀和水沙的耦合关系

基于1980—2020 年年降水量（P）、年土壤侵蚀量（M）、年径流量（Q）和年输沙量（S）数据，通过PLS-SEM 分析，得到长江源区1980—2020 年降水、土壤侵蚀与水沙变化耦合关系图（图10）。计算得知，模型NFI=0.982（NFI>0.9），SRMR=0.03（SRMR<0.12），表明模型对潜变量的拟合较好。1980—2020年长江源区各潜变量之间的主要关系如下：（1）以降水为驱动因素，降水对土壤侵蚀和河道径流均具有直接的显著正向效应［β=0.727（P<0.01）和β=0.837（P<0.01）］，并进一步通过径流和侵蚀产沙间接影响输沙量变化。（2）径流和土壤侵蚀对输沙量变化具有直接正向效应［β=0.885（P<0.01）和β=0.062］，径流对输沙量变化的贡献度明显大于土壤侵蚀。通过间接效应分析发现，降水对输沙量具有显著正向影响。降水、径流和土壤侵蚀对输沙量变化的解释度达72.5%，其影响路径主要为降水➝径流➝输沙（表2）。

表2 1980—2020年长江源区降水、土壤侵蚀、径流和输沙量路径效应Tab.2 Route effects of precipitation,soil loss,runoff,and sediment load in the headwaters of Yangtze River during 1980-2020

图10 长江源区1980—2020年降水、土壤侵蚀与水沙变化耦合关系Fig.10 Coupling relationship between precipitation,soil loss,runoff,and sediment load in the headwaters of Yangtze River during 1980-2020

3 讨论

3.1 长江源区土壤侵蚀对泥沙变化的影响

本研究利用CSLE模型估算1980—2020年长江源区多年平均土壤侵蚀模数的结果与Liu等［28］和魏梦美等［29］的分析结果较为一致。但是，显著高于Teng 等［15］估算的2.76 t·hm-2·a-1。产生差异的主要原因是计算相关因子值的数据来源及估算方法不同，特别是R因子和K因子。由于所选气象站点数量和位置的不同，计算得到的降雨侵蚀力具有一定差异，本研究中年均R值（269.25 MJ·mm·hm-2·h-1·a-1）与于文竹［30］研究的降雨侵蚀力范围（25.72～448.35 MJ·mm·hm-2·h-1·a-1）略有不同，但在可接受变化范围内。同时，区域土壤侵蚀量和降雨侵蚀力的空间变化均为由西北向东南增加。另外，本文对K值栅格数据集进行重采样及空间插值后，发现源区侵蚀主要集中在草甸土，与全国水土流失与生态安全考察结果一致，这也使得降雨入渗量降低进而导致径流量增加［31］。

源区多年平均土壤侵蚀模数呈显著增加趋势，然而，多年输沙量和含沙量无明显变化趋势。一方面，可能是源区水力侵蚀主要分布在源区东南部，该区域以高山为主，且LS值较高，进而使泥沙容易在坡面淤积；另一方面，虽然新的土壤侵蚀仍在增加，但源区水土保持工作取得了一定成效，源区输沙量整体处于一个相对平衡的状态［9］。另外，20 世纪80 年代源区人为无序开采沙金导致源区的侵蚀面积和侵蚀强度都较大［32］，整体侵蚀情况较为严重；20 世纪90 年代源区植被生态发展趋势良好，没有出现退化区域［33］，生态环境总体趋势良好，因此，水土流失状况得到明显改善；2000—2004 年，近70%的地区出现植被轻微退化的现象［32］，且长江源区降雨侵蚀力增加，导致影响水土流失的外部动力条件增强［34］，进一步加剧了水土流失状况。但是，得益于近些年三江源自然生态保护规划的实施［35］，2010年后侵蚀增加速度有所放缓（图7）。

尽管近年来源区进行了大量的生态保护与恢复工作，但是对比前期水土流失遥感调查数据发现，土壤侵蚀仍在加剧，对生态环境的威胁依然存在，长江源区的水土保护及生态修复工作仍需继续开展。

3.2 不同时期降水、土壤侵蚀与水沙的响应关系的差异

由1980—2020 年长江源区年均径流量累积距平图可知（图11），2004年是年径流量累积距平值发生显著变化的一年，在2004 年达到最低，随后升高。结合滑动t检验法，进一步确定了2004 年为径流量的突变点。年输沙量累积距平值在1988 年和2006 年发生变化，但两个突变点均未通过检验，年含沙量在1994年发生的突变也未通过检验。因此，根据径流量突变点将水文序列划分为1980—2004年、2005—2020年两个阶段。

图11 长江源区1980—2020年径流量累积距平变化Fig.11 Cumulative anomaly changes of runoff in the headwaters of Yangtze River during 1980-2020

源区直门达水文站输沙量较为稳定，径流量在2004年出现拐点，与刘彦等［11］和张永勇等［36］的研究结果一致。径流量发生突变的原因主要是21 世纪初降水的增加［37］与三江源生态保护项目的实施［34］。因此，将2004 年作为分界点，对1980—2004 年、2005—2020年两个时期进行耦合关系分析可知（图12），降水对径流量变化的解释度上升［β=0.690（P<0.01）和β=0.869（P<0.01）］，结合图10分析，3个时期中模型均显示降水与河道径流呈显著正相关，且降水对径流变化贡献度达70%左右，降水是径流量变化的主导因素，与李凯等［37］的研究结果较为一致。

图12 长江源区1980—2004年和2005—2020年降水、土壤侵蚀与水沙变化耦合关系Fig.12 Coupling relationship between precipitation,soil loss,runoff,and sediment load in the headwaters of Yangtze River during 1980-2004 and 2005-2020

通过对图10 和图12 中3 个时期耦合关系结果对比发现降水对土壤侵蚀的正向影响也略有上升，且均达到显著水平。结合降雨侵蚀力（图3b）与土壤侵蚀变化，认为降水是引起源区侵蚀的主要因素。但是，源区仍面临由于人类活动导致草场退化及土地沙化的严峻形势［38］，未来有可能会加剧源区侵蚀并影响源区径流与泥沙变化。

目前大多数研究认为区域土壤侵蚀对输沙量有较大影响［39］，然而，本文研究结果显示土壤侵蚀与输沙量变化没有显著相关关系，源区土壤侵蚀与含沙量、输沙量变化趋势也明显不同。一方面，水文特性和流域的形态等因素也会引起输沙量变化；另一方面由于源区侵蚀泥沙在坡面淤积且源区内河湖、湿地等存储效应的影响［18］，从而导致输沙量变化较为稳定，所以源区输沙量的变化与侵蚀并没有表现出必然联系［40］。

3.3 土地利用及植被覆盖变化对流域水沙的影响

土地利用/覆被变化对流域水文活动起着重要的调节作用［9］，其对流域产流产沙的影响是流域水文学和生态水文学研究的关键问题之一。

21世纪初期长江源区人类扰动较显著，2010年后较为平稳［41］。在本研究时段内（1980—2020 年）源区草地多年平均占比约在68%～70%，其次是未利用地和水域，近40 a 土地利用格局相对稳定（图2）。基于1998—2009年、2010—2018年土地利用和1 km 空间分辨率的归一化植被指数（NDVI）数据与水沙进行相关性分析发现（表3），土地利用和植被指数变化与水沙变化整体几乎呈现不显著相关关系，与郭帅等［42］和Zhang 等［5］的研究结果相同。同时，Naveed 等［12］和Ji 等［43］认为长江源区径流量变化主要由气候变化引起，与土地利用相比，气候变化的影响要大得多。另外，在本研究中土地利用的变化已囊括于土壤侵蚀的计算中，因此，本文重点探讨了降水对水沙变化的影响，不再单独讨论土地利用和植被覆盖变化对流域侵蚀与水沙变化的影响。随着人类活动的增加，加之气候变化的影响和长江源区各项生态项目的不断推进，未来长江源区土地利用和植被覆盖因子的变化不可忽视，应注重未来源区生态保护及水源涵养［37］。

表3 1998—2009年、2010—2018年土地利用、植被覆盖变化与径流、泥沙相关性分析Tab.3 Correlation analysis of landuse,NDVI,runoff,and sediment load from 1998-2009,2010-2018

4 结论

本文基于中国土壤流失方程和结构方程模型分析长江源区水土流失与水沙变化特征及其与降水的响应关系，得出以下结论：

（1）长江源区1980—2020 年土壤侵蚀强度整体较低，以轻度侵蚀为主，约占总侵蚀面积的70.27%，年均土壤侵蚀模数为4.71 t·hm-2·a-1，呈显著增加趋势（P<0.05），但在2000年左右增加趋势最为明显，主要是由于降水的影响。侵蚀强度分布具有明显空间异质性，强度较大的地方主要分布在源区东南部。

（2）1980—2020年直门达水文站多年径流量呈显著增加趋势，年均径流量为142.03×108m3，输沙量和含沙量均无明显变化趋势，水沙变化具有较好的一致性，径流量在2004 年发生突变，突变后呈现增加趋势。

（3）降水对土壤侵蚀、径流和输沙均具有显著正向效应，且降水是土壤侵蚀发生的主导因素，降水主要通过影响径流对输沙变化产生正向影响。未来应及时关注降水带来的土壤侵蚀影响，加强未来源区生态保护和恢复工作。