信息速览

数据集名称五十年同一套网格 — Landsat 陆地影像档案

数据类型栅格影像,多光谱,地表反射率,质量评估波段,云掩膜

规模全球陆地 / 无个体对象

接入方式https://landsat.gsfc.nasa.gov/

AI 就绪度

做遥感模型的人手边通常有两种食材。第一种是切好的小图：一个瓦片、一个类别标签、几万到几十万张，下载下来就能塞进 DataLoader，跑出来的数字也好看。第二种是连续观测：同一个地方每隔若干天被同一类传感器拍一次，拍了几十年，像元对得像元，年份对得上年份。前者的麻烦在别处——瓦片是从更大的影像上随机切下来的，同一块地的不同瓦片会同时落进训练集和测试集，模型记住的是"这一带的纹理"而不是"这一类地物"，成绩里有相当一部分是地理泄漏给的。后者的麻烦在于它大得没法下载，且一个标签也没有。

Landsat 属于第二种的全部，第一种的一点都没有。它是 NASA 与国家地质调查局（USGS）共同维持的陆地卫星计划，Landsat 陆地影像档案就是这条计划的产物，1972 年 7 月 23 日第一颗星入轨至今没有断过档；档案存量超过一千万景，每天新增一千五百景上下。它的独门之处可以压缩成一句话：五十年里，同一块地的像元始终落在同一套网格上。这句话听起来平淡，做起来极难——它要求九颗卫星的传感器互相标定、几何配准控制在十米级、大气校正用同一套算法重跑整条历史。USGS 用"Collections"（集合）这套分级管理制度来做这件事：一次全档案重处理定一个集合号，集合之间像素级可比，集合内部不再变。当前在用的一代叫 Collection 2，2020 年 12 月发布，Level-2 第一次把全球地表反射率与地表温度做成了常备库存。

它的代价同样直白。三十米的像元看不清一栋房子，十六天的重访抓不住一场洪水的涨落过程，而云会直接作废掉一部分观测；档案整体在 PB 量级，不是"下载下来训"的东西，而是"在云上算"的东西；最关键的是，它没有任何人工标注——分类标签得用户自己造。所以它不是"数据集"这个词的通常含义，而是一条需要你自己定义监督信号的观测记录。一屏速览见下表。

维度	内容
全名	Landsat 系列陆地卫星计划及其影像档案（现行产品体系为 Collection 2）
发布方	NASA 负责仪器、平台与发射校验；USGS 负责运营、接收、归档、产品生成与分发（EROS 中心）
首次发布	1972-07-23 Landsat 1 入轨；2008 年整档转为免费开放
最近改动	Collection 2 于 2020 年 12 月发布；Collection 1 于 2022-12-30 停止下载
档案规模	USGS EROS 档案逾一千万景；Landsat 8 与 9 合计每日约一千五百景入档
时间跨度	Level-1：1972 年至今（1–9 号星）；Level-2 全球：1982 年至今（4–9 号星）
空间分辨率	MSS 80 米（产品重采样到 60 米）；TM / ETM+ / OLI 30 米；全色 15 米
波段数	MSS 4 个；TM 7 个；ETM+ 8 个；OLI 9 个加 TIRS 2 个
重访周期	单星 16 天，8 号与 9 号联合 8 天；1–3 号星为 18 天
单景覆盖	约 170 km 南北乘 183 km 东西；Landsat 8 约 6300 乘 6000 个三十米像元
原始来源	WRS-1 / WRS-2 全球轨道网格，USGS EROS 与国际地面站网
标注情况	无人工标注；逐像元质量评估波段给出云、云影、雪、水、饱和与地形遮挡
文件格式	Cloud Optimized GeoTIFF，附 MTL / XML / JSON 元数据与 ANG 角度系数文件
许可	USGS 声明产品无使用限制；科学数据目录以 CC0-1.0 公共领域奉献标注
典型用途	土地覆盖与变化检测、植被指数时序、城市扩张、灾害前后对比、蒸散与地表温度
获取方式	EarthExplorer、M2M API、LandsatLook Viewer、AWS S3（us-west-2，请求方付费）、STAC

§1 它要解决的具体麻烦

1.1 两种食材的分工，以及各自的软肋

遥感算法论文里最常见的形态，是"一堆切好的瓦片加一个类别标签"。这个形态之所以流行，是因为它把问题简化成了普通图像分类：输入是三通道或十三通道的小图，输出是十类土地覆盖。它上手快、可复现、榜单漂亮，但它把两件真正困难的事藏起来了。

第一件是地理泄漏。瓦片是从有限的若干景原始影像上切下来的，相邻瓦片共享同一片天空、同一套大气条件、同一天的光照角度，甚至共享同一块地的同一片树冠。随机划分的时候，训练集和测试集会同时包含来自同一景的兄弟瓦片，模型只要学会"这一带的色调偏绿"就能拿分。这个问题在遥感圈被反复讨论，根源不在模型，而在切分口径——按地理块划分和按随机划分得到的数字能差出一大截。

第二件是时间冻结。瓦片数据集给出的是"某一时刻的地物类别"，它能回答"这是不是工业用地"，回答不了"这块地什么时候变的、变了几轮"。而现实里的监测需求——森林退化、城市扩张、矿区复垦、灾后重建——问的恰恰是后者。

Landsat 提供的是回答后者所需要的那件东西：同一块地在五十年里被反复测量的记录，且这些测量之间可以逐像元相减。它的软肋正好是瓦片数据集的长处：没有标签、体量巨大、需要云掩膜和大气校正的前置知识。两者不是替代关系，而是上下游关系——绝大多数遥感分类基准的原始影像，最终都能追到这一类星载档案。

1.2 三条常见的错误期待

第一条，期待它自带标签。它一个人工标注都没有。Level-2 产品给你的是物理量（地表反射率、地表温度）和逐像元质量标记，不是类别。想做监督学习，标签要自己造：要么配外部的土地覆盖产品，要么自己标注，要么用官方的 Level-3 产品（水体范围、火烧迹地、雪盖比例）当弱标签。这是一条硬边界，很多人是在下载完几百 GB 之后才发现的。

第二条，期待三十米能看清地物。三十米意味着一个像元边长相当于半个足球场，一栋独栋住宅只占几个像元，道路在多数情况下窄于一个像元而只表现为混合光谱。它适合看地块、看林分、看农田、看城市边界，不适合看车辆、看建筑轮廓、看单体树冠。需要米级或亚米级细节的任务应该去找商业或政府的高分影像，而不是在这里耗。

第三条，期待"免费开放"等于"零门槛"。数据本身确实不收费、无使用限制，但通道各有门槛：EarthExplorer 与 M2M API 需要注册 USGS 账号；AWS 上的 usgs-landsat 桶是请求方付费（requester pays）桶，需要 AWS 账号并承担请求与出网流量费用；Google Earth Engine 需要单独申请。费用不在数据上，在计算与传输上。

1.3 它能给到的与给不到的

你的需求	它给的	够不够
长时序土地覆盖变化检测	1982 年至今全球三十米地表反射率，像元级可比	够，这是它最不可替代的用法
植被指数时序与物候	红与近红波段齐备，十六天一个观测	够做年际与季节尺度，不够做日尺度
城市扩张与不透水面	三十米反射率加质量波段	够做市县级，不够做街区级
地表温度与热岛	1982 年至今地表温度产品，开尔文	够，注意它是地表温度而非气温
灾害前后对比	十六天重访，新数据六小时内可检索	够做事后评估，不够做应急实时
语义分割标注数据	无任何人工标注	不够，需自建或外部配准
目标检测（车、船、建筑）	三十米混合像元	不够，分辨率不足
小时级或日级过程监测	单星十六天，双星八天	不够，应看静止轨道或极轨高频传感器

1.4 你的处境与它能帮到的事

你的处境	它能帮你什么	下一步该做什么
你训好一个土地覆盖模型，怀疑成绩里有地理泄漏	按 WRS-2 路径行号或按地理块划分，做空间留一验证	先按路径行分组再划集合，别按随机像元划（§5.6）
你要做多年变化检测，手上的影像来自不同年份不同传感器	同一集合内跨传感器已做交叉辐射标定与几何配准	只用 Collection 2，且优先 Tier 1（§4.2）
你要算植被指数，担心大气把年际差异吃掉	Level-2 地表反射率已做大气校正	用 SR 波段而不是 Level-1 的表观反射率（§2.2）
你要做热岛或干旱监测	地表温度产品与多种中间波段	用 ST 波段反算开尔文，同时用 ST_QA 看不确定度（§5.5）
你要给模型造预训练语料	全球覆盖、量大、无版权负担	先定"以地理块为单位"的采样策略，再谈数据规模
你要复现已发表结果	处理日期写在产品名里，可精确锁定版本	记录完整产品标识，含采集日期与处理日期（§2.3）

1.5 一句话判断

它不是拿来直接训练的"数据集"，而是拿来造数据集的"底片库"：给你的是可比、可追溯、可全球铺开的物理量，不是标签。

1.6 不适合谁：四条边界

不适合想拿现成标签的人。 全档没有一条人工标注。所有监督信号都要自建或外接，这一步的工作量与不确定性往往超过模型本身。

不适合小团队全量落地。 档案在 PB 量级，几十景就是几十 GB。没有云上计算安排而打算"先下载下来看看"的团队，会在带宽和磁盘上先耗尽预算。

不适合需要高频观测的任务。 单星十六天、双星八天，汛期一周的过程、作物关键期的逐日胁迫、短时强对流前后，它都覆盖不到。这类需求应选更高时间分辨率的传感器。

不适合需要亚米级形态的任务。 三十米像元里混着屋顶、树冠与路面，任何依赖形状与纹理细节的目标级任务都超出了它的能力边界。

§2 数据构成与规格

2.1 一个产品包由哪些文件组成

Collection 2 的 Level-2 产品以"景"为单位分发，一景对应 WRS-2 网格上的一个路径行号与一个采集日期。同一景的所有文件共享一个产品标识前缀，靠后缀区分：

后缀	内容	说明
`_SR_B1` 至 `_SR_B7`	地表反射率波段	Landsat 8-9 为 1–7 波段；TM / ETM+ 为 1–5 与 7
`_ST_B10`	地表温度	TM / ETM+ 为 `_ST_B6`
`_QA_PIXEL`	逐像元质量	位打包，云、云影、雪、水、卷云及各项置信度
`_QA_RADSAT`	辐射饱和与地形遮挡	位打包，逐个波段一个饱和位
`_SR_QA_AEROSOL`	气溶胶质量	大气反演质量分级
`_ST_ATRAN` / `_ST_EMIS` / `_ST_EMSD`	大气透过率、地表比辐射率及其标准差	温度反演的中间量
`_ST_TRAD` / `_ST_URAD` / `_ST_DRAD`	热辐射、上行与下行辐射	温度反演的中间量
`_ST_QA` / `_ST_CDIST`	温度不确定度、到最近云的距离	单位分别为开尔文与千米
`_MTL.txt` / `_MTL.xml` / `*_MTL.json`	元数据	ODL、XML、JSON 三种并存
`_ANG.txt`	角度系数文件	逐像元太阳与观测角度
`*_thumb_large.jpeg` 等	快视图	仅供预览，不能用于分析

这个清单值得注意的一点是：温度产品附带了六七个中间量。它们不是冗余，而是让使用者能判断"这个温度值为什么是这样"——比如 ST_CDIST 告诉你这个像元离最近被标为云的像元有多远，云边若干千米内的温度反演是已知不可靠的，官方文献明确记载了这一点。

2.2 波段与取值的完整字典

下表把 Landsat 8-9 的 Level-2 波段逐条列出。需要强调：这里的"字段"不是表格列，而是栅格波段与它们的标定点。

波段文件	名称	数据类型	单位	填充值	乘系数	加常数	有效值区间
`_SR_B1`	1 波段反射率	16 位无符号整型	反射率	0	0.0000275	−0.2	0.0–1.0
`_SR_B2`	2 波段反射率	16 位无符号整型	反射率	0	0.0000275	−0.2	0.0–1.0
`_SR_B3`	3 波段反射率	16 位无符号整型	反射率	0	0.0000275	−0.2	0.0–1.0
`_SR_B4`	4 波段反射率	16 位无符号整型	反射率	0	0.0000275	−0.2	0.0–1.0
`_SR_B5`	5 波段反射率	16 位无符号整型	反射率	0	0.0000275	−0.2	0.0–1.0
`_SR_B6`	6 波段反射率	16 位无符号整型	反射率	0	0.0000275	−0.2	0.0–1.0
`_SR_B7`	7 波段反射率	16 位无符号整型	反射率	0	0.0000275	−0.2	0.0–1.0
`_ST_B10`	10 波段地表温度	16 位无符号整型	开尔文	0	0.00341802	149.0	150–359
`_ST_QA`	温度不确定度	16 位有符号整型	开尔文	−9999	0.01	—	0–327.67
`_ST_TRAD`	热辐射	16 位有符号整型	瓦每平方米球面度微米	−9999	0.001	—	0–22
`_ST_ATRAN`	大气透过率	16 位有符号整型	无量纲	−9999	0.0001	—	0–1
`_ST_EMIS`	地表比辐射率	16 位有符号整型	无量纲	−9999	0.0001	—	0–1
`_ST_CDIST`	到最近云的距离	16 位有符号整型	千米	−9999	0.01	—	0–240
`_QA_PIXEL`	逐像元质量	16 位无符号整型	位索引	1	—	—	位打包
`_SR_QA_AEROSOL`	气溶胶质量	8 位无符号整型	位索引	1	—	—	位打包
`_QA_RADSAT`	饱和与遮挡	16 位无符号整型	位索引	0	—	—	位打包

Collection 1 与 Collection 2 在这张表上的差别是最容易踩的坑：Collection 1 的填充值是 −9999、数据类型是有符号 16 位、反射率系数是 0.0001 且无加常数；Collection 2 全部改成了无符号整型、填充值 0、带加常数。官方解释这么改是为了与 Level-1 一致并更高效利用无符号整型的取值空间。混用两代的反射率图，会得到一张"看着正常但整体偏移"的镶嵌图。

2.3 一景的完整解剖

产品标识的格式是 LXSS_LLLL_PPPRRR_YYYYMMDD_yyyymmdd_CC_TX。以一个真实出现在官方教程里的标识为例：

LC08_L2SP_140041_20130503_20200911_02_T1

从左往右读：L 表示 Landsat；C 表示传感器为 OLI 与 TIRS 组合（O 为仅 OLI，T 为仅 TIRS 或 TM，E 为 ETM+，M 为 MSS）；08 表示 8 号星；L2SP 是处理级别，L2SP 表示反射率与温度都有，L2SR 表示只有反射率；140041 是 WRS-2 路径 140、行 41；20130503 是采集日期；20200911 是 Level-2 处理日期；02 是集合号；T1 是集合类别。

这个命名里藏着三条对使用者很关键的信息。第一，采集日期与处理日期是分开的——同一景可能被重处理，处理日期一变，产品标识就变，你手上的文件可能与别人论文里的不同。第二，路径行号直接写在名字里，这既是定位手段，也是划分数据集时最方便的分组键。第三，L2SP 与 L2SR 的区别意味着一部分景没有温度产品：辅助大气数据缺失或太阳天顶角不满足小于 76 度约束时，温度反演做不出来，只剩反射率。

读取与反算的标准写法：

import rasterio
import numpy as np

#-- 直接读 AWS 上的 COG，无需整景下载
base = 'LC08_L2SP_140041_20130503_20200911_02_T1'
url  = ('https://landsatlook.usgs.gov/data/collection02/'
        'level-2/standard/oli-tirs/2013/140/041/'
        f'{base}/{base}_SR_B4.TIF')

with rasterio.open(url) as src:
    dn = src.read(1)                      #-- 16 位无符号整型
    profile = src.profile
    win = src.read(1, window=rasterio.windows.Window(
        3000, 3000, 512, 512))            #-- 只取一块窗口

sr = dn.astype('float32') * 0.0000275 - 0.2
sr[dn == 0] = np.nan                      #-- 填充值置为缺测
print(dn.dtype, np.nanmin(sr), np.nanmax(sr))
print(profile['dtype'], profile['crs'], profile['width'], profile['height'])

上面这段的关键是最后两步：整型必须先乘系数加常数变成物理量，填充值 0 必须在运算前剔掉。把 0 当成一个"很黑的像元"直接喂进模型，会在景的边缘（条带外区域全是 0）制造出一批伪样本。

2.4 从景到时间序列：网格、层与集合类别

Landsat 的组织方式是三层嵌套：WRS 网格决定一景覆盖哪里，**集合类别（Tier）**决定这一景的几何质量，处理级别决定这一景是数字量还是物理量。

层	取值	含义	使用建议
处理级别	L1TP	精密地形校正，用地面控制点与高程模型正射	时序分析的首选
处理级别	L1GT	系统级地形校正	可用，几何精度弱于 L1TP
处理级别	L1GS	仅系统级几何校正	谨慎，仅当无替代时使用
处理级别	L2SP	反射率加温度都有	需要温度时认准这个
处理级别	L2SR	只有反射率	温度缺失，不要硬凑
集合类别	T1	径向均方根误差不超过 12 米	时序堆叠的默认选择
集合类别	T2	径向均方根误差超过 12 米	同一辐射标准，几何不达标
集合类别	RT	实时，采集后数小时内发布	应急用途，会被后续重处理取代
网格	WRS-1	1–3 号星使用，18 天重访	仅 1972–1983 年
网格	WRS-2	4–9 号星使用，16 天重访	覆盖 1982 年至今

实时层到 Tier 的转换有明确的时延：Landsat 8 约 14 至 16 天，Landsat 7 约 26 天，Landsat 9 约 3 至 4 天。这意味着"昨天刚过境的影像"在 Tier 里还查不到，做业务化系统时要把这段时延算进设计。新采集的 8 号与 9 号星数据在获取后 6 小时内即可检索与下载，走的是实时层。

2.5 与同类资源的规格差异

对比项	Landsat Collection 2	Sentinel-2（欧空局哥白尼计划）	MODIS（NASA Terra / Aqua）	EuroSAT 一类瓦片集
空间分辨率	30 米（全色 15 米）	10 至 60 米	250 至 1000 米	取决于源影像
重访周期	单星 16 天，双星 8 天	单星 10 天，双星 5 天	每天 1 至 2 次	不适用
时间起点	1972 年（Level-2 为 1982 年）	2015 年	2000 年前后	通常为单时相
波段数	4 至 11	13	36	13 或更少
数据量级	单景数百 MB 至 1 GB 以上	单景数百 MB	单景数十 MB	数 GB 全集
是否带标签	无	无	无	有
几何基准	Sentinel-2 全球参考影像	自有	自有	继承源影像

这张表要说清的只有一件事：Landsat 在"时间深度"这一列上唯一性突出，在"重访频率"与"分辨率"两列上都不占优。它与 Sentinel-2 的关系不是竞争而是接续——Collection 2 特意把几何基准重新对齐到由 Sentinel-2 生成的全球参考影像上，官方的说法是此前两颗星的像元在某些地区能差到 20 米以上，对齐后显著改善但没有彻底消除。做联合使用时，这一点要写进方法节。

2.6 各波段的物理含义

用法上的差异，根子在波段各自的波长位置。把 Landsat 8-9 的波段按用途排好，很多"为什么某个指数用这两个波段"的问题就有了答案：

波段	名称	波长范围（微米）	分辨率	主要用途
B1	海岸 / 气溶胶	0.43–0.45	30 米	水体、气溶胶反演
B2	蓝	0.45–0.51	30 米	真彩色合成、水体
B3	绿	0.53–0.59	30 米	真彩色合成、植被峰
B4	红	0.64–0.67	30 米	叶绿素吸收、归一化植被指数
B5	近红	0.85–0.88	30 米	生物量、归一化植被指数
B6	短波红外 1	1.57–1.65	30 米	含水量、云与雪区分
B7	短波红外 2	2.11–2.29	30 米	岩性、火烧迹地
B8	全色	0.50–0.68	15 米	锐化、形态判读
B9	卷云	1.36–1.38	30 米	卷云检测
B10	热红外 1	10.6–11.2	100 米（产品重采样到 30 米）	地表温度主用
B11	热红外 2	11.5–12.5	100 米（产品重采样到 30 米）	分裂窗，Collection 2 未用于温度反演

最后一行值得留意：Landsat 8-9 有两个热红外波段，设计意图是支持分裂窗温度反演，但 Collection 2 的地表温度走的是单通道算法，只用 B10。B11 的定标问题不在本文范围内，但使用者应当知道"两个热波段不等于两个温度产品"。

§3 它是谁的数据：一条五十年没有断过的记录

3.1 贡献者与机构分工

Landsat 是 NASA 与 USGS 的联合计划，两家分工明确：NASA 负责技术投入、仪器研制、发射前定标与在轨校验、空间段以及科学与应用；USGS 负责用户需求、在轨运行与数据获取、发射后定标、地面段（下传、处理、归档、管理与分发）与数据产品。数据分发由南达科他州苏福尔斯的 EROS 中心承担。

这个分工解释了一件看上去矛盾的事：为什么"卫星是 NASA 的，下载要去 USGS"。NASA 把造出来的仪器和平台交出去，USGS 从接手那一刻起负责把它变成可长期依赖的产品流。可持续陆地成像（SLI）计划把这件事写成了跨部门伙伴关系，目标是保证全球陆地成像数据免费、无歧视地获取，且时间上可比。

3.2 时间线

时间	事件	意义
1972-07-23	Landsat 1 入轨（原名 ERTS-1）	民用对地观测卫星起点
1975 / 1978	Landsat 2 / 3 入轨	MSS 时代，18 天重访
1982 / 1984	Landsat 4 / 5 入轨，搭载 TM	30 米与短波红外加入
1993-10-05	Landsat 6 未能入轨	计划里唯一一次失败
1999-04-15	Landsat 7 入轨，搭载 ETM+	15 米全色与 60 米热红外
2003-05-31	Landsat 7 扫描行校正器失效	此后单景缺失约 22%
2008	USGS 宣布整档免费开放	下载量出现百倍级增长
2013-02-11	Landsat 8 入轨，OLI 加 TIRS	首个业务化推扫陆地成像仪
2016	引入集合（Collections）管理制度	全档案重处理成为常规动作
2020-12	Collection 2 发布	全球 Level-2 成为常备库存
2021-09-27	Landsat 9 入轨	与 8 号联合 8 天重访
2022-12-30	Collection 1 停止下载	存量脚本需迁移
2030 年代初期	Landsat Next 计划入轨	三颗同构星、26 波段、6 天重访

Landsat 7 的退役时间在不同官方页面上口径不一：Earthshots 专题页给出的退役日期是 2025 年 6 月 4 日，而产品页说明其成像持续到 2022 年、此后卫星降轨。做长时序研究时以产品页的成像时段为准，退役日期只影响"什么时候彻底不再有新数据"。

3.3 采集场景：按网格而不是按需

Landsat 的采集是计划性的。卫星沿太阳同步轨道运行，轨道高度 705 公里，每 99 分钟绕地球一圈、每天约 14.5 圈，幅宽 185 公里，按 WRS-2 网格把陆地一块块扫过去。长期采集计划（LTAP）决定两颗在轨星如何分配观测优先级。

这个采集方式带来三个对使用者很现实的后果。第一，你不能在需要的时候点播一次观测——你只能在既有的过境日历里挑，挑不到就是没有。第二，观测时刻基本固定，上午约 10 点前后过境，这意味着太阳高度角与阴影形态在多年之间高度一致，对时序可比性是好事，对需要多角度观测的任务是限制。第三，云是硬约束，一次过境被云盖住，这次观测就作废，只能等十六天后的下一次；在多云地区，一年能拿到的有效观测次数会显著少于理论值。

3.4 数据来源的四个层次

理解这份数据要先分清四层：Level-1（辐射定标与正射校正后的数字量与表观量，1972 年至今）、Level-2（大气校正后的地表反射率与地表温度，1982 年至今，全球）、Level-3（由 Level-2 派生的专题产品，覆盖美国）、U.S. ARD（按统一投影切好的分析就绪瓦片，覆盖美国）。

四层的覆盖范围并不一致，这一点常被忽略：Level-2 是全球的，Level-3 与 ARD 只有美国（本土、阿拉斯加、夏威夷）。做中国、欧洲或非洲的研究时，能拿到的是 Level-1 与 Level-2，"分析就绪瓦片"这层要自己做。

§4 它最独特的东西：可比，而不是清晰

4.1 五十年同一套网格为什么稀缺

长时序遥感的最大敌人不是分辨率不够，而是不可比。两个年份的两张影像，如果大气校正算法不同、传感器定标不同、几何配准差了几十米，那么相减得到的"变化"里有很大一部分是处理差异而不是地面变化。这是为什么很多长时序研究宁可用粗糙但一致的传感器，也不用精细但不一致的。

USGS 用"集合"来根治这件事：一次全档案重处理定一个集合号，集合内所有景用同一套算法、同一套辅助数据、同一套几何基准重跑一遍；在下次重处理之前，已发布的像素不再变动。官方的说法很朴素——正是这种逐像素的确定性，让更高层级的时序分析成为可能。

所以它的稀缺性不在"拍得清楚"，而在"五十年里没有换过规矩"。这个前提一旦成立，很多本来做不了的研究就变得可做：冰川退缩的速率、海岸线的年际摆动、城市边界的逐年推进、森林从砍伐到复绿的全过程。

4.2 同一地点不同年份，像元能不能直接相减

能，但有前提：两景都要是 Tier 1，且来自同一集合。Tier 1 的门槛是径向均方根误差不超过 12 米，也就是像元中心位置的不确定性控制在三分之一像元以内。Tier 2 的辐射标准与 Tier 1 相同，几何不达标——差在哪儿？多半是轨道参数不够准（老卫星常见）、云太多导致匹配不上足够的控制点、或者地面控制不足。

你想做的事	需要哪一层	理由
逐年镶嵌、像元级时序	Tier 1，同一集合	亚像元级配准才有可比性
单时相分类、目视解译	Tier 1 或 Tier 2 均可	几何误差不影响单景判读
灾害应急、越快越好	实时层	接受几何待修正
跨传感器（TM 与 OLI）拼接	同一集合 Level-2	已做传感器间交叉定标
与 Sentinel-2 联合使用	Collection 2	几何基准已对齐到全球参考影像

最后一行要补一句限定：官方明确表示这次对齐"没有彻底解决问题，只是显著改善"。在两者混用做像元级分析时，仍建议先做一次局部配准检查。

4.3 案例：官方给出的三张对照图

USGS 在产品页上给了一组对照，正好把 Level-1 与 Level-2 的差别讲清楚：中国株洲 1999 年 9 月 7 日的 Landsat 5 影像，并排给出表观反射率、地表反射率与地表温度三张。同一天的同一块地，中间的差别全部来自大气——表观反射率里混着气溶胶散射与水汽吸收，看起来发灰发白；地表反射率把这些扣掉之后，植被与水体的对比立刻拉开；地表温度则是完全不同的物理量，量纲是开尔文，颜色条从 288.2 K 到 313.1 K。

另一组对照是美国华盛顿州布雷默顿，路径 47 行 27，Landsat 5 于 2010 年 10 月 6 日采集，同样是表观反射率与地表反射率并排。第三组是尼泊尔巴拉拉附近萨普塔科西河，路径 140 行 41，Landsat 8 于 2013 年 5 月 3 日采集，三张并排：Level-1 表观反射率、Level-2 地表反射率、Level-2 地表温度。

用代码把这件事复现一遍，比看说明文字有用：

import rasterio
import numpy as np

base = 'LC08_L2SP_140041_20130503_20200911_02_T1'
root = ('https://landsatlook.usgs.gov/data/collection02/'
        'level-2/standard/oli-tirs/2013/140/041/' + base + '/')

def read(suffix, scale=None, offset=None, fill=None):
    with rasterio.open(root + base + suffix + '.TIF') as src:
        a = src.read(1).astype('float32')
    if scale is not None:
        a = a * scale + offset
        a[a == fill * scale + offset] = np.nan
    return a

red  = read('_SR_B4', 0.0000275, -0.2, 0)      #-- 红
nir  = read('_SR_B5', 0.0000275, -0.2, 0)      #-- 近红
st   = read('_ST_B10', 0.00341802, 149.0, 0)   #-- 开尔文
qa   = read('_QA_PIXEL')                        #-- 位打包质量

#-- 归一化植被指数
ndvi = (nir - red) / (nir + red)

#-- 从 QA_PIXEL 拆出云与云影：bit3 云、bit4 云影
cloud     = (qa.astype('uint16') >> 3) & 1
shadow    = (qa.astype('uint16') >> 4) & 1
clear     = (cloud == 0) & (shadow == 0) & np.isfinite(ndvi)

print('NDVI 在晴空像元上的范围',
      round(float(np.nanmin(ndvi[clear])), 3),
      round(float(np.nanmax(ndvi[clear])), 3))
print('地表温度（摄氏度）范围',
      round(float(np.nanmin(st[clear]) - 273.15), 1),
      round(float(np.nanmax(st[clear]) - 273.15), 1))

这段代码里有两个容易被忽略的细节。一是 QA 波段必须按位拆，不能直接当数值比较——它的值 21824 以上代表某些填充与云的位组合，直接比较大小毫无意义。二是 fill * scale + offset 这个写法：填充值 0 经过线性变换之后不再等于 0，而是等于加常数，所以剔除填充必须在变换之后按原值剔除，或者在变换之前按 0 剔除。上面写成先变换再按变换后的填充值剔除，两种写法都对，但不能混。

4.4 质量波段的粒度：它替你做的与没替你做的

Landsat 的"标注"全部集中在质量波段上，把它们摊开看，才知道它替你做了什么、没做什么：

位	含义	谁用它
0	填充（条带外、无数据）	所有分析必须首先剔除
1	膨胀云	保守掩膜，云边外扩若干像元
2	卷云	薄云，反射率会被抬升
3	云	基本掩膜
4	云影	反射率被压暗，易误判为水体
5	雪 / 冰	与云在可见光上易混
6	晴空	该像元通过大气反演
7	水	水体掩膜
8–9	云置信度	低 / 中 / 高
10–11	云影置信度	低 / 中 / 高
12–13	雪冰置信度	低 / 中 / 高
14–15	卷云置信度	低 / 中 / 高

它替你做的是：把云、云影、雪、水、饱和、地形遮挡这些"物理上无法反演"的像元标出来，并给出置信度分级，让你能按任务松紧选择掩膜强度。它没替你做的是：任何语义上的类别。“这是耕地还是林地”“这是不是火烧迹地”——这些一概没有。把它当成带标签的数据集，是这条档案最常被误用的地方。

4.5 边界条件：什么情况下它不再适用

条件	是否适用	判断依据
做十年以上尺度的土地覆盖变化	适用，且几乎无可替代	Level-2 覆盖 1982 年至今
做 1972–1982 年的变化	谨慎	该时段只有 MSS Level-1，无 Level-2
做美国境内的专题分析	适用，可直接用 Level-3	水体范围、火烧迹地、雪盖比例已有官方产品
做美国之外的专题分析	适用 Level-2，专题需自建	Level-3 与 ARD 仅覆盖美国
做日尺度过程监测	不适用	重访 16 天（双星 8 天）
做亚米级目标识别	不适用	像元 30 米
需要夜间观测	不适用	太阳天顶角需小于 76 度，夜间无温度产品
做 Landsat 7 的 2003 年后镶嵌	谨慎	扫描行校正器失效，单景缺失约 22%

中间两行最容易被跳过。1972–1982 这一段只有 MSS 的表观反射率，官方明确说 MSS 的地表反射率要等到"业务化大气补偿算法成熟"才会提供——也就是说这一段至今没有 Level-2，任何跨 1982 年前后的长时序研究都会在这里出现算法断层。

§5 怎么拿到：从目录树到可训练张量

5.1 目录结构与下载面

云上的目录是按年、路径、行三层组织的，路径与行各补零到三位：

路径	内容
`s3://usgs-landsat/collection02/level-2/standard/oli-tirs/2022/026/039/<产品标识>/`	Landsat 8-9 Level-2 一景的全部文件
`s3://usgs-landsat/collection02/level-2/standard/etm/2005/136/041/<产品标识>/`	Landsat 7 Level-2
`s3://usgs-landsat/collection02/level-2/standard/tm/1995/136/041/<产品标识>/`	Landsat 4-5 Level-2
`s3://usgs-landsat/collection02/level-1/standard/oli-tirs/2023/...`	Level-1
`s3://usgs-landsat/collection02/level-2/standard/ard/...`	美国分析就绪瓦片

桶位于 AWS 的 us-west-2（俄勒冈）区域，是请求方付费桶：数据本身不收费，请求次数与出网流量由请求方承担。所有产品都配有 STAC 元数据，可通过 STAC API 检索而不触碰数据本体，这一点在筛选阶段能省下大量费用。

STAC 集合标识	内容
`landsat-c2l1`	Level-1
`landsat-c2l2-sr`	Level-2 地表反射率
`landsat-c2l2-st`	Level-2 地表温度
`landsat-c2ard-sr` / `-ta` / `-bt` / `-st`	美国 ARD 各分量
`landsat-c2l3-dswe` / `-ba` / `-fsca`	Level-3 三项专题产品

5.2 批量下载

#-- 先按 STAC 检索，只拿资产链接，不下载数据
python3.11 - <<'PY'
from pystac_client import Client
#-- USGS 的 STAC 端点无需凭据即可检索
cat = Client.open('https://landsatlook.usgs.gov/api/stac')
search = cat.search(
    collections=['landsat-c2l2-sr'],
    bbox=[113.0, 27.8, 114.0, 28.8],      #-- 株洲一带
    datetime='1999-01-01/2013-12-31',
    query={'platform': {'eq': 'landsat-5'}, 'eo:cloud_cover': {'lt': 20}})
items = list(search.items())
print(len(items), items[0].id if items else None)
with open('urls.txt', 'w') as f:
    for it in items:
        for k in ('red', 'nir'):
            f.write(it.assets[k].href + '\n')
PY

#-- 再按清单取对象；--request-payer 是请求方付费桶的必需参数
while read -r u; do
  aws s3 cp "$u" ./scenes/ --request-payer requester
done < urls.txt

第一段里 eo:cloud_cover 这个筛选字段来自景级元数据，能在下载之前把云量高的景剔掉，这是省钱最有效的一步。第二段逐条 aws s3 cp 而不是 aws s3 sync 整个目录，是因为后者会遍历你并不需要的大量景。

5.3 用窗口读取代替整景下载

COG（Cloud Optimized GeoTIFF）的用途就在这里：文件内部带金字塔与分块索引，可以只取一块窗口而不落地整景。

import rasterio
from rasterio.windows import Window

url = ('https://landsatlook.usgs.gov/data/collection02/level-2/'
       'standard/oli-tirs/2013/140/041/'
       'LC08_L2SP_140041_20130503_20200911_02_T1/'
       'LC08_L2SP_140041_20130503_20200911_02_T1_SR_B4.TIF')

with rasterio.open(url) as src:
    print(src.width, src.height, src.crs)       #-- 约 6300 乘 6000，UTM
    a = src.read(1, window=Window(3000, 3000, 256, 256))
    b = src.read(1, out_shape=(1, 128, 128))    #-- 概览层，快速预览
print(a.shape, b.shape)

训练阶段如果只需要小块样本，这种读法能把传输量降两个数量级。代价是每次取窗口都要发起一次或多次 HTTP 请求，样本数很大时请求数会成为瓶颈——这时候应该改成"整景拉到云主机上再切"，而不是"一块一块远程取"。

5.4 转成通用的表结构与掩膜口径

import numpy as np
import pandas as pd

#-- 假定 red / nir / qa 已按 §4.3 的方式读到同一窗口上
bits = qa.astype('uint16')
fill   = (bits >> 0) & 1
cloud  = (bits >> 3) & 1
shadow = (bits >> 4) & 1
snow   = (bits >> 5) & 1
water  = (bits >> 7) & 1

#-- 晴空陆地：无填充、无云、无云影、非雪非水
good = (fill == 0) & (cloud == 0) & (shadow == 0) & (snow == 0) & (water == 0)

ndvi = (nir - red) / (nir + red)
ndvi[~good] = np.nan

rows = []
h, w = ndvi.shape
for i in range(0, h, 64):
    for j in range(0, w, 64):
        blk = ndvi[i:i+64, j:j+64]
        rows.append({'row': i, 'col': j,
                     'ndvi_mean': float(np.nanmean(blk)),
                     'valid': int(np.isfinite(blk).sum())})
df = pd.DataFrame(rows)
print(df.shape, df['ndvi_mean'].describe())

这段代码只演示了从栅格到表的口径。上正式实验时有三处必须改：一是掩膜口径要写进方法节，"晴空陆地"的定义不同，结果就完全不同；二是分块大小要与任务匹配，六十四像元一格约合两公里，对城市任务太粗；三是这里只算了归一化植被指数，做分类模型时应保留全部波段而不是只留指数。

5.5 标定反算与温度换算

import numpy as np

def to_reflectance(dn, fill=0, scale=0.0000275, offset=-0.2):
    a = dn.astype('float32')
    a[dn == fill] = np.nan
    return a * scale + offset

def to_celsius(dn, fill=0, scale=0.00341802, offset=149.0):
    a = dn.astype('float32')
    a[dn == fill] = np.nan
    return a * scale + offset - 273.15

#-- 哨兵式的自检：反算后的物理值应落在文档给出的有效区间内
sr = to_reflectance(red_dn)
st = to_celsius(st_dn)
print('反射率 ', np.nanmin(sr), np.nanmax(sr))     #-- 应落在 0.0 至 1.0
print('摄氏度 ', np.nanmin(st), np.nanmax(st))     #-- 应落在 -123 至 86
print('不确定度', np.nanmin(st_qa * 0.01))         #-- ST_QA 另用 0.01 系数

最后一行是三个波段三种系数组的一个缩影：反射率用 0.0000275 与 −0.2，温度用 0.00341802 与 149.0，温度不确定度又用 0.01。三套系数不能通用，写反了不会报错，只会给出"看着合理、其实差十倍"的结果。

5.6 五个最容易出错的地方

位置	错误写法	后果	正确做法
数据集划分	随机划分像元或瓦片	地理泄漏，成绩虚高	按 WRS-2 路径行或地理块分组
反射率换算	把整型直接当反射率	数值差四个数量级	乘 0.0000275 减 0.2
填充值	把 0 当暗像元	条带外区域变成伪样本	反算前按 0 剔除，置为缺测
版本混用	同时用 Collection 1 与 2	前者填充 −9999、系数 0.0001	全链统一到 Collection 2
掩膜	只用云位不用膨胀云位	云边残留进入统计	保守任务叠加膨胀云位

§6 与其他数据集的关系

6.1 关系图谱

把遥感与环境类数据按"时间深度 × 空间粒度"两个轴摆开，Landsat 的位置很清楚：

数据集	时间深度	空间粒度	与 Landsat 的关系
EuroSAT	单时相瓦片	取决于源影像	下游：带标签的分类基准，与它构成"标签—原料"关系
多模态卫星公共卫生数据集	三年 156 周	卫星影像加地面指标	同类：用卫星时序配社会指标
逐日海温 OISST	四十余年逐日	四分之一度格网	互补：海面与陆地的两套长期记录
ERA5 单层再分析	数十年逐小时	数十公里格网	上游辅助：常作为大气与气象协变量
全球逐日气候观测 GHCNd	数十年逐日	地面站点	互补：点位真值，用于验证卫星反演
全球土壤栅格 SoilGrids	单时相或少数版本	250 米栅格	互补：静态属性层，常与它叠加
Overture Maps	持续更新	矢量要素	互补：矢量底图，用于样本框定
空气质量系统 AQS	数十年	地面站点	互补：点位观测，用于暴露建模

6.2 检索路径

你要找的东西	检索词	会落到哪里
带标签的遥感分类基准	`EuroSAT`、`land cover classification`	§6.1 表与 §8.2
更密的时间采样	`Sentinel-2`、`MODIS`、`VIIRS`	§2.5 规格对比
更高分辨率的影像	`high-resolution`、`commercial imagery`	§1.6 分辨率边界
大气与气象协变量	`ERA5`、`MERRA-2`、`GEOS-5`	§2.1 与 §7.3
地面的真值点位	`GHCNd`、`AQS`、`flux tower`	§7.4 外部验证
分析就绪的切块	`ARD`、`CEOS ARD`、`analysis ready`	§3.4 四层结构

6.3 与其他数据集的分工判据

任务	主粮仓该用什么	Landsat 的位置	原因
土地覆盖分类的预训练	带标签的瓦片集	提供原始波段，不提供标签	标签需自建
长时序变化检测	无可替代	主力	五十年同一套网格
高频过程监测	静止轨道或极轨高频传感器	不适用	重访不足
地表温度与热环境	可用	主力之一	1982 年至今的三十米温度
暴露与环境的协变量建模	地面站点加再分析	提供空间连续的自变量	点位与格网互补
目标级识别	亚米级影像	不适用	像元过大

这张表的第一行是大多数人的真实处境：他们需要一个带标签的基准来跑分类，而 Landsat 给的是没有标签的原料。被高估的用法是把它当分类基准直接训练，被低估的用法是用它构造"时间维度的对照"。

§7 质量评估与局限性

7.0 阅读本节前的一个前置校正

先纠正一个常见误判：把"没有标签"当成"质量差"。Landsat 不提供标签，是因为它是一个观测系统而不是一个标注工程；它的质量体现在辐射与几何的可追溯性上——每一个像元都能追到哪颗星、哪次过境、哪套算法、哪次重处理。评价它不能用"有没有标签"这把尺子，而应该问：在"可比"这件事上，它做到了什么程度？

7.1 五维评分

按常用维度打分，每项 1 至 5 分，分数针对"作为遥感算法研究原料"这个用途：

维度	分数	理由
数据规模	5	逾千万景，全球覆盖，五十年连续
标注质量	2	无语义标注，只有逐像元质量标记
任务适配	4	时序与变化检测几乎无可替代，分类任务需自建标签
文档完整	5	产品指南、格式控制书、数据字典、常见问题齐全
可获取性	3	免费无限制，但需注册账号，云上取数需付费账户

五项合计 19 分。分数结构说明它不是"开箱即训"的材料：规模与文档拉满，标注拖后腿，获取有门槛。

7.2 质量的两层

第一层是几何。Tier 1 的门槛是径向均方根误差不超过 12 米，也就是三分之一像元；Tier 2 超过这个数。做时序堆叠必须只用 Tier 1，且必须来自同一集合——跨集合比较等于跨处理版本比较。

第二层是辐射。反射率由 LaSRC 1.5.0（Landsat 8-9）与 LEDAPS 3.4.0（Landsat 4-7）生成，温度由单通道算法 1.3.0 生成。两套算法输入与方法都不相同：LEDAPS 用 6S 辐射传输模型、水汽与臭氧来自外部再分析、并把热波段用于云估计；LaSRC 用内部辐射传输模型、用海岸气溶胶波段做气溶胶反演、在反射率处理中不使用热波段。这意味着跨传感器拼接时，TM / ETM+ 与 OLI 的反射率差异里有一部分来自算法而不是地面。

7.3 可枚举的偏差

偏差	表现	影响
云与云影	一部分观测直接不可用	多云地区有效样本数远低于理论值
传感器与算法断层	1982 年前无 Level-2，MSS 至今无 Level-2	跨 1982 年的时序出现算法断点
Landsat 7 条带缺失	2003 年 5 月后单景缺失约 22%	镶嵌需多景填补或接受空缺
大气辅助数据切换	2023 年 2 月 17 日起改用新版气溶胶气候态	长时序中隐含一次辅助数据变更
温度在云边的已知误差	云与云影附近温度反演有系统性误差，文献已记载	热分析须结合 ST_CDIST 与 ST_QA 剔除
几何对齐不彻底	与 Sentinel-2 的像元偏移已大幅改善但未消除	联合使用需局部配准检查
过境时刻固定	上午约 10 点过境	无法反映日内变化，也无法按需点播
空间分辨率	三十米混合像元	目标级与形态级任务不适用

7.4 外部验证矩阵

验证角色	做法	讲得通的地方	风险
与地面站点对照	用地面观测校验反演值	点位真值独立于卫星	点位代表性与像元尺度不匹配
与再分析对照	用再分析场做协变量或交叉检验	空间连续、时间完整	再分析本身是模式输出，不是观测
与更高分辨率影像对照	用高分影像判读结果做样本	细节更丰富	高分影像时相不同，地物可能已变
与另一传感器对照	用同日过境的其他卫星交叉比对	独立观测链	波段与角度不同，需先做光谱转换

外部验证里最容易犯的错误是拿"不同日期的影像"当真值。地面变了，比对结果就既包含卫星误差也包含实际变化，两者分不开。可行的做法是锁定同日或相邻日期的过境，并把日期差写进误差预算。

7.5 它为何能持续：稀缺的不是分辨率，是连续性

它没有因为分辨率不如商业影像而被取代，原因是它占着一个别人没法短期重建的坑：五十年的连续记录。新增一条同等长度的记录需要再等五十年；新增一颗更高分辨率的卫星只需要几年。只要"长期可比"这个需求还在，它作为底片库的地位就不会被新传感器替代。Landsat Next 的设计也印证了这一点——它的卖点写的是"延续五十余年的记录"加上更高的时间、空间与光谱分辨率，而不是另起炉灶。

§8 生态位与历史

8.1 历史沿革

Landsat 的起点是 1960 年代中期：受无人遥感卫星在行星探测上的成功刺激，内政部、NASA 与农业部一起推动第一颗民用对地观测卫星，1972 年 7 月 23 日以"地球资源技术卫星"（ERTS-1）的名义发射，后来改称 Landsat 1。此后五十年，卫星换了九颗，传感器从 MSS 到 TM、ETM+、OLI / TIRS，分辨率从 80 米到 30 米，波段从 4 个到 11 个，唯一没变的是"同一套网格、同一套规矩"这个承诺。

2008 年是一个分水岭。此前影像按成本收费，一次订购的门槛足以劝退大多数研究；USGS 在 2008 年宣布整档免费，下载量随即出现百倍级增长。这个决定把 Landsat 从一个专业圈子内的资源变成了公共基础设施，也直接催生了此后十余年的遥感算法繁荣。

8.2 与基准数据集的关系

主流遥感分类基准，包括以 Sentinel-2 影像为源的十类土地覆盖瓦片集，本质上都是"从星载档案上切下来的、带标签的小图"。它们与 Landsat 是下游与上游的关系：前者解决"有没有标签"，后者解决"能不能追时间"。

一个常被混淆的点：做这类基准拿到好分数，不等于模型能处理长时序。基准里的瓦片是单时相、经过挑选的晴空样本，模型没见过云雾、没见过条带缺失、没见过跨传感器的辐射差异。把基准成绩直接外推到 Landsat 时序上，会遇到一批在基准里根本不存在的问题。

8.3 生态位：底片库而不是标签库

角色	标签库	底片库
典型用法	直接训练与评测分类模型	构造时序、变化检测、自监督语料
需要的前提	标签覆盖任务类别	可比性、可追溯性、全球覆盖
代表	各类土地覆盖瓦片集	Landsat、Sentinel-2、MODIS
失败模式	类别不全、地理泄漏	无标签、云掩膜与配准的前置工作量被低估

Landsat 坐在右边这一列。把它当标签库用，会立刻发现一个标签也没有；把标签库当底片库用，会立刻发现时间只有一层、地理泄漏无法回避。两者不能互换。

8.4 谁在用它，用来做什么

从公开材料看，使用者集中在四类：做土地覆盖与变化检测的研究组（把它当主数据源）、做农业与水资源管理的业务机构（用它算蒸散与作物耗水，美国农业部就是典型用户）、做灾害评估的应急机构（灾前灾后对比）、做机器学习的研究组（用它构造自监督与时序预训练语料）。第四类增长最快，因为无标签恰恰是自监督预训练最欢迎的形态。

§9 使用指引

9.1 三种角色的路径

角色	第一天	第一周	注意事项
算法研究者	读产品指南的波段与标定两章，跑通 §2.3	建立按路径行分组的划分，构造一个时序基线	先定掩膜口径再跑实验（§5.4）
数据工程师	用 STAC 检索并按云量筛选，跑通 §5.2	搭好窗口读取与缓存，估算请求费用	请求方付费，先算请求数再算数据量
业务分析人员	在 LandsatLook 上按地点看快视图，确认可用时相	用 Level-3 或自建指数产出专题图	美国之外没有 Level-3（§3.4）

9.2 一条最短的可用流程

#-- 1. 用 STAC 检索某地上一年云量低于 20% 的反射率景
python3.11 - <<'PY'
from pystac_client import Client
cat = Client.open('https://landsatlook.usgs.gov/api/stac')
s = cat.search(collections=['landsat-c2l2-sr'],
               bbox=[113.0, 27.8, 114.0, 28.8],
               datetime='2024-01-01/2024-12-31',
               query={'eo:cloud_cover': {'lt': 20}})
for it in list(s.items())[:5]:
    print(it.id, it.properties.get('eo:cloud_cover'))
PY

#-- 2. 取单景的红与近红，算归一化植被指数（见 §4.3 的 Python 段）
#-- 3. 记录完整产品标识，含采集日期与处理日期，便于复现

import rasterio
import numpy as np

#-- 复现的关键：把产品标识完整写进实验记录
pid = 'LC08_L2SP_140041_20130503_20200911_02_T1'
root = ('https://landsatlook.usgs.gov/data/collection02/'
        'level-2/standard/oli-tirs/2013/140/041/' + pid + '/')

def sr(band):
    with rasterio.open(f'{root}{pid}_SR_B{band}.TIF') as src:
        dn = src.read(1)
    out = dn.astype('float32') * 0.0000275 - 0.2
    out[dn == 0] = np.nan
    return out

ndvi = (sr(5) - sr(4)) / (sr(5) + sr(4))
print(pid, float(np.nanmean(ndvi)))

第三步看起来最琐碎，实际最要紧。产品标识里同时含采集日期与处理日期，重处理会改变后者；论文里只写"Landsat 8 影像"而不写完整标识，别人就无法确定你用的是哪一版。

9.3 引用方式

USGS 要求按数据集级别引用，并为三个传感器组分别给了 DOI：

引用对象	内容
Landsat 8-9 OLI / TIRS Level-2, Collection 2	EROS Center (2020)，DOI 10.5066/P9OGBGM6
Landsat 7 ETM+ Level-2, Collection 2	EROS Center (2020)，DOI 10.5066/P9C7I13B
Landsat 4-5 TM Level-2, Collection 2	EROS Center (2020)，DOI 10.5066/P9IAXOVV
Landsat 1-5 MSS Level-1, Collection 2	EROS Center (2020)，DOI 10.5066/P9AF14YV

写论文时建议同时注明：使用的集合号（Collection 2）、层级（Level-2）、集合类别（T1 或 T2）、以及是否使用美国 ARD 或 Level-3。官方对引用只提"请以如下方式引用"，未附加使用限制。

§10 AI 使用声明

项目	声明
文本生成	本文由 AI 辅助撰写，编辑把关结构与事实
数字来源	规模、结构、边界三类数字来自 NASA 与 USGS 官方页面、产品指南与科学数据目录
文件核验	官方产品页、产品指南、数据字典、科学数据目录与 EROS 活动页经检索核验；USGS 主站对自动请求返回人机验证，改由检索结果正文核对
不一致说明	Landsat 7 退役日期在不同官方页面口径不一，已在 §3.2 并列说明；档案景数在 2021 年与 2024 年官方材料上分别为九百万与逾一千万，已按时间顺序标注
论文信息	仅参考公开摘要与官方引用的文献条目，未读全文
未核验内容	各路径行的可用景数、逐景云量的具体数值、各年度重处理批次的完整清单
图像	无配图；封面为生成式提示词，不含真实卫星影像
免责	本文不构成业务或工程建议；数据使用须遵守官方引用规范与云上取数的计费规则

常见误读（避坑清单）

说法	实际	证据
「它是带标签的遥感分类数据集」	无任何人工标注，只有逐像元质量波段	产品波段清单里只有反射率、温度与各类 QA
「免费等于免注册」	EarthExplorer 与 M2M 需 USGS 账号，AWS 桶需账号且请求方付费	官方商业云访问页说明桶为请求方付费
「下载下来就能训」	档案在 PB 量级，应按云上计算或窗口读取组织	官方云访问指南与 COG 格式说明
「30 米能看清房子」	一个像元约合半个足球场，建筑只占数个像元	分辨率规格表
「1972 年就有地表反射率」	Level-2 从 1982 年起，MSS 至今没有 Level-2	Level-2 覆盖范围说明
「反射率整型可以直接用」	需乘 0.0000275 减 0.2，填充值 0 须先剔除	官方标定系数常见问题页
「Collection 1 还能下」	2022 年 12 月 30 日起停止下载	官方停服公告
「随机划分像元就行」	相邻像元共享同一景的天时地利，会造成地理泄漏	集合类别与网格结构说明
「Landsat 7 的影像是完整的」	2003 年 5 月 31 日后单景缺失约 22%	扫描行校正器失效影响评估文献
「地表温度就是气温」	是地表温度，量纲开尔文，与百叶箱气温不是一回事	Level-2 温度产品说明
「两个热红外波段有两个温度产品」	Collection 2 走单通道算法，只用第 10 波段	产品指南波段表
「美国之外的 Level-3 也有」	Level-3 与美国 ARD 只覆盖本土、阿拉斯加、夏威夷	ARD 与 Level-3 覆盖范围说明

FAQ

Q1：Landsat 到底有没有标签？
没有。它提供的是物理量（地表反射率、地表温度）与逐像元质量标记，不含任何地物类别。想做监督学习，标签要自建或外接其他土地覆盖产品；也可以退一步用官方的 Level-3 专题产品（水体范围、火烧迹地、雪盖比例）当弱标签，但那三项只覆盖美国。

Q2：为什么我下载的反射率图看起来一片黑？
多半是把整型当成了反射率。原始存储是 16 位无符号整型，必须乘 0.0000275 再减 0.2 才得到 0 到 1 之间的反射率。另一种可能是把填充值 0 当成了有效像元——条带外的区域全是 0，直接显示就是黑边。

Q3：1982 年之前的数据能不能用？
能用，但只有 Level-1。1972 至 1982 年是 MSS 传感器，官方至今没有提供 MSS 的地表反射率产品，理由是业务化大气补偿算法尚未满足计划要求。跨 1982 年前后的长时序研究要在这里处理一次算法断层，不能把 Level-1 的表观反射率与 Level-2 的地表反射率混在同一条曲线上。

Q4：训练集和测试集该怎么分？
按 WRS-2 的路径行号分组，或者按地理块分组，不要按随机像元分。同一个路径行的相邻像元来自同一次过境，共享同一片大气与同一套光照，随机划分会把它们分到两侧，成绩里会掺进记忆而非泛化。产品标识里就带着路径行号，是最方便的分组键。

Q5：云太多怎么办？
两条路。一是放宽时间窗口，把多年同期影像合成一期中值影像，用中值而不是单景；二是换更密的传感器。Landsat 单星十六天、双星八天，理论年观测次数有限，在多云地区实际可用次数会更少，这是设计约束而不是数据缺陷。

Q6：地表温度怎么换算成摄氏度？
先用 0.00341802 乘、加 149.0 得到开尔文，再减 273.15。填充值 0 要先剔除。同时建议读 _ST_QA（系数 0.01，单位开尔文）看不确定度，并读 _ST_CDIST 看该像元到最近云的距离——云边附近的温度反演有已知的系统性误差，文献已记载。

Q7：Landsat 和 Sentinel-2 能不能混用？
能，但要知道 Collection 2 做了一件具体的事：把几何基准重新对齐到由 Sentinel-2 生成的全球参考影像上。官方说这显著改善了两者的像元对齐，但"没有彻底解决问题"。做像元级联合分析前，建议先在局部做一次配准检查。

Q8：Collection 1 还能不能用？
2022 年 12 月 30 日起 USGS 停止提供 Collection 1 下载。存量脚本要迁移到 Collection 2，迁移时注意三处变化：填充值从 −9999 变成 0，数据类型从有符号变成无符号，反射率系数从 0.0001 变成 0.0000275 并加了 −0.2 的加常数。

Q9：档案有多大？能一次下完吗？
不能。USGS EROS 档案存量逾一千万景，每天新增约一千五百景，整体在 PB 量级。可行做法是先用 STAC 按时间、地点与云量筛选出需要的几十到几百景，再决定是拉到云主机上处理，还是用 COG 窗口读取只取需要的块。

Q10：商业使用需要授权吗？
官方声明产品无使用限制，USGS 科学数据目录以 CC0-1.0 公共领域奉献标注。但"无使用限制"指的是数据本身，不覆盖取数通道的费用：AWS 桶是请求方付费，EarthExplorer 与 M2M 需要注册账号。商业项目应把这两类成本分开核算。

Q11：Landsat Next 什么时候能用？
官方定位是 2030 年代初期入轨，三颗同构卫星、26 个波段、空间分辨率 10 至 60 米、三星联合 6 天重访、采用新的 WRS-3 网格。网格一换，与既有档案的像元对应关系就要重新定义，这是未来十年长时序研究必须提前考虑的一件事。

术语表

术语	解释
WRS-1 / WRS-2	全球参考系统，Landsat 的轨道网格编号，由路径号与行号定位一景
集合（Collection）	USGS 的档案分级管理制度，一次全档案重处理定一个集合号
Tier 1 / Tier 2 / RT	集合类别，分别对应几何误差不超过 12 米、超过 12 米、实时未定级
Level-1	辐射定标与正射校正后的产品，含数字量与表观反射率
Level-2	大气校正后的地表反射率与地表温度
Level-3	由 Level-2 派生的专题产品，含水体范围、火烧迹地、雪盖比例
U.S. ARD	美国分析就绪数据，按统一等积投影切好的瓦片
COG	云优化 GeoTIFF，内部带金字塔与分块索引，支持窗口读取
地表反射率	地物反射太阳辐射的比例，已扣除大气影响，无量纲
表观反射率	传感器入瞳处测到的反射率，含大气贡献
地表温度	地表辐射温度，单位开尔文，不等于气温
LaSRC / LEDAPS	两套大气校正算法，分别用于 Landsat 8-9 与 Landsat 4-7
太阳天顶角约束	Level-2 生成的前置条件，需小于 76 度
STAC	时空资产目录，用于检索卫星数据的元数据规范
请求方付费	AWS 桶的计费方式，数据免费但请求与流量由请求方承担
扫描行校正器	Landsat 7 上补偿卫星前进运动的部件，2003 年失效
MSS / TM / ETM+ / OLI / TIRS	五代传感器：多光谱扫描仪、专题制图仪、增强型专题制图仪、陆地成像仪、热红外传感器
地理泄漏	同一地理区域同时出现在训练与测试集，导致成绩虚高
归一化植被指数	近红与红波段之差除以二者之和，植被丰度的常用指标
EROS	地球资源观测与科学中心，USGS 的归档与分发机构
Landsat Next	下一代计划，三颗同构星、26 波段、6 天重访

附录 A：核心数字核对表

项	值	来源
档案存量	逾 1000 万景	NASA Landsat Next 报告（2024）；2021 年官方口径为逾 900 万
每日新增	约 1500 景（8 号与 9 号合计）	USGS 卫星概览说明书
每日覆盖面积	约 4000 万平方公里	NASA Landsat Next 报告
累计下载	2020 年突破 1 亿景	NASA 五十周年专题
免费开放起始	2008 年	NASA 与 USGS 多处说明
Level-1 起始	1972 年 7 月 23 日	Landsat 1 入轨日期
Level-2 起始	1982 年	USGS Collection 2 说明
MSS 数据区间	1972-07-25 至 2013-01-07	USGS 科学数据目录
空间分辨率	30 米（全色 15 米）	USGS 卫星概览说明书
MSS 分辨率	80 米采集，产品重采样到 60 米	USGS 卫星概览说明书
单景覆盖	约 170 km 南北乘 183 km 东西	USGS Landsat 7 产品页
单景像元数（8 号）	约 6300 乘 6000	USGS ARD 瓦片说明
重访周期	单星 16 天，双星 8 天	USGS 卫星概览说明书
轨道高度	705 公里（1–3 号星为 920 公里）	USGS 卫星概览说明书
Tier 1 几何门槛	径向均方根误差不超过 12 米	USGS 集合说明书与数据字典
反射率标定	乘 0.0000275 减 0.2，填充 0	USGS 常见问题页与产品指南
温度标定	乘 0.00341802 加 149.0，填充 0	USGS 常见问题页与产品指南
温度有效区间	150 至 359 开尔文	产品指南波段表
反射率算法版本	LaSRC 1.5.0 / LEDAPS 3.4.0	USGS 地表反射率页
温度算法版本	单通道算法 1.3.0	USGS Level-2 科学产品页
ARD 瓦片大小	5000 乘 5000 个 30 米像元	USGS ARD 事实说明书
ARD 瓦片数	本土 422、阿拉斯加 110、夏威夷 8	USGS ARD 产品页
SLC 失效	2003 年 5 月 31 日，单景缺失约 22%	USGS SLC 影响评估文献
大气辅助切换日	2023 年 2 月 17 日	产品指南
实时层转 Tier	8 号 14–16 天，7 号 26 天，9 号 3–4 天	USGS ARD 事实说明书
新数据可检索	获取后 6 小时内	USGS Level-1 产品页
太阳天顶角约束	小于 76 度	USGS Level-2 说明
许可	无使用限制；CC0-1.0	USGS 产品页与科学数据目录
Collection 2 发布	2020 年 12 月	USGS EROS 活动页
Collection 1 停服	2022-12-30	USGS 产品页公告
Landsat Next	三颗星、26 波段、6 天重访、WRS-3	USGS Landsat Next 事实说明书

附录 B：数字的来源分层

第一层：官方页面与官方产品文档。 档案存量、分辨率、重访、轨道、单景尺寸、标定系数、有效区间、算法版本、瓦片规格、集合类别门槛、覆盖范围、许可与日期。这些数字来自 NASA 与 USGS 的项目页、产品指南、数据字典、事实说明书与科学数据目录记录。

第二层：官方引用文献的摘要。 扫描行校正器失效的影响评估、温度在云边的已知误差，这两项来自 USGS 与 NASA 官方页面上引用的技术文献，本文只引用其结论性表述，未做量化复算。

第三层：未核实。 各路径行的可用景数、逐景云量的具体数值、各年度重处理批次的完整清单、中国境内各年份的有效观测次数。这些需要按路径行逐个查询统计，超出了本篇的范围，正文中凡是涉及的地方都只做定性表述。

附录 C：完整参考文献

编号	文献	链接	用途
1	Landsat 计划主页（NASA）	https://landsat.gsfc.nasa.gov/	计划定位、波段、连续性与免费开放起始
2	Landsat Collection 2 Level-2 科学产品页（USGS）	https://www.usgs.gov/land-resources/nli/landsat/level-2-science-products	算法版本、引用 DOI、无使用限制声明
3	Landsat 8-9 C2 L2SP 产品指南 LSDS-1619	https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/media/files/LSDS-1619_Landsat8-9-Collection2-Level2-Science-Product-Guide-v5.pdf	波段表、标定系数、有效区间、命名规则
4	Landsat Collection 2 事实说明书（USGS，2021）	https://pubs.usgs.gov/fs/2021/3002/fs20213002.pdf	Collection 1 与 2 的逐项差异
5	Landsat 集合说明书（USGS，2018）	https://pubs.usgs.gov/fs/2018/3049/fs20183049.pdf	Tier 定义与 12 米门槛
6	标定系数与填充值常见问题（USGS）	https://www.usgs.gov/faqs/why-are-fill-values-and-scaling-factors-landsat-collection-2-level-2-products-different-those	反射率与温度的系数组
7	Landsat Collection 2 数据字典（USGS）	https://www.usgs.gov/centers/eros/science/landsat-collection-2-data-dictionary/	产品标识与景标识字段定义
8	地球观测卫星概览说明书（USGS，2015）	https://pubs.usgs.gov/fs/2015/3081/fs20153081.pdf	各卫星时段、传感器、分辨率、SLC 失效
9	Landsat 项目历史与卫星退役表（USGS Earthshots）	https://eros.usgs.gov/earthshots/the-landsat-program	九颗星的发射与退役日期
10	扫描行校正器失效影响评估（USGS）	https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/atoms/files/SLC Failure Image Impact Paper V1.1.pdf	失效日期与 22% 缺失比例
11	美国 ARD 事实说明书（USGS，2023）	https://pubs.usgs.gov/fs/2023/3015/fs20233015.pdf	瓦片规格、覆盖范围、处理时延
12	商业云数据访问页（USGS）	https://www.usgs.gov/landsat-missions/landsat-commercial-cloud-data-access	S3 桶位置、请求方付费、STAC 与 SNS
13	Landsat 8-9 OLI/TIRS Level-2 科学数据目录记录	https://data.usgs.gov/datacatalog/data/USGS:EROS5e83d14f2fc39685	覆盖范围、CC0-1.0 标注、DOI 10.5066/P9OGBGM6
14	Landsat 1-5 MSS Level-1 科学数据目录记录	https://data.usgs.gov/datacatalog/data/USGS:EROS5e83d0e09a752cff	MSS 数据区间与 DOI 10.5066/P9AF14YV
15	Landsat Next 事实说明书（USGS，2024）	https://pubs.usgs.gov/fs/2024/3005/fs20243005.pdf	三星、26 波段、6 天重访、WRS-3
16	Landsat 五十周年专题（NASA）	https://www.nasa.gov/earth-and-climate/landsat-legacy-nasa-usgs-program-observing-earth-from-space-turns-50/	2008 年免费开放与累计下载量

附录 D：任务代码

任务一：按路径行分组的时序构造骨架

import re
import numpy as np
import pandas as pd

#-- 从产品标识解析路径行与日期
PAT = re.compile(r'^L(?P<sensor>[CTEOM])(?P<sat>\d{2})_'
                 r'(?P<lvl>L\d\w{2})_(?P<path>\d{3})(?P<row>\d{3})_'
                 r'(?P<acq>\d{8})_(?P<proc>\d{8})_(?P<coll>\d{2})_(?P<tier>\w+)$')

def parse(pid):
    m = PAT.match(pid)
    if not m:
        raise ValueError(pid)
    d = m.groupdict()
    return {'pid': pid, 'path': d['path'], 'row': d['row'],
            'pathrow': d['path'] + d['row'],
            'acq': d['acq'], 'proc': d['proc'], 'tier': d['tier']}

ids = ['LC08_L2SP_140041_20130503_20200911_02_T1',
       'LE07_L2SP_140041_20120424_20200908_02_T1',
       'LT05_L2SP_140041_20111008_20200903_02_T1']
meta = pd.DataFrame([parse(i) for i in ids])

#-- 只保留 Tier 1，并按路径行分组做留一验证
meta = meta[meta['tier'] == 'T1']
for pr, grp in meta.groupby('pathrow'):
    print('路径行', pr, '可用景数', len(grp),
          '采集日期', sorted(grp['acq'].tolist()))

任务二：多年同期的晴空合成

import numpy as np

#-- 假定 scenes 是一组已反算的 (日期, red, nir, qa) 三元组
stack = []
for date, red, nir, qa in scenes:
    bits = qa.astype('uint16')
    ok = ((bits >> 0) & 1) == 0     #-- 非填充
    ok &= ((bits >> 3) & 1) == 0    #-- 非云
    ok &= ((bits >> 4) & 1) == 0    #-- 非云影
    ndvi = (nir - red) / (nir + red)
    ndvi[~ok] = np.nan
    stack.append(ndvi)

arr = np.stack(stack)
med = np.nanmedian(arr, axis=0)     #-- 取多年中值，抵消残余云与噪声
cnt = np.isfinite(arr).sum(axis=0)  #-- 每个像元的可用年数

print('合成后有效像元比例', float(np.isfinite(med).mean()))
print('可用年数的最小值', int(cnt.min()), '最大值', int(cnt.max()))

第二段代码里的 cnt 比 med 更值得先看一眼。合成图上每一个像元的"可用年数"差别很大，常年有云的山谷可能只有两三年可用，而干旱地区可能有十几二十年。把这个计数图一起输出，才能让下游知道哪些像元的时序是可靠的。

附录 E：如果只记五件事

序号	记住什么	为什么
1	1972 年至今连续，Level-2 从 1982 年起，三十米，十六天	这四个数字定义了一切用法与限制
2	一个人工标注都没有	它是底片库不是标签库，标签要自建
3	整型要乘 0.0000275 减 0.2，填充值 0 须剔除	不做这一步，所有反射率都错
4	按路径行分组划分，不要随机分像元	地理泄漏会让结果虚高
5	档案在 PB 量级，用 STAC 筛选加 COG 窗口读取	它在云上算，不在本地下

编辑说明

本文的核验时间为 2026 年 10 月 29 日。核对方式为直接访问 NASA Landsat 计划页（返回 200）与 USGS 的产品页、产品指南、数据字典、事实说明书与科学数据目录；USGS 主站对自动请求返回人机验证页，相关页面改由检索结果正文核对。两处官方口径不一致已在正文标出：一是 Landsat 7 的退役日期（Earthshots 专题页为 2025 年 6 月 4 日，产品页称成像持续到 2022 年），二是档案景数（2021 年为逾九百万，2024 年为逾一千万，本文按时间顺序分别标注）。未覆盖的内容包括：各路径行的可用景数、逐景云量的具体数值、各年度重处理批次的完整清单、以及中国境内各年份的有效观测次数。这些需要按路径行逐个查询统计，本文按要求只写规模、结构、边界三类数字。站内自链接以 __RID__ 形式留空，发布时由流水线注入真实编号。