CN109255781B - 一种面向对象的多光谱高分辨率遥感影像变化检测方法 - Google Patents

Abstract

本发明公开了一种面向对象的基于CRF的多时相多光谱遥感影像变化检测方法，包括以下步骤：输入同一区域、不同时相的两幅高分辨率光学遥感影像；利用遥感软件ENVI对X₁和X₂进行影像配准，分为粗校正和精校；利用多元变化检测MAD方法对X₁和X₂进行辐射归一化校正；基于SLIC算法进行分割；计算多时相影像的差异影像，利用FCM算法进行聚类；根据SLIC的分割图像构建CRF的多邻域系统，构建CRF的二元项能量；根据循环信度传播LBP优化算法对CRF能量项U进行优化，获取最终的变化检测结果。本发明使得变化检测的结果更加可靠，也更加具有稳健性。

Description

一种面向对象的多光谱高分辨率遥感影像变化检测方法

技术领域

本发明涉及一种面向对象的基于条件随机场CRF的多时相多光谱高分辨率遥感影像非监督变化检测方法，属于光学遥感影像变化检测技术领域。

背景技术

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

遥感影像的变化检测就是从不同时期的遥感数据中，定量地分析和确定地表变化的特征与过程。各国学者从不同的角度和应用研究提出了许多有效的检测算法，总体来说，根据检测过程中是否需要训练样本，可以将变化检测分为三个大类：非监督变化检测算法、半监督变化检测算法和监督变化检测算法。

非监督变化检测一般先通过构造差分影像，进而再通过预先设定的阈值将差分影像分为变化和非变化区域。传统的非监督变化检测主要有算法运算法、CVA法及其改进的算法、主成分变换(Principle Component Analysis，PCA)法、影像回归法等。上述传统的变化检测方法大多是基于像素的方法，适用于中-低分辨率的遥感影像，且这些方法都是暗含像素在空间上具有独立性。针对高分辨率影像，由于其结构、纹理等信息更加突出，面向对象的变化检测更有利于知识的结合利用，也更能有效地利用高分辨率影像所具有的多特征优势。

另外一方面，变化检测建模过程中，像素和其邻域具有高度相关性，MRF和CRF两种概率图模型均可以刻画这种邻域关系。由于CRF能够同时考虑观察场与标记场的空间上下文信息，比MRF建模更加灵活，其在多时相变化检测领域应用具有优势。如李民等提出的基于非平稳分析与条件随机场的SAR图像变化检测方法(2015，中国发明专利，西安电子科技大学)，但是该发明专利技术是监督类变化检测，实际应用中需消耗大量的人力成本等去构建训练样本。Cao等提出的基于CRF的非监督变换检测方法应用于多光谱多时相遥感图像检测过程中，提高了检测的精度。但是该方法的是基于像素的，检测的效率和精度有一定的局限性。B.Desclée等提出的面向对象的基于卡方变换的检测方法(OBCT)，可以克服基于像素的变化检测的局限性，但是该面向对象的方法没有考虑对象之间的相关性。

发明内容

针对现有技术存在的不足，本发明目的是提供一种面向对象的基于条件随机场CRF的多时相多光谱高分辨率遥感影像非监督变化检测方法，本发明首先利用简单线性迭代聚类SLIC算法对多时相遥感影像进行分割，然后对差异影像的均值应用FCM算法算法获取CRF的一元能量项；其次根据图像的每一个节点(影像的对象)的邻域和差异影像的幅值，获取CRF的二元能量项；最后通过循环信度传播LBP算法对CRF的能量最小化获取最终的变化检测结果，使得变化检测的结果更加可靠，也更加具有稳健性。

为了实现上述目的，本发明是通过如下的技术方案来实现：

本发明的一种面向对象的基于CRF的多时相多光谱遥感影像变化检测方法，包括以下步骤：

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

步骤2：利用遥感软件ENVI对X₁和X₂进行影像配准，分为粗校正和精校；

步骤3：利用多元变化检测MAD方法对X₁和X₂进行辐射归一化校正；

步骤4：对经过预处理的多时相高分辨率影像进行叠合，然后基于SLIC算法进行分割；

步骤5：计算多时相影像的差异影像，对所述差异影像的每一个波段，计算每一个对象的均值；记各波段各对象均值所组成的影像为X_M，继而利用FCM算法进行聚类；

步骤6：根据SLIC的分割图像构建CRF的多邻域系统，并结合步骤5的差异影像，构建CRF的二元项能量；

步骤7：根据循环信度传播LBP优化算法对CRF能量项U进行优化，获取最终的变化检测结果。

步骤2中，对于所述粗校正，具体步骤如下：

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

2.2)采集地面控制点GCPs，所述GCPs均匀分布在整幅图像内，所述GCPs的数目至少大于等于25；

2.3)计算均方根误差；

2.4)选择多项式校正模型；

2.5)采用双线性插值进行重采样输出。

上述双线性插值具体方法如下：

求未知函数f在坐标点P＝(x,y)的值，x是横坐标，y是纵坐标；基于预设的已知函数f在Q₁₁＝(x₁,y₁),Q₁₂＝(x₁,y₂),Q₂₁＝(x₂,y₁)及Q₂₂＝(x₂,y₂)四个坐标点的值(x₁，x₂是横坐标，y₁，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 (1)

步骤2中，所述精校正即将经过几何粗校正的多光谱遥感影像数据，利用自动匹配与三角剖分法进行几何精校正。

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

步骤3中，所述多元变化检测MAD相对辐射归一化校正具体方法如下：

首先，通过典型相关分析找到两期影像各波段亮度值的一个线性组合，得到变化信息增强的差异影像，通过对归一化的MAD变量平方和取阈值确定变化和未变化区域；

然后，通过所述未变化区域对应的两时相像元对的映射方程，完成相对辐射校正。

步骤4，具体实现步骤如下：

4.1)初始化种子点即聚类中心：按照设定的超像素个数，在图像内均匀的分配种子点；基于预设的图片总共有N个像素点，预分割为K个相同尺寸的超像素，那么每个超像素的大小为N/K，则相邻种子点的距离即步长为

4.2)在种子点的z×z邻域内重新选择种子点；

4.3)在每个种子点周围的邻域内为每个像素点分配类标签，搜索范围限制为2S×2S；

4.4)对于每个搜索到的像素点，分别计算它和种子点的颜色距离和空间距离；取最小值对应的种子点作为像素点的聚类中心；

4.5)对步骤4.1)-步骤4.4)进行迭代直到误差收敛。

步骤4.2)中，所述种子点具体选择方法为：

计算邻域内所有像素点的梯度值，将种子点移到邻域内梯度最小的地方。

步骤5中，基于预设的加权指数m＝2，聚类数目C＝2，具体方法如下：

5.1)基于预设的构建FCM：

式中，C是聚类数目，K是对象样本的总数，k表示样本的序号，j表示类别中心的序号，X_M表示多时相差异影像各波段所有对象均值所组成的影像，v_j表示第j聚类中心，d_kj表示第k个样本和第j聚类中心的欧式距离，

表示第k样本对于第j类聚类中心的模糊隶属度，m为隶属度的加权指数，u_jk∈[0,1]且

5.2)公式(2)的目标函数最小化用下述的公式交替进行：

式中，d_ki表示第k个样本和第i聚类中心的欧式距离、d_kj表示第k个样本和第j聚类中心的欧式距离；

5.3)FCM算法收敛后，获得CRF的一元能量项U_j＝u_jk。

步骤6)中，构建CRF的二元项能量如下：

式中，i,j表示图像相邻的两个节点即相邻的对象(本发明中，将每一个对象看作图模型中的一个节点)；c_i,c_j表示两个节点i,j处的类别标签；两种类别为变化和非变化；D_M(i)表示对象i的内部所有差异影像的幅值均值；L(i，j)表示两个对象的公共边界的长度；[c_i≠c_j]表示0-1指示函数；

最后构建的CRF能量项如下：

式中，U表示CRF的总能量项，U_j、U_i,j分别表示对象j或节点j的一元能量项、相邻对象j和i间的二元能量项，S表示节点的集合，E表示相邻节点的集合，参数β是正则化参数，用于控制一元能量项和二元能量项间的比重。

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

(1)利用面向对象技术实现变化检测，可以消除变化检测结果中的“椒盐”噪声，使得变化检测的结果更加可靠，也更加具有稳健性。

(2)利用CRF对图像的相邻对象空间关系进行刻画，提高了变化检测的精度。

附图说明

图1是本发明面向对象的基于CRF的非监督多时相多光谱高分辨率遥感影像变化检测方法的实现流程示意图；

图2(a)是本发明所采用的2007年1月的沙特阿拉伯Mina地区高分辨率IKONOS图像第3波段示意图；

图2(b)是本发明所采用的2007年12月的沙特阿拉伯的Mina地区高分辨率IKONOS图像第3波段示意图；

图2(c)是变化检测参考(Ground Truth)图像；

图3(a)EM-MRF算法检测结果图像；

图3(b)CRF算法检测结果图像；

图3(c)OBCT算法的检测结果图像；

图3(d)本发明算法检测结果图像。

具体实施方式

为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。

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

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

步骤2：利用遥感软件ENVI对X₁和X₂进行影像配准，分为粗校正和精校正两个步骤：

对于几何粗校正，利用ENVI4.8软件中的相关功能实现，具体操作步骤为：(1)显示基准影像和待校正影像；(2)采集地面控制点GCPs；GCPs应均匀分布在整幅图像内，GCPs的数目至少大于等于25。(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 (1)

对于几何精校正，将经过几何粗校正的多光谱遥感影像数据，利用自动匹配与三角剖分法进行几何精校正。

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

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

步骤4：对经过预处理的多时相高分辨率影像进行叠合，然后基于SLIC算法进行分割。具体实现步骤如下：

4.1)初始化种子点(聚类中心)：按照设定的超像素个数，在图像内均匀的分配种子点。假设图片总共有N个像素点，预分割为K个相同尺寸的超像素，那么每个超像素的大小为N/K，则相邻种子点的距离(步长)近似为

4.2)在种子点的n×n邻域内重新选择种子点(一般取n＝3)。具体方法为：计算该邻域内所有像素点的梯度值，将种子点移到该邻域内梯度最小的地方。

4.3)在每个种子点周围的邻域内为每个像素点分配类标签，搜索范围限制为2S×2S。

4.4)对于每个搜索到的像素点，分别计算它和该种子点的颜色距离和空间距离。取最小值对应的种子点作为该像素点的聚类中心。

4.5)对步骤4.1-4.4进行迭代直到误差收敛。

4.6)对出现多连通情况、超像素尺寸过小，单个超像素被切割成多个不连续超像素等，通过增强连通性解决。

步骤5：计算多时相影像的差异影像，对差异影像的每一个波段，计算每一个对象的均值。记各波段各对象均值所组成的影像为X_M，继而利用FCM算法进行聚类，设定隶属度的加权指数m＝2，聚类数目C＝2，具体实施如下：

3)首先构建FCM的模型如下：

式中，C是聚类数目，K是对象样本的总数，k表示样本的序号，j表示类别中心的序号，X_M表示多时相差异影像各波段所有对象均值所组成的影像，v_j表示第j聚类中心，d_kj表示第k个样本和第j聚类中心的欧式距离，

表示第k样本对于第j类聚类中心的模糊隶属度，m为隶属度的加权指数，u_jk∈[0,1]且

2)式(2)的目标函数最小化可以用下述的公式交替进行：

式中，d_ki表示第k个样本和第i聚类中心的欧式距离、d_kj表示第k个样本和第j聚类中心的欧式距离。

3)FCM算法收敛后，获得CRF的一元能量项U_j＝u_jk。

步骤6：根据SLIC的分割图像构建CRF的4邻域系统，并结合步骤5的差异影像，构建CRF的二元项能量如下：

式中，i,j表示图像相邻的两个节点(对象)；c_i,c_j表示两个节点i,j处的类别标签(两种类别：变化和非变化)；D_M(i)对象i的内部所有差异影像的幅值均值；L(i，j)表示两个对象的公共边界的长度；[c_i≠c_j]表示0-1指示函数。最后构建的CRF能量项如下：

式中S表示节点的集合，E表示相邻节点的集合，参数β是正则化参数，用于控制一元能量项和二元能量项间的比重。

步骤7：根据优化算法循环信度传播(Loopy Belief Propagation,LBP)算法对CRF能量项U进行优化，获取最终的变化检测结果。

本发明的实验数据为沙特阿拉伯的Mina地区的多时相IKNOS高分辨影像数据，图像大小为700×950，使用B1、B2和B3三个波段。为了验证本发明的有效性，将本发明变化检测方法与下述变化检测方法进行比对：

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

(2)Cao等所提的CRF检测(CRF)方法[Guo Cao,Xuesong Li&LicunZhou.Unsupervised change detection in high spatial resolution remote sensingimages based on a conditional random field model.European Journal of RemoteSensing,2016,49:225-237.]

(3)B.Desclée等提出的面向对象的基于卡方变换的检测方法(OBCT)[B.Desclée,P.Bogaert,and P.Defourny,“Forest change detection by statistical object-basedmethod,”Remote Sens.Environ,2006,102(1–2),1–11.]

(4)本发明方法。

检测性能用错检数FP、漏检数FN、总错误数OE和Kappa系数四个指标来衡量。FP、FN和OE越接近于0、Kappa系数越接近于1，表明变化检测方法的性能越好。检测结果如表1所示。

表1 Mina地区的多时相IKONOS影像变化检测结果比较

由表1可见，本发明所提的检测方法Kappa系数最大，比其他三种检测算法相比更接近于1。另外，本发明的总错误数OE在对比算法中是最小的，更接近于0。综上，本发明变化检测算法的性能优于其他三种检测方法，这表明本发明所提的变化检测方法是有效的。

图2(a)是Mina地区的前一时相多光谱IKONOS影像，图2(b)是Mina地区的后一时相多光谱IKONOS影像，图2(c)是变化检测的参考图。图3(a)是EM-MRF算法的变化检测结果，图3(b)是CRF算法的变化检测结果，图3(c)是OBCT算法的变化检测结果，图3(d)是本发明方法的变化检测结果。从图2(c)的参考图和图3(a)-(d)的对比来看，目视效果上，本发明所提算法的检测效果是最好的。

以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (10)

1.一种面向对象的基于CRF的多时相多光谱遥感影像变化检测方法，其特征在于，包括以下步骤：

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

步骤2：利用遥感软件ENVI对X₁和X₂进行影像配准，分为粗校正和精校；

步骤3：利用多元变化检测MAD方法对X₁和X₂进行辐射归一化校正；

步骤4：对经过预处理的多时相高分辨率影像进行叠合，然后基于SLIC算法进行分割；

步骤5：计算多时相影像的差异影像，对所述差异影像的每一个波段，计算每一个对象的均值；记各波段各对象均值所组成的影像为X_M，继而利用FCM算法进行聚类；

步骤6：根据SLIC的分割图像构建CRF的多邻域系统，并结合步骤5的差异影像，构建CRF的二元项能量；

步骤7：根据循环信度传播LBP优化算法对CRF能量项U进行优化，获取最终的变化检测结果。

2.根据权利要求1面向对象的基于CRF的多时相多光谱遥感影像变化检测方法，其特征在于，步骤2中，对于所述粗校正，具体步骤如下：

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

2.2)采集地面控制点GCPs，所述GCPs均匀分布在整幅图像内，所述GCPs的数目至少大于等于25；

2.3)计算均方根误差；

2.4)选择多项式校正模型；

2.5)采用双线性插值进行重采样输出。

3.根据权利要求2面向对象的基于CRF的多时相多光谱遥感影像变化检测方法，其特征在于，所述双线性插值具体方法如下：

求未知函数f在坐标点P＝(x,y)的值，x和y分别表示横坐标和纵坐标；基于预设的已知函数f在Q₁₁＝(x₁,y₁),Q₁₂＝(x₁,y₂),Q₂₁＝(x₂,y₁)及Q₂₂＝(x₂,y₂)四个坐标点的值，x₁，x₂是横坐标，y₁，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 (1)。

4.根据权利要求1面向对象的基于CRF的多时相多光谱遥感影像变化检测方法，其特征在于，步骤2中，所述精校正即将经过几何粗校正的多光谱遥感影像数据，利用自动匹配与三角剖分法进行几何精校正。

5.根据权利要求4面向对象的基于CRF的多时相多光谱遥感影像变化检测方法，其特征在于，所述三角剖分法即采用逐点插入法构建Delaunay三角网，对每一个三角形，利用其三个顶点的行列号与其对应的基准影像同名点的地理坐标来确定该三角形内部的仿射变换模型参数，对待校正影像进行纠正，得到校正后的遥感影。

6.根据权利要求1面向对象的基于CRF的多时相多光谱遥感影像变化检测方法，其特征在于，步骤3中，所述多元变化检测MAD相对辐射归一化校正具体方法如下：

首先，通过典型相关分析找到两期影像各波段亮度值的一个线性组合，得到变化信息增强的差异影像，通过对归一化的MAD变量平方和取阈值确定变化和未变化区域；

然后，通过所述未变化区域对应的两时相像元对的映射方程，完成相对辐射校正。

7.根据权利要求1面向对象的基于CRF的多时相多光谱遥感影像变化检测方法，其特征在于，步骤4，具体实现步骤如下：

4.1)初始化种子点即聚类中心：按照设定的超像素个数，在图像内均匀的分配种子点；基于预设的图像总共有N个像素点，预分割为K个相同尺寸的超像素，那么每个超像素的大小为N/K，则相邻种子点的距离即步长为

4.2)在种子点的z×z邻域内重新选择种子点；

4.3)在每个种子点周围的邻域内为每个像素点分配类标签，搜索范围限制为2S×2S；

4.4)对于每个搜索到的像素点，分别计算它和种子点的颜色距离和空间距离；取最小值对应的种子点作为像素点的聚类中心；

4.5)对步骤4.1)-步骤4.4)进行迭代直到误差收敛。

8.根据权利要求7面向对象的基于CRF的多时相多光谱遥感影像变化检测方法，其特征在于，步骤4.2)中，所述种子点具体选择方法为：

计算邻域内所有像素点的梯度值，将种子点移到邻域内梯度最小的地方。

9.根据权利要求1面向对象的基于CRF的多时相多光谱遥感影像变化检测方法，其特征在于，步骤5中，基于预设的加权指数m＝2，聚类数目C＝2，具体方法如下：

5.1)基于预设的构建FCM：

式中，C是聚类数目，K是对象样本的总数，k表示样本的序号，j表示类别中心的序号，X_M表示多时相差异影像各波段所有对象均值所组成的影像，v_j表示第j聚类中心，d_kj表示第k个样本和第j聚类中心的欧式距离，

表示第k个样本对于第j类聚类中心的模糊隶属度，m为隶属度的加权指数，u_jk∈[0,1]且

5.2)公式(2)的目标函数最小化用下述的公式交替进行：

其中，d_ki表示第k个样本和第i聚类中心的欧式距离、d_kj表示第k个样本和第j聚类中心的欧式距离；

5.3)FCM算法收敛后，获得CRF的一元能量项U_j＝u_jk。

10.根据权利要求1面向对象的基于CRF的多时相多光谱遥感影像变化检测方法，其特征在于，步骤6)中，构建CRF的二元项能量如下：

式中，i,j表示图像相邻的两个对象或者节点；c_i,c_j表示两个节点i,j处的类别标签；两种类别为变化和非变化；D_M(i)表示对象i的内部所有差异影像的幅值均值；L(i，j)表示两个对象的公共边界的长度；[c_i≠c_j]表示0-1指示函数；

最后构建的CRF能量项如下：

式中，U表示CRF的总能量项，U_j、U_i,j分别表示对象j或节点j的一元能量项、相邻对象j和i间的二元能量项，S表示节点的集合，E表示相邻节点的集合，参数β是正则化参数，用于控制一元能量项和二元能量项间的比重。

Images