基于哨兵数据和类胡萝卜素时序特征的玉米制图方法

1.本发明涉及农作物生长监测技术领域，具体涉及一种基于哨兵数据和类胡萝卜素时序特征的玉米制图方法。


背景技术：

2.快速准确获取大范围农作物种植分布数据，对于确保粮食安全至关重要。传统的农业抽样调查方法，数据获取成本高、耗时长并且难以避免实现研究区域全覆盖。遥感影像具有时效性强、覆盖范围大的特点，在农田调查特别是农作物分布信息提取中越来越发挥重要作用。但由于不同类型的农作物光谱具有一定的相似性，而同一种农作物光谱在不同区域和不同年份具有很强的异质性，给常规基于光谱的农作物制图方法带来挑战。基于时序遥感影像的农作物制图已经成为农作物遥感的前沿发展方向。
3.玉米作为全球种植范围最广、产量最高的粮食作物，在全球农业生产中具有重要地位。快速准确获取大范围的玉米种植面积及其空间分布，对于科学预测全球粮食安全态势具有重要作用。由于农作物生长期普遍均在半年以内，因此对农作物遥感监测所需的时序遥感数据提出了较高的要求。modis遥感数据具有每日覆盖全球、时效性高的优势，因此在大尺度农作物遥感监测中得到很好的应用。在农作物遥感监测领域，已有学者提出了基于modis时序遥感数据的系列大宗农作物遥感监测技术方法。如基于生长盛期nmdi增减比值指数的玉米自动制图方法以及基于分蘖-抽穗关键物候期水体与植被指数变化比值指数的水稻制图方法。这些农作物制图方法通过巧妙锁定关键物候期提取不同农作物独特的时序变化特征，能有效地实现大范围多年分农作物分布信息自动提取。但不足之处在于modis数据的分辨率较低，混合像元问题比较突出。
4.sentinel-2数据，具有10米空间分辨率、5天时间分辨率，并且免费提供公众使用。哨兵数据特有的三个红光边沿波段，在农作物遥感监测方面具有巨大潜力。但哨兵数据作为一种较为新型的遥感时序数据，基于哨兵时序数据的农作物制图技术方法相对匮乏。如果完全照搬借用基于modis时序影像的农作物制图方法，其不足之处在于：其一，哨兵影像数据时间分辨率毕竟不如modis时序数据，特定关键物候期如生长盛期哨兵数据可能仅仅少数几期影像，难以全面刻画描述特定农作物特征；其二，未能很好地利用哨兵影像数据的优势，未能有效地利用哨兵红光边沿波段信息，充分发挥其在农作物制图领域的潜力。


技术实现要素：

5.有鉴于此，本发明的目的在于提供一种基于哨兵数据和类胡萝卜素时序特征的玉米制图方法，克服中高分辨率遥感影像数据可获得性有限导致的农作物制图方法难以大范围推广应用的技术瓶颈，确保玉米制图方法的适用性与鲁棒性。
6.为实现上述目的，本发明采用如下技术方案：一种基于哨兵数据和类胡萝卜素时序特征的玉米制图方法，包括以下步骤：步骤s1:基于哨兵波段的反射率数据，构建研究区植被指数、类胡萝卜素、叶绿素
指数时序数据集；步骤s2:根据得到的研究区植被指数时序数据集，逐像元获取农作物生长期；步骤s3:根据生长期类胡萝卜素指数时序数据集,构建农作物生长期类胡萝卜素变化特征指标;步骤s4:根据生长期叶绿素指数时序数据集,提取农作物生长期叶绿素极值指标；步骤s5:结合农作物生长期类胡萝卜素变化特征指标、生长期叶绿素指数极值指标，构建玉米制图判别模型；步骤s6:根据玉米制图判别模型，获得研究区玉米空间分布图。
7.进一步的，所述步骤s1具体为：步骤s11:基于哨兵波段的反射率数据,按照时间顺序逐景依次计算植被指数、类胡萝卜素、叶绿素指数，获得该年份研究区植被指数、类胡萝卜素、叶绿素指数时序数据集;步骤s12:基于无云日的研究区植被指数、类胡萝卜素、叶绿素指数时序数据集，采用whittaker smoother数据平滑方法，逐像元构建研究年份内研究区逐日连续平滑的植被指数、类胡萝卜素、叶绿素指数时序数据集。
8.进一步的，所述步骤s2具体为：步骤s21:根据植被指数evi2年内时序曲线局部最大值，确定农作物生长峰值期;步骤s22:逐像元求算evi2时序曲线中数值大于预设值的局部最大值，如果同时找到两个或两个以上的局部evi2最大值，依据相邻两个局部最大值的时间间隔进一步筛选，保留相邻间隔大于预设天数以上的evi2局部最大值；步骤s23:设置约束条件，在出现多个evi2局部最大值的情况下，删除与邻近局部最小值差异小于预设值的局部最大值,获得满足条件的evi2局部最大值；步骤s24:将这些evi2局部最大值记录为pn，其中n表示n季节作物，evi2局部最大值pn所在时刻，分别为该生长期所对应的农作物生长峰值期；步骤s25:根据农作物生长峰值期，提取农作物生长期。
9.进一步的，所述步骤s25具体为：预设生长峰值期前50天定义为农作物生长开始期，生长峰值期后30天定义为农作物生长结束期，而农作物生长开始期至生长结束期的这段时间，记录为农作物生长期。
10.进一步的，所述步骤s3具体为:步骤s31:逐像元获取农作物生长期类胡萝卜素指数cri时序数据集，分别计算农作物生长期类胡萝卜素指数cri时序数据集的最大值cmax、最小值cmin以及标准差csd;步骤s32:根据下式，计算得到农作物生长期类胡萝卜素变化特征指标criticriti=(cmax-cmin)*(csd）。
11.进一步的，所述步骤s4具体为:逐像元获取农作物生长期叶绿素指数mtci时序数据集，计算农作物生长期叶绿素指数mtci极值指标，将农作物生长期叶绿素指数mtci的最大值，记录为mmax。
12.进一步的，所述步骤s5具体为： 根据农作物生长期类胡萝卜素变化特征指标criti、生长期叶绿素指数mtci极值指标mmax，建立玉米制图判别模型，判别条件为：如果criti》ω1 & mmax》ω2，则该像元为玉米，否则为其他农作物。
13.本发明与现有技术相比具有以下有益效果：
1、本发明通过挖掘基于红光边沿波段构建的植被色素指数在农作物生长期内的时序特征，构建基于农作物生长期色素时序特征的光谱指数，有助于摆脱仅仅依赖植被指数时序数据难以应对光谱类内异质性问题；2、本发明通过综合整个生长期农作物植被色素的变异性特征，而非某期或几期影像光谱特征，有助于克服中高分辨率遥感影像数据可获得性有限导致的农作物制图方法难以大范围推广应用的技术瓶颈，确保玉米制图方法的适用性与鲁棒性。
附图说明
14.图1是本发明方法流程图；图2是本发明一实施例中玉米、大豆的evi2、cri指数的时序信号图；图3是本发明一实施例中玉米、大豆的evi2、mtci指数的时序信号图；图4是本发明一实施例中农作物生长期示意图；图5是本发明一实施例中农作物生长前期类胡萝卜素变化特征指标示意图；图6是本发明一实施例中农作物生长期叶绿素指数mtci极值指标示意图；图7是本发明一实施例中玉米制图技术流程方法；图8是本发明一实施例中研究区玉米空间分布图。
具体实施方式
15.下面结合附图及实施例对本发明做进一步说明。
16.请参照图1，本发明提供一种基于哨兵数据和类胡萝卜素时序特征的玉米制图方法，包括以下步骤：步骤s1:基于哨兵波段的反射率数据，构建研究区植被指数、类胡萝卜素、叶绿素指数时序数据集；步骤s2:根据得到的研究区植被指数时序数据集，逐像元获取农作物生长期；步骤s3:根据生长期类胡萝卜素指数时序数据集,构建农作物生长期类胡萝卜素变化特征指标;步骤s4:根据生长期叶绿素指数时序数据集,提取农作物生长期叶绿素极值指标；步骤s5:结合农作物生长期类胡萝卜素变化特征指标、生长期叶绿素指数极值指标，构建玉米制图判别模型；步骤s6:根据玉米制图判别模型，获得研究区玉米空间分布图。
17.在本实施例中，基于哨兵2号（sentinel-2）的反射率数据，具有宽幅和高时空分辨率以及免费共享等多重优势。
18.在本实施例中，所述步骤s1具体为：步骤s11:基于哨兵波段的反射率数据,按照时间顺序逐景依次计算植被指数、类胡萝卜素、叶绿素指数，获得该年份研究区植被指数、类胡萝卜素、叶绿素指数时序数据集;步骤s12:基于无云日的研究区植被指数、类胡萝卜素、叶绿素指数时序数据集，采用whittaker smoother数据平滑方法，逐像元构建研究年份内研究区逐日连续平滑的植被指数、类胡萝卜素、叶绿素指数时序数据集。参考图2，为构建的玉米、大豆的evi2、cri指数的时序信号图。所构建的玉米、大豆的evi2、mtci指数的时序信号图，若图3所示。
[0019] 优选的，植被指数是表征植被生长状态以及空间分布密度的因子。常见的植被指数有ndvi和evi2。ndvi为归一化植被指数，全称为normalized difference vegetation index。evi2为增强型植被指数，全称为enhanced vegetation index。evi2指数的计算公式为：其中ρ
red
, ρ
nir
分别为哨兵影像的红光、近红外波段的反射率，红光波段、近红外波段的范围分别为650-680nm、785-899nm。
[0020]
类胡萝卜素指数（carotenoid reflectance index，cri），类胡萝卜素是一种植物中广泛存在的重要色素。类胡萝卜素含量及其变化能很好地反应不同农作物生长特性及其生长变化规律。类胡萝卜素指数的计算公式为：其中，表示哨兵影像蓝光、第一红光边缘波段的反射率，蓝光波段。
[0021]
叶绿素指数 (meris terrestrial chlorophyll index, mtci)，叶绿素是植物进行光合作用的重要色色。叶绿素吸收大部分的红光和紫光但反射绿光，在光合作用的光吸收中起核心作用。叶绿素指数mtci的计算公式为：其中，表示哨兵影像红光、第一红光边缘波段、第二红光边缘波段的反射率。
[0022]
在本实施例中，所述步骤s2具体为：步骤s21:根据植被指数evi2年内时序曲线局部最大值，确定农作物生长峰值期;步骤s22:逐像元求算evi2时序曲线中数值大于预设值的局部最大值，如果同时找到两个或两个以上的局部evi2最大值，依据相邻两个局部最大值的时间间隔进一步筛选，，保留相邻间隔大于80天以上的evi2局部最大值；步骤s23:设置约束条件，在出现多个evi2局部最大值的情况下，删除与邻近局部最小值差异小于0.15的局部最大值,获得满足条件的evi2局部最大值；步骤s24:将这些evi2局部最大值记录为pn，其中n表示n季节作物，n可以取值为1，2，3，分别对应单季作物、双季作物以及三季作物，evi2局部最大值pn所在时刻，分别为该生长期所对应的农作物生长峰值期；步骤s25:根据农作物生长峰值期，提取农作物生长期。
[0023]
优选的，预设生长峰值期前50天定义为农作物生长开始期，生长峰值期后30天定义为农作物生长结束期，而农作物生长开始期至生长结束期的这段时间，记录为农作物生长期。
[0024]
参考图5，在本实施例中，所述步骤s3具体为:步骤s31:逐像元获取农作物生长期类胡萝卜素指数cri时序数据集，分别计算农
作物生长期类胡萝卜素指数cri时序数据集的最大值cmax、最小值cmin以及标准差csd;步骤s32:根据下式，计算得到农作物生长期类胡萝卜素变化特征指标criticriti=(cmax-cmin)*(csd）。
[0025]
参考图6，在实施例中，步骤s4具体为:逐像元获取农作物生长期叶绿素指数mtci时序数据集，计算农作物生长期叶绿素指数mtci极值指标，将农作物生长期叶绿素指数mtci的最大值，记录为mmax。
[0026]
在本实施例中，所述步骤s5具体为： 根据农作物生长期类胡萝卜素变化特征指标criti、生长期叶绿素指数mtci极值指标mmax，建立玉米制图判别模型，判别条件为：如果criti》ω1 & mmax》ω2，则该像元为玉米，否则为其他农作物。
[0027]
以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。