基于多源环境变量和随机森林模型的江西省耕地土壤pH值空间预测

钟骁勇, 李洪义, 郭冬艳, 谢模典, 赵婉如, 胡碧峰

(1.江西财经大学财税与公共管理学院,南昌 330013; 2.中国自然资源经济研究院,北京 101149; 3.江西财经大学旅游与城市管理学院,南昌 330013; 4.江西财经大学财经数据科学重点实验室,南昌 330013)

0 引言

耕地资源是人类赖以生存的重要物质基础和生产资料,土壤pH值作为耕地土壤关键属性之一,不仅是影响农作物生长发育的重要因素,也是影响土壤养分有效性和土壤重金属元素迁移转换的关键因子[1]。但是,随着人类社会工业化进程加快,酸性气体(SO2)、颗粒污染物排放以及氮肥的大量施用等,使得土壤酸化在许多地方成为一个普遍的趋势[2-4]。张福锁院士团队在Science上发表的研究结果表明,从20世纪80年代—21世纪初,中国主要农作物产区土壤pH值明显下降[5]。土壤酸化会增加土壤重金属的活性,不利于农作物健康生长,进而威胁中国的粮食生产安全,这一问题在中国南方地区尤为突出[6-7]。因此,获得准确的土壤pH值空间分布信息对于遏制土壤酸化趋势、保护土壤环境质量和保障粮食生产安全具有十分重要的意义和价值。

国内外许多研究学者使用地统计学等空间插值方法用于土壤属性空间信息的预测,主要有普通克里格插值法(ordinary Kriging,OK)、反距离权重插值法(inverse distance weighted,IDW)、普通最小二乘法(ordinary least squares,OLS)、地理加权回归(geographically weighted regression,GWR)等[6,8-12]。近年来,随着计算机技术、卫星遥感和土壤近地传感器等技术的快速发展,机器学习算法应用于数字土壤制图(digital soil mapping,DSM)成为一个新的研究热点[13-21]。张甘霖等[22]指出: DSM未来有从地统计学走向机器学习为主导模型的趋势。随机森林(random forest,RF)作为一种以数据为驱动的机器学习方法[23],具有训练速度快、对多元线性不敏感、易于处理缺失数据等优点,还可避免对模型的过度拟合[24-26]。RF能够获取和表达土壤属性空间变化与环境变量的复杂非线性关系,并根据这种关系推测土壤属性空间分布,还能以相对较少的样本量取得较好的预测精度[27]。目前已有研究将RF模型用于土壤有机碳、土壤类型和土壤质地等属性空间预测和制图[28-31,26]。但是目前将RF应用于土壤pH值空间预测制图的研究还较为少见,尤其是国内运用RF模型开展大尺度耕地pH值空间预测分析的研究鲜有报道。

因此,本文以土壤酸化现象显著的江西省为例,使用RF算法和环境协变量信息开展耕地土壤pH值的预测研究。同时将预测结果与使用OK方法插值结果进行对比,验证使用RF算法开展大尺度耕地土壤pH值预测研究的可行性和可靠性,为江西省土壤酸化问题的治理和农业生产管理提供数据基础和依据。

1 研究区概况及数据源

1.1 研究区概况

江西省位于中国东南部,长江中下游南岸,位于E113°34′36″～118°28′58″,N24°29′14″～30°04′41″之间(图1),属于亚热带季风性湿润气候,是我国重要的粮食、油料、棉花、蔬菜和水果生产基地。江西地貌类型以山地、丘陵为主,山地占全省面积的36%,丘陵占42%,平原占12%,水域占10%。主要土壤类型为偏酸性的红壤,面积达93 111.32 km2,约占全省总面积的56%,具有较为明显的地带性和地域分布规律。

图1 研究区土壤采样点分布

1.2 数据源及预处理

采集土壤样点16 582个,同时选取样点的土壤属性以及植被指数、地形变量、气候变量和耕地利用条件作为耕地土壤pH值空间预测的环境变量。其中土壤属性包括土壤类型、有机质、有效磷、速效钾、全氮、全磷、全钾、阳离子交换量、成土母质和耕层质地,数据来源于2018年江西省耕地质量等别数据库、1980年和2010年江西省农业测土配方项目数据。归一化植被指数、地形变量的地貌类型和高程以及气象因子的年均气温和年降水量均通过资源环境科学数据中心(http: //www.resdc.cn/)网站下载获取,影像空间分辨率均为30 m。为了更好地反映出江西省地形条件对土壤属性变化的影响,选择高程、坡度、坡向、剖面曲率、平面曲率、沟谷深度、地形湿度指数、垂直到沟谷距离和多分辨率谷底平坦度等地形因子,分别借助ArcGIS 10.2和SAGA 7.6.2软件从数字高程模型(digital elevation model,DEM)数据中提取。耕地利用条件包括灌溉保证率、排水条件、氮肥用量、磷肥P2O5用量、钾肥K2O用量、秸秆还田比例和秸秆还田量,数据来源于2018年江西省耕地质量等别数据库(表1)。

表1 土壤pH值空间预测的环境变量及数据来源

2 研究方法

2.1 环境变量的组合和筛选

为全面分析土壤pH值空间变异规律以及探索不同环境变量对模型预测精度的贡献程度,将表1中30个环境变量分成2种不同的组合: ①以地形变量、气候变量和植被指数等遥感数据为基础,增加土壤属性和耕地利用条件等实测数据进行模型训练(RF-A); ②采用传统的地形变量、气候变量和植被指数等遥感数据训练模型(RF-B)。

2.2 RF模型构建

基于采集到的16 582个土壤样点数据,按照8∶2,7∶3和6∶4的比例分为训练集和验证集。RF模型采用bootstrap方法对样本进行放回抽样,没有被抽取的样本会自动生成一个对照集,因此不需要进行交叉验证。RF模型中2个可调参数决策树数量ntree和节点分裂次数mtry决定了模型的配置。

2.3 精度分析与比较

RF模型精度评价选用皮尔逊相关系数r、平均误差ME、绝对平均误差MAE和均方根误差RMSE等4个统计指标,其计算公式分别为:

,

(1)

,

(2)

,

(3)

,

(4)

3 结果与分析

3.1 土壤pH值的描述性统计分析

通过江西省耕地土壤pH值描述性统计结果可知,土壤pH值范围在4.00～8.30之间,平均值为5.20,中位数为5.10,标准差为4.91,变异系数为9.43%,表明江西省耕地土壤pH值呈弱变异。参考通用标准将江西省土壤pH值划分为7个等级: 极强酸性(pH<4.5)、强酸性[4.5,5.5)、酸性[5.5,6.5)、中性[6.5,7.5)、碱性[7.5,8.5)、强碱性[8.5,9.5)和极强碱性(pH≥9.5)。结果显示,江西省耕地土壤pH值样点中,极强酸性样本数67个,占比0.43%; 强酸性样本数12 718个,占比76.69%; 酸性样本数3 301个,占比19.90%; 中性样本数435个,占比2.62%; 碱性样本数61个,占比0.36%; 没有强碱性样本。综上可知,江西省耕地土壤总体呈现明显的酸性特征。

3.2 RF模型最优参数分析

通过逐次试验,确定RF模型中mtry和ntree参数的最优值。固定mtry(分别设为2,4,6和8),逐次调整ntree(分别设为300,600,900,1 200和1 500),进行4组20次试验。为避免过拟合问题,通过比较训练集和验证集的相关系数(r值),根据两者最为接近的结果确定最优预测模型的参数。结果表明(表2): 当mtry=8,ntree=900时,pH值预测模型的训练集和验证集r值最为接近,表明此时的模型稳定性最佳。

表2 RF模型中节点分裂次数和决策树数量的筛选

3.3 不同训练集的RF和OK模型预测精度比较

基于不同训练集运行模型所获得的预测结果有所不同(表3)。从相关系数r看,不同训练集的RF-A和RF-B模型预测的土壤pH值与实测值的r较高(均大于0.96)且差别不大,而OK模型r处于0.65左右,预测精度明显低于RF模型; 从误差看,RF-A模型的ME,MAE和RMSE值最小,而OK模型的误差明显更高。3种模型的拟合能力随训练样本量呈现不同的变化趋势。总体上看,RF-A模型略微优于RF-B模型,OK模型明显弱于前2个模型,并且训练集和验证集抽样比在7∶3时RF-A模型精度最佳。在验证数据集中,除了抽样比7∶3时,OK模型比RF-A模型精度高外,其他情况下RF-A模型的r更高。这表明环境变量增加耕地土壤实测数据有利于提高模型预测精度。因此在土壤样点数据充足的情况下,应优先选择RF-A模型用于土壤pH值的建模预测制图。

表3 使用不同训练集时RF和OK模型预测精度比较

3.4 环境变量重要性分析

使用RF模型定量计算出不同变量对于土壤pH值预测的重要性(图2)。由图2可知,气候因素如年均气温和年降水量是决定土壤pH值的最重要因素,这一结果与Chen等[32]、卢宏亮等[33]和王世航等[34]的研究结果基本一致。此外,地形地貌因子如沟谷深度、高程和多分辨率谷底平坦度也是决定土壤pH值的重要因子。说明江西省土壤pH值分布与地形地貌特征高度相关。土壤理化性质对于土壤pH值也有着不可忽视的影响。人为因素如化肥施用量、排水和灌溉条件等也对土壤pH值有较大影响。但是部分因子如土壤类型、秸秆还田比例的重要性低于预期,需要后续进一步研究加以分析。

图2 RF-A模型变量相对重要性

3.5 耕地土壤pH值空间预测制图

研究区耕地土壤pH值空间分布如图3所示,整体上使用OK模型、RF-A和RF-B模型获得的研究区土壤pH值空间分布特征大体相同,其中土壤pH低值区主要分布在江西省中部鄱阳湖平原的南昌市、抚州市和吉安市等地区,高值区域主要分布于江西省西部和东部山地丘陵区域,如九江市、鹰潭市、萍乡市和上饶市等地区。结合图1高程信息可知,这一分布特征与江西省地形地势较为相似,表明江西省土壤pH值分布与地形地貌因素存在较大的关联性,这一结果与图2结果相互印证。本文预测结果与周宏冀[35]采用OK模型、协同克里格和BP神经网络所获得的预测结果也基本一致。

(a) RF-A模型预测结果 (b) RF-B模型预测结果 (c) OK模型预测结果

4 讨论

4.1 模型精度分析

根据表3评价结果所示,当训练集和验证集划分比例为8∶2时,RF-A模型对于pH值的预测结果精度最高,其r为0.599,RMSE为0.392,其预测值与实测值散点图见图4。模型精度略低于齐雁冰等[36]的研究结果。主要原因可能是江西省地貌类型复杂,山地丘陵地貌区所占比重较大,山区耕地分布面积小而破碎,土壤pH值预测难度相比地势平坦、地形地貌类型单一的小区域难度更大[37]。虽然本文使用RF模型预测精度仍需进一步提高,但是预测结果基本符合预期,仍然可以为江西省耕地资源管理提供宝贵的数据基础。

图4 RF-A模型预测值与实测值散点图

4.2 pH值空间变异的影响因素分析

通过运用R Studio软件开展pH值空间变异的影响因素分析,结果显示江西省耕地土壤pH值与年均气温和年降水量关系密切(图2),这一结论与Chen等[32]、卢宏亮等[33]和王世航等[34]的研究结果基本一致。此外,沟谷深度、多分辨率谷底平坦度、高程、坡度和坡向等地形地貌因子也可以通过影响温度、光照、水分和热量的再分配以及母岩形成分布等影响土壤的发育,从而间接影响土壤pH值[38]。因此,地形因子对于预测土壤pH值也具有很高的重要性(图2)。土壤pH值与土壤有机质含量呈现极其显著负相关关系,而土壤pH值与阳离子交换量呈现极其显著的正相关关系,这一结果与戴万宏等[39]、许亚琪[40]、康婷等[41]和王晖等[42]的研究结果一致。耕层质地和土壤类型对于土壤pH值也有明显的影响,轻壤pH值明显高于其他耕层质地类型,黏土次之,不同土壤类型pH值差异也十分显著,其中草洲沙泥土、鳝泥棕红壤和棕色石灰土土壤pH值相对最高,这一结果与前人研究基本一致[32]。需要注意的是,排水条件和灌溉保证率对于耕地土壤pH值也有不可忽视的影响,本研究结果表明,耕地排水能力不满足的情况下耕地土壤pH值最高,在排水能力充分满足和基本满足的耕地土壤pH值相对较低。同时灌溉能力不满足的耕地土壤pH值也高于灌溉能力基本满足和充分满足的耕地土壤,这一结果与朱丹妹等[43]和Ding等[44]的研究结果较为相符(图5)。

(a) pH与有机质含量的相关性 (b) pH与阳离子交换量的相关性 (c) pH与土壤类型的相关性

5 结论

综合使用地形、气候、植被指数、土壤属性和耕地利用条件等自然环境协变量和人为管理因素数据,基于大量土壤实测样点,分别训练使用2种不同变量组合、不同模型参数设置和不同训练抽样比数据集合的随机森林模型,并预测江西省耕地土壤pH值,再将预测结果与普通克里格方法进行对比。得到如下结论:

1)增加土壤属性和耕地利用条件作为环境变量的RF-A模型预测耕地土壤pH值的精度优于以地形、气候、植被属性作为环境变量的RF-B模型和普通克里格的预测结果。

2)环境变量重要性分析结果显示,气候因素如年均气温和年降水量是决定土壤pH的最重要因素,地形因子、化肥施用量、排水和灌溉条件对土壤pH值变异的影响次之,分析结果与前人研究基本一致,表明相关结论较为可靠,对开展大尺度耕地土壤pH值预测具有一定的理论借鉴和实际应用价值。

但是,对于土壤pH值空间预测精度仍有待提高,今后将考虑增加环境协变量因子,提高基础数据空间分辨率等,进一步提高随机森林模型预测土壤pH值的精度。