CN103871039B - 一种sar图像变化检测差异图生成方法 - Google Patents

Abstract

本发明公开了一种的SAR图像变化检测差异图生成方法，主要解决现有SAR图像变化检测技术中差异图生成的问题。其实现步骤为：(1)输入两幅不同时相的SAR图像；(2)对两幅图全图逐像素遍历，确定待恢复像素点，初始化待恢复像素点的搜索窗口和相似性测量图块；(3)计算每个待恢复像素点与其搜索窗口中其他像素点的相似性测量图块之间对应位置像素点灰度比值的概率密度分布函数；(4)计算每个待恢复像素点与其搜索窗口中其他像素点的比值相似度；(5)计算待恢复像素点与其搜索窗口中其他像素点的相似性权值；(6)用相似度权值和搜索窗口中像素点进行加权求和，重构待恢复像素点的非局部均值；(7)得到基于比值相似度的非局部均值的非局部均值比差异图。

Description

一种SAR图像变化检测差异图生成方法

技术领域

本发明属于SAR图像变化检测领域，涉及SAR图像变化检测中的差异图生成技术。

背景技术

近年来，基于同一场景不同时间的SAR图像变化检测技术在遥感领域具有大量应用。SAR图像虽然具有其数据复杂性及斑点噪声的存在，但SAR传感器可以在所有天气条件下正常工作，不受大气和光照条件的影响，使得SAR图像变化检测很具有吸引力。

SAR图像变化检测一般可分为三个步骤：（1）图像预处理。（2）差异图生成。（3）差异图分析。预处理步骤主要包括配准、几何校正及图像去噪等。在第二步，两幅图像将进行比较并生成差异图，在差异图中变化区域的像素和非变化区域的像素在灰度上将呈现出明显的差异。在第三步，通常对差异图进行阈值或聚类分割就可以得到被分为变化类和不变化类的最终二值检测结果。由此看出，差异图的精度高低直接影响着变化检测的性能。

经典的差异图生成方法有差值法和比值法。差值法主要思想是将配准后、不同时间获取的两幅影像中对应像素的灰度值相减，得到差异图像。但差值法对SAR图像的乘性噪声不鲁棒。比值法是将两幅影像对应像素点进行除法运算。两幅图像均做了相对辐射校正等预处理，产生的比值法差异图使变化区域的信息得到了增强。

国内外学者在比值法基础上提出了目前流行的对数比和局部均值比差异图产生方法。对数比方法对比值差异图像取对数，能使SAR图像的乘性噪声转化为便于处理的加性噪声。局部均值比方法对两幅原图进行局部均值滤波后再做比值运算，能够一定程度上抑制SAR图像噪声。

非局部均值算法在图像去噪领域非常流行，它通过整幅图的像素加权平均值来估计一个像素点，而权重则通过测量其他像素点邻域图块与此像素点邻域图块的相似度来得到。非局部均值利用了一幅图像中的所有相似冗余图块来进行去噪。同时，一个纹理或边缘图块可通过图像中其他相似的纹理或边缘图块得以完好的保存甚至修复。但是，传统非局部均值基于高斯距离的图块相似度测量方法不适用于SAR图像块的测量。

目前最流行的差异图生成方法中，对数比方法容易漏检变化区域，同时该方法对图像的配准精度要求很高。局部均值比方法使用局部均值滤波时，容易过度平滑SAR图像，丢失边缘和纹理细节，误检率较高。而传统非局部均值的图块相似度测量方法不适用于SAR图像。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种基于比值相似度非局部均值的SAR图像变化检测差异图生成方法，以实现在SAR图像变化检测中既能抑制乘性斑点噪声的影响，又能较好地保持图像边缘及纹理细节，使得差异图质量更好，变化检测结果精度更高。

本发明提出一种新的SAR图像块相似度测量方法，其更加鲁棒于SAR图像。本发明的技术方案如下：

步骤1.通过合成孔径雷达获取两幅不同时相同样大小的SAR图像，对这两幅SAR图像进行预处理，预处理不改变图像的大小，获得预处理后的适合变化检测使用的SAR图像X¹和X²；

步骤2.分别对图像X¹和图像X²进行全图逐像素遍历，选取第i个位置的像素点X_i作为待恢复像素点，i＝1,2,…,m×n，这里X_i泛指两幅图像X¹和X²第i个位置的像素点，其中m，n分别是步骤1中获取的两幅不同时相同样大小的SAR图像的长度和宽度，选取以X_i为中心的正方形W_i ^r作为搜索窗口，选取以X_i为中心的正方形图块作为X_i的相似性测量图块，的大小应小于搜索窗口W_i ^r的大小；

步骤3.在所述两幅图像X¹和X²中，对于以X_i为中心的搜索窗口W_i ^r中其他像素点X_j,X_j∈W_i ^r，选取以X_j为中心的正方形图块作为像素点X_j的相似性测量图块，正方形图块与正方形图块的大小相等，计算以像素点X_i为中心的图块和以Xj为中心的图块这两个相似性测量图块对应的矩阵对应位置处的像素点灰度比值r_i,j的概率密度分布函数p(r_i,j)；

步骤4.对概率密度分布函数p(r_i,j)进行高斯核函数映射得到以像素点X_i为中心的图块和以X_j为中心的图块这两个相似性测量图块对应的矩阵对应位置处的像素点比值相似度R(i,j)；

步骤5.在两幅图像X¹和X²中，分别计算待恢复像素点X_i与其搜索窗口W_i ^r中所有其他像素点X_j的相似性权值ω(i,j)；

步骤6.对于像素点X_i，用相似性权值ω(i,j)与搜索窗口W_i ^r中第j位置的像素点X_j的灰度值加权相乘，对加权乘积求和，得到两幅图像X¹和X²第i个位置的待恢复像素点X_i的基于比值相似度的非局部均值

步骤7.使用步骤6中得到的基于比值相似度的非局部均值NL[X_i]，得到非局部均值比差异图中对于像素点X_i的像素点

其中，为SAR图像X¹和X²的差异图中第i位置像素点，与分别表示SAR图像X¹与X²的第i位置像素点，与表示SAR图像X¹和X²的第i位置像素点灰度值的非局部均值。

本发明提出了一种比值距离用来测量受斑点噪声影响的SAR图像块间的距离，这种比值距离的概率密度函数被用来将其映射为相似度值。然后将这种新的非局部均值引入局部均值比算子，提出了一种新的差异图生成方法。本发明与现有技术相比具有以下优点：

第一，本发明提出了一种鲁棒于SAR图像的比值相似度测量方法改进了非局部均值。更详细地说，提出了一种比值距离用来测量受斑点噪声影响的SAR图像块间的距离，这种比值距离的概率密度函数被用来将其映射为相似度值。

第二，本发明采用了一种新的非局部均值方法对SAR图像进行差异图生成，这种非局部均值方法不仅能很好的抑制SAR图像的噪声，还能利用全图的边缘纹理信息，很好的保持检测的边缘纹理信息。

附图说明

图1为本发明的流程图；

图2为本发明提出的比值概率密度函数随图像视数变化情况图；

图3为本发明所用的Ottawa地区两幅SAR影像和变化参考图像；

图4为本发明采用的方法和其它三种比较方法对Ottawa地区SAR影像的检测结果图；

图5为本发明对Ottawa数据集的ROC曲线和其它方法的ROC曲线比较图；

图6为本发明所用的黄河地区两幅SAR影像和变化参考图像；

图7为本发明采用的方法和其它三种比较方法对黄河地区SAR影像的检测结果图；

图8为本发明对黄河数据集的ROC曲线和其它方法的ROC曲线比较图。

具体实施方式

下面结合附图，对本发明做进一步的描述。

参照附图1，本发明具体实施步骤如下：

步骤1.通过合成孔径雷达获取两幅不同时相同样大小的SAR图像，对这两幅SAR图像进行预处理，预处理不改变图像的大小，获得预处理后的适合变化检测使用的SAR图像X¹和X²。

所谓两幅不同时相的SAR图像，是指同一时间、不同场景下拍摄的两幅SAR图像。

所述对这两幅SAR图像进行预处理包括对两幅SAR图像进行图像配准和几何校正。

图像配准是对不同时间、不同视角的同一地区的两幅图像利用已知信息进行匹配、校准的过程。

几何校正是根据图像坐标和地面坐标之间的数学模型，确定校正公式，对原始影像进行几何变换计算。

具体实现时，可以采用输入到Matlab2010图像处理库和C++计算机视觉库软件的方式。

步骤2.分别对图像X¹和图像X²进行全图逐像素遍历，选取第i个位置的像素点X_i作为待恢复像素点，i＝1,2,…,m×n，这里X_i泛指两幅图像X¹和X²第i个位置的像素点，其中m，n分别是步骤1中获取的两幅不同时相同样大小的SAR图像的长度和宽度。选取以X_i为中心的正方形W_i ^r作为搜索窗口。对于每个待恢复像素点X_i，选取以其为中心的正方形图块作为X_i的相似性测量图块，的大小应小于搜索窗口W_i ^r的大小。

之所以称为待恢复像素点，是因为此像素点原本可能被噪声污染，我们对其进行非局部均值运算，可以对其进行恢复。

理论上，在小于整幅图像大小的前提下，搜索窗口W_i ^r越大越好。但考虑到算法效率，根据试验结果，大小为21×21的搜索窗口已足够满足搜索性能。因此，搜索半径可根据(21-1)/2计算得到r＝10。

考虑到算法效率与性能，本发明设置相似性测量图块大小为7×7。

步骤3.在所述两幅图像X¹和X²中，对于以X_i为中心的搜索窗口W_i ^r中其他像素点X_j,X_j∈W_i ^r，选取以X_j为中心的正方形图块与的大小相等）作为像素点X_j的相似性测量图块。计算以像素点X_i为中心的图块和以X_j为中心的图块这两个相似性测量图块对应的矩阵对应位置处的像素点灰度比值r_i,j＝X_i/X_j，接着根据SAR图像符合的伽马（gamma）分布假设推导像素点灰度比值的概率密度分布函数。

伽马分布假设形式如下：

其中，Y代表SAR图像的像素点灰度的观测值，p_Y(Y)表示伽马分布的概率密度分布函数，SAR图像的像素点灰度的观测值的分布符合此概率密度分布函数，L代表视数，X代表该SAR图像的雷达反射值，Γ(·)为伽马函数，其表达式为

通常我们处理两幅视数均为L的SAR图像Y₁和Y₂时，假设这两幅图像间的复杂相关度是0，则像素点灰度比值r_i,j＝Y₁/Y₂的概率密度分布函数如下：

从而计算出以像素点X_i为中心的图块和以X_j为中心的图块这两个相似性测量图块对应的矩阵对应位置处的像素点灰度比值r_i,j的概率密度分布函数p(r_i,j)。

研究概率密度函数p(r_i,j)，我们发现其最大值在r^*随着视数的增大越来越接近于1。图2显示了r^*随视数变化的情况。

步骤4.对概率密度分布函数p(r_i,j)进行高斯核函数映射得到以像素点X_i为中心的图块和以X_j为中心的图块这两个相似性测量图块对应的矩阵对应位置处的像素点比值相似度R(i,j)：

其中，公式中（./）符号表示矩阵点除运算，即矩阵对应位置元素相除；M为相似性测量图块和中像素点的个数；G与G(k)表示高斯核函数，形式为G(k)＝1/[2(k-3)+1]²；分别表示图块和的第k个对应位置的像素的灰度值，k＝1,2，…，M；R(i,j)表示图块和的相似度，也就是图块和中心两个像素点X_i和X_j的相关性。R(i,j)越大表示图块和图块越相似。

步骤5.在两幅图像X¹和X²中，分别计算待恢复像素点X_i与其搜索窗口W_i ^r中其他像素点X_j,X_j∈W_i ^r的相似性权值：

其中，exp表示指数运算，Z(i)＝Σexp(-R(i,j)/h²)表示标准化常数，其作用是保证相似性权值ω(i,j)满足约束条件0≤ω(i,j)≤1和Σω(i,j)＝1。h是控制指数函数衰减速度的滤波参数，我们令σ²表示斑点噪声的方差，则h与σ的线性关系可以被估计为h＝fσ，f为常数。对于L视数的SAR图像，我们可得知因此h可以被固定赋值为一些文献中将f固定为2，而在本发明中对于真实SAR图像将其取1.4效果较好。

步骤6.对于像素点X_i，用相似性权值ω(i,j)与搜索窗口W_i ^r中第j位置的像素点X_j的灰度值加权相乘，即ω(i,j)X_j。类似的，将搜索窗口W_i ^r中其他所有像素点进行加权相乘，然后对以上所有加权乘积求和，即可得到两幅图像X¹和X²第i个位置的待恢复像素点X_i的灰度值即待恢复像素点的基于比值相似度的非局部均值，其中W_i ^r为以像素X_i为中心的搜索窗口。对两幅图像X¹和X²进行全图逐像素遍历，计算每个像素点的基于比值相似度的非局部均值。

例如，像素点X_i的搜索窗口W_i ^r是一个21*21的正方形图块，包含了21*21=442个像素点，搜索窗口W_i ^r中第1个位置像素点X₁的灰度值需要与相似性权值ω(i,1)相乘。搜索窗口W_i ^r中第2个位置像素点X₂的灰度值需要与相似性权值ω(i,2)相乘。以此类推，搜索窗口W_i ^r中第442个位置像素点X₄₄₂的灰度值需要与相似性权值ω(i,442)相乘。对以上442个乘积求和，就可以得到像素点X_i的非局部均值。

步骤7.使用步骤6中得到的基于比值相似度的非局部均值NL[X_i]，得到非局部均值比差异图，其形式为其中，为SAR图像X¹和X²的差异图中第i位置像素点的灰度。与分别表示SAR图像X¹与X²的第i位置像素点的灰度值。与表示图像X¹与X²的第i位置像素点灰度值的非局部均值。

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

1.仿真条件：

本发明是在CPU为Intel Core22.3GHz，内存2G，Windows7系统上使用Matlab2010进行仿真。

2.仿真内容：

我们采用三种对比方法，包括经典的局部均值比算子(MRD)和对数比算子(LRD)，其中对数比差异图形式为我们使用提出的非局部均值对其进行改进，得到第三种对比方法非局部对数比算子(NLLRD)：我们采用四种差异图生成方法分别对带参考图的Ottawa地区和黄河地区真实SAR图像数据进行了实验。

仿真1：

参照图3，图3为Ottawa地区真实SAR影像，(a)表示1997年5月Ottawa地区的地貌信息。(b)表示1997年8月Ottawa地区的地貌信息。(c)表示变化检测的变化对比图。其中黑色为水域，灰色为河岸。

参照图4，图4为本发明采用的方法和其它三种比较方法对Ottawa地区SAR影像生成的差异图；

(a)为MRD差异图，边缘模糊，细节保持不好，图像噪声较多。

(b)为LRD差异图，没有使用邻域信息的LRD方法也存在很大的斑点噪声。

(c)和(d)为本发明的NLMRD方法和NLLRD方法，无论在变化区域还是非变化区域都更加清晰有效。相比较于MRD差异图，这两幅差异图的边缘更加清楚并且图中的斑点噪声也得到了很好的抑制，主要原因是我们引入的非局部均值在SAR图像去噪及细节保持方面都比MRD使用的局部均值比效果更好。相比较于没有使用邻域信息的LRD方法，我们提出的差异图构造方法生成的差异图背景明显光滑的多。

参照图5，图5为本发明对Ottawa数据集的ROC曲线和其它方法的ROC曲线比较图。ROC曲线下面积越大表示差异图效果越好。从图5我们可计算出MRD、LRD、NLMRD和NLLRD方法的ROC曲线下面积分别为0.9867,0.9406,0.9969和0.9962。通过精确的数字我们可以看到引入我们所提出的非局部均值的NLMRD和NLLRD两种方法明显优于两种经典方法，而且NLMRD的效果还稍微优于NLLRD方法。

仿真2：

参考图6，图6为本发明所用的黄河地区两幅SAR影像和变化参考图像；(a)表示2008年6月18日黄河地区地貌信息。(b)表示2009年6月19日黄河地区地貌信息。(c)表示变化检测的变化对比图。中间方形部分为农田，周围灰白色为河岸，下方黑色为水域。

参考图7，图7为本发明采用的方法和其它三种比较方法对黄河地区SAR影像的检测结果图；

(a)为MRD差异图，在背景上都存在大量杂乱的斑点并且边缘也不是很清晰。

(b)为LRD差异图，包含很大的斑点噪声，变化区域边缘不够清晰。

(c)和(d)为本发明的NLMRD方法和NLLRD方法，从视觉观察来看，NLMRD和NLLRD差异图比其他差异图更有效的突出了变化区域和非变化区域。很好的反映了农田的形状和纹理信息，而且也明显对背景部分的噪声有很大抑制。而LRD和MRD差异图在非变化区域仍然包含着很严重的噪点，这在之后的差异图分析阶段会产生很严重的影响。

参考图8，图8为本发明对黄河数据集的ROC曲线和其它方法的ROC曲线比较图。ROC曲线下面积越大表示差异图效果越好。MRD、LRD、NLMRD和NLLRD差异图的ROC曲线下面积分别是0.9181,0.9014,0.9317和0.9276。根据下面积指标我们可以看到引入非局部均值的两种算法明显表现出了更好的性能。

综上所述，MRD算法和LRD算法由于在计算过程中没有针对SAR图像的噪声特性对乘性相干斑噪声采取足够的抑制措施，效果不佳。而本发明提出的算法，全图信息的支持相比于局部窗口信息更有效的抑制了乘性噪声，并且更好的保存了细节信息，实验结果也进一步证明了本方法的有效性。

Claims (7)

1.一种SAR图像变化检测差异图生成方法，包括以下步骤：

步骤1.通过合成孔径雷达获取两幅不同时相同样大小的SAR图像，对这两幅SAR图像进行预处理，预处理不改变图像的大小，获得预处理后的适合变化检测使用的SAR图像X¹和X²；

步骤2.分别对图像X¹和图像X²进行全图逐像素遍历，选取第i个位置的像素点X_i作为待恢复像素点，i＝1,2,…,m×n，这里X_i泛指两幅图像X¹和X²第i个位置的像素点，其中m，n分别是步骤1中获取的两幅不同时相同样大小的SAR图像的长度和宽度，选取以X_i为中心的正方形作为搜索窗口，选取以X_i为中心的正方形图块作为X_i的相似性测量图块，的大小应小于搜索窗口的大小；

步骤3.在所述两幅图像X¹和X²中，对于以X_i为中心的搜索窗口中其他像素点X_j,选取以X_j为中心的正方形图块作为像素点X_j的相似性测量图块，正方形图块与正方形图块的大小相等，计算以像素点X_i为中心的图块和以X_j为中心的图块这两个相似性测量图块对应的矩阵对应位置处的像素点灰度比值r_i,j的概率密度分布函数p(r_i,j)；

步骤4.对概率密度分布函数p(r_i,j)进行高斯核函数映射得到以像素点X_i为中心的图块和以X_j为中心的图块这两个相似性测量图块对应的矩阵对应位置处的像素点比值相似度R(i,j)；

步骤5.在两幅图像X¹和X²中，分别计算待恢复像素点X_i与其搜索窗口中所有其他像素点X_j的相似性权值ω(i,j)，其中，

$\mathrm{\ω}(i,j)=\frac{\exp (-\frac{R(i,j)}{h^2})}{Z\left(i\right)}$

exp表示指数运算，表示标准化常数，h是控制指数函数衰减速度的滤波参数，σ²表示斑点噪声的方差，h与σ的线性关系可以被估计为h＝fσ，f为常数，对于L视数的SAR图像，h可以被固定赋值为

步骤6.对于像素点X_i，用相似性权值ω(i,j)与搜索窗口中第j位置的像素点X_j的灰度值加权相乘，对加权乘积求和，得到两幅图像X¹和X²第i个位置的待恢复像素点X_i的基于比值相似度的非局部均值

步骤7.使用步骤6中得到的基于比值相似度的非局部均值NL[X_i]，得到非局部均值比差异图中对于像素点X_i的像素点

其中，为SAR图像X¹和X²的差异图中第i位置像素点，与分别表示SAR图像X¹与X²的第i位置像素点，与表示SAR图像X¹和X²的第i位置像素点灰度值的非局部均值。

2.根据权利要求1所述的SAR图像变化检测差异图生成方法，其中，所述对这两幅SAR图像进行预处理包括对两幅SAR图像进行图像配准和几何校正。

3.根据权利要求1所述的SAR图像变化检测差异图生成方法，其中，搜索窗口的大小为21×21。

4.根据权利要求1所述的SAR图像变化检测差异图生成方法，其中，相似性测量图块大小为7×7。

5.根据权利要求1所述的SAR图像变化检测差异图生成方法，其中，当SAR图像X¹和X²的视数均为L，如果这两幅SAR图像间的复杂相关度是0，则以像素点X_i为中心的图块和以X_j为中心的图块这两个相似性测量图块对应的矩阵对应位置处的像素点灰度比值r_i,j的概率密度分布函数p(r_i,j)如下：

Γ(·)为伽马函数。

6.根据权利要求1所述的SAR图像变化检测差异图生成方法，其中，

$R(i,j)=\left|\right|p(X_{N_i}./X_{N_j})|{|}_{2,G}^2=\underset{k=1}{\overset{M}{\Σ}}G\left(k\right){\[p\left(\frac{X_{N_i}\left(k\right)}{X_{N_j}\left(k\right)}\right)\]}^2$

公式中(./)符号表示矩阵点除运算，M为相似性测量图块和中像素点的个数，G与G(k)表示高斯核函数，形式为G(k)＝1/[2(k-3)+1]²，分别表示图块和的第k个对应位置的像素的灰度值，k＝1,2，…，M。

7.根据权利要求1所述的SAR图像变化检测差异图生成方法，其中，f＝1.4。