基于植被指数时空间统计特征的农业灾害信息遥感提取方法与流程
本发明涉及基于植被指数时空间统计特征的农业灾害信息遥感提取方法，属于农业灾害信息获取技术领域。
背景技术：
目前大多数农业灾害遥感监测的研究中，对小面积地区或某一时期、某一特定类型灾害的研究较多，对于大尺度区域或长时间序列的灾害监测与评价研究较少，不能全面且宏观地掌握和分析农业灾害的时空分布规律；由于不同作物类型、生长地区、生长阶段的植被指数存在差异，利用单时相影像所获得的植被指数，仅能反映在该地区、该时间内作物长势的相对优劣，进行灾害监测不具有普遍性。
MODIS全称为中分辨率成像光谱仪(moderate-resolution imaging spectroradiometer)。MODIS是搭载在Terra和Aqua卫星上的一个重要的传感器，是卫星上唯一将实时观测数据通过x波段向全世界直接广播，并可以免费接收数据并无偿使用的星载仪器，全球许多国家和地区都在接收和使用MODIS数据。MODIS用于对陆表、生物圈、固态地球、大气和海洋进行长期全球观测。中分辨率成像光谱仪归一化植被指数MODIS-NDVI时序数据空间分辨率较低，但时间分辨率高，覆盖范围广，适用于大尺度范围、长时间序列监测。
高分一号卫星影像和环境减灾卫星影像空间分辨率高，已被证实可以满足灾害监测。高分一号GF-1卫星搭载了两台2m分辨率全色/8m分辨率多光谱相机，四台16m分辨率多光谱相机。高分一号卫星的宽幅多光谱相机幅宽达到了800公里。环境与灾害监测预报小卫星星座A、B、C星HJ-1A/B/C包括两颗光学星HJ-1A/B和一颗雷达星HJ-1C，可以实现对生态环境与灾害的大范围、全天候、全天时的动态监测。
技术实现要素：
本发明目的是为了解决现有农业灾害遥感监测不能用于对大尺度区域或区域内长时间序列的灾害进行监测与评价，监测方法不具有普遍性的问题，提供了一种基于植被指数时空间统计特征的农业灾害信息遥感提取方法。
本发明所述基于植被指数时空间统计特征的农业灾害信息遥感提取方法，它包括以下步骤：
步骤一：采集被监测区域作物生长期的MODIS反射率产品MOD09Q1数据和MYD09Q数据；采集被监测区域基准年全年的植被指数产品MOD13Q1数据；采集被监测区域基准年内典型受灾区域所需的HJ-1A/1B_CCD影像数据和GF-1/WFV影像数据；
步骤二：对所述MOD09Q1数据和MYD09Q数据转投影，分别计算获得归一化植被指数MOD09Q1-NDVI和MYD09Q1-NDVI，将二者拼接结合，获得被监测区域作物生长期的MOD_MYD-NDVI时间序列；
对所述MOD13Q1数据进行转投影，提取被监测区域基准年全年的23期归一化植被指数MOD13Q1-NDVI时间序列；
步骤三：对所述HJ-1A/1B_CCD影像数据和GF-1/WFV影像数据分别进行辐射定标、大气校正、正摄校正和自动匹配项目的预处理，提取获得所述典型受灾区域的基准年植被覆盖指数HJ-NDVI；用已知的该典型受灾区域的矢量文件裁剪所述基准年植被覆盖指数HJ-NDVI，并对裁剪后获得的HJ-NDVI数据进行分级处理；将分级后的HJ-NDVI数据与原始的HJ-1A/1B_CCD影像数据和GF-1/WFV影像数据对比，去除分级后的HJ-NDVI数据中非受灾的等级，保留受灾的等级作为灾害区域HJ-NDVI数据，用灾害区域HJ-NDVI数据裁剪所述的分级后的HJ-NDVI数据，再重新分级获得最终灾害监测基准区；
步骤四：对归一化植被指数MOD13Q1-NDVI时间序列进行平滑重构，在被监测区域内，利用基准年全年的植被指数产品MOD13Q1数据，根据作物的物候特征，进行S-G滤波平滑，提取11个物候参数；对11个物候参数进行标准化处理，统一数值范围，然后提取每个物候参数的有效成分，再对被监测区域进行分区获得物候分区；
步骤五：将被监测区域与各物候分区的矢量文件分别进行交汇处理，获得包含物候信息及作物信息的矢量文件，为该矢量文件添加包含物候信息和作物类型信息的字段，获得包含物候信息和作物类型信息的矢量文件；
步骤六：对步骤五中获得的包含物候信息和作物类型信息的矢量文件进行提取，获得各个物候分区各类作物的归一化植被指数中值NDVIm与标准差STD；
步骤七：对步骤三中所述最终灾害监测基准区的受灾年份与未受灾年份的NDVI灰度直方图进行分析，获得作物受灾的判断依据：当某时相某区域作物的实际NDVI值小于NDVIm-xSTD，并且在连续的下两个时相内仍然小于NDVIm-xSTD时，认定该区域作物受灾；由此，建立作物受灾监测模型：NDVI-(NDVIm-xSTD)；式中x为受灾监测模型的阈值；
步骤八：由步骤二中的被监测区域作物生长期的MOD_MYD-NDVI时间序列和步骤七中作物受灾监测模型，在0～1范围内尝试不同阈值，通过空间分析法提取被监测区域的灾害范围，将被监测区域分为6月下旬至7月中旬、7月中旬至8月上旬、8月上旬至八月下旬三个时段，每个时段分为三个时相，计算连续三个时相都受灾的栅格数据；
步骤九：将连续三个时相都受灾的栅格数据转成矢量数据，计算受灾面积，将大于三个像元面积的区域认作受灾区域的中心，再将受灾区域向外扩展一个像元的大小，获得新的矢量文件，用新的矢量文件裁剪之前得到的连续三个时相都受灾的栅格数据，得到不同阈值对应的灾害监测结果；
步骤十：将不同阈值对应的灾害监测结果与相应的已知最终灾害监测基准区的监测结果进行对比，将每一对两种不同分辨率监测结果分别转成矢量数据，进行交汇处理，得到同时包含两种分辨率监测结果的矢量文件；取阈值＝0.5作为最优阈值，获得最终的灾害监测模型为NDVI-(NDVIm-0.5STD)；
步骤十一：将最终的灾害监测模型NDVI-(NDVIm-0.5STD)用于被监测区域，获得被监测区域的农业灾害信息。
本发明的优点：本发明方法能够快速、准确的获取农业灾害信息，它在考虑到不同物候区、作物和生长阶段的植被指数存在差异的基础上，提取已知灾害所在物候区内的各像元NDVI以及相同作物的NDVI平均值NDVIm和标准差STD，根据受灾前后NDVI灰度直方图特征，利用统计学特征分析各参数之间的关系，建立灾害监测模型，提取农业灾害。该方法考虑到了因生长地区、不同作物和生长阶段造成的干扰因素，提高了监测结果精度。
本发明提取作物的NDVI、NDVIm、STD，通过已知典型灾害，根据受灾前后灾害所在物候区的NDVI灰度直方图，分析验证各参数之间的关系，建立灾害监测模型，并确定了模型的最优阈值，得出最终的灾害监测模型，基于不同物候区与作物分类进行灾害信息提取，解决了因为地域不同造成的作物长势差异和不同作物的长势差异等影响因素对提取结果造成的误差，提高了精度。将每年作物关键生长期分为177-193、193-209、209-225三个时段，分别对黑龙江省2013年至2015年投保地块提取灾害空间信息，验证了最终的灾害监测模型的通用性，实现了大尺度范围、长时间序列农业灾害监测。
附图说明
图1是2013年6月末四方山涝灾HJ灾害监测图，其中1级为受灾最重的等级，逐级变轻；
图2是2013年6月末四方山涝灾MODIS灾害监测图；
图3是2015年6月24日洪河农场雹灾HJ灾害监测图，其中1级为受灾最重的等级，逐级变轻；
图4是2015年6月24日洪河农场雹灾MODIS灾害监测图；
图5是对2014年黑龙江省投保地块采用本发明方法进行MODIS灾害监测获得的监测图；图中图例由上至下依次为黑龙江省省界、黑龙江省县界、177-193时段(6月下旬至7月中旬)MODIS农业灾害监测结果、193-209时段(7月中旬至8月上旬)MODIS农业灾害监测结果、209-225时段(8月上旬至8月下旬)农业灾害监测结果；
图6是对2015年黑龙江省投保地块采用本发明方法进行MODIS灾害监测获得的监测图。图中图例由上至下依次为黑龙江省省界、黑龙江省县界、177-193时段(6月下旬至7月中旬)MODIS农业灾害监测结果、193-209时段(7月中旬至8月上旬)MODIS农业灾害监测结果、209-225时段(8月上旬至8月下旬)农业灾害监测结果。
具体实施方式
具体实施方式一：下面结合图1至图6说明本实施方式，本实施方式所述基于植被指数时空间统计特征的农业灾害信息遥感提取方法，它包括以下步骤：
步骤一：采集被监测区域作物生长期的MODIS反射率产品MOD09Q1数据和MYD09Q数据；采集被监测区域基准年全年的植被指数产品MOD13Q1数据；采集被监测区域基准年内典型受灾区域所需的HJ-1A/1B_CCD影像数据和GF-1/WFV影像数据；
步骤二：对所述MOD09Q1数据和MYD09Q数据转投影，分别计算获得归一化植被指数MOD09Q1-NDVI和MYD09Q1-NDVI，将二者拼接结合，获得被监测区域作物生长期的MOD_MYD-NDVI时间序列；
对所述MOD13Q1数据进行转投影，提取被监测区域基准年全年的23期归一化植被指数MOD13Q1-NDVI时间序列；
步骤三：对所述HJ-1A/1B_CCD影像数据和GF-1/WFV影像数据分别进行辐射定标、大气校正、正摄校正和自动匹配项目的预处理，提取获得所述典型受灾区域的基准年植被覆盖指数HJ-NDVI；用已知的该典型受灾区域的矢量文件裁剪所述基准年植被覆盖指数HJ-NDVI，并对裁剪后获得的HJ-NDVI数据进行分级处理；将分级后的HJ-NDVI数据与原始的HJ-1A/1B_CCD影像数据和GF-1/WFV影像数据对比，去除分级后的HJ-NDVI数据中非受灾的等级，保留受灾的等级作为灾害区域HJ-NDVI数据，用灾害区域HJ-NDVI数据裁剪所述的分级后的HJ-NDVI数据，再重新分级获得最终灾害监测基准区；
步骤四：对归一化植被指数MOD13Q1-NDVI时间序列进行平滑重构，在被监测区域内，利用基准年全年的植被指数产品MOD13Q1数据，根据作物的物候特征，进行S-G滤波平滑，提取11个物候参数；对11个物候参数进行标准化处理，统一数值范围，然后提取每个物候参数的有效成分，再对被监测区域进行分区获得物候分区；
步骤五：将被监测区域与各物候分区的矢量文件分别进行交汇处理，获得包含物候信息及作物信息的矢量文件，为该矢量文件添加包含物候信息和作物类型信息的字段，获得包含物候信息和作物类型信息的矢量文件；
步骤六：对步骤五中获得的包含物候信息和作物类型信息的矢量文件进行提取，获得各个物候分区各类作物的归一化植被指数中值NDVIm与标准差STD；
步骤七：对步骤三中所述最终灾害监测基准区的受灾年份与未受灾年份的NDVI灰度直方图进行分析，获得作物受灾的判断依据：当某时相某区域作物的实际NDVI值小于NDVIm-xSTD，并且在连续的下两个时相内仍然小于NDVIm-xSTD时，认定该区域作物受灾；由此，建立作物受灾监测模型：NDVI-(NDVIm-xSTD)；式中x为受灾监测模型的阈值；
步骤八：由步骤二中的被监测区域作物生长期的MOD_MYD-NDVI时间序列和步骤七中作物受灾监测模型，在0～1范围内尝试不同阈值，通过空间分析法提取被监测区域的灾害范围，将被监测区域分为6月下旬至7月中旬、7月中旬至8月上旬、8月上旬至八月下旬三个时段，每个时段分为三个时相，计算连续三个时相都受灾的栅格数据；
步骤九：将连续三个时相都受灾的栅格数据转成矢量数据，计算受灾面积，将大于三个像元面积的区域认作受灾区域的中心，再将受灾区域向外扩展一个像元的大小，获得新的矢量文件，用新的矢量文件裁剪之前得到的连续三个时相都受灾的栅格数据，得到不同阈值对应的灾害监测结果；
步骤十：将不同阈值对应的灾害监测结果与相应的已知最终灾害监测基准区的监测结果进行对比，将每一对两种不同分辨率监测结果分别转成矢量数据，进行交汇处理，得到同时包含两种分辨率监测结果的矢量文件；取阈值＝0.5作为最优阈值，获得最终的灾害监测模型为NDVI-(NDVIm-0.5STD)；
步骤十一：将最终的灾害监测模型NDVI-(NDVIm-0.5STD)用于被监测区域，获得被监测区域的农业灾害信息。
步骤四中的11个物候参数为作物生长起始期、作物生长结束期、振幅、NDVI的平均值、生长期长度、NDVI的积分、NDVI最大值、左侧上升曲线间斜率、右侧下降曲线间斜率、整个时期的中间点和整个时期NDVI的积分。
一个像元的面积为250×250m2。
本发明使用的是30m空间分辨率、2d时间分辨率的HJ-1A/1B_CCD影像，16m空间分辨率、4d时间分辨率的GF-1/WFV影像、250m空间分辨率、8d时间分辨率的MODIS影像、250m空间分辨率、16d时间分辨率的MODIS影像。
步骤四中，获取11个物候参数采用的S-G滤波平滑手段可通过Timesat软件实现；对11个物候参数统一数值范围是为了减少较大数值在综合评价中的影响作用，对其进行主成分分析，提取有效成分，避免数据冗余，使影像颜色更加突出明显。物候分区是采用面向对象的方法进行基于物候信息的多尺度分割，经多次试验，得到2015年最优尺度为70，该尺度下分割效果良好且不破碎，能够充分体现物候影像的差异性，故70为最优尺度。，得到物候分区结果，避免不同物候区造成的气候、土壤、水分等干扰因素。
步骤五中，为矢量文件添加包含物候信息和作物类型信息的字段的具体方法是，在文件的属性表中新添加一列字段，通过field calculator使这列新的字段同时包含物候信息和作物类型信息。
步骤六中，获得各个物候分区各类作物的归一化植被指数中值NDVIm与标准差STD，能够排除因为不同物候区和不同作物造成的NDVI差异的影响。
步骤七中，在ARCGIS软件中分析已知典型灾害所在物候区受灾年份与未受灾年份的NDVI灰度直方图，发现未受灾年份的NDVI灰度直方图大致以NDVIm为中心呈正太分布，受灾年份的NDVI灰度直方图呈偏态分布，且偏向NDVI值较小的方向，故当某时相某区域作物的实际NDVI值小于NDVIm-xSTD，并且在连续的下两个时相仍然小于NDVIm-xSTD时，认定作物受灾。
步骤八中，对被监测区域划分的关键生长期三个时段的对应关系分别为177-193对应于6月下旬至7月中旬、193-209对应于7月中旬至8月上旬、209-225对应于8月上旬至八月下旬三个时段，利用步骤七的受灾监测模型，利用ARCGIS软件的空间分析功能分别提取2013年黑龙江177-225不同时期全省主要作物：玉米、水稻、大豆和麦类等农业灾害空间信息。将相同作物的每个时段三个时相的受灾像元在ARCGIS软件中利用calculator功能得到连续三个时相都受灾的栅格数据。利用raster to polygon功能转成矢量文件。
步骤九中，由于MODIS影像的空间分辨率较低，会出现混合像元的干扰，因此认定至少大于3个MODIS像元的面积才算受灾；将连续三个时相都受灾的栅格数据转成矢量数据，计算面积，导出大于三个像元，即2502×3＝127500m2面积的区域，认作受灾范围的中心，再将受灾区域向外buffer一个像元，250m的大小，得到新的矢量文件。
步骤十中，根据与30个已知典型灾害对比分析，再在全黑龙江省投保地块范围内随机选取30个非受灾点，利用窗口变步长法不断搜索阈值，因为作物受灾后其NDVI变小，故阈值x在0～1的范围内，计算不同阈值对应的kappa系数PA％和总体精度Pc。最终当阈值x为0.5时，此时的PA％为80％，Pc为0.9，确定了最终的灾害监测模型。阈值＝0.5时，包含两种分辨率监测结果的矢量文件中包含两种分辨率监测结果的典型灾害最多，因此作为最优阈值。
具体实施例如下：
步骤一：下载某研究区所需的MODIS反射率产品MOD09Q1和MYD09Q1，对作物生长季的MOD09Q1和MYD09Q1数据进行转投影，利用空间分析功能计算NDVI，二者通过拼接处理结合使用，消除云影响。得到研究区2013年5月初到10月末的归一化植被指数NDVI时间序列。
步骤二：下载研究区2013年已知典型灾害区域的HJ-1A/1B_CCD或GF-1/WFV无云影像，进行辐射定标、大气校正、正摄校正、自动匹配等预处理，提取NDVI，分级处理，保留灾害区域，得到最终的高分辨率灾害监测结果。
步骤三：利用物候数据进行信息提取，将11个物候数据进行主成分分析及标准化处理。对处理后的数据采用面向对象的方法进行基于物候信息的多尺度分割，确定最优尺度，得到物候分区结果，避免不同物候区造成的气候、土壤、水分等干扰因素。
步骤四：将耕地与物候区进行交汇，使新得到的矢量文件属性表中既含有物候又含有作物信息。
步骤五：利用Arcgis中分区统计功能提取不同物候区各类作物的NDVI中值NDVIm与标准差STD。
步骤六：通过分析已知典型灾害所在物候区受灾前后的NDVI灰度直方图，得出当某时相某区域作物的实际NDVI值小于NDVIm-xSTD，并且在连续的下两个时相仍然小于NDVIm-xSTD时，认定作物受灾。依据这个原理，建立初步监测模型。
步骤七：通过窗口变步长法，调整步长搜寻最优阈值，得出Kappa系数和总体精度最高时即为最优阈值。表1为窗口变步长最优阈值搜寻结果，得出最优阈值为0.5。
表1 窗口变步长最优阈值搜寻结果
步骤八：通过步骤六确定的最优阈值，得出最终的灾害监测模型。利用此模型提取2014年和2015年全黑龙江省农业灾害空间信息，每年分为177-193、193-209、209-225三个时段，与14个已知典型灾害对比分析，表2为2014年和2015年上报灾害详情及提取结果。总体精度Pc为0.93，Kappa系数PA％为85.9％，详情见表3，精度较高，证明了该监测模型的通用性，可适用于大尺度范围、长时间序列灾害空间信息监测。
表2 2014年和2015年上报灾害详情及提取结果
注：表2中提取结果为是或否，表示MODIS灾害监测结果是否提取出已知典型灾害。
表3 灾害监测模型精度评价表
- 还没有人留言评论。精彩留言会获得点赞！