简介：这是一份“遥感+Python绝对辐射定标”配套工程资源，面向遥感、GIS及Python学习者，解决如何利用GDAL库将遥感图像原始数字值转换为有物理意义的辐射值这一核心问题。压缩包共9个文件、约16KB，包含一个Python定标脚本、一份JSON辐射定标参数文件、若干XML工程配置及Git忽略规则等，结构精简，便于对照教程逐行理解定标流程。内容以绝对辐射定标为主线，覆盖数据准备、元数据与定标参数读取、各波段Band对象遍历、辐射校正值计算，以及按需结合大气校正模型的处理思路；代码框架清晰，可快速改写成多波段批量处理工具，或作为独立模块嵌入其他遥感流程。已有903人学习下载，适合作为课程设计、科研预处理或GDAL二次开发的入门实操样板，尤其有助于理解从DN值到反射率/辐射值的换算链条。

1. 绝对辐射定标：为什么DN转成物理量才是定量分析的前提

做遥感定量分析的人，多半被同一件事坑过：两期影像在同一块农田上算NDVI，趋势一会儿升一会儿降，换了传感器之后所有结论都反了。原因在于DN不是一个物理量，它只是传感器量化输出的一串整数，增益、偏置、太阳高度角不同都会让这串数跑偏。遥感+Python绝对辐射定标要做的，就是把DN转换成辐亮度或表观反射率，让不同时期、不同传感器、不同光照条件下的影像能在同一物理尺度上互相比较。这个操作看起来只是一次乘加和一次余弦改正，但系数从哪里读、单位怎么统一、输出怎么保存，每个环节都有细节。这篇笔记会沿着一条完整链路走：读TIF、解析元数据、批量定标写回，再把验证方法和踩坑记录一次讲清，新老手都能直接照着改。

2. 定标链路：增益、偏置与太阳参数如何把DN变成表观反射率

2.1 绝对辐射定标与相对辐射定标的分工

遥感里说的“定标”其实有两层。相对辐射定标只是把不同影像的DN分布对齐，比如做直方图匹配，让两期影像看起来亮度一致，但输出仍然是DN，不是物理量。绝对辐射定标则必须把DN换算成传感器入瞳处的辐亮度，或者进一步换算成表观反射率。这个区别在做时序分析和定量反演时非常关键：NDVI差值分析用相对定标勉强能过，但要反演叶绿素浓度、水体悬浮物含量，输入必须是辐亮度或反射率这类有单位的物理量，否则经验模型里的系数全部错位。

严格说起来，绝对辐射定标只是整个辐射校正链路的前半段。DN经乘加运算得到辐亮度，再结合太阳角度和日地距离得到表观反射率，这一层一般不涉及大气。真正的地表反射率还需要大气校正，把大气散射和吸收的影响去掉。我在处理业务数据时，会把DN到表观反射率这一段统称为绝对辐射定标，把大气校正作为单独一步处理，这样职责清楚，也方便排查问题。下面所有代码默认输出到表观反射率这一层。

2.2 直读式定标系数：MTL与XML里的Mult和Add

现在主流光学遥感卫星的L1级产品，官方元数据文件里会直接给出定标系数，用户不需要自己去查传感器的绝对定标参数表。以常见的30米多光谱数据为例，它的元数据文件通常叫MTL.txt，里面有一组字段长这样：

RADIANCE_MULT_BAND_1 = 0.0211 RADIANCE_ADD_BAND_1 = -5.2311 REFLECTANCE_MULT_BAND_1 = 0.000023 REFLECTANCE_ADD_BAND_1 = -0.099 SUN_ELEVATION = 56.52

细看会发现，这一份文件里同时给了两套系数：一套是辐亮度定标系数，公式是L = ML * DN + AL；另一套是表观反射率定标系数，公式是ρ = (Mρ * DN + Aρ) / sin(θe)，其中 θe 是太阳高度角。两套系数不能混用，ML与Mρ相差上千倍，代入错的系数，结果立刻乱套。

另有一类10米级多光谱数据的L1C产品，处理方式更简单，量化值DN除以10000就是表观反射率，不需要逐波段查找增益和偏置。但这条规则只适用于该产品的L1C级别，如果拿到的是L2A地表反射率产品，那已经是大气校正后的结果，再做一次“定标”反而画蛇添足。所以动手之前先确认产品级别，这比选Python库还重要。

对于一些科学卫星数据，定标系数不是写在文本里，而是存在HDF文件属性中。比如MODIS的L1B辐亮度产品，每个波段有两个属性，一个叫radiance_scales，一个叫radiance_offsets，读取方式变成查属性表而不是查文本。HDF属性读取用Python的h5py或pyhdf也很快，几十行能搞定，但这里先不展开，后面提到的代码逻辑是一样的。

2.3 传统公式里的日地距离与ESUN：两种反射率算法对比

如果拿到的是旧一代卫星数据，元数据里没有直读式反射率系数，就需要走经典公式计算表观反射率：

ρ = π * L * d² / (ESUN * cos θz)

其中L是辐亮度，d是日地距离修正因子，ESUN是波段平均太阳辐照度，θz是太阳天顶角。d不是一个天文单位距离的绝对值，而是一个以年为周期的比例系数，常见近似公式是：

d = 1 - 0.01673 * cos(0.0172 * (DOY - 2))

DOY是影像获取日期的年积日，1月1日为第1天。这个近似公式的精度已经足够，误差在0.1%以内，完全不影响后续反演。ESUN则是各波段的太阳等效辐照度，在数据产品的官方文档里能查到，不同卫星、不同波段的ESUN都不一样，不能跨数据源混用。

把直读式和经典式放在一起对比，各自的适用情况就很清楚了。

我一般优先用直读式系数，因为它少传两个参数，少两个出错点。经典式保留在代码里做兜底，处理老数据时再用。

3. 用Python跑通最小定标流程：读DN、算辐亮度、写回GeoTIFF

3.1 最小命令：rasterio读DN并完成辐射定标

先给一个能直接跑的最小示例。假设你已经从影像产品中拿到了一个波段TIF和对应的定标系数，第一步就是把DN读进来做乘加运算。

import rasterio import numpy as np with rasterio.open("scene_band1.tif") as src: dn = src.read(1).astype("float64") profile = src.profile gain = 0.0211 # 示例值，必须替换成元数据里的RADIANCE_MULT_BAND_x offset = -5.2311 # 示例值，必须替换成元数据里的RADIANCE_ADD_BAND_x radiance = gain * dn + offset print("radiance range:", radiance.min(), radiance.max())

这段代码取的是单波段。点开TIF时注意波段序号的起始位置，rasterio的read(1)读的是第一个波段，对应元数据里的BAND_1。如果影像的波段序号从0开始的产品，就必须先确认波段映射关系，否则定标完的波段和实际物理波段对不上。这是最容易在批量处理时翻车的细节。

把DN转成float64是刻意的。很多光学卫星影像是uint16存储，DN范围能到0到65535，增益乘进去之后数值可能跨好几个数量级，用uint16中间存储会溢出，用float32虽然精度够，但在极端高值和低值像元混合的场景下有舍入风险。定标是逐像元线性变换，我用float64跑一遍不心疼，反正后面还会转回float32写出。

3.2 自动解析MTL文件：别再手动抄系数

手动把gain和offset填进代码，能跑通一个波段，但处理一景影像十几个波段时效率太低，而且很容易把BAND_5的系数贴到BAND_6上。更好的办法是直接从MTL文件里读取系数。

import re def parse_mtl(mtl_path): params = {} with open(mtl_path, "r", encoding="utf-8") as f: for line in f: line = line.strip() if "=" not in line: continue key, value = line.split("=", 1) params[key.strip()] = value.strip() return params mtl = parse_mtl("scene_MTL.txt") band = 1 refl_mult = float(mtl[f"REFLECTANCE_MULT_BAND_{band}"]) refl_add = float(mtl[f"REFLECTANCE_ADD_BAND_{band}"]) sun_elev = float(mtl["SUN_ELEVATION"]) # 直读式反射率公式 toa = (refl_mult * dn + refl_add) / np.sin(np.radians(sun_elev))

这个parse_mtl函数看起来很简陋，但对MTL这类“缩进由空格构成、分组用GROUP/END_GROUP包裹”的文本，逐行按等号切分反而比正则更稳。正则匹配如果要处理GROUP嵌套、键名重复，容易写出体积膨胀且难维护的模式。逐行split的唯一注意点是，某些字段值可能带引号，如果引号存在，需要再处理一次：

def clean_value(v): v = v.strip() if v.startswith('"') and v.endswith('"'): return v[1:-1] return v

解析完成后，mtl就是一个以字段名为键的字典，接下来无论你是做单波段测试还是批量处理，都可以从同一个字典取系数。取值时用f-string拼波段号，循环里非常方便，而且波段号写错会在第一次运行时直接报KeyError，比静默读到错误系数好查得多。

3.3 批量定标多波段并写回GeoTIFF：保留地理信息还能压缩

有了解析函数，下一步就是批量处理一景影像的所有波段，并把结果写成新的GeoTIFF。如果还是用老办法——先转成numpy数组计算，再用rasterio新建文件写CRS和transform，很容易把地理信息写丢。正确姿势是直接从输入文件的profile复制，再修改必要字段。

import rasterio import numpy as np def calibrate_toa_band(src_path, dst_path, refl_mult, refl_add, sun_elev): with rasterio.open(src_path) as src: dn = src.read(1).astype("float32") profile = src.profile.copy() toa = (refl_mult * dn + refl_add) / np.sin(np.radians(sun_elev)) profile.update( dtype="float32", nodata=-9999.0, compress="deflate", predictor=3, ) with rasterio.open(dst_path, "w", **profile) as dst: dst.write(toa, 1)

逐行说三个参数。nodata设置成-9999，是因为定标后有效反射率理论上在-0.1到1.2之间，-9999足够远离有效范围，后续统计时一眼能识别。compress="deflate"是GeoTIFF里通用性最好的压缩方式，不要用LZW，defrate在大多数平台兼容性更好而且压缩比相当。predictor=3是浮点数据专用的预测器，它按浮点数的指数部分做差分预测，对反射率这种取值相对连续的浮点数组，文件体积能比不设预测器再小20%到40%。

批量处理时，我会建一个波段映射表，把输入文件名列表和波段序号对齐：

band_files = [ "scene_B2.TIF", "scene_B3.TIF", "scene_B4.TIF", "scene_B5.TIF", ] for idx, src_path in enumerate(band_files, start=1): band_idx = idx + 1 # 根据产品实际波段号调整 mult = float(mtl[f"REFLECTANCE_MULT_BAND_{band_idx}"]) add = float(mtl[f"REFLECTANCE_ADD_BAND_{band_idx}"]) dst_path = f"TOA_B{band_idx}.TIF" calibrate_toa_band(src_path, dst_path, mult, add, sun_elev) print("done:", dst_path)

这里最值得注意的就是波段映射。不同产品的文件名后缀与实际波段号的关系并不固定，有的用B2、B3、B4，有的用Band1、Band2。批量循环中如果不对应，程序不会报错，只会默默产出错误的反射率图。所以我常在循环里加一个断言：

assert refl_mult > 0 and abs(refl_add) < 1.0, "系数疑似读错"

多光谱反射率定标系数的增益通常是个很小的正值，偏置通常在零附近。如果读出来的是负数或上百，十有八九是波段号对应错了。

4. 定标路上翻车最多的地方：5个常见问题排查

4.1 定标结果全黑，是失败还是显示设置问题

现象：定标完的浮点影像在GIS软件里打开，整幅图是全黑的，看不到任何地物轮廓。

原因：反射率浮点影像的数值集中在0到0.5之间，而大多数软件默认按全范围线性拉伸显示。如果影像里有一个异常的极大值，比如云顶反射率接近1，或者有一个负的无效值-9999，显示系统会把有效值压缩到极窄的灰度区间，看起来就是全黑。

解决：先别怀疑定标失败，用统计量确认数值分布。

with rasterio.open("TOA_B4.TIF") as src: arr = src.read(1) mask = arr < -1000 # 过滤nodata valid = arr[~mask] print(np.nanpercentile(valid, [2, 50, 98]))

如果2%分位数和98%分位数落在0到0.6之间，定标结果基本正常。显示时手动拉伸到这两个分位数，或者用软件的2%线性拉伸就能正常看。我习惯把定标代码里顺手输出一组统计量，处理完第一眼就能判断有没有低级错误。

4.2 辐亮度大了1000倍，先检查单位是μm还是nm

现象：定标出来的辐亮度数值比经验估算大了1000倍，或者某个波段的辐射值明显不合理。

原因：辐亮度的单位在不同数据产品里不完全一致。光学遥感常用的单位是W/(m²·sr·μm)，但部分传感器产品在文档里给出的是W/(m²·sr·nm)，两者差了1000倍。如果按μm的单位读入数据，却套用nm系数，整幅图所有像元都会放大1000倍。这个错误不报错、不影响显示拉伸，但下游反演模型会彻底崩溃。

解决：拿到数据先看波段中心波长的量级。中心波长在几百nm量级的传感器，辐亮度如果写成nm单位，数值会非常小，接近0.001级别；μm单位的数值则在1到10之间。根据量级判断该用哪套系数。更稳妥的做法是定标后输出一个中位数，水体在红波段辐亮度应该在个位数水平，如果中位数变成几千，几乎可以断定单位错位。

4.3 天顶角与高度角混用，反射率直接超过1

现象：计算表观反射率后，大量像元值超过1，甚至出现几十倍的极端值。

原因：太阳参数有两种表示方式，太阳天顶角是垂直方向与太阳光线的夹角，太阳高度角是水平面与太阳光线的夹角，两者互余，差90度。很多元数据文件里写的是SUN_ELEVATION，也就是高度角，但计算反射率时有人会把它当成天顶角，用余弦去算。天顶角50度对应的余弦是0.64，高度角50度对应的正弦是0.77，混用之后反射率会被放大或缩小约20%。更极端的案例是把高度角直接代入余弦，天顶角代入正弦，结果反射率超过2。

解决：明确记住一句口诀——元数据给的是高度角就用正弦，给的是天顶角就用余弦。下面的代码同时支持两种模式：

def to_toa_reflectance(dn, refl_mult, refl_add, sun_elevation=None, sun_zenith=None): if sun_elevation is not None: denominator = np.sin(np.radians(sun_elevation)) elif sun_zenith is not None: denominator = np.cos(np.radians(sun_zenith)) else: raise ValueError("需要提供太阳高度角或天顶角") return (refl_mult * dn + refl_add) / denominator

额外提醒：角度转弧度这步不能省，numpy的三角函数默认按弧度计算，直接传入50度，结果会非常荒谬。np.radians就是用来做这个转的，很多人漏掉它。

4.4 多景影像共用一份定标系数，时序分析数据全部作废

现象：批量处理一个月内的多景影像时，个别景的NDVI时序出现突变，但原始影像直方图看起来正常。

原因：光学卫星的定标系数不是永久不变的。卫星发射后会定期更新绝对定标参数，同一天的L1产品可能因为元数据版本更新，增益系数发生微小变化。如果整个批次处理时只从某一景的MTL里读了一次系数，全批次共用了，那其他几景影像都会带上系统偏差。偏差小的时候看不出来，一旦做NDVI阈值分类或变化检测，边界就出来了。

解决：每一景影像单独解析它自己的MTL文件，哪怕来自同一目录。批量处理时，脚本里的循环应该是“一景一读”，而不是“一批共读一份”：

import pathlib scene_list = pathlib.Path("scenes").glob("*_MTL.txt") for mtl_path in scene_list: mtl = parse_mtl(str(mtl_path)) sun_elev = float(mtl["SUN_ELEVATION"]) # 再按波段处理该景的TIF，文件前缀与MTL同名

我一般还会把处理日期和MTL文件路径记录到日志文件里，这样复查数据的时候能快速定位是哪一份MTL产生了问题。

4.5 写入float影像时把nodata当成有效值

现象：定标后的影像在查询统计量时，最小值总是-9999，或者融合分析时边缘出现大片异常低值。

原因：原始影像的边界外像元本身是无效值，DN可能是0或者某个填充值。定标公式只做乘加和除法，不会自动识别这些无效值，结果边缘外像元也参与了运算，变成负值或极大值。

解决：把掩膜信息从输入影像带出来，在输出前标记无效值。

with rasterio.open(src_path) as src: dn = src.read(1).astype("float32") mask = src.read_masks(1) # 有效像元为255，无效为0 profile = src.profile.copy() toa = (mult * dn + add) / np.sin(np.radians(sun_elev)) toa = np.where(mask == 0, -9999.0, toa) profile.update(dtype="float32", nodata=-9999.0, compress="deflate", predictor=3) with rasterio.open(dst_path, "w", **profile) as dst: dst.write(toa, 1)

read_masks(1)读取的是波段的有效区域掩膜，返回0到255的数组，只有像元有效时是255。np.where把无效位置替换成nodata。这是最容易忽略的一步，尤其当原始影像四周有黑边时，不定标看不出来，定标后黑边变成巨大的负值，会让后续所有统计和可视化都受到污染。

5. 定标精度验证：没有地面真值也能做的三个检查

5.1 用官方示例像元回算，验证公式与系数无误

大多数公开遥感产品的文档里，会给出至少一个示例像元的DN值以及对应的辐亮度或反射率。这是一种非常廉价的验证方式。先把你的定标函数写成一个纯函数，然后把官方示例DN代入：

def dn_to_toa(dn, refl_mult, refl_add, sun_elev): return (refl_mult * dn + refl_add) / np.sin(np.radians(sun_elev)) official_dn = 8200 official_toa = 0.2513 # 假设文档给出的参考值 my_toa = dn_to_toa(official_dn, refl_mult, refl_add, sun_elev) print("relative error:", abs(my_toa - official_toa) / official_toa * 100)

相对误差在0.1%以内，说明公式和系数都没问题。如果误差超过1%，先查有没有把Mρ和ML写反，再查太阳角度单位。这个验证步骤最好写进批处理脚本的开头，每次换数据源先跑一遍，能省掉后面排查的半天时间。

5.2 用典型地物的反射率范围做定性校验

没有官方示例也没关系，遥感图像灰度直方图里隐含了地物信息。定标后的表观反射率有个经验范围：清洁深水体在近红外波段通常低于0.05，植被在近红外波段通常高于0.3，裸土在红波段大约0.15到0.3，云和雪可以冲到0.8以上。

拿到定标影像后，我习惯直接统计几个已知地物点的反射率，比如一片湖泊、一片农田、一个裸土区域：

import rasterio points = { "water": (420000, 3260000), "crop": (421000, 3261000), "bare": (422000, 3262000), } with rasterio.open("TOA_B4.TIF") as src: for name, (x, y) in points.items(): row, col = src.index(x, y) value = src.read(1, window=((row, row + 1), (col, col + 1)))[0, 0] print(name, value)

近红外波段水体反射率如果超过0.1，或者裸土接近0，说明定标有问题。用这种方式做不了精确验证，但能在十分钟内筛掉80%的低级错误。注意这些阈值适用于表观反射率，如果做了大气校正得到地表反射率，水体值会更低，裸土上限也有变化。

5.3 多期影像一致性检验，定位定标异常

时序分析场景下，还有一个更实用的验证逻辑：同一地物在多期影像里的表观反射率虽然会有波动，但不会出现成倍变化。如果同一块农田在A影像里近红外反射率是0.35，B影像里变成0.8，排除极端天气和作物收割因素后，第一个怀疑对象就是定标系数。

def check_consistency(toa_a, toa_b, thresh=0.5): ratio = toa_b / (toa_a + 1e-6) bad = np.abs(ratio - 1.0) > thresh return np.mean(bad) * 100 # 返回异常像元比例

这个检验函数适合在两期影像覆盖重叠区域时跑，异常比例超过5%就要回去查太阳角度或者系数是否取自正确的波段。我自己处理几十年土地变化数据时，这套方法抓出过两期MTL的SUN_ELEVATION字段读取顺序错误的问题，比肉眼盯直方图快得多。

6. 把这套流程封装成命令行工具：让下一次定标不再重新踩坑

我习惯把解析MTL、定标、写GeoTIFF封装进一个独立脚本，用命令行参数驱动。这样每次处理新数据不用再打开Jupyter手动改路径，批处理和多景并行也能复用同一套逻辑。一个精简版的结构如下：

# toa_calib_cli.py import argparse import json import pathlib import rasterio import numpy as np def parse_mtl(mtl_path): params = {} with open(mtl_path, "r", encoding="utf-8") as f: for line in f: line = line.strip() if "=" not in line: continue key, value = line.split("=", 1) params[key.strip()] = value.strip() return params def calibrate_file(src_path, dst_path, mult, add, sun_elev): with rasterio.open(src_path) as src: dn = src.read(1).astype("float32") mask = src.read_masks(1) profile = src.profile.copy() toa = (mult * dn + add) / np.sin(np.radians(sun_elev)) toa = np.where(mask == 0, -9999.0, toa) profile.update(dtype="float32", nodata=-9999.0, compress="deflate", predictor=3) with rasterio.open(dst_path, "w", **profile) as dst: dst.write(toa, 1) def main(): parser = argparse.ArgumentParser(description="TOA绝对辐射定标工具") parser.add_argument("--mtl", required=True, help="MTL/XML元数据文件路径") parser.add_argument("--band", type=int, required=True, help="波段号") parser.add_argument("--src", required=True, help="输入TIF路径") parser.add_argument("--dst", required=True, help="输出TIF路径") args = parser.parse_args() mtl = parse_mtl(args.mtl) mult = float(mtl[f"REFLECTANCE_MULT_BAND_{args.band}"]) add = float(mtl[f"REFLECTANCE_ADD_BAND_{args.band}"]) sun_elev = float(mtl["SUN_ELEVATION"]) calibrate_file(args.src, args.dst, mult, add, sun_elev) with open(pathlib.Path(args.dst).with_suffix(".json"), "w", encoding="utf-8") as f: json.dump({ "src": args.src, "band": args.band, "mult": mult, "add": add, "sun_elev": sun_elev, }, f, indent=2, ensure_ascii=False) if __name__ == "__main__": main()

命令行运行方式：

python toa_calib_cli.py --mtl scene_MTL.txt --band 4 --src scene_B4.TIF --dst TOA_B4.TIF

脚本会在输出TIF旁边生成一个同名JSON，记录这次定标用的系数和太阳高度角。此举不是为了好看，是为了下次查数据时不用再翻原始MTL——看一眼JSON就知道某幅影像的数值是怎么来的。我在这套小工具上修修补补了两年，翻车最多的地方从来不是numpy，而是系数有没有取自对应波段、太阳角度有没有搞错。现在把关键参数落在JSON里，每次处理完顺手核对一眼，基本没再出过低级错误。希望帮到你。

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