植被指数与美豆单产的相关性探讨

北美大豆的收获正在进行中，虽然定产预估需等待明年1月USDA报告公布，但是新作产量的波动对市场的影响已在盘面得到反映。2022年北美干旱问题贯穿始终，8月下旬Pro Farmer的田野调查显示多数州的大豆、玉米单产值都显著低于3年平均水平，数据公布当周CBOT大豆、玉米、小麦合约均表现出较大程度的涨幅。单产的变化对市场的影响无疑，因此利用遥感数据来判断单产预估将是本文探讨的重点。 在上一期的遥感专题《另眼看风景——以植被指数为例初探遥感数据在农产品研究中的应用》中，我们谈到了关于NDVI、EVI的基本概念、季节性特征，以及利用植被指数进行土地覆盖类型分类的问题。在此基础之上，本文将研究聚焦到具体的品种（美国大豆），选取了9个主产州（Iowa, Illinois, Indiana, Ohio, Missouri, Minnesota, Nebraska, South Dakota, North Dakota）作为对象，进一步探讨植被指数和单产的相关性。 首先介绍文中使用的数据来源。第一，美国大豆种植区域的数据来自USDA NASS CDL（Cropland Data Layer），该数据覆盖时间自1997年至2021年，结合卫星数据和实地数据，利用监督学习的方式给出某自然年度内不同农作物的种植区域（相关介绍参见https://www.nass.usda.gov/Research_and_Science/Cropland/sarsfaqs2.php）。因为2022年数据尚未提供，故分析时以2021年数据替代。第二，NDVI、EVI数据来自LP DAAC MOD13Q1（6.1版本）250米16天 合 成 数 据 （ 相 关 介 绍 参 见https://lpdaac.usgs.gov/products/mod13q1v061/）。第三，美国各州单产数据来自USDA NASS数据库。 图 1 2021年美国大豆种植区域（黑色部分） 数据来源：NASS CDL，国投安信期货 大豆的生长期NDVI曲线特征 以2021年美国大豆主产州的NDVI时间序列为例，通过NDVI的变化来反映大豆生长期的过程。从图一中能总结出两个较明显的特点：一、从时间上看，自5月开始NDVI的增长速度变快，其最大值一般发生在8、9月份，这代表了大豆从开始种植到结荚鼓粒成熟的时期。从9月之后，NDVI开始下降，且9-10月的下降速度较快，反映大豆已进入收获期。到11月收获近结束时，NDVI下降速率已然趋缓。二、从空间上看，在2021年，位置偏北的各州，例如North Dakota，South Dakota，Minnesota等的NDVI整体上低于位置偏南的各州，例如Indiana, Ohio, Missouri等，且其到达峰值、开始下降的时间点更早一些。 从图3各州的NDVI季节性图上则能反映出更多细节问题。 首先，NDVI值是否对单产的预测具备意义。诸多文献资料已经针对NDVI在预测单产上的作用方面有过研究并得出积极的结论。这里我们定性地观察一下，从North Dakota的NDVI季节性图上找到一条明显低于其他年份的曲线，即2021年的NDVI。而2021年North Dakota的单产为25.5蒲式耳/英亩，是自2005年以来的最低值，这一点符合NDVI的表现。但这样较极端的情况比较少见，且有可能是异常点，因此我们将9个主产州2000-2021年的单产和NDVI数据进行了相关性分析，具体讨论见下一节。 其次，生长期内不同阶段NDVI的重要性。图3中多个主产州在2019年的NDVI曲线都明显向右偏移，说明该年度大豆种植的时间处于过去6年的偏晚位置。结合2019年大豆生长期的crop progress数据（图4下），可以看到，2019年从种植到收获的进度都是明显低于2018年和过去5年平均值，这一点与NDVI的表现非常吻合。回顾2019年美国大豆的种植情况，降雨在其中的影响非常显著。从图5可以看出，从5月中下旬开始，2019年的降雨量明显处于近5年高位，在种植窗口期多个地区的降雨量几乎翻倍。降雨对作物单产的影响可以从优良率反映出来（图4上），2019年整个生长期的优良率都位于5 年低位。再结合NDVI的曲线，可以发现，尽管曲线反映出生长期向后推迟的特征，但NDVI最大值的数值变化却不太明显，所以只使用生长期内最大NVDI和单产做相关性分析还不足以解释单产的变化。因此我们将生长期内各个月份的NDVI分别和单产做相关性和一元线性回归，试图找出关键时期的NDVI，具体讨论见下一节。 图 2 2021年美国大豆产区NDVI时间序列 数据来源：LP DAAC MOD13Q1，国投安信期货 图 3 美国大豆主产州NDVI季节性 数据来源：LP DAAC MOD13Q1，国投安信期货 图 4 2019年美国大豆生长优良率（上图）与进展（下图） 数据来源：USDA NASS 图 5 大豆主产区降雨情况 数据来源：路透，国投安信期货 植被指数与大豆单产的相关性 接下来讨论植被指数与大豆单产的关系。数据的时间为2000年至2021年，范围包括9个大豆主产州，选取大豆生长期内（4月-10月）的植被指数和当年单产为研究对象。因为数据量较为有限，此处将不同州、不同年份的数据混置在一起研究，并未考虑不同地区或者不同年份的影响。 从NDVI与单产的散点图（图6）中，可以发现7月27日至8月12日、8月13日至8月28日这两段时间的NDVI与单产具有较为明显的正向相关关系，分别为0.606、0.648。若以EVI作为观察对象，也得到类似的结果，只是相关系数较高的时间段多了7月12日-27日。具体的统计结果见表1。 表格 1 不同时间段植被指数与单产相关系数 时间段 相关系数（NDVI） 相关系数（EVI） 4/7 – 4/22 -0.205 -0.114 4/22 – 5/7 -0.260 -0.102 5/9 – 5/24 -0.374 -0.183 5/24-6/9 -0.274 -0.088 6/10-6/25 0.015 0.153 6/25-7/10 0.358 0.438 7/12-7/27 0.548 0.602 7/27-8/12 0.606 0.627 8/13-8/28 0.648 0.668 8/28-9/13 0.358 0.422 9/14-9/29 0.021 0.100 9/29-10/14 -0.153 -0.155 数据来源：国投安信期货 图 6 不同时段大豆NDVI与单产散点图 数据来源：国投安信期货 在和单产的相关性比较中，NDVI和EVI的表现一致，同时EVI与单产的相关性略强于NDVI。接下来，我们利用2000年至2021年7月27日-8月12日以及8月13日-28日的NDVI、EVI分别与当年的单产建立一元回归模型，并将2022年的数据作为测试集。从结果来看，以8/13-8/28时间段内的植被指数为解释变量的训练集拟合效果最佳，但从测试集结果来看，7/27-8/12时间段内的效果略好，且以NDVI作为解释变量的效果好于EVI。从图7测试集的真实值和拟合值来看，有一个点的拟合效果不佳（由于数据集 有限，无法判断该点是否是异常值）。该点代表的植被指数显著低于其他数据点，但是单产却没有明显的降低，由此仅依靠植被指数来拟合单产无法很好地解释这类数据。从测试集的效果来看，2022年单产的拟合值有一定的解释度，但明显比真实值偏低。因为2022年的单产非定产数据，且CDL沿用2021年的数据，因此和训练集的处理方式有一定差异，这可能为结果带来了偏差。另外，回归的效果和训练集和测试集的分组也有关系，此处将2022年的数据作为测试组，2000-2021年的数据作为训练组，得到的可能并不是最好的回归效果，可以考虑使用k-fold 交叉验证进一步优化。 表格 2 回归评价 解释变量 时间段 相关系数 训练集R2 测试集R2 NDVI 7/27-8/12 0.606 0.367 0.273 8/13-8/28 0.648 0.420 0.269 EVI 7/27-8/12 0.627 0.393 0.165 8/13-8/28 0.668 0.447 0.144 数据来源：国投安信期货 图 7 以7/27-8/12 NDVI为解释变量的拟合效果（左：训练集，右：测试集） 数据来源：国投安信期货 结论 7月末-8月末的NDVI、EVI均与美豆单产有较强的正相关关系，其中EVI的相关系数高于NDVI。在以NDVI、EVI分别和美豆单产建立的一元回归模型中，7/27-8/12 时间段内的NDVI表现最优，且对单产变化具有一定的解释力度。因此在进行单产预测时，可以考虑将7月末-8月末这一时间段内的NDVI或EVI作为因子之一。而简单的一元线性回归拟合度仍有提升空间，且出现预测值普遍偏低的问题，因此在下一次的期市月月谈中，我们将拓展单产预测模型，并进行更多的讨论。