MODIS,NDSI产品去云算法及最优阈值选择研究

王晓艳， 陈思勇， 郭慧， 谢佩瑶， 王建， 郝晓华

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

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

3. 中国科学院西北生态环境资源研究院, 兰州 730000;

4. 江苏省地理信息资源开发与利用协同中心, 南京 210023

积雪是全球气候变化的指示器，对气候系统具有显著的正反馈作用，是全球气候变化研究的关键变量（秦大河 等，2007；
Choi等，2010）。积雪在可见光和近红外波段的高反照率特征会影响地表能量平衡（Robinson 等，1986；
Huang 等，2017），积雪融化消耗能量会有效延缓温室效应（Robinson和Kukla，1985）。与此同时，积雪融化导致地表反照率下降，使得地面吸收更多的热量，促进温度升高（杨梅学 等，2000）。在干旱和半干旱地区，季节性融雪是一种主要的水源，有超过世界六分之一的人口依赖冰川和积雪融水维持生命（Barnett等，2005；
车涛和李新，2005；
王建和李硕，2005；
赵晓艳 等，2022）。积雪范围的定量化研究对气候变化、地表能量平衡和水文研究都具有重要意义（王建 等，2018）。

传统的气象站和积雪测线可以提供精确的地表积雪信息，但是由于站点的有限性并不能用于大范围积雪监测。卫星遥感技术可以提供大范围的积雪信息面上监测，弥补了常规观测在时间和空间上的不连续问题。中等分辨率成像仪MODIS（Moderate-resolution Imaging Spectroradiometer） 包含积雪识别所需的可见光和短波红外波段，并且具有很好的时间和空间分辨率，在积雪检测中得到广泛的应用（Hall等，2002；
Frei等，2012）。

大量研究表明，基于MODIS的积雪产品在晴空下具有较高的精度（Klein和Barnett，2003；
Huang等，2011；
Yang等，2015；
Simic等，2004；
Parajka和Blöschl，2006）。然而，云遮挡导致MODIS产品大量的数据缺失并限制了MODIS 产品的应用，如何高精度地恢复云下像元、生成时空连续的高精度MODIS 积雪产品成为近年来研究的热点问题。目前MODIS 积雪产品去云方法可以归纳为空间滤波方法、时间滤波方法、时空联合滤波方法和多源数据融合方法4 大类（Li等，2019a）。空间滤波是当前常用的去云方法，该方法通常在空间邻域上选择无云像元来估计云覆盖像元，对于大片的云遮挡区域不能获得满意的结果（Tran 等，2019；
Gafurov等，2015）；
采用Terra和Aqua一天合成以及多天合成的时间滤波方法，在连续多天的云覆盖情况下无法恢复云下数据，同时在积雪快速变化的阶段，多天合成会导致很大的误差（Wang等，2014；
Hall等，2010；
Dozier等，2008）；
时空联合的方法综合利用积雪的时空特征，在不同区域开展的实验表明，时空联合可有效提高积雪产品去云的精度（Li等，2017；
Da Ronco和De Michele，2014；
Parajka和Blöschl，2008；
邱玉宝 等，2017）；
多源融合方法充分利用光学观测、微波观测和地面观测等不同数据源的互补信息，也是一种有效的云去除方法。光学与微波数据相结合，可以用于较粗分辨率的积雪产品，但是产品的精度依赖于输入的遥感产品精度（Li等，2019b；
Huang等，2016；
Yu等，2016；
Deng等，2015；
Gao等，2010；
Liang等，2008；
Chen等，2018）。近年来，出现了一些基于新技术的去云算法。Dong 和Menzel（2016）在实地雪深观测的基础上，利用条件概率插值对MODIS 积雪图上的剩余云量进行重新分类，对德国西南部云覆盖和积雪变化较高地区的积雪产品进行了去云处理；
Huang等（2018）采用隐形马尔科夫随机场框架内集成MODIS光谱信息、时空上下文信息和环境关联的技术，对2006年—2008年积雪季上里约热内卢盆地的积雪产品去云；
Hou等（2019）采用时空滤波结合机器学习的技术，生成了2001 年—2016 年中国北疆地区逐日无云MODIS 积雪面积比例产品。以上的MODIS 去云算法综合MODIS 数据的时空特性并结合多源数据的优势特征，取得了较好的效果。

当前的积雪研究中常用的MODIS 积雪产品包括版本V005和V006，以上对积雪产品去云的研究主要针对V005 版本提供的二值积雪面积以及积雪面积比例产品。V005 版本二值积雪面积产品是基于归一化差值积雪指数NDSI（Normalized Difference Snow Index）的特定阈值生成的，积雪面积比例产品也是通过与NDSI 建立线性回归关系产生（Hall等，2002；
Hall和Riggs，2007；
Salomonson和Appel，2004）。然而，采用统一算法生成的MODIS积雪产品在不同地形、不同地表条件下的精度不尽相同，许多地区的积雪产品并不能达到令人满意的精度（黄晓东 等，2007；
曹云刚 等，2007；
唐志光 等，2013；
于小淇 等，2017；
王晓艳 等，2017；
Riggs等，2017）。NASA（National Aeronautics and Space Administration）于2016 年发布的MODIS 积雪产品V006 版本不再提供二值积雪面积和积雪面积比例产品，而是直接提供归一化差值积雪指数NDSI 产品。与传统的二值积雪产品和积雪面积比例产品相比，NDSI 产品可以保留更多的信息，从而更加灵活地应用于水文和生态过程模拟（Riggs 和Hall，2015）。然而，由于云的遮挡，当前的NDSI 数据集存在大量的数据缺失，影响了数据的使用。因此，亟需开展相关研究，生成高精度逐日无云NDSI数据集以更好地满足积雪遥感监测的需求。

Jing等（2019）提出了一种通过高斯核函数和误差校正的融合框架生成无云MODIS NDSI 产品的方法，在青藏高原地区取得了较高的精度。该方法重构的NDSI 值为0 和10—100 的整数，并未考虑0—10的NDSI低值。然而在森林地区，由于冠层的遮挡，森林积雪通常具有较低的NDSI 值（Rittger等，2013；
Hall 等，2007）。因而，恢复云下像元的NDSI原始值，包括低值的NDSI是非常必要的。

本文发展了一种基于邻近相似像元的云去除方法。在NDSI 影像中，对于t0时刻某一个云覆盖的目标像元，依次在前后日期无云的影像上寻找到与目标像元NDSI 值相近的像元作为邻近相似像元，然后用t0时刻的这些邻近相似像元进行加权计算得到目标像元的NDSI 值。本文以东北积雪区为实验区，选用一个积雪季（2017-10-01—2018-04-30）的NDSI 影像进行去云实验，并且将去云的NDSI 序列与气象站点雪深序列进行对比用于确定二值积雪制图中最优的NDSI阈值。

2.1 研究区概况

东北是中国3大主要积雪区之一，包括黑龙江、吉林、辽宁和内蒙古东部地区，这里冬季漫长寒冷，积雪丰富，研究区包括森林、草原、耕地、居民地等多种地表类型（图1）。东北平原被大兴安岭山脉、长白山脉和小兴安岭山脉所环绕，平均海拔为443 m。该区域的植被分布类型主要为温带针阔叶混交林、寒温带针叶林和温带草原3 类，植被覆盖度平均在70%以上，其中森林覆盖率达到40%（Che 等，2016）。选取该区域为研究区，获取逐日无云NDSI 时间序列，结合气象站点的雪深数据，可以揭示不同地表类型下NDSI 序列与积雪分布的变化规律，从而增强NDSI 在积雪识别中的应用。

图1 研究区及气象站点分布图Fig. 1 Map of study area and the location of the meteorological stations

2.2 数据及预处理

2.2.1 MODIS积雪产品

为了减小积雪监测的低估误差和高估误差，以便在全球尺度上获得精确的积雪分布，MODIS积雪产品V006（MOD10A1，MYD10A1）不再提供V005 中的二值积雪面积和积雪面积比例产品，而是采用保守的积雪监测方法，只提供MODIS NDSI和MODIS NDSI SNOW COVER 产品。其中，MODIS NDSI 的中像元的值在-10000—10000 之间，是将像元的原始NDSI 值扩大10000 倍的结果；
MODIS NDSI SNOW COVER 中保留了可能是积雪的像元的NDSI 值（0—100），并且已经识别出云（250）、内陆水体（237）、海洋（239）等地物类型，同时还标识出了缺数据（200）、无精度（202）、夜间（211） 和探测器饱和（254） 等（Riggs 和Hall，2015）。MODIS NDSI SNOW COVER 中，NDSI 值在0—0.1 时均被置为0，10—100 为原始NDSI 值扩大100倍的结果（Riggs等，2016）。

东北地区森林覆盖度高达40%，由于冠层的遮挡作用，森林积雪通常具有较低的NDSI 值（王晓艳 等，2017）。由于MODIS NDSI SNOW COVER将0—0.1 的NDSI 值均置为0，可能会丢失一些森林积雪的信息。因此我们选择MODIS NDSI 序列进行去云处理，以获取森林地区NDSI 的时序变化规律。本文选用了2017-10-01—2018-04-30 的MODIS Terra 和MODIS Aqua 的每日积雪产品（MOD10A1和 MYD10A1 V006），数据下载于美国国家雪冰数据 中 心（http：//nsidc.org/［2020-06-03］）。使 用MRT 工具将数据转换为WGS84 投影。采用MODIS NDSI SNOW COVER 数据层的云类别对无云的MODIS NDSI 进行云掩膜，再将去云结果与MODIS NDSI实际值进行对比来评价该去云算法的精度。

2.2.2 气象站点雪深数据

本文选取了研究区内均匀分布的24 个国家气象局气象站点，站点分布见图1。以站点的雪深数据作为积雪分布的真值，与同时期的NDSI 进行对比分析。站点数据采用以cm 为单位的整数记录了各个站点的积雪厚度，气象站点的雪深数据存在个别异常值，在使用前进行了剔除。表1为各站点的位置信息及对应的地表类型。

表1 气象站点位置及对应的地表类型Table1 Location and land cover type of the meteorological stations

2.2.3 MODIS土地覆盖产品

为了确定不同地表覆盖类型下NDSI 的最优阈值，本研究采用2017 年空间分辨率为500 m 的土地覆盖数据MCD12Q1。该产品是年度土地覆盖分类产品，有5 种分类标准，本文选用Land_cover_type5 分类标准，在使用前将土地覆盖类型合并林地、草地、耕地、水体和其他。该数据来源于美国航空航天局（https：//lpdaac.usgs.gov［2020-06-03］）。

为了充分利用无云的MODIS NDSI 数据信息，在使用基于邻近相似像元去云方法之前，首先使用MOD10A1/MYD10A1 结合和邻近时间数据滤波去除部分云像元。MOD10A1/MYD10A1 结合和邻近时间数据滤波的去云操作已经被广泛采用（Paudel和Andersen，2011；
Gao等，2010；
Hall 等，2010），因此接下来的部分将重点阐述基于邻近相似像元的去云方法。

3.1 基于邻近相似像元去云算法

地表的地物分布具有一定的相关性，在空间上并不是相互独立的，通常邻近像元具有很大的光谱相似性，被称为邻近相似像元，这些邻近相似像元具有两个基本特征：（1）近似相等的光谱特征；
（2）在多时相影像上具有相同的改变（Liu等，2019；
Zhu 等，2010；
Cheng 等，2017）。

如图2 所示，假设M为t0时刻云遮挡的目标像元，该像元在无云的t0+d时刻记为M′，假定SM′为M′的邻近相似像元，对应位置的SM为M的邻近相似像元。那么，

图2 邻近相似像元示意图Fig.2 Sketch map of similar pixels.

根据对邻近相似像元的假设，从t0到t0+d，目标像元和邻近相似像元的改变应该是一致的，也就是说Δ和Δ′相等。于是可以得到，

因此，选择已知的M′、SM和SM′，即可预测出云覆盖像元M 的NDSI 值。然而，由于一段时间内地表类型的改变，可能会使得某些邻近相似像元的变化并不相同。选取当前窗口下n个邻近相似像元，并对每个邻近相似像元赋予权重来计算目标像元的值，可以有效避免由于单个邻近相似像元的突变导致的较大预测误差。这里，与目标像元的NDSI 差值最小的前n个像元被认为是邻近相似像元。那么，目标像元的计算公式如下：

式中，xi，yi为第i个邻近相似像元的x，y坐标值，x0，y0为目标像元的x，y坐标值，L为当前窗口的大小。式（6）将Di约束在1到1+ 2 范围。

3.2 窗口大小与邻近相似像元个数的选择

为了确定式（4）中邻近相似像元的个数n和式（6）中窗口的大小L，本文首先将无云的NDSI影像进行云掩模，然后采用不同L和n值进行云去除，采用本文3.3 节的去云算法预测云下的NDSI值，然后计算预测的NDSI值与NDSI实际值的相关系数，结果如图3 所示。当L=15，n=20 时，相关系数r达到95.73%，继续增加L和n，r并无明显变化，综合考虑预测精度和计算耗时，本文选择的窗口大小为15×15，邻近相似像元个数为20。

图3 预测的NDSI与真值的相关系数变化曲面Fig. 3 The surface diagram of the correlation coefficient between predicted NDSI and true NDSI

3.3 去云流程

步 骤1，对MOD10A1 和MYD10A1 的NDSI 数据层进行合成得到MOYD。优先采用MOD10A1 的有效数据，即如果MOD10A1 为无云像元，则取该值作为MOYD 的值；
若MOD10A1 为云遮挡，MYD10A1 无云，则采用MYD10A1 为MOYD 的值；
若均为云遮挡，则该像元仍然为云像元，等待下一步去云操作。早期的MYD10A1 使用B7 代替B6进行积雪判别具有较低的精度（Riggs 等，2016）。Gladkova 等（2012）发展的QIR（Quantitative Image Restoration）算法被用来恢复Aqua MODIS 的B6 波段用于积雪监测，结果表明积雪监测的精度与Terra MODIS 的MOD10A1相当。目前，QIR 算法已集成到V006 的MYD10A1 中（Riggs 等，2016）。因此，联合V006 中MOD10A1 和MYD10A1 可以充分利用每天的无云信息。

步骤2，利用临近时间的数据进行前后时间段的去云处理。本文选用了前后两天的数据进行融合，结果记为ATC（Adjacent Temporal Composite）。如果选择的时间间隔太长，可能会由于积雪的变化导致去云结果存在较大的误差。如果仅采用前后一天的数据，则剩余的云量太多导致最后一步去云运算量过大。东北地区积雪相对稳定，采用前后两天数据融合是合理的。具体融合规则如下：对于当前的云遮挡像元，若前一天和后一天该位置像元均无云，则用这两天的像元平均值来替代云遮挡像元；
若只有一天为无云，则采用该无云像元的值替代云遮挡像元；
若前一天和后一天该位置像元均被云遮挡，则进一步考虑前两天和后两天的影像，云像元的替换规则与前后各一天的替换规则一致。

1.2 家庭农场在现代农业生产的优势 家庭农场的内在市场化运作方式和家庭经营的稳定性决定了在市场经济下有较高的竞争力，在发达国家，以家庭农场为经营方式的已存在近200年，今天依然为农业的最重要经营方式。我国家庭农场是在家庭承包经营基础上发展起来的，它保留了家庭承包经营的稳定性特点，同时又吸纳了现代农业市场化运作方式，其优势有以下几点。

步骤3，采用基于邻近相似像元的加权运算，对剩余的云遮挡像元进行去云处理，过程如下：

（1）对于t0时刻云遮挡的目标像元M，如果该目标像元在t0+d时刻是无云的，即M′已知，那么以该目标像元为中心，在15×15 窗口范围内选取前20 个与目标像元具有最小NDSI 差异的像元SM′作为它的邻近相似像元，同时还需满足SM也是无云的，即M′、SM′和SM均为无云像元。

（2）依据式（6）计算每个邻近相似像元到中心像元的距离Di，然后由式（5）计算出该邻近相似像元的权重Wi，最后根据式（4）计算出云覆盖的目标像元M的NDSI值。

（3）d 从1 开始，即先分别在前后一天的影像上寻找满足条件的邻近相似像元，如果前后一天都能找到满足条件的邻近相似像元，那么计算出的两个M平均值作为最终的预测值；
若只有一天满足，则用这一天的M作为最终值；
如果d=1时没找到符合条件的20 个邻近相似像元，则d=d+1 继续以上操作，直到找到满足条件的邻近相似像元，计算出M。

（4）移动窗口，重复（1）—（3），对下一个有云像元进行去云处理，直至整个影像的云覆盖率为0。

3.4 去云结果与精度验证

本研究对东北地区2017-10-01—2018-04-30的MODIS NDSI 数据进行了去云处理，得到了无云的NDSI 序列。图4 给出了部分影像逐步去云的结果，图中CF为云覆盖比例。

图4 最后一列从上到下分别为2018-01-01、2018-02-01 和2018-03-01 云去除后的NDSI 影像，影像上NDSI 的分布符合东北积雪分布规律。在这些日期，大兴安岭、小兴安岭以及长白山地区都被积雪覆盖，而森林地区积雪NDSI 较低，在0—0.3 之间；
高纬度地区的草地耕地居民地积雪具有较高的NDSI 值，在0.5 以上；
西南部草原耕地少雪，NDSI为负值。

图4 MODIS NDSI去云影像结果Fig. 4 Cloud removal results for the MODIS NDSI images

为了进一步定量验证去云的精度，本文采用“云假设”的方法。“云假设”的思想是在无云的影像上添加云掩膜，然后对添加云掩膜的影像进行去云处理，再将去云结果与真实的影像进行对比来评价算法的好坏。为了使得添加的云接近真实的云分布，本文选取了每个月中云量最小的两景影像作为NDSI 影像的真值，将本月中云量最大的云覆盖到该NDSI 影像上，这里采用的是MODIS NDSI SNOW COVER中的云类别（值为250的像元）对MODIS NDSI进行云掩膜。

MOD10A1和MYD10A1的合成以及ATC 已被广泛采用并证明具有较高的精度（Paudel和Andersen，2011；
Dong和Menzel，2016）。这里只讨论基于邻近相似像元的去云结果，表2中的样本比例即为采用基于邻近相似像元进行去云的像元占总像元的比例，时间窗口是指完全去除云遮挡需要的天数。使用相关系数r、均方根误差RMSE和平均偏差MAE这几个指标（式（7）—（9））对去云结果进行精度评价。

表2 基于“云假设”的去云结果评价Table 2 Cloud removal result assessment based on “cloud assumption”

表2 给出了10 景实验影像在基于邻近相似像元去云过程中云像元的比例、时间窗口的大小以及r、RMSE、和MAE 的值。较高的r值和较低的RMSE 以及较低的MAE 代表反演的结果更加接近真值。10景实验影像去云结果中，相关系数r的平均值为0.95，均方根误差RMSE 的均值为0.08，平均绝对误差MAE 的均值为0.05。时间窗口从22 天到61 天不等，时间窗口的大小与云覆盖比例并无显著相关，对去云结果的精度也无明显影响。

4.1 NDSI时序变化

图5 为一些典型气象站点所在像元的去云NDSI序列与站点雪深序列对应图（站点编号见图1和表1）。站点所在的像元包括耕地、居民地、林地和混合像元，序列的日期从2017 年10 月1 日—2018 年4 月30 日共212 天，图中横轴day1 代表2017年10月1日。

图5 典型地表的NDSI序列与雪深变化曲线（横轴的序列日期从2017-10-01—2018-04-30；

第一天代表2017-10-01）Fig.5 The curve of NDSI and snow depth at typical surface （The date on the x-axis is from October 1， 2017 to April 30， 2018， and the first day is October 1， 2017）

图5（a）为气象站点1 所在像元的NDSI 及雪深序列，该像元为纯净的居民地像元。稳定积雪时期该像元的NDSI 值在0.6 附近，雪深在5 cm 以上。2017 年11 月10 日气象站点观测到雪深由0 变为5 cm，对应的NDSI 由负值上升为正值，之后一直在0.6 上下浮动，直到2018 年3 月16 日融雪开始。融雪过程自2018 年3 月16 日持续到2018 年3 月28 日，期间气象站点所测雪深由15 cm 逐渐减小，NDSI 值也逐渐降低。当雪深降为0 时，即积雪完全融化，此时NDSI 变为负值，融雪过程持续了12 d。

图5（b）对应的气象站点2所在像元为纯耕地像元，该像元的NDSI 时序变化与气象站点所测的雪深序列也有很好的对应关系。气象站点观测到雪深由0 变为15 cm 的同时，NDSI 由负值骤升到0.6。之后漫长的冬季，NDSI 一直保持在0.6 左右，站点观测的雪深在10—20 cm。2018 年3 月18 日开始，气象站点所测雪深逐渐降低，NDSI 值也随之降低。至2018 年3 月26 日积雪完全融化，雪深降为0，NDSI 同时降为负值，积雪融化过程大约持续8天。

图5（c）对应的气象站点4所在像元为耕地居民地的混合像元，该像元的NDSI 变化与雪深变化也有很好的一致性。在稳定的积雪季，像元NDSI值在0.6上下波动，雪深在5—10 cm。从2018年3月26日开始至2018年3月29日，雪深由5 cm 降为0，NDSI值也由0.6降为负值。因雪深较浅，融雪过程仅持续3天。

图5（d）对应气象站点19，该站点所在的像元为大兴安岭的落叶林。林地NDSI 的变化与站点所测雪深具有很好的对应关系。2017年10月初有一次5 cm 降雪，因为积雪快速融化，NDSI 值有响应但并未能升至大于0。接下来10月下旬和11月初的降雪NDSI值都有相应的变化响应。2017年11月10日至2018年3月18日为稳定的积雪期，雪深在15 cm左右，此时NDSI 的值在0.2—0.4 之间，低于积雪期耕地和居民地的NDSI值。从2018年3月9日开始NDSI值持续下降，至2018年3月25日NDSI降为负值，雪深由10 cm降为0，该融雪过程持续16天。

图5（e）对应的气象站点22 位于小兴安岭的针阔叶混交林。自2017年11月10日站点观测到的雪深由0变为5 cm，接下来的积雪季，雪深逐步增加，2018年3月份雪深增至25 cm，而NDSI的值始终在0—0.2之间。融雪开始于2018年3月28日，至2018 年3 月31 日雪深变为0，同时NDSI 降至小于0。因该像元融雪开始时间较晚，较高的气温使得融雪过程仅持续了3 天时间。

图5（f）为纬度较低的气象站点15所在像元，从该站点的积雪观测数据来看，在2018年3月12日有一场明显的降雪，NDSI 的时序变化图也准确地反映了此次降雪过程。

图5（a）—5（e）站点均位于高纬度地区，冬季漫长寒冷，有明显的积雪稳定期。NDSI 与气象站点所测雪深有很好的对应关系。在气象站点观测的雪深由0 变为非零时，NDSI 会有相应的跃升。随着积雪融化，观测的雪深逐渐减小为0，相应的NDSI 值也会变为负值。图5（f）表明，NDSI 也能准确反映短暂的降雪过程。

4.2 积雪提取中NDSI的最优阈值选择

图6 表明，NDSI 时序变化与观测的雪深具有很好的对应关系，因此NDSI 是进行积雪识别的有效指数。接下来讨论NDSI 的阈值与积雪识别精度的关系。

图6 积雪识别精度随NDSI阈值的变化曲线Fig. 6 The curve of snow recognition accuracy with NDSI threshold

通常，站点雪深数据大于等于1 cm 时，认为该站点所在的像元为有雪像元；
站点雪深为0，则认为该站点所在的像元为无雪像元（Parajka和Blöschl，2006；
Ault 等，2006）。那么根据站点雪深数据，可以将站点所对应的像元分为积雪像元和无雪像元，作为积雪分布的真值。选取一定的NDSI 阈值，同样可以将像元分为积雪像元（NDSI>给定阈值）和无雪像元（NDSI<给定阈值），得到二值积雪分布。以雪深的积雪分类作为真值，对采用不同NDSI阈值制作的二值积雪分布图进行精度评价，评价指标包括总体精度OA （Overall Accuracy）、低 估误差UE （Underestimation Error）和OE 高估误 差（Overestimation Error）（Parajka 和Blöschl，2008），见式（10）—（12），其中，a、b、c、d含义如表3 所示。因为森林和非森林地区的NDSI值有较大差异，这里分别对森林和非森林地区进行评价。

表3 基于雪深和NDSI的二值积雪识别的混淆矩阵Table 3 Confusion matrix comparing the snow depth（SD） and NDSI binary snow cover

（1）非森林像元。采用了图1中的非森林站点1—17，雪深大于等于1cm 则记为snow。随着NDSI阈值的变化，总体精度的变化如图6（a）：当NDSI=0.1时总体精度最高，达到95.6%，此时，高估误差和低估误差分别为1.3%和4.7%。

（2）森林的像元。通常站点安置在森林附近的道路旁或者居民区，站点所在的像元为混合像元。为了统计纯森林像元的NDSI 变化与雪深的关系，参照Google Earth 选取站点附近（<2 km）的纯森林像元，用站点的雪深数据代表周围森林像元的雪深。本文选取了站点18—24 周围的18 个纯森林像元，对比森林像元的NDSI 序列与站点雪深的变化规律，统计不同NDSI 阈值下积雪识别的精度。在森林地区，随着NDSI 阈值的变化，森林像元的积雪识别的精度变化如图6（b）：当NDSI=0时，精度最高可达到93.5%，低估误差和高估误差均为3.8%。选取的NDSI 阈值大于0.1，低估误差会急剧增加，总体精度也相应快速降低。

本文提出了一种基于邻近相似像元的去云方法，对MODIS 积雪产品NDSI V006 进行去云处理，得到逐日无云NDSI 序列。选取中国东北地区2017-10-01—2018-04-30 一 个 积 雪 季 的MODIS NDSI 产品进行去云实验并采用“云假设”方法进行精度评价。结果表明，本文提出的基于邻近相似像元的去云方法可以很好的恢复云下像元的NDSI 值，预测值与真值的相关系数达到0.95，是一种有效的去云方法。

无论在林地还是非林地、高纬度还是低纬度地区，无云的NDSI 序列与气象站点观测到的雪深序列都具有很好的一致性，NDSI 可以准确地反映站点所在像元的积雪的变化。无雪时的林地和非林地像元的NDSI 基本为负值，在有雪的情况下，林地像元和非林地像元的NDSI 存在差异。非林地积雪像元具有较高的NDSI值，约为0.6。由于森林冠层的影响，林地积雪的NDSI 值较低，大兴安岭落叶林具有较高的透过率，林地积雪的NDSI 值在0.3 左右，而小兴安岭的针阔叶混交林透过率较低，林地积雪的NDSI值在0—0.2之间。

将雪深大于等于1 cm 记为有雪，作为地表积雪分布的真值，统计不同NDSI 阈值下积雪识别的精度。结果表明，在森林和非森林地区的NDSI 的阈值分别取0 和0.1 时，可以得到最高的积雪识别精度95.6%和93.5%。Zhang 等（2019）基于中国境内的200多个站点的日雪深数据，推荐中国境内对MODIS V6 产品使用NDSI 阈值0.1 来判定积雪。由于采用的森林站点数量有限，Zhang 等得出的NDSI=0.1 的阈值主要是指在非森林地区，本研究的结论与此一致。

Hall 等（2002）给出的MODIS 全球积雪面积产品NDSI阈值为0.4；
郝晓华等（2008）选取祁连山中部山区的3 个积雪子区进行积雪识别并采用Landsat 数据进行精度评价，结果表明，NDSI 平均阈值为0.33 时积雪识别的总精度达到最高。本文给出的NDSI 阈值与以上研究的结果有较大差异，主要是因为去云之后的NDSI 产品，在积雪提取中不需要再考虑云雪混淆问题，只需区分有雪和无雪的情形，故NDSI 阈值适当降低，减少低估误差的同时不会引入将云错分为雪的高估误差。

图6 表明，非森林像元的NDSI 阈值从-0.1 到0.4，积雪识别的精度均在90%以上，NDSI=0 时的积雪识别精度为93.2%，略低于NDSI=0.1 时的95.6%。而森林像元的NDSI 阈值大于0.1 之后，低估误差急剧增加，积雪识别的总精度明显下降。因此，如果在积雪识别过程中存在地表分类图，则可以在森林和非森林地区分布采用NDSI 阈值为0 和0.1，若没有相应的地表分类图，则可以统一使用NDSI=0作为积雪识别的阈值。

本文提出的MODIS NDSI 产品去云方法，可以高精度地恢复云下像元，获取逐日无云的NDSI 产品，使其更好地应用于积雪识别中。基于站点雪深选取的东北地区森林和非森林地区NDSI 最优阈值，为NDSI 在积雪识别中的应用提供了有力的支撑。

猜你喜欢 雪深气象站积雪 珠峰上架起世界最高气象站环球时报(2022-05-05)2022-05-05我们福建文学(2019年12期)2019-08-06心灵气象站趣味(语文)(2019年3期)2019-06-12小二沟地区1971～2010年积雪的特征分析现代农业(2018年5期)2018-05-23大粮积雪 谁解老将廉颇心炎黄地理(2017年10期)2018-01-31积雪少年文艺·开心阅读作文(2017年1期)2017-02-24铁路防灾雪深图像采集的设计和实现铁道通信信号(2016年5期)2016-06-01铁路雪深监测系统数据过滤算法研究铁路计算机应用(2016年10期)2016-02-16青藏高原积雪深度时空分布与地形的关系自然资源遥感(2015年4期)2015-12-25自动气象站异常记录分析及处理方法Advances in Meteorological Science and Technology(2015年5期)2015-12-10