CN106097290A - 基于nmf图像融合的sar图像变化检测方法 - Google Patents

Abstract

本发明提出了一种基于NMF图像融合的SAR图像变化检测方法，用于解决现有的基于图像融合的SAR图像变化检测方法中存在的准确率低且物理意义不明确的技术问题，实现步骤如下：输入同地域不同时间的两幅原始SAR图像，经去噪处理后构造两基本差异图X_l和X_m；分别对X_l和X_m进行邻域特征描述得到特征描述矩阵X_L和X_M；利用两特征描述矩阵构造变化检测特征矩阵X；对X进行NMF非负矩阵分解，得到分解后的系数矩阵H；利用系数矩阵得到融合后差异图；利用聚类算法将融合后差异图聚成两类，得到变化检测结果。本发明更好的利用了原始图像的细节信息，且该图像融合方法具有更显著的物理意义，可用于SAR图像变化信息的分析和标注。

Description

基于NMF图像融合的SAR图像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及一种SAR图像变化检测方法，具体涉及一种基于NMF图像融合的SAR图像变化检测方法，可用于SAR图像变化信息的分析和标注。

背景技术

随着一系列光学和雷达卫星的运行服务的发展，积累了大量的遥感数据，通过对这些图像数据进行进一步的分析、处理及解译，就可以将经遥感技术获得的图像数据转化为有用的信息，进而服务于国民经济的各个领域。

凡是记录各种地物电磁波大小的图片都称为遥感图像，其中合成孔径雷达SAR(Synthetic Aperture Radar)图像是遥感图像的一个重要分支，它是一种先进的主动式微波对地观测设备，具有全天候、全天时工作的特点，对地面植被有一定穿透能力，能获得类似光学照片的目标图像，在国民经济和军事应用领域有着十分重要的作用。

SAR图像变化检测技术是遥感图像研究的重要组成部分，它是对同一地域不同时刻拍摄的多幅遥感图像，采用图像处理和模式识别等手段进行对比分析，根据图像之间的差异来得到人们所需要的地物或目标随时间发生的变化信息。可以应用于军事研究领域、遥感气象服务、灾害动态监测、环境污染检测以及城市规划等方面。

Singh.A将图像变化检测方法分为分类比较法和直接比较法，分类比较法的变化检测结果的精度容易受分类精度的影响，初始分类的错误将造成难以弥补的误检或漏检现象；直接比较法的算法相对简单，容易实现，但难以克服大气条件和传感器噪声等对检测结果带来的误差，受图像的预处理精度的影响严重。图像融合能有效的解决该问题，图像融合是指将多源信道所采集到的关于同一目标的图像数据经过图像处理和计算机技术等，最大限度的提取各自信道中的有利信息，最后综合成高质量的图像，以提高图像信息的利用率、改善计算机解译精度和可靠性、提升原始图像的空间分辨率和光谱分辨率，利于监测。如公茂果等2012年在国际电子与电气工程师协会地球科学与遥感杂志上发表的名为基于小波融合的SAR图像变化检测方法文章中，发表了一种基于三级小波融合的SAR图像变化检测方法，该方法是将对数比值差异图和均值比值差异图通过小波变换和逆变换得到融合差异图，融合规则是基于低频信息的加权平均和高频信息的最小标准差。该文章采用了三级小波变换和逆变换，得到小波融合差异图，最终实现对融合差异图的变化信息的分析和标记。该方法为图像融合提供了新的思路，但物理意义不够明确，并且三级小波运算及逆运算导致运算效率并不高。

随着图像质量的增强，以及图像数量的增多，行之有效的处理大量遥感影像数据成为亟待解决的问题。矩阵分解在降维、低秩和稀疏方面有着重要应用，且考虑图像数据本身的特点，为变化检测问题的解决提供了更清晰的物理意义。吸引了很大一批学者的兴趣。

发明内容

本发明的目的在于克服上述现有技术存在的缺陷，提供了一种基于NMF图像融合的SAR图像变化检测方法，用于解决现有的基于图像融合的SAR图像变化检测方法中存在的准确率低且物理意义不明确的技术问题。

本发明的技术关键是基于像素的原始图像特征描述与基于非负矩阵分解的图像融合方法。特征描述过程充分考虑了SAR图像的成像特点，使用像素的邻域信息来表征该像素点。同时图像融合过程利用NMF非负矩阵分解融合方法，有效利用了基本差异图像的互补信息，在有效抑制图像背景信息的同时大大增强感兴趣区域，能有效的提高检测准确率。NMF非负矩阵分解方法充分考虑了物理图像信息的非负特点，物理意义更明确。并且分解过程能实现矩阵的降维操作，提高了效率。

根据上述技术思路，实现本发明目的采用的技术方案，包括如下步骤：

1.对同一地域不同时间的两幅SAR原始图像分别进行去噪处理，得到去噪后图像I₁和去噪后图像I₂；

2.分别利用对数比算子和均值比算子，对得到的去噪后图像I₁和去噪后图像I₂进行运算，得到对数比值差异图X_l和均值比值差异图X_m，计算公式如下：

$X_l=\left|log\frac{I_2+1}{I_1+1}\right|=|log(I_2+1)-log(I_1+1)|;$

$X_m=1-min(\frac{{\mathrm{\μ}}_1}{{\mathrm{\μ}}_2},\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1});$

其中，μ₁，μ₂分别表示I₁和I₂的局部均值；

3.对得到的对数比值差异图X_l和均值比值差异图X_m进行特征描述，得到对应的特征描述矩阵X_L和特征描述矩阵X_M，实现步骤如下：

3a)分别提取得到的对数比值差异图X_l和均值比值差异图X_m中每个像素点的l层邻域像素特征矩阵，其中l≥3，且l为奇数；

3b)将提取的l层邻域特征像素矩阵展开成一列，得到该像素点特征描述列向量；

3c)分别将对数比值差异图X_l和均值比值差异图X_m中所有像素点的特征描述列向量作为不同的列依次排列，得到的特征描述矩阵X_L和特征描述矩阵X_M。

4.根据像素位置对应的原则，对得到的特征描述矩阵X_L和特征描述矩阵X_M按列拼接，得到变化检测特征矩阵：

5.采用乘性迭代的分解方式，对得到的变化检测特征矩阵X进行NMF非负矩阵分解，该NMF非负矩阵分解的矩阵表达式为X_m×n≈(WH)_m×n＝W_m×rH_r×n，其中W为基矩阵，H为系数矩阵，r为矩阵分解的秩，分解步骤如下：

5a)设矩阵分解的秩r＝1，并对基矩阵W和系数矩阵H进行随机初始化，其中基矩阵W≥0，系数矩阵H≥0；

5b)设迭代次数Iterator的初始值为0，初始化迭代误差上限为ε，迭代次数上限为Iter；

5c)验证迭代条件是否满足：误差E＝|X-WH|满足||E||≥ε并且Iterator≤Iter，若是，执行步骤5d)，否则，迭代结束，输出非负矩阵分解后的系数矩阵H，并执行步骤(6)；

5d)对基矩阵W和系数矩阵H分别进行迭代，单次迭代完成执行步骤5c)，迭代按如下步骤实现：

$W_{ia}\←W_a\underset{\mathrm{\μ}}{\Σ}\frac{X_{i\mathrm{\μ}}}{{\left(WH\right)}_{i\mathrm{\μ}}}{H}_{a\mathrm{\μ}};$

$W_{ia}\←\frac{W_{ia}}{\underset{j}{\Σ}{W}_{ja}};$

$H_{a\mathrm{\μ}}\←H_{a\mathrm{\μ}}\underset{i}{\Σ}{W}_{ia}\frac{X_{i\mathrm{\μ}}}{{\left(WH\right)}_{i\mathrm{\μ}}};$

Iterator++；

E＝|X-WH|；

6.将输出的非负矩阵分解后的系数矩阵H，按照对数比值差异图X_l或均值比值差异图X_m中对应像素的排列顺序，恢复为原始图像的尺寸，得到融合后差异图X_FUSION；

7.将得到的融合后差异图X_FUSION聚类成两类，得到同一地域不同时间的两幅SAR图像的变化检测结果图CM，并输出该结果图CM。

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

1.本发明由于在获取融合图像时，采用乘性迭代的分解方式，对得到的变化检测特征矩阵X进行NMF非负矩阵分解，充分考虑了物理图像信息的非负特点，同时在获取融合图像之前增加了对像素特征描述的步骤，能更好的利用原始图像的信息，与现有技术采用小波融合进行图像融合的方法相比，有效地提高了SAR图像变化检测的准确率，且物理意义明确。

2.本发明利用NMF非负矩阵分解方法进行图像融合的同时，实现了数据的降维操作，与现有技术中采用的将数据降维和其他操作分开处理的方法相比，简化了运算步骤，提高了运算效率。

附图说明

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

图2为本发明所使用Ottawa数据集参考图和对该数据集进行SAR图像变化检测的结果对比图；

图3为本发明所使用Berne数据集参考图和对该数据集进行SAR图像变化检测的结果对比图；

图4为本发明对Ottawa数据集进行NMF融合，和现有技术对该数据集分别进行均值比运算、对数比运算和小波融合的差异对比图以及对应差异图的变化检测结果；

图5为本发明对Berne数据集进行NMF融合，和现有技术对该数据集分别进行均值比运算、对数比运算和小波融合的差异对比图以及对应差异图的变化检测结果；

具体实施方式

以下结合附图和实施例，对本发明作进一步详细说明。

参照图1，本发明的实现步骤如下：

步骤1，获取同一地域不同时间的两幅相同大小的SAR图像，对其进行去噪处理，常见的去噪处理方法有PPB去斑去噪算法、均值滤波算法和非局部均值算法。本实施例采用均值滤波算法,均值滤波是典型的线性滤波算法，其采用的主要方法为邻域平均法。考虑到SAR图像成像原理：各像素均受到邻域像素点的散射影响，可看作是其邻域像素点在该点一定程度的叠加，因此对待处理的当前像素点，采用其邻域均值来代替该像素点，并且均值滤波算法实现简单，因此本实施例中采用均值滤波算法对两幅原始SAR图像分别进行去噪。均值滤波的基本原理是用均值代替原图像中的各个像素值，即对待处理的当前像素点(x,y)，选择一个模板，该模板由其近邻的若干像素组成，求模板中所有像素的均值，再把该均值赋予当前像素点(x,y)，作为处理后图像在该点上的灰度值g(x,y)，即m为该模板中包含当前像素在内的像素总个数，f(x,y)为原始图像中像素灰度值，去噪后图像记为I₁和I₂。

步骤2，分别利用对数比算子和均值比算子，对得到的去噪后图像I₁和去噪后图像I₂进行运算，得到对数比值差异图X_l和均值比值差异图X_m，计算公式如下：

$X_l=\left|log\frac{I_2+1}{I_1+1}\right|=|log(I_2+1)-log(I_1+1)|;$

$X_m=1-min(\frac{{\mathrm{\μ}}_1}{{\mathrm{\μ}}_2},\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1});$

其中，μ₁，μ₂分别表示I₁和I₂的局部均值。

步骤3，分别提取得到的对数比值差异图X_l和均值比值差异图X_m中的每个像素点的l层邻域像素特征矩阵，其中l≥3,且l为奇数；然后将提取的邻域特征像素矩阵展开成列，得到该像素点特征描述列向量；再分别将对数比值差异图X_l和均值比值差异图X_m中所有像素点的特征描述列向量作为不同的列依次排列，得到对数比值差异图X_l和均值比值差异图X_m对应的特征描述矩阵X_L和特征描述矩阵X_M；本发明中使用每个像素点的l＝3邻域像素点以及本身来描述该像素点。并将该特征向量展开排列成一列，逐像素展开得到图像I₁和I₂对应的特征描述矩阵X_L和X_M。若原始图像大小为n×m，则对应特征描述矩阵X_L和X_M的大小均为l²×nm，其中nm为像素点个数。

步骤4，根据像素位置对应的原则，对得到的特征描述矩阵X_L和特征描述矩阵X_M按列拼接，得到变化检测特征矩阵：X的大小为2l²×nm。

步骤5，采用乘性迭代的分解方式，对得到的变化检测特征矩阵X进行NMF非负矩阵分解，分解公式为其中W为基矩阵，H为系数矩阵，r为分解矩阵的秩，分解步骤如下：

步骤5a)设矩阵分解的秩r＝1，并对基矩阵W和系数矩阵H进行随机初始化，其中基矩阵W≥0，系数矩阵H≥0；

步骤5b)设迭代次数Iterator的初始值为0，初始化迭代误差上限为ε，迭代次数上限为Iter；

步骤5c)验证迭代条件是否满足：误差E＝|X-WH|满足||E||≥ε并且Iterator≤Iter，若是，执行步骤5d)，否则，迭代结束，输出非负矩阵分解后的系数矩阵H，并执行步骤6；

步骤5d)针对变化检测特征矩阵X，其非负矩阵分解形式中基矩阵W和系数矩阵H的求解公式如下：

$W_{ia}\←W_{ia}\underset{\mathrm{\μ}}{\Σ}\frac{X_{i\mathrm{\μ}}}{{\left(WH\right)}_{i\mathrm{\μ}}}{H}_{a\mathrm{\μ}};$

$W_{ia}\←\frac{W_{ia}}{\underset{j}{\Σ}{W}_{ja}};$

$H_{a\mathrm{\μ}}\←H_{a\mathrm{\μ}}\underset{i}{\Σ}{W}_{ia}\frac{X_{i\mathrm{\μ}}}{{\left(WH\right)}_{i\mathrm{\μ}}};$

Iterator++；

E＝|X-WH|；

步骤6，将系数矩阵H按照原始像素排列顺序恢复为原始图像尺寸，得到融合后差异图X_FUSION，其大小分别为2l²×r和r×nm，将系数矩阵H按照原始像素排列顺序对应恢复成原始图像尺寸，得到大小为n×m的融合后差异图X_FUSION。

步骤7，根据k-means聚类算法，将得到的融合后差异图X_FUSION聚为变化和未变化两类，得到变化检测结果CM。该算法接受参数k，以空间中k个点为中心进行聚类，对最靠近他们的对象归类。通过迭代的方法，逐次更新各聚类中心的值，直至得到最好的聚类结果。假设要把样本集分为k个类别，算法描述如下：

7.1适当选择k个类的初始中心；

7.2在第i次迭代中，对任意一个样本，求其到k个中心的距离，将该样本归到距离最短的中心所在的类；

7.3利用均值等方法更新该类的中心值；

7.4对于所有的k个聚类中心，如果利用7.2、7.3的迭代法更新后，值保持不变，则迭代结束，否则继续迭代。

如原始数据的集合为(x₁,x₂,…,x_n)，并且每个x_i为d维的向量，k-means聚类的目的就是，在给定分类组数k的条件下，其中k≤n，将原始数据分成k类S＝{S₁,S₂,…,S_k}，其目标函数如下：

以下结合仿真实验，对本发明的技术效果作进一步说明：

1.仿真条件

本实例在AMD Athlon(tm)P340Dual-Core Processor 2.20GHz Windows 8.1系统下，Matlab 2014a运行平台上，完成本发明仿真实验。

2.仿真实验内容

1)本发明所使用的两个数据集参考图和对该数据集进行SAR图像变化检测的结果对比；两数据集结果分别如图2和图3所示。

2)本发明对两个数据集进行NMF融合，和现有技术对该数据集分别进行均值比运算、对数比运算和小波融合的差异对比图以及对应差异图的变化检测结果；如图4和图5所示。

本实验选择了Ottawa和Berne两个数据集做变化检测结果的验证，Ottawa数据集中图像的大小为290×350，Berne数据集图像的大小为301×301。

3.仿真实验结果及分析

参照图2，其中图2(a)是Ottawa数据集的参考图，图2(b)是本实施例得到的Ottawa数据集变化检测结果，可以发现本实施例得到的Ottawa数据集变化检测结果能很好的标注出变化和未变化部分；参照图3，其中图3(a)是Berne数据集的参考图，图3(b)是本实施例得到的Berne数据集变化检测结果，可以发现，出来少部分的漏检之外，本实施例得到的Berne数据集变化检测结果也较准确的标注了变化和未变化部分，实现了变化检测的目的。

参照图4，其中图4(a)-图4(d)分别表示Ottawa数据集对应的均值比值差异图、对数比值差异图、小波融合差异图以及NMF非负矩阵分解融合差异图，图4(e)-图4(h)分别表示上述四幅图差异图像对应的变化检测结果。对比图4(a)、图4(b)和图4(c)、图4(d)可以发现，图4(c)和图4(d)融合差异图对比图4(a)和图4(b)有效的抑制了图4(a)图中均值比值差异图的背景噪声，同时增强了图4(b)图中变化部分。从图4(e)-图4(h)的对比可以看出，图像融合有效的增强了变化检测的结果，同时对比图4(g)和图4(h)可以发现，图4(h)比图4(g)增加了更多的细节信息，变化检测结果更准确。

参照图5，其中图5(a)-图5(d)分别表示Berne数据集对应的均值比值差异图、对数比值差异图、小波融合差异图以及NMF非负矩阵分解融合差异图，图5(e)-图5(h)分别表示上述四幅图差异图像对应的变化检测结果。对比图5(a)、图5(b)和图5(c)、图5(d)可以发现，图5(c)和图5(d)融合差异图与图5(a)和图5(b)相比，有效的抑制了图5(a)图中均值比值差异图的背景噪声，同时增强了图5(b)图中变化部分。从图5(e)-图5(h)的对比可以看出，图像融合有效的增强了变化检测的结果，同时对比图5(g)和图5(h)可以发现，图5(h)进一步抑制了图5(g)中的背景噪声，同时增加了更多的变化区域细节信息,即NMF融合方法相比小波融合方法，能更有效的实现变化检测这一目的。

表1和表2分别给出了两个数据集上对比试验结果。进一步验证了我们的分析。

表1

Ottawa数据集上基于四种不同差异图和k-means聚类方法变化检测结果

差异图	虚警数	漏警数	准确率	一致性评价
					均值比差异图	2555	365	0.9712	0.8976
对数比差异图	930	1957	0.9716	0.8903
					小波融合差异图	949	894	0.9818	0.9319
NMF融合差异图	687	409	0.9891	0.9594

表2

Berne数据集上基于四种不同差异图和k-means聚类方法变化检测结果

差异图	虚警数	漏警数	准确率	一致性评价
					均值比差异图	15552	8	0.8283	0.107
对数比差异图	363	329	0.9924	0.703
					小波融合差异图	503	77	0.9936	0.784
NMF融合差异图	100	214	0.9950	0.8460

通过对比表1和表2可以发现，本发明中所采用方法的变化检测漏警数和虚警数相对基本的均值比差异图、对数比差异图都有着明显的降低，并且相对已有的基于小波融合的变化检测结果来讲，准确率和一致性评价参数也有着较大幅度的提升。

根据以上结果可知，本发明的方法操作更简单，对于SAR图像变化检测更有效。总体来说，本发明方法明显提高了SAR图像变化检测的准确率，可实现SAR图像变化信息的分析和标注。

Claims (3)

1.一种基于NMF图像融合的SAR图像变化检测方法，包括如下步骤：

(1)对同一地域不同时间的两幅SAR原始图像分别进行去噪处理，得到去噪后图像I₁和去噪后图像I₂；

(2)分别利用对数比算子和均值比算子，对得到的去噪后图像I₁和去噪后图像I₂进行运算，得到对数比值差异图X_l和均值比值差异图X_m，计算公式如下：

$X_l=\left|log\frac{I_2+1}{I_1+1}\right|=|log(I_2+1)-log(I_1+1)|;$

$X_m=1-min(\frac{{\mathrm{\μ}}_1}{{\mathrm{\μ}}_2},\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1});$

其中，μ₁，μ₂分别表示I₁和I₂的局部均值；

(3)对得到的对数比值差异图X_l和均值比值差异图X_m进行特征描述，得到对应的特征描述矩阵X_L和特征描述矩阵X_M，实现步骤如下：

3a)分别提取得到的对数比值差异图X_l和均值比值差异图X_m中每个像素点的l层邻域像素特征矩阵，其中l≥3，且l为奇数；

3b)将提取的l层邻域特征像素矩阵展开成一列，得到该像素点特征描述列向量；

3c)分别将对数比值差异图X_l和均值比值差异图X_m中所有像素点的特征描述列向量作为不同的列依次排列，得到的特征描述矩阵X_L和特征描述矩阵X_M。

(4)根据像素位置对应的原则，对得到的特征描述矩阵X_L和特征描述矩阵X_M按列拼接，得到变化检测特征矩阵：

(5)采用乘性迭代的分解方式，对得到的变化检测特征矩阵X进行NMF非负矩阵分解，该NMF非负矩阵分解的矩阵表达式为X_m×n≈(WH)_m×n＝W_m×rH_r×n，其中W为基矩阵，H为系数矩阵，r为矩阵分解的秩，分解步骤如下：

5a)设矩阵分解的秩r＝1，并对基矩阵W和系数矩阵H进行随机初始化，其中基矩阵W≥0，系数矩阵H≥0；

5b)设迭代次数Iterator的初始值为0，初始化迭代误差上限为ε，迭代次数上限为Iter；

5c)验证迭代条件是否满足：误差E＝|X-WH|满足||E||≥ε并且Iterator≤Iter，若是，执行步骤5d)，否则，迭代结束，输出非负矩阵分解后的系数矩阵H，并执行步骤(6)；

5d)对基矩阵W和系数矩阵H分别进行迭代，单次迭代完成执行步骤5c)，迭代按如下步骤实现：

$W_{ia}\←W_{ia}\underset{\mathrm{\μ}}{\Σ}\frac{X_{i\mathrm{\μ}}}{{\left(WH\right)}_{i\mathrm{\μ}}}{H}_{a\mathrm{\μ}};$

$W_{ia}\←\frac{W_{ia}}{\underset{j}{\Σ}{W}_{ja}};$

$H_{a\mathrm{\μ}}\←H_{a\mathrm{\μ}}\underset{i}{\Σ}{W}_{ia}\frac{X_{i\mathrm{\μ}}}{{\left(WH\right)}_{i\mathrm{\μ}}};$

Iterator++；

E＝|X-WH|；

(6)将输出的非负矩阵分解后的系数矩阵H，按照对数比值差异图X_l或均值比值差异图X_m中对应像素的排列顺序，恢复为原始图像的尺寸，得到融合后差异图X_FUSION；

(7)将得到的融合后差异图X_FUSION聚类成两类，得到同一地域不同时间的两幅SAR图像的变化检测结果图CM，并输出该结果图CM。

2.根据权利要求书1所述的基于NMF图像融合的SAR图像变化检测方法，其特征在于，步骤(1)中所述的去噪处理，采用均值滤波算法其中，(x,y)表示当前像素点，g(x,y)表示处理后图像在该点上的灰度值，m为该模板中包含当前像素在内的像素总个数，f(x,y)为原始图像中像素灰度值。

3.根据权利要求书1所述的基于NMF图像融合的SAR图像变化检测方法，其特征在于，步骤(7)中所述的聚类，采用k-means聚类算法。