基于AVHRR影像的北半球积雪识别算法

王京达， 郝晓华， 和栋材， 王 建，3， 李弘毅， 赵 琴

（1.太原理工大学矿业工程学院，山西太原 030024； 2.中国科学院西北生态环境资源研究院，甘肃兰州 730000；3.江苏省地理信息资源开发与利用协同创新中心，江苏南京 210023）

0 引言

积雪作为冰冻圈的重要组成部分，在地球表面分布较广，大约98%的季节性积雪位于北半球，面积最大可达4.7×107km2[1］，且变化较为活跃，是影响气候变化的重要因子[2-4］。在全球水循环中，积雪的积累和消融过程调节着水的年内再分配，是干旱、半干旱地区春季最重要的淡水资源。积雪是气候变化的指示器[5-6］，积雪范围、积雪深度及雪水当量等信息是水文模型和气候模式的重要输入参数[7-8］，全球尺度下的积雪参数监测对于气象、水循环、气候变化以及能量平衡领域的研究具有重要的意义[9］。因此，准确获取北半球的积雪信息可为寒区水文循环、水资源管理、雪灾预警等重大研究提供必要可靠的参考资料[10］。

北半球积雪范围和气候变化研究需要长时间序列、高时空分辨率、高精度的积雪范围数据。地面观测积雪范围资料是目前研究积雪长期变化特征最广泛、最可靠的资料，具有时间尺度长、连续性好的优点，但地面观测资料无法准确获取积雪范围[11-12］。可用于北半球研究的积雪范围产品有微波产品和光学产品，微波积雪范围产品主要来源于SMMR、SSM/I[13-14］和AMSR-E[15-16］传感器。微波系列产品时间序列较长，但空间分辨率较低，目前可用的长时间序列微波积雪范围产品的空间分辨率为25 km。 光学积雪范围产品主要来源于AVHRR[17］传感器、MODIS[18-19］传感器，以及Landsat系列卫星[20］和哨兵（Sentinel）系列卫星[21］。MODIS积雪范围产品（2000 年至今）时空分辨率较高，精度也较高，广泛用于全球的水文和气候变化研究，但对于气候变化研究来说，时间序列仅有约20 年，较短的时间序列无法满足研究；Landsat 时间序列长，空间分辨率高，但时间分辨率低（16 d），无法准确监测季节性积雪的变化；Sentinel系列卫星虽然空间分辨率高，但时间序列短，时间分辨率较低，也无法满足气候变化研究需求。AVHRR 产品具有时间序列长（1981 年至今）、时空分辨率较高的特点（时间分辨率为1 d，空间分辨率为5 km），是研究北半球气候变化的理想数据。

一些学者利用AVHRR 发展了积雪范围产品。Hüsler 等[22］利用1 km 分辨率的AVHRR 原始数据集，开发了一种改进的积雪范围检索算法来识别欧洲阿尔卑斯山的积雪范围；Zhou 等[23］利用决策树算法基于1 km 分辨率的AVHRR LAC 数据和HRPT数据，以中亚地区为例开发了一种长时间序列积雪范围的识别算法。但此类研究主要集中在小区域的积雪识别，算法不具有全球普适性。日本宇航局（Japan Aerospace Exploration Agency，JAXA）利用AVHRR 和MODIS 结合发展的北半球1978—2015年5 km 分辨率的JASMES（JAXA Satellite Monitor⁃ing for Environmental Studies）积雪范围产品，弥补了大尺度长时间序列逐日积雪范围产品数据的缺失，但是该产品经验证发现在青藏高原精度较低[24］。Hao等[25］在中国区域发展了多级决策树积雪识别算法，并生成了一套长时间序列积雪范围产品，但算法中的阈值对于全球产品是否适用仍需要探索。

本研究的目的是利用AVHRR 数据，获取多级决策树积雪识别算法中适合北半球的阈值，生成1981—2000 年北半球历史积雪范围数据集，从而为气候变化研究提供更长时间序列、更高精度的逐日积雪范围数据集。研究主要解决两方面问题：①如何通过有限的6个波段，结合研究区给定合理阈值；②如何填补AVHRR 在高纬度地区数据完全缺失和低纬度地区数据部分缺失问题，即如何填补数据空缺值。主要研究内容包括：①以高空间分辨率的无云Landsat-5 TM 参考积雪范围影像作为参考真值，获取北半球积雪识别算法合理阈值；②填充积雪范围数据中的空缺值，制作北半球AVHRR 逐日积雪范围数据集；③利用北半球气象台站的雪深数据和Landsat-5 TM 作为参考真值，验证数据集精度。该数据集将为北半球积雪变化研究提供重要支撑。

1 数据来源

1.1 AVHRR地表反射率数据

NOAA-AVHRR 地表反射率数据第四版（Sur⁃face Reflectance data Version 4，SR V4）被用作本研究的主要数据，来源于美国国家环境信息中心（https：//www. ncdc. noaa. gov/），它 是 由NOAA GAC（Global Area Coverage）L1 数据经过一系列的定标、辐射校正、大气校正及几何校正后得到的地表反射率数据和亮温数据。该数据时间范围为1981年6月24日至2019年5月16日，时间分辨率为1 d，空间分辨率为5 km。本研究使用1981 年6 月24日至2000年2月27日时间范围内的数据，数据被转换投影坐标为椭球体WGS84 Albers 投影坐标，包含各个波段的地表反射率、亮温和质量控制标记（quality control bit flags assessment，QA）等，具体如表1所示。

表1 NOAA-AVHRR地表反射率数据Table 1 NOAA-AVHRR surface reflectivity data

1.2 Landsat-5 TM数据

Landsat-5 TM 专 题 成 像 仪（thematic mapper，TM）是搭载在Landsat-5陆地观测卫星的传感器，它于1984 年3 月发射成功直到2013 年6 月退役。本研究通过SNOWMAP 算法[18，26］获取Landsat-5 TM二值积雪范围，即雪和非雪。SNOWMAP 算法的核心是通过归一化雪被指数（normalized difference snow index，NDSI）区分积雪和非雪，计算公式为

式中：NDSIL为Landsat-5 TM 的归一化雪被指数；Band2 为Landsat-5 TM 影像的第二波段（0.52~0.60 μm）；Band5 为第五波段（1.55~1.75 μm）。算法设定NDSIL≥0.4 来区分雪和非雪，同时为了排除水体的干扰，设定第四波段（0.76~0.90 μm）的阈值大于0.11。利用该算法获取Landsat-5 TM 空间分辨率为30 m 的二值积雪范围后，通过最近邻法将空间分辨率30 m 的Landsat-5 TM 二值积雪范围影像重采样为5 km 分辨率，转换过程中设定5 km 像元中30 m 积雪像元占比≥50%作为积雪像元，否则为非雪像元，投影与AVHRR 保持一致，最终将获取的5 km 分辨率Landsat-5 TM 二值积雪范围影像作为参考真值。

Landsat-5 TM 参考积雪范围影像分为两组：第一组参考真值用作训练AVHRR 积雪识别算法阈值，选取1986年、1987年、1996年、1997年和1998年积雪季（11 月1 日开始至次年3 月31 日）的无云Landsat-5 TM 作为参考积雪范围影像，共1 094 景，训练数据共包含积雪像元数为156 661 个，非雪像元数为696 296 个；第二组参考真值用于评估AVHRR 积雪范围数据集，筛选1988 年、1989 年和1990 年积雪季Landsat-5 TM 无云且积雪占比大于50%的影像作为验证数据集，共939景，验证数据分布如图1所示。

图1 北半球气象台站和Landsat-5 TM验证数据集的空间分布Fig. 1 Spatial distributions of meteorological stations and Landsat-5 TM validation data in the Northern Hemisphere

1.3 全球长时间序列逐日雪深数据集

本研究利用微波雪深数据集来填补AVHRR 中的空缺值，使用Che 等[27］采用SMMR 和SSM/I2 生成的全球长时间序列逐日雪深数据集。该微波雪深数据集，采用动态亮温梯度算法，通过建立不同频率亮温差与实测雪深在空间和季节上的动态关系[28-29］获取雪深。该数据集来源于时空三极环境大数据平台（http：//poles. tpdc. ac. cn/zh-hans/）。本研究使用的时间范围为1981 年6 月24 日至2000 年2月27日，数据集空间分辨率为25 km，利用最近邻法将该数据集重采样为5 km，并转换为与AVHRR 相一致的投影。

1.4 地面雪深观测数据

雪深数据用于对AVHRR 积雪识别算法和数据进行评估，本研究使用的北半球雪深数据来自于NOAA 的国家环境信息中心网（http：//www. ncdc.noaa.gov/），是由世界气象组织（World Meteorologi⁃cal Organization，WMO）提供的气象台站数据。本研究使用1981—1999 年2 546 个气象台站积雪季的实测雪深数据，使用的气象台站分布如图1所示。

1.5 辅助数据

本研究共使用两种辅助数据，一种是数字高程模型，另一种是气候再分析数据-地表温度产品。数字高程模型数据（Digital Elevation Model，DEM）为SRTM（Shuttle Radar Topography Mission）90 m数字高程数据库第4 版本数据，是由美国国家航天局联合美国国防部国家地理情报局通过航天飞机雷达测绘的数据产品（http：//srtm. csi. cgiar. org/），本研究中将STRM 数据用作云雪区分和积雪识别算法辅助数据，也用来作为消除积雪伪像元的辅助数据。

ERA5 气候再分析数据-地表温度逐日产品[30］是由欧盟资助欧洲中距离天气预报中心（European Center for Medium-Range Weather Forecast，EC⁃MWF）执行哥白尼气候变化服务推出的适用于气候研究的数据集。本研究中将ERA5数据集也作为辅助数据用来消除可能的积雪伪像元，空间分辨率为10 km，同样采用最近邻法将上述数据重采样为5 km，并转换为与AVHRR相一致的投影。

2 方法

本研究AVHRR 北半球积雪识别算法流程如图2所示，主要分为以下四个步骤：

图2 北半球AVHRR积雪识别算法流程Fig. 2 Flowchart of AVHRR snow discrimination algorithm in the Northern Hemisphere

（1）数据预处理

根据AVHRR SR V4 中自带的质量控制标记剔除反射率无效值。

（2）云识别算法

针对冰云和积雪容易混淆的问题，结合云雪混淆区分技术，提高了积雪识别精度。

（3）积雪识别算法

利用大量的样本进行训练，针对北半球积雪研究区的特性，获取多级决策树积雪识别算法的最优阈值，建立适合北半球的积雪识别算法，通过云、雪识别算法生成AVHRR L1级积雪范围数据集。

（4）空缺值插补

针对AVHRR 高纬度地区数据完全缺失和低纬度地区部分数据缺失问题，利用微波雪深数据集填充空缺数据。从而生成了一套北半球AVHRR 1981—1999 年逐日积雪范围数据集，即AVHRR L2 级积雪范围数据集。

2.1 数据预处理

根据AVHRR SR V4 中自带的质量控制标记，如表2 所示，剔除AVHRR 各个通道的反射率无效值。有效值被作为下一步云、雪、非雪的区分，无效值被认为是空缺值。

表2 AVHRR SR V4的质量控制标记描述Table 2 Descriptions of quality control bit flags assessment of AVHRR SR V4

2.2 云识别算法

关于云的识别算法，很多科学家也做了研究[10，31］，Hao 等[25］针对AVHRR SR V4 产品发展了一套云检测算法，该算法可以有效解决云雪混淆问题。因此本研究结合Hao 等[25］生成的AVHRR 云雪混淆区分技术进行云检测，算法如表3所示。

表3 云检测及其阈值Table 3 Cloud detection tests and their threshold values

2.3 积雪识别算法

由于Hao 等[25］开发的多级决策树积雪识别算法是面对中国地区来提取积雪，在北半球算法的阈值还有待改进。因此本研究采用该算法的方案（图3），结合Landsat-5 TM 参考积雪范围影像作为参考真值，来获取北半球积雪识别算法最优阈值。根据前人对AVHRR 积雪识别的分类研究[10，32］，涉及雪阈值算法的变量主要有Ref1、（Ref3/Ref2）和BT11三个一级阈值变量，NDVI 和（Ref3−Ref2）两个二级阈值变量，NDSI 一个三级阈值变量，其中NDSI 的区分度最好。NDSI 通常使用绿色（大约0.50 μm波长）和短波红外波段（大约1.60 μm 波长）计算。由于在AVHRR SR V4 中没有1.60 μm 附近的短波红外观测，参考Hori 等[10］使用3.75 μm 波长的Ref3用 于NDSI 的 计 算，NDSI、NDVI 可 由 下 式 计 算得到。

（1）一级指标阈值

如图3 所示，一级指标为Ref1、（Ref3/Ref2）、DEM 与BT11。一级指标确立原则是排除大量的非雪数据，使尽可能多的雪被识别进而再用下一级指标进行区分。一级指标由于仅仅是初级筛选，因此一级指标阈值采用Hao 等[25］算法的阈值。表4 展示了一级指标和对应的阈值，即当Ref1＞0.14，（Ref3/Ref2）＜0.32、DEM＜1 300 m 并且BT11＜280 K 或者DEM≥1 300 m 并且BT11＜281.5 K 时，则认为可能为积雪，否则为非雪，进行下一级别区分。

表4 北半球积雪识别算法最优阈值选取Table 4 Selection of optimal threshold for snow discriminant algorithm in the Northern Hemisphere

图3 多级决策树积雪识别算法流程Fig. 3 Flowchart of snow discriminant algorithm based on multi-level decision tree

（2）二级指标阈值

如图3所示，二级指标是NDVI和（Ref3−Ref2）。将一级指标筛选后可能为积雪的数据进行分段统计，将非雪样本1%各指标对应的反射率值作为二级指标阈值，通过计算：获取NDVI最优阈值-0.11，即满足NDVI＜-0.11 设定为第一类积雪（雪1）；获取（Ref3−Ref2）最优阈值-0.96，即满足（Ref3−Ref2）＜-0.96 设定为第二类积雪（雪2），不符合条件的被认为可能为积雪，进行下一级别区分。

（3）三级指标阈值

如图3 所示，三级指标是NDSI。利用二级指标筛选后可能为积雪的数据，通过计算最佳总体精度的方法得出三级指标，即逐步长方式计算最佳总体精度（步长设为0.01），最佳总体精度对应的NDSI值即为三级指标。最优阈值的确定方法如图4 所示。图4（a）展示了积雪和非雪的NDSI 频率分布，选择具有最高总体精度OA（88%）时对应的NDSI值作为最优阈值（0.76）。计算结果显示NDSI＞0.76，符合条件的设定为第三类积雪（雪3），不符合条件的被判别为非雪。北半球积雪识别算法的最优阈值如表4所示。

图4 三级指标NDSI积雪和非雪频率分布（a）和总体精度（b）Fig. 4 Frequency distribution of third-level indicator NDSI with and without snow（a）and overall accuracy（b）

2.4 空缺值插补

微波数据具有不受云干扰的特点，可对AVHRR 的空缺值进行插补。利用Che 等[27］、Dai等[28］生成的全球被动微波数据逐日雪深数据集进行插补。插补分为两个步骤：

（1）将25 km 空间分辨率的雪深数据利用重采样方法（最近邻法）降尺度到5 km，投影和AVHRR一致。

（2）根据Che 等[27］研究结果，插补规则如下：①纬度处于0°~60° N 时，若AVHRR 的云像元对应当日雪深≥2 cm，则该云像元被分类为积雪像元；若雪深＜2 cm，则该云像元被分类为非雪像元。②纬度处于60°~90°N时，若所有像元对应当日雪深≥2 cm，则该像元被分类为积雪像元；若雪深＜2 cm，则该像元被分类为非雪像元。

3 结果与精度评估

3.1 结果

AVHRR L1 级和L2 级逐日积雪范围数据集示意图如图5 所示。图5（a）展示了通过云、雪识别后的北半球AVHRR L1 级逐日积雪范围数据集示意图，其中雪像元赋值为1，云像元赋值为255，陆地赋值为0。图5（b）是北半球AVHRR L2级逐日积雪范围数据集示意图，其中雪像元赋值为1，云像元赋值为250，陆地赋值为0。

图5 L1级（a）和L2级（b）数据集示意图（选取日期为1998-01-13）Fig. 5 Schematic diagram of datasets at Level 1（a）and Level 2（b）（The selected date is 1998-01-13）

3.2 精度评估方法

3.2.1 基于气象台站雪深

利用1981—1999 年间北半球2 546 个气象台站的雪深数据分别对AVHRR L1 级和L2 级积雪范围数据集进行精度评估。设定气象台站雪深大于等于1 cm 时认为该气象台站所在的像元为积雪像元，雪深小于1 cm 时认为该气象台站所在的像元为非雪像元[33］。将气象台站积雪信息与AVHRR 各级积雪范围数据集相应位置的积雪信息进行比对验证。利用气象台站雪深数据转化为雪和非雪的二值数据，并作为参考真值，然后分别构建混淆矩阵分析总体精度及误差。具体的SS、SN、NS 和NN 的定义参见表5，其中SS 指AVHRR 积雪范围数据集和参考真值两者都判别为雪的像元个数，NS 指AVHRR积雪范围数据集判别为有雪且参考真值判别为非雪的像元个数，SN 指AVHRR 积雪范围数据集判别为非雪且参考真值判别为雪的像元个数，NN 指两者都判别为非雪的像元个数。

表5 精度验证混淆矩阵的指标含义Table 5 Index implication of confusion matrix for accuracy validation

分析数据集的总体精度（OA）、用户精度（UA）、制图精度（PA）和F-Score（FS）[34］。其中，OA 为雪和非雪像元被正确分类的总体精度；UA 为气象台站雪深数据被分类为积雪的概率，作为用户精度；PA为AVHRR 积雪范围数据集被分类为积雪的概率，作为制图精度；FS 为UA 和PA 的调和平均值，是综合考虑了UA 和PA 的结果，当FS较高时说明UA 和PA较为均衡。

计算公式为

3.2.2 基于Landsat-5 TM

基于Landsat-5 TM 的验证同样使用3.2.1 节中构建混淆矩阵的方法，进一步评估AVHRR 积雪范围数据集（L2 级）的精度。将AVHRR 积雪范围数据集分别在北半球和中国、北美洲、欧亚大陆、西伯利亚、北极地区这五个北半球重要积雪区[35］进行精度评估，分析数据集的OA、UA、PA 和FS。（Landsat-5 TM参考积雪范围被当作参考真值）

3.3 验证结果与分析

3.3.1 基于气象台站雪深的精度评估

利用北半球1981—1999 年的共2 546 个实测雪深气象台站，分别对AVHRR L1 级和L2 级积雪范围数据集进行精度评估，计算结果如表6 所示，AVHRR L1 级数据集的总体精度为81.8%，AVHRR L2 级数据集的总体精度为82.2%，总体验证精度较高。如图6 所示，总体精度年际变化比较稳定，都在80%左右，其中1996 年精度最高达到85.9%。AVHRR L1 级数据集的UA 为83.7%，PA为81.7%，AVHRR L2级数据集的UA为83.8%，PA为84.2%，表明AVHRR L1 级数据集漏分误差略大于错分误差，尤其是在1993 年，由于AVHRR 大面积的数据空缺，导致漏分误差较大，而AVHRR L2级数据集，制图精度有所提高，表明雪深插值在保证精度的情况下，不仅可以填补空缺数据，而且可以有效降低漏分误差，预测正确的积雪像元个数在AVHRR 识别的积雪像元中和气象台站记录的积雪信息中的比例较高。FS 是平衡UA 和PA 的指数，AVHRR L1级和L2级数据集的FS 分别为82.9%和85.1%，FS 较高，表明漏分误差和错分误差都比较均衡，算法可有效解决积雪识别问题。

图6 L1级（a）和L2级（b）数据集的OA、UA和PA年际变化Fig. 6 Interannual variations of OA，UA and PA of datasets at Level 1（a）and Level 2（b）

表6 气象台站验证评估结果Table 6 Verification and evaluation results of meteorological stations

3.3.2 基于Landsat-5 TM的精度评估

在空间上，利用1.2 节中的筛选方法，获取了939 景Landsat-5 TM 参考积雪范围影像作为参考真值，其中中国155 景、北美洲382 景、欧亚大陆455景、西伯利亚174 景及北极地区63 景。使用3.2.2节中的验证方法，对待检验的AVHRR 积雪范围数据集（L2级）进行评估。

表7 展示了北半球和五个北半球重要积雪区（中国、北美洲、欧亚大陆、西伯利亚及北极地区）的验证评估结果。结果表明，AVHRR 积雪范围数据集的总体验证结果OA、UA 和PA 分别为90.3%、90.2%和99.1%，该北半球积雪范围数据集的总体精度、制图精度和用户精度都较高，数据具有一定可靠性。其中PA 最高达到99.1%，说明漏分误差较低，预测正确的积雪像元个数在AVHRR 识别的积雪像元中和Landsat-5 TM 识别的积雪像元中的比例较高。FS 为93.4%，是UA 和PA 的平衡，较高则说明漏分误差和错分误差较为均衡，积雪识别算法比较有效。对中国、欧亚大陆和西伯利亚地区的总体验证结果OA、UA、PA 和FS 都在90.0%以上，北美洲地区总体精度OA 为85.2%，总体精度在中纬度地区较高，而在北极地区，OA 为73.7%，UA 为73.7%，PA 为99.9%，说明错分误差大于漏分误差，且漏分误差为0.1%，这是由于AVHRR 在高纬度地区数据的完全缺失，使用2.4 节中的数据和方法对其进行插补，该数据集是通过微波数据和气象台站数据融合的数据集，说明该数据集错分误差较高，漏分误差较低。

表7 Landsat-5 TM验证评估结果Table 7 Verification and evaluation results of Landsat-5 TM

5 景获取的Landsat-5 TM 验证样本的详细信息如表8 所示，OA、UA、PA 和FS 的平均值分别为93.4%、93.7%、99.2%和96.2%，数据总体验证精度较高，具有一定可靠性。每景验证样本对应的Landsat-5 TM 参考积雪范围、AVHRR 积雪范围，以及AVHRR 和TM 的两者对比如图7 所示。说明在大范围覆盖的积雪区域，AVHRR 积雪范围与Land⁃sat-5 TM 参考积雪范围是较吻合的，而在积雪边界区，AVHRR 积雪识别算法有未能识别出的Landsat-5 TM 参考积雪范围影像的积雪像元，也有更多识别出的积雪像元，这可能是因为AVHRR 积雪范围数据集对较低积雪范围比例的像元识别能力较不准确。总体而言，在积雪范围识别上算法存在一定的漏分和错分情况，但算法在保持漏分误差和错分误差较为均衡的前提下，总体精度较高，数据具有可靠性。

图7 Landsat-5 TM参考积雪范围影像与AVHRR积雪范围数据集对比Fig. 7 Comparison between Landsat-5 TM snow map and AVHRR snow cover extent datasets：China（a），North America（b），Eurasia（c），Siberia（d）and Arctic region（e）

表8 Landsat-5 TM验证样本信息Table 8 Verification sample information of Landsat-5 TM

4 结论

本研究基于AVHRR 地表反射率数据，获取了基于多级决策树积雪识别算法的各指标最优阈值，并结合云雪混淆区分技术，生成了北半球AVHRR 1981—1999 年L1 级逐日积雪范围数据集。利用微波雪深数据插值填补了AVHRR 的空缺数据，生成了北半球AVHRR 1981—1999 年L2 级逐日积雪范围数据集。最后利用气象台站实测雪深数据对L1级和L2级数据集分别进行精度验证，利用高空间分辨率遥感影像对L2级数据集进一步进行精度评估。

通过气象台站雪深数据对L1 级积雪范围数据集验证，评估结果表明：数据集总体精度为81.8%，用户精度为83.7%，生产者精度为81.7%，数据集总体精度较高，漏分误差比错分误差略高。FS 为82.9%，表明漏分误差和错分误差比较均衡，整体上，北半球积雪识别算法精度较高。

通过气象台站雪深数据对L2 级积雪范围数据集验证，评估结果表明：数据集总体精度为82.2%，用户精度为83.8%，生产者精度为84.2%，FS 为85.1%，L2 级数据集总体精度较高，漏分误差和错分误差相近。同时结合Landsat-5 TM 参考积雪范围影像，在面尺度上对L2 级数据集进行进一步验证，结果表明：L2 级数据集总体精度为90.3%，用户精度为90.2%，生产者精度为99.1%，在面尺度上，AVHRR L2 级积雪范围数据集展现了更高的精度，表明插值方法可以有效解决AVHRR 在高纬度地区数据的完全缺失和低纬度地区数据的部分缺失问题。该算法将为北半球积雪范围提取提供新的方案，L1 级和L2 级数据集也为北半球气候变化研究提供了重要数据支撑。