任彦润,张耀南,康建芳
1. 中国科学院西北生态环境研究院,科学大数据中心,兰州 730070
2. 中国科学院大学,北京 10004
数据库(集)基本信息简介
在全球气候变暖的背景下,近几十年我国冰川面积总体呈现加速退缩的趋势。对比两次冰川编目数据可知,近 30年来我国冰川面积约减少 17%[1],区域水文过程和冰湖变化受冰川环境影响较大。中巴经济走廊范围内地势险峻,地质构造复杂,水文特征多变,固态降水丰沛,现代冰川发育充分,冰川变化十分活跃,冰川阻塞湖分布广泛,冰川的进退与消融极易引发各类冰川灾害,冰湖溃决引发的洪水灾害是研究区典型的灾害类型之一。
冰湖是由冰川作用形成的湖泊或以冰川融水为主要补给源的湖泊[2],冰湖的发育和变化过程与气候环境的变化、冰川物质平衡、冰川的温度状态及水力特征等有着密切的关系。现代冰川是本区内高山冰雪融水和冰湖溃决洪水的主要水源,其中气候变化导致的冰湖上游来水增多、冰坝自身消融变形及冰滑坡、冰崩、滑坡、泥石流进入湖泊,都有可能造成堤坝溃决,形成冰湖溃决型洪水。冰湖溃决引发的洪水及其次生灾害是高山冰川区常见的灾害之一,冰川阻塞湖和冰碛阻塞湖是最常发生溃决的两类冰湖。
中巴公路大多路段均位于河谷之中,两岸高山及沟谷之中冰川、多年积雪密布,冰川活动是中巴经济走廊地质灾害的主要激发因素之一。针对冰川活动引起冰湖灾害及其可能对中巴经济走廊形成的危害,本文利用来源于美国地质勘探局数据中心的Landsat 8 OLI遥感影像,基于地理信息系统平台的支持,采用面向对象的分类方法,对中巴经济走廊范围,包括洪扎河流域、努布拉流域、盖孜河流域、克勒青河流域等多个典型流域进行冰川及冰湖分布的分类和提取。鉴于冰川与冰湖存在共生关系,有时两者的空间位置还会发生重叠,再加上冰湖与山体阴影具有相近的光谱特征,精确提取冰川和冰湖的边界存在一定困难。高山冰湖的提取由于自身理化特征的复杂性及周边背景的影响,其提取工作大多数采用传统野外监测结合人工解译的方式实现,需要进行大量、长时间的分析和处理。本文基于面向对象的分类方法,以多尺度分割算法对高分辨率影像进行分割,以雪盖指数法、归一化水体指数法对典型地物进行分类和提取,剔除干扰因素,实现了遥感影像中冰川、冰湖信息的自动提取。对研究区内冰川和冰湖的空间分布及变化特征进行长期、定时的监测可为中巴经济走廊的进一步建设提供数据支持,对区域水资源变化、冰湖溃决危险性评估的科学决策也具有重要意义。
中国科学院寒区旱区工程与研究所组织和参加了多次冰川、冰湖的科考调查,参与编制了大量冰川编目数据,如 2006年开始的中国第二次冰川编目工作、2015年中国喜马拉雅山区冰湖遥感调查与编目等,并将其共享在了寒区旱区科学数据中心。在中巴经济范围走廊内,基于交通、水利等建设工程的喀喇昆仑山区冰川灾害科考也已开展了大量调查和研究工作。
我国的冰湖研究和编目工作大都集中在西藏喜马拉雅山区,西藏是我国冰湖分布最多的地区。1987年我国科学家便联合尼泊尔进行了该区域冰湖溃决洪水的考察及冰湖编目工作,并基于此预测了冰湖溃决的规模和波及范围,成果显著。该区域冰湖类型主要为冰碛阻塞湖。冰川阻塞湖溃决引发的溃决型洪水在我国境内主要发生在喀喇昆仑山北坡,该地发生的洪水对叶尔羌河的水资源变化和研究有重大影响,但尚未见针对本区域的完整冰湖数据,本次冰川和冰湖数据集的制备一定程度上弥补了冰湖编目的空缺。
Landsat 8是Landsat系列于2013年2月发射的最新卫星,携带了TIRS(Thermal Infrared Sensor,热红外传感器)和OLI(Operational Land Imager,陆地成像仪)两个传感器,在原有基础上新增了一些波段,并对部分波段进行了调整。Landsat 8遥感影像有像幅面积大(成像宽幅185 km×185 km)、获取周期短(16天)及波段信息丰富(9个波段)的特点,可以通过不同波段组合形式进行影像分析。Landsat 8 OLI影像全色波段最大空间分辨率为15 m,其他波段空间分辨率为30 m。通过对多光谱影像不同波段的组合可以凸显用户感兴趣的地物类型,根据冰雪和水体对不同波段的反射率特征,可选择第2、5、6波段的RGB组合进行判读和目视解译。本次冰川冰湖分布数据集的制备主要基于Landsat 8 OLI遥感影像,共选取从2013-2017年共80幅影像,基本覆盖了中巴经济走廊周边绝大部分冰川分布区,数据集所用Landsat影像的详细信息见表1,Landsat影像在整个研究区中的分布见图1。
图1 2015年数据集所用遥感影像在研究区中的分布图
在遥感影像的解译过程中,由于“同物异谱”“同谱异物”现象,高山冰川和冰湖的提取会受到云雪和山体阴影的影响。为排除云、雾、雪、山体阴影等对遥感影像的解译精度造成影响,优先选取年内 7-10月中云量较少的遥感影像,此时段影像占总影像总数的70%。由于冰湖受年内气温和降水的影响较大,故选取影像时将备选影像时段扩展为研究区气温较高、降水较为丰富的6-11月。此时段内冰湖面积相对较大、水位较高,便于其提取和识别。USGS(United States Geological Survey,美国地质勘探局)网站提供的所有 Landsat 影像均已经过系统辐射校正、地面控制点几何校正及基于DEM的地形校正,因此在本次数据预处理阶段中没有再进行影像校正工作。
表1 数据集用到的Landsat影像信息列表
本次冰川冰湖分布数据集的研究范围如图2所示。根据流域规模、水文变化特征、冰川空间分布等特点选取喀喇昆仑山区洪扎河流域、努布拉流域、盖孜河流域和克勒青河流域为典型研究对象。
洪扎河流域位于喀喇昆仑山脉和兴都库什山脉交界处,地势险峻,冰川密布,其中80%为山岳冰川[3],流域及其周边分布有帕苏-慕士塔格山活动性冰川群,包括古尔米特冰川、固尔金冰川、帕苏冰川、巴托拉冰川等多条北西-南东走向的冰川,由冰川变化导致的冰湖溃决及冰川泥石流等危害较大。据统计在洪扎河谷分布有 110个面积动态变化的冰湖[4],在对洪扎河流域冰川的观测中还发现相邻冰川在相同时期进退不同步的情况[5],可见对其进行定期动态的监测有助于对冰川进退规律的认识和冰川灾害的防治。洪扎河谷附近的帕苏-慕士塔格冰川群是中巴公路沿线活动性冰川群中最为活跃的部分,我国自 1974年就开始了对喀喇昆仑山南坡洪扎河谷的冰川泥石流灾害的考察和研究。近期发生的冰川灾害中,以2010年由于巴基斯坦北部山区地震引发山体滑坡导致洪扎河被堵,进而形成了世界上最大的堰塞湖最为严重。
努布拉流域属于印度河流域,位于喀喇昆仑山中部,平均海拔高,面积大约 5×103平方千米,受西风环流的影响降水丰沛,冰川覆盖率接近50%,河流主要由冰川融水补给,流域内分布有被称为“克什米尔地理中枢”的锡亚琴冰川,周边分布有狮泉河、塔什库尔干等气象观测站。每年都有大量的喀喇昆仑山冰川融水流入印度河支流,努布拉河就是其中之一,锡亚琴冰川的融水是努布拉河的主要补给来源。现有的努布拉流域灾害研究相对较少,对努布拉流域冰川的长期监测对区域水资源利用及灾害研究有实际意义。
图2 研究区范围地理位置示意图
盖孜河流域地处新疆西南部,帕米尔高原东缘,受西风环流影响,冰川规模与地貌发育充分。盖孜河由慕士塔格、公格尔等高山的冰雪融水形成,是喀什噶尔河水系的源流之一,设有喀拉库里水文站、克勒克水文站、维塔克水文站及喀什、塔什库尔干等气象观测站。1984年我国就对帕米尔东缘盖孜河谷的冰川泥石流进行了相关科学调查。流域内分布有克拉牙依拉克、其木干等多条冰川,其中克拉牙依拉克冰川是公格尔山北坡最大的现代冰川,2015年曾被观测到冰川跃动事件。
克勒青河流域位于塔里木河河源区,坐落在喀喇昆仑山分水岭北侧,处于西风环流控制区[6],冰川分布高度密集,克勒青河由冰川融水形成,是叶尔羌河的重要支流之一,流域内分布有我国境内最长冰川音苏盖提冰川、特拉木坎力冰川等多条大型冰川,曾多次被记录和观测到冰川跃动现象[7-8],是我国冰坝湖溃决洪水的主要研究区域,周边分布有塔什库尔干、吐尔尕特和乌恰三个气象观测站。1985年我国就对开始对喀喇昆仑山北坡克勒青河谷冰川洪水进行考察科考调查,为叶尔羌河水资源的开发和利用提供了资料支持。
1.3.1 数据处理方法
为了快速、准确地提取目标类别地物,选择使用面向对象分类方法对影像进行处理。按照分类过程基本单元的不同,遥感影像的分类分为基于像元的分类方法和面向对象的分类方法。基于像元的分类方法主要以光谱特征为分类依据,难以避免由于“同物异谱”和“异物同谱”带来的错分,分类精度不高。面向对象的分类方法以对象为影像分类的基本单位,充分利用了影像的光谱信息,通过形状、纹理、空间结构等特征及上下文信息进行类别划分[9-11],分类效果较好。
遥感影像的处理基于影像光谱信息,由于影像间生成时间不同导致的光照强度差异、卫星系统偏差等问题,同一地物的光谱在不同影像中可能存在差异。为避免此种情况,选择在各幅影像完成所有数据处理和分类工作后再行拼接。
数据处理流程如图3,具体步骤如下:
图3 数据处理流程
(1)研究区的选定:根据中巴经济走廊的覆盖范围及冰川分布的空间位置确定研究区范围,获取相应Landsat影响的条带号和行编号。
(2)数据的选取和下载:研究数据下载于美国地质勘探局网站,选取夏、秋季少云的影像,减少干扰因素。
(3)遥感影像的预处理:由于选取的Landsat 8 OLI遥感影像产品级别为L1,故省略大气校正和几何精校正的步骤,直接对影像的全色波段和多光谱影像进行融合。
(4)影像融合:融合多光谱和全色波段,在保持多光谱特性的基础上提升空间分辨率,达到提高地物提取精度的目的。
(5)分类标准的确定:选取需要提取的地物类型并对其类别特征值进行统计和分析,确定具体的地物分类标准、调试阈值范围。
(6)面向对象分类方法:结合分类标准确定最优的分割参数组合,在此基础上进行面向对象的分类,实现对遥感影像的自动解译;由于冰川和冰湖在理化特性上的特异性,无法确定一个统一提取标准,在分类和提取过程中采用全局阈值提取和局部优化结合的方法。
(7)结果导出和验证:结合DEM数据,在分类结果导出时剔除坡度不符合水体特征的部分(坡度大于5°),避免了山体阴影和水体错分的情况;根据冰湖的定义,以本次冰川的分类结果做缓冲区分析并结合水系分布图,去除不符合冰湖空间分布规律的水体和天然水系;导出分类结果;通过验证样本对分类结果进行评价。
(8)结果订正:对分类结果不佳的地区结合多源数据以目视解译的方法对结果进行订正和修改。
1.3.2 山体阴影和云雪影响的去除
山体阴影和积雪范围在夏季影像中的面积最小,选择夏季的遥感影像可在最大程度上消除山体阴影和薄雪覆盖带来的影响,是文中目标地物提取的最佳影像来源。由于综合考虑了云量等影响因素(云量信息见图4),本文的遥感影像选取范围扩展到了6-11月,结合太阳高度角提取冰川边界可以最大程度剔除山体阴影的影像,对其中非夏季受积雪影像较大的影像,在进行后续地物分类和提取过程中将结合目视解译来进行冰川和积雪的边界处理,将积雪的影响降到最小。ENVI提供的去云插件haze tool对去除分类过程中少量云给影像带来的影响效果较好。
图4 文中所用遥感影像的云量信息图
结合Landsat 8 OLI影像获取时的太阳入射角度进行分析,可消除山体阴影对目标地物提取带来的影响。山体阴影在不同太阳高度角下的面积和分布存在差异,读取遥感影像成像时的太阳角和方位角,结合地形坡度图生成地形晕渲图,可获取山体阴影可能的分布范围(山体阴影一般存在于巨大山体山脊线的阴面,地形晕渲值一般为0);Landsat影像中水体与山体阴影存在光谱特征上的相似性,在结果导出时可通过由DEM数据生成的坡度图来进行山体阴影的剔除,当坡度大于5°便认定其为山体阴影并进行删除操作。
研究区的大部分冰川存在较厚的上覆积雪,体现在光谱特征上为光谱特征单一且光谱范围较窄,对冰川边界的提取有利。为了去除季节性积雪对冰川提取的影响,对本文中非夏季时段获取的数据会结合Google Earth原始影像对部分区域进行人工勾绘。中分辨率单一时相的影像对冰川和多年积雪的区分存在难度,在气象资料充足的情况下,结合气象数据可确定遥感影像中不受积雪影响的冰川边界,为最大程度区分积雪和冰川范围提供依据。
1.3.3 冰川和水体的提取方法
提取雪冰覆盖地表的常用方法有比值阈值法、雪盖指数法(Normalized Difference Snow/Ice Index,简称NDSI)等,比值阈值法基于雪冰对短红外波段的强吸收及对可见光波段的强反射特性进行地物区分,雪盖指数法将雪冰地表对可见光波段和近红外波段反射率做归一化处理提取目标地物对象。本文选取归一化雪盖指数法进行雪冰地表和非雪冰地表的区分,选取Landsat 8 OLI影像的可见光绿色波段和短红外波段参与计算,计算公式为(1):
式中,NDSI为积雪指数,BandGreen为 Landsat 8 OLI影像的可见光绿色波段,BandSWIR1为Landsat 8 OLI影像的短红外波段。
按徐演的说法,这个版本本来要作为“云南民族民间文学丛书之一”出版,但因为这时由宣传部领头的,具体由一群学生组成的调查又已经启动,所以出版暂停。徐嘉瑞甚至“拿出这份由出版社已经打印成校样的整理稿,以完全无私的精神,无条件地全部交给了学生们”。
提取水体的常用方法有单波段阈值法、多波段阈值法(水体指数法、谱间关系法)等,单波段阈值法比较简单易行但精度有限,本文中选取以水体与其他地物在各波段光谱值之间的差异为分类依据的水体指数法,通过归一化水体指数(Normalized Difference Water Index,简称NDWI)提取雪冰地表中的水体,由于水体在近红外波段的反射率几乎为零,选用Landsat 8 OLI影像的可见光绿波段和近红外波段参与计算,计算公式为(2):
式中,NDWI为水体指数,BandGreen为Landsat 8 OLI影像的可见光绿色波段,BandNIR为Landsat 8 OLI影像的近红外波段。
对比雪盖指数法和归一化水体指数法计算所得灰度图的统计特征,确定了最佳的阈值提取范围。NDSI大于0.05则为雪冰覆盖区,在雪冰地表中NDWI大于0.15划分为水体。由于本文选用的遥感影像中地形地貌条件、生成时间存在差异,仅使用NDSI、NDWI等指数作为提取依据并给定统一的阈值无法有效处理所有的数据,因此在具体处理过程中针对部分信息提取效果不佳的影像对阈值和方法进行了相应调整,以期获取研究区内冰川、冰湖相对较为准确的边界信息。
1.3.4 分割与分类
面向对象的分类主要流程包括影像分割和分类两部分。影像分割是面向对象分类的第一步,通过比较试验效果,选定多尺度分割为影像分割方法。作为面向对象分类中较为常用的分割方法,多尺度分割综合考虑影像的光谱、形状、纹理信息进行自下而上的区域合并,以保证对象间异质性最小和对象内同质性最大为前提,合并异质性小于设定阈值的相邻像元或“小对象”[12-14],分割形成的多边形不仅包含了被合并像元原有的光谱信息,还形成了形状、纹理及空间位置等信息。
面向对象的分类过程技术路线如图5,具体步骤如下:
图5 面向对象的分类过程技术路线
(1)多尺度分割及参数选取:对融合后的影像进行分割。结合经验和分割实验结果建立合适的分割规则,调整影像的分割尺寸并选定分割参数。分割尺度和分割参数的选择会直接影响分类的精度,根据Landsat影像的分辨率和本次分类的精度要求,确定分割尺度(Scale)为200;在参照前人经验的基础上,多次尝试后确定分割参数:紧凑度(Compactness)为0.5、形状参数(Shape)为0.2,得到了较好的分类结果。
(3)计算机自动分类:根据确定的分类规则按顺序进行初步分类和提取,雪冰地表的分类规则为NDSI>0.05,冰川为NDSI>0.05且NDWI≤0.15,冰湖为NDSI>0.05且NDWI>0.15。根据以上初始阈值进行的全局阈值分割会提取研究区中所有可能的水体,其中包括部分融化的冰川、山体阴影等,以近红外波段NIR<0.15,短波红外波段SWIR<0.05为条件可去除结果中的融化冰川。进行局部优化,剔除误判的水体,以各冰湖为局部分割单元,根据NDWI值的统计特征,冰湖的直方图为双峰分布,以此为判别标准。如果是单峰则为非冰湖,认定该单元为非冰湖并剔除。如果直方图符合双峰分布准则,则认定该单元为冰湖并确定冰湖边界。统计各个冰湖单元内的坡度及晕渲值均值,剔除坡度均值不在水体坡度值范围内的冰湖单元,剔除晕渲值在山体阴影晕渲值范围内的冰湖单元。
(4)合并、去除小图斑:将相邻同类别对象合并为一个整体;将面积较小的图斑与邻近的同一类别的图斑进行合并,如果没有同一类别则与其他地物类别的对象合并;参照前人在基于Landsat等中等分辨率遥感数据制备冰川编目时的标准并结合影像像元实际面积,选取以 0.01km2作为最小冰川面积阈值,以0.001 km2作为保留冰湖的最小面积阈值。
(5)分类后的评价:验证和评价分类结果。
研究范围内包含冰川、积雪、水体等多种高寒山区典型地表类型,以冰川和冰湖为目标地物类别进行提取,本数据集中的地物类型和范围采用GCS_WGS_1984投影坐标,参照现有冰川和冰湖的编目规范,以国际通用的冰川编目方法进行冰川边界的提取和冰川属性的计算,文件中包括 FID、形状类型(Shape)、GLIMS 编码(GLIMS_ID)、冰川名称(Glc_Name)、面积(area)、周长(perimeter)、类别(Class_name)、质心位置坐标(x-coordina、y-coordinate)等信息。2015年中巴经济走廊典型区域冰川冰湖分布如图6。
在数据制备过程中,通过ARCGIS自带的工具对属性信息中的面积和周长进行统计计算;按照GLIMS编码方法对冰川和冰湖进行编码,并补充了部分冰川名称。具体以冰川或冰湖质心的经纬度坐标,按照以下方式编码:
其中G代表冰川,GL代表冰湖,n和m可通过质心小数格式得经纬度乘以1000后四舍五入得到,N与S分别代表南半球和北半球。
面向对象分类方法在本次数据制备过程中适用性较好。由于面向对象分类过程中考虑了纹理等信息,所以对形状规整的地物分类效果尤佳,在水体地物的提取中表现较好,它在对冰川的提取中可有效抑制云雪干扰,提取精度较高,分类结果的完整性和均质性有了提升[15-16]。文中所用方法总体分类效果较好,但是也存在雪盖指数法中将部分山体阴影、河道或水体错分为雪冰地表的情况,归一化水体指数法也并不能完全将水体与雪冰地表区分开,部分水体被错分为冰川。
图6 2015年中巴经济走廊重点区域冰川冰湖分布图
非监督分类可以在缺少先验地物样本、没有人为干涉的条件下自动进行识别和分类,是遥感影像自动提取处理的主流方法。监督分类通过选取具有代表性和确定性的地表点为训练样本,依据样本的特征参数建立判别函数并以此对像元进行类别判断,后续还可以通过检验和评估控制样本来提高分类精度,是地物分类中公认精度较高的方法。为了评价面向对象方法在分类过程中的适应性和精确程度,以监督分类中最大似然法的分类结果结合目前公认分类精度最高的目视解译分类结果形成标准数据,通过计算混淆矩阵(表2)、生产精度、用户精度、总体精度和Kappa系数等进行精度评价,结果如表 3。从表中可发现面向对象的分类方法对冰川、水体的分类精度较高,满足用户对Landsat遥感影像自动解译的需求,总体分类精度0.9072165,Kappa系数0.8088042,在Landsat 8遥感影像的地物分类过程中适用性较好。本文研究区部分与中国第二次冰川编目数据重合,故选用同区域已有冰川分布数据进行针对性的冰川类别精度评价,结果如表 4。从表中可发现本文的分类方法对冰川的分类精度较高,与冰川编目的一致性较好。由于冰湖季节性变化强烈,且不同类别冰湖的变化规律存在年际、季节性差异,本文选用夏秋两季遥感影像进行拼接虽然适用于冰川数据,但对冰湖数据的提取存在较大影响,在接下来的工作中将结合重访周期较短的遥感影像,以多源数据进行冰湖的信息提取。
表2 面向对象分类的混淆矩阵
表3 单一类别精度及总体精度分析
表4 以中国第二次冰川编目数据进行的冰川类别精度分析
本文研究范围内冰川灾害分布较广,而且类型众多、危害性较大,对当地冰川和冰湖进行空间位置和分布情况的定时监测对冰川灾害的研究和预防具有重要意义。喀喇昆仑山在全球冰川气候变暖的大背景下,在部分地区存在冰川萎缩停滞甚至前进的现象,冰川运动活跃,选取典型区进行研究对了解喀喇昆仑山区的气候变化有重要意义。中巴经济走廊附近地质灾害类型多、分布广泛,且大都为冰川活动激发的灾害,相关观测和研究可为中巴公路的防灾减灾工作提供理论和数据支持。在西风气流和印度洋季风同时作用下的喀喇昆仑山和帕米尔东段,冰川消积量较大,洪水、冰川泥石流等冰川灾害频发,但由于地形复杂且环境恶劣,实地观测数据较少,通过遥感手段对其进行检测十分必要。
Landsat影像的地面分辨率已经难以满足现在对地表精细化研究的不断发展需求,但作为地表类型数据制备中的参考数据依旧具有不可替代的地位。本文通过影像的选择、多尺度分割、地物提取函数的确定等自动解译了冰川、水体等地物类别。相比传统方法,面向对象的分类方法在保证解译时效性的基础上提高了解译精度。制备的数据集针对中巴经济走廊冰川分布范围进行了基于面向对象分类方法的特定类型的地物提取,利用优化后的分割参数及分类标准,得出了总体分类精度0.9072165,Kappa系数为0.8088042的分类结果。
面向对象的分类方法中,分割参数的选取对分类结果的影响较大,不同类型遥感影像、目标地物所对应的最优参数组合是不同的,由于缺乏对分割效果统一的评价指标,通过人工目视判定分割效果给分类带来了一定的主观影响,接下来应着重关注分割效果评价指标的建立,提升面向对象分类方法的应用范围。由于地形地貌以及季节性变化特征的差异,没有一个统一的方法可以处理所有的数据并获取好的结果,冰川形态和特征上的差异决定了单一阈值无法在大范围内获取较好的提取结果,在具体处理过程中除了使用NDSI、NDWI等指数作为提取依据,针对某些影像需要使用阈值分割与分类算法结合的信息提取方法。
验证分类精度的最好办法是将实地考察点作为样点,但本次验证样本的地物类别来源于基于目视解译的监督分类结果,难免存在一定偏颇。
本数据集基于Landsat 8 OLI影像获取了2013-2017年中巴经济走廊典型区域冰川和冰湖的分布数据,在面向对象方法的基础上进行地物分类,为冰川变化监测及当地水文变化研究提供了基础数据,是当地冰川灾害研究重要的参考数据。数据集保存为矢量SHP格式,ArcGIS、QGIS、ENVI、ERDAS等常用的GIS与遥感软件均支持该数据的读取和操作。
致 谢
感谢美国地质勘探局数据中心USGS提供的Landsat数据,感谢寒区旱区科学数据中心提供的中国帕米尔高原第二次冰川编目数据和巴基斯坦冰川编目数据。