基于高程变化及遥感影像的高亚洲地区跃动冰川数据集

吕明阳，郭华东，闫世勇，李冠宇,3，蒋迪,3，张豪磊，张子彦

1.中国科学院空天信息创新研究院，数字地球重点实验室，北京 100094

2.中国矿业大学，环境与测绘学院，江苏徐州 221116

3.中国科学院大学，北京 101408

引 言

高亚洲地区是指亚洲中部包括青藏高原及周围的高山高原地区在内的海拔在 4000-5000 m以上的高海拔区域，南北方向以喜马拉雅山南麓和天山山脉为边界，东西以青藏高原东端和兴都库什山脉为边界，这些山脉和高原环绕着中国新疆维吾尔自治区和西藏自治区，包括了尼泊尔和不丹全境，以及部分巴基斯坦和印度地区(图1)。21世纪以来，全球气候变化致使高亚洲地区山地冰川退缩剧烈，引起各国学者的广泛关注[1-4]。高亚洲地区的冰川融水是亚洲内陆地区的重要淡水来源，尤其对于中国西北干旱区而言，该地区冰冻圈的变化会直接影响区域内部及下游的生态环境[5-6]。除此之外，冰川的急剧变化也会诱发许多冰川相关灾害，如山洪、滑坡、泥石流、冰川垮塌、冰川跃动等，是各流域内居民生产生活的巨大安全隐患[7-8]。冰川跃动作为一类冰川相关灾害，因其事件突发性、地点隐蔽性、过程难探测性等特点，长期以来是冰冻圈灾害研究的难点。

冰川跃动是介于正常冰川流动和冰崩之间的一种特殊的周期性冰川运动现象。在跃动过程中，冰川泄冰区积累的大量冰会在短时间内迅速传输至冰川下部的冰补给区，多数情况下会引起冰川末端的前移。跃动通常会持续数月至数年，其间冰川的流速可以达到平静期流速的数十至数百倍[9]。虽然跃动冰川仅占全世界冰川数量的1%左右，但它们对于研究冰川作用过程、流动不稳定性等方面有着至关重要的作用[10]。相关研究表明全球的冰川跃动趋于集中发生在一些特定地区，Sevestre and Benn[10]通过总结已有文献记录，获得了包括全球跃动冰川分布情况，形成了基于文献整理的跃动冰川数据集，并以属性信息的形式收录于Randolph Glacier Inventory (RGI) 6.0版[11]。值得注意的是，由于斯瓦尔巴群岛和阿拉斯加-育空等环北极地区研究历史悠久、相关文献较丰富，RGI的跃动记录在这些区域具有很高的可信度。但是高亚洲地区冰川跃动研究起步较晚，相关文献不甚详尽，多数区域仍可能存在未被探明的跃动冰川。对于喀喇昆仑地区和帕米尔地区，有学者进行了细致的跃动冰川的空间观测识别工作[12-14]，但仅为区域性研究成果。许艾文等对喀喇昆仑地区乔戈里峰附近的跃动冰川进行了全面的遥感监测，并提出了导致跃动发生的可能原因[15]。针对一些大型且造成灾害的跃动事件也有学者开展了详细的研究工作，例如2015年的帕米尔地区克拉牙依拉克冰川跃动[16]、2014-2016年喀喇昆仑地区的Hisper冰川跃动[17]。国内相关团队已经开展了高亚洲地区的跃动冰川的识别工作，并已取得相关成果[18]，但目前尚未有覆盖高亚洲地区的跃动冰川数据集发布。

为了填补这一研究数据空白，本团队通过配准及差分两批覆盖高亚洲地区的数字高程模型(DEM)，得到了一期冰川表面的高程变化数据，并结合2000-2016年冰川表面高程变化结果[15]以及20世纪70年代至今的历史光学遥感影像，识别了该地区的跃动冰川。经过重新绘制跃动冰川的最大轮廓，补充和完善名称、地理位置、面积、高程、跃动判断依据等冰川属性信息，最终获得了高亚洲地区的跃动冰川数据集。相比RGI基于文献整理的跃动记录信息，本研究对于跃动冰川的筛选标准是统一的，而且本数据集提供详实的冰川跃动判断依据，以及最大跃动冰川轮廓。本数据集可作为深入研究高亚洲地区冰川跃动的基础数据，也可作为该地区公共基础设施规划和建设的参考资料，还可为预防冰川跃动造成居民生命财产损失提供可靠研究支撑。

基于研究的便利性考虑，依据冰川分布地区地理地形的差异，并结合前人的研究分区情况[2]，本研究将高亚洲拆为8个分区：天山山脉、帕米尔地区、兴都库什山脉、喀喇昆仑山脉、昆仑山脉、青藏高原内部地区、喜马拉雅地区、念青唐古拉山脉(图1)。

1 数据采集和处理方法

1.1 数据源

本数据集的生成主要使用了3种数据类型：光学遥感影像、DEM、冰川编目数据。光学遥感影像主要用于冰川边缘轮廓提取、历史跃动事件的影像验证；DEM主要用于冰川表面高程信息及高程变化提取；冰川编目数据主要用于冰川边界参照、DEM配准辅助。

1.1.1 光学遥感影像

本研究收集了1972-2019年共2904景Landsat系列卫星影像，具体数据情况见表1，影像的筛选原则为在确保数据时间连续性的前提下，尽量选择云量较少、季节性积雪覆盖较少的影像。Landsat系列影像是通过USGS的地球资源观测科学中心(EROS，http://glovis.usgs.gov/)免费获取，为经过辐射校正和几何校正的Level 1级产品。谷歌地球高分辨率历史影像主要用于辅助Landsat影像识别冰川边缘轮廓和验证历史冰川跃动事件。

表1 Landsat系列卫星影像参数及使用的数据量Table 1 Parameters and image numbers of used Landsat satellites

1.1.2 DEM

本研究中使用了覆盖全高亚洲地区的SRTM DEM和AW3D30 DEM两套DEM数据，以及一期DEM差分数据。

SRTM DEM是现存发布最早的覆盖全球陆表的DEM产品，是由2000年2月11日至22日奋进号航天飞机获取的C波段雷达数据处理而得，第三版SRTM DEM最优空间分辨率达1弧度秒(30 m)，其水平精度和垂直精度均优于10 m[19]。本研究选择以SRTM DEM代表2000年高亚洲地区冰川表面高程信息。可从美国宇航局(NASA)的Earthdata网站(https://earthdata.nasa.gov/)获得SRTM DEM。

AW3D30 DEM是基于ALOS PRISM传感器在2006-2011年拍摄的多视影像（约300万景）通过立体像对技术处理后得到的全球陆表DEM产品，空间分辨率为1弧度秒(30 m)，垂直精度约为4.4 m。本研究使用的为2.1版本，相比于ASTER GDEM，AW3D30 DEM的数据精度更高，而且用于生成DEM的影像时间范围(2006-2011)确定，因此本研究选取其与SRTM DEM配准差分，获得2000-2006～11高亚洲地区冰川表面高程变化数据。本数据集中相关高程、坡度、坡向信息也是由AW3D30 DEM 计算而来。AW3D30 DEM 可以通过日本宇宙航空研究开发机构(JAXA)网站(https://www.eorc.jaxa.jp/ALOS/en/aw3d30/index.htm)申请获得。

法国学者Fanny Brun博士利用Ames Stereo Pipeline(ASP)工具处理ASTER Level 1A级的5万多景影像，得到了2000-2016年高亚洲地区时序DEM数据集，通过配准及差分，得到了2000-2016高亚洲地区冰川表面高程变化数据[20]。该DEM差分数据(FB DEM差分数据)的分辨率为30 m，基本能够覆盖高亚洲地区。

本研究将上述2000-2006～11和2000-2016两期DEM差分数据作为识别高亚洲地区跃动冰川的基础数据，用以初步筛选跃动冰川。使用的DEM数据量详情见表2。

表2 高亚洲地区跃动冰川识别使用的DEM数据量(幅)Table 2 DEM used for surge-type glaciers identification in the High Mountain Asia

1.1.3 冰川编目数据

本研究使用的冰川编目数据为RGI 6.0版本。RGI中的冰川轮廓的勾绘主要依据1999-2010年的遥感数据，其首要目的是确保覆盖地球上所有冰川，约198000条冰川总面积为726800±34000平方千米(km2)，面积误差主要来自对单一冰川面积的不确定性评估，以及与其他冰川编目的对比结果，对季节性积雪和冰碛物覆盖的误判是导致误差的主要原因[21]。RGI可在网站http://www.glims.org/RGI免费获取。

由于冰川表面高程变化数据仅用于定性识别泄冰区和冰补给区的冰量得失，因此直接使用 RGI的冰川轮廓辅助DEM配准工作，不需要对冰川轮廓进行修正。而本数据集提供的冰川轮廓矢量文件是以RGI冰川轮廓为基础，依据跃动冰川在历史影像中的最大轮廓重新勾绘而成。

1.2 数据采集方法

本研究识别跃动冰川并生成数据集的步骤如下：(1)以 SRTM DEM 代表 2000年地表高程、AW3D30 DEM代表2006～11年地表平均高程，通过DEM配准得到2000-2006～11年冰川表面DEM差分结果。结合2000-2006～11年和2000-2016年冰川表面DEM差分结果(FB DEM差分数据)，挑选表面存在异常升降现象的冰川为可能的跃动冰川。(2)针对这些潜在的跃动冰川，通过目视解译冰川表面特征变化和末端的演化情况，逐一识别冰川跃动事件，进而确认跃动冰川。(3)最终利用ArcGIS软件勾绘跃动冰川的最大轮廓，计算位置、面积、高程、坡度、坡向等冰川属性信息，合并RGI中相关属性信息，得到高亚洲地区跃动冰川数据集(图2)。

1.2.1 跃动冰川空间观测识别标准

20世纪中期，Meier and Post[22]在对北美洲西部冰川跃动的研究中指出跃动过程中冰川流速会急剧增加，可高出其平静期流速的十倍至百倍；跃动事件后，流速远低于周围非跃动冰川。Meier and Post[22]将冰川末端的快速剧烈前移也作为冰川跃动事件的判断标准之一，但末端位置的变化仅可作为判断参考，因为在后期学者的研究中发现并不是所有冰川跃动都会导致末端的前移[9,14]。此外，因流速快慢变化而形成的地貌特征也可作为冰川跃动识别标准之一，例如水滴状表面冰碛物特征、纵横交错的冰裂隙、翻卷褶皱的冰碛物、蛇形丘等[10,23-24]。以上标准即为Sevestre and Benn[10]在筛选全球冰川跃动记录所用的判断标准。需要注意的是，由于可用文献有限，其冰川跃动记录中并不是所有冰川均包含流速异常的记录。

冰川跃动会将上部泄冰区的大量冰体传输至下部的冰补给区，这种冰川系统的冰量重新分布会直接导致泄冰区高程降低、冰补给区高程增加，部分情况下冰补给区会冲出原有冰川末端边界，导致冰川末端前移和冰川面积增加。跃动事件结束后，泄冰区开始新一轮冰量积累、高程增加，而冰补给区在跃动期间接收的冰体会在平静期内逐渐消融、高程降低。在相同的气候条件下，临近区域正常冰川的冰量变化和流速变化不会出现明显的差异，因此，某一时段内的冰川表面高程变化不仅可用于识别该时段内发生的跃动事件，也可用于识别可能在该时段前发生跃动的冰川。本研究通过对比不同时间泄冰区和冰补给区的高程变化差异，初步挑选可能的跃动冰川。考虑到跃动前后，冰川表面高程变化存在相同的可能，而且不同DEM在不同地区数据质量不同，仅使用一期冰川表面高程数据来判别跃动冰川存在遗漏的可能。因此本研究使用了2000-2006～11年和2000-2016年两期DEM差分结果，以确保跃动冰川初步筛选的完备性。图3a、b和图4a、b给出了kk149035_32和pm150033_17跃动冰川的两期表面高程变化，可见冰川跃动会导致冰川表面高程发生明显变化。目前已有部分区域性研究采用通过DEM差分结果识别跃动冰川，并得到了较好的识别效果[13]。

针对初步筛选的跃动冰川，本研究通过目视解译1972-2019年的遥感影像中与跃动相关的表面特征变化，确认发生过的冰川跃动事件。跃动相关的图像特征变化具体是指与临近冰川相比，目标冰川末端是否经历过快速前移、是否存在由流动快慢振荡导致的冰川表面特征变化(图3a、b；图4c、d)。冰川流速的异常变化也常被作为跃动的识别标准之一[22,25]。在高亚洲地区的冰川跃动过程中，冰川流速急剧增加，会高出正常流速的十倍甚至百倍，并在跃动后明显降低[8,14]，这是唯一可以量化的判别跃动冰川的标准，但考虑到现有覆盖高亚洲地区的冰川流速产品较少，且其时空分辨率尚不能较好满足跃动过程的识别，因此本数据集的生成并未涉及相关的冰川流速数据。对于跃动事件开始和结束时间的判断，则是通过解译历史遥感影像中冰川表面裂隙、冰碛物覆盖、冰川末端等冰川特征开始前移和停止前移的时间而确认。若日后加入较理想的冰川流速产品，满足识别冰川跃动的时空分辨率，将有助于完善本数据集。

1.2.2 数据预处理

SRTM DEM和AW3D30 DEM有不同的投影坐标系统，空间信息的不匹配会使得DEM配准过程失败。因而，在配准前，需将这两套DEM数据以三次卷积内插的方式重采样至30 m，并投影于WSG84 World Mercator坐标系。用于目视解译的Landsat系列卫星影像均为经过正射校正的影像。为凸显冰川边缘，在后期对图像进行了增强处理。

1.2.3 DEM配准

由于传感器不稳定、DEM处理方法局限、实地勘测情况不充分、后期 DEM处理步骤差异等，会导致不同源DEM数据之间存在偏差，主要体现在三个方面：地理定位的偏差、高程扭曲的偏差、卫星轨道模式相关的偏差[26]。本研究采用一种三步修正法，逐步修正AW3D30 DEM和SRTM DEM之间的偏差[26]。在地形配准工作前，冰川覆盖地区的DEM需要使用RGI将其排除在外，批量配准工作通过Ames Stereo Pipeline的dem_align工具完成[27]。

2 数据样本描述

高亚洲地区跃动冰川数据集包含 5个数据文件，其中 3个矢量文件(STG_HMA.shp、STGlikely_HMA.shp、HMA_regionshape.shp)、2个表格文件(Attributes of STG.xlsx、Attributes of STGlikely.xlsx)。

HMA_regionshape.shp为高亚洲地区各分区矢量文件，可在其属性表中查看各分区名称。

每个跃动冰川的最大轮廓均以多边形的形式记录于STG_HMA.shp文件中，并带有相应的属性信息(表 3)。需要说明的是 STG_HMA_ID为本数据集对各跃动冰川的命名，命名规则为XXZZZZZZ_YY。XX为该冰川所处研究分区名称缩写；ZZZZZZ为覆盖该冰川位置的Landsat影像的行列号，即全球参考系World Wide Reference System(WRS)行列号；YY为该行列号区域中跃动冰川的序号。例如kk148035_04代表在喀喇昆仑地区148列35行Landsat影像中识别的第4条跃动冰川。Attributes of STG.xlsx为STG_HMA.shp对应的跃动冰川属性信息表格文件，其中包含对各跃动事件的判断依据。

STGlikely_HMA.shp中记录可能的跃动冰川轮廓。由于在历史影像中未观测到跃动事件，因此其属性信息中未包含表 3中的 Ini_date、Ter_date和 Surge_dur。对于可能的跃动冰川命名规则为likelyXXZZZZZ_YY，具体规则同上。Attributes of STGlikely.xlsx为STGlikely_HMA.shp对应的属性信息表格文件，其中包含判断其为可能的跃动冰川的依据。

表3 STG_HMA.shp文件中的属性信息类别Table 3 Attribute information category of STG_HMA.shp

本数据集确认了高亚洲地区发育的457条跃动冰川，总的跃动冰川面积为42906.46 km2。在研究期内发生过跃动的跃动冰川有362条，其中天山地区32条，帕米尔地区134条，兴都库什地区6条，喀喇昆仑地区127条，昆仑山地区20条，青藏高原内部地区35条，喜马拉雅地区7条，念青唐古拉地区1条。这362条跃动冰川从20世纪70年代至今发生了421次冰川跃动事件。另外还有95条可能的跃动冰川，其中天山地区5条，帕米尔地区5条，兴都库什地区8条，喀喇昆仑地区37条，昆仑山地区13条，青藏高原内部地区23条，喜马拉雅地区4条(图5)。

3 数据质量控制和评估

由于本研究使用的DEM数据和光学遥感影像分辨率较低，对于面积大于5 km2的跃动冰川具有较好的识别效果，但对于面积小于 5 km2的跃动冰川识别效果较差。本研究尽可能多地获取遥感影像，以控制跃动事件起止年份的识别误差在一年之内，但由于可用遥感影像有限和人工识别的误差，对于少部分跃动事件起止年份的识别可能会超过一年。

本研究使用的SRTM DEM水平精度和垂直精度均优于10 m[19]，AW3D30 DEM垂直精度约为4.4 m。在高亚洲的高山高原地区，现有DEM产品的总体精度会相对低于其他地区，但已有研究表明通过配准AW3D30 DEM和SRTM DEM来识别跃动冰川是可行的[13]。为了去除AW3D30 DEM和SRTM DEM之间的偏差，本研究采用Nuth and Kaab的方法，对两套DEM之间的地理定位的偏差、高程扭曲的偏差、卫星轨道模式相关的偏差配准进行配准校正[26]。配准结果的精确与否直接决定了冰川高程变化结果的可信度。假设非冰川区域在研究期内应无明显变化，那么配准后非冰川区域的偏差将直接反映配准的质量。图6为AW3D30 DEM - SRTM和FB DEM差分数据的两期数据在非冰川区域的对比图。经统计，这两期DEM差分结果在非冰川区域的标准差均接近0值。虽然在不同DEM数据在山地区域的配准结果存在较明显偏差，但并不影响对个体冰川表面高程变化的定性识别(图 3、图 4)。

冰川末端的冰碛物覆盖会严重影像冰川边界的识别，而且手动勾绘冰川轮廓也会增加结果的不一致性和不可重复性[28]。本研究使用轮廓测量误差评估方法[29]，分别针对同一冰川和不同大小不同区域的多条冰川，结合Landsat影像和谷歌地球高分辨率影像的特征辨别冰川末端的变化，多次重复描绘冰川边缘轮廓，测量冰川末端变化（表4）及冰川面积（表5），并将结果的标准差平均值作为相应冰川属性信息的误差。最终得到冰川末端变化测量的误差为15.70 m，人工修正冰川轮廓后冰川面积测量的误差近于1.42%。

表4 冰川末端变化测量误差统计表Table 4 Error statistics of glacial terminal change measurement

表5 冰川面积测量误差统计表Table 5 Error statistics of glacial area measurement

为了确保数据集的准确性，我们将本数据集与高亚洲地区相关的跃动冰川数据集进行了对比。需要注意的是，由于不同研究所用方法、数据、研究期跨度存在差异，使得对比工作存在一定困难[30]，而且不同研究记录跃动事件及识别跃动冰川轮廓的方式不同，因此为了便于对比统计，我们整理了多个基于遥感研究得到的跃动冰川数据集，并将其与RGI冰川轮廓对应，统计各研究识别的跃动冰川数量，这其中包括RGI 6.0中记录的跃动冰川，以及多个区域性跃动冰川数据集(表6)。表6中本数据集的跃动冰川数量也为对应的 RGI冰川数量。RGI 6.0中的跃动冰川信息均来自 Sevestre and Benn[10]，其中记录了两类跃动冰川信息：可能的(1 Possible和2 Probable)、确认的(3 Observed)[11]。在高亚洲地区，RGI未记录可能的跃动冰川，仅记录了102条确认的跃动冰川，远低于本研究中记录385条确认的跃动冰川和95条可能的跃动冰川。Goerlich等[13]和Mukherjee等[31]使用了与本研究相似的识别跃动冰川的方法，即判别冰川末端是否前移和冰川表面高程变化，分别识别了帕米尔地区的202条跃动冰川和天山地区的39条跃动冰川。Mukherjee等[31]在天山地区的研究结果与本研究结果(33+5)相似，Goerlich等[13]在帕米尔地区的研究结果多出本研究结果(150+5)约47条，这可能是因为该研究使用的数据分辨率较高、时间跨度较长，对于一些小型跃动冰川和持续时间较长的冰川跃动有较好的识别度。Bhambri等[12]基于1840s至2017年的实地考察数据和卫星数据，通过判别冰川末端、表面流速和表面特征的变化，识别了喀喇昆仑地区的172条跃动冰川，与本研究得到的研究结果(137+37)相似，2组跃动冰川数据在空间分布上有些许出入，这可能是因为Bhambri等[12]的研究时间跨度较长，对于获取卫星数据之前发生的跃动有较好的识别性，而本研究使用的遥感数据较多样，可以较好地识别近年来发生的跃动事件。Yasuda and Furuya[32]通过研究昆仑山部分地区冰川末端和冰川表面流速的变化识别了9条跃动冰川，少于本研究识别的数量(20+13)。以上对比结果表明本数据集有较高的质量和可信度。

表6 与其他数据集识别的跃动冰川数量对比Table 6 Comparison of surge-type glacier numbers between this study and other inventories

4 数据价值

本研究通过空间观测手段，依据冰川表面高程和表面特征变化情况，识别了高亚洲地区的跃动冰川。本数据集可作为深入研究高亚洲地区冰川跃动机理和山地冰川演化过程的基础数据，也可作为该地区公共基础设施规划和建设的参考资料，还可为预防冰川跃动造成当地及下游居民的生命财产损失提供可靠的支撑材料。

5 数据使用方法和建议

本数据集中矢量文件的地理坐标系均为WGS-1984，可用ArcGIS等GIS软件查阅使用。研究人员可根据需要灵活选择部分地区或全部地区的跃动冰川数据，以补充其现有研究；也可针对数据集中记录的跃动事件，开展跃动机理的深入研究。

致 谢

本项目受中国科学院战略性先导科技专项（A类）（XDA19070202）资助。感谢中国科学院空天信息创新研究院刘广研究员提供的宝贵意见和建议，感谢USGS的EROS提供的Landsat系列卫星影像，感谢NASA和JAXA提供的SRTM DEM和AW3D30 DEM，感谢Fanny Brun博士提供的2000-2016高亚洲地区冰川表面高程变化数据，感谢GLIMS提供的6.0版RGI数据。

数据作者分工职责

吕明阳（1991—），男，吉林省松原市人，博士，助理研究员，研究方向为冰冻圈遥感。主要承担工作：数据收集与处理，跃动冰川筛选与识别，跃动冰川最大轮廓勾绘，论文撰写。

郭华东（1950—），男，江苏省徐州市人，博士，研究员，研究方向为遥感理论与数字地球。主要承担工作：研究方案设计，论文修改。

闫世勇（1982—），男，江苏省徐州市人，博士，副教授，研究方向为雷达遥感应用研究。主要承担工作：研究方案设计，数据处理，论文撰写及修改。

李冠宇（1999—），男，山东省滕州市人，硕士研究生，研究方向为极地遥感。主要承担工作：冰川属性信息整理，数据集信息校对。

蒋迪（1997—），女，陕西省西安市人，博士研究生，研究方向为极地遥感。主要承担工作：跃动冰川最大轮廓勾绘，数据集信息校对。

张豪磊（1998—），男，浙江省绍兴市人，硕士研究生，研究方向为雷达遥感。主要承担工作：冰川属性信息整理，数据集信息校对。

张子彦（1999—），男，江苏省常州市人，硕士研究生，研究方向为雷达遥感。主要承担工作：冰川属性信息整理，数据集信息校对。