1960—2020年高加索山地冰川变化研究

2022-09-14侯姗姗王宁练李志杰

冰川冻土 2022年3期

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

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

0 引言

冰川对气候变化极为敏感。随着气候加速变暖，近几十年来全球大多数地区的冰川均处于持续退缩状态，且21世纪以来的退缩幅度更加显著［1-2］。当前的全球升温已打破了冰川与气候环境的平衡，因此即便升温停止，全球冰川也将进一步退缩［3］。冰川作为全球各地的“水塔”［4-5］，对于区域生态平衡、社会经济发展有着巨大影响［6］。高加索山脉是全球山地冰川集中分布的地区之一，冰川融水是下游地区农业灌溉、工业发电的重要水源［7］，也是诱发地缘冲突的潜在因素之一。此外，由于气候变暖导致冰川稳定性减弱，该地区冰川灾害发生频率显著增加［8-9］，严重影响当地社会稳定发展。

18世纪至19世纪末，早期探险家就高加索地区冰川开展考察，留下了一些观察资料［10-11］。20世纪初期，科学家们致力于冰川编目工作［12-13］，例如Podozerskiy等［10］依据1881—1910年地形图完成了高加索地区的首套冰川编目（Glaciers of the Caucasus Mountain Range，PGI；1 329条，1 967.4 km2）。此后，苏联科学家依据1950—1960年航空照片对高加索冰川编目进行了更新（Catalog of Glaciers of the USSR，UGI；2 080条，1 427.12 km2）［13］；Gobejishvili等［14］和Tielidze等［15］则利用同样的数据源对UGI进行了校正，获取了格鲁吉亚高加索（786条，563.1 km2）和高加索整体（2 349条，1 674.9 km2）的冰川边界。21世纪初期，基于卫星遥感数据的冰川编目方法日渐成熟。Khromova等［16］利用1999—2004年的Landsat ETM+和ASTER影像解译了2000年时高加索山脉的1 706条，覆盖面积1 174.52 km2，但其冰川编目并不完整［17］；Tielidze等［18］基于1999—2002年的Landsat影像解译得到高加索2 186条冰川［（1 381.5±58.2）km2］，但其用于解译的遥感影像存在大量的季节性积雪、山体阴影以及云雾（图1）；Shahgedanova等［19］利用1999—2001年的ASTER和Landsat影像解译提取了2000年高加索山脉中部和西部共计407.3 km2的冰川面积。此外，Lur’e等［20］的研究表明1895—2011年北高加索山地冰川面积退缩速率为0.45%·a-1；Shahgedanova等［19］的研究结果表明1987—2010年包含高加索山中西段的冰川面积萎缩率为0.20%·a-1；Tielidze等［15］对高加索冰川面积变化的研究结果显示，1960—2014年该地区冰川面积的整体萎缩率为0.53%·a-1，2000—2020则增加为1.16%·a-1［18］，但其使用的冰川边界却存在较多的错误与误差（图1）。总体上，当前对于高加索冰川分布与变化的认识尚存在争议，冰川变化的时空格局也有待深入。

图1 Tielidze等［18］解译冰川边界时使用的部分Landsat ETM+/TM影像（a）～（d）及其编目中存在的部分错误（e）～（f）Fig.1 Some of the Landsat ETM+/TM images used in the interpretation of glacier boundaries by Tielidze et al［18］（a）～（d）and some inaccurate glacier boundaries（e）～（f）

高加索是“一带一路”沿线山地水塔的重要组成部分，探究区域气候对于冰川变化时空格局的影响对预测冰川未来演化趋势至关重要［8-9，21-24］。高加索山脉的Djankuat和Garabashi冰川均是全球42条参照冰川之一，物质平衡监测时长均超过30年［25］，为该地区冰川变化研究提供了难以取代的观测资料。本文基于美国Landsat系列影像以及欧空局Sentinel-1影像生成的干涉图，采用目视解译的方法获取高加索山脉2000年和2020年冰川边界，并结合全球陆地冰空间测量计划（Global Land Ice Measurements from Space，GLIMS）中高加索山脉1960年和1986年的冰川编目，分析了近60年来高加索地区冰川变化的时空格局。在此基础上，结合世界冰川监测服务处（World Glacier Monitoring Service，WGMS）提供的Djankuat和Garabashi两条冰川的物质平衡监测资料及英国CRU（Climate Research Unit）气温和降水资料，探究了近60年来高加索地区气候变化对冰川变化的影响。

1 研究区概况

高加索山位于欧亚大陆交界，黑海与里海之间（40°～44° N，40°～49° E）。山脉海拔由西向东先升高后降低，降水量则表现为由东向西逐渐增加，最西端地区的降水量（3 200 mm）是东部地区（1 000 mm）的三到四倍［26］。受气候和地形影响，在高山区降水形式以降雪为主，尤其是受来自黑海潮湿气团影响的高加索西部迎风坡［27］，山脉北坡相较于南坡更加干燥寒冷，年降水量减少500～1 000 mm、夏季气温低1～2℃［28］，山脉内的冰川集中分布在山脉中段，且北坡的冰川规模远大于南坡［29］。

图2 高加索山脉地理位置（冰川资料来源于GLIMS所提供的1960年高加索山地冰川边界；CRU资料来源于英国东英吉利大学的气候研究中心）Fig.2 Geographic location of the Caucasus Mountains（glacier information from GLIMS glacier database；CRU data from University of East Anglia，UK）

2 数据与方法

2.1 数据来源

2.1.1 遥感影像

本研究采用43景Landsat光学遥感影像（附表1），6景IW成像模式下的Sentinel-1（S-1）单视复数影像（SLC，附表2），对高加索山脉内2000年和2020年的冰川边界进行解译。其中，Landsat卫星影像均为L1T级产品，来源于美国地质调查局（United States Geological Survey，USGS；https：//glovis.usgs.gov），S-1影像来源于欧洲航天局（European Space Agency，ESA；https：//scihub.esa.int）。

附表1 Landsat影像列表Attached table 1 List of Landsat data used in this study

附表2 S-1影像数据列表Attached table 2 List of S-1 data used in this study

为减少季节性积雪、云雾、山地阴影等的影响，缩小冰川边界的获取时间跨度，本研究尽可能选取2000年和2020年消融季末期的影像。当特定年份影像质量无法满足解译需求时，采用相邻年份的影像进行补充（图3）。此外，在冰川边界的解译过程中，还参考了高分辨率的Google Earth影像。

图3 2000年和2020年冰川编目所使用的Landsat遥感影像的空间分布、数量及覆盖冰川面积Fig.3 Spatial distribution，number and area of glaciers covered by Landsat images used for glacier inventory in 2000 and 2020

2.1.2 GLIMS编目

GLIMS是一个全球性的冰川编目数据库，本研究选用了其中高加索山脉1960年和1986年的冰川编目（http：//www.glims.org/）。其中1960年冰川编目是根据1950—1960年1∶50 000航测地形图解译而成，1986年冰川编目则是基于1985—1987年Landsat TM影像解译完成［29］。

2.1.3 数字高程数据集

本研究选用美国国家航空航天局（National Aeronautics and Space Administration，NASA）于2020年2月发布的NASADEM提取冰川属性信息，该数据可通过IPDACC平台下载（https：//lpdaac.usgs.gov/products/nasadem_hgtv001/）。NASADEM是基于SRTM DEM，通过高程控制、空洞填充等手段处理改进而生成的全球DEM数据，空间分辨率为30 m。相比于原始的SRTM DEM，NASADEM的精度更高，数据空洞也获得了有效填充［30］。

2.1.4 冰川物质平衡观测资料

Djankuat和Garabashi冰川拥有长时期的物质平衡观测资料［25］，是全球42条参照冰川之一，能够有效代表高加索地区的冰川物质平衡变化。监测冰川物质平衡资料来源于WGMS冰川物质平衡通报（https：//wgms.ch/products_ref_glaciers/），其中Djankuat冰川观测时段为1967—2020年，Garabashi冰川为1983—2020年［25］。

2.1.5 气象数据

本研究的气温和降水资料来源于东英吉利大学的气候研究中心发布的CRU V4.05气候再分析资料（http：//www.cru.uea.ac.uk/）。CRU数据集是目前全球使用较为广泛的近地面气候数据集之一，包含气温、降水、湿度、风速和蒸发等气象要素，为1901—2022年间逐月格点数据（0.5°×0.5°）。本研究选取了其中1961—2020年间高加索地区的平均气温与降水量月值数据。

2.2 方法

2.2.1 冰川边界提取

本研究利用ENVI软件中的Gram-Schmidt光谱锐化方法［31-32］，对Landsat OLI/ETM+影像进行全色波段融合，将其空间分辨率由30 m提升为15 m，并基于Landsat假彩色影像中冰川区与非冰川区色彩、纹理差异解译冰川边界［图4（a）］。对于高加索山大量存在的表碛覆盖型冰川，本文借助S-1的相干系数图辅助判识冰川边界［图4（b）］。S-1干涉可以捕捉到地表微小的形变，其核心是估算复相关系数的振幅，即相干性。相干性值的范围在1到0之间，地物形变越快，相干性越弱。由于冰川具有运动和消融特性，冰川区相较于周边的稳定基岩存在明显的变化特征，因此该方法已被广泛地应用于冰川边界解译中［33-34］。最后，本研究使用高分辨率的Google Earth影像对冰川边界进行检验［图4（c）］，并将冰川面积的最小阈值设置为0.01 km2。

图4 使用相干系数图像辅助解译表碛覆盖冰川示例Fig.4 Example of using coherence images to assist in the interpretation of debris-covered glacier：Landsat OLI image（bands 654 as RGB）（a）；debris-covered glacier shown on S-1 coherence image（b）；Google Earth（Maxar imagery，0.30 m）（c）

目视解译虽然耗时较长，但却是当前精度较高的冰川边界解译方法之一，可在很大程度上避免计算机提取时阴影、云、积雪等造成的错误，且适用于各类冰川边界的提取。本研究参照中国第二次编目中采用的方法，对冰川分冰岭、冰川平均朝向、中值高度等进行提取［35］。此外，通过冰川边界叠加DEM，统计不同海拔高度带内的冰川面积及其变化；通过冰川的平均朝向统计不同坡向范围内的冰川面积及其变化。

2.2.2 冰川面积不确定性评估

为了确保观测到的冰川变化为实际变化，而不是由解译误差引起，必须准确评估冰川边界解译的不确定性。目前，对于冰川解译不确定性的评估主要通过缓冲区的方法进行［36-38］，但这一方法更适用于自动提取的冰川边界［35］。在目视解译的过程中，即使是由同一人对同一区域的冰川进行多次解译的结果都会有所不同，且表碛区冰川边界和裸冰区冰川边界误差往往存在差异，因此采用目视解译的冰川边界误差很大程度上取决于解译者的经验和影像分辨率。

鉴于以上原因，本研究结合Paul等［39］提出的对冰川多次数字化及Guo等［35］提出的以Google Earth高分辨率影像得出的冰川边界为真值的两种方法，分别计算裸冰区和表碛区冰川边界与“真值”之间的平均距离，并将结果作为冰川边界解译误差，从而计算得到冰川面积变化误差（图5）。本研究使用的1960年和1986年冰川编目的误差则沿用了数据提供者的计算结果［29］。在对冰川面积变化误差的计算上，本文综合考虑前后两期冰川边界的面积误差，基于误差传递公式来进行计算［40］。

图5 基于高分辨率遥感影像的冰川面积不确定性评估：基于Google EarthMaxar影像手工数字化的冰川边界（a）；基于Landsat OLI假彩色影像多次手工数字化的冰川边界示例（b）Fig.5 Glacier area uncertainty assessment based on high-resolution imagery：manual digitization of glacier boundaries based on Google EarthMaxar imagery（a）；example of glacier boundaries based on multiple manual digitization of Landsat OLI composite false colour imagery（b）

经测算，Landsat TM和OLI影像中冰川裸冰区和表碛区边界的解译误差分别为20 m、21 m和11 m、9 m。由于冰川内部分冰岭不会导致区域冰川面积的误差，因此忽略了分冰岭的影响。本研究采用公式（1）对冰川边界解译误差进行计算［35］，采用公式（2）对不同阶段冰川面积变化误差进行计算，在计算时仅考虑冰川边界发生变化的部分［40］：

式中：Δσ指冰川面积的不确定性；σpc和σpd分别指Landsat影像上裸冰区和表碛区的冰川边界解译误差；Lc和Ld分别指裸冰区和表碛区边界的长度；σAB指冰川在某一时段内冰川面积变化的不确定性；σA、σB分别指A时期和B时期冰川面积的不确定性。经计算，2000年和2020年高加索山地冰川边界解译的不确定性分别为±138.37 km2和±66.44 km2，占同期冰川面积的10.29%和6.11%。1960—2020年高加索山地冰川面积变化的不确定性为±98.66 km2（±5.89%）。

3 结果

3.1 1960—2020年高加索冰川面积变化

2020年高加索山共分布有山地冰川1 912条，总面积（1 087.36±66.44）km2。1960—2020年高加索山脉的冰川面积减少了（587.36±98.66）km2，面积萎缩率为（35.07±5.89）%，年均萎缩率为（0.58±0.10）%。从不同阶段来看，1960—1986年，冰川面积共减少（193.51±85.54）km2，年均萎缩率为（0.44±0.20）%·a-1；1986—2000年，冰川面积共减少（136.20±159.40）km2，年均萎缩率增加为（0.66±0.77）%·a-1；2000—2020年，冰川面积共减少（257.66±82.18）km2，年均萎缩率进一步增加为（0.96±0.31）%·a-1（表1）。就空间变化格局而言（图6），1960—2020年高加索山东段冰川面积萎缩率明显高于中西段，这体现了冰川规模对冰川变化存在直接影响，小规模冰川相较于大规模冰川的面积退缩更为显著。

表1 1960—2020年高加索山地冰川面积变化特征Table 1 The Caucasus Mountains glacier area changes during 1960—2020

图6 0.25°×0.25°格网中显示的1960年冰川面积分布和1960—2020年冰川面积变化（以冰川中心点所在格网为依据统计）Fig.6 Distribution of glacier area in 1960 and change in glacier area from 1960 to 2020 shown in the 0.25°×0.25°grid（statistics based on the grid where the glacier centroid is located）

统计表明，高加索山地冰川在数量构成上以面积≤1 km2小规模冰川为主，在面积构成上则以＞1 km2的较大规模冰川为主，其中，2～5 km2区间的冰川面积占比最大（图7）。1960—2020年间大于20 km2的冰川面积减少最多（156.99 km2）；其次是2～5 km2的冰川（112.79 km2）；1～2 km2，0.1～0.5 km2，5～10 km2和0.5～1 km2等区间的冰川面积减少量基本相当（64.05～83.16 km2）；10～20 km2区间的冰川面积减少14.09 km2；而小于0.1 km2的冰川面积基本保持不变。在高加索山冰川持续退缩的同时，规模≤0.1 km2的冰川数量却持续增加，面积占比也在1960—1986年和2000—2020年分别增加了1.48 km2、0.09 km2，这主要是更大规模冰川消融分裂所致；10～20 km2区间内的冰川面积和数量1986—2000年增加同样是更大规模冰川分裂所致。

图7 1960—2020年不同时期各规模等级冰川的面积和数量值Fig.7 Area and number of glaciers of different sizes for different periods during 1960—2020

以200 m高差为间距，统计了高加索冰川在不同高程带内的分布与变化。结果表明，高加索山地冰川面积随海拔的变化呈正态分布（图8）。2020年高加索冰川分布在海拔1 381～5 627 m之间，其中超过90%（1 026.88 km2）的冰川面积集中分布在海拔2 400～4 400 m之间；随着海拔的增加，冰川面积萎缩速率逐渐减小，在海拔2 400 m以下时冰川面积萎缩速率为1.03%·a-1，但当海拔上升到4 600 m以上时，冰川面积萎缩速率仅为0.02%·a-1。1960—2020年海拔2 600～3 800 m内的冰川面积共计减少514.01 km2，占总萎缩量的87.51%。随着低海拔地区冰川面积的退缩，高加索冰川的平均中值面积高度也由3 295 m上升至3 320 m。

图8 1960—2020年不同时期冰川面积随高程带的分布及其变化情况Fig.8 Distribution of glacier area with elevation and its variation in different periods from 1960 to 2020

从冰川平均坡向分布来看（图9），高加索山地冰川呈现显著的北多南少的特点，其中以正北朝向的冰川面积和冰川数量比重最大，接近三分之一。1960—2020年间各朝向冰川面积均呈减少趋势，N、NE朝向的冰川面积萎缩率明显大于其他朝向（143.56 km2、122.98 km2），分别占总萎缩量的24.44%和20.94%，其次是SE朝向（-94.16 km2）；NW和SW朝向的面积萎缩量分别为77.71 km2和75.77 km2，E和SE朝向的冰川面积萎缩率分别为33.59 km2和32.83 km2，W朝向冰川面积萎缩量最少，仅为6.76 km2；从各朝向面积萎缩率看，SE和SW朝向的冰川面积退缩幅度最大，分别为1%·a-1、0.99%·a-1，W朝向幅度最小，仅为0.13%·a-1，其余朝向则介于（0.47～0.58）%·a-1之间。各朝向冰川面积萎缩率的不同主要与区域气候有关，西风水汽对W朝向的冰川补给充足，导致其面积萎缩率较低；南朝向冰川接收的太阳辐射偏多，面积萎缩率普遍偏大。

图9 1960—2020年不同坡向的冰川面积百分比（a）、数量百分比（b）Fig.9 Proportion of glacier aspect by percentage of area（a）and number（b）during 1960—2020

3.2 参照冰川物质平衡变化

Djankuat（G042761E43193N）冰川位于高加索山北坡，1967—2020年面积萎缩率为0.44%·a-1，Garabashi（G042466E43313N）冰川位于高加索最高峰——厄尔布鲁士山峰南坡，1983—2020年面积萎缩率为0.36%·a-1。Djankuat和Garabashi冰川物质平衡年际变化幅度较大，且以负平衡年份为主（图10）。1967—2020年，Djankuat冰川年物质平衡介于-2 010～1 540 mm w.e.，平均为-276 mm w.e.·a-1，负平衡年份与正平衡年份之比为3：2。1983—2020年，Garabashi冰川年物质平衡介于-1 510～410 mm w.e.，平均为-393 mm w.e.·a-1，负平衡年份与正平衡年份之比为2：1。

图10 Djankuat和Garabashi冰川的物质平衡和累积物质平衡变化过程Fig.10 Change of mass balance and cumulative mass balance of Djankuat（a）and Garabashi（b）Glacier during 1967—2020 and 1983—2020

如图10所示，Djankuat和Garabashi冰川的累积物质平衡总体呈持续下降趋势。1967—2020年Djankuat冰川累积物质平衡达到-13 790 mm w.e.，其中1967—2000年累积物质平衡下降速率为56 mm w.e.·a-1，2000—2020年剧增为517 mm w.e.·a-1，达到了1968—2000年的约9倍；1983—2020年Garabashi冰川的累积物质平衡为-14 148 mm w.e.，其中1983—2000年累积物质平衡下降速率为97 mm w.e.·a-1，2000—2020年剧增为571 mm w.e.·a-1。此外，Djankuat和Garabashi冰川的物质平衡线高度（Equilibrium Line Atitude，ELA）也在不断上升，积累区比率（Accumulation Area Ratio，AAR）在不断减小。

4 讨论

4.1 与其他研究对比

4.1.1 与其他编目对比

目前能够完全覆盖高加索山脉的冰川编目主要有Randolph Glacier Inventory（RGI）6.0、GLIMS以及Tielidze等［18］进行的编目。其中，RGI 6.0中高加索区域的冰川编目获取时间集中于2000年，GLIMS提供的最新的高加索冰川编目集中在2014年，Tielidze等公布的冰川编目则分别是2000年和2020年。考虑到不同编目之间的可对比性，本文对RGI 6.0、Tielidze等发布2000年冰川编目与本研究中解译得到的2000年冰川编目进行对比（表2）。

表2 与高加索不同冰川编目结果对比Table 2 Comparison with the different glacier inventory in the Caucasus

通过表2可知，3种冰川编目在冰川数量和面积上都存在一定差异，这些差异的原因主要与数据来源、解译标准以及解译方法的不同有关。首先，RGI 6.0是基于Landsat影像采用自动提取加后期人工修正的方法获得，而Tielidze等的编目和本研究的冰川边界都是直接通目视解译获得；其次，RGI 6.0在高加索山东段和西段部分地区存在数据缺失［17］，而Tielide等用于解译冰川边界的Landsat遥感影像质量不佳，存在大量积雪、云雾和山体阴影（图1），导致解译得到的冰川边界与实际不符；最后，编目之间的冰川数量和面积的差异在一定程度上也与最小冰川面积的阈值不同有关。

4.2.2 与不同时段研究结果对比

Koromova等［41］的研究结果显示1910—1965年高加索山地冰川面积萎缩率为0.45%·a-1，Tielidze等［15］估算的1960—2014年高加索山地冰川面积萎缩速率为0.53%·a-1，本研究结果表明1960—2020年高加索冰川面积萎缩率达到（0.58±0.01）%·a-1。因此，近百年来高加索山地冰川发生了持续退缩，且退缩速率在不断增加。

4.2.3 与其他山系对比

与其他山系的对比表明（表3），近几十年以来高加索山冰川面积萎缩率明显低于欧洲阿尔卑斯山和亚洲阿尔泰山，与加拿大西部的海岸山脉基本持平，但明显高于亚洲的天山、祁连山、帕米尔高原、喀喇昆仑、西昆仑和喜马拉雅山等地区。冰川退缩主要受气候影响，阿尔卑斯山气温较高、阿尔泰山冰川发育规模则普遍偏小，因此冰川面积萎缩率偏高；天山、祁连山、喜马拉雅、念青唐古拉山的平均海拔超过4 000 m，气候高寒，因此冰川退缩相对缓慢［42］；特别是帕米尔高原、喀喇昆仑山和西昆仑山冰川发育规模巨大，近几十年来总体保持稳定［36］。因此高加索山冰川面积萎缩速率低于阿尔卑斯山和阿尔泰山，但普遍高于亚洲高山区。

表3 与其他山系冰川面积变化的对比Table 3 Comparison of glacier area change with other mountains

4.3 气候对冰川变化的影响

1961—2020年高加索山脉平均气温以0.24℃·（10a）-1的速度持续升高，降水量变化则基本保持稳定。高加索山气温和降水变化可以2000年为界分为两个阶段（图11），1961—2000年，气温和降水波动平缓，无显著变化；2000—2020年，气温以0.58℃·（10a）-1的速率快速上升，降水则-55 mm·（10a）-1的速率快速下降。

图11 1961—2020年高加索山地冰川区年均气温（a）和年均降水量变化（b）Fig.11 Variations of annual mean air temperature（a）and annual precipitation（b）in glaciated areas of the Caucasus Mountains from1961 to 2020

气温和降水变化直接影响冰川物质的消融和积累［52］，但山地冰川的长期变化主要受气温控制［53］，因此本文重点分析了高加索地区气温的长期变化对冰川变化的影响。1961—2020年，高加索地区气温变化存在显著的空间差异（图12），总体上升温速率由西向东逐渐增加，其中在厄尔布鲁士山附近的升温速率最低［0.21℃·（10a）-1］，在山脉东段升温速率最高，可达到0.27℃·（10a）-1。近几十年来，高加索山冰川变化的空间格局与气温变化的空间格局高度一致（图12），主要表现为在升温速率较低的中段地区，冰川面积萎缩率最小，而升温率最为显著的东段地区冰川面积萎缩率最大。

图12 1961—2020年高加索山地冰川面积萎缩率与气温之间关系Fig.12 Relationship between the glacier area shrinkage and air temperature in the Caucasus Mountains during 1961—2020

冰川物质平衡受到区域水热条件的影响，对气候变化的响应迅速［54-55］。Djankuat与Garabashi冰川所在地区气温均呈上升趋势，在2000年以前气温缓慢上升，两条冰川的正负物质平衡年份交替出现，冰川总体呈微弱的物质亏损状态；2000年之后，气温上升速率显著加快，两条监测冰川的物质负平衡不断增强，表明冰川消融在不断加速（图13）。气温的持续快速上升已打破了监测冰川与气候环境的平衡，因此即使高加索地区升温停止，冰川也会持续处于物质亏损状态［56］。