CN110111300B - 一种图像变化检测方法 - Google Patents

一种图像变化检测方法 Download PDF

Info

Publication number
CN110111300B
CN110111300B CN201910205137.3A CN201910205137A CN110111300B CN 110111300 B CN110111300 B CN 110111300B CN 201910205137 A CN201910205137 A CN 201910205137A CN 110111300 B CN110111300 B CN 110111300B
Authority
CN
China
Prior art keywords
difference image
pixel
image
energy
classification result
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
CN201910205137.3A
Other languages
English (en)
Other versions
CN110111300A (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.)
Xidian University
Original Assignee
Xidian 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 Xidian University filed Critical Xidian University
Priority to CN201910205137.3A priority Critical patent/CN110111300B/zh
Publication of CN110111300A publication Critical patent/CN110111300A/zh
Application granted granted Critical
Publication of CN110111300B publication Critical patent/CN110111300B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/29—Graphical models, e.g. Bayesian networks
    • G06F18/295—Markov models or related models, e.g. semi-Markov models; Markov random fields; Networks embedding Markov models
    • 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
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • 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/20—Special algorithmic details
    • G06T2207/20081—Training; Learning

Landscapes

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

Abstract

本发明公开了一种图像变化检测方法，包括：获取两时相图像，包括：第一图像和第二图像；利用两时相图像得到第一差异图像，并对第一差异图像训练后得到第二差异图像；利用MRF对第二差异图像分类后得到初始分类结果，其中，初始分类结果包括概率值和第一分类结果；将概率值输入一元势函数得到第一能量，将第一分类结果输入二元势函数得到第二能量，并对第一能量、第二能量求和得到总能量；利用CRF方法对总能量分割后得到第二分类结果；利用迭代法对第二分类结果循环得到最终结果。本发明对第一差异图像训练后得到第二差异图像，将第二差异图像作为MRF随机场的初始检测图像，抑制了噪声的干扰，且MRF随机场结合CRF随机场，提升了检测结果的准确度。

Description

一种图像变化检测方法
技术领域
本发明属于图像处理技术领域，具体涉及一种图像变化检测方法。
背景技术
遥感图像变化检测主要通过对某一相同区域不同时刻所对应的遥感图像进行分析，从而获取图像之间的差异，再利用相关技术，进一步的识别出地物状态信息的变化，随着科学技术的发展，变化检测方法有很多：阈值法、还有基于超像素和基于语义的方法；还有基于随机场的方法，其中：MRF(Markov Random Field，马尔科夫随机场)和CRF(Conditional Random Field，条件随机场)是目前较为常用的两个方法。基于MRF和CRF的方法是将图像看成一个随机场，随机场中每一类像素的变化只与其邻域像素的变化相关，MRF是根据联合概率分布对图像进行建模，CRF是根据条件概率分布对图像进行建模。
第一种方法，首先用FCN(Fully Convolutional Network)模型生成粗略的包含飞机的候选图，将该候选图作为下一步的多层MRF方法的初始分类标签，然后使用多层MRF去进行优化，生成RP(regional proposal，候选区域)最后使用基于CNN(ConvolutionalNeural Network，卷积神经网络)的分类器去获取最终的结果；第二种方法，首先将SAR图像分割成两个分平稳区域，利用两个区域和SVM去定义CRF模型中的UP(Unary Potential，一元势函数)；然后，使用ROEWA(Ratio Of Exponentially Weighted Averages，指数加权平均比)算子从对数比值图像中提取边缘信息来定义标记场中邻域位置的空间相关性，进而定义对应的PP(Pairwise Potential，二元势函数)，使用ICE(Iterative ConditionalEstimation，条件迭代估计)对此步中的参数进行估计；最后使用ICM(IterationConditional Model，条件迭代模型)迭代优化，获取最终的变化检测结果图。
由于第一种方法中的FCN是对各个像素进行分类，未能充分的考虑像素之间的空间相关性，第二种方法考虑的是局部邻域信息之间的关系，两种方法的检测结果不是很理想，即检测结果的准确度不高。
发明内容
为了解决现有技术中存在的上述问题，本发明提供了一种图像变化检测方法。本发明要解决的技术问题通过以下技术方案实现：
本发明提供了一种图像变化检测方法，包括：
获取两时相图像，包括：第一图像和第二图像；
利用所述两时相图像得到第一差异图像，并对所述第一差异图像训练后得到第二差异图像；
利用MRF对所述第二差异图像进行分类后得到初始分类结果，其中，所述初始分类结果包括概率值和第一分类结果；
将所述概率值输入一元势函数得到第一能量，将所述第一分类结果输入二元势函数得到第二能量，并对所述第一能量、所述第二能量求和得到总能量；
利用CRF方法对所述总能量分割后得到第二分类结果；
利用迭代法对所述第二分类结果进行循环得到最终结果。
在本发明的一个实施例中，所述第一差异图像包括第一对数差异图像和/或第一均值差异图像；
所述第二差异图像包括第二对数差异图像和/或第二均值差异图像。
在本发明的一个实施例中，利用所述两时相图像得到第一差异图像，包括：
分别计算所述第一图像、所述第二图像上第一相同位置的像素8邻域像素值的和得到第一像素和、第二像素和；
对所述第二像素和除以所述第一像素和的商求对数得到第一对数值，并将所述第一对数值作为差异图像第一相同位置上的像素值，以得到所述第一对数差异图像。
在本发明的一个实施例中，利用所述两时相图像得到第一差异图像，包括：
分别计算所述第一图像、所述第二图像上第一相同位置的像素8邻域像素值的和得到第一像素和、第二像素和；
分别计算所述第一像素和、所述第二像素和的均值，得到第一均值、第二均值；
分别计算所述第一均值与所述第二均值的比值，将二者中的最小值作为差异图像第一相同位置上的像素值，以得到所述第一均值差异图像。
在本发明的一个实施例中，利用所述两时相图像得到第一差异图像，并对所述第一差异图像训练后得到第二差异图像，包括：
对所述第一对数差异图像训练后得到所述第二对数差异图像；和/或
对所述第一均值差异图像训练后得到所述第二均值差异图像。
在本发明的一个实施例中，利用MRF对所述第二差异图像进行分类后得到初始分类结果，包括：
将所述第二差异图像进行初始化；
计算初始化后的所述第二差异图像的均值和方差，再通过所述均值和方差计算出能量函数；
利用迭代法对所述能量函数进行最小化求解后得到所述初始分类结果。
在本发明的一个实施例中，将所述第二差异图像进行初始化，包括：
计算所述第二差异图像的期望值；
将所述第二差异图像的期望值进行最大化求解后完成初始化。
在本发明的一个实施例中，计算初始化后的所述第二差异图像的均值和方差，再通过所述均值和方差计算出能量函数，包括：
分别计算所述第二差异图像的均值和方差；
通过所述第二差异图像的均值和方差计算出所述第二差异图像中像素的分类能量和像素先验能量，并将所述分类能量和所述像素先验能量求和得到所述能量函数。
在本发明的一个实施例中，所述CRF方法的计算公式如下：
Figure GDA0003020881180000041
上式中：
y表示像素对应的分类类别；
xd表示差异图像中对应的像素值；z表示归一化常数；
S表示第二差异图像中像素的位置集合；
ym表示第m个像素对应的分类；
wp表示惩罚系数；δi表示像素的8邻域；θ表示参数。
本发明的有益效果：
本发明对所述第一差异图像进行模型训练后得到第二差异图像，将所述第二差异图像作为MRF随机场的初始检测图像，在一定程度上抑制了噪声的干扰，并且MRF随机场结合CRF随机场，对图像的空间信息和光谱信息进行了合理的利用，提升了检测结果的准确度。
以下将结合附图及实施例对本发明做进一步详细说明。
附图说明
图1是本发明实施例提供的一种图像变化检测方法的流程示意图；
图2是本发明实施例提供的一种模型训练的流程示意图；
图3是本发明实施例提供的一种输出初始分类结果的流程示意图；
图4是本发明实施例提供的另一种模型训练的流程示意图；
图5是本发明实施例提供的另一种输出初始分类结果的流程示意图。
具体实施方式
下面结合具体实施例对本发明做进一步详细的描述，但本发明的实施方式不限于此。
实施例一
请参见图1，图1是本发明实施例提供的一种图像变化检测方法的流程示意图。
本发明提供一种图像变化检测方法，包括：
步骤一、获取两时相图像，包括：第一图像和第二图像；
具体地，由卫星获取两幅图像作为两时相图像，包括：第一图像和第二图像，分别命名为I1和I2。
步骤二、利用所述两时相图像得到第一差异图像，并对所述第一差异图像训练后得到第二差异图像；
对所述第一对数差异图像训练后得到所述第二对数差异图像；和/或
对所述第一均值差异图像训练后得到所述第二均值差异图像。
具体地，所述第一差异图像是指第一对数差异图像，所述第二差异图像是指第二对数差异图像。
具体地，分别计算所述第一图像、所述第二图像上第一相同位置的像素8邻域像素值的和得到第一像素和、第二像素和；对所述第二像素和除以所述第一像素和的商求对数得到第一对数值，并将所述第一对数值作为差异图像第一相同位置上的像素值，以得到所述第一对数差异图像。
进一步地，所述第一相同位置为所述两时相图像的所有相同位置，例如位置(i，j)，其中，i表示行数；j表示列数。
进一步地，分别计算两幅图像即所述两时相图像I1和I2中相同位置(i，j)上像素的8邻域像素值的和，分别计作S1和S2，然后对所述S1和S2取对数，即：
Figure GDA0003020881180000061
用这个对数值作为差异图像(i，j)位置上的像素值，对两幅图像相同位置上的所有像素均采取上述操作，就可以得到第一对数差异图像。
进一步地，对应公式如下所示：
Figure GDA0003020881180000062
上式中：
(i，j)表示图像中像素的坐标；i表示行数；j表示列数；
(k，l)表示图像像素8邻域内像素的坐标；Ωi,j表示像素邻域；
I1(k，l)表示图像I1中某一像素的8邻域中(k，l)位置对应的像素；
I2(k，l)表示图像I2中某一像素的8邻域中(k，l)位置对应的像素；
Dlog表示第一对数差异图像；
进一步地，所述8邻域是指(i，j)位置的上、下、左、右、左上、左下、右上、右下。
进一步地，利用编码器对所述第一对数差异图像进行训练后得到所述第二对数差异图像。
请参见图2，图2是本发明实施例提供的一种模型训练的流程示意图。
具体地，将所述第一对数差异图像作为训练模型的训练样本，利用均值方差计算模块计算出所述第一对数差异图像中每一个像素的均值和方差，再通过所述第一对数差异图像中每一个像素的均值和方差计算出每一个像素对应的正态分布，对变量进行采样，经过生成器对所述第一对数差异图像进行重构后生成所述第二对数差异图像。
其中，所述编码器可以采用降噪自编码器、栈式自编码器，不同编码器对应不同的训练网络。
优选地，本发明实施例采用变分自编码器对所述第一对数差异图像进行训练。
进一步地，所述第二对数差异图像中的分类结果图只包含两类，变化为一类，记作1，不变化为另一类，记作2；将分类结果进行显示的时候，需要将变化类的1转换为真值1，不变化类2转换为真值0，最终显示的图像为一个二值图像，即只包括0和1。
本发明利用变分自编码器对所述第一对数差异图像进行训练后得到所述第二对数差异图像，将所述第二对数差异图像作为下一步MRF随机场的初始检测图像，在一定程度上抑制了噪声的干扰，提升了最终检测结果的准确度。
步骤三、利用MRF对所述第二差异图像进行分类后得到初始分类结果，其中，所述初始分类结果包括概率值和第一分类结果；
请参见图3，图3是本发明实施例提供的一种输出初始分类结果的流程示意图。
具体地，将所述第二对数差异图像进行参数初始化；计算所述第二对数差异图像的均值和方差，再通过所述均值和方差计算出能量函数；利用迭代法对所述能量函数进行最小化求解后得到所述初始分类结果。
进一步地，所述迭代是用于对所述均值和方差进行数据更新，使所述第二对数差异图像的像素数据更接近平稳值。
具体地，将所述第二差异图像进行初始化，包括：计算所述第二差异图像的期望值；将所述第二差异图像期望值进行最大化求解后完成初始化。
具体地，使用EM(Expectation Maximization，期望最大化)算法获取初始的分割结果；
进一步地，所述EM算法分为两步，对应公式如下：
E步(求期望)：
Figure GDA0003020881180000081
M步(最大化)：
Figure GDA0003020881180000082
上式中：
Qm表示z的某种分布(例如：高斯分布)；
z(m)表示像素对应的分类；
x(m)表示第m个像素值；θ表示参数。
步骤3.2、MRF随机场求解；
具体地，计算初始化后的所述第二差异图像的均值和方差，再通过所述均值和方差计算出能量函数，包括：分别计算所述第二差异图像的均值和方差；通过所述第二差异图像的均值和方差计算出所述第二差异图像中像素的分类能量和像素先验能量，并将所述分类能量和所述像素先验能量求和得到所述能量函数。
进一步地，根据贝叶斯公式简化得到，对应公式如下：
L＝argmax P(L|X)＝argmax{P(L)P(X|L)} (4)
上式中：
P(L)表示图像中每个像素类别标签的先验概率；
P(X|L)表示图像中每一个像素属于分类标签的条件概率密度函数；
根据Hammersly-Clifford定理可知：
Figure GDA0003020881180000091
上式中：
z表示归一化常数；
exp表示自然对数e为底的指数函数。
通过上述公式(5)计算可得所述概率值；
因此，将求最大概率问题转换成求最小能量函数问题。
能量函数计算公式如下所示：
minL[U(L|X)]＝minL[U(X|L)+U(L)] (6)
上式中：
U(X|L)表示一个像素属于某一分类的能量，即根据已知的分类推算哪些像素属于这一类；
U(L)表示像素属于哪一个分类的先验能量；
U(X|L)，U(L)对应公式如下：
Figure GDA0003020881180000101
Figure GDA0003020881180000102
上式中：
μk表示第k类中像素对应的均值；
Figure GDA0003020881180000103
表示第k类中像素对应的方差；
β表示惩罚系数；Nm表示像素m的邻域；
e以自然对数e为底的指数函数；
σk表示标准差(方差开平方根后的值)；
xm表示第m个像素；
xn表示第n个像素对应的8邻域内的像素。
进一步地，ρ(l(xm),l(xn))表示Potts模型，具体公式如下：
Figure GDA0003020881180000104
上式中：
l(xm)表示第m个像素的分类；
l(xn)表示第n个像素对应的8邻域内的像素分类。
综上所述，将公式(9)带入公式(8)计算得出U(L)的先验能量，再将公式(7)和公式(8)带入公式(6)求出能量函数，即所述公式(6)的输出即能量函数。将所述公式(6)输出的能量函数带入公式(5)得到所述概率值。
步骤3.2、利用迭代法对所述能量函数进行最小化求解后得到第一分类结果；
具体地，利用ICM迭代法对所述能量函数进行最小化求解后得到第一分类结果。
进一步地，利用设置其迭代次数控制其循环迭代的程度，设置其迭代次数为大于20次，当完成迭代次数后，输出的分类结果即所述第一分类结果。
优选地，所述迭代次数为50次。
具体地，所述第一分类结果的分类结果图中只包含两类，变化为一类，记作1，不变化为另一类，记作2；将分类结果进行显示的时候，需要将变化类的1转换为真值1，不变化类2转换为真值0，最终显示的图像为一个二值图像，即只包括0和1。
进一步地，所述初始分类结果包括所述概率值和所述第一分类结果。
步骤四、将所述概率值输入一元势函数得到第一能量，将所述第一分类结果输入二元势函数得到第二能量，并对所述第一能量、所述第二能量求和得到总能量；
步骤4.1、计算第一能量；
具体地，将所述概率值输入到一元势函数得到第一能量，计算公式如下：
Figure GDA0003020881180000111
上式中：
ym表示第m个像素；L表示像素分类类别；l表示第几类；
x表示像素x；θ表示参数；
Figure GDA0003020881180000112
表示由上述公式(5)得到的概率值，是图像中像素隶属于某一类的概率；
Figure GDA0003020881180000113
表示图像中像素和对应的分类标签之间的关系，即第一能量。
具体地，将公式(5)得到的概率值带入上述公式(10)中得到第一能量。
步骤4.2、计算第二能量；
具体地，将所述第一分类结果输入到二元势函数得到第二能量，计算公式如下：
Figure GDA0003020881180000121
上式中：
d(xm,xn)表示图像中像素之间的欧式距离；
σ表示d(xm,xn)的均方差；
ym表示第m个像素的分类结果；
yn表示第n个像素对应的8邻域内像素的分类结果；
x表示像素x；θ表示参数；
Figure GDA0003020881180000122
表示图像中的像素、其对应的分类标签及其领域像素对应的分类标签之间的关系，即第二能量。
步骤4.3、计算总能量；
具体地，总能量为所述第一能量与所述第二能量之和，计算公式如下：
Figure GDA0003020881180000123
上式中：
Figure GDA0003020881180000124
表示总能量。
步骤五、利用CRF方法对所述总能量分割后得到第二分类结果；
具体地，所述CRF方法是直接对后验概率建模，求取图像对应的分类结果的一种方法，其对应的目标函数如下所示：
Figure GDA0003020881180000131
上式中：
y表示像素对应的分类类别，即1和2两个类；
xd表示差异图像中对应的像素值；z表示归一化常数；
exp表示以自然对数e为底的指数函数；x表示像素x；
S表示差异图像中像素的位置集合(即像素在图像中的第几行第几列这样的一个位置集合)；
ym表示第m个像素对应的分类；
wp表示惩罚系数(用来约束
Figure GDA0003020881180000132
)；
δm表示像素的8邻域；θ表示参数。
进一步地，将所述总能量带入上述公式(13)后得到所述第二分类结果。
步骤六、利用迭代法对所述第二分类结果进行循环得到最终结果；
具体地，利用ICM迭代法对所述第二分类结果进行最小化求解后得到最终结果。
进一步地，利用设置其迭代次数控制其循环迭代的程度，设置其迭代次数为大于20次，当完成迭代次数后，输出的分类结果即所述最终结果。
优选地，所述迭代次数为50次。
本发明的有益效果：
本发明对所述第一差异图像进行模型训练后得到第二差异图像，将所述第二差异图像作为MRF随机场的初始检测图像，在一定程度上抑制了噪声的干扰，并且MRF随机场结合CRF随机场，对图像的空间信息和光谱信息进行了合理的利用，提升了检测结果的准确度。
实施例二
请再次参见图1，本发明提供了另一种图像变化检测方法，包括：
步骤一、获取两时相图像，包括：第一图像和第二图像；
具体地，由卫星获取两幅图像作为两时相图像，包括：第一图像和第二图像，分别命名为I1和I2。
步骤二、利用所述两时相图像得到第一差异图像，并对所述第一差异图像训练后得到第二差异图像；
对所述第一对数差异图像训练后得到所述第二对数差异图像；和/或
对所述第一均值差异图像训练后得到所述第二均值差异图像。
具体地，所述第一差异图像是指第一均值差异图像，所述第二差异图像是指第二均值差异图像。
具体地，分别计算所述第一图像、所述第二图像上第一相同位置的像素8邻域像素值的和得到第一像素和、第二像素和；分别计算所述第一像素和、所述第二像素和的均值，得到第一均值、第二均值；分别计算所述第一均值与所述第二均值的比值，将二者中的最小值作为差异图像第一相同位置上的像素值，以得到所述第一均值差异图像。
进一步地，所述第一相同位置为所述两时相图像的所有相同位置，例如位置(i，j)，其中，i表示行数；j表示列数。
进一步地，分别计算两幅图像即所述两时相图像I1和I2中相同位置(i，j)上像素的8邻域像素值的和，分别计作S1和S2，然后求出二者对应的均值，再求出两个均值对应的比值，将二者中的最小值作为图像(i，j)位置上的像素值，重复上述过程即可得到所述第一均值差异图。
进一步地，对应公式如下所示：
Figure GDA0003020881180000151
上式中：
Figure GDA0003020881180000152
表示图像I1的邻域灰度均值；
Figure GDA0003020881180000153
表示图像I2的邻域灰度均值；
Dmean表示第一均值差异图像。
进一步地，所述8邻域是指(i，j)位置的上、下、左、右、左上、左下、右上、右下。
进一步地，利用编码器对所述第一均值差异图像进行训练后得到所述第二均值差异图像，也可以对所述实施例一中的第一对数差异图像和所述第一均值差异图像训练后得到所述第二均值差异图像。
请参见图4，图4是本发明实施例提供的另一种模型训练的流程示意图。
具体地，将所述第一均值差异图像作为训练模型的训练样本，利用均值方差计算模块计算出所述第一均值差异图像中每一个像素的均值和方差，再通过所述第一均值差异图像中每一个像素的均值和方差计算出每一个像素对应的正态分布，对变量进行采样，经过生成器对所述第一均值差异图像进行重构后生成所述第二均值差异图像。
其中，所述编码器可以采用降噪自编码器、栈式自编码器，不同编码器对应不同的训练网络。
优选地，本发明实施例采用变分自编码器对所述第一均值差异图像进行训练。
进一步地，所述第二均值差异图像中的分类结果图只包含两类，变化为一类，记作1，不变化为另一类，记作2；将分类结果进行显示的时候，需要将变化类的1转换为真值1，不变化类2转换为真值0，最终显示的图像为一个二值图像，即只包括0和1。
本发明利用变分自编码器对所述第一均值差异图像进行训练后得到所述第二均值差异图像，将所述第二均值差异图像作为下一步MRF随机场的初始检测图像，在一定程度上抑制了噪声的干扰，提升了最终检测结果的准确度。
步骤三、利用MRF对所述第二差异图像进行分类后得到初始分类结果，其中，所述初始分类结果包括概率值和第一分类结果；
请参见图5，图5是本发明实施例提供的另一种输出初始分类结果的流程示意图。
具体地，将所述第二均值差异图像进行参数初始化；计算所述第二均值差异图像的均值和方差，再通过所述均值和方差计算出能量函数；利用迭代法对所述能量函数进行最小化求解后得到所述初始分类结果。
进一步地，所述迭代是用于对所述均值和方差进行数据更新，使所述第二均值差异图像的像素数据更接近平稳值。
具体地，将所述第二差异图像进行初始化，包括：计算所述第二差异图像的期望值；将所述第二差异图像期望值进行最大化求解后完成初始化。
具体地，使用EM算法获取初始的分割结果；
进一步地，所述EM算法分为两步，对应公式如下：
E步(求期望)：
Qm(z(m))：＝p(z(m)|x(m)；θ) (15)
M步(最大化)：
Figure GDA0003020881180000161
上式中：
Qm表示z的某种分布(例如：高斯分布)；
z(m)表示像素对应的分类；
x(m)表示第m个像素值；
θ表示参数。
步骤3.2、MRF随机场求解；
具体地，计算初始化后的所述第二差异图像的均值和方差，再通过所述均值和方差计算出能量函数，包括：分别计算所述第二差异图像的均值和方差；通过所述第二差异图像的均值和方差计算出所述第二差异图像中像素的分类能量和像素先验能量，并将所述分类能量和所述像素先验能量求和得到所述能量函数。
具体地，根据贝叶斯公式简化得到，对应公式如下：
L＝argmax P(L|X)＝argmax{P(L)P(X|L)} (17)
上式中：
P(L)表示图像中每个像素类别标签的先验概率；
P(X|L)表示图像中每一个像素属于分类标签的条件概率密度函数。
根据Hammersly-Clifford定理可知：
Figure GDA0003020881180000171
上式中：
z表示归一化常数；
exp表示自然对数e为底的指数函数。
通过上述公式(18)计算可得所述概率值；
因此，将求最大概率问题转换成求最小能量函数问题。
能量函数计算公式如下所示：
minL[U(L|X)]＝minL[U(X|L)+U(L)] (19)
上式中：
U(X|L)表示一个像素属于某一分类的能量，即根据已知的分类推算哪些像素属于这一类；
U(L)表示像素属于哪一个分类的先验能量。
U(X|L)，U(L)对应公式如下：
Figure GDA0003020881180000181
Figure GDA0003020881180000182
上式中：
μk表示第k类中像素对应的均值；
Figure GDA0003020881180000183
表示第k类中像素对应的方差；
β表示惩罚系数；Nm表示像素m的邻域；
e以自然对数e为底的指数函数；
σk表示标准差(方差开平方根后的值)；
xm表示第m个像素；
xn表示第m个像素对应的8邻域内的像素。
进一步地，ρ(l(xm),l(xn))表示Potts模型，具体公式如下：
Figure GDA0003020881180000184
上式中：
l(xm)表示第m个像素的分类；
l(xn)表示第n个像素对应的8邻域内的像素分类。
综上所述，将公式(22)带入公式(21)计算得出U(L)的先验能量，再将公式(20)和公式(21)带入公式(19)求出能量函数，即所述公式(19)的输出即能量函数。
将所述公式(19)输出的能量函数带入公式(18)得到所述概率值。
步骤3.2、利用迭代法对所述能量函数进行最小化求解后得到第一分类结果；
具体地，利用ICM迭代法对所述能量函数进行最小化求解后得到第一分类结果。
进一步地，利用设置其迭代次数控制其循环迭代的程度，设置其迭代次数为大于20次，当完成迭代次数后，输出的分类结果即所述第一分类结果。
优选地，所述迭代次数为50次。
具体地，所述第一分类结果的分类结果图中只包含两类，变化为一类，记作1，不变化为另一类，记作2；将分类结果进行显示的时候，需要将变化类的1转换为真值1，不变化类2转换为真值0，最终显示的图像为一个二值图像，即只包括0和1。
进一步地，所述初始分类结果包括所述概率值和所述第一分类结果。
步骤四、将所述概率值输入一元势函数得到第一能量，将所述第一分类结果输入二元势函数得到第二能量，并对所述第一能量、所述第二能量求和得到总能量；
步骤4.1、计算第一能量；
具体地，将所述概率值输入到一元势函数得到第一能量，计算公式如下：
Figure GDA0003020881180000201
上式中：
ym表示第m个像素；
L表示像素分类类别；l表示第几类；
x表示像素x；θ表示参数；
Figure GDA0003020881180000202
表示由上述公式(18)得到的概率值，是图像中像素隶属于某一类的概率；
Figure GDA0003020881180000203
表示图像中像素和对应的分类标签之间的关系，即第一能量。
具体地，将公式(18)得到的概率值带入上述公式(23)中得到第一能量。
步骤4.2、计算第二能量；
具体地，将所述第一分类结果输入到二元势函数得到第二能量，计算公式如下：
Figure GDA0003020881180000204
上式中：
d(xm,xn)表示图像中像素之间的欧式距离；
σ表示d(xm,xn)的均方差；
ym表示第m个像素的分类结果；
yn表示第n个像素对应的8邻域内像素的分类结果；
x表示像素x；θ表示参数；
Figure GDA0003020881180000205
表示图像中的像素、其对应的分类标签及其领域像素对应的分类标签之间的关系，即第二能量。
步骤4.3、计算总能量；
具体地，总能量为所述第一能量与所述第二能量之和，计算公式如下：
Figure GDA0003020881180000211
上式中：
Figure GDA0003020881180000212
表示总能量。
步骤五、利用CRF方法对所述总能量分割后得到第二分类结果；
具体地，所述CRF方法是直接对后验概率建模，求取图像对应的分类结果的一种方法，其对应的目标函数如下所示：
(26)
上式中：
y表示像素对应的分类类别，即1和2两个类；
xd表示差异图像中对应的像素值；z表示归一化常数；
exp表示以自然对数为e底的指数函数；x表示像素x；
S表示差异图像中像素的位置集合(即像素在图像中的第几行第几列这样的一个位置集合)；
ym表示第m个像素对应的分类；
wp表示惩罚系数(用来约束
Figure GDA0003020881180000213
)；
δm表示像素的8邻域；θ表示参数。
进一步地，将所述总能量带入上述公式(26)后得到所述第二分类结果。
步骤六、利用迭代法对所述第二分类结果进行循环得到最终结果；
具体地，利用ICM迭代法对所述第二分类结果进行最小化求解后得到最终结果。
进一步地，利用设置其迭代次数控制其循环迭代的程度，设置其迭代次数为大于20次，当完成迭代次数后，输出的分类结果即所述最终结果。
优选地，所述迭代次数为50次。
本发明的有益效果：
本发明对所述第一差异图像进行模型训练后得到第二差异图像，将所述第二差异图像作为MRF随机场的初始检测图像，在一定程度上抑制了噪声的干扰，并且MRF随机场结合CRF随机场，对图像的空间信息和光谱信息进行了合理的利用，提升了检测结果的准确度。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (4)

1.一种图像变化检测方法，其特征在于，包括：
获取两时相图像，包括：第一图像和第二图像；
利用所述两时相图像得到第一差异图像，并对所述第一差异图像训练后得到第二差异图像，所述第一差异图像包括第一对数差异图像，所述第二差异图像包括第二对数差异图像；
利用所述两时相图像得到第一差异图像，包括：
分别计算所述第一图像、所述第二图像上第一相同位置的像素8邻域像素值的和得到第一像素和、第二像素和，所述第一相同位置为所述两时相图像的所有相同位置；
对所述第二像素和除以所述第一像素和的商求对数得到第一对数值，并将所述第一对数值作为差异图像第一相同位置上的像素值，以得到所述第一对数差异图像，所述第一对数差异图像的计算公式为：
Figure FDA0003020881170000011
上式中：
(i，j)表示图像中像素的坐标；i表示行数；j表示列数；
(k，l)表示图像像素8邻域内像素的坐标；Ωi,j表示像素邻域；
I1(k，l)表示图像I1中某一像素的8邻域中(k，l)位置对应的像素；
I2(k，l)表示图像I2中某一像素的8邻域中(k，l)位置对应的像素；
Dlog表示第一对数差异图像；
对所述第一差异图像训练后得到第二差异图像，包括：
将所述第一对数差异图像作为训练模型的训练样本，利用均值方差计算模块计算出所述第一对数差异图像中每一个像素的均值和方差，再通过所述第一对数差异图像中每一个像素的均值和方差计算出每一个像素对应的正态分布，对变量进行采样，经过生成器对所述第一对数差异图像进行重构后生成所述第二对数差异图像；
利用MRF对所述第二差异图像进行分类后得到初始分类结果，其中，所述初始分类结果包括概率值和第一分类结果，所述第一分类结果的分类结果图中只包含两类，变化为一类，不变化为另一类；
将所述概率值输入一元势函数得到第一能量，将所述第一分类结果输入二元势函数得到第二能量，并对所述第一能量、所述第二能量求和得到总能量，所述第一能量表示图像中像素和对应的分类标签之间的关系，所述第二能量表示图像中的像素、其对应的分类标签及其邻域像素对应的分类标签之间的关系；
利用CRF方法对所述总能量分割后得到第二分类结果，所述CRF方法的计算公式如下：
Figure FDA0003020881170000021
上式中：
y表示像素对应的分类类别；
xd表示差异图像中对应的像素值；Z(x；θ)表示归一化常数；
exp表示自然对数e为底的指数函数；
S表示第二差异图像中像素的位置集合；
ym表示第m个像素对应的分类结果；yn表示第m个像素对应的8邻域内像素的分类结果；
Figure FDA0003020881170000022
表示图像中像素和对应的分类标签之间的关系，即第一能量；
Figure FDA0003020881170000031
表示图像中的像素、其对应的分类标签及其领域像素对应的分类标签之间的关系，即第二能量；
wp表示惩罚系数；δi表示像素的8邻域；θ表示参数；x表示像素x；
利用ICM迭代法对所述第二分类结果进行最小化求解得到最终结果。
2.根据权利要求1所述的方法，其特征在于，利用MRF对所述第二差异图像进行分类后得到初始分类结果，包括：
将所述第二差异图像进行初始化；
计算初始化后的所述第二差异图像的均值和方差，再通过所述均值和方差计算出能量函数；
利用迭代法对所述能量函数进行最小化求解后得到所述初始分类结果。
3.根据权利要求2所述的方法，其特征在于，将所述第二差异图像进行初始化，包括：
计算所述第二差异图像的期望值；
将所述第二差异图像的期望值进行最大化求解后完成初始化。
4.根据权利要求2所述的方法，其特征在于，计算初始化后的所述第二差异图像的均值和方差，再通过所述均值和方差计算出能量函数，包括：
分别计算所述第二差异图像的均值和方差；
通过所述第二差异图像的均值和方差计算出所述第二差异图像中像素的分类能量和像素先验能量，并将所述分类能量和所述像素先验能量求和得到所述能量函数。
CN201910205137.3A 2019-03-18 2019-03-18 一种图像变化检测方法 Active CN110111300B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910205137.3A CN110111300B (zh) 2019-03-18 2019-03-18 一种图像变化检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910205137.3A CN110111300B (zh) 2019-03-18 2019-03-18 一种图像变化检测方法

Publications (2)

Publication Number Publication Date
CN110111300A CN110111300A (zh) 2019-08-09
CN110111300B true CN110111300B (zh) 2021-06-15

Family

ID=67484408

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910205137.3A Active CN110111300B (zh) 2019-03-18 2019-03-18 一种图像变化检测方法

Country Status (1)

Country Link
CN (1) CN110111300B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111008644B (zh) * 2019-10-30 2023-04-21 西安电子科技大学 基于局部动态能量函数fcn-crf模型的生态变化监测方法
CN112200137B (zh) * 2020-10-29 2022-11-25 内蒙古工业大学 图像的识别方法及相应装置、存储介质及电子设备
CN112308025B (zh) * 2020-11-23 2023-01-31 内蒙古工业大学 图像变化的识别方法、装置、存储介质及电子设备

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102169584A (zh) * 2011-05-28 2011-08-31 西安电子科技大学 基于分水岭和treelet的遥感图像变化检测方法
CN104680549A (zh) * 2015-03-24 2015-06-03 西安电子科技大学 基于高阶邻域tmf模型的sar图像变化检测方法
CN106127236A (zh) * 2016-06-17 2016-11-16 西安电子科技大学 基于狄利克雷mrf混合模型的极化sar图像分类方法
CN109191503A (zh) * 2018-08-23 2019-01-11 河海大学 基于条件随机场的遥感影像变化检测方法及系统

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7139409B2 (en) * 2000-09-06 2006-11-21 Siemens Corporate Research, Inc. Real-time crowd density estimation from video
CN102496154B (zh) * 2011-10-24 2013-10-30 华中科技大学 一种基于马尔科夫随机场的多时相遥感图像变化检测方法
CN107228867A (zh) * 2017-06-21 2017-10-03 同方威视技术股份有限公司 安检图像显示方法、设备和安检系统

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102169584A (zh) * 2011-05-28 2011-08-31 西安电子科技大学 基于分水岭和treelet的遥感图像变化检测方法
CN104680549A (zh) * 2015-03-24 2015-06-03 西安电子科技大学 基于高阶邻域tmf模型的sar图像变化检测方法
CN106127236A (zh) * 2016-06-17 2016-11-16 西安电子科技大学 基于狄利克雷mrf混合模型的极化sar图像分类方法
CN109191503A (zh) * 2018-08-23 2019-01-11 河海大学 基于条件随机场的遥感影像变化检测方法及系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
DSSRM级联分割的SAR图像变化检测;张建龙等;《遥感学报》;20170430(第04期);第2.3节 *
Unsupervised change detection based on hybrid conditional random field model for high spatial resolution remote sensing imagery;pengyuan lv et al.;《IEEE transactions on geosciences and remote sensing》;20180516;第56卷(第7期);第4002-4015页 *
基于空-谱先验条件随机场的高分辨率遥感影像变化检测方法;吕鹏远等;《南京信息工程大学学报（自然科学版）》;20180228;第10卷(第01期);第1-2节 *

Also Published As

Publication number Publication date
CN110111300A (zh) 2019-08-09

Similar Documents

Publication Publication Date Title
CA3091035C (en) Systems and methods for polygon object annotation and a method of training an object annotation system
CN109118564B (zh) 一种基于融合体素的三维点云标记方法和装置
CN108537102B (zh) 基于稀疏特征与条件随机场的高分辨sar图像分类方法
Liu et al. Learning converged propagations with deep prior ensemble for image enhancement
Wilson et al. Structural matching by discrete relaxation
CN105701502B (zh) 一种基于蒙特卡罗数据均衡的图像自动标注方法
CN110111300B (zh) 一种图像变化检测方法
CN112287839B (zh) 一种基于迁移学习的ssd红外图像行人检测方法
CN109033978B (zh) 一种基于纠错策略的cnn-svm混合模型手势识别方法
CN111931816A (zh) 一种视网膜图像平行处理方法及装置
CN111160553A (zh) 一种新的领域自适应学习方法
CN111008644B (zh) 基于局部动态能量函数fcn-crf模型的生态变化监测方法
CN107358172B (zh) 一种基于人脸朝向分类的人脸特征点初始化方法
CN106600610B (zh) 一种fcm图像分割方法及装置
Suzuki et al. Adversarial transformations for semi-supervised learning
CN115861333A (zh) 基于涂鸦标注的医学图像分割模型训练方法、装置及终端
CN115311449A (zh) 基于类重激活映射图的弱监督图像目标定位分析系统
CN113627240B (zh) 一种基于改进ssd学习模型的无人机树木种类识别方法
CN107292268A (zh) 快速脊波反卷积结构学习模型的sar图像语义分割方法
Katartzis et al. A hierarchical Markovian model for multiscale region-based classification of vector-valued images
CN113313179A (zh) 一种基于l2p范数鲁棒最小二乘法的噪声图像分类方法
CN109191503A (zh) 基于条件随机场的遥感影像变化检测方法及系统
CN104933410A (zh) 一种高光谱图像光谱域与空间域联合分类方法
CN106504260B (zh) 一种fcm图像分割方法和系统
CN115393631A (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