周 柽，邵东国，顾文权，姚明磊

（武汉大学水资源与水电工程科学国家重点实验室，武汉 430072）

0 引 言

受地形、冻融、风力、水力、重力等多种侵蚀营力影响，青藏高原土壤侵蚀严重［1］，致使土层变薄、生态失衡、环境退化［1-3］。近年来，青藏高原大面积的退耕或开垦等土地利用变化在改变局部微地形的同时，会对土壤侵蚀动力造成影响，减弱或诱发土壤侵蚀［4］。因此，探究土地利用变化对土壤侵蚀的影响及其对土壤减蚀的贡献，不仅对揭示该地区土地利用对土壤侵蚀的扰动效应具有重要的科学意义，而且对合理开发利用水土资源、防治土壤侵蚀具有实用前景。

青藏高原作为我国重要的生态屏障，其生态环境亟待保护，备受中国政府与学者的关注。国务院于2011年印发了《关于印发青藏高原区域生态建设与环境保护规划（2011-2030年）的通知》，以加强青藏高原的生态建设与环境保护。多年来，我国学者就青藏高原的土壤侵蚀特征［5，6］、强度［7，8］及其影响因素［9，10］开展了大量研究。胡玲等［11］采用野外冲刷试验研究了青藏高原公路边坡在不同植被防护措施下的产流产沙特征，揭示了初始产流时间、流量和植被盖度的关系，边坡产流量、产沙量与植被盖度的关系，以及不同植被措施的边坡防护效果；康琳琦等［12］基于RUSLE 模型分析了青藏高原的土壤侵蚀强度时空变化特征，发现不同时期内土壤侵蚀变化不大，侵蚀强度由南向北逐渐减弱；以及许多学者在局部地区开展了地质地貌［13］、降水［14］、植被［15］与土壤侵蚀的关系研究，发现拉萨河流域坡度为15°～25°的区域土壤侵蚀最为严重；三江源地区在降雨次数和降雨量增多的情况下，降雨侵蚀力在1961-2012年间不断增大；三江源地区的植被覆盖度随物种多样性增大而减小，土壤侵蚀随之增大。现有青藏高原土地利用的研究成果，主要聚焦于估算单一年份不同土地利用类型的土壤侵蚀强度［13］，从空间维度分析不同土地利用的侵蚀特征，以此指导土地利用规划，而时空维度上的分析欠缺量化，即土地利用变化对土壤侵蚀的驱动作用，在青藏高原地区尚缺乏量化研究，需要深入研究同一地块上年际间土地利用变化对土壤侵蚀的驱动作用机制，揭示土地利用变化对当下土壤侵蚀的影响与减蚀贡献度，为未来土地利用发展规划提供科学依据。

本文以青藏高原农耕区为研究对象，基于CSLE 模型，解析了农耕区2015年土壤侵蚀状态，探讨了土壤侵蚀强度变化综合指数及区域减蚀贡献度，量化了不同耕地变化模式对自身土壤侵蚀的影响及其对区域土壤减蚀的贡献，以期为青藏高原各个地区的耕地开发利用和保护、土壤侵蚀防治等提供参考。

1 材料与方法

1.1 研究区概况

青藏高原位于中国西南部地区，涉及六个省级行政区（图一），南邻喜马拉雅山脉南缘，北邻昆仑山、阿尔金山和祁连山北缘，西邻帕米尔高原和喀喇昆仑山脉，东邻秦岭山脉西段和黄土高原，介于北纬26°00'～39°47'，东经73°19'～104°47'之间，面积达257.24 万km2。降水季节性显著，主要分布在5-9月份，区域性显著，自东南向西北逐渐减少；地域辽阔，水土资源丰富，其中水资源总量占全国22.71%，土地资源主要用于牧业，其次用于林业，耕地稀少，仅占0.9%，主要分布在一江两河流域与湟水流域，其他区域分布较为破碎，集约化程度较低。

图1 青藏高原省级行政区划图Fig.1 The provincial administrative divisions of the QTP

1.2 数据来源

本研究选用的数据为：①降雨数据，利用由气象站点日观测数据而生成的年降水量数据集，获取R因子的空间分布；②土壤数据，利用世界土壤数据库（HWSD）中的土壤质地和土壤有机碳，获取K因子的空间分布；③地形数据：利用DEM 数据生成坡度及坡长数据，获取LS因子的空间分布；④土地利用及植被数据：利用NDVI 生成植被覆盖度图，结合土地利用类型图，获取B因子的空间分布；土地利用类型图，结合地形数据，获取T因子的空间分布。基础数据、分辨率及来源如表1所示。

表1 本文使用的基础数据Tab.1 Basic data used in this paper

1.3 研究方法

1.3.1 土壤侵蚀模数计算

Liu 等［16］建立了中国土壤流失方程（CSLE），该方程相比USLE 模型，充分考虑中国坡面的土壤侵蚀实际特征和土壤侵蚀植物覆盖与生物防治措施，且已被众多学者广泛应用于中国地区的土壤侵蚀定量评价［17，18］，因此本文采用中国土壤流失方程CSLE 定量评价青藏高原农耕区的土壤侵蚀状况，其基本结构形式如下：

式中，A为单位公顷的年土壤流失量，即土壤侵蚀模数，t/(hm2·a)；R为降雨侵蚀力因子，MJ·mm/(hm2·h·a)；K为土壤可蚀性因子，t·hm2·h/(hm2·MJ·mm)；L为坡长因子，无量纲；S为坡度因子，无量纲；B为植被覆盖与生物措施因子，无量纲；E为工程措施因子，无量纲；T为耕作措施因子，无量纲。其各项因子的计算方法如下所述。

（1）降雨侵蚀力因子R。R因子的计算方法主要分为经典算法和简易算法，前者计算复杂、数据要求高，难以推广至大规模的土壤侵蚀研究。在简易算法中，雨量模式由于数据的易获取性而备受青睐。由于现有数据的限制，本文采用年雨量估算模型［19］计算R因子，其计算模型如下：

式中：R为降雨侵蚀力因子，MJ·mm/(hm2·h·a)；P为年降水量，mm。

（2）土壤可蚀性因子K。由于侵蚀-生产力计算模型（EPIC模型）所需数据较易获取，本文中K因子采用EPIC 模型进行计算，其模型如下［20］：

式中：Sa为砂粒（0.05～2 mm）百分含量，%；Si为粉砂（0.002～0.05）百分含量，%；Ci为黏粒（＜0.002 mm）百分含量，%；C为有机碳含量，% ；EPIC模型计算K值为美制单位，t·arce·h/(100 arce·ft·tf·in)，需将其再乘以0.131 7 方转化为国际制单位，t·hm2·h/(hm2·MJ·mm)。

（3）坡长因子L和坡度因子S。某一点的坡长在土壤侵蚀计算中是流域顶部至该点的坡长距离，即该点的上坡长，坡长因子L和坡度因子S采用刘宝元算法计算，公式如下：

式中：L为坡长因子；S为坡度因子；λ为坡长值，m；θ为坡度值，°；m为坡长系数。

（4）植被覆盖与生物措施因子B。由于不同的土地利用类型和植被覆盖度会致使植被对土壤侵蚀的抑制作用程度不同，故本文参考已有研究成果［21］，通过土地利用类型和植被覆盖度确定B因子，如表2所示。

表2 B因子赋值Tab.2 The B factor values of different land-use types on the QTP

根据农村土地调查，坡度≤2°的耕地为平耕地，坡度＞2°的耕地为坡耕地。植被覆盖度由NDVI计算得来，见下式［22］：

式中：C1为植被覆盖度；NDVI为归一化植被指数；NDVImax与NDVImin分别为研究区内NDVI的最大值与最小值。

（5）工程措施因子E。E因子依赖于工程措施的收集，但在大尺度研究区上，难以进行收集以致没有基础数据，故本文暂未考虑该因子对土壤侵蚀的影响。

（6）耕作措施因子T。现有研究中［23］，一般是根据不同的地形坡度下，等高耕作抑制土壤侵蚀的效果得到耕作措施因子T，参考相关研究成果［24］，耕地中与该因子的关系见表3。

表3 青藏高原耕地在不同坡度条件下的T因子Tab.3 The T factor values of the different slope of farmland on the QTP

1.3.2 土壤侵蚀强度变化综合指数及减蚀贡献度计算

穆天龙等［25］构建了土壤侵蚀强度综合指数，用以研究不同土地利用类型土壤侵蚀程度，得到了广泛认可和应用［26，27］。但其公式采用分级值计算，将所有处于某一级侵蚀的地块按照同一侵蚀模数处理，与真实计算结果不符，易产生较大的系统误差。因此，在穆天龙公式的基础上，以单元值代替分级值，提出了土壤侵蚀强度变化综合指数和区域土壤侵蚀强度变化综合指数的概念及计算公式，以此评估土地利用变化对自身和农耕区土壤侵蚀的影响。为评估不同土地利用变化对区域水土流失防治的贡献，提出了区域土壤减蚀贡献度的概念及计算公式，以量化土地利用变化对区域土壤减蚀的贡献。

（1）土壤侵蚀强度变化综合指数，指某种土地利用变化类型引起的土壤侵蚀强度变化加权值，反映不同土地利用变化的自身土壤侵蚀的变化程度。其计算式为：

式中：j为土地利用变化类型；i为土地利用变化单元；n为土地利用变化单元数；ΔIj为j类土地利用变化下，土壤侵蚀强度变化综合指数，t/(hm2·a)；ΔAi为第i个单元的土壤侵蚀模数变化，t/(hm2·a)；Sji为第i个单元的j类土地利用变化面积，hm2；Sj为j类土地利用变化总面积，hm2。

其中，ΔAi的计算如下：

通过CSLE 模型可知，在短时间内，土地利用变化实质上是通过引起植被覆盖与生物措施因子B、工程措施因子E和耕作措施因子T的变化，进而引起土壤侵蚀的变化。因此，本文将2015年的降雨侵蚀力因子R代入2000年的CSLE 计算式［式（9）］，获得在2000-2015年之间土地利用类型不变的假定情景下，2015年的土壤侵蚀模数情景估算值，将情景估算值与2015年的实际值作比较，差值即为土地利用变化引起的土壤侵蚀模数变化（式（10））。

式中：A估为土地利用不变情景下，2015年的土壤侵蚀模数情景估算值，t/(hm2·a)；R2015为2015年的R因子；B2000、T2000为2000年的B和T因子；其余因子为常数，同式（1）。

式中：ΔA为土地利用变化引起的土壤侵蚀模数变化，t/(hm2·a)；A2015为2015年的土壤侵蚀模数实际值，t/(hm2·a)。

（2）区域土壤侵蚀强度变化综合指数，指区域内的土地利用变化引起的土壤侵蚀强度变化加权值，反映出不同土地利用变化类型引起的区域土壤侵蚀变化，以及该区域因土地利用变化形成的土壤侵蚀变化程度。其计算式为：

式中：m为土地利用变化类型数；ΔI为区域土壤侵蚀强度变化综合指数，t/(hm2·a)；S为区域总面积，hm2。

（3）区域土壤减蚀贡献度，指某种土地利用变化所引起的土壤侵蚀变化，占区域内各类土地利用变化引起的土壤减蚀总和的比例，能反映不同土地利用变化类型对研究区土壤减蚀的贡献程度。计算式为：

式中：k为土地利用变化类型，k∈[0，j]；Cj为j类土地利用变化类型对区域土壤减蚀的贡献度；ΔIk为k类土地利用变化下，土壤侵蚀强度变化综合指数，其中ΔIk＜0，t/(hm2·a)；Sk为k类土地利用变化总面积，hm2。

2 结果与分析

2.1 土壤侵蚀空间分布

（1）各因子的空间分布。CSLE 各因子的空间分布如图2，年降雨侵蚀力因子R如图2（a）所示，R的均值为1 137.87 MJ·mm/(hm2·h·a)，最小的R因子(R= 41.63MJ·mm/(hm2·h·a))位于农耕区的西北部地区，主要是由于喜马拉雅山和喀喇昆仑山阻挡了暖湿气流，致使全年降雨极少［28］；最大的R因子(R=2 737.75 MJ·mm/(hm2·h·a)）位于农耕区的东南部地区，这是因为雅鲁藏布江下游的峡谷向南展开的喇叭口，便于孟加拉湾暖湿空气进入青藏高原东南部，形成降水丰沛的“舌状地带”［28］。

土壤可蚀性因子K如图2（b）所示，K的均值为0.034 8 t·hm2·h/(hm2·MJ·mm)，可蚀性最低的土壤［K=0.025 51 t·hm2·h/(hm2·MJ·mm)]位于柴达木盆地，富含碳酸钙、黏胶状物质，使得土壤颗粒胶结在一起，不易分离。最容易被侵蚀的土壤(K=0.039 18 t·hm2·h/[hm2·MJ·mm)]位于阿里山地，冷漠土广布，植被稀疏，常见裸地，易被侵蚀。K的分布和Teng 等人［29］、刘斌涛等人［20］的成果相似。

坡度坡长因子LS如图2（c）所示，LS的均值为11.53，LS的最小值(LS= 0.56)位于柴达木盆地，海拔低且地势平坦。LS的最大值(LS= 29.20)位于川西藏东高山深谷及云贵高原，处于中国地势第二级阶梯上，呈波状起伏，山高谷深，落差大。

植被覆盖与生物措施因子B如图2（d）所示，B的最小值(B= 0.122)位于东喜马拉雅南翼的热带雨林中，植被茂密，此外，在川西高山和云南高原里的热带雨林，B因子较小。B的最大值(B= 0.484)位于昆仑山北翼，植被覆盖度低至15.60%～20.73%之间。

耕作措施因子T如图2（e）所示，反映的是人类耕作措施对耕地侵蚀的减少，集中分布在农田，在一定范围内，地势越平坦，耕作措施效果越好，T因子越小。

图2 CSLE各因子空间分布Fig.2 Distribution of the CSLE factors of the QTP

（2）土壤侵蚀空间分布特征。将CSLE 各因子的图层相互乘积，计算得到2015年的青藏高原农耕区的土壤侵蚀模数，根据《土壤侵蚀分类分级标准》（SL190-2007）（表5），得到土壤侵蚀强度等级图（图3）。

图3 2015年土壤侵蚀强度等级空间分布Fig.3 Spatial distribution of soil erosion intensity grades in 2015

结果表明，2015年，青藏高原农耕区的平均土壤侵蚀模数为195.61t/(hm2·a)，潜在土壤流失量达4.465 1 亿t/a。农耕区土壤侵蚀呈明显的空间分异特征，东南部地区侵蚀相对强烈，西北部地区侵蚀相对较轻，这与降雨、地形因子的空间分布规律相似，根据地理探测器的核心思想，即当自变量对于因变量的影响至关重要时，自变量与因变量的空间分布理应存在相似性，因此，在青藏高原农耕区，降雨、地形因素是土壤侵蚀的重要影响因素。

未利用地的土壤侵蚀模数最高（表4），但由于面积有限，造成的土壤侵蚀量并不大。耕地的土壤侵蚀模数相对较高，侵蚀量占总量的99.52%。草地土壤侵蚀模数相对较小，林地的土壤侵蚀模数最低，这与其他地区的相关研究结果相似［4，30］，说明土地利用对土壤侵蚀的影响规律在不同地区具有一定可复用性。

表4 不同土地利用类型的土壤侵蚀情况Tab.4 Soil erosion of different land use types

微度侵蚀的面积最大，占49.04%（表5），主要分布在西北部，强烈侵蚀面积最小，占4.07%，主要分布在中部，剧烈侵蚀分布在东南部，占19.12%，土壤侵蚀强度由西北向东南逐渐增大，这很可能是受降雨和地形因子的影响。

表5 研究区的各土壤侵蚀等级面积及占比Tab.5 Areas and proportions of different soil erosion grades

（3）结果合理性分析。将本文计算得到的2015年土壤侵蚀估算值与国家青藏高原科学数据中心的数据集［31］作比较，不同土壤侵蚀等级面积的重合率为90.91%，不重合面积来源于本文对剧烈等级面积的估算较大，其原因是本文采用年降雨量模型进行计算，未区分降雨和非降雨时段，一定程度上放大了土壤侵蚀强度，剧烈等级面积的差异比例为7.42%，但对剧烈及以下侵蚀等级的估算结果良好，其土壤侵蚀量相对误差为0.48%。说明模型估算结果是合理的。

2.2 不同土地利用变化对土壤侵蚀的影响

基于公式（8）～（10），得到农耕区不同土地利用变化下的土壤侵蚀强度变化综合指数，计算结果如图4 所示。土壤侵蚀强度变化综合指数大于0，表明该类土地利用变化引起增蚀，数值越大，增蚀作用越强；小于0，引起减蚀，数值越小，减蚀作用越强；绝对值越大，土壤侵蚀强度变化越大，说明该土地利用变化对土壤侵蚀的影响越大。

图4 不同土地利用变化的土壤侵蚀强度变化综合指数Fig.4 Comprehensive index of soil erosion intensity change for different land-use changes

就研究区而言，综合指数大小排序为：退耕还林＜退耕还草＜未利用地转耕地＜草地转耕地＜林地转耕地＜耕地转未利用地，反映出退耕还林、退耕还草和未利用地转耕地减弱了土壤侵蚀，减蚀作用依次降低；草地、林地开垦为耕地和耕地退化为未利用地增大了土壤侵蚀，增蚀作用依次增大，该规律在不同地域基本表现一致，说明土地利用是影响土壤侵蚀的关键因素。

综合指数绝对值大小排序为：退耕还林＞耕地转未利用地＞退耕还草＞林地转耕地＞草地转耕地＞未利用地转耕地，其中，退耕还林的综合指数绝对值最大，原因一是土地利用类型之间的自身土壤保持能力差异，林地最高，草地次之，耕地和未利用地分列三四，其中林地和耕地之间的差别最大；另一个原因是降水、地形的区域性显著，退耕还林区有66.35%的面积分布在四川和云南，领先其他变化类型，而这些区域降雨量最大且地形有利于降水的冲刷，因此，当下垫面有了更强的降雨植物截留能力，能更大程度减轻降雨对土地的冲刷作用，从而引起了更大的减蚀。其他土地利用变化的综合指数绝对值排序的原因亦然，主要是源于自身土壤保持能力差异和降水、地形的地区性差异。

与此同时，降水、地形及植被覆盖的地区性差异使得土地利用变化对土壤侵蚀的影响具有一定的区域性特征：

（1）不同土地利用变化在不同区域的表现稍有差异，如在青海和西藏，退耕还草的综合指数绝对值比退耕还林更大。

（2）相同土地利用变化在不同区域的指数不同，如，退耕还林在各省的综合指数都不相等。

（3）甚至局部地区的土壤侵蚀变化方向与研究区不一致，如青海海西蒙古族藏族自治州和西藏阿里地区的草地转耕地区，其综合指数均小于0，减弱了土壤侵蚀，主要原因是原有草地的植被覆盖度仅31.37%和20.48%，草地抗侵蚀能力大大下降，转为耕地后，由于坡度较缓，仅1.2°和2.2°，耕作措施起到了水土保持作用。

2.3 不同土地利用变化对区域土壤侵蚀的影响及减蚀贡献

基于2.2 节的土壤侵蚀强度变化综合指数，按公式（11）计算农耕区不同土地利用变化类型与区域的土壤侵蚀侵蚀强度变化综合指数，计算结果如图5 所示。土壤侵蚀强度变化区域综合指数大于0，表明该类土地利用变化引起区域的增蚀，数值越大，增蚀作用越强；小于0，引起减蚀，数值越小，减蚀作用越强。

图5 不同土地利用变化的区域土壤侵蚀强度变化综合指数Fig.5 Comprehensive index of regional soil erosion intensity change for different land-use changes

从区域角度分析，由研究区的区域土壤侵蚀强度变化综合指数可知，农耕区在2000-2015年间的土地利用变化，使得研究区的土壤侵蚀强度减少245.54 t/(hm2·a)，较土地利用不变时下降55.66%，且各省农耕区的区域土壤侵蚀强度变化综合指数均小于0，说明不论是在整体上还是局部地区，近年来农耕区土地利用规划有助于水土流失的减少和生态环境的保护，其中四川省的区域变化综合指数最小，这是因为退耕还林和退耕还草的土壤侵蚀强度变化综合指数最小，即减蚀效果最佳，且四川省拥有45.74%的退耕还林和29.12%的退耕还草，因此四川省的减蚀幅度最大。

基于区域土壤侵蚀强度变化综合指数，按公式（12）计算得不同土地利用变化的研究区减蚀贡献度如表6所示。从土地利用变化角度分析，对研究区而言，不同土地利用变化对农耕区土壤减蚀的贡献排序为：退耕还草＞退耕还林＞未利用地转耕地＞林地转耕地＞耕地转未利用地＞草地转耕地，该次序主要是受土壤侵蚀强度变化综合指数和变化类型的面积的双重影响，如草地相比于林地，缺少繁茂的树冠、树干、枝叶和发达的植物根系，抗侵蚀能力稍逊一筹，因此在同面积下，退耕还草的减蚀效果弱于退耕还林，因此退耕还草的土壤侵蚀强度变化综合指数绝对值小于退耕还林，但退耕还草的面积是退耕还林的4.88倍，综合计算，退耕还草贡献了比退耕还林更多的土壤减蚀量。

表6 不同土地利用变化的研究区减蚀贡献度Tab.6 Contribution of regional erosion reduction of different land-use changes

由前述可知，不论是对自身还是农耕区而言，退耕还林和退耕还草都起到了较大的减蚀作用，林地、草地开垦为耕地和耕地转未利用地引起了6.53%土壤减蚀的流失，因此，为减少农耕区的水土流失，应侧重于开展退耕还林、退耕还草，重点防止耕地退化，尽量避免将林地、草地开垦为耕地，而应多关注未利用地，科学论证未利用地开发为农用地的可行性。

3 结 论

本文研究了年际间土地利用变化对土壤侵蚀的驱动作用，与单一年份土地利用的土壤侵蚀特征相比［14，17］，揭示了土地利用变化对当下土壤侵蚀的影响，提出了土壤侵蚀强度变化综合指数和区域减蚀贡献度的概念和计算公式，量化分析了2000-2015年间土地利用变化对土壤侵蚀的驱动效果以及对区域土壤减蚀的贡献度。主要结论如下：

（1）2015年，青藏高原农耕区的平均土壤侵蚀模数为195.61 t/(hm2·a)，潜在土壤流失量达4.460 1 亿t/a。农耕区以微度侵蚀为主，受降雨、地形因子的影响，土壤侵蚀呈明显的空间分异特征，由西北向东南逐渐增大。

（2）土地利用变化是引起土壤侵蚀变化的重要因素，退耕还林、退耕还草和未利用地转耕地的土壤侵蚀强度变化综合指数分别为-378.24、-118.54 和-5.44 t/(hm2·a)，减蚀作用依次降低；草地、林地开垦为耕地和耕地退化为未利用地的土壤侵蚀强度变化综合指数分别为16.53、32.34和200.61 t/(hm2·a)，增蚀作用依次增大，效果的差异可能主要来源于土地利用自身的土壤保持能力差异和降水、地形的地区性差异，并且这些差异使得土地利用变化对土壤侵蚀的影响具有一定的区域性特征。

（3）青藏高原农耕区近年来的土地利用规划有助于水土流失的防治，使得青藏高原农耕区的土壤侵蚀强度减少245.54 t/(hm2·a)，较土地利用不变时下降55.66%。不同土地利用变化对农耕区土壤减蚀的贡献排序为：退耕还草＞退耕还林＞未利用地转耕地＞林地转耕地＞耕地转未利用地＞草地转耕地。□