CN109740570B

CN109740570B - 一种基于地形指数和物候差异的梯田遥感识别方法

Info

Publication number: CN109740570B
Application number: CN201910056183.1A
Authority: CN
Inventors: 周祖昊; 韩宁; 严子奇; 姜欣彤; 刘佳嘉; 顾世祥; 陈刚; 梅伟; 浦承松
Original assignee: YUNNAN PROVINCE WATER RESOURCES AND HYDROPOWER SURVEY AND DESIGN INSTITUTE; China Institute of Water Resources and Hydropower Research
Current assignee: YUNNAN PROVINCE WATER RESOURCES AND HYDROPOWER SURVEY AND DESIGN INSTITUTE; China Institute of Water Resources and Hydropower Research
Priority date: 2019-01-18
Filing date: 2019-01-18
Publication date: 2020-01-03
Anticipated expiration: 2039-01-18
Also published as: CN109740570A

Abstract

本发明公开了一种基于地形指数和物候差异的梯田遥感识别方法，涉及水土保持遥感监测技术领域。该方法首先利用易于获取的中等分辨率的Landsat遥感影像提取不同物候期的NDVI差值，然后结合包括海拔、坡度和坡度变率在内的地形因子信息用以有效识别梯田的面积和位置，以此提取梯田分布信息，解决了大范围、多时相、高效率、低成本提取梯田的难题。

Description

一种基于地形指数和物候差异的梯田遥感识别方法

技术领域

本发明涉及水土保持遥感监测技术领域，尤其涉及一种基于地形指数和物候差异的梯田遥感识别方法。

背景技术

梯田是治理坡面水土流失的一项行之有效的水土保持工程措施，具有明显的保水、保土与保肥效用，在我国黄土高原地区的水土流失防治工作中发挥着不可替代的作用。梯田调查是水土保持动态监测和评价工作中不可或缺的环节，准确快速的梯田位置和面积等信息提取为估算区域尺度水土流失量及厘定梯田在水沙减少中的作用提供理论依据和技术支持。

目前，通常采用外业调查的统计法确定梯田面积，这种梯田调查方法，速度慢、工期长而且耗费大量人力和物力，不能满足现阶段较大空间尺度上高效获取梯田分布面积和位置的需求。此外，高精度遥感影像如5m以内精度的SPOT影像、IKONOS影像和Quick-Bird卫星数据等图像上地物的形状和纹理信息突出，可有效辨别地物类型，也多用于纹理和光谱特征明显的梯田提取。但由于高精度遥感影像价格昂贵，数据获取困难，难以实现大范围、多时相的梯田面积和位置的遥感解译，即使在小流域范围内进行梯田面积提取，也多采用面向对象法利用遥感影像基于梯田纹理和光谱特征信息对梯田进行提取，对遥感影像的精度依赖过高，或仅仅考虑梯田的地形因子特征尝试梯田提取，提取结果往往精度较差。

发明内容

本发明的目的在于提供一种基于地形指数和物候差异的梯田遥感识别方法，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明采用的技术方案如下：

一种基于地形指数和物候差异的梯田遥感识别方法，包括如下步骤：

S1，实地调查研究区梯田作物种植类型及其物候条件，根据调查结果选取对应类型作物收割前后的中等分辨率Landsat8遥感影像；

S2，对选取的遥感影像进行包括辐射定标、大气校正、全色和多光谱波段精度融合、研究区剪裁在内的预处理步骤，得到处理后的遥感影像；

S3，基于处理后的遥感影像计算月时间尺度、空间分辨率为15m的不同物候期的植被指数NDVI值，并通过数据滤波后求差值△NDVI，区域△NDVI值越大，说明该区域为耕地的概率越大；

S4，根据DEM数据，获取研究区的的坡度Slope、坡度变率PC和海拔信息；

S5，参考研究区高精度Google地图，分别选取提取梯田的△NDVI值、DEM值、Slope值和PC值的阈值范围；并根据该范围进行提取，分别得到最终的△NDVI值、DEM值、Slope值和PC值的提取图层；

S6，在ARCGIS中使用Intersect叠加分析工具提取△NDVI、DEM、Slope和PC图层的相交部分，并把相交结果确定为研究区最终的梯田分布。

优选地，S3中，所述植被指数NDVI值按照如下公式进行计算：

NDVI＝(NIR-R)/(NIR+R)

其中，NIR表示近红外波段，R表示可见光红波段；

对于landsat8数据：NDVI＝(B5-B4)/(B5+B4)，在ENVI5.3软件中通过NDVI工具或者波段计算。

优选地，S4中，所述研究区的的坡度Slope、坡度变率PC和海拔信息按照如下方法获取：在ARCGIS软件中利用空间分析工具包提取。

优选地，S5中，所述ΔNDVI值、DEM值、Slope值和PC值的阈值范围按照如下方法进行选取：首先参照梯田分布的调查统计数据和相关著作给出的范围作为粗略的经验范围，然后将ΔNDVI值、DEM值、Slope值和PC值图层叠加到相同时期的Google地图里，采用目视解译的方法不断的尝试缩小提取范围使之与Google地图里的梯田分布最为接近，从而确定为合适的阈值范围。

本发明的有益效果是：本发明提供的基于地形指数和物候差异的梯田遥感识别方法，首先利用易于获取的中等分辨率的Landsat遥感影像提取不同物候期的NDVI差值，然后结合包括海拔、坡度和坡度变率在内的地形因子信息用以有效识别梯田的面积和位置，以此提取梯田分布信息，解决了大范围、多时相、高效率、低成本提取梯田的难题。

附图说明

图1是本发明提供的基于地形指数和物候差异的梯田遥感识别方法流程示意图；

图2是具体实施例中，研究区8月份与10月份NDVI差值分布图；

图3是具体实施例中，研究区DEM分布图；

图4是具体实施例中，研究区坡度Slope分布图；

图5是具体实施例中，研究区坡度变率PC分布图；

图6是具体实施例中，阈值范围内△NDVI分布图；

图7是具体实施例中，阈值范围内DEM分布图；

图8是具体实施例中，阈值范围内坡度Slope分布图；

图9是具体实施例中，阈值范围内坡度变率PC分布图；

图10是具体实施例中，研究区Google earth图；

图11是具体实施例中，研究区梯田分布图；

图12是具体实施例中，梯田分布区域放大图。

具体实施方式

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

本发明中涉及到的术语解释如下：

NDVI(Normalized Difference Vegetation Index，归一化植被指数):反映植被长势和营养信息的重要参数之一，由近红外波段与红光波段的反射率值计算得到。

物候：指植物在一年的生长中，随着气候的季节性变化而发生萌芽、抽枝、展叶、开花、结果及落叶、休眠等规律性变化的现象，称之为物候或物候现象。

坡度(Slope)：指地表单元陡缓程度，通常把坡面的垂直高度和水平距离的比值叫做坡度。

坡度变率(profile curvature，PC)：坡度沿地面最大坡降方向的地面高程变化率，是依据坡度的计算原理，在所提取的坡度值的基础上对地面每一点再求算一次坡度，即坡度之坡度。坡度变率在一定程度上可以很好地反映坡面曲率信息。

如图1所示，本发明实施例提供了一种基于地形指数和物候差异的梯田遥感识别方法，包括如下步骤：

S3中，所述植被指数NDVI值可以按照如下公式进行计算：

NDVI＝(NIR-R)/(NIR+R)

其中，NIR表示近红外波段，R表示可见光红波段；

S4中，所述研究区的坡度Slope、坡度变率PC和海拔信息按照如下方法获取：在ARCGIS软件中利用空间分析工具包提取。

S5中，所述△NDVI值、DEM值、Slope值和PC值的阈值范围按照如下方法进行选取：首先参照梯田分布的调查统计数据和相关著作给出的范围作为粗略的经验范围，然后将△NDVI值、DEM值、Slope值和PC值图层叠加到相同时期的Google地图里，采用目视解译的方法不断的尝试缩小提取范围使之与Google地图里的梯田分布最为接近，从而确定为合适的阈值范围。

具体实施例：

选取甘肃省兰州市榆中县部分区域作为研究实例，按照本发明提供的上述方法对梯田进行识别，具体按照如下步骤进行实实施：

步骤一，首先实地考察研究区梯田作物种植类型及其物候条件，调查发现，该地区主要种植春小麦、冬小麦和夏玉米，其中夏玉米在8月中旬达到NDVI值的巅峰，而一般在10月上旬全部收割完毕，NDVI值达到最小值，而该时期非作物的林、草地分布区，其NDVI值一般变化不大。

根据上述结果，本实施例中，需要分别下载8月份中旬和同年10月份上旬的Landsat 8遥感影像，通过筛选选取了质量较高的2015年10月上旬的Landsat8遥感影像：LC81300352015283LGN00和2015年8月中旬的Landsat 8遥感影像：LC81300352015235LGN00。

步骤二，对选取的两个Landsat 8遥感影像数据，在ENVI软件中进行预处理，包括辐射定标、大气校正、图像裁剪和融合处理。借助ENVI遥感软件的Gram-Schmidt PanSharpening工具对两时相的图像分别实现多光谱和全色光谱的波段融合，得到15m精度的两期图像，然后用ENVI遥感软件的NDVI工具分别计算两时期的NDVI值，反演出了月时间尺度、空间分辨率为15m的NDVI，并利用数据滤波工具Median Filter对NDVI进行滤波处理，最后求差值△NDVI，进而得到研究区作物收割前后NDVI差值分布，区域△NDVI值越大，说明该区域受人类活动影响较大，可能正处于作物收割前后，其为耕地的概率越大，一般将大于0的区域初步判断为耕地。结果如图2所示。从图2可以看出研究区8月份与10月份△NDVI分布在大于0.03的范围内呈明显的斑块连接状，初步判断为主要耕地区。

步骤三，借助于10m精度的DEM数据，在提取研究区海拔的基础上，计算相应坡度(Slope)和坡度变率(profile curvature，PC)地形指数，提取结果分别如图3所示，梯田依山而建，主要分布在山坡上，从图3可以看出，梯田所处地势相对较高。

在ARCGIS中打开工具箱，使用Spatial Analyst Tools-Surface-Slope工具提取DEM的坡度，得到Slope分布图，如图4所示，从图4可以看出在坡度在低于4°的区域地连续大片且分布清晰分明，判断为非坡耕地和非梯田的坝地或者其它分布。

在ARCGIS中打开工具箱，使用Spatial Analyst Tools-Surface-Slope工具提取Slope的坡度，得到坡度变率PC值，如图5所示。梯田是在丘陵山坡地上沿等高线方向修筑的条状阶台式或波浪式断面的田地，其坡度变率一般较大，而坡耕地的坡度均匀变化，因此坡度变率较小。从图5可以看出坡度变率分布规律较杂乱，可尝试通过去除较小的坡度变率去除坡耕地的影响。

步骤四，参考研究区高精度Google地图，在ARCGIS中经过多次试验，分别针对不同图层利用栅格计算器的工具(Spatial analyst——Raster calculator)选取合适的△NDVI值、DEM值、Slope值和PC值的阈值范围，如表1所示。

表1

步骤五，根据步骤四中确定的阈值范围，对梯田进行提取，得到阈值范围内的△NDVI分布图(如图6所示)、DEM分布图(如图7所示)、坡度Slope分布图(如图8所示)、坡度变率PC分布图(如图9所示)，从图6-9中可以看出按照表1设置的提取范围处理后，各个提取的图层结果显示非白色的区域是有梯田分布的区域，而图6-9的存在的公共部分可以认定为梯田分布。

在ARCGIS中打开工具箱，使用Analysis Tools——Overly——Intersect工具提取△NDVI、DEM、Slope和PC图层的相交部分，得到梯田解译图层，最后，将上述提取结果通过参考同期高精度Google影像(如图10所示)和相应年份土地利用分类图进行目视解译和手动修正，得到最终的梯田分布结果，可如图11和图12所示，从图中可以看出根据上述步骤实现的梯田提取结果(由线条勾勒的区域)与同时期的高精度Google影像上目视可见的梯田分布区域一致，从而可确定本方法的可靠性与准确性。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：本发明提供的基于地形指数和物候差异的梯田遥感识别方法，首先利用易于获取的中等分辨率的Landsat遥感影像提取不同物候期的NDVI差值，然后结合包括海拔、坡度和坡度变率在内的地形因子信息用以有效识别梯田的面积和位置，以此提取梯田分布信息，解决了大范围、多时相、高效率、低成本提取梯田的难题。

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

Claims

1.一种基于地形指数和物候差异的梯田遥感识别方法，其特征在于，包括如下步骤：

S6，在ARCGIS中使用Intersect叠加分析工具提取△NDVI、DEM、Slope和PC图层的相交部分，并把相交结果确定为研究区最终的梯田分布；

2.根据权利要求1所述的基于地形指数和物候差异的梯田遥感识别方法，其特征在于，S3中，所述植被指数NDVI值按照如下公式进行计算：

NDVI＝(NIR-R)/(NIR+R)

其中，NIR表示近红外波段，R表示可见光红波段；

在ENVI5.3软件中通过NDVI工具或者波段计算。

3.根据权利要求1所述的基于地形指数和物候差异的梯田遥感识别方法，其特征在于，S4中，所述研究区的的坡度Slope、坡度变率PC和海拔信息按照如下方法获取：在ARCGIS软件中利用空间分析工具包提取。