CN101937079A - 基于区域相似度的遥感影像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于区域相似度的遥感影像变化检测方法，包括：在监测区域内选择两景不同时相的遥感影像，所述遥感影像的影像数据包括多光谱数据和全色波段数据；对选择的遥感影像进行图像预处理；将预处理后的两期遥感影像的所有波段组成一景影像，对一景影像进行分割，建立两期所述遥感影像的对象的一一映射关系；基于分割结果，构建每个对象的特征集描述；基于所述对象的光谱特征，计算所述对象间的区域相似度，根据所述对象的纹理特征、NDVI特征进行代数计算；设定一个区域相似度阈值，结合代数计算的计算结果，建立变化区域检测准则，提取所述两期遥感影像的变化区域。本发明能有效消除“椒盐”现象，提高遥感影像变化区域提取精度。

Description

基于区域相似度的遥感影像变化检测方法

技术领域

本发明涉及遥感影像处理技术领域，特别涉及一种基于区域相似度的遥感影像变化检测方法。

背景技术

遥感信息变化检测是遥感应用的一个重要研究课题，用同一地区、不同时相的卫星遥感影像进行变化检测，是对土地等资源开展动态监测、资源变更调查、环境监测、基础地理数据库更新等对地观测应用中的关键技术，具有广泛的应用领域。从现有的变化检测方法来看，针对中低分辨率的遥感影像，主要还停留在像元级数据导引的方法，缺少知识导引的特征级变化检测方法。

近年来，技术的发展使得拍摄到遥感影像的分辨率越来越高，空间分辨率的提高有利于反映土地利用类型差异，有利于表征同一地物内部组成要素丰富的细节信息，但地物类别内部的光谱响应差异增大，导致地类的光谱特征不稳定。面向对象影像的分析方法是一种适用于高分辨率遥感影像的分析与信息提取方法，该方法结合影像的光谱特征、形状特征以及纹理特征等，将影像划分为若干个有意义的同质区域(对象)，将对象取代像元作为影像分析的基本单元，分析人员不仅仅可以利用光谱信息进行分析，还可以选择使用影像对象(多边形)有关的多种形状测度，对象纹理特征，以及对象空间相关特征等进行分析。一旦获得对象的这些特征，从理论上讲，基于对象的所有特征信息都可以有效地用于遥感信息提取过程，这就使影像分析更加灵活和稳定。目前该方法已成功应用于土地利用和土地覆盖分类等研究中。

将面向对象影像的分析方法应用于变化检测研究，是进行高分辨率影像变化检测的一种行之有效的途径。现有技术常用的是分类后比较变化检测方法，该方法对两个时相的遥感影像分别进行分割和分类，然后对两幅分类后的图逐像元进行比较。然而，单时相影像中的任何分类错误在最终的变化检测图中都会表现出来，称之为“椒盐”现象，因此该方法的主要缺点是变化检测的精度依赖于两个独立的分类结果的精度，存在分类误差累积现象。

发明内容

(一)主要解决的技术问题

本发明要解决的技术问题是提高遥感影像变化检测的精度，消除检测结果“椒盐”现象。

(二)技术方案

为此，本发明提供了一种基于区域相似度的遥感影像变化检测方法，包括：

步骤10、在监测区域内选择两景不同时相的遥感影像，所述遥感影像的影像数据包括多光谱数据和全色波段数据；

步骤20、对步骤10选择的遥感影像进行图像预处理；

步骤30、将预处理后的两期遥感影像的所有波段组成一景影像，对所述一景影像进行分割，建立两景所述遥感影像的对象的一一映射关系；

步骤40、基于分割结果，构建每个对象的特征集描述，所述特征集描述包括光谱特征、纹理特征以及经过变换的二次特征；

步骤50、基于所述对象的光谱特征，计算所述两期遥感影像对象间的区域相似度，根据所述对象的纹理特征、归一化植被指数(NDVI)特征进行代数计算，计算结果作为区域变化检测的参考值；

步骤60、设定一个区域相似度阈值，结合所述代数计算的计算结果建立变化区域检测准则，提取所述遥感影像的变化区域。

其中，所述步骤20中的所述预处理包括辐射校正、数据融合、精确几何配准以及影像降噪处理。

所述预处理具体为：采用回归分析归一化法对所述遥感影像进行相对辐射校正；采用Gram-Schmidt融合方法对辐射校正后的遥感影像进行数据融合；采用多项式或有理函数法对数据融合后的遥感影像进行配准；采用基于形态学的混合开闭重建滤波运算方法对配准后的遥感影像进行降噪处理。

所述步骤30中采用基于标记控制的分水岭分割方法对所述一景影像进行多波段、多尺度分割。

所述步骤50中根据下式计算所述区域相似度，

$\mathrm{RSIM}=\frac{4{\mathrm{\σ}}_{12}\stackrel{\‾}{{\mathrm{BV}}_1}\stackrel{\‾}{{\mathrm{BV}}_2}}{({\mathrm{\σ}}_1^2+{\mathrm{\σ}}_2^2)[{\left(\stackrel{\‾}{{\mathrm{BV}}_1}\right)}^2+{\left(\stackrel{\‾}{{\mathrm{BV}}_2}\right)}^2]}$

其中：

$\stackrel{\‾}{{\mathrm{BV}}_1}=\frac{1}{B}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{\mathrm{BV}}_{i1}$

$\stackrel{\‾}{{\mathrm{BV}}_2}=\frac{1}{B}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{\mathrm{BV}}_{i2}$

${\mathrm{\σ}}_1^2=\frac{1}{B-1}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{({\mathrm{BV}}_{i1}-\stackrel{\‾}{{\mathrm{BV}}_1})}^2$

${\mathrm{\σ}}_2^2=\frac{1}{B-1}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{({\mathrm{BV}}_{i2}-\stackrel{\‾}{{\mathrm{BV}}_2})}^2$

${\mathrm{\σ}}_{12}=\frac{1}{B-1}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}({\mathrm{BV}}_{i1}-\stackrel{\‾}{{\mathrm{BV}}_1}){({\mathrm{BV}}_{i2}-\stackrel{\‾}{{\mathrm{BV}}_2})}^2$

RSIM为区域相似度，B为经过预处理后的遥感影像的波段数，BV_i1表示某对象区域内，位于时相1波段i中所有像素值的平均值；BV_i2表示某对象区域内，位于时相2波段i中所有像素值的平均值；

表示某对象区域内，位于时相i中所有波段所有像素值的平均值。

所述步骤50中：

所述代数计算包括差值计算或比值计算，计算结果包括：差值纹理Entropy特征图或比值纹理Entropy特征图，以及差值NDVI特征图或比值NDVI特征图。

所述步骤40中采用面向对象的灰度共生矩阵计算对象的纹理特征，所述灰度共生矩阵的窗口是所述对象本身，所述窗口的大小是所述对象的大小；

所述纹理特征包括：纹理均值和纹理熵，

$\mathrm{Mean}=\underset{i=0}{\overset{{\mathrm{quant}}_k}{\mathrm{\Σ}}}\underset{j=0}{\overset{{\mathrm{quant}}_k}{\mathrm{\Σ}}}p(i,j)*i$

$\mathrm{Entropy}=-\underset{i=0}{\overset{{\mathrm{quant}}_k}{\mathrm{\Σ}}}\underset{j=0}{\overset{{\mathrm{quant}}_k}{\mathrm{\Σ}}}p(i,j)*\ln p(i,j)$

其中，Mean表示纹理均值，Entropy表示纹理熵，quant_k表示原图像灰度压缩的等级，即灰度共生矩阵的维数，p(i，j)表示共生矩阵中元素的值，也即概率值。

所述经过变换的二次特征包括归一化植被指数(NDVI)，根据下式计算所述归一化植被指数，

$\mathrm{NDVI}=\frac{{\mathrm{\ρ}}_{\mathrm{nir}}-{\mathrm{\ρ}}_{\mathrm{red}}}{{\mathrm{\ρ}}_{\mathrm{nir}}+{\mathrm{\ρ}}_{\mathrm{red}}}$

其中：NDVI表示归一化植被指数，ρ_nir表示近红外波段的反射率，ρ_red表示红波段的反射率。

所述步骤60包括：

根据下式计算所述区域相似度、纹理特征和NDVI特征作为变化区域的参考阈值，

T＝Mean±nR

其中：T为参考阈值；Mean为区域相似度(或差值图像)的均值；n为常数；R为区域相似度(或差值图像)的标准差，n取0.5；

通过实验调整所述参考阈值得到变化阈值，将各个变化阈值结合形成检测准则，根据所述检测准则提取所述遥感影像的变化区域。

(三)有益效果

本发明提供的技术方案具有如下有益效果：对遥感影像进行分割得到对象，通过区域相似度来评价两个遥感影像的相似程度，结合对象纹理特征、NDVI等特征建立变化区域检测准则，最终提取到两个时相的遥感影像的变化区域。本发明提高了对遥感影像变化区域的检测精度，消除了现有技术检测结果中的图斑“椒盐”现象。

附图说明

图1是本发明实施例的基于区域相似度的遥感影像变化检测方法流程图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

如图1所示，是本发明实施例的基于区域相似度的遥感影像变化检测方法流程图，本实施例包括如下步骤：

步骤10：在监测区域内选择两景不同时相的遥感影像，该遥感影像的影像数据包括多光谱数据和全色波段数据；

对遥感影像的变化检测需要考虑两方面重要因素：遥感影像的自身特征和环境特征。理想情况下，用于变化检测的遥感影像的分辨率应该与空间、光谱和辐射的分辨率一致；而且在变化检测过程中，应该尽量保持大气条件、土壤湿度状况、物候周期等环境参数相近，否则不一致的环境因素将降低变化检测的精度。本实施例选择同一区域、同一传感器所采集的，季节差异不大的两个时相的高分辨率遥感影像，进行高分辨率遥感影像变化检测。

步骤20：对选择的遥感影像进行图像预处理；

本实施例的预处理包括辐射校正、数据融合、精确几何配准以及影像降噪处理；包括如下步骤：

步骤21、辐射校正预处理：绝对的辐射校正所需参数难以获取，本实施例采用回归分析归一化法进行两时相遥感影像的相对辐射校正；以一个时相的遥感影像作为参考影像，选择辐射稳定、光谱稳定的地物作为辐射控制点，即伪不变特征点(Pseudo Invariant Features，PIF)；基于伪不变特征点，采用回归分析求解两时相遥感影像之间的光谱特性关系参数；根据该光谱特性关系参数，得到与参考影像具有相同或相近辐射值的结果影像，完成对另一时相的遥感影像的相对辐射校正；

步骤22、数据融合预处理：分别对经过辐射校正后的两期遥感影像的全色影像和多光谱影像进行融合处理，得到分辨率与全色波段一致、波段数与多光谱波段相同的融合后的遥感影像；

本实施例采用的是Gram-Schmidt融合方法，该方法能够在保证原始影像高分辨率特性的基础上，很好的保持多波段影像的光谱信息；

步骤23、运用多项式或有理函数法对步骤22得到的遥感影像进行精确配准，具体地，通过将两时相的遥感影像中的同名点进行匹配，最终使得同名地物出现在配准后的影像的相同位置，以检测两期遥感影像之间的变化；将配准后的遥感影像作为最终进行变化检测的遥感数据；

步骤24、采用基于形态学的混合开闭重建滤波运算(Morphological Reconstruction)方法，分别对步骤23得到的遥感影像的各波段数据进行降噪处理，以减弱分水岭分割的“过分割”现象。

步骤30：将步骤20得到的两期遥感影像的所有波段进行整合，组成一景影像，选择适宜的尺度对该一景影像进行影像分割；

本实施例采用基于标记控制(Marker-controlled)的分水岭分割算法对步骤20得到的遥感影像进行初始分割，然后基于区域邻接图(Region Adjacency Graph，RAG)将初始分割得到的小区域合并为一景影像，合并判断准则需要综合考虑小区域的光谱特征和形状特征。获得两期遥感影像一一对应的分割对象是该变化检测方法的关键步骤，本发明提供两种方案建立两期遥感影像分割对象的一一映射关系。

方案1包括如下步骤：

步骤31’、选择一个时相的遥感影像作为变化检测的基准影像；

步骤32’、选择适宜的尺度对步骤31’选择的影像进行分割，尺度选择采用多次实验结合目视的方式来确定，分割结束后获得所有对象的边界信息；

步骤33’、将得到的边界信息转换为面状矢量文件，作为另外一个时相的影像分割的专题图层；

步骤34’、在步骤33’获得的矢量专题图层的约束下，对另一个时相的影像进行分割，最终获得的分割结果在数目和位置都与对上一时相的影像的分割结果一致，即建立了两期影像对象的一一映射关系。

该方案在进行第一次分割时，根据分水岭分割算法，分割结果中的对象认为是同质的，当以该分割对象获得的专题图层为约束条件对另一时相的影像进行划分时，将获得部分不同质的区域，这样的区域是由于发生了变化而造成的。在实际操作中，适当减小第一次分割尺度，能够保证第二次分割产生对象的同质性。

方案2：

本方案中，两个对象之间满足的理想情况是：

1、分别属于两个时相的遥感影像的同一区域，满足对应关系；

2、两时相的遥感影像分割后对应的对象都是同质区域。

基于以上两点，本发明实施例采用该技术方案建立两个时相的对象的一一映射关系，步骤如下：

步骤31、依次将步骤20预处理后的遥感影像的所有波段组合在一起，产生一景影像；

步骤32、基于标记控制的分水岭分割算法对一景影像进行多波段、多尺度影像分割。

这种分割方式考虑到了两时相的所有波段数据，因此产生的对象是近似的同质区域，而且所有对象对于两时相的遥感影像数据来说都是一一对应的。实验证明当分割尺度选择较小的时候，有利于提高分割对象的同质性，同时也能保证最终的变化检测精度。

下面本实施例将基于两时相遥感影像的分割结果进行区域相似度计算。

在现有技术中，Zhou Wang和C.Bovik等人于2002年首次提出了结构信息的概念，结构相似度(Structural Similarity Index，SSIM)适用于很多图像处理应用，提供了一种评价图像失真的定量表达参数。本发明实施例基于SSIM和影像分割区域，提出区域相似度(Region Similarity，RSIM)的概念，并以RSIM作为评价两个区域间相似程度的参数，形成一种新的面向对象高分辨率遥感影像的变化检测方法，用来检测和提取两时相遥感影像间的变化区域。

假设时相1某一个区域表示为：X＝{x_i|i＝1，2，...，N}；时相2中与之对应的区域表示为：Y＝{y_i|i＝1，2，...，N}，则RSIM可以表示为：

$\mathrm{RSIM}=\frac{4{\mathrm{\σ}}_{\mathrm{xy}}\stackrel{\‾}{x}\stackrel{\‾}{y}}{({\mathrm{\σ}}_x^2+{\mathrm{\σ}}_y^2)[{\left(\stackrel{\‾}{x}\right)}^2+{\left(\stackrel{\‾}{y}\right)}^2]}---\left(1\right)$

其中

$\stackrel{\‾}{x}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_i,\stackrel{\‾}{y}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{y}_i$

${\mathrm{\σ}}_x^2=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^2,{\mathrm{\σ}}_y^2=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(y_i-\stackrel{\‾}{y})}^2$

${\mathrm{\σ}}_{\mathrm{xy}}=\frac{1}{N-1}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(x_i-\stackrel{\‾}{x})(y_i-\stackrel{\‾}{y})$

RSIM的变化范围为[-1，1]，仅当x_i＝y_i，i＝1，2，...，N时，RSIM＝1；当，i＝1，2，...，N时，RSIM＝-1。区域相似度参数由对象间的相关性、亮度和对比度三部分组合而成。为了更好的理解，将RSIM改写为：

$\mathrm{RSIM}=\frac{{\mathrm{\σ}}_{\mathrm{xy}}}{{\mathrm{\σ}}_x{\mathrm{\σ}}_x}\×\frac{2\stackrel{\‾}{x}\stackrel{\‾}{y}}{{\left(\stackrel{\‾}{x}\right)}^2+{\left(\stackrel{\‾}{y}\right)}^2}\×\frac{2{\mathrm{\σ}}_x{\mathrm{\σ}}_x}{{\mathrm{\σ}}_x^2+{\mathrm{\σ}}_y^2}---\left(2\right)$

公式(2)的第一部分

表示时相1和时相2的遥感影像数据间的相关系数，为两数据间的线性相关系数测度，其变化范围为[-1，1]，当y_i＝ax_i+b，i＝1，2，...，N且a＞0时，该线性相关系数测度为1，尽管两时相数据线性相关，但是仍可能存在相对失真，失真的信息将用公式(2)的第二和第三部分进行评价，第二部分

变化范围为[0，1]，评价了两时相的数据间的亮度相似性，为亮度相似性参数，仅当

时，亮度相似性参数值达到1；标准差通常用来评价图像的对比度，所以公式(2)的第三部分

表示两时相的遥感影像间对比度，为对比度参数，变化范围为[0，1]，仅当σ_x＝σ_x时，对比度参数达到1。

对于本实施例，经过预处理后的两时相的遥感影像，波段数都为B，将B个波段重新组合成一景影像进行分割，得到分割区域，新的影像具有2B个波段。针对分割后的某一个对象，分别对于两期影像的波段数据进行统计，变量说明如下：

$\stackrel{\‾}{{\mathrm{BV}}_1}=\frac{1}{B}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{\mathrm{BV}}_{i1}---\left(3\right)$

$\stackrel{\‾}{{\mathrm{BV}}_2}=\frac{1}{B}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{\mathrm{BV}}_{i2}---\left(4\right)$

BV_i1表示某对象区域内，位于时相1波段i中所有像素值的平均值；BV_i2表示某对象区域内，位于时相2波段i中所有像素值的平均值；

表示某对象区域内，位于时相i中所有波段所有像素值的平均值，

则区域相似度(RSIM)可以表示为：

$\mathrm{RSIM}=\frac{4{\mathrm{\σ}}_{12}\stackrel{\‾}{{\mathrm{BV}}_1}\stackrel{\‾}{{\mathrm{BV}}_2}}{({\mathrm{\σ}}_1^2+{\mathrm{\σ}}_2^2)[{\left(\stackrel{\‾}{{\mathrm{BV}}_1}\right)}^2+{\left(\stackrel{\‾}{{\mathrm{BV}}_2}\right)}^2]}---\left(5\right)$

其中：

${\mathrm{\σ}}_1^2=\frac{1}{B-1}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{({\mathrm{BV}}_{i1}-\stackrel{\‾}{{\mathrm{BV}}_1})}^2$

${\mathrm{\σ}}_2^2=\frac{1}{B-1}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{({\mathrm{BV}}_{i2}-\stackrel{\‾}{{\mathrm{BV}}_2})}^2$

${\mathrm{\σ}}_{12}=\frac{1}{B-1}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}({\mathrm{BV}}_{i1}-\stackrel{\‾}{{\mathrm{BV}}_1}){({\mathrm{BV}}_{i2}-\stackrel{\‾}{{\mathrm{BV}}_2})}^2$

本实施例中区域相似度参数比传统的相关系数评价系数方法更加稳定，更能准确表达两个对象之间的相似性，而不再简单的表达两个对象间的线性相关性。

经过方案2建立起了时相1和时相2对象间的一一映射对应关系，针对一一映射对应的对象计算区域相似度，最终得到区域相似度影像图，作为建立检测变化区域准则的最主要特征。

步骤40：基于分割结果，对步骤30得到的对象，构建每个对象的特征集描述，对象的特征集描述主要包括光谱特征、纹理特征以及经过变换的二次特征；其中经过变换的二次特征如归一化植被指数(NDVI)等；

其中对象的光谱特征是指该对象在某一波段上所有像素的均值。对象的纹理特征的计算采用面向对象的灰度共生矩阵(GLCM)计算，在计算面向对象的GLCM时，基窗口是对象本身，窗口的大小即对象的大小。

传统的GLCM计算过程中，一般都是先把图像的灰度级别进行量化压缩，比如压缩到16、32或64灰度级，然后再进行GLCM的计算。然而，对整个图像进行灰度级别压缩之后，一些原本具有纹理特征的区域将会减弱甚至丢失纹理信息。本实施例首先统计对象内所有灰度级别，然后基于这些灰度级别再进行GLCM的计算。通过上述方式来计算，不仅减少了GLCM的运算量，提高了GLCM纹理特征的计算速度，而且没有损失对象的纹理信息，增强了具有不同纹理特征的对象的可分性。

本实施例根据公式(6)和公式(7)计算纹理均值(Mean)和纹理熵(Entropy)两个特征：

$\mathrm{Mean}=\underset{i=0}{\overset{{\mathrm{quant}}_k}{\mathrm{\Σ}}}\underset{j=0}{\overset{{\mathrm{quant}}_k}{\mathrm{\Σ}}}p(i,j)*i---\left(6\right)$

$\mathrm{Entropy}=-\underset{i=0}{\overset{{\mathrm{quant}}_k}{\mathrm{\Σ}}}\underset{j=0}{\overset{{\mathrm{quant}}_k}{\mathrm{\Σ}}}p(i,j)*\ln p(i,j)---\left(7\right)$

其中quant_k表示变化检测影像灰度压缩的等级，即共生矩阵的维数，p(i，j)表示共生矩阵中元素的值，也就是概率值。

对象NDVI特征根据公式(8)进行计算：

$\mathrm{NDVI}=\frac{{\mathrm{\ρ}}_{\mathrm{nir}}-{\mathrm{\ρ}}_{\mathrm{red}}}{{\mathrm{\ρ}}_{\mathrm{nir}}+{\mathrm{\ρ}}_{\mathrm{red}}}---\left(8\right)$

其中ρ_nir表示近红外波段的反射率，ρ_red表示红波段的反射率。

经过公式(8)的比值处理，可以减少两时相遥感影像多个波段中很多形式的乘性噪声，有利于提高变化检测区域的提取精度。

步骤50：基于步骤30，获得了区域相似度影像图，在该图中，区域相似度值越大，说明两时相的对象相似度越高，即代表该对象实际发生变化的概率越小，区域相似度值越小代表对象实际发生变化的概率越高。同时利用区域对象的纹理、NDVI等特征值进行代数计算，计算结果作为区域变化检测的参考值；

步骤40中对两时相的对象分别计算光谱特征、纹理特征、NDVI形成特征影像图，即特征集描述，然后步骤50基于特征影像图，对应对象间的特征进行代数运算，如差值运算或比值运算，最终形成差值纹理Mean特征图或比值纹理Mean特征图、差值纹理Entropy特征图或比值纹理Entropy特征图，以及差值NDVI特征图或比值NDVI特征图，这两类特征图将作为建立检测变化区域准则的辅助特征，即参考值。

步骤60：设定某一个区域相似度阈值，结合纹理特征、NDVI等特征建立变化区域检测准则，提取两期遥感影像的变化区域；

由于面向对象的变化检测方法中评价的单元是对象，如果仅以光谱特征计算区域相似度，容易将光谱统计特征相似、纹理特征不同的区域误认为未发生变化的区域；另一方面，如果两时相遥感影像存在季节差异性，植被在不同生长期的光谱特征也就存在较大差异，单一基于区域光谱特征相似度的变化检测容易将上述两类对象误认为发生变化的区域。基于以上考虑，本发明实施例以区域相似度为主要判别特征，辅以纹理Mean、纹理Entropy、NDVI等特征差值，在综合考虑上述特征基础上，建立进行变化区域提取的判别准则。

针对单个特征，如何选择合适的阈值，本实施例进行变化检测的关键。虽然很多研究人员提出了一些阈值确定的方法，但只能作为参考，目前阈值的边界仍无法准确地预知，只能通过经验，或者结合目视实验确定。

本发明提出的基于区域相似度的面向对象的遥感影像变化检测方法，首先采用经验阈值法获得各个特征的参考阈值，然后通过实验、调整的策略确定各个特征的变化阈值，最后将各个特征的变化阈值结合起来形成最终的变化区域判别准则。

经验阈值法是一种人工交互式的阈值确定方法，最常用的是n倍标准差法。这种方法首先假设研究区只有少量区域的覆盖类型发生了变化，并且差异图像呈正态分布，在此假设下，阈值可以根据下式确定：

T＝Mean±nR                                  (9)

其中：T为阈值；Mean为区域相似度(或差值图像)的均值；n为常数；R为区域相似度(或差值图像)的标准差。实验证明，随着研究区的差异以及研究对象和目的的不同，n的取值皆有不同。本实施例中n取0.5。

根据上述建立的变化区域提取准则，逐一对每个对象进行检测，最终提取得到变化区域，将变化区域矢量化，更新至地理信息数据库存储。

现有技术中的变化检测通常是基于像素，或者是某个大小的窗口。本发明提出的基于区域相似度(RSIM)的高分辨率遥感影像变化检测方法，通过计算两个时相对应对象(图斑)间的相似度，并辅以纹理、NDVI等特征来确定变化的区域，虽然变化区域判别准则依赖于研究者的经验知识，但该方法提出的RSIM能够很好的表达对象之间的相似性；对图像噪声有很好的鲁棒性，使算法更加稳定；而且提取的变化区域避免了板块破碎现象，易于矢量化，在土地利用遥感动态监测和全国“一张图”工程有着广阔的应用潜力。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变型，这些改进和变型也应视为本发明的保护范围。

Claims (9)

1.一种基于区域相似度的遥感影像变化检测方法，其特征在于，包括：

步骤10、在监测区域内选择两景不同时相的遥感影像，所述遥感影像的影像数据包括多光谱数据和全色波段数据；

步骤20、对步骤10选择的遥感影像进行图像预处理；

步骤30、将预处理后的两期遥感影像的所有波段组成一景影像，对所述一景影像进行分割，建立两期所述遥感影像的对象的一一映射关系；

步骤40、基于分割结果，构建每个对象的特征集描述，所述特征集描述包括光谱特征、纹理特征以及经过变换的二次特征；

步骤50、基于所述对象的光谱特征，计算所述两期遥感影像对象间的区域相似度，根据所述对象的纹理特征、归一化植被指数(NDVI)特征进行代数计算，计算结果作为区域变化检测的参考值；

步骤60、设定一个区域相似度阈值，结合所述代数计算的计算结果建立变化区域检测准则，提取所述遥感影像的变化区域。

2.如权利要求1所述的基于区域相似度的遥感影像变化检测方法，其特征在于，所述步骤20中的所述预处理包括辐射校正、数据融合、精确几何配准以及影像降噪处理。

3.如权利要求2所述的基于区域相似度的遥感影像变化检测方法，其特征在于，所述预处理具体为：采用回归分析归一化法对所述遥感影像进行相对辐射校正；采用Gram-Schmidt融合方法对辐射校正后的遥感影像进行数据融合；采用多项式或有理函数法对数据融合后的遥感影像进行配准；采用基于形态学的混合开闭重建滤波运算方法对配准后的遥感影像进行降噪处理。

4.如权利要求1所述的基于区域相似度的遥感影像变化检测方法，其特征在于，所述步骤30中采用基于标记控制的分水岭分割方法对所述一景影像进行多波段、多尺度分割。

5.如权利要求1所述的基于区域相似度的遥感影像变化检测方法，其特征在于，所述步骤50中根据下式计算所述区域相似度，

$\mathrm{RSIM}=\frac{4{\mathrm{\σ}}_{12}\stackrel{\‾}{{\mathrm{BV}}_1}\stackrel{\‾}{{\mathrm{BV}}_2}}{({\mathrm{\σ}}_1^2+{\mathrm{\σ}}_2^2)[{\left(\stackrel{\‾}{{\mathrm{BV}}_1}\right)}^2+{\left(\stackrel{\‾}{{\mathrm{BV}}_2}\right)}^2]}$

其中：

$\stackrel{\‾}{{\mathrm{BV}}_1}=\frac{1}{B}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{\mathrm{BV}}_{i1}$

$\stackrel{\‾}{{\mathrm{BV}}_2}=\frac{1}{B}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{\mathrm{BV}}_{i2}$

${\mathrm{\σ}}_1^2=\frac{1}{B-1}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{({\mathrm{BV}}_{i1}-\stackrel{\‾}{{\mathrm{BV}}_1})}^2$

${\mathrm{\σ}}_2^2=\frac{1}{B-1}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}{({\mathrm{BV}}_{i2}-\stackrel{\‾}{{\mathrm{BV}}_2})}^2$

${\mathrm{\σ}}_{12}=\frac{1}{B-1}\underset{i=1}{\overset{B}{\mathrm{\Σ}}}({\mathrm{BV}}_{i1}-\stackrel{\‾}{{\mathrm{BV}}_1}){({\mathrm{BV}}_{i2}-\stackrel{\‾}{{\mathrm{BV}}_2})}^2$

RSIM为区域相似度，B为经过预处理后的遥感影像的波段数，BV_i1表示某对象区域内，位于时相1波段i中所有像素值的平均值；BV_i2表示某对象区域内，位于时相2波段i中所有像素值的平均值；

表示某对象区域内，位于时相i中所有波段所有像素值的平均值。

6.如权利要求1所述的基于区域相似度的遥感影像变化检测方法，其特征在于，所述步骤50中代数计算包括差值计算或比值计算，计算结果包括：差值纹理Entropy特征图或比值纹理Entropy特征图，以及差值NDVI特征图或比值NDVI特征图。

7.如权利要求1所述的基于区域相似度的遥感影像变化检测方法，其特征在于，所述步骤40中采用面向对象的灰度共生矩阵计算对象的纹理特征，所述灰度共生矩阵的窗口是所述对象本身，所述窗口的大小是所述对象的大小；

所述纹理特征包括：纹理均值和纹理熵，

$\mathrm{Mean}=\underset{i=0}{\overset{{\mathrm{quant}}_k}{\mathrm{\Σ}}}\underset{j=0}{\overset{{\mathrm{quant}}_k}{\mathrm{\Σ}}}p(i,j)*i$

$\mathrm{Entropy}=-\underset{i=0}{\overset{{\mathrm{quant}}_k}{\mathrm{\Σ}}}\underset{j=0}{\overset{{\mathrm{quant}}_k}{\mathrm{\Σ}}}p(i,j)*\ln p(i,j)$

其中，Mean表示纹理均值，Entropy表示纹理熵，quant_k表示原图像灰度压缩的等级，即灰度共生矩阵的维数，p(i，j)表示共生矩阵中元素的值，也即概率值。

8.如权利要求1所述的基于区域相似度的遥感影像变化检测方法，其特征在于，所述经过变换的二次特征包括归一化植被指数，根据下式计算所述归一化植被指数，

$\mathrm{NDVI}=\frac{{\mathrm{\ρ}}_{\mathrm{nir}}-{\mathrm{\ρ}}_{\mathrm{red}}}{{\mathrm{\ρ}}_{\mathrm{nir}}+{\mathrm{\ρ}}_{\mathrm{red}}}$

其中：NDVI表示归一化植被指数，ρ_nir表示近红外波段的反射率，ρ_red表示红波段的反射率。

9.如权利要求1所述的基于区域相似度的遥感影像变化检测方法，其特征在于，所述步骤60包括：

根据下式计算所述区域相似度、纹理特征和NDVI特征作为变化区域的参考阈值，

T＝Mean±nR

其中：T为参考阈值；Mean为区域相似度(或差值图像)的均值；n为常数；R为区域相似度(或差值图像)的标准差，n取0.5；

通过实验调整所述参考阈值得到变化阈值，将各个变化阈值结合形成检测准则，根据所述检测准则提取所述遥感影像的变化区域。

Images