1987—2016年雀儿山冰川变化及其对气候变化的响应

欧健滨，许刘兵，蒲 焘

（1.华南师范大学地理科学学院，广东广州510631；2.中国科学院西北生态环境资源研究院冰冻圈科学国家重点实验室玉龙雪山冰川与环境观测研究站，甘肃兰州730000）

0 引言

作为地球冰冻圈的重要组成部分，冰川是气候变化最为敏感的指示器，它已成为了解气候变化过程和趋势的关键指标［1］。过去数十年间，学者们借助多种手段对青藏高原及其周边山地现代冰川的状态及其变化过程开展了大量工作。总体来看，青藏高原地区现代冰川呈现广泛退缩趋势，但高原上不同地区、不同类型冰川的退缩速率有较大差异，而帕米尔高原和喀喇昆仑山等地冰川甚至出现了前进现象［2-3］。关于青藏高原冰川差异性变化的原因，有研究认为可能与西风环流的加强和印度季风的减弱［2］、大气环流季节性变化作用［4］以及区域气候的干湿程度［5］等因素有关。此外，青藏高原内部和西部的大陆型和极大陆型冰川，与高原东南部的海洋型温冰川对气温与降水变化的敏感性有所不同［6-7］。所以，要全面了解青藏高原现代冰川对气候变化的响应规律，还需选取不同气候区的典型冰川开展系统监测研究。

实地观测是研究现代冰川变化最为精确的方法，其观测数据是理解现代冰川过程的重要基础［8］。但对于青藏高原广袤的冰川覆盖范围而言，现有实测研究点的数量仍非常有限［2］。近年来，国内外众多学者基于遥感数据对青藏高原冰川变化及其气候响应［9-11］、高原内陆冰川消融差异［12-13］以及冰川编目［14-15］等方面开展了大量的研究工作。但这些研究主要聚焦于喜马拉雅山、天山等地的大陆型或极大陆型冰川分布地区［16］，对藏东南地区海洋型温冰川的研究还不够系统。同时，在冰川遥感中，比值指数、NDSI（normalized difference snow index，归一化雪盖指数）及监督与非监督分类等基于光谱特征的传统分类方法，容易导致水体的混分以及阴影区域冰川的漏分［17-18］，如何提高遥感信息提取冰川的精度，也是相关研究的重点［19］。

本文以青藏高原东南部地区的雀儿山为例，尝试将地形阴影模拟结果引入冰川分类过程，通过光谱与阴影模拟信息的综合，提取了雀儿山地区1987—2016年的多年冰川边界，并结合1979—2016年的气象数据，对研究区的冰川变化及其气候响应机制进行探讨。这对认识青藏高原东南部地区海洋型冰川变化及其与气候变化之间的关系、了解不同类型冰川对气候变化的响应方式具有重要意义。

1 研究区概况

雀儿山（98°47′～99°35′E，31°34′～32°6′N）位于青藏高原东南缘，呈西北—东南走向（图1），是青藏高原较典型的海洋型冰川分布区。该区冰川平衡线高度为5 100～5 200 m，冰川区降水可达1 000 mm以上［20-21］。据第二次冰川编目数据，研究区内共有冰川113条，总面积为103.9 km2，主要分布在海拔4 700～6 000 m范围内，其中大部分冰川面积小于1 km2（表1）［22］，附近的德格（海拔3 280 m）与甘孜（海拔3 380 m）气象站年平均气温分别为6.82℃与5.89℃，年平均降水量为638 mm与652 mm。

图1 研究区位置示意图Fig.1 Location of the study area

表1 研究区不同规模冰川的数量及面积Table 1 The glacier number and glacierized area of various sizes in the study area

2 数据与方法

2.1 数据来源

本文所采用的遥感影像以及数字高程数据（DEM）来自于美国地质调查局（http：//glovis.usgs.gov/），具体包括Landsat TM（Thematic Mapper，专题绘图仪）、Landsat 7 ETM+（Enhanced Thematic Mapper Plus，增强型专题绘图仪）和Landsat 8 OLI（Operational Land Imager，陆地成像仪）三种传感器的L1T级产品以及30 m空间分辨率的ASTER GDEM数据。为减少云层和积雪对冰川边界提取精度的影响，主要选择云量较低，且没有云层覆盖冰川的7—9月遥感影像，合计13景，摄于1987—2016年（表2）。本文采集了研究区附近的德格、甘孜、色达、石渠以及新龙五个气象站点1979—2016年的月均气象数据（站点位置见图1），来源于中国气象科学数据共享服务网站的中国地面气候资料数据集（http：//cdc.cma.gov.cn/）。

2.2 冰川提取方法

指数法［23］、比值法［24］、面向对象分类［25］等基于光谱信息特征的方法是目前主要采用的冰川计算机分类方法。这些方法在图像质量较好的情况下能相对精确地识别和提取冰川信息，但在阴影、云、雪等因素的干扰下，这些方法在分类过程中容易出现错分、漏分现象［17-18］。为解决这一问题，本文在比值法、波段阈值法等方法的基础上，尝试将基于DEM的山地阴影模拟结果纳入到遥感分类过程中，从而提高冰川信息提取的精度。

表2 遥感影像及其参数Table 2 The parameters of selected remote sensing images

2.2.1 基于比值法与波段阈值法的冰川初步提取

比值法操作相对简单，精度较高，因而在遥感冰川信息提取中，它是采用最多的一种方法［19］。在TM图像中，考虑到阴影区域冰川的光谱特征变化，常选用红波段和中红外波段进行比值计算来突出冰川信息［18］。本文选取TM 3和TM 5进行比值计算，以2.0为阈值来初步区分冰川与非冰川区域。但由于冰雪与水体信息在这两个波段之间的比例关系较为接近，因此在使用比值法提取冰川时，往往需要进行人工修正或结合其他信息排除水体信息的干扰。

在遥感冰川信息提取中，阈值分割法可以较为方便高效地区分冰川和某些特定地物［26］。上述比值法的冰川提取结果混有了部分云和水体信息，而部分处于阴影当中的冰雪区却被漏分。因此，本文基于2016年的OLI影像，采用样本点像元统计方法，对冰川、云、阴影和水体四部分进行样本点DN值统计分析（图2）。由于冰川在中红外与短波红外波段的DN值明显小于云，因此本文选择上述两个红外波段，通过阈值分割的方法，对云的混分信息进行排除。

2.2.2 基于山体阴影模拟的冰川提取方法

在阴影的干扰下，被阴影覆盖冰川的光谱辐射特征与水体十分接近，难以依靠计算机自动分类方法对两者进行区分。因此，学者们通常采用目视解译［14］、蓝波段阈值分割［15］或通过人工调试比值阈值［17-18］等方法来减少阴影的影响。但在本研究区，水体与阴影区域冰川在蓝波段的光谱特征非常接近（图2），且不同部位的最优波段比值阈值有较大差别，因此需要借助遥感上的其他特征进行阴影冰川的提取。

图2 不同波段样本点DN值统计Fig.2 The DN（digital number）curves of glacier，shadow，cloud and water

Li等［27］将DEM数据结合到遥感分类过程中，较好地解决了高原上山体阴影与冰湖水体的混分问题，这为解决阴影区域冰川与水体混分问题提供了参考。在遥感影像中，山体阴影与山体的坡度、坡向以及遥感影像成像时刻太阳高度角和方位角存在着定量关系［28-29］。因此，本文利用经配准的DEM以及遥感图像成像时刻太阳高度角与方位角数据，通过GIS中的表面阴影分析得到阴影二值图像。随后基于决策树分类方法，将阴影模拟结果、坡度图像以及剩余的阴影冰川-水体混分信息进行综合并通过目视解译对冰川边界进行修正，实现冰川的整体识别。提取流程与结果见图3～4。

2.3 冰川分冰岭提取方法

通过识别冰川复合体中的分冰岭，可以将其划分为多条独立冰川。由于分冰岭的实质是山脊线，因此多数研究主要通过提取冰川所在小流域的边界作为山脊线，随后通过相交操作实现对冰川复合体的分割。但该方法需要确定所要提取最小流域面积的参数，而该参数的准确性直接影响分冰岭提取的精度［15］。由于雀儿山地区的冰川复合体数量较少，为方便起见，本文基于全局分析法（global process）原理，利用GIS中的表面分析和栅格计算器等功能实现山脊线（分冰岭）的提取。具体过程为：①求出正、反地形的坡向变率，进一步得出无误差DEM的坡向变率；②基于原始DEM计算邻域为11×11的平均值数据层；③将原始DEM减去平均值数据层，得到正负地形分布区域；④求出正地形下无误差DEM坡向变率大于70的区域，即为山脊线区域；⑤将生成的栅格结果转换为矢量弧段，进行平滑与人工修正，得到分冰岭提取结果。

图3 冰川提取流程Fig.3 Glacier image extraction process

图4 冰川提取结果示例Fig.4 Examples of glacier image extraction results：TM 7，6，5 false color remote sensing image（a）and identified image（b）

2.4 不确定性分析

本文在冰川自动提取的初步结果上，进一步采用目视解译对冰川边界进行修整。冰川与非冰川之间存在着包含两类信息的混合像元，需要通过统计冰川轮廓线长度与像元大小来计算冰川面积信息的不确定性［15，30］。

式中：EA为冰川面积信息提取的误差；LC为提取冰川的边界长度；EPc为边界的定位精度；E为某一时段内冰川面积变化的不确定性；EA1与EA2分别为该时段起始和终止时的冰川面积提取误差。计算可知，本文各个时段年均冰川面积变化不确定性为±5.85%～±7.71%。

3 结果与讨论

3.1 研究区冰川面积及其变化

图5 与表3为研究区冰川面积变化情况及其多项式拟合结果。近30 a来，雀儿山地区冰川面积从1987年 的（103.63±29.31）km2减 少 到2016年 的（52.82±13.84）km²，共退缩了（48.97±25.29）km2（47.25%），远超其他地区的平均水平（18%）［22］。根据冰川面积变化率的不同，可以将研究区1987—2016年间冰川面积的变化分为三个阶段：①1987—1993年的快速退缩阶段，此时冰川面积减少了（21.72±29.14）km2，变化率达（-3.49±4.69）%·a-1；快速退缩，并在2009—2011年呈现突变下降趋势，该时段冰川年变化率达（-6.67±6.46）%·a-1，但在2011年后冰川面积退缩趋势放缓，年变化率减小为（-2.19±5.73）%·a-1。

图5 1987—2016年雀儿山冰川面积变化Fig.5 Variation of the glacierized area and its fitting curve in the Que’er Mountains，1987—2016

表3 1987—2016年部分年份雀儿山冰川面积及其面积变化量和变化率Table 3 Glacierized area at a certain year，and its changing amount and changing ratio，in the Que’er Mountains，1987—2016

3.2 不同规模冰川的变化特征

如图6所示，面积较小的冰川变化幅度较大，其中面积＜0.5 km²且面积损失大于50%的冰川占雀儿山冰川总数量的53.60%，面积＜1 km²的冰川面积共减少25.91 km²，占冰川损失面积的51.10%；而面积≥5 km2的冰川退缩幅度为8.81km²，占冰川损失面积的17.38%。因此，无论是相对变化面积还是绝对变化面积，＜1 km²的小型冰川都占绝对主导地位。从1987—2016年雀儿山地区冰川面积减少比例来看（图7），由于位于研究区中部的新路海地区冰川面积最大，因此其面积减少比例最小，为19.4%，其次为雀儿山西北部，为24.8%。

图6 1987—2016年雀儿山不同规模冰川的面积变化Fig.6 Changed glacierized area changing with glacierized area size（a），percentage of changed glacierized area changing with glacierized area size（b），in the Que’er Mountains from 1987 to 2016

图7 1987—2016年雀儿山冰川面积减少比例Fig.7 Decreasing proportion of glacierized area in the Que’er Mountains from 1987 to 2016

3.3 冰川随坡向分布特征及其变化

为统计研究区冰川在不同坡向的分布及其变化特征，本文以22.5°方位角为间隔对冰川朝向进行分类。如图8所示，在过去30 a间，雀儿山地区北向、南向、东向和西向的冰川年变化率分别为-1.71%·a-1、-2.16%·a-1、-1.74%·a-1和-2.14%·a-1，其中以朝向202.5°～225°之间的冰川年变化率最大，达到了-2.26%·a-1。这说明雀儿山地区的冰川退缩速度南向大于北向，西向大于东向，与冰川面积分布特征相反。这种冰川变化差异可能与北半球不同坡向冰川所接受的太阳辐射量、降水的差异［31］，以及与冰川规模等因素有关。

图8 冰川变化坡向特征Fig.8 Rose graphs showing glacierized areas of various aspects in 1987 and 2016（a），and rose graphs showing annual mean shrinkage ratio of glacierized area from 1987 to 2016（b）

3.4 冰川垂直分布特征及其变化

为研究雀儿山地区冰川变化垂直特征，本文以100 m为间距，将1987年与2016年冰川分布的高程范围划分为11个梯度（图9）。总体而言，研究区各时期冰川面积随海拔升高呈现先增加后减小的趋势，且在5 101～5 200 m处达到最大。从冰川变化率来看：冰川变化率先随海拔的升高而增大，在5 101～5 200 m处年变化率达到最大，为-2.52%·a-1，之后随海拔进一步上升，变化率逐渐减小，在海拔5 400 m以上年变化率稳定在-0.47%·a-1～-0.66%·a-1区间内。

图9 不同海拔带冰川的面积变化Fig.9 Profiles of per 100 m elevation glacierized area in 1987 and 2016，and their mean shrinkage ratio

另外，1987—2016年间，雀儿山冰川末端平均高程从5 181.4 m上升到5 244.7 m，年变化率为-2.26 m·a-1，是青藏高原主要山系中冰面高程上升速度最快的地区之一［16］。如图10所示，1987—2016年间，新路海冰川末端在水平方向上退缩了近500 m，末端平均高程则从4 343 m上升到4 688 m，年变化率高达-11.33 m·a-1。

图10 1987—2016年新路海冰川末端变化Fig.10 The terminus peripheries of the Xinluhai Glacier in 1987，1993 and 2016

值得注意的是，海拔越低的地区气温越高，其冰川退缩速度理论上也应随之增大。但在雀儿山地区，冰川退缩速率最大的区域并非海拔最低区域，这可能与冰川所处的坡向以及地形遮蔽程度有关。在雀儿山海拔4 900 m以下地区，朝向偏西的冰川占该海拔范围全部冰川比例大于海拔4 901～5 000 m地区（图11）。由于该地区偏西方向的冰川变化率大于偏东方向，由此推断，冰川的朝向差异对雀儿山地区不同海拔高度的冰川退缩速率有一定的影响。此外，本文将1987年冰川提取结果与阴影模拟图像与叠加分析，发现在夏季接近正午时刻，海拔4 900 m以下处于阴影区域的冰川像元数占冰川总像元数的12.18%，而在海拔4 901～5 000 m处为11.87%，即低海拔区域的冰川被阴影遮蔽的程度更为显著。因此，阴影区范围的差异也可能是影响不同海拔地区冰川退缩速率的一个因素。

图11 不同朝向冰川面积分布情况Fig.11 Rose graphs of glacierized area distribution ratio，below 4 900 m a.s.l.and between 4 901 m and 5 000 m a.s.l.

4 讨论

4.1 气候变化与冰川变化对比分析

冰川变化主要受太阳辐射、降水、地形、朝向、及冰川规模等因素的影响，但相对而言，夏季气温和冬季降水量对冰川变化的影响更为直接［20］。同时，冰川的变化通常滞后于气候变化。根据丁永建［32］的研究，长度小于5 km的冰川对气候变化的响应时间约为2～3 a，而长度大于5 km的冰川其响应时间大约为8～9 a。由于研究区中的冰川规模较小，故本文以3 a的滞后期来分析冰川对气候变化的响应。本文的显著性检验系数包括R2值和回归P值，其中P≤0.01为非常显著，0.01＜P≤0.05为显著，P＞0.05为不显著。

4.1.1 研究区年降水量变化特征

本文以研究区附近甘孜、色达、德格等六个气象站点（图1）上一年度9月至本年度8月的降水量总和为年均降水量，对研究区的降水变化趋势进行分析（图12）。结果显示，年均降水与3年滑动年均降水的线性拟合检验系数分别为R2=0.03，P=0.90以及R2=0.08，P=0.12，两者均没有通过显著性检验，说明雀儿山地区降水变化呈现不显著的变化趋势，此外，3年滑动年均降水数据也未与冰川面积变化存在明显的相关性。总体而言，研究区年均降水波动较大，但无明显变化趋势，从3年滑动年均降水拟合结果来看，研究区降水变化增幅约为0.76 mm·a-1。

图12 年均降水变化趋势Fig.12 Variation trend of annual mean precipitation，measured and simulated

4.1.2 研究区夏季平均气温变化特征

在分析冰川变化特征时，通常将气象数据按月平均值是否高于0℃来划分冬、夏两季［33］。由于研究区附近气象站点月平均气温观测值仅在6—9月高于0℃，故本文以上述站点在6—9月的平均气温算术均值作为年度夏季平均气温值进行统计分析。

图13 为雀儿山地区不同时间段的年度夏季平均气温以及3年滑动夏季平均气温的线性拟合结果（P＜0.01）。总体而言，1987—2016年间，雀儿山地区的年度夏季平均气温呈现显著上升的趋势（P＜0.01），上升速率约为0.30～0.31℃·a-1，是全球平均升温率0.15℃·（10a）-1的2倍［34］。由于冰川面积与3年滑动夏季平均气温数据呈现显著负相关，其相关系数达-0.76，因此认为，年度夏季平均气温的上升是引起雀儿山地区冰川30 a以来持续退缩的主要原因。

图13 夏季平均气温变化Fig.13 Variations of summer mean air temperature，measured and simulated

从不同时段来看，1984—1993年、1993—2007年以及2007—2016年三个时段的年度夏季平均气温增幅分别为0.26℃·（10a）-1、0.35℃·（10a）-1和0.46℃·（10a）-1，对应均值分别为11.64℃、11.66℃和12.43℃，三个时段的线性拟合R2＞0.56，P＜0.05，通过显著性检验。虽然三个时段的温度上升幅度逐步递增，但第二时段的夏季平均气温相对第一时段无明显上升，这可能是1993—2007年间雀儿山地区冰川面积变化幅度相对较小的主要原因。而在2007—2011年间，雀儿山地区冰川年变化率达（-6.67±6.46）%·a-1，是其他时段的数倍。为探究波动的真实性及其发生原因，本文对五个站点的夏季平均气温数据进行Mann-Kendall突变分析（图14）。结果显示，在2004—2005年间，研究区的夏季平均气温出现了快速升高的状况，该突变点通过滑动t检验，期间研究区的降水并未发生明显变化。由于研究区冰川对气候变化响应的滞后期约为2～3 a，由此推测，当地气温的突变上升，可能是2007—2011年冰川突然快速退缩的主要因素。

图14 雀儿山地区1979—2016年夏季平均气温M-K趋势变化Fig.14 Variation in Mann-Kendall trend of summer air temperature in the Que’er Mountains，1979—2016

4.2 与青藏高原部分地区冰川变化的对比

影响冰川退缩速度主要有冰川性质、冰川规模以及气候变化等因素。由表4可知，雀儿山地区是青藏高原冰川退缩最剧烈的地区之一，仅次于岗日嘎布山地区冰川。由于雀儿山地区大部分冰川小于2 km²，且大部分冰川分布于海拔5 300 m以下，加之海洋型冰川的性质决定了它对温度变化十分敏感，因而在夏季平均气温上升速率高达0.30～0.31℃·a-1的背景下，雀儿山地区冰川的退缩速率超过了青藏高原上大多数山地。

4.3 冰川提取结果的精度分析

4.3.1 不同分类方法的精度对比

本文参考相关高清遥感影像图以及调查资料，采用分层随机采样，在11个不同时期的图像中分别获取了1 018～1 994个地物类样本点作为地面真实值，对所得冰川分类结果与不同分类方法的结果进行对比（图15）。在多种分类结果中，本文结果的精度相比其他分类方法提高了2.8%～28.2%，平均kappa系数达到0.97，意味着结合地形阴影的冰川提取方法能更精确地进行冰川信息识别与提取。

图15 多种分类方法精度统计比较Fig.15 Comparison of analysis precision of different classification methods

4.3.2 与中国第二次冰川编目数据的对比

中国第二次冰川编目是以Landsat TM/ETM+图像数据为主要数据源，主要通过波段比值阈值分割方法，结合人工修订形成最终的冰川边界数据。

图16 （a）和（b）为本文提取的2003年雀儿山地区部分冰川边界与中国第二次冰川编目边界叠加在不同遥感影像的结果，可见两者的冰川边界有细小的差异。本文以10 m为间隔，从冰川提取结果边界取点，得出各点与中国第二次冰川编目对应冰川边界平均距离误差为12.48 m，不足像元边长的一半，两组数据的冰川边界差异可能是计算机分类结果的平滑、人工修正等随机误差所造成的。此外，两组冰川的面积差异仅为2.29%（表5），说明本文方法提取的冰川边界是可靠的。

4.3.3 与GAMDAM冰川编目数据的对比

GAMDAM冰川编目的主要数据源为Landsat TM/ETM+影像数据和Google Earth高分辨率图像，主要采用人工目视解译提取冰川信息［14］。

图16 （c）和（d）为本文提取的1993年冰川边界与GAMDAM冰川编目边界叠加在不同遥感影像的结果，其中图16（c）为新路海冰川末端的边界对比。由于GAMDAM冰川编目的冰川范围包含了被表碛覆盖或是分布于陡坡上的冰川，而这部分冰川难以在30 m空间分辨率遥感影像中进行识别，加之在云影响下部分小冰川也未能被识别［图16（d）］，因此本文提取的冰川范围相对较小。虽然如此，本文提取的1993年冰川面积与GAMDAM冰川编目数据相比（表6），差异仅为-6.93%，而两组冰川的边界距离误差也仅为8.17 m，考虑到随机误差，认为本文所提取的冰川边界具有较高的精度。

图16 本文冰川边界提取结果与不同冰川编目的对比Fig.16 Comparison of glacial periphery between this study and the second Chinese glacier inventory［（a）and（b）］，between this study and the GAMDAM glacier inventory［（c）and（d）］

表5 本研究与中国第二次冰川编目结果（雀儿山2003年冰川面积）的对比Table 5 Comparison of extracted glacierized area in the Que’er Mountains between this study（in 2003）and the second Chinese glacier inventory（in 2003）

5 结论

根据多期遥感数据和DEM数据，将地形阴影模拟结果引入到冰川分类方法中，提取了雀儿山地区1998—2016年间的冰川面积信息，并在此基础上结合气象数据对冰川变化及其对气候变化的响应进行了分析，得到以下结论：

（1）从冰川变化与气候响应关系上看，年度夏季平均气温的上升，是雀儿山地区冰川退缩的主要原因。1987—2016年期间，雀儿山地区冰川呈现明显的退缩趋势，冰川面积合计减少了（48.97±25.29）km2，整体年变化率为（-1.69±0.87）%·a-1，是青藏高原主要山系中冰川退缩速度最快的地区之一。

（2）雀儿山地区冰川损失面积主要来源于＜1 km²的小型冰川。从坡向上看，雀儿山地区西、南朝向的冰川退缩相比其他方向更为明显，这与冰川面积大小分布特征相反，可能与太阳辐射能量和热量等分布差异有关。此外，随海拔的上升，雀儿山地区冰川退缩率呈现“先上升后减少”的趋势，并在4 901～5 000 m处达到最大，这可能与不同海拔冰川的坡向及其阴影覆盖范围差异有关。

（3）在遥感冰川分类方法上，将传统基于光谱特征的冰川提取方法与基于GIS的阴影模拟方法进行结合，能较好地对山体阴影下的冰川进行提取，提高冰川信息识别提取的精度。这为利用遥感技术提取冰川信息的方法改进提供了一种新的思路。