CN113628234B - 基于综合邻域信息的显著性极化sar图像变化检测方法 - Google Patents
基于综合邻域信息的显著性极化sar图像变化检测方法 Download PDFInfo
- Publication number
- CN113628234B CN113628234B CN202110936927.6A CN202110936927A CN113628234B CN 113628234 B CN113628234 B CN 113628234B CN 202110936927 A CN202110936927 A CN 202110936927A CN 113628234 B CN113628234 B CN 113628234B
- Authority
- CN
- China
- Prior art keywords
- pixel
- information
- map
- 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.)
- Active
Links
- 230000008859 change Effects 0.000 title claims abstract description 63
- 238000000034 method Methods 0.000 title claims abstract description 54
- 230000010287 polarization Effects 0.000 title claims description 24
- 238000001514 detection method Methods 0.000 claims abstract description 88
- 238000009499 grossing Methods 0.000 claims abstract description 32
- 230000011218 segmentation Effects 0.000 claims abstract description 22
- 239000011159 matrix material Substances 0.000 claims description 21
- 230000003044 adaptive effect Effects 0.000 claims description 16
- 238000000354 decomposition reaction Methods 0.000 claims description 16
- 238000012937 correction Methods 0.000 claims description 15
- 238000010586 diagram Methods 0.000 claims description 13
- 238000004364 calculation method Methods 0.000 claims description 11
- 238000001914 filtration Methods 0.000 claims description 11
- 230000015572 biosynthetic process Effects 0.000 claims description 9
- 238000003786 synthesis reaction Methods 0.000 claims description 9
- 238000004422 calculation algorithm Methods 0.000 claims description 7
- 238000003708 edge detection Methods 0.000 claims description 6
- 230000016507 interphase Effects 0.000 claims description 6
- 238000012163 sequencing technique Methods 0.000 claims description 6
- 238000005829 trimerization reaction Methods 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 4
- 230000000295 complement effect Effects 0.000 claims description 3
- 238000005192 partition Methods 0.000 claims description 3
- 238000000926 separation method Methods 0.000 claims description 3
- 239000013589 supplement Substances 0.000 claims description 3
- 230000000153 supplemental effect Effects 0.000 claims description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 11
- 230000007547 defect Effects 0.000 description 8
- 238000004088 simulation Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 4
- 230000014759 maintenance of location Effects 0.000 description 3
- 230000000717 retained effect Effects 0.000 description 3
- 230000001629 suppression Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000006870 function Effects 0.000 description 2
- 230000005764 inhibitory process Effects 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 239000002131 composite material Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000012938 design process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000005571 horizontal transmission Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000036039 immunity Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000008447 perception Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000005570 vertical transmission Effects 0.000 description 1
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通过交替的水平发射和接收以及垂直发射和接收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对应的三聚类集合其中,/>分别表示CVm的三聚类结果中的高、中、低同质度类别;
(3b)根据变差系数图CVm的三聚类结果中的高同质度类别与中同质度类别/>对应区域之间的分界值/>以及中同质度类别/>和低同质度类别/>对应区域之间的分界值/>计算中同质度类别/>对应区域内和低同质度类别/>对应区域内的变差系数信息动态值/>
其中是中同质度类别/>对应区域内和低同质度类别/>对应区域内Em(x,y)为1的像素点个数占对应区域内像素点总数的比例;
(3c)通过变差系数图CVm中所有像素点的变差系数信息的最大值中同质度类别/>和低同质度类别/>对应区域内的变差系数动态值/>对初始边缘信息图Em中所有位置为(x,y)的像素的边缘信息Em(x,y)进行修正,得到修正后的边缘信息图集合E′={E′1,E′2,...,E′m,...,E′M},其中Em(x,y)的修正公式为:
其中,CVm(x,y)表示变差系数图CVm所在位置为(x,y)的像素的变差系数信息,E′m(x,y)表示Em(x,y)的修正结果;
(4)通过修正后的边缘信息对变差系数信息进行指导:
(4a)通过变差系数图CVm中高同质度类别对应区域内所有像素点的变差系数信息的平均值/>以及修正后的边缘信息图E′m中所有像素点的边缘信息的最大值/>对每个时相的变差系数图CVm中所有位置为(x,y)的变差系数CVm(x,y)进行指导,得到指导后的初始综合邻域权重信息图集合CV′={CV′1,CV′2,...,CV′m,...,CV′M},其中CVm(x,y)的修正公式为:
其中,E′m(x,y)表示修正后的边缘信息图E′m所在位置为(x,y)的像素的边缘信息,CV′m(x,y)表示初始综合邻域权重信息图CV′m所在位置为(x,y)的初始综合邻域权重信息;
(4b)对每个时相的初始综合邻域权重信息图CV′m进行归一化,得到范围在修正后的边缘信息图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中在高同质度类别对应区域/>中所有像素点的平均值/>
(5c)通过综合邻域权重信息图CV″m中所有像素点对应的综合邻域权重信息的最大值/>对每个时相的初始平滑因子图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)))
其中,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)的像素点的差异度值计算公式为:
A=min(Fp(x,y),Fq(x,y))
B=max(Fp(x,y),Fq(x,y))
其中,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进行超像素分割,得到超像素分割区域集合并采用结合超像素分割的显著性检测方法,通过Areap,q对DIp,q进行显著性检测,得到DIp,q的显著性检测图/>其中,U≥2,表示第u个超像素分割子区域;
(7b)采用模糊C均值聚类算法,对DIp,q与进行二聚类,并按照聚类结果中每个类别均值由大到小的顺序对两个聚类类别进行排序,得到DIp,q对应的聚类集合以及/>对应的聚类集合/>其中/>表示初始变化类,对应于DIp,q中的初始变化区域/>表示的初始未变化类,对应于DIp,q中的初始未变化区域/>表示非主体显著类,对应于/>中的非主体显著区域/> 表示主体显著类,对应于/>中的主体显著区域/>将/>中的所有像素点的像素值置1,同时将/>中的所有像素点的像素值置0;
(7c)将初始变化区域区域内值为1的像素点且主体显著区域/>内值不为1的像素点对应的区域设置为显著性遗失区域/>统计每个超像素分割子区域/>中包含像素值的个数/>并对显著性遗失区域/>在每个对应/>内值为1的像素点个数进行统计,得到统计个数Numu,将Numu大于/>的对应/>内的所有像素值设定为1,得到显著性补充区域/>
(7d)对主体显著区域中像素值为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的封闭子区域划分至主体显著性区域中,将下标小于或等于nmax-1的封闭子区域划分至孤立区域/>中;
(7f)将孤立区域与显著性补充区域/>进行求与操作,得到完整显著性补充区域/>对显著性检测图/>进行自适应阈值化操作,得到/>对应的完整主体显著性区域/>将/>与/>进行求与操作,得到综合邻域信息差异图DIp,q的显著性二值检测图/>
(8)获取显著性差异图:
将综合邻域信息差异图DIp,q与显著性二值检测图同一位置像素的值进行相乘,实现对DIp,q的修正,得到显著性差异图/>其中/>表示DIp,q和/>同一位置像素的值相乘;
(9)获取极化SAR图像的变化检测结果:
利用阈值二值化方法对进行二值化分割,得到极化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分解:
其中,[Sm(x,y)]表示极化SAR散射矩阵,表示接收方式为i′、发射方式为j′的散射矩阵,H表示水平方式,V表示垂直方式,/>表示Pauli基,a′,b′,c′和d′都是复数,当介质满足互易条件时,可将这些量写为向量Km的形式:
Pauli基是完备正交基,具有一定的抗噪性,|a′m|2,|b′m|2和|c′m|2,对应着明显的物理散射机制,所以利用Pauli分解各系数,可以合成任意时相RGB图像Pm:
其中,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)的计算公式为:
其中,||·||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的邻域窗口本实施例中,sizeF=7,并计算Fm(x,y)在对应邻域窗/>内的变差系数CVm(x,y):
其中,Stdm(x,y)是Fm(x,y)在其对应邻域窗内的标准差;μm(x,y)是Fm(x,y)在其对应邻域窗/>内的均值;
步骤3通过变差系数信息对初始边缘信息进行修正:
(3a)采用模糊C均值聚类算法,根据变差系数图CVm中每个像素点的变差系数信息对CVm进行三聚类,并按照每个像素簇内变差系数信息均值由小到大的顺序对CVm的三个聚类类别进行排序,得到CV对应的三聚类集合其中,/>分别表示CVm的三聚类结果中的高、中、低同质度类别;
(3b)根据变差系数图CVm的三聚类结果中的高同质度类别与中同质度类别/>对应区域之间的分界值/>以及中同质度类别/>和低同质度类别/>对应区域之间的分界值计算中同质度类别/>对应区域内和低同质度类别/>对应区域内的变差系数信息动态值/>
其中是中同质度类别/>对应区域内和低同质度类别/>对应区域内Em(x,y)为1的像素点个数占对应区域内像素点总数的比例;
(3c)通过变差系数图CVm中所有像素点的变差系数信息的最大值中同质度类别/>和低同质度类别/>对应区域内的变差系数动态值/>对初始边缘信息图Em中所有位置为(x,y)的像素的边缘信息Em(x,y)进行修正,得到修正后的边缘信息图集合E′={E′1,E′2,...,E′m,...,E′M},其中Em(x,y)的修正公式为:
其中,CVm(x,y)表示变差系数图CVm所在位置为(x,y)的像素的变差系数信息,E′m(x,y)表示Em(x,y)的修正结果;
步骤4通过修正后的边缘信息对变差系数信息进行指导:
(4a)通过变差系数图CVm中高同质度类别对应区域内所有像素点的变差系数信息的平均值/>以及修正后的边缘信息图E′m中所有像素点的边缘信息的最大值/>对每个时相的变差系数图CVm中所有位置为(x,y)的变差系数CVm(x,y)进行指导,得到指导后的初始综合邻域权重信息图集合CV′={CV′1,CV′2,...,CV′m,...,CV′M},其中CVm(x,y)的修正公式为:
其中,E′m(x,y)表示修正后的边缘信息图E′m所在位置为(x,y)的像素的边缘信息,CV′m(x,y)表示初始综合邻域权重信息图CV′m所在位置为(x,y)的初始综合邻域权重信息;
(4b)对每个时相的初始综合邻域权重信息图CV′m进行归一化,得到范围在修正后的边缘信息图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)的计算公式为:
其中,gk是每个时相极化SAR Frobenius范数图Fm图的灰度值级数,是灰度共生矩阵图Gm所在位置为(x,y)的像素点的灰度共生矩阵值Gm(x,y)中第l行第s列的元素值;本实例中,gk=16;
在得到Homm之后,进一步统计Homm中在高同质度类别对应区域/>中所有像素点的平均值/>
(5c)通过综合邻域权重信息图CV″m中所有像素点对应的综合邻域权重信息的最大值/>对每个时相的初始平滑因子图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)))
其中,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;
上述步骤中,非局部均值滤波器的具体设计方法为:
其中,B′(ac,rr)表示一个以ac为中心,大小为(2rr+1)×(2rr+1)像素的邻域,Norm(ac)表示归一化因子,u′是滤波前的值,是滤波后的值,本实施例中,rr=3;
权重w′(ac,a′c)取决于分别以ac和a′c为中心的(2rr+1)×(2rr+1)大小的彩色patch块的欧氏距离dis2(B′(ac,rr),B′(a′c,rr))的平方:
利用指数核函数并结合自适应平滑因子h′来计算w′(ac,a′c):
式中,σ是噪声的标准差,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)的像素点的差异度值计算公式为:
A=min(Fp(x,y),Fq(x,y))
B=max(Fp(x,y),Fq(x,y))
其中,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进行超像素分割,得到超像素分割区域集合表示第u个超像素分割子区域,本实例中超像素分割方法采用效果优良且应用广泛的简单线性迭代聚类(SLIC)方法,其中初始分割子区域的窗口尺寸为rslic=27,并采用结合超像素分割的显著性检测方法,通过对DIp,q进行显著性检测,实现过程为:
统计Areap,q中每个超像素分割子区域内所有像素值的众数/>以及像素值的总个数/>并遍历/>中的每一个像素,得到每个/>中第一个像素值等于/>的像素点/>通过显著性检测方法对每一个/>进行显著性检测,得到每个/>对应的显著性检测值/>将对应/>中的其余/>个像素点的显著检测值设定为/>得到DIp,q对应的显著性检测结果/>
本实施例中,显著性检测方法选择了应用广泛、效果良好的上下文感知显著性检测(CA)方法,具体设计过程为:
首先获得一个衡量显著性的综合距离,上下文感知显著性检测方法对图像块之间在Lab颜色空间做对比,若某一块与其它块差距大,则说明是显著性特征。考虑空间距离,与显著性区域相似的区域一般都离得比较近,据此可以得到综合距离:
其中,Patchci、Patchcj分别表示中心像素为ci、cj的像素块,dcolor(·)代表两点间的Lab颜色距离,dposition(·)代表两点间的空间距离,一对图像块的综合距离越大,差异度就越大,若对于任意的块,得到的差异值都很大,则可判定该块为显著块;本实施例中,para=3;
进一步确定显著性公式,为降低算法的计算复杂度,只计算前K′个与最相似的块即可,得到的显著性公式为:
其中,ci,k′分别像素块的中心像素,r′表示尺度值,/>表示像素ci处尺度为r′的显著性值,本实施例中,K′=65;
将上式的单个尺度所得到的显著性值引申为求多个尺度下的显著性值的平均值来进一步提高显著与非显著区域的对比度,即:
式中,表示像素ci处的多尺度平均显著值,M′为尺度的个数,R′={r′1,...,r′M};本实施例中,M′=4;
进一步加入上下文修正,设定一个显著性阈值TCA,并从显著图中提取注意力最大的局部区域,在注意力区域之外的像素显著性值由与其最近的注意力像素之间的欧式距离加权,得到新的显著性值:
式中,dfoci(ci)为离像素点ci最近的注意力像素之间的欧式距离,表示像素ci处的多尺度平均显著值,/>表示像素ci处的最终显著性值;本实施例中,TCA=0.8;
(7b)采用模糊C均值聚类算法,对DIp,q与进行二聚类,并按照聚类结果中每个类别均值由大到小的顺序对两个聚类类别进行排序,得到DIp,q对应的聚类集合/>以及/>对应的聚类集合/>其中/>表示初始变化类,对应于DIp,q中的初始变化区域/>表示的初始未变化类,对应于DIp,q中的初始未变化区域/>表示非主体显著类,对应于/>中的非主体显著区域/> 表示主体显著类,对应于/>中的主体显著区域/>将/>中的所有像素点的像素值置1,同时将/>中的所有像素点的像素值置0;
(7c)将初始变化区域区域内值为1的像素点且主体显著区域/>内值不为1的像素点对应的区域设置为显著性遗失区域/>统计每个超像素分割子区域/>中包含像素值的个数/>并对显著性遗失区域/>在每个对应/>内值为1的像素点个数进行统计,得到统计个数Numu,将Numu大于/>的对应/>内的所有像素值设定为1,得到显著性补充区域/>
(7d)对主体显著区域中像素值为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的封闭子区域划分至主体显著性区域中,将下标小于或等于nmax-1的封闭子区域划分至孤立区域/>中;
(7f)将孤立区域与显著性补充区域/>进行求与操作,得到完整显著性补充区域/>对显著性检测图/>进行自适应阈值化操作,得到/>对应的完整主体显著性区域/>实现步骤为:
(7f1)设定阈值TCA的初始值为0,以Δ为步长增加TCA;本实施例中,Δ=0.01;
(7f2)将显著性检测图中显著性值大于TCA的像素点的像素值设定为1,小于TCA的像素点的像素值设定为0,得到/>对应的显著性阈值图/>其中,0≤TCA≤1;
(7f3)将主体显著性区域中的像素点与/>中对应像素点进行求与操作,并统计求与结果中像素值为1的像素点个数numone,同时统计/>中的像素点个数num′;
(7f4)判断numone<num′是否成立,若是,则以Δ为步长增加TCA,并执行步骤(7f2),否则,得到TCA对应的中止阈值Tend,将显著性检测图中显著性值大于Tend的像素点的像素值设定为1,小于TCA的像素点的像素值设定为0,得到/>对应的完整主体显著性区域
之后,将与/>进行求与操作,得到综合邻域信息差异图DIp,q的显著性二值检测图/>
步骤8获取显著性差异图:
将综合邻域信息差异图DIp,q与显著性二值检测图同一位置像素的值进行相乘,实现对DIp,q的修正,得到显著性差异图/>其中/>表示DIp,q和/>同一位置像素的值相乘;
步骤9获取极化SAR图像的变化检测结果:
利用阈值二值化方法对进行二值化分割,得到极化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
参照表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,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对应的三聚类集合其中,/>分别表示CVm的三聚类结果中的高、中、低同质度类别;
(3b)根据变差系数图CVm的三聚类结果中的高同质度类别与中同质度类别/>对应区域之间的分界值/>以及中同质度类别/>和低同质度类别/>对应区域之间的分界值计算中同质度类别/>对应区域内和低同质度类别/>对应区域内的变差系数信息动态值/>
其中是中同质度类别/>对应区域内和低同质度类别/>对应区域内Em(x,y)为1的像素点个数占对应区域内像素点总数的比例;
(3c)通过变差系数图CVm中所有像素点的变差系数信息的最大值中同质度类别和低同质度类别/>对应区域内的变差系数动态值/>对初始边缘信息图Em中所有位置为(x,y)的像素的边缘信息Em(x,y)进行修正,得到修正后的边缘信息图集合E′={E′1,E′2,...,E′m,...,E′M},其中Em(x,y)的修正公式为:
其中,CVm(x,y)表示变差系数图CVm所在位置为(x,y)的像素的变差系数信息,E′m(x,y)表示Em(x,y)的修正结果;
(4)通过修正后的边缘信息对变差系数信息进行指导:
(4a)通过变差系数图CVm中高同质度类别对应区域内所有像素点的变差系数信息的平均值/>以及修正后的边缘信息图E′m中所有像素点的边缘信息的最大值/>对每个时相的变差系数图CVm中所有位置为(x,y)的变差系数CVm(x,y)进行指导,得到指导后的初始综合邻域权重信息图集合CV′={CV1′,CV′2,…,CV′m,...,CV′M},其中CVm(x,y)的修正公式为:
其中,E′m(x,y)表示修正后的边缘信息图E′m所在位置为(x,y)的像素的边缘信息,CV′m(x,y)表示初始综合邻域权重信息图CV′m所在位置为(x,y)的初始综合邻域权重信息;
(4b)对每个时相的初始综合邻域权重信息图CV′m进行归一化,得到范围在修正后的边缘信息图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中在高同质度类别对应区域/>中所有像素点的平均值/>
(5c)通过综合邻域权重信息图CV″m中所有像素点对应的综合邻域权重信息的最大值/>对每个时相的初始平滑因子图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)))
其中,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时相的极化SARFrobenius范数图Fp、Fq,得到任意p、q时相间的综合邻域信息差异图DIp,q,其中DIp,q中位置为(x,y)的像素点的差异度值计算公式为:
A=min(Fp(x,y),Fq(x,y))
B=max(Fp(x,y),Fq(x,y))
其中,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进行超像素分割,得到超像素分割区域集合并采用结合超像素分割的显著性检测方法,通过Areap,q对DIp,q进行显著性检测,得到DIp,q的显著性检测图/>其中,U≥2,表示第u个超像素分割子区域;
(7b)采用模糊C均值聚类算法,对DIp,q与进行二聚类,并按照聚类结果中每个类别均值由大到小的顺序对两个聚类类别进行排序,得到DIp,q对应的聚类集合/>以及/>对应的聚类集合/>其中/>表示初始变化类,对应于DIp,q中的初始变化区域/>表示的初始未变化类,对应于DIp,q中的初始未变化区域/>表示非主体显著类,对应于/>中的非主体显著区域/>表示主体显著类,对应于中的主体显著区域/>将/>中的所有像素点的像素值置1,同时将/> 中的所有像素点的像素值置0;
(7c)将初始变化区域区域内值为1的像素点且主体显著区域/>内值不为1的像素点对应的区域设置为显著性遗失区域/>统计每个超像素分割子区域/>中包含像素值的个数/>并对显著性遗失区域/>在每个对应/>内值为1的像素点个数进行统计,得到统计个数Numu,将Numu大于/>的对应/>内的所有像素值设定为1,得到显著性补充区域/>
(7d)对主体显著区域中像素值为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的封闭子区域划分至主体显著性区域中,将下标小于或等于nmax-1的封闭子区域划分至孤立区域/>中;
(7f)将孤立区域与显著性补充区域/>进行求与操作,得到完整显著性补充区域/>对显著性检测图/>进行自适应阈值化操作,得到/>对应的完整主体显著性区域/>将/>与/>进行求与操作,得到综合邻域信息差异图DIp,q的显著性二值检测图/>
(8)获取显著性差异图:
将综合邻域信息差异图DIp,q与显著性二值检测图同一位置像素的灰度值进行相乘,实现对DIp,q的修正,得到显著性差异图/>其中DotM(·)表示矩阵间同一位置像素值进行相乘运算;
(9)获取极化SAR图像的变化检测结果:
利用阈值二值化方法对进行二值化分割,得到极化SAR变化检测结果CDp,q。
2.根据权利要求1中所述的基于综合邻域信息的显著性极化SAR图像变化检测方法,其特征在于,步骤(2b)中所述的计算配准后的每个时相的极化SAR图像T’m的Frobenius范数,其中T’m中位置为(x,y)的像素点的Frobenius范数Fm(x,y)的计算公式为:
其中,||·||F表示求Frobenius范数操作,|·|表示求绝对值操作,t′ij(x,y)表示T’m(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的邻域窗口sizeF≥2,并计算Fm(x,y)在对应邻域窗/>内的变差系数CVm(x,y):
其中,Stdm(x,y)是Fm(x,y)在其对应邻域窗内的标准差;μm(x,y)是Fm(x,y)在其对应邻域窗/>内的均值。
4.根据权利要求1中所述的基于综合邻域信息的显著性极化SAR图像变化检测方法,其特征在于,步骤(5b)中所述的求解每个时相的灰度共生矩阵图Gm的同质度特征图Homm,其中Homm中位置为(x,y)的像素点的同质度特征值Homm(x,y)的计算公式为:
其中,gk是每个时相极化SAR Frobenius范数图Fm图的灰度值级数,是灰度共生矩阵图Gm所在位置为(x,y)的像素点的灰度共生矩阵值Gm(x,y)中第l行第s列的元素值。
5.根据权利要求1中所述的基于综合邻域信息的显著性极化SAR图像变化检测方法,其特征在于,步骤(7a)中所述的采用结合超像素分割的显著性检测方法,通过Areap,q对DIp,q进行显著性检测,实现过程为:
统计Areap,q中每个超像素分割子区域内所有像素值的众数/>以及像素值的总个数numAreap,q,并遍历/>中的每一个像素,得到每个/>中第一个像素值等于的像素点/>通过显著性检测方法对每一个/>进行显著性检测,得到每个对应的显著性检测值/>将对应/>中的其余/>个像素点的显著检测值设定为/>得到DIp,q对应的显著性检测结果/>
6.根据权利要求1中所述的基于综合邻域信息的显著性极化SAR图像变化检测方法,其特征在于,步骤(7f)中所述的对显著性检测图进行自适应阈值化操作,实现步骤为:
(7f1)设定阈值TCA的初始值为0,以Δ为步长增加TCA;
(7f2)将显著性检测图中显著性值大于TCA的像素点的像素值设定为1,小于TCA的像素点的像素值设定为0,得到/>对应的显著性阈值图/>其中,0≤TCA≤1;
(7f3)将主体显著性区域中的像素点与/>中对应像素点进行求与操作,并统计求与结果中像素值为1的像素点个数numone,同时统计/>中的像素点个数num′;
(7f4)判断numone<num′是否成立,若是,则以Δ为步长增加TCA,并执行步骤(7f2),否则,得到TCA对应的中止阈值Tend,将显著性检测图中显著性值大于Tend的像素点的像素值设定为1,小于TCA的像素点的像素值设定为0,得到/>对应的完整主体显著性区域
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 CN113628234A (zh) | 2021-11-09 |
CN113628234B true 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) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114898224B (zh) * | 2022-05-13 | 2023-04-21 | 北京科技大学 | 一种基于物理散射机制的变化检测方法 |
CN115797299B (zh) * | 2022-12-05 | 2023-09-01 | 常宝新材料(苏州)有限公司 | 一种光学复合膜的缺陷检测方法 |
CN117333468B (zh) * | 2023-10-17 | 2024-05-24 | 南京北斗创新应用科技研究院有限公司 | 面向多模式时序PolSAR影像的洪涝灾害监测方法 |
Citations (4)
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图像变化检测模糊聚类分析方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104966065B (zh) * | 2015-06-23 | 2018-11-09 | 电子科技大学 | 目标识别方法及装置 |
-
2021
- 2021-08-16 CN CN202110936927.6A patent/CN113628234B/zh active Active
Patent Citations (4)
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图像变化检测模糊聚类分析方法 |
Non-Patent Citations (1)
Title |
---|
基于多特征的SAR影像溢油暗斑提取;金杰;吴雅男;康仲林;;测绘与空间地理信息(第02期);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN113628234A (zh) | 2021-11-09 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN113628234B (zh) | 基于综合邻域信息的显著性极化sar图像变化检测方法 | |
Chierchia et al. | Multitemporal SAR image despeckling based on block-matching and collaborative filtering | |
Cross et al. | Segmentation of remotely-sensed images by a split-and-merge process+ | |
Khazai et al. | An approach for subpixel anomaly detection in hyperspectral images | |
CN110927706A (zh) | 基于卷积神经网络的雷达干扰检测识别方法 | |
Zhang et al. | Regions of interest detection in panchromatic remote sensing images based on multiscale feature fusion | |
CN110097101B (zh) | 一种基于改进可靠性因子的遥感图像融合与海岸带分类方法 | |
CN105718957A (zh) | 非下采样轮廓波卷积神经网络的极化sar图像分类方法 | |
Modava et al. | Hierarchical coastline detection in SAR images based on spectral‐textural features and global–local information | |
Xiang et al. | Adaptive statistical superpixel merging with edge penalty for PolSAR image segmentation | |
CN112508963B (zh) | 一种基于模糊c均值聚类的sar图像分割方法 | |
CN106295498A (zh) | 光学遥感图像目标区域检测装置与方法 | |
CN109360191B (zh) | 一种基于变分自编码器的图像显著性检测方法 | |
Wang et al. | The PAN and MS image fusion algorithm based on adaptive guided filtering and gradient information regulation | |
CN107392211B (zh) | 基于视觉稀疏认知的显著目标检测方法 | |
CN107464247B (zh) | 一种基于g0分布的随机梯度变分贝叶斯sar图像分割方法 | |
Hao et al. | A novel change detection approach for VHR remote sensing images by integrating multi-scale features | |
CN108509835B (zh) | 基于DFIC超像素的PolSAR图像地物分类方法 | |
Kusetogullari et al. | Unsupervised change detection in landsat images with atmospheric artifacts: a fuzzy multiobjective approach | |
Chen et al. | Change detection of multispectral remote-sensing images using stationary wavelet transforms and integrated active contours | |
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 | |
CN115937302A (zh) | 结合边缘保持的高光谱图像亚像元定位方法 | |
CN108764145A (zh) | 一种面向扎龙湿地遥感图像密度峰值聚类方法 | |
CN109697466A (zh) | 一种自适应区间型空间模糊c均值的地物分类方法 | |
CN110751652B (zh) | 基于巴氏距离和纹理模式度量的sar图像分割方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |