基于随机森林算法的土壤侵蚀影响因子研究——以赣江上游流域为例

朱 青， 国佳欣， 郭 熙， 韩 逸， 刘士余， 陈 蕾， 罗海玲

(江西农业大学 国土资源与环境学院 江西省鄱阳湖流域农业资源与生态重点实验室， 江西 南昌 330045)

土壤侵蚀是指土壤及其母质在内营力和外营力的共同作用下，发生风化、剥蚀、搬运和堆积的过程[1]。土壤侵蚀不仅会引起土地退化、土壤肥力下降，还会造成水体富营养化、水库河道淤塞、旱涝灾害加剧等一系列生态环境问题[2-4]，因此科学认识和定量评估区域土壤侵蚀，分析土壤侵蚀空间分布特征及其影响因素对水土流失防治工作和生态环境恢复具有重要意义[5-6]。近年来，修正的通用土壤流失方程(revised universal soil loss equation, RUSLE)[7]在各区域尺度下的土壤侵蚀定量研究方面得到了广泛应用，该方程全面考虑了影响土壤侵蚀的因素[6,8]，具有较强的代表性和实用性[9]。众多学者将RUSLE模型与空间信息技术(RS和GIS)相结合对土壤侵蚀进行了评估与分析[2,6,10]，但目前对土壤侵蚀及其影响因素的研究多基于单一流域[5,8,11]，缺乏对影响流域内不同子流域土壤侵蚀因素的定量分析，这对关键区域水保措施的开展具有一定的局限性，无法提供针对性的措施建议。因此，准确掌握影响子流域土壤侵蚀的主要因素，对于提升水土保持工作的精准化具有现实指导意义。此外，前人研究对土壤侵蚀影响因素的探究多停留在区间统计分析[6,12]、相关性分析[13]和经典回归模型[8,14]等传统分析方法，难以揭示因子间的复杂过程和非线性关系[15]。而随机森林是一种基于统计学的机器学习算法[16]，可以有效地解决过度拟合、多元共线性等问题，结果具有可解释性，能够很好地量化自变量对因变量的重要程度，有效地弥补了传统统计分析方法的不足[15]，在各行各业中均有广泛应用[17]。

赣江上游流域为鄱阳湖流域最大支流—赣江的水源区，该地区的土壤侵蚀不仅会影响到中下游流域的饮水质量和粮食安全，还会影响到鄱阳湖入湖泥沙量，进而造成一系列生态环境问题，因此该区域生态环境状况是维系鄱阳湖一湖清水，促进鄱阳湖生态经济区和绿色生态江西建设的重要保障[18]。但该区域生态环境较为脆弱，部分区域的土壤侵蚀在自然因素和人为因素的双重作用下越来越剧烈，水土流失已成为当地重要的生态环境问题之一[19]，引起了学者们的高度重视，诸多有关土壤侵蚀的研究在该地区相继开展[12,18]；但在量化影响因子与土壤侵蚀之间的研究尚不多见。因此，本文以赣江上游流域为研究区，基于2015年Landsat 8遥感影像、MODIS NDVI数据、数字高程模型(DEM)、土壤类型和雨量站点等数据，在RS和GIS技术的支持下，采用RUSLE模型对赣江上游流域土壤侵蚀进行定量分析；并将各子流域作为一个独立单元，识别研究区内的重点防治区域，同时借助随机森林算法系统地分析影响因子在不同子流域间的重要程度，旨在科学准确地界定影响因子与子流域土壤侵蚀的量化关系，并在划分土壤侵蚀重点治理区时，为决策者科学规划实施有针对性的水土保持措施提供参考依据。

1 研究区概况

以赣江上游流域为案例区，该区地处江西省南部，南岭之北，位于113°54′—116°38′E，24°29′—27°09′N。流域面积为32 948.72 km2，约占赣江流域面积的39.46%，占赣州市全域面积的83.67%，涉及章贡区、赣县区、兴国县等16个县区。该区属亚热带湿润季风气候，热量丰富、降水充沛，多年平均气温为18.9 ℃，1月蒸发量最小，7月蒸发量最大；多年平均降雨量为1 584.14 mm。流域内地势呈四周高中间低，地形复杂多样，多山地丘陵，坡度在5°～25°的区域占流域面积的64.40%；地貌类型分为西部中、低山构造侵蚀地貌，南部低山、丘陵构造剥蚀地貌，中部丘陵、河谷侵蚀堆积地貌和东北部低山、丘陵构造剥蚀地貌和溶蚀地貌。成土母质以易产生水土流失的花岗岩类、泥质岩类和酸性结晶岩类风化物为主。土壤类型以红壤为主，约占总面积的82%，可蚀性较高[4]，是典型的南方红壤丘陵分布区。流域内植被覆盖率约为64%，植被群落结构单一，多为马尾松和湿地松林，林下植被匮乏，林下流现象严重[20]。

2 研究数据与技术流程

2.1 数据来源

研究所用主要数据包括： ①降雨数据。赣江上游流域15个雨量站1981—2015年逐日降雨量数据，来源于中国气象科学数据共享服务网(http:∥data.cma.cn/)； ②土壤数据。世界土壤数据库(HWSD)中的中国土壤数据集与江西省第二次土壤普查数据集； ③高程数据。数字高程模型(DEM)数据，空间分辨率为30 m，来源于地理空间数据云(http:∥www.gscloud.cn/)； ④植被数据。由最大值合成的2015年MODIS MOD13 Q1 NDVI产品数据，空间分辨率为250 m，时间分辨率为16 d，来源于地理空间数据云(http:∥www.gscloud.cn/)； ⑤土地利用类型数据。选用2015年Landsat 8 OLI/TIRS 5景遥感影像(空间分辨率为30 m)，结合研究区2015年1∶10 000土地利用变更调查数据库，通过人机交互解译共生成8大地类：水田、旱地、草地、园地、林地、建设用地、未利用地和水域； ⑥2015年年均人口密度栅格数据。分辨率为1 km，来自于中国科学院资源环境数据中心(http:∥data.cma.cn/)； ⑦赣江上游流域及其子流域矢量边界数据。由Arc SWAT软件中的水文分析模块(watershed delineator)提取而来。

2.2 技术流程

研究主要分为两个流程：①在空间信息技术RS和GIS的支持下，利用RUSLE模型对赣江上游流域土壤侵蚀特征进行定量评估，得到研究区2015年土壤侵蚀强度空间分布图。为保证各模型因子计算的准确性，本文统一采用30 m×30 m栅格单元大小、投影坐标系统一为WGS_1984_UTM_Zone_50N； ②为减少数据量和运算量，选用最邻近法对土壤侵蚀栅格图重采样至250 m，将其与MODIS NDVI数据的空间分辨率保持一致。在R软件中调用随机森林算法程度包RandomForest对赣江上游子流域中的影响因子(R，K，LS，C，P)进行重要性分析，以确定影响流域中土壤侵蚀的关键因子。

3 研究方法

3.1 RUSLE模型

本文选用修正的通用土壤流失方程(RUSLE)[7]对赣江上游流域土壤侵蚀量进行估算，数学表达式为：

A=R·K·L·S·C·P

(1)

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

3.1.1 降雨侵蚀力因子(R)R是指降雨过程对土壤流失存在的潜在影响力[9]。本文选用Renard等[21]提出的降雨侵蚀力算法模型，将赣江上游流域15个雨量站点数据运用如下公式进行计算：

(2)

式中：Pa为年均降雨量(mm)。将研究区15个雨量站的降雨数据计算得到的R值(表1)进行克里金插值[22]。

表1 赣江上游流域雨量站及其降雨侵蚀力

3.1.2 土壤可蚀性因子(K) 土壤可蚀性因子K是表征土壤性质对侵蚀敏感程度的指标[7]。本文应用侵蚀—生产力评价模型(EPIC)[23]进行估算，计算公式为：

(3)

式中：K为土壤可蚀性因子〔t·hm2·h/(hm2·MJ·mm)〕；Wd为砂粒含量(%)；Wi为粉粒含量(%)；Wt为黏粒含量(%)；Wc为有机碳含量(%)。

3.1.3 坡长坡度因子(LS)LS是反映地形地貌特征对土壤侵蚀影响的指标[7]。

考虑到研究区坡度大于10°的区域占总面积的59.30%，本文采用RUSLE模型计算的坡度因子公式不适宜南方丘陵区实际情况，因此本文借鉴Liu等[24]和McCool等[25]提出的坡度坡长计算方法对LS因子进行提取：

(4)

式中：S为坡度因子；θ为坡度(°)；L为坡长因子；λ为水平投影坡长(m)； 22.13为标准小区坡长；m为可变坡长指数，即坡面角度的正切值。

3.1.4 植被覆盖与管理因子(C)C为侵蚀动力的抑制因子，对水土保持具有积极作用[7]，其与土地利用类型、植被覆盖度密切相关[26]。根据野外考察发现，研究区林下流现象较为严重，主要发生在海拔500 m以下的林地类型中，且在海拔200 m以下大多为马尾松和湿地松林，林下流普遍存在。因此，本文在咨询相关专家意见和研究报道[10,22]的基础上，结合植被覆盖度和海拔的影响，对不同的土地利用类型赋予不同的C值(表2)。植被覆盖度fc计算公式为：

(5)

式中：fc为植被覆盖度，fc值介于0～1之间，fc越趋近于1，植被覆盖度越高。本研究以2%和98%的累计百分比为置信度区间，读取对应的像元值为研究区内NDVI的最小值和最大值。

表2 赣江上游流域植被与覆盖管理因子C取值

3.1.5 水土保持措施因子(P)P是指实施水土保持措施下土壤的流失量与顺坡种植时的土壤流失量之比[7]，范围在0～1之间，p=0表示在水土保持措施后不发生侵蚀，P=1表示未采取相应水土保持措施。赣江上游流域为我国重点水土流失综合治理区，全域范围内开展了多种水保措施，主要集中在经果林开发、基本农田改造等。鉴于这种情况，本文在参考相关学者研究成果[4-5,12]的基础上，结合流域内水土保持措施在不同生态系统的强弱程度，对不同土地利用类型的P因子进行赋值，结果详见表3。

表3 赣江上游流域水土保持措施因子P取值

3.2 基于随机森林的重要性分析

随机森林(random forest, RF)是一种基于分类树的机器学习算法[16]，其对多元共线性不敏感，无需对数据进行归一化处理，结果对缺失和非平衡的数据有较好的容忍度。变量的重要性评估是RF算法的一个重要特点，它通过利用Bootstrap方法，从原始数据集中进行有放回的随机抽取，形成训练数据集。在使用Bootstrap对原始数据集进行抽取时，每个数据未被抽取的概率为(1-N-1)N，其中N为原始数据集中数据的总量。这些未被抽中的数据称为袋外数据。模型主要是利用袋外数据进行估计，其随自变量的变换而产生预测误差，RF算法可以根据误差的变化来计算自变量的重要性，重要性值越大，表明对土壤侵蚀的影响越大[15,17]。本研究基于R软件中的Random Forest程序包，采用回归算法来计算各影响因子的重要性评分，其中决策树的数量与分割节点的个数分别设置为500，2。

4 结果与分析

4.1 各因子特征

1981—2015年平均降雨侵蚀力因子R为9 092.31 (MJ·mm)/(hm2·h·a)，变幅在7 569.60～11 621.35 (MJ·mm)/(hm2·h·a)之间，标准差为1 036.69 (MJ·mm)/(hm2·h·a)，这与马良等[27]对江西省约50 a来降雨侵蚀力变化范围研究结果基本一致。由图1可以看出，流域内降雨侵蚀力总体呈四周高，中间低的分布趋势。其中，降雨侵蚀力的高值区主要集中分布在东北部的宁都县、石城县气象站附近，低值区位于北部的赣县区气象站附近。土壤可蚀性K值范围在0～0.02 (t·hm2·h)/(hm2·MJ·mm)之间，均值为0.018 (t·hm2·h)/(hm2·MJ·mm)，标准差为0.003 (t·hm2·h)/(hm2·MJ·mm)。

注:地貌类型分图中，1为低海拔平原; 2为低海拔台地; 3为低海拔丘陵; 4为小起伏低山; 5为小起伏中山; 6为中海拔丘陵; 7为中起伏低山; 8为中起伏中山。

图1 赣江上游流域因子空间分布

从图1可以看出，土壤可蚀性高值区域主要呈条带状和块状分布在中部低海拔台地和西部的中海拔丘陵处，而海拔相对较高处的山地、丘陵土壤可蚀性较小，朱成刚等[28]对犁河谷土壤理化性质及可蚀性特征研究时也得到类似结果。研究区LS因子范围在0.038～48.98之间，均值为5.51，标准差为4.60，该结果与李新艳等[29]对赣南地区LS因子计算范围相一致。由图1可知，在赣江上游流域不同地貌类型中，低海拔平原的LS因子最小，均值为0.66；中起伏中山最大，均值为14.91。整个流域四周的LS因子较大，中部较小，且北部总体小于南部，这与赣南四周高山环绕，中部丘陵起伏，平原和小盆地散布、河流水系呈辐辏状向中心—章贡区汇集的地形特征相吻合(图1)。赣江上游流域C因子均值为0.018，标准差为0.043。从图1可以看出，低值区分布较为广泛，这与该区域内植被生长状况良好，植被覆盖度总体达63.54%有关。水土保持措施因子P主要依据不同土地利用类型赋值而来，均值为0.79，标准差为0.33。高值区主要呈片状分布在林地、草地类型中；低值区则主要呈块状分布在建设用地和水田类型中。

4.2 模型结果验证

本文将RUSLE模型估算出的2015年研究区土壤侵蚀总量与2001年江西省水利厅第三次土壤侵蚀遥感调查数据结果[30]进行对比。该模型估算出的赣江上游流域土壤侵蚀总量为3.45×107t/a；调查结果记录赣州市的土壤侵蚀总量为4.87×107t/a，且研究区约占赣州市的83.67%，因此RUSLE模型估算出的结果按比例换算约为公报的84.52%(相对误差在20%以内)，结果在合理范围内，具有较高的可信度。另外，李恒凯等[12]利用RUSLE模型对2013年赣州市土壤侵蚀进行估算时得出轻度侵蚀以上所占面积为16 177.17 km2，与本文赣江上游流域轻度侵蚀以上面积为15 399.70 km2相接近。综上可知，该模型的土壤侵蚀估算结果较为准确。

4.3 流域土壤侵蚀总体特征

根据水利部最新颁布的南方红壤土壤侵蚀的技术标准(SL657-2014)[31]，将研究区土壤侵蚀强度划分为微度〔<500 t/(km2·a)〕、轻度〔500～1 500 t/(km2·a)〕、中度〔1 500～3 000 t/(km2·a)〕、强烈〔3 000～5 000 t/(km2·a)〕、极强烈〔5 000～10 000 t/(km2·a)〕和剧烈〔>10 000 t/(km2·a)〕6种等级。从图2可以看出，赣江上游流域主要以微度和轻度侵蚀为主，空间分布格局明显，侵蚀强度由东南向西北逐渐加剧；局部地区出现强烈以上的侵蚀强度，主要分布在研究区西北部等海拔相对较低与植被覆盖度低的低海拔台地、丘陵地区；而微度侵蚀的区域则分散在植被覆盖度高的小起伏低山和低海拔平原等地区。

图2 赣江上游流域土壤侵蚀强度空间分布

从表4可以看出，流域微度侵蚀面积所占比重最大，占流域面积的53.26%；中度侵蚀的面积为3 943.95 km2，占11.97%；强烈以上的侵蚀面积占总面积的10.19%，占比相对较高，因此须重点加强该区域土壤侵蚀的防治工作，以有效减少水土流失。从侵蚀量上来看，流域内以中度和强烈侵蚀为主，两者所占比重相当，各占侵蚀总量的25.00%，24.27%；剧烈侵蚀占比最小，仅为3.69%。整体而言，赣江上游流域土壤侵蚀总量为3.45×107t/a，平均土壤侵蚀模数为1 046.38 t/(km2·a)，低于陆建忠等[22]观测到的2000年赣江流域平均土壤侵蚀模数〔1 321.09 t/(km2·a)〕，表明近15 a来研究区土壤侵蚀状况呈现好转趋势，但总体上仍远大于南方红壤丘陵区土壤允许流失量500 t/(km2·a)的标准，因此流域内水土流失治理任务依然艰巨。

表4 赣江上游流域土壤侵蚀统计分析

4.4 子流域土壤侵蚀特征

本文利用Arc SWAT软件，根据研究区地形特征进行子流域的划分，经过水文分析模块计算后共得到26个子流域(图3)。从图4可以看出，各子流域平均土壤侵蚀模数范围在645.59～1 715.83 t/(km2·a)之间，流域间的差异较大。其中子流域9，11，15处于中度侵蚀级别，平均土壤侵蚀模数分别为1 672.66，1 715.83，1 565.36 t/(km2·a)；其余子流域均为轻度侵蚀，子流域14最小，为645.59 t/(km2·a)。从土壤侵蚀量来看(图4)，子流域15最大，其数值达到了4.18×106t/a，其次为子流域20；最小的是子流域9，为8.38×104t/a。

图3 研究区子流域空间分布

综上可知，子流域9，11，15是赣江上游流域水土流失的关键区域，为重点防治区。其中子流域9，11位于低缓河谷区，在流域总出水口附近；子流域15位于研究区西部低海拔丘陵地带。由表5可以看出，3个关键区域的红壤类型占比均达到了70%，降雨量均在1 400 mm以上。且从子流域15来看，其现有土地利用类型主要是以植被覆盖度较高的林地为主，占比接近90%，坡度坡长因子较大，分别为18.17°和87.98 m，高程为506.33 m。与之相比，子流域9，11建设用地面积占比相对较大，人口密度分别为256.16和636.55人/km2，林地、园地受人类活动干扰严重，植被覆盖度较低，分别为37.48%和40.28%；坡度均在15°以下，坡长适中，高程分别为149.92和201.30 m。

图4 子流域土壤侵蚀统计

表5 重点防治区域地理环境因素状况

4.5 影响土壤侵蚀因子的重要性分析

为进一步厘清土壤侵蚀影响因子的重要性，本文在子流域划分的基础上，通过250 m网格化采样后，在R软件中采用随机森林中的回归算法对R，K，LS，C，P因子进行重要性评分，然后将其换算为相对重要程度。从图5可以看出，各子流域回归系数R2均在0.7以上，拟合效果较好；子流域间的土壤侵蚀受C因子和LS因子影响较大，这与前人研究结果相一致[2,14,32]，且两者重要程度分别在30%和20%以上，其中C因子重要程度最高可达48.54%，表明C因子在很大程度上影响着该区域的土壤侵蚀，体现了流域内整体的相似性[11]；但在子流域9，11，21中，LS因子重要程度比C因子更大，说明在这三个子流域中，地形因素对土壤侵蚀起着主导作用；P因子的重要程度均在10%以上；R因子和K因子对土壤侵蚀的重要程度偏低，均未超过10%，表明这两种因子在整体上的差异性相对较小，各子流域间的变化起伏不大。此外，在重点防治区域中，C因子对子流域15土壤侵蚀的重要程度达到了42.88%，其次为LS因子，重要程度为21.88%；P和K因子重要程度分别为18.87%，7.2%；R因子最低，仅为6.04%。而在子流域9，11中，LS因子重要程度分别为36.15%，38.15%；C因子次之，为32.86%和28.00%；R因子最低，分别为2.83%，5.52%。

图5 子流域土壤侵蚀影响因子重要性分析

5 讨 论

赣江上游流域土壤侵蚀主要受到C因子和LS因子影响，且C因子对子流域土壤侵蚀的重要程度较大，均在30%以上。这与Kosmas C等[33]研究结果相一致，其认为植被覆盖和土地利用是最重要的影响因素，并且在一定程度上超过了降雨强度和坡度的影响。但胡刚等[14]在基于RUSLE模型估算卧虎山水库子流域土壤侵蚀时发现LS因子为主要因素，究其原因，本文研究区地处南方红壤丘陵区，林下流现象较为严重[20]，因此C因子在赋值时兼顾到了海拔的影响，涵盖了部分地形因素。但C因子在各子流域中的重要程度具有差异性，这主要是因为不同土地利用类型间的差异性较大[32]，C因子赋值相差甚远，且植被覆盖度在不同土地利用类型上的变化对土壤侵蚀具有放大效应[34]。

由重点防治区域结果分析可知(表5)，其红壤类型占比均达到了70%以上，该土壤土层薄、酸性强、黏重板结、有机质含量低，保肥保水能力较差[4,19]，受到侵蚀的敏感程度高；在降雨量较大的情况下，容易引起土壤发生分离和搬运，对土壤侵蚀具有潜在影响[9]；且从子流域15来看，虽然区域内植被覆盖度较高，达到63.81%，但植被群落结构单一，林下植被匮乏，“远看青山在，近看水土流”是研究区林下流的真实写照[20]，尤其是海拔500 m以下的林下水土流失较为严重，因此在对C因子赋值时未将植被覆盖度考虑在内；从图5可以看出，C因子对子流域15土壤侵蚀的重要程度最高，为42.88%，LS因子次之。因此，针对子流域15可采取增加林下覆盖度，形成层次分明的针阔混交林地，以减少地表径流量的生物措施为主，以改变坡长坡度、控制径流的工程措施为辅。与之相比，子流域9，11则位于低缓河谷区，受到人类干扰活动较大。陈思旭等[3]研究发现在南方丘陵山区，土地利用率较高的地方多集中在缓坡和海拔小于500 m的人口密集区，其多以种植柑橘，油茶和桉树等经济作物为主，植被覆盖度相对较低，本文研究结果与此相保持一致。LS因子对子流域9，11的土壤侵蚀重要程度分别为36.15%，38.15%，C因子次之。因此在制定水保措施时，一方面应考虑实施合理的工程措施进行微地形生态修复；另一方面可以改善土地利用类型的使用情况，注重植被覆盖度的提高，防止人类活动对生态系统造成进一步破坏。

本研究基于RUSLE模型，分析了国内水土流失重点区—赣江上游流域2015年土壤侵蚀的空间分布特征，对影响土壤侵蚀因子(R，K，LS，C，P)的重要性进行了定量化分析，标识了易发生水土流失的关键子流域，明确了各影响因子在流域之间的重要程度，这对于选择合适的水土保持措施来缓解关键区域的土壤侵蚀状况具有重要意义。本文下一步将缩小研究尺度，着重开展重点防治区域(子流域9，11，15)土壤侵蚀的内在机理研究，掌握小尺度下的径流输沙机制，以更好地完善子流域的水土保持工作。但本文还存在一定的不足之处，如在结果验证方面，只是借鉴已有的相关文献做参考，后期会在流域内建立径流试验小区进行模型参数的修正，使计算结果更为精确。另外，评估模型中的因子数据精度不一，因此在栅格计算过程中难免会产生一定的误差；且未能考虑林地类型及其垂直分层结构对土壤侵蚀的影响，而是在实地考察和询问相关专家意见的基础上，以海拔作为C因子对林下流赋值的依据。

6 结 论

赣江上游流域总体处于轻度侵蚀水平，空间分布格局明显，侵蚀强度由东南向西北逐渐加剧；各子流域平均土壤侵蚀模数范围在645.59～1 715.83 t/(km2·a)之间，流域间的差异较大，其中子流域9，11，15是水土流失的重点防治区域；C因子和LS因子是各流域内土壤侵蚀的主要因素，R因子和K因子的重要程度偏低，均未超过10%。研究考虑了不同子流域的土壤侵蚀情况，明确了各子流域科学管理与治理的优先顺序及其主控因子，可为决策者提供有针对性的水保措施建议。