CN109635713A

CN109635713A - 高原山地阴影区冰川识别方法

Info

Publication number: CN109635713A
Application number: CN201811492914.9A
Authority: CN
Inventors: 季漩; 陈云芳; 罗贤; 李运刚
Original assignee: Yunnan University YNU; Yunnan Normal University
Current assignee: Yunnan University YNU; Yunnan Normal University
Priority date: 2018-12-07
Filing date: 2018-12-07
Publication date: 2019-04-16

Abstract

本发明公开了一种高原山地阴影区冰川识别方法，获取Landsat 8OLI高原山地影像，然后对其进行辐射定标预处理，并利用FLAASH模型对定标后的影像进行大气校正；基于DEM数据和遥感成像时的太阳方位角、高度角信息，利用地理信息系统软件模拟出相应的阴影区域，进行阴影区识别；依据阴影区冰川与非冰川的波段反射特征，确认蓝光波段为阴影区冰川信息识别优势波段；利用阴影区冰川信息识别优势波段，通过增强型识别方法计算得到ENDSI灰度图，再依据计算结果直方图选择合适的阈值，对ENDSI灰度图进行分割，得到阴影区冰川信息。采用优势波段的增强雪被指数法有助于提高整体高原山地阴影区冰川识别工作效率。

Description

高原山地阴影区冰川识别方法

技术领域

本发明属于图像识别技术领域，特别是涉及一种基于Landsat 8 OLI影像的高原山地阴影区冰川识别方法。

背景技术

冰川变化是气候变化的产物，同时也是全球气候变化的敏感指示器，尤其是处于中纬度地区的高原山地冰川对气候变化响应极为敏感，其变化对区域气候、生态与环境、水文水资源等均有重要影响。因此，山地冰川变化一直是地学研究的热点。而无论是冰川对气候变化响应还是对区域水资源的影响，对山地冰川的提取、变化监测都是面临的首要问题。

山地冰川通常地处偏远，地形复杂，实地监测成本巨大，很难大范围开展变化监测工作。随着遥感技术的发展，冰川变化研究可以获得及时、准确的数据源，一定程度上解决了高山区冰川信息资料稀缺等问题。

基于遥感影像的冰川信息提取方法已有大量研究，总体上可以分为目视解译和计算机辅助分类。目视解译方法虽然精度较高，但目视解译工作量大、周期较长、耗时耗力。计算机辅助分类方法精度整体上不如人工目视解译，但能快速获取大尺度冰川区域的信息，优越性十分明显。

常用的计算机辅助分类方法有比值阈值法、雪盖指数法等。对于早期冰川编目工作而言，大多选用Landsat数据，基于冰雪在可见光波段高反射率、近红外波段低反射率的原理，进行波段运算处理进而获取冰川信息。随着遥感技术的成熟，面向对象的信息提取方法、雷达INSAR干涉技术等也陆续应用于冰川的提取与监测中。虽然基于遥感的冰川信息提取方法有很多，但由于冰川类型、影像空间分辨率、影像质量及研究区域的不同，每种方法也都有自己的适用性，很难找到一种方法适用于所有地区冰川的提取。对于高原山地区域而言，遥感影像通常会有较大面积的山体阴影。阴影使地物目标反映的信息量有所损失或受到干扰，在遥感影像数据上表现为DN值偏低，难以判读。因此，基于遥感影像的山体阴影区冰川识别也成为一个技术难点。许多研究结果显示，监督和非监督分类方法不能有效识别阴影下的冰川，NDSI指数(归一化差异雪被指数)效果也较差；相对于NDSI，波段比值法对阴影区冰川的识别表现尚可，但是对于选用哪些波段效果更好，结论不一。

发明内容

本发明的目的在于提供一种高原山地阴影区冰川识别方法，基于landsat 8OLI数据，从分析阴影区冰川与非冰川光谱特征入手，提出阴影区冰川信息增强算法，对冰川信息进行有效提取，以解决山体阴影区冰川识别问题。

本发明所采用的技术方案是，高原山地阴影区冰川识别方法，具体步骤如下：

步骤S1、获取高原山地影像，并对其进行预处理；

步骤S2、对预处理后的影像，基于DEM数据和遥感成像时的太阳方位角、高度角信息，利用地理信息系统软件模拟出相应的阴影区域，进行阴影区识别；

步骤S3、依据阴影区冰川与非冰川的波段反射特征，确认阴影区冰川信息识别优势波段；

步骤S4、利用阴影区冰川信息识别优势波段，通过增强型识别方法提取阴影区冰川信息。

进一步的，所述高原山地影像是Landsat 8 OLI影像。

进一步的，所述步骤S1中对高原山地影像进行预处理，具体为先对其进行辐射定标预处理，再利用FLAASH模型对定标后的影像进行大气校正。

进一步的，所述步骤S3确定的优势波段为波长小于0.515μm的蓝光波段。

进一步的，所述步骤S4具体实现过程如下：

步骤S41、利用阴影区冰川信息识别优势波段，通过增强型雪被指数法，计算增强型雪被指数：

增强型雪被指数法计算公式为：

其中，BandX为阴影区冰川信息识别第X优势波段的DN值或反射率；BandY为阴影区冰川信息识别第Y优势波段的DN值或反射率，ENDSI为增强型雪被指数，形成ENDSI灰度图；

步骤S42、对ENDSI灰度图的像元值进行频率统计，得到计算结果直方图；

步骤S43、依据计算结果直方图选择合适的阈值，并利用阈值对ENDSI灰度图进行分割，即可得到阴影区冰川信息。

进一步的，所述步骤S43中依据计算结果直方图选择合适的阈值，是选择能够明显区分阴影冰川区与非冰川区的值，即选择计算结果直方图中两个明显的峰值之间的最低值，作为合适的阈值。

进一步的，所述步骤S42中依据阈值对ENDSI灰度图进行分割依据为：当ENDSI灰度图的某像素点的灰度值大于或等于阈值时，该点为阴影冰川区；当ENDSI灰度图的某像素点的灰度值小于阈值时，该点为阴影非冰川区。

本发明的有益效果是：本发明基于Landsat 8 OLI影像数据，分析了山体阴影区冰川与非冰川的波段反射特征，得出波长较短(小于0.515μm)的蓝光波段因具有更高的散射强度，是阴影区冰川识别的优势波段；在此基础上提出针对山体阴影区冰川信息提取的增强指数算法，并与常规的冰川信息提取方法进行效果对比，结果表明增强指数方法得到的直方图分割阈值更为明显；且从冰川信息提取结果来看，无论是空间分布还是面积误差比例，采用优势波段的增强雪被指数法效果最好。在高原山地区域进行大规模冰川提取时，采用本文提出的山体阴影区冰川信息增强指数算法，有助于提高整体高原山地阴影区冰川识别工作效率，且能够实时、准确的获得冰川变化数据。

附图说明

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

图1是2015年10月研究区位置及其中的冰川与阴影Landsat 8 OLI影像652波段假彩色合成示意图；

图2是基于DEM的阴影区提取结果及其中目视解译的冰川范围示意图；

图3是波段1～5在阴影区冰川部分的DN值直方图；

图4是波段1～5在阴影区非冰川部分的DN值直方图；

图5是阴影区六种提取方法的计算结果直方图；

图6是六种提取方法对阴影区冰川范围识别效果对比图。

具体实施方式

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

选择青藏高原中部唐古拉山的岗加曲巴冰川及其周边冰川群区域为实验研究区，具体位置如图1所示。实验区共涉及大小冰川37条，总面积约305km²。区域最高海拔6571m，山体高大，地形较为复杂，有较多冰川长时间位于山体阴影区。

本研究主要涉及的基础数据包括DEM数据和两期Landsat 8 OLI数据。

DEM数据选用30m空间分辨率的ASTER GDEM(Version 2)。数据覆盖范围包括了全球南北纬83°之间的区域。其Version2版本于2011年10月发布，相对于Version 1精度有了显著提升，能保证本研究对于山体阴影的提取准确度。

Landsat 8 OLI影像数据来自美国地质调查局网站(http://www.glovis.usgs.gov)。相较于之前的Landsat ETM+数据，OLI(Operational LandImager，陆地成像仪)传感器增加了两个波段：海蓝波段(band 1：0.433–0.453μm)主要应用海岸带观测；短波红外波段，又称卷云波段(band 9：1.360–1.390μm)包含水汽强吸收特征，主要用于云检测。该产品经过了系统辐射校正和地面控制点几何校正，并且经过数字高程模型(DEM)校正处理。本实验研究选用了两期OLI影像，具体信息见表1。其中2015年10月的影像质量很好，由于太阳高度角较大，山体阴影区很小。2016年12月的影像，太阳高度角较小，产生了较大范围的山体阴影区域。本研究对这两期影像进行了辐射定标等预处理，用FLAASH模型对定标后的影像进行了大气校正。

表1选用的LandSat 8 OLI影像数据信息

成像日期	条带号	云量	太阳方位角	太阳高度角
					2015年10月2日	138/37	0.13％	150.30882071	49.25647732
2016年12月23日	138/37	0.99％	157.44656004	29.88404772

阴影区的识别

本研究针对阴影区进行冰川识别，首先要提取准确的阴影区域。近年来，许多光谱信息遥感影像阴影检测方法被提出。但是这些方法还是会因“同物异谱”或“异物同谱”存在一定程度错分。

本发明基于DEM数据和遥感成像时的太阳方位角、高度角信息，利用地理信息系统软件模拟出相应的阴影区域，如图2阴影区域所示。这种方式相对于基于光谱信息的阴影检测方法，更方便快捷，在大规模冰川信息处理工作中，具有更高效率。

采用增强型识别方法提取阴影区冰川信息

(1)确定阴影区冰川信息识别优势波段：

根据图2目视解译的冰川边界，实验区阴影分为冰川区和非冰川区。分别按照阴影冰川区和阴影非冰川区统计各波段直方图，如图3和图4所示。由于遥感影像阴影区域的反射率随波长增加急剧减小，第6、7波段在阴影区DN极小，提供的信息十分有限，因此本文仅针对前5波段进行分析。

根据图3和图4统计直方图可知，不论是阴影区冰川部分，还是非冰川部分，随着波长增加，其DN值减小。

对于阴影区非冰川部分，1-5波段分别有97.1％、98.5％、99.1％、99.3％及99.8％的像元DN值在5以下。与此形成明显对比的是，阴影区冰川部分，1-5波段DN值小于5的像元比例分别有15.7％、37.2％、86.1％、91.1％及96.5％。

波段3、4、5无论在阴影区冰川部分还是非冰川部分，像元值重合比例过高，不利于区分出阴影区的冰川信息。而第1、2两个蓝光波段对于阴影区的冰川具有较高的辨识度，是优势波段。尽管如此，若采用波段阈值法对第1、2波段进行阈值分割提取，还是很难准确分离出阴影区冰川信息，因为仍有部分像元(15.7％、37.2％)与其他地物难以区分。因此需要对优势波段进一步增强处理。

(3)阴影区冰川信息增强型识别方法

对确定的阴影区，对冰川信息识别优势波段进行增强，形成灰度图：

波段比值法计算公式：

式中，Band_C表示遥感影像第X波段的反射率或DN值，Band_D表示遥感影像第Y波段的反射率或DN值，计算结果为Ratio灰度值图像；

归一化雪被指数(NDSI)计算公式：

式中，计算结果为NDSI灰度值图像；

增强型比值法计算公式：

增强型雪被指数法计算公式：

式中，Band_X为阴影区冰川信息识别第X优势波段的DN值或反射率；Band_Y为阴影区冰川信息识别第Y优势波段的DN值或反射率，计算结果ERatio和ENDSI分别为阴影区增强型波段比值灰度图和增强型雪被指数灰度图。

对ERatio灰度图或ENDSI灰度图选择合适的阈值，并对其进行分割即可得到冰川面积。

图2中底图为2016年12月23日的Landsat 8 OLI数据652波段假彩色合成影像。可以看出基于DEM提取出的山体阴影区与影像阴影区匹配良好，对于容易误分的水体和阴影，也能很好的区分开来。图2中阴影区冰川范围是通过对图1目视解译得到。根据该结果统计可知，实验区山体阴影总面积为34.42km²，其中冰川面积29.34km²，非冰川面积5.08km²。

由于两期影像属于相邻年份不同月份的数据，因此对其冰川分布范围的一致性进行了检验。对2016年12月的影像中非阴影区冰川及薄阴影区冰川进行目视解译，将解译结果与2015年10月质量良好的影像在相同范围内进行了对比，冰川面积一致性达到99.1％。在此基础上，假定该区域冰川在两期影像成像时间内并无显著变化。认为2015年10月无阴影的OLI影像目视解译冰川范围可作为计算机辅助解译结果的对比标准。

对Landsat 8 OLI影像第1、2波段，根据公式(1)～(4)分别对整个阴影区进行计算；此外，为了对比传统做法，对第1波段(可见光波段)和第5波段(近红外波段)也分别计算波段比值和雪被指数。根据计算结果，对于整个阴影区做直方图分析，如图5所示。

可以看出，对于优势波段B1和B2，使用增强指数计算的结果，其直方图“双峰”效果明显，分割阈值非常易于识别(ENDSI＝0.64，ERatio＝4.57)；而对于优势波段B1和B2使用常规的比值法和雪被指数法，其直方图总体呈“单峰”形态，对于阴影区的两种冰川与非冰川情况并无明显分离情况，阈值可辨识度较低。对于常规做法使用第1、5波段计算的比值和雪被指数，其直方图都呈现出一个“主峰”和一个不明显的“次峰”，阈值的可辨识度仍然不够明确。总体而言，对优势波段采用增强指标计算后，阴影区图像可分割性较强，阈值可辨识度相较于其他常规方法要高。

从图5可知优势波段增强算法的阈值非常明显，而常规比值法和雪被指数法从直方图上难以明确分割阈值，因此在使用这种方式提取冰川时，分别采用多次试验，找到最佳的分割阈值，优势波段B1和B2的NDSI阈值为0.11；优势波段B1和B2的Ratio阈值为1.25；B1和B5波段的NDSI阈值为0.82；B1和B5波段的Ratio阈值为10.2。根据最佳阈值对灰度图进行分割，当像素点灰度值大于最佳阈值时，判断该像素点为阴影区冰川范围，否则为阴影区非冰川范围，依次对整个灰度图像进行分割，提取出阴影区冰川范围，结果如图6所示。

图6中对比了不同方法提取的阴影区冰川范围，并与目视解译冰川范围作了比较。可以看出，采用优势波段增强指数提取的冰川效果较好；优势波段常规算法效果次之；而采用第1、5波段做比值计算和雪被指数计算得到的结果最差。使用优势波段增强算法提取的阴影区纯冰川范围效果极好，对于阴影区内冰川与非冰川混杂的情况，略有遗漏。第1、5波段做常规比值法和雪被指数法结果对于阴影区纯冰川情况的提取也有较多遗漏，而对于非冰川区又有较多区域被识别为冰川，说明采用第1、5波段做常规比值法和雪被指数对于阴影区冰川识别能力不高。

对于几种冰川信息提取方法的效果，采用面积误差比例进行评价，计算公式为：

式中A_τ为实际提取的阴影区冰川面积；A_M为通过目视解译法得到的阴影区冰川面积。计算结果P为误差率。

表2中列出了各种算法提出冰川的面积以及相对于目视解译面积的误差。可以看出基于优势波段的增强雪被指数提取结果误差最小。其次是，优势波段增强比值法和常规比值法，二者效果相近。而基于第1、5波段提取的面积误差较大。

表2不同方法提取的阴影区冰川面积比较

解译方式	非冰川面积(km<sup>2</sup>)	冰川面积(km<sup>2</sup>)	冰川面积误差率(％)
				ENDSI(B1&B2)	5.63	28.80	-1.83
ERatio(B1&B2)	5.72	28.70	-2.16
				NDSI(B1&B2)	5.92	28.50	-2.84
Ratio(B1&B2)	5.71	28.72	-2.11
				NDSI(B1&B5)	3.88	30.55	4.13
Ratio(B1&B5)	3.63	30.79	4.97
				目视解译	5.08	29.34	0.00

综合图6与表2的面积误差比例，基于优势波段增强的两种方法具有明显优势。

以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围内。

Claims

1.高原山地阴影区冰川识别方法，其特征在于，具体步骤如下：

步骤S1、获取高原山地影像，并对其进行预处理；

2.根据权利要求1所述高原山地阴影区冰川识别方法，其特征在于，所述高原山地影像是Landsat 8OLI影像。

3.根据权利要求1所述高原山地阴影区冰川识别方法，其特征在于，所述步骤S1中对高原山地影像进行预处理，具体为先对其进行辐射定标预处理，再利用FLAASH模型对定标后的影像进行大气校正。

4.根据权利要求1所述高原山地阴影区冰川识别方法，其特征在于，所述步骤S3确定的优势波段为波长小于0.515μm的蓝光波段。

5.根据权利要求1～4任意一项所述高原山地阴影区冰川识别方法，其特征在于，所述步骤S4具体实现过程如下：

增强型雪被指数法计算公式为：

其中，Band_X为阴影区冰川信息识别第X优势波段的DN值或反射率；Band_Y为阴影区冰川信息识别第Y优势波段的DN值或反射率，ENDSI为增强型雪被指数，形成ENDSI灰度图；

6.根据权利要求5所述高原山地阴影区冰川识别方法，其特征在于，所述步骤S43中依据计算结果直方图选择合适的阈值，是选择能够明显区分阴影冰川区与非冰川区的值，即选择计算结果直方图中两个明显的峰值之间的最低值，作为合适的阈值。

7.根据权利要求5所述高原山地阴影区冰川识别方法，其特征在于，所述步骤S42中依据阈值对ENDSI灰度图进行分割依据为：当ENDSI灰度图的某像素点的灰度值大于或等于阈值时，该点为阴影冰川区；当ENDSI灰度图的某像素点的灰度值小于阈值时，该点为阴影非冰川区。