一种集成像元与对象的水稻自动变化检测方法
【专利摘要】本发明涉及一种集成像元与对象的水稻自动变化检测方法，包括以下步骤：建立基于MODIS水体指数时序数据集及MODIS水体指数时序曲线；逐像元计算起始与结束年份MODIS水体指数时序曲线的相似度；基于MODIS水体指数时序曲线的相似性，判断疑似变化区域；综合起始年份多期Landsat影像，通过影像分割获得影像对象；将疑似变化区域与影像对象叠加，获得研究区疑似变化对象；综合结束年份多期Landsat影像，将研究区疑似变化对象进一步细分割；对细分割后的研究区疑似变化对象，利用Landsat遥感影像特征，分别开展起始和结束年份的水稻制图，从而获得研究区起始到结束年份水稻动态变化分布图。本发明具有自动化程度高、简单易用、鲁棒性好以及分类精度高等特点。
【专利说明】
一种集成像元与对象的水稻自动变化检测方法
技术领域
[0001 ]本发明涉及一种集成像元与对象的水稻自动变化检测方法。
【背景技术】
[0002]在遥感应用研究领域，早期的研究成果主要集中在土地利用及其变化监测方面。 21世纪以来，随着遥感影像的时间与空间分辨率的不断提高，遥感应用的深度与广度不断 拓展。近年来，利用遥感影像开展农作物监测的研究备受关注。但是，由于不同农作物之间 的光谱相似性，以及同种农作物在不同区域内受物候变化、品种差异等各种因素导致的光 谱不一致性，农作物遥感监测面临挑战。常规的遥感监测方法难以满足实时快速自动监测 农作物种植面积变化的实际需要，急需建立新型高效的自动半自动变化检测方法。
[0003]针对农作物遥感监测面临的难题，目前主导方向与实施策略是大力发展基于时序 遥感影像的分类方法。即：通过连续多期的遥感影像，监测农作物在播种、发芽、开花、结实 乃至成熟的整个生长周期内的变化，从而达到农作物遥感分类的目的。近年来，基于时序遥 感影像的分类技术发展迅速，并且在农作物遥感监测中取得了一定的成效。相关研究成果 主要集中在针对大宗农作物如水稻、小麦等的种植面积遥感监测方面。
[0004] 在水稻遥感监测方面，比较常用的方法有Xiangming Xiao等（2005)提出的基于植 被指数与水体指数的差值的方法。该方法的判断依据是，水稻在移栽期间植被指数通常比 较低，而水体指数受到移栽期灌水的影响急剧增高，如果水体指数与植被指数之差大于一 定的阈值，判断为水稻。该方法简便有效，在东南亚、中国等水稻主产区得到了很好的应用。 邱炳文等(2015)针对水体指数易受到降水等因素干扰的特点，提出了一种基于特定物候期 的水体与植被指数变化比值指数的水稻自动制图方法。该方法充分利用特定物候期内水体 指数与植被指数变化量的比值，具备较好的抗干扰能力。
[0005] 诸多研究表明，基于时序遥感影像的分类方法具有很好的应用前景。基于时序遥 感影像的分类方法的应用前提条件是，能获得研究时段内比较理想的连续的时序遥感影像 数据。虽然近年来遥感影像数据源日益丰富，但遥感影像的时间、空间分辨率不可兼得的矛 盾将在很长一段时间内存在。一方面，高时间分辨率的卫星遥感数据的空间分辨率通常比 较低，如常见的M0DIS遥感影像数据;另一方面，中等或较高空间分辨率的卫星遥感影像数 据，通常时间分辨率偏低，如常用的Landsat卫星数据。因此在建立中等或较高空间分辨率 遥感影像的时序数据集的过程中，通常苦于雨雾霾等天气因素的影响，难以建立数据时间 步长比较短、影像质量比较好的、完整的时序数据集。
[0006] 构建较高的时间与空间分辨率的时序数据集中所遇到的种种障碍，一定程度上阻 碍了基于时序遥感影像的分类方法及其应用的深入发展。即使某些遥感影像数据兼具较高 的空间与时间分辨率的特点，其影像处理与计算的工作量倍增，如能有效地将不同分辨率 遥感影像进行结合不失为一种有效途径。
[0007] 依据分类对象的不同，可以将遥感影像分类方法分为两种:基于像元以及面向对 象的分类方法。对于比较粗分辨率的遥感影像，通常采用基于像元的分类方法。但基于像元 的遥感影像分类方法运用到较高分辨率遥感影像时，却遇到一系列问题。比如，随着遥感影 像空间分辨率的提高，像元与完整地物之间很难一一对应，通常多个像元组合起来才对应 一种地物。因此利用面向对象分类方法，通过将遥感影像分割为一系列对象的基础上，结合 形状和纹理等信息进行分类，不失为较高分辨率遥感影像分类的合适方法。如果将面向对 象分类方法运用到时序遥感影像，将难以处理利用某一期遥感影像分割建立的对象内部的 某些像元的类别已经发生了变化这类问题。因此基于时序遥感影像的分类方法通常采用基 于像元的分类方法。
【发明内容】
[0008] 有鉴于此，本发明的目的在于提供一种集成像元与对象的水稻自动变化检测方 法，适用于水稻变化连续监测，具有自动化程度高、简单易用、鲁棒性好以及分类精度高等 特点。
[0009] 为实现上述目的，本发明采用如下技术方案:一种集成像元与对象的水稻自动变 化检测方法，其特征在于，包括以下步骤： 步骤S01:构建研究区的起始年份与结束年份的M0DIS水体指数时序数据集，然后逐像 元分别建立起始年份与结束年份的M0DIS水体指数时序曲线； 步骤S02:逐像元计算起始年份与结束年份的M0DIS水体指数时序曲线的相似度； 步骤S03:基于所述起始年份与结束年份的M0DIS水体指数时序曲线的相似度，判断疑 似变化区域； 步骤S04:综合起始年份的多期Landsat影像，通过影像分割获得影像对象； 步骤S05:将所述疑似变化区域与所述影像对象进行叠加，获得研究区疑似变化对象； 步骤S06:综合结束年份的多期Landsat影像，将所述研究区疑似变化对象进一步细分 割； 步骤S07:对细分割后的研究区疑似变化对象分别进行水稻制图，获得起始年份与结束 年份的研究区疑似变化对象的水稻分布图； 步骤S08:获得研究区起始年份至结束年份的水稻动态变化分布图。
[0010] 本发明与现有技术相比具有以下有益效果： 1、本发明充分利用低分辨率遥感影像时序数据比较完整，因而能很好地揭示研究区变 化区域的特点，基于低分辨率遥感影像时序数据提取研究区疑似变化区域，能有效地缩小 变化监测范围，避免多次分类带来的累加误差。
[0011] 2、本发明采用动态时间弯曲距离，计算起始和结束年份水体指数时序曲线的相似 度，通过水体指数时序曲线的相似性程度，获得研究区疑似变化区域。水稻与非水稻的水体 指数时序曲线存在本质上的差异，并且动态时间弯曲距离能很好地度量时序曲线的相似 性，因此采用本发明获得的研究区疑似变化区域具有科学性与合理性。
[0012] 3、本发明通过基于起始年份遥感影像分割获得影像对象的前提下，将基于像元检 测得到的研究区疑似变化区域与之叠加，从而获得研究区疑似变化对象，进一步在疑似变 化对象中综合结束年份的遥感影像特征做细分割，有效地避免了分割对象与实际变化区域 的不匹配，真正实现了基于像元与对象的结合。
【附图说明】
[0013] 图1是本发明本发明的实现流程图。
[0014] 图2是本发明一实施例的研究区起始年份与结束年份水稻指数的DTW距离的空间 分布图。
[0015] 图3是本发明一实施例的研究区起始年份与结束年份水体指数的DTW距离的直方 图。
[0016] 图4是本发明一实施例的疑似变化区域空间分布图。
[0017] 图5是本发明一实施例的分割后形成的影像对象分布图。
[0018] 图6A是本发明一实施例的基于像元的水体指数与增强型植被指数变化比值指数 的空间分布图。
[0019] 图6B是本发明一实施例的基于像元的纹理特征空间分布图。
[0020]图7是本发明一实施例的研究区2014-2015年水稻变化监测结果图。
【具体实施方式】
[0021]下面结合附图及实施例对本发明做进一步说明。
[0022]请参照图1，本发明提供一种集成像元与对象的水稻自动变化检测方法，其特征在 于，包括以下步骤： 步骤S01:构建研究区的起始年份与结束年份的M0DIS水体指数时序数据集，然后逐像 元分别建立起始年份与结束年份的M0DIS水体指数时序曲线； 利用8天最大化合成的M0DIS波段反射率数据，计算水体指数LSWK1年共有46期8天最 大化合成的M0DIS数据。逐期计算基于M0DIS水体指数，获得年内MODIS LSWI时序数据。 [0023]基于年内MODIS LSWI时序数据，在排除有云干扰的观测值的情况下，利用线性插 值方法，逐像元建立研究区逐日LSWI时序数据。然后采用Whittaker smoother等数据平滑 方法，逐像元构建研究区年内逐日MODIS LSWI时序数据集。依据此方法，分别建立起始年份 与结束年份的逐日MODIS LSWI时序数据集。在此基础上，逐像元分别建立开始年份与结束 年份的M0DIS水体指数时序曲线。
[0024]步骤S02:采用动态时间弯曲距离逐像元计算起始年份与结束年份的M0DIS水体指 数时序曲线的相似度，动态时间弯曲距离(Dynami c time warping，简称DTW)采用动态规划 的思想，通过构建邻接矩阵，寻找起始与结束年份水体指数时序曲线之间的最短距离;具体 包括以下步骤： 步骤S21:设起始年份与结束年份的水体指数时序曲线，分别为长度是1*365的两个时 间序列LSWI_A{al，a2，a3, · · ·，a365}和LSWI_B{bl，b2，b3, · · ·，b365}; 步骤S22:计算1^11_4和1^11_8中每个元素之间的距离，得到一个365*365的距离矩阵 M;使用欧式距离计算每个元素之间的距离，即:M(i，j) = (ai_bj)2,其中ai，bj分别属于起始 年份和结束年份水体指数时间序列1^11_4和1^11_8, ; 步骤S23:查找距离矩阵Μ的最短路径D{d(l)，d(2)，. . .，d(x)，. . .，d(k)}，其中， M5HO29, Sfc ;在最短路径 D{d(l)，d(2)，· · ·，d(x)，· · ·，d(k)}中，当中 x=l 时，d (l)=M(l，l)，iX>Ut，gd(X)=M(i，j)dlJd(X+l)=min{M(i+l，j)，M(i，j+l)，M(i+l，j+l)}; 步骤S24:累加最短路径D{ d(l)，d(2)，...，d(x)，...，d(k)}，最终获得结束年份LSWI_ B和起始年份LSWI_A水体指数时序曲线的动态时间弯曲距离DTW;
其中：d(x)为距离矩阵Μ的最短路径D{d(l)，d(2)，. . .，d(x)，. . .，d(k)}中第X个数的 取值。
[0025] 采用动态时间弯曲距离（Dynamic Time Warping distance，简称DTW距离），计算 起始年份与结束年份水体指数的相似度。在所有农作物中，水稻与非水稻最大的区别，在于 水稻在整个生长期内土壤都比较湿润。尤其是水稻在移栽期前，通常需要漫灌，秧苗移栽后 很长一段时间内稻田处于淹水状态。水体指数时序曲线，能很好地表征水稻生长期内地表 水分的动态变化情况。水稻与非水稻的水体指数时序曲线存在非常大的差距。因此，通过评 估两个年份水体指数的相似性，能有效地监测水稻的动态变化情况。在本实施例中，基于 MODIS LSWI时序数据，获得的研究区2014-2015年水体指数时序曲线的动态时间弯曲距离 的空间分布图见图2。
[0026]步骤S03:基于所述起始年份与结束年份的M0DIS水体指数时序曲线的相似度，判 断疑似变化区域;其中将起始年份与结束年份的MODI S水体指数时序曲线的动态时间弯曲 距离超过阈值的像元判断为疑似变化像元，所有的疑似变化像元的汇总即为所述疑似变化 区域。
[0027]某像元起始与结束年份MODI S水体指数时序曲线的相似性越大，这两个年份内发 生变化的可能性越小。反之，如果某像元起始与结束年份MODI S水体指数时序曲线的相似性 很小，则极有可能已经发生了变化。起始与结束年份MODI S水体指数时序曲线的动态时间弯 曲距离，能很好地指示这两个年份的水体指数时序曲线的相似性。因此，可以基于起始与结 束年份M0DIS水体指数时序曲线的动态时间弯曲距离，判断研究区疑似变化区域。当该像元 起始与结束年份M0DIS水体指数时序曲线的动态时间弯曲距离超过一定的阈值时，判断为 研究区疑似变化像元。
[0028]在本实施例中，从研究区2014-2015年水体指数时序曲线的动态时间弯曲距离的 直方图（图3)可以看出，其基本符合正态分布特征。依据显著性检验，选取正态分布中右侧 单尾95%置信度区间的值，作为设定的阈值，作为该像元的M0DIS水体指数时序曲线的动态 时间弯曲距离是否超出正常的依据。如本实施例中确定的阈值为0.14。逐像元依据上述步 骤进行判断，汇总研究区所有的疑似变化像元，获得的研究区疑似变化区域分布图（图4)。 [0029]步骤S04:综合起始年份最能体现水稻生长特征的多期Landsat影像主要波段的第 一主成分数据，进行影像分割获得影像对象。
[0030] 综合起始年份最能体现水稻生长特征的多期遥感影像，如水稻移栽期、分蘖期、抽 穗期以及成熟期等多个时期的Landsat遥感影像，作为遥感影像分割的数据基础。可以将这 些时期Landsat遥感影像的一些重要波段数据，如近红外、短波红外波段，通过主成分变换 然后选取第一主成分数据进行遥感影像分割。利用Canny、Rob ert等分割算子，将研究区 Landsat影像进行分割，获得Landsat影像对象。如综合上述步骤获得的研究区2014年 Landsat影像对象图见图5。
[0031] 步骤S05:将所述疑似变化区域与所述影像对象进行叠加，获得研究区疑似变化对 象;所述研究区疑似变化对象的判断依据为:将所述疑似变化区域与所述影像对象进行叠 加，只要影像对象中存在步骤S03中获得的疑似变化区域，则判断该影像对象为研究区疑似 变化对象。
[0032]步骤S06:综合结束年份最能体现水稻生长特征的多期Landsat影像主要波段的第 一主成分数据，对步骤S05获得的研究区疑似变化对象进一步细分割，获得综合起始年份与 结束年份的多期Landsat影像细分割后的研究区疑似变化对象。
[0033] -般而言，发生变化的区域通常比较小。并且，通常存在起始年份的影像对象部分 发生改变的情况。因此，有必要综合结束年份的Landsat影像，对研究区疑似变化对象做进 一步的细分割。
[0034] 综合结束年份最能体现水稻生长特征的多期遥感影像，如水稻移栽期、分蘖期、抽 穗期以及成熟期等多个时期的Landsat遥感影像，作为研究区疑似变化对象进一步细分割 的数据基础。
[0035] 可以选取结束年份多期Landsat影像的一些重要波段数据(如近红外、短波红外波 段），通过主成分变换，选取第一主成分数据开展研究区疑似变化对象细分割。利用Canny、 Robert等分割算子，将步骤S05获得的研究区疑似变化对象进一步细分割，获得综合起始与 结束年份细分割后的研究区疑似变化对象。
[0036]步骤S07:对细分割后的研究区疑似变化对象分别进行水稻制图，获得起始年份与 结束年份的研究区疑似变化对象的水稻分布图；综合多期Landsat影像，建立基于像元的 Landsat影像特征，所述Landsat影像特征包括分蘖期影像纹理，水稻分蘖期到抽穗期内水 体指数的最小值、水体指数与增强型植被指数变化比值指数;然后采用中位数方法汇总到 Landsat影像对象，建立面向对象的Landsat影像特征;基于面向对象的Landsat影像特征， 针对所述步骤S06细分割后生成的研究区疑似变化对象，分别进行水稻制图，获得起始年份 与结束年份的研究区疑似变化对象的水稻分布图，具体如下： 基于结束年份多期Landsat遥感影像，从分蘖期影像纹理、水稻分蘖期到抽穗期内水体 指数的最小值、水体指数与增强型植被指数变化比值指数等方面，建立基于像元的Landsat 影像特征。如计算获得的水稻分蘖期到抽穗期水体指数与增强型植被指数变化比值指数的 空间分布图见图6A。如选取水稻分蘖期Landsat遥感影像的第一主成分数据，通过三层离散 小波分解，母小波基为db4,邻域窗口为7X7,所获得的研究区纹理指标的空间分布图见图 6B〇
[0037]将基于像元的Land sat遥感影像特征，通过取中位数的方法，汇总到研究区 Landsat影像对象，从而建立面向对象的Landsat遥感影像特征。考虑到如果影像对象中个 别像元的数据质量不理想，将会对该像元的Landsat遥感影像特征带来干扰。因此，选取中 位数而非平均值的方法，建立面向对象的Landsat遥感影像特征。针对步骤S06细分割后生 成的研究区疑似变化对象，基于分蘖期影像纹理、水稻分蘖期到抽穗期内水体指数的最小 值、水体指数与增强型植被指数变化比值指数等遥感影像特征，进行结束年份的水稻制图， 获得结束年份研究区疑似变化对象的水稻分布图。
[0038] 水稻的判断规则为:分蘖期影像纹理大于1800,水稻分蘖期到抽穗期水体指数的 最小值大于〇. 1，并且水体指数与增强型植被指数变化比值指数小于〇. 6。
[0039] 按照上述步骤，开展起始年份研究区疑似变化对象的水稻制图。首先基于起始年 份多期Landsat遥感影像，建立基于像元的Landsat影像特征，然后采用中位数方法汇总到 Landsat影像对象，进而建立研究区起始年份面向对象的Landsat遥感影像特征。基于分蘖 期影像纹理、水稻分蘖期到抽穗期内水体指数的最小值、水体指数与增强型植被指数变化 比值指数等遥感影像特征，对步骤S06细分割后生成的研究区疑似变化对象，开展起始年份 的水稻制图，获得起始年份研究区疑似变化对象的水稻分布图。
[0040]步骤S08:获得研究区起始年份至结束年份的水稻动态变化分布图。将所述步骤 S07获得的起始年份与结束年份的研究区疑似变化对象的水稻分布图进行叠加，以起始年 份为非水稻而结束年份为水稻的区域为水稻增加区，以起始年份为水稻而结束年份为非水 稻的区域为水稻减少区，获得研究区起始年份至结束年份的水稻动态变化分布图，综上述 步骤，本实施例中研究区哈尔滨市及其附近区域2014-2015年水稻变化区域分布图见图7。 [0041]以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与 修饰，皆应属本发明的涵盖范围。
【主权项】
1. 一种集成像元与对象的水稻自动变化检测方法，其特征在于，包括W下步骤： 步骤SOI:构建研究区的起始年份与结束年份的MODIS水体指数时序数据集，然后逐像 元分别建立起始年份与结束年份的M0DIS水体指数时序曲线； 步骤S02:逐像元计算起始年份与结束年份的M0DIS水体指数时序曲线的相似度； 步骤S03:基于所述起始年份与结束年份的M0DIS水体指数时序曲线的相似度，判断疑 似变化区域； 步骤S04:综合起始年份的多期Landsat影像，通过影像分割获得影像对象； 步骤S05:将所述疑似变化区域与所述影像对象进行叠加，获得研究区疑似变化对象； 步骤S06:综合结束年份的多期Landsat影像，将所述研究区疑似变化对象进一步细分 割； 步骤S07:对细分割后的研究区疑似变化对象分别进行水稻制图，获得起始年份与结束 年份的研究区疑似变化对象的水稻分布图； 步骤S08:获得研究区起始年份至结束年份的水稻动态变化分布图。2. 根据权利要求1所述的集成像元与对象的水稻自动变化检测方法，其特征在于:所述 步骤S02中，采用动态时间弯曲距离计算起始年份与结束年份的MODIS水体指数时序曲线的 相似度，具体包括W下步骤： 步骤S21:设起始年份与结束年份的水体指数时序曲线，分别为长度是1*365的两个时 间序列LSWI_A{al，a2，a3,...，曰365}和LSWI_B{bl，b2，b3,. . .，b365}; 步骤S22:计算LSWI_A和LSWI_B中每个元素之间的距离，得到一个365*365的距离矩阵 M;使用欧式距离计算每个元素之间的距离，即：1(1^)=(曰1-叫）*(曰1-叫），其中曰1，叫分别 属于起始年份和结束年份水体指数时间序列LSWI_A和LSWI_B，i.<i <365，！< |<說5 ; 步骤S23:查找距离矩阵Μ的最短路径D{d(l)，d(2)，. . .，d(x)，. . .，d(k)}，其中， 3财 <k < 巧!^ J k ;在最短路径D{d(l)，d(2)，...，d(x)，...，d(k)}中，当中x=l时，d (1)=M(1，1)，当x〉l时，若d(x)=M(i,j)，则d(x+l)=min{M(i+l，j)，M(i,j+1)，M(i+l，j+1)}; 步骤S24:累加最短路径D{ d(l)，d(2)，...，d(x)，...，d(k)}，最终获得结束年份LSWI_ B和起始年份LSWI_A水体指数时序曲线的动态时间弯曲距离DTW;其中：d(x)为距离矩阵Μ的最短路径D{d(l)，d(2)，. . .，d(x)，. . .，cKk)}中第X个数的 取值。3. 根据权利要求1所述的集成像元与对象的水稻自动变化检测方法，其特征在于:所述 步骤S03中，将起始年份与结束年份的M0DIS水体指数时序曲线的动态时间弯曲距离超过阔 值的像元判断为疑似变化像元，所有的疑似变化像元的汇总即为所述疑似变化区域。4. 根据权利要求1所述的集成像元与对象的水稻自动变化检测方法，其特征在于:所述 步骤S04中，综合起始年份最能体现水稻生长特征的多期Landsat影像主要波段的第一主成 分数据，进行影像分割获得影像对象。5. 根据权利要求1所述的集成像元与对象的水稻自动变化检测方法，其特征在于:所述 步骤S05中，所述研究区疑似变化对象的判断依据为:将所述疑似变化区域与所述影像对象 进行叠加，只要影像对象中存在步骤S03中获得的疑似变化区域，则判断该影像对象为研究 区疑似变化对象。6. 根据权利要求1所述的集成像元与对象的水稻自动变化检测方法，其特征在于:所述 步骤S06中，综合结束年份最能体现水稻生长特征的多期Landsat影像主要波段的第一主成 分数据，对步骤S05获得的研究区疑似变化对象进一步细分割，获得综合起始年份与结束年 份的多期Landsat影像细分割后的研究区疑似变化对象。7. 根据权利要求1所述的集成像元与对象的水稻自动变化检测方法，其特征在于:所述 步骤S07中，综合多期Landsat影像，建立基于像元的Landsat影像特征，所述Landsat影像特 征包括分葉期影像纹理，水稻分葉期到抽穗期内水体指数的最小值、水体指数与增强型植 被指数变化比值指数；然后采用中位数方法汇总到Landsat影像对象，建立面向对象的 Landsat影像特征;基于面向对象的Landsat影像特征，针对所述步骤S06细分割后生成的研 究区疑似变化对象，分别进行水稻制图，获得起始年份与结束年份的研究区疑似变化对象 的水稻分布图。8. 根据权利要求1所述的集成像元与对象的水稻自动变化检测方法，其特征在于:所述 步骤S08中，将所述步骤S07获得的起始年份与结束年份的研究区疑似变化对象的水稻分布 图进行叠加，W起始年份为非水稻而结束年份为水稻的区域为水稻增加区，W起始年份为 水稻而结束年份为非水稻的区域为水稻减少区，获得研究区起始年份至结束年份的水稻动 态变化分布图。
【文档编号】G06K9/62GK106096630SQ201610385818
【公开日】2016年11月9日
【申请日】2016年6月3日 公开号201610385818.9, CN 106096630 A, CN 106096630A, CN 201610385818, CN-A-106096630, CN106096630 A, CN106096630A, CN201610385818, CN201610385818.9
【发明人】邱炳文, 卢迪非, 齐文, 刘哲
【申请人】福州大学