CN104318270A - 一种基于modis时间序列数据的土地覆盖分类方法 - Google Patents

Abstract

一种基于MODIS时间序列数据的土地覆盖分类方法，本发明涉及土地覆盖分类领域，本发明要解决传统方法用时长、植被指数的负偏差以及SG重建结果准确性降低的问题，而提出的一种基于MODIS时间序列数据的土地覆盖分类方法，该方法具体是按照以下步骤进行的：1、建立原始曲线；2、对原始曲线进行滤波拟合成初始曲线；3、建立初始曲线像元的无云影像二维数组；4、设置为阈值T，其中，Y_i≠y_i；5、处理过的原始曲线；6、得到重建后的NDVI年变化曲线；7、提取植被生长季参数组成特征影像；8、决定最终投票分类结果等步骤进行的；本发明应用于基于MODIS时间序列数据的土地覆盖分类领域。

Description

一种基于MODIS时间序列数据的土地覆盖分类方法

技术领域

本发明涉及土地覆盖领域，特别涉及基于MODIS时间序列数据的土地覆盖分类方法领域；

背景技术

目前利用数据统计理论方法结合人工解译仍为在大尺度内进行遥感分类的主导方法。显然这种方法具有算法成熟、充分利用人机交互和影响等特点，然而其用时长，对参与解译分析的人员依赖性强，很大程度上不具备可重复性等。这些局限性影响了迅速、准确、客观地获取大面积土地覆盖类型信息。

尽管SG(Savitzky和Golay的滤波方法)在拟合过程中一定程度上较为客观的反应真实地物的归一化植被指数(NDVI)值，但在重建过程中仍然存在两个主要问题：即：1、由于大气影响通常引起植被指数的负偏差，所以均匀的权重分布对于年际动态变化应用是不适合的，经试验与分析，曲线下面的点应该比上面的点得到更小的权重。2、经过SG重建结果低于上包络线，使其峰值降低，准确性降低。

发明内容

本发明的目的是为了解决传统方法用时长、植被指数的负偏差以及SG重建结果准确性降低的问题，而提出的一种基于MODIS时间序列数据的土地覆盖分类方法。

上述的发明目的是通过以下技术方案实现的：

步骤一、将一年中的原始MODIS NDVI时间序列影像中的无云影像设置为n+1景，儒略日为X，NDVI值为Y，建立了一个二维数组即(X₀，Y₀)，(X₁，Y₁)，…(X_n，Y_n)即为原始曲线；

步骤二、采用C5科学数据集中的VI质量评价数据QA来设置对应像素的权重，利用该权重采用SG方法对原始曲线进行滤波拟合成初始曲线；其中，C5第五代MODIS植物指数科学数据集；SG的全称为Savitzky-Golay平滑滤波器；

步骤三、将初始曲线的一个像元的无云影像设置为n+1景，儒略日为x，NDVI值为y，建立了一个初始曲线像元的无云影像二维数组(x₀，y₀)，(x₁，y₁)，…(x_n，y_n)；

步骤四、将原始曲线的峰值与初始曲线峰值的绝对值差设置为阈值T＝min{(Y₀-y₀),(Y₁-y₁),...,(Y_i-y_i)}，其中，Y_i≠y_i；

步骤五、若M超过该阈值T(M>T)，则原始曲线上的点被初始曲线上的点所取代；若M<T，则保留原始曲线的点；其中，M＝|Y_i-y_i|；

步骤六、根据阈值长度设置权重采用B样条曲线对经过步骤五处理过的原始曲线进行拟合，然后转到步骤四计算新阈值，依次带入步骤五进行计算，直到前后两次经过B样条曲线拟合的曲线无区别之后则停止计算，完成对MODIS NDVI时间序列重建得到重建后的NDVI年变化曲线；

步骤七、采用和Eklundh提出的动态阈值法，对步骤六中得到的重建后的NDVI年变化曲线中提取植被生长季参数组成特征影像；其中，光谱信息为红波段、近红外波段、重建后NDVI的最大值、最小值、平均值和标准偏差；

步骤八、在特征影像上随机采样并以谷歌地球为参考影像，根据制定的分类体系进行目视解译获得训练样本集即根据提取的生长季参数作为分类器的输入参数，利用随机森林方法的分类算法，将训练样本即提取的生长季参数构建决策树分类器得到每个决策树的分类结果，根据每个决策树的分类结果投票，决定最终投票分类结果；其中，随机森林(Random Forests)是一种基于分类与回归决策树的组合分类算法；即完成了一种基于MODIS时间序列数据的土地覆盖分类方法。

发明效果

本发明建立一套宏观尺度土地覆盖分类方法模式，以MODIS的归一化植被指数(NDVI)时间序列为主要数据源，从重建后的MODIS NDVI时间序列中反演物候特征作为主要特征，并与光谱特征信息以及坡度坡向信息结合后参与分类，采用随机森林组合分类器的方法对进行宏观尺度土地覆盖分类研究。挖掘MODIS遥感数据的优势，为宏观土地覆盖类型监测及土地资源调查提供服务。

本发明建立的一套基于MODIS时间序列数据的宏观尺度土地覆盖制图的方法流程。采用6组特征集的分类结果。利用单独的物候特征集比NDVI统计特征集的总体分类精度高5.6％，说明物候特征对于土地覆盖分类的识别度高于NDVI，且具有更实际的物理意义；但其分类精度比加入下表的第3组光谱特征的NDVI特征集低2％，在此基础上加入物候特征时精度又提高了0.8％。考虑到地形特征时，加入下表的第6组物候特征比下表的第5组特征集的总体精度高0.5％。从分析可得，利用下表的第6组特征集进行分类能达到最佳的分类效果，可见物候特征在土地覆盖类型分类中具有重要的指导作用。同时，从表中可以得出，无混合类的总体精度比有混合类的精度要高，且在总体精度最高为下表的第6组数据集中，无混合植物比有混合植物别的精度高2.3％。但考虑到MODIS数据本身混合像元较严重的问题，同时希望本研究的土地覆盖制图方法有很好的推广型。所以，本研究在分类系统中定义的混合植物是有必要的；如以下特征集评价表所示：

图2～6是林地样本的原始序列曲线和五种重构方法重构NDVI序列的结果。从图2～6为林地重构前后的对比分析图，可得出中，时间序列谐波分析法(HANTS)对林地滤波后使得曲线的结果整体偏低，和上包络线偏离较远。与HANTS相比，非对称性高斯函数拟合法(AG)、双Logistic函数拟合法(DL)对于林地的时间序列重建的相对较好，但没有峰值的逼近度不够且峰值的位置偏移，这对于后期物候特征的反演有一定程度的影响。之所以会出现这样的现象，是因为一些强噪声使得算法也很容易将其误判为提取目标而获得虚假区间，使得局部拟合出现偏移现象。SG的方法模拟较好的同时，使得峰值几乎没有位置偏移，提高后期物候特征提取的精度如图5，但仍然对于峰值没有较好的逼近，这样导致峰值信息以及相关信息偏低。使用基于加权样条曲线的SG方法，较好的逼近了包络线，且峰值与峰值的位置都有较好的保真性如图6。

附图说明

图1是具体实施方式一提出的一种基于MODIS时间序列数据的土地覆盖分类方法流程图；

图2是具体实施方式一提出的采用原始处理方法与HANTS处理方法在不同儒略日下得到NDVI的曲线对比图；其中，HANTS处理方法为时间序列谐波分析法；；

图3是具体实施方式一提出的采用原始处理方法与AG处理方法在不同儒略日下得到NDVI的曲线对比图；其中，AG处理方法为非对称性高斯函数拟合法；

图4是具体实施方式一提出的采用原始处理方法与DL处理方法在不同儒略日下得到NDVI的曲线对比图；其中，DL处理方法为双Logistic函数拟合法；

图5是具体实施方式一提出的采用原始处理方法与SG滤波处理方法在不同儒略日下得到NDVI的曲线对比图；其中，SG处理方法为Savitzky-Golay平滑滤波器；

图6是具体实施方式一提出的采用原始处理方法与BWISG处理方法在不同儒略日下得到NDVI的曲线对比图；其中，BWSG为本发明方法；

图7是具体实施方式一提出的基于NDVI(归一化植被指数)时间序列曲线提取物候指标原理图；

图8是具体实施方式一提出的随机森林原理示意图；

图9是实施例一提出的河北省2010年9月2日MODIS假彩色合成影像图；

图10是实施例一提出的河北省坡度数据示意图；

图11是实施例一提出的河北省坡向数据示意图；

图12是实施例一提出的本发明分类体系基本地物八大类型示意图；其中，a为水田类型，b为旱地类型，c为林地类型，d为草地类型，e为水体类型，f为建筑用地类型，g为裸地类型，h为混合地物类型；

图13(a)是实施例一提出的应用本发明方法对河北省土地覆盖分类的制图结果示意图；

图13(b)是实施例一提出的MCD12Q1土地覆盖产品制图结果示意图；

图14是实施例一提出的本发明制图方法与MCD12Q1土地覆盖类型产品精度对比示意图；

图15是实施例一提出的随机森林(Random forest)分类精度示意图；

图16是实施例一提出的MCD12Q1土地覆盖产品精度示意图；

图17具体实施方式一提出的基于加权样条曲线拟合的S-G滤波算法流程图。

具体实施方式

具体实施方式一：本实施方式的一种基于MODIS时间序列数据的土地覆盖分类方法，具体是按照以下步骤制备的：

步骤一、将一年中的原始MODIS NDVI时间序列影像中的无云影像设置为n+1景，儒略日为X，NDVI值为Y，建立了(无云影像的)一个二维数组即(X₀，Y₀)，(X₁，Y₁)，…(X_n，Y_n)即为原始曲线；

步骤二、采用C5科学数据集中的VI质量评价数据QA来设置对应像素的权重，利用该权重采用SG方法对原始曲线进行滤波拟合成初始曲线；其中，C5第五代MODIS植物指数科学数据集；SG的全称为Savitzky-Golay平滑滤波器；

步骤三、将初始曲线的一个像元的无云影像设置为n+1景，儒略日为x，NDVI值为y，建立了一个初始曲线像元的无云影像二维数组(x₀，y₀)，(x₁，y₁)，…(x_n，y_n)；

步骤四、将原始曲线的峰值与初始曲线峰值的绝对值差设置为阈值T＝min{(Y₀-y₀),(Y₁-y₁),...,(Y_i-y_i)}，其中，Y_i≠y_i；

步骤五、若M超过该阈值T(M>T)，则原始曲线上的点被初始曲线上的点所取代；若M<T，则保留原始曲线的点；其中，M＝|Y_i-y_i|；

步骤六、根据阈值长度设置百分比作为对应点的权重采用B样条曲线对经过步骤五处理过的原始曲线进行拟合，然后转到步骤四计算新阈值，依次带入步骤五进行计算，直到前后两次经过B样条曲线拟合的曲线无明显区别之后则停止计算，完成对MODISNDVI时间序列重建得到重建后的NDVI年变化曲线如图7即通过这种反复进行的计算过程实现图像的重构形成基于加权样条曲线拟合的Savitzky-Golay平滑滤波方法；

本发明基于加权样条曲线拟合的Savitzky-Golay平滑滤波方法；其基本思想为：在基于加权SG滤波后，对被忽略的真值点，采用样条曲线并加入已知真值点对已经滤波结果进行修正；发明提出改进的SG方法的优点在于，克服了SG滤波结果低于包络线的问题，同时，充分利用MODIS数据的优势，采用QA评价数据提高精度；同时，自适应的选择阈值进行计算，使重构的曲线逼近真值，具体流程如图17所示；

步骤七、采用和Eklundh提出的动态阈值法，对步骤六中得到的重建后的NDVI年变化曲线中提取植被生长季参数即提取物候特征；将物候特征、光谱信息以及坡度坡向辅助信息共25维特征向量组成特征影像(表1)；其中，光谱信息为红波段、近红外波段、重建后NDVI的最大值、最小值、平均值和标准偏差；

表1特征组合

特征序号	特征属性	特征序号	特征属性
				1	重构后NDVI的最大值	14	生长季结束NDVI
2	重构后NDVI的最小值	15	整个生长季的积累值
				3	重构后NDVI的平均值	16	基线值
4	重构后NDVI的标准差	17	整个生长季NDVI的中值
				5	红波段的最大值	18	生长季NDVI的峰值
6	红波段的最小值	19	峰值与基线值之间的差值

7	红波段的平均值	20	生长季开始时的比率
				8	红波段的标准差	21	生长季结束时的比率
9	近红外波段的最大值	22	面积整体的积累值
				10	近红外波段的最小值	23	峰值与基线的积累值
11	近红外波段的平均值	24	坡度
				12	近红外波段的标准差	25	坡向
13	生长季开始期NDVI的值

步骤八、在特征影像上随机采样并以谷歌地球为参考影像，根据制定的分类体系进行目视解译获得训练样本集即根据提取的生长季参数作为分类器的输入参数，利用随机森林方法的分类算法(Random Forests)，基本思想是将训练样本即提取的生长季参数构建决策树分类器得到每个决策树的分类结果，根据每个决策树的分类结果投票，决定最终投票分类结果(如图8)从而对土地覆盖进行分类和识别；为了保证MODIS土地覆盖类型制图精度，本研究分别对特征集、样本数量和随机森林分类器进行评价与分析；其中，随机森林(Random Forests)是一种基于分类与回归决策树(Classification And Regression Tree，CART)的组合分类算法；基本思想是根据训练样本即提取的生长季参数构建决策树分类器得到每个决策树的分类结果，根据每个决策树的分类结果投票，决定最终投票分类结果即完成了一种基于MODIS时间序列数据的土地覆盖分类方法如图1。

本实施方式效果：

本实施方式建立一套宏观尺度土地覆盖分类方法模式，以MODIS的归一化植被指数(NDVI)时间序列为主要数据源，从重建后的MODIS NDVI时间序列中反演物候特征作为主要特征，并与光谱特征信息以及坡度坡向信息结合后参与分类，采用随机森林组合分类器的方法对进行宏观尺度土地覆盖分类研究。挖掘MODIS遥感数据的优势，为宏观土地覆盖类型监测及土地资源调查提供服务。

本实施方式建立的一套基于MODIS时间序列数据的宏观尺度土地覆盖制图的方法流程。采用6组特征集的分类结果。利用单独的物候特征集比NDVI统计特征集的总体分类精度高5.6％，说明物候特征对于土地覆盖分类的识别度高于NDVI，且具有更实际的物理意义；但其分类精度比加入表2的第3组光谱特征的NDVI特征集低2％，在此基础上加入物候特征时精度又提高了0.8％。考虑到地形特征时，加入表2的第6组物候特征比表2的第5组特征集的总体精度高0.5％。从分析可得，利用表2的第6组特征集进行分类能达到最佳的分类效果，可见物候特征在土地覆盖类型分类中具有重要的指导作用。同时，从表中可以得出，无混合类的总体精度比有混合类的精度要高，且在总体精度最高为表2的第6组数据集中，无混合植物比有混合植物别的精度高2.3％。但考虑到MODIS数据本身混合像元较严重的问题，同时希望本研究的土地覆盖制图方法有很好的推广型。所以，本研究在分类系统中定义的混合植物是有必要的如表2所示：

表2特征集评价

图2～6是林地样本的原始序列曲线和五种重构方法重构NDVI序列的结果。从图2～6为林地重构前后的对比分析图，可得出中，时间序列谐波分析法(HANTS)对林地滤波后使得曲线的结果整体偏低，和上包络线偏离较远。与HANTS相比，非对称性高斯函数拟合法(AG)、双Logistic函数拟合法(DL)对于林地的时间序列重建的相对较好，但没有峰值的逼近度不够且峰值的位置偏移，这对于后期物候特征的反演有一定程度的影响。之所以会出现这样的现象，是因为一些强噪声使得算法也很容易将其误判为提取目标而获得虚假区间，使得局部拟合出现偏移现象。SG的方法模拟较好的同时，使得峰值几乎没有位置偏移，提高后期物候特征提取的精度如图5，但仍然对于峰值没有较好的逼近，这样导致峰值信息以及相关信息偏低。使用基于加权样条曲线的SG方法，较好的逼近了包络线，且峰值与峰值的位置都有较好的保真性如图6。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤二中采用权重为C5科学数据集中的VI质量评价数据QA为质量总评0～3，将设定对应像素值的权重为100％、60％、20％和0，如果质量评价为0，权重就为100％；如果质量评价为1，权重就为60％；如果质量评价为2，权重就为20％；如果质量评价为3，权重就为0％。其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤六中B样条曲线表达式为

$C\left(u\right)=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}{P}_i{N}_{i,k}\left(u\right)---\left(1\right)$

其中，P_i为节点N处属于i类样本数占总样本数的频度；N_i,k(u)是调和函数，也称为基函数，按照递归公式可定义为：

$N_{i,k}\left(u\right)=\frac{(u-t_i)N_{i,k-1}\left(u\right)}{t_{i+k}-t_i}+\frac{(t_{i+k+1}-u)N_{i+1,k-1}\left(u\right)}{t_{i+k+1}-t_{i+1}}---\left(3\right)$

其中t_i是节点值，T＝[t₀,t₁,...,t_L+2k+1]构成了k次B样条函数的节点矢量，节点沿着参数轴是均匀等距分布，t_i+1-t_i＝α。其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤七中采用和Eklundh提出的动态阈值法是一种动态比值形式，即给定像元和年份的植被指数(vegetation index,VI)值与当年VI振幅的比即动态阈值；动态阈值在时间域和空间域上比绝对阈值和差值阈值以及植被指数VI值具有更好的适用性；给定像元和年份的植被指数(vegetation index,VI)值是任意给定像元和年份的植物指数；提取植被生长季参数即提取物候特征；将物候特征、光谱信息以及坡度坡向辅助信息共25维特征向量组成特征影像，其中，提取植物生长季参数(物候特征)的应用的植物物候软件TIMESAT来实现的。其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤七中提取的物候特征包括：(1)、生长季开始期NDVI的值；(2)、生长季结束时NDVI的值；(3)、在生长季开始期生长季开始期NDVI的增长比率；(4)、在生长季结束时的应该是生长季结束时的NDVI减少比率；(5)、在整个生长季中NDVI的峰值；(6)、在整个生长季中NDVI的基线值；(7)、重建后的NDVI年变化曲线的峰值与基线值之间的差值；(8)、整个生长季中的NDVI的中值；(9)、整个生长季期间NDVI的累计值，NDVI曲线到基线值之间的面积；(10)、在整个生长季过程中，NDVI的累计面积和；(11)、在整个生长季过程中，峰值与基线值之间的积分结果；生长季开始期为NDVI增长达到当年NDVI振幅20％的时刻；生长季结束时为NDVI降低到当年NDVI振幅20％的时刻；整个生长季为从植被生长季开始到结束所需要的时间；积分结果为采用基于加权样条曲线的Savitzky-Golay平滑滤波方法对MODIS NDVI时间序列重建后的曲线的峰值和基线的结果。其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：步骤八中根据提取的生长季参数作为分类器的输入参数，利用随机森林方法的分类算法(Random Forests)，基本思想是将训练样本M即提取的生长季参数构建决策树分类器得到每个决策树的分类结果，根据每个决策树的分类结果投票，决定最终投票分类结果具体过程：

(1)随机重采样bootstrap技术对全体训练样本M进行随机放回抽样M次，将此抽样过程重复N次，得到S1，S2…，Sn作为N棵决策树的训练样本；

(2)N棵决策树根据各自的训练样本，采用节点随机分裂技术从全体属性特征T中随机选取t个属性特征(t≤T)作为分裂该棵树的属性特征集；

(3)根据方差不纯度指标对t个属性特征进行建树得到节点的不纯度，其计算公式为：

$f\left(N\right)=\underset{f\&NotEqual;j}{\mathrm{\&Sigma;}}p\left({\mathrm{\&omega;}}_f\right)p\left({\mathrm{\&omega;}}_j\right)=1-\underset{j}{\mathrm{\&Sigma;}}{p}^2\left({\mathrm{\&omega;}}_j\right)---\left(4\right)$

其中：ω_j为第j类的属性；P(ω_j)为节点N处属于ω_j类样本数占总样本数的频度；f为方差不纯度；

(4)根据分支停止准则预先设定一个不纯度下降差的阈值；当分支使得节点的不纯度的下降差小于这个阈值时，停止分支；完成对N棵决策树的构建，最终N棵决策树构成一个随机森林；当进行分类时，将所有分类树的分类结果进行综合采用投票方式得到最终投票分类结果；其中，不纯度的下降差是指决策树左右两个节点；

随机森林优点主要体现为随机性和多分类器投票性；其中，随机性表现为随机选取训练样本集，目的为了扩大树的棵树、随机选取分裂属性集，为了增加每棵树之间的差异度，从而提高森林的泛化误差以及所有的树都自然生长，不进行剪枝；随机森林的理论依据为大数定律且用边缘函数最大化准则；由大数定律和树的结构可知：

$\underset{k\&RightArrow;\&infin;}{\lim }{\mathrm{PE}}^{*}=P_{\mathrm{XY}}(P_{\mathrm{\&Theta;}}(h(X,\mathrm{\&Theta;})=Y)-\underset{j\&NotEqual;Y}{\max }{P}_{\mathrm{\&Theta;}}(h(X,\mathrm{\&Theta;})=j)<0)---\left(5\right)$

其中，PE^*为泛化误差，k为森林中树的数目；公式(5)表明随着树的增加，泛化误差PE将趋向一个上界，分类精度会提高；

应用N棵分类树集(h₁(x),h₂(x),...,h_n(x))以及按照随机向量X分布获取的训练集和正确的分类向量Y，定义边缘最大化函数为：

$\mathrm{mg}(X,Y)={\mathrm{av}}_{k}I(h_k\left(X\right)=Y)-\underset{j\&NotEqual;Y}{\max }{\mathrm{av}}_{k}I(h_k\left(X\right)=j)---\left(6\right)$

其中av_k为平均得票数，I(.)为指示器函数；该边缘函数刻画了对向量X正确分类Y的平均得票数超过其它任何类得票数的程度；显然，边缘函数越大，分类精度越高；j为不正确的分类向量。其它步骤及参数与具体实施方式一至五之一相同。

采用以下实施例验证本发明的有益效果：

实施例一：

本实施例一种基于MODIS时间序列数据的土地覆盖分类方法，具体是按照以下步骤制备的：

1、数据收集与处理

(1)基础数据的收集及处理

本研究采用的MODIS数据产品包括植被指数数据(MOD13Q01)和土地覆盖类型数据(MCD12Q1)(见表3所示)，其中植被指数数据(MOD13Q1)作为土地覆盖类型制图研究的基础数据，土地覆盖产品(MCD12Q1)为土地覆盖分类结果精度评价的重要参照数据集。研究区域为河北省，MOD13Q1数据是16天合成产品，一年共23幅，河北省一期数据需要4幅拼接而成，研究区一年需92景数据；本研究主要采用MCD12Q1数据产品中IGBP分类体系下的土地覆盖类型产品。

表3MODIS数据集详细信息表

本研究采用MRT软件(MODIS Reprojection Tool)对研究区2005～2010年的552景影像批处理，将投影方式从原始的Sinusoidal投影转换为UTM_ZONE_50N(WGS84坐标系)，采用双线性内插算法将原始231.7m像元重采样成250m，最后进行镶嵌和裁切，得到河北省2005～2010年每年MODIS NDVI、EVI时间序列以及光谱反射率数据(蓝波段、红波段、近红波段、中红外波段)以及VI的像元质量和像素质量，并将结果保存为TIFF数据格式。河北省2010年9月2日MODIS假彩色(近红外、红、绿波段)合成影像示意图如图9。

(2)高程数据的收集及处理

SRTM数据为90m空间分辨率的数字高程模型(DEM)数据(httP://srtm.esi.egiar.org/)，数据格式为GeoTIFF。将其经过转投影、重采样、拼接和裁切后，得到与MOD13Q1格式相同的DEM数据，并采用ACRGIS从DEM数据中提取坡度和坡向信息(如图10和图11)。

(2)验证数据的收集及整理

本研究收集的样本数据是以Google Earth中高空间分辨率影像作为参考。Google Earth遥感数据的优势为便于解译、覆盖区域大，而且可以采用相同的规则集。本研究主要考虑样本的数量和纯度。

(1)样本数量与空间形状确定。本研究中应用Hawth’s Tools(应用抽样模拟软件)，采用简单随机抽样的方法，在研究区内随机抽取3500个像素，为了保证样本点的质量，在MOD13Q1的QA(图像像元质量)数据层中选取QA＝1的样本点。采用ARCGIS软件将样本的栅格像素点转成单位大小为250m×250m的矢量格网。

(2)样本属性确定。本研究采用Google Earth高分辨率遥感影像确定样本的属性。原因在于河北省地貌复杂多样，样本分布广，且样本的网格较大，野外调查无法实现，而中分辨率遥感影像，如：TM、HJ遥感数据等，存在光谱差异、尺度转换等问题对目视判读的结果也存在较大误差。因此，采用Google Earth的高分辨率遥感影像作为研究区的参考影像进行目视判读有利于提高其准确程度且效率较高。采用方法为，将已选定的样本格网单元叠加到Google Earth的高分辨率遥感影像中，进行目视判读，将地物类型多于70％的类别标定为该样本的地物类别。图12为本研究分类体系的土地覆盖基本类型图，红框为一个250m MODIS像素大小的范围。

2、分类体系的制定

本研究目的在于建立一套MODIS时间序列数据辅助下的宏观尺度的土地覆盖分类方法流程，使其具有一定的推广性。研究者可以根据其研究需要制定相对应的分类体系。所以，本研究以河北省为例，参考中国科学院“国家资源环境遥感宏观调查与动态研究”制定的土地资源分类系统，将其一级类定义为本研究的分类体系，并对其进行适当的扩充。同时考虑到MODIS数据存在混合像元的问题，在林地、草地或耕地交界处，在操作时无法具体判定地物类型。为了提高可操作性，将分类体系进行调整，增加定义植被混合类。此外，混合类别的优点是提高林地、草地和耕地的各类别的分类精度。

表4本研究采用的土地覆盖类型分类体系

3、精度验证

(1)混淆矩阵精度评价

随机森林(Random Forests)对试验影像进行分类，结果如图13(a)所示。采用测试样本对各个分类结果进行评价，并且对随机森林分类方法建立混淆矩阵，生成各类的总体分类精度和Kappa系数，结果如表5及图15所示。实验结果表明，随机森林分类法分类精度为84.30％，Kappa系数为0.79。

表5Random forest分类结果

(2)与MODIS土地覆盖类型产品对比分析

通过将MCD12Q1产品的分类体系统一到本研究分类体系中，对利用本研究提出的土地覆盖类型分类法所得的制图结果与MCD12Q1土地覆盖产品进行对比，结果如图13(a)和13(b)。图14为河北省本研究土地覆盖制图方法和MCD12Q1土地覆盖产品数据的方法对比。如图13(a)和13(b)所示，本研究制图结果对实验区的主要地物类型的提取较图16更为准确，没有椒盐现象，提取结果的整体景观结构与研究区基本情况相符合优于MCD12Q1土地覆盖产品数据。

采用检验样本，分别对本研究土地覆盖类型制图结果和MCD12Q1土地覆盖类型产品进行精度评价。如图14所示，本研究分类结果的分类精度为84.3％，Kappa系数为0.79，而MCD12Q1产品的分类精度为57.32％，Kappa系数为0.42。本研究制图方法比MOD12Q1产品的总体精度和Kappa系数分别提高了26.98％和0.37。

通过如图13(a)和13(b)与表5对比可知，除混合类别外，本研究的分类制图方法中每一类的用户精度均在80％以上，均优于MCD12Q1产品数据的用户精度。植被和耕地用户精度高的同时制图精度均较高，而MCD12Q1数据则波段较大，错分漏分现象较为明显。本方法考虑到对植被植物不敏感的地物类型，如：建筑用地和旱地。均取得较好的效果，但MCD12Q1产品数据在此方面效果不佳。虽说本研究方法对于水域、水田等含水量较高的地物精度比其他植被类型地物相比较低，但与MCD12Q1产品数据已经得以提高。混合类别的制图精度和用户精度均较低，原因在于受混合像元的影响，由植被类型的多样性引起的光谱和NDVI值差异较大。但由于混合类的使用，使得旱地、有林地和草地的精度分别提高了2.3％、3.5％和2.7％。基于此，本研究土地覆盖制图方法流程整体上优于MCD12Q1土地覆盖产品数据。此外，将其空间分辨率从500m提高到250m。

(3)本实施例建立的一套基于MODIS时间序列数据的宏观尺度土地覆盖制图的方法流程。采用6组特征集的分类结果。利用单独的物候特征集比NDVI统计特征集的总体分类精度高5.6％，说明物候特征对于土地覆盖分类的识别度高于NDVI，且具有更实际的物理意义；但其分类精度比加入下表的第3组光谱特征的NDVI特征集低2％，在此基础上加入物候特征时精度又提高了0.8％。考虑到地形特征时，加入下表的第6组物候特征比下表的第5组特征集的总体精度高0.5％。从分析可得，利用下表的第6组特征集进行分类能达到最佳的分类效果，可见物候特征在土地覆盖类型分类中具有重要的指导作用。同时，从表中可以得出，无混合类的总体精度比有混合类的精度要高，且在总体精度最高为下表的第6组数据集中，无混合植物比有混合植物别的精度高2.3％。但考虑到MODIS数据本身混合像元较严重的问题，同时希望本研究的土地覆盖制图方法有很好的推广型。所以，本研究在分类系统中定义的混合植物是有必要的；如以下特征集评价表所示：

整个流程可作为河北省土地利用/覆盖类型实时监测的方法参照模版，为后续河北省宏观土地覆盖类型监测及土地资源评价提供服务。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (6)

1.一种基于MODIS时间序列数据的土地覆盖分类方法，其特征在于：一种基于MODIS时间序列数据的土地覆盖分类方法具体是按照以下步骤进行的：

步骤一、将一年中的原始MODIS NDVI时间序列影像中的无云影像设置为n+1景，儒略日为X，NDVI值为Y，建立了一个二维数组即(X₀，Y₀)，(X₁，Y₁)，…(X_n，Y_n)即为原始曲线；

步骤二、采用C5科学数据集中的VI质量评价数据QA来设置对应像素的权重，利用该权重采用SG方法对原始曲线进行滤波拟合成初始曲线；其中，C5第五代MODIS植物指数科学数据集；SG的全称为Savitzky-Golay平滑滤波器；

步骤三、将初始曲线的一个像元的无云影像设置为n+1景，儒略日为x，NDVI值为y，建立了一个初始曲线像元的无云影像二维数组(x₀，y₀)，(x₁，y₁)，…(x_n，y_n)；

步骤四、将原始曲线的峰值与初始曲线峰值的绝对值差设置为阈值T＝min{(Y₀-y₀),(Y₁-y₁),...,(Y_i-y_i)}，其中，Y_i≠y_i；

步骤五、若M超过该阈值T(M>T)，则原始曲线上的点被初始曲线上的点所取代；若M<T，则保留原始曲线的点；其中，M＝|Y_i-y_i|；

步骤六、根据阈值长度设置权重采用B样条曲线对经过步骤五处理过的原始曲线进行拟合，然后转到步骤四计算新阈值，依次带入步骤五进行计算，直到前后两次经过B样条曲线拟合的曲线无区别之后则停止计算，完成对MODIS NDVI时间序列重建得到重建后的NDVI年变化曲线；

步骤七、采用和Eklundh提出的动态阈值法，对步骤六中得到的重建后的NDVI年变化曲线中提取植被生长季参数组成特征影像；其中，光谱信息为红波段、近红外波段、重建后NDVI的最大值、最小值、平均值和标准偏差；

步骤八、在特征影像上随机采样并以谷歌地球为参考影像，根据制定的分类体系进行目视解译获得训练样本集即根据提取的生长季参数作为分类器的输入参数，利用随机森林方法的分类算法，将训练样本即提取的生长季参数构建决策树分类器得到每个决策树的分类结果，根据每个决策树的分类结果投票，决定最终投票分类结果；其中，随机森林是一种基于分类与回归决策树的组合分类算法；即完成了一种基于MODIS时间序列数据的土地覆盖分类方法。

2.根据权利要求1所述一种基于MODIS时间序列数据的土地覆盖分类方法，其特征在于步骤二中采用权重为C5科学数据集中的VI质量评价数据QA为质量总评0～3。

3.根据权利要求1所述一种基于MODIS时间序列数据的土地覆盖分类方法，其特征在于步骤六中B样条曲线表达式为

$C\left(u\right)=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}{P}_i{N}_{i,k}\left(u\right)---\left(1\right)$

其中，P_i为节点N处属于i类样本数占总样本数的频度；N_i,k(u)是调和函数，也称为基函数，按照递归公式可定义为：

$N_{i,k}\left(u\right)=\frac{(u-t_i)N_{i,k-1}\left(u\right)}{t_{i+k}-t_i}+\frac{(t_{i+k+1}-u)N_{i+1,k-1}\left(u\right)}{t_{i+k+1}-t_{i+1}}---\left(3\right)$

其中t_i是节点值，T＝[t₀,t₁,...,t_L+2k+1]构成了k次B样条函数的节点矢量，节点沿着参数轴是均匀等距分布，t_i+1-t_i＝α。

4.根据权利要求1所述一种基于MODIS时间序列数据的土地覆盖分类方法，其特征在于步骤七中采用和Eklundh提出的动态阈值法是一种动态比值形式，即给定像元和年份的植被指值与当年VI振幅的比即动态阈值；给定像元和年份的植被指数值是任意给定像元和年份的植物指数；提取植被生长季参数组成特征影像，其中，提取植物生长季参数的应用的植物物候软件TIMESAT来实现的。

5.根据权利要求1所述一种基于MODIS时间序列数据的土地覆盖分类方法，其特征在于步骤七中提取的物候特征包括：(1)、生长季开始期NDVI的值；(2)、生长季结束时NDVI的值；(3)、在生长季开始期生长季开始期NDVI的增长比率；(4)、在生长季结束时的应该是生长季结束时的NDVI减少比率；(5)、在整个生长季中NDVI的峰值；(6)、在整个生长季中NDVI的基线值；(7)、重建后的NDVI年变化曲线的峰值与基线值之间的差值；(8)、整个生长季中的NDVI的中值；(9)、整个生长季期间NDVI的累计值，NDVI曲线到基线值之间的面积；(10)、在整个生长季过程中，NDVI的累计面积；(11)、在整个生长季过程中，峰值与基线值之间的积分结果；其中，生长季开始期为NDVI增长达到当年NDVI振幅20％的时刻；生长季结束时为NDVI降低到当年NDVI振幅20％的时刻；整个生长季为从植被生长季开始到结束所需要的时间；积分结果为采用基于加权样条曲线的Savitzky-Golay平滑滤波方法对MODIS NDVI时间序列重建后的曲线的峰值和基线的结果。

6.根据权利要求1所述一种基于MODIS时间序列数据的土地覆盖分类方法，其特征在于步骤八中根据提取的生长季参数作为分类器的输入参数，利用随机森林方法的分类算法将训练样本M即提取的生长季参数构建决策树分类器得到每个决策树的分类结果，根据每个决策树的分类结果投票，决定最终投票分类结果具体过程：

(1)随机重采样bootstrap技术对全体训练样本M进行随机放回抽样M次，将此抽样过程重复N次，得到S1，S2…，Sn作为N棵决策树的训练样本；

(2)N棵决策树根据各自的训练样本，采用节点随机分裂技术从全体属性特征T中随机选取t个属性特征(t≤T)作为分裂该棵树的属性特征集；

(3)根据方差不纯度指标对t个属性特征进行建树得到节点的不纯度，其计算公式为：

$f\left(N\right)=\underset{f\&NotEqual;j}{\mathrm{\&Sigma;}}p\left({\mathrm{\&omega;}}_f\right)p\left({\mathrm{\&omega;}}_j\right)=1-\underset{j}{\mathrm{\&Sigma;}}{p}^2\left({\mathrm{\&omega;}}_j\right)---\left(4\right)$

其中：ω_j为第j类的属性；P(ω_j)为节点N处属于ω_j类样本数占总样本数的频度；f为方差不纯度；

(4)根据分支停止准则预先设定一个不纯度下降差的阈值；当分支使得节点的不纯度的下降差小于这个阈值时，停止分支；当进行分类时，将所有分类树的分类结果进行综合采用投票方式得到最终投票分类结果；其中，不纯度的下降差是指决策树左右两个节点。