帕米尔中部North Kyzkurgan冰川跃动变化遥感监测

李志杰， 王宁练，3， 侯姗姗

（1.陕西省地表系统与环境承载力重点实验室，陕西西安 710127； 2.西北大学城市与环境学院地表系统与灾害研究院，陕西西安 710127； 3.中国科学院青藏高原地球科学卓越创新中心，北京 100101）

0 引言

跃动冰川指周期性发生快速运动的冰川，当冰川发生跃动时其运动速度可达到平时的数倍乃至上百倍，进而引起冰川物质的快速转移和重新分布［1-2］。尽管跃动冰川数量仅占全球冰川的约1%，但其重要影响却不可忽视［3］。一方面，探究冰川发生跃动的原因机制有助于拓展当前对于冰川变化的认识；另一方面，监测跃动冰川可为防灾减灾提供重要支撑［4］。冰川跃动可能直接导致冰湖溃决洪水、冰川融水泥石流、高山堰塞湖等一系列自然灾害，对下游地区产生重大威胁［5］。例如，2015 年5 月新疆公格尔九别峰北坡的克拉牙依拉克冰川发生跃动，造成部分草场和房屋被冰体淹没冲毁［6］。因此，关注和监测跃动冰川具有十分重要的科学和现实意义。

近年来“帕米尔-喀喇昆仑异常”现象受到了广泛关注，多项研究表明自1970s 以来，帕米尔、喀喇昆仑、西昆仑等地区的冰川接近或处于正平衡状态，基本保持稳定甚至有扩张现象［7-9］。跃动冰川的大量存在是“帕米尔-喀喇昆仑异常”现象的突出特征，由此蕴含的冰川灾害风险也在不断增大［4］。帕米尔跃动冰川主要分布在高原西北部和东部，Goerlich 等［10］利用多源DEM 和Landsat 影像对帕米尔跃动冰川进行了编目，认为该地区有多达186 条冰川发生过跃动现象，其中也包括本研究中的North Kyzkurgan冰川，但Goerlich等未深入分析这些冰川的跃动过程，其编目结果也存在一定的不确定性。Kotlyakov 等［11］利用多源影像和实地调查资料对西帕米尔冰川进行监测，发现1960—2003 年有19 条冰川发生跃动。Lü 等［12］利用Landsat ETM+/OLI、ASTER 影像对东帕米尔昆盖山冰川进行监测，发现1999—2016 年有9 条冰川发生跃动。此外，Wendt等［13］、张震等［14］还分别报道了Fedchenko 冰川的Bivachny 支冰川、昆盖山5Y663L0023 冰川的跃动现象，深入分析了冰川跃动的过程与机理。关于2015年5 月的东帕米尔克拉牙依拉克冰川跃动事件，Shangguan 等［15］、张震等［6］、冯力力等［16］分别进行了报道。因此在总体上，当前对于帕米尔跃动冰川的认识多集中在局部区域或特定冰川，而对多数跃动冰川缺乏深入研究，对其跃动时间、过程、周期以及最新的动态变化等尚缺乏基本了解。

跃动冰川具有明显的表面特征，如末端在数月或数年内突然前进、冰面裂隙发育和破碎化、积蓄区和接收区高程的剧烈波动、冰流速突然加快至平时的数倍到数十倍等，通过这些特征可有效的识别和监测跃动冰川［1，17］。遥感和GIS技术的快速发展，为山地冰川的监测研究提供了有效的技术手段。通过多源遥感影像和DEM 等对山地冰川开展联合监测，可有效识别出冰川跃动导致的面积、高程、流速变化［18］。本文以帕米尔中部近50 年来面积、高程、流速变化极为显著的North Kyzkurgan 冰川为研究对象，利用1973 年以来的Landsat 影像、ASTER立体像对和ITS_LIVE 数据产品，揭示了该冰川的完整跃动发生过程，并对跃动周期、跃动原因等进行了探讨，以深化对于帕米尔跃动冰川的认识。

1 研究区概况

帕米尔高原指位于高亚洲西部，阿赖山以南，兴都库什山以北的广阔高原山地，面积超过12×104km2（图1）。帕米尔高原地形高峻，西部发育了一系列呈西南—东北走向的平行山脉和山间河谷，使得湿润的西风气流可以顺坡爬升，带来丰富的降水，因此发育了面积广阔的山地冰川［19］。North Kyzkurgan 冰川（G072447E38852N）位于帕米尔高原中部，是一条典型的复式山谷冰川，末端有少量表碛覆盖。 2020 年时冰川面积为（26.27±0.50）km2，最大长度约9.4 km，平均高程4 960 m，平均坡度约16°。若以雪线高度（4 700 m）划分［20］，该冰川的积累区面积比率超过0.8，冰川作用正差约1 000 m。Fedchenko气象站位于研究冰川以西约20 km，该站的观测资料表明1935—1990 年该地区年均气温约为-6.9 ℃，年均降水量约为1 200 mm且以降雪为主，极低的气温和丰富的降雪为该冰川提供了充足的物质补给。

图1 North Kyzkurgan冰川的地理位置Fig. 1 Geographic location of the North Kyzkurgan Glacier

2 数据与方法

2.1 数据来源

2.1.1 遥感影像

本研究使用了1973—2020 年的20 景Landsat MSSTMETM+OLI 影像，这些影像均为L1T 级产品，由USGS（United States Geological Survey）进行了系统辐射校正和几何校正，并结合DEM 进行了地形校正（表1）。由于Landsat 系列影像有相当高的正射校正精度，达到1/2 像元左右，因此本研究直接利用Landsat 影像开展冰川变化研究［17，21］。本研究还选用了获取于2002 年、2009 年、2014 年、2020年的4 对ASTER 立体像对，以提取不同时期的DEM，探究冰川跃动前、跃动中和跃动后的表面高程变化（表1）。为避免季节性积雪、云、山地阴影的影响，本研究尽可能选取了消融季末期，冰川区无云、阴影覆盖的影像。Landsat 影像来源于USGS（http：//glovis. usgs. gov），ASTER 立体像对来源于NASA（https：//search.earthdata.nasa.gov/）。

表1 遥感影像数据列表Table 1 List of the satellite data used in this study

2.1.2 冰川流速数据

ITS_LIVE（Inter-mission time series of land ice velocity and elevation）产品包含了1985—2018 年高亚洲主要冰川的流速数据，来源于NASA 的MEa-SUREs 项目（https：//its-live. jpl. nasa. gov/）。该产品基于Landsat 4、5（1～4 波段）和7、8（全色波段）影像，采用Gardner 等［22］提出的自动裂缝特征跟踪处理链方法，通过进行局部归一化、过采样、特征跟踪等提取冰川运动速度。该数据具有良好的鲁棒性，减少了因图像匹配、冰川跃动等因素引起的误差，由于Landsat 8 空间分辨率和波谱分辨率的提高，2013 年以后该产品的误差显著降低，数据质量得到进一步提升［23-24］。ITS_LIVE 产品已广泛应用于高亚洲冰川运动的相关研究中，充分验证了其可靠性和准确性［24-25］。本研究中选用了2000—2018 年的单年流速及误差数据，空间分辨率为240 m，数据版本为最新的V01。

2.1.3 气象数据

本研究使用了塔吉克斯坦Fedchenko 气象站1935—1990 年的气温和降水观测数据（https：//nsidc. org/data/G02174），该气象站位于Fedchenko冰川中游海拔4 179 m 处（38.83° N，72.22° E），距North Kyzkurgan冰川中心仅约20 km（图1）。

2.2 方法

2.2.1 冰川边界解译

本研究采用目视解译方法提取冰川边界。由于North Kyzkurgan 冰川末端仅存在少量表碛，因此本研究直接在Landsat假彩色合成影像中，根据冰川区与非冰川区显著的色彩差异进行冰川边界的数字化［26-27］。在实际解译过程中，本研究首先完成2020 年的冰川边界提取，在此基础上结合其他年份的Landsat 影像，保持冰川上部积累区不变，只修改冰川末端发生变化的部分，从而获取其他年份的冰川边界。

2.2.2 DEM空间匹配

本研究基于ENVI 软件中的“DEM Extraction”模块利用ASTER 立体像对提取DEM。我们采用相对定向法提取DEM，在每对立体像对中自动生成连接点并进行修正，确保每对影像中连接点不少于80个且均匀分布，最终连接点的Y方向最大视差控制在一个像元以内［28］。为保证DEM 提取的精度，将提取的ASTER DEM 空间分辨率设置为30 m 并全部投影到WGS84 UTM zone 43 N坐标系下。

我们采用统计学方法对DEM 空间匹配误差进行校正，但多源DEM 之间的高程差是数据匹配偏差与冰川变化共同作用的结果，因此在非冰川区开展校正工作以排除冰川变化的影响。本研究基于DEM 的空间匹配偏差与坡向、坡度等地形因子的相关关系，对不同时期ASTER DEM 之间的空间匹配偏差进行了计算和校正（表2）［29］。在校正过程中我们采用±100 m 作为剔除非冰川区高程差中异常值的阈值，还剔除了地面坡度小于5°的区域，以提高匹配精度［30］。

表2 不同DEM之间的空间偏移量Table 2 The X-Y-Z shift vector between two DEMs

2.3 不确定性分析

2.3.1 冰川边界提取的不确定性

本研究采用统计冰川轮廓线经过的像元数量的方法，评价目视解译提取冰川边界的不确定性［31］。如公式（1）：

式中：EA为解译结果的不确定性，N为冰川轮廓线经过的像元数量，λ为像元面积（MSS影像为3 600 m2，TM为900 m2，ETM+/OLI为225 m2）。

冰川面积变化的不确定性（EAC）由公式（2）计算得出：

式中：EA1和EA2分别为两个时期冰川边界提取的不确定性。

2.3.2 冰面高程变化的不确定性

经过空间匹配后，尽管DEM 之间在非冰川区的高程差已趋近于0，但仍会存在一定的高程残差，因此必须对DEM 的空间匹配精度进行评价。如公式（3）、（4）所示［30］：

式中：SE 和STDV 为非冰川区DEM 高程差的均方根误差和标准差；n为去空间自相关处理后的像元数量（本研究采用600 m作为去空间自相关距离）；σ为DEM 间的高程变化误差；MED 为非冰川区DEM间的平均高程差。表3 为不同时期DEM 空间匹配校正前后的误差特征分布，可见校正后DEM 间在非冰川区的高程差明显减小，匹配精度显著提高。

表3 DEM高程差不确定性Table 3 The uncertainties of DEM differences

3 结果

3.1 冰川面积和长度变化

基于Landsat 影像的冰川边界解译结果表明（图2），1973—2020 年North Kyzkurgan 冰川的面积和长度都经历了显著变化。1973—2006 年，North Kyzkurgan 冰川处于退缩状态，末端累计退缩约1.31 km，冰川面积也由（26.63±2.00）km2减少为（24.27±0.48）km2。2006—2011 年，North Kyzkurgan 冰川整体处于退缩状态，面积减少至（24.01±0.93）km2，但值得注意的是在此期间，在冰川末端退缩的同时还发生了膨胀变宽的现象［图3（a），黑色箭头处］。总体上，1973—2011 年North Kyzkurgan 冰川处于退缩状态且退缩速率较为稳定，冰川面积累计减少（2.62±0.41）km2（9.84%±1.54%）。2011—2016 年，North Kyzkurgan 冰川末端突然大幅度前进了1.17 km，冰川面积也急剧扩张了（2.49±0.14） km2（10.37%±0.58%），达到了（26.50±0.50）km2。这6 年内North Kyzkurgan 冰川的长度和面积变化极为显著，几乎相当于该冰川自1973—2011 年来近40 年的长度和面积退缩量。2016 年，North Kyzkurgan 冰川面积和长度达到了1970s以来的最大值，此后又开始退缩。

图2 1973—2020年North Kyzkurgan冰川面积和长度变化Fig. 2 Glacier area and length change for North Kyzkurgan Glacier during 1973—2020

本研究进一步分析了2011—2016 年，North Kyzkurgan 冰川末端的前进和表面特征变化过程。如图3 所示，2006 年8 月—2011 年10 月North Kyzkurgan 冰川末端出现了膨胀变宽现象［图3（a），黑色箭头处］，但冰川表面依然平整，表明冰川物质已开始挤压转移，但规模较为有限。2011 年10 月—2013 年10 月，冰川面积发生了迅速扩张，末端累计前进1.08 km，冰川表面裂隙广泛发育，剧烈破碎化，表明冰川物质由积蓄区向接收区迅速转移［图3（b）和3（c）］。2013年10月—2016年9月，冰川面积进一步扩张，末端继续前进，但前进的幅度已明显减弱，在2016 年9 月后末端前进停滞，表明冰川物质的转移过程逐渐减弱直至结束［图3（d）］。

3.2 冰面高程变化

本研究选择影像质量较优的2002 年、2009 年、2014 年和2020 年4 期ASTER 立体像对提取DEM，对North Kyzkurgan 冰川的表面高程变化状况进行了计算和分析。计算结果表明：2002—2009 年冰川跃动发生前［图4（a）］，North Kyzkurgan 冰川积蓄区上部表面高程基本保持稳定，末端表面高程降低，积蓄区下部高程却略有升高。同时期冰川的相应位置也发生了膨胀变宽，但冰川表面仍基本保持平整［图3（a）］，表明积蓄区下部冰体因受到了强烈挤压而膨胀变形［2］。2009—2014 年跃动发生期间［图4（b）］，North Kyzkurgan 冰川表面高程变化非常显著，冰川末端高程明显升高，最大升高幅度达到了近180 m，冰川中游高程则明显降低，表明在跃动过程中大量物质由积蓄区转移到了接收区，因此也导致了冰川末端的大幅度前进。2014—2020 年冰川跃动结束后［图4（c）］，North Kyzkurgan 冰川的表面高程普遍降低，表明由于跃动造成的冰面剧烈破碎化，导致冰川末端处于强烈消融状态。

图3 2011—2016年North Kyzkurgan冰川末端位置变化Fig. 3 Terminus positions at different time spans for North Kyzkurgan Glacier during 2011—2016

图4 2002—2009年、2009—2014年以及2014—2020年North Kyzkurgan冰川表面高程变化Fig. 4 Elevation difference of the North Kyzkurgan Glacier during 2002—2009，2009—2014 and 2014—2020

此外，Zhou 等［20］利用KH-9 DEM 和SRTM DEM 计算帕米尔中部冰川表面高程变化的结果表明，1975—1999年North Kyzkurgan冰川积蓄区增厚速率为（0.50±0.31）m·a-1，接收区减薄速率高达（1.87±0.31）m·a-1，而同时期帕米尔中部冰川的整体物质平衡水平却仅为（0.03±0.24）m w. e. ·a-1。这表明1975—1999 年在帕米尔中部冰川基本保持稳定的同时，North Kyzkurgan 冰川的表面高程却发生了显著变化，冰川末端处于强烈退缩状态。

3.3 冰面流速变化

本研究基于ITS_LIVE 产品计算和分析了2000—2018 年North Kyzkurgan 冰川表面的年平均［图5（a）］和年最大流速［图5（b）］变化。如图5 所示，在2000—2011 年，North Kyzkurgan 冰川运动速度十分缓慢且稳定，年平均流速介于7.38～15.47 m·a-1，年最大流速介于34.37～105.12 m·a-1。即便是冰川下游发生膨胀变宽的2006—2011 年，冰川流速也仅有微弱的增加。2011年之后伴随着冰川末端的快速前进，冰川流速也显著加快，2013年冰川表面最大流速达到最高约（400.60±1.91）m·a-1，2014 年平均流速达到最高约（39.39±0.54）m·a-1，是跃动发生前的数十倍。冰川最大流速出现的位置集中在中下游，并伴随着流速的增加向下游转移。总体上，North Kyzkurgan冰川表面流速在2013年前后达到顶峰，此后迅速减缓，在2016年基本恢复到了跃动发生前的水平。

图5 2000—2018年North Kyzkurgan冰川的年平均流速（a）和年最大流速（b）Fig. 5 Annual average and maximum surface flow velocities of the North Kyzkurgan Glacier from 2000 to 2018

4 讨论

4.1 冰川跃动过程

冰川跃动的发生具有周期性，从十几年到上百年不等，通常根据运动状态划分跃动周期内的恢复阶和跃动阶［1，32］。由于冰川跃动会引发面积、长度、高程、流速等一系列的变化，因此可以根据这些变化识别跃动冰川的恢复阶和跃动阶［10］。就North Kyzkurgan 冰川而言，2011 年以前，在整体上冰川面积持续退缩，表面高程缓慢降低，运动速度基本稳定；2011—2016 年，冰川末端迅速前进，面积大幅度扩张，积蓄区和接收区高程剧烈波动，冰川流速也显著加快，同时冰川表面裂隙广泛发育，剧烈破碎化；2016 年之后，冰川末端前进停滞，表面高程降低，流速也恢复到跃动发生前的水平。这些综合变化表明，North Kyzkurgan 冰川是一条典型的跃动冰川，1973—2011 年处于恢复阶，2011—2016 年处于跃动阶，2016 年之后重新进入了恢复阶，因此该冰川的跃动周期至少在40年以上。

此外，值得注意的是2006—2011 年North Kyzkurgan 冰川下游出现了膨胀变宽的现象［图3（a）］，这种现象的出现主要与冰川跃动的发生过程有关。冰川跃动是冰体内部应力集中释放的结果，由于积蓄区不断增厚，下部冰体受到的应力逐渐增加，当达到冰体所能承受的极限时，应力获得集中释放，即引发冰川跃动［33-34］。因此跃动过程先是冰川积蓄区下部挤压冰体产生缓慢变形，然后才引起上下游的快速运动［2，17，32］。图3和图4显示冰川膨胀变宽的位置出现在积蓄区下部向接收区的过渡处，且此时冰川表面仍基本保持平整，即表明这种现象的出现是上游物质积累挤压所致。在高亚洲其他冰川跃动事件中也观察到了同类现象，如喀喇昆仑山的Balt Bare 冰川和木孜塔格的鱼鳞川冰川［2，10，17］。在2006—2011 年North Kyzkurgan 冰川下游膨胀变宽的同时，冰川流速却基本保持稳定，同样表明此时尚处于快速跃动前。

帕米尔是全球跃动冰川集中分布的地区之一，现有的跃动冰川编目和相关研究表明，帕米尔跃动冰川存在较大的异质性，跃动阶从数月到十几年不等，跃动周期则从十几年到数十年不等［10-11］。总体而言，西帕米尔和东帕米尔跃动冰川的跃动阶和跃动周期相对较短，分别集中在数年和十几年，而帕米尔中部跃动冰川的跃动阶和跃动周期相对较长。西帕米尔的Medvezhy冰川分别在1937年、1951年、1963 年、1973 年、1989 年、2001 年、2011 年发生跃动，每次跃动阶长度不到1 年，跃动周期在10～14 年左右［35］。Oshanina 冰川分别在1961—1962 年、1983—1984 年、2010—2011 年发生跃动，其跃动周期稳定在25 年左右［10-11］。东帕米尔昆盖山的G074348E 39282N 冰川跃动阶为3 年左右，两次跃动的间隔也仅为10 年［12］，5Y663L0023 冰川跃动阶为4 年，恢复阶最短为15 年左右［14］。本研究中位于帕米尔中部的North Kyzkurgan 冰川跃动阶为6 年，跃动周期则达到了40 年以上，而同样位于帕米尔中部的Bivachny冰川分别于1975—1978年、2011—2014年发生跃动，两次跃动间隔也在35 年左右［11，13］，因此从东、西帕米尔向帕米尔中部，跃动冰川的跃动周期呈现增加趋势。这种现象在一定程度上与地形因素有关，东、西帕米尔山地，如科学院山脉、昆盖山等地区的跃动冰川上游往往存在大面积的积雪陡壁，在积蓄区接受降雪直接补给的同时也接受雪崩补给，因此物质积累过程迅速，跃动周期相对较短［4，36］。而研究冰川所在的帕米尔中部腹地，地形平缓，缺乏雪崩补给（图1），只能依赖降雪的缓慢积累，因此跃动周期较长。

4.2 跃动发生的潜在因素

冰川跃动的发生机理较为复杂，可以将跃动冰川分为温冰川和多温型冰川［37］。温冰川的跃动机制主要由水文模型解释，夏季冰川消融产生的大量融水进入冰体和冰下，水量增加导致静水压力增加，排水系统发生扩张，当融水减少、静水压力下降时引发排水通道坍塌，在冰川底部融水的润滑作用下触发冰川跃动，因此温冰川的跃动周期短、运动速度快、通常在冬季开始夏季结束［38-39］。多温型冰川的跃动机制主要由热力模型解释，冰川上部物质不断积累、应力不断增加，导致底部达到压力熔点，产生大量融水，而冰川下游的冰体和冻土阻碍了融水流失，使融水在冰床不断蓄积，在融水的润滑和顶托作用下，冰川发生快速滑动即冰川跃动，因此多温型冰川的跃动周期长、运动速度慢，在任何季节都可开始与结束［40］。

North Kyzkurgan 冰川所在的帕米尔中部地区气候寒冷，降水量丰富且以降雪为主，Fedchenko 气象站1935—1990年的气象观测资料显示，该地区年平均气温为-6.9 ℃［图6（a）］，年均降水量高达1 200 mm［图6（b）］。而近三十年来的气象观测资料以及CRU 再分析资料都表明，帕米尔中部气温呈上升趋势［0.1 ℃·（10a）-1］，降水量则基本保持稳定［41-42］。与此同时，North Kyzkurgan 冰川的积累区面积比率超过0.8，冰川作用正差达到1 000 m，因此极低的气温和丰富的降雪为该冰川提供了充足的物质补给。跃动冰川的恢复阶长短一般取决于降雪填平积蓄区所需的时间［43］，而North Kyzkurgan冰川缺乏雪崩补给，物质积累几乎完全依赖降雪直接补给，因此跃动周期漫长。基于DEM 和遥感影像的观测结果显示，在1973—2011 年长达近40 年的时间里，North Kyzkurgan 冰川积蓄区表面高程不断升高，表明这一时期该冰川获取了充分的物质积累。在积蓄区物质不断增多，冰川底部受到的应力不断增大以及近几十年气温升高的共同影响下，冰川底部达到了压力熔点，融水不断增加，最终触发了冰川跃动。综上所述，North Kyzkurgan 冰川跃动的发生是物质长期积累的结果，主要由热力学因素所致。

图6 1935—1990年研究区年平均气温（a）和年降水量变化（b）Fig. 6 Mean annual air temperature and precipitation for 1935—1990 in study region

5 结论

本文利用1973 年以来Landsat 影像、ASTER 立体像对和ITS_LIVE 数据产品，对帕米尔中部North Kyzkurgan 冰川开展联合监测，从冰川面积、高程、流速等方面分析了该冰川跃动前、跃动中和跃动后的变化特征，揭示了该冰川的完整跃动过程。研究结果表明：

North Kyzkurgan 冰川是一条典型的跃动冰川，在1973—2011年处于恢复阶，2011—2016年处于跃动阶，2016 年之后重新进入恢复阶。North Kyzkurgan 冰川所在的帕米尔中部地区气候寒冷，降雪量丰富，同时该冰川的积累区面积比率超过0.8，冰川作用正差近1 000 m，因此冰川补给物质充足。随着积蓄区不断增厚，冰川底部达到压力熔点，融水不断增多，在融水的润滑和顶托作用下，最终导致冰川发生跃动。因此，North Kyzkurgan 冰川跃动的发生主要与热力学因素有关。