基于随机森林模型的卡马兰河流域土壤有机碳储量空间分布研究

Research on the Spatial Distribution of Soil Organic Carbon Storage in the Kamalan River Basin Based on the Random Forest Model

摘要: 大兴安岭林区土壤有机碳储量对于区域碳源汇变化及气候响应研究具有重要意义。本文基于土壤剖面数据，结合多种环境变量，采用随机森林(RF)模型，估计了大兴安岭卡马兰河流域0~30 cm深度的土壤有机碳储量(SOCS)及其空间分布特征。结果表明：RF模型(R² = 0.75, RMSE = 6.81)对小尺度区域的模拟细节表现较好，研究区0~30 cm的SOCS主要分布在地势相对平坦的河流沿岸，与河流走向基本保持一致。研究结果得出了相对精确的卡马兰河流域表层土壤有机碳储量及空间分布特征，该成果能够丰富大兴安岭多年冻土区土壤有机碳储量的认识，并为生态过程相关模拟研究提供数据支持。

Abstract: Soil organic carbon stocks in the forested areas of Great Khingan are of great significance for the study of regional carbon source and sink changes and climate response. In this paper, based on soil profile data and combining multiple environmental variables, we estimated the soil organic carbon stock (SOCS) and its spatial distribution characteristics at 0~30 cm depth in the Kamalan River watershed of Great Khingan by using the random forest (RF) model. The results showed that the RF model (R² = 0.75, RMSE = 6.81) performed well in simulation details for small-scale areas, and the SOCS at 0~30 cm in the study area was mainly distributed along the river with relatively flat topography, which was basically consistent with the river course. The results of the study yielded a relatively accurate characterisation of the surface soil organic carbon stock and spatial distribution in the Kamalan River basin, which can enrich the understanding of soil organic carbon stock in the perennial permafrost region of the Great Khingan and provide data support for simulation studies related to ecological processes.

1. 引言

土壤作为陆地生态系统的重要组成部分，承载着地球上最为庞大且周转时间最长的有机碳库[1]-[3]。研究表明，全球1米深土壤内约储存了1500 Pg (1 Pg = 10¹⁵ g)的有机质形态的碳，几乎是陆地植被有机碳储量的3倍[2] [4]，表明土壤在碳储存中的核心地位。北半球陆地面积的25%被多年冻土覆盖[5]，储存了超过全球一半的土壤有机碳[6] [7]，多年冻土对维持全球碳平衡具有十分重要的作用。然而，随着全球变暖趋势的加剧，多年冻土区正经历着前所未有的退化[8]，这不仅影响了区域生态系统的稳定性，也对土壤有机碳库及其碳循环过程构成了重大挑战。IPCC预测指出，21世纪全球平均气温将持续上升，特别是在高纬度和高海拔地区，增温幅度可能更为显著[9]，这将进一步加剧多年冻土退化和土壤碳储量变化的不确定性。因此，了解高纬度多年冻土区土壤有机碳储量及空间分布特征，对调节全球碳收支、应对气候变化具有重要意义。

随着卫星遥感技术和机器学习技术的快速发展，基于土壤剖面数据和多环境因子的模型方法为土壤有机碳的估算带来了新的机遇[10]-[13]。例如，ElKamil Tola等[14]使用多元回归技术和源自Landsat影像的植被指数研究了耕作方式对沙特阿拉伯超干旱条件下农田土壤有机碳含量的长期影响；Mohammad Akbari等[15]基于遥感数据和机器学习算法推导出半干旱地区干湿月份的土壤有机碳预测模型；Kandrika Sreenivas等[16]使用RF模型，以气候、土地覆盖、岩石类型、土壤类型、多年NDVI、灌溉状况为独立输入变量，预测了印度250米空间分辨率的碳密度空间分布；袁玉琦等[17]利用RF模型模拟了福建闽东南地区的耕地表层土壤有机碳含量与环境变量之间的复杂非线性关系，评估土壤属性对SOC空间预测精度的影响，并与普通克里金(Ordinary Kriging, OK)插值法进行对比，研究发现RF模型精度明显优于OK插值法。众多研究表明，RF模型凭借其能够处理大量环境变量、捕捉复杂非线性关系以及提供高预测精度的能力，在估算土壤有机碳储量空间分布研究方面展现出了显著的优势。与传统的插值方法(如OK插值法)相比，RF算法不仅能够更准确地反映SOCS与环境变量之间的复杂关系，还能在更细的空间尺度上提供可靠的估计结果[18] [19]。

大兴安岭多年冻土区，作为我国第二大的多年冻土区，仅次于青藏高原多年冻土区[20]，这一区域不仅蕴藏着丰富的土壤有机碳储量，成为碳循环的重要“储库”，并且其多年冻土的稳定性对于维持碳储量的长期封存至关重要。针对该区域土壤有机碳储量的系统性研究仍显不足，特别是就大兴安岭卡马兰河流域这一地理区域单元而言，该地区的SOCS分布状况仍不是很清楚，亟需我们对其进行更为精确的研究。

2. 材料与方法

2.1. 研究区概况

卡马兰河流域位于大兴安岭腹地，全长约106千米，流域面积约为2207平方千米，大致介于东经122˚30′~123˚40′，北纬51˚50′~52˚30′。流域以山地地貌为主，海拔在400~1400 m之间，多年平均气温为−2℃~3.5℃，多年平均年降水量在450~550 mm之间。该区域属于寒温带针叶林区，林型以针叶林和针阔混交林为主[21]，土壤类型主要有棕色针叶林土、暗棕壤、草甸土、腐殖质沼泽土和泥炭土。

2.2. 数据获取

2.2.1. 土壤剖面数据

综合考虑研究区道路通达性和采样均匀性，使用标准环刀获取土壤剖面样品，每层取3个重复样，并同步收集坡度、海拔、地势和森林类型等环境要素数据。采集的土壤样品在野外及时烘干处理，计算每个样点不同深度土壤容重(BD)。在实验室内对采样土壤进行研磨过筛等处理，用于测量土壤有机碳含量(SOC)，并利用公式(1)计算土壤有机碳密度(SOCD)。

本研究所采用的有机碳密度SOCD_i (kg/m²)为[22]：

${\text{SOCD}}_i=\frac{\left({\text{BD}}_i\ast {\text{SOC}}_i\ast d_i\ast \left(1-V_i\right)\right)}{100}$ (1)

式中：${\text{SOC}}_i$ 代表第i层土壤有机碳含量(g/kg)，${\text{BD}}_i$ 代表第i层土壤容重(g/cm³)，$d_i$ 代表第i层土壤(cm)，$V_i$ 代表第i层土壤中大于2 mm砾石的体积百分比(%)。

若某一土壤剖面由j层采样土层组成，那么此剖面的有机碳储量SOCS (kg)为[22]：

$\begin{array}{c}\text{SOCS}={\displaystyle \sum }_{i=1}^j{S}_i\ast {\text{SOCD}}_i\end{array}$ (2)

式中；SOCS代表剖面有机碳储量(kg)，S_i代表i层土壤中与有机碳密度一一对应的实地面积，j代表该剖面的采样层数。

2.2.2. 环境数据

数字高程模型(DEM)来自地理空间数据云(https://www.gscloud.cn/home)的ASTER GDEM V3类型30 m分辨率数据集。根据DEM提取高程(Elvation)、坡度(Slope)、坡向(Aspect)、地形湿度指数(Topographic Wetness Index, TWI)、表面曲率(Surface Curvature, SC)、剖面曲率(Profile Curvature, PC)数据以及地理坐标X和Y。

遥感影像使用Landsat8 OLI/TIRS C2 L2，该数据来自于美国航空航天局(NASA)和美国地质调查局(USGS)的官方网站(https://ers.cr.usgs.gov/)。本文的波段数据选择Landsat8 OLI传感器的前7个波段(B1~B7)的地表反射率数据，并基于Landsat8遥感影像提取归一化植被指数(NDVI)和比值植被指数(RVI)。植被覆盖类型(Vegetation coverage type, VCT)从MODIS土地覆被类型产品(MCD12Q1)中提取(https://search.earthdata.nasa.gov/)，并由三个数字编码变量代表：1 = 针叶林、2 = 阔叶林和3 = 针阔混交林。纹理特征的提取是利用ENVI 5.3软件基于Landsat8 OLI遥感影像的第一主成分(PCA)进行的，提取出常用的8个纹理特征统计量(T1~T8)分别是均值(Mean)、方差(Var)、协同性(Hom)、对比度(Con)、相异性(Dis)、信息熵(Ent)、二阶性(SM)和相关性(Cor)。

年均温(MAT)和年降雨(MAP)气象数据根据彭守璋[23]的中国1 km分辨率逐月降水量和气温数据集(https://data.tpdc.ac.cn/)获得，分辨率为1 km。土壤类型(Soil Type, ST)数据根据世界土壤数据库中国1:100万土壤数据获取(http://vdb3.soil.csdb.cn/)，分辨率为1 km。

2.3. 研究方法

2.3.1. 随机森林模型

随机森林(RF)算法是一种非参数方法，它由大量训练样本生成的个体决策树组成[24]，通过将所有单个决策数的结果聚合实现预测。该方法还可以基于袋外预测的回归预测误差对输入变量的相对重要性进行排序。RF最初主要通过R语言中的随机森林包、栅格包等实现[25]，2019年ArcGIS Pro在其2.3版本中增加了这一功能模块，既可对采样点数据建模，也可直接对矢量或栅格数据建模，通过加入相应的栅格或矢量变量可直接获得预测结果的空间分布，并对多个变量的重要性进行排序，比编程更便捷高效，且对尺度区域的细节表现更好。

本研究使用ArcGIS Pro 2.3中的RF模块估计SOCS空间分布特征，参数设置为：1700棵树量，5个叶子节点，和4~18树深范围，并设置每棵树可用65%的数据量，以总数据的10%作为验证数据集。

2.3.2. 随机森林建模因子选择

本次实验得到的相关性分析结果均采用皮尔逊相关系数(Pearson correlation coefficient)。两个变量之间的协方差与标准差之商，即为两者间的皮尔逊相关系数，其公式如下：

(3)

式中：$X_i$ 和$Y_i$ 代表变量X和Y的第i个观测值，$\overline{X}$ 和$\overline{Y}$ 代表变量X和Y的均值。$R_{\left(X,Y\right)}$ 在−1到1之间取值，R = 0表示不具有相关性，值为正表示正相关，反之则为负相关。因此R的绝对值越大，表明变量之间相关性越强，反之则越弱。

方差膨胀因子(VIF)是用于评估多元线性回归模型中多重共线性程度的指标。它是存在多重共线性时解释变量方差与不存在多重共线性时解释变量方差的比值[26]。计算公式如下：

$VI{F}_i=\frac{1}{1-R_i^2}$ (4)

式中：$VI{F}_i$ 表示第i个自变量的VIF值，$R_i^2$ 表示该自变量在其他自变量线性组合中的R方。VIF值从1开始，没有上限。VIF的值为1表示给定解释变量与模型中的任何其他解释变量之间不存在相关性。1到5之间的值表示给定解释变量与模型中的其他解释变量之间存在中等相关性，但通常没有严重到需要特别注意的程度。大于5的值表示给定解释变量与模型中的其他解释变量之间可能存在严重相关性。在这种情况下，回归结果中的系数估计和p值可能不可靠[27]。

2.3.3. 精度评价指标

为了验证模型预测结果的准确性，采用决定系数(R²)、均方根误差(RMSE)来评价精度[28]，计算公式如下：

$\begin{array}{c}{\text{R}}^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(p_i-\overline{O}\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(O_i-\overline{O}\right)}^2}\end{array}$(5)

$\text{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left[P_i-O_i\right]}^2}$ (6)

其中$\overline{O}$ 为土壤有机碳密度的平均值，$p$ 为土壤有机碳密度的预测值，$\begin{array}{c}{O}_i\end{array}$ 为土壤有机碳密度的实测值，n为样点数量。

3. 结果与讨论

3.1. 土壤有机碳统计特征分析

统计结果表明土壤表层的SOC含量范围为5.81至415.69 g/kg之间，平均值为68.19 g/kg。BD是直接反映土壤松紧程度的一项物理指标，在理论及生产实践中具有多方面的实用意义，其数值大小与土壤质地、结构、松紧和有机质含量等有关。本研究区表层腐殖质含量高，土壤容重较低，土壤孔隙度大，导致表层土壤在干湿交替过程中，水分条件有利于从物理变化方面抑制土壤有机碳的分解[29]。经分析计算该流域土壤容重数值介于0.12至1.83 g/cm³之间，均值为0.72 g/cm³，接近大兴安岭林区土壤容重(0.9~1.0 g/cm³) [30]，低于柴华[31]等人关于中国土壤容重的研究结果(均值1.32 ± 0.21 g/cm³)。由公式(1)可知，本文的SOCD是根据SOC、BD数据以及剔除大于2 mm砾石计算得出，其值介于1.39至35.98 kg/m²之间，均值为9.83 kg/m²，变异系数(CV)用于描述研究变量的变异程度。相关研究认为，变异系数小于10%表示变异性低，10%~90%表示变异性中等，而变异系数大于90%表示变异性高[32]，本文SOC的变异系数为1.18，表明土壤有机碳空间分布变异性很强，BD和SOCD为中等强度变异(表1)。

Table 1. Descriptive statistical characteristics of variables

表1. 变量的描述性统计特征

变量名	均值	最小值	最大值	标准差	变异系数
SOC (g/kg)	68.19	5.81	415.69	80.27	1.18
BD (g/cm^3)	0.72	0.12	1.83	0.30	0.42
SOCD (kg/m^2)	9.83	1.39	35.98	6.91	0.70

3.2. SOCS空间分布特征

RF模型结合多源辅助变量预测表层土壤有机碳储量具有较高精度，模型拟合R²为0.75，RMSE为6.81，模型能够适用于环境复杂的大兴安岭多年冻土区。通过随机森林回归建模处理分析，得到表层土壤有机碳储量空间分布数据。

土壤有机碳储量整体呈现随坡度的增加而降低，阴坡高于阳坡的趋势。地势平坦且靠近河岸的区域土壤有机碳含量较高，这是由于河流沿岸土壤含水量较高，有助于土壤有机碳的固存[33]。河流沿岸的土壤类型有机碳含量更高的是松软潜育土(GLm)，其他区域有机碳含量较低的主要土壤类型为简育高活性淋溶土(LVh)。B. Witkowska-Walczak [34]的研究也表明，GLm在水分保持、水分传导和蒸发抑制方面表现更优于LVh，且GLm的土壤有机碳含量明显高于LVh。

4. 结论

本研究使用RF模型对卡马兰河流域多年冻土区的土壤碳储量进行预测，并分析了环境因子对土壤有机碳储量的影响。主要发现如下：

1) 通过对采样点数据汇总，得到SOC含量介于5.81至415.690 g/kg之间，平均值为68.19 g/kg，属于较高水平。研究区的SOC空间变异性很强且分布不均，这主要是土壤物理、化学和地形等环境变量的差异性导致的。

2) 影响研究区SOCS空间分布的主要环境因子是地形因子。海拔相对较高、坡度较大的区域的SOCS较低，且阴坡明显高于阳坡。

3) 预测精度方面，RF模型(R² = 0.75, RMSE = 6.81)对小尺度区域的细节表现较好，很好地展现了小尺度流域SOCS分布的精细特征。

这些发现不仅深化了对多年冻土区碳循环机制的理解，也为区域生态环境保护、资源可持续利用及应对气候变化提供了科学依据。本研究为政策制定者、环境管理者及科研人员提供了宝贵的参考信息，有助于推动区域可持续发展战略的实施，也将为卡马兰河流域多年冻土区的区域性碳循环研究提供重要的参考。