USLE土壤流失方程是经验统计模型中出现最早且影响最广的模型，Wischmeier归纳出通用土壤流失方程USLE^[1]，被广泛应用于土壤侵蚀预测、土壤侵蚀时空变化、评估区域土壤侵蚀风险等领域^[2-3]。北京市房山区现有的土壤侵蚀研究主要集中在水土保持效益评价^[4]以及不同雨型对坡面土壤侵蚀的影响分析^[5]，欠缺针对房山区土壤侵蚀对典型地形因子响应的研究；因此，笔者采用通用土壤流失方程USLE，结合GIS与RS技术对北京市房山区2013—2018年土壤水力侵蚀进行模拟计算，定量评估北京市房山区土壤侵蚀分布格局与动态特征，探索土壤侵蚀与地形因子的关系，并为研究区土壤侵蚀防治的重点布局决策提供参考。

北京市房山区地处华北平原和太行山交界地带，介于E 115°25′~116°15′，N 39°30′~39°55′之间，总面积2 019 km²。房山区气候属暖温带半湿润大陆性季风气候，四季特征分明，春季干旱多风，夏季炎热多雨，多大到暴雨且雨量集中；多年平均气温10.8~11.6 ℃，多年平均降雨量为602.5~645.2 mm。房山区地貌西北高，东南低，由西北向东南依次为中山、低山、丘陵、岗台地、冲(洪)积平原，其中66%的区域为山地丘陵；海拔区间为-8~2 014 m，平均381.92 m。根据南京土壤所提供的1∶100万土壤数据分析得知，房山区土壤类型以褐土为主，占全区面积的51.91%，集中分布于研究区南部及西部，土壤质地以壤质为主，其中轻壤质占78.80%，中壤质和砂壤质分别占全区总面积的6.14%和9.42%。房山区植物茂盛，种类众多，植物种类共137科、545属、1 064种。区域内土地利用类型以农用地为主，其中林地占农用地总面积的63.24%，分布于全区的西部和北部山区，耕地占农用地面积的25.50%，分布于东部平原地区。根据水利部办公厅关于印发《全国水土保持规划国家级水土流失重点预防区和重点治理区复核划分成果》的通知(办水保[2013]188号文)，房山区为太行山国家级水土流失重点治理区；房山区属于全国水土保持规划(2015—2030)中北方土石山区(北方山地丘陵区)太行山东部山地丘陵水源涵养保土区，水土流失以微度和轻度水力侵蚀为主。

笔者利用GIS，结合多源遥感数据建立土壤侵蚀方程，计算研究区土壤侵蚀模数及侵蚀量，并划分土壤侵蚀等级，分析房山区土壤侵蚀状况。主要数据及数据来源见表 1所示。

表 1(Tab. 1)

表 1 数据及来源 Tab. 1 Data and their sources 数据类别 Data category 来源 Sources 精度 Precision 气象数据 Meteorological data 中国气象数据共享网以及国家气象中信息中心 China Meteorological Data Sharing Network and National Meteorological Information Center 2013—2018年降水日值数据 Daily rainfall data from 2013 to 2018 植被指数数据 Vegetation index data 中科院资源环境数据中心 Chinese Academy of Sciences Resources and Environment Data Center 500 m旬合成数据 500 m resolution 10-day synthetic data 土地利用数据 Land use data 资源环境科学数据中心 Resource and Environment Science Data Center 空间分辨率1 000 m Spatial resolution 1 000 m 土壤属性数据 Soil attribute data 南京土壤所提供的1∶100万土壤数据 1∶1 million soil data provided by Nanjing Institute of Soil Science, 空间分辨率1 000 m Spatial resolution 1 000 m 数字高程模型 Digital elevation model (DEM) 资源环境数据云平台 Resource environment data cloud platform 空间分辨率30 m Spatial resolution 30 m

表 1 数据及来源 Tab. 1 Data and their sources

通用土壤流失方程由6个影响因子组成^[1]，其中降雨侵蚀力因子R由气象数据计算并插值后得出；土壤可蚀性因子K先根据土壤属性划分出研究区各土壤类型，然后根据土壤类型计算得出；坡长、坡度因子LS由数字高程模型计算；作物覆盖和管理因子C由植被指数数据计算得出；水土保持措施因子P由土地利用数据查表得出；最后将6项因子相乘得到土壤侵蚀模数A，其中各因子含义、计算公式与单位见表 2。

表 2(Tab. 2)

表 2 各因子含义、计算公式与单位 Tab. 2 Corresponding meanings, calculation formulas and units of factors 因子 Factor 含义 Meaning 公式 Formula 单位 Unit A 土壤侵蚀模数 Soil erosion modulus A=R×K×L×S×C×P t/(hm2·a)(or t/(km2·a)) R 降雨侵蚀力因子[6] Rainfall erosivity factor $R=0.698 \sum\limits_{j=1}^{k}\left(p_{j}\right)^{1.671} $ (MJ·mm)/(hm2·h·a) K 土壤可蚀性因子[7] Soil erodibility factor $ \begin{gathered} K_{a}=\left\{0.2+0.3 \exp \left[0.0256 \mathrm{SAN} \times\left(1-\frac{\mathrm{SIL}}{100}\right)\right]\right\} \\ \left(\frac{\mathrm{SIL}}{\mathrm{CAL}+\mathrm{SIL}}\right)^{0.3}\left(1.0-\frac{0.25 C}{C+\exp (3.72-2.95 C)}\right) \times \\ \left(1.0-0.7 \mathrm{SN}_{1}\right) /\left[\mathrm{SN}_{1}+\exp \left(-5.51+22.9 \mathrm{SN}_{1}\right)\right] \\ \mathrm{SN}_{1}=1-\mathrm{SAN} / 100 \\ K=K_{\mathrm{a}} \times 0.1317 \end{gathered}$ (t·hm2·h)/(hm2·MJ·mm) L、S 坡长因子、坡度因子[8]，用于表述研究区地形 Slope length factor, slope factor. Used for describing the topography of the study area $\begin{gathered} S=\left\{\begin{array}{l} 10.8 \sin \theta+0.03 \quad \theta \leqslant 5^{\circ} \\ 16.8 \sin \theta-0.05 & 5^{\circ}<\theta \leqslant 10^{\circ} \\ 21.9 \sin \theta-0.96 & \theta \geqslant 10^{\circ} \end{array}\right. \\ n= \begin{cases}0.2 & \text { if }~~ \theta \leqslant 1^{\circ} \\ 0.3 & \text { if }~~ 1^{\circ}<\theta \leqslant 3^{\circ} \\ 0.4 & \text { if }~~ 3^{\circ}<\theta \leqslant 5^{\circ} \\ 0.5 & \text { if }~~ \theta>5^{\circ}\end{cases} \\ L=\left(\frac{\lambda}{22.13}\right)^{n} \end{gathered} $ 1 C 作物覆盖和管理因子[9] Crop coverage and management factors $ \begin{array}{l} G = \frac{{N - {N_{{\rm{soil }}}}}}{{{N_{{\rm{veg }}}} - {N_{{\rm{soil }}}}}}\\ \left\{ {\begin{array}{*{20}{l}} {C = 1}&{G = 0}\\ {C = 0.6508 - 0.3436{{\log }_{10}}G}&{0 < G < 78.3\% }\\ {C = 0}&{G \ge 78.3\% } \end{array}} \right. \end{array}$ 1 P 水土保持措施因子[10] Soil and water conservation measures factor 1 注：pj为日降水量，mm；SAN、SIL、CAL分别为土壤砂粒、粉粒、黏粒含量，%；Ka为土壤侵蚀力因子(美制单位)；θ为坡度，(°)；n为坡长指数，量纲1；λ为坡长，m；G为植被覆盖度，%；Nsoil为裸地NDVI值，Nveg为高纯度植被像元NDVI值。Notes：pj is the daily rainfall amount, mm; SAN, SIL, and CAL are the soil sand, silt, and clay content, respectively, %; Ka is the soil erosivity factor (US unit); θ is the slope, °; n is the slope length index; λ is slope length, m; G is vegetation coverage, %; Nsoil is the NDVI value of bare land, and Nveg is the NDVI value of high-purity vegetation pixels.

表 2 各因子含义、计算公式与单位 Tab. 2 Corresponding meanings, calculation formulas and units of factors

根据土壤流失方程计算得出2013—2018年房山区土壤侵蚀量，各年土壤侵蚀总量逐渐减少，从2013年的最高70.694万t降低至2018年的29.082万t，其余各年土壤侵蚀量为：2014年53.138万t、2015年42.699万t、2016年48.504万t、2017年43.509万t(图 1)。根据SL190—2007《土壤侵蚀分类分级标准》将土壤水力侵蚀划分为微度、轻度、中度、强烈、极强烈、剧烈6个不同的等级，微度侵蚀等级各年平均占比64.13%，其次为轻度侵蚀占比29.99%，其余各更高侵蚀强度等级面积占比 < 6%，且各年微度侵蚀面积占比逐步上升，由2013年的55.74%至2018年的71.70%。综上2013年至2018年间土壤侵蚀强度的变化，表明房山区土壤侵蚀状况逐步好转。

图 1(Fig. 1)

图 1 房山区不同年度土壤侵蚀量 Fig. 1 Soil erosion in different years in Fangshan district

高程对于土壤侵蚀量具有显著影响，不同高程带上地貌、降雨、植被类型特征均直接影响土壤侵蚀量^[11]。根据研究区的高程范围(21~2 019 m)，参考类似研究^[12]的做法，以区分高程梯度为目的，将房山区所涉及高程均分 < 400 m、[400，800) m、[800，1 200) m、[1 200，1 600) m、≥1 600 m 5个范围。并与各年房山区土壤侵蚀模数分布图进行叠加分析，得到不同高程带的土壤侵蚀特征。

从平均土壤侵蚀模数来看(图 2)，各年[400, 800) m高程的土壤侵蚀模数均最高，其次为 < 400 m高程范围，[400，1600) m高程范围平均土壤侵蚀模数随高程升高而下降。

图 2(Fig. 2)

图 2 各年各高程带土壤平均侵蚀量模数与土壤侵蚀量 Fig. 2 Soil erosion modulus and soil erosion amount of each elevation grade in each year

土壤侵蚀量方面(图 3)，各高程带土壤侵蚀量随着高程升高而减小，40%以上侵蚀量都发生在 < 400 m高程一带，且除 < 400 m高程带外，各高程带的土壤侵蚀量占比总体上逐渐减少，说明≥400 m高程带的土壤侵蚀状况恢复较明显。

图 3(Fig. 3)

图 3 各年各高程带土壤侵蚀量占比 Fig. 3 Ratio of soil erosion amount at each elevation grade in each year

按照GB/T 15772—2008《水土保持综合治理规划通则》，将地表坡度划分为微坡[0°，5°)、较缓坡[5°，8°)、缓坡[8°，15°)、较陡坡[15°，25°)、陡坡[25°，35°)、极陡≥35°6个等级并与房山区土壤侵蚀模数分布图进行叠加分析，得到不同坡度带的土壤侵蚀特征。

从土壤侵蚀模数上来看(图 4和图 5)，在[0°，25°)坡度等级，平均土壤侵蚀模数随着坡度增加而增大，坡度越陡的地区在难以形成植被覆盖的条件下，潜在的土壤侵蚀量也较大。但在≥25°的坡度等级下，土壤侵蚀模数出现了一定程度的下降，其原因是坡度过陡时，降水在这一坡度等级地点停留时间缩短，快速形成径流转移至其他地点。在2013年较陡坡平均土壤侵蚀模数达到各年各坡度等级最高的841.40 t/(km²·a)，但随着时间推移，各坡度带上的土壤侵蚀模数都出现逐步下降的情况。

图 4(Fig. 4)

图 4 各年各坡度等级土壤侵蚀模数与土壤侵蚀量变化图 Fig. 4 Soil erosion modulus and soil erosion amount of each slope grade in each year

图 5(Fig. 5)

图 5 各年各坡度等级土壤侵蚀量变化图 Fig. 5 Soil erosion amount of slope grade in each year

从土壤侵蚀量上来看，主要的侵蚀量分布于[0°，5°)的微坡等级区域，该等级的土壤侵蚀量占到总侵蚀量>73%，微坡、较缓坡、缓坡等级在各年的侵蚀量分布都相对均匀。从各坡度带侵蚀量占比年际变化上来看，各年各坡度等级侵蚀量分布较为均衡，但2018年微坡等级侵蚀量占比较2013年有一定降低。总体上由于微坡等级面积较大，则微度等级土壤侵蚀量也较大；其他坡度等级因面积较小，尽管土壤侵蚀模数较微坡等级大，但其土壤侵蚀量小于总量的30%。

1) 北京市房山区的年度土壤侵蚀总量，随着时间推移逐渐减少，自2013年至2018年土壤侵蚀总量减少58.86%，研究区的土壤侵蚀总量明显降低。

2) 在不同侵蚀强度面积中，房山区以微度侵蚀面积占比最大，各年度微度侵蚀面积占比均在55%以上；随时间推移微度侵蚀面积逐渐增加，自2013年55.74%增至2018年71.70%，则研究区的土壤侵蚀面积明显减少。

3) 在不同高程范围内，以[400，800) m高程区间的土壤侵蚀模数较高，以[0，800) m高程区间的土壤侵蚀量占比较大，各年占比均>60%。

4) 在不同坡度范围内，以[15°，25°)较陡坡等级的土壤侵蚀模数较高，以[0°，5°)微坡等级的土壤侵蚀量占比较大，各年微坡等级的土壤侵蚀量占比>73%。

本研究结论还局限于基于1 000 m卫星影像获取土地利用等数据的计算结果，研究区在更高分辨率(如2 m等)影像下表现的土壤侵蚀分布特征及其与本结果的差别还需要进一步研究。