2015年东帕米尔高原克拉牙依拉克冰川跃动数据集 中华文库

摘要：2015年5月东帕米尔高原克拉牙依拉克冰川（38°35′6″N–38°44′48″N，75°7′47″E–75°22′29″E）发生跃动，给当地牧民生产、生活带来一定的损失。加强对跃动冰川的监测和分析，对冰川跃动机理和灾害预警预报研究具有重要意义。本研究收集了冰川跃动期间Landsat 8影像，基于特征匹配法提取了克拉牙依拉克冰川2015年4月13日至7月11日表面运动速度数据集（GeoTIFF格式，32位浮点型），空间分辨率为240 m，冰川运动速度误差最高为±0.42 m d-1。根据跃动前后的2013年和2015年两期质量较好的ASTER立体像对数据生成的DEM，提取了克拉牙依拉克冰川高程变化数据（GeoTIFF格式，32位浮点型），空间分辨率为30 m，经非冰川区域的高程残差统计，冰川高程变化误差为±0.7 m。本数据集可以作为跃动冰川的本底调查资料，也可以作为冰川动力学模型的输入，为冰川跃动研究提供基础数据支撑。另外，跃动冰川的监测结果有助于冰川灾害预警预报研究，为中巴经济走廊的经济社会发展提供数据支撑。

关键词：帕米尔；跃动冰川；冰川表面运动；冰川高程变化；Landsat；ASTER；克拉牙依拉克冰川

Abstract: Karayaylak Glacier (38°35′6″N–38°44′48″N, 75°7′47″E–75°22′29″E) is located in the northern slope of Jiubie Peak of Kongur Mountains, eastern Pamir, China. Its surge in May 2015 brought property losses to local residents. To enhance the monitoring and analysis of the surged glacier has important practical significance for exploring the glacier surge mechanism and for establishing a disaster early-warning mechanism. This study obtained the surface velocities of Karayaylak Glacier from April 13, 2015 to July 11, 2015 from Landsat images by using feature-tracking method. With a spatial resolution of 240 m, these data were stored in GeoTiff format in 32-bit floating point. A statistical analysis indicates that the results have an overall accuracy of less than ± 0.42 m d-1. The elevation change between 2013 and 2015 was also obtained from ASTER-derived DEM. With a spatial resolution of 240 m, the data were stored in GeoTiff format in 32-bit floating point. A statistical analysis of the residual error in non-glacial areas indicates that the results have an overall accuracy of about ± 0.7 m. The dataset can serve as the basis for investigating surged glaciers in this region, or as input data for building glacier dynamic models. It provides statistical support for glacier surge research. In addition, the monitoring data not only support further research in risk assessment and the early warning of glacier disasters, but also ensure the economic and social development of China–Pakistan Economic Corridor.

Keywords: Pamir, glacier surge; glacier surface velocity; glacier elevation change; Landsat; ASTER; Karayaylak Glacier

跃动冰川是具有间歇性和周期性快速运动的冰川，其在跃动期间运动速度突然加快并引起冰川物质快速转移和冰体重新分布。冰川跃动时不仅可能淹没下游森林、道路、村庄、草场等，而且可能带来河道阻塞、冰湖溃决、泥石流、冰崩等灾害，严重威胁生命和财产安全^[1]。世界主要冰川作用区均有跃动冰川的报道，帕米尔和喀喇昆仑山地区冰川跃动较为频繁^[2][3]。2015年5月，东帕米尔高原公格尔山九别峰北坡克拉牙依拉克冰川西支（图1）突然发生跃动，造成阿克陶县61户牧民房屋受损，上百头牲畜失踪，部分草场被淹没。该跃动事件不仅受到新闻媒体和科研工作者的广泛关注，而且发生在中巴经济走廊穿越区、丝绸之路经济带上重要交通区域，赋予了本研究重要的现实意义。为认识常态冰川和跃动冰川的控制机理和准确识别跃动冰川过程，需要针对已发生跃动的冰川建立动力学模型。在建模的过程中，冰川运动速度和冰川高程变化数据是不可或缺的。由于冰川所处环境恶劣，野外观测困难，并且缺乏定位观测，因此遥感成为监测跃动冰川的主要手段^[4][5]。基于遥感影像的冰川表面运动速度是指两期或者多期影像的同名特征点在单位时间内的位移变化。基于遥感方法获取冰川运动速度和冰川高程变化均在冰川研究中取得了广泛应用^[6][7]。

克拉牙依拉克冰川位于东帕米尔高原公格尔山九别峰（海拔7530 m）北坡，总长度约20.3 km。根据我国第二次冰川编目数据^[8]，该冰川面积115.6 km2，是东帕米尔高原最大的冰川。该冰川覆盖大量表碛，表碛覆盖率约为22.3%，冰川雪线高度4220 m，末端海拔2817 m。本文利用Landsat 8影像提取了克拉牙依拉克冰川跃动期间冰川表面平均运动速度数据，并且利用ASTER立体像对提取了跃动前后的高程变化，采用方法均已得到大量验证^[6][7]，获取数据可靠，可为进一步研究跃动冰川提供必要支撑。

图1 克拉牙依拉克冰川地理位置示意图

1.1 数据源与数字高程模型（DEM）的提取

本研究收集了跃动前后（2015年4月13日至2015年7月11日）Landsat 8卫星影像（表1），用于提取冰川表面运动速度。采用2013年6月11日和2015年6月17日获取的ASTER L1A立体像对数据，分别在ENVI中提取数字高程模型（DEM）。平面控制坐标来源于Landsat影像，高程坐标来源于地形图DEM，坐标系设置为WGS1984 UTM 43N，空间分辨率设为30 m。

表1 研究采用的数据

1.2 数据处理方法和流程

1.2.1 冰川运动速度数据处理方法和流程

冰川表面运动速度是通过不同时相影像之间冰川表面特征的位移量计算获得的。对于光学影像来说，其位移量可以通过加州理工学院开发的COSI-CORR（Co-registration of Optically Sensed Images and Correlation）软件包获取。COSI-CORR基于IDL开发并集成于ENVI软件中，在ENVI 5.0版本中支持ENVI classic（ENVI经典界面）。本研究采用2014年10月更新的版本（下载地址：http://www.tectonics.caltech.edu/slip_history/spot_coseis/download_software.html）。COSI-CORR根据光学遥感图像的特征匹配及互相关计算位移信息，其原理是基于频域互相关算法求影像的偏移量^[9]。所有的相位相关依赖于傅立叶变换理论：一对影像的相对偏移量通过傅里叶变换的相位差进行提取。用i1 和i2 代表两景影像，x和y为坐标值，△x 和△y 代表偏移量，那么：

\({i}_{2}\left(x,y\right)={i}_{1}\left(x-{△}_{x},y-{△}_{y}\right)\) (1)

用I1 和I2 分别代表通过傅立叶移位定理得到的傅里叶变换，可以表示为：

\({I}_{2}\left({\mathrm{\omega }}_{x},{\omega }_{y}\right)={I}_{1}\left({\omega }_{x},{\omega }_{y}\right){e}^{-j\left({\omega }_{x}{△}_{x}+{\omega }_{y}{△}_{y}\right)}\) (2)

ωx 和ωy 分别代表行和列上的频率变量，e为自然常数，j表示虚数形式。因此图像i1 和i2 标准互功率谱（C）可以表示为：

\({C}_{{i}_{1}{i}_{2} }\left({\omega }_{x},{\omega }_{y}\right)=\frac{{I}_{1}\left({\omega }_{x},{\omega }_{y}\right){I}_{2}^{*}\left({\omega }_{x},{\omega }_{y}\right)}{\left|{I}_{1}\left({\omega }_{x},{\omega }_{y}\right){I}_{2}^{*}\left({\omega }_{x},{\omega }_{y}\right)\right|}={e}^{j\left({\omega }_{x}{△}_{x}+{\omega }_{y}{△}_{y}\right)} \) (3)

其中*表示复共轭。通过将标准互功率谱进行傅里叶逆变换（\({\mathcal{F}}^{-1}\)），便可得到狄拉克δ（Dirac delta）函数\(\delta \left(x+{△}_{x},y+{△}_{y}\right)\)：

\({\mathcal{F}}^{-1}\left\{{e}^{j\left({\omega }_{x}{△}_{x}+{\omega }_{y}{△}_{y}\right)}\right\}=\delta \left(x+{△}_{x},y+{△}_{y}\right)\) (4)

通过找到\(\delta \)函数的峰值位置来确定△x 和△y 。当偏移量为亚像元级别时，这个峰值就不是狄拉克δ（Dirac delta）函数了，而是下采样的狄利克雷（Dirichlet kernel）函数。因此COSI-CORR的相位相关方法可分为两个步骤：（1）通过计算图像的互功率谱的线性相位重构图像的相对偏移；（2）通过严格的相关峰值的位置确定图像的相对偏移量。在计算相对偏移量前需要对其中的一景影像做正射校正处理，然后以正射影像作为基准影像对另外一景影像进行配准处理。本研究采用的Landsat数据均从美国地质调查局（USGS）获取，所有数据均由USGS采用的GLS2005全球陆地控制点系统进行了正射校正。郭万钦等^[10]对其校正精度进行了验证，结果显示USGS提供的正射校正后的Landsat遥感影像具有相当高的校正精度，大部分影像校正精度都在半个像元左右，部分影像甚至达到1/6～1/10左右的精度。因此，认为Landsat影像的正射校正精度满足表面运动速度分析。本研究采用Landsat OLI第8波段（全色波段，空间分辨率为15 m），在COSI-CORR计算相关系数时采用频率域算法，参考窗口设为128，搜索窗口设为32，获得的位移数据有3个图像组成：东西向位移（EW）、南北向位移（NS）以及信噪比（SNR）。本研究采用SNR≥0.9部分并剔除云和阴影覆盖部分，最后即可得到冰川表面运动场的信息。

1.2.2 高程变化数据处理方法和流程

冰川高程变化由跃动前后两景ASTER DEM求差值即可获取。不同ASTER DEM可能存在空间匹配误差，求差运算前按照图2的流程消除该部分误差。理论上空间匹配造成的高程误差（dh）可采用公式（5–6）进行描述^[7]：

\(\frac{dh}{\mathit{tan}\left(\alpha \right)}=a*\mathit{cos}\left(b-\phi \right)+c\) （5）

\(c=\frac{\stackrel{-}{dh}}{tan\left(\stackrel{-}{\alpha }\right)}\) （6）

其中α为坡度，\(\phi \)为坡向，\(\stackrel{-}{dh}\)为不同DEM整体高程差异，代表垂直偏移量，\(\stackrel{-}{\mathrm{\alpha }}\)为基准DEM提取的平均坡度。参数a、b、c可以通过回归分析获取。不同DEM数据间X、Y、Z方向上的偏移校正量可采用公式（7–9）获取：

\(X=a*sin\left(b\right)\) （7）

\(Y=a*cos\left(b\right)\) （8）

\(Z=c*tan\left(\stackrel{-}{\alpha }\right)\) （9）

选择无冰区稳定地形对如上述过程进行迭代，当dh的标准差减小幅度小于2％或者（X2+Y2 ）小于0.25即可完成迭代。在校正的过程中，选择高程偏差的5％和95％分位数以剔除异常值的影响^[11]。

由于受积累区云雪影响，数据结果不理想，并且跃动发生在消融区，因此本研究只提取了消融区部分的数据。本数据集的数据均存储为GeoTIFF格式，数据类型为32位浮点型，分为两个文件夹存储，分别为冰川运动速度（Velocity）和冰川高程变化（ElevationChange）。前者文件夹中含有5个文件，命名规则为YYYYMMDD_YYYYMMDD_velocity.tif，表示两个时间间隔冰川表面运动速度，单位m d-1，结果显示如图3。后者文件夹中只有一个文件（ElevationChange.tif），该文件为2013–2015年冰川表面高程变化，单位为m，数据结果显示如图4（图3和图4均是通过ArcGIS渲染得到的）。

图2 DEM空间匹配校正主要流程示意图[9]

图3 冰川运动速度（底图为SRTM生成的山体阴影，冰川边界来源于我国第二次冰川编目数据集[13]）

图4 冰川高程变化（底图为SRTM生成的山体阴影，冰川边界来源于我国第二次冰川编目数据集[13]）

表面运动速度的误差可用公式（10）和（11）评估：

\(E_{d}=\sqrt {E_{1}^{2}+E_{2}^{2}}\) (10)

\(E_{v}=E_{d}/d\) (11)

式中Ed 和Ev 分别为两时间间隔的位移与速度的误差；E1 和E2 分别是东西向和南北向位移误差，d为两日期相隔的天数。本研究对遥感结果中无冰区稳定地形区域东西向和南北向位移偏移进行统计（表2），发现无冰区稳定地形的东西向和南北向偏移平均值均接近于0。本文假定无冰区稳定地形东西向和南北向理论上不发生偏移，其实际的偏移量可作为位移变化残差，并认为冰川区域与无冰区残差相似。因此本文将无冰区东西向和南北向偏移的标准差作为E1 和E2 ，估计冰川区表面运动速度的误差（表2）。对于SNR＜0.9的部分，用空值表示。

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

配准后的两期DEM差值结果在无冰区是趋向于0的，本文假定理论上无冰区高程变化值为0，其实际结果可认为是DEM高程差的残差。高程变化的误差（σ）可以根据经空间匹配校正后两期DEM的无冰区平均高程差（MED）和标准差（STDV无冰区 ）进行评估^[12]：

(12)

(13)

其中SE为引入的中间变量即标准平均误差（Standard error of the mean），N为包含的像元数，为了消除空间自相关影响，按600 m进行采样，N值为采样的像元数。最终结果显示，无冰区残差平均值（MED）减少到0.39 m，高程变化的误差（σ）为±0.70 m。受云的影像，部分像元值为空值。

冰川运动速度遥感监测主要分为基于光学影像或者SAR图像的特征匹配法和雷达干涉测量法。受卫星重访周期的限制，冰川区获取高相干图像较为困难。特征匹配法在山地冰川上多有研究，被证明了可靠性，其所得结果与野外花杆测量一致^[13]。冰川高程变化遥感监测最为成功的方法是本研究使用的大地测量法。大地测量法所需要的DEM主要来源于光学立体像对和雷达干涉测量的方法。雷达干涉测量限制于相干性和对冰雪可能的穿透，光学立体像对需要视觉对比度较优，受积雪等因素限制。本研究只考虑消融区部分，视觉对比度较优，所得结果可靠。因此，本研究结果在一定程度上为缺乏野外观测的跃动冰川研究提供数据支撑。

冰川表面运动速度数据和高程变化数据均是研究冰川物质平衡、冰川动力学模型、冰川灾害预警预测的基础数据。本研究利用Landsat 8影像提取了冰川跃动期间的运动速度数据，反映了冰川跃动前后的运动特征；并利用ASTER立体像对数据提取了跃动前后的高程变化数据，反映了冰体重新分布的特征。二者为跃动冰川动力学模型的发展提供基础数据。

2015年东帕米尔高原跃动冰川数据集所有数据存储格式均为GeoTIFF格式，空间坐标系为WGS1984 UTM 43N。ArcGIS、SuperMap、ENVI、ERDAS等常用的GIS与遥感软件均支持该数据的读取和操作。冰川运动速度日变化值以影像的像元值表示，单位m d-1，空间分辨率为240 m；冰川高程年变化值以影像的像元值表示，单位为m，空间分辨率为30 m。本研究数据集的空洞建议采用空间插值的方法进行处理，可将数据集应用于冰川动力学模型，研究冰川跃动规律。

感谢USGS提供的Landsat、SRTM 1数据，感谢NASA EARTHDATA提供ASTER数据，感谢寒区旱区科学数据中心提供中国第二次冰川编目数据集。

张震, 刘时银, 魏俊锋. 2015年东帕米尔高原克拉牙依拉克冰川跃动数据集[DB/OL]. Science Data Bank, 2018. (2018-10-31). DOI: 10.11922/sciencedb.572.