CN106469452A - 基于空间约束的卡方变换的多时相遥感影像变化检测方法 - Google Patents

Abstract

本发明公开了基于空间约束的卡方变换的多时相遥感影像变化检测方法，首先计算多时相高分辨率影像的多波段差异影像，在此基础上，计算影像上每一个点的CST值，根据致信水平获取阈值，得到初步的变化检测结果，然后再对该初步结果进行众数滤波(嵌入空间信息)，并根据滤波结果重新计算非变化区域的均值和方差矩阵。重复上述过程直到检测结果没有变化为止，最后获取最终的变化检测结果。本发明解决了多时相多光谱遥感影像背景信息复杂、噪声干扰严重的问题。

Description

基于空间约束的卡方变换的多时相遥感影像变化检测方法

技术领域

本发明属于遥感影像处理技术领域，特别涉及了基于空间约束的卡方变换的多时相遥感影像变化检测方法。

背景技术

随着多时相高分辨率遥感数据的不断积累以及空间数据库的相继建立，如何从这些遥感数据中提取和检测变化信息已成为遥感科学和地理信息科学的重要研究课题。根据同一区域不同时相的遥感影像，可以提取城市、环境等动态变化的信息，为资源管理与规划、环境保护等部门提供科学决策的依据。我国“十二五”将加大拓展实施“十一五”已启动实施的高分辨率对地观测工程，关注包括高分辨率遥感目标与空间环境特征分析及高可靠性自动解译等基础理论与关键技术研究，正在成为解决国家安全和社会经济发展重大需求的研究焦点。

遥感影像的变化检测就是从不同时期的遥感数据中，定量地分析和确定地表变化的特征与过程。各国学者从不同的角度和应用研究提出了许多有效的检测算法，如变化矢量分析法(Change Vector Analysis，CVA)、基于Fuzzy C-means(FCM)的聚类方法等。其中，传统的基于卡方变换(Chi-Squared Transform,CST)的多时相光学遥感变化检测，先计算差异影像的均值和方差矩阵，然后再基于置信水平，确定变化检测的阈值，进而得到变化检测结果。该类技术中，使用CST的不足是仅使用多时相高分辨率差异影像的光谱信息，没有利用空间信息。另外，在计算差异影像的均值和方差矩阵时，估计的精度不高。

针对上述问题，有必要研究新的高分辨率可见光遥感图像变化检测技术来有效克服上述难点。

发明内容

为了解决上述背景技术提出的技术问题，本发明旨在提供基于空间约束的卡方变换的多时相遥感影像变化检测方法，采用基于空间约束的卡方变换的变化检测方法，克服了多时相多光谱遥感影像背景信息复杂、噪声干扰严重的问题。

为了实现上述技术目的，本发明的技术方案为：

基于空间约束的卡方变换的多时相遥感影像变化检测方法，包括以下步骤：

(1)输入两时相的高分辨率光学遥感影像，分别记为X₁和X₂；

(2)对X₁和X₂进行影像配准；

(3)利用多元变化检测方法分别对X₁和X₂进行辐射归一化校正；

(4)计算多时相差异影像D_X＝X₁-X₂；

(5)确定D_X上的非变化区域，计算非变化区域的均值和方差矩阵，并计算D_X上每个点的卡方值；初始化时，将整个D_X作为非变化区域；

(6)给定置信水平，计算检测阈值，并根据该阈值进行卡方检测，获得初步的检测结果M₀；

(7)给定滤波窗口尺寸，对M₀进行众数滤波，即滤波窗口内，如果变化的像素数目大于非变化的，则窗口中心点像素为变化的，反之，窗口中心点像素为非变化的，记众数滤波的结果为M₁；

(8)判断M₁中的非变化区域相较于M₀中的非变化区域是否有改变，如果两者没有改变，则将M₁作为最终的变化检测结果；如果两者改变了，则将M₁中的非变化区域作为新的非变化区域，返回步骤(5)，循环迭代。

进一步地，在步骤(7)中，给定滤波窗口的尺寸为3×3或者5×5。

进一步地，在步骤(6)中，采用下式计算检测阈值：

上式中，1-α是置信水平，C_ij表示D_X在坐标(i，j)处的卡方值，为检测阈值。

进一步地，在步骤(2)中，对X₁和X₂进行影像配准包括几何粗校正和几何精校正，所述几何粗校正的过程：

(A)选择X₁和X₂分别作为基准影像和待校正影像；

(B)在基准影像和待校正影像上分别采集地面控制点，地面控制点的数量大于等于9，且地面控制点均匀分布在影像上；

(C)计算基准影像和待校正影像各地面控制点处的均方误差；

(D)采用多项式纠正法对待校正影像进行纠正；

(E)采用双线性插值法对待校正影像进行重采样；

所述几何精校正是将经过几何粗校正的遥感影像，利用自动匹配与三角剖分法进行校正。

进一步地，步骤(3)的具体过程如下：

首先找到两期影像X₁和X₂各波段亮度值的一个线性组合，得到变化信息增强的差异影像，根据设定的阈值确定该差异影像的变化区域和未变化区域，然后通过未变化区域对应的两时相像元对的映射方程，完成相对辐射校正。

采用上述技术方案带来的有益效果：

(1)本发明将众数滤波引入到CST检测结果之后，使得使用CST变换具有空间约束能力；

(2)本发明采用迭代的方法，估计非变化区域的均值和标准差，克服了仅仅利用整个多波段影像Dx计算均值和方差矩阵的不足，可以使变化检测的结果更加可靠，也更加具有稳健性。

附图说明

图1是本发明的基本流程图；

图2(a)、2(b)分别是2007年1月的沙特***Mina地区高分辨率IKONOS图像第3波段示意图、2007年12月的沙特***的Mina地区高分辨率IKONOS图像第3波段示意图；

图3是变化检测的参考图像；

图4(a)、4(b)、4(c)、4(d)分别是EM-CVA算法、ICST算法、RCST算法、本发明算法的检测结果示意图。

具体实施方式

以下将结合附图，对本发明的技术方案进行详细说明。

如图1所示，基于空间约束的卡方变换的多时相遥感影像变化检测方法，具体步骤如下：

步骤1：输入同一区域、不同时相的两幅高分辨率光学遥感影像，分别记为X₁和X₂。

步骤2：对X₁和X₂进行影像配准，分为粗校正和精校正两个步骤：

对于几何粗校正，利用ENVI4.8软件中的相关功能实现，具体操作步骤为：

(1)显示基准影像和待校正影像；

(2)采集地面控制点GCPs，GCPs应均匀分布在整幅图像内，GCPs的数目至少大于等于9；

(3)计算均方误差；

(4)采用多项式纠正法对待校正影像进行纠正；

(5)采用双线性插值进行重采样输出，若求未知函数f在点P＝(x,y)的值，假设已知函数f在Q₁₁＝(x₁,y₁),Q₁₂＝(x₁,y₂),Q₂₁＝(x₂,y₁),及Q₂₂＝(x₂,y₂)四个点的值，如果选择一个坐标***使得这四个点的坐标分别为(0,0)、(0,1)、(1,0)和(1,1)，那么双线性插值公式就可以表示为：

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

对于几何精校正，将经过几何粗校正的多光谱遥感影像数据，利用自动匹配与三角剖分法进行几何精校正。采用逐点***法构建Delaunay三角网，对每一个三角形，利用其三个顶点的行列号与其对应的基准影像同名点的地理坐标来确定该三角形内部的仿射变换模型参数，对待校正影像进行纠正，得到校正后的遥感影。

步骤3：利用多元变化检测(Multivariate Alteration Detection，MAD)方法对X₁和X₂进行辐射归一化校正。首先找到X₁和X₂各波段亮度值的一个线性组合，得到变化信息增强的差异影像，通过阈值确定变化和未变化区域，然后通过未变化区域对应的两时相像元对的映射方程，完成相对辐射校正。

步骤4：对输入的多时相高分辨率影像X₁和X₂，计算多时相差异影像D_X：

D_X＝X₁-X₂

步骤5：确定D_X上的非变化区域(初始化时，将整个D_X作为非变化区域)，计算非变化区域的均值和方差矩阵，并计算D_X上每个点的卡方值：

C_ij＝(x_ij-m)^TΣ^-1(x_ij-m)～χ²(b)

上式中，C_ij表示D_X在(i，j)坐标点的卡方值，其服从自由度为b的卡方分布；

x_ij表示D_X在(i，j)坐标点的矢量值；Σ^-1表示方差矩阵的逆矩阵；b表示D_X的波段数目。

步骤6：给定置信水平，计算检测阈值，并根据该阈值进行卡方检测，获得初步的检测结果M₀。检测阈值的计算公式如下：

当置信水平为1-α时，C_ij的值大于的概率为α。如果α取值较小，则大于的C_ij可以视为出界点(outlier)或者变化点，由此确定检测阈值为

步骤7：给定滤波窗口尺寸(一般设置为3×3或者5×5的窗口)，对M₀进行众数滤波，即滤波窗口内，如果变化的像素数目大于非变化的，则窗口中心点像素为变化的，反之，窗口中心点像素为非变化的，记众数滤波的结果为M₁。

步骤8：判断M₁中的非变化区域相较于M₀中的非变化区域是否有改变，如果两者没有改变，则将M₁作为最终的变化检测结果；如果两者改变了，则将M₁中的非变化区域作为新的非变化区域，返回步骤(5)，循环迭代。

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

本发明的实验数据为沙特***的Mina地区的多时相IKNOS高分辨影像数据，图像大小为700×950，使用B1、B2和B3三个波段，图2(a)和图2(b)为两时像B3波段的遥感影像。

为了验证本发明的有效性，将本发明变化检测方法与下述变化检测方法进行比对：

(1)基于CVA的EM方法(CVA-EM)[意大利的Bruzzone L.等在文章“Automaticanalysis of difference image for unsupervised change detection”(IEEETransactions on Geoscience and Remote Sensing,2000,38(3):1171-1182.)中所提的检测方法]。

(2)基于迭代的CST检测(ICST)方法[B.Desclée,P.Bogaert,and P.Defourny在文章“Forest change detection by statistical object-based method”(Remote Sensingof Environment,2006,102(1-2):1-12.)中所提的方法]。

(3)基于鲁棒估计的CST检测(RCST)方法[Aiye Shi等在文章“Unsupervisedchange detection based on robust chi-squared transform for bitemporalremotely sensed images”(International of Remote Sensing,2006,102(1-2):1-12.)中所提的方法]。

(4)本发明方法。

图3为变化检测的参考图像。上述四种方法的检测结果如图4(a)、4(b)、4(c)、4(d)所示，检测性能用错检数FP、漏检数FN、总错误数OE和Kappa系数四个指标来衡量。FP、FN和OE越接近于0、Kappa系数越接近于1，表明变化检测方法的性能越好。检测结果如表1所示。由表1可见，本发明所提的检测方法性能优于其他三种检测方法，这表明本发明所提的变化检测方法是有效的。

表1

方法	FP	FN	OE	k
					CVA-EM	49065	1456	50521	0.556
ICST	11236	868	12104	0.850
					RCST	2192	4346	6538	0.906
本发明方法	4837	1602	6439	0.913
					理想	0	0	0	1

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (5)

1.基于空间约束的卡方变换的多时相遥感影像变化检测方法，其特征在于，包括以下步骤：

(1)输入两时相的高分辨率光学遥感影像，分别记为X₁和X₂；

(2)对X₁和X₂进行影像配准；

(3)利用多元变化检测方法分别对X₁和X₂进行辐射归一化校正；

(4)计算多时相差异影像D_X＝X₁-X₂；

(5)确定D_X上的非变化区域，计算非变化区域的均值和方差矩阵，并计算D_X上每个点的卡方值；初始化时，将整个D_X作为非变化区域；

(6)给定置信水平，计算检测阈值，并根据该阈值进行卡方检测，获得初步的检测结果M₀；

(7)给定滤波窗口尺寸，对M₀进行众数滤波，即滤波窗口内，如果变化的像素数目大于非变化的，则窗口中心点像素为变化的，反之，窗口中心点像素为非变化的，记众数滤波的结果为M₁；

(8)判断M₁中的非变化区域相较于M₀中的非变化区域是否有改变，如果两者没有改变，则将M₁作为最终的变化检测结果；如果两者改变了，则将M₁中的非变化区域作为新的非变化区域，返回步骤(5)，循环迭代。

2.根据权利要求1所述基于空间约束的卡方变换的多时相遥感影像变化检测方法，其特征在于：在步骤(7)中，给定滤波窗口的尺寸为3×3或者5×5。

3.根据权利要求1所述基于空间约束的卡方变换的多时相遥感影像变化检测方法，其特征在于：在步骤(6)中，采用下式计算检测阈值：

$P(C_{ij}<{\mathrm{\&chi;}}_{1-\mathrm{\&alpha;}}^2(b\left)\right)=1-\mathrm{\&alpha;}$

上式中，1-α是置信水平，C_ij表示D_X在坐标(i，j)处的卡方值，为检测阈值。

4.根据权利要求1所述基于空间约束的卡方变换的多时相遥感影像变化检测方法，其特征在于：在步骤(2)中，对X₁和X₂进行影像配准包括几何粗校正和几何精校正，所述几何粗校正的过程：

(A)选择X₁和X₂分别作为基准影像和待校正影像；

(B)在基准影像和待校正影像上分别采集地面控制点，地面控制点的数量大于等于9，且地面控制点均匀分布在影像上；

(C)计算基准影像和待校正影像各地面控制点处的均方误差；

(D)采用多项式纠正法对待校正影像进行纠正；

(E)采用双线性插值法对待校正影像进行重采样；

所述几何精校正是将经过几何粗校正的遥感影像，利用自动匹配与三角剖分法进行校正。

5.根据权利要求1所述基于空间约束的卡方变换的多时相遥感影像变化检测方法，其特征在于：步骤(3)的具体过程如下：

首先找到两期影像X₁和X₂各波段亮度值的一个线性组合，得到变化信息增强的差异影像，根据设定的阈值确定该差异影像的变化区域和未变化区域，然后通过未变化区域对应的两时相像元对的映射方程，完成相对辐射校正。