CN107346549A - 一种利用遥感影像多特征的多类别变化动态阈值检测方法 - Google Patents

Abstract

本发明公开了一种利用遥感影像多特征的多类别变化动态阈值检测方法，首先分别提取两时期影像的光谱变化特征和光谱梯度差异特征，按波段进行叠加，并对其进行主成分分析，选择前两个主成分分量作为最终的变化特征，计算变化强度和变化方向。然后，在极坐标框架下计算变化强度阈值，将像素分为未变化和变化两部分，对变化像素的变化方向设置阈值，区分变化类别。在此基础上，对初始变化检测结果进行分析，对所有变化类别依次生成仅包含单一变化的像素集，再获取针对各变化类别的最优变化强度阈值，从而生成基于动态阈值的变化检测结果。最后，利用多相水平集算法对检测结果进行后处理，提高了多类别变化的检测精度。

Description

一种利用遥感影像多特征的多类别变化动态阈值检测方法

技术领域

本发明涉及遥感影像多类别检测方法，特别涉及一种利用遥感影像多特征的多类别变化动态阈值检测方法，属于遥感影像多类别检测技术领域。

背景技术

遥感影像变化检测技术已成为一种有效的地表变化监测方法，广泛应用于土地利用、土地覆盖、环境监测、灾害评估和城镇发展监测等领域。

现有的多类别变化检测方法可分为监督和非监督两类，监督变化检测是利用分类器对不同时期的影像进行分类后，通过比较分类结果发现变化。该类方法需要参考数据进行分类器的训练，而在实际应用中采集参考数据需要耗费大量的人力物力，且变化检测的精度受分类精度影响严重。多类别变化的非监督检测方法在通过分析不同时期地物光谱的变化，实现变化的检测及不同变化类别的区分，无需要参考数据，且检测精度较高。然而，现有的多类别变化的非监督检测方法仅以用光谱变化特征，受光照、大气辐射和季节等因素影响无法准确描述所有类别变化，设置单一阈值进行多类别变化的检测，也因忽略了各变化类别的特性而无法精确检测所有类别的变化。

发明内容

本发明所要解决的技术问题是：提供一种利用遥感影像多特征的多类别变化动态阈值检测方法，利用多种特征组合，减弱光照、大气等因素的影响，并充分考虑不同变化类别的特性，提出一种动态阈值的多类别变化检测模型，准确检测所有类别的变化，提高多类别变化的检测精度。

本发明为解决上述技术问题采用以下技术方案：

一种利用遥感影像多特征的多类别变化动态阈值检测方法，包括如下步骤：

步骤1，对前、后两时期的遥感影像进行预处理，包括几何校正、配准、辐射校正；

步骤2，计算经步骤1预处理后的前、后两时期遥感影像的光谱变化特征和光谱梯度差异特征，将光谱变化特征和光谱梯度差异特征按波段进行叠加，利用主成分分析方法对叠加后的特征进行主成分分析，选取前两个主成分分量作为变化特征，根据变化特征计算变化强度和变化方向；

步骤3，在极坐标框架下，对变化强度和变化方向进行空间分析，针对变化强度计算变化强度阈值，根据变化强度阈值检测未变化和变化类别；针对变化类别的变化方向设置阈值，根据阈值区分不同的变化类别，生成初步多类别变化检测结果；

步骤4，分析初步多类别变化检测结果，针对各变化类别，构建仅包含当前变化类别的像素集，并设置动态的变化强度阈值，重新判定各变化类别中的变化像素与未变化像素，生成基于动态阈值的多类别变化检测结果。

作为本发明的一个进一步方案，该方法还包括步骤5，利用多相水平集算法，对步骤4生成的基于动态阈值的多类别变化检测结果进行处理，通过集成空间信息同时去除虚检测和漏检测像素。

作为本发明的一个优选方案，步骤2所述光谱变化特征和光谱梯度差异特征计算公式分别为：

X_D＝X₂-X₁

G＝G₂-G₁

其中，X_D表示光谱变化特征，X₁、X₂分别表示经预处理后的前、后时期遥感影像，G表示光谱梯度差异特征，G₁、G₂分别表示经处理后的前、后时期遥感影像光谱梯度。

作为本发明的一个优选方案，步骤2所述变化强度和变化方向计算公式分别为：

其中，ρ、α分别表示变化强度、变化方向，P₁、P₂表示经主成分分析后选取前两个主成分分量。

作为本发明的一个优选方案，步骤4所述针对各变化类别，构建仅包含当前变化类别的像素集具体为：

获取各变化类别在极坐标空间中的分布，构建仅包含单一变化类别的像素集，针对变化类别k，找出该类别中变化方向的最小值T_α,k和最大值T_α,k+1，将极坐标空间中变化方向值在(T_α,k,T_α,k+1)内的变化像素集S_k和未变化像素集合并，生成仅包含变化类别k的像素集S：

其中，α(i,j)表示极坐标空间中变化像素或未变化像素的变化方向值。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

本发明针对现有多类别变化非监督检测方法存在的不足，提出一种利用遥感影像多特征的多类别变化动态阈值检测方法，旨在利用多种特征组合，减弱光照、大气等因素的影响，并充分考虑不同变化类别的特性，提出一种动态阈值的多类别变化检测模型，准确检测所有类别的变化，提高多类别变化的检测精度。

附图说明

图1是本发明方法的流程图。

图2是本发明方法中光谱梯度差异计算示意图。

图3是本发明方法中第k类变化在极坐标空间中的表示。

图4是本发明方法中基于水平集的影像分割示意图；其中，(a)为单相水平集，(b)为二相水平集。

图5是本发明实例影像与参考变化；其中，(a)为前时期影像第4波段，(b)为后时期影像第4波段，(c)为参考变化。

图6是本发明基于不同特征的变化强度及变化方向在极坐标空间的分布；其中，(a)为光谱变化特征，(b)为光谱梯度差异特征，(c)为叠加特征的前两个主成分分量。

图7是传统方法与本发明方法的检测结果对比；其中，(a)为光谱变化特征的检测结果，(b)为光谱梯度差异的检测结果，(c)为本发明的检测结果。

具体实施方式

下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

如图1所示，一种利用遥感影像多特征的多类别变化动态阈值检测方法，包括以下步骤：

步骤1：对两时期遥感影像进行几何校正、配准、辐射校正等预处理。

步骤2：分别计算两时期影像的光谱变化特征及光谱梯度差异特征，并通过波段叠加生成多特征。在此基础上，利用主成分分析方法对叠加后的特征进行主成分分析，选取合适的主成分作为最终的变化特征，并计算变化强度和变化方向，进行多类别变化检测。

采用影像差分技术，根据式(1)从前时期影像X₁和后时期影像X₂生成光谱变化特征X_D：

X_D＝X₂-X₁ (1)

利用式(2)计算前时期影像光谱梯度G₁，如图2所示。

G₁＝{G_1,b+1(i,j)-G_1,b(i,j)|1≤i≤I,1≤j≤J,1≤b≤B-1} (2)

其中，I、J和B分别为影像行、列数和波段数，G_1,b(i,j)和G_1,b+1(i,j)分别表示前时期影像第b和b+1波段在(i,j)位置上的像素灰度值。同理，计算后时期影像光谱梯度G₂，在此基础上根据下式计算两时期影像光谱梯度差异G：

G＝G₂-G₁ (3)

将光谱变化特征与光谱梯度差异特征按波段进行叠加，并利用主成分分析算法对叠加特征进行主成分分析，本发明选择前两个主成分分量作为最终的变化特征P＝{P₁,P₂}，用于变化检测。

然后利用式(4)，针对提取的主成分分量计算变化强度ρ和变化方向α：

步骤3：生成初步多类别变化检测结果。在极坐标框架下，对变化强度和方向进行空间分析，针对变化强度设置阈值，检测未变化和变化类别；针对变化类别的变化方向设置阈值，区分不同类别的变化。

计算变化强度阈值时，假设影像中未变化和变化像素的变化强度满足混合高斯分布，可用下式表示：

p(ρ)＝P(ω_n)p(ρ/ω_n)+P(ω_c)p(ρ/ω_c) (5)

其中，P(ω_n)、P(ω_c)和p(ρ/ω_n)、p(ρ/ω_c)分别是未变化类别ω_n和变化类别ω_c的先验概率和概率密度函数，且先验概率密度函数可用高斯函数表示。在此基础上，利用最大期望算法估计变化和未变化类别的均值和方差。以未变化类别为例，利用下式计算未变化类别的后验概率P(ω_n/ρ)：

然后，利用如下公式重新估计均值μ_n、方差和先验概率：

其中，下标t和t+1分别表示当前和下次迭代次数。重复式(6)至(9)，直到收敛，同理用上述方法估计变化类别的均值μ_c和方差根据Bayesian最小错误准则，依据下式计算变化强度阈值T_ρ，并区分变化和未变化类别。

然后，针对变化像素的变化方向，利用模糊C均值聚类(FCM)算法，通过在极坐标框架中像素的空间分布，确定聚类数目，对变化像素进行聚类，区分变化类别。FCM聚类算法采用各个样本与所在类均值的差值平方和最小准则，通过迭代和更新隶属度矩阵U和聚类中心V，使目标函数F达到最小，实现最优聚类，目标函数如下：

其中，u_wk(u_wk∈[0,1]且)表示数据x_w对于类别k的隶属度，U＝{u_wk}是隶属度矩阵，x_w表示变化像素w的变化方向，V＝{v₁,v₂,…,v_K}表示聚类中心点集，q∈[1,+∞)为加权指数，用来控制聚类结果的模糊程度，当q＝1时，模糊聚类变为传统的C均值聚类，通常情况下q＝2时计算简单且效果较理想，故本发明中设q＝2。FCM聚类方法的具体步骤为：

步骤1)：设定聚类类别数K，聚类中心V＝{v₁,v₂,…,v_K}初始化；

步骤2)：按下式计算模糊隶属度矩阵：

步骤3)：按下式更新聚类中心：

步骤4)：判定式(11)是否收敛。当||V^t+1-V^t||<ε(t表示迭代次数，ε>0是算法停止阈值)时，式(11)收敛，算法停止迭代；否则重复步骤2)和3)，直到式(11)收敛。最后，经过FCM处理后，可得到各像素的隶属度信息，按最大隶属度原则，将变化像素分到相应的类别，实现变化类别区分。

步骤4：建立动态阈值模型，进行多类别变化检测。分析步骤3所得的初步变化检测结果，针对各变化类别，构建新的仅包含一类变化的像素集，并设置动态的变化强度阈值，重新判定各变化类别中的变化与未变化像素。最终，生成基于动态阈值的多类别变化检测结果。

通过分析步骤3获取的各初始变化类别在极坐标空间中的分布，构建仅包含单一类别变化的像素集，如图3所示。针对变化类别k，T_α,k和T_α,k+1分别表示该类别中变化方向的最小和最大值，将变化方向值在此范围内的变化像素集S_k(图3中T_α,k和T_α,k+1之间的白色部分)和未变化像素集(图3中将T_α,k和T_α,k+1向圆心延伸到达圆心后形成的灰色部分)合并，生成仅包含变化类别k的像素集S(白色与灰色部分合并成扇形部分)。

假设像素集S的变化强度ρ为混合高斯分布，根据步骤3中式(5)至式(10)，计算针对变化类别k的最优阈值T_k。同理，针对所有变化类别，逐一生成仅包含单一变化类别的像素集，并针对变化强度计算变化阈值。最终针对所有变化类别，根据下式，实现基于动态阈值的多类别变化非监督检测，提高变化检测的精度和自动化程度。

步骤5：利用多相水平集算法，对步骤4中生成的基于动态阈值的多类别变化检测检测进行后处理，通过集成空间信息同时去除虚检和漏检测像素，进一步提高多类别变化非监督检测的精度。

首先，以二相水平集为例，可将影像分为两部分，如下式所示：

其中，Ω_in，Ω_out和c₁，c₂分别表示轮廓曲线C内部和分外区域及其均值，μ为调节轮廓曲线C的常数。上式可用水平集φ表示：

其中，H表示Heaviside函数，如下式：

水平集φ可定按下式定义：

多相水平集可将影像分为多个类别，以二相水平集φ₁和φ₂为例，可按下式将影像分为四个类别，如图4的(a)、(b)所示。

其中，c₁₁表示同时位于水平集φ₁和φ₂内部的像素均值，c₁₀表示同时位于φ₁内部和φ₂外部像素的均值，c₀₁表示同时位于φ₁外部和φ₂内部像素的均值，c₀₀同时位于水平集φ₁和φ₂外部的像素均值。根据e＝log₂(K+1)确定包含变化类和未变化类的所有类别数K+1与水平集相数e的对应关系，最终利用拉格朗日算法求解能量最小解，得到最优的水平集，最终实现多类别变化的非监督检测。

本发明分别提取两时期影像的光谱变化特征和光谱梯度差异特征，按波段进行叠加，并利用主成分分析方法对其进行主成分分析，进行相关性抑制和异质性增强，最终选择前两个主成分分量作为最终的变化特征，计算变化强度和变化方向。然后，在极坐标框架下利用Bayesian最小错误准则计算变化强度阈值，将像素分为未变化和变化两部分，利用FCM聚类方法对变化像素的变化方向设置阈值，区分变化类别。在此基础上，在极坐标框架下对初始变化检测结果进行分析，对所有变化类别，依次生成仅包含单一变化的像素极，再利用Bayesian最小错误准则获取针对各变化类别的最优变化强度阈值，从而生成基于动态阈值的变化检测结果。最后，利用多相水平集算法对基于动态阈值的多类别变化检测结果进行后处理，通过集成空间信息同时去除虚检和漏检像素。本发明从多特征融合、动态阈值策略和集成空间信息的多相水平集三个方面克服现在多类别变化非监督变化检测的不足，可有效提高多类别变化检测的精度。

以下结合附图和实例详细说明本发明的技术方案。

实例用两个时期获取的Landsat-5TM卫星数据进行实验，实验区数据大小为412×300像素，分辨率为30米，获取时间分别是1995年9月和1996年7月，采用除热红外波段的6个波段光谱信息进行变化检测，如图5的(a)、(b)和(c)所示。实验区包括三种变化，分别为由于采矿引起的变化、由于森林火灾引起的变化和湖泊水位变化引起的变化。

一种利用遥感影像多特征的多类别变化动态阈值检测方法，包括以下步骤：

步骤1：对两时期多波段遥感影像进行几何校正、配准、辐射校正；

步骤2：分别计算两时期影像的光谱变化特征(6波段)及光谱梯度差异特征(5波段)，并通过波段叠加生成多特征(11波段)。在此基础上，利用主成分分析方法对叠加后的特征进行主成分分析，选取前两个主成分分量作为最终的变化特征，在此基础上，计算变化强度和变化方向，进行多类别变化检测。

步骤3：在极坐标框架下，对变化强度和变化方向进行空间分析，图6的(a)、(b)和(c)分别为光谱变化特征、光谱梯度差异特征和主成分分量特征在极坐标空间中的分布，可以看出本发明采用的主成分分量特征最大程度的区分不同类别的变化。假设变化强度为高斯混合分布，利用Bayesian最小错误准则计算变化强度阈值，与此同时，通过视觉分析可得出变化类别K＝3，利用FCM对变化像素的变化方向进行聚类，区分不同类别的变化，最终阈值设置为：S₁＝{ρ,α:ρ≥54，0°≤α≤59°∪278°≤α≤360°}，S₂＝{ρ,α:ρ≥54，73°≤α<145°}和S₃＝{ρ,α:ρ≥54，145°≤α<278°}，分别对应采矿变化、火灾变化和湖泊变化。

步骤4：分析步骤3所得的初步变化检测结果，针对各变化类别，构建新的仅包含一类变化的像素集，利用步骤3中提取的Bayesian最小错误准则对各像素集设置最优的变化强度阈值，最终得到的动态阈值分别为：S₁＝{ρ,α:ρ≥58，0°≤α≤59°∪278°≤α≤360°}，S₂＝{ρ,α:ρ≥46，73°≤α<145°}和S₃＝{ρ,α:ρ≥63，145°≤α<278°}，生成基于动态阈值的多类别变化检测结果。

步骤5：通过步骤4中的分析，本实例中共有K+1＝4种类别，即未变化、采矿变化、火灾变化和湖泊变化，因此采用二相水平集算法，对步骤4中生成的基于动态阈值的多类别变化检测检测进行后处理，设置轮廓长度控制参数μ＝0.12。最终通过二相水平集的演化，集成空间信息同时去除虚检和漏检测像素，进一步提高多类别变化非监督检测的精度。

图7的(a)、(b)和(c)分别为基于光谱变化特征的固定阈值、光谱梯度差异特征的固定阈值和本发明方法得到的多类别变化非监督检测结果，通过与参考数据的比较可以看出，前两种传统方法的结果中存在较多的虚检和漏检错误，而本发明提出的方法可以生成与参考数据更相似的检测结果。表1、2和3分别为上述三种方法得到变化检测结果定量评价精度。从表中可以看出：本发明方法对在总体提高了各类别的生产者精度和用户精度，总分类精度和Kappa系数都有很大提高，极大提高了遥感影像多类别变化非监督检测的精度。

表1：基于光谱变化特征的固定阈值检测结果精度

表2：基于光谱梯度差异特征的固定阈值检测结果精度

表3：本发明方法得到的变化检测结果精度

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (5)

1.一种利用遥感影像多特征的多类别变化动态阈值检测方法，其特征在于，包括如下步骤：

步骤1，对前、后两时期的遥感影像进行预处理，包括几何校正、配准、辐射校正；

步骤2，计算经步骤1预处理后的前、后两时期遥感影像的光谱变化特征和光谱梯度差异特征，将光谱变化特征和光谱梯度差异特征按波段进行叠加，利用主成分分析方法对叠加后的特征进行主成分分析，选取前两个主成分分量作为变化特征，根据变化特征计算变化强度和变化方向；

步骤3，在极坐标框架下，对变化强度和变化方向进行空间分析，针对变化强度计算变化强度阈值，根据变化强度阈值检测未变化和变化类别；针对变化类别的变化方向设置阈值，根据阈值区分不同的变化类别，生成初步多类别变化检测结果；

步骤4，分析初步多类别变化检测结果，针对各变化类别，构建仅包含当前变化类别的像素集，并设置动态的变化强度阈值，重新判定各变化类别中的变化像素与未变化像素，生成基于动态阈值的多类别变化检测结果。

2.根据权利要求1所述利用遥感影像多特征的多类别变化动态阈值检测方法，其特征在于，该方法还包括步骤5，利用多相水平集算法，对步骤4生成的基于动态阈值的多类别变化检测结果进行处理，通过集成空间信息同时去除虚检测和漏检测像素。

3.根据权利要求1所述利用遥感影像多特征的多类别变化动态阈值检测方法，其特征在于，步骤2所述光谱变化特征和光谱梯度差异特征计算公式分别为：

X_D＝X₂-X₁

G＝G₂-G₁

其中，X_D表示光谱变化特征，X₁、X₂分别表示经预处理后的前、后时期遥感影像，G表示光谱梯度差异特征，G₁、G₂分别表示经处理后的前、后时期遥感影像光谱梯度。

4.根据权利要求1所述利用遥感影像多特征的多类别变化动态阈值检测方法，其特征在于，步骤2所述变化强度和变化方向计算公式分别为：

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>&amp;rho;</mi> <mo>=</mo> <msqrt> <mrow> <msubsup> <mi>P</mi> <mn>1</mn> <mn>2</mn> </msubsup> <mo>+</mo> <msubsup> <mi>P</mi> <mn>2</mn> <mn>2</mn> </msubsup> </mrow> </msqrt> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>&amp;alpha;</mi> <mo>=</mo> <mi>arctan</mi> <mrow> <mo>(</mo> <mfrac> <msub> <mi>P</mi> <mn>1</mn> </msub> <msub> <mi>P</mi> <mn>2</mn> </msub> </mfrac> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced>

其中，ρ、α分别表示变化强度、变化方向，P₁、P₂表示经主成分分析后选取前两个主成分分量。

5.根据权利要求1所述利用遥感影像多特征的多类别变化动态阈值检测方法，其特征在于，步骤4所述针对各变化类别，构建仅包含当前变化类别的像素集具体为：

获取各变化类别在极坐标空间中的分布，构建仅包含单一变化类别的像素集，针对变化类别k，找出该类别中变化方向的最小值T_α,k和最大值T_α,k+1，将极坐标空间中变化方向值在(T_α,k,T_α,k+1)内的变化像素集S_k和未变化像素集合并，生成仅包含变化类别k的像素集S：

<mrow> <mi>S</mi> <mo>=</mo> <msub> <mi>S</mi> <mi>k</mi> </msub> <mo>&amp;cup;</mo> <msubsup> <mi>S</mi> <mi>k</mi> <mo>*</mo> </msubsup> <mo>=</mo> <mo>{</mo> <msub> <mi>T</mi> <mrow> <mi>&amp;alpha;</mi> <mo>,</mo> <mi>k</mi> </mrow> </msub> <mo>&lt;</mo> <mi>&amp;alpha;</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>&lt;</mo> <msub> <mi>T</mi> <mrow> <mi>&amp;alpha;</mi> <mo>,</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>}</mo> </mrow>

其中，α(i,j)表示极坐标空间中变化像素或未变化像素的变化方向值。