Landsat 8/9 影像从下载到预处理完整流程：辐射定标、大气校正、裁剪一篇搞定

一、Landsat 8/9 数据体系与预处理逻辑起点

1.1 传感器谱系与数据产品层级

Landsat 8 于 2013 年 2 月 11 日发射，Landsat 9 于 2021 年 9 月 27 日发射，两者均搭载 OLI（Operational Land Imager）和 TIRS（Thermal Infrared Sensor）两台载荷。OLI 提供 9 个光谱波段，涵盖可见光、近红外、短波红外与卷云波段，空间分辨率 30 米（全色波段 15 米）；TIRS 提供两个热红外波段，原始分辨率 100 米，重采样至 30 米分发。Landsat 9 的 OLI-2 在辐射分辨率上从 12 位提升至 14 位，这一变化直接影响了暗目标区域的信噪比表现。

USGS 将 Landsat 产品划分为多个处理层级：L0 为原始遥测数据，L1 为经过辐射定标和几何校正的辐射亮度或表观反射率产品，L2 为经过大气校正的地表反射率产品，L3 为更高层次的地球物理产品。自 2020 年起，USGS 开始为 Landsat 8/9 提供 Collection 2 Level-2 产品，这意味着用户可以直接获取大气校正后的地表反射率数据，无需自行运行校正算法。

1.2 预处理链条的“保真度”主线

如果将遥感影像的定量应用比作一条信息传递链，那么预处理就是这条链上信噪比衰减最剧烈的环节。本文提出以“数据保真度”（Data Fidelity）作为贯穿全文的分析主线，其核心内涵是：在每一个预处理步骤中，尽可能减少对原始辐射信息的非必要改变，同时将必要的物理转换做到可追溯、可验证、可复现。

保真度的损失主要来自三个渠道：一是定标系数引入的系统误差，二是大气校正模型对气溶胶参数的敏感性，三是重采样与裁剪过程中的空间信息损失。这三个渠道分别对应本文的第三、四、六章。理解这一主线，有助于读者在每一步操作中建立“误差预算”意识，而非机械地执行流程。

二、数据下载：从 EarthExplorer 到云原生接口

2.1 传统下载路径与数据格式

USGS EarthExplorer 是最经典的数据获取入口，支持按路径/行号（WRS-2）、经纬度、行政区划、云量阈值等条件检索。下载得到的 L1 产品为 GeoTIFF 格式，每个波段一个文件，附带 MTL.txt 元数据文件和 BQA 质量波段。自 Collection 2 起，L1 产品还包含 ANG.txt 角度文件，记录了太阳天顶角、方位角以及传感器观测天顶角、方位角，这些角度信息是后续大气校正和地形校正的必要输入。

下载时需注意几个关键参数：云量阈值建议设置在 10%–20% 之间，过低会导致可用影像稀少，过高则增加后续云掩膜的工作量；数据级别选择 L1TP（Precision Terrain Corrected）而非 L1GT（Geometric Terrain Corrected），前者使用地面控制点和 DEM 进行精校正，几何精度通常在 12 米以内。

2.2 云原生接口与批量获取

近年来，USGS 与 NASA 联合推出了 Landsat 数据云原生访问方案。通过 AWS 的公开数据集（Registry of Open Data on AWS），用户可以直接在 S3 存储桶中访问 Landsat 8/9 的 Collection 2 数据，无需逐景下载。这一方案的核心优势在于：数据以 Cloud Optimized GeoTIFF（COG）格式存储，支持 HTTP Range 请求，可以只读取影像的特定窗口而非整景下载。

对于需要批量处理的用户，USGS 提供了 Machine-to-Machine（M2M）API，支持通过编程方式提交检索请求、获取下载链接。此外，Google Earth Engine 平台也集成了 Landsat 8/9 的 L1 和 L2 数据，用户可以直接在平台上进行云端计算，避免了本地存储和算力的瓶颈。

三、辐射定标：DN 值到表观反射率的物理转换

3.1 辐射定标的物理基础

辐射定标是将传感器记录的原始数字量化值（DN）转换为具有物理意义的辐射亮度或表观反射率的过程。其理论基础是传感器响应线性假设：在动态范围内，传感器输出的 DN 值与入射辐射亮度呈线性关系。USGS 为每个波段提供了两个关键定标系数：RADIANCE_MULT_BAND_x（增益）和 RADIANCE_ADD_BAND_x（偏移），以及 REFLECTANCE_MULT_BAND_x 和 REFLECTANCE_ADD_BAND_x。

辐射亮度转换公式为：Lλ = ML × Qcal + AL，其中 Lλ 为传感器处辐射亮度（W/(m²·sr·μm)），ML 为增益，AL 为偏移，Qcal 为 DN 值。表观反射率转换公式为：ρλ' = Mρ × Qcal + Aρ，再除以 sin(θSE) 进行太阳高度角校正，其中 θSE 为太阳高度角。

# 辐射定标核心代码（Python + rasterio）
import rasterio
import numpy as np

def dn_to_radiance(dn_path, mult, add):
    with rasterio.open(dn_path) as src:
        dn = src.read(1).astype(np.float32)
        radiance = mult * dn + add
    return radiance

def dn_to_reflectance(dn_path, mult, add, sun_elevation):
    with rasterio.open(dn_path) as src:
        dn = src.read(1).astype(np.float32)
        rho = (mult * dn + add) / np.sin(np.radians(sun_elevation))
    return rho

3.2 定标系数的时间演变

Landsat 8/9 的定标系数并非一成不变。USGS 会定期根据传感器在轨表现更新定标参数，以补偿探测器老化带来的响应衰减。以 Landsat 8 OLI 为例，自 2013 年发射以来，USGS 已发布多次定标系数更新，其中 2016 年和 2020 年的更新幅度较大。如果使用过时的定标系数处理近期影像，可能导致辐射亮度出现 2%–5% 的系统偏差。

MTL.txt 文件中记录的定标系数是处理该景影像时的“正确”系数，因此直接读取 MTL 文件是最稳妥的做法。但需要注意的是，如果研究需要长时间序列的一致性，则应使用 USGS 发布的最新定标系数重新处理所有影像，而非依赖各景 MTL 中的历史系数。

四、大气校正：从理论模型到工程实现

4.1 大气校正的物理机制

大气校正是预处理链条中误差贡献最大的环节。太阳辐射在到达地表之前，会经历瑞利散射、气溶胶散射和水汽吸收等过程；地表反射后，辐射再次穿过大气到达传感器。大气校正的目标就是消除这些大气效应，反演出地表反射率。

辐射传输方程可以简化为：LTOA = Lpath + (ρ × T↓ × T↑ × Edown) / (π × (1 - ρ × S))，其中 Lpath 为路径辐射，T↓ 和 T↑ 分别为下行和上行透过率，Edown 为下行辐照度，S 为球面反照率。大气校正的核心挑战在于：Lpath、T 和 Edown 都依赖于大气状态参数（气溶胶光学厚度、水汽含量、臭氧含量等），而这些参数在空间和时间上都是变化的。

4.2 主流校正方法对比

目前主流的大气校正方法可分为三类：基于辐射传输模型的方法（如 FLAASH、6S、MODTRAN）、基于暗像元的方法（如 DOS、LaSRC）和基于统计的方法（如 QUAC）。FLAASH 集成在 ENVI 平台中，使用 MODTRAN 辐射传输代码，支持多种气溶胶模型和大气模式；6S 模型是开源辐射传输代码，计算效率高，适合批量处理；LaSRC 是 USGS 官方 L2 产品使用的校正器，基于 6S 模型改进，利用 MODIS 气溶胶产品辅助估算。

4.3 6S 模型工程化配置

6S（Second Simulation of the Satellite Signal in the Solar Spectrum）模型是开源大气校正的基石。其输入参数包括几何参数（太阳天顶角、方位角、观测天顶角、方位角）、大气模式（热带、中纬度夏/冬、亚北极夏/冬）、气溶胶类型（大陆型、海洋型、城市型、沙漠型）和气溶胶光学厚度（550nm）。

工程化配置的关键在于气溶胶光学厚度的获取。对于无实测数据的区域，可以从 AERONET 站点获取，或使用 MODIS 的 MOD04 气溶胶产品，或采用暗像元法从影像自身估算。笔者在实践中发现，对于中国东部地区，使用 MODIS 气溶胶产品结合暗像元法进行交叉验证，可以将气溶胶光学厚度的不确定性控制在 ±0.05 以内。

五、几何精校正与地形校正

5.1 几何精校正的必要性

Landsat L1TP 产品已经过系统几何校正和地形校正，理论几何精度在 12 米以内。但在实际应用中，由于 DEM 精度限制和控制点分布不均，局部区域仍可能存在几何偏差。对于变化检测、影像融合等应用，亚像元级的几何配准是必要的。

几何精校正的流程包括：选择基准影像、选取控制点、计算变换模型、重采样。控制点的选择应遵循“均匀分布、特征明显、数量充足”的原则，一般每景影像需要 20–30 个控制点。变换模型可选择多项式模型（一阶、二阶）或有理多项式系数（RPC）模型。重采样方法中，最近邻法不改变像元值但会产生锯齿，双线性内插法平滑但会改变光谱值，三次卷积法在两者之间取得平衡。

5.2 地形校正方法

在山区，地形起伏导致太阳入射角和观测角发生剧烈变化，同一地物类型在不同坡向上的反射率差异可达 2–3 倍。地形校正的目标就是消除这种地形效应，使地表反射率更接近真实值。常用方法包括 Cosine 校正、C 校正、Minnaert 校正和 SCS 校正。

Cosine 校正最为简单，假设地表为朗伯体，将反射率除以太阳入射角余弦。但该方法在阴影区域会产生过度校正。C 校正引入经验参数 C 来调节校正强度，对阴影区域更为稳健。Minnaert 校正引入非朗伯体假设，适用于植被覆盖区域。SCS 校正则同时考虑太阳入射角和传感器观测角，适用于复杂地形。

本文评述认为，地形校正方法的选择应基于研究区的地形特征和地物类型。对于地形起伏较小的平原地区，地形校正的收益有限；对于高山峡谷区域，SCS 校正通常表现最优，但其对 DEM 精度的要求也最高。笔者建议在应用地形校正前，先使用统计分析（如反射率与太阳入射角的回归分析）判断地形效应的显著程度，避免“为校正而校正”。

六、影像裁剪与掩膜处理

6.1 裁剪的两种策略

影像裁剪看似简单，实则涉及“保真度”的微妙权衡。裁剪策略可分为两类：按矢量边界裁剪和按矩形范围裁剪。按矢量边界裁剪可以精确提取研究区，但边界像元会因部分覆盖而产生混合像元问题；按矩形范围裁剪则保留了完整的像元，但包含研究区外的区域。

对于定量遥感应用，笔者建议采用“先矩形裁剪、后矢量掩膜”的两步策略：先用矩形范围裁剪减少数据量，再用矢量边界生成掩膜，将研究区外的像元设为 NoData。这样既保证了边界像元的完整性，又避免了混合像元对统计分析的干扰。

6.2 云掩膜与阴影处理

Landsat Collection 2 产品提供了 QA_PIXEL 质量波段，其中包含了云、云阴影、雪、水体等标记信息。QA_PIXEL 采用位运算编码，需要按位解析。例如，第 3 位为云标记，第 4 位为云阴影标记，第 5 位为雪标记。

# QA_PIXEL 位运算解析（Python）
import numpy as np

def parse_qa(qa_band):
    cloud = (qa_band & (1 << 3)) > 0
    shadow = (qa_band & (1 << 4)) > 0
    snow = (qa_band & (1 << 5)) > 0
    water = (qa_band & (1 << 7)) > 0
    return cloud, shadow, snow, water

近年来，基于深度学习的云检测方法发展迅速。以 Landsat 8/9 为对象的云检测模型（如 CloudNet、SegNet 变体）在精度上已超过传统阈值方法，尤其是在薄云和碎云检测方面。但深度学习方法需要大量标注数据，且模型泛化能力受训练区域限制。本文评述认为，对于常规应用，QA_PIXEL 波段已足够；对于云污染严重的区域，可考虑结合 Fmask 算法或深度学习模型进行补充。

七、批量处理工程化：从脚本到流水线

7.1 自动化脚本设计

单景影像的预处理流程可以通过 GUI 软件完成，但面对数十景甚至数百景的批量处理任务，自动化脚本是唯一可行的方案。一个健壮的批量处理脚本应具备以下特征：参数可配置、错误可捕获、进度可追踪、结果可验证。

以 Python 为例，可以使用 rasterio 进行影像读写，使用 py6s 调用 6S 模型，使用 geopandas 处理矢量边界。脚本的核心逻辑是：遍历影像列表 → 读取 MTL 元数据 → 辐射定标 → 大气校正 → 几何裁剪 → 输出结果。每一步都应记录日志，便于排查问题。

7.2 并行计算与云平台

对于大规模批量处理，单机串行效率低下。可以使用 Python 的 multiprocessing 模块进行多进程并行，或使用 Dask 进行分布式计算。在云平台上，AWS Batch、Google Cloud Run 等服务可以按需启动计算资源，处理完成后自动释放。

Google Earth Engine 提供了另一种思路：将预处理算法部署到云端，用户只需提交处理请求，平台自动完成计算。GEE 内置了 Landsat 8/9 的 L1 和 L2 数据，以及辐射定标、大气校正（使用 LaSRC）、云掩膜等函数，用户可以直接调用。但 GEE 的封闭性也意味着用户无法完全控制校正参数，对于需要精细调参的研究，仍需本地处理。

八、质量评价与精度验证

8.1 定性评价方法

预处理结果的质量评价可以从定性和定量两个维度进行。定性评价包括：目视检查影像是否清晰、有无明显条带或噪声、水体是否呈现深色、植被是否呈现红色（假彩色合成）、云影是否被有效掩膜。定量评价则包括：统计反射率的动态范围、检查有无异常值（如负值或超过 1 的值）、与官方 L2 产品进行对比。

8.2 定量验证方法

定量验证的“金标准”是与地面实测反射率进行对比。地面实测可使用 ASD FieldSpec 等光谱仪，在影像过境前后 1 小时内完成测量。验证点应选择均质、平坦、无云的区域，如大面积水体、裸土或草地。

在没有地面实测数据的情况下，可以使用交叉验证方法：将自行校正的结果与 USGS 官方 L2 产品进行逐像元对比，计算均方根误差（RMSE）和平均绝对误差（MAE）。根据笔者经验，在清洁大气条件下，自行校正结果与官方 L2 产品的 RMSE 通常在 0.01–0.02 之间；在气溶胶负载较高的区域，RMSE 可能达到 0.03–0.05。

九、前沿趋势与范式预判

9.1 云原生与按需计算

遥感预处理正在经历从“本地工作站”向“云原生”的范式迁移。这一迁移的核心驱动力是数据量的爆炸式增长和云计算成本的持续下降。以 Landsat 8/9 为例，全球每年新增约 150 万景影像，单景 L1 产品约 1GB，总数据量超过 1.5PB。在本地处理如此规模的数据已不现实。

云原生预处理的核心特征是：数据存储在对象存储中（如 AWS S3），计算任务以容器化方式部署（如 Docker + Kubernetes），处理流程以工作流引擎编排（如 Apache Airflow）。用户只需定义处理逻辑，平台自动完成资源调度和容错处理。这一模式的优势在于弹性伸缩和按需付费，但也对用户的工程能力提出了更高要求。

9.2 深度学习与物理模型的融合

近年来，深度学习在遥感预处理中的应用日益广泛，尤其是在云检测、大气参数反演和超分辨率重建方面。但纯数据驱动的方法存在可解释性差、泛化能力弱的问题。未来的趋势是物理模型与深度学习的融合：用物理模型生成训练样本，用深度学习加速参数反演，再用物理约束保证结果的合理性。

本文评述认为，这种“物理引导的深度学习”（Physics-Informed Deep Learning）范式有望在未来五年内成为大气校正领域的主流方法。其核心思想是将辐射传输方程作为正则化项嵌入神经网络的损失函数，使网络在拟合数据的同时满足物理约束。这既保留了深度学习的计算效率，又避免了“黑箱”模型的不可解释性。

9.3 标准化与互操作性

随着遥感数据源的多样化（Landsat、Sentinel、MODIS、高分系列），预处理流程的标准化和互操作性成为迫切需求。Open Data Cube、STAC（SpatioTemporal Asset Catalog）等标准的推广，正在推动不同数据源之间的无缝集成。未来，用户可能不再需要关心数据来自哪个传感器，而是通过统一的接口获取经过标准化预处理的地表反射率产品。

十、参考文献与数据来源