青藏高原不同土地覆盖类型下积雪面积判别算法优化

doi:10.7522/j.issn.1000-0240.2022.0364

青藏高原不同土地覆盖类型下积雪面积判别算法优化

谢佩瑶^,, 韩超, 欧阳志棋, 王晓艳^,

兰州大学资源环境学院，甘肃兰州 730000

Optimization of snow area discrimination algorithm under different land cover types in Qinghai-Tibet Plateau

XIE Peiyao^,, HAN Chao, OUYANG Zhiqi, WANG Xiaoyan^,

College of Earth and Environmental Sciences，Lanzhou University，Lanzhou 730000，China

通讯作者: 王晓艳，副教授，主要从事积雪遥感研究. E-mail： wangxiaoy@lzu.edu.cn

收稿日期: 2022-01-05 修回日期: 2022-05-25

基金资助:

国家自然科学基金面上项目“森林冠层降雪截留遥感监测方法研究”. 42271373

Received: 2022-01-05 Revised: 2022-05-25

作者简介 About authors

谢佩瑶，硕士研究生，主要从事积雪遥感研究.E-mail：xiepy19@lzu.edu.cn , E-mail：xiepy19@lzu.edu.cn

摘要

MODIS V006版本数据仅提供了归一化积雪指数（NDSI），而用户往往关心的是直观的积雪分类，包括积雪范围或积雪覆盖率。美国国家冰雪数据中心推荐全球积雪范围最佳的NDSI阈值为0.4，但是青藏高原地形复杂多样，积雪斑块化特征明显，单一阈值并不能精确地判识不同下垫面上的积雪。青藏高原被称为地球的第三极，是中国三大稳定积雪区之一，蕴藏了大量的淡水资源。随着全球气候变暖，青藏高原地区积雪融化时间提前，冰川融水增加，影响河流水量，造成洪涝灾害，进而影响人类正常生产生活，因此通过确定不同下垫面阈值，改善传统阈值的积雪高估低估现象，提高积雪识别精度，进而更准确地探究青藏高原积雪状况，显得尤为迫切。本文以青藏高原为研究对象，首先生成MODIS逐日无云NDSI序列并进行验证；其次对应站点雪深数据与NDSI序列，证实在下垫面为林地和非林地的区域，去云NDSI序列与站点雪深均有良好的对应关系，确定不同下垫面最优阈值范围；最后在最优阈值范围内通过混淆矩阵确定最优阈值。计算得出，林地NDSI=0.03时，总体精度最高为94.02%，在该NDSI之下，高估误差OE和低估误差UE分别为1.21%和4.60%；非林地NDSI=0.26时，总体精度OA最高为94.27%，在该NDSI之下，高估误差OE和低估误差UE分别为0.51%和5.03%。因此选取优化后林地阈值为NDSI=0.03，非林地阈值为NDSI=0.26。为避免地面常规观测资料尺度上的局限性，本文采用高精度的Landsat 8 OLI卫星数据识别结果，作为“真值”对优化后阈值的判别结果进行“像元—像元”级别的验证。在定量验证中，优化后 NDSI阈值对MOD10A1 V006积雪判别结果的总体精度OA为84.21%，高估误差OE为5.33%，低估误差UE为10.46%；传统阈值对MOD10A1 V006积雪判别结果的总体精度OA为82.86%，高估误差OE为1.48%，低估误差UE为15.66%。可以看出在定量验证中，优化后阈值的积雪判别精度更高。同时在定性验证中，积雪大面积集中的区域，新的阈值与传统阈值提取效果均相对较好；积雪相对分散破碎的区域，优化后阈值能提取出大量积雪，传统阈值则不能。这表明考虑不同土地覆盖类型下的NDSI阈值优化可以有效地提高青藏高原积雪判别精度，为NDSI在积雪识别中的应用提供有力的支撑，有助于更准确地了解该地区积雪分布状况。

关键词： STAGFM ; 青藏高原 ; 林地 ; 积雪提取

Abstract

The data provided by MODIS V006 version is the Normalized snow cover Index （NDSI）， but the majority of users are often concerned with the intuitive snow cover classification results， including snow cover extent or snow cover rate. The National Snow and Ice Data Center （NSIDC） recommended 0.4 is the best NDSI threshold for the global snow cover. However， the Qinghai-Tibet Plateau has complex and diverse terrain and obvious snow patch characteristics， so a single NDSI threshold of 0.4 cannot accurately distinguish the snow cover on different underlying surfaces. The Qinghai-Tibet Plateau， known as the third pole of the earth， is one of the three stable snow areas in China and contains a large amount of fresh water resources. With global warming， the Tibetan plateau ahead of time， the snow is melting glaciers， increase， affect the rivers of water， causing floods， and thus affect the normal production of human life， so determining different underlying surface threshold， improve the traditional threshold value of snow overestimated underestimate phenomenon， improve the identification accuracy of snow， and then more accurate study of the Tibetan plateau snow conditions， Is particularly urgent. In this study， This paper takes the Qinghai-Tibet Plateau as the research object. Firstly， MODIS daily cloud-free NDSI sequence is generated and its reliability is verified. Secondly， the underlying surface is forested and non-forested areas， and the NDSI sequence of cloud removal has a good corresponding relationship with the snow depth at the site. NDSI can accurately reflect the snow melting phenomenon of the pixel where the station is located. Determine the optimal threshold range of different underlying surface； Finally， the optimal threshold was determined by confusion matrix within the optimal threshold range. When NDSI=0.03， the highest overall accuracy was 94.02%. Under this NDSI， the overestimation error OE and underestimation error UE were 1.21% and 4.6%， respectively. When NDSI=0.26 for non-forestland， overall accuracy （OA） is 94.27%. Under this NDSI， he overestimates error （OE） and underestimate error （UE） are 0.51% and 5.03%， respectively. Therefore， the optimized threshold of forestland is NDSI=0.03， and that of non-forestland is NDSI=0.26. Because snow cover is a large scale phenomenon， the conventional observation data are mostly point scale observations. In order to avoid the limitations in the scale of conventional ground observation data， this paper uses the high-precision Landsat 8 OLI satellite data identification results as the “truth value” and the snow discrimination results of the optimized threshold and the snow discrimination results of the traditional threshold to verify the “pixel to pixel” level. In quantitative verification， the overall accuracy （OA） of the optimized NDSI threshold to MOD10A1 V006 snow discrimination results is 84.21%， the overestimation error OE is 5.33%， and the underestimation error （UE） is 10.46%. The overall accuracy （OA） of the traditional threshold for the snow discrimination results of MOD10A1 V006 is 82.86%， the overestimation error （OE） is 1.48%， and the underestimation error （UE） is 15.66%. It can be seen that in the quantitative verification， the snow discrimination accuracy of the optimized threshold is higher. At the same time， it can be seen from the qualitative verification that the new threshold and the traditional threshold are relatively good at snow recognition in the area with large area of concentrated snow. In the region with relatively scattered and broken snow， the optimized threshold can identify a large number of snow， while the traditional threshold cannot identify the same number of threshold. These results indicate that NDSI threshold optimization considering different land cover types can effectively improve the accuracy of snow discrimination on the Qinghai-Tibet Plateau， which provides a strong support for the application of NDSI in snow recognition. It is helpful to understand the snow distribution in this area more accurately.

Keywords： STAGFM ; Qinghai-Tibet Plateau ; forest ; snow extraction

PDF (3815KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

谢佩瑶, 韩超, 欧阳志棋, 王晓艳. 青藏高原不同土地覆盖类型下积雪面积判别算法优化[J]. 冰川冻土, 2023, 45(3): 1168-1179 doi:10.7522/j.issn.1000-0240.2022.0364

XIE Peiyao, HAN Chao, OUYANG Zhiqi, WANG Xiaoyan. Optimization of snow area discrimination algorithm under different land cover types in Qinghai-Tibet Plateau[J]. Journal of Glaciology and Geocryology, 2023, 45(3): 1168-1179 doi:10.7522/j.issn.1000-0240.2022.0364

0 引言

冰冻圈泛指地球表层中水体以固态保存的圈层，包括冰川、积雪、多年冻土及河流湖泊中的淡水冰、海冰、地下冰等^［1］。积雪作为陆地表面覆盖物中的重要组成部分，不仅在冰冻圈和气候系统扮演着不可或缺的角色，同时也是表征局部地区物候最活跃的要素之一，在气象、水文和生态过程中发挥着重要作用^［2-3］。青藏高原是北半球中低纬度海拔最高、积雪覆盖最大的地区^［4］。积雪是青藏高原季节性变化最大的下垫面和重要的生态环境组成要素，冰雪融水是高原及其下游地区重要的水资源^［5-6］。同时，高原积雪作为一种重要的陆面强迫因子，与东亚、南亚季风以及长江中下游的旱涝等灾害性天气紧密相关，是短期气候预测的重要指示因子和全球气候变化最为敏感的响应因子之一^［7］。青藏高原地区平均海拔超过4 000 m，环境恶劣，交通条件极为不便，仅靠实测手段和摄影测量难以实现高效快速的大范围积雪监测^［8-9］。遥感影像具有大面积同步观测，实时性等特点，对于青藏高原地区积雪提取具有巨大优势^［10］。国内外利用不同传感器（Landsat、SPOT^［11］、SMMR^［12］、SSM/I^［13］、AVHRR^［14］、MODIS^［15］、AMSR-E^［16］）生产了多种积雪产品，MODIS是地球观测系统（Earth Observation System，EOS）系统中主要的传感器之一^［17-18］，其积雪范围产品具有较高的时空分辨率，是目前最广泛使用的积雪覆盖产品。MODIS数据使用SNOMAP算法制备积雪范围产品，其中归一化积雪指数（normalized difference snow index， NDSI）阈值是最主要的积雪判识条件^［19-20］。

NDSI阈值法是最早针对积雪范围进行智能化制图的探索，通过设置NDSI大于0.4制备二值积雪产品^［21］。但NDSI大于0.4在森林地区的积雪识别精度明显低于非森林地区，Klein等^［22］利用不同指数差异提取雪像元。通过对不同森林类型有雪和无雪时NDSI值及NDVI值的变化统计，综合使用NDSI和NDVI阈值法来提取森林积雪。将SNOMAP算法对森林积雪的识别规则定义为NDVI≥0.1时，将0.1≤NDSI≤0.4的像元识别为积雪，通过降低NDSI阈值，可以有效减小在森林区积雪像元的低估，但NDSI的降低同时也会将非积雪像元高估为积雪，增加高估误差。

青藏高原土地覆盖包含草地、稀疏植被、耕地、针叶林、阔叶林等多个类型^［23］。Wang等^［24］基于Landsat数据，使用近红外波段代替NDSI计算公式中的可见光波段，构建归一化差值林地积雪指数NDFSI（normalized difference forest snow index）。当NDSI小于0.4时，以NDFSI阈值0.4进行森林积雪的识别，对黑河上游青海云杉林地积雪的识别精度可以达到93.92%。Zhou等^［25］验证了中国地区的MODIS V6最优NDSI阈值为0.1，郝晓华等^［26］验证了祁连山地区的MODIS数据最优NDSI为0.33，王雪璐等^［27］验证了青海地区的MODIS最优NDSI为0.37，王建^［28］在北疆对AVHRR数据进行积雪提取的最优NDSI阈值范围为0.2~0.5，在黑河对MODIS数据的最优NDSI阈值范围为0.35~0.45，在天山站对Landsat TM数据的最优NDSI阈值范围为0.57~0.72。高扬等^［10］在青藏高原针对草地，稀疏植被和其他下垫面地表类型计算了MODIS积雪范围产品的最优NDSI阈值，分别为0.33、0.40、0.47。王玮^［29］采用最大似然法分类获取TM数据的积雪分类结果，并将其作为“真实”的雪盖，将青藏高原牧区积雪识别的MODIS的NDSI阈值确定为0.35。由于青藏高原地区地形复杂多样，积雪覆盖范围在空间上分布不连续，受混合像元、地形等多种因素的影响使其精度较差，其精度明显低于其他地区^［30］。使用单一NDSI阈值进行积雪识别判别，未完全考虑积雪的物理特性、土地覆盖类型、大气状况等自然条件的影响，必然会降低积雪识别精度^［31］。同时青藏高原东南部地区有茂密林地，由于其在青藏高原地区总体占比小于百分之10%，在各类研究中均没有对该地区进行深入讨论。因此根据区域实际情况设置符合当地情况的NDSI阈值是十分必要和迫切的。

本文以MODIS V006版本逐日积雪范围产品MOD10A1/MYD10A1为基础，生成逐日无云NDSI序列，以气象站点雪深数据作为参考真值，获取不同土地覆盖类型（林地与非林地）的NDSI最优阈值，生成MODIS积雪二值产品。使用Landsat 8 OLI数据的验证结果表明，优化后的阈值有效提高了MODIS积雪产品的精度。

1 研究区与数据

1.1　研究区概况

青藏高原是中国最大、世界海拔最高的高原，西起帕米尔高原，东至横断山脉，南自喜马拉雅山脉南缘，北至昆仑山—祁连山北侧，范围为26°00′12″~39°46′50″ N，73°18′52″~104°46′59″ E，被称为“世界屋脊”和“第三极”，同时该地区也是我国的三大稳定积雪区之一^［32］。青藏高原的冰川积雪大都处于高海拔低纬度地区^［33］。积雪资源丰富，是长江、黄河、雅鲁藏布江等诸多河流的发源地。春季积雪融水是河流的主要补给来源，直接影响江河流量，甚至引发春汛。图1中标注为红色的点为典型气象站点。用于确定阈值范围，其中A~H共8个点为非林地区域典型气象站点，I、J点为林地区域典型气象站点，标注为黑色的点为典型气象站点用于计算总体精度OA，高估误差OE和低估误差UE，确定最优阈值。其中1~15共15个点为非林地气象站点，均为2017年数据，16~24共9个点为林地气象站点，点16为2013年气象站点数据，点17~19为2014年气象站点数据，点20~22为2015年气象站点数据，点23~24为2016年气象站点数据。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 研究区、气象站点和Landsat 8检验数据以及土地覆盖类型的空间分布

Fig. 1 Spatial distribution of land cover types in the study area， meteorological stations， and Landsat 8 test data

1.2　数据来源

为了便于积雪特征的研究，资料选取前一年9月1日至当年8月31日，作为一个水文年数据，同时选用前一年11月1日至当年3月31日作为一个积雪季^［34］。本文采用MODIS积雪产品数据利用STAGFM方法（时空自适应去云算法）进行逐日无云积雪数据的合成，利用气象站点数据获取不同类型下垫面的阈值^［35］。本文采用具有较高可信度的Landsat 8 OLI卫星数据判识结果来做“像元—像元”级的验证，相对于地面常规观测资料的局限，能够更好地验证MODIS卫星数据积雪判识结果的精度。

1.2.1　MODIS积雪产品数据

本文采用MODIS积雪产品V006版本^［17-18］，包括MOD10A1数据与MYD10A1数据，V006版本包含七层数据：2个NDSI数据、2个质量评估数据、1个反照率数据和2个轨道信息（表1）。该版本属性定义如下表所示。本研究采用NDSI积雪范围（NDSI Snow Cover）和NDSI（归一化积雪指数）进行积雪识别^［36］。

表1 MODIS V006版NDSI Snow Cover属性定义

Table 1 MODIS V006 NDSI Snow Cover attribute definition

数值	属性
0~100	NDSI值
200	缺失数据
201	未确定
211	夜晚
225	陆地
237	内陆水域
239	海洋
250	云
254	探测器饱和
255	填充

新窗口打开| 下载CSV

MODIS V006积雪产品（MOD10A1/MYD10A1）与MODIS V005版本相比进行了很大的改进，改进了大气校准方法，更好地区分云和积雪，但是，由于积雪和云具有相似的光谱反射特性，使得当前的MODIS积雪产品在空间和时间上表现出质量不高且不完整的现象。

1.2.2　气象站点雪深数据

收集了青藏高原2017年11月1日到2018年3月31日共计151天的240个地面气象台站地面雪深数据集，以及青藏高原2013—2017年各年11月1日到下一年3月31日共计151天的28个林地地面气象台站地面雪深数据集，包括台站号、经纬度、海拔及雪深等信息。数据来源于中国国家气象科学数据共享服务平台（http：//data.cma.cn/）。参考于小淇等^［37］的方法将台站雪深“>0.5 cm”的数据表示为有雪；“=0 cm”表示无雪，即为陆地；32700表示微量积雪，将其赋为陆地，即无雪；32766表示缺测数据（图1）。

1.2.3　Landsat 8 OLI积雪产品数据

Landsat 8数据主要包含OLI的9个波段和TIRS的2个波段，共11个波段，目前，SNOMAP算法是用于积雪遥感监测最普遍的技术手段，其核心是利用固定阈值的归一化差分积雪指数（normalized difference snow index， NDSI）进行积雪识别。因此在下垫面属于非林地的地区使用SNOMAP方法来提取积雪范围^［38］。在下垫面为林地的区域使用多指标的积雪监测算法分类提取积雪范围，进行“真值”计算。对Landsat 8 OLI影像进行升尺度计算，使Landsat 8 OLI影像与MODIS影像空间分辨率一致，便于数据整合与分析。选取多景影像，覆盖林地与非林地两个下垫面，与优化后阈值与传统阈值进行对比分析，验证优化后阈值的积雪判别精度。

1.2.4　土地覆盖数据

土地覆盖数据选用MCD12Q1 V006产品数据。数据来源于NSIDC网站（https：//search.earthdata.nasa.gov/search）该产品使用MODIS Terra和Aqua反射数据进行监督分类，按年间隔（2001—2016年）提供全球土地覆盖类型。利用年度国际地圈生物圈计划（IGBP）分类，共分为17个类，包括：常绿针叶林、常绿阔叶林、落叶针叶林、落叶阔叶林、混交林、稠密灌丛、稀疏灌丛、稀疏草原、热带稀疏草原、草地、永久湿地、耕地、城市和建筑区、自然植被、稀疏植被、雪和冰、水体。通过对各土地覆盖类型所占比例分析，所占比例最大的类别为草原，占比超过50%；其次为稀疏植被，即植被不足10%的区域；其余类别占比均少于10%，本文使用2016年度MCD12Q1数据将青藏高原地区东南部包含常绿针叶林、常绿阔叶林、落叶针叶林、落叶阔叶林、混交林、稠密灌丛、稀疏灌丛、自然植被、稀疏植被9种林地类型，统一归属于林地类型；其余地物类型，归属于非林地类别^［39］（图1）。

2 研究方法

图2为本文NDSI阈值优化算法流程图，主要分为以下几个步骤：

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 NDSI阈值优化算法流程图

Fig. 2 Flow chart of NDSI threshold optimization algorithm

①利用GEE云平台获取地表覆盖类型MCD12Q1数据^［40］、MODIS V006版本下MOD10A1、MYD10A1的NDSI和NDSI Snow Cover Class数据，合成逐日无云产品；获取地面气象站点雪深数据；

②利用地面气象站点雪深数据和逐日无云产品数据通过NDSI与雪深对应时序图确定阈值范围，在阈值范围内计算总体精度确定林地和非林地的最优阈值；

③使用确定最优阈值制作二值积雪分布图；

④利用Landsat 8数据提取“真值”，对二值积雪分布图进行精度验证。

2.1　逐日无云NDSI影像生成

采用Chen等^［35］提出的STAGFM方法生产逐日无云影像。具体去云过程如下：

首先对MOD10A1和MYD10A1的NDSI数据层进行合成得到MOYD。优先采用MOD10A1的有效数据，即如果MOD10A1为无云像元，则取该值作为MOYD的值；若MOD10A1为云遮挡，MYD10A1无云，则采用MYD10A1为MOYD的值；若均为云遮挡，则该像元仍然为云像元，等待下一步去云操作。

然后利用临近时间的数据进行前后时间段的去云处理。本文选用了前后两天的数据进行融合，结果记为ATC（adjacent temporal composite）。如果选择的时间间隔太长，可能会由于积雪的变化导致去云结果存在较大的误差。如果仅采用前后一天的数据，则剩余的云量太多导致最后一步去云运算量过大。

最后采用基于相似像元的加权运算，对剩余的云遮挡像元进行去云处理。该算法主要采用上下午星结合；临近日数据判别，临近日选取前后各两天数据，最后采用空间相似像元进行加权运算，得到逐日无云影像。

2.2　去云精度验证

图2中最后一列为相似像元加权去云影像；为了进一步定量验证去云的精度，本文采用“云掩膜”的方法。“云掩膜”方法是指在无云的影像上添加云掩膜，然后对添加云掩膜的影像进行去云处理，再将去云结果与真实的影像进行对比来评价算法的好坏。为了使得添加的云接近真实的云分布，本文选取了一个水文年中云量最小的五景影像作为NDSI影像的真值，将本月中云量最大的云覆盖到该NDSI影像上，这里采用的是MODIS NDSI SNOW COVER中的云类别（250）对MODIS NDSI进行云掩膜。利用相关系数r，均方根误差RMSE，和平均偏差MAE这几个指标对去云结果进行精度评价。

2.3　利用典型气象站点数据确定阈值

2.3.1　阈值范围确定

选取了11个典型气象站点，见图1中的A~H。另外，图1中A~H为非林地，I~J为林地。建立典型气象站点（林地与非林地）去云NDSI序列与站点雪深对应图，其中A~H站点数据为2017年11月1日至2018年3月31日共151天数据。图中横轴day1代表该年11月1日。

图3（a）为气象站点A所在像元的NDSI序列与站点雪深对应图，该像元位于中国西藏自治区日喀则市亚东县，为纯净的裸地像元。可以看出在积雪稳定期，即2017年12月30日—2018年1月15日，NDSI均大于0.2；且在2017年10月10日左右出现了明显降雪，雪深为4 cm时，NDSI大于0.2，出现明显峰值。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 典型地表的NDSI序列与雪深变化曲线（A~J站点见图1）

Fig. 3 NDSI series and snow depth variation curve of typical surface： the meteorological station A~J （a）~（j）， and A~J see the Fig.1

图3（b）为气象站点B所在像元的NDSI序列与站点雪深对应图，该像元位于中国青海省海南藏族自治州同德县，为纯净的耕地像元。可以看出在出现积雪时，NDSI均出现峰值，峰值NDSI多大于0.2。可以看出2018年3月20日前后均有积雪出现，2018年3月20日无积雪出现，NDSI出现了明显的低谷，在2018年3月20日时，NDSI为-0.4，与积雪存在日期的NDSI值区分显著。

图3（c）为气象站点C所在像元的NDSI序列与站点雪深对应图，该像元位于中国西藏自治区山南市，为雅鲁藏布江与居民地混合像元。该像元NDSI与雪深同样具有良好对应关系，出现降水的三个时期，NDSI的峰值分别为0.2、0.4、0.4。

图3（d）为气象站点D所在像元的NDSI序列与站点雪深对应图，该像元位于中国青海省果洛藏族自治州达日县，为纯净的居民地像元。在2017年11月25日至2017年12月5日之间出现稳定积雪时，NDSI均大于0.2。

图3（e）为气象站点E所在像元的NDSI序列与站点雪深对应图，该像元位于中国四川省甘孜藏族自治州石渠县，为山地与居民地混合像元。在出现积雪时，NDSI值多大于0.2。2018年3月20日附近，积雪为1 cm，NDSI没有显著变化，由于青藏高原积雪太浅，探测不够显著。

图3（f）为气象站点F所在像元的NDSI序列与站点雪深对应图，该像元位于中国西藏自治区昌都市左贡县，为裸地与草地混合像元。该气象站点分别于2017年11月3日、2018年1月15日、2018年3月22日共出现三次积雪。在三次积雪时间内，NDSI均处于0.2~0.4之间。

图3（g）为气象站点G所在像元的NDSI序列与站点雪深对应图，该像元位于中国西藏自治区那曲市嘉黎县，像元下垫面类型均为山地。在2018年1月30日的短时间积雪与2018年2月10日至2018年2月30日之间的长时间积雪中，NDSI与雪深均有良好的对应关系，NDSI的值大于0.2。

图3（h）为气象站点H所在像元的NDSI序列与站点雪深对应图，该像元位于中国四川省阿坝藏族羌族自治州小金县，为混合像元主要下垫面类型为居民地。2018年2月10日出现积雪，NDSI与其对应较好，NDSI值大于0.2。

图3（i）为气象站点I所在像元的NDSI序列与站点雪深对应图，该像元位于中国云南省迪庆藏族自治州香格里拉市，为林地与居民地混合像元，为了更合理地统计森林像元NDSI与雪深的对应关系，选取距离原站点距离附近的纯森林像元（距离小于5 km），站点I在2018年1月10日至2018年1月20日出现连续积雪时，NDSI均大于0。

图3（j）为气象站点J所在像元的NDSI序列与站点雪深对应图，该像元位于中国四川省甘孜藏族自治州稻城县，为林地与居民地混合像元，为了更合理地统计森林像元NDSI与雪深的对应关系，选取距离原站点距离附近的纯森林像元（距离小于5 km），在积雪季中出现四次积雪雪深，其中三次与NDSI对应效果较好，此时NDSI大于0。

因此在包括下垫面为裸地、耕地、河流、居民地、草地、山地的各类非林地区域，NDSI与雪深均有较好的对应关系。在出现积雪的时期内，NDSI多大于0.2，因此在非林地区域，本文选取0.2以上作为积雪时期时NDSI的阈值范围。

在选取的各气象点附近选取的纯森林像元，NDSI与雪深均有较好的对应关系。在出现积雪的时期内，NDSI所在区间范围多大于0，因此在林地区域，本文选取0以上作为积雪时期时NDSI的阈值范围。

2.3.2　最优阈值确定

由图3确定不同下垫面阈值选取范围，即林地阈值范围为NDSI大于0，非林地阈值范围为NDSI大于0.2；在阈值范围内，将0.01作为分界点，计算每个分界点的总体精度［式（1）］、低估误差［式（2）］和高估误差［式（3）］以确定最优阈值。

O A = a + d / a + b + c + d

（1）

U E = b / a + b + c + d

（2）

O E = c / a + b + c + d

（3）

式中：a表示NDSI与雪深均为有雪像元；b表示雪深为有雪像元，NDSI为无雪像元；c表示NDSI为有雪像元，雪深为无雪像元；d表示NDSI与雪深均为无雪像元。

（1）非森林像元

非森林像元选取为编号为1~15共15个典型气象站点数据（图1），雪深大于等于1 cm则记为有雪。随着NDSI的变化，总体精度OA，高估误差OE，低估误差UE的变化，如图4所示。当NDSI等于0.26时，总体精度OA最高，为94.46%；在该NDSI之下，高估误差和低估误差分别为0.51%和5.03%。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 非林地积雪识别精度随NDSI阈值的变化曲线

Fig. 4 Variation curve of non-forest snow recognition accuracy with NDSI threshold： non-forest OA （a）， non-forest OE （b） and non-forest （c）

（2）森林像元

森林像元选取编号为16~24共9个典型气象站点数据（图1），其中点编号16为2013年森林像元；点编号17~19为2014年森林像元；点编号20~22为2015年森林像元，点编号23~24为2016年森林像元。由于所选取典型气象站点下垫面为混合像元，且站点的雪深数据可以代表周围森林像元的雪深。因此我们通过Google Earth目视解译选取在该雪深监测点附近（距离小于5 km）查找纯像元下垫面。对其进行最优阈值的选取计算，如图5所示。当NDSI等于0.03时，总体精度OA最高，为94.23%；在该NDSI之下，高估误差和低估误差分别为1.21%和4.6%。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 林地积雪识别精度随NDSI阈值的变化曲线

Fig. 5 Variation curve of forest snow recognition accuracy with NDSI threshold： forest OA （a）， forest OE （b） and forest （c）

3 积雪面积数据精度验证

3.1　积雪产品精度评估

利用GEE云平台对Landsat 8 OLI“真值”影像对MODIS产品数据进行精度验证。由于Landsat 8 OLI影像分辨率30 m，MODIS影像分辨率为500 m。因此对Landsat 8 OLI影像进行重采样至500 m，方便进行精度对比。以积雪面积比例≥50%提取积雪分类精度。如图所示，其中图6（a）、图7（a）分别为林地、非林地区域利用本文所选取精度提取的积雪，即森林地区选用0.03作为阈值，非森林地区选用0.26作为阈值，选取的Landsat 8 OLI数据的空间位置如图1所示。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 林地区域Landsat 8数据与对应V6影像对比

Fig. 6 Comparison of Landsat 8 data and corresponding V6 images in woodland area

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 非林地区域Landsat 8数据与对应V6影像对比

Fig. 7 Comparison of Landsat 8 data and corresponding V6 images in non-forested areas

利用6景Landsat 8 OLI数据，其中包括1景林地区域数据，5景非林地区域数据（图1 Landsat 8定量精度验证），进行验证得到，在传统阈值下总体精度OA为82.86%，高估误差OE为1.48%，低估误差UE为15.66%；在新的阈值下总体精度OA为84.21%，高估误差OE为5.33%，低估误差UE为10.46%。在新的阈值下低估误差UE下降明显，高估误差OE由于阈值选取值变低，有了一定程度升高，总体精度OA得到了提高，证明新的阈值提高了在青藏高原地区积雪的提取精度。

图6中Landsat 8 OLI为Landsat 8时相为2018年1月11日行列号为131/037的影像“真值”图像（图6中左图）。图6（a）为林地区域（蓝框区域）新的阈值下的积雪提取，图6（b）为林地区域（蓝框区域）Landsat 8 OLI“真值”，图6（c）为林地区域（蓝框区域）传统阈值0.4下的积雪提取。结果表明，与传统阈值相比，新阈值下的积雪提取更接近真值。传统阈值与“真值”对比可以看出，传统阈值更易低估积雪区域为非积雪区域（图6蓝框区域）。新的阈值虽然也存在一定情况的高估低估现象，但相对传统阈值有了明显提高，可以看出新的阈值能更有效地识别积雪。

图7中Landsat 8 OLI为Landsat 8时相为2018年1月10日行列号为138/036的影像“真值”图像（图7中左图）。图7（a）为非林地区域（蓝框区域）新的阈值下的积雪提取，图7（b）为非林地区域（蓝框区域）Landsat 8 OLI“真值”，图7（c）为非林地区域（蓝框区域）传统阈值0.4下的积雪提取。结果表明，与传统阈值相比，新阈值下的积雪提取更接近真值。图中部积雪大范围集中，新的阈值与传统阈值提取效果均相对较好；图像右下角区域，积雪相对分散破碎，新阈值能提取出大量积雪，传统阈值则不能。在对传统阈值进行调整后，新的阈值虽然也存在一定情况的高估低估现象，但相对传统阈值有了明显提高，可以看出新的阈值能更有效地识别积雪。

4 结论与讨论

本文主要考虑不同下垫面对青藏高原地区积雪识别中NDSI阈值选取的影响。通过对林地、非林地区域的NDSI阈值优化，提高不同土地覆盖类型下积雪面积判别精度，主要结论如下：

（1）本文采用的MODIS数据，并对生成的逐日无云积雪数据进行精度验证，在5景影像的去云结果中，相关系数r均大于0.8，相关系数均值大于0.86。均方根误差RMSE均小于0.15，平均绝对误差MAE均小于0.1，证明反演结果具有可靠性和合理性。积雪季数据去云NDSI序列与站点雪深在下垫面为林地和非林地的区域均有良好的对应关系。NDSI可以准确反映站点所在像元的降雪、融雪现象。

（2）以实测站点数据与MODIS V006遥感影像进行对照训练，其中林地采用2013—2017年积雪季数据，非林地采用2017年积雪季数据，在青藏高原的林地与非林地下垫面的最优NDSI阈值，分别为0.03、0.26。

（3）土地覆盖类型条件是影响NDSI阈值选取的主要因素。以Landsat 8 OLI数据作为“真值”进行“面”的评估定量验证时，得到MOD_0.4_V6、MOD_OPT_V6在整个研究区的OA分别为82.86%、84.21%，优化后的NDSI阈值在一定程度上改进了积雪识别精度。在以Landsat 8 OLI数据作为“真值”进行定性验证时，在积雪大范围集中的区域内，新的阈值与传统阈值提取效果均相对较好；而在积雪相对分散破碎，新阈值能提取出大量积雪，传统阈值则不能。可以看出在对传统阈值进行调整后，新的阈值虽然也存在一定情况的高估低估现象，但相对传统阈值有了明显提高，可以看出新的阈值能更有效的识别积雪。

本文主要考虑不同下垫面对青藏高原地区积雪识别中NDSI阈值选取的影响。通过对林地、非林地区域的NDSI阈值优化，与传统阈值相比，优化后阈值的积雪判别精度得到了提高，为NDSI在积雪识别中的应用提供了有力的支撑。

尽管如此，利用MODIS V006版本数据所选取的阈值的整体流程还有许多待完善的地方。首先，在NDSI逐日无云影像生成的过程中，NDSI值的填补可能出现误差。其次，由于气象站点数据为点尺度数据，而遥感影像对地面进行了栅格化，直接和地面点的数据进行验证会存在尺度上的偏差。最后本文选取的林地气象站点海拔均为3 000~4 000 m之间，但由于林地气象站点多修建于混合像元区域，为保证下垫面的均一性，在Google Earth进行了目视判别确定选取的NDSI区域，保证了林地下垫面的均一性。同时由于林地像元位置与气象站点位置并不完全一致，林地的保温作用，可能会导致在积雪较浅且不稳定的时期内，气象站点位置融雪快于林地像元位置融雪，增大高估误差，影响总体精度。因此需再使用精度更高的遥感数据进行提取验证，当然，尝试进行亚像元尺度的积雪覆盖信息提取，也将是改善积雪有效监测的另一种途径。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Cao

Meisheng

， Li

Xin

， Chen

Xianzhang

， et al. Cryosphere remote sensing［M］. Beijing： Science Press， 2006.