CN104008376A - 基于可能性中心点聚类的多光谱遥感影像混合像元分解方法 - Google Patents

Abstract

本发明属于遥感图像处理技术领域，具体为一种基于可能性中心点聚类的多光谱遥感影像混合像元分解方法。本发明根据可能性理论和中心点聚类方法的基本原理，构建可能性C中心点聚类方法的代价函数式和可能性隶属度函数，运用启发式算法对离散的解空间进行搜索，通过迭代运算更新聚类中心及可能性隶属度矩阵，获得最终的光谱端元和组分比例，即地表覆被的盖度。本方法不仅对噪声和孤立点数据有良好的鲁棒性，而且可有效处理重叠聚类问题，在聚类间彼此距离较近时也能准确地识别聚类中心，因此，可以在噪声环境下获得高精度的地表覆被盖度和端元光谱信息。本发明在采用多光谱遥感影像进行高精度地物分类、目标检测方面具有重要的应用价值。

Description

基于可能性中心点聚类的多光谱遥感影像混合像元分解方法

技术领域

本发明属于遥感图像处理技术领域，具体涉及一种基于可能性中心点聚类的从多光谱遥感影像中准确提取地表覆被盖度的混合像元分解新方法。

背景技术

高精度的遥感解译信息是各类地学空间分析和过程演化研究的重要数据源，如何从遥感数据中获取高精度的解译信息是遥感研究和应用中的基本问题。遥感影像像元记录的是探测系统的瞬时视场角所对应的地面范围内目标的辐射能量总和[1]。如果瞬时视场(IFOV)仅包含同一类性质的目标，则该像元记录的是同一性质地面目标的辐射能量的总和，这样的像元称为纯像元；如果IFOV中包含不止一种土地覆盖类型或要素，则该像元记录的是多类不同性质的地面目标的辐射能量总和，这样的像元称为混合像元。遥感影像通常是由纯像元和混合像元组合而成，其组合程度取决于IFOV和地面特征的空间复杂性。在通常情况下，将每个混合像元赋予单一地表覆盖类型的硬分类方法往往会发生误分类，并且由于该方法易于忽视小尺寸地物、无法准确识别地类边界，因而难以提供较为精确的地表覆盖信息。与传统的硬分类方法相比，混合像元分解方法可以预测得到遥感影像中每一像元内不同组分端元的盖度信息，从而有效降低遥感影像的分类误差和不确定性，提高地表覆盖目标的识别精度。因此，混合像元分解已成为遥感影像解译与定量研究中的重要方向之一。

现有的混合像元分类技术，如光谱混合模型[2]、模糊分类[3]、人工神经网络[4]、支持向量机[5]等不断涌现，它们通过建立不同的模型对每一像元内的地表覆盖类型信息进行解混以预测该像元内各类型组分的比例。在这些分类方法中，光谱混合模型和模糊分类方法由于被集成于各种遥感图像处理软件中而得到了广泛地应用。光谱混合分析使用线性混合模型，假设图像中的光谱变异是由有限个地表覆盖类型的混合所引起的，即地面上被观测到的区域光谱响应模式是区域内各种地表覆盖类型组分(也称端元组分)的光谱信号的线性混合。线性光谱混合模型没有考虑某些诸如多光谱反射的因素，这种反射会导致光谱混合过程中复杂的非线性过程。也就是说，从像元中观测到的光谱信号可能包括来自各种端元组分光谱信号的混合物，但同时可能也包括地物要素间额外的多时反射信号，例如，植被和土壤之间的多次交叉反射、散射作用等。对于这种情况就需要应用非线性光谱混合模型[6]。

模糊分类方法通过计算每个像元属于各光谱端元的隶属度值来获取每一端元组分在像元内所占的比例，并要求像元中的所有可能端元组分之和必须为1。在各种模糊分类方法中，模糊C均值(FCM)算法作为一种模糊非监督分类法被较多地应用于遥感影像的混合像元分类。根据FCM算法的约束条件[3]，其计算得到的模糊隶属度值无法准确表示数据点对聚类中心的归属度(belonging)和兼容度(compatibility)，而是代表它在类别间的共享(sharing)程度，并且对于噪声较多的数据集，该算法的分类效果较差。为了克服FCM算法的缺陷，Krishnapuram和Keller[7]在可信度及可能性理论[8][9]的基础上，通过构建可能性隶属度函数和目标函数，提出了可能性C均值(PCM)聚类法，由该法所得的可能性隶属度值能很好反映数据点对聚类中心的归属度和兼容度，并对噪声和孤立点数据具有一定的鲁棒性。但是，PCM对于聚类中心的初值十分敏感，并且在聚类间彼此距离较近、出现一定程度重叠时，该算法往往会收敛于相同的聚类中心。为了避免产生一致的聚类中心，出现了一些改进的PCM算法，然而相关实验表明，这些方法并不能很好地解决这一问题[10][11][12]。因此，上述方法都有成功应用的一面，但由于这些算法本身存在的固有缺陷，其分类精度和效率与实际应用要求仍存在一定差距。

发明内容

为了克服现有技术的不足，本发明提出一种多光谱遥感影像混合像元分解方法，本发明可以在噪声环境下获得高精度的地表覆被分类结果和端元光谱信息，为遥感影像混合像元分解提供了一种有效的新方法。

本发明根据可能性理论和中心点聚类方法的基本原理，构建可能性C中心点聚类(PCRMDD)方法的代价函数和可能性隶属度函数；为提高算法的效率，运用启发式算法对离散的解空间进行搜索，通过迭代运算更新聚类中心及可能性隶属度矩阵，将代价函数达到最小值时的聚类中心和隶属度归一化结果作为光谱端元和组分比例。由该方法所得的可能性隶属度值能很好反映数据点对聚类中心的归属度(belonging)和匹配度(compatibility)，并对噪声和孤立点数据有良好的鲁棒性。此外，本发明由于可以有效处理重叠聚类问题，在聚类间彼此距离较近时也能准确地识别聚类中心，其分类结果则能正确地反映数据集的聚类结构，本发明提出的遥感影像混合像元分解方法，基本内容如下：

1.初始聚类中心的提取

首先利用减法聚类法(Subtractive Clustering)获取初始聚类中心。其核心思想是根据各个数据点周围的数据点密度来计算该点作为聚类中心的可能性，被选为聚类中心的数据点具有最高的密度值，并通过修正所有数据点的密度值，排除聚类中心附近一定范围内的其他点作为聚类中心的可能，再以剩下的数据点重复这一过程，直到所有剩余数据点的密度值小于某一阈值或满足终止准则为止[13]。减法聚类法是一种简单而高效的单次(one-pass)算法，由该算法得到的聚类估计可用于初始化那些基于重复优化过程的聚类方法，其具体步骤如下：

考虑M维空间的n个数据点(x₁,x₂,…,x_n)，并假定所有数据点已归一化到一个单位超立方体中。该法将每个数据点作为聚类中心的候选者，并构造一个山峰函数来表示各数据点一定邻域以内数据点的密集程度，数据点x_i处的密度指标，即山峰函数的高度定义为：

$D_i=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\exp [-\frac{{\left|\right|x_i-x_j\left|\right|}^2}{{(r_a/2)}^2}]i=\mathrm{1,2},...,n---\left(1\right)$

其中，r_a是一个正常数，定义了该点的一个邻域，该邻域以外的数据点对该点的密度指标贡献甚微。显然，如果一个数据点有多个邻近的数据点，则该数据点具有高密度值，其作为聚类中心的可能性也越大。常量r_a同时影响山峰函数的高度和光滑程度，当数据集有足够的规模、且能很好地聚类时，聚类结果通常不会对r_a的值敏感，但是在实际应用中，必需根据数据规模和输入维数来选择r_a的取值。

在计算得到每个数据点的密度指标后，选择具有最高密度指标的数据点作为第一个聚类中心，令为选中的点，为其密度指标。然后，为了消除第一个聚类中心密度指标的影响，对每个数据点x_i的密度指标用公式(2)进行修正。

$D_i=D_i-D_{C_1}\exp [-\frac{{\left|\right|x_i-x_{C_1}\left|\right|}^2}{{(r_b/2)}^2}]---\left(2\right)$

其中，r_b是一个正常数，定义了一个密度指标显著减小的邻域，r_b通常大于r_a，以避免出现相距很近的聚类中心。显然，靠近第一个聚类中心的数据点的密度指标将显著减小，这样使得这些点不太可能被选为下一个聚类中心。根据文献[13]所建议，一般取r_b＝1.5r_a。

在修正了每个数据点的密度指标后，继续采用类似的方法选定下一个聚类中心，再次修正所有数据点的密度指标。该过程不断重复，直至产生足够多的聚类中心或满足终止准则，即(δ为一较小分数)，δ的取值会影响到聚类个数，δ的值越小，聚类数就越多，反之，聚类数越少。

2.可能性C中心点聚类(PCRMDD)方法

为有效处理重叠聚类问题，并在噪声环境下提高可能性聚类的精度，本发明基于可能性理论和中心点聚类方法的基本原理，建立了可能性C中心点(Possibilistic CRepulsive Medoids，简称PCRMDD)算法。

本发明定义

$J_3(U,V)=\underset{i=1}{\overset{c}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{u}_{\mathrm{ij}}^m{\left|\right|x_j-v_i\left|\right|}^2+\underset{i=1}{\overset{c}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\frac{{\mathrm{\η}}_i}{m}(u_{\mathrm{ij}}^m\ln u_{\mathrm{ij}}^m-u_{\mathrm{ij}}^m)---\left(3\right)$

作为PCRMDD算法的目标函数。其中V＝{v₁,v₂,…,v_c}为聚类中心，c(c≥2)为类别数，x_j为数据对象，x_j∈R^s，j＝1,2,…,n，m为模糊加权指数，||x_j-v_i||²为x_j和v_i之间的欧几里得距离平方，η_i为一正常数，它与第i个聚类的形状和大小有关，决定了特定点对第i个聚类中心估计的影响程度，U是一个c×n的二维隶属度矩阵，u_ij表示x_j属于第i个聚类的可能性，满足如下约束条件：

u_ij∈[0,1],i＝1,...,c,j＝1,...,n

$0<\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}\&le;n,$ i＝1,...,c

且 $\underset{i}{\max }{u}_{\mathrm{ij}}>0,j=1,...,n---\left(4\right)$

因此，隶属度u_ij只取决于x_j和第i个聚类中心，而与x_j属于其他聚类的可能性无关。

由于本发明是将数据集中的若干个数据对象(即中心点)作为PCRMDD算法中的聚类中心，即该算法的解空间是离散的，因而无法通过计算目标函数(3)对中心点V的偏导数来获取算法的最优解。为了提高算法的效率，运用启发式算法[14]对数据集X进行搜索，得到J₃(U,V)取最小值时的中心点V，此时，划分矩阵U的表达式为

$u_{\mathrm{ij}}=\exp (-\frac{{\left|\right|x_j-v_i\left|\right|}^2}{{\mathrm{\&eta;}}_i})---\left(5\right)$

将隶属度函数表达式(5)代入目标函数(3)式中，得到

$J_3(U,V)=-\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&eta;}}_i}{m}\exp (-\frac{m}{{\mathrm{\&eta;}}_i}\&times;{\left|\right|x_j-v_i\left|\right|}^2)=\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{h}_j---\left(6\right)$

式中， $h_j=-\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&eta;}}_i}{m}\exp (-\frac{m}{{\mathrm{\&eta;}}_i}\&times;{\left|\right|x_j-v_i\left|\right|}^2).$

本发明方法的具体步骤如下：

(1)确定聚类个数C、模糊加权指数m和最大迭代次数max_iter，运用减法聚类法从数据集X中选择C个初始聚类中心V＝{v₁,v₂,…,v_c}，将迭代次数iter设为0；

(2)对所有数据对象计算h_j(j＝1,…,n)，并对h_j进行升序排列构成h_j:n，选取对应于前s个h_j:n的s个数据对象，其中，h_j:n代表将h_j进行升序排列后得到的第j项；

(3)根据式(5)对这s个数据对象计算其隶属度值u_ij:n(i＝1,…,C；j＝1,…,s)；

(4)储存当前聚类中心V_current＝V；

(5)计算新的聚类中心v_i＝x_q(i＝1,…,C)，

其中， $q=\underset{1\&le;k\&le;s}{\arg \min }[\underset{j=1}{\overset{s}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}:n}^m{\left|\right|x_{j:n}-x_{k:n}\left|\right|}^2+\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&eta;}}_i}{m}(u_{\mathrm{ij}:n}^m\ln u_{\mathrm{ij}:n}^m-u_{\mathrm{ij}:n}^m)];$

(6)将iter增1，即iter＝iter+1。如果聚类中心不再发生变化，或者iter＝max_iter，那么终止迭代过程，即可获得可能性隶属度矩阵U与聚类中心V。否则，转向第(2)步。

如果上述步骤中的h_j和q分别改为

$h_j={\left(\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{\left|\right|x_j-v_i\left|\right|}^{\frac{2}{1-m}}\right)}^{1-m}---\left(7\right)$

$q=\underset{1\&le;k\&le;s}{\arg \min }\left[\underset{j=1}{\overset{s}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}:n}^m{\left|\right|x_{j:n}-x_{k:n}\left|\right|}^2\right]---\left(8\right)$

并根据式(9)

$u_{\mathrm{ij}:n}=\frac{{\left|\right|x_{j:n}-v_i\left|\right|}^{-2/(m-1)}}{\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}{\left|\right|x_{j:n}-v_k\left|\right|}^{-2/(m-1)}},i=1,...,C,j=1,...,s---\left(9\right)$

计算模糊隶属度值u_ij:n，则可能性C中心点(PCRMDD)算法将转变为模糊C中心点(FCTMDD)算法[15]。

本发明中，在具体计算过程中，估计η_i值的方法是使用模糊C中心点(FCTMDD)算法，在该算法收敛后得到η_i估计值为

${\mathrm{\&eta;}}_i=\frac{\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{\left|\right|x_j-v_i\left|\right|}^2}{\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m}---\left(10\right)$

或

${\mathrm{\&eta;}}_i=\frac{\underset{u_{\mathrm{ij}}>\mathrm{\&alpha;}}{\mathrm{\&Sigma;}}{{\left|\right|x}_j-v_j\left|\right|}^2}{\underset{u_{\mathrm{ij}}>\mathrm{\&alpha;}}{\mathrm{\&Sigma;}}1}---\left(11\right)$

其中α为合适的阈值，通常在[0.1,0.4]区间内取值。

此外，在迭代运算中，我们始终考虑对应于前s个h_j:n的s个数据对象，若取s＝80％×n，则其余占总数20％、对应于较大h_j的数据对象被视为噪声数据，不参与运算，这样我们仅以前80％的数据求取目标函数的最小值来获得算法的最优解，因而，PCRMDD算法对于噪声数据具有很好的鲁棒性。

为了将PCRMDD算法应用于遥感影像混合像元分解，本发明对可能性隶属度u_ij进行归一化处理，并将噪声数据的隶属度值设为0，其计算公式为

$u_{\mathrm{ij}}^{\&prime;}=\left\{\begin{array}{cc}\frac{u_{\mathrm{ij}:n}}{\underset{i=1}{\overset{C}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}:n}}& j=1,...,s\\ 0& j=s+1,...,n\end{array}\right.---\left(12\right)$

本发明的优点在于：本发明提供了一种基于可能性C中心点聚类(PCRMDD)算法的多光谱遥感影像混合像元分解新方法。基于可能性理论和中心点聚类方法的基本原理，本发明提出了可能性聚类的代价函数式，计算得到的可能性隶属度值能很好反映数据点对聚类中心的归属度和匹配度，在噪声较多的情况下可取得较高的分类精度。同时，本方法克服了模糊聚类算法易于陷入局部最优解的缺点，并且在聚类间彼此距离较近、出现一定程度重叠时也能准确地识别出聚类中心。模拟数据实验发现，所提出的算法对噪声和孤立点数据具有良好的鲁棒性，与现有的混合像元分解方法相比，PCRMDD算法的性能最优。对于实际遥感影像实验的结果表明，本方法能在噪声环境下获得精度较高的地表覆被盖度和端元光谱信息。本发明在采用多光谱遥感影像进行高精度地物分类、目标检测方面具有重要的应用价值。

附图说明

图1为Bridge数据的聚类结果，(a)为FCM算法的硬分类结果，(b)为PCM、FCTMDD和PCRMDD算法的硬分类结果。

图2为含噪声数据的Bridge数据集的聚类结果，(a)为FCM算法的硬分类结果，(b)为PCM、FCTMDD和PCRMDD算法的硬分类结果。

图3为Iris数据集及其聚类结果，(a)为Iris数据集，(b)为PCM算法[7]的硬分类结果，(c)为改进的PCM算法[10]的硬分类结果，(d)为FCTMDD算法的硬分类结果，(e)为PCRMDD算法的硬分类结果。

图4为研究区遥感影像混合像元分解结果，(a)为低反照度建筑，(b)为高反照度建筑，(c)为水体，(d)为植被，(e)为土壤。

图5为研究区各地类端元反射光谱曲线。

图6为均方根误差图像(a)及其频率分布图(b)。

图7为高反照度建筑的盖度估算结果与实际面积比率之间的差异情况。

具体实施方式

下面，分别以模拟数据和实际遥感影像数据为例说明具体的实施方式。

1.模拟数据实验

为检验PCRMDD算法的结果及其有效性，采用两组模拟数据为实验数据，通过引入划分系数(PC)、划分熵系数(PE)、FS指数和XB指数四种指标作为聚类有效性测度标准，将PCRMDD算法与现有的混合像元分解方法(如FCM、PCM和FCTMDD)的聚类效果进行比较分析来验证PCRMDD算法的性能。

聚类有效性测度

(1)划分系数(Partition Coefficient)[16]，定义为

$\mathrm{PC}\left(c\right)=\frac{1}{n}\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^2---\left(13\right)$

PC的取值范围是[1/c,1]，PC的值越接近1，聚类结果越硬，即不同聚类间X中的共享向量部分越小；而当所有的u_ij值相等时，得到的PC值最小为1/c，因而PC的值越接近于1/c，聚类结果越模糊，这意味着或者X不具有聚类结构，或者采用的聚类算法不能识别出X所具有的聚类集合。

(2)划分熵系数(Partition Entropy)[17]，定义为

$\mathrm{PE}\left(c\right)=-\frac{1}{n}\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}{\log }_2{u}_{\mathrm{ij}}---\left(14\right)$

PE的取值范围是[0,log₂c]，PE的值越接近0，聚类结果越硬；另一方面，PE的值越接近于log₂c，聚类结果越模糊，同前者一样，这意味着或者X不具有聚类结构，或者采用的聚类算法不能识别出X所具有的聚类集合。

PC和PE指标可以用来测量聚类间的重叠程度，它们仅使用了隶属度值，而与数据的几何结构无关，同时，它们对于模糊度m也颇为敏感。而下面描述的FS指数和XB指数则同时考虑了隶属度值和数据的几何特征。

(3)FS指数(Fukuyama-Sugeno)[18]，定义为

$\begin{array}{c}\mathrm{FS}\left(c\right)=\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{\left|\right|x_j-v_i\left|\right|}^2-\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{\left|\right|v_i-\stackrel{\&OverBar;}{v}\left|\right|}^2\\ =J_{\mathrm{FCM}}(u,v)+K(u,v)\end{array}---\left(15\right)$

其中，FS指数可以分解为J_FCM(u,v)和K(u,v)两项，前一项是FCM算法的代价函数，用于测量聚类的致密度，而后一项用于测量分离度，FS值越小，聚类的性能越优，即算法获得的聚类致密且分离性较好。

(4)XB指数(Xie-Beni)[19]，定义为

$\begin{array}{c}\mathrm{XB}\left(c\right)=\frac{\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{\left|\right|x_j-v_i\left|\right|}^2}{n{\min }_{i,j=1,...,m,i\&NotEqual;j}{\left|\right|v_i-v_j\left|\right|}^2}\\ =\frac{J_{\mathrm{FCM}}(u,v)/n}{\mathrm{Sep}\left(v\right)}\end{array}---\left(16\right)$

其中，J_FCM(u,v)是聚类致密度的测度，Sep(v)是分离度的测度，与FS指数类似，XB值越小，聚类的性能越优，即算法获得的聚类致密且分离性较好。同样，由于XB和FS指数对于模糊度m敏感，在实际计算中，将m值限制为2。

1)实验一：Bridge数据

在该实验中，利用文献[7]中的Bridge数据，对FCM、PCM、FCTMDD和PCRMDD算法的运行结果进行比较。由图1所示，该数据集形成了两个明显的聚类，各类分别包含7个数据点。图1(a)为运行FCM后获得的硬分类结果，图1(b)为PCM、FCTMDD和PCRMDD的运行结果，图中，黑色的圆形和矩形分别代表两个聚类中心，根据计算得到的划分矩阵(表1和2)将每个数据点分配到最大隶属度值所对应的类中，从而得到硬分类结果。由图1和表1可见，PCM、FCTMDD和PCRMDD得到了基本一致的聚类中心，与可能性聚类法(PCM和PCRMDD)相比，模糊聚类法(FCM和FCTMDD)获得隶属度值接近于1或0，因此，这类算法难以区分各个类中成员距聚类中心的远近，相反，可能性聚类法则能提供数据点的位置信息，即数据点离开聚类中心越远，其隶属度值就越小；同时，模糊聚类法的划分矩阵表现出一种非对称性，并且当聚类彼此不断靠近时，这种非对称性则越为明显，因此，FCM得到的聚类中心出现向另一类靠近的现象。此外，由聚类中心的位置坐标可以看出，中心点聚类法(FCTMDD和PCRMDD)的结果分离性较好，在聚类间彼此距离较近时也能准确地识别聚类中心。当在图1的数据集中加入连接点A和孤立点B时(图2)，FCM由于受到新加入数据点的影响，其聚类中心向两点靠近，并受初始化影响，FCM分别将A和B分配至不同的类中；而PCM、FCTMDD和PCRMDD的硬分类结果(图2(b))与图1(b)所示情况基本相同，由此说明，这3种算法对于噪声数据和孤立点数据具有一定的鲁棒性，尤其是PCRMDD算法所获得的聚类中心和划分矩阵完全相同。另一方面，表2显示了可能性聚类法和模糊聚类法在隶属度值上的差异，应用模糊聚类法，A点和B点对于每个类均得到0.5的隶属度值，因而最终影响到聚类中心的估计；而在PCM和PCRMDD中，A和B点分别得到了不同的较小的隶属度值，因此，所得的聚类中心未发生变化，各数据点的隶属度值也未出现明显改变。

表1运用4种方法所得的Bridge数据集的聚类中心和划分矩阵

表2运用4种方法所得的含噪声数据的Bridge数据集的聚类中心和划分矩阵

2)实验二：Iris数据

在该实验中，采用Iris数据集[20][21]中的二维数据向量(图3(a))来测试聚类算法的分类效果，Iris数据集由150个四维向量所构成的3个聚类集合组成，每一类包含50个数据点，并且其中有2个聚类具有较大程度的重叠。图3(b)-(e)显示了运用PCM算法、FCTMDD算法和PCRMDD算法所获得的聚类结果，其中图3(b)和(c)分别为文献[7]提出的PCM算法及其改进算法[10]的硬分类结果。从图3中可以看到，这4种方法得到的硬分类结果和聚类中心各不相同。对于具有较大重叠的两个类别，PCM算法均产生了近于相同的聚类中心，这表明采用改进的目标函数[10]仍不能处理重叠聚类问题，而早期的PCM算法[7]的分类结果发生了较大错误。与此相反，FCTMDD和PCRMDD算法获得了较为精确的分类结果和聚类中心，而PCRMDD算法则能更准确地识别出聚类中心。表3显示了运用这4种方法所得的Iris数据集的聚类有效性测度指标的计算结果，其中PCM[7]和PCM[10]分别代表文献[7]提出的PCM算法及其改进算法[10]。由该表可见，PCM算法的划分系数PC接近于1/2，表明该算法无法准确识别出Iris数据集所具有的聚类集合，并且较大的XB指数值进一步证实其所得聚类之间的重叠程度很高；而FCTMDD和PCRMDD算法得到的分类结果则能正确地反映数据集的聚类结构，并根据FS和XB指数值可得出PCRMDD算法的性能最优，其获得的聚类致密且分离性最好。

表3Iris数据集的聚类有效性测度计算结果

2.实际遥感影像数据实验

本实验采用上海市中心城区2002年11月11日Landsat7ETM+影像的第1～5和7波段数据。在进行混合像元分解之前对遥感影像进行了辐射校正和几何校正。首先对遥感数据进行辐射定标，将各波段像元值(DN值)转化为反射率值。进而，以上海地方坐标系为基准，采用全球横轴墨卡托投影(UTM)，对其进行几何精校正和投影变换，并运用最邻近法将反射率图像重采样至30m空间分辨率，经上述方法配准的图像，其均方根误差(RMSE)被控制在0.5个像元以内。根据研究区地表覆盖的主要特征、遥感影像的空间分辨率水平及动态监测的需要，将中心城区土地利用类型分为高反照度建筑(H)、低反照度建筑(L)、水体(W)、植被(R)和土壤(S)5大类，这5种地类在经校正后的4、3、2波段合成影像上的判读依据如表4所示。

表4经校正后的4、3、2波段合成影像上研究区域用地类型的判读依据

运用本发明提出的PCRMDD算法，我们对研究区反射率影像进行混合像元分解，得到了各个端元组分的空间分布图(图4)以及聚类中心向量。根据各个聚类中心的反射率值，对每一类选取在反射率特征空间中最邻近聚类中心的50个像元作为影像光谱端元，对选取的像元计算其在各个波段上的反射率均值，从而绘制成端元反射光谱曲线(图5)。

为评价研究区域遥感影像混合像元分解的准确性，首先采用均方根(RMS)误差对总体精度进行统计运算；此外，利用高分辨率航空影像并结合实地调查，通过选取一定数量样本区，对解混结果进行检验。图6为运用PCRMDD算法对反射率影像进行混合像元分解后的RMS误差图像，其平均RMS值为0.0114，最大RMS值为0.1093，对照RMS的频率直方图发现绝大部分RMS都小于0.03，其中小于0.025的像元占到总数的96％，可见PCRMDD算法解混的总体精度较高。城区中大面积水域、绿地等地表性质单一的区域均方根误差都较小，RMS误差较大的像元则基本上分布于反射率非常高的区域，这些像元对应于一些采用玻璃、金属及其他特殊装饰材料的建筑物表面，其反射率要远远高于其它高反照度材料。

进而，我们在配准后的上海市2002年高分辨率航空影像上均匀提取150个样区并结合实地调查对端元盖度进行统计检验。为保证样本检验质量，选取的矩形样区基本上要涵盖2种以上地物类型，并且其面积为150m×150m。对这些样区内的地物端元进行手工数字化，并参照实地调查结果来获取每个样区内各类端元实际的面积比率，再将这些样区和各类端元盖度图像进行叠置，统计每一样区内各类端元盖度值，通过比较发现，应用PCRMDD算法得到的端元盖度值与实际面积比率的平均误差为0.0663，其中最高为高反照度建筑占主体的样区，最低为水面占主体的样区。以高反照度建筑为例，图7中给出了地物端元的估算结果与实际面积比率之间的差异情况以反映PCRMDD算法解混的误差特征，如图7所示，对于包含高反照度建筑的83个样本区，解混误差绝大部分在±0.1以内，并且误差的绝对值随着高反照度建筑盖度的增大呈现出减小的趋势。因而，尽管在航空影像解译过程中存在着判读误差及数字化误差，但是根据上述使用航空影像选择样本区的检验结果表明，运用本发明提出的PCRMDD算法所获得的混合像元分类结果与上海市中心城区地表覆盖的实际分布十分吻合。

参考文献：

[1]朱述龙,张占睦.遥感图像获取与分析.北京:科学出版社,2000.

[2]F.J.Garcia-Haro,M.A.Gilabert,J.Melia.Linear spectral mixture modelling toestimate vegetation amount from optical spectral data.International Journal of RemoteSensing,1996,17:3373-3400.

[3]J.C.Bezdek.Pattern recognition with fuzzy objective function algorithms.NewYork:Plenum Press,1981.

[4]G.M.Foody.Hard and soft classifications by a neural network with anonexhaustively defined set of classes.International Journal of Remote Sensing,2002,23:3853-3864.

[5]C.F.Lin,S.D.Wang.Fuzzy support vector machines.IEEE Transactions onNeural Networks,2002,13(2):464-471.

[6]C.C.Borel,S.A.W.Gerstl.Nonlinear spectral mixing models for vegetative andsoil surfaces.Remote Sensing of Environment,1994,47:403-416.

[7]R.Krishnapuram,J.M.Keller.A possibilistic approach to clustering.IEEETransactions on Fuzzy Systems,1993,1(2):98-110.

[8]D.Dubois,H.Prade.Possibility theory:An approach to computerized processingof uncertainty.New York:Plenum Press,1988.

[9]G.Shafer.A mathematical theory of evidence.Princeton,NJ:Princeton UniversityPress,1976.

[10]R.Krishnapuram,J.M.Keller.The possibilistic c-means algorithm:Insights andrecommendations.IEEE Transactions on Fuzzy Systems,1996,4(3):385-393.

[11]H.Timm,R.Kruse.A modification to improve possibilistic fuzzy cluster analysis.In:The2002IEEE International Conference on Fuzzy Systems,Honolulu,Hawaii,2002:1460-1465.

[12]M.Yang,K.Wu.Unsupervised possibilistic clustering.Pattern Recognition,2006,39:5-21.

[13]S.L.Chiu.Fuzzy model identification based on cluster estimation.Journal ofIntelligent and Fuzzy Systems,1994,2:267-278.

[14]K.S.Fu.Syntactic pattern recognition and applications.San Diego,CA:Academic Press,1982.

[15]R.Krishnapuram,A.Joshi,L.Yi.A fuzzy relative of the k-medoids algorithmwith application to web document and snippet clustering.In:Proceedings of IEEEInternational Conference on Fuzzy Systems-FUZZ IEEE(99),Seoul,Korea,1999:1281-1286.

[16]J.C.Bezdek.Cluster validity with fuzzy sets.Journal of Cybernetics,1974,3(3):58-73.

[17]J.C.Bezdek.Mathematical models for systematics and taxonomy.In:Proceedings of8th International Conference on Numerical Taxonomy,Freeman,SanFrancisco,1975:143-166.

[18]Y.Fukuyama,M.Sugeno.A new method of choosing the number of clusters forthe fuzzy c-means method.In:Proceedings of the5th Fuzzy Systems Symposium,Kobe,Japan,1989:247-250.

[19]X.L.Xie,G.Beni.A validity measure for fuzzy clustering.IEEE Transactionson Pattern Analysis and Machine Intelligence,1991,13(8):841-847.

[20]E.Anderson.The irises of the Gaspe Peninsula.Bulletin of the American IrisSociety,1935,59:2-5.

[21]J.C.Bezdek,J.M.Keller,R.Krishnapuram,L.I.Kuncheva,N.R.Pal.Will thereal iris data please stand up？IEEE Transactions on Fuzzy Systems,1999,7(3):368-369。

Claims (2)

1.一种基于可能性中心点聚类的多光谱遥感影像混合像元分解方法，其特征在于，根据可能性理论和中心点聚类方法的基本原理，构建可能性中心点聚类方法的代价函数式和可能性隶属度函数，运用启发式算法对离散的解空间进行搜索，通过迭代运算更新聚类中心及可能性隶属度矩阵，将代价函数达到最小值时的聚类中心和隶属度归一化结果作为光谱端元和组分比例，即各类地表覆被的盖度；具体步骤包括：

1)确定聚类个数C、模糊加权指数m和最大迭代次数max_iter，运用减法聚类法从数据集X中选择C个初始聚类中心V＝{v₁,v₂,…,v_c}，将迭代次数iter设为0；

2)对所有数据对象计算(i＝1,…,C；j＝1,…,n)，并对h_j进行升序排列构成h_j:n，选取对应于前s个h_j:n的s个数据对象，其中，h_j:n代表将h_j进行升序排列后得到的第j项；η_i为运行模糊C中心点算法收敛后所得到的估计值：

${\mathrm{\&eta;}}_i=\frac{\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{\left|\right|x_j-v_i\left|\right|}^2}{\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m}$

或

${\mathrm{\&eta;}}_i=\frac{\underset{u_{\mathrm{ij}}>\mathrm{\&alpha;}}{\mathrm{\&Sigma;}}{{\left|\right|x}_j-v_j\left|\right|}^2}{\underset{u_{\mathrm{ij}}>\mathrm{\&alpha;}}{\mathrm{\&Sigma;}}1}$

其中α为合适的阈值，在[0.1,0.4]区间内取值；

3)对这s个数据对象计算其隶属度值(i＝1,…,C；j＝1,…,s)；

4)储存当前聚类中心V_current＝V；

5)计算新的聚类中心v_i＝x_q(i＝1,…,C)，

其中， $q=\underset{1\&le;k\&le;s}{\arg \min }[\underset{j=1}{\overset{s}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}:n}^m{\left|\right|x_{j:n}-x_{k:n}\left|\right|}^2+\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&eta;}}_i}{m}(u_{\mathrm{ij}:n}^m\ln u_{\mathrm{ij}:n}^m-u_{\mathrm{ij}:n}^m)];$

6)将iter增1，即iter＝iter+1；如果聚类中心不再发生变化，或者iter＝max_iter，那么终止迭代过程，即可获得可能性隶属度矩阵U与聚类中心V；否则，转向第(2)步；

7)对算法收敛后得到的可能性隶属度u_ij:n进行归一化处理，并将噪声数据的隶属度值设为0，其计算公式为

$u_{\mathrm{ij}}^{\&prime;}=\left\{\begin{array}{cc}\frac{u_{\mathrm{ij}:n}}{\underset{i=1}{\overset{C}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}:n}}& j=1,...,s\\ 0& j=s+1,...,n\end{array}\right.$

2.按权利要求1所述的多光谱遥感影像混合像元分解方法，其特征在于，在迭代运算中，对应于前s个h_j:n的s个数据对象，通过取s＝80％×n，求取目标函数的最小值来获得算法的最优解。