杨映，唐忠喜，邢德，侯金亮

1.甘肃洮河国家级自然保护区管护中心，甘肃甘南藏族自治州 747600

2.甘肃省地质矿产勘查开发局第三地质矿产勘查院，兰州 730050

3.国防科技大学海洋气象学院，长沙 050024

4.中国科学院西北生态环境资源研究院甘肃省遥感重点实验室，兰州 730000

引 言

2019年9月18日，习近平总书记在主持召开黄河流域生态保护和高质量发展座谈会时，将黄河流域生态保护和高质量发展提升为重大国家战略，特别强调了黄河流域的水源涵养问题。黄河源区位于青藏高原东北部，平均海拔为4026 m，海拔3000 m以上的区域占整个研究区的95.80%，年平均径流量约为 204.0亿立方米，占黄河流域年径流量的三分之一，是黄河流域的主要产水区和水源涵养区，被誉为“黄河水塔”[1-2]。而青藏高原是北半球中低纬度海拔最高、积雪覆盖最大的地区，既是气候变化的敏感区，又对水资源系统产生重要影响[3-5]。已有研究表明，青藏高原东部的积雪变化最显著，其积雪主导了整个青藏高原积雪的变化[6-8]。因此，获取黄河源区准确的积雪覆盖变化信息，具有十分重要的意义。

卫星遥感是持续监测积雪时空变化的有力工具[9]。光学传感器利用积雪与云以及其他地物在可见光（VIS, 0.4-0.7 μm）和近红外（NIR, 1-4 μm）波段的较大反射率差异特性，在全球、区域尺度上的积雪监测中发挥了重要作用[10]。目前，国内外学者利用甚高分辨率辐射计（Advanced Very High Resolution Radiometer, AVHRR）、中分辨率成像光谱仪（Moderate Resolution Imaging Spectroradiometer,MODIS）、可见光红外成像辐射仪（Visible infrared Imaging Radiometer Suite, VIIRS）和风云三号气象卫星（FY-3）上搭载的中分辨率光谱成像仪（Medium Resolution Spectral Imager, MERSI）等光学传感器研制了不同时间分辨率和空间分辨率的多种全球和区域尺度的积雪覆盖产品[11-14]。这些积雪产品不仅可以用于监测积雪覆盖的时空变化，还可以作为各种大气模式和陆面/水文过程模型的输入参数。其中，AVHRR和FY-3积雪面积产品的空间分辨率较粗，在积雪斑块化分布、消融变化剧烈以及地形复杂山区等情况下，难以精准反映积雪的空间分布细节特征。VIIRS和MODIS积雪产品非常相似，但VIIRS时间序列较短（2012年起）。简而言之，已有的各种积雪覆盖产品中，MODIS系列产品因其较优的质量、较高的时空分辨率（可达到500m，逐日）、全球覆盖、长时间序列（2000年起）、免费使用等优势[9,15-16]，在积雪研究领域应用最广。全球不同区域上的大量验证结果表明，MODIS积雪产品具有很好的精度，晴空条件下的总体精度一般在 90%-96%之间[17-19]。然而，由于云覆盖的影响，导致MODIS积雪产品中存在大量的数据缺失。已有研究表明，逐日MODIS积雪产品中由于云覆盖造成的数据缺失在 39%-70%之间[20-21]。在所有天气条件下，MODIS积雪产品的总体平均精度仅为34%-50%[22]。显然，MODIS积雪产品中的数据缺失严重阻碍了其在水文模拟、气候变化以及雪灾预警预报等方面的应用。

为提高 MODIS积雪产品的监测精度，大量学者对 MODIS积雪产品中的云覆盖进行了重建研究，如黄晓东等[23]、YU等[24]和邱玉宝等[25]均研制了青藏高原区域的逐日无云二值积雪产品，HAO等[26]基于MODIS发展了中国逐日无云二值积雪产品。唐志光等[27]发展了青藏高原区域的MODIS积雪面积比例（FSC）产品。自2016年起，新版本的MODIS积雪产品不再直接提供二值积雪产品和FSC产品，而是提供归一化积雪指数（NDSI）产品，因为利用NDSI值可以很方便地推算出二值积雪覆盖和FSC值[13]。目前，已有研究基于相似像元/块的思想，对MODIS NDSI积雪产品中的云覆盖进行了重建[28-30]，但受青藏高原云系旺盛、积雪斑块化严重以及积雪消融较快等因素影响，在该区域搜索到合适的相似像元/块并非易事，有时甚至根本就不存在合理的相似像元。XING等[31]利用深度学习技术深入挖掘时空特征，充分利用MODIS NDSI积雪产品中可用的时空信息，构建了基于部分卷积神经网络（Partial Convolution Neural Network，PCNN）的MODIS NDSI积雪产品云像元重建模型，并在黄河源区进行了验证，结果表明PCNN模型重建后的MODIS NDSI产品达到了可靠的较高精度。

本数据集利用2000-2021年逐日的MODIS NDSI产品，使用上述基于PCNN的MODIS NDSI云像元重建模型，生成了时空连续的MODIS NDSI数据产品；其次，再利用NASA FSC产品的标准算法，制备了黄河源区2000-2021年无云逐日MODIS FSC遥感监测数据集，将为黄河源区的积雪分布、雪水储量估计、积雪变化分析和雪灾风险评估等研究工作提供数据支撑。

1 数据采集和处理方法

1.1 数据生产流程

黄河源区MODIS FSC数据集的制作流程主要分为5个步骤：获取MODIS积雪产品、数据预处理、MODIS NDSI云像元重建（包括：MODIS/Terra和MODIS/Aqua数据合成、临近三天时间滤波和基于PCNN的云像元时空重建）、FSC估计和精度验证，具体流程如图1所示。

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

1.2 数据获取

本研究使用的是MODIS/Terra和MODIS/Aqua逐日500 m积雪覆盖L3级产品（MOD10A1和MYD10A1, V6）。从NASA的雪冰数据中心（https://nsidc.org/）免费下载了研究区2000-2021年的数据。覆盖整个黄河源区需要二个瓦片（即：h25v05和h26v06），影像数据均为HDF文件存储格式、Sinusoidal投影。

雪深验证数据来自中国气象局（http://data.cma.cn）地面台站观测的2000-2020年每日地面气候积雪资料数据集，雪深观测以厘米（cm）为单位，取整数，雪深小于1 cm的记录为无雪。黄河源区上6个气象站记录的雪深资料被用来对积雪覆盖比例数据集进行基于离散值的精度评估。

1.3 数据预处理

利用MODIS MRT批处理工具对两个瓦片进行拼接、重投影、坐标转换、重采样等相关预处理，将原始数据转换为地理坐标系，0.005°空间分辨率。本研究选取Raw NDSI和NDSI snow cover作为研究数据。在Raw NDSI的值在-10000-10000，转换系数为0.0001。NDSI snow cover数据的编码为整数（0，10-100），其他类型包括云、内陆水体、海洋、缺失值、无决策、夜间和传感器饱和数据。由于通过Raw NDSI的值无法判断出一个像元是否为有效像元，本研究利用NDSI snow cover数据中的分类，将NDSI snow cover数据中所有其他类型都重分类为云，再使用重分类后的NDSI snow cover数据对Raw NDSI数据进行云掩膜，并将Raw NDSI实际值小于0的像元的NDSI值设置为0，将处理后的Raw NDSI数据作为基础数据层。

1.4 MODIS NDSI积雪产品云像元重建

1.4.1 MOD10A1和MYD10A1数据的逐日合成

根据云会移动的特点，对上午星的Terra/MODIS和下午星的Aqua/MODIS的NDSI数据进行逐日合成处理，采用“最大化”合成准则，具体是将一日内两幅NDSI影像进行对比，每个像元NDSI值取两者中的较大值，合成后的数据简称为MOYD产品。

1.4.2 邻近时间滤波

根据积雪降落之后，会在地表保留一段时间的特点，对逐日合成后的MOYD积雪产品进行时间上的滤波分析。考虑到黄河源区存在大量斑块状积雪、消融变化快的特点，选择最严格的时间窗口限制，即邻近三天时间滤波。具体方法是：当一个像元被云覆盖时，若该像元前一天和后一天均未被云覆盖，取前、后两天NDSI的平均值，作为该云像元的NDSI值。

1.4.3 基于PCNN的云像元重建

使用XING等[31]提出的基于端到端深度学习的三维PCNN的NDSI预测模型，利用时空信息估计云像元的地表覆盖信息。首先，选择块的大小为 64×64，在研究区上遍历，寻找不包含云的完整NDSI块，研究区上2000-2021年间21个积雪季（即前一年11月1日至后一年4月30日）上共搜索到103630个完整的块。其次，随机选择103630个有云的块（云覆盖率0%-10%，……，90%-100%各占1/10），以七天为时间窗口（即前、后各三天），提取完整块七天的数据，形成时空信息块组，并用对应有云块组对完整块组进行云掩膜。再次，将103 630个块组，按7:3划分为训练集和测试集（即72 541个训练块组，31 089个测试块组），在训练集上训练PCNN模型，独立测试集上的数据用于验证模型的性能。最后，将所构建的模型应用到整个黄河源区上，重建所有缺失像元的NDSI值，得到黄河源区2000-2021年时空连续的NDSI数据集。

1.5 MODIS FSC估计

使用NASA第五版MODIS FSC的标准算法[15]，根据像元的NDSI值推算其FSC值，公式如下：

根据公式（1）分别利用Terra/MODIS原始NDSI产品、Aqua/MODIS原始NDSI产品、双星合成后的NDSI数据、临近三天时间滤波的NDSI数据以及PCNN模型重建的时空连续NDSI数据计算FSC值，得到对应的FSC产品，分别简称为MOD、MYD、MOYD、ATF和PCNN产品。

1.6 精度验证

1.6.1 基于站点雪深观测值的验证

利用黄河源区上6个气象观测站点（玛多、达日、河南、久治、若尔盖、红原）2000年1月1日至2020年3月31日观测的雪深值为“真值”，对计算出的MODIS FSC进行评估。提取6个站点对应像元的MODIS FSC值，与站点实测的SD值进行比较，构建混淆矩阵，如表1。

表1 MODIS FSC和站点观测雪深的对比混淆矩阵Table 1 Comparative confusion matrix between MODIS FSC value and snow depth observations from stations

3个验证指标总体精度（OA）、高估积雪事件（MO）和低估积雪事件（MU）分别定义如下：

1.6.2 基于云假设的验证

本研究进一步基于云假设，即用“假想”的云覆盖，对原始晴空 NDSI影像进行掩膜，再利用PCNN方法对云像元重建后，对FSC的估计结果进行精确的定量评价，定义两个定量评价指标平均绝对偏差（MAE）和相关系数（R）如下：

其中，N是总的云像元数；xis是第i个假设“云”像元FSC的估计值；xiT是对应的参考“真实”FSC值；̅和̅是所有xis和xiT的平均值。

2 数据样本描述

黄河源区2000-2021年无云积雪覆盖比例数据集包含从2000-2001积雪季开始（即2000年11月1日）至2020-2021积雪季结束（即2021年4月30日）的7485个逐日的栅格数据文件。为便于数据在通用遥感软件下处理与展示，将数据存储为tif格式，投影方式为WGS84地理坐标系，文件命名规则为：YYYYDDD.tif，其中 YYYY代表年，DDD代表儒略日（001-365），空间分辨率为0.005°，属性值FSC为像元上积雪覆盖的百分比，覆盖范围为整个黄河源区。积雪覆盖比例产品值的取值范围均为0-100，填充值为200（即黄河源区范围之外像元上的值），单位：无。注意，由于MODIS 原始数据的缺失（2000301、2001166-2001184、2001221、2001310、2001339、2002079-2002087以及2002105共33天无数据），故2000年11月1日至2021年4月30日的7518天中，共有7485天有数据。

图2展示了2018-2019积雪季上6个样例日MODIS FSC的空间分布状况。可以看出，本无云MODIS FSC数据集从视觉上看，积雪分布变化平滑，黄河源区的积雪主要分布在中部的阿尼玛青山和南部的巴彦喀拉尔山主峰；同时，黄河源区上分布着大量斑块状的较浅积雪，这与黄河源区上的实际情况相符合。积雪季上不同时间的 FSC空间分布很好地展示出了积雪的积累-消融变化过程，一些积雪事件（如2018年11月15日以及2019年2月底的新降雪），MODIS FSC产品都有很好的及时响应。

图2 MODIS FSC分布图。(a) 2018年10月15日；(b) 2018年11月15日；(c) 2018年12月15日；(d) 2019年1月15日；(e) 2019年2月15日；(f) 2019年4月15日Figure 2 MODIS FSC distribution map.(a) October 15, 2018; (b) November 15, 2018; (c) December 15, 2018; (d)January 15, 2019; (e) January 15, 2019; (f) April 15, 2019

3 数据质量控制和评估

3.1 基于站点雪深观测值的评估

基于站点雪深观测的验证，是最直接的验证方式。利用站点实测雪深值的验证结果如表2所示。表中平均云覆盖日数百分比表示6个站点平均云覆盖日数占总日数的百分比。可以看出，由于云覆盖造成的严重数据缺失，导致MODIS原始产品超过一半的时间被云覆盖。双星合成和临近时间窗口滤波，均可以将 MODIS FSC产品的云覆盖降低 10%左右，PCNN重建所有云像元后，时空连续MODIS FSC产品的总体精度显著提高，可以达到94%，并且跟原始产品晴空像元条件下比较，并没有显著增加高估积雪事件和低估积雪事件发生的概率（高估和低估均增加1%）。

表2 基于站点观测雪深值的验证结果Table 2 Performance evaluation results based on snow depth observations from stations

3.2 基于云假设的评估

受黄河源区地形复杂、气象站点稀少且站点主要分布在低海拔平地（只有玛多站位于海拔4000 m以上的区域）等因素影响，导致站点雪深对积雪信息描述存在较大不确定性。基于站点雪深值的验证，可能并不能准确代表整个研究区上的情况。基于云掩膜的验证，是对基于 PCNN的 MODIS FSC估计模型的精确定量评价，结果如图3所示。可以看出，PCNN模型重建的MODIS FSC值的平均MAE为10.43%，平均R约为0.93，模型在FSC值40%-70%时，性能稍差一些，主要是因为此时雪层较浅，积雪分布不稳定。

此外，我们还进一步选择了研究区上5天无云的影像（2011314、2013320、2014084、2016306、2017348）作为真值，再假设每景都被20%、50%、80%的云覆盖，使用本研究提出的基于PCNN的方法去云后，将去云的影像与原始真值进行了比较，性能定量评价结果如表3。以2017328这一天为例，原始影像、假设的云覆盖情况以及去云后的影像如图4所示。结果表明：PCNN方法的鲁棒性较好，即使云覆盖率达到80%，基于PNCC的云像元重建方法也可以有效利用时空信息，将云像元的积雪覆盖信息高质量地重建。

图3 基于云假设的性能定量评价结果。(a) MAE；(b) RFigure 3 Performance quantitative evaluation results based on cloud assumption.(a) MAE; (b) R

表3 5个晴空日上基于云假设的验证结果Table 3 Performance evaluation results based on cloud assumption for five selected cloudless days

图4 2017年12月14日（即2017328）的FSC分布图。(a1)真值；(b1)-(b3)分别指20%、50%、80%的云覆盖（白色区域）掩膜真值影像；(c1)-(c3) 分别指PCNN方法去云后的影像Figure 4 FSC distribution map on 14 December 2017 (i.e.2017328).(a1) Truth FSC image on December 14,2017; (b1)-(b3) indicate that the 20%, 50% and 80% of the original truth image are covered by artificial clouds (white area); (c1)~(c3) indicate the reconstructed FSC distribution map of (b1)-(b3) by PCNN method.

4 数据价值

本数据以 MODIS V6版积雪产品为基础，对 MODIS NDSI产品中的云像元重建（包括：MODIS/Terra和MODIS/Aqua数据合成、临近三天时间滤波和基于PCNN的云像元时空重建）的基础上，进一步计算对应的FSC，制备了2000-2021年黄河源区积雪覆盖比例数据集，可以为黄河源区积雪状况评估、积雪演变特征分析、水资源综合管理、雪灾风险评估及预测等研究工作提供数据支撑。