一种基于历史数据挖掘的高分遥感影像土地覆盖变化检测方法

摘要

本发明涉及一种基于历史数据挖掘的高分遥感影像土地覆盖变化检测方法，包括：获取历史土地覆盖数据、对应历史期高分遥感影像和现时高分遥感影像，并在历史土地覆盖数据约束下，采用图像分割算法分割高分遥感影像，得到内部光谱、纹理匀质的对象；选取部分不变对象，对两期高分遥感影像进行相对辐射校正，并以辐射校正后的中误差为阈值，将其余对象分为不变对象和疑似变化对象；提取历史期影像中对象的光谱、纹理特征，并使用机器学习算法，获取分类规则；使用分类规则对现时高分遥感影像中的疑似变化对象分类，从而得到变化检测结果。本发明在使用高分遥感影像自动化更新土地覆盖数据、变化检测或者分类方面，具有重大的应用前景。

说明书

技术领域

本发明涉及高分遥感影像变化检测技术领域，特别涉及一种基于历史数据挖掘的高分遥 感影像土地覆盖变化检测方法。

背景技术

土地覆盖的空间分布反映着人类社会经济活动过程，决定着地表的水热和物质平衡，其 变化直接影响到生物地球化学循环，改变着陆地－大气的水分、能量和碳循环，以至引起气 候变化。科学准确地测定全球地表覆盖的空间分布与动态变化，对于研究地球系统的能量平 衡、碳循环及其他生物地球化学循环、气候变化等有着十分重要的意义(参见文献：陈军,陈 晋,宫鹏,等.全球地表覆盖高分辨率遥感制图[J].地理信息世界,2011,9(2):12-14.)。

遥感技术作为一种非接触式获取信息技术，长久以来，被视为大范围获取土地覆盖数据 最有效的技术(参见文献：Chen J,Chen J,Liao A,et al.Global land cover mapping at 30m  resolution:A POK-based operational approach[J].ISPRS Journal of Photogrammetry and Remote  Sensing,2014.)。传统土地覆盖变化检测方法主要包括代数运算法，变化向量分析法，空间 变换法，分类比较法(参见文献：(1)Singh A.Review article digital change detection techniques  using remotely-sensed data[J].International journal of remote sensing,1989,10(6):989-1003. (2)Lu,D.,et al.(2004)."Change detection techniques."International Journal of Remote Sensing 25(12):2365-2401.)。代数运算法主要指影像差值法，即利用两期遥感影像对应波段的差值与 选定阈值的大小进行比较，以确定变化区域；变化向量分析法与插值法相比，能够给出变化 类别的信息；空间变换法是利用两期影像波段信息的冗余性，将波段重新组合，对新的波段 进行分析。总体而言，以上方法都存在阈值难以确定的问题。分类比较法主要包括直接分类 法和分类后比较法两种。直接分类法是将两期影像联合，直接参与分类，获取变化检测结果； 分类后比较法是对两期影像分别分类，然后对分类结果进行比较得出变化检测结果。在土地 覆盖变化检测方面，分类比较法是一个比较活跃的算法，但由于时代等原因，传统分类比较 法所采用的处理手段与当前技术相比，自动化程度和检测精度都较低。

随着遥感技术的发展，可用于土地覆盖变化检测的卫星遥感影像出现许多新的特点，主 要包括卫星数据源增多，空间分辨率增加，时间分辨率提高，这些特点为遥感卫星影像变化 检测的精度、速度以及自动化方面提出了新的挑战。一方面，以面向对象为代表并融合机器 学习算法的高分遥感影像解译技术，为变化检测的智能化带来了新的机遇。一批以对象代替 像元或者融合智能分类算法的变化检测技术诞生，提高了变化检测的精度和自动化水平(参 见文献：Chen G,Hay G J,Carvalho L M T,et al.Object-based change detection[J].International  Journal of Remote Sensing,2012,33(14):4434-4457.)。但现有方法的对象生成策略主要有两 种，独立分割方法和复合分割方法，其中独立分割方法存在对象边界不匹配的问题，复合分 割方法则存在变化对象边界不准确的问题。另一方面，多数智能分类算法都存在样本库构建 的问题，这也限制了变化检测的自动化水平。随着遥感土地覆盖数据的积累，通过挖掘历史 数据提高当前的变化检测结果，成为遥感土地覆盖变化检测技术的一个新方向。一些研究人 员挖掘历史土地覆盖数据，以解决样本库构建的问题(参见文献：(1)Walter V.Object-based  classification of remote sensing data for change detection[J].ISPRS Journal of photogrammetry  and remote sensing,2004,58(3):225-238.；(2)Xian G,Homer C,Fry J.Updating the 2001 National Land Cover Database land cover classification to 2006by using Landsat imagery change  detection methods[J].Remote Sensing of Environment,2009,113(6):1133-1147.)。Walter(2004) 等仅利用了历史土地覆盖数据的边界来获取对象，且未挖掘历史期高分遥感影像中对象的信 息。Xian(2009)等挖掘了历史土地覆盖数据和历史期高分遥感影像中的信息，但仍是基于像 元且面向中等分辨率的遥感影像土地覆盖变化检测。如何针对新时期高分遥感影像的特点， 综合面向对象分类技术以及机器学习算法，挖掘历史土地覆盖数据和历史期高分遥感影像的 信息，对提高土地覆盖变化检测精度和自动化水平具有重要意义。

发明内容

针对上述问题，本发明提出了一种基于历史数据挖掘的高分遥感影像土地覆盖变化检测 方法，该方法以历史土地覆盖数据约束分割高分遥感影像，并利用机器学习算法挖掘土地覆 盖数据以及对应历史期高分遥感影像的分类规则信息，对现时高分遥感影像自动分类，从而 实现变化检测。

本发明技术解决方案：本发明提供了一种基于历史数据挖掘的高分遥感土地覆盖变化检 测方法，其特征在于，该方法包括：

步骤1：利用历史土地覆盖数据约束分割高分遥感影像，获取光谱、纹理匀质的对象， 该步骤主要包括：

步骤1-1：获取研究区域两期时相一致的高分遥感影像和历史期高分遥感影像对应的历 史土地覆盖数据；

步骤1-2：配准两期高分遥感影像；

步骤1-3：利用所述历史土地覆盖数据的边界范围作为约束条件，分割历史期高分遥感 影像；

步骤1-4：以历史期高分遥感影像分割结果为约束条件，分割现时高分遥感影像，从而 得到光谱、纹理匀质的对象。

步骤2：以相对辐射校正中误差为阈值，提取疑似变化对象，该步骤主要包括：

步骤2-1：选择若干未变化对象，以历史期影像为准，对现时高分影像进行相对辐射校 正，并计算校正后的中误差；

步骤2-2：计算两期高分影像对应波段的差值，并以对象为单元统计差值的均值特征；

步骤2-3：以相对中误差为阈值，将差值的均值特征超过阈值的对象设为疑似变化对象， 阈值以内的设为未变化对象。

步骤3：挖掘历史数据，获取土地覆盖分类规则信息，该步骤主要包括：

步骤3-1：计算历史期影像中对象的光谱、纹理特征，构成对象的特征向量；

步骤3-2：根据历史土地覆盖数据的类别信息，使用机器学习算法，获取土地覆盖分类 规则信息。

步骤4：对疑似变化对象分类，获得变化检测结果，该步骤主要包括：

步骤4-1：计算现时高分遥感影像中对象的光谱、纹理特征；

步骤4-2：使用步骤3获取的分类规则信息，对现时遥感影像中疑似变化对象分类；

步骤4-3：将疑似变化对象分类结果与历史土地覆盖数据中不一致的对象设为变化区域， 其余为未变化区域，从而得到最终的变化检测结果。

进一步地，步骤1中，所述研究区域两期时相一致的高分遥感影像指用于变化检测的高 分遥感影像分辨率基本相同，高分影像的成像时相相同；所述历史期高分遥感影像对应的历 史土地覆盖数据指用于变化检测的两期高分遥感影像中，历史时期的高分遥感影像成像时间 在历史土地覆盖数据所反应的土地覆盖信息有效时间范围内；所述历史土地覆盖数据的边界 范围为约束条件指分割历史期高分遥感影像时，对象的边界与历史土地覆盖数据的边界不存 在交叉现象；所述光谱、纹理匀质的对象指对象无论在历史期高分遥感影像中，还是在现时 高分遥感影像中，对象内部均保持光谱或者纹理的匀质性。

进一步地，其特征在于，所述步骤2-1中相对辐射校正指以历史期高分影像为准，对现 时高分遥感影像进行线性回归校正；相对中误差指相对辐射校正的中误差，具体公式如下：

DN₁₂(λ)_＝K·DN₂(λ)+B

$K=\frac{\underset{i=1}{\overset{n}{\Sigma }}(D{N}_2^i\left(\lambda \right)-\mathrm{mD}{N}_2\left(\lambda \right))(D{N}_1^i\left(\lambda \right)-\mathrm{mD}{N}_1\left(\lambda \right))}{\underset{i=1}{\overset{n}{\Sigma }}(D{N}_2^i\left(\lambda \right)-\mathrm{mD}{N}_2\left(\lambda \right))}$

B＝mDN₁(λ)-K·mDN₂(λ)

$\sigma \left(\lambda \right)=\sqrt{\frac{\underset{i=1}{\overset{n}{\Sigma }}{(D{N}_{12}^i\left(\lambda \right)-D{N}_1^i\left(\lambda \right))}^2}{n-2}}$

其中，DN₁(λ)、DN₂(λ)分别指历史期和现时高分遥感影像波段λ所对应的影像灰度值， DN₁₂(λ)指现时高分遥感影像波段λ相对辐射校正后的影像灰度值，指第 i个样本对象内的历史期和现时高分遥感影像波段λ所对应的影像灰度值的均值、 mDN₁(λ)、mDN₂(λ)指样本对象在历史期和现时高分遥感影像波段λ所对应的影像灰度值均 值特征的均值，K和B指相对辐射校正时线性回归的斜率和截距，σ(λ)指波段λ相对辐射 校正后的中误差。

步骤3是本文发明创新的关键部分，因此，下面对该步骤进一步展开说明。

该步骤的目的是挖掘历史土地覆盖数据以及与历史期影像相结合的信息，获取现时高分 遥感影像的分类规则。步骤3-1中的光谱特征指对象所在区域的波段均值、标准差，纹理特 征指基于灰度共生矩阵的角二阶矩、熵、逆差矩，公式如下：

$D{N}_1\left(\lambda \right)=\frac{1}{\#P}\underset{(x,y)\in P}{\Sigma }D{N}_1(\lambda ,x,y)$

${\sigma }_1\left(\lambda \right)=\sqrt{\frac{1}{\#P}\underset{(x,y)\in P}{\Sigma }{(D{N}_1(\lambda ,x,y)-D{N}_1\left(\lambda \right))}^2}$

$\mathrm{ASM}\left(\lambda \right)=\underset{i=1}{\overset{k}{\Sigma }}\underset{j=1}{\overset{k}{\Sigma }}{\left(G(i,j)\right)}^2$

$\mathrm{ENT}\left(\lambda \right)=-\underset{i=1}{\overset{k}{\Sigma }}\underset{j=1}{\overset{k}{\Sigma }}G(i,j)\log G(i,j)$

$\mathrm{IDM}\left(\lambda \right)=-\underset{i=1}{\overset{k}{\Sigma }}\underset{j=1}{\overset{k}{\Sigma }}\frac{G(i,j)}{1+{\left(G(i,j)\right)}^2}$

其中，DN₁(λ)指对象在历史期高分遥感影像波段λ的均值，σ₁(λ)指对象在历史期高分 遥感影像波段λ的标准差，指对象的像素集，#P指对象的像素个数，DN₁(λ,x,y)指历史 期高分遥感影像波段λ(x,y)的灰度值，G(i,j)是灰度共生矩阵的值，ASM(λ)指角二阶矩， ENT(λ)指熵，IDM(λ)指逆差矩。

灰度共生矩阵是对影像上保持某距离的两象素分别具有某灰度的状况进行统计得到的 矩阵。如果灰度共生矩阵中的值集中在某一块(比如对连续灰度值影像，值集中在对角线； 对结构化的影像，值集中在偏离对角线的位置)，则ASM有较大值，若G(i,j)中的值分布较 均匀(如噪声严重的影像)，则ASM有较小的值。若灰度共生矩阵值分布均匀，也即影像近 于随机或噪声很大，ENT会有较大值。如果影像纹理的不同区域间缺少变化，局部非常均匀， IDM就会取较大的值。

计算对象的光谱、纹理特征后，可将其联合，构成对象的特征向量，如下：

V＝(DN₁(λ₁),σ₁(λ₁),ASM(λ₁),ENT(λ₁),IDM(λ₁),

…,DN₁(λ_n),σ₁(λ_n),ASM(λ_n),ENT(λ_n),IDM(λ_n))^T

其中，λ₁,…,λ_n指波段1,2,…n。

步骤3-2中使用机器学习算法，获取高分遥感影像分类规则指利用决策树、支持向量机、 神经网络等智能分类算法，由历史土地覆盖数据提供的类别信息和步骤3-1中计算的特征， 训练分类器，得到分类规则(机器学习算法的原理与实现请参考模式识别、数据挖掘、机器 学习等相关领域的书籍或软件)。

进一步地，步骤4中所述的现时高分遥感影像中对象的光谱、纹理特征与步骤3-1中所 述特征一致，但是计算特征的数据是现时高分遥感影像。步骤4-2中的分类规则指步骤3中 学习得到的分类规则，且分类的对象是步骤2中提取的疑似变化对象。步骤4-3将疑似变化 对象分类结果与历史土地覆盖数据中不一致的对象设为变化区域，且其类别信息为步骤4-2 中通过分类得到的类别信息，而未变化对象的类别信息则保持为历史土地覆盖数据中的类别 信息，由此可以得到现时高分遥感影像的分类结果，最后对比历史土地覆盖数据和现时高分 遥感影像的分类结果，即可得到最终的变化检测结果。

本发明与现有技术相比的优点在于：

(1)本发明方法以历史土地覆盖数据约束分割高分遥感影像，并利用机器学习算法挖 掘历史土地覆盖数据以及与对应历史期高分遥感影像相结合的分类规则等信息，对现时高分 遥感影像自动分类，从而实现变化检测。

(2)本发明方法实现了高分遥感影像变化检测的高精度自动化，在使用高分遥感影像 自动更新土地覆盖数据、变化检测或者分类方面，可以大大节省成本，有利于土地利用的管 理、评估以及科学研究。

附图说明

图1是基于历史数据挖掘的高分遥感影像土地覆盖变化检测方法流程图；

图2是本方法实施例所使用的影像、数据及变化检测结果图，其中(a)为历史期高分 遥感影像，(b)为现时高分遥感影像，(c)为对应的历史土地覆盖数据，(d)为本实施例的 变化检测结果图。

具体实施方式

下面结合沿海某局部地区进行实验，对本发明的基于历史数据挖掘的高分遥感影像土地 覆盖变化检测方法进行详细描述。

在本发明实施例中，所采用的历史期高分遥感影像和现时高分遥感影像分别是2003年 10月23日和2010年11月9日的SPOT5多光谱影像(见图2中的(a)和(b))，包括四个波段： 绿波段(0.49-0.61um)、红波段(0.61-0.68um)、近红外(0.78-0.89um)和短波红外(1.58-1.78um)， 其中短波红外的分辨率为20米，其余皆为10米。此外，本实施例中采用的历史土地覆盖数 据为基于2003年SPOT5影像人工目视解译并综合地面调查得到的结果(见图2中的(c))。

本实施例分成以下四个主要部分：

(一)利用历史土地覆盖数据约束分割高分遥感影像，获取光谱、纹理匀质的对象。该 步骤主要包括：

步骤1-2：配准两期高分遥感影像。影像配准方法可参考ENVI、PCI、ERDAS等遥感 处理软件，配准精度在1个像素以内。

步骤1-3：在eCogniton软件中，以历史期高分遥感影像四个波段为影像图层，并且以历 史土地覆盖矢量数据为专题图层，采用多尺度分割算法(Multiresolution Segmentation)进行分 割，创建分割图层Date1Seg-level，分割参数设置为Scale＝10,Shape＝0.1,Compactness＝ 0.5。

步骤1-4：以现时高分遥感影像四个波段为影像图层，在分割图层Date1Seg-level中继 续以参数Scale＝10,Shape＝0.1,Compactness＝0.5分割，并创建分割图层Date2Seg-level。 分割图层Date2Seg-level即为光谱、纹理匀质的对象所在图层。

(二)以相对辐射校正中误差为阈值，提取疑似变化对象。该步骤主要包括：

步骤2-1：选择若干未变化对象，以历史期影像为准，对现时高分影像进行相对辐射校 正。本实施例选择了8个未变化对象，计算未变化对象光谱的均值，以均值作线性回归，计 算相对辐射校正系数和中误差，具体计算公式如下，回归系数和中误差结果见表1：

$K=\frac{\underset{i=1}{\overset{n}{\Sigma }}(D{N}_2^i\left(\lambda \right)-\mathrm{mD}{N}_2\left(\lambda \right))(D{N}_1^i\left(\lambda \right)-\mathrm{mD}{N}_1\left(\lambda \right))}{\underset{i=1}{\overset{n}{\Sigma }}(D{N}_2^i\left(\lambda \right)-\mathrm{mD}{N}_2\left(\lambda \right))}$

B＝mDN₁(λ)-K·mDN₂(λ)

$\sigma \left(\lambda \right)=\sqrt{\frac{\underset{i=1}{\overset{n}{\Sigma }}{(D{N}_{12}^i\left(\lambda \right)-D{N}_1^i\left(\lambda \right))}^2}{n-2}}$

表1 相对辐射校正系数与中误差

由以上计算得到的结果，对现时SPOT5影像分波段进行相对辐射校正，校正公式如下：

DN₁₂(λ)＝K·DN₂(λ)+B

步骤2-2：计算相对辐射校正后的SPOT5影像与历史期SPOT5影像分波段的差值，并 以对象为单元统计差值的均值。

步骤2-3：以相对辐射校正中误差σ(λ)为阈值，将对象差值绝对值的均值超过阈值的对 象设为疑似变化对象，阈值以内的设为未变化对象。

(三)使用智能分类算法，学习土地覆盖分类规则信息。本实施例选择的智能分类算法 为C5.0决策树(实现方法，可参见eCognition使用手册)。该步骤主要包括：

步骤3-1：由2003年SPOT5影像计算对象光谱、纹理特征，构成对象的特征向量：

V＝(DN₁(λ₁),σ₁(λ₁),ASM₁(λ₁),ENT₁(λ₁),IDM₁(λ₁),

…,DN₁(λ₄),σ₁(λ₄),ASM₁(λ₄),ENT₁(λ₄),IDM₁(λ₄))^T

步骤3-2：根据土地覆盖数据的类别信息，使用C5.0决策树算法获取土地覆盖分类规则 信息。

(四)对疑似变化对象分类，获得变化检测结果，该步骤主要包括：

步骤4-1：计算2010年SPOT5影像中疑似变化对象的光谱、纹理特征，构成特征向量：

V＝(DN₁₂(λ₁),σ₁₂(λ₁),ASM₁₂(λ₁),ENT₁₂(λ₁),IDM₁₂(λ₁),

…,DN₁₂(λ₄),σ₁₂(λ₄),ASM₁₂(λ₄),ENT₁₂(λ₄),IDM₁₂(λ₄))^T

步骤4-2：使用步骤3获取的分类规则信息，对2010年SPOT5影像中疑似变化对象分 类；

步骤4-3：将2010年SPOT5影像中疑似变化对象分类结果与2003年土地覆盖数据中类 别不一致的对象设为变化区域，并将分类结果作为2010年土地覆盖数据，而其余未变化区 域，保留其2003年的土地覆盖类别，作为2010年土地覆盖数据。由此，得到2010年SPOT5 的总体分类结果，同时，变化区域的分类结果为2003年到2010年土地覆盖的变化检测结果 (见图2中的(d))。

虽然上面结合本发明的优选实施例对本发明的原理进行了详细的描述，本领域技术人员 应该理解，上述实施例仅仅是对本发明的示意性实现方式的解释，并非对本发明包含范围的 限定。实施例中的细节并不构成对本发明范围的限制，在不背离本发明的精神和范围的情况 下，任何基于本发明技术方案的等效变换、简单替换等显而易见的改变，均落在本发明保护 范围之内。