1992—2009 年格拉丹东冰川及冰前湖面积变化的遥感研究

周文明，李志伟，李佳，崔志勇，汪长城

(中南大学 地球科学与信息物理学院，湖南 长沙，410083)

世界范围内冰川的加速退缩是全球气候不断变暖的最明显信号[1]。目前，全球气候观测系统(global climate observing system，GCOS)将冰川作为研究气候变化的基本监测对象之一[2]。同时，冰川也是地球上重要的固体淡水资源，与人类的生产生活息息相关。处于世界“第三极”的青藏高原蕴藏了丰富的冰川资源。在全球变暖的大背景下，青藏高原大部分冰川呈明显的退缩趋势[3]。青藏高原冰川的变化监测对于研究全球变化和评估区域内径流量的变化都具有重要意义。位于长江发源地的格拉丹东雪山发育了大规模冰川，这些冰川能深刻影响长江源头生态环境和中下游社会发展，且规模较大，是研究气候变化的理想标本。但是，由于该地区环境恶劣，地势险峻，常规实地测量不能满足大范围、高空间分辨率、经济和安全的冰川变化监测要求，而遥感技术凭借其面上的、无接触式的观测特点可以满足这些要求[4]。鲁安新等[5]利用人工目视判读遥感分类方法研究了格拉丹东地区冰川在小冰期、1969 年和2000 年这3 个时期的面积变化，但没有分析引起冰川变化的原因；Ye 等[6]利用非监督遥感分类和人工修整方法得到该地区冰川在1969—2002 期间的分布变化，并分析了冰川变化与沱沱河站和安多站的气象数据之间的联系，但其采用的2002年ETM+影像有季节性积雪覆盖，虽然冰川边界经过了人工处理，但结果仍然受到一定的影响，并且没有研究该地区冰前湖的时空变化规律及与冰川变化之间的关系。上述研究只涉及了该地区冰川2002 年之前的变化，而未研究2002 年之后的变化情况。完善和验证目前的冰川变化预测模型以及探究区域的气候变化机制迫切需要最新的冰川变化参数[7]。针对目前格拉丹东冰川的研究状况，本文作者综合利用遥感和地理信息系统平台，采用分类精度更高的方法得到格拉丹东冰川近20 a 来的变化规律，并分析引起这些变化的原因，同时揭示冰前湖的演化规律，发现冰前湖的演化与冰川变化密切相关，而且由于采用的分类方法不同，本文揭示的冰川退缩速度略大于文献[6]中的结果。

1 研究区与数据源

1.1 研究区概况

格拉丹东雪山位于青藏高原腹地的唐古拉山西段山区(90.0°—91.5°N，33.0°—30.7°E)，海拔5 200～6 621 m，是长江源区的最高峰[5]。该地区冰川规模较大而且集中，冰川的分布形式多为星状分布，冰川面积大约为1 000 km2(图1 中虚线矩形框标注)，冰川类型为冷性大陆性冰川。格拉丹东雪山处于青藏高原腹地，受季风气候影响很小，具有气温低和降水少的特点。该地区气象站资料显示年平均气温约为-4 ℃，每年5～9月份在0 ℃以上；年降水量在500 mm 以内，90%集中在5～9 月份。

1.2 数据源

图1 研究区域Fig.1 Location of study area

Landsat 卫星的TM/ETM+影像覆盖整个地球的时间序列较长，并且影像数据的空间分辨率和光谱分辨率均适合冰川监测，所以，TM/ETM+影像成为监测冰川变化的主要数据源，并被广泛应用于世界冰川编目中[1-2,7-10]。通常作为冰川研究的遥感影像要满足2 方面要求：1) 无云覆盖或影像中的冰川区无云覆盖；2) 影像获取时间处于或接近冰川的消融末期[1,8-9]。

综合考虑云、阴影及季节性积雪的影响，本文选取了1992 年、1999 年、2006 年和2009 年中多幅TM/ETM+数据，见表1。影像来源于美国地质调查局(United States Geological Survey, USGS)及全球土地覆盖数据数据库(global land cover facility,GLCF)。所有影像均为正射校正影像[7]。

表1 研究区域的TM/ETM+影像Table 1 TM/ETM+images of the study area

2 研究方法

目前基于光学遥感影像的冰川边界提取方法有人工目视判读、监督与非监督分类、波段比值法、NDSI(雪盖指数法)[11-16]。对于表面无冰碛物覆盖的冰川，波段比值法被认为是精确、快速、稳健的冰川提取方法，并广泛应用于冰川编目中[2,7-10]。波段比值法的基本原理是利用冰和雪在可见光近红外波段(波长为0.4～1.2 μm)的高反射特性和在短波红外(波长为1.4～2.5 μm)的强吸收特性[10]。对于Landsat TM/ETM+影像，波段比值法主要应用f(TM3)/f(TM5)和f(TM4)/(TM5)这2 个参数(其中f(TM3)，f(TM4)和f(TM5)分别代表波段3，4 和5 的原始灰度)。在阴影区和薄岩屑覆盖的冰川区域，f(TM3)/f(TM5)比f(TM4)/f(TM5)识别冰川能力更强，所以，比值f(TM3)/f(TM5)在冰川目录的编译中应用更广泛[1-2,7-8]。但是，采用波段比值法对于冰川识别来说有2 个缺陷：冰前湖在遥感图像上的光谱反射特性与冰雪的反射特性相似，使冰前湖容易被误分为冰川；冰面上有岩屑覆盖的部分与周围裸露的岩石的光谱特性相似，使得这部分冰川不能识别[8]。

本文首先通过波段比值法得到二值图，对二值图进行矢量化，得到每个时期冰川的初始边界矢量图。有的影像无季节性积雪覆盖，但部分区域受到云的影响。本文通过当年中的多期影像相互辅助，人工修整被云覆盖的区域；通过将冰川矢量边界叠加在假彩色合成的遥感影像上后人工修整冰川的边界，消除冰前湖对冰川解译结果影响[8,10]。另外，结合研究区影像的假彩色合成图，采用人工目视判读的方法提取该地区冰川湖边界。在遥感分类精度评价中，一般采用高分辨率或相同分辨率人工解译的结果作为参考数据。由于该地区缺少成像时间相同的高分辨率影像，因此，本文没有定量评价冰川和冰川湖信息提取的精度。但是，本文采用波段比值和人工目视判读相结合的方法提取冰川信息，采用人工目视判读提取冰川湖信息，不仅兼顾了时间效率，同时也保证了冰川信息提取的精度。图像处理的基本流程如下：

1) 将4 个时期的影像配准，本次实验选取1992年TM 影像为主影像，其他影像都配准到这幅影像上，配准误差都在1 个像元以内(低于30 m)。

图2 1992 年研究区假彩色图及冰川分类结果Fig.2 TM false-color composite and glacier mask of study area in 1992

2) 选取TM/ETM+影像的部分波段合成假彩色图(R,G 和B 分布对应第5，3 和1 波段)，如图2(a)所示。根据前面提到的波段比值f(TM3)/f(TM5)得到影像二值图，二值图中冰川像元值为1，非冰川区像元值为0。

3) 对冰川二值图进行窗口为3×3 的中值滤波处理，这样可以减少阴影中的噪声，同时可以消除一些季节性雪块[2,9,13]。滤波完成后将二值图(图2(b))进行矢量化处理，得到冰川的矢量边界。

4) 将冰川的矢量边界叠加在同时期的假彩色合成遥感图像上，通过人工判读检查冰川边界。修整阴影区的冰川边界，同时删除被误分为冰川的水体(尤其是冰前湖)(见图2(a)中的c 和f)；有云覆盖的区域可以通过当年不同时期的遥感影像综合判别来修改其边界；在当年不同时期的遥感影像中选择最小的冰川面积作为当年的统计冰川面积，可以最大限度地消除季节性雪的影响。

5) 统计各个时相冰川的面积参数，并计算该地区冰川1992—2009 年的分布变化数据。

3 结果与分析

3.1 冰川变化及分析

依据上述步骤获得格拉丹东雪山冰川的变化情况，见表2。从表2 可以看出：1992—2009 年，格拉丹东冰川的总面积减少了66.68 km2(减小率为7.37%)。其中，在1992—1999 年，该地区冰川面积共减少19.90 km2(减小率为2.20%)，平均每年减少2.84 km2；在1999—2006 年，冰川总面积退缩了39.95 km2(减小率为4.52%)，平均每年减小5.71 km2，与1992—1999 年相比，每年减小的面积增加2.87 km2，冰川的退缩速度加快；在2006—2009 年，总面积减少6.83 km2(减小率为0.81%)，每年减小量为2.28 km2，冰川的退缩速度有所减缓。通过上述数据可知，自1992年以来，格拉丹东冰川总体处于连续退缩之中。

表2 1992—2009 年格拉丹东地区冰川的面积变化Table 2 Glacier area change in Geladandong from 1992 to 2009

文献[5]和[6]获取了格拉丹东雪山1969—2002 年的冰川面积变化信息。本文增加该地区2006 年和2009年的冰川面积参数，同时部分覆盖了文献[5]和[6]中的监测时段，因此，可以更好地揭示该地区冰川的时空变化规律。本文综合波段比值法与人工解译获取的1999 年面积与文献[5]中手动解译的2000 年面积一致，表明本文解译的结果精度很高。由于各自采用的分类方法对处于特殊位置的积雪与冰川的辨别能力不同，本文得到的1992—1999 年时间段内的冰川年面积变化率与文献[6]得到的1992—2002 年时间段内的面积变化率略有差别，但并不影响对冰川变化趋势的判定[7]。

根据该地区冰川空间分布特点，本文将研究区划分为A 区和B 区(如图2(a)所示)。通过计算A 和B 这2 个区域面积参数(见表3)发现：A 和B 这2 个区域冰川均出现退缩，但B 区的冰川比A 区退缩更快。其原因可能是B 区的冰川面积基数更大；A 区3 个时段的冰川面积年平均变化率比较稳定，约为0.45 km2/a，B区则波动较大；1999—2006 年期间退缩速度最快，2006—2009 年期间退缩速度有所减缓。

从整体来看，冰川退缩现象广泛存在，而不是局限于某几条冰川，冰川退缩主要集中于冰舌区(消融区)。雪山东侧的冰川比西侧的退缩更严重，局部地形以及冰川中心线的走向可能是导致变化不一致的原因。其中坡度会影响冰川流速，在冰川系统中主要靠冰川运动将积累区的冰雪运送到消融区，每条冰川冰下地形的坡度与冰川流速有直接联系，进而影响冰川的物质平衡，使得冰川出现不同的退缩状况；冰川中心线的走向直接影响冰川接收太阳辐射量，而冰川接收太阳辐射量与冰川表面温度直接相关，进而影响冰川消融速度。另外，冰川大多位于山谷，当冰川两边的山峰高度较低时，冰川接收更多的太阳光能量，消融会更加剧烈。冰川的加速消融也导致某些复合冰川的末端已经开始分裂。图3 所示为图2(a)中e 冰川的末端在4 个时期的变化情况。对比图3(a)～(d)中的黑色圆圈部分，很清晰地看出在1992—2009 年间冰川消融加快，末端严重退缩，面积不断减少。图4 所示为图2(a)中a，b，c 和d 这4 条冰川的末端在不同时间的位置线。对比这4 条冰川末端位置可以发现：a，b和c 这3 条冰川处于逐年退缩状态；1992—1999 年，d 冰川末端有所前进，而在1999—2009 年间又退缩，但2009 年的冰川面积仍然比1992 年的大，冰川的加速运动或者某些年份的降水增多可能是导致这种现象的原因。

表3 1992—2009 年格拉丹东雪山A 区和B 区冰川的面积变化Table 3 Glacier area change of part A and part B in Geladandong from 1992 to 2009

3.2 冰前湖变化及分析

当冰川融化加快时，大量的冰融水以及脱落的冰块就会在冰川的前端、冰面或冰下汇集形成冰川湖，如图4(c)所示。该冰前湖在1992 年、1999 年、2006年和2009 年这4 个时期内的面积分别为0.19，0.23，0.28 和0.29 km2，呈连续扩大的趋势。图5 所示为研究区4 个时期的冰川和冰川湖面积参数。由于研究区能满足冰川监测要求的遥感数据数量有限，本文只获取研究区4 个时期的面积。由于实验数据过少，当利用函数拟合冰川和冰川湖的面积变化时，会导致拟合函数存在很大的误差，无法统计拟合误差。但是，通过分析图5 发现：冰川面积不断减少，而冰川湖面积在增加，冰川和冰川湖的面积随时间变化呈负相关，并且冰川面积和冰川湖面积呈反比例函数关系，因此，认为冰前湖面积与冰川消融面积呈非线性正相关。另外，冰川消融和冰川前端的冰块脱落速度加快是导致冰前湖水量增加的主要因素。夏季温度升高，冰川消融加强，融水注入冰川湖后使湖水上涨，伸入湖中的冰川前端底部压力和浮力随之增大，冰川冰开始脱落形成浮冰，从而使冰川物质转移至冰川湖中。为了补充冰舌损失的物质，冰体脱落和冰舌消融反过来又会加速冰川运动，使更多的上游物质抵达至此，导致冰川湖面积持续扩大。冰川湖水量激增至冰坝所能承受的最大值后就引发冰湖溃决，冰川洪水会严重威胁流域内人们的生命财产安全。另外，湖水增加会影响当地的生态平衡和农牧业，所以，冰川湖变化监测十分重要。

图3 冰川末端在4 个不同时期的退缩变化Fig.3 Examples of recession of glacier tongue in four different periods

图4 4 个时期冰川边界变化(底图为2009 年TM 影像)Fig.4 Terminus boundary of a glacier(TM images in 2009 is used as base map)

图5 研究区冰川和冰川湖在1992 年、1999 年、2006 年和2009 年的面积Fig.5 Areas of glaciers and glacial lake in study area in 1992,1999,2006 and 2009

3.3 冰川变化与气候的关系

冰川的变化与冰川所在地区的气温和降水有密切联系[3-4,17]。本文依据距离格拉丹东雪山最近的沱沱河(34°13′N,92°25′E)气象站资料[6]，整理1992—2009 年中夏季平均气温和夏季降水量，如图6 所示。从图6可见：在1992—1999 年、1999—2006 年和2006—2009年这3 个时间段内，该地区夏季月平均气温分别为5.3，5.6 和5.8 ℃，夏季降水量分别为230，300 和332 mm。夏季月平均气温一直处于升高状态。在冰川区域，夏季气温的升高会导致3 个方面的变化[17]：1) 增加外部环境融化冰和雪的能量；2) 增加液态水的比例，从而减少了冰川积累区湿雪或干雪的体积；3) 降低了冰川表面的反照率，增加冰川吸热量，加强冰、湿雪和干雪的消融。当冰川周围环境温度升高时，会导致冰川内部的温度缓慢升高，从而引起冰川内部消融。这些变化最终导致冰川表面高程降低，冰川末端的后退和雪线位置升高，使冰川出现负的物质平衡。格拉丹东地区冰川属于“夏季积累型冰川”，降水主要集中在夏季(5～9 月份)。数据显示，这3 个时间段内夏季降水量呈现增加的趋势。在冰川系统中，气温主要影响冰川的消融，降水主要影响冰川的积累。若冰川消融量与积累量相等，则冰川的物质平衡不会发生变化，冰川的总面积也不会发生变化。对于格拉丹东冰川，1992—2009 年期间冰川总面积一直减少，冰川一直处于负的物质平衡状态，说明该冰川的消融量大于积累量。虽然该区域夏季的降水量增加，使得冰川的物质积累增加，但冰川处于负的物质平衡状态。夏季升温导致冰川消融量大于积累量，因此，认为夏季降水量增加量并不能弥补夏季气温升高造成的消融增加量[17-19]。通过上述冰川的面积变化和气象数据，认为夏季气温的升高可能是导致格拉丹东冰川消融加快、面积减少的主要因素。

图6 1992—2009年沱沱河气象站的夏季(5～9月份)平均气温和降水量变化Fig.6 Annual mean summer(May to September)temperate and precipitation at Tuotuohe Station in 1992—2009

4 结论

1) 基于1992 年、1999 年、2006 年和2009 年的LandsatTM/ETM+影像，采用波段比值法和人工解译结合的分类策略，提取格拉丹东地区冰川与冰前湖面积参数并矢量化，然后对当年中多期遥感影像对应的矢量边界进行综合评判消除了阴影、云和季节性积雪的影响，以提高分类的精度。1992—2009 年，冰川面积共减少66.68 km2(减小率为7.37%)，其中在1992—1999 年、1999—2006 年和2006—2009 年分别减少19.90 km2(减小率为2.84 km2/a)，39.95 km2(减小率为5.71 km2/a)和6.83 km2(减小率为2.28 km2/a)，处于连续退缩状态，而冰前湖面积共增加了0.1 km2(增大率为52.6%)，在相应3 个时段内分别增加了0.04 km2(增大率为0.006 km2/a)，0.05 km2(增大率为0.007 km2/a)和0.01 km2(增大率为0.003 km2/a)，处于连续扩张状态。冰川与冰川湖两者紧密相连，其面积呈非线性负相关。

2) 得到的1992—1999 年时间段内的冰川年面积变化率与文献[6]中1992—2002 年时间段内的冰川年面积变化率略有差别，但并不影响对冰川变化趋势的判定。冰川退缩主要集中于冰舌区(消融区)，研究区域内几乎所有冰舌都在退缩，但雪山东侧的冰舌比西侧退缩更严重。沱沱河气象站资料显示上述3 个时间段内夏季气温和降水量都呈现上升的趋势。对于夏季积累型的格拉丹东冰川来说，夏季降水量增加量并不能弥补夏季持续气温上升造成的消融增加量，因此，认为夏季持续升温可能是导致格拉丹东冰川连续退缩的主要因素。

3) 对于大范围的冰川变化监测，在冰川末端的冰碛区，还可考虑高精度的DEM 和遥感热红外波段信息综合解译冰川冰碛区的边界变化。目前格拉丹东冰川的实地监测资料缺乏，要更深入地了解该地区冰川的物质平衡变化，还需要利用其他手段监测冰川体积变化和冰川的运动变化。在以后的研究中还将使用其他高精度的光学遥感影像、SAR 影像等综合研究格拉丹东以及其他区域冰川变化情况，以便更全面地掌握冰川与气候变化之间的关系。

致谢：

感谢美国地质调查局(USGS)及全球土地覆盖数据数据库(GLCF)提供LandSat 影像；感谢中国气象局提供的气象数据；感谢中科院寒旱所张华伟博士对本文给予的大量帮助。

[1] Paul F, Kääb A. Perspectives on the production of a glacier inventory from multispectral satellite data in Arctic Canada：Cumberland Peninsula, Baffin Island[J]. Annals of Glaciology,2005,42：59-66.

[2] Andreassen L M,Paul F,Kääb A, et al.The new landsat-derived glacier inventory for Jotunheimen, Norway, and deduced glacier changes since the 1930s[J]. The Cryosphere, 2008, 2(2)：131-145.

[3] 蒲健辰, 姚檀栋, 王宁练, 等.近百年来青藏高原冰川的进退变化[J]. 冰川冻土,2004,26(5)：517-522.PU Jianchen,YAO Tandong,WANG Ninglian,et al.Fluctuations of the glaciers on the Qinghai—Tibetan Plateau during the past century[J]. Journal of Glaciology and Geocryology, 2004, 26(5)：517-522.

[4] 叶庆华, 陈锋, 姚檀栋, 等. 近30 年来喜马拉雅山脉西段纳木那尼峰地区的冰川变化的遥感监测研究[J]. 遥感学报,2007,11(4)：511-520.YE Qinghua, CHEN Feng,YAO Tandong, et al. Tupu of glacier variations in the Mt. Naimona’nyi Region, Western Himalayas,in the last three decades[J]. Journal of Remote Sensing, 2007,11(4)：511-520.

[5] 鲁安新, 姚檀栋, 刘时银, 等. 青藏高原格拉丹东地区冰川变化的遥感监测[J]. 冰川冻土,2002,24(5)：559-562.LU Anxin, YAO Tandong, LIU Shiyin, et al. Glacier change in the Geladandong area of the Tibetan Plateau monitored by remote sensing[J]. Journal of Glaciology Geocryology, 2002,24(5)：559-562.

[6] Ye Q H, Kang S C, Chen F, et al. Monitoring glacier variations on Mt. Geladandong, Central Tibetan Plateau,from 1969 to 2002 using remote sensing and GIS technologies[J]. Journal of Glaciology,2006,52(179)：537-545.

[7] Bolch T,Yao T,Kang S,et al.A glacier inventory for the western Nyainqentanglha Range and the Nam CoBasin,Tibet,and glacier changes 1976—2009[J].The Cryosphere,2010,4(3)：419-433.

[8] Bolch T, Menounos B, Wheate R. Landsat-based inventory of glaciers in western Canada, 1985—2005[J]. Remote Sensing of Environment,2010,114(1)：127-137.

[9] Paul F, Kääb A, Maisch M, et al. The new remote-sensing derived glacier inventory： I Methods[J]. Annals of Glaciology,2002,34：355-361.

[10] Racoviteanu A E, Paul F, Raup B, et al.Challenges and recommendations in mapping of glacier parameters from space：results of the 2008 Global Land Ice Measurements from Space(GLIMS) workshop, Boulder, Colorado, USA[J]. Annals of Glaciology,2009,50：53-69.

[11] Gjermundsen E, Mathieu R, Kaab A, et al. Assessment of multispectral glacier mapping methods and derivation of glacier area changes, 1978—2002, in the central Southern Alps, New Zealand, from ASTER satellite data, field survey and existing inventory data[J]. Journal of Glaciology, 2011, 57(204)：667-683.

[12] Paul F.Evaluation of different methods of glacier mapping using landsat TM[C]// Proceedings of EARSeL SIG-Work-shop Land Ice and Snow.Dresden,2000：16-17.

[13] Raup B, Kääb A, Kargel J S, et al. Remote sensing and GIS technology in the global land ice measurements from space(GLIMS)Project[J].Computes&Geosciences,2007,33(1)：104-125.

[14] 宋波, 何元庆, 庞洪喜, 等. 基于遥感和GIS 的我国季风海洋型冰川区冰碛物覆盖型冰川边界的自动识别[J]. 冰川冻土,2007,29(3)：456-462.SONG Bo, HE Yuanqing, PANG Hongxi, et al. Identifying automatically the debris-covered glaciers in China’s monsoonal temperate-glacier regions based on remote sensing and GIS[J].Journal of Glaciology and Geocryology,2007,29(3)：456-462.

[15] 王建. 卫星遥感积雪制图方法对比与分析[J]. 遥感技术与应用,1999,14(4)：29-36.WANG Jian. Comparison and analysis on methods of snow cover mapping by using satellite remote sensing data[J]. Remote Sensing Technology and Application,1999,14(4)：29-36.

[16] 张世强, 卢健, 刘时银. 利用TM 高光谱图像提取青藏高原喀喇昆仑山区现代冰川边界[J]. 武汉大学学报(信息科学版),2001,26(5)：435-440.ZHANG Shiqiang, LU Jian, LIU Shiyin. Deriving glacier border information on Qinghai Tibet by TM high spectrum image[J].Geomatics and Information Science of Wuhan University, 2001,26(5)：435-440.

[17] Narama C, Kääb A, Duishonakunov M, et al. Spatial variability of recent glacier area changes in the Tien Shan Mountains,Central Asia,using Corona(～1970),Landsat(～2000),and ALOS(～2007) satellite data[J]. Global and Planetary Change, 2010,71(1/2)：42-54.

[18] Khromova T E, Osipova G B, Tsvetkov D G, et al. Changes in glacier extent in the eastern Pamir, Central Asia, determined from historical data and ASTER imagery[J]. Remote Sensing of Environment,2006,102(1/2)：24-32.

[19] Paul F, Huggel C, Kääb A. Combining satellite multispectral image data and a digital elevation model for mapping debris-covered glaciers[J]. Remote Sensing of Environment,2004,89(4)：510-518.