土壤侵蚀不仅破坏土地资源，引起生态环境恶化，带来淤积、干旱、洪涝等自然灾害，更会直接影响区域生态环境和社会环境的可持续发展，给民众生存与经济发展带来影响^[1]；因此，研究土壤侵蚀的影响因素对区域水土保持规划的制订和施行是十分重要的。区域尺度的土壤侵蚀和降雨过程、土地利用、土壤、地形和地貌等地表状况相关^[2-3]；但是对于区域尺度而言，降雨、土壤类型和地貌在较短的时间尺度内不会发生显著的变化，这就意味着人为活动、土地利用方式的改变是区域土壤侵蚀变化的主要影响因素^[4]；因此，研究土地利用/覆被变化对土壤侵蚀的影响对探索区域内土壤侵蚀规律和开展水土保持等方面工作具有重要的现实意义。

已经有许多研究表明区域土壤侵蚀变化和土地利用/覆被的变化高度相关。例如，郝芳华等^[5]在黄河流域使用SWAT模型发现林地的减少会导致流域产沙增加，耕地的增加则导致产沙增加。刘世梁等^[6]则通过NDVI发现植被是影响区域内土壤侵蚀的重要因素之一。Ciampalini等^[7]在研究Aksum区域的土壤侵蚀对多年农业水土保持措施变化的响应时发现大量的耕地转化成牧场是土壤侵蚀量显著增加的主要原因。徐晓雅等^[8]在澜沧县使用通用水土流失方程发现草地和耕地转变为林地和桉树人工林导致植被覆盖恢复，土壤侵蚀减缓。但是，目前对于中国北方土石山区土地利用/覆被变化对土壤侵蚀的影响的研究相对较少，而对于北方土石山区，土地利用方式改变也是土壤侵蚀的重要影响因素。

本研究以北方土石山区典型区域，北京市境内潮白河水系，包括潮白河及其支流流经的区域作为研究区。潮白河是密云水库的主要入库河流，对密云水库的水质有重要影响，而密云水库是北京市最大的地表饮用水水源地，是首都的生态屏障，其生态环境状况对北京的和谐发展有重要意义^[9]；因此，笔者基于RUSLE，借鉴径流小区的实验思路，利用研究区1980—2010年的多年年均降雨量数据作为多期降雨数据，假设在研究区的时空尺度内土壤类型和地貌条件保持不变的情况下，得到不同时期、不同土地利用/覆被下北京市境内潮白河流域这一区域尺度的土壤侵蚀量；分析25年间区域内土地利用方式的转变和植被覆盖度的变化及特征，模拟区域内不同时期不同地类土壤侵蚀量的变化规律，着重探讨土地利用/覆被变化对区域内土壤侵蚀的影响；期望为实现北京市山区土地资源的优化配置和有效控制土壤侵蚀提供科学指导的决策性依据。

1 研究区概况

潮白河流域位于华北平原北部，地处E 115°34′59.1″~116°57′52.7″，N 39°47′2.6″~41°22′57.9″。北京市境内的潮白河及其支流组成潮白河水系，位于E 116°56′0.46″~116°9′10.28″，N 39°48′13.10″~40°39′40.42″，水系流经延庆区、密云区、怀柔区、顺义区、通州区5大行政区，流域面积为5 688 km²。流域的地形主要为山地和丘陵。流域气候属温带半湿润季风气候，多年平均降水量为488.9 mm，流域内主要降水集中在6—8月，降水量约为300~400 mm，占全年降水量的65%~75%。地带性土壤为褐土，并有山地草甸土、潮土及棕壤等^[10]。植被以暖温带落叶阔叶林为主，树种类型丰富，如华山松(Pinus armandi)、油松(Pinus tabulaeformis)、白皮松(Pinus bungeana)、麻栎(Quercus acutissima)，栓皮栎(Quercus variabilis)和侧柏(Platycladus orientalis)等。农田利用类型多为旱地，耕作制度一般为两年三熟或者一年一熟，农业结构以粮食种植为主，主要的粮食作物有小麦(Triticum aestivum)、玉米(Zea mays)、大豆(Glycine max)等，经济类作物有芝麻(Sesamum indicum)、花生(Arachis hypogaea)、烟草(Nicotiana tabacum)等^[11]。

2 数据来源与方法 2.1 数据来源

本研究所使用的数据(表 1)包括遥感数据、DEM、土壤物理属性以及气象数据。遥感数据下载于地理空间数据云(http://www.gscloud.cn/)(1990、1995、2000、2005、2010、2015年)的Landsat数据，其空间分辨为30 m，据此得到相应年份的植被覆盖度以及土地利用方式。DEM数据来自于地理空间数据云(http://www.gscloud.cn/)的ASTER GDEM V2数据，空间分辨率为30 m。土壤物理属性数据下载于寒区旱区科学数据中心(http://westdc.westgis.ac.cn)，来源于联合国粮农组织(FAO)和维也纳国际应用系统研究所(IIASA)所构建的世界土壤数据库(Harmonized World Soil Database version 1.1)(HWSD)，中国境内数据源为第2次全国土地调查南京土壤所所提供的1:100万土壤数据^[12]。研究区内及附近的气象站点数量偏少，且部分气象站点存在每日、每月降水量的数据记录缺失的情况；因此，本研究采用的气象数据来自国家气象科学数据共享服务平台(http://data.cma.cn/)中国地面累年值年值数据集(1981—2010年)，选取在研究区内及附近的顺义、怀柔、延庆、通州、密云、密云上甸子等11个气象站点，数据类型为累年年平均降水量。

2.2 修正的通用土壤流失方程

本研究采用的是修正的土壤流失方程RUSLE：A=RKLSCP。式中A为年土壤侵蚀模数，t/(hm²·a)，R、K、LS、C和P的计算如下。

1) 降雨侵蚀力因子R。笔者采用徐丽等^[13]提出的北京市年降雨侵蚀力简易估算公式计算：

$ R_{年}=0.440 P_{年}^{1.463}。$

(1)

式中：R_年为年降雨侵蚀力，MJ·mm/(hm²·h·a)；P_年为年降雨量，mm。根据式(1)计算得到北京市潮白河流域内及周边水文站的年降雨侵蚀力。因降雨受地形影响较大，所以结合研究区域的数字高程模型，利用Jin Qiutong等^[14]提出的方法对研究区内及周边的降雨侵蚀力进行多元线性克里格插值，得到1990—2015年北京市潮白河流域年平均降雨侵蚀力因子栅格图。

2) 土壤可蚀性因子K。由于研究区范围较大，且考虑到数据的可获得性、计算的方便性以及前人的研究成果，采用Williams^[15]在EPIC模型中提出的公式，并对其修正以符合我国的土壤特性^[16]。

$ \begin{aligned} K_{\mathrm{EPIC}} &=\left\{0.2+0.3 \exp \left[-0.025 \;6 S_{\mathrm{d}}(1-\right.\right.\\ & S_{\mathrm{i}} / 100 ) ] \}\left[S_{\mathrm{i}} /\left(S_{\mathrm{i}}+C_{\mathrm{i}}\right)\right]^{0.3}\cdot \end{aligned} \\ \begin{array}{r}{\{1.0-0.25 C /[C+\exp (3.72-2.95 C)]\}} \cdot \\ {\left\{1.0-0.7 S_{\mathrm{n}} /\left[S_{\mathrm{n}}+\exp \left(-5.51+22.9 S_{\mathrm{n}}\right)\right]\right\}};\end{array} $

(2)

$ S_{\mathrm{n}}=1-S_{\mathrm{d}} / 100; $

(3)

$ K=-0.013 \;83+0.515\;75 K_{\mathrm{EPIC}}。$

(4)

式中：S_d为砂粒质量分数，%；S_i为粉粒质量分数，%；C_i为黏粒质量分数，%；C为有机质质量分数，%；S_n为粉粒和砂粒质量分数，%。

3) 坡长坡度因子LS。本研究采用的DEM分辨率为30 m，由其直接提取的坡度坡长不能准确反应实际的坡度坡长，因此采用汤国安等^[17]提出的模型计算坡度误差，并据此对GIS软件计算得出的坡度进行校正。

在此基础上，由于国内坡度小区观测资料的缺乏，根据国内外的相关研究，采用McCool^[18]等建立的公式计算缓坡坡度因子S(θ＜10°)，Liu Baoyuan等^[19]建立的公式计算陡坡坡度因子S(θ≥10°)，见式(5)。为了使计算的坡长更加接近实际坡长，坡长因子根据符素华等^[20]采用的由Foster等^[21]提出的分段坡坡长因子公式计算，见式(6)。坡长指数则根据刘宝元等^[22]的研究计算，见式(7)。

$ \left\{\begin{array}{ll}{S=10.8 \sin \theta+0.03, } & {\theta<5^{\circ}} ；\\ {S=16.8 \sin \theta-0.50, } & {5^{\circ} \leqslant \theta<10^{\circ}} ；\\ {S=20.91 \sin \theta-0.96, } & {\theta \geqslant 10^{\circ} 。}\end{array}\right. $

(5)

$ L_{i}=\left(\lambda_{\text { out }}^{m+1}-\lambda_{\text { in }}^{m+1}\right) /\left[\left(\lambda_{\text { out }}-\lambda_{\text { in }}\right) 22.13^{m}\right]。$

(6)

$ \left\{\begin{array}{ll}{m=0.2, } & {\theta<0.5^{\circ}}；\\ {m=0.3, } & {0.5^{\circ} \leqslant \theta<1.5^{\circ}}；\\ {m=0.4, } & {1.5^{\circ} \leqslant \theta<3^{\circ}} ；\\ {m=0.5, } & {\theta \geqslant 3^{\circ} 。}\end{array}\right. $

(7)

式中：L_i为第i个栅格的坡长因子；m为坡长指数；λ_out、λ_in分别为栅格出口及入口的坡长，m。

4) 植被覆盖和管理因子C。植被覆盖与管理因子是植被的覆盖程度对水土流失的抑制程度，起着保持水土的作用，其值在0~1之间。考虑到研究区的雨季集中在6—8月，以及影像的质量状况，选择1990、1995、2000、2005、2010和2015年6月的Landsat影像，计算NDVI(normalized difference vegetation index)^[23]。

根据研究区内的实际植被覆盖状况，采用像元二分的方法，计算植被覆盖度：

$ c=\frac{\mathrm{NDVI}-\mathrm{NDVI}_{\mathrm{min}}}{\mathrm{NDVI}_{\max }-\mathrm{NDVI}_{\min }}。$

(8)

采用蔡崇法等^[24]建立的C因子值与植被覆盖度之间关系方程式计算C值：

$ \left\{\begin{array}{ll}{C=1, } & {c=0} ；\\ {C=0.650 \;8-0.343 \;6 \mathrm{lgc}, } & {0<c \leqslant 78.3 \%} ；\\ {C=0, } & {c>78.3 \% 。}\end{array}\right. $

(9)

式中c为植被覆盖度，%。

5) 水土保持措施因子P。水土保持措施因子P是采取专门措施后的土壤流失量之比，取值范围在0~1之间，0值代表不会发生土壤侵蚀，1值代表未采取任何水土保持措施^[25]。本研究因尺度较大，故结合实地考察资料给不同土地利用类型赋值。将土地利用类型重分类为耕地、林地、草地、水域、城镇居民用地5种，并参考国内学者的研究^[26-28]，将它们的P因子分别赋值为0.400、1.000、1.000、0和0.001。

2.3 土壤侵蚀强度指数

土壤侵蚀强度指数考虑了土地利用类型所占的面积，可更有效地对各地表覆被类型土壤侵蚀强度的贡献情况进行把握^[29-30]；因此本研究将土壤侵蚀强度指数作为土地利用/地表覆被变化对土壤侵蚀影响的量化指标。其计算公式如下：

$ E_{j}=100 \sum\limits_{i=1}^{n}\left(C_{i} A_{i} / S_{j}\right)。$

(10)

式中：E_j为第j单元的土壤侵蚀强度指数；C_i为第j单元第i类型土壤侵蚀强度分级值；A_i为第j单元第i类型土壤侵蚀所占面积，km²；S_j为第j单元所占的土地面积，km²；n为第j单元土壤侵蚀的类型总数。按照水利部颁布的SL190—2007《土壤侵蚀分类分级标准》^[31]，土壤侵蚀强度被划分为微度、轻度、中度、强烈、极强烈和剧烈侵蚀，据此，本研究将各等级土壤侵蚀强度分级赋值为1、2、3、4、5和6。

3 结果与分析 3.1 土地利用变化特征

本研究基于多期Landsat数据，采用目视解译方法对潮白河流域1990—2015年土地利用状况进行提取，验证结果显示解译精度达到90%以上。解译结果显示，1990—2015年研究区土地利用类型以林地、耕地、草地为主(图 1)，三者占研究区总面积的90%以上，其中林地是最主要的土地覆被类型，占总面积的50%以上，主要分布在潮白河流域的北部，这25年间呈先增加后减少再保持不变的变化趋势，尤其是1994年之后的防护林建设和退耕还林等一系列政策的执行，致使初期潮白河林地面积迅速增加。

综合图 1对多期土地利用数据进行叠加分析，统计1990—2015年间不同土地利用类型的面积，并据此获得25年间研究区内各土地利用类型间的转移矩阵(表 2)。耕地在1990—2015年间总体减少14.10%(179.79 km²)，并且流向其他土地利用类型，1995年相对于1990年减少28.20%(361.24 km²)，但随后到2000年期间却增加了26.47%(243.47 km²)，在之后的2000—2005年期间又开始少量减少3.90%(44.82 km²)，2005—2015年期间少量减少，可以看做基本保持不变。林地则在1990—1995年显著增加13.79%，之后逐渐减少，直至2005年后基本保持不变。草地在1990—1995年显著减少27.20%(203.56 km²)，之后缓慢增加，至2015年总体减少4.01%(29.78 km²)。水域面积则增加30.42%(59.25 km²)，其面积的增加主要是在1990—1995年期间，1995年之后基本保持不变。城镇居民用地的变化最为显著，1990—2015年期间，其面积持续增加，直至2015年总共增加110.16%(171.06 km²)。

整体看来，耕地转向林地的面积最多，为242.70 km²，其次转化为城镇，225.47 km²；林地和草地的相互转化最为频繁，二者其次转化为耕地；城镇面积的增加主要是由耕地转入。

3.2 植被覆盖度时空演化特征

笔者参照水利部颁布的SL190—2007《土壤侵蚀分类分级标准》，按照≤30%，>30%~45%，>45%~60%，>60%~75%，＞75%的标准，将1990—2015年潮白河流域的植被覆盖度分为低、中低、中、中高、高度植被覆盖5个等级(图 2)。结果显示，1990—2015年研究区在不同年份的雨季内，植被覆盖度呈先下降后上升的趋势，主要是因为在防护林工程建设初期，以破坏为主，并且树种单一，郁闭度不高，但是在政策的长期实施中，植被在不断恢复，生态环境不断改善，2015年的植被覆盖度达到最佳。

笔者将植被盖度和土地利用数据进行叠加，分析不同地类的植被盖度变化特征(表 3)。1990年，耕地以中低和中等覆盖度为主，林地和草地则是以中高和高植被覆盖度为主；1995年，耕地的植被覆盖度等级降低为以低和中低覆盖度为主，林地和草地仍以中高和高植被覆盖度为主，但是在这2种植被覆盖等级下的林地面积有所增加，草地则有所减少；2000年，耕地的植被覆盖状况恢复至以中低和中等覆盖度为主，林地和草地仍然保持以中高和高植被覆盖度为主。

2005年，这三者的植被覆盖状况均发生退化，耕地的植被覆盖度等级虽然仍以中低和中为主，但是低覆盖耕地面积明显增加，而林地和草地覆盖度等级下降为以中和中高为主。2010年，低覆盖耕地面积进一步增加，但是林地和草地的植被覆盖状况得到小幅恢复，二者在中等覆盖状况的面积减少，高等覆盖状况的面积增加。2015年，这三者均恢复到以中高和高覆盖度为主。

将耕地、林地和草地三者在25年间的变化规律结合土地利用类型的变化特征(表 2)来看，北京市潮白河流域的土地利用/覆被格局持续承受着快速城市化需求及社会经济发展需要与生态环境保护及建设之间彼此矛盾的巨大压力。

3.3 土壤侵蚀强度估算模型的参数因子

根据前文所述公式计算得到因子结果如图 3。研究区降雨侵蚀力因子和地形因子的分布情况有相似之处，东部和中部较高。位于研究区西北部的延庆、怀柔山区，受到海拔和植被覆盖程度的影响，降雨侵蚀力因子值较低，而研究区的东南部为平原地区，其地形因子值较低。研究区的土壤可蚀性因子的空间分布呈片状和点状，由于研究区的主要土壤类型为褐土；因此研究区大部分的土壤可蚀性因子值相差不大，但西部和北部的部分地区值较高，南部的值则相对较低。

3.4 土壤侵蚀强度的时空变化特征

笔者采用RUSLE估算空间尺度的潮白河流域土壤侵蚀强度(图 4)，并分析不同等级土壤侵蚀强度的面积变化(表 4)。结果显示，1990—2015年北京市潮白河流域无剧烈侵蚀，流域整体处于微度至强烈侵蚀，25年间土壤侵蚀强度呈先加重后恢复的趋势，2005年为侵蚀最为严重的年份，之后直至2015年侵蚀得到控制。受到植被覆盖度变化的影响，土壤侵蚀强度的变化趋势与植被覆盖等级先下降后上升的变化趋势密切相关。土壤侵蚀强度的空间分布则与图 3的降雨、土壤和地形条件密切相关，整体呈现南部平原较低，西、北部山区及东部海拔较高地区的土壤侵蚀强度较大。

1990—2015年研究区的侵蚀强度以微度和轻度侵蚀为主，几乎无极强烈及以上侵蚀。1990年微度侵蚀面积占研究区土壤侵蚀面积的67.96%(3 667.13 km²)，轻度侵蚀面积占30.26%(1 632.73 km²)，中度侵蚀和强烈侵蚀面积仅占研究区土壤侵蚀面积的1.78%(96.40 km²)。1995年与1990年的侵蚀状况相差不大。2000年微度侵蚀面积少幅减少，占研究区总面积的59.75%(3 224.27 km²)，轻度侵蚀面积所占比例少幅增长，为37.39%(1 648.27 km²)，中度侵蚀和强烈侵蚀面积所占比例增长至2.86%(154.25 km²)。2005年是土壤侵蚀强度最严重的一年，轻度侵蚀成为主要侵蚀强度，其面积占研究区土壤侵蚀面积的58.39% (3 151.06 km²)，微度侵蚀面积所占比例减少至34.17%(1 843.82 km²)，中度侵蚀和强烈侵蚀所占比例也增至7.44%(401.51 km²)。2010年则恢复到微度土壤侵蚀为主，其面积所占比例为65.06%(3 510.90 km²)，轻度侵蚀面积所占比例下降到28.74%(1 550.67 km²)，中度和强烈侵蚀面积所占比例也降至6.21%(334.83 km²)。2015年的侵蚀强度得到控制，微度侵蚀面积为占研究区土壤侵蚀面积的77.72%(4 193.85 km²)，轻度侵蚀面积占20.91%(1 128.13 km²)，中度和强烈侵蚀面积占比仅1.38%(74.38 km²)。

3.5 土地利用变化对土壤侵蚀强度的影响

笔者估算了不同用地类型的土壤侵蚀强度指数(表 5)。结果显示：耕地的土壤侵蚀强度指数小于林地和草地。对于林地和草地的土壤侵蚀强度指数变化来说，在2005年达到最大分别为201.57和190.10，在2005年之前，林地和草地的土壤侵蚀强度指数变化相对较大且呈现增大趋势，2010年略微下降，2015年却是急剧下降，总体来看，相较于1990年，林地的土壤侵蚀强度指数降幅为8.40%，草地降幅为10.20%。

通过研究区的自然条件(图 3)和土地利用方式的空间分布(图 1)上来看，研究区的中部和南部地形平坦，地势较低，土地利用类型以耕地和城镇为主，因此在1990—2015年间这一区域的侵蚀等级多为微度。轻度及以上的侵蚀则发生在研究区的北部，西部及东部山区的林地和草地，这些地区海拔高，坡度陡，地形复杂，是北京市水土流失重点治理区。发生在水域和城镇居民用地这两种土地利用类型的土壤侵蚀强度指数为100.00，耕地的土壤也因此侵蚀强度指数小于林地和草地，这是由于土壤侵蚀强度指数综合考虑了不同土壤侵蚀强度下不同土地利用类型的面积。为了进一步探讨耕地、林地和草地的相互转移对土壤侵蚀的影响，利用GIS软件将这3种土地利用类型与上述植被盖度分级结果进行叠加之后，考虑相同盖度等级条件下的类型转移，结果如图 5。

1990—2000年，这3种土地利用类型的相互转移比较剧烈，耕地转移至林地面积为132.44 km²，耕地转移至草地面积为51.42 km²，林地转移至草地面积为549.16 km²。1990—2015年，在相同植被盖度等级下，耕地—林地转化区、耕地—草地转化区的土壤侵蚀模数均值呈减少趋势，分别为362.32和259.54 t/(km²·a)，而林地—草地转化区的土壤侵蚀模数均值呈增加的趋势，为230.94 t/(km²·a)。这一结果与众多对于土地利用变化对土壤侵蚀的影响的研究观点一致，认为林地的水土流失量少于耕地、草地和牧草地的^[32-34]。

4 讨论与结论 4.1 讨论

此次研究采用的是研究区30年的累年年平均降雨量作为计算降雨侵蚀力因子的数据源，在计算具体年份降雨侵蚀力因子的精度上有所降低，而且RUSLE的因子间具有很强的相关性，2005年和2010年北京市潮白河流域的植被覆盖程度整体下降与降雨量之间的关笔者未作考虑。

RUSLE是国际上和国内使用较多的土壤侵蚀预报模型，但是作为经验模型，其应用有较大的局限性；而比如SWAT模型等机理模型，往往存在着数据获取困难等问题。从国内外的现状看，依然处于研究型阶段，尚未实现应用，因此模拟流域及更大尺度的土壤侵蚀的模型有待进一步的发展，以对土壤侵蚀过程、变化的研究提供更为科学有力的支撑^[35]。

4.2 结论

1) 1990—2015年，北京市潮白河流域以微度侵蚀和轻度侵蚀为主，二者占侵蚀面积的90%以上，微度侵蚀多发生在城镇居民用地、耕地、中高覆盖度及以上的林地和草地，轻度及以上程度的侵蚀则多发生在中覆盖度及以下程度的草地和林地。

2) 在相同植被盖度等级下，耕地—林地转化区、耕地—草地转化区的土壤侵蚀模数呈减少趋势，均值分别减少362.32和259.54 t/(km²·a)，而林地—草地转化区的土壤侵蚀模数呈增加的趋势，为230.94 t/(km²·a)。

3) 但总体来看，北京市潮白河流域内林地和草地的土壤侵蚀强度指数大于耕地，并在2005年到达最大值，分别为201.57和190.10。这主要是林草地质量退化致使土壤侵蚀等级较高区域的面积较大，对于山区敏感的生态环境而言，这无疑是水土流失加重的主要原因；所以，应该继续施行天然林保护，强化封山育林，充分发挥山区环境的自然修复能力，同时，在继续重视坡耕地的水土流失防治的基础上，加强对覆盖度低、垂直结构不完整的低效人工林和退耕初期、弃耕之后无人管理的草地的土壤侵蚀监管。