简介：这份资源面向人工智能与机器学习方向的学习者，尤其是需要处理遥感影像的初学者与中级开发者，聚焦Landsat8影像的批量预处理问题。内容涵盖数据清洗、缺失值与异常值处理、辐射校正、大气校正、云层去除，以及波段组合、PCA、光谱指数计算等特征工程方法，并借助Python生态中的rasterio、numpy、pandas、scikit-image等库完成自动化批处理，最终输出保留空间信息的TIFF或GeoTIFF数据，可直接用于土地覆盖分类、植被监测、灾害检测等建模任务。资源包共16个文件，以xml配置、py脚本、zbak备份、tif影像及md说明为主，压缩后约46.93MB，目录结构包含源码、测试与数据模块，便于按流程查阅。目前已有151人学习下载，适合希望掌握遥感预处理完整链路、提升模型输入数据质量的读者参考实践。
1. 从一堆原始波段到能用的地表反射率：Landsat8 预处理到底在做什么
手里攒了几十景 Landsat8 的 Level-1 数据，打开一看全是分离的单波段 GeoTIFF，还带着一层薄薄的云和大气散射带来的灰蒙蒙的雾感，直接拿去做 NDVI 或者分类，结果总是不对劲。这不是数据坏了，而是从原始 DN 值到真正能用的地表反射率之间，隔着一整套辐射定标、大气校正、云掩膜和波段堆叠的流程。Landsat8 影像数据预处理方法，说白了就是把卫星传感器记录的原始数字量化值，一步步还原成地表真实反射率的过程，同时把云、云阴影、水体这些干扰因素标出来，最后整理成分析软件能直接读的格式。这套流程单景手动做一遍大概二十分钟，但如果你有几十上百景，手动操作不仅耗时，还容易在参数设置上前后不一致。支持批量处理的意义就在这里：把定标系数、大气模型、气溶胶参数这些固定下来，用脚本一次性跑完整个时间序列，保证每一景的处理逻辑完全一致。适合谁看？做土地利用变化监测的、搞植被长势分析的、需要长时间序列遥感指数的从业者，只要你的工作流里出现过 Landsat8，这套预处理就是绕不开的第一步。
2. 辐射定标与大气校正：从 DN 值到地表反射率的关键两步
2.1 为什么不能跳过辐射定标直接做大气校正
Landsat8 的 Level-1 产品里，像素值是以 16 位无符号整数存储的 DN 值，范围 0 到 65535。这个数值本身没有物理意义，它只是传感器记录的电信号强度。辐射定标的作用，就是利用元数据文件里的乘系数和加系数，把 DN 值转换成传感器处的辐射亮度或者表观反射率。很多人觉得大气校正软件能直接读 DN 值，就跳过定标这一步，结果出来的反射率要么溢出要么偏低，这就是典型的翻车现场。常见做法是用 USGS 提供的定标参数，在元数据 MTL 文件里找到 RADIANCE_MULT_BAND_x 和 RADIANCE_ADD_BAND_x 这两组系数，按波段逐一计算。公式很简单：辐射亮度 = DN 值 × 乘系数 + 加系数。这一步做完，数据就从无量纲的整数变成了有物理单位的辐射亮度值，单位是瓦特每平方米每球面度每微米。
用 Python 的 rasterio 库做批量定标，核心代码大概长这样：
import rasterio import numpy as np import os import glob def radiance_calibration(dn_path, mult, add, out_path): """ 对单波段 DN 影像做辐射定标 dn_path: 输入 DN 影像路径 mult: 乘系数，从 MTL 文件读取 add: 加系数，从 MTL 文件读取 out_path: 输出辐射亮度影像路径 """ with rasterio.open(dn_path) as src: dn = src.read(1).astype(np.float32) profile = src.profile.copy() # 定标公式：辐射亮度 = DN * mult + add radiance = dn * mult + add # 将无效值（通常为 0）设为 NaN radiance[dn == 0] = np.nan profile.update(dtype=rasterio.float32, nodata=np.nan) with rasterio.open(out_path, 'w', **profile) as dst: dst.write(radiance, 1) return out_path # 批量处理示例 mtl_params = { 'B1': {'mult': 0.012353, 'add': -61.765}, 'B2': {'mult': 0.012353, 'add': -61.765}, # 其他波段参数从 MTL 文件解析 } for band, params in mtl_params.items(): dn_file = f'LC08_L1TP_123032_20230101_B{band[1:]}.TIF' radiance_calibration(dn_file, params['mult'], params['add'], f'rad_{band}.TIF')这段代码的逻辑很直接：读入 DN 影像，按波段对应的系数做线性变换，把 0 值区域标记为无效。参数说明上，乘系数的量级通常在 0.01 左右，加系数是负数，用来修正传感器的暗电流偏移。注意不同波段的系数完全不同，必须从 MTL 文件里逐波段读取，不能用一个系数套所有波段。定标完成后，建议抽查几个像素的辐射亮度值，热红外波段应该在 0 到 20 之间，可见光近红外波段在 0 到 500 之间，如果超出这个范围，大概率是系数用错了。
2.2 大气校正的三种路线与参数选择
辐射定标之后，数据是传感器处的辐射亮度，但大气层在太阳和地表之间来回散射和吸收，导致卫星看到的信号里混入了大气程辐射。大气校正的目标就是去掉这部分影响，还原地表反射率。目前主流的路线有三条：基于辐射传输模型的 6S 和 MODTRAN，基于暗像元法的简化校正，以及基于查找表的快速校正。6S 模型精度高，但需要输入大气模式、气溶胶类型、水汽含量等参数，跑一景大概三到五分钟。暗像元法速度快，适合批量处理，但精度依赖影像里有没有干净的暗目标。我一般会这样选：如果做定量遥感，比如反演叶面积指数，用 6S 或者 MODTRAN；如果只是做 NDVI 时间序列，暗像元法够用，而且批量跑起来快得多。
用 Py6S 做大气校正的典型配置：
from Py6S import SixS, AtmosProfile, AeroProfile, Geometry, Wavelength import numpy as np def atmospheric_correction_6s(radiance, solar_zenith, solar_azimuth, view_zenith, view_azimuth, day_of_year): """ 基于 6S 模型计算地表反射率 radiance: 辐射亮度数组 solar_zenith: 太阳天顶角，度 view_zenith: 观测天顶角，度 day_of_year: 年积日，用于确定大气模式 """ s = SixS() # 设置大气模式，中纬度夏季常用 s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer) # 气溶胶类型，大陆型适合内陆地区 s.aero_profile = AeroProfile.PredefinedType(AeroProfile.Continental) # 几何参数 s.geometry = Geometry.User() s.geometry.solar_z = solar_zenith s.geometry.solar_a = solar_azimuth s.geometry.view_z = view_zenith s.geometry.view_a = view_azimuth # 波长设置为对应波段 s.wavelength = Wavelength(0.66) # 红光波段中心波长 s.run() # 获取大气校正系数 xa = s.outputs.coef_xa xb = s.outputs.coef_xb xc = s.outputs.coef_xc # 反射率 = (xa * radiance - xb) / (1 + xc * (xa * radiance - xb)) rho = (xa * radiance - xb) / (1 + xc * (xa * radiance - xb)) return np.clip(rho, 0, 1)这里的关键参数是大气模式和气溶胶类型。大气模式根据纬度和季节选，中纬度夏季、中纬度冬季、热带、极地等，选错了反射率整体会偏。气溶胶类型有大陆型、海洋型、城市型、沙漠型，内陆地区一般用大陆型，城市周边用城市型。几何参数从 MTL 文件里的 SUN_ELEVATION 和 SUN_AZIMUTH 换算，观测天顶角近似为 0 因为 Landsat8 是近星下点观测。波长要对应每个波段的中心波长，蓝光 0.48 微米、绿光 0.56 微米、红光 0.66 微米、近红外 0.86 微米。跑完 6S 会得到三个系数 xa、xb、xc，代入公式就能算出反射率。注意反射率要截断在 0 到 1 之间，负值通常是水汽吸收或者云边缘造成的，直接置零。
2.3 批量处理的目录组织与并行策略
批量处理最容易出问题的地方不是算法，而是文件管理。几十景影像，每景有 11 个波段，加上元数据、质量评估波段，文件数量轻松上千。我习惯按景组织目录，每景一个文件夹，里面放原始波段、中间结果和最终产品。中间结果包括辐射亮度、反射率、云掩膜，最终产品是堆叠后的多波段反射率影像。并行策略上，如果机器有 8 核，可以同时跑 4 景，每景内部按波段串行，这样内存占用可控。用 Python 的 multiprocessing 池做并行，注意每个进程要独立打开文件，不能共享文件句柄。
from multiprocessing import Pool import os def process_scene(scene_dir): """处理单景影像的完整流程""" # 解析 MTL 文件获取定标系数和几何参数 # 逐波段辐射定标 # 逐波段大气校正 # 生成云掩膜 # 波段堆叠 return scene_dir if __name__ == '__main__': scene_dirs = [d for d in os.listdir('raw_data') if d.startswith('LC08')] with Pool(processes=4) as pool: results = pool.map(process_scene, scene_dirs)并行数不要超过物理核心数，否则磁盘 IO 会成为瓶颈。每个进程处理完一景后，把中间文件清理掉，只保留最终产品，不然磁盘很快爆满。如果做时间序列分析，还要注意不同景之间的几何配准，Landsat8 的 Level-1 产品已经做了几何精校正，但不同时相之间还是可能有半个像素的偏移，批量堆叠前最好做一次相对配准。
3. 云掩膜与波段堆叠：让数据真正可用的最后一步
3.1 利用 QA 波段生成云和云阴影掩膜
Landsat8 提供了一个质量评估波段，叫 BQA 或者 QA_PIXEL，里面用比特位存储了每个像素的质量信息。云、云阴影、雪、水体、卷云都有对应的比特位。很多人忽略这个波段，直接用反射率阈值去云，结果要么漏掉薄云，要么把亮地表误判成云。正确做法是解析 QA 波段的比特位，把云和云阴影的像素标记出来。QA_PIXEL 的比特定义是：第 3 位是云，第 4 位是云阴影，第 5 位是雪，第 6 位是晴空，第 7 位是水体，第 8 位是卷云。用位运算提取这些标志。
import rasterio import numpy as np def generate_cloud_mask(qa_path, out_mask_path): """ 从 QA_PIXEL 波段生成云和云阴影掩膜 qa_path: QA 波段路径 out_mask_path: 输出掩膜路径，1 表示云或云阴影，0 表示晴空 """ with rasterio.open(qa_path) as src: qa = src.read(1) profile = src.profile.copy() # 提取云比特位（第 3 位） cloud = (qa >> 3) & 1 # 提取云阴影比特位（第 4 位） cloud_shadow = (qa >> 4) & 1 # 提取卷云比特位（第 8 位） cirrus = (qa >> 8) & 1 # 合并掩膜 mask = (cloud | cloud_shadow | cirrus).astype(np.uint8) profile.update(dtype=rasterio.uint8, nodata=255) with rasterio.open(out_mask_path, 'w', **profile) as dst: dst.write(mask, 1) return out_mask_path位运算的逻辑是右移对应位数再和 1 做与运算，这样就能把单个比特位提取出来。云和云阴影的掩膜合并后，值为 1 的像素在后续分析中应该被剔除或者插值。注意卷云比特位对薄云很敏感，如果做地表反射率反演，建议把卷云也去掉。掩膜生成后，可以用形态学膨胀操作把云的边缘扩大几个像素，因为云的边缘往往有过渡带，直接按比特位剔除可能会留下亮边。
3.2 波段堆叠与无效值处理
预处理完的各个波段是分开的 GeoTIFF，分析软件通常需要一个多波段文件。波段堆叠就是把蓝、绿、红、近红外、短波红外这些波段按顺序写到一个文件里，同时把云掩膜作为单独的一个波段或者外部文件关联。堆叠时要注意波段顺序，常见的排列是 B2、B3、B4、B5、B6、B7，对应蓝、绿、红、近红外、短波红外 1、短波红外 2。无效值处理上，云覆盖的像素在反射率影像里应该设为 NaN，这样后续做统计或者时序分析时能自动忽略。
def stack_bands(band_paths, out_path, cloud_mask_path=None): """ 将多个波段堆叠为一个多波段文件 band_paths: 波段文件路径列表，按顺序排列 out_path: 输出文件路径 cloud_mask_path: 云掩膜路径，可选 """ import rasterio import numpy as np # 读取第一个波段获取元数据 with rasterio.open(band_paths[0]) as src: profile = src.profile.copy() height, width = src.height, src.width # 波段数等于输入波段数，如果有云掩膜则加 1 count = len(band_paths) + (1 if cloud_mask_path else 0) profile.update(count=count, dtype=rasterio.float32, nodata=np.nan) with rasterio.open(out_path, 'w', **profile) as dst: for i, band_path in enumerate(band_paths, start=1): with rasterio.open(band_path) as src: data = src.read(1).astype(np.float32) # 如果有云掩膜，将云像素设为 NaN if cloud_mask_path: with rasterio.open(cloud_mask_path) as mask_src: mask = mask_src.read(1) data[mask == 1] = np.nan dst.write(data, i) # 写入云掩膜作为最后一个波段 if cloud_mask_path: with rasterio.open(cloud_mask_path) as mask_src: dst.write(mask_src.read(1).astype(np.float32), count) return out_path堆叠后的文件可以直接拖进 GIS 软件或者用 Python 读取做分析。注意反射率数据用 float32 存储，文件大小会比原来的 16 位整数大不少，一景大概 500MB 到 1GB。如果做长时间序列，建议按研究区裁剪后再堆叠，能省很多存储空间。云掩膜作为最后一个波段写入，方便后续按波段索引读取。
3.3 批量处理中的元数据一致性检查
批量跑几十景，最怕的是某一景的元数据缺失或者参数异常，导致整批结果里混入一景坏数据。我一般会在处理前做一次元数据扫描，检查每景的 MTL 文件是否完整，定标系数是否存在，太阳高度角是否在合理范围。太阳高度角低于 15 度的影像，大气校正的误差会很大，建议直接标记为低质量。检查脚本可以输出一个 CSV 报告，列出每景的检查结果，人工确认后再批量跑。
import os import csv def check_metadata(scene_dirs, report_path): """检查每景元数据的完整性""" with open(report_path, 'w', newline='') as f: writer = csv.writer(f) writer.writerow(['scene', 'mtl_exists', 'sun_elevation', 'has_calib', 'status']) for scene in scene_dirs: mtl_file = os.path.join(scene, [f for f in os.listdir(scene) if f.endswith('_MTL.txt')][0]) # 解析 MTL 文件，提取太阳高度角和定标系数 # 这里省略具体解析代码 sun_elev = 45.0 # 示例值 has_calib = True status = 'ok' if sun_elev > 15 and has_calib else 'check' writer.writerow([scene, True, sun_elev, has_calib, status])这个检查步骤花不了几分钟，但能避免跑完一整批才发现某景参数不对的尴尬。太阳高度角从 MTL 文件里的 SUN_ELEVATION 字段读取，定标系数检查 RADIANCE_MULT_BAND_x 是否存在。状态标记为 check 的景，人工看一眼再决定要不要跑。
4. 避坑与排查：批量预处理中最容易翻车的五个地方
4.1 定标系数用错波段导致反射率整体偏移
现象：跑完大气校正后，植被的 NDVI 值普遍偏低，水体反射率出现负值。原因：辐射定标时把不同波段的乘系数搞混了，比如用蓝光波段的系数去定标红光波段。Landsat8 的 MTL 文件里，每个波段的系数是分开的，RADIANCE_MULT_BAND_2 对应蓝光，RADIANCE_MULT_BAND_4 对应红光，名字里的数字就是波段号。解决：写脚本时按波段号动态读取系数，不要硬编码。跑完定标后抽查几个像素，蓝光波段的辐射亮度应该比红光波段高，因为蓝光在大气散射中更强。
4.2 大气模式选错导致季节间反射率不可比
现象：同一区域夏季和冬季的影像，做时间序列分析时发现反射率有系统性差异，夏季整体偏高。原因：大气校正时用了固定的大气模式，比如全年都用中纬度夏季，冬季影像的水汽含量低，用夏季模式会高估大气程辐射。解决：根据影像的成像月份和纬度选择大气模式，北半球中纬度地区，5 到 9 月用中纬度夏季，10 到 4 月用中纬度冬季。如果做长时间序列，最好逐月检查大气模式设置。
4.3 云掩膜过于激进把亮地表误删
现象：城市建成区或者沙地的影像，云掩膜把大量像素标记为云，导致有效数据所剩无几。原因：QA 波段的云检测算法对亮地表敏感，水泥地、沙地、雪地的反射率接近云，容易被误判。解决：不要只依赖 QA 波段，结合反射率阈值做二次判断。云在蓝光波段的反射率通常高于 0.2，而亮地表在蓝光波段反射率低于 0.15，可以用这个差异做区分。另外，云掩膜生成后做一次连通域分析，面积小于 10 个像素的云斑块可能是误判，可以剔除。
4.4 并行处理时内存溢出导致进程被杀
现象：批量跑的时候，跑到一半某个进程突然消失，日志里没有报错。原因：并行数设得太多，每个进程都要读入整景影像，内存不够被系统杀掉。一景 Landsat8 的 11 个波段，每个波段约 7000×7000 像素，float32 存储约 200MB，11 个波段就是 2GB 多，如果同时跑 8 个进程，内存需求超过 16GB。解决：根据机器内存设置并行数，8GB 内存的机器最多跑 2 个进程，16GB 跑 4 个。另外，处理完一个波段就释放内存，不要把所有波段都读进来再处理。
4.5 波段堆叠顺序错误导致色彩异常
现象：堆叠后的影像在 GIS 里打开，植被显示为蓝色，水体显示为红色。原因：波段顺序搞反了，把近红外波段放到了蓝光的位置。常见的 RGB 显示顺序是红、绿、蓝，对应 Landsat8 的 B4、B3、B2。如果堆叠时按 B2、B3、B4 的顺序写入，GIS 默认按写入顺序显示，就会把蓝光当红光。解决：堆叠时按分析软件的习惯顺序排列，或者在使用时手动指定波段组合。做 NDVI 的话，近红外和红光的波段号要记牢，B5 是近红外，B4 是红光，NDVI = (B5 - B4) / (B5 + B4)。
5. 用 Python 把整条链路串起来：一个可复现的批量处理脚本框架
5.1 主控脚本的模块划分与参数配置
把前面几章的步骤串成一个完整的批量处理脚本，主控脚本负责遍历景目录、调用各个处理模块、记录日志。参数配置单独放在一个 YAML 文件里，包括大气模式、气溶胶类型、并行数、输出路径这些。这样换一个研究区，只需要改配置文件，不用动代码。
import yaml import os import logging from multiprocessing import Pool def load_config(config_path): """加载 YAML 配置文件""" with open(config_path, 'r', encoding='utf-8') as f: config = yaml.safe_load(f) return config def setup_logging(log_path): """配置日志输出""" logging.basicConfig( level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s', handlers=[ logging.FileHandler(log_path), logging.StreamHandler() ] ) def process_scene(scene_dir, config): """单景处理流程""" try: # 步骤 1：解析 MTL 文件 # 步骤 2：辐射定标 # 步骤 3：大气校正 # 步骤 4：云掩膜生成 # 步骤 5：波段堆叠 logging.info(f'完成处理：{scene_dir}') return {'scene': scene_dir, 'status': 'success'} except Exception as e: logging.error(f'处理失败：{scene_dir}，错误：{str(e)}') return {'scene': scene_dir, 'status': 'failed', 'error': str(e)} if __name__ == '__main__': config = load_config('config.yaml') setup_logging(config['log_path']) scene_dirs = [os.path.join(config['raw_dir'], d) for d in os.listdir(config['raw_dir']) if d.startswith('LC08')] with Pool(processes=config['parallel']) as pool: results = pool.starmap(process_scene, [(s, config) for s in scene_dirs]) # 输出处理报告 success = sum(1 for r in results if r['status'] == 'success') logging.info(f'处理完成：成功 {success} 景，失败 {len(results) - success} 景')配置文件里，大气模式用字符串表示，比如 midlatitude_summer，气溶胶类型用 continental，并行数根据机器内存设。日志同时输出到文件和控制台，方便排查。处理报告记录每景的成功或失败状态，失败的景可以单独重跑。
5.2 处理结果的快速验证方法
批量跑完后，不能直接就用，得做一次快速验证。我一般会抽三景：一景晴空、一景有薄云、一景有厚云，分别检查反射率范围、NDVI 分布和云掩膜覆盖比例。晴空影像的反射率应该在 0 到 0.6 之间，NDVI 植被区在 0.4 到 0.9 之间，云掩膜覆盖比例低于 5%。有薄云的影像，云掩膜覆盖比例在 10% 到 30% 之间，云边缘的反射率过渡自然。厚云影像的云掩膜覆盖比例超过 50%，有效像素主要集中在云间隙。
def validate_result(stack_path, mask_path): """验证单景处理结果""" import rasterio import numpy as np with rasterio.open(stack_path) as src: # 假设波段顺序为 B2, B3, B4, B5, B6, B7, mask red = src.read(3) # B4 nir = src.read(4) # B5 # 计算 NDVI ndvi = (nir - red) / (nir + red + 1e-10) # 统计有效像素的 NDVI 分布 valid_ndvi = ndvi[~np.isnan(ndvi)] print(f'NDVI 均值：{np.mean(valid_ndvi):.3f}，' f'标准差：{np.std(valid_ndvi):.3f}') # 检查反射率范围 print(f'红光反射率范围：{np.nanmin(red):.3f} - {np.nanmax(red):.3f}') with rasterio.open(mask_path) as src: mask = src.read(1) cloud_ratio = np.sum(mask == 1) / mask.size print(f'云掩膜覆盖比例：{cloud_ratio:.2%}')验证脚本输出几个关键指标，和预期范围对比。如果 NDVI 均值低于 0.2，可能是大气校正过度，反射率被压低了。如果云掩膜覆盖比例超过 80%，这景影像基本不可用，直接标记为废弃。验证通过后，把结果按时间序列整理，就可以进入后续分析了。
5.3 时间序列分析前的相对辐射归一化
批量预处理完的影像，虽然每一景都做了绝对辐射校正，但不同时相之间还是可能有微小的辐射差异，来源包括传感器衰减、太阳角度变化、大气校正残差。做时间序列分析前，最好做一次相对辐射归一化，选一景高质量的作为参考，把其他景的反射率调整到和参考景一致。常用方法是直方图匹配或者线性回归，用不变目标（比如深水体、裸土）建立回归关系。
def relative_normalization(target_path, reference_path, out_path): """相对辐射归一化，用线性回归调整目标影像""" import rasterio import numpy as np from sklearn.linear_model import LinearRegression with rasterio.open(reference_path) as ref_src: ref_data = ref_src.read().astype(np.float32) profile = ref_src.profile.copy() with rasterio.open(target_path) as tgt_src: tgt_data = tgt_src.read().astype(np.float32) # 逐波段做线性回归 normalized = np.zeros_like(tgt_data) for i in range(tgt_data.shape[0]): ref_band = ref_data[i].flatten() tgt_band = tgt_data[i].flatten() # 剔除无效值 valid = ~np.isnan(ref_band) & ~np.isnan(tgt_band) if np.sum(valid) < 100: normalized[i] = tgt_data[i] continue model = LinearRegression() model.fit(tgt_band[valid].reshape(-1, 1), ref_band[valid]) normalized[i] = model.predict(tgt_data[i].reshape(-1, 1)).reshape(tgt_data[i].shape) profile.update(dtype=rasterio.float32, nodata=np.nan) with rasterio.open(out_path, 'w', **profile) as dst: dst.write(normalized) return out_path归一化的关键是选好不变目标，深水体的反射率在时间序列上最稳定，裸土次之。回归的斜率接近 1、截距接近 0 说明两景辐射条件接近，如果斜率偏离 1 超过 10%，说明大气校正可能有问题，得回头检查。归一化后的影像再做 NDVI 时间序列，曲线会平滑很多，物候特征也更清晰。
这套流程我跑了三年多，从最开始手动一景一景点，到后来脚本批量跑，最大的教训就是参数一定要统一。有一次偷懒，大气模式用了默认值没改，结果半年的影像反射率整体偏高，NDVI 曲线完全看不出物候规律，白白浪费了一周时间重跑。现在我的习惯是，每批数据跑之前，先拿一景做测试，确认参数无误再批量。希望帮到你。
本文还有配套的精品资源，点击获取
