CN102629377A - 基于显著性度量的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于显著性度量的遥感图像变化检测方法，主要解决现有变化检测方法容易受到噪声干扰影响导致检测结果中虚警较高的问题。其实现过程是：输入两时相遥感图像，首先对每幅图像分别进行中值滤波，得到两时相滤波后图像；然后分别计算两时相滤波后图像的灰度差异图和基于Chi-square距离的差异图，并对灰度差异图进行显著性度量，得到一幅显著图；接着分别对基于chi-square距离的差异图和显著图进行模糊c均值分割，最后对显著图的分割结果进行修正，生成最终的变化检测结果。实验结果表明，本发明能够有效抑制噪声的干扰，提高检测结果的精度，可用于对自然灾害的分析、城市规划建设等领域。

Description

基于显著性度量的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，具体来说是一种基于显著性度量的遥感图像变化检测方法，可应用于自然灾害检测、城市规划建设、道路监测等诸多领域。

背景技术

遥感图像的变化检测是指利用多时相获取的、覆盖同一地表区域的遥感图像来确定和分析地表变化。由于遥感技术是变化检测最主要的技术手段之一，变化检测成为目前遥感应用方法研究中的热点问题。它被广泛于应用于自然灾害评估、城市规划建设、道路检测、水系检测及环境保护等诸多领域。

目前出现的变化检测算法中，利用不同时相图像的特征进行差值比较的方法因其操作简单，且能取得较好的检测结果得到了广泛研究。在进行差值计算之前，有时需要对多时相遥感影像进行精确配准，否则会导致许多虚假变化信息的生成，因此，徐丽华等提出一种顾及配准误差的遥感影像变化检测方法，该方法能在一定程度上减轻阴影和投影差对变化检测结果的影响，变化区域的边缘保持较为理想，但虚警率较高，变化区域周围存在较多孤立的伪变化信息。Sicong Liu等提出一种基于光谱特征差异的变化检测方法，分别对提取的五组光谱特征差异图进行支撑矢量机分类，最后对各组结果进行融合得到最终结果，该方法是一种有监督的检测方法，在进行支撑矢量机分类时需要通过实地考察和目测来选取训练样本，这在实际中往往难以做到。MinghuiTian等提出一种基于视觉注意模型的遥感图像变化检测方法，首先分别提取两时相图的强度、纹理、方向等多种底层图像特征，对应差值得到强度、纹理、方向差异图，接着计算各个特征通道的显著图，最后通过线性组合的方法将其整合成最终的显著图，得到变化检测的结果。该方法对低对比度的图像效果较好，但是当图像的对比度较高时，检测结果不理想，虚警率较高。

发明内容

本发明的目的在于针对上述已有的遥感图像变化检测技术的不足，提出一种基于显著性度量的遥感图像变化检测方法，以使得对低对比度和高对比度的图像均能取得较好的检测结果，并在有效抑制虚警的同时提高检测结果的整体精度。

为实现上述目的，本发明的检测方法包括如下步骤：

(1)分别对输入的两幅已配准的大小均为m×n的两时相遥感图像I¹和I²进行中值滤波，得到两时相遥感图像I¹和I²的滤波后图像X¹和X²；

(2)对滤波后图像X¹和X²作差得到灰度差异图D1；

(3)计算滤波后图像X¹和X²基于Chi-square距离的差异图D2：

$D2=\left[\begin{array}{cccc}D2\left(\mathrm{1,1}\right)& D2\left(\mathrm{1,2}\right)& ...& D2(1,n)\\ D2\left(\mathrm{2,1}\right)& ...& ...& ...\\ ...& ...& D2(i,j)& ...\\ D2(m,1)& D2(m,2)& ...& D2(m,n)\end{array}\right]$

$D2(i,j)=\underset{k=1}{\overset{n1}{\mathrm{\Σ}}}{\left(\frac{X^1(i,j)-X^2(i,j)}{{\mathrm{\σ}}_{D1}}\right)}^2$

式中i，j分别表示当前像素点的坐标索引，σ表示标准差，∑表示求和，n1表示以当前像素点为中心，所取滑动窗口内像素点的总个数，本发明中取n1＝9，即滑动窗口大小为3×3；

(4)对灰度差异图D1进行显著性度量，得到显著图D3；

(5)对差异图D2和显著图D3分别进行模糊c均值聚类分割，得到两组分割结果，依次表示为Z1和Z2，其中分割结果图Z1和初始变化检测结果图Z2中被标记为变化类的像素点赋值为1，被标记为非变化类的像素点赋值为0；

(6)利用分割结果图Z1对初始变化检测结果图Z2依据不同情况进行修正，去除分割结果图Z2中的孤立小噪声点，生成变化检测结果图Z。

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

(1)本发明由于对灰度差异图进行显著性度量，将得到的显著图作为差异影像进行分割，可以改善灰度差异图中变化区域周围干扰较多的情况，能有效抑制虚警；

(2)本发明由于对初始变化检测结果进行修正，能进一步减少伪变化信息，对低对比度和高对比度的图像均能取得较好的检测效果，提高了检测结果的整体精度。

附图说明

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

图2是本发明使用的两时相遥感图像及对应参考图；

图3是用本发明及对比方法对模拟遥感数据集实验得到的变化检测结果图；

图4是用本发明及对比方法对第一组真实遥感数据集实验得到的变化检测结果图；

图5是用本发明及对比方法对第二组真实遥感数据集实验得到的变化检测结果图。

具体实施方式

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

步骤1，对输入的两幅已配准且大小均为m×n的两时相遥感图像I¹和I²，如2(a)和图2(b)所示，分别进行窗口大小为5×5的中值滤波，得到滤波后的图像X¹和X²。

步骤2，对滤波后的图像X¹和X²作差得到灰度差异图D1；

D1＝|X¹-X²|。

步骤3，计算滤波后图像X¹和X²基于Chi-square距离的差异图D2：

$D2=\left[\begin{array}{cccc}D2\left(\mathrm{1,1}\right)& D2\left(\mathrm{1,2}\right)& ...& D2(1,n)\\ D2\left(\mathrm{2,1}\right)& ...& ...& ...\\ ...& ...& D2(i,j)& ...\\ D2(m,1)& D2(m,2)& ...& D2(m,n)\end{array}\right]$

$D2(i,j)=\underset{k=1}{\overset{n1}{\mathrm{\Σ}}}{\left(\frac{X^1(i,j)-X^2(i,j)}{{\mathrm{\σ}}_{D1}}\right)}^2$

其中σ表示标准差，∑表示求和，n1表示以当前像素点为中心，所取滑动窗口内像素点的总个数，n1＝9，即滑动窗口大小为3×3，i，j分别为当前像素点的坐标索引，即第i行第j列，i的取值范围为[1，m]，j的取值范围为[1，n]，m，n为图像大小。

步骤4，对灰度差异图D1进行显著性度量，得到显著图D3，具体按照以下步骤进行：

(4a)将灰度差异图D1作为初始图像DS⁰，即DS⁰＝D1；

(4b)当迭代次数h＝1时，对初始图像DS⁰中的每个像素点按照下式计算，得到第一次迭代后的图像DS¹＝DS^h：

$D{S}^h=\left[\begin{array}{cccc}D{S}^h\left(\mathrm{1,1}\right)& D{S}^h\left(\mathrm{1,2}\right)& ...& D{S}^h(1,n)\\ D{S}^h\left(\mathrm{2,1}\right)& ...& ...& ...\\ ...& ...& D{S}^h(i,j)& ...\\ D{S}^h(m,1)& D{S}^h(m,2)& ...& D{S}^h(m,n)\end{array}\right]$

$D{S}^h(i,j)=\left[\underset{u=1}{\overset{m}{\mathrm{\Σ}}}\underset{v=1}{\overset{n}{\mathrm{\Σ}}}\right|D{S}^0(i,j)-D{S}^0(u,v)\left|\right]/(m\×n)$

式中i和j分别表示当前像素点的坐标索引，即第i行，第j列，i的取值范围为[1，m]，j的取值范围为[1，n]，m，n表示图像大小，u和v分别为图像中除当前像素点外其它像素点的坐标索引，即第u行，第v列，DS⁰(i，j)表示初始图像DS⁰中坐标为(i，j)的元素值，DS⁰(u，v)表示初始图像DS⁰中坐标为(u，v)的元素值；

(4c)当迭代次数h＝2时，对第一次迭代后的图像DS¹中的每个像素点按照下式计算，得到第二次迭代后的图像DS²＝DS^h：

$D{S}^h=\left[\begin{array}{cccc}D{S}^h\left(\mathrm{1,1}\right)& D{S}^h\left(\mathrm{1,2}\right)& ...& D{S}^h(1,n)\\ D{S}^h\left(\mathrm{2,1}\right)& ...& ...& ...\\ ...& ...& D{S}^h(i,j)& ...\\ D{S}^h(m,1)& D{S}^h(m,2)& ...& D{S}^h(m,n)\end{array}\right]$

$D{S}^h(i,j)=\left[\underset{u=1}{\overset{m}{\mathrm{\Σ}}}\underset{v=1}{\overset{n}{\mathrm{\Σ}}}\right|D{S}^1(i,j)-D{S}^1(u,v)\left|\right]/(m\×n)$

式中DS¹(i，j)表示第一次迭代后图像DS¹中坐标为(i，j)的元素值，DS¹(u，v)表示第一次迭代后图像DS¹中坐标为(u，v)的元素值；

(4d)当达到最高迭代次数即h＝3时，按照下式计算，得到的DS³就是灰度差异图D1经过显著性度量得到的显著图D3＝DS³，DS³＝DS^h：

$D{S}^h=\left[\begin{array}{cccc}D{S}^h\left(\mathrm{1,1}\right)& D{S}^h\left(\mathrm{1,2}\right)& ...& D{S}^h(1,n)\\ D{S}^h\left(\mathrm{2,1}\right)& ...& ...& ...\\ ...& ...& D{S}^h(i,j)& ...\\ D{S}^h(m,1)& D{S}^h(m,2)& ...& D{S}^h(m,n)\end{array}\right]$

$D{S}^h(i,j)=\left[\underset{u=1}{\overset{m}{\mathrm{\Σ}}}\underset{v=1}{\overset{n}{\mathrm{\Σ}}}\right|D{S}^2(i,j)-D{S}^2(u,v)\left|\right]/(m\×n)$

式中DS²(i，j)表示第二次迭代后图像DS²中坐标为(i，j)的元素值，DS²(u，v)表示第二次迭代后图像DS²中坐标为(u，v)的元素值。

步骤5，分别对差异图D2和显著图D3采用模糊c均值聚类算法进行分割，每幅图像均被分为变化类和非变化类这两类，将变化类的像素点赋值为1，非变化类的像素点赋值为0，得到分割结果图Z1和初始变化检测结果图Z2，其中Z1为差异图D2的分割结果，Z2为显著图D3的分割结果。

步骤6，利用分割结果图Z1对初始变化检测结果图Z2依据不同情况进行修正，去除初始变化检测图Z2中的孤立小噪声点，生成变化检测结果图Z，具体按照以下步骤进行：

对初始变化检测结果图Z2中的每个像素点，分别按以下三种情况进行处理：

(1)当初始变化检测结果图Z2中的元素Z2(i，j)＝1时，如果分割结果图Z1中以像素点(i，j)为中心，窗口大小为7×7邻域内存在被标记为变化类的像素点，即

则令变化检测结果图Z中的元素Z(i，j)＝Z2(i，j)，i，j表示当前像素点的坐标索引，即第i行，第j列；

(2)当初始变化检测结果图Z2中的元素Z2(i，j)＝1时，如果分割结果图Z1中以像素点(i，j)为中心，窗口大小为7×7邻域内不存在被标记为变化类的像素点，即

则令变化检测结果图Z中的元素Z(i，j)＝0；

(3)当初始变化检测结果图Z2中的元素Z2(i，j)＝0时，令变化检测结果图Z中的元素Z(i，j)＝Z2(i，j)。

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

1.实验数据

本发明的实验数据为一组模拟遥感数据集和两组真实遥感数据集，如图2所示。

模拟数据集的原始图像为470×335像素大小的位于英国Feltwell村庄农区的图像，配准误差为1.5个像素左右，模拟变化图像是通过模拟地球的天气变化和电磁波的辐射特性等因素影响并人工嵌入一些变化区域得到的，分别如图2(a)和图2(b)所示。

第一组真实遥感图像数据集原始图像及参考变化图如图2(d)和(e)所示。该组真实遥感数据集由1995年9月和1996年7月在意大利撒丁岛Mulargia湖泊区域的Landsat5TM第5波段光谱图像组成。两幅图像大小均为300×412像素，灰度级为256，它们之间的变化为湖水水位上涨引起的，发生变化区域的参考图如图2(f)所示，其中包括115974个非变化像素和7626个变化像素，图中白色像素区表示发生变化的区域。

第二组真实遥感图像数据集原始图像及参考变化图如图2(g)和(h)所示，该组真实遥感数据集是分别于1994年8月和1994年9月在意大利Elba岛西部地区的两时相Landsat-5TM第4波段多光谱图像组成，两幅图像大小均为326×414，灰度级为256，它们之间发生的变化是由于森林火灾破坏了大量植被所致，如图2(g)中左上角较暗的区域，其参考变化图如图2(i)所示，图中包含2415个变化像素和132549个非变化像素，白色像素区域表示变化的区域。

2.实验方法

方法1：徐丽华等在文章“顾及配准误差的遥感影像变化检测”中提出的一种基于相对辐射校正和区域变化率的变化检测方法；

方法2：Minghui Tian等学者在文章“A novel approach for change detection inremote sensing image based on saliency map.”提出的一种视觉注意模型的图像变化检测方法；

方法3：本发明方法。

3.评价指标

本发明采用常见的变化检测评价指标对实验结果进行定量评价，包括：虚警数、漏检数和总错误数。虚景数指实际没发生变化却被检测为变化的像素点个数之和，漏检数指实际发生变化却未被检测出的像素点个数之和，总错误数为虚景数与漏检数之和。

4.实验内容与结果分析

实验1，用不同方法对如图2(a)和图2(b)所示的两幅不同时相的模拟数据集进行变化检测，变化检测参考图如图2(c)所示，对图2(a)和图2(b)的模拟遥感数据集进行变化检测得到的结果如图3所示，其中图3(a)为现有方法1得到的变化检测结果图，图3(b)为现有方法2得到的变化检测结果图，图3(c)为本发明方法得到的变化检测结果图，从图3可以看出，现有方法1和现有方法2结果图中存在少量较明显的伪变化信息，本发明方法变化区域边缘保持较好，伪变化信息很少，比较接近参考图2(c)。

实验2，用不同方法对如图2(d)和图2(e)所示的两幅不同时相的真实遥感数据集进行变化检测，变化检测参考图如图2(f)所示，对图2(d)和图2(e)的真实遥感数据集进行变化检测得到的结果如图4所示，图4(a)为现有方法1得到的变化检测结果图，图4(b)为现有方法2得到的变化检测结果图，图4(c)为本发明方法得到的变化检测结果图，从图4可以看出，现有方法1结果图中都存在许多明显的噪声区域，现有方法2中的伪变化信息有所减少，但是变化区域存在明显的漏检，区域一致性保持不是很好，本发明方法在抑制虚景的同时能有效减少漏检数，结果图和参考图2(f)最为接近，变化检测结果的正确率最高。

实验3，用不同方法对如图2(g)和图2(h)所示的两幅不同时相的真实遥感数据集进行变化检测，变化检测参考图如图2(i)所示，对图2(g)和图2(h)的真实遥感数据集进行变化检测得到的结果如图5所示，图5(a)为现有方法1得到的变化检测结果图，图5(b)为现有方法2得到的变化检测结果图，图5(c)为本发明方法得到的变化检测结果图，从图5可以看出，现有方法1结果图中虽然变化区域保持较完整，但是噪声过多，现有方法2中变化区域的边缘信息丢失较多，本发明方法能较好的保持变化区域的几何信息，抑制伪变化信息的生成，结果图和参考图2(i)最为接近，变化检测结果的正确率最高。

对实验变化检测结果进行定量分析，即比较变化检测结果图与标准参考图，计算总错误数、虚景数、漏检数等性能指标，结果如表1所示。

表1.三组遥感图像数据变化检测结果评价指标

从表1上可以更为直观的看出，本发明方法能有效稳定地的获得遥感图像的变化区域，与其它两种现有方法相比较变化检测结果的精度最高。

Claims (3)

1.一种基于显著性度量的遥感图像变化检测方法，包括如下步骤：

(1)分别对输入的两幅已配准的大小均为m×n的两时相遥感图像I¹和I²进行中值滤波，得到两时相遥感图像I¹和I²的滤波后图像X¹和X²；

(2)对滤波后图像X¹和X²作差得到灰度差异图D1；

(3)计算滤波后图像X¹和X²基于Chi-square距离的差异图D2：

$D2=\left[\begin{array}{cccc}D2\left(\mathrm{1,1}\right)& D2\left(\mathrm{1,2}\right)& ...& D2(1,n)\\ D2\left(\mathrm{2,1}\right)& ...& ...& ...\\ ...& ...& D2(i,j)& ...\\ D2(m,1)& D2(m,2)& ...& D2(m,n)\end{array}\right]$

$D2(i,j)=\underset{k=1}{\overset{n1}{\mathrm{\Σ}}}{\left(\frac{X^1(i,j)-X^2(i,j)}{{\mathrm{\σ}}_{D1}}\right)}^2$

式中i，j分别表示当前像素点的坐标索引，σ表示标准差，∑表示求和，n1表示以当前像素点为中心，所取滑动窗口内像素点的总个数，本发明中取n1＝9，即滑动窗口大小为3×3；

(4)对灰度差异图D1进行显著性度量，得到显著图D3；

(5)对差异图D2和显著图D3分别进行模糊c均值聚类分割，得到两组分割结果，依次表示为Z1和Z2，其中分割结果图Z1和初始变化检测结果图Z2中被标记为变化类的像素点赋值为1，被标记为非变化类的像素点赋值为0；

(6)利用分割结果图Z1对初始变化检测结果图Z2依据不同情况进行修正，去除分割结果图Z2中的孤立小噪声点，生成变化检测结果图Z。

2.根据权利要求1所述的遥感图像变化检测方法，其中步骤(4)所述的对灰度差异图D1进行显著性度量，按以下步骤执行：

(4a)将灰度差异图D1作为初始图像DS⁰，即DS⁰＝D1；

(4b)当迭代次数h＝1时，对初始图像DS⁰中的每个像素点按照下式计算，得到第一次迭代后的图像DS¹＝DS^h：

$D{S}^h=\left[\begin{array}{cccc}D{S}^h\left(\mathrm{1,1}\right)& D{S}^h\left(\mathrm{1,2}\right)& ...& D{S}^h(1,n)\\ D{S}^h\left(\mathrm{2,1}\right)& ...& ...& ...\\ ...& ...& D{S}^h(i,j)& ...\\ D{S}^h(m,1)& D{S}^h(m,2)& ...& D{S}^h(m,n)\end{array}\right]$

$D{S}^h(i,j)=\left[\underset{u=1}{\overset{m}{\mathrm{\Σ}}}\underset{v=1}{\overset{n}{\mathrm{\Σ}}}\right|D{S}^0(i,j)-D{S}^0(u,v)\left|\right]/(m\×n)$

式中i和j分别表示当前像素点的坐标索引，u和v分别为图像中除当前像素点外其它像素点的坐标索引，m，n表示图像大小，DS⁰(i，j)表示初始图像DS⁰中坐标为(i，j)的元素值，DS⁰(u，v)表示初始图像DS⁰中坐标为(u，v)的元素值；

(4c)当迭代次数h＝2时，对第一次迭代后的图像DS¹中的每个像素点按照下式计算，得到第二次迭代后的图像DS²＝DS^h：

$D{S}^h=\left[\begin{array}{cccc}D{S}^h\left(\mathrm{1,1}\right)& D{S}^h\left(\mathrm{1,2}\right)& ...& D{S}^h(1,n)\\ D{S}^h\left(\mathrm{2,1}\right)& ...& ...& ...\\ ...& ...& D{S}^h(i,j)& ...\\ D{S}^h(m,1)& D{S}^h(m,2)& ...& D{S}^h(m,n)\end{array}\right]$

$D{S}^h(i,j)=\left[\underset{u=1}{\overset{m}{\mathrm{\Σ}}}\underset{v=1}{\overset{n}{\mathrm{\Σ}}}\right|D{S}^1(i,j)-D{S}^1(u,v)\left|\right]/(m\×n)$

式中DS¹(i，j)表示第一次迭代后图像DS¹中坐标为(i，j)的元素值，DS¹(u，v)表示第一次迭代后图像DS¹中坐标为(u，v)的元素值；

(4d)当达到最高迭代次数即h＝3时，按照下式计算，得到的DS³就是灰度差异图D1经过显著性度量得到的显著图D3＝DS³，DS³＝DS^h：

$D{S}^h=\left[\begin{array}{cccc}D{S}^h\left(\mathrm{1,1}\right)& D{S}^h\left(\mathrm{1,2}\right)& ...& D{S}^h(1,n)\\ D{S}^h\left(\mathrm{2,1}\right)& ...& ...& ...\\ ...& ...& D{S}^h(i,j)& ...\\ D{S}^h(m,1)& D{S}^h(m,2)& ...& D{S}^h(m,n)\end{array}\right]$

$D{S}^h(i,j)=\left[\underset{u=1}{\overset{m}{\mathrm{\Σ}}}\underset{v=1}{\overset{n}{\mathrm{\Σ}}}\right|D{S}^2(i,j)-D{S}^2(u,v)\left|\right]/(m\×n)$

式中DS²(i，j)表示第二次迭代后图像DS²中坐标为(i，j)的元素值，DS²(u，v)表示第二次迭代后图像DS²中坐标为(u，v)的元素值。

3.根据权利要求1所述的遥感图像变化检测方法，步骤(6)所述的利用分割结果图Z1对初始变化检测结果Z2依据不同情况进行修正，生成变化检测结果图Z，按照以下步骤进行：

对初始变化检测结果图Z2中的每个像素点，分别按以下三种情况进行处理：

(1)当初始变化检测结果图Z2中的元素Z2(i，j)＝1时，如果分割结果图Z1中以像素点(i，j)为中心，窗口大小为7×7邻域内存在被标记为变化类的像素点，即

则令变化检测结果图Z中的元素Z(i，j)＝Z2(i，j)，i，j表示当前像素点的坐标索引，即第i行，第j列；

(2)当初始变化检测结果图Z2中的元素Z2(i，j)＝1时，如果分割结果图Z1中以像素点(i，j)为中心，窗口大小为7×7邻域内不存在被标记为变化类的像素点，即

则令变化检测结果图Z中的元素Z(i，j)＝0；

(3)当初始变化检测结果图Z2中的元素Z2(i，j)＝0时，令变化检测结果图Z中的元素Z(i，j)＝Z2(i，j)。

Images