1 引言

地表反照率是指地表反射的太阳辐射与入射太阳辐射之比，反映地表对太阳辐射的反射能力，是陆面过程模式和气候模式中的重要参数之一(Dickinson, 1983; Collins et al., 2004)。地表反照率的细微变化会显著影响地—气系统的辐射能量收支，从而影响局地、区域乃至全球的气候变化（肖登攀等，2011；孟宪红等，2015）。研究和分析地表反照率的时空分布及动态变化特征，对于改进陆面过程模拟、促进全球气候变化和地—气相互作用研究具有重要意义。

近年来，地表反照率的时空变化特征一直倍受关注，如王鸽等（王鸽，2010；王鸽等，2011）利用1982~1998年NOAA/AVHRR (Advanced Very High Resolution Radiometer)数据反演的地表反照率研究了中国地区地表反照率的时空变化特征；Zhang et al. (2003)利用1985~1998年AVHRR数据反演得到的地表反照率研究了阿拉斯加北斜坡的地表反照率变化特征；王艺等（2011）利用2000~2009年MODIS (MODerate resolution Imaging Spectroradiometer)地表反照率气候网格产品MCD43C3研究了中国地区地表反照率的分布特征及年际、季节变化趋势；Atlaskina et al. (2015)利用2000~2013年MODIS地表反照率气候网格产品MCD43C3研究了北半球春季雪盖地区地表反照率的时空变化特征。然而，文献表明针对青藏高原地表反照率时空分布及其动态变化特征的研究仍然较为缺乏。

青藏高原平均海拔约为4500 m，有“世界屋脊”和“第三极”之称，不仅范围广阔、地表类型复杂多样，而且广泛分布着冰川、积雪和冻土等高反照率地表类型，其大范围热动力作用和地—气相互作用对我国和亚洲地区的灾害性天气，甚至全球气候及生态环境变化均有重大影响（叶笃正和高由禧，1979；赵勇等，2013；岑思弦等，2014）。然而，气候和陆面过程模拟研究中仍然广泛采用参数化方案获得地表反照率（张强等，2009；蔡福等，2011）。因此，研究青藏高原地表反照率的时空分布及动态变化特征，对于全球气候变暖背景下更加准确地认识青藏高原的地—气相互作用过程及其气候效应、验证和改进气候和陆面过程模拟的研究具有重要的科学意义。

美国国家航空航天局NASA将Terra和Aqua卫星的MODIS数据联合，反演制作的全球MODIS双向反射分布函数BRDF (Bi-directional Reflectance Distribution Function)和地表反照率产品MCD43已累积了长达18年的产品。研究表明：MCD43产品的反演精度较高，在青藏高原的绝对精度可以达到0.02~0.05，完全可以满足气候和陆面过程模式的精度要求(Wang et al., 2004; 陈爱军等，2012)，而且时空分布连续性较好(Schaaf et al., 2011; Wang et al., 2012, 2014; 陈爱军等，2016)，对于研究全球地表反照率的时空分布及动态变化特征具有重要的应用价值。

本文采用2000~2016年MODIS地表反照率气候网格产品MCD43C3 (Collection005)，分析青藏高原地表反照率的时空分布及动态变化特征，这将为改进气候和陆面过程模拟研究提供参考，并有助于准确地理解青藏高原地—气相互作用过程及其气候效应。

2 数据与方法 2.1 研究数据

（1）MODIS地表反照率气候网格产品MCD43C3

MODIS地表反照率产品MCD43采用半经验的线性RossThick-LiSparse-R核驱动的BRDF模型反演，时间间隔为8 d，输入数据为过去16 d期间经过大气校正的Terra和Aqua双星MODIS数据。根据输入的多角度晴空观测数据（不少于7个），采用最小二乘法获得模拟地表各向异性散射特性最佳的3个BRDF参数，即各向同性散射、RossThick体散射核和LiSparse-R几何光学核的权重系数。然后，通过积分运算获得黑空反照率(Black-Sky Albedo, BSA)和白空反照率(White-Sky Albedo, WSA)(Schaaf et al., 2002, 2011)。MCD43以500 m、1 km和0. 05°气候网格分辨率（分别编号A、B和C）提供BRDF参数、BRDF反演质量、反照率和BRDF校正后的天底反射率NBAR(Nadir BRDF Adjusted Reflectance)等4种产品（分别编号1、2、3和4），其中：MCD43A和MCD43B采用正弦投影，MCD43C采用等经纬度投影(Schaaf et al., 2002, 2011)。

本文采用0.05°气候网格分辨率、等经纬度投影的地表反照率产品MCD43C3。该产品提供MODIS 1~7通道的黑空反照率和白空反照率，以及可见光（0.3~0.7 μm）、近红外（0.7~5.0 μm）、短波（0.3~5.0 μm）等3个宽波段的BSA和WSA。此外，该产品还提供了积雪覆盖率、当地正午时分的太阳天顶角及数据质量信息。本文所用MCD43C3数据由美国地质调查局USGS (United States Geological Survey)地球资源观测与科学EROS (Earth Resource Observation and Science)数据中心免费提供，数据质量完全可以满足科学研究的精度要求（https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mcd43c3 [2017-01-03]）。

（2）MODIS植被指数气候网格产品MOD13C1

MODIS植被指数数据集MOD13C1由Terra/MODIS数据生成，每16天合成一次，采用等经纬度投影，网格分辨率为0.05°，提供归一化植被指数(Normalized Differential Vegetation Index, NDVI)、增强型植被指数(Enhanced Vegetation Index, EVI)，以及质量信息和其他数据信息。与NDVI相比，EVI的敏感性更高，可以监测从稀疏到茂密的植被（王正兴等，2003）。本文采用MOD13C1数据提供的EVI分析植被对地表反照率变化的影响，该数据同样由EROS数据中心免费提供，数据质量完全可以满足科学研究的精度要求。

2.2 研究方法

（1）真实地表反照率计算

根据MCD43C3数据集提供的当地正午时分的太阳天顶角，通过计算获得当地天空散射比因子的近似值(Stokes and Schwartz, 1994)。然后，利用MCD43C3数据集提供的短波波段（0.3~5.0 μm）的BSA和WSA，采用线性加权平均法获得短波波段的真实地表反照率(Lewis and Barnsley, 1994)。

（2）年际变化趋势计算

本文采用一元线性回归法拟合青藏高原地表反照率的年际变化速率，即倾向值b（魏凤英，1999），其计算公式为

$ b = \frac{{n\sum\limits_{i = 1}^n {\left({i{a_i}} \right) - \left({\sum\limits_{i = 1}^n i } \right)\sum\limits_{i = 1}^n {{a_i}} } }}{{n\sum\limits_{i = 1}^n {{i^2} - {{\left({\sum\limits_{i = 1}^n i } \right)}^2}} }}, $

(1)

其中，n为研究时段的总年数；α_i为第i年的平均地表反照率值。如果b＞0，表明地表反照率呈增大的趋势，反之呈减小的趋势，b的绝对值大小反映了地表反照率变化的速率。然后，利用相关系数法对每个像元的变化趋势进行显著性检验，分析变化趋势的可信度。

（3）变率计算

为了定量分析高原地表反照率的年际变率和异常变化情况，计算了年和季节平均地表反照率标准差。地表反照率标准差的绝对值表示地表反照率异常和偏离平均值的程度，其空间分布的高值区是地表反照率异常变化的敏感区（除多，2016），其计算公式为

$ \sigma = \sqrt {\frac{1}{{n{\rm{ - }}1}}\sum\limits_{i = 1}^n {{{\left({{a_i} - \bar a} \right)}^{^2}}} }, $

(2)

其中α为多年平均地表反照率。

（4）自然带划分

青藏高原地形复杂，地势总体呈西北高、东南低的特征，高山大川密布，除祁连山平均海拔约为4 km外，昆仑山、冈底斯山、巴颜喀拉山、喀喇昆仑山等山脉海拔均在5 km以上，山脉之间分布着地表起伏小、面积广阔的高原和盆地（图 1a）。MODIS IGBP (International Geosphere-Biosphere Program)地表分类数据显示：从东南至西北，青藏高原依次分布着森林、草地、荒漠等植被类型，喜马拉雅山、念青唐古拉山、横断山、昆仑山、祁连山等山脉地区分布着冰川积雪（图 1b）。

为了更好地分析青藏高原不同区域地表反照率的变化及对气候响应的区域差异，本文参考郑度（1996）的青藏高原自然地带划分法，将青藏高原按温度带、地表干湿状态分为11个自然地带（图 1、表 1），逐一分析各自然地带的地表反照率时空分布及变化情况。

此外，研究还应用MCD43C3数据集提供的积雪覆盖率数据和植被指数产品MOD13C1数据集提供的增强型植被指数数据分析积雪和植被变化对地表反照率的影响。

3 结果与分析 3.1 高原地表反照率多年平均特征

（1）空间分布特征

多年平均地表反照率能够反映区域地表反照率的整体情况。分析青藏高原2000~2016年平均地表反照率的空间分布，可以看出：青藏高原近17年平均地表反照率的空间分布存在明显的差异，东南部较低、西北部较高。喜马拉雅山、昆仑山、祁连山、巴颜喀拉山、唐古拉山以及念青唐古拉山等地是地表反照率高值中心，一般超过0.4，喀喇昆仑山脉的乔戈里峰甚至超过0.7。这与当地属于高海拔地区，常年分布着大量具有强反射能力的冰川、冻土和积雪密不可分。青藏高原西北部地表反照率较高，一般为0.2~0.3，这与当地属干旱、半干旱区，大多为高寒荒漠地带有关。青藏高原东南部则是地表反照率低值区，一般为0.1~0.2，这主要是因为当地属湿润、半湿润区，大多为草原森林地带。至于青海湖、阿牙克库湖等湖泊分布区，地表反照率一般不足0.1（图 2a），说明青藏高原地表反照率的空间分布受地形和地表覆盖影响较大。分析年均地表反照率标准差的空间分布，结果表明：年均地表反照率标准差较大的区域主要分布在祁连山、念青唐古拉山、唐古拉山、横断山喜马拉雅山以及果洛那曲等地区，有的甚至超过0.06。这些地区海拔高、地形复杂，地表反照率容易受到降雪、冰川伸缩、植被覆盖等因子的影响，因此标准差较大，是地表反照率异常变化的敏感区。羌塘高原、柴达木盆地、藏南和东喜马拉雅南翼地区年均地表反照率标准差小于0.01。这可能与当地的地表覆盖类型较为均一稳定有关（图 2b）。

进一步分析不同季节青藏高原多年平均地表反照率的空间分布特征，结果表明：春季，念青唐古拉山、巴颜喀拉山、阿玛尼卿山、喀喇昆仑山脉的乔戈里峰等地是高原地表反照率高值中心，一般超过0.6；高原西部、北部地区的灌丛、草甸、荒漠地带的地表反照率也比较高，一般为0.2~0.3，东南部和南部地区森林草原地带的地表反照率则比较低，一般为0.1~0.2（图 3a）。夏季，随着温度升高，地表积雪融化，土壤湿度升高，植被生长，只有喀喇昆仑山脉的乔戈里峰等高海拔地区的地表反照率超过0.6，柴达木、昆仑北翼等荒漠、高寒草原地带的地表反照率较高，大约0.2~0.25，其余大部分地区的地表反照率均小于0.2（图 3b）。秋季，随着降雪的发生和植被的枯萎，果洛那曲和青南等高寒灌丛草甸地带地表反照率值增至0.3~0.4，念青唐古拉山、祁连山等山脉地表反照率也增至0.4左右（图 3c）。冬季，高原大部分地区地表反照率都比较高，一般为0.2~0.4，喜马拉雅山南部、唐古拉山、念青唐古拉山以及祁连山等高海拔地区是地表反照率的高值中心，一般都超过0.4，尤其是祁连山和青海湖在冬季地表反照率达到一年中的最高值（图 3d）。

总体而言，高原地表反照率四季的空间分布变化明显，冬季的地表反照率平均值最高，其次是春季和秋季，夏季的地表反照率平均值最低。喜马拉雅山、昆仑山、唐古拉山、念青唐古拉山、巴颜喀拉山、阿玛尼卿山、祁连山等高海拔山脉地区以及果洛那曲和青南等高寒灌丛草甸地带都是高原地表反照率年内变化较为剧烈的地区。

图 4给出了青藏高原不同季节平均地表反照率标准差的空间分布，可以看出：春季，平均地表反照率标准差大于0.07的主要分布在喜马拉雅山、念青唐古拉山、唐古拉山、巴颜喀拉山和阿玛尼卿山等地区，表明这些地区春季地表反照率年际差异较大。其余大部分地区平均地表反照率标准差均小于0.05，春季地表反照率年际差异较小。夏季，高原整体平均地表反照率标准差小于0.01，说明高原夏季地表反照率年际差异小；东喜马拉雅南翼地区受西南暖湿气流和地形抬升作用的影响，多云雨天气，无法获得地表反照率反演结果，因而标准差为无效值。与夏季相比，秋季平均地表反照率的标准差明显增大，尤其在祁连山、阿玛尼卿山、巴颜喀拉山以及昆仑高寒地区，这些地区地表反照率标准差甚至超过0.1。冬季，祁连山、喜马拉雅山以及果洛那曲地区平均地表反照率标准差均大于0.09，说明高海拔山脉和中部的高寒灌丛草甸地带的地表反照率年际差异大。上述分析结果表明：高海拔地区始终是地表反照率异常变化的敏感区。

（2）时间变化特征

分析青藏高原17年平均地表反照率的年变化特征，可以发现：高原地表反照率年变化一般介于0.19~0.26，在一定程度上表现为“双峰单谷”型：峰值位于冬末春初（年积日第49天）和秋季中期（年积日第297天），平均地表反照率值分别为0.26和0.24，谷值位于夏末（年积日第241日），平均地表反照率值约0.19。和积雪覆盖率的时间变化进行对比，可以看出：青藏高原地表反照率与积雪覆盖率的时间变化正相关。春季，随着积雪的消融，地表反照率逐渐减小，并在夏末降到最低值。进入秋季，伴随着降雪的发生，地表反照率逐渐增大并在秋季中期达到峰值。冬季，随着高原积雪覆盖率的不稳定变化，地表反照率也有小幅振荡，并在冬末春初达到峰值（图 5a）。另外，将地表反照率的周期处理成16 d，与植被指数的时间变化进行对比，结果表明：青藏高原地表反照率与植被指数的时间变化呈较好的反相关。春季，随着植被的萌芽生长，植被指数逐渐增大，地表反照率逐渐减小，在夏末二者达到最低值。进入秋季，伴随着植被的枯萎，植被指数逐渐减小，地表反照率则逐渐增大（图 5b）。由此可知：青藏高原地表反照率的年变化与地表覆盖类型的季节变化密切相关。

3.2 高原地表反照率年际变化特征

（1）空间变化趋势

采用一元线性回归法逐像元拟合青藏高原年均地表反照率的年际变化速率，结果表明：高原年均地表反照率的年际变化存在明显的空间差异，大部分地区的年均地表反照率呈减小趋势，减小速率约为1.0×10^-3 a^-1，其面积约占高原总面积的66%。年均地表反照率减小较快的地区主要分布在高原南部、东北部的祁连山地区、阿玛尼卿山地区以及念青唐古拉山等地区，这些地区年均地表反照率的减小速率约为6.0×10^-3 a^-1，其中川西—藏东针叶林带的西南部地区年均地表反照率减小得最快，减小速率超过1.0×10^-2 a^-1。在青藏高原西南部、中部和东北部等地区年均地表反照率有缓慢增加趋势，增加速率约为7.0×10^-4 a^-1（图 6a）。

显著性分析结果表明：年均地表反照率年际变化通过0.1显著性水平检验的地区约占高原总面积的10%，青东祁连山中南部地区、柴达木边缘地区、昆仑北翼局部地区和川西藏东等地区年均地表反照率年际变化速率为负且通过0.1显著性水平检验，说明17年间这些地区的年均地表反照率显著减小。柴达木盆地中部、青南、羌塘和藏南的零星地区年均地表反照率年际变化速率为正且通过0.1显著性水平检验，说明17年间这些地区年均地表反照率显著增加（图 6b）。

土地利用/覆被变化(LUCC)是气候变化研究的重要内容，LUCC直接改变地表反照率，从而对气候变化产生影响。本文通过分析17年间青藏高原积雪覆盖率和植被指数的年际变化情况来研究其对地表反照率的影响，可以发现，积雪覆盖率的年际变化趋势空间分布和地表反照率的年际变化趋势空间分布具有很好的对应关系：高原西部边缘、东北部以及中部地区积雪覆盖率增加，地表反照率也增加；祁连山地区、阿玛尼卿山地区、唐古拉山以及念青唐古拉山等地区积雪覆盖率减小较快，尤其是察隅和墨脱地区，相应的这些地区地表反照率也减小较快（图 7a）。这说明青藏高原地表反照率的年际变化和积雪覆盖率的变化密切相关。植被指数的年际变化趋势空间分布和地表反照率的年际变化趋势空间分布对应关系也较好：高原西部边缘和中部部分地区植被覆盖变化呈退化趋势，地表反照率呈增加趋势，除个别地区外，大部分地区植被覆盖呈改善趋势，地表反照率呈减小趋势（图 7b）。这表明高原植被覆盖整体改善是引起地表反照率减小的一个重要因素。

（2）时间变化特征

分析17年间青藏高原整体和各自然带年均地表反照率的年际变化情况，结果表明：青藏高原绝大多数自然带年均地表反照率的年际变化呈缓慢波动且减小的趋势，波动范围约为0.01~0.02；只有藏南(ⅡC1)、阿里(ⅡD1)和青南(ⅠC1)地区年均地表反照率的年际变化呈微弱增加的趋势；果洛那曲(ⅠB1)和阿里(ⅡD1)地区的年均地表反照率年际变化振幅较大，波动范围约为0.04~0.07（图 8a）。统计分析各自然带年均地表反照率的年际变化速率，结果表明：青藏高原整体的年均地表反照率年际变化速率约为－0.4×10^-3 a^-1，藏南(ⅡC1)、阿里(ⅡD1)和青南(ⅠC1)地区年均地表反照率增加速率较小，不足0.2×10^-3 a^-1，其余自然带年均地表反照率的年际变化速率约为－0.4×10^-3~ －1.1×10^-3 a^-1，其中川西藏东(ⅡAB1)地区减小得最快，年际变化速率达－1.1×10^-3 a^-1（图 8b）。进一步统计分析各自然带年均地表反照率的标准差，结果表明：阿里(ⅡD1)、果洛那曲(ⅠB1)、昆仑(ⅠD1)和青南(ⅠC1)等地区年均地表反照率标准差高于0.02，与其他自然带相比，这些自然带地表反照率年际变化波动较大，这可能是因为这些自然带多为高寒灌丛草甸和荒漠地带，地形复杂，地表反照率容易受到降雪、植被覆盖等因子的影响。东喜马拉雅南翼(OA1)地区的地表反照率年际标准差最小，不足0.011，这与该地区多高密度的森林覆盖区，地表均一稳定密不可分（图 8c）。

分析青藏高原整体的地表反照率与积雪覆盖率的年际变化，结果表明：2000~2016年高原整体积雪覆盖率与地表反照率年际变化趋势十分吻合（图 9a），二者相关系数达到0.93，且通过99%的显著性检验。同样，对地表反照率年际变化振幅较大的果洛那曲（ⅠB1）地区的分析结果也表明：当地地表反照率与积雪覆盖率的变化趋势也十分吻合，两者的相关系数高达0.98，且通过99%的显著性检验（图 9b）。以上分析进一步说明，青藏高原积雪覆盖率对地表反照率有着显著影响，近年来青藏高原地表反照率的减小与高原冰川消融和积雪减少密切相关。

4 小结与讨论

本文应用MODIS地表反照率气候网格产品MCD43C3数据，结合青藏高原自然带数据、积雪覆盖率和植被指数数据，统计分析了2000~2016年青藏高原地表反照率的时空分布及动态变化特征，结果表明：

（1）青藏高原地表反照率空间分布存在明显的地域差异，整体呈东南部低、西北部高的分布特征，受地形和地表覆盖影响较大。

（2）青藏高原地表反照率四季的空间分布变化明显，冬季值最高，其次是春季和秋季，夏季值最低，高海拔山脉地区、高寒灌丛草甸地区是高原地表反照率年内变化较为剧烈的地区，也是地表反照率异常变化的敏感区。

（3）青藏高原地表反照率的年内变化介于0.19~0.26，一定程度上表现为“双峰单谷”型，峰值位于冬末春初和秋季中期，谷值位于夏末，地表反照率年内变化与地表覆盖类型的季节变化密切相关。

（4）青藏高原年均地表反照率年际变化整体呈缓慢波动且减小的趋势，平均变率约为－0.4×10^-3 a^-1，减小区域约占高原总面积的66%，川西藏东针叶林带的西南部是地表反照率减小最快的地区，减小速率超过1.0×10^-2 a^-1，果洛那曲和阿里地区地表反照率年际波动较大。

（5）青藏高原年均地表反照率年际变化整体呈缓慢波动减小与季节性积雪减少和冰川消融密切相关，高原植被覆盖改善也是一个重要因素。

本研究利用长时间序列、高精度的卫星遥感数据研究了青藏高原地表反照率时空分布及其动态变化特征，一方面对青藏高原地表反照率的认知水平有了提高，另一方面，青藏高原地表反照率改变后必然会影响辐射、感热和潜热通量等，造成地—气之间能量交换的差异，进而影响区域及全球气候，因此这对气候模式中地表反照率的参数化方案有了新的要求。将地表反照率的时空变化应用到模式中，改进和提高模式的准确性，这有待进一步研究。

致谢: 感谢中国科学院地理科学与资源研究所张镱锂研究员提供青藏高原边界矢量数据；感谢鄱阳湖湿地与流域研究教育部重点实验室丁明军博士提供青藏高原自然带区划矢量数据。