利用先进的遥感技术估算氮素胁迫及其对作物产量的影响：以苏丹盖齐拉灌溉区为例

《International Journal of Applied Earth Observation and Geoinformation》：Estimating nitrogen stress and its impact on crop yield using advanced remote sensing approach: A case study of Gezira irrigation scheme, Sudan

农业遥感与机器学习结合估算作物氮含量及产量的方法研究及验证。

　　在当今世界，确保粮食安全已经成为一项至关重要的任务。随着人口增长和农业需求的增加，对粮食和淡水资源的需求持续上升，同时，气候变化的不确定性也给农业带来了更多挑战。粮食安全与水资源之间的相互依赖性已经被联合国粮食及农业组织（FAO）所强调，因为在一个领域采取的措施可能会对另一个领域产生深远的影响。因此，提升水资源和土地的利用效率对于改善人类福祉、减少贫困以及促进可持续发展具有重要意义。为了实现这一目标，持续提高水和土地的生产率是关键，这不仅对当前世代至关重要，也对未来世代具有深远影响。

农业作为人类生存的基础，对经济结构和国家发展有着深远的影响。因此，对作物在其生长周期中的状态进行细致监测变得尤为重要。这种监测需要考虑多个因素，包括作物的营养状况、土壤水分含量、植物活力以及受到生物和非生物因素影响的植物压力。在作物生长的不同阶段，如果某些关键参数偏离了最佳值，可能会对作物的发育产生负面影响，从而导致产量和生产效率的下降。在这一背景下，氮元素的作用尤为关键。氮是植物、动物和土壤生态系统中不可或缺的成分，其在土壤中的存在形式多种多样，且经常在不同形式之间进行转化。当氮含量不足时，会对植物生理学产生显著影响，包括叶片叶绿素含量和光合作用速率，进而影响植物的发育和最终的生物量产出。此外，氮含量不足导致的光合作用速率下降主要源于气孔导度的减少，而不是叶片结构内部碳固定能力的下降。

精准农业（PA）因其能够提升农民对田间作物和土壤状况的理解而受到高度重视，它在提高生产效率和优化投入方面展现出巨大潜力。随着精准农业的发展，变量率技术（VRT）逐渐成为其中的重要组成部分，它通过利用特定区域的数据，使农民能够以不同的速率在田间精准施用投入品，如肥料和水。进一步发展，变量率施肥（VRF）作为一种先进的精准农业技术，结合了嵌入式高速计算机、遥感（RS）技术、地理信息系统（GIS）、精确全球定位系统（GPS）接收器、土壤图以及电子传感器和执行器的实时作物特性测量和预测能力。VRF方法有望提高肥料使用效率，减少对环境的影响，并改善作物质量和产量。

遥感技术在精准农业和变量率施肥中的发展起到了关键作用。遥感产品通过使用先进的传感器和分析工具，为农业决策支持系统提供了重要的信息。近年来，基于遥感的精准农业技术已被集成到商业农业实践中。将变量率施肥与处方图相结合，被认为是实现精准施肥的一种有前景的策略，从而提高氮的施用效率。为了生成变量率图，可以通过使用各种主题图层作为输入，结合经典的或无监督的聚类算法，从而划分出具有特定管理需求的区域。

遥感技术提供了一种灵活的方法来估算各种参数，这些参数能够反映作物的氮素状况。在植物中，氮素含量通常与叶绿素（Chl）相关的波长表现出一定的相关性。值得注意的是，随着水分含量的减少，反射率会增加，而水分含量的关键信息通常编码在特定波长范围内，如1450 nm、1940 nm和2500 nm。在新鲜叶片的反射率光谱中，由于水分的显著吸收，可能掩盖了与植物生物化学成分（如氮素吸收特征）相关的较弱吸收信号。

为了模拟电磁辐射在大气和植被冠层中的传播过程，已经开发了多种模型，这些模型考虑了作物的结构和内在特性，以准确预测辐射的传输情况。辐射传输模型（RTM）是一种基于物理原理的模型，已被广泛应用于氮素的估算。RTMs能够解释光子与植物生物物理和生物化学属性之间的相互作用，并在遥感领域有广泛的应用，包括从卫星波段反演参数、开发参数化和非参数化回归技术，以及将机器学习（ML）回归方法与RTM模拟相结合的混合回归方法。

尽管传统的估算方法，如使用NDVI等指数和现场作物生物物理数据开发的回归模型，已被广泛应用于氮素估算，但这些方法在不同时间和空间背景下的可迁移性往往受到限制。近期的研究突显了这些挑战，并强调了需要更稳健的替代方案。因此，本研究聚焦于一种新颖的方法，该方法利用基于RTM的物理模型，这些模型在估算作物生物物理和生物化学参数方面已被证明更加稳健、可迁移且准确。这种方法的主要局限性包括“不适定问题”，即不同的生物物理和生物化学参数组合可能产生相似的反射率值，以及与其它方法相比计算速度较慢的问题。

本研究旨在利用遥感数据和RTMs来评估整个种植季节中农田的氮素浓度。这种技术的结合提供了一种新颖的方法，用于估算与土壤肥力和作物营养相关的关键变量。研究的具体目标包括利用卫星数据捕捉关键变量的潜力、通过训练有监督的机器学习模型评估选定作物的氮素水平和生物量估算的准确性、评估氮素水平对作物最终产量的影响，以及测试所生成模型的时空可迁移性。通过实现这些目标，本研究希望为遥感和机器学习在精准农业中的应用提供有价值的见解，弥合这些技术在优化作物管理实践方面的现有知识空白。鉴于盖兹拉计划作为全球最大的灌溉项目之一，其历史贡献包括自1980年至1990年期间占苏丹总小麦产量的60%，本研究对于推动大规模灌溉系统中的可持续实践具有特别重要的意义。

为了测试训练模型在不同时间和空间环境下的可迁移性，本研究在2022-23年冬季进行了评估。评估涉及在盖兹拉计划中的一个不同子区——瓦德希尔尔区域（Wad-Hilal）的八个三级单位（Nimras）中运行训练好的GPR模型，以预测氮素浓度和fAPAR。随后，计算了预测氮素浓度与模型和实际记录的产量之间的相关性。

在本研究中，使用了Sentinel-2（S2）数据进行分析。S2任务以其频繁的重访周期（每5天一次）、精细的空间分辨率（10米）、广泛的全球覆盖以及开放获取政策而著称。这些特性使S2成为开发近实时农业服务的重要资源。本研究中，使用了S2多光谱仪器（MSI）的L2A产品，以获取特定作物季节的S2任务瓷砖数据。L2A数据通过Sen2Cor处理器进行了大气校正，从而消除了额外校正步骤的需要，并提高了反射率值的可靠性。所有图像的云覆盖比例均低于20%的阈值。随后，获取的L2A产品被裁剪以匹配研究区域，并通过QGIS中的半自动分类插件将数据从数字数转换为10米分辨率的反射率值。这些10个波段将作为输入，用于预测氮素含量和fAPAR。这些波段的选择基于其与植物氮素含量和生物量估算的强相关性。值得注意的是，已有研究强调了某些波段在这些估算中的重要性。例如，Guerif等人（2007）展示了利用可见光和红边波段的光谱特征评估氮素浓度的可行性，这些波段能够反映叶绿素含量。

同时，Berger等人（2020）主张利用短波红外（SWIR）光谱域，通过蛋白质作为氮素的代理来估算氮素水平。结合敏感性分析的结果和现有文献的见解，确定了用于估算氮素浓度和生物量的最佳波段。

在本研究中，采用了Python版本的PROSAIL-PRO模型，生成一个查找表（LUT），作为训练数据库。这一耦合模型能够模拟冠层尺度的反射率，其预测结果基于一组多样化的生物物理参数（如叶面积指数（LAI）和平均叶倾角）以及叶片生物化学输入（包括叶绿素-a+b含量、等效水厚和叶黄质含量）。氮素含量（N）的直接计算从叶蛋白含量（Cp）进行，利用蛋白质到氮素的转换因子4.43，该因子由Yeoh和Wee（1994）在研究了90种植物后确定。这一转换因子的使用基于以下公式：N = Cp / 4.43（g/cm2）。这一过程使我们能够从叶片水平推导出冠层水平的氮素含量。

在从叶片水平过渡到冠层尺度的过程中，LAI作为尺度参数发挥了重要作用。通过使用LUT中的输入点，结合LAI和叶蛋白含量，可以计算出冠层水平的氮素含量。这种方法推荐使用Berger等人（2020）提出的方法，公式如下：冠层氮素含量（g/m2）= N（g/cm2） × LAI（m2/m2） × 10,000。

在冠层辐射传输领域，fAPAR的计算涉及多个因素。这包括冠层对辐射的直接吸收（a1）以及背景反射后被植被吸收的分数（a2）。因此，绿色冠层的总能量吸收可以表示为：a(λ) = a1(λ) + a2(λ)。这里，λ表示波长。a1的计算基于背景反射率（rg）可以忽略不计的假设。在此背景下，瞬时fAPAR可以表示为：fAPAR = ∫0.4 μm到0.7 μm a(λ) dλ。这一概念框架借鉴自Fan等人（2014）的研究。在这一框架中，计算fAPAR的关键参数，无论是直接计算还是作为SAIL模型的次级量，都至关重要。SAIL模型中的参数skyl，即漫反射的太阳辐射分数，考虑了太阳天顶角（SZA）（Fran?ois等人，2002），该参数基于中纬度地区的平均大气条件。skyl的计算公式如下：skyl = 0.847 - 1.61 sin(90° - SZA) + 1.04 sin2(90° - SZA)。

模拟过程从计算直接到达冠层的太阳辐射（edir）和漫反射的辐射（edif）开始，其中es和ed分别表示直接和漫反射太阳辐射的光谱。随后，利用4SAIL生成的参数和参数dn（表示与背景土壤的相互作用），计算了αs和αd，即冠层对太阳和漫反射辐射的吸收率。最终，通过应用公式计算了fAPAR：fAPAR = ∑(αs × edir[i] + αd × edif[i]) / ∑(edir[i] + edif[i])，其中i从0到300。该切片因子定义在0到300之间，用于计算这些值，同时仅考虑400到700 nm之间的波长。

为了评估作物产量，使用了以下公式：NPP = Sc × Rs × fAPAR × SM × ?lue × ?p × ?t × ?co2 × ?AR。其中，Sc表示从干物质生产（DMP）到NPP的缩放因子（0.045），Rs表示总短波入射辐射（GJ/ha/day），来自ERA5-Land Hourly的“surface_solar_radiation_downwards”。?p表示气候效率（0.48），SM表示土壤水分压力减少因子（等于1），?lue表示最佳条件下的光能利用效率（对于C3作物如小麦为2.7 kgDM/GJPA），?AR表示自动呼吸后的保留分数（0.5），而?co2表示归一化的二氧化碳效应。二氧化碳效应的估计使用了Veroustraete（1994）提出的函数，其中二氧化碳混合比（CO2）为参考年份1833年的281 ppm，而2021年和2022年的CO2值分别为414.7 ppm和417.07 ppm。

随后，使用简单的线性插值估算NPP在每月开始和结束时的累积值，以确定整个生长季节结束时的总NPP。接着，利用Mul和Bastiaanssen（2019）提出的公式估算产量：产量 = (NPP × (AOT × 22.222 × fc)) / (1 - Mc)。其中，AOT表示上述总生物量（0.85），fc表示光能利用效率的相关因子，Mc表示新鲜生物量中的水分含量（0.15用于小麦），HI表示收获指数（取自盖兹拉地区的水利研究机构，设为0.4）。

在模型优化过程中，我们采用了ADAM优化器，并将其学习率设置为0.001。ADAM优化器因其在训练过程中动态调整每个参数的学习率而受到青睐。Kingma和Ba（2014）详细阐述了ADAM优化器在处理高维数据和参数空间时的计算效率、低内存占用、对梯度对角线缩放的鲁棒性以及其适应性。在优化过程中，我们定义了一个负对数似然（NLL）损失函数，该函数包含了高斯过程（GP）的关键要素，包括均值函数、核函数、索引点（在此情境中代表反射率）和观测噪声方差。NLL损失函数旨在最小化模型的负对数似然，从而优化模型以做出精确预测，同时考虑到观测值中的固有噪声。

这些优化步骤共同促进了高斯过程回归（GPR）模型在后续分析中的性能和可靠性。随后，模型进入训练循环，其中计算并存储了批量NLL损失值（batch_nlls），同时定期评估了完整的NLL值以全面评估模型性能（full_ll）。训练过程持续了10,000次迭代（nb_iterations），并在每次20次迭代时记录更新（log_interval）。

为了保存训练好的GP模型，我们存储了模型的训练参数，包括核函数的超参数和观测噪声方差，以及包含反射率值和相应参数值的训练数据集。该数据集有助于在对新数据进行预测时重新计算均值函数和估计标准差（STD）。

在进行后训练预测时，我们使用了tfp.Distribution GPR类。GPR模型使用训练好的核函数超参数、观测噪声方差和特定于每个模型的训练数据进行初始化。随后，这些GPR模型用于预测每个在研究区域中获取的Sentinel-2（S2）波段对应参数（N和fAPAR）及其相关的标准差（不确定性）。这些行级预测随后被准确地映射到预测地图（Predicted_parameter_map和Prediction_stdev_map）中的相应位置。

在结果和讨论部分，我们分析了PROSAIL-PRO模型的输入数据集。这些数据集由21,000个输入点组成，每个点包含14个生物物理和生物化学变量（BV）值，代表PROSPECT和SAIL参数。在这些输入点中，我们观察到LAI与叶结构参数（n）、叶绿素含量（Cab）、叶黄质含量（Car）、叶水含量（Cw）、叶干物质（Cm）和土壤亮度（rs）之间存在清晰的线性关系。LAI的分布遵循正态分布，其中大部分值集中在LAI=2附近，标准差为2。LAI的范围设定在1到8之间，因为小麦的预期LAI趋势从0到4.5。这一范围的设定有助于限制在低LAI值下进行冠层反射率模拟的潜在问题，从而影响训练好的GPR模型在作物生长早期阶段估算N和fAPAR的能力。同时，较高的LAI值对于GPR模型的训练也很重要，以避免模型在训练过程中被高值饱和，并在预测极端情况时表现更佳。

在这些数据分布中，我们观察到了一些异常值，如Cm均匀分布的初始下降，以及n和rs的分布趋势与LAI相关。这些异常值的出现可以归因于与高LAI值相关的共分布约束。同样，Cab的高斯分布的下降趋势也源于这些约束。

在生成输入点后，我们使用了第2.4.1节中的公式计算了N的值。这导致了N值的范围从2.26到45.15 g/m2。图5展示了这些N值的分布。在该直方图中，每个条形代表一个特定的N值范围，而条形的高度表示每个范围内出现的频率。最高的条形对应于最频繁的N值，这提供了数据分布的清晰视觉表示。

总体来看，N值主要集中在5到20 g/m2之间，这与小麦生长季节中叶片氮素含量的正常变化范围（2.4到13.2 g/m2）相一致。此外，N值的分布符合正态分布，几乎以均值为中心并均匀分布，这为GPR模型的训练提供了良好的基础。

在本研究中，我们还评估了模型的时空可迁移性。在2022-23年冬季，我们对瓦德希尔尔子区的八个三级单位进行了分析，这些单位覆盖了214个小麦种植农场。所有图像的云覆盖比例均低于20%的阈值。随后，获取的L2A产品被裁剪以匹配研究区域，并通过QGIS中的半自动分类插件将数据从数字数转换为10米分辨率的反射率值。这些10个波段将作为输入，用于预测氮素含量和fAPAR。这些波段的选择基于其与植物氮素含量和生物量估算的强相关性。

在模型训练过程中，我们采用了ADAM优化器，并将其学习率设置为0.001。ADAM优化器因其在训练过程中动态调整每个参数的学习率而受到青睐。在优化过程中，我们对数据集进行了划分，其中14,000个数据点用于训练，7,000个数据点用于验证。模型选择过程关注于通过验证集评估模型的性能，从而减少过拟合的风险，并确保模型对未知数据的稳健泛化能力。

N GPR模型在训练过程中表现出最强的性能。这可能是由于N值的正态分布（如图5所示），其几乎以均值为中心并均匀分布。相比之下，fAPAR值的分布（如图7所示）是右偏的。这表明，N GPR模型在预测过程中表现更为稳健，而fAPAR GPR模型在训练过程中可能出现过拟合。

在模型训练过程中，我们观察到负NLL值的出现，这可能表明模型在训练数据上出现了过拟合。负NLL值与概率密度函数中的似然值相关。这一现象在fAPAR（如图7所示）中尤为明显，因为fAPAR值通常在0到1之间。为了评估这一现象，我们考虑了训练开始和结束时的相对差异。交叉验证利用验证数据点来确认这一差异的重要性。

N GPR模型和fAPAR GPR模型的RMSE和R2值分别为3.56 g/m2和0.66，以及0.044和0.82。这些指标表明，所有模型都没有过拟合数据，并在各自的任务中表现出相对较强的性能。

通过图12中的条形图，我们可以比较实际和预测的产量值。预测的产量值在三个农场中分别为3.86吨/公顷、3.67吨/公顷和3.91吨/公顷。产量估算的RMSE为0.54吨/公顷，这在使用简单的线性插值进行累积计算的背景下可以被认为是一个良好的结果。这些估算进一步验证了fAPAR-GPR模型在预测作物产量方面的性能。条形图显示， Farms 1和2的产量被高估，而Farm-3的产量则被低估。由于Farm-2的产量与实际值之间的差距最大，可以推测这些差距可能与Farm-2边缘区域的高不确定性（STD）值有关（如图11（a）所示）。值得注意的是，这些区域在N估算（如图10所示）中也表现出高STD值，这表明这些像素可能代表了表现不佳的区域，即使在生长季节末期，它们的NDVI值也低于整个田块的平均水平。由于这些像素未能被准确预测fAPAR值，导致了产量的低估和实际与预测值之间的最大差距。

在氮素稀释曲线的分析中，我们发现随着生物量的增加，稀释曲线的斜率降低，而临界氮素百分比则相应提高。根据Yao等人（2021）报告的曲线，当生物量约为2吨/公顷时，临界氮素百分比可以近似为3%。在考虑这一假设的情况下，我们计算了所有生物量像素的瞬时临界氮素含量（以千克/公顷为单位），并与预测的冠层氮素含量进行了比较，这些预测值是在2021年12月28日估算的。这一比较提供了关于最终尿素施用（发生在2022年1月7日）影响的宝贵见解。

值得注意的是，尽管三个农场的面积相对较小，但利用S2高分辨率数据，本研究成功检测到了这些有限区域内的空间变化。这些建议性见解，即使不具有极高的精确度，也可以为实现更高效的资源利用和经济收益提供支持。

本研究的结论表明，将遥感数据与机器学习相结合，利用基于冠层辐射传输模型（PROSAIL-PRO）进行作物氮素状况和生物量的估算，具有显著价值。通过模拟反射率数据和Sentinel-2影像，我们训练了高斯过程回归（GPR）模型，以估算氮素浓度和生物量，其RMSE值分别为7.9千克/公顷和0.54吨/公顷。

我们的分析强调了多个因素对估算准确度的重要性，包括模拟输入数据的质量和分布、GPR中的核函数和均值函数选择，以及外部条件如作物生长阶段和数据获取时间。这些模型能够有效区分耕地和非耕地区域，但准确的作物类型识别（例如小麦）对于农业区域内的生物物理估算至关重要。

为了应对GPR模型在实际应用中的可扩展性挑战，我们建议探索更计算高效的替代方案，如线性或多项式回归，特别是在有足够的训练数据的情况下。

相关性分析揭示了氮素动态与最终产量之间的强联系，这强调了将氮素作为预测变量纳入作物生产模型以提高产量估算价值的重要性。此外，我们的工作流程在盖兹拉灌溉计划中的实际应用展示了其在使用PROSAIL驱动的GPR模型和Sentinel-2影像识别氮素应激区域的潜力。通过氮素稀释曲线映射营养缺乏区域，为更高效的氮素资源管理和改进精准农业实践提供了有前景的路径。

在本研究中，我们还开发了一个名为“nret”的Python包，该包整合了从输入数据点生成到将其分布与叶面积指数（LAI）联系起来的所有功能和脚本，并通过PROSAIL-PRO模型运行。所有用于训练GPR模型和对Sentinel-2波段进行预测的脚本和函数已公开可用。该存储库简化了整个氮素检索过程，为研究人员和农业实践者提供了一个用户友好且高效的工具。此外，它也为那些希望复制或扩展我们的研究，促进科学界透明度和协作的人们提供了宝贵资源。