韩超，沈言龙，欧阳志棋，谢佩瑶，郭慧，陈思勇，王晓艳*

1.兰州大学资源环境学院，兰州 730000

2.南京大学地理与海洋科学学院，南京 210008

引 言

积雪是地表重要的覆盖物，也是重要的淡水资源[1]，其较高的反照率和较低的导热率特性，对地表辐射平衡起到了重要的影响，同时积雪时空变化对全球生态系统及区域气候的变化起到至关重要的作用[2-7]。中国积雪主要分布于新疆，青藏高原和东北地区[2]。其中东北地区是中国重要的粮食生产基地。积雪对春季土壤墒情和作物生长的影响十分重要，会直接影响作物实际产量[8]。因而，研究东北地区积雪分布及其时空变化特性具有重要意义。

随着遥感和地理信息技术的发展，人们逐步开展了一系列有关积雪产品的研发。MODIS积雪产品具有较高的时间分辨率，其空间分辨率为500 m，近年来被广泛应用于积雪时空变化分析研究。其中，归一化差值积雪指数（Normalized Difference Snow Index，NDSI）在MODIS积雪面积制图中被广泛应用[9-12]。我国东北地区地处高纬度，地形地貌相对复杂，森林分布广泛[13-15]。受冬季太阳光照条件以及林地冠层的影响，当前 MODIS积雪产品 V6版本在该地区存在过度云掩膜问题，MOD/MYD10A1的NDSI_Snow_Cover数据层存在大量的数据空缺，尤其在大兴安岭及小兴安岭森林地区该问题较为严重。本文对位于大兴安岭地区的一个站点（编号为50247，经纬度为123.5579°E，52.0285°N）进行了云覆盖累计天数统计，发现该站点所处位置在MODIS积雪产品中月内云覆盖累计天数在1、2、3、11和12月份都高达26天以上。这导致该地区整个积雪季可用的NDSI数据极少，很大原因是MODIS积雪产品中将积雪林地错分为了云，从而限制了MODIS积雪产品在中国东北地区的应用[16]。当前的 MODIS积雪产品去云算法主要采用时空邻域的像元信息进行插值填补云下像元，这些方法对于大范围连续的云覆盖去云效果并不理想[17-23]。CHEN等[24]提出的时空自适应的方法可以生成时空连续的 NDSI序列，但是当存在过度云掩膜的情况下，该算法的效率和精度都会大大降低。另外，由于MODIS积雪产品中采用的low_NDSI screen，使得许多具有较低NDSI的森林积雪像元被判为无雪像元，造成森林积雪的漏分。

本文以MODIS V6版本MOD/MYD10A1和MOD09GA数据集为数据源，采用上下午星合成、决策树判断和时空自适应去云算法，显著提高了产品精度和计算效率，生成了2000-2020年中国东北地区积雪季逐日无云 NDSI序列产品。本数据可用于支持中国东北地区积雪分布时空变化、积雪物候等后续科研工作。

1 数据采集和处理方法

1.1 数据源

在MODIS V6积雪产品中，MOD10A1/MYD10A1产品包括Raw_NDSI、NDSI_snow_cover等数据集。两个数据层分辨率一致，像元位置一一对应，其中，Raw_NDSI数据集的数值范围为-10000～10000，该值是每个像元真实NDSI值的10000倍。NDSI_snow_cover数据集在Raw_NDSI数据集的基础上进行了属性分类，属性编码内容如下表所示[25]。NDSI_snow_cover数据集中10-100的数值代表可能的积雪像元，该值是像元真实NDSI值的100倍。MOD09GA产品包括地面7个波段的地表反射率数据，本研究选取了MOD/MYD10A1中的Raw_NDSI和NDSI_snow_cover以及MOD09GA中的 sur_refl_bo4(band4)作为数据源（表1）。上述产品数据来源于 Google Earth Engine（https://code.eartengine.google.com）。空间分辨率为500 m，时间分辨率为1 d，数据格式为*.tiff。

表1 MODIS积雪产品V6版属性定义Table1 Attribute definitions of MODIS Snow Product V6

1.2 数据处理方法

1.2.1 影像裁剪

在Google Earth Engine上首先对研究区积雪产品数据集进行裁剪。由于所需数据集的空间分辨率都是500 m，故不用重采样。数据统一处理后的空间范围为北纬37°-54°30′，东经114°-136°之间，范围包含中国东北积雪区（如图1），地理坐标系为WGS84。

图1 b2、b4和b6波段合成结果（2019-02-12）Figure 1 Synthetic results of Band2, Band4 and Band6（Feb.12, 2019）

1.2.2 影像合成

结合MOD/MYD10A1中的Raw_NDSI和NDSI_snow_cover两个数据集，进行上下午星合成，填充一些云像元，恢复其NDSI值[26]。算法规则为若原产品上午定义有云，下午定义无云，则将下午星NDSI值赋值给当天；若上下午星都有云，则未被恢复，依然保留该像元为云像元；其余情况则采用上午星的NDSI值。经过此步骤处理，可以一定程度地减小云覆盖比例。以2018年为例，平均每景影像中云量占比由45%降至36%。

1.2.3 决策树判别

采用决策树判别方法，恢复MOD10A1中错分为云的像元NDSI值。此步骤只针对MOD10A1中的云像元进行操作，决策树流程见图2。研究发现，在有雪季节，东北地区的云雪混淆现象非常明显，主要表现为将积雪或者森林积雪错分为云。故此步骤针对每年11月至次年2月的积雪季数据实施。NSIDC（National Snow and Ice Data Center，美国国家冰雪数据中心）发布的MODIS全球积雪面积产品采用的NDSI阈值为0.4，可以有效进行云雪的区分[27-28]。本文对东北地区大量的云样本进行了统计，结果表明云像元NDSI最大值不超过0.4。因此，首先对Raw_NDSI进行均值滤波，算法规则是像元周围9个像元的值直接平均赋值给中心像元，得到新的数据层Mean_NDSI。若MOD10A1中的云像元对应的Mean_NDSI大于0.4，则认为该像元不是云像元，恢复该像元的Raw_NDSI值。

图2 决策树流程图Figure 2 The flow chart of decision tree

当 NDSI在0-0.4之间的像元，主要包含了森林积雪、积雪边界和云像元。在假彩色合成影像（RGB=b2/b4/b6）中，积雪、林地积雪和云表现出不同的色彩，可以通过目视解译直接分辨，如图1所示。本文分别选取森林积雪像元和云像元进行统计分析，发现二者第四波段地表反射率（sur_refl_bo4），即绿光波段具有明显的差异。森林积雪像元在该波段值地表反射率大多分布在0-0.4之间；而云像元在绿光波段的反射率大于0.4，如图3所示。故采用0.4作为阈值进行判别，若绿光波段地表反射率小于或等于0.4，则该像元为非云像元，恢复其Raw NDSI值；若该波段地表反射率大于0.4，则该像元仍标记为云。统计2018年数据，经过此步骤图中云量占比大约降至20%。此步骤消除了东北地区积雪产品中的过度云掩膜现象，可以大大提高后续计算效率，而且恢复的Raw NDSI为像元的NDSI真值，保证了生成的时空连续的NDSI序列的精度。

图3 不同地表下b4波段地表反射率频率曲线分布图Figure 3 The curve distribution diagram of Band 4 reflectance frequency under different land surfaces

1.2.4 临近日合成

一般来说，积雪的累积和消融会持续一段时间，故可以采用临近日合成的方法，填补部分云下NDSI值。为了保证精度，本文只采用前后两天的数据作为填补数据[30,31]。算法规则为若目标像元前后天都是非云像元，则将前后天 NDSI值平均后赋值给当天；若前后天中有一天为非云像元，则直接使用该非云像元的NDSI值赋值给当天；若前后天都为云像元，则需寻找前后两天像元继续判别，算法规则同上；若前后两天同样也都是云像元，则目标像元仍未得到恢复，计算结束，等待下一步去云操作。

1.2.5 时空自适应去云

采用CHEN等[24]提出的时空自适应去云算法，进行全图层彻底去云。其主要思路是判断相似像元，基于相似像元的NDSI在自适应时空中赋予相应的权重去云。

其中，M(x0,y0,t0)是t0时刻待填充像元NDSI值，M(xk,yk,t0)是t0时刻相似像元NDSI值，M(x0,y0,t0+s)和M(xk,yk,t0+s)是t0+s时刻M(x0,y0)目标像元和M(xk,yk)相似像元的NDSI值。Δt0和Δ(t0+s)表示两个时刻相似像元和目标像元各自的差异变化。在一定时期内，我们可以认为两者差异变化相同，则得到公式（3）。

基于上述理论，通过时空相邻相似像元的NDSI值插补云像元，可将云像元NDSI值得到恢复。在w×w局部窗口中选择前n个NDSI差异最小的像元作为相似像元，其中相似像元差异定义为公式（4），Sk越小，表明两者差异越小。公式（5）是填充目标云像元的加权平均计算公式，公式（6）则为相似像元相应的求权重公式，此处采用的是反距离权重算法。公式（7）中，Dk表示目标像元和相似像元之间的相对距离，其中将Dk归一化到1-w/2（w是窗口大小）。本论文的参数设置为n=20，w=15。

1.2.6 精度评价

采用云假设方法对最终去云结果进行精度评价。未被云层覆盖的真实地表 NDSI值是已知的，故可以作为验证数据。使用云假设将验证样本进行云掩膜并重新计算，对其计算得到的数据与验证数据进行比较，评估去云算法的性能和数据产品质量。本文采用2018年作为验证数据集，将积雪季各月云覆盖率最小的图像作为真实影像，使用各月云覆盖率最大的图像作为掩膜影像进行验证，通过相关系数（r）、均方根误差（RMSE）和平均绝对误差（MAE）来评估产品数据的精度[24]。

2 数据样本描述

本数据集包含2000-2020年中国东北地区积雪季（每年9月1日至次年4月30日）逐日无云的NDSI时间序列产品，均为栅格数据，空间分辨率为500 m。文件并命名为ChinaNE_ndsi_xxxxxxxx.tif（其中xxxxxxxx表示日期，例如20200101是指2020年1月1日）。其中，数据储存格式为双精度类型。数据集中数值范围为-1～1，该值为去云后每个像元的NDSI值。图4给出了数据集中2018年各月15日的东北积雪区无云NDSI的分布。需要注意的是，数据集中未对水体进行掩膜，因此在使用该产品进行积雪制图时应考虑水体要素。

3 数据质量控制和评估

MODIS V6版本从数据生产的每个环节进行质量控制，与之前V5相比，有了很大的改进，对水、云、气溶胶等进行了大气校正[25,31]。因此，本研究数据来源可靠。

从原始数据到最终产品的计算中，每个处理步骤都确保数据准确性。上下午星合成去云是积雪产品去云常用的方法，主要是考虑同一天的地表具有稳定性；决策树中的阈值也是基于文献或者统计特征选取，并且选取了大量执行决策树之后的结果与假彩色合成影像进行了对比，该步骤可以有效恢复NDSI_snow_cover数据集中一些由森林积雪误分成云的问题。最终采用CHEN等的方法进行完全去云，该方法在原文中经过云假设检验证明了其精度，去云结果与原始数据的平均相关系数r为0.94，平均RMSE为0.12，平均MAE为0.09[24]。本文对最终产品也经过云假设检验，平均相关系数r、平均RMSE和平均MAE分别为0.95、0.10和0.06。最后对2018年逐日的数据产品均进行了目视检验，结果表明该数据集NDSI分布合理，可以作为东北地区积雪研究的基础数据源。

图4 积雪年2018年各月15日NDSI产品展示图Figure 4 NDSI product display on the 15th day of each month in 2018 Snow Cover year

4 数据价值

积雪是全球气候变化的指示器，是全球气候变化研究的关键变量，因此探究积雪时空变化受到众多学者的广泛关注。在现阶段遥感积雪信息提取中，仍主要使用的是归一化差值积雪指数（NDSI）方法。但当前MODIS 积雪产品V6版本在东北地区，尤其是在东北森林地区存在过度云掩膜，使该地区有效的NDSI指数存在大范围和长时间的数据值空缺，严重影响了MODIS积雪产品在该地区的实际应用。本文针对MOD/MYD10A1的Raw_NDSI数据层中大量云像元带来的数据空缺，采用多种有效的方法进行数据恢复，生产了2000-2020年中国东北地区500 m分辨率积雪季逐日无云时间序列 NDSI数据。较已有的积雪产品，本数据集在去云和林区积雪识别方面有更好的精度，且有较长的时间跨度和较高的时间连续性，可为后续积雪二值提取，积雪混合像元分解等后续积雪产品的研究和生产工作提供数据源支撑。同时本数据集对东北地区积雪时空变化、森林水文过程、农业生态、碳循环过程和气候环境变化等研究具有重要的科学意义。

5 数据使用方法和建议

中国东北长时间序列逐日无云归一化差值积雪指数（NDSI）数据集，保存为 TIFF格式，能够在ArcGIS、ENVI及MATLAB等相关软件中对数据进行读取、查看、编辑及后续的一系列统计分析工作。需要说明的是本数据集并未对水体进行掩膜。此外，考虑到本数据集服务于后续的积雪产品生产工作，本文保留了东北地区边界外的数值，即并未使用东北积雪区矢量边界对数据集进行裁剪工作。数据集中图像所覆盖的范围完全包含了东北积雪区。由于本文采用的自适应时空去云算法，需要用到前后日期的影像数据，因此导致本数据集中积雪季起始时期9月初和结束时期次年4月末的数据质量存在一定问题，但是考虑9月初和4月末中国东北地区基本无积雪分布，因此不影响本数据集的使用。