2025, Vol. 57
2. 安徽理工大学空间信息与测绘工程学院, 安徽淮南 232001;
3. 中国科学院大学, 北京 100049
黑土是世界上最肥沃的土壤，具有高可耕性，土壤肥力极高[1]。东北黑土区是我国最重要的商品粮供应地之一[2]。但由于我国黑土区耕地长期高强度、超负荷使用，引起土壤质量下降、功能衰退，导致产生“变薄、变瘦、变硬”等退化问题，严重影响了农业和环境可持续发展[3]。土壤pH、土壤有机质(SOM)和土壤全氮(TN)作为关键评价土壤质量和土壤肥力的属性特征，在土壤养分和植物生长发育过程中具有重要地位，对土壤质量和关键功能产生重要影响[4]。因此，为探究黑土地土壤质量，最大程度利用土地资源，提高农业生产效益，刻画黑土区土壤属性的空间分布特征，深入了解黑土地土壤状况，显得尤为重要。
近年来，数字土壤制图技术广泛应用于土壤学领域，已成为一种新兴的、高效表达土壤空间分布的方法[5]。随着3S等信息技术的快速发展，国内外许多学者提出了多种数字土壤制图的方法，用以描述土壤属性的空间分布规律。目前，常见的方法包括地统计学、回归分析法、机器学习和混合模型等。其中基于地统计学方法的普通克里格法(Ordinary kriging，OK)最具代表性。徐尚平等[6]和Duan等[7]采用OK法分别对内蒙古土壤pH与SOM和湖南省土壤有机碳(SOC)的空间变异规律进行了分析。结果表明，在不同尺度预测不同土壤属性空间变异规律中，OK法均能取得较好的预测结果。然而，OK法只对单一指标进行插值，并未考虑环境变量(例如气候因子、地形因子和遥感指数等)对土壤属性的影响[7]。以结合辅助环境变量的回归克里格(Regression kriging，RK)相比OK法，能更好地揭示土壤属性的空间变异规律。张素梅等[8]和杨顺华等[9]分别采用RK法对吉林省农安县SOM、TN和宜都市SOM进行了空间分布预测。结果表明，结合辅助环境变量的RK法能够有效提高土壤属性的预测精度。此外，由于土壤发育过程中影响土壤属性的环境变量具有空间非平稳性，传统的地统计学方法不能很好地体现土壤属性的局部特征[10]。作为一种有效处理回归分析中空间非平稳性和空间依赖性的局部模型，地理加权回归模型(Geographically weighted regression，GWR)近年来在土壤属性制图中得到广泛应用[11-12]。例如，罗梅等[11]和Costa等[12]分别采用GWR等多种模型对中国和巴西东南部山区SOC空间分布进行预测，并对不同模型的预测性能进行分析与评价。
上述克里格法和回归分析法均是对变量进行线性估计的方法，然而实际中土壤属性与环境变量间存在着复杂的非线性关系[13-15]。因此，许多学者将机器学习算法如决策树[13]、随机森林(Random forest，RF)[14-15]、支持向量机[16]和回归树[17]等应用于土壤属性制图。其中，RF不仅能够有效克服过度拟合和不稳定的缺陷，而且还能很好地处理变量间复杂的非线性关系，近年来在土壤属性预测中应用广泛[15，18]。例如Liu等[19]应用RF生成了我国第一版90 m空间分辨率国家土壤信息网格数据集(包含pH、SOC和TN等)；Song等[20]采用RF对东北松嫩平原SOC和TN含量进行了空间分辨率为90 m的空间预测；庞龙辉等[21]运用RF对青海省pH、SOC和TN等属性进行了空间分辨率为1 km的预测，决定系数(R2)均高于0.47。此外，一些学者将RF模型的预测残差与预测结果相结合，形成混合模型——随机森林-普通克里格法(Random forest-ordinary kriging，RF-OK)，该方法也能够表现出良好的预测性能[22-24]。Guo等[22]采用RF和RF-OK对海南岛橡胶林的SOM进行空间分辨率为30 m的预测制图，研究表明，RF-OK具有较低的均方根误差(RMSE)和较高的R2，在预测SOM空间分布方面优于RF。Mousavi等[23]利用RK和RF-OK对伊朗马什哈德平原的SOC进行了空间分辨率为30 m的预测，其模型R2分别为0.33和0.86，RMSE分别为0.48和0.24。Zhang等[24]运用RF和RF-OK对我国西南部复杂地形的SOC进行了空间分辨率为30 m的预测，其RMSE分别为28.09% 和27.61%。综上所述，当前没有一种数字土壤制图方法能够适用于所有不同种类土壤属性的预测，每种方法都有其适用的环境，而且不同土壤属性对应的最佳预测方法也不同。因此，选择最佳的数字土壤制图模型对于提高区域土壤属性的制图精度非常重要。
友谊农场地处东北典型黑土区，位于黑龙江省友谊县境内，是我国首个大型机械化国营农场，有“中国第一农场”称号。然而，当前黑土区农场尺度土壤属性数据比较粗略，其空间分辨率仅为90 m，难以满足现代精准农业的迫切需求。因此，本研究以友谊农场核心示范区为研究区，选用OK、GWR、RF和RF-OK 4种代表性数字土壤制图方法，对示范区土壤pH、SOM和TN含量进行空间分布预测及不确定分析。主要目标包括：①采用空间分辨率为12.5 m的DEM和10 m的哨兵影像数据作为主要环境变量，以提高土壤属性制图精度；②对比分析4种模型对示范区土壤pH、SOM和TN的预测性能，选择最优模型阐明示范区土壤属性的空间分布格局。研究结果将为友谊农场核心示范区的高效、智能和精准农业管理提供数据支持与科学参考。
1 材料与方法 1.1 研究区概况研究区友谊农场位于黑龙江、松花江和乌苏里江三江平原腹地，是我国重要的商品粮基地之一和现代化农业示范窗口。友谊农场核心示范区地处友谊县、富锦市和宝清县三地，地理位置为131°27′E ~ 132°17′E，46°35′N ~ 46°55′N(图 1)，东西长约64 km，南北宽约17.5 km，耕地面积约为974.5 km2，占总区域面积的87% 左右。示范区西南为丘陵，东北为低洼地，由西南向东北倾斜，海拔为71 ~ 413 m。年均温约为3.2 ℃，年均降水量506.03 mm，气候适宜，为农业生产提供了优越的条件。土地利用类型主要以耕地为主，作物主要有水稻、大豆和玉米。土壤类型以草甸土、沼泽土、白浆土为主。友谊农场作为我国现代化农业的典型代表，对于我国的粮食安全和农业现代化建设具有重要的意义。
|
图 1 研究区概况与样点分布 Fig. 1 Overview of the study area and distribution of sample sites |
通过综合考虑研究区的土壤类型、土地利用和道路可达性等因素，在遵循代表性和空间分布均匀性原则的基础上，采用近似网格方法共布设土壤样点152个(图 1)。于2022年10月中旬进行野外采集，使用手持GPS仪确定样点地理位置，并采用梅花形五点采样法获取5份表层(0 ~ 20 cm)土壤样品后进行充分混合。随后使用四分法采集约1.5 kg土壤样品，经风干和研磨，并通过20目和100目筛后测定样点土壤pH、SOM和TN含量。其中，土壤pH使用1∶5土水质量比浸提电位法测定，SOM含量采用重铬酸钾氧化-硫酸消化法测定，TN含量采用半微量凯氏定氮法测定[25]。
1.3 环境变量获取及预处理基于SCORPAN函数模型[26]，遵循相关性和可获得性原则，选取土壤因子、地形因子和遥感生物指数作为预测示范区土壤属性的环境变量。其中土壤因子来源于全国尺度空间分辨率为90 m的土壤质地栅格数据[27]。地形因子直接影响着地表物质的能量循环以及土壤的发生与演变，是土壤制图中常用的环境变量[5，28]。采用ALOS-12.5 m DEM地形数据，来源于ALOS的PALSAR传感器(https://search.asf.alaska.edu)，在DEM数据基础上，通过SAGA-GIS软件提取出地形衍生因子。生物因子间接由遥感影像的特征波段及不同组合有效反映土壤属性形成的地表状况和植被景观特征[5，28]。采用由欧洲航天局(ESA，https://scihub.copernicus.eu)公开数据平台下载的空间分辨率为10 m的Sentinel-2A卫星数据，获取影像时间为2022年8月，云量较少。随后对获取到的影像数据进行了格式转换、拼接、投影变换和重采样等预处理。所选环境变量见表 1。
|
表 1 环境变量基本信息 Table 1 Basic information on environmental variables |
为便于后续建模，使用ArcGIS 10.8软件将上述所选环境变量空间分辨率重采样至10 m，并统一其空间范围和坐标系。在后续建模中，将上述29个辅助因子均作为RF和RF-OK模型的输入变量进行属性预测。而考虑到GWR模型作为一种局部线性回归方法，为避免变量共线性对模型预测精度的影响，在建模之前，调用R语言car包中的vif函数计算了所有环境变量的方差膨胀因子(VIF)，并剔除VIF > 10的变量[10]，仅保留剔除后的部分变量进行建模。图 2展示了研究区部分环境变量。
|
图 2 研究区部分环境变量 Fig. 2 Selected environmental variables in the study area |
普通克里格法(OK)：又称空间局部估计或空间局部插值法，是利用输入数据和半方差函数的结构，在有限区域内对区域化变量进行无偏最优估计的一种方法[29]。其基于土壤属性的空间连续性，通过邻近的相关观测点权重均值来预测采样点位置的土壤属性值，因其简单易操作被广泛应用于土壤理化性质及土壤属性的空间制图中[30]。
地理加权回归(GWR)：是一种适用于“空间非平稳性”地理现象的空间分析方法，为典型的局部回归模型。通过建立空间范围内每个点处的局部回归方程，来探索研究对象在某一尺度下的空间变化及相关驱动因素，并可用于对未来结果的预测[10]。可表示为：
| $ {y_i} = {\beta _0}\left( {{u_i}, {v_i}} \right) + \sum\limits_{k = 1}^{\text{m}} {{\beta _k}\left( {{u_i}, {v_i}} \right){x_{ik}}} + {\varepsilon _i} $ | (1) |
式中：yi是在位置i处的因变量值；xik(k=1，2，…，m)为位置i处的自变量值；(ui，vi)为样点i的坐标；β0((ui，vi)为截距项；βk(ui，vi)为采样点i上的第k个自变量x的回归系数，由自适应函数bi-square函数确定，用改进的赤池信息量准则(AICc)作为选择标准确定最优带宽；
随机森林(RF)：是一种强大的非参数方法，由多棵分类和回归树(CART)模型集成的机器学习算法，CART模型是一个相对复杂的模型，通过训练多个决策树来提高模型的准确性[23，31]。每个决策树使用从训练数据集中随机选择的样本进行训练，利用随机重采样技术bootstrap来创建决策树节点；多棵回归树集成随机森林，并选取所有决策树预测值的平均值作为最终结果。在RF模型中，决策树的数量ntree和每棵决策树节点分裂时选用的随机特征数量mtry是两个关键参数[15，19，23]。本研究设定mtry和ntree两参数分别为5和500。
随机森林–普通克里格(RF-OK)：是一种将随机森林和普通克里格结合的方法。在该方法中，随机森林输出的土壤属性预测值被认为是趋势项，训练集的属性实测值减去趋势项得到残差项；接着，针对残差项进行普通克里格插值；最后，将趋势项和进行插值的残差项相加，得到研究区随机森林–普通克里格混合模型土壤属性空间预测结果[22-24，32]。其输出结果可表示为：
| $ f\left( x \right) = m\left( x \right) + \varepsilon \left( x \right) $ | (2) |
式中：f(x)为因变量在点x处的预测值；m(x)为随机森林模型拟合的趋势项；ε(x)为残差普通克里格插值的误差项。
地统计学分析、地理加权回归及随机森林等建模过程均在R-4.2.3编程软件中，分别运用gstat、spgwr、randomForest和caret等包完成，并运用SPSS软件对数据进行统计分析，模型预测得到的土壤属性空间预测制图在ArcGIS10.8软件中完成。
1.5 模型精度评价与不确定性评估为了评估不同预测模型的性能，利用R语言caret包中的createDataPartition函数从152个样点中随机选取76% 的样点(116个)作为建模集参与建模，其余24% 的数据(36个)作为验证集(图 1)，用于评估不同模型的预测精度。建模集采用十折交叉验证，并对每种模型重复100次，取其平均值作为模型的预测结果。同时，逐像元计算重复迭代100次的模型预测值的标准差(SD)，并生成标准差图来评估模型预测的不确定性[32]。采用均方根误差(RMSE)、一致性相关系数(CCC)、决定系数(R2)和偏差(Bias)4个常用指标来验证模型的预测精度[23]。计算公式如下：
| $ {R^2} = 1 - \frac{{\sum\limits_{i = 1}^n {{{\left( {{{\hat y}_i} - {y_i}} \right)}^2}} }}{{\sum\limits_{i = 1}^n {{{\left( {{{\bar y}_i} - {y_i}} \right)}^2}} }} $ | (3) |
| $ {\text{RMSE}} = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {{{\left( {{{\hat y}_i} - {y_i}} \right)}^2}} } $ | (4) |
| $ {\text{CCC}} = \frac{{2\rho {\sigma _y}{\sigma _{\hat y}}}}{{\left( {\bar Y - \bar {\hat Y}} \right) + {\sigma _y} + {\sigma _{\hat y}}}} $ | (5) |
| $ {\text{Bias}} = \frac{1}{n}\sum\limits_{i = 1}^n {\left( {{y_i} - {{\hat y}_i}} \right)} $ | (6) |
式中：n为样点数量，i为第i个样点，yi为实测值，
友谊农场核心示范区总样本、建模集和验证集的土壤属性描述性统计特征如表 2所示。其中，示范区152个样点pH介于5.26 ~ 8.42，平均值为6.63；SOM含量的范围在19.41 ~ 109.17 g/kg，平均值为42.26 g/kg，但其标准差较大，说明示范区局部区域间土壤SOM含量具有不同程度波动；TN含量的范围在0.94 ~ 5.20 g/kg，平均值为1.94 g/kg。pH、SOM和TN的偏度值分别为0.36、1.72和1.76，峰度值分别为–1.12、6.48和5.80。经K-S正态性检验，pH(P=0.11)和SOM含量(P=0.16)均呈正态分布(P > 0.05)，而TN含量(P=0.03)不符合正态分布特征，因此对其进行了对数转换。变异系数(CV)的大小可反映土壤养分指标的变异程度[33]，示范区样点土壤pH、SOM和TN的变异系数分别为13.67%、29.50% 和31.98%，均属于中等程度的变异性。此外，验证集与建模集土壤样点的各项描述性统计特征相似(表 2)。
|
|
表 2 示范区土壤属性基本统计特征 Table 2 Basic statistical characteristics of soil properties in the demonstration area |
基于建模集和验证集数据，4种模型对示范区土壤pH、SOM和TN含量预测精度评价指标见表 3。总体来看，不同模型对不同属性的预测精度表现各异。对于建模集，4种模型对土壤pH、SOM和TN预测中，R2范围分别为0.60 ~ 0.81、0.43 ~ 0.71和0.28 ~ 0.71，表明模型对3种属性的解释率分别在60%、43% 和28% 以上；CCC值分别在0.75、0.62和0.43以上，说明4种模型对土壤属性具有不同程度的预测效果；RMSE值分别在0.41 ~ 0.57、7.06 ~ 10.87 g/kg和0.20 ~ 0.55 g/kg。综合来看，建模集中RF模型预测土壤pH(R2=0.81，CCC=0.87，RMSE=0.41)、SOM (R2=0.71，CCC=0.79，RMSE=7.06 g/kg)和TN(R2=0.71，CCC=0.80，RMSE=0.35 g/kg)效果最为突出。
|
|
表 3 示范区土壤属性预测模型精度比较 Table 3 Comparison of the accuracies of soil property prediction models in demonstration areas |
在验证集土壤pH预测中，OK模型(R2=0.84，CCC=0.86，RMSE=0.39)和RF-OK模型(R2=0.83，CCC=0.84，RMSE=0.41)预测性能均表现最优，其次是RF模型(R2=0.79，CCC=0.82，RMSE=0.45)，GWR模型预测性能最低(R2=0.65，CCC=0.78，RMSE= 0.56)。由于OK法只对单一指标进行插值，并不考虑环境变量对土壤属性的影响，其次模型验证集的RF-OK和OK法预测的R2仅差0.01、RMSE仅差0.02。因此，结合土壤pH验证集与建模集各评价指标及是否考虑环境变量等，综合考虑将预测残差作为误差项加入RF模型的RF-OK模型作为对预测土壤pH的最佳制图模型。土壤SOM预测中，RF-OK模型预测效果最突出(R2=0.72，CCC=0.68，RMSE=7.36 g/kg)，模型预测效果依次为RF-OK > OK > RF > GWR。土壤TN预测中，验证集和建模集RF模型预测效果均表现最佳(验证集：R2=0.59，CCC=0.68，RMSE= 0.36 g/kg)。在大多数情况下，3种土壤属性建模集和验证集的Bias均接近0，表明模型的拟合程度整体均表现较好。
综上，基于验证集模型的表现性能，不同模型对不同土壤属性的预测评价指标结果各有差异，其中RF-OK模型对示范区土壤pH和SOM含量的预测性能表现最佳，RF模型对土壤TN含量预测性能表现最佳。
示范区验证集土壤属性实测值与4种模型预测值关系散点图如图 3 ~ 图 5所示。由最优模型实测值与预测值的散点图可知，RF-OK模型得到的pH和SOM含量的实测值与预测值之间的R2分别达到0.83和0.72，RF模型得到的TN含量的实测值与预测值之间的R2为0.59，表明预测值与实测值拟合程度较好，表现出较强的相关性。其中pH的实测值与预测值的散点图相较SOM和TN含量的离散程度较为明显，但整体基本接近在1∶1对角线附近，与精度评价效果一致。虽然预测模型存在高值区被低估、低值区被高估的现象，但整体效果较好。
|
图 3 示范区验证集土壤pH实测值与4种模型预测值关系散点图 Fig. 3 Scatter plots of measured soil pH in the validation set in relation to the predicted values of the four models |
|
图 4 示范区验证集土壤SOM实测值与4种模型预测值关系散点图 Fig. 4 Scatter plots of measured SOM in the validation set in relation to the predicted values of the four models |
|
图 5 示范区验证集土壤TN实测值与4种模型预测值关系散点图 Fig. 5 Scatter plots of measured soil TN in the validation set in relation to the predicted values of the four models |
研究区3种土壤属性的空间分布表现出较强的空间异质性(图 6 ~ 图 8)，4种模型预测的土壤pH、SOM和TN空间分布的整体变化趋势基本一致，均表现出东北部高、西南部低的空间格局。图 9为基于最优模型所预测的3种属性的不确定性分布图。
|
图 6 示范区不同模型土壤pH空间分布图 Fig. 6 Spatial distributions of soil pH predicted by different models in the demonstration area |
|
图 7 示范区不同模型土壤SOM空间分布图 Fig. 7 Spatial distributions of SOM predicted by different models in the demonstration area |
|
图 8 示范区不同模型土壤TN空间分布图 Fig. 8 Spatial distributions of soil TN predicted by different models in the demonstration area |
|
图 9 基于最优模型示范区土壤属性空间分布不确定性图 Fig. 9 Uncertainty maps of spatial distributions of soil properties in the demonstration area based on the optimal model |
不同属性虽整体分布趋势基本一致，但局部分布仍存在差异。从最优模型预测的3种土壤属性空间分布和不确定性结果来看，示范区土壤pH表现出东北部主要呈中性土壤，而西南部主要为酸性土壤；其不确定性范围为0.01 ~ 0.58，平均值为0.16，标准差为0.05。土壤SOM和TN含量分布表现为东北部的含量水平明显高于西南部，不确定性范围分别为0.15 ~ 14.65 g/kg和0.01 ~ 1.04 g/kg，平均值分别为1.91 g/kg和0.15 g/kg，标准差分别为0.98 g/kg和0.09 g/kg。其中，pH的不确定性整体最低；而SOM和TN含量则表现出在示范区为湿地保护区且未布设样点的东北部和地势相对较高的西南部拥有较高的不确定性，而低不确定性相对集中在中部。总体而言，基于最优模型的土壤属性预测结果表现稳定且具有较低的不确定性。
土壤pH、SOM和TN是评价土壤肥力状况的重要指标。本研究根据第二次全国土壤普查土壤属性分级标准和前人研究[33]，将3种土壤属性指标划分为6个级别(表 4)。表 5展示了友谊农场核心示范区土壤属性的最优模型空间预测结果描述性统计。从表中可以看出，示范区土壤pH预测范围为5.41 ~ 8.15，平均值为6.64，土壤酸碱度等级为Ⅲ级标准，表明示范区整体呈现中性土壤，具有较高的肥料利用率；土壤SOM含量的预测范围为22.65 ~ 106.11 g/kg，平均值为43.95 g/kg，属于Ⅰ级标准，表明示范区表层土壤SOM含量整体上较高；土壤TN含量预测范围为1.21 ~ 3.79 g/kg，平均值为1.99 g/kg，属于Ⅱ级标准。从变异系数来看，示范区土壤pH属弱变异性，变异系数为9.49%；土壤SOM和TN含量变异系数分别为17.95% 和20.10%，均属中等程度变异性。总体而言，示范区的pH在东北部主要呈碱性或中性，在西南部主要呈酸性，土壤酸碱度主要处于Ⅱ ~ Ⅳ级；SOM和TN含量的空间分布呈现东北部含量水平整体高于西南部，并且养分指标丰缺度等级为Ⅰ ~ Ⅲ级。与示范区调查样点的属性统计结果相比，最优模型预测结果的均值在不同程度上有所增加，最值差距明显缩小，变异系数明显降低。由此进一步体现出不同的预测模型对不同土壤属性的预测性能存在差异。
|
|
表 4 土壤属性指标分级 Table 4 Classification of soil property indicators |
|
|
表 5 示范区土壤属性空间预测结果描述性统计 Table 5 Descriptive statistics of spatial prediction results of soil properties in the demonstration area |
示范区土壤属性的最优模型预测中，不同环境变量对空间预测结果的重要性存在差异，如图 10所示。在RF-OK模型预测土壤pH中，CNBL、DEM、MrVBF、Sand和VD等变量对模型的贡献率较高，特别是地形因子CNBL和DEM对示范区表层土壤pH的空间变异具有重要性影响。而在RF-OK模型预测土壤SOM含量中，Sand、MrVBF、RSP、CNBL和Clay等环境变量对模型的贡献率较高，尤其是土壤质地因子Sand和地形因子MrVBF对示范区土壤SOM含量的空间预测影响最为突出。另外，在RF模型预测土壤TN含量中，MrVBF、CNBL、CI、Silt和Sand等是对模型贡献率较高的变量，特别是地形因子MrVBF和CNBL对示范区土壤TN含量的空间预测影响最为显著。总体来说，地形因子和土壤因子是影响示范区土壤属性预测结果的重要环境变量。其中，地形因子在反映区域环境的同时，也体现了水文地理风貌对土壤属性分布的影响；土壤因子则直接反映了土壤自身性质对土壤属性分布的影响；而遥感数据指标则能够提供更多的区域信息，反映不同土地利用和覆盖类型对土壤属性分布的影响。由此可见，综合利用这些环境变量，能够提高土壤属性预测模型对区域土壤属性分布的解释能力和预测精度。
|
（pH和SOM采用RF-OK模型；TN采用RF模型） 图 10 基于最优模型土壤属性变量重要性 Fig. 10 Importance of soil property variables based on the optimal model |
土壤pH、SOM和TN是影响土壤肥力、作物生长状况和生态安全的重要土壤属性。本研究表明，在示范区内，不同模型对不同土壤属性的预测性能表现各异。其中，RF-OK模型在预测土壤pH和SOM方面表现较好，而RF模型则在预测土壤TN方面表现较好。该结果证实了RF模型能有效处理土壤属性与环境变量间复杂的非线性关系。同时，结合RF和OK两种方法的RF-OK模型能够进一步解释土壤-景观关系，减小RF回归残差的空间变异，表现出较好的预测性能。这与前人的研究结果相吻合。例如，蒙莉娜等[34]和赵彦锋等[35]运用GWR和RF等方法分别预测了渭干河-库车河绿洲土壤盐分和南阳市1 m土体SOC密度的空间分布，发现RF模型拥有较高的预测精度。此外，Mousavi等[23]和Zhang等[24]的研究验证了将RF模型的预测残差经克里格插值后与预测结果相结合的RF-OK模型能有效提高预测精度，并且表现出较好的预测性能。然而，也有研究指出RF-OK模型在预测不同土壤属性时并没有明显提高预测精度。例如，马重阳等[36]采用OK、RF和RF-OK对许昌市的SOM含量及其空间分布进行了预测，结果显示3种模型的性能表现没有明显差距。周洋等[32]利用RF和RF-OK对寻乌县的TN空间分布进行了预测，发现RF-OK模型并未提高其预测精度。这进一步说明目前还没有一种制图方法能够对各种土壤属性的预测都具有较好的性能，未来还需要更多的研究来改进和优化，以寻找更有效的方法来处理各种土壤属性。
由不同模型的土壤pH、SOM和TN的空间分布预测(图 6 ~ 图 8)可知，示范区3种土壤属性的整体空间分布趋势基本一致，表现为东北部含量水平明显高于西南部。这主要是由于：①地形因子是影响土壤属性空间分布的重要因素[24，37]。示范区西南部为丘陵地带，且靠近林地地区坡度相对较高，而东北部为低洼地，坡度较缓，整体地势由西南向东北倾斜。因此，示范区土壤pH随地势的升高而降低，而土壤SOM和TN含量在地势较低的东北低洼地积累。②土地利用与土壤属性的空间分布具有相关关系[38]。调查发现，示范区东北部主要为水田，而西南部主要是旱地，以种植大豆和玉米为主。水田要求地势较低，以利于灌溉和水稻种植。因为水田处于滞水状态，土壤SOM分解速率较低，有利于累积。研究还发现示范区水田的pH主要呈中性，而旱地则以酸性土壤为主，这与之前的研究结果[39-40]基本一致。所以，以上因素共同作用导致了示范区土壤属性含量水平在空间上的差异性，表现为东北部高于西南部的趋势。本研究对于示范区进一步优化土壤管理和合理利用土地具有重要的指导意义。
此外，本研究在示范区土壤属性预测中存在一定的局限性，这些限制与所考虑的环境变量有关。由于示范区地势平坦且范围较小，导致降水量与温度等气候因子对其土壤的影响表现相对均一，因此在本研究中没有考虑[5]。然而，有学者指出，引入基于遥感获取的地表动态反馈信息可以有效地揭示平原或地形平缓地区土壤的空间差异[41-42]。另外，示范区作为我国北方旱粮现代化生产试验基地，土壤受到人为因素的影响非常严重，如施肥、轮作、有机肥和秸秆管理等人类活动对土壤属性含量具有重要影响[43]。因此，在选择示范区环境变量时，除了考虑易获取和相对廉价的常用辅助变量外，还应考虑引入地表温度以及人类活动等多源辅助变量，将有助于更准确地预测示范区土壤属性，进而可为北方旱粮现代化生产提供有效的支持。
4 结论1) 友谊农场核心示范区pH介于5.26 ~ 8.42，均值为6.63；土壤SOM含量范围为19.41 ~ 109.17 g/kg，均值为42.26 g/kg；TN含量为0.94 ~ 5.20 g/kg，均值为1.94 g/kg。3种属性变异系数分别为13.67%、29.50% 和31.98%，均属中等程度的空间变异。
2) 友谊农场核心示范区不同模型对不同土壤属性的预测精度表现各异，RF-OK模型对示范区土壤pH和SOM含量的预测性能表现最佳，RF模型对土壤TN含量预测性能表现最佳。
3) 友谊农场核心示范区土壤pH、SOM含量和TN含量整体空间分布预测趋势基本一致，其中，土壤pH在示范区东北部主要以中性土壤为主，碱性土壤次之，西南部主要呈酸性土壤；土壤SOM含量和TN含量空间分布趋势表现为东北部低洼区含量水平明显高于西南部丘陵区。
| [1] |
韩长赋. 加强东北黑土地保护推进农业绿色发展[N]. 人民日报, 2018-02-05 (007).
(  0) |
| [2] |
韩晓增, 邹文秀. 东北黑土地保护利用研究足迹与科技研发展望[J]. 土壤学报, 2021, 58(6): 1341-1358 ( 0) |
| [3] |
李保国, 刘忠, 黄峰, 等. 巩固黑土地粮仓保障国家粮食安全[J]. 中国科学院院刊, 2021, 36(10): 1184-1193 ( 0) |
| [4] |
张江周, 李奕赞, 李颖, 等. 土壤健康指标体系与评价方法研究进展[J]. 土壤学报, 2022, 59(3): 603-616 ( 0) |
| [5] |
朱阿兴, 杨琳, 樊乃卿, 等. 数字土壤制图研究综述与展望[J]. 地理科学进展, 2018, 37(1): 66-78 ( 0) |
| [6] |
徐尚平, 陶澍, 曹军. 内蒙古土壤pH值、粘粒和有机质含量的空间结构特征[J]. 土壤通报, 2001, 32(4): 145-148 ( 0) |
| [7] |
Duan L X, Li Z W, Xie H X, et al. Large-scale spatial variability of eight soil chemical properties within paddy fields[J]. Catena, 2020, 188: 104350 DOI:10.1016/j.catena.2019.104350 ( 0) |
| [8] |
张素梅, 王宗明, 张柏, 等. 利用地形和遥感数据预测土壤养分空间分布[J]. 农业工程学报, 2010, 26(5): 188-194 ( 0) |
| [9] |
杨顺华, 张海涛, 郭龙, 等. 基于回归和地理加权回归Kriging的土壤有机质空间插值[J]. 应用生态学报, 2015, 26(6): 1649-1656 ( 0) |
| [10] |
赵明松, 刘斌寅, 卢宏亮, 等. 基于地理加权回归的地形平缓区土壤有机质空间建模[J]. 农业工程学报, 2019, 35(20): 102-110 ( 0) |
| [11] |
罗梅, 郭龙, 张海涛, 等. 基于环境变量的中国土壤有机碳空间分布特征[J]. 土壤学报, 2020, 57(1): 48-59 ( 0) |
| [12] |
Costa E M, Tassinari W S, Pinheiro H S K, et al. Mapping soil organic carbon and organic matter fractions by geographically weighted regression[J]. Journal of Environmental Quality, 2018, 47(4): 718-725 DOI:10.2134/jeq2017.04.0178 ( 0) |
| [13] |
Mousavi S R, Sarmadian F, Omid M, et al. Three-dimensional mapping of soil organic carbon using soil and environmental covariates in an arid and semi-arid region of Iran[J]. Measurement, 2022, 201: 111706 DOI:10.1016/j.measurement.2022.111706 ( 0) |
| [14] |
蒙莉娜, 丁建丽, 张振华. 耦合植被物候的新疆绿洲土壤盐分信息挖掘[J]. 土壤, 2022, 54(3): 629-636 DOI:10.13758/j.cnki.tr.2022.03.025 ( 0) |
| [15] |
Zeraatpisheh M, Ayoubi S, Jafari A, et al. Digital mapping of soil properties using multiple machine learning in a semi-arid region, central Iran[J]. Geoderma, 2019, 338: 445-452 DOI:10.1016/j.geoderma.2018.09.006 ( 0) |
| [16] |
王欣怡, 王昌昆, 马海艺, 等. 基于双时相卫星遥感光谱指数估算土壤有机质含量[J]. 土壤, 2023, 55(5): 1106-1113 DOI:10.13758/j.cnki.tr.2023.05.021 ( 0) |
| [17] |
Lu Q K, Tian S, Wei L F. Digital mapping of soil pH and carbonates at the European scale using environmental variables and machine learning[J]. Science of the Total Environment, 2023, 856: 159171 DOI:10.1016/j.scitotenv.2022.159171 ( 0) |
| [18] |
Yang H, Lim H, Moon H, et al. Simple optimal sampling algorithm to strengthen digital soil mapping using the spatial distribution of machine learning predictive uncertainty: A case study for field capacity prediction[J]. Land, 2022, 11(11): 2098 DOI:10.3390/land11112098 ( 0) |
| [19] |
Liu F, Wu H Y, Zhao Y G, et al. Mapping high resolution national soil information grids of China[J]. Science Bulletin, 2022, 67(3): 328-340 DOI:10.1016/j.scib.2021.10.013 ( 0) |
| [20] |
Song X D, Yang F, Ju B, et al. The influence of the conversion of grassland to cropland on changes in soil organic carbon and total nitrogen stocks in the Songnen Plain of Northeast China[J]. Catena, 2018, 171: 588-601 DOI:10.1016/j.catena.2018.07.045 ( 0) |
| [21] |
庞龙辉, 刘峰, 赵霞, 等. 青海省表层土壤属性数字制图[J]. 土壤通报, 2019, 50(3): 505-513 ( 0) |
| [22] |
Guo P T, Li M F, Luo W, et al. Digital mapping of soil organic matter for rubber plantation at regional scale: An application of random forest plus residuals Kriging approach[J]. Geoderma, 2015, 237: 49-59 ( 0) |
| [23] |
Mousavi A, Karimi A, Maleki S, et al. Digital mapping of selected soil properties using machine learning and geostatistical techniques in Mashhad plain, northeastern Iran[J]. Environmental Earth Sciences, 2023, 82(9): 234 DOI:10.1007/s12665-023-10919-x ( 0) |
| [24] |
Zhang W C, Wan H S, Zhou M H, et al. Soil total and organic carbon mapping and uncertainty analysis using machine learning techniques[J]. Ecological Indicators, 2022, 143: 109420 DOI:10.1016/j.ecolind.2022.109420 ( 0) |
| [25] |
鲍士旦. 土壤农化分析(3版)[M].
中国农业出版社, 北京, 2000: 1-495
( 0) |
| [26] |
Minasny B, McBratney A B, Malone B P, et al. Digital mapping of soil carbon[M]//Advances in Agronomy. Amsterdam: Elsevier, 2013: 1–47.
( 0) |
| [27] |
Song X D, Wu H Y, Ju B, et al. Pedoclimatic zone-based three-dimensional soil organic carbon mapping in China[J]. Geoderma, 2020, 363: 114145 DOI:10.1016/j.geoderma.2019.114145 ( 0) |
| [28] |
张振华, 丁建丽, 王敬哲, 等. 集成土壤-环境关系与机器学习的干旱区土壤属性数字制图[J]. 中国农业科学, 2020, 53(3): 563-573 ( 0) |
| [29] |
马冉. 流域尺度土壤特性空间分布及影响因素研究——以三峡库区草堂河流域为例[D]. 重庆: 西南大学, 2019.
( 0) |
| [30] |
杜挺, 杨联安, 张泉, 等. 县域土壤养分协同克里格和普通克里格空间插值预测比较——以陕西省蓝田县为例[J]. 陕西师范大学学报(自然科学版), 2013, 41(4): 85-89 ( 0) |
| [31] |
Breiman L. Random forests[J]. Machine Learning, 2001, 45: 5-32 DOI:10.1023/A:1010933404324 ( 0) |
| [32] |
周洋, 赵小敏, 郭熙. 基于多源辅助变量和随机森林模型的表层土壤全氮分布预测[J]. 土壤学报, 2022, 59(2): 451-460 ( 0) |
| [33] |
王世豪, 徐新良, 黄麟, 等. 1980s—2010s东北农田土壤养分时空变化特征[J]. 应用生态学报, 2023, 34(4): 865-875 ( 0) |
| [34] |
蒙莉娜, 丁建丽, 王敬哲, 等. 基于环境变量的渭干河–库车河绿洲土壤盐分空间分布[J]. 农业工程学报, 2020, 36(1): 175-181 ( 0) |
| [35] |
赵彦锋, 李豪杰, 陈杰, 等. 基于土壤变异解释力的几种土壤制图方法的对比研究——以南阳市1 m土体土壤有机碳密度制图为例[J]. 土壤学报, 2018, 55(1): 43-53 ( 0) |
| [36] |
马重阳, 孙越琦, 巫振富, 等. 基于不同模型的区域尺度耕地表层土壤有机质空间分布预测[J]. 土壤通报, 2021, 52(6): 1261-1272 ( 0) |
| [37] |
马冉, 刘洪斌, 武伟. 三峡库区草堂河流域土壤pH空间分布预测制图[J]. 长江流域资源与环境, 2019, 28(3): 691-699 ( 0) |
| [38] |
陈金萍, 李奕, 李祥余, 等. 不同农业土地利用方式对土壤养分含量分布影响及养分等级评价[J]. 萍乡学院学报, 2022, 39(3): 111-116 ( 0) |
| [39] |
Li Y T, Xie X N, Zhu Z J, et al. Land use driven change in soil organic carbon affects soil microbial community assembly in the riparian of Three Gorges Reservoir Region[J]. Applied Soil Ecology, 2022, 176: 104467 ( 0) |
| [40] |
Wei L, Ge T D, Zhu Z K, et al. Comparing carbon and nitrogen stocks in paddy and upland soils: Accumulation, stabilization mechanisms, and environmental drivers[J]. Geoderma, 2021, 398: 115121 ( 0) |
| [41] |
Liu F, Geng X Y, Zhu A X, et al. Soil texture mapping over low relief areas using land surface feedback dynamic patterns extracted from MODIS[J]. Geoderma, 2012, 171: 44-52 ( 0) |
| [42] |
Zeng C Y, Zhu A X, Liu F, et al. The impact of rainfall magnitude on the performance of digital soil mapping over low-relief areas using a land surface dynamic feedback method[J]. Ecological Indicators, 2017, 72: 297-309 ( 0) |
| [43] |
韩晓增, 李娜. 中国东北黑土地研究进展与展望[J]. 地理科学, 2018, 38(7): 1032-1041 ( 0) |
2. School of Geomatics, Anhui University of Science and Technology, Huainan, Anhui 232001, China;
3. University of Chinese Academy of Sciences, Beijing 100049, China