简介面向遥感图像处理初学者与需用 Python 完成辐射定标的从业者这份资料聚焦基于 GDAL 的绝对辐射定标实践帮助把原始数字值转换为具有物理意义的辐射量并衔接后续大气校正与分析流程。资源包共 9 个文件、约 16KB以 5 个 XML 配置文档、1 个 JSON 定标参数表、1 个 Python 主脚本及 IDE 模块与忽略配置文件组成JSON 用于承载分波段校正系数等参数Python 脚本演示读取影像、逐波段计算与写出结果XML 和项目文件则维护运行环境与工程结构。已有 903 人学习下载适合作为小型可运行示例参考。读者可从中获得按波段组织的定标代码骨架、参数外置思路、结果写出与元数据更新要点也能据此理解增益、暗电流等参数如何落到 GDAL Band 操作中便于迁移到自己的遥感数据与传感器参数上减少从公式到工程实现的试错成本。1. 从 DN 到物理量为什么绝对两个字值一次重跑同一景影像两个人跑出来的反射率差了近三倍最后追查下去不是公式写错而是定标系数表的表头写的是 L DN / gain脚本里却当成了乘。绝对辐射定标这件事公式只有一行失败点几乎全在公式之外DN 的填充值是多少、增益和偏移的符号怎么带、太阳天顶角用的是场景中心值还是逐像元角度、辐亮度量纲是 per μm 还是 per sr。遥感影像的绝对定标用 Python 加 GDAL 把 DN 一路推到辐亮度和大气顶反射率这条链路本身不长长的是参数核对。这份资源面向手里已经有 L1 级产品、需要把像素值变成物理量再去做反演或时序对比的人覆盖系数解析、分块读写、掩膜处理和精度复核。2. 定标公式拆开看增益、偏移和太阳几何的四个输入项定标链路一般是两段DN 到辐亮度辐亮度到大气顶反射率。前一段靠传感器定标系数后一段靠太阳几何和波段平均太阳辐照度。两段的输入项加起来其实只有四个gain、offset、太阳天顶角、日地距离或者说年积日。看起来简单但每个输入项都有至少两种常见写法写错哪一种结果都不报错只是数值悄悄偏掉。我一般的习惯是先在 5 个像素的小数组上把方向跑通再上全图。全图跑一遍动辄十几分钟方向错了等于白等。2.1 DN→辐亮度gain 和 offset 落在哪一步DN 是传感器记录的无量纲整数值辐亮度 L 是有量纲的物理量常见单位是 W·m⁻²·sr⁻¹·μm⁻¹。两者之间的线性关系是定标系数的全部内容但系数表的给出方式并不统一。多数产品文档给的是乘加形式也有不少给的是除的形式还有一些只给绝对定标系数 K公式写成 L DN / K。这三种写法在数值上能差一到两个数量级靠肉眼是看不出来的。import numpy as np # 五个像素先定方向跑全图之前必须过这一关 dn np.array([0, 100, 500, 1000, 4095], dtypenp.float32) gain, offset 0.0123, -0.45 # 写法 AL gain * DN offset L_a gain * dn offset k_div 0.0123 # 写法 BL DN / K L_b dn / k_div print(乘加:, np.round(L_a, 4)) print(除法:, np.round(L_b, 5))上面两组数的量级完全不同判断哪种才对的办法是查产品文档里给出的辐亮度动态范围或者找一景已知典型地物的影像看算出来的辐亮度落不落在文档声明的区间内。offset 的符号也要留意有些产品把偏移写成正数但公式里是减号等价于 L gain * DN - offset这种坑我在模拟项目里遇到过不止一次。量纲是第二个容易翻车的地方。辐亮度的分母里可能带 μm也可能带 nm分子可能是 W也可能是 μW。单位和数值是一起给的改单位就得改数值不能只改注释。我在脚本里习惯把 gain、offset 写进配置文件的同一行后面加一个 unit 字段强迫自己每次都读一眼。2.2 辐亮度→大气顶反射率太阳天顶角与日地距离辐亮度是传感器视角看到的量和成像时的太阳几何强相关同一地物在冬天和夏天成像辐亮度能差不少。所以做时序对比或者跨传感器对比之前通常要再推到大气顶反射率把太阳几何的影响归一化掉。公式是ρ π · L · d² / (ESUN · cos θs)其中 d 是日地距离单位天文单位ESUN 是波段平均太阳辐照度必须和 L 同量纲才能约掉θs 是太阳天顶角。这里有两个细节。第一很多产品的元数据给的是太阳高度角Sun Elevation天顶角等于 90 减去高度角这一步漏掉的话在低纬度地区误差还能忍到了中高纬度就是十几个百分点的偏差。第二d 不是一个常数一年之内在 0.983 到 1.017 之间变化平方之后约 3.4% 的浮动做精细对比时不能省。import numpy as np def sun_earth_distance(doy): 年积日 - 日地距离天文单位 AU常用的一阶近似 g 2.0 * np.pi * (doy - 1) / 365.0 return 1.00014 - 0.01671 * np.cos(g) - 0.00014 * np.cos(2 * g) def toa_reflectance(L, esun, sun_zenith_deg, doy, cos_min0.05): 辐亮度 - 大气顶反射率返回无量纲的 0~1 数值 cos_s np.cos(np.deg2rad(sun_zenith_deg)) cos_s np.maximum(cos_s, cos_min) # 挡住高纬低太阳高度角时的分母趋零 d sun_earth_distance(doy) return (np.pi * L * d * d) / (esun * cos_s)cos_min这个截断是必须的。太阳天顶角超过 85 度时 cos 值小于 0.09分母很小反射率会被放大到几十甚至上百看起来像是算法爆了。实际上那些像元要么是晨昏线附近的无效观测要么本来就该被掩掉。我的做法是先把 cos 截断到 0.05同时在掩膜里把天顶角大于 80 度的像元标记为无效两道保险。太阳天顶角还要确认是一景一个值还是逐像元一个值。小范围影像用场景中心值误差可以接受跨几个纬度的条带就必须用逐像元角度层很多 L1 产品会单独给一个角度波段别忘了读。2.3 元数据落在哪头文件字段和脚本变量的映射定标系数和太阳几何通常散落在产品自带的 XML 头文件里字段名各产品不同但语义就那么几个。我一般先做一个字段映射表把元数据名和脚本里的变量名对齐写进配置后面所有处理都从配置读不碰硬编码。元数据语义脚本变量单位 / 形式常见坑波段增益gain浮点可能是科学计数法表头写的是倒数形式波段偏移offset浮点可能带负号公式里是减号还是加号波段平均太阳辐照度esun与辐亮度同量纲单位是 W 还是 μW太阳高度角sun_elev度要转成天顶角成像日期doy年积日给的是 UTC 还是地方时有效 DN 范围dn_min / dn_max整数填充值常是 0 或 65535提示映射表里除了系数把 DN 的有效范围也一起记下来。后面写掩膜时直接引用不用每次翻文档。3. GDAL 读写骨架分块、掩膜与 dtype 的三个决定很多人写定标脚本第一版都是一次性ReadAsArray()读整景本地跑小图没问题换成几千行几千列的全色或者多光谱就内存溢出。GDAL 的分块读写不难难的是三个决定块开多大、掩膜什么时候套、输出用什么 dtype。这三个决定会直接影响脚本能不能跑完以及跑完的结果还能不能被下游软件正确读进去。3.1 打开数据集先确认五件事拿到一个影像文件别急着读数组先把这几项打出来看一眼。行列数决定分块方式数据类型决定中间计算要不要升到 float32nodata 决定掩膜条件投影和仿射变换决定输出能不能对齐原始影像。from osgeo import gdal gdal.UseExceptions() # 打开后出错会抛异常而不是悄悄返回 None ds gdal.Open(input_l1.tif, gdal.GA_ReadOnly) band ds.GetRasterBand(1) print(size:, ds.RasterXSize, ds.RasterYSize) print(dtype:, gdal.GetDataTypeName(band.DataType)) print(nodata:, band.GetNoDataValue()) print(gt[0], gt[1]:, ds.GetGeoTransform()[:2]) # 左上角 x 和像元宽度 print(proj:, ds.GetProjection()[:60])gdal.UseExceptions()是我每个脚本第一行都会加的。GDAL 默认行为是失败返回 None后面再调方法才报一个莫名其妙的 AttributeError排查起来很绕。打开异常之后文件路径写错、驱动不支持这类问题会直接指出来。nodata 打印出来是 None 的情况也很常见意思是产品没有声明填充值。这时候不能当成没有无效像元而是要按产品文档手工指定或者退一步用 DN 为 0 作为无效条件。这一步判断错了边缘一圈填充值会被当成真实观测参与定标。3.2 分块流式处理窗口、行数与内存上限分块读的核心是ReadAsArray(xoff, yoff, xsize, ysize)这四个参数。按行切是最省事的做法每次读一个固定高度的条带处理完直接写出去内存占用只跟块大小有关跟影像多大无关。def iter_blocks(band, block_h512): 按行分块yield (yoff, rows, array) ysize band.YSize for y in range(0, ysize, block_h): rows min(block_h, ysize - y) # 最后一块可能不满 arr band.ReadAsArray(0, y, band.XSize, rows) yield y, rows, arrrows这一行必须用min截一下。我见过有人直接用block_h最后一块多读了越界的行GDAL 不报错返回的数组却比预期短写回去的时候位置整体错位成品图上出现一条横向的错位带非常难查。块高度取多少合适512 行对于 8 位或 16 位单波段大概几 MB 到几十 MB多波段就相应往下调。这个值不用太精确关键是别一次读全图。还有一种情况是影像带金字塔或者被分块存储TILEDYES这时按存储块对齐读会快一些。但对定标这种逐像元的线性变换来说IO 不是瓶颈按行读就够用没必要为了对齐去增加代码复杂度。3.3 写出 Float32 GeoTIFF压缩、nodata 与缩放输出格式的选择有两个取向。一是直接写 Float32反射率保留 0~1 的小数下游做反演最方便二是缩放成 UInt16比如乘 10000 存整数文件小一半和很多成型的反射率产品格式一致但下游必须记得除回来。drv gdal.GetDriverByName(GTiff) dst drv.Create( out_reflectance.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32, options[TILEDYES, COMPRESSLZW, PREDICTOR3, BIGTIFFIF_SAFER], ) dst.SetGeoTransform(ds.GetGeoTransform()) dst.SetProjection(ds.GetProjection()) ob dst.GetRasterBand(1) ob.SetNoDataValue(-9999.0)PREDICTOR3是浮点数据的专用预测器配合 LZW 压缩能把浮点影像压到原始大小的三到四成比不设预测器明显。整数数据则用PREDICTOR2。BIGTIFFIF_SAFER是给大影像留的后路超过 4GB 自动切 BigTIFF省得跑到一半失败。nodata 的设置要和实际写入的填充值对上。写 Float32 时用 -9999写 UInt16 缩放时用 0两者不能混。我见过有人 nodata 设成 0但填充值写的是 -9999结果下游软件读完把整幅图当成有效区统计出来的均值完全不对。4. 一份能直接跑的定标脚本参数、计算与输出把前面的碎片拼起来就是一份完整脚本。我习惯把系数放 JSON、几何参数放命令行脚本本身不出现任何硬编码的传感器数值这样换一景、换一个波段只改参数不改代码。下面的系数值全部是模拟项目里随手编的占位真实处理时换成产品文档里的表。4.1 系数配置和命令行参数配置按波段组织每段三个字段。把单位也写进配置里虽然脚本不读但下次翻回来能立刻确认量纲。{ sensor: sim-optical-a, unit: W·m-2·sr-1·um-1, doy: 168, sun_zenith_deg: 27.4, bands: { 1: {gain: 0.0123, offset: -0.45, esun: 1958.0}, 2: {gain: 0.0131, offset: -0.52, esun: 1850.0} } }doy和sun_zenith_deg放在顶层是因为它们和波段无关。如果产品给了逐像元角度层这两个字段就退化成默认值脚本里加一个可选的角度文件参数覆盖掉。命令行参数用 argparse必填项只留源文件、目标文件、系数文件和年积日其余给默认值。这样批量跑的时候命令行短手工调试的时候也能临时改块大小和缩放系数。4.2 核心计算循环与掩膜处理#!/usr/bin/env python3 # calib.py 单波段绝对辐射定标DN - 辐亮度 - 大气顶反射率 import argparse import json import numpy as np from osgeo import gdal gdal.UseExceptions() def read_coef(cfg_path, band_id): with open(cfg_path, r, encodingutf-8) as f: cfg json.load(f) b cfg[bands][str(band_id)] return float(b[gain]), float(b[offset]), float(b[esun]) def sun_earth_distance(doy): g 2.0 * np.pi * (doy - 1) / 365.0 return 1.00014 - 0.01671 * np.cos(g) - 0.00014 * np.cos(2 * g) def main(): p argparse.ArgumentParser() p.add_argument(--src, requiredTrue) p.add_argument(--dst, requiredTrue) p.add_argument(--band, typeint, default1) p.add_argument(--coef, requiredTrue) p.add_argument(--doy, typeint, requiredTrue) p.add_argument(--sza, typefloat, requiredTrue, help太阳天顶角90-高度角) p.add_argument(--scale, typeint, default0, help0 时缩放并写 UInt16) p.add_argument(--block, typeint, default512) a p.parse_args() gain, offset, esun read_coef(a.coef, a.band) # 把反射率公式里的常数项合并成一个系数循环里只做一次乘加 d2 sun_earth_distance(a.doy) ** 2 cos_s max(np.cos(np.deg2rad(a.sza)), 0.05) k np.pi * d2 / (esun * cos_s) src gdal.Open(a.src, gdal.GA_ReadOnly) band src.GetRasterBand(a.band) nodata band.GetNoDataValue() dtype gdal.GDT_UInt16 if a.scale 0 else gdal.GDT_Float32 drv gdal.GetDriverByName(GTiff) dst drv.Create( a.dst, src.RasterXSize, src.RasterYSize, 1, dtype, options[TILEDYES, COMPRESSLZW, PREDICTOR2 if a.scale 0 else PREDICTOR3, BIGTIFFIF_SAFER], ) dst.SetGeoTransform(src.GetGeoTransform()) dst.SetProjection(src.GetProjection()) ob dst.GetRasterBand(1) valid 0 for y in range(0, src.RasterYSize, a.block): rows min(a.block, src.RasterYSize - y) dn band.ReadAsArray(0, y, src.RasterXSize, rows).astype(np.float32) mask np.ones(dn.shape, dtypebool) if nodata is not None: mask dn ! nodata mask dn 0 # 0 一般是填充值按文档调整 rho (gain * dn offset) * k # 辐亮度那一步已经合进 k 里 out np.where(mask, rho, np.nan) valid int(mask.sum()) if a.scale 0: arr np.where(np.isnan(out), 0, np.clip(out * a.scale, 0, 65535)).astype(np.uint16) else: arr np.where(np.isnan(out), -9999.0, out).astype(np.float32) ob.WriteArray(arr, 0, y) ob.SetNoDataValue(0 if a.scale 0 else -9999.0) ob.FlushCache() dst None print(fdone, valid px {valid}) if __name__ __main__: main()几个参数值得单独说。--sza是太阳天顶角不是高度角写命令行之前先算一下。--scale设成 10000 就输出 UInt16 的反射率下游读的时候除以 10000设成 0 就输出 Float32不做缩放。--block是分块行数普通机器 512 够用内存紧张就调到 256。k这个合并系数是把 π、d²、ESUN、cos θs 提前算好的结果循环里每个像元只做一次乘加。这种合并在数学上完全等价但要注意它把 cos θs 固定成了一个常数如果要做逐像元角度层就得把 k 拆开成数组参与计算不能在循环外算死。valid计数器看着不起眼实际很有用。跑完打出来和影像总像元数一比就知道掩掉了多少正常应该在 90% 以上低得离谱就要回去查填充值设错了没有。4.3 跑完先看三个统计量脚本跑完不要直接拿去做下一步先做一次快速统计。反射率的物理范围是 0 到 1 之间水体、阴影、云这些会偏低裸地植被偏高但正常景的 99 分位数一般不会超过 0.8。如果 p99 跑到 3 或者 5说明要么系数方向错了要么 cos 截断没生效要么掩膜漏了大量填充值。gdalinfo -stats out_reflectance.tif | tail -20-stats会算最小、最大、均值、标准差几秒钟出结果比写脚本快。重点看最小值是不是 -9999说明 nodata 生效最大值有没有超过 1 太多均值在不在 0.1 到 0.4 这个典型区间。这三项一起看八成的问题能定位到。5. 避坑与排查五条踩过的坑定标脚本的 bug 有个共同特征就是不会报错。数组照常算完影像照常写出只是数值不对。所以排查要靠对比和统计不能靠看代码。下面五条是我自己或者身边开发者真实翻过车的记录按现象、原因、解决写。5.1 反射率整景出现一块特别大的数值现象输出的反射率最大值到几万最小值还是正常的 0.05 左右直方图右侧拖一条很长的尾巴。原因分母里的 cos θs 用了太阳高度角而不是天顶角高度角几十度的时候 cos 值还大于 0.5但有一批像元被算成了负数天顶角cos 变成负的或者接近零除法直接把值放大。解决确认元数据字段语义sun_zenith 90 - sun_elevation然后在脚本里给 cos 加下限截断同时把天顶角大于 80 度的像元整体掩掉。改完重跑p99 立刻回到 0.6 以下。5.2 成品图上出现横向错位带现象大图缩略看正常放大到某几行会有一条内容错开的横带位置大概在每个分块的接缝处。原因分块写回时WriteArray的 yoff 用了块序号乘块高但循环里用的是实际偏移 y两者在最后一块对不上还有一种情况是rows没截断多读的行被写到了下一块的位置。解决统一用循环变量 y 作为写回的偏移rows一定用min算。改完跑一遍把源和结果各自切同一行做差分差值应该是一个常数乘子不是错位。5.3 输出全是 nodata 或者有效像元极少现象valid px打出来只有总像元的百分之几影像大面积空白。原因掩膜条件写成了dn ! 0但这景产品的填充值其实是 65535 或者 -9999也可能产品本身声明了 nodata脚本又额外加了dn 0两道条件叠加把正常像元筛掉了。解决先不带掩膜跑一小块把 DN 的唯一值统计出来看填充值到底是什么。掩膜条件只保留一条或者用np.isin(dn, [fill_list])的方式显式排除。5.4 和成熟软件的结果差一个数量级现象同一景影像脚本算出来的辐亮度和商业软件对比数值差了正好 10 倍或者 100 倍。原因量纲不一致。辐亮度的单位有的是 W·m⁻²·sr⁻¹·μm⁻¹有的是 μW·cm⁻²·sr⁻¹·nm⁻¹后者换算到前者要乘一个 10 左右的因子ESUN 的单位也可能跟着不匹配两边不同源就会整体偏移。解决把 gain、offset、ESUN 三个值的单位写进配置计算前统一换算到一套单位。我一般统一到 W·m⁻²·sr⁻¹·μm⁻¹ 和对应的辐照度换算完在日志里打一行确认。5.5 跑到一半内存暴涨然后失败现象前几块跑得很快跑到中间内存持续上涨最后进程被系统结束。原因循环里创建了中间数组却没释放比如每一步都np.stack累积结果或者ReadAsArray没带窗口一次读了整景分块循环形同虚设。解决检查读数组那行有没有带四个参数中间结果用np.where就地生成不要累积循环结束把块变量置空。改完内存曲线应该是一条平的直线。6. 精度复核与批量落地抽样对比和降采样预跑脚本能跑通只是第一步数值对不对得单独验。我一般做两层复核。第一层是抽样点对比在景内挑几类均匀地物比如深水体、裸沙、密植被各取一个 5×5 的小窗口把原始 DN 均值、算出来的辐亮度均值、反射率均值列出来看同类地物在不同景之间是否一致。第二层是和产品自带的大气顶反射率层有些 L1 产品会附一层做逐像元差分算 RMSE 和偏差判定标准可以照下面这张表。复核项检查方式合理区间不通过时先查反射率范围p1 与 p990.01 ~ 0.8系数方向、cos 截断有效像元占比valid / 总像元大于 90%填充值、掩膜条件与参考层偏差均值差绝对值小于 0.02单位换算、ESUN 来源均匀地物一致性同类窗口均值相对差小于 5%太阳天顶角、日地距离直方图形态目视单峰、无长尾分母趋零、错位写入还有一个我自己养成的习惯全分辨率跑之前先用gdal_translate -outsize 6.25% 6.25%生成一张 1/16 的缩略图拿它把参数跑通。缩略图几秒钟就能出结果统计量能立刻看出方向对不对等确认无误再上全分辨率省下的时间远超过这一步的开销。批量处理多景的时候我会写一个 bash 循环把每景的年积日和天顶角从文件名或者 XML 里带进去跑完统一出一份统计汇总哪一景的 p99 异常就单独回去看。从那以后我每次上批量之前都强制先跑一遍缩略图加统计宁可多花五分钟也不想第二天早上看到几百景全废。希望帮到你。本文还有配套的精品资源点击获取