CN116682003A - 一种植被异常的遥感即时探测方法和装置 - Google Patents

一种植被异常的遥感即时探测方法和装置 Download PDF

Info

Publication number
CN116682003A
CN116682003A CN202310524502.3A CN202310524502A CN116682003A CN 116682003 A CN116682003 A CN 116682003A CN 202310524502 A CN202310524502 A CN 202310524502A CN 116682003 A CN116682003 A CN 116682003A
Authority
CN
China
Prior art keywords
vegetation
month
pixel
historical
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.)
Pending
Application number
CN202310524502.3A
Other languages
English (en)
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.)
Beijing Normal University
Original Assignee
Beijing Normal University
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 Beijing Normal University filed Critical Beijing Normal University
Priority to CN202310524502.3A priority Critical patent/CN116682003A/zh
Publication of CN116682003A publication Critical patent/CN116682003A/zh
Pending legal-status Critical Current

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
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/10—Image acquisition
    • G06V10/12—Details of acquisition arrangements; Constructional details thereof
    • G06V10/14—Optical characteristics of the device performing the acquisition or on the illumination arrangements
    • G06V10/143—Sensing or illuminating at different wavelengths
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/764—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using classification, e.g. of video objects
    • G06V10/765—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using classification, e.g. of video objects using rules for classification or partitioning the feature space
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/766—Arrangements for image or video recognition or understanding using pattern recognition or machine learning using regression, e.g. by projecting features on hyperplanes

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Multimedia (AREA)
  • General Health & Medical Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Computing Systems (AREA)
  • Databases & Information Systems (AREA)
  • Evolutionary Computation (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Image Processing (AREA)

Abstract

本发明公开一种植被异常的遥感即时探测方法和装置，包括：将目标区域在研究时间段内的地表反射率影像划分成历史时期和监测期分别进行预处理，将历史时期中各年各月的所有数据合成为各年逐月数据；对预处理后的各年逐月数据和待探测影像逐像元计算多种植被指数；对各像元的各年逐月数据的植被指数逐月份采用两次线性回归依次去除噪声和植被生长趋势信息；将各像元的植被指数在各时间点的历史残差作为样本数据进行核密度估计以构建各像元各月的植被正常模式，进而根据待探测影像所属月份找到对应月的植被正常模式，结合各像元的植被指数在对应月的预测值和植被正常模式的边界，逐像元判别植被是否发生异常。本发明可提高植被异常即时探测的普适性。

Description

一种植被异常的遥感即时探测方法和装置
技术领域
本发明涉及遥感图像处理技术领域，具体涉及一种植被异常的遥感即时探测方法和装置，可应用于植被遥感的动态变化监测研究中。
背景技术
森林等植被作为地表重要组分之一，是地球碳水循环和能量交换的重要场所，对维持全球生态平衡具有不可替代的作用。当前，日益活跃的人类活动或自然因素如砍伐、火灾、病虫害等时有发生，其引起的植被损失属于植被异常。对植被异常进行即时探测，及时掌握植被覆盖变化、介入不可持续的人类活动，有助于维护生态系统的稳定。
植被异常探测归属于变化检测，遥感数据因其高频、大覆盖的优势是变化检测的主流工具，且越来越多的遥感历史影像在探测和量化地表覆盖变化方面发挥着重要作用。近十年来，利用陆地卫星Landsat时间序列数据，发展了众多参数、非参数植被异常探测算法，如基于陆地卫星的干扰和再生趋势检测算法LandTrendr(Landsat-based detectionof Trends in Disturbance and Recovery)、季节和趋势分量断点检测算法BFAST(BreaksFor Additive Seasonal and Trend)、土地覆盖连续变化检测和分类模型CCDC(Continuous Change Detection and Classification)、植被异常变化检测算法AVOCADO(Anomaly Vegetation Change Detection)等，但这些算法主要针对历史异常的探测，即时性不强。
现今为了更快捕捉植被变化，即时探测需求在增加。大多数现有即时探测方法仅针对热带地区森林砍伐所发展，大概可分为两类：第一类基于条件概率，主要思想为将某植被特征值分为过去、现在和未来，对每个时刻特征值计算其为非森林的条件概率，若t时刻的概率大于0.5则认为该时刻为潜在异常点，进而使用t-1、t以及未来的t+i时刻的观测值迭代贝叶斯更新计算，以确认或拒绝t时刻处的异常事件；第二类基于植被物候参数模型，主要思想为利用某植被特征的历史时间序列数据，根据正弦、余弦或双逻辑函数等，拟合正常物候并预测待探测t时刻的观测值，将实际观测与预测值之间的偏差作为扰动指示。
然而，现有即时探测方法应用场景较为单一，还未有应用于不同植被异常类型的普适性方法；且物候参数模型算法数学理论较为复杂，探测结果受参数拟合效果影响；另外，在算法使用的特征方面，大多数方法仅使用某一种植被指数，可能会忽略掉其他波段的信息。
发明内容
有鉴于此，本发明的主要目的在于提供了一种植被异常的遥感即时探测方法和装置，用于提高植被异常即时探测的普适性。
依据本发明的第一方面，提供了一种植被异常的遥感即时探测方法，包括：
收集目标区域在研究时间段内所有可用的Landsat系列的地表反射率影像；
将收集的地表反射率影像划分成历史时期和监测期，其中历史时期数据用于构建各月的植被正常模式，监测期数据用于植被异常的即时探测；
分别对历史时期数据和监测期中的待探测影像进行预处理，并将历史时期中各年各月的所有数据通过均值合成方式处理为历史各年逐月数据；
分别对预处理后的历史各年逐月数据和待探测影像逐像元计算多种植被指数；
对各像元的历史各年逐月数据的植被指数逐月份采用两次时间序列线性回归依次去除噪声和植被生长趋势信息，并根据第二次线性回归结果得到各像元的植被指数在各历史时间点的实际值与回归值之差的历史残差，以及得到各像元的植被指数在未来对应月的预测值；
将各像元的植被指数的历史残差作为样本数据进行非参数的核密度估计以构建各像元各月的植被正常模式，并根据设定概率阈值确定出各像元各月的植被正常模式的边界；
根据待探测影像所属月份找到对应月的植被正常模式，结合各像元的植被指数在未来对应月的预测值和植被正常模式的边界，逐像元判别目标区域的植被是否发生异常。
依据本发明的第二方面，提供了一种植被异常的遥感即时探测装置，包括：
影像收集单元，用于收集目标区域在研究时间段内所有可用的Landsat系列的地表反射率影像；
时期划分单元，用于将收集的地表反射率影像划分成历史时期和监测期，其中历史时期数据用于构建各月的植被正常模式，监测期数据用于植被异常的即时探测；
预处理单元，用于分别对历史时期数据和监测期中的待探测影像进行预处理，并将历史时期中各年各月的所有数据通过均值合成方式处理为历史各年逐月数据；
植被指数计算单元，用于分别对预处理后的历史各年逐月数据和待探测影像逐像元计算多种植被指数；
时间序列处理单元，用于对各像元的历史各年逐月数据的植被指数逐月份采用两次时间序列线性回归依次去除噪声和植被生长趋势信息，并根据第二次线性回归结果得到各像元的植被指数在各历史时间点的实际值与回归值之差的历史残差，以及得到各像元的植被指数在未来对应月的预测值；
正常模式构建单元，用于将各像元的植被指数的历史残差作为样本数据进行非参数的核密度估计以构建各像元各月的植被正常模式，并根据设定概率阈值确定出各像元各月的植被正常模式的边界；
异常判别单元，用于根据待探测影像所属月份找到对应月的植被正常模式，结合各像元的植被指数在未来对应月的预测值和植被正常模式的边界，逐像元判别目标区域的植被是否发生异常。
依据本发明的第三方面，提供了一种电子设备，包括：处理器和存储计算机可执行指令的存储器，所述可执行指令在被所述处理器执行时，实现前述的植被异常的遥感即时探测方法。
依据本发明的第四方面，提供了一种计算机可读存储介质，所述计算机可读存储介质存储一个或多个程序，所述一个或多个程序当被处理器执行时，实现前述的植被异常的遥感即时探测方法。
本发明的有益效果是：
本发明面向各植被典型异常类型，同时使用从不同侧度反映植被状态的多种植被指数，结合非参数核密度估计，可以普遍适用于不同区域不同植被异常的探测；而且本发明的整个异常探测过程逐像元进行，按月构建植被正常模式以及逐像元处理可以进一步提高本发明在不同时间点以及影像上不同空间位置处进行异常探测的普适性，从而不限区域和植被类型，一旦有可用影像时，就能够实现植被异常的即时探测。
附图说明
通过阅读下文优选实施方式的详细描述，各种所有的优点和益处对于本领域普通技术人员将变得清楚明了。附图仅用于示出优选实施方式的目的，而并不认为是对本发明的限制。而且在整个附图中，用相同的参考符号表示相同的部件。在附图中：
图1示出了本发明一个实施例的植被异常的遥感即时探测方法的流程示意图；
图2为对应图1所示方法的植被异常的遥感即时探测的流程框架图；
图3示出了本发明一个实施例的NDVI、NBR指数两次线性回归示意图；
图4示出了本发明一个实施例的核密度估计后的概率密度示意图；
图5示出了本发明一个实施例的像元异常判别示意图；
图6示出了本发明一个实施例的植被异常的遥感即时探测装置的框图；
图7示出了本发明一个实施例的电子设备的框图；
图8示出了本发明实施例1的哥伦比亚亚马孙森林的Landsat影像数据；
图9示出了本发明实施例1的各待即时探测日期的短波红外2波段图及对应的即时探测结果图；
图10示出了本发明实施例2的内蒙古根河镇大兴安岭的Landsat影像数据；
图11示出了本发明实施例2的各待即时探测日期的近红外波段图及对应的即时探测结果图。
具体实施方式
下面将参照附图更详细地描述本发明的示例性实施例。提供这些实施例是为了能够更透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。虽然附图中显示了本发明的示例性实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。
图1示出了本发明一个实施例的植被异常的遥感即时探测方法的流程示意图，图2为对应图1所示方法的植被异常的遥感即时探测的流程框架图。
参见图1和图2，本发明实施例的方法包括如下步骤S110至步骤S160：
步骤S110，收集目标区域在研究时间段内所有可用的Landsat系列的地表反射率影像。
Landsat(陆地资源卫星系列)数据因其有长期连续观测记录、较高空间分辨率、免费共享等优点，是研究不同类型土地覆盖变化的重要数据源之一。本发明对目标区域使用研究时间段内所有可用的Landsat系列collection2 level2地表反射率产品，卫星涉及Landsat4、Landsat5、Landsat7、Landsat8，可在https://earthexplorer.usgs.gov/下载。
步骤S120，将收集的地表反射率影像划分成历史时期和监测期，其中历史时期数据用于构建各月的植被正常模式，监测期数据用于植被异常的即时探测。
将收集的Landsat系列的地表反射率影像的时间序列数据分为历史时期和监测期，历史时期中的数据应较稳定，没有发生明显异常，用于构建各月的植被正常模式；监测期通常始于历史时期结束的下一年，监测期中的数据用于植被异常的即时探测。
步骤S130，分别对历史时期数据和监测期中的待探测影像进行预处理，并将历史时期中各年各月的所有数据通过均值合成方式处理为历史各年逐月数据。
由于收集的地表反射率影像数据涉及多个卫星传感器，波段设置不同，且同一区域的每景影像在空间位置上存在一定差异，故预处理步骤中需对所有数据进行波段统一和空间配准，处理后的数据共有7个波段：蓝、绿、红、近红外、短波红外1、短波红外2和QA(Quality Assessment)质量控制波段，并拥有相同的空间范围。
另外，云、云影、积雪等会影响地表真实反射率，故需结合地表反射率产品提供的QA波段将其掩膜，使其不参与后续处理。关于QA波段的具体信息可查询
https://www.usgs.gov/landsat-missions/landsat-collection-2-surface- reflectance中的产品手册。
历史时期数据和监测期中的待探测影像均需要进行波段统一、空间配准以及云、云影、积雪等去除的预处理，除此之外对历史时期中各年各月的所有数据还进一步通过均值合成方式处理为历史各年逐月数据，便于后续构建各月的植被正常模式。
步骤S140，分别对预处理后的历史各年逐月数据和待探测影像逐像元计算多种植被指数。
本步骤S140分别对预处理后的历史各年逐月数据和待探测影像逐像元计算的植被指数包括但不限于：归一化差异植被指数NDVI(Normalized Difference VegetationIndex)和归一化燃烧比指数NBR(Normalized Burn Ratio)；其中，
NDVI指数放大了近红外波段与红波段的反射率差异，可表征植被覆盖程度；NBR指数放大了近红外波段与短波红外波段的反射率差异，可表征植被叶片缺水程度。这两个指数分别从植被绿度和含水量两个方面反映植被状态，具体计算公式如下：
其中的RNIR、Rred、RSWIR2分别为近红外、红以及短波红外2波段反射率。
可以理解，除了逐像元计算NDVI和NBR两种植被指数外，本发明还可以根据需要计算其他植被指数，例如增强植被指数EVI(Enhanced Vegetation Index)、水分胁迫指数MSI(Moisture Stress Index)等。
步骤S150，对各像元的历史各年逐月数据的植被指数逐月份采用两次时间序列线性回归依次去除噪声和植被生长趋势信息，并根据第二次线性回归结果得到各像元的植被指数在各历史时间点的实际值与回归值之差的历史残差，以及得到各像元的植被指数在未来对应月的预测值。
本步骤S150主要针对历史时期数据，也即对预处理后的历史各年逐月数据，逐月份采用两次线性回归对噪声和植被生长趋势信息依次进行去除。
由于QA波段掩膜并不彻底，影像上还会存在残余的云、云影等噪声，为后续植被正常模式构建造成误差，故需去除。另外，植被正常模式构建基于的是时间序列数据，包含有植被生长趋势信息，同样需要去除，可见本步骤S150也能在一定程度上避免不同植被类型的生长差异，从而增加本发明的普适性。
云信息在NDVI指数的时间序列上常表现为一个向下的“尖峰”，针对每个像元各月份数据，利用NDVI指数进行第一次时间序列线性回归，若某个历史时间点的实际值低于回归值，且两者之间的差异大于某阈值，则认为该点为噪声点需要去除。参与异常探测的NDVI、NBR指数在保留下来的历史时间点上进行第二次时间序列线性回归，通过将各历史时间点的实际值与回归值相减的方式可以去除植被生长趋势信息。
因此本步骤S150中，对各像元的历史各年逐月数据的植被指数逐月份采用两次时间序列线性回归依次去除噪声和植被生长趋势信息，具体可以包括：
针对每个像元的历史各年逐月数据，利用NDVI指数进行第一次时间序列线性回归，若某个时间点实际值低于回归值，且两者之间的差异大于设定剔除阈值，则将该时间点作为噪声点进行剔除，同时NBR指数中对应时间点的数据也被剔除；
对NDVI指数和NBR指数保留下来的时间点进行第二次时间序列线性回归，得到各时间点的实际值与回归值之差以去除植被生长趋势信息。
图3示出了本发明一个实施例的NDVI、NBR指数两次线性回归示意图。图3中“1th_regression”、“2th_regression”依次表示第一次线性回归和第二次线性回归，“data tobe removed”表示被剔除的数据，“final data”表示保留下来的数据。在图3所示的实施例中，第一次线性回归设定的噪声点剔除阈值为0.1，对NDVI指数进行第一次时间序列线性回归，不满足该阈值的数据认为是残留云信息，被剔除(图3上图中的圆点)，同时NBR指数中对应位置的数据也被剔除(图3下图中的圆点)；对NDVI和NBR指数保留下来的数据(图3上下图中的菱形点)进行第二次时间序列线性回归以去除趋势信息。
根据第二次线性回归结果可以得到各像元的植被指数在各历史时间点的实际值与回归值之差的历史残差，该历史残差进一步用于构建植被正常模式。
根据第二次线性回归结果还可以得到各像元的植被指数在未来对应月的预测值。预测是连接历史时期和监测期的桥梁。第二次线性回归结果提供了各像元植被指数在某月的变化趋势，据此趋势可预测得到待探测影像所处月份的指数预测值，进而得到与对应月的待探测影像的指数实际值间的预测值残差，该预测值残差用于后续判别植被是否发生异常。
步骤S160，将各像元的植被指数的历史残差作为样本数据进行非参数的核密度估计以构建各像元各月的植被正常模式，并根据设定概率阈值确定出各像元各月的植被正常模式的边界。
本步骤S160将上一步骤得到的各像元的植被指数的历史残差作为样本数据进行核密度估计，以逼近每个像元各个月份数据的正常分布模式。核密度估计属于非参数估计方法，完全利用数据本身信息，避免人为主观带入的先验知识，根据样本集合来拟合总体数据的分布密度函数，从而相对于参数估计能够对样本数据进行最大程度的近似。
核密度估计可以调用python sklearn KernelDensity代码包实现。其原理是将一个带宽为h的核函数居中放置在每个观测点上，并对所有核函数取平均值，直到得到最终的概率密度估计。具体地，
定义X＝(X1，X2，......，Xn)T为样本数据，每个样本包含了参与构建的N个植被指数的历史残差：
Xi＝(Index1i，Index2i，......，IndexNi)，i＝1，......，n
则N维高斯核密度估计表达式为：
其中：x为N维植被指数空间中的点，h为核函数带宽，K()代表高斯核函数，||x-Xi||代表N维空间中x与样本点的距离。
需要说明的是，N为参与正常模式构建的植被指数数量，如果N为1，那么核密度估计就是一维的，如果为2，则是二维的，以此类推。每个植被指数相当于一个维度，在实际计算时根据N维核密度估计的公式一起计算。并且由于各个植被指数的历史残差都可认为是近似的高斯分布，本发明将各核函数按一致的情况进行计算，N维核密度估计的表达式也都是高斯核函数。
带宽的确定方法有经验公式法、交叉验证法、嵌入式方法(plug-in method)等，其受样本数量以及核函数维度的影响。在核密度估计为N维时，可以采用经验公式h＝1/(N+3)n(-1/(N+4))来确定核函数带宽。该经验公式法相比于其他带宽确定方法更为简单，应用起来方便快速，且该经验公式的参数系数已经过反复实验，核密度估计效果较好。
图4示出了本发明一个实施例的核密度估计后的概率密度示意图。其中，图4(a)为二维核密度估计后的概率密度示意图，图4(b)为图4(a)在x-o-y平面上的投影。图4中“NDVI_Diff”、“NBR_Diff”分别表示经过两次线性回归后，NDVI、NBR指数的历史残差，“Kernel Density Estimation”表示核密度估计。显然通过核密度估计，数据集中分布区域的概率密度较高。若以一个平行于x-o-y的平面截取概率密度曲面，那么截面以上区域的积分可用来表示正常的概率，据此根据一定概率阈值能提取出植被正常模式分布的大致范围，此时植被正常模式的边界为截面对应的概率密度。如下式所示：
其中：为植被正常模式的边界，即截面对应的概率密度，可由/>的/>分位数计算得到，/>为截面以上所包含的区域。由上式可知，该区域至少包含正常的概率为一旦设定概率阈值即可确定/>从而确定植被正常模式分布的边界。例如当概率阈值取95％时，代表至少包含了95％的正常，超出该范围的点为正常的概率小于5％，此时边界为核密度估计结果5％分位数处的概率密度。
步骤S170，根据待探测影像所属月份找到对应月的植被正常模式，结合各像元的植被指数在未来对应月的预测值和植被正常模式的边界，逐像元判别目标区域的植被是否发生异常。
步骤S150中得到了各像元的植被指数在未来对应月的预测值，步骤S160中通过核密度估计构建起了植被正常模式与其边界，本步骤S170是逐像元将待探测影像的植被指数的实际值与对应月的预测值相减得到预测值残差，然后将预测值残差带入对应月的植被正常模式，可以得到相应的概率密度，通过与对应边界的概率密度相对比来判别目标区域内的植被是否发生异常。
本步骤S170具体包括：根据待探测影像所属月份找到对应月的植被正常模式，逐像元将待探测影像的植被指数的实际值与对应月的预测值之差的预测值残差带入对应月的植被正常模式，将得到的概率密度与对应边界的概率密度相对比，若待探测影像的某个像元的概率密度小于对应边界的概率密度，则判别与该像元对应区域的植被异常，否则判别与该像元对应区域的植被正常。
图5示出了本发明一个实施例的像元异常判别示意图。其中图5(a)、(b)为不同待判别像元的异常判别示意图，“samples”圆点为参与核密度估计的样本点，“new data”叉形点为待判别是否发生异常的点。图5(a)、(b)中若干等高线从外至内分别代表至少包含的正常概率为95％、90％、75％和50％时所对应的正常边界。显然，若将95％作为概率阈值，那么图5(a)中叉形点识别为正常，而图5(b)中叉形点识别为异常。
综上所述，本发明通过收集目标区域在研究时间段内所有可用的Landsat系列的地表反射率影像，将收集的地表反射率影像划分成历史时期和监测期，其中历史时期数据以月为单位，用于非参数核密度估计构建各月的植被正常模式，监测期数据用于植被异常的即时探测；然后分别对历史时期数据和监测期中的待探测影像进行预处理，并将历史时期中各年各月的所有数据通过均值合成方式处理为历史各年逐月数据，接着分别对预处理后的历史各年逐月数据和待探测影像逐像元计算多种植被指数，这些植被指数可从不同侧度反映植被状态；之后对各像元的历史各年逐月数据的植被指数采用两次时间序列线性回归依次去除噪声和植被生长趋势信息，并根据第二次线性回归结果得到各像元的植被指数在各历史时间点的实际值与回归值之差的历史残差，以及得到各像元的植被指数在未来对应月的预测值，进而将各像元的植被指数的历史残差作为样本数据进行非参数的核密度估计以构建各像元各月的植被正常模式，并根据设定概率阈值确定出各像元各月的植被正常模式的边界，最后根据待探测影像所属月份找到对应月的植被正常模式，结合各像元的植被指数在未来对应月的预测值和植被正常模式的边界，逐像元判别目标区域的植被是否发生异常。
可见，本发明面向各植被典型异常类型，同时使用从不同侧度反映植被状态的多种植被指数，结合非参数核密度估计，可以普遍适用于不同区域不同植被异常的探测；而且本发明的整个异常探测过程逐像元进行，按月构建植被正常模式以及逐像元处理可以进一步提高本发明在不同时间点以及影像上不同空间位置处进行异常探测的普适性，从而实现不限区域和植被类型，一旦有可用影像时，就能够进行植被异常的即时探测。
与前述的植被异常的遥感即时探测方法同属于一个技术构思，本发明还提供了一种植被异常的遥感即时探测装置。图6示出了本发明一个实施例的植被异常的遥感即时探测装置的框图，参见图6，本发明的植被异常的遥感即时探测装置600包括：
影像收集单元610，用于收集目标区域在研究时间段内所有可用的Landsat系列的地表反射率影像；
时期划分单元620，用于将收集的地表反射率影像划分成历史时期和监测期，其中历史时期数据用于构建各月的植被正常模式，监测期数据用于植被异常的即时探测；
预处理单元630，用于分别对历史时期数据和监测期中的待探测影像进行预处理，并将历史时期中各年各月的所有数据通过均值合成方式处理为历史各年逐月数据；
植被指数计算单元640，用于分别对预处理后的历史各年逐月数据和待探测影像逐像元计算多种植被指数；
时间序列处理单元650，用于对各像元的历史各年逐月数据的植被指数逐月份采用两次时间序列线性回归依次去除噪声和植被生长趋势信息，并根据第二次线性回归结果得到各像元的植被指数在各历史时间点的实际值与回归值之差的历史残差，以及得到各像元的植被指数在未来对应月的预测值；
正常模式构建单元660，用于将各像元的植被指数的历史残差作为样本数据进行非参数的核密度估计以构建各像元各月的植被正常模式，并根据设定概率阈值确定出各像元各月的植被正常模式的边界；
异常判别单元670，用于根据待探测影像所属月份找到对应月的植被正常模式，结合各像元的植被指数在未来对应月的预测值和植被正常模式的边界，逐像元判别目标区域的植被是否发生异常。
在本发明的一个实施例中，植被指数计算单元640分别对预处理后的历史各年逐月数据和待探测影像逐像元计算的植被指数包括但不限于：归一化差异植被指数NDVI和归一化燃烧比指数NBR；其中，
NDVI指数放大了近红外波段与红波段的反射率差异，可表征植被覆盖程度，其计算公式为：
NBR指数放大了近红外波段与短波红外波段的反射率差异，可表征植被叶片缺水程度，其计算公式为：
其中的RNIR、Rred、RSWIR2分别为近红外、红以及短波红外2波段的反射率。
在本发明的一个实施例中，时间序列处理单元650对各像元的历史各年逐月数据的植被指数逐月份采用两次时间序列线性回归依次去除噪声和植被生长趋势信息，具体包括：
针对每个像元的历史各年逐月数据，利用NDVI指数进行第一次时间序列线性回归，若某个时间点实际值低于回归值，且两者之间的差异大于设定剔除阈值，则将该时间点作为噪声点进行剔除，同时NBR指数中对应时间点的数据也被剔除；
对NDVI指数和NBR指数保留下来的时间点进行第二次时间序列线性回归，得到各时间点的实际值与回归值之差以去除植被生长趋势信息。
在本发明的一个实施例中，正常模式构建单元660将各像元的植被指数的历史残差作为样本数据进行非参数的核密度估计以构建各像元各月的植被正常模式，具体包括：
定义X＝(X1，X2，......，Xn)T为样本数据，每个样本包含了参与构建的N个植被指数的历史残差：
Xi＝(Index1i，Index2i，......，IndexNi)，i＝1，......，n
则N维核密度估计表达式为：
其中：x为N维植被指数空间中的点，h为核函数带宽，K()代表高斯核函数，||x-Xi||代表N维空间中x与样本点的距离。
其中，在核密度估计为N维时，采用经验公式h＝1/(N+3)n(-1/(N+4))来确定核函数带宽。
在本发明的一个实施例中，正常模式构建单元660根据设定概率阈值确定出各像元各月的植被正常模式的边界，具体包括：
利用下式确定出各像元各月的植被正常模式的边界：
其中：为植被正常模式的边界，即截面对应的概率密度，可由/>的/>分位数计算得到，/>为截面以上所包含的区域；由上式可知，该区域至少包含正常的概率为一旦设定概率阈值即可确定/>从而确定植被正常模式分布的边界。
在本发明的一个实施例中，异常判别单元670具体用于：
根据待探测影像所属月份找到对应月的植被正常模式，逐像元将待探测影像的植被指数的实际值与对应月的预测值之差的预测值残差带入对应月的植被正常模式，将得到的概率密度与对应边界的概率密度相对比，若待探测影像的某个像元的概率密度小于对应边界的概率密度，则判别与该像元对应区域的植被异常，否则判别与该像元对应区域的植被正常。
与前述的植被异常的遥感即时探测方法同属于一个技术构思，本发明还提供了一种电子设备。图7示出了本发明一个实施例的电子设备的框图，参见图7，在硬件层面，本发明的电子设备包括：处理器和存储器。可选地还包括接口模块、通信模块等。存储器可以包括内存，例如高速随机存取存储器(Random-Access Memory，RAM)，还可以包括非易失性存储器(non-volatile memory)，例如至少一个磁盘存储器等。当然，该电子设备还可能包括其他业务所需要的硬件。
处理器、接口模块、通信模块和存储器可以通过内部总线相互连接，该内部总线可以是ISA(Industry Standard Architecture，工业标准体系结构)总线、PCI(PeripheralComponent Interconnect，外设部件互连标准)总线或EISA(Extended Industry StandardArchitecture，扩展工业标准结构)总线等。总线可以分为地址总线、数据总线、控制总线等。为便于表示，图7中仅用一个双向箭头表示，但并不表示仅有一根总线或一种类型的总线。
存储器，用于存放计算机可执行指令。存储器通过内部总线向处理器提供计算机可执行指令。
处理器，执行存储器所存放的计算机可执行指令，并具体用于实现以下操作：
收集目标区域在研究时间段内所有可用的Landsat系列的地表反射率影像；
将收集的地表反射率影像划分成历史时期和监测期，其中历史时期数据用于构建各月的植被正常模式，监测期数据用于植被异常的即时探测；
分别对历史时期数据和监测期中的待探测影像进行预处理，并将历史时期中各年各月的所有数据通过均值合成方式处理为历史各年逐月数据；
分别对预处理后的历史各年逐月数据和待探测影像逐像元计算多种植被指数；
对各像元的历史各年逐月数据的植被指数逐月份采用两次时间序列线性回归依次去除噪声和植被生长趋势信息，并根据第二次线性回归结果得到各像元的植被指数在各历史时间点的实际值与回归值之差的历史残差，以及得到各像元的植被指数在未来对应月的预测值；
将各像元的植被指数的历史残差作为样本数据进行非参数的核密度估计以构建各像元各月的植被正常模式，并根据设定概率阈值确定出各像元各月的植被正常模式的边界；
根据待探测影像所属月份找到对应月的植被正常模式，结合各像元的植被指数在未来对应月的预测值和植被正常模式的边界，逐像元判别目标区域的植被是否发生异常。
上述如本发明图6所示实施例揭示的植被异常的遥感即时探测装置执行的功能可以应用于处理器中，或者由处理器实现。处理器可能是一种集成电路芯片，具有信号的处理能力。在实现过程中，上述方法的各步骤可以通过处理器中的硬件的集成逻辑电路或者软件形式的指令完成。上述的处理器可以是通用处理器，包括中央处理器(CentralProcessing Unit，CPU)、网络处理器(Network Processor，NP)等；还可以是数字信号处理器(Digital Signal Processor，DSP)、专用集成电路(Application Specific IntegratedCircuit，ASIC)、现场可编程门阵列(Field－Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等，用于实现或者执行本发明各实施例公开的方法、步骤及流程框架图。通用处理器可以是微处理器或者任何常规的处理器。结合本发明所公开的方法步骤可以直接体现为硬件译码处理器执行完成，或者用译码处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器，闪存、只读存储器，可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器，处理器读取存储器中的信息，结合其硬件完成本发明方法的各步骤。
该电子设备还可执行图1中植被异常的遥感即时探测方法执行的步骤，并实现图2所示植被异常的遥感即时探测的流程框架图的功能，本发明在此不再赘述。
本发明还提出了一种计算机可读存储介质，该计算机可读存储介质存储一个或多个程序，该一个或多个程序当被处理器执行时，实现前述的植被异常的遥感即时探测方法的各个实施例，本发明在此也不再赘述。
下面是两个应用实施例：
实施例1
实施例1位于哥伦比亚亚马孙森林，大致中心经纬度1°56′9″N，73°28′34″W，主要植被类型为热带雨林，主要异常类型为森林砍伐。图8示出了本发明实施例1的哥伦比亚亚马孙森林的Landsat影像数据，其中历史时期1985年至2017年，监测期2018年，“LT04”、“LT05”、“LT07”、“LT08”分别代表Landsat 4、5、7、8。选择了该实施例在2018年云量较少的日期进行即时探测实验，分别为20180129、20180222、20180326、20180716、20181231。
图9示出了本发明实施例1的各待即时探测日期的短波红外2波段图及对应的即时探测结果图。其中图9(a-e)为各待即时探测日期的短波红外2波段图，图9(f-j)为各日期对应的即时探测结果图。结果图中background背景值有两种情况：一是该像元为已去除的云像元，不再参与异常探测判别；二是该像元处的历史时间序列数据数量不足以进行两次线性回归，无法进行后续的核密度估计，故也不参与异常探测判别。
从结果来看，即时探测结果与短波红外2波段图所呈现情况基本一致。具体来说，1月29日至3月26日，砍伐迹地出现扩张；到了7月16日，半年内部分被砍伐区域由于发生次生演替出现了浅层植被，在本发明的方法中不再认为其为异常；12月31日，在原始砍伐迹地上出现了小区域的砍伐扩张，曾经的被砍伐区域几乎也出现了浅层植被的生长。
实施例2
实施例2位于内蒙古根河镇大兴安岭，大致中心经纬度51°15′33″N，121°50′47″E，主要指被类型为兴安落叶松、白桦与落叶松混交林等，在2003年5月5日发生了森林火灾。图10示出了本发明实施例2的内蒙古根河镇大兴安岭的Landsat影像数据，其中历史时期1986年至2002年，监测期2003年。选择的即时探测日期为20030126、20030526、20030611、20030814。
图11示出了本发明实施例2的各待即时探测日期的近红外波段图及对应的即时探测结果图，其中图11(a-d)为各待即时探测日期的近红外波段图，图11(e-h)为各日期对应的即时探测结果图。背景值的情况如上述实施例1森林砍伐结果中所述。
从结果来看，即时探测结果与近红外波段所呈现情况也基本一致。具体来说，2003年1月26日，部分区域因为数据不足无法参与异常检测判别，其余区域均判别为正常，与实际相符；5月26日是火灾发生后Landsat获取的第一幅影像，通过与近红外波段图对比，虽然存在漏检与错检，该方法能探测到绝大部分的火灾异常；6月11日至8月14日，部分火烧迹地出现了次生演替且随时间推移增加，在本发明的方法中不再探测为异常。
以上两个应用实施例表明，本发明的方法对不同区域不同异常类型有一定普适性，能够实现在有可用影像时，对植被异常进行即时探测。

Claims (10)

1.一种植被异常的遥感即时探测方法，其特征在于，包括：
收集目标区域在研究时间段内所有可用的Landsat系列的地表反射率影像；
将收集的地表反射率影像划分成历史时期和监测期，其中历史时期数据用于构建各月的植被正常模式，监测期数据用于植被异常的即时探测；
分别对历史时期数据和监测期中的待探测影像进行预处理，并将历史时期中各年各月的所有数据通过均值合成方式处理为历史各年逐月数据；
分别对预处理后的历史各年逐月数据和待探测影像逐像元计算多种植被指数；
对各像元的历史各年逐月数据的植被指数逐月份采用两次时间序列线性回归依次去除噪声和植被生长趋势信息，并根据第二次线性回归结果得到各像元的植被指数在各历史时间点的实际值与回归值之差的历史残差，以及得到各像元的植被指数在未来对应月的预测值；
将各像元的植被指数的历史残差作为样本数据进行非参数的核密度估计以构建各像元各月的植被正常模式，并根据设定概率阈值确定出各像元各月的植被正常模式的边界；
根据待探测影像所属月份找到对应月的植被正常模式，结合各像元的植被指数在未来对应月的预测值和植被正常模式的边界，逐像元判别目标区域的植被是否发生异常。
2.根据权利要求1所述的方法，其特征在于，分别对预处理后的历史各年逐月数据和待探测影像逐像元计算的植被指数包括但不限于：归一化差异植被指数NDVI和归一化燃烧比指数NBR；其中，
NDVI指数放大了近红外波段与红波段的反射率差异，可表征植被覆盖程度，其计算公式为：
NBR指数放大了近红外波段与短波红外波段的反射率差异，可表征植被叶片缺水程度，其计算公式为：
其中的RNIR、Rred、RSWIR2分别为近红外、红以及短波红外2波段的反射率。
3.根据权利要求2所述的方法，其特征在于，所述对各像元的历史各年逐月数据的植被指数逐月份采用两次时间序列线性回归依次去除噪声和植被生长趋势信息包括：
针对每个像元的历史各年逐月数据，逐月份利用NDVI指数进行第一次时间序列线性回归，若某个时间点实际值低于回归值，且两者之间的差异大于设定剔除阈值，则将该时间点作为噪声点进行剔除，同时NBR指数中对应时间点的数据也被剔除；
对NDVI指数和NBR指数保留下来的时间点进行第二次时间序列线性回归，得到各时间点的实际值与回归值之差以去除植被生长趋势信息。
4.根据权利要求1所述的方法，其特征在于，所述将各像元的植被指数的历史残差作为样本数据进行非参数的核密度估计以构建各像元各月的植被正常模式，包括：
定义X＝(X1，X2，......，Xn]T为样本数据，每个样本包含了参与构建的N个植被指数的历史残差：
Xi＝(Index1i，Index2i，......,IndexNi),i＝1,......,n
则N维核密度估计表达式为：
其中：x为N维植被指数空间中的点，h为核函数带宽，K()代表高斯核函数，||x-Xi||代表N维空间中x与样本点的距离。
5.根据权利要求4所述的方法，其特征在于，在核密度估计为N维时，采用经验公式h＝1/(N+3)n(-1/(N+4))来确定核函数带宽。
6.根据权利要求4所述的方法，其特征在于，所述根据设定概率阈值确定出各像元各月的植被正常模式的边界，包括：
利用下式确定出各像元各月的植被正常模式的边界：
其中：为植被正常模式的边界，即截面对应的概率密度，可由/>的/>分位数计算得到，/>为截面以上所包含的区域；
由上式可知，该区域至少包含正常的概率为一旦设定概率阈值即可确定/>从而确定植被正常模式分布的边界。
7.根据权利要求1所述的方法，其特征在于，所述根据待探测影像所属月份找到对应月的植被正常模式，结合各像元的植被指数在未来对应月的预测值和植被正常模式的边界，逐像元判别目标区域的植被是否发生异常，包括：
根据待探测影像所属月份找到对应月的植被正常模式，逐像元将待探测影像的植被指数的实际值与对应月的预测值之差的预测值残差带入对应月的植被正常模式，将得到的概率密度与对应边界的概率密度相对比，若待探测影像的某个像元的概率密度小于对应边界的概率密度，则判别与该像元对应区域的植被异常，否则判别与该像元对应区域的植被正常。
8.一种植被异常的遥感即时探测装置，其特征在于，包括：
影像收集单元，用于收集目标区域在研究时间段内所有可用的Landsat系列的地表反射率影像；
时期划分单元，用于将收集的地表反射率影像划分成历史时期和监测期，其中历史时期数据用于构建各月的植被正常模式，监测期数据用于植被异常的即时探测；
预处理单元，用于分别对历史时期数据和监测期中的待探测影像进行预处理，并将历史时期中各年各月的所有数据通过均值合成方式处理为历史各年逐月数据；
植被指数计算单元，用于分别对预处理后的历史各年逐月数据和待探测影像逐像元计算多种植被指数；
时间序列处理单元，用于对各像元的历史各年逐月数据的植被指数逐月份采用两次时间序列线性回归依次去除噪声和植被生长趋势信息，并根据第二次线性回归结果得到各像元的植被指数在各历史时间点的实际值与回归值之差的历史残差，以及得到各像元的植被指数在未来对应月的预测值；
正常模式构建单元，用于将各像元的植被指数的历史残差作为样本数据进行非参数的核密度估计以构建各像元各月的植被正常模式，并根据设定概率阈值确定出各像元各月的植被正常模式的边界；
异常判别单元，用于根据待探测影像所属月份找到对应月的植被正常模式，结合各像元的植被指数在未来对应月的预测值和植被正常模式的边界，逐像元判别目标区域的植被是否发生异常。
9.一种电子设备，包括：处理器和存储计算机可执行指令的存储器，所述可执行指令在被所述处理器执行时，实现权利要求1-7任一项所述的植被异常的遥感即时探测方法。
10.一种计算机可读存储介质，所述计算机可读存储介质存储一个或多个程序，所述一个或多个程序当被处理器执行时，实现权利要求1-7任一项所述的植被异常的遥感即时探测方法。
CN202310524502.3A 2023-05-11 2023-05-11 一种植被异常的遥感即时探测方法和装置 Pending CN116682003A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310524502.3A CN116682003A (zh) 2023-05-11 2023-05-11 一种植被异常的遥感即时探测方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310524502.3A CN116682003A (zh) 2023-05-11 2023-05-11 一种植被异常的遥感即时探测方法和装置

Publications (1)

Publication Number Publication Date
CN116682003A true CN116682003A (zh) 2023-09-01

Family

ID=87782704

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310524502.3A Pending CN116682003A (zh) 2023-05-11 2023-05-11 一种植被异常的遥感即时探测方法和装置

Country Status (1)

Country Link
CN (1) CN116682003A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117688505A (zh) * 2024-02-04 2024-03-12 河海大学 一种植被大范围区域化负异常的预测方法及系统

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117688505A (zh) * 2024-02-04 2024-03-12 河海大学 一种植被大范围区域化负异常的预测方法及系统
CN117688505B (zh) * 2024-02-04 2024-04-19 河海大学 一种植被大范围区域化负异常的预测方法及系统

Similar Documents

Publication Publication Date Title
Hawbaker et al. The Landsat Burned Area algorithm and products for the conterminous United States
Verbesselt et al. Evaluating satellite and climate data-derived indices as fire risk indicators in savanna ecosystems
Berberoglu et al. Assessing different remote sensing techniques to detect land use/cover changes in the eastern Mediterranean
Li et al. Satellite-based mapping of Canadian boreal forest fires: evaluation and comparison of algorithms
Al-Doski et al. NDVI differencing and post-classification to detect vegetation changes in Halabja City, Iraq
Brown et al. Multitemporal, moderate-spatial-resolution remote sensing of modern agricultural production and land modification in the Brazilian Amazon
Shang et al. Near-real-time monitoring of land disturbance with harmonized Landsats 7–8 and Sentinel-2 data
CN116682003A (zh) 一种植被异常的遥感即时探测方法和装置
CN116091497A (zh) 遥感变化检测方法、装置、电子设备和存储介质
CN116844053B (zh) 一种小麦种植区域识别方法、系统、电子设备及存储介质
CN112364681A (zh) 基于二维表的植被覆盖度估算方法及装置
Maynard et al. Effect of spatial image support in detecting long-term vegetation change from satellite time-series
CN113887324A (zh) 基于卫星遥感数据的火点检测方法
CN113052153B (zh) 遥感反射率影像的检验方法、装置、电子设备及存储介质
MacLean et al. Applicability of multi-date land cover mapping using Landsat-5 TM imagery in the Northeastern US
Collins et al. Multi-temporal analysis of landsat data to determine forest age classes for the mississippi statewide forest inventory~ preliminary results
CN115795237A (zh) 一种基于四元Copula的综合遥感生态指数并行计算方法
CN114494281A (zh) 基于增强型燃烧指数的中小型火烧迹地自动提取方法
Ngcofe et al. The South African land cover change detection derived from 2013_2014 and 2017_2018 land cover products
Sempio et al. Assessment of different image transformation methods on diwata-1 smi images using structural similarity measure
Ahmad et al. Haze effects on satellite remote sensing imagery and their corrections
CN113640226B (zh) 绿潮空间覆盖率的反演方法、装置及电子设备
Cal et al. Automatic Classification of Agricultural Summer Crops in Uruguay
Flores-Romero Improving ocean chlorophyll estimation in satellite hyperspectral images using ensemble machine learning
CN114299385B (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