CN113628234A - 基于综合邻域信息的显著性极化sar图像变化检测方法 - Google Patents

基于综合邻域信息的显著性极化sar图像变化检测方法 Download PDF

Info

Publication number
CN113628234A
CN113628234A CN202110936927.6A CN202110936927A CN113628234A CN 113628234 A CN113628234 A CN 113628234A CN 202110936927 A CN202110936927 A CN 202110936927A CN 113628234 A CN113628234 A CN 113628234A
Authority
CN
China
Prior art keywords
information
region
pixel
value
neighborhood
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
CN202110936927.6A
Other languages
English (en)
Other versions
CN113628234B (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 CN202110936927.6A priority Critical patent/CN113628234B/zh
Publication of CN113628234A publication Critical patent/CN113628234A/zh
Application granted granted Critical
Publication of CN113628234B publication Critical patent/CN113628234B/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/10—Segmentation; Edge detection
    • G06T7/13—Edge detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • 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/12—Edge-based segmentation
    • 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
    • G06T2207/10044—Radar image

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Remote Sensing (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Electromagnetism (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)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明提出了一种基于综合邻域信息的显著性极化SAR图像变化检测方法，用于解决现有技术中存在的检测总体精度与区域内部一致性较低的技术问题。实现步骤为：对多时相极化SAR图像进行配准；获取多时相极化SAR图像的初始边缘信息和变差系数信息；通过变差系数信息对初始边缘信息进行修正；通过修正的边缘信息对变差系数进行指导；设计平滑因子自适应于综合邻域权重信息的非局部均值滤波器；获取任意两时相的综合邻域信息差异图；获取基于综合邻域信息差异图的显著性二值检测图；获取综合邻域信息差异图的显著性差异图；利用阈值二值化方法对显著性差异图进行二值化分割，得到变化检测结果。

Description

基于综合邻域信息的显著性极化SAR图像变化检测方法
技术领域
本发明属于雷达图像处理技术领域，涉及一种极化SAR图像变化检测方法，具体涉及一种基于综合邻域信息和显著性修正的极化SAR图像变化检测方法，可用于滑坡地质灾害的监测以及城市化建设的调查等极化SAR图像变化检测任务。
背景技术
极化SAR通过交替的水平发射和接收以及垂直发射和接收4种组合方式来获得含有丰富极化特性的极化散射矩阵，以此来对不同的地物目标进行衡量。因此，相较于合成孔径雷达(Synthetic Aperture Radar，SAR)，极化合成孔径雷达(Polarimetric SyntheticAperture Radar，PolSAR)拥有更丰富的地物目标信息和更全面的极化信息。极化SAR图像变化检测是指利用极化SAR对同一地区不同时间的两幅图像之间是否存在变化进行判断，进而得到变化检测结果图，其在监测毁林、灾害评估、城市规划等邻域发挥着巨大作用。极化SAR图像变化检测主要分为有监督的极化SAR图像变化检测方法和无监督的极化SAR图像变化检测方法。近年来，无监督极化SAR图像变化检测的主流框架包含三个部分：1)预处理：包括图像配准、相干斑抑制等操作；2)构造变化检测特征量并生成差异图(DifferenceImage，DI)；3)提取变化区域，得到变化检测结果。由于极化SAR特有的成像特性，在获取极化SAR图像过程中会引入相干斑噪声，这限制了利用极化SAR数据进行变化检测的可操作性。针对极化SAR图像变化检测任务易受到噪声影响的问题，在变化检测前通常会进行适当的滤波预处理，但这也会在一定程度上破坏最终变化检测结果的局部细节信息。综上所述，极化SAR图像变化检测方法的难点是在抑制相干斑噪声的同时保持对局部边缘细节信息的保留。
为了解决上述问题，Mahdavi,S等人在其发表的论文“A PolSAR ChangeDetection Index Based on Neighborhood Information for Flood Mapping”(RemoteSensing,2019,11(16):1854-)中，公开了一种基于极化散射总功率比值与邻域变差系数权重因子的极化SAR变化检测方法。该方法在产生两时相极化SAR图像的差异图时，引入了一个基于极化散射总功率比值并顾及了邻域信息的极化权重因子，该权重因子通过邻域内标准差与均值的比值(即变差系数)来实现，可以用于衡量中心像素与其对应邻域内其它像素的重要程度，进而对差异图求解算子中的原值项与滤波项各自对应的权重进行自适应调节。该方法考虑了图像局部区域的纹理特性，进而可以针对不同区域有区别地对原始图像进行滤波，所以产生的差异图可以在抑制噪声的同时也保留了图像的边缘细节信息，在一定程度上提高了变化检测的总体精度与内部一致性。但是其存在的缺陷在于，该方法忽略了变差系数在分布上的不均匀性，这导致差异图求解算子中的滤波项权重始终处在一个较高范围内，产生的差异图因过度均值滤波而无法保留完整边缘信息，且该方法未考虑异质区域中包含信息的复杂性，在保留异质区域中局部边缘细节信息的同时也保留了局部微小斑点噪声，最终导致变化检测结果的总体精度与区域内部一致性较低。
发明内容
本发明的目的在于针对上述现有技术存在的不足，提出一种基于综合邻域信息和显著性修正的极化SAR图像变化检测方法，旨在提高变化检测的总体精度与区域内部一致性。
为实现上述目的，本发明采取的技术方案包括如下步骤：
(1)对同一地区M个时相的极化SAR图像进行配准：
以同一地区M个时相的极化SAR图像T＝{T1,T2,...,Tm,...,TM}中任意一个时相的极化SAR图像作为参考图像，并将其它时相的极化SAR图像与参考图像进行配准，得到配准后的极化SAR图像集合T′＝{T′1,T′2,...,T′m,...,T′M}，其中，M≥2，Tm表示像素大小为w×d的第m个时相的极化SAR图像，w≥200，d≥200，T’m表示Tm的配准结果；
(2)获取配准后的每个时相的极化SAR图像对应的初始边缘信息和变差系数信息：
(2a)对配准后的每个时相的极化SAR图像T’m进行极化分解，得到T′对应的极化SAR图像分解合成图集合P＝{P1,P2,...,Pm,...,PM}，并对每个时相的极化SAR图像的分解合成图Pm进行边缘检测，得到T′对应的初始边缘信息图集合E＝{E1,E2,...,Em,...,EM}，其中，Pm、Em表示Tm对应的分解合成图、初始边缘信息图；
(2b)计算配准后的每个时相的极化SAR图像T’m的Frobenius范数，得到T′对应的Frobenius范数图集合F＝{F1,F2,...,Fm,...,FM}，Fm表示T’m对应的Frobenius范数图；
(2c)计算每个Frobenius范数图Fm的变差系数图，得到T′对应的变差系数图集合CV＝{CV1,CV2,...,CVm,...,CVM}，其中，CVm表示T’m对应的变差系数图；
(3)通过变差系数信息对初始边缘信息进行修正：
(3a)采用模糊C均值聚类算法，根据变差系数图CVm中每个像素点的变差系数信息对CVm进行三聚类，并按照每个像素簇内变差系数信息均值由小到大的顺序对CVm的三个聚类类别进行排序，得到CV对应的三聚类集合
Figure BDA0003213587370000031
其中，
Figure BDA0003213587370000032
分别表示CVm的三聚类结果中的高、中、低同质度类别；
(3b)根据变差系数图CVm的三聚类结果中的高同质度类别
Figure BDA0003213587370000033
与中同质度类别
Figure BDA0003213587370000034
对应区域之间的分界值
Figure BDA0003213587370000035
以及中同质度类别
Figure BDA0003213587370000036
和低同质度类别
Figure BDA0003213587370000037
对应区域之间的分界值
Figure BDA0003213587370000038
计算中同质度类别
Figure BDA0003213587370000039
对应区域内和低同质度类别
Figure BDA00032135873700000310
对应区域内的变差系数信息动态值
Figure BDA00032135873700000311
Figure BDA00032135873700000312
其中
Figure BDA00032135873700000313
是中同质度类别
Figure BDA00032135873700000314
对应区域内和低同质度类别
Figure BDA00032135873700000315
对应区域内Em(x,y)为1的像素点个数占对应区域内像素点总数的比例；
(3c)通过变差系数图CVm中所有像素点的变差系数信息的最大值
Figure BDA00032135873700000316
中同质度类别
Figure BDA00032135873700000317
和低同质度类别
Figure BDA00032135873700000318
对应区域内的变差系数动态值
Figure BDA00032135873700000319
对初始边缘信息图Em中所有位置为(x,y)的像素的边缘信息Em(x,y)进行修正，得到修正后的边缘信息图集合E′＝{E′1,E′2,...,E′m,...,E′M}，其中Em(x,y)的修正公式为：
Figure BDA00032135873700000320
其中，CVm(x,y)表示变差系数图CVm所在位置为(x,y)的像素的变差系数信息，E′m(x,y)表示Em(x,y)的修正结果；
(4)通过修正后的边缘信息对变差系数信息进行指导：
(4a)通过变差系数图CVm中高同质度类别
Figure BDA0003213587370000041
对应区域内所有像素点的变差系数信息的平均值
Figure BDA0003213587370000042
以及修正后的边缘信息图E′m中所有像素点的边缘信息的最大值
Figure BDA0003213587370000043
对每个时相的变差系数图CVm中所有位置为(x,y)的变差系数CVm(x,y)进行指导，得到指导后的初始综合邻域权重信息图集合CV′＝{CV′1,CV′2,...,CV′m,...,CV′M}，其中CVm(x,y)的修正公式为：
Figure BDA0003213587370000044
Figure BDA0003213587370000045
其中，E′m(x,y)表示修正后的边缘信息图E′m所在位置为(x,y)的像素的边缘信息，CV′m(x,y)表示初始综合邻域权重信息图CV′m所在位置为(x,y)的初始综合邻域权重信息；
(4b)对每个时相的初始综合邻域权重信息图CV′m进行归一化，得到范围在
Figure BDA0003213587370000046
修正后的边缘信息图E′m指导的综合邻域权重信息图CV″m，则CV′对应的综合邻域权重信息图集合为CV″＝{CV″1,CV″2,...,CV″m,...,CV″M}；
(5)设计平滑因子自适应于综合邻域权重信息的非局部均值滤波器：
(5a)对极化SAR Frobenius范数图集合F中每个时相的Fm求解灰度共生矩阵，得到F对应的灰度共生矩阵图集合G；
(5b)求解G中每个时相的灰度共生矩阵图Gm的同质度特征图Homm，并统计Homm中在高同质度类别
Figure BDA0003213587370000047
对应区域
Figure BDA0003213587370000048
中所有像素点的平均值
Figure BDA0003213587370000049
(5c)通过
Figure BDA00032135873700000410
综合邻域权重信息图CV″m中所有像素点对应的综合邻域权重信息的最大值
Figure BDA0003213587370000051
对每个时相的初始平滑因子图hm中所有位置为(x,y)的初始平滑因子值hm(x,y)进行更新，得到更新后的自适应平滑因子图h′m，并设计h′m自适应于综合邻域权重信息图CV″m的非局部均值滤波器Filterm，以及Filterm对应的平滑因子自适应于综合邻域权重信息的非局部均值滤波器集合Filter＝{Filter1,Filter2,...,Filterm,...,FilterM}，其中hm(x,y)的更新公式为：
h′m(x,y)＝hm(x,y)(1-δ(x,y)×ε(δ(x,y)))
Figure BDA0003213587370000052
Figure BDA0003213587370000053
其中，h′m(x,y)表示更新后的自适应平滑因子图h′m所在位置为(x,y)的像素的自适应平滑因子值，log(·)表示求以自然常数e为底的对数，CV″m(x,y)表示综合邻域权重信息图CV″m所在位置为(x,y)的综合邻域权重信息；
(6)获取综合邻域信息差异图：
通过极化SAR Frobenius范数图集合F＝{F1,F2,...,Fm,...,FM}中任意p、q时相的极化SAR Frobenius范数图Fp、Fq，得到任意p、q时相间的综合邻域信息差异图DIp,q，其中DIp,q中位置为(x,y)的像素点的差异度值计算公式为：
Figure BDA0003213587370000054
A＝min(Fp(x,y),Fq(x,y))
B＝max(Fp(x,y),Fq(x,y))
Figure BDA0003213587370000055
其中，A′、B′表示A、B当前所在的极化SARFrobenius范数图的时相值，FilterA′[·]、FilterB′[·]表示利用第A′、B′个时相的平滑因子自适应于综合邻域权重信息的滤波器FilterA′、FilterB′对A、B进行滤波操作，Fp(x,y)、Fq(x,y)表示第p、q个时相的极化SARFrobenius范数图Fp、Fq所在位置为(x,y)的Frobenius范数值，CV″p(x,y)、CV″q(x,y)表示第p、q个时相综合邻域权重信息图CV″p、CV″q所在位置为(x,y)的综合邻域权重信息，CV″p,q(x,y)是CV″p(x,y)与CV″q(x,y)的平均值；
(7)获取基于综合邻域信息差异图的显著性二值检测图：
(7a)利用超像素分割方法，对DIp,q进行超像素分割，得到超像素分割区域集合
Figure BDA0003213587370000061
并采用结合超像素分割的显著性检测方法，通过Areap,q对DIp,q进行显著性检测，得到DIp,q的显著性检测图
Figure BDA0003213587370000062
其中，U≥2，
Figure BDA0003213587370000063
表示第u个超像素分割子区域；
(7b)采用模糊C均值聚类算法，对DIp,q与
Figure BDA0003213587370000064
进行二聚类，并按照聚类结果中每个类别均值由大到小的顺序对两个聚类类别进行排序，得到DIp,q对应的聚类集合
Figure BDA0003213587370000065
以及
Figure BDA0003213587370000066
对应的聚类集合
Figure BDA0003213587370000067
其中
Figure BDA0003213587370000068
表示初始变化类，对应于DIp,q中的初始变化区域
Figure BDA0003213587370000069
表示的初始未变化类，对应于DIp,q中的初始未变化区域
Figure BDA00032135873700000610
表示非主体显著类，对应于
Figure BDA00032135873700000611
中的非主体显著区域
Figure BDA00032135873700000612
Figure BDA00032135873700000613
表示主体显著类，对应于
Figure BDA00032135873700000614
中的主体显著区域
Figure BDA00032135873700000615
将
Figure BDA00032135873700000616
中的所有像素点的像素值置1，同时将
Figure BDA00032135873700000617
中的所有像素点的像素值置0；
(7c)将初始变化区域
Figure BDA00032135873700000618
区域内值为1的像素点且主体显著区域
Figure BDA00032135873700000619
内值不为1的像素点对应的区域设置为显著性遗失区域
Figure BDA00032135873700000620
统计每个超像素分割子区域
Figure BDA00032135873700000621
中包含像素值的个数
Figure BDA00032135873700000622
并对显著性遗失区域
Figure BDA00032135873700000623
在每个对应
Figure BDA00032135873700000624
内值为1的像素点个数进行统计，得到统计个数Numu，将Numu大于
Figure BDA00032135873700000625
的对应
Figure BDA00032135873700000626
内的所有像素值设定为1，得到显著性补充区域
Figure BDA00032135873700000627
(7d)对主体显著区域
Figure BDA00032135873700000628
中像素值为1的封闭区域进行统计，得到封闭区域集合L＝{L1,L2,...,Ln,...,LN}，并统计每一个封闭子区域Ln中包含的像素点个数，按照个数进行降序排序，得到降序排序后的封闭子区域Ln中的像素数量集合Num＝{Num1,Num2,...,Numn,...,NumN}，计算相邻两个封闭子区域Ln-1、Ln内像素数量Numn-1、Numn间的比值集合Rio＝{Rio1,Rio2,...,Rion,...,RioN-1}，其中，N≥2，Ln表示第n个封闭区域，Numn表示第n个封闭区域的像素数量，Rion表示Numn-1与Numn的比值；
(7e)统计上述比值集合Rio中的最大值Riomax对应的下标nmax，以下标nmax的前一个下标nmax-1为分隔下标，将L中下标大于nmax-1的封闭子区域划分至主体显著性区域
Figure BDA0003213587370000071
中，将下标小于或等于nmax-1的封闭子区域划分至孤立区域
Figure BDA0003213587370000072
中；
(7f)将孤立区域
Figure BDA0003213587370000073
与显著性补充区域
Figure BDA0003213587370000074
进行求与操作，得到完整显著性补充区域
Figure BDA0003213587370000075
对显著性检测图
Figure BDA0003213587370000076
进行自适应阈值化操作，得到
Figure BDA0003213587370000077
对应的完整主体显著性区域
Figure BDA0003213587370000078
将
Figure BDA0003213587370000079
与
Figure BDA00032135873700000710
进行求与操作，得到综合邻域信息差异图DIp,q的显著性二值检测图
Figure BDA00032135873700000711
(8)获取显著性差异图：
将综合邻域信息差异图DIp,q与显著性二值检测图
Figure BDA00032135873700000712
同一位置像素的值进行相乘，实现对DIp,q的修正，得到显著性差异图
Figure BDA00032135873700000713
其中
Figure BDA00032135873700000714
表示DIp,q和
Figure BDA00032135873700000715
同一位置像素的值相乘；
(9)获取极化SAR图像的变化检测结果：
利用阈值二值化方法对
Figure BDA00032135873700000716
进行二值化分割，得到极化SAR变化检测结果CDp,q。
本发明与现有的技术相比，具有以下优点：
第一，本发明通过任意两个时相极化SAR图像的边缘信息指导的综合邻域权重信息及平滑因子自适应于该综合邻域权重信息的非局部均值滤波器，构建了任意两个时相间的综合邻域信息差异图，避免了现有技术获取的差异图因过度均值滤波而无法保留完整边缘信息的缺陷以及因未对异质区域包含的复杂信息进行有效区分而保留了局部微小斑点噪声的缺陷，提高了对边缘细节信息的保留能力与局部微小噪声的抑制能力，有效地提升了极化SAR变化检测的总体精度与区域内部一致性。
第二，本发明通过面向综合邻域信息差异图的自适应显著性二值检测图来对综合邻域信息差异图进行修正，构建了对应的显著性差异图来突出变化区域并抑制非变化区域，避免了现有技术中产生的差异图因未对异质区域包含的复杂信息进行有效区分而保留了局部微小斑点噪声的缺陷，提高了对微小噪声的抑制能力，有效地提升了极化SAR变化检测的准确率与区域内部一致性。实验结果表明，本发明能够在极化SAR变化检测中获得更高的总体精度与区域内部一致性。
附图说明
图1为本发明的实现流程图；
图2为本发明和现有技术应用于两时相太湖水域1数据的极化SAR图像的变化检测结果对比图；
图3为本发明和现有技术应用于两时相太湖水域2数据的极化SAR图像的变化检测结果对比图。
具体实施方式
下面结合附图和具体实施例，对本发明作进一步详细描述。
参照图1，本发明包括如下步骤：
步骤1对同一地区M个时相的极化SAR图像进行配准：
以同一地区M个时相的极化SAR图像T＝{T1,T2,...,Tm,...,TM}中任意一个时相的极化SAR图像作为参考图像，并将其它时相的极化SAR图像与参考图像进行配准，得到配准后的极化SAR图像集合T′＝{T′1,T′2,...,T′m,...,T′M}，其中，M≥2，Tm表示像素大小为w×d的第m个时相的极化SAR图像，w≥200，d≥200，T’m表示Tm的配准结果；以两时相太湖水域1数据为例，本实施例中，M＝2，w＝272，d＝240；
步骤2获取配准后的每个时相的极化SAR图像对应的初始边缘信息和变差系数信息：
(2a)对配准后的每个时相的极化SAR图像T’m进行极化分解，得到T′对应的极化SAR图像分解合成图集合P＝{P1,P2,...,Pm,...,PM}，本实施例中选择了应用范围广、效果良好的Pauli极化分解方法进行极化分解，具体实现步骤为：
将T’m中的每一个位置(x,y)处的极化SAR数据T’m(x,y)转换为极化SAR散射矩阵Sm(x,y)，并对极化SAR数据进行Pauli分解：
Figure BDA0003213587370000091
其中，[Sm(x,y)]表示极化SAR散射矩阵，
Figure BDA0003213587370000092
表示接收方式为i′、发射方式为j′的散射矩阵，H表示水平方式，V表示垂直方式，
Figure BDA0003213587370000093
表示Pauli基，a′，b′，c′和d′都是复数，当介质满足互易条件时，可将这些量写为向量Km的形式：
Figure BDA0003213587370000094
Pauli基是完备正交基，具有一定的抗噪性，|a′m|2，|b′m|2和|c′m|2，对应着明显的物理散射机制，所以利用Pauli分解各系数，可以合成任意时相RGB图像Pm：
Figure BDA0003213587370000095
其中，Red、Green和Blue分别代表合成图像中的红、绿蓝三通道信息；
在得到极化SAR图像分解合成图集合P后，进一步对每个时相的极化SAR图像的分解合成图Pm进行边缘检测，得到T′对应的初始边缘信息图集合E＝{E1,E2,...,Em,...,EM}，其中，Pm、Em表示Tm对应的分解合成图、初始边缘信息图；本实施例中，边缘检测方法采用Canny边缘检测方法，具体实现步骤为：
将每个时相的Pm进行灰度化处理，得到Pm对应的灰度图Graym，进一步先后对Graym进行高斯滤波器滤波、计算各像素梯度强度与方向、非极大值抑制杂散响应、双阈值检测以及抑制孤立弱边缘，得到Canny边缘检测结果Em；
(2b)计算配准后的每个时相的极化SAR图像T’m的Frobenius范数，得到T′对应的Frobenius范数图集合F＝{F1,F2,...,Fm,...,FM}，Fm表示T’m对应的Frobenius范数图，其中T’m中位置为(x,y)的像素点的Frobenius范数Fm(x,y)的计算公式为：
Figure BDA0003213587370000101
其中，||·||F表示求Frobenius范数操作，|·|表示求绝对值操作，t′ij(x,y)表示T’m(x,y)中第i行第j列的元素值；
(2c)计算每个Frobenius范数图Fm的变差系数图，得到T′对应的变差系数图集合CV＝{CV1,CV2,...,CVm,...,CVM}，CVm表示T’m对应的变差系数图，每个变差系数图CVm中位置为(x,y)的像素点的变差系数CVm(x,y)的实现步骤为：
以Fm中每个位置为(x,y)的像素点的Frobenius范数Fm(x,y)为中心像素，将Fm划分成设定大小为sizeF×sizeF的邻域窗口
Figure BDA0003213587370000102
本实施例中，sizeF＝7，并计算Fm(x,y)在对应邻域窗
Figure BDA0003213587370000103
内的变差系数CVm(x,y)：
Figure BDA0003213587370000104
其中，Stdm(x,y)是Fm(x,y)在其对应邻域窗
Figure BDA0003213587370000105
内的标准差；μm(x,y)是Fm(x,y)在其对应邻域窗
Figure BDA0003213587370000106
内的均值；
步骤3通过变差系数信息对初始边缘信息进行修正：
(3a)采用模糊C均值聚类算法，根据变差系数图CVm中每个像素点的变差系数信息对CVm进行三聚类，并按照每个像素簇内变差系数信息均值由小到大的顺序对CVm的三个聚类类别进行排序，得到CV对应的三聚类集合
Figure BDA0003213587370000107
其中，
Figure BDA0003213587370000108
分别表示CVm的三聚类结果中的高、中、低同质度类别；
(3b)根据变差系数图CVm的三聚类结果中的高同质度类别
Figure BDA0003213587370000109
与中同质度类别
Figure BDA00032135873700001010
对应区域之间的分界值
Figure BDA0003213587370000111
以及中同质度类别
Figure BDA0003213587370000112
和低同质度类别
Figure BDA0003213587370000113
对应区域之间的分界值
Figure BDA0003213587370000114
计算中同质度类别
Figure BDA0003213587370000115
对应区域内和低同质度类别
Figure BDA0003213587370000116
对应区域内的变差系数信息动态值
Figure BDA0003213587370000117
Figure BDA0003213587370000118
其中
Figure BDA0003213587370000119
是中同质度类别
Figure BDA00032135873700001110
对应区域内和低同质度类别
Figure BDA00032135873700001111
对应区域内Em(x,y)为1的像素点个数占对应区域内像素点总数的比例；
(3c)通过变差系数图CVm中所有像素点的变差系数信息的最大值
Figure BDA00032135873700001112
中同质度类别
Figure BDA00032135873700001113
和低同质度类别
Figure BDA00032135873700001114
对应区域内的变差系数动态值
Figure BDA00032135873700001115
对初始边缘信息图Em中所有位置为(x,y)的像素的边缘信息Em(x,y)进行修正，得到修正后的边缘信息图集合E′＝{E′1,E′2,...,E′m,...,E′M}，其中Em(x,y)的修正公式为：
Figure BDA00032135873700001116
其中，CVm(x,y)表示变差系数图CVm所在位置为(x,y)的像素的变差系数信息，E′m(x,y)表示Em(x,y)的修正结果；
步骤4通过修正后的边缘信息对变差系数信息进行指导：
(4a)通过变差系数图CVm中高同质度类别
Figure BDA00032135873700001117
对应区域内所有像素点的变差系数信息的平均值
Figure BDA00032135873700001118
以及修正后的边缘信息图E′m中所有像素点的边缘信息的最大值
Figure BDA00032135873700001119
对每个时相的变差系数图CVm中所有位置为(x,y)的变差系数CVm(x,y)进行指导，得到指导后的初始综合邻域权重信息图集合CV′＝{CV′1,CV′2,...,CV′m,...,CV′M}，其中CVm(x,y)的修正公式为：
Figure BDA00032135873700001120
Figure BDA0003213587370000121
其中，E′m(x,y)表示修正后的边缘信息图E′m所在位置为(x,y)的像素的边缘信息，CV′m(x,y)表示初始综合邻域权重信息图CV′m所在位置为(x,y)的初始综合邻域权重信息；
(4b)对每个时相的初始综合邻域权重信息图CV′m进行归一化，得到范围在
Figure BDA0003213587370000122
修正后的边缘信息图E′m指导的综合邻域权重信息图CV″m，则CV′对应的综合邻域权重信息图集合为CV″＝{CV″1,CV″2,...,CV″m,...,CV″M}；
步骤5设计平滑因子自适应于综合邻域权重信息的非局部均值滤波器：
(5a)对极化SAR Frobenius范数图集合F中每个时相的极化SAR Frobenius范数图Fm求解灰度共生矩阵，得到F对应的灰度共生矩阵图集合G，G中每个时相灰度共生矩阵Gm的求解过程为：
在Fm中任意一点(x,y)及偏离它的一点(x+α,y+β)(其中α、β为整数)构成点对，设该点对的灰度值为(gr1,gr2),假设图像的最大灰度级为gk,则gr1与gr2的组合共有gk×gk种。对于整幅Fm，统计每一种(gr1,gr2)值出现的次数，排列成一个方阵，再用(gr1,gr2)出现的总次数将它们归一化为出现的概率，得到灰度共生矩阵Gm；
(5b)求解G中每个时相的灰度共生矩阵图Gm的同质度特征图Homm，其中Homm中位置为(x,y)的像素点的同质度特征值Homm(x,y)的计算公式为：
Figure BDA0003213587370000123
其中，gk是每个时相极化SAR Frobenius范数图Fm图的灰度值级数，
Figure BDA0003213587370000124
是灰度共生矩阵图Gm所在位置为(x,y)的像素点的灰度共生矩阵值Gm(x,y)中第l行第s列的元素值；本实例中，gk＝16；
在得到Homm之后，进一步统计Homm中在高同质度类别
Figure BDA0003213587370000125
对应区域
Figure BDA0003213587370000126
中所有像素点的平均值
Figure BDA0003213587370000131
(5c)通过
Figure BDA0003213587370000132
综合邻域权重信息图CV″m中所有像素点对应的综合邻域权重信息的最大值
Figure BDA0003213587370000133
对每个时相的初始平滑因子图hm中所有位置为(x,y)的初始平滑因子值hm(x,y)进行更新，得到更新后的自适应平滑因子图h′m，并设计h′m自适应于综合邻域权重信息图CV″m的非局部均值滤波器Filterm，以及Filterm对应的平滑因子自适应于综合邻域权重信息的非局部均值滤波器集合Filter＝{Filter1,Filter2,...,Filterm,...,FilterM}，其中hm(x,y)的更新公式为：
h′m(x,y)＝hm(x,y)(1-δ(x,y)×ε(δ(x,y)))
Figure BDA0003213587370000134
Figure BDA0003213587370000135
其中，h′m(x,y)表示更新后的自适应平滑因子图h′m所在位置为(x,y)的像素的自适应平滑因子值，log(·)表示求以自然常数e为底的对数，CV″m(x,y)表示综合邻域权重信息图CV″m所在位置为(x,y)的综合邻域权重信息；本实施例中，时相1的初始平滑因子图h1中所有位置为(x,y)的初始平滑因子值均设定为11.3，时相2的初始平滑因子图h2中所有位置为(x,y)的初始平滑因子值均设定为10.1；
上述步骤中，非局部均值滤波器的具体设计方法为：
Figure BDA0003213587370000136
Figure BDA0003213587370000137
其中，B′(ac,rr)表示一个以ac为中心，大小为(2rr+1)×(2rr+1)像素的邻域，Norm(ac)表示归一化因子，u′是滤波前的值，
Figure BDA0003213587370000138
是滤波后的值，本实施例中，rr＝3；
权重w′(ac,a′c)取决于分别以ac和a′c为中心的(2rr+1)×(2rr+1)大小的彩色patch块的欧氏距离dis2(B′(ac,rr),B′(a′c,rr))的平方：
Figure BDA0003213587370000141
利用指数核函数并结合自适应平滑因子h′来计算w′(ac,a′c)：
Figure BDA0003213587370000142
式中，σ是噪声的标准差，h′为平滑参数，e为自然常数，控制着高斯函数的衰减程度。h′越大高斯函数变化越平缓，去噪水平越高，但同时这也会导致图像越模糊；反之，h′越小，边缘细节成分保持得越多，但会残留过多的噪声点；
步骤6获取综合邻域信息差异图：
通过极化SAR Frobenius范数图集合F＝{F1,F2,...,Fm,...,FM}中任意p、q时相的极化SAR Frobenius范数图Fp、Fq，得到任意p、q时相间的综合邻域信息差异图DIp,q，其中DIp,q中位置为(x,y)的像素点的差异度值计算公式为：
Figure BDA0003213587370000143
A＝min(Fp(x,y),Fq(x,y))
B＝max(Fp(x,y),Fq(x,y))
Figure BDA0003213587370000144
其中，A′、B′表示A、B当前所在的极化SARFrobenius范数图的时相值，FilterA′[·]、FilterB′[·]表示利用第A′、B′个时相的平滑因子自适应于综合邻域权重信息的滤波器FilterA′、FilterB′对A、B进行滤波操作，Fp(x,y)、Fq(x,y)表示第p、q个时相的极化SARFrobenius范数图Fp、Fq所在位置为(x,y)的Frobenius范数值，CV″p(x,y)、CV″q(x,y)表示第p、q个时相综合邻域权重信息图CV″p、CV″q所在位置为(x,y)的综合邻域权重信息，CV″p,q(x,y)是CV″p(x,y)与CV″q(x,y)的平均值；本实施例中，p＝1，q＝2；
步骤7获取基于综合邻域信息差异图的显著性二值检测图：
(7a)利用超像素分割方法，对DIp,q进行超像素分割，得到超像素分割区域集合
Figure BDA0003213587370000151
表示第u个超像素分割子区域，本实例中超像素分割方法采用效果优良且应用广泛的简单线性迭代聚类(SLIC)方法，其中初始分割子区域的窗口尺寸为rslic＝27，并采用结合超像素分割的显著性检测方法，通过
Figure BDA0003213587370000152
对DIp,q进行显著性检测，实现过程为：
统计Areap,q中每个超像素分割子区域
Figure BDA0003213587370000153
内所有像素值的众数
Figure BDA0003213587370000154
以及像素值的总个数
Figure BDA0003213587370000155
并遍历
Figure BDA0003213587370000156
中的每一个像素，得到每个
Figure BDA0003213587370000157
中第一个像素值等于
Figure BDA0003213587370000158
的像素点
Figure BDA0003213587370000159
通过显著性检测方法对每一个
Figure BDA00032135873700001510
进行显著性检测，得到每个
Figure BDA00032135873700001511
对应的显著性检测值
Figure BDA00032135873700001512
将对应
Figure BDA00032135873700001513
中的其余
Figure BDA00032135873700001514
个像素点的显著检测值设定为
Figure BDA00032135873700001515
得到DIp,q对应的显著性检测结果
Figure BDA00032135873700001516
本实施例中，显著性检测方法选择了应用广泛、效果良好的上下文感知显著性检测(CA)方法，具体设计过程为：
首先获得一个衡量显著性的综合距离，上下文感知显著性检测方法对图像块之间在Lab颜色空间做对比，若某一块与其它块差距大，则说明是显著性特征。考虑空间距离，与显著性区域相似的区域一般都离得比较近，据此可以得到综合距离：
Figure BDA00032135873700001517
其中，Patchci、Patchcj分别表示中心像素为ci、cj的像素块，dcolor(·)代表两点间的Lab颜色距离，dposition(·)代表两点间的空间距离，一对图像块的综合距离越大，差异度就越大，若对于任意的块，得到的差异值都很大，则可判定该块为显著块；本实施例中，para＝3；
进一步确定显著性公式，为降低算法的计算复杂度，只计算前K′个与最相似的块即可，得到的显著性公式为：
Figure BDA0003213587370000161
其中，ci,k′分别像素块
Figure BDA0003213587370000162
的中心像素，r′表示尺度值，
Figure BDA0003213587370000163
表示像素ci处尺度为r′的显著性值，本实施例中，K′＝65；
将上式的单个尺度所得到的显著性值引申为求多个尺度下的显著性值的平均值来进一步提高显著与非显著区域的对比度，即：
Figure BDA0003213587370000164
式中，
Figure BDA0003213587370000165
表示像素ci处的多尺度平均显著值，M′为尺度的个数，R′＝{r′1,...,r′M}；本实施例中，M′＝4；
进一步加入上下文修正，设定一个显著性阈值TCA，并从显著图中提取注意力最大的局部区域，在注意力区域之外的像素显著性值由与其最近的注意力像素之间的欧式距离加权，得到新的显著性值：
Figure BDA0003213587370000166
式中，dfoci(ci)为离像素点ci最近的注意力像素之间的欧式距离，
Figure BDA0003213587370000167
表示像素ci处的多尺度平均显著值，
Figure BDA0003213587370000168
表示像素ci处的最终显著性值；本实施例中，TCA＝0.8；
(7b)采用模糊C均值聚类算法，对DIp,q与
Figure BDA0003213587370000169
进行二聚类，并按照聚类结果中每个类别均值由大到小的顺序对两个聚类类别进行排序，得到DIp,q对应的聚类集合
Figure BDA00032135873700001610
以及
Figure BDA00032135873700001611
对应的聚类集合
Figure BDA00032135873700001612
其中
Figure BDA00032135873700001613
表示初始变化类，对应于DIp,q中的初始变化区域
Figure BDA00032135873700001614
表示的初始未变化类，对应于DIp,q中的初始未变化区域
Figure BDA00032135873700001615
表示非主体显著类，对应于
Figure BDA00032135873700001616
中的非主体显著区域
Figure BDA00032135873700001617
Figure BDA00032135873700001618
表示主体显著类，对应于
Figure BDA00032135873700001619
中的主体显著区域
Figure BDA00032135873700001620
将
Figure BDA00032135873700001621
中的所有像素点的像素值置1，同时将
Figure BDA00032135873700001622
中的所有像素点的像素值置0；
(7c)将初始变化区域
Figure BDA00032135873700001623
区域内值为1的像素点且主体显著区域
Figure BDA00032135873700001624
内值不为1的像素点对应的区域设置为显著性遗失区域
Figure BDA0003213587370000171
统计每个超像素分割子区域
Figure BDA0003213587370000172
中包含像素值的个数
Figure BDA0003213587370000173
并对显著性遗失区域
Figure BDA0003213587370000174
在每个对应
Figure BDA0003213587370000175
内值为1的像素点个数进行统计，得到统计个数Numu，将Numu大于
Figure BDA0003213587370000176
的对应
Figure BDA0003213587370000177
内的所有像素值设定为1，得到显著性补充区域
Figure BDA0003213587370000178
(7d)对主体显著区域
Figure BDA0003213587370000179
中像素值为1的封闭区域进行统计，得到封闭区域集合L＝{L1,L2,...,Ln,...,LN}，并统计每一个封闭子区域Ln中包含的像素点个数，按照个数进行降序排序，得到降序排序后的封闭子区域Ln中的像素数量集合Num＝{Num1,Num2,...,Numn,...,NumN}，计算相邻两个封闭子区域Ln-1、Ln内像素数量Numn-1、Numn间的比值集合Rio＝{Rio1,Rio2,...,Rion,...,RioN-1}，其中，N≥2，Ln表示第n个封闭区域，Numn表示第n个封闭区域的像素数量，Rion表示Numn-1与Numn的比值；
(7e)统计上述比值集合Rio中的最大值Riomax对应的下标nmax，以下标nmax的前一个下标nmax-1为分隔下标，将L中下标大于nmax-1的封闭子区域划分至主体显著性区域
Figure BDA00032135873700001710
中，将下标小于或等于nmax-1的封闭子区域划分至孤立区域
Figure BDA00032135873700001711
中；
(7f)将孤立区域
Figure BDA00032135873700001712
与显著性补充区域
Figure BDA00032135873700001713
进行求与操作，得到完整显著性补充区域
Figure BDA00032135873700001714
对显著性检测图
Figure BDA00032135873700001715
进行自适应阈值化操作，得到
Figure BDA00032135873700001716
对应的完整主体显著性区域
Figure BDA00032135873700001717
实现步骤为：
(7f1)设定阈值TCA的初始值为0，以Δ为步长增加TCA；本实施例中，Δ＝0.01；
(7f2)将显著性检测图
Figure BDA00032135873700001718
中显著性值大于TCA的像素点的像素值设定为1，小于TCA的像素点的像素值设定为0，得到
Figure BDA00032135873700001719
对应的显著性阈值图
Figure BDA00032135873700001720
其中，0≤TCA≤1；
(7f3)将主体显著性区域
Figure BDA00032135873700001721
中的像素点与
Figure BDA00032135873700001722
中对应像素点进行求与操作，并统计求与结果中像素值为1的像素点个数numone，同时统计
Figure BDA00032135873700001723
中的像素点个数num′；
(7f4)判断numone＜num′是否成立，若是，则以Δ为步长增加TCA，并执行步骤(7f2)，否则，得到TCA对应的中止阈值Tend，将显著性检测图
Figure BDA0003213587370000181
中显著性值大于Tend的像素点的像素值设定为1，小于TCA的像素点的像素值设定为0，得到
Figure BDA0003213587370000182
对应的完整主体显著性区域
Figure BDA0003213587370000183
之后，将
Figure BDA0003213587370000184
与
Figure BDA0003213587370000185
进行求与操作，得到综合邻域信息差异图DIp,q的显著性二值检测图
Figure BDA0003213587370000186
步骤8获取显著性差异图：
将综合邻域信息差异图DIp,q与显著性二值检测图
Figure BDA0003213587370000187
同一位置像素的值进行相乘，实现对DIp,q的修正，得到显著性差异图
Figure BDA0003213587370000188
其中
Figure BDA0003213587370000189
表示DIp,q和
Figure BDA00032135873700001810
同一位置像素的值相乘；
步骤9获取极化SAR图像的变化检测结果：
利用阈值二值化方法对
Figure BDA00032135873700001811
进行二值化分割，得到极化SAR变化检测结果CDp,q；本实施例中，阈值二值化方法采用大津法。
本发明效果可以通过以下实验进一步证实：
1.仿真条件和内容：
仿真实验环境为：MATLAB R2016a，Intel(R)Core(TM)i7-8700 CPU 3.20GHz，Window 10。
采用两时相不同数据对本发明与现有的基于极化散射总功率比值与邻域变差系数权重因子的极化SAR变化检测方法的变化检测总体精度和内部一致性进行验证，其结果如仿真结果分析中的表1所示，其中内部一致性是通过Kappa系数进行评价的。
仿真1，对本发明和现有技术应用于两时相太湖水域1数据的极化SAR图像的变化检测结果进行对比仿真，其结果如图2所示。
仿真2，对本发明和现有技术应用于两时相太湖水域2数据的极化SAR图像的变化检测结果进行对比仿真，其结果如图3所示。
2.仿真结果分析：
参照图2，其中图2(a)是太湖水域1数据现有方法的检测结果，图2(b)是太湖水域1数据本发明方法的变化检测结果，图2(c)是太湖水域1数据变化检测结果参考图，图中白色区域表示变化区域，黑色区域表示未变化区域，分别将图2(a)、图2(b)与图2(c)进行对照，可以看出图2(b)在对局部边缘信息的保留能力以及局部微小噪声的抑制能力相较于图2(c)更强，与参考图给出结果更加接近；
参照图3，其中图3(a)是太湖水域2数据现有方法的检测结果，图3(b)是太湖水域2数据本发明方法的变化检测结果，图3(c)是太湖水域2数据变化检测结果参考图，图中白色区域表示变化区域，黑色区域表示未变化区域，分别将图3(a)、图3(b)与图3(c)进行对照，可以看出图3(b)在对局部边缘信息的保留能力以及局部微小噪声的抑制能力相较于图3(c)更强，与参考图给出结果更加接近；
表1
Figure BDA0003213587370000191
参照表1，本发明在检测精度上优于现有技术。这是由于本发明通过任意两个时相极化SAR图像的边缘信息指导的综合邻域权重信息及平滑因子自适应于该综合邻域权重信息的非局部均值滤波器，构建了任意两个时相间的综合邻域信息差异图，避免了现有技术获取的差异图因过度均值滤波而无法保留完整边缘信息的缺陷以及因未对异质区域包含的复杂信息进行有效区分而保留了局部微小斑点噪声的缺陷，此外，本发明通过面向综合邻域信息差异图的自适应显著性二值检测图来对综合邻域信息差异图进行修正，构建了对应的显著性差异图来突出变化区域并抑制非变化区域，避免了现有技术产生的差异图因未对异质区域包含的复杂信息进行有效区分而保留了局部微小斑点噪声的缺陷，因此本发明的变化检测总体精度与区域内部一致性明显高于现有技术。

Claims (6)

1.一种基于综合邻域信息的显著性极化SAR图像变化检测方法，其特征在于，包括如下步骤：
(1)对同一地区M个时相的极化SAR图像进行配准：
以同一地区M个时相的极化SAR图像T＝{T1,T2,...,Tm,...,TM}中任意一个时相的极化SAR图像作为参考图像，并将其它时相的极化SAR图像与参考图像进行配准，得到配准后的极化SAR图像集合T′＝{T1′,T′2,...,T′m,...,T′M}，其中，M≥2，Tm表示像素大小为w×d的第m个时相的极化SAR图像，w≥200，d≥200，Tm’表示Tm的配准结果；
(2)获取配准后的每个时相的极化SAR图像对应的初始边缘信息和变差系数信息：
(2a)对配准后的每个时相的极化SAR图像Tm’进行极化分解，得到T′对应的极化SAR图像分解合成图集合P＝{P1,P2,...,Pm,...,PM}，并对每个时相的极化SAR图像的分解合成图Pm进行边缘检测，得到T′对应的初始边缘信息图集合E＝{E1,E2,...,Em,...,EM}，其中，Pm、Em表示Tm对应的分解合成图、初始边缘信息图；
(2b)计算配准后的每个时相的极化SAR图像Tm’的Frobenius范数，得到T′对应的Frobenius范数图集合F＝{F1,F2,...,Fm,...,FM}，Fm表示Tm’对应的Frobenius范数图；
(2c)计算每个Frobenius范数图Fm的变差系数图，得到T′对应的变差系数图集合CV＝{CV1,CV2,...,CVm,...,CVM}，其中，CVm表示Tm’对应的变差系数图；
(3)通过变差系数信息对初始边缘信息进行修正：
(3a)采用模糊C均值聚类算法，根据变差系数图CVm中每个像素点的变差系数信息对CVm进行三聚类，并按照每个像素簇内变差系数信息均值由小到大的顺序对CVm的三个聚类类别进行排序，得到CV对应的三聚类集合
Figure FDA0003213587360000021
其中，
Figure FDA0003213587360000022
分别表示CVm的三聚类结果中的高、中、低同质度类别；
(3b)根据变差系数图CVm的三聚类结果中的高同质度类别
Figure FDA0003213587360000023
与中同质度类别
Figure FDA0003213587360000024
对应区域之间的分界值
Figure FDA0003213587360000025
以及中同质度类别
Figure FDA0003213587360000026
和低同质度类别
Figure FDA0003213587360000027
对应区域之间的分界值
Figure FDA0003213587360000028
计算中同质度类别
Figure FDA0003213587360000029
对应区域内和低同质度类别
Figure FDA00032135873600000210
对应区域内的变差系数信息动态值
Figure FDA00032135873600000211
Figure FDA00032135873600000212
其中
Figure FDA00032135873600000213
是中同质度类别
Figure FDA00032135873600000214
对应区域内和低同质度类别
Figure FDA00032135873600000215
对应区域内Em(x,y)为1的像素点个数占对应区域内像素点总数的比例；
(3c)通过变差系数图CVm中所有像素点的变差系数信息的最大值
Figure FDA00032135873600000216
中同质度类别
Figure FDA00032135873600000217
和低同质度类别
Figure FDA00032135873600000218
对应区域内的变差系数动态值
Figure FDA00032135873600000219
对初始边缘信息图Em中所有位置为(x,y)的像素的边缘信息Em(x,y)进行修正，得到修正后的边缘信息图集合E′＝{E′1,E′2,...,E′m,...,E′M}，其中Em(x,y)的修正公式为：
Figure FDA00032135873600000220
其中，CVm(x,y)表示变差系数图CVm所在位置为(x,y)的像素的变差系数信息，E′m(x,y)表示Em(x,y)的修正结果；
(4)通过修正后的边缘信息对变差系数信息进行指导：
(4a)通过变差系数图CVm中高同质度类别
Figure FDA00032135873600000221
对应区域内所有像素点的变差系数信息的平均值
Figure FDA0003213587360000031
以及修正后的边缘信息图E′m中所有像素点的边缘信息的最大值
Figure FDA0003213587360000032
对每个时相的变差系数图CVm中所有位置为(x,y)的变差系数CVm(x,y)进行指导，得到指导后的初始综合邻域权重信息图集合CV′＝{CV1′,CV′2,...,CV′m,...,CV′M}，其中CVm(x,y)的修正公式为：
Figure FDA0003213587360000033
Figure FDA0003213587360000034
其中，E′m(x,y)表示修正后的边缘信息图E′m所在位置为(x,y)的像素的边缘信息，CV′m(x,y)表示初始综合邻域权重信息图CV′m所在位置为(x,y)的初始综合邻域权重信息；
(4b)对每个时相的初始综合邻域权重信息图CV′m进行归一化，得到范围在
Figure FDA0003213587360000035
修正后的边缘信息图E′m指导的综合邻域权重信息图CV″m，则CV′对应的综合邻域权重信息图集合为CV″＝{CV1″,CV″2,...,CV″m,...,CV″M}；
(5)设计平滑因子自适应于综合邻域权重信息的非局部均值滤波器：
(5a)对极化SAR Frobenius范数图集合F中每个时相的Fm求解灰度共生矩阵，得到F对应的灰度共生矩阵图集合G；
(5b)求解G中每个时相的灰度共生矩阵图Gm的同质度特征图Homm，并统计Homm中在高同质度类别
Figure FDA0003213587360000036
对应区域
Figure FDA0003213587360000037
中所有像素点的平均值
Figure FDA0003213587360000038
(5c)通过
Figure FDA0003213587360000039
综合邻域权重信息图CV″m中所有像素点对应的综合邻域权重信息的最大值
Figure FDA00032135873600000310
对每个时相的初始平滑因子图hm中所有位置为(x,y)的初始平滑因子值hm(x,y)进行更新，得到更新后的自适应平滑因子图h′m，并设计h′m自适应于综合邻域权重信息图CV″m的非局部均值滤波器Filterm，以及Filterm对应的平滑因子自适应于综合邻域权重信息的非局部均值滤波器集合Filter＝{Filter1,Filter2,...,Filterm,...,FilterM}，其中hm(x,y)的更新公式为：
h′m(x,y)＝hm(x,y)(1-δ(x,y)×ε(δ(x,y)))
Figure FDA0003213587360000041
Figure FDA0003213587360000042
其中，h′m(x,y)表示更新后的自适应平滑因子图h′m所在位置为(x,y)的像素的自适应平滑因子值，log(·)表示求以自然常数e为底的对数，CV″m(x,y)表示综合邻域权重信息图CV″m所在位置为(x,y)的综合邻域权重信息；
(6)获取综合邻域信息差异图：
通过极化SAR Frobenius范数图集合F＝{F1,F2,...,Fm,...,FM}中任意p、q时相的极化SAR Frobenius范数图Fp、Fq，得到任意p、q时相间的综合邻域信息差异图DIp,q，其中DIp,q中位置为(x,y)的像素点的差异度值计算公式为：
Figure FDA0003213587360000043
A＝min(Fp(x,y),Fq(x,y))
B＝max(Fp(x,y),Fq(x,y))
Figure FDA0003213587360000044
其中，A′、B′表示A、B当前所在的极化SAR Frobenius范数图的时相值，FilterA′[·]、FilterB′[·]表示利用第A′、B′个时相的平滑因子自适应于综合邻域权重信息的滤波器FilterA′、FilterB′对A、B进行滤波操作，Fp(x,y)、Fq(x,y)表示第p、q个时相的极化SARFrobenius范数图Fp、Fq所在位置为(x,y)的Frobenius范数值，CV″p(x,y)、CV″q(x,y)表示第p、q个时相综合邻域权重信息图CV″p、CV″q所在位置为(x,y)的综合邻域权重信息，CV″p,q(x,y)是CV″p(x,y)与CV″q(x,y)的平均值；
(7)获取基于综合邻域信息差异图的显著性二值检测图：
(7a)利用超像素分割方法，对DIp,q进行超像素分割，得到超像素分割区域集合
Figure FDA0003213587360000051
并采用结合超像素分割的显著性检测方法，通过Areap,q对DIp,q进行显著性检测，得到DIp,q的显著性检测图
Figure FDA0003213587360000052
其中，U≥2，
Figure FDA0003213587360000053
表示第u个超像素分割子区域；
(7b)采用模糊C均值聚类算法，对DIp,q与
Figure FDA0003213587360000054
进行二聚类，并按照聚类结果中每个类别均值由大到小的顺序对两个聚类类别进行排序，得到DIp,q对应的聚类集合
Figure FDA0003213587360000055
以及
Figure FDA0003213587360000056
对应的聚类集合
Figure FDA0003213587360000057
其中
Figure FDA0003213587360000058
表示初始变化类，对应于DIp,q中的初始变化区域
Figure FDA0003213587360000059
Figure FDA00032135873600000510
表示的初始未变化类，对应于DIp,q中的初始未变化区域
Figure FDA00032135873600000511
Figure FDA00032135873600000512
表示非主体显著类，对应于
Figure FDA00032135873600000513
中的非主体显著区域
Figure FDA00032135873600000514
Figure FDA00032135873600000515
表示主体显著类，对应于
Figure FDA00032135873600000516
中的主体显著区域
Figure FDA00032135873600000517
将
Figure FDA00032135873600000518
Figure FDA00032135873600000519
中的所有像素点的像素值置1，同时将
Figure FDA00032135873600000520
Figure FDA00032135873600000521
中的所有像素点的像素值置0；
(7c)将初始变化区域
Figure FDA00032135873600000522
区域内值为1的像素点且主体显著区域
Figure FDA00032135873600000523
内值不为1的像素点对应的区域设置为显著性遗失区域
Figure FDA00032135873600000524
统计每个超像素分割子区域
Figure FDA00032135873600000525
中包含像素值的个数
Figure FDA00032135873600000526
并对显著性遗失区域
Figure FDA00032135873600000527
在每个对应
Figure FDA00032135873600000528
内值为1的像素点个数进行统计，得到统计个数Numu，将Numu大于
Figure FDA00032135873600000529
的对应
Figure FDA00032135873600000530
内的所有像素值设定为1，得到显著性补充区域
Figure FDA00032135873600000531
(7d)对主体显著区域
Figure FDA00032135873600000532
中像素值为1的封闭区域进行统计，得到封闭区域集合L＝{L1,L2,...,Ln,…,LN}，并统计每一个封闭子区域Ln中包含的像素点个数，按照个数进行降序排序，得到降序排序后的封闭子区域Ln中的像素数量集合Num＝{Num1,Num2,…,Numn,…,NumN}，计算相邻两个封闭子区域Ln-1、Ln内像素数量Numn-1、Numn间的比值集合Rio＝{Rio1,Rio2,…,Rion,…,RioN-1}，其中，N≥2，Ln表示第n个封闭区域，Numn表示第n个封闭区域的像素数量，Rion表示Numn-1与Numn的比值；
(7e)统计上述比值集合Rio中的最大值Riomax对应的下标nmax，以下标nmax的前一个下标nmax-1为分隔下标，将L中下标大于nmax-1的封闭子区域划分至主体显著性区域
Figure FDA0003213587360000061
中，将下标小于或等于nmax-1的封闭子区域划分至孤立区域
Figure FDA0003213587360000062
中；
(7f)将孤立区域
Figure FDA0003213587360000063
与显著性补充区域
Figure FDA0003213587360000064
进行求与操作，得到完整显著性补充区域
Figure FDA0003213587360000065
对显著性检测图
Figure FDA0003213587360000066
进行自适应阈值化操作，得到
Figure FDA0003213587360000067
对应的完整主体显著性区域
Figure FDA0003213587360000068
将
Figure FDA0003213587360000069
与
Figure FDA00032135873600000610
进行求与操作，得到综合邻域信息差异图DIp,q的显著性二值检测图
Figure FDA00032135873600000611
(8)获取显著性差异图：
将综合邻域信息差异图DIp,q与显著性二值检测图
Figure FDA00032135873600000612
同一位置像素的值进行相乘，实现对DIp,q的修正，得到显著性差异图
Figure FDA00032135873600000613
其中
Figure FDA00032135873600000614
表示DIp,q和
Figure FDA00032135873600000615
同一位置像素的值相乘；
(9)获取极化SAR图像的变化检测结果：
利用阈值二值化方法对
Figure FDA00032135873600000616
进行二值化分割，得到极化SAR变化检测结果CDp,q。
2.根据权利要求1中所述的基于综合邻域信息的显著性极化SAR图像变化检测方法，其特征在于，步骤(2b)中所述的计算配准后的每个时相的极化SAR图像Tm’的Frobenius范数，其中Tm’中位置为(x,y)的像素点的Frobenius范数Fm(x,y)的计算公式为：
Figure FDA0003213587360000071
其中，||·||F表示求Frobenius范数操作，|·|表示求绝对值操作，t′ij(x,y)表示Tm’(x,y)中第i行第j列的元素值。
3.根据权利要求1中所述的基于综合邻域信息的显著性极化SAR图像变化检测方法，其特征在于，步骤(2c)中所述的计算每个Frobenius范数图Fm的变差系数图，其中每个变差系数图CVm中位置为(x,y)的像素点的变差系数CVm(x,y)的实现步骤为：
以Fm中每个位置为(x,y)的像素点的Frobenius范数Fm(x,y)为中心像素，将Fm划分成设定大小为sizeF×sizeF的邻域窗口
Figure FDA0003213587360000072
sizeF≥2，并计算Fm(x,y)在对应邻域窗
Figure FDA0003213587360000073
内的变差系数CVm(x,y)：
Figure FDA0003213587360000074
其中，Stdm(x,y)是Fm(x,y)在其对应邻域窗
Figure FDA0003213587360000075
内的标准差；μm(x,y)是Fm(x,y)在其对应邻域窗
Figure FDA0003213587360000076
内的均值。
4.根据权利要求1中所述的基于综合邻域信息的显著性极化SAR图像变化检测方法，其特征在于，步骤(5b)中所述的求解每个时相的灰度共生矩阵图Gm的同质度特征图Homm，其中Homm中位置为(x,y)的像素点的同质度特征值Homm(x,y)的计算公式为：
Figure FDA0003213587360000081
其中，gk是每个时相极化SAR Frobenius范数图Fm图的灰度值级数，
Figure FDA0003213587360000082
是灰度共生矩阵图Gm所在位置为(x,y)的像素点的灰度共生矩阵值Gm(x,y)中第l行第s列的元素值。
5.根据权利要求1中所述的基于综合邻域信息的显著性极化SAR图像变化检测方法，其特征在于，步骤(7a)中所述的采用结合超像素分割的显著性检测方法，通过Areap,q对DIp,q进行显著性检测，实现过程为：
统计Areap,q中每个超像素分割子区域
Figure FDA0003213587360000083
内所有像素值的众数
Figure FDA0003213587360000084
以及像素值的总个数
Figure FDA0003213587360000085
并遍历
Figure FDA0003213587360000086
中的每一个像素，得到每个
Figure FDA0003213587360000087
中第一个像素值等于
Figure FDA0003213587360000088
的像素点
Figure FDA0003213587360000089
通过显著性检测方法对每一个
Figure FDA00032135873600000810
进行显著性检测，得到每个
Figure FDA00032135873600000811
对应的显著性检测值
Figure FDA00032135873600000812
将对应
Figure FDA00032135873600000813
中的其余
Figure FDA00032135873600000814
个像素点的显著检测值设定为
Figure FDA00032135873600000815
得到DIp,q对应的显著性检测结果
Figure FDA00032135873600000816
6.根据权利要求1中所述的基于综合邻域信息的显著性极化SAR图像变化检测方法，其特征在于，步骤(7f)中所述的对显著性检测图
Figure FDA00032135873600000817
进行自适应阈值化操作，实现步骤为：
(7f1)设定阈值TCA的初始值为0，以Δ为步长增加TCA；
(7f2)将显著性检测图
Figure FDA00032135873600000818
中显著性值大于TCA的像素点的像素值设定为1，小于TCA的像素点的像素值设定为0，得到
Figure FDA00032135873600000819
对应的显著性阈值图
Figure FDA00032135873600000820
其中，0≤TCA≤1；
(7f3)将主体显著性区域
Figure FDA00032135873600000821
中的像素点与
Figure FDA00032135873600000822
中对应像素点进行求与操作，并统计求与结果中像素值为1的像素点个数numone，同时统计
Figure FDA0003213587360000091
中的像素点个数num′；
(7f4)判断numone＜num′是否成立，若是，则以Δ为步长增加TCA，并执行步骤(7f2)，否则，得到TCA对应的中止阈值Tend，将显著性检测图
Figure FDA0003213587360000092
中显著性值大于Tend的像素点的像素值设定为1，小于TCA的像素点的像素值设定为0，得到
Figure FDA0003213587360000093
对应的完整主体显著性区域
Figure FDA0003213587360000094
CN202110936927.6A 2021-08-16 2021-08-16 基于综合邻域信息的显著性极化sar图像变化检测方法 Active CN113628234B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110936927.6A CN113628234B (zh) 2021-08-16 2021-08-16 基于综合邻域信息的显著性极化sar图像变化检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110936927.6A CN113628234B (zh) 2021-08-16 2021-08-16 基于综合邻域信息的显著性极化sar图像变化检测方法

Publications (2)

Publication Number Publication Date
CN113628234A true CN113628234A (zh) 2021-11-09
CN113628234B CN113628234B (zh) 2024-04-19

Family

ID=78385771

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110936927.6A Active CN113628234B (zh) 2021-08-16 2021-08-16 基于综合邻域信息的显著性极化sar图像变化检测方法

Country Status (1)

Country Link
CN (1) CN113628234B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114898224A (zh) * 2022-05-13 2022-08-12 北京科技大学 一种基于物理散射机制的变化检测方法
CN115797299A (zh) * 2022-12-05 2023-03-14 常宝新材料(苏州)有限公司 一种光学复合膜的缺陷检测方法
CN117333468A (zh) * 2023-10-17 2024-01-02 南京北斗创新应用科技研究院有限公司 面向多模式时序PolSAR影像的洪涝灾害监测方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20060007816A (ko) * 2004-07-22 2006-01-26 학교법인 중앙대학교 영상 합성 방법
CN101982835A (zh) * 2010-11-12 2011-03-02 西安电子科技大学 Sar图像机场道路边缘检测水平集方法
CN102208109A (zh) * 2011-06-23 2011-10-05 南京林业大学 X射线图像和激光图像的异源图像配准方法
CN102938071A (zh) * 2012-09-18 2013-02-20 西安电子科技大学 基于非局部均值的sar图像变化检测模糊聚类分析方法
US20160379053A1 (en) * 2015-06-23 2016-12-29 University Of Electronic Science And Technology Of China Method and Apparatus for Identifying Object

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20060007816A (ko) * 2004-07-22 2006-01-26 학교법인 중앙대학교 영상 합성 방법
CN101982835A (zh) * 2010-11-12 2011-03-02 西安电子科技大学 Sar图像机场道路边缘检测水平集方法
CN102208109A (zh) * 2011-06-23 2011-10-05 南京林业大学 X射线图像和激光图像的异源图像配准方法
CN102938071A (zh) * 2012-09-18 2013-02-20 西安电子科技大学 基于非局部均值的sar图像变化检测模糊聚类分析方法
US20160379053A1 (en) * 2015-06-23 2016-12-29 University Of Electronic Science And Technology Of China Method and Apparatus for Identifying Object

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
金杰;吴雅男;康仲林;: "基于多特征的SAR影像溢油暗斑提取", 测绘与空间地理信息, no. 02 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114898224A (zh) * 2022-05-13 2022-08-12 北京科技大学 一种基于物理散射机制的变化检测方法
CN115797299A (zh) * 2022-12-05 2023-03-14 常宝新材料(苏州)有限公司 一种光学复合膜的缺陷检测方法
CN115797299B (zh) * 2022-12-05 2023-09-01 常宝新材料(苏州)有限公司 一种光学复合膜的缺陷检测方法
CN117333468A (zh) * 2023-10-17 2024-01-02 南京北斗创新应用科技研究院有限公司 面向多模式时序PolSAR影像的洪涝灾害监测方法
CN117333468B (zh) * 2023-10-17 2024-05-24 南京北斗创新应用科技研究院有限公司 面向多模式时序PolSAR影像的洪涝灾害监测方法

Also Published As

Publication number Publication date
CN113628234B (zh) 2024-04-19

Similar Documents

Publication Publication Date Title
CN113628234A (zh) 基于综合邻域信息的显著性极化sar图像变化检测方法
CN107833220B (zh) 基于深度卷积神经网络与视觉显著性的织物缺陷检测方法
CN110930416B (zh) 一种基于u型网络的mri图像前列腺分割方法
CN108765465B (zh) 一种无监督sar图像变化检测方法
CN106875395B (zh) 基于深度神经网络的超像素级sar图像变化检测方法
CN108549891A (zh) 基于背景与目标先验的多尺度扩散显著目标检测方法
CN110309781B (zh) 基于多尺度光谱纹理自适应融合的房屋损毁遥感识别方法
Zhang et al. Regions of interest detection in panchromatic remote sensing images based on multiscale feature fusion
CN109978848A (zh) 基于多光源颜色恒常模型检测眼底图像中硬性渗出的方法
CN104102928B (zh) 一种基于纹理基元的遥感图像分类方法
CN109712149B (zh) 一种基于小波能量和模糊c均值的图像分割方法
CN107169962B (zh) 基于空间密度约束核模糊聚类的灰度图像快速分割方法
CN109146890A (zh) 基于滤波器的高光谱图像的异常目标检测方法
CN117764864B (zh) 基于图像去噪的核磁共振肿瘤视觉检测方法
Wang et al. The PAN and MS image fusion algorithm based on adaptive guided filtering and gradient information regulation
CN111291615A (zh) 一种多时相遥感影像变化监测方法
CN107392211B (zh) 基于视觉稀疏认知的显著目标检测方法
CN115690086A (zh) 一种基于对象的高分辨率遥感影像变化检测方法及系统
CN111626380A (zh) 一种基于超像素和卷积网络的极化sar图像分类方法
Wu et al. A modified fuzzy dual-local information c-mean clustering algorithm using quadratic surface as prototype for image segmentation
CN107944497A (zh) 基于主成分分析的图像块相似性度量方法
CN107346549B (zh) 一种利用遥感影像多特征的多类别变化动态阈值检测方法
CN112651951A (zh) 一种基于dce-mri的乳腺癌分级方法
CN110211106B (zh) 基于分段Sigmoid带宽的均值漂移SAR图像海岸线检测方法
Shen et al. Cropland extraction from very high spatial resolution satellite imagery by object-based classification using improved mean shift and one-class support vector machines

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