CN104834942A - 基于掩膜分类的遥感影像变化检测方法及系统 - Google Patents

Abstract

本发明提供了一种基于掩膜分类的遥感影像变化检测方法及系统，包括：步骤1，多波段遥感影像的预处理；步骤2，对预处理后的多波段遥感影像进行单波段分离；步骤3，构造两时相单波段遥感影像的单波段差异影像；步骤4，构造变化掩膜差异影像，基于大津法分割变化掩膜差异影像；步骤5，将变化掩膜与两时相多波段遥感影像之一进行乘积融合，获得变化区域；步骤6，利用模糊C均值法对变化区域进行聚类；步骤7，根据变化区域的聚类结果，结合两时相多波段遥感影像之二的地物类型先验知识，确定变化区域的变化类型。本发明抗噪声干扰强，可有效缓解局部最优；可兼顾影像变化过程的定性和定量描述，变化检测精度高且可靠性强。

Description

基于掩膜分类的遥感影像变化检测方法及系统

技术领域

本发明属于摄影测量与遥感影像应用技术领域，具体涉及一种基于掩膜分类的遥感影像变化检测方法及系统，适用于实时、自动(或半自动)的检测遥感影像变化。

背景技术

遥感影像变化检测是通过分析同一地区不同时相的遥感图像，检测该地区的地物随时间变化的信息。随着航空航天技术的飞速发展，遥感观测数据以实时、快速、覆盖范围广、高时空分辨率等特点得到越来越广泛的应用。如何有效地提取海量数据中的变化信息，并用于环境、农业、生态系统以及人类面临的各种自然灾害的预防，已经成为目前遥感应用研究中的热点问题。近年来，国内外众多学者相继开展了许多变化检测的理论研究与系统开发工作，并取得一系列成果。

目前，常用的遥感影像变化检测方法有直接比较法和分类后比较法，前者操作简单、速度快，但只能定量描述目标区是否发生变化，很难确定变化性质；后者可提供变化类型信息，但必须进行两次影像分类且制定统一的分类标准，检测精度受到单独分类的误差传播影响。此外，中高分辨率遥感影像存在混合像元的现象，用传统“硬”分类法无法获得较高的精度，而模糊C均值聚类(FCM)是一种软聚类算法，利用隶属度使类与类之间没有明显界限，用来处理混合像元行之有效，但存在过分依赖初始聚类中心、类别数难以自动确定、对孤立点噪声数据敏感等缺陷。

文中涉及如下参考文献：

[1]吴芳,刘荣,田维春等.遥感变化检测技术及其应用综述[J].地理空间信息,2007,5(4):57-60.

[2]Lu D,M ausel P,Brondzio E,etal.Change detection techniques[J].Internationaljournal of remote sensing,2004,25(12):2365-2401.

[3]Robe M.Haralick,K.Shanmugam，Textural features for image classification[J].Systems,Man and Cybernetics,IEEE Transactions on,1973(6):610-621.

[4]Mas J F.Monitoring land-cover changes:a comparison of change detectiontechniques.International Journal of Remote Sensing，1999，20(1)：l39—152.

[5]李敏.基于模糊C均值算法的遥感图像变化检测的研究[D].湖南大学,2009.

发明内容

针对现有技术中存在的问题，本发明提供了一种可防止信息丢失、提高检测精度的基于掩膜分类的遥感影像变化检测方法及系统。

本发明将模糊聚类的思想引入到遥感影像变化检测中，在利用改进FCM进行分类处理基础上，采用一种掩膜分类变化检测法，来防止有用信息丢失、提高变化检测精度。

本发明主要包括两大阶段，第一阶段：构造变化掩膜；第二阶段：利用改进的模糊C均值聚类法进行掩膜分类，以检测影像变化。

第一阶段中，由于直接比较法中常用的比值法能增强变化信息、抑制背景信息，但有时会过于夸大部分变化；差值法直接且便于操作，但不能完全反应出地物辐射能量的变化。考虑到比值法和差值法各自的优势及局限性，本发明在差值和比值影像乘积变换融合的基础上，采用一种新的复合法构造差异影像，为后续的变化检测提供数据源。该复合法引入偏差因子，对灰度层级有调整作用，能增强差异影像的振幅及结构信息，更有利于变化检测。此外，由于遥感影像分辨率的不断提高，地物在多波段影像的不同通道存在波段间的敏感性差异，给变化检测带来虚检、漏检等影响。本发明综合考虑各个波段对变化检测的贡献度，采用邻域熵权法对单波段差异影像进行调权处理。信息熵度量某随机事件包含的信息量，该波段像元的邻域信息量越大或不确定性越高，该波段对应的权值越大、对变化检测的贡献度也越大。所以，在该阶段，分离R、G、B单波段后，采用差值/比值复合法分别构造R、G、B单波段的差异影像。再采用邻域熵权法对R、G、B单波段差异影像进行调权求和处理，权值大小对应各波段的贡献度，并结合Ostu阈值分割法确定二值变化掩膜。

第二阶段中，由于分类前比较法操作简单，检测执行速度快，能定量描述目标区是否发生了变化，但很难确定目标区域变化的性质。而分类后比较法可提供变化区域的变化类型信息，但精度受到误差传播的影响。因此，本发明综合考虑两者的优劣性，将二值变化掩膜与时相二多波段遥感影像乘积融合确定变化区域，并利用改进的模糊C均值聚类法(FCM)对变化区域聚类。最后参考时相一多波段遥感影像的地物类型先验知识，确定变化过程，即“从什么”变化“到什么”。掩膜分类法有效防止有用信息丢失、抑制噪声、提高变化检测精度。

由于模糊C均值聚类法是一种典型的“软聚类”，算法简单、收敛速度快，隶属度的引进使得处理混合像元行之有效。但该方法存在一些不足：如过分依赖初始聚类中心、类别数不能自动确定、忽略空间信息导致局部最优等。本发明对模糊C均值聚类法做两点改进，其一：基于区间划分思想，以各区间像元灰度均值作为初始聚类中心；其二：引入邻域像素的位置和属性信息得到相似度加权隶属度和加权空间函数，使模糊C均值聚类法具有更强的抗噪能力，有效避免局部最优，改善分类效果。

本发明的技术方案如下：

一、基于掩膜分类的遥感影像变化检测方法，包括步骤：

步骤1，多波段遥感影像的预处理，即多波段遥感影像的滤波和增强处理；

步骤2，对预处理后的多波段遥感影像进行单波段分离，获得单波段遥感影像；

步骤3，构造两时相单波段遥感影像的单波段差异影像；

步骤4，构造变化掩膜差异影像基于大津法分割变化掩膜差异影像，获得变化掩膜；

上述k表示单波段编号，k＝1,2,...N，N为波段数；Y_ij表示变化掩膜差异影像中(i,j)像元的灰度值；表示单波段k的差异影像中(i,j)像元的邻域熵权，H_ij-k表示单波段k的差异影像中(i,j)像元的邻域信息熵；C_ij-k表示单波段k的差异影像中(i,j)像元灰度值；

步骤5，将变化掩膜与两时相多波段遥感影像之一(记为：时相一多波段遥感影像)进行乘积融合，融合后灰度值与融合前灰度值保持一致的区域即变化区域；

步骤6，以变化区域像元灰度值为样本，基于区间划分思想对样本进行分区，以各区样本均值作为初始聚类中心；以邻域像元隶属度调权求和作为中心像元隶属度，各邻域像元对中心像元的邻域相似权利用模糊C均值法对变化区域进行聚类；

上述中，t、r分别表示邻域像元和中心像元在邻域窗口内的位置编号；x_r表示中心像元r的灰度值，x_t表示邻域像元t的灰度值；

步骤7，根据变化区域的聚类结果，结合两时相多波段遥感影像之二(记为：时相二多波段遥感影像)的地物类型先验知识，确定变化区域的变化类型。

步骤2具体为：

采用加权求和法将多波段遥感影像灰度化，当R、G、B单波段某一分量的权值为1，其余为0时，经加权求和即获得R、G、B单波段遥感影像。

步骤3中结合差值法和比值法构造单波段差异影像，具体为：

根据公式构造两时相单波段遥感影像的单波段差异影像，其中：

C_ij-k表示单波段k的差异影像中(i,j)像元灰度值；

R_ij-k表示单波段k的比值影像中(i,j)像元灰度值；

D_ij-k表示单波段k的差值影像中(i,j)像元灰度值；

a为比值影像所有像元在第k波段的灰度均值和差值影像中所有像元在第k波段的灰度均值之差；

b为比值影像所有像元在第k波段的灰度值均方根和差值影像所有像元在第k波段的灰度值均方根之差。

步骤4中，所述的基于大津法分割变化掩膜差异影像，具体为：

基于大津法单阈值分割法，将使变化掩膜差异影像的背景像元灰度和目标像元灰度的类间方差最大的阈值作为分割阈值；采用分割阈值分割变化掩膜差异影像，获得变化掩膜。

步骤6进一步包括子步骤：

6.1根据经验和目测观测结果设定类别数C，根据经验设定模糊度m和收敛值ε，随机设定初始隶属度矩阵，隶属度矩阵由变化区域中各像元对聚类中心的隶属度构成；令初始迭代步数l＝0；

6.2以变化区域像元灰度值为样本，基于区间划分思想将样本划分为C个子区间，以各子区间样本均值作为初始聚类中心；

6.3将邻域像元t当前隶属度的加权和作为中心像元r隶属度即邻域像元t的当前隶属度根据当前隶属度矩阵U获得，其权值 $w_{rt}=\frac{2x_r}{x_r+x_t};$

6.4采用公式 $\left\{\begin{array}{c}{u}_{rk}^{(l+1)}=\frac{1}{{\Σ}_{y=1}^C{\left(\frac{u_{rk}^{\left(l\right)}\×\left|\right|v_k^{\left(l\right)}-x_r\left|\right|}{u_{ry}^{\left(l\right)}\×\left|\right|v_y^{\left(l\right)}-x_r\left|\right|}\right)}^{\frac{2}{m-1}})},\∀y,\∀k\\ v_k^{(l+1)}=\frac{{\Σ}_{r=1}^n\[{\left(u_{rk}^{\left(l\right)}\right)}^m{x}_r\]}{{\Σ}_{r=1}^n{\left(u_{rk}^{\left(l\right)}\right)}^m}\end{array}\right.$ 修正隶属度矩阵及聚类中心，其中，t、r分别表示邻域像元和中心像元在邻域窗口内的位置编号，k、y为聚类中心编号；和分别表示第(l+1)和l次迭代中、中心像元r对第k类聚类中心的隶属度；表示第l次迭代中、中心像元r对第y类聚类中心的隶属度；x_r表示中心像元r的灰度值；和表示第l次迭代获得的第k类和第y类聚类中心，y在[1,C]范围内遍历；表示第(l+1)次迭代获得的第k类聚类中心；n为邻域窗口内像元数；

6.5比较隶属度矩阵范数||U^(l+1)-U^(l)||，U^(l+1)、U^(l)分别表示第(l+1)次迭代、第l次迭代获得的隶属度矩阵；若||U^(l+1)-U^(l)||≤ε，结束迭代，按照当前聚类中心和隶属度矩阵进行聚类；否则，令l＝l+1，执行步骤6.4。

步骤7进一步包括：

7.1对两时相多波段遥感影像之二进行监督分类，获得地物类型先验知识；

7.2结合两时相多波段遥感影像之一的聚类结果，采用人工方式确定变化区域的变化类型。

二、基于掩膜分类的遥感影像变化检测系统，包括：

预处理模块，用来多波段遥感影像的预处理，即多波段遥感影像的滤波和增强处理；

单波段分离模块，用来对预处理后的多波段遥感影像进行单波段分离，获得单波段遥感影像；

单波段分离模块，用来构造两时相单波段遥感影像的单波段差异影像；

变化掩膜构造模块，用来构造变化掩膜差异影像基于大津法分割变化掩膜差异影像，获得变化掩膜；

上述k表示单波段编号，k＝1,2,...N，N为波段数；Y_ij表示变化掩膜差异影像中(i,j)像元的灰度值；表示单波段k的差异影像中(i,j)像元的邻域熵权，H_ij-k表示单波段k的差异影像中(i,j)像元的邻域信息熵；C_ij-k表示单波段k的差异影像中(i,j)像元灰度值；

乘积融合模块，用来将变化掩膜与两时相多波段遥感影像之一进行乘积融合，融合后灰度值与融合前灰度值保持一致的区域即变化区域；

模糊C均值聚类模块，用来以变化区域像元灰度值为样本，基于区间划分思想对样本进行分区，以各区样本均值作为初始聚类中心；以邻域像元隶属度调权求和作为中心像元隶属度，各邻域像元对中心像元的邻域相似权利用模糊C均值法对变化区域进行聚类；

上述中，t、r分别表示邻域像元和中心像元在邻域窗口内的位置编号；x_r表示中心像元r的灰度值，x_t表示邻域像元t的灰度值；

变化检测模块，用来根据变化区域的聚类结果，结合两时相多波段遥感影像之二的地物类型先验知识，确定变化区域的变化类型。

传统的遥感影像变化检测法大多从单像素、单决策的角度出发，对空间结构复杂、细节丰富的中高分辨率遥感影像容易出现“同物异谱”、“异物通谱”所导致的波段敏感性等问题，且单决策变化检测法往往不能兼顾定性、定量描述变化区域和变化类型。

和传统的遥感影像变化检测方法相比，本发明具有如下优点和有益效果：

1、抗噪声干扰强，可有效缓解局部最优。

(1)在单波段差异影像构造环节，采用差值/比值复合法，并引入偏差因子对灰度层级进行调整，增强了遥感影像的变化信息、且能反应出地物辐射能量的变化。另外，考虑到地物在不同波段的敏感性差异给变化检测带来的影响，采用邻域熵权法对单波段差异影像进行调权求和处理，差异影像的振幅及结构信息也增强，更有利于变化检测。

(2)在采用改进的FCM机进行遥感影像变化检测中，基于区间划分思想，以各区间像元灰度均值作为初始聚类中心，并引入邻域空间信息得到相似度加权隶属度和加权空间函数，从而增强抗噪能力，有效避免局部最优，改善分类效果，提高变化检测精度。

2、可兼顾影像变化过程的定性和定量描述，变化检测精度高且可靠性强。

分类前比较法可定量描述目标区是否发生变化，但很难确定目标区域的变化性质。分类后比较法则可提供变化区域的变化类型信息，但精度受到误差传播影响。本发明利用变化掩膜定性描述可靠的变化区域，再单独对变化区域进行聚类，从而可在保证检测精度和可靠性的前提下，快速、有效地获得目标区域的变化类型信息。

附图说明

图1为本发明方法的具体流程示意图；

图2为改进的模糊C均值聚类流程图。

具体实施方式

本发明主要设计两部分内容，即：(1)单波段差异影像与二值变化掩膜的构造和(2)掩膜分类变化检测。本发明技术方案主要包括差值/比值复合法构造单分波段差异影像、Ostu阈值分割构造变化掩膜、基于模糊C均值聚类法和掩膜分类的影像变化检测。本发明可实现定性描述变化区域、定量描述变化类型，同时可缓解同物异谱、异物同谱引起的单波段敏感性差异对变化检测虚检、漏检的影响，能有效抑制孤立噪声干扰、避免局部最优等，还具有较优的时效性和精度。

下面将结合附图和具体实施方式进一步说明本发明，具体步骤如下：

步骤1：多波段遥感影像的预处理，即多波段遥感影像的滤波和增强处理。

采用高斯滤波器和直方图均衡化增强法对多波段遥感影像进行滤波及增强处理，以增大影像对比度、抑制噪声干扰、改善视觉效果。

步骤2：预处理后多波段遥感影像的单波段分离，获得单波段遥感影像。

为了能适应单波段对比熵权修正多波段联合互信息量的波段修正变化检测处理，具体实施中，采用加权求和法将多波段遥感影像灰度化，当R、G、B单波段某一分量的权值为1，其余为0时，加权求和即实现R、G、B单波段分离，获得R、G、B单波段遥感影像。

步骤3：采用差值/比值复合法构造单波段差异影像。

为增强遥感影像的变化信息、抑制背景信息，同时反应出地物辐射能量变化，本具体实施中采用差值/比值复合法分别构造R、G、B三个单波段的差异影像。

单波段差异影像的常规构造方法中，比值法能增强影像变化信息、抑制背景信息，但有时会过于夸大部分变化；差值法直接且便于操作，但不能完全反应地物辐射能量的变化。考虑到两者的优势及局限性，本发明在差值和比值影像乘积融合的基础上，引入偏差因子，对灰度层级有调整作用，能增强差异影像的振幅及结构信息，更有利于变化检测。

差值/比值复合法的原理如下：

$C_{ij-k}=a+b\×\sqrt{D_{ij-k}*R_{ij-k}}---\left(1\right)$

$R_{ij-k}=255\×\frac{\min \[x_{ij-k}\left(t_1\right),x_{ij-k}\left(t_2\right)\]}{\max \[x_{ij-k}\left(t_1\right),x_{ij-k}\left(t_2\right)\]}---\left(2\right)$

D_ij-k＝|x_ij-k(t₂)-x_ij-k(t₁)|                  (3)

$a=|u_{D_{ij-k}}-u_{R_{ij-k}}|---\left(4\right)$

$b=\frac{{\mathrm{\σ}}_{D_{ij-k}}}{{\mathrm{\σ}}_{R_{ij-k}}}---\left(5\right)$

公式(1)为差值/比值复合法原理公式。公式(2)为比值法原理公式，公式(3)为差值法原理公式，比值法和差值法均基于多时相单波段遥感影像的像元点对进行处理。

公式(1)～(5)中：

R_ij-k表示单波段k的比值影像中(i,j)像元灰度值，其绝对值在[0，255]范围内；

x_ij-k(t₁)、x_ij-k(t₂)分别表示T₁时相、T₂时相下单波段k的遥感影像中(i,j)像元的灰度值；

D_ij-k表示单波段k的差值影像中(i,j)像元灰度值，其值在[0，255]范围内；

C_ij-k表示单波段k的差异影像中(i,j)像元灰度值，由比值影像R_ij-k、差值影像D_ij-k的乘积均方根复合后，并引入偏差因子a、b调整得到；

a为比值影像和差值影像的像元灰度均值差，表示比值影像所有像元在第k波段的灰度均值，表示差值影像所有像元在第k波段的灰度均值；

b为比值影像和差值影像的像元灰度均方根之比，表示比值影像所有像元在第k波段的灰度值均方根，表示差值影像所有像元在第k波段的灰度值均方根。

步骤4，变化掩膜的构造。

考虑到地物在不同波段的敏感性差异给影像变化检测带来的影响，本发明采用邻域熵权法对R、G、B三个单波段的差异影像进行调权求和处理，权值对应各单波段的贡献度，并结合Ostu阈值分割法确定变化掩膜。

由于地物在多波段遥感影像的不同波段存在波段间的敏感性差异，给遥感影像变化检测带来虚检、漏检等影响，本发明综合考虑各波段对变化检测的贡献度，采用邻域熵权法对单波段差异影像进行调权处理。信息熵度量某随机事件包含的信息量，该波段像元的邻域信息量越大或不确定性越高，该波段对应的权值越大、对变化检测的贡献度也越大。

变化掩膜差异影像Y_ij和邻域熵权w_ij-k计算公式如下所示：

$\left\{\begin{array}{c}{Y}_{ij}={\Σ}_{k=1}^N(w_{ij-k}\×C_{ij-k})\\ w_{ij-k}=\frac{H_{ij-k}}{{\Σ}_{k=1}^N{H}_{ij-k}}\\ H_{ij-k}=-{\Σ}_{t=1}^c{\Σ}_{s=1}^c\left(p_{ts}{\mathrm{logp}}_{ts}\right)\end{array}\right.---\left(6\right)$

式(6)中：

k表示单波段编号，k＝1,2,...N，N为波段数；

Y_ij表示变化掩膜差异影像中(i,j)像元的灰度值；

w_ij-k表示单波段k的差异影像中(i,j)像元的邻域熵权；

C_ij-k表示单波段k的差异影像中(i,j)像元灰度值，可由公式(1)获得；

H_ij-k表示单波段k的差异影像中(i,j)像元的邻域信息熵；

p_ts表示单波段k的差异影像中(i,j)像元邻域窗口内灰度值(该灰度值取值范围为[0,255])出现的概率或频数，t、s分别代表邻域窗口内像元的行号和列号，t、s取值范围[1,c]，c×c为邻域窗口大小。

变化掩膜差异影像的像元灰度取值范围为0～255，要获得变化掩膜还需对变化掩膜差异影像进行阈值分割。本发明采用Ostu法(大津法，又称最大类间方差法)分割获得变化掩膜。其流程为：确定分割阈值t^*，使变化掩膜差异影像中背景和目标两个像素类的类间方差最大。采用分割阈值t^*分割变化掩膜差异影像获得变化掩膜，变化掩膜中像元根据灰度值分成A类和B类，A类像元由灰度值为[0,t^*]的像元组成，B类像元由灰度值为及[t^*+1,L-1]的像元组成，L表示灰度层级。

分割阈值t^*由下列准则函数确定：

$\left\{\begin{array}{c}{w}_0={\mathrm{\Σ}}_{f=0}^t{P}_f,w_1={\mathrm{\Σ}}_{f=t}^{L-1}{P}_f\\ {\mathrm{\μ}}_T={\mathrm{\Σ}}_{f=0}^{L-1}\left({\mathrm{fP}}_{fi}\right),{\mathrm{\μ}}_t={\mathrm{\Σ}}_{f=0}^t\left({\mathrm{fP}}_f\right)\\ {\mathrm{\μ}}_0=\frac{{\mathrm{\μ}}_t}{w_0},{\mathrm{\μ}}_1=\frac{{\mathrm{\μ}}_T-{\mathrm{\μ}}_t}{w_1},{\mathrm{\σ}}_0^2=\frac{{\mathrm{\Σ}}_{f=0}^t{(f-{\mathrm{\μ}}_0)}^2}{w_0}1,{\mathrm{\σ}}_1^2=\frac{{\mathrm{\Σ}}_{f=t}^{L-1}{(f-{\mathrm{\μ}}_1)}^2}{w_1}\\ {\mathrm{\σ}}_W^2=w_0{\mathrm{\σ}}_0^2+w_1{\mathrm{\σ}}_1^2\\ {\mathrm{\σ}}_B^2=w_0{({\mathrm{\μ}}_0-{\mathrm{\μ}}_T)}^2+w_1{({\mathrm{\μ}}_1-{\mathrm{\μ}}_T)}^2=w_0{w}_1{({\mathrm{\μ}}_0-{\mathrm{\μ}}_1)}^2\\ {\mathrm{\σ}}_T^2={\mathrm{\σ}}_W^2+{\mathrm{\σ}}_B^2\end{array}\right.---\left(7\right)$

式(7)中：

P_f表示变化掩膜差异影像中灰度为f的像元数；

w₀表示变化掩膜差异影像中灰度不大于t的像元数，w₁表示变化掩膜差异影像中灰度不小于t的像元数；

μ_T表示变化掩膜差异影像中所有像元灰度值之和，μ_t表示变化掩膜差异影像中灰度不大于t的像元的灰度值之和；

μ₀、分别表示变化掩膜差异影像中灰度不大于t的像元的灰度均值和灰度方差；

μ₁、分别表示变化掩膜差异影像中灰度小于t的像元的灰度均值和灰度方差；

分别代表类内方差、类间方差、总体方差；

L表示灰度层级，一般取256。

步骤5，将变化掩膜与两时相多波段遥感影像之一(记为时相二多波段遥感影像)进行乘积融合，确定变化区域。

变化掩膜影像各像元灰度值为0或1，分别对应非变化类和变化类。将时相二多波段影像与变化掩膜进行乘积融合，灰度值为0的区域即非变化区域，灰度值与融合前保持一致的区域即变化区域，即变化区域的像元灰度值仍为时相二多波段遥感影像的灰度值。

步骤6，利用改进的模糊C均值聚类法对时相二多波段遥感影像的变化区域进行聚类。

变化掩膜能有效确定变化区域，但变化类型信息还需要借助分类算法实现。本发明对模糊C均值聚类法做两点改进，其一：基于区间划分思想，以各区间灰度均值作为初始聚类中心；其二：引入当前像元的邻域像素位置和属性信息，得到当前像元的相似度加权隶属度和加权空间函数。

区间划分思想如下：

类别相同的样本具有类似特征，若两样本的各子向量近似相等，则该两样本必属于同一类别；另外，类别中心必然聚集更多样本，即类别中心样本密度较大。

设影像样本X＝{x₁,x₂…x_N}^T，向量x_n样本中最大值记为x_n ⁺，最小值记为x_n ^-，本发明中每一维向量即变化区域的像元灰度集合。

样本空间为：

W＝W₁×W₂×…W_N＝[x₁ ^-,x₁ ⁺]×[x₂ ^-,x₂ ⁺]×…[x_N ^-,x_N ⁺]           (8)

W₁、W₂、……W_N分别表示向量x₁、x₂、…….x_N的样本空间。

设类别数为C，根据公式(9)将每维输入空间等分为C个子区间，以各子区间样本均值分别作为C个初始聚类中心。

$\left\{\begin{array}{c}{W}_{nq}=\[x_n^{-}+\frac{(q-1)(x_n^{+}-x_n^{-})}{C},x_n^{-}+\frac{q(x_n^{+}-x_n^{-})}{C}\]\\ Z_{nq}=\frac{1}{m_{nq}}\underset{s=1}{\overset{m_{nq}}{\Σ}}{X}^s,n=1,2...N;q=1,2...C\end{array}\right.---\left(9\right)$

式(9)中，

n表示波段号，n∈[1,2…N]，N为波段数；

q表示区间号，即类别号，q∈[1,2…C]，C为类别数；

W_nq表示第n维子空间的第q个子区间；

x_n ⁺、x_n ^-分别表示第n波段遥感影像的灰度最大值和最小值；

Z_nq表示第n维子空间第q个子区间的聚类中心，则第q个聚类中心v_q＝(Z_1q,Z_2q,...Z_Nq)；

m_nq表示第n维子空间第q个子区的像素数；

s表示第n维子空间第q个子区的像素编号，取值为s∈[1,2…m_nq]。

变化区域中像元的相似度加权隶属度如下：

$\left\{\begin{array}{c}{w}_{rt}=1-\frac{x_r-x_t}{x_r+x_t}=\frac{2x_r}{x_r+x_t}\\ u_{rk}^{\left(l\right)}={\mathrm{\Σ}}_{t=r-\frac{c}{2}}^{r+\frac{c}{2}}(w_{rt}\×u_{tk}^{\left(l\right)})\end{array}\right.---\left(10\right)$

式(10)中：

w_rt表示邻域像元t对中心像元r的邻域相似权，中心像元r即邻域窗口的中心像元，邻域像元t即邻域窗口内像元，其取值范围为t、r分别表示邻域像元和中心像元在邻域窗口内的位置编号，位置编号即从邻域窗口第一行第一列像元开始，按照从左向右、从上向下的方式依次编号，例如，对于第i行第j列的像元，其位置编号为i*c+j；

x_r表示中心像元r的灰度值，x_t表示邻域像元t的灰度值；

表示第l次迭代中、中心像元r对第k类聚类中心的隶属度，由其邻域像元调权求和得到；

表示第l次迭代中、邻域像元t对第k类聚类中心的隶属度，其初始值随机获得；

k＝1,2，邻域窗口大小为c×c。

改进的模糊C均值法对变化区域进行变化检测的流程见图2，包括步骤：

(1)根据经验和目测观测结果设定类别数C；设定模糊度m和收敛值ε，模糊度m和收敛值ε为经验值，本具体实施中，令m＝2，ε＝0.01；随机设定初始隶属度矩阵U，隶属度矩阵U由辩护区域中各像元对聚类中心的隶属度；设定初始迭代步数l＝0。

(2)基于区域划分思想，以C个子区间的像元灰度均值分别作为C个初始聚类中心V＝[v₁,v₂,…v_C]，参见公式(9)。

(3)根据公式(10)获得变化区域中邻域像元t对像元r的邻域相似权w_rt，并调整像元r对聚类中心V的隶属度

(4)采用公式(11)修正隶属度矩阵U及聚类中心V：

$\left\{\begin{array}{c}{u}_{rk}^{(l+1)}=\frac{1}{{\Σ}_{y=1}^C{\left(\frac{u_{rk}^{\left(l\right)}\×\left|\right|v_k^{\left(l\right)}-x_r\left|\right|}{u_{ry}^{\left(l\right)}\×\left|\right|v_y^{\left(l\right)}-x_r\left|\right|}\right)}^{\frac{2}{m-1}})},\∀y,\∀k\\ v_k^{(l+1)}=\frac{{\Σ}_{r=1}^n\[{\left(u_{rk}^{\left(l\right)}\right)}^m{x}_r\]}{{\Σ}_{r=1}^n{\left(u_{rk}^{\left(l\right)}\right)}^m}\end{array}\right.---\left(11\right)$

式(11)中：

t、r分别表示邻域像元和中心像元在邻域窗口内的位置编号，k、y为聚类中心编号，C为类别数；

和分别表示第(l+1)和l次迭代中、中心像元r对第k类聚类中心的隶属度；

表示第l次迭代中、中心像元r对第y类聚类中心的隶属度；

x_r表示中心像元r的灰度值；

和表示第l次迭代获得的第k类和第y类聚类中心，y在[1,C]范围内遍历，可以与k相等；

表示第(l+1)次迭代获得的第k类聚类中心；

m表示模糊度，n为邻域窗口内像元数量。

(5)比较隶属度矩阵范数||U^(l+1)-U^(l)||，U^(l+1)、U^(l)分别表示第(l+1)次迭代、第l次迭代下的隶属度矩阵；若||U^(l+1)-U^(l)||≤ε，结束迭代，按照当前聚类中心和隶属度矩阵进行聚类；否则，令l＝l+1，执行步骤(4)。

步骤7，根据时相二遥感影像变化区域的聚类结果，参考时相一遥感影像的地物类型先验知识，确定变化区域的变化类型，即“从什么”变化“到什么”。

本步骤采用人机半自动操作实现。首先，借助ERDAS软件的监督分类模块，获取时相一遥感影像的地物类型先验知识；然后，利用ERDAS软件读取时相二遥感影像模糊聚类结果的属性表，参考时相一遥感影像的地物类型属性值，确定变化区域的变化类型。

这里需指出：第一，时相一遥感影像采用监督法分类，借助人的经验知识来识别地物，具有一定真实度、准确性；第二，“从什么”变化“到什么”的变化类型确定只针对变化区域，该区域由前面变化掩膜(分类前比较法)得到，结果相对可靠。

Claims (7)

1.基于掩膜分类的遥感影像变化检测方法，其特征在于，包括步骤：

步骤1，多波段遥感影像的预处理，即多波段遥感影像的滤波和增强处理；

步骤2，对预处理后的多波段遥感影像进行单波段分离，获得单波段遥感影像；

步骤3，构造两时相单波段遥感影像的单波段差异影像；

步骤4，构造变化掩膜差异影像基于大津法分割变化掩膜差异影像，获得变化掩膜；

上述k表示单波段编号，k＝1,2,...N，N为波段数；Y_ij表示变化掩膜差异影像中(i,j)像元的灰度值；表示单波段k的差异影像中(i,j)像元的邻域熵权，H_ij-k表示单波段k的差异影像中(i,j)像元的邻域信息熵；C_ij-k表示单波段k的差异影像中(i,j)像元灰度值；

步骤5，将变化掩膜与两时相多波段遥感影像之一进行乘积融合，融合后灰度值与融合前灰度值保持一致的区域即变化区域；

步骤6，以变化区域像元灰度值为样本，基于区间划分思想对样本进行分区，以各区样本均值作为初始聚类中心；以邻域像元隶属度调权求和作为中心像元隶属度，各邻域像元对中心像元的邻域相似权利用模糊C均值法对变化区域进行聚类；

上述中，t、r分别表示邻域像元和中心像元在邻域窗口内的位置编号；x_r表示中心像元r的灰度值，x_t表示邻域像元t的灰度值；

步骤7，根据变化区域的聚类结果，结合两时相多波段遥感影像之二的地物类型先验知识，确定变化区域的变化类型。

2.如权利要求1所述的基于掩膜分类的遥感影像变化检测方法，其特征在于：

步骤2具体为：

采用加权求和法将多波段遥感影像灰度化，当R、G、B单波段某一分量的权值为1，其余为0时，经加权求和即获得R、G、B单波段遥感影像。

3.如权利要求1所述的基于掩膜分类的遥感影像变化检测方法，其特征在于：

步骤3中结合差值法和比值法构造单波段差异影像，具体为：

根据公式构造两时相单波段遥感影像的单波段差异影像，其中：

C_ij-k表示单波段k的差异影像中(i,j)像元灰度值；

R_ij-k表示单波段k的比值影像中(i,j)像元灰度值；

D_ij-k表示单波段k的差值影像中(i,j)像元灰度值；

a为比值影像所有像元在第k波段的灰度均值和差值影像中所有像元在第k波段的灰度均值之差；

b为比值影像所有像元在第k波段的灰度值均方根和差值影像所有像元在第k波段的灰度值均方根之差。

4.如权利要求1所述的基于掩膜分类的遥感影像变化检测方法，其特征在于：

步骤4中，所述的基于大津法分割变化掩膜差异影像，具体为：

基于大津法单阈值分割法，将使变化掩膜差异影像的背景像元灰度和目标像元灰度的类间方差最大的阈值作为分割阈值；采用分割阈值分割变化掩膜差异影像，获得变化掩膜。

5.如权利要求1所述的基于掩膜分类的遥感影像变化检测方法，其特征在于：

步骤6进一步包括子步骤：

6.1根据经验和目测观测结果设定类别数C，根据经验设定模糊度m和收敛值ε，随机设定初始隶属度矩阵，隶属度矩阵由变化区域中各像元对聚类中心的隶属度构成；令初始迭代步数l＝0；

6.2以变化区域像元灰度值为样本，基于区间划分思想将样本划分为C个子区间，以各子区间样本均值作为初始聚类中心；

6.3将邻域像元t当前隶属度的加权和作为中心像元r隶属度即邻域像元t的当前隶属度根据当前隶属度矩阵U获得，其权值 $w_{\mathrm{rt}}=\frac{2x_r}{x_r+x_t};$

6.4采用公式 $\left\{\begin{array}{c}{u}_{\mathrm{rk}}^{(l+1)}=\frac{1}{{\mathrm{\Σ}}_{y=1}^C{\left(\frac{u_{\mathrm{rk}}^{\left(l\right)}\×\left|\right|v_k^{\left(l\right)}-x_r\left|\right|}{u_{\mathrm{ry}}^{\left(l\right)}\×\left|\right|v_y^{\left(l\right)}-x_r\left|\right|}\right)}^{\frac{2}{m-1}})},\∀y,\∀k\\ v_k^{(l+1)}=\frac{{\mathrm{\Σ}}_{r=1}^n\left[{\left(u_{\mathrm{rk}}^{\left(l\right)}\right)}^m{x}_r\right]}{{\mathrm{\Σ}}_{r=1}^n{\left(u_{\mathrm{rk}}^{\left(l\right)}\right)}^m}\end{array}\right.$ 修正隶属度矩阵及聚类中心，其中，t、r分别表示邻域像元和中心像元在邻域窗口内的位置编号，k、y为聚类中心编号；和分别表示第(l+1)和l次迭代中、中心像元r对第k类聚类中心的隶属度；表示第l次迭代中、中心像元r对第y类聚类中心的隶属度；x_r表示中心像元r的灰度值；和表示第l次迭代获得的第k类和第y类聚类中心，y在[1,C]范围内遍历；表示第(l+1)次迭代获得的第k类聚类中心；n为邻域窗口内像元数；

6.5比较隶属度矩阵范数||U^(l+1)-U^(l)||，U^(l+1)、U^(l)分别表示第(l+1)次迭代、第l次迭代获得的隶属度矩阵；若||U^(l+1)-U^(l)||≤ε，结束迭代，按照当前聚类中心和隶属度矩阵进行聚类；否则，令l＝l+1，执行步骤6.4。

6.如权利要求1所述的基于掩膜分类的遥感影像变化检测方法，其特征在于：

步骤7进一步包括：

7.1对两时相多波段遥感影像之二进行监督分类，获得地物类型先验知识；

7.2结合两时相多波段遥感影像之一的聚类结果，采用人工方式确定变化区域的变化类型。

7.基于掩膜分类的遥感影像变化检测系统，其特征在于，包括：

预处理模块，用来多波段遥感影像的预处理，即多波段遥感影像的滤波和增强处理；

单波段分离模块，用来对预处理后的多波段遥感影像进行单波段分离，获得单波段遥感影像；

单波段分离模块，用来构造两时相单波段遥感影像的单波段差异影像；

变化掩膜构造模块，用来构造变化掩膜差异影像基于大津法分割变化掩膜差异影像，获得变化掩膜；

上述k表示单波段编号，k＝1,2,...N，N为波段数；Y_ij表示变化掩膜差异影像中(i,j)像元的灰度值；表示单波段k的差异影像中(i,j)像元的邻域熵权，H_ij-k表示单波段k的差异影像中(i,j)像元的邻域信息熵；C_ij-k表示单波段k的差异影像中(i,j)像元灰度值；

乘积融合模块，用来将变化掩膜与两时相多波段遥感影像之一进行乘积融合，融合后灰度值与融合前灰度值保持一致的区域即变化区域；

模糊C均值聚类模块，用来以变化区域像元灰度值为样本，基于区间划分思想对样本进行分区，以各区样本均值作为初始聚类中心；以邻域像元隶属度调权求和作为中心像元隶属度，各邻域像元对中心像元的邻域相似权利用模糊C均值法对变化区域进行聚类；

上述中，t、r分别表示邻域像元和中心像元在邻域窗口内的位置编号；x_r表示中心像元r的灰度值，x_t表示邻域像元t的灰度值；

变化检测模块，用来根据变化区域的聚类结果，结合两时相多波段遥感影像之二的地物类型先验知识，确定变化区域的变化类型。