2002-2021年青海省积雪日数遥感监测数据集

李晓东，李林，史飞飞,4*，苏文将，肖建设，李红梅,5

1.青海省气象科学研究所，西宁 810008

2.青海省防灾减灾重点实验室，西宁 810008

3.青海省气象局，西宁 810008

4.青海师范大学地理科学学院，西宁 810008

5.青海省气候中心，西宁 810008

引 言

第六次IPCC评估报告指出，随着全球变暖的加剧，强降水、高温热浪等高影响气候事件将更为频繁[1]。被誉为世界“第三极”的青藏高原，对全球气候变化的反馈具有一定的“前瞻性”。近年来青藏高原部分地区异常降雪过程增多，积雪状况也发生了显著而又独特的变化[2-3]。青藏高原作为全球典型的高海拔积雪分布区以及中国和周边诸多河流的发源地，积雪消融及其变化状况对该地区乃至全球尺度的水循环、气候变化进程将产生较大影响，因而准确获积雪变化信息具有十分重要意义[4-6]。遥感具有宏观、高效和高频次等优势，能够快速实现积雪的动态监测。在近五十余年的积雪遥感监测过程中，国内外学者利用 SMMR、SSM、AMSR-E等微波传感器，以及 AVHRR、MODIS、Suomi NPP和FY-3 MERSI等光学传感器研制了多种全球和区域尺度的积雪覆盖、积雪深度和雪水当量产品[7-10]，但受青藏高原云系旺盛、积雪斑块化严重以及积雪消融较快等因素影响，现有的全球积雪产品精度在高原地区普遍不够理想[11]。为提高积雪产品在青藏高原的监测精度，黄晓东等[12]、于金媛等[13]和邱玉宝等[14]均基于MODIS研制了逐日无云积雪产品，唐志光等[15]发展了MODIS积雪面积比例产品，上述积雪产品的监测精度均得到较大程度提高并在积雪评估研究中得到广泛应用。

青藏高原雪层普遍偏浅且呈非连续性分布，在反映积雪类型分区和分布特征方面积雪日数相比其他积雪监测数据更具有代表性[16-17]，但目前积雪日数多以地面气象站的积雪观测资料为主，受青藏高原地形复杂和气象站点稀少等因素影响其对积雪信息描述存在较大不确定性。本数据集利用2002-2021年MODIS V6版积雪产品中的归一化积雪指数（Normalized Difference Snow Index，NDSI）等数据，通过逐日积雪覆盖提取、积雪覆盖范围校正、雪季内积雪日数合成、冰川区积雪日数修正和气象站点数据校正等步骤，制备了青海省积雪日数遥感监测数据集，将为青藏高原地区的积雪分布及雪水储量评估、积雪变化分析和雪灾风险区划等研究工作提供数据支撑。

1 数据采集和处理方法

1.1 数据生产流程

青海省积雪日数数据集的制作流程主要分为8个步骤：获取MODIS积雪产品、积雪覆盖提取、积雪覆盖范围校正（包括：上下午数据雪盖合成、临近三天雪盖合成）、雪季内积雪日数合成、冰川区积雪日数修正、基于气象站点的积雪日数校正、精度验证和构建积雪日数数据集，具体流程如图1所示。其中输入的原始数据为2002年10月至2021年5月的MODIS V6版积雪数据，精度验证采用地面气象台站的人工观测数据，最终输出各积雪季内青海省的积雪日数数据。

1.2 数据获取

积雪日数数据集的构建选用 MODIS积雪产品，该产品在青藏高原地区的积雪研究中广泛使用[4]，其500 m的空间分辨率以及逐日的观测频次能够满足高原地区积雪特性监测的实际需求。MODIS积雪产品已更新至第六版，本文结合Terra和Auqa双星积雪数据，其产品编号分别为MOD10A1和MYD10A1，在谷歌地球引擎（Google Earth Engine，GEE）中获取2002-2021年青海省的逐日积雪数据。

图1 数据处理流程Figure 1 Data processing flow

1.3 积雪覆盖提取方法

此前在MODIS V5版积雪产品中是利用NDSI进行积雪覆盖的二值判识，其阈值设定为0.4，但在MODIS V6版积雪产品中不再提供积雪覆盖的直接判识结果，转而提供NDSI指数。ZHANG等[18]和高杨等[19]研究指出在对青藏高原进行积雪覆盖判识时将NDSI的阈值设定为0.4会存在严重漏判，而将NDSI的阈值设定为0.1较0.4更为合理。本文在此研究基础上，通过前期对MODIS影像中积雪覆盖目视判识和进行野外观测对比，获取了青海省积雪判识时 NDSI的经验性阈值，并最终采用分段阈值的方法开展青海省积雪覆盖范围提取，具体方法见公式（1）。在GEE平台下通过编写脚本实现了逐日积雪覆盖提取，同时利用MODIS积雪产品中的云（250）、内陆水体（237）和无数据（200）等质量控制编码值对每日MOD/MYD10A1数据中的非积雪信息进行掩膜。

式中：B2为0.55-0.57 um可见光波段的反射率，B1为1.62-1.65 um短波红外波段的反射率，NDSIth为判定阈值，通过实验将青海省内10-12月的NDSIth设定为0.1，1-3月的NDSIth设定为0.13，4-5月的NDSIth设定为0.08。

1.4 积雪覆盖范围校正

上下午数据雪盖合成是指对MODIS/Terra上午星积雪数据（MOD10A1）判识的积雪覆盖范围和MODIS/Aqua下午星积雪数据（MYD10A1）判识的积雪覆盖范围进行最大合成，具体是将一日内两幅积雪判识后的二值栅格图像进行叠加，将每个栅格单元值取两幅积雪判识结果中的最大值，最后合成一日内具有最大积雪覆盖范围的二值数据。

临近三天雪盖合成是指在较短时间内积雪未发生快速消融，在当天积雪监测结果中有栅格像元受云影响无法判识是否存在积雪时，如前天和后天该像元上有积雪存在，则认定当天依旧有雪覆盖。

1.5 积雪日数合成

积雪覆盖日数合成是将每日的积雪覆盖提取结果依据使用要求进行旬/月/季/年或任意时段的合成，具体是将每日内积雪覆盖的二值栅格数据进行空间叠加，将每个栅格单元值进行累加后，便可生成统计时段内的累计积雪日数栅格数据。本数据集是按积雪季进行合成，其中当年积雪季一般从上年10月开始至当年5月结束。

1.6 冰川区积雪日数修正

青藏高原冰川分布区海拔高且常年有积雪覆盖，但上述区域经常受云遮挡导致部分地区积雪日数的遥感监测结果偏低，针对该问题选用中国第二次冰川数据集[20]获取青海省内冰川分布区，并对上述区域内的积雪日数进行校正，即整个雪季内冰川区均有积雪覆盖。

1.7 站点积雪日数校正

利用 MODIS积雪数据初步获取的积雪日数会受到云的干扰造成一定误差，可利用地面观测的积雪日数对其进行校正[21]。其中积雪日数观测数据来源于地面气象台站，从中国气象局综合气象信息共享平台（China Integrated Meteorological Information Sharing System，CIMISS）获取并已进行质量控制，观测员于每日8时定时开展观测，如在后期有降雪过程则在14时和20时进行补测，规定在积雪观测场内积雪覆盖范围超过一半，且多次测量的积雪深度均≥2 cm时则记录当天积雪日数为1。在校正方法的具体实施时，首先按3:2比例随机将全省气象台站划分为建模站点（31个）和验证站点（20个），其次需获取整个积雪季内（上年10月至当年5月）建模气象台站观测的积雪日数和与气象台站空间位置对应处MODIS提取的遥感积雪日数，对上述两组时序数据按公式(2)进行拟合，以获得校正系数a和b。最后当其拟合效果较好时则按照公式(3)对MODIS积雪数据提取的积雪日数进行逐栅格像元校正。

式中：X为利用MODIS积雪数据监测的积雪日数，Y为站点观测的积雪日数，单位均为天（d），a、b为校正系数，Y1为MODIS积雪数据监测的积雪日数，Y2为校正后的MODIS积雪数据监测的积雪日数。

1.8 精度验证

将积雪日数数据集与验证站点（20个）的积雪日数观测数据进行对比，以评价数据的准确性，按下列公式计算决定系数(R2)和相对误差(Yr)，一般R2越接近1，Yr越小则认为制备的积雪日数数据精度较高。

式中：Xi为站点积雪覆盖日数，Yi为遥感监测的积雪覆盖日数，Y̅为遥感监测的积雪覆盖日数平均值，X̅为站点积雪覆盖日数的平均值，n为验证站点个数。

2 数据样本描述

青海省积雪日数遥感监测数据集共由19个雪季年的栅格数据文件组成，为便于数据在通用遥感软件下处理与展示，将数据存储为tif格式，其中文件名以雪季年（上年10月至当年5月）命名，空间分辨率为500 m，属性值为雪季内的积雪日数，单位为天（d），覆盖范围为青海省，图2展示了2017和2019雪季内青海省积雪日数的空间分布状况。

图2 2017雪季（a）和2019雪季（b）积雪日数空间分布图Figure 2 The spatial distribution of snow cover days in 2017 (a) and 2019 (b) snow season

3 数据质量控制和评估

以参与验证的地面台站人工观测的积雪日数作为“真值”，逐积雪季对本文制作的积雪日数数据集进行精度检验，结果表明在各个雪季内遥感提取的积雪日数与台站观测相关程度较高，其R2在0.75-0.93，平均相对误差25%。其中2017和2019雪季的积雪日数及其精度状况如图2、图3所示，积雪日数高值区主要分布在三江源腹地、可可西里地区和祁连山地区，而积雪日数低值区主要分布在柴达木盆地，其较好地展示了积雪在空间尺度上的分布差异，并与实际积雪分布现状较为相符；同时，针对异常降雪年份例如2019年，强降雪日数异常偏多、气温普遍较低、积雪持续时间长，在青海省内玛多县、杂多县、甘德县和玉树市等地雪灾等级为历史最重，本文积雪日数数据也能较好刻画积雪覆盖在空间和时间尺度上的积累状况，并对雪灾发生的风险也具有很好的指示意义。

4 数据价值

本数据以MODIS V6版积雪产品为基础，通过逐日积雪覆盖提取、积雪覆盖范围校正、雪季内积雪日数合成和气象站点数据校正等步骤制作了2003-2021雪季内青海省的积雪日数数据集，将对开展青藏高原积雪状况评估、积雪演变特征分析、雪灾风险评估及预测等研究工作提供数据支撑。

图3 2017雪季（a）和2019雪季（b）MODIS积雪日数与站点积雪日数相关性图Figure 3 The diagram of correlation between snow cover days in MODIS and site snow days in 2017 (a) and 2019 snow seasons (b)

5 数据使用方法和建议

本数据集为tif文件格式，可使用ArcGIS和ENVI等图像处理软件进行查看、统计和分析等，各数据以积雪季命名，数据起始时间为2003雪季（即2002年10月至2003年5月），数据结束时间为2021雪季（即2020年10月至2021年5月）。