播种前裸土光谱代理变量与结荚期冠层指数相辅相成，助力基于无人机影像的田块尺度大豆产量预测

《iScience》：Pre-season bare-soil spectral proxies complement pod-filling canopy indices for field-scale soybean-yield prediction from UAV imagery

　　

**要点**
- 季前裸土光谱代理变量与结荚期冠层指数形成互补
- 加入裸土代理变量后，外部测试预测决定系数从0.240提升至0.325
- 在全部10次重复中，引入完整裸土代理变量池均使平均绝对误差下降
- 在未见的坡度单元上性能下降，支持了外部验证的必要性

**摘要**
基于冠层影像的田间尺度大豆产量模型可能遗漏土壤地表变异。我们评估了一种结合季前裸土和结荚期多光谱影像的双架次无人机框架，在中国东北四个坡度单元的101个地点开展研究。在嵌套特征选择和重复外部测试评估下，比较了五种随机森林情景。将裸土代理变量池加入结荚期植被指数后，外部测试平均预测决定系数从0.240提升至0.325，平均绝对误差从573.06降至518.76千克/公顷；该误差在全部10次重复中均有下降。模型归因分析显示冠层和裸土预测变量均有贡献。预测值对低产量高估、对高产量低估。在该组合情景下，当逐一将每个坡度单元排除出模型开发时，合并预测决定系数降至0.161。因此，季前裸土影像在结荚期采集之外提供了适度的互补信息，但该框架仍属于田间尺度概念验证，需要多年份和多地点的验证。

**图解摘要**

**引言**
大豆（*Glycine max* (L.) Merrill）是全球主要的植物油和植物蛋白来源。中国东北是国家级大豆主产区。及时且空间明确的 大豆产量预测对于作物评估和农业决策具有重要意义。然而，作物产量估算本质上十分复杂。传统的田间产量评估劳动密集型，且提供的空间细节有限。基于过程的模型，如WOFOST、DSSAT和ALMANAC，可以从机制上表征作物发育，但需要详细的作物、土壤、气象和管理输入数据。遥感可以通过在收获前提供空间分布的观测数据来补充这些方法。

配备多光谱传感器的无人机可在田间尺度提供灵活的高空间分辨率作物冠层观测。由这些观测衍生出的植被指数已被用于刻画生长期冠层变异并预测大豆产量。将冠层光谱信息与结构、热红外或纹理特征相结合，可以为大豆产量预测提供互补信息，尽管多传感器或多时相方案也需要额外的数据获取和处理工作。然而，冠层衍生预测变量主要描述影像采集时刻所表达的作物状况，可能无法充分表征土壤和根系环境中持续存在的空间变异。大豆发育和产量也会随土壤质地、土壤水分和地形位置的变化而不同。一旦冠层郁闭，从上方直接观测土壤地表信号变得困难；因此，季前裸土采集为作物出苗前捕捉互补的土壤地表光谱变异提供了机会。可见光和近红外土壤反射率受多种土壤和地表条件的共同影响，因此由裸土影像衍生的指数应被解释为与土壤相关的光谱代理变量，而非单个土壤属性的直接测量值。这些代理变量是否能提供超越及时冠层采集的预测信息，在田间尺度大豆产量建模中尚未得到充分检验。

随机森林（RF）是由随机化回归树组成的集成模型，能够处理非线性关系和预测变量间的交互作用。由于此类集成模型的预测无法直接归因于单个输入变量，Shapley加性解释（SHAP）可用于量化每个预测变量对拟合预测的加性贡献，并汇总跨预测的贡献。然而，模型解释并不能纠正有偏的验证设计。如果预测变量选择或超参数调优利用了后续用于评估的数据，所得误差估计可能过于乐观；因此，这些模型开发步骤应嵌套在训练数据内。在具有空间结构的田间数据中，随机划分也可能将相邻或环境相似的观测同时放入训练集和测试集，而预测未见过的空间单元则需要与该任务相匹配的空间分离验证设计。因此，预测性能应使用排除在模型开发之外的样本进行评估，空间迁移能力应在解释模型之前单独评估。

本研究评估了一种结合季前裸土多光谱影像和单次结荚期冠层采集的双架次无人机田间尺度大豆产量预测框架，在中国东北两个站点开展（图1）。比较了五种预测变量情景，以检验在将预测变量选择和超参数调优限定于训练数据的嵌套验证设计下，与土壤相关的光谱代理变量是否能在结荚期植被指数之外提供增量预测信息。通过留出法坡度单元（LOSO）验证检验了向完全未见坡度单元的迁移能力，并使用SHAP描述了所选冠层和裸土预测变量对排除在模型开发之外的样本预测的贡献方式。由于数据集仅代表一个生长季、两个站点和101个观测值，本研究定位为田间尺度概念验证，而非广泛验证的区域预测框架。

**图1** 地理位置与研究田间

(A) 中国东北地区位置；(B) 中国东北大豆分布，数据源自Liu等；(C) 海伦和北安的位置；(D) 蓝点代表样本点分布；黑色调影像代表裸土阶段，绿色调影像代表结荚期阶段。面板(A)–(C)中的比例尺表示公里距离。

**结果**

**嵌套预测变量选择与土壤相关光谱代理变量的外部测试标定**

在101个观测值中，大豆产量平均为3438.23 ± 842.29千克/公顷，范围为1305.82至5184.32千克/公顷。土壤有机碳（SOC）、砂粒、粉粒和黏粒含量的平均值分别为25.50 ± 9.72克/千克、12.24% ± 2.76%、53.48% ± 6.28%和34.28% ± 5.83%（表S1）。

嵌套特征选择在各10次外部重复中产生了不同的预测变量集（表S3）。在仅植被光谱指数（VSI）情景中，结荚期LCI和结荚期NDRE在10次重复中的各7次被选中。由于每个外部重复中选出的VSI子集在组合情景中被保留，这两个指数在VSI+SOC、VSI+texture和VSI+土壤光谱指数（SSI）情景中具有相同的选中频率。在VSI+SSI情景中，裸土NCMI在9次重复中被选中，裸土RI在3次重复中被选中，而裸土IV3和裸土RDVI各在2次重复中被选中。裸土DVI、裸土LCI、裸土TCARI/OSAVI和裸土CIrededge各被选中1次。没有预测变量在全部10次重复中均被保留，表明所选预测变量集随外部训练数据组成而变化。

外部测试标定在四个实测土壤目标之间有所不同（表S4）。SOC相关代理变量模型实现了0.543 ± 0.060的平均预测决定系数和0.592 ± 0.085的平均Pearson r2，均方根误差（RMSE）为6.40 ± 0.67克/千克，相对均方根误差（rRMSE）为25.18% ± 1.90%。相比之下，质地相关代理变量模型的平均预测决定系数分别为砂粒?0.031 ± 0.193、粉粒?0.128 ± 0.973和黏粒?0.072 ± 1.094。相应的Pearson r2值分别为0.086 ± 0.099、0.424 ± 0.207和0.507 ± 0.206，RMSE值分别为2.74 ± 0.36、5.34 ± 2.35和4.53 ± 2.71个百分点。仅SOC相关代理变量模型实现了正的平均预测决定系数。因此，裸土指数被作为与土壤相关的光谱代理变量处理，而非SOC、砂粒、粉粒或黏粒的直接估计值。

**裸土光谱代理变量为外部测试产量预测提供了适度的增量价值**

所有五种预测变量情景在10次重复划分中表现出差异性的外部测试性能（图2；表1）。仅VSI情景实现了0.240 ± 0.251的平均预测决定系数，均方根误差为695.92 ± 108.75千克/公顷，平均绝对误差（MAE）为573.06 ± 81.12千克/公顷，rRMSE为20.21% ± 3.21%。仅SSI情景的平均预测决定系数最低（0.071 ± 0.255），平均均方根误差最高（771.44 ± 98.18千克/公顷）。在跨重复平均值层面，三个组合情景的预测决定系数均高于仅VSI情景，且均方根误差、平均绝对误差和rRMSE均低于仅VSI情景。VSI+SOC、VSI+texture和VSI+SSI的平均预测决定系数分别为0.302 ± 0.226、0.321 ± 0.204和0.325 ± 0.199，对应的均方根误差分别为665.10 ± 91.99、658.50 ± 99.38和656.15 ± 90.30千克/公顷。其平均绝对误差分别为532.64 ± 56.82、524.27 ± 70.63和518.76 ± 65.62千克/公顷。因此，VSI+texture和VSI+SSI表现出非常接近的平均外部测试性能。VSI+SOC和VSI+texture标签标识了用于构建相应裸土代理变量候选池的实测土壤目标；其产量预测性能并不代表SOC或土壤质地估计的精度。

**图2** 五种预测变量情景的外部测试性能

10次重复外部划分中(A)预测决定系数、(B)均方根误差和(C)平均绝对误差的分布。点代表各外部重复，箱体表示中位数和四分位距，须延伸至四分位距1.5倍范围内的极端值，菱形表示跨重复平均值。面板(A)中的水平虚线表示预测决定系数=0。五种情景使用相同的外部划分进行评估。均方根误差和平均绝对误差以千克/公顷表示。VSI表示结荚期植被光谱指数；SSI表示裸土光谱指数预测变量组，解释为与土壤相关的光谱代理变量；SOC表示土壤有机碳。

**表1** 重复嵌套验证下五种预测变量情景的外部测试性能

数值为10次外部测试评估的平均值 ± 标准差。所有特征选择和超参数调优仅使用相应的外部训练数据完成；外部测试样本仅用于最终模型评估。预测决定系数计算为1 ? SSE/SST，Pearson r2计算为观测产量与预测产量之间Pearson相关系数的平方。均方根误差和平均绝对误差以千克/公顷表示，rRMSE以百分比表示。VSI表示从结荚期影像衍生的植被光谱指数。在情景标签中，SSI表示裸土光谱指数预测变量组；该组中的指数被解释为与土壤相关的光谱代理变量。SOC，土壤有机碳；SD，标准差；SSE，误差平方和；SST，总平方和；RMSE，均方根误差；MAE，平均绝对误差；rRMSE，相对均方根误差。

| 情景 | 预测决定系数 | Pearson r2 | RMSE/千克·公顷?1 | MAE/千克·公顷?1 | rRMSE/% |
|------|------------|-----------|----------------|----------------|---------|
| 仅VSI | 0.24±0.25 | 0.35±0.17 | 695.92±108.75 | 573.06±81.12 | 20.21±3.21 |
| 仅SSI | 0.07±0.25 | 0.19±0.14 | 771.44±98.18 | 597.55±72.22 | 22.46±3.50 |
| VSI+SOC | 0.30±0.23 | 0.40±0.20 | 665.10±91.99 | 532.64±56.82 | 19.30±2.62 |
| VSI+texture | 0.32±0.20 | 0.39±0.17 | 658.50±99.38 | 524.27±70.63 | 19.14±3.06 |
| VSI+SSI | 0.33±0.20 | 0.38±0.17 | 656.15±90.30 | 518.76±65.62 | 19.07±2.87 |

内部折外与外部测试误差估计在五种情景中总体相似（表S5）。对于VSI+SSI，内部折外均方根误差和平均绝对误差分别为648.14 ± 16.18和516.68 ± 14.18千克/公顷，与相应的外部测试值656.15 ± 90.30和518.76 ± 65.62千克/公顷密切吻合。内部折外预测决定系数为0.403 ± 0.033，而外部测试评估为0.325 ± 0.199。由于外部测试样本未参与特征选择、超参数调优或模型拟合，所有预测性能的主要比较均基于外部测试结果。

仅VSI和VSI+SSI在10次相同外部划分中的重复级预测决定系数和均方根误差值，以及各重复内的差值，列于表S6。基于这些匹配的重复级差值的配对比较表明，加入完整裸土代理变量池相对于仅VSI提供了适度的增量价值（图3；表S7）。对于VSI+SSI，预测决定系数平均增加0.0857，在10次外部重复中的8次得到改善（95%自助法置信区间：0.0185至0.1488；p = 0.037）。均方根误差平均减少39.77千克/公顷，在8次重复中改善（95%自助法置信区间：?66.90至?8.84千克/公顷；p = 0.041）。平均绝对误差平均减少54.29千克/公顷，在全部10次重复中均改善（95%自助法置信区间：?73.82至?34.66千克/公顷；p = 0.002）。平均rRMSE差值为?1.134个百分点，在8次重复中改善（95%自助法置信区间：?1.918至?0.238个百分点；p = 0.039）。因此，平均绝对误差的下降是最一致的改善，而预测决定系数、均方根误差和rRMSE在大多数但非所有外部重复中得到改善。VSI+texture在10次外部重复中的9次也在预测决定系数、均方根误差、平均绝对误差和rRMSE方面表现出有利差异，而VSI+SOC的增益则不那么一致（表S7）。

**图3** VSI+SSI与仅VSI之间的配对外部测试性能差异

10次相同外部划分中(A)预测决定系数、(B)均方根误差和(C)平均绝对误差的差异。差异计算为VSI+SSI减去仅VSI。正值表示预测决定系数的有利差异，负值表示均方根误差和平均绝对误差的有利差异。实心点表示有利差异，空心点表示不利差异，垂直线从零延伸至观测到的配对差异值。

虚线水平线表示各情景之间没有差异。RMSE和MAE差异以kg ha?1表示。VSI代表结荚期植被光谱指数；SSI代表从裸土影像中提取的土壤相关光谱代理变量。

**预测向均值偏移以及在未见坡面上的性能下降**

在10次重复试验中累积的200条外测试预测记录中，预测偏差在由相应外训练数据定义的产量层次中表现出一致的方向性模式（表S8）。对于VSI+SSI，低产量层次的平均偏差为+460.76 kg ha?1，中产量层次为+164.02 kg ha?1，高产量层次为?571.88 kg ha?1。仅VSI也表现出相同的方向性模式，相应偏差分别为+527.86、+143.69和?630.73 kg ha?1。因此，低产量通常被高估，而高产量则被低估。与仅VSI相比，VSI+SSI降低了低产量和高产量层次的偏差幅度，尽管中产量层次的正偏差略有增加，且预测向均值回归的总体趋势仍然存在。

留一坡面法（LOSO）验证为101个观测值中的每一个产生了唯一的外测试预测，其汇总性能低于重复随机划分（图S1A和S1B）。对于VSI+SSI，预测R2从重复随机划分下的均值0.325降至LOSO下的汇总值0.161，而RMSE从656.15增至767.64 kg ha?1。在LOSO下，三个组合情景表现出非常接近的汇总性能：VSI+SOC的预测R2为0.158，VSI+texture为0.167，VSI+SSI为0.161。VSI+texture的预测R2数值最高且RMSE最低（765.10 kg ha?1），但其性能与其他组合情景非常接近。在两种验证方案下，没有任何单一组合情景在数值上表现最佳。

VSI+SSI在四个完全留出的坡面上的性能差异显著（图S1C和S1D）。预测R2值分别为BAA的?0.272、BAB的?0.019、HLA的0.040和HLB的0.238，相应的RMSE值分别为694.26、1016.86、536.40和495.31 kg ha?1。因此，BAB的RMSE最高，而HLB的预测R2最高且RMSE最低。在LOSO下，所选预测变量的分布差异也更大（图S2）。在重复随机划分下VSI+SSI的34个"特征×重复"记录中，中位绝对标准化均差（|SMD|）为0.058，中位Kolmogorov-Smirnov（KS）统计量为0.140。在它们被选中的各个LOSO折中，裸土NCMI和裸土TCARI/OSAVI的|SMD|值分别为4.071和3.599，KS统计量分别为1.000和0.973。这些较大的分布差异出现在预测性能较低且在留出的坡面间存在差异的同一LOSO分析中。综合来看，这些结果表明，采样坡面内的性能并不能均匀地推广到未见坡面。

**SHAP归因显示冠层和裸土预测变量均对外测试预测有贡献**

选择加权的全局SHAP重要性在各预测变量情景间有所不同，但组合情景表现出一种共同模式：结荚期冠层指数仍然显著，而裸土预测变量提供了额外的贡献（图S3）。在仅VSI中，结荚期LCI和结荚期NDRE的选择加权平均绝对SHAP值分别为240.5和225.5 kg ha?1。在仅SSI中，裸土NCMI排名第一，为172.4 kg ha?1，而排名第二的预测变量裸土RDVI和裸土RI的值显著较低，分别为56.0和54.0 kg ha?1。在每个组合情景中，结荚期LCI和结荚期NDRE仍然是排名最高的两个预测变量。排名最高的裸土预测变量分别为VSI+SOC中的裸土GNDVI（62.2 kg ha?1）、VSI+texture中的裸土NCMI（104.0 kg ha?1）和VSI+SSI中的裸土NCMI（126.4 kg ha?1）。

在VSI+SSI中，结荚期LCI、结荚期NDRE和裸土NCMI的选择加权平均绝对SHAP值最高，分别为195.4、167.5和126.4 kg ha?1（图4A）。这些预测变量在10次外重复中的7次、7次和9次中被选中，而其余每个裸土预测变量最多仅在3次重复中被选中，且选择加权平均绝对SHAP值为26.0 kg ha?1或以下。领先预测变量之间的自助法置信区间存在重叠，表明尽管相同的三个预测变量反复显著，但其确切排序存在不确定性。有符号SHAP分布显示了领先预测变量中贡献方向的趋势（图4B）。结荚期LCI和结荚期NDRE的高值更倾向于出现在正SHAP侧，而其低值则强烈延伸至负贡献区域。裸土NCMI的高值更常与负SHAP贡献相关联，而其低值倾向于接近零或位于正侧。由于特征值在每个预测变量内进行了单独的标准化，颜色模式应在预测变量内部而非跨预测变量进行解读。相应的观测值范围和汇总统计量见表S9。

**图4 VSI+SSI情景外测试产量预测的全局SHAP归因**

显示完整标题 | 图形查看器

（A）在10次外重复中保留的预测变量的选择加权全局重要性。点表示选择加权平均绝对SHAP值，水平误差条表示相应的95%自助法置信区间，括号中的值表示每个预测变量被选中的外重复次数。未选中某预测变量的重复对选择加权均值贡献为零。

（B）在外测试预测中实际选中预测变量的有符号SHAP贡献。每个点代表一个外测试预测中的一个选中预测变量，颜色表示该特征内标准化后的预测变量值。正的SHAP值使模型预测相对于其基线增加，而负值则使其减少。SHAP值以kg ha?1表示，解释的是模型预测而非对观测产量的因果效应。

VSI，植被光谱指数；SSI，土壤相关光谱代理变量；SHAPs，SHapley加性解释。在（A）中，蓝色和柔和的陶土色符号分别表示结荚期VSI和裸土SSI预测变量。在（B）中，蓝色至赭石色渐变表示每个预测变量内从低到高的标准化值。

局部SHAP分解说明了冠层和裸土预测变量如何共同形成个体的外测试预测（图S4）。此处"局部"指的是单个预测的加性分解，而非地理或空间尺度。对于北安低产量案例BAB34，最大的负贡献与结荚期LCI和结荚期NDRE相关，每个使预测减少约590 kg ha?1，而裸土NCMI和裸土TCARI/OSAVI分别贡献了约?114和?17 kg ha?1。相比之下，对于北安中产量案例BAA3，裸土NCMI、结荚期LCI、结荚期NDRE和裸土RI分别使预测增加了约248、241、231和150 kg ha?1。在六个案例中，冠层和裸土预测变量均以正负方向做出贡献，尽管某些个体预测主要受单一方向的贡献支配。这六个案例说明了拟合模型如何构建个体预测；它们未被用于推断特定产量层次的普遍关系，且SHAP贡献描述的是模型行为而非对观测产量的因果效应。

**讨论**

**裸土光谱代理变量在结荚期冠层指数之上提供了适度的增量价值**

主要发现是结荚期VSI为大豆产量预测提供了更强的独立信号，而将裸土光谱代理变量添加到仅VSI基线中提供了适度的额外预测信息。与仅VSI相比，VSI+SSI将预测R2提高了0.0857，并将RMSE、MAE和rRMSE分别降低了39.77 kg ha?1、54.29 kg ha?1和1.134个百分点（图3；表S7）。在全部10次外重复中，MAE均有改善，而预测R2、RMSE和rRMSE在8次重复中有改善。对于VSI+SSI与仅VSI的比较，预测R2和RMSE的配对增益得到了精确配对符号翻转置换检验的支持（p = 0.037和p = 0.041，分别对应），且相应的95%自助法置信区间排除了零。VSI+texture在9次重复中也改善了所有四个指标，而仅SSI的预测能力低于仅VSI（图2；表1）。综合来看，这些结果表明裸土光谱代理变量补充了结荚期VSI已经捕获的产量信息，而非替代冠层信号。

表2提供了紧凑的跨情景比较。三个组合情景的平均外测试RMSE和rRMSE均低于仅VSI，而仅SSI的误差最高。VSI+texture和VSI+SSI的均值误差水平非常相似。预测变量选择频率在组合情景内有所不同，但结荚期LCI在所有三个组合情景中始终是排名最高的选择加权SHAP特征。综合来看，这些模式表明裸土预测变量池提供了增量信息，而并未替代结荚期冠层信号的主导贡献。

| 情景 | 平均外测试RMSE (kg ha?1) | 平均外测试rRMSE (%) | 预测变量选择频率范围 (n/10) | 排名最高的SHAP特征 |
|---|---|---|---|---|
| 仅VSI | 695.92 | 20.21 | 7–7 | 结荚期LCI |
| 仅SSI | 771.44 | 22.46 | 1–9 | 裸土NCMI |
| VSI+SOC | 665.10 | 19.30 | 1–7 | 结荚期LCI |
| VSI+texture | 658.50 | 19.14 | 1–7 | 结荚期LCI |
| VSI+SSI | 656.15 | 19.07 | 1–9 | 结荚期LCI |

**表2 外测试误差、预测变量选择频率和SHAP归因的紧凑跨情景汇总**

RMSE和rRMSE为10次外测试评估的跨重复均值；相应的标准差见表1。预测变量选择频率范围表示在给定情景中至少被选中一次的预测变量被保留的重复次数的最小值和最大值。排名最高的SHAP特征是各情景内选择加权平均绝对SHAP值最高的预测变量，如图S3所述。VSI代表从结荚期影像中提取的植被光谱指数；SSI代表裸土光谱指数预测变量组，解释为土壤相关光谱代理变量；SOC，土壤有机碳；RMSE，均方根误差；rRMSE，相对均方根误差；SHAPs，SHapley加性解释。

在新标签页中打开表格

VSI基线的相对强度与冠层观测时机的选择一致。本研究的冠层航飞于2024年8月下旬的田间确定的结荚期进行。在Fu等人对黑龙江大豆进行的县尺度分析中（跨越4月至9月），8月被确定为最佳收获前预测月份，8月在整体误差指标方面表现最佳，且MAPE显著低于7月。最近从黑龙江的田块尺度无人机观测中获得了更直接的物候证据：在R3、R6、R7和R8阶段中，灌浆满阶段（R6）的单阶段预测精度最高。Johnson同样发现仲夏NDVI与大豆产量呈正相关，且季节内遥感变量在收获前仍保留了有用的预测信息。这些研究支持晚季生殖期冠层观测作为有信息量的产量预测基线，但日历月份和物候阶段在不同品种、地点和年份间不应被视为可以互换。

然而，即使是时机合适的冠层影像，也无需完全代表根际环境的持续性空间变异。冠层指数描述的是观测日期所表达的作物状况，而质地、水分和地形位置等方面的土壤相关梯度可在全季节影响植物发育。Falco等人通过将无人机和地球物理测量相结合，将大豆发育和产量的空间差异与土壤黏土含量、土壤水分、灌溉和地形控制联系起来。Franz等人同样发现，精细尺度的产量模式反映了植被状况、土壤含水量、地形和地块位置的综合作用。Pejak等人的田间大豆研究也将Sentinel-2观测和土壤参数纳入同一预测框架。在此背景下，从仅VSI到组合情景的改善与裸土影像提供了8月下旬冠层指数未能完全表征的空间背景信息是一致的。然而，由于输入是光谱代理变量，这些额外信息不应被解释为某一特定实验室测定的土壤性质的直接影响。

增益幅度适中，且在重复外划分中总体上是有利的。到8月下旬，冠层状况已经整合了部分全季作物对土壤和天气的响应，因此VSI与裸土信息之间的部分冗余是可以预期的。以往的基于无人机的豆研究同样表明，结合非冗余的光谱、结构、热或纹理信息可以改善产量预测，但效益取决于附加信息源的信息含量和验证方案。在本研究中，VSI+texture和VSI+SSI的外测试性能非常接近，且在重复随机划分和LOSO验证两种方案下，没有任何组合情景在数值上表现最佳。

因此，现有证据支持在两次飞行框架内利用裸土光谱代理变量所提供的信息——即在结荚期VSI基线之上增加有用信息——而不是声称某一特定代理变量池具有普遍优越性。

**冠层指数与土壤相关光谱代理变量代表了预测信号的不同组成部分**

SHAP结果显示，两个影像获取时期为预测信号贡献了不同的组成部分。

在VSI+SSI中，结荚期LCI和结荚期NDRE的选择加权平均绝对SHAP值最高，分别为195.4和167.5 kg ha?1，且各自在10次外循环中有7次被保留（图4A）。这两个指数较高的值倾向于更频繁地出现在SHAP正值一侧，而较低的值则强烈向负贡献方向延伸。NDRE对比了近红外和红边反射率，而本研究中使用的LCI公式结合了近红外、红边和红光反射率。这些指数可能受到叶片和冠层叶绿素、氮素状态、绿叶面积及冠层结构变化的影响。以往研究已证明可见光、红边和近红外反射率对叶绿素评估具有价值，而红边指数已被证明可用于作物冠层叶绿素和氮素估算。近期对大豆的无人机观测同样将红边、近红外及与绿色度相关的变量识别为叶片叶绿素含量的重要预测因子。,，因此，SHAP模式与以下观察结果一致：即绿色度和叶绿素相关的冠层信号较强的观测点被拟合模型赋予了更高的产量预测值。这一解释涉及在结荚期冠层中表达的作物状态，并不意味着LCI或NDRE直接导致了产量差异。

裸土NCMI代表了另一个不同但界定清晰的预测信号组成部分。它在10次VSI+SSI外循环中的9次被保留，且具有第三高的选择加权平均绝对SHAP值，为126.4 kg ha?1。NCMI在数学上被定义为近红外反射率与红绿反射率之和之间的归一化对比。它最初是为作物监测而开发的，因为其在水稻物候发育和管理处理过程中对这三个光谱区域变化的敏感性。然而，当从季前裸土影像中计算时，其生物物理解释发生了变化。在本研究中，裸土NCMI应被解释为一种综合性的土壤表面光谱对比，它概括了裸露土壤反射率的空间差异，而非植被管理指标。其频繁被选中且具有相对较大的SHAP重要性表明，这种光谱对比被拟合模型反复使用，并在重复的外循环中包含了与产量预测相关的信息。较高的裸土NCMI值倾向于与负SHAP贡献相关联，表明拟合模型将这种光谱构型与较低的预测产量联系在一起。这是一个有意义的预测关联，尽管其本身并不能确定导致该关联的物理土壤因子。

无法将裸土NCMI归因于某一个具体的土壤特性，并不意味着该指数缺乏物理意义。可见光和近红外土壤反射率共同受到有机质、矿物成分、铁氧化物、质地、颗粒排列、土壤含水量、表面粗糙度、作物残留物以及照明或观测几何因素的影响。这些影响在光谱上相互重叠，因此多光谱指数可以概括可重复的土壤表面状态，而无需作为某一实验室特性的选择性测量。仅土壤水分就能够在可见光和近红外反射率中产生显著的非线性变化，而田间土壤光谱学的综述强调，光谱与SOC、质地、矿物学和含水量之间的关系在很大程度上取决于校准域和共变的土壤条件。,，这一区别有助于解释为何SOC相关代理变量模型取得了正的外循环测试预测R2，而沙、粉沙和黏土相关模型则没有，同时质地相关代理变量池在添加到VSI中时仍能改善产量预测。估算实验室测量的土壤特性与检测与产量相关的复合土壤表面光谱模式是两种不同的预测任务。因此，裸土NCMI可能通过表征土壤和地表条件的混合状态来提供有用的产量信息，尽管其信号无法唯一归属于SOC、质地、含水量或其他单一测量特性。这就是它被描述为土壤相关光谱代理变量而非直接土壤特性估计值的原因。

然而，SHAP结果仍应被解释为模型归因，而非确定性的生理或土壤机制等级。选择加权重要性结合了预测变量被保留的频率及其被保留时的贡献大小。结荚期LCI、结荚期NDRE和裸土NCMI的反复突出表现比其精确数值排序更具信息量，特别是因为其自举置信区间相互重叠。LCI和NDRE还共享相关的光谱信息，且Shapley归因可以在具有依赖性或者相关的预测变量之间重新分配，具体取决于所假设的背景分布。更根本地，SHAP解释了预测变量如何对拟合模型输出做出贡献；它并不能独立建立对观测产量的因果效应。,因此，六个局部分解应被视为个体预测如何组装的示例，而非某一预测变量控制产量分层内所有观测值的证据（图S4）。所支持的解释是，结荚期冠层指数代表了季节后期表达的作物状态，而裸土NCMI及其他土壤相关光谱代理变量则提供了关于季前土壤表面空间变化的额外信息，这些信息被拟合模型所使用。

**预测向均值偏差以及对未见坡面单元的有限迁移能力**

产量分层残差显示，预测值被拉向观测产量分布的中部：低产量被高估，而高产量被低估。仅VSI和VSI+SSI均高估了低产量分层并低估了高产量分层，而中产量分层的偏差较小。低、中、高产量分层的累计记录数相近——分别为64、67和69——表明方向性偏差并非三个分层中外部测试记录数严重不平衡的人工伪迹（表S8）。更合理的解释是极端产量代表不足与回归森林平均化的组合。随机森林回归在终端节点和树之间汇总训练响应值；当连续响应分布尾部附近的观测值在训练数据中缺乏相近的类似项时，其预测值可能被拉向更密集表示的中心值。更普遍地，不平衡回归研究表明，连续目标分布中的稀有区域通常预测精度低于中心范围，即使总体误差指标看起来是可接受的。,，因此，观测到的偏差应被解释为代表完整产量范围的局限性，而非均匀加性校准误差的证据。

添加完整的裸土代理变量池减少但未能消除预测值向产量分布中心移动的趋势。与仅VSI相比，VSI+SSI将低产量分层的高估幅度从527.86降至460.76 kg ha?1，将高产量分层的低估幅度从630.73降至571.88 kg ha?1。然而，中产量偏差从143.69略微增加至164.02 kg ha?1。因此，额外的裸土信息帮助模型区分了产量分布两端的差异部分，但并未在所有产量水平上均匀改善校准。这一结果与所报告的VSI+SSI总体增量收益并不矛盾：配对RMSE和MAE概括了完整外部测试集上的平均误差，而分层分析则识别了残差误差仍然集中的位置。对于那些依赖区分异常低产和高产地点的应用而言，这种压缩很重要，因为它可能降低田间最具农学差异性部分之间表面上的对比。

在LOSO验证下的较低表现也应结合预测任务来理解。在重复的坡面分层随机划分中，每个坡面单元都同时出现在外部训练和外部测试划分中。因此，这些结果估计的是在已采样坡面单元集合内对新观测值的预测。LOSO施加了一个不同的且更具挑战性的任务：预测整个坡面单元，而其样本被排除在特征选择、调优和模型拟合之外。对于VSI+SSI，预测R2从重复随机划分下的0.325降至LOSO下的0.161，而RMSE从656.15增至767.64 kg ha?1（图S1A和S1B）。方法论研究反复表明，当训练数据和测试数据中存在空间或环境相似的观测值时，随机验证可能会夸大地理迁移能力，而空间分离验证在部署涉及新地点时提供了更相关的估计。,，作物产量研究也报告了同样的总体模式：基于无人机的大豆产量模型的空间验证效果远低于普通随机交叉验证，而近期的田间子区分析同样发现常规交叉验证比空间替代方法更为乐观。,相反，完全阻断设计评估了比已采样坡面单元内新观测值之间的内插更严格的迁移任务。因此，随机外部测试和LOSO结果是互补而非竞争的估计：前者描述研究内的预测，后者描述向未见坡面单元迁移的敏感性。

LOSO折叠之间的变异性进一步表明，迁移难度取决于特定的坡面单元。VSI+SSI的RMSE范围为从HLB的495.31 kg ha?1到BAB的1016.86 kg ha?1，预测R2范围为从BAA的-0.272到HLB的0.238（图S1C和S1D）。一些被选中的预测变量在LOSO下的训练-测试分布差异也比在重复随机划分下大得多；例如，裸土NCMI和裸土TCARI/OSAVI在被选中的折叠中标准化均值差的绝对值分别为4.071和3.599，KS统计量分别为1.000和0.973（图S2）。这些观察结果与向其余坡面弱代表的协变量空间部分进行预测的情况一致。空间预测研究强调，交叉验证误差在训练数据所代表的预测变量空间域内最为可靠，而当新地点环境不相似时可能不适用。本分析并未证明预测变量偏移是LOSO较差结果的唯一原因，因为只有四个坡面折叠可用，且各折叠中选择的特征集也不同。然而，它表明了合并后的关系并未在坡面单元之间均匀迁移。合并四个坡面满足了估计在采样研究域上共同预测关系的研究目标，而LOSO量化了将该关系迁移到完全未见坡面单元的额外难度。

**研究的局限性**

本研究是基于两个站点、四个坡面单元、单一生长季中101个观测值进行的田间尺度概念验证。重复的外循环量化了对不同训练-测试划分的敏感性，但并未增加独立年份、站点或空间单元的数量。被选中预测变量的变异性、跨循环表现、SHAP置信区间和LOSO结果因此反映了有限数据集所固有的真实不确定性。重复交叉验证可以描述这种不确定性，但无法替代额外的独立观测值，且小样本的性能估计可能仍存在相当大的误差。因此，所报告的结果应被解释为在采样研究域内的证据，而非东北中国大豆产量预测的性能基准。

此外，模型仅代表了两个无人机获取时期，并未包括产量的所有决定因素，如天气、土壤水分、品种、地形和管理。这些遗漏的因素可能对剩余预测误差以及低产高估和高产低估有所贡献。此外，裸土指数是光谱代理变量而非直接土壤测量。其对产量预测的贡献无法归属于某一单一实验室测量特性，尤其是因为质地相关代理变量模型的外循环测试校准较差。SHAP解释同样依赖于拟合模型、被选中的预测变量和观测数据。全局排名不应被视为确定性的机制等级，六个局部案例说明的是个体预测而非群体层面的效应或因果关系。

尽管如此，结果支持在结荚期VSI基线上添加季前裸土光谱代理变量的研究内价值。无泄漏嵌套工作流程旨在通过将土壤相关候选池构建、面向产量的特征选择和超参数调优限制在对应的训练数据内，来降低乐观偏差的风险。

外部测试集平均预测R2仅为0.325，且在使用留一观测单位（LOSO）验证时进一步下降至0.161，这进一步凸显了在将该框架考虑用于实际业务化应用之前，亟需开展独立的多年份、多站点重复验证。此类验证还应纳入额外的空间单元、不同品种和管理条件，以及更广泛的天气与土壤水分环境。预测因子的选择应嵌套在模型开发流程内部，性能评估应使用独立的空间和时间数据集，而不只是对相同观测样本进行新的随机划分。预测结果也应限定在训练数据具有充分代表性的条件下，因为交叉验证误差可能不适用于模型预测因子空间适用范围之外的区域。尽管存在这些局限性，本研究仍证明了紧凑的两次飞行设计在田间尺度大豆产量预测中具有研究内可行性。其主要贡献在于证明，在无信息泄露的外部测试评估条件下，季前裸土多光谱影像所包含的预测信息在大豆鼓粒期冠层单次采集之后依然具有实用价值。现有结果为进一步外部验证提供了明确基础，以确定该紧凑采集设计在不同环境和生长季节之间的迁移能力。

资源可用性

主要联系人

有关进一步信息和资源的请求应提交给主要联系人Zhiyuan Tian（[e-mail]），并由其履行。

材料可用性

本研究未产生新的独特材料。

数据与代码可用性

• 支持本研究的一个样本级数据集已存档于Figshare，并公开发布于https://doi.org/10.6084/m9.figshare.33130349。
• 与本研究所报告分析相关的R脚本已存档于Figshare，并公开发布于https://doi.org/10.6084/m9.figshare.33130349。
• 重新分析本文所报告数据所需的任何附加信息，可在提出合理请求后向主要联系人获取。

致谢

本研究由国家重点研发计划（项目编号：2024YFD1501101和2021YFD1500801）以及中国科学院战略先导重点研究计划（项目编号：XDA28010202）提供资金支持。我们衷心感谢Sitan Wang和Dongbo Li在调查和数据收集方面提供的帮助。

作者贡献

Y.Z.：调查、方法、软件、形式化分析、可视化、撰写——初稿撰写与撰写——审阅与编辑；Z.T.：资金获取、概念化、调查、撰写——初稿撰写与撰写——审阅与编辑；R.M.：调查与撰写——审阅与编辑；Y.L.：资金获取、监督与项目管理。

利益声明

作者声明不存在利益冲突。

方法（STAR★Methods）

关键资源表

试剂或资源 | 来源 | 标识符

存档数据 | 样本级数据集和R分析代码 | 本研究 Figshare：https://doi.org/10.6084/m9.figshare.33130349

软件与算法 | R（版本4.5.1） | R Core Team https://www.r-project.org/
ranger（版本0.17.0） | CRAN https://cran.r-project.org/package=ranger
kernelshap（版本0.8.0） | CRAN https://cran.r-project.org/package=kernelshap
terra（版本1.8–29） | CRAN https://cran.r-project.org/package=terra
sf（版本1.0–21） | CRAN https://cran.r-project.org/package=sf
exactextractr（版本0.10.1） | CRAN https://cran.r-project.org/package=exactextractr
ggplot2（版本3.5.2） | CRAN https://cran.r-project.org/package=ggplot2
DJI Pilot 2 | DJI https://www.dji.com/downloads/djiapp/dji-pilot-2
DJI Terra | DJI https://enterprise.dji.com/dji-terra

其他 | DJI Mavic 3 Multispectral | DJI https://enterprise.dji.com/mavic-3-m

方法细节

研究区域与田间实测数据采集

研究区域位于中国东北黑龙江省的两个地点：海伦（126°49′42″E，47°21′23″N）和北安（126°38′52″E，48°19′08″N）（图1）。该区域位于典型的黑土区，属大陆性季风气候，冬季寒冷干燥，夏季炎热多雨。年平均气温在?5°C至5°C之间，年平均降水量为400–650毫米。依据《中国土壤系统分类》，土壤类型为黑土；依据美国农业部土壤分类系统（USDA Soil Taxonomy），则归类为Mollisol（软土/黑钙土类）。研究区域的主要土地利用类型为小农经营的耕地，实行统一的玉米—大豆轮作。大豆通常于5月初播种，10月初收获，每年收获一季。

每个地点选取了两个坡面单元：海伦的HLA（n = 25）和HLB（n = 20），北安的BAA（n = 19）和BAB（n = 37），共101个采样点（图1D）。所选田块在区域气候与地貌背景方面具有可比性，并采用同一集中化管理系统。在每个坡面单元上，采样点以10米的间隔设置，用于土壤采集和产量测定（n = 101）。由于耕作效应主要局限于0–20厘米土层，所有土壤分析均针对该深度进行。土壤采样分别在2023年4月的海伦和2024年4月的北安进行。2023年采集的海伦土壤样本对应于2024年无人机观测和产量测定所使用的相同地理配准采样点和田块。在每个采样点均记录经纬度坐标。随后以该点为中心建立1米×1米的样方，并在样方内随机选取三个位置采集0–20厘米土层土壤子样本。这些子样本混合后充分混匀，形成每个样方的复合样本。每个复合样本风干后分为两份：一份通过2毫米筛网，采用吸管法测定砂、粉砂和黏粒组分；另一份研磨至通过0.149毫米筛网，采用重铬酸钾氧化—外加热法测定土壤有机碳（SOC）。选择土壤质地和SOC作为指标，是因为它们在较短时间尺度上保持稳定，能够在一定程度上可靠反映土壤状况。

大豆产量测定于2024年9月下旬完成。根据记录的地理坐标，在每个采样点中心设置一个产量样地，样地包含两行作物，行长1–1.44米，宽1米，平均样地面积为1.28平方米。每个样地内所有豆荚均手工采收，脱粒获得洁净籽粒并称取鲜重。籽粒在75°C下烘干至恒重后再次称重，以确定含水量。大豆产量按13%标准含水量折算，并以kg ha?1表示。实测土壤性质与大豆产量的描述性统计结果见表S1。

无人机数据采集与处理

多光谱影像由DJI Mavic 3 Multispectral无人机平台采集（深圳市大疆创新科技有限公司，中国广东深圳）。该平台配备一台2000万像素RGB相机和四个500万像素多光谱传感器。所采集的多光谱波段包括绿光（G，560 nm ± 16 nm）、红光（R，650 nm ± 16 nm）、红边（RE，730 nm ± 16 nm）和近红外（NIR，860 nm ± 26 nm）。集成的RTK模块为每张影像的曝光中心提供厘米级定位。集成的日照传感器会在影像文件中记录入射太阳辐照度，从而在正射影像重建过程中实现光照补偿。

共进行了两次无人机采集：一次在2024年4月下旬的季前裸土期（大豆出苗之前），另一次在2024年8月下旬、田间识别确定的大豆鼓粒期。以往研究表明，该生育期大豆的表型响应与最终产量关联最为密切。飞行作业在当地时间09:00至11:00之间进行，且要求风力和天气条件平稳、无云遮挡。飞行航线使用DJI Pilot 2软件进行规划与执行。飞行高度为30–50米，航向和旁向重叠度均设为80%，飞行速度为2.6–4.4 m s?1，相应空间分辨率约为2.5 cm·pixel?1。多光谱影像在DJI Terra中处理生成地理配准的正射影像，供后续光谱指数计算使用。目视检查表明，鼓粒期采集时大豆冠层几乎完全覆盖土壤表面；因此未对鼓粒期影像实施显式的土壤背景掩膜。

光谱指数计算与预测因子情景

最终候选池包含36个光谱指数。候选池仅限于已发表、可由DJI Mavic 3 Multispectral平台提供的绿光、红光、红边和近红外波段计算得出，且在作物冠层或裸土多光谱研究中已有使用或有其光谱依据的指数。36个指数的缩写、全称、计算公式及文献来源见S2表。

通过一个R脚本，将同样的36个公式分别应用于季前裸土和鼓粒期正射影像。对于每个地理配准的采样点，以采样点为中心建立轴对齐的1米×1米方形窗口。每个指数首先在像元水平上由绿光、红光、红边和近红外栅格值计算得出。随后将像元级指数值在采样窗口内取平均，其中与窗口边界相交的像元按其覆盖比例加权。该流程为101个观测中的每一个产生36个鼓粒期植被光谱指数和36个裸土光谱指数。由鼓粒期影像衍生出的指数此后记为植被光谱指数（VSI）；SSI表示裸土光谱指数预测因子组，这些指数被解释为土壤相关的谱学替代指标，而非土壤性质的直接测量值。

共评估了五种预测因子情景。VSI-only最初包含由鼓粒期影像衍生的全部36个指数。为构建土壤相关候选池，基于模型训练数据，使用随机森林排列重要性对36个裸土指数分别针对实测土壤有机碳（SOC）、砂粒、粉砂和黏粒进行排序。对于每一种实测土壤性质，保留排名前两位的指数以形成针对该指标的候选池。SSI-only使用四种土壤性质所识别候选者的并集。VSI+SOC将所选VSI子集与SOC相关候选池相结合；VSI+Texture使用砂、粉砂和黏粒相关候选池的并集；VSI+SSI使用所有四种实测土壤性质相关候选者的并集。在多个针对指标的候选池中重复出现的指数仅保留一次。实测土壤性质仅用作构建和评估相应裸土替代指标候选池的目标变量，并不直接作为大豆产量模型中的预测因子。

用于大豆产量建模的预测因子选择采用随机森林递归特征消除（RF-RFE）方法。在模型开发过程中，VSI-only所选的VSI子集在相应的VSI+SOC、VSI+Texture和VSI+SSI情景中保持不变，随后从相应的候选池中选取额外的裸土预测因子。所有土壤性质指导下的候选池构建和面向产量的预测因子选择仅使用模型训练数据完成，不涉及测试集信息。详细的重采样设计、预测因子排序流程、RF-RFE实现、子集规模选择及超参数调优方法将在下文"嵌套特征选择与随机森林建模"部分进行说明。

嵌套特征选择与随机森林建模

嵌套重采样与预测因子选择

模型开发采用10个固定的、按坡面分层的80:20训练—测试划分。每个划分包含81个外部训练观测和20个外部测试观测，其中外部测试集由BAA的4个观测、BAB的7个、HLA的5个和HLB的4个组成。所有候选池构建、预测因子选择、超参数调优及模型拟合均限于外部训练数据，外部测试数据则留作最终评估使用。此处将单一划分中的完整模型开发与评估流程统称为一次外部重复。

在每个外部训练集内，5折交叉验证重复进行两次。观测样本在每个坡面单元内被分别分配到5折中。在每个折中，四组作为内部训练数据，剩余一组作为内部验证数据。此处将每次重复与每一折所形成的训练—验证划分称为一次内部重采样，每次外部重复共产生10次内部重采样。

土壤相关的SSI候选池在面向产量的预测因子选择之前构建。对于每次内部重采样，仅使用内部训练数据，分别拟合随机森林以从36个裸土指数预测SOC、砂粒、粉砂和黏粒。采用排列重要性对指数进行排序，该指标量化了某一预测因子数值被随机重排后预测性能的下降幅度。每个筛选森林包含500棵树，每次分裂考虑6个预测因子，最小终端节点大小为7，并对内部训练观测按80%比例有放回抽样以构建每棵树。重要性值在每个森林内部进行归一化，并在10次内部重采样中进行汇总。对于每种土壤目标变量，指数按其在10次内部重采样中作为前两位高排名预测因子出现次数的多少进行排序，随后依据其平均归一化重要性及平均排名。每种目标保留排名前两位的指数并用于构建如"光谱指数计算与预测因子情景"中所述的候选池。每种土壤目标采用固定的"两位"截断标准，以在小样本嵌套设计中既保留多个高排名谱学候选项，又限制候选池维度；该截断标准统一应用于四种土壤目标，以避免针对不同目标对候选池规模进行专门调优。

随后采用随机森林递归特征消除（RF-RFE）进行面向产量的预测因子选择，以排列重要性作为消除判据。

在每个内部重采样中，首先对当前预测变量集拟合随机森林。在每个缩减步骤中，重新拟合森林，重新计算置换重要性，并移除排名靠后的预测变量以形成下一个更小的子集。使用内部验证RMSE评估允许的子集。RF-RFE森林包含500棵树，使用当前预测变量数量的平方根向上取整作为每次分裂时考虑的预测变量数，终端节点最小大小为7，以有放回方式抽取80%的内部训练观测值。仅VSI限制为1–3个预测变量，仅SSI限制为2–4个预测变量。对于每个组合场景，同一外部重复中为仅VSI所选的VSI子集被锁定，即在特征消除过程中保留这些预测变量，同时评估额外的SSI预测变量。这种锁定策略确保在每个外部重复中，每个组合场景都在相同的外部划分和相同的选定VSI基线下进行评估，从而使比较针对添加的裸土预测变量的增量贡献。VSI+SOC允许额外1–2个SSI预测变量，VSI+纹理允许1–3个，VSI+SSI允许2–4个，任何组合模型中不超过6个预测变量。

最终子集大小使用一个标准误规则确定。保留最小平均内部验证RMSE不超过最小平均RMSE加一个标准误的最小子集。确定子集大小后，预测变量按它们在10次内部重采样中的选择频率排名，然后按平均归一化置换重要性排名。选择排名最高的预测变量，在组合场景中还包括锁定的VSI预测变量。

随机森林调优与最终模型拟合

选定预测变量子集后，使用相同的10次内部重采样调优随机森林超参数。mtry（每次树分裂时考虑的预测变量数）的候选值为1、选定预测变量数量的平方根向上取整、以及完整的选定预测变量数量，移除重复值。评估的终端节点最小大小为3、7和15。每个调优模型包含400棵树。选择平均内部验证RMSE最低的参数组合。并列时依次按RMSE标准差更低、mtry更小、终端节点最小尺寸更大解决。为选定组合生成内部折外（inner OOF）预测，即观测值被排除在相应内部模型拟合之外时生成的预测。由于重复了两次5折交叉验证，对每个外部训练观测值的两个内部OOF预测取平均。随后使用选定的预测变量和调优后的超参数，对所有81个外部训练观测值拟合包含2000棵树的最终随机森林模型，并对相应的20个外部测试观测值应用一次。每个树使用以有放回方式抽取的80%外部训练观测值。光谱指数预测变量以其原始提取单位使用，未经标准化。

分别拟合SOC、砂、粉砂和粘土的土壤代理标定模型，每个目标保留两个裸土指数。这些模型遵循相同的内部重采样、调优和最终拟合程序，但预测相应的测量土壤性质而非大豆产量。其外部测试预测仅用于评估代理标定；测量的或预测的土壤性质值均未输入大豆产量模型。

量化与统计分析

模型验证与空间敏感性分析

性能评估与配对场景比较

使用内部折外（inner OOF）和外部测试预测评估模型性能，以外部测试性能作为主要估计。计算五个指标：预测R2（定义为1?SSE/SST，相对于相应评估集中的观测均值）；Pearson r2（定义为观测值与预测值之间Pearson相关系数的平方）；均方根误差（RMSE）；平均绝对误差（MAE）；以及相对RMSE（rRMSE，计算为RMSE除以观测均值，以百分比表示）。当预测误差超过使用评估集均值时的误差时，预测R2可能为负值。产量RMSE和MAE以kg ha?1表示。土壤代理标定模型使用相同指标；SOC的RMSE和MAE以g kg?1表示，砂、粉砂和粘土以百分点表示。结果总结为10次外部重复的平均值和标准差。统计分析在R 4.5.1版本中进行，随机森林模型使用ranger包（版本0.17.0）拟合。研究级别和分析特定的样本量以及汇总统计的定义在METHOD DETAILS、结果部分和相应的表格与图注中提供。

由于所有五个预测变量场景使用相同的外部划分，每个组合场景在同一外部重复中与仅VSI进行比较。配对差值计算为组合场景减去仅VSI；正差值表示预测R2改善，负差值表示RMSE、MAE和rRMSE改善。报告有利的外部重复次数。每个平均差值的95%置信区间通过对10个配对差值进行10000次非参数Bootstrap重采样获得。使用精确的配对符号翻转置换检验计算双侧p值，该检验在零假设（无系统性配对差值）下评估所有21?种可能的符号分配。效应量、置信区间和有利重复次数与p值一起考虑。

产量分层残差与空间敏感性

检查仅VSI和VSI+SSI在低、中、高产层的残差模式。在每个外部重复中，使用外部训练集中观测产量的33.3百分位和66.7百分位数定义分层，并将所得阈值应用于相应的外部测试观测值。偏差计算为预测减去观测，正值表示过预测。每个层总结偏差、RMSE和MAE。计数代表跨10次外部重复合累的预测记录；由于划分重复，同一观测值可能贡献多个记录。

使用留一坡度单元（LOSO）验证评估空间敏感性。在四次LOSO折叠中，每次排除一个完整的坡度单元（BAA、BAB、HLA或HLB），并使用其余三个坡度重复候选池构建、RF-RFE、超参数调优和模型拟合。为每个留出的坡度计算性能，并跨101个唯一的LOSO预测汇总。重复的坡度分层外部划分评估采样坡度单元内新观测值的预测，而LOSO评估向完全未见坡度单元的迁移。因此两种验证设计不视为配对估计。

对于随机划分和LOSO模型中选定的预测变量，使用绝对标准化均值差（|SMD|）和Kolmogorov–Smirnov（KS）统计量描述训练–测试分布差异。前者衡量训练和测试均值相对于其合并变异的差异，后者衡量其经验累积分布之间的最大分离。较大的值表示较大的分布差异。这些统计量仅作描述性使用，无假设检验p值。随机划分结果总结跨选定特征–重复记录，LOSO结果总结为每个预测变量在其被选中的折叠中。

外部测试预测的SHAP分析

使用SHapley加性解释（SHAP）描述选定预测变量如何对拟合的随机森林预测做出贡献。SHAP值表示一个预测变量使预测高于或低于基线预测的程度。对所有五个预测变量场景跨10次外部重复进行SHAP分析，共50个最终模型。LOSO模型被排除，因为LOSO用于空间敏感性分析而非主要模型解释框架。

对每个模型，使用相应外部重复中选定的预测变量和超参数重构最终随机森林。重构的外部测试预测与存储的正式预测进行验证。81个外部训练观测值用作SHAP背景数据，定义基线预测和计算预测变量贡献的参考分布，SHAP值仅对相应的20个外部测试观测值计算。因为所有最终模型不超过6个预测变量，所以计算精确SHAP值。使用kernelshap包（版本0.8.0）实现精确KernelSHAP用于两个或更多预测变量的模型，对单预测变量模型使用等效的直接分解。通过确认基线加所有预测变量特定SHAP值重现模型预测来验证可加性。SHAP值以kg ha?1表示，与预测产量单位相同。

在每个外部重复中，全局重要性计算为每个选定预测变量在外部测试观测值上的平均绝对SHAP值。由于预测变量集在外部重复间不同，使用选择加权平均绝对SHAP值总结重要性。对于该度量，预测变量被选定时贡献其平均绝对SHAP值，未被选定时贡献零，然后跨所有10次外部重复取平均。该总结因此同时纳入选择频率和贡献幅度。其95%置信区间使用重复级值的10000次Bootstrap重采样估计。

对于VSI+SSI，使用蜂群图显示带符号SHAP分布，每个点代表一个外部测试预测中的一个选定预测变量。正SHAP值增加相对于基线的预测，负值减少。预测变量值在每个预测变量内单独标准化，仅用于定义点颜色；因此颜色在预测变量内解释，而非跨预测变量。

为六个VSI+SSI外部测试预测准备局部SHAP分解，代表在Beian和Hailun两个地点的低、中、高产层。此处"局部"指单个预测的分解，而非空间尺度。在每个地点内，使用客观规则从每个产量层选择一个独特样本，该规则最小化从层特定中位绝对预测误差的综合距离；平局按样本ID和外部重复解决。每个分解显示预测变量特定SHAP值如何将模型基线转移到最终预测。SHAP值描述对拟合模型预测的贡献，未被解释为对观测产量的因果效应。

补充信息

（1）PDF（685.10 KB）文档S1。图S1–S4和表S1–S9