密云水库流域土壤侵蚀模数时空动态监测

郭 晴，任永康，陈 涛*

（1. 中国地质大学（武汉） 地球物理与空间信息学院，湖北 武汉 430074）

土壤侵蚀对土壤分级、农业生产、水文系统和环境的不利影响，长期以来被认为是影响人类可持续性发展的严重问题[1]。截至2019 年，我国土壤侵蚀总面积已达到271.08 万km2，成为世界上土壤侵蚀最严峻的国家之一。其中，水力侵蚀面积为113.47 万km2，占总面积的41.86%；风力侵蚀面积为157.61 万km2，占总面积的58.14%[2]。土壤侵蚀会引起土壤肥力衰退，进而威胁农业生产安全。径流泥沙搬移的污染源将影响侵蚀邻近区生态和经济的发展，且侵蚀泥沙的搬运使得土壤碳、氮、磷等元素的含量与组分产生变化，从而影响全球生源要素的循环，甚至成为全球气候变化的驱动要素之一[3]。

土壤侵蚀研究通常采用实地测算和构建土壤侵蚀模型来开展。由于实地观测难以直接监测空间大范围土壤侵蚀情况，且易受主观因素和监测技术的限制，因此不能满足土壤侵蚀治理工作对时效性、快速性与准确性的需求[4]；而土壤侵蚀模型不仅可以解决实地观测中的矛盾，还保障了土壤侵蚀预测和评价的准确性。根据研究方法，土壤侵蚀模型可分为物理模型和经验统计模型[5]，物理模型中最具代表性的是水蚀预报模型[6]，虽对土壤侵蚀的描述十分全面，但所需数据量大且难以在大区域中实现[7]；经验统计模型中认可度最高、应用范围最广的是Kenneth G R[8]等提出的土壤流失方程（USLE）和修正后的通用土壤流失方程（RUSLE）[6]。随着土壤侵蚀研究的不断深入，许多学者结合研究区的真实状况对模型中各因子进行了修正，并与GIS 和遥感技术进行集成，这样不仅具有了RUSLE 模型结构简单、所需数据量少的特点[9]，而且具有了GIS 和遥感技术快速获取、有效管理并实时分析时空大数据的优势，为土壤侵蚀提供了科学、快速、准确的监测和评估。例如，刘宝元[10]等根据我国国情，结合GIS 和遥感技术，辅以遥感图像，修正了RUSLE 中的各个因子；杨佳佳[11]等基于遥感和GIS 技术，利用修正后的RUSLE模型对黑龙江典型黑土区土壤侵蚀进行了监测；王丹媛[12]等利用遥感数据修改了桂南沿海诸河流域USLE 模型中的因子，并评价了土壤侵蚀强度。近年来，众多学者将RUSLE模型与GIS相结合以研究大尺度土壤侵蚀和水土保持格局与规律[13-14]，如陈峰[15]等利用RUSLE 模型研究了滇南山区土壤侵蚀的时空演变情况；王志杰[16]等运用RUSLE模型探究10 a间贵阳市土壤侵蚀空间演变规律。

作为北京市最大的饮用水源供应地，密云水库担负着北京生产生活用水的重要任务。然而，密云水库流域生态环境较脆弱，土壤侵蚀不利于流域内的水循环，进而将威胁北京市民的饮用水安全[17]。因此，研究密云水库的土壤侵蚀情况，为改进水库周边水土保持工作提供科学性的指导，对于保护其生态环境具有重要意义。目前，针对密云水库流域土壤侵蚀的研究主要为单一年份，缺乏时空动态变化的分析，在一定程度上不利于决策和管理，阻碍了流域的可持续发展。本文基于遥感和GIS技术，利用RUSLE模型对密云水库流域2001—2020年土壤侵蚀情况进行了动态监测，分析了时空动态变化特征及其影响因素，为密云水库流域的水土保持提供了科学合理的依据。

1 研究区概况

密云水库流域位于北京市城北 13 km2（41°14′~41°05′N、116°07′~117°30′E），面积为180 km2， 包括北京市的密云、怀柔、延庆以及河北省的滦平、丰宁、赤城等10个区县（图1）；坐落在潮河以及白河中游偏下，流域面积为15 788 km2。水库流域年均气温为9～10 ℃，年均降水量为488.9 mm，主要集中在汛期。流域内土壤主要包括褐土、棕壤、草甸土和粟钙土，植被以阔叶混交杂木林、油松、刺槐和落叶松为主。

图1 密云水库概要图

2 数据来源与研究方法

2.1 数据来源与预处理

本文采用的遥感数据包括2001 年、2005 年、2009 年 、 2011 年 、 2015 年 和 2020 年 6 期 共 24 景Landsat8/TM 遥感影像。由于部分年份遥感影像存在云量较多、噪声明显和数据条带丢失等问题，因此在等间隔的基础上从研究年份中选择数据精度较高的6 期数据作为研究对象；由于研究区遥感影像数据具有多时相、多光谱、多传感器等特点，因此利用ENVI 5.3对原始影像进行辐射定标、大气校正等预处理，形成标准化的多源遥感数据集。本文采用的地形数据为ASTER GDEM 30M产品；降水量数据包括2001—2020年流域内丰宁、赤城、沽源、崇礼、承德、怀来、滦平、怀柔、密云、兴隆、延庆和张家口地区的年均降水量；土地利用数据包括2001年、2005年、2009年、2011 年、2015 年和 2020 年 6 期 1 km 分辨率的中国土地利用数据（包含耕地、林地、草地、水域、建设用地和未利用地6个一级地类）；土壤数据为土壤科学数据库中的中国第二次普查数据，包括土壤质地砂砾、粉粒和有机碳含量等数据；各区县生产总值来源于当地政府官方网站数据；人口数据来源于中国统计信息网（表1）。

表1 数据列表

2.2 研究方法

基于GIS和遥感技术，本文利用RUSLE模型计算得到土壤侵蚀模数，并利用转换矩阵分析了2001—2020年土壤侵蚀强度的空间动态变化情况。RUSLE模型的计算公式为：

式中， A 为土壤侵蚀模数，单位为t/(hm2·a)；R 为降雨侵蚀力因子， 单位为(MJ·mm)/(hm2·h·a)； K 为土壤可蚀性因子 ，单位为(t·hm2·h)/(hm2·MJ·mm) ； LS 为坡度坡长因子，无量纲；C 为植被覆盖与管理因子，无量纲；P 为水土保持措施因子，无量纲。

1）降雨侵蚀力因子表示降雨引起土壤侵蚀的潜在能力。本文基于监测区域代表站点连续30 min的最大降雨强度和每次降雨的降雨量，利用基于日降雨量的拟合模型[18-19]计算降雨侵蚀因子。为了得到流域内降雨侵蚀力因子的连续数据，在各测站降雨侵蚀力值和空间分布的基础上，利用GIS 的样条插值实现。其计算公式为：

式中，Pf为汛期雨量；I30B为侵蚀性降雨最大30 min雨强的年代表值，单位为mm；i 为侵蚀性降雨次数；Pt为汛期月份的月降雨总量。

2）土壤可蚀性因子是指土壤对侵蚀的敏感性，与降雨、土壤自身的理化特性等相关。采用Wisch⁃meier W H[20]等提出的公式，根据土壤质地、土壤中的有机质含量等主要土壤性质数值指标估算土壤可蚀性因子。其计算公式为：

式中， M 为优势粒径组成的乘积；a 为有机质含量百分比；b 为土壤结构等级；c 为土壤渗透等级。

3）坡度坡长因子。坡长因子是指任意坡长的单位面积土壤流失量与标准坡长单位面积土壤流失量之比；坡度因子是指任意坡度下的单位面积土壤流失量与标准坡度下单位面积土壤流失量之比。由于坡度从20%增加到40%和60%时，坡长指数无变化[6]，因此当坡度<21%时，采用式（5），当坡度＞21%时，采用式（6）[21]。

4）植被覆盖与管理因子是指在土壤、坡度和降雨条件相同的情况下，具有特定植被覆盖的土地的土壤流失量与裸地的土壤流失量的比值。本文采用蔡崇法[22]等提出的利用植被覆盖度与定量估算管理因子的回归方程，即

式中， f 为植被覆盖度；C 为定量估算管理因子

植被覆盖度计算采用李苗苗[23]等提出的改进的像元二分模型，则有：

式中， NDVI 为归一化植被指数值，无量纲；NDVImax、NDVImin为遥感图像中给定置信度置信区间内的最大值与最小值，本文采用的置信区间为[5%，95%]。

5）水土保持措施因子是指采取一定的水保措施后的土壤流失量与标准状况下土壤流失量的比值。水土保持措施因子目前主要通过布设天然小区试验得到。由于难以获得水土保持措施因子，在考虑财力物力人力等综合因素后，结合前人的研究结果[24-25]，本文将水土保持措施因子设定为1。

2.3 转移矩阵

土壤侵蚀转移矩阵反映了研究区在某一时段期初和期末不同侵蚀等级面积之间的转化信息，包括期初各类面积的转出量和期末各类面积的转入量[26]。本文计算得到2001 年和2020 年不同侵蚀等级土壤侵蚀模数的转移矩阵，并探究了侵蚀等级的变化特征。

2.4 相关性分析

相关性分析能反映因变量与自变量之间相互关系的密切程度。本文将土壤侵蚀模数分别与年均降雨量与平均植被覆盖度进行统计回归分析[27]，以探究年均降雨量、平均植被覆盖度与土壤侵蚀模数的相关性。

3 研究结果分析

3.1 土壤侵蚀模数的时间变化特征

密云水库流域年均土壤侵蚀模数呈“快速减小、小幅增长、又迅速回落”的趋势，20 a间研究区土壤侵蚀状况整体得到明显控制，2009—2015年土壤侵蚀呈现恶化趋势。根据国家水利部发布的SL190-2007《土壤侵蚀分类分级标准》[28]，本文将研究区土壤侵蚀分为微度、轻度、中度、强烈、极强烈和剧烈侵蚀6个等级，并统计分析不同时期密云水库流域土壤侵蚀强度的面积比例。20 a间研究区土壤侵蚀以轻度和微度侵蚀为主，二者占总面积的70%以上，且呈波浪式变化；其次为中度侵蚀；除强烈侵蚀和极强烈侵蚀面积以外，其他侵蚀强度面积呈总体下降趋势。2001—2020年各侵蚀强度等级的年变化率（降幅）均在5%~10%，但2001—2020年强烈和极强烈侵蚀增幅明显（图2）。

图2 2001—2020年不同等级平均土壤侵蚀模数折线图

由面积转移矩阵可知（表2）， 2001—2020 年研究区土壤侵蚀以微度侵蚀少量转出、其他侵蚀等级不同程度转入为主要特征；侵蚀强度等级的转出方向以向低等级侵蚀转化为主要特征，微度侵蚀转入了14.69 km2，占总微度侵蚀转出面积的67.35%，轻度侵蚀转入微度侵蚀的面积为2.91 km2，占轻度侵蚀转出面积的56.7%，中度侵蚀转入微度和轻度侵蚀的面积为0.43 km2，占总中度侵蚀转出面积的67.2%，剧烈侵蚀几乎全部转向低等级侵蚀；但极强烈侵蚀大部分面积转入为更高等级，强烈侵蚀转向低等级侵蚀的面积占比为44.7%，因此流域内20 a 间土壤侵蚀呈“整体好转、局部恶化”的趋势。

表2 2001—2020年研究区土壤侵蚀强度面积转移矩阵/km2

3.2 土壤侵蚀模数的空间变化特征

3.2.1 研究区不同行政区内土壤侵蚀变化分析

本文叠置分析了密云水库流域研究期土壤侵蚀强度分布图与研究区行政区划图，获得流域内研究期间各县区土壤侵蚀强度的面积比例（表3），可以看出，流域内土壤侵蚀最严重的区域为密云区和滦平县，均值都在1 500 t/(km2·a)以上；密云、崇礼和滦平县的年均侵蚀模数降幅最大，分别为53.9%、63.5%和42.5%；崇礼县水土保持情况最好，平均土壤侵蚀模数接近 700 t/(km2·a)，且降幅较大。

表3 密云水库流域2001—2020年各县区侵蚀强度面积/km2

3.2.2 区域内不同海拔下土壤侵蚀变化分析

本文叠置分析了土壤侵蚀分布图与数字高程模型分级图，得到密云水库流域2001—2020年不同海拔下土壤侵蚀变化情况（表4、5），可以看出，随着海拔的升高，研究区土壤侵蚀模数整体呈“高—低—高”的特点，海拔500 m 以下以及海拔1 500 m 以上的地区侵蚀相对严重。在低海拔地带，由于地形平坦有利于开发，导致人口分布密集、城市化程度较高，进而对地表植被造成严重破坏；而在1 500~3 000 m高海拔地区，土地类型以高山裸土区和荒草地为主，地势陡峻，高强度的降雨易产生滑坡、泥石流等自然灾害。

表4 密云水库流域2001—2009年不同海拔下土壤侵蚀变化分析

表5 密云水库流域2011—2020年不同海拔下土壤侵蚀变化分析

3.2.3 区域内不同土地利用类型下土壤侵蚀变化分析

本文统计分析了不同土地利用类型的侵蚀状况，结果如表6所示。各种土地利用类型中，耕地的平均土壤侵蚀强度普遍较高，而城乡居民用地和水域的平均土壤侵蚀强度较低。2001年、2009年耕地和林地的侵蚀强度接近于未利用地，处于较高水平；2015年耕地、林地和草地的侵蚀强度均大于1 500 t/（km2·a），林地的侵蚀强度最大，未利用地的侵蚀强度未达到1 500 t/（km2·a）；2020年耕地、林地和草地的侵蚀强度均有所回落，与逐步治理耕地以及大力推行退耕还林政策密不可分。

表6 密云水库流域2001—2020年不同土地利用类型下土壤侵蚀变化/（t/（km2·a））

3.3 影响因素分析

3.3.1 自然因素

本文对密云水库流域研究期间的平均土壤侵蚀模数与平均植被覆盖变化、平均降雨侵蚀模数进行回归分析，结果表明，平均土壤侵蚀模数与平均植被覆盖变化存在明显的负相关性，相关系数为-0.77（图3）；而与平均降雨侵蚀模数呈正相关性，相关系数为0.71（图4）。由于平均土壤侵蚀模数与植被覆盖度的决定系数大于其与降雨侵蚀模数的决定系数，因此平均土壤侵蚀模数与植被覆盖度的相关程度更高。

图3 植被覆盖度与土壤侵蚀模数关系散点图

图4 平均降雨侵蚀模数和土壤侵蚀模数关系散点图

3.3.2 人为因素

分析研究区各县市人口数据（图5）和土壤侵蚀强度的面积比例统计表可知，研究期内密云区、兴隆县、滦平县和怀柔区的平均侵蚀面积均在1 000 km2以上，且人口逐年递增，其中密云区和怀柔区的人口基数大、增速也较快，人口的迅速增长促使农民大量开荒种粮，导致土地垦荒率越来越高，加剧了水土流失；后期（2015—2020年）崇礼县、沽源县和延庆区等地区人口出现零增长甚至负增长的情况，人口对环境的压力下降了，水土流失也有所缓解。

图5 研究区各县区常驻人口变化图

分析各区县GDP 变化情况（图6）和土壤侵蚀强度的面积比例统计表可知，土壤侵蚀与社会经济发展具有相关性，研究区初期和中期（2001—2011年）各行政区的GDP增速较大，如怀柔区、密云区和滦平县，其土壤侵蚀面积相对较大，原因在于其地理位置离北京主城区较近，经济社会活动较频繁，经济的发展伴随着开发区建设等土建工程的急剧增加，使得大量草地、林地遭到了破环，加剧了水土流失；后期各行政区GDP增速逐渐放缓，政府践行绿色发展理念，加大对水土保持工作的投入，使得土壤侵蚀面积较前期有所回落。

图6 研究区各县区GDP变化图

4 结 语

本文利用RUSLE 模型对密云水库流域2001—2020 年的土壤侵蚀情况进行了定量研究。结果表明，2001—2020年密云水库流域土壤侵蚀模数大致经历了先下降后回升又大幅下降的变化过程，最小值出现在2009 年，2020 年比2001 年的土壤侵蚀面积有所减少，说明整体情况在好转；流域内侵蚀等级主要由轻度转为微度，但强烈侵蚀和剧烈侵蚀主要集中在滦平县和密云区；随着海拔的上升，流域内侵蚀模数先升高后降低又升高；土壤侵蚀模数与土地利用类型相关性程度较高；植被覆盖度升高，土壤侵蚀模数下降，降雨侵蚀加剧，土壤侵蚀模数降低，但前者的影响程度较大；人口与社会经济的发展也会对土壤侵蚀模数造成影响。