基于时间序列重要点分析的高光谱遥感图像波段选择方法

摘要

本发明公开了一种基于时间序列重要点分析的高光谱遥感图像波段选择方法，本发明方法通过基于DBI的K-means聚类方法将高光谱遥感图像数据样本聚类成K个类别，并使用小波分析进行去噪处理，再基于时间序列分析进行重要点的提取，以实现特征波段的选择。相比现有技术，本发明方法具有计算复杂度低、实现过程方便快捷的优点，并为高维数据的降维提供了一种全新的思路。

说明书

技术领域

本发明涉及图像处理领域，尤其是涉及一种高光谱遥感图像波段选择方法。

背景技术

随着科学技术的快速发展和航空航天技术及遥感科学的飞速进步，可获取的高光谱遥感数据正以惊人的速度增长，如何处理与利用这样丰富的光谱信息，又希望提高学习效率和速度，是人们研究的热门课题。处理高光谱遥感数据一个最重要的前提是要对波段数目进行降维，常用的降维方式分为特征提取与特征选择，此时特征选择即波段选择。通常特征提取会改变光谱波段的原始物理意义，不利于地物反演，因此许多研究者为了充分利用高光谱遥感的丰富波段信息，积极研究波段选择方法。

目前的波段选择方法较多，如常用的搜索算法，有穷举搜索法、顺序前进或后退法、分支定界法等。这类方法用于维数不多的特征选择时，计算量与计算复杂度都还可以接受，但是用于成百上千的高光谱图像波段选择时，由于计算效率较低而不实用。刘春红等人[刘春红，赵春晖，张凌雁：一种新的高光谱遥感图像降维方法.中国图象图形学报，2005，10(2)：218-222]便提出了自适应波段选择方法(ABS，Adaptive Band Selection)。该方法充分考虑了高光谱图像的空间相关性和谱间相关性，并通过计算各个波段的指数来选择信息量大并且与其他波段相关性小的波段，效果较好，其计算复杂度也大大降低，但是这种方法选择出的波段分布不均匀，有信息损失，不能充分利用高光谱图像丰富的光谱信息。还有人[Chavez P S，Berlin G L，Sowers L B：Statistical method for selectingLandsat MSS ratios.Journal of applied photographic engineering，1982，1(8)：23-30.]提出基于信息量的最佳波段选择方法，通过计算高光谱图像的熵与联合熵、组合波段的协方差矩阵行列式或者最佳指数(OIF，Optimum Index Factor)，并将这些反应信息量的指标值按照从大到小的顺序排序，根据需要选择最优组合波段。这些计算指标通常都能找到信息量最大的波段组合，但是都没有考虑到不同地物类别之间在组合波段上区分的难易度。而刘建平等人[刘建平，赵英时，孙淑玲：高光谱遥感数据最佳波段选择方法试验研究.遥感信息，2001，1(1)：7-13.]针对这个问题提出了基于类间可分性的最佳波段组合选择方法。实验证明此类方法能够找到最佳波段组合，实现维数降低。但是在实际应用中，存在一定的局限性，并不能总是找到最优组合，而且不适用于高维数据。王立国等人[王立国，谷延锋，张晔：基于支持向量机和子空间划分的波段选择方法.系统工程与电子技术，2005，27(6)：974-977.]提出了基于支持向量机和子空间划分的波段选择方法，该方法通过对支持向量机判决函数进行敏感度分析，同时对数据源进行子空间划分，然后结合敏感度分析结果和子空间划分结果来实现有效的波段选择。此方法有效去除了冗余，并大幅度的提高了分类精度，可行有效，但计算过程复杂。

综上所述，现有的高光谱遥感图像波段选择方法均存在计算复杂度高的问题，对图像处理系统的软、硬件配置也提出了极高的要求。

发明内容

本发明的目的在于提供一种计算复杂度低、实现过程方便快捷的高光谱遥感图像波段选择方法。

本发明是通过将时间序列分析方法引入高光谱遥感图像的波段选择来实现上述目的。下面先对时间序列分析方法的原理和内容进行简单介绍：

时间序列分析方法是一种动态数据处理的统计方法。从统计的意义上讲，所谓时间序列，就是根据某一个指标在不同的时间上顺序记录的一系列的数据。严格按照上述定义来说，时间序列数据都是按照时间顺序变化的，但是王振龙[王振龙：时间序列分析.北京：中国统计出版社，2000.]提出只要观测量按照一定的顺序排列而成均可认为是时间序列。那么引申而出，一定的顺序可以是时间顺序，也可以是其他物理意义的物理量，例如长度，温度，速度等。因此，时间序列只强调顺序的重要性，而没有强调必须要用时间顺序排列，因此完全可以将时间序列与高光谱遥感图像数据联系起来。以Washington DC Mall数据集的为例，下表是其中第一类中一个样本点的数据：

波段号           1            2    …         190         191

波段值    0.757223     0.824438    …    0.479944    0.490868

对于上表中的数据，如果看成是样本的波段值按照波段号这个递增的顺序排列而成，那么它就是一个时间序列，这样，不仅符合了时间序列的定义，而且为分析高光谱图像数据提供了更加简单直观的工具，依次类推，高光谱图像的数据集实际上就是由多个时间序列组成的。如Washington DC Mall数据集分成训练集与测试集，训练集样本总共有4428个，也就是训练集中有4428个时间序列。而每个时间序列都是191个“时间”点，有191个数值，属于离散时间序列。

一般来说，时间序列往往具有某种趋势性或者呈现出一定的周期性。那么如何根据这些时间序列找出其隐含的趋势性、周期性与内在统计规律，则是时间序列分析要解决的问题。时间序列分析在很多领域都有极其重要的作用，而在高光谱图像处理中，则很少有人引用时间序列分析。

肖辉和胡运发[肖辉，胡运发：基于分段时间弯曲距离的时间序列挖掘，计算机研究与发展，2005，42(1)：72-77]提出在数据挖掘领域，时间序列的极值点通常带有比较多的信息，具有重要分析意义。Perng c s等人[Perng c s，Wang h，Zhang s r，et al：Landmarks：a new model for similarity-based pattern querying in time series databases.Proceedings of the 16th International Conference on Data Engineering，San Diego：IEEE，2000：33-42.]则根据时间序列的关键点(与极值点类似)进行相似性匹配，大大加速了搜索速度，证明了关键点的重要处理价值。周黔、吴铁军[周黔，吴铁军：基于重要点的时间序列趋势特征提取方法.浙江大学学报，2007，41(11)：1782-1787.]提出时间序列重要点的概念，并将其应用于时间序列趋势特征提取中，效果较好。关键点与重要点都与时间序列的极值点非常类似，它们都带有时间序列本身重要的信息，对分析时间序列有重要作用，因此这些点可认为是重要信息点。时间序列的重要点描述了时序的趋势变化特征，本身带有重要信息，因此完全可以在波段选择中引入重要点分析，达到降低数据维数，且保持分类精度的目的。通过实际试验观察发现，高光谱遥感图像样本数据的时间序列曲线趋势特征复杂，上下波动频繁，噪声较多，不利于后期的重要点的提取，因此需要进行去噪预处理，本发明选择小波变换的方法来进行去噪预处理。

综上所述，可以得到本发明的技术方案如下：

一种基于时间序列重要点分析的高光谱遥感图像波段选择方法，其特征在于，包括以下各步骤：

步骤1)通过K-means聚类方法将原始高光谱遥感图像数据聚类成K个类别，并保留每个类别的聚类中心；

步骤2)将步骤1得到的所有聚类中心看成一系列的时间序列，用小波变换的方法对每一时间序列进行平滑去噪处理；

步骤3)分别提取步骤2得到的平滑去噪处理后的时间序列中的重要点，然后将得到的所有重要点进行合并，合并后保留的重要点所对应的波段即为选择的波段；其中所述重要点的定义为：

对于时间序列n是时间序列X的长度，t_i表示第i个时间坐标，X(t_i)表示时间序列在第i个时间坐标的取值；定义X的第q个重要点其中P_q∈{1，2，…，n}表示第q个重要点在时间序列中的位置，为满足以下条件的数据点：

或者

或者

或者

本发明方案中采用的小波变换方法是现有技术，具体内容可参见文献[Daubechies，I.1992.Ten Lectures on Wavelets.Society for Industrial and Applied Mathematics，Philadelphia，PA，USA]。在使用小波变换方法进行去噪处理时，小波基函数的选择和分解层次的确定直接影响到去噪的效果，经大量实验发现，选择sym4小波基函数进行分解，分解层数为4效果最佳。

上述技术方案中所使用的K-means聚类方法是图像处理领域常用的方法，主要目的是将n个数据对象最后分成k个聚类，并且使得每个类别中的数据点到该聚类中心的平方和最小。它属于无监督的聚类方法，无需事先知道样本的输入顺序，因此它不仅可以处理大的数据集，而且也是相对可伸缩的和高效率的。K-means聚类方法的具体内容可参见文献[J.B.McQueen，Some methods of classification and analysis of multivariateobservations，Proc.5^th Berkeley Symp.on Mathematical Statistics and Probability，1967，vol.1：281-297.]，此处不再赘述。在使用K-means聚类方法时，需要用户事先确定聚类数目K的取值，往往K的值需要根据个人经验多次尝试，才能找到最佳的聚类个数K，如果用户给出的K的取值不恰当，会影响到最终聚类的效果，进而导致波段选择的准确性下降。

针对K-means聚类算法本身的一些缺陷，为了实现聚类过程完全自动化，本发明提出在K-means算法的基础之上使用DBI(Davies-Bouldin Index)有效性指标来进行最佳聚类数目K的确定。DBI指标是一种聚类有效性的指标，主要是关于同一类中所有样本的紧密程度与不同类样本之间分散程度的一个函数。DBI指标主要就是利用几何原理，通过分别计算同一类别样本之间的相似性及不同类别之间的相异性，再取它们的比值，然后根据此比值(即DBI值)来评价聚类的效果。当同类样本之间的相似性大，不同类别之间的相异性大的时候，DBI值较小，表示聚类效果较好。

DBI的计算公式：

$\mathrm{DBI}=\frac{1}{K}\underset{i=1}{\overset{K}{\Sigma }}\underset{1<=j<=K,j\ne i}{\max }\left\{\frac{S_i+S_j}{d_{i,j}}\right\}$

其中，

$S_i=\frac{1}{\left|C_i\right|}\underset{x\in C_i}{\Sigma }\left|\right|x-v_i||,$

d_i，j＝‖v_i-v_j‖

S_i表示第i类样本之间的相似度；S_j表示第j类样本之间的相似度；d_i，j表示第i类样本与第j类聚类中心之间的距离；K表示聚类数目；C_i表示属于第i类的样本个数；v_i表示第i个聚类的类别中心；v_j表示第j个聚类的类别中心；x表示属于第i类的各个样本；i和j均为大于0且小于等于K的整数；公式中所用的距离计算方法都是指欧几里德距离。

因此可以在K-means聚类过程中，通过不断计算DBI的值，找到最小DBI对应的聚类，这个聚类结果就是最佳聚类结果，相应的K值就是最佳聚类数目。可以根据具体数据集设定一个K的初始值，一直不断的增加到某一预先设定的阈值，本发明的技术方案中，K的初始值取2，考虑到计算复杂度和计算效率，建议预先设定的阈值范围为10-50之间的整数；具体的聚类按照如下各步骤：

步骤101)设定聚类数目K的初始值为2；

步骤102)判断K是否小于或等于一个预先设定的阈值，如是，则在样本数据中等间隔的取K个聚类中心，继续执行步骤103；如否，则转而执行步骤107；

步骤103)分别计算各个样本点到各个聚类中心的距离，并且将它归类到距离最近的类别中；

步骤104)将所有样本点都归类完之后，再重新计算K个类别的聚类中心；

步骤105)比较新计算的聚类中心与前一次计算的聚类中心，如果聚类中心有改变，则转向步骤103；否则，转向步骤106；

步骤106)按照以下公式计算此时的DBI的值，然后设置K＝K+1，转至步骤102：

$\mathrm{DBI}=\frac{1}{K}\underset{i=1}{\overset{K}{\Sigma }}\underset{1<=j<=K,j\ne i}{\max }\left\{\frac{S_i+S_j}{d_{i,j}}\right\}$

其中，

$S_i=\frac{1}{\left|C_i\right|}\underset{x\in C_i}{\Sigma }\left|\right|x-v_i\left|\right|,$

d_i，j＝‖v_i-v_j‖

S_i表示第i类样本之间的相似度；S_j表示第j类样本之间的相似度；d_i，j表示第i类样本与第j类聚类中心之间的距离；K表示聚类数目；C_i表示属于第i类的样本个数；v_i表示第i个聚类的类别中心；v_j表示第j个聚类的类别中心；x表示属于第i类的各个样本；i和j均为大于0且小于等于K的整数；公式中所用的距离计算方法都是指欧几里德距离；

步骤107)比较所有的DBI的值，找出最小DBI值所对应的K，并输出该K值时的聚类结果。

本发明通过基于DBI的K-means聚类方法将高光谱遥感图像数据样本聚类成K个类别，并使用小波分析进行去噪处理，再基于时间序列分析进行重要点的提取，以实现特征波段的选择。相比现有技术，本发明方法具有计算复杂度低、实现过程方便快捷的优点，并为高维数据的降维提供了一种全新的思路。

附图说明

图1为本发明方法的流程图。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

此处仍使用HYDICE光谱仪所获取的Washington DC Mall地区的公共测试图像数据为例来具体说明本发明的具体实施方式。该图像数据的波长范围为0.40μm至2.40μm，共包含210个连续的波段，去除了无效波段后，剩余191个有效波段。该图像中包含草地、屋顶以及道路等七种类别，已经手工标出137个已知地面真实类别的区域，为保证训练和测试数据集的不重叠，取其中序号为奇数的标注区域作为训练样本，序号为偶数的标注区域作为测试样本。

使用本发明的方法对该图像进行波段选择，如附图1所示，具体按照以下各步骤：

步骤1)通过基于DBI的K-means聚类方法将原始高光谱遥感图像数据聚类成K个类别，并保留每个类别的聚类中心；本步骤具体包括以下各步骤：

步骤101)设定聚类数目K的初始值为2；

步骤102)判断K是否小于或等于一个预先设定的阈值，如是，则在样本数据中等间隔的取K个聚类中心，继续执行步骤103；如否，则转而执行步骤107；

本具体实施方式中，所述预先设定的阈值取值为30；

步骤103)分别计算各个样本点到各个聚类中心的距离，并且将它归类到距离最近的类别中；

步骤104)将所有样本点都归类完之后，再重新计算K个类别的聚类中心；

步骤105)比较新计算的聚类中心与前一次计算的聚类中心，如果聚类中心有改变，则转向步骤103；否则，转向步骤106；

步骤106)按照以下公式计算此时的DBI的值，然后设置K＝K+1，转至步骤102：

$\mathrm{DBI}=\frac{1}{K}\underset{i=1}{\overset{K}{\Sigma }}\underset{1<=j<=K,j\ne i}{\max }\left\{\frac{S_i+S_j}{d_{i,j}}\right\}$

其中，

$S_i=\frac{1}{\left|C_i\right|}\underset{x\in C_i}{\Sigma }\left|\right|x-v_i\left|\right|,$

d_i，j＝‖v_i-v_j‖

S_i表示第i类样本之间的相似度；S_j表示第j类样本之间的相似度；d_i，j表示第i类样本与第j类聚类中心之间的距离；K表示聚类数目；C_i表示属于第i类的样本个数；v_i表示第i个聚类的类别中心；v_j表示第j个聚类的类别中心；x表示属于第i类的各个样本；i和j均为大于0且小于等于K的整数；公式中所用的距离计算方法都是指欧几里德距离；

步骤107)比较所有的DBI的值，找出最小DBI值所对应的K，并输出该K值时的聚类结果。

在本具体实施方式中，K的取值与相应的DBI之间的关系如下表：

  K  7  8  9  10  11  12  13  14  DBI  0.874982  0.85155  0.907315  0.770845  0.751972  0.808533  1.057282  0.763676  K  15  16  17  18  19  20  21  22  DBI  1.023653  0.812052  0.974334  0.914839  0.992921  1.027819  1.136527  1.005891  K  23  24  25  26  27  28  29  30  DBI  1.166676  1.208495  1.112705  1.080183  1.155581  1.161763  1.215406  1.166322

从上表中可以发现，在K＝11时，获得最小的DBI的值，是0.751972，也就是说，Washington DC Mall数据集的训练样本在聚成11类的时候，各类之间相异性大，同类之间相似性大，这样聚类最合理，因此最佳聚类个数为11。

步骤2)将步骤1得到的所有聚类中心看成一系列的时间序列，用小波变换的方法对每一时间序列进行平滑去噪处理；

本具体实施方式中，选择sym4小波基函数进行分解，分解层数为4；

步骤3)分别提取步骤2得到的平滑去噪处理后的时间序列中的重要点，然后将得到的所有重要点进行合并，合并后保留的重要点所对应的波段即为选择的波段；其中所述重要点的定义为：

对于时间序列n是时间序列X的长度，t_i表示第i个时间坐标，X(t_i)表示时间序列在第i个时间坐标的取值；定义X的第q个重要点其中P_q∈{1，2，…，n}表示第q个重要点在时间序列中的位置，为满足以下条件的数据点：

或者

或者

或者

由于各时间序列可能存在相同的趋势，即存在同样的重要点，因此可将重复的波段号去掉；同时，由于高光谱图像数据本身具有高度相关性，各类样本提取出来的重要点还可能存在非常接近的现象，因此，需要对所有重要点进行合并处理。在本具体实施方式中，采用如下方法进行合并：

首先对所有重要点按照从小到大进行排序；然后从最小的重要点开始进行判断是否有其他的重要点与其距离小于5，如有，则保留该重要点，并将与其距离小于5的其它重要点去除；如果没有，则保留该重要点，并从下一个距离最近的重要点开始判断；直到所有重要点均判断完毕，此时剩余的重要点即为合并后的结果。

上述的数值5是根据本具体实施方式中所处理的图像数据选择的，当然，针对不同的图像，可以选取其它合适的数值。

经过合并后，得到了22个重要点，其所对应的22个波段即为需要得到的信息量最大并且与其他波段相关性小的波段。

为了验证本发明的效果，采用支持向量机方法分别对上述图像的原始数据以及本发明所得到的数据进行分类，并比较两者的分类准确率。实验中，采用台湾大学林智仁(LinChih-Jen)副教授等开发设计的LIBSVM2.86版本的软件包进行分类。分类结果如下表所示：

  波段数  分类准确率  原始数据  191  99.5344％  本发明方法的数据  22  99.5892％

通过实验结果发现，采用本发明方法得到的波段数据进行支持向量机分类，准确率高达99.5892％，比原始数据的分类精度还要高。因此，这22个重要信息波段保留了非常重要的有利于分类的原始信息，可以代表原始的191维数据。采用本发明的基于时间序列分析的遥感图像波段选择方法所选择出的波段数据作为原始数据的重要波段来进行后期的应用(比如分类)，可以达到减少计算量与降低计算复杂度的目的。

本发明方法可通过编程与计算机系统结合，从而实现自动化的图像处理与分析。