1901—2016年黄土高原土壤侵蚀格局演变及其驱动机制

穆兴民，李朋飞，刘斌涛，赵广举，高 鹏，孙文义

(1.西北农林科技大学 黄土高原土壤侵蚀与旱地农业国家重点实验室，陕西 杨凌 712100;2.中国科学院 水利部 水土保持研究所，陕西 杨凌 712100;3.西安科技大学 测绘科学与技术学院，陕西 西安 710054;4.中国科学院 水利部 成都山地灾害与环境研究所，四川 成都 610041）

土壤侵蚀对陆地生态系统安全构成极大威胁[1-2］,是我国头号生态环境问题。黄土高原地形复杂、沟壑纵横、土质松散、降雨集中,是我国乃至世界水土流失危害最严重的地区之一[3-5］。气候变化和人类活动从正反两个方面影响土壤侵蚀的发生和发展。全新世中期以前,黄土高原土壤侵蚀主要受气候变化影响,而近1000 a来,人类活动是黄土高原土壤侵蚀不断加剧之主导因素[6］,人口增长导致农作物种植面积扩大与土地利用强度提高,进而加剧水土流失(面积增加、强度提高)[7-8］；近100 a来,在气候变化和人类活动双重作用下,黄土高原土地利用及地表覆盖显著变化,土壤侵蚀模数剧烈波动[9-11］。探明黄土高原土壤侵蚀模数时空变化特征,对黄土高原土壤侵蚀演变及其驱动机制的认识意义重大。水土保持作为新兴学科,历史观测数据不足制约了对历史时期黄土高原土壤侵蚀时空动态变化及其驱动机制认识的深化。

土壤侵蚀模数获取方法有控制模拟试验[12-14］、野外定位观测[15-16］、遥感影像分析[17-18］、土壤侵蚀模型模拟[19-21］等。控制模拟试验、野外定位观测费时费力且成本较高,适用于小尺度研究；遥感影像分析仅能获取20世纪80年代以后的土壤侵蚀变化特征[22］；土壤侵蚀模型模拟虽因目前对土壤侵蚀发生发展的机理尚不非常清楚而存在诸多问题,但仍不失为土壤侵蚀时空演变评估的有效方法。已有较多学者采用上述方法对黄土高原进行了不同时间尺度的土壤侵蚀调查、监测与模拟,但研究时段多集中于20世纪80年代以后[23-24］,鲜有针对百年尺度土壤侵蚀模数的模拟评估。

本文基于修正通用土壤流失方程RUSLE(Revised Universal Soil Loss Eguation)构建了1901—2016年逐年黄土高原土壤侵蚀模数栅格数据集,并采用人类活动干扰较小的22个黄河支流水文站的输沙量实测数据对其进行验证,基于双累积曲线法分析了不同时期土壤侵蚀模数变化的主要驱动因素,以期解决黄土高原土壤侵蚀模数历史数据缺乏问题,为黄土高原土壤侵蚀演变过程研究及水土保持生态建设规划实施等提供方法与数据支撑。

1 资料与方法

1.1 黄土高原概况

黄土高原地处我国西北部,总面积约64万km2,地势西北高、东南低,大致分为高塬沟壑区、农灌区、沙地沙漠区、丘陵沟壑区、土石山区、河谷平原区等6个土壤侵蚀地貌类型区(见图1)。黄土高原属干旱和半干旱气候区[25］,多年平均降水量150～700 mm[8］,降水量远低于潜在蒸散发能力[26］,土壤含水率低,植被类型从南部的森林过渡到北部的荒漠草原[27］。由于黄土高原地区降雨主要发生在6—9月且具有降雨历时短、强度大的特点[28］,加之地形复杂、沟壑纵横、土质疏松,极易发生严重的水土流失[4,10-11］。20世纪以来,黄土高原人类活动日益加剧,如农业生产、战争、森林砍伐、植被恢复、水土保持措施与水利工程建设、城市化等,对区域土壤侵蚀的发生发展产生重要的影响[5］,加之气候变化导致的极端天气事件,使黄土高原土壤侵蚀表现出显著的时空分异特征。

图1 黄土高原土壤侵蚀地貌类型分区

1.2 RUSLE

RUSLE根据降雨、土壤、地形、植被覆盖、水土保持措施等因子计算土壤侵蚀模数,其概念清晰,所需数据较少且易于获取,计算方便[29］,是应用广泛的土壤侵蚀模型之一。RUSLE形式为

式中:A为土壤侵蚀模数,t/(hm2·a)；R为降雨侵蚀力,MJ·mm/(hm2·h·a)；K为土壤可蚀性因子,t·hm2·h/(hm2·MJ·mm)；LS为地形因子(其中L为坡长因子、S为坡度因子)；C为植被覆盖因子；P为水土保持措施因子。

(1)降雨侵蚀力的确定。已有大量研究提出了我国不同区域的降雨侵蚀力计算方法,本研究采用以月平均雨量估算侵蚀力的简易算法[30］计算黄土高原地区降雨侵蚀力:

其中

式中:F为与年降雨量的季节分布有关的指标,mm；p、pi分别为年降雨量、第i个月份降雨量(降雨量数据源自高分辨率地表气候格网数据集,网址为https://chelsa-climate.org/chelsacruts/),mm。

(2)土壤可蚀性因子的确定。基于EPIC模型估算土壤可蚀性因子,并采用Zhang等[31］提出的方法对其进行修正:

式中:KEPIC为基于EPIC模型估算的土壤可蚀性因子；San、Sil、Cla分别为土壤的沙粒、粉粒、黏粒含量,%；OC为土壤有机碳含量,%；Sn为常数,Sn＝1-San/100。土壤类型和属性数据从中国土壤科学数据库(http://vdb3.soil.csdb.cn/extend/jsp/introduction)获取。

(3)坡长和坡度因子的确定。坡长和坡度因子基于数字高程模型(DEM)进行计算,DEM由分辨率为30 m的全球数字高程模型ASTER GDEM(https://nordpil.com/blog/astergdem/)经重采样得到。计算坡长和坡度因子的传统方法适用于坡度≤18%的区域[32］,而黄土高原大部分区域地表坡度＞18%[33］,Liu等[34］发现坡度为9%～55%的坡面土壤侵蚀强度与坡面坡度的正弦呈线性关系,江忠善等[35］根据我国不同地区土壤侵蚀研究成果提出了幂函数形式的坡度坡长因子计算方法,本研究综合考虑上述研究成果,采用如下公式计算坡度和坡长因子:

其中

式中:θ为坡度,%；r为坡长,m；m为指数。

(4)植被覆盖因子的确定。植被覆盖因子根据土地利用类型和植被特征来确定。对源于全新世土地利用数据集HYDE3.2.1(https://archaeology.datastations.nl/dataset.xhtml？persistentId＝doi:10.17026/dans-25g-gez3)的土地利用数据进行重采样,获得分辨率为1 km的栅格数据[36］,将土地利用类型分为森林、草地、耕地、居民地、水体,把水体、居民地的植被覆盖因子分别设定为0.001、0.1,森林、草地、耕地的植被覆盖因子根据2011年全国土壤侵蚀调查数据[37］与1 km分辨率NDVI之间的关系进行计算:

根据遥感影像可以获取1982年之后的NDVI数据,但1982年以前因遥感影像资料缺乏而难以确定NDVI。笔者分析NDVI与月降水量的关系并构建了二者的线性回归方程(其相关系数为0.77、显著性水平为0.001),据此方程,结合由CHELSA cruts数据集得出的月降水量构建了1901—2016年NDVI数据集,并依此估算历史时期的C值。

(5)水土保持措施因子的确定。水土保持措施因子根据HYDE3.2.1和土地利用数据集CCI＿LC(https://www.esa-landcover-cci.org/)来确定,其中CCI＿LC提供了灌溉农田和雨养农田的空间分布[38］。将HYDE3.2.1数据重采样至300 m,与CCI＿LC数据进行叠加,提取各项水土保持措施和耕地(把耕地分为梯田和坡耕地,其中梯田包括灌溉农田和缓坡雨养农田)。参照已有研究成果[24］,把梯田、居民地的P值分别设定为0.2、0.01,把森林、草地、水体、坡耕地的P值均设定为1。

1.3 RUSLE模拟计算精度验证

区域土壤侵蚀模数是难以测度的。本研究基于RUSLE计算土壤侵蚀模数的空间分辨率为1 km,无法根据小区实测数据进行验证,故采用中小河流实测输沙量对RUSLE计算结果进行验证。黄土高原地区在无坝库拦蓄情况下中小河流泥沙输移比接近于1[8］,在1970年以前水土保持措施尤其淤地坝数量很少,中小河流输沙量接近坡面侵蚀产沙量,因此本研究选取有实测资料的22个黄河支流水文站1919—1969年实测输沙量(换算为输沙模数)对RUSLE模拟计算的对应支流土壤侵蚀模数进行验证,验证结果用二者的相关系数、Nash系数及显著性来反映。

1.4 土壤侵蚀强度分级

根据水利部发布的《土壤侵蚀分类分级标准》(SL 190—2007),依据土壤侵蚀模数将土壤侵蚀强度划分为6级:微度,侵蚀模数＜1000 t/(km2·a)；轻度,侵蚀模数为1000～2500 t/(km2·a)；中度,侵蚀模数为2500～5000 t/(km2·a)；强烈,侵蚀模数为5000～8000 t/(km2·a)；极强烈,侵蚀模数为8000～15000 t/(km2·a)；剧烈,侵蚀模数＞15000 t/(km2·a)。

1.5 驱动因素贡献评估

双累积曲线是有效评估气候和人类活动对河川径流量及输沙量变化影响的有效方法之一[39］。本研究采用土壤侵蚀模数与降水量双累积曲线确定土壤侵蚀模数变化的各驱动因素的贡献(如图2所示),在不受人类活动影响的情况下,某一河流的多年累积土壤侵蚀模数和累积降水量二者之比为常数,双累积曲线为线段OT,TT′为其延长线；当流域土壤侵蚀受人类活动影响剧烈时,双累积曲线变为TB1或TB2(TB1为人类活动加剧侵蚀,TB2为人类活动减缓侵蚀),其与直线OT′的偏差为人类活动对土壤侵蚀模数变化的贡献量。

图2 土壤侵蚀模数与降水量双累积曲线示意

2 结 果

2.1 RUSLE模拟计算精度验证结果

分别计算22个支流水文站实测输沙量与RUSLE模拟计算的土壤侵蚀模数的相关系数和Nash系数,以验证RUSLE模拟计算结果的精度。分析表1可知:1957—1969年选取的16个支流水文站实测输沙量与计算的土壤侵蚀模数相关系数为0.48～0.80,其中有10个支流的相关系数大于0.60、有9个支流的相关性达到极显著水平；除白家川和旧县外,其余支流模拟结果的Nash系数接近1,其中有9个支流的Nash系数大于等于0.4。此外,1919—1953年享堂、静乐、兰村、赵城、头、南河川等水文站共有实测输沙资料22 a,与对应年份RUSLE计算的土壤侵蚀模数相关系数为0.75、相关性达到极显著水平。验证结果说明,RUSLE模拟计算结果总体可信,对部分支流模拟计算的可信度较高。

表1 1957—1969年RUSLE模拟计算精度验证结果

2.2 黄土高原土壤侵蚀模数时空变化情况

黄土高原地区1901—2016年多年平均土壤侵蚀模数为5056.86 t/(km2·a),不同区域土壤侵蚀模数存在明显差异(见图3),其中:黄土高原西北部沙地沙漠区与农灌区土壤侵蚀模数大都小于1000 t/(km2·a)；黄土高原中部的丘陵沟壑区、高塬沟壑区平均侵蚀模数分别达8570.06、5781.82 t/(km2·a),黄土高原土壤侵蚀模数大于8000 t/(km2·a)的极强烈和剧烈侵蚀主要集中在这2个类型区,因此本研究把这2个类型区作为分析的重点区域(图3中红线所围区域)。

图3 黄土高原年平均土壤侵蚀模数空间变化情况

1901—2016年黄土高原及重点区域土壤侵蚀模数年际变化情况见图4。1901—1909年,黄土高原平均土壤侵蚀模数为4160.55 t/(km2·a),维持在相对较低水平；1910—1919年,黄土高原平均土壤侵蚀模数与1901—1909年相比有一定的上升,达到4928.70 t/(km2·a)；1920—1929年平均土壤侵蚀模数上升至5663.00 t/(km2·a)；1930—1979年,黄土高原平均土壤侵蚀模数为6195.69 t/(km2·a),强烈以上等级的土壤侵蚀面积急剧增长,有13 a平均侵蚀模数超过8000 t/(km2·a)；1980—2009年,黄土高原的平均土壤侵蚀模数显著下降,平均侵蚀模数为3401.30 t/(km2·a),比侵蚀严重的1930—1979年平均值下降了45.1%；2010—2016年,黄土高原平均土壤侵蚀模数有所反弹。

图4 黄土高原及其重点区域土壤侵蚀模数年际变化过程

表3 1901—2016年丘陵沟壑区各级土壤侵蚀面积占比变化情况 %

1901—2016年黄土高原及2个重点区域各级土壤侵蚀面积占比变化情况见表2～表4。黄土高原中度及以下土壤侵蚀面积占比较大(为54.56%～81.31%),尤其1930年以前与1980年以后中度及以下土壤侵蚀面积占比较大、强烈及以上土壤侵蚀面积占比较小,而1930—1979年强烈及以上土壤侵蚀面积占比相对较大(为43.80%～45.44%)；丘陵沟壑区强烈及以上土壤侵蚀面积占比1979年以前为56.21%～72.16%、1980年以后下降到36.24%～55.96%,高塬沟壑区强烈及以上土壤侵蚀面积占比1920—1979年为50.24%～53.40%、1980年后为24.38%～36.05%；2010—2016年黄土高原及2个重点区域强烈及以上土壤侵蚀面积占比均有所反弹。

表2 1901—2016年黄土高原各级土壤侵蚀面积占比变化情况 %

表4 1901—2016年高塬沟壑区各级土壤侵蚀面积占比变化情况 %

根据黄土高原气候变化及人类活动情况及其对土壤侵蚀的影响,可将研究时段(1901—2016年)大致划分为自然状态时期(1901—1929年)、灾害破坏时期(1930—1949年)、初步治理时期(1950—1969年)、全面治理时期(1970—1998年)及高质量生态恢复时期(1999—2016年)等5个时期(见图5),各时期黄土高原年均土壤侵蚀模数分别为4943.49、6415.55、6022.78、4575.70、3431.79 t/(km2·a)。从自然状态时期到灾害破坏时期,黄土高原土壤侵蚀模数呈上升态势,在灾害破坏时期以后随着人类活动(各类水土保持措施的实施)土壤侵蚀模数呈降低趋势。

图5 黄土高原各时期土壤侵蚀模数变化趋势

黄土高原土壤侵蚀模数与汛期(6—9月)降雨量关系密切,二者波动变化趋势基本一致,尤其1920—1980年丘陵沟壑区二者的波动一致性更为明显,见图6。极端暴雨往往造成极端侵蚀,把黄土高原平均土壤侵蚀模数超过10000 t/(km2·a)的年份称为极端侵蚀年份,则研究时段(1901—2016年)有4 a为极端侵蚀年份(见图5),分别为1917年、1925年、1949年、2013年,这4个年份土壤侵蚀模数分别为11717、15202、10853、10968 t/(km2·a)。

图6 土壤侵蚀模数与汛期降雨量逐年变化过程

2.3 不同时期土壤侵蚀模数变化情况及驱动机制

黄土高原土壤侵蚀以水蚀为主,其主要影响因素可分为气候(主要是降水)和人类活动影响下的土地利用、地表覆盖及微地形改变两类。采用土壤侵蚀模数和降水量双累积曲线(见图7)分析黄土高原及重点区域土壤侵蚀模数变化的阶段性特征及驱动机制:黄土高原及2个重点区域土壤侵蚀模数均发生了2次突变,其中黄土高原及高塬沟壑区土壤侵蚀模数2次突变年份均为1924年和1981年、丘陵沟壑区土壤侵蚀模数2次突变年份为1931年和1981年。按突变年份可把研究时段(1901—2016年)分为3个阶段,把第一次突变前、后分别称为基准期、变化期(分为第一变化期和第二变化期),黄土高原和高塬沟壑区的基准期为1901—1923年、丘陵沟壑区的基准期为1901—1930年,按照前述方法基于双累积曲线计算降水和人类活动造成变化期土壤侵蚀模数的变化量及其占比,结果见表5。

图7 土壤侵蚀模数与降水量双累积曲线

分析表5可知:第一变化期,黄土高原及2个重点区域土壤侵蚀模数较基准期均有所上升,黄土高原地区土壤侵蚀模数阶段平均值由4367.6 t/(km2·a)上升至6317.9 t/(km2·a),其中人类活动是土壤侵蚀模数上升的主要驱动因素(造成的变化量为1652.7 t/(km2·a),贡献占比为84.7%)；第二变化期,黄土高原及2个重点区域土壤侵蚀模数较基准期均显著降低,黄土高原土壤侵蚀模数阶段平均值降至3475.9 t/(km2·a),相对于基准期、第一变化期的降幅分别为20.4%、45.0%,人类活动是土壤侵蚀模数降低的驱动因素,降水对土壤侵蚀模数降低的贡献占比为负值,尽管该阶段降水量较之前有所增大、理论上应增大土壤侵蚀模数,但大规模的水土保持措施和水利工程建设等人类活动仍然使土壤侵蚀模数大幅度降低。

表5 降水与人类活动对侵蚀模数的影响

3 讨 论

3.1 近百年来黄土高原土壤侵蚀模数演变特征

我国区域土壤侵蚀评估研究多聚焦于20世纪80年代以后,而河流水沙变化研究始于20世纪20年代。本研究通过区域土壤侵蚀模拟计算,首次构建了黄土高原长时段逐年土壤侵蚀模数数据集,并探讨了百年时间尺度的黄土高原土壤侵蚀模数变化情况,以期为黄土高原历史时期侵蚀模数研究提供方法与数据支撑。黄土高原20世纪初土壤侵蚀模数相对较低、30—70年代土壤侵蚀模数明显上升、80年代后呈下降趋势,21世纪初处于历史最低水平,2010年以后土壤侵蚀模数有所反弹。与黄土高原土壤侵蚀模数的年代际变化相对应,黄河输沙量变化趋势为在20世纪30年代之前相对较小、30—60年代持续增大、80年代后急剧减小[40-41］。

3.2 人类活动对土壤侵蚀模数变化的影响

人类活动是黄土高原百年尺度土壤侵蚀模数变化的主导因素,这与以往研究结论基本一致[42-44］,上述第一变化期不合理的人类活动导致黄土高原土壤侵蚀模数上升,第二变化期大规模水土保持措施的实施有效降低了土壤侵蚀模数。

20世纪30—40年代我国处于战争时期,人口迁移导致黄土高原中部人口大量增加[41］、森林和草地被大面积砍伐和垦殖,加之当时人们对生态环境保护认识不足,因而使土壤侵蚀迅速加剧。

新中国成立后的20世纪50—70年代,国家实施“以粮为纲”的方针,黄土高原地区人口快速增长,经济发展仍以种植农业为主,为了增加耕地面积和解决农村薪柴问题,森林被进一步砍伐、荒草地被开垦种植[45］,导致土壤侵蚀进一步加剧。尽管这一时期已开始加强水土保持工作,但其效果并不理想[8］。

20世纪80年代以后国家进一步重视水土保持工作,黄土高原地区推广旱地农业增产技术和以小流域为单元的山水林田路综合治理,在保障农业生产的同时,使土壤侵蚀模数有所降低、入黄泥沙量显著减少。生态系统是人与自然生命共同体之基,是人类经济社会存在和发展之基。随着改革开放的推进,全国社会经济快速发展,黄土高原乃至西北地区水土流失严重、生态环境脆弱问题与经济社会发展不适应的矛盾日益突出,为了再造一个山川秀美的西北地区,从20世纪90年代末大规模退耕还林(草)工程的实施,以及黄土高原化石能源开发利用根本上改变了农村能源结构,黄土高原植被覆盖率明显提高,积极贯彻落实生态文明建设基本国策、科学实施各类水土保持项目使水土保持工作逐步向高质量发展转变,同时城市化吸引农村人口大规模向城市迁移,因而有效缓解了黄土高原地区的人地矛盾[46］,使黄土高原地区土壤侵蚀环境得到系统性改变,土壤侵蚀模数显著降低(尤其丘陵沟壑区和高塬沟壑区降低更为显著)[47］、入黄泥沙量显著减少。深思细想,各项水土保持措施是黄土高原土壤侵蚀模数降低的驱动因素,而退耕还林(草)等正确的水土保持政策是控制黄土高原土壤侵蚀的根本驱动力。

3.3 极端降雨及其侵蚀灾害的积极应对

黄土高原时段平均土壤侵蚀模数在2000—2009年处于近百年最低水平的情况下,于2010—2016年出现小幅反弹,其与2013年7月降雨异常偏多造成极端侵蚀有关。随着社会经济发展,在全球极端暴雨频发的背景下,未来黄土高原水土保持工作须重点考虑强监管、重维护、提标准,以提高防御极端降雨引发极端侵蚀的能力。长期以来,对水土保持措施存在重建轻管的问题,使得各项水土保持措施因缺乏维护维修资金而老化失修,制约了其效应的持续有效发挥,如大量老旧的小型淤地坝几近淤满,淤地坝和梯田的管理维护不到位,存在垮塌致灾风险,这是极端暴雨导致土壤侵蚀模数及河流输沙量突增的主要原因之一,因此应加强“后水土保持”监管工作。经过20余a来退耕还林(草)政策引导和持续治理,黄土高原已基本呈现沟坡林灌成荫、梁峁芳草铺地的景观特点,水土流失得到有效遏制,但如果对刚刚恢复的仍显脆弱的林草“资源”进行不合理开发利用,那么将很快使林草植被遭到破坏、入黄泥沙量“回归”高位,因此必须巩固水土流失防治和退耕还林(草)成果。此外,植被恢复中重林轻草、盲目还“林”导致的人工林退化以及土壤干燥化问题比较突出,亟须在巩固和提升退耕还林(草)成果基础上进一步优化各项水土保持措施配置,提升生态系统保持水土的功能；现有水土保持措施投资标准及建设标准规范多产生于20世纪80年代,据此建设的水土保持措施难以抵御全球气候变化背景下的极端暴雨及极端侵蚀事件,应尽快修订和完善水土保持措施设计、施工、验收等标准规范体系,以满足新时代水土保持高质量发展的需求。

4 结论

基于RUSLE建立1901—2016年黄土高原土壤侵蚀模数逐年序列1 km分辨率栅格数据集,评估百年尺度的黄土高原土壤侵蚀模数时空变化情况,并探讨了不同时期黄土高原土壤侵蚀模数变化的驱动机制,结论如下:

(1)用22个黄河支流水文站在人类活动影响较小时期的实测输沙量对RUSLE模拟计算的土壤侵蚀模数进行验证,结果表明二者具有较好的一致性,采用RUSLE模拟计算黄土高原土壤侵蚀模数的精度较高。

(2)百余年来黄土高原土壤侵蚀模数在1924年左右和1981年左右发生了2次突变,20世纪30—70年代土壤侵蚀面积持续增大、侵蚀模数持续上升,20世纪80年代—21世纪初黄土高原土壤侵蚀模数持续降低,至21世纪初降为百余年来的最低值,2010—2016年土壤侵蚀模数出现小幅反弹。

(3)人类活动是黄土高原土壤侵蚀模数变化的主要驱动因素。人口增长导致森林过度砍伐、耕地扩张及过度放牧是20世纪30—70年代土壤侵蚀加剧的主要原因,20世纪80年代及以后各项水土保持措施建设及退耕还林(草)政策的实施使黄土高原土壤侵蚀模数大幅降低。

(4)丘陵沟壑区和高塬沟壑区是黄土高原水土流失严重区域和治理的重点区域。

(5)极端降雨事件是2010—2016年黄土高原土壤侵蚀模数发生反弹的主要原因,鉴于全球气候变化背景下未来极端暴雨事件或将更加频繁,亟须创新黄土高原土壤侵蚀防治策略,提高防御极端暴雨造成严重水土流失灾害的能力。