一种土地覆被制图方法

本发明涉及土地覆盖信息领域，尤其涉及一种土地覆被制图方法。

背景技术：

土地利用与土地覆被是人类社会经济活动与自然生态过程的纽带，直接表征着人类活动对地表生态系统的影响。精确的土地覆被数据是全球变化、陆地表层过程模拟、生态文明建设及区域可持续发展等研究的重要基础数据，在土地资源可持续利用以及政府决策中发挥重要作用。遥感数据成为大范围绘制土地覆被数据的主要手段，由于已有产品在应用于区域或者全球尺度时精度不高且存在明显差异，采用数据融合方法可以有效克服目前土地覆被数据集在兼容性、可比性和准确性方面的局限性。

目前常用的数据融合方法包括地理加权回归和建立融合决策规则等。国内外学者对此有一定研究，如应用地理加权回归方法融合多种已有森林产品，提取了全球范围的森林覆盖；基于贝叶斯理论定量地描述每个像元的土地覆被类型，融合生成了多种土地覆被产品全球土地覆被图；基于证据理论建立融合决策规则，完成了整合多源产品的中国土地覆被图；基于模糊逻辑思想融合不同土地分类系统，实现了多源土地覆被产品的综合；通过建立置信图描述了不同产品的分类误差，并协同了多种土地覆被产品。

然而，基于地理加权回归方法的准确性受到训练样本的数量和采样密度的影响；而基于融合决策规则的方法需要大量先验知识和土地利用统计数据，并且原始土地覆盖产品的质量将会显著影响最终融合产品的准确性。此外，通过以上方法生成的融合产品的空间分辨率较低，难以满足土地科学研究的需求。现有土地覆被产品在遥感数据源、时空分辨率、土地覆被分类系统和分类方法等方面的一致性较差，导致不同地类的准确性差异显著。在异质性地区，尤其是混合地类，现有土地覆盖产品的准确性和一致性较差，限制了其在区域和全球尺度上的应用。目前缺乏通过数据融合的方法实现高精度土地覆被制图的技术，随着土地覆盖产品数量的增加和研究区域的扩大，现有的数据融合方法可能非常耗时且难以实施。

技术实现要素：

本发明提供了一种土地覆被制图方法，以解决目前土地覆被数据精度不高和一致性较差、数据融合算法应用在大尺度上不足的技术问题的一者或多者。

为了解决上述技术问题，本发明实施例提供了一种土地覆被制图方法，包括：

获取目标区域的多源土地覆被产品，对所述多源土地覆被产品进行重分类，得到多类土地覆被产品，以及分别对每一类土地覆被产品进行预处理，得到格式一致的土地覆被产品；

获取目标区域的遥感影像数据，对所述遥感影像数据中的低质量像元进行掩膜，提取所述遥感影像数据中的高质量多波段数据并计算植被指数、建筑指数、裸土指数、燃烧指数、雪指数和水体指数；

根据目标区域中的气候分区和行政范围边界，对目标区域进行划分，得到多个子区域；

对所述格式一致的土地覆被产品进行空间叠置处理，提取所有产品中为同一地类的栅格单元作为完全一致区域，并对完全一致区域进行面向对象的分割进一步筛选可靠区域，剩余的区域作为不一致区域；

对每个子区域中的完全一致区域根据地类进行分层随机抽样，提取初始的遥感影像解译训练样本，以及对获取的训练样本进行提纯，得到每个子区域的遥感影像解译训练样本集；

根据遥感影像多波段合成值、系列遥感指数和地形辅助数据，生成土地覆被解译特征集合，通过机器学习算法对所述每个子区域的遥感影像解译训练样本集进行构建分类器，对不一致区域进行影像解译；其中，所述系列遥感指数包括所述植被指数、建筑指数、裸土指数、燃烧指数、雪指数和水体指数；

将目标区域中不一致区域进行影像解译后的土地覆被结果与完全一致区域的土地覆被进行镶嵌，生成土地覆被修正数据。

作为优选方案，所述遥感影像数据为landsatsr影像和sentinel-1a影像，所述高质量多波段数据包括landsat5/7的第1波段至第5波段、第7波段共6个波段，sentinel-1a的vv、vh和vv-vh共3个波段，以及landsat8第2波段至第7波段共6个波段；其中，所述计算植被指数、建筑指数、裸土指数、燃烧指数、雪指数和水体指数的公式具体如下：

植被指数：ndvi＝(nir-red)/(nir+red)；

evi＝2.5×(nir-red)/(nir+6×red-7.5×blue+1)；

建筑指数：ndbi＝(mir-nir)/(mir+nir)；

裸土指数：bsi＝((swir2+red)-(nir+blue))/((swir2+red)+(nir+blue))；

燃烧指数：nbr＝(swir1-swir2)/(swir1+swir2)；

雪指数：ndsi＝(green-swir1)/(green+swir1)；

水体指数：mndwi＝(green-swir1)/(green+swir1)；

式中，blue，green，red，nir，swir1，swir2分别为landsat5tm/7etm+影像的波段1-蓝，波段2-绿，波段3-红，波段4-近红外，波段5-短波红外1，波段7-短波红外2；landsat8oli影像的波段2-蓝，波段3-绿，波段4-红，波段5-近红外，波段6-短波红外1，波段7-短波红外2。

作为优选方案，在所述对所述多源土地覆被产品进行重分类之前，包括：对每个多源土地覆被产品进行分幅的合并镶嵌和掩膜提取；

所述对所述遥感影像数据中的低质量像元进行掩膜的步骤中，包括：对所述遥感影像数据中的低质量像元进行掩膜处理，其中，所述低质量像元包括云、雪和云阴影像元；其中针对sentinel-1a影像，对每一景原始影像的散斑噪声进行滤波处理，并重新采样与landsat影像进行配准。

作为优选方案，所述分别对每一类土地覆被产品进行预处理的步骤中，包括：

对所述土地覆被产品的空间范围和参照坐标系进行统一处理；

对所述土地覆被产品的格式统一为栅格格式；

选择所述土地覆被产品中的中等空间分辨率作为标准，根据最优面积原则对土地覆被产品进行重采样，对剩余土地覆被产品的空间分辨率进行统一处理。

作为优选方案，所述根据目标区域中的气候分区和行政范围边界，对目标区域进行划分，得到多个子区域的步骤中，包括：

根据-geiger气候分区的一级类和国家行政边界，对目标区域进行空间合并提取，划分成多个子区域。

作为优选方案，所述对每个子区域中的完全一致区域根据地类进行分层随机抽样，提取初始的遥感影像解译训练样本的步骤中，包括：

对每一地类的完全一致区域根据地类使用简单非迭代聚类算法进行面向对象的分割，依据1σ准则和分割区域的多波段合成值提取各地类可靠的训练区；根据筛选后完全一致区域的数量比例进行分层随机抽样，以栅格为中心建立90m×90m的正方形样本，作为初始的遥感影像解译训练样本。

作为优选方案，所述对获取的训练样本进行提纯，得到每个子区域的遥感影像解译训练样本集的步骤中，包括：

提取训练样本内植被指数的均值和标准差，计算训练样本内植被指数的变异系数，对变异系数大于1的训练样本进行剔除，其中，变异系数＝标准差/均值；

根据从低到高排序建立直方图，在直方图首尾两端分别去除2.5％的训练样本，保留中心95％的训练样本；

其中，每个地类的训练样本数不低于1000个，若不足，则从相邻的子区域中进行随机抽取补充。

作为优选方案，所述遥感影像多波段合成值包括所述遥感影像数据在一年内4个时期的各波段中值合成及全年所有影像的波段数值标准差，其中，4个时期分别为1-3月、4-6月、7-9月、10-12月。

作为优选方案，所述的面向对象分割算法为简单非迭代聚类算法。

作为优选方案，所述地形辅助数据包括海拔和坡度数据。

作为优选方案，所述机器学习算法为随机森林算法。

相比于现有技术，本发明实施例具有如下有益效果：

1、本发明技术方案建立了多源土地覆被产品的高精度融合方法，既有效提取了已有土地覆被产品的有效信息，又对已有产品的不确定性进行了有效的修正，解决了目前多源土地覆被数据一致性程度差、精度低的问题；

2、基于空间分析和统计筛选方法定义了多源土地覆被产品的完全一致区域和不一致区域，建立了快速提取大范围可靠训练样本的方法，可以有效解决了谷歌地球引擎等云计算平台训练样本不足的难点；

3、本技术方案无需实地土地覆被解译参考样本点以及先验知识，可快速、准确地生成大范围区域较高空间分辨率的土地覆被数据融合产品，能够有效实现高精度和高效率的土地覆被制图，具有可移植性。

附图说明

图1：为本发明实施例提供的一种土地覆被制图方法的步骤流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例

请参照图1，为本发明实施例提供的一种土地覆被制图方法的步骤流程图，包括步骤101至步骤107，各步骤具体如下：

步骤101，获取目标区域的多源土地覆被产品，对所述多源土地覆被产品进行重分类，得到多类土地覆被产品，以及分别对每一类土地覆被产品进行预处理，得到格式一致的土地覆被产品。

在本实施例中，在所述对所述多源土地覆被产品进行重分类之前，包括：对每个多源土地覆被产品进行分幅的合并镶嵌和掩膜提取。

在本实施例中，所述分别对每一类土地覆被产品进行预处理的步骤中，包括：对所述土地覆被产品的空间范围和参照坐标系进行统一处理；对所述土地覆被产品的格式统一为栅格格式；选择所述土地覆被产品中的中等空间分辨率作为标准，根据最优面积原则对土地覆被产品进行重采样，对剩余土地覆被产品的空间分辨率进行统一处理。

具体地，所述的土地覆被数据可包括多地类的综合土地覆被数据和某一地类的专题土地覆被数据，所述海拔数据为srtmdem数据；所述土地覆被重分类是建立统一的分类体系，对所有土地覆被产品进行重分类；所述的空间配准预处理是统一所有产品的空间范围及参照坐标系，所述的数据格式预处理为将所有土地覆被产品格式统一为栅格格式，所述的空间分辨率统一预处理选择各产品的中等空间分辨率作为标准，依据最优面积原则进行土地覆被产品重采样，统一其余产品的数据分辨率，重采样后的栅格土地覆被类型依据该栅格单元内最大面积的地类来确定。

本实施方式中，获取覆盖研究区的多种土地覆被产品，进行各产品分幅的合并镶嵌和掩膜提取，本例使用的产品包括cgls-lc、cci-lc、from-glc和mcd12q1等综合土地覆被产品，以及耕地产品gfsad-30、林地产品palsar、水体产品gswd、城市用地产品ghs-build等专题土地覆被产品，以及获取海拔数据，并以此计算坡度。

本实施例建立的统一土地分类系统为耕地、林地、草地、灌丛、水体、城市用地、裸地、冰川积雪、苔原和湿地等，依据最优面积将土地覆被产品进行重分类。其中，统一各土地覆被产品的空间坐标系，本例选择为wgs84坐标系统，通过比较选取了各土地覆被产品分辨率居中的cci-lc产品的分辨率作为参照，对各土地覆被产品进行栅格重采样。

步骤102，获取目标区域的遥感影像数据，对所述遥感影像数据中的低质量像元进行掩膜，提取所述遥感影像数据中的高质量多波段数据并计算植被指数、建筑指数、裸土指数、燃烧指数、雪指数和水体指数。

在本实施例中，所述对所述遥感影像数据中的低质量像元进行掩膜的步骤中，包括：对所述遥感影像数据中的低质量像元进行掩膜处理，其中，所述低质量像元包括云、雪和云阴影像元。

在本实施例中，所述影像来自谷歌地球引擎平台上landsatsr影像和sentinel-1a影像，所述的低质量像元掩膜是通过landsatsr影像自带的qa波段完成；具体地，基于谷歌地球引擎平台获取landsatsr数据和海拔数据srtmdem，依据landsat影像的qa波段信息，掩膜遥感影像中的云、雪和云阴影等质量差的像元，并应用海拔数据计算坡度。其中针对sentinel-1a影像，对每一景原始影像的散斑噪声进行滤波处理，并重新采样与landsat影像进行配准。

所述高质量多波段数据包括landsat5/7的第1波段至第5波段、第7波段共6个波段，sentinel-1a的vv、vh和vv-vh共3个波段，以及landsat8第2波段至第7波段共6个波段；其中，所述计算植被指数、建筑指数、裸土指数、燃烧指数、雪指数和水体指数的公式具体如下：

植被指数：ndvi＝(nir-red)/(nir+red)；

evi＝2.5×(nir-red)/(nir+6×red-7.5×blue+1)；

建筑指数：ndbi＝(mir-nir)/(mir+nir)；

裸土指数：bsi＝((swir2+red)-(nir+blue))/((swir2+red)+(nir+blue))；

燃烧指数：nbr＝(swir1-swir2)/(swir1+swir2)；

雪指数：ndsi＝(green-swir1)/(green+swir1)；

水体指数：mndwi＝(green-swir1)/(green+swir1)；

式中，blue，green，red，nir，swir1，swir2分别为landsat5tm/7etm+影像的波段1-蓝，波段2-绿，波段3-红，波段4-近红外，波段5-短波红外1，波段7-短波红外2；landsat8oli影像的波段2-蓝，波段3-绿，波段4-红，波段5-近红外，波段6-短波红外1，波段7-短波红外2。

步骤103，根据目标区域中的气候分区和行政范围边界，对目标区域进行划分，得到多个子区域。

在本实施例中，所述根据目标区域中的气候分区和行政范围边界，对目标区域进行划分，得到多个子区域的步骤中，包括：根据-geiger气候分区的一级类和国家行政边界，对目标区域进行空间合并提取，划分成多个子区域。

具体地，分区是考虑研究区范围大、区域影像特征差异明显划定的，若研究区范围较小，则不作分区；而划定多源土地覆被产品融合与修正的子区，首先提取-geiger气候分类图的一级类，本实施例包括了赤道气候带、干燥气候带、暖温气候带、冷温气候带和极地气候带，并结合研究区的国家行政边界范围，进行空间合并提取，将气候带在个别国家过小的斑块划入国家行政边界，最终对研究区划分为11个子区域。

步骤104，对所述格式一致的土地覆被产品进行空间叠置处理，提取所有产品中为同一地类的栅格单元作为完全一致区域，并对完全一致区域进行面向对象的分割进一步筛选可靠区域，剩余的区域作为不一致区域。

具体地，对选用的所有土地覆被产品进行空间叠置，并进行一致性分析，提取所有产品均认为是同一地类的栅格单元，划分为多源土地覆被产品完全一致区域，其余的区域划分为不一致区域。

步骤105，对每个子区域中的完全一致区域根据地类进行分层随机抽样，提取初始的遥感影像解译训练样本，以及对获取的训练样本进行提纯，得到每个子区域的遥感影像解译训练样本集。

具体地，所述对每个子区域中的完全一致区域使用简单非迭代聚类算法进行面向对象的分割进一步筛选；根据筛选后完全一致区域每一地类的数量比例分层随机抽样，提取初始的遥感影像解译训练样本的步骤中，包括：在各子区域内，应用提取的完全一致区域，使用简单非迭代聚类算法进行面向对象的分割，依据1σ准则和分割区域的多波段合成值进一步筛选各地类可靠的训练区，其中设置简单非迭代聚类算法的参数compactness为0，connectivity为8，neighborhoodsize为64，size为3；

在筛选后的完全一致性区域内依据每一地类的数量比例进行分层随机抽样，以栅格为中心建立90m×90m的正方形样本，作为初始的遥感影像解译训练样本。

具体地，所述对获取的训练样本进行提纯，得到每个子区域的遥感影像解译训练样本集的步骤中，包括两个部分：

第一部分：将初步提取的样本导入谷歌地球引擎平台，获取计算的样区范围内的植被指数ndvi的均值和标准差，计算样区内ndvi数值的变异系数，即变异系数＝标准差/均值，剔除变异系数大于1的样本；第二部分：建立直方图统计筛选方法对训练样本进行提纯，将上述步骤的样本ndvi依据从低到高排序建立数据直方图，在直方图首尾两端分别去除2.5％的样本，保留中心95％的样本；为保证每个地类的训练样本数不低于1000个，本例中，灌丛、建设用地、冰川积雪和苔原等地类在部分区域样本数不足，从相邻的子区域中进行随机抽取补充。

步骤106，根据遥感影像多波段合成值、系列遥感指数和地形辅助数据，生成土地覆被解译特征集合，通过机器学习算法对所述每个子区域的遥感影像解译训练样本集进行构建分类器，对不一致区域进行影像解译；其中，所述系列遥感指数包括所述植被指数、建筑指数、裸土指数、燃烧指数、雪指数和水体指数。

在本实施例中，所述遥感影像多波段合成值包括所述遥感影像数据在一年内4个时期的各波段中值合成及全年所有影像的波段数值标准差，其中，4个时期分别为1-3月、4-6月、7-9月、10-12月。所述地形辅助数据包括海拔和坡度数据。所述机器学习算法为随机森林算法。

具体地，基于谷歌地球引擎平台，生成用于遥感影像土地覆被解译特征集合，包括遥感影像多波段合成值、系列遥感指数和地形等辅助数据；将一年内的所有landsat影像分为4个时段的影像，分别是1-3月、4-6月、7-9月以及10-12月等四个时段，生成各时段的中值合成影像，并生成landsat影像各波段和遥感指数在一年中的标准差合成影像，共生成74个子类的影像集合。

应用随机森林算法，将完全一致区域内提纯的样本作为算法的训练样本集合，设置参数treenumber为300，每棵树选取的特征数量为波段和遥感指数总数的平方根，每棵树输入的训练数量为总样本量的63％，进行不一致区域进行修正。

步骤107，将目标区域中不一致区域进行影像解译后的土地覆被结果与完全一致区域的土地覆被进行镶嵌，生成土地覆被修正数据。

具体地，基于最优面积准则将不一致区域的修正结果重新采样至完全一致区域的空间分辨率并进行空间镶嵌，最终的土地覆被修正数据。

本发明技术方案提供了一种可以快速而准确地进行大区域范围的土地覆被制图技术，通过建立多源数据空间一致性识别及样本统计筛选方法，提取已有土地覆被产品的有效信息，并基于遥感云计算平台和海量遥感影像，应用机器学习算法对多种土地覆被产品间的不一致信息进行修正，实现了大范围内多源土地覆被产品的融合和修正，具有简易快速以及分类精度高的特点。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步的详细说明，应当理解，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围。特别指出，对于本领域技术人员来说，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。