全球变暖背景下, 与大气中CO2浓度递增相关联的“碳问题”, 备受人们关注[1~3].农田作为陆地生态系统的重要组成部分, 受自然环境与人类活动因素的双重影响[4].在农田土壤碳库中, 土壤有机碳密度(soil organic carbon density, SOCD)表征单位面积一定深度的土层中土壤有机碳(SOC)的储量[5, 6], 它不仅是衡量土壤肥力和土壤质量的重要指标[7], 而且其微小变化会引起农田不同组分之间碳通量的变化, 进而改变农田生物地球化学循环过程[8~10].因此, 理解农田SOCD的变化特征对全球变暖背景下农田温室气体减排以及地球系统模式优化等具有重要意义.

近年来, 许多学者通过多源数据(如野外采集[11, 12]、文献整理[13, 14]和国家土壤普查数据[15~17]等), 对我国不同区域或站点尺度的农田SOC含量以及SOCD的变化特征, 进行了深入分析[18~20].在国家尺度上, 虽然一些研究表明我国农田SOC含量呈增加趋势[9, 21, 22], 但不同结果在增加速率及空间格局等方面存在较大的差异.农田SOCD是SOC含量、土壤容重和土层深度的综合反映, 许泉等[23]基于第二次土壤普查数据, 发现我国农田SOCD表现出明显的区域差异, 并在一定程度上受水热因素的影响.由于我国气候类型、下垫面条件及农业生产复杂多样, 加之不同来源的监测数据往往受空间覆盖面小、采样时间不统一、分析方法不一致等问题的制约[24], 难以系统认识农田SOCD长期且连续的变化规律.此外, 在农田SOCD影响因素的分析方面, 一些研究常采用相关分析的方法[25], 衡量SOCD与影响因素之间的线性关系.如果相关系数不显著, 则表明二者之间没有明显的线性关系, 但这并不能说明它们之间没有非线性关系.因此, 如何定量解析不同因素(如地理环境、气候因素、作物类型和土壤因素等)对农田SOCD的影响也是现阶段研究亟待解决的重要问题之一.近年来, 地理探测器方法的提出, 为分析各种复杂现象背后的驱动力以及多影响因子之间的交互作用[26], 提供了重要工具.目前, 该方法已广泛应用于地理学、生态学和环境科学等领域[27~29], 但涉及农业数据的大尺度空间应用研究还相对较少.

为此, 本研究以中国生态系统研究网络(CERN)典型农田生态站为例, 基于长期(2005~2015年)的SOCD数据, 分析农田SOCD的空间分异特征; 并利用地理探测器方法, 揭示农田SOCD空间分异的影响因素, 以期为农田温室气体减排以及地球系统模式优化等提供科学依据.

本研究所用的2005~2015年农田表层(0~20 cm)逐年SOCD数据来源于CERN的19个典型农田生态站.如表 1所示, 这些站点基本涵盖了我国不同气候区、不同土壤类型的农田生态系统, 主要位于松嫩平原、下辽河平原、太行山前平原、黄淮海平原、太湖平原、南方丘岗区、洞庭湖平原、川中丘陵区、黄土高原和荒漠-绿洲区等[24].各站点结合不同作物类型, 在每季作物监测1次的频度(收获期调查)上, 利用SOCD计算公式[式(1)], 得到各站点的SOCD.

式中, SOCD为土壤有机碳密度(kg·m-2), OM为土壤有机质含量(g·kg-1), γ为土壤容重(g·cm-3), H为土层厚度(cm).

为了探测不同因素对农田SOCD空间格局的影响作用, 本研究选取的影响因素包括地理环境、气候因素、作物类型和土壤理化性状等.具体如下：经度(Lon)、纬度(Lat)和高程(DEM); 气温(Tem)和降水量(Pre); 作物类型(Ctyp); 土壤类型(Styp); 土壤砂粒(Sand)、粉粒(Silt)和黏粒(Clay)含量; 土壤pH(SpH); 土壤碱解氮(AN)、速效磷(AP)和速效钾(AK)含量.其中, DEM信息一般不随时间发生变化, 通过地理空间数据云(http://www.gscloud.cn)的数字高程模型, 提取19个典型农田生态站点的DEM信息; 各站点作物生长季期间的平均气温和总降水量来源于国家气候中心的逐日气象栅格化资料, 从中提取19个典型农田生态站点的气温和降水量信息; 各站点的经纬度、作物类型和土壤数据来源于CERN, 作物类型为每年记录1次, 土壤数据一般每5a监测1次, 取平均值用于本研究的分析.对于连续型变量数据, 本研究借助ArcGIS软件的重分类功能进行离散化处理, 处理采用分位数分类原则, 共分为9类.

如前所述, 地理探测器是一种探测空间分异性, 并揭示其背后驱动力的统计学方法.与传统的统计学方法相比, 地理探测器具有明确的物理含义, 没有线性假设, 能客观地识别出自变量对因变量的解释程度, 它包括因子探测、交互作用探测、生态探测等功能.具体如下.

(1) 因子探测探测某因子对农田SOCD空间分异的解释力, 用q值度量[27], 它的计算见式(2).

式中, q为某指标的空间异质性; N为研究区全部样本数; σ2为指标的方差; h=1, 2, …, L, h为标识分区, L为分区数目.q值的值域介于0~1, q值越大, 表明该因子的解释力越强, 反之亦然.

(2) 交互作用探测用于判断两因子共同作用时是否会增加或减弱对农田SOCD的解释力[27].根据单因子作用时的q值与双因子交互作用时的q值, 将交互作用类型分为5类(表 2).

(3) 生态探测用于比较两因子对农田SOCD的空间分布影响是否有显著差异, 以F统计量衡量[27].本研究利用RStudio软件和地理探测器的R语言包(https://cran.r-project.org/web/packages/geodetector/index.html), 从地理环境、气候因素、作物类型、土壤因素等方面, 揭示农田SOCD空间分异的影响因素.

我国不同站点SOCD表现出明显的空间差异性(表 3), SOCD的平均值变化范围为0.83~4.97 kg·m-2.如表 3所示, 约有31.58%的站点SOCD在3.0 kg·m-2以上.其中, 海伦站SOCD最高, 介于4.89~5.21 kg·m-2之间, 平均值为4.97 kg·m-2; 其次为常熟站和桃源站, 分别为4.82 kg·m-2和3.70 kg·m-2.约有52.63%的站点SOCD在2.0 kg·m-2以下, 特别是策勒站SOCD最低, 平均值仅为0.83 kg·m-2.整体而言, 湿润季风区的SOCD[(2.96±1.16)kg·m-2]明显高于干旱半干旱地区[(1.58±0.39)kg·m-2], 特别是东北和华东地区的SOCD相对较高, 但同一区域内不同站点的SOCD仍表现出较大的差异性.

表 4显示了不同土地利用方式下典型农田SOCD的统计特征.从中可知, 水稻田SOCD最高, 其SOCD介于1.78~5.57 kg·m-2之间, 平均值为3.31 kg·m-2; 其次为大豆田, SOCD介于1.77~5.0 kg·m-2之间, 平均值为3.14 kg·m-2; 此外, 小麦田和玉米田SOCD也相对较高, 二者的平均值分别为2.44 kg·m-2和2.66 kg·m-2.相较于4种粮食作物而言, 种植经济作物的棉花田和花生田SOCD较低, 二者的SOCD的平均值分别为1.39 kg·m-2和1.30 kg·m-2.

2005~2015年我国不同区域SOCD变化曲线, 如图 1(a)所示.整体而言, 不同区域SOCD呈波动增加的趋势, 湿润季风区SOCD的增加速率明显高于干旱半干旱地区, 特别是华东地区SOCD呈显著增加趋势(P < 0.05), 增加速率为0.09 kg·(m2·a)-1.

进一步分析表明, 2005~2015年我国各站点SOCD也均表现出不同程度的增加趋势[图 1(b)], SOCD的增加速率介于0~0.33 kg·(m2·a)-1之间.如图 1(b)所示, 约有52.63%的站点SOCD呈显著增加趋势(P < 0.05).其中, 桃源站、常熟站和拉萨站SOCD增加速率相对较高, 分别为0.33、0.17和0.14 kg·(m2·a)-1.约有31.58%的站点SOCD变化趋势较小, 不足0.02 kg·(m2·a)-1, 主要分布在干旱半干旱地区, 特别是长武站和临泽站SOCD增加速率在0.01 kg·(m2·a)-1以下.

图 1(c)显示了不同土地利用方式下典型农田SOCD的变化趋势.从中可知, 不同作物类型的农田SOCD均呈波动增加的趋势.水稻田SOCD呈显著增加的趋势(P < 0.05)且增加速率最高, 为0.13 kg·(m2·a)-1.其次, 小麦田和玉米田SOCD增加速率也相对较高, 分别为0.08 kg·(m2·a)-1和0.05 kg·(m2·a)-1.在种植经济作物的农田中, 棉花田SOCD的增加速率要明显高于花生田, 二者的变化速率分别为0.04 kg·(m2·a)-1和0.01 kg·(m2·a)-1.

基于地理探测器方法, 分析不同因素(Lon、Lat和DEM; Tem和Pre; Ctyp; Styp; Sand、Silt和Clay; SpH; AN、AP和AK)对典型农田SOCD空间分异的影响程度.由图 2可知, 在因子探测中, 不同因素对SOCD空间分布的影响存在一定的差异性.其中, AN和SpH对农田SOCD空间分布的影响较大, q值均超过了0.75, 分别为0.78和0.76; 其次为Styp、Pre和Clay, 它们的q值分别为0.68、0.64和0.60, 并且这5种因素也都通过了显著性检验(P < 0.05).其他因素由于q统计量未通过显著性检验, 而不存在显著的解释能力, 这表明土壤理化性状和降水量是影响农田SOCD空间分布的重要因素.

在实际环境条件下, 农田SOCD的空间格局往往是由多种因素共同作用的结果.如表 5所示, 在交互作用探测中, 任意两个因素交互后对农田SOCD的空间分布影响均表现为协同效应, 即交互后对农田SOCD空间分布的解释能力明显高于各单因素的作用, 表现出双因子增强或非线性增强的特征.从影响农田SOCD交互作用的强度看, AN与Lat的交互作用影响最强(0.99), AN与Styp的交互作用影响次之, 为0.98.此外, AN与Pre、AN与SpH和SpH与Lat的交互作用影响也均达到了0.97.

为理解不同因素对农田SOCD空间分布的影响是否存在显著差异, 在生态探测中, 如果任意两个因素的影响具有显著差异(P < 0.05), 则标记为Y, 否则就记为N.如表 6所示, 除AN外, 其他因素两两之间的相互作用对农田SOCD的空间分布影响并未表现出一致的显著性差异, 而AN与任意一个因素结合对农田SOCD空间分布的影响具有显著差异, 这表明AN在各因素中表现得最为突出.

基于2005~2015年19个CERN典型农田生态站SOCD数据发现, 不同站点SOCD表现出明显的空间差异性, 整体表现为湿润季风区的SOCD高于干旱半干旱地区, 主要是因为干旱半干旱地区植被净初级生产力较低, 使新增的SOC输入量较少, 而土壤中较多的砂粒含量又使SOC矿化速率较高, 进而导致SOCD较低[6, 9, 30].与其他学者的研究结果对比发现, 本研究得出的典型农田SOCD空间分布特征与李金全等[9]通过文献整理方法, 得到的我国SOC空间分布格局基本相似; 但与许泉等[23]通过第二次国家土壤普查资料, 得到的我国SOCD空间分布特征存在一定的差异, 主要表现在东北地区SOCD偏低和西北地区SOCD偏高, 这主要与国家土壤普查资料在上述地区的取样点分布不均有关, 如西北地区的取样点并未涉及风沙土、灰漠土等主要土壤类型, 从而使该地区SOCD偏高.

与前期通过文献整理方法或国家土壤普查资料的结果相比, 本研究基于逐年SOCD数据发现, 不同区域SOCD的变化曲线呈波动增加的趋势, 这在之前的研究中较少被报道[31].进一步分析发现, 湿润季风区SOCD的增加速率明显高于干旱半干旱地区, 特别是华东地区SOCD呈显著增加趋势(P < 0.05).此外, 各站点SOCD也均表现出不同程度的增加趋势, 增加速率介于0~0.33 kg·(m2·a)-1之间.

水田和旱地是现阶段我国农田常见的土地利用方式.其中, 水稻田SOCD明显高于其他类型的农田, 这与通过文献整理方法[9, 31]或国家土壤普查资料[23]的分析结果相似, 反映了水田具有较强的SOC储存能力.这主要是因为水田一般长时间处于淹水状态, 相较于旱地, 其SOC矿化速率较慢, 有利于SOC的积累, 进而使SOCD较高[32].并且, 水稻田SOCD呈显著增加的趋势(P < 0.05), 增加速率明显高于其他类型的农田.

前期研究证实, 我国SOCD的空间分布特征大体与Styp的地带性分布相一致[33].由地理探测器方法可知, 包括Styp在内的一些土壤理化性状, 如AN、SpH、Clay以及降水量是影响农田SOCD空间分布格局的重要因素, 特别是AN在各因素中表现得最为突出, 反映了土壤中C元素与N元素之间存在较强的耦合关系.土壤中N素的富集可以通过影响微生物活动和SOC矿化速率等途径, 从而有利于SOC积累[30], 即土壤N素积累与SOC的关系密切.此外C∶N比作为土壤质量的重要参数, 较高的C∶N比也有利于SOC贮存[34, 35], 从而提高SOCD.SpH一般与SOC含量表现为负相关关系, 酸性土壤会降低微生物和酶的活性, 影响SOC的溶解性; 反之在中性或弱碱性土壤中, 土壤微生物活性增强, 使得SOC分解加快[36], 从而降低SOCD.Moharana等[37]基于大田试验结果也发现, 通过合理施用N肥, 并调节SpH, 可以提高SOC的积累[37].除AN与SpH外, 由于Styp本身包含气温、降水与成土母质等因素, 因而在前期研究中被认为是影响SOCD空间分布的重要因素之一[33].Clay含量直接影响土壤的通气性, 并会与SOC形成有机-无机复合体.当Clay含量较高时, 它对SOC的物理性保护作用就越强[38].在实际环境条件下, AN分别与Lat、Styp、Pre及SpH等因素交互后影响力显著增强, 反映了农田SOCD空间分布在一定程度上是AN与其他因素综合作用的结果.

近年来Meta分析研究也表明, 农田表层SOCD可能还受到耕作方式、秸秆还田等因素的影响[31], 但受限于相关数据获取以及指标定量化等因素的制约, 本研究目前未能充分考虑这些因素在农田SOCD空间分布中的具体影响.后续将整合涵盖不同地理区域的多源农田SOCD数据, 利用机器学习与地学信息图谱方法[39, 40], 定量分析自然与人为因素对农田SOCD的影响效应, 并形成长期且连续的农田SOCD时空变化图谱, 为我国农田温室气体减排能提前实现“2030年联合国可持续发展目标”, 提供重要的科学依据.

(1) 我国不同站点农田SOCD表现出明显的空间差异性, 整体表现为湿润季风区的SOCD高于干旱半干旱地区; 在不同土地利用方式下, 水稻田SOCD明显高于其他类型的农田.

(2) 2005~2015年各站点SOCD均表现出不同程度的增加趋势, 但干旱半干旱地区SOCD增加速率相对较小.水稻田SOCD呈显著增加的趋势(P < 0.05), 其增加速率高于其他类型的农田.

(3) 土壤理化性状和降水量是影响农田SOCD空间分布的重要因素, 特别是AN在各因素中表现得最为突出, 且AN与其他因素交互后的解释力显著增强.

致谢: 感谢中国生态系统研究网络提供的19个典型农田生态站土壤有机碳密度的监测数据.