黄土高原典型流域土壤侵蚀对退耕还林土地利用变化的响应

汪滨，张志强

（1. 北京林业大学水土保持学院，水土保持与荒漠化防治教育部重点实验室，北京 100083；2. 山西吉县森林生态系统国家野外科学观测研究站，临汾 042200）

0 引言

由于土壤抗蚀能力较弱以及人类长期的不合理土地开发利用，黄土高原属于世界上土壤侵蚀最严重的地区[1-2]，为了从根本上有效控制土壤侵蚀，中国于1999年起在黄土高原实施了大规模的退耕还林工程。该工程以改变不合理的土地利用、遏制土壤侵蚀、减少入黄泥沙并兼顾生态、经济和社会效益为主要目的，以坡耕地退耕还林还草、荒山造林、封山育林和梯田建设等为主要措施[3-4]，实施后土地利用结构明显优化，取得了显著的综合效益尤其是水土保持效益[4-6]。退耕还林工程历时较长，不同阶段实施了不同的主要措施，直接驱动土地利用持续发生变化，深入研究这一变化对土壤侵蚀的影响，对于促进退耕还林成果巩固和推动退耕还林高质量发展具有重要意义。

土壤侵蚀的发生是降雨、土地利用、地形、土壤等各种因素相互影响与制约的综合结果[7]，其中降雨和土地利用是2个主要影响因素[8-9]。黄土高原有关试验与观测研究表明，降雨量、降雨强度和土地利用都对土壤侵蚀产生了明显影响[10-12]。运用模型情景模拟（如RUSLE（revised universal soil loss equation）[13]、SWAT（soil and water assessment tool）[14]、MMF（Morgan,Morgan and Finney）[15]）和数理统计（如双累积曲线[16]、线性回归[17]、弹性系数[18]）等方法，黄土高原已开展了降雨和土地利用对退耕还林期间土壤侵蚀变化影响程度方面的研究，得出的共同结论是土地利用变化在一定程度上减轻了土壤侵蚀，但在影响程度方面存在着2 种完全不同的结论，一些研究观点认为土地利用变化的影响程度相对较大[16,18-20]，另一些研究的结论则相反[13,17,21]。退耕还林工程的具体实施对象是相关的土地利用类型（以下简称地类），实施的结果是这些地类中部分地区（或地块）发生了地类的转换（即类别变更），如坡耕地通过退耕还林转至林地、草地通过工程整地转至耕地[14,22]，与此同时这些地类内未发生转换的地区（或地块）也在发生变化，如坡耕地实施了梯改、林草地实施了封禁，且从已有的研究来看未发生转换的面积较大[14,22-23]。土地利用是人类对土地自然属性的利用方式和状况[24]，地类转换地区所发生的变化和未转换地区所发生的变化通过改变侵蚀环境，如耕地转草地后植物由作物转变为草丛、林草地的封禁提高了植被覆盖度，都会对土壤侵蚀产生直接影响，故同时研究这2 种变化对土壤侵蚀的影响尤其是剔除降雨变化后的影响，可完整地揭示土地利用变化对土壤侵蚀的影响过程。目前黄土高原已开展了大量有关这方面的研究，其中以对地类或区域土壤侵蚀整体变化影响的研究为主[25-27]，也有一些研究侧重转换对侵蚀的影响[28-30]，但前者没有区分这2 种变化，后者没有关注未转换地区的变化，且均叠加了降雨变化的影响，不能直接反映土地利用变化对土壤侵蚀的影响，关于剔除降雨变化后这2 种变化影响方面的研究少有报道。实施退耕还林工程后，一些地类对土壤侵蚀的影响最大，黄土高原已有的研究表明影响最大的地类主要有旱地和中低覆盖草地[31]、林地和草地[32]等，有的地区建设用地也发挥了重要作用[33]，这些研究主要依据降雨变化下地类整体或其转换地区侵蚀变化的重要性，关于剔除降雨变化后地类转换地区和未转换地区侵蚀变化重要性方面的研究也少有报道。

本研究以山西省吉县清水河流域为对象，应用RUSLE模型分析流域2000—2020 年土壤侵蚀强度的变化特征，通过情景模拟判别土地利用变化和降雨变化对土壤侵蚀的影响程度，将土地利用变化分解为土地利用转换和土地利用改造2 种形式，并在剔除降雨变化影响的基础上，分析退耕还林工程实施过程中这2 种变化形式对土壤侵蚀的影响及其在影响中的相互关系和各自作用，识别影响流域土壤侵蚀变化的主要地类，以完整地反映土地利用变化对土壤侵蚀的影响程度和过程，为黄土高原退耕还林成果巩固及高质量发展有效措施的制定提供科学依据。

1 研究流域及其退耕还林工程实施概况

清水河发源于山西省吉县高天山西麓，经州川河向西南汇入黄河，是黄河中游河口镇至龙门区间左岸的一级支流。吉县水文站控制的清水河流域地理位置为110°36′47″E～110°56′0″E，36°2′18″N～36°16′23″N（图1）。

图1 清水河流域区位及水文站和雨量站分布图Fig.1 Location of the Qingshuihe watershed and distributions of hydrologic and rainfall stations

跨吉县、乡宁两县，总土地面积436 km2，其中吉县境内占94.67%。为黄土残塬沟壑区典型流域，土壤类型以褐土为主。属暖温带半湿润大陆性季风气候，多年平均气温10.2 ℃，多年平均降水量541.5 mm，年内降水不均，根据吉县水文站2000—2020 年的观测结果，汛期5—10 月降雨量占全年降水量的86.58%。

清水河流域从2000 年起属于黄土高原退耕还林工程的重点实施地区，经历了不同的发展时期：2000 年下半年开始实施，2005 年完成退耕和造林任务；从2008 年起进入退耕还林成果巩固期，流域所在的吉县在继续实施荒山造林和封山育林的同时，大力开展坡改梯与基本农田建设；2010 年起吉县发挥地方优势，将发展苹果产业作为成果巩固的重要举措；2014 年下半年起开始实施新一轮工程。为了完整地反映流域退耕还林工程的实施过程，本研究以2000、2005、2010、2014 和2020 年5个标志年份为节点，将流域已进行的退耕还林实施时期（2000—2020 年）分为2000—2005、2005—2010、2010—2014 和2014—2020 年4个时段，并依据工程实施期间主要措施的不同，将其对应的实施阶段分别分为第1 阶段（2000—2005 年）、第2 阶段（2006—2009 年）、第3阶段（2010—2014 年）和第4 阶段（2015—2020 年）。

2 数据与方法

2.1 数据来源与处理

2.1.1 遥感影像与土地利用分类

本研究选取2000 年5 月21 日Landsat 7、2005 年6月12 日Landsat 5、2011 年6 月29 日Landsat 5、2014 年5 月20 日Landsat 8 和2020 年6 月5 日Landsat 8 影像分别作为5个标志年份的代表影像（因2009 年G22 临吉高速公路的开工建设对2010 年沿线地区土地利用的影响较大，故用2011 年影像替代，相应的时段分别调整为2005—2011 和2011—2014 年，对应的实施阶段分别调整为第2 阶段（2006—2010 年）和第3 阶段（2011—2014 年）），这些影像均源自中国科学院计算机网络信息中心地理空间数据云（http://www.gscloud.cn/）。以中国土地利用现状遥感监测数据库中的土地利用分类系统为基础，并参照土地利用现状分类国家标准（GB/T 21010—2017），将土地利用类型划分为耕地、园地、有林地、疏林地、灌木林地、草地、建设用地和水域8 类。在ERDAS 9.2 中经监督分类并借助Google Earth 和91卫图的高空间分辨率影像，进行人工目视判读与修正[34]后得到5 期土地利用空间分布图。利用野外实地采点数据、地块利用历史实地访问调查资料等对5 期分类结果进行精度检验，总体解译精度均达到90%以上。

2.1.2 降雨数据与侵蚀性降雨标准

流域土壤侵蚀主要发生在汛期，流域内及周边共10个站点的汛期逐日降雨数据来源于《中华人民共和国水文年鉴黄河流域水文资料》。依据谢云等[35]基于黄土高原试验小区实测数据的研究结果，将12 mm/d 作为侵蚀性降雨标准。

2.1.3 土壤类型的分布和属性数据

流域土壤类型的分布数据来源于《山西省土壤图集》，属性数据来源于《山西土种志》。

2.1.4 DEM 和NDVI 数据

DEM 数据为空间分辨率30 m 的ASTER GDEM 影像，源自中国科学院计算机网络信息中心地理空间数据云（http://www.gscloud.cn/）。NDVI 数据中2000、2005、2011 和2014 年数据为5—10 月共24 期分辨率500 m 的MODIS NDVI 月合成影像，同样源自地理空间数据云（http://www.gscloud.cn/）。因地理空间数据云无2020 年分辨率500 m 的月合成影像，2020 年数据选取美国航空航天局（NASA）地球观测数据与信息系统（EOSDIS）（https://search.earthdata.nasa.gov/）4 月22 日至10 月30日共12 期分辨率500 m 的MODIS NDVI 16 日合成影像。

2.2 研究方法

2.2.1 RUSLE 模型的应用

黄土高原沟蚀问题突出[36]，其中浅沟侵蚀分布较广[37]，模拟试验结果表明浅沟侵蚀量占坡面侵蚀总量的26.6%～59.7%[38]。流域属于残塬沟壑区，沟蚀严重，参考黄土高原地区的有关研究[39-41]，其侵蚀强度的估算采用增加了浅沟侵蚀影响因子G的改进RUSLE 模型[42]：

式中A为年平均土壤侵蚀模数，t/(hm2·a)；R为降雨侵蚀力因子，MJ·mm/(hm2·h·a）；K为土壤可蚀性因子，t·hm2·h/(hm2·MJ·mm）；L和S分别为坡长和坡度因子，无量纲；G为浅沟侵蚀影响因子，无量纲；C为覆盖管理因子，无量纲；P为水土保持措施因子，无量纲。本研究根据数据的可获取情况，选用了各因子的适当算法。

1）降雨侵蚀力因子R。降雨侵蚀力反映了降雨引起土壤侵蚀的潜在能力[43]，采用章文波等[44]基于日降雨资料提出的算法：

式中Ri为第i个半月的降雨侵蚀力值，MJ·mm/(hm2·h）；Pj为该半月内第j日的侵蚀性降雨量，mm；k为该半月内的侵蚀性降雨日数，d；α和β为参数；Pd和Py分别为研究时段内侵蚀性日降雨量和侵蚀性年降雨量的多年平均值，mm。采用泰森多边形法为各雨量站所控制的区域赋值。

2）土壤可蚀性因子K。土壤可蚀性反映了土壤性质对侵蚀的敏感程度[45]，采用EPIC 模型[46]中的算法以及张科利等[47]基于中国土壤属性提出的因子值修正算法：

式中Sa为砂粒质量分数，%；Si为粉粒质量分数，%；Cl为黏粒质量分数，%；Co为有机碳质量分数，%；Sn=1—Sa/100；K为基于中国土壤属性对KEPIC的修正值。

3）坡长因子L和坡度因子S。坡长和坡度是诱发土壤侵蚀的重要地形因素，L因子和10°以下坡度的S因子直接应用RUSLE 用户手册[48]中的算法，10°以上的S因子采用LIU等[49]基于实测数据提出的算法：

式中λ为水平坡长，m；m为坡长指数；γ为细沟侵蚀和细沟间侵蚀的比值；θ为坡度，（°）。

4）浅沟侵蚀影响因子G。在无植被覆盖的黄土陡坡条件下，浅沟侵蚀发生的临界坡度为15°，采用江忠善等[42]基于实测数据提出的算法：

5）覆盖管理因子C。植被是土壤侵蚀的重要抑制因素，采用蔡崇法等[50]基于实测数据提出的算法：

式中f为植被覆盖度，%。通过采用像元二分模型[51]对NDVI 数据进行反演获取。

6）水土保持措施因子P。水土保持措施也是土壤侵蚀的重要抑制因素。采用谢怡凡等[52]关于耕地、建设用地和水域的赋值方法，耕地（不包括梯田）按坡度分级赋值，坡度≤5°、＞5°～10°、＞10°～15°、＞15°～20°、＞20°～25°和＞25°分别赋值0.100、0.221、0.305、0.575、0.705 和0.800，建设用地和水域按不发生侵蚀赋值0。采用秦伟等[53]关于梯田、园地、林地和草地的赋值方法，梯田赋值0.084，园地参照鱼鳞坑造林地一并赋值0.187，未实施水土保持措施的林地和草地赋值1。

2.2.2 土地利用变化形式的分解及数据提取

在人类活动驱动下土地利用与土地覆被变化关系的分析中，TURNER Ⅱ等[54]分解并定义了不同土地利用下土地覆被的2 种变化形式：一是土地覆被的转换，即一种覆被类型完全转至另一种覆被类型，如森林皆伐改种牧草后转至草地，二是土地覆被的改造，即土地覆被属性的改变，如草地过牧后成为退化草地，并且无论是转换还是改造，均需对土地覆被进行维护，如梯田维修、森林抚育等。本研究应用这一土地覆被变化形式分解原理，将土地利用变化分解为土地利用转换和土地利用改造，两者分别指同一空间位置上土地利用类别发生变更和未发生变更的变化。转换区和改造区的面积数据提取基于土地利用转移矩阵。

2.2.3 土地利用和降雨变化对土壤侵蚀影响程度的判别

为了判别土地利用变化和降雨变化各自对土壤侵蚀的影响程度，本研究采用RUSLE 模型情景模拟方法，提出了依据侵蚀模数升降变化量差值的算法：首先以土地利用与降雨共同变化实际情景下研究时段初的侵蚀模数为参照，分别计算研究时段末较时段初土地利用与降雨共同变化实际情景下侵蚀模数的升降变化量（ΔA）、仅土地利用变化模拟情景下侵蚀模数的升降变化量（ΔAl）和仅降雨变化模拟情景下侵蚀模数的升降变化量（ΔAr），然后计算ΔAl和ΔAr分别与ΔA的差值（D），最后将这一差值在总差值（均取绝对值）中的比重作为影响程度（E，%）的判别依据。仅土地利用变化和仅降雨变化对土壤侵蚀影响程度（分别为El和Er）的计算方法如下：

式中Dl和Dr分别为仅土地利用变化模拟情景和仅降雨变化模拟情景下侵蚀模数升降变化量的差值（t/(hm2·a)），Ab和Ae分别为研究时段初和时段末土地利用与降雨共同变化实际情景下的侵蚀模数，Al和Ar分别为研究时段末仅土地利用变化模拟情景和仅降雨变化模拟情景下的侵蚀模数，t/(hm2·a)。

在RUSLE 模型因子中，K、L、S和G因子短时期内相对稳定，R和C因子逐年发生变化[17]，因水土保持措施贯穿于整个退耕还林工程实施过程之中，P因子也会在短时期内发生变化，参考黄土高原地区的有关研究，将C和P作为反映土地利用变化的因子[19-20]。Ab、Ae、Al和Ar的计算如下：

以上各因子的计算方法与式（1）相同。仅土地利用变化模拟情景下各地类及流域转换区和改造区时段初、末的平均侵蚀模数均分别依据式（20）和式（22）计算。

2.2.4 影响土壤侵蚀变化的主要土地利用类型识别

为了兼顾地类侵蚀的强度和面积，避免地类侵蚀强度变化大而面积变化小或侵蚀强度变化小而面积变化大2 种极端情况并剔除降雨变化的影响，本研究依据仅土地利用变化情景的模拟结果，将侵蚀量变化的大小作为主要影响地类的识别指标。由于各地类转换区与改造区侵蚀量的变化可能有增有减，且这样的增和减均是地类变化影响土壤侵蚀的结果，本研究提出了依据侵蚀量增减变化绝对值来识别主要影响地类的算法，计算式如下：

式中Wi为地类i侵蚀量合计变化量的贡献率，%；Mi为地类i侵蚀量的合计变化量，万t；Mci和Mmi分别为地类i转换区和改造区侵蚀量的增减变化量，万t。

3 结果与分析

3.1 土壤侵蚀变化基本特征

根据水力侵蚀强度分级国家行业标准[55]，将流域土壤侵蚀强度划分为微度、轻度、中度、强烈、极强烈和剧烈6个级别。在ArcGIS 9.3 中计算RUSLE 模型，得到流域2000、2005、2011、2014 和2020 年的土壤侵蚀强度图，对侵蚀强度进行分级，得到土壤侵蚀强度分级空间分布图（图2），叠加分析土壤侵蚀强度图与土地利用图，得到流域各地类的平均土壤侵蚀模数（表1）。

表1 清水河流域2000—2020 年土地利用类型的平均土壤侵蚀模数Table 1 Average soil erosion moduli of land use types in the Qingshuihe watershed from 2000 to 2020

图2 清水河流域 2000—2020 年土壤侵蚀强度分级空间分布Fig.2 Spatial distribution of soil erosion intensity grades in the Qingshuihe watershed from 2000 to 2020

由图2 可知，微度侵蚀主要分布在高天山和人祖山山区、东部和中部大部残塬地区以及清水河河谷地带，与2000 年相比较，2005、2011、2014 和2020 年其分布范围依次明显扩大；轻度、中度和强烈侵蚀主要分布在山地周边、部分残塬和河谷两侧地区，2005、2011 和2014 年其分布范围稳定缩小，2020 年明显缩小；极强烈和剧烈侵蚀主要分布在沟壑地区，2005、2011、2014 和2020 年其分布范围总体上呈依次缩小趋势。由表1 可知，2000 年流域平均侵蚀模数36.21 t/(hm2·a)，侵蚀强度属于中度级别，2005 年升至41.02 t/(hm2·a)，接近强烈级别，2011 年大幅降至24.93 t/(hm2·a)，进入轻度级别，2014 年继续降至23.72 t/(hm2·a)，仍保持轻度级别，2020年大幅降至8.24t/(hm2·a)，进入微度级别。

由表1 还可以看出：2000 年草地平均侵蚀模数高达53.04 t/(hm2·a)，属于强烈级别，而园地只有5.43 t/(hm2·a)，属于微度级别，表明除建设用地和水域外，退耕还林初各地类之间的平均侵蚀模数差异较大；与前一年份相比较，2005 年除草地上升外，其他地类均为下降，而2014 年除有林地下降外，其他地类均为上升，2011 和2020 年各地类则均为下降，2005、2011 和2014 年各地类平均侵蚀模数自高至低的排序依次为草地、疏林地、灌木林地、耕地、有林地和园地，2020 年的排序则依次为草地、疏林地、耕地、灌木林地、园地和有林地，表明通过不同阶段退耕还林措施的实施，各地类的平均侵蚀模数总体处于下降趋势，但地类之间仍存在明显差异。

3.2 土地利用变化对土壤侵蚀的影响程度

在土地利用变化和降雨变化对土壤侵蚀影响程度的情景模拟过程中，将侵蚀性降雨量作为降雨变化的计算依据。流域1995—2020 年侵蚀性降雨量的变化见图3，其中退耕还林实施前（1995—1999 年）与4个阶段（分别为 2000—2005、2006—2010、2011—2014 和2015—2020 年）的年均侵蚀性降雨量分别为269.0、333.6、317.9、434.6 和331.3 mm。TFPW—MK 趋势性检验结果为：显著性检验统计量Z=1.76＞1.64，表明在显著性水平0.10 下1995—2020 年的侵蚀性降雨量呈显著增加趋势（P=0.078＜0.10，年均变化率3.7 mm）。

图3 清水河流域1995—2020 年侵蚀性降雨量的变化Fig.3 Change of annual erosive rainfall in the Qingshuihe watershed from 1995 to 2020

依据式（13）～（23）模拟了流域土地利用和降雨变化各自对土壤侵蚀的影响程度，结果见表2。

表2 清水河流域土地利用和降雨变化对土壤侵蚀影响程度的模拟结果Table 2 Simulations of effect degrees of land use and rainfall change on soil erosion in the Qingshuihe watershed

由表2 可知第1 阶段、第3 阶段末较阶段初仅土地利用变化情景下的流域平均侵蚀模数分别下降5.54 和7.34 t/(hm2·a)，仅降雨变化情景下的流域平均侵蚀模数分别上升14.94 和8.67 t/(hm2·a)，第1 阶段仅土地利用变化和仅降雨变化对土壤侵蚀的影响程度分别为49.46%和50.54%，差距较小，表明土地利用变化和降雨变化在流域土壤侵蚀的变化中起了同等重要的作用，第3 阶段的影响程度分别为61.71%和38.29%，差距明显，表明土地利用变化起了主导作用。

第2 阶段、第4 阶段末较阶段初仅土地利用变化情景和仅降雨变化情景下的流域平均侵蚀模数均为下降，第2 阶段分别下降11.90 和6.04 t/(hm2·a)，第4 阶段分别下降13.92 和12.95 t/(hm2·a)，这2个阶段仅土地利用变化的影响程度分别为70.58%和61.86%，明显大于仅降雨变化的影响程度（分别为29.42%和38.14%），表明土地利用变化在流域土壤侵蚀的变化中都起了主导作用。

综上可知，仅土地利用变化情景下的流域侵蚀模数阶段平均下降9.67 t/(hm2·a)，仅降雨变化情景下侵蚀模数阶段平均上升1.15 t/(hm2·a)，两者对土壤侵蚀的影响程度分别为75.23%和24.77%，前者明显大于后者，表明土地利用变化在流域土壤侵蚀的变化中总体上起了主导作用。

3.3 土地利用变化对土壤侵蚀的影响过程

3.3.1 土地利用转换和改造过程对土壤侵蚀的影响

为了分析流域退耕还林不同阶段土地利用转换和改造过程对土壤侵蚀的影响，依据了土地利用转移矩阵和仅土地利用变化情景下阶段初和阶段末各地类转换区与改造区平均侵蚀模数的模拟结果。第1 阶段退耕还林工程实施的主要措施为坡耕地退耕、荒山造林（包括补植，配套了鱼鳞坑整地等）和封山育林，未退耕地和园地加强了相关水土保持措施，如坡耕地梯改和水保耕作、园地整治等，并且这些措施实施后草地、灌木林地、疏林地和有林地内部及之间植被的恢复与更新演替开始见效，由此驱动的土地利用转换和改造过程影响了各地类转换区和改造区土壤侵蚀强度的变化。耕地转草地、草地转3 类林地是此阶段的主要转换过程，地类类别的变更（改变了转换区的土壤侵蚀环境）及所实施主要相关措施的差异直接影响流域土地利用转换区侵蚀强度的变化。以坡耕地退耕为例，地类类别从耕地主要转至草地，所实施的主要相关措施由水保耕作转为草地植被封禁，由于地类的侵蚀强度本身已经蕴涵了地类类别和实施措施的影响，因而这一过程主要属于侵蚀强度较低地类向侵蚀强度较高地类转换的过程，耕地转换区的平均侵蚀模数由17.03 t/(hm2·a)升至21.48 t/(hm2·a)（图4a）。耕地、灌木林地和草地的改造是此阶段的主要改造过程，所实施的主要相关措施直接影响流域土地利用改造区侵蚀强度的变化。以灌木林地为例，补植和封禁实施后改造区植被的恢复与更新演替开始见效，平均侵蚀模数由16.58 t/(hm2·a)降至6.90 t/(hm2·a)（图4b）。

图4 清水河流域仅土地利用变化情景下土地利用类型转换区和改造区平均土壤侵蚀模数的升降变化Fig.4 Variations of average soil erosion moduli in conversion regions and modification regions of land use types under scenario of only land use change in the Qingshuihe watershed

第2 阶段为退耕还林成果巩固期，在继续实施第1阶段造林和封禁措施的同时，地方政府开展了新果园建设和荒坡垦殖，以提高退耕后农民的生活水平并解决粮食问题，期间草地和3 类林地内部及之间植被的恢复与更新演替效果明显。耕地转园地、疏林地和灌木林地转有林地、草地转耕地和3 类林地是此阶段的主要转换过程，以疏林地为例，在主要转至有林地的过程中，主要实施措施均为补植和封禁，但地类类别发生变更，转换过程主要属于侵蚀强度较高地类向侵蚀强度较低地类的转换，转换区的平均侵蚀模数由阶段初的27.11 t/(hm2·a)降至阶段末的20.41 t/(hm2·a)（图4c）。有林地、灌木林地和草地的改造是此阶段的主要改造过程，以有林地为例，补植和封禁后植被的恢复与更新演替明显见效，改造区的平均侵蚀模数由8.51 t/(hm2·a)降至6.61 t/(hm2·a)（图4d）。

第3 阶段在延续第2 阶段主要措施的同时，坡耕地梯改持续推进，地方政府为了发展苹果产业，进一步加大了果园整治和新果园建设力度，且期间草地和3 类林地内部及之间植被的恢复与更新演替已开始显著见效。此阶段的主要转换过程为耕地转园地、灌木林地转有林地、草地转园地和3 类林地，以灌木林地为例，主要转至有林地，同上阶段的疏林地一样也属于侵蚀强度较高地类向侵蚀强度较低地类的转换，转换区的平均侵蚀模数由阶段初的10.95 t/(hm2·a)降至阶段末的6.33 t/(hm2·a)（图4e）。此阶段的主要改造过程为园地、有林地、灌木林地和草地的改造，以草地为例，通过严格封禁，植被的恢复与更新演替使得其改造区平均侵蚀模数由55.37 t/(hm2·a)降至50.25 t/(hm2·a)（图4f）。

第4 阶段为新一轮退耕还林工程实施期，地方政府围绕生态林业、民生林业两大主体，突出了第1 轮退耕还林成果巩固、植被恢复和生态果园建设等工作，由此果园整治和新果园建设稳步推进，草地和3 类林地内部及之间植被的恢复与更新演替已显著见效。此阶段的转换过程主要为疏林地转有林地、灌木林地转有林地和疏林地、草地转园地和3 类林地，以草地为例，因其自身的侵蚀强度在各地类中最高，故在地类类别变更及所实施主要措施（草地实施的措施主要为封禁）差异的影响下，转至任何其他地类都属于侵蚀强度较高地类向侵蚀强度较低地类的转换，转换区的平均侵蚀模数由阶段初的42.11 t/(hm2·a)降至阶段末的12.61 t/(hm2·a)（图4g）。此阶段的主要改造过程仍为园地、有林地、灌木林地和草地的改造，以园地为例，整治措施的实施直接降低了侵蚀强度，改造区的平均侵蚀模数由4.83 t/(hm2·a)降至3.23 t/(hm2·a)（图4h）。

3.3.2 土地利用转换和改造在影响土壤侵蚀中的相互关系

由于各地类内的侵蚀强度在空间上并非均匀分布，故转换区和改造区的起始侵蚀强度存在差异：如果某地类内侵蚀强度相对较高的部分地区（或地块）转至其他地类，其转换区侵蚀强度变化的起点就相应较高（反之较低），此时剩余地区（或地块）的侵蚀强度也就相对较低，改造区侵蚀强度变化的起点也就相应较低（反之较高），两者之间存在着此消彼长的关系。

为了分析这一关系，对各阶段初仅土地利用变化情景下流域土地利用转换区和改造区的平均侵蚀模数（表3）进行了比较：流域4个阶段转换区的起始侵蚀模数分别为27.07、24.02、21.84 和23.01 t/(hm2·a)，分别较改造区的起始侵蚀模数（37.59、45.21、26.35 和23.90 t/(hm2·a)）低27.99%、46.87%、17.12%和3.72%，改造区的阶段平均起始侵蚀模数为33.63 t/(hm2·a)，较转换区的阶段平均起始侵蚀模数（23.44 t/(hm2·a)）高43.47%。这一结果表明，各阶段土地利用转换主要是从侵蚀强度相对较轻的地区开始的，而土地利用改造则正好相反，改造区水土流失综合治理的难度总体上大于转换区。

表3 清水河流域仅土地利用变化情景下土地利用转换区和改造区各时段初与末的平均土壤侵蚀模数Table 3 Average soil erosion moduli in land use conversion regions and modification regions under scenario of only land use change at the beginning and the end of each period in the Qingshuihe watershed

3.3.3 土地利用转换和改造在影响土壤侵蚀中的各自作用

土地利用转换是在不同地类之间进行的，而改造是在同一地类内进行的，在起始侵蚀强度差异的基础上，转换和改造过程的不同直接导致转换区和改造区之间侵蚀强度的升降变化量出现差异。为了分析这一差异，利用表3 对仅土地利用变化情景下流域土地利用转换区和改造区各阶段末较阶段初平均侵蚀模数的升降变化量进行了比较：流域4个阶段转换区侵蚀模数的下降量分别为10.11、12.05、14.15 和14.73 t/(hm2·a)，分别较改造区的下降量（4.86、11.86、4.21 和13.71 t/(hm2·a)）高108.02%、1.60%、236.10%和7.44%，转换区侵蚀模数的阶段平均下降量为13.18 t/(hm2·a)，较改造区的阶段平均下降量（8.74 t/(hm2·a)）高50.80%。这一结果表明土地利用转换拥有较强的侵蚀模数调节能力，在降低其实施地区的侵蚀强度方面发挥了重要作用。

土地利用转换和改造在土壤侵蚀变化中的作用还与各自面积变化的大小有关。从表3 可以看出，流域4个阶段转换面积均明显小于改造面积，这样的面积差异与上述各阶段末较阶段初侵蚀强度升降变化量差异的结合，直接导致这2 种变化区之间侵蚀量的增减变化量也随之出现差异。为了分析这一差异，利用表3进行了侵蚀量增减变化量的计算，结果显示：流域4个阶段转换区侵蚀量的减少量分别为5.76、10.39、19.40 和13.12 万t，改造区侵蚀量的减少量分别为18.42、41.49、12.58 和47.57 万t，除第3 阶段转换区是改造区的1.54 倍外，其他3个阶段改造区分别是转换区的3.20、3.99 和3.63 倍，4个阶段改造区侵蚀量的合计减少量120.06 万t，是转换区阶段合计减少量（48.67 万t）的2.47 倍，占流域阶段合计减少量（168.73 万t）的71.16%。这一结果表明，由于实施面积较大，土地利用改造拥有较强的侵蚀量调节能力，在减少流域侵蚀总量方面发挥了重要作用。

3.3.4 对流域土壤侵蚀变化影响最大的土地利用类型

依据仅土地利用变化情景下各地类转换区和改造区侵蚀量增减变化的模拟结果，利用式（24）～（25）计算各地类侵蚀量合计变化量的贡献率，结果见表4。

表4 清水河流域仅土地利用变化情景下各土地利用类型土壤侵蚀量合计变化量的贡献率Table 4 Contribution rates of total variation of soil erosion amount of land use types under scenario of only land use change in the Qingshuihe watershed

由表4 可知：草地侵蚀量合计变化量（12.12、47.68、22.68 和 39.01万t）及贡献率（46.37%、91.31%、69.32%和63.71%）在4个阶段均位居各地类之首，4个阶段合计变化量为121.49 万t，对流域侵蚀量阶段合计变化量（172.31 万t）的贡献率为70.51%，在4个阶段草地合计变化量中，改造区的阶段合计变化量（81.90万t）占67.41%。这一结果表明在各地类中草地对流域土壤侵蚀变化的影响最大，并主要源自其改造的影响。

4 结论与建议

2000、2005、2011、2014 和2020 年清水河流域平均侵蚀模数分别为36.21、41.02、24.93、23.72 和8.24 t/(hm2·a)，侵蚀强度从中度降至微度级别，地类之间侵蚀强度差异明显，土地利用和降雨变化对土壤侵蚀的阶段平均影响程度分别为75.23%和24.77%，退耕还林土地利用变化是土壤侵蚀变化的主导影响因素。通过将土地利用变化分解为土地利用转换和改造2 种形式并通过情景模拟剔除降雨变化的影响，得到如下结论：

1）流域土地利用转换区侵蚀强度的变化直接受转换过程中地类类别变更及所实施主要相关措施差异的影响，而改造区侵蚀强度的变化直接受改造过程中所实施主要相关措施的影响。

2）流域土地利用改造区的阶段平均起始侵蚀模数为33.63 t/(hm2·a)，较转换区的阶段平均起始侵蚀模数（23.44 t/(hm2·a)）高43.47%，其水土流失综合治理的难度总体上大于转换区。

3）流域土地利用转换区的侵蚀模数阶段平均下降13.18 t/(hm2·a)，较改造区的阶段平均下降量（8.74 t/(hm2·a)）高50.80%，土地利用转换在降低其实施地区侵蚀强度方面发挥了重要作用，而改造区的侵蚀量阶段合计减少120.06 万t，占流域阶段合计减少量（168.73 万t）的71.16%，土地利用改造因实施面积较大在减少流域侵蚀总量方面发挥了重要作用。

4）草地转换区和改造区侵蚀量的阶段合计变化量为121.49 万t，占流域阶段合计变化量（172.31 万t）的70.51%，在草地阶段合计变化量中，改造区的阶段合计变化量（81.90 万t）占67.41%，草地是对流域土壤侵蚀变化影响最大的地类，并主要源自其改造的影响。

流域总体侵蚀强度已降至微度级别，但草地和部分林地的侵蚀仍较严重，建议后续退耕还林工程的实施及成果维护应将提高林草生态系统的质量作为重点，尤其应注重草地改造措施的实施，以期获得更大的综合效益。本研究中土壤侵蚀强度的估算主要基于RUSLE 模型，降雨变化的情景模拟主要针对降雨侵蚀力R因子，没有剔除降雨变化对覆盖管理C因子的影响，建议深入开展这方面的相关研究。

致谢：非常感谢北卡罗来纳大学教堂山分校宋从和教授与自然资源部国土卫星遥感应用中心李宪文研究员的指导以及吉县林业局在实地考察中所提供的帮助。