基于自适应权值图像融合的SAR图像变化检测方法

摘要

本发明公开了一种基于自适应权值图像融合的SAR图像变化检测方法，主要解决现有技术中单一类型差异图检测精度低的问题。其实现步骤为：1.读入不同时刻同一地区的两幅SAR图像，分别进行PPB滤波得到滤波后的图像X1和X2；2.根据滤波后的图像，得到差值图Ds与对数比值图Dl；3.分别对Ds和Dl进行均值滤波，得到均值滤波后的差值图Ds'与中值滤波后的比值图Dl'；4.计算融合参数η(i,j)；5.根据融合参数对差值图Ds'与比值图Dl'进行融合产生差异图D；6.将差异图D聚成不同的两类，得到变化检测结果。本发明具有操作简单、抗噪性好和检测精度高的优点，可应用于环境监测、灾害预测。

说明书

技术领域

本发明属于图像处理技术领域，涉及SAR图像变化检测方法，可用于城区规划、 自然灾害评估、预测灾害、土地利用和土地覆盖的动态监测等领域。

背景技术

合成孔径雷达图像又叫SAR图像，属于一种微波雷达成像，与可见光和红外遥 感相比，微波遥感具有无可比拟的优点：一、微波能穿透云雾、雨雪，具有全天候、 全天时的工作能力。二，微波对地物有一定穿透能力。三，采用侧视方式成像，覆盖 面积大。由于SAR图像这些优点，使得SAR图像日益成为变化检测的重要数据源。 遥感图像变化检测是指通过对同一地区不同时期的两幅或多幅图像进行比较分析，根 据图像之间的差异来获取地物的变化信息。变化检测技术可检测出不同时期遥感图像 辐射值及局部纹理的变化信息，该特性使得遥感图像在资源和环境监测，如土地利用 和覆盖变化，森林和植被变化，城市扩展等，地理空间数据更新，农业作物生长监测， 自然灾害监测与评估，军事领域研究等方面得到的广泛应用。

在变化检测研究中，学者们将现有的多种变化检测方法归为不同的类别，其中最 常见的是基于差异图像分析的变化检测方法。该方法通常包括3个关键的步骤：(1) 图像预处理；(2)差异图的构造；(3)变化信息的提取。其中：

差异图构造的好坏对变化检测的结果有着重要的影响，差值法和比值法作为两种 最基本的方法，分别通过对校正后的两时相遥感图像逐像素相减和相除运算来获得差 异图。差值法构造差异图时，差值为0或接近0的认为是不变区域，不为0的认为是 变化区域，该方法理论简单、直接，容易掌握，可以较好地保留细节信息，但缺点是 易受SAR图像质量和相干斑影响，产生“伪变化”信息。SAR图像中由于相干斑噪 声的影响，多采用比值法构造差异图，且比值法对校准误差不敏感。目前在用比值法 构造差异图时主要采用对数比或均值比的形式。Dekker提出使用对数比的方法构造 差异图，对数比法就是对比值法取对数，它除了可以将乘性噪声转化为加性噪声外， 还对比值图像的变化范围进行了压缩，因此削弱了高像素值区域的变化，部分真实变 化区域轮廓的像素灰度值与非变化类像素的灰度值比较接近,不利于保留变化区域的 轮廓和提高变化检测精度。均值比法就是取对应像素的邻域均值再求比值，它对噪声 具有更强的鲁棒性，但是对背景区域也就是未变化区域带来过多的伪变化信息。因此， 使用单一类型的差异图进行变化检测会存在检测精度低以及适用范围窄等问题。

因此将单一类型的差异图进行融合在变化检测中起到至关重要的作用，现有的将 单一类型的差异图进行融合的方法有基于参数权重的差异图融合和离散小波融合等 方法。其中：

离散小波变换融合的方法为：对均值比图和对数比值图分别进行离散小波分解, 得到一系列子频带图像，即低频分量，水平方向高频分量，竖直方向高频分量，对角 线方向高频分量，在不同的尺度上对各分解层进行融合处理,从而形成融合图像的小 波分解表示，最后通过离散小波逆变换得到融合图像。该方法的不足之处是融合的过 程比较复杂。

基于参数权重的差异图融合方法，其过程为：对两幅待融合的图像分别得到他们 的差值图与对数比值图，提取差值图D_s中保留细节信息较好的部分与对数比值图D_l中的抗噪性能较好的部分，结合两类差异图的优点，按照不同的权重参数α进行融合， 融合后的差异图为D，融合公式为D＝α*D_s+(1-α)*D_l，将差值图与对数比值图根 据不同的权重进行融合，取α的值从0-1范围进行多次试验，选取融合效果最好的值 作为最终α的值得到最终融合后的差异图。该方法的不足之处为：(1)在没有标准参 考图的情况下，无法判断融合结果的好坏，其适用范围窄。(2)在参数α的选择时需 要进行多次试验来对比融合的效果，参数α不能自适应，需要手动调整参数，操作繁 琐。

发明内容

本发明目的在于提出一种基于自适应权值图像融合的SAR图像变化检测方法， 以提高检测的精度，同时解决现有技术在图像融合时参数难以选取，参数不能自适应 的问题。

为了实现上述目的，本发明实现的步骤包括如下：

(1)读入在不同时刻同一地区获取的两幅SAR图像I₁和I₂，并对这两幅SAR图 像I₁和I₂分别进行基于块的权重概率PPB滤波，得到滤波后的两幅图像X₁和X₂；

(2)计算滤波后的两幅图像X₁和X₂所对应像素灰度值的差，并进行归一化，得 到差值图D_s；

(3)计算滤波后的两幅图像X₁和X₂对应像素灰度值的商，并进行归一化，得到 对数比值图D_l；

(4)将差值图D_s进行11×11窗口的均值滤波，消除差值图D_s中的孤立像素点， 得到均值滤波后的差值图D_s'；

(5)将对数比值图D_l进行3×3窗口的中值滤波，抑制对数比值图D_l中的孤立像 素点，得到中值滤波后的对数比值图D_l'；

(6)分别计算差值图D_s'中每一个像素点x(i,j)的3×3邻域Ω_x像素的均值 μ_x(i,j)与方差σ_x(i,j)，得到融合参数：

$\eta (i,j)=\frac{{\sigma }_x(i,j)}{{\mu }_x(i,j)},$

该融合参数η(i,j)随着像素变化而变化，用于体现测量像素点x(i,j)所处的3*3邻域 Ω_x的异质性，即匀质区域或异质区域；匀质区域是指图像中平滑的区域，异质区域 是指图像中的噪声部分或边缘部分，像素点x(i,j)处于异质区域时的融合参数η(i,j) 大于其处于匀质区域时的融合参数η(i,j)；

(7)根据得到的融合参数η(i,j)将均值滤波后的差值图D_s'与中值滤波后的对数 比值图D_l'进行图像融合，得到融合后的差异图D，该差异图D中每一个横纵坐标在 i，j处的像素点为：

D(i,j)＝η(i,j)×D_s′(i,j)+(1-η(i,j))×D_l′(i,j)，

其中，D'_s(i,j)为均值滤波后差值图中横纵坐标分别为i，j的点，D′_l(i,j)为中值滤波 后对数比值图中横纵坐标分别为i，j的点；

(8)用k-means聚类方法将融合后的差异图D聚类成不同的两类，分别计算这两 个不同类别的均值，定义均值较大的那一类为变化类，均值较小的那一类为非变化类， 得到最终的变化检测结果。

本发明与现有技术相比具有以下优点：

1.本发明在对均值滤波后差值图D'_s(i,j)与中值滤波后对数比值图D′_l(i,j) 融合时利用D'_s(i,j)中每一个像素3*3邻域Ω_x的方差σ_x(i,j)与均值μ_x(i,j)的比值 η(i,j)作为融合参数，使原有的融合参数能根据不同的像素点自动调整，无需进行多 次试验手动调整融合参数，操作简单。

2.本发明采用的融合参数η(i,j)考虑了邻域像素的灰度信息与空间信息， 提高了抗噪性能与保留图像细节的能力，进而提高变化检测精度。

仿真实验表明，本发明对多幅SAR图像都能得到满意的变化检测结果，具有较好 的鲁棒性。

附图说明

图1是本发明的实现流程图；

图2是本发明仿真使用的黄河SAR图像数据集图；

图3是用本发明和对比方法对图2进行仿真实验的结果图；

图4是本发明仿真使用的渥太华地区SAR图像数据集图；

图5是本发明和对比方法对图4仿真实验的结果图；

图6是本发明仿真使用的伯尔尼地区SAR图像数据集图；

图7是本发明和对比方法对图6仿真实验的结果图。

具体实施方式

参照图1，本发明的具体步骤如下：

步骤1：读入在不同时刻获取同一地区的两幅SAR图像I₁和I₂。

步骤2：对两幅SAR图像I₁和I₂分别进行PPB滤波得到滤波后的图像X₁和X₂。

对第一幅SAR图像I₁进行probabilistic patch-based weights,即PPB滤波，得到 滤波后的图像X₁；

对第二幅SAR图像I₂进行PPB滤波得到滤波后的图像X₂；

所述PPB滤波，是一种基于噪声分布模型，具有更普遍的统计相似性准则的滤 波方法。PPB滤波时，去除噪声的过程被表示成一个加权最大似然估计问题，这些权 重在迭代过程中基于噪声块的相似性与之前计算的块的相似性进行更新，这种迭代过 程提高了去噪的性能。

步骤3：获取两幅SAR图像的差值图D_s。

3a)对第一幅SAR图像X₁在坐标(i,j)处的灰度值X₁(i,j)与第二幅SAR图像X₂在坐标(i,j)处的灰度值X₂(i,j)进行差值运算，得到差矩阵D在坐标(i,j)处的差值： D(i,j)＝|X₁(i,j)-X₂(i,j)|，进而得到差矩阵：D＝{D(i,j)}；

3b)归一化差矩阵，得到差值图D_s在坐标(i,j)处的值D_s(i,j)：

$D_s(i,j)=\frac{D(i,j)-\min \left(D\right)}{\max \left(D\right)-\min \left(D\right)},$进而得到差值图：D_s＝{D_s(i,j)}。

步骤4：计算滤波后的两幅图像X₁和X₂对应像素灰度值的商，并进行归一化， 得到对数比值图D_l。

4a)通过如下公式计算计算滤波后的两幅图像X₁和X₂对应像素灰度值的商 R(i,j)：

$R(i,j)=\left|\log \frac{X_2(i,j)+1}{X_1(i,j)+1}\right|,$

其中X₁(i,j)为第一幅SAR图像X₁在横纵坐标为i，j处的像素点，X₂(i,j)为第 二幅SAR图像X₂在横纵坐标为i，j处的像素点；

4b)通过如下公式对像素灰度值的商进行归一化，得到每一个像素灰度值的商归 一化后的结果：

$D_l(i,j)=\frac{R(i,j)-\min \left(R\right)}{\max \left(R\right)-\min \left(R\right)},$其中R＝{R(i,j)}；

4c)用每一个像素灰度值的商归一化后的结果得到对数比值图：D_l＝{D_l(i,j)}。

步骤5：将差值图D_s进行11×11窗口的均值滤波，消除差值图D_s中的孤立像素 点，使得区域更完整，局部区域保持较好的连续性，图像更加平滑，最终得到均值滤 波后的差值图D_s'。

步骤6：将对数比值图D_l进行3×3窗口的中值滤波，抑制对数比值图D_l中的孤 立像素点，以保留边缘信息，最终得到中值滤波后的对数比值图D_l'。

步骤7：计算融合参数η(i,j)。

7a)利用如下公式计算差值图D_s'中每一个像素点x(i,j)的3×3邻域Ω_x像素的均 值μ_x(i,j)

${\mu }_x(i,j)=\frac{\underset{n=1}{\overset{9}{\Sigma }}{x}_n(i,j)}{9},$

其中，为3*3邻域Ω_x内每一个像素点x(i,j)的灰度值之和；

7b)根据均值μ_x(i,j)，利用如下公式计算差值图D_s'中每一个像素点x(i,j)的3×3 邻域Ω_x像素的方差σ_x(i,j)：

${\sigma }_x(i,j)=\frac{\underset{n=1}{\overset{9}{\Sigma }}{(x_n(i,j)-{\mu }_x(i,j))}^2}{9},$

其中，为3*3邻域Ω_x内每一个像素点x(i,j)灰度值减去邻域 Ω_x内均值μ_x(i,j)结果的平方之和；

7c)根据均值μ_x(i,j)和方差σ_x(i,j)计算融合参数η(i,j)：

$\eta (i,j)=\frac{{\sigma }_x(i,j)}{{\mu }_x(i,j)},$

该融合参数η(i,j)随着像素变化而变化，其变化范围在0-1之间，用于体现测量 像素点x(i,j)所处的3*3邻域Ω_x的异质性，即匀质区域或异质区域；匀质区域是指 图像中平滑的区域，异质区域是指图像中的噪声部分或边缘部分，像素点x(i,j)处于 异质区域时的融合参数η(i,j)大于其处于匀质区域时的融合参数η(i,j)。

步骤8：融合均值滤波后的差值图D_s'与中值滤波后的对数比值图D_l'，得到融合 后的差异图D，

根据步骤7得到的融合参数η(i,j)，将均值滤波后的差值图D_s'与中值滤波后的 对数比值图D_l'进行图像融合，得到融合后的差异图D，该差异图D中每一个横纵坐 标在i，j处的像素点为：

D(i,j)＝η(i,j)×D_s'(i,j)+(1-η(i,j))×D_l'(i,j)，进而得到融合后的差异图 D＝{D(i,j)}，

其中，D'_s(i,j)为均值滤波后差值图中横纵坐标分别为i，j的点，D′_l(i,j)为中值 滤波后对数比值图中横纵坐标分别为i，j的点；由于系数η(i,j)是测量x(i,j)邻域区 域异质性的参数，当η(i,j)值小时x(i,j)处于匀质区域，D_l'在产生差异图时占主导地 位，更好地抑制图像的斑点噪声，当η(i,j)值大时x(i,j)处于异质区域，D_s'在产生差 异图时占主导地位，更好地保留图像细节信息。

步骤:9：对差异图D进行k-means聚类，得到最终变化检测结果。

9a)从差异图D中随机选择2个像素点作为初始聚类中心；

9b)分别计算差异图D中每一个像素点x(i,j)与两个聚类中心的距离，并根据最 小距离重新对相应像素点进行划分，把x(i,j)归到与它距离最近的那个类中去；

9c)重新计算每个类的均值作为新的聚类中心；

9d)重复步骤9b)-9c)直到每个聚类不再发生变化为止，得到两个类别的均值；

9e)将两个类别中均值小的那一类像素转化为0即为未变化类，均值大的那一类 像素转化为255即为变化类，输出得到的结果图，由此得到最终的变化检测结果。

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

1.实验条件

实验环境为：windows XP，SPI，CPU Pentium(R)4，基本频率2.4GHZ，软 件平台为MatlabR2010a。

第一个数据集为黄河的SAR图像数据集的一部分截取，如图2(a)和2(b)所 示，在山东省东营被采集，分辨率为8m*8m。其中2(a)为2008年6月的图像，2 (b)为2009年6月的图像。图像大小均为257×289个像素，图像灰度级为256。 图2(c)为标准变化检测结果图，其中包括13423个变化像素点和60841个未变化 像素点。

第二个数据集为渥太华地区SAR图像数据集，它由两幅在不同时刻拍摄的 Radarsat SAR图像组成，变化区域主要是由洪水灾情引起的，如图4(a)和4(b)所示， 其中图4(a)为1997年5月的图像，图4(b)为1997年8月的图像，图像大小均为290 ×350像素，灰度级为256。图4(c)为该数据集的标准变化检测图，包括85451个未 变化像素点和16049个变化像素点。

第三个数据集为伯尔尼地区SAR图像数据集，它由两幅在不同时刻拍摄的 Radarsat SAR图像组成，如图6(a)和6(b)所示，其中图6(a)为1999年4月 图像，图6(b)为1999年5月的图像，图像大小均为301×301像素，灰度级为256。 图6(c)为该数据集标准变化检测结果图，包括1269个变化像素点，89332个未变 化像素点。

2.实验评价指标

实验使用的评价指标是漏检数FN、误检数FP、总错误数OE、正确检测率PCC 和Kappa系数。其中，漏检数为没有检测出来的实际发生了变化的像素的总和，误检 数为实际没有发生变化但被当作变化检测出来的像素的总和，总错误数是漏检数和误 检数之和，正确检测率＝1-总错误数/总像素个数。总像素数为N，总变化点数为Nc， 总未变化点数为Nu，正确检测出的变化点为TP，TP＝总变化点-漏检数，正确检测出 的未变化点为TN，TN＝总未变化点-误检数，令得到如果变化图和参考图完全一致，则Kappa系数为1，否则 为0。Kappa系数包含了更多分类的细节信息，所以是比PCC更精确的评价标准。

3.实验内容

用本发明方法和现有的4种变化检测方法对3个图像数据集进行变化检测，该4 种对比方法分别为：采用PPB滤波后的图像X₁和X₂进行差值运算构造差异图然后用 K-means聚类的方法记作差值法；采用PPB滤波后的图像X₁和X₂进行对数比运算构 造差异图然后用K-means聚类的方法记做对数比值法；采用PPB滤波后的图像X₁和 X₂进行均值比运算构造差异图然后用K-means聚类的方法记作均值比法； Yaoguozheng在文章“Using Combined Difference Image and K-means Clustering for  SAR Image Change Detection”提出的方法记做CDI-K法。

实验1

用本发明方法和现有的4种变化检测方法对第一个数据集图2(a)和2(b)进行变化 检测，2(c)为参考图，实验结果如图3，其中3(a)是差值法的变化检测结果图，3(b) 是对数比值法的变化检测结果图，3(c)是均值比法的变化检测结果图，3(d)是CDI-K 方法的变化检测结果图，3(e)是本发明的变化检测结果图。

实验2

用本发明方法和现有的4种变化检测方法对第二个数据集图4(a)和4(b)进行变化 检测，4(c)为参考图，实验结果如图5，其中5(a)是差值法的变化检测结果图，5(b) 是对数比值法的变化检测结果图，5(c)是均值比法的变化检测结果图，5(d)是CDI-K 方法的变化检测结果图，5(e)是本发明的变化检测结果图。

实验3

用本发明方法和现有的4种变化检测方法对第三个数据集图6(a)和6(b)进行变化 检测，6(c)为参考图，实验结果如图7，其中7(a)是差值法的变化检测结果图，7(b) 是对数比值法的变化检测结果图，7(c)是均值比法的变化检测结果图，7(d)是CDI-K 方法的变化检测结果图，7(e)是本发明的变化检测结果图。

4.实验结果分析

从图3、图5、图7可以看出，传统的单一类型的差异图的构造抗噪性能差，总 错误数大，检测精度低，直观效果不好。CDI-K方法利用参数α将均值滤波后的差值 图与中值滤波后的对数比值图融合，生成差异图，但是在融合的过程中参数没有自适 应，参数α需要从0-1经过多次实验选择变化检测最优结果来确定参数α的值，此过 程需要手动调整参数，较为繁琐，而且只考虑了像素的灰度信息，没有考虑邻域的空 间信息，所以抗噪性能不好，检测精度低。本发明在产生差异图时使得参数自适应， 无需手动调整参数，大大简化了实验步骤，使得方法简便且高效，而且在参数计算时 考虑了像素灰度信息也考虑了邻域的空间信息，使得抗噪能力增强，总错误数减少， 提高了变化检测的精度。

将本发明方法和现有的4种变化检测方法对3个数据集进行变化检测的数值进行 统计，结果如表1所示。

表1  本发明与现有四种方法的变化检测数值比较

由表1可看出，对上述三幅SAR图像本发明的总错误数均少于其他方法，PCC和 Kappa系数优于其他方法。