中国区域Landsat 地表反射率产品

彭燕，何国金，张兆明，尹然宇

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

2. 中国科学院大学，北京 100049

3. 海南省地球观测重点实验室，海南三亚 572029

引 言

Landsat 系列卫星自1972 年发射以来，已获取了大量中空间分辨率的卫星影像，这些影像记录着人类活动和自然变化，成为最长时间序列的星载陆地观测数据集[1]。地表反射率常用于土地覆盖及变化研究，是很多地表地球物理参数（如叶面积指数、叶绿素和生物量）反演的关键输入参数。同时，随着对地观测领域进入大数据时代，卫星遥感数据应用用户和研究者不仅希望遥感数据在几何位置上具有一致性，同时也对遥感数据的辐射一致性提出了更高的要求，以更好地进行遥感数据应用和信息挖掘分析。由于卫星遥感数据的大气校正是一项繁琐而专业的工作，因此，研究生产高质量中国区域的Landsat 地表反射率产品并公开共享，具有重要的现实意义。

本文提供30 m 空间分辨率的中国区域Landsat 系列的地表反射率产品，给出了中国区域Landsat数据地表反射率产品生产的技术流程：首先采集覆盖中国区域的Landsat 正射影像，然后经过辐射定标，获取6S（Second Simulation of Satellite Signal in the Solar Spectrum）辐射传输模型所需要的几何参数和大气参数并进行大气校正，生成中国区域Landsat 地表反射率产品。

1 数据采集和处理方法

1.1 数据采集方法

本文所生产的中国区域地表反射率产品主要是在Landsat 系列正射影像数据的基础上进行生产的，1980s-2012 年为Landsat 5 数据，2000-2003 年期间为无条带的Landsat 7 数据，2013-2019 年采用Landsat 8 数据，这些数据主要由中国遥感卫星地面站接收。对于太阳天顶角大于76°的Landsat 数据不进行大气校正处理。

1.2 数据处理方法

中国区域地表反射率产品基于6S 辐射传输模型进行大气校正获得。6S 模型是目前比较完善的大气校正模型之一，由Vermote 等人在5S 模型的基础上改进而来[2]，适用于0.25-4 μm 波长范围内电磁波的大气辐射传输的模拟。假设气溶胶类型为大陆型，并将地表视为平面朗伯体，则其地表反射率可根据表观反射率计算而得到。其公式为[3]：

式中，ρs为地表反射率；ρTOA为表观反射率，通过辐射定标可计算得到；Tg(OG)为O2、O3、CO2、NO2、CH4气体总的透过率，Tg(OG)=Tg(O2)×Tg(O3)×Tg(CO2)×Tg(NO2)×Tg(CH4)；ρR+A为瑞利和气溶胶的反射率；TR+A为瑞利和气溶胶透过率；Tg(H2O)为水汽的透过率；SR+A为瑞利和气溶胶球面反照率。而Tg(OG)、ρR+A、TR+A、Tg(H2O)、SR+A等大气相关系数均可通过调用6S 模型计算得到。而基于6S 模型进行大气校正的关键和难点在于大气参数的获取。本数据集是基于公式（1）的计算原理和6S 辐射传输模型[1,3-4]进行生产的，首先需要进行辐射定标，然后获取6S 模型所需要的大气参数、DEM 以及几何参数进行大气校正，其技术路线如图1 所示。

1.2.1 辐射定标

表观反射率与大气顶层进入卫星传感器的光谱辐射亮度、日地间距离、大气顶层的平均太阳光谱辐照度，以及太阳的天顶角有关。Landsat 5/7 的表观反射率计算公式为：

式中，ρ为表观反射率，d为日地距离，θz为太阳天顶角（单位为度，与元数据文件中给出的太阳高度角互为余角），ESUN（Mean solar exoatmospheric spectral irradiances）是大气层顶平均太阳光谱辐照度，Lλ为光谱辐射亮度，单位为w/(m2∙μm∙sr)，利用定标系数进行线性计算得到。

图1 地表反射率产品生产技术路线图

Landsat 8 的表观反射率计算公式为：

式中，Mρ和Aρ可从元数据文件中获取，分别为Mρ= REFLECTANCE_MULT_BAND_x，Aρ= REFLECTANCE_ADD_BAND_x，x 为波段号。

1.2.2 模型参数获取

6S 模型的输入参数主要包括太阳天顶角、方位角和观测天顶角、观测方位角等观测几何参数以及水汽、臭氧、气溶胶、气压等大气参数和海拔高度信息。在本文中，观测几何参数可根据元数据文件获取，海拔高度信息则是利用0.05°空间分辨率的GCM DEM（Global Climate Model Digital Elevation Model）[5]得到的（https://search.earthdata.nasa.gov/search?q=ASTER&ok=ASTER）。对于Landsat 5/7 数据而言，大气校正所需要的水汽数据以及气压数据均来自于NCEP（NOAA National Centers for Environmental Prediction）再分析数据（http://dss.ucar.edu/），空间分辨率为2.5°×2.5°，臭氧 数 据 来 自 TOMS （ Total Ozone Mapping-Spectrometer ） 数 据（https://ozoneaq.gsfc.nasa.gov/data/ozone/#），空间分辨率为1.25°×1.00°，气溶胶光学厚度则是利用暗目标法反演得到的[3,6]。而对于Landsat 8 数据而言，大气校正所需要的水汽数据、臭氧数据、气溶胶光学厚度（AOT）来自于0.05°空间分辨率的MODIS09CMA（MODIS Aerosol Optical Thickness Climate Modeling Grid）和MODIS09CMG（MODIS surface reflectance Climate Modeling Grid）数据[7]（https://earthdata.nasa.gov/earth-observation-data/near-real-time/download-nrt-data/modis-nrt）。

2 数据样本描述

本数据集是以景为单位存放文件夹，由每个波段的地表反射率产品、质量文件（Quality Assessment, QA）、元数据、缩略图组成。

（1）文件夹的命名规则为卫星-传感器-path-row-成像日期-LSR，如L5-TM-115-026-19840418-LSR。

（2）每个波段的地表反射率产品：Landsat 5/7 包括波段1，2，3，4，5，7（即分别为蓝、绿、红、近红、中红外1、中红外2），Landsat 8 包括波段1，2，3，4，5，6，7（即分别为深蓝、蓝、绿、红、近红、中红外1、中红外2）。命名规则为卫星-传感器-path-row-成像日期-LSR-BX（X 表示第几个波段）.TIF，如L5-TM-115-026-19840418-LSR-B1.TIF。影像的空间分辨率为30 m，投影坐标系是WGS84 UTM。为了降低存储空间，地表反射率结果由原本的0-1 范围内的浮点型均乘以10000变成16 位整型，背景填充值为-9999，并进行了“LZW”的无损压缩。

（3）质量文件（QA）：包括在原始数据的基础上生成的QA，为PIXEL-QA，主要是对填充值（Fill）、晴空（Clear）、云（Cloud）、云置信度（Cloud Confidence）、云阴影（Cloud Shadow）、冰雪（Snow/Ice）以及水（Water）等信息进行标识，命名规则为：卫星-传感器-path-row-成像日期-PIXELQA.TIF。和在地表反射率的基础上生成的QA，主要是对气溶胶相关的信息进行标识，Landsat 5/7 为LSR-CLOUD-QA，命名规则为卫星-传感器-path-row-成像日期-LSR-CLOUD-QA.TIF；Landsat 8 为LSR-AEROSOL，命名规则为卫星-传感器-path-row-成像日期-LSR-AEROSOL.TIF。Landsat 5/7 和Landsat 8 的QA 的属性分别如表1-4 所示。

表1 Landsat 5/7 的PIXEL-QA 属性表

表2 Landsat 5/7 的LSR-CLOUD-QA 的属性表

表3 Landsat 8 的PIXEL-QA 属性表

表4 Landsat 8 的LSR-AEROSOL 属性表

（4）元数据文件：包括原始数据的元数据文件（MTL.txt 文件），命名规则为卫星-传感器-pathrow-成像日期-MTL.txt；和地表反射率产品的XML 元数据文件，命名规则为卫星-传感器-path-row-成像日期-LSR.xml。XML 元数据文件描述了地表反射率产品的基本信息以及各个波段产品的相关信息。

（5）缩略图：包括512 像素大小和1024 像素大小的缩略图，命名规则分别为卫星-传感器-pathrow-成像日期-LSR-THUMB.JPG 和卫星-传感器-path-row-成像日期-LSR-BROWSER.JPG。

图2 为2005 年江西省的地表反射率结果示意图，其中图2（a）为大气校正前的影像，图2（b）为地表反射率结果，图2（c）为每景数据的成像时间，如0420 表示2005 年4 月20 日的Landsat5 数据。从图2 中可以看出，经过大气校正后能较大地消除大气的影响，具有较好的辐射一致性，尤其是2005 年6 月23 日这一景数据，其地表反射率结果明显消除了大部分薄云薄雾的影响。

图2 江西省地表反射率产品示意图

3 数据质量控制和评估

为了评估本数据集的可靠性，在全国范围内均匀选取了12 个轨道号如图3 所示，每个轨道号分别选取了Landsat5/7/8 地表反射率（具体的时间见图3 的列表），同时每景影像均随机选取50 个样本点，与USGS 发布的地表反射率产品进行对比验证。结果显示本数据与USGS 所发布的地表反射率产品的相关性系数R2均在0.99 以上，均方根误差（RMSE）均小于1%，这表明本数据集与USGS所发布的地表反射率产品具有很好的一致性，如图4 所示。Feng 等人将USGS 发布的全球2000 年和2005 年的Landsat 5/7 地表反射率产品分别与相应的MODIS 地表反射率产品进行比较分析，结果表明Landsat 5 TM 的均方根差（RMSD）约为2.2%-3.5%，Landsat 7 ETM+的均方根差约为1.3%-2.8%[8]。同时利用2014 年6 月11 日获取的南京地区的实测光谱数据对Landsat 8 数据地表反射率产品进行验证，结果显示其地表反射率产品的RMSD 约为2%-4%之间。并对USGS 发布的Landsat8地表反射率产品进行验证，结果显示USGS 发布的Landsat8地表反射率结果的RMSD分别为4.07%、4.37%、3.49%和5.3%，而本产品的RMSD 分别为2.71%、3.1%、2.71%、2.08%，本产品的RMSD均低于USGS 发布的地表反射率产品[1]。

4 数据价值

本文推出的1986-2019 年中国区域Landsat 系列卫星长时序地表反射率产品，采用目前国际上较为成熟的地表反射率算法，精度较高，未来将持续更新，在森林、水资源、气候变化等领域长时序信息挖掘分析方面具有重要的应用价值。

图3 对比验证的轨道号分布示意图

5 数据使用方法和建议

本数据集可通过地球大数据科学工程（ CASEarth ） Databank 在线服务网址（http://databank.casearth.cn）获取数据。用户注册成功并登录系统后，进入平台产品查询界面，产品类型选择“LSR”，然后在所需要的数据下点击下载，随后会弹出一个下载对话框，该框中列出了各波段的地表反射率结果、QA 文件和元数据文件，根据需要下载数据即可。用户可根据行政区、地图选择以及行列号等方式查询所需要的数据。目前共享的产品主要包括1980s-2012 年云量小于50%的Landsat 5、2000-2003 年云量小于20%的Landsat 7 以及2018-2019 年的全部Landsat 8 的地表反射率产品，后续作者将持续生产我国区域云量大于50%的Landsat 5、云量大于20%的Landsat 7 和2019年以后的Landsat 8 地表反射率产品，以提供更好的、持续的数据共享服务。若平台系统中暂时缺少用户所需的数据或者有与本数据相关的其他数据需求，可通过咨询本文作者进行申请。

图4 与USGS 地表反射率产品对比验证结果

致 谢

衷心感谢王桂周和龙腾飞在产品生产时遇到的存储与效率问题给予的建设性意见。

数据作者分工职责

彭燕（1988—），女，湖南省郴州市人，在读博士，工程师，研究方向为遥感图像智能处理。主要承担工作：算法集成程序编写，数据生产流程设计，论文撰写。

何国金（1968—），男，福建省龙岩市人，博士，研究员，研究方向为遥感数据智能处理与信息挖掘。主要承担工作：总体思路与方案设计，论文修改。

张兆明（1980—），男，河南省郑州市人，博士，副研究员，研究方向为遥感数据智能处理与信息提取。主要承担工作：技术指导，论文修改。

尹然宇（1996—），男，山东省临沂市人，在读博士，研究方向为遥感图像智能处理。主要承担工作：数据挑选与整合。