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

Abstract

本发明公开了一种面向对象的多光谱高分辨率遥感影像变化检测方法，包括以下步骤：输入同一区域、不同时相的两幅高分辨率光学遥感影像；利用遥感软件ENVI对X₁和X₂进行影像配准，分为粗校正和精校正；利用多元变化检测MAD方法对X₁和X₂进行辐射归一化校正；基于SLIC算法进行分割；计算多时相影像的差异影像，利用K‑Means算法进行聚类，依据聚类结果并结合SLIC的分割对象确定CRF一元能量项；根据SLIC的分割图像构建CRF的邻域***，并根据差异影像和对象的空间坐标，构建CRF的二元项能量；对CRF能量项U进行优化，获取最终的变化检测结果。本发明不仅考虑相邻对象的光谱差异，还考虑相邻对象的空间位置关系，提高了变化检测的精度。

Description

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

技术领域

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

背景技术

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

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

基于像素的非监督变化检测方法，适用于中-低分辨率的遥感影像，且这些方法都是暗含像素在空间上具有独立性。针对高分辨率影像，由于其结构、纹理等信息更加突出，利用传统的基于像元光谱统计的处理方式，其效率及其获得的结果信息都是十分有限的，且处理结果往往会存在许多小斑块，难以对其分析和描述。而面向对象的变化检测更有利于知识的结合利用，也更能有效地利用高分辨率影像所具有的多特征优势。

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

发明内容

针对现有技术存在的不足，本发明目的是提供一种面向对象的多光谱高分辨率遥感影像变化检测方法，将基于像素和面向对象技术二者结合，使得变化检测结果更稳健，在构建CRF二元能量项时，不仅考虑相邻对象的光谱差异，还考虑相邻对象的空间位置关系，提高了变化检测的精度。

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

本发明的一种面向对象的多光谱高分辨率遥感影像变化检测方法，包括以下步骤：

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

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

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

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

步骤5：计算多时相影像的差异影像，记各波段差异影像所组成的影像为X_M，然后利用K-Means算法进行聚类，根据聚类结果，确定CRF一元能量项；

步骤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)的值，基于预设的已知函数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)。

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

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

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

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

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

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

先将经过预处理的两时相多光谱高分辨率遥感影像进行叠合。然后执行如下步骤：

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

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

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

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

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

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

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

步骤5中，基于预设的聚类数目K＝2，所述CRF一元能量项具体确定方法如下：

5.1)随机选取K个聚类中心；

5.2)将每个数据点归类到离它最近的那个中心点所代表的类别中；

5.3)重新计算已经得到的各个类的中心；

5.4)迭代步骤5.2)～步骤5.3)，直至新的中心与原中心相等或小于设定阈值；

5.5)K-Means算法收敛后，统计每一个对象内部属于不同类别的个数，记对象j属于非变化类的个数为N_j1，属于变化类的个数为N_j2，则对象j的一元能量项为：

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

式中，i,j表示图像相邻的两个对象或者节点(本发明中，将每一个对象看作图模型中的一个节点)；c_i,c_j表示两个节点i,j处的类别标签；两种类别为变化和非变化；L(i，j)表示两个对象的公共边界的长度；[c_i≠c_j]表示0-1指示函数；F_M(i)、F_M(j)分别表示由对象i、j内部差异影像值的均值、坐标的均值构成的矢量，其中F_M(i)的构造如下：

F_M(i)＝[Mean(Dif_i)Mean(A_ix)Mean(A_iy)] (4)

式中，Mean(Dif_i)表示差异影像的幅值在对象i内的均值，Dif表示对象i内差异影像值的均值，Mean(A_ix)表示对象i内所有像素的横坐标的均值，Mean(A_iy)表示对象i内所有像素的纵坐标的均值；

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

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

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

(1)构建CRF一元能量项时，首先利用基于像素的K-means聚类获取各像素的类别。然后基于对象，统计对象内个像素的出现频次，进而构建鲁棒的一元能量项，将基于像素和面向对象技术二者结合，使得变化检测结果更稳健。

(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₂进行辐射归一化校正。该方法首先通过典型相关分析找到两期影像各波段亮度值的一个线性组合，得到变化信息增强的差异影像，通过对归一化的MAD变量平方和取阈值确定变化和未变化区域，然后通过未变化区域对应的两时相像元对的映射方程，完成相对辐射校正。

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

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

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

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

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

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

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

步骤5：计算多时相影像的差异影像，记各波段差异影像所组成的影像为X_M，继而利用K-Means算法进行聚类，根据聚类结果，确定CRF一元能量项。设定聚类数目K＝2，具体实施如下：

5.1)从随机选取K个聚类中心；

5.2)将每个数据点归类到离它最近的那个中心点所代表的类别中；

5.3)重新计算已经得到的各个类的中心。

5.4)迭代5.2)～5.3)步直至新的中心与原中心相等或小于指定阈值，算法结束。

5.5)K-Means算法收敛后，统计每一个对象内部属于不同类别的个数，记对象j属于非变化类的个数为N_j1，属于变化类的个数为N_j2，则对象j的一元能量项为：

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

式中，i,j表示图像相邻的两个节点或者对象(本发明中，将每一个对象看作图模型中的一个节点)；c_i,c_j表示两个节点i,j处的类别标签(两种类别：变化和非变化)；L(i，j)表示两个对象的公共边界的长度；[c_i≠c_j]表示0-1指示函数，F_M(i)、F_M(j)分别表示由对象i、j内部差异影像值的均值、坐标的均值构成的矢量，F_M(i)的构造如下：

F_M(i)＝[Mean(Dif_i)Mean(A_ix)Mean(A_iy)] (4)

式中，Mean(Dif_i)表示差异影像在对象i内的均值，Dif表示差异影像的幅值，Mean(A_ix)表示对象i内所有像素的横坐标的均值，Mean(A_iy)表示对象i内所有像素的纵坐标的均值。

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

式中U_j、U_i,j分别表示对象j或节点j(本发明中，将每一个对象看作图模型中的一个节点)的一元能量项、相邻对象j和i间的二元能量项，S表示节点的集合，E表示相邻节点的集合，参数β是正则化参数，用于控制一元能量项和二元能量项间的比重。

步骤10：根据优化算法循环信度传播(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所示。

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

方法	FP	FN	OE	k
					EM-MRF	13373	1367	14740	0.821
CRF	650	3913	4563	0.934
					OBCT	8813	780	9593	0.878
本发明方法	1624	2139	3763	0.947
					理想	0	0	0	1

由表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.一种面向对象的多光谱高分辨率遥感影像变化检测方法，其特征在于，包括以下步骤：

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

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

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

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

步骤5：计算多时相影像的差异影像，记差异影像为X_M，然后利用K-Means算法进行聚类，根据聚类结果，确定CRF一元能量项；

步骤6：根据SLIC的分割图像构建CRF的邻域***，并根据步骤5的差异影像和对象的空间坐标，构建CRF的二元项能量；

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

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

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

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

2.3)计算均方根误差；

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

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

3.根据权利要求2所述的面向对象的多光谱高分辨率遥感影像变化检测方法，其特征在于，求未知函数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所述的面向对象的多光谱高分辨率遥感影像变化检测方法，其特征在于，步骤2中，所述精校正即将经过几何粗校正的多光谱遥感影像数据，利用自动匹配与三角剖分法进行几何精校正。

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

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

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

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

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

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

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

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

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

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

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

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

9.根据权利要求1所述的面向对象的多光谱高分辨率遥感影像变化检测方法，其特征在于，步骤5中，基于预设的聚类数目K＝2，所述CRF一元能量项具体确定方法如下：

5.1)随机选取K个聚类中心；

5.2)将每个数据点归类到离它最近的那个中心点所代表的类别中；

5.3)重新计算已经得到的各个类的中心；

5.4)迭代步骤5.2)～步骤5.3)，直至新的中心与原中心相等或小于设定阈值；

5.5)K-Means算法收敛后，统计每一个对象内部属于不同类别的个数，记对象j属于非变化类的个数为N_j1，属于变化类的个数为N_j2，则对象j的一元能量项为：

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

式中，i,j表示图像相邻的两个节点即相邻的对象；c_i,c_j表示两个节点i,j处的类别标签；两种类别为变化和非变化；L(i，j)表示两个对象的公共边界的长度；[c_i≠c_j]表示0-1指示函数；F_M(i)、F_M(j)分别表示由对象i、j内部差异影像值的均值、坐标的均值构成的矢量，其中F_M(i)的构造如下：

F_M(i)＝[Mean(Dif_i) Mean(A_ix) Mean(A_iy)] (4)

式中，Mean(Dif_i)表示差异影像的幅值在对象i内的均值，Dif表示对象i内差异影像值的均值，Mean(A_ix)表示对象i内所有像素的横坐标的均值，Mean(A_iy)表示对象i内所有像素的纵坐标的均值；

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

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

Images