本发明涉及遥感信息处理领域，具体涉及一种基于多时序指标变化趋势的退耕还林自动识别方法。
背景技术：
万物土中生。耕地是人类赖以生存和发展最基本的自然资源。随着人口的增长，人类向土地不断索取食物的需求加剧。山地丘陵区由于耕地稀少，盲目开垦播种坡地带来严重的水土流失问题。干旱区沙地耕种，造成干旱、沙尘暴等自然灾害频发，严重威胁了生态安全。因此通过投入大量资金，有计划、分步骤地实施退耕还林政策。
据报道，自退耕还林政策实施10多年来，取得了显著的成效。但退耕还林究竟发生在哪些区域，有必要开展退耕还林实施区域大范围实时监测，而且大范围快速自动监测退耕还林的空间分布，具有重要意义。传统的人工调查方法不仅费时费力，而且难以实现全覆盖，容易受到主观因素的干扰。遥感影像具有时效性强、覆盖范围大的特点，在土地变化检测中发挥着重要作用。然而，由于不同土地覆盖光谱的“类内异质性、类间相似性”的干扰，以及遥感分类人员的经验、分类器、相关数据的丰富程度以及研究区特点等众多因素的影响，导致土地覆盖自动制图精度难以保证。
传统的土地覆盖变化检测多基于分类后变化检测的方法，这种方法的不足之处在于：多次分类带来累计误差，并且发生土地覆盖变化区域本身比较小，很容易造成错分或遗漏现象。因此直接提取变化区域的土地覆盖变化检测方法受到极大关注。本发明所提出的方法，通过设计多个时序指标，刻画耕地范围内种植的单季农作物、多季农作物以及森林等多种植被类型，进而依据多个时序指标的变化趋势，揭示单季农作物或多季农作物转变为森林的变化过程，最终建立一种基于多时序指标变化趋势的退耕还林自动识别方法。
技术实现要素：
本发明的目的在于提供一种基于多时序指标变化趋势的退耕还林自动识别方法，以克服现有技术中存在的缺陷。
为实现上述目的，本发明的技术方案是：一种基于多时序指标变化趋势的退耕还林自动识别方法，按照如下步骤实现：
步骤s1：逐像元建立研究区历年增强型植被指数年内逐日时序数据集；
步骤s2：设定丰度指标、时序离散度指标以及生长期长度指标；
步骤s3：设定高分位持续度指标、低分位持续度指标；
步骤s4：建立上述五个时序指标的多年时序数据集；
步骤s5：计算所述五个时序指标的多年变化趋势q；
步骤s6：对所述五个时序指标的多年变化趋势q的显著性进行检验；
步骤s7：建立退耕还林自动识别流程图；
步骤s8：获得研究区退耕还林空间分布图。
在本发明一实施例中，在所述步骤s1中，所述增强型植被指数年内逐日时序数据集按照如下步骤获取：
步骤s11：通过基于每日的modis遥感影像波段数据逐日计算或通过基于8天最大化合成的modis遥感影像波段数据计算以及线性插值获取原始年内时序数据集；
步骤s12：利用whittakersmoother平滑方法，对所述原始年内时序数据集进行平滑处理，获取所述增强型植被指数年内逐日时序数据集。
在本发明一实施例中，在所述步骤s2中，计算所述增强型植被指数年内逐日时序数据集中增强型植被指数年内逐日时序数据的第50分位值，并提取其中大于或等于第50分位值的所有数据，并计算该所有数据的平均值，作为所述丰度指标。
在本发明一实施例中，在所述步骤s2中，计算所述增强型植被指数年内逐日时序数据集中增强型植被指数年内逐日时序数据的第75分位值，逐像元依次提取所述增强型植被指数年内逐日时序数据集中大于或等于第75分位值的所有数据，建立高值区域；将所述高值区域的极差和标准差的乘积，作为所述时序离散度指标。
在本发明一实施例中，在所述步骤s2中，计算所述增强型植被指数年内逐日时序数据集中增强型植被指数年内逐日时序数据的最大值max和最小值min，计算该最大值max和最小值min的极差(max-min)，将按时间顺序依次获取的大于或等于(min+(max-min)/2)的数据对应的持续时间累加，作为生长期长度指标。
在本发明一实施例中，在所述步骤s3中，计算所述增强型植被指数年内逐日时序数据集中增强型植被指数年内逐日时序数据的第65分位值，将经逐像元按时间顺序依次检索大于或等于第65分位值的持续时间后获取的最大值作为高分位持续度指标。
在本发明一实施例中，在所述步骤s3中，计算所述增强型植被指数年内逐日时序数据集中增强型植被指数年内逐日时序数据的第50分位值，将经逐像元按时间顺序依次检索大于或等于第50分位值的持续时间后获取的次大值作为低分位持续度指标。
在本发明一实施例中，在所述步骤s5中，通过采用sen氏斜率法获取所述多年变化趋势q。
在本发明一实施例中，在所述步骤s6中，所述多年变化趋势q的显著性检验包括空间显著性检验以及时序显著性检验且将显著性检验的结果记为：显著正趋势、无趋势和显著负趋势。
在本发明一实施例中，在所述步骤s7中，通过如下方式建立所述退耕还林自动识别流程图：
(1)当丰度和生长期长度均具有显著正趋势，而时序离散度呈现出显著负趋势，并且高分位持续度、低分位持续度均无趋势时，则该像元由单季农作物转变为森林，即该像元为退耕还林区域；
(2)当高分位持续度、生长期长度均具有显著正趋势，而低分位持续度呈现出显著负趋势，并且丰度、时序离散度均无趋势时，则该像元由多季农作物转变为森林，即该像元为退耕还林区域。
相较于现有技术，本发明具有以下有益效果：
(1)通过采用基于植被指数时序数据分位数，划分若干植被生长阶段。由于单季农作物、多季农作物、森林等植被在不同植被生长阶段呈现出各自的特点。基于若干植被生长阶段，通过从覆盖时间、平均状态、变化幅度等侧面，设计丰度、时序离散度、连续度等多个时序指标，用于刻画不同植被类型的多方位特征，具有很好的灵活性和拓展性。
(2)通过采用基于多个时序指标的总体变化趋势，对某些年份数据出现错误或异常等情况，具有更好的容错能力。
(3)不依赖多次分类，直接自动识别退耕还林这种特定的变化方式，具有简洁明了的特点。
(4)可以不借助已知训练数据，不需要人机交互，不依赖监督分类或机器学习方法，实现退耕还林自动识别。
附图说明
图1为本发明一实施例中的流程图。
图2为本发明一实施例中单季农作物2001-2016年退耕还林点位的evi时序曲线图。
图3为本发明一实施例中双季农作物2001-2016年退耕还林点位的evi时序曲线图。
图4为本发明一实施例中单季农作物2001-2016年退耕还林点位的五个时序指标的变化趋势图。
图5为本发明一实施例中双季农作物2001-2016年退耕还林点位的五个时序指标的变化趋势图。
图6为本发明一实施例中2001-2016年研究区五个时序指标变化趋势的空间分布图。
图7为本发明一实施例中退耕还林自动识别流程图。
图8为本发明一实施例中研究区退耕还林空间分布图。
具体实施方式
下面结合附图，对本发明的技术方案进行具体说明。
本发明提供一种基于多时序指标变化趋势的退耕还林自动识别方法，如图1所示：首先建立研究区每个栅格单元的多年逐日增强型植被指数时序数据集，然后逐像元基于增强型植被指数时序数据，逐年提取丰度、时序离散度、生长期长度、高分位持续度、低分位持续度五个时序指标，逐像元依次检测上述五个时序指标的多年变化趋势，在此基础上基于五个时序指标的变化趋势，建立退耕还林自动识别流程图，最终达到退耕还林自动识别的目标。
为了让本领域技术人员进一步了解本发明提出的基于多时序指标变化趋势的退耕还林自动识别方法，下面结合具体步骤以及实施例进行说明。
在本实施例中，具体包括以下步骤：
步骤s01：逐像元建立研究区多年逐日植被指数时序数据。
在本实施例中，分别建立研究区年内逐日evi时序数据集。modisevi时序数据的建立，可以基于每日的modis遥感影像波段数据逐日计算，也可以基于8天最大化合成的modis遥感影像波段数据进行计算的基础上通过线性插值获得。利用whittakersmoother平滑方法，对原始年内时序数据集进行平滑处理，从而获得研究区平滑的年内逐日evi时序数据集，作为退耕还林自动识别的基础。如基于年内逐日evi时序数据集，建立的研究区退耕还林的evi时序曲线图见图2和图3。
在本实施例中，modis数据为中分辨率成像光谱仪数据，全称为moderateresolutionimagingspectroradiometer。
植被指数是表征植被生长状态以及空间分布密度的因子。常见的植被指数有ndvi和evi。ndvi为归一化植被指数，全称为normalizeddifferencevegetationindex。evi为增强型植被指数，全称为enhancedvegetationindex。evi指数的计算公式为：
步骤s02：设计丰度、时序离散度以及生长期长度指标
基于年内逐日modisevi时序数据，设定体现植被生长状态的若干时序指标，包括：丰度、时序离散度以及生长期长度。首先提取modisevi年内时序数据的极值以及重要分位值，如最大值、最小值、第50分位值、第65分位值、第75分位值。然后依据这些极值和分位值将modisevi年内时序数据划分为体现不同植被生长状态的若干区间，在此基础上综合时间和植被指数值域分布情况两方面设计建立丰度、时序离散度以及生长期长度指标。
进一步的，基于modisevi年内时序数据的第50分位值，建立丰度指标，其计算流程如下：提取增强型植被指数大于或等于第50分位值的所有数据，求取其平均值，作为丰度的数值。
进一步的，基于modisevi年内时序数据的第75分位值，提取时序离散度指标，其计算流程如下：逐像元依次提取增强型植被指数大于或等于第75分位值的所有数据，对应为该年份的modisevi年内时序数据的高值区域p。基于该年份增强型植被指数时序数据的高值区域p，计算该高值区域p内数据的极差和标准差sdm。依据植被指数高值区域内数据的极差和标准差sdm，构建时序离散度(temporaldispersion)指标td。
其中，td为植被指数中高值区域内数据的极差与标准差的乘积，其计算公式为：td＝(max-quan3)(sdm)。其中，max为对应年份增强型植被指数时序数据的最大值，quan3为该年份modisevi年内时序数据的第三分位数，即第75分位值，sdm为该年份modisevi年内时序数据高值区域p内的标准差。
进一步的，基于modisevi年内时序数据的最大值max和最小值min，提取生长期长度，其计算流程如下：逐像元提取年内modisevi时序数据的最大值max和最小值min，求取极差(max-min)，将数值达到极差的一半作为植被开始生长的依据，并将按时间顺序依次检索获取的大于或等于(min+(max-min)/2)的数据对应的持续时间累加，获得生长期长度。
步骤s03：设计高分位持续度、低分位持续度指标
进一步的，基于年内逐日modisevi时序数据的第65分位值，提取高分位持续度指标，其计算流程如下：逐像元按时间顺序依次检索出现大于或等于第65分位值的持续时间，依次记录为t1，t2，……，tn。求算所有的高值持续时间{t1，t2，……，tn}的最大值，获得高分位持续度指标。
进一步的，基于年内逐日modisevi时序数据的第50分位值，提取低分位持续度指标，其计算流程如下：逐像元按时间顺序依次检索出现大于或等于第50分位值的持续时间，依次记录为t1，t2，……，tn。求算所有的高值持续时间{t1，t2，……，tn}的次大值，获得低分位持续度指标。
步骤s04：建立五个时序指标的多年时序数据集
进一步的，基于2001-2016年多年逐日modisevi时序数据，逐像元逐年依次提取丰度、时序离散度、生长期长度、高分位持续度、低分位持续度指标。针对每个像元，建立以年份为步长的包括丰度、时序离散度、生长期长度、高分位持续度、低分位持续度在内五个时序指标的时序数据集。
在本实施例中，基于该时序数据集，以a(单季农作物变为森林)和b(多季农作物变为森林)为例，所形成的2001-2016年丰度、时序离散度、生长期长度、高分位持续度、低分位持续度指标的时序曲线图见图4以及图5所示。相比而言，单季农作物时序离散度和高分位持续度均较高，但植被丰度、生长期长度以及低分位持续度均较低。而多季农作物的植被丰度、低分位持续度均较高，而其他三个时序指标较低。森林植被的植被丰度、高分位持续度以及生长期长度均较高，并且其他两个时序指标较低。因此，在耕地植被(单季农作物、多季农作物)转变为森林植被时，其相应的五个时序指标呈现出对应的变化趋势。
步骤s05：计算五个时序指标的多年变化趋势q。
在本实施例中，利用sen氏斜率法，逐像元检测丰度、时序离散度、生长期长度、高分位持续度、低分位持续度等指标的多年变化趋势。sen氏斜率法能很好地避免数据异常值、数据空间分布形态对结果的干扰，对于较短的时间序列趋势分析更具可靠性。因此采用sen氏斜率法，逐像元分别计算丰度、时序离散度、生长期长度、高分位持续度、低分位持续度等多年时序数据的变化趋势q。
步骤s06：开展五个时序指标变化趋势q的显著性检验。
在本实施例中，当时序指标多年变化趋势q>0时，表示该时间序列呈一定的上升趋势；当q<0时，表示该时间序列具有一定的下降趋势。变化趋势q的绝对值越大，表示变化幅度越大。为了判断每个像元对应的丰度、时序离散度、生长期长度、高分位持续度、低分位持续度等时序指标的变化趋势q是否显著，需要进一步从空间和时序两个方面进行显著性检验。
进一步的，在空间显著性检验方面，其目的是为了检验丰度、时序离散度、生长期长度、高分位持续度、低分位持续度等时序指标的变化趋势q在所处研究区域内是否显著。
以丰度为例进行说明。假设研究区域内丰度的变化趋势q满足正态分布n(μ，σ2)，在给定的α显著性水平上，如果|q|≥q1-α/2，则该像元的丰度指标的变化趋势q在研究区域内显著。
进一步的，在时序显著性检验方面，其目的是检验丰度、时序离散度、生长期长度、高分位持续度、低分位持续度这五个时序指标的变化趋势q在研究时段内是否显著。
在本实施例中，时序显著性检验采用mann-kendall方法。mann-kendall方法具有对时序数据序列分布无要求，对异常值不敏感等优点，因此采用mann-kendall方法能很好地检验时序数据序列的变化趋势。
进一步的，综合sen氏斜率法计算获得的变化趋势q，综合空间显著性和时序显著性检验结果，将五个时序指标的变化趋势分为三类：显著正趋势(显著的上升趋势)、无趋势(不变)和显著负趋势(显著的下降趋势)。逐像元记录五个时序指标的多年变化趋势检测结果，建立相应的变化趋势q的空间分布图，如图6所示。
步骤s07：建立退耕还林自动识别流程图。
进一步的，如图7所示，基于每个像元的丰度、时序离散度、生长期长度、高分位持续度、低分位持续度的变化趋势特征，建立退耕还林自动识别流程图。
相对耕地农作物植被而言，森林的丰度较高，时序离散度较低，高分位持续度较高而低分位持续度较低，生长期长度较高。所建立退耕还林自动识别流程图如下：
(1)如果丰度和生长期长度均具有显著正趋势，而时序离散度呈现出显著负趋势，并且高分位持续度、低分位持续度均无趋势，则该像元由单季农作物转变为森林，即该像元为退耕还林区域。
(2)如果高分位持续度、生长期长度均具有显著正趋势，而低分位持续度呈现出显著负趋势，并且丰度、时序离散度均无趋势，则该像元由多季农作物转变为森林，即该像元为退耕还林区域。
步骤s08：获得研究区退耕还林空间分布图。
基于所建立的退耕还林自动识别流程，逐像元提取退耕还林区域，最终生成研究区退耕还林空间分布图。依据上述流程，可实现较精确的退耕还林自动提取。进一步的，如图8所示，为以山东省为例，综合研究区退耕还林自动提取结果，获得的退耕还林空间分布图。
以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。