Landsat 8批量预处理全流程：从MTL解析到TOA反射率与云掩膜

简介：陆地卫星8号（Landsat8）影像的批量预处理是遥感数据建模的关键环节。这套资源面向机器学习与数据预处理领域的技术人员，针对影像去云、辐射校正、大气校正、波段组合和光谱指数计算等常见需求，梳理出一套可复用的Python处理流程，重点解决批处理效率与特征工程质量问题。压缩包共16个文件，大小约46.93MB，除核心Python脚本（PreprocessL8.py、base.py）外，还包含项目配置文件（xml/iml/json）、数据备份（zbak）、示例GeoTIFF以及README说明，整体结构贴近真实工程目录。目前已有148人学习浏览。读者可以从可直接运行的脚本、附赠内容与说明文档中，快速掌握从影像加载、云遮挡处理、特征构建到尺度归一化的完整操作思路，便于在此基础上扩展土地覆盖分类、植被监测等机器学习应用，也可根据自身数据替换参数并二次开发。

1. 为什么说 Landsat 8 预处理的第一步不是写代码

Landsat 8 的 L1TP 产品不是"打开就能出图"的数据。它存储的是 DN 值，没有经过大气校正，也没有统一量纲；不同场景、波段的 DN 范围差异很大，直接拿去算 NDVI，会同时受到太阳高度角、大气散射和传感器增益干扰。预处理就是把 DN 还原成有物理意义的反射率，把云和云影挑出来，再裁剪到目标区域。这个流程单看每一步都不难，难在要重复几十上百次：一景压缩包解压后有十几个波段，时序分析常常一次跑二三十景，手动在 ENVI 里点一轮要十几分钟。所以这个标题真正要解决的是让预处理变得可重复、可并行、可断点续跑。适合已经手动跑通一景、想把它产品化的从业者。

2. 读懂 Landsat 8 L1TP 结构与 MTL 元数据，再写批量预处理脚本

2.1 L1TP 压缩包里装了哪些文件，预处理用到哪几个

Landsat 8 最常见的下载形态是 Collection 2 L1TP 的.tar.gz压缩包，解压后包含 MTL 元数据文件、十多个波段 TIF、角度文件和 QA 文件。L1TP 已经做了系统几何校正，所以整个预处理流程不涉及几何纠正，重心全部放在辐射、大气和掩膜上。批处理脚本真正要读的只有 MTL 和下面这张表里的波段，其余文件可以不解压，按文件名模式定向读取：

B8 全色波段是 15 m，通常不参与多光谱预处理，只在融合阶段用；B1 在陆表应用里噪声偏大，我一般不会默认纳入批量输出。热红外波段的定标系数和其他波段不同，批量脚本里要单独处理。QA_PIXEL 是 Collection 2 的命名，Collection 1 里叫 BQA，位定义也不一样，写脚本前先确认数据版本。

2.2 MTL 文件中决定预处理精度的五类参数

预处理脚本从 MTL 里读字段，但读哪些字段决定了精度和通用性。我通常把这些字段提取出来，作为每个场景的元数据字典：

常见的误用是把不同波段的系数搞混。每个波段在 MTL 里是独立一行，从 RADIANCE_MULT_BAND_1 到 RADIANCE_MULT_BAND_11，脚本不能用"取第一个"或"全部相同"的写法。另一个容易漏的是 SUN_ELEVATION：Collection 2 的 TOA 反射率系数公式是 ρ = (Mρ × DN + Aρ) / sin(SUN_ELEVATION)，漏掉这一步，中纬度地区反射率会整体偏高 15%–30%，叠进时序分析就是系统性偏差。所以批量脚本的第一步不是读影像，而是扫目录、读 MTL、列清单，把一景场景当作一个自包含的任务单元。

2.3 解析 MTL 并生成批量任务清单的最小脚本

MTL 正文是形如KEY = VALUE的键值对，值可能是数字或带引号的字符串。用正则解析比引入 XML 解析器更省事，代码也更容易让别人接手：

import re from pathlib import Path def parse_mtl(mtl_path: Path) -> dict: text = mtl_path.read_text(encoding="utf-8", errors="ignore") fields = {} for key, val in re.findall(r'(\w+)\s*=\s*"?([^"\n]*)"?', text): fields[key.strip()] = val.strip() return fields scene_root = Path("L8_raw") tasks = [] for tar_dir in sorted(scene_root.glob("LC08*")): mtl = next(tar_dir.glob("*MTL.txt")) meta = parse_mtl(mtl) if "REFLECTANCE_MULT_BAND_2" not in meta: print(f"跳过 {tar_dir.name}：缺少反射率系数，可能是旧版本产品") continue tasks.append({ "scene_id": meta.get("PRODUCT_ID"), "mtl": mtl, "dir": tar_dir, **{k: meta[k] for k in ["SUN_ELEVATION", "EARTH_SUN_DISTANCE"] if k in meta} }) print(f"有效任务数：{len(tasks)}")

正则里(\w+)匹配键名，"?([^"\n]*)"?匹配等号右侧的值并去掉外层引号，遇到注释行也能跳过。把 MTL 解析结果和场景路径绑成一条任务记录，而不是每步处理时重新解析文件，好处有三点：单景失败时按 scene_id 就能单独重跑；预处理前就能筛掉损坏或版本不符的场景，避免跑到第 20 景才中断；后续要加字段（比如云量统计）只改一处解析逻辑。任务清单还可以落成 CSV 人工检查。

3. 用 Python 批量完成 Landsat 8 辐射定标与 TOA 反射率计算

3.1 辐射定标的两个公式和一条分界线

辐射定标解决的是把传感器记录的 DN 变成物理量。Landsat 8 需要两个层次的物理量，公式不一样：

• 辐亮度，热红外和定量遥感用：Lλ = RADIANCE_MULT_BAND_x × DN + RADIANCE_ADD_BAND_x
• 表观反射率 TOA，多光谱波段最常用：ρλ = (REFLECTANCE_MULT_BAND_x × DN + REFLECTANCE_ADD_BAND_x) / sin(SUN_ELEVATION)

分界线在于：做陆表分类、植被指数的，直接用 TOA 反射率就够；做地表温度、水体定量反演的，才需要先走辐亮度再进大气传输模型。很多教程把两套系数混着用，最常见的问题是手动用辐亮度公式算完，又去乘 cos(太阳天顶角)，等于把太阳角校正做两遍。用 MTL 提供的反射率系数时，日地距离已经由 USGS 在定标参数里折算进去了，脚本里不要再乘 EARTH_SUN_DISTANCE。对 NDVI 这类比值型指数，大气影响在红和近红外的贡献方向相反，TOA 反射率算出的 NDVI 与地表 NDVI 相关性很强，因此快速监测流程通常直接拿 TOA 当中间产物。

3.2 rasterio 加 numpy 的批量 TOA 反射率实现

多光谱预处理只取 B2–B7 六个波段，覆盖了绝大多数陆表应用。实现时用 rasterio 读取原始 DN，numpy 做数组运算，最后统一写成一个多波段 GeoTIFF：

import numpy as np import rasterio from pathlib import Path BANDS = ["B2", "B3", "B4", "B5", "B6", "B7"] OUT_DIR = Path("L8_toa") OUT_DIR.mkdir(exist_ok=True) def toa_reflectance(task: dict) -> Path: meta = parse_mtl(task["mtl"]) sin_se = np.sin(np.radians(float(meta["SUN_ELEVATION"]))) # 用 B2 的 profile 作为输出模板，保证投影、分辨率、原点一致 with rasterio.open(next(task["dir"].glob("*_B2.TIF"))) as ref: profile = ref.profile arrays = [] for band in BANDS: key_m = f"REFLECTANCE_MULT_BAND_{band[1:]}" key_a = f"REFLECTANCE_ADD_BAND_{band[1:]}" mult = float(meta[key_m]) add = float(meta[key_a]) with rasterio.open(next(task["dir"].glob(f"*_{band}.TIF"))) as src: dn = src.read(1).astype("float32") arrays.append((mult * dn + add) / sin_se) stacked = np.stack(arrays) profile.update(count=len(BANDS), dtype="float32", compress="deflate", nodata=-9999.0) out_path = OUT_DIR / f"{task['scene_id']}_TOA.tif" with rasterio.open(out_path, "w", **profile) as dst: dst.write(stacked) return out_path

(mult * dn + add) / sin_se这个顺序不要拆开写，多光谱系数和太阳角校正必须在同一个反射率域内完成。用float32而不是uint16，因为反射率是 0–1 的小数，整型会把小数截断；单波段 float32 在 30 m 场景下约 240 MB，六波段用 DEFLATE 压缩后体积在 400–500 MB，属于可接受范围。nodata=-9999.0在写 GeoTIFF 时同时写入元数据，后续裁剪和镶嵌不需要另外约定。输出波段顺序固定为 B2–B7，这是和下游算法（NDVI、分类器特征矩阵）的约定，不要跟着 MTL 里的顺序走。

3.3 批量脚本最容易出错的三个位置

多景出错往往不是算法问题，而是这三处。

第一，太阳高度角单位。MTL 里 SUN_ELEVATION 单位是度，np.sin默认接收弧度，漏写np.radians会算出完全不可用的反射率。这个错非常隐蔽，单景肉眼不一定看得出来，批量均值对比时会发现整体偏大。

第二，热红外波段混入多光谱栈。B10、B11 没有 REFLECTANCE 系数，一旦循环写成"遍历所有 B 开头的文件"，跑到热红外就会 KeyError。处理办法是把热红外和卷云波段从多光谱栈中显式排除，单独走辐射定标分支。

第三，内存峰值。单个 float32 波段约 240 MB，六波段叠加 1.4 GB，再叠加大气校正中间数组，单场景内存需求接近 4 GB。并行时要按这个基准倒推 worker 数量，不要上来就填满 CPU 核数。

4. 大气校正、QA 云掩膜与矢量裁剪的批量衔接

4.1 大气校正选型：批量场景下先想清楚要不要自己做

大气校正把 TOA 反射率变成地表反射率。完整的辐射传输模型（FLAASH、6S）精度高，但批量场景下代价很大。下面是常见做法的比较：

提示：如果项目最终只是算 NDVI、做土地利用分类，直接下载 Collection 2 L2 的 SR 产品，连大气校正都不用自己做；只有只能拿到 L1 时才需要自行校正。

如果只能在本地批量处理 L1，DOS 暗目标法是成本最低的路，简化实现甚至只有一行：对每个波段做 1% 分位截断，再整体夹逼到 0。

sr_band = np.maximum(toa_band - np.percentile(toa_band, 1), 0.0)

它假设波段最暗像元基本是大气程辐射贡献的，减去这个分位值就近似去掉了气溶胶影响。局限在于暗目标区域必须真实存在，全图都是高反射地物时会扣过头，出现负值，所以截断到 0 是必要的。不要把 DOS 的结果当成精确地表反射率，它只够支撑分类和时序趋势分析。

4.2 用 QA_PIXEL 位掩膜批量生成云、云影与水体

Collection 2 的 QA_PIXEL 每个像元是一个 16 位整数，用位标记不同属性。批量脚本最关心下面几位：

批量解码脚本：

def decode_qa_pixel(qa_path: Path, out_mask: Path): with rasterio.open(qa_path) as src: qa = src.read(1) profile = src.profile cloud = ((qa & (1 << 3)) > 0) | ((qa & (1 << 2)) > 0) dilated = (qa & (1 << 1)) > 0 shadow = (qa & (1 << 4)) > 0 fill = (qa & (1 << 0)) > 0 # 坏像元 = 云 + 膨胀云 + 云影，剔除填充区 bad = (cloud | dilated | shadow) & ~fill profile.update(dtype="uint8", count=1, compress="deflate", nodata=255) with rasterio.open(out_mask, "w", **profile) as dst: dst.write(bad.astype("uint8"), 1) return bad

位与运算qa & (1 << n)判断第 n 位是否为 1，|把云、膨胀云、云影合并成一个坏像元掩膜。掩膜输出用 uint8 而不是布尔类型，是为了兼容下游的分类器和统计工具。这里有个容易踩的坑：Collection 1 的 BQA 位定义和 Collection 2 的 QA_PIXEL 不同，位 3 和位 4 的含义可能颠倒，脚本要根据 PRODUCT_ID 里的 C1/C2 标识分支处理。

4.3 按矢量边界批量裁剪并把掩膜同步裁出来

预处理流水线的最后一步是把 TOA 反射率和坏像元掩膜裁到 AOI 范围，减少后续存储和计算量。用 rasterio.mask 的 crop 参数可以同时拿到裁剪数组和变换矩阵：

import geopandas as gpd from rasterio.mask import mask as rio_mask import rasterio def crop_scene(toa_path: Path, qa_path: Path, aoi_path: Path, out_dir: Path): aoi = gpd.read_file(aoi_path) with rasterio.open(toa_path) as src: aoi_geom = aoi.to_crs(src.crs).geometry arr, tf = rio_mask(src, aoi_geom, crop=True, nodata=-9999.0) prof = src.profile.copy() prof.update(height=arr.shape[1], width=arr.shape[2], transform=tf) out_toa = out_dir / f"{toa_path.stem}_crop.tif" with rasterio.open(out_toa, "w", **prof) as dst: dst.write(arr) # 掩膜用同一套 AOI 裁剪，保证像元位置对齐 with rasterio.open(qa_path) as src: aoi_geom = aoi.to_crs(src.crs).geometry mask_arr, mask_tf = rio_mask(src, aoi_geom, crop=True, nodata=255) prof = src.profile.copy() prof.update(height=mask_arr.shape[1], width=mask_arr.shape[2], transform=mask_tf) out_qa = out_dir / f"{qa_path.stem}_crop.tif" with rasterio.open(out_qa, "w", **prof) as dst: dst.write(mask_arr) return out_toa, out_qa

注意三点。第一，AOI 矢量必须先转到影像的 CRS 再传给 rio_mask，直接用 WGS84 的 shp 去裁 UTM 影像会得到空结果。第二，TOA 和 QA 必须用同一个aoi_geom裁剪，否则两幅输出像元错位，后面统计坏像元比例时会出问题。第三，多景覆盖同一 AOI 时这个环节不做镶嵌，保留场景粒度，到做时序合成时再按日期和云量选像元，这样每条处理链路都保持独立可重跑。

5. Landsat 8 预处理产物的三个数值校验与多进程提速

5.1 批量完成后先校验三个数值

批量预处理最怕的不是报错，而是不报错但数值错了。每批跑完我会先做三个校验：

1. 反射率分布范围。正常植被场景 B5 中位数在 0.2–0.4，1% 分位接近 0，99% 分位不超过 1.2。如果中位数小于 0.05 或最大值大于 2，先怀疑太阳高度角 sin() 那一步。
2. 坏像元比例和 MTL 云量对得上。统计 QA 掩膜坏像元比例，和 MTL 里 CLOUD_COVER 字段对比，偏差超过 20% 说明 QA 解码位顺序用错了。
3. 相邻场景重叠区一致性。两景重叠条带各自算中位数，差值大于 0.03 时检查是否某景大气校正参数异常。

def quick_check(toa_path: Path, qa_path: Path): with rasterio.open(toa_path) as src: b5 = src.read(src.indexes[3]) # B2-B7 中 B5 是第 4 个波段 with rasterio.open(qa_path) as src: bad = src.read(1) bad_ratio = float((bad > 0).mean()) percentiles = np.nanpercentile(b5, [1, 50, 99]) return {"b5_pct": percentiles, "bad_ratio": bad_ratio}

5.2 用进程池提速，并按内存倒推并发数

单景六波段 float32 栈约 1.4 GB，DOS 和掩膜处理的中间数组会再翻倍，每个 worker 预留 3–4 GB 比较稳妥。16 GB 内存的机器 max_workers=4 是安全值；核数再多，瓶颈也在内存不在 CPU。进程级并行而不是线程级并行，因为 rasterio 底层读写在多线程下会争用文件句柄，而且场景之间天然独立：

from concurrent.futures import ProcessPoolExecutor, as_completed def process_one(task): toa = toa_reflectance(task) qa_path = next(task["dir"].glob("*QA_PIXEL.TIF")) mask_path = OUT_DIR / f"{task['scene_id']}_bad.tif" decode_qa_pixel(qa_path, mask_path) crop_scene(toa, mask_path, AOI_PATH, OUT_DIR) return task["scene_id"] with ProcessPoolExecutor(max_workers=4) as pool: futures = [pool.submit(process_one, t) for t in tasks] for fut in as_completed(futures): print("完成:", fut.result())

Windows 下要把入口代码放进if __name__ == "__main__":，否则进程池会递归导入主模块导致死锁。跑批时我会把每个场景状态写进 log.csv，每行包含 scene_id、状态、异常信息、B5 中位数和坏像元比例；中断后只把状态为 fail 的行转成 task 子集重跑，断点续跑比整批重来省的时间往往比预处理本身还多。

本文还有配套的精品资源，点击获取