CN110276746B - 一种鲁棒性遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种鲁棒性遥感图像变化检测方法，包括以下步骤：采集获取目标地区两个时相的图像信息；对两个时相的遥感图像进行预处理，获取图像特征；根据Relief算法对图像特征进行选择和筛选；运用PCA算法进行特征降维，得到最终优选的特征变化矢量；利用NR‑Kmeans方法进行分类，得到变化检测结果。本发明改进了传统的kmeans聚类算法，检测精度更高，准确率高。

Description

一种鲁棒性遥感图像变化检测方法

技术领域

本发明属于图像变化技术领域，尤其涉及一种鲁棒性遥感图像变化检测方法。

背景技术

遥感，即远距离感知。遥感一词最早于1960年出现在美国的一项军事科研计划之中，在1962年密执安大学等单位举行的第一届环境遥感学术研讨会上被正式采用。如今，遥感技术应用于越来越多的领域，比如：森林或植被的动态变化监测、对自然灾害的灾后分析及评估、对土地利用的变化分析、对农田进行监控、对城镇变化实时监测、分析农作物生长状况、对军事战略目标，比如机场、道路，进行动态监视等领域，极大地促进了经济和社会的发展。

遥感图像变化检测属于遥感图像处理领域，用于分析处理同一地点不同时期的遥感图像而获得变化信息。国内外学者已经对遥感图像变化检测问题进行了大量的研究，提出了各种各样的变化检测方法。其中基于特征的变化检测方法是利用图像中的光谱特征、纹理特征、空间结构特性，对遥感图像进行变化检测。基于特征的变化检测方法考虑了图像丰富的特征信息，具有较好的噪声鲁棒性和抗光照、辐射干扰能力。

由于遥感图像受光照、大气环境等影响，而采集的图像特征也受其影响。不仅如此，由于噪声干扰，传统的聚类算法受稀疏数据的影响分类的结果不够精确。选择具有代表性的特征并精确分类成为遥感图像变化检测技术的关键。

发明内容

发明目的：针对以上问题，本发明提出一种鲁棒性遥感图像变化检测方法。该方法是一种基于特征的遥感图影变化检测方法，根据Relief算法对图像特征进行选择和筛选，然后利用PCA算法降维，最后使用NR-Kmeans算法，排除噪声干扰点对检测结果的影响，最后输出结果。本发明排除了噪声干扰点，检测更加精确。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种鲁棒性遥感图像变化检测方法，包括如下步骤：

(1)采集获取目标地区T₁和T₂时刻的遥感图像，所述图像满足分辨率需求并且包含清晰的纹理特征信息；

(2)采用ENVI软件实现图像数据预处理，消除拍摄误差；采用eCogntion软件平台对T₁和T₂时刻的预处理后的图像进行多尺度分割，将两幅图像分别划分成N个图像块，即分别得到N个对象；分别获取T₁和T₂时刻的图像的光谱特征和纹理特征，构建特征向量，得到变化矢量集合M_d，M_d＝[m₁,m₂,...m_N]，其中，m_i表示第i个对象所对应的变化矢量；将M_d中元素按照T₁或T₂时刻图像中N个对象的排列方式排列，构成一幅差值图；

(3)根据Relief特征选择算法，获得变化矢量集合M_d的特征子空间M_l中每个对象的特征的权重，选取e个权重最大的变化特征，e取值小于对象包含的特征总数，进而分别在N个对象中选取所述e个变化特征，结合变化矢量集合M_d，得到矩阵X；

(4)利用PCA算法对矩阵X进行特征降维，得到变换之后的变化特征矩阵Y，Y＝[y₁，y₂，...,y_N]；

(5)使用NR-Kmeans算法对变化特征矩阵Y的元素y_i进行分类，从变化特征矩阵Y中剔除被噪声干扰的稀疏数据，得到密集点集合Y′；

(6)将密集点集合Y′分为未变化的类和变化的类，未变化的类元素对应的差值图中的对象的灰度值转化为255，变化的类中元素对应的差值图中的对象的灰度值转化为0，根据对象的灰度值输出变化检测结果图，图中标出了图像变化的区域。

进一步，步骤(2)中，采用ENVI软件实现图像数据预处理，消除拍摄误差；采用eCogntion软件平台对T₁和T₂时刻的预处理后的图像进行多尺度分割，将两幅图像分别划分成N个图像块，即分别得到N个对象；分别获取T₁和T₂时刻的图像的光谱特征和纹理特征，构建特征向量，得到变化矢量集合，步骤如下：

(2.1)采用ENVI软件分别对T₁和T₂时刻的图像进行几何校正、图像配准、大气校正操作，实现图像数据预处理；

(2.2)采用eCogntion软件平台对T₁和T₂时刻的预处理后的图像进行多尺度分割，将两幅图像分别划分成N个图像块，即分别得到N个对象；

(2.3)求取预处理后的图像的Mean-Std特征，作为图像的光谱特征；采用灰度共生矩阵提取预处理后的图像的纹理特征；

(2.4)根据预处理后的两幅图像提取的光谱特征和纹理特征，构建特征向量，获取变化矢量集合M_d。

进一步，步骤(2.3)所述求取图像的Mean-Std特征，作为图像的光谱特征；方法如下：

对图像的每个对象，按照下列公式求取光谱特征中的均值Mean和标准差Std：

式中，Mean为对象中像素点灰度的均值，std为对象中像素点灰度的标准差，A表示对象中像素点数，c_i代表对象中第i个像素点的灰度大小。

进一步，步骤(2.3)所述采用灰度共生矩阵提取图像的纹理特征；方法如下：

(2.3.1)将图像中任意一像素点(τ₁，τ₂)及另一像素点(τ₁+Δ₁，τ₂+Δ₂)构成点对，其中Δ₁，Δ₂为整数；设像素点(τ₁，τ₂)的灰度值为f1，(τ₁+Δ₁，τ₂+Δ₂)的灰度值为f2，则该点对的灰度值为(f1，f2)，设图像的最大灰度级为L，则f1与f2的组合共有L*L组；

(2.3.2)统计图像中每一组不同的(f1，f2)值出现的次数，然后排列成一个矩阵，其中，位于矩阵位置(L，L)上的矩阵元素值就是两个灰度值均为L的点对出现的次数；

(2.3.3)根据(f1，f2)出现的总次数将每一组不同的(f1，f2)值出现的次数归一化为出现的概率g(f1，f2)，归一化之后的的方阵称为灰度共生矩阵；

(2.3.4)灰度共生矩阵提供了14种统计量作为纹理特征，本发明选取6种统计量作为纹理特征，所述6种统计量为对比度、熵、能量、均等性、不相似性、相关性；对比度反映了图像的清晰度和纹理沟纹深浅的程度；熵表示了图像中纹理的非均匀程度或复杂程度；能量反映了图像灰度分布均匀程度和纹理粗细度；均等性反映图像纹理的同质性，度量图像纹理局部变化的多少；相关性反应了图像纹理的一致性；这6种统计量都从各个方面表示了图像的特性，故选择这6种统计量作为纹理特征；

(2.3.5)设灰度共生矩阵的位置(i，j)上的矩阵元素值为g(i，j)，L为图像的最大灰度级，则所述纹理特征表示如下：

式中，Con表示对比度，Ent表示熵，Ene表示能量，Hom表示均等性，Dis表示不相似性，Cor表示相关性；

进一步，步骤(2.4)所述根据预处理后的两幅图像提取的光谱特征和纹理特征，获取变化矢量集合M_d；方法如下：

(2.4.1)根据步骤(2.3)得到的T₁和T₂时刻图像的光谱特征和纹理特征，分别构建特征向量，对两幅图像的相对应的特征向量计算差值；两幅图像相对应的第i个对象的第j个特征向量的差值M_d(i，j)，表示如下：

M_d(i，j)＝M_d，1(i，j)-M_d，2(i，j)

式中，M_d，1(i，j)为T₁时刻图像第i个对象的第j个特征向量，M_d，2(i，j)为T₂时刻图像第i个对象的第j个特征向量；i≤N，N为图像中对象的总数；j≤S_feature，S_feature为对象包含的特征总数；S_feature＝S_band×8，S_band为对象总波段数；

根据差值M_d(i，j)，得到第i个对象所对应的变化矢量m_i，表示如下：

m_i＝(M_d(i，1)，M_d(i，2)，...，M_d(i，S_feature))

(2.4.2)重复步骤(2.4.1)，计算两幅图像中相对应的每一个对象的特征向量的差值，得到每个对象所对应的变化矢量，构成变化矢量集合M_d，表示如下：

M_d＝[m₁，m₂，…m_N]

将M_d中元素按照T₁和T₂时刻图像中N个对象的排列方式排列，构成一幅差值图。

进一步，步骤(3)所述根据Relief特征选择算法，获得变化矢量集合M_d的特征子空间M_i中每个对象的特征的权重，选取e个权重最大的变化特征，e取值小于对象包含的特征总数，进而分别在N个对象中选取所述e个变化特征，结合变化矢量集合M_d，得到矩阵X；方法如下：

(3.1)变化矢量集合M_d也称为特征空间，提取特征空间M_d＝[m₁，m₂，m₃，…m_N]的子空间M_l＝[m₁,m₂,m₃,…m_n]作为训练集，其中n＜N；使用人工解译将特征子空间M_l分为M_l1与M_l2两个数集，其中M_l1包含l1个元素，M_L2包含l2个元素，l1+l2＝n，M_l1∈M_l，M_l2∈M_l；

(3.2)计算M_l中所有元素两两之间变化矢量的欧式距离，组成距离矩阵，即

其中，d表示距离矩阵，

表示M_l中元素两两之间变化矢量的欧式距离；

(3.3)在特征子空间M_l中，对于对象i的变化矢量m_i，i＝1，2，3，...，n，若m_i∈M_l1，选取M_l1内与其欧式距离最小的s个矢量，并选取M_l2内与其欧氏距离最大的s个矢量；若m_i∈M_l2，选取M_l2内与其欧式距离最小的s个矢量，并选取M_l1内与其欧氏距离最大的s个矢量；即对m_i求取近邻矩阵Q_i＝[q₁，q₂，...，q_s，...，q_2s]，其中，[q₁，q₂，...，q_s]为所述欧式距离最小的s个矢量，[q_s+1,q_s+2，...，q_2s]为所述欧氏距离最大的s个矢量，s＝min(l1，l2)；

(3.4)根据猜错近邻相关统计量和猜中近邻相关统计量，计算对象i的每个特征对应的权重，即

式中，

是特征f_j对应的权重，diffr_j是特征f_j的猜错近邻相关统计量，diffw_j是特征f_j的猜中近邻相关统计量；diffr_j和diffw_j的计算表达式如下：

式中，s＝min(l1，l2)，i表示对象i，j表示第j个特征，当a＝1，2，3，...，s时，Q_i(a，j)表示[q₁，q₂，...，q_s]，当a＝s+1，s+2，...，2s时，Q_i(a，j)表示[q_s+1,q_s+2，...，q_2s]；当b＝1，2，3，...，s时，Q_i(b，j)表示[q₁,q₂,...，q_s]，当b＝s+1，s+2，...，2s时，Q_i(b，j)表示[q_s+1,q_s+2，...，q_2s]；

(3.5)根据步骤(3.4)计算得到的特征子空间M_l中每个对象的特征的权重，选取e个权重最大的变化特征，进而分别在N个对象中选取所述e个变化特征，结合原变化矢量集合M_d，得到矩阵X，表示如下：

式中，X_i＝[x_i1 x_i2 ... x_ie]，i＝1，2，...，N；矩阵元素x_ij表示对象i的e个变化特征中的第j个特征，e取值小于对象包含的特征总数。

进一步，步骤(4)所述利用PCA算法对矩阵X进行特征降维，得到变换之后的变化特征矩阵Y；方法如下：

(4.1)根据步骤(3)得到的矩阵X，计算矩阵每一行的均值μ_i，i＝1，2，...，N：

式中，j＝1，2，...，e，E(X_i)表示数学期望；

(4.2)根据均值μ_i，计算矩阵X的协方差矩阵C，即：

式中，i＝1，2，...，N，r＝1，2，...，N，j＝1，2，...，e；由上述公式可知，协方差矩阵C为实对称矩阵；

(4.3)计算协方差矩阵C的特征值λ₁，λ₂，...，λ_N，对应的特征向量为v₁，v₂，...，v_N，即λ_iv_i＝Cv_i，i＝1，2，...，N；

(4.4)将特征值λ₁，λ₂，...，λ_N进行降序排序，得到λ′₁，λ′₂，...，λ′_N，对应的特征向量为v′₁，v′₂，...，v′_N，取排序后的特征值构成一个对角矩阵P；

(4.5)计算矩阵X在新的特征向量v′₁，v′₂，...，v′_N上的投影，得到降维后的变化特征矩阵Y，即：

Y＝X*P＝[y₁，y₂，...，y_N]

其中，y₁，y₂，...，y_N表示对象的变化特征。

进一步，所述步骤(5)中使用NR-Kmeans算法对变化特征矩阵Y的元素y_i进行分类，从变化特征矩阵Y中剔除被噪声干扰的稀疏数据，得到密集点集合Y′；步骤如下：

(5.1)在矩阵Y＝[y₁，y₂，...，y_N]中，对于每一个元素点y_α，y_α∈Y，给定邻域半径δ，得到符合

的元素点

Q_b(y_α)表示包含元素点y_α在内的y_α的最近邻元素点集合；其中，

表示元素点y_α、

的距离，b为y_α在邻域范围内最近邻元素点个数，f＝1，2，...，b；

(5.2)计算元素点y_α的密度函数值；元素点y_α的密度函数值表示邻域半径δ范围内所有最近邻元素点对其影响的函数之和；y_α的密度函数值计算公式如下：

式中，

为元素点y_α、

的距离，f＝1，2，...，b；

(5.3)重复步骤(5.1)～(5.2)，直到矩阵Y中所有元素点的密度函数值计算完成；

(5.4)当y_β的密度函数值大于等于矩阵Y中所有元素点平均密度值时，则将元素点y_β视为可用数据放入密集点集合Y′中，即y_β符合：

式中，DF(y_β)表示y_β的密度函数值，DF(y_α)为数据点y_α的密度函数值；

(5.5)根据步骤(5.4)得到矩阵Y中所有的可用数据，所述可用数据构成密集点集合Y′，Y′为剔除了被噪声干扰的稀疏数据后的元素点集合。

进一步，步骤(6)所述将密集点集合Y′分为未变化的类和变化的类，未变化的类元素对应的差值图中的对象的灰度值转化为255，变化的类中元素对应的差值图中的对象的灰度值转化为0，根据对象的灰度值输出变化检测结果图；方法如下：

(6.1)从密集点集合Y′＝[y′₁，y′₂，...，y′_η]，η为Y′中元素点个数，η≤N，y′_η表示剔除了被噪声干扰的数据后的对象的特征；选取密度值最大的元素点，作为第一个初始聚类中心C₁，计算Y′中其他元素点的特征值与聚类中心C₁的特征值的差值，选取与C₁的特征值相差最大的元素点，即差值绝对值最大的点作为第二个初始聚类中心C₂；C₁，C₂相应的特征值为c₁，c₂，聚类中心C₁和C₂对应类D₁和D₂；

(6.2)计算密集点集合Y′中的元素点y′_θ与两个聚类中心的特征值的差值，即

d_θ1＝|y′_θ-c₁|

d_θ2＝|y′_θ-c₂|

d_θ1表示元素点y′_θ与聚类中心C₁的特征差值的绝对值，d_θ2表示元素点y′_θ与聚类中心C₂的特征差值的绝对值，θ＝1，2，...，η；将密集点集合Y′中的元素点分到差值绝对值较小的聚类中心所对应的类；

(6.3)计算类D₁，D₂内相应的元素点特征值的平均值c′₁，c′₂，即

式中，c′₁，c′₂分别为类D₁，D₂相应的元素点特征值的平均值；|D₁|表示类D₁中元素点的个数，|D₂|表示类D₂中元素点的个数，y′对应类内元素点特征值；如果c′₁≠c₁，将类D₁的聚类中心C₁的特征值c₁更新为c′₁；如果c′₂≠c₂，将类D₂的聚类中心C₂的特征值c₂更新为c′₂；

(6.4)重复步骤(6.2)～(6.3)，直至两个聚类中心的特征值都不再发生改变；

(6.5)将两个类D₁，D₂分为未变化的类和变化的类，其中，平均值较小的为未变化的类，平均值较大的为变化的类；未变化的类元素对应的差值图中的对象的灰度值转化为255，变化的类中元素对应的差值图中的对象的灰度值转化为0；

(6.6)根据对象的灰度值输出检测结果图，图中标出了图像变化的区域。

有益效果：与现有技术相比，本发明的技术方案具有以下有益的技术效果：

(1)为了减小噪声对于数据的干扰，使用密度分布函数排除密度稀疏点，剔除一些噪声干扰点。

(2)结合最远距离原则，在高密度数据集合中选取初始聚类中心，避免了原始Kmeans算法随机选取初始聚类中心的缺点。

附图说明

图1是本发明的流程图；

图2是NR-Kmeans算法流程图。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明。

本发明所述的一种鲁棒性遥感图像变化检测方法，如图1所示，包括以下几个步骤：

(1)采集获取目标地区T₁和T₂时刻的遥感图像，所述图像满足分辨率需求并且包含清晰的纹理特征信息；选择合肥市某区域的图像，包含了农田作物面积覆盖情况；

(2)采用ENVI软件实现图像数据预处理，消除拍摄误差；采用eCogntion软件平台对T₁和T₂时刻的预处理后的图像进行多尺度分割，将两幅图像分别划分成N个图像块，即分别得到N个对象；分别获取T₁和T₂时刻的图像的光谱特征和纹理特征，构建特征向量，得到变化矢量集合M_d，M_d＝[m₁,m₂,...m_N]，其中，m_i表示第i个对象所对应的变化矢量；将M_d中元素按照T₁或T₂时刻图像中N个对象的排列方式排列，构成一幅差值图；

(3)根据Relief特征选择算法，获得变化矢量集合M_d的特征子空间M_l中每个对象的特征的权重，选取e个权重最大的变化特征，e取值小于对象包含的特征总数，进而分别在N个对象中选取所述e个变化特征，结合变化矢量集合M_d，得到矩阵X；

(4)利用PCA算法对矩阵X进行特征降维，得到变换之后的变化特征矩阵Y，Y＝[y₁,y₂,...,y_N]；

(5)使用NR-Kmeans算法对变化特征矩阵Y的元素y_i进行分类，从变化特征矩阵Y中剔除被噪声干扰的稀疏数据，得到密集点集合Y'；

(6)将密集点集合Y'分为未变化的类和变化的类，未变化的类元素对应的差值图中的对象的灰度值转化为255，变化的类中元素对应的差值图中的对象的灰度值转化为0，根据对象的灰度值输出变化检测结果图，图中标出了图像变化的区域。

步骤(2)中，采用ENVI软件实现图像数据预处理，消除拍摄误差；采用eCogntion软件平台对T₁和T₂时刻的预处理后的图像进行多尺度分割，将两幅图像分别划分成N个图像块，即分别得到N个对象；分别获取T₁和T₂时刻的图像的光谱特征和纹理特征，构建特征向量，得到变化矢量集合，步骤如下：

(2.1)采用ENVI软件分别对T₁和T₂时刻的图像进行几何校正、图像配准、大气校正操作，实现图像数据预处理；两幅图像大小均为512*512；

(2.2)采用eCogntion软件平台对T₁和T₂时刻的预处理后的图像进行多尺度分割，将两幅图像分别划分成N个图像块，即分别得到N个对象；

(2.3)求取预处理后的图像的Mean-Std特征，作为图像的光谱特征；采用灰度共生矩阵提取预处理后的图像的纹理特征；

(2.4)根据预处理后的两幅图像提取的光谱特征和纹理特征，构建特征向量，获取变化矢量集合M_d。

步骤(2.3)所述求取图像的Mean-Std特征，作为图像的光谱特征；方法如下：

对图像的每个对象，按照下列公式求取光谱特征中的均值Mean和标准差Std：

式中，Mean为对象中像素点灰度的均值，std为对象中像素点灰度的标准差，A表示对象中像素点数，c_i代表对象中第i个像素点的灰度大小。

步骤(2.3)所述采用灰度共生矩阵提取图像的纹理特征；方法如下：

(2.3.1)将图像中任意一像素点(τ₁，τ₂)及另一像素点(τ₁+Δ₁，τ₂+Δ₂)构成点对，其中Δ₁，Δ₂为整数；设像素点(τ₁，τ₂)的灰度值为f1，(τ₁+Δ₁，τ₂+Δ₂)的灰度值为f2，则该点对的灰度值为(f1，f2)，设图像的最大灰度级为L，则f1与f2的组合共有L*L组；本实施例设L＝16，则f1与f2的组合共有16*16组；

(2.3.2)统计图像中每一组不同的(f1，f2)值出现的次数，然后排列成一个矩阵，其中，位于矩阵位置(L，L)上的矩阵元素值就是两个灰度值均为L的点对出现的次数；

(2.3.3)根据(f1，f2)出现的总次数将每一组不同的(f1，f2)值出现的次数归一化为出现的概率g(f1，f2)，归一化之后的的方阵称为灰度共生矩阵；

(2.3.4)灰度共生矩阵提供了14种统计量作为纹理特征，本发明选取6种统计量作为纹理特征，所述6种统计量为对比度、熵、能量、均等性、不相似性、相关性；对比度反映了图像的清晰度和纹理沟纹深浅的程度；熵表示了图像中纹理的非均匀程度或复杂程度；能量反映了图像灰度分布均匀程度和纹理粗细度；均等性反映图像纹理的同质性，度量图像纹理局部变化的多少；相关性反应了图像纹理的一致性；这6种统计量都从各个方面表示了图像的特性，故选择这6种统计量作为纹理特征；

(2.3.5)设灰度共生矩阵的位置(i，j)上的矩阵元素值为g(i，j)，L为图像的最大灰度级，则所述纹理特征表示如下：

式中，Con表示对比度，Ent表示熵，Ene表示能量，Hom表示均等性，Dis表示不相似性，Cor表示相关性；

步骤(2.4)所述根据预处理后的两幅图像提取的光谱特征和纹理特征，获取变化矢量集合M_d；方法如下：

(2.4.1)根据步骤(2.3)得到的T₁和T₂时刻图像的光谱特征和纹理特征，分别构建特征向量，对两幅图像的相对应的特征向量计算差值；两幅图像相对应的第i个对象的第j个特征向量的差值M_d(i，j)，表示如下：

M_d(i，j)＝M_d，1(i，j)-M_d，2(i，j)

式中，M_d，1(i，j)为T₁时刻图像第i个对象的第j个特征向量，M_d，2(i，j)为T₂时刻图像第i个对象的第j个特征向量；i≤N，N为图像中对象的总数；j≤S_feature，S_feature为对象包含的特征总数；S_feature＝S_band×8，S_band为对象总波段数；本实施例中图像有4个波段，则S_feature＝32；

根据差值M_d(i，j)，得到第i个对象所对应的变化矢量m_i，表示如下：

m_i＝(M_d(i，1)，M_d(i，2)，...，M_d(i，32))

(2.4.2)重复步骤(2.4.1)，计算两幅图像中相对应的每一个对象的特征向量的差值，得到每个对象所对应的变化矢量，构成变化矢量集合M_d，表示如下：

M_d＝[m₁，m₂，...m_N]

将M_d中元素按照T₁和T₂时刻图像中N个对象的排列方式排列，构成一幅差值图。

步骤(3)所述根据Relief特征选择算法，获得变化矢量集合M_d的特征子空间M_l中每个对象的特征的权重，选取e个权重最大的变化特征，e取值小于对象包含的特征总数，进而分别在N个对象中选取所述e个变化特征，结合变化矢量集合M_d，得到矩阵X；方法如下：

(3.1)变化矢量集合M_d也称为特征空间，提取特征空间M_d＝[m₁，m₂，m₃，...m_N]的子空间M_l＝[m₁,m₂,m₃,...m_n]作为训练集，其中n＜N；使用人工解译将特征子空间M_l分为M_l1与M_l2两个数集，其中M_l1包含l1个元素，M_l2包含l2个元素，l1+l2＝n，M_l1∈M_l，M_l2∈M_l；

(3.2)计算M_l中所有元素两两之间变化矢量的欧式距离，组成距离矩阵，即

其中，d表示距离矩阵，

表示M_l中元素两两之间变化矢量的欧式距离；

(3.3)在特征子空间M_l中，对于对象i的变化矢量m_i，i＝1，2，3，...，n，若m_i∈M_l1，选取M_l1内与其欧式距离最小的s个矢量，并选取M_l2内与其欧氏距离最大的s个矢量；若m_i∈M_l2，选取M_l2内与其欧式距离最小的s个矢量，并选取M_l1内与其欧氏距离最大的s个矢量；即对m_i求取近邻矩阵Q_i＝[q₁，q₂，...，q_s，...，q_2s]，其中，[q₁，q₂，...，q_s]为所述欧式距离最小的s个矢量，[q_s+1,q_s+2，...，q_2s]为所述欧氏距离最大的s个矢量，s＝min(l1，l2)；

(3.4)根据猜错近邻相关统计量和猜中近邻相关统计量，计算对象i的每个特征对应的权重，即

式中，

是特征f_j对应的权重，diffr_j是特征f_j的猜错近邻相关统计量，diffw_j是特征f_j的猜中近邻相关统计量；diffr_j和diffw_j的计算表达式如下：

式中，s＝min(l1，l2)，i表示对象i，j表示第j个特征，当a＝1，2，3，...，s时，Q_i(a，j)表示[q₁，q₂，...，q_s]，当a＝s+1，s+2，...，2s时，Q_i(a，j)表示[q_s+1,q_s+2，...，q_2s]；当b＝1，2，3，...，s时，Q_i(b，j)表示[q₁,q₂,...，q_s]，当b＝s+1，s+2，...，2s时，Q_i(b，j)表示[q_s+1,q_s+2，...，q_2s]；

(3.5)根据步骤(3.4)计算得到的特征子空间M_l中每个对象的特征的权重，选取e个权重最大的变化特征，进而分别在N个对象中选取所述e个变化特征，结合原变化矢量集合M_d，得到矩阵X，表示如下：

式中，X_i＝[x_i1 x_i2 ... x_ie]，i＝1，2，...，N；矩阵元素x_ij表示对象i的e个变化特征中的第j个特征，e取值小于对象包含的特征总数。

步骤(4)所述利用PCA算法对矩阵X进行特征降维，得到变换之后的变化特征矩阵Y；方法如下：

(4.1)根据步骤(3)得到的矩阵X，计算矩阵每一行的均值μ_i，i＝1，2，...，N：

式中，j＝1，2，...，e，E(X_i)表示数学期望；

(4.2)根据均值μ_i，计算矩阵X的协方差矩阵C，即：

式中，i＝1，2，...，N，r＝1，2，...，N，j＝1，2，...，e；由上述公式可知，协方差矩阵C为实对称矩阵；

(4.3)计算协方差矩阵C的特征值λ₁，λ₂，...，λ_N，对应的特征向量为v₁，v₂，...，v_N，即λ_iv_i＝Cv_i，i＝1，2，…，N；

(4.4)将特征值λ₁，λ₂，...，λ_N进行降序排序，得到λ′₁，λ′₂，...，λ′_N，对应的特征向量为v′₁，v′₂，...，v′_N，取排序后的特征值构成一个对角矩阵P；

(4.5)计算矩阵X在新的特征向量v′₁，v′₂，...，v′_N上的投影，得到降维后的变化特征矩阵Y，即：

Y＝X*P＝[y₁，y₂，...，y_N]

其中，y₁，y₂，...,y_N表示对象的变化特征。

所述步骤(5)中使用NR-Kmeans算法对变化特征矩阵Y的元素y_i进行分类，从变化特征矩阵Y中剔除被噪声干扰的稀疏数据，得到密集点集合Y′；步骤如下：

(5.1)在矩阵Y＝[y₁，y₂，...，y_N]中，对于每一个元素点y_α，y_α∈Y，给定邻域半径δ，得到符合

的元素点

Q_b(y_α)表示包含元素点y_α在内的y_α的最近邻元素点集合；其中，

表示元素点y_α、

的距离，b为y_α在邻域范围内最近邻元素点个数，f＝1，2，...，b；

(5.2)计算元素点y_α的密度函数值；元素点yx的密度函数值表示邻域半径δ范围内所有最近邻元素点对其影响的函数之和；y_α的密度函数值计算公式如下：

式中，

为元素点y_α、

的距离，f＝1，2，...，b；

(5.3)重复步骤(5.1)～(5.2)，直到矩阵Y中所有元素点的密度函数值计算完成；

(5.4)当y_β的密度函数值大于等于矩阵Y中所有元素点平均密度值时，则将元素点y_β视为可用数据放入密集点集合Y′中，即y_β符合：

式中，DF(y_β)表示y_β的密度函数值，DF(y_α)为数据点y_α的密度函数值；

(5.5)根据步骤(5.4)得到矩阵Y中所有的可用数据，所述可用数据构成密集点集合Y′，Y′为剔除了被噪声干扰的稀疏数据后的元素点集合。

步骤(6)所述将密集点集合Y′分为未变化的类和变化的类，未变化的类元素对应的差值图中的对象的灰度值转化为255，变化的类中元素对应的差值图中的对象的灰度值转化为0，根据对象的灰度值输出变化检测结果图；方法如下：

(6.1)从密集点集合Y′＝[y′₁，y′₂，...，y′_η]，η为Y′中元素点个数，η≤N，y′_η表示剔除了被噪声干扰的数据后的对象的特征；选取密度值最大的元素点，作为第一个初始聚类中心C₁，计算Y′中其他元素点的特征值与聚类中心C₁的特征值的差值，选取与C₁的特征值相差最大的元素点，即差值绝对值最大的点作为第二个初始聚类中心C₂；C₁，C₂相应的特征值为c₁，c₂，聚类中心C₁和C₂对应类D₁和D₂；

(6.2)计算密集点集合Y′中的元素点y′_θ与两个聚类中心的特征值的差值，即

d_θ1＝|y′_θ-c₁|

d_θ2＝|y′_θ-c₂|

d_θ1表示元素点y′_θ与聚类中心C₁的差值的绝对值，d_θ2表示元素点y′_θ与聚类中心C₂的差值的绝对值，θ＝1，2，...，η；将密集点集合Y′中的元素点分到差值绝对值较小的聚类中心所对应的类；

(6.3)计算类D₁，D₂内相应的元素点特征值的平均值c′₁，c′₂，即

式中，c′₁，c′₂分别为类D₁，D₂相应的元素点特征值的平均值；|D₁|表示类D₁中元素点的个数，|D₂|表示类D₂中元素点的个数，y′对应类内元素点特征值；如果c′₁≠c₁，将类D₁的聚类中心C₁的特征值c₁更新为c′₁；如果c′₂≠c₂，将类D₂的聚类中心C₂的特征值c₂更新为c′₂；

(6.4)重复步骤(6.2)～(6.3)，直至两个聚类中心的特征值都不再发生改变；

(6.5)将两个类D₁，D₂分为未变化的类和变化的类，其中，平均值较小的为未变化的类，平均值较大的为变化的类；未变化的类元素对应的差值图中的对象的灰度值转化为255，变化的类中元素对应的差值图中的对象的灰度值转化为0；

(6.6)根据对象的灰度值输出检测结果图，图中标出了图像变化的区域。

本发明改进了传统的kmeans聚类算法，检测精度更高，准确率高。

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

Claims (8)

1.一种鲁棒性遥感图像变化检测方法，其特征在于：该方法包括如下步骤：

(1)采集获取目标地区T₁和T₂时刻的遥感图像，所述图像满足分辨率需求并且包含清晰的纹理特征信息；

(2)采用ENVI软件实现图像数据预处理，消除拍摄误差；采用eCogntion软件平台对T₁和T₂时刻的预处理后的图像进行多尺度分割，将两幅图像分别划分成N个图像块，即分别得到N个对象；分别获取T₁和T₂时刻的图像的光谱特征和纹理特征，构建特征向量，得到变化矢量集合M_d，M_d＝[m₁,m₂,...m_N]，其中，m_N表示第N个对象所对应的变化矢量；将M_d中元素按照T₁或T₂时刻图像中N个对象的排列方式排列，构成一幅差值图；

(3)根据Relief特征选择算法，获得变化矢量集合M_d的特征子空间M_l中每个对象的特征的权重，选取e个权重最大的变化特征，e取值小于对象包含的特征总数，进而分别在N个对象中选取所述e个变化特征，结合变化矢量集合M_d，得到矩阵X；

(4)利用PCA算法对矩阵X进行特征降维，得到变换之后的变化特征矩阵Y，Y＝[y₁,y₂,...,y_N]；

(5)使用NR-Kmeans算法对变化特征矩阵Y的元素y_i进行分类，从变化特征矩阵Y中剔除被噪声干扰的稀疏数据，得到密集点集合Y'；步骤如下：

(5.1)在矩阵Y＝[y₁,y₂,...,y_N]中，对于每一个元素点y_α，y_α∈Y，给定邻域半径δ，得到符合

的元素点

Q_b(y_α)表示包含元素点y_α在内的y_α的最近邻元素点集合；其中，

表示元素点y_α、

的距离，b为y_α在邻域范围内最近邻元素点个数，f＝1,2,...,b；

(5.2)计算元素点y_α的密度函数值；元素点y_α的密度函数值表示邻域半径δ范围内所有最近邻元素点对其影响的函数之和；y_α的密度函数值计算公式如下：

式中，

为元素点y_α、

的距离，f＝1,2,…,b；

(5.3)重复步骤(5.1)～(5.2)，直到矩阵Y中所有元素点的密度函数值计算完成；

(5.4)当y_β的密度函数值大于等于矩阵Y中所有元素点平均密度值时，则将元素点y_β视为可用数据放入密集点集合Y'中，即y_β符合：

式中，DF(y_β)表示y_β的密度函数值，DF(y_α)为数据点y_α的密度函数值；

(5.5)根据步骤(5.4)得到矩阵Y中所有的可用数据，所述可用数据构成密集点集合Y'，Y'为剔除了被噪声干扰的稀疏数据后的元素点集合；

(6)将密集点集合Y'分为未变化的类和变化的类，未变化的类元素对应的差值图中的对象的灰度值转化为255，变化的类中元素对应的差值图中的对象的灰度值转化为0，根据对象的灰度值输出变化检测结果图，图中标出了图像变化的区域。

2.根据权利要求1所述的一种鲁棒性遥感图像变化检测方法，其特征在于：步骤(2)中，采用ENVI软件实现图像数据预处理，消除拍摄误差；采用eCogntion软件平台对T₁和T₂时刻的预处理后的图像进行多尺度分割，将两幅图像分别划分成N个图像块，即分别得到N个对象；分别获取T₁和T₂时刻的图像的光谱特征和纹理特征，构建特征向量，得到变化矢量集合，步骤如下：

(2.1)采用ENVI软件分别对T₁和T₂时刻的图像进行几何校正、图像配准、大气校正操作，实现图像数据预处理；

(2.2)采用eCogntion软件平台对T₁和T₂时刻的预处理后的图像进行多尺度分割，将两幅图像分别划分成N个图像块，即分别得到N个对象；

(2.3)求取预处理后的图像的Mean-Std特征，作为图像的光谱特征；采用灰度共生矩阵提取预处理后的图像的纹理特征；

(2.4)根据预处理后的两幅图像提取的光谱特征和纹理特征，构建特征向量，获取变化矢量集合M_d。

3.根据权利要求2所述的一种鲁棒性遥感图像变化检测方法，其特征在于：步骤(2.3)求取图像的Mean-Std特征，作为图像的光谱特征；方法如下：

对图像的每个对象，按照下列公式求取光谱特征中的均值Mean和标准差Std：

式中，Mean为对象中像素点灰度的均值，std为对象中像素点灰度的标准差，A表示对象中像素点数，c_p代表对象中第p个像素点的灰度大小。

4.根据权利要求2所述的一种鲁棒性遥感图像变化检测方法，其特征在于：步骤(2.3)所述采用灰度共生矩阵提取图像的纹理特征；方法如下：

(2.3.1)将图像中任意一像素点(τ₁,τ₂)及另一像素点(τ₁+△₁,τ₂+△₂)构成点对，其中△₁,△₂为整数；设像素点(τ₁,τ₂)的灰度值为f1，(τ₁+△₁,τ₂+△₂)的灰度值为f2，则该点对的灰度值为(f1,f2)，设图像的最大灰度级为L，则f1与f2的组合共有L*L组；

(2.3.2)统计图像中每一组不同的(f1,f2)值出现的次数，然后排列成一个矩阵，其中，位于矩阵位置(L，L)上的矩阵元素值就是两个灰度值均为L的点对出现的次数；

(2.3.3)根据(f1,f2)出现的总次数将每一组不同的(f1,f2)值出现的次数归一化为出现的概率g(f1,f2)，归一化之后的的方阵称为灰度共生矩阵；

(2.3.4)选取灰度共生矩阵所提供的6种统计量作为纹理特征，所述6种统计量为对比度、熵、能量、均等性、不相似性、相关性；

(2.3.5)设灰度共生矩阵的位置(z,t)上的矩阵元素值为g(z,t)，L为图像的最大灰度级，则所述纹理特征表示如下：

式中，Con表示对比度，Ent表示熵，Ene表示能量，Hom表示均等性，Dis表示不相似性，Cor表示相关性；

5.根据权利要求2所述的一种鲁棒性遥感图像变化检测方法，其特征在于：步骤(2.4)所述根据预处理后的两幅图像提取的光谱特征和纹理特征，获取变化矢量集合M_d；方法如下：

(2.4.1)根据步骤(2.3)得到的T₁和T₂时刻图像的光谱特征和纹理特征，分别构建特征向量，对两幅图像的相对应的特征向量计算差值；两幅图像相对应的第i个对象的第j个特征向量的差值M_d(i,j)，表示如下：

M_d(i,j)＝M_d,1(i,j)-M_d,2(i,j)

式中，M_d,1(i,j)为T₁时刻图像第i个对象的第j个特征向量，M_d,2(i,j)为T₂时刻图像第i个对象的第j个特征向量；i≤N,N为图像中对象的总数；j≤S_feature,S_feature为对象包含的特征总数；S_feature＝S_band×8,S_band为对象总波段数；

根据差值M_d(i,j)，得到第i个对象所对应的变化矢量m_i，表示如下：

m_i＝(M_d(i,1),M_d(i,2),...,M_d(i,S_feature))

(2.4.2)重复步骤(2.4.1)，计算两幅图像中相对应的每一个对象的特征向量的差值，得到每个对象所对应的变化矢量，构成变化矢量集合M_d，表示如下：

M_d＝[m₁,m₂,...m_N]

将M_d中元素按照T₁和T₂时刻图像中N个对象的排列方式排列，构成一幅差值图。

6.根据权利要求1-5任一所述的一种鲁棒性遥感图像变化检测方法，其特征在于：步骤(3)所述根据Relief特征选择算法，获得变化矢量集合M_d的特征子空间M_l中每个对象的特征的权重，选取e个权重最大的变化特征，e取值小于对象包含的特征总数，进而分别在N个对象中选取所述e个变化特征，结合变化矢量集合M_d，得到矩阵X；方法如下：

(3.1)变化矢量集合M_d也称为特征空间，提取特征空间M_d＝[m₁,m₂,m₃,...m_N]的子空间M_l＝[m₁,m₂,m₃,...m_n]作为训练集，其中n＜N；将特征子空间M_l分为M_l1与M_l2两个数集，其中M_l1包含l1个元素，M_l2包含l2个元素，l1+l2＝n，M_l1∈M_l,M_l2∈M_l；

(3.2)计算M_l中所有元素两两之间变化矢量的欧式距离，组成距离矩阵，即：

其中，d表示距离矩阵，

表示M_l中元素两两之间变化矢量的欧式距离；

(3.3)在特征子空间M_l中，对于对象i的变化矢量m_i，i＝1,2,3,…,n，若m_i∈M_l1，选取M_l1内与其欧式距离最小的s个矢量，并选取M_l2内与其欧氏距离最大的s个矢量；若m_i∈M_l2，选取M_l2内与其欧式距离最小的s个矢量，并选取M_l1内与其欧氏距离最大的s个矢量；即对m_i求取近邻矩阵Q_i＝[q₁,q₂,...,q_s,...,q_2s]，其中，[q₁,q₂,...,q_s]为所述欧式距离最小的s个矢量，[q_s+1,q_s+2,...,q_2s]为所述欧氏距离最大的s个矢量，s＝min(l1,l2)；

(3.4)根据猜错近邻相关统计量和猜中近邻相关统计量，计算对象i的每个特征对应的权重，即

w_fj＝diffr_j-diffw_j

式中，w_fj是特征f_j对应的权重，diffr_j是特征f_j的猜错近邻相关统计量，diffw_j是特征f_j的猜中近邻相关统计量；diffr_j和diffw_j的计算表达式如下：

式中，s＝min(l1,l2)，i表示对象i，j表示第j个特征，当a＝1,2,3,…,s时，Q_i(a,j)表示[q₁,q₂,...,q_s]，当a＝s+1,s+2,…,2s时，Q_i(a,j)表示[q_s+1,q_s+2,...,q_2s]；当b＝1,2,3,…,s时，Q_i(b,j)表示[q₁,q₂,...,q_s]，当b＝s+1,s+2,…,2s时，Q_i(b,j)表示[q_s+1,q_s+2,...,q_2s]；

(3.5)根据步骤(3.4)计算得到的特征子空间M_l中每个对象的特征的权重，选取e个权重最大的变化特征，进而分别在N个对象中选取所述e个变化特征，结合原变化矢量集合M_d，得到矩阵X，表示如下：

式中，X_i＝[x_i1 x_i2...x_ie]，i＝1,2,...,N；矩阵元素x_ij表示对象i的e个变化特征中的第j个特征，e取值小于对象包含的特征总数。

7.根据权利要求1-5任一所述的一种鲁棒性遥感图像变化检测方法，其特征在于：步骤(4)所述利用PCA算法对矩阵X进行特征降维，得到变换之后的变化特征矩阵Y；方法如下：

(4.1)根据步骤(3)得到的矩阵X，计算矩阵每一行的均值μ_i,i＝1,2,...,N：

式中，j＝1,2,...,e，E(X_i)表示数学期望；

(4.2)根据均值μ_i，计算矩阵X的协方差矩阵C，即：

式中，i＝1,2,...,N，r＝1,2,...,N，j＝1,2,...,e；

(4.3)计算协方差矩阵C的特征值λ₁,λ₂,...,λ_N，对应的特征向量为ν₁,ν₂,...,ν_N，即λ_iv_i＝Cv_i,i＝1,2,…,N；

(4.4)将特征值λ₁,λ₂,...,λ_N进行降序排序，得到λ'₁,λ'₂,...,λ'_N，对应的特征向量为ν'₁,ν'₂,...,ν'_N，取排序后的特征值构成一个对角矩阵P；

(4.5)计算矩阵X在新的特征向量ν'₁,ν'₂,...,ν'_N上的投影，得到降维后的变化特征矩阵Y，即：

Y＝X*P＝[y₁,y₂,...,y_N]

其中，y₁,y₂,...,y_N表示对象的变化特征。

8.根据权利要求1-5任一所述的一种鲁棒性遥感图像变化检测方法，其特征在于：步骤(6)所述将密集点集合Y'分为未变化的类和变化的类，未变化的类元素对应的差值图中的对象的灰度值转化为255，变化的类中元素对应的差值图中的对象的灰度值转化为0，根据对象的灰度值输出变化检测结果图；方法如下：

(6.1)从密集点集合Y'＝[y'₁,y'₂,...,y'_η]，η为Y'中元素点个数，η≤N，y'_η表示剔除了被噪声干扰的数据后的对象的特征；选取密度值最大的元素点，作为第一个初始聚类中心C₁，计算Y'中其他元素点的特征值与聚类中心C₁的特征值的差值，选取与C₁的特征值相差最大的元素点，即差值绝对值最大的点作为第二个初始聚类中心C₂；C₁,C₂相应的特征值为c₁,c₂，聚类中心C₁和C₂对应类D₁和D₂；

(6.2)计算密集点集合Y'中的元素点y'_θ与两个聚类中心的特征值的差值，即

d_θ1＝|y'_θ-c₁|

d_θ2＝|y'_θ-c₂|

d_θ1表示元素点y'_θ与聚类中心C₁的特征差值的绝对值，d_θ2表示元素点y'_θ与聚类中心C₂的特征差值的绝对值，θ＝1,2,...,η；将密集点集合Y'中的元素点分到差值绝对值较小的聚类中心所对应的类；

(6.3)计算类D₁,D₂内相应的元素点特征值的平均值c'₁,c'₂，即

式中，c'₁,c'₂分别为类D₁,D₂相应的元素点特征值的平均值；|D₁|表示类D₁中元素点的个数，|D₂|表示类D₂中元素点的个数，y'对应类内元素点特征值；如果c'₁≠c₁，将类D₁的聚类中心C₁的特征值c₁更新为c'₁；如果c'₂≠c₂，将类D₂的聚类中心C₂的特征值c₂更新为c'₂；

(6.4)重复步骤(6.2)～(6.3)，直至两个聚类中心的特征值都不再发生改变；

(6.5)将两个类D₁,D₂分为未变化的类和变化的类，其中，平均值较小的为未变化的类，平均值较大的为变化的类；未变化的类元素对应的差值图中的对象的灰度值转化为255，变化的类中元素对应的差值图中的对象的灰度值转化为0；

(6.6)根据对象的灰度值输出检测结果图，图中标出了图像变化的区域。

Images