CN104915674B

CN104915674B - Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法

Info

Publication number: CN104915674B
Application number: CN201410572829.9A
Authority: CN
Inventors: 张锦水; 潘耀忠; 谢登峰; 孙佩军; 袁周米琪
Original assignee: Beijing Normal University
Current assignee: Beijing Normal University
Priority date: 2014-10-24
Filing date: 2014-10-24
Publication date: 2018-12-14
Anticipated expiration: 2034-10-24
Also published as: CN104915674A

Abstract

本发明公开了一种Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法，主要包括数据准备、数据预处理、利用STDFA模型构建数据、数据滤波处理、秋粮识别特征向量选择、秋粮识别特征组合和SVM分类获得识别结果。本发明可以有效区分秋粮作物的类型，在一定程度上可以解决秋粮识别过程中高分辨率数据不足的问题。根据实验结果可知，对秋粮识别的各项精度均达到了较高的识别精度，也在一定程度上验证了物候数据有助于秋粮识别。

Description

Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法

技术领域

本发明涉及一种利用遥感数据时空融合方法(Spatial Temporal Data FusionApproach，STDFA)构建高时间、空间分辨率的遥感影像数据。并以此为基础，构建分类数据集，然后利用支持向量机(Support Vector Machine，SVM)进行秋粮作物识别的方法，验证不同维度的数据集进行秋粮作物识别的适用性。属于遥感影像应用领域。

背景技术

秋粮作物是我国粮食作物的重要组成部分，准确地获取秋粮作物种植面积对制定农业政策和管理具有十分重要的意义。利用遥感技术进行农业监测是其重要的一个应用方向，无论在面积监测还是遥感估产等方面发挥着重要作用。然而，秋粮作物分布广泛，种植结构复杂，同期作物生长，造成“异物同谱”现象，这对用于作物识别的遥感数据要求比较高；秋粮作物生长关键期内云雨天气较多，增加了有效遥感数据获取的难度，对利用遥感技术进行秋粮作物的识别带来了一定的挑战。

搭载在Terra和Aqua卫星上的中分辨率成像光谱仪(Moderate-resolutionImaging Spectroradiometer，MODIS)是一种广泛应用于农作物监测和识别的传感器。MODIS传感器具有很高的时间分辨率，重返周期为1-2天，能很好地记录农作物生长的时间信息，应用于地表植被的生长监测和识别，这对大尺度、种植结构单一的作物(如小麦、水稻)具有很大的优势。利用MODIS 250米分辨率的NDVI时间序列产品数据对大面积种植的水稻进行识别，能够获得较高的识别精度。但我国秋粮作物、种植结构复杂(玉米、水稻、棉花、大豆和花生等交错种植)、种植地块破碎，再加上自然植被(树木等)的影响，这使得MODIS数据的空间分辨率限制了破碎地块小面积种植的秋粮的识别。中高分辨率的卫星具有比较适中的空间分辨率，如美国的陆地卫星(Landsat)和法国SPOT卫星，在作物识别中发挥着重要的作用。Landsat卫星空间分辨率为30米，回返周期为16天，对于大面积种植的作物具有较高的识别精度。但该卫星的时间分辨率会限制植被生长关键期遥感数据的获取，再加上我国的雨季与秋粮作物的生长期重叠，这使得获取作物生长关键期高质量的遥感影像数据(云量＜10％)非常困难(获取概率＜10％)。

发明内容

本发明的目的是针对现有技术中存在的问题，以遥感数据时空融合方法(SpatialTemporal Data Fusion Approach，STDFA)构建高时间、空间分辨率的遥感影像数据。以此为基础，构建15种30米分辨率分类数据集，然后利用支持向量机(Support VectorMachine，SVM)进行秋粮作物识别，验证不同维度的数据集进行秋粮作物识别的适用性。

本发明通过以下技术方案解决上述问题，本发明提供一种Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法，主要包括以下步骤：

步骤一、数据准备；

步骤二、数据预处理；

步骤三、利用STDFA模型构建数据；

步骤四、数据滤波处理；

步骤五、秋粮识别特征向量选择；

步骤六、秋粮识别特征组合；

步骤七、SVM分类获得识别结果。

优选的，上述步骤一～步骤七具体为：

步骤一、数据准备，选取目标区域，获取同一区域的Landsat8影像数据、 MODIS数据和航拍数据；

步骤二、数据预处理，包括对步骤一获取的Landsat8影像数据进行几何校正、对MODIS数据进行格式转换和重采样，对航拍数据进行坐标系转换、解译、矢量化和属性赋值并形成栅格数据；

步骤三、利用STDFA模型构建数据，即以步骤二的数据为基础，利用数据时空融合STDFA模型构建红波段(RED)和近红外波段(NIR)高时空分辨率数据，再利用RED和NIR数据计算归一化植被指数(NDVI)数据；

步骤四、数据滤波处理，对上述步骤三构建的RED、NIR和NDVI数据进行滤波处理，并对NDVI时间序列数据滤波提取的多种物候指标组成物候数据(Phenology)；

步骤五、秋粮识别特征向量选择，根据航拍数据，选取时间序列中能区分地物的波段用于后续的秋粮识别；

步骤六、秋粮识别特征组合，以步骤四滤波后的Red、NIR和NDVI时间序列数据以及重组的物候指标数据为基础，进行所有可能的类型组合用于秋粮识别；

步骤七、SVM分类获得识别结果，即利用对高维数据分类具有明显优势的支持向量机(SVM)的分类方法分别步骤六的数据组合类型进行分类，得到秋粮作物种类的识别结果。

优选的，上述方法还包括步骤八、采用航拍数据评价分类结果，对选取目标区域拍影像矢量化数据转化的栅格数据进行重分类，使重分类后的分类结果与所述步骤七的分类数据进行像元对像元的叠加比较，得到用于分类结果精度评价的混淆矩阵。

优选的，上述步骤二数据预处理中对MODIS地表反射率产品首先转化为反射率数据，然后利用MODIS重投影转换成与Landsat8影像一致的坐标系，将MODIS产品的数据进行格式转换，并重采样到Landsat8数据分辨率的整数倍的空间分辨率，以便进行后续的MODIS混合像元分解。

优选的，上述步骤三种以MOD09Q1和Landsat8为基础数据利用数据时空融合STDFA模型构建三种高时空分辨率数据(Red、NIR和NDVI)，所述NDVI数据是利用时空数据融合STDFA模型构建的Red和NIR数据，通过下列公式计算得到：

式中，ρ_NIR与ρ_Red分别表示近红外(NIR)波段和红波段(Red)的反射率。

优选的，上述步骤四是以Savitzky-Golay(S-G)拟合法为核心算法对时间序列数据进行重构建。

优选的，上述的航拍数据为高分辨率的无人机航拍数据，拍摄多个无人机航拍样方，对多幅样方进行拼接后转换坐标系，并且目视解译矢量化及地块属性赋值，以面积占优的类型转换成一定分辨率的栅格数据。

本发明提供的Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法可以有效区分秋粮作物的类型，在一定程度上可以解决秋粮识别过程中高分辨率数据不足的问题。根据实验结果可知，对秋粮识别的各项精度均达到了较高的识别精度，也在一定程度上验证了物候数据有助于秋粮识别。

附图说明

图1为本发明实施例研究区位置及其四期Landsat8影像图；

图2为本发明实施例无人机样方及分布图；

图3(a)-3(1)为本发明实施例融合影像与真实Landsat8影像比较及对应波段的相关系数图；(a)为融合的2013-7-20Red影像图；(b)为真实的2013-7-26Landsat8Red波段影像图；(c)为融合的2013-07-20Red波段影像与真实2013-07-26Red影像散点图及相关系数；(d)为融合的2013-7-20NIR影像图；(e)为真实的2013-7-26Landsat8NIR波段影像图；(f)为融合的2013-07-20NIR波段影像与真实2013-07-26NIR影像散点图及相关系数；(g)为融合的2013-08-05Red影像图；(h)为真实的2013-08-11Landsat8Red波段影像图；(i)为融合的2013-08-05Red波段影像与真实的2013-08-11Red影像散点图及相关系数；(j)为融合的2013-08-05NIR影像图；(k)为真实的2013-08-11Landsat NIR波段影像图；(1)为融合的2013-08-05NIR波段影像与真实的2013-08-11NIR影像散点图及相关系数；

图4为本发明方法流程示意图；

图5为本发明实施例不同数据典型地物曲线图；(a)、(b)、(c)表示三种构建数据的不同地物时间序列曲线；(d)图表示提取的不同地物物候特征，横坐标表3中的物候指标，纵坐标表示物候指标归一化值；

图6为本发明实施例矢量数据与分类结果子区比较图。

图中标记：1-森林；2-水稻；3-玉米；4-芦苇；5-居民区。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及具体实施方式对本发明作进一步的详细描述。

如图1所示，本发明公开的Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法，结合MODIS数据和中高分辨率数据各自的优点进行秋粮作物识别。

融合MODIS数据与Landsat数据构建高时空遥感数据是综合二者优势的一个有效方法。近些年，国内外学者提出了几种高时间分辨率数据和高空间分辨率数据的融合方法，构建高时空分辨率遥感数据应用于不同方面的研究工作。这些方法大部分都是基于线性模型分解粗分辨率混合像元。在考虑像元反射率的环境影响时，不同学者改进了线性分解模型，提出相应的融合模型。Maselli(2001)提出了一种基于像元反射率在一定邻域范围内不会发生剧烈的变化的线性分解模型集成高低时间分辨率的NDVI数据用于监测植被的研究；Gao等人(2006)不仅考虑到像元间的距离和光谱差异，还考虑了像元时间上的差异而提出了一种自适应遥感图像融合模型(spatial and temporal adaptive reflectance fusionmodel，STARFM)用于破碎地块地表覆盖类型的识别；Hilker等人(2009)提出了一种提取反射率变化的时空自适应融合算法(spatial temporal adaptive algorithm for mappingreflectance change，STAARCH)用于监测森林覆盖的研究；Walker等人(2012)利用MODIS与Landsat的融合数据用于干旱区森林物候的分析研究；Wu等人(2012)提出了一种基于时间变化特征的时空数据融合模型(spatial temporal data fusion approach，STDFA)用于水稻面积的提取等。

融合MODIS数据与Landsat数据构造高时空分辨率遥感数据，生成高维度的遥感数据，选择适当的高维度遥感数据进行秋粮作物识别，是保证秋粮识别精度关键。本发明基于STDFA时空融合模型，构造高时空间分辨率的数据(红波段数据、近红外波段数据、NDVI数据)，利用Timesat软件对NDVI时间序列数据进行滤波提取物候特征数据(Phenology)。针对这4种类型的数据，分析数据对典型地物(玉米、水稻)的可分性并确定分类数据的组合类型，用于探讨秋粮识别的可行性及最佳数据组合类型，为利用高时空遥感数据在大范围进行秋粮作物的遥感识别积累一定的实验基础。

为了验证本发明提出的方法有效性，本发明选择了实地区域进行了实验研究，研究区位于辽宁省锦州市、盘锦市和鞍山市境内，其范围为：纬度40°58′47"N-41°46′51″N，经度121°26′0"E-122°27′02"E(如图1)。该区域地势西北高、东南低，从海拔高的山区，向东南逐渐降到海拔20米以下的海滨平原。研究区位于中纬度地带，属于暖温带大陆性半湿润季风气候，常年温差较大，全年平均气温8℃～9℃，年降水量平均为540～640毫米，无霜期达180天。该区域的自然条件决定了农作物为一年一熟，秋粮作物以玉米和水稻为主，兼有少量的花生等。西北部是山区生长着森林，其它地区是平原，主要种植着水稻和玉米，并兼分布着零散的居民地，其中南部地区是大面积居民地盘锦市区，同时在河流附近生长着大面积的水生植物芦苇(图1(a))。

一、数据与预处理

1、Landsat 8数据

选取研究区秋粮作物生长物候期5月至10月份内的4景Landsat 8影像数据(如图1)，其投影坐标系为UTM-WGS84Zone 51N，获取时间分别为2013-05-23，2013-07-26，2013-08-11和2013-09-12。4景Landsat8陆地成像仪(OLI)各波段数据在ENVI5.0sp3软件中利用辐射校正功能转换为反射率数据，选取实验用的红、绿、蓝、近红和两个短波红外共6个波段。然后，对Landsat8数据的红、绿和蓝波段采用暗目标减法(Dark Object Subtraction，DOS)进行大气校正，而近红外波段因其波长比较大，可以忽略大气对该波段的影响。由于Landsat 8数据做过基于地形数据的几何校正，一般情况下可以直接使用而不需要做几何校正。因此，本实验中Landsat8数据不再进行几何校正。

2、MODIS数据

由于秋粮作物生长季节云雨天气偏多，MODIS每天地表反射率产品数据受到“云污染”的影响较大，导致数据无法使用。因此，本实验采用MODIS 8天合成的250m分辨率的地表反射率产品数据(MOD09Q1)，该数据只包含的红波段(Red)和近红外波段(NIR)，与其对应的Landsat8(OLI)波段如表1所示。MOD09Q1产品数据的标注日期是合成日期中的第一天。因此，选择与秋粮作物生长周期相对应的2013-05-17至2013-09-22期间的MODIS数据。

表1 用于实验的Landsat 8和MODIS影像数据波段信息

MODIS地表反射率产品首先乘以0.0001转化为-1～1之间的反射率数据，然后利用MODIS重投影工具(MODIS Re-projection Tool，MRT)转换成与Landsat8影像一致的UTM-WGS84坐标系，同时把MODIS产品的数据格式转换成Geo-tif格式，并采用最近邻域法重采样到240m空间分辨率，为Landsat8数据30m分辨率的整数倍，以便进行后续的MODIS混合像元分解。用于本实验中Landsat8和MODIS数据主要特征如表2所示。

表2 用于实验的Landsat 8和MODIS影像数据的主要特征

3、验证数据

本实验采用高分辨率的无人机航拍数据用于评价秋粮作物的识别精度，研究区内共有13个无人机航拍样方(如图2)，大小为1500m×1800m，约3KM2。在ArcGIS软件中把拼接后的无人机航拍数据转换成UTM-WGS84坐标系，然后目视解译矢量化及地块属性赋值，以面积占优的类型转换成30m分辨率的栅格数据，用于后续秋粮分类结果的精度评价。

二、基于STDFA模型构建高时空间分辨率数据

本发明根据时空数据融合方法(spatial temporal data fusion approach，STDFA)和MOD09Q1数据产品的特点构建红波段和近红外波段数据。由于Wu等利用STDFA模型开展实验研究时，选择的Landsat影像是无云的高质量遥感数据，而现实应用中，尤其是秋粮作物识别的遥感影像上云是经常出现的。本实验中的四期Landsat8影像(如图1(a)-(d))除2013-05-23影像没有受到“云污染”外，其余三期影像均有少量的云覆盖。因无法消除云对STDFA模型中地表变化聚类的影响，在此提出假设：云在各期遥感影像的分布位置是不尽相同的，这在一定程度上可以减少云对地表变化聚类的影响。具体实现方法是：用2013-05-23遥感影像分别与2013-07-26、2013-08-11、2013-09-12三期影像进行ISO-Data聚类，提取地表覆盖的变化类型，获得三幅具有相同变化类型数量的地表变化聚类图C1、C2和C3；以C3为基准，把C1、C2中非云和阴影变化类像元的属性与C3中相应像元的属性相统一；然后，把C3中云和阴影类的像元值用C1或C2中相应的非云和阴影类像元值进行替换。运算后，得到一幅最终的地表变化类聚类图C，用于STDFA模型中的MODIS混合像元分解。

统计MODIS像元内的地表变化类及每种变化类面积占该MODIS像元面积的比例，即得到地表变化类c的丰度f_c(i，c)。利用全约束的混合像元线性分解模型分别对多时相的MODIS地表反射率产品中的红波段(Red)数据和近红外波段(NIR)数据按公式(1)进行分解，得到红波段和近红波段t_j时期的地表变化类c的平均反射率

约束条件：且f_c(i，c)≥0

式中：R(i，t_j)为t_j时期i位置的MODIS混合像元的反射率；f_c(i，c)为i位置的MODIS混合像元内地表变化类c像元占该混合像元的面积比；为t_j时期地表变化类c的平均反射率；ξ(i，t_j)为残差；k为研究区内地表变化类型的数量。

假设在时间段t₀至t_n及该时间段两端有限外延的时间范围内，地表同一变化类c像元的反射率变化趋势是一致的。在这个假设下，利用公式(1)和公式(2)可以构建t_j时期的高分辨率的数据r(c，t_j)，从而构建出高分辨率的时间序列数据。

式中：和分别是利用公式(1)求出的t_j时期和t₀时期的地表变化类c像元的平均反射率；r(c，t₀)是初期t₀时期的Landsat8遥感影像相对应c类像元的地表反射率。

因此，由8天合成的地表反射率产品数据(MOD09Q1)构建2013-07-20和2013-08-05两期Red和NIR波段的高分辨率影像，选取研究区2013-07-26和2013-08-11两期Landsat8影像质量较好的子区域与相应区域的融合影像进行对比分析(图3)。

从图3融合结果与真实Landsat影像的比较及对应波段的相关系数可以看出，融合影像与真实影像的目视效果比较接近并且两者也达到了较高的相关系数，其散点图分布在1∶1对角线的周围，融合影像在一定程度上能够反映同时期Landsat影像的光谱信息，可以用于秋粮作物的识别，但融合影像的质量及其与真实影像之间的相关系数还是受到其它一些因素的影响：

(1)融合的08-05影像与08-11真实影像对应红波段和近红外波段的相关系数分别达到0.86和0.81，而融合的07-20影像与07-26真实影像对应红波段和近红外波段的相关系数也分别达到0.78和0.73。单从相关系数上看，前者的融合效果比后者好，在图3(c)、(f)散点图上出现许多真实地表反射率与融合结果相差很大的异常像元，这很可能是由2013-07-26的Landsat8遥感影像上少量云和阴影(圆圈所示)造成的，也可能是造成后者相关系数较低的一个因素。

(2)MOD09Q1是8天合成的地表反射率数据，融合构建的数据也会受到合成期内其它日期地表反射率的影响。而与之进行比较的真实Landsat8影像只是合成期间内某一天的数据，两者存在着时相上的差异，这对融合影像与真实影像之间的相关性存在一定的影响。但融合的近红外影像反射率要比真实Landsat影像的近红外影像反射率高(图3(f)、(1))，而红波段影像则相反，造成这一现象可能是由于8月份左右是植被快速生长阶段，近红外波段反射率高，红波段反射率低，变化快，融合影像和真实影像之间是存在时相差异等造成的。

(3)在构建融合数据时，使用初期MODIS影像和初期Landsat影像存在时相的差异，这也会造成融合影像与同期真实影像之间存在一定的偏差。

三、秋粮作物识别

本实验是以MOD09Q1和Landsat8为基础数据利用数据时空融合STDFA 模型构建三种高时空分辨率数据(Red、NIR和NDVI)，在对NDVI时间序列数据进行滤波时提取作物物候特征数据，然后分析数据对秋粮作物的可分性，选择有效的识别特征向量组成不同分类数据集进行秋粮作物识别，具体流程如图4。

1、植被指数反演及滤波

植被指数是反映植被在可见光、近红外波段反射与土壤背景之间差异的指标，在一定条件下能用来定量说明植被的生长状况。而归一化植被指数(NormalizedDifferential Vegetation Index，NDVI)是目前应用广泛的植被指数，不同的NDVI值对应不同的土地覆被类型，可以进行土地覆被方面的研究。为此，利用时空数据融合STDFA模型构建的Red和NIR数据，通过公式(3)计算NDVI数据。

式(3)中，ρ_NIR与ρ_Red分别表示近红外(NIR)波段和红波段(Red)的反射率。

由于构建的Red、NIR和NDVI数据具有较多的噪音，为减少噪音对植被信息提取的干扰，需要对其进行滤波处理，本实验利用TIMESAT软件对构建的时间序列数据进行滤波，以Savitzky-Golay(S-G)拟合法为核心算法对时间序列数据进行重构建。

对构建的17期Red、NIR和NDVI时间序列数据分别复制两份组成三个周期的时间序列数据，利用与时间序列数据曲线形状比较接近的S-G滤波方法进行数据重构建，然后选取重构建数据的中间17期数据用于秋粮识别。典型地物在重构建数据中的时间序列曲线如图5(a)、(b)、(c)。

2、物候特征提取

作物的物候特征定量刻画了作物的生长过程，能够提高低分辨率遥感影像对土地覆盖类型的识别精度。利用多时相的MODIS数据与物候数据可以提高大尺度区域的荒废农用地的识别精度。为了验证物候数据是否也有助于秋粮识别，本实验对NDVI时间序列数据滤波提取的11种物候指标组成物候数据(Phenology)用于秋粮识别(如表3)。典型地物的物候指数曲线如图5(d)

表3 提取的物候指标

3、秋粮识别特征向量选择

由STDFA时空融合模型构建的时间序列数据可以有多种分类特征的组合，为了确保分类特征具有代表性，减少数据组合的冗余，需进行分类特征的选择。根据无人机航片数据，选取研究区内种植面积比较大的玉米和水稻以及居民区、水体、树木和芦苇等6种典型地物来分析4种时间序列数据对典型地物的可分性，选取时间序列中能区分地物的波段用于后续的秋粮识别，典型地物的时间序列曲线如图5。

在图5(a)、(b)、(c)中，由于树木和芦苇在五月中下旬时已长出大量的叶片，对红光的反射率比较低，对近红外反射率高，具有较高的NDVI值；而研究区内的水稻是在5月下旬插秧，随着水稻的生长，该区域的地表覆盖类型也经历了裸土-水体-植被的变化过程，表现在近红外和NDVI时间序列曲线上会有先下降后上升的趋势，在8月中上旬达到峰值；玉米一般在4月下旬和5月上旬种植，早于水稻，因此在时间序列中起始点的红波段的反射率是下降的，近红外波段和NDVI数据是上升的；水体里会因生长水草等使其在红、近红和NDVI时间序列曲线中会有一定的波动；对于居民区和道路，由于其周围树木的影响也会对其时间序列曲线有相应的波动。综上分析，构建的时间序列数据是符合典型地物的变化特征。而对于物候特征图曲线(图5(d))，绿色植被在f和g(生长期的中期时刻、NDVI峰值)物候特征上具有很高的相似性而不能有效地区分，因而选择除f和g特征之外的其它9个物候特征进行依次重组作为本实验的物候数据(Phenology)。

综上分析，滤波后的Red数据、NIR数据、NDVI及重组后的物候数据(Phenology)可以用于秋粮作物识别。

四、秋粮识别特征组合及SVM分类

以滤波后的Red、NIR和NDVI时间序列数据以及重组的物候指标数据等4种基本分类特征数据为基础，进行所有可能的类型组合用于秋粮识别，数据组合类型特征如表4：

结合Landsat8遥感影像目视解译、无人机航拍数据选取玉米、水稻、居民区、水体、树木和芦苇等6种典型地物训练样本。其他少量或不易划分类别的地物包含在具有相似光谱特征的典型地物之中：居民区包括居民地、裸地和道路；树木包括除芦苇之外的其它非作物绿色植被。利用对高维数据分类具有明显优势的支持向量机(Support VectorMachine，SVM)的分类方法分别对表4中的数据组合类型进行分类。

表4 用于作物分类的数据组合类型特征

数据集类型	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)	(14)	(15)
																波段数量	18	18	18	9	36	36	27	36	27	27	45	54	45	45	63

注：(1)Red(2)NIR(3)NDVI(4)Phenology(5)Red+NIR(6)Red+NDVI(7)Red+Phenology(8)NIR+NDVI(9)NIR+Phenology(10)NDVI+Phenology(11)Red+NDVI+Phenology(12)Red+NIR+NDVI(13)Red+NIR+Phenology(14)NIR+NDVI+Phenology(15)Red+NIR+NDVI+Phenology

五、结果与分析

1、精度评价

本实验采用无人机航拍数据作为分类结果精度评价的验证数据。对研究区内无人机航拍影像矢量化数据转化的栅格数据进行重分类，使重分类后的地物为水稻、玉米和其它；本实验15种分类数据的分类结果也分别进行相同的重分类；然后，使重分类后的分类结果与验证数据进行像元对像元的叠加比较，得到用于分类结果精度评价的混淆矩阵。混淆矩阵如表5，分类结果精度比较如表6。

表5 不同分类数据的分类结果混淆矩阵

表6 不同分类数据的分类结果精度比较

从表6分类结果精度比较中可以看出：

(1)每种组合类型的分类结果，秋粮识别的各项精度均达到了较高的识别精度，水稻的制图精度和用户精度均分别达90％和83％左右；玉米的制图精度和用户精度均分别达到85％和75％左右；二者的总体精度也达到85％左右。而每种组合类型对水稻的识别效果都好于玉米，其原因可能在于：

a)由于水稻和玉米的种植特点不同，水稻种植相对集中连片，而玉米种植周围往往会有其它小面积的作物。因在选择玉米训练样本时忽略其周围种植的小面积作物，这些作物会被分到光谱相似的典型地物(玉米)中，而无人机影像7-8cm的分辨率可以识别小面积作物，在采样为30米分辨率的验证数据时会有一定数量的像元表示这些小面积作物，精度验证时会造成一定的影响。

b)无人机航拍样方在水稻种植区域少且集中，在玉米种植区多且分布广泛，这也可能是作物分类精度差异的一个原因。

c)分类数据与验证数据存在地理配准误差，有像元错位现象。以像元对象元的方式进行精度评价时对种植不是集中连片的作物会有较大的影响。

(2)分类结果精度与组合分类数据集的数据类型数量之间并非正相关，即并不是分类数据组合的数据类型数量越多分类精度越高。本实验中，水稻和玉米的制图精度、用户精度以及两者的总体精度的最大值分别出现在Red+NDVI、NIR、Red、NDVI+Phenology和Red数据的分类结果中，分别达到了92.23％、85.44％、89.15％、76.21％和89.42％；各项精度的最小值都出现在单一数据组合的分类结果中，这说明仅使用单一的数据进行分类不是最佳的选择。但物候数据在一定程度上可以提高单一数据组合分类结果中的部分项的(Red、NDVI和NIR)。另外，物候数据对玉米识别的制图精度最低，这可能是玉米与周围其它地物具有相似的物候特征。Red+Phenology组合的分类结果中具有三项精度指标达到第二高值，其它两项指标也达到了较高的精度。综合分析，Red+Phenology组合对本次实验的适用性最好，可以作为本次试验的最佳分类组合，同时也在一定程度上说明了物候数据有助于秋粮识别。

(3)当分类数据集由3种或3种以上的数据组合时，水稻和玉米的制图精度、用户精度以及总体精度均趋向于稳定，不同分类数据的相同精度评价项之间的差异比较小，表明只有达到一定数据的数据组合时才能得到比较稳定的分类结果精度。

2、空间特征分析

在15种分类数据集的分类结果中，各分类精度项差别不是很大，为了突出不同数据集对玉米和水稻分类结果的影响，分别选取最大和最小分类精度项所在分类结果中的子区域与矢量化地块进行比较分析(图6)。

图6(a1)、(a2)是NIR数据和Red+NDVI数据的分类结果，Red+NDVI数据对水稻的分类效果好于NIR数据，前者水稻被漏分为玉米的像元少于后者，造成NIR水稻的制图精度低于Red+NDVI；同理，对图6(b1)、(b2)分类结果，Red数据对玉米的分类结果好于Phenology数据，相比前者，后者有部分玉米像元漏分为其它，造成Phenology数据对玉米的制图精度低于Red。对于图6(c1)、(c2)分类结果，NDVI数据错分其它类为水稻的像元多于NIR数据，这会造成前者水稻的用户精度低于后者；同理，图6(d1)、(d2)分类结果中，NIR数据错分其它类为玉米的像元多于NDVI+Phenology数据，造成前者玉米的用户精度低于后者。

本方法以辽宁省盘锦市和锦州市的部分区域为实验区，基于STDFA数据时空融合模型利用Landsat8和MODIS数据构建高时空分辨率数据(Red、NIR和NDVI)和通过TIMESAT软件从NDVI时间序列数据提取物候数据，然后分析各时间序列数据对秋粮作物的可分性，从中选取具有良好区分性的特征数据组成15种不同数据集用于秋粮识别，得出以下结论：

(1)基于STDFA模型构建的Red、NIR和NDVI数据及利用NDVI数据滤波提取的Phenology数据等4种基本类型数据可以有效区分秋粮作物的类型(图5)，在一定程度上可以解决秋粮识别过程中高分辨率数据不足的问题。由4种基本类型数据组成15种分类数据集对秋粮识别的各项精度均达到了较高的识别精度(表6)，水稻的制图精度和用户精度均分别达90％和83％左右；玉米的制图精度和用户精度也均分别达到85％和75％左右；二者的总体精度也达到85％左右。对表6中不同分类结果的各项精度进行综合分析，Red+Phenology数据组合对水稻和玉米识别的总体效果最好，水稻和玉米的制图精度、用户精度以及两者的总体识别精度分别达到了91.76％、82.49％、85.80％、74.97％和86.90％。同时，也在一定程度上验证了物候数据有助于秋粮识别。

(2)在分类数据的组合方式上，从分类结果精度中可以看出，秋粮作物的识别精度并不是随着分类数据的数据组合数量的增多而增加，二者之间并非呈正相关。但可以看出：随着分类数据的数据集组合数量的增多，同一精度评价项之间的差值越趋于稳定，当分类数据由三种或三种以上数据组合时，作物的各项识别精度不再发生大的变化，这说明不同数据的组合有助于秋粮作物识别的稳定性。

以上所述，仅是用以说明本发明的具体实施案例而已，并非用以限定本发明的可实施范围，举凡本领域熟练技术人员在未脱离本发明所指示的精神与原理下所完成的一切等效改变或修饰，仍应由本发明权利要求的范围所覆盖。

Claims

1.一种Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法，其特征在于所述方法包括以下步骤：

步骤一、数据准备，选取目标区域，获取同一区域的Landsat8影像数据、MODIS数据和航拍数据；

2.根据权利要求1所述的Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法，其特征在于：所述方法还包括步骤八、采用航拍数据评价分类结果，对选取目标区域拍影像矢量化数据转化的栅格数据进行重分类，使重分类后的分类结果与所述步骤七的分类数据进行像元对像元的叠加比较，得到用于分类结果精度评价的混淆矩阵。

3.根据权利要求1所述的Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法，其特征在于：所述步骤二数据预处理中对MODIS地表反射率产品首先转化为反射率数据，然后利用MODIS重投影转换成与Landsat8影像一致的坐标系，将MODIS产品的数据进行格式转换，并重采样到Landsat8数据分辨率的整数倍的空间分辨率，以便进行后续的MODIS混合像元分解。

4.根据权利要求3所述的Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法，其特征在于：所述步骤三中以MOD09Q1和Landsat8为基础数据利用数据时空融合STDFA模型构建三种高时空分辨率数据Red、NIR和NDVI，所述NDVI数据是利用时空数据融合STDFA模型构建的Red和NIR数据，通过下列公式计算得到：

5.如权利要求1所述的Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法，其特征在于：所述步骤四中在滤波处理之后还包括以Savitzky-Golay(S-G)拟合法为核心算法对时间序列数据进行重构建。

6.如权利要求1-5之一所述的Landsat8和MODIS融合构建高时空分辨率数据识别秋粮作物的方法，其特征在于：所述的航拍数据为高分辨率的无人机航拍数据，拍摄多个无人机航拍样方，对多幅样方进行拼接后转换坐标系，并且目视解译矢量化及地块属性赋值，以面积占优的类型转换成一定分辨率的栅格数据。