简介面向遥感图像处理初学者与Python开发者的一套轻量实践资料围绕绝对辐射定标这一关键预处理环节演示如何借助GDAL将影像原始数字值转换为具有物理意义的辐射量并衔接后续反射率、亮温及大气校正等分析场景。压缩包共9个文件、约16KB以5个XML元数据或工程配置、1个Python示例脚本、1个JSON定标参数文件为主体另含gitignore与iml工程文件结构紧凑便于在本地IDE中直接查看和改造。已有903人学习下载。资料提供了可复用的脚本骨架与参数组织方式覆盖影像波段读取、定标系数加载、逐波段校正计算及结果写出等环节并保留工程配置与元数据示例适合对照自身传感器参数做二次开发也可作为遥感辐射定标流程的学习参考与代码排错起点。1. 从DN到辐亮度一景影像进定量反演前必须跨过的门槛上个月有同学拿来两景同一片农田的影像问为什么隔了半个月近红外波段的数值差了一大截。我翻到元数据一看两景的RADIANCE_MULT_BAND_4差了将近一倍——DN值本身没有物理意义直接跨影像比较等于拿两台没校准过的秤读数。遥感影像要做定量分析绝对辐射定标就是那道绕不开的门槛用 Python 加 GDAL 把传感器输出的整数 DN 还原成辐亮度或表观反射率让不同时相、不同传感器、不同增益下的数据能放在同一个尺度上说话。GDAL 负责分块读写和地理参考numpy 负责逐像元运算一套脚本同时能吃 Landsat 的 MTL、Sentinel-2 的 XML 和国产高分系列的定标文件。这篇写给已经能打开影像、但定标这一步还靠手工点鼠标的人。2. 定标公式与系数来源把MTL、XML和定标文件对齐动手写代码之前得先弄清楚自己在算哪个物理量。遥感里常说的定标其实是一条链DN 到辐亮度辐亮度到表观反射率表观反射率再经过大气校正才是地表反射率。绝对辐射定标只负责前两环第三环是大校正的事混在一起做结果对不上账。这一章把公式和系数来源拆开后面写代码时才知道每个数字该从哪个字段取。2.1 物理量链条DN、辐亮度、表观反射率各是什么DN 是传感器记录的数字量化值无量纲本质上是一个 ADC 读数。同一颗卫星在不同时期可能调整增益同一景影像不同波段的量化范围也不一样所以 DN 之间没有可比性。辐亮度的单位是 W·m⁻²·sr⁻¹·μm⁻¹描述的是传感器入瞳处单位面积、单位立体角、单位波长宽度内接收到的辐射能量这是一个有物理量纲的值可以跨传感器比较。表观反射率再进一步把辐亮度除以入射太阳辐照度并做太阳天顶角余弦修正归一到 0 到 1 之间好处是消除了太阳高度角和日地距离的季节变化。常用的两个公式分别是Lλ ML × Qcal AL和ρλ (Mρ × Qcal Aρ) / sin(θSE)。第一式中 ML、AL 是辐亮度乘系数与加系数从元数据里直接读第二式中 Mρ、Aρ 是反射率乘系数与加系数θSE 是太阳高度角。很多传感器只提供辐亮度系数反射率系数要自己用 ESUN 反推公式是ρ π × L × d² / (ESUN × cosθs)d 是日地距离修正因子。我一般优先用元数据里现成的反射率系数除非确认它缺失因为自己推的 ESUN 版本不同结果会有千分之几的偏差做时序分析时这千分之几会被放大。还有一个容易搞混的地方表观反射率不等于地表反射率。前者包含大气散射和吸收的影响天空晴朗时蓝波段的表观反射率会比地表高出好几倍这是大气程辐射造成的不是定标算错了。拿到定标结果直接做 NDVI 或者水体指数问题不大因为比值型指数对乘性因子不敏感但要做地表温度反演、生物量回归这类吃绝对值的模型就必须再接一步大气校正。2.2 Landsat 8/9 的MTL字段与国产卫星XML的对应关系Landsat 系列的好处是元数据齐全一景 Collection 2 Level-1 产品的_MTL.txt里每个波段都成对给出了辐亮度和反射率的乘加系数外加太阳高度角、方位角、日地距离。字段名规律性很强RADIANCE_MULT_BAND_1、RADIANCE_ADD_BAND_1、REFLECTANCE_MULT_BAND_1、REFLECTANCE_ADD_BAND_1把波段号替换掉就行不需要为每个波段写一遍解析逻辑。Sentinel-2 L1C 走的是另一套产品里存的是放大了一万倍的反射率整数定标就是把 DN 除以 10000再除以 sin(太阳高度角)。它的元数据在MTD_MSIL1C.xml里用 XPath 取QUANTIFICATION_VALUE和每个波段的SOLAR_IRRADIANCE不同处理基线版本的字段层级会变硬编码路径容易翻车建议先打印一份 XML 结构再写解析。国产高分系列GF-1、GF-2 这类通常不随影像附带 XML定标系数单独发布在定标参数文件里一份文件对应一个传感器一个年份。这类系数一般只有一组增益值公式是L DN / Gain没有加系数。表里整理了三种来源的差异写代码时按这个分支处理数据源系数载体辐亮度公式反射率是否现成Landsat 8/9 L1MTL.txtL ML×DN AL现成含太阳角Sentinel-2 L1CMTD_MSIL1C.xmlρ DN / 10000 / sinθ现成国产高分系列独立定标文件L DN / Gain需用 ESUN 自算2.3 乘法还是除法定标系数的方向判断与样本核算这是国内数据上最容易翻车的一步。Landsat 的系数是乘系数量级在 10⁻² 左右DN 乘上去得到几十到几百的辐亮度。国产高分的定标文件里写的是绝对定标系数本质是增益单位是 DN/(W·m⁻²·sr⁻¹·μm⁻¹)数值动辄几百上千公式是 DN 除以它。两种写法长得像方向反了结果能差三到四个数量级。判断方向不能只看文件名要看三样东西系数文件里的单位说明、系数数量级、以及拿一景影像试算后的量级。我常用的做法是先取一小块样本两种方向各算一遍和人眼能判断的物理常识对一下。植被在近红外波段的辐亮度一般在几十的量级如果算出来是 0.00x 或者上万方向基本就错了。import numpy as np def check_direction(dn_sample: np.ndarray, coef: float) - None: 用一小块样本核算定标方向避免整景跑完才发现反了 dn dn_sample.astype(np.float64) dn dn[dn 0] # 先剔掉背景 0 值 if dn.size 0: print(样本全为 0换一块区域) return as_mult dn * coef as_div dn / coef print(fDN 均值: {dn.mean():.2f}) print(f按乘法: 均值 {as_mult.mean():.4f}, 范围 [{as_mult.min():.4f}, {as_mult.max():.4f}]) print(f按除法: 均值 {as_div.mean():.4f}, 范围 [{as_div.min():.4f}, {as_div.max():.4f}]) # 近红外辐亮度合理量级在 10^1 附近落在 10^-3 或 10^5 的那个方向就是错的这段代码不做任何写盘只用来在正式跑之前做一次量级体检。coef传定标文件里的原始值不要预先取倒数两种方向让它自己算。判据很简单近红外波段植被区辐亮度落在 10 到 100 之间算正常蓝波段因为大气散射会高一些但也在同一量级。如果两个方向的输出离这个区间都很远那多半是系数文件本身对应错了传感器或者年份得回去查定标文件的适用说明。3. GDAL分块流水线写一个能吃下GB级影像的定标脚本把公式落成代码难点不在数学在工程。一景 Landsat 8 全波段产品解压后接近 1GB国产高分成像幅宽更大整幅读进内存再运算普通办公机直接爆。GDAL 的分块读写接口就是为这个场景设计的配合 numpy 的布尔掩膜可以在常数内存下把整景处理完。这一章的脚本按能直接抄去改的标准写参数都留了口子。3.1 环境准备与GDAL的三种安装方式GDAL 是这套流程里唯一有安装门槛的依赖。pip 直接装gdal经常因为缺少系统级 libgdal 编译失败尤其是 Windows 上。我给的建议按成功率排序conda 最稳conda install -c conda-forge gdal numpy其次是 pip 装预编译轮子pip install GDAL$(gdal-config --version)这要求本机已经装了 GDAL 命令行工具最不推荐源码编译除非你有特定版本需求。# 推荐路径conda 环境里一次装齐 conda create -n calib python3.10 -y conda activate calib conda install -c conda-forge gdal numpy -y # 验证安装输出应大于等于 3.4 python -c from osgeo import gdal; print(gdal.__version__)版本上没有硬性要求但建议 GDAL 3.4 以上因为它对 COG 和 BIGTIFF 的支持更完整处理超过 4GB 的输出文件时不用额外配置。numpy 用 1.24 以上即可。装完之后先跑一遍上面那行验证ImportError多半是 conda 环境和系统 Python 混用了which python确认一下路径。3.2 元数据解析别用硬编码波段号元数据解析的原则是波段号从文件名或循环变量来不要写死在函数里。一景 Landsat 有 11 个波段每个波段单独一个 TIF如果函数里写死BAND_4换一景就得改代码。我一般把解析和计算拆成两个函数解析负责把 MTL 变成一个扁平字典计算函数只接收mult、add两个浮点数这样国产卫星的定标文件也能塞进同一个计算函数。from pathlib import Path def parse_mtl(mtl_path: str) - dict: 把 Landsat MTL.txt 解析成扁平字典GROUP 行和注释行自动跳过 meta {} text Path(mtl_path).read_text(encodingutf-8, errorsignore) for raw in text.splitlines(): line raw.strip() if not line or not in line: continue key, val line.split(, 1) meta[key.strip()] val.strip().strip() return meta def landsat_coef(meta: dict, band: int, kind: str radiance): 按波段号取系数kind 决定取辐亮度还是反射率那一组 if kind radiance: mult float(meta[fRADIANCE_MULT_BAND_{band}]) add float(meta[fRADIANCE_ADD_BAND_{band}]) return mult, add, None mult float(meta[fREFLECTANCE_MULT_BAND_{band}]) add float(meta[fREFLECTANCE_ADD_BAND_{band}]) sun_elev float(meta[SUN_ELEVATION]) return mult, add, sun_elev if __name__ __main__: m parse_mtl(LC08_L1TP_XXX_MTL.txt) print(landsat_coef(m, band4, kindradiance)) print(landsat_coef(m, band4, kindreflectance))split(, 1)里的1不能省MTL 里个别字段的值本身带等号不限制分割次数会拆出三段后面float()直接报错。strip()处理带引号的字符串字段对数值字段无副作用。反射率分支一定要把SUN_ELEVATION一起返回后面做太阳角修正时用得上漏掉这一步反射率会整体偏大太阳高度角 45 度时偏差接近 41%。3.3 分块读写窗口、数据类型与压缩参数主函数的骨架是两层循环加ReadAsArray(x, y, cols, rows)。三个细节决定脚本能不能跑通大文件输出数据类型必须是GDT_Float32块大小取 512 或 1024压缩用DEFLATE PREDICTOR3。PREDICTOR 是浮点预测器对定标这类平滑变化的浮点数据压缩率能提升一倍左右而且是无损的不改变数值。import numpy as np from osgeo import gdal gdal.UseExceptions() def calibrate(src_path, dst_path, mult, add, sun_elevNone, scaleNone, nodata_in0, nodata_out-9999.0, block512, clip_rangeNone): src_path : 输入单波段 DN 影像 mult/add : 定标乘系数与加系数国产数据传 mult1.0/gain, add0 sun_elev : 太阳高度角度传 None 表示只做辐亮度 scale : 额外缩放因子如 Sentinel-2 传 10000 clip_range : 反射率裁剪区间如 (0.0, 1.0) src gdal.Open(src_path, gdal.GA_ReadOnly) if src is None: raise RuntimeError(f打不开输入影像: {src_path}) band src.GetRasterBand(1) w, h src.RasterXSize, src.RasterYSize dst gdal.GetDriverByName(GTiff).Create( dst_path, w, h, 1, gdal.GDT_Float32, options[TILEDYES, BLOCKXSIZE512, BLOCKYSIZE512, COMPRESSDEFLATE, PREDICTOR3, BIGTIFFIF_SAFER] ) dst.SetGeoTransform(src.GetGeoTransform()) dst.SetProjection(src.GetProjection()) out_band dst.GetRasterBand(1) out_band.SetNoDataValue(nodata_out) sun_sin None if sun_elev is not None: sun_sin np.sin(np.deg2rad(sun_elev)) if sun_sin 0: raise ValueError(f太阳高度角异常: {sun_elev}) for y in range(0, h, block): rows min(block, h - y) for x in range(0, w, block): cols min(block, w - x) dn band.ReadAsArray(x, y, cols, rows) if dn is None: continue dn dn.astype(np.float32) valid dn ! nodata_in out np.full(dn.shape, nodata_out, dtypenp.float32) val dn[valid] * mult add if scale is not None: val val / scale if sun_sin is not None: val val / sun_sin if clip_range is not None: np.clip(val, clip_range[0], clip_range[1], outval) out[valid] val out_band.WriteArray(out, x, y) out_band.FlushCache() dst None src None return dst_pathmin(block, h - y)是防止最后一块越界的唯一正确写法直接用block会让ReadAsArray返回一个比预期小的数组WriteArray时形状不匹配落盘结果右下角出现黑边或者错位。valid dn ! nodata_in先把无效像元摘出来只在有效位置做运算无效位置统一填nodata_out这样 NoData 不会参与乘加后变成一个看似合理的假值。clip_range只对反射率有意义辐亮度不要裁因为地物辐亮度的物理上限不好定裁了反而引入人为误差。3.4 反射率分支太阳高度角与0-1裁剪反射率分支在辐亮度基础上多一步除法参数从元数据里来代码里不做任何硬编码。Landsat 的SUN_ELEVATION单位是度np.deg2rad必须显式调用numpy 的三角函数只认弧度。裁剪区间传(0.0, 1.0)物理上表观反射率不会超过 1超过的部分要么是云要么是系数配错了裁掉能防止后续指数计算出无穷大。if __name__ __main__: meta parse_mtl(LC08_L1TP_XXX_MTL.txt) for b in range(1, 8): mult, add, sun landsat_coef(meta, bandb, kindreflectance) calibrate( src_pathfLC08_L1TP_XXX_B{b}.TIF, dst_pathfLC08_L1TP_XXX_B{b}_SRF.tif, multmult, addadd, sun_elevsun, nodata_in0, nodata_out-9999.0, clip_range(0.0, 1.0), )这段循环把 1 到 7 波段一次性处理完输出文件名带_SRF后缀和原始 DN 文件区分开。波段 6 是热红外定标公式不同不能套这套系数代码里单独跳过。波段 8 是全色分辨率和其他波段不一致单独处理更省事。国产数据的调用方式是把mult换成1.0 / gainsun_elev传None因为多数国产定标文件给的系数直接对应辐亮度不带太阳角修正。4. 反射率与质量掩膜让定标结果能直接进反演定标跑完不等于能用。云覆盖、云影、水体边缘的混合像元这些在原始 DN 上看着正常定标之后会变成反射率里的极值直接喂给反演模型会把回归系数带偏。Landsat 和 Sentinel-2 都带了质量波段用位运算解析出来做成掩膜是让结果能进下游流程的最后一道工序。批量处理也得在这一步解决不然一景一景点几十景数据能耗掉一整天。4.1 QA_PIXEL位运算云、云影、水体怎么掩Landsat Collection 2 的QA_PIXEL是一个 16 位整数每一位代表一类质量标记。常用的几位是 bit 1膨胀云、bit 3云、bit 4云影、bit 5雪、bit 6晴空、bit 7水体。判云不能只判 bit 3膨胀云标记覆盖范围更保守做定量分析时连膨胀云一起掩掉更安全。def build_qa_mask(qa: np.ndarray, cloudTrue, shadowTrue, snowFalse) - np.ndarray: 返回 True 表示该像元需要被掩掉 bits 0 if cloud: bits | (1 1) | (1 3) # 膨胀云 云 if shadow: bits | (1 4) if snow: bits | (1 5) if bits 0: return np.zeros(qa.shape, dtypebool) return (qa bits) ! 0位运算的关键是先把要掩的位或起来再去和 QA 值做与运算比逐位判断快得多。(qa bits) ! 0返回布尔数组直接拿来索引反射率数组赋nodata_out就行。雪要不要掩看研究区做积雪反演的当然不能掩做植被时序的必须掩雪在可见光的反射率接近 0.9和云的干扰差不多。掩膜应用的时候有个顺序问题先掩膜再算指数还是先算指数再掩膜。我的习惯是先掩膜因为指数计算里一旦掺进云的反射率比值会出现大于 1 或者负值虽然最终也会被掩掉但中间过程出现的 inf 和 nan 会污染统计量。如果非要在后处理阶段掩记得用np.errstate把除零警告压掉。4.2 批量调度进程池、命名规范与断点续跑几十景数据串行跑一景三分钟加起来就是两三个小时。GDAL 的读写是 IO 密集型多进程能明显提速一般开到 CPU 核数的一半开满反而因为磁盘争抢变慢。进程池里每个任务传独立的文件路径不要传已打开的数据集对象GDAL 的数据集句柄不能跨进程共享。from concurrent.futures import ProcessPoolExecutor, as_completed def run_batch(tasks, workers4): ok, fail [], [] with ProcessPoolExecutor(max_workersworkers) as ex: futures {ex.submit(calibrate, **t): t[src_path] for t in tasks} for fut in as_completed(futures): name futures[fut] try: fut.result() ok.append(name) print(fOK {name}) except Exception as e: fail.append((name, str(e))) print(fFAIL {name}: {e}) print(f完成 {len(ok)}失败 {len(fail)}) return ok, failtasks是一个字典列表每个字典的键必须和calibrate的形参名对上这样**t展开才不会报错。断点续跑的逻辑放在任务生成阶段先扫一遍输出目录把已经存在的目标文件对应的任务剔除再提交给进程池。失败的任务单独记录下来人工查完原因重新跑不要自动重试失败通常是数据本身有问题重试也是白跑。命名规范建议固定成{传感器}_{日期}_{波段}_{处理级别}.tif处理级别用RAD表示辐亮度、TOA表示表观反射率。这样下游做镶嵌或者时序分析时用 glob 模式一批就能筛出来不用打开每个文件读元数据。我有一次偷懒没带处理级别后缀结果把辐亮度和反射率混在一个目录里做了均值统计数字出来怎么都对不上回头排查花了两个小时。4.3 结果自检统计量、直方图与异常值分布每景跑完都该看一眼统计量不用精细但必须看。方法是读输出文件排除 NoData 后算均值、分位数和超出物理范围的比例。反射率正常应该在 0 到 1 之间如果超过 1 的像元占比超过 5%要么是云没掩干净要么是系数配错。def quick_stats(path: str, nodata-9999.0) - dict: ds gdal.Open(path) arr ds.GetRasterBand(1).ReadAsArray() arr arr[arr ! nodata] if arr.size 0: return {n: 0} return { n: int(arr.size), mean: float(arr.mean()), p01: float(np.percentile(arr, 1)), p50: float(np.percentile(arr, 50)), p99: float(np.percentile(arr, 99)), max: float(arr.max()), over1: float((arr 1.0).mean()) if arr.max() 1 else 0.0, }p01和p99比min和max更可靠因为单像元的坏点会把min/max拉爆。p99接近 1、p50在 0.1 到 0.3 之间是植被区影像的常见形态如果p50接近 0检查一下是不是把 DN 值直接当反射率输出了。over1只在最大值超过 1 时才计算省一次全数组比较。5. 避坑与排查五条在定标现场翻车的记录定标脚本写出来第一次能跑通不代表结果是对的。下面五条是我和身边人实际遇到过的每一条都按现象、原因、解决三段写对照着排查能省掉大量试错时间。5.1 输出写成UInt16定标结果只剩0和1现象输出影像在 GIS 里打开一片黑统计出来只有 0、1 两个值和原始 DN 影像的纹理完全对不上。原因Create时数据类型写成了gdal.GDT_UInt16。反射率在 0 到 1 之间乘加之后转成整型小数部分全部截断0.3 变成 00.95 也变成 0只有极少数超过 1 的被截成 1。解决输出类型固定为gdal.GDT_Float32创建之后再写SetNoDataValue。如果下游软件只认整型把反射率乘以 10000 转成整型存储同时在元数据里标注缩放因子别在原地截断。5.2 国产系数当乘法用量级差出三个数量级现象辐亮度结果都是 0.001 量级或者反过来几万和同区域的 Landsat 产品对比差了上千倍。原因国产高分的定标文件给的是增益公式为L DN / Gain误当成了L DN × Gain。系数本身在几百的量级乘和除差了系数平方倍。解决先用 2.3 节的样本核算脚本跑一遍看两个方向的量级哪个落在物理区间。代码里把mult1.0/gain显式写出来并在注释里标注系数来源文件和公式别让下一个人再猜。5.3 NoData0参与运算边缘出现鬼影现象影像有效区外围出现一圈均匀的低值数值刚好等于加系数 AL或者反射率等于太阳角修正后的 Aρ/sinθ。原因没有先掩膜就做乘加。DN 为 0 的无效像元代入DN×ML AL乘项归零剩下的就是 AL看起来像一个合理的辐亮度实际上没有任何物理意义。解决计算前先做valid dn ! nodata_in只在有效位置运算无效位置统一填-9999。填的值不能是 0因为 0 在反射率里是合法值深水、阴影都可能接近 0用 0 当 NoData 会让下游分不清真假。5.4 太阳高度角没转弧度反射率整体飘高现象反射率均值比预期高出三到五成植被近红外到了 0.6 以上水体蓝波段超过 0.3。原因np.sin(sun_elev)直接把角度值传进去了。sin(45) 的弧度计算结果是 0.85而 sin(45°) 是 0.707差了两成太阳高度角越小偏差越大冬季高纬度影像能差出一倍。解决np.sin(np.deg2rad(sun_elev))并且在函数入口加一次if sun_sin 0: raise。太阳高度角为负说明是极夜这种影像本来就不该做光学定标。5.5 整幅ReadAsArray4GB影像直接把内存吃满现象小影像跑得好好的换成一景宽幅国产影像进程直接被系统杀掉日志里没有任何 Python 异常。原因band.ReadAsArray()不带参数读整幅一个 30000×30000 的 UInt16 波段就是 1.8GB转成 Float32 翻倍到 3.6GB中间再产生几个临时数组内存峰值轻松超过 10GB。解决用 3.3 节的两层循环分块读写块大小 512 或 1024内存峰值控制在几十 MB。同时给Create加上BIGTIFFIF_SAFER输出超过 4GB 时自动切 BigTIFF不然写到一半报错。6. 交叉验证拿官方产品和直方图给自己兜底定标结果对不对光看代码逻辑看不出来得有一个外部参照。最直接的办法是拿同一景的官方 L2 产品做交叉验证但要注意 L2 是地表反射率已经做过了大气校正和 TOA 反射率直接比会有系统性偏差正确做法是比较两者的空间格局和地物间的相对关系而不是逐像元差值。更实用的参照是物理区间表把结果按波段统计出来对照经验范围判断。波段类型辐亮度经验区间TOA 反射率经验区间越界时的常见原因蓝、绿水体30900.030.12气溶胶厚或系数方向反了红植被15500.030.10太阳角未修正近红外植被20800.200.45云未掩干净拉高均值短波红外5300.100.30增益配错或用了错误年份系数区间只是量级参考不同大气条件、不同物候下浮动很大别当成硬阈值卡。真正好用的自检是同一传感器不同时相的直方图对比找两景同一区域、不同日期的影像分别定标后画反射率直方图两条曲线的峰位应该接近形状差异来自物候和大气。如果峰位差了 0.1 以上八成是其中一景的系数用错了。还有一个低成本的办法是选稳定目标做时序核查。深水水体在近红外波段的反射率长期接近 0沙漠在可见光波段的反射率长期稳定在 0.3 到 0.4拿这些已知目标点几个样本半年内的定标结果波动超过 10% 就说明流程里有环节不稳定。我自己是在批量脚本跑完之后加一步随机抽三景各取一万个有效像元算中位数和上一批次的结果对比偏差超过阈值就停下来查不往下游送。从那以后我每次做新传感器的定标都强制先跑一遍样本核算和统计自检两关都过了才允许进批量。定标这一步出错下游所有的反演、分类、时序分析全是在错误的基座上盖楼返工成本比多花半小时自检高得多。希望帮到你。本文还有配套的精品资源点击获取