一种玉米收割进度高时空分辨率监测方法

技术领域

本发明涉及一种玉米高时空分辨率监测方法。

背景技术

东北地区是我国玉米的主产区，其玉米产量居主要农作物之首。玉米是我国主要的粮食作物之一，对于我国粮食安全具有重要意义。对玉米的收割进度进行及时的监测，有利于及时调整农业管理决策，保障玉米的品质和产量。

然而，现在进行农作物监测的主要方法仍然依靠人工采样调查。该种方法需要大量的时间、人力和物力，不仅不能实现大面积监测，而且监测周期长，不能满足数据的时效性要求。遥感技术作为一种“遥远的观测”技术，能够在较短周期内获取大面积地物的时空动态变化数据。因此，遥感技术能够实现大面积动态制图，是一种较好的监测玉米的收割进度的方法。

目前，利用遥感技术进行农作物收割进度监测，主要是利用光学影像对冬小麦的收割进度进行监测。主要的方法分为两类，一类是基于NDVI阈值分割的方法对单时相数据进行分类，另一类是对多时相的数据进行分类。根据多期数据构建NDVI时间序列曲线，由冬小麦前后NDVI时间序列曲线的波动寻找阈值，进行分类。但是NDVI极易受到植被含水量的影响，在玉米收割后的变化不稳定，因此其监测精度有待提高。而且，光学影像在玉米收割进度监测方面也存在一定的局限性。

光学数据受到云、雾、雪等天气状况的影响较大，极易造成影像的缺失，因此容易造成玉米收割进度的监测周期较长。雷达数据作为遥感技术中的另一类影像，不易受到天气状况的影响，能够实现玉米收割进度全天时和全天候的监测。但是，雷达影像没有丰富的光谱信息，并且单一的雷达数据其重访周期较长，无法满足玉米收割进度的监测周期。

发明内容

本发明提供了一种玉米收割进度高时空分辨率监测方法，是一种协同Sentinel-1和Sentinel-2数据且基于时序误分校正算法的玉米收割进度高时空分辨率监测方法。

本发明玉米收割进度高时空分辨率监测方法按照以下步骤进行：

一、确定待监测地区玉米收割时间，获取收割前后的Sentinel-1和Sentinel-2影像；

二、对Sentinel-1数据和Sentinel-2数据进行预处理；

三、根据玉米收割前和收割后在影像上的变化进行目视解译，提取样本点数据；

四、根据样本点数据，提取并分析Sentinel-1影像的极化特征和Sentinel-2光谱反射率曲线；

五、对Sentinel-1影像和Sentinel-2影像构建初步的特征集；

六、对Sentinel-1影像和Sentinel-2影像特征集进行选择，形成特征集；

七、利用阈值分割方法和随机森林方法基于特征集对Sentinel-1和Sentinel-2数据进行分类，并比较精度；

八、Sentinel-1影像和Sentinel-2影像对分类方法进行选择；

九、对分类后的影像进行一般的分类后处理；

十、进行时间序列误分校正的特殊分类后处理方法，对分类精度进行优化；

十一、绘制玉米收割进度的时空动态分布图，并统计累计玉米累计收割面积和收割比例；即实现了玉米收割进度高时空分辨率监测。

本发明联合了Sentinel-1雷达影像和Sentinel-2光学影像对玉米收割进度进行监测。在缩短监测周期的同时，结合了后向散射系数信息和丰富的光谱信息，大大提高玉米收割进度监测的精度。本发明以黑龙江省友谊农场为例，对友谊农场的玉米收割进度进行了动态监测，其Sentinel-1雷达影像选择阈值分割的方法，Sentinel-2光学影像选择随机森林方法，总体精度(Kappa)分别可以达到0.97(0.94)和0.99(0.99)。

本发明对于农业生产决策具有极其重要的作用。

附图说明

图1为实施例1中时间序列误分校正方法；

图2为实施例1玉米累计收割面积和收割比例。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

具体实施方式一：本实施方式玉米收割进度高时空分辨率监测方法按照以下步骤进行：

一、确定待监测地区玉米收割时间，获取收割前后的Sentinel-1和Sentinel-2影像；

二、对Sentinel-1数据和Sentinel-2数据进行预处理；

三、根据玉米收割前和收割后在影像上的变化进行目视解译，提取样本点数据；

四、根据样本点数据，提取并分析Sentinel-1影像的极化特征和Sentinel-2光谱反射率曲线；

五、对Sentinel-1影像和Sentinel-2影像构建初步的特征集；

六、对Sentinel-1影像和Sentinel-2影像特征集进行选择，形成特征集；

七、利用阈值分割方法和随机森林方法基于特征集对Sentinel-1和Sentinel-2数据进行分类，并比较精度；

八、Sentinel-1影像和Sentinel-2影像对分类方法进行选择；

九、对分类后的影像进行一般的分类后处理；

十、进行时间序列误分校正的特殊分类后处理方法，对分类精度进行优化；

十一、绘制玉米收割进度的时空动态分布图，并统计累计玉米累计收割面积和收割比例；即实现了玉米收割进度高时空分辨率监测。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤二中预处理的方法：对Sentinel-1数据采用7*7的窗口大小进行Lee滤波，抑制斑点噪声；对Sentinel-2数据，进行去云处理；然后分别对Sentinel-1雷达数据和Sentinel-2光学数据根据待监测地区农作物分布数据进行裁剪处理。其他步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一不同的是：步骤四中提取并分析Sentinel-1影像的极化特征：VV极化方式、VH极化方式、VH/VV的后向散射系数和后向散射系数的比值。其他步骤及参数与具体实施方式一相同。

具体实施方式四：本实施方式与具体实施方式一不同的是：步骤四中Sentinel-2光谱反射率曲线：Sentinel-2影像的B2(Blue),B3(Green),B4(Red),B8(NIR),B11(SWIR1)and B12(SWIR2)的光谱反射率曲线。其他步骤及参数与具体实施方式一相同。

具体实施方式五：本实施方式与具体实施方式一不同的是：步骤五中构建初步特征集：Sentinel-1光学影像上，选择B2(Blue)、B3(Green)、B4(Red)、B8(NIR)、B11(SWIR1)和B12(SWIR2)的光谱反射率和NDVI、NDTI、STI、NDI7、VTI、SVTI指数作为初步特征集；在Sentinel-2雷达影像上，选择VV、VH和RI(VH/VV)作为初步特征集。其他步骤及参数与具体实施方式一相同。

具体实施方式六：本实施方式与具体实施方式一不同的是：形成特征集的选择方法：通过对初步特征集绘制直方图和特征重要性排序，选择出能够区分玉米收割前和玉米收割后的特征集，分别为Sentinel-2光学影像的B2(Blue)、B3(Green)、B4(Red)、B8(NIR)、B11(SWIR1)和B12(SWIR2)的光谱反射率和NDVI、NDTI、STI、NDI7、VTI、SVTI指数以及Sentinel-1雷达影像的VH极化方式的后向散射系数。其他步骤及参数与具体实施方式一相同。

本实施方式形成特征集是通过绘制直方图和特征重要性排序，选择出区分玉米收割前和玉米收割后较好的特征集。

具体实施方式七：本实施方式与具体实施方式一不同的是：步骤七中分别对Sentinel-1雷达影像和Sentinel-2光学影像进行阈值分割方法进行分类和随机森林分类，比较两种分类方法对于雷达影像和光学影像的适用性。其他步骤及参数与具体实施方式一相同。

具体实施方式八：本实施方式与具体实施方式一不同的是：步骤八中对分类方法进行选择：Sentinel-1雷达影像选择阈值分割方法进行玉米收割进度监测，Sentinel-2光学影像选择随机森林分类方法进行玉米收割进度监测。其他步骤及参数与具体实施方式一相同。

具体实施方式九：本实施方式与具体实施方式一不同的是：步骤九中一般的分类后处理是进行去除小图斑(Majority and Clump)的分类后处理。其他步骤及参数与具体实施方式一相同。

实施例1 2021年在黑龙江省友谊农场，利用本发明的方法对其玉米收割进度进行高时空分辨率监测，具体步骤如下：

一、确定待监测地区玉米收割时间，获取Sentinel-1和Sentinel-2影像；

二、对Sentinel-1数据和Sentinel-2数据进行预处理；

三、根据玉米收割前和收割后在影像上的变化进行目视解译，提取样本点数据；

四、根据样本点数据，提取并分析Sentinel-1影像的极化特征和Sentinel-2光谱反射率曲线；

五、对Sentinel-1影像和Sentinel-2影像构建后向散射系数和指数等初步的特征集；

六、对Sentinel-1影像和Sentinel-2影像的特征集进行选择，形成最优的特征集；

七、利用阈值分割方法和随机森林方法基于特征集对Sentinel-1和Sentinel-2数据进行分类，并比较精度；

八、对分类方法进行选择；

九、进行一般的分类后处理；

十、进行时间序列误分校正的特殊分类后处理方法，对分类精度进行优化；

十一、绘制玉米收割进度的时空动态分布图，并统计累计玉米累计收割面积和收割比例；即实现了玉米收割进度高时空分辨率监测。

其中，步骤一中，黑龙江省友谊农场玉米的收割时间集中于9月末到11月初，从Google Earth Engine平台上获取友谊农场2021年9月末至11月初的Sentinel-1数据和Sentinel-2数据；

步骤二中，对Sentinel-1数据采用7*7的窗口大小进行Lee滤波，抑制斑点噪声；对Sentinel-2数据，进行去云处理；然后分别对Sentinel-1雷达数据和Sentinel-2光学数据根据友谊农场农作物分布数据进行裁剪处理；

步骤三中，玉米收割后，由成熟的玉米变为粉碎的秸秆，再进行翻土处理。在光学遥感的真彩色合成影像上，会先变亮后变暗；玉米收割后，在雷达影像上，亮度会变暗；根据这些特征，进行目视解译，选择玉米收割前与玉米收割后的样本点数据；

步骤四中，根据样本点数据，提取并分析Sentinel-1影像的VV极化方式、VH极化方式、VH/VV的后向散射系数和后向散射系数的比值，提取并分析Sentinel-2影像的B2-Blue、B3-Green、B4-Red、B8-NIR、B11-SWIR1 and B12-SWIR2的光谱反射率曲线；

步骤五中，通过分析，在光学影像上，选择B2-Blue，B3-Green，B4-Red，B8-NIR，B11-SWIR1和B12-SWIR2的光谱反射率和NDVI、NDTI、STI、NDI7、VTI、SVTI指数作为初选特征；在雷达影像上，选择VV、VH和RI-VH/VV作为初选特征；其中，Sentinel-2各个波段介绍表1所示，光谱指数公式的介绍如表2所示，其中NIR为近红外波段反射率，SWIR1为短波红外1波段反射率，SWIR2为短波红外2波段反射率。

步骤六中，通过绘制直方图和特征重要性排序，选择出区分玉米收割前和玉米收割后最优的特征集，最后选择B2-Blue)、B3-Green-、B4-Red、B8-NIR、B11-SWIR1和B12-SWIR2的光谱反射率和NDVI、NDTI、STI、NDI7、VTI、SVTI指数和VH作为特征集；

步骤七中，分别对Sentinel-1雷达影像和Sentinel-2光学影像进行阈值分割方法进行分类和随机森林分类，比较两种分类方法对于雷达影像和光学影像的适用性；

步骤八中，经过实验，雷达影像选择阈值分割方法进行玉米收割进度监测，光学影像选择随机森林分类方法进行玉米收割进度监测；

步骤九中，对分类后的影像进行去除小图斑(Majority and Clump)的分类后处理；

步骤十中，因为光谱反射率易受到外界条件，如太阳位置、传感器位置、地理位置、地形、季节、气候变化、地面湿度变化、地物本身的变异、大气状况的影响，所以加入光谱反射率作为特征，就需要对其监测结果的不确定性进行进一步消除；为了提高进一步玉米收割进度监测的精度，构造了一种时间序列的误分校正方法对其进行纠正；对于分类结果，假设收获的像素值为0，未收获的像素的值为255；基于分类结果，提取并分析了100个样本像素；获得了四种情景：全部收获，全部未收获，未收获-收获和收获-未收获，具有多重波动；修正方法侧重于第四类；以使其成为前三类之一；因为三个连续0或三个连续255的可能性很小，所以使用三个连续周期作为阈值，可以执行以下算法来修改像素值：

1、如果0或255连续出现三次或三次以上；当0在连续出现三次或三次以上时，0之后的所有值都被赋值为0；当255在连续出现三次以上，255之前的所有值均被赋值为255；然后，继续搜索连续三个0或连续三个255；

2、基于第一步的结果，再次搜索两个连续0或两个连续255；当0在一行中出现两次时，0之后的所有值都被赋值为0；当255在一行出现两次后，255之前的所有值被赋值为255；然后，继续搜索两个连续的0或两个连续255；

3、根据第二步的结果，需要找到第一次出现0的位置，并将0之后的所有值赋值为0；时间序列误分校正方法如图1所示。

步骤十一中，利用ArcGIS软件，绘制玉米收割进度的时空动态分布图，并统计累计玉米累计收割面积和收割比例，如图2所示。

表1Sentinel-2各个波段介绍

表2光谱指数介绍

/>