青藏高原MODIS逐日无云积雪范围产品精度验证

李 诺，韩其飞，马 英，黄晓东

（1.南京信息工程大学 地理科学学院，江苏 南京 210041；2.兰州大学 草地农业生态系统国家重点实验室 草地农业科技学院，甘肃 兰州 730000）

0 引言

积雪作为冰冻圈的重要组成成分之一，对气候变化十分敏感。根据政府间气候变化专门委员会（Intergovernmental Panel on Climate Change，IPCC）于2021 年8 月发布的第六次评估报告，冰冻圈正呈现加速萎缩状态，北半球积雪范围也在加速减小［1］。积雪因为其特定的高反照率和低导热率，对地球能量平衡及水循环有着深刻的影响［2-3］。青藏高原地处中纬度地区，作为我国三大稳定积雪区之一，该地区的积雪影响东亚地区的大气环流及天气系统，其冰雪融水更是农业灌溉及周边众多河流的重要来源［4］。因此准确获取青藏高原地区的积雪信息对研究青藏高原的径流变化及气候演变有着重要意义。

遥感技术相较于传统的气象台站观测数据具有大尺度、获取信息方便快捷、不受地理条件限制等特征，成为大范围积雪监测的关键手段［5］。搭载在Terra 和Aqua 卫星的中分辨率成像光谱仪（Moderate Resolution Imaging Spectroradiometer，MODIS）因其在时间、空间和光谱分辨率的较大优势，在区域及半球尺度积雪监测中应用广泛。MODIS积雪产品第5 版本（V5）提供了二值积雪范围产品和积雪覆盖率产品［6］，两种产品的算法均基于归一化差值积雪指数（Normalized Difference Snow Index，NDSI）。二值积雪产品通过设置NDSI 阈值对每个像元进行分类，即一个MODIS 像元如果NDSI 值大于等于0.4 且波段2 反射率大于11%，则该像元被分类为积雪［7］，并引入归一化植被指数（Normalized Difference Vegetation Index，NDVI）用于提高MODIS 在林区的积雪分类精度［8］。评估结果表明，MODIS 逐日二值积雪产品（V5）在晴空状态下的总体分类精度介于85%～99%［9］。积雪覆盖率（Fractional Snow Cover，FSC）产品通过建立Landsat 获取的FSC 与NDSI 之间的线性关系生成［10］。2016年发布的MODIS 第6版本积雪产品（V6）不再提供积雪二值产品与积雪覆盖率产品，仅提供网格的NDSI数据。V6的产品结合高程重新定义了地表温度对高海拔地区积雪的分类阈值。另外，采用定量图像恢复算法重建了Aqua/MODIS 短波红外波段（波段6），替代了V5 版本使用波段7 计算NDSI的策略［11］。V6版本较V5版本的积雪产品，质量和精度均得到了有效的提升［12］。

云是影响MODIS 逐日积雪遥感产品准确获取地表积雪覆盖范围的主要因素之一。虽然MODIS逐日积雪产品在晴空状态下具有可靠的精度，但由于云的影响，导致该产品在区域积雪监测中受到极大的限制。因此针对MODIS 逐日积雪产品的云下信息恢复，一直是领域内学者关注的一个热点问题。国内外许多学者基于V5 产品针对云下信息恢复开展了大量的研究。唐志光等［13］利用三次样条函数插值法对MODIS 逐日积雪产品MOD10A1 进行去云处理，得到青藏高原MODIS 逐日无云的积雪覆盖率产品，该产品获取的积雪日数（Snow-Covered Days，SCD）与地面观测值得到的SCD 平均一致性为87.0%。侯小刚等［14］通过融合MOD10A1 与MYD10A1 数据，发展了一套适合新疆地区的逐日无云积雪范围数据集，去云后产品与气象台站积雪观测的一致性为88.1%，总体精度达到90.6%。黄晓东等［15］利用MODIS 上下午星获取的逐日积雪产品和被动微波数据AMSR-E 雪水当量产品，对MODIS 逐日积雪图像进行上下午星合成和邻近日像元分析，并结合数字高程数据利用雪线算法进行云下信息恢复，获取了青藏高原时空连续的MODIS 积雪逐日产品，总体分类精度达到90.7%。Gafurov等［16］在阿富汗东北部的Kokcha盆地，综合了双星合成、时间滤波、空间滤波等五种方法组成连续六步去云算法，生成的MODIS 逐日积雪产品，其总体分类精度也达到90%。Huang 等［17］提出一种基于隐马尔可夫随机场（Hidden Markov Random Field，HMRF）的时空滤波模型，充分利用了立方体内的时空信息对MODIS 逐日积雪产品进行去云处理，总体分类精度达88.0%，云量减少到1%。Hall等［18］提出一种基于MOD10C1积雪覆盖产品的时间滤波去云方法，生成MODIS Cloud-Gap-Filled（CGF）逐日云填补积雪产品。随着V6 版本积雪产品的发布，陆续有学者基于V6 积雪产品开展了类似研究。邱玉宝等［19］于2021 年11 月在科学数据银行（Science Data Bank，ScienceDB）更新了基于V6 版本的青藏高原MODIS 逐日无云积雪面积数据集，该数据集共使用八个连续步骤进行去云。Muhammad 等［20］于2021 年3 月在Data Publisher for Earth & Environmental Science（PANGAEA）发布了基于MODIS V6产品的亚洲高山区逐日无云二值积雪数据集。郝晓华［21］于2020 年11 月在国家冰川冻土沙漠科学数据中心发布中国逐日无云500 m 积雪面积产品数据集（2000—2020 年），该数据集基于MODIS V6 逐日地表反射率产品重建NDSI，通过训练不同土地覆盖类型下NDSI 识别阈值，然后利用隐马尔可夫随机场时空滤波模型进行了去云处理。上述产品自发布后，在青藏高原还未得到系统验证。

青藏高原作为我国三大积雪分布地区之一，海拔高且辐射强，导致积雪消融迅速、积雪破碎化严重，使得目前不同积雪产品在该地区的积雪监测精度均不够理想。因此，从应用角度出发，对MODIS逐日无云积雪产品在青藏高原进行系统评估，研究其在积雪监测中的精度及影响因素尤为关键。鉴于此，本文利用中高分辨率的Landsat影像数据作为参考值，对上述三套MODIS 逐日无云积雪产品进行验证，并分析土地覆盖类型和积雪覆盖率对MODIS积雪产品精度的影响，为评估青藏高原积雪范围及其时空动态变化提供参考依据。

1 研究区及数据

1.1 研究区概况

青藏高原位于我国西南部，平均海拔4 000 m左右，被称为“世界屋脊”“第三极”［22］（图1）。青藏高原西起帕米尔高原向东至横断山脉，自南喜马拉雅山脉南缘延伸至北昆仑山脉北缘和祁连山脉，总面积约2.57×106km2［23］。青藏高原以冰川、积雪、冻土等形态储存了巨大的水资源，我国的黄河、长江、怒江、澜沧江、雅鲁藏布江等均发源于青藏高原，有“中华水塔”和“亚洲水塔”之称［24］。青藏高原积雪水储量关系着所在区域及周边区域的生活和灌溉用水，影响高原及周边区域的植被生长，该地区的积雪具有重要的水文、生态及气候意义［25］。

图1 青藏高原土地覆盖类型及Landsat-8影像分布Fig.1 Land cover types in Tibetan Plateau and Landsat-8 images spatial distribution

1.2 MODIS逐日无云积雪产品

（1）中国逐日无云500 m 积雪面积产品数据集（MODIS CGF SCE）：来源于国家冰川冻土沙漠科学数据中心（http：//www.ncdc.ac.cn/）。该产品基于MODIS反射率产品MOD09GA/MYD09GA，利用高分辨率Landsat TM 数据作为真值，结合MODIS土地覆盖分类产品，确定林区和非林区下积雪判别的指标阈值，利用MODIS 积雪反演算法获取初级产品，通过隐马尔科夫随机场时空滤波算法进行云下信息恢复。该产品以HDF5 文件格式存储，每个HDF5文件包含18个数据要素，其中包括数据值、数据起始日期、经纬度等。本文采用无云积雪数据，为二值影像。其编码有0、1、2、3、4和255，对应的编码含义为陆地、影像识别积雪、去云插补积雪、雪深插补积雪、水体和填充值［21］。

（2）青藏高原MODIS 逐日无云积雪面积数据集（MODIS_Dysno_Cloudfree）：来源于ScienceDB（https：//www.scidb.cn/）。该产品以MODIS V6 积雪产品MOD10A1 和MYD10A1 为基础，通过上下午星合成、三天合成、“长时间”积雪和陆地法、邻近四像元法、高程滤波法、修改阴影区错误分类和最大积雪陆地范围掩膜这7 个连续步骤，获得MODIS少云积雪产品，而后将研究区根据坡度划分为“印度平原”“高原北部”“高原腹地”“藏东南山区”“帕米尔高原”“天山山脉”以及“喜马拉雅山脉”7 个区域，针对7 个区域采用拟合预期雪线方法去除全部云污染，获得MODIS 逐日无云积雪产品。本次所用产品为该数据集中青藏高原MODIS 逐日无云积雪产品面积数据集第二次更新版（2002—2021 年MODIS积雪产品C6.1版，http：//www.csdata.org/p/15/）。该产品为二值产品，影像分类代码与MOD10A1 一致，其编码有25、37、40、100、150 和200，对应含义为无积雪覆盖的陆地、湖泊、湖泊不确定、湖冰、湖冰不确定和积雪。该产品利用地面台站积雪深度观测数据进行验证，结果表明其积雪分类精度为78.4%，去掉雪深≤3 cm 的样本后，精度达到89.0%［19］。

（3）亚洲高山区MODIS 逐日无云积雪数据集（M*D10A1GL06）：数据来源于PANGAEA（https：//www.pangaea.de/）。该产品以MODIS 积雪产品MOD10A1 和MYD10A1 为基础，使用8d MOYDGL06*产品（双星8 天融合产品）分别代替Terra 和Aqua 雪产品中的云像素。然后与全球陆地冰空间测量计划（Global Land Ice Measurements from Space，GLIMS）发布的RGI（Randolph Glacier Inventory）6.0 冰川编目数据产品进行融合，最终达到99.99%的去云效果［20］。本产品为二值产品，其编码含义如表1所示。

表1 M*D10A1GL06积雪产品属性定义Table 1 Attribute definition of M*D10A1GL06 snow product

1.3 Landsat-8

前期研究大多基于站点雪深观测数据对MODIS积雪产品进行精度验证，但由于二者空间尺度的不匹配，导致基于点-面验证的结果存在较大的不确定性。因此，本研究采用空间分辨率为30 m的Landsat-8 陆地成像仪（Operational Land Imager，OLI）获取的影像作为验证数据。OLI 影像首先要选择晴空状态的数据，以避免云雪混淆造成验证数据集的不确定性；影像应分布在青藏高原积雪分布的主要区域，且能包含高原主要的土地覆盖类型；另外，为了避免破碎化积雪对OLI积雪制图的影响，验证时段选择积雪分布较为稳定的时期。因此，选取时间范围为2015 年12 月2 日至2016 年1 月23 日共34 景影像作为验证MODIS 逐日积雪产品的参考数据。

采用混合像元分析法首先获取OLI像元尺度的积雪面积比例（Fractional Snow Cover，FSC）。其中端元提取的规则如表2所示。将30 m的OLI积雪面积比例数据聚合成500 m，定义升尺度后的亚像元积雪面积比例阈值0.5，生成二值积雪数据。其中喜马拉雅山中东段的18 景Landsat 积雪数据来源于ScienceDB：“2013—2020 年喜马拉雅山中段和东段Landsat 8 积雪覆盖范围数据”。该数据集采用支持向量机（Support Vector Machine，SVM）分类方法，选取不同地形、阴影等条件的积雪特征训练样本进行积雪分类，结合地表水体等辅助数据及空间邻域分析进行分类后处理。通过对比Sentinel-2 高分辨率积雪分类数据，在900 m×900 m的网格内，其积雪覆盖率相关系数在0.95 以上，均方根误差约0.1%。该数据集为二值影像，其编码有1、2、100、200 和300，对应的编码含义分别为积雪、陆地、河流、冰湖和湖泊［27］。

表2 Landsat-8 OLI端元提取规则［26］Table 2 Endmember extraction rules for Landsat-8 OLI［26］

1.4 土地覆盖数据

MODIS 土地覆盖产品（MCD12Q1）用来比较不同土地覆盖类型下的MODIS 逐日无云积雪产品的精度。采用国际地圈生物圈计划（International Geosphere-Biosphere Program，IGBP）分类方案确定的17个土地覆盖类别，其中包括11个自然植被类别，3个开发和镶嵌土地类别，以及3 个非植被土地覆盖类别［28］。为了避免类型过多造成评价结果的复杂性，本文根据将IGBP 分类方案重新分类为七大类：森林、灌丛、草原、耕地、城市建设用地，冰川和裸地（图1）［29］。

2 精度验证方法

使用高分辨率遥感产品验证低分辨率产品的精度是遥感产品验证的常用手段［30-31］。本研究利用高分辨率的Landsat积雪图对MODIS两个版本的数据进行对比验证，采用精度评估指标包括总体分类精度（overall accuracy，OA）、低估误差（underestimated error，UE）、高估误差（overestimated error，OE），积雪分类精度P和Kappa系数。其中Kappa一致性检验是评价两幅图像之间一致性的方法［32］。当Kappa 系数为0～0.20 时，表明两幅图像的一致性极低；当Kappa 系数为0.20～0.40 时，表明两幅图像一致性一般；当Kappa 系数为0.40～0.60 时，表明两幅图像一致性中等；当Kappa 系数为0.60～0.80 时，表明两幅图像一致性极好；当Kappa 系数为0.80～1.00时，表明两幅图像高度一致［33］。

精度评价指标计算公式如下所示：

式中：a为Landsat 和MODIS 均记录为积雪的样本数；b为Landsat记录有雪而MODIS被分类为非雪的样本数，即漏判数；c为Landsat 和MODIS 均为非雪的样本数；d为Landsat记录为非雪而MODIS被分类为积雪的样本数，即为误判数；Pe表示偶然一致性，即在偶然机会下预测结果与真实结果一致的概率。

3 结果与讨论

3.1 精度评价

验证结果表明，MODIS CGF SCE、M*D10A1 GL06 和MODIS_Dysno_Cloudfree 产品的总体精度（OA）分别为88.2%、88.2%和89.6%，表明三种积雪产品有着相对较高的精度。其中MODIS CGF SCE低估误差（UE）最小，仅为9.4%，M*D10A1GL06和MODIS_Dysno_Cloudfree 积雪产品UE 分别为14.7%和18.6%。三套产品的高估误差（OE）分别为12.2%、8.9%和11.2%。由于MODIS CGF SCE 产品积雪分类精度（P）最高，达到90.6%，其他两种产品分别仅为85.3%和81.4%。Kappa 系数统计结果表明三种积雪数据与Landsat积雪数据的一致性均表现为极好，分别为0.642、0.629 和0.652，MODIS CGF SCE 和MODIS_Dysno_Cloudfree 与Landsat 积雪数据一致性最为接近，M*D10A1GL06产品略差。

M*D10A1GL06 和MODIS_Dysno_Cloudfree 产品均是基于MODIS V6 版本生产的去云积雪产品，两种产品的用户精度理论上是一致的，造成积雪分类误差不一致的原因主要是云重分类误差导致的，去云后M*D10A1GL06 的积雪分类精度优于MODIS_Dysno_Cloudfree，说明后者采用的云下信息恢复算法存在较大的误分类现象，与站点验证结论也有一定的差距［20］。MODIS CGF SCE 产品是基于大气校正的地表反射率数据生成的，并且分土地覆盖类型设置最优NDSI 阈值，在积雪分类时利用雪深等辅助数据进行插值，因此积雪分类精度最高；而MODIS V6产品是经过辐射校正的反射率数据计算NDSI。验证结果表明，基于大气校正的地表反射率数据在积雪信息提取方面更具优势，积雪分类精度优于另外两种去云产品。

图2 MODIS逐日无云积雪产品验证结果Fig.2 Validation results for three kinds of MODIS daily cloud-free snow cover products

3.2 土地覆盖对MODIS积雪产品精度的影响

表3中统计了MODIS逐日无云积雪产品在不同土地覆盖类型下的积雪分类精度。其中草原、裸地和森林的验证样本（N）分别占比67%（N=2 391 609），20%（N=729 769）和10%（N=370 443），其余验证样本分布在其他土地覆盖类型，因样本占比较小，未统计精度评价结果。结果显示，三种去云积雪产品的总体分类精度在不同土地覆盖类型下均大于80%，其中MODIS CGF SCE 在林区的总体分类精度最优，达到89.6%，MODIS_Dysno_Cloudfree在草原和裸地的总体分类精度最优，分别达到91.6%和84.0%。三种产品在裸地的总体分类精度均相对较差，介于81%～84%之间。MODIS CGF SCE 产品在草原的漏判和误判误差相当（～10%），在裸地存在严重的高估误差（20.9%），但在林区漏判误差较大（21.5%）。M*D10A1GL06 产品在裸地高估误差大，在林区漏判误差大，分别达到20.0%和27.2%。MODIS_Dysno_Cloudfree 产品仅在林区存在较为严重的漏判误差（21.3%）。MODIS CGF SCE 产品的积雪分类精度在草原和裸地类型下均优于另外两种产品，林区的积雪识别精度与MODIS_Dys-no_Cloudfree 产品相当，分别为78.2%和78.7%，M*D10A1GL06 产品在林区的积雪识别精度最低，为72.8%。与OLI 积雪数据相比，三种产品在草原和裸地类型下，一致性均较高，林区的一致性较差，但MODIS CGF SCE 产品在林区与OLI 一致性最高，为0.249。虽然该产品针对不同土地覆盖类型分别设置NDSI 阈值，并结合了林区积雪指数（Normalized Difference Forest Snow Index，NDFSI）用于提高林区积雪分类精度，但验证结果显示，相比基于MODIS 标准逐日积雪产品生成的去云产品，对林区积雪识别精度的改进并不明显，依然存在较高的漏判误差。

表3 土地覆盖类型对MODIS积雪产品精度的影响Table 3 Assessment of MODIS cloud-free snow products under different land cover types

3.3 FSC对MODIS积雪产品精度的影响

为了评估积雪覆盖率（FSC）对MODIS 逐日无云产品积雪分类精度的影响，本研究基于OLI 获取的FSC，对像元尺度的FSC 进行分级，验证不同FSC级别MODIS 逐日无云积雪产品的积雪分类精度。结果如图3 所示，三种积雪分类精度随着积雪覆盖率的增大在逐渐升高。其中MODIS CGF SCE 逐日无云积雪产品积雪分类精度在NDSI 大于10（原始值/100 代表实际值）均大于其他两种产品，Dysno_Cloudfree 产品积雪分类精度相对较差。进一步说明基于大气校正的地表反射率数据计算NDSI，进而对积雪信息进行提取，相比MODIS 标准逐日积雪产品具有一定的优势。

图3 不同积雪覆盖率下的MODIS积雪分类精度Fig.3 Accuracy of snow classification in graded fractional snow cover

4 结论

近年来，MODIS积雪产品在水文、水资源管理、气象和气候变化等众多领域得到了广泛的应用。在青藏高原地区，受限于该地区较强的太阳辐射，积雪消融迅速且破碎化严重，使得MODIS 积雪产品在该地区的准确性明显低于高纬度稳定积雪区。本文选取了青藏高原基于MODIS 数据发展的三套逐日无云积雪产品，其中MODIS CGF SCE 是基于MODIS 地表反射率数据生产的，其他两套产品（M*D10A1GL06 和MODIS_Dysno_Cloudfree）是基于美国雪冰数据中心（National Snow and Ice Data Center，NSIDC）发布的MODIS 逐日积雪产品进行去云处理后获取的。本文利用分辨率较高的Landsat-8 OLI 积雪影像数据，对比验证了三套逐日无云积雪产品的精度，并分析了土地覆盖类型及积雪覆盖率对产品精度的影响，得出以下结论：

（1）三套无云积雪产品的总体精度均优于80%，MODIS_Dysno_Cloudfree 积雪产品总体分类精度最高，达到89.6%，Kappa 系数为0.652。MODIS CGF SCE 低估误差最小，为9.4%，积雪分类精度也是最高的，达到90.6%。M*D10A1GL06 高估误差最低（8.9%）。

（2）土地覆盖类型对MODIS 积雪产品积雪分类的精度影响较大，尤其在林区。三种积雪产品在林区的积雪分类精度均较差，存在较高的漏判误差。虽然MODIS CGF SCE 产品针对林区进行了NDSI阈值优化，但效果甚微。

（3）随着积雪覆盖率的增加，三套无云积雪产品的积雪分类精度也随之增加，其中MODIS CGF SCE 产品在不同NDSI 等级（NDSI＞10）均优于其他两套产品。

综上，基于MODIS 地表反射率产品生产的逐日无云产品（MODIS CGF SCE），针对不同土地覆盖类型进行NDSI 阈值优化，在积雪识别精度方面提升较为明显，但是林区积雪识别误差较大的问题依旧没有得到有效解决。当积雪覆盖率较低时，积雪识别精度相比未经过大气校正的MODIS 标准积雪产品，精度相对较好，但提升不显著。因此，在积雪破碎化严重的青藏高原地区，除了考虑NDSI 阈值的优化，应进一步结合该地区的积雪特性，包括环境特征，进一步发展适合青藏高原的积雪制图算法。