2000–2020 年喀喇昆仑山洪扎河流域冰川高程变化数据集

吴坤鹏，刘时银，朱钰

1.中国科学院空天信息创新研究院可持续发展大数据国际研究中心，北京 100094

2.云南大学国际河流与生态安全研究院，昆明 650091

3.中国科学院西北生态环境资源研究院国家冰川冻土沙漠科学数据中心，兰州 730000

引 言

作为冰冻圈系统的重要组成部分，冰川与区域气候、自然环境演化、海平面变化、地球表面过程等有着密切联系[1-6]。20 世纪90 年代以来，亚洲高山区冰川不断退缩，冰川持续消融[7-8]，但空间差异显著，主要表现为喜马拉雅地区东部和念青唐古拉山脉的冰川强烈消融，喀喇昆仑山–帕米尔地区冰川呈稳定平衡状态，甚至物质有所增加，称之为“喀喇昆仑异常”[9-10]。喀喇昆仑山冰川变化关系到下游地区的工农业、水力发电等的正常运作，其中西喀喇昆仑山是中国政府推动的“一带一路”陆上丝绸之路关键区之一——中巴经济走廊的必经之地，冰川表面高程变化关系到该地区水资源利用、基础设施建设、灾害防治预警，因此对该地区的冰川变化研究具有重要的社会意义。

当前对冰川变化研究运用最为广泛的方法有传统冰川学法、遥感监测法、模型模拟法等。传统冰川学法通过采集足够数量的花杆和雪坑数据估算冰川物质平衡，是监测单条冰川运用最为广泛的方法；但冰川多发育在边远山区，自然环境恶劣，能够开展野外实地监测的冰川有限[11]。模型模拟法可以模拟年际时间尺度的冰川物质平衡，并预测未来气候变化情景下冰川变化趋势[12-13]，但该方法需要详细的冰川和气象观测资料驱动。遥感监测法通过比较不同时期遥感数据获取的冰川规模、表面高程、运动速度等，分析冰川对气候变化的响应[14-16]。相较于传统冰川学法和模型模拟法，遥感能够开展长时间序列、大空间尺度的冰川变化监测。随着遥感技术的发展和海量数据的积累，遥感监测法已成为研究冰川变化的主要方法。

当前对喀喇昆仑山地区冰川变化研究大多数依赖国外遥感卫星数据，如KH9、ASTER、SPOT等，且与冰川高程变化相关的开源数据集主要为基于ASTER 提取的冰川高程变化[8-10]，包括2000–2016 年高亚洲冰川高程变化数据集[10]（空间分辨率30 m）和2000–2020 年全球冰川高程变化数据集[8]（空间分辨率100 m）。为凸显我国遥感数据在科学研究中的实际价值，本文选用2020 年1–2 月获取的资源3 号卫星立体像对数据，基于大地测量法[17]，通过DEM 误差校正、冰雪穿透估算等[18]，获取了洪扎河流域冰川表面高程变化，并以稳定地形区残余高程差为统计特征，判别和验证数据集产品的可靠性，为“喀喇昆仑异常”研究及该地区水资源利用、冰川灾害防治等提供更精准的数据支撑。

1 数据采集和处理方法

1.1 数据源

2000–2020 年喀喇昆仑山洪扎河流域冰川高程变化数据集主要基于三种数据：2000 年2 月的数字高程模型（SRTM DEM）、2020 年1–2 月获取的资源三号卫星立体像对数据和冰川编目数据（表1、图1）。

表1 资源三号立体像对参数信息Table 1 Specifications of ZY3 stereo images

SRTM DEM 为美国地质调查局提供的、未填补“空洞”的 SRTM1 C 波段数据产品（https://earthexplorer.usgs.gov/），获取时间为2000年2月，空间分辨率30 m，坐标系WGS84/EGM96。在90%的置信区间内，SRTM 垂直精度优于16 m，其高程精度受地形的影响较大，在平坦区域的高程精度可达10 m，且随地表坡度变化[19]。SRTM 有两套合成孔径雷达系统，分别采用C 波段（波长5.6 cm，频率2.2 Hz）和X 波段（波长3.1 cm，频率8.8 Hz），为有效评估C 波段对积雪的穿透，本文还选取了覆盖喀喇昆仑山地区的SRTM X 波段数据，通过C 波段和X 波段差值比较，估算C 波段对积雪的穿透深度[20]。

资源三号是我国首颗民用三线阵立体测绘卫星，于2012 年1 月9 日成功发射。资源三号卫星搭载1 台多光谱和3 台全色线阵光学相机。全色相机包括前视、正视和后视组成三线阵立体相机，其中前视和后视地面分辨率为3.5 m[21-22]。

冰川编目数据由国际冰雪数据中心（https://nsidc.org/data/glims）的RGI 6.0（Randolph Glacier Inventory 6.0）[23]，主要用于冰川区与非冰川区的划分。由于喀喇昆仑山洪扎河流域存在大量跃动前进冰川，为了确保冰川表面高程变化的完整性，本文还采用了RGI 6.0[23]、中国第二次冰川编目[24]及GAMDAM（Glacier Area Mapping for Discharge from the Asian Mountains）[25]冰川编目的集合，以保证在研究时段冰川边界的最大范围。

1.2 数据处理方法和流程

冰川表面高程变化获取主要由2 部分组成：基于立体像对提取2020 年数字高程模型，基于两期数字高程模型获取2000–2020 年的冰川表面高程变化。本文采用ENVI、ArcGIS 和Python 软件平台进行数据的处理和分析，其中ENVI 用于DEM 提取，ArcGIS 用于DEM 编辑与镶嵌，Python 用于DEM 配准，主要处理流程如图2。

图2 ZY3 DEM 的提取及冰川表面高程变化生成流程示意图Figure 2 The flow chart of ZY3 DEM extraction and glacier height changes

基于资源三号立体像对提取数字高程模型总体分为6 步，分别是输入立体像对、定义地面控制点、定义连接点、设定DEM 提取参数、输出并编辑DEM、镶嵌DEM[26]。（1）输入立体像对：左影像加载正视影像，右影像加载后视/前视影像，加载影像后，ENVI 软件能够自动识别RPC 文件；（2）定义地面控制点：ENVI 软件提供不定义、交互式定义和读取控制点文件3 种定义地面控制点的方式，其中不定义控制点提取的是相对高程，水平面为WGS84 基准面。由于获取冰川表面高程变化需要进行两期DEM 配准，因此提取DEM 可采用不定义控制点方式；（3）定义连接点：ENVI 软件提供自动寻找、交互式定义和读取连接点文件3 种定义连接点的方式，本文先选择自动寻找连接点，然后进行交互式编辑，确保连接点误差达到最小；（4）设定DEM 提取参数：设定投影参数为WGS84 UTM Zone 43N，输出像元大小为10 m，背景值为0；（5）输出并编辑DEM：生成DEM 并检查由云量等造成的DEM 错误，通过手动编辑修改高程数据；（6）镶嵌DEM：利用ArcGIS 软件将提取的DEM 进行镶嵌，得到研究区相对DEM（以下称ZY3 DEM）。

两期DEM 数据（SRTM DEM 和ZY3 DEM）之间的配准采用大地测量法。同一区域不同时间的两幅DEM 数据，在没有配准的情况下直接高程相减得到高程差值图，而高程差值与坡度坡向之间存在如下余弦函数关系：

式中：dh为不同DEM 数据之间各个像元的高程差，a是平移矢量的大小，b是平移矢量的方向， 和φ分别为像元对应的坡度和坡向， ̅ℎ̅̅̅为不同DEM 数据之间的整体高程差异（垂直矢量的大小）。从两幅DEM 对应的像元属性中可以获得dh、坡度 和坡向φ。不同DEM 数据之间的高程差存在异常值，影响余弦曲线的拟合。本文对高程差样本进行统计，认为低于5%或高于95%分位数的高程差是异常值，进行剔除。坡度的大小对DEM 数据匹配的精度有很大影响，因此要剔除地形平坦的区域和坡度太大的区域。本文剔除了地面坡度小于5°和大于80°的区域。然后通过最小二乘法拟合余弦曲线可以得到系数a、b、c，不同DEM 数据在X、Y 和Z 方向的平移量可表示为：

式中：α̅为DEM 数据的平均坡度。求出shiftx，shifty和shiftz，并按照在X、Y 和Z 方向的平移量进行平移，但一次平移后的DEM 并不一定能够完全消除空间匹配错位造成的误差，需要将平移后的DEM 再次迭代平移，重复上述步骤计算dh标准差和平移矢量。当dh的标准差减小幅度小于2%或平移距离小于0.5 m，即达到了两幅DEM 水平方向的配准要求[17]。最后利用非冰川区高程差与最大曲率的关系，校正由不同分辨率造成的冰川区的高程差异残差。最终得到的高程差值即为研究区高程变化信息。

2 数据样本描述

本数据集的数据存储格式为 GeoTIFF 格式，数据类型为 32 位浮点型，命名为DH_Karakoram_Hunza_20002020.tif，其中DH 表示为Difference Height，Karakoram_Hunza 表示为喀喇昆仑山洪扎河流域，20002020 表示为2000 年至2020 年时间段，数据单位是米，地理坐标系WGS84 UTM Zone 43N，空间分辨率30 m。样本展示如图3。

图3 喀喇昆仑山洪扎河流域冰川表面高程变化Figure 3 The glacier height changes in Hunza Basin

3 数据质量控制和评估

由于SRTM 的C 波段能够穿透积雪，在测定积雪覆盖的冰川表面地形时得不到真实的地面信息。穿透深度随着地面状况的变化而不同，以新雪和粒雪的穿透深度最大，表碛物覆盖的冰面穿透深度最小[27]，穿透深度能介于0–10 m 之间[20,28-29]。相比SRTM 的C 波段，X 波段对积雪的穿透效应要小[20]。利用配准和校正后的SRTM X 波段DEM 和SRTM C 波段DEM 相比较，得到SRTM C波段在喀喇昆仑山地区平均穿透1.6±0.05 m。SRTM DEM 和ZY3 DEM 在高海拔地区均存在一定“空洞”，但通过大地测量法获取的冰川表面高程变化数据，在不同海拔均有变化信息，数据空洞占比较小，对冰川表面高程变化统计分析不会产生显著影响。

基于SRTM DEM 和ZY3 DEM 获取地面高程变化后，可利用稳定地形区高程差的均方根误差（Root Mean Square Error,RMSE）或标准差（Standard Deviation,STDV）作为高程差不确定性的初步估计（图4），但是会高估高程变化的不确定性，因为误差在一个较大区域内平均后会减小[30-31]。根据空间邻近相似性原理，相邻的像元间存在自相关，因此需要除去自相关造成的影响。为此，引入标准平均误差（Standard Error of the mean,SE）对高程差不确定性进行评估。

图4 无冰区高程差及其直方图Figure 4 Elevation changes and the histogram in off-glacier areas

式中：n为去自相关后的像元个数。除去自相关造成的影响，需要根据DEM 的分辨率计算空间自相关的长度。空间自相关距离可以通过自定义的方式确定，也可以通过莫兰指数（Moran’s I）计算。研究表明，30 m 分辨率选取600 m 的自相关长度，10–20 m 分辨率选取400 m 的自相关长度较为合适[30]。因此，本文采用自定义方式，自相关长度选择400 m。最终SRTM DEM 和ZY3 DEM 数据高程差的精度（σ）用无冰区高程差的均值MED和标准平均误差SE计算：

结果显示，误差校正后，SRTM DEM 和ZY3 DEM 数据高程差的误差减小到±1.1 m。

4 数据价值

相较于开源的冰川表面高程变化数据集，本数据集比2000–2016 年高亚洲冰川高程变化数据集[10]有更长的时间跨度、比2000–2020 年全球冰川高程变化数据集[8]具有更高的空间分辨率。更长的时间跨度与更高空间分辨率，能够为喀喇昆仑山地区冰川变化研究提供更充实的基础数据。

“喀喇昆仑异常”是当前研究热点，主要表现为轻微的物质正平衡及大量跃动冰川的分布[15]。开源跃动冰川分布数据集多是根据冰川末端进退及冰川运动速度变化进行制备[32-33]，对冰川跃动过程中冰量的输送研究较为欠缺，难以充分开展冰川跃动过程及其机理研究。本文基于SRTM DEM 和ZY3 DEM 获取了喀喇昆仑山洪扎河流域21 世纪初冰川表面高程变化，为更完整的跃动冰川识别及典型冰川的跃动控制机理研究奠定了数据基础。与此同时，洪扎河流域是中巴经济走廊的必经之地，全面了解冰川表面高程变化、估算冰川物质平衡，为该地区水资源利用及冰川灾害防治提供必要的基础信息。

5 数据使用方法和建议

喀喇昆仑山洪扎河流域2000–2020 年冰川表面高程变化数据集通过国家冰川冻土沙漠数据中心NCDC（National Cryosphere Desert Data Center）网站提供全部数据下载服务。本数据集数据存储格式为32 位浮点型GeoTIFF 栅格数据，地理坐标系WGS84 UTM Zone 43N。可以在ArcGIS、SuperMap等GIS 平台软件直接调用，也可以使用MATLAB 等数据分析软件编码进行批量处理和分析。研究人员可利用该数据与其他时期冰川物质平衡相比较，或选择该数据集作为冰川模型的输入参数/验证参数，开展冰川对气候变化响应研究工作。

致 谢

感谢美国地质调查局（USGS）提供SRTM1 DEM 数据，感谢国际冰雪数据中心及国家冰川冻土沙漠科学数据中心提供冰川编目数据集。

数据作者分工职责

吴坤鹏（1990—），男，安徽省安庆市人，博士，副研究员，研究方向为冰川遥感。主要承担工作：数据筛选，数据处理，论文撰写。

刘时银（1963—），男，河南省信阳市人，博士，研究员，研究方向为冰川水文。主要承担工作：研究思路设计。

朱钰（1992—），男，甘肃省平凉市人，博士，助理研究员，研究方向为冰川遥感。主要承担工作：数据处理。