基于Sentinel-2遥感与涡度协方差通量估算的农田碳输入监测方法研究及其对温室气体清单的改进意义

《Carbon Management》：Improving agricultural carbon monitoring with Sentinel-2 and eddy-covariance-based plant productivity estimates

本刊推荐：本研究创新性融合Sentinel-2卫星遥感数据与涡度协方差（EC）测量技术，构建了基于红边叶绿素指数（CIrededge）和光合有效辐射（PAR）的初级生产力（GPP）估算模型，显著提升农田土壤碳输入量的估算精度。通过揭示现行国家温室气体清单中碳预算与实测数据的差异，为整合地球观测与通量测量数据提供了关键技术路径，对优化农业碳汇评估和政策制定具有重要科学价值。

　　

摘要

本研究通过整合哨兵2号（Sentinel-2）卫星遥感数据和基于涡度协方差（EC）的植被生产力估算，提出了一种改进农业碳监测的方法。该方法利用红边叶绿素指数（CI_rededge）和光合有效辐射（PAR）构建统计回归模型，以估算每日总初级生产力（GPP），并借助五个EC站点的数据进行模型训练。留一站点交叉验证表明，模型对年GPP的预测均方根误差约为均值的10%。此外，预测累积GPP在135个农田生物量测量中解释了60%–70%的观测变异。最终，该方法应用于芬兰三个区域，通过将净初级生产力（NPP，假设为GPP的50%）与收获碳输出之差作为土壤碳输入量，结果显示与当前国家温室气体清单中的同化方法相比，小麦碳输入量差异在5%–30%之间，而大麦、绿休耕和饲草分别高出50%–100%、40%–50%和100%–160%。这些结果突显了当前清单与EC测量碳预算之间的差异，为将地球观测和通量数据整合到温室气体清单和报告系统提供了路径。

引言

可靠的土壤有机碳（SOC）储量监测与预测方法对于理解不同土地管理实践的气候影响至关重要。 deforestation、农业扩张和传统耕作方式导致土壤碳储存下降，许多农业土壤仍是二氧化碳（CO₂）的源。然而，再生农业实践的成功实施为逆转SOC损失提供了机会。许多国家和公司依赖独立的土壤碳模型来评估土壤碳储量的演变，这些模型的预测依赖于作物残茬、死根和根系分泌物等土壤碳输入的准确估算。当前，土壤碳输入常通过同化方程从收获产量中估算，但不同方程可能导致不一致的结果，从而引入SOC估算的不确定性。卫星遥感通过持续监测植物光合作用，可能为土壤碳输入提供自上而下的约束。本研究定量比较了卫星驱动方法与现有同化方法在区域尺度上估算土壤碳输入的差异。

遥感监测作物生产力的方法主要集中在估算光合作用速率（GPP）和生物量密度。地上生物量密度或产量可直接通过遥感冠层反射率与机器学习结合生物量测量来估算，但这些方法可能无法捕捉同化碳的完整周转。更全面的方法是将遥感衍生的生物物理变量（如LAI）同化到过程模型中，但这些模型需要详细的土壤和管理信息，且计算量大，可扩展性受限。相比之下，使用遥感植被指数直接或通过光利用效率（LUE）模型估算植物光合作用是一种简单方法。在LUE模型中，GPP是光合有效辐射（PAR）、吸收PAR比例（fAPAR）和LUE的乘积，其中fAPAR通常通过植被指数估算，LUE则基于生物群落类型和环境因素（如温度和蒸汽压赤字）估算。校准植被指数以确定fAPAR和生物群落特定LUE是一项挑战。通过利用fAPAR、LUE与冠层叶绿素含量之间的紧密联系，可将GPP建模为PAR和叶绿素敏感植被指数的函数，该方法在本研究中被采用。

卫星遥感衍生的GPP已用于通过累积生长季GPP来估算产量，使用MODIS和Landsat数据。高分辨率的Sentinel-2数据能够在更细空间尺度（如单个田块）上估算GPP，但目前尚不清楚这些方法预测的土壤碳输入与传统同化估算的比较。本研究旨在从Sentinel-2多光谱仪器（MSI）测量的田块尺度植物生产力出发，估算土壤碳输入。基于早期工作，我们开发了一个非线性回归模型来描述每日光合作用，使用CI_rededge和PAR作为解释变量，并通过EC站点的碳通量估算进行校准。通过插值缺失卫星数据的CI_rededge值，我们能够预测季节和年尺度的GPP，并通过与EC和生物量数据的比较验证估算。最后，我们将卫星衍生的植物生产力估算与作物产量统计结合，估算芬兰三个区域的土壤碳输入，并与国家温室气体清单中使用的产量同化方法进行比较。

我们的方法将土壤碳输入确定为净碳同化与收获产量之差，这使得估算对碳分配变化的敏感性低于典型同化方法。这可能提高准确性，特别是对于地下碳输入（如根系周转），这些通常发现仅受收获产量的弱约束。通过提高土壤碳输入估算的稳健性，我们的方法有助于改进碳模型，并为增强碳汇的农业政策提供信息。

材料与方法

总体设计

本研究总体设计概述如图1所示。主要步骤包括预处理芬兰五个站点的EC通量测量、拟合回归模型，以及使用留一站点交叉验证和与40个额外采样站点的生物量测量比较来评估模型结果。EC和生物量采样站点的位置如图2所示。最后，我们通过在区域水平上评估从国家土地田块识别系统中随机选择的作物田块的净初级生产力（NPP）来应用该方法。区域NPP估算进一步与基于调查的产量统计结合，以估算区域尺度的土壤碳输入。

我们的方法强调不确定性的 consistent 传播。我们使用贝叶斯技术拟合回归模型和卫星数据的时间插值，并通过采样后验概率分布来量化预测不确定性。不确定性估算考虑了模型训练数据和卫星时间序列中的数据差距带来的不确定性。贝叶斯混合效应模型用于估算预测和采样误差对区域GPP和碳输入估算的联合影响。

通量测量

涡度协方差测量与数据处理

使用EC技术在芬兰五个不同农田站点测量净CO₂通量，站点列表如表1所示。四个站点位于芬兰南部的Uusimaa（SMEAR-Agri Haltiala和SMEAR-Agri Viikki）、西南芬兰（Qvidja）和Kanta-H?me（Hauho）地区，一个站点位于北部的北博滕区（NorPeat, Ruukki）。Ruukki为泥炭土，Haltiala、Viikki和Qvidja为细质地粘土，Hauho站点为更粗质地的粉壤土。测量期在2018年5月至2023年12月之间，各站点从1.5年到5.5年不等。测量年间，Qvidja种植饲草（提摩西草、草地羊茅和三叶草），Hauho种植燕麦，其余三个站点为混合轮作，包括饲草和谷物。部分数据集已发布（Qvidja和Ruukki），本研究重新计算以使数据后处理在所有站点间一致，重新计算的通量与之前发布的数据集无显著变化。

EC仪器在所有站点相似。CO₂混合比使用封闭路径CO₂/H₂O气体分析仪（LI?7200）测量，三维风使用超声风速计（u-Sonic3 Scientific）记录，采样率为10 Hz。测量高度从Qvidja的2.3 m到Ruukki的3.3 m不等。样本空气通过加热管道以12 L min^?1的流速传输至气体分析仪。辅助气象测量包括气温和相对湿度、土壤温度、土壤湿度和光合有效辐射，数据以30分钟平均值存储。

半小时湍流CO₂通量使用Eddypro软件计算，遵循标准质量控制和处理协议。首先，对10 Hz原始数据进行筛选，包括去尖峰、振幅分辨率和绝对限制，之后允许30分钟原始数据最多缺失1%，否则丢弃该时段。接下来，应用坐标系统双旋转以沿当地风流线对齐超声风速计，并块平均每个30分钟时段以去除线性趋势。气体分析仪和超声风速计信号之间的时间滞后通过自动时间滞后优化确定。水蒸气和温度波动对LI-7200测量的影响在分析仪内部和足够长的采样线中考虑。块平均和垂直风速与CO₂混合比同谱中最高频率衰减引起的频谱损失通过应用方法校正。30分钟通量通过相对平稳性和积分湍流特性筛选，如果这些测试的组合标志指示最低质量，则丢弃通量。弱湍流时期通过应用摩擦速度（u^*）阈值去除，该阈值通过移动点转换方法估算。最后，进行源区分析以过滤不代表目标区域的通量数据，通过应用微气象足迹模型计算从EC塔到田间边缘的跨风积分足迹分布，如果累积足迹小于70%，则丢弃数据。各站点的数据覆盖率见表S3。

间隙填充与分区

数据使用神经网络深度集成进行间隙填充。我们使用五模型集成，每个基础模型为具有两个隐藏层（50个节点）和非线性ReLU激活的顺序模型。每个模型用于预测拉普拉斯分布的平均值和尺度参数对数，以考虑EC数据中的异方差噪声。拉普拉斯负对数似然用作损失函数。权重使用HeNormal初始化随机初始化，使用Adam优化器进行训练，10%数据用作验证集用于早停，耐心为20轮。预测因子包括PAR、气温、土壤温度和土壤湿度、蒸汽压赤字、距上次收获的天数以及描述一天时间和季节的两个循环函数。实现使用TensorFlow库。

当估算CO₂通量或不稳定生态系统交换（NEE）的不确定性时，我们考虑了与测量、间隙填充和u过滤相关的不确定性。测量不确定性通过通量大小与残差标准偏差之间的线性关系从模型残差估算。间隙填充不确定性通过结合模型预测的偶然和认知不确定性估算，这些使用拉普拉斯分布的预测均值和尺度以及模型预测的方差计算。与u阈值相关的不确定性通过使用不同阈值过滤数据、填充替代数据并计算NEE的方差传播到NEE中。

间隙填充数据被划分为总生态系统呼吸（TER）和GPP。TER假设遵循Arrhenius型模型，其中TER = R₀ × e^{E(1/T₀ ? 1/(T?T₁))}，R₀为283.15 K时的生态系统呼吸，E为呼吸的温度敏感性，T₀ = 56.02 K，T₁ = 227.13 K，T为实测气温。呼吸的温度敏感性从夜间（PAR < 20 mol m^?2 s^?1）数据使用15天移动窗口（每拟合至少20个观测且温度范围至少5 K）确定。E作为成功拟合（标准误差<50%）的中位数确定。R₀使用7天移动窗口拟合，同时保持E固定。在推导两个参数估算后，TER使用气温建模并从NEE中减去以获得GPP。NEE的不确定性通过使用不同u^*阈值划分NEE数据和间隙填充集成的个体模型预测传播到GPP。仅使用最多70%间隙填充数据的每日GPP值用于拟合S2-based GPP模型。

模型输入数据

植被指数

开发的GPP模型基于CI_rededge植被指数，该指数对叶绿素含量敏感，在先前使用LUE-based模型估算GPP的研究中表现良好，并在测试六个植被指数时产生最佳拟合模型（补充图S2）。

CI_rededge源自Copernicus Sentinel-2卫星数据。Sentinel-2由两颗卫星组成，在芬兰的重访时间为两到三天。我们使用Google Earth Engine（GEE）提供的底层大气反射率产品（level-2A）。包含在level-2A产品中的场景分类数据用于确定图像获取日期，当田间无云、云阴影、雪或水覆盖时。每个田块的反射率数据以10 m空间分辨率检索。对于区域级分析，我们从芬兰食品局提供的地理空间援助应用（GSAA）数据中检索田块边界。否则，使用已知田块边界。对于EC站点，卫星图像提取与典型EC足迹一致的目标区域。CI_rededge植被指数使用Sentinel-2反射波段B5（705 nm）和B7（783 nm）计算，作为田间所有像素的平均值。

太阳辐射

为估算GPP，植被指数与每日PAR通量结合。辐射通量从Copernicus大气监测服务（CAMS）全天辐射网格数据产品获得，该产品由结合非洲和欧洲0.1度地理网格上卫星产品衍生的辐射通量组成。为评估每日PAR，我们提取每个田块质心的15分钟全球辐照度时间序列，聚合到每日时间分辨率，并将值乘以PAR分数0.5。

建模初级生产力

我们假设GPP的平均响应μ_GPP是植被指数CI_rededge和光合有效辐射的函数，由μ_GPP = (a × CI_rededge × PAR) / (b + CI_rededge × PAR)给出，其中a定义为全局分量a_glob和站点效应a_i的乘积，b为全局参数。这种层次结构允许我们捕捉GPP的全局模式，同时考虑站点特定差异。由于数据限制，未纳入作物特定效应，尽管研究站点种植了不同作物，且三个站点的作物类型逐年变化。作物类型对GPP的可能影响在结果部分讨论。

为考虑不确定性，GPP的标准偏差σ定义为σ = √(σ_EC² + σ_model²)，其中σ_EC代表EC-derived GPP的不确定性，σ_model代表回归模型未解决的每日变异性的残差误差。观测GPP值GPP_obs假设服从均值为μ_GPP、标准偏差为σ的正态分布。使用贝叶斯推断，站点间观测变异性也纳入模型，允许将该变异性传播到预测的不确定性估算中。使用马尔可夫链蒙特卡罗（MCMC）采样估计全局和站点特定参数的后验分布。使用numpyro包实现的No U-Turn Sampler。模型参数的先验总结在表2中。

CI_rededge的时间插值

尽管Sentinel-2A和2B重访相对频繁，但由于云污染，生长季数据缺口常见，植被指数需要在时间上插值以评估每日GPP。我们使用高斯过程（GP）回归插值植被指数，这是一种非参数贝叶斯技术，其中未观测时间的值从以观测值为条件的多元正态分布中采样。高斯过程技术的一个关键优势是其内置的不确定性量化。预测不确定性自动校准到输入数据，且当两个观测之间的时间间隔增加时，局部预测误差增大。

GP给出的插值由均值函数和协方差核控制。均值函数规定数据中的确定性分量，这里设置为零。因此，插值趋向于零（随着方差增加）当时间到最近观测增加时。这对于在北方条件下插值植被指数是可接受的，因为无数据的延长时期通常与雪或裸土条件相关。

协方差核参数化GP的随机分量。由于植被指数通常随时间平滑变化，我们使用带有对角线项的平方指数核，通常用于模拟受噪声污染的光滑信号。核依赖于三个超参数，描述时间相关分量的幅度、不相关分量（噪声）的幅度和相关时间尺度。超参数为每个站点和年份单独估计。我们使用tinygp包实现高斯过程回归，并使用numpyro包进行超参数推断。技术细节见S1.1。

为减少冬季插值过程中的不确定性，当田间被雪覆盖时，我们将植被指数设为零。这是基于S2数据提供的场景分类层完成的。然而，由于分类偶尔在夏季指示雪覆盖，实施了质量控制步骤以确定雪覆盖的合理性。首先，排除雪覆盖时间插值植被指数，并采样GP的100个实现。然后，对于每个指示雪覆盖的观测，执行独立双面t检验，针对植被指数为零的零假设。p值低于0.01的雪分类被拒绝为不可能。最后，植被指数插值，并在合理雪分类处插入零值。

平滑变化植被指数的假设在草地收获时被违反。为适应这一点，我们实现了修改的平方指数核，描述具有有限数量不连续性的分段平滑函数。通过此修改，插值仅由同一间隔内的观测通知，该间隔由收获划定。收获日期通常可用于实验站点，但对于随机采样的田块不可用，对于这些站点，不连续性的存在作为推断过程的一部分估算，如Sect. S1.1.3所述。

交叉验证

为评估模型在不同站点间的泛化能力，我们进行了留一站点交叉验证。这里，一个站点的数据被排除并用作验证集，而模型使用剩余站点的数据拟合。此过程重复五次，使得每个站点的数据被用作验证集一次。

在每次迭代中，排除站点的站点特定缩放效应（a_i）不可用，预测通过从后验预测分布中采样站点特定效应进行，使用尺度参数的后验分布。因此，对于留出站点，预测的总不确定性包括站点特定缩放效应a_i的变异性贡献，以及全局参数（a_glob和b_glob）和残差不确定性（σ_model）的变异性。

模型性能使用均方根误差（RMSE）、平均偏差（MB）和决定系数（R²）评估，使用EC-derived GPP值和排除站点的预测值计算。此外，我们通过比较观测和预测概率评估预测不确定性的校准，通过计算落在指定预测区间内的观测数据点比例完成。这里，我们将EC-derived GPP的不确定性添加到预测不确定性中，以考虑观测误差。

生物量测量

为补充使用通量测量的留一站点交叉验证，我们测试了预测累积GPP与40个采样站点测量的地上生物量之间的关系。虽然我们的方法不直接估算植物生物量，但地上生物量的增长不能超过总碳同化，因此我们期望作物生物量与累积GPP之间存在单调关系。

生物量测量从Carbon Action数据集中提取，该数据集包含参与碳农业实验的20个芬兰农场的数据。数据集包括约2公顷的年作物和草田块，位于常规和有机管理农场。对于每个站点，数据包括一个对照田块和一个处理田块，测试碳农业实践；本研究中未区分这两个田块。生物量通过从每个田块切割4-6个生物量样方（33×33 cm）测量，围绕三个GPS定位采样点。由于田间定期耕作和重新播种，生物量的空间变异性低，样方大小被认为足以捕捉田块的代表性估算。生物量切割至土壤表面1 cm。湿生物量在现场7月测量，并取生物量样本进行烘箱干燥以确定干物质含量。

对于每个生物量测量，我们评估了相应田块的每日GPP之和，要么从年初开始（对于收获前的年作物和草地），要么从上次收获或割草开始（对于收获后的草地）。放牧田块由于量化放牧牲畜引起的生物量变化的不确定性而被省略。我们还省略了由留田的茬或割草植物生物量组成的测量。遵循这些标准，从原始182个测量中保留了135个。

区域GPP和土壤碳输入

使用卫星based GPP方法，我们估算了2020–2023年间芬兰三个行政区域（县）为小麦、大麦、饲草和绿休耕管理的农田的平均年GPP。估算的GPP进一步转换为NPP，并与区域报告的收获产量结合，以估算年土壤碳输入，即NPP中未以收获输出的部分。选择四种作物以覆盖一系列对比管理实践，并包括对各区域经济重要的作物。三个区域（西南芬兰、南博滕和北萨沃）气候条件不同，包括主要谷物种植区（西南芬兰）、草地主导区（北萨沃）和混合土地利用区（南博滕）。

区域GPP通过调查采样估算。对于每年、区域和作物，我们评估了从国家土地田块识别系统中采样的n=100个农业田块的GPP。田块按大小（表面积）比例概率抽样，无放回。如讨论所示，这种策略减少了在种群内估算面积加权平均的不确定性，这是评估区域上 aggregate 碳通量所需的量。田块每年重新采样，因此总共处理4800个站点年以评估区域GPP。

区域GPP估算的不确定性是采样误差和GPP模型中随机和系统误差的函数。模型不确定性通过蒙特卡罗采样估算。对于每个站点，我们采样100个高斯过程实现以确定馈送到GPP模型的植被指数的不确定性，同时，我们采样回归模型的后验参数分布以包括参数不确定性的影响。站点特定模型参数为每个站点重新采样，而站点独立模型参数仅采样一次，并统一应用于给定作物和区域的所有站点和年份。

然后，通过将贝叶斯混合效应模型拟合到样本并遵循概述的调查方法，估算全种群中的区域平均GPP。这种方法分两步进行：首先，我们估计一组参数，定义种群内GPP的分布。然后，我们使用第一步的参数估计在全种群田块上采样GPP。第一步中估计的种群参数的不确定性传播到第二步。对采样GPP的初步检查显示田块大小与GPP之间存在弱但显著的相关性（图S4），因此我们的模型包括田块大小效应的校正项。

混合效应模型假设每个采样田块的年GPP可分解为加性项：y_i = x_i + β(s_i ? s?) + ε + δ_i，其中y_i为第i个采样田块的观测GPP，x_i为无田块大小效应时的田块GPP，δ_i和ε为站点特定和站点独立误差。观测值y_i取等于回归模型后验预测分布的平均值。项β(s_i ? s?)允许田块大小s_i（给定为与样本平均田块面积s?的偏差）与田块GPP之间的线性关联。我们假设未调整GPP x_i正态分布，均值为μ，方差为σ²，这些与斜率β一起估计。

模型使用MCMC方法求解（类似于建模初级生产力），假设先验分布列于Sect. S1.2。这导致μ、σ²和β的后验样本，进而用于采样未观测田块的后验预测分布。然后，每个MCMC样本中的值在空间上聚合，以评估全种群田块上面积加权平均的后验分布。如果由β参数化的田块大小效应显著，则该值不同于x_i的平均值。然而，由于田块按大小比例概率抽样，差异最小化，因为与样本平均田块大小的偏差平均而言对加权平均贡献最大的田块较小。

土壤碳输入

农田土壤碳输入主要包括收获残茬、茬、凋落物、死根和根系分泌物，以及可能的外源输入如进口粪肥或堆肥。外源C输入本研究未考虑，由于大田作物缺乏长寿木质生物量库，年水平上的碳输入几乎等于净初级生产（GPP减去自养呼吸）与收获生物量中去除的碳之差。

我们使用区域GPP估算通过将GPP缩放为NPP并减去芬兰自然资源研究所统计数据库报告的收获产量来估算区域尺度的土壤C输入：F_C = r_NPP × GPP ? f_CY_DM，其中F_C表示碳输入，f_C为植物生物量中的碳质量比例，Y_DM为干物质作物产量，r_NPP为净与总初级生产之比。芬兰的收获残茬通常留田或用于动物垫料，遵循研究，我们仅考虑谷物和饲草收获从生态系统中输出。

作物产量统计基于年度调查，并与我们的GPP在同一区域水平给出；此聚合水平也用于芬兰国家排放清单中的CO₂排放计算。报告的作物产量转换为干物质，假设谷物和干草水分含量14%，青贮水分含量66%。所有作物的碳含量f_C为0.45 g C g DM^-1。

不同生态系统中NPP:GPP比的估算通常范围在0.2至0.8之间，但中心约0.5，这里采纳为r_NPP的中心估计。我们模拟比r_NPP的不确定性为 normally 分布随机变量，标准偏差0.1，与研究一致，采样r_NPP作为站点独立随机效应，并按照前节描述GPP的步骤进行混合效应模型。

我们将我们的C输入估算与当前用于芬兰国家温室气体清单的方法估算进行比较，该方法应用同化系数和规定的根和茎生物量组合，基于报告的作物产量评估年C输入。类似于方法，该方法将土壤C输入划分为地上残茬、根生物量和根沉积。地上残茬从产量使用收获指数计算，除绿休耕外，对其使用规定的生物量输入。根生物量使用谷类和其他年作物的茎:根比计算，对于草作物，使用规定的根生物量。根沉积与根生物量成比例计算，使用固定根周转率。我们使用列出的系数计算这些同化based C输入，用于与我们Sentinel-2 based估算相同的区域和年份。

结果

模型拟合和留一交叉验证

EC-derived GPP与CI_rededge × PAR之间的关系，以及拟合模型及其预测不确定性，如图3所示。所有站点的GPP随CI_rededge × PAR增加而增加，但对于较高CI_rededge × PAR值，速率开始 plateau。在较高辐射水平和增加植被覆盖下，不同站点的GPP估算变异性更大。饲草和谷物的GPP似乎有些 distinct，但差异也可能归因于站点效应（补充图S5），因为Qvidja的EC-derived GPP估算始终低于使用所有站点数据拟合模型的S2-based估算。然而，大多数每日GPP总和落在预测不确定性范围内。

留一站点交叉验证结果表明，GPP模型的平均RMSE在日尺度上为1 g C m^?2 d^?1，月尺度上为31 g C m^?2 mo^?1（图4），站点间变异小（1.2–1.7 g C m^?2 d^?1和28–34 g C m^?2 mo^?1）。EC-derived