CN102968790A - 基于图像融合的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于图像融合的SAR图像变化检测方法，主要解决现有技术中细节信息丢失和检测精度低的问题。其实现过程是：提取变化前后图像的差值差异图和比值差异图，并提取差值图像的方差矩阵；用KI阈值分割法，得到最佳阈值T；根据最佳阈值T和方差矩阵将差值图与比值图融合，得到融合后的最终差异图D；对所得的最终差异图D的灰度值进行类别数为2的模糊C均值聚类分割，将聚类中心取值较大的一类作为变化类，另一类为非变化类，得到变化监测结果。本发明克服了差值图对噪声敏感和比值图边缘信息丢失的缺点，降低了错误率，更好的保存了细节信息。

Description

基于图像融合的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及遥感图像变化检测，即根据不同时相的SAR图像检测地表产生的变化，可用于地物覆盖与利用、自然灾害监控，城区规划，地图更新等领域。

背景技术

遥感图像的变化检测是指通过分析在不同时间来自同一地区的两幅或多幅遥感图像，检测出该地区的地物随时间发生的变化信息。遥感图像的变化检测已经广泛地应用于如森林资源的动态监测、土地覆盖和利用的变化监测、城市规划布局、环境监测分析、自然灾害评估、以及在军事中的人造目标监测和地面武装部署分析等许多领域。

在变化检测研究中，学者们将现有的多种变化检测方法归为不同的类别，其中最常见的是基于差异图像分析的变化检测方法。该方法通常包括3个关键的步骤：1)图像预处理；2)构造差异图像；3)分析差异图像提取变化检测结果。图像预处理包括图像去噪、配准、辐射校正等。构造差异图像是基于差异图像分析方法的重要一步，尤其是构造视觉效果较好的差异图像难度比较大。其中，图像差值法和比值法是最常见的差异图提取算子，图像差值法对SAR图像的变化区域非常敏感，有利于保留变化信息的几何边缘，从而有利于精确地提取变化区域，但由于SAR图像相干斑噪声的影响，差值法并不适用于分析受噪声影响较大的SAR图像。而图像对数比虽然能有效地抑制SAR图像的相干斑噪声，但却对变化区域的边缘造成一定的模糊，特别是存在大量的微弱变化和小区域变化信息的图像中，不利于变化区域轮廓和细节信息的保持。

发明内容

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

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

（1）输入变化前后的SAR图像P₁和P₂，分别计算出两幅图对应像素灰度值的差并归一化，得到差值图像S的灰度矩阵H_S，计算两幅图像P₁和P₂对应素灰度值的商对数值并归一化，得到比值图R的灰度矩阵H_R；

（2）计算差值图S每一个像素灰度值S^ab的局部方差值V^ab，得到整幅差值图的方差矩阵V，并对该方差矩阵V进行KI阈值分割，得到最优阈值T*,其中a,b为像素点在图像中的空域坐标变量；

（3）根据方差矩阵V和最佳阈值T^*，对差值差异图S和比值差异图R进行融合，得到融合后的最终差异图D在每一个像素点的灰度值：

$D^{\mathrm{ab}}=\left\{\begin{array}{cc}{S}^{\mathrm{ab}}& V^{\mathrm{ab}}\&GreaterEqual;T^{*}\\ R^{\mathrm{ab}}& V^{\mathrm{ab}}<T^{*}\end{array}\right.,$

其中a,b为像素点的空域坐标变量，S^ab，R^ab，D^ab分别为差值图S，比值图R和融合后的最终差异图D在像素点(a,b)的灰度值，V^ab为差值差异图S^ab的在像素点(a,b)的局部方差值；

（4）对融合后的最终差异图D的灰度值进行类别数为2的模糊C均值聚类分割，聚类结束后，聚类中心值较大的一类为变化类，另一类为非变化类。

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

1.本发明主要以差值图和比值图为基础，操作简单，运算速度快；

2.本发明由于将方差信息对差值图与比值图进行融合，因而克服了差值图对噪声敏感和比值图易丢失边缘信息的缺点，降低了变化检测的错误率，更好地保持了细节信息。

附图说明

图1是本发明的流程图；

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

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

图4是用本发明和现有差值法，比值法对图2的变化检测结果图；

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

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

图7是用本发明和现有差值法，比值法对图5的变化检测结果图。

具体实施方式

参照图1，本发明图像融合的遥感图像变化检测方法，包括如下步骤：

步骤1:输入变化前后的SAR图像P₁和P₂，计算出两幅图对应像素灰度值的差，并归一化得到差值图S的灰度矩阵H_R：

1a)计算输入图像P₁在像素点(a,b)的灰度值与输入图像P₂在像素点(a,b)的灰度值的差Z_S ^ab：

${Z_S}^{\mathrm{ab}}=|P_1^{\mathrm{ab}}-P_2^{\mathrm{ab}}|,$

其中，

和

分别为输入图像P₁在像素点(a,b)的灰度值和输入图像P₂在像素点(a,b)的灰度值，a,b为像素点在图像中的空域坐标；

1b)归一化Z_S ^ab，得到差值图S在每一个像素点(a,b)的灰度值S^ab，

$S^{\mathrm{ab}}=255\&times;\frac{Z_S^{\mathrm{ab}}-\min \left(Z_s\right)}{\max \left(Z_S\right)-\min \left(Z_S\right)},$

其中，Z_s为输入图像P₁与P₂对应像素灰度相减产生的矩阵；

1c)根据每一个像素点(a,b)的灰度值S^ab，得到差值图S的灰度矩阵H_S

H_s={S^ab}。

步骤2：计算两幅图对应像素灰度值商的对数值并归一化，得到比值图R的灰度矩阵H_R：

2a)计算输入图像P₁在像素点(a,b)的灰度值与输入图像P₂在像素点(a,b)的灰度值的商对数Z_R ^ab：

${Z_R}^{\mathrm{ab}}=\left\{\begin{array}{cc}\ln \left(\frac{P_1^{\mathrm{ab}}}{P_2^{\mathrm{ab}}}\right)& \text{if}(P_1^{\mathrm{ab}}\&GreaterEqual;P_2^{\mathrm{ab}})\\ \ln \left(\frac{P_2^{\mathrm{ab}}}{P_1^{\mathrm{ab}}}\right)& \mathrm{else}\end{array}\right.,$

其中，

和

分别为输入图像P₁在像素点(a,b)的灰度值和输入图像P₂在像素点(a,b)的灰度值；

2b)归一化Z_R ^ab，得到比值图R在每一个像素点(a,b)的灰度值R^ab，

$R^{\mathrm{ab}}=255\&times;\frac{Z_R^{\mathrm{ab}}-\min \left(Z_R\right)}{\max \left(Z_R\right)-\min \left(Z_R\right)},$

其中，Z_R为输入图像P₁与P₂对应像素灰度取商的对数得到的矩阵；

2c)根据每一个像素点(a,b)的灰度值R^ab，得到差值图S的灰度矩阵H_S

H_R={R^ab}。

步骤3:计算差值图S每一个像素灰度值S^ab的局部方差值V^ab：

$V^{\mathrm{ab}}=\frac{1}{3\&times;3}\underset{m=a-1}{\overset{a+1}{\mathrm{\&Sigma;}}}\underset{n=b-1}{\overset{b+1}{\mathrm{\&Sigma;}}}{(S^{\mathrm{mn}}-\frac{1}{3\&times;3}\underset{m=a-1}{\overset{a+1}{\mathrm{\&Sigma;}}}\underset{n=b-1}{\overset{b+1}{\mathrm{\&Sigma;}}}{S}^{\mathrm{mn}})}^2,$

其中，a,b为像素点的在图像中的空域坐标变量，m,n为局部窗口内像素点的坐标变量，m=a-1,a,a+1，n=b-1,b,b+1，窗口大小为3x3，S^mn为差值图S在局部窗口内像素点(m,n)的灰度值。

步骤4：对差值图S每一个像素灰度值S^ab的局部方差值V^ab进行KI阈值分割，得到最优阈值T^*：

1a）在给定灰度等级x和阈值T情况下，分别计算变化类w_c的后验概率p(w_c|x,T)和非变化类w_u的后验概率p(w_u|x,T)：

$p\left(w_c\right|x,T)=\frac{p\left(w_c\right)p\left(x\right|w_c,T)}{p\left(w_u\right)p\left(x\right|w_u,T)+p\left(w_c\right)p\left(x\right|w_c,T)},$

$p\left(w_u\right|x,T)=\frac{p\left(w_u\right)p\left(x\right|w_u,T)}{p\left(w_u\right)p\left(x\right|w_u,T)+p\left(w_c\right)p\left(x\right|w_c,T)},$

其中，x=1,2,…255，T＝1,2,…255，p(x|w_c,T)为变化类的对数正态似然概率，p(x|w_u,T)为非变化类的对数正态似然概率，p(w_c)为变化类的先验概率，p(w_u)为非变化类的先验概率；

1b）计算分类误差J(T)

$J\left(T\right)=\underset{x=0}{\overset{255}{\mathrm{\&Sigma;}}}h\left(x\right)c(x,T)$

这里， $c(x,T)=\left\{\begin{array}{cc}-2\ln P\left(w_u\right|x,T)& \mathrm{ifx}\&le;T\\ -2\ln P\left(w_c\right|x,T)& \mathrm{ifx}>T\end{array}\right.,$ h(x)为灰度值为x的像素点在图像中出现的概率；

2c）将分类误差J(T)达到最小值时的T，作为最佳分类阈值T^*：

T^*＝argminJ(T)，T＝0,1,....255。

步骤5，根据方差矩阵V和最佳阈值T^*，对差值差异图S和比值差异图R进行融合，得到融合后的最终差异图D在每一个像素点的灰度值：

$D^{\mathrm{ab}}=\left\{\begin{array}{cc}{S}^{\mathrm{ab}}& V^{\mathrm{ab}}\&GreaterEqual;T^{*}\\ R^{\mathrm{ab}}& V^{\mathrm{ab}}<T^{*}\end{array}\right.,$

其中，a,b为像素点的空域坐标变量，S^ab，R^ab，D^ab分别为差值图S，比值图R和融合后的最终差异图D在像素点(a,b)的灰度值，V^ab为差值差异图S^ab的在像素点(a,b)的局部方差值。

步骤6:对所得的最终差异图D的灰度值进行类别数为2的模糊C均值聚类分割，得到聚类中心c₁，c₂：

6a)将待聚类图像的灰度矩阵转化为一个列向量X=[x_i]，x_i为该向量第i个像素的灰度值，i=1,2,3,…N，N为总像素点数，将第i个像素属于第k类的隶属度记为u_ki，u_ki∈[0,1]且

k表示两类的类标，k=1,2，随机初始化隶属度矩阵U^b=[u_ki]的每一个元素，设定循环计数器b=0；

6b)按如下公式计算k类的聚类中心c_k：

$c_k=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ki}}^2{x}_i}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ki}}^2};$

6c)按如下公式计算隶属度矩阵U^(b+1）:

U^（b+1）={u_ki}，

其中， $u_{\mathrm{ki}}=\frac{1}{{\left(\frac{x_i-c_k}{x_i-c_1}\right)}^2+{\left(\frac{x_i-c_k}{x_i-c_2}\right)}^2};$

6d)设定聚类的停止条件ε为0.001，如果max{U^(b)-U^(b+1)}<ε则停止迭代，否则设置b=b+1,转至步骤4b)继续迭代，直到满足条件为止，聚类结束，聚类中心c₁和c₂；

6e)通过比较c₁与c₂的大小确定变化类与非变化类，如果c₁≥c₂，则为以c₁为聚类中心的那一类为变化类，以c₂为聚类中心的那一类为非变化类；如果c₁<c₂，为以c₂为聚类中心的那一类为变化类，以c₁为聚类中心的一类为非变化类。

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

1.实验条件

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

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

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

2.实验内容和实验结果

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

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

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

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

试验结果数据表

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

从上表可以看出，由于本发明将方差信息对差值图与比值图进行融合，与所述的两种传统变化检测方法相比，可以获得较少的漏检数和误检数，以及最少的总错误数，提高了变化检测的检测精度。

Claims (5)

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

（1）输入变化前后的SAR图像P₁和P₂，分别计算出两幅图对应像素灰度值的差并归一化，得到差值图像S的灰度矩阵H_S，计算两幅图像P₁和P₂对应素灰度值的商对数值并归一化，得到比值图R的灰度矩阵H_R；

（2）计算差值图S每一个像素灰度值S^ab的局部方差值V^ab，得到整幅差值图的方差矩阵V，并对该方差矩阵V进行KI阈值分割，得到最优阈值T^*,其中a,b为像素点在图像中的空域坐标变量；

（3）根据方差矩阵V和最佳阈值T^*，对差值差异图S和比值差异图R进行融合，得到融合后的最终差异图D在每一个像素点的灰度值：

$D^{\mathrm{ab}}=\left\{\begin{array}{cc}{S}^{\mathrm{ab}}& V^{\mathrm{ab}}\&GreaterEqual;T^{*}\\ R^{\mathrm{ab}}& V^{\mathrm{ab}}<T^{*}\end{array}\right.,$

其中a,b为像素点的空域坐标变量，S^ab，R^ab，D^ab分别为差值图S，比值图R和融合后的最终差异图D在像素点(a,b)的灰度值，V^ab为差值差异图S^ab的在像素点(a,b)的局部方差值；

（4）对融合后的最终差异图D的灰度值进行类别数为2的模糊C均值聚类分割，聚类结束后，聚类中心值较大的一类为变化类，另一类为非变化类。

2.根据权利要求书1中所述的方法，其中步骤(1)所述的计算出两幅图对应像素灰度值的差并归一化，得到差值图像S的灰度矩阵H_S，按如下步骤进行：

1a)计算输入图像P₁在像素点(a,b)的灰度值与输入图像P₂在像素点(a,b)的灰度值的差Z_S ^ab：

${Z_S}^{\mathrm{ab}}=|P_1^{\mathrm{ab}}-P_2^{\mathrm{ab}}|,$

其中，

和

分别为输入图像P₁在像素点(a,b)的灰度值和输入图像P₂在像素点(a,b)的灰度值，a,b为像素点在图像中的空域坐标；

1b)归一化Z_S ^ab，得到差值图S在每一个像素点(a,b)的灰度值S^ab:

$S^{\mathrm{ab}}=255\&times;\frac{Z_S^{\mathrm{ab}}-\min \left(Z_s\right)}{\max \left(Z_S\right)-\min \left(Z_S\right)},$

其中，Z_s为输入图像P₁与P₂对应像素灰度相减产生的矩阵；

1c)根据每一个像素点(a,b)的灰度值S^ab，得到差值图S的灰度矩阵H_S

H_s={S^ab}。

3.根据权利要求书1中所述的方法，其中步骤(1)中计算两幅图像P₁和P₂对应像素灰度值的商对数值并归一化，得到差值图像R的灰度矩阵H_R，按如下步骤进行：

1.1)计算输入图像P₁在像素点(a,b)的灰度值与输入图像P₂在像素点(a,b)的灰度值的商对数Z_R ^ab：

${Z_R}^{\mathrm{ab}}=\left\{\begin{array}{cc}\ln \left(\frac{P_1^{\mathrm{ab}}}{P_2^{\mathrm{ab}}}\right)& \mathrm{if}(P_1^{\mathrm{ab}}\&GreaterEqual;P_2^{\mathrm{ab}})\\ \ln \left(\frac{P_2^{\mathrm{ab}}}{P_1^{\mathrm{ab}}}\right)& \mathrm{else}\end{array}\right.,$

其中，

和

分别为输入图像P₁在像素点(a,b)的灰度值和输入图像P₂在像素点(a,b)的灰度值；

1.2)归一化Z_R ^ab，得到比值图R在每一个像素点(a,b)的灰度值R^ab，

$R^{\mathrm{ab}}=255\&times;\frac{Z_R^{\mathrm{ab}}-\min \left(Z_R\right)}{\max \left(Z_R\right)-\min \left(Z_R\right)},$

其中，Z_R为输入图像P₁与P₂对应像素灰度取商的对数得到的矩阵；

1.3)根据每一个像素点(a,b)的灰度值R^ab，得到差值图S的灰度矩阵H_S

H_R={R^ab}。

4.根据权利要求书1中所述的方法，其中步骤(2)中所述的计算差值图S每一个像素灰度值S^ab的局部方差值V^ab，按如下公式计算：

$V^{\mathrm{ab}}=\frac{1}{3\&times;3}\underset{m=a-1}{\overset{a+1}{\mathrm{\&Sigma;}}}\underset{n=b-1}{\overset{b+1}{\mathrm{\&Sigma;}}}{(S^{\mathrm{mn}}-\frac{1}{3\&times;3}\underset{m=a-1}{\overset{a+1}{\mathrm{\&Sigma;}}}\underset{n=b-1}{\overset{b+1}{\mathrm{\&Sigma;}}}{S}^{\mathrm{mn}})}^2,$

其中，a,b为像素点的在图像中的空域坐标变量，m,n为局部窗口内像素点的坐标变量，m=a-1，a,a+1，n=b-1,b,b+1，窗口大小为3x3，S^mn为差值图S在局部窗口内像素点(m,n)的灰度值。

5.根据权利要求书1所述的方法，其中步骤(4)中对融合后的最终差异图D的灰度值进行类别数为2的模糊C均值聚类分割，按照如下步骤实现：

4a)将待聚类图像的灰度矩阵转化为一个列向量X=[x_i]，x_i为该向量第i个像素的灰度值，i＝1,2,3,…N，N为总像素点数，将第i个像素属于第k类的隶属度记为u_ki，u_ki∈[0,1]且

k表示两类的类标，k=1,2，随机初始化隶属度矩阵U^b=[u_ki]的每一个元素，设定循环计数器b=0；

4b)按如下公式计算k类的聚类中心c_k：

$c_k=\frac{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ki}}^2{x}_i}{\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ki}}^2};$

4c)按如下公式计算隶属度矩阵U^(b+1）:

U^（b+1）={u_ki}，

其中， $u_{\mathrm{ki}}=\frac{1}{{\left(\frac{x_i-c_k}{x_i-c_1}\right)}^2+{\left(\frac{x_i-c_k}{x_i-c_2}\right)}^2};$

4d)如果max{U^(b)-U^(b+1)}<0.001则停止迭代，否则设置b=b+1,转至步骤4b)继续迭代，直到满足条件为止；

4e)停止迭代后，通过比较c₁与c₂的大小确定变化类与非变化类，如果c₁≥c₂，则为以c₁为聚类中心的那一类为变化类，以c₂为聚类中心的那一类为非变化类；如果c₁<c₂，为以c₂为聚类中心的那一类为变化类，以c₁为聚类中心的一类为非变化类。

Images