一种土地覆盖分类模型的构建方法

技术领域

本发明涉及信息技术与遥感交叉领域，具体来说，涉及基于MODIS遥感数据进行土地覆盖分类领域，更具体地说，涉及基于MODIS数据进行全球土地覆盖分类模型的构建方法以及全球土地覆盖分类方法。

背景技术

搭载在Terra和Aqua两颗卫星上的中分辨率成像光谱仪(MODIS)，是美国地球观测系统(EOS)计划中用于观测全球生物和物理过程的重要仪器。它具有36个中等分辨率水平(0.25-1μm)的光谱波段，每1-2天对地球表面观测一次，获取陆地和海洋温度、初级生产率、陆地表面覆盖、云、气溶胶、水汽、火情等目标的图像。MODIS第1-2波段分辨率为250m，3-7波段分辨率为500m，其他波段分辨率为1000m，MODIS数据是对全球土地进行分类所使用的重要数据。但是，现在基于MODIS数据做全球或者全国尺度的长时间序列的研究目前正面临着海量数据、处理复杂等问题。近年来，已有许多基于MODIS数据进行土地覆盖方式分类的方法和策略。2015年有研究者将不同时间节点NDVI时间序列作为输入变量，采用最大似然法，光谱角填图法和最小距离法对不同输入变量进行土地覆盖分类，将黑龙江省土地覆盖类型分为耕地、林地、草地、水域、居住用地和未利用地等6类，从而提取耕地范围。2018年又有研究者利用平滑后的NDVI时序数据进行支持向量机(SVM)分类，得到农用地等分类信息。但是以上的方法一般适用某一区域，分类精度不高，在实际制作全球农地分类时各区域间的光谱数据信息存在差异较大、分类更细等特点，因此上述方法在全球细粒度农地分类图的制作上存在一定的局限性，处理实际数据时往往存在比较严重的错分情况，与此同时，由于MODIS数据量的庞大，现有技术无法解决高效处理海量数据的问题，设计适用于海量MODIS数据的并行化处理方法也是急需解决的问题之一。

发明内容

因此，本发明的目的在于克服上述现有技术的缺陷，提供一种新的土地覆盖分类模型的构建方法以及基于构建的模型对土地进行分类的方法。

根据本发明的第一方面，提供一种土地覆盖分类模型的构建方法，包括：S1、获取土地历史数据，对其进行聚类处理以形成多个区域，每个区域中包含多个网格点数据；S2、对每一个区域中的网格点数据提取每个网格点的指数特征数据，并将网格点指数特征数据插值到土地覆盖数据对应的站点，获得多个站点数据；S3、从每个区域的站点数据中选取样本组成该区域对应的训练集和测试集；S4、用每个区域对应的训练集训练多个随机森林模型，并用验证集验证每个随机森林模型的分类准确率，然后用每个区域对应的测试集验证训练后的该区域对应的随机森林模型的分类准确率，从中选出分类准确率最高的随机森林模型作为该区域的分类模型。

优选的，所述土地历史数据是全球土地历史数据。

在本发明的一些实施例中，在所述步骤S1中，采用AP聚类方法对土地历史数据进行聚类处理，其中，土地历史数据是MODIS数据，包括多个网格点数据，每个网格点数据包含多种波段数据，所述步骤S1包括：S11、获取土地历史数据中每个网格点的不同时刻的不同波段数据组成波段数据特征向量；S12、根据网格点的波段数据特征向量计算网格点之间的相似度矩阵，并以相似度矩阵的均值初始化每个网格点的参考度矩阵；S13、初始化网格点之间的吸引度矩阵和归属度矩阵，按照预设的迭代次数和衰减系数多次计算网格点之间的吸引度矩阵和归属度矩阵；S14、按照网格点之间最终的吸引度矩阵和归属度矩阵对网格点进行聚类，聚类后每个网格点对应一个聚类标签，相同的聚类标签对应的网格点组成一个区域。其中，所述每个网格点的不同波段数据包括归一化植被指数、增强型植被指数、植被指数、红光波段反射率、近红外波段反射率、蓝光波段反射率、中红外波段反射率、像素质量控制码。

在本发明的一些实施例中，所述步骤S2中，土地覆盖数据是联合国粮食及农业组织发布的全球土地覆盖数据，包括多个站点，站点类别包括农田、森林、草地、灌丛、水体、城市、裸地、冰雪，且所述步骤S2包括针对每一个站点进行如下操作：S21、提取每个网格点的指数特征，其中，指数特征包括：归一化植被指数、增强型植被指数、植被指数、红光波段反射率、近红外波段反射率、蓝光波段反射率、中红外波段反射率、像素质量控制码、归一化水体指数；S22、计算站点与其邻近网格点之间的欧式距离；S23、根据计算出的欧式距离计算站点每个临近网格点的权重；S24、基于每个邻近网格点的权重，将每个邻近网格点对应的每种指数特征插值到站点。其中，优选的，在所述步骤S21中，采用多种主从并行方法提取网格点的指数特征，其中，所述主从并行方法包括：主从加速并行方法、主从协同并行方法、主从异步并行方法、主从动态并行方法。

优选的，所述步骤S23中，通过如下方式将站点的邻近网格点对应的指数特征插值到站点：

F

其中，F

在本发明的一些实施例中，在所述步骤S3中，对从站点数据中选取样本进行如下加工处理：S31、进行错误样本的甄别、修复或去除；S32、对不同类别样本比重进行均匀化处理；S33、将样本的合理衍生数据参与运算处理。

在本发明的一些实施例中，在所述步骤S4中，针对每个区域执行如下操作：S41、将步骤S3中获取当前区域的训练集作为初始训练集，对其进行多次有放回的抽样，得到多个新的决策树训练集；S42、针对每个决策树训练集，计算该决策树训练集中每个类别对应的基尼系数，将类别按照基尼系数从小到大进行排列，依次将类别作为决策树的第一分类依据、第二分类依据，依次类推，以构成决策树，使得每个决策树训练集对应一个决策树最终获得多个决策树；S43、将步骤S42获得的决策树按照预设的深度范围和树的个数范围，组成多个随机森林模型，用当前区域验证集去验证每个随机森林模型的分类准确率，选出分类准确率最高的随机森林模型作为当前区域的分类模型。优选的，决策树训练集中每个类别对应的基尼系数通过如下方式计算：

其中，D代表决策树训练集，k代表决策树训练集D中的样本类别数，w

根据本发明的第二方面，提供一种土地覆盖分类方法，包括：T1、获取当前土地数据，对其进行聚类处理以形成多个区域，每个区域中包含多个网格点数据；T2、获取如本发明第一方面一所述的方法构建的土地覆盖分类模型，用每个区域对应的分类模型对步骤T1中每个区域的网格点数据进行土地覆盖分类。

与现有技术相比，本发明的优点在于：本发明能够提高分类的准确度、细粒度以及运行效率，采用本发明的方法可以将土地覆盖类型从海量MODIS数据中快速的提取出来，并实现快速、准确的土地覆盖分类。

附图说明

以下参照附图对本发明实施例作进一步说明，其中：

图1为根据本发明实施例的一种全球土地覆盖分类模型的构建方法流程示意图；

图2为根据本发明实施例的一种全球土地覆盖分类模型的构建方法中数据并行处理方案示意图；

图3为根据本发明实施例的一种全球土地覆盖分类模型的构建方法中网格点数据插值到站点的示意图；

图4为根据本发明实施例的一种全球土地覆盖分类模型的构建方法中随机森林模型构建示意图；

图5为根据本发明实施例的一种全球土地覆盖分类结果示意图。

具体实施方式

为了使本发明的目的，技术方案及优点更加清楚明白，以下通过具体实施例对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

鉴于现有技术存在的问题，本发明的目的是为了解决全球土地覆盖分类问题，同时实现全球土地海量数据的高效并行化处理，以将土地覆盖类型从海量MODIS数据中快速的提取出来，并实现快速准确的土地覆盖分类。其中，MODIS数据中的数据是标准的网格点数据，按照经纬度标准划分。

需要说明的是，本发明中土地覆盖分类采用联合国粮食及农业组织发布的全球土地覆盖数据，包括多个站点，站点类别包括农田、森林、草地、灌丛、水体、城市、裸地、冰雪，当然，实际应用中可以不局限于这些，本发明仅以这些类别举例进行说明,本发明实施例中以全球土地数据为例。

根据本发明的一个实施例，本发明提供一种土地覆盖分类模型的构建方法，如图1所示，包括步骤S1、S2、S3、S4，下面结合附图详细说明每个步骤。

在步骤S1中，进行区域聚类，通过获取全球土地历史数据，对其进行聚类处理以形成多个区域，每个区域中包含多个网格点数据。本发明中的全球土地历史数据进行聚类时可以是以MODIS遥感指数为要素进行聚类，也可以是以经纬度、降水、气温等要素进行聚类。根据本发明的一个实施例，以遥感指数为例进行聚类说明，采用AP聚类方式对全球历史原始遥感数据资料进行处理，实现特征区域划分。

首先说明一下AP聚类的原理：AP聚类算法是一种新的无监督聚类算法，聚类的目标是使数据点与其类代表点之间的距离达到最小化，因此选用欧氏距离作为相似度的测量指标，即任意两个点p(x

用代表矩阵(也叫吸引度矩阵)r(p,q)(responsiblity)和适选矩阵(也叫归属度矩阵)a(p,q)(availability)来表示数据点之间的两类信息，其中r(p,q)是从数据点p(x

AP算法的输入是N个数据点之间的相似度矩阵S，以矩阵S对角线上的数据作为该点成为聚类中心的评判标准，称之为参考度(P)，初始时候，所用点的参考度设为相同的P值(通常取S的均值)。其中r(p,q)和a(p,q)按照式如下方式进行计算：

同时为了避免振荡，AP算法更新信息时引入了衰减系数λ，每条信息被设置为它前次迭代更新的λ倍加上本次信息更新的1-λ倍，其中衰减系数为λ∈[0,1]的实数，根据所需衰减速度设置具体的衰减系数，即第t+1次更新后的吸引度矩阵r

r

a

采用AP聚类方式对本发明的全球土地历史数据进行聚类包括如下步骤：

a)获取全球土地历史数据中每个网格点的历史MODIS数据中的不同时刻的不同波段数据组成波段数据特征向量，本实施例中使用的波段主要包括归一化植被指数、增强型植被指数、植被指数、红光波段反射率、近红外波段反射率、蓝光波段反射率、中红外波段反射率、像素质量控制码等；

b)根据波段数据特征向量计算网格点之间的相似度矩阵S，以相似度矩阵S的均值初始化每个网格点的参考度矩阵；

c)根据上述公式(1)和公式(2)初始化吸引度和归属度矩阵r(p,q)和a(p,q)，设置最大迭代次数T，衰减系数λ，其中最大迭代次数和衰减系数根据实验确定，默认情况下衰减系数为0.5；

d)对于全球土地历史数据中所有的网格点[row,col]，其中row为全球网格化之后的行数(纬度范围)，col为全球网格化之后的列数(经度范围)，按照以下方式迭代计算T次：

d1、按照公式(1)和公式(2)计算r(p,q)和a(p,q)，p＝1，2，…row；q＝1,2，…col。

d2、按照公式(3)和公式(4)更新r

e)按照网格点之间最终的吸引度矩阵和归属度矩阵对网格点进行聚类，聚类后每个网格点对应一个聚类标签，相同的聚类标签对应的网格点组成一个区域，最后获得聚类后的结果图，每一个网格点都对应一个聚类标签g，g＝1…n。

聚类后，具有相同类别标签的网格点，具有相似的植被状态，不同类别网格点之间的植被状态差异较大。因为具有相同聚类标签的网格点，植被状态接近(例如华东区域、华南区域等)，针对这个类别单独训练这个类别的分类模型，可以更好的进行有监督的分类。

在步骤S2中，按照聚类的区域进行并行化数据预处理，是指对每一个区域中的网格点数据进行预处理以提取每个网格点的指数特征数据，并将网格点指数数据插值到全球土地覆盖数据对应的站点，获得多个站点数据；

根据本发明的一个实施例，所示步骤S2包括:

首先，提取每个网格点的指数特征，指数特征包括：归一化植被指数、增强型植被指数、植被指数、红光波段反射率、近红外波段反射率、蓝光波段反射率、中红外波段反射率、像素质量控制码、归一化水体指数；其中，归一化植被指数、增强型植被指数、植被指数、红光波段反射率、近红外波段反射率、蓝光波段反射率、中红外波段反射率、像素质量控制码等数据均可以直接从MODIS数据中提取得到，归一化水体指数(NDWI)需要再次计算得出，其计算公式为：

其中，NIR为MODIS数据中提取的近红外波段反射率，WIR为MODIS数据中提取的短波红光波段反射率，进一步计算水体指数的目的是为了增加特征，便于更好的区分水体和其它类别的土地覆盖。特征提取完成以后根据MODIS自带的像素质量控制码进行数据清洗，保留质量控制码为0-3的网格点。

由于MODIS遥感数据是海量数据，数据量非常庞大，为了提高数据处理的效率，根据本发明的一个实施例，设计了多种并行化处理方法，用于实现指数特征的快速提取、数据清洗和计算等操作。并行化方式为土地覆盖分类中的数据预处理提供高性能计算的多种主从并行方法，根据本发明的一个实施例，本发明通过设计主从加速并行、主从协同并行、主从异步并行、主从动态并行方法实现多粒度计算加速的提升，为遥感数据的计算和服务提供自适应的加速，如图2所示，将海量遥感数据放到超级计算机中进行并行加速计算，其中，利用主从加速并行进行MODIS遥感数据的提取，此时主核用于应用程序的通信、I/O和部分串行代码的计算，从核用于加速计算；利用主从协同并行进行NDWI数据的计算，此时主核按能力负载分配完成计算，从核按计算能力负载分配任务完成计算；利用主从异步并行进行数据的质量控制，此时主核用于应用程序的通信、I/O和部分串行代码的计算，从核用于加速计算；此外，在这里本发明还设计了一个主从动态并行化方法，主核用于任务分配，从核用于任务计算机返回计算结果，该方法主要用于后面的分类模型训练。在加速计算过程中，每条记录里面，为遥感光谱及类别标签，按照聚类的类别进行存储。并行加速方法主要是针对数据处理过程的加速，通过超级计算机来实现，此处不再赘述。

然后，对土地类别进行定义，主要是提取全球土地覆盖空间分布信息。类别的定义需要满足，同级类别之间的概念不能有重叠，全部类别应尽量覆盖研究区域内所有的对象特征。根据FAO(联合国粮食及农业组织)发布的全球土地覆盖数据，这里面类别包括农田、森林、草地、灌丛、水体、城市、裸地、冰雪。但FAO发布的土地覆盖数据是不规则的站点数据，MODIS遥感获取的数据是规则的网格点数据，在进行多遥感指数融合时候，需要把规则的格点插值到站点上，根据本发明的一个实施例，采用反距离权重的方式把格点数据和站点数据进行融合，距离越近权重越大。如图3所示，以4近邻为例将站点附近4个邻近格点插值到站点，包括:

(1)站点的设四个近邻的网格点的经纬度为(也可以用8个或16个近邻网格点)Q

(2)计算每个网格点与站点的距离

(3)计算每个网格点的权重

(4)计算网格点中每种指数特征插值到站点的值

网格点和站点融合后，站点数据不仅具有了分类，还具有了指数特征。

需要说明的是，根据实际聚类后的网格点数据特点，插值可以是4近邻、8近邻、16近邻等。

在步骤S3中，构建训练集和测试集，即从每个区域的站点数据中选取样本组成该区域对应的训练集和测试集；其中，从站点数据构建样本，基于样本构成的训练集要满足以下的条件：训练集要有代表性；训练集中不能有错误的样本；训练集要尽量完备。受数据获取精度及数据集完备性的限制，实际获得的原始数据未必都能完全满足上述条件。为了使分类结果免受质量不高数据的影响，需要对实际获得的数据做前期的加工和处理，主要包括以下方法：错误样本的甄别、修复或去除，不同类别样本比重的均匀化，样本的合理衍生数据参与运算处理。

在步骤S4中，用每个区域对应的训练集训练多个随机森林模型，用该区域对应的测试集验证训练后的随机森林模型，从中选出分类准确率最高的随机森林模型作为该区域的分类模型。由于FAO发布的土地覆盖数据，数据量很少，有效的数据只有几万条，没法实现全球的土地覆盖分类，如果想要实现全球的精细化覆盖分类如(250m分辨率，或者更高的分辨率)，需要利用已有的标签数据，建立遥感光谱信息和土地覆盖之间关系，然后把模型进行固化，利用固化的模型，对土地覆盖进行分类。为了实现快速准确的土地覆盖分类，根据本发明的一个实施例，采用随机森林模型，搭建分类模型进行训练及预测，包括：

Step1：利用前面得到的数据集，针对每个聚类的类别，训练一个分类模型，对于一个包含m个训练样本的数据集，采用有放回的抽样，组成含有m*α(其中0<α<0.8)个样本的新的决策树训练集。

Step2：重复Step1进行T遍，得到T个决策树训练集，针对每一个新构建的决策树训练集D，根据样本属性，按照公式(5)计算该决策树训练集中每个类别的基尼系数：

k代表决策树训练集D中的样本类别数，w

Step3：计算按照某属性划分后的决策树训练集的基尼系数，选择基尼系数最小的那个类别作为第一划分依据，然后继续选择第二小的属性，以此类推，完成单个决策树的构建；

Step4：在新构建的T个决策树训练集基础上上，根据step2、step3的流程，独立的训练出T个决策树，组成森林，示意如图4所示。

根据本发明的实施例，随机森林数的个数取值范围5-20，树的深度的取值范围是4-9。

通过上述步骤，可以构建多个随机森林模型。此处的随机森林指的是利用多棵树对样本进行训练并预测的一种分类器。在机器学习中，随机森林是一个包含多个决策树的分类器，并且其输出的类别是由决策树输出的类别的众数而定。样本数据很少能直接获取，往往需要将多种数据源综合运用，得到最终满足分类器训练要求的数据集。数据的准备还包括依据分类器要求，转换数据格式、进行数据格式上的准备等。将准备好的训练数据集导入训练分类器，开始分类器的训练。在训练过程中，使用主从动态并行方式进行训练，主核分配任务给从核进行计算，从核计算完损失函数和参数更新之后损失函数和准确率返回主核，主核依据所有从核的反馈结果计算某一轮训练的整体损失和准确率。将训练后的模型进行固化，并利用训练集按照如下公式来进行分类结果评价以进行模型土地覆盖的分类准确率验证：

验证完成之后，保存每个区域中准确率最高的模型参数(树的个数和深度)，得到全球土地覆盖数据的分类模型，在最终训练好的模型中输入需要分类的遥感数据就可以生成农地分类结果图。

例如，将2019年的全球MODIS遥感数据输入本发明构建的分类模型，可以得到如图5所示的全球土地覆盖分类结果图。

本发明能够提高分类的准确度、细粒度以及运行效率，采用本发明的方法可以将土地覆盖类型从海量MODIS数据中快速的提取出来，并实现快速、准确的土地覆盖分类。

需要说明的是，虽然上文按照特定顺序描述了各个步骤，但是并不意味着必须按照上述特定顺序来执行各个步骤，实际上，这些步骤中的一些可以并发执行，甚至改变顺序，只要能够实现所需要的功能即可。

本发明可以是系统、方法和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于使处理器实现本发明的各个方面的计算机可读程序指令。

计算机可读存储介质可以是保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以包括但不限于电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意合适的组合。计算机可读存储介质的更具体的例子(非穷举的列表)包括：便携式计算机盘、硬盘、随机存取存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、静态随机存取存储器(SRAM)、便携式压缩盘只读存储器(CD-ROM)、数字多功能盘(DVD)、记忆棒、软盘、机械编码设备、例如其上存储有指令的打孔卡或凹槽内凸起结构、以及上述的任意合适的组合。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。