中喀喇昆仑布拉尔杜冰川近期跃动分析

王振峰， 蒋宗立， 刘时银， 杨婧睿， 马致远， 王雯清

（1.湖南科技大学测绘遥感信息工程湖南省重点实验室，湖南 湘潭 411201；2.湖南科技大学地球科学与空间信息工程学院，湖南 湘潭 411201；3.云南大学国际河流与生态安全研究院，云南 昆明 650500）

冰川跃动是冰川准周期性的运动，在短暂的跃动阶段（通常以高流速为特征）和漫长的静止阶段（以末端停滞或退缩为特征）之间交替［1-2］。在跃动阶段，大量的冰川物质从上游积蓄区快速向下游接收区迁移，导致冰川表面高程发生剧烈变化。大规模的冰川物质再分配可能会导致冰川末端前进、冰川湖溃决洪水（GLOF）等一系列灾害，从而威胁流域内居民的生命和财产安全。例如，2015年东帕米尔地区的克拉牙依拉克冰川发生跃动，引发的雪崩造成大量的牧场和牧民的房屋被破坏［3］；2016 年西藏阿鲁错冰川跃动/冰崩造成牧民和牲畜死亡［4］。因此，研究冰川跃动的过程及控制机理对于预防和减轻跃动带来的灾害具有重要理论和实践意义。

跃动型冰川并未均匀的分布在世界冰川区，它们往往聚集在少数地区［5-6］，如阿拉斯加及加拿大育空、斯瓦尔巴德群岛、格陵兰岛以及喀喇昆仑、天山、帕米尔等中亚部分山系［7-9］，俄罗斯的高加索和新西伯利亚群岛等地区也曾出现过关于冰川跃动的报道［10］。虽然跃动型冰川仅约占全球冰川总数的1%，但是其重要影响却不可忽视［5］。

目前存在两种普遍接受的假说解释冰川跃动：冰下水文机制模型和冰下热转换机制模型［11-12］。冰下水文机制模型是在阿拉斯加杂色（Variegated）冰川1982—1983年跃动研究的基础上发展起来的［13］，水文控制的跃动主要由于冰下孔隙水压力的改变而导致的不稳定，通常伴随着冰下排水通道的重组，这类跃动通常在排水效率低下的冬季开始，结束于排水效率高的夏季。热转换机制模型是在斯瓦尔巴德地区冰川的研究中发展起来的［14］，热转换机制的跃动是由于冰川积蓄区的物质不断积累致使冰川底部温度达到融点，冰川底部静水压力的上升致使底部快速滑动，当底部融水全部流失后跃动结束。由于冰川底碛特性很难观测，大部分地区的跃动机制仍不明确。

在全球变暖的气候背景下，全球冰川整体上呈现不同程度的持续退缩，而喀喇昆仑地区由于异常的区域物质平衡［15］和大量跃动型冰川［16］倍受关注。研究表明喀喇昆仑的跃动机理存在异质性，一方面认为喀喇昆仑地区的跃动是受热转换机制控制，Quincey 等［17］利用卫星遥感对喀喇昆仑地区5 条冰川的跃动过程进行研究，量化了这5 条冰川流速的时空变化，并认为喀喇昆仑地区近期的冰川跃动是由热转换机制控制的，且与该地区降水增加和冰川积累模式导致的高海拔地区变暖相一致；另一方面对单个冰川跃动的事件观测表明冰下水文条件的变化是主要触发机制，如Bhambri 等［18］基于多源遥感数据对喀喇昆仑地区Shispare冰川跃动前后的流速、高程和表面特征的变化进行分析，认为Shispare冰川的跃动是由冰下水文机制控制。

目前喀喇昆仑地区详细的冰川跃动过程描述仍然缺乏，至今仅有Shispare、Hispar、Kyagar、Khurdopin 以及Baltoro 等冰川的跃动过程有相对详细的遥感观测［18-22］。随着遥感技术快速的发展使得对跃动冰川的识别和研究更加快捷，利用获得的冰川跃动过程及特征可以推断冰川的跃动控制机制。本研究收集了布拉尔杜冰川跃动期间Sentinel-1 数据，基于雷达影像偏移追踪法提取了此冰川表面运动速度，此外还基于TerraSAR-X/TanDEM-X 数据生成的数字高程模型（DEM）对冰川表面高程变化进行量化，描绘其跃动特征，进而对其控制机理进行探讨，以期对喀喇昆仑布拉尔杜冰川跃动的进一步认识，为喀喇昆仑冰川跃动研究提供更多的实例。

1 研究区概况

布拉尔杜冰川位于喀喇昆仑中部，其融水汇入什玛莎河，上游东侧与音苏盖提冰川相接。喀喇昆仑地区主要受亚洲季风影响，夏季降水量占总降水量的80%。布拉尔杜冰川编号为RGI60-14.04411，中心位置为36°07′N，75°52′E（图1），此冰川主体呈现自南向北的流向，是由多条支冰川组成的树枝状山谷冰川，面积约为176 km2［23］。冰川末端被大量表碛所覆盖，表面还育有少量冰面湖。

图1 布拉尔杜冰川位置Fig.1 Location of the Braldu Glacier

2 数据与方法

2.1 数据来源

2014 年发射的Sentinel-1A 卫星搭载C 波段传感器，具有全球采集、重访周期短和良好的空间分辨率等特点，使其成为研究冰川运动的重要数据源。本研究利用Sentinel-1A 干涉宽幅（IW）模式的影像获取布拉尔杜冰川2014—2018年的流速变化，并结合ITS_LIVE 项目（https://nsidc.org/apps/itslive/）中的高亚洲地区冰川的年均流速数据进行分析。

TerraSAR-X/TanDEM-X（TSX/TDX）系统是由德国宇航局（DLR）先后发射的两颗X波段SAR卫星组成的单轨双天线双星分布式系统。两颗卫星具有相似的结构，且采用近距离螺旋结构运行，重复周期为11 d。该系统在双基站模式下，可以同时采集干涉对的主从影像，具有“0 时间基线”的特点。SRTM是在2000年由美国国家航空航天局和国家影像与测绘局联合进行的航天飞机雷达地形测绘任务，SRTM 的观测宽度为225 km，其生成了全球约80%陆地表面的数字高程模型。本文使用TerraSAR-X/TanDEM-X 影像提取新的DEM，与SRTM DEM 共同用于计算布拉尔杜冰川跃动前后表面高程的变化，同时SRTM DEM也作为参考DEM对SAR数据进行基于地形的精配准。

Landsat数据通过美国地质调查局（USGS，https://earthexplorer.usgs.gov/）进行处理和分发。本文选取了夏秋两季云雾较少的影像，主要用于监测布拉尔杜冰川跃动期间冰川末端和内侧冰碛的变化。本研究中使用数据见表1。

表1 使用的数据源Tab.1 Datasets used in this study

2.2 冰川表面流速与高程提取

2.2.1 冰川表面流速提取由于冰流速度过快而导致的冰面失相干使得D-InSAR难以获得高质量的干涉条纹图，而偏移量跟踪法中的强度追踪法却可以克服冰面失相干的问题，被认为是获取冰川流速的可靠手段［24-25］。强度追踪的原理是根据图像的纹理（条纹或后向散射强度）的最大相关系数来计算目标的偏移量。本文采用瑞士GAMMA合成孔径雷达软件中的Offset Tracking 模块对Sentinel-1A 影像进行偏移量估计，具体流程如图2。在偏移量追踪计算时需要对搜索窗口的大小进行设定，窗口尺寸过小会导致单个窗口不稳定，可能无法在其他图像中找到对应的窗口。经试验，将匹配窗口设置为256×64效果最佳。

图2 偏移量追踪方法获取冰川表面流速流程Fig.2 Schematic diagram for acquirement of the surface velocity of glaciers using offset tracking method

本文利用SAR 影像偏移量追踪技术获取的冰川流速，其误差主要来自图像配准、图像对的偏移估计及系统误差。为了对冰川表面流速的精度进行评估，通常假设非冰川区为静止区域，静止区域包含上述所有误差源。因此通过非冰川区残差位移对冰川流速结果进行评估，具体公式［26］为：

式中：eoff为非冰川区偏移量的误差；MED为非冰川区偏移量的平均值；SE为非冰川区偏移量的标准误差；SD 为非冰川区偏移量平均值的标准差；Neff为去自相关后的像元个数；Ntotal为非冰川区像元总数；R为像元分辨率；D为去除空间自相关距离。根据上述公式可得非冰川区偏移量误差为0.016 m·d-1，远小于冰川速度的观测值。

2.2.2 冰川表面高程提取及变化计算在瑞士GAMMA软件平台下对TSX/TDX SAR数据进行差分干涉处理，在差分干涉处理的过程中引入1 rad·s-1SRTM DEM 作为参考DEM，从而获取TSX/TDX DEM。主要处理步骤：首先将SRTM DEM 模拟转换为SAR 成像数据，其次TSX/TDX 与SRTM 模拟的SAR 数据进行配准并共轭相乘生成干涉图［27］，通过差分处理，去除SRTM DEM 模拟的地形相位和基线误差导致的线性相位，从而得到冰川表面高程变化的残差相位，将得到的残差相位加上SRTM 模拟的地形相位，最后经过地理编码，即可生成高精度的TSX/TDX DEM［28］，具体流程如图3。将生成的TSX/TDX DEM与SRTM DEM进行配准后，采用大地测量法计算2000—2014年、2014—2018年布拉尔杜冰川的表面高程变化。因SRTM-X 未完全覆盖研究区，故本文选取SRTM-C进行研究。由于微波对冰雪具有一定的穿透性，利用InSAR 得到的冰川表面高程可能低于实际高程。考虑到TSX/TDX与SRTM-X的载频几乎相同，因此在进行冰川高程变化计算时，可利用SRTM-X DEM对C波段的冰雪穿透深度进行估算和校正［29］，经计算得出SRTM-C 在喀喇昆仑地区的穿透深度约为2.4 m［30］。

图3 TSX/TDX差分干涉获取DEM流程Fig.3 Flow chart for DEM generation from TSX/TDX differential interferometry

冰川表面高程变化的不确定性评估方法与上文冰川表面速度的不确定性评估类似，假定非冰川区高程无变化，利用非冰川区的高程差对DEM数据间相对误差进行评估。最终不同DEM 数据间高程误差（σ）可用下列公式表示［26］：

式中：N为去自相关后的像元个数，经计算得到2000—2014 年和2014—2018 年非冰川区表面高程变化的总体误差分别为0.85 m和2.41 m（图4）。

图4 非冰川区表面高程变化Fig.4 Elevation change in non-glacier area

3 结果与分析

3.1 冰川表面流速变化

利用偏移量追踪的方法获取了布拉尔杜冰川2014年10月—2018年1月的时序表面流速数据，为了进一步地观察冰川流速的时间变化，本研究沿冰川中心线提取速度剖面图（图5）。从图中可以看到冰川中上游在2014 年10 月已具有较高的流速，最高流速达3.6 m·d-1，此时布拉尔杜冰川处于高速运动期。随后冰川的跃动前锋保持较高的流速向下游推进，在此过程中最高流速达4.9 m·d-1。到2015年8月底，冰川表面流速急剧下降，跃动前锋的推进随着冰川流速的下降而终止，冰川进入平静期。然而到2016 年1 月初，冰川表面流速突然增加，最高流速达到3.9 m·d-1，这种高流速的状态维持了一个月，随后冰川进入平静期。结合ITS_LIVE 项目（https://nsidc.org/apps/itslive/）中的高亚洲年均流速数据对布拉尔杜冰川进行分析（图6），发现从2006年起，布拉尔杜冰川主干中上游的流速逐年升高，具有明显的跃动前锋并向下游缓慢推进，此时冰川已经显示出跃动特征趋势。布拉尔杜冰川在2012—2015年处于快速运动期，在此期间最高流速可达400 m·a-1。布拉尔杜冰川由多个支冰川汇流而成，图7为布拉尔杜冰川详细的冰川流速数据，从图7a～d、图7g中可以观察到在跃动期间不仅冰川主干具有较高的流速，支冰川流速也相对较高；而在平静期，布拉尔杜冰川主干和其支冰川均处于缓慢流速的状态（图7e～f、图7 h～i）。

图5 布拉尔杜冰川沿主流线日均流速变化Fig.5 Variation of daily average velocity along the mainstream of the Braldu Glacier

图6 1989—2018年布拉尔杜冰川沿主流线流速变化Fig.6 Velocity variation along the mainstream of the Braldu Glacier during 1989—2018

3.2 冰川表面高程变化

为了定量分析冰川跃动前后表面高程的变化，本文基于大地测量法获得2000—2014 年、2014—2018 年布拉尔杜冰川的表面高程变化（图8），由于2014 年的TSX/TDX 数据没有完全覆盖布拉尔杜冰川，因此导致图8 部分区域缺失。2000—2014 年3、4 号支冰川向主干输送大量的冰川物质，致使冰川主干中上游有明显隆起，1、2、5、6 号支冰川与主干相接的位置也明显增厚，而冰川主干的下游以消融为主，其高程明显降低。2014—2018年布拉尔杜冰川主干中上游以及各支冰川有不同程度的减薄，而冰川主干的接收区明显增厚，最大增厚可达120 m。

沿冰川中心线提取布拉尔杜冰川2000—2014年和2014—2018 年的表面高程变化（图9），结果表明在2000—2014 年，冰川的下游（0～13 km）以消融为主，最大减薄89 m，而冰川主干中上游（距末端13～26 km）有不同程度增厚；2014—2018年，冰川主干的接收区（6～13 km）高程显著增加，冰川主干中上游减薄5～20 m不等。

图9 布拉尔杜冰川沿主流线表面高程及高程变化Fig.9 Elevation and elevation change profile of the Braldu Glacier along the mainstream

4 讨论

ITS_LIVE 年均流速（图6）表明此冰川从2006年开始缓慢加速，并逐渐向下游推进。根据Landsat影像发现5、6 号支冰川的冰碛从2006 年起逐渐向外扩展，结合图8a 中的高程变化可以发现5、6号支冰川与冰川主干相接位置有明显的增厚，故认为5、6 号支冰川从2006 年起便处于跃动状态，这与Gardelle 等［31］研究的结果相一致，但他们的研究中认为1 号支冰川也处在跃动期，而本文仅通过2000—2014年高程和表碛的变化无法判别1号支冰川是否从2006年开始跃动。从2009年开始，3、4号支冰川的冰碛（距末端23 km处）逐渐扩大并沿主干向下游推进，结合冰川表面流速的结果（图5、图6）发现布拉尔杜冰川从2011 年开始处于高速运动的状态，随着支冰川和主干上游的冰川物质向下游推进，冰川主干中部的表碛受到挤压发生扭曲并向下游继续推进。而布拉尔杜冰川表面冰碛的变化是由于表面速度的不均匀性所导致，同时也表明冰川处于跃动阶段［32］。

图8 布拉尔杜冰川表面高程变化Fig.8 Change of the Braldu Glacier surface elevation

通过冰川表面流速和冰碛的变化发现布拉尔杜各支冰川跃动的起始时间并不一致，根据2000—2014 年布拉尔杜冰川高程变化结果（图8a）可以推断布拉尔杜冰川的各分支在2014 年之前均已开始跃动。图8a 的结果显示在2014 年之前仅有3、4 号支冰川的物质向下游传递，致使冰川主干的中上游（13～20 km）明显增厚；1、2、5、6 号支冰川虽在此之前已经发生跃动，但来自这四条支冰川上游的物质大多聚集在与冰川主干相接的位置，故此认为3、4号支冰川的跃动引起主干冰川的第一阶段跃动。而在2014年之后，冰川主干的中上游和各支冰川的冰川物质沿冰川主干向下游传播，导致冰川主干的接收区（6～13 km）的高程显著增加，最大增厚约120 m，而1、2、5、6号支冰川明显减薄，由此认为1、2、5、6号支冰川的跃动引起第二阶段的跃动。

由于缺乏2006—2013年详细的流速数据，无法获取布拉尔杜冰川跃动初期流速变化的特征。Quincey等［33］利用光学影像特征匹配法获取2014年布拉尔杜冰川夏季表面流速数据，发现8 月的流速与7 月相比大幅下降，且通过现有的流速数据也发现此冰川在2015年8月底突然减速，在保持一段低速运动后，于2015—2016 年冬季突然加速，随后逐渐减速。布拉尔杜冰川流速的季节性变化比较符合冰下水文控制冰川的跃动特征，速度在夏季急剧下降，在冬季会上升［34-35］。受水文控制的冰川跃动，流速与冰下水压直接相关，流速的大小取决于冰下排水系统的状态［36］。在夏季，大量融水的输入致使排水渠道高效的运转，冰下水压下降，导致流速下降。由于冰体的持续压缩运动和可塑性，排水通道在秋冬季节随着融水的减少而变得低效运转，随后冰下水压升高，最终冬季冰流加速。

与喀喇昆仑其他跃动冰川相比，布拉尔杜冰川在跃动期间并未表现出较高的流速，最高流速仅4.9 m·d-1。而与布拉尔杜冰川相邻的Hispar冰川在2015 年春季的最高流速可达14 m·d-1，此冰川在2015 年夏季被观察到流速突然下降，却在2015—2016年冬季再次加速，其跃动也被认为是与冰下水文状况有关，且Hispar 冰川此次跃动主要受其支流影响［19］。Quincey 等［33］通过对一些喀喇昆仑冰川跃动过程研究发现，喀喇昆仑地区的冰川跃动并不完全符合经典的水文控制模型或热控制模型，地形和表碛物也是影响跃动的重要因素。布拉尔杜冰川表面覆盖着厚厚的表碛和死冰，这些物质可能会阻碍冰川物质向下游输送。

此外，张祥松［37］曾观察到布拉尔杜冰川的末端在1968—1973 年前进了约900 m，并且前进的冰舌与韦斯米亚兹河口的冲积扇相接；Rankl 等［38］利用SRTM、TerraSAR-X、ALOS PALSAR 和ENVISAT ASAR等多源雷达影像对喀喇昆仑冰川变化进行研究，发现布拉尔杜冰川上次跃动发生在1976年之前，因此推断布拉尔杜冰川的跃动间隔约为40 a。

5 结论

（1）根据ITS_LIVE年均流速数据发现，从2006年开始，布拉尔杜冰川年均流速逐渐增加，展现出跃动趋势。

（2）结合布拉尔杜冰川高程变化数据发现，在2014 年之前，此冰川的各支冰川均已开始跃动，但仅有3、4 号支冰川持续向下游进行物质传输，而其余支冰川在2014 年之后开始向下游进行物质的传输，可认为由于3、4 号支冰川的跃动引起主干冰川的第一阶段跃动，1、2、5、6 分支在2014—2018 年物质向下游迁移引起第二阶段跃动。

（3）根据布拉尔杜冰川2014—2018 年详细的流速数据，发现此冰川流速的变化特征比较符合水文控制冰川的跃动特征。结合现有的文献及影像数据，推断布拉尔杜冰川的跃动间隔约为40 a。