本发明涉及罂粟监测领域，特别涉及一种基于遥感的罂粟提取方法及系统。
背景技术：
我国面临的禁毒形势严峻复杂，毒品犯罪多发高发，禁毒工作任务十分艰巨。精确识别毒品原植物是禁毒工作的关键，而罂粟作为一种主要毒品原植物，对其进行识别至关重要。非法罂粟种植的位置一般十分隐蔽，利用传统人工方式进行排查不仅费时费力，而且效率十分低下。遥感技术的发展为大范围的罂粟种植监测提供了一条有效途径。
从遥感影像上提取罂粟的方法主要可分为监督分类与非监督分类两种。基于监督分类方法提取罂粟的研究，如联合国毒品与犯罪问题办公室(The United Nations Office On Drugs and Crime)使用监督分类方法对高空间分辨率遥感影像进行分类，从而监测主要罂粟种植区域的罂粟种植情况。该方法需要人工选取罂粟样本，其过程不仅费时费力，而且罂粟样本通常很难获得，因而限制了该方法的应用。基于非监督分类方法提取罂粟的研究，如采用一种多基元混合像元分解的非监督分类方法从高光谱影像中识别罂粟，但该方法效率很低，针对该问题又提出了一种混合调谐匹配滤波的非监督分类方法从高光谱影像中检测罂粟。现有的方法都是针对高光谱影像进行分析的，而高光谱影像虽然光谱分辨率很高，但其覆盖率很低同时幅宽小，因此采用高光谱影像实现大范围的罂粟检测也不太符合实际情况。
总体来看，以上这些方法都是基于单时相遥感影像进行分析的，未充分利用不同植被物候特征的差异进行分析，在有云的情况会导致影像成像质量较差，会导致分类结果出现很大误差，而且通常罂粟的生长期较短，而高光谱影像重返周期较长且地面覆盖率不高，会由于数据缺失导致该区域无法进行罂粟监测，因此现有的在遥感影像上提取罂粟的方法，无法实现罂粟种植的有效监测。
技术实现要素：
本发明的目的是，为了实现罂粟种植的有效监测，提供一种基于遥感的罂粟提取方法及系统。
为实现上述目的，本发明提供了如下方案：
一种基于遥感的罂粟提取方法，所述提取方法包括如下步骤：
获取罂粟不同生长阶段的研究区域的多幅遥感影像；
计算多幅所述遥感影像的特征值，所述特征值包括归一化植被指数、增强型植被指数、近红外波段反射率和红波段反射率；
将所述特征值按照研究区域的罂粟生长阶段组成时间序列影像文件，所述时间序列影像文件包括归一化植被指数时间序列影像文件、增强型植被指数时间序列影像文件、近红外波段反射率时间序列影像文件和红波段反射率时间序列影像文件；
根据所述时间序列影像文件选取少量罂粟和研究区域的其它类型植被的样本点，生成所述罂粟和其它各个类型植被对应的时序变化曲线，所述时序变化曲线包括归一化植被指数时序变化曲线、增强型植被指数时序变化曲线、近红外波段反射率时序变化曲线和红波段反射率时序变化曲线；
根据所述时序变化曲线建立罂粟提取模型；
设定罂粟提取模型的阈值，根据所述时间序列影像文件对研究区域内的罂粟进行提取，获得提取结果。
可选的，所述计算多幅所述遥感影像的特征值，所述特征值包括归一化植被指数、增强型植被指数、近红外波段反射率和红波段反射率；具体包括：
对多幅所述遥感影像进行预处理，得到预处理后的遥感影像；
提取各预处理后的遥感影像的近红外波段反射率、红波段反射率和蓝波段反射率，获得所述遥感影像的近红外波段反射率和红波段反射率；
根据所述近红外波段反射率和红波段反射率计算各预处理后的遥感影像的归一化植被指数，获得所述遥感影像的归一化植被指数；
根据所述近红外波段反射率、红波段反射率和蓝波段反射率计算各预处理后的遥感影像的增强型植被指数，获得所述遥感影像的增强型植被指数。
可选的，所述根据所述近红外波段反射率和红波段反射率计算各预处理后的遥感影像的归一化植被指数，获得所述遥感影像的归一化植被指数，具体包括：
根据所述近红外波段反射率和红波段反射率，利用公式(1)计算各预处理后的遥感影像的归一化植被指数，获得所述遥感影像的归一化植被指数；
其中，ρNIR和ρRED分别表示近红外波段反射率和红波段反射率。
可选的，所述根据所述近红外波段反射率、红波段反射率和蓝波段反射率计算各预处理后的遥感影像的增强型植被指数，获得所述遥感影像的增强型植被指数，具体包括：
根据所述近红外波段反射率、红波段反射率和蓝波段反射率，利用公式(2)计算各预处理后的遥感影像的增强型植被指数，获得所述遥感影像的增强型植被指数；
其中，ρNIR、ρRED和ρBLUE分别表示近红外波段反射率、红波段反射率和蓝波段反射率。
可选的，所述设定罂粟提取模型的阈值，根据所述时间序列影像文件对研究区域内的罂粟进行提取，获得提取结果之后，还包括：
采用小斑消除算法对所述提取结果进行后处理，将零碎图斑去除，获得最终的提取结果。
一种基于遥感的罂粟提取系统，所述提取系统包括：
遥感影像获取模块，用于获取罂粟不同生长阶段的研究区域的多幅遥感影像；
特征值计算模块，用于计算多幅所述遥感影像的特征值，所述特征值包括归一化植被指数、增强型植被指数、近红外波段反射率和红波段反射率；
时间序列影像建立模块，用于将所述特征值按照研究区域的罂粟生长阶段组成时间序列影像文件，所述时间序列影像文件包括归一化植被指数时间序列影像文件、增强型植被指数时间序列影像文件、近红外波段反射率时间序列影像文件和红波段反射率时间序列影像文件；
时序变化曲线建立模块，用于根据所述时间序列影像文件选取少量罂粟和研究区域的其它类型植被的样本点，生成所述罂粟和其它各个类型植被对应的时序变化曲线，所述时序变化曲线包括归一化植被指数时序变化曲线、增强型植被指数时序变化曲线、近红外波段反射率时序变化曲线和红波段反射率时序变化曲线；
提取模型建立模块，用于根据所述时序变化曲线建立罂粟提取模型；
提取结果获取模块，用于设定罂粟提取模型的阈值，根据所述时间序列影像文件对研究区域内的罂粟进行提取，获得提取结果。
可选的，所述特征值计算模块，具体包括：
预处理子模块，用于对多幅所述遥感影像进行预处理，得到预处理后的遥感影像；
近红外波段反射率和红波段反射率获取子模块，用于提取各预处理后的遥感影像的近红外波段反射率、红波段反射率和蓝波段反射率，获得所述遥感影像的近红外波段反射率和红波段反射率；
归一化植被指数计算子模块，用于根据所述近红外波段反射率和红波段反射率计算各预处理后的遥感影像的归一化植被指数，获得所述遥感影像的归一化植被指数；
增强型植被指数计算子模块，用于根据所述近红外波段反射率、红波段反射率和蓝波段反射率计算各预处理后的遥感影像的增强型植被指数，获得所述遥感影像的增强型植被指数。
可选的，所述归一化植被指数计算子模块，具体包括：
归一化植被指数计算单元，用于根据所述近红外波段反射率和红波段反射率，利用公式(1)计算各预处理后的遥感影像的归一化植被指数，获得所述遥感影像的归一化植被指数；
其中，ρNIR和ρRED分别表示近红外波段反射率和红波段反射率。
可选的，所述增强型植被指数计算子模块，具体包括：
增强型植被指数计算单元，用于根据所述近红外波段反射率、红波段反射率和蓝波段反射率，利用公式(2)计算各预处理后的遥感影像的增强型植被指数，获得所述遥感影像的增强型植被指数；
其中，ρNIR、ρRED和ρBLUE分别表示近红外波段反射率、红波段反射率和蓝波段反射率。
可选的，所述提取系统还包括：
后处理模块，用于采用小斑消除算法对所述提取结果进行后处理，将零碎图斑去除，获得最终的提取结果。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明公开了一种基于遥感的罂粟提取方法及系统，该方法，首先，生成时间序列影像文件，然后，通过选取少量的罂粟和研究区域的其它类型的背景植被的样本点，生成罂粟和每个类型的背景植被对应的时序变化曲线，并根据时序变化曲线建立罂粟提取模型，最后，根据所述罂粟提取模型，从所述时间序列影像文件中提取研究区域内种植的罂粟，获得提取结果。本发明通过时间序列影像对罂粟与背景植被的物候及光谱变化特征进行分析，自动提取罂粟，实现了罂粟种植的有效监测。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明提供的一种基于遥感的罂粟提取方法的流程图；
图2为本发明提供的一种基于遥感的罂粟提取系统的结构框图；
图3为本发明提供的一个具体的实施例的基于遥感的罂粟提取方法的流程图；
图4为本发明提供的一个具体的实施例的归一化植被指数时序变化曲线；
图5为本发明提供的一个具体的实时例的增强型植被指数时序变化曲线；
图6为本发明提供的一个具体的实时例的近红外波段反射率时序变化曲线；
图7为本发明提供的一个具体的实施例的红波段反射率时序变化曲线。
具体实施方式
本发明的目的是，为了实现罂粟种植的有效监测，提供一种基于遥感的罂粟提取方法及系统。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对发明作进一步详细的说明。
如图1所示，本发明提供了一种基于遥感的罂粟提取方法，所述提取方法包括如下步骤：
步骤101，获取罂粟不同生长阶段的研究区域的多幅遥感影像。
步骤102，计算多幅所述遥感影像的特征值，所述特征值包括归一化植被指数、增强型植被指数、近红外波段反射率和红波段反射率；具体包括：
对多幅所述遥感影像进行预处理，得到预处理后的遥感影像；提取各预处理后的遥感影像的近红外波段反射率、红波段反射率和蓝波段反射率，获得所述遥感影像的近红外波段反射率和红波段反射率；根据所述近红外波段反射率和红波段反射率计算各预处理后的遥感影像的归一化植被指数，获得所述遥感影像的归一化植被指数；根据所述近红外波段反射率、红波段反射率和蓝波段反射率计算各预处理后的遥感影像的增强型植被指数，获得所述遥感影像的增强型植被指数。
所述根据所述近红外波段反射率和红波段反射率计算各预处理后的遥感影像的归一化植被指数，获得所述遥感影像的归一化植被指数，具体包括：
根据所述近红外波段反射率和红波段反射率，利用公式(1)计算各预处理后的遥感影像的归一化植被指数，获得所述遥感影像的归一化植被指数；
其中，ρNIR和ρRED分别表示近红外波段反射率和红波段反射率。
所述根据所述近红外波段反射率、红波段反射率和蓝波段反射率计算各预处理后的遥感影像的增强型植被指数，获得所述遥感影像的增强型植被指数，具体包括：
根据所述近红外波段反射率、红波段反射率和蓝波段反射率，利用公式(2)计算各预处理后的遥感影像的增强型植被指数，获得所述遥感影像的增强型植被指数；
其中，ρNIR、ρRED和ρBLUE分别表示近红外波段反射率、红波段反射率和蓝波段反射率。
步骤103，将所述特征值按照研究区域的罂粟生长阶段组成时间序列影像文件，所述时间序列影像文件包括归一化植被指数时间序列影像文件、增强型植被指数时间序列影像文件、近红外波段反射率时间序列影像文件和红波段反射率时间序列影像文件。
步骤104，根据所述时间序列影像文件选取少量罂粟和研究区域的其它类型植被的样本点，生成所述罂粟和其它各个类型植被对应的时序变化曲线，所述时序变化曲线包括归一化植被指数时序变化曲线、增强型植被指数时序变化曲线、近红外波段反射率时序变化曲线和红波段反射率时序变化曲线。
步骤105，根据所述时序变化曲线建立罂粟提取模型。
步骤106，设定罂粟提取模型的阈值，根据所述时间序列影像文件对研究区域内的罂粟进行提取，获得提取结果。
步骤106所述设定罂粟提取模型的阈值，根据所述时间序列影像文件对研究区域内的罂粟进行提取，获得提取结果之后，还包括：采用小斑消除算法对所述提取结果进行后处理，将零碎图斑去除，获得最终的提取结果。
如图2所示，本发明还提供一种基于遥感的罂粟提取系统，所述提取系统包括：
遥感影像获取模块201，用于获取罂粟不同生长阶段的研究区域的多幅遥感影像；
特征值计算模块202，用于计算多幅所述遥感影像的特征值，所述特征值包括归一化植被指数、增强型植被指数、近红外波段反射率和红波段反射率；所述特征值计算模块202，具体包括：预处理子模块，用于对多幅所述遥感影像进行预处理，得到预处理后的遥感影像；近红外波段反射率和红波段反射率获取子模块，用于提取各预处理后的遥感影像的近红外波段反射率、红波段反射率和蓝波段反射率，获得所述遥感影像的近红外波段反射率和红波段反射率；归一化植被指数计算子模块，用于根据所述红外波段反射率和红波段反射率计算各预处理后的遥感影像的归一化植被指数，获得所述遥感影像的归一化植被指数；增强型植被指数计算子模块，用于根据所述近红外波段反射率、红波段反射率和蓝波段反射率计算各预处理后的遥感影像的增强型植被指数，获得所述遥感影像的增强型植被指数。
所述归一化植被指数计算子模块，具体包括：归一化植被指数计算单元，用于根据所述近红外波段反射率和红波段反射率，利用公式(1)计算各预处理后的遥感影像的归一化植被指数，获得所述遥感影像的归一化植被指数。
其中，ρNIR和ρRED分别表示近红外波段反射率和红波段反射率。
所述增强型植被指数计算子模块，具体包括：增强型植被指数计算单元，用于根据所述近红外波段反射率、红波段反射率和蓝波段反射率，利用公式(2)计算各预处理后的遥感影像的增强型植被指数，获得所述遥感影像的增强型植被指数。
其中，ρNIR、ρRED和ρBLUE分别表示近红外波段反射率、红波段反射率和蓝波段反射率。
时间序列影像建立模块203，用于将所述特征值按照研究区域的罂粟生长阶段组成时间序列影像文件，所述时间序列影像文件包括归一化植被指数时间序列影像文件、增强型植被指数时间序列影像文件、近红外波段反射率时间序列影像文件和红波段反射率时间序列影像文件。
时序变化曲线建立模块204，用于根据所述时间序列影像文件选取少量罂粟和研究区域的其它类型植被的样本点，生成所述罂粟和其它各个类型植被对应的时序变化曲线，所述时序变化曲线包括归一化植被指数时序变化曲线、增强型植被指数时序变化曲线、近红外波段反射率时序变化曲线和红波段反射率时序变化曲线。
提取模型建立模块205，用于根据所述时序变化曲线建立罂粟提取模型。
提取结果获取模块206，设定罂粟提取模型的阈值，根据所述时间序列影像文件对研究区域内的罂粟进行提取，获得提取结果。
所述提取系统还包括：后处理模块，用于采用小斑消除算法对所述提取结果进行后处理，将零碎图斑去除，获得最终的提取结果。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明公开了一种基于遥感的罂粟提取方法及系统，该方法首先生成时间序列影像文件，然后通过选取少量的罂粟和研究区域的其它类型的背景植被的样本点，生成罂粟和每个类型的背景植被对应的时序变化曲线，并根据时序变化曲线建立罂粟提取模型，最后，根据所述罂粟提取模型，从所述时间序列影像文件中提取研究区域内种植的罂粟，获得提取结果。本发明通过时间序列影像对罂粟与背景植被的物候及光谱变化特征进行分析，自动提取罂粟，实现了罂粟种植的有效监测。
结合一个示例性实施例对本发明技术方案进行详细说明，该实施例的具体提取方法流程图，如图3所示。研究区域位于甘肃省金昌市国家罂粟种植基地，通过调研及查阅相关文献得知，该区域罂粟播种时间为3月末，出苗时间为5月初，开花时间为6月末至7月初，收割时间为7月中下旬。
步骤1，根据该区域罂粟的生长阶段分为出苗、拔节、开花、结果、收割，提取该区域的不同生长阶段的遥感影像，获取Landsat8OLI影像，具体如表1所示：
表1
步骤2：对各景Landsat8OLI影像进行辐射定标、大气校正、几何校正等预处理操作。
步骤3：对预处理后的各景影像计算归一化植被指数NDVI(Normalized Differential Vegetation Index)，计算公式如下所示：
其中，ρNIR和ρRED分别表示近红外波段反射率和红波段反射率。
归一化植被指数(NDVI)：反映土地覆盖植被状况的一种遥感指标，定义为近红外波段与红波段反射率之差与之和的商。
步骤4：对预处理后的各景影像计算增强型植被指数EVI(Enhanced Vegetation Index)，计算公式如下所示：
其中，ρNIR、ρRED和ρBLUE分别表示近红外波段反射率、红波段反射率和蓝波段反射率。
增强型植被指数(EVI)：通过加入蓝波段以增强植被信号，矫正土壤背景和气溶胶散射的影响。
步骤5：将影像的第N1、N2、N3、N4、N5时相所对应的NDVI按顺序组合成归一化植被指数时间序列影像文件，即N1、N2、N3、N4、N5时相的NDVI分别对应时间序列影像文件的波段1、波段2、波段3、波段4、波段5。接着按照上述过程分别组合生成增强型植被指数时间序列影像文件、红波段时间序列影像文件、近红外波段时间序列影像文件。
步骤6：分析研究区域内包含的主要植被类型，然后在4种时间序列影像文件上分别选取少量罂粟与各植被类型的样本点生成对应罂粟与各植被类型的NDVI时序变化曲线、EVI时序变化曲线、红波段反射率时序变化曲线以及近红外波段反射率时序变化曲线。
经分析，罂粟生长阶段内研究区域主要种植作物包括小麦、大麦、苜蓿等8种类型，其次还包括菜地等零散分布的农作物。根据高分辨率遥感影像在各时间序列影像文件上选取罂粟与各作物类型的样本点，然后绘制与其对应的时序变化曲线，结果如图4-7所示。
步骤7：对罂粟与其他植被的时序变化曲线进行分析，找出可分性最大的特征，然后对该特征中各时相对应特征值进行数学运算构建罂粟提取模型(如两个时相对应特征的比值；两个时相对应特征值大小的比较；单个时相对应特征值的大小等)，最后对罂粟提取模型设定阈值，实现罂粟的提取。
对图4-7进行分析可知，NDVI时序变化曲线与EVI时序变化曲线中poppy、crop1、crop2、crop4在N1至N2之间增长迅速，其他几种作物增长缓慢，并且crop1在N2至N3之间减少，因此可利用该特征将poppy、crop2、crop4与其他作物进行区分。在近红外波段反射率时序变化曲线与红波段反射率时序变化曲线中crop4与poppy、crop2在N4至N5内变化趋势相反，同时poppy与crop2在N5时段反射率不同，并且红波段反射率时序变化曲线中poppy与crop2的差异比近红外波段反射率时序变化曲线的差异大，因此可以利用红波段反射率时序变化曲线对poppy、crop2与crop4进行区分。因为NDVI在植被生长旺盛期容易达到饱和，而EVI则能够克服这一现象，因此选用EVI时序变化曲线结合红波段反射率时序变化曲线对罂粟提取模型进行研究。经过对阈值的反复试验，最终确定利用EVI时序变化曲线中(EVIN2/EVIN1)>5且EVIN3>EVIN2对poppy、crop2、crop4与其他农作物进行区分，同时利用红波段反射率时序变化曲线中ref(R)N4<ref(R)N5且ref(R)N5<0.12对poppy、crop2、crop4进行区分，即罂粟提取模型为(EVIN2/EVIN1)>5&EVIN3>EVIN2&ref(R)N4<ref(R)N5&ref(R)N5<0.12。
步骤8：采用小斑消除算法对罂粟提取结果进行后处理，从而将零碎图斑去除获得最终的罂粟提取结果。
本发明采用时间序列遥感影像对罂粟进行提取，可以综合考虑罂粟与背景植被的光谱变化特征及物候特征差异，从而提高罂粟提取的精度；通过对区分性大的特征对应的时序变化曲线建立罂粟提取模型，然后对罂粟提取模型设定阈值实现罂粟的自动提取，该过程不需要分类，只需通过简单的波段运算然后进行阈值分割即可实现罂粟的提取，因此效率大幅提高。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。