基于nmf图像融合的sar图像变化检测方法
【专利摘要】本发明提出了一种基于NMF图像融合的SAR图像变化检测方法，用于解决现有的基于图像融合的SAR图像变化检测方法中存在的准确率低且物理意义不明确的技术问题，实现步骤如下：输入同地域不同时间的两幅原始SAR图像，经去噪处理后构造两基本差异图Xl和Xm；分别对Xl和Xm进行邻域特征描述得到特征描述矩阵XL和XM；利用两特征描述矩阵构造变化检测特征矩阵X；对X进行NMF非负矩阵分解，得到分解后的系数矩阵H；利用系数矩阵得到融合后差异图；利用聚类算法将融合后差异图聚成两类，得到变化检测结果。本发明更好的利用了原始图像的细节信息，且该图像融合方法具有更显著的物理意义，可用于SAR图像变化信息的分析和标注。
【专利说明】
基于NMF图像融合的SAR图像变化检测方法
技术领域
[0001] 本发明属于图像处理技术领域，涉及一种SAR图像变化检测方法，具体涉及一种基 于NMF图像融合的SAR图像变化检测方法，可用于SAR图像变化信息的分析和标注。
【背景技术】
[0002] 随着一系列光学和雷达卫星的运行服务的发展，积累了大量的遥感数据，通过对 这些图像数据进行进一步的分析、处理及解译，就可以将经遥感技术获得的图像数据转化 为有用的信息，进而服务于国民经济的各个领域。
[0003] 凡是记录各种地物电磁波大小的图片都称为遥感图像，其中合成孔径雷达SAR (Synthetic Aperture Radar)图像是遥感图像的一个重要分支，它是一种先进的主动式微 波对地观测设备，具有全天候、全天时工作的特点，对地面植被有一定穿透能力，能获得类 似光学照片的目标图像，在国民经济和军事应用领域有着十分重要的作用。
[0004] SAR图像变化检测技术是遥感图像研究的重要组成部分，它是对同一地域不同时 刻拍摄的多幅遥感图像，采用图像处理和模式识别等手段进行对比分析，根据图像之间的 差异来得到人们所需要的地物或目标随时间发生的变化信息。可以应用于军事研究领域、 遥感气象服务、灾害动态监测、环境污染检测以及城市规划等方面。
[0005] Singh.A将图像变化检测方法分为分类比较法和直接比较法，分类比较法的变化 检测结果的精度容易受分类精度的影响，初始分类的错误将造成难以弥补的误检或漏检现 象;直接比较法的算法相对简单，容易实现，但难以克服大气条件和传感器噪声等对检测结 果带来的误差，受图像的预处理精度的影响严重。图像融合能有效的解决该问题，图像融合 是指将多源信道所采集到的关于同一目标的图像数据经过图像处理和计算机技术等，最大 限度的提取各自信道中的有利信息，最后综合成高质量的图像，以提高图像信息的利用率、 改善计算机解译精度和可靠性、提升原始图像的空间分辨率和光谱分辨率，利于监测。如公 茂果等2012年在国际电子与电气工程师协会地球科学与遥感杂志上发表的名为基于小波 融合的SAR图像变化检测方法文章中，发表了一种基于三级小波融合的SAR图像变化检测方 法，该方法是将对数比值差异图和均值比值差异图通过小波变换和逆变换得到融合差异 图，融合规则是基于低频信息的加权平均和高频信息的最小标准差。该文章采用了三级小 波变换和逆变换，得到小波融合差异图，最终实现对融合差异图的变化信息的分析和标记。 该方法为图像融合提供了新的思路，但物理意义不够明确，并且三级小波运算及逆运算导 致运算效率并不高。
[0006] 随着图像质量的增强，以及图像数量的增多，行之有效的处理大量遥感影像数据 成为亟待解决的问题。矩阵分解在降维、低秩和稀疏方面有着重要应用，且考虑图像数据本 身的特点，为变化检测问题的解决提供了更清晰的物理意义。吸引了很大一批学者的兴趣。
【发明内容】
[0007] 本发明的目的在于克服上述现有技术存在的缺陷，提供了一种基于NMF图像融合 的SAR图像变化检测方法，用于解决现有的基于图像融合的SAR图像变化检测方法中存在的 准确率低且物理意义不明确的技术问题。
[0008] 本发明的技术关键是基于像素的原始图像特征描述与基于非负矩阵分解的图像 融合方法。特征描述过程充分考虑了SAR图像的成像特点，使用像素的邻域信息来表征该像 素点。同时图像融合过程利用NMF非负矩阵分解融合方法，有效利用了基本差异图像的互补 信息，在有效抑制图像背景信息的同时大大增强感兴趣区域，能有效的提高检测准确率。 NMF非负矩阵分解方法充分考虑了物理图像信息的非负特点，物理意义更明确。并且分解过 程能实现矩阵的降维操作，提高了效率。
[0009] 根据上述技术思路，实现本发明目的采用的技术方案，包括如下步骤：
[0010] 1.对同一地域不同时间的两幅SAR原始图像分别进行去噪处理，得到去噪后图像 II和去噪后图像工2 ;
[0011] 2.分别利用对数比算子和均值比算子，对得到的去噪后图像h和去噪后图像12进 行运算，得到对数比值差异图心和均值比值差异图X m，计算公式如下：
[0014] 其中，μι，μ2分别表不Ιι和12的局部均值；
[0015] 3.对得到的对数比值差异图Xi和均值比值差异图Xm进行特征描述，得到对应的特 征描述矩阵Xl和特征描述矩阵Xm，实现步骤如下：
[0016] 3a)分别提取得到的对数比值差异图Xi和均值比值差异图Xm*每个像素点的1层邻 域像素特征矩阵，其中1多3，且1为奇数；
[0017] 3b)将提取的1层邻域特征像素矩阵展开成一列，得到该像素点特征描述列向量；
[0018] 3c)分别将对数比值差异图Xi和均值比值差异图Xm*所有像素点的特征描述列向 量作为不同的列依次排列，得到的特征描述矩阵Xl和特征描述矩阵Xm。
[0019] 4.根据像素位置对应的原则，对得到的特征描述矩阵4和特征描述矩阵XM按列拼 接，得到变化检测特征矩阵I;
[0020] 5.采用乘性迭代的分解方式，对得到的变化检测特征矩阵X进行NMF非负矩阵分 解，该NMF非负矩阵分解的矩阵表达式为X mXn~(WH)mXn = WmXrHrXn，其中W为基矩阵，Η为系数 矩阵，r为矩阵分解的秩，分解步骤如下：
[0021] 5a)设矩阵分解的秩r=l，并对基矩阵W和系数矩阵Η进行随机初始化，其中基矩阵 W彡0,系数矩阵Η彡0;
[0022] 5b)设迭代次数Iterator的初始值为0,初始化迭代误差上限为ε，迭代次数上限为 Iter；
[0023] 5c)验证迭代条件是否满足：误差E= |X_WH|满足| |E| |彡ε并且Iterator彡Iter， 若是，执行步骤5d)，否则，迭代结束，输出非负矩阵分解后的系数矩阵H，并执行步骤(6);
[0024] 5d)对基矩阵W和系数矩阵Η分别进行迭代，单次迭代完成执行步骤5c)，迭代按如 下步骤实现： LU〇28」 lterator++;
[0029] E= |X-ffH| ；
[0030] 6.将输出的非负矩阵分解后的系数矩阵H，按照对数比值差异图Xi或均值比值差 异图对应像素的排列顺序，恢复为原始图像的尺寸，得到融合后差异图X FUSI〇N;
[0031] 7.将得到的融合后差异图XFUSIQN聚类成两类，得到同一地域不同时间的两幅SAR图 像的变化检测结果图CM，并输出该结果图CM。
[0032] 本发明与现有技术相比，具有如下优点：
[0033] 1.本发明由于在获取融合图像时，采用乘性迭代的分解方式，对得到的变化检测 特征矩阵X进行NMF非负矩阵分解，充分考虑了物理图像信息的非负特点，同时在获取融合 图像之前增加了对像素特征描述的步骤，能更好的利用原始图像的信息，与现有技术采用 小波融合进行图像融合的方法相比，有效地提高了SAR图像变化检测的准确率，且物理意义 明确。
[0034] 2.本发明利用NMF非负矩阵分解方法进行图像融合的同时，实现了数据的降维操 作，与现有技术中采用的将数据降维和其他操作分开处理的方法相比，简化了运算步骤，提 高了运算效率。
【附图说明】
[0035] 图1为本发明的实现流程图；
[0036]图2为本发明所使用Ottawa数据集参考图和对该数据集进行SAR图像变化检测的 结果对比图；
[0037]图3为本发明所使用Berne数据集参考图和对该数据集进行SAR图像变化检测的结 果对比图；
[0038]图4为本发明对Ottawa数据集进行NMF融合，和现有技术对该数据集分别进行均值 比运算、对数比运算和小波融合的差异对比图以及对应差异图的变化检测结果；
[0039]图5为本发明对Berne数据集进行匪F融合，和现有技术对该数据集分别进行均值 比运算、对数比运算和小波融合的差异对比图以及对应差异图的变化检测结果；
【具体实施方式】
[0040]以下结合附图和实施例，对本发明作进一步详细说明。
[0041 ]参照图1，本发明的实现步骤如下：
[0042] 步骤1，获取同一地域不同时间的两幅相同大小的SAR图像，对其进行去噪处理，常 见的去噪处理方法有PPB去斑去噪算法、均值滤波算法和非局部均值算法。本实施例采用均 值滤波算法，均值滤波是典型的线性滤波算法，其采用的主要方法为邻域平均法。考虑到 SAR图像成像原理:各像素均受到邻域像素点的散射影响，可看作是其邻域像素点在该点一 定程度的叠加，因此对待处理的当前像素点，采用其邻域均值来代替该像素点，并且均值滤 波算法实现简单，因此本实施例中采用均值滤波算法对两幅原始SAR图像分别进行去噪。均 值滤波的基本原理是用均值代替原图像中的各个像素值，即对待处理的当前像素点(x，y)， 选择一个模板，该模板由其近邻的若干像素组成，求模板中所有像素的均值，再把该均值赋 予当前像素点(x，y)，作为处理后图像在该点上的灰度值g(x，y)，即
为该模板中包含当前像素在内的像素总个数，f(x，y)为原始图像中像素灰度值，去噪后图 像记为Ιι和12。
[0043]步骤2,分别利用对数比算子和均值比算子，对得到的去噪后图像和去噪后图像 12进行运算，得到对数比值差异图h和均值比值差异图X?，计算公式如下：
[0046] 其中，μι，μ2分别表不Ιι和12的局部均值。
[0047] 步骤3,分别提取得到的对数比值差异图h和均值比值差异图Xm中的每个像素点的 1层邻域像素特征矩阵，其中1多3,且1为奇数;然后将提取的邻域特征像素矩阵展开成列， 得到该像素点特征描述列向量;再分别将对数比值差异图Χι和均值比值差异图1中所有像 素点的特征描述列向量作为不同的列依次排列，得到对数比值差异图Χι和均值比值差异图 X m对应的特征描述矩阵4和特征描述矩阵Xm;本发明中使用每个像素点的1 = 3邻域像素点 以及本身来描述该像素点。并将该特征向量展开排列成一列，逐像素展开得到图像IdPI2对 应的特征描述矩阵Xl和Xm。若原始图像大小为nXm，则对应特征描述矩阵4和乂(《的大小均为 1 2 X nm，其中nm为像素点个数。
[0048] 步骤4,根据像素位置对应的原则，对得到的特征描述矩阵4和特征描述矩阵XM按 列拼接，得到变化检测特征矩阵，X的大小为21 2Xnm。 _Λ Μ _
[0049] 步骤5,采用乘性迭代的分解方式，对得到的变化检测特征矩阵X进行NMF非负矩阵 分解，分解公式为胃= ?χ~，其中W为基矩阵，Η为系数矩阵，r为分解矩阵的 a -1 秩，分解步骤如下：
[0050] 步骤5a)设矩阵分解的秩r=l，并对基矩阵W和系数矩阵Η进行随机初始化，其中基 矩阵W>0,系数矩阵Η>0;
[0051] 步骤5b)设迭代次数Iterator的初始值为0,初始化迭代误差上限为ε，迭代次数上 限为Iter;
[0052] 步骤5c)验证迭代条件是否满足：误差E= |X_WH|满足| |E| |彡ε并且Iterator彡 Iter，若是，执行步骤5d)，否则，迭代结束，输出非负矩阵分解后的系数矩阵H，并执行步骤 6；
[0053]步骤5d)针对变化检测特征矩阵X，其非负矩阵分解形式中基矩阵W和系数矩阵Η的 求解公式如下：
[0058] E= |X-ffH| ；
[0059] 步骤6,将系数矩阵Η按照原始像素排列顺序恢复为原始图像尺寸，得到融合后差 异图Xfusiqn，其大小分别为212Xr和rXnm，将系数矩阵Η按照原始像素排列顺序对应恢复成 原始图像尺寸，得到大小为η X m的融合后差异图XFUSI0N。
[0060] 步骤7，根据k-means聚类算法，将得到的融合后差异图XFUSIQN聚为变化和未变化两 类，得到变化检测结果CM。该算法接受参数k，以空间中k个点为中心进行聚类，对最靠近他 们的对象归类。通过迭代的方法，逐次更新各聚类中心的值，直至得到最好的聚类结果。假 设要把样本集分为k个类别，算法描述如下：
[0061] 7.1适当选择k个类的初始中心；
[0062] 7.2在第i次迭代中，对任意一个样本，求其到k个中心的距离，将该样本归到距离 最短的中心所在的类；
[0063] 7.3利用均值等方法更新该类的中心值；
[0064] 7.4对于所有的k个聚类中心，如果利用7.2、7.3的迭代法更新后，值保持不变，则 迭代结束，否则继续迭代。
[0065] 如原始数据的集合为(X1，X2，···，Xn)，并且每个Xi为d维的向量，k-means聚类的目 的就是，在给定分类组数k的条件下，其中k彡η，将原始数据分成k类S = {Si，&，…，Sk}，其目 标函数如下：
[0066] 以下结合仿真实验，对本发明的技术效果作进一步说明：
[0067] 1.仿真条件
[0068]本实例在AMD Athlon(tm)P340Dual_Core Processor 2 · 20GHz Windows 8 · 1系统 下，Mat lab 2014a运行平台上，完成本发明仿真实验。
[0069] 2.仿真实验内容
[0070] 1)本发明所使用的两个数据集参考图和对该数据集进行SAR图像变化检测的结果 对比;两数据集结果分别如图2和图3所示。
[0071] 2)本发明对两个数据集进行NMF融合，和现有技术对该数据集分别进行均值比运 算、对数比运算和小波融合的差异对比图以及对应差异图的变化检测结果；如图4和图5所 不。
[0072] 本实验选择了 Ottawa和Berne两个数据集做变化检测结果的验证，Ottawa数据集 中图像的大小为290 X 350，Berne数据集图像的大小为301 X 301。
[0073] 3.仿真实验结果及分析
[0074]参照图2，其中图2(a)是Ottawa数据集的参考图，图2(b)是本实施例得到的Ottawa 数据集变化检测结果，可以发现本实施例得到的Ottawa数据集变化检测结果能很好的标注 出变化和未变化部分;参照图3,其中图3(a)是Berne数据集的参考图，图3(b)是本实施例得 到的Berne数据集变化检测结果，可以发现，出来少部分的漏检之外，本实施例得到的Berne 数据集变化检测结果也较准确的标注了变化和未变化部分，实现了变化检测的目的。
[0075]参照图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)增加了更多的细节信息，变化检测结果更准确。
[0076]参照图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融合方法相 比小波融合方法，能更有效的实现变化检测这一目的。
[0077]表1和表2分别给出了两个数据集上对比试验结果。进一步验证了我们的分析。
[0078] 表 1
[0079] Ottawa数据集上基于四种不同差异图和k-means聚类方法变化检测结果
[0082] Berne数据集上基于四种不同差异图和k-means聚类方法变化检测结果
[0084]通过对比表1和表2可以发现，本发明中所采用方法的变化检测漏警数和虚警数相 对基本的均值比差异图、对数比差异图都有着明显的降低，并且相对已有的基于小波融合 的变化检测结果来讲，准确率和一致性评价参数也有着较大幅度的提升。
[0085]根据以上结果可知，本发明的方法操作更简单，对于SAR图像变化检测更有效。总 体来说，本发明方法明显提高了 SAR图像变化检测的准确率，可实现SAR图像变化信息的分 析和标注。
【主权项】
1. 一种基于NMF图像融合的SAR图像变化检测方法，包括如下步骤： (1) 对同一地域不同时间的两幅SAR原始图像分别进行去噪处理，得到去噪后图像Ii和 去噪后图像12; (2) 分别利用对数比算子和均值比算子，对得到的去噪后图像Ii和去噪后图像12进行运 算，得到对数比值差异图Xl和均值比值差异图Xm，计算公式如下：其中，μι，μ2分别表不Ii和12的局部均值； (3) 对得到的对数比值差异图Xl和均值比值差异图Xm进行特征描述，得到对应的特征描 述矩阵Xl和特征描述矩阵Xm，实现步骤如下： 3a)分别提取得到的对数比值差异图Xl和均值比值差异图Xm中每个像素点的1层邻域像 素特征矩阵，其中1 > 3，且1为奇数； 3b)将提取的1层邻域特征像素矩阵展开成一列，得到该像素点特征描述列向量； 3c)分别将对数比值差异图Xl和均值比值差异图Xm中所有像素点的特征描述列向量作 为不同的列依次排列，得到的特征描述矩阵Xl和特征描述矩阵Xm。 (4) 根据像素位置对应的原则，对得到的特征描述矩阵Xl和特征描述矩阵Xm按列拼接， 得到变化检测特征矩阵：(5) 采用乘性迭代的分解方式，对得到的变化检测特征矩阵X进行NM巧自负矩阵分解，该 非负矩阵分解的矩阵表达式为XmXn>(WH)mXn = WmXrHrXn，其中W为基矩阵，Η为系数矩阵，r 为矩阵分解的秩，分解步骤如下： 5a)设矩阵分解的秩r=l，并对基矩阵W和系数矩阵Η进行随机初始化，其中基矩阵 0，系数矩阵11>0; 5b)设迭代次数Iterator的初始值为0,初始化迭代误差上限为ε，迭代次数上限为 Iter； 5c)验证迭代条件是否满足:误差E= |X-WH|满足I |E| I 并且Iterator^Iter，若是， 执行步骤5d)，否则，迭代结束，输出非负矩阵分解后的系数矩阵H，并执行步骤(6); 5d)对基矩阵W和系数矩阵Η分别进行迭代，单次迭代完成执行步骤5c)，迭代按如下步 骤实现：Iter3tor++; E= |X-WH| ； (6) 将输出的非负矩阵分解后的系数矩阵H，按照对数比值差异图Xi或均值比值差异图 Xm中对应像素的排列顺序，恢复为原始图像的尺寸，得到融合后差异图Xfusion; (7) 将得到的融合后差异图化USI0N聚类成两类，得到同一地域不同时间的两幅SAR图像 的变化检测结果图CM，并输出该结果图CM。2. 根据权利要求书1所述的基于NMF图像融合的SAR图像变化检测方法，其特征在于，步 骤(1)中所述的去噪处理，采用均值滤波算其中，（x，y)表示当前像 素点，g(x，y)表示处理后图像在该点上的灰度值，m为该模板中包含当前像素在内的像素总 个数，f(x，y)为原始图像中像素灰度值。3. 根据权利要求书1所述的基于NMF图像融合的SAR图像变化检测方法，其特征在于，步 骤(7)中所述的聚类，采用k-means聚类算法。
【文档编号】G06T5/50GK106097290SQ201610390224
【公开日】2016年11月9日
【申请日】2016年6月3日 公开号201610390224.7, CN 106097290 A, CN 106097290A, CN 201610390224, CN-A-106097290, CN106097290 A, CN106097290A, CN201610390224, CN201610390224.7
【发明人】张建龙, 李巧, 高新波, 周晓鹏
【申请人】西安电子科技大学