2024, Vol. 61
表1(Table 1)
2. 中国科学院大学, 北京 100049;
3. 黑龙江省黑土保护利用研究院, 哈尔滨 150086
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. Heilongjiang Academy of Black Soil Conservation and Utilization, Harbin 150086, China
东北地区土壤黑土层特指富含腐殖质的暗沃表土层[1]，是植被生长根基及营养重要来源，其厚度成为土壤质量及资源禀赋的重要标志[2-6]。
土壤剖面挖掘是黑土层厚度调查的经典方法。早在20世纪二三十年代起，国内外学者利用土壤剖面挖掘调查我国的黑土层厚度[7-11]。1925年和1926年俄国Gordeef[7]先后调查了中东铁路沿线土壤，将典型剖面分为腐殖质层和底土层，其中腐殖质层包括A层和B层。1935年美国土壤学家Pendleton和Hou[8]与中国学者在哈尔滨、海伦、安达、双城等地进行剖面调查，得出典型剖面表土层厚度为25 ~ 100 cm，一般厚度在40 ~ 50 cm。1958—1960年，黑龙江省农业厅[9]在全国第一次土壤普查工作中以耕地为主，典型剖面记录的黑土层厚度一般在30 ~ 50 cm。1979—1994年，全国第二次土壤普查[10-11]中，以典型剖面A层为黑土层，《中国土种志》[12]等686个典型剖面黑土层厚度一般小于60 cm。
近些年，通过土壤剖面分析黑土层变薄问题也得到关注[13-15]。陆继龙[13]对比分析了黑龙江省和吉林省第一次土壤普查（简称“一普”）与第二次土壤普查（简称“二普”）的土壤剖面统计资料，得出黑土层厚度在此期间每年减少0.4 ~ 0.5 cm。沈波等[14]的研究表明，东北黑土区从黑土地开垦之初到开垦20年（“一普”）、再到开垦40年（“二普”），黑土层厚度每年减少0.7 ~ 1 cm。张之一[3]对比分析了黑龙江省北部地区和内蒙古呼伦贝尔盟“一普”与“二普”黑土层厚度剖面统计资料，指出并非所有黑土层厚度均变薄，变薄原因主要是坡耕地面蚀和耕地土壤由松到紧的压实过程。刘凯等[15]通过61个“二普”剖面点的对应点位剖面调查，点对点比较得出东北黑土区的黑土层厚度平均每年减少0.32 cm。
然而，野外土壤剖面调查方法，不仅存在调查样点类型代表性问题，也因人力消耗及成本费用高导致土壤剖面样本量少，严重影响了黑土层厚度及其变化调查分析的准确性。现有调查报告也多基于点位数据统计，缺乏连续性空间分布特征分析。仅有的某小流域黑土层厚度空间预测[16]，区域覆盖范围较小（1.39 km2）也不能满足需求。因此，针对土壤剖面调查样本量少的问题，需建立黑土层厚度快速、简易调查获取方法；针对连续空间分析不足问题，需建立高性能黑土层厚度空间预测方法。
目前土层厚度快速调查方法也在不断尝试，如土钻、雷达探测方法等。2013年和2014年张福平等[17]在甘肃省敦煌市利用土钻法调查了124个样点1 m深土层厚度信息。张世文等[18]采用500 MHZ探地雷达调查北京市昌平区耕地土层厚度，得出土层厚度最佳探测深度小于45 cm，否则探测精度得不到保证。探地雷达接受雷达信号波，经过模型处理可以快速获取土层厚度信息，一定深度范围内模型精度高，但显然，目前雷达探测方法多数仍处于试验阶段，探测精度和稳定性问题依然没有得到解决，大范围区域应用也受到限制。而土钻快速获取方法，也还存在表层土壤松散滑落以及深层土壤压缩等问题；相对于土壤剖面调查断面，土钻调查穴位空间极小，缺乏多穴位重复调查，样点微域空间变异性势必影响调查精度。相较于探地雷达等方法存在问题，土钻法更易于满足调查急迫需求，但其调查规范性等问题需得到彻底解决。
黑土层厚度空间预测方法主要有经典统计[19-20]、地统计学[21-22]和机器学习[23-24]。经典统计方法往往判定变量之间需具有较好的线性相关性，但对于变量和表达式的选择有一定主观性，不能很好拟合非线性数据[25]；地统计学方法更侧重土壤空间自相关，忽视土壤与环境之间的关系，较少考虑成土因素影响，在大面积复杂景观尺度下的空间插值结果的真实性较弱[25]；而机器学习在处理复杂景观以及非线性变量方面具有优势。黑土层厚度土钻快速调查方法，虽可满足机器学习对大样本量的数据要求，但对于大面积复杂景观地区，受限于道路通达性等条件，所能采集的样本数量往往还不够充分。在样本量受限的条件下，随机森林方法得到改进，产生了与分位数相结合的分位数回归森林方法[26]，但这种模型预测的性能及稳定性，是否随样本数据集不同而发生改变还需进一步论证。
因此，考虑到野外土壤剖面调查方法费时费力、探地雷达精度未解决等问题，土钻法更易于满足需求。为解决土钻法调查规范性问题，同时考虑微域空间变异性，本文将构建调查样点黑土层厚度多穴位土钻快速获取方法，通过优化回归森林模型参数开展空间预测，分析不同样本对优化模型预测性能及稳定性的影响，为黑土层厚度高精度快速调查与预测提供新经验、新途径。
1 材料与方法 1.1 研究区概况黑龙江省地处中国东北部（43°26′ ~ 53°33′ N，121°11′ ~ 135°05′ E），辖区总面积47.3万km2，地势西北、北部和东南部高，东北、西南部低，属寒温带与温带大陆性季风气候，年均温4.2 ℃，年均降雨量608.5 mm[27]。土壤类型主要有黑土、黑钙土、栗钙土、草甸土、暗棕壤、白浆土、沼泽土、泥炭土、水稻土等；成土母质主要有黄土状黏质土、坡积物、冲积物等。研究区土地利用可概括为“五山一水一草三分田”，耕地总面积为17.19万km2，约占全国耕地总面积的13%，集中于黑土、黑钙土、草甸土、白浆土、暗棕壤等土壤类型。
1.2 样点布设与野外调查充分考虑土地利用、土壤及母质类型与地形因素，并兼顾坡耕地完整性，于2022年夏季布设和调查311个耕地土壤及15个林草地对照样点（图 1）。其中，旱地样点283个、水田样点28个；不同土壤类型样点分别为黑土101个、草甸土74个、暗棕壤71个、黑钙土28个、沼泽土23个、白浆土20个、水稻土8个、棕色针叶林土1个。
|
图 1 研究区位置与采样点分布 Fig. 1 Study area location and sampling point distribution |
在每个调查样点2 m×3 m的空间范围内，随机选择3个土壤调查穴位A、B、C，采用“浅层挖掘与深层土钻”相结合方法开展黑土层厚度快速调查。在每个穴位，首先挖掘0 ~ 40 cm简易土壤剖面（宽40 cm×长50 cm），如果黑土层厚度≤40 cm，可直接判断测量黑土层厚度，避免采用土钻方法因表层土壤松散沿土钻凹槽上窜或下滑带来的观测误差。当黑土层厚度≥40 cm时，则在相应穴位中心，采用长100 cm、直径3 cm的荷兰Eijkelkamp钻（0401SB）取土，判定黑土层厚度；为避免土壤被挤压带来观测误差，按20 cm等间距分别取40 ~ 100 cm完整土样。若黑土层厚度≥100 cm时，以100 cm计。
野外调查中记录各调查样点经纬度、坡度、坡向、坡位、土地利用类型、地貌类型、土壤类型、母质类型、水土保持措施以及各穴位的表层土壤颜色、土壤质地等基本信息，同时采集表层（0 ~ 20 cm）、亚表层（20 ~ 40 cm）3个穴位的混合土壤样品各2 kg，带回实验室测定其理化性质。以3个穴位黑土层厚度的平均值作为该调查样点的黑土层厚度（图 2）。
|
图 2 野外调查样点各观测穴位分布示意图 Fig. 2 Map of the distribution of observation points in the field survey sample sites |
考虑野外调查可操作性，参考美国土壤系统分类[1]及中国土壤系统分类[28]关于暗沃表土层判定标准，本次黑土层调查判定指标及标准为：①土壤颜色润态明度≤3.5，润态彩度≤3.5；②土壤结构呈粒状、小角块状和小亚角块状；③石块体积占比小于1/2[29]。在实际操作中，先判定样点是否满足土壤结构和石块比例条件，再用芒塞尔比色卡比对土壤颜色。为减少土壤颜色判定上的主观误差，每个穴位均由3~4人集体判定。
2023年7月在黑龙江省嫩江市示范区，又布设31个耕地土壤样点，涉及土壤类型主要有黑土（11个）、草甸土（7个）、暗棕壤（12个）、沼泽土（1个），分别同时采用剖面法和上述多穴位土钻快速调查法获取各样点黑土层厚度数据，并以剖面调查黑土层厚度数据为真实值，以验证多穴位黑土层厚度快速勘察与预测方法的精确度和可行性。
1.3 样点黑土层厚度观测误差由于土壤调查样点不可避免地存在微域空间上的变异性，单个穴位调查难以获取调查样点的真实值，通过在一定空间（标准剖面挖掘空间）范围内测定多个观测穴位的黑土层厚度，取其平均值作为该调查样点的观测值。观测穴位重复越多，平均值越接近真实值，但样点平均黑土层厚度计算误差依然存在。考虑到重复穴位越多采样所消耗的成本也就越大，在每个调查点设立了3个重复调查穴位（A、B、C）。设YT为样点黑土层厚度的真实值，YS为黑土层厚度观测值，观测误差方差[30]：
| δ2(s)=var(YS−YT) | (1) |
依据第i个样点的3个穴位黑土层厚度实际观测值获取该样点平均值
| δ(s)=√∑ni=1(YS(i)−ˉYT(i))2n | (2) |
遥感因子、地形因子、气候因子、土壤因子、人为因子共同构成黑土层厚度预测的环境协变量指标体系（表 1）。其中，DEM数据、2022—2023年遥感影像数据从USGS获取，土壤属性数据从世界土壤数据库（HWSD）获取，气候数据从国家地球系统科学数据中心（http://www.geodata.cn）获取并计算年平均值，土地利用数据从GlobeLand30（http://www.globallandcover.com）和资源环境科学与数据中心（https://www.resdc.cn）下载。将2022年和2023年卫星影像通过中值合成后提取20个遥感因子指数值（表 1）。
|
表 1 预测黑土层厚度的环境协变量 Table 1 Environmental covariates for predicting black soil layer thickness |
随机森林模型（Random Forest Algorithm，RF）是一种用于分类和回归的机器学习方法，是随机分类与回归树集合的组合。利用Boruta特征指标选择，将原始75个变量的特征矩阵随机生成影子特征矩阵，并与原始变量特征矩阵组合成新的特征矩阵训练RF模型，计算每个变量的相对重要性。选取影子特征矩阵相对重要性的最大值，当原始变量特征矩阵相对重要性大于这个最大值则被标记为重要变量[31]。
对随机森林模型进行训练时，仅利用2/3样本作为训练数据，此外1/3作为袋外数据（out-of-bag，OOB）检验预测精度，以评估随机森林性能[32]。将袋外数据重复运行了100次，计算袋外数据预测误差对模型进行校准。
随机森林模型校准主要对三项模型参数进行调整。其中，参数nodesize为拆分二叉树的最小节点大小；mtry是指定节点用于二叉树的变量个数；ntree是形成森林的树木数量，默认值500足以达到稳定预测[33]。例如以穴位均值作为样点观测值，针对mtry（2 ~ 15）和nodesize（3 ~ 10）参数值范围，设置构成112个参数组合方案，并采用10折交叉验证方法，选择RMSE值最小的参数组合方案，最终确定模型最优参数值。
单个样点或栅格像元的随机森林预测值为其所有决策树加权线性组合[26]的平均值，表示为：
| ˆY(i)=∑nm=1wm(i)×ym(i),(i=1,2,3...) | (3) |
式中，
在随机森林预测过程中，基于样点均值以及各穴位A、B、C观测值分别进行特征指标选择、模型评估、参数优化调整，获得最优随机森林模型进行黑土层厚度空间预测；比较四种样本预测结果，分析不同穴位随机观测对环境变量重要性、预测空间分布、模型精度的影响。
1.6 模型预测精度检验及不确定性分析分别采用10折（T1）、留一（T2）以及5倍10折（T3）三种交叉验证方法，同时检验四种样本数据集下随机森林模型精度，评价随机森林回归预测潜能；采用的评价指标为平均绝对误差（MAE）、均方根误差（RMSE）和决定系数（R2）[34]。
| R2=∑ni=1(Yi−ˉYi)2∑ni=1(Yi−ˆYi)2 | (4) |
| MAE=∑ni=1(Yi−ˆYi)n | (5) |
| RMSE = √∑ni=1(Yi−ˆYi)2n | (6) |
式中，
利用同种交叉验证方法获得的不同样本集MAE、RMSE、R2值之间比较，分析样点多穴位及单穴位对模型预测性能的影响。利用不同交叉验证算法对同一样本数据集的MAE、RMSE、R2值差异，分析模型预测潜能；当这3种评价指标值变化较小时，表明模型稳定性较高，具备较好空间预测潜能。
空间预测不确定性采用定量分析方法。应用分位数回归森林模型获取研究区100 m空间分辨率各栅格单元5%、50%、95%分位数的黑土层厚度预测值；将每个像元90% 预测区间设定为5%与95%分位数黑土层厚度预测值之间的差值，其与50%分位数黑土层厚度预测值的比值，即为该像元不确定性指数[35]。
| UC(i)=ˆY0.95(i)−ˆY0.05(i)ˆY0.5(i), (i=1,2,3...) | (7) |
式中，
同时，利用分位数回归森林计算的预测区间覆盖率（PICP），对栅格单元不确定性分析的准确性进行评估[36]。例如，当预测区间为90%，记为PICP90，表示在5%与95%分位数之间验证数据的覆盖比例。将观测值作为验证数据，设黑土层厚度观测值为TS，α为预测区间，
| PICP=1n∑ni=1Ci | (8) |
式中，n为样点数。
为直观反映模型预测不确定性分析的准确性水平，绘制了不同观测穴位样点预测值不确定性分析的精度图。其中，横、纵坐标分别代表预测区间α及其覆盖率PICP。当精度图的点线越靠近1︰1线，则说明空间预测不确定性评估良好，高于或者低于1︰1线则表示不确定性分别被高估或者低估。
2 结果 2.1 黑土层厚度样点统计特征及勘察精度检验各调查样点的观测穴位A、B、C的黑土层厚度平均值分别为55.52、55.29、56.86 cm，调查样点平均值为55.56 cm（表 2）。样点平均黑土层厚度的标准差最小，其数据离散程度均低于观测穴位A、B、C。样点观测穴位A、B、C及其平均的变异系数均介于10% ~ 100%之间，属于中等变异性；偏度均大于0，说明数据呈正偏态，数据分布主体集中在左侧；峰度均小于0，说明数据呈矮胖分布。此外，样点观测穴位A、B、C在95%置信度下的均值误差均高于其样点平均的均值误差。不同土壤类型下，黑土的黑土层厚度测定结果最高为68.06 cm，白浆土的黑土层厚度测定结果最低为34.40 cm（表 3）。
|
|
表 2 土壤样点黑土层厚度描述性统计 Table 2 Descriptive statistics of black soil layer thickness at soil sample sites(n=357) |
|
|
表 3 主要土壤类型土壤样品黑土层厚度快速调查结果统计特征 Table 3 Statistical characteristics of the results of a rapid survey of black soil layer thickness in soil samples of major soil types |
本调查获得的各样点平均黑土层厚度小于Zhang等[16]在小流域平原地带调查的平均黑土层厚度（80 cm），但高于“二普”调查的平均厚度（40 cm）[12]，原因在于研究区以及黑土层厚度判定方法均存在差异。本研究基于暗沃表层分类标准判定黑土层，而“二普”《中国土种志》确定土壤发生分类A层为黑土层，前者已被国内外广泛接收和认可[1，28]。
观测穴位A、B、C的观测误差
多穴位黑土层厚度快速勘察方法与剖面法获取的黑土层厚度基本一致（y=0.9517x+ 3.1252），两者间的皮尔森相关性高达0.9744，决定系数R2为0.9494，说明该方法精确度较高，具有可行性，可以替代剖面调查方法[37]。
2.2 随机森林模型参数优化与检验以穴位均值作为样点观测值，采用boruta方法优选出15个用于黑土层厚度预测的环境因子协变量，由于随机森林模型预测受共线性影响较小，较多的环境协变量可提升预测能力[30]。其中，遥感因子5个、地形因子7个、气候因子2个、土壤因子1个。人为因子未被选中，是因为耕地主要集中在黑土、黑钙土、草甸土和白浆土4种土壤类型上（图 1），土地利用类型信息已经被土壤类型所解释。此外，观测穴位A、B、C分别获取了17、15、13个环境因子协变量。
以穴位均值作为样点观测值，袋外数据RMSE介于25.25 ~ 27.96 cm之间，平均值为26.71 cm；MAE介于21.28 ~ 24.36 cm之间，平均值为23.09 cm。袋外数据RMSE和MAE的平均值占研究区黑土层厚度平均值的48.11%和41.59%；Zhang等[16]构建随机森林绘制的1.39 km2东北小流域黑土层厚度图，袋外数据RMSE占流域黑土层平均厚度46.85%，均说明随机森林对测试数据的模型构建表现良好。观测穴位A、B、C袋外数据RMSE平均值分别为26.87 cm、27.37 cm、26.92 cm，MAE平均值分别为23.19 cm、23.72 cm、23.10 cm。
以穴位均值作为样点观测值，采用的随机森林最优模型参数值（mtry = 15和nodesize = 4）与默认模型参数（mtry = 2和nodesize = 5）相比，最优模型参数值将模型性能R2、RMSE、MAE分别提高了5.22%、4.60%、5.42%，表明构建模型时参数调优，能够提高模型结构的预测性能。在最优模型中，10倍交叉验证R2为0.60，即对黑土层厚度空间变化解释力达到60%，略高于Liu等[25]对祁连山地区14万平方千米的土壤深度的解释力（R2=0.59）。此外，观测穴位A、B、C随机森林最优模型参数值分别为A（mtry = 17和nodesize = 3）、B（mtry = 13和nodesize = 4）、C（mtry = 13和nodesize = 3）。
2.3 黑土层厚度空间分布及其不确定性土壤样点黑土层厚度预测值按土壤类型统计结果表明，暗棕壤、白浆土、黑钙土、棕色针叶林土黑土层厚度平均预测值较测定值略高，草甸土、黑土、水稻土、沼泽土的平均预测值较测定值略低，但总体差异不大（表 3—4）。不同土壤类型下，黑土的黑土层厚度平均预测值最高，白浆土的黑土层厚度平均预测值最低，与实际观测结果具有一致性。
|
|
表 4 主要土壤类型样点黑土层厚度空间预测结果统计特征 Table 4 Statistical characteristics of the results of spatial prediction of black soil layer thickness at sample sites of major soil types |
张之一 [3]在全国第二次土壤普查中，将黑土层厚度≤30 cm定义为薄层黑土，介于30 ~ 60 cm之间为中层黑土，> 60 cm为厚层黑土。研究区以穴位均值空间预测的黑土层厚度均值为51.73 cm，其中耕地平均黑土层厚度为53.42 cm，分布格局与其地形地貌特征吻合（图 3a），平原地区黑土层厚度较高，山地黑土层厚度较低。厚层黑土主要分布在松嫩平原和三江平原，中层黑土分布在大兴安岭、小兴安岭以及东南山地，薄层黑土分布在小兴安岭东南部地区以及东南山地西部地区（图 3d）。嫩江市处于毗邻内蒙古的嫩江东岸平原，土地肥沃，黑土地资源相对丰富，黑土层厚度较高。此外，松花江流域地势相对较低，加上河流汇集泥沙堆积，黑土层厚度也相对较高。
|
图 3 黑土层厚度预测及不确定性空间分布 Fig. 3 Spatial distribution of black soil layer thickness prediction and uncertainty |
黑土层厚度分布的局部细节表明，从山顶到山底黑土层厚度逐步变厚（图 3b）；从山地到平原，由南至北黑土层厚度也逐步变厚（图 3c）。
黑土层厚度空间分布不确定性指数UC平均值仅为2.13，远低于不确定性指数中位数4，表明以穴位均值预测的不确定性指数总体较低（图 3e）。高不确定性指数主要分布在西北大兴安岭、北部小兴安岭以及东南山地，一方面由于山地黑土层厚度较薄，空间变异性高，预测误差相对较高；另一方面是由于山地可达性较差，样点数目稀少。低不确定性指数分布在三江平原、松嫩平原、松花江流域以及嫩江地区，这些区域地势开阔可达性好，耕地分布集中连片，样点调查密度较高。结合黑土层厚度不确定性空间信息，黑土层厚度空间分布特征更为清晰。
3 讨论 3.1 样点观测穴位对环境协变量重要性影响利用样点观测值，通过boruta特征提取方法获取15个环境协变量，不同环境变量对黑土层厚度预测的重要性有明显差异（图 4a）。而基于穴位A、B、C观测分别获取17、15、13个环境协变量（图 4b~图 4d），其中共同的环境变量有13个，重要性也因穴位不同发生变化。
|
图 4 环境协变量的重要性排名：（a）样点平均；（b）观测穴位A；（c）观测穴位B；（d）观测穴位C Fig. 4 Importance ranking of environmental covariates: (a)sample point average; (b)observation point A; (c)observation point B; (d)observation point C |
以穴位均值作为样点观测值，黑土层厚度最重要的协变量有GNDI、PRE、MRV、BAI、CND、MCA、SOIL。其中，GNDI、BAI代表遥感因子，分别反映了黑龙江省植被指数、耕地和非耕地的分布状况，BAI燃烧指数代表秸秆燃烧程度反映种植作物类型、生物量[38]。燃烧明显的耕地区域黑土层较厚；未燃烧区的林地黑土层厚度较薄，Zhao等[38-39]研究很好论证了BAI等遥感指数是黑土层厚度主要协变量。PRE代表气候因子，决定了黑土成土过程的水热条件，不仅影响土壤母质风化，也控制着植物和微生物生长，影响有机质合成、分解和迁移[40]，进一步影响黑土层厚度。SOIL代表土壤因子，反映土壤类型属性及其剖面结构特征，与黑土层厚度直接相关，黑土的黑土层较厚、暗棕壤的黑土层较薄[41]。MRV、CND、MCA代表地形因子，主要改变了黑土层成土母质、水分和热量等条件再分配，水分、母质在较低地形部位积累，在较高地形部位流失，热量的再分配表现为南坡大于北坡[25，40]。
遥感因子、地形因子、气候因子、土壤因子对黑土层厚度的重要性值分别为37.37%、38.52%、17.97%、6.14%，各因子单个变量重要性平均值分别为7.47%、5.50%、8.99%、6.14%，说明地形因子与单变量平均气候因子解释黑土层厚度能力最高（图 4a）。但各调查样点不同观测穴位的空间变异性对变量重要性值也产生了明显影响（图 4b~图 4d），穴位A的遥感因子、地形因子、气候因子、土壤因子的重要性值变化量分别为2.12%、8.14%、7.13%、1.11%；穴位B的各因子变化量分别为13.67%、11.19%、8.94%、6.46%；穴位C各因子变化量分别为12.45%、7.7%、0.92%、5.67%。其中，穴位B产生的变量重要性值变化最大，这与观测穴位随机选择有关。
3.2 样点观测穴位对黑土层厚度空间预测分布的影响基于样点不同穴位A、B、C观测值与平均值预测的黑土层厚度空间分布（图 5a~图 5c；图 3a）具有明显差异，表明基于调查样点单穴位观测的随机性，会影响黑土层厚度空间分布预测。
|
图 5 基于各观测穴位与调查样点的黑土层厚度空间预测分布差异 Fig. 5 Differences in the spatially predicted distribution of black soil layer thickness based on each observation site and the survey sample sites |
以调查样点多穴位均值预测的黑土层厚度空间分布（图 3a）作为对照底图，分析样点不同观测穴位A、B、C对空间分布预测的影响。对照结果表明，观测穴位A预测的区域黑土层平均厚度增加了0.64 cm，在松嫩平原四周以及三江平原增加最明显，在大兴安岭、小兴安岭以及东南山地也被严重高估，而在松嫩平原核心区被显著低估（图 5a，图 5d）。对观测穴位B而言（图 5b，图 5e），黑土层平均厚度降低了0.95 cm，在松嫩平原四周、小兴安岭被显著高估，在松嫩平原、大兴安岭和东南山地则被明显低估。基于观测穴位C预测研究区平均厚度值（图 5c）提高了0.14 cm，在松嫩平原、大兴安岭和小兴安岭被显著高估，在松嫩平原四周、三江平原和东南山地被显著低估（图 5f）。
总之，样点不同穴位调查与预测的黑土层厚度存在差异，黑土层厚度被高估或低估具有随机性。
3.3 样点观测穴位对黑土层厚度空间预测不确定性的影响对比分析样点不同观测穴位A、B、C观测值与平均值定量分析的黑土层厚度空间分布不确定性指数（图 6a~图 6c；图 3e）之间的差异，表明样点单穴位观测也会影响空间预测不确定性。以调查样点多穴位均值预测的黑土层厚度空间不确定性指数分布（图 3e）作为对照底图，穴位A预测的区域黑土层厚度不确定性指数平均值降低了0.18（图 6d），说明以穴位A预测的黑土层厚度空间分布总体上可解释性增加，主要体现在大兴安岭和小兴安岭调查样点较为稀疏区域，而在调查样点稠密的松嫩平原黑土层厚度空间分布解释性却显著降低。同样，基于穴位B预测研究区黑土层厚度不确定性指数平均值（图 6b）下降了0.22（图 6e），降低了松嫩平原黑土层厚度空间分布解释性，增加了大兴安岭、小兴安岭以及东南山地的黑土层厚度空间分布解释性。对于穴位C而言（图 6c，图 6f），黑土层厚度不确定性指数平均值下降了0.16，总体不确定性呈下降趋势，在山地和平原黑土层厚度不确定性指数分别具有增加和下降趋势。
|
图 6 基于各观测穴位与调查样点的黑土层厚度空间不确定性指数分布差异 Fig. 6 Difference in the distribution of spatial uncertainty indices of the thickness of the black soil layer based on each observation point and the survey sample sites |
总之，样点不同穴位调查与预测的黑土层厚度不确定性指数之间存在差异，不确定性增加和下降具有较大的随机性。
从空间预测不确定性评估的准确性水平看，当预测值在预测区间的占比曲线越贴近1︰1线，即预测区间α与预测占比PICP值相当，表明黑土层厚度空间预测不确定性的准确性水平越高。基于调查样点平均观测值的预测区间α介于0 ~ 0.5之间，略微低于1︰1线，说明空间不确定性的准确性水平被低估，即预测区间宽度相比于预测值略窄；在预测区间α介于0.5 ~ 1之间相对贴合1︰1线，空间不确定性的准确性水平较高。穴位A在预测区间α介于0 ~ 1之间离1︰1线稍远，空间不确定解释能力稍弱，其准确性水平被明显低估。而穴位B离1︰1线最远，表明其空间不确定性的准确性水平最低，即预测区间宽度太窄，不确定性能力被低估。穴位C在预测区间α介于0 ~ 1之间总体略低于1︰1线，空间不确定性的准确性水平也被低估。
总之，基于调查样点平均观测值的空间预测不确定性的准确性水平相对较高，而基于单穴位观测值预测的空间不确定性的准确性水平均被低估，总体表现为穴位C优于穴位A，进一步优于穴位B。
3.4 样点观测穴位对黑土层厚度空间预测精度的影响通过对比分析T1、T2、T3三种交叉验证方法、4种样本数据集下的随机森林空间预测精度差异，表明样点多穴位观测较任何单穴位观测的模型预测精度高（表 5），从空间预测角度证明了样点多穴位观测的必要性。如在T1交叉验证方法下，观测穴位A、B、C相比于样点平均RMSE值分别提高了30.11%、32.19%、30.80%，MAE值升高明显，R2降幅均较大，极大降低了模型预测精度。同样，在T2、T3交叉验证方法下均得出相似结果。显然，多穴位快速调查获取样点均值法，能够有效提高模型预测精度，得到更准确的黑土层厚度空间分布。
|
|
表 5 三种交叉验证方法获取4种样本随机森林模型精度值 Table 5 Three cross-validation methods to obtain four-sample random forest model accuracy values |
多穴位快速调查获取样点均值法，能代表微域范围内黑土层厚度整体水平，而单穴位观测则只能反映微域范围内某个点的实际水平，不能代表微域整体水平。因为栅格化的环境是样点所在栅格的平均值，均值与均值的配对更能反映实际存在的规律。所以，在模型预测精度上，多穴位快速调查获取样点均值法预测精度较高。
3.5 基于随机森林模型的黑土层厚度空间预测潜能稳定性基于样点均值数据集在不同交叉验证方法下的模型拟合精度值差异较小（表 5），说明模型相对稳定，即模型具备较好的空间预测潜能。以RMSE为例，样点平均从10折交叉验证（T1）到留一交叉验证（T2）再到5倍10折交叉验证（T3），其值从21.59 cm到21.58 cm再到21.57 cm，基本稳定。同理，观测穴位A、B、C在不同交叉验证方法下的预测精度具有较好稳定性。观察T1、T2、T3三种交叉验证方法在4种样本下的预测值与实测值之间的散点图（图 7），点散射均未明显出现；同一样本下不同交叉验证方法的散点图格局也大致保持一致。上述特征均表明随机森林模型具备稳定的空间预测潜能。
|
注：基于样点平均观测值：（a）T1，（b）T2，（c）T3；基于穴位A观测值：（d）T1，（e）T2，（f）T3；基于穴位B观测值：（g）T1，（h）T2，（i）T3；基于穴位C观测值：（j）T1，（k）T2，（l）T3。 Note: Based on sample mean observations(a)T1, (b)T2, (c)T3;based on point A observations: (d)T1, (e)T2, (f)T3;based on point B observations: (g)T1, (h)T2, (i)T3;based on point C observations: (j)T1, (k)T2, (l)T3. 图 7 三种不同方法交叉验证散点图 Fig. 7 Cross-validation of scatterplots by three different methods |
（1）建立了调查样点黑土层厚度多穴位“浅层挖掘+深层土钻”快速勘察新方法，检验结果表明，与剖面法获取的黑土层厚度基本一致，精度达到94.94%，可以替代剖面调查方法。（2）基于样点的空间预测证明了多穴位观测的科学性及必要性。原因在于土壤样点微域空间依然存在变异性，单穴位观测随机性改变了模型变量重要性值，也影响了黑土层厚度空间分布预测；相比于调查样点多穴位平均值的模型预测，基于单穴位观测值的预测精度显著下降，不确定性评估的准确性也降低。（3）构建及优化的随机森林模型，具有稳定的黑土层厚度空间预测潜能。通过特征选择、模型性能评估、模型参数调优等步骤，确定了随机森林最优模型。随机森林最优模型与默认模型相比，将模型性能指标R2提高了5.22%，RMSE、MAE分别降低了4.60%、5.42%，黑土层厚度空间变异解释力R2达到60%。对于不同样本集在相同规模下的黑土层厚度空间预测，随机森林模型具有极强稳定性。（4）基于调查样点与空间预测统计的当前黑龙江省平均黑土层厚度分别为55.56cm和51.73 cm，预测耕地平均黑土层厚度为53.42 cm。其中，厚层黑土主要分布在东北三江平原和西南松嫩平原，中层黑土分布在大兴安岭、小兴安岭以及东南山地，薄层黑土分布在小兴安岭东南部地区以及东南山地西部地区。黑土层厚度空间分布格局，可为研究揭示黑土层变薄特征、强化黑土资源保护等提供科学依据。
| [1] |
Nachtergaele F. Soil taxonomy—a basic system of soil classification for making and interpreting soil surveys: Second edition, by Soil Survey Staff, 1999, USDA–NRCS, Agriculture Handbook number 436, Hardbound[J]. Geoderma, 2001, 99: 336-337. DOI:10.1016/S0016-7061(00)00097-5
(  3) |
| [2] |
Sadras V O, Calviño P A. Quantification of grain yield response to soil depth in soybean, maize, sunflower, and wheat[J]. Agronomy Journal, 2001, 93(3): 577-583. DOI:10.2134/agronj2001.933577x
( 1) |
| [3] |
Zhang Z Y. The thickness changes of Ah horizon after the phaeozems cultivated (In Chinese)[J]. Journal of Heilongjiang August First Land Reclamation University, 2010, 22(5): 1-3. DOI:10.3969/j.issn.1002-2090.2010.05.001 [张之一. 黑土开垦后黑土层厚度的变化[J]. 黑龙江八一农垦大学学报, 2010, 22(5): 1-3.]
( 2) |
| [4] |
Duan X W, Xie Y, Ou T H, et al. Effects of soil erosion on long-term soil productivity in the black soil region of northeastern China[J]. Catena, 2011, 87(2): 268-275. DOI:10.1016/j.catena.2011.06.012
( 0) |
| [5] |
Gu Z J, Xie Y, Gao Y, et al. Quantitative assessment of soil productivity and predicted impacts of water erosion in the black soil region of northeastern China[J]. Science of the Total Environment, 2018, 637/638: 706-716. DOI:10.1016/j.scitotenv.2018.05.061
( 0) |
| [6] |
Pieters A, Baruch Z. Soil depth and fertility effects on biomass and nutrient allocation in jaraguagrass[J]. Journal of Range Management, 1997, 50(3): 268-273. DOI:10.2307/4003728
( 1) |
| [7] |
Gordeef T P. Description of the soils and rocks amid which has been found a mammoth - task. The society for the study of Manchuria, Harbin Bul. 6, 1926.
( 2) |
| [8] |
Pendleton R L, Hou K. Reconnaissance soil survey of the Harbin region[J]. Soil Bulletin, 1935, 11: 42-106.
( 1) |
| [9] |
Heilongjiang Province Agriculture Department. Heilongjiang soil[M]. Printed by the Department of Agriculture of Heilongjiang Province, 1962: 15—20.[黑龙江省农业厅. 黑龙江土壤[M]. 黑龙江省农业厅印, 1962: 15—20.]
( 1) |
| [10] |
Heilongjiang Provincial Land Management Bureau, Heilongjiang Provincial Soil Census Office. Heilongjiang soil (In Chinese). Beijing: Agriculture Press, 1992. [黑龙江省土地管理局, 黑龙江省土壤普查办公室. 黑龙江土壤[M]. 北京: 农业出版社, 1992.]
( 1) |
| [11] |
National Soil Census Office. Soil of China (In Chinese). Beijing: China Agriculture Press, 1998. [全国土壤普查办公室. 中国土壤[M]. 北京: 中国农业出版社, 1998.]
( 2) |
| [12] |
National Soil Census Office. Soil species of China (In Chinese). Beijing: China Agriculture Press, 1995. [全国土壤普查办公室. 中国土种志[M]. 北京: 中国农业出版社, 1995.]
( 2) |
| [13] |
Lu J L. Phaeozem degradation and sustainable agriculture in China (In Chinese)[J]. Journal of Soil and Water Conservation, 2001, 15(2): 53—55, 67. DOI:10.3321/j.issn:1009-2242.2001.02.014 [陆继龙. 我国黑土的退化问题及可持续农业[J]. 水土保持学报, 2001, 15(2): 53—55, 67.]
( 2) |
| [14] |
Shen B, Fan J R, Pan Q B, et al. General situation of pilot project of comprehensive prevention and control of soil and water loss of chernozem region in northeast of China (In Chinese)[J]. Soil and Water Conservation in China, 2003(11): 7-8. DOI:10.3969/j.issn.1000-0941.2003.11.004 [沈波, 范建荣, 潘庆宾, 等. 东北黑土区水土流失综合防治试点工程项目概况[J]. 中国水土保持, 2003(11): 7-8.]
( 1) |
| [15] |
Liu K, Wei M H, Dai H M, et al. Spatiotemporal variation of black soil layer thickness in black soil region of northeast China (In Chinese)[J]. Geology and Resources, 2022, 31(3): 434—442, 394. [刘凯, 魏明辉, 戴慧敏, 等. 东北黑土区黑土层厚度的时空变化[J]. 地质与资源, 2022, 31(3): 434—442, 394.]
( 2) |
| [16] |
Zhang S, Liu G, Chen S L, et al. Assessing soil thickness in a black soil watershed in northeast China using random forest and field observations[J]. International Soil and Water Conservation Research, 2021, 9(1): 49-57. DOI:10.1016/j.iswcr.2020.09.004
( 3) |
| [17] |
Zhang F P, Gao Z, Ma Q Q, et al. Construction of minimum data set for soil quality assessment in the Dunhuang oasis (In Chinese)[J]. Chinese Journal of Soil Science, 2017, 48(5): 1047-1054. [张福平, 高张, 马倩倩, 等. 面向敦煌市绿洲土壤质量评价的最小数据集构建研究[J]. 土壤通报, 2017, 48(5): 1047-1054.]
( 1) |
| [18] |
Zhang S W, Song X X, Wang Y, et al. Estimation of effective soil layer thickness of newly increased cultivated land based on GPR envelope detection (In Chinese)[J]. Journal of Anhui Science and Technology University, 2022, 36(5): 54-64. [张世文, 宋孝心, 王阳, 等. 基于探地雷达包络检波的新增耕地有效土层厚度估算[J]. 安徽科技学院学报, 2022, 36(5): 54-64.]
( 1) |
| [19] |
Ziadat F M. Prediction of soil depth from digital terrain data by integrating statistical and visual approaches[J]. Pedosphere, 2010, 20(3): 361-367. DOI:10.1016/S1002-0160(10)60025-2
( 1) |
| [20] |
Gessler P E, Moore I D, McKenzie N J, et al. Soil-landscape modelling and spatial prediction of soil attributes[J]. International Journal of Geographical Information Systems, 1995, 9(4): 421-432. DOI:10.1080/02693799508902047
( 1) |
| [21] |
Kuriakose S L, Devkota S, Rossiter D G, et al. Prediction of soil depth using environmental variables in an anthropogenic landscape, a case study in the Western Ghats of Kerala, India[J]. Catena, 2009, 79(1): 27-38. DOI:10.1016/j.catena.2009.05.005
( 1) |
| [22] |
Scarpone C, Schmidt M G, Bulmer C E, et al. Modelling soil thickness in the critical zone for Southern British Columbia[J]. Geoderma, 2016, 282: 59-69. DOI:10.1016/j.geoderma.2016.07.012
( 1) |
| [23] |
Lacoste M, Mulder V L, Richer-de-Forges A C, et al. Evaluating large-extent spatial modeling approaches: A case study for soil depth for France[J]. Geoderma Regional, 2016, 7(2): 137-152. DOI:10.1016/j.geodrs.2016.02.006
( 1) |
| [24] |
Lu Y Y, Liu F, Zhao Y G, et al. An integrated method of selecting environmental covariates for predictive soil depth mapping[J]. Journal of Integrative Agriculture, 2019, 18(2): 301-315. DOI:10.1016/S2095-3119(18)61936-7
( 1) |
| [25] |
Liu F, Yang F, Zhao Y G, et al. Predicting soil depth in a large and complex area using machine learning and environmental correlations[J]. Journal of Integrative Agriculture, 2022, 21(8): 2422-2434. DOI:10.1016/S2095-3119(21)63692-4
( 4) |
| [26] |
Meinshausen N, Ridgeway G. Quantile regression forests[J]. Journal of Machine Learning Research, 2006, 7(6): 983-999.
( 2) |
| [27] |
Zhang Z Y, Zhai R C, Cai D L. Introduction to soil series in Heilongjiang Province (In Chinese). Harbin: Harbin Cartographic Publishing House, 2006. [张之一, 翟瑞常, 蔡德利. 黑龙江土系概论[M]. 哈尔滨: 哈尔滨地图出版社, 2006.]
( 1) |
| [28] |
Gong Z T, Chen Z C, Shi X Z, et al. Chinese Soil Taxonomy: Theory, methodology and practice (In Chinese). Beijing: Science Press, 1999. [龚子同, 陈志诚, 史学正, 等. 中国土壤系统分类: 理论·方法·实践[M]. 北京: 科学出版社, 1999.]
( 2) |
| [29] |
Xie Y, Gao Y, Gu Z J, et al. Potential hazard assessment of slope farmland soil erosion in black soil region of northeastern China (In Chinese)[J]. Science of Soil and Water Conservation, 2020, 18(6): 105-114. [谢云, 高燕, 顾治家, 等. 东北黑土区坡耕地水土流失危险程度评价[J]. 中国水土保持科学, 2020, 18(6): 105-114.]
( 1) |
| [30] |
Takoutsing B, Heuvelink G B M. Comparing the prediction performance, uncertainty quantification and extrapolation potential of regression kriging and random forest while accounting for soil measurement errors[J]. Geoderma, 2022, 428: 116192. DOI:10.1016/j.geoderma.2022.116192
( 2) |
| [31] |
Chen Y, Ma L X, Yu D S, et al. Comparison of feature selection methods for mapping soil organic matter in subtropical restored forests[J]. Ecological Indicators, 2022, 135: 108545. DOI:10.1016/j.ecolind.2022.108545
( 1) |
| [32] |
Breiman L. Random forests[J]. Machine Learning, 2001, 45(1): 5-32. DOI:10.1023/A:1010933404324
( 1) |
| [33] |
Kuhn M. Building predictive models in R using the caret package[J]. Journal of Statistical Software, 2008, 28(5): 1-26.
( 1) |
| [34] |
Armstrong J S, Collopy F. Error measures for generalizing about forecasting methods: Empirical comparisons[J]. International Journal of Forecasting, 1992, 8(1): 69-80. DOI:10.1016/0169-2070(92)90008-W
( 1) |
| [35] |
Tian Z Y, Liu F, Liang Y, et al. Mapping soil erodibility in southeast China at 250 m resolution: Using environmental variables and random forest regression with limited samples[J]. International Soil and Water Conservation Research, 2022, 10(1): 62-74. DOI:10.1016/j.iswcr.2021.06.005
( 1) |
| [36] |
Malone B P, McBratney A B, Minasny B. Empirical estimates of uncertainty for mapping continuous depth functions of soil attributes[J]. Geoderma, 2011, 160(3/4): 614-626.
( 1) |
| [37] |
Zhang G L, Li D C. Manual of soil description and sampling (In Chinese). Beijing: Science Press, 2016. [张甘霖, 李德成. 野外土壤描述与采样手册[M]. 北京: 科学出版社, 2016.]
( 1) |
| [38] |
Yang G Y, Zhao H M, Tong D Q, et al. Impacts of post-harvest open biomass burning and burning ban policy on severe haze in the Northeastern China[J]. Science of the Total Environment, 2020, 716: 136517. DOI:10.1016/j.scitotenv.2020.136517
( 2) |
| [39] |
Zhao H M, Yang G Y, Tong D Q, et al. Interannual and seasonal variability of greenhouse gases and aerosol emissions from biomass burning in Northeastern China constrained by satellite observations[J]. Remote Sensing, 2021, 13(5): 1005.
( 1) |
| [40] |
Geng Z C, Dai W. Agrology (In Chinese). Beijing: Science Press, 2011. [耿增超, 戴伟. 土壤学[M]. 北京: 科学出版社, 2011.]
( 2) |
| [41] |
Bei'an Soil Census Office. Bei'an soil[M]. Printed by Bei'an Soil Census Office, 1985.[北安市土壤普查办公室. 北安土壤[M]. 北安市土壤普查办公室印, 1985.]
( 1) |