CN109856056B - 一种沙化土地遥感快速识别方法 - Google Patents

Abstract

本发明提供一种沙化土地遥感快速识别方法，属于遥感监测技术领域；所述方法包括，获取数据，并对数据进行预处理；在所述进行预处理的数据中提取非植被类型信息；对提取的所述非植被类型信息进行处理，包括，提取水体、建筑用地类型信息；提取裸地类型信息；分析、识别沙化土地。本发明的数据源时间序列长、时间分辨率高；监测周期短，以一年为周期，可连续每年掌握沙化土地的发展或逆转。

Description

一种沙化土地遥感快速识别方法

技术领域

本发明属于遥感监测技术领域，特别涉及一种沙化土地遥感快速识别方法。

背景技术

土地沙化是指因气候变化和人类活动所导致的天然沙漠扩张和沙质土壤上植被破坏、沙土裸露的过程。土地沙化，是指主要因人类不合理活动所导致的天然沙漠扩张和沙质土壤上植被及覆盖物被破坏，形成流沙及沙土裸露的过程；本法所称的沙化土地包括已经沙化的土地和具有明显沙化趋势的土地。是否会发生沙化，决定的因素在于土壤中含有多少水分可供植物吸收、利用，并通过植物叶面而蒸发。任何破坏土壤水分的因素都会最终导致土壤沙化，土地沙化的大面积蔓延就是荒漠化，是最严重的全球环境问题之一。截至2014年，全国沙化土地总面积17211.75万hm²，占国土总面积的17.93％，分布于除上海市、台湾省和香港、澳门特别行政区外的30个省的920个县。

遥感具有大面积同步观测、时效性高以及成本低廉等特点，可以相对便捷地获得大范围、全覆盖的沙化土地空间分布信息，已被广泛应用于沙化土地的提取。目前我国沙化监测是每5年一次，对于沙化来说，因其对气候变化比较敏感，年际间甚或更短时期的气候波动都可能导致土地沙化情况加剧或减弱，如何以一年为周期或更短时期快速识别沙化土地的发展或逆转，是现阶段沙化监测需要探索的一个问题。

目前对于沙化土地变化监测的研究，大多数是对多时相影像分别进行分类的基础上，进行不同时期沙化土地的对比分析。虽然此技术可实现沙化土地监测且方法成熟，但是由于数据量大导致此方法计算过程繁琐、处理速度慢，且由于监测时段不连续、过程孤立导致一年时间周期或更短周期内土地沙化情况无法及时掌握。

发明内容

针对上述问题，本发明提供一种沙化土地遥感快速识别方法；所述方法包括：

获取数据，并对数据进行预处理；

在所述进行预处理的数据中提取非植被类型信息；

对提取的所述非植被类型信息进行处理，包括，

提取水体、建筑用地类型信息；

提取裸地类型信息；

分析、识别沙化土地。

进一步地，所述获取数据是获取连续时间序列的AVH13 NDVI影像数据。

进一步地，所述数据进行预处理是利用ArcGIS Python对获取的数据影像，进行批量拼接、投影、裁剪的预处理工作。

进一步地，所述进行预处理的数据中提取非植被类型信息为对预处理后的数据影像进行像元统计工作，统计一年内NDVI值保持在一定范围内不变的像元，从而提取非植被类型空间分布信息。

进一步地，所述提取水体、建筑用地类型信息具体操作为，

NDVI月值数据合成；

多波段合成；

主成分分析；

非监督分类。

进一步地，所述NDVI月值数据采用国际通用的最大化合成法获得。

进一步地，所述多波段合成是将全年12个月NDVI最大值合成影像进行叠加，把12个月影像进行多通道合并，使其成为一幅具有12个层的复合影像。

进一步地，所述非监督分类的方法包括ISODATA迭代聚类算法、决策树分类法、NDVI阈值分类法中的一种。

进一步地，所述裸地类型信息是从非植被类型中排除水体、建筑用地类型，从筛选中得出裸地类型信息。

进一步地，所述分析、识别沙化土地是在时间序列上进行土地利用变化分析，得到一定时间范围内裸地面积增加的区域。

本发明的数据源时间序列长、时间分辨率高；监测周期短，以一年为周期，可连续每年掌握沙化土地的发展或逆转。本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所指出的结构来实现和获得。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1示出了根据本发明实施例的整体流程示意图；

图2示出了根据本发明实施例的操作流程示意图。

具体实施方式

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

图1示出了根据本发明实施例的整体流程示意图。如图1所示，一种沙化土地遥感快速识别方法，所述方法包括，获取数据，并对数据进行预处理；在所述进行预处理的数据中提取非植被类型信息；对提取的所述非植被类型信息进行处理，包括，提取水体、建筑用地类型信息；提取裸地类型信息；分析、识别沙化土地。

以一年内，某地区的土地沙化识别为例，进行操作：

步骤一，获取数据，并对数据进行预处理；

具体的，所述获取数据是获取连续时间序列的AVH13 NDVI(AVH13，AVHRR地表反射率产品；AVHRR，Advanced Very High Resolution Radiometer，就是NOAA系列气象卫星传感器；NDVI，Normalized Difference Vegetation Index，归一化植被指数)影像数据，其获得影像数据空间分辨率为8KM，时间分辨率为日。

具体的，所述数据进行预处理包括对研究区范围内AVH13 NDVI影像数据进行拼接、裁剪、投影。示例性的，投影坐标使用UTM ZONE 48 North(UTM，Unified ThreatManagement，通用横轴墨卡托投影北纬48度带)，WGS84(WGS84，World GeodeticSystem1984，1984年世界大地测量系统)。

另外需要说明的是，由于AVH13 NDVI的影像数据对瑞利散射和臭氧吸收作了校正，又通过最大值合成法(MAC)对其他噪音作了校正，可以直接使用；其中瑞利散射是光学散射的一种；光的散射以及臭氧对获取的数据有影响，因而采用AVH13 NDVI可消除这些影响。

具体的，影像批量拼接、投影、裁剪等预处理工作可利用ArcGIS Python(ArcGIS软件脚本语言)实现。

步骤二，在所述进行预处理的数据中提取非植被类型信息；

具体的，利用步骤一处理后的数据进行像元统计工作，统计一年内NDVI(Normalized Difference Vegetation Index，归一化植被指数)值几乎保持不变的像元，即统计一年内NDVI值保持在一个范围内不变的像元，从而提取非植被类型空间分布信息。

其中，AVH13 NDVI数据时间分辨率为日，可利用ArcGIS Python调用像元统计工具，对一年内日NDVI数据进行批处理；需要说明的是，将NDVI值小于一个标准差的像元均视为一年内NDVI值几乎保持不变的像元，准确的描述就是，确定一个NDVI的标准差值，只要检测得到的一年内的NDVI变化值小于这个标准差值，即可视为NDVI值几乎保持不变的像元；NDVI值一年内几乎保持不变即可认为其非植被类型，因为植被类型一年内随着季节的改变，其形态如植物到了秋天叶子枯萎、落叶，春天长新芽都会发生改变，因此其影像数据也会发生改变，一年内的NDVI变化值会大于这个标准差值，视为植被类型；而水体、建筑用地、裸地其一年内形态的变化小于这个标准差值，则可视为几乎不变，因此将在一个小于标准差内NDVI值几乎保持不变的提出，即可实现提取非植被类型。

提取非植被类型信息具体的步骤为，

(1)分类体系的构建

分类体系是将土地划分为不同的类型，主要可以划分为自然地物、人工建造的地物或者划分为植被类型、非植被类型；此处根据需求将地物类型划分为植被类型、非植被类型进行讨论。

示例性的，设定甘南藏族自治州的地物类型为植被、非植被类型，其中非植被类型划分为裸地、水体以及建筑用地。

(2)非植被类型信息提取

利用ArcGIS Python调用像元统计工具，统计NDVI值在一年中的变化小于一个标准差的像元，符合上述条件的像元，可以看作是一年内NDVI值几乎保持不变的像元，而这样的像元在很大意义上可以归为非植被类型的像元。

示例性的，进行连续时间序列上像元NDVI值批量统计的Python(计算机程序设计语言)核心代码如下：

在经过叠加求和后的影像图中，像元值为0的像元为符合条件的非植被类型像元；利用ArcGIS10.2(全面的GIS平台)软件对上述结果进行重分类，像元值为0的像元赋值为0，其他值赋值为1，从而得到包含植被、非植被两类型的空间分布。

步骤三，所述对提取的所述非植被类型信息进行处理，包括，

第一：提取水体、建筑用地类型信息；其提取步骤为，

(1)NDVI月值数据合成

NDVI月值数据采用国际通用的最大化合成法(MVC)获得，具体的，此法假设每月中NDVI值最大的那一天的天气晴朗，不受云层的影响，进一步消除大气、云和太阳高度角的部分干扰，以使NDVI更为贴近实际地反映每月地表植被覆盖状况；通过这个处理过程，得到连续几年中(如2013-2015年)每年每月的NDVI最大值，作为该月的植被指标。

示例性的，通过ArcGIS Python调用像元统计工具，批量统计叠加栅格数据，最终合成结果为所有栅格的最大值，其核心代码如下：

currentDir＝"F:\\GISdata\\gsdata\\max\\data\\2015max\\aug15\\"

IMGFileList＝glob.glob(currentDir+"/*.tif")

outCellStats＝CellStatistics(IMGFileList，"MAXIMUM"，"DATA")

outCellStats.save("F:\\GISdata\\gsdata\\max\\process\\2015m\\201508\\ma x1508.tif")

print('OK')

(2)多波段合成

将全年12个月NDVI最大值合成影像进行叠加，把12个月影像进行多通道合并，使其成为一幅具有12个层的复合影像；示例性的，此操作可以通过ArcGIS10.2软件中数据管理工具/栅格/栅格处理/波段合成工具实现。

(3)主成分分析

为了提高运算速度，减少图像数据信息的含量，对全年12个月NDVI复合影像进行主成分分析。

具体的，主成分分析(Principal component analysis，PCA)是减低数据维数的一种很有效的方法，属于光谱增强处理方法的一种，可以去除波段间的多余信息、将多波段图像信息压缩到比原波段更有效的少数几个转换波段；这种方法常被用于多时相NDVI数据处理，对12个月的NDVI合成影像，在用于分类之前作主成分分析，计算出各波段的统计特征和主成分分析的各主成分及其信息载荷。

示例性的，利用ENVI5.1(遥感图像处理平台)中forward PCA Rotation NewStatistics and Rotate工具对全年12个月NDVI合成影像进行主成分分析，分析结果第一、第二、第三分量信息载荷达到90％以上，能基本反应图像信息的含义；可用前三分量作为新的影像的三个波段，对研究区的土地类型进行分类。

(4)非监督分类

具体的，非监督分类采用ISODATA(迭代自组织数据分析算法)，其分类速度快、精度较好，现阶段常用的分类方法如决策树分类、NDVI阈值分类等方法亦可达到分类要求。

示例性的，对全年12个月NDVI综合数据经过主成分分析(PCA处理后，取前3个主要分量作为新的影像的三个波段，对其实施非监督分类，采用ISODATA，迭代次数设为20次，结合实验区的实际土地覆被情况，将实验区分为20类，最终合并为6类，从中提取水体及建筑用地类型。

第二：提取裸地类型信息；其提取步骤为，

从非植被类型中排除水体、建筑用地类型，从筛选中得出裸地类型信息；示例性的，首先，利用ArcGIS10.2按属性提取工具，可提取时序NDVI分类图中的水体、建设用地类型；之后，利用ArcGIS10.2中栅格数据裁剪工具，筛选出非植被类型中的裸地类型信息。

步骤四，分析、识别沙化土地。

由于土地沙化是一个渐进的过程，在时间序列上进行土地利用变化分析，可得到一定时间范围内裸地面积增加的区域；其中一段时间范围即为研究时间段内的时间范围。

示例性的，以提取某地实验区2013年至2014年间沙化土地为例，在ArcGIS10.2软件中，首先对2013年裸地空间分布图进行重分类，将裸地属性信息赋值为2，其他地类赋值为1；同时，对2014年裸地空间分布图进行重分类，裸地类型信息赋值为0，其他地类赋值为3；随后，利用栅格计算器对两期影像进行减法运算，裸地类型信息不变的像元值为2，其他地类改变为裸地的像元值为1，即为出现土地沙化的区域。

示例性的，利用全国土地利用数据产品(空间分辨率1KM×1KM)对所述方法的结果进行精度评价，评价公式如下：

P＝S_沙S_总×100％

其中，P为模型识别正确率，S_沙是正确识别为沙化土地的面积，S_总是识别的沙化土地总面积。

最终验证结果，此方法识别沙化土地精度达到88.2％，能够较好识别出试验区域沙化土地空间分布。

综上所述，在进行沙化土地识别时，例如进行2013-2014年某地区的沙化土地的识别，首先需要对获取数据，获取数据时主要是通过获取这一年连续时间序列的AVH13 NDVI影像数据，通过数据主要记录了需要进行沙化土地识别地区的样貌变化；数据获取后通过ArcGIS Python对AVH13 NDVI影像数据进行拼接、裁剪、投影，处理后的数据进行像元统计工作，统计一年内NDVI值保持在一定范围内不变的像元，一年内保持不变的像元即可视为非植被类型。

非植被类型确定后，在非植被类型中删除水体、建筑类型用地，剩余的即可视为裸地类型的空间，裸地类型的空间内2013年-2014年裸地变化的区域，即增加的区域即可视为沙化的区域。具体的在非植被类型中删除水体、建筑类型用地，主要是通过首先利用最大合成法确定NDVI月值数据，而后利用多段合成法，得到NDVI合成数据，最后通过主题分析以及非监督分类的方法从中提取水体及建筑用地类型；在水体及建筑用地类型中通过像元的筛选得到裸地空间分布，最后通过时序的分析识别本年内沙化土地的增加量。

本发明所使用的AVH13 NDVI影像数据，图像分辨率为8km×8km，时间分辨率以日为单位，具有时间序列长、覆盖范围广、较强地植被动态变化表征能力等特点；高时间分辨率对于沙化土地的监测具有重要意义，干旱区生长大量短命植物，植物利用早春雨水和融化的雪水提供的土壤湿度和一定的温度生长发育，并在夏季干热季节来临之前短短几周到2个月时间里迅速完成生活周期，随后整个植株或地上部分干枯死亡；利用高时间分辨率影像监测土地沙化相对于周期较长的MODIS NDVI(MODIS植被指数产品，其中MODIS为搭载在terra和aqua卫星上的一个重要的传感器)、Landsa tETM+(Landsat卫星为美国国家航空航天局陆地卫星，ETM+为传感器名称)等具有缩短监测周期、减少误差等优势。

尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (7)

1.一种沙化土地遥感快速识别方法，其特征在于，所述方法包括，

获取数据，并对数据进行预处理；其中，所述获取数据是获取连续时间序列的AVH13NDVI影像数据，其获得影像数据空间分辨率为8KM，时间分辨率为日；所述数据进行预处理是利用ArcGIS Python对获取的数据影像，进行批量拼接、投影、裁剪的预处理工作；

在所述进行预处理的数据中提取非植被类型信息；其中，所述进行预处理的数据中提取非植被类型信息为对预处理后的数据影像进行像元统计工作，统计一年内NDVI值保持在一定范围内不变的像元，从而提取非植被类型空间分布信息；

对提取的所述非植被类型信息进行处理，包括，

提取水体、建筑用地类型信息；

提取裸地类型信息；

分析、识别沙化土地。

2.根据权利要求1所述的快速识别方法，其特征在于，所述提取水体、建筑用地类型信息具体操作为，

NDVI月值数据合成；

多波段合成；

主成分分析；

非监督分类。

3.根据权利要求2所述的快速识别方法，其特征在于，所述NDVI月值数据采用国际通用的最大化合成法获得。

4.根据权利要求2所述的快速识别方法，其特征在于，所述多波段合成是将全年12个月NDVI最大值合成影像进行叠加，把12个月影像进行多通道合并，使其成为一幅具有12个层的复合影像。

5.根据权利要求2所述的快速识别方法，其特征在于，所述非监督分类的方法包括ISODATA迭代聚类算法、决策树分类法、NDVI阈值分类法中的一种。

6.根据权利要求1-5中任意一项所述的快速识别方法，其特征在于，所述裸地类型信息是从非植被类型中排除水体、建筑用地类型，从筛选中得出裸地类型信息。

7.根据权利要求1-5中任意一项所述的快速识别方法，其特征在于，所述分析、识别沙化土地是在时间序列上进行土地利用变化分析，得到一定时间范围内裸地面积增加的区域。

Images