CN108647719A - 用于大数据量光谱遥感图像分类的非监督聚类方法 - Google Patents

Abstract

本发明公开了一种用于大数据量光谱遥感图像分类的非监督聚类方法。将原始数据划分为若干数据块，通过峰值密度搜索方法得到各数据子块的聚类中心；将聚类中心重新分成若干数据块，通过峰值密度搜索方法再次聚类，减少聚类中心数；重复分块‑聚类过程至可以用一个二维矩阵表征任意两个聚类中心的相似度，进而得到最终的分类结果。本发明方法的优势在于：适用性好，不仅可以用于谱段数较多的高光谱遥感图像分类，也适用于谱段数较少的多光谱遥感图像或谱段选择后的高光谱遥感图像分类；运算效率较高，分块处理减少了相似度矩阵的计算冗余，且由于各数据块的聚类处理相互独立，可采用并行处理加快分类速率。

Description

用于大数据量光谱遥感图像分类的非监督聚类方法

技术领域

本发明涉及光谱遥感图像地物分类的技术领域，具体涉及一种用于大数据量光谱遥感图像分类的非监督聚类方法。

背景技术

高光谱与多光谱遥感图像记录了同一区域在不同观测谱段内的辐射特性。由于植被、土壤、建筑物、水体等各类地表物体的光谱辐射特性存在显著差异，通过分析光谱遥感数据，能够获取不同类型地物的空间分布信息及光谱辐射特征信息。这些分类结果在地表植被分布研究、土壤与地质勘探、城市覆盖物调查及水体质量监测等领域有重要应用。

非监督分类是一种在不预先知晓地物种类数、地物光谱特性等先验信息的情况下，对图像进行分类的方法。方法中主要包含选择分类特征，确定分类依据，制定分类准则并实现算法三个部分。选取两个谱段数据实现光谱遥感图像分类的过程如图1所示。

为了便于分析，将二维图像数据转换为由像元组成的一维列向量，则原始的三维光谱遥感图像转换为二维光谱遥感数据。假设原始的光谱遥感图像N个像元，L个谱段，则原始图像数据Y可以表示为：

Y_N×L＝{y₁,y₂,...,y_L}_N×L＝{x₁,x₂,...,x_N}^T (1)

式中，y_i(i＝1,2,…,L)为第i谱段的图像数据，为N×1的列向量；x_j(j＝1,2,…,N)为第j像元的光谱数据，为L×1的列向量。

选择分类特征是从原始图像数据Y中选取所需谱段的数据y_i建立分类特征空间的过程，又可称为谱段选择(参见文献[1]C.Chang,Q.Du,T.Sun,and M.L.G.Althouse.AJoint Band Prioritization and Band-Decorrelation Approach to Band Selectionfor Hyperspectral Image Classification[J].IEEE Transactions on Geoscience andRemote Sensing,1999,37(6):2631-2641.,文献[2]C.Chang and S.Wang.ConstrainedBand Selection for Hyperspectral Imagery[J].IEEE Transactions on Geoscienceand Remote Sensing,2006,44(6):1575-1585.及文献[3]Q.Du and H.Yang.Similarity-Based Unsupervised Band Selection for Hyperspectral Image Analysis[J].IEEEGeoscience and Remote Sensing Letters,2008,5(4):564-568.)。由于仅选取了两个谱段，故图1中构建的特征空间为二维平面，各个像元在特征空间的坐标取决于像元在两个谱段的辐射强度。若分类中使用的谱段数多于三个，则需构造高维特征空间(参见文献[4]王永革，藤岩梅，贾超华，冯伟杰.应用泛函分析[M].北京：北京航空航天大学出版社.2012.)。

分类依据用于判断仍以两个像元之间的相似度，定义为：

d_ij＝||x_i-x_j|| (2)

式中，d_ij为像元i与像元j之间的相似度；‖·‖为范数符号，可以采用不同的计算方法，如相似系数、距离空间及相对熵等(参见文献[5]童庆禧，张兵，郑兰芬.高光谱遥感——原理、技术与应用[M].北京：高等教育出版社.2006.)。计算任意两个像元相似度后，可以生成相似度矩阵，即：

式中，D为N×N的二维相似度矩阵。

分类准则通常与使用的分类算法有关，最核心的问题为确认图像中地物的种类数及各类地物的光谱特性。理论上而言，遥感图像中的同类地物在相同的观测条件下应该具有相同的光谱特性。然而，由于各个像元的地表地形、大气环境、光照条件等因素存在差异，且受传感器的空间分辨率限制，单一像元中往往包含多类地物，因此遥感图像中的同类地物表现出不同光谱特性，而不同地物之间的差别变得模糊。例如图1所示的分类特征空间中，任意两类地物的分类边界出现不确定像元(六边形所示)，分类结果产生误差。

对于非监督分类方法而言，确定地物的种类数及各类地物的光谱特性需要利用式(3)中定义的相似度矩阵。然而，对于大数据量遥感图像，由于像元数很多，生成二维相似度矩阵容易超出计算机的物理内存。此外，由于像元i至像元j的距离与像元j至像元i的距离相等，计算所得的相似度矩阵包含冗余数据。

与本发明相关的现有技术介绍如下：

1.1现有技术一的技术方案

正交子空间投影方法既可以用于选择分类特征(参见文献[3])也可以用于高光谱遥感图像分类(参见文献[6]C.-I Chang.Orthogonal Subspace Projection(OSP)Revisited:A Comprehensive Study and Analysis[J].IEEE Transactions onGeoscience and Remote Sensing,2005,43(3):502-518.)。正交子空间投影方法假定能够找到一组像元，它们光谱特性的线性组合能够表示其他任意像元的光谱特性。通过前向搜索方法寻找一组最大线性无关的像元组合，此时选取的像元数即为地物种类数，选取像元的光谱特性即为各类地物的光谱特性。由于该方案不需要计算式(3)中的相似度矩阵，故采用该方案对光谱遥感图像进行分类不会超出计算机的物理内存限制。

1.2现有技术一的缺点

根据线性泛函理论(参见文献[4])，由于构成高维空间的向量数应当等于维度数，故选择的像元个数不应大于谱段数。参考文献[2]中的研究结论表明，选取高光谱遥感图像中少数几个谱段的数据就能够取得较准确的分类结果。由于正交子空间投影方法可区分的最大地物种类数取决于选取的谱段数，该方案不适用于谱段选择后的遥感图像分类。

2.1现有技术的技术方案二

监督分类方法常用于高光谱遥感图像处理(参见文献[7]F.Melgani andL.Bruzzone.Classification of Hyperspectral Remote Sensing Images with SupportVector Machines[J].IEEE Transactions on Geoscience and Remote Sensing,2004,42(8):1778-1790.和文献[8]F.Ratle,G.C.-Valls,and J.Weston.Semisupervised NeuralNetworks for Efficient Hyperspectral Image Classification[J].IEEETransactions on Geoscience and Remote Sensing,2010,48(5):2271-2282.)。该方法利用实地勘探或人工鉴别等方法获取先验信息，以此确定遥感图像中地物的种类数及各类地物的典型光谱特性。进而根据设定的分类依据，计算实测数据与各类先验地物的相似度，并确定实测数据的类别。由于不需要从实测遥感图像中获取地物种类数与各类地物的光谱特性，因此不需要生成式(3)中的相似度矩阵。

2.2现有技术二的缺点

星载或机载传感器获取的光谱遥感图像，其条带通常覆盖较广的地理区域。通过实地勘探获取某一观测区域的先验信息需要花费大量的人力、物力。而且在其他观测区域中，很可能存在不同于先验信息的地物类别。因此，采用监督分类方案处理光谱遥感图像，一些像元的类别归属无法确定，分类结果不完整。

3.1现有技术三的技术方案

聚类方法是一类常用的数据分类方法，其基本思路是将相似度接近的两类数据归为一类，实现方法包括K均值算法、模糊聚类方法、子空间聚类方法等(参见文献[9]A.K.Jain.Data Clustering:50Years beyond K-means[J].Pattern RecognitionLetters,2010,31:651-666.)。采用不同的聚类方法处理光谱遥感图像，得到的分类结果大体相同。此外，聚类方法对谱段数不敏感，适用于任意谱段数的光谱遥感图像分类。

3.2现有技术三的缺点

聚类方法需要计算任意两个像元之间的相似度，利用式(3)生成相似度矩阵。对于大数据量遥感图像，由于像元数很多，生成相似度矩阵容易超出计算机的物理内存。此外，由于相似度矩阵存在冗余数据，该方案的计算效率低。

发明内容

本发明所要解决的技术问题为：针对以上三种分类方案各自的缺点，为解决大数据量遥感图像的地物分类问题，本发明提出用于大数据量光谱遥感图像分类的非监督聚类方法。先将原始光谱遥感数据分为许多数据块，采用密度峰值搜索方法寻找各个数据块的聚类中心。接着，采用密度峰值搜索方法对将各数据块选取的聚类中心进行再次聚类，从而得到最终的聚类结果。若第一次聚类后的聚类中心数仍然过多，则采用相似思路将聚类中心分为多组并分别聚类，从而减少聚类中心数。重复“分块-聚类”过程，直到可以采用一个二维矩阵描述任意两个聚类中心的相似度。该方案适用于任意谱段数的光谱遥感图像分类，且分类过程中不需要知晓地物的先验信息。

本发明采用的技术方案为：一种用于大数据量光谱遥感图像分类的非监督聚类方法，包括如下步骤：

步骤1：谱段选择；

步骤2：按像元数n划分数据块；

步骤3：搜寻各数据块的聚类中心；

步骤4：判断聚类中心数是否大于步骤2中设定的像元数n；

步骤5：增加像元数n；

步骤6：选取局域密度较高的聚类中心；

步骤7：光谱遥感数据按聚类中心分类。

与现有光谱遥感图像分类方法相比，本发明提出的基于数据分块的光谱遥感图像非监督聚类方法具有以下优点：

(1)本发明方法具有普适性。所提出的分类方法适用于大数据量光谱遥感图像的地物非监督分类处理，不仅可以用于谱段数较多的高光谱遥感图像分类，也适用于谱段数较少的多光谱遥感图像或谱段选择后的高光谱遥感图像分类；

(2)本发明方法具有较高的运算效率。分块处理减少了相似度矩阵的计算冗余，且由于各数据块的聚类处理相互独立，可采用并行处理加快分类速率。

附图说明

图1为光谱遥感图像分类过程示意图；

图2为基于图像数据分块的非监督聚类流程图；

图3为高光谱遥感数据分块聚类示意图；

图4为华盛顿DC高光谱遥感数据分类结果，其中，图4(a)为实测场景的伪彩图，图4(b)为全谱段数据分类结果图，图4(c)为谱段选择后数据分类结果图；

图5为MODIS的A2014001号数据谱段1至谱段7分类结果，其中，图5(a)为平均反射率图像，图5(b)为全谱段分类结果图；

具体实施方式

下面结合附图以及具体实施方式进一步说明本发明。

本发明的技术方案框图如图2所示，基本技术原理如下。

步骤1：谱段选择。

高光谱遥感图像包含上百个连续谱段，数据量很大。在图像处理过程中，常采用谱段选择方法选取最优的特征谱段，牺牲一些分类精度以大幅提升分类处理的效率。各类谱段选择方法的复杂度不同，这里简要介绍文献[1]中的主成分分析方法。

如式(1)所示，假设原始图像数据Y包含N个像元，L个谱段，定义该数据的相关矩阵为：

式中，R为L×L的相关矩阵；x_i为第i像元的光谱数据，为L×1的列向量；m为所有像元光谱数据的均值，为L×1列向量。

由于相关矩阵R正定，故对其进行特征值分解可得特性值向量{λ₁,λ₂,…,λ_L}及特征向量矩阵v，即：

定义特征因子为：

式中，ρ_l表示第l谱段的数据方差。

假定选取了m个谱段，则这m个谱段图像所占的功率比率为：

利用式(7)可以计算得到各谱段图像所占的功率比率。谱段选择过程中认为功率比率越高的谱段优先级越高，故式(7)定义了各谱段的优先级。

主成分分析认为同一像元、不同谱段的数据之间应当具有相关性，可以采用互熵(又称为Kullback-Leibler距离)衡量两个谱段的相关度，即：

D(p,q)＝L(p；q)+L(q；p) (8)

式中，D为两幅图像的互熵；p与q为两组数据，谱段选择中定义为两个谱段的辐射图像；L为相对熵，有计算式：

主成分分析方法先通过式(6)与式(7)确定各谱段的优先级，然后利用式(8)计算优先级低的谱段与优先级高的谱段的相关度，若计算所得的互熵大于预设门限，则认为该谱段包含重要光谱特性信息并纳入被选谱段。因此，主成分分析方法是一种“排序-去冗余”的前向谱段选择方法。

与传统的主成分分析方法不同，本发明中采用聚类后各聚类中心的光谱辐射廓线作为输入数据进行谱段选择，而非高光谱实测图像。由于聚类中心数目远远小于像元数，故采用廓线进行主成分分析能够提升谱段选择的效率。

应当指出，不论是高光谱遥感数据还是多光谱遥感数据均可以采用所有谱段数据进行地物分类处理。故提出的分类方案，应该既适用于全谱段的图像分类处理，又适用于分类特征选择后的图像分类处理。

步骤2：按像元数n划分数据块。

图3所示为高光谱遥感数据三次聚类示意图。如图3所示，经过维度变换后，原始的三维高光谱遥感图像转变为二维高光谱遥感数据，高光谱数据失去了图像的纹理信息。像元数n的初始值为主观设定。若n太小，则所需的聚类次数很多，需要较长的分类时间。若n太大，则相似度矩阵中存在大量的冗余数据。

步骤3：搜寻各数据块的聚类中心。

非监督聚类方法有许多种，这里采用文献[10](A.Rodriguez andA.Laio.Clustering by Fast Search and Find of Density Peaks[J].Science,2014,344(6191):1492-1496.)中的密度峰值搜索方法对光谱遥感数据进行分类。该聚类方法已用于数据量较少时的分类处理，如二维数据点分类、航班线路规划、高光谱谱段选择等，但尚未用于数据量较大时的分类处理或是光谱遥感图像地物分类处理。再者，该聚类方法尚未推广至数据分块分类或分级分类的情形。

定义数据块中某一像元的局域密度为：

式中，ρ_i为像元i的局域密度；d_ij为像元i与像元j之间的相似度，由式(2)定义；d_c为相似度门限，其大小将影响分类结果；χ为判决函数，有定义式：

如图1所示，由式(11)与式(12)可知，某一像元i的局域密度是指在分类特征空间内，以像元i为中心、相似度门限d_c为半径围成的区域内包含的像元个数。

采用欧式距离定义任意两像元i与j之间的相似度为：

式中，x_ik与x_jk分别为像元i与像元j的光谱数据；k为谱段编号。由式(13)可知，两个像元之间的欧氏距离越小，相似度越高。

在各数据块中，利用式(11)计算所有像元的局域密度，并从大到小排序。定义某像元的欧式距离最小值为：

δ_i＝min(d_ij),j∈ρ_j＞ρ_i (14)

式中，δ_i为像元i至高于其局域密度像元的距离最小值；min为求最小值符号。由式(14)可知，像元的距离最小值代表了相似度最大值。

式(11)与式(14)提供了判断聚类中心的两个标准。一方面，在分类特征空间中，聚类中心应该具有高的局域密度，即聚类中心附近的像元很多；另一方面，新聚类中心应该与已有聚类中心的相似度低，即新聚类中心与已有聚类中心的光谱特性不同。因此，聚类中心的局域密度ρ应高于阈值，且距离最小值δ亦高于阈值。

考虑到分块方法破坏了光谱遥感数据的完整性，故在搜寻各数据块的聚类中心时，仅采用距离阈值作为聚类中心的判断标准。局域密度标准仅用于最终聚类中心的选取，将在步骤5中介绍。

步骤4：判断聚类中心数是否大于步骤2中设定的像元数n。

按步骤3进行处理，得到各数据块的聚类中心后，需要进一步确定光谱遥感图像的聚类中心。对于非监督分类方法，由于不知晓各数据块聚类结果的相关性，故使用密度峰值搜索方法对所得的聚类中心进行再次聚类，进而获取最终的聚类中心。此外，考虑到像元数n的设置较为主观，若第一次聚类所得的聚类中心数过多，无法直接计算相似度矩阵，则将所得的聚类中心划分为若干数据块进行再次聚类。重复该过程，直至聚类中心数满足要求。

对于再次聚类，类似于式(11)，可以定义新的局域密度为：

式中，ρ_i′为再次聚类的局域密度；d_ij为聚类中心i与聚类中心j之间的相似度；d_c为相似度门限；χ′为判决函数，有定义式：

式中，ρ_j为前一次聚类所得的局域密度。

由式(15)与式(16)可知，再次聚类的局域密度是前一次聚类的局域密度之和。因此，对于相同的光谱遥感数据，不论如何设定像元数n、进行多少次聚类，采用该方法得到的聚类中心理应相同。

步骤5：增加像元数n。

所需的聚类次数取决于像元数n。适当增加像元数能够提高所提出方法的分类效率。具体操作时，设定每次聚类的像元数增加1.2倍

步骤6：选取局域密度较高的聚类中心。

经过步骤2至步骤5的处理，可以认为最终聚类中心的局域密度表示了某类地物在光谱遥感图像中的丰富程度。尽管，理论上而言，遥感图像中的同类地物在相同的观测条件下应该具有相同的光谱特性，然而受观测条件和传感器系统噪声的影响，实测光谱遥感图像中往往存在奇异点。由于地物分类处理更加关注成分丰富的地物类别，而对异常像元的关注度较低。具体操作时，认为某类地物的像元数应该大于遥感图像总像元数的千分之一。

步骤7：光谱遥感数据按聚类中心分类。

经过步骤2至步骤6处理得到光谱遥感图像的聚类中心后，需要对各像元的类别归属进行判决。判决过程可以采用分类器、分类判别算法或是神经网络^[5]。为提高分类效率，采用最小距离分类判断各像元的类别，即认为像元类别与距离最近的聚类中心相同。

假定遥感数据包含N个像元，最终的聚类中心有n个。与式(14)类似可定义最小距离：

δ_i′＝min(d_ij),1≤i≤N,1≤j≤n (17)

式中，δ_i′为像元i至聚类中心的最小距离；d_ij为距离函数，采用式(13)定义的欧式距离计算。

受传感器空间分辨率的限制，光谱遥感图像中有些像元的光谱特性是多类地物的混合，最小距离分类判断将混合像元划归为混合比例较大的类别。

具体实施举例：

下面采用实测光谱遥感数据对本发明做进一步说明。使用的数据为HYDICE高光谱传感器测量得到的华盛顿DC Mall数据图像(参见文[11]D.Landgrebe.HyperspectralImage Data Analysis[J].IEEE Signal Processing Magazine,2002,1:17-28.)及MODIS传感器谱段1至谱段7的多光谱实测数据(参见文献[12]C.O.Justice,E.Vermote,J.R.G.Townshend,et al.The Moderate Resolution Imaging Spectroradiometer(MODIS):Land Remote Sensing for Global Change Research[J].IEEE Transactionson Geoscience and Remote Sensing,1998,36(4):1228-1249.)。

(1)采用发明方法对HYDICE高光谱传感器数据进行验证。

步骤1：利用式(4)至式(10)进行谱段选择。

HYDICE高光谱传感器原始测量数据大小为1280×307×210，即一副辐射图像包含1280×307个像元，全部数据包含210个谱段。由于某些谱段内的大气吸收效应很强，测量数据的信噪比很低，故全谱段数据选取了191个信噪比较高的谱段数据用于地物分类处理。对于该数据而言，主成分分析方法选取的谱段为：谱段40(579nm)，谱段104(1435nm)及谱段116(1595nm)。

步骤4：判断聚类中心数是否大于步骤2中设定的像元数n。

设定数据分块的初始像元数n为1200，随后每次聚类的像元数增加1.2倍。由于图像数据为量化后的灰度值，设定距离门限d_c为2000。在确定最终聚类中心时认为某类地物的像元数应该大于总像元数的万分之一，即39个像元。

步骤7：光谱遥感数据按聚类中心分类。

图4所示为采用本发明方法得到的数据分类结果图。图4(a)为由谱段17(459nm)、谱段27(503nm)及谱段60(758nm)合成的伪彩图；图4(b)为采用全谱段数据得到的分类结果；图4(c)为主成分分析方法选取谱段40(579nm)，谱段104(1435nm)及谱段116(1595nm)三个谱段数据得到的分类结果。图4(a)的亮度代表辐射强度，容易辨别图片中道路两旁的树木，而图4(b)与图4(c)的亮度代表地物类别编号，受混合像元影响，树木的形状变得模糊。由图可知，不论是采用全谱段还是采用谱段选择后的数据，分类结果中均可以明显区分水域、植被、道路和建筑物等主要地物类别，并且能够区分不同种类的建筑物材料。

对比图4(b)与图4(c)可知，高光谱遥感数据的谱段相关度很高，数据存在大量冗余，故采用几个谱段数据即能够区分遥感数图像中的主要地物。然而，地物分类的精细程度及准确度均与选用的谱段数有关。理论上而言，选取的谱段数目越多，分类的精细程度越高。另外，降低聚类过程中的局域密度门限，增加最终的聚类中心数，也能够提升分类的精细程度。

(2)为了验证本发明的通用性，采用MODIS传感器实测数据进行验证。

步骤1：利用式(4)至式(10)进行谱段选择。

MODIS的Level 3数据产品MOD09A1中提供了八天的平均地表反射率。验证过程中选取了A2014001号数据，该数据包含2400×2400个像元。由于仅包含7个谱段的数据，故不作谱段选择处理。

步骤4：判断聚类中心数是否大于步骤2中设定的像元数n。

数据分块的初始像元数n仍为1200，距离门限d_c的设置与传感器信噪比有关，取值为2.4×10^-4。

步骤7：光谱遥感数据按聚类中心分类。

图5所示为采用本发明方法得到的数据分类结果图。图5(a)为该数据的平均反射率图像；图5(b)为采用本发明方法得到的全谱段分类图。由图可知，聚类方法可以区分的地物种类数可以多于分类中采用的谱段数，且分类图中可以区分平均反射率相近但光谱特性不同的像元。此外，由分类图5(b)可知，场景中云层的光谱特性存在较大差异，故云层对地物分类结果有较大影响。

本发明还可以采用的替代方案同样能完成发明目的：分类特征选择(步骤1：谱段选择)可采用除主成分分析以外的其他选择方法替代，例如约束能量优化方法、正交子空间投影方法等；像元类别判决过程(步骤7)可以采用除最小距离以外的其他类别判决方法替代，如最大似然判别法、模糊分类等。

Claims (2)

1.一种用于大数据量光谱遥感图像分类的非监督聚类方法，其特征在于：包括如下步骤：

步骤1：谱段选择；

步骤2：按像元数n划分数据块；

步骤3：搜寻各数据块的聚类中心；

步骤4：判断聚类中心数是否大于步骤2中设定的像元数；

步骤5：增加像元数n；

步骤6：选取局域密度较高的聚类中心；

步骤7：光谱遥感数据按聚类中心分类。

2.根据权利要求1所述的一种用于大数据量光谱遥感图像分类的非监督聚类方法，其特征在于：步骤2至步骤6描述的基于密度峰值搜索方法的数据分块-聚类方法，其核心过程如下：

对于第一次分块-聚类过程而言，定义数据块中某一像元的局域密度为：

式中，ρ_i为像元i的局域密度；d_ij为像元i与像元j之间的相似度；d_c为相似度门限，其大小将影响分类结果；χ为判决函数，有定义式：

采用欧式距离定义任意式(11)中两像元i与j之间的相似度为：

式中，x_ik与x_jk分别为像元i与像元j的光谱数据；k为谱段编号；

在各数据块中，利用式(11)计算所有像元的局域密度，并从大到小排序，定义某像元的欧式距离最小值为：

δ_i＝min(d_ij),j∈ρ_j＞ρ_i (14)

式中，δ_i为像元i至高于其局域密度像元的距离最小值；min为求最小值符号，由式(14)可知，像元的距离最小值代表了相似度最大值；

在搜寻各数据块的聚类中心时，仅采用距离阈值d_c作为聚类中心的判断标准，局域密度标准仅用于选取最终的聚类中心；

对于第二次及后续分块-聚类过程而言，将前一次聚类所得聚类中心划分为若干数据库并聚类，定义新的局域密度为：

式中，ρ_i′为再次聚类的局域密度；χ′为判决函数，有定义式：

式中，ρ_j为前一次聚类所得的局域密度；

对于最后一次聚类过程：通过计算聚类中心之间的相似度，合并距离小于距离阈值d_c的聚类中心；剔除局域密度较小的聚类中心，从而得到最终的聚类结果。