CN104867150A - 遥感影像模糊聚类的波段修正变化检测方法及*** - Google Patents

Abstract

本发明提供了一种遥感影像模糊聚类的波段修正变化检测方法及***，包括步骤：步骤1，多波段遥感影像的预处理，即多波段遥感影像的滤波和增强处理；步骤2，对预处理后的多波段遥感影像进行单波段分离，获得单波段遥感影像；步骤3，构造两时相单波段遥感影像的单波段差异影像；步骤4，基于预处理后的两时相多波段遥感影像构造多波段联合影像，并获得多波段联合影像各像元的邻域互信息量；步骤5，修正多波段联合影像像元的邻域互信息量，得到修正影像；步骤6，利用模糊C均值法对修正影像进行变化检测。本发明检测精度高，抗噪声干扰强、自动化程度高。

Description

遥感影像模糊聚类的波段修正变化检测方法及***

技术领域

本发明属于摄影测量与遥感影像应用技术领域，具体涉及一种遥感影像模糊聚类的波段修正变化检测方法及***。

背景技术

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

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

文中涉及如下参考文献：

[1]Song C,Woodcock CE,Seto KC,Lenney MP,Macomber SA.Classification andchange detection using Landsat TM data:when and how to correct atmospheric effects[J].Remote sensing of Environment.2001,75(2):230-44.

[2]Mas J-F.Monitoring land-cover changes:a comparison of change detectiontechniques[J].International journal of remote sensing.1999,20(1):139-52.

[3]Jensen JR.Introductory digital image processing:a remote sensing perspective:Prentice-Hall Inc.1996.

[4]申邵洪,赖祖龙,万幼川.基于融合的高分辨率遥感影像变化检测[J].测绘通报.2009(3):16-9.

[5]Jensen JR.Introductory digital image processing:a remote sensing perspective:Prentice-Hall Inc.；1996.

[6]Xie XL,Beni G.A validity measure for fuzzy clustering[J].IEEE Transactions onpattern analysis and machine intelligence.1991,13(8):841-7.

发明内容

针对现有技术中存在的问题，本发明提供了一种高实时性、高检测精度、高自动化程度的遥感影像模糊聚类的波段修正变化检测方法。

本发明将模糊聚类的思想引入遥感影像变化检测中，在利用改进FCM进行分类处理基础上，采用一种单波段对比熵权修正多波段联合互信息量的波段修正法来提高变化检测精度。

本发明主要包括两大阶段，第一阶段：单波段对比熵权修正多波段联合互信息量；第二阶段：利用改进的模糊C均值聚类法对修正影像进行变化检测。

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

另一方面，为了综合利用单波段与多波段丰富的光谱及空间结构信息，且有效降低由于不同传感器或不同时相影像存在的地理和散射差异(灰度层级整体偏离)导致对分类和变化检测的影响，本发明以两个时相的多波段影像的像元点对为基础，计算邻域互信息量作为多波段联合处理结果。互信息量基于一对像元灰度值的统计特征，不需要选择基准点或提取图像的纹理特征，且不必考虑图像中的灰度值的线性关系，因此在两图像相似度的比较及变化信息的检测上具有更大的灵活度和准确性。目前常用的MI(互信息量)类型包括归一化互信息量NMI、图像梯度互信息量GMI、区域互信息量RMI。传统的邻域互信息量虽然以邻域为基础，但目的是寻找最优的配准变换参数，使两幅图像间的互信息最大，属于整图幅范围影像配准范畴。而变化检测往往发生在某些像元或小区域，变化的程度与像元点对的互信息量大小紧密相关，具体表现为互信息量越大，像元点对的相关性越大，对应的变化程度越小。

第二阶段，可以看成特殊的分类过程，即将波段修正影像聚类为变化、非变化两类。模糊C均值聚类算法(FCM)是一种典型的“软聚类”，算法简单、收敛速度快，隶属度的引进使得处理混合像元行之有效。但该算法存在一些不足：如算法性能过分依赖初始聚类中心、类别数不能自动确定、忽略空间信息导致局部最优等。本发明从两个方面对FCM算法进行改进：其一，结合大津法单阈值分割思想，求出满足类间方差最大的最优阈值t^*，并以该阈值分割下两类像元的灰度均值作为初始聚类中心。其二，引入邻域像元的位置和属性信息得到相似度加权隶属度，并进一步得到加权空间函数，使模糊C均值法具有更强的抗噪能力，有效避免局部最优，改善分类效果。最后利用改进的模糊C均值算法对修正影像进行变化检测。

本发明的技术方案如下：

一、遥感影像模糊聚类的波段修正变化检测方法，包括步骤：

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

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

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

步骤4，基于预处理后的两时相多波段遥感影像构造多波段联合影像，并获得多波段联合影像各像元的邻域互信息量；

步骤5，根据公式修正多波段联合影像像元的邻域互信息量，得到修正影像；其中：k表示单波段编号，k＝1,2,...N，N为波段数；M_ij表示修正影像(i,j)像元灰度值；w_ij-k为表示单波段k的差异影像中(i,j)像元的邻域熵权，H_ij-k表示单波段k的差异影像中(i,j)像元的邻域信息熵；C_ij-k表示单波段k的差异影像中(i,j)像元灰度；RMI_ij表示多波段联合影像(i,j)像元的邻域互信息量；Mod_ij表示(i,j)像元对应N个单波段的熵值模，

步骤6，以修正影像中像元灰度值为样本，令类别c＝2，以修正影像阈值分割后两类像元的灰度均值作为初始聚类中心；以邻域像元隶属度调权求和作为中心像元隶属度，各邻域像元对中心像元的邻域相似权利用模糊C均值法对修正影像进行变化检测；

上述的阈值分割中，基于大津法单阈值分割，将使修正影像的背景像元灰度和目标像元灰度的类间方差最大的阈值作为分割阈值；

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

步骤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波段的灰度值均方根之差。

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

6.1设定类别C＝2，根据经验设定模糊度m和收敛值ε，随机设定初始隶属度矩阵，隶属度矩阵由修正影像中各像元对聚类中心的隶属度构成；令初始迭代步数l＝0；

6.2基于大津法单阈值分割法，将使修正影像的背景和目标像元灰度的类间方差最大的阈值作为分割阈值，以修正影像阈值分割后两类像元的灰度均值作为初始聚类中心

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}})}\text{,\∀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为聚类中心编号，C为类别数；和分别表示第(l+1)和l次迭代中、中心像元r对第k类聚类中心的隶属度；表示第l次迭代中、中心像元r对第y类聚类中心的隶属度；x_r表示中心像元r的灰度值；和表示第l次迭代获得的第k类和第y类聚类中心，y在[1,C]范围内遍历，可以与k相等；表示第(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。

二、遥感影像模糊聚类的波段修正变化检测***，包括：

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

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

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

邻域互信息量获得模块，用来基于预处理后的两时相多波段遥感影像构造多波段联合影像，并获得多波段联合影像各像元的邻域互信息量；

修正影像获得模块，用来根据公式修正多波段联合影像像元的邻域互信息量，得到修正影像；其中：k表示单波段编号，k＝1,2,...N，N为波段数；M_ij表示修正影像(i,j)像元灰度值；w_ij-k为表示单波段k的差异影像中(i,j)像元的邻域熵权，H_ij-k表示单波段k的差异影像中(i,j)像元的邻域信息熵；C_ij-k表示单波段k的差异影像中(i,j)像元灰度；RMI_ij表示多波段联合影像(i,j)像元的邻域互信息量；Mod_ij表示(i,j)像元对应N个单波段的熵值模，

模糊C均值聚类模块，用来以修正影像中像元灰度值为样本，令类别c＝2，以修正影像阈值分割后两类像元的灰度均值作为初始聚类中心；以邻域像元隶属度调权求和作为中心像元隶属度，各邻域像元对中心像元的邻域相似权利用模糊C均值法对修正影像进行变化检测；

上述的阈值分割中，基于大津法单阈值分割，将使修正影像的背景像元灰度和目标像元灰度的类间方差最大的阈值作为分割阈值；

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

传统的遥感影像变化检测法大多从单像元、单决策的角度出发，对空间结构复杂、细节丰富的中高分辨率遥感影像容易出现“同物异谱”、“异物同谱”所导致的波段敏感性等问题。如何快速、准确地提取中高分辨率遥感影像的变化信息则显得至关重要。

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

1、检测精度高，抗噪声干扰强。

(1)在单波段对比差异影像构造环节，本发明采用差值/比值复合法，在差值、比值乘积融合的基础上，引入偏差因子，对灰度层级进行调整。另外，考虑到地物在不同波段的敏感性差异给变化检测带来的影响，采用邻域熵权法、邻域互信息量法分别对单波段差异影像和多波段影像对进行调权处理，权值表示该波段的贡献度。经过上述处理，可增强遥感影像的变化信息、且能反应地物辐射能量变化，复合差异影像的振幅及结构信息也增强，更有利于变化检测。

(2)采用改进的模糊C均值聚类法对修正影像进行变化检测。其中，改进的模糊C均值聚类法根据类间方差最大最优原理，确定初始聚类中心，且引入邻域空间信息得到相似度加权隶属度和加权空间函数，使模糊C均值聚类法具有更强的抗噪能力，从而有效避免局部最优，改善分类效果，提高变化检测精度。

2、自动化程度高，全过程自动检测。

整个影像变化检测过程无需人工干预，能达到自动检测。

附图说明

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

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

具体实施方式

本发明主要设计两部分内容，即(1)单波段对比熵权修正多波段联合互信息量和(2)基于模糊C均值聚类法的影像变化检测。本发明技术方案主要包括差值/比值复合法构造单波段对比差异影像、单波段对比熵权修正多波段互信息量、基于模糊C均值聚类法的影像变化检测。本发明可缓解同物异谱、异物同谱引起的单波段敏感性差异对变化检测虚检、漏检的影响，能有效抑制孤立噪声干扰、避免局部最优等，具有较优的时效性和精度。

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

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

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

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

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

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

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

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

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

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

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

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

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

$b=\frac{{\mathrm{\σ}}_{D_{\mathrm{ij}-k}}}{{\mathrm{\σ}}_{R_{\mathrm{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三个单波段的差异影像进行调权处理，权值对应各单波段的贡献度。

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

邻域熵权计算公式如下：

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

式(6)中：

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

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

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

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

步骤5：针对预处理后的两时相多波段遥感影像构造多波段联合影像，获得邻域互信息量。

互信息量基于一对像元灰度值的统计特征，不需要选择基准点或提取影像的纹理特征，且不必考虑影像灰度值的线性关系，因此在两影像相似度的比较及变化信息的检测上具有更大的灵活度和准确性。变化的程度与像元点对的互信息量大小紧密相关，具体表现为互信息量越大，像元点对的相关性越大，对应的变化程度越小。

因此，为充分利用图像的灰度和空间信息，有效降低由不同传感器或不同时相影像存在的地理和散射差异(即，灰度层级整体偏离)导致对分类和变化检测的影响。本发明以两个时相多波段遥感影像的像元点对为基础，构造当前像元的邻域像元的灰度值向量P₁、P₂，并计算当前像元的邻域互信息量作为多波段联合处理结果。

邻域互信息量的计算方法如下：

$\left\{\begin{array}{c}{\mathrm{RMI}}_{\mathrm{ij}}=H_{\mathrm{ij}}\left(C_{X1}\right)+H_{\mathrm{ij}}\left(C_{X2}\right)-H_{\mathrm{ij}}\left(C_{X12}\right)\\ C_{X1}=\frac{1}{d}{(P_1-m_1)}^{^2},C_{X2}=\frac{1}{d}{(P_2-m_2)}^{^2},C_{X12}=\frac{1}{d}(P_1-m_1)\×(P_2-m_2)\\ m_1=\frac{1}{d}{\mathrm{\Σ}}_d{P}_1,m_2=\frac{1}{d}{\mathrm{\Σ}}_d{P}_2\\ H_{\mathrm{ij}}\left(\mathrm{\Σd}\right)=\log \left({\left(2\mathrm{\πe}\right)}^{d/2}\det {\left(2\mathrm{\πe}\right)}^{d/2}\right)\end{array}\right.---\left(7\right)$

式(7)中：

RMI_ij表示多波段联合影像中(i,j)像元的邻域互信息量；

H_ij(·)表示求邻域信息熵的运算；

P₁为第一时相多波段遥感影像中(i,j)像元的邻域像元灰度值x_1i按像元位置顺序组合的向量[x₁₁..x_1i..x_1d]；

P₂为第二时相遥感影像中(i,j)像元的邻域像元灰度值x_2i按像元位置顺序组合的向量[x₂₁..x_2i..x_2d]；

C_X1、C_X2分别为第一时相遥感影像和第二时相遥感影像的方差矩阵，C_X12为协方差矩阵；

m₁、m₂分别向量P₁、P₂中元素均值，d＝c×c，表示为邻域大小；

e表示自然对数。

步骤6：采用单波段邻域熵权和单波段差异影像修正多波段联合影像像元的邻域互信息量，得到修正影像。

函数修正关系如下：

$\left\{\begin{array}{c}{M}_{\mathrm{ij}}=\frac{{\mathrm{RMI}}_{\mathrm{ij}}}{{\mathrm{Mod}}_{\mathrm{ij}}}\×{\mathrm{\Σ}}_{k=1}^N(w_{\mathrm{ij}-k}\×C_{\mathrm{ij}-k})\\ {\mathrm{Mod}}_{\mathrm{ij}}=\sqrt{{\mathrm{\Σ}}_{k=1}^N{\left(H_{\mathrm{ij}-k}\right)}^2}\end{array}\right.---\left(8\right)$

式(8)中：

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

M_ij表示修正影像(i,j)像元灰度值；

RMI_ij表示多波段联合影像(i,j)像元的邻域互信息量；

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

w_ij-k为C_ij-k的邻域熵权，表示单波段k的贡献度；

Mod_ij表示(i,j)像元对应N个单波段的熵值模；

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

步骤7：利用改进的模糊C均值聚类法对修正影像进行变化检测。

遥感影像变化检测可看成特殊的分类，即将修正影像分为变化、非变化两大类。模糊C均值聚类引入隶属度的概念，是一种软聚类算法，对混合像元具有很好的处理效果。

本发明对模糊C均值聚类法做两点改进，其一：结合大津法单阈值分割思想，求出满足类间方差最大的分割阈值t^*，并根据分割阈值t^*对修正影像进行分割，并将分割后两类像元的灰度均值作为初始聚类中心。其二：引入当前像元的邻域像元的位置和属性信息得到当前像元的相似度加权隶属度，并进一步得到加权空间函数，使模糊C均值聚类法具有更强的抗噪能力，有效避免局部最优，改善分类效果。

阈值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(f{P}_{\mathrm{fi}}\right),{\mathrm{\μ}}_t={\mathrm{\Σ}}_{f=0}^t\left(f{P}_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},{\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(9\right)$

式(9)中：

P_f表示修正影像中灰度为f的像元数；

w₀表示修正影像中灰度不大于t的像元数，w₁表示修正影像中灰度不小于t的像元数；

μ_T表示修正影像中所有像元灰度值之和，μ_t表示修正影像中灰度不大于t的像元的灰度值之和；

μ₀、分别表示修正影像中灰度不大于t的像元的灰度均值和灰度方差；

μ₁、分别表示修正影像中灰度小于t的像元的灰度均值和灰度方差；

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

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

使类间方差最大的阈值t即分割阈值t^*。

修正影像中像元的相似度加权隶属度如下：

$\left\{\begin{array}{c}{w}_{\mathrm{rt}}=1-\frac{x_r-x_t}{x_r+x_t}=\frac{2x_r}{x_r+x_t}\\ u_{\mathrm{rk}}^{\left(l\right)}={\mathrm{\Σ}}_{t=r-\frac{c}{2}}^{r+\frac{c}{2}}(w_{\mathrm{rt}}\×u_{\mathrm{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＝2，分别对应变化类和非变化类；设定模糊度m和收敛值ε，模糊度m和收敛值ε为经验值，本具体实施中，令m＝2，ε＝0.01；随机设定初始隶属度矩阵U，隶属度矩阵U由修正影像中各像元对聚类中心的隶属度；设定初始迭代步数l＝0。

(2)利用Ostu阈值分割法获得分割阈值t^*，以分割阈值t^*分割下两类像元的灰度均值作为初始聚类中心V＝[v₁,v₂]。

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

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

式(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)。

Claims (5)

1.遥感影像模糊聚类的波段修正变化检测方法，其特征在于，包括步骤：

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

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

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

步骤4，基于预处理后的两时相多波段遥感影像构造多波段联合影像，并获得多波段联合影像各像元的邻域互信息量；

步骤5，根据公式修正多波段联合影像像元的邻域互信息量，得到修正影像；其中：k表示单波段编号，k＝1,2,...N，N为波段数；M_ij表示修正影像(i,j)像元灰度值；w_ij-k为表示单波段k的差异影像中(i,j)像元的邻域熵权，H_ij-k表示单波段k的差异影像中(i,j)像元的邻域信息熵；C_ij-k表示单波段k的差异影像中(i,j)像元灰度；RMI_ij表示多波段联合影像(i,j)像元的邻域互信息量；Mod_ij表示(i,j)像元对应N个单波段的熵值模，

步骤6，以修正影像中像元灰度值为样本，令类别c＝2，以修正影像阈值分割后两类像元的灰度均值作为初始聚类中心；以邻域像元隶属度调权求和作为中心像元隶属度，各邻域像元对中心像元的邻域相似权利用模糊C均值法对修正影像进行变化检测；

上述的阈值分割中，基于大津法单阈值分割，将使修正影像的背景像元灰度和目标像元灰度的类间方差最大的阈值作为分割阈值；

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

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所述的遥感影像模糊聚类的波段修正变化检测方法，其特征在于：

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

6.1设定类别C＝2，根据经验设定模糊度m和收敛值ε，随机设定初始隶属度矩阵，隶属度矩阵由修正影像中各像元对聚类中心的隶属度构成；令初始迭代步数l＝0；

6.2基于大津法单阈值分割法，将使修正影像的背景和目标像元灰度的类间方差最大的阈值作为分割阈值，以修正影像阈值分割后两类像元的灰度均值作为初始聚类中心

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|}{)}^{\frac{2}{m-1}}\right)},\∀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为聚类中心编号，C为类别数；和分别表示第(l+1)和l次迭代中、中心像元r对第k类聚类中心的隶属度；表示第l次迭代中、中心像元r对第y类聚类中心的隶属度；x_r表示中心像元r的灰度值；和表示第l次迭代获得的第k类和第y类聚类中心，y在[1,C]范围内遍历，可以与k相等；表示第(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。

5.遥感影像模糊聚类的波段修正变化检测***，其特征在于，包括：

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

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

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

邻域互信息量获得模块，用来基于预处理后的两时相多波段遥感影像构造多波段联合影像，并获得多波段联合影像各像元的邻域互信息量；

修正影像获得模块，用来根据公式修正多波段联合影像像元的邻域互信息量，得到修正影像；其中：k表示单波段编号，k＝1,2,...N，N为波段数；M_ij表示修正影像(i,j)像元灰度值；w_ij-k为表示单波段k的差异影像中(i,j)像元的邻域熵权，H_ij-k表示单波段k的差异影像中(i,j)像元的邻域信息熵；C_ij-k表示单波段k的差异影像中(i,j)像元灰度；RMI_ij表示多波段联合影像(i,j)像元的邻域互信息量；Mod_ij表示(i,j)像元对应N个单波段的熵值模，

模糊C均值聚类模块，用来以修正影像中像元灰度值为样本，令类别c＝2，以修正影像阈值分割后两类像元的灰度均值作为初始聚类中心；以邻域像元隶属度调权求和作为中心像元隶属度，各邻域像元对中心像元的邻域相似权利用模糊C均值法对修正影像进行变化检测；

上述的阈值分割中，基于大津法单阈值分割，将使修正影像的背景像元灰度和目标像元灰度的类间方差最大的阈值作为分割阈值；

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