0 引言

青藏高原是全球中低纬度冰川的主要分布区， 发育现代冰川36 793条， 面积49 873 km2， 冰储量4 561 km3， 为我国及周边国家重要的淡水资源库[1]， 具有“亚洲水塔”之称， 其雪冰融水对周边地区14亿居民的生产生活具有重要影响[2]。短波净辐射是冰川融化最主要的能量源， 而冰川反照率直接影响冰面对太阳短波辐射的吸收， 故成为影响冰川融化的重要因素之一[3]。

目前， 在野外实地观测冰川反照率还存在一定的困难。随着高空间分辨率、 高光谱遥感技术的发展， 国内外学者通过遥感影像进行陆面反照率的反演工作逐渐开展。王建等[4]通过对比地面实测与遥感数据， 认为中分辨率成像光谱仪图像反演获取的积雪反射特性可有效反映积雪的变化过程， 而且同类积雪在可见光波段的反射光谱基本一致。Stroeve等[5-6]通过格陵兰气候监测网络(Greenland Climate Network， GC-Net)在格陵兰岛开展冰雪反照率地面观测工作， 并用MODIS(Moderate Resolution Imaging Spectrpradiometer)卫星反照率16日产品(MOD43)和逐日产品(MOD10A1)进行了验证， 结果显示Terra反演结果略优于Aqua的反演结果。He等[7]用GLASS数据产品反演获取了格陵兰地区1981-2012年冰川反照率变化情况， 结果显示这31年冰川表面反照率呈现下降趋势， 从2000年开始反照率加速降低， 特别是在海拔1 000～1 500 m 的地区冰川反照率的降低更加明显。Tedesco等[8]通过MODIS数据对格陵兰2004-2010年夏季(5-8月)冰川反照率的反演结果进一步表明格陵兰西南地区反照率明显降低。Wright等[9]在格陵兰冰川应用ASD光谱仪测量的雪冰反照率为0.78～0.90， 并将MODIS产品和BSRN(Baseline Surface Radiation Network)的观测值进行比较验证， 结果表明MODIS产品可以很好地反演雪冰反照率。

青藏高原地区雪冰反照率的MODIS卫星反演研究已经取得了相当的进展， 发现大部分冰川反照率呈下降趋势， 但也有少数冰川反照率增大。王杰等[10]对冬克玛底冰川的反演得到2000-2009年间反照率呈微弱降低趋势。Qu等[11]反演的扎当冰川2001-2010年反照率呈下降趋势， 并与2006-2012年扎当冰川物质平衡的实测结果进行了对比， 发现反照率降低而冰川融化增强， 反照率最高的年份(2008年)冰川处于积累状态。Ming等[12]对喜马拉雅山脉地区南北坡11条冰川2000-2009年的反照率进行了反演， 其结果表明， 10年间绝大部分冰川反照率均降低， 仅有PT1冰川反照率升高。毛瑞娟等[13]反演的青藏高原木孜塔格冰川反照率(2000-2009年)不论是在消融区还是积累区都表现出增大趋势。目前， 整个青藏高原面上的冰川反照率变化及其对冰川融化影响的研究还相对缺乏。本文拟通过MODIS卫星MOD10A1产品反演获取2000-2013年间青藏高原11条冰川反照率的时空变化特征， 并结合已发表的11条冰川反照率研究结果， 深入分析2000年以来青藏高原地区冰川反照率时空变化特征及其对冰川融化的影响。

1 数据与方法

1.1 MOD10A1数据产品

MOD10A1是NASA陆地产品组开发的逐日积雪数据产品， 投影方式为正弦曲线投影， 空间分辨率500 m， 采用EOS-HDF(Hierarchical Data Format， HDF)格式进行分层存储。该数据由四部分组成： 分别是积雪分类(Snow cover)、 积雪反照率(Snow albedo)、 积雪覆盖率(Fractional snow cover)和质量控制(Quality assessment)。2000年3月以来全球范围的逐日积雪数据集由美国国家冰雪数据中心(National Snow and Ice Data Center， NSIDC)网站(http://nsidc.org/data)提供开源下载。

MODIS雪冰反照率反演算法由Klein等[14]提出， 该算法首先筛选出满足雪冰覆盖且无云条件的像元， 将符合条件的像元从MOD09反射率产品中进一步提取大气校正后的窄波段光谱反射率， 然后应用离散坐标法求解辐射传输方程并生成雪冰地物的双向反射模型， 并对非森林覆盖像元进行各向异性散射纠正， 最后基于MODIS传感器的窄-宽波段转换公式生成积雪反照率产品[15-16]。详细工作流程参见图1。

1.2 MOD10A1反照率数据提取

本文研究选取11条冰川(图2)位于4幅遥感影像中(逐日)， 产品编码分别是h24v05、 h25v04、 h25v05、 h26v06， 选取的遥感影像时间为2000年3月25日-2013年12月31日。利用MODIS投影转换工具(MODIS Reprojection Tool， MRT)将投影转换为地理坐标投影， 椭球体为1984年世界大地坐标系(WGS84)， 采用最邻近法对图像进行重采样， 图像格式转换为GeoTIFF格式， 像元大小为500 m， 然后用MATLAB软件程序对数据进行提取， 得到研究区反照率数据； 其中面积小的冰川选取了一个像元， 而面积大的冰川选取的像元数相对增多， 详细像元数见表1。图2中所选11条冰川分别是： ①费德琴科， ②木吉， ③慕士塔格， ④古里雅， ⑤藏色岗日， ⑥仲巴， ⑦普若岗日， ⑧羌塘1号， ⑨敦德， ⑩来古， 作求普； 所选他人已研究冰川11条： 1木孜塔格， 2冬克玛底， 3扎当， 4枪勇， 5东绒布， 6抗物热， 7杰玛央宗， 8纳木那尼， 9 PT3， 10 PT4， 11 PT5。

2 结果

2.1 冰川反照率及变化速率的空间分布

研究区内冰川反照率的空间分布特征见图3。

图3(a)和图3(b)的对比发现青藏高原冰川年平均反照率明显大于夏季平均反照率， 同时年平均与夏季平均的冰川反照率变化范围较大， 其变化范围分别是0.42(枪勇冰川)～0.75(PT5冰川)和0.45(PT2冰川)～0.69(东绒布和古里雅冰川)。此外， 不论是夏季还是年平均冰川反照率在南部喜马拉雅山脉地区大小差异较为明显， 而北部边缘地区的冰川反照率则相对较稳定， 但是冰川反照率的空间分布并没有体现出明显的规律性。冰川反照率的变化速率在空间分布上规律性较为明显： 南部变化速率较大， 往北逐渐降低， 北部变化反照率增大； 而且夏季变化速率明显大于年变化速率[图3(c)、 3(d)]。也就是说， 青藏高原北缘反照率变化较小， 而且有一定的增加趋势； 冰川反照率相对稳定的中部地区变化速率存在小幅度降低； 南部的喜马拉雅山脉地区冰川反照率的降低速率整体较大， 是冰川融化最为剧烈的地区。

注： 海拔范围是由谷歌地球测量的大致结果。

2.2 冰川反照率及变化速率的变化趋势

2000年以来， 研究区内冰川反照率的年平均值和夏季平均值变化呈现一致性。年际尺度上， 除个别冰川外， 总体表现出波动降低的趋势； 季节尺度上， 绝大部分冰川年平均反照率大于夏季平均反照率(图4、 5)。夏季平均反照率和年平均反照率的降低速率最大值都出现在枪勇冰川， 分别是-0.019 a-1和-0.015 a-1。同时少数冰川(木吉冰川、 木孜塔格冰川)的反照率呈现波动上升趋势， 木孜塔格冰川和木吉冰川反照率的年变化速率分别是0.003 a-1和0.001 a-1， 两条冰川夏季变化速率都是0.001 a-1。从图4木吉冰川反照率变化可以看出， 2012年冰川反照率明显高于其他年份， 将2012年反照率数据去掉之后得出的年变化速率由0.001 a-1变为-0.002 a-1， 夏季变化速率由0.001 a-1变为-0.001 a-1， 可见2012年的反照率高值导致整个时间序列上变化速率呈增加趋势。木孜塔格冰川反照率升高主要是由于该地区气温降低， 降水增多造成的[16]。

在藏东南地区的来古冰川和作求普冰川， 年平均反照率与夏季平均反照率规律一致， 但差值较大(图4)。结合藏东南地区的冰川物质平衡观测结果可知， 夏季冰川严重融化导致冰川反照率降低较大。

在整个研究区域内， 冰川反照率大部分呈波动降低趋势(图4、 5)， 年平均和夏季平均反照率变化速率的波动范围分别是-0.015～0.003 a-1和-0.019～0.001 a-1。变化速率呈现增加趋势的冰川主要集中在青藏高原北缘， 青藏高原中部地区变化速率相对较稳定， 而南部的喜马拉雅山脉地区变化速率差异最大[图3(c)、 3(d)]。

3 讨论

3.1 青藏高原冰川反照率降低的影响因素

在全球变暖的大背景下， 青藏高原同样处在升温过程中[17]， 气温升高是导致冰川融化的主要原因， 同时冰川反照率降低也是冰川融化的重要影响因素。2000年以来青藏高原大部分冰川反照率降低明显， 同时大部分夏季平均冰川反照率低于年平均冰川反照率。夏季气温升高， 雪冰变质和融化作用加强， 冰川融化强烈， 且粉尘、 黑碳等污化物在表雪层富集， 使冰川反照率逐步降低[18-19]， 低反照率使冰川吸收能量增多， 又促使冰川融化进一步增强； 这说明冰川融化与冰川反照率两者相互反馈。当冰川积雪融化， 冰面裸露， 至出现覆盖冰碛物和融水漫流现象时， 冰川反照率降至夏秋季节的最低值[20]。但少数冰川出现夏季平均反照率高于年平均反照率的情况(慕士塔格、 东绒布、 纳木那尼和枪勇冰川)。因极地地区温度较低， 降水主要以雪的形式降到地面； 在温度更低的冬春季节， 降水稀少， 污化物不断在冰川表面富集， 雪反照率下降[21]。慕士塔格、 东绒布、 纳木那尼和枪勇四条冰川海拔较高， 从谷歌地球测量得到最高海拔均超过6 500 m， 温度较低。由于夏季降水较多， 所以杂质难以富集， 而冬季气温更低， 降水偏少， 杂质不断在冰川表面富集， 与极地情况类似， 所以夏季平均反照率高于年平均反照率。

沙尘沉降到雪冰表面降低反照率， 同时加速雪冰融化[18,22]。Gautam等[23]通过MODIS捕捉到发生在青藏高原西南部的一次强沙尘事件， 发现大气中的沙尘同样可造成反照率的降低， 当其沉降到雪冰表面后不仅增加了对雪冰反照率的降低作用而且还延长了起作用的时间。而自二十世纪末期以来青藏高原整体气溶胶指数(AI)呈现增加趋势[24]， 说明青藏高原上空粉尘含量增多， 这也是冰川反照率降低的一个可能因素。黑碳对太阳辐射的吸收作用比粉尘更强， 从青藏高原不同地点的五根冰芯分析结果来看， 藏东南地区黑碳浓度在持续增大， 特别是二十世纪末， 黑碳浓度在迅速增大[19]。除气温升高， 冰川表面吸光性物质的增多也是青藏高原冰川反照率降低的重要原因。

3.2 青藏高原冰川反照率空间变化影响因素

从空间分布来看， 青藏高原冰川反照率空间分布没有明显规律， 但是变化速率的空间分布规律明显： 南高北低； 从局部来看， 在南部喜马拉雅山脉地区冰川反照率差异较大， 同时反照率降低速率较大， 这主要与喜马拉雅山脉地区的复杂地形有着密切的关系， Ming等[25]的研究结果显示兴都库什-喀喇昆仑-喜马拉雅地区(Hindu Kush-Karakoram-Himalaya， HKH)北部冰川反照率要大于南部冰川反照率， 冰川反照率与海拔关系密切： 3 000 m以下平均冰川反照率为0.47， 5 500～6 000 m平均冰川反照率升高至0.58， 7 000 m以上平均冰川反照率反而降低至0.56； 而中部及北部地区冰川反照率及变化速率相对比较稳定。此外， 除少数冰川外， 年平均反照率均大于夏季平均反照率。夏季气温降低、 年平均降水量增加则是青藏高原部分冰川反照率增大的重要原因[13]。1979-2010年间青藏高原降水变化表现为： 青藏高原南部尤其是喜马拉雅山脉地区降水明显减少， 青藏高原西北部与帕米尔高原东部地区降水显著增加， 而青藏高原中部地区降水比较稳定[26]。这与本文研究的青藏高原冰川反照率变化速率基本一致， 表明降水在一定程度上对冰川反照率具有重要影响。

3.3 冰川反照率与物质平衡的关系

本文统计了部分冰川物质平衡观测结果(表2)， 其中用帕隆94号冰川代表地理位置相近的来古和帕隆4号冰川， 可见除西北部的慕士塔格冰川外， 其他冰川均处在融化状态， 其中纳木那尼冰川2000-2009年冰川厚度的减薄变化速率为(-0.63±0.32) m·a-1， 而且有加速趋势[27]。从扎当冰川物质平衡与反照率关系看， 两者变化趋势基本一致， 即反照率低的年份冰川融化越严重， 冰川反照率高的年份冰川融化较弱甚至呈积累状态， 例如2008年扎当冰川具有较高的反照率(0.72)而冰川物质平衡观测结果显示出现相应的正积累[11]。物质平衡观测结果显示慕士塔格冰川在2005-2010年呈正平衡， 这期间冰川反照率微弱降低， 但是慕士塔格冰川反照率整体较高(0.65以上)， 特别是夏季反照率高于年平均反照率(图4)， 夏季又是该地区降水最多的季节， 高反照率有利于冰川的积累。

注： 扎当冰川物质平衡数据引自文献[11]， 其余冰物质平衡数据引自文献[26]。

不论是卫星还是实测数据都体现出青藏高原大部分冰川正处在退缩状态， 其中南部的喜马拉雅山脉地区退缩较为明显， 藏东南地区最为严重； 中部地区退缩并不明显； 青藏高原的北缘及西北地区冰川则处在积累状态[26,28]， 这与本文中反照率变化速率结果相一致。除了温度和降水两个主要的因素[29]， 反照率也可以作为指示冰川物质平衡变化的一个指标。

4 结论

结合已有文献资料和MOD10A1数据， 分析青藏高原2000-2013年冰川反照率的时空变化特征， 得出以下结论： 青藏高原冰川整体表现出反照率波动降低的趋势； 绝大部分冰川夏季反照率低于年平均反照率， 而夏季反照率的变化速率高于年变化速率； 此外， 冰川反照率的变化速率从南往北减小， 南部的喜马拉雅山脉地区反照率的变化速率最大， 是冰川变化最为剧烈的地区， 青藏高原北缘冰川反照率变化较小， 个别冰川有一定的增大趋势。

致谢： 感谢美国国家冰雪数据中心(National Snow and Ice Data Center， NSIDC)提供的数据支持， 感谢中国科学院青藏高原研究所张宏波、 韩孟磊、 张潇文在数据处理方面提供的帮助， 感谢中国气象科学研究院明镜副研究员提供的喜马拉雅山脉地区冰川反照率数据。

参考文献：

[1] Immerzeel W W, van Beek L P H, Bierkens M F P. Climate change will affect the Asian water towers[J]. Science, 2010, 328(5984): 1382-1385.

[2] Liu Shiyin, Yao Xiaojun, Guo Wanqin, et al. The contemporary glaciers in China based on the Second Chinese Glacier Inventory[J]. Acta Geographica Sinica, 2015, 70(1): 3-16. [刘时银, 姚晓军, 郭万钦, 等. 基于第二次冰川编目的中国冰川现状[J]. 地理学报, 2015, 70(1): 3-16.]

[3] Hock R. Temperature index melt modelling in mountain areas[J]. Journal of Hydrology, 2003, 282(1): 104-115.

[4] Wang Jian, Chen Zidan, Li Wenjun, et al. Research on spectral reflectance properties of snow using moderate resolution imaging spectroradiometer data[J]. Journal of Glaciology and Geocryology, 2000, 22(2): 165-170. [王建, 陈子丹, 李文君, 等. 中分辨率成像光谱仪图像积雪反射特性的初步分析研究[J]. 冰川冻土, 2000, 22(2): 165-170.]

[5] Stroeve J, Box J E, Gao Feng, et al. Accuracy assessment of the MODIS 16-day albedo product for snow: Comparisons with Greenland in situ measurements[J]. Remote Sensing of Environment, 2005, 94(1): 46-60.

[6] Stroeve J C, Box J E, Haran T. Evaluation of the MODIS (MOD10A1) daily snow albedo product over the Greenland ice sheet[J]. Remote Sensing of Environment, 2006, 105(2): 155-171.

[7] He Tao, Liang Shulin, Yu Yunyue, et al. Greenland surface albedo changes in July 1981-2012 from satellite observations[J]. Environmental Research Letters, 2013, 8(4): 044043. DOI:10.1088/1748-9326/8/4/044043.

[8] Tedesco M, Fettweis X, van den Broeke M R, et al. The role of albedo and accumulation in the 2010 melting record in Greenland[J]. Environmental Research Letters, 2011, 6(1): 014005. DOI:10.1088/1748-9326/6/1/014005.

[9] Wright P, Bergin M, Dibb J, et al. Comparing MODIS daily snow albedo to spectral albedo field measurements in Central Greenland[J]. Remote Sensing of Environment, 2014, 140: 118-129.

[10] Wang Jie, He Xiaobo, Ye Baisheng, et al. Variations of albedo on the Dongkemadi glacier, Tanggula range[J]. Journal of Glaciology and Geocryology, 2012, 34(1): 21-28. [王杰, 何晓波, 叶柏生, 等. 唐古拉山冬克玛底冰川反照率变化特征研究[J]. 冰川冻土, 2012, 34(1): 21-28.]

[11] Qu Bin, Ming Jing, Kang Shichang, et al. The decreasing albedo of the Zhadang glacier on western Nyainqentanglha and the role of light-absorbing impurities[J]. Atmospheric Chemistry and Physics, 2014, 14(20): 11117-11128.

[12] Ming Jing, Du Zhencai, Xiao Cunde, et al. Darkening of the mid-Himalaya glaciers since 2000 and the potential causes[J]. Environmental Research Letters, 2012, 7(1): 014021. DOI:10.1088/1748-9326/7/1/014021.

[13] Mao Ruijuan, Wu Hongbo, He Jianqiao, et al. Spatiotemporal variation of albedo of Muztagh glacier in the Kunlun mountains and its relation to dust[J]. Journal of Glaciology and Geocryology, 2013, 35(5): 1133-1142. [毛瑞娟, 吴红波, 贺建桥, 等. 昆仑山木孜塔格冰川反照率变化特征及其与粉尘的关系[J]. 冰川冻土, 2013, 35(5): 1133-1142.]

[14] Klein A G, Stroeve J. Development and validation of a snow albedo algorithm for the MODIS instrument[J]. Annals of Glaciology, 2002, 34(1): 45-52.

[15] Tekeli A E, Sensoy A, Sorman A, et al. Accuracy assessment of MODIS daily snow albedo retrievals with in situ measurements in Karasu basin, Turkey[J]. Hydrological processes, 2006, 20(4): 705-721.

[16] Hall D K, Riggs G A. Accuracy assessment of the MODIS snow products[J]. Hydrological processes, 2007, 21(12): 1534-1547.

[17] Duan Jianping, Li Lun, Fang Yongjie. Seasonal spatial heterogeneity of warming rates on the Tibetan Plateau over the past 30 years[J]. Scientific Reports, 2015, 5: 11725. DOI:10.1038/srep11725.

[18] Warren S G, Wiscombe W J. A model for the spectral albedo of snow II: Snow containing atmospheric aerosols[J]. Journal of the Atmospheric Sciences, 1980, 37(12): 2734-2745.

[19] Xu Baiqing, Cao Junji, Hansen J, et al. Black soot and the survival of Tibetan glaciers[J]. Proceedings of the National Academy of Sciences, 2009, 106(52): 22114-22118.

[20] Perovich D K, Grenfell T C, Light B, et al. Seasonal evolution of the albedo of mutiyear Arctic sea ice[J]. Journal of Geophysical Research, 2002, 107(C10). DOI:10.1029/2000JC000438.

[21] Doherty S J, Warren S G, Grenfell T C, et al. Light-absorbing impurities in Arctic snow[J]. Atmospheric Chemistry and Physics, 2010, 10(23): 11647-11680.

[22] Painter T H, Deems J S, Belnap J, et al. Response of Colorado River runoff to dust radiative forcing in snow[J]. Proceedings of the National Academy of Sciences, 2010, 107(40): 17125-17130.

[23] Gautam R, Hsu N C, Lau W K M, et al. Satellite observations of desert dust-induced Himalayan snow darkening[J]. Geophysical Research Letters, 2013, 40(5): 988-993.

[24] Gong Xiaoqian, Wu Guangjian, Zhang Chenglong, et al. Dust change over the Tibetan Plateau in recent years using ice core records and satellite remote sensing data[J]. Journal of Glaciology and Geocryology, 2012, 34(2): 257-266. [宫晓倩, 邬光剑, 张成龙, 等. 基于冰芯记录与遥感数据的近期青藏高原粉尘变化研究[J]. 冰川冻土, 2012, 34(2): 257-266.]

[25] Ming Jing, Wang Yaqiang, Du Zhencai, et al. Widespread albedo decreasing and induced melting of Himalayan snow and ice in the early 21st century[J]. Plos One, 2015, 10(6): e0126235. DOI:10.1371/journal.pone.0126235.

[26] Yao Tandong, Thompson L, Yang Wei, et al. Different glacier status with atmospheric circulations in Tibetan Plateau and surroundings[J]. Nature Climate Change, 2012, 2(9): 663-667.

[27] Zong Jibiao, Ye Qinghua, Tian Lide. Recent Naimona′Nyi Glacier surface elevation changes on the Tibetan Plateau based on ICESat/GLAS, SRTM DEM and GPS measurements[J]. Chinese Science Bulletin, 2014, 59(21): 2108-2118. [宗继彪, 叶庆华, 田立德. 基于ICESat/GLAS, STRM DEM和GPS观测青藏高原纳木那尼冰面高程变化(2000-2010年)[J]. 科学通报, 2014, 59(21): 2108-2118.]

[28] Neckel N, Kropáek J, Bolch T, et al. Glacier mass changes on the Tibetan Plateau 2003-2009 derived from ICESat laser altimetry measurements[J]. Environmental Research Letters, 2014, 9(1): 014009. DOI:10.1088/1748-9326/9/1/014009.

[29] Shi Yafeng. Concise glacier inventory of China[M]. Shanghai: Shanghai Popular Science Press, 2008.