企业官网建设流程全解析

简介：这是一份面向遥感图像处理初学者与测绘、地信方向学生和科研人员的实操型教程资源，围绕 Python 结合 GDAL 完成遥感图像绝对辐射定标展开，帮助读者把原始灰度值转换为具有物理意义的辐射值，为后续定量遥感分析打基础。压缩包共 9 个文件，约 16KB，体量轻便，包含 1 个 Python 主脚本承载定标流程、1 个 JSON 辐射校正参数表存放增益与偏置等系数，另有若干 XML 与 iml 等项目配置文件，便于在 IDE 中直接导入工程并对照阅读。内容覆盖元数据读取、波段遍历、定标参数加载、线性校正计算及结果写出等环节，代码结构清晰，可按实际传感器参数调整系数与输出设置。目前已有 903 人学习下载，适合用来快速理解绝对辐射定标的实现思路并作为自有项目的脚本起点。

1. 一景影像两个季节对不上，问题往往出在 DN 值

一个做农业遥感的朋友，把两景相隔三个月的影像丢进脚本直接算 NDVI，同一块地阈值一刀切，分类结果差了三成，他第一反应是云没去干净。抽了几个纯净像元一看，云确实没有，问题是这两景的 DN 值根本不在同一个物理量纲上。DN 只是传感器输出的原始计数值，随增益设置、积分时间、日地距离和太阳高度角一起漂移，夏天和秋天的同一个地物，DN 可以差出几倍。绝对辐射定标干的就是把 DN 拉回到有物理意义的量纲上：要么变成辐亮度（W·m⁻²·sr⁻¹·μm⁻¹），要么变成表观反射率（无量纲，0~1）。这份基于 python + gdal 的定标实现，做的就是这件事，输出带地理参考的 float32 GeoTIFF，多时相、跨传感器都能直接对比。适合需要自己做预处理流水线、又不想在离线环境里装一整套商业软件的人。

2. 项目骨架与技术选型：GDAL 管 IO、NumPy 管算、配置管传感器

2.1 目录结构与模块职责划分

拿到一份定标代码，我最先看的不是算法，是目录。算法写得再漂亮，如果读取、计算、写盘三件事揉在一个main.py里，想加第二个传感器就得重写。这个项目的分层很清楚：

profile.py和configs/是这套结构的关键。传感器之间的差异，本质上只有三样：波段编号怎么对应、每个波段的乘加系数是多少、量化系数（DN 到物理量之前要不要先除 10000）。把这三样抽出去，calibrate.py里就只剩下纯数学，一份代码能同时跑好几类 L1 级产品。我见过不少同学把系数直接硬编码在循环里，if sensor == 'xxx'一路往下堆，到第五个传感器的时候函数已经两百行，改一个系数要在文件里搜三遍。

reader.py单独拎出来是因为元数据格式不统一。带 MTL 文本的那类产品，字段是KEY = VALUE的等号分隔；欧空局那类 L1C 产品，元数据在 XML 里，还得处理命名空间。解析逻辑放在一处，上层拿到的永远是统一的 dict。

2.2 为什么是 GDAL + NumPy 这套组合

选 GDAL 不是因为情怀。定标输出的影像动辄七八千像素见方、十几个波段，还要保留投影、仿射变换、NoData 和内部瓦片，这几个需求摆在一起，能一次满足的库不多。PIL 存 16 位多波段就已经别扭，投影信息得自己想办法塞进 TIFF tag；rasterio 的接口更 Pythonic，但底层还是 GDAL，某些版本在 window 读写叠加 creation options 的时候行为不一致，我在离线环境里踩过一次，输出影像的 block size 被悄悄改回默认值，压缩率掉了一大截。

NumPy 负责的是逐像元运算。定标本身是纯线性变换，dn * mult + add，用循环写 Python 能跑到天荒地老：一景 8000×8000 的影像有六千四百万像元，纯 Python 循环按每次 100 纳秒算也要六秒多，实际比这慢一个数量级。NumPy 的向量化把这个操作压到内存带宽的水平，同一个波段通常零点几秒就完事。真正花时间的是磁盘 IO，所以后面写盘那一节的参数比算法更值得琢磨。

提示：如果目标机器上没有 GDAL 的 Python 绑定，装gdal包时务必让包的版本和系统 libgdal 的版本对齐，否则 import 阶段就会报符号缺失，这种错误信息很难看懂。

2.3 用 SensorProfile 抽象传感器差异

先把参数结构定下来，后面的代码才有着力点。我一般用一个 dataclass 描述单波段参数，再用一个描述整景：

from dataclasses import dataclass, field from typing import Dict, Optional @dataclass class BandCal: """单个波段的定标参数""" rad_mult: float # 辐亮度增益，DN -> 辐亮度 rad_add: float # 辐亮度偏置 ref_mult: Optional[float] = None # 反射率增益，部分产品直接给 ref_add: Optional[float] = None # 反射率偏置 esun: Optional[float] = None # 大气顶太阳辐照度 W/(m2*um) wavelength: Optional[float] = None # 中心波长，做大气校正时要用 @dataclass class SensorProfile: """一景影像的整体定标配置""" name: str band_map: Dict[int, int] # 源文件波段序号 -> 传感器波段号 bands: Dict[int, BandCal] # 传感器波段号 -> 参数 quantize: float = 1.0 # DN 到量化值的除数，L1C 常见 10000 sun_elev_key: str = "SUN_ELEVATION" # 元数据里太阳高度角的键名 date_key: str = "DATE_ACQUIRED" # 成像日期，算日地距离要用

band_map这个字段看起来多余，实际是血泪经验。有些产品文件里只存了可见光近红外几个波段，波段 1 对应的是蓝，不是海岸带；有些产品把全色波段排在最后。如果代码里写死band_idx + 1当波段号，换一景数据出来的反射率就整体错位，而且数值看起来完全正常——这种错最难查，因为你不会怀疑一个 0.03 到 0.4 之间的反射率是错的。

quantize的存在是因为不同产品的 DN 定义不同。有一类 L1C 产品，DN 除以 10000 才是大气顶反射率，系数不用再乘；而带 MTL 的那类产品，DN 直接乘RADIANCE_MULT就是辐亮度，不需要先除。把这件事做成字段，calibrate.py里就只需要一行q = dn / profile.quantize，不用写分支。

esun只在那条「先算辐亮度再转反射率」的路径上用得到。大气顶太阳辐照度是个跟波段中心波长绑定的常数，各波段的典型值差异很大，短波红外波段的量级只有可见光波段的几十分之一。这个表我一般直接内置进 JSON，不动态查询，因为查询需要联网光谱库，离线环境跑不通。

2.4 MTL 元数据解析：别用正则，用 split

元数据解析这段代码很短，但写法差异会直接影响可维护性。我见过用正则一把梭的版本，遇到字段名里带下划线、值里带空格的行就开始出问题。

def parse_mtl(path: str) -> dict: """解析 KEY = VALUE 格式的元数据文本""" meta = {} with open(path, "r", encoding="utf-8", errors="ignore") as f: for raw in f: line = raw.strip() # 跳过空行、分组标记行和注释 if not line or line.startswith("#") or line.startswith("GROUP"): continue if "=" not in line: continue key, value = line.split("=", 1) # 只切第一个等号，值里可能还有等号 meta[key.strip()] = value.strip().strip('"') return meta

逻辑上就三件事：按行读、按第一个等号切、去掉两侧引号。split("=", 1)这个1很重要，某些字段的值是带等号的字符串，不限制切分次数会把值切碎。strip('"')处理的是部分产品给字符串字段加引号的情况，不加这一步，后面按日期解析就会失败。

拿到 dict 之后，按前缀抽系数：

def build_profile_from_mtl(meta: dict) -> SensorProfile: bands = {} for key, value in meta.items(): if not key.startswith("RADIANCE_MULT_BAND_"): continue b = int(key.rsplit("_", 1)[1]) # 从键名末尾拿波段号 bands[b] = BandCal( rad_mult=float(value), rad_add=float(meta[f"RADIANCE_ADD_BAND_{b}"]), ref_mult=float(meta.get(f"REFLECTANCE_MULT_BAND_{b}", "nan")), ref_add=float(meta.get(f"REFLECTANCE_ADD_BAND_{b}", "nan")), ) band_map = {i + 1: b for i, b in enumerate(sorted(bands))} return SensorProfile(name="mtl-l1", band_map=band_map, bands=bands)

rsplit("_", 1)[1]是取键名最后一段，比正则省事也更快。band_map用enumerate按排序后的波段号生成，前提是源文件里的波段顺序和波段号升序一致——这个前提在绝大多数产品里成立，但我在cli.py里留了一个--band-map参数，遇到不成立的情况可以手动覆盖。

数值转换用float()，不包 try。元数据里的系数偶尔会出现nan字符串，float("nan")是合法的，会得到一个 NaN，后面判断ref_mult is not None时 NaN 会走进去，算出来的整幅影像全是 NaN。所以更稳妥的写法是判断math.isnan，这一点在排查章节还会提到。

3. 核心定标算法拆解：辐亮度、反射率与太阳几何校正

3.1 辐亮度定标：一次线性变换，但顺序有讲究

辐亮度是最基础的一步，公式本身没有任何玄学：

L_λ = M_L × Q_cal + A_L

M_L是绝对辐射定标增益（Radiance Mult），A_L是偏置，Q_cal是量化后的 DN 值。落到 NumPy 上就一行，但这一行之前必须处理数据类型：

import numpy as np def dn_to_radiance(dn: np.ndarray, cal) -> np.ndarray: """DN -> 辐亮度，单位 W/(m2*sr*um)""" # 必须先转 float32：uint16 直接参与浮点运算会先升到 float64， # 内存翻倍；若走整数路径还会在减法时回绕 dn = dn.astype(np.float32, copy=False) return dn * np.float32(cal.rad_mult) + np.float32(cal.rad_add)

astype(np.float32, copy=False)里的copy=False是让 NumPy 在类型已经匹配时直接复用原数组，避免多一次内存拷贝。对 8000×8000 的 uint16 数组来说，一次全量拷贝就是 128 MB，多波段叠加起来很可观。把系数显式转成np.float32是为了防止 NumPy 的类型提升规则把整个结果推成 float64——Python 的float是双精度，float32 数组 * Python float会得到 float64 数组。

这一步的输出单位是 W·m⁻²·sr⁻¹·μm⁻¹，量级通常在 0 到几百之间。如果算出来的最大值只有 0.0x，八成是rad_mult读错了，很多产品的增益是2.0000E-05这种科学计数法，解析时按字符串截断就会变成 2。

3.2 表观反射率：把太阳高度角补回去

反射率有两条路径，走哪条取决于元数据里给了什么系数。

第一条是官方直接给反射率系数的路径，先算未校正的表观反射率：

ρ' = M_ρ × Q_cal + A_ρ

第二条是只给了辐亮度系数的路径，用辐亮度反推：

ρ = π × L_λ × d² / (ESUN_λ × cos θ_SZ)

这里的d是日地距离（天文单位），θ_SZ是太阳天顶角。两条路径最后都要做太阳高度角归一化，ρ = ρ' / sin(θ_SE)。最容易翻车的地方在于：走第二条路径时，cos θ_SZ已经等于sin θ_SE，再除一次就等于除了两遍，反射率整体偏大。我第一次写这个函数的时候就是这么错的，出来的影像肉眼看着还挺正常，直到拿官方 L2 产品做对比才发现系统性偏高。

def earth_sun_distance(doy: int) -> float: """由一年中的第几天算日地距离因子，精度足够定标使用""" return 1.0 - 0.0167 * np.cos(2.0 * np.pi * (doy - 3.0) / 365.0) def dn_to_reflectance(dn: np.ndarray, cal, sun_elev_deg: float, quantize: float = 1.0, doy: int = 180) -> np.ndarray: """DN -> 表观反射率，结果无量纲，正常范围 0~1""" q = dn.astype(np.float32, copy=False) / np.float32(quantize) sin_elev = np.float32(np.sin(np.deg2rad(sun_elev_deg))) if cal.ref_mult is not None and not np.isnan(cal.ref_mult): # 路径一：官方反射率系数，得到的是未做太阳角校正的 rho' rho_uncorr = q * np.float32(cal.ref_mult) + np.float32(cal.ref_add) return rho_uncorr / sin_elev # 路径二：辐亮度反推，cos(天顶角) 已等价于 sin(高度角)，不再重复除 radiance = dn_to_radiance(q, cal) d = np.float32(earth_sun_distance(doy)) denom = np.float32(cal.esun) * sin_elev return (np.pi * radiance * d * d / denom)

参数说明：sun_elev_deg是太阳高度角，单位度，直接从元数据的SUN_ELEVATION字段读，别用天顶角，两者相加等于 90；quantize默认 1.0，L1C 那类产品传 10000；doy是成像日期的年内序号，1 月 1 日算 1，用datetime.strptime(date_str, "%Y-%m-%d").timetuple().tm_yday拿。

sin_elev有个边界情况：极区冬季的影像太阳高度角可能接近 0 甚至为负，除下来会得到天文数字。我一般加一道保护，sun_elev_deg < 3的影像直接跳过并打日志，这种数据本身也没法做定量分析。

3.3 分块读写与 NoData 掩膜

大影像不能整幅读进来。一景 8000×8000、11 个波段的 uint16 数据是 1.4 GB，读成 float32 就是 2.8 GB，如果代码里再产生几个中间数组，16 GB 内存的机器直接 OOM。分块处理是标配：

def process_band(src_ds, src_idx: int, dst_ds, dst_idx: int, fn, block: int = 512, nodata_out=np.nan): """按块读取 -> 定标 -> 写回，整幅影像内存占用只跟 block 有关""" src_band = src_ds.GetRasterBand(src_idx) src_nodata = src_band.GetNoDataValue() dst_band = dst_ds.GetRasterBand(dst_idx) dst_band.SetNoDataValue(float(nodata_out)) width, height = src_ds.RasterXSize, src_ds.RasterYSize for y in range(0, height, block): rows = min(block, height - y) for x in range(0, width, block): cols = min(block, width - x) arr = src_band.ReadAsArray(x, y, cols, rows) if arr is None: continue # 无效值掩膜：官方 NoData 和常见的 0 填充都要盖掉 invalid = (arr == 0) if src_nodata is not None: invalid |= (arr == src_nodata) out = fn(arr).astype(np.float32, copy=False) out[invalid] = nodata_out # 用 NaN 标记，后续统计会自动跳过 dst_band.WriteArray(out, x, y)

block取 512 是因为大多数 GeoTIFF 的内部瓦片是 256 或 512，块大小跟瓦片对齐时，GDAL 的一次ReadAsArray正好落在一到两个瓦片上，不会触发额外的重排。取 1024 或 2048 也能跑，内存占用线性上升，但 IO 次数下降，实际速度差别不大，因为瓶颈在解压。

invalid |= (arr == 0)这一行值得单独说。绝对定标里 DN=0 几乎总是填充值，但有些产品的有效地物 DN 也可能是 0（深水体、阴影），一刀切会丢掉一点真实数据。稳妥做法是只在元数据明确写了FILL_VALUE = 0的时候才掩，或者在cli.py里加--keep-zero开关。我默认掩掉，因为不掩的后果更严重：0 DN 算出来的反射率是一个负偏置值，混进统计里会把整幅影像的最小值拉到负数。

3.4 输出写盘：压缩参数比算法更影响体验

输出这一段的参数选择，直接决定了文件是 800 MB 还是 200 MB：

import osgeo.gdal as gdal def create_output(path: str, ref_ds, count: int) -> gdal.Dataset: options = [ "TILED=YES", # 内部瓦片，配合后面的分块随机读 "BLOCKXSIZE=512", "BLOCKYSIZE=512", "COMPRESS=DEFLATE", # 无损，定标结果不建议用有损压缩 "PREDICTOR=3", # 浮点专用预测器，压缩率提升明显 "ZLEVEL=6", # 压缩级别，6 是速度与体积的平衡点 "BIGTIFF=IF_SAFER", # 超过 4GB 自动切 BigTIFF "NUM_THREADS=ALL_CPUS", # 压缩阶段多线程 ] driver = gdal.GetDriverByName("GTiff") dst = driver.Create(path, ref_ds.RasterXSize, ref_ds.RasterYSize, count, gdal.GDT_Float32, options=options) dst.SetGeoTransform(ref_ds.GetGeoTransform()) dst.SetProjection(ref_ds.GetProjection()) return dst

PREDICTOR=3是最容易被忽略的一项。浮点数据的相邻像元高位往往相同，加这个预测器之后再 DEFLATE，反射率影像的压缩率通常能到 3:1 以上；不加的话，float32 的随机低位很难压，文件体积会是前者的两三倍。BIGTIFF=IF_SAFER是防止多波段输出超过 4 GB 时 GeoTIFF 直接写失败，这种失败在老版本 GDAL 上不会报错，只会在关闭数据集时静默截断。

注意：dst用完必须显式dst = None或者dst.FlushCache()，否则进程退出前数据可能没落盘。我在批处理脚本里因为漏了这一句，跑完三百景发现最后几景全是 0 字节。

4. 命令行与批处理：把单景脚本变成能跑的流水线

4.1 参数设计：一个入口覆盖单景和批量

脚本能跑和能天天跑是两回事。前者手改路径就行，后者必须有稳定的命令行接口：

import argparse def build_parser(): p = argparse.ArgumentParser(description="L1 产品绝对辐射定标") p.add_argument("-i", "--input", required=True, help="输入影像路径") p.add_argument("--mtl", help="元数据文件路径，不给则从 input 同名推断") p.add_argument("-o", "--outdir", default="./out", help="输出目录") p.add_argument("--mode", choices=["radiance", "reflectance"], default="reflectance", help="定标模式") p.add_argument("--bands", default="all", help="要处理的波段，逗号分隔，如 2,3,4；all 表示全部") p.add_argument("--profile", help="传感器配置 JSON，优先于 MTL 自动解析") p.add_argument("--block", type=int, default=512, help="分块大小") p.add_argument("--overwrite", action="store_true", help="覆盖已有输出") p.add_argument("--workers", type=int, default=1, help="并发进程数") return p

--mode用choices限定取值，比自由字符串好，输错了 argparse 会直接提示。--bands支持只定标部分波段，实际业务里很常用：做植被指数只需要红和近红外两个波段，把十个波段全定标一遍纯属浪费磁盘。--profile的优先级高于 MTL 自动解析，这是为了应对元数据缺字段的情况——有些二次分发的产品会把 MTL 精简掉一部分。

--overwrite默认关闭是有意为之。批处理跑一半中断，重跑时跳过已完成的输出，能省很多时间。要强制重算就显式加这个开关，比每次都去手动删目录可靠。

4.2 单景主流程：从打开到落盘的完整链路

把前面几节的零件串起来：

def run_single(args) -> str: src = gdal.Open(args.input, gdal.GA_ReadOnly) if src is None: raise RuntimeError(f"无法打开输入影像: {args.input}") # 1) 拿定标参数：优先读显式配置，否则从元数据构造 if args.profile: profile = load_profile_json(args.profile) else: mtl_path = args.mtl or guess_mtl_path(args.input) profile = build_profile_from_mtl(parse_mtl(mtl_path)) # 2) 拿太阳几何与成像日期 meta = parse_mtl(args.mtl or guess_mtl_path(args.input)) sun_elev = float(meta.get(profile.sun_elev_key, "nan")) doy = parse_doy(meta.get(profile.date_key)) out_path = os.path.join(args.outdir, os.path.basename(args.input)) if os.path.exists(out_path) and not args.overwrite: return out_path # 断点续跑 # 3) 选出要处理的波段对 pairs = [(src_idx, b) for src_idx, b in sorted(profile.band_map.items()) if args.bands == "all" or b in parse_band_list(args.bands)] dst = create_output(out_path, src, len(pairs)) for out_idx, (src_idx, band_no) in enumerate(pairs, start=1): cal = profile.bands[band_no] if args.mode == "radiance": fn = lambda a, c=cal: dn_to_radiance(a, c) else: fn = lambda a, c=cal: dn_to_reflectance( a, c, sun_elev, profile.quantize, doy) process_band(src, src_idx, dst, out_idx, fn, block=args.block) dst.GetRasterBand(out_idx).SetDescription(f"B{band_no}") dst.FlushCache() dst = None src = None return out_path

lambda a, c=cal: ...里的c=cal是 Python 闭包的经典坑：循环里创建的 lambda 如果不绑定默认参数，等到真正调用时cal已经指向最后一个波段，所有波段的定标结果都会用错的系数。这个 bug 不会报错，只会让输出影像的波段间关系彻底错乱。

SetDescription(f"B{band_no}")把原始波段号写进波段描述。QGIS 里打开影像，图层名会显示成B4而不是Band 1，做后续分析的时候不用再翻元数据对照。

pairs的排序按源文件波段序号来，保证输出波段的顺序和输入一致。如果中间用--bands 4,3,2过滤，输出就是三个波段，顺序按源文件排，不按你输入的逗号顺序——这个行为要写进--help，不然用户会以为能靠参数顺序调整波段排列。

4.3 多进程批处理与日志

单景跑通之后，批量就是加一层调度。定标是 CPU 密集型，多进程比多线程有效：

from concurrent.futures import ProcessPoolExecutor, as_completed import logging def run_batch(file_list, args): logging.basicConfig( level=logging.INFO, format="%(asctime)s | %(levelname)s | %(message)s", handlers=[logging.FileHandler("radcal.log"), logging.StreamHandler()], ) ok, failed = 0, [] with ProcessPoolExecutor(max_workers=args.workers) as pool: futures = {pool.submit(run_single, args, f): f for f in file_list} for fut in as_completed(futures): name = os.path.basename(futures[fut]) try: out = fut.result() ok += 1 logging.info("完成 %s -> %s", name, out) except Exception as exc: failed.append(name) # 单景失败不中断整批，把原因记下来后继续 logging.exception("失败 %s: %s", name, exc) logging.info("批量结束：成功 %d，失败 %d", ok, len(failed)) return failed

as_completed而不是按提交顺序取结果，是为了让快的先落日志，用户能实时看到进度。单景异常用logging.exception记完整堆栈，不往外抛，不然一景数据损坏就会让整批几百景中断，重跑还得从头开始。

workers不建议超过物理核数。定标过程中内存占用是workers × 单景峰值内存，八个进程同时处理 8000×8000 的十一波段影像，32 GB 内存也会吃紧。我一般按min(核数, 内存GB / 4)来估。

5. 常见问题排查：从数值异常到 IO 卡顿

5.1 反射率出现大于 1 的值，是不是算错了

现象：输出影像的最大值到了 1.5 甚至 3，直方图右侧拖了一条长尾。

原因：表观反射率本来就可以超过 1。云顶、雪面、水面镜面反射、以及太阳耀斑区域，DN 值接近量化上限，算出来的反射率超过 1 是正常物理现象。真正的错误是整幅影像的均值在 0.8 以上，那才是系数或太阳高度角出了问题。

解决：先看直方图分布，如果 95% 分位数在 0.05~0.5 之间，长尾只是少数亮目标，不用管。如果整体偏高，检查是不是重复除了sin(θ_SE)——把sun_elev临时设成 90 度再跑一遍，结果不变说明公式里有多余的除法。

5.2 输出影像全黑或者全 NaN

现象：QGIS 打开只有一片黑，或者统计信息显示 min/max 都是 NaN。

原因通常有三个：一是quantize没设对，L1C 产品没除 10000，DN 上千，乘 1e-4 量级的系数是对的，但如果元数据里的系数本身已经是量化后的，再除就小了四个数量级；二是 NoData 掩膜把整幅影像盖掉了，比如src_nodata读出来是 0，而invalid = (arr == 0)又叠加了一次，正常数据不可能是 0，但如果某景的填充值恰好是 65535，而有效值范围也覆盖到了，就会误伤；三是输出时dst没置 None，缓冲区没落盘。

解决：先用gdalinfo -stats看输入影像的 DN 分布，再单独取一个已知地物（比如城区）的像元，手算一遍公式对数值。这一步花五分钟，比盲改代码一小时有用。

5.3 分块处理的耗时比整幅读取还长

现象：设了block=512之后，跑一景要十几分钟，以前整幅读只要两分钟。

原因：ReadAsArray的窗口和 TIFF 内部瓦片不对齐，每次读取都触发一次完整的瓦片解码再裁剪，同一个瓦片被反复解码几十次。

解决：先用gdalinfo看Tile Size，把block设成它的整数倍；如果影像是条带式（BLOCKYSIZE很大），先gdal_translate -co TILED=YES转一次再处理。这条我在处理一批二次分发的产品时踩过，那批数据被重新编码成 strip 格式，分块读的代价比顺序读高一个数量级。

5.4 多波段输出波段顺序错乱

现象：输出影像的第三个波段在 QGIS 里显示出来是近红外，跟输入的波段 3（红）对不上。

原因：band_map构造时用了sorted(bands)，而某些产品的元数据里波段编号不是从 1 连续开始的，比如缺了波段 8，enumerate生成的序号就和源文件里的波段顺序错开了。

解决：band_map不要靠enumerate猜，直接从元数据里读FILE_NAME_BAND_x的顺序，或者用--band-map 1:1,2:2,3:3,4:4手动指定。写完输出后拿gdalinfo -nomd核对一遍波段数和描述，比事后发现好。

5.5 批处理跑到一半内存被打满

现象：并行 8 个进程，跑到第七景的时候系统开始 swap，速度断崖式下跌。

原因：每个进程的峰值内存是block² × 4 字节 × 中间数组个数，加上 GDAL 自己的瓦片缓存，实际占用比估算高一截。GDAL 默认的块缓存是物理内存的 5%，八个进程就吃掉 40%。

解决：在程序启动时调gdal.SetCacheMax(256 * 1024 * 1024)把每个进程的缓存压到 256 MB，workers按内存 / 单景峰值保守取值。我现在的习惯是先在单景上跑一遍看htop的峰值，再定并发数，不靠拍脑袋。

6. 二次开发方向：从能跑到好用之间要补的几块

这套代码的骨架不大，改造空间主要在三处，按投入产出排一下：

新增传感器是最划算的。把build_profile_from_mtl复制一份改成build_profile_from_xml，用xml.etree.ElementTree解析，注意命名空间{http://www.esa.int/safe/sentinel-1.0}这类前缀，用iterparse或者findall(".//{*}TagName")通配命名空间最省事。系数映射关系写进 JSON，剩下的代码一行不动。

输出 COG 值得单独提一句，如果定标结果要上传到对象存储给在线分析用，普通 GeoTIFF 的读取效率很差。转换很简单：

def to_cog(src_path: str, dst_path: str) -> None: """把普通 GeoTIFF 转成 Cloud Optimized GeoTIFF，带金字塔""" gdal.Translate( dst_path, src_path, format="COG", creationOptions=[ "COMPRESS=DEFLATE", "PREDICTOR=3", "OVERVIEWS=IGNORE_EXISTING", # 没有金字塔时自动生成 "BLOCKSIZE=512", ], )

OVERVIEWS=IGNORE_EXISTING会让 GDAL 在转换时按 2 的幂次生成概览层，缩放到小比例尺时只读概览不读原图，在线浏览的响应速度差别很明显。BLOCKSIZE保持 512，和前面的分块参数对齐。

验证改造是否正确的办法，我一般用抽样对比而不是目视。找一景同产品已经发布的大气顶反射率产品，用gdal.Translate抽出同样的几百个像元坐标，两边相减算 RMSE。城市、水体、植被各取一批，RMSE 在 0.01 以内就算通过。目视对比不靠谱，显示拉伸会把 5% 的系统偏差完全抹掉。

最后说个习惯。这个项目我一开始图快，把系数解析、计算、写盘全塞在cli.py里，单景跑得挺好，等到要加第二个传感器的时候，改一处崩三处，最后花了两天重构回现在的分层。从那以后我每次新开一个定标脚本，都强制先写profile.py和calibrate.py两个文件，跑通一个波段的单元测试再碰 IO。定标这件事，公式本身没有难度，难的是数据源之间的差异和边界情况，把差异挡在配置层外面，剩下的路就好走了。希望帮到你。

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