基于NDVI的2017年阿根廷大豆分布及长势数据分析

孙 伟 张 峭 齐亚楠 窦红涛

（1. 中国农业科学院农业信息研究所农业部农业大数据重点实验室，北京 100081；2. 中国人民大学环境学院，北京 100872）

基于NDVI的2017年阿根廷大豆分布及长势数据分析

孙 伟1张 峭1齐亚楠2窦红涛2

（1. 中国农业科学院农业信息研究所农业部农业大数据重点实验室，北京 100081；2. 中国人民大学环境学院，北京 100872）

采用遥感手段提取作物分布并监测其长势对动态感知粮食安全具有重要意义。基于MODIS-NDVI产品，结合阿根廷大豆物候特征，快速提取阿根廷大豆分布及长势情况，并使用农用地连片系数FCI校正由尺度效应造成的面积失真。在此基础上，使用距平分析监测大豆长势，计算DSI指数监测干旱过程。结果表明：提取的近3年阿根廷大豆种植面积与美国农业部公布数据变化趋势一致，主产区种植比例高度吻合。2017年阿根廷大豆整体长势较差，副产区长势优于主产区，各主产省种植空间分布呈现离散或聚集变化。大豆种植面积及分布、长势的变化与干旱过程在时间和空间上具有高度耦合关系。提出的“物候匹配-LUCC降尺度作物提取法”具有准确快速的特点，通过改良和修正可用于其他作物的有效提取与监测。

植被指数；大豆种植面积；长势监测；阿根廷；距平分析；干旱严重指数

1 引言

遥感技术的快速发展为监测全球大尺度农作物种植及生长状况提供了可能，是动态感知粮食安全的主要手段。近年来，Landsat、SPOT、MODIS以及CBERS等卫星的影像数据已经被成功用于宏观尺度上小麦、玉米、水稻等作物的种植区提取和长势监测且得到了满意的精度[1]，但很少用于对大豆种植的分布提取和长势监测，而面对国内大豆及其产品供应偏紧的市场格局，国际大豆种植情况的变化将直接影响我国大豆进口价格。因此，研究开发用于大豆种植遥感识别和长势监测的方法与技术流程对于分析全球大豆市场行情具有重要意义。

目前，常见作物的遥感识别的方法主要有物候识别、光谱识别和物候-光谱识别[2]。物候识别容易受到生长期相同的其他作物或植被的干扰[3]，而光谱识别则存在同物异谱和异物同谱的问题[4]。因此，作物识别精度受地区农业种植结构复杂程度的影响较大。阿根廷是全球大豆主产区之一，且其作物种植结构比较单一，大豆占主导地位，采用中低分辨率的遥感影像提取种植分布具备较大的可行性。此外，作物识别精度取决于选取的遥感信息源，在国家尺度上的大豆种植监测中，对于时空分辨率的选择应更侧重时间分辨率，使数据能够覆盖作物生长期的各个时段，避免过高光谱分辨率带来的时相不佳、数据不全等系列问题[5-8]，由于降低了对空间分辨率的要求则导致监测结果因尺度效应与区域统计数据产生一定偏离，因此需要进行后期修正。本文利用MODIS植被指数产品（NDVI）时相丰富的特点，采用物候法识别大豆，并采用具有较高分辨率的农用地分布数据进行降尺度处理，从而得到精度较高的阿根廷大豆空间分布数据。基于大豆分布数据，使用距平分析法对大豆长势逐月监测，同时利用干旱指数DSI监测生长期内干旱发生情况，全面解析2017年阿根廷大豆种植及长势情况。

2 研究区与数据源

阿根廷位于南美洲南部，地理范围为：53°35′W-73°37′W，21°46′S-55°07′S。东临南大西洋，西接智利，南与南极洲隔海相望，北接巴拉圭、玻利维亚。国土面积278.1万km2，是拉丁美洲第二大国，仅次于巴西。阿根廷地形西高东低，西部是以安第斯山脉为主体的山地；东部和中部的潘帕斯草原是著名的农牧区；北部主要是格二查科平原，沼泽、森林较多；南部是巴塔哥尼亚高原。阿根廷农业用地主要分布于中东部地区，北部也有少量分布，总面积达35.32万km2。其农牧区大都位于温带和亚热带地区，天气温暖，降水量较大，湿度较高，是世界最适宜农业生产的区域之一。阿根廷农作物种植结构较为单一，大豆、玉米和小麦是阿根廷的三大主产作物，并且大豆种植面积占据绝对优势。阿根廷大豆主产区包括布宜诺斯艾利斯省、科尔多瓦省、圣菲省，副产区包括恩特雷里奥斯省、圣地亚哥-德尔埃斯特罗省、拉潘帕省等。图1为研究区位置示意图。

研究使用的数据包括联合国粮农组织各国粮食简报中提供的阿根廷作物生长期数据（http：//www.fao.org/giews/country-analysis/countrybriefs）；美国国家航空航天局地球科学数据系统项目提供的MODIS Terra星16天合成NDVI产品、8天合成净蒸散（ET）产品，空间分辨率均为500m（https：//search.earthdata.nasa.gov/search）；中国科学院地理科学与资源研究所资源环境数据中心的全球1：100万国家边界和省级行政边界数 据（http：//www.resdc.cn/Default.aspx）； 阿 根廷2010年30米土地覆盖数据，来源于全球30米土地覆盖数据集（http：//www.globeland30.org/GLC30Download/index.aspx）。在时间范围上包括2011—2017年11月至次年4月NDVI产品，2014—2017年11月至次年3月ET产品，数据较为完整地覆盖了阿根廷大豆生长周期。

图1 研究区位置示意图

3 研究方法

3.1 大豆物候特征

不同作物及同一作物在不同生长发育时期，具有不同的生理特征，这些特征在一定程度上可以通过其光谱反射和吸收特性表现出来，因此，光谱反射率的差异及其组合可以作为作物类型识别的重要依据[9]。植被指数利用了绿色植物叶片在不同光谱段的反差组合来定量描述植被特征，是到目前为止用于区分不同地物类型或植被类型以及监测植被生长状态的最有效方法[10]。而在众多的植被指数中，NDVI在农业遥感中的应用最为普遍，被认为是作物生长状态的最佳指示因子[11]。NDVI的计算方法如式（1）所示：

在大豆出苗、开花、结荚、鼓粒和成熟的整个发育过程中，NDVI会表现出相应的增减变化。当作物的叶片面积增大时，NDVI随之增大，叶面积最大时植被指数也达到最大值，而当作物进入成熟期，光合作用下降，叶绿素降低，NDVI值也相应降低，作物成熟后随着收割活动的开展会造成植被覆盖的急剧下降。阿根廷大豆播种从11月上旬开始，一直持续到12月下旬，其中11月是最重要的种植时间；1月到3月上旬为大豆的关键生长期，大部分种植区的大豆在2月下旬到3月上旬是结荚灌浆期；从3月下旬开始进入收获期，3月下旬到5月上旬大豆陆续成熟收割。一般而言，在结荚灌浆期大豆叶面积达到最大，此时NDVI也达到了峰值[12]，而在4月因作物成熟和农户的收获活动，NDVI会急剧下降。根据这一物候特征，可以利用时间序列NDVI产品提取大豆种植区。图2为大豆多年种植区11月上旬到来年4月下旬平均NDVI值时间动态变化曲线。

图2 大豆生长期平均NDVI变化曲线

3.2 大豆种植区提取流程

大豆种植区的提取包括数据预处理和物候特征匹配，其中数据预处理包括对NDVI产品进行图幅拼接、投影转换和影像裁剪以及去负值等预处理工作，坐标系统由原始的Sinusoidal转换为在南美地区变形较小的WGS_1984_UTM_19S投影。物候特征匹配主要通过最大合成NDVI与大豆结荚灌浆期时间的匹配和4月份NDVI明显下降两个特征匹配。式（2）为匹配技术原理。

在式（2）中，Rpotential表示潜在大豆种植区，当满足match中的条件时，值取1，表示是大豆潜在种植区，否则值为0，表示非大豆种植区。Max（NDVIi）为大豆生长期12个旬里NDVI最大合成值，NDVI为第i日的NDVI值，取值分别为某年第305天至次年第113天以16天为间隔的12个旬的NDVI，如NDVI065表示次年第65天的NDVI，-0.05表示后一期NDVI比前一期下降至少5%。

由于MODIS空间分辨率较低，无法区分细碎地物，使得在潜在大豆种植区中包含了其他地物。一个较为快捷和准确的方法是使用农用地连片系数FCI （Farmland Connectivity Index）进行校正，从而得到能够表示大豆种植率的空间网格Rsoybean。FCI借助具有较高分辨率的30m土地利用，可以通过统计栅格单元中农用地的面积比例得到，优点是通过高精度的土地利用/覆盖产品实现了降尺度的效果，过滤了存在于非农用地中的提取误差。而Rsoybean的优点是直观表达大豆潜在种植区混合像元的异质程度，且便于统计有效种植面积。计算方法如式（3）和式（4）所示。

在式（3）和式（4）中，Aa为某栅格单元中农用地的面积，Ag为栅格单元面积，含义同式（2），Rsoybean是FCI在潜在大豆种植区的表达，综合了大豆的物候特征和农用地分布情况，直观反映大豆空间分布及其栅格单元尺度上的种植比例，取值范围与FCI一致，介于0和1之间，而FCI则表示一个栅格单元中可用于大豆种植的概率。图3是FCI与Rsoybean的对比图。

图3 FCI与Rsoybean对比图

3.3 长势与干旱指数

本文对于大豆长势的判断是通过NDVI距平分析得到。距平分析是指通过将某时段的指标值与多年同期平均水平进行差值分析判断该时段的指标状况。NDVI的距平计算方法如式（5）。

其中，Andvi为NDVI距平值，NDVIi含义同式（2），为NDVI多年同期平均值。

干旱严重指数（DSI）是通过卫星遥感数据监测地球表面的干旱程度，它通过综合MODIS卫星数据产品NDVI和ET/PET（蒸散/潜在蒸散）来表达作物的光合作用状态和水分胁迫情况。其中，ET是植被蒸腾和土壤蒸发的总和值，PET是指在一系列最佳条件下土壤水分供应充足时的蒸散量。ET基于Penman-Monteith公式[13]计算所得，PET基于Penman公式计算所得。ET与PET的比值通常用作陆地水分可利用度和相关湿度或干旱程度的指标。DSI提供了近乎实时的干旱监测能力，不仅能够用来监测气象干旱，而且能够有效监测农业干旱[14]，从而帮助决策者进行区域干旱评估和减灾工作。相关的计算如下：

如今，当我怀揣着文学梦，亦步亦趋来到北大之后，每天我只有珍惜，只有感恩，有时，我会担心第二天一觉醒来，发现自己不在北大。我现在能做的，就是对身边的每一个人好，用我的笔勤奋地耕耘下去，纵然艰难，但身处这样一个和谐的环境内，不对立就是进步，就是人生最大的成功。

式（6）—式（10）中，NDVI、ET/PET、DSI分别为某年某旬的归一化植被指数、蒸散与潜在蒸散比以及干旱严重指数；ZNDVI、ZET/PET、ZDSI分别为标准差标准化后某年某旬的NDVI、ET/PET以及DSI；分别为多年同期平均NDVI、多年同期平均ET/PET以及多年同期平均Z值；σNDVI、σET/PET、σZ分 别 为NDVI、ET/PET以及DSI的标准差；σ为变量X的标准差，为变量X的平均值，n为X的个数。本文中n为数据的旬数。关于ET/PET的计算方法见http：//www.unc.edu/courses/2008spring/geog/577/001/www/Mu07-RSE.pdf[15]。

4 结果与分析

（1）大豆种植区空间分布及面积统计

根据大豆种植率空间网格计算得到阿根廷近3年大豆种植面积。结果表明，2015-2017年阿根廷大豆种植面积呈小幅下降趋势，2017年阿根廷大豆总的种植面积为1695.77万公顷，相比2016年和2015年分别减少8.79万公顷、43.79万公顷。图4为本研究提取面积与美国农业部USDA公布的阿根廷大豆种植面积对比结果，变化趋势一致。

图4 USDA-PSD大豆种植面积与本研究提取面积对比

2017年阿根廷大豆种植仍然以布宜诺斯艾利斯省、科尔多瓦省和圣菲省最多。这3个省的种植面积占总种植面积的70.51%。其次是恩特雷里奥斯省、圣地亚哥-德尔埃斯特罗省、拉潘帕省，种植面积占总种植面积的18.27%。因此，全国近90%的大豆都种植在上述6个省份里，是阿根廷大豆的主副产区。其他省份仅有少量种植。

除了大豆种植总面积发生小幅波动外，各省的种植面积年度变化特征也各异。2015—2017年，布宜诺斯艾利斯省大豆种植面积在数量上呈现小年→大年→小年，在分布上呈现由离散→东部、南部局部聚集→西部聚集的动态特征；科尔多瓦省大豆种植面积在数量上呈现大年→小年→大年，在分布上表现为由中部带状聚集转为中东部、南部片状离散；圣菲省大豆种植面积逐年增加，在分布上表现为整体离散，中部较集中的特征。而在副产区，种植面积主要表现为北部的圣地亚哥-德尔埃斯特罗省逐年增加、中东部的恩特雷里奥斯省先小幅增加后大幅减少以及中南部的拉潘帕省先减后增等特征，在分布上这3个省分别表现出东部集中→东部离散、中部集中→中北部离散以及东部离散→东部相对集中的分布态势。查科省、圣路易斯省、图库曼省种植面积逐年下降，分布上存在离散化趋势；萨尔塔省2017年种植面积与2015年和2016年相比有较大幅度下降，分布明显离散化。其他省份种植面积很少。种植面积与分布的变化如图5所示。

（2）阿根廷大豆长势情况距平分析

长势监测通过对当年作物生长各阶段长势与多年平均水平的比较，得出作物长势好坏的相对信息，为作物估产提供基础信息，并有助于间接识别旱涝灾害和作物病虫害等气象、农业灾害。

距平分析表明，2016—2017年度阿根廷大豆长势前期较差，后期较好，整体较差，尤其是主产区布宜诺斯艾利斯省、科尔多瓦省和圣菲省。从空间分布来看，2016年11月长势较差的地区分布在布宜诺斯艾利斯省中部和南部，至12月上旬扩展至科尔多瓦省、圣菲省中南部。12月下旬至2017年1月中旬布宜诺斯艾利斯省及其同纬度的拉潘帕省长势渐好，但科尔多瓦省长势持续较差。2017年2月长势较差的区域主要集中在圣菲省。3月上旬除科尔多瓦省和圣菲省少数区域长势较差外，其余地区长势较好。图6展现了阿根廷主副产区大豆生长期旬尺度的长势情况。

图5 阿根廷大豆种植分布及省级种植面积统计

图6 阿根廷大豆生长期旬尺度的长势状况示意图

统计结果表明，2017年阿根廷大豆长势整体较差，尤其表现为大豆种植期和出苗期主产区长势较差，主产区长势差于多年平均水平的比例，均在10%左右，而副产区基本低于10%并且多保持在5%及以下；长势差于或者略差于多年平均水平的比例主产区大都保持在35%以下而副产区大部分在30%以下；长势好于或者略好于多年平均水平的比例主产区保持在35%左右，而副产区保持在45%～65%，尤其是圣地亚哥-德尔埃斯特罗省和恩特雷里奥斯省长势很好；主产区长势好于多年平均水平的多低于10%（图7）。

（3）干旱过程监测及对种植面积和长势影响

利用MODIS的NDVI产品和ET产品，计算干旱严重指数DSI，监测阿根廷大豆生长期干旱过程。受NDVI和ET数据可得性的影响，计算得到的DSI指数没有覆盖整个大豆生长期（缺少2017年1月上旬—2月上旬和2017年3月下旬—4月下旬，其中，后一个时段为收获期，干旱对作物产量影响很小，可以不考虑）。因此，本文对2016—2017年大豆干旱情况分析仅限于2016年11月上旬—12月中旬以及2017年2月下旬—2017年3月上旬。

干旱监测结果表明，在大豆主副产区，2016年11月下旬干旱严重，并呈现逐渐加重趋势，直至2016年12月中旬；2017年2月下旬至2017年3月上旬主产区干旱缓解，干旱转移至西部非大豆种植区。在大豆种植期和出苗期的严重干旱主要发生在布宜诺斯艾利斯省，而科尔多瓦省、圣菲省以及恩特雷里奥斯省也较为干旱，其他个别省份也有部分地区发生干旱，但大豆种植面积不多。大豆结荚灌浆期主副产区长势良好，较为干旱的地区分布在布宜诺斯艾利斯省南缘以及其北部与科尔多瓦省、圣菲省交界处且面积较小，程度较轻（图8）。

综合大豆生长期干旱过程及种植面积变化、长势状况可以得出：大豆种植期和出苗期主产区的干旱过程造成了布宜诺斯艾利斯省大豆种植面积的减少，且12月上旬以前主产区长势较差，并呈现出种植区由东部、南部聚集转为西部聚集的变化。在大豆结荚灌浆期，布宜诺斯艾利斯省北部和科尔多瓦省、圣菲省交界处的轻度干旱也造成了该区域大豆长势不佳。由干旱指数分布图与年际间大豆面积变化可以看出，阿根廷大豆种植面积的减少和空间分布的变化与干旱区域分布存在高度耦合关系。

5 结论与讨论

本文提出了一种适于阿根廷大豆种植分布识别的“物候匹配-LUCC（土地利用）降尺度”提取法，利用MODIS的NDVI、ET等产品及土地利用数据，结合阿根廷大豆物候特征，对2017年阿根廷大豆种植的面积提取、长势监测和干旱监测进行了研究。结果表明，本文使用的提取流程和技术方法能够快速、准确地识别阿根廷大豆种植区，利用30m分辨率土地利用构建的FCI可有效减少尺度效应的影响从而提高种植范围的识别精度。采用美国农业部预测的PSD（产量、供给和分布）数据进行验证，本文提取出的近3年大豆种植面积均略低于该预测数据，约为其预测值的90%，具有一定的一致性与可信度。分析误差原因在于，虽然本文在一定程度上考虑了阿根廷因南北跨度而存在的物候期差异，但由于缺乏实地考察求证，在大豆物候期判定上难免存在误差，进而导致大豆面积提取存在少许偏离。在省级区域的空间分布上，本文区分出的大豆主副产区及其种植比（88.78%）与阿根廷农业部USDA公布的大豆主副产区及种植比（90%）高度一致。在一定程度上，本文可佐证美国农业部每年发布的全球主要国家农业监测数据，通过对比也证实了该方法的合理性和可靠性。应用该方法持续监测大豆主产区种植分布的空间变化、长势的时序变化、干旱发生等情况，可为全球粮食安全动态感知与预判提供基础信息。

图7 阿根廷大豆长势状况统计

充分利用低分辨率数据源的多时相信息以及较高空间分辨率遥感数据产品的地类信息，具有宏观、高效、准确的特点，非常适宜国家尺度的大片种植作物的快速提取，具有广阔的应用前景。但应注意的是，该方法并不能直接用于像我国种植结构复杂、种植面积不集中、物候差异显著的国家进行作物识别和长势监测。通过添加和修改物候匹配的约束条件可在一定程度上提高该方法的适用性。未来的研究可从物候信息校准、NDVI产品分辨率提高、土地利用时效性及精度提升、作物异种程度区分等方面进一步深入探讨。

[1]吴炳方, 蒙继华, 李强子. 国外农情遥感监测系统现状与启示[J]. 地球科学进展, 2010, 25(10)： 1003-1012.

[2]张喜旺, 刘剑锋, 秦奋, 等. 作物类型遥感识别研究进展[J]. 中国农学通报, 2014, 30(33)： 278-285.

[3]胡宗辰. 基于MODIS的中国主要粮食作物种植时空分布信息提取方法研究[D].成都： 电子科技大学,2013.

[4]朱钟正, 苏伟. 基于局部空间统计分析的SPOT 5影像分类[J]. 遥感学报, 2011, 15(5)： 957-972.

[5]YANG Wei, CHEN Jin, MATESUSHITA Bundei, et al.基于相关系数匹配的混合像元分解算法[J].遥感学报 , 2008, 12(3)： 454-461.

[6]潘耀忠, 李乐, 张锦水, 等.基于典型物候特征的MODIS-EVI时间序列数据农作物种植面积提取方法： 小区域冬小麦实验研究[J]. 遥感学报, 2011,15(3)： 578-594.

[7]陆俊, 黄进良, 王立辉, 等. 基于时空数据融合的江汉平原水稻种植信息提取[J].长江流域资源与环境,2017, 26(6)： 874-881.

[8]王玉. 基于时序光谱库的棉花种植面积信息提取研究[D].北京： 中国地质大学(北京), 2013.

[9]盛永伟, 陈维英, 肖乾广, 等.利用气象卫星植被指数进行我国植被的宏观分类[J].科学通报, 1995, 40(1)：68-71.

[10]陈晓苗. 基于MODIS-NDVI的河北省主要农作物空间分布研究[D]. 石家庄： 河北师范大学, 2010.

[11]刘真真. 基于改进CASA模型的小区域冬小麦遥感估产研究[D]. 开封： 河南大学, 2016.

[12]ZHAO H, YANG Z, DI L, et al. Crop phonology date estimation based on NDVI derived from the reconstructed MODIS daily surface re fl ectance data[C]//International Conference on Geoinformatics. IEEE, 2009：1-6.

[13]黄健熙, 张洁, 刘峻明, 等. 基于遥感DSI指数的干旱与冬小麦产量相关性分析[J]. 农业机械学报, 2015,46(3)： 166-173.

[14]MONTEITH J L. Evaporation and environment[M]//FOGG B D. The state and movement of water in living organism, symposium of the society of experimental biology. Cambridge： Cambridge University Press, 1965：205-234.

[15]MU Q, HEINSCH F A, ZHAO M, et al. Development of a global evapotranspiration algorithm based on MODIS and global meteorology data[J]. Remote Sensing of Environment, 2007, 111(4)： 519-536.

Condition Analysis on NDVI-Based Soybean Spatial Distribution and Growth in Argentina 2017

SUN Wei1, ZHANG Qiao1, QI Yanan2, DOU Hongtao2
(1.Agricultural Information Institute, Chinese Academy of Agricultural ScienceKey Laboratory of Agricultural Big Data, Ministry of Agriculture, Beijing 100081; 2.School of Environment and Natural Resources, Renmin University of China, Beijing 100872)

Timely grasp crop distribution and monitor its growth condition is of great signi fi cance to dynamic perception of food security. Based on the MODIS-NDVI product and the phonology characteristics of soybean in Argentina, we extract the distribution of soybean and monitor its growth condition, then use farmland contiguous index(FCI) to rectify the area distortion causing by scale effect. The result shows that the change trend of soybean planting area in Argentina for past three years is in line with the trend of data released by the USDA. The planting ratio in the main producing areas is also highly consistent with of ficial data. The overall growth of soybean is poor than normal in 2017, and the soybean harvest area decrease than the past two years.The sub-region is better than the main producing area, and the spatial distribution of the main producing provinces presents a discrete or aggregated change. The soybean distribution and growth condition are highly coupled with the drought degree in time and space. In this study, the method named “phonological matching-LUCC downscaling crop extraction” is comparatively accurate and rapid, it is could be used for effective extraction of soybeans in other areas or other staple crop crops.

NDVI, crop growth monitoring, growth monitoring, Argentina, distance average analysis, DSI

P237

A

10.3772/j.issn.1674-1544.2017.06.010

孙伟（1982—），女，中国农业科学院农业信息研究所副研究员，研究方向：农业遥感；张峭（1962—），男，中国农业科学院农业信息研究所研究员，研究方向：农业风险管理（通讯作者）；齐亚楠（1995—），女，中国人民大学环境学院硕士研究生，研究方向：环境规划和3S应用；窦红涛（1992—），男，中国人民大学环境学院硕士研究生，研究方向：环境规划和3S应用。

项目名称：中国农业科学院科技创新工程“全球典型农作物生产遥感监测与风险评估”（CAAS-ASTIP-2016-AII）；科技基础性工作专项重点项目“科技基础性工作数据资料集成与规范化整编”（2013FY110900）。

2017年8月23日。