1989–2020 年基于Landsat 的东帕米尔高原冰川运动速度数据集

张震，黄丹妮，陆艺杰，张莎莎

1.安徽理工大学空间信息与测绘工程学院，安徽淮南 232001

引 言

冰川不仅是重要的水资源，而且是气候变化的敏感指示器。在气候变化的影响下，冰川发生快速响应，主要体现在形态或厚度变化上，众多研究者对此开展了大量研究[1-3]。总体来说，大多数冰川处于减薄或退缩状态，但是冰川变化具有差异性，即不同地区或者同一地区不同性质、不同规模的冰川对气候变化的响应时间和程度不同[3]，具体表现为同一地区相同时段内冰川退缩程度不一，少数冰川甚至处于稳定或者前进状态，而不同地区的冰川变化差异性更大[4]。然而当前对于冰川变化差异性机理的认识还不够清晰，为了量化冰川响应差异，需要进行一定的动力学模拟，但是在动力学模拟过程中又缺乏区域尺度、较长时间的冰川运动速度、物质平衡等数据，传统的定位观测又很难获取这些数据。遥感技术由于具有方便快捷和空间覆盖范围大等特点，近些年在区域尺度冰川参数获取方面得到了大量的应用研究。但研究大多数集中在冰川面积、高程变化和物质平衡上，对冰川运动速度的研究，尤其是对流域尺度或者山系尺度的冰川运动速度的研究总体上还比较薄弱。冰川运动控制冰量传输变化，进而影响冰川长度、面积、物质平衡及厚度变化，冰川运动速度是认识冰川性质、量化冰川变化差异性的关键因素，也是冰川动力学模型的重要参数[5]。此外，冰川运动数据为冰川跃动、冰崩等冰冻圈灾害预警及气候变化研究提供重要依据。因此，基于遥感技术获取区域尺度的冰川运动速度数据集，对冰川研究至关重要。

东帕米尔高原（图1）是我国西部巨大山系，北依天山，西邻塔吉克斯坦，南靠昆仑山，东向塔克拉玛干沙漠，境内有慕士塔格山、公格尔山、昆盖山、外阿赖山和萨雷阔勒岭等山脉，是中巴经济走廊的穿越区。由于其极高海拔而发育众多规模较大的冰川，是下游喀什地区水资源源区，对未来社会经济发展、生态文明建设都具有重要意义。根据中国第二次冰川编目数据[6]，东帕米尔高原共有冰川1265 条，面积2054 km2。近几十年来，该地区冰川变化空间差异显著[7]，具体表现在不同流域的冰川变化差异性显著，即使是同一流域或同一山区，甚至相邻的两条冰川变化呈相反的趋势[4]。冰川变化时间上也具有一定差异性：部分研究者认为该地区冰川近些年呈现正物质平衡现象[4]，此现象被称之为“帕米尔–喀喇昆仑山(冰川)异常”[8]，而2000 年以前却呈现物质亏损状态。虽然对于整个地区是否都存在类似的状态目前还存在一定的争议，但是部分冰川物质亏损减缓或者转为正平衡的趋势是显著的[9]。因此，东帕米尔地区是认识冰川变化差异性研究的典型区域，该地区冰川运动速度的研究对于冰川变化差异性的认识具有重要作用。此外，东帕米尔地区还发现了多条跃动冰川，比如2015 年克拉牙依拉克冰川发生跃动[10]，给人类生产和生活带来巨大破坏，并且引起了新闻媒体、政府和科研工作者的高度重视。冰川运动速度是冰川跃动等相关灾害模拟和预测的重要参数，对于冰川跃动机理认识及预警体系建立具有重要意义。区域尺度冰川运动速度一般基于遥感技术实现，可使用基于光学影像或者SAR 影像的特征匹配法和雷达干涉测量法，其中雷达干涉测量法只适合短时间内冰川运动速度监测，主要应用于冰盖运动速度的提取[11-12]，对于山地冰川，因其规模小和地形复杂导致的失相干问题严重而应用较少。特征匹配法可以获取长时间段冰川运动速度，并且在山地冰川上得到了很多尝试性研究[13-15]。在SAR 领域，特征匹配法又称之为偏移量追踪法。SAR具有全天候的特点，但是历史资料数据不如光学遥感影像丰富。Landsat 历史资料丰富，并且可以免费获取，在冰川运动速度研究中取到较好应用。因此，本研究收集1989–2020 年Landsat TM、ETM+、OLI 影像，基于特征匹配法提取了东帕米尔高原冰川运动数据，为东帕米尔高原冰川研究提供参考。

图1 东帕米尔高原地理位置示意图

1 数据采集和处理方法

1.1 数据源

本研究以影像无云少雪为原则，从美国地质调查局（USGS）网站（https://earthexplorer.usgs.gov/）上选择收集了1989–1999 年和2008–2009 年的Landsat TM、1999–2003 年的Landsat ETM+和2013–2020 年Landsat OLI 遥感影像数据（轨道号149/33 和150/33，表1）。选择两期间隔一年左右的影像用于冰川运动速度的提取。Landsat TM 主要使用红光波段，分辨率为30 m；Landsat ETM+和OLI 影像主要使用全色波段，分辨率为15 m[16]。此外，中国第二次冰川编目数据[6,17]用于裁剪提取冰川区域的运动速度。

表1 研究采用的数据

1.2 数据处理方法和流程

基于光学遥感影像的冰川运动速度的提取通过追踪不同时相影像之间冰川表面特征的位移量获得[18]。其核心思想是，先对两景影像进行精确配准；然后在两景配准后的影像上，根据合适大小的窗口和步长，在一定搜索区域内移动窗口进行互相关计算。当窗口内相关系数最大时的两个点即为同名点，计算同名点之间的偏移量。重复以上步骤，即可获取两幅图像重叠区域在间隔时间内的位移。位移分为东西向（EW）和南北向（NS）位移，从而获取总位移，根据间隔时间即可获取运动速度[19]。特征匹配的方法有很多，其中基于COSI-CORR 软件包的频率域变换的影像互相关方法，能够实现比较精确的亚像元配准和影像互相关，在阴影区和积雪区也能取得较好效果[20]，目前广泛应用于冰川运动速度的研究[13-15]。

本研究使用的Landsat 数据由USGS 采用的GLS2005 全球陆地控制点系统进行了正射校正，校正精度为半个像元左右，部分影像甚至达到1/6～1/10 左右的精度[20]。本研究采用Landsat TM 的第3 波段（红光波段，空间分辨率为30 m）、Landsat ETM+的第7 波段（全色波段，空间分辨率为15 m）和Landsat OLI 第8 波段（全色波段，空间分辨率为15 m），在COSI-CORR 计算相关系数时采用频率域算法，参考窗口设为128，搜索窗口设为32，获得的位移数据有3 个图像组成：东西向位移（EW）、南北向位移（NS）以及信噪比（SNR）。本研究采用SNR≥0.9 部分并剔除云和阴影覆盖部分[10]，得到整幅图像的运动速度（V），然后利用中国第二次冰川编目数据对V 进行裁剪，从而获得冰川运动速度数据（glacV）（图2）。

图2 数据处理流程图

2 数据样本描述

本数据集的数据均存储为GeoTIFF 格式，数据类型为32 位浮点型，分为23 个文件夹存储，文件夹命名规则为YYYY_YYYY，YYYY 表示年份，前者为前一年，后者为后一年；每个文件夹至少有一组数据，L5 代表基于Landsat 5 TM 影像提取，L7 代表基于Landsat 7 ETM+提取，L8 代表基于Landsat 8 OLI 影像提取。YYYYMMDD_YYYYMMDD 代表两个时间间隔，MM 代表月，DD 代表日，后缀SNR 代表信噪比，后缀EW 代表整幅图像东西向位移（单位m），后缀NS 代表整幅图像南北向位移（单位m），后缀V 代表整幅图像的运动速度（单位md-1，为了剔除噪声影响，本结果剔除了信噪比小于0.9 的数据[10]，如果读者需要修改信噪比阈值，可以根据我们提供的SNR 数据进行重新处理），后缀glacV 代表冰川运动速度（根据中国第二次冰川编目数据对V 进行裁剪获取，单位md-1），冰川运动速度结果样本显示如图3–5。图3 为2015 年4 月13 日至2016 年5 月1 日慕士塔格山区冰川运动速度示例图，图4 为2016 年5 月1 日至2017 年5 月20 日公格尔山区冰川运动速度示例图。图5 为科克萨依冰川沿主流线运动速度变化情况。

图3 慕士塔格区域冰川运动速度

图4 公格尔区域冰川运动速度

图5 科克萨依冰川沿主流线运动速度变化情况（左图显示主流线位置，右图为沿主流线冰川运动速度变化情况）

3 数据质量控制和评估

在天山托木尔峰科其喀尔冰川野外实测中证明了本研究方法的可靠性[13,21]。由于东帕米尔高原野外实测困难，本研究使用间接的方法评估。冰川表面运动一般明显快于无冰区稳定地形表面的构造运动，因此本研究假设无冰区稳定地形无位移变化，即运动速度真值为0，那么通过遥感影像提取的无冰区稳定地形的误差可当作冰川区域运动速度的误差[22-23]（表2），利用公式（1）–（3）评估[24]：

式（1）中，eoff为无冰稳定地形区域的运动速度误差，MED表示无冰稳定地形区域运动速度平均值，SE为无冰稳定地形区域的运动速度标准误差，SE计算公式为[24]：

式（2）中，STDV为无冰稳定地形区域运动速度的标准差，Neff代表无冰稳定地形区域去除自相关影响的有效像元的个数，其计算公式为[24]：

式（3）中，Ntotal代表无冰区稳定地形区域的总像元个数，PS表示为像元分辨率。D为去除自相关影响的空间自相关距离，一般是像元分辨率的20 倍[25]。误差统计结果（表2）表明，误差在可接受范围内。Yan 等[26]基于ALOS/PALSAR 提取的科克萨依冰川2009 年夏季和冬季运动速度分别为0.09 md-1和0.15 md-1，Yang 等[27]基于1998–2010 年Landsat 数据提取的科克萨依冰川上游运动速度为0.14 md-1，Zhou 等[28]认为科克萨依冰川2009 年冬季运动速度最大为0.11 md-1，Holzer 等[4]则认为该冰川2011 年最大运动速度为0.15 md-1，本研究结果（图5）也与上述结果接近。

表2 冰川运动速度误差统计表

4 数据价值

冰川运动与冰川物质平衡、冰川厚度、温度状态、冰川水力特征等密切相关。在气候变化的影响下，冰川运动速度控制冰量输送变化，从而引起冰川长度、面积、厚度变化的响应。冰川运动速度数据不仅是冰川研究的基础数据，更是冰川动力学模型不可或缺的参数数据，在研究冰川变化及预测、冰川资源利用及灾害防治措施中都能发挥重要作用。

5 数据使用方法和建议

1989–2020 年基于Landsat 的东帕米尔高原冰川运动数据集所有数据存储格式均为栅格数据（GeoTIFF 格式），空间坐标系为WGS 1984 UTM 43N。ArcGIS、SuperMap、ENVI、ERDAS 等常用GIS 与遥感软件均支持本数据的读取和操作。基于Landsat 5 的所有数据分辨率为240 m，基于Landsat 7 和Landsat 8 的所有数据分辨率为480 m，信噪比（SNR）无单位，东西向位移（EW）单位为m，南北向位移（NS）单位为m，运动速度（整幅图像的V 或者冰川边界内的glacV）单位为md-1。

致 谢

Landsat 数据来源于USGS（https://earthexplorer.usgs.gov/），中国第二次冰川编目数据集来源于“国家青藏高原科学数据中心”（http://data.tpdc.ac.cn），谨致谢忱。