CN105718936A

CN105718936A - 一种森林动态变化模式自动提取方法

Info

Publication number: CN105718936A
Application number: CN201610069963.6A
Authority: CN
Inventors: 邱炳文; 刘哲
Original assignee: Fuzhou University
Current assignee: Fuzhou University
Priority date: 2016-02-02
Filing date: 2016-02-02
Publication date: 2016-06-29

Abstract

本发明涉及一种森林动态变化模式自动提取方法，建立研究区的多年逐日增强型植被指数以及积雪指数的时序数据集，逐像元基于增强型植被指数时序数据，逐年提取植被覆盖强度、离散度、持续度等指标，并逐年提取亮度指标，逐像元依次检测上述指标的多年变化趋势，在此基础上基于植被覆盖强度、离散度、持续度、亮度等指标的变化趋势，建立森林动态变化模式识别流程图，最终达到森林动态变化自动监测的目的。该方法分别从覆盖时间、平均状态、变化幅度等侧面，设计若干指标，基于这些指标的多年总体变化趋势，有效地识别森林动态变化模式，不需要依赖已知样区的训练数据，不人机交互，具有鲁棒性好、分类精度高、自动化程度与抗干扰能力强等特点。

Description

一种森林动态变化模式自动提取方法

技术领域

本发明涉及遥感分类技术领域，特别是一种森林动态变化模式自动提取方法。

背景技术

森林作为地球之肺，源源不断地输送给我们所需的氧气和大量的有机物。人类的生存发展与森林生态系统息息相关。然而，近现代以来，由于人类不节制的砍伐，直接导致森林大量减少。我国是森林贫乏的国家，人均森林面积不足世界的1/4。在21世纪我国经济快速发展时期，对木材的需求量将持续增长，因此高效快速监测森林的动态变化对于自然资源的合理利用与保护至关重要。传统的森林调查方法，难以满足现势性以及时空连续监测等各方面的需求，因此基于遥感时序数据的森林动态监测已成为主流发展方向。

在森林遥感动态监测方面，比较有名的方法有LandTrendr(Landsat-basedDetectionofTrendsinDisturbanceandRecovery)和BFAST(BreakDetectionForAdditiveandTrend)方法。LandTrendr方法利用Landsat遥感影像时序影像，通过时间序列分割与重构，实现包括林火、病虫害等因素带来的森林干扰监测。在BFAST方法中，首先将植被演变时间序列分解为渐进的线性变化趋势、季节变化以及突变等三方面，并分别赋予一定的含义，如将前两者映射为干扰(火灾和虫灾)和植被物候变化(如土地覆盖变化)，然后通过迭代检测突变点，实现森林动态变化监测。这些基于时间序列分析的森林变化监测方法，为实现森林的时空连续变化监测提供了新思路。众多因素表明，基于遥感影像时空序列数据分析技术的森林动态监测，将成为未来的主流发展方向。现有的相关方法虽然展露出此类方法的特色优势，但在以下几个方面仍存在一些不足之处：(1)算法比较复杂，如通常涉及时间序列的分解和重构；(2)算法自动化程度有待提高；(3)多以森林增加或减少为主题，无法有效地识别森林动态变化模式。本发明拟从植被覆盖演变趋势分析的角度，自动提取森林动态变化模式，达到森林高效快速监测的目的。

发明内容

本发明的目的在于提供一种森林动态变化模式自动提取方法，适用于森林动态监测或其他植被类型动态变化的遥感自动监测领域，以克服现有技术中存在的缺陷。

为实现上述目的，本发明的技术方案是：一种森林动态变化模式自动提取方法，包括如下步骤，

S01：逐像元建立研究区多年逐日增强型植被指数以及积雪指数时序数据集；

S02：逐像元逐年分别提取植被覆盖强度指标、离散度指标以及持续度指标；

S03：逐像元逐年提取亮度指标；

S04：逐像元检测植被覆盖强度指标、离散度指标、持续度指标以及亮度指标的多年变化趋势；

S05：建立森林动态变化模式识别流程图；

S06：实现森林变化模式自动提取，获得研究区森林动态变化分布图。

在本发明一实施例中，在所述步骤S01中，利用MODIS波段反射率数据，通过如下方式建立研究区每个栅格单元的多年逐日增强型植被指数MODISEVI以及归一化差分积雪指数MODISNDSI时序数据，计算MODISEVI2以及NDSI时序数据，分别为：

E V I 2 = \frac{2.5 * (N i r - Re d)}{(N i r + Re d + 1)},

N D S I = \frac{G r e e n - {SWIR}_{6}}{G r e e n + {SWIR}_{6}},

其中，Nir，Red，Green，SWIR6分别为MODIS的近红外、红光、绿光以及中红外6波段的反射率；基于无云日计算得到MODISEVI以及NDSI时序数据，利用线性插值方法获得研究区多年的逐日MODISEVI以及NDSI时序数据，采用Whittakersmoother数据平滑方法，逐像元构建多年时空连续的逐日MODISEVI和NDSI时序数据集。

在本发明一实施例中，在所述步骤S02中，逐像元逐年基于植被指数时序数据的高值区域的均值，构建该像元每年的植被覆盖强度指标；逐像元逐年基于植被指数时序数据的中高值区域的极差以及标准差，构建该像元每年的植被覆盖离散度指标，其中，DM＝(Max-Quan₂)(SD_m)，DM为植被指数中高值区域内数据的极差与标准差的乘积，Max为对应年份增强型植被指数时序数据的最大值，Quan₂为对应年份增强型植被指数时序数据的二分位数；逐像元逐年基于植被指数时序曲线上大于等于预设阈值的持续时间，构建该像元每年的植被覆盖持续度指标。

在本发明一实施例中，在所述步骤S03中，逐像元逐年基于归一化差分积雪指数的年均值，构建该像元每年的亮度指标。

在本发明一实施例中，在所述步骤S04中，基于Sen氏斜率法和Mann-Kendall方法，逐像元检测植被覆盖强度指标、离散度指标、持续度指标以及亮度指标的变化斜率，结合空间分布显著性检验结果以及时间序列趋势显著性检验结果，将多年变化趋势的检测结果分为三类：显著的正趋势、没有趋势以及负趋势。

在本发明一实施例中，在所述步骤S05中，依据植被覆盖强度指标、离散度指标、持续度指标以及亮度指标的变化趋势特征，建立森林动态变化模式识别流程图。

在本发明一实施例中，在所述森林动态变化模式识别流程图中，结合所述步骤S04的检测结果，判断植被覆盖强度与离散度的变化趋势是否一致；如果植被覆盖强度以及离散度均存在显著的正趋势，由此判断该像元由其他植被转换为森林，森林增加，并进一步判断为以下两种情况：(1)如果植被覆盖亮度以及持续度均没有显著的变化趋势，则该像元由非植被转换为森林；(2)如果植被覆盖持续度不具有显著的变化趋势，而亮度具有显著的下降趋势，则该像元由稀疏植被转换为森林。

在本发明一实施例中，在所述森林动态变化模式识别流程图中，结合所述步骤S04的检测结果，判断植被覆盖强度与离散度的变化趋势是否一致；如果植被覆盖强度以及离散度均存在显著的负趋势，由此判断该像元由森林转换为其他植被类型，森林减少，并进一步判断为以下两种情况：(1)如果植被覆盖持续度以及亮度均不没有显著的变化趋势，则该像元由森林转换为稀疏植被；(2)如果植被覆盖亮度呈显著的上升趋势，而植被覆盖持续度不具有显著的变化趋势，则该像元由森林转换为非植被。

在本发明一实施例中，在所述森林动态变化模式识别流程图中，结合所述步骤S04的检测结果，判断植被覆盖强度与离散度的变化趋势是否一致；如果植被覆盖强度与离散度的变化趋势不一致，且在植被覆盖强度以及亮度均不具备显著的变化趋势的前提下，进一步分为以下四种情况：(1)植被覆盖离散度不变，而植被覆盖持续度呈显著的上升趋势，则该像元由多季农作物转换为森林；(2)植被覆盖持续度不具备显著的变化趋势，而植被覆盖离散度呈显著的下降趋势，则该像元由单季农作物转换为森林；(3)植被覆盖持续度不具有显著的变化趋势，而离散度呈显著的上升趋势，则该像元由森林转换为单季农作物；(4)植被覆盖离散度不具有显著的变化趋势，而植被覆盖持续度呈显著的下降趋势，则该像元由森林转换为多季农作物。

在本发明一实施例中，基于所述森林动态变化模式识别流程，逐像元提取森林动态变化模式，生成研究区森林动态变化模式分布图。

相较于现有技术，本发明具有以下有益效果：

(1)基于植被指数时序数据，分别从覆盖时间、平均状态、变化幅度等侧面，设计测植被覆盖强度、离散度、持续度等指标，用于表征每年植被覆盖多方位特征，一定程度上避免了时间序列分析与重构的繁琐。

(2)基于归一化差分积雪指数时序数据，设计亮度指标，增加了植被变化监测的维度，更进一步弥补了单一植被指数信息量不足的局限性。

(3)基于植被覆盖强度、离散度、连续度、亮度等指标的总体变化趋势，对某些年份数据出现错误或异常等情况，具有更好的容错能力。

(4)能有效地识别森林变化模式，不仅检测出森林出现变化两种情况(增加或减少)，同时能找到森林增加或减少的来源或去向。

(5)可以不借助已知森林变化样区的训练数据，不需要人机交互，不依赖监督人类或机器学习方法，顺利地实现森林变化模式自动识别。

附图说明

图1为本发明中森林动态变化模式自动提取方法的实现流程图。

图2为本发明一实施例中2001-2014年植被MODISEVI和NDSI指数时序曲线图。

图3为本发明一实施例中2001-2014年植被覆盖强度、离散度、持续度、亮度指标的变化趋势图。

图4为本发明一实施例中2001-2014年研究区植被覆盖强度、离散度、持续度、亮度指标变化趋势空间分布图。

图5本发明中森林动态变化模式自动提取方法的森林动态变化模式识别流程图。

图6为本发明一实施例中研究区森林动态变化分布图。

具体实施方式

下面结合附图，对本发明的技术方案进行具体说明。

本发明一种森林动态变化模式自动提取方法，如图1所示：首先建立研究区每个栅格单元的多年逐日增强型植被指数以及积雪指数的时序数据集，然后逐像元基于增强型植被指数时序数据，逐年提取植被覆盖强度、离散度、持续度等指标，同时基于多年逐日积雪指数的时序数据集逐年提取亮度指标，逐像元依次检测上述指标的多年变化趋势，在此基础上基于植被覆盖强度、离散度、持续度、亮度等指标的变化趋势，建立森林动态变化模式识别流程图，最终达到森林动态变化自动监测的目的。

进一步的，在本实施例中，一种森林动态变化模式自动提取方法，包括以下步骤：

步骤S01：建立研究区2001-2014年的多年逐日增强型植被指数MODISEVI和归一化差分积雪指数MODISNDSI时序数据。

利用MODIS波段反射率数据，计算MODISEVI2和NDSI时序数据，计算公式分别为：

E V I 2 = \frac{2.5 * (N i r - Re d)}{(N i r + Re d + 1)}

N D S I = \frac{G r e e n - {SWIR}_{6}}{G r e e n + {SWIR}_{6}}

其中Nir，Red，Green，SWIR6分别为MODIS的近红外、红光、绿光以及中红外6波段的反射率。基于无云日计算得到的MODISEVI和NDSI时序数据，利用线性插值方法获得研究区多年的逐日MODISEVI和NDSI时序数据，然后采用Whittakersmoother等数据平滑方法，逐像元构建多年时空连续的逐日MODISEVI和NDSI时序数据集。基于该数据集，以A，也即森林种植与抚育、农作物变为森林，和B，也即森林减少、森林变为城市，为例，所形成的多年植被MODISEVI和NSDI指数时序曲线图见图2。

进一步的，在本实施例中，MODIS数据：中分辨率成像光谱仪数据，全称为ModerateResolutionImagingSpectroradiometer。

植被指数：植被指数是表征植被生长状态以及空间分布密度的因子。常见的植被指数有NDVI和EVI。NDVI为归一化植被指数，全称为NormalizedDifferenceVegetationIndex。EVI为增强型植被指数，全称为EnhancedVegetationIndex。

归一化差分积雪指数(NormalizedDifferenceSnowIndex,NDSI)：能指示地表亮度变化，积雪指数可用于识别积雪和冰等高亮度地物。

多年逐日植被指数时序数据/时序曲线：从起始年份的元旦开始，按时间顺序，一直到结束年份的年末，逐日记录所形成的植被指数的数据序列/时序曲线。

步骤S02：基于多年逐日MODISEVI时序数据，逐年提取植被覆盖强度、离散度、持续度等指标。

逐像元逐年依次提取MODISEVI时序数据的重要分位值，进一步的，在本实施例中，可选择第65分位值、第75分位值。

基于多年逐日MODISEVI时序数据，逐像元基于每年MODISEVI时序数据的第75分位值，建立植被覆盖强度指标，其计算流程如下：逐像元逐年依次提取增强型植被指数大于或等于第75分位值的所有数据，逐年求取其年平均值，作为每年的植被覆盖强度数据。

基于多年逐日MODISEVI时序数据，逐像元逐年基于每年MODISEVI时序数据的第65分位值，提取植被覆盖离散度指标，其计算流程如下：逐像元逐年依次提取增强型植被指数大于或等于第65分位值的所有数据，分别对应为每个年份的增强型植被指数时序数据的中高值区域M。基于每个年份的增强型植被指数时序数据的中高值区域M，逐年分别计算计算植被指数中高值区域内数据的极差和标准差SD_m，从而构建每年的植被覆盖离散度(Dispersion)指标DM。DM为植被指数中高值区域内数据的极差与标准差的乘积，其公式为：DM＝(Max-Quan₂)(SD_m)，其中Max为对应年份增强型植被指数时序数据的最大值，Quan₂为对应年份增强型植被指数时序数据的二分位数。

基于多年逐日MODISEVI时序数据，逐像元逐年基于每年MODISEVI时序数据的第65分位值，提取植被覆盖持续度指标，其计算流程如下：逐像元按时间顺序依次检索出现大于或等于第65分位值的持续时间，依次记录为t₁，t₂，……，t_n。求算所有的高值持续时间{t₁，t₂，……，t_n}的最大值，获得植被覆盖持续度指标。

步骤S03：基于多年逐日积雪指数的时序数据，逐年提取亮度指标。

基于多年逐日积雪指数的时序数据，逐像元逐年依次计算积雪指数的年均值，作为该年份亮度指标的数据。分别对植被覆盖强度、离散度、持续度以及亮度指标，逐像元逐年依次提取其数据，建立以年份为步长的植被覆盖强度、离散度、持续度以及亮度时序数据集。基于该数据集，以A，也即森林种植与抚育、农作物变为森林，和B，也即森林减少、森林变为城市，为例，所形成的2002-2014年植被覆盖强度、离散度、持续度、亮度指标的时序曲线图见图3。

步骤S04：检测每个像元的植被覆盖强度、离散度、持续度以及亮度等指标的多年变化趋势。

将Sen氏斜率法和Mann-Kendall方法相结合，检测每个像元的植被覆盖强度、离散度、持续度以及亮度等指标的多年变化趋势。Sen氏斜率法能很好地避免数据异常值、数据空间分布形态对结果的干扰，对于较短的时间序列趋势分析更具可靠性，但需要结合Mann-Kendall方法对时间序列变化趋势做进一步的显著性检验。Mann-Kendall方法具有对时序数据序列分布无要求，对异常值不敏感等优点，因此采用Mann-Kendall方法能很好的检验时序数据序列的变化趋势。

首先采用Sen氏斜率法，逐像元分别计算植被覆盖强度、离散度、持续度以及亮度等时序数据的变化斜率Q。当Q>0时，表示该时间序列呈一定的上升趋势；当Q<0时，表示该时间序列具有一定的下降趋势。变化斜率Q的绝对值越大，表示变化幅度越大。为了判断每个像元对应的植被覆盖强度、离散度、持续度以及亮度等时序数据的变化斜率Q是否明显，需要进一步从空间分布和时间序列趋势这两个方面进行显著性检验。在变化斜率Q的空间分布方面，接下来需要进一步检测植被覆盖强度、离散度、持续度以及亮度等时序数据的变化斜率Q在研究区域内是否显著。以植被覆盖强度为例进行说明，假设研究区域内植被覆盖强度的变化斜率Q满足正态分布N(μ，σ²)，在给定的α显著性水平上，进一步判断每个像元对应植被覆盖强度指标的变化斜率Q所处数值范围，如果|Q|≥Q_1-α/2，则该像元的植被覆盖强度指标的变化斜率Q在研究区域内显著。

在时间序列趋势显著性方面，基于Mann-Kendall方法，逐像元分别依次检测植被覆盖强度、离散度、持续度以及亮度等指标的时序数据的变化趋势是否显著。综合Sen氏斜率法计算获得的变化斜率Q，结合空间分布和时间序列趋势两个方面获得的显著性检验结果，将多年变化趋势的检测结果，分为三类：显著的正趋势、没有趋势和显著的负趋势，且分别对应显著的上升趋势、不变以及显著的下降趋势。将每个像元中植被覆盖强度、离散度、持续度以及亮度等指标的多年变化趋势检测结果保存在数据库中。基于该数据集，所形成的研究区植被覆盖强度、离散度、持续度以及亮度的变化趋势的空间分布图如图4所示。

步骤S05：建立森林动态变化模式识别流程图。

基于每个像元的植被覆盖强度、离散度、持续度以及亮度的变化趋势特征，建立森林动态变化模式识别流程图，如图5所示。相对其他植被而言，森林的植被覆盖持续度较高，离散度较低。相对于非植被而言，森林的植被覆盖强度明显偏高。所建立森林动态变化模式识别流程图如下：

(1)如果植被覆盖强度和离散度均具有正趋势，并且植被覆盖亮度和持续度均不具备显著的变化趋势，则该像元由非植被转换为森林；如果植被覆盖强度和离散度均具有正趋势，植被覆盖持续度不存在显著的变化趋势，亮度具有下降的变化趋势，则该像元由稀疏植被转换为森林。

(2)如果植被覆盖强度和离散度均具有负趋势，并且植被覆盖持续度和亮度均不具有显著的变化趋势，则该像元由森林转换为稀疏植被；如果植被覆盖强度和离散度均具有负趋势，植被覆盖亮度呈上升趋势而植被覆盖持续度不具备显著的变化趋势，则该像元由森林转换为非植被。

(3)如果植被覆盖强度和亮度均没有显著的变化趋势，并且植被覆盖离散度不变而植被覆盖持续度呈上升趋势，则该像元由多季农作物转换为森林；如果植被覆盖强度、亮度和持续度均没有显著的变化趋势，而植被覆盖离散度呈下降趋势，则该像元由单季农作物转换为森林；如果植被覆盖强度、亮度和持续度均没有显著的变化趋势，而离散度呈上升趋势，则该像元由森林转换为单季农作物；如果植被覆盖强度、亮度和离散度均没有显著的变化趋势，而植被覆盖持续度呈下降趋势，则该像元由森林转换为多季农作物。

步骤S06：实现森林变化模式自动提取，获得研究区森林动态变化分布图。

基于所建立的森林动态变化模式识别流程，逐像元提取森林动态变化模式，最终生成研究区森林动态变化模式分布图。依据上述流程，可实现较精确的森林变化模式自动提取。以河南省为例，综合研究区森林动态变化模式自动提取结果，获得的森林动态变化分布图，如图6所示。

以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。

Claims

1.一种森林动态变化模式自动提取方法，其特征在于：包括如下步骤，

S03：逐像元逐年提取亮度指标；

S05：建立森林动态变化模式识别流程图；

2.根据权利要求1所述的一种森林动态变化模式自动提取方法，其特征在于：在所述步骤S01中，利用MODIS波段反射率数据，通过如下方式建立研究区每个栅格单元的多年逐日增强型植被指数MODISEVI以及归一化差分积雪指数MODISNDSI时序数据，计算MODISEVI2以及NDSI时序数据，分别为：

E V I 2 = \frac{2.5 * (N i r - Re d)}{(N i r + Re d + 1)},

N D S I = \frac{G r e e n - {SWIR}_{6}}{G r e e n + {SWIR}_{6}},

3.根据权利要求1所述的一种森林动态变化模式自动提取方法，其特征在于：在所述步骤S02中，逐像元逐年基于植被指数时序数据的高值区域的均值，构建该像元每年的植被覆盖强度指标；逐像元逐年基于植被指数时序数据的中高值区域的极差以及标准差，构建该像元每年的植被覆盖离散度指标，其中，DM＝(Max-Quan₂)(SD_m)，DM为植被指数中高值区域内数据的极差与标准差的乘积，Max为对应年份增强型植被指数时序数据的最大值，Quan₂为对应年份增强型植被指数时序数据的二分位数；逐像元逐年基于植被指数时序曲线上大于等于预设阈值的持续时间，构建该像元每年的植被覆盖持续度指标。

4.根据权利要求1所述的一种森林动态变化模式自动提取方法，其特征在于：在所述步骤S03中，逐像元逐年基于归一化差分积雪指数的年均值，构建该像元每年的亮度指标。

5.根据权利要求1所述的一种森林动态变化模式自动提取方法，其特征在于：在所述步骤S04中，基于Sen氏斜率法和Mann-Kendall方法，逐像元检测植被覆盖强度指标、离散度指标、持续度指标以及亮度指标的变化斜率，结合空间分布显著性检验结果以及时间序列趋势显著性检验结果，将多年变化趋势的检测结果分为三类：显著的正趋势、没有趋势以及负趋势。

6.根据权利要求1所述的一种森林动态变化模式自动提取方法，其特征在于：在所述步骤S05中，依据植被覆盖强度指标、离散度指标、持续度指标以及亮度指标的变化趋势特征，建立森林动态变化模式识别流程图。

7.根据权利要求6所述的一种森林动态变化模式自动提取方法，其特征在于：在所述森林动态变化模式识别流程图中，结合所述步骤S04的检测结果，判断植被覆盖强度与离散度的变化趋势是否一致；如果植被覆盖强度以及离散度均存在显著的正趋势，由此判断该像元由其他植被转换为森林，森林增加，并进一步判断为以下两种情况：(1)如果植被覆盖亮度以及持续度均没有显著的变化趋势，则该像元由非植被转换为森林；(2)如果植被覆盖持续度不具有显著的变化趋势，而亮度具有显著的下降趋势，则该像元由稀疏植被转换为森林。

8.根据权利要求6所述的一种森林动态变化模式自动提取方法，其特征在于：在所述森林动态变化模式识别流程图中，结合所述步骤S04的检测结果，判断植被覆盖强度与离散度的变化趋势是否一致；如果植被覆盖强度以及离散度均存在显著的负趋势，由此判断该像元由森林转换为其他植被类型，森林减少，并进一步判断为以下两种情况：(1)如果植被覆盖持续度以及亮度均不没有显著的变化趋势，则该像元由森林转换为稀疏植被；(2)如果植被覆盖亮度呈显著的上升趋势，而植被覆盖持续度不具有显著的变化趋势，则该像元由森林转换为非植被。

9.根据权利要求6所述的一种森林动态变化模式自动提取方法，其特征在于：在所述森林动态变化模式识别流程图中，结合所述步骤S04的检测结果，判断植被覆盖强度与离散度的变化趋势是否一致；如果植被覆盖强度与离散度的变化趋势不一致，且在植被覆盖强度以及亮度均不具备显著的变化趋势的前提下，进一步分为以下四种情况：(1)植被覆盖离散度不变，而植被覆盖持续度呈显著的上升趋势，则该像元由多季农作物转换为森林；(2)植被覆盖持续度不具备显著的变化趋势，而植被覆盖离散度呈显著的下降趋势，则该像元由单季农作物转换为森林；(3)植被覆盖持续度不具有显著的变化趋势，而离散度呈显著的上升趋势，则该像元由森林转换为单季农作物；(4)植被覆盖离散度不具有显著的变化趋势，而植被覆盖持续度呈显著的下降趋势，则该像元由森林转换为多季农作物。

10.根据权利要求1一种森林动态变化模式自动提取方法，其特征在于：基于所述森林动态变化模式识别流程，逐像元提取森林动态变化模式，生成研究区森林动态变化模式分布图。