基于多源地理数据和随机森林模型的土壤类型模拟预测研究——以宁洱县为例

卢加华

(云南省地矿测绘院有限公司,云南 昆明 650218)

0 引言

土壤普查是查明土壤类型、理化性质及其空间分布规律,查清土壤资源数量和质量等的重要手段,普查成果可以为土壤的科学分类、改良利用、保护管理等提供科学依据,也可以为国民经济建设和地方政府重大政策的制定提供决策依据。我国先后于1958年和1979年开展了两次全国土壤普查。现如今距第二次全国土壤普查(简称二普)已经有40余年,二普的普查成果数据的现势性已经不能全面如实地反映当前农用地的土壤质量情况。2022年2月16日,国务院印发《关于开展第三次全国土壤普查的通知》,决定自2022年起开展第三次全国土壤普查。当前,全国上下正在紧锣密鼓地开展第三次全国土壤普查(简称三普)工作。通过普查,可以对耕地、园地等农用地和未利用地的土壤类型、性状、立地条件等进行“全面体检”,摸清土壤质量家底。

土壤类型模拟预测,是第三次土壤普查成果形成的核心工作。土壤属性图、土壤农业适宜类型评价、土壤志等成果形成,都需要基于土壤型图来开展,土壤类型模拟预测的重要性不言而喻。本研究在宁洱县土壤三普试点工作的基础上开展,旨在研究土壤类型模拟预测实用方法和存在的问题,以指导后续土壤三普工作按时、按质、按量完成。

1 数据来源及研究方法

1.1 区域概况

宁洱县位于云南省南部,在普洱市中部,距省会昆明市370 km,南距市政府驻地思茅区33 km。地处东经100°43′至101°37′,北纬22°41′至23°36′之间。北部与镇沅县相接,东北部和东部沿把边江与墨江县分界,东南与江城县山水相连,西南与思茅区接壤,西北沿小黑江与景谷县为邻。全县总面积3 670 km2,其中,县城建成区面积7.6 km2,山区面积占96.8%。全县辖6镇3乡、89个村(居)民委员会,总人口16.27万人,城镇化率45.3%,少数民族人口占55.5%。全境山区较多,地质构造属喀斯特地貌,地貌纵横交错复杂,气候属于南亚热带山地季风气候,并兼有南温带、热带、中亚热带等气候类型。

1.2 数据来源

数据类型涵盖了地质、地形地貌、遥感、土地利用类型和土壤三普外业采样等成果。其中,1∶25万地质图来源于云南省地质矿产勘查开发局,行政区划图、数字高程模型(DEM)、土地利用类型等数据来源于宁洱县自然资源局,1∶5万土壤名称校准后土壤图和样点制备数据来源于云南省第三次土壤普查办公室下发给宁洱县第三次土壤普查办公室的数据,高分多光谱遥感影像从地理空间数据云(https：∥geocloud.cgs.gov.cn/)网上获取。

1.3 随机森林模型

经典的随机森林(Random Forest,简称RF),是由Breiman提出的一种基于分类回归树的非参数机器学习模型,能从数据中学习复杂的非线性关系和交互特征[1]。它通过建立多个决策树来完成分类或回归任务,并通过集成这些决策树的结果来作出最终的预测,具有较高的预测准确性和稳定性。

随机森林模型中存在大量需要调优的超参数,其中一部分来自于基决策树,一部分来自于最终的随机森林模型自身[2]。主要参数包括基决策树棵树n(ntree)和特征数量m(mtry)。基决策树棵树n(ntree)是指随机森林所包含的决策树数量,默认为500。若ntree太小容易欠拟合,太大容易过拟合,可根据误差大小调整决策树的数量,默认情况下取数据集变量个数的1/3。通常情况下,需要进行反复训练来确定最佳值。

2 模拟预测结果

2.1 成土环境因素确定和预处理

环境变量的选取原则是基于土壤发生学理论,综合考虑研究区域的土壤景观特点和成土环境条件,选取与土壤类型形成与演变相关或协同的环境因素变量,结合宁洱县实际情况及已有资料完整性来分析。本研究选取母岩母质、植被覆盖指数、土地利用现状、高程、坡度、坡向、坡位、平面曲率、剖面曲率、地形湿度指数共10种成土环境因素。确定成土环境因素后,利用GIS软件对10个成土环境因素进行处理,转换为格式统一栅格数据,使之能被R语言识别,并进行回归分析,预测土壤类型。

1)母质(MZ),采用校正后的二普土壤类型图与样点数据进行空间连接并赋值,得到母质分布图。为便于R语言识别,将母质类型按数字连续编号处理,研究区域的范围为1～11。

2)归一化植被指数(Normalized Difference Vegetation Index,NDVI),是一种常用的遥感指数,用于评估地表植被状况。它通过计算红外波段和可见光波段的差异来反映植被的活力和繁茂程度。归一化是指将NDVI的值范围进行统一化,使其在0到1之间,使得不同数据集之间的植被指数可以进行比较和分析。研究区域的归一化植被指数范围为-0.21～0.69。

3)土地利用类型(TDLY),土地利用方式也是影响土壤养分分布的重要因素,但土地利用方式为类别变量,不能直接用于回归分析。本研究采用哑变量方法处理土地利用类型,哑变量处理后,对数值做连续对应关系处理,处理后研究区域的土地利用类型取值为2～14。

4)地形湿度指数(Topographic Wetness Index,TWI)是一个地表湿度评估指标,用于描述地表的排水性和潮湿程度。地形湿度指数根据地形因素计算得出,包括地形斜率和流域贡献面积。地形斜率越大,流域贡献面积越大,则地形湿度指数越高,表示该区域可能有较高的湿度和排水性较差。研究区域的地形湿度指数取值范围为2.56～23.64。

5)高程数据(GC)、坡度(PD)、坡向(PX)、坡位(PW),均由数字高程模型(DEM)生成。对坡度重分类为5级,坡位重分类为上、中、下3级,对应数值为1、2、3。

6)平面曲率(PMQL)及剖面曲率(POMQL)代表地形表面的凹凸程度,其间接影响土壤发育。平面曲率值为正值时,表明地形表面为向上凸,负值表明地形向下凹,值为0时表面为水平。剖面曲率与平面曲率相反,正值表明地面向下凹,负值表明向上凸[2]。

主要成土环境因素的空间分布特征见图1。

图1 主要成土环境因素的空间分布特征Fig.1 Spatial distribution characteristics of main soil forming environmental factors

2.2 随机森林模型训练

2.2.1 模型训练及参数调优

首先,基于土壤三普剖面样点、表层样点制作样点数据集,并以此来训练随机森林模型。土种为因变量,10种成土环境因子为自变量。预测结果表明,仅选择三普采样点作为样本训练模型,可解释度(Mean of squared residuals：81.78)和整体解释率(% Var explained：9.54)均不高,且预测出的土种类型缺失,共缺失10个土种。

其次,提取一定数量的典型虚点来参与随机森林模型的构建。从二普土壤图上拾取土壤类型典型点(虚点,非实际调查观测点)作为补充性样本点。通过人机交互的方式,筛选典型虚点数量和空间分布的合理性。对土种样点数量少于10个、图斑面积较少或空间分布较零碎的土种进行空间插值处理。全县共提取1 190个典型虚点,典型虚点主要分布于样点少于10个、土种面积较小和空间分布零碎的水稻土、灰泡土等土种。

融合典型虚点和三普外业采集的样点数据形成样点数据集,来训练随机森林模型。经多次训练结果对比显示,可解释度(Mean of squared residuals：98.25)和整体解释率(% Var explained：23.32)分别提高了16.47和13.78个百分点,且各成土环境因子重要性得分均有提高,其中母质重要性提高了约15个百分点。增加典型虚点前后的重要性对比如图2所示。

图2 原始样点与增加典型虚点重要性对比图Fig.2 Comparison of importance between original samples and added typical virtual points

确定样点数据集后,对随机森林模型中的超参数经过不断的调试挑选,最终确定n(ntree)最佳取值为1 000,m(mtry)最佳取值为4。

2.2.2 环境变量重要性分析

对利用土壤类型与环境变量关系,进行土壤类型与环境变量之间的相关性分析,保证两者之间存在显著相关性,以判断哪些环境变量可以保留在模型中,并去除环境变量之间的共线性。采用成土环境因子重要性得分(%IncMSE)来评价相关性,是通过对每一个成土环境因子进行随机赋值,如果该成土环境因子更为重要,那么它被随机替换后模拟预测的误差会增大。因此,该值越大表示该变量的重要性越大。从表1可以看出,成土环境因子母质的重要性得分最高,归一化植被指数次之,坡位的影响最小。

2.3 土壤类型空间推测

将最终筛选样点数据集、训练好的随机森林模型、成土环境因子栅格数据集成到R语言中。依次读入预处理好的成土因子栅格变量数据,将栅格转为矩阵和向量,再将向量合并为数据框,然后,利用Predict函数预测土种类型的栅格分布。研究区域共43个土种,用连续数字分别进行编码,预测出38个土种,预测完整率为88.37%。其中,有5种土种因图斑面积太小未被预测出来,分别为黑香面土、红砂土、棕末香土、黄胶泥田和灰泡土,通过外业校核来完善土壤类型图。 R语言模拟预测土种图见图3。经GIS处理校准后的土种图见图4。

图3 R语言模拟预测土种图Fig.3 Soil type map of simulated prediction using R language

图4 经GIS处理校准后的土种图Fig.4 Soil type map after processing and calibration by GIS

2.4 土壤图校核与更新

2.4.1 土壤类型图校核

土壤图校核工作包括内业和野外校核。内业校核结合第三次全国国土调查数据和成土环境因子进行边界调整,调整依据为地形地貌、母质、植被、土地利用等在景观上的明显变异点。

利用GIS软件提取土壤类型名称或边界疑似发生改变区域和模拟预测不出来的土种图斑,经人工筛选优化处理后,形成的校核图斑。全县共提取54个图斑,分三条路线开展野外校核工作。野外校核工作组由土壤调查分类、土壤制图专家和熟悉当地土壤类型的专家组成。通过打钻和专家经验现场判别土种类型,利用GNSS记录校核点的经纬度坐标、景观部位和土壤利用情况等信息,验证模拟预测的土种名称和土壤类型图斑边界正确性。外业结束后,内业根据外业校核反馈的结果来修正模拟预测结果。

通过野外校核,发现土壤名称发生改变的图斑有3个,分析其原因主要有两方面：一是二普时候为水田,后长期水改旱,导致耕作层已经不具备水稻土特征;二是二普有两个名称,对应到三普只有一个名称,干扰机器学习结果。绝大部分图斑边界均需要调整,调整原因多数为土地利用类型发生改变,导致土壤类型边界改变。

2.4.2 土壤类型图更新

通过野外校核和内业边界调整,获取土壤类型改变区代表性图斑的土壤类型变化情况,经过归纳整理,形成县域内土地利用变更等原因导致土壤类型变化的知识规则,根据这些知识规则对土壤类型改变区进行土壤类型和边界更新。将土壤类型改变区更新图斑与土壤类型未改变区更新图斑在ArcGIS软件中进行合并和融合,生成新的土壤三普土壤类型图。

经过分析发现,引起宁洱县土壤类型发生改变的知识规则主要有以下几种情况：一是水田改为旱地、园地、林地、草地,导致耕作层已不具备水稻土特征;二是旱地、林地、草地等改为水田,耕作层具备水稻土特征;三是通过占补平衡措施,采取覆土、填埋等方式建成的新增耕地;四是潜育化土壤因水分条件变化脱潜;五是表土层因土壤侵蚀导致表土层变薄或表土层消失;六是水土流失、酸雨等其他原因。

3 结论

宁洱县第三次土壤普查土壤类型,经过GIS软件处理多源地理数据,选取母岩母质、植被覆盖指数、土地利用现状、高程、坡度、坡向、坡位、平面曲率、剖面曲率、地形湿度指数10个成土环境因素,在R语言中模拟预测土壤类型生成土壤类型栅格图。经过内业和野外校核,得出以下结论：

1)基于多源地理数据和随机森林模型的土壤类型预测方法具有较高的准确性和可靠性,预测的土壤类型边界与实际吻合度较高。

2)提高外业采样点分布的合理性,可以提高预测成果的准确性。

3)提取合理的典型虚点作为样本数据集,可以提高预测成果的准确性。

4)预测土种有一定的准确率,但需要结合外业采样点和第二次土壤普查成果的土壤类型图来校准土种名称。

此外,通过对宁洱县土壤类型模拟预测研究发现,宁洱县第三次全国土壤普查试点工作还存在以下问题：

一是表层样采样样点布局不合理。表层样样点过分集中于某几种土种,少部分土种样点数据过多,过度拟合;大部分土种样点数量过少,欠拟合。从而导致样点少的土种不能被预测出来,过度拟合的土种边界与实际情况不吻合,加大了外业校核工作量。

二是土壤类型图边界校核工作前后倒置。受各种因素的影响,三普土壤类型图斑边界校核工作滞后。外业取样工作完成后才开展土壤类型图边界校核,致使部分区域外业工作重复。

三是利用低影响高精度底图来校核高影响低精度地图。如：重要性最高的母质图比例尺为1∶25万,而土地利用类型数据精度为1∶1万,重要性排名为第四位。用低影响高精度去校核高影响低精度比例尺地图,需要大量的野外校核工作来核实确认边界的正确性。

因此,结合本研究成果和试点成果形成工作,对土壤类型图制作提出几点建议：

一是加大外业采样样点布局分布合理性检查。在样点校核布局的时候,各级三普办应加大各土种表层样取样的数量和空间分布的合理性检查,对样点数量少于10个的土种进行加密,建议每个土种图斑不少于3个样点。

二是将土壤类型图边界校核工作前置。建议将土壤类型图边界校核工作和表层样采样工作同步开展,减少外业工作量。通过外业判别对土壤类型边界发生改变区域进行标注,大幅提高土壤类型图边界校核的准确性和真实性。

三是采取全图斑校核方法,提高土壤类型图精度。建议依托地方政府和基层专家的力量,采取沿图斑边界发生改变的区域开展全图斑野外校核工作,以提高土壤类型图边界的正确性和准确性。