CN111832506B - 一种基于长时序植被指数的重建植被遥感判别方法 - Google Patents

一种基于长时序植被指数的重建植被遥感判别方法 Download PDF

Info

Publication number
CN111832506B
CN111832506B CN202010701602.5A CN202010701602A CN111832506B CN 111832506 B CN111832506 B CN 111832506B CN 202010701602 A CN202010701602 A CN 202010701602A CN 111832506 B CN111832506 B CN 111832506B
Authority
CN
China
Prior art keywords
vegetation
ndvi
year
fitting
data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010701602.5A
Other languages
English (en)
Other versions
CN111832506A (zh
Inventor
米家鑫
马占元
赵晨德
张绍良
常晓华
路喜
王玉明
魏川博
岳雪莲
杨永均
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Mining and Technology CUMT
Datong Coal Mine Group Co Ltd
Original Assignee
China University of Mining and Technology CUMT
Datong Coal Mine Group Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Mining and Technology CUMT, Datong Coal Mine Group Co Ltd filed Critical China University of Mining and Technology CUMT
Priority to CN202010701602.5A priority Critical patent/CN111832506B/zh
Publication of CN111832506A publication Critical patent/CN111832506A/zh
Application granted granted Critical
Publication of CN111832506B publication Critical patent/CN111832506B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/188—Vegetation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/02—Agriculture; Fishing; Forestry; Mining

Landscapes

  • Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • Strategic Management (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Economics (AREA)
  • Theoretical Computer Science (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Educational Administration (AREA)
  • Development Economics (AREA)
  • General Business, Economics & Management (AREA)
  • Health & Medical Sciences (AREA)
  • Tourism & Hospitality (AREA)
  • General Health & Medical Sciences (AREA)
  • Marketing (AREA)
  • Quality & Reliability (AREA)
  • Operations Research (AREA)
  • Game Theory and Decision Science (AREA)
  • Multimedia (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Agronomy & Crop Science (AREA)
  • Animal Husbandry (AREA)
  • Marine Sciences & Fisheries (AREA)
  • Mining & Mineral Resources (AREA)
  • Primary Health Care (AREA)
  • Image Processing (AREA)

Abstract

一种基于长时序植被指数的重建植被遥感判别方法，适用于矿山生态环境监测技术领域。通过植被的长时序植被指数，在进过噪声识别与预测值计算后基于傅里叶级数与最小二乘法构建植被生长曲线，进一步得到区域植被在时序年份内的年度植被指数均值，利用线性拟合计算每5年间植被指数的变化趋势，通过某一时点过去与未来5年间变化趋势的差异度作为判别指标，当差异度超过设置的阈值时则认定为重建植被且差异度最大的年份为植被重建起始时间。其有效反映了植被完整的生长过程从而减少了单一时点数据的片面性，噪声识别与预测值计算提高了时序数据的精确度，而植被生长趋势差异度作为判别指标减少了单一年份的异常数据对植被整体生长趋势的影响。

Description

一种基于长时序植被指数的重建植被遥感判别方法
技术领域
本发明涉及一种基于长时序植被指数的重建植被遥感判别方法，尤其适用于矿山生态环境监测技术领域使用。
背景技术
重建植被的快速判别是矿山土地复垦与生态修复的重要工作。过去20年来，矿山生态修复受到广泛关注，全国各地矿区开展了大量的植被重建工程。为了掌握植被重建工作落实情况，首先需要对矿区内的各种植被进行判别。然而传统的野外调查难以满足高效、准确、实时的监测要求，从而导致监管部门无法及时地评估矿山植被重建的有效性。
针对大范围的重建植被，遥感技术是实现判别的主要方法，目前的判别方法主要是通过重建植被在遥感影像中的光谱特征为判别依据，从而与其他类型植被进行区分然后进行判别和提取，主要包括遥感影像分类法、土地利用变化检测法、植被指数阈值法。遥感影像分类法主要是利用无人机或卫星遥感影像，依靠重建植被在光谱特征或形态结构上的特征，通过监督分类或非监督分类算法确定重建植被的面积和分布；土地利用变化检测法则是在遥感影像分类法的基础上，通过对多期遥感影像进行分类，确定植被覆盖发生变化的区域并将其认定为重建植被区；植被指数阈值法则是利用遥感影像计算区域内的植被指数，并确定重建植被的植被指数变化范围，然后设置植被指数阈值从而提取出区域内的重建植被。
目前的重建植被遥感判别方法十分依赖重建植被与其他类型植被的光谱特征差异，然而植被之间的光谱特征差异并不明显，此外单一时点的遥感影像常常存在数据噪声、异常气候、人类扰动造成的影响，无法保证判别结果的可靠性。
发明内容
针对上述技术的不足之处，提供一种步骤简单，辨别效果精准，满足矿山生态修复监管和评价的需要的基于长时序植被指数的重建植被遥感判别方法。
为解决上述问题，本发明基于长时序植被指数的重建植被遥感判别方法，包括以下步骤：
S1.利用卫星采集需要识别的研究区域的植被信息，生成具有明确边界的制备矢量图层遥感影像；
S2.获取研究区域的历年内所有时序植被指数遥感影像图，统计每期遥感影像图中拟识别植被区域的归一化植被指数NDV均值，利用时序NDVI均值建立时序数据集；
S3.利用均值迭代法对研究区域的归一化植被指数进行噪声识别，去除遥感影像中冰雪覆盖及云雾遮挡的数据噪声，并利用相邻时点数据均值作为预测值替代从而剔除噪声；
S4.对剔除噪声后的时序NDVI数据基于傅里叶级数与普通最小二乘法进行拟合，获得研究区域植被在不同年份的生长曲线；
S5.基于植被生长曲线计算研究区域植被年度均值，并从不同年份的生长曲线中提取植被关键参数，包括5年NDVI均值线性拟合方程的斜率及决定系数，并基于植被关键参数提出植被生长趋势差异度以反映植被的持续变化情况；
S6.将植被生长趋势差异度作为阈值判别植被类型，从而将遥感影像图中人工种植的重建植被与演替发育的自然植被进行区分，并通过植被生长趋势差异度的最大值获得重建植被的种植时间。
所述步骤S1具体方法为：通过调查时点的卫星遥感影像或植被覆盖现状图，在ArcMap或Super Map地理信息软件中人工选择或者导入拟判别植被的边界范围，将选择的植被的边界范围保存为具有坐标系的矢量图层，植被区域的矢量图层为UTM坐标系。
所述步骤S2的具体方法为：利用Landsat系列遥感影像的归一化植被指数NDVI时序数据，包括Landsat4/5TM，Landsat7 ETM及Landsat8 OLC影像的8天间隔NDVI数据集，在ArcMap或ENVI软件中通过批量裁剪功能提取并输出研究区域植被的NDVI数据均值，分别对获取的历年所有遥感影像数据进行统计后从而建立判别区域的时序植被指数数据集。
所述步骤S3的具体方法为：比较第i期NDVI数据的观测值NOi，i＞2且为自然整数，与其前后相邻两期植被指数中最大值的均值进行比较，若小于其均值的50％则认为观测值为噪声，并用其均值作为预测值NPi进行代替，反之则保留其观测值，然后持续重复该步骤，直到消除所有的数据噪声，具体筛选规则为
ELSE NPi＝NOi
其中NOi为第i期NDVI数据的观测值，NPi为第i期NDVI数据预测值，MAX()为识别输入数据中的最大值并将其输出的函数，NOi-2，NOi-1，NOi+1和NOi+2分别为第i-2，i-1，i+1，和i+2期NDVI数据的观测值。
所述步骤S4的具体方法为：将剔除噪音后的时序NDVI数据根据所在年份进行分类，并将其日期转化为年纪日形式(n/365)，然后分别对各个年份的NDVI数据通过傅里叶级数与最小二乘法进行拟合，其具体公式为：
其中T为函数的周期，根据植被的生长周期将其设置为365天，x为第i期NDVI数据点年纪日，a，b1，b2为待确定的拟合方程常数项；
通过带入获取的NDVI观测值确定拟合方程的常数项，从而得到各个年份的NDVI拟合曲线，通过拟合方程与原始观测值计算拟合曲线的决定系数R2，具体公式为：
其中yi为第i期遥感影像的NDVI观测值，fi为第i期NDVI的拟合值，为该年份内NDVI观测值的均值；
R2反映了拟合方程对NDVI的拟合效果，R2如果过低则说明方程无法反映植被的真是变化，反之则说明效果较好，对于拟合曲线的决定系数R2低于0.8的年份进行二次检验，重复S3、S4步骤，直至获得的拟合曲线精度达到要求。
所述步骤S5的具体方法为：通过步骤S4中确定的不同年份的NDVI拟合曲线分别计算该年度的NDVI均值，
首先将年纪日形式的1至365天分别带入NDVI拟合曲线，求得该年份内所有的NDVI拟合值并计算该年份的NDVI均值NAj；
按照此方法依次求得所有年份的NDVI均值后，开始计算每5年NDVI均值的线性拟合方程并提取关键植被参数，包括5年NDVI均值线性拟合方程的斜率及决定系数：首先对第m年起当年至5年的NDVI拟合曲线均值进行线性拟合，其中5＜m＜n-5且为自然整数，n为NDVI时序数据的跨越年度，公式为：
r＝FTmX+b1 III
其中X为年份，Y为该年份的拟合曲线计算得到NDVI均值，FTm为方程斜率，b1为该线性拟合方程的常数项；
记录其线性拟合方程的斜率FTm以及决定系数FTR2’，决定系数FTR2的计算公式III如下所示：
其中NAi为第j年份的NDVI的均值，Yj为第j年份NDVI的线性拟合值，为参与该线性拟合方程的所有年份的NDVI均值；
类似的，对第m年之前的5年NDVI数据，进行线性拟合求得其拟合线性方程的斜率PTm与决定系数PTR2，公式为：
Y＝PTmX+b2 V
其中X为年份，Y为拟合曲线计算得到该年份的NDVI均值，PTm为方程斜率，b2为该线性拟合方程的常数项；
其中NAi为第j年份的NDVI的均值，Yj为第j年份NDVI的线性拟合值，为参与该线性拟合方程的所有年份的NDVI均值；
在获得过去5年与未来5年的NDVI均值线性拟合方程斜率之后，开始结算植被生长趋势差异度；首先通过前后5年NDVI均值线性拟合方程斜率的比值，并将两段拟合方程的决定系数均值作为权重，计算得到第m年前后植被生长趋势差异度Dm，公式如下所示：
由于NDVI拟合曲线均值的线性拟合方程斜率存在正负差异，因此在计算生长趋势差异度Dm时加入逻辑判断语句：当PTm小于0时说明该区域植被在未来5年处于下降趋势，明显不属于重建植被，从而对直接对其Dm赋值为-1；当FTm小于0而PTm大于0是说明该区域植被从第m年起NDVI从过去5年的下降趋势转变为增长趋势，具有较大的可能为重建植被，其可能性取决于拟合方程的准确度，因此将其Dm值赋值为当FTm与PTm均大于0时，则按照公式VII计算；逻辑判断式为：
IF PTm＜0
THEN Dm＝-1
ELSE IF FTm＜0
其中PTm为线性拟合公式III的斜率，FTm为线性拟合公式V的斜率，FTR2为拟合方程III的决定系数，PTR2为拟合公式V的决定系数，Dm为植被生长趋势差异度，
所述步骤S6的具体方法为：考虑到植被生长收到光照、降水、物候等多因素影像下，不同年度中NDVI存在一定的自然波动，因此将判别是否为重建植被的Dm阈值设定为1.5，即在未来5年植被增长的斜率较过去5年增长斜率提高50％，当判别区域植被在任意5年时间段中出现2个年度的Dm大于1.5时则认定这一区域植被为重建植被；反之，当Dm小于1.5时则认定这一区域植被不是重建植被，同时，Dm最大值出现的年份则被认为是植被重建的起始年份。
有益效果：
1)在时序植被指数拟合曲线过程中增加了噪声识别和预测值计算，从而提高了拟合方程的精度，同时避免了剔除噪声信息带来的植被数据缺失；
2)提出植被生长趋势差异度作为判别指标，量化了植被在多年期的趋势变化，减少了某一年度植被的异常变化对植被遥感判别结果的影响；
3)通过植被生长趋势差异度最大值的出现时间，准确确定监测区域的重建植被的栽种时间。
附图说明
图1为本发明基于长时序植被指数的重建植被遥感判别方法的流程示意图。
图2为本发明植被指数观测值经过噪声识别与预测值计算后的对比图。
图3为本发明植被指数的最小二乘法(OLS)拟合示意图。
图4为本发明植被指数生长趋势变化示意图。
具体实施方式
下面结合附图及实例对本发明做进一步说明。
如图1所示，本发明的基于长时序植被指数的重建植被遥感判别方法，包括以下步骤：
S1.利用卫星采集需要识别的研究区域的植被信息，生成具有明确边界的制备矢量图层遥感影像；通过调查时点的卫星遥感影像或植被覆盖现状图，在ArcMap或Super Map地理信息软件中人工选择或者导入拟判别植被的边界范围，将选择的植被的边界范围保存为具有坐标系的矢量图层，植被区域的矢量图层为UTM坐标系。
S2.获取研究区域的历年内所有时序植被指数遥感影像图，统计每期遥感影像图中拟识别植被区域的归一化植被指数NDV均值，利用时序NDVI均值建立时序数据集；利用Landsat系列遥感影像的归一化植被指数NDVI时序数据，包括Landsat4/5TM，Landsat7 ETM及Landsat8 OLC影像的8天间隔NDVI数据集，在ArcMap或ENVI软件中通过批量裁剪功能提取并输出研究区域植被的NDVI数据均值，分别对获取的历年所有遥感影像数据进行统计后从而建立判别区域的时序植被指数数据集。
S3.利用均值迭代法对研究区域的归一化植被指数进行噪声识别，去除遥感影像中冰雪覆盖及云雾遮挡的数据噪声，并利用相邻时点数据均值作为预测值替代从而剔除噪声；比较第i期NDVI数据的观测值NOi，i＞2且为自然整数，与其前后相邻两期植被指数中最大值的均值进行比较，若小于其均值的50％则认为观测值为噪声，并用其均值作为预测值NPi进行代替，反之则保留其观测值，然后持续重复该步骤，直到消除所有的数据噪声，具体筛选规则为
ELSE NPi＝NOi
其中NOi为第i期NDVI数据的观测值，NPi为第i期NDVI数据预测值，MAX()为识别输入数据中的最大值并将其输出的函数，NOi-2，NOi-1，NOi+1和NOi+2分别为第i-2，i-1，i+1，和i+2期NDVI数据的观测值。
S4.对剔除噪声后的时序NDVI数据基于傅里叶级数与普通最小二乘法进行拟合，获得研究区域植被在不同年份的生长曲线；将剔除噪音后的时序NDVI数据根据所在年份进行分类，并将其日期转化为年纪日形式(n/365)，然后分别对各个年份的NDVI数据通过傅里叶级数与最小二乘法进行拟合，其具体公式为：
其中T为函数的周期，根据植被的生长周期将其设置为365天，x为第i期NDVI数据点年纪日，a，b1，b2为待确定的拟合方程常数项；
通过带入获取的NDVI观测值确定拟合方程的常数项，从而得到各个年份的NDVI拟合曲线，通过拟合方程与原始观测值计算拟合曲线的决定系数R2，具体公式为：
其中yi为第i期遥感影像的NDVI观测值，fi为第i期NDVI的拟合值，为该年份内NDVI观测值的均值；
R2反映了拟合方程对NDVI的拟合效果，R2如果过低则说明方程无法反映植被的真是变化，反之则说明效果较好，对于拟合曲线的决定系数R2低于0.8的年份进行二次检验，重复S3、S4步骤，直至获得的拟合曲线精度达到要求。
S5.基于植被生长曲线计算研究区域植被年度均值，并从不同年份的生长曲线中提取植被关键参数，包括5年NDVI均值线性拟合方程的斜率及决定系数，并基于植被关键参数提出植被生长趋势差异度以反映植被的持续变化情况；
确定的不同年份的NDVI拟合曲线分别计算该年度的NDVI均值：
首先将年纪日形式的1至365天分别带入NDVI拟合曲线，求得该年份内所有的NDVI拟合值并计算该年份的NDVI均值NAj；
按照此方法依次求得所有年份的NDVI均值后，开始计算每5年NDVI均值的线性拟合方程并提取关键植被参数，包括5年NDVI均值线性拟合方程的斜率及决定系数：首先对第m年起当年至5年的NDVI拟合曲线均值进行线性拟合，其中5＜m＜n-5且为自然整数，n为NDVI时序数据的跨越年度，公式为：
Y＝FTmX+b1 III
其中X为年份，Y为该年份的拟合曲线计算得到NDVI均值，FTm为方程斜率，b1为该线性拟合方程的常数项；
记录其线性拟合方程的斜率FTm以及决定系数FTR2’，决定系数FTR2的计算公式III如下所示：
其中NAi为第j年份的NDVI的均值，Yj为第j年份NDVI的线性拟合值，为参与该线性拟合方程的所有年份的NDVI均值；
类似的，对第m年之前的5年NDVI数据，进行线性拟合求得其拟合线性方程的斜率PTm与决定系数PTR2，公式为：
Y＝PTmX+b2 V
其中X为年份，Y为拟合曲线计算得到该年份的NDVI均值，PTm为方程斜率，b2为该线性拟合方程的常数项；
其中NAi为第j年份的NDVI的均值，Yj为第j年份NDVI的线性拟合值，为参与该线性拟合方程的所有年份的NDVI均值；
在获得过去5年与未来5年的NDVI均值线性拟合方程斜率之后，开始结算植被生长趋势差异度；首先通过前后5年NDVI均值线性拟合方程斜率的比值，并将两段拟合方程的决定系数均值作为权重，计算得到第m年前后植被生长趋势差异度Dm，公式如下所示：
由于NDVI拟合曲线均值的线性拟合方程斜率存在正负差异，因此在计算生长趋势差异度Dm时加入逻辑判断语句：当PTm小于0时说明该区域植被在未来5年处于下降趋势，明显不属于重建植被，从而对直接对其Dm赋值为-1；当FTm小于0而PTm大于0是说明该区域植被从第m年起NDVI从过去5年的下降趋势转变为增长趋势，具有较大的可能为重建植被，其可能性取决于拟合方程的准确度，因此将其Dm值赋值为当FTm与PTm均大于0时，则按照公式VII计算；逻辑判断式为：
IF PTm＜0
THEN Dm＝-1
ELSE IF FTm＜0
其中PTm为线性拟合公式III的斜率，FTm为线性拟合公式V的斜率，FTR2为拟合方程III的决定系数，PTR2为拟合公式V的决定系数，Dm为植被生长趋势差异度，
S6.将植被生长趋势差异度作为阈值判别植被类型，从而将遥感影像图中人工种植的重建植被与演替发育的自然植被进行区分，并通过植被生长趋势差异度的最大值获得重建植被的种植时间。
考虑到植被生长收到光照、降水、物候等多因素影像下，不同年度中NDVI存在一定的自然波动，因此将判别是否为重建植被的Dm阈值设定为1.5，即在未来5年植被增长的斜率较过去5年增长斜率提高50％，当判别区域植被在任意5年时间段中出现2个年度的Dm大于1.5时则认定这一区域植被为重建植被；反之，当Dm小于1.5时则认定这一区域植被不是重建植被，同时，Dm最大值出现的年份则被认为是植被重建的起始年份。
实施例一、
S1.确定植被范围。
在本实例中，首先利用GoogleEarth软件定位至山西省大同市南郊区，参照最新的卫星遥感影像通过目视解译分别选择山西省大同市南郊区两处乔木覆盖区域：林地A与林地B，然后利用软件自带的工具描绘两处植被的大致边界，并保存为UTM坐标系的shp格式文件。
S2.获取时序植被指数数据。
在本实例中，通过USGS(https://glovis.usgs.gov/)获取所有覆盖山西省大同市南郊区的Landsat系列遥感影像文件合成的归一化植被指数(NDVI)数据，得到包括Landsat4/5 TM，Landsat7 ETM及Landsat8OLC影像的8天间隔NDVI数据集，总计829期数据，时间序列为1986年5月29日至2017年4月23日。在ENVI软件中通过批量处理插件按照S1中获取的植被边界对NDVI数据集进行裁剪，并统计和输出两片植被区域的NDVI均值，从而得到拟判别植被区域1986至2017年的时序NDVI数据集。
S3.去除时序植被指数数据噪声。
如图2所示，在本实例中，利用EXCEL对两组NDVI数据分别进行噪声消除，对时序数据最初以及最后的四个数据根据其获取时间所处的季节进行人工判断，若存在明显错误则将其剔除，否则使其预测值等于观测值。然后对其余所有数据依次与前后相邻的4个数据进行比较，若第i期观测值NOi小于相邻时点中最大值均值的50％时，则将第i期观测值存在较大偏差，并用相邻时点中最大值均值进行替换。具体的逻辑判断语言如下所示：
ELSE NPi＝NOi
由于遥感影像中NDVI的大小取决于地表在红光波段与近红外波段反射率的差异，绿色植被由于对红光具有较强的吸收能力从而有着较高的NDVI。当存在冰雪覆盖或云雾遮挡时，原有的植被区域将反射大部分的红光从而导致NDVI下降，因此大部分的NDVI数据噪声表现为观测值出现的异常下降。另一方面由于植被的生长具有明显的周期性与季节性规律，植被总是处于生长或衰落的过程中，其NDVI由此处于缓慢上升或下降中，某一时点的NDVI在理想条件下总是处于相邻两个时点NDVI的均值附近。因此在识别噪音与计算预测值时，应当将低于相邻时点NDVI均值的异常值作为噪音进行剔除，同时将其相邻时点NDVI均值作为预测值，从而避免直接剔除异常值带来的数据量较少。此外为了消除由于相邻时点的数据存在噪声带来的噪音识别与预测值计算误差，通过分别选择前后两个时点的最大值进行均值计算能够有效改善数据的准确性。
S4.重建植被生长曲线。
在本实例中，首先将剔除噪音后的829期NDVI数据根据其所在年份进行分类，并将其日期转化为年纪日形式(n/365)。由于不同年代、不同卫星在数据完整度上存在差异，导致不同年度获取的遥感数据也有所不同，本实例涵盖的时间序列中数据量最少的年份为1988年仅有7条，最多为2014年的43条。分别对各个年份的NDVI数据进行最小二乘法拟合，由于植被生长具有明显的周期性变化规律，且植被生长-凋落的周期为1年，因此建立正弦函数与余弦函数叠加的拟合曲线方程，其具体公式为：
如图3所示，黑色圆点为2016年度中噪声识别与预测后得到的NDVI预测值，将预测值及对应的年纪日时间带入拟合曲线方程求得方程各常数项后，得到该年度的NDVI最小二乘法(OLS)拟合方程(图3中曲线)。然后对获得各个年份的NDVI拟合曲线进行精度检验，对于拟合曲线的决定系数R2低于0.8的年份进行二次检验，重复S3、S4步骤，直至获得的拟合曲线精度达到要求。最终本实例中林地A与林地B的NDVI拟合曲线的决定系数R2均值分别为0.82与0.86，且各年份拟合取钱的决定系数均大于0.8。
S5.计算植被年度均值和生长趋势差异度。
在本实例中，首先分别利用林地A与林地B在不同年份的NDVI拟合曲线计算1987-2017年间各年度的NDVI均值。然后从第6年开始比较相邻的两个5年时间段内NDVI的变化趋势，首先分别对1987-1991年的5年间NDVI均值进行线性拟合，从而得到拟合方程的斜率与决定系数，用于反映两片植被在时间序列的前5年的生长趋势，斜率与决定系数分别用PT6与表示；然后对1992-1995年的5年间NDVI均值进行线性拟合，从而得到拟合方程的斜率与决定系数FT6与/>表示。然后将过去5年与未来5年的NDVI均值线性拟合方程斜率相除，并将两段拟合方程的决定系数均值作为权重，从而得到第6年前后植被生长趋势差异度D6，如下公式所示：
其中m为参与评价的年度序数。
如图4所示，其中黑色圆点为1999年至2008年各年份的NDVI均值，两条直线依据两个方程斜率的比值可以计算得到这10年间的植被生长趋势差异度。由于NDVI拟合曲线均值的线性拟合方程斜率存在正负差异，因此在计算生长趋势差异度Dm时加入逻辑判断语句。当PTm小于0时说明该区域植被在未来5年处于下降趋势，明显不属于重建植被，从而对直接对其Dm赋值为-1；当FTm小于0而PTm大于0是说明该区域植被从第m年起NDVI从过去5年的下降趋势转变为增长趋势，具有较大的可能为重建植被，其可能性取决于拟合方程的准确度，从而将其Dm值赋值为当FTm与PTm均大于0时，则按照上式计算。逻辑判断语句为：
IF PTm＜0
THEN Dm＝-1
ELSE IF FTm＜0
由此获得了整个时间序列中第6年时的植被生长趋势差异度，然后依次对剩余年份进行计算，从而获得林地A与林地B在1987-2017年间的植被生长趋势差异度，如表1所示。值得注意的是，由于1987-1991年与2014-2017年作为时间序列的始末年际段，所有这些年份的植被生长趋势将作为背景数据，不参与植被生长趋势差异度计算。
表1实施例中林地A与林地B与1987-2017年间植被生长趋势差异度
S6.通过阈值判别植被类型。
在本实例中，考虑到植被生长收到光照、降水、物候等多因素影像下，不同年度中NDVI存在一定的自然波动，因此将重建植被的Dm阈值设定为1.5，即未来5年中植被增长的速率较过去5年增长速率提高50％。当判别植被在任意5年时间段中出现2个年度的Dm值大于1.5时则认定这一区域植被为重建植被，且Dm值最大时的年份认定为重建植被的起始年份。分析数据可以发现，林地A在2004-2011年的8年间Dm值大于1.5的年份有5年，且最大值出现在2004年，而林地B仅在2004年Dm值大于1.5。因此认定林地A为重建植被区域，且植被重建开始于2004年，而林地B为自然植被区域。

Claims (7)

1.一种基于长时序植被指数的重建植被遥感判别方法，其特征在于包括以下步骤：
S1.利用卫星采集需要识别的研究区域的植被信息，生成具有明确边界的制备矢量图层遥感影像；
S2.获取研究区域的历年内所有时序植被指数遥感影像图，统计每期遥感影像图中拟识别植被区域的归一化植被指数NDVI均值，利用时序NDVI均值建立时序数据集；
S3.利用均值迭代法对研究区域的归一化植被指数进行噪声识别，去除遥感影像中冰雪覆盖及云雾遮挡的数据噪声，并利用相邻时点数据均值作为预测值替代从而剔除噪声；
S4.对剔除噪声后的时序NDVI数据基于傅里叶级数与普通最小二乘法进行拟合，获得研究区域植被在不同年份的生长曲线；
S5.基于植被生长曲线计算研究区域植被年度均值，并从不同年份的生长曲线中提取植被关键参数，包括5年NDVI均值线性拟合方程的斜率及决定系数，并基于植被关键参数提出植被生长趋势差异度以反映植被的持续变化情况；
S6.将植被生长趋势差异度作为阈值判别植被类型，从而将遥感影像图中人工种植的重建植被与演替发育的自然植被进行区分，并通过植被生长趋势差异度的最大值获得重建植被的种植时间。
2.根据权利要求1所述的基于长时序植被指数的重建植被遥感判别方法，其特征在于：所述步骤S1具体方法为：通过调查时点的卫星遥感影像或植被覆盖现状图，在ArcMap或Super Map地理信息软件中人工选择或者导入拟判别植被的边界范围，将选择的植被的边界范围保存为具有坐标系的矢量图层，植被区域的矢量图层为UTM坐标系。
3.根据权利要求1所述的基于长时序植被指数的重建植被遥感判别方法，其特征在于：所述步骤S2的具体方法为：利用Landsat系列遥感影像的归一化植被指数NDVI时序数据，包括Landsat4/5TM，Landsat7ETM及Landsat8 OLC影像的8天间隔NDVI数据集，在ArcMap或ENVI软件中通过批量裁剪功能提取并输出研究区域植被的NDVI数据均值，分别对获取的历年所有遥感影像数据进行统计后从而建立判别区域的时序植被指数数据集。
4.根据权利要求1所述的基于长时序植被指数的重建植被遥感判别方法，其特征在于：所述步骤S3的具体方法为：比较第i期NDVI数据的观测值NOi，i>2且为自然整数，与其前后相邻两期植被指数中最大值的均值进行比较，若小于其均值的50％则认为观测值为噪声，并用其均值作为预测值NPi进行代替，反之则保留其观测值，然后持续重复该步骤，直到消除所有的数据噪声，具体筛选规则为
IF
THEN
ELSE NPi＝NOi
其中NOi为第i期NDVI数据的观测值，NPi为第i期NDVI数据预测值，MAX()为识别输入数据中的最大值并将其输出的函数，NOi-2，NOi-1，NOi+1和NOi+2分别为第i-2，i-1，i+1，和i+2期NDVI数据的观测值。
5.根据权利要求1所述的基于长时序植被指数的重建植被遥感判别方法，其特征在于：所述步骤S4的具体方法为：将剔除噪音后的时序NDVI数据根据所在年份进行分类，并将其日期转化为年纪日形式，即n/365，然后分别对各个年份的NDVI数据通过傅里叶级数与最小二乘法进行拟合，其具体公式为：
其中T为函数的周期，根据植被的生长周期将其设置为365天，x为第i期NDVI数据点年纪日，a,b1,b2为待确定的拟合方程常数项；
通过带入获取的NDVI观测值确定拟合方程的常数项，从而得到各个年份的NDVI拟合曲线，通过拟合方程与原始观测值计算拟合曲线的决定系数R2，具体公式为：
其中yi为第i期遥感影像的NDVI观测值，fi为第i期NDVI的拟合值，为该年份内NDVI观测值的均值；
R2反映了拟合方程对NDVI的拟合效果，R2如果过低则说明方程无法反映植被的真实变化，反之则说明效果较好，对于拟合曲线的决定系数R2低于0.8的年份进行二次检验，重复S3、S4步骤，直至获得的拟合曲线精度达到要求。
6.根据权利要求5所述的基于长时序植被指数的重建植被遥感判别方法，其特征在于所述步骤S5的具体方法为：通过步骤S4中确定的不同年份的NDVI拟合曲线分别计算该年度的NDVI均值，
首先将年纪日形式的1至365天分别带入NDVI拟合曲线，求得该年份内所有的NDVI拟合值并计算该年份的NDVI均值NAj；
按照此方法依次求得所有年份的NDVI均值后，开始计算每5年NDVI均值的线性拟合方程并提取关键植被参数，包括5年NDVI均值线性拟合方程的斜率及决定系数：首先对第m年起当年至5年的NDVI拟合曲线均值进行线性拟合，其中5<m<n-5且为自然整数，n为NDVI时序数据的跨越年度，公式为：
Y＝FTmX+b1 Ⅲ
其中X为年份，Y为该年份的拟合曲线计算得到NDVI均值，FTm为方程斜率，b1为该线性拟合方程的常数项；
记录其线性拟合方程的斜率FTm以及决定系数FTR2，决定系数FTR2的计算公式Ⅳ如下所示：
其中NAj为第j年份的NDVI的均值，Yj为第j年份NDVI的线性拟合值，为参与该线性拟合方程的所有年份的NDVI均值；
类似的，对第m年之前的5年NDVI数据，进行线性拟合求得其拟合线性方程的斜率PTm与决定系数PTR2，公式为：
Y＝PTmX+b2 Ⅴ
其中X为年份，Y为拟合曲线计算得到该年份的NDVI均值，PTm为方程斜率，b2为该线性拟合方程的常数项；
其中NAj为第j年份的NDVI的均值，Yj为第j年份NDVI的线性拟合值，为参与该线性拟合方程的所有年份的NDVI均值；
在获得过去5年与未来5年的NDVI均值线性拟合方程斜率之后，开始结算植被生长趋势差异度；首先通过前后5年NDVI均值线性拟合方程斜率的比值，并将两段拟合方程的决定系数均值作为权重，计算得到第m年前后植被生长趋势差异度Dm，公式如下所示：
由于NDVI拟合曲线均值的线性拟合方程斜率存在正负差异，因此在计算生长趋势差异度Dm时加入逻辑判断语句：当PTm小于0时说明该区域植被在未来5年处于下降趋势，明显不属于重建植被，从而对直接对其Dm赋值为-1；当FTm小于0而PTm大于0是说明该区域植被从第m年起NDVI从过去5年的下降趋势转变为增长趋势，具有较大的可能为重建植被，其可能性取决于拟合方程的准确度，因此将其Dm值赋值为当FTm与PTm均大于0时，则按照公式Ⅶ计算；逻辑判断式为：
IF PTm<0
THEN Dm＝-1
ELSE IF FTm<0
THEN
ELSE
其中PTm为线性拟合公式Ⅴ的斜率，FTm为线性拟合公式Ⅲ的斜率，FTR2为拟合方程Ⅲ的决定系数，PTR2为拟合公式Ⅴ的决定系数，Dm为植被生长趋势差异度。
7.根据权利要求6所述的基于长时序植被指数的重建植被遥感判别方法，其特征在于所述步骤S6的具体方法为：考虑到植被生长受到光照、降水、物候多因素影响下，不同年度中NDVI存在一定的自然波动，因此将判别是否为重建植被的Dm阈值设定为1.5，即在未来5年植被增长的斜率较过去5年增长斜率提高50％，当判别区域植被在任意5年时间段中出现2个年度的Dm大于1.5时则认定这一区域植被为重建植被；反之，当Dm小于1.5时则认定这一区域植被不是重建植被，同时，Dm最大值出现的年份则被认为是植被重建的起始年份。
CN202010701602.5A 2020-07-20 2020-07-20 一种基于长时序植被指数的重建植被遥感判别方法 Active CN111832506B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010701602.5A CN111832506B (zh) 2020-07-20 2020-07-20 一种基于长时序植被指数的重建植被遥感判别方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010701602.5A CN111832506B (zh) 2020-07-20 2020-07-20 一种基于长时序植被指数的重建植被遥感判别方法

Publications (2)

Publication Number Publication Date
CN111832506A CN111832506A (zh) 2020-10-27
CN111832506B true CN111832506B (zh) 2023-10-17

Family

ID=72924435

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010701602.5A Active CN111832506B (zh) 2020-07-20 2020-07-20 一种基于长时序植被指数的重建植被遥感判别方法

Country Status (1)

Country Link
CN (1) CN111832506B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112861810B (zh) * 2021-03-23 2021-11-02 中国科学院、水利部成都山地灾害与环境研究所 一种基于时序遥感观测数据的人工林种植时间自动检测方法
CN113516077B (zh) * 2021-07-14 2023-05-26 生态环境部卫星环境应用中心 地表环境变化监测方法、装置、可读存储介质和设备
CN114092831B (zh) * 2021-12-02 2023-03-24 中国科学院东北地理与农业生态研究所 一种草本沼泽植被物候信息提取方法
CN116305875B (zh) * 2023-02-17 2023-08-29 中国科学院地理科学与资源研究所 数值天气预报模式中叶面积指数的数据处理方法及装置
CN117407838B (zh) * 2023-11-07 2024-07-12 生态环境部华南环境科学研究所(生态环境部生态环境应急研究所) 一种绿色植被适应气候变化恢复力确定方法及系统
CN117789024B (zh) * 2023-12-26 2024-07-23 成都理工大学 一种基于机器学习的黑土滩演变过程遥感监测方法
CN117689959B (zh) * 2024-01-30 2024-06-14 中交第二公路勘察设计研究院有限公司 一种融合植被生命周期特征的遥感分类方法
CN117688505B (zh) * 2024-02-04 2024-04-19 河海大学 一种植被大范围区域化负异常的预测方法及系统

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106022224A (zh) * 2016-05-12 2016-10-12 北京师范大学 一种冬小麦识别方法
CN106803059A (zh) * 2016-12-02 2017-06-06 浙江工业大学 一种遥感植被指数时间序列森林监测方法
CN108982369A (zh) * 2018-04-28 2018-12-11 中国农业大学 融合gf-1wfv和modis数据的地块尺度作物长势监测方法
CN109840516A (zh) * 2019-03-06 2019-06-04 福州大学 一种基于时序遥感影像的水体变化自动识别方法
CN110169336A (zh) * 2019-06-21 2019-08-27 昆山千亿圆生物科技有限公司 一种大棚种植灌溉系统及方法
CN110598514A (zh) * 2019-06-10 2019-12-20 南京大学 一种土地整治项目区地块尺度农作物播种面积监测方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106022224A (zh) * 2016-05-12 2016-10-12 北京师范大学 一种冬小麦识别方法
CN106803059A (zh) * 2016-12-02 2017-06-06 浙江工业大学 一种遥感植被指数时间序列森林监测方法
CN108982369A (zh) * 2018-04-28 2018-12-11 中国农业大学 融合gf-1wfv和modis数据的地块尺度作物长势监测方法
CN109840516A (zh) * 2019-03-06 2019-06-04 福州大学 一种基于时序遥感影像的水体变化自动识别方法
CN110598514A (zh) * 2019-06-10 2019-12-20 南京大学 一种土地整治项目区地块尺度农作物播种面积监测方法
CN110169336A (zh) * 2019-06-21 2019-08-27 昆山千亿圆生物科技有限公司 一种大棚种植灌溉系统及方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
NDVI-Based Long-Term Vegetation Dynamics and Its Response to Climatic Change in the Mongolian Plateau;Gang Bao等;《Remote Sens.》;第6卷(第9期);8337-8358 *
The long-term effects of underground mining on the growth of tree, shrub, and herb communities in arid and semiarid areas in China;Jiaxin Mi等;《Land Degradation and Development》;第32卷(第3期);1-28 *
半干旱区井工矿山地表形变对植被的长期影响研究;米家鑫;《中国博士学位论文全文数据库_基础科学辑》;A006-278 *
基于多源遥感数据的竹林LAI多尺度同化及在碳循环模拟中的应用;邢璐琪;《中国优秀硕士学位论文全文数据库_农业科技辑》;D049-111 *
青藏高原典型植被生长季遥感模型提取分析;常清等;《地球信息科学学报》;第16卷(第05期);815-823 *

Also Published As

Publication number Publication date
CN111832506A (zh) 2020-10-27

Similar Documents

Publication Publication Date Title
CN111832506B (zh) 一种基于长时序植被指数的重建植被遥感判别方法
CN112070234B (zh) 复杂场景下水体叶绿素和藻蓝素陆基遥感机器学习算法
CN108595414B (zh) 基于源汇空间变量推理的土壤重金属企业污染源识别方法
CN110008948B (zh) 基于变分自编码网络的高光谱图像目标检测方法
CN112395808A (zh) 一种结合随机森林和协同克里金的生物量遥感制图方法
CN112464920B (zh) 基于极端随机树的fy-3d红外高光谱云检测方法
CN108896021B (zh) 基于航空摄影测量点云提取人工林林分结构参数的方法
CN115437036A (zh) 一种基于葵花卫星的对流初生预报方法
CN110765962A (zh) 一种基于三维点云轮廓分维值的植物识别分类方法
CN110263412B (zh) 副热带急流或极锋急流强度和径向位置协同变化的表征方法
CN110888186A (zh) 基于gbdt+lr模型的冰雹和短时强降水预报方法
CN114120137B (zh) 一种基于时序植被遥感影像的湿地要素时空演变监测方法
CN115512223B (zh) 一种融合多种变化检测算法的红树林动态监测方法
CN114048944A (zh) 一种暴雨诱发地质灾害下应撤离人口及损毁房屋的预估方法
CN113281716A (zh) 一种光子计数激光雷达数据去噪方法
CN115907574A (zh) 一种流域性暴雨洪涝承灾体重置成本遥感模拟方法
CN113361782A (zh) 基于改进mkpls的光伏发电功率短期滚动预测方法
CN117035174A (zh) 一种木麻黄单木地上生物量的估算方法与系统
CN117132894A (zh) 一种基于时序遥感影像的露天煤矿区生态损伤区识别方法
CN112033937A (zh) 一种水体提取精度的评价方法
CN112380984B (zh) 一种基于遥感的盐沼植被缓流能力空间评价方法
CN110208876B (zh) 副热带急流和极锋急流径向位置协同变化的表征方法
CN114971041A (zh) 一种基于残差网络的海冰面积预测方法及系统
CN114722909A (zh) 一种基于低维卷积神经网络的太阳耀斑时间序列分类方法
CN113128523A (zh) 一种基于时间序列遥感影像自动提取珊瑚礁的方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB02 Change of applicant information

Address after: 221116 Research Institute of China University of Mining and Technology, 1 University Road, Xuzhou, Jiangsu

Applicant after: CHINA University OF MINING AND TECHNOLOGY

Applicant after: Datong Coal Mine Group Co.,Ltd.

Address before: 037003 xinpingwang mining area, Datong City, Shanxi Province

Applicant before: Datong Coal Mine Group Co.,Ltd.

Applicant before: CHINA University OF MINING AND TECHNOLOGY

CB02 Change of applicant information
CB03 Change of inventor or designer information

Inventor after: Mi Jiaxin

Inventor after: Yang Yongjun

Inventor after: Ma Zhanyuan

Inventor after: Zhao Chende

Inventor after: Zhang Shaoliang

Inventor after: Chang Xiaohua

Inventor after: Lu Xi

Inventor after: Wang Yuming

Inventor after: Wei Chuanbo

Inventor after: Yue Xuelian

Inventor before: Ma Zhanyuan

Inventor before: Yang Yongjun

Inventor before: Mi Jiaxin

Inventor before: Zhao Chende

Inventor before: Zhang Shaoliang

Inventor before: Chang Xiaohua

Inventor before: Lu Xi

Inventor before: Wang Yuming

Inventor before: Wei Chuanbo

Inventor before: Yue Xuelian

CB03 Change of inventor or designer information
GR01 Patent grant
GR01 Patent grant