中国农业管理方式改变了中国不同轮作制度下的表土有机碳动态

《Soil and Tillage Research》：Agricultural management reshaped topsoil organic carbon dynamics across diverse crop rotation systems in China

　　•对国家历史土壤调查数据集的重新分析。•对表土有机碳动态的时空非实验性评估。•在环境约束下农业管理对碳分布的影响。•过量施氮与轮作导致大量碳积累。•当氮肥施用量超过约15克/平方米·年时，其固碳效益会下降。

1. 引言
耕地中的土壤有机碳是调控全球碳循环与气候变化的关键因素，也是保障粮食安全的基石（Bossio等人，2020；Ma等人，2023；Wang等人，2023）。由于农业土壤有从历史上的碳源转变为重要碳汇的潜力（Ruehr等人，2023），它们为减缓气候变化与改善土壤健康提供了“双赢”的解决方案（Frank等人，2024）。然而，土壤有机碳的动态变化受人为干扰与气候、地形、土壤条件等自然因素的复杂相互作用影响。虽然自然因素决定了微生物生物量产生（Wang等人，2025）和分解过程（Davidson和Janssens，2006）的基准分布（Tan等人，2024），但农业集约化趋势日益削弱了这些自然控制因素的作用。这一现象在中国尤为明显，自20世纪80年代以来，中国农业集约化程度达到前所未有的水平（Jiao等人，2018；Zhao等人，2018），表现为大量化学肥料的使用（Li等人，2013）、秸秆还田政策的实施（Li等人，2018；Lin等人，2024；Qu等人，2012）、种植结构的调整（Chen等人，2025；Yu等人，2024）以及灌溉条件的改善（McDermid等人，2023；Zhang等人，2022）。通过研究长期管理措施如何影响土壤有机碳的空间分布，可以了解不同地区和轮作体系下农业活动与土壤有机碳之间的关联程度，以及这些关系在不同管理方式下的变化情况。这有助于区分农业活动与自然因素在土壤有机碳重新分配中的作用。此外，明确管理措施与土壤有机碳时间变化之间的关系，有助于估算耕地土壤的碳固存能力，并评估未来农业策略的有效性。

准确理解农业管理在调控土壤有机碳动态中的作用，对于构建精确的地球系统模型和制定政策至关重要。但目前的研究方法存在明显局限。长期的实地试验虽能提供微观层面的见解（Chen等人，2024；Gocke等人，2023），但缺乏空间代表性。元分析虽可从不同研究中提炼新发现（Lin等人，2023；Mo等人，2024；Ren等人，2024），但难以消除因实验差异及纳入/排除标准不同而带来的偏差。另一种方法是利用数字土壤制图技术（Yin等人，2024；Zhang等人，2024）、基于过程的生物地球化学模型（Attia等人，2024；Lin等人，2024），或结合这两种方法（Xie等人，2024），来推断有限的原始土壤数据，并识别环境因素与土壤有机碳之间的关联。不过，这类推断往往会使空间异质性趋于平滑，或无法捕捉农业活动时间演变与土壤有机碳变化之间的非线性、时变关系，比如土壤有机碳对气候及气候变化的不同响应。目前仍存在一个关键缺口，即难以明确判断土壤有机碳的时空变化在多大程度上是由农业活动引起的，也难以确定管理措施发挥有效作用或变得多余的临界点。

为解决这一问题，我们整合了两个历史时期的全国土壤剖面数据集（20世纪80年代和21世纪10年代），旨在重建中国耕地土壤有机碳的时空演变规律，具体目标包括：（1）描述不同轮作类型下的土壤有机碳时空变化；（2）区分管理措施（化肥施用、秸秆还田和灌溉）与自然因素对土壤有机碳变化的贡献；（3）为制定针对不同轮作体系的土壤管理策略提供新的依据。通过架起农业与生物地球化学之间的桥梁，本研究为依据区域种植系统发展合理农业技术提供了数据驱动的非实验性证据。

2. 材料与方法

2.1 数据收集

2.1.1 土壤数据
我们收集了中国20世纪80年代和21世纪10年代的两个历史土壤剖面数据集。80年代的数据集是根据《中国土壤志》（国家土壤调查办公室，1993）的记录数字化得到的，该志书基于1979年至1984年的第二次全国土壤调查。21世纪10年代的数据集则来自《中国土壤系列》（Zhang，2017），这是2009年至2019年全国土壤系列调查之后编纂的一系列书籍。这两个数据集中的观测值均为野外土壤剖面数据，而非实验数据，其采样是基于相应土壤类型和自然环境条件进行的。本研究仅考虑耕地区域，只保留那些明确记载了耕地用途和作物类型的土壤剖面数据。我们最终获得了包含作物类型、轮作类型以及土壤性质描述的有效耕地观测数据，这些土壤性质包括土壤有机碳含量（SOC，单位：克/千克）、黏土含量（百分比）和土壤pH值。部分观测数据记录的是土壤有机质含量（SOM）而非SOC（21世纪10年代的数据中有约30%如此，80年代的所有数据均为SOM）。为便于数据统一处理，考虑到获取动态转换因子的不便（Pribyl，2010），我们采用传统的固定转换系数1.724（Wolff，1864）将SOM转换为SOC（单位：克/千克）。在本研究中，黏土含量和土壤pH值被用作反映土壤条件的指标。

土壤性质是在土壤剖面的每一层中测量的。由于各层厚度不一致，我们采用线性变换方法将土壤性质标准化到表土层（0–20厘米深度），具体公式如下：（1）、（2）、（3）。在标准化之前，土壤pH值需先转换为H+离子浓度。

（1）SOC20=∑i=1nSOMdi×di+SOMdj×20?dj20×1.724
（2）Clay20=∑i=1nClaydi×di+Claydj×20?dj20
（3）pH20=?lg∑i=1n10?pHdi×di+10?pHdj×20?dj20

其中，SOC20、SOM20、Clay20和pH20分别表示0–20厘米深度处的标准化土壤性质。di为某一土壤层的厚度，若下一层深度大于20厘米，则dj为最上层土壤的深度。n为完全位于0–20厘米深度范围内的土壤层数量。

2.1.2 环境与管理指标
本研究采用年平均降水量（MAP）、年平均气温（MAT）和干旱指数（AI）作为气候指标（Peng，2024；Peng，2019a；Peng，2019b；Peng等人，2019）。这些指标的空间分辨率为1公里，数据基于1980年代（1971–1980年）和21世纪10年代（2001–2010年）采样前的10年平均值。海拔高度和复合地形指数（CTI）则用作地形指标。海拔数据来源于NASA Earthdata提供的Shuttle Radar Topography Mission（SRTM）全球数字高程模型（DEM，NASA，2013），该模型已预处理为250米分辨率。CTI数据则来自OpenTopography在2020年发布的Geomorpho90m地形产品，其分辨率为90米，可反映地形湿润程度。

在本研究中，我们以化肥氮输入量作为管理指标。该数据来源于Yu等人的研究（Yu等人，2022a；Yu等人，2022b），其中包括仅含氮肥和复合化学肥料带来的氮输入量（已转换为元素氮含量，单位：克/平方米·年），有机肥料未被计入。化肥氮输入量的空间分辨率为5公里，数据同样基于1980年代（1971–1980年）和21世纪10年代（2001–2010年）采样前的10年平均值。

秸秆还田可以将收获后的一部分碳重新返还到土壤中，包括机械直接还田以及通过焚烧将富含矿物质的秸秆灰分还田。我们根据以往的研究（Liu和Li，2017；Ran等人，2021；Zhao等人，2018）以及谷物产量数据，估算了直接还田带来的秸秆碳输入量。为此，我们使用了Cao等人在2025年和2024年开发的GlobalCropYield5min数据集，该数据集提供了水稻、小麦、玉米和大豆的5弧分钟分辨率的谷物产量信息，用于在采样前计算网格级别的秸秆碳输入量。对于每个网格单元，年度秸秆碳输入量是这四种作物输入量的总和。每种作物的秸秆碳输入量则通过公式（4）计算：

（4）StrawCinput=yield×(1?wc)×r×α×β×0.45

其中，yield为作物产量，wc为水分含量，r为不同地区的典型秸秆与谷物比例，α为可收集系数，表示实际可收集的秸秆资源量与理论资源量之比，β为直接还田率，这一数值会因地区和时间的不同而有所变化。0.45这一转换系数用于将作物生物量转换为SOC含量。秸秆碳输入量指标采用了两种不同的平均时间跨度：由于80年代的数据较为有限，采用5年平均值（1982–1986年），而21世纪00年代则采用10年平均值（2001–2010年）。各项参数的具体数值详见补充表1和表2。

我们还通过公式（5）计算了化肥氮输入量和秸秆碳输入量的平均时间梯度，以此反映这些指标在不同时间段内的变化幅度：

（5）G?=1nx1?x0+∑t=1n?2xt+1?xt?12+xn?1?xn?2

其中，G?为某指标的平均时间梯度，n为年份数，x0到xn?1分别为各年份的指标数值。在本研究中，管理指标的绝对水平指的是采样前该指标的平均值，而时间梯度则代表该时期指标的平均变化速率，因此能够反映管理活动的动态变化，而非单纯的输入量大小。

灌溉也是本研究中的一个管理指标。每个网格单元的灌溉面积占比数据来自CIrrMap250数据集（Zhang等人，2024a；Zhang等人，2024b），该数据集包含了中国耕地灌溉情况的年度地图。由于80年代缺乏空间上明确的灌溉数据，我们使用2000年的灌溉地图来代表80年代的灌溉状况，2010年的地图则用于表示21世纪10年代的灌溉情况。所有地图都通过双线性插值处理为1公里分辨率，以减少在将灌溉数值分配到各个数据点时因空间偏移而产生的误差。

2.2 数据处理

2.2.1 数据清洗与分组
我们总结了中国的10种主要轮作体系（见表1）。只有那些明确记录了作物类型且并非仅用于种植蔬菜的观测数据才会被保留。任何包含一个或多个指标缺失值的观测数据都会被剔除。经过筛选后，80年代共有1234个观测数据被保留，21世纪10年代则有1763个观测数据被保留。土壤观测数据的分布情况如图1所示。对于每一个观测数据，我们还增加了9个二元变量来表示相应的轮作体系，“1”表示该观测数据属于该轮作体系，其他轮作体系对应的变量则为“0”。

表1 本研究中总结的10种主要轮作体系的特征。

| 轮作体系 | 80年代观测数据数量 | 21世纪10年代观测数据数量 | 特征描述 |
|------------------|--------------------|------------------------|--------------------------------------------------------------------------|
| 单一旱地（西北） | 1782 | 2237 | 以种植谷物为主的单一旱地，主要分布在中国的西北部。 |
| 单一旱地（东北） | 1802 | 2376 | 以种植谷物为主的单一旱地，主要分布在中国的东北部。 |
| 旱地轮作 | 2573 | 3971 | 如小麦-玉米或小麦-大豆之类的轮作旱地，主要分布在中国华北平原。 |
| 单一水稻 | 6113 | 1613 | 以种植单季水稻为主的稻田，广泛分布于中国南部，北部也有少量分布。 |
| 单一水稻-旱地（西部）| 1007 | 7710 | 如水稻-小麦、水稻-玉米或水稻-油菜之类的轮作稻田，主要分布在中国西南部。 |
| 单一水稻-旱地（东部）| 7125 | 536 | 如水稻-小麦、水稻-玉米或水稻-油菜之类的轮作稻田，主要分布在中国长江中下游地区。 |
| 双季水稻轮作 | 2162 | 5664 | 用于种植双季水稻或与其他谷物作物及油菜轮作的稻田，广泛分布于中国南部。 |
| 谷物/根茎类作物 | 1301 | 1498 | 用于种植谷物、根茎类作物或二者轮作的旱地，主要分布在中国西南部。 |
| 糖蔗/根茎类作物 | 2737 | 673 | 用于种植根茎类作物、甘蔗或二者轮作的旱地，主要分布在中热带和热带地区。 |
| 高原地区 | 1111 | 12 | 主要分布在中国青藏高原地区，以种植青稞、玉米和油菜为主。 |
| 总计 | 1234（含3个未明确类别） | 1763（含6个未明确类别） | |

下载：下载高分辨率图像（754KB）
下载：下载全尺寸图像

图1 中国主要轮作体系下土壤观测数据的空间分布：（a）20世纪80年代的数据，浅绿色点，数量为1234；（b）21世纪10年代的数据，X标记，数量为1763；（c–k）不同轮作体系下的虚拟土壤样本对，深绿色点，数量为411。边界线为中国农业区划划分。

2.2.2 生成虚拟土壤样本对以模拟重复采样
由于缺乏重复采样数据，为模拟土壤有机碳的变化情况，我们采用了基于地理相似性原则的样本配对方法（Cui等人，2025；McBratney等人，2003；Zhu等人，2015；Zhu等人，2020；Zhu和Turner，2022）。根据这一原则，地理环境越相似，目标变量的数值也应越接近。对于土壤数据集而言，该方法的基本思路是根据一系列环境协变量的相似性，从其他时间点的真实土壤数据集中筛选出环境特征相似的样本。经过滤波后的样本的土壤属性会根据其与真实样本的整体相似度进行加权，从而生成伪土壤属性。我们选择了8个环境协变量，包括海拔、CTI、坡度、年平均温度、年平均降水量、干旱指数、经度和纬度，用以计算地理相似度。坡度数据则来自Geomorpho90m产品。具体流程如下：首先，所有协变量栅格数据都被重新投影到相同的坐标系统，然后重采样为5公里分辨率。每个网格单元包含2010年代双线性插值得到的协变量值向量。1980年代数据集中的每个数据点都会与其对应的网格单元进行空间关联，代表2010年代环境条件的单元值向量被作为伪协变量值赋予该点。其次，我们使用公式(6)、(7)、(8)、(9)来计算2010年代数据集中的点与协变量网格单元之间的环境相似度。(6)X=x1,x2,…,xm(7)ESi,j=1m∑k=1mexp?xk,i?xk,j22SDxk×SDxkRMSDxk,j2(8)SDxk=1n∑i=1nxk,i?μk2(9)RMSDxk,j=1n∑i=1nxk,i?xk,j2其中X是环境协变量向量，m为协变量的数量。ESi,j表示点i与点j之间的环境相似度。xk,i和xk,j分别是点i和点j的第k个协变量值。SDxk是所有网格单元中第k个协变量的标准差，μk是其均值。RMSDxk,j是所有数据点与点j之间第k个协变量的均方根偏差。n为网格单元的数量。第三，我们提取了与1980年代数据集中的点进行空间关联的网格单元的相似度向量。基于1980年代的伪协变量值和2010年代的实际协变量值，构建了两个土壤数据集的成对相似度矩阵。最后，那些环境相似度超过0.90的点被用来根据公式(10)计算1980年代数据集的伪土壤属性。我们选择0.90作为保守的阈值，以在保证环境可比性的同时尽量保留更多样本，不过也要认识到不同的阈值可能会影响配对结果和模型表现。(10)Pj=∑i=1n′ESi,j×Pi∑i=1n′ESi,j其中Pj是点j的伪土壤属性值。n′是与环境相似度超过0.90的点j相关的点数。ESi,j是所选点i与点j之间的相似度。Pi是点i的实际土壤属性值。为避免异常值的影响，只有那些所用点数在3到20之间的伪样本被选中用于计算伪土壤属性。经过筛选后，共有411对伪土壤样本被生成用于进一步分析（见图1）。这些配对代表的是推断出的匹配结果，并非同一地点的真实重复观测数据。伪有机碳含量、伪黏土含量以及伪土壤pH值被用来计算与1980年代实际观测到的土壤属性相比的时间变化情况。需要说明的是，环境和管理指标的时间变化是根据原始分辨率的栅格数据计算的，而非用于计算环境相似度的重采样栅格网格。生成伪土壤样本对的流程如图2所示。下载：下载高分辨率图像（397KB）下载：下载全尺寸图像图2. 基于环境相似度生成伪土壤样本对并计算伪土壤属性的流程。2.3. 机器学习建模及解析影响有机碳含量和有机碳变化的各因素我们运用随机森林模型来确定气候、地形、土壤环境以及农业管理指标在1980年代和2010年代有机碳含量回归分析以及有机碳变化中的相对贡献。随机森林具有对多重共线性不敏感且不易过拟合的优点，非常适合用于处理环境数据。为更好地评估模型稳定性，我们采用了5折交叉验证方法。随机森林模型在训练集上得到训练，在验证集上进行评估。对于1980年代和2010年代的数据集，交叉验证共进行了10次，以此根据验证集的预测均值生成评估指标。随机森林回归分析是使用Python中的scikit-learn包实现的。我们还运用SHapley加性解释方法来解读随机森林回归的结果。SHAP的核心理念源自于联盟博弈论中的Shapley值。考虑到所有可能的玩家组合，某个玩家的Shapley值是指其在所有可能组合中对预测结果的平均边际贡献。因此，Shapley值并非简单地等于从模型中移除某个特征后预测值的差异，因为那样做会忽略特征之间的交互效应。SHAP引入了新的方法来估算用于解读机器学习模型的Shapley值。SHAP以数据集的预测均值作为基准，单个数据点的最终预测值等于基准值加上所有特征的正面和负面Shapley值之和。此外，整个数据集中所有特征的绝对Shapley值平均值能够反映某个特征的整体重要性。SHAP分析也是使用Python中的shap包来实现的。基于Shapley交互指数概念，该概念不仅考虑单个玩家，还考虑所有玩家组合之间的贡献分配，SHAP能够生成Shapley交互值，从而在随机森林中区分交互效应与普通Shapley值。我们对所有特征都提取了成对的Shapley交互值。在交互作用分析中，重点关注了两个管理变量：化肥氮输入量和秸秆碳输入量。对于每一组管理措施与轮作方式的组合，我们将属于目标轮作方式（轮作变量值为1）的样本与其余所有轮作方式（轮作变量值为0）的样本进行比较。我们通过为目标轮作方式以及其余轮作方式绘制估计的LOWESS曲线，来总结二者之间Shapley交互值的变化差异。在每个十年里，我们还会使用Welch双样本t检验来比较目标轮作方式与其余轮作方式之间的Shapley交互值。3. 结果3.1. 1980年代、2010年代不同轮作方式下表土有机碳的特征以及有机碳变化1980年代和2010年代数据集的有机碳分布情况分别如图3a和3b所示。1980年代的有机碳含量在0.7到57.7克/千克之间，中位数为10.0克/千克。2010年代的有机碳含量则在0.3到49.6克/千克之间，中位数为12.0克/千克。各轮作方式子组的有机碳中位数以及25%至75%的百分位数分别展示在图3d和3e中。总体而言，2010年代各轮作方式下的有机碳中位数都要高于1980年代。我们对不同轮作方式子组之间进行了非参数的Dunn检验，并应用Bonferroni校正法，结果显示存在显著差异，表明不同轮作方式下的有机碳含量存在差别。旱地土壤的有机碳含量通常低于水田土壤。旱地轮作和单一旱地轮作（西北地区）的有机碳含量最低，而单一旱地轮作（东北地区）的有机碳含量最高。水田土壤的有机碳含量则更为均匀。双季水稻轮作和单季水稻与旱地轮作（西部）的有机碳含量最高，而单季水稻与旱地轮作（东部）的有机碳含量最低。下载：下载高分辨率图像（377KB）下载：下载全尺寸图像图3. 不同时期（a–c）以及不同轮作方式（d–f）下有机碳的含量及其变化情况。粗黑色的垂直线表示有机碳的中位数。图3d–f中的实心点代表各轮作方式的中位数有机碳含量，而水平线则表示四分位数范围（25%至75%的百分位数）。括号内标注的是轮作方式组别之间存在显著差异的情况（经过Bonferroni校正的Dunn检验，P值小于0.01）。为便于查看，图中仅展示了P值小于0.0001的显著差异结果（见图3d–e）。图3f中的星号表示基于双侧Student’s t检验，有机碳变化与0之间存在显著差异：*表示P值小于0.05；**表示P值小于0.01；***表示P值小于0.001。从1980年代的1234个观测数据中，共生成了411对伪土壤样本。这些伪样本对应的有机碳变化分布情况如图3c所示。有机碳变化范围在-15.68到19.59克/千克之间，中位数为2.77克/千克。各轮作方式子组的有机碳变化中位数以及25%至75%的百分位数展示在图3f中。总体来看，大多数轮作方式下的有机碳含量都呈现上升趋势，近70%的样本对显示出有机碳含量的增加。其中，旱地轮作方式下的有机碳含量增幅最为显著，在Student’s t检验中其平均变化量与0存在显著差异。值得注意的是，尽管这些地区的初始有机碳含量较低，类似情况也出现在中国西北地区的单季旱地轮作系统中，但依然出现了有机碳含量的上升。相比之下，单季水稻轮作系统以及非轮作系统的有机碳变化中位数均低于0，这表明在这些有机碳含量原本较高的系统中可能存在净有机碳损失的情况。而双季水稻轮作系统则显示出明显的有机碳含量上升趋势，反映出以水田为基础的种植系统之间存在差异。对于那些分布在气候较为干燥或半干燥的中国西北地区的单季谷物旱地系统，其有机碳含量也有了显著提升。相反，分布在气候较为湿润的中国东北地区的此类系统则没有出现明显的有机碳含量增加。3.2. 有机碳与各类指标之间的关系为了构建1980年代和2010年代数据集的随机森林模型，我们共使用了12个连续型指标和9个二元轮作变量。随机森林模型的详细评估指标详见补充表3和表4。图中用平均绝对Shapley值来表示各指标在随机森林模型中的重要性，如图4所示。所有连续型指标对应的Shapley值变化情况则展示在图5中。对于1980年代的数据集，干旱指数、黏土含量以及土壤pH值是最为重要的指标。从SHAP总结图（图4a）可以看出，较高的干旱指数和黏土含量值对应着正的Shapley值，这意味着它们会对预测的有机碳含量产生正向影响。相反，土壤pH值则会对预测的有机碳含量产生负向影响。旱地轮作是其中最重要的农业相关变量，它通常会降低预测的有机碳含量。而双季或单季水稻轮作变量则往往会使预测的有机碳含量上升。化肥氮输入量及其随时间的变化趋势也对有机碳含量的预测有着重要影响，其重要性程度与年平均温度、年平均降水量以及复合地形指数相当。较高的化肥氮输入量会对预测的有机碳含量产生负向影响，其转折点大约在2克氮/平方米·年左右（见图5ai）。同样，化肥氮输入量的变化趋势也会使预测的有机碳含量下降，其转折点大约在0.3克氮/平方米·年左右（见图5aj）。与化肥氮输入量相比，秸秆碳输入量及其随时间的变化趋势对有机碳含量的影响要小一些。秸秆碳输入量与有机碳含量呈负相关关系。而秸秆碳输入量的变化趋势则与有机碳含量呈正相关关系，其转折点大约在0吨碳/公顷·年左右。下载：下载高分辨率图像（550KB）下载：下载全尺寸图像图4. 1980年代（a）、2010年代（b）以及有机碳变化情况（c）下，各类环境和管理指标的Shapley值，这些指标是按照平均绝对Shapley值从高到低排序的（d–f）。数据点的颜色是根据该数据点对应的特征值来确定的。AI：干旱指数；MAT：年平均气温；MAP：年平均降水量；CTI：复合地形指数。下载：下载高分辨率图像（857KB）下载：下载全尺寸图像图5. 1980年代（aa–al）、2010年代（ba–bl）以及有机碳变化情况（ca–co）下，各类连续型环境和管理指标对应的Shapley值变化情况。主图中的颜色代表了数据点的密度，而边缘图则展示了各指标值以及Shapley值的频率分布情况。到了2010年代，各指标的重要性排名与1980年代大致相同。黏土含量依然是最重要的指标，其次是年平均气温、干旱指数、年平均降水量以及旱地轮作。土壤pH值的重要性有所下降，海拔的重要性也同样有所降低。化肥氮输入量及其变化趋势的重要性也有轻微下降，不过它们与有机碳含量的负相关关系依然存在（见图5bi、bj）。秸秆碳输入量以及灌溉措施对有机碳含量的预测贡献则有所增加。灌溉强度越高，也就是每个网格单元中处于灌溉状态的面积比例越大，那么预测的有机碳含量就会越高。关于更详细的空间分布特征，补充图1展示了1980年代和2010年代各网格级别上各类指标的主导类型分布情况。两个随机森林模型都表明，土壤属性和气候指标始终具有很高的重要性，而化肥氮输入量以及旱地轮作等农业管理措施也对有机碳含量的预测有着显著的贡献。各个指标与有机碳含量之间关系的方向——无论是正向还是负向——在两个时期都基本保持一致。不过，从1980年代到2010年代，农业管理指标的重要性总体上有所下降，这一现象需要谨慎解读。它有可能是管理水平整体提升以及这三十年间管理指标的空间变异程度逐渐降低所导致的，也可能是数据分辨率存在限制所致。3.3.SOC变化与各指标之间的关系
SOC的时间变化是通过1980年代的伪样本值与实际样本值之间的差异来计算的。因此，它代表的是基于伪样本对推导出的变化，而非直接观测到的重复测量比率。共有15个连续型指标和1个二元轮作变量（旱地轮作）被用于为伪样本对数据集构建随机森林模型。随机森林模型的详细评估指标列于补充表5中。各指标的重要性在图4中显示，而所有连续型指标的SHAP值响应情况则展示在图5中。初始黏土含量对SOC变化具有负向影响。旱地轮作是影响SOC变化的第二大因素，其后是气候指标。尽管采用旱地轮作的农田中SOC含量较低，但与其他轮作系统相比，这些样本的SOC变化幅度通常更大（图3d–f）。秸秆碳输入的初始水平及其时间变化，以及化肥氮输入的变化，是最重要的管理指标。初始秸秆碳输入对SOC变化的负向影响与初始黏土含量类似，但程度较轻。初始输入量过高的负面影响可理解为观测数据集中的空间关联性，而非直接证据表明更大输入量会降低某地的SOC变化。化肥氮输入的变化对预测的SOC变化有正向影响，其转折点约为15克氮/平方米·年（图5cl）。其他管理指标，包括1971年至2010年的化肥氮输入梯度以及初始化肥氮输入量，也对SOC变化预测有显著贡献。气候指标同样具有重要性，其对SOC变化预测的影响方向十分明确。1980年代的MAP、MAT、AI值（即1971–1980年的平均值）及其变化均对预测的SOC变化有正向影响。总体而言，土壤属性和气候指标是影响SOC变化预测的主要因素。不过，农业管理措施也发挥着重要作用，尤其是旱地轮作。

3.4 轮作系统对SOC与管理指标之间关系的影响
SOC的轮作特异性特征表明，SOC与管理指标之间的关系因轮作系统不同而存在差异。管理指标与轮作变量的SHAP交互值进一步揭示了这些关系。图6展示了管理指标与SHAP交互值之间的斯皮尔曼相关系数，以及部分显示从1980年代到2010年代明显特征或显著变化的依赖图。在中国北方的单季作物旱地地区，半干旱和干旱地区的化肥氮输入呈现正的SHAP交互值（图6b、g），表明单一作物种植与氮肥施用之间存在协同效应。相比之下，在中国东北著名的黑土区，单季旱地（东北）系统的正交互值在2010年代减弱甚至转为负值（图6h）。旱地轮作系统在两个十年间均表现出化肥氮输入的持续正交互值（图6d、i），说明在这种集约化种植系统中，氮肥施用与SOC积累仍存在持续协同效应。双季水稻轮作系统也出现了类似模式，但在1980年代至2010年代间，以水稻为主的系统的这种效应有所减弱（图6e、j）。在中国西南部的轮作系统中，化肥氮输入的SHAP交互值为负，即单季水稻-旱地（西部）系统以及谷物/根茎类作物系统（图6f、k）。

下载：下载高分辨率图像（1017KB）
下载：下载全尺寸图像

图6. 管理指标通过SHAP交互值对预测SOC的轮作特异性影响。（a）管理指标与其SHAP交互值之间的斯皮尔曼相关系数（b–u为R值）。韦尔奇t检验用于判断热力图中每个网格单元对应的目标轮作系统与其他轮作系统之间是否存在显著的SHAP交互值差异。（b–q）SHAP交互值与管理指标的依赖图。彩色线和灰色线分别为包含/排除某一轮作系统的95%置信区间内的LOWESS回归线。*P<0.05；**P<0.01；***P<0.001。

关于秸秆碳输入，中国北方的旱地轮作系统呈现出正-负-正的SHAP交互值曲线（图6n、s）。而单季旱地（东北）系统则显示出正的SHAP交互值（图6m、r）。双季水稻轮作系统在2010年代秸秆碳输入的交互值逐渐变为负值（图6t），表明在这种集约化系统中，过量施用秸秆可能对SOC产生不利影响。中国西南部的谷物和根茎类作物旱地系统中，秸秆碳输入的交互值总体呈下降趋势（图6p、u），说明在这种主要分布在山区的系统中，秸秆对SOC的影响较为有限。

SHAP交互分析进一步揭示了管理指标在不同轮作系统中的随机森林模型表现，以及管理措施与轮作方式之间是否存在协同或拮抗效应。

4. 讨论

4.1 环境调控下表层土壤SOC空间变化与管理实践之间的关联
我们的分析表明，中国农田中SOC的空间变化是由自然环境条件与数十年来集约化农业管理共同作用形成的复杂格局。我们的研究结果证实了这样一种公认观点：大规模的SOC分布格局本质上受环境因素控制（Wiesmeier等人，2019）。无论是在1980年代还是2010年代的模型中，黏土含量、AI值和MAT值都具有重要意义（图4），这表明气候和土壤质地决定了SOC储存的基本边界条件。这些因素决定了生物量生产力与微生物分解作用之间的平衡——前者决定碳输入量，后者则控制碳的输出量。土壤pH值与其SHAP值之间呈负相关，这一现象可能反映出稻田中较低的pH值限制了微生物的酶活性（Averill和Waring，2018）。

关于SOC空间变化最关键且反直觉的发现是，无论是哪个年代，化肥氮输入与SOC变化之间都存在持续的负相关关系（图5ai、bi）。这一现象挑战了“化肥通过提高作物产量从而增加碳输入，必然导致更高的SOC储存”这一简单假设。这一矛盾现象可从不同角度加以解释。一方面，中国平原地区的农田在从天然森林或草原转变为农田后，经历了显著的碳损失（Amelung等人，2020；Yue等人，2020），这一过程是由长期耕作和粮食生产驱动的。然而，这些地区也是最早受益于农业技术进步的地区，自20世纪70年代起，这里开始大量使用合成氮肥（Jiao等人，2018；Li等人，2013）。另一方面，这一关系也可能与“启动效应”的长期表现有关——该效应指的是添加易分解的基质，如施肥刺激产生的新鲜根系分泌物，会改变原有稳定SOC的分解速率（Bernard等人，2022；Fontaine等人，2007）。虽然氮肥的添加有时可以通过缓解微生物的氮需求来抑制启动效应（即负启动效应），但在那些长期从事集约化农业且氮输入量较高的生态系统中，微生物可能已不再受氮限制，尤其是在中国这种氮盈余极为严重的地区（Schulte-Uebbing等人，2022）。在氮素充足的条件下，进一步增加氮肥用量，再加上植物生长增强带来的易分解碳源，可能会加速微生物活动，从而导致原有SOC的更快分解，因为可能存在其他限制微生物活动的养分因素（Chen等人，2014；Zheng等人，2024）。这种持续的正启动效应降低了SOC含量的稳定状态，使得氮肥投入较多的农田，如中国北方的旱地轮作系统，SOC平衡水平更低。图5bi中确定的约8克氮/平方米·年这一转折点，可能代表着一个临界阈值，超过这一阈值后，启动效应的负面作用开始占据主导地位，超过合成氮肥带来的益处，进而对预测的SOC产生负面影响。这与化肥氮输入梯度对预测SOC的负面影响相似（图5aj、bj），因为化肥氮输入与其时间梯度之间存在很强的正相关性（见补充图3），这意味着氮肥投入量较高的农田，其氮添加速率始终更高。

秸秆碳输入与SOC变化之间也呈负相关关系（图5ak、bk），这与人们认为作物残体回归土壤能增加土壤碳的预期相悖。然而，秸秆碳输入的时间梯度与其对SOC影响之间的正相关关系，为理解这一做法的时间动态提供了重要线索。自20世纪80年代起，中国逐步提倡直接将秸秆还田（Liu和Li，2017）。与合成氮肥的施用类似，最早实施秸秆还田政策的地区往往是农业发展水平和粮食产量较高，但SOC含量较低的地区。由于20世纪80年代的秸秆碳输入量普遍较低，其直接结果是，高产粮食带来的高秸秆碳输入反而对预测的SOC产生了负面影响，导致SHAP值整体为负。经过三十年的粮食产量提升以及持续过量的氮肥施用，这种负相关关系变得更为复杂。根据先前的研究（Li等人，2018），秸秆还田的效果与时间密切相关——连续3–15年进行秸秆还田可以显著提高SOC含量。因此，这种负相关关系可能是由协同启动效应造成的：持续将富含碳但氮含量低的新鲜秸秆施入氮素充足的土壤中，为微生物的代谢活动和物质更新提供了充足能量，这可能会加速新增秸秆和原有SOC的分解，以满足微生物的营养需求（Fontaine等人，2011；Zheng等人，2024）。不过，如果整个研究期间秸秆碳输入总量呈上升趋势，即正的秸秆碳输入梯度，那么它会对预测的SOC产生正向影响（图5bl），这说明持续的秸秆还田可能抵消了启动效应的负面影响，使得SOC积累量高于那些秸秆还田较少的地区，这一结果与中国农业政策的走向一致。

总之，分析不同时间点上农田SOC的空间变化，反映了在环境调控背景下，随着农业集约化的发展，表层土壤碳分布格局发生了长期变化。各种管理措施通过对土壤化学性质产生空间上不均匀且非线性的影响，对表层农田SOC的重新分布起到了重要作用。

4.2 农业集约化背景下与管理相关的表层土壤SOC积累
为了进一步探讨这一问题，我们生成了伪土壤样本对，以识别SOC变化的关键驱动因素，即SOC变化的整体时间趋势。在本研究涵盖的三十年间，我们的研究结果表明，中国农田的表层土壤SOC呈现明显的累积趋势。这一发现与众多全国性研究的结果一致（Lin等人，2023；Ren等人，2021；Zhao等人，2018；Zhou等人，2025），这些研究均指出，在农业集约化发展的过程中，中国农田的表层土壤SOC总体上有所增加，具备成为净碳汇的潜力。我们的分析不仅揭示了这一总体趋势，还阐明了推动这一积累过程的具体管理机制。

农业管理有望提升土壤有机碳库的容量，尤其是活性颗粒有机碳的比例（Mitchell等人，2018）。与化肥氮输入对预测SOC的负面影响不同，我们的研究结果显示，研究期间化肥氮输入的变化显著推动了SOC的积累（图4c和图5cl）。考虑到正启动效应和碳饱和缺陷理论的存在，这一发现表明，氮素可用性的持续增加在空间上非线性地提升了SOC含量，从而克服了此前讨论过的合成肥料带来的负面影响。一方面，氮素可用性的增加直接促进了作物生产力的提升（Xinyue Zhang等人，2025），使得通过根系生物量返回土壤的碳量大幅增加。再加上秸秆还田政策的推动，更多的新碳源成为SOC积累的主要动力，尤其是在SOC含量较低的土壤中。另一方面，先前的研究表明，氮肥的添加会导致土壤酸化（D. Zhou等人，2025），而这一现象在中国的主要农田中十分普遍（Guo等人，2010）。氮素引起的pH值下降促进了土壤有机碳的积累（Lu等人，2022年），并通过降低POC与MAOC的比值来减轻氮素添加对土壤有机碳稳定性的负面影响（Xinsheng Zhang等人，2025年）。再加上秸秆还田导致的土壤酸化加剧现象（Liang等人，2023年），进一步强化了氮肥施用对土壤有机碳随时间变化的正向影响。此外，SHAP值对氮肥投入变化的非线性响应（图5cl）可能表明存在一个阈值（约15克氮/平方米·年），超过该阈值后氮素促进碳封存的效果开始减弱，说明由于氮素过剩导致施肥效率下降。土壤有机碳的初始水平通常被视为调控其动态变化的重要因素，这一观点得到了土壤碳饱和度缺陷理论的支持。该理论认为土壤具有有限的稳定和储存土壤有机碳的能力，而这种能力受到与矿物质结合的有机碳饱和度的限制，这一最大容量完全由土壤矿物基质的化学和物理性质决定（Georgiou等人，2025年）。从理论上讲，MAOC的含量受其饱和度缺陷的限制（King和Sokol，2025年），并且可以通过管理措施和气候条件得到提升（Georgiou等人，2025年）。我们构建了一个随机森林模型，将初始土壤有机碳水平作为特征变量，并进行了SHAP分析，结果见补充图4。分析表明初始土壤有机碳是预测土壤有机碳变化的最重要因素，且与其SHAP值呈负相关。然而，由于土壤有机碳变化是通过伪土壤有机碳与初始土壤有机碳的差值计算得出的，二者之间的负相关关系可能是由数学上的回归到均值偏差所导致的（Slessarev等人，2023年）。我们进行了回归到均值偏差诊断（补充图5），发现两者之间的关系虽然有所减弱，但仍然显著为负，这为国家尺度上碳饱和度缺陷理论提供了验证。此外，在排除初始土壤有机碳之后，初始黏土含量成为预测土壤有机碳变化的最重要因素。初始黏土含量与其SHAP值在土壤有机碳变化方面的负相关关系（图5ca）进一步支持了该理论，尤其是考虑到MAOC与黏土及裂隙含量之间存在正相关关系。

4.3 中国不同轮作制度下表土土壤有机碳的空间和时间变化

通过分析不同轮作制度下的研究结果，我们可以看出农业管理方式如何在中国形成了不同的区域土壤有机碳变化轨迹。图7展示了20世纪80年代、21世纪10年代不同轮作制度下的氮肥投入量和秸秆碳投入量以及这些数值在此期间的变化情况。由于当地经济状况、地形、气候条件的差异，不同轮作制度下的管理方式存在明显区别，进而影响了耕作策略。SHAP交互作用分析的结果进一步揭示了管理差异对土壤有机碳动态的影响，表明不同轮作制度是在何种程度上增强了或抵消了管理措施对土壤有机碳空间和时间变化的影响。

下载：下载高分辨率图像（319KB）
下载：下载全尺寸图像
图7. 20世纪80年代（a）、21世纪10年代（b）不同轮作制度下的氮肥投入量与秸秆碳投入量，以及这两者之间的变化（c）。实心点表示各轮作制度的中位数，横线和竖线则表示四分位数范围（第25百分位到第75百分位）。

华北地区的主要轮作制度是以小麦-玉米或小麦-大豆的密集轮作为例，属于典型的碳封存不足型轮作模式。在这一制度下，研究开始时的土壤有机碳水平较低，但在接下来的三十年里土壤有机碳含量出现了显著上升（图3d、f）。由于该地区基础肥力低、地形平坦且人口密度高，因此存在较大的碳饱和度缺陷，但在粮食生产中起着重要作用。因此，该地区对政策推动的农业集约化带来的大量碳和氮输入反应十分敏感。与其他轮作制度相比，氮肥投入量与其SHAP交互作用值之间存在正相关关系（图6d、i），这说明土壤碳的不足在充足的氮素和碳输入作用下能够促进碳封存，且其效果超过了负向的初始效应。不过，由于SHAP交互作用值与氮肥投入量变化之间存在负相关关系（补充图6b），这可能表明该地区的氮素已经过剩。此外，20世纪80年代秸秆碳投入量的SHAP交互作用值总体为正（图6n和补充图6c），说明在该碳缺乏的轮作制度中，碳输入有助于土壤有机碳的积累。关于华北平原农田土壤有机碳的积累趋势以及秸秆还田和氮肥施用的作用，之前的研究也有相关论述（Han等人，2018年）。中国西北部的单季水稻轮作制度也呈现出类似的趋势，尽管20世纪80年代的初始土壤有机碳水平较低，但由于农业集约化程度相对较低且气候更为干旱，其土壤有机碳的上升幅度虽然显著，但低于华北地区的轮作制度。

中国南方的双季水稻轮作制度则属于高投入、高保持型的轮作模式。这一制度在20世纪80年代就已经具有较高的土壤有机碳水平，而且令人惊讶的是，其土壤有机碳含量还在持续增加（图3d、f）。较高的初始土壤有机碳水平是由于淹水水稻种植造成的长期厌氧环境抑制了有机物的分解速度。整个研究期间，该地区的氮肥投入量一直较高（图7a、b），再加上轮作制度的协同作用（图6a），表明施肥对作物生长和土壤环境维护有着积极作用。此外，尽管该地区的秸秆还田率相对较低，且管理指标的变化幅度较小（图7），但在这类集约化管理模式下，合理的氮素和碳输入依然能够在土壤有机碳初始储量较高的水稻系统中实现碳封存潜力。通过冬季休耕和种植绿肥等策略，进一步提升了碳封存潜力。不过，氮肥投入量与轮作制度之间的相关性逐渐减弱，同时秸秆碳投入量与轮作制度之间的相关性也由正转负，这可能表明施肥和秸秆还田存在最佳效果的最佳阈值。相比之下，单季水稻-旱地（东部）制度尽管在管理上投入较多，但并未带来显著的土壤有机碳积累，这一有限的效果可能意味着农业投入过度。

4.4 对农田土壤碳管理和估算的启示

本研究的结果为制定针对性策略、提升中国不同农业区域的土壤碳汇能力提供了有力的数据支持。不同轮作制度下土壤有机碳变化轨迹的差异，要求我们采取针对性的政策干预和农场管理措施。对于那些土壤有机碳含量较低的旱地系统，研究结果表明，这些地区是中国实现农业碳封存的最大且最紧迫的机遇。因此，政策制定应优先考虑采用能够减少碳损失并避免土壤酸化等不良后果的管理措施。对于以水稻种植为主的系统，单季水稻系统土壤有机碳积累有限的现状表明，需要调整耕作策略，以在农业投资与维持农田肥力之间找到平衡，例如像双季水稻轮作制度那样，在休耕期延长冬季覆盖作物种植或施用绿肥。

本研究强调了氮肥施用和秸秆还田在调节表土土壤有机碳变化中的重要性，表明这两种措施在农田碳管理中具有显著效果。尽管这两种措施在保障粮食生产和实现碳封存方面都发挥着作用，但它们也带来了诸多环境问题，如水污染、空气污染、过度施用氮素导致的土壤酸化和退化，以及直接掩埋作物残体所带来的病虫害问题。在评估农业管理政策时，应综合考虑经济成本、作物产量、土壤环境效益以及对区域生态和人类生计的潜在影响。

本研究最重要的启示在于，土壤碳建模中必须纳入动态管理数据。通过引入反映农业管理实践随时间变化的指标，本研究明确了影响土壤有机碳变化的驱动因素，表明这些指标并非仅仅是补充信息，而是解释农业生态系统土壤有机碳动态变化的关键要素。因此，有必要投资于收集、整合高分辨率、具有空间明确性且能反映时间动态的农业管理实践数据的相关系统。

4.5 研究局限性

本研究存在若干局限性，这些局限性可能会影响研究结果。伪样本配对方法虽然克服了全国土壤调查中非重复采样的弊端，有助于了解土壤有机碳积累的驱动因素，但为了更准确地模拟农田土壤有机碳的动态变化，尤其是后续的全国土壤调查数据，应结合更详细的原位采样以及管理条件记录。目前所使用的氮肥投入量和秸秆碳投入量数据来自分辨率较低的统计数据，而且氮肥数据中未包含有机肥的投入量。由于早期缺乏空间数据，基准期的灌溉情况是根据2000年的数据估算的。这一假设在2000年之前灌溉面积大幅扩张的情况下会引入不确定性，也可能削弱对基准期灌溉效应的判断。在解读所报告的关联强度和机制时，应考虑这些局限性。不过，它们不会显著影响土壤有机碳动态变化的整体趋势和方向。利用遥感技术获取的其他详细网格化数据能够提供更多关于农业活动特征的信息，从而提高研究的可解释性，避免土壤样本与指标之间出现不匹配的情况，比如用于估算作物残体的光谱指数（Dong等人，2024年），这些都将进一步提升我们对农业在过去、现在和未来土壤有机碳动态变化中作用的理解。

CRediT作者贡献说明
岳浦：写作——审阅与编辑，写作——初稿撰写，可视化，软件应用，方法论，正式分析。
杨林：写作——审阅与编辑，监督，方法论，概念构建。
崔文凯：写作——审阅与编辑，方法论。
李翔：写作——审阅与编辑。
沈飞雪：写作——审阅与编辑。
周成虎：监督。