亚洲高山区积雪物候时空动态及其对气候变化的响应

唐志光， 邓 刚， 胡国杰， 王 欣，， 蒋宗立， 桑国庆

（1.湖南科技大学测绘遥感信息工程湖南省重点实验室，湖南湘潭 411201； 2.湖南科技大学地理空间信息技术国家地方联合工程实验室，湖南湘潭 411201； 3.中国科学院西北生态环境资源研究院，甘肃兰州 730000）

0 引言

积雪是冰冻圈的重要组成要素，是地球表层上地域分布最广泛、季节变化最明显、对气候异常变化最敏感的覆盖物［1］。在气候变化研究中，积雪的季节变化能导致地表反照率变化，进而引起地气能量收支平衡和区域水平热力差异。在全球水循环过程中，积雪的积累和消融过程起到水的年内再分配作用，是干旱半干旱地区春季最重要的淡水资源［2］。积雪影响着地球陆地上中高纬度地带和高海拔地区的生态地理环境，关系着农牧业生产、生态系统和社会经济的发展。

亚洲高山区（又称高海拔亚洲）是地球上除南北两极之外的第三大冰冻圈，被称为地球的第三极，是全球气候变化研究的热点区，也被称为“亚洲水塔”［3］，是诸多大江大河的发源地，也是“一带一路”的核心区。在全球气候变暖背景下，亚洲高山区的冰冻圈（冰川、积雪等）正经历着快速的变化，但这种变化及其对全球变化的响应在不同时空尺度表现出明显差异［2，4-12］。亚洲高山区冰冻圈的变化将引发水资源减少、极端天气事件增多、灾害频发等生态和环境问题，深刻影响相关国家和地区的经济社会可持续发展，已受到各国的广泛关注［7，13-14］。亚洲高山区积雪的时空变化是全球气候变化的敏感指示器，同时也影响着全球气候变化。开展亚洲高山区积雪变化的时空差异性研究，可为全球气候变化研究、区域水资源的管理和经济社会可持续发展等提供科学依据。

冰冻圈遥感技术最大的优势在于获取大范围冰冻圈要素的时空分布信息，已成为“三极”大科学计划的重要技术手段［15-16］。近年来，不少学者利用遥感技术的优势，分析了不同时空尺度下的积雪变化特征，主要包括积雪范围、雪深以及积雪物候的变化。积雪物候信息是衡量积雪变化特征的重要指标，是寒区水文模型中融雪径流模拟的重要参数，能为干旱区河川径流预测提供重要数据支撑。同时积雪物候变化对陆地生态系统产生重要影响，尤其是显著影响着高寒植被的物候期［17-18］。Peng等［19］发现1980—2006 年欧亚大陆大部分气象站的积雪结束日期提前，开始日期推后，积雪持续时间降低，并且这些变化与气温升高密切相关。Ke等［20］和Ma［21］分别采用1952—2010 年和1970—2014 年的中国气象台站数据，得出了中国的积雪物候总体上呈现结束日期提前，开始日期推后的结论。汪箫悦等［22］采用MODIS 积雪产品分析了2002—2012 年青藏高原积雪物候的时空变化特征，发现18.1%的地区积雪开始期明显提前（推迟占8.5%），23.2%的地区积雪结束期明显推迟（提前占6.9%）。乔德京等［23］利用1980—2009 水文年的被动微波遥感雪深数据从20 世纪80 年代、90 年代和21 世纪初三个时间段分别分析了青藏高原积雪物候变化特征和异常分布。

本文基于2000—2020 年的逐日MODIS 积雪产品，在对其进行去云处理的基础上提取了亚洲高山区每一水文年的积雪物候参数（积雪日数、积雪开始日期、积雪结束日期和积雪持续日数）。并对研究区积雪物候变化的时空差异及其对气温、降水变化的响应关系进行了分析。

1 研究区及数据

1.1 研究区概况

亚洲高山区（High Mountain Asia）位于亚洲中部（25°～51° N、65°～105° E），是以青藏高原为核心的亚洲高海拔区域，由各大山脉（阿尔泰山、天山、希萨尔－阿莱山、昆仑山、祁连山、横断山和喜马拉雅山等）及周边高原区组成［6］（图1）。该地区地形和气候模式复杂，拥有多种不同类型的独特地貌和生态系统。作为地球上除两极之外最大的冰雪资源分布区，亚洲高山区是许多国际性大江大河的发源地，也被称为“亚洲水塔”，其冰雪变化关系着该地区及周边区域的生活和灌溉用水，是生态环境变化中的关键环节［24-26］。本研究以亚洲高山区（海拔大于1 500 m 的区域）作为研究区，总面积约为5.14×106km2，并按照全球冰川编目（Randolph Glacier Inventory，RGI 6.0［27］）的分区开展区域性差异研究（图1）。

图1 研究区概况图Fig. 1 Sketch map of the study area

1.2 数据源及预处理

1.2.1 MODIS积雪数据及去云处理

MODIS 积雪产品（MOD10A1，V006 版）来源于美国冰雪数据中心（https：//nsidc. org/data/MOD10A1/），空间分辨率为500 m，时间分辨率为1 d，文件格式为HDF，投影方式为正弦曲线投影（sinusoidal projection）。V006 版积雪产品较V005 版具有质量和精度的提升，并采用归一化积雪指数（Normalized Snow Index，NDSI）代替了V005 版的积雪覆盖率（Fractional_Snow_Cover）。利用MRT工具（MODIS Reprojection Tools）［28］对研究区2000—2020 年的逐日积雪产品进行预处理，主要包括：影像拼接、投影转换、数据裁剪、格式转换等；转换后投影为Geographic Lat/Lon，输出格式为Geo-TIFF。并对极少数缺失数据采用其最邻近日的数据进行替代。

云的存在极大地限制了MODIS 积雪产品的积雪监测效果。不少学者基于MODIS 二值（雪或非雪）积雪产品开展去云研究，主要方法包括：Terra和Aqua（上下午星）产品融合［29-31］，临近日数合成［31-33］，临近像元空间插值［29，31，34］，雪线高度去云［30-31，35］，融合被动微波雪深［35-36］或者其他光学遥感数据［37］，以及基于时空建模技术去云［38-40］等。积雪覆盖率产品相对于二值积雪产品更能有效的反映像元内积雪覆盖的渐变（累积/消退）特征，Tang 等［40-42］发展了适用于积雪覆盖率产品的三次样条函数插值去云算法。本研究采用三次样条函数插值去云算法对MODIS NDSI 产品进行去云处理，最终得到2000—2020 年亚洲高山区逐日无云的MODIS NDSI 数据集。

1.2.2 ERA5气候再分析资料

ERA5 数据是欧洲中期天气预报中心（European Centre for Medium-Range Weather Forecasts，ECMWF）的全球气候第五代大气再分析资料，涵盖了从1979 年1 月1 日到接近实时的时间段，且与其他气候再分析资料相比在青藏高原地区具有更好的适用性［43-46］。本文利用ERA5 逐月数据集“ERA5 monthly averaged data on single levels from 1979 to present”（https：//cds. climate. copernicus. eu/cdsapp#！/search？type=dataset&text=ERA5，空间分辨率为25 km）来提取2000—2020 年亚洲高山区的月降水量和月平均气温，以探究积雪物候与气候因素之间的相关性。

1.2.3 其他数据

SRTM（Shuttle Radar Topography Mission）数字高程数据（Digital Elevation Model，DEM）用于积雪物候的空间分布特征分析。目前，SRTM V4 DEM空间分辨率为90 m，以WGS84 地理坐标系为基准，可直接在线获得（http：//srtm.csi.cgiar.org/）。为了在分析过程与MODIS 积雪数据的空间分辨率保持一致，将DEM数据重采样成500 m分辨率。

2 研究方法

2.1 去云算法精度评价

本文采用“云假设”的精度验证方法［40-41］，对三次样条函数插值法得到的云覆盖下的NDSI 精度进行验证。首先选取研究区2020 年云覆盖较少的4日（第238、261、289 和315 d）的MODIS NDSI 数据作为“真值”影像，影像中晴空条件下所有像元的NDSI值即为“真实值”；然后假设这4幅“真值”影像完全被云覆盖（即影像像元值为250），并利用去云算法对“云假设”后的影像进行去云处理，得到相应的无云NDSI（即“计算值”）和云持续日数（即为云覆盖像元最邻近的前后两次无云覆盖之间所持续的日数）；最后将“计算值”与“真值”进行对比，统计不同云持续时间条件下的平均绝对误差（MAE）和均方根误差（RMSE），计算公式如下：

式中：n为所有参与比较的像元个数；ci为MODIS 产品的NDSI 值（“计算值”）；si为对应的MODIS NDSI“真值”。

2.2 积雪物候参数提取

本文提取的积雪物候参数包括积雪日数（snowcovered days，SCD）、积雪开始日期（snow onset date，SOD）、积雪结束日期（snow end date，SED）和积雪持续日数（snow duration days，SDD），其具体定义与提取方法如下。

积雪日数（SCD）：一个水文年内观测到的积雪覆盖日数之和；

积雪开始日期（SOD）：查找一个水文年内首次出现连续5 天为积雪覆盖的事件，若存在，则取这5天中第一天的日期为SOD；

积雪结束日期（SED）：查找一个水文年内最后出现连续5 天为积雪覆盖的事件，若存在，则取这5天中最后一天的日期为SED；

积雪持续日数（SDD）：从积雪开始日期到积雪结束日期的持续日数，即两者的差值。

原MODIS V005 版的二值积雪产品以NDSI 阈值0.40 作为积雪判识条件，但已有研究［47-50］表明该阈值在青藏高原地区的积雪判识中存在较大误差，并证明了NDSI 阈值0.29 更适合于青藏高原地区的积雪制图［47］。因此，本研究在积雪物候参数的提取过程中，将NDSI≥0.29 的像元视为积雪像元，水文年定义为当年9 月1 日至次年8 月31 日。为了便于分析积雪物候参数的变化特征，若SOD 和SED出现在水文年中的第一个自然年，SOD 和SED 值以儒略日计数；若SOD 和SED 出现在该水文年中的第二个自然年内，其值的计算则以儒略日+365（或366）计数，单位记为day of year（DOY）。例如SOD 等于365、370，分别代表积雪开始日期出现在第一个自然年的12 月31 日、第二个自然年的1月5日。

2.3 趋势分析法

采用Sen’s斜率估计法计算亚洲高山区近20年的积雪物候参数的变化趋势，并结合Mann-Kendall检验方法对变化趋势进行显著性检验。Sen’s 斜率估计是一种广泛应用于气象学和水文学中的非参数统计检验方法，其优点在于可以处理非正态分布的数据，用它代替线性回归，是因为它限制了离群值对斜率的影响［51-52］。Sen’s 斜率β的计算公式如下：

式中：Median为取中值函数；β为时间序列的趋势斜率。当β＞0 时，序列呈上升趋势；当β=0 时，序列趋势不明显；当β＜0时，序列呈下降趋势。

Mann-Kendall 检验是一种单调趋势的非参数方法，已广泛应用于检测时间序列的趋势［53-55］。时间趋势的统计显著性由Z值来假设，当|Z|≥1.28、1.64、2.32 时，分别通过了置信度为90%、95%、99%的显著性检验。对于时间序列变量x1，x2，…，xn，n为时间序列长度，定义检验统计量S：

其中，sgn为符号函数

式中：当n≥10 时，统计量S近似服从正态分布。将S标准化得到Z，利用统计检验值Z进行显著性检验，计算公式如下：

其中，

3 结果分析

3.1 去云算法精度分析

统计2020 年逐日MODIS NDSI 原始影像中云覆盖像元的云持续时间，得到研究区云覆盖像元按照云持续时间的频率分布图［图2（a）］，可以看出，随着云持续时间的增加，研究区云覆盖像元出现的频率逐渐降低，云持续时间在5 d 内的云覆盖像元频率可达63.7%。图2（b）为去云算法在不同云持续时间条件下的MAE 和RMSE 分布。去云算法的MAE和RMSE 随着云持续时间的增加而增加，即云覆盖像元的云持续时间越短，去云精度越高。以云覆盖像元随云持续日数的分布频率为权重，得到NDSI 去云算法的总体MAE 和RMSE 分别为0.041和0.108。可见，本文的基于三次样条函数插值的NDSI去云算法具有较高的精度。

图2 去云算法精度检验Fig. 2 Accuracy verification of cloud removal methodology

3.1 积雪范围年内变化

图3 为近20 年研究区积雪范围的平均年内波动特征（S1～S7 七个子区域是根据亚洲高山区地形地貌的相似性和复杂程度，由16个RGI分区进行邻域合并而来）。可以看出，研究区积雪范围季节变化显著，有着明显的积雪累积期和融雪期。从整个研究区来看［图3（a）］，积雪从9月开始积累，直至次年1 月中旬积雪范围达到峰值（45.5%）；随后积雪开始融化，至8 月份积雪范围达到谷值（3.0%）；此外，较大的标准差分布表明积雪范围年际波动较大。

从图3（b）可以看出，不同子区域积雪范围的年内波动差异明显。其中，高纬度的阿尔泰－萨彦岭（S1）、天山/希萨尔－阿莱山（S2），积雪范围峰值明显高于其他子区域，且在融雪期4 月至6 月的积雪范围下降较快；较高纬度的帕米尔/兴都库什/喀喇昆仑/西昆仑（S3）和高海拔的喜马拉雅山脉（S7）年内积雪范围也均处于较高水平；而东昆仑/祁连山（S4）、青藏高原内部（S5）和青藏高原东南部/横断山（S6），年内的积雪范围较低。此外不同区域的积雪范围峰值出现的时间表现出明显差异。

图3 积雪范围水文年内变化Fig.3 The annual cycle of snow cover extent during a hydrological year

3.2 积雪物候空间分布

积雪日数（SCD）大于60 d的地区被认为是稳定积雪区，也是雪水资源的主要来源地［56］。研究区的稳定积雪区主要分布在天山、喀喇昆仑山、昆仑山、祁连山、喜马拉雅山以及横断山等高海拔山脉，以及高纬度的阿尔泰-萨彦岭［图4（a）］，占研究区总面积的44.7%；其中，SCD 在120 d 以上的区域占研究区总面积的28.0%。此外，SCD 低值区（SCD＜20 d，占30.7%）主要分布在青藏高原周边的低海拔区以及高原内部区域。青藏高原内部地区主要是由于受到四周高山地形的阻挡，水汽输入量少，因而降雪很少。

研究区积雪开始日期（SOD）、积雪结束日期（SED）和积雪持续日数（SDD）也体现出较强的空间差异性［图4（b）～4（d）］。整体上，SOD主要集中在9—12 月，SED 主要集中在3—6 月，SDD 与SCD 在数值上虽有所差异，但两者呈现出极其相似的空间分布格局。在天山、喀喇昆仑山、西昆仑、喜马拉雅山脉和横断山等高海拔山脉以及高纬度的阿尔泰－萨彦岭，SOD 普遍较早，9 月初就开始出现积雪；在空间分布上，SOD 有着明显的海拔梯度，海拔越低，SOD 越晚（图5）；低海拔区的SOD 多出现在11月及以后。从SED 的空间分布来看，在低海拔区SED 出现在3 月份及以前；SED 同样呈现出主要受海拔影响的垂直地带性分布规律，海拔越高，SED 越晚（图5）；高海拔山区则推迟到6月份及以后。

图4 2000—2019水文年亚洲高山区平均积雪物候Fig.4 Average snow phenology in the High Mountain Asia during the hydrological year of 2000—2019

图5 不同子区域积雪物候随海拔变化Fig. 5 Variations of snow phenology with elevation in different subdivisions

3.3 积雪物候年际变化

对研究区2000—2019 水文年的SOD、SED 和SDD 进行变化趋势分析和显著性检验（图6），并统计了不同区域趋势显著性所占面积比例（表1）。总体而言，大部分地区的SOD 无显著变化趋势［图6（a）、6（b）］；SOD 在8.2%的区域显著提前（P＜0.1），提前速率约为3 d·a-1，主要分布在阿尔泰山、东天山以及青藏高原东南部的部分区域；SOD 在11.4%的区域表现为显著推迟（P＜0.1）的趋势，推迟速率约为1～5 d·a-1，主要分布在萨彦岭南部、西天山、喀喇昆仑、中喜马拉雅和东喜马拉雅。从SED 的变化趋势和显著性空间分布来看［图6（c）、6（d）］，15.8%的区域SED 显著提前（P＜0.1），主要集中在阿尔泰－萨彦岭、天山以及横断山西部，提前速率为1～6 d·a-1；而显著推迟（P＜0.1）的区域占比仅为8.8%，主要分布在兴都库什南部、西喜马拉雅，推迟速率在3 d·a-1左右。研究区13.5%的区域SDD 呈现显著缩短（P＜0.1）趋势，其中，西天山、中喜马拉雅和东喜马拉雅、横断山西部、萨彦岭等地区缩短了1～6 d·a-1，兴都库什北部、喀喇昆仑、东天山以及阿尔泰－萨彦岭部分区域缩短了约3 d·a-1；7.4%的区域SDD 显著延长（P＜0.1），主要集中分布于兴都库什南部、西喜马拉雅以及青藏高原东南部，大约延长3 d·a-1［图6（e）、6（f）］。从积雪物候显著变化面积占比来看（表1），整个研究区SDD以显著缩短为主（显著缩短多于显著延长）、SED 以显著提前为主（显著提前多于显著推迟）、SOD 以显著推迟为主（显著推迟多于显著提前）。

图6 2000—2019水文年亚洲高山区积雪物候变化趋势与显著性分布图Fig. 6 Trends and significance of snow phenology in the High Mountain Asia during the hydrological year of 2000—2019

表1 积雪物候显著变化面积占比Table 1 The area proportion of significantly changed snow phenology

3.4 积雪物候变化和气候因子相关性分析

积雪物候变化是一个周期性、持续性的动态过程，对气候环境变化非常敏感。为了分析积雪物候对气候变化的响应关系，本文从不同空间区域和时间段计算了SOD 和SED 与降水和气温之间的相关关系（表2）；其中，累积期为9—12 月、融雪期为3—6 月。气温和降水对研究区积雪物候的影响表现出较强的空间差异性。但总体而言，SOD 与累积期降水呈负相关，与累积期气温呈正相关；几乎整个研究区（除东喜马拉雅以外）的SOD 均与累积期气温的相关性更大，累积期气温是影响SOD 年际变化的主导因子。SED 与融雪期降水和年降水均呈正相关，与融雪期气温和年平均气温均呈负相关；阿尔泰－萨彦岭、希萨尔－阿莱山、天山、帕米尔以及兴都库什的SED 都与融雪期气温的相关性更大（呈显著负相关），融雪期气温是其影响SED 年际变化的主导因子；而喀喇昆仑的SED 主要受融雪期降水的影响，西喜马拉雅的SED主要受年降水量的影响。

表2 不同分区积雪开始日期、结束日期与降水、温度的相关性Table 2 Correlation of snow onset date and snow end date with precipitation and temperature in different subdivisions

4 结论

本文在对2000—2020 年逐日MODIS 积雪产品进行了去云处理的基础上，提取了研究区20个水文年的SCD、SOD、SED 和SDD，从像元尺度以及不同空间分区对其年际变化趋势进行了分析，并分析了积雪物候对气候变化的响应关系，主要结论包括：

（1）基于三次样条函数插值的MODIS 积雪覆盖率去云算法，能够有效地获取云遮蔽像元的NDSI 信息。云覆盖像元的云持续时间越短，该算法的去云精度越高；该去云算法的总体平均绝对误差和均方根误差分别为0.041和0.108。

（2）亚洲高山区积雪物候空间差异性较大，并呈现出主要受海拔影响的垂直地带性分布规律。稳定积雪区（SCD＞60 d）主要分布在天山、喀喇昆仑山、昆仑山、祁连山、喜马拉雅山以及横断山等高海拔山脉，以及高纬度的阿尔泰－萨彦岭。研究区SOD 主要集中在9—12 月，高海拔山脉及高纬度地区的SOD 出现较早，而低海拔区的SOD 多出现在11 月及以后；SED 主要集中在3—6 月，在低海拔区SED出现在3月份及以前，而高海拔山区则推迟到6月份及以后。

（3）20 年间，研究区SDD 主要呈缩短趋势，在13.5%的区域呈现显著缩短趋势，而仅有7.4%的区域显著延长；SED 主要呈提前趋势，在15.8%的区域为显著提前，而显著推迟的区域仅占8.8%；SOD 主要以推迟趋势为主，在11.4%的区域表现为显著推迟，在8.2%的区域为显著提前。

（4）亚洲高山区积雪物候年际变化对气候变化有着明显的响应关系。SOD 与累积期降水呈负相关，与累积期气温呈显著正相关，累积期气温是影响SOD 年际变化的主导因子；SED 与融雪期降水和年降水均呈正相关，与融雪期气温和年平均气温均呈负相关，融雪期气温是影响SED 年际变化的主导因子。