2014—2020年科其喀尔巴西冰川运动速度年际与季节变化特征分析

2022-02-12赵晋彪许杨杨王荣军蒋宗立

冰川冻土 2022年6期

赵晋彪，张震，，许杨杨，王荣军，蒋宗立

（1.安徽理工大学空间信息与测绘工程学院，安徽淮南 232001；2.安徽理工大学矿山采动灾害空天地协同监测与预警安徽普通高校重点实验室，安徽淮南 232001；3.安徽理工大学矿区环境与灾害协同监测煤炭行业工程研究中心，安徽淮南 232001；4.中国科学院西北生态环境资源研究院冰冻圈科学国家重点实验室，甘肃兰州 730000；5.湖南科技大学地球科学与空间信息工程学院，湖南湘潭 411201）

0 引言

运动是冰川的固有属性，是冰川区别于其他冰体的主要标志［1］。冰川运动在其动力学系统的控制下，通过改变物质的分布，影响冰川几何形状、物质平衡、水热及边界的变化［1］。冰川运动速度作为冰川动力学模型的重要参数，不仅能响应物质平衡变化，同时也能反应冰川的动态信息［2-3］。在全球变暖的背景下，研究冰川运动及其变化规律，对区域经济发展、资源合理利用和自然灾害预警有着重要意义。

冰川运动速度的提取方法主要包括实地测量和遥感监测［4-5］。实地测量是在冰面上布设花杆，通过花杆的位置变化测量冰川运动速度。尽管实地测量能获得精度高、连续性好的数据，但由于山地环境恶劣和成本的限制，无法进行长周期、大范围的观测［6-7］。遥感监测具有获取速度快、监测范围广的优势，可弥补实地测量的不足［8-9］。基于遥感影像的冰川运动速度提取方法主要有雷达差分干涉测量、SAR 偏移量跟踪和光学影像互相关［10］。由于山地冰川地形复杂且消融较为强烈，获取高相干影像较为困难，因此很少采用雷达差分干涉测量［11-12］。基于SAR 数据的偏移量跟踪方法受天气影响小，但长时间冰川消融导致较大表面差异，较适合短时间尺度的冰川运动速度提取［5，13］。而光学影像互相关方法虽受天气影响，但历史资料丰富，已被广泛应用于年尺度的冰川运动速度研究中［14-15］。

随着遥感数据时空分辨率的不断提高，基于遥感数据提取冰川运动速度已成为当前研究的热点［16-18］。部分学者利用遥感数据对大区域尺度的冰川运动特征进行分析，比如Sam 等［19］基于GoLIVE数据分析了喜马拉雅山冰川运动特征；Dehecq 等［20］基于Landsat 等数据分析了亚洲高山区冰川运动变化特征。也有部分研究者针对单条冰川进行时空变化特征分析，比如Das 等［21］对米亚冰川1992—2019 年的运动速度时空分布特征及影响因素进行分析；周中正等［22］基于Sentinel-1A 影像分析了2016年岗纳楼冰川表面运动速度变化。

科其喀尔巴西冰川表面分布有大量表碛，对其运动速度进行研究，有助于了解表碛覆盖型冰川对气候变化的复杂响应。以往的研究更侧重于冰川运动速度的空间分布以及区域性变化［23-24］，对于冰川运动速度年内特征及年际变化特征的认识仍然不足。因此，本研究利用Landsat 8影像提取2014—2020 年科其喀尔巴西冰川的运动速度，分析其年内波动特征和年际变化特征，为冰川动力学模型及物质平衡变化提供基础支撑。

1 研究区概况

科其喀尔巴西冰川（冰川编号：5Y674A0005），位于新疆温宿县，属亚大陆型冰川。冰川表面高程范围为3 020～6 342 m，冰川总长约25.1 km，面积83.56 km2，冰储量15.79 km3，物质平衡线海拔高度为4 350 m［25］。该冰川消融期较长［26］，从4 月中下旬开始消融，9 月中下旬消融结束。冰川区受潮湿气流影响［27］，降水丰富，5—9 月降水量约占全年降水量的75%。据冰川编目资料显示［25］，该冰川消融区分布有大量表碛，表碛面积为17.5 km2，占整个冰川面积的21%，其表碛厚度随海拔降低逐渐增加，在冰川末端厚度达到200 cm［28］。

2 数据与方法

2.1 数据

2.1.1 Landsat 8

本研究使用Landsat 8 L1T 影像数据，该级别影像是经过辐射校正和几何精校正处理的数据。根据影像质量及数据可获取性，结合表面积雪状况及研究时段间隔，主要采用冰川表面积雪较少，冰川纹理清晰的数据用于获取冰川运动速度。共选取2014—2021 年间的16景Landsat 影像（行列号：147，31），用于冰川运动速度提取。并按照影像成像时间间隔组成影像对（表1），时段1～7为年运动速度影像对，时段8～21为季节运动速度影像对。数据获取网站（https：//earthexplorer.usgs.gov/）。

图1 科其喀尔巴西冰川位置示意图Fig.1 Location of Koxkar Baxi Glacier

表1 用于提取科其喀尔巴西冰川运动速度的Landsat 8影像对Table 1 The Landsat 8 image pairs for extracting Koxkar Baxi Glacier flow velocity

2.1.2 ASTER GDEM

ASTER GDEM（ASTER Globe Digital Elevational Model）由美国NASA 和日本METI 共同研制并面向全球免费发布。该产品基于ASTER 近红外影像数据生成，其水平空间分辨率为30 m［29］，数据覆盖了全球99%的陆地表面（83° N～83° S）。本研究使用的是第三版数据（ASTER GDEM V3），主要用于获取研究区高程信息。该产品的发布网站为NASA EARTH DATA（https：//search.earthdata.nasa.gov/search）。

2.1.3 冰川编目数据

冰川边界数据来源于中国第二次冰川编目数据集（http：//www.ncdc.ac.cn），该数据是以Landsat TM/ETM+和SRTM V4 为数据源，采用高效波段比值阈值分割方法，并结合人工修订所获得。

2.1.4 其他数据

气候数据是来自欧洲中期天气预报中心（ECMWF）的第五代大气再分析产品（ERA5），是一种综合性的再分析数据，本文中利用大气再分析产品获取气温和降水数据（数据获取网站为https：//cds.climate.copernicus.eu/#！/home）。

ITS_LIVE（The inter-mission Time Series of Land Ice Velocity and Elevation）数据集［30］提取自Landsat 4、5、7、8 卫星影像，涵盖了所有面积大于5 km2的陆地冰区，时间跨度为1985—2020年。该数据有120 m和240 m 两种分辨率，本文使用由光学卫星图像对生成的120 m 分辨率的速度数据（该数据获取网站为https：//its-live.jpl.nasa.gov/）。

GoLIVE（Global Land Ice Velocity Extraction from Landsat 8）数据集［31］是将Landsat 8影像在空间域中进行归一化互相关，进而提取的冰川运动速度，经过误差调整后的数据精度可达到亚像素级别。该数据集覆盖了82° N～82° S 之间的所有陆地冰川，空间分辨率为300 m，并提供多种时间跨度的数据（16 d、32 d、48 d、64 d、80 d和96 d）。本文利用多个时间跨度为32 d 的数据进行计算，获取年尺度冰川运动速度（该数据集获取网站为https：//nsidc.org/data/）。

Sentinel-1 冰川运动速度数据集［32］是通过强度偏移跟踪技术，从Sentinel-1 SAR 影像中提取的冰川运动速度，该数据集涵盖了极地冰盖以外的12个主要冰川区域，并提供分辨率为200 m 的月度和年度速度数据。该数据集可从埃尔朗根-纽伦堡大学冰川入口网站获取（http：//retreat.geographie.uni-erlangen.de.）。

实测数据是通过观测冰面上布设的花杆所获得，2015 年5 月至7 月在冰川表面海拔3 100～3 900 m区域布设了40根花杆（共布设40根花杆，其中10、12、39号杆丢失），2015年7月3日、9月28日、9月29日、9月30日分别对每根花杆采用2台水平精度和垂直精度分别为5 mm+1 ppm 和10 mm+2 ppm的X20单频静态接收机进行GPS测量。实测中将一台GPS 接收机作为基准站（位于科其喀尔巴西冰川大本营内），而另一台GPS接收机作为流动站在冰川表面进行流动观测。将观测数据导入华测GPS数据后处理软件中，进行求差解算坐标，从而获取每根花杆的位置信息，即经纬度坐标，最后利用花杆坐标位移量以及测量时间间隔即可计算得到速度信息。

2.2 研究方法

提取冰川运动速度的流程主要包括两部分，获取冰川的位移量以及运动速度的计算。本研究使用COSI-Corr（Co-registration of Optically Sensed Images and Correlation）软件包来获取位移信息，该软件包是加州理工学院基于IDL（Interactive Data Language）语言开发的，最初被应用于检测同震形变［30］，后被广泛用于处理Landsat、ASTER、SPOT 和Quickbird卫星影像。通过COSI-Corr软件对影像进行配准和互相关，其结果精度可达到亚像元级别。在积雪区域，该方法仍能获取精度较高的位移信息，已被证实是基于光学影像监测冰川运动速度的可靠方法［21，33-34］。

本研究使用的Landsat 8 L1T 产品是经过几何精校正的数据，各景数据间的匹配精度较高，无需再次配准，可直接用于位移信息的提取。将参考影像和搜索影像进行傅里叶变换后，通过设置合适大小的窗口和步长，滑动计算两景影像上的最大相关系数，进而获得整个区域在影像成像期间的相对位移［10］。该方法相关参数设置如下，参考窗口设为128，搜索窗口设为32，步长设为8。此过程生成3个图层：东西向位移（EW）、南北向位移（NS）、信噪比（SNR）。信噪比为从0到1的数值，值越大表示计算出的位移量越可信［2，35］。为保证结果的准确性，选择SNR≥0.9的部分。受云和阴影的影响，部分值不能准确反映冰川表面的位移信息，对其进行剔除后，即可得到冰川表面位移量。

利用冰川表面位移量计算运动速度Vi：

式中：Di为每对影像之间测量的位移；Ti是两景影像间的时间间隔。

2.3 精度评估

假设无冰区没有运动，那么可将基于遥感图像提取的无冰区运动速度看作误差值［36］。误差评估公式为：

式中：E为冰川运动速度误差；MED为无冰区运动速度的平均值；SE为无冰区运动速度的标准误差，其具体计算公式为：

式中：STDV为无冰区运动速度的标准差；Neff为去除空间自相关影响的有效像元个数，其计算公式为：

式中：Ntotal为无冰区像元的总数；R为像元分辨率；D为去除空间自相关距离，本研究使用Bolch 等［37］采用的保守估计，取空间去相关距离为像元分辨率的20倍。由此可得冰川运动速度误差为0.007 m·d-1，远小于冰川的实际运动速率。无冰区的误差分布直方图如图2所示。

图2 无冰区残余运动频率分布Fig.2 The frequency distribution of residual velocity in non-glacierized regions

3 结果分析

3.1 冰川运动速度空间分布特征

科其喀尔巴西冰川运动速度的空间分布如图3所示，2014—2020 年科其喀尔巴西冰川运动速度的空间分布较为相似。科其喀尔巴西冰川运动速度沿中流线分布特征为中上游快、下游慢。冰川末端及下游（0～10 km）冰川运动较为缓慢（0.01 m·d-1）；冰川中上游（10～18 km）运动速度随着海拔升高逐渐加快并达到最大值，最大速度为0.17～0.20 m·d-1；之后冰川运动速度随着海拔升高逐渐减小至0.01 m·d-1。以上符合冰川运动速度的分布规律，即从积累区上游至物质平衡线附近冰川运动逐渐加快，物质平衡线附近冰川运动最快，从物质平衡线附近至末端逐渐减小。

图3 2014—2020年科其喀尔巴西冰川运动速度空间分布Fig.3 Spatial distribution of flow velocity on the Koxkar Baxi Glacier from 2014 to 2020

科其喀尔巴西冰川运动速度沿横向剖面线的分布情况如图4（b）所示。冰川横剖面X1 运动速度明显呈现“中间快两边慢”的特点，中间最大运动速度为0.10～0.15 m·d-1。冰川横剖面X2 和X3 特征虽不及X1 显著，但也是中间快，中间最大运动速度分别为0.06～0.08 m·d-1和0.02～0.03 m·d-1，边缘最小运动速度小于0.01 m·d-1。

图4 冰川运动速度剖面图Fig.4 Glacial velocity profile：the velocity along the centreline（a）；transects of surface velocity（b）

3.2 冰川运动速度季节性变化特征

为了探讨冰川运动速度的季节性变化，本研究根据气温、降水的变化定义了暖季和冷季，将每年5—9 月作为暖季、9 月—次年5 月作为冷季。并提取该冰川沿中流线的运动速度对季节性变化特征进行分析，结果显示（图5），该冰川运动速度具有季节性波动特征，暖季运动速度均值为0.07 m·d-1左右，冷季运动速度均值为0.06 m·d-1左右，暖季比冷季快16.67%。

图5 沿中流线季节运动速度箱线图Fig.5 Box diagram of seasonal velocity along the center line

3.3 冰川运动速度年际变化特征

2014—2020 年科其喀尔巴西冰川年平均运动速度为0.04～0.05 m·d-1左右，运动速度的峰值均出现在距离冰川末端约18 km处（图6），最大运动速度在0.17～0.20 m·d-1之间。对科其喀尔巴西冰川中流线上年均运动速度进行统计学分析，结果显示（图7），线性回归模型R2为0.284，显著性P值为0.218，2014—2020 年该冰川年际运动速度变化特征不显著，但整体呈微弱减小趋势，平均运动速度减小约0.01 m·d-1。其中2014—2016 年冰川运动速度呈微弱减慢趋势（减小约21.3%），2016—2018 年冰川运动速度呈加快趋势（加快约25.6%），2018—2020 年冰川运动速度呈减小趋势（减小约27.6%）。

图6 沿中流线冰川年运动速度Fig.6 Annual glacier velocity along the center line

图7 2014—2020年均运动速度变化趋势Fig.7 The variation trend of flow velocity from 2014 to 2020

4 讨论

4.1 不确定性分析

由于实测数据有限，基于遥感方法的冰川运动速度误差很难评估［10，38-39］。一般假设无冰区在研究时间内保持稳定且未发生运动，那么无冰川区域的测量值可以代表冰川运动速度的误差。结果表明，冰川运动速度误差值为0.007 m·d-1，误差结果整体呈正态分布（图2）。

为了进一步评估本研究的数据质量，以2015年冰川运动速度数据为例，将其与有限的花杆实测数据、ITS_LIVE、GoLIVE 以及Sentinel-1 冰川运动速度数据进行比较（数据获取日期见表2）。结果表明（图8），本研究利用COSI-corr 提取的冰川运动速度在总体空间分布上与其他数据产品具有较好的一致性。但由于数据源、数据时间间隔以及提取方法的不同，导致各结果之间存在差异。在冰川中上游（6～18 km）本研究结果与ITS_LIVE 数据最为接近，平均运动速度比ITS_LIVE 数据高1.07%，但在靠近冰川末端处（0～3 km）两者结果相差较大，可能是因为ITS_LIVE 运动速度是一年内多期数据融合的结果；Sentinel-1 冰川运动速度数据虽然整体较低，但其不确定性为±0.08 m·d-1，在误差允许范围内，与本研究结果一致。GoLIVE 数据在冰川末段区域（0～8 km）与本研究的提取结果较为接近，但在冰川消融区中上游的差异较大，可能是由于选用间隔32天的GoLIVE 数据取平均，导致数据误差叠加，使最终结果产生较大差异。花杆实测数据在冰川末端区域（0～8 km）与本研究结果较为接近，但由于实测数据的获取时间为夏季，较高的温度和冰川融水导致冰的变形和滑动增大，从而使得末端的一些实测数据高于本研究结果。

图8 不同方法的冰川运动速度比较Fig.8 Comparison of glacier flow velocity by different methods

表2 研究所使用的主要数据Table 2 List of the data used in this study

4.2 冰川运动速度季节变化特征影响因素

科其喀尔巴西冰川运动速度呈现显著的季节变化，该变化趋势与气温和降水变化趋势一致（图9）。冰川运动是冰的变形、底部滑动和底部岩屑层运动等共同作用的结果，与冰川物质平衡、冰温、冰川的水力特征等密切相关［40-41］。冰川运动速度的季节性变化通常与表面融化引起的底部润滑有关［42］，冰川在暖季虽会因消融变薄，但暖季（5—9 月）冰温高，冰川表面消融强烈，大量冰川融水通过冰川裂隙流入冰川底部，减少冰川底部摩擦力，进而冰川运动速度增大［43-44］。降雨和融雪产生的水也可能与地表融水叠加，导致更大的融水输入，促使冰川底部滑动，从而影响冰川速度［45-46］。此外，在底部滑动开始时，冰川与底部基岩之间的摩擦产生了大量的热量，为冰川底部润滑提供了额外的融水［47］，进一步削弱了冰川的抗滑力，导致冰川速度增加。

图9 2014—2021年研究区月均气温及降水Fig.9 Monthly mean air temperature and precipitation in the study area from 2014 to 2021

4.3 冰川运动速度年际变化特征影响因素

科其喀尔巴西冰川年际运动速度呈现微弱下降趋势（图7），该现象可归因于冰川厚度减薄。近几十年来全球气候变暖导致冰川呈现负物质平衡状态，进而造成冰川运动呈现减速趋势。联合国政府间气候变化专门委员会（IPCC）第六次评估报告显示，近几十年全球平均升温速度为（0.2±0.1）℃·（10a）-1，山地冰川区域升温速率达到（0.3±0.2）℃·（10a）-1［48］，而气温是控制冰川物质损失的主要因素。Brun 等［49］对亚洲高山区冰川物质平衡的研究结果表明，2000—2016 年间天山山区年平均物质平衡为（-0.28±0.20）m w.e·a-1，累计物质平衡始终为负，这与本研究中冰川年际运动特征呈减小趋势的结论相一致。此外，液态降水也是加速冰川融化的重要驱动因素，当携带热量的雨滴落在冰体时，雨滴可能会侵蚀冰体［42，50］，导致冰川进一步消融。气温和液态降水是导致冰川物质平衡改变的重要因素，当它们共同促进冰川物质亏损时，就会使得冰川运动速度减慢。