植物碳通量在气候变化条件下影响着土壤有机碳的动态变化：机器学习揭示了关键的全球植物生产力（GPP）阈值

《Ecological Indicators》：Plant carbon fluxes govern soil organic carbon dynamics under climate change: Machine learning reveals critical GPP thresholds

　　黄飞宁|张永坤|卢星星|尚冠伟|李 Qingliang|秦章才|魏中旺|袁欢|卢莉|戴永久|

**广东省气候变化与自然灾害研究重点实验室，中山大学大气科学学院，珠海519082**

**摘要**
土壤有机碳（SOC）的动态受到气候变化的深刻影响，然而植物碳通量、气候强迫和土地利用变化的贡献仍然缺乏定量研究。在此，我们开发了一个机器学习框架，用于评估1980年至2100年间中国0-20厘米和20-100厘米深度土壤中SOC动态的驱动因素。该模型通过10折交叉验证进行了验证，对于100厘米和20厘米土层分别获得了0.41和0.48的R2值。此外，在两个不同时间段进行的独立时间验证也显示了该模型能够准确捕捉SOC存量的空间分布和时间动态。根据四种共享社会经济路径（SSP）情景下的未来预测，到2100年，0-100厘米深度的全国SOC将增加4.8至7.7Pg C。敏感性分析表明，植物碳通量是控制SOC变异性的主要因素。Shapley Additive Explanations（SHAP）分析进一步揭示了总初级生产力（GPP）与SOC之间的非线性关系，这在一定程度上抵消了由气候变化引起的SOC损失。这种关系具有两个关键阈值：首先，在GPP约为4 gC m-2 d-1时，GPP对SOC的边际贡献（SHAP值）从负值转变为正值，表明生态系统从净碳亏损状态转变为盈余状态，这与等降水量线400毫米的生态敏感性区域相吻合；低于这一阈值时，不足的输入会促进已有SOC的分解（可能是启动效应）。其次，在GPP约为7 gC m-2 d-1时，GPP-SOC的正敏感性达到平台期，表明随着生产力的进一步提高，碳封存效果减弱。在高排放情景下，这一转变区向西北方向扩展到历史上干旱地区。

**1. 引言**
土壤有机碳（SOC）是最大的陆地碳库，其储存的碳量超过大气和全球植被中的总碳量（Hengl等人，2017；Jobbágy和Jackson，2000；Todd-Brown等人，2013）。因此，即使SOC储量的微小波动也会显著影响区域碳通量，并放大全球气候变化的影响（Jungkunst等人，2022；Potter等人，1993）。作为土壤碳封存的核心，SOC的组成、来源和稳定机制对于理解陆地碳汇的功能以及制定有效的气候变化缓解策略至关重要（Lal，2004；Jackson等人，2017；Bossio等人，2020；Bai和Cotrufo，2022）。在当前的气候变化下，人类活动（包括直接干预如造林/毁林（Fujisaki等人，2017；Shao等人，2019）和间接效应如温度驱动的植被生长（Li等人，2024a；de Mello等人，2025）通过改变植被状况深刻影响了SOC动态。总初级生产力（GPP）衡量的是植被通过光合作用每单位面积和时间的总碳固定量，它是进入陆地生态系统的最大碳通量，并在碳循环中起关键作用（Beer等人，2010；Pandey等人，2024；Amir等人，2025；Xu等人，2025）。因此，GPP-SOC之间的联系有效地说明了生态植被变化如何影响SOC的储存。这一联系对中国尤为重要，因为那里复杂的植被模式（Su等人，2022）和大量的SOC储量因造林政策（Hong等人，2020）而经历了显著的GPP变化。这些变化有可能重新配置SOC动态，对国家碳预算和全球气候反馈产生连锁影响。

GPP-SOC关系通过几个相互关联的生物地球化学过程体现：包括通过植物凋落物（Fan等人，2023）和根系分泌物（Chari等人，2024）的碳输入、植被介导的C/N比率变化影响分解速率（Mo等人，2022）、以及森林扩张导致的木质素含量增加（Yang等人，2021），这些都有助于长期SOC的储存。然而，气候变化同时通过温度、降水、大气CO?浓度和土地利用的变化重塑了这些联系。虽然升温可以增加某些地区的GPP（He等人，2022）和碳输入，但也会加速微生物活动（Crowther等人，2016；García-Palacios等人，2021）和SOC分解（尤其是在永久冻土中，He等人，2023；Liu等人，2024），这可能抵消部分收益（Beillouin等人，2023；Griscom等人，2017）。降水变化（表现为湿干对比度的加剧）进一步复杂化了这一平衡：极端降雨事件可能引发土壤侵蚀和SOC损失（Jung等人，2014；Su等人，2024），而干旱则可能促进微生物碳输入对SOC的贡献（Wang等人，2024）。同样，二氧化碳浓度升高可以增加根系分泌物和地下碳的分配（Norby等人，2024），但也可能激活微生物分解现有SOC（Li等人，2022a）。这些相互竞争的机制引入了关于未来气候驱动的GPP变化下SOC变化方向的根本不确定性，突出了需要综合、过程意识的建模框架来确信陆地生态系统在未来几十年是否会继续作为碳汇或成为碳源。

当前挑战阻碍了我们确定气候驱动的GPP变化下SOC变化净方向的能力。传统的基于过程的模型（PBM）由于对SOC形成、稳定和分解过程表示不一致，常常产生不同的未来SOC预测（Abramoff等人，2018；Georgiou等人，2021）。此外，它们依赖于高分辨率、空间明确的强迫数据进行校准，这也构成了另一个障碍（Pierson等人，2022）。PBM理论上通过明确模拟机械路径来量化植被碳输入与SOC之间的耦合，这些路径可以通过操纵实验进行验证。然而，这些关系通常是基于PBM内的理论参数化制定的，而不是基于观测数据的。近年来，机器学习（ML）算法在SOC建模中取得了显著进展，特别是在数字土壤制图（DSM）框架内（Bernardini等人，2024；Shi等人，2024；Li等人，2024b；Peng等人，2025）。基于DSM的未来SOC预测通常依赖于空间对时间的替代假设（Pickett，1989），该假设通过空间梯度推断时间趋势，这一策略已应用于欧洲（Yigini和Panagos，2016）、中国（Zhang等人，2023）和阿根廷（Heuvelink等人，2021）。在ML算法中，随机森林（RF）已经成为最流行的SOC预测方法，证明了其在这一领域的有效性（Lamichhane等人，2019）。为了量化GPP-SOC联系，通过SHAP值解释的RF在统计归因研究中最为常用，能够提供对非线性和交互效应的透明、可解释的洞察（Huang等人，2023；Li等人，2023）。虽然PBM理论上可以表示植被碳输入与SOC之间的耦合，但其功能形式通常是基于理论参数化制定的，缺乏在大陆尺度上的经验基础。同样，现有的ML研究大多将GPP视为另一种预测因子，而没有明确模拟动态的非线性GPP–SOC联系或解决历史观测与未来ESM预测之间的结构不一致性问题。

为了克服这些限制，我们引入了一个新颖的、基于观测的框架，从三个基本方面推进了植被-土壤碳耦合的诊断。首先，我们建立了基于观测准确性的实证GPP–SOC基线。通过将高质量土壤剖面与通量塔约束的GPP相结合，我们得出了一个物理上可解释的关系，规避了理论模型的结构不确定性。与以往仅基于历史数据拟合模型（面临外推风险）或使用未经校正的ESM输出（引入结构偏差）的ML研究不同，我们的研究引入了一个受观测约束的、动态一致的框架来预测气候变化下的SOC。其核心创新在于对CMIP6 ESM碳通量应用了非线性的、保留异常的偏差校正。这一程序将未来的GPP轨迹锚定在中国土壤剖面和通量塔观察到的实证GPP–SOC关系上，从而确保了历史与未来之间的物理连续性。因此，我们的ML模型不仅仅是从静态协变量预测SOC；它模拟了在新气候状态下植被生产力和土壤碳之间动态耦合的演变。其次，随后的基于SHAP的分析（Shapley Additive Explanations）探究了这个物理上连贯的系统，提供了关于阈值和转变的稳健见解。这样，我们的方法弥合了数据驱动的灵活性和基于过程的真实性之间的差距——这是对传统PBM和现有ML评估的关键改进。除了将GPP视为简单的预测因子外，我们还利用SHAP值来量化GPP–SOC耦合的数学和生物学特性：（1）我们将功能联系强度定义为每个驱动因素对SOC变异性的相对贡献，通过其平均绝对SHAP值的比例来量化；（2）利用历史基线，我们将SHAP值的符号作为局部碳预算偏差的指标。负SHAP值表示“碳输入亏损”，即某个驱动因素相对于全国平均水平而言起到生物物理瓶颈的作用；相反，正SHAP值表示“碳输入盈余”，表明该变量使SOC储存超过大陆平均水平；（3）特征值与其SHAP组分之间的依赖性提供了SOC边际响应的稳健表示。这使我们能够明确模拟动态的非线性GPP–SOC联系，识别出如400毫米等降水量线附近的敏感气候阈值，这些阈值决定了碳亏损和盈余状态之间的转变。

**2. 材料与方法**
**2.1. 研究区域**
中国位于欧亚大陆东部，西太平洋沿岸，拥有广阔且地形复杂的领土，涵盖了显著的纬度和经度梯度。明显的陆海对比和剧烈的地形变化产生了多样的气候类型，包括热带季风、亚热带季风、温带季风、温带大陆性和高山高原气候。这些气候梯度塑造了独特的植被模式和生物群落，从而对区域SOC储量产生了强烈影响。因此，SOC表现出与纬度和经度一致的系统空间模式（Chen等人，2024）。寒冷湿润的气候通常通过抑制微生物分解促进SOC积累，这种模式在高海拔地区尤为明显，那里低温和充足的水分减缓了分解速率。因此，中国最高的SOC浓度通常出现在东北部山区森林和青藏高原东部，而西北部的干旱和半干旱地区则表现出较低的SOC浓度（Ai等人，2025；Liu等人，2022；Zhuo等人，2022）。这种环境异质性产生了跨气候、土壤和生物条件的SOC储量综合自然梯度。因此，中国成为研究SOC动态驱动因素以及开发整合多种环境因素的预测模型的理想自然实验室。

**2.2. 数据与处理**
**2.2.1. 观测数据**
为了使用观测历史数据构建SOC模型，我们收集了包括土壤剖面、静态变量和动态变量在内的多种数据。历史数据覆盖了1977-2015年的期间。所有数据集的完整概述，包括其来源、时间覆盖范围、空间分辨率以及针对每个数据集应用的特定预处理方法，见表1。第二次全国土壤调查（SNSS）代表了一项具有里程碑意义的全国性努力，其中专业技术人员团队与当地农民一起，在2444个县、312个国营农场和44个林业中心系统地研究了土壤。这项行动产生了超过25,000个具有代表性的土壤剖面图，其中包含了每个土壤层关键参数的详细测量数据。为了解决这些历史剖面图中缺乏精确地理坐标的问题，我们应用了一种强大的地理参考校正程序（Shi等人，2024年），将每个剖面图记录的土地利用、海拔和坡度信息与高分辨率的空间协变量层对齐。从这个丰富的档案中，我们提取了土壤有机质数据，并将其转换为SOC（土壤有机碳），作为我们分析的历史基准。

图1. 数字高程图（DEM）叠加了第二次全国土壤调查（1979-1984年，蓝点）和中国陆地生态系统碳密度数据集（2000-2014年，红点）的土壤剖面采样位置的空间分布。较新的CTEC数据集通过提供包括森林、草地、农田、湿地和灌木丛在内的主要生态系统类型的广泛覆盖范围，补充了这些历史记录。它记录了植被碳密度（地上和地下）以及两个深度区间（0-20厘米和0-100厘米）的SOC密度。该数据集结合了原始的实地测量数据和经过严格审查的文献数据。对于本研究，我们使用了其7,683条SOC密度记录的子集，其中包括4,536个0-20厘米层的测量值和3,147个0-100厘米层的测量值。

(2) 静态变量
本研究使用地形、母质和土壤因素作为静态变量，因为它们在时间上的相对稳定性（表1）。具体来说，地形和母质在长时间内基本保持不变，而土壤性质如含沙量在世纪时间尺度上的变化很小。数字高程模型（DEM）来源于去除了多种误差的改进地形DEM（Yamazaki等人，2017年），而其他地形参数如坡度则是使用SAGA GIS工具从这个基础DEM计算得出的。对于土壤母质的表征，我们利用了美国地质调查局地球资源观测与科学中心的全球岩性数据（Hartmann和Moosdorf，2012年）。土壤参数结合了中国特定的土壤类型分类数据和基岩深度地图，并补充了从Harmonized World Soil Database（Wieder，2014年）中提取的土壤性质信息。这个综合数据集为历史时期和未来预测时期的这些静态环境控制因素提供了一个一致的框架，确保了长期环境条件的稳健表示。

(3) 动态变量
气候数据（表1），包括每月最高和最低温度以及降水量，来源于WorldClim2数据库（Fick和Hijmans，2017年），为基准时期（1970-2000年）提供了高分辨率（约1公里）的历史平均值。土地利用数据来自CMIP6土地利用协调2（LUH2）数据集（Hurtt等人，2020年），空间分辨率为0.25°×0.25°，包含了14种不同土地利用类别的面积百分比的历史时间序列。这些数据涵盖了原始森林和非森林土地、潜在的次生森林和非森林区域、管理的牧场、天然草地以及C3和C4一年生/多年生作物、C3固氮作物，还包括次生植被的年龄和地上生物量。具体来说，这些变量不是直接测量得到的，而是通过结合卫星观测到的森林覆盖损失数据和基于FAO（粮食及农业组织）林业数据的簿记模型得出的。历史植物碳通量数据包括两个关键变量：GPP（总初级生产力）和NEP（净生态系统生产力）（Shangguan等人，2023年）。GPP表示单位面积和时间内植被光合作用固定的总大气碳，而NEP表示在考虑了植被的自养呼吸和土壤及其他生物的异养呼吸后的剩余碳。这些数据集为我们时空分析提供了必要的气候、土地利用和生物地球化学输入。为了确保与ESMs（地球系统模型）在结构上的一致性，我们建模框架中使用的是LUH2数据集中的土地利用数据，该数据集为CMIP6 SSPs（共享社会经济路径）提供了标准化的土地利用驱动因素。为了将其纳入我们的建模框架，这些数据被转换成了随时间变化的分类变量：对于每个网格单元和时间段（例如，1980-1999年、2000-2015年、2021-2040年）。这允许RF模型学习特定于土地利用的SOC动态，并捕捉随时间变化的土地利用效应。

为了使环境驱动因素与土壤观测结果一致，所有协变量数据都根据土壤剖面采集周期进行了时间同步。具体来说，1980-1999年和2000-2015年间的气候、土地利用和碳通量变量分别与SNSS（1979-1985年）和CTEC（2000-2014年）的相应土壤数据集进行了匹配，确保了建模输入的时间一致性。

2.2.2. 地球系统模型数据
对于未来的预测，我们从参与第六阶段耦合模型比较项目（CMIP6）的四个地球系统模型（ESMs）中获取了数据，分别是ACCESS-ESM1-5（Ziehn等人，2020年）、EC-Earth3-Veg（D?scher等人，2022年）、IPSL-CM6A-LR（Boucher等人，2020年）和MPI-ESM1-2-LR（Mauritsen等人，2019年）。选择这些特定的ESMs是基于最近的基准研究，这些研究表明它们在模拟历史土壤碳储量和初级生产力方面的高保真度（Hu等人，2022年；Varney等人，2022年；Lu等人，2023年）。模型输出涵盖了四个共享社会经济路径（SSPs）：SSP126、SSP248、SSP370和SSP585，时间跨度为2021年至2100年的二十年时间段。这些模拟得出的动态协变量包括气候变量、土地利用模式和碳通量。这组ESM-SSP组合提供了在不同未来气候和社会经济轨迹下SOC动态的多方面评估。

2.2.3. 数据处理
所有协变量都被重新采样到统一的2.5角分分辨率，以确保输入数据集在空间上的一致性。

(1) 土壤有机碳
SNSS记录了土壤有机质（SOM）浓度（重量百分比）。我们使用经验系数0.58将SOM转换为SOC。SOC储量（SOCS，克/平方厘米）通过以下公式计算：
SOCS = SOM × 0.58 × D × BD × 1?G
其中D表示土壤层厚度（厘米），BD是容重（克/立方厘米），G是砾石含量（大于2毫米，体积百分比）。BD和G的数据来源于Shangguan等人（2023年）。

为了实现跨剖面比较和建模的一致性，我们将原本由不同厚度层（通常范围从5到100厘米）和最大深度（集中在20厘米和100厘米；见图S1）组成的所有SNSS测量值，使用等面积二次样条插值（Odgers等人，2012年）统一到两个标准深度区间（0-20厘米和0-100厘米）。为了减少相邻层次之间突变引起的异常，当连续层次之间的绝对比率超过1.225时，我们插入了1厘米的过渡层，从而防止重建的剖面出现不现实的波动。在计算SOC储量之后，我们应用了3σ原则来识别和移除统计异常值，假设数据呈正态分布，并移除了偏离平均值三个标准差以上的338个数据点，以确保后续分析的基准测量结果的稳健性。

(2) 未来碳通量
所有碳通量数据（GPP和NEP）被汇总成季节性的季度数据集（春季、夏季、秋季和冬季），以捕捉生态系统碳动态的年内变化。GPP表示通过光合作用流入的碳的总量，NEP表示包括异养呼吸在内的碳同化和生态系统整体损失的结果，两者共同提供了生态系统碳平衡的全面表示。

为了处理CMIP6 ESMs之间历史和未来碳通量预测之间的显著差异（Lu等人，2023年；Luo等人，2012年），这可能导致转换点出现人为的不连续性，我们实施了一种基于delta-change（异常保留）方法的偏差校正程序，并进行了非线性偏差校正的调整（Xu等人，2021年）：
CFi,j* = CFhistory + CFi,j?CFhistory,j
其中CFi,j*是来自ESM j的未来时期i的校正后的碳通量，CFhistory是来自Shangguan等人（2023年）的观测约束下的历史碳通量，CFi,j是来自同一ESM j的未来时期i的原始预测通量，CFhistory,j是同一ESM j在2000-2014年期间的历史输出。这种调整通过将未来预测基于一致的观测基准，同时保留每个模型内部的气候驱动趋势和年际变异性，从而消除了系统特定的偏差。这种校正考虑了GPP和NEP在模型表示上的差异性，确保了我们协变量时间序列的时间连续性。这减少了碳循环动态中的虚假跳跃，特别是在观测限制数据过渡到模型预测未来的关键点。

2.3. SOC模型构建和预测
我们开发了RF模型来使用两个历史参考时期的数据预测SOC：1980-1999年和2000-2015年。回归矩阵是通过将土壤剖面测量数据与相应的时空协变量配对构建的。我们提取了空间上对应的协变量值来构建包括土地利用协变量在内的统一回归矩阵。然后根据重要性得分对变量进行排名，对于高度相关的变量对（Pearson的|r| > 0.85；见图S2），通过保留重要性较高的变量进行修剪。经过迭代重新排序和模型运行的平均后，选择了最具影响力的前30个协变量（表S1）用于最终建模。模型的准确性通过十折交叉验证进行评估，使用了三个指标：决定系数（R2）、均方根误差（RMSE）和平均绝对误差（MAE）。为了测试我们用于SOC映射的模型的稳健性，我们应用了时间独立的验证方法，即在每次折叠验证中测试SNSS样本和CTEC样本。

我们的实验涵盖了六个时间框架：两个历史时期1980-1999年和2000-2015年，以及四个未来的二十年时间段（2021-2040年、2041-2060年、2061-2080年和2081-2010年）。遵循空间换时间替代策略（Blois等人，2013年；Liu等人，2020年），我们将这六个时期的动态协变量作为顺序输入到模型中。这种方法使我们能够重建历史SOC分布，并在不同ESMs和SSPs之间预测未来的SOC，同时保持模型结构和预测器空间的一致性。

2.3.1. SOC模型构建
为了评估我们建模框架的时间稳健性，并确定最适合SOC时空预测的策略，我们使用不同的训练数据集开发并比较了三个独立的RF模型：（i）仅使用1980-1999年的数据，（ii）仅使用2000-2015年的数据，以及（iii）结合两个时期的数据（1980-2015年）。RF是一种集成学习方法，结合了装袋算法和决策树，通过在回归任务中平均多棵树的输出来汇总预测。它在处理高维数据时表现出稳健性能，具有固有的抗过拟合能力，并内置了处理缺失值的机制。

由于计算限制，如果输入全部162个协变量（动态变量：36个气候变量、12个碳通量变量、26个土地利用变量），将需要超过1TB的内存，因此我们进行了特征选择。对于每个历史土壤采样时期（1980-1999年和2000-2015年），我们提取了空间上对应的协变量值来构建统一的回归矩阵。然后根据重要性得分对变量进行排名，对于高度相关的变量对（Pearson的|r| > 0.85），通过保留重要性较高的变量进行修剪。经过迭代重新排序和模型运行的平均后，选择了最具影响力的前30个协变量（表S1）用于最终建模。模型超参数通过R语言中的caret包进行了十折交叉验证进行优化，分别针对三个训练数据集进行了优化。我们对关键调整参数mtry（每次树分割时考虑的随机采样候选变量的数量）进行了网格搜索，使用R2作为优化标准。0-20厘米SOC模型的最佳mtry值为5，0-100厘米模型的最佳mtry值为11，分别在每个模型深度和训练数据集配置之间达到了模型复杂性和泛化能力之间的有效平衡。

2.3.2. 模型不确定性量化
与第2.3.1节中描述的基于十折交叉验证的模型性能评估不同，我们的不确定性量化直接来自预测模型本身。DSM中的误差可能来源于数据流程的多个阶段，包括土壤测量、数字化、数据输入、解释、分类、泛化和插值。此外，不确定性还会通过建模选择传播，例如算法偏差、参数化和输入协变量的不准确性，进一步影响预测的可靠性。在许多情况下，特别是在农业管理和环境政策制定中，量化预测不确定性与进行预测本身一样重要。因此，不确定性图对于做出明智的决策至关重要。在DSM中，不确定性分析对于评估预测土壤图是否足够可靠以用于实际应用至关重要（Nenkam等人，2024年）。此外，明确表征不确定性是迈向模型透明度和可解释性的关键步骤。我们特别关注预测区间，即可能包含真实或未观测到的值的概率范围。然而，只有大约30%的研究量化了预测的不确定性（Wadoux等人，2020年）。为了解决这一差距，我们采用了分位数随机森林（quantile RF）进行不确定性分析，以获得第10、50和90百分位的预测结果（Schmidinger和Heuvelink，2023年；Vaysse和Lagacherie，2017年）。根据这些结果，我们使用以下公式计算了相对不确定性指标：(3)相对不确定性 = (Q90 - Q10) / (Q50) × 100%，其中Q90、Q50和Q10分别是第10、50和90百分位的预测值。这种标准化的度量方法有助于在不同土壤有机碳（SOC）量级的区域之间进行不确定性比较。

2.4. 将分组环境因素归因于SOC动态
为了量化主要环境驱动因素对未来SOC变化的相对贡献，我们将影响分为三个组：气候、碳通量和土地利用。我们采用了因子扰动方法来隔离每个驱动因素组的总体效应，而不是单个变量，以避免破坏每个组内紧密耦合的变量之间的内部共变。具体来说，对于每个驱动因素组，我们进行了一系列受约束的模拟，在这些模拟中，该组内的所有变量都被固定在其历史基线值（2000-2015年的平均值），而所有其他驱动因素则根据未来情景动态变化。这有效地将选定的驱动因素组转换为静态输入，模拟了一个“无变化”的反事实情景。根据第2.3.1节建立的建模框架，我们在每种受约束的配置下生成了不同的SOC预测。然后，每个驱动因素组的边际贡献被计算为全模型预测（所有驱动因素动态变化）与相应固定驱动因素预测之间的差异。为了便于解释，我们反转了这些差异的符号，使得正值表示对SOC储存的增强效应，反之亦然。这些差异量化了每个变量类别对预测SOC储量的净影响，提供了一个严谨的归因框架，区分了在不同社会经济路径（SSPs）下气候、生物化学和土地利用变化在控制未来陆地碳动态方面的相对重要性。

2.5. GPP-SOC关系的量化
2.5.1. 随机森林重要性
我们使用RFI解释方法量化了预测因子对SOC的相对影响。该指标通过测量在故意破坏给定预测因子与响应变量之间的关联时模型性能的下降来评估特征的重要性。这种方法通过以下步骤操作：(1) 重排：随机打乱每个预测变量的值，同时保持所有其他变量不变，从而打破被扰动特征与目标变量（SOC）之间的原始关系；(2) 影响评估：使用重排后的数据与原始数据集相比，测量模型性能的下降幅度。较大的增加表示更高的预测相关性；(3) 基于方差的计算：具体通过R语言中的ranger包和impurity参数设置实现，该参数计算由于随机排列干扰引入的预测响应的方差；(4) 坚健的平均：为了考虑排列过程中固有的随机变异性，重要性分数是在10次独立重复实验上的平均值，从而得到稳健的排名。

2.5.2. 随机森林的Tree-SHAP
为了进一步解释环境协变量与SOC之间学习到的关系，我们应用了SHAP（Shapley Additive Explanations），这是一个基于合作博弈论的模型不可知的框架（Lundberg和Lee，2017年）。对于每个土壤剖面或网格单元，特定特征$i$的边际贡献计算为所有可能的特征子集S?N?i的加权平均值：φi=∑S?N?iS!M?S?1!M!fxS∪i?fxS，其中Mis是特征的总数，fxS=EfxxS表示在特征子集S条件下的预期模型输出。这种方法满足了可加性属性，允许将模型的最终预测fx分解为基线值和各个特征贡献的总和：fx=EfX+∑i=1Mφi，这里EfX是整个训练数据集上的SOC预测的预期（平均值），作为我们碳预算分析的历史基线。由于我们的建模框架基于随机森林，我们使用了Tree-SHAP，这是一种优化算法，它利用决策树的内部结构来精确高效地计算SHAP值（Lundberg等人，2020年）。Tree-SHAP通过递归跟踪训练样本通过树节点的流动来估计条件期望EfxxS。对于一个分割特征i?S的节点j，算法将其计算为其子节点的加权和：EfxxSnodej=∑k∈LRrkrjEfxxSchildkj，其中rj表示到达节点j的训练样本数量，L和R分别代表左子节点和右子节点。这种结构集成允许算法捕捉高阶特征相互作用，计算复杂度为OT·L·D2，其中T、L和D分别代表树的数量、叶子和最大深度。我们使用R语言的iml包计算Shapley值，保持与我们SOC映射模型中使用的协变量的一致性。鉴于SHAP计算的高计算需求，我们实施了空间聚合策略：将研究区域划分为10×10的网格块，并在每个块内对空间平均的协变量计算SHAP值。这保留了主导的区域模式，同时实现了国家尺度上的可行计算。得到的SHAP图显示了关键驱动因素对SOC的方向性和相对影响。然而，需要注意的是，SHAP值反映的是边际效应而非机制关系。因此，应结合领域知识来解释它们。

3. 结果
3.1. 模型性能和验证
为了验证我们建模框架在时空SOC预测方面的稳健性，我们首先比较了三种不同的训练策略的性能（图2），然后基于最佳模型配置，我们通过空间不确定性分析对模型进行了验证（图3）。我们评估了三种不同的配置：(i) 仅在1980-1999年数据上训练的模型（“1980”模型），(ii) 仅在2000-2015年数据上训练的模型（“2000”模型），以及(iii) 综合时期的模型（“1980+2000”模型）。

首先，该框架表现出了一定的时间稳定性。单个时期的模型显示了不同的错误特征，反映了它们各自数据集的性质。历史上的1980-1999年模型表现出较低的精确度（较低的R2和较高的RMSE/MAE）以及较低的标准差，这可能表明与SNSS的方法论限制有关。相比之下，2000-2015年模型显示出较小的偏差，但随机误差较大。通过整合这两个不同的时期，综合时期的模型取得了显著的性能提升；数据集的互补性质允许部分抵消特定时期的偏差，从而显著降低了RMSE和MAE，并稳定了误差分布。关键的是，综合时期的模型表现出显著的时间稳定性和一致性。独立的时间验证（图2d，h）显示，综合模型分别在1980年代和2000年代的保留样本上进行了测试，表明回归斜率和准确性指标在两个不同的十年间保持基本一致。这种跨时期的稳定性表明，RF框架成功地识别了超越特定时间间隔的稳态生物物理关系，尽管SOC空间模式和环境协变量发生了变化。通过纳入这两个时期的更广泛的环境状态，模型在未来气候和植被条件下的外推能力得到了根本加强，这些条件可能会与历史平均值有显著偏差，从而最小化了未来预测中的潜在系统偏差。鉴于其在偏差和方差方面的优越平衡，我们采用了综合模型作为SOC映射和深度区间未来预测的基础。综合模型（图2中的“1980+2000”）在0-100cm预测上的R2为0.41，RMSE为0.30 g C/cm2，MAE为0.22 g C/cm2；在0-20cm层上的R2为0.48，RMSE为0.15 g C/cm2，MAE为0.10 g C/cm2。这些准确性指标处于最近SOC映射研究中报告的范围内（表S2；Li等人，2022a；Liang等人，2019a；Song等人，2020a；Zhang等人，2024b）。这反映了我们庞大的国家数据集（SNSS和CETC）中固有的高环境方差与模型在大陆尺度上解决SOC动态主导生物物理驱动因素的能力之间的稳健权衡。

其次，空间不确定性分析进一步确认了该框架的可靠性。图3显示了使用第2.3节描述的分位数回归森林（QRF）量化出的两个历史时期和两个深度区间（0-20cm和0-100cm）内建模不确定性的空间模式。中国西北部的不确定性相对较高，特别是在青藏高原。这种模式可能是由于土壤采样密度稀疏和协变量空间中高山土壤过程表示有限所致。值得注意的是，0-20cm层的不确定性通常高于0-100cm层。这可能是因为地表SOC更容易受到短期环境扰动的影响，如土地利用变化、植被动态和气候极端事件。尽管存在这些区域热点，中国大部分地区的预测不确定性仍低于10%，这是一个广泛认为可以接受的阈值（Sun等人，2012；Vaysse和Lagacherie，2017）。这表明我们基于QRF的不确定性估计足够可靠，可以支持农业管理和碳政策制定中的明智决策。

3.2. 历史SOC映射和未来预测
为了验证我们建模框架的稳健性，我们量化了两个历史时期（1980-1999年和2000-2015年）0-100 cm范围内SOC储量的变化幅度和空间模式（图4a，b）。结果显示了时期间的变化很小，分别为88.17（1980-1999年）和88.72（2000-2015年）Pg C。这些估计与独立的历史平均SOC储量估计值非常接近，例如89.03 Pg C（Yang等人，2023）和84.62 Pg C（Zhang等人，2023）。从空间上看，SOC的分布在这两个时期之间基本保持一致。青藏高原东部和内蒙古东北部的SOC较高区域保持稳定，而通常碳含量较低的内蒙古地区则显示出可测量的增加，这与自1990年代末以来为对抗荒漠化和土壤退化而实施的大规模生态恢复计划的影响一致。在我们的0-20 cm估计中观察到了类似的空间模式（图S3），这进一步证实了我们的建模方法能够可靠地捕捉中国SOC动态的绝对幅度和地理异质性。

展望未来，图4c-f展示了1980-1999年基线时期与21世纪末（2081-2100年）之间0-100 cm SOC储量的预测变化（?SOC），这些变化是在四种共享社会经济路径（SSPs）下平均得出的。图4c-f中的SSP126、SSP245、SSP370和SSP585分别表示。插图（c-f）展示了来自四个单独地球系统模型（ACCESS-ESM1-5 [红色]、EC-Earth3-Veg [蓝色]、MPI-ESM1-2-LR [紫色]、IPSL-CM6A-LR [黄色]）及其多模型平均值（黑色）的预测变化的核密度估计（KDE）分布曲线。这些结果表明，在所有情景下，尤其是高排放情景下，青藏高原南部出现了显著的SOC积累（高达0.6 g C/cm2），这与新的研究发现一致（Jiang等人，2025；Zhang等人，2024d）。这种增长很可能反映了由气候驱动的植被扩张，这种扩张得益于全球变暖引发的冰川退缩、生长季节延长以及净初级生产力的提高。相比之下，中国东北部高SOC（土壤有机碳）地区在所有排放情景下都表现出SOC的持续损失。这些损失可能是由于温度上升导致的分解加速造成的。尽管不同模型的预测结果在极端变化的幅度和位置上存在一些不确定性（见图4c-f和图S7-13），但它们揭示了一个明显的空间极值：寒冷高海拔地区的碳储量增加（尤其是在青藏高原），而中国东北部的温带耕地和森林则出现了碳储量减少。模型间相对较小的不确定性范围表明，对于这些相互对立的区域趋势存在广泛的共识。这增加了我们基于多模型预测的可靠性，并强调了其在指导中国生态多样性地区的气候适应和土壤碳管理策略方面的实用性。

3.3. 环境因素对SOC动态的影响
图5显示了四个SSP（共享社会路径）下0-100厘米深度SOC储量的预期趋势，以及一系列敏感性实验，这些实验依次隔离了特定驱动因素组（气候变量、土地利用变化和碳通量）的影响。0-20厘米层的结果呈现出类似的模式（见图S14）。排除ACCESS-ESM1-5在内的多模型集合显示，在所有SSP下SOC储量都呈一致增加趋势，尤其是在SSP370情景下预计到2081-2100年将超过95Pg C，但在最高排放情景SSP585下则不是这样。这种非单调响应表明，较高的排放量并不一定导致更大的SOC损失或增加，而是反映了气候强迫、植被生产力、凋落物输入和SOC稳定过程之间的复杂相互作用。

值得注意的是，ACCESS-ESM1-5模型与其他模型存在差异，它在所有情景下都预测SOC会减少；这一异常值突显了模型在表示生物地球化学反馈方面的特殊性。这些差异（空间分布见图S6-12）可能与该模型对植被和GPP（全球初级生产力）的系统低估有关，这导致了未来SOC的减少（见图S14-17；Gier等人，2024年）。到2080-2100年期间，ACCESS-ESM1-5模型在中国西北部大部分地区的预测结果与其他模型存在较大偏差（见图S18），这主要是由于其对GPP的低估。尽管如此，其余模型之间的一致性支持了一个核心结论：适度的排放路径（例如SSP245）可能为中国陆地生态系统创造了最佳的净SOC封存条件，这可能是由于变暖增强了生产力并且分解速率在可控范围内。

为了区分关键驱动因素的相对作用，我们进行了因子敏感性模拟，其中一组驱动因素被固定在2000-2015年的历史平均值上，而其他因素则根据每个SSP情景进行演变（见图5，第2-4行）。碳通量被证明是主要控制因素，多模型集合预测SOC储量将呈负增长，从SSP126到SSP585情景下到2100年减少了约20 Pg C（空间分布见图S7和S8）。这种日益不利的效应反映了未来碳输入的减少，而这种减少本可以在动态植被响应下发生（例如CO?施肥、变暖增强的生产力）。随着时间的推移，这种赤字的幅度会放大，并且在较高排放情景下变得更加明显，尽管模型间的差异仍然存在。相比之下，气候变量（温度、降水量等）的影响较为温和（<+7 Pg C）。这表明，允许气候根据SSP情景演变通常会抑制SOC的积累，这可能是由于变暖加速了分解过程。

土地利用变化在SSP126、SSP245和SSP370下始终产生负面的SOC影响：当土地覆盖保持在2015年的条件下时，模拟的SOC储量低于完全动态模拟的结果，这意味着这些情景下预期的未来土地利用变化将导致SOC的净增加。SSP126和SSP245下的净SOC增加与森林面积的增长相符（见Hou等人，2022年）。对于SSP370，当森林面积减少时，SOC意外增加（见图S19中的次级平均生物量密度示例）。这可能与SSP370下耕地面积的扩大有关（见Hou等人，2022年）。C3植被在预测0-100厘米深度SOC方面起着重要作用（见图S20），耕地面积的增加增强了C3植被的影响，从而抵消了森林损失导致的SOC减少。然而，在SSP585下，这种效应在2041-2060年后发生了逆转，这可能反映了大规模的造林或高排放情景下的牧场恢复。

总体而言，这些因子实验证实了生态系统碳通量的变化是长期SOC轨迹的主要决定因素。随着本世纪的发展，由于碳输入受到抑制而产生的赤字逐渐超过了气候和土地覆盖变化的较温和影响。这种层次结构强调了未来SOC储存的能力在很大程度上取决于陆地碳循环的完整性，特别是生态系统在全球变化下维持或增强碳输入的能力。

3.4. GPP与SOC联系的相对强度
我们提供了GPP的历史平均值和空间明确的相对重要性，以量化中国境内GPP与SOC储存之间的功能联系强度（见图6）。通过划分这些基于SHAP（Shapley Value）的贡献，我们确定了GPP-SOC耦合的主要区域和幅度，阐明了植物碳通量如何调节SOC的变化。

变量重要性分析表明，GPP是SOC分布的主要控制因素。对于100厘米深度的土壤剖面，GPP解释了模型可解释性的13.5%，仅次于最高温度（Tmax_04，见图6a）。在表土层（0-20厘米）中，GPP占据了9.1%的预测权重（见图6d）。这种高联系强度证实了植物碳通量与全国范围内SOC的长期稳定功能上的紧密耦合，是预测框架的基石。

这种联系的空间分布显示了中国西北部干旱和半干旱地区的GPP影响最为显著，那里SOC的变化对植被生长的波动最为敏感。绝对SHAP地图（见图6c和S8）进一步通过局部碳核算的角度描述了这种影响。我们发现中国西北部存在明显的碳输入不足（表现为负SHAP值），那里的生产力水平显著低于全国平均水平。在这种情况下，我们发现了一个生物物理瓶颈，即生产力不足成为将当地SOC储量拉低到大陆基线的根本限制因素。相反，中国南部和东北部的湿润地区表现出碳输入过剩（正SHAP值），那里丰富的生产力积极推动了SOC的积累。这种空间分解表明，负的集合平均值反映了高敏感性区域的资源限制，强化了GPP作为正向碳封存驱动因素的作用，其强度受到水文气候因素的地理限制。

3.5. GPP与SOC关系的非线性
图7展示了SOC对夏季GPP（GPP_S2）的非线性依赖性，这是通过RF模型得出的条件SHAP值来量化的。在所有SSP情景和历史时期，GPP的SHAP值都表现出一个非常一致的阈值响应：在低GPP水平下贡献为负值，当达到约4 g C m-2 d-1的临界阈值时转变为正值，然后在超过约7 g C m-2 d-1时趋于稳定。这种符号反转和这些拐点的稳定性在多个模型和不同土壤深度上都是可靠的（见图S20-23），证实了内在的GPP-SOC联系是生态系统的一个结构特性，对外部气候强迫具有不变性。这种符号反转在多个模型（见图S21）和表土层（0-20厘米；见图7b）中都是稳定的，证实了非线性不是模型结构或深度整合的产物。

3.6. GPP-SOC关系的整体趋势和空间转变
GPP-SOC联系的时间演变进一步强调了未来碳增加的外在性质。如图8a和8c所示，无论SSP轨迹如何，SHAP份额（%）——GPP相对联系强度的代理指标——在历史和未来时期（1980-2100年）几乎保持不变。这种稳定性表明，模型对GPP的内部敏感性随时间得到了保持。因此，绝对SHAP值从负值稳步上升至正值（见图8b和d）完全是由GPP幅度的上升趋势驱动的。

图9展示了这种动态的空间异质性。在中国大部分地区，GPP对SOC的贡献在整个时期要么持续为负，要么为正。有趣的是，出现了一个明显的过渡区，在这里GPP-SOC关系从负转为正（以红色表示）。关键的是，这种符号反转区域的范围随着排放量的增加而扩大。此外，这种转变在高排放情景下的发生时间逐渐推迟。在SSP585下，表现出从负转正的像素比例显著高于SSP126（见图9中的插入直方图），特别是在中国北部的生态过渡带。这些过渡区域集中在400毫米等降水量线沿线，这是中国湿润/农业区与半干旱/牧区之间的生态边界。这突显了半干旱生态带作为未来SOC封存的关键前沿，在这里由变暖引起的生产力增加最终可能克服历史上的分解主导地位，将边缘土地转变为新兴的碳汇。

总之，这些结果揭示了一个稳健的非线性GPP-SOC联系，其中净碳贡献受到相对于定义的生物物理阈值的生产力绝对大小的控制。

3.7. GPP-SOC关系的总体趋势和空间转变
GPP-SOC联系的时间演变进一步强调了未来碳增加的外在性质。如图8a和8c所示，无论SSP轨迹如何，SHAP份额（%）——GPP相对联系强度的代理指标——在历史和未来时期（1980-2100年）几乎保持不变。这种稳定性表明，模型对GPP的内部敏感性随着时间的推移而保持。因此，绝对SHAP值从负值稳步上升到正值（见图8b和d）完全是由GPP幅度的上升趋势驱动的。蓝色（PtN）表示从正方向变为负方向，而红色（NtP）表示从负方向变为正方向，周期则表明了变化发生的时间。插图中的直方图（计数百分比）量化了这些转变的相对频率。

4. 讨论
本研究探讨了在全球气候变化背景下，全球初级生产力（GPP）与土壤有机碳（SOC）之间的联系如何塑造中国的未来SOC动态。我们首先利用机器学习（ML）框架，基于历史观测数据构建了一个随机森林（RF）模型，以在历史条件和四种不同的气候变化路径（SSPs）下预测SOC的变化（模型性能：R2 = 0.41）。接着进行了因子敏感性分析，系统地分离了不同环境因素的贡献，以揭示全球变化影响SOC的途径。最后，通过应用模型不可知的解释工具（如SHAP和RFI），我们展示了GPP-SOC关系的强度、非线性以及时空转变特征。分析表明，碳通量是控制SOC变化的主要因素，其影响超过了气候和土地利用变化的直接影响。尽管最高温度的重要性稍高一些，但GPP-SOC关系在各种模型和情景中都表现出极高的稳健性，在驱动因素层级中始终排在第三位。关键的是，这种关系具有非线性特征，存在一个约为4克碳每平方米每天（g C m-2 d-1）的阈值。这导致在年降水量为400毫米的等雨量线上形成了一个过渡区，在该区域GPP的作用下，局部碳预算可以从负值变为正值。这些结果突显了植被驱动的碳封存能力在调节气候变化对SOC响应方面的潜力。对于中国多样的生物群落，尤其是其脆弱的半干旱生态过渡带，管理景观以跨越关键的生产力阈值，可以将边缘土地转变为持久的碳汇，为国家的碳中和策略提供可行的路径。

我们研究结果的稳健性得到了多种证据和分析框架的支持。首先，SOC预测RF模型使用涵盖两个关键时期（1979-1985年和2000-2014年）的全国性土壤数据集进行了训练和验证（图1），并在SOC预测方面取得了可接受的性能（R2 = 0.41）。其次，重建的历史SOC储量与近期研究的结果一致（Li等人，2022b；Liang等人，2019b；Liu等人，2022；Song等人，2020a；Yang等人，2023a；Zhang等人，2023a；表S2）。同时，我们的模型还预测了在所有气候变化路径下全国SOC的净增加（图S4和S5），这与之前的中国研究结果一致（Wu等人，2024a）。此外，我们还再现了中国东北地区的SOC减少现象（Wang等人，2023a）。我们的模型能够协调这些模式，说明它有效地整合了来自多源预测因子的信息，从而能够反映复杂的SOC响应。第三，不确定性量化显示了稳定的误差范围（相对不确定性约为10%，图3），进一步证实了该模型的稳健性优于其他方法（Kakhani等人，2024；Ogle等人，2010）。第四，我们采用的分组变量分解技术（Wu等人，2024b；Zhang等人，2024a）能够清晰地区分SOC对不同驱动因素的响应，从而分离出植物碳通量与直接气候影响的贡献。最后，SHAP值的应用为GPP-SOC关系提供了可靠的量化（图S20-23），相比传统的基于相关性的方法具有两个主要优势：（i）对预测因子多重共线性的鲁棒性（Huang等人，2023）；（ii）与基于排列的RFI结果的高度一致性（图S25）。这些交叉验证共同证实了我们方法论的可靠性和主要发现的科学可信度。

GPP-SOC关系代表了气候变化下碳循环中一个关键但研究不足的维度。尽管大多数现有研究，尤其是那些使用ML模型的研究，主要关注温度升高和降水模式变化如何通过影响微生物分解速率来影响SOC（García-Palacios等人，2021；Wang等人，2023b），但往往忽略了另一个重要途径：二氧化碳（CO?）施肥作用导致的植物生产力增加以及随后土壤中碳的积累（Wang等人，2025）。我们的分析表明，碳通量对SOC动态的控制作用比直接的气候变量更强且更持久（图5）。特别是，GPP作为CO?施肥效应的替代指标比大气中的CO?浓度本身更为有效（Chen等人，2022）。当我们明确将CO?浓度作为协变量纳入RF模型时，模型未能捕捉到预期的施肥效应。这种差异源于CO?浓度的空间和季节性变化相对于未来变化的幅度较小，因此基于时间替代的空间方法无效。相比之下，我们提出的利用碳通量作为协变量的方法对于在ML模型中模拟CO?施肥效应非常有用。GPP直接反映了植被的碳吸收能力及其对气候和升高CO?水平的响应（Liu等人，2025；Manzoni等人，2018），使其成为潜在SOC输入的机制性指标。变量重要性分析也显示，GPP的重要性超过其他所有变量，只有温度除外（图6）。这些发现强调了植被介导的碳输入可能在一定程度上抵消由升温引起的SOC损失。

我们的分析揭示了一个稳健的、非线性的GPP–SOC关系，这种关系在不同的时间尺度、地球系统模型（ESMs）配置和气候变化路径（SSPs）中都普遍存在（图7；图S20–S23），这一一致性强烈表明所识别的阈值反映了土壤-植被系统的固有生物物理特性，而非模型偏差。在未来的气候变化驱动下，GPP的时空重新分配通过三种不同的机制从根本上重塑了全国的SOC平衡。在低生产力区域（例如中国西北部的温带草原），由于基准GPP通常低于4克碳每平方米每天（g C m-2 d-1）且土壤碳氮比（C/N）较低，即使GPP略有增加也无法增加SOC。这一反直觉的结果源于微生物的启动效应：易分解的碳输入（如根系分泌物）会刺激微生物代谢，加速已有稳定SOC库的分解（Fontaine等人，2003；Ding等人，2025）。这种机制对水热变化的敏感性在低GPP下的不同气候变化路径响应中得到了体现（图7a-b），进一步证实了其对升温引起的碳损失的脆弱性，这与在高温下碳封存作用增强的培养实验结果一致（Schiedung等人，2023；von Fromm等人，2024）。这种碳赤字机制与传统植被碳封存范式形成鲜明对比（Tian等人，2015）。然而，其稳健性得到了三方面的证据支持：（i）经过偏差校正后的12个CMIP6 ESMs结果的一致性（图S21）；（ii）来自干燥生态系统实地观察的结果，即在水分供应不足的情况下植树造林会导致SOC下降（Hoffmann等人，2017；Kunkel等人，2022）；（iii）历史重建的可重复性，其中GPP持续低于约4克碳每平方米每天的区域显示出亚平衡的SOC储量（图7a，b）。

当GPP超过这个阈值时，即年降水量达到400毫米的等雨量线（图9），系统转变为碳盈余状态。此时，增强的光合作用输入（通过凋落物、根系周转和分泌物）超过了分解速度，使得SOC在物理和生化上趋于稳定（例如矿物结合、团聚体封闭；Schrumpf等人，2013；Kramer和Chadwick，2018）。在高排放情景（SSP5-8.5）下，CO?施肥作用推动这一正向耦合区向西北方向扩展，将边际土地转化为新的碳汇（图9d）。关键的是，这种转变取决于水分的可用性：只有在水分同时支持持续的生产力和有机矿物形成时，GPP才能转化为持久的SOC增益（Védère等人，2022）。在高生产力区域（中国南部/东部），GPP–SOC反应在约7克碳每平方米每天（7 g C m-2 d-1）附近达到饱和，表明土壤碳饱和度达到极限（Shi等人，2021）。虽然持续的高碳负荷可能会重新激活微生物启动效应（Sitch等人，2015），但这些生态系统仍保持净碳盈余，证实了它们作为稳定碳库的作用。

这些结果颠覆了“更多植物=更多土壤碳”的简单观点。相反，它们揭示了一个由生产力阈值控制的三分法框架（赤字、盈余和饱和），这些阈值受水分可用性的调节。模型、情景、历史数据和独立文献之间的一致性为这些非直观的动态提供了强有力的验证，为预测和管理变化气候下的未来SOC轨迹提供了精细的机制基础。在我们的预测中，青藏高原（TP）在未来的气候情景下会出现广泛的SOC积累，这一结果与许多升温模拟实验和早期模型预测的SOC下降相矛盾。为了解释这一差异，我们首先考察了历史基准数据（图4b）：在1980-2015年间，基于观测数据的模型表明，在TP的干旱西部地区，SOC储量基本保持稳定或略有下降，这与寒冷且水分受限条件下碳输入有限的实际观察结果一致。然而，CMIP6多模型集合一致预测TP在所有排放情景下会经历“升温-湿润”的趋势（图4c-f），不仅温度上升，降水量和生长季节长度也显著增加（Dang等人，2014）。这种水文气候变化从根本上改变了碳平衡：虽然升温确实加速了异养呼吸作用，但也解除了对植被生产力的关键限制，从而显著增加了GPP和凋落物碳输入——尤其是在东南部和中部TP，那里水分增加最为明显。我们的框架捕捉到了这种非线性响应：一旦GPP超过当地阈值（高山草甸约为3-4克碳每平方米每天），新的碳输入开始超过分解损失，驱动净SOC积累。这一机制得到了基于过程的最新建模研究的支持（Jiang等人，2025；Zhang等人，2024b），这些研究也将预测的TP碳汇归因于CO?施肥效应与增强的水分可用性的共同作用。因此，与早期实验研究的结果相比，我们的发现强调了在考虑实际未来气候情景时需同时考虑碳-水耦合动态的必要性。我们的结果表明，TP作为碳源或碳汇的命运不仅取决于温度，还取决于降水是否足够增加以激活生产力的响应，从而抵消呼吸作用造成的损失。

尽管我们的研究为中国未来的SOC动态提供了重要见解，但在解释结果时也应认识到几个局限性。首先，土壤样本的空间分布不均，特别是在干旱的西北地区和高海拔地区（图1），这在模型构建中引入了不确定性，并可能使数据贫乏区域的SOC估计产生偏差。其次，我们的预测结果包含了来自ESMs的固有不确定性。ESMs在模拟碳通量、气候极端事件和陆地-大气反馈方面的差异会传递到我们的ML框架中，反映了表示复杂生物地球化学过程的广泛挑战。例如，ACCESS-ESM-1-5模型显示GPP值系统性地偏低（图S14-17），特别是在SSP585中，导致未来SOC的低估（图5）。这突显了继续进行模型比较、过程级改进以及整合多源观测数据的必要性。第三，数据驱动方法的一个根本局限性是它们依赖于历史数据的静态性，这种假设认为SOC与其驱动因素之间的关系随时间保持不变。特别是LUH2中的土地利用数据可能存在固有偏差，从而影响建模结果。然而，在快速气候变化的情况下，新的环境条件（如前所未有的温度-湿度组合、CO?水平或干扰机制）可能会以过去数据无法捕捉的方式改变生态系统功能（Pickett，1989）。因此，我们的模型在超出历史协变量范围时的外推能力，尤其是在高排放情景下，仍然受到限制。最后，尽管SHAP和基于排列的可解释性方法为我们提供了关于驱动因素重要性和非线性的宝贵见解，但它们并非不受方法学偏差的影响，例如对特征相关性的敏感性和采样偏差。在从观测模型推断因果机制时，这些因素都应予以考虑。

总体而言，这项研究提供了基于观测的证据，表明植被驱动的碳封存可以显著抵消中国各地由温度升高引起的SOC损失。我们确定了年降水量达到400毫米的等雨量线作为控制国家尺度SOC动态的关键生物地理阈值，为优化中国的重大生态项目提供了直接的科学指导，包括长期植树造林和 Grain-for-Green 计划（Song等人，2014）、水资源再分配（Lv等人，2019）、沙尘暴缓解（Mahmoodi，2025）和生态恢复（Ding等人，2022）等。这一划分标识了从碳赤字向碳盈余制度的转变，为实现中国的碳中和目标，实现了从简单的绿化向以碳封存效率为导向的土地管理的战略转型。具体的空间分析突出了两个关键方面：（1）在年降水量400毫米线附近的过渡区，由于总初级生产力（GPP）与土壤有机碳（SOC）的关系从碳输入赤字转变为盈余，应该优先进行有针对性的生态恢复，以最大化SOC的增加；（2）当GPP超过某个临界阈值（约7克碳/平方米·天）时，SOC的增加可能会达到饱和状态，这意味着在年降水量已经超过800毫米的地区（例如中国南部），植树造林对土壤碳封存的贡献可能会逐渐减少。因此，我们的研究结果为在中国的主要生态工程项目中重新分配资源提供了定量依据，将资源从生态饱和区转移到这些对气候敏感的过渡前沿地区。此外，我们对GPP与SOC之间关系的量化分析强调了改进地球系统模型（ESM）中植被-土壤碳反馈表示的必要性。然而，这种向西北方向的扩展带来了水分与碳之间的权衡。虽然二氧化碳的施肥作用在干旱地区创造了新的碳封存机会，但追求最高的GPP可能会超出当地的生态水文承载能力。为了确保可持续性，恢复策略必须优先考虑水分与碳之间的协同作用，而不是仅仅关注生物量积累，以防止有限水资源的枯竭同时稳定土壤碳含量。将这些特定于生物群的阈值纳入土地利用规划对于将边际土地转化为有效的碳汇至关重要，同时避免在已经饱和的生态系统中出现收益递减的情况，从而直接支持中国的碳中和目标。

**5. 结论**
我们基于数据驱动的机器学习方法揭示了气候变化下碳通量是控制SOC动态的关键因素，其中GPP-SOC关系成为仅次于温度的第二大驱动因素。通过对不同气候情景（SSPs）的全面分析，我们量化了GPP-SOC关系的时空演变，并发现了两个关键特征：（1）当GPP约为4克碳/平方米·天时，与历史基线相比，GPP从碳输入赤字转变为盈余的非线性阈值；（2）存在一个与年中国400毫米等降水量线对齐的敏感气候过渡区，在高排放情景下，这一转变趋势逐渐向西北方向移动。这些发现对气候智能型土地管理具有重要意义。首先，应将植被介导的碳输入纳入SOC预测中，以考虑其抵消气候变暖引起的SOC损失的能力。其次，植树造林和生态恢复政策应基于碳循环反馈来制定，优先考虑那些GPP提升能够可靠超过所确定的生产力阈值的地区，从而引发积极的SOC响应。第三，400毫米降水量线是一个关键的干预前沿，通过有针对性的植被恢复，可以将GPP-SOC关系从碳源转变为持久的碳汇。通过将植被生产力作为SOC管理的着力点，我们的工作为地球系统建模提供了可靠的框架，并为基于自然的气候解决方案提供了可行的建议。未来的研究应集中在不同土壤类型、土地利用历史和管理实践下进一步完善这些生物地理阈值，以确保实施的可行性和针对性。

**作者贡献声明**
Feini Huang：撰写、审阅与编辑、验证、方法论、调查、数据分析、概念化。
Yongkun Zhang：撰写初稿、可视化、软件应用、资源管理、方法论、数据分析、数据整理。
Xingjie Lu：撰写初稿、验证、调查、概念化。
Wei Shangguan：撰写、审阅与编辑、项目监督、资金筹措、数据分析。
Qingliang Li：验证、方法论。
Zhangcai Qin：撰写初稿、验证、项目监督。
Zhongwang Wei：验证、项目监督、概念化。
Huan Yuan：撰写初稿、可视化、验证、项目监督、概念化。
Lu Li：撰写初稿、验证、数据分析、概念化。
Yongjiu Dai：项目监督、资金筹措。