CN106918816A - 农作物长势监测系统及其数据处理方法和模块 - Google Patents

Abstract

本发明公开了农作物长势监测系统及其数据处理方法和模块，涉及遥感图像处理和农业监测技术领域。所述数据处理方法包括：获取监测用遥感数据和修正用遥感数据；分别对所述监测用遥感数据和修正用遥感数据进行预处理，得到对应的植被指数数据；根据修正用植被指数数据获得修正参数；根据所述修正参数修正所述监测用植被指数数据。根据本发明所述修正后的植被指数数据来监测农作物长势，削弱了由于作物轮作和耕地种植区年际变化造成作物长势监测结果的不确定性，提高了农作物长势监测的精准度。

Description

农作物长势监测系统及其数据处理方法和模块

技术领域

本发明涉及遥感图像处理和农业监测技术领域，特别涉及一种农作物长势监测系统及其数据处理方法和模块。

背景技术

近年，全球各国，无论是发展中国家还是发达国家，无一例外的越来越多的依赖及时、准确、透明的农情监测结果，以应对在全球尺度、地区尺度、国家尺度以及省州级尺度所面临的粮食安全问题。研究表明，与作物收割之后统计获得的产量数据相比较，在作物生长阶段尽早的获取作物生长状况信息，更有利于应对突发的粮食供应短缺以及粮食进出口相关政策的制定。当粮食产量预期可能出现下降的情况下，稳定国际粮价对于粮食进口国和出口国均具有重要意义。

目前许多国家或组织纷纷构建了各自的农业监测系统，用以实现国家尺度、地区尺度或是全球尺度的作物生长形势的监测，以及产量变化形势的预测。尽管这些数量众多的农业监测系统采用监测方法可能各有相同，但作物长势监测始终是监测系统中最核心的模块之一。

大多数现有的农情监测系统均采用某种遥感指数来实现作物生长形势的监测。在这些遥感指数中，归一化植被指数(The Normalized Difference Vegetation Index，简称NDVI)，因为其能够直观反映出作物光合作用的强弱，并在一定程度上削弱植物冠层的背景影响，如土壤、潮湿地面、雪、枯叶、粗糙度等，同时又能在一定程度上反映出植被是否受到水分、温度等外在条件的胁迫，是最常用的遥感指数之一。

根据获得的遥感图像数据，利用年度间遥感指数的变化，分析人员一方面能够将当前时期的遥感指数与其他年份同期的数据进行对比，另一方面也可以对当前生长季的遥感指数的动态变化与其他年份遥感指数的动态变化进行对比，进而评估作物生长状况并在早期进行产量预测。

根据前述的对比获得了作物长势的监测结果，所述的监测结果既包含真实的作物长势信息，也包含由作物物候、作物轮作等因素造成的伪变化信息。另外，年际间耕地上同一像元播种情况的空间差异更是直接反应在实时作物长势监测和作物生长过程监测结果上。

然而，现有的农情监测系统采用的长势监测方法均没有考虑作物轮作、作物播种状况变化造成的耕地利用状况的变化，这些年际间的作物轮作造成的耕地利用状况的变化均会在作物长势监测的过程中引入很大的不确定性。当作物轮作以及耕地种植区发生显著变化时，作物长势信息将会失去原有的指示意义。

发明内容

针对现有技术中的不足，本发明提供了一种农作物长势监测系统及其数据处理方法和模块，用以消除耕地利用状况的变化在作物长势监测过程中引入的不确定性，从而获得更加准确的作物长势监测结果。

为了解决上述技术问题，本发明提供了一种农作物长势监测数据处理方法，其中，包括：

获取监测用遥感数据和修正用遥感数据；

分别对所述监测用遥感数据和修正用遥感数据进行预处理，得到对应的植被指数数据；

根据修正用植被指数数据获得修正参数；

根据所述修正参数修正所述监测用植被指数数据。

优选地，所述根据修正用植被指数数据获得修正参数的步骤包括：

选择识别用关键时期；

基于所述修正用植被指数数据，在所述识别用关键时期，识别各像元对应的耕地为未种植耕地或种植耕地；

对所述修正用遥感数据按照监测用遥感数据的图像分辨率重新采样；

获取重新采样后各像元所覆盖的采样前未种植耕地和种植耕地的像元数量，根据公式(2-1)计算得到重新采样后各像元中未种植耕地的比例，所述未种植耕地的比例为所述的修正参数；

优选地，所述选择识别用关键时期的步骤包括：

采用类间Jeffries-Matusita距离对所述修正用植被指数数据进行量化；

对应于多个不同的两个时期的植被指数组合，分别计算对应的未种植耕地和种植耕地两种不同类别的Jeffries-Matusita距离；

比较多个所述Jeffries-Matusita距离的大小，获取Jeffries-Matusita距离最大的植被指数组合；

分别获取Jeffries-Matusita距离最大的植被指数组合对应的两个时期，将其作为用于识别的关键时期。

优选地，基于所述修正用植被指数数据，在所述识别用关键时期，识别各像元对应的耕地为未种植耕地或种植耕地的步骤包括：

读取与两个所述关键时期对应的植被指数值；

计算所述两个植被指数值的差值，并将植被指数值大的关键时期设为高峰时期；

逐像元比较高峰时期的植被指数值与预设的种植区提取阈值的大小，如果高峰时期的植被指数值大于或等于所述种植区提取阈值，则将所述像元对应的耕地识别为种植耕地；

如果高峰时期的植被指数值小于所述种植区提取阈值，逐像元比较所述两个所述关键时期对应的植被指数值的差值与预设的未种植区提取阈值的大小，如果所述差值大于或等于所述未种植区提取阈值，则将所述像元对应的耕地识别为种植耕地，如果所述差值小于所述未种植区提取阈值，则将所述像元对应的耕地识别为未种植耕地。

优选地，根据所述修正参数修正所述监测用植被指数数据的步骤包括：

根据红波段和近红外波段的反射率在种植作物端元和未种植作物端元的关系，获得以下二元线性方程：

其中，所述ρ_R和ρ_NIR为MODIS像元尺度的端元红波段与近红外波段反射率；

所述UALR为对应像元的未种植耕地比例；

ρ_Rcropped为种植作物耕地端元的红波段反射率；

ρ_NIRcropped为种植作物耕地端元的近红波段反射率；

ρ_Runcropped为未种植作物耕地端元的红波段反射率；

ρ_NIRuncropped为未种植作物耕地端元的近红外波段反射率；

其中，未种植作物耕地端元的红波段反射率ρ_Runcropped、近红外波段反射率ρ_Runcropped通过统计未种植耕地比例大于98％的像元对应波段的均值确定；

根据公式(2-2)获得种植作物耕地端元的红波段反射率ρ_Rcropped和近红外波段反射率ρ_NIRcropped；

根据公式(2-3)逐像元计算出对应的植被指数值，从而获得修正后用于农作物长势监测的植被指数数据；

本发明还提供了一种农作物长势监测数据处理模块，其中，包括：

数据获取子模块，用于获取监测用遥感数据和修正用遥感数据，并对监测用遥感数据和修正用遥感数据进行预处理，获得对应的植被指数数据；

数据预处理子模块，用于分别对所述监测用遥感数据和修正用遥感数据进行预处理，得到对应的植被指数数据；

修正参数计算子模块，用于根据修正用遥感数据获得修正参数；和

数据修正子模块，用于根据所述修正参数修正所述监测用植被指数数据。

优选地，所述修正参数计算子模块包括：

关键时期选择单元，用于选择识用别关键时期；

耕地种植状况识别单元，用于基于所述修正用植被指数数据，在所述识别用关键时期，识别各像元对应的耕地为未种植耕地或种植耕地；

采样单元，用于将修正用遥感数据按照监测用遥感数据的图像分辨率重新采样，并获取重新采样后各像元所覆盖的采样前未种植耕地和种植耕地的像元数量；

修正参数计算单元，用于根据重新采样后各像元所覆盖的采样前未种植耕地和种植耕地的像元数量计算得到未种植耕地的比例，所述未种植耕地的比例为所述的修正参数。

优选地，所述数据修正子模块包括：

第一计算单元，用于根据公式(2-2)提供的二元线性方程，及未种植作物的耕地端元的红波段反射率ρ_Runcropped和近红外波段反射率ρ_Runcropped计算得到种植作物的耕地端元的红波段反射率ρ_Rcropped和近红外波段反射率ρ_NIRcropped；

其中，所述ρ_R和ρ_NIR为MODIS像元尺度的红波段与近红外波段反射率；

所述UALR为对应像元的未种植耕地比例；

和

第二计算单元，用于根据公式(2-3)，针对监测用植被指数数据，逐像元计算出对应的植被指数值，从而获得修正后用于农作物长势监测的植被指数数据；

本发明还提供了一种农作物长势监测系统，其中，包括前面所述的数据处理模块和监测模块，所述监测模块以所述数据处理模块修正后的植被指数数据作为监测数据，进行作物长势的过程监测和/或作物长势的实时监测。

优选地，所述农作物长势监测系统还包括数据输出模块，其与所述监测模块相连接，用于输出监测结果数据。

本发明以高分辨率遥感数据为数据源，识别不同年份间耕地利用状况，并计算未种植耕地的比例。通过所述未种植耕地比例来修正低分辨率遥感植被指数数据，从而还原耕地种植区的低分辨率遥感植被指数的真实数据，削弱了由于作物轮作和耕地种植区年际变化造成的农作物长势监测不确定性，因而在开展年度间作物长势的对比分析及实时作物长势监测时，提高了耕地种植状况动态变化情况下作物长势监测的精准度。

附图说明

通过参照以下附图对本发明实施例的描述，本发明的上述以及其它目的、特征和优点将更为清楚，在附图中：

图1为本发明实施例一提供的一种农作物长势监测数据处理方法流程示意图；

图2为本发明实施例一中获得修正参数的流程示意图；

图3为本发明实施例一利用植被指数阈值建立决策树、在关键时期识别各像元农田种植状况的具体流程示意图；

图4为本发明实施例二所述农作物长势监测数据处理模块的结构原理框图；

图5为本发明实施例三所述农作物长势监测系统的结构原理框图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。附图中的流程图、框图图示了本发明实施例的系统、方法、装置的可能的体系框架、功能和操作，流程图和框图上的方框可以代表一个模块、程序段或仅仅是一段代码，所述模块、程序段和代码都是用来实现规定逻辑功能的可执行指令。也应当注意，所述实现规定逻辑功能的可执行指令可以重新组合，从而生成新的模块和程序段。因此附图的方框以及方框顺序只是用来更好的图示实施例的过程和步骤，而不应以此作为对发明本身的限制。

实施例一

参见图1，为本发明实施例一提供的一种农作物长势监测数据处理方法的流程示意图。所述方法包括：

步骤S1，获取监测用遥感数据和修正用遥感数据。

步骤S2，对所述监测用遥感数据和修正用遥感数据进行预处理，得到对应的植被指数数据。其中，所述植被指数数据包括由多个特定时期顺序组成的时间序列植被指数数据。

步骤S3，根据修正用植被指数数据获得修正参数。

步骤S4，根据所述修正参数修正所述监测用植被指数数据。

具体地，在步骤S1中，所述监测用遥感数据和修正用遥感数据分别采用不同分辨率的遥感数据。其中，监测用遥感数据为中低分辨率的遥感数据，来源于中分辨率成像光谱仪(MODerate resolution Imaging Spectro-radiometer，简称MODIS)，具体为MODIS上午星(Terra)和下午星(Aqua)16日合成归一化植被指数(The Normalized DifferenceVegetation Index，简称NDVI)产品。例如采用产品编号为MOD13Q1和MYD13Q1的三级(L3)产品，其空间分辨率为250m，时间范围为2009-2010年度和2010-2011年度冬季作物完整生长期，共包括五个不同影像块(h26v04，h26v05，h27v04，h27v05和h28v05)的时间序列数据。所述数据是从美国国家宇航局(NASA)的Reverb网络工具(http://reverb.echo.nasa.gov)上下载获取。

所述的修正用遥感数据为高分辨率的遥感数据。例如，在本实施例中，依据研究区内的主要作物的物候期，并结合地面观测时间，申请2010年2月中旬至2010年5月下旬、2011年2月中旬至2011年5月下旬每半月一期的环境1A、1B星(HJ-1A/1B)的CCD数据，研究区需要7景环境星以实现全覆盖，本实施例共获取了112景数据质量较好的HJ-1A/B CCD数据，其分辨率为30m。

在步骤S2中，分别对前述两种遥感数据进行预处理。对于中低分辨率的检测用遥感数据，其预处理过程包括数据拼接、重投影、数据类型转换和数据格式转换等过程。例如，利用MRT软件(lpdaac.usgs.gov/tools/modis_reprojection_tool)完成重投影，投影参数使用统一的Albers110大地坐标系。然后，根据产品提供的标定参数，提取植被指数产品中的NDVI数据，或者是EVI数据。再将数据换算成真实植被指数，并转存成GEO-TIFF格式，利用研究区的矢量边界对图像进行裁剪。最后得到监测用遥感图像和植被指数数据，所述植被指数数据是指由多个特定时期顺序组成的时间序列植被指数数据。

对于修正用遥感数据，即HJ-1CCD数据，在使用之前，需要进行预处理，主要包括几何精校正和辐射校正，而大气纠正是辐射校正中重要的一步，也是影响冠层反射率质量、植被指数计算及陆表类型分类精度的重要环节，因此本发明将辐射校正分为两个步骤进行，第一步将DN值转换为辐亮度，第二步进行大气校正，具体过程如下：

根据每景环境星影像自带的配置文件中标明的各波段绝对定标系数(包含在*.xml文件中)将CCD图像各波段的DN值转换为辐亮度图像，其计算公式如(1-1)所示：

式中A为绝对定标系数增益，L₀为绝对定标系数偏移量，转换后辐亮度单位为W·m^-2·sr^-1·μm^-1；

本发明通过ENVI软件提供的FLAASH(FAST Line-of-sight AtmosphericAnalysis of Spectral Hypercube)大气纠正模块对影像进行大气纠正。在FLAASH模型中，依据影像获取时间将大气模型选择夏天或冬天中纬度模型，气溶胶模型则选择农村模型(Rural)；由于空气能见度与气溶胶厚度存在着很好的转换关系，利用转换公式(1-2)(1-3)将气溶胶厚度转换成FLAASH所需要的能见度。气溶胶厚度则来自于MODIS的气溶胶产品(MOD04_V5.1，http://ladsweb.nascom.nasa.gov/data/search.html)。

春夏：

秋冬：

式中，t₅₅₀表示550nm波长范围的气溶胶光学厚度，V代表水平气象视距(km)。

关于环境星的几何纠正，通常需要人工选取大量的地面控制点进行校正。由于本实施例中所用到的环境星数量众多，如果全部采用人工方式进行处理将耗费大量的时间和精力。因此本实施例采用基于AROP程序包的类Landsat遥感影像配准与正射纠正软件包基础上，升级实现了环境星CCD数据并行化预处理系统进行，并对HJ-1CCD数据进行几何校正和正射校正，其中参考影像为预处理后的Landsat TM数据，DEM数据为90m分辨率的SRTM数据(https://lta.cr.usgs.gov/SRTM2)。系统输出的经过处理的影像采用UTM投影，WGS84椭球体，UTM投影参数是北半球50区。几何纠正的平均误差基本控制在1个像元内，部分山区地带的影像几何偏差稍大。

经过几何纠正、正射校正、辐射定标和大气纠正之后，就获取了各景影像蓝波段、绿波段、红波段和近红外波段的反射率，利用植被指数计算公式可以计算各景影像对应的植被指数(NDVI)数值，并获取了各时期影像对应的NDVI影像。

在步骤3中，根据修正用植被指数数据获得修正参数的过程具体如图2所示，首先，在步骤S31，选择识别用关键时期，所述识别用关键时期为确定修正参数时用到的关键时期。具体地选择的原理及过程如下：

在获得时间序列植被指数数据后，将区分未种植耕地与种植耕地的关键时期选择问题转换为最佳波段组合问题，本发明采用Jeffries-Matusita距离(简称J-M距离)以指示未种植耕地和种植耕地两种类别间可分性，基于修正用遥感数据，采用未种植耕地和种植耕地两种类别间J-M距离进行量化。

将监测用时间序列植被指数数据的不同时期的数据视为不同波段的遥感数据，针对多个不同的两个时期波段组合，分别计算对应的未种植耕地和种植耕地两种不同类别的J-M距离；具体如公式(1-4)(1-5)所示：

JM_ij＝2(1-e^-B) (1-4)

其中，μ_i，μ_j分别为第i，j类的均值矢量，Σ_i，Σ_j分别为第i，j类的协方差矩阵。由于本发明仅有未种植耕地和种植耕地两个类别，因此i类和j类分别代表未种植耕地和种植耕地两个类别。

而后比较多个所述J-M距离的大小，得到J-M距离最大的两个特定波段组合，该波段组合即为区分未种植耕地与有作物种植耕地的最佳波段组合。所述最佳波段组合对应的两个时期即为未种植耕地提取的关键时期，选择该时期对应的环境星数据，即修正用的植被指数数据，来进行未种植耕地的识别。

在得到了识别用关键时期之后，在步骤S32，基于所述修正用遥感数据，在所述识别用关键时期，识别各像元对应的耕地为未种植耕地或种植耕地。

在一个具体实施例中，选择经过预处理后的识别用关键时期覆盖的某个区域(如黄淮海平原)的环境星影像。首先利用2010年黄淮海平原的土地利用数据对环境星CCD影像中的非耕地区域进行去除，并计算出两个识别用关键时期环境星CCD影像对应的植被指数。可以预见的是，在这两个关键时期中，其中一个关键时期一定处于作物生长高峰时期，因为作物生长高峰阶段未种植耕地的植被指数将显著低于植被区，因而，在这里将植被指数值大的时期设为高峰时期T_高峰，则另一时期设为T₂。

针对环境星CCD遥感数据中的耕地像元，逐像元来识别耕地像元的种植状况，即有作物种植，还是没有作物种植。具体如图3所示：

步骤S321，获取T_高峰时期的像元i的植被指数值NDVI_高峰。

步骤S322，比较所述像元i的植被指数值NDVI_高峰与预设阈值X的大小，判断像元i的植被指数值NDVI_高峰是否大于或等于阈值X，如果大于或等于，则确认当前像元i为种植耕地，并执行步骤S326；如果像元i的植被指数值NDVI_高峰小于阈值X，转到步骤S323。

步骤S323，获取T₂时期的像元i的植被指数值NDVI₂。

步骤S324，计算两个时期对应的植被指数差值NDVI_差值，其中，NDVI_差值＝NDVI_高峰-NDVI₂。

步骤S325，比较NDVI_差值与预设阈值Y的大小，判断NDVI_差值是否大于等于阈值Y,如果NDVI_差值大于等于阈值Y，则确认当前像元i为种植耕地，并执行步骤S326；如果NDVI_差值小于阈值Y，则确认当前像元i为未种植耕地，并执行步骤S326。

步骤S326，判断是否识别完所有的像元，如果识别完了所有的像元，则识别完成，结束该流程；如果没有识别完所有的像元，在步骤S327，取下一个像元，并返回步骤S321，重复上述步骤。

其中，以地面采样过程中获取的种植作物耕地区训练样本，统计T_高峰时期的植被指数，绘制种植作物耕地区高峰植被指数的直方图，截取频率为2％处的植被指数值作为种植区提取阈值，记为X。

以地面采样过程获取的未种植耕地区训练样本，计算出T_高峰时期和T₂时期的植被指数差值，并统计绘制出未种植耕地区训练样本区域的植被指数差值直方图，截取累积频率为98％处的植被指数差值作为未种植耕地识别的阈值之一，记为Y。

通过上述步骤，识别完了环境星CCD遥感数据中耕地种植状况。

在本实施例中，由于识别用的环境星CCD遥感数据与监测用的遥感数据(MODISNDVI)的分辨率不同，像元大小不同，为了获得环境星CCD遥感数据与MODIS NDVI像元大小一致的未种植耕地比例(UALR)数据，本发明利用最邻近方法，在步骤S33将30米分辨率的环境星CCD遥感数据中的种植耕地与未种植耕地分布图重新采样到250米分辨率上，根据采样后各像元的覆盖区域中采样前未种植耕地像元个数(Nuc)和有作物种植耕地像元个数(Nc)，根据公式(2-1)推导计算出采样后的各像元中未种植耕地比例UALR，具体公式如下：

本发明假定在冬小麦生长季内各像元未种植耕地比例UALR数值基本保持不变。

在步骤S4，本实施例将MODIS影像各耕地像元考虑为仅由有作物种植和无作物种植的两种端元组成的混合像元，其他端元忽略不计，并假定混合像元的地表反射率由两种端元的地表反射率按照各自丰度构成的线性混合。由于需要计算MODIS像元尺度的NDVI数据，因此本实施例中仅对红波段和近红外波段的反射率进行线性混合像元分解。各混合像元内，“有作物种植”端元的红波段、近红外波段的地表反射率可以通过求解下列的二元线性方程组获得：

其中ρ_R和ρ_NIR为某一MODIS像元尺度的端元红波段与近红外波段反射率，UALR代表对应像元的未种植耕地比例，ρ_Rcropped、ρ_NIRcropped和ρ_Runcropped、ρ_NIRuncropped分别代表在种植作物的耕地端元红波段反射率、近红外波段反射率以及未种植耕地端元红波段反射率、近红外波段反射率，其中未种植耕地端元红波段的反射率和近红外波段的反射率通过统计未种植耕地比例大于98％的像元对应波段的均值确定，因此式(2-2)中仅有有作物种植耕地端元红波段、近红外波段反射率(ρ_Rcropped、ρ_NIRcropped)两个未知数，通过求解式(2-2)中的二元线性方程组即可获得ρ_Rcropped、ρ_NIRcropped的值。

由此，根据所述ρ_Rcropped、ρ_NIRcropped的值，根据公式(2-3)，以监测用数据，即16天合成的时间序列MODIS数据为数据源，计算出各MODIS像元内有作物种植耕地部分的NDVI值，计算公式见式(2-3)。

基于前述步骤S4中修正后得到的NDVI数据，进行作物长势监测，包括作物长势实时监测与作物长势过程监测两种方法。

对于作物长势实时监测，其反映的是当前监测时段与前一年同时段的NDVI变化特征，识别出哪些区域在当前时期作物生长状况优于去年，哪些区域不及去年，哪些区域与去年作物长势相当。实时监测突出的是两年间空间上的差异性。而对作物长势过程监测，记录的是某一区域生长季内NDVI的时间序列变化，通常利用各像元耕地比例作为权重统计某一区域各期NDVI的均值，并与往年NDVI动态变化曲线进行对比分析，作物长势过程监测突出的是作物生长状况的动态变化并能反映出不同年份之间不同生长阶段的作物长势差异。

实施例二

如图4所示，为本发明提供的一种农作物长势监测数据处理模块，包括：数据获取子模块1、数据预处理子模块2、修正参数计算子模块3和数据修正子模块4。另外，为了完整地体现出本发明的技术方案，在图4中还包括数据存储子模块5，如各种存储介质。具体地，所述数据获取子模块1获取的数据包括遥感数据，即监测用遥感数据和修正用遥感数据。还包括各种参数，如阈值、公式中涉及到必要的参数等，这类参数数据，在获取到后，直接存储到数据存储子模块5中，并根据相应模块或单元的读取指令，向其提供相应的参数。获取到遥感数据后，将其发送给数据预处理子模块2。其中，监测用遥感数据为中低分辨率的遥感数据，例如通过访问相关网站下载或通过数据接口从外部存储装置中读取数据等方式，获得本发明所需要的数据。如前述实施例一中所述的监测用中低分辨率的MODISNDVI数据、修正用的环境1A、1B星的CCD数据及用于获取相应阈值及验证的地面观测得到的样本数据等。另外，数据获取子模块1还包括获取一些必要的计算参数、公式等。如未种植耕地端元红波段、近红外波段的反射率ρ_Runcropped、ρ_NIRuncropped和前面实施例一中涉及到的各种公式。

所述数据预处理子模块2收到遥感数据后，分别对所述监测用遥感数据和修正用遥感数据进行预处理，例如，对所述MOD13Q1和MYD13Q1的三级(L3)产品的数据进行括数据拼接、重投影、数据类型转换和数据格式转换等。对所述的环境星CCD数据进行几何精校正和辐射校正等操作。由于这些操作在实施例一中已进行较为详细地说明，因而在此不再重复说明。经过上述的预处理后，得到对应的植被指数数据。

所述修正参数计算子模块3具体包括关键时期选择单元31、耕地种植状况识别单元32、采样单元33和修正参数计算单元34，数据获取子模块1为上述各单元提供必要的参数数据。其中，所述关键时期选择单元31用于选择关键时期。所述的关键时期为用于识别耕地种植状况的关键时期。在本发明中，将时期选择的问题转换为最佳波段组合问题，将两个最佳波段组合对应的时期确定为关键时期。具体如实施例一所述，在此不再重复说明。

耕地种植状况识别单元32基于所述修正用植被指数数据，在所述关键时期，识别各像元对应的耕地为未种植耕地或种植耕地。例如，耕地种植状况识别单元32按照图3所示的流程图逐像元地识别出耕地类别。在具体识别过程中，利用了预先设置的两个阈值，即种植区提取阈值X和未种植耕地识别阈值Y。这两个数据可以预先计算好后输入到本模块中，或根据本模块根据存储的地面采样训练样本数据计算获得。

所述采样单元33将修正用遥感数据按照监测用遥感数据的图像分辨率重新采样，例如，将原来分辨率为30m的遥感数据按照250m重新采样，从而得到新的像元。对应于采样前的像元，重新采样后一个像元则覆盖了多个采样前的像元，并且，这些像元的类别已确定，因而可以得到采样后的一个像元中未种植耕地的像元数量和种植耕地的像元数量。

所述修正参数计算单元34根据公式(2-1)、重新采样后各像元所覆盖的采样前未种植耕地和种植耕地的像元数量，计算得到未种植耕地的比例，所述未种植耕地的比例为所述的修正参数。

所述数据修正子模块4包括两个计算单元：第一计算单元41和第一计算单元42。其中，所述第一计算单元41根据公式(2-2)的二元线性方程，及未种植作物的耕地所处像元的红波段、近红外波段的反射率ρ_Runcropped和ρ_Runcropped，计算得到种植作物的耕地端元的红波段、近红外波段反射率ρ_Rcropped、ρ_NIRcropped；

其中，所述ρ_R和ρ_NIR为MODIS像元尺度的端元红波段与近红外波段反射率；

所述UALR为对应像元的未种植耕地比例；

ρ_Rcropped为种植作物耕地端元的红波段的反射率；

ρ_NIRcropped为种植作物耕地端元的近红波段的反射率；

ρ_Runcropped为未种植作物耕地端元的红波段的反射率；

ρ_NIRuncropped为未种植作物耕地端元的近红外波段的反射率。

所述第二计算单元42根据公式(2-3)，针对监测用植被指数数据，逐像元计算出对应的植被指数值，从而获得修正后用于农作物长势监测的植被指数数据；

通过上述农作物长势监测数据处理模块，实现了对中低分辨率MODIS数据的修正。

实施例三

本实施例三提供了一种农作物长势监测系统，其中，包括前面实施例二中提供的数据处理模块1a和监测模块2a，通过数据处理模块1a得到经过修正的中低分辨率MODIS数据。所述监测模块2a以所述修正后的植被指数数据作为监测数据，进行作物长势的过程监测和/或作物长势的实时监测。

为了便于查看，所述监测系统还包括数据输出模块3a，所述数据输出模块3a为打印模块或/和显示模块等。用于将监测数据以特定的形式输出，如以图表、图形的方式显示或打印。

应用实施例：

采用本实施例中的数据处理方法和模块，对取得的遥感数据进行修正，以修正前和修正后的数据作为数据源进行作物长势监测，对监测结果进行对比分析如下：

首先利用分类获取的黄淮海地区30米分辨率的作物种植区与未种植区耕地分布数据，按前述的未种植耕地比例的统计方法，重采样获得了250米分辨率的黄淮海平原冬小麦生长季内的未种植耕地比例分布图。显而易见的，未种植耕地比例分布图与中的在耕地与未种植耕地分布具有非常相似的分布特征。2011年渤海湾西海岸以及黄河三角洲附近区域的未种植耕地比例明显高于2010年；同样的，2011年河南省西部地区的未种植耕地比例较2010年也有所增加，尽管这些区域的未种植耕地比例与渤海湾西海岸附近区域相比总体较低。

而后，关于作物长势实时监测，利用2010年与2011年两年的NDVI数据进行对比分析，根据获得的研究区的作物长势实时图，经过对比，两幅数据源不同的作物长势图具有相似的空间分布态势，但在个别区域还是存在一定差异的。统计结果显示，利用修正前MODISNDVI数据生成的作物长势图中，研究区约有8％的耕地作物长势好于2010年，约9.2％的耕地作物长势比2010年5月上旬差，而其余大部分地区的作物长势与2010年持平。而利用修正后的NDVI数据导出的作物长势图中，约有13.4％的有作物种植的耕地长势好于2010年同期，而14.3％的在农作物长势与2010年相比偏差，且某些区域的作物长势明显偏差，而原MODIS NDVI数据生成的作物长势图中，该地区的作物长势仅比2010年同期稍差。这一点与未种植耕地比率分布图的分布态势相似，未种植耕地比例监测结果显示黄淮海平原南部区域(主要包括河南最南部和安徽北部)在2009-2010年度冬小麦生长季内保留了较2010-2011年度更多的未种植耕地，因此直接利用MODIS NDVI数据开展作物长势监测时有一部分作物长势较差的信号被未种植耕地的变化所掩盖，而利用修正后的NDVI数据开展作物长势监测时，恰恰恢复了真实的作物长势状况。

在以本发明提供的监测系统对河南省最北部和河北省最南部地区进行监测时，如果以修正前的MODIS NDVI数据监测作物长势，其结果显示2010年和2011年作物长势基本持平，但应用修正之后NDVI数据时，其监测结果却显示该区域作物长势较2010年稍差，这些长势好坏统计比例以及空间分布态势的变化进一步凸显了利用未种植耕地比例进行NDVI修正来完善作物长势监测的重要性。

关于作物长势过程监测的对比，以对黄淮海平原以及其覆盖的五个省级行政单元开展长势监测得到结果为例，2009-2010年度与2010-2011年度的植被指数变化范围差异较大，2009-2010年度冬小麦生长季内，修正前MODIS NDVI数据的数值范围为0.19-0.66，而2010-2011年度的修正前MODIS NDVI数据的数值范围为0.31-0.63；相应时间段的修正后NDVI数据的数值范围分别为0.2-0.68及0.32-0.63。从上述数据可以得知，经过修正后的NDVI值与修正前MODIS NDVI值均有较大差异。

因而，与现有技术相比，本发明具有以下明显的技术优势：

(1)耕地种植状况归一化修正后的NDVI数据因其结合了年度间耕地利用状况(有作物种植耕地与未种植耕地)的空间变化，解决了年度间耕地利用状况变化以及作物轮作对作物长势的影响。更适合用于作物长势监测，由此生成的作物长势监测结果较以原始NDVI数据为数据源的作物长势监测结果更加准确。

(2)本发明提出的农作物长势监测系统仅需要区分出有作物种植的耕地与没有作物种植的耕地两类，不需要精确识别出不同作物类别，便可对原始NDVI数据进行修正，从而开展更加准确的作物长势监测，这种对种植耕地与未种植耕地的识别更加适用于遥感数据的处理。

(3)本发明采用混合像元分解策略，有效规避了植被指数的非线性混合效应，在一定程度上减小了修正后的纯净种植区NDVI数据的误差。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (10)

1.一种农作物长势监测数据处理方法，其中，包括：

获取监测用遥感数据和修正用遥感数据；

分别对所述监测用遥感数据和修正用遥感数据进行预处理，得到对应的植被指数数据；

根据修正用植被指数数据获得修正参数；

根据所述修正参数修正所述监测用植被指数数据。

2.如权利要求1所述的农作物长势监测数据处理方法，其中，所述根据修正用植被指数数据获得修正参数的步骤包括：

选择识别用关键时期；

基于所述修正用植被指数数据，在所述识别用关键时期，识别各像元对应的耕地为未种植耕地或种植耕地；

对所述修正用遥感数据按照监测用遥感数据的图像分辨率重新采样；

获取重新采样后各像元所覆盖的采样前未种植耕地和种植耕地的像元数量，根据公式(2-1)计算得到重新采样后各像元中未种植耕地的比例，所述未种植耕地的比例为所述的修正参数；

$UALR=\frac{N_{uc}}{N_{uc}+N_c}---(2-1).$

3.如权利要求2所述的农作物长势监测数据处理方法，其中，所述选择识别用关键时期的步骤包括：

采用类间Jeffries-Matusita距离对所述修正用植被指数数据进行量化；

对应于多个不同的两个时期的植被指数组合，分别计算对应的未种植耕地和种植耕地两种不同类别的Jeffries-Matusita距离；

比较多个所述Jeffries-Matusita距离的大小，获取Jeffries-Matusita距离最大的植被指数组合；

分别获取Jeffries-Matusita距离最大的植被指数组合对应的两个时期，将其作为用于识别的关键时期。

4.如权利要求3所述的农作物长势监测数据处理方法，其中，基于所述修正用植被指数数据，在所述识别用关键时期，识别各像元对应的耕地为未种植耕地或种植耕地，包括：

读取与两个所述关键时期对应的植被指数值；

计算所述两个植被指数值的差值，并将植被指数值大的关键时期设为高峰时期；

逐像元比较高峰时期的植被指数值与预设的种植区提取阈值的大小，如果高峰时期的植被指数值大于或等于所述种植区提取阈值，则将所述像元对应的耕地识别为种植耕地；

如果高峰时期的植被指数值小于所述种植区提取阈值，逐像元比较所述两个所述关键时期对应的植被指数值的差值与预设的未种植区提取阈值的大小，如果所述差值大于或等于所述未种植区提取阈值，则将所述像元对应的耕地识别为种植耕地，如果所述差值小于所述未种植区提取阈值，则将所述像元对应的耕地识别为未种植耕地。

5.如权利要求2所述的农作物长势监测数据处理方法，其中，根据所述修正参数修正所述监测用植被指数数据，包括：

根据红波段和近红外波段的反射率在有种植作物端元和未种植作物端元的关系，获得以下二元线性方程：

$\left\{\begin{array}{c}{\mathrm{\ρ}}_R={\mathrm{\ρ}}_{R_{cropped}}*\left(1-UALR\right)+{\mathrm{\ρ}}_{R_{uncropped}}*UALR\\ {\mathrm{\ρ}}_{NIR}={\mathrm{\ρ}}_{{\mathrm{NIR}}_{cropped}}*\left(1-UALR\right)+{\mathrm{\ρ}}_{{\mathrm{NIR}}_{uncropped}}*UALR\end{array}\right.---\left(2-2\right)$

其中，所述ρ_R和ρ_NIR为MODIS像元尺度的端元红波段与近红外波段反射率；

所述UALR为对应像元的未种植耕地比例；

ρ_Rcropped为种植作物耕地端元的红波段反射率；

ρ_NIRcropped为种植作物耕地端元的近红波段反射率；

ρ_Runcropped为未种植作物耕地端元的红波段反射率；

ρ_NIRuncropped为未种植作物耕地端元的近红外波段反射率；

其中，未种植作物耕地端元的红波段的反射率ρ_Runcropped和近红外波段反射率ρ_Runcropped通过统计未种植耕地比例大于98％的像元对应波段的均值确定；

根据公式(2-2)计算获得修正后的种植作物耕地端元的红波段反射率ρ_Rcropped和近红外波段反射率ρ_NIRcropped；

根据公式(2-3)逐像元计算出对应的植被指数值，从而获得修正后用于农作物长势监测的植被指数数据；

${\mathrm{NDVI}}_{cropped}=\frac{{\mathrm{\ρ}}_{{\mathrm{NIR}}_{cropped}}-{\mathrm{\ρ}}_{R_{cropped}}}{{\mathrm{\ρ}}_{{\mathrm{NIR}}_{cropped}}+{\mathrm{\ρ}}_{R_{cropped}}}---(2-3).$

6.一种农作物长势监测数据处理模块，其中，包括：

数据获取子模块，用于获取监测用遥感数据和修正用遥感数据；

数据预处理子模块，用于分别对所述监测用遥感数据和修正用遥感数据进行预处理，得到对应的植被指数数据；

修正参数计算子模块，用于根据修正用遥感数据获得修正参数；和

数据修正子模块，用于根据所述修正参数修正所述监测用植被指数数据。

7.如权利要求6所述农作物长势监测数据处理模块，其中，所述修正参数计算子模块包括：

关键时期选择单元，用于选择识用别关键时期；

耕地种植状况识别单元，用于基于所述修正用植被指数数据，在所述识别用关键时期，识别各像元对应的耕地为未种植耕地或种植耕地；

采样单元，用于将修正用遥感数据按照监测用遥感数据的图像分辨率重新采样，并获取重新采样后各像元所覆盖的采样前未种植耕地和种植耕地的像元数量；

修正参数计算单元，用于根据重新采样后各像元所覆盖的采样前未种植耕地和种植耕地的像元数量，按照公式(2-1)计算得到未种植耕地的比例，所述未种植耕地的比例为所述的修正参数：

$UALR=\frac{N_{uc}}{N_{uc}+N_c}---(2-1).$

8.如权利要求6所述农作物长势监测数据处理模块，其中，所述数据修正子模块包括：

第一计算单元，用于根据公式(2-2)提供的二元线性方程，及未种植作物耕地端元的红波段反射率ρ_Runcropped和近红外波段反射率ρ_Runcropped计算得到种植作物耕地端元红波段反射率ρ_Rcropped和近红外波段反射率ρ_NIRcropped；

$\left\{\begin{array}{c}{\mathrm{\ρ}}_R={\mathrm{\ρ}}_{R_{cropped}}*\left(1-UALR\right)+{\mathrm{\ρ}}_{R_{uncropped}}*UALR\\ {\mathrm{\ρ}}_{NIR}={\mathrm{\ρ}}_{{\mathrm{NIR}}_{cropped}}*\left(1-UALR\right)+{\mathrm{\ρ}}_{{\mathrm{NIR}}_{uncropped}}*UALR\end{array}\right.---\left(2-2\right)$

其中，所述ρ_R和ρ_NIR为MODIS像元尺度的耕地端元的红波段与近红外波段反射率；

所述UALR为对应像元的未种植耕地比例；

和

第二计算单元，用于根据公式(2-3)，针对监测用植被指数数据，逐像元计算出对应的植被指数值，从而获得修正后用于农作物长势监测的植被指数数据；

${\mathrm{NDVI}}_{cropped}=\frac{{\mathrm{\ρ}}_{{\mathrm{NIR}}_{cropped}}-{\mathrm{\ρ}}_{R_{cropped}}}{{\mathrm{\ρ}}_{{\mathrm{NIR}}_{cropped}}+{\mathrm{\ρ}}_{R_{cropped}}}---(2-3).$

9.一种农作物长势监测系统，其中，包括如权利要求6-8任一所述的数据处理模块；和

监测模块，以所述数据处理模块修正后的植被指数数据作为监测数据，进行作物长势的过程监测和/或作物长势的实时监测。

10.如权利要求9所述的农作物长势监测系统，其中，还包括：数据输出模块，与所述监测模块相连接，用于输出监测结果数据。