机器学习用于耕地土壤有机碳空间预测对比研究*——以亚热带复杂地貌区为例

任必武,陈瀚阅,张黎明,聂祥琴,邢世和,范协裕

(1.福建农林大学资源与环境学院/福建省土壤生态系统健康与调控重点实验室 福州 350002;2.福建农林大学公共管理学院 福州 350002)

耕地土壤有机碳(SOC)是影响土壤养分和理化性质的重要因素,也是衡量土壤质量和肥力的重要指标[1-3]。获取SOC含量及空间变化对提升土壤结构、保证粮食安全和缓解全球气候变化具有重要意义。传统的土壤制图方法因在大量采样点的采集和分析上存在困难,限制了其在大范围复杂地形地貌区的应用[4]。数字土壤制图方法,如传统的普通克里格法(Ordinary Kriging,OK)因简单、插值效果显著[5],以及具备良好的空间自相关性[6],而得到广泛应用[7],但其未能考虑土壤属性与环境因子间复杂的非线性关系[8],限制其在复杂地形地貌等SOC可能产生剧烈变化区域的应用。

随着3S技术的发展,基于机器学习算法(Machine Learning,ML)建立的辅助环境变量与土壤属性之间的预测模型已越来越多地被用于区域土壤属性空间预测。如随机森林算法(Random Forest,RF)因其处理多元非线性数据方面的优势而在SOC空间预测研究中表现较为突出[9];支持向量机算法(Support Vector Machine,SVM)在解决高维问题和非线性问题等方面表现出良好的泛化能力和预测性能[10],也逐渐被用于SOC空间预测研究。如:Wiesmeier等[11]利用RF模型成功预测半干旱草原生态系统SOC的空间分布(R2=0.76);Sreenivas等[12]基于RF模型算法利用气候、土壤、植被等辅助变量模拟并预测印度地区SOC空间分布特征(R2=0.82);Emadi等[13]利用SVM模型成功预测了伊朗地区森林环境下SOC含量(R2=0.55)。在不同方法预测SOC对比研究中,Were等[14]得出东非森林土壤环境下SVM模型表现(R2=0.64,RMSE=14.88)优于RF模型表现(R2=0.53,RMSE=17.57);而Siewert[15]则发现冻土环境下RF模型预测效果(R2=0.736,RMSE=14.13)优于SVM模型预测效果(R2=0.726,RMSE=14.9);杨煜岑等[16]对西安市SOC进行研究也发现随机森林模型表现最好(r=0.78)。

由此可见,RF和SVM等ML算法具备提高复杂地貌区SOC空间模拟精度的潜力,但不同ML算法用于区域SOC模拟的研究仍处于起步阶段,尤其在不同算法的对比研究方面仍需有所补充。不同算法可能适应不同变量组合和不同环境下的SOC预测,因此需要对不同算法的预测精度进行对比分析。此外,国内多基于单一模型对简单地貌类型的小尺度区域进行预测研究,且空间分辨率多为中低分辨率[17],不能准确反映复杂地貌区零碎耕地图斑上SOC含量,而高空间分辨率更能捕捉SOC空间分异特征。针对福建省小而破碎的耕地图斑,10 m空间分辨率更有利于捕获SOC空间分异特征,因此有必要对高分辨率影像数据SOC应用效果进行评价。

目前,全国性的土壤普查已积累了大量的土壤样点及其属性数据,可为高空间分辨率SOC预测模型对比研究提供极为便利的研究条件。基于此,本研究以典型的亚热带复杂地貌区为例,选取较易获取的地形因子、遥感植被因子、气候因子等辅助变量作为模型输入,分别通过RF模型和SVM模型预测SOC,并与普通克里格插值方法进行比较,找出适合耕地SOC预测的方法,以期为复杂地形地貌区SOC空间分布预测提供理论基础。

1 材料与方法

1.1 研究区概况

研究区位于福建省东北部(24°59′～27°4′N,118°08′～120°44′E),包括宁德、福州和莆田3个地级市(图1),南北长约880 km,东西宽约490 km,区域总面积29 829.26 km2,其中耕地面积4747.5 km2,地跨福建省最大的水系——闽江,流域面积超过60 000 km2,约占全省面积的1/2。该区位于鹫峰山脉、太姥山脉和戴云山脉之间,海拔高度400～1500 m,地形地貌复杂,以山地、丘陵、盆地和平原为主,其中山地(＞500 m)面积高达70.18%,丘陵(300～500 m)面积达20.01%。全年受到亚热带季风气候和地形影响,年平均气温为17～21 ℃,平均降雨量为1400～2000 mm,主要集中在3—8月,雨量充沛,光照充足,适合农作物生长。全区耕地以水稻土为主,占耕地总面积的79.51%。

1.2 数据来源与处理

研究区SOC数据来源于国家农业农村部2017年末测土配方调查样点数据,共计1128个(图1),每个样点均按照代表性、均匀性和适当性原则进行采样,采样深度为0～20 cm,采样的同时记录采样点GPS位置信息,最后对采集样品进行风干、筛选备用,其中土壤有机质含量采用重铬酸钾氧化-外加热方法测定。其他辅助数据包括1∶50 000福建省土地利用现状数据库和1∶50 000土壤类型数据库,分别来源于福建省国土资源厅和农业农村厅,并于ArcGIS 10.2中空间叠加后提取耕地图斑作为研究区评价底图。本研究选取与SOC空间分异密切相关的遥感植被指数、气候因子和地形因子3类辅助因子作为预测环境变量(表1)。

1.2.1 遥感植被因子

植被指数基于Sentinel-2卫星影像数据计算得到,空间分辨率为10 m,数据来源于欧洲航天局(European Space Agency,https://scihub.copernicus.eu)。影像选取植被生长旺盛的季节以有效反映有机质状况,且云量低于10%,保证数据质量。利用ArcGIS 10.2对影像数据进行大气校正、镶嵌和裁剪等预处理,统一坐标系为WGS_1984_UTM_Zone_50N,获得研究区地表反射率,最后通过波段运算(Band Math)计算比值植被指数(RVI)和归一化植被指数(NDVI)以反映植被生长状况,计算公式如下:

式中:RNIR为近红外波段,Rred为红光波段。

1.2.2 气候因子

研究区气候因子包括月最高温、月最低温和月降水量等栅格数据,空间分辨率约为4.6 km,来自世界气象数据库(WorldClim Database,http://www.worldclim.org/)。年降水量通过波段运算对月降水量进行求和计算,月最高温和月最低温分别求和,平均得到年最高温和年最低温。利用ArcGIS 10.2对3类气候因子进行镶嵌、裁剪、掩膜等预处理,并统一坐标,最后通过最邻近算法(NEAREST)重采样为10 m空间分辨率得到结果。

1.2.3 地形因子

研究区地形因子中DEM数据由福建省ASTERGDEM数据通过镶嵌、裁剪等预处理,统一坐标并通过最邻近算法将30 m空间分辨率重采样为10 m空间分辨率得到,来源于地理空间数据云(Geospatial Data Cloud,http://www.gscloud.cn/),相关衍生因子(坡度、坡向、曲率等)在ArcGIS 10.2通过Spatial Analyst模块计算得到结果。

表1 影响土壤有机碳的环境变量Table 1 Environmental variables affecting soil organic carbon

1.3 SOC空间预测模型构建

1.3.1 环境变量筛选

影响SOC含量的环境变量众多,模型训练之前需对模型变量进行筛选。机器模型因被称为“黑匣子”而不能揭示环境变量与目标变量之间的函数关系,因此将每个变量依次排除在模型之外,根据RMSE增减对变量进行筛选,其中RMSE增加则变量保留,反之剔除。

1.3.2 随机森林模型

随机森林(RF)是一种组合分类器的学习算法,在回归树模型基础上发展而来,它是利用bootstrap重抽样的方法[18]从原始样本中抽取多个样本,对每个抽取样本进行决策树建模,然后组成多棵决策树进行预测,通过投票得出最终预测结果。模型用于回归预测时,取所有回归树预测的平均值作为最后的输出结果[19]。其中n_estimators和max_depth两个参数最为重要,分别代表决策树的数量和决策树的最大深度。根据RMSE筛选后的环境变量和SOC实测值为自变量和因变量参与模型计算,通过网格搜索[14](Grid Search CV)工具设置参数选项n_estimators值(500、600、700、800、900和1000)和max_depth值(15、16、17、18、19和20)进行逐次参数组合计算,根据训练集和验证集R2最为接近为原则,确定决策树数量(n_estimators=600)和决策树的最大深度(max_depth=16)为最优参数对SOC进行预测。RF模型的构建和预测通过python中Random Forest Regressor模块实现。

1.3.3 支持向量机

支持向量机(SVM)是在统计学理论的基础上,以结构风险最小化为原则建立起来的机器学习模型,通过生成最优分离平面,控制参数、调节模型结构,实现经验风险和结构风险的最小化[20]。SVM回归预测基于不敏感函数及核函数算法进行计算,针对非线性回归,常通过非线性映射核函数(Φ)把数据映射到高维空间进行线性回归处理[21],其中惩罚系数(C)和不敏感损失函数(ε)两个参数最为重要,用于平衡误差和调整模型复杂程度[22]。根据RMSE筛选后的环境变量和SOC实测值作为模型输入,通过网格搜索工具(Grid Search CV)设置参数选项C值(5、8、10、12、15、18和20)和ε值(0.01、0.005、0.001、0.0005和0.0001)进行逐步参数组合计算,根据训练集和验证集R2最为接近为原则,确定惩罚系数(C=10)和损失函数(ε=0.001)为模型外推的最优参数进行SOC的空间预测。SVM模型的构建和预测通过python中SVR模块实现。

1.3.4 普通克里格

普通克里格(OK)基于SOC的空间自相关性,进行内插或外推来达到预测SOC的目的,其实质就是对SOC实测数据进行线性无偏最优估计[23]。OK模型基于SOC实测样点于ArcGIS 10.2中通过3D Analysis模块进行空间插值得到预测结果。

1.4 模型对比

1.4.1 模型精度对比

针对3种模型(RF、SVM、和OK模型)整体精度,采用80%的训练样本和20%的验证样本进行SOC的模拟与验证,通过均方根误差(RMSE)、平均绝对误差(MAE)、决定系数(R2)和相关系数(r)定量化模型精度并进行比较,其中RMSE、MAE越小,R2和r值越接近1,表明模型预测精度越高;不同范围SOC模型精度通过不同比例5次随机抽样计算RMSE平均值得到(抽样比例分别为30%、40%、50%、60%、70%、80%、90%),避免样点在不同范围数量较少带来误差的不准确。

1.4.2 模型变量重要性对比

去除某一环境变量计算RMSE与全部环境变量RMSE的差值代表该变量重要性大小,差值越大表明该因子相对重要性程度越高,以此比较同一机器模型不同变量的重要性大小以及不同模型反演SOC环境变量组合间的差异。

1.4.3 模型预测SOC空间分布对比

比较RF、SVM模型和OK模型预测SOC空间分布图,对差异性较大的区域放大对比,用于评价SOC预测模型空间异质性的优劣。

1.5 区域SOC反演

基于研究区耕地图斑对预处理好的辅助变量进行裁剪得到耕地对应的遥感植被因子、气候因子、地形因子等环境变量,并作为机器模型的输入对象,最后得到不同机器模型预测SOC空间分布格局;OK模型预测结果基于SOC实测样点于ArcGIS 10.2进行插值并重采样为10 m得到。ArcGIS 10.2中完成SOC空间分布图的制作。

2 结果与讨论

2.1 不同模型预测SOC精度对比

基于SOC实测样点对不同模型性能进行统计,结果如图2。根据RMSE、MAE、R2和r4个模型指标比较:RF模型得到最低的预测误差和最高的r值(RMSE=2.004,R2=0.880,r=0.897),模型表现最优,能解释SOC 88%的空间变异性;相较于RF模型,SVM模型和OK模型中RMSE分别增加159.0%和128.1%,误差较大。在SOC所有含量范围的误差中RF模型均低于其他模型(图2d),表现最好;在SOC为10～15 g·kg−1和15～20 g·kg−1范围内SVM模型误差低于OK模型误差,而在其余SOC含量范围SVM模型误差高于OK模型。

如表2,不同方法预测SOC中值和平均值都与实测结果接近,约为15 g·kg−1,但变异系数和标准偏差均变小,预测结果范围被明显压缩,这与Siewert[15]预测亚北极SOC结论一致,主要因为算法回归趋于平均所导致;相较于RF和OK模型,SVM模型这一现象更为明显,压缩范围高达56.03%;从标准偏差(SD)和变异系数(CV)发现,3种算法中,RF模型预测结果的SD和CV与实测值最为接近,更能表征SOC的动态变化范围。总之,3种预测模型中,RF表现最佳,其次OK模型,SVM相对最弱。

表2 不同模型土壤有机碳(SOC)预测值与实测值对比Table 2 Comparison of predicted by different models and measured soil organic carbon(SOC)values

2.2 RF和SVM模型变量重要性对比

图3所示,RF模型中重要性占主导地位的变量主要包括最低温度(Mint,87)、高程(DEM,75)和降水(Rainfall,48);SVM模型中重要性占主导地位的变量主要包括降水(Rainfall,93)、高程(DEM,28)和地形起伏度(Rel,22),其中Rainfall和DEM在山区SOC的贡献率表现突出[24]。Mza等[25]认为地形和降水在一定程度上改变了土壤性质促进SOC的积累,杨煜岑等[16]和任丽等[17]在SOC研究中也发现DEM是影响SOC最重要的变量之一;且张厚喜等[26]和钟兆全[27]分别运用不同模型预测福建省SOC,发现高程是影响SOC含量的重要因子,且SOC含量随海拔的升高而增加,与本研究观点一致。

两种模型中遥感植被指数NDVI和RVI对SOC重要性都低于高程和降水,但都作为贡献因子参与模型预测,对SOC空间预测不可或缺。Guo等[28]预测海南橡胶园SOC,齐雁冰等[29]针对陕西省SOC研究中也发现遥感植被因子重要性均低于高程,且卢宏亮等[30]和马冉[23]在SOC空间预测研究中均得到相同结论。SVM模型中坡向(Aspect)因子被排除不参与SOC建模;RF模型中地形起伏度(Rel)、地形湿度指数(TWI)和平面曲率(Plan)3个因子被排除不参与模型计算,其余因子[剖面曲率(Profile)、坡度(Slope)、年最高气温(Maxt)]对模型的贡献率较低。年最低气温(Mint)在RF模型中贡献率较高而在SVM模型中贡献率较低。

2.3 土壤有机碳空间分布

图4展示了不同模型预测的SOC空间分布状况。总体上看3种模型预测的SOC空间分布相似,呈现出北部、中部高于南部地区,西部高于东部沿海地区的分布态势。统计分析发现(图5),SOC空间分布特征与3类环境变量代表性因子(DEM,Rainfall,NDVI)呈现明显的相关性。针对亚热带复杂的高山地貌环境,一般海拔越高,降水量越高,越有利于SOC的积累,这与本研究中SOC与高程和降水呈正相关的结论一致(图5a,b),且这一结论被多位学者所证明[31-33]。植被是SOC含量的重要来源,控制着土壤有机质含量的输入,故遥感植被指数NDVI与SOC含量呈现明显的正相关关系(图5c),这与Shi等[34]观点相一致。研究区北部、中部和西部地区(SOC高值区)高海拔高降水的土壤环境有利于SOC的积累和保持,形成较高的SOC分布格局;东部沿海与南部低海拔地区SOC含量低于北部、中部和西部高海拔地区,可能是SOC积累速率较低导致。

选取不同模型预测SOC分布差异较大的子区域(图4Ⅰ-Ⅲ,a,b框)对SOC空间异质性进行比较发现:机器模型(RF模型和SVM模型)所表达空间异质性更为精细,而OK模型空间表达相对粗糙,主要归因于OK模型仅考虑SOC空间自相关因素,而缺乏对影响SOC多种环境因子的考量[35]。图2d不同范围SOC误差对比得出:RF模型在SOC所有范围误差最小,故RF模型预测SOC空间分布更为平滑,预测结果更为精细(图4Ⅰ);SVM模型在高值区(＞25 g·kg−1)和低值区(＜5 g·kg−1)存在较大的预测误差(RMSE＞9.2,图2d),造成SOC预测值过于趋向平均,其SOC空间分布上表现为高值低值区较少(图4Ⅱ),不能完全反映SOC的动态变化范围。最后,无论从模型精度还是空间异质性表达,RF模型表现最优,可用于亚热带复杂地貌环境的SOC预测。以下仅对最优模型预测的SOC空间分布状况进行描述。

RF模型预测SOC含量均值为15.33±4.07 g·kg−1,范围是3.57～29.91 g·kg−1,其中SOC中高值区(＞20 g·kg−1)主要分布在鹫峰山脉和戴云山脉的高海拔地区,仅占研究区总面积的4.13%;SOC中低值区(＜10 g·kg−1)主要分布于福州平原和兴化平原等低海拔的沿海地区,占研究区总面积的5.15%;10～20 g·kg−1区间所占面积超过90%,位于海拔高低过渡带之间。

2.4 局限性分析

福建省地形地貌环境复杂,本研究选取了容易获取的地形因子、气候因子以及遥感植被因子参与模型构建并预测SOC,但缺乏人类活动对SOC影响的研究。实践证明农业活动[36](如:轮作、灌溉、施肥)等人为要素对SOC尤其表层SOC含量产生重要影响,如:尹萍[37]认为农业活动会改变表层SOC含量,影响气候等自然环境变量与SOC的关系;田慎重等[38]研究发现保护性轮作措施能有效提高表层SOC含量。因此,寻找更多SOC相关性强的辅助变量以及能代表人类活动的替代因子作为模型输入将是提高SOC预测准确度的重要途径之一。

本研究仅针对两种机器模型(RF模型和SVM模型)预测SOC,并与OK进行比较,机器模型应用较少,而未来SOC研究需要更多的机器学习方法进行对比研究,如:分类树[39](Classification Tree models)、朴素贝叶斯[40](Naive Bayesian Classifier)、人工神经网络[41](Artificial Neural Network)等,随着机器模型的发展,深度神经网络(Deep Neural Networks,DNN)也被应用于数字土壤制图,如:Emadi等[13]研究伊朗北部地区SOC发现DNN模型优于其他机器模型,表现最好。因此,多种机器模型的评估和比较更有利于对研究区SOC空间分布状况进行整体而全面的把握。

3 结论

本研究基于福建省复杂地貌区大量实测样点和10 m空间分辨率Sentinel-2卫星影像数据,选择植被遥感变量、地形变量以及气候变量作为模型输入,重点比较RF模型和SVM模型在模型精度、变量重要性和空间分布的差异,并与OK模型进行对比。结果显示:1)从模型误差和相关系数比较,RF模型表现最好(RMSE=2.004,R2=0.880,r=0.897),OK模型次之,SVM模型表现相对最差(RMSE=5.190,R2=0.193,r=0.431);2)RF模型和SVM模型重要性变量选择上,高程和降水最为重要,与SOC呈正相关关系,而遥感植被因子重要性低于高程,也与SOC呈正相关关系;3)3种模型预测SOC空间分布趋势总体一致,表现为:北、中、西部高海拔地区SOC含量高于南部和东部低海拔地区,相较于SVM和OK模型,RF模型能体现更多的空间变异性信息,空间异质性表达更为精确,可作为复杂地貌区SOC含量预测的高效方法。