CN102663740B - 基于图切的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于图切的SAR图像变化检测方法，主要解决现有技术中检测精度低的问题。其实现过程是：将变化前后的SAR图像用图切的方法进行分块预处理；然后用分块后的变化前后的SAR图像，分别构造变化前邻接图和变化后邻接图；连接两个邻接图的对应顶点，将两个邻接图融合得到变化前后SAR图像的二维图；最后用NJW谱聚类的方法对二维图进行切割，得到变化后图像各块的标签，将变化后图像中表示了变化区域的一类数据块提取出来，得到最终的变化检测结果图。本发明避免了已有技术中构造差异图对图像进行降唯，所造成的信息丢失和累积误差，提高了变化检测的检测精度。

Description

基于图切的SAR图像变化检测方法

技术领域

本发明涉及图像处理技术领域，特别是SAR图像的处理，主要可应用于生态和环境的监控以及自然灾害的评测与预防。

背景技术

由于合成口径雷达SAR不受大气条件和云层覆盖等相关条件的影响，因此，多时相SAR图像技术在人们的生产生活中发挥着越来越重要的作用。而SAR图像变化检测技术，在相关应用中正起着至关重要的作用，例如生态和环境的监控以及自然灾害的评测与预防，都在起着越来越重要的作用。图像变化检测就是研究同一场景不同时段的两幅或多幅图像之间发生的变化，不同时段同一场景的SAR图像的获得，使通过SAR图像变化检测技术提供目标与场景实时动态信息成为可能，使得SAR图像变化检测技术在实际应用中变得越来越重要。

在SAR图像中寻找“非变化”和“变化”两分类的研究目前还处于初步阶段，大致有两个路线，一条路线是分类后比较方法，也称后分类比较法，该方法先对两个时相的图像进行独立分类，再对两幅分类图像进行逐像素的比较，最后得到变化检测图；另外一条路线是差异图分类方法，该方法先对两个时相的图像进行逐像素的比较，比如逐像素点的差值、比值、CVA等，再对比较得出的差异图像进行进一步地多种变换、概率分布等处理以达到两分类，最后得到变化检测图。后分类比较法能够减少由于数据获取平台和环境的不同引起的伪变化信息，不需要数据的辐射校正等复杂的预处理，目前较多的研究是沿着差异分类模型这个路线进行的。差异图分类方法简单直观，得到的变化细节较为显著。

目前变化检测方法的缺陷主要有：对于差异图分类法而言，差异图可分性的高低直接影响着变化检测结果，而高可分性差异图的构建至今是一个难题，同时在构造差异图的过程中不可避免的造成了图像信息的丢失；对于分类后比较法而言，各时相分类后的效果对后续的分类结果起着比较至关重要的作用，而图像分类本身就是个比较困难的问题，分类的精度很难保证，在比较的过程中存在着误差累积的问题。

发明内容

本发明的目的在于针对上述已有技术的不足，提出了一种基于图切的SAR图像变化检测方法，以减小信息丢失和累积误差，提高SAR图像变化检测的检测精度。

为实现上述目的，本发明基于图切的SAR图像变化检测方法，包括如下步骤：

(1)输入变化前后的SAR图像P₁和P₂，并选用变化前后SAR图像的变化检测标准结果图q₀；

(2)将变化前后的SAR图像P₁和P₂进行分块预处理，得到m块以及预处理后的变化前图像I₁和变化后图像I₂，其中预处理后的变化前图像I₁中的块标记为d_1i，i＝1，2，…，m，预处理后的变化后后图像I₂中的块标记为d_2i，i＝1，2，…，m；

(3)用分块预处理后的变化前图像I₁构造变化前邻接图：G₁＝[V₁，S₁]，变化前邻接图G₁的顶点V₁表示分块预处理后的变化前图像I₁中的各个小块，S₁表示变化前邻接图G₁的加权边集合：S₁＝{s₁(i，j)}i，j＝1，2，…，m，

$s_1(i,j)=\exp (-\frac{{\left|\right|d_{1i}-d_{1j}\left|\right|}^2}{2{\mathrm{\σ}}^2}),$

其中，σ是高斯核的尺度参数，exp是自然指数；

(4)用分块预处理后的变化后图像I₂构造变化后邻接图：G₂＝[V₂，S₂]，变化后邻接图G₂的顶点V₂表示分块预处理后的变化后图像I₂中的各个小块，S₂表示变化后邻接图G₁的加权边集合：S₂＝{s₂(i，j)}i，j＝1，2，…，m，

$s_2(i,j)=\exp (-\frac{{\left|\right|d_{2i}-d_{2j}\left|\right|}^2}{2{\mathrm{\σ}}^2});$

(5)将变化前邻接图G₁和变化后邻接图G₂中对应顶点连接，得到融合图：G＝[V₀，S₀]，融合图G的顶点V₀表示分块预处理后的变化前图像I₁和变化后图像I₂中的各个小块，融合图G中的加权边的集合S₀＝S₁∪S₂∪V，V表示变化连接变化前邻接图G₁和变化后邻接图G₂对应顶点的加权边集合：

V＝{v(i)}i＝1，2，…，m，

$v\left(i\right)=1-\exp (-\frac{{\left|\right|d_{1i}-d_{2i}\left|\right|}^2}{2{\mathrm{\σ}}^2});$

(6)计算融合图G的差异矩阵W＝[w(i，j)]，i，j＝1，2，…，m，

$w(i,j)=\left\{\begin{array}{cc}|s_1(i,j)-s_2(i,j)|+v\left(i\right)+v\left(j\right)& \mathrm{ifi}\≠j\\ 3\·v\left(i\right)& \mathrm{otherwise}\end{array}\right.$

其中，v(i)表示分块预处理后的变化前图像I₁和变化后图像I₂中块d_1i和d_2i之间的不相似度，s₁(i，j)表示分块预处理后的变化前图像I₁中块d_1i与块d_1j之间的相似度，s₂(i，j)表示分块预处理后的变化后图像I₂中块d_2i与块d_2j之间的相似度；

(7)对差异矩阵W进行NJW谱聚类，类别数k＝3，得到差异矩阵W的第i行w_i的类别标签，i＝1，2，…，m；

(8)根据差异矩阵W中行w_i的类别标签，给对应块d_2i指派类别标签，将分块处理后的变化后图像I₂聚为k类，得到第h类的数据集合A_h，1≤h≤k；

(9)根据标准结果图q₀，找到聚类后的k类变化后图像I₂中代表变化区域的一类数据块A_h，将从分块处理后的变化后图像I₂中提取出来的A_h作为最终的变化检测结果图。

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

1.本发明以构造图的形式将变化前后的SAR图像之间的相关信息体现在一个二维图中，通过对图的切割和后续处理，完成了变化前后SAR图像中变化区域和非变化区域的划分，实现了图像的变化检测；

2.本发明克服了已有技术中对数据维数的限制，避免了已有技术中首先对图像进行降维所造成的信息丢失和累积误差问题，降低了错误率；

3.本发明中采用NJW谱聚类的方法实现对图的切割，使得本发明不会陷入局部最优以及较充分的利用了类内和类间信息，提高了变化检测的检测精度。

附图说明

图1是本发明的流程图；

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

图3是现有的对Feltwell SAR图像数据集变化检测的标准结果图；

图4是用本发明和传统变化检测方法对图2的变化检测结果图；

图5是本发明仿真所使用的撒丁岛SAR图像数据集；

图6是现有的对撒丁岛SAR图像数据集变化检测的标准结果图；

图7是用本发明和传统变化检测方法对图5的变化检测结果图。

具体实施方式

参照图1，本发明基于图切的SAR图像变化检测方法，包括如下步骤：

步骤1：输入变化前的SAR图像P₁和变化后的SAR图像P₂，以及选用变化前后SAR图像的变化检测标准结果图q₀，用图切的方法对变化前后的SAR图像P₁和p₂进行分块预处理。

1a)将变化前的SAR图像P₁和变化后的SAR图像P₂相减得到差值图：I＝|P₁-P₂|；

1b)用图切的方法将差值图I分割成若干小块，并分别记录各小块标记e，e∈{1，…，m}，m为块的个数；

1c)根据差值图像I的块标记，将变化前后的SAR图像P₁和P₂分割成与差值图像I具有相同划分的预处理后的变化前图像I₁和变化后图像I₂，分别得到块d_1i和d_2i，i＝1，2，…，m。

步骤2：分别用分块预处理后的变化前图像I₁和变化后图像I₂，构造变化前邻接图G₁和变化后邻接图G₂。

2a)用分块预处理后的变化前图像I₁构造变化前邻接图：G₁＝[V₁，S₁]，变化前邻接图G₁的顶点V₁表示分块预处理后的变化前图像I₁中的各个小块，S₁表示变化前邻接图G₁的加权边集合：S₁＝{s₁(i，j)}i，j＝1，2，…，m，

$s_1(i,j)=\exp (-\frac{{\left|\right|d_{1i}-d_{1j}\left|\right|}^2}{2{\mathrm{\σ}}^2}),$

其中，σ是高斯核的尺度参数，exp是自然指数。

2b)用分块预处理后的变化后图像I₂构造变化后邻接图：G₂＝[V₂，S₂]，变化后邻接图G₂的顶点V₂表示分块预处理后的变化后图像I₂中的各个小块，S₂表示变化后邻接图G₁的加权边集合：S₂＝{s₂(i，j)}i，j＝1，2，…，m，

$s_2(i,j)=\exp (-\frac{{\left|\right|d_{2i}-d_{2j}\left|\right|}^2}{2{\mathrm{\σ}}^2}).$

步骤3：连接变化前邻接图G₁和变化后邻接图G₂的对应顶点，将变化前邻接图G₁和变化后邻接图G₂融合，得到二维图：G＝[V₀，S₀]，二维图G的顶点V₀表示分块预处理后的变化前图像I₁和变化后图像I₂中的各个小块，二维图G中的加权边集合S₀＝S₁∪S₂∪V，V表示连接变化前邻接图G₁和变化后邻接图G₂的加权边集合：V＝{v(i)}i＝1，2，…，m，

$v\left(i\right)=1-\exp (-\frac{{\left|\right|d_{1i}-d_{2i}\left|\right|}^2}{2{\mathrm{\σ}}^2}).$

步骤4：计算二维图G的差异矩阵W＝[w(i，j)]，i，j＝1，2，…，m，

$w(i,j)=\left\{\begin{array}{cc}|s_1(i,j)-s_2(i,j)|+v\left(i\right)+v\left(j\right)& \mathrm{ifi}\≠j\\ 3\·v\left(i\right)& \mathrm{otherwise}\end{array}\right.,$

其中，v(i)表示分块预处理后的变化前图像I₁和变化后图像I₂中块d_1i和d_2i之间的不相似度，s₁(i，j)表示分块预处理后的变化前图像I₁中块d_1i与块d_1j之间的相似度，s₂(i，j)表示分块预处理后的变化后图像I₂中块d_2i与块d_2j之间的相似度。

步骤5：对差异矩阵W进行NJW谱聚类，类别数k＝3，得到差异矩阵W的行数为w_i的标签，i＝1，2，…，m。

5a)定义对角矩阵为D，将差异矩阵W的w_i行元素之和作为对角矩阵D对角上的元素： $D(i,i)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{w}_{\mathrm{ij}};$

5b)构造拉普拉斯矩阵：L＝D^-1/2WD^-1/2；

5c)对L进行SVD奇异值分解，得到它的前k个最大特征值所对应的特征向量v¹，v²，...，v^k，k＞1且k∈N并构造特征向量矩阵U＝[v¹，v²，...，v^k]，其中v_i，i＝1，2，…，k为列向量；

5d)归一化特征向量矩阵U的行向量，得到归一化后的矩阵Y，其中

将归一化后的矩阵Y的每一行看成是R^k空间内的一点，用模糊C均值的万法将归一化后的矩阵Y聚为k类。

步骤6：根据差异矩阵W的行数为w_i的标签，给对应块d_2i指派类别标签，将分块处理后的变化后图像I₂聚为k类，即如果差异矩阵W的w_i行为第h类，1≤h≤k且h∈N，就将原数据块d_2i也划分到第h类，得到第h类的数据集合A_h。

步骤7：将分块处理后的变化后图像I₂中k类的每一类数据块A_h分别从分块处理后的变化后图像I₂中提取出来对其赋值255，在剩余区域赋值0，得到与标准结果图大小相同的类别图像q_h，1≤h≤k，将k幅类别图像q₁，…，q_h，…，q_k分别与标准结果图q₀相减，得到k幅比较图像c₁，…，c_h，…c_k，

c_h＝|q_h-q₀|，

其中，标准结果图q₀中变化区域的灰度值为255，非变化区域的灰度值为0。

步骤8：在k幅比较图像c₁，…，c_h，…c_k中查找非零元素最少的比较图像c_h，将对应的类别图像q_h作为最终的变化检测结果图。

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

1.实验条件

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

仿真选用的第一个数据集为英格兰Feltwell村庄SAR图像数据集，如图2所示，其中图2(a)为英格兰Feltwell村庄变化前的SAR图像，图2(b)为英格兰Feltwell村庄变化后的SAR图像。该图像位于英国Feltwell村庄的一个农田区，模拟变化图像通过模拟地球的天气变化和电磁波的辐射特性等因素影响并人工地嵌入一些变化区域得到，图像大小均为470×335像素，灰度级为256，包括153214个非变化像素和4236个变化像素。检测的标准结果图采用如图3所示的对Feltwell SAR图像数据集变化检测的结果图。在实验过程中的高斯核的尺度参数σ选取的是0.5。

第二个数据集为意大利撒丁岛SAR图像数据集，如图5所示，其中图5(a)为意大利撒丁岛变化前的SAR图像，图5(b)为意大利撒丁岛变化后的SAR图像。该组真实遥感数据集由1995年9月和1996年7月Landsat-5卫星TM传感器接收的两幅多光谱图像构成。图像大小均为300×412像素，灰度级为256。包括115974个非变化像素和7626个变化像素。检测的标准结果图采用如图6所示的对撒丁岛SAR图像数据集变化检测的结果图。在实验过程中的高斯核的尺度参数σ选取的是0.3。

2.实验内容和实验结果

实验一：用本发明方法和两种变化检测传统方法：ID差值法和KI阈值分割法，对图2进行变化检测。实验结果如图4所示，其中4(a)为现有ID差值法对图2(b)进行变化检测的结果图，4(b)为现有KI阈值分割法对图2(b)进行变化检测的结果图，4(c)为本发明方法对图2(b)进行变化检测的结果图。

实验二：用本发明方法和两种变化检测传统方法：ID差值法和KI阈值分割法，对图5进行变化检测。实验结果如图7所示，其中7(a)为现有ID差值法对图5(b)进行变化检测的结果图，7(b)为现有KI阈值分割法对图5(b)进行变化检测的结果图，7(c)为本发明方法对图5(b)进行变化检测的结果图。

从图4(a)和图7(a)中可以看出，ID差值法对图像进行变化检测的实验结果中有很多的杂点；从图4(b)和图7(b)中可以看出，KI阈值分割法对图像进行变化检测的实验结果产生了很多的细节丢失，边缘很模糊，图4(b)丢失了大部分的变化区域；从图4(c)和图7(c)中可以看出，本发明方法能够在较好保持边缘细节的情况下，减少杂点的个数。

本发明方法和所述两种传统变化检测方法，对图2和图5进行变化检测的结果数据，如表1所示：

表1试验结果数据

表1中列出了三种评价指标：分别为虚警数FA，漏检数MA和总错误数OE，其中，虚警数FA为实际没有发生变化但被当作变化检测出来的像素，漏检数MA为没有检测出来的实际发生了变化的像素，总错误数OE＝FA+MA。

从表1可以看出，本发明与所述的两种传统变化检测方法相比，可以获得较少的虚警数FA和漏检数MA，以及最少的总错误数OE，提高了变化检测的检测精度。

Claims (4)

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

（1）输入变化前后的SAR图像P₁和P₂，并选用变化前后SAR图像的变化检测标准结果图q₀；

（2）将变化前后的SAR图像P₁和P₂进行分块预处理，得到m块以及预处理后的变化前图像I₁和变化后图像I₂，其中预处理后的变化前图像I₁中的块标记为d_1i，i＝1,2,...，m，预处理后的变化后图像I₂中的块标记为d_2i，i＝1,2,...,m；

（3）用分块预处理后的变化前图像I₁构造变化前邻接图：G₁=[V₁,S₁]，变化前邻接图G₁的顶点V₁表示分块预处理后的变化前图像I₁中的各个小块，S₁表示变化前邻接图G₁的加权边集合：S₁={s₁(i,j)}，i,j＝1,2,...,m，

$s_1(i,j)=\exp (-\frac{{\left|\right|d_{1i}-d_{1j}\left|\right|}^2}{2{\mathrm{\σ}}^2}),$

其中，σ是高斯核的尺度参数，exp是自然指数；

（4）用分块预处理后的变化后图像I₂构造变化后邻接图：G₂=[V₂,S₂]，变化后邻接图G₂的顶点V₂表示分块预处理后的变化后图像I₂中的各个小块，S₂表示变化后邻接图G₁的加权边集合：S₂={s₂(i,j)}，i,j＝1,2,...，m，

$s_2(i,j)=\exp (-\frac{{\left|\right|d_{2i}-d_{2j}\left|\right|}^2}{2{\mathrm{\σ}}^2});$

（5）将变化前邻接图G₁和变化后邻接图G₂中对应顶点连接，得到融合图：G=[V₀,S₀]，融合图G的顶点V₀表示分块预处理后的变化前图像I₁和变化后图像I₂中的各个小块，融合图G中的加权边的集合S₀=S₁∪S₂∪V，V表示变化前邻接图G₁和变化后邻接图G₂对应顶点的加权边集合：V={v(i)}，i＝1,2,...,m，

$v\left(i\right)=1-\exp (-\frac{{\left|\right|d_{1i}-d_{2i}\left|\right|}^2}{{2\mathrm{\σ}}^2});$

（6）计算融合图G的差异矩阵W=[w(i,j)]，i,j＝1,2,...,m，

$w(i,j)=\left\{\begin{array}{cc}|s_1(i,j)-s_2(i,j)|+v\left(i\right)+v\left(j\right)& \mathrm{ifi}\≠j\\ 3\·v\left(i\right)& \mathrm{otherwise}\end{array}\right.$

其中，v(i)表示分块预处理后的变化前图像I₁和变化后图像I₂中块d_1i和d_2i之间的不相似度，s₁(i,j)表示分块预处理后的变化前图像I₁中块d_1i与块d_1j之间的相似度，s₂(i,j)表示分块预处理后的变化后图像I₂中块d_2i与块d_2j之间的相似度；

（7）对差异矩阵W进行NJW谱聚类，类别数k=3，得到差异矩阵W的第i行w_i的类别标签，i＝1,2,...,m；

（8）根据差异矩阵W中行w_i的类别标签，给对应块d_2i指派类别标签，将分块处理后的变化后图像I₂聚为k类，得到第h类的数据集合A_h，1≤h≤k；

（9）根据标准结果图q₀，找到聚类后的k类变化后图像I₂中代表变化区域的一类数据块A_h，得到最终的变化检测结果。

2.根据权利要求1所述的方法，其中步骤（2）所述的将变化前后的SAR图像P₁和P₂进行分块预处理，按如下步骤进行：

2a)将变化前的SAR图像P₁和变化后的SAR图像P₂相减得到差值图像I;

2b)用图切的方法将差值图像I分割成若干小块e,e∈{1,...,m}，m为块的个数；

2c)根据差值图像I的块标记，将变化前后的SAR图像P₁和P₂分割成与差值图像I具有相同划分的预处理后的变化前图像I₁和变化后图像I₂，分别得到块d_1i和d_2i，i＝1,2,...,m。

3.根据权利要求1所述的方法，其中步骤（7）所述的对差异矩阵W进行NJW谱聚类，按如下步骤进行：

7a)定义对角矩阵为D，将差异矩阵W的w_i行元素之和作为对角矩阵D对角上的元素： $D(i,i)=\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{w}_{\mathrm{ij}};$

7b)构造拉普拉斯矩阵：L＝D^-1/2WD^-1/2；

7c)对拉普拉斯矩阵L进行SVD奇异值分解，得到它的前k个最大特征值所对应的特征向量v₁,v₂,...,v_k，k＞1且k∈N，并构造特征向量矩阵：U＝[v₁v₂...v_k]∈R^n×k，其中v_i,i＝1,2,...,k为列向量；

7d)归一化特征向量矩阵U的行向量，得到归一化后的矩阵Y，将归一化后的矩阵Y的每一行看成是R^k空间内的一点,用模糊C均值FCM的方法将归一化后的矩阵Y聚为k类。

4.根据权利要求1所述的方法，其中步骤（9）所述的根据标准结果图q₀，找到聚类后的k类变化后图像I₂中表示变化区域的一类数据块A_h，按如下步骤进行：

9a)将分块处理后的变化后图像I₂中k类的每一类数据块A_h分别从分块处理后的变化后图像I₂中提取出来对其赋值255，在剩余区域赋值0，得到与标准结果图大小相同的类别图像q_h，1≤h≤k，将k幅类别图像q1,...,q_h,...,q_k分别与标准结果图q₀相减，得到k幅比较图像c₁,...,c_h,...c_k，

c_h=|q_h-q₀|，

其中，标准结果图q₀中变化区域的灰度值为255，非变化区域的灰度值为0；

9b)在k幅比较图像c₁,...,c_h,...c_k中查找非零元素最少的比较图像c_h，将对应的类别图像q_h作为最终的变化检测结果图。

Images