1 引言
青藏高原是世界上海拔最高的高原， 被誉为“世界屋脊”、 “地球第三极”， 其动热力效应对高原周边甚至全球天气和气候变化产生深远影响（叶笃正等， 1979； Qiu， 2008； Wu et al， 2011； 乔钰等， 2014； Yao et al， 2019）。而高原积雪作为重要的下垫面因子， 具有“高反照率”“高发射率”以及“低热传导率”的特点， 其通过热力和水文效应， 改变地表辐射以及土壤水分属性（柯长青和李培基， 1998； Zhao et al， 2007； 朱玉祥等， 2007； 王顺久， 2017； 罗江鑫等， 2020； 王梓月等， 2022）， 在天气和气候变化中起着至关重要的作用。特别是青藏高原积雪表现出显著的年际变化特征， 对中国东部夏季降水预测具有一定指示意义（Yeh et al， 1983； Brown and Mote， 2009）。鉴于高原积雪年际变化对东亚地区短期气候预测的重要性， 过去对积雪异常对亚洲夏季风和降水变化影响和机理已开展了大量研究（陈乾金等， 2000； 张顺利和陶诗言， 2001； Wu and Qian， 2003； 穆松宁和周广庆， 2010； Xu et al， 2012； Liu et al， 2014； 张人禾等， 2016； Xiao and Duan， 2016； Wang et al， 2017； 保云涛等， 2018； 段安民等， 2018； 车涛等， 2019； 曹言超和王晓春， 2022）。在全球尺度上， 高原积雪可以通过改变背景风场以及湿度场来影响北半球夏季季节内振荡（BSISO）， 从而影响印度、 孟加拉湾、 西北太平洋的对流活动（曹言超和王晓春， 2022）。高原的雪灾频数还与赤道中东太平洋和印度洋海温异常模态相联系（刘彩红等， 2021）。而对于我国， 高原冬春积雪与来年中国南部至长江中下游地区降水具有较好的正相关（Wu and Qian， 2003； 高涛和谢立安， 2014； 冯文霞和蔡晓军， 2018； 宋敏红等， 2020）。高原积雪偏多年， 积雪高反照率使得其上方对流层温度降低， 高原地表加热减弱， 高原垂直运动受抑制， 减小海陆热力梯度对比（Liu et al， 2014）， 同时， 伴随着大气环流的扰动还可以激发出异常波列（朱玉祥等， 2009）， 导致副热带高压、 西风急流和南支气流异常以及夏季风强度减弱位置偏南（Xiao and Duan， 2016）， 长江流域夏季降水偏多（张顺利和陶诗言， 2001）， 华南与华北夏季降水偏少（Wu and Qian， 2003； Wang et al， 2017）。而在最近， 罗红羽等（2023）回顾了高原冬春季积雪作为前期信号对我国旱区降水的预测研究， 并指出其内在的影响机制仍需更进一步的探讨。
由于特殊的地势， 高原上气象观测站点密度低且分布不均匀， 具体表现为高原西部站点分布较少， 中东部站点分布较多， 而中西部地区由于自然环境限制几乎没有台站， 即使东部有限的站点资料也存在缺测值较多的情况， 且开始观测的年代不统一。目前观测资料不足以及资料的代表性和精度仍是影响和制约积雪变化及其天气气候影响机理研究的一个重要因素（曹梅盛， 1995； You et al， 2008； 段安民等， 2018）。随着同化技术和卫星遥感等观测技术的发展， 再分析数据和卫星反演积雪数据产品精度和空间分辨率得到了提高， 在全球雪冰监测领域应用越来越广泛， 也很大程度上弥补了气象站点（尤其高原中西部）资料稀缺的不足（车涛等， 2004； Wu et al， 2011； 唐志光等， 2017； 周婵， 2017； 李文杰等， 2018； 黄晓东等， 2019； 闵文彬等， 2021）。但是由于不同数据之间的获取以及处理方式存在差异， 这些再分析数据和卫星反演的高原积雪数据无论在精度还是分辨率上都存在一定差距。实际研究中， 学者们往往仅选取一种或两种积雪资料分析不同季节青藏高原积雪异常的气候效应， 致使分析结果可能依赖于不同资料的选择， 存在较大的不确定性。因此， 不同类型的积雪资料（如再分析数据和卫星反演的积雪数据）是否能够较好描述高原积雪时空变化特征亦是要回答的一个重要科学问题。
为了探究不同来源积雪数据间的差异， 国内学者进行了很多相关研究。例如， 高文良等（2002）通过对NCEP/NCAR和SMM/I两种积雪资料的对比分析， 表明两种积雪资料反映出青藏高原的积雪变化特征比较一致， NCEP/NCAR积雪资料具有一定的可信度。韦志刚等（2002）利用青藏高原地面站雪深、 NOAA/NESDIS 卫星观测的雪盖面积资料以及美国宇航局微波遥感积雪水当量等资料进行了对比分析， 指出青藏高原地面站雪深资料能较好地反映出青藏高原地区积雪量的年际变化， 且其反映的冬春积雪与我国夏季降水有很好的相关。陈栋等（2006） 利用美国冰雪中心（NSIDC）提供的北半球冰雪盖周资料和欧洲数值预报中心（ECMWF）的ERA40积雪深度再分析资料研究分析青藏高原积雪的时间和空间的分布特征。王芝兰等（2013）对比1978-2010 年中国雪深长时间序列数据集与台站观测资料后指出， 对于积雪的显著性特征及气候尺度变化特征来说， 遥感资料与台站资料具有很好的一致性， 但前者对积雪的季节内波动等细节表现较差。陈爱军等（2018）利用台站雪深观测资料分别对青藏高原IMS雪冰产品精度进行了验证， 结果表明两者的总体准确率在85%以上， 指出利用IMS雪冰产品监测高原积雪范围是可靠的。除多等（2021）利用站点雪深资料对NOAA/IMS雪冰产品在青藏高原的精度和适用性进行了验证和评估， 研究发现IMS分辨率从4 km 到1 km 总体精度平均提高了2.9%， 积雪分类精度平均提高了0.9%。上述研究为青藏高原再分析和卫星反演的积雪资料使用提供了很好的参考。但需要指出的是， 已有研究多关注于积雪资料精度的评估检验， 而对更具有短期气候预测的高原积雪年际变化特征的检验尚不多。即使稍有研究关注， 但也受制于研究资料时段较短、 资料品种相对单一等条件的限制， 还缺少对多种资料综合比较结果的认知。
因此， 本文基于青藏高原台站观测和不同再分析（ERA5 和NOAA-V3）、 MODIS卫星反演以及IMS雪盖产品等多源积雪资料， 采用偏差、 均方根误差以及相关分析等多元统计方法， 通过比较不同积雪资料的时空分布和变化区域及其在积雪异常年份的空间分布特征， 重点检验了多源高原积雪数据对高原积雪年际变化尺度特征的刻画能力， 以期提升多源高原积雪资料不确定性的认知， 并为相关后继研究可提供有意义的参考。
2 数据来源和方法介绍
2.1 积雪数据
本文使用的高原积雪站点观测来自国家气象中心整理的积雪专题数据集， 包含全国2400多个站点1980 -2017年的积雪天气现象、 积雪深度、 积雪雪压、 积雪增量等数据。图1为青藏高原积雪观测站点分布， 可见50个测站多位于高原中东部， 而高原西部及南部地区站点分布较少， 因此本文在对比不同积雪数据的年际变化时， 先确定研究范围为高原东部（90°E -104°E， 28°N -38°N）。
为验证站点观测资料的准确性， 选取两种再分析资料提供的积雪数据， 一种为美国国家海洋和大气管理局（National Oceanic and Atmospheric Administration， 简写为NOAA）提供的Twentieth Century Reanalysis （V3版本）再分析数据集， 该数据下载网址https： //psl.noaa.gov/data/gridded/data.20thC_ReanV3.monolevel.html， 另外一种数据为欧洲中期天气预报中心（European Centre for Medium-Range Weather Forecasts， ECMWF）提供的ERA5再分析数据集， 数据下载地址https： //cds.climate.copernicus.eu/cdsapp#！/dataset/reanalysis-era5-land-monthly-means？tab=form。20CR（Twentieth Century Reanalysis）采用集合滤波资料同化的方法， 每三小时评估一次全球大气的状态， V3为最新版本。ERA5地表逐月数据是通过对ERA5气候资料的地表数据而产生的， 再分析资料结合了模式资料和观测资料。
三类卫星遥感雪盖数据分别选取罗格斯大学全球积雪实验室提供的北半球逐周雪盖变量（数据下载网址https： //climate.rutgers.edu/snowcover/docs.php？target=datareq）、 美国国家冰雪数据中心提供的IMS（the Interactive Multisensor Snow and Ice Mapping System）逐日雪盖产品（数据下载网址： https： //nsidc.org/）和MODIS（moderate-resolution imaging spectroradiometer）反演的积雪数据（数据下载网址https： //modis.gsfc.nasa.gov/）： （1）IMS雪冰产品是多种光学数据与微波数据集合的融合产品， 提供逐日无云的北半球积雪覆盖产品， 随着更多卫星遥感数据的融合， 产品精度和空间分辨率得到了提高， 在全球雪冰监测领域应用越来越广泛。这里采用的IMS冰雪产品空间分辨率为24 km、 该数据是由而成积雪产品， 时间为1999 -2012年逐日数据。（2）MODIS传感器全称中分辨率成像光谱仪， 搭载于EOS卫星上， 其主要有36个光谱波段， 覆盖从可见光（0.4 µm）到热红外（14.4 µm）的全光谱范围。MODIS多波段覆盖的特点可以同时反应大气的多种特征信息， 是大范围积雪和雪灾监测最为有效的手段， 也是青藏高原积雪业务化监测的主要卫星遥感数据。MODIS利用积雪在可见光波段的反射率较高、 而近红外波段反射率较低。本文使用的MODIS雪盖产品时间为2001 -2012年， 水平分辨率为0.05°×0.05°； （3）还使用了罗格斯大学全球积雪实验室提供的北半球雪盖数据为逐周数据， 时间为1966年至今， 水平分辨率为25 km， 原始的积雪数据来自AVHRR、 GOES等可见光波段的卫星数据。
2.2 分析方法
2.2.1 偏差及均方根误差的计算
本文分别计算了NOAA-V3和ERA5再分析雪深数据与站点雪深之间的偏差（BIAS）与均方根误差（RMSE， Root Mean Square Error）， 具体计算公式为：
式中： 与 分别为再分析雪深数据和站点雪深数据； BIAS直接反映了再分析雪深与站点雪深之间的数值差异； RMSE评估了再分析雪深数据和站点雪深数据之间的偏离程度。
为了将格点数据与站点数据进行匹配， 将每个站点周围的四个格点的平均作为与站点相对应的值。由此分别计算每个站点与其周围四个格点的偏差与均方根误差， 再对求得的值进行平均得到与每个站点对应的偏差值和均方根误差。分别求出每个站点对应的偏差与均方根误差之后再进行多年平均得到最终的结果。
2.2.2 卫星雪盖数据处理
本文选取了三种卫星遥感积雪数据与站点积雪数据作对比， 除了罗格斯大学提供的北半球雪盖为逐周数据外， MODIS雪盖和IMS雪盖都为逐日数据。通过计算每个格点上有积雪的日数占该月总日数的百分比， 可以得到该月的积雪覆盖率。同样的， 将罗格斯大学提供的北半球雪盖周数据计算成逐日数据后， 再计算每个格点上有积雪存在日数的占比来表示该月的积雪覆盖率。由于罗格斯大学提供的北半球雪盖与IMS雪盖都为等面积可伸缩地球网格数据， 需要将其插值为规则的经纬度网格。根据格点数的多少分别将其插值为分辨率为0.5°×0.5°和0.25°×0.25°的规则经纬度网格数据后再进行统计。
2.2.3 相关分析法
相关分析法主要通过计算两个变量之间的相关系数来判断两者之间的相关关系。本文主要通过计算皮尔逊（Pearson）相关系数来表征两个场变量之间的相关关系。皮尔逊（Pearson）相关系数是描述两个随机变量线性相关的统计量， 一般用r表示， 计算公式如下：
式中： 相关系数r的取值在-1.0~1.0之间。当r的值大于0时， 表示两变量之间呈现正相关关系， 且越接近1.0， 正相关越显著； 当r的值小于0时， 表示两变量之间呈现负相关关系， 且越接近-1.0， 负相关越显著。
3 结果与讨论
3.1 再分析资料中雪深数据的误差估计
在对积雪数据的时间变化特征进行比较之前， 先对NOAA-V3和ERA5两种再分析雪深数据与站点观测雪深之间的偏差进行分析。分别计算NOAA-V3和ERA5雪深数据和站点观测雪深之间的偏差（BIAS）和均方根误差（RMSE）来判断再分析积雪数据与观测数据之间的差异及偏离程度。由于再分析积雪为逐月数据， 因此分别计算对应年份冬春季平均再分析雪深和站点雪深的偏差和均方根误差， 再将多年的值平均得到最后结果。如图2所示， NOAA-V3和ERA5对高原中东部冬季积雪大部分站点表现为一定高估， 部分站点平均偏差大于15 cm。NOAA-V3在高原南部表现相对较好， 平均偏差小于5 cm， 除了高原中东部的部分站点外， ERA5在其他区域都表现较好。整体来看， ERA5雪深数据与站点观测雪深之间的偏差较NOAA-V3更小， 且部分站点的平均偏差都在1 cm以内。
NOAA-V3和ERA5再分析雪深的多年平均均方根误差如图3所示， 与平均偏差的空间分布类似， NOAA-V3和ERA5在高原中东部部分站点的均方根误差相对较大。NOAA-V3的平均均方根误差分布在整个高原范围内较为相近， 而ERA5的平均均方根误差在不同区域相差较大， 具体表现为在高原东南部部分站点较大而在高原北部和西南部较小。
上述偏差分析结果可见， 虽然NOAA-V3和ERA5再分析雪深数据与站点观测雪深在部分站点存在一定的偏差， 但大部分站点的平均偏差与平均均方根误差较小。其中ERA5雪深数据与站点观测雪深之间具有更好的一致性。ERA5数据虽然存在个别站点偏离程度较大的情况， 但整体来看其与站点雪深之间的平均偏差相对较小。NOAA-V3数据的平均偏差虽然存在普遍高估， 但整体来看各个站点较为稳定， 无偏离程度异常高的站点存在。
3.2 高原积雪年际及年代际变化特征
上文分析了NOAA-V3和ERA5再分析雪深数据和站点观测雪深间的偏差， 虽然部分站点存在一定的偏离， 但整体来看再分析数据与站点观测雪深之间存在较好的一致性。除了再分析雪深数据外， 为了检验积雪数据的不确定性， 本文还选取了对应的再分析雪盖数据以及卫星反演雪盖数据， 首先对站点观测数据和再分析积雪数据以及卫星反演积雪数据进行比较， 以检验不同积雪数据在描述年际及年代际变化特征方面的一致性。
图4（a）和（d）分别给出了高原站点观测雪深数据以及NOAA-V3和ERA5两种大气再分析积雪资料的雪深及雪盖数据标准化时间序列， 并定义站点观测雪深的标准化时间序列为高原冬春雪深指数SDI。计算结果显示， NOAA-V3再分析雪深和雪盖与高原50个站点雪深之间的相关系数分别为0.5和0.52； ERA5雪深和雪盖与站点观测雪深之间的相关关系分别为0.46和0.52， 两种大气再分析积雪数据和站点观测积雪的相关性均通过0.01的显著性统计检验， 这说明两种大气再分析积雪数据与站点积雪都呈现较好的对应关系， 表明站点雪深和再分析积雪数据年际变化特征较为一致。
图4 青藏高原站点雪深（黑色）、 NOAA-V3再分析雪深（红色）及雪盖（蓝色）年际变化时间序列（a~c）及青藏高原站点雪深（黑色）、 ERA5再分析雪深（红色）及雪盖（蓝色）年际变化时间序列（d~f）r 1表示站点积雪和再分析雪深的相关关系； r 2表示和雪盖的相关关系； 虚线为5年滑动平均时间序列 Fig.4 Interannual variation of snow depth observed by the 50 stations over the Qinghai-Xizang （Tibetan） Plateau （black） and snow depth by NOAA-V3 （red） snow cover by NOAA-V3 （blue） （a~c）， Interannual variation of snow depth observed by the 50 stations over the Qinghai-Xizang （Tibetan） Plateau （black） and snow depth by ERA5 （red） snow cover by ERA5 （blue） （d~f）， respectively.The data are standardized， where r 1 represents the correlation with snow depth and r 2 represents the correlation with snow cover， and the dotted line is the time series of five-year moving average |
分别对站点积雪数据和两种大气再分析积雪数据的时间序列做五年滑动平均处理以过滤高频振荡， 从而反映高原积雪的年代际变化特征。图4（b）、 （c）、 （e）、 （f）显示， 站点观测雪深在1996年前后出现相反的变化趋势， 在1996年以前高原雪深为增长趋势， 1996年之后转为减少趋势。NOAA-V3雪深在1994年之前呈现短暂的下降趋势随后转为上升趋势， 而1994年之后再次呈现下降趋势； NOAA-V3雪盖在1994年之前无明显上升或下降趋势， 而在1994年之后呈现下降趋势。ERA5雪深和雪盖的年代际变化趋势与站点观测雪深的年代际变化趋势较为一致， 两者均在1996年之前呈现上升趋势， 1996年之后呈现下降趋势。
比较NOAA-V3和EAR5两种再分析积雪数据和站点观测雪深的年际及年代际变化特征时， 我们选择的是再分析数据位于高原东部（28°N -38°N， 90°E -104°E）的平均值与高原站点观测数据平均值的对比。由于是格点数据和站点数据的对比， 可能会存在误差。为了减小误差， 分别筛选出再分析积雪数据对应每个站点位置周围的格点， 由于每个站点的经纬度不是整齐排列的数据， 因此分别选出每个站点对应经向和纬向两个方向的共4个格点， 总计200个格点的数值， 再将平均之后得到的数值与高原站点观测雪深平均后的值做比较。由此得到图5所示的时间序列。
图5 青藏高原站点雪深（黑色）、 NOAA-V3再分析雪深（红色）及雪盖（蓝色）年际变化时间序列（a~c）及青藏高原站点雪深（黑色）、 ERA5再分析雪深（红色）及雪盖（蓝色）年际变化时间序列（d~f）r 1表示站点积雪和再分析雪深的相关关系； r 2表示和雪盖的相关关系； 虚线为5年滑动平均时间序列； 再分析数据为对应站点周围格点的平均值 Fig.5 Interannual variation of snow depth observed by the 50 stations over the Qinghai-Xizang （Tibetan） Plateau （black） and snow depth by NOAA-V3 （red） snow cover by NOAA-V3 （blue） （a~c）， Interannual variation of snow depth observed by the 50 stations over the Qinghai-Xizang （Tibetan） Plateau （black） and snow depth by ERA5 （red） snow cover by ERA5 （blue） （d~f）， respectively.The data are standardized， where r 1 represents the correlation with snow depth and r 2 represents the correlation with snow cover， the dotted line is the time series of five-year moving average， and the reanalysis data is the average of the grid points around the corresponding site |
将再分析数据的数值由高原东部平均换为对应站点周围4个格点平均值后， 比较再分析积雪数据和站点观测雪深之间的相关关系发现， NOAA-V3雪深与站点观测雪深的相关系数从0.50变为0.47， 雪盖与站点雪深的相关系数从0.52变为0.53； 同样地， ERA5雪深与站点观测雪深的相关系数由0.46变为0.49， 雪盖与站点雪深的相关系数由0.52变为0.57。由此可见， 虽然上述相关分析结果都通过了0.01的显著性检验， 但进一步精确比较范围后， ERA5再分析积雪数据与站点雪深的相关性都相对增强， 而NOAA-V3积雪数据的变化不大。
同样对筛选出的再分析积雪数据做五年滑动平均处理来比较其年代际变化特征。与高原东部范围的积雪变化一致， NOAA-V3再分析雪深在1994年之前呈现先下降后上升趋势， 在1994年之后呈现下降趋势； 雪盖在1994年之后呈现下降趋势。ERA5雪深和雪盖都在1996年前后表现为相反的变化趋势， 1996年之前表现为上升而1996年之后为下降趋势， 此结果与站点雪深的年代际变化趋势一致。
通过对比高原东部再分析积雪的变化特征和对应站点的再分析积雪的变化特征发现， 两种方法处理的再分析积雪数据与站点观测积雪的相关性都较为一致， 且各自的年际及年代际变化趋势也大致相同， 此结果可以说明高原50站点积雪深度的变化在一定程度上能够代表高原东部积雪的变化。
车涛等（2004）研究指出， 青藏高原西部山脉崎岖， 站点分布较少， 卫星遥感手段被大量用于高原积雪的观测中。因此， 我们同样选取三种卫星反演雪盖数据来比较其与站点观测积雪数据的年际及年代际变化特征的一致性。由于卫星反演雪盖数据统一将有积雪覆盖的格点用特定的数字表示， 因此将整个区域内有积雪覆盖的格点累加作为该区域内的积雪覆盖的总格点数， 并且通过计算研究时段里目标区域内冬春季节有积雪的总格点数来表示高原冬春雪盖的年际变化序列。且卫星反演的积雪数据时间序列都相对较短， IMS提供的卫星雪盖数据的时段为1999 -2012年， MODIS传感器雪盖数据的时间范围为2001 -2012年， 只有北半球雪盖与本文研究时间段一致。
比较北半球雪盖数据和站点观测雪深时间序列时发现， 这两种积雪数据的变化趋势是大致相同的。但通过计算发现两者的相关系数仅为0.06， 结果显示站点雪深在1996年的标准差呈现正的最大值， 而北半球雪盖在1996年的标准差却为负值。于是对1996年的积雪数据进行剔除处理， 发现剩下的37年序列， 北半球雪盖与站点雪深的相关性由0.06提高为0.17。进一步将研究范围变为1980 - 1995年这16年， 发现北半球雪盖与站点雪深之间的相关系数为0.42， 而加上1996年后这17年间两者之间的相关系数减小为0.23。同样计算1997 - 2017年21年里北半球雪盖和站点雪深之间的相关关系， 得到两者间的相关系数为0.62（通过95%的信度检验）， 而加上1996年之后的22年里， 两者间的相关系数则降低到0.15。总体而言， 北半球雪盖与青藏高原站点观测雪深在不包含1996年的时间序列相关性更高（如表1粗体所示）。由此对比结果， 我们猜测高原50站点雪深数据在1996年的结果对整个时间序列的影响较大， 而卫星反演雪盖数据在1996年并未有明显的正异常甚至表现为负异常， 因此在包含了1996年的时间序列里， 卫星反演雪盖数据和站点雪深数据的相关性不强。建议在使用两种数据时应该尤其注意这一年的异常情况， 多方面进行对比力争得到最合理的结果。本文通过比较站点观测雪深数据与再分析积雪数据的结果， 以及与前人对积雪异常年份的判定结果确认1996年为典型积雪偏多年份。
表1 北半球雪盖与青藏高原站点观测雪深的相关关系Table 1 Correlation coefficient between northern hemisphere snow cover and snow depth at station on the eastern Plateau |
| 相关关系 | 高原50站点雪深 |
|---|---|
| 北半球雪盖（1980 -2017年） | 0.06 |
| 北半球雪盖（1980 -2017年）（剔除1996年） | 0.17 |
| 北半球雪盖（1980 -1996年） | 0.23 |
| 北半球雪盖（1980 -1995年） | 0.42 |
| 北半球雪盖（1996 -2017年） | 0.15 |
| 北半球雪盖（1997 -2017年） | 0.62 |
加粗数值表示表示不含1996年的时间序列 |
接下来， 我们比较三种卫星雪盖的年代际变化趋势， 由图7所示， IMS雪盖和MODIS雪盖在对应的时间段内无明显增加或者减少趋势， 相应的， 站点雪深在21世纪后也未表现出明显的变化趋势。而北半球雪盖在2000年之前表现为缓慢上升趋势， 2000年之后表现为下降趋势， 此种趋势在站点观测雪深的年代际变化中也能体现出来。
图7 IMS雪盖（a）、 MODIS传感器雪盖（b）、 北半球雪盖（c）以及青藏高原站点观测雪深滑动平均的年际变化时间序列对应各数据时间段， 虚线为五年滑动平均时间序列 Fig.7 Interannual variation of snow cover derived from IMS （a）， MODIS （b）， the Northern Hemisphere （c） and the five-year moving average snow depth observed by the 50 stations over the Qinghai-Xizang （Tibetan） Plateau.For their time periods， and the dotted line is the time series of five-year moving average |
3.3 积雪异常年选取结果比较
通过对高原上不同积雪数据年际及年代际变化特征的分析， 我们发现站点观测积雪数据和再分析积雪数据以及卫星反演积雪数据之间存在较好的对应关系， 尤其在年际变化特征方面呈现较高的一致性。该结果表明由站点观测雪深的年际变化时间序列来筛选积雪异常年份是可行的。
因此， 选取站点观测的1980 -2017年冬春积雪深度数据， 在多雪年和少雪年的选取数量相对较一致的情况下， 根据高原冬春雪深指数时间序列以0.9个单位为阈值判定多雪年和少雪年， 如图8所示， 多雪年为1988/1989年、 1992/1993年、 1994/1995年、 1995/1996年、 1997/1998年、 2007/2008年共6年； 少雪年为1983/1984年、 1998/1999年、 2002/2003年、 2009/2010年、 2015/2016年共5年。
在以往的研究中， 除了用雪深划分积雪异常年份外， Xiao and Duan（2016）还利用积雪日数划分了多雪年与少雪年， 以及以雪深和积雪日数的变化共同划分（罗江鑫等， 2020）。通过与前人划分的积雪异常年份的对比（陈乾金等， 2000； Xiao and Duan， 2016； 罗江鑫等， 2020）， 同样说明本文所选取的多雪年和少雪年都为典型积雪异常年份。
3.4 高原积雪空间分布特征
前文通过比较不同积雪数据的时间分布特征发现， 站点观测积雪数据与再分析积雪数据以及卫星反演积雪数据在年际变化特征上具有较好的一致性， 此结果也佐证了由站点观测积雪数据筛选的积雪异常年的可靠性。进一步比较各积雪数据的空间分布特征， 结果显示， NOAA-V3再分析雪深和ERA5雪深多年平均分布一致（图略）， 冬春季青藏高原雪深的大值区主要位于喜马拉雅山脉附近， 以及高原的中东部和北部。这与沈鎏澄等（2019）得到的结果一致。接着根据高原站点雪深挑选的积雪异常年份， 分别对比不同再分析数据在多雪年和少雪年的空间分布。如图9所示， NOAA-V3雪深与ERA5雪深在多雪年和少雪年的分布整体基本一致， 其中NOAA-V3雪深大值区位于高原中东部以及西北部昆仑山区域； ERA5雪深在积雪异常年份的大值区域主要为高原中东部以及南部喜马拉雅山脉区域。可以看出NOAA-V3雪深与ERA5雪深数据在表现高原雪深大值区的时候出现细微差异， 高原中东部的大值中心两者皆能体现， 而NOAA-V3雪深更加突出高原西北部昆仑山地区的积雪大值区， 而ERA5数据在表现高原南部喜马拉雅山脉地区雪深大值区的效果更为显著。
如图9（c）、 （f）所示， 在多雪年和少雪年的差值分布中， NOAA-V3再分析雪深增加的大值区主要位于高原中部地区， 高原北部以及东南部在多雪年主要为减少趋势； ERA-5再分析雪深在高原中部以及南部喜马拉雅山脉地区为积雪增加的大值区， 高原东南部为积雪减少的大值区。
分析了NOAA-V3和EAR5再分析雪深数据在积雪异常年份的分布及其差值分布， 发现两者对于青藏高原中东部雪深大值区的表征都较为显著， 差异在于NOAA-V3雪深数据更好地凸显出高原西北部昆仑山地区的积雪深度大值区， 但对高原南部喜马拉雅山脉地区的雪深大值区未能较好体现； 相反， ERA5雪深数据能较好地体现出高原南面喜马拉雅山脉附近积雪深度的大值区， 但未能较好地体现出高原西北部昆仑山地区的积雪。两者在积雪异常年份的差值分布大致相同， 在多雪年， 高原中部都表现为显著的雪深增加的大值区， 高原北部及东部为雪深减少的大值区。
同样比较NOAA-V3和ERA5再分析雪盖数据， 结果显示两种再分析雪盖数据的多年平均空间分布与对应的雪深分布基本一致（图略）， 雪盖的大值区分布在高原西部、 中东部， 以及北部部分区域。
如图10所示， 不同再分析雪盖数据与雪深数据在积雪异常年份的空间分布大致相同。NOAA-V3雪盖的大值区主要位于高原中东部以及西北部昆仑山脉地区， 此结果与NOAA-V3雪深的分布一致； ERA5雪盖的大值区主要位于高原中东部和南部喜马拉雅山脉地区， 同样与ERA5雪深的分布一致。另外， NOAA-V3雪盖在多雪年主要在高原中部以及高原东南角和东北角呈增加趋势， ERA5雪盖除了高原中部、 东南以及东北角外， 高原的南部也为雪盖增加的大值区。而两种再分析雪盖在高原的北部边缘以及高原中东部的变化也较为一致， 都表现出了高原雪盖在多雪年的负异常变化。
整体上看， NOAA-V3和ERA5两种再分析数据都能较好地体现出积雪异常年份的空间分布特征。无论是雪深还是雪盖， 两种再分析数据在积雪异常偏多年对应增减的区域较为一致， 都能较好地表征高原中东部的积雪大值区， 说明站点雪深数据和选取的两种再分析积雪数据在一定程度上是相对应的， 进一步证实了由站点雪深数据筛选出的积雪异常年份的典型性。但也应该看到， NOAA-V3和ERA5再分析积雪数据在积雪异常年份的积雪分布存在差异， 具体来说， NOAA-V3对于体现高原西北部昆仑山地区的积雪大值区有较好的效果， 而ERA5更能较好地表征高原南部喜马拉雅山脉地区的积雪大值区， 此种差异可能与不同再分析数据的处理方式以及原始数据的选取不同有关。除此之外两者之间的空间分布大致相同。
图11给出了罗格斯大学提供的北半球雪盖在积雪异常年份的分布及其差值分布、 MODIS传感器雪盖在积雪异常年份的分布及其差值分布以及IMS逐日雪盖在积雪异常年份的分布及其差值分布。结果显示， 三种卫星反演雪盖在多雪年和少雪年的分布大致相同， 积雪大值区主要位于高原西南部喜马拉雅山脉附近， 以及高原中东部和北部山区， 该结果与再分析积雪数据的分布也基本一致。
三种雪盖数据在多雪年和少雪年的差值分布也大致相同， 负异常的大值区均位于高原东南部， 其余大部分地区表现为正异常。多雪年， 北半球雪盖除了在高原东南缘减少外， 高原北部也有一定的减少， 这与再分析积雪数据的变化较一致。MODIS传感器雪盖和IMS雪盖在多雪年和少雪年的差值分布基本一致， 除了高原东南沿雪盖减少外， 其余大部分地区雪盖增加， 其中高原中东部和东北部为雪盖增加的大值区。
总体来说， 卫星反演雪盖数据在积雪异常年份的积雪增加区域与再分析数据的增加区域大致相同， 高原的中东部及南部边缘都为积雪增加的大值区， 但多雪年再分析数据在高原北部边缘体现出的负异常却未能在MODIS雪盖和IMS雪盖数据上体现， 这可能与不同积雪数据的获取以及处理方式不同有关， 具体原因有待更进一步研究。
4 结论
本文利用青藏高原台站观测、 再分析（ERA5 和NOAA-V3）和卫星反演（MODIS雪盖以及IMS雪盖）的多源积雪资料， 采用偏差分析、 均方根误差以及相关分析等多元统计方法， 通过比较不同积雪资料的时空分布和变化特征， 重点检验了多源高原积雪数据在描述积雪年际变化特征方面的不确定性， 以期提升多源高原积雪资料适用性的认知， 并为相关研究提供有意义的参考。通过分析得到如下结论：
（1） 与NOAA-V3雪深相比， ERA5雪深对高原站点观测雪深的表现效果更好。除了高原中东部分站点外， ERA5雪深数据的平均偏差和平均均方根误差均较小， 而NOAA-V3雪深数据的平均偏差和均方根误差在整个高原范围内均存在一定程度的高估。
（2） NOAA-V3和ERA5再分析积雪数据和高原站点雪深在年际变化特征上具有较好的一致性； 长期变化趋势上， ERA5再分析积雪数据与高原站点观测雪深数据的变化较为类似； MODIS雪盖以及IMS雪盖在对应的年份与高原观测站点雪深的年际变化特征较为一致， 而北半球雪盖与站点观测雪深的年代际变化趋势更加类似。
（3） 同一种再分析资料的雪深和雪盖数据在高原区域分布差别较小， 而不同再分析资料间略有区别， NOAA-V3再分析雪深和雪盖在多雪年和少雪年的大值区都位于高原南部、 北部边缘以及高原中东部； 而ERA5再分析雪深和雪盖在多雪年和少雪年的大值区主要分布在高原的中东部、 西部以及南部边缘。两种再分析数据在积雪异常年份的差值分布较为一致， 无论是雪深还是雪盖， 在多雪年增加的大值区都位于高原中部、 西部和南部边缘， 而高原北部的负异常分别在两种再分析数据中都有体现。
（4） 三种卫星数据在积雪异常年份分布结果也大致相同， 雪盖在多雪年大值区都位于高原中东部、 南部以及西部地区。另外， 不同卫星数据在积雪异常增多年份雪盖的增加分布也较为一致， 但在再分析数据中高原北部边缘的负异常， 并未在卫星雪盖上体现出来， 此结果可能与卫星反演积雪数据的涵盖时间短以及处理方式不同有关。
通过对比不同积雪数据的时空变化特征， 我们对积雪数据的不确定性进行检验， 虽然站点观测、 再分析积雪数据以及卫星反演雪盖在描述高原积雪变化时存在一定的偏差， 但高原站点观测的雪深变化和其他积雪数据具有类似的变化趋势， 特别在年际变化尺度上具有较好的一致性， 进一步通过比较不同积雪数据在多雪年和少雪年的空间分布亦能说明各积雪数据整体的一致性特征。另外本文所选取的卫星反演积雪资料覆盖年限较短， 我们下一步工作将会选取更长时间的数据来进行不确定性检验以便结果更具有典型性和代表性。