CN107481235A - 一种数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法 - Google Patents

一种数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法 Download PDF

Info

Publication number
CN107481235A
CN107481235A CN201710733914.2A CN201710733914A CN107481235A CN 107481235 A CN107481235 A CN 107481235A CN 201710733914 A CN201710733914 A CN 201710733914A CN 107481235 A CN107481235 A CN 107481235A
Authority
CN
China
Prior art keywords
mrow
msub
msup
msubsup
mfrac
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
CN201710733914.2A
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.)
Hohai University HHU
Original Assignee
Hohai University HHU
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 Hohai University HHU filed Critical Hohai University HHU
Priority to CN201710733914.2A priority Critical patent/CN107481235A/zh
Publication of CN107481235A publication Critical patent/CN107481235A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Quality & Reliability (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，首先计算多时相多光谱影像的差异影像，在此基础上，计算影像上每一个点的CST值，根据致信水平获取阈值，得到初步的变化检测结果，然后再对该初步结果进行众数滤波(嵌入空间信息)，并根据滤波结果重新计算非变化区域的均值和方差矩阵。重复上述过程直到检测结果没有变化为止。其中检测过程中置信水平的选择是通过伪训练样本集来选择的。在最优的置信水平基础上，获取最终的变化检测结果。本发明在迭代估计卡方变换的统计参数过程中，加入数学形态学来对检测结果进行空间约束，通用性好，且提高了检测精度。

Description

一种数学形态学滤波结合卡方变换的多时相遥感影像变化检 测方法
技术领域
本发明涉及一种数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，属于多时相光学遥感影像变化检测技术领域。
背景技术
随着多时相遥感数据的不断积累以及空间数据库的相继建立，如何从这些遥感数据中提取和检测变化信息已成为遥感科学和地理信息科学的重要研究课题。根据同一区域不同时相的遥感影像，可以提取城市、环境等动态变化的信息，为资源管理与规划、环境保护等部门提供科学决策的依据。
遥感影像的变化检测就是从不同时期的遥感数据中，定量地分析和确定地表变化的特征与过程。各国学者从不同的角度和应用研究提出了许多有效的检测算法，如变化矢量分析法(Change Vector Analysis，CVA)、基于Fuzzy C-means(FCM)的聚类方法等。其中，传统的基于卡方变换(Chi-Squared Transform,CST)的多时相光学遥感变化检测，先计算差异影像的均值和方差矩阵，然后再基于置信水平，确定变化检测的阈值，进而得到变化检测结果。该类技术中，常规的CST的不足是仅使用多时相高分辨率差异影像的光谱信息，没有利用空间信息。Shi等提出的基于空间约束的CST结合MRF模型(SCCSTMRF)的非监督变换检测方法通过众数(mode)滤波将空间信息嵌入到检测过程中，提高了检测的精度[AiyeShi,Chao Wang,Shaohong Shen,Fengchen Huang,and Zhenli Ma.Unsupervised changedetection of multispectral images based on spatial constraint chi-squaredtransform and Markov random field model.Journal of Applied Remote Sensing,2016,10(4),046028:1-18.]。该空间约束的CST方法的不足是(1)容易引起较大的虚检率。(2)计算较复杂，对CST的结果须要采用MRF模型进一步精炼。
发明内容
针对现有技术存在的不足，本发明目的是提供一种数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，在迭代估计卡方变换的统计参数过程中，加入数学形态学来对检测结果进行空间约束，通用性好，且提高了检测精度。
为了实现上述目的，本发明是通过如下的技术方案来实现：
本发明的一种数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，包括以下几个步骤：
步骤1：输入两时相的高分辨率光学遥感影像，分别记为：X1和X2；
步骤2：利用遥感软件对X1和X2进行影像配准，分为几何粗校正和几何精校正；
步骤3：利用多元变化检测方法对X1和X2进行辐射归一化校正；
步骤4：对输入的多时相高分辨率光学遥感影像分别计算(是现有的计算方法，即简单的波段相减，此处不再赘述)每一个波段的差异影像，最后将每一个波段的差异影像进行组合(就是简单的组合，为现有技术。比如将独立的3个波段影像，组合为一个3波段的单独影像，此处不再赘述)，获取多时相差异影像Dx，DX＝X1-X2；
步骤5：计算Dx的模值，记为XM；并利用Bayes原理，基于最大期望算法获取最优分割阈值T，设定Dx的模值动态范围，并依据动态范围的百分比确定伪训练样本集；
步骤6：设定置信水平1-α的搜索范围和搜索步长；
步骤7：按照迭代步长增加置信水平；
步骤8：针对设定的置信水平，计算非变化区域的均值矢量m和方差矩阵∑，并计算卡方值；
步骤9：对每一个置信水平计算相应的阈值，并根据阈值获取初步的检测结果M0；
步骤10：设定数学形态学滤波结构单元B的形状和大小，对步骤9的检测结果M0进行数学形态学开操作滤波，记滤波的结果为M1：
式中，表示图像的开操作；
步骤11：利用步骤10的结果M1，确定影像的变化区域和非变化区域；
步骤12：判读本次迭代的检测结果与前一次检测结果是否有改变，如果没有改变，回到步骤7；如果有改变，计算DX中相对应于非变化区域的均值矢量和方差矩阵，回到步骤8；
步骤13：判读置信水平是否达到最大值；如果是，终止迭代；如果否，回到步骤7；
步骤14：针对每一个置信水平，判断(判断方法为现有技术，精度的指标可以选为分类的正确率)伪训练样本集的精度，选择精度最高的置信水平，并输出对应的变化检测结果。
步骤2中，对于所述几何粗校正，利用ENVI4.8软件实现，具体操作步骤为：
步骤2.1：显示基准影像和待校正影像；
步骤2.2：采集地面控制点GCPs，所述地面控制点GCPs应均匀分布在整幅图像内，地面控制点GCPs的数目至少大于等于9；
步骤2.3：按照计算同名控制点的均方根误差，其中P1 i，分别是前一时相和后一时相影像的同名控制点，S表示同名控制点的数目；
步骤2.4：选择两时相影像的同名控制点间的匹配模型为二次多项式模型；
步骤2.5：采用双线性插值法进行重采样输出。
上述双线性插值法具体的方法如下：
若求未知函数f在点P＝(x,y)的f(x，y)值，假设已知函数f在Q11＝(x1,y1),Q12＝(x1,y2),Q21＝(x2,y1)及Q22＝(x2,y2)四个点的值；
如果选择一个坐标系统使得这四个点的坐标分别为(0,0)、(0,1)、(1,0)和(1,1)，那么双线性插值公式就可以表示为：
f(x,y)≈f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy (1)。
步骤2中，所述几何精校正具体的方法如下：
将经过所述几何粗校正的多光谱遥感影像数据，利用三角剖分法进行几何精校正；所述三角剖分法为，采用逐点插入法构建Delaunay三角网，对每一个三角形，利用其三个顶点的行列号与其对应的基准影像同名点的地理坐标，来确定(具体的确定方法为现有技术，此处不再赘述)该三角形内部的仿射变换模型参数，对待校正影像进行纠正，得到校正后的遥感影。
步骤3中，具体的方法如下：
首先，找到两期影像各波段亮度值的一个线性组合，得到变化信息增强的差异影像，通过分割阈值确定变化区域和未变化区域；
然后，通过未变化区域对应的两时相像元对的映射方程，完成相对辐射校正。
步骤5中，将|XM-T|≤δ的区域作为伪训练样本集，其中δ的选择为XM动态范围的15％，最优分割阈值T的计算过程及伪训练样本集的构建如下：
1)假设XM影像上未变化类ωn和变化类ωc服从如下的高斯分布，即：
未变化类的均值和方差为mn和σn，变化类的均值和方差为mc和σc，X表示差异影像的模值，ωi表示类别，可以为非变化类ωn，也可以为变化类ωc；
2)采用最大期望算法估计mn、σn、mc和σc这四个参数，
其中，I和J分别表示影像的行数和列数，t表示迭代次数，pt+1(ωn)表示第t+1次迭代后，类别ωn的先验概率，表示非变化类在t次迭代后的均值，表示非变化类在t+1次迭代后的方差，公式中的p为均值概率；
3)依据Bayes最小误差准则，求解变化矢量幅值图像XM的分割阈值T，在高斯分布的情况下，等同于求解下式：
4)根据最大期望算法所估计的阈值T，伪训练样本集的构建分为如下两个部分：
其中，未变化类伪训练集样本为：
变化类伪训练样本集为：
步骤6中，所述置信水平1-α的搜索范围为0.95-0.999，搜索步长为0.001。
步骤8中，第一次迭代时，将整个Dx视为非变化的，
Cij＝(xij-m)TΣ-1(xij-m)～χ2(b) (10)
其中，Cij表示差异影像(i，j)坐标点的卡方值，其服从自由度为b的卡方分布；xij表示差异影像(i，j)坐标点的矢量值；Σ-1表示方差矩阵的逆矩阵；b表示差异影像的波段数目，χ表示卡方分布。
步骤9中，利用下式计算阈值
当置信水平为1-α时，Cij的值大于的概率为α；如果α取值小(一般0.01-0.05)，则大于的Cij视为出界点outlier或者变化点，由此确定阈值为并根据该阈值获得初步的检测结果M0。
本发明与现有的技术相比具有以下优点：
(1)在基于CST的变化检测中，将数学形态学滤波引入到CST检测结果之后，使得使用CST变换具有空间约束能力；
(2)变化检测中，采用迭代的方法，估计非变化区域的均值和标准差，克服了仅仅利用整个多波段影像Dx计算均值和方差矩阵的不足；可以使得变化检测的结果更加可靠，也更加具有稳健性；
(3)CST检测中最优置信水平是基于伪训练样本集获取，使得该技术通用性好。
附图说明
图1是本发明基于数学形态学结合空间约束CST的多时相遥感影像变化检测方法的实现流程示意图；
图2(a)是本发明所采用的2007年1月的沙特阿拉伯Mina地区高分辨率IKONOS图像第3波段示意图；
图2(b)本发明所采用的2007年12月的沙特阿拉伯的Mina地区高分辨率IKONOS图像第3波段示意图；
图2(c)Ground Truth图像；
图3(a)EM-CVA算法检测结果图像；
图3(b)ICST算法检测结果图像；
图3(c)SCSTMRF算法的检测结果图像；
图3(d)本发明算法检测结果图像(结构单元为方形、尺寸为5)。
具体实施方式
为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。
参照图1，本发明的实现步骤如下：
步骤1：输入同一区域、不同时相的两幅高分辨率光学遥感影像，分别记为：X1和X2。
步骤2：利用ENVI遥感软件对X1和X2进行影像配准，分为粗校正和精校正两个步骤：
对于几何粗校正，利用ENVI4.8软件中的相关功能实现，具体操作步骤为：(1)显示基准影像和待校正影像；(2)采集地面控制点GCPs；GCPs应均匀分布在整幅图像内，GCPs的数目至少大于等于9。(3)计算误差；(4)选择多项式模型；(5)采用双线性插值进行重采样输出。
双线性插值法，若求未知函数f在点P＝(x,y)的值，假设我们已知函数f在Q11＝(x1,y1),Q12＝(x1,y2),Q21＝(x2,y1),及Q22＝(x2,y2)四个点的值。如果选择一个坐标系统使得这四个点的坐标分别为(0,0)、(0,1)、(1,0)和(1,1)，那么双线性插值公式就可以表示为：
f(x,y)≈f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy (1)
对于几何精校正，将经过几何粗校正的多光谱遥感影像数据，利用自动匹配与三角剖分法进行几何精校正。
三角剖分法为，采用逐点插入法构建Delaunay三角网，对每一个三角形，利用其三个顶点的行列号与其对应的基准影像同名点的地理坐标来确定该三角形内部的仿射变换模型参数，对待校正影像进行纠正，得到校正后的遥感影。
步骤3：利用多元变化检测(Multivariate Alteration Detection，MAD)方法对X1和X2进行辐射归一化校正。该方法首先找到两期影像各波段亮度值的一个线性组合，得到变化信息增强的差异影像，通过阈值确定变化和未变化区域，然后通过未变化区域对应的两时相像元对的映射方程，完成相对辐射校正。
步骤4：对输入的多时相高分辨率影像X1和X2，计算多时相差异影像DX
DX＝X1-X2 (2)
步骤5：计算DX的模值，计为XM，利用Bayes原理，并基于最大期望(Expectation-Maximization，EM)算法获取最优分割阈值T。将|XM-T|≤δ的区域作为伪训练样本集。其中δ的选择为XM动态范围的15％，最优分割阈值T的计算过程及伪训练样本集的构建如下：
3)假设XM影像上未变化类ωn和变化类ωc服从如下的高斯分布，即：
未变化类的均值和方差为mn和σn，变化类的均值和方差为mc和σc。
4)采用EM算法可估计mn、σn、mc和σc这四个参数，下面仅以未变化类的参数估计为例进行说明，变化类参数估计类似。
其中，I和J分别表示影像的行数和列数，t表示迭代次数。
3)依据Bayes最小误差准则，求解变化矢量幅值图像XM的分割阈值T，在高斯分布的情况下，等同于求解下式：
4)根据EM算法所估计的阈值T，伪训练样本集的构建分为如下两个部分：
其中未变化类伪训练集样本为：
变化类伪训练样本集为：
步骤6：给定置信水平1-α的搜索范围为0.95-0.999，搜索步长为0.001。
步骤7：按照迭代步长增加置信水平。
步骤8：计算非变化区域其均值矢量m和方差矩阵∑，并计算差异影像每一个点的卡方值(开始迭代时，将整个DX视为非变化的)：
Cij＝(xij-m)TΣ-1(xij-m)～χ2(b) (10)
其中，Cij表示差异影像(i，j)坐标点的卡方值，其服从自由度为b的卡方分布；xij表示差异影像(i，j)坐标点的矢量值；Σ-1表示方差矩阵的逆矩阵；b表示差异影像的波段数目。
步骤9：利用下式计算阈值
当置信水平为1-α时，Cij的值大于的概率为α。如果α取值较小，则大于的Cij可以视为出界点(outlier)或者变化点，由此确定阈值为并根据该阈值获得初步的检测结果M0。
步骤10：给定数学形态学滤波结构单元B的形状、大小，一般设置结构单元为方形结构，结构单元大小设置为5×5，对M0进行数学形态学开操作滤波，记滤波的结果为M1：
式中，符合表示图像的开操作。
步骤11：根据步骤10的结果，确定M1中的非变化区域。
步骤12：判读M1中的标志相较于M0(变化和非变化标记)是否有改变，如果没有改变，回到步骤7。否则，计算DX中相对应于非变化区域(步骤11的结果)的均值矢量和方差矩阵。回到步骤8。
步骤13：判读置信水平是否达到最大值0.999。如果是，终止迭代；如果否，回到步骤7。
步骤14：针对每一个置信水平，判断伪训练样本集的精度，选择精度最高的置信水平，在此基础上，输出对应的变化检测结果。
本发明的效果可通过以下实验结果与分析进一步说明：
本发明的实验数据为沙特阿拉伯的Mina地区的多时相IKNOS高分辨影像数据，图像大小为700×950，使用B1、B2和B3三个波段。为了验证本发明的有效性，将本发明变化检测方法与下述变化检测方法进行比对：
(1)基于CVA的EM方法(CVA-EM)[意大利的Bruzzone L.等在文章“Automaticanalysis of difference image for unsupervised change detection”(IEEETransactions on Geoscience and Remote Sensing,2000,38(3):1171-1182.)中所提的检测方法]。
(2)基于迭代的CST检测(ICST)方法[B.Desclée,P.Bogaert,and P.Defourny在文章“Forest change detection by statistical object-based method”(Remote Sensingof Environment,2006,102(1-2):1-12.)中所提的方法]
(3)基于空间约束的CST结合MRF模型的检测(SCCSTMRF)方法[Aiye Shi,ChaoWang,Shaohong Shen,Fengchen Huang,and Zhenli Ma.Unsupervised change detectionof multispectral images based on spatial constraint chi-squared transform andMarkov random field model.Journal of Applied Remote Sensing,2016,10(4),046028:1-18]。
(4)本发明方法。
检测性能用错检数FP、漏检数FN、总错误数OE和Kappa系数四个指标来衡量。FP、FN和OE越接近于0、Kappa系数越接近于1，表明变化检测方法的性能越好。检测结果如表1所示。
表1 Mina地区的多时相IKONOS影像变化检测结果比较
方法 FP FN OE k
CVA-EM 49065 1456 50521 0.556
ICST 11236 868 12104 0.850
SCCSTMRF 4654 1296 5950 0.920
本发明方法 2258 2227 4485 0.937
理想 0 0 0 1
由表1可见，本发明所提的检测方法性能优于其他三种检测方法，这表明本发明所提的变化检测方法是有效的。
图2(a)是Mina地区的前一时相多光谱IKONOS影像，图2(b)是Mina地区的后一时相多光谱IKONOS影像，图2(c)是变化检测的参考图。图3(a)是CVA-EM算法的变化检测结果，图3(b)是ICST算法的变化检测结果，图3(c)是SCCSTMRF算法的变化检测结果，图3(d)是本发明方法的变化检测结果。
以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (9)

1.一种数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，其特征在于，包括以下几个步骤：
步骤1：输入两时相的高分辨率光学遥感影像，分别记为：X1和X2；
步骤2：利用遥感软件对X1和X2进行影像配准，分为几何粗校正和几何精校正；
步骤3：利用多元变化检测方法对X1和X2进行辐射归一化校正；
步骤4：对输入的多时相高分辨率光学遥感影像分别计算每一个波段的差异影像，最后将每一个波段的差异影像进行组合，获取多时相差异影像Dx，DX＝X1-X2；
步骤5：计算Dx的模值，记为XM；并利用Bayes原理，基于最大期望算法获取最优分割阈值T，设定Dx的模值动态范围，并依据动态范围的百分比确定伪训练样本集；
步骤6：设定置信水平1-α的搜索范围和搜索步长；
步骤7：按照迭代步长增加置信水平；
步骤8：针对设定的置信水平，计算非变化区域的均值矢量m和方差矩阵∑，并计算卡方值；
步骤9：对每一个置信水平计算相应的阈值，并根据阈值获取初步的检测结果M0；
步骤10：设定数学形态学滤波结构单元B的形状和大小，对步骤9的检测结果M0进行数学形态学开操作滤波，记滤波的结果为M1：
式中，表示图像的开操作；
步骤11：利用步骤10的结果M1，确定影像的变化区域和非变化区域；
步骤12：判读本次迭代的检测结果与前一次检测结果是否有改变，如果没有改变，回到步骤7；如果有改变，计算DX中相对应于非变化区域的均值矢量和方差矩阵，回到步骤8；
步骤13：判读置信水平是否达到最大值；如果是，终止迭代；如果否，回到步骤7；
步骤14：针对每一个置信水平，判断伪训练样本集的精度，选择精度最高的置信水平，并输出对应的变化检测结果。
2.根据权利要求1所述的数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，其特征在于，步骤2中，对于所述几何粗校正，利用ENVI4.8软件实现，具体操作步骤为：
步骤2.1：显示基准影像和待校正影像；
步骤2.2：采集地面控制点GCPs，所述地面控制点GCPs应均匀分布在整幅图像内，地面控制点GCPs的数目至少大于等于9；
步骤2.3：按照计算同名控制点的均方根误差，其中P1 i，分别是前一时相和后一时相影像的同名控制点，S表示同名控制点的数目；
步骤2.4：选择两时相影像的同名控制点间的匹配模型为二次多项式模型；
步骤2.5：采用双线性插值法进行重采样输出。
3.根据权利要求2所述的数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，其特征在于，所述双线性插值法具体的方法如下：
若求未知函数f在点P＝(x,y)的f(x，y)值，假设已知函数f在Q11＝(x1,y1),Q12＝(x1,y2),Q21＝(x2,y1)及Q22＝(x2,y2)四个点的值；
如果选择一个坐标系统使得这四个点的坐标分别为(0,0)、(0,1)、(1,0)和(1,1)，那么双线性插值公式就可以表示为：
f(x,y)≈f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy(1)。
4.根据权利要求1所述的数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，其特征在于，步骤2中，所述几何精校正具体的方法如下：
将经过所述几何粗校正的多光谱遥感影像数据，利用三角剖分法进行几何精校正；
所述三角剖分法为，采用逐点插入法构建Delaunay三角网，对每一个三角形，利用其三个顶点的行列号与其对应的基准影像同名点的地理坐标，来确定该三角形内部的仿射变换模型参数，对待校正影像进行纠正，得到校正后的遥感影。
5.根据权利要求1所述的数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，其特征在于，步骤3中，具体的方法如下：
首先，找到两期影像各波段亮度值的一个线性组合，得到变化信息增强的差异影像，通过分割阈值确定变化区域和未变化区域；
然后，通过未变化区域对应的两时相像元对的映射方程，完成相对辐射校正。
6.根据权利要求1所述的数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，其特征在于，步骤5中，将|XM-T|≤δ的区域作为伪训练样本集，其中δ的选择为XM动态范围的15％，最优分割阈值T的计算过程及伪训练样本集的构建如下：
1)假设XM影像上未变化类ωn和变化类ωc服从如下的高斯分布，即：
<mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>X</mi> <mo>|</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <msqrt> <mrow> <mn>2</mn> <mi>&amp;pi;</mi> </mrow> </msqrt> <msub> <mi>&amp;sigma;</mi> <mi>i</mi> </msub> </mrow> </mfrac> <mi>exp</mi> <mo>{</mo> <mo>-</mo> <mfrac> <msup> <mrow> <mo>(</mo> <mi>X</mi> <mo>-</mo> <msub> <mi>m</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mrow> <mn>2</mn> <msubsup> <mi>&amp;sigma;</mi> <mi>i</mi> <mn>2</mn> </msubsup> </mrow> </mfrac> <mo>}</mo> <mo>,</mo> <msub> <mi>&amp;omega;</mi> <mi>i</mi> </msub> <mo>&amp;Element;</mo> <mo>{</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>,</mo> <msub> <mi>&amp;omega;</mi> <mi>c</mi> </msub> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
未变化类的均值和方差为mn和σn，变化类的均值和方差为mc和σc，X表示差异影像的模值，ωi表示类别，可以为非变化类ωn，也可以为变化类ωc；
2)采用最大期望算法估计mn、σn、mc和σc这四个参数，
<mrow> <msup> <mi>p</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>&amp;Element;</mo> <msub> <mi>X</mi> <mi>M</mi> </msub> </mrow> </munder> <mfrac> <mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <mi>X</mi> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>|</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <mi>X</mi> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mfrac> </mrow> <mrow> <mi>I</mi> <mi>J</mi> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <msubsup> <mi>m</mi> <mi>n</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msubsup> <mo>=</mo> <mfrac> <mrow> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>&amp;Element;</mo> <msub> <mi>X</mi> <mi>M</mi> </msub> </mrow> </munder> <mfrac> <mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <mi>X</mi> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>|</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <mi>X</mi> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mfrac> <mi>X</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> </mrow> <mrow> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>&amp;Element;</mo> <msub> <mi>X</mi> <mi>M</mi> </msub> </mrow> </munder> <mfrac> <mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <mi>X</mi> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>|</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <mi>X</mi> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mfrac> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
<mrow> <msup> <mrow> <mo>(</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>n</mi> <mn>2</mn> </msubsup> <mo>)</mo> </mrow> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <mfrac> <mrow> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>&amp;Element;</mo> <msub> <mi>X</mi> <mi>M</mi> </msub> </mrow> </munder> <mfrac> <mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <mi>X</mi> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>|</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <mi>X</mi> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mfrac> <mo>&amp;lsqb;</mo> <mi>X</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>-</mo> <msubsup> <mi>m</mi> <mi>n</mi> <mi>t</mi> </msubsup> <mo>&amp;rsqb;</mo> </mrow> <mrow> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>X</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>&amp;Element;</mo> <msub> <mi>X</mi> <mi>M</mi> </msub> </mrow> </munder> <mfrac> <mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <mi>X</mi> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>|</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> </mrow> <mrow> <msup> <mi>p</mi> <mi>t</mi> </msup> <mrow> <mo>(</mo> <mi>X</mi> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>)</mo> </mrow> </mfrac> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
其中，I和J分别表示影像的行数和列数，t表示迭代次数，pt+1(ωn)表示第t+1次迭代后，类别ωn的先验概率，表示非变化类在t次迭代后的均值，表示非变化类在t+1次迭代后的方差，公式中的p为均值概率；
3)依据Bayes最小误差准则，求解变化矢量幅值图像XM的分割阈值T，在高斯分布的情况下，等同于求解下式：
<mrow> <mo>(</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>n</mi> <mn>2</mn> </msubsup> <mo>-</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>c</mi> <mn>2</mn> </msubsup> <mo>)</mo> <msup> <mover> <mi>T</mi> <mo>^</mo> </mover> <mn>2</mn> </msup> <mo>+</mo> <mn>2</mn> <mo>(</mo> <msub> <mi>m</mi> <mi>n</mi> </msub> <msubsup> <mi>&amp;sigma;</mi> <mi>c</mi> <mn>2</mn> </msubsup> <mo>-</mo> <msub> <mi>m</mi> <mi>c</mi> </msub> <msubsup> <mi>&amp;sigma;</mi> <mi>n</mi> <mn>2</mn> </msubsup> <mo>)</mo> <mover> <mi>T</mi> <mo>^</mo> </mover> <mo>+</mo> <msubsup> <mi>m</mi> <mi>c</mi> <mn>2</mn> </msubsup> <msubsup> <mi>&amp;sigma;</mi> <mi>n</mi> <mn>2</mn> </msubsup> <mo>-</mo> <msubsup> <mi>m</mi> <mi>n</mi> <mn>2</mn> </msubsup> <msubsup> <mi>&amp;sigma;</mi> <mi>c</mi> <mn>2</mn> </msubsup> <mo>-</mo> <mn>2</mn> <msubsup> <mi>&amp;sigma;</mi> <mi>n</mi> <mn>2</mn> </msubsup> <msubsup> <mi>&amp;sigma;</mi> <mi>c</mi> <mn>2</mn> </msubsup> <mi>l</mi> <mi>n</mi> <mo>&amp;lsqb;</mo> <mfrac> <mrow> <msub> <mi>&amp;sigma;</mi> <mi>c</mi> </msub> <mi>p</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>n</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>&amp;sigma;</mi> <mi>n</mi> </msub> <mi>p</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>c</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>&amp;rsqb;</mo> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
4)根据最大期望算法所估计的阈值T，伪训练样本集的构建分为如下两个部分：
其中，未变化类伪训练集样本为：
变化类伪训练样本集为：
7.根据权利要求1所述的数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，其特征在于，步骤6中，所述置信水平1-α的搜索范围为0.95-0.999，搜索步长为0.001。
8.根据权利要求1所述的数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，其特征在于，步骤8中，第一次迭代时，将整个Dx视为非变化的，
Cij＝(xij-m)TΣ-1(xij-m)～χ2(b) (10)
其中，Cij表示差异影像(i，j)坐标点的卡方值，其服从自由度为b的卡方分布；xij表示差异影像(i，j)坐标点的矢量值；Σ-1表示方差矩阵的逆矩阵；b表示差异影像的波段数目，χ表示卡方分布。
9.根据权利要求1所述的数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法，其特征在于，步骤9中，利用下式计算阈值
<mrow> <mi>P</mi> <mrow> <mo>(</mo> <msub> <mi>C</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> </msub> <mo>&lt;</mo> <msubsup> <mi>&amp;chi;</mi> <mrow> <mn>1</mn> <mo>-</mo> <mi>&amp;alpha;</mi> </mrow> <mn>2</mn> </msubsup> <mo>(</mo> <mi>b</mi> <mo>)</mo> <mo>)</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>-</mo> <mi>&amp;alpha;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>
当置信水平为1-α时，Cij的值大于的概率为α；如果α取值小，则大于的Cij视为出界点outlier或者变化点，由此确定阈值为并根据该阈值获得初步的检测结果M0。
CN201710733914.2A 2017-08-24 2017-08-24 一种数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法 Pending CN107481235A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710733914.2A CN107481235A (zh) 2017-08-24 2017-08-24 一种数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710733914.2A CN107481235A (zh) 2017-08-24 2017-08-24 一种数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法

Publications (1)

Publication Number Publication Date
CN107481235A true CN107481235A (zh) 2017-12-15

Family

ID=60602458

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710733914.2A Pending CN107481235A (zh) 2017-08-24 2017-08-24 一种数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法

Country Status (1)

Country Link
CN (1) CN107481235A (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109191503A (zh) * 2018-08-23 2019-01-11 河海大学 基于条件随机场的遥感影像变化检测方法及系统
CN109242832A (zh) * 2018-08-23 2019-01-18 河海大学 一种多时相多光谱遥感影像变化检测方法及系统
CN109523516A (zh) * 2018-10-19 2019-03-26 中国科学院遥感与数字地球研究所 一种基于双重约束条件的对象级土地覆盖变化检测方法
CN110827330A (zh) * 2019-10-31 2020-02-21 河海大学 一种时序集成的多光谱遥感图像变化检测方法及系统
CN110837787A (zh) * 2019-10-31 2020-02-25 河海大学 一种三方生成对抗网络的多光谱遥感图像检测方法及系统
CN111127393A (zh) * 2019-11-14 2020-05-08 苏州中科天启遥感科技有限公司 雷达影像变化检测的样本制作方法及系统、存储介质、设备
CN111340792A (zh) * 2020-03-05 2020-06-26 宁波市测绘设计研究院 一种遥感影像变化检测方法
CN114936992A (zh) * 2022-07-25 2022-08-23 北京数慧时空信息技术有限公司 一种遥感影像可用域的建立方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105405133A (zh) * 2015-11-04 2016-03-16 河海大学 一种遥感影像变化检测方法
CN106447653A (zh) * 2016-09-09 2017-02-22 河海大学 基于空间约束的卡方变换的多时相遥感影像变化检测方法
CN106469452A (zh) * 2016-09-09 2017-03-01 河海大学 基于空间约束的卡方变换的多时相遥感影像变化检测方法
CN106709500A (zh) * 2015-11-13 2017-05-24 国网辽宁省电力有限公司检修分公司 一种图像特征匹配的方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105405133A (zh) * 2015-11-04 2016-03-16 河海大学 一种遥感影像变化检测方法
CN106709500A (zh) * 2015-11-13 2017-05-24 国网辽宁省电力有限公司检修分公司 一种图像特征匹配的方法
CN106447653A (zh) * 2016-09-09 2017-02-22 河海大学 基于空间约束的卡方变换的多时相遥感影像变化检测方法
CN106469452A (zh) * 2016-09-09 2017-03-01 河海大学 基于空间约束的卡方变换的多时相遥感影像变化检测方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
杨金刚: "基于数学形态学的遥感图像边缘信息提取技术研究", 《中国优秀硕士学位论文全文数据库 信息科技辑（月刊）》 *

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109191503A (zh) * 2018-08-23 2019-01-11 河海大学 基于条件随机场的遥感影像变化检测方法及系统
CN109242832A (zh) * 2018-08-23 2019-01-18 河海大学 一种多时相多光谱遥感影像变化检测方法及系统
CN109242832B (zh) * 2018-08-23 2021-08-27 河海大学 一种多时相多光谱遥感影像变化检测方法及系统
CN109523516A (zh) * 2018-10-19 2019-03-26 中国科学院遥感与数字地球研究所 一种基于双重约束条件的对象级土地覆盖变化检测方法
CN110827330A (zh) * 2019-10-31 2020-02-21 河海大学 一种时序集成的多光谱遥感图像变化检测方法及系统
CN110837787A (zh) * 2019-10-31 2020-02-25 河海大学 一种三方生成对抗网络的多光谱遥感图像检测方法及系统
CN110827330B (zh) * 2019-10-31 2022-08-12 河海大学 一种时序集成的多光谱遥感图像变化检测方法及系统
CN111127393A (zh) * 2019-11-14 2020-05-08 苏州中科天启遥感科技有限公司 雷达影像变化检测的样本制作方法及系统、存储介质、设备
CN111340792A (zh) * 2020-03-05 2020-06-26 宁波市测绘设计研究院 一种遥感影像变化检测方法
CN111340792B (zh) * 2020-03-05 2022-04-12 宁波市测绘和遥感技术研究院 一种遥感影像变化检测方法
CN114936992A (zh) * 2022-07-25 2022-08-23 北京数慧时空信息技术有限公司 一种遥感影像可用域的建立方法

Similar Documents

Publication Publication Date Title
CN107481235A (zh) 一种数学形态学滤波结合卡方变换的多时相遥感影像变化检测方法
Chen et al. An automated approach for updating land cover maps based on integrated change detection and classification methods
CN103400151B (zh) 一体化的光学遥感影像与gis自动配准与水体提取方法
CN105405133B (zh) 一种遥感影像变化检测方法
CN103810704B (zh) 基于支持向量机和判别随机场的sar图像变化检测方法
CN107067405B (zh) 基于尺度优选的遥感影像分割方法
CN109255781B (zh) 一种面向对象的多光谱高分辨率遥感影像变化检测方法
CN105551031B (zh) 基于fcm和证据理论的多时相遥感影像变化检测方法
CN105574534A (zh) 基于稀疏子空间聚类和低秩表示的显著性目标检测方法
CN107689055A (zh) 一种多时相遥感影像变化检测方法
CN105389817B (zh) 一种两时相遥感影像变化检测方法
CN106447653B (zh) 基于空间约束的卡方变换的多时相遥感影像变化检测方法
CN106295124A (zh) 利用多种图像检测技术综合分析基因子图相似概率量的方法
Huang et al. Automatic building change image quality assessment in high resolution remote sensing based on deep learning
Yue et al. Texture extraction for object-oriented classification of high spatial resolution remotely sensed images using a semivariogram
CN106469452B (zh) 一种遥感影像的变化检测处理方法
CN103679734A (zh) 基于svm和pde的有眼台风二维表面风场反演方法
CN102651132A (zh) 一种基于交叉视觉皮质模型的医学图像配准方法
CN102779346A (zh) 基于改进c-v模型的sar图像变化检测方法
CN110458192A (zh) 基于视觉显著性的高光谱遥感图像分类方法及系统
CN113723492A (zh) 一种改进主动深度学习的高光谱图像半监督分类方法及装置
CN106650571B (zh) 一种基于自适应卡方变换的多时相遥感影像变化检测方法
CN109300115B (zh) 一种面向对象的多光谱高分辨率遥感影像变化检测方法
CN109242832B (zh) 一种多时相多光谱遥感影像变化检测方法及系统
Durduran Automatic classification of high resolution land cover using a new data weighting procedure: The combination of k-means clustering algorithm and central tendency measures (KMC–CTM)

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
RJ01 Rejection of invention patent application after publication

Application publication date: 20171215

RJ01 Rejection of invention patent application after publication