简介：围绕遥感影像绝对辐射定标，这份资料面向具备一定 Python 基础、希望用 GDAL 将原始 DN 值转换为物理辐射量的遥感学习者与图像处理人员。内容以脚本实践为主，核心是一个 Python 处理程序，配合 JSON 定标参数组织各波段的增益、暗电流等信息，可据此理解参数读取、逐波段校正、结果写回及元数据更新的实现路径。压缩包共 9 个文件，以 5 个 XML 工程配置、1 个 Python 脚本和 1 个 JSON 参数文件为主，另含工程辅助文件，整体约 16KB，结构轻量，适合快速查看代码框架并移植到自己的数据流程中。目前已有 903 人学习下载。借助其中的代码骨架与参数示例，读者可掌握 GDAL 波段读写、线性辐射校正计算和输出数据集创建的基本方法，也能据此排查波段对应、参数匹配与文件保存等环节的问题，为后续大气校正和定量遥感分析打基础。
1. DN 直算 NDVI 差多少：定标这道坎绕不过去
把一张多波段遥感影像读进 numpy，减一减、除一除算出 NDVI，图一画看着挺像回事——这是很多人用 python 和 gdal 起步时踩的第一块地雷。DN 值只是传感器记录的整数计数，同一块地不同日期、不同传感器拍出来完全不同，直接拿去算指数、做时序对比，结论在物理上站不住脚。绝对辐射定标就是把这些整数拉回有量纲的辐亮度、或者无量纲的大气顶反射率，把增益、偏移、太阳高度角这些成像条件统一掉。它不改变地物的相对分布，改变的是你说这句话时的底气：这块水体反射率 0.03、那片植被近红外 0.28，是能写进报告、能和别的影像放一起比对的数字。凡是跨时相、跨传感器、要把结果交给别人复现的活，这一步都省不得。
2. 从 DN 到 TOA 反射率：三个物理量与两条计算路径
2.1 DN、辐亮度、反射率各自的量级感
DN 是整型，L1 产品常见 uint16，12bit 有效数据落在 0～4095 或 0～65535，取决于产品封装方式。它没有单位，只表示探测器输出的相对强弱。辐亮度的单位是 W/(m²·sr·μm)，量级通常在 0～700 之间，蓝天、厚云、裸岩能拉开好几倍差距，出问题的时候一眼能看出来。TOA 反射率的量纲是 1，正规产品会落在 0～1 之间（部分产品存成 0～10000 的整数），水体在可见光通常 0.02～0.08，健康植被近红外能到 0.3～0.5，裸土可见光 0.15～0.30。
很多新手会问：既然最后都要用反射率，为什么不跳过辐亮度？因为辐亮度是传感器真正测到的量，所有的定标系数、交叉定标、场地替代定标都是围绕辐亮度做的。做传感器性能分析、做绝对定标精度评估，必须停在辐亮度这一层；做地物识别、时序分析、植被指数，才继续往上走。两条路径的中间产物不一样，用途也不一样，混着用就是事故的起点。
2.2 两组乘加系数：元数据里到底给了什么
主流 L1 产品的元数据一般会给两组线性系数，形式都是 y = mult × DN + add。第一组把 DN 变成辐亮度，第二组把 DN 直接变成反射率。第二组的系数在生成时其实已经把波段平均太阳辐照度和日地距离揉进去了，所以你拿到手就是「差不多可以用的反射率」，只剩太阳高度角需要你自己除一下。
这里有个非常隐蔽的坑：既然有两组系数，就一定有人拿辐亮度系数去套反射率公式。算出来的数量级能差几十倍，图上还是花花绿绿，肉眼一时半会看不出来，等到拿去和标准产品比对才发现全错。判断方法很简单，看系数本身的数量级：反射率系数的 mult 通常是 1e-4 到 1e-5 这个量级，加了之后正好落在 0～1；辐亮度系数会大一些，具体取决于波段。看到 1e-2 这种量级的 mult，多半是辐亮度系数。
2.3 太阳高度角和日地距离：哪条路要除、哪条路不用
从 DN 走反射率系数这条路径，公式是：
ρ = (Mρ × DN + Aρ) / sin(θ_SE)
θ_SE 是元数据里的太阳高度角，单位是度，代码里必须转弧度再取 sin。这个除法是很多人的翻车点——用 cos 的、用天顶角的、忘了转弧度的，都能让结果整体偏移几十个百分点。太阳高度角是「地面与太阳的仰角」，不是天顶角，两者加起来是 90 度，元数据字段名里带 ELEVATION 的才是高度角。
从辐亮度走反射率这条路径，公式多一项日地距离：
ρ = π × L × d² / (ESUN × cos(θ_S))
d 是日地距离，单位天文单位 AU，随一年中的日序变化，在 0.983 到 1.017 之间摆动，d² 带来的差异最大约 7%。ESUN 是波段的平均太阳光谱辐照度，从传感器文档附表里抄，不同来源能差 1%～2%，对 TOA 反射率影响很小，可以接受。这两条路径的差别不在精度，在于「谁已经把校正做完了」：走系数路径就不要额外乘 d²，否则等于重复校正；走辐亮度路径就必须带上 d²，漏了在冬季高纬会明显偏低。
日地距离不需要查表，Spencer 近似就够：
import math def earth_sun_distance(doy): """Spencer 近似，输入年积日(1-366)，返回日地距离(AU)""" g = 2.0 * math.pi * (doy - 1) / 365.0 return (1.00011 - 0.034221 * math.cos(g) - 0.001280 * math.sin(g) - 0.000719 * math.cos(2 * g) - 0.000077 * math.sin(2 * g))参数说明：doy 从元数据里的采集日期算出来，比如datetime.strptime(date_str, "%Y-%m-%d").timetuple().tm_yday。这个近似式在全年范围内的误差小于 0.01%，做 TOA 级别的定标绰绰有余，不用上更复杂的星历。
提示：走系数路径时，元数据给的反射率系数已经内含日地距离归一化，代码里再乘 d² 会把 7% 的季节性偏差重新引进来，表现为「夏天的影像普遍偏亮」。
3. GDAL 落地：元数据解析、分块读写与批量输出
3.1 把元数据解析成字典并校验完整性
文本型元数据是最常见的一种，一行一个KEY = VALUE，直接按等号切就行。这里不建议用正则硬抠字段名，因为不同产品版本会加前后缀，键名统一转小写再查，兼容性最好。
from pathlib import Path def parse_kv_meta(path): """解析 KEY = VALUE 形式的文本元数据，键统一转小写""" meta = {} for raw in Path(path).read_text(encoding="utf-8", errors="ignore").splitlines(): line = raw.strip() if not line or line.startswith("#") or "=" not in line: continue k, _, v = line.partition("=") meta[k.strip().lower()] = v.strip() return meta读完之后别急着算，先校验。定标系数缺一个，后面全白干，而且往往是某几个波段缺，整景跑完才发现就晚了。
REQUIRED_TPL = ( "radiance_mult_band_{b}", "radiance_add_band_{b}", "reflectance_mult_band_{b}", "reflectance_add_band_{b}", "sun_elevation", ) def build_band_calib(meta, band, route="reflectance"): """取出指定波段的系数，缺失直接抛错，避免静默算错""" missing = [t.format(b=band) for t in REQUIRED_TPL if t.format(b=band) not in meta] if missing: raise KeyError(f"波段 {band} 缺少元数据字段: {missing}") if route == "reflectance": return { "mult": float(meta[f"reflectance_mult_band_{band}"]), "add": float(meta[f"reflectance_add_band_{band}"]), "sun_elev": float(meta["sun_elevation"]), "route": "reflectance", } return { "mult": float(meta[f"radiance_mult_band_{band}"]), "add": float(meta[f"radiance_add_band_{band}"]), "sun_elev": float(meta["sun_elevation"]), "route": "radiance", }逻辑说明：build_band_calib把「取系数」和「选路径」绑在一起，调用方只能声明自己要反射率还是辐亮度，不能自己拼字段名。这样以后换产品，只改这一个函数，主流程不动。参数route只有两个合法值，其它值直接走 else 分支返回辐亮度，所以我一般会在函数开头加一句断言把非法值挡掉，避免手滑传错字符串悄悄走错路。
3.2 分块读写：别把整景读进内存
一景中分辨率影像单波段动辄上亿像素，uint16 读进来 200MB，转 float32 再翻倍，同时开三四个波段就顶到内存天花板了。gdal 的ReadAsArray支持窗口读取，按行分块是最省心的做法，块高取 512 或 1024，兼顾 IO 效率和内存占用。
from osgeo import gdal import numpy as np gdal.UseExceptions() def calibrate_reflectance(src_path, dst_path, band, mult, add, sun_elev_deg, chunk=1024, src_nodata=0): src = gdal.Open(src_path, gdal.GA_ReadOnly) w, h = src.RasterXSize, src.RasterYSize dst = gdal.GetDriverByName("GTiff").Create( dst_path, w, h, 1, gdal.GDT_Float32, options=["TILED=YES", "COMPRESS=DEFLATE", "BLOCKXSIZE=256", "BLOCKYSIZE=256", "BIGTIFF=IF_SAFER"]) dst.SetProjection(src.GetProjection()) dst.SetGeoTransform(src.GetGeoTransform()) sin_se = math.sin(math.radians(sun_elev_deg)) if sin_se <= 0: raise ValueError(f"太阳高度角异常: {sun_elev_deg}") src_band = src.GetRasterBand(band) out_band = dst.GetRasterBand(1) out_band.SetNoDataValue(float("nan")) for y in range(0, h, chunk): rows = min(chunk, h - y) dn = src_band.ReadAsArray(0, y, w, rows).astype(np.float32) rho = (dn * mult + add) / sin_se rho[dn == src_nodata] = np.nan # 填充值必须置 NaN out_band.WriteArray(rho, 0, y) del dn, rho out_band.FlushCache() dst = None src = None逻辑说明与参数说明：
gdal.UseExceptions()让失败变成异常而不是返回 None，省掉一堆if ds is None。COMPRESS=DEFLATE对浮点影像的压缩比一般，但能把体积压到三分之一左右；如果你的下游工具对压缩敏感，换成LZW。BIGTIFF=IF_SAFER是给大区域镶嵌留的后路。chunk=1024是个折中值。块太小 IO 次数多，块太大内存吃紧，1024 行在 30 米分辨率、7000 像素宽的景上大约是 30MB 的 float32 临时数组，很安全。src_nodata=0要按产品实际情况改。有的产品填充值是 0，有的是 65535，还有的是 -9999，不确定就先跑一句src_band.GetNoDataValue()看一眼，别想当然。- 输出统一 float32，不要用 float64，体积翻倍但精度对你没有意义；也不要输出 uint16，会丢掉负值和小于 1 的小数。
3.3 批量处理整景目录与输出命名
单景跑通之后就是批量。目录结构每家都不一样，我的做法是先扫目录找出元数据文件，再按元数据里登记的波段文件名去定位影像，而不是靠文件名猜。这样换产品、换命名规则都不用改代码。
from pathlib import Path def batch_calibrate(scene_dir, out_dir, bands=(4, 5), route="reflectance"): scene_dir, out_dir = Path(scene_dir), Path(out_dir) out_dir.mkdir(parents=True, exist_ok=True) mtl = next(scene_dir.glob("*_MTL.txt")) meta = parse_kv_meta(mtl) for b in bands: # 优先用元数据里登记的波段文件名，找不到再退回通配 key = f"file_name_band_{b}" if key in meta: src = scene_dir / meta[key] else: src = next(scene_dir.glob(f"*_B{b}.TIF")) cal = build_band_calib(meta, b, route=route) if cal["route"] == "reflectance": calibrate_reflectance(src, out_dir / f"B{b}_toa.tif", 1, cal["mult"], cal["add"], cal["sun_elev"]) else: calibrate_radiance(src, out_dir / f"B{b}_rad.tif", 1, cal["mult"], cal["add"])逻辑说明：bands=(4, 5)默认只出红和近红外，够算 NDVI 了，全波段跑一遍没必要。out_dir单独放，不要写回原目录，否则下次扫目录会把定标结果当成原始影像再处理一遍。命名里带上toa或rad后缀，半年后回来看一眼就知道这是什么层级的产物。
4. 多产品差异：MTL 文本与 XML 元数据的系数处理
4.1 文本型产品：系数已含太阳几何归一化
文本型元数据的产品（一景一个*_MTL.txt，波段系数全在里面）比较好处理，因为反射率系数是逐波段独立给的，你只需要除一个 sin(太阳高度角)。但要注意两点。
第一，太阳高度角是整景一个值，不是逐像元。地形起伏大的区域，这个近似会带来明显误差，山区做定量分析得上 DEM 做地形校正，那是另一个话题。第二，部分产品的反射率系数是「按波段中心波长算的平均值」，全色波段和热红外波段不提供反射率系数，只提供辐亮度系数。所以你写批量脚本时，波段列表里混进全色波段会直接抛 KeyError，得提前过滤。
热红外波段的处理更特殊，它的物理量是亮温不是反射率，用辐亮度反演亮温需要中心波长和两个常数 K1、K2，这套流程和可见光近红外完全不是一回事，别塞进同一个循环里。
4.2 XML 型产品：量化值和负偏移
XML 型元数据的产品结构不一样，整景一个*_MTD_MSIL1C.xml，每个波段一个<Band>节点，属性band_id从 0 开始编号，而波段号从 1 开始。这个错位一不留神就把 B4 的系数套到 B3 上，算出来的 NDVI 全是错的，而且因为数值范围看着正常，很难第一眼发现。
另一个差异是量化。DN 先被存成整数，元数据给一个QUANTIFICATION_VALUE（常见是 10000），真实反射率是 DN 除以它。近几年的处理基线还给了一个负的偏移量，较早的标签叫BOA_ADD_OFFSET，新基线的 L1C 里叫RADIO_ADD_OFFSET，含义都是在除法之前先加回来。公式变成：
ρ = (DN + offset) / quantification
import xml.etree.ElementTree as ET def parse_xml_meta(xml_path): root = ET.parse(xml_path).getroot() quant = float(root.findtext(".//QUANTIFICATION_VALUE")) offsets = {} # 兼容新旧两种标签，有哪个用哪个 for tag in ("RADIO_ADD_OFFSET", "BOA_ADD_OFFSET"): for node in root.iter(tag): offsets[int(node.attrib["band_id"])] = float(node.text) if offsets: break return quant, offsets逻辑说明：findtext(".//QUANTIFICATION_VALUE")用.//做全树搜索，因为不同产品把这个节点放在不同层级。offset 部分做了两轮尝试，先找新标签，找不到再退回旧标签，两轮都空就返回空字典，调用方按 offset=0 处理——老产品确实没有这个字段，不能因为找不到就报错。参数band_id是 0 基的，取系数时务必用band_id = band_no - 1。
注意：新基线产品在未减偏移的情况下，水体等暗目标的反射率会被整体抬高约 0.1，做水体提取会大面积漏提，而且这种错误在图上看不出来，只能靠统计和比对发现。
4.3 抽成统一参数表：一份配置跑两种产品
两种产品的差异集中在「系数从哪来、要不要减偏移、要不要除 sin」三件事上。把这三点收进一个结构体，主流程就只剩「读 → 算 → 写」，后面加第三种产品也只改解析层。
from dataclasses import dataclass @dataclass class BandCalib: mult: float = 1.0 add: float = 0.0 sun_elev_deg: float = 90.0 physical: str = "reflectance" # reflectance / radiance def apply(self, dn): import numpy as np out = dn.astype(np.float32) * self.mult + self.add if self.physical == "reflectance": out /= np.sin(np.radians(self.sun_elev_deg)) return out这样写的好处是，apply里只有一个公式，两种产品的差异被压到构造BandCalib的那几行里。文本型产品的 mult 直接用元数据给的，XML 型产品的 mult 是1/quantification、add 是offset/quantification，再配合太阳高度角，出来的就是同一个物理量。参数physical用来区分反射率和辐亮度两种输出，辐亮度路径直接跳过除法。
5. 避坑与排查：五条定标结果异常的现场记录
5.1 反射率整体大几十倍
现象：算出来的反射率最大值到了 30 多，直方图整体右移，但图还是能出。
原因：拿辐亮度系数套了反射率公式。两组系数的数量级差两三个量级，代码里字段名只差radiance和reflectance几个字母，复制粘贴最容易错。
解决：在build_band_calib里加一道数量级断言，反射率系数 mult 大于 1e-2 就直接抛错。同时把取字段的逻辑收进函数，不让调用方自己拼字段名。
5.2 影像边缘出现一条 0 值黑边
现象：定标后的影像左上角或四周有一圈 0，算 NDVI 时变成无效值，统计出来地物比例全乱。
原因：原始 DN 的填充值是 0，套上公式后变成add的值（一个很小的正数），没有置成 NoData。下游工具把它当有效像元参与统计。
解决：在分块循环里显式做rho[dn == src_nodata] = np.nan，并且给输出波段设SetNoDataValue(float("nan"))。注意比较要在转 float32 之前做，或者用dn == 0的整型比较，浮点比较容易因为精度漏掉。
5.3 跨年份比对出现系统性偏移
现象：同一块农田，2020 年的 NDVI 均值 0.72，2023 年的 0.61，物理上不可能差这么多。
原因：新基线的产品引入了负偏移，处理脚本还在按老公式算，暗目标被整体抬高，NDVI 被压低。元数据里的值是有的，代码里没读。
解决：把 offset 解析做成必选项，读不到就打印警告并记录到日志，同时在输出文件名里带上处理基线的版本号。跨年份分析前，先挑一景两个年份都有的稳定地物（大型水库水面）做对比，两个年份的反射率差异超过 0.02 就要回头查。
5.4 跑到一半内存爆掉
现象：单景能跑，批量跑到第三景被系统杀掉，日志停在ReadAsArray。
原因：整景读入ReadAsArray()不带参数会一次性读全图，一景 8000×8000 的 uint16 转 float32 就是 500MB，加上中间数组和 numpy 的临时副本，峰值能到 2GB 以上。
解决：按行分块读，块高 1024，循环里及时del掉中间数组。另外gdal.Open拿到的数据集和 band 对象要在函数末尾显式置 None，Python 的引用计数在 gdal 对象上不是立刻回收的，循环里创建几百次会累积。
5.5 反射率出现大于 1 和负数
现象：直方图两端有长尾，最大值 1.3，最小值 -0.05。
原因：分三种情况。太阳高度角用错（度当弧度用、用成天顶角），会让整景系统性偏移；云和雪的高亮目标确实能超过 1，属于物理正常；阴影区加上大气散射，出现小的负值也正常。
解决：先看太阳高度角的量级，度制下常见 20～70，弧度制会是 0.3～1.2 的小数，一眼就能分辨。超过 1 的像元统计一下占比，云雪区 5% 以内算正常，全图 20% 以上就要怀疑系数。负值在 -0.05 以内不必处理，超过就检查是不是 offset 减多了。
6. 结果自检：三条量程红线与一次跨传感器比对
定标跑完不算结束，得有一套能当场判死的自检。我习惯跑三条量程红线，全都是几行 numpy 的事，但能挡住八成以上的错误。
第一条，全图分位数。健康植被的近红外 TOA 反射率 95 分位通常在 0.35～0.55，可见光红波段中位数在 0.05～0.15。如果你的红波段中位数跑到 0.4，基本可以确定是把辐亮度当反射率了。第二条，超量程占比。反射率大于 1.2 或小于 -0.08 的像元，在无云场景下应当低于 1%，超过就回头查 offset 和太阳高度角。第三条，波段间关系。植被像元的近红外减红波段的差值分布应当明显右偏，如果 NDVI 的中位数落在 0 附近，说明两个波段的系数张冠李戴了。
import numpy as np def sanity_check(rho_red, rho_nir): valid = np.isfinite(rho_red) & np.isfinite(rho_nir) r, n = rho_red[valid], rho_nir[valid] ndvi = (n - r) / (n + r + 1e-6) report = { "red_p50": float(np.percentile(r, 50)), "nir_p95": float(np.percentile(n, 95)), "over_range": float(np.mean((n > 1.2) | (r > 1.2))), "ndvi_p50": float(np.percentile(ndvi, 50)), } assert 0.02 < report["red_p50"] < 0.20, f"红波段中位数异常: {report}" assert 0.25 < report["nir_p95"] < 0.60, f"近红外 95 分位异常: {report}" assert report["over_range"] < 0.05, f"超量程占比过高: {report}" return report参数说明：分位数的阈值不是死的，拿几景你自己确认过没问题的影像先跑一遍，把实际值记下来当基线，以后每换一批数据就跑一次对比。断言里的上下限要留出季节余量，冬季高纬的植被近红外会低到 0.25 附近，硬套夏季阈值会误报。
更硬的一招是交叉比对：找同一区域、同一天、另一颗卫星的 L1 产品，各自定标后取重叠区，在稳定地物（大型水库、裸岩、机场跑道）上比反射率。两颗卫星的波段响应函数不同，允许 0.02～0.04 的差异，超过 0.1 就一定是某一侧算错了。这个比对不需要大气校正，因为两边都是 TOA 层，条件对等。我一般会挑一块大面积深水，水体在近红外几乎不反射，任何一侧的系数错了都会在这个目标上暴露得最明显。
血泪教训是有一年做长时序分析，脚本跑了大半年没人管，后来发现某一批产品的元数据字段名多了一个后缀，解析函数静默返回了默认系数，整批影像的反射率全都偏了 15%，而因为偏得「均匀」，图上看不出任何异常，直到和另一颗卫星比对才暴露。从那以后我每次加新的数据源，第一件事就是拿一景已知结果的数据跑通sanity_check，把三项指标打印出来存进日志，再开始批量。写脚本的人和用脚本的人常常不是同一个，日志里留下这几个数，半年后回来查也有据可依。希望帮到你。
本文还有配套的精品资源，点击获取
