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

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

背景技术：

随着多时相遥感数据的不断积累以及空间数据库的相继建立，如何从这些遥感数据中提取和检测变化信息已成为遥感科学和地理信息科学的重要研究课题。根据同一区域不同时相的遥感影像，可以提取城市、环境等动态变化的信息，为资源管理与规划、环境保护等部门提供科学决策的依据。

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

针对上述问题，有必要在检测过程中加入空间约束来提高检测精度，另一方面，对关键参数的自适应选择也是检测技术中的重点。

技术实现要素：

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

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

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

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

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

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

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

(5)计算D_X的模值X_M，利用贝叶斯原理，并基于最大期望算法获取最优分割阈值T，将|X_M-T|≤δ的区域作为伪训练样本集，其中δ取X_M动态范围的百分比；

(6)设定置信水平的搜索范围和搜索步长；

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

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

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

(10)判断M₁中的非变化区域相较于M₀中的非变化区域是否有改变，如果两者没有改变，将给定的置信水平A加上步骤(6)设定的搜索步长作为新的置信水平，并返回步骤(7)；如果两者改变了，则将M₁中的非变化区域作为新的非变化区域，返回步骤(7)；

(11)当置信水平达到步骤(6)设定的搜索范围的上界，终止循环迭代；

(12)针对每一个置信水平，计算伪训练样本集的精度，选择精度最高的伪训练样本集对应的置信水平，并在此基础上，输出最终的变化检测结果。

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

(a)假设X_M影像上的非变化类ω_n和变化类ω_c服从如下的高斯分布：

ω_i∈{ω_n,ω_c}，m_i∈{m_n,m_c}，σ_i∈{σ_n,σ_c}上式中，m_n和σ_n分别为非变化类的均值和方差，m_c和σ_c分别为变化类的均值和方差；

(b)采用最大期望算法算法估计m_n、σ_n、m_c和σ_c这四个参数；

(c)依据贝叶斯最小误差准则，求解X_M的分割阈值T，在高斯分布的情况下，等同于求解下式：

上式中，为分割阈值T的估计值，p(ω_n)、p(ω_c)分别为非变化类先验概率和变化类先验概率；

(d)根据估计的分割阈值，构建伪训练样本集：

非变化区域伪训练样本集：

变化区域伪训练样本集：

其中，δ取X_M动态范围的15％。

进一步地，在步骤(6)中，设定置信水平的搜索范围为0.95-0.999，搜索步长为0.001。

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

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

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

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

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

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

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

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

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

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

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

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

(3)在本发明中，最优置信水平是基于伪训练样本集获取，使得该技术通用性较好。

附图说明

图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的模值X_M，利用贝叶斯原理，并基于最大期望算法获取最优分割阈值T，将|X_M-T|≤δ的区域作为伪训练样本集，其中δ取X_M动态范围的百分比。具体过程如下：

(a)假设X_M影像上的非变化类ω_n和变化类ω_c服从如下的高斯分布：

ω_i∈{ω_n,ω_c}，m_i∈{m_n,m_c}，σ_i∈{σ_n,σ_c}上式中，m_n和σ_n分别为非变化类的均值和方差，m_c和σ_c分别为变化类的均值和方差；

(b)采用最大期望算法算法估计m_n、σ_n、m_c和σ_c这四个参数，下面以非变化类的参数估计为例进行说明，变化类参数估计类似：

上式中，I和J分别表示X_M的行数和列数，上标t表示迭代次数；

(c)依据贝叶斯最小误差准则，求解X_M的分割阈值T，在高斯分布的情况下，等同于求解下式：

上式中，为分割阈值T的估计值，p(ω_n)、p(ω_c)分别为非变化类先验概率和变化类先验概率；

(d)根据估计的分割阈值，构建伪训练样本集：

非变化区域伪训练样本集：

变化区域伪训练样本集：

其中，δ取X_M动态范围的15％。

步骤6：设定置信水平的搜索范围为0.95-0.999，以及置信水平的搜索步长0.001。

步骤7：确定D_X中的非变化区域(初始化时，将整个D_X作为非变化区域)，计算非变化区域的均值矢量m和方差矩阵Σ，并计算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的波段数目。

步骤8：在给定的置信水平A的基础上，计算检测阈值，并根据该阈值进行卡方检测，获得初步的检测结果M₀。检测阈值的计算方法如下：

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

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

步骤10：判断M₁中的非变化区域相较于M₀中的非变化区域是否有改变，如果两者没有改变，将给定的置信水平A加上步骤6设定的搜索步长作为新的置信水平，并返回步骤7；如果两者改变了，则将M₁中的非变化区域作为新的非变化区域，返回步骤7。

步骤11：当置信水平达到0.999，终止循环迭代。

步骤12：针对每一个置信水平，计算伪训练样本集的精度，选择精度最高的伪训练样本集对应的置信水平，并在此基础上，输出最终的变化检测结果。

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

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

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

(1)基于CVA的EM方法(CVA-EM)[意大利的Bruzzone L.等在文章“Automatic analysis of difference image for unsupervised change detection”(IEEE Transactions 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 Sensing of Environment,2006,102(1-2):1-12.)中所提的方法]。

(3)基于鲁棒估计的CST检测(RCST)方法[Aiye Shi等在文章“Unsupervised change detection based on robust chi-squared transform for bitemporal remotely 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

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