一种基于Sentinel-2遥感影像的大蒜种植信息自动化提取方法与流程

本发明涉及农业遥感技术领域，特别涉及一种基于sentinel-2遥感影像的大蒜种植信息自动化提取方法，利用大蒜生长过程中独特的物候特征，在不依赖于地面数据的情况下实现大蒜种植信息的自动化提取。

背景技术：

农作物种植信息的实时监测和精细制图是农业遥感的基础，详尽的作物类型分布信息对粮食安全、政策制定以及社会经济等有着重要的意义。世界粮食安全首脑会议宣称，预计到2050年，世界总人口将增长到近100亿，在经济温和增长的情况下，农业需求将比2013年增长约50％。粮食产量的增加必然伴随着对农田的可持续管理，使得能够停止或至少减缓对水土资源质量和数量、土地退化、生物多样性以及温室气体排放等多方面的负面影响。鉴于此，在不同的地点和环境背景下，对农作物生长状况的实时监测则显得意义非凡。

由于遥感影像具有覆盖范围广、更新速度快、长时序数据监测等优点，所以能够提供一种高效的农作物种植信息获取方法。随着遥感技术的不断发展，已经涌现出许多性能更加优越的传感器，空间分辨率、时间分辨率以及光谱分辨率均得到了很大的改进和提升。特别是最近发射的sentinel-2系列卫星，不仅拥有较高的10m空间分辨率，而且还具有幅宽大、重访周期短、数据免费、波段信息丰富以及较多的红边波段等特点，这些特点在作物识别和面积提取方面具有很大的优势。目前，在基于遥感数据的农作物种植信息监测方面，已有的研究大多都聚焦于冬小麦、水稻、玉米等大面积种植的大众作物，对大蒜种植信息的监测却鲜有研究。

大蒜是一种百合科葱属多年生草本植物，日常生活中被人们作为普遍食用的传统蔬菜和重要调味品；除此之外，多年的科学研究证明，大蒜具有显著的杀菌、抗癌等药用价值，是一种珍贵的保健食品。近年来，大蒜的价格受其产量等方面的影响，波动变化剧烈，2010年网络上更是一度有着“蒜你狠”的调侃热词，究其根本是因为产量信息的不合理性。中国是迄今为止世界上最大的大蒜生产国，产量占据全球的70％以上。因此，及时而精确的掌握我国大蒜的种植及产量等信息对于指导市场经济有着重要的意义。

目前，在农业遥感领域，对农作物种植信息的提取主要依赖于机器学习等监督分类算法，该类方法及其依赖于训练数据，而数据的获取通常会耗费大量的人力物力；所以，如何寻求其它不依赖训练数据的作物自动提取算法是本研究领域技术人员亟待解决的问题。本发明旨在提出一种自动的大蒜提取算法。需要注意的是，大蒜是一种冬季作物，而研究区中通常会存在物候特征和大蒜相近的另一种冬季作物冬小麦，因此，算法的关键在于如何避免冬小麦对大蒜提取的干扰。

技术实现要素：

本发明目的是提供一种基于sentinel-2遥感影像的大蒜种植信息自动化提取方法，sentinel-2号系列卫星具有幅宽大、重访周期短、数据免费、波段信息丰富以及较多的红边波段等特点，这些特点在作物识别和面积提取方面具有很大的优势。为了减少对人工标记样本数据的依赖，极大地解放人力物力，同时尽可能地使我们提出的方法具有良好的可迁移性，本发明引入了一种和叶面积指数(lai)变化密切相关的绿色叶绿素植被指数(greenchlorophyllvegetationindex,gcvi)，并结合经典的归一化植被指数(normalizeddifferencevegetationindex,ndvi)，充分探讨大蒜和冬小麦的叶面积指数及光谱特征之间的差异，构建相关的提取模型，实现大蒜种植信息的自动化提取。

为达到上述目的，本发明采用的技术方案是：一种基于改进gcvi指数的大蒜种植信息自动化提取模型，包括以下步骤：

步骤1，获取研究区大蒜生长周期内的所有可用sentinel-2遥感影像；

步骤2，对获得的所有sentinel-2遥感影像进行预处理操作，包括大气校正、重采样以及影像裁剪拼接等；

步骤3，利用步骤2中经过预处理后的sentinel-2数据，生成覆盖整个生长周期的ndvi时间序列和仅覆盖越冬后到成熟收获这段时期的gcvi时间序列曲线；

步骤4，利用步骤3中获得的gcvi时间序列，结合大蒜和冬小麦两种作物之间的光谱特征差异，生成改进后的gcvi时间序列；

步骤5，利用步骤3中获得的ndvi时间序列，基于大蒜和冬小麦独特的变化特征，结合ndvi最大值以及曲线二次拟合后的系数建立相应的阈值判断规则，生成冬季作物掩膜，获得过滤掉其它地物的大蒜和冬小麦分布图；

步骤6，利用步骤4中获得的改进gcvi时间序列，通过对曲线积分值的判断生成冬小麦结果，并结合步骤5中获得的大蒜和冬小麦分布图，得到最终的大蒜提取结果。

而且，步骤1中，获取的所有sentinel-2影像均分布在大蒜一个生长周期内。大蒜有六个独特的物候阶段，分别为发芽期、幼苗期、越冬期、花芽和鳞芽分化期、鳞茎膨大期以及最后的生理休眠期，每个阶段都要有对应的一景或多景sentinel-2影像；同时影像云覆盖应低于50％，对于部分云覆盖严重的影像，要事先剔除。

而且，步骤2中，预处理操作包括以下具体步骤：利用sen2cor插件对获取的所有sentinel-2原始大气顶端反射率(toa)l1c级数据进行大气校正，获得大气底层反射率(boa)l2a级数据；并对大气校正之后的数据进行重采样，使得波段分辨率统一；最后利用研究区边界矢量数据对影像进行相应的裁剪拼接操作。

而且，步骤3中，对所有经过预处理的sentinel-2数据均进行波段运算操作，分别生成ndvi数据和gcvi数据，并按时间顺序分别合成ndvi时间序列和gcvi时间序列。ndvi时间序列针对的是大蒜整个生长周期内的所有影像，其中大蒜在整个周期内的ndvi值略低于冬小麦；而gcvi时间序列仅针对越冬后到成熟收获这段时期的影像，其中该时期内大蒜的gcvi值远小于冬小麦。

ndvi的表达式如下所示：

其中，nir表示近红外波段反射率，red表示红光波段反射率。

gcvi的表达式如下所示：

其中，nir表示近红外波段反射率，green表示绿光波段反射率。

而且，步骤4中，在越冬后到成熟收获这段时期内的大蒜和冬小麦的光谱曲线中，始终存在着这样的差异，即在可见光b2(蓝光)、b3(绿光)和b4(红光)等波段范围内，大蒜的光谱反射率要高于冬小麦；而在近红外b8等波段范围内，大蒜的光谱反射率则低于冬小麦。基于这些差异，我们可以构造出这样的改进后的gcvi指数，即用原始gcvi值减去红光波段和近红外波段的比值与一个缩放系数α的乘积，这样做不仅可以进一步放大大蒜和冬小麦之间的差异，而且可以有效地避免云阴影带来的误差。因为引入了缩放系数α，所以原始gcvi表达式中的常数项可以省略掉，最终的改进gcvi指数表达式如下所示：

其中，nir、green和red分别代表近红外波段、绿光波段和红光波段的反射率，α为缩放系数。当α取值合理时，便可以有效地放大大蒜和冬小麦之间的差异，使得可以更好地对大蒜进行提取。经过多次实验探究，α的取值为15时比较合理，实际应用中可根据情况略微调整。

而且，步骤5中，利用构建好的覆盖整个生长周期的ndvi时间序列，设计相应的冬季作物掩膜，提前滤除一些非冬季作物的其它地物。在整个ndvi时间序列中，大蒜和冬小麦的最大值均可以超过0.5，而其它的一些地物在整个生长周期内的ndvi值最大都不会高于这个值，比如水体、建筑物等。因此，可基于这一特征滤除掉一些其它地物。同时，在整个ndvi时间序列中，大蒜和冬小麦会呈现独特的“双峰”特征，即在越冬前和越冬后的ndvi时间序列中分别会有一个小峰值，其中第二个峰值要大于第一个，而其它地物不会出现该特征。所以可分别对包含两个峰值的ndvi时间序列进行二次拟合，通过拟合系数的判定进一步滤除掉其它地物，最后只留下大蒜和冬小麦。具体实施步骤如下：首先，通过ndvi最大值的判断，滤除掉一些不可能为冬季作物的地物，判定公式如下所示：

其中，ndvimax表示ndvi时间序列中的最大值；ti表示sentinel-2影像获取的日期；i表示所有可用的sentinel-2影像按成像日期先后顺序排列的次序。在表达式中，我们将判定阈值设置为0.5。当ndvi最大值小于该阈值时，我们将其判定为其它地物，当ndvi最大值大于该阈值时，我们将其判定为可能冬季作物，此时，需进一步对其进行判定识别。

其次，通过对包含两个峰值的ndvi时间序列分别进行二次拟合，根据二次拟合系数的值进一步滤除掉其它所有地物，拟合公式如下所示：

y1＝a1t²+b1t+c1(5)

y2＝a2t²+b2t+c2(6)

其中，y1表示由包含第一个峰值的第一部分时间序列拟合得到的函数；y2表示由包含第二个峰值的第二部分时间序列拟合得到的函数；t表示自变量；a1和a2分别为拟合得到的两条曲线的二次系数，b1和b2分别为拟合得到的两条曲线的一次系数，c1和c2分别为常数项。当两个二次系数a1和a2同时小于0时，才可以将其判定为冬季作物，否则为其它地物。

而且，步骤6中，对缩放系数定为15的改进gcvi时间序列进行积分操作，可以对冬小麦进行识别提取。其中求积分的表达式如下所示：

其中，modified_gcviint表示改进gcvi时间序列积分值；t0和tn分别表示越冬后到成熟收获这段时期第一景和最后一景sentinel-2影像的日期；f(x)表示改进gcvi时间序列多项式函数。在求得的积分结果中，只有冬小麦的积分值大于0，而其它地物得值均小于0。因此，依据这一点可以获得冬小麦的提取结果。最后结合步骤5得到的冬季作物提取结果，即可获得最终的大蒜提取结果。

相较于现有技术，本发明有以下有益效果是：

1、针对目前农业遥感中缺乏大蒜种植信息提取相关研究的问题，本发明提出了一种大蒜种植信息提取方法，并通过相关实验证明该方法可以有效地对大蒜主产区的种植信息进行提取，在一定程度上填补了农业遥感领域中的空白。

2、本发明中提出的大蒜种植信息提取方法为自动化的提取算法，相较于传统的利用机器学习等监督分类的算法，本方法不需依赖训练样本，可以有效地节省人工标定训练样本的时间；同时，方法具有良好的可迁移性，能够应用于大尺度范围内的大蒜种植信息提取。

附图说明

图1是整个大蒜生长周期内的ndvi时间序列；

图2是越冬后到成熟收获这段时期的gcvi时间序列；

图3是大蒜和冬小麦的光谱曲线。

图4是大蒜提取总体流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图对本发明进一步说明。应当理解，此处所描述的具体实施例仅仅用于解释本发明，并不用于限定本发明。

步骤1，获取研究区大蒜生长周期内的所有可用sentinel-2遥感影像。

本发明主要以我国山东省金乡县、河南省杞县和中牟县、河北省永年县以及江苏省邳州市和丰县六个大蒜主产区作为研究区，获取2018年10月到2019年5月这一大蒜生长周期内的所有可用sentinel-2影像。大蒜有六个独特的物候阶段，分别为发芽期、幼苗期、越冬期、花芽和鳞芽分化期、鳞茎膨大期以及最后的生理休眠期，每个阶段都要有对应的一景或多景sentinel-影像；同时需要对获取的影像进行筛选，有些影像因为传感器等问题，会出现残缺，此时影像不能使用，需剔除；同时影像云覆盖应低于50％，对于部分云覆盖严重的影像，同样要剔除。

步骤2，对获得的所有sentinel-2遥感影像进行预处理操作，包括大气校正、重采样以及影像裁剪拼接等。

其具体操作为：首先进行大气校正，获取的sentinel-2影像为经过辐射校正和几何校正的标准l1c级大气顶端反射率(toa)数据，不能直接用于地物分类提取，需要结合欧空局的sentinel应用平台sentinelapplicationplatform(snap)，利用sen2cor插件，将原始的l1c级数据进行大气校正，处理为大气底层反射率(boa)l2a级数据；其次进行重采样操作，经过大气校正之后的l2a级数据由13张jpeg-2000图像组成，对应13个sentinel-2光谱波段，具有三种不同的空间分辨率，分别位10m、20m和60m，所以为了便于后续数据的统一处理，需要在snap中以b2波段为依据重采样到10m分辨率，并转换成envi格式；最后进行影像裁剪拼接操作，利用各研究区的矢量数据，对处理后的sentinel-2影像进行裁剪，得到仅覆盖研究区的影像。其中邳州市不能由一景影像完全覆盖，所以首先需要对覆盖邳州市的所有影像进行无缝拼接操作，然后再裁剪获得能完全覆盖邳州市的拼接影像。

步骤3，利用步骤2中经过预处理后的sentinel-2数据，生成覆盖整个生长周期的ndvi时间序列和仅覆盖越冬后到成熟收获这段时期的gcvi时间序列曲线。

在envi中对所有经过预处理的sentinel-2数据均进行bandmath波段运算操作，分别生成相应的ndvi数据和gcvi数据，并在envi中利用layerstacking工具，对覆盖大蒜整个生长周期的所有sentinel-2影像合成ndvi时间序列，如图1所示，可以发现整体而言大蒜的ndvi值要略低于冬小麦；同时合成越冬后到成熟收获这段时期的gcvi时间序列，如图2所示，可以发现在这段时期大蒜的值要远小于冬小麦。

ndvi的表达式如下所示：

其中，nir表示近红外波段反射率，red表示红光波段反射率。

gcvi的表达式如下所示：

其中，nir表示近红外波段反射率，green表示绿光波段反射率。

步骤4，利用步骤3中获得的gcvi时间序列，结合大蒜和冬小麦两种作物之间的光谱特征差异，生成改进后的gcvi时间序列。

如图3所示，在越冬后到成熟收获这段时期内的大蒜和冬小麦光谱曲线中，始终存在着这样的差异，即在可见光b2(蓝光)、b3(绿光)和b4(红光)等波段范围内，大蒜的光谱反射率要高于冬小麦；而在近红外b8等波段范围内，大蒜的光谱反射率则低于冬小麦。基于这些差异，我们可以构造出这样的改进后的gcvi指数，即用原始gcvi值减去红光波段和近红外波段的比值与一个缩放系数α的乘积，这样做不仅可以进一步放大大蒜和冬小麦之间的差异，而且可以有效地避免云阴影带来的误差。因为引入了缩放系数α，所以原始gcvi表达式中的常数项可以省略掉，最终的改进gcvi指数表达式如下所示：

其中，nir、green和red分别代表近红外波段、绿光波段和红光波段的反射率，α为缩放系数。当α取值合理时，便可以有效地放大大蒜和冬小麦之间的差异，使得可以更好地对大蒜进行提取。经过多次实验探究，α的取值为15时比较合理，实际应用中可根据情况略微调整。

步骤5，利用步骤3中获得的ndvi时间序列，基于大蒜和冬小麦独特的变化特征，结合ndvi最大值以及曲线二次拟合后的系数建立相应的阈值判断规则，生成冬季作物掩膜，获得过滤掉其它地物的大蒜和冬小麦分布图。

利用构建好的覆盖整个生长周期的ndvi时间序列，设计相应的冬季作物掩膜，提前滤除一些非冬季作物的其它地物。如图1所示，在整个ndvi时间序列中，大蒜和冬小麦的最大值均可以超过0.5，而其它的一些地物在整个生长周期内的ndvi值最大都不会高于这个值，比如水体、建筑物等。因此，可基于这一特征滤除掉一些其它地物。同时，在整个ndvi时间序列中，大蒜和冬小麦会呈现独特的“双峰”特征，即在越冬前和越冬后的ndvi时间序列中分别会有一个小峰值，其中第二个峰值要大于第一个，而其它地物不会出现该特征。所以可分别对包含两个峰值的ndvi时间序列进行二次拟合，通过拟合系数的判定进一步滤除掉其它地物，最后只留下大蒜和冬小麦。具体实施步骤如下：首先，通过ndvi最大值的判断，滤除掉一些不可能为冬季作物的地物，判定公式如下所示：

其中，ndvimax表示ndvi时间序列中的最大值；ti表示sentinel-2影像获取的日期；i表示所有可用的sentinel-2影像按成像日期先后顺序排列的次序。在表达式中，我们将判定阈值设置为0.5。当ndvi最大值小于该阈值时，我们将其判定为其它地物，当ndvi最大值大于该阈值时，我们将其判定为可能冬季作物，此时，需进一步对其进行判定识别。

其次，通过对包含两个峰值的ndvi时间序列分别进行二次拟合，根据二次拟合系数的值进一步滤除掉其它所有地物，拟合公式如下所示：

y1＝a1t²+b1t+c1(5)

y2＝a2t²+b2t+c2(6)

其中，y1表示由包含第一个峰值的第一部分时间序列拟合得到的函数；y2表示由包含第二个峰值的第二部分时间序列拟合得到的函数；t表示自变量；b1和b2分别为拟合得到的两条曲线的一次系数，c1和c2分别为常数项，a1和a2分别为拟合得到的两条曲线的二次系数。当这两个系数同时小于0时，才可以将其判定为冬季作物，否则为其它地物。

步骤6，利用步骤4中获得的改进gcvi时间序列，通过对曲线积分值的判断生成冬小麦结果，并结合步骤5中获得的大蒜和冬小麦分布图，得到最终的大蒜提取结果。

对缩放系数定为15的改进gcvi时间序列进行积分操作，可以对冬小麦进行识别提取。其中求积分的表达式如下所示：

其中，modified_gcviint表示改进gcvi时间序列积分值；t0和tn分别表示越冬后到成熟收获这段时期第一景和最后一景sentinel-2影像的日期；f(x)表示改进gcvi时间序列多项式函数。在求得的积分结果中，只有冬小麦的积分值大于0，而其它地物得值均小于0。因此，依据这一点可以获得冬小麦的提取结果。最后结合步骤5得到的冬季作物提取结果，即可获得最终的大蒜提取结果。整体的大蒜提取流程如图4所示。

以上所述仅为本发明的较佳实施例，但本发明不局限于所描述的实施方式。对于本领域的技术人员而言，在不脱离本发明的原理和精神的情况下对实施方式进行的变化、修改、替换和变形，皆应属本发明的保护范围内。