CN102360503B - 基于空间贴近度和像素相似性的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于空间贴近度和像素相似性的SAR图像变化检测方法，属于SAR图像变化检测领域。其实现过程为：首先，对同一地区不同时间的两幅SAR图像进行几何校正和配准；其次，运用双边相似性构造差异影像图，最后，利用阈值方法对差异影像图进行阈值处理，得到变化检测结果图。本发明通过空间贴近度和像素相似性的差异影像图构造方法解决了SAR图像变化检测的差异影像图构造受SAR图像斑点噪声影响严重的问题，从而提高了变化检测的精度，本发明可用于环境监控，农业研究，城市区域研究，森林监控等领域。

Description

基于空间贴近度和像素相似性的SAR图像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及SAR图像变化检测，可以用来解决SAR图像变化检测领域中检测精度不高的问题，适用对多种复杂SAR图像提高SAR图像变化检测的精度。

背景技术

合成孔径雷达SAR是现代遥感领域的一项突破性技术，成为获取信息的一个重要手段。随着SAR技术的不断发展和成熟，及SAR系统的分辨率的不断提高，使得SAR系统具有了全天候、全天时及覆盖面积大等优点。不同时间同一场景的SAR图像的获得，使得通过SAR图像变化检测技术提供目标与场景实时动态信息成为可能。

SAR图像变化检测是遥感技术的主要应用之一。通过对不同时间同一场景的各种图像进行对比分析，根据图像之间的差异来得到人们所需要的地物或目标的变化信息。变化检测技术可以检测出不同时期图像灰度值或局部纹理之间的变化，在此基础上获得需要的目标在形状、位置、数量、及其它属性的变化情况。

然而，SAR图像在具有上述优点的同时，斑点噪声影响是其应用的一个主要瓶颈，如何从SAR图像自身特点出发，有效抑制斑点噪声是SAR图像实际应用中必须考虑的一个问题。

随着计算机和遥感成像技术的不断发展，SAR图像变化检测技术也在不断提高和改进中，在各国学者的努力下提出了以下几种经典的变化检测算法：

图像差值法：其主要过程是将两幅不同时间相同地域的SAR图像对应像素的灰度值相减得到差异影像图。理论上，在得到的差值图像上，差值为0或接近0的认为是不变区域，不为0的认为是变化区域。这种方法的优点在于理论相对简单，容易理解和掌握。但是缺点也比较明显，由于该方法主要通过分析地物光谱值改变的灰度差值图像来实现变化检测，但在某些情况下，仅仅利用光谱特征的差值图像难以实现地物的变化检测，易产生“伪变化”信息，这样会使得检测结果的错误率表较高。

图像比值法：其主要过程是计算两幅不同时间相同地域的SAR图像对应像素灰度值的比值以得到差异影像图，其中像素比值为1或近似于1的认为是未发生变化的区域，反之，远大于或远小于1的认为是发生变化的区域。该方法虽然可以克服乘性噪声的干扰，但该方法没有考虑邻域信息，使得检测错误率高。

分类后比较法：其主要过程是采用相同的分类体系分别对两幅不同时间相同地域的SAR图像进行独立的分类，然后比较分类结果，从中提取变化信息。如果对应像素的分类类别相同，则认为该像素没有发生变化，反之，发生了变化。该方法虽然可以回避所用多时相遥感图像数据应获取环境条件不同和传感器不同所带来的辐射归一化问题，并可直接获取变化的类型，数量，位置。但却存在着工作量很大，精度低的缺陷。

上面几种方法是较为常用的变化检测方法，另外还有变化向量分析法，主分量分析法等等。但是这些方法均受SAR图像斑点噪声影响较大，必须预先对两时相SAR图像进行滤波去噪，这样就有可能对原始图像数据进行改变，另外这些方法构造的差异影像图中变化类与非变化类区别都不够明显，使分割后检测精度不高。

发明内容

本发明的目的在于针对上述已有技术的不足，提出了一种基于空间贴近度和像素相似性的SAR图像变化检测方法，以克服现有方法的检测错误率高检测精度低的问题。

本发明的技术方案是：基于空间贴近度和像素相似性的原理，结合了比值方法与均值比值方法的优点，提出了一种双边相似性构造差异影像图的方法，根据双边相似性构造同一区域不同时间的两幅SAR图像的差异影像图DI，再利用KI阈值得到变化检测结果图。其具体实现步骤如下：

(1)对同一区域不同时间的两幅SAR图像A和B进行几何校正和配准；

(2)对校正和配准后的两幅SAR图像A和B构造差异影像图DI：

2a)取SAR图像A在位置(i，j)上的像素灰度值邻域集合y′₁(i，j)和SAR图像B在位置(i，j)上的像素灰度值邻域集合y′₂(i，j)：

$y_1^{\′}(i,j)=y_2^{\′}(i,j)=\underset{i=0}{\overset{N}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}{w}_d(i,j)\·w_s(i,j)\·X(i,j)$

其中：N表示邻域集合大小，N∈{3，5，7，9}，(i，j)表示SAR图像横轴为i、纵轴为j的像素点，X(i，j)表示在位置(i，j)上的像素值；

w_d(x，y)表示空间贴近度， $w_d(x,y)=\exp (-{|{\left(\frac{x}{\mathrm{sx}}\right)}^2+{\left(\frac{y}{\mathrm{sy}}\right)}^2|}^{\frac{1}{2}})$

x，y分别表示在邻域集合大小为N×N中，以中心像素点为坐标原点时，对应的像素点的横轴与纵轴的坐标值，sx，sy分别表示横轴和纵轴邻域集合大小的一半；

w_s(i，j)表示像素的相似性， $w_s(i,j)=\exp (-\frac{{|X(i,j)-X(x,y)|}^2}{{2\mathrm{\σ}}^2})$

X(x，y)表示邻域内中心点的像素值，X(i，j)表示邻域内位置(i，j)上的像素值；σ表示SAR图像A和B的标准差的和；

y′₂(i，j)的相似性，得到差异影像图DI在位置(i，j)上的相似度值DI(i，j)：

$\mathrm{DI}(i,j)=1-\min \begin{array}{}\end{array}\{\frac{y_1^{\′}(i,j)}{y_2^{\′}(i,j)},\frac{y_2^{\′}(i,j)}{y_1^{\′}(i,j)}\}$

其中，集合y′₁(i，j)与y′₂(i，j)越相似，则DI(i，j)值越小，即DI(i，j)的值越接近于0，SAR图像A和B在位置(i，j)上的像素点属于非变化类的可能性越大，反之，属于变化类的可能性越大；

2c)对相似度值乘以255，得到位置(i，j)的像素灰度值；

2d)对SAR图像A和B的每一位置(i，j)从左到右，从上到下重复步骤2a)-2c)，得到差异影像图DI；

(3)对差异影像图DI进行阈值处理，得到变化检测结果图。

本发明与现有技术相比具有如下特点：

1、本发明采用空间贴近度和像素相似性的SAR图像变化检测方法，抑制了SAR图像的部分斑点噪声，不需要对原始SAR图像滤波去噪，很好的保持原始图像的信息；同时有效地减少了过渡带的信息，而且将差异影像图的直方图分布调整为双峰分布，有利于对差异影像图的分割，提高了变化检测的精度；

2、仿真结果表明，本发明采用的双边相似性构造差异影像图方法较比值方法，对数比值方法，均值比值方法构造差异影像图方法的检测精度更高。

附图说明

图1是本发明方法的流程图；

图2是本发明仿真使用的Ottawa两幅SAR影像图和标准图；

图3是本发明方法和现有三种差异影像图构造方法的对比试验结果图；

图4是本发明仿真使用的Berma两幅SAR影像图和标准图；

图5是本发明方法和现有三种差异影像图构造方法的对比试验结果图。

具体实施方式

参照图1，本发明包括如下步骤：

步骤一：对同一区域不同时间的两幅SAR图像A和B进行几何校正和配准，最常应用的校正和配准方法是控制点位法和灰度插值的最近邻法，参见胡东红“两种图像校正算法在实际应用中的比较”2009，计算机工程与应用。

实际中，已有经过几何校正和配准好的图像，本发明只对几何校正和配准后的图像进行仿真实验。

步骤二：对校正和配准后的两幅SAR图像A和B构造差异影像图DI：

2a)取SAR图像A在位置(i，j)上的像素灰度值邻域集合y′₁(i，j)和SAR图像B在位置(i，j)上的像素灰度值邻域集合y′₂(i，j)：

$y_1^{\′}(i,j)=y_2^{\′}(i,j)=\underset{i=0}{\overset{N}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}{w}_d(i,j)\·w_s(i,j)\·X(i,j)$

其中：N表示邻域集合大小，N∈{3，5，7，9}，(i，j)表示SAR图像横轴为i、纵轴为j的像素点，X(i，j)表示在位置(i，j)上的像素值；

w_d(x，y)表示空间贴近度， $w_d(x,y)=\exp (-{|{\left(\frac{x}{\mathrm{sx}}\right)}^2+{\left(\frac{y}{\mathrm{sy}}\right)}^2|}^{\frac{1}{2}})$

x，y分别表示在邻域集合大小为N×N中，以中心像素点为坐标原点时，对应的像素点的横轴与纵轴的坐标值，sx，sy分别表示横轴和纵轴邻域集合大小的一半；

w_s(i，j)表示像素的相似性， $w_s(i,j)=\exp (-\frac{{|X(i,j)-X(x,y)|}^2}{{2\mathrm{\σ}}^2})$

X(x，y)表示邻域内中心点的像素值，X(i，j)表示邻域内位置(i，j)上的像素值；σ表示SAR图像A和B的标准差的和；

2b)比较两个邻域集合y′₁(i，j)和y′₂(i，j)的相似性，得到差异影像图DI在位置(i，j)上的相似度值DI(i，j)：

$\mathrm{DI}(i,j)=1-\min \begin{array}{}\end{array}\{\frac{y_1^{\′}(i,j)}{y_2^{\′}(i,j)},\frac{y_2^{\′}(i,j)}{y_1^{\′}(i,j)}\}$

其中，集合y′₁(i，j)与y′₂(i，j)越相似，则DI(i，j)值越小，即DI(i，j)的值越接近于0，SAR图像A和B在位置(i，j)上的像素点属于非变化类的可能性越大，反之，属于变化类的可能性越大；

2c)对相似度值DI(i，j)乘以255，得到位置(i，j)的像素灰度值；

2d)对SAR图像A和B的每一位置(i，j)从左到右，从上到下重复步骤2a)-2c)，得到差异影像图DI；

步骤三：对差异影像图DI进行阈值处理，得到变化检测结果图。

在SAR图像变化检测中，常用的阈值处理方法有KI阈值法和Ostu阈值法：

A)采用KI阈值方法进行，其步骤如下：

3a)设差异影像图DI的灰度值为0到L-1级，差异影像图上的变化类w_c和非变化类w_u，灰度值X_l的像素个数为h(X_l)；

3b)计算非变化类w_u的总概率P_u(T)、均值m_u(T)和方差

如下：

$P_u\left(T\right)=\underset{X_l=0}{\overset{T}{\mathrm{\Σ}}}h\left(X_l\right),$ $m_u\left(T\right)=\frac{1}{P_u\left(T\right)}\underset{X_l}{\overset{T}{\mathrm{\Σ}}}{X}_{l}h\left(X_l\right),$ ${\mathrm{\σ}}_u^2\left(T\right)=\frac{1}{P_u\left(T\right)}\underset{X_l}{\overset{T}{\mathrm{\Σ}}}{(X_l-m_u\left(T\right))}^{2}h\left(X_l\right)$

3c)计算变化类w_c的总概率P_c(T)、均值m_c(T)和方差

如下：

P_c(T)＝1-P_u(T)， $m_c\left(T\right)=\frac{1}{P_c\left(T\right)}\underset{X_l=T+1}{\overset{L-1}{\mathrm{\Σ}}}{X}_{l}h\left(X_l\right),$ ${\mathrm{\σ}}_c^2\left(T\right)=\frac{1}{P_c\left(T\right)}\underset{X_l=T+1}{\overset{L-1}{\mathrm{\Σ}}}{(X_l-m_c\left(T\right))}^{2}h\left(X_l\right)$

3d)计算两类间的准则函数如下：

J_KI(T)＝1+2[P₁(T)lnσ₁(T)+P₂(T)lnσ₂(T)]-2[P₁(T)lnP₁(T)+P₂(T)lnP₂(T)]

从1～L-1之间改变T，当J_KI(T)为最小值时的T值为最佳阈值TMax，对差异影像图DI按照最佳阈值进行分割，小于TMax的像素设为0，大于TMax的像素设为255，得到变化检测结果图。

B)采用Ostu阈值方法进行，其步骤如下：

3a)设差异影像图DI的灰度值为0到L-1级，像素灰度值i的像素个数为n_i，此时得到的像素总数为各值概率为

然后用阈值T将其分成两组，其中第一组C₀＝{0～T-1}，第二组C₁＝{T～L-1}；

3b)计算第一组C₀产生的概率w₀和平均值μ₀：

$w_0=\underset{i=0}{\overset{T-1}{\mathrm{\Σ}}}{p}_i,$ ${\mathrm{\μ}}_0=\underset{i=0}{\overset{T-1}{\mathrm{\Σ}}}\frac{{\mathrm{ip}}_i}{w_0},$

3c)计算第二组C₁产生的概率w₁和平均值μ₁：

$w_1=\underset{i=T}{\overset{L-1}{\mathrm{\Σ}}}{p}_i$ ${\mathrm{\μ}}_1=\underset{i=T}{\overset{L-1}{\mathrm{\Σ}}}\frac{{\mathrm{ip}}_i}{w_1}$

3d)计算第一组和第二组间的方差如下式：

δ²(T)＝w₀(μ₀-μ)²+w₁(μ₁-μ)²

其中，

是差异影像图DI的灰度平均值；

从1～L-1之间改变T，取δ²(T)为最大值时的T值为最佳阈值TMax，对差异影像图DI按照最佳阈值进行分割，小于TMax的像素设为0，大于TMax的像素设为255，得到变化检测结果图。

本发明的效果可以通过以下仿真进一步说明：

1、仿真参数

对于具有参考图的两组实验仿真图进行结果分析，其定量的评判准则参数漏检个数FN、错检个数FP和总检测错误数OE如下：

①计算漏检个数：统计实验结果图中发生变化区域的像素个数，与参考图中变化区域的像素个数进行对比，把参考图中发生变化但实验结果图中检测为未变化的像素个数，称为漏检个数FN；

②计算错检个数：统计实验结果图中未发生变化区域的像素个数，与参考图中未变化区域的像素个数进行对比，把参考图中未发生变化但实验结果图中检测为变化的像素个数，称为错检个数FP；

③总检测错误数OE等于漏检个数和错检个数的和。

2、仿真内容

①将Ottawa图像数据用双边相似性方法构造的差异影像图，以及比值方法、对数比值方法、均值比值方法构造的差异影像图，运用KI阈值得到结果图后，用定性、定量的方法对实验结果进行了分析；

②将Berma图像数据用双边相似性方法构造的差异影像图，以及比值方法、对数比值方法、均值比值方法构造的差异影像图，运用KI阈值得到结果图后，用定性、定量的方法对实验结果进行了分析。

3、仿真实验结果及分析

①带参考图的Ottawa区域SAR图像数据

应用比值法、对数比值法、均值比值法和双边相似性差异影像图构造方法在真实Ottawa区域SAR图像数据上。实验原图及参考图如图2所示。其中图2(a)表示1997年5月该地区地貌信息，图2(b)表示1997年8月该地区地貌信息，图2(c)表示变化检测的标准结果图。

比值法、对数比值法、均值比值法和双边相似性方法构造的差异影像图的变化检测结果如图3所示，检测精度的比较结果列于表1中。其中图3(a)表示比值方法的结果图，图3(b)表示对数比值方法的结果图，图3(c)表示均值比值方法的结果图，图3(d)表示双边相似性方法结果图。

表1  不同差异影像图构造方法的变化检测结果分析

从变化检测结果图中可以看出，本发明方法对于杂点的处理效果比较明显。这正是由于本发明方法不仅考虑了邻域空间信息对象素的影响，还注意到了邻域中像素的影响不能与中心像素等同看待，另外，利用像素相似性对变化类和非变化类像素进行更新，加大了两类的区分度。虽然，均值比值也考虑了邻域信息，但是，它将邻域像素与中心像素等同对待，这样会使未发生变化的像素被误认为变化像素；而比值法和对数比值法直接对象素作比值，忽略了邻域信息，差异图受斑点噪声的影响会使变化像素被漏检。

从表1中可以看出，本发明方法在变化检测总错误数上，比其它对比实验方法都要少，而且在误检数与漏检数上相对都比较均衡。同时从参考图3中可看出，本发明在视觉效果上也较好，其它方法要么噪声较大要么轮廓不清晰，而本发明方法所得的结果更加接近于参考图。总体来看，通过表1性能结果和参考图3的视觉效果说明了本发明方法对变化类与非变化类的精确检测，从而保证了阈值准确的确定，提高了最终变化检测的精确度，使得总错误数更少，误检数与漏检数上相对都比较均衡。

②带参考图的Berma区域SAR图像数据

应用比值法、对数比值法、均值比值法和双边相似性差异影像图构造方法在真实Berma区域SAR图像数据上。实验原图及参考图如图4所示。其中图4(a)表示1999年4月该地区地貌信息，图4(b)表示1999年5月该地区地貌信息，图4(c)表示变化检测的标准结果图。

比值法、对数比值法、均值比值法和双边相似性方法构造的差异影像图的变化检测结果如图5所示，检测精度的比较结果列于表2中。其中图5(a)表示比值方法的结果图，图5(b)表示对数比值方法的结果图，图5(c)表示均值比值方法的结果图，图5(d)表示双边相似性方法结果图。

表2  不同差异影像图构造方法的变化检测结果分析

从表2所列检测精度统计结果可以看出，对Berma数据的实验中，双边相似性方法的错检个数和漏检个数的都低于其它三种方法，总的错误数甚至低于对数比值方法4倍；通过检测Ottawa数据和Berma数据，可以看出，双边相似性方法效果优于现有的方法，这使得SAR图像变化检测的精度更加准确。也可以看出本发明方法有着很好的稳定性和适用性。

Claims (3)

1.一种基于空间贴近度和像素相似性的SAR图像变化检测方法，包括如下步骤： 

(1)对同一区域不同时间的两幅SAR图像A和B进行几何校正和配准； 

(2)对校正和配准后的两幅SAR图像A和B构造差异影像图DI： 

2a）取SAR图像A在位置(i,j)上的像素灰度值邻域集合y′₁(i,j)和SAR图像B在位置(i,j)上的像素灰度值邻域集合y′₂(i,j)： 

其中：N表示邻域集合大小，N∈{3,5,7,9}，(i,j)表示SAR图像横轴为i、纵轴为j的像素点，X(i,j)表示在位置(i,j)上的像素值； 

w_d(x,y)表示空间贴近度，

x，y分别表示在邻域集合大小为N×N中，以中心像素点为坐标原点时，对应的像素点的横轴与纵轴的坐标值，sx，sy分别表示横轴和纵轴邻域集合大小的一半； 

w_s(i,j)表示像素的相似性，

X_c表示邻域内中心点的像素值，X(i,j)表示邻域内位置(i,j)上的像素值；σ表示SAR图像A和B的标准差的和； 

2b）比较两个邻域集合y₁(i,j)和y′₂(i,j)的相似性，得到差异影像图DI在位置(i,j)上的相似度值DI(i,j)： 

其中，集合y₁(i,j)与y′₂(i,j)越相似，则DI(i,j)值越小，即DI(i,j)的值越接近于0，SAR图像A和B在位置(i,j)上的像素点属于非变化类的可能性越大，反之， 属于变化类的可能性越大； 

2c）对相似度值乘以255，得到位置(i,j)的像素灰度值； 

2d）对SAR图像A和B的每一位置(i,j)从左到右，从上到下重复步骤2a）-2c），得到差异影像图DI； 

(3)对差异影像图DI进行阈值处理，得到变化检测结果图。 

2.根据权利要求1所述的SAR图像变化检测方法，其中步骤(3)所述的对差异影像图DI进行阈值处理，采用KI阈值方法进行，其步骤如下： 

3a）设差异影像图DI的像素灰度值为0到L-1级，差异影像图上的变化类w_c和非变化类w_u，像素灰度值X_l的像素个数为

此时得到的像素总数为各像素灰度值概率为

3b）计算非变化类w_u的总概率P_u(T)、均值m_u(T)和方差

如下： 

3c）计算变化类w_c的总概率P_c(T)、均值m_c(T)和方差

如下： 

P_c(T)＝1-P_u(T)，

3d）计算两类间的准则函数如下： 

J_KI(T)＝1+2[P₁(T)lnσ₁(T)+P₂(T)lnσ₂(T)]-2[P₁(T)lnP₁(T)+P₂(T)lnP₂(T)] 

从1～L-1之间改变T，当J_KI(T)为最小值时的T值为最佳阈值TMax，对差异影像图DI按照最佳阈值进行分割，小于TMax的像素设为0，大于TMax的像素设为255，得到变化检测结果图。 

3.根据权利要求1所述的SAR图像变化检测方法，其中步骤(3)所述的对差异影像图DI进行阈值处理，进一步采用Ostu阈值方法进行，其步骤如下： 

3.1）设差异影像图DI的灰度值为0到L-1级，像素灰度值i的像素个数为n_i，此时得到的像素总数为

各值概率为

然后用阈值T将其分成两组，其中第一组C₀＝{0~T-1}，第二组C₁＝{T~L-1}； 

3.2）计算第一组C₀产生的概率w₀和平均值μ₀： 

3.3）计算第二组C₁产生的概率w₁和平均值μ₁： 

3.4）计算第一组和第二组间的方差如下式： 

δ²(T)＝w₀(μ₀-μ)²+w₁(μ₁-μ)²

其中，

是差异影像图DI的灰度平均值； 

从1～L-1之间改变T，取δ²(T)为最大值时的T值为最佳阈值TMax，对差异影像图DI按照最佳阈值进行分割，小于TMax的像素设为0，大于TMax的像素设为255，得到变化检测结果图。 

Images