CN105389817A - 一种两时相遥感影像变化检测方法 - Google Patents

Abstract

本发明公开了一种两时相遥感影像变化检测方法，包括：获取两时相的原始光学遥感影像；对所述原始光学遥感影像进行影像配准；利用直方图调整方法对影像配准后的遥感影像进行辐射归一化校正；对辐射归一化校正后的遥感影像，计算每个相应波段的波段差；对每一个波段差，分别计算扩展属性剖面；将每一个波段差作为FCM的输入，并将扩展属性剖面作为空间邻域信息嵌入到FCM模型中，进而确定影像的变化区域和非变化区域。本发明可以提高变化检测的精度。

Description

一种两时相遥感影像变化检测方法

技术领域

本发明涉及遥感影像处理技术领域，尤其涉及一种两时相遥感影像变化检测方法。

背景技术

遥感影像的变化检测就是从不同时期的遥感数据中，定量地分析和确定地表变化的特征与过程。各国学者从不同的角度和应用研究提出了许多有效的检测算法，如变化矢量分析法(ChangeVectorAnalysis，CVA)、基于FuzzyC-means(FCM)的聚类方法等。其中，传统的基于FCM的多时相光学遥感变化检测，多先进行CVA变换，然后对变化矢量的幅值进行FCM聚类，进而得到变化检测结果。该类技术中，使用FCM的不足是仅仅使用变化矢量的幅值，使得原始多波段信息没有得到充分的挖掘。

基于FCM的变化检测算法中，由于存在空间相关性，图像中的像素点所属类别(变化类与非变化类)不仅与当前像素点的灰度值有关，还受到其邻域内像素点灰度值的影响。目前在FCM算法的基础上，目前多以局部空间邻域均值或中值作为空间约束条件嵌入到FCM的目标函数中。然而这些方法使用固定的窗口得到空间信息，导致了尺寸比例的选择问题。另外，由于多时相光学遥感影像检测环境的复杂化以及目标先验信息匮乏等，导致这些算法都存在着一定的局限性。为此，有必要研究新的光学遥感图像变化检测技术来有效克服上述难点。

发明内容

本发明所要解决的技术问题在于，提供一种两时相遥感影像变化检测方法，基于扩展属性滤波(ExtendedAttributeProfile，EAP)和FCM(FuzzyC-Means)，提高了变换检测精度。

为了解决上述技术问题，本发明提供了一种两时相遥感影像变化检测方法，包括：

获取两时相的原始光学遥感影像；

对所述原始光学遥感影像进行影像配准；

利用直方图调整方法对影像配准后的遥感影像进行辐射归一化校正；

对辐射归一化校正后的遥感影像，计算每个相应波段的波段差；

对每一个波段差，分别计算扩展属性剖面；

将每一个波段差作为FCM的输入，并将扩展属性剖面作为空间邻域信息嵌入到FCM模型中，确定影像的变化区域和非变化区域。

实施本发明，具有如下有益效果：

1)采用EAP提取影像的纹理特征，其滤波窗口尺寸的选择不是固定的，而是自适应变化的。

2)在FCM目标函数中，利用EAP滤波影像作为空间约束，提高了影像的变化检测精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明提供的两时相遥感影像变化检测方法的一个实施例的流程示意图；

图2是原始光学遥感影像；

图3是本发明方法与其他方法的处理结果对比图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1是本发明提供的两时相遥感影像变化检测方法的一个实施例的流程示意图，如图1所示，包括：

S101、获取两时相的原始光学遥感影像X₁和X₂。

其中，X₁、X₂是同一区域不同时相的两幅光学遥感影像。

S102、对所述原始光学遥感影像X₁和X₂进行影像配准。

具体的，步骤S102具体包括步骤：

S1021、采用ENVI14.8遥感软件对光学遥感影像X₁和X₂进行几何粗校正。

几何粗校正具体操作步骤为：(1)显示基准影像和待校正影像；(2)采集地面控制点GCPs；GCPs应均匀分布在整幅图像内，GCPs的数目至少大于等于9；(3)计算误差；(4)选择多项式模型；(5)采用双线性插值进行重采样输出。

其中的双线性差值法为：若求未知函数f1在点P＝(x,y)的值，假设我们已知函数f1在Q₁₁＝(x₁,y₁),Q₁₂＝(x₁,y₂),Q₂₁＝(x₂,y₁),及Q₂₂＝(x₂,y₂)四个点的值。如果选择一个坐标系统使得这四个点的坐标分别为(0,0)、(0,1)、(1,0)和(1,1)，那么双线性插值公式就可以表示为：

f1(x,y)≈f1(0,0)(1-x)(1-y)+f1(1,0)x(1-y)+f(0,1)(1-x)y+f1(1,1)xy。

S1022、利用自动匹配与三角剖分法对几何粗校正后的遥感影像X₁和X₂进行几何精校正。

其中，三角剖分法为，采用逐点插入法构建Delaunay三角网，对每一个三角形，利用其三个顶点的行列号与其对应的基准影像同名点的地理坐标来确定该三角形内部的仿射变换模型参数，对待校正影像进行纠正，得到校正后的遥感影像。

S103利用直方图调整方法对影像配准后的遥感影像进行辐射归一化校正。

具体的，步骤S103包括步骤：将影像配准后的遥感影像的像素灰度值分别调整为从而实现辐射归一化校正，

其中，影像配准后的遥感影像各个波段服从高斯分布，μ_f、σ_f分别为参考影像的像素均值和标准差，g_r、μ_r、σ_r分别为影像配准后的遥感影像的像素灰度值、均值和标准差。

本步骤中按照公式分别对影像配准后的遥感影像X₁和X₂进行辐射归一化校正。

S104、对辐射归一化校正后的遥感影像，计算每个相应波段的波段差。

其中，计算到的波段差为D_i，i＝1,2,…,B，B表示影像的波段数目。

S105、对每一个波段差，分别计算扩展属性剖面。

具体的，步骤S105包括步骤：

S1051、对每一个波段差D_i，设置其灰度级为L+1，其层级定义为X_l(D_i)＝{p∈D_i|D_i(p)≥l}，对影像D_i在每一灰度级上进行二值化，得到一系列二值图像X_l(D_i)_，p表示像素点，l表示影像的灰度值；

S1052、对于D_i上的每一个像素点p，根据不同的属性门限T_λ做属性开变换和属性闭变换，得到p点处的属性开剖面Π((γ^T(D_i))(p))和属性闭剖面其中，

$\mathrm{\Π}\left(\right({\mathrm{\γ}}^T\left(D_i\right)\left)\right(p\left)\right)=\{{\mathrm{\Π\γ}}^{T_{\mathrm{\λ}}}:{\mathrm{\Π\γ}}^{T_{\mathrm{\λ}}}=\left({\mathrm{\γ}}^{T_{\mathrm{\λ}}}\right(D_i\left)\right)\left(p\right)\},\∀\mathrm{\λ}\∈\[0,1,...,n\]$

其中，属性开变换和闭变换的计算步骤为：

对每一个二值图像I，其属性开变换定义为： 属性闭变换定义为Ψ^T(I)＝(Γ^T(I^c))^c。

因此，D_i的属性开变换为：(γ^T(D_i))(p)＝max{k:p∈Γ^T[X_l(D_i)]}

D_i的属性闭变换为：

S106、将每一个波段差作为FCM的输入，并将扩展属性剖面作为空间邻域信息嵌入到FCM模型中，确定影像的变化区域和非变化区域。

具体的，步骤S106包括步骤：

S1061、构建FCM结合空间邻域信息的模型为：

$J_m=\underset{j=1}{\overset{C}{\mathrm{\Σ}}}\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{u}_{jk}^m\left|\right|D\left(k\right)-v_j|{|}^2+\mathrm{\α}\[\underset{j=1}{\overset{C}{\mathrm{\Σ}}}\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{u}_{jk}^m||open\left(k\right)-v_j|{|}^2+\underset{j=1}{\overset{C}{\mathrm{\Σ}}}\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{u}_{jk}^m\left|\right|close\left(k\right)-v_j|{|}^2\]$

式中，C是聚类数目，C为2，N是样本的总数，表示第k样本对于第j类聚类中心的模糊隶属度，m为隶属度的加权指数，m为2，u_jk∈[0,1]且参数α控制惩罚效果，

D(k)＝[D₁(k)D₂(k)…D_B(k)]，

open(k)＝[Π((γ^T(D₁))(k))Π((γ^T(D₂))(k))…Π((γ^T(D_B))(k))]，

S1062、采用式 $u_{jk}=\frac{{\left(\right|\left|D\right(k)-v_j|{|}^2+\mathrm{\α}\left(\left|\right|open\left(k\right)-v_j|{|}^2+||close\left(k\right)-v_j|{|}^2\right))}^{-1/(m-1)}}{\underset{j=1}{\overset{C}{\Σ}}{\left(\right|\left|D\right(k)-v_j|{|}^2+\mathrm{\α}\left(\left|\right|open\left(k\right)-v_j|{|}^2+||close\left(k\right)-v_j|{|}^2\right))}^{-1/(m-1)}}$ 更新模糊划分矩阵；

S1063、采用式 $v_j=\frac{\underset{k=1}{\overset{N}{\Σ}}{u}_{jk}^m\left(D\right(k)+\mathrm{\α}(open\left(k\right)+close\left(k\right)\left)\right)}{(1+\mathrm{\α})\underset{k=1}{\overset{N}{\Σ}}{u}_{jk}^m}$ 更新聚类中心；

S1064、重复更新模糊划分矩阵和聚类中心，直至相邻两次迭代的聚类中心聚类小于阈值ε；

S1065、根据模糊划分矩阵u_jk确定最终的变化检测图，确定影像的变化区域和非变化区域。

本发明的实验数据为LandsatTM数据，位于巴西的亚马逊森林地区的2幅遥感影像，获取时间分别为2000年7月和2006年7月，选用前4个波段，实验区大小为320像元×320像元，图2(a)和(b)分别为两个时相的真彩色遥感影像。变化参考图如图2(c)所示，共有16,826个变化像素。

(1)基于CVA的EM方法(CVA-EM)[意大利的BruzzoneL.等在文章“Automaticanalysisofdifferenceimageforunsupervisedchangedetection”(IEEETransactionsonGeoscienceandRemoteSensing,2000,38(3):1171-1182.)中所提的检测方法]。

(2)FCM结合空间邻域信息的分类方法(FCM-S)[Chensongchan等在文章“RobustImageSegmentationUsingFCMWithSpatialConstraintsBasedonNewKernel-InducedDistanceMeasure”(IEEETransactionsonSystems,Man,andCybernetics—PartB:Cybernetics,2004,34(4):1907-1916.)中所提的方法]

(3)本发明方法(采用标准差属性，阈值设为30)。

检测性能用错检数FP、漏检数FN、总错误数OE和Kappa系数四个指标来衡量。FP、FN和OE越接近于0、Kappa系数越接近于1，表明变化检测方法的性能越好。检测结果如表1所示。由表1可见，本发明所提的检测方法性能优于其他两种检测方法，这表明本发明所提的变化检测方法是有效的。

表1各算法变化检测结果比较

实施本发明，具有如下有益效果：

1)采用EAP提取影像的纹理特征，其滤波窗口尺寸的选择不是固定的，而是自适应变化的。

2)在FCM目标函数中，利用EAP滤波影像作为空间约束，提高了影像的变化检测精度。

需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者装置不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者装置所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者装置中还存在另外的相同要素。

专业人员还可以进一步意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。软件模块可以置于随机存储器(RAM)、内存、只读存储器(ROM)、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM、或技术领域内所公知的任意其它形式的存储介质中。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (5)

1.一种两时相遥感影像变化检测方法，其特征在于，包括：

获取两时相的原始光学遥感影像；

对所述原始光学遥感影像进行影像配准；

利用直方图调整方法对影像配准后的遥感影像进行辐射归一化校正；

对辐射归一化校正后的遥感影像，计算每个相应波段的波段差；

对每一个波段差，分别计算扩展属性剖面；

将每一个波段差作为FCM的输入，并将扩展属性剖面作为空间邻域信息嵌入到FCM模型中，确定影像的变化区域和非变化区域。

2.如权利要求1所述的两时相遥感影像变化检测方法，其特征在于，所述对所述原始光学遥感影像进行影像配准，具体包括：

采用ENVI14.8遥感软件对原始光学遥感影像进行几何粗校正；

利用自动匹配与三角剖分法对几何粗校正后的遥感影像进行几何精校正。

3.如权利要求1所述的两时相遥感影像变化检测方法，其特征在于，所述利用直方图调整方法对影像配准后的遥感影像进行辐射归一化校正，具体包括：

将影像配准后的遥感影像的像素灰度值分别调整为从而实现辐射归一化校正，

其中，影像配准后的遥感影像各个波段服从高斯分布，μ_f、σ_f分别为参考影像的像素均值和标准差，g_r、μ_r、σ_r分别为影像配准后的遥感影像的像素灰度值、均值和标准差。

4.如权利要求1所述的两时相遥感影像变化检测方法，其特征在于，所述对每一个波段差，分别计算扩展属性剖面，具体包括：

对每一个波段差D_i，设置其灰度级为L+1，其层级定义为X_l(D_i)＝{p∈D_i|D_i(p)≥l}，对影像D_i在每一灰度级上进行二值化，得到一系列二值图像X_l(D_i)，其中，i＝1,2,…,B，B表示影像的波段数目，p表示像素点，l表示影像的灰度值；

对于D_i上的每一个像素点p，根据不同的属性门限T_λ做属性开变换和属性闭变换，得到p点处的属性开剖面Π((γ^T(D_i))(p))和属性闭剖面其中，

$\mathrm{\Π}\left(\right({\mathrm{\γ}}^T\left(D_i\right)\left)\right(p\left)\right)=\{{\mathrm{\Π\γ}}^{T_{\mathrm{\λ}}}:{\mathrm{\Π\γ}}^{T_{\mathrm{\λ}}}=\left({\mathrm{\γ}}^{T_{\mathrm{\λ}}}\right(D_i\left)\right)\left(p\right)\},\∀\mathrm{\λ}\∈\[0,1,...,n\]$

5.如权利要求4所述的两时相遥感影像变化检测方法，其特征在于，所述将每一个波段差作为FCM的输入，并将扩展属性剖面作为空间邻域信息嵌入到FCM模型中，确定影像的变化区域和非变化区域，具体包括：

构建FCM结合空间邻域信息的模型为：

$J_m=\underset{j=1}{\overset{C}{\Σ}}\underset{k=1}{\overset{N}{\Σ}}{u}_{jk}^m\left|\right|D\left(k\right)-v_j|{|}^2+\mathrm{\α}\[\underset{j=1}{\overset{C}{\Σ}}\underset{k=1}{\overset{N}{\Σ}}{u}_{jk}^m||open\left(k\right)-v_j|{|}^2+\underset{j=1}{\overset{C}{\Σ}}\underset{k=1}{\overset{N}{\Σ}}{u}_{jk}^m\left|\right|close\left(k\right)-v_j|{|}^2\]$

式中，C是聚类数目，C为2，N是样本的总数，表示第k样本对于第j类聚类中心的模糊隶属度，m为隶属度的加权指数，m为2，u_jk∈[0,1]且参数α控制惩罚效果，

D(k)＝[D₁(k)D₂(k)…D_B(k)]，

open(k)＝[Π((γ^T(D₁))(k))Π((γ^T(D₂))(k))…Π((γ^T(D_B))(k))]，

采用式 $u_{jk}=\frac{{\left(\right|\left|D\right(k)-v_j|{|}^2+\mathrm{\α}\left(\right|\left|open\right(k)-v_j|{|}^2+\left|\right|close\left(k\right)-v_j\left|{|}^2\right))}^{-1/(m-1)}}{\underset{j=1}{\overset{C}{\Σ}}{\left(\right|\left|D\right(k)-v_j|{|}^2+\mathrm{\α}\left(\right|\left|open\right(k)-v_j|{|}^2+\left|\right|close\left(k\right)-v_j\left|{|}^2\right))}^{-1/(m-1)}}$ 更新模糊划分矩阵；

采用式 $v_j=\frac{\underset{k=1}{\overset{N}{\Σ}}{u}_{jk}^m\left(D\right(k)+\mathrm{\α}(open\left(k\right)+close\left(k\right)\left)\right)}{(1+\mathrm{\α})\underset{k=1}{\overset{N}{\Σ}}{u}_{jk}^m}$ 更新聚类中心；

重复更新模糊划分矩阵和聚类中心，直至相邻两次迭代的聚类中心聚类小于阈值ε；

根据模糊划分矩阵u_jk确定最终的变化检测图，确定影像的变化区域和非变化区域。