基于时序相似性的植被变化发生时间检测方法

摘要

本发明涉及一种基于时序相似性的植被变化发生时间检测方法，该方法首先建立研究区多年时空连续的植被指数时序数据，然后逐像元逐年计算历年与起始年份植被指数时序曲线的JM距离，生成历年与起始年份JM距离的时序曲线；利用logistic模型拟合历年与起始年份JM距离的时序曲线，从logistic模型参数中获取时间参数，实现植被变化时间的自动提取。该方法利用历年与起始年份植被指数时序曲线的JM距离指示时序相似性，从逐年时序相似性的变化规律中获取植被发生变化时间。该方法能有效地检测时序曲线在幅度、频率等各方面的变化，避免了将原始光谱指数时序数据进行分解的繁琐程序，解决了难以直接从原始光谱指数时序数据中提取指标来全面表征植被变化的难题。

说明书

技术领域

本发明数据挖掘技术领域，特别是涉及一种基于时序相似性的植被变化发生时间检测方法。

背景技术

植被为我们提供氧气和食物，在生态系统平衡中发挥着非常重要的作用。植被变化与全球变化关系密切，因此备受关注。目前植被变化监测集中在森林植被变化方面。在植被遥感动态监测方面，比较常用的方法有LandTrendr(Landsat-based Detection of Trends in Disturbance and Recovery)和BFAST(Break Detection For Additive and Trend)方法。这些基于遥感影像时间序列数据的变化监测方法，为植被变化时空连续监测提供了新的发展导向。但这些方法一般基于光谱指数时序数据，通常需要开展时间序列的分解和重构，通过阈值的设置判断植被发生突变或干扰情况。由于遥感影像原始波段反射率数据往往受到大气条件、太阳高度角变化等各种因素的影响，从而导致在此基础上计算的光谱指数数据不可避免存在一定的不确定性。因此，在基于光谱指数的植被变化监测方法应用过程中，不可避免存在一定的问题。本发明拟从历年与起始年份时序相似性变化的角度，通过建立时序相似性变化曲线，揭示历年与起始年份相比是否发生了变化，达到自动获取植被变化发生时间的目的。

发明内容

有鉴于此，本发明的目的是提供一种基于时序相似性的植被变化发生时间检测方法，该方法适用于大范围快速监测的需求，具有自动化程度高、简单易用、鲁棒性好以及分类精度高等特点。

本发明采用以下方案实现：一种基于时序相似性的植被变化发生时间检测方法，包括以下步骤：

步骤S01：逐像元建立多年时空连续植被指数时序曲线；

步骤S02：基于植被指数时序曲线，逐年依次计算其他年份与起始年份的JM距离；

步骤S03：依次按时间顺序生成其他各年份与起始年份JM距离的时序曲线；

步骤S04：基于多年JM距离时序曲线，进行logistic模型拟合；

步骤S05：从logistic模型拟合结果中获取植被变化发生时间；

步骤S06：对植被变化发生时间进行自动提取，获得研究区植被变化时间分布图。

特别地，该方法逐像元基于多年植被指数时序曲线，通过逐年计算与开始年份的Jeffries-Matusita(JM)距离，衡量历年与开始年份的时序相似性，进一步基于多年时序相似性曲线的变化规律，实现植被变化发生时间的自动提取。

进一步地，在所述步骤S02中，逐期基于遥感影像近红外、红光波段反射率数据，计算植被指数；按照时间顺序，生成原始植被指数时序数据；然后采用Whittaker smoother数据平滑方法，逐像元构建多年时空连续的植被指数时序数据集，在此基础上逐像元建立多年时空连续植被指数时序曲线。

进一步地，在所述步骤S02中，逐年依次计算该年份与起始年份植被植被指数时序曲线的JM距离；通过两个年份植被指数时序曲线的JM距离揭示其时序相似性变化；时序曲线的JM距离能很好地揭示不同时序曲线在频率、幅度以及物候等各方面的差异；JM距离越大表示两者的相似程度越小，反之，JM距离越小表示两个年份的时序相似性越强。

进一步地，在所述步骤S03中，基于历年与起始年份植被植被指数时序曲线的JM距离，依次按时间顺序生成JM距离时序曲线，用以指示历年与起始年份的时序相似性变化规律。

进一步地，在所述步骤S04中，对历年与起始年份的JM距离时序曲线，进行logistic模型拟合；从logistic模型拟合参数中获取植被变化类型以及变化时间。

进一步地，在所述步骤S04中，logistic模型的公式如下所示：

其中：f(x)表示其他年份与起始年份的JM距离，自变量x为时间，用年份表示；其中参数a代表了研究时段内JM距离的变化量；参数b代表变化速率；参数c指示变化发生的时间；参数d表示变化发生前的JM距离。

进一步地，logistic模型参数b代表变化速率同时也指示变化类型；其中变化速率接近1表示为渐变型；设变化速率b的值域处在[0.9,1.1]区段内，则像元的变化类型为渐变型；变化速率b大于1.1或者小于0.9，则像元的变化类型为突变型。

进一步地，在所述步骤S05中，对于突变型的像元，从logistic模型参数c中获取植被变化发生时间。

特别地，多年植被指数时序数据/时序曲线为从起始年份的元旦开始，按时间顺序，一直到结束年份的年末，按照一定的时间步长，如每8天或逐日记录记录所形成的植被指数的数据序列/时序曲线。

进一步地，在所述步骤S01-S06中，通过计算历年与起始年份时序曲线的时序相似性的变化规律，获得植被变化发生时间。

特别地，该方法在城市化、退耕还林、耕地抛荒、田园园林化、林地采伐与矿山开采引起植被变化的发生时间自动检测领域中的应用。

与现有技术相比，本发明的显著优点在于：

1、基于时序相似度，而非原始光谱指数时序数据，一方面避免了将原始光谱指数时序数据分解为趋势、季节和扰动项等繁琐程序，另一方面也解决了难以直接从原始光谱指数时序数据中提取指标来全面表征植被变化的难题。

2、可以不借助已知训练数据，不需要人机交互，不依赖监督学习或机器学习方法，简便地实现植被变化时间的自动获取。

附图说明

图1为本发明实施例的实现流程图。

图2为2001-2015年MODIS OSAVI时序曲线。

图3为2002-2015年与2001年逐年JM距离的时序曲线图。

图4为logistic模型及其模型参数对应的含义图。

图5为研究区植被变化发生时间的空间分布图。

图6为2002-2015年历年植被发生变化的面积的直方图。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

请参照图1，本实施例提供一种基于时序相似性的植被变化发生时间检测方法，包括以下步骤：

步骤S01：建立2001-2015年土壤调节型植被指数MODIS OSAVI时序曲线。

利用MODIS波段反射率数据，计算MODIS OSAVI时序数据，计算公式为：

其中NIR,Red分别为MODIS的近红外、红光波段的反射率。依据上述公式，逐期基于遥感影像波段数据，计算植被指数。按照时间顺序，生成原始MODIS OSAVI时序数据。数据的时间步长为8天。然后采用Whittaker smoother等数据平滑方法，逐像元构建多年时空连续的MODIS OSAVI时序数据集。在此基础上，逐像元建立2001-2015年土壤调节型植被指数MODIS OSAVI时序曲线。以宁夏回族自治区灵武市某区域实施退耕还林为例，所形成的2001-2015年8天合成的MODIS OSAVI时序曲线如图2所示。

其中，MODIS数据为中分辨率成像光谱仪数据，全称为Moderate Resolution Imaging Spectroradiometer。OSAVI为土壤调节植被指数，全称为Optimized Soil Adjusted Vegetation Index。用于表征植被生长状态以及空间分布密度。

步骤S02：基于MODIS OSAVI时序曲线，逐年计算其他年份与起始年份的JM距离

逐年计算2002-2015年与起始年份2001年MODIS OSAVI时序曲线的Jeffries-Matusita(JM)距离：

其中：(p(x|w_i))^1/2为条件概率密度。JM_i,j的值在0～2之间，其大小指示时序曲线之间的相似程度。JM距离的数值越大表示两条时序曲线的时序相似性越小。比如当0<JM_i,j<1.0时，两条时序曲线比较相似；1.0<JM_i,j<1.5时，两条时序曲线具有一定的时序相似性；1.5<JM_i,j<2.0时，两条时序曲线的时序相似性比较小。

步骤S03：生成2002-2015年与起始年份2001年的JM距离的时序曲线图；

基于逐年计算的该年份与起始年份2001年MODIS OSAVI时序曲线的JM距离，生成2002-2015年与起始年份2001年的JM距离的时序曲线。该时序曲线指示2002-2015年与2001年植被指数的时序相似性变化规律。以图2为例，生成的2002-2015年与2001年逐年的JM距离的时序曲线图见图3。从2002到2007年这段时间内，各年份与起始年份2001年的JM距离都非常小，在0.6以下。从2008年开始略有增大，但仍在1.0以下。2009年与2001年的JM距离迅速增大到1.4以上。从2010年开始，各年份与起始年份2001年的JM距离均在1.7以上，接近2.0。2002-2015年与起始年份2001年的JM距离的时序曲线图，能很好地指示历年与起始年份2001年的植被指数时序曲线的时序相似度的变化规律。从2002到2015年，历年与起始年份2001年的时序相似度经历了一个从比较强突变到比较弱的过程。

步骤S04：基于2002-2015年JM距离时序曲线，进行logistic模型拟合

logistic模型的公式如下所示：

其中：f表示其他年份与2001年的JM距离，自变量x为时间，用年份表示。logistic模型中有四个参数，a,b,c,d分别具有一定的指示意义。其中参数a代表了研究时段内的变化量，变化量越大表示变化的程度越高；b代表变化速率；c指示变化发生的时间；d表示变化发生前的JM距离。logistic模型及其模型参数对应的含义图见图4。

步骤S05：从logistic模型拟合结果中获取植被变化发生时间；

logistic模型拟合的四个参数中，参数b代表变化速率同时也指示了变化类型(突变或渐变)。其中变化速率接近1表示为渐变型，本实施例中，如果变化速率的值域处在[0.9,1.1]区段内，则设定该像元的变化类型为渐变型。如果变化速率大于1.1或者小于0.9，则设定该像元的变化类型为突变型。对于突变型的像元，进一步从logistic模型的参数c中获取植被变化发生时间。对于渐变型的像元，从logistic模型的参数c所获得到的变化发生时间一般落在研究时段范围之外。在本实施例中不予考虑。

步骤S06：实现植被变化时间自动提取，获得研究区植被变化时间分布图；

基于上述所建立的植被变化发生时间检测流程与方法，逐像元提取植被变化发生时间，最终生成研究区植被变化发生时间分布图。依据上述流程，可实现植被变化发生时间的快速自动提取。以中国三北防护林地区为例，获得研究区植被变化发生时间的空间分布图见图5。2002-2015年历年植被发生变化的面积直方图见图6。由图可见，2006-2007年这两年内较大面积的植被发生了变化。

以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。