根据云像元数量分治的遥感图像时间序列聚类方法

摘要

本发明为了解决遥感图像时间序列聚类过程中广泛存在的云像元的干扰，提出了一种根据时间序列中云像元数量的多少分而治之的聚类方法。大致步骤为：(1)根据含云量对时间序列分组(2)对含云较少的序列，排除云像元并用动态时间规整距离进行K-means聚类(3)利用上一步得到的聚类中心作为训练样本，利用欧式距离，且不计算云像元，对含云量中等的序列进行最近邻分类，(4)对含云量很多的序列，参考其空间邻域的情况确定类别。本发明充分考虑了云对聚类结果的影响，并利用适当的方法进行了合理的规避，最大化地利用了数据，对各类传感器数据都有良好的适应能力，并能取得良好的聚类精度。

说明书

技术领域

本发明涉及遥感图像时间序列的聚类技术，具体的说，是一种根据每个遥感图像时间序列中云像元数 量的多少而采取相应策略的图像聚类方法，属于遥感图像处理技术领域。

背景技术

对地球表面变化的监测是遥感学科的主要任务之一，利用卫星遥感影像对地表覆盖进行分类或聚类是 监测地表情况的重要手段。以往的地表覆盖聚类都是基于单幅影像的，但是由于单幅影像提供的信息量不 足，且影像质量千差万别，基于单幅影像的方法始终不能取得很好的效果。卫星遥感影像到目前已经有了 几十年的数据积累，而且新发射的卫星时间分辨率越来越高，遥感图像时间序列的处理技术已经具备成熟 的发展条件。遥感图像时间序列是若干幅图像按照时间排列而成的(如图2所示)它能提供对地表覆盖类 型随时间变化情况的连续观测，因此提供关于地表覆盖类别的更完备的信息。

聚类是非监督的分类，不需要人工指定分类样本和对应的类别。基于遥感图像时间序列的聚类工作， 背景技术主要包括时间序列的相似性度量(距离)和聚类算法两大部分。目前在相关的科研领域，专门针 对遥感图像时间序列聚类的相关技术的研究不是很充分，大多数技术都是迁移自通用的时间序列处理领 域。

法国的Petitjean等人在《IEEE Transactions on Geoscience and Remote Sensing》(IEEE地 球科学与遥感)Volume.50，page.3081-3095上发表的“Satellite Image Time Series Analysis Under Time  Warping”(基于时间规整的卫星图像时间序列分析)中提出了排除遥感图像时间序列中的云像元， 然后利用动态时间规整距离可以处理长度不同的序列的特点进行聚类的方法。这种方法首次将动态时 间规整距离引入遥感图像时间序列处理领域，相当有启发性。但后来经过实验发现，这种方法在面对 部分序列云像元很多情况下效果不好。在面对不同传感器，不同空间尺度时，比较恶劣的云分布的情 况下，我们仍然缺少一套通用的遥感图像处理聚类方法。

发明内容

为了解决在大量云干扰下，有效进行遥感图像时间序列的聚类的问题，本发明提出了一套根据每个时 间序列中云像元数量的多少而采取相应策略的聚类方法。本发明的特点是先对所有时间序列按照含有云像 元的数量分组，不同组的时间序列，采取不同的距离度量和聚类方式，合理的规避云像元对聚类精度的影 响，既最大化的利用了数据，又保证了聚类的效果。

本发明是通过以下技术方案实现的，主要步骤包括：

第一步，根据时间序列中含云像元的数量将时间序列分组，共分为三组，第一组时间序列含云像元数 量少含云量小于或等于20％，第二组时间序列含云像元数量中等，含云量大于20％且小于80％，第三组时间 序列含云像元数量多，含云量大于80％，分组的分界点不限于20％和80％，但至少要保证含云量少的第一组 时间序列中的序列数量大于总数的60％；

第二步，对第一组时间序列排除序列中的云像元，然后采用基于动态时间规整距离进行K-means聚类；

第三步，对第二组时间序列采用第一组的聚类中心作为训练样本，对于序列中非云的数据点用欧氏距 离进行最近邻分类；

第四步，对第三组时间序列采用空间领域聚类。

进一步，第三组时间序列的聚类方式具体为：每个序列空间邻域中哪种类别的序列数量最多，此序列 就属于哪种类别。

进一步，对第一组时间序列进行聚类之前，先排除了序列中的云像元。

进一步，对第一组时间序列进行聚类时，使用的时间序列相似性度量为动态时间规整距离。

进一步，对第二组时间序列进行聚类时，使用的时间序列相似性度量为非云像元的数据点间的欧氏距 离。

本发明与现有的技术相比，不需要放弃那些云比较多的图像，充分的利用了所有可能的数据资源。本 发明对各种云覆盖条件下的，各种类型的传感器的数据，都有很强的适应能力。因此能满足全国乃至全球 等大空间尺度的地表聚类任务。

附图说明

图1根据含云像元数量分治的遥感图像时间序列聚类方法整体流程图

图2遥感图像时间序列概念图

图3时间序列动态时间规整距离示意图

图4只计算非云像元的时间序列欧式距离示意图

图5邻域扩展示意图

具体实施方式

下面结合附图，描述本发明的具体实施方式，图1描述了整体的流程。

(1)本方法的输入：

本方法需要的输入有两个，一是多幅不同时相的遥感图像按时间顺序排列而组成的遥感图像时间序 列；二是每幅图像相应的云掩膜数据，标记了每幅图像上的每一个像素点是否是云像元。这两个输入缺一 不可。

(2)步骤1：根据含云像元数量，对所有时间序列进行分组：

本方法的特点是根据每条遥感图像时间序列中所含有的云像元的多少，而采取不同的聚类策略，分而 治之。第一步是对所有序列按照含云像元的数量进行分组，总共分为三组，第一组序列含云像元数量小于 或等于时间序列长度的20％；第二组序列含云像元数量介于时间序列长度的20％至80％之间；第三组序列的 含云像元数量大于或等于时间序列长度的80％。除了20％和80％这两个分界点，也可以根据数据的具体情况 选择其他合适的分界点，比如10％和90％，但是一定要保证含云量少的第一组序列中的序列数量足够多(大 于总数的60％)。一般20％和80％是比较理想的分界点，不必变化。根据图2所示，每条时间序列是跟图像 上的一个像素位置一一对应的，时间序列的长度即是图像的幅数。根据输入的每幅图像的云掩膜数据，计 算出每条时间序列中云像元的数量，然后再根据上述的分组标准执行分组。

(3)步骤2：对第一组遥感图像时间序列基于动态时间规整距离进行K-means聚类：

动态时间规整距离是一种相似性度量手段，其最大的特点是可以度量长度不相同的时间序列之间的相 似性。图3给出了动态时间规整距离的示意图。但是如果两条时间序列的长度相差太大，或者数据过于破 碎，相似性度量的效果也不理想。云像元是遥感图像时间序列中不可避免的问题，云作为噪声，遮挡了真 实的数据，而且云像元的时间与空间分布是随机的，这都给时间序列的相似性度量带来了严重的问题。

利用动态时间规整距离可以处理不同长度时间序列的特点，我们排除序列中所有的云像元，得到长短 不一的序列，然后进行相似性度量。但是因为序列中的云像元数量差别很大，删除云像元之后剩余的时间 序列要么过于破碎，要么长短差别过于悬殊，都会影响相似性度量的结果。因此，我们先选出含云像元较 少的时间序列，排除它们之中的云像元，所得到的时间序列的长度差距不会很大，对它们利用动态时间规 整距离进行K-means聚类。K-means是最常用的数据聚类方式，鲁棒性很强，精度也比较高，速度较快。

含云像元量少于或等于20％的遥感图像时间序列的数量一般较多，一般占65％以上，且空间分布较均 匀，因此对于这一部分的序列进行聚类，已经可以反映整体的聚类结构，得到的聚类中心具有良好的代表 性。经过这一步，较大部分的数据确定了所属的类别标记，同时我们也得到了全体序列的大致聚类结构和 各类别的代表性序列。

动态时间规整距离的计算是一个动态规划问题，先计算两条序列所有数据点对之间的距离，再利用动 态规划算法计算完整的两条序列之间的距离。K-means是最简单最普遍的聚类方法，对算法本身不再赘述， 但是在K-means的计算过程中，需要计算序列的平均值，在动态时间规整距离下的序列平均值，采用的是 DBA(Dynamic Time Warping Barycenter Averaging)算法。DBA算法需要给定一个初始序列，然后计算 所有待平均序列与初始序列的动态时间规整距离。对于初始序列中的每一个点，根据动态时间规整距离中 的数据点对应关系，找到与之对应的所有待平均序列中的点，然后求这些点的平均值，以更新初始序列中 的点。当初始序列中的点全部更新之后，DBA算法完成一次循环，按照同样的方法继续优化初始序列，直 至初始序列不再有变化，则此时的初始序列就是待平均序列在动态时间规整距离定义下的平均序列。初始 序列的选择是从待平均序列中选择任一最长的序列。

(4)步骤3：用第一组数据的聚类中心作为训练样本，对第二组数据进行最近邻分类：

对于含云像元量中等的序列，因为无效数据较多，如果让它们参与整体聚类结构的确定，会影响聚类 精度。因此，在经过上一步得到了各类别的代表性序列之后，我们计算第二组的数据与各类别代表性序列 间的相似性，如果某序列与某一类别的代表性序列最相似，则该序列属于此类别。

在这一步中，因为序列中的云像元较多，删除云像元会导致序列支离破碎，动态时间规整距离的效果 不佳，因此，我们选择只对序列中非云像元的数据进行欧式距离度量。欧式距离要求两条时间序列的长度 相同。遥感图像时间序列在删除云像元之前的长度都是相同的，在第一组数据的处理中，我们暂时屏蔽了 其中的云像元，导致参与计算的时间序列长度不同。但是在计算每一个类别的平均序列的过程中，即得到 代表性序列的过程中，我们总是选择最长的序列最为初始序列，这样得到的平均序列的长度也是最长的， 可以说在绝大多数情况下，这个最长的序列的长度就是原始的未删除云的序列的长度，因为每一类中几乎 总存在不含云的序列。如果遇到特殊情况，则我们先对初始序列在两端进行复制延拓，得到长度与原始序 列长度相同的序列，再进行平均值的计算。在两端进行复制延拓的过程，具体的说就是如果序列的长度与 原始长度差距为1，就复制序列的第一个数据点，放在序列最前面，如果长度差距为2，就再复制序列的 最后一个数据点，放在序列最后面，循环以上过程直至长度补齐至原始长度。注意这里的平均值的概念是 定义在动态时间规整距离之下的，参与计算的序列的长度和平均值序列的长度可以是不同的。这样一来， 每一类代表性序列的长度和未删除云像元的序列的长度就是相同的，我们可以应用欧式距离度量它们之间 的相似性。计算相似性的过程中，第二组数据中被标记成云的数据点是不考虑的，只计算非云的数据点之 间的距离，如图4所示。

经过这一步，绝大部分的数据都确定了所属的类别。

(5)步骤4：基于前两组数据的分类结果，根据空间邻域对含云量较多的第三组数据进行分类。

第三组数据的含云量很多，即是说无效的数据占了绝大多数，我们认为这种数据如果依然按照时间序 列之间的相似性作为分类的依据，则会产生很大的误差，降低分类的精度。因此，我们考虑采用基于空间 邻域的分类依据。具体的说，就是对于一个遥感图像时间序列，如果它周围的大多数序列都属于同一个类 别，那么这个时间序列就被认为属于这个类别。

第三组时间序列中序列的个数一般很少，前两组序列一般涵盖了所有时间序列数量的90％以上。经过 了对前两组序列的分类之后，绝大多数的序列都已经被标记了所属的类别。因此第三组序列中的某序列的 空间邻域之中大多数都是已经有类别标记的数据。

在具体的实现当中，一个像素位置的直接邻域有八个像素点，可能会出现其中4个点属于同一类，而 另外4个点属于另一类的情况，因而无法判断所属类别。在这种情况下，则扩大邻域范围直至邻域中某一 种类别占多数，如图5所示。

(6)经过以上4个步骤，遥感图像时间序列中的所有序列都获得了所属类别的标记，完成了聚类过程。 整个过程充分考虑到了云对聚类结果的影响，并利用适当的方法进行了合理的规避，在最大化利用数据的 同时取得了良好的聚类效果。

本发明实施例已经在PC平台上实现，经过不同传感器的遥感图像时间序列数据(TM，MODIS传感器的 归一化植被指数数据)验证，取得了比传统方法更好的聚类效果。

应当指出，以上所述具体实施方式可以使本领域的技术人员更全面地理解本发明，但不以任何方式限 制本发明。因此，本领域技术人员应当理解，仍然可以对本发明进行修改或者等同替换；而一切不脱离本 发明的精神和技术实质的技术方案及其改进，其均应涵盖在本发明专利的保护范围当中。