基于RUSLE 和GIS 的博罗县土壤侵蚀特征分析

陈燕奎，张 豫，朱猷亮

（1.嘉应学院 地理科学与旅游学院，广东 梅州 514015;2.深圳航天宏图信息技术有限公司,广东 深圳 518000）

为了全面认识土壤侵蚀的过程，理解全球变化与区域土壤侵蚀关系[1]，Wischmeier 和Smish 分析了全美各地径流泥沙观测资料，提出通用土壤流失方程（Universal Soil Loss Equation,USLE）[2].1985 年美国科学家对USLE 进行修正，研发出修正通用土壤流失方程（Revised Universal Soil Loss Equation,RUSLE）[3].由文献可知，学者利用RUSLE 方程研究水土流失土壤侵蚀量估算是有意义的，主要围绕RUSLE 方程的研究进程主要体现在侧重土壤侵蚀量估算[4-6]、方程因子改进与增加方程参数估算[7]与GIS、RS 等技术结合做土壤侵蚀的时空分析[8-9]以及在全球不同地区、不同尺度、不同时序等方面深度剖析土壤侵蚀方程的应用效益[10-12]；区域研究多以省、市层面，而南方山区县级层面的研究还比较少.本文以广东省博罗县域为例，引入基于修正的通用土壤流失方程（RUSLE），综合运用GIS 和RS 技术，试探索定量评价博罗县流域土壤侵蚀现状及空间分布，并分别估算不同坡度及土地利用类型下的土壤侵蚀量，为本区域的生态环境保护和政府决策提供参考依据.

1 区域概况与数据源

1.1 区域概况

博罗县位于广东省珠江三角洲东北部，地理坐标为北纬23°03′50″～23°43′20″，东经113°49′50″～114°45′50″，全县总面积约2 858 平方千米，属东江水系，地貌类型为北部山地丘陵、间有山谷平原，中部丘陵台地，南部沿东江自东向西的冲积平原等3 个地带，海拔约为208～368 m，属亚热带季风气候，年平均气温22.1℃，多年平均降水量约为1 918 mm，日照时数约为1 871.5 h.全域地形起伏比较大，海拔-10～1 239 m，土壤则以粉砂土为主，粉质黏粘土为辅，另外森林覆盖率较高，土地利用类型多样，植被物种丰富.

1.2 数据源

所用数据来源如表1 所示，在ArcMap10.7 和PIE-Basic6.0 等软件功能的支持下，建立博罗县域的土壤流失方程各类因子空间数据库，并将所有因子图层设立统一坐标系（CGCS2000 坐标系）和划分评价单元：30 m×30 m 的栅格单元.此外，土地利用数据类型按照用地性质初步划分为耕地、有林地、水域、建设用地、草地、裸地以及未利用地7 大类等.

表1 博罗县土壤侵蚀所需基础数据

2 数据处理与方法

2.1 土壤侵蚀方程

基于国内外现有成果和博罗县土壤侵蚀发生的特点，建立了坡面土壤侵蚀模型，即博罗土壤流失方程BLSLE（BoLuo Soil Loss Equation），其基本形式为

式中：A为土壤流失量，t·(hm-2·a-1)；R为降雨侵蚀力因子，MJ·mm·(hm-2·a-1)，1MJ=106J；K为土壤可蚀性因子，t·hm2·h/(hm2·MJ·mm)；L为坡长因子，无量纲；S为坡度因子，无量纲；C为地表植被覆盖与管理因子，无量纲；P为水土保持措施因子，无量纲.

根据式（1），博罗县土壤流失方程构建关键在于6 个区域性参数取值的确定.

2.2 RUSLE 模型中各因子值的确定

2.2.1 降雨侵蚀力因子（R）

在USLE 方程中，降雨侵蚀力因子（R）被定义为降雨动能和最大30 分钟降雨强度的乘积.引用文献[13]、[14]基于日降雨量提出计算广东地区降雨侵蚀力的简易算法模型：

式中，iR表示第i个半月时段的侵蚀力值（MJ·mm·hm-2·-h·-a）；k表示该半月时段内的天数；pj表示半月时段内第j天的日雨量，要求日雨量≥12 mm，否则以0 计算，12 mm 与侵蚀性降雨标准对应；α和β为模型待定参数.模型参数α和β值的计算方程为

式中，Pd12表示日降雨量≥12 mm 的日平均雨量（mm），Pd12表示日降雨量≥12 mm 的年平均雨量.

文献[13]指出，在年降雨量较小的地区，模型的决定系数较小、估算多年平均侵蚀力的相对误差较大，而在年降雨量比较丰富的地区相反，模型的决定系数较大、估算多年平均侵蚀力的相对误差较小.而本研究区年降雨量比较丰富，20 年日降雨量≧12 mm 的年平均雨量为1 615 mm，所以该公式在本地区可行.收集到博罗县及其附近增城、龙门、河源、惠阳、惠东共6 个站点2000～2021 年的日降雨量数据计算博罗县降雨侵蚀力的空间分布，具体地理位置见图1.小部分缺失数据，通过最邻近法空间插值方法，保持数据序列的完整和连续性.

图1 气象站点位置

首先，根据式（3）、（4）和收集到的各气象站点日降雨数据，计算各气象站点的模型参数α和β值（见表2），再根据式（2）计算各气象站点半月时段的降雨侵蚀力值并累加得出各年降雨侵蚀力值，取各站点2000～2021 年平均降雨侵蚀力值为该站点的降雨侵蚀力值，应用样条函数法插值生成研究区降雨侵蚀力分布如图2 所示.

图2 博罗县降雨侵蚀力分布

表2 各站点模型参数与年平均降雨侵蚀力值

2.2.2 土壤可蚀性因子（K）

结合博罗县的具体调查确定土壤机械颗粒组成，选择Williams 等提出的EPIC 模型，计算土壤可蚀性因子的值.按USDA 土壤类型分类标准对1∶100 万土壤类型进行分类，将博罗县区域USDA 土壤分类为粉质粘土、粉质粘壤土、粉土、壤土和砂土，各种土壤在区域所占比例如图3 所示.

图3 博罗县区域USDA土壤分类

根据博罗县USDA 土壤类型分类图，以及各种土壤类型中粘粒、粉粒、砂粒和有机碳含量，运用侵蚀-生产力影响模型EPIC（Erosion Productivity Impact Calculation）中土壤可蚀性因子估算公式（5）计算各种土壤类型的土壤可蚀性[15]，土壤可蚀性值计算结果如表3 和图4 所示.

图4 土壤可蚀性值分布

图5 L·S 乘值的空间分布特征

表3 土壤组成及土壤可蚀性值

式中，KEPIC为EPIC 模型计算的土壤可蚀性因子（t·hm2·h·hm-2·MJ-1·mm-1），SN、SIL和CLA分别为美制土壤粒级分类标准中的砂粒、粉粒和黏粒含量，SN1=1-SN/100，Cα为有机碳含量（%）.

2.2.3 坡度因子（S）与坡长因子（L）

坡度和坡长均是土壤侵蚀量变化的加速因子，主要体现在坡度高低和坡面长度对土壤沙土侵蚀速度、搬运量和长度的影响.同等条件下，坡度越大，坡面的冲刷能力越强，而且搬运量将随着波长（L）的增大而增大.对于波长和坡度研究，小尺度研究可根据实测值代入分析，对于大区域的估算，可从DEM 数据提取坡度和波长.修正土壤流失方程（RULSE）的坡长坡度因子计算方法为：

式中，S为坡度因子；θ为DEM 提取的坡度值；L为坡长因子；λ为DEM 提取的坡长值：22.13 m 为标准小区的坡长；α为坡度坡长因子.

该模式是学者在美国耕地坡度现状基础上建立的，因各地区的坡度坡长变化不用，该模式并非适合直接套用到其他研究地区，只能作为基础模型范本，由此，实际使用该模式需要根据地方实际进行模型参数修正.

（1）坡度因子修正.美国耕地坡度普遍小于20%（11.3°），博罗县大部分耕地分布于山地丘陵的坡地上，坡度值大于15°，坡地利用在博罗县丘陵山地较为普遍，因此借鉴文献[16]对坡度在9%～55%的陡坡侵蚀的研究对坡度因子进行修正，坡度小于14°沿用原有的坡度因子计算方式，其余坡度值，重新提出坡度因子计算方式：

因此将修正后的坡度因子进行合并，则得到

（2）波长因子修正.坡长因子在ArcMap10.7 软件下进行修正计算，通过PIE-Basic6.0 等软件DEM 提取坡度数值图，可以根据坡度值推算出α坡度坡长因子，公式如下：

式中β为细沟侵蚀和面蚀的比值.

2.2.4 地表植被覆盖与管理因子（C）

植被覆盖率与土壤水土保持具有很强的相关性，一般而言植被覆盖度越大，其植被叶、树干对雨水冲刷地表起到降速与阻隔作用和植被根系对土壤起到固定作用就越大，从而使地表土壤在受到雨水作用时越不轻易被剥蚀.学者们对C值研究，定义了C值取值介于0～1 之间，植被覆盖度（ƒ）越高，其C取值越小，表示植被对土壤侵蚀的抑制作用越大.采用蔡崇法等[17]的分析方法，用于分析植被覆盖度（f）与对应C值的关系.

（1）计算植被覆盖度（f）.利用PIE-Basic6.0 软件将2019-2021 年的1 月，4 月，7 月，10 月的TM-L8 数据集，提取涵盖博罗县范围的影像数据，利用红外与近红外波段，运用ENVI 软件的波段运算功能计算博罗县植被归一化指数（NDVI）年度NDVI 平均值：

式中，f为植被覆盖度（%）；NDVI为归一化植被指数，其中NDVImax和NDVImin为归一化植被指数中的最大值和最小值.

（2）计算植被覆盖与管理因子（C），公式如下：

根据文献[17]，当植被覆盖度为0 时，最易发生土壤侵蚀，C取值为1，覆盖度大于78.3%时，发生土壤侵蚀概率小，C取值为0.在PIE-Basic6.0 软件的支持下，利用栅格计算器功能，得出研究区植被覆盖与管理因子空间分布图，如图6 所示.

图7 水土保持措施因子空间分布

2.2.5 水土保持措施因子（P）

水土保持措施因子（P）主要以采取水土保持措施后的土壤侵蚀量与对应的未采取水土保持措施的情况下顺坡种植时的土壤侵蚀量的比值.取值范围0～1 之间，丘陵山区的随着地表利用方式的不同，其P值也有差异，从而土壤侵蚀的作用也有不同，P值越小，表示水土保持能力越好，土壤越不易被侵蚀.参考文献[18]、[19]对不同的土地利用类型的P进行赋值，以博罗县2021 年的土地利用现状数据为基础（见表4），运用ArcMap 分析工具得到P 值栅格图层.

表4 不同土地利用类型P值

3 结果与分析

3.1 研究区土壤侵蚀空间分布特征

将上述各因子的初始数据进行坐标匹配，CGC2000 坐标系，划分栅格单元格大小30 m×30 m，利用ArcMap10.7 软件的栅格计算功能将博罗县RUSLE 模型各因子计算结果进行栅格叠加计算，得到博罗县的土壤侵蚀模数分布图（见图8），博罗县年平均土壤侵蚀模数为9721.24 t·km-2·a-1，年侵蚀量892.61×104t.

图8 博罗县土壤侵蚀模数空间分布

从空间分布特征分析，将博罗县土壤侵蚀分布按照镇级行政区划进行统计分析，得到的各行政区土壤侵蚀情况统计如表5 所示，按照南方红壤丘陵侵蚀等级划分，石湾镇及园洲镇属微度，土壤侵蚀较为严重的行政区排名前五的是罗浮山、横河镇、公庄镇、柏塘镇及泰美镇，土壤侵蚀总体发生区域分布在西北部及东中部山地丘陵区，该区域海拔相对较高，山地多，降雨充沛，可以看出雨水冲刷和地形起伏对土壤侵蚀的影响比较显著，为土壤侵蚀提供了有利条件.

表5 博罗县各行政单元土壤侵蚀量

博罗县属于南方红壤丘陵区，分为6 个侵蚀等级，其大部分地区年平均土壤侵蚀模数属于轻微度水平，按照该标准分类，得到博罗县各等级土壤侵蚀面积如表6 所示，各等级土壤侵蚀空间分布如图9 所示.各侵蚀强度的面积大小依次是微度＞剧烈＞轻度＞强度＞极强度＞中度，分别占47.45%、15.36%、14.28%、10.81%、7.77%、4.34%.虽然博罗县土壤侵蚀等级极强度以上比重只占23.13%，但土壤侵蚀量却占到82.51%，其中剧烈等级的平均侵蚀模数高达46 759.54 t·km-2·a-1，远高于当地容许土壤流失标准500 t·km-2·a-1，微弱侵蚀面积仅为47.45%，而剧烈侵蚀年土壤侵蚀量却为微度侵蚀年土壤侵蚀量的63 倍.

图9 博罗县土壤侵蚀分级

表6 博罗县土壤侵蚀分级

3.2 土壤侵蚀与坡度空间分布特征分析

地形是土壤侵蚀的主要影响因素，利用博罗县DEM 数据生成坡度图，按照水利部土壤侵蚀分类分级标准将坡度范围划分为6 个等级区，统计得出各等级下的土壤侵蚀分析如表7 所示.土壤侵蚀面积最大的为坡度0°～5°，占土壤侵蚀面积比例的51.12%，其土壤侵蚀量的比例仅占总量的10.24%，其平均土壤侵蚀模数为3 265.26 t/（km2·a），处于中度以下侵蚀级别，也是6 个坡度等级中侵蚀模数最低的.8°～15°和15°～25°两个坡度等级是仅次于0°～5°土壤侵蚀面积，分别占到土壤侵蚀面积的19.13%和15.56%；两者的土壤侵蚀量比例分别为25.11%和36.81%，为6 个坡度等级中侵蚀量比例最大的两个值，另外25°以上两个坡度等级，虽然其侵蚀面积较少，仅占土地总面积的2.82%，但其土壤侵蚀量却占到侵蚀总量的19%，而且是6 个坡度等级中平均侵蚀模数最大的两个，远超过南方红壤侵蚀最高等级的基数（15 000 t/（km2·a））.

表7 博罗县不同坡度下土壤侵蚀特征分析

3.3 土壤侵蚀与土地利用空间分布特征

利用土地利用类型图和土壤侵蚀强度等级图进行叠加分析分析如表8 所示.各土地利用类型下的土壤侵蚀量来看，有林地范围内的土壤侵蚀量占到总侵蚀量的45.08%，耕地和园地范围的土壤侵蚀量占到总侵蚀量的25.16%和13.43%，另外虽然裸地范围仅占土地总面积的0.45%，但其土壤侵蚀量却占到总土壤侵蚀量的6.41%，属于激烈侵蚀等级.

表8 博罗县不同土地利用下土壤侵蚀特征分析

3.4 植被覆盖度与土壤侵蚀的影响

利用植被覆盖度（C）和土壤侵蚀空间分布图层进行叠加分析，划分为5 个等级：＜20%、20%～35%、35%～50%、50%～65%、＞65%.统计各覆盖度等级下的土壤侵蚀情况如表9 所示.植被覆盖度小于20%和大于65%的区域土壤侵蚀面积高达26.78%和24.57%，其他植被覆盖度级别下的土壤侵蚀面积变化不大，植被覆盖度小于20%的平均侵蚀模数最大，为30 045.57 t/（km2·a），虽然植被覆盖度低的区域为城镇村建设用地和水域地区，但其还包含裸地、稀疏林等，其雨水冲刷和地形坡度的影响下，使得其平均侵蚀模数较大.而植被覆盖度高的地区，其平均侵蚀模数较小，均为2 240.67 t/（km2·a），侵蚀量也最小，仅为62.98 t/（a），虽然其地形起伏大，因其林木高而密，林木种质多样，林下层级多样，对雨水冲刷起到一定减缓作用.另外20%～65%三个覆盖度等级，其土壤侵蚀变化不大，但其侵蚀量则占到土壤侵蚀总量的81.23%，这些地区多为山区周围，地形起伏变化大，而且种植了大量的人工林、稀疏林木、次生林，特别是桉树类，林分结构较为单一，林下植被较为匮乏，在雨季，特别是大雨、暴风雨等，雨水的冲刷力大，植被对土壤保持起到的作用较弱，其土壤侵蚀情况不容忽视.

表9 博罗县不同被覆盖度下土壤侵蚀分析

4 结论

本文利用修正的通用土壤流失方程RUSLE 模型，在GIS 技术支持下，定量估算了博罗县土壤侵蚀强度，并结合遥感影像、DEM、土地利用等数据，通过空间分析和数理统计方法分析了博罗县内土壤侵蚀的空间分布特征及影响因素，得出如下结论.

（1）2021 年，博罗县土壤平均侵蚀模数为9 721.24 t·km-2·a-1，表明水土流失较为严重，年侵蚀量892.61×104t.6 个侵蚀强度等级的面积大小依次是微度＞剧烈＞轻度＞强度＞极强度＞中度，其中剧烈等级侵蚀面积虽然只占到15.36%，但土壤侵蚀量却占到总侵蚀量的54.7%.

（2）石湾镇及园洲镇属微度，而西北部及东中部山地丘陵区的罗浮山、横河镇、公庄镇、柏塘镇及泰美镇等区域，该区域海拔相对较高，山地多，地形起伏大、降雨充沛，雨水冲刷大，土壤侵蚀较为严重.

（3）博罗县土壤侵蚀量与坡度存在明显的正相关，每增加一个等级，平均侵蚀模数呈增大趋势，而且增大幅度较大；不同的土地利用类型会直接影响土壤侵蚀的量，有林地的面积为52.64%，其土壤侵蚀量比重占45.08%，而裸地虽然面积比重仅占0.45%，但其土壤侵蚀量却占到总土壤侵蚀量的6.41%，属于激烈侵蚀等级，因此裸土地是土壤侵蚀易发地区，需要进一步加强水土保持治理力度，做到治理与预防并重.

（4）植被覆盖度小于20%和大于65%的区域土壤侵蚀面积高达26.78%和24.57%，其他植被覆盖度级别下的土壤侵蚀面积变化不大；植被覆盖度20%～65%之间的土壤侵蚀面积为48.64%，但其土壤侵蚀总量比重达81.23%，该区域多为山区丘陵，地形起伏大，分布着大量的人工林、次生林等林分结构单一、林下植被较为匮乏的经济林木，雨水冲刷严重，不利于水土保持，该区域土壤侵蚀情况不容小觑.