基于随机森林模型的呼中冻土区土壤碳空间分布特征

关键词：GIS； 土壤碳储量； 冻土； RF模型； 大兴安岭

中图分类号：S159 文献标识码：A DOI：10. 7525/j. issn. 1006-8023. 2025. 01. 008

0引言

冻土区是重要的土壤碳库，由于多年冻土长期处于冷湿的环境，抑制了有机质的分解，从而累积着丰富的有机碳，仅北半球多年冻土中储存的有机碳约占全球土壤有机碳50%以上，其储量相当于大气碳库的2倍［1］。气候变化和人为活动的极易影响冻土环境，导致多年冻土区对人为活动和气候变化更为敏感［2］。大兴安岭多年冻土区是中国高纬度多年冻土的核心分布区，也是欧亚大陆高纬度多年冻土区的南缘，在全球变暖趋势下，北半球高纬度地区的增温幅度更大，使其成为最敏感的地区之一［3-4］。随着气温的升高多年冻土正逐渐消融，冻土层逐渐退化，导致封存在冻土中的碳以温室气体的形式向外缓慢释放，因此，在全球变暖的大背景下，研究土壤碳储量的空间变异及其分布特征和环境因子的关系，对于了解冻土生态系统，丰富我国冻土土壤碳储量数据库十分必要。

随着科技的发展，利用遥感数据与实地测量相结合，通过机器学习构建相关关系模型反演生态系统碳储量的方法，在生态系统碳储量的研究中得到了大量的应用。张扬［5］利用最小二乘法回归（ordinaryleast squares，OLS）、地理加权回归（geographicallyweighted regression，GWR）和地理加权克里格回归（GWR kriging，GWRK）3种回归模型分别对大兴安岭呼玛河流域不同优势树种森林地上碳储量的空间分布进行预测。魏雪梅［6］通过森林资源二类调查等数据结合遥感数据，估算了安徽省金寨县森林地上生物量，得到了较高精度的估算结果。卢宏亮等［7］利用随机森林算法（random forest，RF）模型，通过地理信息系统（GIS）和遥感影像（RS）技术，采集相关的因子及数据，对安徽省土壤有机碳进行空间分布预测。Chagas等［8］比较了RF模型和多元线性回归方法的空间预测制图的效果，认为RF模型可以避免过拟合且预测精度更高。Dharumarajan等［9］通过预测印度南部半干旱热带地区的土壤有机碳（SOC）、pH等属性，认为RF模型可以提高土壤属性空间预测的精度。RF模型合适进行土壤属性空间的预测，研究RF模型在土壤属性空间预测中的应用对数字土壤制图具有一定的意义。基于此，本研究以大兴安岭典型冻土区呼中区为研究区域，采用RF模型，基于大量土壤有机碳实测样点分别训练 2种不同变量组合数据集驱动的RF模型，进行冻土区土壤碳储量空间分布格局的预测，以期为冻土区土壤碳储量的准确估算提供理论基础。

1研究区域与方法

1. 1研究区概况

呼中区位于大兴安岭主脉伊勒呼里山的北坡，呼玛河中上游地区，地理坐标为51°14′51. 0″～52°25′28. 6″N，122°36′58. 7″～124°15′46. 9″E。地貌类型属于大兴安岭北部石质中低山山地，具有老年期的明显特征，坡度平缓，一般坡度在15°以下，局部的阳坡较陡，可达35°以上，河谷宽而平坦。海拔多在610～1 400 m，相对高差一般在200～400 m［10］。土壤以棕色针叶林土为主，分布面积最大，其次有草甸土、沼泽土和石质土等。该局气候属于寒温带大陆性季风气候。地处纬度较高，是我国最寒冷的地区之一。冬季受西伯利亚蒙西高原气团影响，严寒而干燥；夏季受海洋气团影响，温暖而多雨。年平均气温-4. 3 ℃，绝对最低气温-53. 2 ℃，被称为高寒禁区，是大兴安岭海拔最高、中国城镇气象记录最低值的地区。

1. 2样品采集和指标测定

根据研究区面积和交通可达性，对研究区选择160个样点的表层土壤为预测的目标变量，如图1所示，图1中红点为采集样点。采集区域为呼中区全域分散随机采集，包含林地、草地、湿地，以及部分城市居民区，采集表层（0～15、15～30 cm）土壤样品。将采集的土样带回室内自然风干，取适量土样研磨过筛备用。土壤总碳含量测定采用TOC分析仪测定，土壤容重采用环刀法测定。

1. 3研究方法

1. 3. 1环境变量提取

本研究基于大兴安岭呼中区Landsat8 OLI影像和ASTER GDEM高程影像得到遥感变量中的植被指数和地形因子，通过ArcGIS、ENVI空间分析模块对遥感影像数据进行指标的提取计算。共计算了4个地形因子：海拔、坡度、坡向、曲率。选取与土壤样品采集同时期的8个植被指数：地面反射率（SR）、归一化植被指数（NDVI）、土壤调节植被指数（SAVI）、优化型土壤调节植被指数（OSAVI）、改良土壤调整植被指数（MSAVI）、增强植被指数（EVI）、耐大气植被指数（ARVI）、绿色植被指数（VIGreen），来反映研究区植被覆盖情况，其中植被指数空间分辨率为250 m。将上述12个因子作为探讨影响土壤有机碳的环境因子。

1. 3.2 RF模型

随机森林（RF）的原理是通过构建多个决策树来组建随机森林［11］，将运算结果用分类回归树进行加权回归，无须筛选变量，相对于线性回归，RF模型具有很强的鲁棒性和准确性，是预测土壤属性和确定因子重要性的首选模型［12-14］。本研究通过MATLAB软件进行随机森林机器学习方法构建模拟模型。采用经典的8∶2 的分配比例将数据分为2 部分，将80%的数据用于模型训练，20%的数据用于模型精度检验。

1. 4数据处理

利用excel和spss软件对数据进行整理及相关性分析；地形因子通过ArcGIS提取、计算；RF模型在MATLAB软件中运行；运用Python3. 10. 1软件的GDAL（geospatial data abstraction library）库对冻土区土壤碳储量进行有空间分布预测。本研究中空间分布图均采用WGS_1984地理坐标系。

2数据采集

2. 1土壤属性的描述性统计特征

对呼中典型冻土区土壤总碳、容重、碳密度进行描述性统计分析，见表1。0～15 cm土层深度，土壤总碳、容重、碳密度的平均值分别为48. 56 g/kg、1. 07 g/cm³、7. 40·kg/m²，极小值分别为5. 90 g/kg、0. 64g/cm³、1.07 kg/m²，极大值为291. 60 g/kg、1. 58 g/cm³、44. 38 kg/m²。15～30 cm 土层深度，土壤总碳、容重、碳密度的平均值分别为48. 39 g/kg、1.06g/cm³、14. 80 kg/m²，极小值分别为6. 18 g/kg、0. 62 g/cm³、2. 33 kg/m²，极大值为230. 00 g/kg、1. 47 g/cm³、69. 78 kg/m²。0～15 cm 土壤总碳的均值、极大值高于15～30 cm，0～15 cm 土壤容重的均值、极大值、极小值高于15～30 cm。0～15 cm土壤碳密度的均值、极大值、极小值小于15～30cm。变异系数可以客观地反映区域内的变异程度，土壤总碳、容重、碳密度的变异系数属于变异水平中等的范畴，0～15 cm的变异系数大于15～30cm。

2. 2数字高程模型（digital elevation model，DEM）数据提取

通过呼中区DEM 数据得到地形特征，获取海拔、坡度、坡向、曲率4个地形因子，最终的高程数据如图2所示。根据提取的遥感影像数据分析其分布的地理环境，由此可以看出，呼中主要分布在海拔400～1 500 m 的范围，分布大致为由南向北逐渐减小，由东向西逐渐减小。对比海拔和坡度的图像可以看出，低海拔地区的地势较平，坡度较小。

3结果与分析

3. 1环境变量的筛选及相关性分析

利用RF 模型对遥感变量地形因子（海拔、坡度、坡向、曲率）和植被指数因子（地面反射率（SR）、归一化植被指数（NDVI）、土壤调节植被指数（SAVI）、优化型土壤调节植被指数（OSAVI）、改良土壤调整植被指数（MSAVI）、增强植被指数（EVI）、耐大气植被指数（ARVI）、绿色植被指数（VIGreen））进行变量相关性排序。由表2可知，对相关性较低的环境变量进行排除后重复建模，选取最优环境变量组合用于预测。最终确定海拔、ARVI、VIGreen、EVI、OSAVI、NDVI这6个环境变量作为自变量预测SOC含量。

3. 2随机森林模型参数设定

本研究通过逐次试验，确定训练模型中trees和leaf的最优参数。设定固定leaf值（分别设为4、5 和6），逐次调整trees 值（分别设为800、900和1000），进行3 组9 次试验。为避免过拟合问题，本研究通过比较训练集和测试集，选择两者最为接近的结果作为最优预测模型，如图3所示。结果表明，当leaf 值为5，trees值为900时，预测模型的训练集和测试集最为接近，表明此时的模型稳定性最好。

3. 3预测精度分析及碳储量空间分布

模型精度评价选用平均相对误差（MBE）、平均绝对误差（MAE）以及决定系数（R²）3个指标，其中MAE和MBE越小表明预测精度越高，训练集R²用于评价建模的拟合精度，测试集R²用于评价预测精度及模型泛化能力［15］。由表3可知，通过计算MBE、MAE、R²参数来进行RF模型的性能评估，0～15 cm土层训练集MBE、MAE、R²分别为0. 18、7. 01、0. 23，测试集MBE、MAE、R² 分别为0. 31、5. 96、0. 35；15～30 cm 土层训练集MBE、MAE、R² 分别为-0. 06、6. 84、0. 43，测试集MBE、MAE、R²分别为-1. 47、8. 93、0. 24，采用经过参数优化后的RF模型作为最终预测模型。基于ArcGIS平台，实现冻土区土壤碳储量的空间模拟。由图4可知，0～15 cm土层土壤碳储量分布大致为由南向北逐渐减小，由西向东逐渐减小。15～30cm土层土壤碳储量分布差异不明显，但南向略高于北向，西向略高于东向，土壤碳储量随土层深度增加而增加。基本上可以利用RF模型对大尺度区域上土壤属性的空间分布进行预测制图。

4讨论

1）本研究为基于Landsat8遥感数据，虽然相较于传统方法需要大量的实测数据，遥感估算方法更为便捷［16］，但由于选择数据的时间段为采集样品的同时期遥感影像，云层覆盖等天气条件会对遥感数据造成影响，该时期遥感影像少部分区域难以完全满足云量低于5%的要求，故采取了辐射定标、大气校正等方法，以消除大气等因素对地物反射的影响，但研究区遥感影像的时相差异性也在一定程度上影响研究结果［17-18］。此外本次预测只考虑了4个地形因子（坡度、坡向、海拔、曲率）、8个植被指数和其他因子（经纬度坐标）等影响，对于其他影响因子如母质、黏粒含量及其他遥感中常用的植被指数如红边叶绿素植被指数（RECL）、归一化差异红边植被指数（NDRE）等因子并没有考虑到模型中。因此实测值与估测值很难保持一致。

2）呼中整体地势呈现大致为由南向北逐渐减小、由东向西为逐渐减小的趋势，其中东南方地势最低。从空间预测分布图上也可看出，土壤碳储量的高值区也主要分布于地势较高地区，可能是因为随着海拔降低热稳定性也逐渐降低，地势较高地区要以多年冻土区为主，且海拔较高的地区冻土层较厚，多年冻土区冻土层常年处于低温冻结的情况，会使得碳的分解速率缓慢，利于碳在土壤中固存。

3）本次研究由于取样点为呼中全区域采集，采集样点分散于林地、草地、湿地及居民区，导致不同采样点区域地貌差异较大，不同冻土区植被类型分布状况不同，使其土壤持水能力、土壤入渗能力和土壤导水率存在差异，形成不同的土壤水文状况，而这些变化都可能影响植被覆盖及土壤基质分布的变化［19-20］，导致不同采样点的土壤总碳、容重之间产生较大的异质性，从而影响土壤碳储量的分布预测，故估算值和实测值会有个别数值超过了3倍标准差。

4）植被根系的分布状态会影响到土壤碳储量，植物根系分泌物可以将碳元素与保护性的矿物质相分离，导致土壤释放碳元素的速度加快。冻土层土壤环境变得相对封闭，有机质的输入量以及微生物的活动减弱，导致土壤碳元素的主要源为植物的根系，故植被根系的密度会在较大程度上影响冻土区的土壤碳储量分布。

5）在15～30cm深度的预测中分布差异较0～15cm土层的差异不明显，冻土退化过程中，含水量降低会导致碳库更容易分解释放，且受多年冻土区中强烈的冻融扰动影响［22-23］，有机物质在土壤中向下迁移，所以有机物质在多年冻土层集聚所以土壤碳含量在垂直方向上出现较复杂的波动。

5结论

本研究基于Landsat 8 OLI影像，利用随机森林（RF）模型机器学习算法建模，选择易获取的遥感变量，基于大量土壤有机碳实测样点构建的RF模型。通过12个遥感变量因子中筛选相关性最高变量，最终确定海拔、ARVI、VIGreen、EVI、OSAVI、NDVI 这6个植被指数作为自变量。结果表明，当leaf值为5，trees值为900时，预测模型的训练集和测试集的R²最为接近，表明此时的模型稳定性最好，区域尺度内能够有效地减少过拟合问题。实测冻土区0～15、15～30 cm土壤总碳的平均值分别为48. 56、48.39g/kg，极小值分别为5. 90、6. 18 g/kg，极大值为291. 60、230. 00 g/kg。0～15、15～30cm土壤碳密度的平均含量分别为7. 40、14. 80 kg/m²，极小值分别为1.07、2. 33 kg/m²，极大值为44. 38、69. 78kg/m²。以上参数预测表层（0～15、15～30cm）土壤碳储量空间分布特征，0～15土层土壤碳储量分布大致为由南向北逐渐减小，由西向东逐渐减小。15～30 cm土层土壤碳储量分布差异不明显。碳储量较高的区域多存在于植被覆盖度较高的山区，该地区森林植被群落结构稳定，人为扰动因素较小，碳储量较低的区域较集中在于居民区等人为活动较多的区域，这说明土壤碳储量极易受人为因素的影响。