CN110070035B - 一种多光谱遥感影像分类方法和系统 - Google Patents

Abstract

一种多光谱遥感影像分类方法和系统。该分类方法采用十进制编码方法对多波段进行压缩，将多个波段的灰度值转换成一个十进制编码值，且编码的位数和波段数相同，采用最小距离法或光谱夹角法对所有十进制编码进行聚类分类，得到聚类中心的参考十进制编码，并根据最小距离法或光谱夹角法计算所有十进制编码和参考十进制编码对应的类别属性，最后，通过索引的方式实现整个多光谱遥感影像的自动分类。本发明的多光谱遥感影像分类方法处理相比于现有的分类方法具有精度更高、处理速度更快的特点。

Description

一种多光谱遥感影像分类方法和系统

技术领域

本发明涉及图像处理技术领域，具体涉及一种多光谱遥感影像分类方法和系统。

背景技术

遥感影像分类方法分为两大类：非监督分类和监督分类。

监督分类是需要学习训练的分类方法，首先从研究区域选取有代表性的训练场地作为样本，根据已知训练区提供的样本，选择特征参数，建立判别函数，依据样本类别的特征来识别非样本像元的归属类别。具体方法有最大似然法，人工神经网络法等。

非监督分类不需要人工采样地面样本点数据，在没有先验类别作为样本的条件下，根据像元间相似度大小进行计算机自动判别，通过聚类的方法实现自动分类，分类后再确定地面类别。主要方法有迭代自组织的数据分析法(ISODATA法)，K均值法等。

非监督分类方法主要有两大类：

一、迭代自组织的数据分析法(ISODATA法)

迭代自组织数据分析算法(Iterative Self-Organizing Data AnalysisTechniques Algorithm,ISODATA)与K均值算法类似，即聚类中心的位置同样是通过样本均值的迭代运算决定。不同的是，这种算法在运算的过程中聚类中心数目不是固定不变的，而是反复进行修改，以得到较合理的类别数K，这种修改通过模式类的合并和分类来实现，合并和分类在一组预先选定的参数指导下进行。ISODATA的特点是计算简单，适用于识别致密聚类。

合并主要发生在某一类样本数量比较少的情况下，或者两类样本聚类中心之间距离太小的情况。为此，需要指定每一类中最少样本数和两类聚类中心之间的最小距离参数。类分类主要发生在某一类的某分量出现类内方差过大的现象时，适合将其分裂成两类，使类内方差比较合理。为此，需要指定类内某个分量方差的参数，用以决定是否需要将某一类分裂成两类。

ISODATA具体计算方法如下：

其中，设有N个模式样本，X₁，X₂，...，X_N；K为期望的聚类中心数；θ_N为每一类中至少应包含的样本数；θ_s为类内样本标准差阈值；θ_c为两个聚类中心之间的最小距离阈值；L为一次迭代中允许合并的最多对数；I为允许迭代的次数。

第一步，任选c个聚类中心z₁(I)，z₂(I)，...，z_c(I)；定义算法参数K，θ_N，θ_s，θ_c，L，I。其中c不要求等于期望的聚类中心数K。

第二步，分配N个样本，按最近邻规则分配到c个聚类中。若||x-z_i||＜||x-z_j||，i＝1，2，...，c，i≠j，则

其中，X_i表示分类到聚类中心z_i的样本子集，N_i为X_i中的样本个数。

第三步，若N_i＜θ_N，则去除X_i，使c＝c-1，即将样本数比θ_N少的样本子集删去。

第四步，修正各聚类中心z_i

第五步，计算X_i中样本与各聚类中心间的平均距离：

第六步，计算总体的平均距离

其中N为样本集中样本总数。

第七步，判断分裂、合并及迭代运算步骤

(1)若迭代已达I次，置θ_c＝0，转到第十一步，算法结束；

(2)若c≤K/2，即聚类中心小于或等于希望数的一半，转到第八步，将已有分类分裂；

(3)若迭代次数是偶数，或c≥2K，即聚类中心数目大于期望数的两倍，则转到第十一步，进行合并处理；

(4)若(2)和(3)不满足则继续，转入第八步。

第八步：计算每个聚类中样本距离的标准差向量。对第X_i类有σ_i＝[σ_i1，σ_i2，...，σ_in]^T，i＝1，2，...，c，

其中各分量为

式中，n是样本维数，即x_i是n维模式向量；x_ij是X_i类样本x_i的第j个分量，z_ij是X_i类聚类中心距离z_i的第j个分量，所以σ_ij是x_i的第i个分量的标准差。

第九步：求出σ_i中最大分量σ_imax，i＝1,2,…,c。

第十步，若σs_imax，i＝1,2,…,c，同时满足以下条件之一：

(1)

和N_i＞2(θ_N+1)，即类内平均距离大于总体平均距离，并且X_i类样本数过大；

(2)c≤K/2，即聚类数小于期望数的一半；则将X_i分成两个新的聚类中心，

和

删去z_i，并使c＝c+1.其中，

这样构成：z_i中对应σ_imax的分量加上ασ_imax；

这样构成：z_i中对应σ_imax的分量减去ασ_imax；

其中，0＜α≤1；选择α的基本要求是，使任意样本到这两个新的聚类中心

和

之间有一个足够可检测的距离差别，但又不能太大。

如果完成分裂，则迭代次数加1，I＝I+1，转到第二步；否则继续进行第十一步。

第十一步，计算全部聚类中心的两两距离d_ij；

d_ij＝||z_i-z_j||，i≠j，i，j＝1，2，...，c

第十二步，比较：如果d_ij>θ_c,转到第十四步；否则如果d_ij<θ_c,将d_ij<θ_c的值升序排列，即d_i1j1＜d_i2＜...＜d_iIjI，I<L

第十三步：从d_i1j1开始，逐对合并，算出新的聚类中心：

删z_iI和z_jI，并使c＝c-1.注意，只允许一对对合并，并且一个聚类中心只能合并一次。

第十四步，迭代处理：若是最后一次迭代，l＝I,算法结束。否则有两种情况：

(1)不修改参数，l＝l+1，转到第二步；

(2)需要人工输入参数，l＝l+1，转到第一步。

以上是ISODATA的具体计算方法步骤。

二、K均值法

K-means算法是硬聚类算法，是典型的基于原型的目标函数聚类方法的代表，它是数据点到原型的某种距离作为优化的目标函数，利用函数求极值的方法得到迭代运算的调整规则。K-means算法以欧式距离作为相似度测度，它是求对应某一初始聚类中心向量V最优分类，使得评价指标J最小。算法采用误差平方和准则函数作为聚类准则函数。

K-means算法是很典型的基于距离的聚类算法，采用距离作为相似性的评价指标，即认为两个对象的距离越近，其相似度就越大。该算法认为簇是由距离靠近的对象组成的，因此把得到紧凑且独立的簇作为最终目标。

k个初始类聚类中心点的选取对聚类结果具有较大的影响，因为在该算法第一步中是随机的选取任意k个对象作为初始聚类的中心，初始地代表一个簇。该算法在每次迭代中对数据集中剩余的每个对象，根据其与各个簇中心的距离将每个对象重新赋给最近的簇。当考察完所有数据对象后，一次迭代运算完成，新的聚类中心被计算出来。如果在一次迭代前后，J的值没有发生变化，说明算法已经收敛。

算法过程如下：

1)从N个文档随机选取K个文档作为质心；

2)对剩余的每个文档测量其到每个质心的距离，并把它归到最近的质心的类；

3)重新计算已经得到的各个类的质心；

4)迭代2～3步直至新的质心与原质心相等或小于指定阈值，计算结束。

具体如下：

输入：k,data[n]；

(1)选择k个初始中心点，例如c[0]＝data[0],…c[k-1]＝data[k-1]；

(2)对于data[0]....data[n]，分别与c[0]…c[k-1]比较，假定与c[i]差值最少，就标记为i；

(3)对于所有标记为i点，重新计算c[i]＝{所有标记为i的data[j]之和}/标记为i的个数；

(4)重复(2)(3)，直到所有c[i]值的变化小于给定阈值。

K-MEANS算法的工作原理及流程：

输入：聚类个数k，以及包含n个数据对象的数据库。

输出：满足方差最小标准的k个聚类。

处理流程：

(1)从n个数据对象任意选择k个对象作为初始聚类中心；

(2)根据每个聚类对象的均值(中心对象)，计算每个对象与这些中心对象的距离；并根据最小距离重新对相应对象进行划分；

(3)重新计算每个(有变化)聚类的均值(中心对象)；

(4)循环(2)到(3)直到每个聚类不再发生变化为止。

k-means算法接受输入量k；然后将n个数据对象划分为k个聚类以便使得所获得的聚类满足：同一聚类中的对象相似度较高；而不同聚类中的对象相似度较小。聚类相似度是利用各聚类中对象的均值所获得一个“中心对象”(引力中心)来进行计算的。

工作过程k-means算法的工作过程说明如下：

首先从n个数据对象任意选择k个对象作为初始聚类中心；而对于所剩下其它对象，则根据它们与这些聚类中心的相似度(距离)，分别将它们分配给与其最相似的(聚类中心所代表的)聚类；然后再计算每个所获新聚类的聚类中心(该聚类中所有对象的均值)；不断重复这一过程直到标准测度函数开始收敛为止。一般都采用均方差作为标准测度函数。k个聚类具有以下特点：各聚类本身尽可能的紧凑，而各聚类之间尽可能的分开。

以上是现有技术方案非监督分类方法中的ISODATA法和K均值法，主要有三个缺陷：

(1)现有方法都是从N个样本中选择k个初始类聚类中心点，而初始类聚类中心的选择直接影响分类的效果。

(2)现有方法分类精度低，效果差。

(3)现有方法针对大图像的处理速度慢。

(4)现有方法人工干预多。

针对上述缺陷，本发明提出一种多光谱遥感影像分类方法和系统，属于一种非监督分类方法。

发明内容

本发明针对现有非监督分类方法人工干预多、分类精度低、处理速度慢等现状，参考人眼识别地物的基本原理，基于十进制编码压缩和多光谱遥感影像索引，提出一种多光谱遥感影像的快速分类方法和系统，并和传统的ISODATA法、K-mean法分类结果进行定量比较，得到多光谱遥感影像的快速高精度自动分类方法和系统。

本发明的第一方面提供了一种多光谱遥感影像分类方法，包括如下步骤：

读取多光谱遥感影像，获取所述影像中每个像素在预定波段的灰度值；

分别统计所述多光谱遥感影像的所有像素在各预定波段的灰度直方图，获取每个灰度值对应的直方图分布频率；

对直方图分布频率大于一预定分布频率的灰度值进行十进制编码压缩处理，建立所述灰度值在预定波段对应的十进制编码；

对所述十进制编码进行聚类分析，获得参考十进制编码；

为所有十进制编码和参考十进制编码分配分类赋值，建立索引查找表；

根据所述索引查找表对所述多光谱遥感影像的灰度值进行分类并保存分类结果。

进一步的，所述对直方图分布频率大于一预定频率分布的灰度值进行十进制编码压缩处理，建立所述灰度值在预定波段对应的十进制编码的步骤包括：

获取压缩系数；

根据所述压缩系数对各预定波段的灰度值进行十进制编码压缩，得到被压缩的灰度值；

对被压缩的灰度值建立所述预定波段对应的十进制编码。

进一步的，按照下式获取压缩系数Gain：

Gain＝max(DN_main)/Scale (1)

其中，DN_main为直方图分布频率大于所述预定分布频率对应的灰度值数据集；max为取最大值函数，Scale为压缩比例，取值设为9。

进一步的，根据所述压缩系数对各预定波段的灰度值进行十进制编码压缩：

DN_new＝round(DN_old/Gain) (2)

其中，DN_old是原始多光谱遥感影像的灰度值，DN_new是经过压缩后的灰度值，取值范围为0-9，round是取整函数。

进一步的，所述对被压缩的灰度值建立所述预定波段对应的十进制编码包括：

假设预定波段共有N个；

十进制编码＝DN_new1*1+DN_new2*10+……DN_newi*10^i-1+……DN_newN*10^N-1；

其中，DN_new1为第1个波段的被压缩的灰度值，DN_new2为第2个波段的被压缩的灰度值，DN_newi为第i个波段的被压缩的灰度值，i∈N,DN_newN为第N个波段的被压缩的灰度值。

进一步的，所述预定波段包括四个波段，分别是蓝波段，绿波段，红波段和近红外波段。

进一步的，对所述十进制编码进行聚类分析，获得参考十进制编码的步骤包括：

对所有的十进制编码进行直方图统计，得到十进制编码对应的直方图分布频率；

提取分布频率大于一预定频率的十进制编码，作为M类初始化聚类中心，M为分布频率大于一预定频率的十进制编码的类别数；

对初始化聚类中心进行聚类分析，得到K类参考十进制编码，K为最终分类后的类别数，且K<M。

进一步的，所述聚类分析的方法包括最小距离法和光谱夹角法；

采用最小距离法进行所述聚类分析，确定参考十进制编码的步骤包括：

将十进制编码转换为N个被压缩后灰度值组成的数组，所述压缩后灰度值为0-9之间的数，N为预定波段的波段数，并记录每个十进制编码对应的分布频率；

依次计算每个十进制编码对应压缩后灰度值数组与其他十进制编码对应压缩后灰度值数组之间的相对距离：

其中，K_a，i是十进制编码a在第i波段对应的压缩后灰度值，K_b，i是十进制编码b在第i波段对应的压缩后灰度值，Dis为两个十进制编码对应的压缩后灰度值数组的相对距离，N为预定波段的波段数；

计算所有压缩后灰度值数组之间的相对距离，找出相对距离最小的两组压缩后灰度值数组，对应的十进制编码分别表示为十进制编码A和十进制编码B；

若十进制编码A对应的分布频率大于十进制编码B对应的分布频率，则将十进制编码B去掉，聚类中心数目由M类减小为M-1类；反之，则将十进制编码A去掉，聚类中心数目由M类减小为M-1类；

判断聚类后的类别数是否达到指定类别数，如果是，返回最终的参考十进制编码；如果否，回到所述计算相对距离，并对最小距离的两类进行合并的步骤；

或

采用光谱夹角法进行所述聚类分析，确定参考十进制编码的步骤包括：

将十进制编码转换为压缩后灰度值组成的数组，所述压缩后灰度值为0-9之间的数，并记录每个十进制编码对应的分布频率；

依次计算每个十进制编码对应压缩后灰度值数组与其他十进制编码对应压缩后灰度值数组之间的光谱夹角：

其中，K_a，i是十进制编码a在第i波段对应的灰度值，K_b，i是十进制编码b在第i波段对应的灰度值，Angle为两个十进制编码对应的灰度值数组的光谱夹角，N为预定波段的波段数；

计算所有压缩后灰度值数组之间的光谱夹角，找出光谱夹角最小的两组压缩后灰度值数组，对应的十进制编码分别表示为十进制编码A和十进制编码B；

若十进制编码A对应的分布频率大于十进制编码B对应的分布频率，则将十进制编码B去掉，聚类中心数目由M类减小为M-1类；反之，则将十进制编码A去掉，聚类中心数目由M类减小为M-1类；

判断聚类后的类别数是否达到指定类别数，如果是，返回最终的参考十进制编码；如果否，回到所述计算光谱夹角，并对最小夹角的两类进行合并的步骤。

进一步的，所述给所有十进制编码和参考十进制编码分配分类赋值，建立索引查找表的步骤包括：

依次计算每个十进制编码对应压缩后灰度值数组与所有参考十进制编码对应的压缩后灰度值数据之间的相似度，相似度可采用最小距离法或光谱夹角法进行计算，

找到相似度最高的参考十进制编码，将该十进制编码赋予与所述相似度最高的参考十进制编码相同的分类赋值，以建立相应的索引查找表；其中，相似度最高为相对距离最小或光谱夹角最小。

本发明的第二方面提供了一种多光谱遥感影像分类系统，包括：

多光谱遥感影像读取单元，用于读取多光谱遥感影像并获取所述影像中每个像素在预定波段的灰度值；

灰度直方图统计单元，分别统计所述多光谱遥感影像的所有像素在各预定波段的灰度直方图，获取每个灰度值对应的直方图分布频率；

十进制编码压缩处理单元，对直方图分布频率大于一预定分布频率的灰度值进行十进制编码压缩处理，建立所述灰度值在预定波段对应的十进制编码；

聚类分析单元，对所述十进制编码进行聚类分析，获得参考十进制编码；

十进制编码查找表建立单元，给所有十进制编码和参考十进制编码分配分类赋值，建立十进制编码查找表；

分类保存单元，根据所述索引查找表对所述多光谱遥感影像的灰度值进行分类并保存分类结果。

综上所述，本发明提供了一种多光谱遥感影像分类方法和系统，与现有技术相比，本发明有如下有益的技术效果：

(1)本发明的分类方法和系统的精度高，分类精度要高于传统的ISODATA和k-means方法。

(2)本发明的处理速度快。针对一景高分一号卫星影像，本发明的处理速度分别是Isodata方法和Kmean方法的5.1倍和3.7倍。

(3)本发明只需要输入指定类别数一个参数，整个分类过程不需要人工干预，可实现多光谱遥感影像分类的批处理。

附图说明

图1是本发明的多光谱遥感影像分类方法的流程示意图；

图2是本发明多光谱遥感影像分类方法中十进制编码压缩灰度值以获得十进制编码的方法流程示意图；

图3(a)是原始灰度值多光谱遥感影像，图3(b)是压缩后的多光谱遥感影像；

图4是本发明多光谱遥感影像分类方法中对十进制编码进行聚类分析以获得参考十进制编码的方法流程示意图；

图5是最小距离法的方法流程示意图；

图6是光谱夹角法的方法流程示意图；

图7是本发明的索引查找表建立示意图；

图8(a)是原始多光谱遥感影像，图8(b)是本发明分类结果多光谱遥感影像，图8(c)是ISODATA方法分类结果多光谱遥感影像，图8(d)是K-mean方法分类结果多光谱遥感影像；

图9是植被区不同非监督分类方法的分类结果图；

图10是非植被区不同非监督分类方法的分类结果图；

图11是本发明的多光谱遥感影像分类系统的结构框图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明了，下面结合具体实施方式并参照附图，对本发明进一步详细说明。应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。

本发明提出一种基于十进制编码和图像索引的多光谱遥感影像非监督分类方法，用于多光谱遥感影像数据的快速自动分类。当前，虽然利用计算机进行遥感图像的分类处理具有处理速度快等优点，在遥感领域得到了广泛的应用。但从分类精度上，对小区域分类时在大多数情况下不如人工目视解译结果。本发明将参考肉眼人工识别的原理，提出一种新的分类方法。

人眼在进行遥感图像分类时，并不具备对图像灰度值的精细识别能力。研究表明，人眼对灰度值的快速识别能力约为10种左右，也就是说，若仅根据图像灰度值特征，最多可将单波段灰度值图像分为10类。四个波段的灰度值图像，最多可分为10000类。而在实际应用中，我们在一幅多光谱影像中最多可分十几类影像。因此，只要将这10000类聚合成十几类，即可完成多光谱遥感影像的非监督分类。

因此，本方法将借鉴人眼对灰度值的识别能力，首先采用十进制编码方法对多光谱遥感影像进行压缩。将多个波段(以下以四个波段为例)的所有灰度值转换成四位数的十进制编码。然后采用最小距离法或光谱夹角法对十进制编码进行聚类分类，得到参考十进制编码。最后，根据最小距离法或光谱夹角法计算所有十进制编码对应的类别属性，通过索引的方式实现整个多光谱遥感影像的自动分类。

如图1所示，本发明的第一方面提供了一种多光谱遥感影像分类方法，详细的流程如下：

步骤S100，读取多光谱遥感影像。读取待分类的灰度值多光谱遥感影像，获取整个多光谱遥感影像在预定波段的灰度值。优选地，该预定波段的数量为4个。

步骤S200，统计灰度直方图。分别统计多光谱遥感影像在各预定波段的灰度直方图，获取每个灰度值对应的直方图分布频率。

步骤S300，十进制编码压缩灰度值以获得十进制编码。对直方图分布频率中大于一预定分布频率的灰度值进行十进制编码压缩处理，建立所述灰度值在预定波段对应的十进制编码。优选地，该预定分布频率为0.01％。具体的，如图2所示，包括如下步骤：

步骤S310，按照下式获取压缩系数Gain：

Gain＝max(DN_main)/Scale (1)

其中，DN_main为直方图分布频率大于所述预定分布频率对应的灰度值数据集；max为取最大值函数，Scale为压缩比例，因为十进制压缩，取值设为9。

步骤S320，根据该压缩系数对各预定波段的灰度值进行十进制编码压缩，得到被压缩的灰度值DN_new，按照下式计算：

DN_new＝round(DN_old/Gain) (2)

其中，DN_old是原始多光谱遥感影像的灰度值，DN_new是压缩后灰度值，取值范围为0-9，round是取整函数。

如图3所示为压缩前后多光谱遥感影像的灰度值比较，其中图3(a)是原始灰度值多光谱遥感影像，图3(b)是压缩后的多光谱遥感影像，可以看出，压缩后多光谱遥感影像保留了原始多光谱遥感影像的绝大部分信息，且经过压缩后，更容易把相似的光谱归为一类。

步骤S330，对被压缩的灰度值建立所述预定波段对应的十进制编码，具体按照如下进行计算：

十进制编码＝DN_new1*1+DN_new2*10+……DN_newi*10ⁱ+……DN_newN*10^N-1；

其中，假设预定波段共有N个，DN_new1为第1个波段的被压缩的灰度值，DN_new2为第2个波段的被压缩的灰度值，DN_newi为第i个波段的被压缩的灰度值，i∈N,DN_newN为第N个波段的被压缩的灰度值。

在一个具体的实施例中，对进行十进制编码压缩后的四个波段多光谱遥感影像灰度值建立十进制编码，即，十进制编码＝DN_new1*1+DN_new2*10+DN_new3*100+DN_new4*1000，得到一个四位数的索引图。为使处理过程更容易理解，举例说明如下：假设在一幅图像中，某一个像元的四个波段灰度值分别为[30，75，66，130]，对应的有效最大值分别为[240，230，250，235]，则压缩系数分别为240/9＝26.67，230/9＝25.56，250/9＝27.78，235/9＝26.11。则四个波段的进行十进制编码压缩处理的灰度值分别为1，3，2，5。对应的十进制编码即为5231。

步骤S400，对十进制编码进行聚类分析，获得参考十进制编码。如图4所示，具体包括如下步骤：

步骤S410，对所有的十进制编码进行直方图统计，得到十进制编码对应的直方图分布频率；

步骤S420，提取分布频率大于一预定分布频率的十进制编码，获取初始化聚类中心，聚类中心数目为M。优选地，该一预定分布频率可取值为0.01％。例如，频率超过0.01％的数目有40个，M即为40。

步骤S430，对初始化聚类中心进行聚类分析，得到聚类后的参考十进制编码。

在聚类分析中，有两种相似度计算方法，分别为最小距离法和光谱夹角法。

最小距离法的计算过程如图5所示，包括如下步骤：

步骤A1，将十进制编码转换为N个压缩后灰度值(0-9之间的数)组成的数组，其中，N为预定波段的波段数。优选地，N取值为4。

步骤A2，依次计算每个十进制编码对应的压缩后灰度值数组与其他压缩后灰度值数组之间的相对距离：

其中，K_a，i是十进制编码a在第i波段的值，K_b，i是十进制编码b在第i波段的值，i∈N，Dis为两个十进制编码的相对距离。

步骤A3，计算所有压缩后灰度值数组之间的相对距离，找出相对距离最小的两组压缩后灰度值数组，对应的十进制编码分别表示为十进制编码A和十进制编码B。若十进制编码A对应的分布频率大于十进制编码B对应的分布频率，则将十进制编码B去掉，聚类中心数目由M类减小为M-1类。反之，则将十进制编码A去掉，聚类中心数目由M类减小为M-1类，并进入下一步。

步骤A4，判断聚类后的类别数是否达到指定类别数，如果是，返回最终的参考十进制编码，即达到指定类别数后剩余的十进制编码；如果否，回到上述判断计算相对距离，并对最小距离的两类进行合并的步骤。

为了更容易理解，给出一个计算十进制编码距离的实例。假设两个十进制编码分别为2378，6214，则对应的数组分别为[8，7，3，2]和[4，1，2，6]。两者的相对距离为：

Dis＝sqrt((8-4)^2+(7-1)^2+(3-2)^2+(2-6)^2)＝sqrt(69)＝8.3066

例如，假设最初有40个聚类样本，进行一次聚类，将40个类别减少至39个。将最相似的两个合并，然后重新进行聚类分析，依次将聚类样本数量减小到38，37……，直到指定的数目后停止。

按照上述步骤，对所有十进制编码进行聚类分析，直到最后类别数达到指定类别数目停止计算，得到最终基于最小距离的参考十进制编码。

光谱夹角法的计算过程如图6所示，包括如下步骤：

步骤B1，将十进制编码转换为N个压缩后灰度值(0-9之间的数)组成的数组，其中，N为预定波段的波段数。优选地，N取值为4。

步骤B2，依次计算每个十进制编码与其他十进制编码之间的光谱夹角：

其中，K_a，i是十进制编码a在第i波段的值，K_b，i是十进制编码b在第i波段的值，Angle为两个十进制编码的光谱夹角。

步骤B3，计算所有压缩后灰度值数组之间的光谱夹角，找出光谱夹角最小的两组压缩后灰度值数组，对应的十进制编码分别表示为十进制编码A和十进制编码B。若十进制编码A对应的分布频率大于十进制编码B对应的分布频率，则将十进制编码B去掉，聚类中心数目由M类减小为M-1类。反之，则将十进制编码A去掉，聚类中心数目由M类减小为M-1类，并进入下一步。

步骤B4，判断聚类后的类别数是否达到指定类别数，如果是，返回最终的参考十进制编码，即达到指定类别数后剩余的十进制编码；如果否，回到所述计算光谱夹角，并对最小夹角的两类进行合并的步骤。

利用光谱夹角法对所有十进制编码进行聚类分析，直到最后类别数达到指定类别数目停止计算，得到最终的基于光谱夹角的参考十进制编码。

步骤S500，给所有十进制编码和参考十进制编码分配分类赋值，建立索引查找表。

如图7所示，依次计算每个十进制编码对应压缩后灰度值数组与所有参考十进制编码对应的压缩后灰度值数据之间的相似度，相似度的计算可采用最小距离法或光谱夹角法，按照上述公式(3)和(4)计算，找出每个十进制编码与所有参考十进制编码中相似度最高的参考十进制编码，将该十进制编码赋予和该相似度最高的参考十进制编码相同的分类赋值，建立相应的索引查找表。其中，相似度最高为相对距离最小或光谱夹角最小。

步骤S600，根据所述索引查找表对所述多光谱遥感影像的灰度值进行分类并保存分类结果。

根据十进制编码和对应的参考十进制编码之间的索引关系，将具有十进制编码对应的像元都给予相应的分类赋值，最终实现整个多光谱遥感影像的类别赋值，最后，保存分类结果。

选择ENVI软件中自动的一景Landsat卫星反射率影像为例，分别采用ENVI提供的ISODADA法、K-means法和本发明提出的本发明进行分类。图8示出了原始多光谱遥感影像和三种非监督分类方法的分类结果多光谱遥感影像，其中，图8(a)是原始多光谱遥感影像，图8(b)是本发明分类结果多光谱遥感影像，图8(c)是ISODATA方法分类结果多光谱遥感影像，图8(d)是K-mean方法分类结果多光谱遥感影像。从图中可以看出，ISODATA法和K-means法分类结果基本相同，本发明得到的结果有较大差异。其中最明显的地方有两处，一是在植被区传统方法分类的植被类型更多，而本发明对植被类型又进行细分；二是在非植被区，传统方法的类型更多，更破碎，而本发明的分类结果更完整。为进一步直观显示，图9和图10分别给出了植被和非植被区的放大图。

此外，为了定量分析三种方法对大数据遥感影像的处理速度，选择一景高分一号卫星PMS影像进行分类。分别计算三种方法的处理速度，如下表所示。

表1三种非监督分类方法处理速度

	Isodata方法	Kmean方法	本发明
				处理一景所需时间	189秒	138秒	37秒

可以看出，本发明的处理速度远高于Isodata方法和Kmean方法，处理时间仅为Isodata方法的19.5％和Kmean方法的26.7％,即相对Isodata方法和Kmean方法，本发明的处理速度分别提高了4.1倍和2.7倍。

本发明的第二方面提供了一种多光谱遥感影像分类系统，如图11所示，包括：

多光谱遥感影像读取单元1，用于读取多光谱遥感影像并获取所述多光谱遥感影像中每个像素在预定波段的灰度值；

灰度直方图统计单元2，分别统计所述多光谱遥感影像的所有像素在各预定波段的灰度直方图，获取每个灰度值对应的直方图分布频率；

十进制编码压缩处理单元3，对直方图分布频率大于一预定分布频率的灰度值进行十进制编码压缩处理，建立所述灰度值在预定波段对应的十进制编码；

聚类分析单元4，对所述十进制编码进行聚类分析，获得参考十进制编码；

索引查找表建立单元5，给所有十进制编码和参考十进制编码分配分类赋值，建立索引查找表；

分类保存单元6，根据所述索引查找表对所述多光谱遥感影像的所有像素对应的灰度值进行分类并保存分类结果。

综上所述，本发明提供了一种多光谱遥感影像分类方法和系统。该多光谱遥感影像分类方法采用十进制编码对多波段进行压缩，将多个波段的灰度值转换成与十进制编码，采用最小距离法或光谱夹角法对所有十进制编码进行聚类分类，得到参考十进制编码，并根据最小距离法或光谱夹角法计算所有十进制编码对应的类别属性，最后，通过索引的方式实现整个多光谱遥感影像的自动分类。本发明的多光谱遥感影像分类方法处理相比于现有的分类方法具有精度更高、处理速度更快的特点。

应当理解的是，本发明的上述具体实施方式仅仅用于示例性说明或解释本发明的原理，而不构成对本发明的限制。因此，在不偏离本发明的精神和范围的情况下所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。此外，本发明所附权利要求旨在涵盖落入所附权利要求范围和边界、或者这种范围和边界的等同形式内的全部变化和修改例。

Claims (10)

1.一种多光谱遥感影像分类方法，其特征在于，包括如下步骤：

读取多光谱遥感影像，获取所述影像中每个像素在预定波段的灰度值；

分别统计所述多光谱遥感影像的所有像素在各预定波段的灰度直方图，获取每个灰度值对应的直方图分布频率；

对直方图分布频率大于一预定分布频率的灰度值进行十进制编码压缩处理，建立所述灰度值在预定波段对应的十进制编码；

对所述十进制编码进行聚类分析，获得参考十进制编码；

为所有十进制编码和参考十进制编码分配分类赋值，建立索引查找表；

根据所述索引查找表对所述多光谱遥感影像的灰度值进行分类并保存分类结果。

2.根据权利要求1所述的多光谱遥感影像分类方法，其特征在于，所述对直方图分布频率大于一预定分布频率的灰度值进行十进制编码压缩处理，建立所述灰度值在预定波段对应的十进制编码的步骤包括：

获取压缩系数；

根据所述压缩系数对各预定波段的灰度值进行十进制编码压缩，得到被压缩的灰度值；

对被压缩的灰度值建立所述预定波段对应的十进制编码。

3.根据权利要求2所述的多光谱遥感影像分类方法，其特征在于，按照下式获取压缩系数Gain：

Gain＝max(DN_main)/Scale (1)

其中，DN_main为直方图分布频率大于所述预定分布频率对应的灰度值数据集；max为取最大值函数，Scale为压缩比例，取值设为9。

4.根据权利要求3所述的多光谱遥感影像分类方法，其特征在于，根据所述压缩系数对各预定波段的灰度值进行十进制编码压缩：

DN_new＝round(DN_old/Gain) (2)

其中，DN_old是原始多光谱遥感影像的灰度值，DN_new是经过压缩后的灰度值，取值范围为0-9，round是取整函数。

5.根据权利要求4所述的多光谱遥感影像分类方法，其特征在于，所述对被压缩的灰度值建立所述预定波段对应的十进制编码包括：

假设预定波段共有N个；

十进制编码＝DN_new1*1+DN_new2*10+……DN_newi*10^i-1+……DN_newN*10^N-1；

其中，DN_new1为第1个波段的被压缩的灰度值，DN_new2为第2个波段的被压缩的灰度值，DN_newi为第i个波段的被压缩的灰度值，i∈N,DN_newN为第N个波段的被压缩的灰度值。

6.根据权利要求1-5任一项所述的多光谱遥感影像分类方法，其特征在于，所述预定波段包括四个波段，分别是蓝波段，绿波段，红波段和近红外波段。

7.根据权利要求1-5任一项所述的多光谱遥感影像分类方法，其特征在于，对所述十进制编码进行聚类分析，获得参考十进制编码的步骤包括：

对所有的十进制编码进行直方图统计，得到十进制编码对应的直方图分布频率；

提取分布频率大于一预定频率的十进制编码，作为M类初始化聚类中心，M为分布频率大于一预定频率的十进制编码的类别数；

对初始化聚类中心进行聚类分析，得到K类参考十进制编码，K为最终分类后的类别数，且K<M。

8.根据权利要求7所述的多光谱遥感影像分类方法，其特征在于，所述聚类分析的方法包括最小距离法和光谱夹角法；

采用最小距离法进行所述聚类分析，确定参考十进制编码的步骤包括：

将十进制编码转换为N个被压缩后灰度值组成的数组，所述压缩后灰度值为0-9之间的数，N为预定波段的波段数，并记录每个十进制编码对应的分布频率；

依次计算每个十进制编码对应压缩后灰度值数组与其他十进制编码对应压缩后灰度值数组之间的相对距离：

其中，K_a,i是十进制编码a在第i波段对应的压缩后灰度值，K_b,i是十进制编码b在第i波段对应的压缩后灰度值，Dis为两个十进制编码对应的压缩后灰度值数组的相对距离，N为预定波段的波段数；

计算所有压缩后灰度值数组之间的相对距离，找出相对距离最小的两组压缩后灰度值数组，对应的十进制编码分别表示为十进制编码A和十进制编码B；

若十进制编码A对应的分布频率大于十进制编码B对应的分布频率，则将十进制编码B去掉，聚类中心数目由M类减小为M-1类；反之，则将十进制编码A去掉，聚类中心数目由M类减小为M-1类；

判断聚类后的类别数是否达到指定类别数，如果是，返回最终的参考十进制编码；如果否，回到所述计算相对距离，并对最小距离的两类进行合并的步骤；

或

采用光谱夹角法进行所述聚类分析，确定参考十进制编码的步骤包括：

将十进制编码转换为压缩后灰度值组成的数组，所述压缩后灰度值为0-9之间的数，并记录每个十进制编码对应的分布频率；

依次计算每个十进制编码对应压缩后灰度值数组与其他十进制编码对应压缩后灰度值数组之间的光谱夹角：

其中，K_a,i是十进制编码a在第i波段对应的灰度值，K_b,i是十进制编码b在第i波段对应的灰度值，Angle为两个十进制编码对应的灰度值数组的光谱夹角，N为预定波段的波段数；

计算所有压缩后灰度值数组之间的光谱夹角，找出光谱夹角最小的两组压缩后灰度值数组，对应的十进制编码分别表示为十进制编码A和十进制编码B；

若十进制编码A对应的分布频率大于十进制编码B对应的分布频率，则将十进制编码B去掉，聚类中心数目由M类减小为M-1类；反之，则将十进制编码A去掉，聚类中心数目由M类减小为M-1类；

判断聚类后的类别数是否达到指定类别数，如果是，返回最终的参考十进制编码；如果否，回到所述计算光谱夹角，并对最小夹角的两类进行合并的步骤。

9.根据权利要求8所述的多光谱遥感影像分类方法，其特征在于，所述为所有十进制编码和参考十进制编码分配分类赋值，建立索引查找表的步骤包括：

依次计算每个十进制编码对应压缩后灰度值数组与所有参考十进制编码对应的压缩后灰度值数据之间的相似度，相似度可采用最小距离法或光谱夹角法进行计算，

找到相似度最高的参考十进制编码，将该十进制编码赋予与所述相似度最高的参考十进制编码相同的分类赋值，以建立相应的索引查找表；其中，相似度最高为相对距离最小或光谱夹角最小。

10.一种多光谱遥感影像分类系统，其特征在于，包括：

多光谱遥感影像读取单元，用于读取多光谱遥感影像并获取所述影像中每个像素在预定波段的灰度值；

灰度直方图统计单元，分别统计所述多光谱遥感影像的所有像素在各预定波段的灰度直方图，获取每个灰度值对应的直方图分布频率；

十进制编码压缩处理单元，对直方图分布频率大于一预定分布频率的灰度值进行十进制编码压缩处理，建立所述灰度值在预定波段对应的十进制编码；

聚类分析单元，对所述十进制编码进行聚类分析，获得参考十进制编码；

索引查找表建立单元，给所有十进制编码和参考十进制编码分配分类赋值，建立索引查找表；

分类保存单元，根据所述索引查找表对所述多光谱遥感影像的灰度值进行分类并保存分类结果。

Images