CN112348750B - 基于阈值融合和邻域投票的sar图像变化检测方法 - Google Patents

基于阈值融合和邻域投票的sar图像变化检测方法 Download PDF

Info

Publication number
CN112348750B
CN112348750B CN202011160475.9A CN202011160475A CN112348750B CN 112348750 B CN112348750 B CN 112348750B CN 202011160475 A CN202011160475 A CN 202011160475A CN 112348750 B CN112348750 B CN 112348750B
Authority
CN
China
Prior art keywords
class
pixel
pixels
threshold
unchanged
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
CN202011160475.9A
Other languages
English (en)
Other versions
CN112348750A (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 CN202011160475.9A priority Critical patent/CN112348750B/zh
Publication of CN112348750A publication Critical patent/CN112348750A/zh
Application granted granted Critical
Publication of CN112348750B publication Critical patent/CN112348750B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • 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
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/40—Image enhancement or restoration using histogram techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • 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
    • 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)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种基于阈值融合和邻域投票的SAR图像变化检测方法，其实现步骤为：对待检测的两个不同时相的SAR图像进行预处理；获得差异图；获得融合阈值；获取预分类结果图；对中间类进行邻域投票分类；更新预分类结果图直至更新后与更新前预分类结果图的中间类像素数量相等；判断更新后的预分类结果图中中间类像素数是否为零，若是则检测结束，否则继续下一步骤；对初步邻域投票图进行扩大邻域投票分类；得到分类结果图。本发明采用阈值融合方法得到中间类像素，结合难分类像素的空间位置信息对预分类后的中间类像素进行邻域投票分类，提高了变化检测的准确率。

Description

基于阈值融合和邻域投票的SAR图像变化检测方法
技术领域
本发明属于图像处理技术领域，更进一步涉及图像分类识别技术领域中的一种基于阈值融合和邻域投票的合成孔径雷达(SAR，Synthetic Aperture Radar)图像变化检测方法。该方法可用于检测两时相合成孔径雷达图像的变化区域，完成对洪涝灾害区域的监测以及农业区域调查等SAR图像变化检测任务。
背景技术
合成孔径雷达传感器对天气条件和光照条件不敏感，可以在各种气候条件和不同时段工作。许多星载SAR传感器的重访周期很短，SAR图像数据越来越容易获取，目前，SAR已广泛应用于资源探测、海洋监测、军事侦察以及科学研究等领域，具有重要的研究价值和应用前景。SAR图像变化检测主要可分为有监督的SAR图像变化检测方法和无监督的SAR图像变化检测方法。因为在实际应用中，有标签的SAR图像训练数据难以获得，因此无监督的SAR图像变化检测方法成为SAR图像变化检测方向的研究重点，无监督SAR图像变化检测方法可以在没有标签样本的情况下完成SAR图像变化检测任务，更加符合实际应用需求。
中国海洋大学2016年在IEEE Geoscience and Remote Sensing Letters上发表的"Automatic Change Detection in Synthetic Aperture Radar Images Based onPCANet,"(参考文献：F.Gao,J.Dong,B.Li and Q.Xu,"Automatic Change Detection inSynthetic Aperture Radar Images Based on PCANet,"in IEEE Geoscience andRemote Sensing Letters,vol.13,no.12,pp.1792-1796,Dec.2016,doi:10.1109/LGRS.2016.2611001.)中提出了一种基于PCANet的SAR图像自动变化检测方法。该方法的实现过程为：首先通过FCM聚类算法对对数比差异图进行预分类，得到含有变化类、未变化类和中间类三类像素的预分类结果，并选择一定数量的变化类和未变化类像素组成训练样本去训练PCANet模型，然后用训练好的模型对待分类的样本进行分类得到最后的SAR图像变化检测二值图。该方法虽然在一定程度上提高了变化检测的准确率，但是，该方法仍然存在的不足之处是，首先简单的对数比差异图像不能准确的反映获取的两个时相SAR图像的变化信息，从而不能为变化检测任务的下一个步骤打好基础，影响变化检测准确率的提升；其次，该方法通过FCM(Fuzzy C-means)预分类方法为PCANet模型的训练提供伪标签训练样本，该预分类方法利用差异图进行预分类得到的伪标签可能不够准确，从而导致变化检测的准确率仍不够高。
Y.Bazi，L.Bruzzone和F.Melgani等人在IEEE Transactions on Geoscience andRemote Sensing上发表的"An unsupervised approach based on the generalizedGaussian model to automatic change detection in multitemporal SAR images,"(参考文献：Y.Bazi,L.Bruzzone and F.Melgani,"An unsupervised approach based on thegeneralized Gaussian model to automatic change detection in multitemporal SARimages,"in IEEE Transactions on Geoscience and Remote Sensing,vol.43,no.4,pp.874-887,April 2005,doi:10.1109/TGRS.2004.842441.)中提出了一种基于广义高斯模型的KI阈值选择(GKI，the modified KI criterion based on generalized Gaussianassumption)算法。该方法的实现过程为：通过自适应迭代滤波方式，对每次滤波后两个时相的SAR图像进行对数比操作得到差异图，然后用广义高斯分布的参数模型对差异图的变化类和未变化类像素的直方图分布进行建模，通过修改后的KI准则标准来自动确定最佳的滤波迭代次数对应的最优阈值，最终根据最优阈值实现差异图的阈值分割得到最终的变化检测结果。该方法虽然能够较好拟合差异图像直方图的分布，但是，该方法仍然存在的不足之处是，由于一些难分和易错分的像素都分布于阈值附近，因此使用单一的阈值对差异图进行分割，如果选择出的阈值不够准确，会导致检测结果的总错误数增多，降低了检测结果的准确率。
发明内容
本发明的目的在于针对上述现有技术的不足，提出一种基于阈值融合和邻域投票的SAR图像变化检测方法，用于解决现有技术存在的检测结果的总错误数增多、检测准确率低的问题。
实现本发明目的的具体思路是，首先对两幅待检测的SAR图像进行预处理，预处理后利用改进的差异图生成方法生成两时相SAR图像的差异图像；然后采用阈值融合策略将差异图像的GKI算法阈值和最大类间方差法(OTSU)算法阈值融合得到融合阈值，利用融合阈值对差异图像进行预分类，得到含有中间类，变化类、未变化类像素的预分类结果；最后对预分类结果图中的每个中间类像素结合该像素空间位置信息，通过邻域投票的方式进行细分，得到最终的变化检测结果图。
本发明的具体实现步骤如下：
(1)对待检测的两个不同时相的SAR图像进行预处理：
对两幅大小相同且时相不同的待检测SAR图像分别进行PPB滤波，对滤波后的图像中的每个像素值进行以10为底的对数操作后再进行10倍放大，得到两幅对数尺度的图像I1和I2；
(2)获得差异图：
利用下式，计算差异图中的每个像素值：
其中，D(m,n)表示差异图中第m行、第n列处的像素值，k表示图像I1或图像I2中与差异图中第m行、第n列对应位置像素为中心的像素邻域边长，∑表示求和操作，(x,y)表示图像I1位于第x行、第y列位置的像素或图像I2中与图像I1位于第x行、第y列对应位置的像素，∈表示属于符号，Ω表示图像I1或图像I2与差异图中第m行、第n列对应位置像素为中心的k×k邻域中所有像素的集合，I1(x,y)表示图像I1位于第x行、第y列位置的像素值，I2(x,y)表示与图像I1位于第x行、第y列对应位置的像素值，log(·)表示以10为底的对数操作；
(3)获得融合阈值：
(3a)分别利用GKI算法和OTSU算法，得到差异图的GKI阈值T1和OTSU阈值T2；
(3b)按照下式，对差异图的GKI阈值T1和OTSU阈值T2进行融合，得到融合阈值T：
T＝αT1+βT2
其中，α、β分别表示差异图的GKI阈值T1和OTSU阈值T2的权重，0≤α≤1，0≤β≤1，α+β＝1；
(4)获取预分类结果图：
将差异图中D(m,n)≤0.9T位置处的所有像素的像素值均设置为‘0’，用于表示未变化类的像素；将差异图中D(m,n)≥2.35T位置处的所有像素的像素值均设置为‘1’，用于表示变化类的像素；将差异图中0.9T＜D(m,n)＜2.35T位置处的所有像素的像素值均设置为‘0.5’，用于表示中间类的像素，得到预分类结果图；
(5)对中间类进行邻域投票分类：
(5a)以预分类结果图中每个中间类像素为中心，构建一个k1行、k2列的邻域图像块，3≤k1≤9，3≤k2≤9；
(5b)比较每个邻域图像块中变化类像素的总数e1、未变化类像素的总数f1、中间类的像素总数g1的大小，若e1为最大，则将中心的中间类像素更新为变化类；若f1最大，则将中心的中间类像素更新为未变化类；若g1最大，则中心的中间类像素仍为中间类；若e1＝f1≥g1，则将中心的中间类像素随机更新为变化类或者未变化类；若e1＝g1＞f1，则将中心的中间类像素随机更新为变化类或者中间类；若f1＝g1＞e1，则将中心的中间类像素随机更新为未变化类或者中间类，得到更新后的预分类结果图；
(6)判断预分类结果图与更新后的预分类结果图的中间类像素数量是否相等，若是，则得到初步邻域投票图后执行步骤(7)，否则，对更新后的预分类结果图再次执行步骤(5)，重新进行邻域投票分类；
(7)判断更新后预分类结果图的中间类像素数量是否为零，若是，则执行步骤(9)，否则，执行步骤(8)；
(8)对初步邻域投票图进行扩大邻域投票分类：
(8a)以初步邻域投票图中每个中间类像素为中心，构建一个l1行、l2列的邻域图像块，l1≥k1，l2≥k2，5≤l1≤13，5≤l2≤13；
(8b)比较每个邻域图像块中变化类像素的总数e2、未变化类像素的总数f2的大小，若e2≥f2，则将该邻域图像块中心的中间类像素更新为变化类；若e2＜f2，则将该邻域图像块中心的中间类像素更新为未变化类后执行步骤(9)；
(9)得到仅包含变化类和未变化类的分类结果图。
本发明与现有技术相比具有以下优点：
第一，由于本发明采用GKI和OTSU融合阈值对差异图像进行预分类，有效的将分布在阈值附近区分度不明显的变化类和未变化类的像素和难分类点选择了出来，克服了现有技术中的检测结果总错误数增多的问题，使得本发明具有针对性的提高难分点像素分类正确率的优点。
第二，由于本发明采用邻域投票的中间类分类方法，对含有中间类的初步分类结果进行再分类，引入像素的空间位置信息，根据初步分类结果中间类像素邻域内的三类像素对分类的投票贡献，实现中间类像素的再分类。克服了现有技术中对中间类像素的分类时没有考虑空间位置信息的问题，使得本发明具有提高变化检测的准确率的优点。
附图说明
图1为本发明的实现流程图；
图2为本发明仿真实验中的旧金山地区SAR图；
图3为本发明仿真实验中的Ottawa地区SAR图；
图4为本发明方法与两种现有技术对旧金山地区的仿真结果图；
图5为本发明方法与两种现有技术对Ottawa地区的仿真结果图。
具体实施方式
下面结合附图对本发明做进一步的详细描述。
结合附图1，对本发明的步骤做进一步的详细描述。
步骤1，对待检测的两个不同时相的SAR图像进行预处理。
对两幅大小相同且时相不同的待检测SAR图像分别进行PPB滤波，对滤波后的图像中的每个像素值进行以10为底的对数操作后再进行10倍放大，得到两幅对数尺度的图像I1和I2。
本步骤中的PPB滤波是引用C.Deledalle、L.Denis和F.Tupin在IEEETransactions on Image Processing上发表的"Iterative Weighted MaximumLikelihoodDenoising With Probabilistic Patch-Based Weights"中的方法(参考文献：C.Deledalle,L.Denis and F.Tupin,"Iterative Weighted Maximum LikelihoodDenoising With Probabilistic Patch-Based Weights,"in IEEE Transactions onImage Processing,vol.18,no.12,pp.2661-2672,Dec.2009,doi:10.1109/TIP.2009.2029593.)。
步骤2，获得差异图。
利用下式，计算差异图中的每个像素值：
其中，D(m,n)表示差异图中第m行、第n列处的像素值，k表示图像I1或图像I2中与差异图中第m行、第n列对应位置像素为中心的像素邻域边长，∑表示求和操作，(x,y)表示图像I1位于第x行、第y列位置的像素或图像I2中与图像I1位于第x行、第y列对应位置的像素，∈表示属于符号，Ω表示图像I1或图像I2与差异图中第m行、第n列对应位置像素为中心的k×k邻域中所有像素的集合，I1(x,y)表示图像I1位于第x行、第y列位置的像素值，I2(x,y)表示与图像I1位于第x行、第y列对应位置的像素值，log(·)表示以10为底的对数操作。
本步骤中获得差异图的方法参考了刘本强等人发表的论文“一种利用邻域相对熵的SAR影像变化检测方法”(参考文献：刘本强,赵争,魏钜杰.一种利用邻域相对熵的SAR影像变化检测方法[J].遥感信息,2018,033(003):91-97)中所提出的方法。
步骤3，获得融合阈值。
分别利用GKI算法和OTSU算法，得到差异图的GKI阈值T1和OTSU阈值T2；
所述的利用GKI算法得到差异图像GKI阈值T1的具体步骤如下：
第一步，对差异图中的每个像素进行归一化，得到归一化后差异图的直方图；
第二步，利用下式，分别计算未变化类像素分布的概率Pu(t)、未变化类像素分布的均值mu(t)、未变化类像素分布的方差变化类像素分布的概率Pc(t)、变化类像素分布的均值mc(t)、变化类像素分布的方差/>
其中，t表示归一化差异图的直方图中以间隔为1依次选取的阈值，t＝1,2···255，X表示归一化差异图的直方图中的像素灰度级，h(X)表示归一化差异图的直方图中X像素灰度级的像素占比；
第三步，按照下式，给定的广义高斯分布模型对小于直方图阈值的未变化类和大于直方图阈值的变化类像素的灰度级分布进行拟合，得到每一个阈值对应的变化类的类条件概率密度p(X|ωc)和未变化类像素的类条件概率密度p(X|ωu)：
其中，p(X|ωc)表示在变化类ωc条件下像素灰度级X的类条件概率密度值，ac(t)表示变化类模型的系数、e(·)表示以自然常数e为底的指数操作，|·|表示绝对值操作，bc(t)表示变化类模型的指数系数，βc(t)表示变化类模型的形状参数，p(X|ωu)表示在未变化类ωu条件下像素灰度级X的类条件概率密度值，au(t)表示未变化类模型的系数、bu(t)表示未变化类模型的指数系数，βu(t)表示未变化类模型的形状参数；
从0到100，每隔0.01取一个值分别给βc(t)和βu(t)赋值，按照下式，计算r(βc(t))和r(βu(t))的函数值，将所有赋值后的βc(t)、βu(t)以及与其对应的r(βc(t))和r(βu(t))的函数值组成查找表：
其中，Γ(·)表示伽玛函数；
按照下式，计算变化类的像素灰度级X绝对值的条件均值E[|X||ωc]和未变化类像素灰度级X绝对值的条件均值E[|X||ωu]：
其中，E[|X||ωc]表示变化类条件下像素灰度级X的绝对值的条件均值，E[|X||ωu]表示未变化类条件下像素灰度级X的绝对值的条件均值；
按照下式，计算变化类的归一化方差值ρc(t)和未变化类的归一化方差值ρu(t)：
在包含βc(t)、βu(t)、r(βc(t))和r(βu(t))的查找表中，找到使|ρc(t)-r(βc(t))|最小的βc(t)值，即为变化类模型的形状参数值βc(t)，找到使|ρu(t)-r(βu(t))|最小的βu(t)值，即为未变化类模型的形状参数值βu(t)；
按照下式，计算变化类的类条件概率密度函数的系数ac(t)、指数系数bc(t)和未变化类的类条件概率密度函数的系数au(t)、指数系数bu(t)：
第四步，按照下式，计算每个阈值对应的代价函数值J(t)：
其中，ln(·)表示以e为底的对数操作；
第五步，从所有阈值对应的代价函数值中找出代价函数值最小的阈值，作为GKI阈值T1。
GKI算法指的是Y.Bazi、L.Bruzzone和F.Melgani在IEEE Transactions onGeoscience and Remote Sensing上发表的论文"An unsupervised approach based onthe generalized Gaussian model to automatic change detection in multitemporalSAR images"中的方法。本步骤参考了论文方法的部分内容。(参考文献：Y.Bazi,L.Bruzzone and F.Melgani,"An unsupervised approach based on the generalizedGaussian model to automatic change detection in multitemporal SAR images,"inIEEE Transactions on Geoscience and Remote Sensing,vol.43,no.4,pp.874-887,April 2005,doi:10.1109/TGRS.2004.842441.)。
OTSU算法又称为最大类间方差法，是一种确定图像二值化分割阈值的算法，本实施例中，OTSU算法采用的是N.Otsu发表在IEEE Transactions on Systems,Man,andCybernetics的"A Threshold Selection Method from Gray-Level Histograms"(参考文献：N.Otsu,"A Threshold Selection Method from Gray-Level Histograms,"IEEETransactions on Systems,Man,and Cybernetics,vol.9,no.1,pp.62-66,1979)里的方法，按照OTSU算法求得的阈值T2进行图像二值化分割后，变化类与未变化类的类间方差最大。
按照下式，对差异图的GKI阈值T1和OTSU阈值T2进行融合，得到融合阈值T：
T＝αT1+βT2
其中，α、β分别表示差异图的GKI阈值T1和OTSU阈值T2权重，0≤α≤1，0≤β≤1，α+β＝1。
步骤4，获取预分类结果图。
将差异图中D(m,n)≤0.9T位置处的所有像素的像素值均设置为‘0’，用于表示未变化类的像素；将差异图中D(m,n)≥2.35T位置处的所有像素的像素值均设置为‘1’，用于表示变化类的像素；将差异图中0.9T＜D(m,n)＜2.35T位置处的所有像素的像素值均设置为‘0.5’，用于表示中间类的像素，得到预分类结果图；
本步骤中阈值的两个系数值0.9和2.35是在多次仿真实验中得出的经验值，适用于多组数据集。差异图像经过此步骤，将得到含有变化类、中间类、未变化类三类的初步分类结果，在本实施例中，中间类像素的分布位置大多在变化区域的边界位置，位于变化类与未变化类交界处的像素值既接近变化类，又接近未变化类，区分度不明显，是一些难分并且容易分错的点，通过预分类，将这些难分的像素点筛选并保留，通过后面步骤再进行分类。
步骤5，对中间类进行邻域投票分类。
以预分类结果图中每个中间类像素为中心，构建一个k1行、k2列的邻域图像块，3≤k1≤9，3≤k2≤9。
比较每个邻域图像块中变化类像素的总数e1、未变化类像素的总数f1、中间类的像素总数g1的大小，若e1为最大，则将中心的中间类像素更新为变化类；若f1最大，则将中心的中间类像素更新为未变化类；若g1最大，则中心的中间类像素仍为中间类；若e1＝f1≥g1，则将中心的中间类像素随机更新为变化类或者未变化类；若e1＝g1＞f1，则将中心的中间类像素随机更新为变化类或者中间类；若f1＝g1＞e1，则将中心的中间类像素随机更新为未变化类或者中间类，得到更新后的预分类结果图。
此步骤结合预分类结果中间类像素的邻域空间的预分类结果，将邻域内的每个类像素的数量看作是每个类的票数，根据少数服从多数的投票原则对中间类像素进行投票分类。这个步骤不仅利用了预分类的像素值信息，又结合了投票分类的空间位置预分类结果信息，这提高了中间类这些难分像素点的分类准确率。
步骤6，判断预分类结果图与更新后的预分类结果图的中间类像素数量是否相等，若是，则得到初步邻域投票图后执行步骤7，否则，对更新后的预分类结果图再次执行步骤5，重新进行邻域投票分类。
步骤7，判断更新后预分类结果图的中间类像素数量是否为零，若是，则执行步骤9，否则，执行步骤8。
如果更新后的预分类结果图的中间类像素数量为零，那么在步骤5中就已经将所有的中间类都进行了投票分类，得到的结果就是变化检测的最终结果。
步骤8，对初步邻域投票图进行扩大邻域投票分类。
以初步邻域投票图中每个中间类像素为中心，构建一个l1行、l2列的邻域图像块，l1≥k1，l2≥k2，5≤l1≤13，5≤l2≤13。
比较每个邻域图像块中变化类像素的总数e2、未变化类像素的总数f2的大小，若e2≥f2，则将该邻域图像块中心的中间类像素更新为变化类；若e2＜f2，则将该邻域图像块中心的中间类像素更新为未变化类后执行步骤9。
本步骤增加邻域大小，以中间类像素为中心，结合这些像素点的更大邻域信息对其进行分类，在更大邻域的投票过程中，不考虑中间类的票数，保证了所有像素都能完成分类。
步骤9，得到仅包含变化类和未变化类的二值结果图。
下面结合仿真实验对本发明的效果做进一步的说明。
1.仿真条件：
本发明仿真实验的运行环境是：Windows10，64位操作系统，本发明的仿真实验是在主频2.60GHZ的Intel Core(TM)i7-4720 CPU，内存8GB的硬件环境和MATLAB R2014b软件环境下完成的。
2.仿真内容：
本发明仿真实验是采用本发明和两种现有技术分别对旧金山地区和Ottawa地区的合成孔径雷达图像进行仿真，得到合成孔径雷达图像变化检测图，并对结果进行比对。
两种现有技术分别是：
PCANet方法是中国海洋大学2016年在IEEE Geoscience and Remote SensingLetters上发表的论文"Automatic Change Detection in Synthetic Aperture RadarImages Based on PCANet"中所述的方法和代码。参考文献：F.Gao,J.Dong,B.Li andQ.Xu,"Automatic Change Detection in Synthetic Aperture Radar Images Based onPCANet,"in IEEE Geoscience and Remote Sensing Letters,vol.13,no.12,pp.1792-1796,Dec.2016,doi:10.1109/LGRS.2016.2611001，参考代码网址：https://github.com/summitgao/SAR_Change_Detection_GarborPCANet。
GKI阈值法是对Yakoub Bazi等人在IEEE Transactions on Geoscience andRemote Sensing上发表的论文“An unsupervised approach based on the generalizedGaussian model to automatic change detection in multitemporal SAR images”中所述方法的复现，参考文献：Y.Bazi,L.Bruzzone and F.Melgani,"An unsupervisedapproach based on the generalized Gaussian model to automatic changedetection in multitemporal SAR images,"in IEEE Transactions on Geoscience andRemote Sensing,vol.43,no.4,pp.874-887,April 2005,doi:10.1109/TGRS.2004.842441。
本发明仿真实验所用数据为真实的合成孔径雷达数据集，第一组数据的第一时相图为2003年8月是通过ERS-2卫星传感器获取的美国旧金山地区的SAR图像，如图2(a)所示。第一组数据的第二时相图为2004年5月通过ERS-2卫星传感器获取的美国旧金山地区的SAR图像，如图2(b)所示。图2(c)为旧金山地区真实的变化参考图。图2中的所有图像尺寸均为256×256像素；第二组数据第一时相图为1997年5月加拿大国防研究与发展部提供的由RADARSAT SAR传感器获得的Ottawa地区的数据集，如图3(a)所示。第二组数据第二时相图为1997年8月加拿大国防研究与发展部提供的由RADARSAT SAR传感器获得的Ottawa地区的数据集，如图3(b)所示。图3(c)为Ottawa地区真实的变化参考图。图3中的所有图像尺寸为290×350像素。
图4(a)为采用现有技术PCANet方法对图2(a)、图2(b)旧金山地区合成孔径雷达图像进行仿真的结果图。图4(b)为采用现有技术GKI阈值法对图2(a)、图2(b)旧金山地区合成孔径雷达图像进行仿真的结果图。图4(c)为采用本发明对图2(a)、图2(b)旧金山地区合成孔径雷达图像进行仿真的结果图。
图5(a)为采用现有技术PCANet方法对图3(a)、图3(b)Ottawa地区合成孔径雷达图像进行仿真的结果图。图5(b)为采用现有技术GKI阈值法对图3(a)、图3(b)Ottawa地区合成孔径雷达图像进行仿真的结果图。图5(c)为采用本发明对图3(a)、图3(b)Ottawa地区合成孔径雷达图像进行仿真的结果图。
将图4中三种方法检测结果和图2(c)旧金山地区真实变化的参考图相互进行比较，可以直观的看出，图4(a)中变化区域的检测明显不完整，有大量的漏检像素。图4(b)中明显检测出大量的虚警像素，并且在变化区域边界位置存在大量的虚警，说明现有技术PCANet方法和GKI阈值法不能准确检测出旧金山地区合成孔径雷达图像的变化信息。相比之下，图4(c)中对于小块的变化区域以及边缘检测都较为准确，更接近于真实的变化参考图，说明采用本发明方法对旧金山地区合成孔径雷达图像中发生变化的区域的检测更为准确。
将图5中三种方法检测结果和图3(c)Ottawa地区真实变化的参考图相互比较，可以直观的看出，图5(a)中在变化区域的中间部分也存在着很多的漏检像素，同时对小块变化区域细节检测不好，导致小块变化区域形状保持不好，说明现有技术PCANet方法不能准确检测出Ottawa地区合成孔径雷达图像的变化信息；图5(b)中变化区域边缘位置检测出许多细小突出的虚警毛刺，检测的边缘不平滑，产生了很多虚警，说明现有技术GKI阈值法对Ottawa地区合成孔径雷达图像中变化区域的边缘检测不准确；相比之下，图5(c)中的细节保持更好，并且更接近于真实的变化参考图，说明采用本发明方法对Ottawa地区合成孔径雷达图像中变化区域的检测更为准确。
3.仿真效果分析：
在本发明与两种现有技术(PCANet方法、GKI阈值法)分别对旧金山地区和Ottawa地区的合成孔径雷达图像进行仿真的过程中，使用虚警像素个数FP、漏检像素个数FN、总错误数OE、总体分类精度PCC以及仿真结果图与真实变化参考图之间的Kappa系数五个评价指标对本发明的仿真实验结果进行评价。Kappa系数越大，仿真实验的结果图与真实变化参考图之间越接近，变化检测结果更好。
本发明仿真实验的三种方法对旧金山地区检测结果的评价指标值列于表1。其中，PCANet表示现有技术的PCANet方法，GKI表示现有技术的GKI阈值法。
表1对旧金山地区仿真实验检测结果指标一览表
算法 FP FN OE PCC(％) Kappa(％)
PCANet 125 883 1008 98.46 87.48
GKI 1206 118 1324 97.98 86.26
本发明 337 412 749 98.86 91.32
本发明仿真实验的三种方法对Ottawa地区检测结果的评价指标值列于表2。
表2对Ottawa地区仿真实验检测结果指标一览表
算法 FP FN OE PCC(％) Kappa(％)
PCANet 726 1112 1838 98.22 93.06
GKI 1572 676 2248 97.79 91.87
本发明 1070 658 1728 98.30 93.67
从表1中可以看出，处理旧金山地区合成孔径雷达图像时，本发明达到了91.32％的Kappa系数值，分别比两种现有技术(PCANet方法、GKI阈值法)的Kappa系数值提高了3.84％和5.06％。而且，本发明与两种现有技术相比，本发明总错误数更低，虚警像素和漏检像素数量更加平衡。
从表2中可以看出，处理Ottawa地区合成孔径雷达图像时，本发明达到了93.67％的Kappa系数值，分别比两种现有技术(PCANet方法、GKI阈值法)的Kappa系数值提高了0.61％和1.8％。同样能看出，本发明与两种现有技术相比，本发明总错误数更低，检测结果更好。
综上所述，相较于两种现有技术(PCANet方法、GKI阈值法)，本发明能够得到与参考图像更为接近、准确率更高的变化检测结果。

Claims (1)

1.一种基于阈值融合和邻域投票的SAR图像变化检测方法，其特征在于，采用GKI和OTSU的融合阈值对差异图像进行预分类，采用邻域投票的方式对预分类结果中的中间类进行细分，该方法的步骤包括如下：
(1)对待检测的两个不同时相的SAR图像进行预处理：
对两幅大小相同且时相不同的待检测SAR图像分别进行PPB滤波，对滤波后的图像中的每个像素值进行以10为底的对数操作后再进行10倍放大，得到两幅对数尺度的图像I1和I2；
(2)获得差异图：
利用下式，计算差异图中的每个像素值：
其中，D(m,n)表示差异图中第m行、第n列处的像素值，k表示图像I1或图像I2中与差异图中第m行、第n列对应位置像素为中心的像素邻域边长，∑表示求和操作，(x,y)表示图像I1位于第x行、第y列位置的像素或图像I2中与图像I1位于第x行、第y列对应位置的像素，∈表示属于符号，Ω表示图像I1或图像I2与差异图中第m行、第n列对应位置像素为中心的k×k邻域中所有像素的集合，I1(x,y)表示图像I1位于第x行、第y列位置的像素值，I2(x,y)表示与图像I1位于第x行、第y列对应位置的像素值，log(·)表示以10为底的对数操作；
(3)获得融合阈值：
(3a)分别利用GKI算法和OTSU算法，得到差异图的GKI阈值T1和OTSU阈值T2；
所述GKI算法的具体步骤如下：
第一步，对差异图中的每个像素进行归一化，得到归一化后差异图的直方图；
第二步，利用下式，分别计算未变化类像素分布的概率Pu(t)、未变化类像素分布的均值mu(t)、未变化类像素分布的方差变化类像素分布的概率Pc(t)、变化类像素分布的均值mc(t)、变化类像素分布的方差/>
其中，t表示归一化差异图的直方图中以间隔为1依次选取的阈值，t＝1,2···255，X表示归一化差异图的直方图中的像素灰度级，h(X)表示归一化差异图的直方图中X像素灰度级的像素占比；
第三步，按照下式，给定的广义高斯分布模型对小于直方图阈值的未变化类和大于直方图阈值的变化类像素的灰度级分布进行拟合，得到每一个阈值对应的变化类的类条件概率密度p(X|ωc)和未变化类像素的类条件概率密度p(X|ωu)：
其中，X|ωc表示在变化类ωc条件下像素灰度级X的类条件概率密度值，ac(t)表示变化类模型的系数、e(·)表示以自然常数e为底的指数操作，|·|表示绝对值操作，bc(t)表示变化类模型的指数系数，βc(t)表示变化类模型的形状参数X|ωu表示在未变化类ωu条件下像素灰度级X的类条件概率密度值，au(t)表示未变化类模型的系数、bu(t)表示未变化类模型的指数系数，βu(t)表示未变化类模型的形状参数；
从0到100，每隔0.01取一个值分别给βc(t)和βu(t)赋值，按照下式，计算r(βc(t))和r(βu(t))的函数值，将所有赋值后的βc(t)、βu(t)以及与其对应的r(βc(t))和r(βu(t))的函数值组成查找表：
其中，Γ(·)表示伽玛函数；
按照下式，计算变化类的像素灰度级X绝对值的条件均值E[|X||ωc]和未变化类像素灰度级X绝对值的条件均值E[|X||ωu]：
其中，|X||ωc表示变化类条件下像素灰度级X的绝对值的条件均值，|X||ωu表示未变化类条件下像素灰度级X的绝对值的条件均值；
按照下式，计算变化类的归一化方差值ρc(t)和未变化类的归一化方差值ρu(t)：
在包含βc(t)、βu(t)、r(βc(t))和r(βu(t))的查找表中找到使|ρc(t)-r(βc(t))|最小的βc(t)值，变化类模型的形状参数值βc(t)；
按照下式，计算变化类的类条件概率密度函数的系数ac(t)、指数系数bc(t)和未变化类的类条件概率密度函数的系数au(t)、指数系数bu(t)：
第四步，按照下式，计算每个阈值对应的代价函数值J(t)：
其中，ln(·)表示以e为底的对数操作；
第五步，从所有阈值对应的代价函数值中找出代价函数值最小的阈值，作为GKI阈值T1；
(3b)按照下式，对差异图的GKI阈值T1和OTSU阈值T2进行融合，得到融合阈值T：
T＝αT1+βT2
其中，α、β分别表示差异图的GKI阈值T1和OTSU阈值T2的权重，0≤α≤1，0≤β≤1，α+β＝1；
(4)获取预分类结果图：
将差异图中D(m,n)≤0.9T位置处的所有像素的像素值均设置为‘0’，用于表示未变化类的像素；将差异图中D(m,n)≥2.35T位置处的所有像素的像素值均设置为‘1’，用于表示变化类的像素；将差异图中0.9T＜D(m,n)＜2.35T位置处的所有像素的像素值均设置为‘0.5’，用于表示中间类的像素，得到预分类结果图；
(5)对中间类进行邻域投票分类：
(5a)以预分类结果图中每个中间类像素为中心，构建一个k1行、k2列的邻域图像块，3≤k1≤9，3≤k2≤9；
(5b)比较每个邻域图像块中变化类像素的总数e1、未变化类像素的总数f1、中间类的像素总数g1的大小，若e1为最大，则将中心的中间类像素更新为变化类；若f1最大，则将中心的中间类像素更新为未变化类；若g1最大，则中心的中间类像素仍为中间类；若e1＝f1≥g1，则将中心的中间类像素随机更新为变化类或者未变化类；若e1＝g1＞f1，则将中心的中间类像素随机更新为变化类或者中间类；若f1＝g1＞e1，则将中心的中间类像素随机更新为未变化类或者中间类，得到更新后的预分类结果图；
(6)判断预分类结果图与更新后的预分类结果图的中间类像素数量是否相等，若是，则得到初步邻域投票图后执行步骤(7)，否则，对更新后的预分类结果图再次执行步骤(5)，重新进行邻域投票分类；
(7)判断更新后预分类结果图的中间类像素数量是否为零，若是，则执行步骤(9)，否则，执行步骤(8)；
(8)对初步邻域投票图进行扩大邻域投票分类：
(8a)以初步邻域投票图中每个中间类像素为中心，构建一个l1行、l2列的邻域图像块，l1≥k1，l2≥k2，5≤l1≤13，5≤l2≤13；
(8b)比较每个邻域图像块中变化类像素的总数e2、未变化类像素的总数f2的大小，若e2≥f2，则将该邻域图像块中心的中间类像素更新为变化类；
若e2＜f2，则将该邻域图像块中心的中间类像素更新为未变化类后执行步骤(9)；
(9)得到仅包含变化类和未变化类的分类结果图。
CN202011160475.9A 2020-10-27 2020-10-27 基于阈值融合和邻域投票的sar图像变化检测方法 Active CN112348750B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011160475.9A CN112348750B (zh) 2020-10-27 2020-10-27 基于阈值融合和邻域投票的sar图像变化检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011160475.9A CN112348750B (zh) 2020-10-27 2020-10-27 基于阈值融合和邻域投票的sar图像变化检测方法

Publications (2)

Publication Number Publication Date
CN112348750A CN112348750A (zh) 2021-02-09
CN112348750B true CN112348750B (zh) 2023-08-18

Family

ID=74358653

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011160475.9A Active CN112348750B (zh) 2020-10-27 2020-10-27 基于阈值融合和邻域投票的sar图像变化检测方法

Country Status (1)

Country Link
CN (1) CN112348750B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114519696B (zh) * 2021-12-31 2022-11-29 扬州盛强薄膜材料有限公司 基于光学智能化的pvc热收缩膜检测方法及系统
CN114754353B (zh) * 2022-04-13 2023-05-30 山西大学 融合邻域粗糙集机器学习的循环流化床锅炉燃烧优化方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2006113582A2 (en) * 2005-04-15 2006-10-26 Mississippi State University Change analyst
CN103020978A (zh) * 2012-12-14 2013-04-03 西安电子科技大学 结合多阈值分割与模糊聚类的sar图像变化检测方法
CN103839257A (zh) * 2013-12-24 2014-06-04 西安电子科技大学 一种广义高斯k&i的sar图像变化检测方法
CN104200472A (zh) * 2014-08-30 2014-12-10 西安电子科技大学 基于非局部小波信息的遥感图像变化检测方法
CN106296655A (zh) * 2016-07-27 2017-01-04 西安电子科技大学 基于自适应权值和高频阈值的sar图像变化检测方法
CN108764119A (zh) * 2018-05-24 2018-11-06 西安电子科技大学 基于迭代最大类间方差的sar图像变化检测方法
WO2018213752A1 (en) * 2017-05-19 2018-11-22 Fundación Deusto Method and system for monitoring and evaluation of pressure ulcer severity
CN110555841A (zh) * 2019-09-10 2019-12-10 西安电子科技大学 基于自注意图像融合和dec的sar图像变化检测方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2948226C (en) * 2014-06-30 2023-09-05 Ventana Medical Systems, Inc. Detecting edges of a nucleus using image analysis

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2006113582A2 (en) * 2005-04-15 2006-10-26 Mississippi State University Change analyst
CN103020978A (zh) * 2012-12-14 2013-04-03 西安电子科技大学 结合多阈值分割与模糊聚类的sar图像变化检测方法
CN103839257A (zh) * 2013-12-24 2014-06-04 西安电子科技大学 一种广义高斯k&i的sar图像变化检测方法
CN104200472A (zh) * 2014-08-30 2014-12-10 西安电子科技大学 基于非局部小波信息的遥感图像变化检测方法
CN106296655A (zh) * 2016-07-27 2017-01-04 西安电子科技大学 基于自适应权值和高频阈值的sar图像变化检测方法
WO2018213752A1 (en) * 2017-05-19 2018-11-22 Fundación Deusto Method and system for monitoring and evaluation of pressure ulcer severity
CN108764119A (zh) * 2018-05-24 2018-11-06 西安电子科技大学 基于迭代最大类间方差的sar图像变化检测方法
CN110555841A (zh) * 2019-09-10 2019-12-10 西安电子科技大学 基于自注意图像融合和dec的sar图像变化检测方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于差异图和块分类的雷达图像变化检测;王平;;长春师范大学学报(第08期);全文 *

Also Published As

Publication number Publication date
CN112348750A (zh) 2021-02-09

Similar Documents

Publication Publication Date Title
Li et al. Automatic pixel‐level multiple damage detection of concrete structure using fully convolutional network
Li et al. Automatic pavement-crack detection and segmentation based on steerable matched filtering and an active contour model
Guo et al. Mining parameter information for building extraction and change detection with very high-resolution imagery and GIS data
Chen et al. An automated approach for updating land cover maps based on integrated change detection and classification methods
CN109871875B (zh) 一种基于深度学习的建筑物变化检测方法
CN111860236B (zh) 一种基于迁移学习的小样本遥感目标检测方法及系统
CN111161229B (zh) 一种基于几何主动轮廓模型和稀疏自编码的变化检测方法
CN105608698A (zh) 一种基于sae的遥感图像变化检测方法
CN112348750B (zh) 基于阈值融合和邻域投票的sar图像变化检测方法
CN111274926B (zh) 图像数据筛选方法、装置、计算机设备和存储介质
CN108447055A (zh) 基于spl和ccn的sar图像变化检测方法
Peng et al. Object-based change detection from satellite imagery by segmentation optimization and multi-features fusion
CN112819821B (zh) 一种细胞核图像检测方法
CN108492298A (zh) 基于生成对抗网络的多光谱图像变化检测方法
Yang et al. Classified road detection from satellite images based on perceptual organization
CN108171119B (zh) 基于残差网络的sar图像变化检测方法
CN104899892A (zh) 一种快速的星图图像星点提取方法
CN103971362B (zh) 基于直方图和精英遗传聚类算法的sar图像变化检测
CN113378642B (zh) 一种对农村违法占地建筑物进行检测的方法
CN101839980A (zh) 一种基于分割窗的无监督遥感图像变化检测方法
Oga et al. River state classification combining patch-based processing and CNN
CN113988222A (zh) 一种基于Faster-RCNN的森林火灾检测与识别方法
CN106951924B (zh) 基于AdaBoost算法的地震相干体图像断层自动识别方法及系统
Saboori et al. Combining multi-scale textural features from the panchromatic bands of high spatial resolution images with ANN and MLC classification algorithms to extract urban land uses
CN114419465B (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