CN113837074B - 结合后验概率和空间邻域信息的遥感影像变化检测方法 - Google Patents

结合后验概率和空间邻域信息的遥感影像变化检测方法 Download PDF

Info

Publication number
CN113837074B
CN113837074B CN202111118251.6A CN202111118251A CN113837074B CN 113837074 B CN113837074 B CN 113837074B CN 202111118251 A CN202111118251 A CN 202111118251A CN 113837074 B CN113837074 B CN 113837074B
Authority
CN
China
Prior art keywords
surface coverage
change
neighborhood information
intensity
pixel point
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
CN202111118251.6A
Other languages
English (en)
Other versions
CN113837074A (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.)
Shandong Jianzhu University
Original Assignee
Shandong Jianzhu 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 Shandong Jianzhu University filed Critical Shandong Jianzhu University
Priority to CN202111118251.6A priority Critical patent/CN113837074B/zh
Publication of CN113837074A publication Critical patent/CN113837074A/zh
Application granted granted Critical
Publication of CN113837074B publication Critical patent/CN113837074B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/21Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214Generating training patterns; Bootstrap methods, e.g. bagging or boosting

Landscapes

  • Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Image Analysis (AREA)
  • Image Processing (AREA)

Abstract

一种结合后验概率和空间邻域信息的遥感影像变化检测方法,利用双时相遥感影像,获取各地表覆盖类型的训练样本和变化强度图像,计算地表覆盖类型变化强度的后验概率。考虑空间邻域信息,对后验概率进行计算,构建地表覆盖类型‑空间邻域信息的空间曲面。根据空间曲面对每个像素自适应地确定阈值,获取变化检测结果图像。对多时相遥感影像变化检测研究而言,能够有效减少农田伪变化,并识别变化强度图像中灰度值较低的变化区域,最大限度地减少了漏报和误报的可能性。

Description

结合后验概率和空间邻域信息的遥感影像变化检测方法
技术领域
本发明涉及多时相遥感影像变化检测领域,具体涉及一种结合后验概率和空间邻域信息的遥感影像变化检测方法。
背景技术
多时相遥感变化检测在资源管理、环境保护、地表动态监测等领域发挥着越来越重要的作用。直接比较法是一种常用且原理简单的变化检测方法,该方法由两个步骤组成:(1)计算变化强度图像;(2)获取变化检测结果。变化强度图像的计算主要是比较T1时刻和T2时刻遥感影像光谱值之间的差异程度。常用的变化强度图像计算方法有:差值法、比值法、变化向量分析法等。变化阈值选取方法是变化检测结果生成中一种常见的方法,也是直接影响检测结果的重要步骤之一。值得注意的是,变化阈值的选择是一个相对复杂的过程,受多种因素的影响。在变化检测中,不同的地表覆盖类型具有不同的变化幅度。例如,草原和林地的光谱变化值相对较小,而农田转为建筑用地的光谱变化值则相对较大。此外,不同的地表覆盖呈现不同的空间异质性,这取决于研究区域和遥感图像的分辨率。例如,当研究区域全是农田时,空间异质性则相对较低。所有这些因素都导致了变化阈值选取的难度。
在过去的三十年里,学者们提出了各种变化阈值选取方法,可以分为全局变化阈值选取和局部变化阈值选取。全局变化阈值选取方法通常根据灰度直方图分布情况,利用一个或多个指标将变化强度图划分为变化区域和未变化区域,例如,Otsu方法、Kapur方法、Kittler方法、expectation maximum(EM)方法,等等。尽管方便有效,但全局变化阈值选取方法是从整体上确定阈值,缺乏对地表覆盖类型和空间异质性的考虑。局部变化阈值选取方法通常设置一个合适的窗口(如:3×3),将变化强度图像划分为相同大小的子区域,并在子区域内进行阈值分割。这种方法在某种程度上考虑了空间异质性,但没有整合地表覆盖类型信息。值得注意的是,如果忽略了地表覆盖类型之间的光谱变化差异,阈值的确定容易产生过检和漏检的问题。近年来,一些学者提出在变化阈值选取方法中引入地表覆盖类型信息。然而大多数的方法对空间异质性的考虑还不够充分,不能有效解决复杂地表覆盖类型的变化阈值选取问题。因此,考虑地表覆盖类型和空间异质性对提高变化检测结果的准确性十分重要。
发明内容
本发明为了克服以上技术的不足,提供了一种结合后验概率和空间邻域信息的遥感影像变化检测方法。
本发明克服其技术问题所采用的技术方案是:
一种结合后验概率和空间邻域信息的遥感影像变化检测方法,包括如下步骤:
a)利用双时相遥感影像,获取各地表覆盖类型的训练样本和变化强度图像,计算地表覆盖类型变化强度的后验概率;
b)利用空间邻域信息,对地表覆盖类型变化强度的后验概率进行计算,构建地表覆盖类型-空间邻域信息的空间曲面;
c)根据地表覆盖类型-空间邻域信息的空间曲面,自适应地选取阈值,确定变化强度图像是否发生变化。
进一步的,步骤a)中利用变化向量分析法获取变化强度图像。
进一步的,步骤a)中通过公式计算得到地表覆盖类型变化强度的后验概率I,式中M为T1时刻和T2时刻的地表覆盖类型数,m是常数,m取值范围在1到M之间,/>为T1时刻第m个地表覆盖类型的类概率,/>为T2时刻的第m个地表覆盖类型的类概率,/>与/>通过支持向量机获得,/>为T1时刻的第m个地表覆盖类型,/>为T2时刻的第m个地表覆盖类型,N为遥感影像的波段数,n是常数,n取值范围在1到N之间,/>为T1时刻第n个波段的遥感影像像素值,/>为T2时刻第n个波段的遥感影像像素值。
进一步的,步骤b)中通过公式计算得到对应像素点p的地表覆盖类型-空间邻域信息的空间曲面值/>式中Ibf为地表覆盖类型-空间邻域信息的空间曲面,式中/>为对应像素点p的规范化权重之和,s为空间域,r为范围域,/>为空间域的高斯递减函数,/>为范围域的高斯递减函数,p为领域中心的像素点,q为与像素点p邻近的像素点,||p-q||为像素点p与像素点q之间的欧式距离,Ip为像素点p的强度值,Iq为像素点q的强度值,|Ip-Iq|为Ip与Iq之间距离的度量。
进一步的,步骤c)中通过公式
计算得到对应像素点p的阈值Tp,式中κ为提高阈值可区分度的系数,/> 为地表覆盖类型-空间邻域信息的空间曲面的标准差,σCMM为变化强度图像的标准差,/>为地表覆盖类型-空间邻域信息的空间曲面的均值,μCMM为变化强度图像的均值,当变化强度图像的像素点大于等于对应的阈值Tp时,则该像素点认定为变化,当变化强度图像的像素点小于对应的阈值Tp时,则该像素点认定为不变化,得到最终的变化检测结果。本发明的有益效果是:利用双时相遥感影像,获取各地表覆盖类型的训练样本和变化强度图像,计算地表覆盖类型变化强度的后验概率。考虑空间邻域信息,对后验概率进行计算,构建地表覆盖类型-空间邻域信息的空间曲面。根据空间曲面对每个像素自适应地确定阈值,获取变化检测结果图像。对多时相遥感影像变化检测研究而言,能够有效减少农田伪变化,并识别变化强度图像中灰度值较低的变化区域,最大限度地减少了漏报和误报的可能性。
附图说明
图1为本发明的方法流程图;
图2为本发明得到的参考变化图,图中(a)为2013年遥感影像,(b)为2018年遥感影像,(c)为参考变化图;
图3为本发明的变化检测结果图;
图4为本发明与传统方法的变化检测结果比较图。
具体实施方式
下面结合附图1、附图2对本发明做进一步说明。
如附图1所示,一种结合后验概率和空间邻域信息的遥感影像变化检测方法,包括如下步骤:
a)利用双时相遥感影像,获取各地表覆盖类型的训练样本和变化强度图像,计算地表覆盖类型变化强度的后验概率;
b)利用空间邻域信息,对地表覆盖类型变化强度的后验概率进行计算,构建地表覆盖类型-空间邻域信息的空间曲面;
c)根据地表覆盖类型-空间邻域信息的空间曲面,自适应地选取阈值,确定变化强度图像是否发生变化。
如附图2所示,我们选取2013年5月21日和2018年5月3日的Landsat8-OperationalLand Imager(OLI)遥感影像,研究区域图像大小为376×350像素,位于中国山东省乐陵市。
利用双时相遥感影像,获取各地表覆盖类型的训练样本和变化强度图像,计算地表覆盖类型变化强度的后验概率。考虑空间邻域信息,对后验概率进行计算,构建地表覆盖类型-空间邻域信息的空间曲面。根据空间曲面对每个像素自适应地确定阈值,获取变化检测结果图像。对多时相遥感影像变化检测研究而言,能够有效减少农田伪变化,并识别变化强度图像中灰度值较低的变化区域,最大限度地减少了漏报和误报的可能性。
实施例1:
步骤a)中利用变化向量分析(CVA)法获取变化强度图像。
实施例2:
步骤a)中通过公式计算得到地表覆盖类型变化强度的后验概率I,式中M为T1时刻和T2时刻的地表覆盖类型数,m是常数,m取值范围在1到M之间,/>为T1时刻第m个地表覆盖类型的类概率,/>为T2时刻的第m个地表覆盖类型的类概率,/>与/>通过支持向量机获得,/>为T1时刻的第m个地表覆盖类型,/>为T2时刻的第m个地表覆盖类型,N为遥感影像的波段数,n是常数,n取值范围在1到N之间,/>为T1时刻第n个波段的遥感影像像素值,/>为T2时刻第n个波段的遥感影像像素值。根据上述公式计算获得研究区域的地表覆盖类型变化强度的后验概率。
实施例3:
通过双边滤波公式,对地表覆盖类型变化强度的后验概率进行空间邻域信息整合,进而构建地表覆盖类型-空间邻域信息的空间曲面。具体的步骤b)中通过公式计算得到对应像素点p的地表覆盖类型-空间邻域信息的空间曲面值/>式中Ibf为地表覆盖类型-空间邻域信息的空间曲面,式中/>为对应像素点p的规范化权重之和,/>s为空间域,r为范围域,/>为空间域的高斯递减函数,/>为范围域的高斯递减函数,p为领域中心的像素点,q为与像素点p邻近的像素点,||p-q||为像素点p与像素点q之间的欧式距离,Ip为像素点p的强度值,Iq为像素点q的强度值,|Ip-Iq|为Ip与Iq之间距离的度量。根据上述公式计算获得地表覆盖类型-空间邻域信息的空间曲面。
实施例4:
进一步的,步骤c)中,变化强度图像标准化为0-255,最终阈值由255减去地表覆盖类型-空间邻域信息的空间曲面乘以系数κ得到。因此,提出方法对变化强度图像中的变化区域能够获取相对较小的阈值,从而更容易检测到变化,对未变化区域能够获取相对较大的阈值,从而减少检测到伪变化的可能性。提出方法能更好地区分变化和未变化的区域,提高变化阈值选取的适应性程度。变化强度图像各像素点对应的最终阈值如公式:计算得到对应像素点p的阈值Tp,式中κ为提高阈值可区分度的系数,是变化强度图像和地表覆盖类型-空间邻域信息的空间曲面的离散系数之和,/> 为地表覆盖类型-空间邻域信息的空间曲面的标准差,σCMM为变化强度图像的标准差,/>为地表覆盖类型-空间邻域信息的空间曲面的均值,μCMM为变化强度图像的均值,当变化强度图像的像素点大于等于对应的阈值Tp时,则该像素点认定为变化,当变化强度图像的像素点小于对应的阈值Tp时,则该像素点认定为不变化,进而获得最终的变化检测结果(如图3所示)。
如附图4所示为了证明提出方法的有效性,与传统变化检测方法进行比较分析。对比方法包括:CVA-Otsu方法、CVA-Kittler方法、CVA-EM方法。利用总体精度(OA)、kappa系数、漏检率(OR)和过检率(CR)来评价变化检测结果。变化检测结果如下表所示:
CVA-Otsu CVA-Kittler CVA-EM 提出方法
OA(%) 83.86 83.86 82.22 91.62
kappa系数 0.65 0.65 0.62 0.82
CR(%) 35.11 35.11 40.44 18.52
OR(%) 1.36 1.36 0.14 0.48
通过实验证明可知,在研究区域中,CVA-Otsu方法和CVA-Kittler方法选取的阈值相同,存在较为明显的过检现象。CVA-EM方法也存在较多的过度检测,但该方法的漏检率最低。虽然提出方法的漏检率并非最低,但结合图4可知,提出方法最接近真实的变化。在多时相遥感影像变化检测研究中,本发明能够有效减少农田伪变化,并识别变化强度图像中灰度值较低的变化区域,最大限度地减少了漏报和误报的可能性。
最后应说明的是:以上所述仅为本发明的优选实施例而已,并不用于限制本发明,尽管参照前述实施例对本发明进行了详细的说明,对于本领域的技术人员来说,其依然可以对前述各实施例所记载的技术方案进行修改,或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (2)

1.一种结合后验概率和空间邻域信息的遥感影像变化检测方法,其特征在于,包括如下步骤:
a)利用双时相遥感影像,获取各地表覆盖类型的训练样本和变化强度图像,计算地表覆盖类型变化强度的后验概率;
b)利用空间邻域信息,对地表覆盖类型变化强度的后验概率进行计算,构建地表覆盖类型-空间邻域信息的空间曲面;
c)根据地表覆盖类型-空间邻域信息的空间曲面,自适应地选取阈值,确定变化强度图像是否发生变化;
步骤a)中通过公式计算得到地表覆盖类型变化强度的后验概率I,式中M为T1时刻和T2时刻的地表覆盖类型数,m是常数,m取值范围在1到M之间,/>为T1时刻第m个地表覆盖类型的类概率,/>为T2时刻的第m个地表覆盖类型的类概率,/>与/>通过支持向量机获得,/>为T1时刻的第m个地表覆盖类型,/>为T2时刻的第m个地表覆盖类型,N为遥感影像的波段数,n是常数,n取值范围在1到N之间,/>为T1时刻第n个波段的遥感影像像素值,/>为T2时刻第n个波段的遥感影像像素值;
步骤b)中通过公式计算得到对应像素点p的地表覆盖类型-空间邻域信息的空间曲面值/>式中Ibf为地表覆盖类型-空间邻域信息的空间曲面,式中/>为对应像素点p的规范化权重之和,/>s为空间域,r为范围域,/>为空间域的高斯递减函数,/>为范围域的高斯递减函数,p为领域中心的像素点,q为与像素点p邻近的像素点,p-q为像素点p与像素点q之间的欧式距离,Ip为像素点p的强度值,Iq为像素点q的强度值,Ip-Iq为Ip与Iq之间距离的度量;
步骤c)中通过公式
计算得到对应像素点p的阈值Tp,式中κ为提高阈值可区分度的系数,/> 为地表覆盖类型-空间邻域信息的空间曲面的标准差,σCMM为变化强度图像的标准差,/>为地表覆盖类型-空间邻域信息的空间曲面的均值,μCMM为变化强度图像的均值,当变化强度图像的像素点大于等于对应的阈值Tp时,则该像素点认定为变化,当变化强度图像的像素点小于对应的阈值Tp时,则该像素点认定为不变化,得到最终的变化检测结果。
2.根据权利要求1所述的结合后验概率和空间邻域信息的遥感影像变化检测方法,其特征在于:步骤a)中利用变化向量分析法获取变化强度图像。
CN202111118251.6A 2021-09-24 2021-09-24 结合后验概率和空间邻域信息的遥感影像变化检测方法 Active CN113837074B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111118251.6A CN113837074B (zh) 2021-09-24 2021-09-24 结合后验概率和空间邻域信息的遥感影像变化检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111118251.6A CN113837074B (zh) 2021-09-24 2021-09-24 结合后验概率和空间邻域信息的遥感影像变化检测方法

Publications (2)

Publication Number Publication Date
CN113837074A CN113837074A (zh) 2021-12-24
CN113837074B true CN113837074B (zh) 2023-08-11

Family

ID=78969738

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111118251.6A Active CN113837074B (zh) 2021-09-24 2021-09-24 结合后验概率和空间邻域信息的遥感影像变化检测方法

Country Status (1)

Country Link
CN (1) CN113837074B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115410096B (zh) * 2022-11-03 2023-01-24 成都国星宇航科技股份有限公司 卫星遥感影像多尺度融合变化检测方法、介质及电子装置

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101694719A (zh) * 2009-10-13 2010-04-14 西安电子科技大学 基于非参数密度估计的遥感图像变化检测方法
CN102169584A (zh) * 2011-05-28 2011-08-31 西安电子科技大学 基于分水岭和treelet的遥感图像变化检测方法
CN103488968A (zh) * 2012-06-14 2014-01-01 株式会社日立制作所 遥感图像的混合像素物质构成精细化分解装置及方法
CN111008644A (zh) * 2019-10-30 2020-04-14 西安电子科技大学 基于局部动态能量函数fcn-crf模型的生态变化监测方法
CN111259784A (zh) * 2020-01-14 2020-06-09 西安理工大学 基于迁移学习和主动学习的sar图像变化检测方法
CN112419266A (zh) * 2020-11-23 2021-02-26 山东建筑大学 一种地表覆盖类别约束的遥感影像变化检测方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101694719A (zh) * 2009-10-13 2010-04-14 西安电子科技大学 基于非参数密度估计的遥感图像变化检测方法
CN102169584A (zh) * 2011-05-28 2011-08-31 西安电子科技大学 基于分水岭和treelet的遥感图像变化检测方法
CN103488968A (zh) * 2012-06-14 2014-01-01 株式会社日立制作所 遥感图像的混合像素物质构成精细化分解装置及方法
CN111008644A (zh) * 2019-10-30 2020-04-14 西安电子科技大学 基于局部动态能量函数fcn-crf模型的生态变化监测方法
CN111259784A (zh) * 2020-01-14 2020-06-09 西安理工大学 基于迁移学习和主动学习的sar图像变化检测方法
CN112419266A (zh) * 2020-11-23 2021-02-26 山东建筑大学 一种地表覆盖类别约束的遥感影像变化检测方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Histogram thresholding for unsupervised change detection of remote sensing images;Swarnajyoti Patra et al.;International Journal of Remote Sensing;全文 *

Also Published As

Publication number Publication date
CN113837074A (zh) 2021-12-24

Similar Documents

Publication Publication Date Title
CN107292310B (zh) 一种圆形指针式表盘视觉定位及自动读数方法
CN102968798B (zh) 一种基于小波变换和otsu阈值的sar图像海陆分割方法
CN116188496B (zh) 基于土地利用类型的遥感影像自适应分割方法
CN111783583B (zh) 基于非局部均值算法的sar图像相干斑抑制方法
CN107123150A (zh) 全局颜色对比度检测和分割显著图的方法
CN108564085A (zh) 一种自动读取指针式仪表读数的方法
CN109741356A (zh) 一种亚像素边缘检测方法及***
CN108830899B (zh) 一种基于几何特征的快速自适应角点结合的圆心检测方法
CN116012364B (zh) Sar图像变化检测方法和装置
CN107256409B (zh) 基于sae和显著性检测的高分辨sar图像变化检测方法
CN113628234B (zh) 基于综合邻域信息的显著性极化sar图像变化检测方法
CN101976436A (zh) 一种基于差分图修正的像素级多聚焦图像融合方法
CN109461132A (zh) 基于特征点几何拓扑关系的sar图像自动配准方法
CN110070545B (zh) 一种城镇纹理特征密度自动提取城镇建成区的方法
CN113837074B (zh) 结合后验概率和空间邻域信息的遥感影像变化检测方法
CN107392887A (zh) 一种基于同质像素点转化的异质遥感图像变化检测方法
CN112308873A (zh) 多尺度Gabor小波PCA融合图像边缘检测方法
CN114677525B (zh) 一种基于二值图像处理的边缘检测方法
CN117911408B (zh) 一种道路路面施工质量检测方法及***
CN117541654B (zh) 一种高分辨率遥感图像细节增强方法
CN106940782A (zh) 基于变差函数的高分sar新增建设用地提取软件
CN108986083B (zh) 基于阈值优化的sar图像变化检测方法
CN107369163B (zh) 一种基于最佳熵双阈值分割的快速sar图像目标检测方法
CN116740579A (zh) 一种国土空间规划数据智能采集方法
CN111080647B (zh) 基于自适应滑动窗口滤波和fcm的sar图像分割方法

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