CN104156955A - 一种高分辨率遥感影像的变化检测方法 - Google Patents

Abstract

本发明涉及一种高分辨率遥感影像的变化检测方法，所述方法包括如下步骤：(1)利用Mean-shift影像分割算法，对不同时相的影像进行多尺度分割，形成不同尺度的地理对象；(2)多尺度地理对象的特征提取与特征向量构建；(3)基于变化矢量法(change vector analysis，CVA)对不同时相多尺度地理对象特征向量进行变化检测，获得最终的变化检测的结果。

Description

一种高分辨率遥感影像的变化检测方法

技术领域

本发明涉及一种高分辨率遥感影像的变化检测方法。

背景技术

在像素水平上的变化检测算法应用到高空间分辨率的影像中时很难检测非均质的面状区域的变化；单一尺度的面相对象的变化检测算法存在过分割或不完全分割的情况，也会影响变化检测的准确性。

发明内容

在以上基础上，本发明提出了采用多尺度的面相对象的变化检测算法，可全自动提取高分辨率影像变化区域，提高变化检测的准确性。

本发明提供一种高分辨率遥感影像的变化检测方法，所述方法包括如下步骤：(1)利用Mean-shift影像分割算法，对不同时相的影像进行多尺度分割，形成不同尺度的地理对象；(2)多尺度地理对象的特征提取与特征向量构建；(3)基于变化矢量法(change vector analysis，CVA)对不同时相多尺度地理对象特征向量进行变化检测，获得最终的变化检测的结果。

附图说明

通过参照附图更详细地描述本发明的示例性实施例，本发明的以上和其它方面及优点将变得更加易于清楚，在附图中：

图1为本发明的一种高分辨率遥感影像的变化检测方法的流程图。

具体实施方式

在下文中，现在将参照附图更充分地描述本发明，在附图中示出了各种实施例。然而，本发明可以以许多不同的形式来实施，且不应该解释为局限于在此阐述的实施例。相反，提供这些实施例使得本公开将是彻底和完全的，并将本发明的范围充分地传达给本领域技术人员。

在下文中，将参照附图更详细地描述本发明的示例性实施例。

对于已经配准好的两个不同时期同一个区域的高分辨率遥感影像，采用多尺度分割的变化检测共分为三个步骤(具体见图1)。(1)利用Mean-shift影像分割算法，对不同时相的影像进行多尺度分割，形成不同尺度的地理对象；(2)多尺度地理对象的特征提取与特征向量构建；(3)基于变化矢量法(change vectoranalysis，CVA)对不同时相多尺度地理对象特征向量进行变化检测，获得最终的变化检测的结果。

步骤1)：多尺度分割

影像分割的目的是将具有相同特征的像素聚集在一起，形成能反映特定对象特征的地理对象。但是对于高分辨率影像来说，当选择不同的分割尺度时，分割的结果不一样。由于地物的结构、纹理、空间特征以及相邻地物的关系等信息十分的丰富，采用单一尺度的分割，必然会出现过分割或者分割不完全等情况，很难完整的反应对象的特征。而通过多尺度的分割，形成了多个地理对象层，可以从不同的角度反映遥感数据的客观表现形态，为描述地理对象的多尺度特征提供了很好的基础。本发明采用Mean-Shift分割算法结合分割后区域合并算法，自下而上，从单像素开始，建立多尺度的地理对象。具体的步骤包括：

(1)采用Mean-Shift分割算法，形成初始分割结果，并对分割结果进行标记，统计邻接关系，得到区域邻接图，获得最初始的地理对象。

(2)对(1)中形成的地理对象，根据区域相似性准则对相邻的地理对象进行合并，形成高一层次的地理对象。区域相似性准则考虑区域的灰度相似性、公共边缘强度以及区域的面积作为相似性度量的依据。具体的公式如下：

$Z(R_i,R_j)=\frac{1}{D_{R_i{R}_j}{E}_{R_i{R}_j}{A}_{R_i{R}_j}+1}---\left(1\right)$

其中D_RiRj表示面积加权灰度均值距离，具体见公式(2)，E_RiRj表示相邻地理对象边缘梯度差值，A_RiRj表示地理对象的邻接关系。

$D_{R_i,R_j}=\frac{N\left(R_i\right)+N\left(R_j\right)}{N\left(R_i\right)+N\left(R_j\right)}\left(\right||{\mathrm{\μ}}_{R_i}-{\mathrm{\μ}}_{R_j}|\left|\right)---\left(2\right)$

$E_{R_i,R_j}=\left|\right|\mathrm{MG}\left(R_i\right)-\mathrm{MG}\left(R_i\right)\left|\right|---\left(3\right)$

其中，N(R_i)和N(R_j)分别表示邻接区域的像素数量，μ表示邻接区域的灰度均值，MG表示区域边缘处的梯度均值。

采用基于最大区域熵的方法对图像区域进行合并，合并停止准则采用图像区域最大熵值，因为从信息论的角度看，灰度变化小的地方其信息量少，局部熵值小；灰度变化大的地方信息量多，局部熵值大。信息熵的计算公式见(5).p_kl表示在地里对象R_i中，地理位置(k，1)上的灰度出现的频率。

$H\left(R_i\right)=-{\mathrm{\Σ}}_{(k,l)\∈R_i}{p}_{\mathrm{kl}}\ln \left(p_{\mathrm{kl}}\right)---\left(5\right)$

$p_{\mathrm{kl}}=\frac{f(k,l)}{{\mathrm{\Σ}}_{(k,l)\∈R_i}f(k,l)}---\left(6\right)$

相邻两区域合并的规则如下：即先计算某一区域的信息熵，然后将其与相邻区域合并后，在重新计算区域熵，如果合并后的区域熵小于合并前的区域熵，则这两个区域是相似的，可以合并，否则这两个区域不能合并。

(3)重复步骤(2)直到对所有的区域遍历，

(4)重复步骤(2)和(3)直到完成整个多尺度地理对象构建。

步骤2)：多尺度地理对象特征提取

通过Mean-Shift分割算法结合分割后区域合并的策略，建立了多个尺度上，地理对象的分割结果。对每一个尺度上的地理对象，需要提取该尺度上地理对象的特征参数来反映该地理对象的特征。本发明选择地理对象中所有像素的灰度均值向量作为特征向量衡量变化前后的差异，因为该特征计算最方便，且反应的变化特征最明显。对每一个像素来说，构建的多尺度地理对象特征向量可以描述为

M(k，l)＝{x(k，l)，μ¹(k，l)，...，μ^L(k，l)}  (7)

其中，

${\mathrm{\μ}}^l=\frac{1}{N_i}{\mathrm{\Σ}}_{(k,l)\∈R_i}x(k,l)---\left(8\right)$

其中，i，j表示影像上的像素点坐标，R_i表示第i个地理对象的所有像素集合，N_i表示第i个地理对象的像素个数，μ表示地理对象的灰度均值，L表示尺度的层数。

采用多尺度的分割算法其目的不是找出地理对象的最佳分割尺度；而是建立了每一个像素不同层次的空间邻域信息。L层数的选择取决于影像的几何分辨率以及需要检测的地物对象大小。

步骤3)：尺度地理对象变化矢量检测

由公式(7)中构建的地理对象的特征向量，实际上代表了每一个像素在不同分割尺度上的灰度均值特征，反映了每个像素的多尺度空间邻域信息。通过计算每一个像素所在地理对象中的特征向量M(i，j)，然后利用变化矢量分析法(change vector analysis，CVA)方法可以直接计算两个不同时期的特征向量的差异，构建面向地理对象的差异影像。这种在像素水平上计算的差异值，实际上反应了该像素所在的空间区域的信息，具体的计算如公式(9)所示。

M_D(i，j)＝|M_t1(i，j)-M_t2(i，j)|  (9)对于形成的差异影像M_D，根据阈值可以判断变化与非变化的区域。阈值的选择方法本发明采用了基于最大期望(EM)的自动阈值选择算法。该算法是对不完整数据问题进行最大似然估计的一种常用算法，它无需任何外来数据和先验知识，仅从观测数据本身得到参数的估计值。

由变化矢量分析法得到的差异影像M_D可认为由未变化的像元类ω_n和变化的像元类ω_c所构成。假设未变化的像元类ω_n和变化的像元类ω_c的概率密度函数p(x/ω_n)，p(x/ω_c)服从高斯密度函数分布，如(10)：

$\begin{array}{cc}p\left(x\right|{\mathrm{\ω}}_i)=\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_i^2}}\exp \{-\frac{{(x-m_i)}^2}{2{\mathrm{\σ}}_i^2}\}& i\∈\{n,c\}---\left(10\right)\end{array}$

未变化的像元类ω_n和变化的像元类ω_c的均值和方差分别用m_n、m_c、来表示，则差异影像的总密度函数p(x)是可表示为：

p(x)＝p(x|ω_n)p(ω_n)+p(x|ω_c)p(ω_c)  (11)

在上述假设条件下，求解差异影像的变化阈值的过程可转化为估算ω_n、ω_c的均值与方差。对m_n、m_c、的估算可采用EM算法实现。EM算法通过循环迭代，每次迭代由求期望值和期望最大化两个步骤组成。前者根据待估计参数的当前值，从观测数据中直接估计概率密度的期望值，后者通过最大化这一期望来更新参数的估计量，这两步在整个迭代过程中依次交替进行，直至迭代过程收敛。采用EM算法来估计未变化像元类别ω_n统计分布参数的计算公式如下：

$p^{t+1}\left({\mathrm{\ω}}_n\right)=\frac{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t(X(i,j)/{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}}{I*J}---\left(12\right)$

$m_n^{\text{t+1}}\left({\mathrm{\ω}}_n\right)=\frac{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t(X(i,j)/{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t(X(i,j)/{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}}---\left(13\right)$

${\left({\mathrm{\σ}}_n^2\right)}^{t+1}\left({\mathrm{\ω}}_n\right)=\frac{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t(X(i,j)/{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}[X(i,j)-m_n^t{]}^2}{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t(X(i,j)/{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}}---\left(14\right)$

上面三式估计的分别是先验概率、均值和标准偏差，式中的t和t+1分别代表了当前和下一次迭代所用的估计值，i，j分别代表了差值影像的行数和列数，X(i，j)表示差异影像中第i行j列的像元值，条件概率p(X(i，j)|ω_n)的计算见(10)，全概率p(X(i，j))的值由式(11)给出。变化类ω_c的先验概率p(ω_c)、均值m_c和标准差σ_c的估计方法同上。当相邻两次迭代计算的p(ω_i)、m_i和σ_i的值小于给定的阈值ε(ε＝10^-8)时迭代终止。

在利用EM算法估算了m_n、m_c、后，根据贝叶斯最小错误率理论，最佳变化阈值应当是满足(15)式时的T值。

p(ω_n)p₁(T|ω_n)＝p(ω_c)p(T|ω_c)  (15)

根据ω_n、ω_c的整体分布特性，将结合(10)式，将m_n、m_c、带入到(15)式可得到变化阈值T的求解公式(16)。

$({\mathrm{\σ}}_n^2-{\mathrm{\σ}}_c^2)T^2+2(m_n{\mathrm{\σ}}_c^2-m_c{\mathrm{\σ}}_n^2)T+m_c^2{\mathrm{\σ}}_n^2-m_c^2{\mathrm{\σ}}_n^2-2{\mathrm{\σ}}_c^2{\mathrm{\σ}}_n^2\mathrm{In}\left[\frac{{\mathrm{\σ}}_{n}p\left({\mathrm{\ω}}_n\right)}{{\mathrm{\σ}}_{c}p\left({\mathrm{\ω}}_c\right)}\right]=0$ $\left(16\right)$

以上所述仅为本发明的实施例而已，并不用于限制本发明。本发明可以有各种合适的更改和变化。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种高分辨率遥感影像的变化检测方法，其特征在于：

所述方法包括如下步骤：(1)利用Mean-shift影像分割算法，对不同时相的影像进行多尺度分割，形成不同尺度的地理对象；(2)多尺度地理对象的特征提取与特征向量构建；(3)基于变化矢量法对不同时相多尺度地理对象特征向量进行变化检测，获得最终的变化检测的结果。

2.如权利要求1所述的一种高分辨率遥感影像的变化检测方法，其特征在于：

所述步骤(1)的具体步骤包括：

(a)：采用Mean-Shift分割算法，形成初始分割结果，并对分割结果进行标记，统计邻接关系，得到区域邻接图，获得最初始的地理对象；

(b)：对(a)中形成的地理对象，根据区域相似性准则对相邻的地理对象进行合并，形成高一层次的地理对象；区域相似性准则考虑区域的灰度相似性、公共边缘强度以及区域的面积作为相似性度量的依据，具体的公式如下：

$Z(R_i,R_j)=\frac{1}{D_{R_i{R}_j}{E}_{R_i{R}_j}{A}_{R_i{R}_j}+1}---\left(1\right)$

其中D_RiRj表示面积加权灰度均值距离，具体见公式(2)，E_RiRj表示相邻地理对象边缘梯度差值，A_RiRj表示地理对象的邻接关系；

$D_{R_i,R_j}=\frac{N\left(R_i\right)+N\left(R_j\right)}{N\left(R_i\right)+N\left(R_j\right)}\left(\right||{\mathrm{\μ}}_{R_i}-{\mathrm{\μ}}_{R_j}|\left|\right)---\left(2\right)$

$E_{R_i,R_j}=\left|\right|\mathrm{MG}\left(R_i\right)-\mathrm{MG}\left(R_i\right)\left|\right|---\left(3\right)$

其中，N(R_i)和N(R_j)分别表示邻接区域的像素数量，μ表示邻接区域的灰度均值，MG表示区域边缘处的梯度均值；

采用基于最大区域熵的方法对图像区域进行合并，合并停止准则采用图像区域最大熵值；信息熵的计算公式见(5).p_kl表示在地里对象R_i中，地理位置(k，1)上的灰度出现的频率；

$H\left(R_i\right)=-{\mathrm{\Σ}}_{(k,l)\∈R_i}{p}_{\mathrm{kl}}\ln \left(p_{\mathrm{kl}}\right)---\left(5\right)$

$p_{\mathrm{kl}}=\frac{f(k,l)}{{\mathrm{\Σ}}_{(k,l)\∈R_i}f(k,l)}---\left(6\right)$

相邻两区域合并的规则如下：即先计算某一区域的信息熵，然后将其与相邻区域合并后，在重新计算区域熵，如果合并后的区域熵小于合并前的区域熵，则这两个区域是相似的，可以合并，否则这两个区域不能合并；

(c)：重复步骤(b)直到对所有的区域遍历，

(d)：重复步骤(b)和(c)直到完成整个多尺度地理对象构建。

3.如权利要求1所述的一种高分辨率遥感影像的变化检测方法，其特征在于：

所述步骤(2)的具体步骤包括：

选择地理对象中所有像素的灰度均值向量作为特征向量衡量变化前后的差异；对每一个像素来说，构建的多尺度地理对象特征向量可以描述为

M(k，l)＝{x(k，l)，μ¹(k，l)，...，μ^L(k，l)}  (7)

其中，

${\mathrm{\μ}}^l=\frac{1}{N_i}{\mathrm{\Σ}}_{(k,l)\∈R_i}x(k,l)---\left(8\right)$

其中，i，j表示影像上的像素点坐标，R_i表示第i个地理对象的所有像素集合，N_i表示第i个地理对象的像素个数，μ表示地理对象的灰度均值，L表示尺度的层数。

4.如权利要求1-3任意一项所述的一种高分辨率遥感影像的变化检测方法，其特征在于：

所述步骤(3)的具体步骤包括：

通过计算每一个像素所在地理对象中的特征向量M(i，j)，然后利用变化矢量分析法方法可以直接计算两个不同时期的特征向量的差异，构建面向地理对象的差异影像；具体的计算如公式(9)所示；

M_D(i，j)＝|M_t1(i，j)-M_t2(i，j)|  (9)

对于形成的差异影像M_D，根据阈值可以判断变化与非变化的区域；阈值的选择方法采用基于最大期望(EM)的自动阈值选择算法；

由变化矢量分析法得到的差异影像M_D可认为由未变化的像元类ω_n和变化的像元类ω_c所构成；假设未变化的像元类ω_n和变化的像元类ω_c的概率密度函数p(x/ω_n)，p(x/ω_c)服从高斯密度函数分布，如(10)：

$\begin{array}{cc}p\left(x\right|{\mathrm{\ω}}_i)=\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_i^2}}\exp \{-\frac{{(x-m_i)}^2}{2{\mathrm{\σ}}_i^2}\}& i\∈\{n,c\}---\left(10\right)\end{array}$

未变化的像元类ω_n和变化的像元类ω_c的均值和方差分别用m_n、m_c、来表示，则差异影像的总密度函数p(x)是可表示为：

p(x)＝p(x|ω_n)p(ω_n)+p(x|ω_c)p(ω_c)  (11)

在上述假设条件下，求解差异影像的变化阈值的过程可转化为估算ω_n、ω_c的均值与方差；对m_n、m_c、的估算可采用EM算法实现；EM算法通过循环迭代，每次迭代由求期望值和期望最大化两个步骤组成；前者根据待估计参数的当前值，从观测数据中直接估计概率密度的期望值，后者通过最大化这一期望来更新参数的估计量，这两步在整个迭代过程中依次交替进行，直至迭代过程收敛；采用EM算法来估计未变化像元类别ω_n统计分布参数的计算公式如下：

$p^{t+1}\left({\mathrm{\ω}}_n\right)=\frac{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t(X(i,j)/{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}}{I*J}---\left(12\right)$

$m_n^{\text{t+1}}\left({\mathrm{\ω}}_n\right)=\frac{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t(X(i,j)/{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}X(i,j)}{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t(X(i,j)/{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}}---\left(13\right)$

${\left({\mathrm{\σ}}_n^2\right)}^{t+1}\left({\mathrm{\ω}}_n\right)=\frac{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t(X(i,j)/{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}[X(i,j)-m_n^t{]}^2}{\underset{X(i,j)\∈M_D}{\mathrm{\Σ}}\frac{p^t\left({\mathrm{\ω}}_n\right)p^t(X(i,j)/{\mathrm{\ω}}_n)}{p^t\left(X(i,j)\right)}}---\left(14\right)$

上面三式估计的分别是先验概率、均值和标准偏差，式中的t和t+1分别代表了当前和下一次迭代所用的估计值，i，j分别代表了差值影像的行数和列数，X(i，j)表示差异影像中第i行j列的像元值，条件概率p(X(i，j)|ω_n)的计算见(10)，全概率p(X(i，j))的值由式(11)给出；变化类ω_c的先验概率p(ω_c)、均值m_c和标准差σ_c的估计方法同上；当相邻两次迭代计算的p(ω_i)、m_i和σ_i的值小于给定的阈值ε(ε＝10^-8)时迭代终止；

在利用EM算法估算了m_n、m_c、后，根据贝叶斯最小错误率理论，最佳变化阈值应当是满足(15)式时的T值；

p(ω_n)p₁(T|ω_n)＝p(ω_c)p(T|ω_c)  (15)

根据ω_n、ω_c的整体分布特性，将结合(10)式，将m_n、m_c、带入到(15)式可得到变化阈值T的求解公式(16)；

$({\mathrm{\σ}}_n^2-{\mathrm{\σ}}_c^2)T^2+2(m_n{\mathrm{\σ}}_c^2-m_c{\mathrm{\σ}}_n^2)T+m_c^2{\mathrm{\σ}}_n^2-m_c^2{\mathrm{\σ}}_n^2-2{\mathrm{\σ}}_c^2{\mathrm{\σ}}_n^2\mathrm{In}\left[\frac{{\mathrm{\σ}}_{n}p\left({\mathrm{\ω}}_n\right)}{{\mathrm{\σ}}_{c}p\left({\mathrm{\ω}}_c\right)}\right]=0$ $\left(16\right)$

。