张海星 ，赵智强 ，王瑞锋 ，刘 静 ，楚天舒 ,4※
（1. 中国农业大学 资源与环境学院,北京 100193；
2. 中国农业大学 工学院,北京 100083；
3. 上海市建设用地和土地整理事务中心,上海 200003；
4. 自然资源部大都市区国土空间生态修复工程技术创新中心,上海 200003）
全球气候变化与极端天气频发对国家安全和经济社会可持续发展带来了严重威胁，已成为全球性挑战。2015 年，联合国可持续发展目标和《巴黎气候协定》共同倡议，各国努力建设可持续低碳经济。2019 年中国温室气体排放的二氧化碳当量为140 亿t，占到全球排放量的27%[1]。从生产部门分析，农业生产温室气体排放量占人为温室气体排放总量的23%，但耕地也是巨大的碳库，其固碳量约占陆地生态系统总固碳量的21%[2]。随着农业绿色发展的大力推进，探索提高耕地土壤有机碳含量的途径，显得至关重要。
近年来，国内外围绕农田土壤有机碳展开了大量研究。HAMID 等[3]研究表明，在长期耕作条件下，种植年限的增加会降低农田总有机碳的可用性，其中2018 年的有机碳含量较2014 年降低了45.32%。而ZHOU 等[4]在中国黄土高原区研究发现，相较于小麦—小麦轮作，1984—2010 年苜蓿地表层土壤有机碳含量增加了24.4%。并且在退耕还草中，长期适度放牧也可以促进有机碳的积累[5]。董林林等[6]通过10 a 定位试验发现，稻麦秸秆还田可增加太湖地区水田土壤有机碳含量0.18～0.46 g/(kg·a)并增强其稳定性。从中国科学院桃源农业生态实验站25 a双季稻定位试验结果可知，氮磷钾平衡施肥处理组的土壤有机碳含量较不施肥处理组高出7.7%[7]。相对于不施用氮肥，长期施用氮肥提高了农田土壤有机碳含量5.2%～12.0%，但其时间响应取决于土壤类型，即在有机碳含量较高的土壤类型中具有放大效应[8]。在西班牙穆尔西亚省果园中，长期采用免耕与绿肥的组合方式可抵消耕作对土壤团聚体破裂的影响，并促进新团聚体的形成，与仅免耕相比，免耕与绿肥组合可提高有机碳含量14.0%[9]。在东北黑土区，相对于传统翻耕模式，保护性耕作能显著增加土壤表层有机碳含量25.8%[10]，并能增加大团聚体比例及团聚体稳定性[11]。总的来说，多种生产措施可促进耕地土壤固碳增汇，而如何系统性评估土壤固碳的相关措施，找出主要影响因素，还有待进一步研究。
长江中下游地区作为中国传统农业产区，如何实现耕地固碳增汇，是该地区农业绿色发展的关键问题。当前缺乏系统性探究该地区耕地土壤有机碳密度变化特征与驱动因素的分析。因此，本研究拟以长江中下游地区为研究区域，收集与整理耕地土壤有机碳密度长期观测数据，量化不同耕地土壤有机碳密度变化率的差异特征，应用随机森林模型探究土壤有机碳密度变化率的主要驱动因素，解析耕地土壤有机碳密度变化率对环境和管理因素的非线性响应关系，以期为耕地固碳增汇提供参考。
1.1 数据来源与筛选方法
本研究选取“有机碳”“有机质”“长期定位”“耕地”等关键词，分别在中国知网（CNKI）和Web of Science文献数据库进行检索，收集与整理了公开发表的有关长江中下游地区耕地土壤有机碳密度变化的研究论文。进而，按照如下标准进行文章筛选：1）研究区域为长江中下游地区，从省级行政区划上涵盖：上海、江苏、浙江、安徽、江西、湖北、湖南；
2）试验方式为大田试验，试验年限不少于2 a，试验地点、时间、基础条件清楚，优选长期定位试验；
3）试验土样为表土（0～20 cm），土壤有机碳密度数据准确易获取。此外，若文章中部分土壤基本信息缺失，则在世界土壤数据库（Harmonized World Soil Database,version 1.2）中按经纬度查询以获取[12]。若文章中气象基本信息缺失，如年降水量和年平均气温等，则在东安格利亚大学气候研究中心气象数据库（Climate Research Unit）查询以获取[13]。
本研究按照上述方法筛选出长期试验（表1），提取试验年限、地点、耕地基本信息（耕地利用类型、试验年限、土壤pH 值、黏粒含量、土壤有机碳含量）、作物基本信息（轮作模式、化肥氮量、有机肥氮量、秸秆氮量）、气象基本信息（年平均降水量和年平均气温），建立长江中下游地区耕地表土有机碳密度变化率数据库。其中，化肥氮量、有机肥氮量、秸秆氮量是由周年施用量乘以氮含量计算获得。
表1 试验基本信息Table 1 Basic information of experiment
1.2 耕地土壤有机碳密度变化率计算方法
若文献未直接给出土壤有机碳密度变化率数据，则依据《2006 年IPCC 国家温室气体清单指南2019 修订版》中方法计算获取，具体如下：
式中ΔCSO为耕地土壤有机碳密度变化率，kg/(hm2·a)；
CSO,t为第t年耕地土壤有机碳密度，kg/(hm2·a)；
CSO,0为初始年耕地土壤有机碳密度，kg/(hm2·a)；
t为试验年限，a；
c为耕地土壤有机碳含量，g/kg；
b为耕地土壤容重，g/cm3；
h为耕层厚度，cm；
100 为单位转化系数。
1.3 随机森林算法
随机森林算法作为典型机器学习算法之一，具有良好的处理非线性问题的能力，广泛应用于农业预测与建模。在本研究中，将耕地基本信息、作物基本信息、气象基本信息作为特征变量，而将表土有机碳密度变化率作为目标变量，采用随机抽样的方法确定训练集与测试集。通过自助法采样获得构建N棵树所需的n个子集，每次未抽中的数据称为袋外数据，用于内部误差估计和特征变量的重要性计算。对测试集的目标变量真值与预测值计算平均绝对误差、均方误差和均方根误差，并进行误差估计。选取对特征变量j随机排列后预测值平均精确率减少值（mean decrease accuracy,MDA）表征特征变量的重要性并对其排序，筛选出影响耕地土壤有机碳密度变化率较大的特征变量作为驱动因素，基础原理详见式（2）。最终构建耕地土壤有机碳密度变化率的预测模型。算法采用Python 语言和scikit-learn 机器学习工具实现。
式中ij表示特征变量j的MDA 值；
s表示原始数据集使用随机森林回归的模型分数；
K表示对数据集中特征变量j数据进行重新随机排列的总次数，sk,j表示在第k次对特征变量j随机排列，原始数据集其他列不变的情况下，使用随机森林回归的模型分数；
R2表示预测的决定系数；
ytrue表示目标变量的实际值；
ypred表示随机森林回归模型基于各特征变量得出的预测值；
ytrue.mean表示目标变量实际值的均值。
1.4 施肥策略优化
2022 年11 月农业农村部发布《到2025 年化肥减量化行动方案》，将减少化肥用量、提高有机肥用量作为目标任务，促进农业绿色转型。江苏省作为国家农业绿色发展先行区，对于长江中下游地区农业绿色低碳发展具有引领示范作用。因此，本文以江苏省为优化对象，选取典型水田和旱地轮作模式，将随机森林模型与线性优化法相结合，以土壤有机碳密度变化率最大为目标函数，以氮肥总量和有机无机肥配施比例为约束条件，仿真分析出适宜的有机肥氮量和化肥氮量，以期协同实现耕地土壤固碳与化肥减量化，具体模型如下：
式中ΔCSO,pred(x,y)为随机森林模型预测不同情况下耕地土壤有机碳密度变化率，kg/(hm2·a)；
x为化肥氮量，kg/hm2；
y为有机肥氮量，kg/hm2；
a为氮肥总量，kg/hm2，参考农业农村部发布的作物科学施肥指导意见与当地推荐施肥量确定；
r为有机肥氮量占氮肥总量比例，%；
rmin为有机肥氮量占氮肥总量比例的最小值，%；
rmax为有机肥氮量占氮肥总量比例的最大值，%。根据大田作物常规施肥策略与作物养分需求规律，若有机肥施用量过高，大田作物容易出现“贪青”现象，作物生殖生长期延后，直接影响作物产量。并且，兼顾长江中下游地区农业面源污染防控需要，参考《NY/T 4 163—2022 稻田氮磷流失综合防控技术指南》，本研究将有机肥氮量占氮肥总量比例的最大值设置为30%。
根据策略优化结果，结合《NY/T 3 877—2021 畜禽粪便土地承载力测算方法》，定量评估江苏省畜禽养殖是否能提供充足的有机粪肥资源支持优化后的施肥策略。
2.1 耕地土壤有机碳密度变化率特征分析
通过对59 处290 组长期定位试验数据进行统计与分析可知，水田和旱地土壤有机碳密度变化率范围分 别 为-1 548.15～3 577.10 kg/(hm2·a)和-261.89～3 245.01 kg/(hm2·a)。从土地利用类型来看，水田与旱地的土壤有机碳密度变化率差异不显著（P=0.85）。从轮作模式来看（表2），针对水田，水稻—小麦、水稻—油菜的土壤有机碳密度变化率较高，但数据离散程度大；
单作水稻的土壤有机碳密度变化率相对较低。针对旱地，蔬菜—蔬菜和油菜—棉花的土壤有机碳密度变化率较高，两者间无显著性差异；
而玉米—玉米的土壤有机碳密度变化率最低。
表2 耕地土壤有机碳密度变化率Table 2 Cultivated land soil organic carbon density change rate (kg·hm-2·a-1)
2.2 随机森林模型的构建
本研究以长期定位试验数据为基础，运用随机森林算法，调节模型生成的决策树数目超参数N=10 000，其余取默认值，获得水田和旱地的土壤有机碳密度变化率模型（图1）。其决定系数（R2）分别为0.63 和0.83，预测值与真实值高度相关[68]，表明随机森林算法在耕地土壤有机碳密度变化率模型的构建上取得了较好的效果。
图1 耕地土壤有机碳密度变化率（ΔCSO）的真实值与预测值散点图Fig.1 Scatter diagrams of the true and predictive value of cultivated land soil organic carbon density change rate(ΔCSO)
2.3 驱动因素分析
2.3.1 驱动因素重要性排序
结合耕地土壤有机碳密度变化率的驱动因素重要性排序与系统聚类结果（表3）可知，对于水田，有机肥氮量、土壤pH 值、化肥氮量、秸秆氮量对有机碳密度变化率的影响较大，而试验年限、黏粒含量、年平均降水量、土壤有机碳含量、年平均气温、水旱轮作、复种指数不是主要影响因素。对于旱地，有机肥氮量对有机碳密度变化率的影响最大，化肥氮量次之，而土壤有机碳含量、年平均降水量、试验年限、土壤pH 值、年平均气温、黏粒含量、复种指数、秸秆氮量的影响较小。总的来说，有机肥氮量和化肥氮量是长江中下游地区耕地土壤有机碳密度变化率的主要驱动因素。
表3 耕地土壤有机碳密度变化率的驱动因素重要性排序Table 3 Sorting of the driving factors of cultivated land soil organic carbon density change rate
2.3.2 驱动因素的边际依赖性分析
就水田而言，随着有机肥氮量的增加，土壤有机碳密度变化率呈现“波动增长—下降”趋势，最大值约为987 kg/(hm2·a)（图2a）。随着土壤pH 值从弱酸性至中性再至弱碱性，土壤有机碳密度变化率呈现出阶梯式增长趋势；
当土壤pH 值为7.8 时，土壤有机碳密度变化率最高可达941 kg/(hm2·a)（图2b）。当不施用化学氮肥时，土壤有机碳密度变化率极低；
随着化肥氮量增加，土壤有机碳密度变化率快速增加，而后在616～850 kg/(hm2·a)范围间呈波动变化（图2c）。另外，随着秸秆氮量的增加，土壤有机碳密度变化率呈现波动增长趋势，最大值约为904 kg/(hm2·a)（图2d）。
图2 水田驱动因素的部分依赖图Fig.2 Partial dependence plots for the paddy field driving factors
科学调控土壤pH 值和秸秆还田量有利于耕地土壤固碳（图2b 和图2d）。在土壤pH 方面，虽然弱碱性土壤有利于固碳，但土壤pH 值是影响农作物产量的重要因素之一，土壤碱化可使农作物减产50%以上[69]。因此，本研究建议长江中下游地区合理调控水田土壤pH 值至中性，可兼顾实现耕地固碳与作物高产目标。在秸秆还田量方面，秸秆虽是重要碳源，但随着秸秆还田量的增加，稻田甲烷排放量却呈现线性增加趋势，其增加的温室效应会抵消土壤固碳的减排效益[70]。因此，未来研究还需要进行大量的田间试验与模型仿真分析，探究适宜的秸秆还田量，协同实现土壤固碳与温室气体减排，以促进农田生态系统绿色低碳发展。
对于旱地而言，当不施用有机氮肥时，土壤有机碳密度变化率极低；
随着有机肥氮量增加，土壤有机碳密度变化率呈现出“增长—平稳”趋势，最高值约为1 425 kg/(hm2·a)（图3a）。对于一年一季的轮作模式来说，随着化肥氮量的提升，土壤有机碳密度变化率也随之增长，最大值在659 kg/(hm2·a)左右；
而对于一年两季的轮作模式来说，化肥氮量呈现“增长—波动下降”趋势，当化肥氮量超过300 kg/hm2后，土壤有机碳密度变化率迅速提高，最高可达824 kg/(hm2·a)（图3b）。
图3 旱地驱动因素的部分依赖图Fig.3 Partial dependence plots for the dry land driving factors
从上述分析可知，有机肥氮量和化肥氮量是长江中下游地区耕地土壤有机碳密度变化率的重要驱动因素。因此，进一步对有机肥氮量和化肥氮量进行双驱动因素边际依赖性分析，评估两者是否协同促进耕地土壤有机碳密度变化率。从双驱动因素部分依赖图结果（图4）可知，当有机肥氮量或化肥氮量过低时，耕地土壤有机碳密度变化率很低，不利于耕地土壤固碳；
当有机肥氮量或化肥氮量过高时，耕地土壤有机碳密度变化率并非最大值，这说明有机肥或化肥过量施用，也不是土壤固碳的最佳途径；
相较于单驱动因素分析结果（图2 和图3），有机肥氮量和化肥氮量的双驱动因素分析中，土壤有机碳密度变化率的最大值更大；
这说明有机肥和化肥可协同提升耕地土壤有机碳密度。在水田中，有机肥氮量和化肥氮量分别为346～397 kg/hm2和422～533 kg/hm2，土壤有机碳密度变化率达到最大值，且较为稳定，约为1 156 kg/(hm2·a)。而在旱地中，有机肥氮量和化肥氮量分别在219～284 kg/hm2和338～394 kg/hm2范围内，土壤有机碳密度变化率稳定在最大值，约为1 654 kg/(hm2·a)。由此可见，有机无机肥配施将有效促进耕地土壤固碳。
图4 有机肥与化肥氮量双驱动因素的部分依赖图Fig.4 Partial dependence plots of driving factors for joint effects of organic fertilizer nitrogen and chemical fertilizer nitrogen
2.4 施肥策略优化
在氮肥总量控制、有机无机肥配施比例约束情况下，采用随机森林模型与线性规划仿真分析可知：随着有机肥氮量占比增加，江苏省水田土壤有机碳密度变化率呈现“增加—波动”趋势；
当有机肥氮量占氮肥总量比例超过10%，水田土壤有机碳密度变化率在1 143～1 335 kg/(hm2·a)间波动（图5a）。旱地土壤有机碳密度变化率随有机肥占氮肥总量比例增加而增加，近似幂函数关系（R2=0.89）；
当有机肥氮量占比超过10%，水田土壤有机碳密度变化率增长趋缓，最高可达1 039 kg/(hm2·a)（图5b）。总的来说，有机肥氮量占氮肥总量比例控制在10%～30%，可稳定地促进耕地土壤固碳，同时助力化肥减量化行动。
图5 耕地土壤有机碳密度变化率的仿真分析结果Fig.5 Simulation results of cultivated land soil organic carbon density change rate
与此同时，当有机肥替代比例为30%时，2020 年江苏省耕地需要有机肥氮量为38.31 万t。核算可知，2020年江苏省畜禽养殖可提供有机肥氮量为42.20 万t，粪肥资源充足。
综合上述研究成果可知，有机无机肥配施是提升长江中下游地区耕地土壤有机碳密度的重要措施，并且当地拥有充足的粪肥资源。因此，建议长江中下游地区大力推广有机无机肥配施，并且有机肥氮量占氮肥总量比例控制在10%～30%，以助力化肥减量化行动，促进农田土壤固碳增汇，实现区域种养循环。
本研究应用随机森林模型与线性规划探究长江中下游地区耕地土壤有机碳密度变化率的驱动因素与提升途径，得到以下结论：
1）本研究整理与分析59 处长期定位试验数据发现，长江中下游地区水田和旱地土壤有机碳密度变化率范围分别为-1 548.15～3 577.10 kg/(hm2·a)和-261.89～3 245.01 kg/(hm2·a)。从土地利用类型来看，水田与旱地的土壤有机碳密度变化率差异不显著（P=0.85）。
2）对于水田而言，有机肥氮量、土壤pH 值、化肥氮量、秸秆氮量对土壤有机碳密度变化率的影响较大，平均精确率减少值（mean decrease accuracy,MDA）分别为0.49、0.32、0.23、0.15，而试验年限、黏粒含量、年平均降水量、土壤有机碳含量、年平均气温、水旱轮作不是主要影响因素。对于旱地而言，有机肥氮量对土壤有机碳密度变化率的影响最大、化肥氮量次之，MDA 值分别为0.98、0.10，土壤有机碳含量、年平均降水量、试验年限、土壤pH 值、年平均气温、黏粒含量、复种指数、秸秆氮量的影响较小。有机肥氮量和化肥氮量是长江中下游地区耕地土壤有机碳密度变化率的主要驱动因素。
3）有机无机肥配施、调控土壤pH 值至中性、秸秆还田有利于水田土壤有机碳密度变化率增加；
有机无机肥配施显著提高旱地土壤有机碳密度变化率；
耕地有机肥氮量占氮肥总量比例10%～30%为宜。