简介：面向人工智能与机器学习领域的遥感数据预处理需求，资源聚焦Landsat8影像的批量处理，适合从事环境监测、农业分析、城市规划等方向的数据科学与GIS学习者，也适合需要提升特征工程能力的Python开发者。资源内含完整的Python预处理脚本与示例数据，覆盖数据清洗、缺失值处理、异常值检测、云遮挡处理、辐射校正与大气校正等关键环节，并演示了波段组合、PCA、光谱指数计算等特征工程方法，为后续建模提供高质量输入。压缩包共16个文件，以py脚本、xml配置、zip附件、tif示例数据及md说明文档为主，整体大小约46.93MB，目录结构清晰，便于按需取用。已有149人学习。通过阅读README并运行脚本，可掌握基于rasterio、geopandas、numpy、pandas等库的批处理流程，包括图像加载、波段校正、特征创建与数据统一保存；同时参考项目中的工程配置，能快速迁移至其他遥感任务，提升机器学习模型预测准确性。
1. Landsat8影像数据预处理：把批量当成一条流水线，而不是一次循环
第一次带二十多景Landsat8影像做批量预处理时，我以为工作量最大的是辐射定标和大气校正的参数选取。做完才发现，真正的门槛是想让每一景在同一条流水线上稳定跑通：元数据读对、无效值掩住、投影对齐、云被标出来，任何一步在一景上翻车，整个批次都要返工。这也是Landsat8影像数据预处理最容易被低估的地方——单景可以慢慢调，批量处理逼着你把过程拆成可复现的模块。这篇笔记按一线工程师落地的方式讲：先拆流程，再给可抄的批处理脚本、参数表和避坑记录，适合准备做时序分析、土地覆盖制图或目视解译底图的人。
2. 预处理流程怎么设计：从L1级产品到可分析影像
2.1 五个标准环节，缺一个后面都要补账
Landsat8 L1级产品（Collection 2 L1TP）已经做过辐射校正，并嵌入了正射信息，但开发者拿到手的数据仍然是DN值，不能直接用于多时相对比。一个能支撑后续分类或时序分析的预处理流程，至少要包含五个环节：元数据核对与数据清洗、辐射定标、大气校正、几何对齐与重采样、云遮罩与无效值处理。前两步属于图像预处理里最容易被一带而过的部分，后两步决定数据能不能在不同日期、不同条带之间直接叠加。
这里要特别提醒Collection版本问题。市面上一大半老教程还是按Collection 1的ESUN公式在写辐射定标，Collection 2改掉了这一套，直接在MTL里给REFLECTANCE_MULT和REFLECTANCE_ADD。如果拿着老脚本去跑新数据，输出的TOA反射率会整体偏移，差值在0.05上下，本来能用的数据直接被污染。所以第一步元数据核对不只是看云量，还应该确认数据版本，否则后面的辐射定标、大气校正全建立在一个错误前提上。
流程里的“数据清洗”也是一步容易被跳过的工程。它做的事很像文本表格里的清洗：把坏行、重复条带、拼错的轨道号、格式不一致的数据段拦在正式处理之前。Landsat8单景数据量不大，手动看几十个MTL还能忍，但一旦进入上百景的批处理，文件名和元数据的错位会让后处理全部白跑。因此批量场景下，数据清洗必须在循环最前面独立成模块。
2.2 批量处理的两条路线：ENVI批处理和脚本化
Landsat8批量预处理的主流落地方式有两条。第一条是ENVI里搭模型：用Radiometric Calibration做辐射定标，接QUAC或FLAASH做大气校正，最后Layer Stacking和几何校正一起挂在流程里，用ENVI的批处理工具把几十个MTL交给同一个模型去跑。这条路的好处是界面友好、参数有下拉框，项目组里有人不熟代码也能接手。坏处是流程相对固定，中途一景失败往往要整批重来，日志太粗，定位问题全靠肉眼。
第二条是Python加GDAL的脚本化路线。GDAL负责读影像、写影像、做gdalwarp重投影，MTL.txt这种文本格式用简单正则就能全部解析出来，循环遍历目录就是“批量”。这条路适合几十景以上、要做可复现处理、之后要接机器学习或时序分析的人。命令行日志能精确到“第几景第几个波段哪一步报错”，加上每景独立的try/except，失败一景不会拖垮整批。
我一般推荐脚本化，并且把每一步设计成独立函数再加一个总入口，而不是一个大循环从头写到尾。下面是两条路线的对比，方便按自己的团队情况选。
| 对比项 | ENVI批处理/模型 | Python+GDAL脚本 |
|---|---|---|
| 适合规模 | 单次几十景，作业出图 | 时序、跨年份、大规模反复处理 |
| 学习成本 | 低，界面操作即可 | 中，需要会Python和命令行 |
| 断点恢复 | 较差，一景失败常要整批重跑 | 每景独立捕获异常，记录fail清单 |
| 扩展性 | 依赖ENVI许可与版本 | 开源，便于长期复用和团队共用 |
| 常见翻车点 | QUAC与工程格式不兼容 | 不同UTM分带没统一、坐标单位混乱 |
2.3 数据清洗：批量前的元数据校验
“数据预处理之数据清洗”放在Landsat8批处理里，最典型的工作是校验三类信息：产品ID与文件名是否一致、云量是否超阈值、数据级别是不是L1TP。因为Landsat8 L1TP里偶尔混进个别L1GT（几何粗校正），这类数据做不了同级别正射产品的对齐，混在一起处理会让镶嵌结果出现几百米的偏移。
下面这段对全目录MTL做快速扫描，可以在进入正式预处理前把问题景筛出来。
for mtl in $(ls ./L8/*MTL.txt); do id=$(grep 'LANDSAT_PRODUCT_ID' "$mtl" | awk '{print $2}' | tr -d '"') cloud=$(grep 'CLOUD_COVER' "$mtl" | awk '{print $2}') type=$(grep 'DATA_TYPE' "$mtl" | awk '{print $2}' | tr -d '"') echo "$id | cloud=$cloud | type=$type" done这段代码不依赖任何地理处理库，只从MTL文件里抓三个关键字段，然后排成一行输出。grep把键对应的值抓到，awk取第二个字段，tr去掉引号，最终形成“产品ID | 云量 | 数据类型”的清单。实际项目里会把结果重定向到csv文件，再按云量大于20%、数据类型不是L1TP两个条件过滤，被过滤掉的景直接移出处理目录，后面的大循环就不会再在脏数据上浪费时间。
另一个容易忽略的清洗点：同一场景目录里如果混入了其他传感器的同名文件，比如高分三号预处理输出的辅助栅格，脚本的globbing会把格式错误的文件带进Landsat8处理流。因此2.3节的扫描脚本也建议把文件扩展名、格式、波段数打出来，确认全部输入都是同一规格再做辐射定标。
3. 用Python+GDAL批量做辐射定标：MTL解析与TOA反射率计算
3.1 辐射定标：为什么不能只用DN值
Landsat8 L1级产品的DN值是数字量化值，不同日期、不同太阳高度角下拍出来的同一地表，DN值可能差出一大截。辐射定标就是把DN换算成大气顶（TOA）反射率，消除太阳位置和日地距离的影响，让不同时相的影像具有可比性。
Collection 2 L1TP的MTL文件里已经直接给出了反射率增益和偏置，即REFLECTANCE_MULT_BAND_x与REFLECTANCE_ADD_BAND_x。换算公式很简单：
ρ = (REFLECTANCE_MULT × DN + REFLECTANCE_ADD) / sin(SUN_ELEVATION)
其中SUN_ELEVATION是MTL里的太阳高度角。老教程里的ESUN公式在Collection 2里已经不需要了，因为美国地质调查局把这套系数直接预置在MTL里，再用老公式反而是画蛇添足，输出的结果还会因为单位换算错误漂移。
需要明白的是：辐射定标只做到TOA反射率，它去掉的是传感器响应和太阳几何作用，没有去掉大气散射吸收的影响。TOA图像里仍然有一层“大气滤镜”，肉眼看起来会偏白偏亮，这就是下一步要交给大气校正处理的。
3.2 可复用脚本：按MTL批量循环
下面是一个直接可跑的批量辐射定标脚本。它遍历目录下所有MTL.txt，逐个解析系数，把每个波段写成TOA反射率GeoTIFF。
# -*- coding: utf-8 -*- """ Landsat8 Collection 2 L1TP 批量辐射定标：DN -> TOA反射率 依赖：gdal>=3.0, numpy """ import glob import os import numpy as np from osgeo import gdal def parse_mtl(mtl_path): """把MTL.txt每一行'KEY = VALUE'解析成dict""" kv = {} with open(mtl_path, 'r', encoding='utf-8', errors='ignore') as f: for line in f: if '=' in line: k, v = [x.strip() for x in line.split('=', 1)] kv[k] = v.strip('"') return kv def toa_refl(dn, mult, add, sun_elev): """单波段DN转TOA反射率，0-1裁剪""" rad = dn.astype(np.float64) * mult + add refl = rad / np.sin(np.radians(sun_elev)) return np.clip(refl, 0, 1) gdal.UseExceptions() mtl_list = sorted(glob.glob('./L8/*MTL.txt')) failed = [] for mtl in mtl_list: try: meta = parse_mtl(mtl) scene = meta['LANDSAT_PRODUCT_ID'] out_dir = f'./output/{scene}_toa' os.makedirs(out_dir, exist_ok=True) sun_elev = float(meta['SUN_ELEVATION']) prefix = os.path.basename(mtl).replace('_MTL.txt', '') for band in range(1, 12): # L1TP波段号1-11 m_key = f'REFLECTANCE_MULT_BAND_{band}' a_key = f'REFLECTANCE_ADD_BAND_{band}' if m_key not in meta or a_key not in meta: continue # 热红外等波段没有反射率系数，安全跳过 dn_path = f'./L8/{prefix}_B{band}.TIF' if not os.path.exists(dn_path): continue ds = gdal.Open(dn_path, gdal.GA_ReadOnly) band_ds = ds.GetRasterBand(1) dn = band_ds.ReadAsArray().astype(np.float64) mult = float(meta[m_key]) add = float(meta[a_key]) refl = toa_refl(dn, mult, add, sun_elev) drv = gdal.GetDriverByName('GTiff') out_ds = drv.Create(out_dir + f'/{scene}_B{band}_toa.tif', ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_band = out_ds.GetRasterBand(1) out_band.WriteArray(refl) out_band.SetNoDataValue(-9999) out_ds.FlushCache() out_ds = None ds = None print(f'[OK] {scene} band {band} ' f'refl range=({refl.min():.3f}, {refl.max():.3f})') except Exception as e: failed.append((mtl, str(e))) print(f'[FAIL] {mtl}: {e}') if failed: print('以下文件处理失败：') for f in failed: print(f)这段脚本的逻辑分三层。parse_mtl负责把MTL文本变成字典；toa_refl实现反射率公式；主循环对每个MTL、每个波段执行“读取DN矩阵 → 套系数 → 写浮点GeoTIFF”。外层套了try/except，失败信息全部收进failed列表，这样某一景某个波段读写异常时，后面几十景继续照跑，不会整批中断。
参数说明里有两个细节值得留意。第一，波段循环从1到11，但只有含反射率系数的波段才会被处理，热红外波段没配REFLECTANCE_MULT就自然跳过，不需要单独维护一个“光学波段白名单”。第二，输出裁剪到0到1之间，是为了方便保存成Float32并直接做快视图；如果后面要严格识别云影，建议把np.clip去掉，因为负反射率本身是异常像元或大气残留的信号，应该在统计时保留而不是提前抹掉。
批量跑完后，fail清单就是唯一的复盘列表。我一般把它导出成文本，连同每一景的输出范围记录一起存档。下次再跑同样流程，直接对比range列表就能发现那一景的系数是不是被MTL改动影响了。
3.3 参数表：Landsat8 MTL里会用到的关键系数
| 参数 | 含义 | 批量处理时关注点 |
|---|---|---|
| REFLECTANCE_MULT_BAND_x | 反射率增益 | Collection 2下直接使用，勿再套老ESUN公式 |
| REFLECTANCE_ADD_BAND_x | 反射率偏置 | 通常是负值，公式里直接相加即可 |
| SUN_ELEVATION | 太阳高度角（度） | 先转弧度，再取sin做分母 |
| EARTH_SUN_DISTANCE | 日地距离（天文单位） | 只在用老辐射亮度公式时才需要平方运算 |
| CLOUD_COVER | 整景云量百分比 | 批量筛选时的硬阈值，不能替代像素级云掩膜 |
这张表读懂后，真正的批量调参工作就变成了“确认所有MTL都在同一Collection版本”。只要有少数几景是Collection 1，跑出来的反射率就会在数值上和同批数据有系统性差异。判断方法很简单：打开MTL看有没有REFLECTANCE_MULT_BAND_2，有就是Collection 2；没有则是老版本，需要单独处理或重新下载。
4. 大气校正与几何校正：精度瓶颈在哪
4.1 三种可行大气校正选型
辐射定标后的TOA反射率，仍然带着大气的“滤镜”。大气校正就是把TOA反射率进一步还原为地表反射率，它是整个Landsat8预处理流程里最容易拉开结果差异的环节。选哪个模型，取决于手上有什么辅助数据、批次规模多大、结果要做什么用途。常见的有三种：
- 6S模型：理论最完整，逐波段模拟大气辐射传输，但参数多，需要气溶胶光学厚度、水汽含量、大气模型类型等输入。批量几十景时，每景都要配套再分析气象数据，参数不齐就报错，适合科研级关键日期的精细校正。
- COST模型：属于DOS暗目标法的简化版，只用影像本身加太阳高度角，把场景中的暗像元（假设接近0反射率）的离差当作大气路径辐射来扣除。它的优点是不依赖外部大气数据，批量稳定性好；缺点是暗目标选不准时，水体或阴影会被抬成灰白色。
- QUAC（ENVI内置）：基于场景内光谱多样性的统计估算方法，不需要输入气溶胶和水汽参数，跑得稳，界面上点一下就出结果。它的代价是物理解释性弱，遇到大面积城市高亮区时统计偏差会比6S明显。
对Landsat8来说，高精度路线其实是直接下载Collection 2 Level-2地表反射率产品（L2），官方LaSRC算法生成，附QA_PIXEL云掩膜波段。但实际工程里拿到的数据经常只有L1TP，或者方案要求自己掌控每步处理，本地大气校正仍然绕不开。我一般在大批量场景优先用QUAC，因为它不会因为缺输入而中断；在需要透明公式说明时用COST；只有单景科研分析才跑6S。
| 选型 | 输入依赖 | 批量稳定性 | 物理透明度 | 适合场景 |
|---|---|---|---|---|
| 6S | 气溶胶、水汽、几何 | 差，缺输入即失败 | 高 | 科研级单景、关键日期 |
| COST | 影像＋太阳高度角 | 稳定 | 中 | 独立批量、时序定性分析 |
| QUAC | 仅影像本身 | 稳定 | 中低 | 土地覆盖、快速出图 |
| L2产品 | 官方已算好 | 高 | 中高 | 时序SR分析，按云量筛选后使用 |
4.2 COST暗目标参数的设置细节
COST模型里最关键的一个参数是暗像元DN阈值。取多大的百分比，直接影响大气校正后反射率整体抬升还是压低。经验做法是取全影像1%到5%累积直方图对应的DN作为暗目标值。Landsat8 OLI信噪比高，干净场景一般取1%到2%累计点就够；如果影像里有大面积深水体、干净云影，暗目标位会偏低，此时改用5%更稳。
批量处理时暗目标阈值必须自动计算，不能逐景手填。常见做法是：对蓝波段、绿波段、红波段分别统计直方图，找累计频率2%的DN；若该DN高于某个绝对门限（比如0~65535里的15），就取该值，否则退回去取直方图里最低且出现频率最高的峰值。整套逻辑可以写进一个函数，随批次的MTL循环自动执行，这样就不会出现同一批30景里某几景因为暗目标选得不当而整体偏色的情况。
COST公式里还有一个容易写反的量：太阳高度角的余角是太阳天顶角，校正公式里用的是cos(天顶角)，实际就是sin(高度角)。如果代码里直接拿SUN_ELEVATION取cos，相当于把0.52左右的分母变成0.87，输出会被系统性拉低0.4倍。这个细节新手经常踩，做批量前最好先拿单景验证一下暗像元的归零效果。
4.3 几何校正与重采样：多景对齐的常见做法
L1TP本身自带精校正与正射结果，单景使用一般不需要再做几何校正。真正需要几何处理的场景有三个：多景镶嵌、与另一源影像配准（比如高分三号预处理后的SAR结果要和Landsat8叠合做联合应用），以及跨UTM分带做统一投影。其中跨分带是最典型的批量坑，因为Landsat8不同轨道会落在不同UTM带里，直接把两景叠起来，同一块地物会错开数十米。
常见做法是选定工作区的中央分带作为目标投影，用gdalwarp把所有影像统一到同一个EPSG再进入后续处理。
gdalwarp -t_srs EPSG:32650 -r bilinear -overwrite \ -tr 30 30 -ot Float32 \ ./L8_UTM49/B2_toa.tif ./warped_B2_UTM50.tif这段命令把原本在UTM 49N的影像重投影到UTM 50N，输出像元30米。gdalwarp会自动读取输入影像原有投影信息做转换，不需要手工干预。-r bilinear表示重采样用双线性；-ot Float32保证反射率数值精度不丢失。批量操作时可以把全部输入拼接成一张列表，逐条执行同样的命令，日志里记录每一景的输入源和目标文件名。
重采样方式选取上，双线性适合绝大多数光谱波段，既能保留地物边界锐度，又不会引入过冲；三次卷积在地表均匀的区域呈现更高清晰度，但在地块边缘会产生振铃效应，植被林缘处尤其明显。分类任务我一般优先双线性，宁可让边界平滑一点，也不要引入伪纹理。全色与多光谱融合时才用最邻近法，保证像元值不被插值污染。
5. Landsat8批量预处理的5个坑：现象、原因、解法
5.1 输出全黑或反射率接近0
现象：批量跑完辐射定标后，某几景输出整幅都是0，统计信息里min=0、max=0，打开影像一片黑。
原因：最常见的是MTL解析没抓到SUN_ELEVATION，sin(0)做了分母，整个矩阵变成无穷或无效值，写盘时被当成0。另一种是代码里搞错键名，例如用BAND_20去查REFLECTANCE_MULT_BAND_2，返回None之后被float()转换抛异常，异常没接住就跳过整景。
解法：parse_mtl函数返回后立刻断言几个关键键存在，缺一个就打印警告并跳过；对SUN_ELEVATION做范围检查，低于5度直接判定处理失败。批量脚本里的try/except不能只吞异常，要把mtl路径和异常信息逐条写进fail清单，跑完后逐个排查。
5.2 批处理中断：投影坐标系不一致
现象：第10景处理得正常，第11景开始报投影不匹配；把所有输出叠到遥感软件里，地面控制点在影像之间能对上，但像元网格错位，量算距离多出几十米。
原因：Landsat8相邻条带可能分属不同UTM分带，L1TP本身带投影信息但分带号不一样。直接把这些影像塞进同一个镶嵌网格，就是不重投影硬对齐，结果必然错位。
解法：批量前先扫描全部MTL或对应TIF的投影信息，统计包含几种UTM分带。超过一种就先把每景重投影到统一工作投影再做辐射定标，或者把几何校正拆成独立步骤放在整个流水线最前面。这个工作顺序上的调整，能省掉后面所有波段的对齐麻烦。
5.3 QUAC把水体变成灰白色
现象：跑QUAC大气校正后，某辖区影像里原本深蓝色的水体变成灰白色，反射率统计值在0.65附近，植被边缘还多出一圈光晕。
原因：QUAC靠场景内像元统计推断大气参数。如果这一景里城市高亮屋顶、裸沙地占比很大，统计过程会高估大气路径辐射，扣除过量后反而把暗像元抬高了。Landsat8 30米分辨率下混合像元多，这个问题比高分影像更容易出现。
解法：给QUAC输入一个经过掩膜的场景，把水体、高反射亮区和云先排除掉再统计；或者换COST模型，把暗目标锁定在干净水体上。如果必须用QUAC，处理完后抽几条水体光谱曲线看趋势，不要只信整景均值。
5.4 CLOUD_COVER字段和实际云量不符
现象：MTL里CLOUD_COVER是2.3，人眼检查却发现有十几块明显云斑，最后交付时被验收方质疑预处理质量。
原因：Landsat8的CLOUD_COVER是场景级粗算值，来自分段云检测算法，山区积雪、亮沙地容易被当成云，薄云又容易漏判，和实际云量不一致并不罕见。
解法：不要依赖这个字段做像素级筛选。批量流程里增加QA_PIXEL波段解析，按Collection 2的位定义把cloud和cirrus像元标出来，生成每景的云掩膜。关键日期再辅以人工抽检快视图，保证进入时序分析的像元确实干净。
5.5 时序统计值偏大：没有处理无效值
现象：同一区域不同月份的反射率统计值整体偏高0.03到0.06，NDVI时间曲线在夏季反而凹陷，人工检查像元后发现云影和薄云被当成了地物参与统计。
原因：Landsat8影像里除了背景外，云、阴影、卷云像元在DN值上没有统一标记，很多像素的DN落在20000以上，转成浮点后仍然在0.8到1.0之间，统计均值被这类高反射噪声拉高。
解法：把QA_PIXEL读取为掩膜，凡是标记为云、阴影、卷云的像元一律在统计前剔除；输出GeoTIFF时给这些像元写NoData。这就是影像侧的数据清洗：不是清理表格重复行，而是清理那些没有物理意义的测量值。
6. 批量结果验证与落地习惯：怎么确定预处理好可用
6.1 快视图检验法
批量产出后先看缩略图再谈精度。把每景的B5、B4、B3合成假彩色快视图，植被呈红色，水体呈黑色，一眼就能看出大气校正过没过头。
for scene in $(ls ./output/*_B4_toa.tif | sed 's/_B4_toa.tif//'); do gdalbuildvrt -separate ./tmp/${scene}_rgb.vrt \ ${scene}_B5_toa.tif ${scene}_B4_toa.tif ${scene}_B3_toa.tif gdal_translate -scale 0 0.6 0 255 -ot Byte -of JPEG \ -outsize 10% 10% ./tmp/${scene}_rgb.vrt ./tmp/${scene}_preview.jpg done这段循环对每个场景先用gdalbuildvrt合成三波段VRT，再用gdal_translate拉伸到0到0.6之间输出JPEG。0.6对应地表反射率的上限，如果某个场景快视图整体发黑，说明大气校正扣除过度或辐射定标系数没生效；整体过曝，则说明系数偏置方向反了。30秒能扫完几十景，比逐波段拉直方图直观得多。
6.2 光谱曲线断点测试
选一块稳定的裸地或深水像元，把同一点在多景影像里的TOA或地表反射率拉出来画曲线。正常情况下曲线应当平滑，同季节不会突然跳变。如果某一天反射率突然抬升0.15以上，先查那景有没有薄云，再看SUN_ELEVATION是不是偏低，最后才怀疑大气校正参数。这个测试不用写复杂代码，导入结果后手工取点即可，但它能快速暴露单景异常。
6.3 存档命名规范
最后把每景输出压缩成一个包，目录名规整为“产品ID_处理日期_校正类型”，包内固定三个子目录：radiometric、surface、mask。radiometric放辐射定标后的TOA，surface放大气校正后的地表反射率，mask放QA掩膜产物。每景附一个处理参数JSON，记录MTL解析出的关键系数、COST暗目标阈值、QUAC版本号。
批量处理这类工作，往往跑一遍不难，难的是半年后再来一批新数据，要能按同一套参数复现。我在项目里被“当时怎么调的来着”坑过三次以后，再也不敢不存参数log和fail清单。现在每次跑完都先把归档目录整理好，再谈结果交付。希望帮到你。
本文还有配套的精品资源，点击获取
