淮河上游土壤侵蚀估算与分析

Estimation and Analysis of Soil Erosion in the Upstream of Huaihe River

摘要: 基于地理信息系统(GIS)和遥感(RS)技术，利用修正的通用土壤流失方程(RUSLE)模型，计算了2000~2015年的淮河上游息县流域的土壤侵蚀，分析了土壤侵蚀的时空变化特征。研究表明：近15年来     植被覆盖管理因子和水土保持措施因子基本呈现减小的趋势。主要受降水侵蚀因子的影响，流域土壤侵蚀模数总体呈减小趋势。土壤侵蚀强度主要以微弱侵蚀为主，占总面积的92.49%，主要发生在高程在300 m以下或坡度在15˚以下的区域。

Abstract: Based on Revised Universal Soil Loss Equation (RUSLE), Geographic Information Systems (GIS)and (Remote Sensing) RS spatial analysis techniques, this study calculated the soil erosion modulus in Xixian Watershed of the upstream of Huaihe River from 2000 to 2015, and analyzed the spatial and temporal changes of soil erosion intensity. The results indicated that vegetation cover management factor soil and water conservation support practice factor showed a decreasing trend in the past 15 years. The soil erosion modulus decreased in the main influence factor of rainfall erosivity. The soil erosion in this region as a whole was dominantly trivial in degree with 92.49%, mainly in the regions below 300 m in elevation or 15˚ in gradient.

1. 引言

土壤侵蚀是土壤及其母质在水力、风力、冻融、重力等外力作用下被破坏、剥蚀、搬运和沉积的过程 [1]，作为全球最大的生态环境问题之一 [2]，不仅导致土地质量退化降低农业生产力，也是导致生态环境恶化的重要原因 [3]。土壤侵蚀研究已成为当今全球环境变化研究的重要内容 [4]。

传统土壤调查方法不具有实效性，特别是在大区域尺度下，难以确定水土流失的动态变化情况。基于RS和GIS技术，结合土壤侵蚀模型是定量评估区域土壤侵蚀、分析土壤侵蚀空间分布特征、防治区域水土流失的有效手段 [5]，如通用土壤流失方程(Universal Soil Loss Equation, USLE) [6] 、修订土壤流失方程RUSLE (revised universal soil loss equation) [7] 、中国土壤流失方程(Chinese Soil Loss Equation, CSLE) [8] 、美国的WEPP [9] 、欧洲的 EUROSEM [10] 等模型。由于很多地区数据资料不够完善，以RUSLE为代表的经验统计模型仍是应用最为广泛的模型 [11] [12] [13]，该模型研究了与土壤侵蚀密切相关的各种自然和人文因子的关系，其结构形式简单，参数易获取且计算方便，被广泛应用于土壤侵蚀定量评价。

淮河流域的水土流失一直都是淮河治理所面对的主要问题，特别是上游桐柏大别山区，由于长期以来对自然资源的不合理开发利用，导致许多地区原始植被破坏，土层瘠薄、抗蚀能力差，随着人口的不断增长和社会发展，人口–资源–环境的矛盾加剧，已被列为国家级水土流失重点预防区 [14]。本文以淮河上游息县水文站以上流域为研究区域，利用2000、2005、2010和2015年的遥感影像和其他资料，采用RUSLE模型研究淮河上游的土壤侵蚀时空分异特征，为区域水土保持措施提供决策支持。

2. 研究区概况

息县流域位于淮河上游区域，流域范围为113˚15'E~114˚50'E，31˚28'N~32˚44'N (图1)，控制流域面积为10,190 km²。流域地处北亚热带过度气候，多年平均气温在15℃左右，淮河以南多年平均降水量为800~1400 mm，而淮河以北区域多年平均降水量为500~900 mm，降雨时空分布不均匀，年际、区间差异分布明显。流域内地形以低平的平原和缓丘为主，小部分为平原洼地。

3. 材料与方法

3.1. 数据来源

本研究需要的数据包括：1) 降雨量数据；2) DEM数据；3) 土壤类型数据；4) 土地利用数据；5) NDVI数据。相关数据说明见表1。

Figure 1. Location of the study areas

图1. 研究区位置

Table 1. Data sources and description

表1. 数据来源与说明

3.2. 土壤侵蚀模型

本文采用目前应用最为广泛的修正土壤流失方程(RUSLE)来进行研究区土壤侵蚀的量化，计算公式为：

$A=R\cdot K\cdot LS\cdot C\cdot P$ (1)

式中：A为预测土壤侵蚀量(t/hm²∙a)；R为降水侵蚀力因子((MJ∙mm)/(hm²∙h∙a))；K为土壤可蚀性因子((t∙hm²∙h)/(MJ∙hm²∙mm))；LS为坡长坡度因子，无量纲；C为植被覆盖与管理因子，无量纲；P为水土保持措施因子，无量纲。各个因子计算公式为：

1) 降水侵蚀力因子R

采用月降水量数据和年降水量数据计算：

$R={\displaystyle {\sum }_{i=1}^{12}1.735\times {10}^{\left(1.5\ast \lg \frac{P_i^2}{P}-0.8188\right)}}$ (2)

式中：P_i为月降雨量(mm)；P为年降雨量(mm)。根据研究区气象站点的月、年的降雨数据计算各个站点的降雨侵蚀力，再用kringing插值法进行空间插值为30 m的栅格数据。

2) 土壤可蚀性因子k

$\begin{array}{c}k=\left\{0.2+0.3\exp \left[-0.0256SAN\left(1-\frac{SIL}{100}\right)\right]\right\}{\left(\frac{SIL}{SIL+CLA}\right)}^{0.3}\\ \times \left(1-\frac{0.25C}{C+\exp \left(3.72-2.95C\right)}\right)\left(1-\frac{0.75SN}{SN+\exp \left(22.9SN-5.51\right)}\right)\end{array}$ (3)

式中，SAN表示砂粒含量(%)；SIL表示粉粒含量(%)；CLA表示黏粒含量(%)；C 表示有机碳含量(%)； $SN=1-SAN/100$。利用土壤数据库的相关数据，计算研究区的土壤可蚀性因子k的空间分布。

3) 植被覆盖管理因子C

本文选用蔡崇法等 [15] 建立的植被覆盖管理因子C与植被覆盖度f_c之间的关系，计算公式为：

$C=\{\begin{array}{l}1f_c=0\\ 0.6508-0.3436\ast \lg \left(f_c\right)0<f_c<78.3\%\\ 0f_c\ge 78.3\%\end{array}$ (4)

$f_c=\frac{\text{NDVI}-{\text{NDVI}}_{\min }}{{\text{NDVI}}_{\max }-{\text{NDVI}}_{\min }}$

式中：f_c为植被覆盖度。NDVI为归一化植被指数，NDVI_max和NDVI_min为归一化植被指数中的最大值和最小值。本研究采用MODIS的16天NDVI产品，计算每年的年均NDVI值，NDVI_max和NDVI_min分别取年均NDVI最大1%和最小1%的值。

4) 坡长坡度因子LS

由DEM数据分别计算坡度因子和坡长因子，计算公式为：

$S=\{\begin{array}{l}10.8\sin \theta +0.03\theta <5^{\circ }\\ 16.8\sin \theta -0.505^{\circ }\le \theta <{10}^{\circ }\\ 21.9\sin \theta -0.96\theta \ge {10}^{\circ }\end{array}$ (5)

$L={\left(\frac{\lambda }{22.13}\right)}^m$ (6)

式中，L为坡长因子；S为坡度因子； $\theta$ 为坡度值(˚)； $\lambda$ 为坡长(m)；m为随坡度变化的量，当 $\theta \le 1^{\circ }$ 时， $m=0.\text{2}$ ；当 $1^{\circ }<\theta \le 3^{\circ }$ 时， $m=0.3$ ；当 $3^{\circ }<\theta \le 5^{\circ }$ 时， $m=0.4$ ；当 $\theta >5^{\circ }$ 时， $m=0.5$。

5) 水土保持措施因子P

本研究结合以往研究成果 [6] [11] [12] 和土地利用情况，对研究区的不同土地类型赋予不同的P因子值(表2)。

Table 2. The P value of different land use types

表2. 不同土地利用类型的P值

3.3. 土壤侵蚀分级标准

根据国家水利部颁布的《土壤侵蚀分类分级标准》(SL190-2007) [16]，将计算得到的土壤侵蚀量分为6个等级：1) 微度侵蚀：小于500 t/km²∙a；2) 轻度侵蚀：500~2500 t/km²∙a；3) 中度侵蚀：2500~5000 t/km²∙a；4) 强度侵蚀：5000~8000 t/km²∙a；5) 极强侵蚀：8000~15000 t/km²∙a；6) 剧烈侵蚀：大于15000 t/km²∙a。

4. 结果与分析

4.1. 各个侵蚀因子计算

根据RUSLE模型的五个因子的计算公式，分别计算2000、2005、2010年和2015年各因子的结果。由于流域的地形特征和土壤特性基本不变，各年计算的坡长坡度因子和土壤可蚀性因子都一致。降水侵蚀因子反映了降水对土壤侵蚀的重要驱动因素。图2给出了2000年至2015年研究区的年降水量和降水侵蚀因子的变化，虽然二者都呈现下降趋势，但在相同的年降水条件下，2005年的降水侵蚀因子明显大于2000年的降水侵蚀因子，2010年的降水侵蚀因子明显大于2015年的降水侵蚀因子，主要原因是2005年和2010年汛期雨量较为集中。表明降水量和降水强度对降水侵蚀因子有重要影响。

Figure 2. Annual rainfall and rainfall erosivity factor

图2. 年降水量与降水侵蚀力因子

植被覆盖管理因子反映的是植被对水土流失的抑制程度，植被覆盖程度越高对水土保持能力越好、对水土流失的效果就越强。图3给出了2000年至2015年研究区NDVI和植被覆盖管理因子的变化。近15年来，研究区年均NDVI基本呈现上升趋势，这与近年来国家对淮河上游生态保护的重视有关，退耕还林、封山育林工程，使植被覆盖状况好转。植被管理因子的变化趋势正好跟NDVI变化趋势相反。

不同土地类型具有不同的水土保持因子。表3给出了研究区2000、2005、2010和2015年土地利用类型面积比例。研究区以耕地和林地为主，分别占到了55.1%和35.6%。近15年来，耕地和草地面积比例变化不大，林地呈现下降趋势，水域和建设用地呈现上升趋势，但变化比例并不大，仅在0.8%左右。进一步计算流域2000、2005、2010和2015年水土保持措施因子分别为：0.542、0.540、0.539和0.529。水土保持措施因子近15年来呈下降趋势，但变化较小。

4.2. 土壤侵蚀空间分布特征

将各个土壤侵蚀因子相乘，得到2000、2005、2010和2015年淮河流域上游息县流域的平均侵蚀模数分别为258.9、279.5、173.8和109.1 t/km²∙a，结合流域面积得出土壤侵蚀总量为263.8、284.8、177.1和111.1万t。邹宏荣等 [17] 指出息县流域在2000年和2005年的实测输沙量为294.5和322.4万t，计算结算与实测比为89.6%和88.3%，说明利用RUSLE计算出的息县流域土壤侵蚀结果较可靠。

Figure 3. NDVI and vegetation cover management factor

图3. NDVI与植被管理因子

Table 3. Land use change from 2000 to 2015

表3. 2000~2015年土地利用变化

根据土壤侵蚀分类分级标准，计算2000~2015年各侵蚀强度所占整个流域的侵蚀面积百分比(表4)。淮河上游流域土壤侵蚀面积随侵蚀强度增大显著减小，主要以微度侵蚀为主，轻度侵蚀为辅，中度及其以上侵蚀的面积较小、仅占1%左右。从2000年至2015年，微度侵蚀面积比例呈增加趋势，其他侵蚀轻度面积比例呈现减小趋势。2015年流域微度侵蚀面积占到了97.77%，中度及其以上侵蚀的面积比例仅为0.33%。降雨、植被覆盖和土地利用是影响土壤侵蚀动态变化的主要因素。植被覆盖管理因子和水土保持措施因子近15年来基本呈现为减小的趋势，但变小的比例并不大。而流域土壤侵蚀量变化趋势与降水侵蚀因子变化趋势一致，2000年的平均侵蚀模数与2015年相比下降了57.9%，而对应的降水侵蚀因子也下降了54.7%。

Table 4. Statistics of soil erosion in study area from 2000 to 2015

表4. 2000~2015年研究区土壤侵蚀统计表

图4给出了各年的土壤侵蚀空间分布图，总体呈现西部和南部大、中东部小的趋势。流域下游平原地区以耕地为主、地形平坦，大部分为微弱侵蚀区域，土壤侵蚀严重地区主要分布在西部桐柏山区和南部大别山区，虽然多为植被覆盖度，近年来植被覆盖度呈现上升趋势，但树龄、林种不合理，以中幼林为主，成熟林较少，林种则以经济林为主，使得保持水土、涵养水源等生态功能脆弱。同时山区地形更加复杂，容易发生土壤侵蚀。2010年和2015年降水侵蚀因子下降，侵蚀严重的地区主要分布地形更为复杂的沟谷地区。刘瑞娟等 [18] 研究也表明淮河流域地形起伏度与土壤侵蚀量呈较显著的正相关。

Figure 4. Spatial distribution of soil erosion in study area from 2000 to 2015

图4. 2000~2015年研究区土壤侵蚀空间分布图

4.3. 土壤侵蚀与地形的关系

将2000、2015、2010和2015年土壤侵蚀量结果计算平均值作为流域多年平均土壤侵蚀量，并重新进行土壤侵蚀分级，分析土壤侵蚀在不同DEM和坡度区间的分布特征。表5给出了在6个海拔高度区间(<100 m，100~200 m，200~300 m，300~500 m，500~700 m，>700 m)的土壤侵蚀面积分布情况。研究区土壤侵蚀强度主要以微弱侵蚀为主，占总面积的92.49%，轻度、中度及中度以上分别为6.94%和0.57%。在不同的高程海拔区间，0~100 m高程区间土壤侵蚀面积最大，占侵蚀总面积的43.41%，其次是100~200 m高程区间，其土壤侵蚀面积占土壤侵蚀总面积的36.17%。高程在300 m以下的面积占土壤侵蚀总面积的90.03%。700 m以上的土壤侵蚀面积仅占土壤侵蚀总面积的0.26%。中度、强度、极强度的土壤侵蚀主要发生在100~300 m高程区间，并以此为中心向上、向下土壤侵蚀面积逐渐减少。

Table 5. Soil erosion area proportions at different elevations (unit: %)

表5. 不同海拔高程带上的土壤侵蚀面积占总侵蚀面积比例(单位：%)

由研究区DEM数据生成坡度图，计算在6个坡度区间(0˚~5˚，15˚~25˚，15˚~25˚，25˚~35˚，35˚~45˚，>45˚)的土壤侵蚀面积分布情况(表6)。0˚~5˚坡度区间的土壤侵蚀面积最大，占土壤侵蚀总面积的60.24%，其次是5˚~15˚坡度区间，其土壤侵蚀面积占土壤侵蚀总面积的26.58%，二者共占土壤侵蚀总面积的86.82%。坡度在35˚以上区间土壤侵蚀面积很小。中度、强度和极强度的土壤侵蚀主要发生在15˚~25˚坡度区间。

Table 6. The proportions of soil erosion area of different slopes to the total soil erosion area (unit: %)

表6. 不同坡度等级中的土壤侵蚀面积占总侵蚀面积比例(单位：%)

5. 结论

本研究采用RS和GIS技术，基于RUSLE 土壤侵蚀模型，估算并分析淮河上游息县流域的土壤侵蚀以及空间分布规律。主要结论如下：

1) 2000~2015年的降水侵蚀因子呈下降趋势，主要受到流域降雨量和降雨强度的影响。植被覆盖管理因子和水土保持措施因子近15年来基本呈现为减小的趋势。

2) 2000~2015年流域土壤侵蚀模数总体呈减小趋势。土壤侵蚀强度主要以微弱侵蚀为主，占总面积的92.49%。西部桐柏山区和南部大别山区土壤侵蚀较为严重。

3) 海拔、坡度等地形因子对土壤侵蚀空间分布有重要影响。高程在300 m以下的土壤侵蚀面积占总面积的90.03%。坡度在15˚以下的土壤侵蚀面积占总面积的86.82%。中度、强度和极强度的土壤侵蚀主要发生在100~300 m高程区间，以及15˚~25˚坡度区间。

基金项目

江苏省高等学校大学生创新创业训练计划项目(201810323090X，201810323034Y，201810323035Y)；国家自然科学基金(41801075)；中国科学院流域地理学重点实验室开放基金(WSGS2017009)；江苏省高等学校自然科学研究项目(18KJB170002)；中国博士后科学基金(2018M642349)。