简介陆地卫星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 和下面这张表里的波段其余文件可以不解压按文件名模式定向读取文件后缀波段区间 (µm)分辨率预处理用途B1 海岸气溶胶0.43–0.4530 m水色、气溶胶反演B2–B4 蓝绿红0.45–0.6830 m真彩色合成、NDVI 辅助波谱B5 近红外0.85–0.8830 mNDVI 核心波段B6–B7 短波红外1.56–2.3030 m地表含水量、云/雪识别B9 卷云1.36–1.3930 m卷云掩膜参考B10–B11 热红外10.6–12.5100 m地表温度只用辐射定标QA_PIXEL—30 m云、云影、雪、水的位掩膜B8 全色波段是 15 m通常不参与多光谱预处理只在融合阶段用B1 在陆表应用里噪声偏大我一般不会默认纳入批量输出。热红外波段的定标系数和其他波段不同批量脚本里要单独处理。QA_PIXEL 是 Collection 2 的命名Collection 1 里叫 BQA位定义也不一样写脚本前先确认数据版本。2.2 MTL 文件中决定预处理精度的五类参数预处理脚本从 MTL 里读字段但读哪些字段决定了精度和通用性。我通常把这些字段提取出来作为每个场景的元数据字典字段用途注意事项RADIANCE_MULT_BAND_x / RADIANCE_ADD_BAND_x辐射定标DN 转辐亮度x 按波段号热红外波段要用REFLECTANCE_MULT_BAND_x / REFLECTANCE_ADD_BAND_xTOA 反射率系数只对 B1–B9 存在SUN_ELEVATION太阳高度角别忘了 sin()EARTH_SUN_DISTANCE日地距离用系数法时已包含无需再用PRODUCT_ID / LANDSAT_SCENE_ID结果命名与去重批量日志主键常见的误用是把不同波段的系数搞混。每个波段在 MTL 里是独立一行从 RADIANCE_MULT_BAND_1 到 RADIANCE_MULT_BAND_11脚本不能用取第一个或全部相同的写法。另一个容易漏的是 SUN_ELEVATIONCollection 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(encodingutf-8, errorsignore) 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 读取原始 DNnumpy 做数组运算最后统一写成一个多波段 GeoTIFFimport 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_okTrue) 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 fREFLECTANCE_MULT_BAND_{band[1:]} key_a fREFLECTANCE_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(countlen(BANDS), dtypefloat32, compressdeflate, 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精度高但批量场景下代价很大。下面是常见做法的比较方案精度批量成本适用场景ENVI FLAASH高每景要人工检查气溶胶和水汽参数脚本化麻烦关键景、论文级处理6S 独立程序高需要逐景准备大气参数集成成本高定量反演研究DOS 暗目标法中纯 numpy 可批量无需外部数据地物分类、NDVI 时序直接用 L2 SR 产品官方下载时选择 Surface Reflectance零处理量大多数应用首选提示如果项目最终只是算 NDVI、做土地利用分类直接下载 Collection 2 L2 的 SR 产品连大气校正都不用自己做只有只能拿到 L1 时才需要自行校正。如果只能在本地批量处理 L1DOS 暗目标法是成本最低的路简化实现甚至只有一行对每个波段做 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 位整数用位标记不同属性。批量脚本最关心下面几位位含义脚本用途bit 0填充判定有效像元时剔除bit 1膨胀云云周围缓冲通常并入云掩膜bit 2卷云薄云检测bit 3云核心云检测bit 4云影与云合并作为坏像元bit 5雪高反射目标单独掩膜bit 6晴空有效像元bit 7水水体掩膜批量解码脚本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(dtypeuint8, count1, compressdeflate, nodata255) 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, cropTrue, nodata-9999.0) prof src.profile.copy() prof.update(heightarr.shape[1], widtharr.shape[2], transformtf) 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, cropTrue, nodata255) prof src.profile.copy() prof.update(heightmask_arr.shape[1], widthmask_arr.shape[2], transformmask_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 批量完成后先校验三个数值批量预处理最怕的不是报错而是不报错但数值错了。每批跑完我会先做三个校验反射率分布范围。正常植被场景 B5 中位数在 0.2–0.41% 分位接近 099% 分位不超过 1.2。如果中位数小于 0.05 或最大值大于 2先怀疑太阳高度角 sin() 那一步。坏像元比例和 MTL 云量对得上。统计 QA 掩膜坏像元比例和 MTL 里 CLOUD_COVER 字段对比偏差超过 20% 说明 QA 解码位顺序用错了。相邻场景重叠区一致性。两景重叠条带各自算中位数差值大于 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 GBDOS 和掩膜处理的中间数组会再翻倍每个 worker 预留 3–4 GB 比较稳妥。16 GB 内存的机器 max_workers4 是安全值核数再多瓶颈也在内存不在 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_workers4) 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 子集重跑断点续跑比整批重来省的时间往往比预处理本身还多。本文还有配套的精品资源点击获取