CN112101159B - 多时相林业遥感影像变化监测方法 - Google Patents

多时相林业遥感影像变化监测方法 Download PDF

Info

Publication number
CN112101159B
CN112101159B CN202010921885.4A CN202010921885A CN112101159B CN 112101159 B CN112101159 B CN 112101159B CN 202010921885 A CN202010921885 A CN 202010921885A CN 112101159 B CN112101159 B CN 112101159B
Authority
CN
China
Prior art keywords
image
land
forest
value
current
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
CN202010921885.4A
Other languages
English (en)
Other versions
CN112101159A (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.)
Central South Investigation Planning And Design Institute Of State Forestry And Grassland Administration
Original Assignee
Central South Investigation Planning And Design Institute Of State Forestry And Grassland Administration
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 Central South Investigation Planning And Design Institute Of State Forestry And Grassland Administration filed Critical Central South Investigation Planning And Design Institute Of State Forestry And Grassland Administration
Priority to CN202010921885.4A priority Critical patent/CN112101159B/zh
Publication of CN112101159A publication Critical patent/CN112101159A/zh
Application granted granted Critical
Publication of CN112101159B publication Critical patent/CN112101159B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/045—Combinations of networks
    • 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/20—Image preprocessing
    • G06V10/26—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion
    • G06V10/267—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion by performing operations on regions, e.g. growing, shrinking or watersheds

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Health & Medical Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Multimedia (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Evolutionary Biology (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种多时相林业遥感影像变化监测方法，包括：监测区域影像转化灰度图；计算目标地类的灰度值的阈值；根据目标地类的灰度值的阈值将遥感影像归一化，得到归一化后的前、后两期影像；分别将不同目标地类对应的归一化后的前、后两期影像进行差值计算得到影像变化图；对不同目标地类的影像变化图分别通过众数滤波和交点统计后进行结果整合，得到初步变化检测的矢量图，剔除破碎斑块后通过深度神经网络模型的分类与筛选，得到最终变化检测结果的矢量图；分析最终变化检测结果的矢量图的变化原因。本发明根据森林资源管理“一张图”年度更新成果对监测区域的地形进行精确划分并进行遥感变化检测，得到准确的变化检测结果。

Description

多时相林业遥感影像变化监测方法
技术领域
本发明涉及遥感变化监测领域，尤其涉及一种多时相林业遥感影像变化监测方法。
背景技术
遥感变化监测是利用多时相的遥感影像，采用多种图像识别方法提取变化信息，并定量分析，确定地表变化的特征与过程。它涉及变化的类型、分布状况与变化量，即需要确定变化前、后的地面类型、界线、及变化趋势，进而分析这些动态变化的特点与形成原因。遥感变化监测按照处理对象的级别分为像素级、特征级和目标级三种方法。其中特征级及目标级变化检测的研究尚未成熟，难以成规模的在林业实践工作中使用。
像素级变化检测包括差值图像的方法和其他像素变化检测方法，相较于差值图像的方法，其他像素变化检测方法的普适性较弱，一般选用差值图像的方法。差值图像的方法中：
阈值化方法选定阈值对差值图像进行阈值化从而区分变化类与非变化类像素，其优点是方法简单明了，缺点是未充分利用图像上下文信息，容易受到传感器噪声、配准误差等因素的影响；
模式分类方法利用分类器对不同样本数据集进行分类，从而得到变化检测结果，其优点是克服简单的阈值化方法的不准确性，缺点是无法利用空间上下文信息，并且人工神经网络和支持向量机的方法须人工干预，提供监督信息；
多尺度分析方法利用多尺度几何分析方法进行变化检测，从而克服传感器噪声、配准误差等因素的影响，其优点是克服传感器噪声、配准误差等因素，缺点是细节丢失、计算时间增加以及多尺度结果处理技术难点尚未解决；
马尔科夫随机场方法利用马尔科夫随机场模型模拟空间上下文信息，从而得到变化检测结果，其优点是克服了差值图像进行参数化概率模型假设的不合理性，利用上下文信息，缺点是针对高分影像计算量太大，计算时长过长。
森林资源二类调查的地类划分标准将林地划分为林地和非林地两大类；其中，林地包括有林地(0110)、疏林地(0120)、灌木林地(0130)、未成林地(0140)、苗圃地(0150)、无林木林地(0160)、宜林地(0170)和林业辅助生产用地(0180)；非林地包括耕地(0210)、牧草地(0220)、水域(0230)、未利用地(0240)和建设用地(0250)。森林资源管理“一张图”年度更新在森林资源二类调查的基础上，每年进行更新，从而及时、准确反映上一年度的森林资源植被覆盖的实际情况。
如何依据森林资源二类调查的地类划分标准，并结合森林资源管理“一张图”年度更新成果所反映的植被覆盖的实际情况对遥感变化监测结果进行处理以得到更加准确的结果值得研究和探讨。
发明内容
本发明要解决的技术问题就在于：针对现有技术所存在的技术问题，本发明提供一种多时相林业遥感影像变化监测方法，根据森林资源管理“一张图”年度更新成果对监测区域的地形进行精确划分并进行遥感变化监测，得到准确的监测结果。
为解决上述技术问题，本发明提出的技术方案为：
一种多时相林业遥感影像变化监测方法，包括以下步骤：
S1)监测区域影像转化灰度图：分别将监测区域的前、后两期影像的红光波段转化为灰度图；
S2)计算目标地类的灰度值的阈值：根据森林资源二类调查的地类划分标准以及森林资源管理“一张图”年度更新成果，分析不同地类对应的红光波段灰度值的大小关系，然后根据目标地类在后期影像中的变化，计算目标地类对应的灰度值的阈值，目标地类包括有植被覆盖的林地以及无植被覆盖的林地；
S3)遥感影像归一化：分别根据有植被覆盖的林地以及无植被覆盖的林地对应的灰度值的阈值处理前、后两期影像，得到有植被覆盖的林地以及无植被覆盖的林地对应的处理后的前、后两期影像，然后对处理后的前、后两期影像进行归一化，得到有植被覆盖的林地对应的归一化后的前期影像A1和后期影像A2，以及无植被覆盖的林地对应的归一化后的前期影像A1’和后期影像A2’；
S4)归一化结果差值计算：针对有植被覆盖的林地对应的归一化后的前期影像A1以及后期影像A2，将归一化后的前、后两期影像、中对应的像元一一进行差值计算得到对应的影像变化图B，针对无植被覆盖的林地对应的归一化后的前期影像A1’以及后期影像A2’，将归一化后的前、后两期影像、中的对应的像元一一进行差值计算得到对应的影像变化图B’；
S5)众数滤波：分别利用众数滤波修复影像变化图B和B’的图形空洞现象，得到有植被覆盖的林地以及无植被覆盖的林地对应的修复后的影像变化图B1和B1’；
S6)焦点统计：分别利用焦点统计合并修复后的影像变化图B1和B1’的图中斑块的临近碎斑现象，得到有植被覆盖的林地以及无植被覆盖的林地对应的二次修复后的影像变化图B2和B2’；
S7)结果整合：将二次修复后的影像变化图B2和B2’依次进行二值化和矢量化后合并得到初步变化检测的矢量图D；
S8)剔除破碎斑块：对初步变化检测的矢量图D求算实际面积，剔除图中实际面积小于预设值的破碎斑块；
S9)深度神经网络模型的分类与筛选：将剔除破碎斑块后的初步变化检测的矢量图D输入深度神经网络模型，以初步变化检测的矢量图D为进一步判定对象范围，通过神经网络模型进行分类和筛选后剔除图中影像辐射误差造成的未发生变化的图斑，得到最终变化检测的矢量图D1；
S10)分析变化原因：提取最终变化检测的矢量图D1中图斑的影像平均灰度值，根据灰度值分析变化原因。
进一步的，步骤S2)中计算目标地类对应的灰度值的阈值的步骤具体包括：
S21)选取有植被覆盖的林地或无植被覆盖的林地作为当前地类，选取前期影像或后期影像作为当前影像，在当前影像中找到当前地类对应的参考地类所在区域，将参考地类所在区域内的面数据转化为点数据，有植被覆盖的林地对应的参考地类为耕地，无植被覆盖的林地对应的参考地类为建设用地或未利用地；
S22)获取所有点数据的灰度值，先通过箱线图分析来进行异常值检测，剔除异常值后再进行K均值聚类得到至少一个聚类中心，选择其中最小的聚类中心作为当前影像中当前地类对应的灰度值的阈值。
进一步的，步骤S3)具体包括以下步骤：
S31)选取有植被覆盖的林地或无植被覆盖的林地作为当前地类，根据后期影像中当前地类对应的灰度值的阈值，确定当前地类对应的处理后的后期影像中所有像元的灰度值，函数表达式如下：
Figure BDA0002667001730000031
上式中，xi为当前地类对应的处理后的后期影像中的像元i的灰度值，E1为后期影像中当前地类对应的灰度值的阈值；
S32)根据前期影像中当前地类对应的灰度值的阈值，确定当前地类对应的处理后的前期影像中所有像元的灰度值，函数表达式如下：
Figure BDA0002667001730000032
上式中，yi为当前地类对应的处理后的前期影像中的像元i的灰度值，E2为前期影像中当前地类对应的灰度值的阈值；
S33)将当前地类对应的处理后的前、后两期影像中像元的灰度值归一化，得到当前地类对应的归一化后的前、后两期影像，函数表达式如下：
Figure BDA0002667001730000041
上式中，X'i是当前地类对应的归一化后的后期影像中像元i的灰度值，Y'i是当前地类对应的归一化后的前期影像中像元i的灰度值，xi是当前地类对应的处理后的后期影像中像元i的灰度值，yi是当前地类对应的处理后的前期影像中像元i的灰度值，ximax为当前地类对应的处理后的后期影像中各像元灰度的最大值，ximin为当前地类对应的处理后的后期影像中各像元灰度的最小值，yimax为当前地类对应的处理后的前期影像中各像元灰度的最大值，yimin为当前地类对应的处理后的前期影像中各像元灰度的最小值。
进一步的，步骤S4)中将归一化后的前、后两期影像中的对应的像元一一进行差值计算的函数表达式如下：
Figure BDA0002667001730000042
上式中，X'i是归一化后的后期影像中的像元i的灰度值，Y'i是归一化后的前期影像中的像元i的灰度值，e是预设的差异性阈值。
进一步的，步骤S5)中众数滤波的具体步骤包括：针对影像变化图B或B’，选取影像变化图中的一个像元作为当前像元，然后以当前像元为中心，设置一个像元为距离的矩形窗口，若矩形窗口中当前像元周边的其他像元中超过八分之五的相邻像元具有相同值，将这些相邻像元的值替换当前像元的值。
进一步的，步骤S6)中焦点统计的具体步骤包括：针对修复后的影像变化图B1或B1’，选取修复后的影像变化图中的一个像元作为当前像元，然后以当前像元为中心，设置两个像元为距离的矩形窗口，统计矩形窗口中所有像元的值的平均值，将所有像元的值的平均值作为当前像元的值。
进一步的，步骤S7)具体包括以下步骤：
S71)分别将影像变化图B和B’进行二值化，得到有植被覆盖的林地以及无植被覆盖的林地对应的初步变化监测结果的二值图C0和C0’；
S72)分别对二值图C0和C0’矢量化，得到矢量图D0和D0’；
S73)将矢量化后的矢量图D0和D0’合并得到初步变化检测的矢量图D。
进一步的，步骤S9)具体包括以下步骤：
S91)将剔除破碎斑块后的初步变化检测的矢量图D输入深度神经网络模型；
S92)以初步变化检测的矢量图D为进一步判定对象范围，通过神经网络模型进行分类和筛选后剔除图中影像辐射误差造成的未发生变化的图斑，得到最终变化检测的矢量图D1。
进一步的，步骤S10)具体包括以下步骤：
S101)将最终变化检测的矢量图D1中图斑的所有面数据转化为点数据；
S102)根据点数据提取对应像元的灰度值；
S103)对所有被提取的像元的灰度值先通过箱线图分析来进行异常值检测，剔除异常值后再进行K均值聚类，将最小聚类中心的灰度值作为聚类分析后的灰度值，根据聚类分析后的灰度值与有植被覆盖的林地以及无植被覆盖的林地的灰度值的阈值的大小关系，分析变化原因。
本发明还提出一种计算机可读存储介质，所述计算机可读存储介质存储有被编程或配置以实现上述的多时相林业遥感影像变化监测方法的计算机程序。
与现有技术相比，本发明的优点在于：
(1)本发明依据森林资源二类调查的地类划分标准，并结合森林资源管理“一张图”年度更新成果，将监测区域进行地类划分为有植被覆盖的林地和无植被覆盖的林地后，分别进行变化检测，相比于现有遥感变化监测方法，单独选取林地类型能够有效的监测无植被覆盖的林地类型的变化以及有植被覆盖的林地不同地类相互之间的转化；
(2)本发明采用设置阈值的归一化法，得到初步的变化检测结果，并引进机器学习的方法，具体利用深度神经网络模型进行分类和筛选，一是有效降低了不同时期遥感影像差异所造成的影响，提高普适性；二是避免了因样本选取过程中人为因素造成的模型的误差，提高精度；三是降低运算时长，提高效率；四是降低发明对硬件条件的强制需求，利于推广；
(3)本发明引入图像形态学中膨胀的方法，利用众数滤波和焦点统计的方式，通过矩阵的思路，有效避免了影像变化图中斑块的临近碎斑及异常的图形空洞的现象。
附图说明
图1为本发明实施例的步骤示意图。
图2为本发明实施例的具体流程图。
图3为众数滤波演示图。
图4为焦点统计演示图。
具体实施方式
以下结合说明书附图和具体优选的实施例对本发明作进一步描述，但并不因此而限制本发明的保护范围。
如图1所示，本发明的多时相林业遥感影像变化监测方法包括以下步骤：
S1)将监测区域影像转化灰度图：为方便基层林业部门使用，国家林业和草原局下发的高分影像为经过预处理的假彩色影像，分别由红外、红、绿三个波段融合而成，其中，红光波段最适合进行林业遥感影像变化检测，因此本实施例中，分别将监测区域的前、后两期影像的红光波段转化为灰度图；
S2)计算目标地类的灰度值的阈值：森林资源二类调查的地类划分标准将地类分为林地和非林地两大类，对于红光波段，灰度值越高的地块植被覆盖率越低，现有技术不区分有植被覆盖的林地和无植被覆盖的林地，根据森林资源管理“一张图”年底更新成果我们发现：
对于非林地而言，不同地类对应的红光波段灰度值从大到小进行排序为建设用地、未利用地、草地、耕地、水域；
将林地和非林地整体考虑，地类对应的红光波段灰度值从大到小进行排序为建设用地、未利用地、无植被覆盖的林地、草地、耕地、有植被覆盖的林地；
此外，根据地类变化的实际情况及红光波段灰度值的差异，有植被覆盖的林地在后期影像中可变为建设用地、未利用地、无植被覆盖的林地、草地、耕地，无植被覆盖的林地在后期影像中可变为建设用地和未利用地；
因此，本实施例中，根据森林资源二类调查的地类划分标准以及森林资源管理“一张图”年底更新数据，分析不同地类对应的红光波段的灰度值的大小关系，然后根据目标地类在后期影像中的地类变化计算目标地类对应的灰度值的阈值，林业变化检测主要针对林地变为非林地或无植被覆盖的林地的情况，本实施例的目标地类包括有植被覆盖的林地以及无植被覆盖的林地，对于有植被覆盖的林地而言，能够转化为非林地或无植被覆盖的林地，根据不同地类对应的红光波段的灰度值的大小关系，有植被覆盖的林地对应的灰度值的阈值为地类为耕地的区域聚类分析后的灰度值，后期影像中有植被覆盖的林地转化后的区域的像元的灰度值应大于影像中地类为耕地的区域聚类分析后的灰度值(即后期影像中有植被覆盖的林地转化为非林地或无植被覆盖的林地)，前期影像中有植被覆盖的林地所在区域的像元的灰度值应小于影像中地类为耕地的区域聚类分析后的灰度值(即前期影像中为有植被覆盖的林地)；对于无植被覆盖的林地而言，无植被覆盖的林地能够转化为建设用地或未利用地，根据不同地类对应的红光波段灰度值的大小关系，无植被覆盖的林地对应的灰度值的阈值为地类为建设用地或未利用地的区域聚类分析后的灰度值，后期影像中无植被覆盖的林地转化后的区域的像元的灰度值应大于影像中地类为建设用地或未利用地的区域聚类分析后的灰度值(即后期影像中无植被覆盖的林地转化为非林地)，前期影像中无植被覆盖的林地所在区域的像元的灰度值应小于影像中地类为建设用地或未利用地的区域聚类分析后的灰度值(即前期影像中为无植被覆盖的林地)；
根据上述内容，本实施例中得到目标地类对应的灰度值的阈值的步骤包括：
S21)选取有植被覆盖的林地或无植被覆盖的林地作为当前地类，选取前期或后期影像作为当前影像，在当前影像中找到当前地类对应的参考地类所在区域，将参考地类所在区域内的面数据转化为点数据，有植被覆盖的林地对应的参考地类为耕地，无植被覆盖的林地对应的参考地类为建设用地或未利用地；
S22)获取所有点数据的灰度值，先通过箱线图分析来进行异常值检测，剔除异常值后再进行K均值聚类得到至少一个聚类中心，选择其中最小的聚类中心作为当前影像中当前地类对应的灰度值的阈值；
根据步骤S21)至S22)即可得到前期影像中有植被覆盖的林地对应的灰度值的阈值、前期影像中无植被覆盖的林地对应的灰度值的阈值、后期影像中有植被覆盖的林地对应的灰度值的阈值、后期影像中无植被覆盖的林地对应的灰度值的阈值；
S3)遥感影像归一化：根据步骤S2)中的分析内容，本实施例中将目标地类对应的灰度值的阈值作为目标地类对应的林地与非林地的分界值，同时本实施例中灰度值的阈值为目标地类对应的参考地类所在区域剔除异常值后最小聚类中心的灰度值，可以忽略林地与非林地之间数据的交叉，如图2所示，得到有植被覆盖的林地以及无植被覆盖的林地对应的灰度值的阈值之后，分别根据有植被覆盖的林地以及无植被覆盖的林地对应的灰度值的阈值处理前、后两期影像，选取有植被覆盖以及无植被覆盖的林地对应的处理后的前、后两期影像的范围，使得处理后的前期影像不含非林地，处理后的后期影像不含林地，从而减小后期计算量，然后对处理后的前期影像和处理后的后期影像进行归一化，得到有植被覆盖的林地对应的归一化后的前期影像A1和后期影像A2，以及无植被覆盖的林地对应的归一化后的前期影像A1’和后期影像A2’；
S4)归一化结果差值计算：如图2所示，针对有植被覆盖的林地，将对应的归一化后的前期影像A1以及后期影像A2按照归一化后的前、后两期影像、中对应的像元一一进行差值计算，得到对应的影像变化图B，针对无植被覆盖的林地，将对应的归一化后的前期影像A1’以及后期影像A2’按照归一化后的前、后两期影像中对应的像元一一进行差值计算得到对应的影像变化图B’，影像变化图B和B’即为有植被覆盖的林地和无植被覆盖的林地的变化检测结果；
S5)众数滤波：如图2所示，针对有植被覆盖的林地对应的影像变化图B以及无植被覆盖的林地对应的影像变化图B’，分别利用众数滤波修复图形空洞现象，得到有植被覆盖的林地对应的修复后的影像变化图B1和无植被覆盖的林地对应的修复后的影像变化图B1’；
S6)焦点统计：如图2所示，针对有植被覆盖的林地对应的修复后的影像变化图B1和无植被覆盖的林地对应的修复后的影像变化图B1’，分别利用焦点统计合并图中斑块的临近碎斑现象得到有植被覆盖的林地对应的二次修复后的影像变化图B2和无植被覆盖的林地对应的二次修复后的影像变化图B2’；
S7)结果整合：如图2所示，将有植被覆盖的林地对应的二次修复后的影像变化图B2和无植被覆盖的林地对应的二次修复后的影像变化图B2’进行二值化，再对二值化后的结果矢量化，矢量化后的结果中的图斑即为图形修复后的有植被覆盖的林地以及无植被覆盖的林地的变化区域，然后将矢量化后的结果合并得到初步变化检测的矢量图D，初步变化检测的矢量图D中的图斑即为本实施例初步得到的林地变为非林地或无植被覆盖的林地的区域；
S8)剔除破碎斑块：如图2所示，对初步变化检测的矢量图D求算实际面积，剔除图中实际面积小于预设值的破碎斑块；
S9)深度神经网络模型的分类与筛选：如图2所示，将剔除破碎斑块后的初步变化检测的矢量图D输入深度神经网络模型，以初步变化检测的矢量图D为进一步判定对象范围，通过神经网络模型进行分类和筛选后剔除图中影像辐射误差造成的未发生变化的图斑，得到最终变化检测的矢量图D1，最终变化检测的矢量图D1中的图斑即为本实施例得到的林地变为非林地或无植被覆盖的林地的精确区域；
S10)分析变化原因：如图2所示，提取最终变化检测的矢量图D1中图斑的灰度值，根据灰度值分析变化原因。
本实施例的步骤S3)具体包括以下步骤：
S31)选取有植被覆盖的林地或无植被覆盖的林地作为当前地类，根据后期影像中当前地类对应的灰度值的阈值，确定当前地类对应的处理后的后期影像中所有像元的灰度值，函数表达式如下：
Figure BDA0002667001730000081
上式中，xi为当前地类对应的处理后的后期影像中的像元i的灰度值，E1为后期影像中当前地类对应的灰度值的阈值；
S32)根据前期影像中当前地类对应的灰度值的阈值，确定当前地类对应的处理后的前期影像中所有像元的灰度值，函数表达式如下：
Figure BDA0002667001730000082
上式中，yi为当前地类对应的处理后的前期影像中的像元i的灰度值，E2为前期影像中当前地类对应的灰度值的阈值；
S33)将当前地类对应的处理后的前、后两期影像中像元的灰度值归一化，得到当前地类对应的归一化后的前、后两期影像，函数表达式如下：
Figure BDA0002667001730000091
上式中，X'i是当前地类对应的归一化后的后期影像中像元i的灰度值，Y'i是当前地类对应的归一化后的前期影像中像元i的灰度值，xi是当前地类对应的处理后的后期影像中像元i的灰度值，yi是当前地类对应的处理后的前期影像中像元i的灰度值，ximax为当前地类对应的处理后的后期影像中各像元灰度的最大值，ximin为当前地类对应的处理后的后期影像中各像元灰度的最小值，yimax为当前地类对应的处理后的前期影像中各像元灰度的最大值，yimin为当前地类对应的处理后的前期影像中各像元灰度的最小值。
本实施例中，根据步骤S31)至步骤S33)即可得到有植被覆盖的林地对应的归一化后的前、后两期影像，以及无植被覆盖的林地对应的归一化后的前、后两期影像：
步骤S31)中，若当前地类为有植被覆盖的林地，作为林地与非林地的分界值的后期影像中当前地类对应的灰度值的阈值为有植被覆盖的林地对应的灰度值的阈值。因此E1为后期影像中耕地的区域的所有点数据的灰度值剔除异常值后最小聚类中心的灰度值，根据式(1)得到有植被覆盖的林地对应的处理后的后期影像A01，若当前地类为无植被覆盖的林地，作为林地与非林地的分界值的后期影像中当前地类对应的灰度值的阈值为无植被覆盖的林地对应的灰度值的阈值，因此E1为后期影像中建设用地或未利用地的区域的所有点数据的灰度值剔除异常值后最小聚类中心的灰度值，根据式(1)得到无植被覆盖的林地对应的处理后的后期影像A01’。
步骤S32)中，若当前地类为有植被覆盖的林地，作为林地与非林地的分界值的前期影像中当前地类对应的灰度值的阈值为有植被覆盖的林地对应的灰度值的阈值，因此E2为前期影像中耕地的区域的所有点数据的灰度值剔除异常值后最小聚类中心的灰度值，根据式(2)得到有植被覆盖的林地对应的处理后的前期影像A02，若当前地类为无植被覆盖的林地，作为林地与非林地的分界值的前期影像中当前地类对应的灰度值的阈值为无植被覆盖的林地对应的灰度值的阈值，因此E2为前期影像中建设用地或未利用地的区域的所有点数据的灰度值剔除异常值后最小聚类中心的灰度值，根据式(2)得到无植被覆盖的林地对应的处理后的前期影像A02’。
步骤S33)中，若当前地类为有植被覆盖的林地，根据式(3)分别将对应的处理后的前期影像A02与处理后的后期影像A01中像元的灰度值归一化，即可得到对应的归一化后的前期影像A1和归一化后的后期影像A2，若当前地类为无植被覆盖的林地，根据式(3)分别将对应的处理后的前期影像A02’与处理后的后期影像A01’中像元的灰度值归一化即可得到对应的归一化后的前期影像A1’和归一化后的后期影像A2’。
本实施例的步骤S4)中将归一化后的前、后两期影像中的对应的像元一一进行差值计算的函数表达式如下：
Figure BDA0002667001730000101
上式中，X'i-Y'i是影像变化图中像元i的值，X'i是归一化后的后期影像中的像元i的灰度值，Y'i是归一化后的前期影像中的像元i的灰度值，e是预设的差异性阈值，差异性阈值主要根据测试结果来定，也可根据各监测区的实际情况，将地形地貌进行细分，本实施例中通过将地形地貌分为平原区和山区两大类，分别设置差异性阈值。
根据步骤S4)的式(4)，针对有植被覆盖的林地，将对应的归一化后的前期影像A1以及后期影像A2中对应的像元一一进行差值计算就可以得到对应的影像变化图B，针对无植被覆盖的林地，将对应的归一化后的前期影像A1’以及后期影像A2’中对应的像元一一进行差值计算就可以得到对应的影像变化图B’，影像变化图B和B’均为二值图。
本实施例引入图像形态学中膨胀的方法，步骤S5)通过利用众数滤波修复异常的图形空洞的现象，步骤S6)通过利用焦点统计合并斑块的临近碎斑现象。
众数滤波的原理如图3所示，对图中值为“-3”的最近相邻的八个相邻像元(3×3窗口)进行分析，其数据的值涉及“7”、“5”和“4”，其个数分别为5、2、1。数值为“7”的数量最多，像元具有相同值并相邻，超过八分之五的已连接像元具有相同值，则将“-3”替换为“7”。因此本实施例的步骤S5)中众数滤波的具体步骤包括：针对影像变化图B或B’，选取影像变化图中的一个像元作为当前像元，然后以当前像元为中心，设置一个像元为距离的矩形窗口，若矩形窗口中当前像元周边的其他像元中超过八分之五的相邻像元具有相同值，将这些相邻像元的值替换当前像元的值。
焦点统计的原理如图4所示，对图中以当前像元为中心，5×5的矩形为窗口进行分析，统计其及周边24个像元的平均值作为该像元的值，例如，图A中心位置的5，其5×5的窗口的的值总和为125，计算领域内像元的平均值为5，则其统计结果为5。利用该种方式达到合并斑块的临近碎斑现象的需求。因此本实施例的步骤S6)中焦点统计的具体步骤包括：针对修复后的影像变化图B1或B1’，选取修复后的影像变化图中的一个像元作为当前像元，然后以当前像元为中心，设置两个像元为距离的矩形窗口，统计矩形窗口中所有像元的值的平均值，将所有像元的值的平均值作为当前像元的值。
本实施例的步骤S7)具体包括以下步骤：
S71)因步骤S6)的焦点统计时计算得到的二次修复后的影像变化图B2和B2’中像元的值可能存在小数，并不符合二值图的要求，为便于后续的计算和分析因此需要再次进行二值化，本实施例中分别将有植被覆盖的林地以及无植被覆盖的林地对应的二次修复后的影像变化图B和B’进行二值化，二值化的差异性阈值主要根据测试结果来定，也可根据各监测区的实际情况，将地形地貌进行细分，本实施例中通过将地形地貌分为平原区和山区两大类，分别设置差异性阈值，对于二次修复后的影像变化图B和B’，若二次修复后的影像变化图中像元的值小于差异性阈值，则将像元的值置0，若若像元的值大于差异性阈值，则将像元的值置1，从而得到有植被覆盖的林地以及无植被覆盖的林地对应的二值图C0和C0’。
S72)分别对有植被覆盖的林地以及无植被覆盖的林地对应的二值图C0和C0’矢量化，得到矢量图D0和D0’；矢量图D0中的图斑即为图形修复后的有植被覆盖的林地的变化区域，矢量图D0’中的图斑即为图形修复后的无植被覆盖的林地的变化区域；
S73)将矢量化后的有植被覆盖的林地对应的矢量图D0以及无植被覆盖的林地对应的矢量图D0’合并得到初步变化检测的矢量图D。
本实施例的步骤S9)中，分类模型为深度神经网络模型(DNN)，训练样本为上一年度森林督查结果。DNN是最广泛使用的机器学习算法之一，现阶段林业遥感影像均为高分辨率影像，利用深度神经网络的方法对影像进行变化检测耗时过长，对设备的配置要求高。本实施例基于阈值的进行初步变化检测后，利用DNN进行分类与筛选，克服简单的阈值化方法的不准确的同时，利用已有训练样本减少人工干预，步骤S9)具体包括以下步骤：
S91)将剔除破碎斑块后的初步变化检测的矢量图D输入深度神经网络模型；
S92)以初步变化检测的矢量图D为进一步判定对象范围，通过神经网络模型进行分类和筛选后剔除图中影像辐射误差造成的未发生变化的图斑，得到最终变化检测的矢量图D1。
本实施例的步骤S10)具体包括以下步骤：
S101)将最终变化检测的矢量图D1中图斑的所有面数据转化为点数据；
S102)根据点数据提取对应像元的灰度值；
S103)对所有被提取的像元的灰度值先通过箱线图分析来进行异常值检测，剔除异常值后再进行K均值聚类，将最小聚类中心的灰度值作为聚类分析后的灰度值，根据聚类分析后的灰度值与有植被覆盖的林地以及无植被覆盖的林地的灰度值的阈值的大小关系，分析变化原因。
根据森林资源管理“一张图”年度更新的要求，变化原因可分为按建设项目占用林地(10)、林木采伐(20)、开垦(30)、灾害等引起地类或林相变化(40)，以及可识别的因造林更新引起的地类或林相变化(50)和其他(60)。步骤S103)中，若聚类分析后的灰度值更接近建设用地聚类分析后的灰度值，则变化原因为建设项目占用林地或灾害等引起地类或林相变化，若聚类分析后的灰度值更接近耕地聚类分析后的灰度值，则变化原因为林木采伐、开垦、可识别的因造林更新引起的地类或林相变化或其他，然后结合森林资源管理“一张图”年度更新的地类信息，给出具体变化原因填写建议，例如，除胸径大于5cm有蓄积的乔木林地外变化原因不可填为林木采伐。
本发明还提出一种计算机可读存储介质，所述计算机可读存储介质存储有被编程或配置以实现上述的多时相林业遥感影像变化监测方法的计算机程序。
上述只是本发明的较佳实施例，并非对本发明作任何形式上的限制。虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明。因此，凡是未脱离本发明技术方案的内容，依据本发明技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均应落在本发明技术方案保护的范围内。

Claims (10)

1.一种多时相林业遥感影像变化监测方法，其特征在于，包括以下步骤：
S1)监测区域影像转化灰度图：分别将监测区域的前、后两期影像的红光波段转化为灰度图；
S2)计算目标地类的灰度值的阈值：根据森林资源二类调查的地类划分标准以及森林资源管理“一张图”年度更新成果，分析不同地类对应的红光波段灰度值的大小关系，然后根据目标地类在后期影像中的变化，计算目标地类对应的灰度值的阈值，目标地类包括有植被覆盖的林地以及无植被覆盖的林地；
S3)遥感影像归一化：分别根据有植被覆盖的林地以及无植被覆盖的林地对应的灰度值的阈值处理前、后两期影像，得到有植被覆盖的林地以及无植被覆盖的林地对应的处理后的前、后两期影像，然后对处理后的前、后两期影像进行归一化，得到有植被覆盖的林地对应的归一化后的前期影像A1和后期影像A2，以及无植被覆盖的林地对应的归一化后的前期影像A1’和后期影像A2’；
S4)归一化结果差值计算：针对有植被覆盖的林地对应的归一化后的前期影像A1以及后期影像A2，将归一化后的前、后两期影像、中对应的像元一一进行差值计算得到对应的影像变化图B，针对无植被覆盖的林地对应的归一化后的前期影像A1’以及后期影像A2’，将归一化后的前、后两期影像中的对应的像元一一进行差值计算得到对应的影像变化图B’；
S5)众数滤波：分别利用众数滤波修复影像变化图B和B’的图形空洞现象，得到有植被覆盖的林地以及无植被覆盖的林地对应的修复后的影像变化图B1和B1’；
S6)焦点统计：分别利用焦点统计合并修复后的影像变化图B1和B1’的图中斑块的临近碎斑现象，得到有植被覆盖的林地以及无植被覆盖的林地对应的二次修复后的影像变化图B2和B2’；
S7)结果整合：将二次修复后的影像变化图B2和B2’依次进行二值化和矢量化后合并得到初步变化检测的矢量图D；
S8)剔除破碎斑块：对初步变化检测的矢量图D求算实际面积，剔除图中实际面积小于预设值的破碎斑块；
S9)深度神经网络模型的分类与筛选：将剔除破碎斑块后的初步变化检测的矢量图D输入深度神经网络模型，以初步变化检测的矢量图D为进一步判定对象范围，通过神经网络模型进行分类和筛选后剔除图中影像辐射误差造成的未发生变化的图斑，得到最终变化检测的矢量图D1；
S10)分析变化原因：提取最终变化检测的矢量图D1中图斑的影像平均灰度值，根据灰度值分析变化原因。
2.根据权利要求1所述的多时相林业遥感影像变化监测方法，其特征在于，步骤S2)中计算目标地类对应的灰度值的阈值的步骤具体包括：
S21)选取有植被覆盖的林地或无植被覆盖的林地作为当前地类，选取前期影像或后期影像作为当前影像，在当前影像中找到当前地类对应的参考地类所在区域，将参考地类所在区域内的面数据转化为点数据，有植被覆盖的林地对应的参考地类为耕地，无植被覆盖的林地对应的参考地类为建设用地或未利用地；
S22)获取所有点数据的灰度值，先通过箱线图分析来进行异常值检测，剔除异常值后再进行K均值聚类得到至少一个聚类中心，选择其中最小的聚类中心作为当前影像中当前地类对应的灰度值的阈值。
3.根据权利要求1所述的多时相林业遥感影像变化监测方法，其特征在于，步骤S3)包括以下步骤：
S31)选取有植被覆盖的林地或无植被覆盖的林地作为当前地类，根据后期影像中当前地类对应的灰度值的阈值，确定当前地类对应的处理后的后期影像中所有像元的灰度值，函数表达式如下：
Figure FDA0002667001720000021
上式中，xi为当前地类对应的处理后的后期影像中的像元i的灰度值，E1为后期影像中当前地类对应的灰度值的阈值；
S32)根据前期影像中当前地类对应的灰度值的阈值，确定当前地类对应的处理后的前期影像中所有像元的灰度值，函数表达式如下：
Figure FDA0002667001720000022
上式中，yi为当前地类对应的处理后的前期影像中的像元i的灰度值，E2为前期影像中当前地类对应的灰度值的阈值；
S33)将当前地类对应的处理后的前、后两期影像中像元的灰度值归一化，得到当前地类对应的归一化后的前、后两期影像，函数表达式如下：
Figure FDA0002667001720000023
上式中，X′i是当前地类对应的归一化后的后期影像中像元i的灰度值，Yi'是当前地类对应的归一化后的前期影像中像元i的灰度值，xi是当前地类对应的处理后的后期影像中像元i的灰度值，yi是当前地类对应的处理后的前期影像中像元i的灰度值，ximax为当前地类对应的处理后的后期影像中各像元灰度的最大值，ximin为当前地类对应的处理后的后期影像中各像元灰度的最小值，yimax为当前地类对应的处理后的前期影像中各像元灰度的最大值，yimin为当前地类对应的处理后的前期影像中各像元灰度的最小值。
4.根据权利要求1所述的多时相林业遥感影像变化监测方法，其特征在于，步骤S4)中将归一化后的前、后两期影像中的对应的像元一一进行差值计算的函数表达式如下：
Figure FDA0002667001720000031
上式中，X′i是归一化后的后期影像中的像元i的灰度值，Yi'是归一化后的前期影像中的像元i的灰度值，e是预设的差异性阈值。
5.根据权利要求1所述的多时相林业遥感影像变化监测方法，其特征在于，步骤S5)中众数滤波的具体步骤包括：针对影像变化图B或B’，选取影像变化图中的一个像元作为当前像元，然后以当前像元为中心，设置一个像元为距离的矩形窗口，若矩形窗口中当前像元周边的其他像元中超过八分之五的相邻像元具有相同值，将这些相邻像元的值替换当前像元的值。
6.根据权利要求1所述的多时相林业遥感影像变化监测方法，其特征在于，步骤S6)中焦点统计的具体步骤包括：针对修复后的影像变化图B1或B1’，选取修复后的影像变化图中的一个像元作为当前像元，然后以当前像元为中心，设置两个像元为距离的矩形窗口，统计矩形窗口中所有像元的值的平均值，将所有像元的值的平均值作为当前像元的值。
7.根据权利要求1所述的多时相林业遥感影像变化监测方法，其特征在于，步骤S7)具体包括以下步骤：
S71)分别将影像变化图B和B’进行二值化，得到有植被覆盖的林地以及无植被覆盖的林地对应的初步变化监测结果的二值图C0和C0’；
S72)分别对二值图C0和C0’矢量化，得到矢量图D0和D0’；
S73)将矢量化后的矢量图D0和D0’合并得到初步变化检测的矢量图D。
8.根据权利要求1所述的多时相林业遥感影像变化监测方法，其特征在于，步骤S9)具体包括以下步骤：
S91)将剔除破碎斑块后的初步变化检测的矢量图D输入深度神经网络模型；
S92)以初步变化检测的矢量图D为进一步判定对象范围，通过神经网络模型进行分类和筛选后剔除图中影像辐射误差造成的未发生变化的图斑，得到最终变化检测的矢量图D1。
9.根据权利要求1所述的多时相林业遥感影像变化监测方法，其特征在于，步骤S10)具体包括以下步骤：
S101)将最终变化检测的矢量图D1中图斑的所有面数据转化为点数据；
S102)根据点数据提取对应像元的灰度值；
S103)对所有被提取的像元的灰度值先通过箱线图分析来进行异常值检测，剔除异常值后再进行K均值聚类，将最小聚类中心的灰度值作为聚类分析后的灰度值，根据聚类分析后的灰度值与有植被覆盖的林地以及无植被覆盖的林地的灰度值的阈值的大小关系，分析变化原因。
10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有被编程或配置以实现权利要求1～9任一所述的多时相林业遥感影像变化监测方法的计算机程序。
CN202010921885.4A 2020-09-04 2020-09-04 多时相林业遥感影像变化监测方法 Active CN112101159B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010921885.4A CN112101159B (zh) 2020-09-04 2020-09-04 多时相林业遥感影像变化监测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010921885.4A CN112101159B (zh) 2020-09-04 2020-09-04 多时相林业遥感影像变化监测方法

Publications (2)

Publication Number Publication Date
CN112101159A CN112101159A (zh) 2020-12-18
CN112101159B true CN112101159B (zh) 2022-05-24

Family

ID=73757343

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010921885.4A Active CN112101159B (zh) 2020-09-04 2020-09-04 多时相林业遥感影像变化监测方法

Country Status (1)

Country Link
CN (1) CN112101159B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112927252B (zh) * 2021-04-12 2023-09-22 二十一世纪空间技术应用股份有限公司 一种新增建设用地监测方法及装置
CN113379620B (zh) * 2021-05-18 2023-10-27 中国资源卫星应用中心 一种光学遥感卫星影像云检测方法
CN113095303B (zh) * 2021-06-04 2021-09-28 成都数之联科技有限公司 模型训练方法及林地变化检测方法及***及装置及介质
CN113421273B (zh) * 2021-06-30 2022-02-25 中国气象科学研究院 林草搭配信息遥感提取方法及装置
CN114119575B (zh) * 2021-11-30 2022-07-19 二十一世纪空间技术应用股份有限公司 一种空间信息变化检测方法及***
CN114817616B (zh) * 2022-06-29 2022-09-16 四川省林业和草原调查规划院(四川省林业和草原生态环境监测中心) 一种森林蓄积量连续监测方法、***及其执行方法
CN115830448B (zh) * 2022-11-30 2024-02-09 广州市地质调查院(广州市地质环境监测中心) 一种基于多视角融合的遥感图像对比分析方法
CN116267382A (zh) * 2023-02-16 2023-06-23 西南大学 一种林下调水-保土-固碳的生态复合利用方法
CN117315491B (zh) * 2023-11-28 2024-02-20 泰安市绿威园林有限公司 一种面向林业工程建设的生态红线预警方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101976437B (zh) * 2010-09-29 2012-10-03 中国资源卫星应用中心 基于自适应阈值分割的高分辨率遥感影像变化检测方法
CN108846832B (zh) * 2018-05-30 2021-06-15 理大产学研基地(深圳)有限公司 一种基于多时相遥感影像与gis数据的变化检测方法及***
CN109657616A (zh) * 2018-12-19 2019-04-19 四川立维空间信息技术有限公司 一种遥感影像土地覆盖自动分类方法
CN110879969A (zh) * 2019-10-21 2020-03-13 北京中科锐景科技有限公司 一种基于遥感影像的林地变化检测方法

Also Published As

Publication number Publication date
CN112101159A (zh) 2020-12-18

Similar Documents

Publication Publication Date Title
CN112101159B (zh) 多时相林业遥感影像变化监测方法
CN111709379B (zh) 基于遥感影像的丘陵区柑橘种植地块监测方法及***
CN108846832B (zh) 一种基于多时相遥感影像与gis数据的变化检测方法及***
CN106951836B (zh) 基于先验阈值优化卷积神经网络的作物覆盖度提取方法
US8615133B2 (en) Process for enhancing images based on user input
CN112164062A (zh) 一种基于遥感时序分析的抛荒地信息提取方法及装置
CN111191628B (zh) 基于决策树与特征优化的遥感影像震害建筑物识别方法
CN104598908A (zh) 一种农作物叶部病害识别方法
CN111398176B (zh) 基于像元尺度特征的水体水色异常遥感识别方法和装置
CN108764284B (zh) 一种对松树病死木的高分辨率影像的分类去噪方法及***
CN108710864B (zh) 基于多维度识别及图像降噪处理的冬小麦遥感提取方法
Ji et al. In-field automatic detection of maize tassels using computer vision
CN114581764B (zh) 基于深度学习算法的地下结构裂纹病害判别方法
CN115272848A (zh) 多云多雾耕地保护区内建筑物智能变化检测方法
CN111476197A (zh) 基于多源卫星遥感影像油棕识别及面积提取的方法和***
CN116559111A (zh) 一种基于高光谱成像技术的高粱品种识别方法
Dang-Ngoc et al. Citrus leaf disease detection and classification using hierarchical support vector machine
CN111046838A (zh) 一种湿地遥感信息的识别方法及装置
CN114596489A (zh) 一种高精度人类居住指数的多源遥感城市建成区提取方法
CN113033386B (zh) 一种基于高分辨率遥感影像的输电线路通道隐患识别方法及***
CN116310913B (zh) 基于小型无人机测量技术的自然资源调查监测方法及装置
CN107657246A (zh) 一种基于多尺度滤波建筑指数的遥感影像建筑物检测方法
Wu et al. A Dense Litchi Target Recognition Algorithm for Large Scenes
CN115512159A (zh) 面向对象的高分辨率遥感影像地表覆盖分类方法及***
CN115115948A (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
GR01 Patent grant
GR01 Patent grant