CN105005784A - 一种基于cd-dtw距离的时间序列遥感影像土地覆盖分类方法 - Google Patents

Abstract

本发明涉及一种基于CD-DTW距离的时间序列遥感影像土地覆盖分类方法，主要步骤包括：输入多时相Landsat8遥感影像，进行相对辐射校正和几何精校正；计算多时相Landsat8遥感影像的NDVI图像；将多时相Landsat8遥感影像的绿、红、近红外波段时间序列与NDVI时间序列组成多特征时间序列；对多特征时间序列进行基于CD-DTW距离的K均值聚类，得到聚类结果；对聚类结果进行类别合并，得到最终的分类结果。本发明通过多特征的优化组合以及对DTW距离的改进能有效地实现时间序列遥感影像土地覆盖分类，并且能够消除云等尖锐噪声对时间序列遥感影像分类的影响。

Description

一种基于CD-DTW距离的时间序列遥感影像土地覆盖分类方法

技术领域

本发明属于遥感影像处理技术领域，主要涉及土地覆盖分类技术，具体涉及一种基于CD-DTW距离的时间序列遥感影像土地覆盖分类方法。

背景技术

土地覆盖是指地球表层的植被覆盖物和人工覆盖物的总称，是自然植被与被自然营造体和人工建造所覆盖的地表诸要素的综合反映。土地覆盖是人类认识自然、掌握自然规律必备的信息，也是各种资源管理和地理信息服务所需要的最基本数据。因此，土地覆盖信息的获取、分析和更新显得极其重要。

遥感影像数据以其宏观性、实时性的特点，一直以来都是土地覆盖检测的重要手段。遥感数据的优势在于其包含了丰富的空间信息，有利于研究地物的空间特性。但是地表事件的发生和地物的演化，是随着时间的推进而进行的，在时间和空间上都表现出了一定的变化规律。因此，针对遥感影像的分类不能只依靠地物的空间特性，应该采取合适的方法对遥感影像地物的时间特性进行分析，力求从遥感影像 的时间特性中提取与之相关的信息和知识。随着遥感传感器资源的日益丰富，遥感数据的时间分辨率有了显著提高，同一区域长期积累的遥感影像数据中蕴含着丰富的时间、空间和光谱信息，为土地覆盖动态监测提供了优良的数据源。因此，基于长期积累的遥感数据，建立遥感影像时间序列数据集并进行土地覆盖分类成为了遥感影像分类的一种发展趋势。

目前，在其它学科，尤其是金融和生物领域，对序列数据的表达、处理和挖掘方法已经相对成熟，这些为遥感影像时间序列的研究提供了坚实的基础。然而，遥感影像时间序列具有一定的特殊性，比如采样间隔不等，固有的时间轴畸变(物候期根据每年的温度等自然条件有所提前或者延迟)，存在云遮挡现象等，使得现有的时间序列处理方法不能完全适应遥感应用的需要，因此有必要对传统的时间序列处理方法进行改进，使之适应时间序列遥感影像分类的要求。

发明内容

(一)发明目的

针对面向时间序列遥感影像的土地覆盖分类，提供一种基于CD-DTW距离的K均值聚类方法。并且通过多特征的优化组合提高时间序列遥感影像分类的精度，通过对DTW距离的改进消除云等尖锐噪声对时间序列遥感影像分类的影响。

(二)技术解决方案

本发明提供了一种基于CD-DTW距离的时间序列遥感影像土地覆盖分类方法，包含以下步骤：

步骤1、对多时相的Landsat8遥感影像进行相对辐射校正和几何精校正；

步骤2、利用经过所述步骤1处理后的影像的多光谱波段，计算不同时相影像的NDVI图像。

步骤3、将多时相的Landsat8遥感影像的绿、红、近红外波段时间序列以及经过步骤2处理得到的NDVI时间序列组成多特征时间序列。

步骤4、将经过所述步骤3处理后得到的多特征时间序列进行基于CD-DTW距离的K均值聚类，得到聚类结果。

步骤5、将经过所述步骤4得到的聚类结果进行合并，生成最终的分类图像。

所述步骤2中的NDVI能够很好的区分植被与非植被，其计算公式如下：

$\mathrm{NDVI}=\frac{{\mathrm{\ρ}}_{\mathrm{nir}}-{\mathrm{\ρ}}_{\mathrm{red}}}{{\mathrm{\ρ}}_{\mathrm{nir}}+{\mathrm{\ρ}}_{\mathrm{red}}}$

其中，ρ_nir和ρ_red分别代表某个像素近红波段和红波段的像素值。

所述步骤3中将不同时相的Landsat8遥感影像的绿、红、近红外波段时间序列以及经过步骤2处理得到的NDVI时间序列组成多特征时间序列。遥感影像的绿、红、近红外波段包含地物信息丰富，能够很好地分辨各种地物，遥感图像的NDVI值能够很好的识别出水域、植被、岩石等地物，因此选择这4种波段作为时间序列的特征。

所述步骤4中的DTW(Dynamic Time Warping)是把时间规整和距离测度计算结合起来的一种非线性规整算法，该算法基于动态规划的思想，将一个复杂的全局最优化问题转化为许多局部最优化问题，一步一步地进行决策。DTW算法采用满足一定条件的时间规整函数描述输入序列和参考序列在时间上的对应关系，求解两者匹配时累计距离最小所对应的规整函数。求解最优路径的过程可以通过构造累计距离矩阵实现，累计距离矩阵计算如下：

假设存在两条序列A＝(a₁，a₂，...，a_m)，B＝(b₁，b₂，...，b_n)，则A和B的累计距离矩阵可以通过如下公式计算得到：

$D(A_i,B_j)=\mathrm{\δ}(a_i,b_j)+\min \left\{\begin{array}{c}D(A_{i-1},B_{j-1})\\ D(A_{i-1},B_j)\\ D(A_i,B_{j-1})\end{array}\right.$

$\mathrm{\δ}(a_i,b_j)=\sqrt[p]{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{\left|\right|a_{\mathrm{ik}}-b_{\mathrm{jk}}\left|\right|}_p}$

其中，D(A_i，B_j)代表序列A的第i个时间点与序列B的第j个时间点之间的DTW距离值，D为m×n矩阵，当要处理的时间序列有多个特征时，a_i，b_j为多维向量，n代表特征个数。||·||_p代表p范数，当p＝1时代表曼哈顿距离，当p＝2时代表欧式距离。

得到累计距离矩阵之后，根据每个点的min()的局部选择过程，从累计距离矩阵的结束点回溯到起始点即可得到一条最优规整路径，这条规整路径对应的就是两条序列时间点的最佳对应关系。并且，最后一个点的DTW值D(A_m，B_n)受全局时间点的影响，可以作为衡量两 条序列相似性的度量值。

本发明将兰氏距离(Canberra Distance)与动态时间规整(Dynamic Time Warping)算法相结合作为不同像素时间序列的相似性度量，本发明称之为CD-DTW距离，两条序列的CD-DTW距离的计算公式如下：

$D(A_i,B_j)=\mathrm{\δ}(a_i,b_j)+\min \left\{\begin{array}{c}D(A_{i-1},B_{j-1})\\ D(A_{i-1},B_j)\\ D(A_i,B_{j-1})\end{array}\right.$

$\mathrm{\δ}(a_i,b_j)=\underset{k}{\overset{n}{\mathrm{\Σ}}}\frac{|a_{\mathrm{ik}}-b_{\mathrm{jk}}|}{\left|a_{\mathrm{ik}}\right|+\left|b_{\mathrm{jk}}\right|}$

其中，|a_ik|是序列A的第i个时间点的第k个特征的绝对值。

兰氏距离是曼哈顿距离的一种变形，可以看作带有权重的曼哈顿距离。兰氏距离具有不受量纲影响的特点，在时间序列的度量中能够消除尖锐噪声对序列相似性的度量。因此，本发明结合兰氏距离与DTW算法的优点，提出一种新的时间序列相似性度量标准，并将其应用在时间序列遥感影像的土地覆盖分类中。通过CD-DTW距离的使用，在发挥DTW比较时间序列相似性的优势的同时，很好地解决了时间序列遥感影像中云等尖锐噪声对土地覆盖分类结果的影响。

所述步骤4中非监督聚类采用K均值聚类，其具体计算步骤如下：

1.选取K个类别的初始中心序列初始中心 序列的选择对聚类结果有一定的影响，初始中心序列的选择一般有以下两种方法：

a)任选K个类别的初始中心序列；

b)将全部数据随机地分为K个类别，每个类别随机选择一条序列作为中心序列，利用DBA算法更新该中心序列，将更新后的中心序列作为K个类的初始中心序列。

2.在第n次迭代中，对任一时间序列X按如下的方法把它划分到K个类别中的某一类。对于所有的i≠j，i＝1，2，...，K，如果 则其中是以为中心序列的类。

3.由第2步利用DBA算法得到类新的类中心

4.对于所有的j＝1，2，...，K，如果则迭代结束，否则转到第二步继续进行迭代，θ为阈值。

所述步骤4中的DBA(DTW Barycenter Averaging)算法是根据中心序列与同一类中其他序列的时间点对应关系，将与中心序列的时间点相对应的所有序列时间点值求平均值，通过多次迭代更新中心序列。每次迭代过程包含以下具体步骤：

1.计算中心序列与同一类中其他序列的DTW值，记录中心序列的时间点与其他序列的时间点的对应关系。

2.利用第1步得到的对应关系更新中心序列，假设第j个中心序列有m个时间点对于新的类中心序列 的第k个时间点

$z_{\mathrm{jk}}^{(l+1)}=\mathrm{barycenter}(x_{k1},x_{k2},...,x_{\mathrm{k\α}})=\frac{x_{k1}+x_{k2}+...+x_{\mathrm{k\α}}}{\mathrm{\α}}$

式中，l为DBA算法的迭代次数，{x_k1，x_k2，...，x_kα}为与中心序列的第k个时间点相对应的其他序列的时间点值，α为所有相对应时间点的数目。

所述步骤5中将经过所述步骤4得到的聚类结果进行合并，通过计算类中心序列相互之间的距离，合并明显属于同一类的聚类，并通过目视解译给聚类图斑赋予正确的地物类型，生成最终的分类图像。

(三)技术效果

本发明与现有的技术方案相比具有如下的优点和有益效果：本发明选取Landsat8的绿、红、近红外波段与NDVI波段组合用于遥感影像聚类，包含地类信息丰富，聚类效果好；本发明采用CD-DTW距离作为遥感影像时间序列间的相似性度量值，使序列的相似性比较更加准确，兰氏距离和DTW距离的结合消除了序列中的尖锐时间点的影响，使时间序列遥感影像的聚类不受云等尖锐噪声的影响，聚类效果比基于传统DTW距离的聚类要好。

附图说明

图1是本发明实施例的基于CD-DTW距离的时间序列遥感影像土地覆盖分类方法流程图；

图2是实验影像中的一幅(北京昌平区的Landsat8影像，拍摄于2013年6月13日)；

图3是实验影像中的一幅(北京昌平区的Landsat8影像，拍摄于 2013年11月4日)；

图4是实验影像中的一幅NDVI图像(源图像为北京昌平区的Landsat8影像，拍摄于2013年11月4日)；

图5是基于CD-DTW距离的时间序列遥感影像K均值聚类结果；

图6是将聚类结果合并后最终的分类图像。

具体实施方式

本发明以9幅北京市昌平区部分区域的Landsat8遥感影像为例，说明时间序列遥感影像地物分类的具体实施方式。9景实验影像的拍摄时间分别为2013年6月13日、2013年7月31日、2013年9月1日、2013年10月3日、2013年10月19日、2013年11月4日、2013年11月20日、2013年12月6日、2013年12月22日，实验影像中两幅源图像如图2、图3所示。下面结合附图对本发明进一步说明。

如图1所示，是本发明实施例的基于CD-DTW距离的时间序列遥感影像土地覆盖分类方法流程图，本实施例包括如下步骤：

步骤1、对多时相的Landsat8遥感影像进行相对辐射校正和几何精校正；

步骤2、利用经过所述步骤1处理后的影像的多光谱波段，计算所有影像的NDVI图像。NDVI能够很好的区分植被与非植被，其计算公式如下：

$\mathrm{NDVI}=\frac{{\mathrm{\ρ}}_{\mathrm{nir}}-{\mathrm{\ρ}}_{\mathrm{red}}}{{\mathrm{\ρ}}_{\mathrm{nir}}+{\mathrm{\ρ}}_{\mathrm{red}}}$

其中，ρ_nir和ρ_red分别代表某个像素近红波段和红波段的像素值。

图4是实验影像中的一幅NDVI图像；

步骤3、将多时相的Landsat8遥感影像的绿、红、近红外波段时间序列以及经过步骤2处理得到的NDVI时间序列组成多特征时间序列。

本发明实施例中选取9幅遥感影像组成时间序列，存储待处理时间序列的数组大小为337*340*4*9，其中337为遥感影像的行数，340为遥感影像的列数，4为时间序列的特征，分别为绿波段、红波段、近红波段、NDVI值，9为时间序列的长度，说明有9幅影像。

步骤4、将经过所述步骤3处理后得到的多特征时间序列进行基于CD-DTW距离的K均值聚类，得到聚类结果。

4.1：兰氏距离是曼哈顿距离的一种变形，可以看作带有权重的曼哈顿距离。兰氏距离具有不受量纲影响的特点，在时间序列的度量中能够消除尖锐噪声对序列相似性度量的影响。因此，本发明结合兰氏距离与DTW算法，提出一种新的时间序列相似性度量标准，本发明称之为CD-DTW距离，CD-DTW距离的计算公式如下：

$D(A_i,B_j)=\mathrm{\δ}(a_i,b_j)+\min \left\{\begin{array}{c}D(A_{i-1},B_{j-1})\\ D(A_{i-1},B_j)\\ D(A_i,B_{j-1})\end{array}\right.$

$\mathrm{\δ}(a_i,b_j)=\underset{k}{\overset{n}{\mathrm{\Σ}}}\frac{|a_{\mathrm{ik}}-b_{\mathrm{jk}}|}{\left|a_{\mathrm{ik}}\right|+\left|b_{\mathrm{jk}}\right|}$

其中，序列A＝(a₁，a₂，...，a_m)，B＝(b₁，b₂，...，b_n)，D(A_i，B_j)代表累计距离矩阵中序列A的第i个时间点与序列B的第j个时间点之间的CD-DTW距离值，D为m×n矩阵，当要处理的时间序列有多个特征时，a_i，b_j为多维向量，n代表特征个数，|a_ik|是序列A的第i个时间点的第k个特征的绝对值。累计距离矩阵结束点D(A_m，B_n)受两条序列所有时间点的影响，其值可以作为衡量两条序列相似性的度量值。

在本发明实施例中，共有9幅遥感影像参与分类，其特征分别是绿、红、近红外波段和NDVI值，因此，时间序列的长度为9，a_i，b_j分别为4维向量。

4.2：非监督聚类采用K均值算法，其具体计算步骤如下：

1.选取K个类别的初始中心序列初始中心序列的选择对聚类结果有一定的影响，初始中心序列的选择一般有以下两种方法：

a)任选K个类别的初始中心序列；

b)将全部数据随机地分为K个类别，每个类别随机选择一条序列作为中心序列，利用DBA算法更新该中心序列，将更新后的中心序列作为K个类的初始中心序列。

本发明实施例需要进行15类聚类，因此选取K＝15，并任选15条序列作为初始聚类中心序列。

2.在第n次迭代中，对任一样本X按如下的方法把它划分到K个类别中的某一类。对于所有的i≠j，i＝1，2，...，K，如果 则其中是以为中心序列的类。

3.由第二步利用DBA算法得到类新的类中心

4.对于所有的j＝1，2，...，K，如果则迭代结束，否则转到第二步继续进行迭代，θ为阈值。

本发明实施例中，令p＝2，即欧式距离，阈值θ＝10。

4.3：DBA算法是根据中心序列与同一类中其他序列的时间点对应关系，将与中心序列的时间点相对应的所有序列时间点值求平均值，通过多次迭代更新中心序列。每次迭代过程包含以下具体步骤：

1.计算中心序列与同一类中其他序列的DTW值，根据动态规整路径记录中心序列的时间点与其他序列的时间点的对应关系。

2.利用第1步得到的对应关系更新中心序列，假设第j个中心序列有m个时间点对于新的类中心序列 的第k个时间点

$z_{\mathrm{jk}}^{(l+1)}=\mathrm{barycenter}(x_{k1},x_{k2},...,x_{\mathrm{k\α}})=\frac{x_{k1}+x_{k2}+...+x_{\mathrm{k\α}}}{\mathrm{\α}}$

式中，l为DBA算法的迭代次数，{x_k1，x_k2，...，x_kα}为与中心序列j的第k个时间点相对应的其他序列的时间点值，α为所有相对应时间点的数目。

本发明实施例中，DBA算法迭代10次更新中心序列。

图5是基于CD-DTW距离的时间序列遥感影像K均值聚类结果，每一类随机赋色，总共有15类。

步骤5、将经过所述步骤4得到的聚类结果进行合并，通过计算类中心序列相互之间的距离，合并明显属于同一类的聚类，并通过目视解译给聚类图斑赋予正确的地物类型，生成最终的分类图像。

本发明实施例根据遥感影像目视解译结果将类别划分成林地、耕地、草地、建设用地、水域、其他用地6大类。

图6是将聚类结果合并后生成的最终的分类图像。

实验结果表明，通过本技术方案，对时间序列遥感影像的土地覆盖分类能够取得较理想的结果。通过多特征的优化组合，提高了时间序列遥感影像土地覆盖分类的精度；通过对DTW的改进，能够消除时间序列中尖锐时间点的影响，相比基于传统DTW距离的聚类能够取得更好的聚类效果。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变型，这些改进和变型也应视为本发明的保护范围。

Claims (3)

1.一种基于CD-DTW距离的时间序列遥感影像土地覆盖分类方法，其特征在于包含以下步骤：

步骤1、对多时相的Landsat8遥感影像进行相对辐射校正和几何精校正；

步骤2、利用经过所述步骤1处理后的影像的多光谱波段，计算不同时相影像的NDVI图像。

步骤3、将多时相的Landsat8遥感影像的绿、红、近红外波段时间序列以及经过步骤2处理得到的NDVI时间序列组成多特征时间序列。

步骤4、将经过所述步骤3处理后得到的多特征时间序列进行基于CD-DTW距离的K均值聚类，得到聚类结果。

步骤5、将经过所述步骤4得到的聚类结果进行合并，生成最终的分类图像。

2.如权利要求1所述的基于CD-DTW距离的时间序列遥感影像土地覆盖分类方法，其特征在于，所述步骤3中，提取多光谱影像的绿、红、近红外波段时间序列和NDVI时间序列组成多特征时间序列，通过多特征的优化组合，提高了时间序列遥感影像分类的精度。

3.如权利要求1所述的基于CD-DTW距离的时间序列遥感影像土地覆盖分类方法，其特征在于，所述步骤4中的CD-DTW是指将兰氏距离(Canberra Distance)与动态时间规整(Dynamic Time Warping)算法相结合的一种时间序列相似度计算方法，利用兰氏距离对尖锐噪声不敏感的特性消除时间序列中尖锐时间点对DTW距离计算影响过大的缺点。通过CD-DTW距离的使用，在发挥DTW比较时间序列相似性的优势的同时，很好地解决了时间序列遥感影像中云等尖锐噪声对分类结果的影响。