中国2000-2020年积雪反照率遥感产品

肖鹏峰，胡瑞，张正，秦棽，李震

1.南京大学地理与海洋科学学院，自然资源部国土卫星遥感应用重点实验室，江苏省地理信息技术重点实验室，南京 210023

2.中国科学院空天信息创新研究院，北京 100094

引 言

地表反照率指地表向各个方向反射的全部光通量与总入射光通量的比。它量化了大气与地表之间辐射的相互作用，是重要的气候参数[1-2]，被广泛用在地表能量平衡、中长期天气预测和全球变化研究中[3]。积雪具有显著的高反照率特性[4]，在局部或全球能量收支平衡中起着重要作用。近几十年来北半球年平均辐射强迫每年增加约0.45 W/m2，其中积雪带来的影响将近一半[5]。积雪反照率调节着积雪能量收支，进一步控制积雪的消融过程，从而影响区域乃至全球的水循环过程。同时，积雪为世界上17%的人口提供饮用水及生产生活用水，在部分山区流域积雪融水占到年径流量的50%以上[6]，是绿洲灌溉的宝贵水源[7]。此外，积雪污化物与环境污染相关，其浓度上升会导致积雪反照率显著下降[8]。因此，积雪污化物的辐射强迫研究对人类污染排放与大气活动有一定的指示作用[9]。

目前，国内外主要有4种遥感反照率产品：美国宇航局16天合成地表反照率产品MCD43A3、逐日积雪反照率产品SAD、欧洲空间局八天合成地表反照率产品GlobAlbedo、北京师范大学八天合成地表反照率产品 GLASS。MCD43A3产品是中分辨率光谱成像仪（Moderate Resolution Imaging Spectroradiometer，MODIS）反照率产品中的一部分，采用的反演模型为半经验核驱动模型（RossThick-LiSparse-Reciprocal，RTLSR），根据体散射核、几何光学核和各向同性核的加权和来估算地表BRDF[10-11]。SAD产品是MODIS积雪产品MOD10A1/MYD10A1中的一部分，在球形雪粒假设下，利用离散纵坐标辐射传输（Discrete Ordinates Radiative Transfer，DISORT）模型构建等效光学粒径为250 μm的各向异性响应函数（Anisotropic Response Function，ARF）查找表，经过窄宽转换计算得到积雪黑空反照率[12]。GlobAlbedo产品的反演模型为半经验核驱动模型RTLSR，经过计算得到地表反照率[13-14]。GLASS产品采用直接估算方法，通过训练建立适用于全球不同区域、不同覆盖类型、不同季节的格网查找表，直接估算得到地表反照率[15-17]。

然而，这些反照率产品应用于积雪研究仍然存在问题。MCD43A3、GlobAlbedo和GLASS这三种地表反照率产品主要是基于植被的辐射传输特性生产的，不能充分反映积雪的辐射传输特性。专门针对积雪的全球反照率产品仅有SAD，但它在反演过程中未考虑雪粒径和污化物浓度参数对积雪各向异性特性的影响，且缺乏白空反照率数据。同时，由于积雪反照率存在快速变化的特点，充分捕获积雪反照率的时间变化特征需要高时间分辨率的积雪反照率产品。当前的大部分产品采用了八天或16天合成的方法，如此低的时间分辨率无法表征积雪积累、消融过程中反照率快速变化的特征，导致这些反照率产品在积雪区域的精度欠佳[12,18]。而且，由于全球反照率产品验证数据主要分布在美洲、欧洲，在中国区域数量较少，因此这些产品在中国区域缺少系统性评价，产品选择和使用的不确定性较大。

本文基于积雪的辐射传输模型发展了适应中国区域特点的反照率反演算法，以 MODIS地表反射率产品作为输入，生产了2000-2020年中国积雪反照率的逐日和八天合成去云的遥感产品。利用地面站点和样方的积雪反照率观测数据，对产品开展了精度验证分析，结果表明产品的均方根误差、相关系数、平均偏差均优于国际同类产品。产品可为我国积雪资源分布调查提供数据支持，为积雪类型划分提供时空参数，为地表辐射研究和陆面过程模拟提供科学支撑，满足我国辐射、气候、环境、生态等相关研究的应用需求。

1 数据采集和处理方法

1.1 数据源

本文主要基于MODIS地表反射率反演积雪反照率。因此，数据源包括MOD09GA地表反射率产品、积雪范围、MCD12Q1土地覆盖产品、SRTM数字高程数据、林区积雪反照率观测数据共5类。

（1）MOD09GA地表反射率产品是对MODIS波段1到波段7的辐亮度数据进行大气校正处理，得到的500 m分辨率的反射率数据集和1 km分辨率的观测几何数据集[19]。本文使用2000-2020年逐日MOD09GA数据，最新版本为Version 6，它改进了气溶胶反演算法，构建新的气溶胶反演查找表，同时对数据内部雪、云以及云影等检测算法进行了改进。

（2）积雪范围来自中国MODIS逐日无云500 m积雪面积产品[20]。它基于MOD09GA地表反射率产品，结合积雪野外调查数据、数字高程数据、土地覆盖数据，通过决策树分类方法生成了中国2000-2020年积雪面积产品，与MODIS积雪产品MOD10A1/MYD10A1相比具有更高的精度。

（3）MCD12Q1土地覆盖产品提供每年的全球土地覆盖类型[21]。它通过对MODIS地表反射率数据进行监督分类得到初始分类结果，利用先验知识和辅助信息进行后处理，得到包含六种不同分类标准的土地覆盖数据。本文使用最新Version 6版本的MCD12Q1数据，因目前仅有2001-2019年的逐年产品，故2000年的土地覆盖用2001年代替，2020年的土地覆盖用2019年代替。选择其中的LC_Type1分类标准，提取森林范围。

（4）数字高程数据采用2000年的SRTM数据[22]。SRTM雷达系统搭载于奋进号航天飞机，采用干涉合成孔径雷达技术生成了全球范围内的数字高程模型。本文使用经过数据空白填补的最新Version 4版本，分辨率为90 m。

（5）林区积雪反照率数据来自长白山和小兴安岭森林生态监测站的塔台观测。长白山森林生态监测站位于吉林省安图县二道白河镇，位置为128.47 E、42.40 N，海拔736 m。小兴安岭森林生态监测站位于伊春市五营区，位置为129.20 E、48.23 N，海拔549 m。两个站都架设有森林观测塔，观测塔下方为长势均一的次生林，与周围情况一致，具有较好的空间代表性。在观测塔30 m高处架设四分量辐射表，自动获取长时间序列的林上反照率连续观测数据，时间范围分别为 2003-2010年2017-2019年。

1.2 数据处理方法和流程

首先基于积雪范围产品对反射率进行掩膜，获得积雪反射率；对积雪反射率进行地形校正、各向异性校正，获得窄波段反照率；经过窄宽转换，得到宽波段积雪反照率；基于土地覆盖数据获得林区范围，对林区积雪反照率进行校正；经过异常值处理，获得逐日积雪反照率产品；进行八天均值合成，获得去云积雪反照率产品。因此，数据处理方法和流程包括预处理、反照率计算、林区校正、异常值处理、任务运行、后处理共六个部分（图1）。

（1）预处理。基于MODIS地表反射率数据获得积雪覆盖范围和云覆盖区域，然后对反射率进行掩膜，获得积雪反射率。将数字高程数据和其他数据的投影转换至MOD09GA的坐标系统，并通过重采样统一空间分辨率。利用数字高程数据，采用C校正模型对积雪反射率进行地形校正，降低山区起伏地形对反照率反演的影响。

图1 中国积雪反照率遥感产品生产流程Figure 1 Production process of remote sensing product of snow albedo over China

（2）反照率计算。采用基于积雪弱吸收特性建立的渐近辐射传输（Asymptotic Radiative Transfer，ART）模型[4]，它可以较好地模拟积雪的双向反射特性[23]。因此，使用ART模型对积雪反射率进行各向异性校正，获得窄波段反照率。计算积雪的窄波段黑空和白空反照率的公式分别为：

式中，rb(θ0)和rw分别为积雪黑空和白空反照率，θ、θ0、φ分别为观测天顶角、太阳天顶角、相对方位角，R(θ,θ0,φ)为积雪反射率，f(θ,θ0,φ)为积雪散射特征函数，K0(θ0)为半无限-无吸收介质中光子溢出归一化角分布函数，R0(θ,θ0,φ)为半无限-无吸收介质的反射函数，计算公式分别如下：

式中，P(Ω)和Ω的计算公式分别为：

最后，采用针对积雪高反照率特性得出的窄宽转换公式[24]，获得宽波段积雪反照率：

式中，α为宽波段反照率，α1到α7分别为波段1到波段7的窄波段反照率。

（3）林区校正。林区积雪存在冠层截留、林下堆积、多次散射等复杂情况，反照率比纯雪低得多。由于缺少实用的森林积雪辐射传输模型，利用线性拟合方式进行林区积雪反射率校正。利用土地覆盖数据获得林区范围，得到下垫面类型为森林的积雪区域。基于森林生态监测站的反照率观测数据，通过格网搜索得到林区积雪反照率校正系数，利用该校正系数对所有的林区积雪反照率像元进行校正。

（4）异常值处理。当反射率产品中的云雪标识错误时，误判为雪的云会造成该像元反照率值显著高于正确识别为雪的像元反照率值。因此，通过格网搜索确定阈值，检测并剔除积雪反照率中异常高值的像元，得到正确的积雪反照率像元。

（5）任务运行。产品生产采取云端（Google Earth Engine）和本地端（Python）交互生产的方式，以最大限度地提高生产效率。在本地端编写脚本，提交生产任务至云端；在云端读取源数据，生产积雪反照率产品，下载至本地进行后处理。

（6）后处理。通过八天均值合成的算法，获得去云的积雪反照率产品。最后将产品由tiff格式转换为hdf格式，以提高存储和使用的灵活性。并生成png格式的产品缩略图，以供预览。

2 数据样本描述

基于上述数据处理方法和流程，生产得到的 2000-2020年中国积雪反照率遥感产品，命名为ChinaSA（China’s Snow Albedo）。数据示例如图2，数据波段信息见表1。

产品的数据命名规则如下：PLT_MODIS_SAB_YYYYMMDD_RESOL_XXXXM_VXX.HDF，其中 PLT表示卫星平台，MODIS表示传感器，SAB表示积雪反照率，YYYYMMDD表示年月日，RESOL表示时间分辨率，XXXXM表示以米为单位的空间分辨率，VXX表示版本号。

图1 中国2020年1月1日积雪反照率遥感产品示意图（审图号：GS(2022)747号）Figure 2 Remote sensing product of snow albedo over China on January 1, 2020

表1 中国积雪反照率遥感产品波段信息Table 1 Band information of remote sensing product of snow albedo over China

3 数据质量控制和评估

3.1 验证方法

利用多种地面站点的积雪反照率观测数据，结合积雪反照率样方观测数据，与遥感反演的积雪反照率产品进行对比，计算均方根误差（RMSE）、相关系数（R）、平均偏差（AD），验证积雪反照率产品的精度。

收集的地面站点观测数据包括：中国气象辐射基本要素逐小时曝幅量数据集、全球协调加强观测计划（CEOP）亚澳季风青藏高原试验研究（CAMP-Tibet）数据集、中国陆地生态系统通量观测研究联盟（ChinaFLUX）数据集、中国典型积雪区积雪站积雪特性观测数据集、森林生态监测站数据集。利用高分辨率遥感影像对站点逐一进行了空间代表性分析，发现拥有反射辐射数据的19个气象站全部位于城区，观测场的实测值高估了混合像元的反照率，因此主要选用另外三种数据进行验证。CAMP-Tibet数据集包含2002-2004年7个站点的观测数据，ChinaFLUX数据集包含2003-2010年10个站点的观测数据，积雪站数据集包含2017-2020年6个站点的观测数据。森林生态监测站数据集包含2003-2019年2个站点的观测数据，随机抽取一半数据做校正，另一半数据做验证。站点数据的筛选规则为：（1）选取北京时间上午10:30左右的地面实测数据，以保证地面实测时间与卫星过境时间一致；（2）地面实测与遥感反演的反照率值皆大于0且小于1，以去除无效数据；（3）站点所在像元当天的积雪覆盖率大于95%，保证观测场的实测值基本能代表整个像元。经过以上筛选，共得到242个站点验证数据。

站点对像元的验证仍然存在尺度效应问题，进一步通过样方观测提高验证的效度。针对中国主要积雪区，根据MODIS像元对应的地面范围，结合高分辨率影像提前确定样方范围，一般选择下垫面较为平坦的农田、戈壁、草原等开阔区域。设计便携式反照率观测装置，分纯雪像元样方和混合像元样方两种情况进行观测。纯雪像元样方的反照率为中心点的实测值，混合像元样方的反照率为不同地物的实测值与无人机获取的面积比例进行加权求和的结果。测得的样方数据包括：2019年1月北疆地区67个样方、2020年1月东北-内蒙古地区73个样方。样方数据的筛选规则为：样方所在像元当天被产品标识为积雪。利用该规则进行筛选，共得到82个样方验证数据。

筛选后的验证数据空间分布如图3所示，覆盖了东北-内蒙古（A）、青藏高原（B）、北疆（C）三大典型积雪区。图中的圆点表示站点验证数据，方形表示样方验证数据，数字表示各站点、样方用于验证的数据记录条数。

图3 中国积雪反照率遥感产品验证数据空间分布（审图号：GS(2022)747号）Figure 3 Spatial distribution of validation data for remote sensing product of snow albedo over China

3.2 验证结果与精度比较

利用324个站点与样方实测数据验证积雪反照率产品ChinaSA的精度，并与美国宇航局16天合成地表反照率产品 MCD43A3、逐日积雪反照率产品 SAD、欧空局八天合成地表反照率产品GlobAlbedo、北京师范大学八天合成地表反照率产品GLASS的精度进行对比，结果如图4所示。图中左下方为林区实测数据验证结果，右侧为非林区实测数据的验证结果，可见ChinaSA产品在非林区和林区都与实测数据具有更好的一致性。进一步计算五种产品的精度指标，如表2所示。ChinaSA产品与实测值的均方根误差为0.114，相关系数为0.822，平均偏差为0.087，优于国际同类产品。

图4 五种反照率产品的积雪反照率与实测值对比图Figure 4 Comparison of snow albedo of the five products and the in-situ data

表2 五种反照率产品的积雪反照率精度Table 2 Snow albedo accuracies of the five products

积雪反照率精度的提升，主要源于以下三方面的贡献：（1）ChinaSA产品采用专门针对积雪弱吸收特性建立的ART模型，能够更好地表达积雪的各向异性反射特性；其他地表反照率产品主要采用针对植被建立的辐射传输模型，在积雪区域的精度较低；积雪反照率产品SAD需要较多的输入参数，限制了模型的适用性；（2）利用我国森林生态监测站获取校正系数，显著提升了森林区域的积雪反照率反演精度；（3）通过地形校正处理，进一步提升了起伏地表的积雪反照率反演精度。

4 数据价值

本数据集为中国2000-2020年积雪反照率遥感产品，可为我国积雪资源分布调查提供数据支持，满足我国辐射、气候、环境、生态等相关研究的应用需求。

5 数据使用方法和建议

数据可通过国家冰川冻土沙漠科学数据中心（www.ncdc.ac.cn）下载。数据下载后可用ArcGIS/ENVI等软件或者相应程序代码读取，数据中的波段与本文表1中的波段相对应。

致 谢

感谢南京大学冯学智教授、张学良副教授和研究生叶李灶、马威、刘豪、盛光伟、马腾耀、宋依娜、李振世等参与北疆、青藏高原、东北-内蒙古等地的积雪反照率野外观测验证以及首席科学家王建研究员等专家对产品生产的指导、支持与帮助。