简介：遥感+Python绝对辐射定标资源面向遥感初学者与地理信息开发人员，演示如何用GDAL库将影像原始数字值转为具有物理意义的辐射值，涉及数据准备、定标参数读取、辐射校正计算、大气校正及结果输出等流程。压缩包共9个文件，16KB，包含1个Python脚本、1份JSON定标参数文件、5个XML工程配置，以及项目配置与版本管理文件，结构精简便于直接查看和复用。已有902人学习，适合对照代码理解辐射定标在环境监测、资源调查等场景中的应用。通过学习可掌握基于GDAL的波段读取、参数加载与写回方法，并结合RadiometricCorrectionParameter.json体会定标参数的组织形式，为后续定量遥感分析打下基础。脚本中预留apply_calibration与create_output_dataset函数结构，便于按传感器参数二次开发，同时涉及大气校正的扩展思路，可帮助理解从DN值到地表反射率的完整链路。
1. 绝对辐射定标到底在解决什么问题：别再拿DN值硬算NDVI
某次做多时相植被监测，三个月内的影像直接算NDVI，结果同一个地块的数值毫无规律，一会儿0.6一会儿0.2，第一反应是大气影响，处理半天发现根子上是没做绝对辐射定标。遥感影像L1级产品里存的大多是DN值，它和真实辐射亮度之间隔着增益和偏移两层系数，不同时相、不同传感器的系数都不相同。绝对辐射定标就是用python把DN还原成辐射亮度，再校正太阳角度和日地距离得到表观反射率，让跨时相影像真正可比。这篇笔记从元数据解析讲起，把每一步计算、参数和容易翻车的地方过一遍，适合正在做定量反演、时序监测的从业者照着复现。
2. 辐射定标之前：先把数据形态和元数据摸清楚
拿到影像就急着套公式，是很多人在绝对辐射定标上交的第一笔学费。定标本身只有两步线性运算，但前置工作决定这两步能不能落对。你需要先弄清楚三个问题：最终要什么物理量、定标系数存在哪、手头产品是不是已经被定标过。这一章把这三件事一次讲完。
2.1 先选通道：辐射亮度还是表观反射率
绝对辐射定标不是一个动作，而是分层递进：第一层把DN转成大气层顶辐射亮度，第二层把辐射亮度结合太阳参数转成表观反射率。很多人不区分这两层，结果要么多算了一步做了无谓校正，要么少算了一步导致跨时相对比失真。
我一般会先问一句：你要比较的是多时相物理量，还是只用单景影像做相对分析？单景影像做分类、分割、目标检测这类任务，直接用DN跑特征通常也能出结果，定标不是必需品。但凡是跨时相、跨传感器的定量活，比如NDVI时序、水体面积年变化、城市热岛对比，就必须把绝对辐射定标做完整。
判断目标后，选择做到哪一层可以参照这个表：
| 应用场景 | 需要做到哪一层 | 原因 |
|---|---|---|
| 单景影像地物分类 | 可不做定标 | 特征相对差异已足够 |
| 热红外温度反演 | DN转辐射亮度 | 反演模型输入为辐射亮度或亮度温度 |
| 植被/水体指数时序分析 | 到表观反射率 | 需消除太阳高度角与日地距离影响 |
| 多传感器一致性拼接 | 到表观反射率并再交叉校正 | 各传感器光谱响应不同 |
选错了层的典型表现是：热红外反演时多做了日地距离校正，温度结果整体偏移；时序NDVI对比时只做了辐射定标没做太阳角度校正，冬季影像反射率明显虚低，看起来像植被在冬天集体退化。
2.2 从元数据里读定标系数：MTL文本解析与字段命名
定标系数的存放位置随传感器而异。常见形态有三种：某个光学卫星的L1级产品带一个MTL文本文件，按GROUP分组列出每个波段的定标乘子和加偏移；另外一些平台给XML格式的元数据，定标系数嵌在radiometric节点下；还有的国产卫星直接交付一张Excel参数表，需要手动核对波段顺序。
MTL文本本质是缩进分组的Key=Value文件，不需要依赖第三方解析库，用ast内置库加一个缩进栈就能完整解析。下面这段代码处理常见MTL结构：
import re from ast import literal_eval def parse_mtl(path): """把MTL缩进文本解析成嵌套字典，便于按组取值""" with open(path, 'r', encoding='utf-8') as f: lines = [ln.strip() for ln in f if '=' in ln or ln.strip().startswith(('GROUP', 'END_GROUP'))] root = {} stack = [root] for ln in lines: if ln.startswith('GROUP'): name = ln.split('=')[1].strip() new_dict = {} stack[-1][name] = new_dict stack.append(new_dict) elif ln.startswith('END_GROUP'): stack.pop() else: key, val = ln.split('=', 1) try: stack[-1][key.strip()] = literal_eval(val.strip()) except Exception: stack[-1][key.strip()] = val.strip() return root meta = parse_mtl('scene_MTL.txt') rescaling = meta['RADIOMETRIC_RESCALING'] gain = rescaling['RADIANCE_MULT_BAND_4'] offset = rescaling['RADIANCE_ADD_BAND_4'] print('增益:', gain, '偏移:', offset)这段代码把MTL里所有GROUP块压成嵌套字典，之后不管是取哪个波段哪个参数，都走meta['组名']['字段名']这条路。literal_eval用来把数值字符串转型成浮点数或整数，比float()更稳，能同时处理整数、小数和科学计数法。
注意字段名在不同传感器里差别很大。比较常见的命名有RADIANCE_MULT_BAND_x加RADIANCE_ADD_BAND_x，也有叫RADIANCE_MULT与RADIANCE_BIAS的，还有直接叫GAIN_x和OFFSET_x的。拿到新数据第一件事是把字段名打印出来，别套旧模板，因为命名差异造成的取错系数，比计算本身更容易出问题。
2.3 确认产品是不是已经做过定标：查量级与元数据字段
有些L1级产品交付时就是TOA反射率，元数据里直接给出REFLECTANCE_MULT这样的字段，或者像素值存成0到1之间的浮点。如果这样的数据再走一遍绝对辐射定标，等于把反射率当DN又乘了一遍增益加偏移，数值会被二次缩放，后面所有反演都失真。
判断是否已定标，我习惯用三步快速排查。第一步看像素值量级：DN通常是0到65535的整数，TOA反射率要么是0到1的小数，要么是0到10000的整数（配合scale factor使用）。第二步查元数据字段：出现REFLECTANCE_MULT、REFLECTANCE_ADD、SCALE_FACTOR之类的字样，说明产品已经做过定标。第三步打开影像做波段统计，均值落在0.1到0.5之间且分布平滑，基本可断定是反射率产品。
提示：拿到新数据先写三行代码打印各波段min、max、mean，再打开元数据grep一遍关键字段。这一步比任何公式都重要。
3. 用python把DN转成辐射亮度：numpy计算与数据类型
确认完产品形态，进入绝对辐射定标第一层：DN转辐射亮度。公式极其简单，但数据类型和水系数的组织方式决定了你算出来的是干净结果还是垃圾结果。这一章给出的实现同时照顾了小影像的直观验证和大影像的内存控制。
3.1 定标方程与最小实现：先乘后加别搞反
辐射亮度定标的标准方程是 L = gain × DN + offset ，其中L是大气层顶辐射亮度，单位是 W/m²/sr/μm 。gain来自元数据，offset通常为负值，用于扣除传感器的暗电流基线。
一个最小实现只需要四行：
import numpy as np def dn_to_radiance(dn, gain, offset): """DN转大气层顶辐射亮度。 dn: numpy数组，取自影像单波段 gain: 该波段的定标增益，浮点数 offset: 该波段的定标偏移，浮点数 """ dn_float = dn.astype(np.float32) radiance = gain * dn_float + offset return radiance先astype再乘是有原因的：如果DN保持uint16，乘上增益后numpy可能溢出，而且溢出不会报错，只会给你一个莫名其妙的整数。float32在这个场景下精度足够，又不至于像float64那样让内存翻倍。
参数说明：gain和offset必须是该波段自己的系数，不能拿蓝波段的增益去乘红波段的DN；运算必须先乘后加，因为分配律在这里不成立——先加后乘等于把offset也乘了一遍gain，结果整体偏移了一个gain倍的常数。这个错误在最终反射率上表现为全图系统性偏移，很难一眼看出来。
3.2 数据类型的坑：uint16溢出与负偏移丢失
运行时最常见的两个翻车现场都与数据类型有关。第一个是乘法溢出，前面提过；第二个是负偏移导致的负值被钳制。DN转辐射亮度时，offset是负值，暗目标（水体、阴影）算出来可能是负数，这是物理上正常的暗电流扣除结果。但如果定标后回存成uint16，numpy会把负数钳到0，暗目标细节全部丢失，后续做水体信息提取时这部分像素就废了。
处理方式是中间结果一律用float32的GeoTIFF保存，不要为了省空间回存成整数。实测中，float32的GeoTIFF体积大约是uint16的两倍，但换来的是后续无需再做类型转换，也避免了再次读取时数值被截断。
下面这个函数把读取、定标、写出连贯起来，同时保留了分块能力：
import rasterio def apply_rescaling(src_path, dst_path, gain, offset): """对整景影像做辐射定标并写出float32 GeoTIFF。""" with rasterio.open(src_path) as src: profile = src.profile.copy() profile.update(dtype='float32') with rasterio.open(dst_path, 'w', **profile) as dst: for block_index, window in src.block_windows(1): dn = src.read(1, window=window).astype(np.float32) radiance = gain * dn + offset dst.write(radiance, 1, window=window)block_windows按波段自动遍历所有分块，每块独立完成读、算、写三步，内存占用恒定，不会随着影像尺寸增长。对十几个GB的大面积影像，这个写法比一次性读入全图要稳得多。如果只需要测试单个波段，也可以手动构造一个Window传给读取函数，但分块循环才是处理生产的常态。
3.3 辐射亮度结果怎么验证：先盯暗目标再看亮目标
定标完别急着做反射率，先抽几个点验证辐射亮度在不在合理区间。做法是找影像里一块已知的水体或阴影，读取定标后的值，看是否接近0附近且允许存在小幅负值；再找一块裸土或混凝土屋顶，看是否落在该传感器常规辐射亮度范围内。
如果水体区域全部是正数且数值偏大，说明offset符号取反了，暗电流没有被扣除，结果整体被抬升。如果整景影像出现大量负值且亮目标数值过低，优先检查gain和offset是不是乘反了位置——确实有人把公式写成dn = gain × L + offset过。这个验证步骤耗时不超过两分钟，但能拦住后续整条反演链路带着系统性偏差跑完。
4. 从辐射亮度到表观反射率：太阳辐照度、天顶角与日地距离校正
辐射亮度是传感器口的物理量，但它同时受太阳高度角和日地距离影响：盛夏正午和严冬清晨，同一个地物的辐射亮度差异很大，可地物本身没有变化。表观反射率通过归一化太阳入射条件消除这部分差异，让多时相影像真正可比较。
4.1 表观反射率公式与参数查表：ESUN到底该用哪个值
TOA反射率的标准公式是 ρ = π × L × d² ／ (ESUN × cos(SZA))。这里L是上一章算出的辐射亮度，d是日地距离（天文单位），ESUN是大气层顶等效太阳辐照度，SZA是太阳天顶角。
ESUN的取值是最容易踩坑的地方，因为不同传感器、不同波段的ESUN差异很大。常见某传感器的典型值大致是这样的量级：
| 波段 | ESUN（W/m²/μm） |
|---|---|
| 蓝 | 约1930 |
| 绿 | 约1826 |
| 红 | 约1514 |
| 近红外 | 约1036 |
这张表只适合用来做量级判断，真正生产时必须用你自己传感器头文件里给出的ESUN表。优先读取元数据中的SOLAR_IRRADIANCE相关字段，元数据里没有再去查该传感器官方发布的技术参数表。特别强调：ESUN的单位是W/m²/μm时，必须与辐射亮度的波长单位一致，若某传感器用nm给出增益系数，需要先统一量纲再做除法，差1000倍直接导致结果完蛋。
4.2 日地距离d的计算：从年积日到天文单位的近似公式
日地距离d随着地球绕太阳公转呈周期性变化，影响幅度约±3.4%。对表观反射率来说，忽略这个校正会产生约3%的系统误差，时序分析里足以被识别成虚假的季节趋势。
用python计算d不需要星历表，一个基于轨道离心率的近似公式就够用：
import math def earth_sun_distance(yday): """根据年积日计算日地距离，返回天文单位。 yday: 影像成像日期的年积日，1月1日为1 """ angle = 2 * math.pi * (yday - 1) / 365.0 d = 1 - 0.01672 * math.cos(angle) - 0.00029 * math.sin(angle) return d distance = earth_sun_distance(150) print('日地距离:', distance)这个近似公式的误差在10⁻³量级，换算成反射率的偏差不超过0.1%，完全满足常规绝对辐射定标需求。年积日的计算很简单：从影像元数据里读成像日期，然后用date.timetuple().tm_yday直接拿到day of year。如果元数据里直接给了DAY_OF_YEAR字段，优先用元数据值。
注意：日地距离校正只作用于表观反射率这一层，热红外辐射亮度定标不需要做，因为热红外能量来自地物自身发射而非太阳反射。
4.3 合并成完整定标函数：DN直接到TOA反射率
把辐射定标、日地距离、太阳角度串成一个函数，调用方只需要传DN数组和四个参数，不用关心内部换算细节：
import numpy as np import math def dn_to_toa_reflectance(dn, gain, offset, esun, sza_deg, yday): """DN直接转表观反射率。 dn: numpy数组，uint16原始DN值 gain/offset: 该波段辐射定标系数 esun: 该波段大气层顶等效太阳辐照度，W/m^2/um sza_deg: 太阳天顶角，度 yday: 年积日 """ dn_float = dn.astype(np.float32) radiance = gain * dn_float + offset sza_rad = math.radians(sza_deg) d = earth_sun_distance(yday) reflectance = math.pi * radiance * d**2 / (esun * math.cos(sza_rad)) return reflectance这段代码把两层定标放进一个入口，适合批量处理。参数说明：sza_deg是太阳天顶角，不是太阳高度角，两者互余，元数据里如果给的是高度角EL，要先做90 − EL的转换；math.cos要求输入弧度制，漏掉radians这一步是高频bug，后面避坑章专门讲。
这个函数设计上故意没有dtype隐式转换，DN输入保持原样，输出由numpy自动推断为float32，同时避免了在函数内部重复astype引入的额外拷贝。
4.4 批处理多时相：定标函数的复用要点
多时相绝对辐射定标的本质就是循环调用dn_to_toa_reflectance，但有几个细节值得注意。先把各景影像的元数据集中读一遍，把每景的增益、偏移、太阳天顶角、年积日整理成一个列表，再统一循环，而不是在循环里反复解析同一个MTL文件。
另外，多时相输出文件的命名建议带上时相标签，如toa_b4_20250115.tif，避免定时后无法对应原始影像。批量处理时要注意每景影像的太阳天顶角都在变化，不能为了省事固定用一个角度值，春秋两季差别最大时反射率能差出5%以上。这段批处理逻辑稳定跑通后，后续所有时序分析的数据准备就自动化了。
5. 绝对辐射定标避坑实录：5个高频问题与排查路径
绝对辐射定标计算不难，难的是数据来源千奇百怪：不同传感器、不同产品级别、不同元数据风格。以下五个问题是平时被问得最多的，全部按现象、原因、解决三步记录，可直接对照排查。
5.1 反射率大于1或分布出现亮斑
现象：完成定标后的影像上，部分像素反射率超过1，甚至到1.5以上，直方图右尾拖得极长。
原因：最常见是ESUN取错，用了一个临近波段的数值代替目标波段；其次是太阳天顶角用的角度不对，cos(θ)偏小导致整体放大；第三是增益系数版本陈旧，和实际传感器定标参数不符。
解决：按波段中心波长逐一核对ESUN，不要凭记忆填；再查太阳角度字段，确认用的是天顶角而非高度角。如果所有像素都呈固定倍数偏高，基本可断定是ESUN或角度问题；如果只有高亮目标异常，检查是否有云或耀斑未掩膜。
5.2 近红外波段反射率反常识地低
现象：植被区域的近红外反射率低于红波段，NDVI算出负值，和常识完全相反。
原因：波段顺序对错了。很多传感器的波段编号和自然光谱顺序不一致，比如波段4可能是红波段而不是近红外；或者读取增益时把某个可见光波段的ESUN套在了近红外上。
解决：先打印元数据里的中心波长列表，按波长重新映射波段顺序，再重新做定标。不要相信文件名里的band编号，要以元数据中的中心波长为唯一依据。这个错法在最终结果上非常隐蔽，因为图像看起来正常，只有数值不合逻辑。
5.3 结果整体差1000倍：um与nm单位错配
现象：全图反射率在0.000x量级，影像几乎全黑，只有最亮的目标隐约可见。
原因：辐射亮度单位是W/m²/sr/μm，但某些传感器的增益偏移系数是按W/m²/sr/nm给出的。直接用nm单位的增益算出的辐射亮度比真实值大1000倍，再除ESUN，反射率就被压缩到原来的千分之一。
解决：统一量纲后再算。两个方案选一个：把增益偏移除以1000换算成μm基线，或者把ESUN乘以1000换算成nm基线，保证分子分母单位一致。更稳的做法是在定标函数里加一个unit参数，默认um，遇到nm时显式换算并打印日志留档。
5.4 太阳高度角被当成天顶角用
现象：冬季影像反射率整体虚高，越靠近两极越严重，同一地块夏季冬季对比异常。
原因：元数据里给的是太阳高度角SUN_ELEVATION，比如冬季只有30度，公式里需要的是天顶角60度。直接把30度塞进cos函数，cos(30°)=0.866，而正确cos(60°)=0.5，反射率被放大了1.7倍还多。
解决：写一个防御性检查，读取元数据字段名，如果字段是SUN_ELEVATION，就自动转成天顶角再传入定标函数；同时输出两个角度到日志，人眼扫一眼确认量级合理。这个坑最危险的地方在于结果看起来和谐，只是数值偏高，容易带着错误一路做完反演。
5.5 拿到的产品本身就是TOA反射率
现象：定标后的影像数值全挤在0附近，直方图严重左偏，NDVI等指数完全失真。
原因：L1级产品交付时有可能已经做了TOA反射率处理，像素本身就在0到1之间或0到10000之间。再套一次gain乘加，等于把反射率二次缩放，数据被压扁或平移。
解决：处理前看元数据里有没有REFLECTANCE_MULT、REFLECTANCE_ADD或SCALE_FACTOR字段，有就直接认定该产品已经过反射率定标，跳过gain/offset步骤，只做日地距离和太阳角度校正，或者直接用产品自带的scale factor解算反射率。这属于流程选择问题，和数据本身无关，但造成的后果比任何参数错误都严重。
6. 用伪不变目标验证定标结果：一个值得养成的定标后检查习惯
定标完不知道结果对不对，等于没定标。我最常用的验证方法不是对着参考产品做绝对值比对，而是选一个伪不变目标看多时相一致性：裸岩、沙漠、深水体这类短期内光学性质稳定的地物，在不同时相的TOA反射率应该基本一致，变异系数在5%以内。
操作分三步。第一步在基准影像上勾一个纯像元区域，避开云、云影、地形阴影和像元混合区域，面积不用大，几百个像元就够。第二步把多个时相定标好的影像裁剪到同一个区域，逐波段计算均值和标准差。第三步算变异系数，也就是标准差除以均值，看是否落在合理区间。
一个精简的验证代码如下：
import numpy as np def evaluate_temporal_stability(stack, roi_mask): """统计伪不变目标多时相反射率稳定性。 stack: (时相数, 行, 列) 的反射率堆叠 roi_mask: 与单景影像同形状的布尔掩膜，True为伪不变目标 """ values = stack[:, roi_mask] # 提取所有时相、所有目标像元 means = values.mean(axis=1) stds = values.std(axis=1) cv = stds / (means + 1e-6) # 变异系数，加极小值防除零 return means, cv means, cv = evaluate_temporal_stability(toa_stack, roi_mask) print('各时相均值:', means) print('变异系数:', cv)这段代码的原理是：对每个时相先取ROI内所有像元反射率均值，再计算这组均值跨时相的波动程度。变异系数含义直观，小于0.05说明定标结果稳定，0.05到0.1之间要警惕是否有未校正的太阳角度残差，大于0.1就得回头查增益偏移和ESUN。这个验证不需要真实地物反射率先验知识，只要求目标在时间序列上光学稳定，实操里裸土场地和深水湖是最省心的选择。
我个人的习惯是把验证脚本和定标脚本放在一起，跑完定标后顺手执行一遍稳定性检查，结果存档到CSV里。多时相序列做多了，这些跨时相的变异系数曲线本身就能反映定标参数的漂移情况，再遇到异常值，不用从头排查，直接看哪一段序列的变异系数跳变就能定位到是哪个时相出了问题。这个习惯帮我挡掉了大量的后期返工，比任何理论推导都直接。
绝对辐射定标看起来就是两步线性运算，但几乎所有定量反演的误差都从这一步来。把影像定标做扎实，后面的反演等于在平地上走路。希望帮到你。
本文还有配套的精品资源，点击获取
