高亚洲地区被动微波遥感雪水当量数据集

邱玉宝，卢洁羽 ，石利娟，解鹏飞，梁雯珊，王星星

1. 中国科学院遥感与数字地球研究所，数字地球重点实验室，北京 100094

2. 湖南科技大学，湖南湘潭 411201

数据库（集）基本信息简介

空间分辨率时间分辨率数据量数据格式（1）N A S A-H M A-S W E产品：0.2 5°等经纬度格点；（2）R A D I-Q T P-S W E产品：0.2 5°等经纬度格点；（3）B l e n d e d-H M A-S W E产品：5 0 0 m。（1）N A S A-H M A-S W E产品：每日、每五日、每月（2）R A D I-Q T P-S W E产品：每日、每五日、每月（3）B l e n d e d-H M A-S W E产品：每日2 0.2 9 G B（1）N A S A-H M A-S W E产品数据格式：*.b i n, *.p n g（2）R A D I-Q T P-S W E产品数据格式：*.b i n, *.t i f, *.h d f, *.p n g, *.q a, *.m e t（3）B l e n d e d-H M A-S W E产品数据格式：*.b l d, *.j p g数据服务系统网址 h t t p://w w w.s c i e n c e d b.c n/d a t a S e t/h a n d l e/6 6 0基金项目数据库（集）组成公益性行业（气象）科研专项（G Y H Y 2 0 1 2 0 6 0 4 0），国家重点研发计划“政府间国际科技创新合作”重点专项（2 0 1 7 Y F E 0 1 1 1 7 0 0），中国科学院国际合作局对外合作重点项目（1 3 1 2 1 1 K Y S B 2 0 1 7 0 0 4 1），国家自然科学基金（4 1 3 7 1 3 5 1）。数据集由3个压缩包组成，文件以产品名及其覆盖时间范围命名：（1）N A S A-H M A-S W E_2 0 0 2-2 0 1 1.z i p是N A S A-H M A-S W E产品数据，数据量约1 8 9 M B；（2）R A D I-Q T P-S W E_2 0 0 3-2 0 1 1.z i p是R A D I-Q T P-S W E产品数据，数据量约1 1 6 M B；（3）B l e n d e d-H M A-S W E_2 0 0 3-2 0 1 1.z i p是B l e n d e d-H M A-S W E产品数据，数据量约1 9.9 9 G B。

引 言

雪水当量是全球气候变化的敏感要素之一，是陆表水文、能量模型及气候模式中的重要参数[1]。在常规地面气象台站观测中通常是直接测量雪深，要获得雪水当量值，则需要通过实测的积雪深度和积雪密度求得。遥感弥补了高亚洲地区传统积雪实测方法中站点少的不足，可以提供持续长时间、大范围积雪监测数据。光学传感器可进行积雪面积、积雪反射率等的监测，但是不能有效地估算积雪深度、雪水当量等参数[2]，而被动微波遥感在积雪雪水当量的估算中可以实现大范围高时间分辨率的监测和估算，是对可见光遥感监测的重要补充。

高亚洲地区北半球中纬度海拔最高、积雪覆盖最大的地区，素有雪域高原之称，特别是青藏高原是我国的大江大河的发源地，冰雪融水是这些河流的主要补给源。由于其地处中纬度、平均海拔高、面积大（占我国总面积的四分之一）、地形地貌复杂，使得该地区具有非常独特的气候单元，影响着其周边区域的气候。高精度的高亚洲陆地季节性积雪信息的获取对于气候变化的评价和预测具有重要的意义。

目前世界上已有很多科研和气象业务部门，均开展积雪反演算法研究，具有业务化数据产品提供能力，如在美国NOAA IMS雪冰服务、雪冰数据中心（NSIDC）[3-4]，欧空局ESA GlobSnow-2[5]、GMES（全球安全和环境观测计划）雪冰服务、EUMESAT（欧洲气象中心）的 H-SAF等，国内科学家也发布了几套专门针对我国的积雪雪水当量产品[6-9]，如西部数据中心发布的一套中国范围的SSM/I中国雪深长时间序列数据集[10]等。目前版本相对比较老，最近几年的数据产品未有更新。国际积雪产品在我国的山区和西部地区因缺少地面数据而疏于算法的校正开发，如NASA的积雪雪水当量产品在高亚洲地区高估严重，而GlobSnow产品在高亚洲地区不做估算等问题。

由于高亚洲地区的地形复杂，积雪覆盖分布不均且日变化快[11]，雪水当量算法的研发存在着缺少地面观测和专有模型，数据产品的大部分为全球产品的子集的缺陷。本数据集采用青藏高原地区不同下垫面算法重新计算了青藏高原地区的雪水当量产品（为RADI算法）。此外，还综合考虑了光学和微波遥感的各自优势，利用MODIS数据判别积雪，综合AMSR-E的SWE数据产品，形成了一套高时空分辨率的积雪综合产品数据集。

本数据集一共包括3种雪水当量产品：一是NASA AMSR-E的雪水当量（SWE）在高亚洲地区的重新计算产品，即NASA-HMA-SWE产品；二是青藏高原地区RADI算法的雪水当量（SWE）产品，即 RADI-QTP-SWE产品；三是高亚洲地区光学和被动微波综合雪水当量产品，即Blended-HMA-SWE产品。

1 数据采集和处理方法

1.1 输入数据

NASA-HMA-SWE产品和Blended-HMA-SWE产品的范围为高亚洲地区，高亚洲是亚洲中部以青藏高原为中心的高海拔区域，地理范围为25°N-45°N，67°E-107°E（图1）。RADI-QTP-SWE产品范围为青藏高原，地理范围为 26°N-40°N，73.45°E-104.6°E。

图1 高亚洲研究区域范围

1.1.1 NASA-HMA-SWE和RADI-QTP-SWE产品输入数据

NASA算法和RADI算法都用了AMSR-E L2A亮温数据，AMSR-E L2A亮温数据来自美国雪冰数据中心NSIDC（http://nsidc.org/），使用的亮温数据包含10.65 GHz、18.7 GHz、36.5 GHz和89.0 GHz 4个频率的水平和垂直极化共8个通道的轨道数据。

NASA算法使用的辅助数据包括：森林覆盖度数据（The 1/120th degree 〔1 km〕 MODIS/Terra Land Cover Type 96-Day L3 Global 〔MOD12Q1〕）、森林密度数据（The 1/120th degree 500 m MODIS Vegetation Continuous Fields）、积雪判识数据（基于MOD12Q1的The possibility of snow mask and the land/ocean/ice mask数据文件）、地面高程数据（GTOPO30）和积雪密度分布数据。

RADI算法使用的辅助数据包括：地面实测数据、四种地物覆盖度权系数辅助数据、水体覆盖度，地物标记和提取地物覆盖度数据、积雪判识数据（基于MOD12Q1的The possibility of snow mask and the land/ocean/ice mask数据文件）和积雪密度分布数据。

其中RADI算法考虑到青藏高原具体地表特征及国内外主流的土地覆盖划分类别，使用每年的MODIS土地覆盖类型数据（IGBP，其空间分辨率为500 m），将青藏高原下垫面类型进行重新分类，分为水体、森林、灌木、草原、裸地和其他共6类。重分类后的青藏高原下垫面分类具有明显的过渡性（图2）。高原的东南边缘部分为森林覆盖区，其西北方向为大片草原区；高原的西北部为大片裸地；在裸地与草原交界的南部为灌木区；高原上零星分布有水体。

图2 青藏高原下垫面重分类图[12]

1.1.2 Blended-HMA-SWE输入数据

采用已发布的2002-2018年青藏高原MODIS 逐日无云积雪面积数据集[13]和NASA-HMA-SWE产品作为输入。

其中MODIS 逐日无云积雪面积产品是通过多个云判别步骤得到的。对NASA的AMSR-E SWE产品进行重采样，使其空间分辨率和 MODIS逐日无云积雪面积产品一致。综合这两种数据，以此获得既能反映雪水当量数值，又能表现较高空间分辨率的积雪面积的图像。

1.2 数据处理算法和计算系统

下面分别介绍3种数据产品的处理流程及批量生产3种数据产品的系统。

1.2.1 NASA-HMA-SWE和RADI-QTP-SWE产品

（1）NASA-HMA-SWE数据处理

NASA算法的具体处理步骤分为以下几步：

第一步，利用AMSR-E/Aqua L2A全球轨道重采样亮温数据，以及辅助数据集，逐点对每一个采样点进行判识，对判识为可能是积雪的采样点进行计算，得到每一个采样点的雪深（SD）值。若采样点被判识为浅雪，则SD的值为5.0 cm；若为中到厚雪，则利用如下计算公式：

其中，SDf为森林覆盖的采样点的雪深值，单位为cm；SDo为无森林覆盖的采样点的雪深值，单位为cm；ff为森林覆盖度。

第二步，将各采样点的SD值投影到25 km EASE-Grid网格单元中。投影时，对多个采样点落入统一网格单元的SD值进行累加，记录累加次数，对网格单元中的SD值求平均值。

第三步，利用求得的SD均值及雪密度文件，计算每一个25 km EASE-Grid网格单元中的SWE值，计算公式如下所示：

最后得到空间分辨率为25 km的NASA算法的积雪雪水当量产品。同时对得到的SWE产品重采样，使其空间分辨率与MODIS一致，统一为500 m，便于后续进行Blended积雪产品的生产。

（2）RADI-QTP-SWE数据处理

RADI算法是以AMSR-E L2A亮温数据为基础，结合相关辅助数据，基于不同下垫面，利用青藏高原改进的雪深、雪水当量反演算法得到青藏高原的日雪水当量产品，五日雪水当量产品和月雪水当量产品。

图3 RADI雪水当量算法流程图

算法将利用辅助数据中的积雪判识数据先对像元进行积雪初步判别，剔除含干扰因子的部分，即水体、不可能为积雪和冰川的像元。仅对判识为可能有积雪覆盖的像元进行下一步计算。计算流程如图3所示，计算算法见表1。

表1 青藏高原各下垫面改进雪水当量算法说明[12]

表 1中，SWE为像元内总的雪水当量，单位为 mm；Ff、Fs、Fg、Fb分别代表森林、灌木、草原、裸地覆盖类型在像元内所占的百分比；SDf、SDs、SDg、SDb分别表示森林、灌木、草原、裸地覆盖类型的雪深，单位为cm；为积雪密度，单位为g/cm3；t1036v=TB10V-TB36V，即垂直极化的10 GHz与36 GHz的差值，也称为亮温梯度；t1937v=TB19V-TB37V，即垂直极化的19 GHz与37 GHz的差值；t1963v=TB19V-TB63V，即垂直极化的19 GHz与63 GHz的差值；pol36、pol89为动态因子，公式分别为pol36=t36v-t36h；pol89=t89v-t89h，分别表示36 GHz和89 GHz的极化差；p1为动态系数，如式（3）。

其中tb36v为36 GHz垂直极化的亮度温度，tb36h为36 GHz水平极化的亮度温度。

（3）青藏高原雪水当量被动微波遥感估算系统

NASA-HMA-SWE产品和 RADI-QTP-SWE产品均采用“青藏高原雪水当量被动微波遥感估算系统”开展计算[14-15]，该系统是基于 Linux平台的被动微波积雪雪水当量生产系统，实现了利用AMSR-E L2A亮温数据批量高效生产每日、五日、月雪水当量产品。该系统不仅可实现积雪制图的批量生产，而且可实现新算法产品的快速更新，结合计算服务器可具备近实时业务化处理能力，完成长时间序列的雪水当量制图工作，便于用户预览和产品的快速发布以及后续的模型和各种应用，该系统生产的产品可较好地满足高亚洲积雪监测的应用需求。

1.2.2 Blended-HMA-SWE产品生产

Blended-HMA-SWE产品采用 NASA-HMA-SWE（雪水当量）数据和 MODIS 逐日无云积雪面积数据进行融合处理，以此获得既能反映雪水当量数值，又能表现较为精确的积雪覆盖面的图像。计算流程如图4所示。

图4 Blended-HMA-SWE产品生产过程

具体步骤如下：第一步，统一两种数据的分辨率，对AMSR-E SWE数据重采样，使其空间分辨率与MODIS一致，都为500 m。第二步，提取AMSR-E SWE数据数值为0-240的数据，赋值为1，其余赋值为0。第三步，提取MODIS SCA数据数值，若像素值为200（雪），则重新赋值为1，其余赋值为4。第四步，用赋值之后的两个影像数据进行减法运算，SWE-SCA，各数字代表含义见表2，计算结果：SWE(1)-SCA(1)=0；SWE(1)-SCA(0)=1；SWE(-1)-SCA(1)=-2; SWE(-1)-SCA(0)=-1。运算结果数字代表含义如表2所示。第五步，利用步骤4中得到的图像，在SWE数据（取0-240为有效值）上建立掩膜文件，最后得出的图像既有较精确的积雪覆盖区域面积，又有雪水当量数值数据。

表2 运算结果数字代表含义

2 数据样本描述

2.1 产品命名规则

2.1.1 NASA-HMA-SWE和RADI-QTP-SWE产品命名规则

NASA生产的高亚洲地区雪水当量产品时间覆盖范围为AMSR_E L2A亮温数据全生命周期，时间从2002年到2011年，分为日产品、五日产品、月产品3种。RADI生产的青藏高原地区雪水当量产品，时间从2003年到2011年，分为日产品、五日产品、月产品3种。NASA-HMA-SWE产品和RADI-QTP-SWE产品命名规范表如表3所示。

NASA-HMA-SWE产品命名约定：

（1）日产品：HMA _RADI_AMSR_E_Aqu_01SWE01_yyyymmdd.bin；

（2）五日产品：HMA _RADI_AMSR_E_Aqu_05SWE01_yyyymmdd.bin；

（3）月产品：HMA _RADI_AMSR_E_Aqu_MOSWE01_yyyymm.bin。

RADI-QTP-SWE产品命名约定：

（1）日产品：QTP_RADI_AMSR_E_Aqu_01SWE02_yyyymmdd.bin；

（2）五日产品：QTP_RADI_AMSR_E_Aqu_05SWE02_yyyymmdd.bin；

（3）月产品：QTP_RADI_AMSR_E_Aqu_MOSWE02_yyyymm.bin。

表3 NASA-HMA-SWE产品和RADI-QTP-SWE产品命名规范说明

说明：日产品命名时，由处理亮温数据的日期决定日产品名称中的日期信息；五日产品命名时，由所处理的日产品中日期最靠前的日产品的日期决定五日产品名称中的日期信息；月产品命名时，由所处理的日产品所在月份决定月产品名称中的日期信息。

2.1.2 Blended-HMA-SWE产品命名规则

经处理得到每日微波和光学Blended积雪产品，产品命名为：Blend_HMA_yyyydoy_MOD-AMS R.bld，具体命名规则见表4。

表4 Blended产品命名规范说明

变量y y y y d o y描述年d a y o f y e a r简称M O D M O D I S去云积雪面积产品A M S R A M S R_E雪水当量产品

2.2 分类描述（数值描述）

计算出来的雪水当量数据集中像素点的可能取值有6种，其对应的描述见表5。

表5 雪水当量数据集中数据取值分布

2.3 部分数据样本示例

使用上述青藏高原雪水当量被动微波遥感估算系统批量自动生产出NASA-SWE和RADI-SWE日产品、五日产品、月产品。

经过 Blend算法得到 Blended-HMA-SWE积雪产品，示例如图 5所示。灰色区域表示NASA-HMA-SWE产品数据范围，蓝色区域表示MODIS 逐日无云积雪面积数据范围，白色部分表示没有积雪覆盖的地区。

图5 2010年1月1日至1月3日Blended雪水当量每日产品示例图

3 数据质量控制和评估

3.1 RADI-QTP-SWE产品与地面实测数据对比验证

所用的地面实测数据来源于气象站台实测记录，地面实测数据记录中的记录内容包括省名、观测站点号、经纬度及实测雪深（单位：cm）。87个地面气象站积雪观测数据将用来验证青藏高原积雪改进算法的精度。如图1所示，黑色圆点为RADI-QTP-SWE产品地面气象台站验证点位置。

对2009年12月至2010年3月共计4个月的地面实测雪深数据进行统计，生成实测数据直方图（图6）。从图中可以明显看到，研究区的站台记录中，雪深基本小于10 cm。研究区为典型的浅雪覆盖区，近二分之一区域的雪深仅为1 cm。

图6 实测数据直方图[12]

利用青藏高原地区4种下垫面（森林、灌木、草地、裸地），对应的4个地面实测站点（聂拉木、嘉梨、玛多、都兰）测得的雪深观测数据，验证各下垫面覆盖类型上改进的雪深算法的精度。验证结果见图7及表6。

图7 青藏高原改进雪深算法在不同下垫面验证散点图[12]

由上表结合已发表的文献[16]可看出改进的算法结果明显优于相同时间相同研究区域选取的其他5种算法，算法精度有一定的提高。

表6 青藏高原改进的雪深算法验证结果[12]

3.2 NASA-HMA-SWE产品与RADI-QTP-SWE产品对比验证

将使用由改进的青藏高原雪水当量算法生成的RADI-QTP-SWE产品与NASA AMSR-E SWE产品进行对比验证，利用地面数据验证两种产品对应网格点的精度误差，并进行统计。精度达标的衡量标准分为两种情况：1）当雪水当量小于等于10 mm时，单点估算误差小于等于4 mm即为精度达标；2）当雪水当量大于10 mm时，单点相对误差小于等于20%即为精度达标[12]。利用地面数据验证由改进的青藏高原雪水当量产品与NASA AMSR-E SWE产品对应网格，产品的单点精度达标率统计结果见表7。

表7 雪水当量数据精度达标率统计表[12]

从表7可以明显看到，改进的算法得到的雪水当量产品的精度达标率优于NASA官方AMSR-E雪水当量产品。

图8 RMSE统计直方图[12]

验证得到的RMSE统计直方图如图8所示，改进算法得到的雪水当量产品的RMSE低于NASA官方AMSR-E雪水当量产品50%以上。

4 数据使用方法和建议

本数据集共包括3种产品：NASA官方AMSR-E雪水当量产品（2002-2011年）；利用青藏高原改进算法得到的雪水当量产品（2003-2011年）；及由AMSR-E的雪水当量数据与MODIS的积雪覆盖数据融合的Blended积雪产品（2002-2011年）。本数据集适用于高亚洲地区水和能量循环、水资源、生态和灾害等相关问题研究，特别是对于高亚洲地区冰雪模型、高亚洲地区气候变化研究、积雪长时间时空分布规律、河流流量变化分析、生态效益、雪灾预测以及未来变化趋势等方面等均具有重要价值，并可为进一步针对被动微波雪水当量的降尺度算法提供重要的参考。