一种基于Landsat遥感数据的多云雨雾地区水稻识别方法与流程

本发明属于农业遥感
技术领域：
，具体涉及一种基于landsat遥感数据的多云雨雾地区水稻识别方法。
背景技术：
：水稻作为世界三大粮食作物之一，是人类最主要的粮食来源，在我国乃至世界的粮食生产结构中占有十分重要的地位。全球超过一半的人口以大米为主食，尤其对亚洲、非洲和拉美地区的发展中国家，粮食的安全生产显得尤为重要。水稻的生产状况与整个世界的粮食安全、社会稳定息息相关。掌握水稻种植面积、长势和产量信息，能够为监测中国水稻生产状况、指导农业生产与宏观调控水稻种植区划、水稻产量的预报和评估、粮食价格的预测和政府部门对粮食生产政策的制定等方面提供依据。长期以来，我国水稻种植面积依靠人工方法，通过实地抽样调查和逐级汇总方式来获取数据，这种方法不仅需要耗费大量的人力、物力，且受到各种主客观因素的影响，精度受到很大限制。随着遥感技术的快读发展，为农作物种植面积快速、准确地实现动态监测提供了新的技术手段。遥感信息具有覆盖范围大、探测周期短、现势性强、成本低等特点，有利于短时间内连续获取大范围的地面信息，实现农作物的种植面积提取。农作物种植面积的遥感提取离不开农作物的识别。而农作物的识别主要是利用绿色植物独特的波谱反射特征，从而将农作物与其他地物区分开。利用遥感技术进行水稻种植面积的估算，国内外已有大量的研究。以往的研究主要是通过单一时相的影响分类法对水田进行精细化监测或者通过时序归一化植被指数(ndvi)的差异监测水稻种植信息。近年来，随着新一代中、高空间分辨率卫星数据的出现，其多时相和多通道的优势，在监测水稻种植面积越来越受到重视。自2008年12月起美国政府决定将landsat全球数据全部共享，向所有用户免费提供。landsat是世界上收集时间最长的中等分辨率土地遥感数据，四十多年来为农业、地质、林业、区域规划、教育、测绘和全球变化研究提供了重要资源。landsat卫星数据的免费发布以及landsat-8和sentinel-2的推出为以更高空间分辨率的碎片化景观绘制水稻图谱提供了前所未有的机会。本发明的研究区为多云雨雾的华南西南地区，此地区的光学遥感数据可用的数量极少(如图3所示)，故本发明利用zhezhu的ccdc算法和时间序列模型来模拟、预测任意时间无云、无条带的landsat反射率数据。landsat卫星数据识别水稻主要的三种特征指数为ndvi(归一化植被指数)、(lswi)陆表水体指数和evi(增强型植被指数)。ndvi可较好的反映植被绿度变化，能够消除影像内部和外部的噪音。lswi是与植被水分含量相关的植被指数，利用对水体敏感的短波红外波段，对于处于蓄水期的水稻监测有较好的效果。evi利用蓝光波段修正大气对红光波段的影响，可提高对高生物量区的敏感度，与ndvi互为补充，现阶段应用卫星数据进行水稻种植面积的监测往往是根据水稻的物候历，确定水稻识别的移栽期、生长期和收获期等关键时期，通过关键时期内的特征来识别水稻。水稻生长过程中有三个重要的时期：一是移栽期；二是生长期；三是收获后。在不同生育期，随着水稻生长状况发生变化，相应地光谱特征也随之变化。目前这三种植被指数广泛应用于水稻遥感监测和估产研究中，landsat影像的最佳时相选取基于水稻不同时期的光谱特征。在移栽期，稻田常常存有2～15cm的水，此时地表是水稻和水体的混合，可通过影像中水体和水稻的混合光谱特征，利用对水体和植被较为敏感的波段或光谱指数来监测ndvi(evi)和lswi变化，识别蓄水和移栽期的水稻，并提取水稻种植面积。利用遥感影像高精度地提取作物种植面积，但是，当前研究仍然存在一些缺陷：(1)由于“同物异谱”和“异物同谱”现象以及混合像元影响了结果的准确性；(2)由于地区差异，不同种植类型水稻的耕作制度不同，不分段地一次性提取种植面积往往无法完全反映实际种植情况，导致种植面积提取结果与实际种植情况有出入。遥感影像分析方法的选取将直接影响水稻种植面积的提取精度，如何利用遥感影像实现高精度提取水稻的种植面积成为农业遥感领域亟待解决的技术问题。技术实现要素：针对现有水稻种植面积提取精度较低，无法反映实际种植情况的问题，本发明提供一种基于长时间序列landsat遥感数据的多云雨雾地区水稻识别方法，利用zhezhu的ccdc算法和时间序列模型来模拟、预测任意时间无云、无条带的landsat反射率数据，通过水稻的耕作制度和农时历作为辅助，分段逐步提取遥感影像中水稻种植面积，达到了提高提取精度的目的。本发明的技术方案具体如下：一种基于长时间序列landsat遥感数据的多云雨雾地区水稻识别方法，其特征在于，包括如下步骤：步骤1：首先从美国地质调查局(usgs)提供的espa网站(https://espa.cr.usgs.gov/)下载实验所需的landsat5/7/8二级产品数据(landsatlevel2)，其中包括地表反射率数据(surfacereflectance)、亮度温度数据(brightnesstemperature)和像元质量评估数据(pixelqa)，所有数据都是从2000年至今的云覆盖量小于80％。步骤2：然后根据云检测算法(fmask3.3，nasa所采用的云和云阴影检测算法)进行云和云阴影检测。然后利用ccdc算法将需要的多个波段的地表反射率数据与云和云阴影检测结果fmask波段组合在一起，生成envi标准格式。进一步将在linuxubuntu16.04操作系统上利用并行方式运行ccdc_changeard13_01(一个可执行文件),这一步ccdc算法的关键步骤——连续变化检测。最后利用zhezhu的时间序列模型获取研究区域内任意时刻的无云、无条带的landsat反射率数据集，此数据集的每一景数据共计8个波段，分别是蓝光波段、绿光波段、红光波段、近红外波段、短波红外1、短波红外2、热红外和fmask(blue,green,red,nir,swir1,swir2,tir,fmask)。此步骤所用的时间序列模型是由zhezhu和woodcock在2015年提出的，该模型的基本思想是：基于所有可用的landsat反射率数据，利用每个像元的至少15个有效观测数据拟合一个曲线方程，将该曲线方程作为时间序列模型进行得到任意时刻的landsat反射率数据(如图6所示)。时间序列模型所采用的曲线方程如下：上述公式(1)中，x表示julian日期(儒略日)，i表示landsat影像的第i个波段，t为变化周期(t＝365),a0,i是第i个波段的均值系数，a1,i，b1,i是年内变化系数，c1,i为年际变化系数，为第i个波段在x时刻的模型预测值，表示第k个时间节点。步骤3：根据步骤2得到的任意时刻的无云、无条带的landsat反射率数据集，计算其光谱指数——ndvi(归一化植被指数,式2)、evi(增强型植被指数、式3)、lswi(陆表水体指数，式4)，利用多种光谱指数的不同组合来识别永久水体(式5)、常绿植被(式6)等其它地物类型将水稻与其他明显地物区分，确定研究区域内潜在的水稻种植区域。具体表达式如下所示：ndvi<0.1&&ndvi<lswi(10/46)(5)ndvi>0.7(20/46)||lswi>0.15(40/46)(6)上述公式2、3、4中nir、red、blue分别代表近红外波段、红光波段、蓝光波段、短波红外1。其中式5中的(10/46)表示一年46景数据中每个像元符合大于等于10即为永久水体；式6中与式5表示意思相同，一年46景数据中每个像元满足ndvi>0.7符合大于等于20或者一年46景数据中每个像元满足lswi>0.15符合大于等于40即为常绿植被。步骤4：利用geoexplore3000仪器、地物光谱仪(psr-3500)、lai-2200植物冠层分析仪获取地面确定的大区域的水稻实测数据，其中实测数据主要包括提取水稻种植区域的经纬度信息、水稻不同生长阶段的光谱曲线以及lai数据。然后根据实测的水稻种植区域的经纬度信息提取对应像元移栽期阶段的步骤3中的光谱指数数据(ndvi、evi、lswi)，然后利用数理统计分析的方法分析移栽期前后水稻的不同光谱指数之间的关系(如图8所示)，结合水稻种植区划数据和水稻的潜在移栽期数据获取的水稻相关信息，确定植被指数阈值模型的参数t，建立植被指数阈值模型。植被指数阈值模型表达式如公式7所示：lswi+t≥evi或lswi+t≥ndvi(7)步骤5：根据步骤4中建立的植被指数阈值模型(式7)对步骤3中确定的研究区域内潜在的水稻种植区域进行逐像元的进行判别，如果某像元在农业部市场与经济司所公布水稻潜在移栽期(如图5所示)前后7～15天中有一次或多次满足植被指数阈值模型，即判定为水稻；否则判定为非水稻。根据上述判定结果，获取最终的水稻种植区域分布图。步骤6：结合国家统计局统计年鉴中水稻种植面积数据(http://data.stats.gov.cn/easyquery.htm？cn＝c01)，对提取的水稻种植面积结果进行对比分析，或者利用谷歌地球提供的高分辨率影像数据(veryhighresolutionimages)选择随机样本点进行对比分析，精度评价。如果精度较低，重复步骤4和步骤5，对参数t的阈值进行相应的调整，确定最为合适的植被指数阈值模型。本发明所述遥感影像产品具体为美国地质调查局(usgs)提供的espa数据，即landsatlevel2(landsat5/7/8二级产品)的产品数据(https://espa.cr.usgs.gov/)。进一步地，所述步骤2中任意时刻的无云、无条带的landsat反射率数据集是利用的zhezhu和woodcock在2014年提出的连续变化检测和分类算法(continuouschangedetectionandclassification，ccdc)以及后期提出的时间序列模型，基于美国地质调查局提供的可用的云覆盖量低于20％的landsat影像数据获取的。进一步地，所述步骤3中将水稻与其他明显地物区分是通过包括lswi、ndvi、evi在内的植被指数之间的差值关系，确定永久水体、常绿植被等其它地物类型，从而获取水稻潜在种植区域。进一步地，所述步骤4中水稻种植区划数据是根据国家统计局自2001年开始出版的《中国统计年鉴》中不同种植类型的水稻面积整理分析而确定。进一步地，所述步骤4中水稻的潜在移栽期数据是根据农业部市场与经济司所公布的各地区的水稻农时信息整理分析而确定。进一步地，所述步骤4中植被阈值模型是依据地面实测水稻种植地区的移栽期间的lswi和evi(ndvi)的散点图(如图8所示)，建立lswi和evi(ndvi)之间的相关关系，进而实现水稻识别。进一步地，所述步骤6中还包括：结合统计年鉴数据和高分辨率影像数据，对提取的不同类别的水稻种植面积进行精度评价，如果精度达不到要求，可重复步骤4和步骤5，对植被指数阈值模型中阈值参数t进行调整。本发明考虑到水稻与其它地物类型的区分，剔除干扰像元的影响，并且选取水稻的移栽期作为水稻识别的关键生长期，移栽期土壤的含水量非常高。因此可根据稻田含水量高的特点，从遥感图像上与其他农作物区别，实现水稻的识别。lswi是与植被水分含量相关的植被指数，evi对高生物量区非常敏感，因此利用通过地面实测水稻移栽期的样本数据来建立lswi和evi(ndvi)之间的相关关系来建立植被阈值模型。同时，考虑到不同种植类型水稻的耕作制度不同，进一步地结合水稻的耕作制度和农时信息，采用植被阈值模型进行分段提取水稻种植面积，最终将不同种植类型的水稻种植面积合并，即可准确反映水稻种植情况，得到相对较为准确的水稻种植面积。下面阐述本发明的构思：为考虑到区分水稻与其它地物类型，本研究选取水稻的移栽期作为识别的关键生长期，识别移栽期的水稻，最终达到水稻识别与种植面积提取的目的。与现有技术相比，本发明提供的一种基于长时间序列landsat遥感数据的多云雨雾地区水稻识别方法，达到如下有益效果：第一，本发明利用光学遥感数据快速实现大面积水稻制图，深度挖掘了光学遥感数据在农业方面的应用前景，也为科学的指导农事提供了可靠的依据。第二，与现有的基于光学遥感的水稻制图的区域相比，本发明主要适用区域为多云雨雾的华南西南地区。第三，与传统依靠野外测量调查水稻种植面积的方法相比，本发明节省了大量的人力物力财力，能够真实反映实际种植情况，提高了种植面积提取的精度。附图说明图1是本发明实施例基于长时间序列landsat遥感数据的多云雨雾地区水稻识别方法的流程示意图。图2是本发明实施例中的研究区域图。图3是本发明实施例中研究区内可用的landsat时间序列清晰像元占比分布图(行/列:126/041)。图4是本发明实施例中研究区域内水稻种植区划信息图。图5是本发明实施例中研究区域内水稻的潜在移栽期图。图6是本发明实施例中研究区域内landsat的第四波段反射率数据的真实值与时间序列模型得到的预测曲线图。图7(a)是本发明实施例中研究区域内有云、有条带的landsat数据，图7(b)是本发明实施例中研究区域内通过时间序列模型得到的无云、无条带的landsat数据。图8是本发明实施例中地面实测稻田灌水移栽期lswi和evi(ndvi)的关系图。图9是本发明实施例中基于长时间序列landsat遥感数据的多云雨雾地区水稻识别方法提取得到的2010年至2018年的水稻种植分布图。具体实施方式通过选择水稻生长阶段的关键时期——移栽期作为水稻识别的突破口，在移栽期，稻田常常存有2～15cm的水，因此地表是水稻和水体的混合，可通过遥感影像中水体和水稻的混合光谱特征，利用对水体和植被较为敏感的波段或植被指数来监测ndvi、evi和lswi变化，识别移栽期的水稻，并提取种植面积。ndvi、evi、lswi是较为常见的植被指数，可通过遥感卫星反射率产品(landsatlevel2)的波段计算获取。下面结合说明书附图和具体实施例对本发明进行详细说明：实施例：一种基于长时间序列landsat遥感数据的多云雨雾地区水稻识别方法，包括如下步骤：步骤1：数据获取：依据美国地质调查局(usgs)提供的遥感卫星反射率产品，获取landsat5/7/8二级反射率产品；依据国家统计局发布的《中国统计年鉴》和农业部市场与经济司所公布的各地区的农作物物候信息，获取水稻种植区划数据和水稻潜在移栽期数据；依据地面实测数据，获取水稻灌水移栽期lswi和evi(ndvi)数据。步骤2：数据处理：对步骤1获取的landsat5/7/8二级反射率产品利用ccdc算法依次进行云和云阴影检测、连续变化检测，然后利用时间序列模型结合变化检测结果得到任意时刻无云、无条带的landsat反射率数据(具体流程如图1所示)；对获取的《中国统计年鉴》中水稻种植面积进行整理、分析，得到水稻种植区划数据，结果如图4所示；对获取农业部市场与经济司所公布的各地区的农时历整理分析得到水稻的潜在移栽期数据，结果如图5所示；对获取的无云、无条带的landsat反射率数据计算所需的归一化植被指数ndvi、增强型植被指数evi、陆表水体指数lswi。对地面实测水稻淹水移栽期的lswi和evi(ndvi)统计分析得到水稻识别模型，如图8所示。此步骤所用的时间序列模型是有zhezhu和woodcock在2015年提出的，该模型的基本思想是：基于所有可用的landsat影像，利用每个像元的多个有效观测数据拟合一个曲线方程，将该曲线方程作为时间序列模型进行得到任意时刻的landsat反射率数据(如图6所示)。时间序列模型所采用的曲线方程如下：上述公式(1)中，x表示julian日期(儒略日)，i表示landsat影像的第i个波段，t为变化周期(t＝365),a0,i是第i个波段的均值系数，a1,i，b1,i是年内变化系数，c1,i为年际变化系数，为第i个波段在x时刻的模型预测值，表示第k个时间节点。步骤3：潜在水稻区域的确定：该步骤的原理具体是通过对包括ndsi、lswi、ndvi等植被指数之间的差值关系确定永久水体、常绿植被、雪等其它地物类型，获取水稻可能或潜在的区域；本发明基于光学遥感数据的多时相特征，以年作为时间尺度，将对2018年landsat反射率数据计算获取lswi、ndvi、evi等指数数据，来识别永久水体(式2)、常绿植被(式3)等其它地物类型。具体表达式如下所示：ndvi<0.1&&ndvi<lswi(10/46)(5)ndvi>0.7(20/46)||lswi>0.15(40/46)(6)式2中(10/46)表示一年46景数据中每个像元符合大于等于10即为永久水体；式2中与式3表示意思相同，一年46景数据中每个像元满足ndvi>0.7符合大于等于20或者一年46景数据中每个像元满足lswi>0.15符合大于等于40即为常绿植被。步骤4：水稻种植面积的提取：(1)水稻识别：在移栽期，稻田的反射光谱通常都是水、土壤、秧苗及稻田的背景沟渠、道路、杂草、防护林和其他作物等的混合光谱，此时，稻田高土壤含水量和低植被覆盖度可以利用lswi和evi检测出来。具体的检测原理如下：如果在灌水移栽期时evi值较高，则说明该像素代表的地物为其他植被,如树木、灌木、草地或者其他作物等，因此就可认为是非水稻区；如果lswi很低，则说明土壤含水量较低的区域，同样可认为是非水稻区；反之，如果lswi较高并且evi较低，那么该像元很有可能就是移栽期的水稻田；本实施例中包括对获取的地面实测水稻样本数据的evi(ndvi)和lswi这两种植被指数进行统计分析，建立植被指数evi(ndvi)和lswi之间的相关关系；为了检测稻田的光谱特性,在研究区域范围内选择35个地面实测点，以2018年的遥感影像数据作为分析依据，根据从以上实测点获取的数据，计算各个实测点覆盖下稻田灌水移栽期的平均evi(ndvi)和lswi，稻田灌水移栽期evi(ndvi)和lswi表现出的特征如图8所示；通过以上实测点获取的资料分析结果可得到水稻的提取算法：lswi+t≥evi或lswi+t≥ndvi(7)(2)水稻种植空间分布获取：以2019年在该研究区地面实测水稻点以及多次随机选取的谷歌地球高分辨率的样本点为参考，进行精度评价。下面是精度对比分析表：样本点总数水稻点个数正确率水稻地面实测点403895％谷歌地球水稻样本点(1)1008383％谷歌地球水稻样本点(2)1008686％谷歌地球水稻样本点(3)1008585％谷歌地球水稻样本点(4)1008484％谷歌地球水稻样本点(5)1008787％上表表明，本方法在多云雨雾地区，30米空间分辨率的水稻制图精度大约85％，这说明了本方法的可行性。本实施例采用上述提出的水稻识别模型进一步对landsat数据进行2010至2018年最终的水稻种植分布图，结果如图9所示。步骤5：调整算法：在推广的使用本发明的过程中，由于应用区域水稻生长状态的不同，反映植被指数阈值模型参数上也略有不同，需要对参数t的阈值进行上下的调整。当前第1页12