CN113129292A - 基于迭代马尔科夫的合成孔径雷达影像变化检测方法 - Google Patents

基于迭代马尔科夫的合成孔径雷达影像变化检测方法 Download PDF

Info

Publication number
CN113129292A
CN113129292A CN202110459965.7A CN202110459965A CN113129292A CN 113129292 A CN113129292 A CN 113129292A CN 202110459965 A CN202110459965 A CN 202110459965A CN 113129292 A CN113129292 A CN 113129292A
Authority
CN
China
Prior art keywords
image
binary image
binary
class
iterative
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.)
Granted
Application number
CN202110459965.7A
Other languages
English (en)
Other versions
CN113129292B (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.)
Shaanxi Normal University
Original Assignee
Shaanxi 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 Shaanxi Normal University filed Critical Shaanxi Normal University
Priority to CN202110459965.7A priority Critical patent/CN113129292B/zh
Publication of CN113129292A publication Critical patent/CN113129292A/zh
Application granted granted Critical
Publication of CN113129292B publication Critical patent/CN113129292B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/136—Segmentation; Edge detection involving thresholding
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/143—Segmentation; Edge detection involving probabilistic approaches, e.g. Markov random field [MRF] modelling
    • 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/28—Quantising the image, e.g. histogram thresholding for discrimination between background and foreground patterns
    • 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/20076—Probabilistic image processing
    • 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/20212—Image combination
    • G06T2207/20221—Image fusion; Image merging

Landscapes

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

Abstract

一种基于迭代马尔科夫的合成孔径雷达影像变化检测方法，由构建对数比值差异图D、确定先验概率和先验二值图、构建特征空间、确定特征数值的量化级、确定似然函数、确定变化类BC的后验概率、确定后验概率对应的二值图、马尔科夫迭代融合、输出最终结果步骤组成。本发明采用对同一地区、不同时刻的两幅合成孔径雷达影像进行变化检测，将一次贝叶斯过程中的先验概率图二值图和后验概率二值图进行迭代融合，更新，用后验概率代替先验概率进行马尔科夫融合，更新变化信息。融合了先验概率图与后验概率图，能检测到变化信息。本发明经对比仿真实验结果表明，本发明具有变化检测结果噪声更少，检测到的变化部分更接近真实结果优点。

Description

基于迭代马尔科夫的合成孔径雷达影像变化检测方法
技术领域
本发明属于遥感影像应用技术领域，具体地涉及到同一地区不同时刻两幅合成孔径雷达影像的变化检测方法。
背景技术
用多时相遥感影像探测来监测地球表面发生的变化是遥感技术最重要的应用之一。其中合成孔径雷达(SAR)因为具有全天候、全天时的特点，可方便地获得同一地区不同时刻的图像，被广泛应用于地震区域的定位和灾害评估、对农作物生长状况的监测、对土地使用的监测等等。SAR图像变化检测也成了它的一个重要的应用领域，通过对不同时期图像的比较分析，根据图像之间的差异得到所需要的变化信息。
由于SAR系统采用合成相干方式成像，乘性相干斑总是伴随SAR影像而存在。因此，如何有效地设计基于SAR影像的变化检测方法得到了越来越多学者的关注。
现有的SAR图像处理存在的主要技术问题是：SAR影像的乘性噪声对后续的构造差异图及差异图分类都有较大影像，会导致变化检测的精度降低，造成错分、误分的现象。
在图像处理技术领域，进行合成孔径雷达影像处理时，要充分考虑相干斑的影响，如何做到既能消除相干斑的影响又能保留真实变化信息不受到损失，是图像处理技术领域当前需迫切解决的一个技术问题。
发明内容
本发明所要解决的技术问题在于提供一种鲁棒性好、正确率高、分类精度高的基于迭代马尔科夫的合成孔径雷达影像变化检测方法。
解决上述技术问题所采用的技术方案由下述步骤组成：
(1)构建对数比值差异图D
用对数比值方法将同一地区不同时刻的合成孔径雷达影像I1和合成孔径雷达影像I2构建对数比值差异图D：
Figure BDA0003042053940000021
(2)确定先验概率和先验二值图
按下式确定先验概率
Figure BDA0003042053940000022
Figure BDA0003042053940000023
其中，l表示当前迭代次数，δ是参数，δ≠0，Df为归一差异图，T为将先验概率
Figure BDA0003042053940000024
经过最大类间方差法进行阈值分割得到的分割阈值，将先验概率
Figure BDA0003042053940000025
中大于阈值T的像素点置为1，小于阈值T的像素点置为0，得到先验二值图Tp。
(3)构建特征空间
采用3×3不重叠分块方法将差异图Df分成n个图像块，将图像块拉成列向量构建成自相关矩阵，并采用奇异值分解方法构建成K维特征空间：
Figure BDA0003042053940000026
其中，λj是自相关矩阵的第j个特征值，k为像素点个数，k∈[1,2，…，8]，n为有限的正整数。
(4)确定特征数值的量化级
将特征空间中的像素点用K-means聚类方法分为变化类C、非变化类UC以及混淆类MC三类，并通过下式将变化类C与非变化类UC分布间的K-L距离
DL(C,UC)＝K(C|UC)+K(UC|C)
其中，K(C|UC)表示从变化类C到非变化类UC的KL散度，K(UC|C)表示从非变化类UC到变化类C的KL散度，DL(C,UC)表示变化类C与非变化UC之间的K-L距离；将K-L距离最大化，获得变化类C与非变化类UC分布间相应特征值的量化级。
(5)确定似然函数
将量化级采用直方图方法得到先验二值图Tp中变化类BC的似然函数Pl(Fi|BC)和非变化类BUC的似然函数Pl(Fi|BUC)：
Figure BDA0003042053940000031
Figure BDA0003042053940000032
Figure BDA0003042053940000033
其中
Figure BDA0003042053940000034
表示像素i处的第k个特征属于变化类BC的条件概率，
Figure BDA0003042053940000035
表示像素i处的第k个特征属于非变化类BUC的条件概率，
Figure BDA0003042053940000036
是像素i处的第k个特征，ωk表示第k维特征的权重,
Figure BDA0003042053940000037
和
Figure BDA0003042053940000038
是第k维特征的变化部分和非变化部分的均值。
(6)确定变化类BC的后验概率
按下式确定差异图Df中各像素属于变化类BC的后验概率Pl(BC|Fo)：
Figure BDA0003042053940000041
(7)确定后验概率对应的二值图
通过最大类间方差阈值法对后验概率Pl(BC|Fi)进行二值化，得到后验概率对应的后验二值图Ta。
(8)马尔科夫迭代融合
利用能量最小化原则，采用迭代方式融合先验二值图Tp和后验二值图Ta，得到融合后二值图Y。
(9)输出最终结果
按下式确定先验概率
Figure BDA0003042053940000042
及先验二值图Tp：
Figure BDA0003042053940000043
Tp＝Y
重复步骤(5)-(8)，直至到最大迭代次数lmax，停止迭代，输出变化检测结果图Y。
在本发明的步骤(8)中，所述的迭代方式融合先验二值图Tp和后验二值图Ta的方法为：
1)用先验二值图Tp和后验二值图Ta，将能量函数U1最小化，得到初始化二值图B0，能量函数U1如下：
Figure BDA0003042053940000044
Figure BDA0003042053940000045
Figure BDA0003042053940000051
其中，bmn表示图像B0中像素点(m,n)的类标，AS(m,n)表示二值图A中像素点(m,n)的邻域，A(p,q)表示当前位置处属于AS(m,n)中的像素点，T表示计算二值图A所使用的阈值，xpq表示二值图A所对应差异图在(p,q)位置处的强度值，γ是一个调节参数，γ＞0。
2)将能量函数U2(Bt,Bt-1,Tp,Ta)最小化得到二值图Bt，能量函数U2如下：
Figure BDA0003042053940000052
Figure BDA0003042053940000053
其中，t表示马尔科夫融合迭代次数，β表示平衡因子，β＞0，
Figure BDA0003042053940000054
表示二值图Bt-1中(m,n)位置处的邻域，
Figure BDA0003042053940000055
表示当前位置处属于
Figure BDA0003042053940000056
中的像素点。
3)重复执行步骤2)，当二值图Bt-1和二值图Bt中不同标记的像素个数小于给定的阈值Td或者达到最大迭代次数tmax时，停止迭代，Td的取值小于总体像素点个数，二值图Bt为融合后二值图Y。
在本发明的(8)步骤的1)步骤中，所述的γ取值为(0，10]；在(8)步骤的2)步骤中，所述的β取值为(0，1]。
在本发明的(8)步骤的1)步骤中，所述的γ取值最佳为5。
在本发明的(8)步骤的3)步骤中，所述的tmax取值为500。
在本发明的步骤(9)中，所述的lmax取值最佳为3。
本发明的有益效果如下：
本发明采用基于迭代马尔科夫融合方法对同一地区、不同时刻的两幅合成孔径雷达影像进行变化检测，马尔科夫模型融合空间信息，可以将一次贝叶斯过程中的先验概率图二值图和后验概率二值图进行迭代融合，再通过迭代更新，利用后验概率代替先验概率继续进行马尔科夫融合，不断更新变化信息。本发明融合了先验概率图与后验概率图，能够最大限度检测到变化信息。本发明方法与PCA+K-means方法、TLFCM方法、MRF方法进行了对比仿真实验，实验结果表明，本发明具有变化检测结果噪声更少，检测到的变化部分更接近真实结果优点。
附图说明
图1是本发明实施例1的流程图。
图2是实施例1变化前的影像。
图3是实施例1变化后的影像。
图4是本发明对图2、图3进行变化检测的结果图。
图5是PCA+K-means方法对图2、图3进行变化检测的结果图。
图6是TLFCM方法对图2、图3进行变化检测的结果图。
图7是MRF方法对图2、图3进行变化检测的结果图。
图8是实施例1真实的变化图。
具体实施方式
下面结合附图和实施例对本发明进进一步详细说明，但本发明不限于下述实施方式。
实施例1
以1997年5月和1997年8月加拿大Ottawa地区的合成孔径雷达影像为例，本实施例的基于迭代马尔科夫的合成孔径雷达影像变化检测方法步骤如下(参见图1)：
(1)构建对数比值差异图D
用对数比值方法将同一地区不同时刻的合成孔径雷达影像I1和合成孔径雷达影像I2构建对数比值差异图D：
Figure BDA0003042053940000071
(2)确定先验概率和先验二值图
按下式确定先验概率
Figure BDA0003042053940000072
Figure BDA0003042053940000073
其中，l表示当前迭代次数，δ是参数，本实施例的δ取值为4×10-3，Df为归一差异图，T为将先验概率
Figure BDA0003042053940000075
经过最大类间方差法进行阈值分割得到的分割阈值，将先验概率
Figure BDA0003042053940000076
中大于阈值T的像素点置为1，小于阈值T的像素点置为0，得到先验二值图Tp。
(3)构建特征空间
采用3×3不重叠分块方法将差异图Df分成n个图像块，将图像块拉成列向量构建成自相关矩阵，并采用奇异值分解方法构建成K维特征空间：
Figure BDA0003042053940000074
其中λj是自相关矩阵的第j个特征值，k为像素点个数，k∈[1,2，…，8]，n取值为12168。
(4)确定特征数值的量化级
将特征空间中的像素点用K-means聚类方法分为变化类C、非变化类UC以及混淆类MC三类，并通过下式将变化类C与非变化类UC分布间的K-L距离
DL(C,UC)＝K(C|UC)+K(UC|C)
其中，K(C|UC)表示从变化类C到非变化类UC的KL散度，K(UC|C)表示从非变化类UC到变化类C的KL散度，DL(C,UC)表示变化类C与非变化UC之间的K-L距离；将K-L距离最大化，获得变化类C与非变化类UC分布间相应特征值的量化级。
(5)确定似然函数
将量化级采用直方图方法得到先验二值图Tp中变化类BC的似然函数Pl(Fi|BC)和非变化类BUC的似然函数Pl(Fi|BUC)：
Figure BDA0003042053940000081
Figure BDA0003042053940000082
Figure BDA0003042053940000083
其中
Figure BDA0003042053940000084
表示像素i处的第k个特征属于变化类BC的条件概率，
Figure BDA0003042053940000085
表示像素i处的第k个特征属于非变化类BUC的条件概率，
Figure BDA0003042053940000086
是像素i处的第k个特征，ωk表示第k维特征的权重,
Figure BDA0003042053940000087
和
Figure BDA0003042053940000088
是第k维特征的变化部分和非变化部分的均值。
(6)确定变化类BC的后验概率
按下式确定差异图Df中各像素属于变化类BC的后验概率Pl(BC|Fi)：：
Figure BDA0003042053940000089
(7)确定后验概率对应的二值图
通过最大类间方差阈值法对后验概率Pl(BC|Fi)进行二值化，得到后验概率对应的后验二值图Ta。
(8)马尔科夫迭代融合
利用能量最小化原则，采用迭代方式融合先验二值图Tp和后验二值图Ta，得到融合后二值图Y。
本实施例迭代方式融合先验二值图Tp和后验二值图Ta的方法为：
1)用先验二值图Tp和后验二值图Ta，将能量函数U1最小化，得到初始化二值图B0，能量函数U1如下：
Figure BDA0003042053940000091
Figure BDA0003042053940000092
Figure BDA0003042053940000093
其中，bmn表示图像B0中像素点(m,n)的类标，AS(m,m)表示二值图A中像素点(m,n)的邻域，A(p,q)表示当前位置处属于AS(m,n)中的像素点，T表示计算二值图A所使用的阈值，xpq表示二值图A所对应差异图在(p,q)位置处的强度值，γ是一个调节参数，γ＞0，本实施例的M取值为313，N取值为352，γ取值为5。
2)将能量函数U2(Bt,Bt-1,Tp,Ta)最小化得到二值图Bt，能量函数U2如下：
Figure BDA0003042053940000094
Figure BDA0003042053940000095
其中，t表示马尔科夫融合迭代次数，β表示平衡因子，β＞0，本实施例的β取值为0.5，
Figure BDA0003042053940000101
表示二值图Bt-1中(m,n)位置处的邻域，
Figure BDA0003042053940000102
表示当前位置处属于
Figure BDA0003042053940000103
中的像素点；
3)重复执行步骤(2)，二值图Bt-1和二值图Bt中不同标记的像素个数小于给定的阈值Td或者达到最大迭代次数tmax时，停止迭代，Td的取值小于总体像素点个数，二值图Bt为融合后二值图Y，本实施例的Td取值为2，tmax取值为500。
(9)输出最终结果
按下式确定先验概率
Figure BDA0003042053940000104
及先验二值图Tp：
Figure BDA0003042053940000105
Tp＝Y
重复步骤(5)-(8)，直至到最大迭代次数lmax，停止迭代，本实施例的lmax取值为3，输出变化检测结果图Y。
实施例2
以1997年5月和1997年8月加拿大Ottawa地区的合成孔径雷达影像为例，本实施例的基于迭代马尔科夫的合成孔径雷达影像变化检测方法步骤如下：
(1)-(7)步骤与实施例1相同。
(8)马尔科夫迭代融合
利用能量最小化原则，采用迭代方式融合先验二值图Tp和后验二值图Ta，得到融合后二值图Y。
本实施例迭代方式融合先验二值图Tp和后验二值图Ta的方法为：
1)用先验二值图Tp和后验二值图Ta，将能量函数U1最小化，得到初始化二值图B0，能量函数U1如下：
Figure BDA0003042053940000111
Figure BDA0003042053940000112
Figure BDA0003042053940000113
其中，bmn表示图像B0中像素点(m,n)的类标，AS(m,n)表示二值图A中像素点(m,n)的邻域，A(p,q)表示当前位置处属于AS(m,n)中的像素点，T表示计算二值图A所使用的阈值，xpq表示二值图A所对应差异图在(p,q)位置处的强度值，γ是一个调节参数，γ＞0，本实施例的M取值为313，N取值为352，γ取值为0.5。
2)将能量函数U2(Bt,Bt-1,Tp,Ta)最小化得到二值图Bt，能量函数U2如下：
Figure BDA0003042053940000114
Figure BDA0003042053940000115
其中，t表示马尔科夫融合迭代次数，β表示平衡因子，β＞0，本实施例的β取值为0.1，
Figure BDA0003042053940000116
表示二值图Bt-1中(m,n)位置处的邻域，
Figure BDA0003042053940000117
表示当前位置处属于
Figure BDA0003042053940000118
中的像素点；
3)重复执行步骤(2)，二值图Bt-1和二值图Bt中不同标记的像素个数小于给定的阈值Td或者达到最大迭代次数tmax时，停止迭代，Td的取值小于总体像素点个数，二值图Bt为融合后二值图Y，本实施例的Td取值为2,tmax取值为500。
其它步骤与实施例1相同，输出变化检测结果图Y。
实施例3
以1997年5月和1997年8月加拿大Ottawa地区的合成孔径雷达影像为例，本实施例的基于迭代马尔科夫的合成孔径雷达影像变化检测方法步骤如下：
(1)-(7)步骤与实施例1相同。
(8)马尔科夫迭代融合
利用能量最小化原则，采用迭代方式融合先验二值图Tp和后验二值图Ta，得到融合后二值图Y。
本实施例迭代方式融合先验二值图Tp和后验二值图Ta的方法为：
1)用先验二值图Tp和后验二值图Ta，将能量函数U1最小化，得到初始化二值图B0，能量函数U1如下：
Figure BDA0003042053940000121
Figure BDA0003042053940000122
Figure BDA0003042053940000123
其中，bmn表示图像B0中像素点(m,n)的类标，AS(m,n)表示二值图A中像素点(m,n)的邻域，A(p,q)表示当前位置处属于AS(m,n)中的像素点，T表示计算二值图A所使用的阈值，xpq表示二值图A所对应差异图在(p,q)位置处的强度值，γ是一个调节参数，γ＞0，本实施例的M取值为313，N取值为352，γ取值为10。
2)将能量函数U2(Bt,Bt-1,Tp,Ta)最小化得到二值图Bt，能量函数U2如下：
Figure BDA0003042053940000124
Figure BDA0003042053940000131
其中，t表示马尔科夫融合迭代次数，β表示平衡因子，β＞0，本实施例的β取值为1，
Figure BDA0003042053940000132
表示二值图Bt-1中(m,n)位置处的邻域，
Figure BDA0003042053940000133
表示当前位置处属于
Figure BDA0003042053940000134
中的像素点；
3)重复执行步骤(2)，二值图Bt-1和二值图Bt中不同标记的像素个数小于给定的阈值Td或者达到最大迭代次数tmax时，停止迭代，Td的取值小于总体像素点个数，二值图Bt为融合后二值图Y，本实施例的Td取值为2,tmax取值为500。
其它步骤与实施例1相同，输出变化检测结果图Y。
为了验证本发明的有益效果，发明人采用本发明实施例1的基于迭代马尔科夫的合成孔径雷达影像变化检测方法与PCA+K-means方法、TLFCM方法、MRF方法进行了对比仿真实验，实验时，选取1997年5月和1997年8月加拿大Ottawa地区的合成孔径雷达影像，图像大小为352×313，如图2、图3。本发明方法对图2、图3进行变化检测的结果图见图4，PCA+K-means方法对图2、图3进行变化检测的结果见图5，TLFCM方法对图2、图3进行变化检测的结果见图6，MRF方法对图2、图3进行变化检测的结果见图7。
4种方法的接测结果与真实的变化图像图8进行比较，按下式计算误检像素总数(OE)、检测准确率(Pcc)、Kappa系数(Kappa)：
OE＝FP+FN
Figure BDA0003042053940000135
Figure BDA0003042053940000136
其中，FP是非变化类被误检为变化类的像素个数，FN是变化类被误检为非变化类的像素个数，TP是正确检测为变化类像素个数，TN是正确检测为非变化类像素个数，N为像素总数，NC为真实的变化像素个数，NU为真实的非变化像素个数。
计算结果见表1、图4-7。
表1本发明方法和三种已知方法对两种图像的检测结果
检测方法 FP(个) FN(个) OE(个) Pcc(％) Kappa
PCA+K-means 956 1941 2897 97.371 0.89083
TLFCM 684 2638 3322 96.985 0.87147
MRF 1376 1499 2875 97.391 0.8941
本发明 1426 1270 2696 97.553 0.90142
由表1可见，本发明获得了OE、Pcc和kappa的最佳值，OE是误检像素个数总数，本发明OE数值较其他三种方法最小，说明本发明误检的像素个数最少。Pcc可以体现分类的正确性，本发明的Pcc较其他三种方法取得了最优值。Kappa系数是能更准确衡量分类准确度的参数，本发明具有较高的Kappa系数，说明在实际的多时相SAR图像的变化检测中具有较好的效果。先用马尔科夫迭代融合先验二值图和后验二值图，再通过贝叶斯继续更新先验二值图，在迭代的过程中融合空间信息具有较好的变化检测结果。本发明与现有的技术相比，考虑了空间信息，具有方法简单、鲁棒性高、准确度高等优点。
在图4中，白色部分为本发明检测出的发生变化区域，黑色部分为本发明检测出的未发生变化区域。与真实的变化部分参考图图8相比，本发明能够更加完整地检测到变化区域。
在图5中，白色部分为PCA+K-means方法检测出的发生变化区域，黑色部分为PCA+K-means方法检测出的未发生变化区域。与真实的变化部分参考图图8相比，PCA+K-means方法能检测出大部分变化区域，但仍存在漏检的情况。
在图6中，白色部分为TLFCM方法检测出的发生变化区域，黑色部分为TLFCM方法检测出的未发生变化区域。与真实的变化部分参考图图8相比，TLFCM方法检测的结果虽然去除了噪声，但也缺失了部分变化信息。
在图7中，白色部分为MRF方法检测出的发生变化区域，黑色部分为MRF方法检测出的未发生变化区域。与真实的变化部分参考图图8相比，MRF方法也存在检测变化区域不完整的情况，另外检测的结果中存在一些噪声，影响变化检测结果的准确性。
以上实验结果表明，实施例1的方法与三种已知方法相比，本发明具有变化检测结果噪声更少，检测到的变化部分更接近真实结果的优点。

Claims (6)

1.一种基于迭代马尔科夫的合成孔径雷达影像变化检测方法，其特征在于由下述步骤组成：
(1)构建对数比值差异图D
用对数比值方法将同一地区不同时刻的合成孔径雷达影像I1和合成孔径雷达影像I2构建对数比值差异图D：
Figure FDA0003042053930000011
(2)确定先验概率和先验二值图
按下式确定先验概率
Figure FDA0003042053930000012
Figure FDA0003042053930000013
其中，l表示当前迭代次数，δ是参数，δ≠0，Df为归一差异图，T为将先验概率
Figure FDA0003042053930000014
经过最大类间方差法进行阈值分割得到的分割阈值，将先验概率
Figure FDA0003042053930000015
中大于阈值T的像素点置为1，小于阈值T的像素点置为0，得到先验二值图Tp；
(3)构建特征空间
采用3×3不重叠分块方法将差异图Df分成n个图像块，将图像块拉成列向量构建成自相关矩阵，并采用奇异值分解方法构建成K维特征空间：
Figure FDA0003042053930000016
其中，λj是自相关矩阵的第j个特征值，k为像素点个数，k∈[1，2，...，8]，n为有限的正整数；
(4)确定特征数值的量化级
将特征空间中的像素点用K-means聚类方法分为变化类C、非变化类UC以及混淆类MC三类，并通过下式将变化类C与非变化类UC分布间的K-L距离
DL(C，UC)＝K(C|UC)+K(UC|C)
其中，K(C|UC)表示从变化类C到非变化类UC的KL散度，K(UC|C)表示从非变化类UC到变化类C的KL散度，DL(C，UC)表示变化类C与非变化UC之间的K-L距离；将K-L距离最大化，获得变化类C与非变化类UC分布间相应特征值的量化级；
(5)确定似然函数
将量化级采用直方图方法得到先验二值图Tp中变化类BC的似然函数Pl(Fi|BC)和非变化类BUC的似然函数Pl(Fi|BUC)：
Figure FDA0003042053930000021
Figure FDA0003042053930000022
Figure FDA0003042053930000023
其中
Figure FDA0003042053930000024
表示像素i处的第k个特征属于变化类BC的条件概率，
Figure FDA0003042053930000025
表示像素i处的第k个特征属于非变化类BUC的条件概率，
Figure FDA0003042053930000026
是像素i处的第k个特征，ωk表示第k维特征的权重，
Figure FDA0003042053930000027
和
Figure FDA0003042053930000028
是第k维特征的变化部分和非变化部分的均值；
(6)确定变化类BC的后验概率
按下式确定差异图Df中各像素属于变化类BC的后验概率Pl(BC|Fi)：：
Figure FDA0003042053930000029
(7)确定后验概率对应的二值图
通过最大类间方差阈值法对后验概率Pl(BC|Fi)进行二值化，得到后验概率对应的后验二值图Ta；
(8)马尔科夫迭代融合
利用能量最小化原则，采用迭代方式融合先验二值图Tp和后验二值图Ta，得到融合后二值图Y；
(9)输出最终结果
按下式确定先验概率
Figure FDA0003042053930000031
及先验二值图Tp：
Figure FDA0003042053930000032
Tp＝Y
重复步骤(5)-(8)，直至到最大迭代次数lmax，停止迭代，输出变化检测结果图Y。
2.根据权利要求1所述的基于迭代马尔科夫的合成孔径雷达影像变化检测方法，其特征在于在步骤(8)中，所述的迭代方式融合先验二值图Tp和后验二值图Ta的方法为：
1)用先验二值图Tp和后验二值图Ta，将能量函数U1最小化，得到初始化二值图B0，能量函数U1如下：
Figure FDA0003042053930000033
Figure FDA0003042053930000034
Figure FDA0003042053930000035
其中，bmn表示图像B0中像素点(m，n)的类标，AS(m，n)表示二值图A中像素点(m，n)的邻域，A(p，q)表示当前位置处属于AS(m，n)中的像素点，T表示计算二值图A所使用的阈值，xpq表示二值图A所对应差异图在(p，q)位置处的强度值，γ是一个调节参数，γ＞0；
2)将能量函数U2(Bt，Bt-1，Tp，Ta)最小化得到二值图Bt，能量函数U2如下：
Figure FDA0003042053930000041
Figure FDA0003042053930000042
其中，t表示马尔科夫融合迭代次数，β表示平衡因子，β＞0，
Figure FDA0003042053930000043
表示二值图Bt-1中(m，n)位置处的邻域，
Figure FDA0003042053930000044
表示当前位置处属于
Figure FDA0003042053930000045
中的像素点；
3)重复执行步骤2)，当二值图Bt-1和二值图Bt中不同标记的像素个数小于给定的阈值Td或者达到最大迭代次数tmax时，停止迭代，Td的取值小于总体像素点个数，二值图Bt为融合后二值图Y。
3.根据权利要求2所述的基于迭代马尔科夫的合成孔径雷达影像变化检测方法，其特征在于：在(8)步骤的1)步骤中，所述的γ取值为(0，10]；在(8)步骤的2)步骤中，所述的β取值为(0，1]。
4.根据权利要求2或3所述的基于迭代马尔科夫的合成孔径雷达影像变化检测方法，其特征在于：在(8)步骤的1)步骤中，所述的γ取值为5。
5.根据权利要求2所述的基于迭代马尔科夫的合成孔径雷达影像变化检测方法，其特征在于：在(8)步骤的3)步骤中，所述的tmax取值为500。
6.根据权利要求1或所述的基于迭代马尔科夫的合成孔径雷达影像变化检测方法，其特征在于：在步骤(9)中，所述的lmax取值为3。
CN202110459965.7A 2021-04-27 2021-04-27 基于迭代马尔科夫的合成孔径雷达影像变化检测方法 Active CN113129292B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110459965.7A CN113129292B (zh) 2021-04-27 2021-04-27 基于迭代马尔科夫的合成孔径雷达影像变化检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110459965.7A CN113129292B (zh) 2021-04-27 2021-04-27 基于迭代马尔科夫的合成孔径雷达影像变化检测方法

Publications (2)

Publication Number Publication Date
CN113129292A true CN113129292A (zh) 2021-07-16
CN113129292B CN113129292B (zh) 2023-04-07

Family

ID=76780224

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110459965.7A Active CN113129292B (zh) 2021-04-27 2021-04-27 基于迭代马尔科夫的合成孔径雷达影像变化检测方法

Country Status (1)

Country Link
CN (1) CN113129292B (zh)

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101694719A (zh) * 2009-10-13 2010-04-14 西安电子科技大学 基于非参数密度估计的遥感图像变化检测方法
CN104751453A (zh) * 2015-03-11 2015-07-01 西安电子科技大学 基于cuda和平稳小波变换的sar图像变化检测方法
CN104778717A (zh) * 2015-05-05 2015-07-15 西安电子科技大学 基于导向差异图的sar图像变化检测方法
CN106203521A (zh) * 2016-07-15 2016-12-07 西安电子科技大学 基于差异图自步学习的sar图像变化检测方法
CN106713935A (zh) * 2017-01-09 2017-05-24 杭州电子科技大学 一种基于贝叶斯决策的hevc块划分快速方法
CN107368781A (zh) * 2017-06-09 2017-11-21 陕西师范大学 基于子空间划分的合成孔径雷达影像变化检测方法
CN107644413A (zh) * 2017-08-25 2018-01-30 西安电子科技大学 基于邻域比值和自步学习的sar图像变化区域检测方法
CN108764119A (zh) * 2018-05-24 2018-11-06 西安电子科技大学 基于迭代最大类间方差的sar图像变化检测方法
CN109697474A (zh) * 2018-12-30 2019-04-30 陕西师范大学 基于迭代贝叶斯的合成孔径雷达影像变化检测方法
US20200065968A1 (en) * 2018-08-24 2020-02-27 Ordnance Survey Limited Joint Deep Learning for Land Cover and Land Use Classification

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101694719A (zh) * 2009-10-13 2010-04-14 西安电子科技大学 基于非参数密度估计的遥感图像变化检测方法
CN104751453A (zh) * 2015-03-11 2015-07-01 西安电子科技大学 基于cuda和平稳小波变换的sar图像变化检测方法
CN104778717A (zh) * 2015-05-05 2015-07-15 西安电子科技大学 基于导向差异图的sar图像变化检测方法
CN106203521A (zh) * 2016-07-15 2016-12-07 西安电子科技大学 基于差异图自步学习的sar图像变化检测方法
CN106713935A (zh) * 2017-01-09 2017-05-24 杭州电子科技大学 一种基于贝叶斯决策的hevc块划分快速方法
CN107368781A (zh) * 2017-06-09 2017-11-21 陕西师范大学 基于子空间划分的合成孔径雷达影像变化检测方法
CN107644413A (zh) * 2017-08-25 2018-01-30 西安电子科技大学 基于邻域比值和自步学习的sar图像变化区域检测方法
CN108764119A (zh) * 2018-05-24 2018-11-06 西安电子科技大学 基于迭代最大类间方差的sar图像变化检测方法
US20200065968A1 (en) * 2018-08-24 2020-02-27 Ordnance Survey Limited Joint Deep Learning for Land Cover and Land Use Classification
CN109697474A (zh) * 2018-12-30 2019-04-30 陕西师范大学 基于迭代贝叶斯的合成孔径雷达影像变化检测方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
F. BOVOLO等: "《A detail-preserving scale-driven approach to change detection in multitemporal SAR images》", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》 *
GABRIELE MOSER等: "《Unsupervised Change Detection From Multichannel SAR Data by Markovian Data Fusion》", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》 *
王纬华等: "《基于高分SAR影像的城市用地变化检测》", 《测绘与空间地理信息》 *
申邵洪等: "《基于马尔科夫随机场的多时相SAR影像变化检测研究》", 《长江科学院院报》 *

Also Published As

Publication number Publication date
CN113129292B (zh) 2023-04-07

Similar Documents

Publication Publication Date Title
CN103353989B (zh) 基于先验和融合灰度与纹理特征的sar图像变化检测方法
CN113160274A (zh) 一种基于YOLOv4的改进DeepSort目标检测跟踪方法
US9105109B2 (en) Method for superpixel life cycle management
US7949162B2 (en) System and method for solid component evaluation in mixed ground glass nodules
CN108717694B (zh) 基于模糊c均值聚类的电阻抗层析成像图像质量评价方法
CN111445484B (zh) 一种基于图像级标注的工业图像异常区域像素级分割方法
CN108550131B (zh) 基于特征融合稀疏表示模型的sar图像车辆检测方法
CN112508963B (zh) 一种基于模糊c均值聚类的sar图像分割方法
Streekstra et al. Analysis of tubular structures in three-dimensional confocal images
CN110415339B (zh) 计算输入三维形体间的匹配关系的方法和装置
Lorette et al. Fully unsupervised fuzzy clustering with entropy criterion
CN112270285B (zh) 一种基于稀疏表示和胶囊网络的sar图像变化检测方法
CN109697474B (zh) 基于迭代贝叶斯的合成孔径雷达影像变化检测方法
Ayech et al. Image segmentation based on adaptive Fuzzy-C-Means clustering
CN113129292B (zh) 基于迭代马尔科夫的合成孔径雷达影像变化检测方法
Irving et al. Classification of targets in synthetic aperture radar imagery via quantized grayscale matching
CN108932520B (zh) 结合先验概率估计的sar影像水体概率制图方法
CN116258877A (zh) 土地利用场景相似度变化检测方法、装置、介质及设备
CN107256399B (zh) 一种基于Gamma分布超像素方法和基于超像素TMF的SAR图像海岸线检测方法
CN106709921B (zh) 一种基于空间Dirichlet混合模型的彩色图像分割方法
McCarter The Kernel Density Integral Transformation
Li et al. Crack Detection and Recognition Model of Parts Based on Machine Vision.
Peng et al. Superpixel-Based Urban Change Detection in SAR Images Using Optimal Transport Distance
Xu et al. Unsupervised Change Detection Based on Iterative Histogram Matching and Bayesian Decision of Thresholding
CN115620030B (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