简介遥感图像的绝对辐射定标是将原始DN值转换为地表反射率或辐射亮度等物理量的关键预处理步骤。这一Python工程面向遥感初学者、地学研究人员及GIS开发人员演示如何利用GDAL库打开多波段影像、从JSON文件中读取增益与暗电流等定标参数并完成逐像元辐射校正同时给出大气校正与元数据整理的思路可结合实际任务灵活调整。压缩包共9个文件核心包括CrossRadiationCalibration.py可执行脚本和RadiometricCorrectionParameter.json参数文件另含工程配置文件与XML说明便于理解整个工程结构整体仅16KB轻量且易于阅读。目前已有896人学习下载适合希望快速掌握遥感辐射定标编程实现、理解传感器定标系数应用逻辑的入门者。借助该资源使用者能梳理从数据准备、参数解析、辐射值计算到结果保存的完整流程并基于示例代码扩展自己的处理脚本有效降低从理论到实践的试错成本。1. 概念拆解为什么绝对辐射定标是遥感定量化的第一道门槛拿到一景L1级影像很多人第一反应是直接拿DN值去算NDVI。我之前也这么干过直到一次做长时间序列的植被指数变化分析发现同一条样带的NDVI曲线在季度交界处出现莫名其妙的跳变。排查到最后问题出在影像之间的辐射基准不一致——说白了DN值只是传感器把接收到的能量量化成的一个整数它跟地物真实反射太阳辐射的能力之间还隔着一层“翻译”过程。这层翻译就叫绝对辐射定标。绝对辐射定标解决的问题是把影像记录的DN值转换成具有物理意义的辐亮度或者反射率。DN值的大小不仅受地物本身反射率影响还跟太阳高度角、传感器增益设置、大气路径辐射等因素有关。如果不做定标两景不同时间、不同角度、甚至不同传感器获取的影像根本没法在数值层面直接对比。这一点在定量遥感里是底线在目视解译里无所谓但在做反演、做光谱分析、做机器学习样本标注时踩不踩这一步结果能差出一个量级。用生活化的方式理解DN值像是一张照片上像素的明暗程度同样的白墙中午拍和傍晚拍亮度完全不一样辐亮度则是把相机换成照度计测的是这面墙实际反射了多少光能量。绝对辐射定标就是帮你把“照片的明暗”换算成“物理上的亮度”后面再做大气校正、反射率反演、地表温度反演才有共同的语言基础。本文面对的读者我默认是遥感或地学相关专业的学生、刚入行做遥感数据处理的技术人员以及需要在Python里批量处理遥感影像的开发者。我会从原理讲到代码实现再到坑点排查全程用Landsat系列和国产高分系列的数据做例子尽量把每一步讲透让你拿到自己的影像后能直接照着做。2. 核心原理与数据准备理解定标公式背后那些参数2.1 绝对辐射定标的标准公式与参数含义先给出最常见的绝对辐射定标公式$$L DN \times Gain Offset$$其中L表示大气顶层的辐亮度单位一般是 W/(m²·sr·μm)DN就是影像原始记录的灰度值Gain是增益系数Offset是偏置系数。这两个系数称为定标系数由传感器在发射前实验室定标确定并在在轨运行期间通过定标场、月球观测、交叉定标等方式不断更新。对Landsat 8/9 OLI传感器来说Gain和Offset不会直接给数字而是以元数据文件里的RADIANCE_MULT_BAND_x和RADIANCE_ADD_BAND_x字段形式存在。x是波段号从1到11包含海岸波段和卷云波段。对国产卫星如高分一号PMS、资源三号等定标系数通常放在影像附属的XML文件里字段名可能是Gain、Offset也可能是CalibrationCoefficient不同载荷差异很大处理前一定要先翻元数据。需要注意的是有些老一辈遥感从业者习惯把公式写成L (DN - Offset) / Gain这是因为不同数据产品的定标系数定义方式不同。用之前先确认你的元数据里到底给的是“乘以再加”还是“先减再除”搞反了结果会差几倍甚至变成负值。2.2 从辐亮度到大气顶层反射率为什么要多走一步辐亮度是物理量但它还依赖太阳照射条件。同一个地物夏天太阳高、辐照度强辐亮度就大冬天太阳低辐亮度就小。这给不同季节影像对比又添了一层麻烦。所以实际应用中一般会把辐亮度进一步转化为大气顶层反射率Top of Atmosphere Reflectance简称TOA反射率也叫表观反射率。计算公式为$$\rho \frac{\pi \cdot L \cdot d^2}{ESUN \cdot cos(\theta)}$$先解释每个变量的意义d是日地距离修正因子单位是天文单位AU用来修正地球绕太阳公转导致的日地距离变化取值通常在0.983到1.017之间ESUN是大气顶层太阳光谱辐照度单位是W/(m²·μm)这个值对同一传感器是固定的在定标文件里可以查到θ是太阳天顶角注意是90度减去元数据里的太阳高度角Sun Elevation别直接用高度角代入余弦。我在实际处理中常规流程是直接输出TOA反射率作为标准产品。原因很简单辐亮度受太阳高度角影响而同一区域的时序影像太阳高度角每天都在变用TOA反射率才能消除这个变量。当然TOA反射率还包含大气散射和吸收的影响真正的地表反射率需要额外做大气校正比如6S、MODTRAN、FLAASH等。但TOA反射率作为绝对辐射定标的输出产物已经能满足大部分相对对比分析需求。2.3 基于Python的工具链与安装环境做定标计算Python生态里最顺手的组合是rasterionumpypandas再加一个matplotlib做结果可视化验证。rasterio负责读写GeoTIFFnumpy处理数组运算pandas整理波段清单和元数据。不同机器的安装我建议直接走conda因为rasterio对GDAL的依赖比较敏感用pip装偶尔会碰到proj库版本冲突conda能把这个麻烦降到最低。安装命令非常简单conda create -n rs_calibration python3.10 conda activate rs_calibration conda install rasterio numpy pandas matplotlib -c conda-forge如果你机器上已经装好了Python环境不打算再建虚拟环境直接用pip安装也可以pip install rasterio numpy pandas matplotlib装完之后可以快速验证一下rasterio是否能正常读取数据。我建议你在正式动手前先跑通环境避免处理到一半才发现GDAL相关依赖有问题。import rasterio print(rasterio.__version__)只要能输出版本号环境就算就绪了。3. 实操流程用Python对Landsat影像做绝对辐射定标3.1 读取元数据与定标系数一般Landsat 8/9 L1级数据包的MTL.txt文件里就存了我们需要的定标参数数据包内结构类似这样LC08_L1TP_123032_20220615_20220615_01_T1/ |-- LC08_L1TP_123032_20220615_20220615_01_T1_MTL.txt |-- LC08_L1TP_123032_20220615_20220615_01_T1_B1.TIF |-- LC08_L1TP_123032_20220615_20220615_01_T1_B2.TIF |-- ...用Python读取MTL文件非常容易不需要额外库正则表达式就能搞定核心字段。我一般这样写import re def read_mtl(mtl_path): with open(mtl_path, r) as f: content f.read() mult_dict {} add_dict {} # 提取每个波段的增益和偏置 for band in range(1, 12): mult_pattern rfRADIANCE_MULT_BAND_{band}\s*\s*([\d.Ee-]) add_pattern rfRADIANCE_ADD_BAND_{band}\s*\s*([\d.Ee-]) mult_match re.search(mult_pattern, content) add_match re.search(add_pattern, content) if mult_match and add_match: mult_dict[band] float(mult_match.group(1)) add_dict[band] float(add_match.group(1)) # 提取太阳高度角和日地距离 sun_elevation float(re.search(rSUN_ELEVATION\s*\s*([\d.Ee-]), content).group(1)) earth_sun_distance float(re.search(rEARTH_SUN_DISTANCE\s*\s*([\d.Ee-]), content).group(1)) return mult_dict, add_dict, sun_elevation, earth_sun_distance这段代码的逻辑是把MTL里每个波段的两个关键系数抽出来同时把太阳高度角和日地距离也一并取到。日地距离这个字段在Landsat 8/9的MTL里叫EARTH_SUN_DISTANCE单位就是AU。如果你的数据是Landsat之前的旧版本如Landsat 5/7可能没有这个字段那需要用经验公式计算后面我会补充。3.2 单波段辐射定标的核心代码拿到定标系数后用rasterio把波段影像读成numpy数组再套公式转换成辐亮度和TOA反射率最后写回GeoTIFF。下面这段代码基本上可以直接复用import numpy as np import rasterio import math def apply_rad_calibration(band_path, output_path, gain, offset, sun_elevation, earth_sun_distance, esun, wavelength_band): with rasterio.open(band_path) as src: dn src.read(1).astype(np.float32) profile src.profile.copy() # 处理填充值Landsat的填充值为0 dn[dn 0] np.nan # 某些数据里填充值是65535或-9999根据实际情况调整 # dn[dn 65535] np.nan # 绝对辐射定标DN - 辐亮度 radiance dn * gain offset # 辐亮度 - TOA反射率 cos_theta math.cos(math.radians(90 - sun_elevation)) reflectance (math.pi * radiance * earth_sun_distance ** 2) / (esun[wavelength_band] * cos_theta) # 将无效值转回0或者填nan取决于后续需求 reflectance[np.isnan(reflectance)] 0 # 反射率理论范围0-1将异常值裁剪 reflectance np.clip(reflectance, 0.0, 1.0) profile.update(dtyperasterio.float32, count1, compresslzw, nodata0) with rasterio.open(output_path, w, **profile) as dst: dst.write(reflectance, 1) return reflectance有几个细节值得单独说明。第一读进来的DN数组必须转成float32再做运算不然整型数组计算会截断小数。第二填充值不处理的话会被算成一个异常大的反射率导出来的图全是麻点。Landsat L1级数据的填充值通常是0但你也可能面对其他处理级别的数据务必先检查直方图确认填充值是多少。第三clip到0-1区间是因为TOA反射率物理上不可能超过1出现超过1的情况通常意味着定标系数用错了或者大气条件极端。ESUN的值不同传感器不一样。Landsat 8 OLI的ESUN在USGS官方文档里有表格Landsat 5 TM、Landsat 7 ETM也各有自己的ESUN表。我习惯把这些值存成一个字典方便多波段循环时查表调用。这里给出Landsat 8 OLI的ESUN参考值单位W/(m²·μm)分别为波段1到7esun_oli { 1: 2067.0, 2: 1893.0, 3: 1603.0, 4: 972.6, 5: 245.0, 6: 79.72, 7: 33.99 }需要强调ESUN这类参数在学术圈有不同版本不同文献给的数值微有差别但对大多数应用场景来说差异极其微小不会对结果产生实质性影响。3.3 日地距离修正因子的两种获取方式Landsat 8/9的新版本MTL里会直接给EARTH_SUN_DISTANCE字段直接用就行。但如果面对Landsat 5/7或一些国产卫星数据元数据没给这个参数就得自己算。一个经典的经验公式是$$d 1 - 0.0167 \cdot \cos(2\pi \cdot (J-3) / 365)$$其中J是儒略日Julian Day也就是这一年在第几天。1月1日是第1天。用Python计算非常简洁from datetime import datetime def julian_day(date_str): dt datetime.strptime(date_str, %Y-%m-%d) return dt.timetuple().tm_yday def earth_sun_distance(jday): return 1 - 0.0167 * np.cos(2 * np.pi * (jday - 3) / 365)这个公式算出来的d误差在0.0001 AU级别对定量分析完全够用。需要小心的是有些国产卫星的元数据文件名里带的日期字符串格式可能是“20220615”而不是“2022-06-15”解析前先做格式转换。3.4 批量处理整个数据包多波段循环与输出管理真实项目中很少只处理单个波段。一个Landsat场景至少要用到可见光到短波红外的六七个波段国产高分系列也类似。如果每个波段都手写一遍上面那些代码不仅重复劳动还容易在复制粘贴时把增益系数搞混。我通常用外层循环一次性批量处理import os from glob import glob def batch_calibrate(data_dir, mtl_path, output_dir, band_list): mult_dict, add_dict, sun_elev, earth_dist read_mtl(mtl_path) esun { 1: 2067.0, 2: 1893.0, 3: 1603.0, 4: 972.6, 5: 245.0, 6: 79.72, 7: 33.99 } os.makedirs(output_dir, exist_okTrue) for band in band_list: band_files glob(os.path.join(data_dir, f*_B{band}.TIF)) if not band_files: print(f波段{band}文件未找到跳过) continue band_path band_files[0] output_path os.path.join(output_dir, fTOA_B{band}.tif) apply_rad_calibration( band_path, output_path, gainmult_dict[band], offsetadd_dict[band], sun_elevationsun_elev, earth_sun_distanceearth_dist, esunesun, wavelength_bandband ) print(f波段{band}定标完成 - {output_path})这段代码中band_list你可以只给需要的波段比如做植被分析就给4、5波段红波段和近红外这样能大幅减少磁盘占用和计算时间。处理完的结果是多个单波段GeoTIFF后续如果需要合成RGB或者计算NDVI直接用rasterio再合并或者做波段运算即可。这里提一句多波段定标不推荐一次性把所有波段读进内存再算尤其是一景Landst 8全波段有11个波段每个1.5GB左右16位整型全部读入内存很容易爆掉。逐个波段流式处理更稳健。4. 常见问题与实战排查技巧4.1 反射率全部为负或超过1这是定标最常翻车的地方。反射率为负大概率是DN值里包含了填充值且没剔除干净或者Offset大于DN乘以Gain的结果——本质是填了无效像元。反射率超过1常见于ESUN查表错误、太阳天顶角用错了角度制或者元数据里的太阳高度角没转换成天顶角就代入了余弦计算。我调试时习惯先画一个直方图看看反射率分布范围。正常的TOA反射率分布应该在0到0.5之间集中水体在0.05以下植被在0.3左右超过1就一定有bug。画图用matplotlibimport matplotlib.pyplot as plt plt.hist(reflectance[reflectance 0].ravel(), bins100, range(0, 1.2)) plt.xlabel(TOA Reflectance) plt.ylabel(Pixel Count) plt.show()如果直方图在0.999附近出现一个尖峰说明大量像元被clip到了1.0往往代表增益系数或太阳高度角的符号出了问题。4.2 批量处理中元数据字段缺失国产卫星数据元数据格式五花八门。有些高分系列的XML文件里增益系数命名可能是CalibrationCoefficient甚至同一个卫星的A/B相机命名还不一样。碰到这种情况我建议直接打开XML文件人工核对字段名再针对性写解析规则。不要赌字段名更不要在没确认字段结构的情况下批量跑几百景数据那样出了问题排查成本极高。另外有些L1级数据处理软件比如某些国产预处理系统会在文件名或附属文件中标注“已定标”或“表观反射率产品”这类数据拿到手就已经是TOA反射率不需要再做一次绝对辐射定标。判断方法很简单打开原文件头信息看单位如果单位是reflectance或者percent reflectance就直接跳过多算一步的环节如果单位还是DN或者count再执行定标。4.3 输出文件精度选择与磁盘空间规划我见过不少人定标完直接把结果存成uint16类型结果就是浮点型的反射率被四舍五入成0和1两个整数整个图像变成黑白二值图。绝对辐射定标的结果必须保存为浮点型建议float32这样能保证反射率的小数精度。还有一点容易被忽略L1级原始数据通常是16位整型单波段文件体积可能在150MB左右定标后转成float32文件体积直接翻倍到300MB左右。如果一景影像做8个波段的批量处理输出就是2.4GB。时间序列动辄几十上百景存储规划要在处理前就想清楚。压缩推荐用LZW无损压缩并且对浮点型遥感数据压缩率很高有时能压到原来的一半以下。输出路径组织方面我习惯这样管理output_dir/ |-- LC08_20220615/ | |-- TOA_B2.tif | |-- TOA_B3.tif | |-- TOA_B4.tif | |-- TOA_B5.tif |-- LC08_20220701/ | |-- TOA_B2.tif | |-- ...一个场景一个文件夹文件夹名带上日期和轨道号后面做时间序列管理会非常省事。4.4 定标结果叠加显示时的投影与范围不一致单波段定标结果如果直接用来做RGB合成要确保每个波段的空间范围和分辨率完全一致。Landsat L1级产品理论上各波段已经配准到一个网格了但实际中偶尔会遇到某景影像部分波段范围偏移半个像元。用rasterio打开两个波段对比一下transform和bounds就能发现。如果有细微偏移用rasterio的reproject或者warp功能重采样对齐统一到参考波段的网格上再做合成。在我自己的流程中定标完成后的第一步永远是先目检结果。用rasterio读入红、绿、蓝三个波段的TOA反射率用matplotlib的imshow拉伸显示一遍看看地物纹理是否自然、水体是否呈暗色、植被是否突出。这一步花不了半分钟能挡住98%的明显错误。别急着往下游流程走等发现问题再回头排查成本会高好几倍。5. 进阶经验从定标到应用的延伸5.1 定标后的NDVI计算示例定标完成后TOA反射率可以直接用于计算各种光谱指数。拿NDVI举例需要红波段Landsat 8对应的波段4和近红外波段波段5代码极短但很能说明问题import rasterio import numpy as np def calc_ndvi(red_path, nir_path, output_path): with rasterio.open(red_path) as src_red: red src_red.read(1).astype(np.float32) profile src_red.profile.copy() with rasterio.open(nir_path) as src_nir: nir src_nir.read(1).astype(np.float32) ndvi (nir - red) / (nir red 1e-10) ndvi np.clip(ndvi, -1, 1) profile.update(dtyperasterio.float32, count1, compresslzw, nodata-9999) with rasterio.open(output_path, w, **profile) as dst: dst.write(ndvi, 1)分母加一个1e-10是为了避免除以零。如果直接用没有定标过的DN值算NDVI比值虽然能约掉一部分增益影响但大气和太阳高度角的影响依然存在做时序对比时误差会非常明显。定标后再算指数不同时相的指数才具备可比性。5.2 Python环境与性能优化心得在处理大范围高分辨率影像时例如高分一号PMS全色波段分辨率为2米一景全色影像可能达到2GB以上。直接用numpy整幅读取有时会内存告急这时候可以用rasterio的窗口读取方式分块处理with rasterio.open(band_path) as src: for ji, window in src.block_windows(1): dn_block src.read(1, windowwindow) # 对block做定标计算 radiance_block dn_block * gain offset # 写回输出文件的对应窗口 dst.write(radiance_block, 1, windowwindow)这样能把内存峰值降下来一个数量级。此外如果每块的计算都是纯数组操作且波段之间没有依赖还可以用concurrent.futures做线程池并行多核CPU环境下提速显著。我个人在使用中发现rasterio本身在读写大文件时性能瓶颈往往在磁盘IO而不在计算。使用SSD存储数据、处理过程中避免同时读写同一块磁盘、及时关闭文件句柄这些看似不起眼的习惯对整体效率影响非常大。建议批处理时先在一个小范围测试场景上跑通全流程再扩展到全场景批处理能少踩很多无谓的坑。5.3 关于定标系数版本更新的提醒定标系数不是永远不变的。卫星传感器在轨运行会逐渐老化探测元件响应度下降导致同样的地物反射输入DN输出值逐渐漂移。所以传感器定标团队会定期更新定标系数比如Landsat 8每次数据处理版本更新定标系数都可能变化。因此建议你从USGS官网下载L1级数据时注意数据包的处理级别和生成时间用MTL里自带的最新系数不要拿着2015年某个教程里的固定系数套用到2024年的数据上。国产卫星数据同样如此各省中心或高校接收的高分数据定标系数通常每季度或每半年更新一次从数据接收方渠道获取最新版本最可靠。此外在做多期影像对比时尽量统一数据来源和处理级别。不要混用同一传感器的不同处理级别比如L1T和L2A混着用否则定标系数差异可能造成虚假的时序变化信号。这几个细节注意好你的定标结果才能经得起审稿人或者项目验收的推敲。6. 写在最后的几点实践体会做遥感数据处理这几年我最大的感受是绝对辐射定标这个步骤看似基础却决定了后面所有分析结果的可信度。很多人在数据集上花了大量时间调模型参数却不舍得花几分钟把影像定标做对最后模型精度上不去还以为是算法问题。如果你正被光谱指数异常、时序曲线跳变这类问题困扰不妨先回头检查一下定标环节。另一个想分享的小技巧是做任何批处理前先做一景完整流程并保存一个验证脚本。把输入影像的行列数、DN值范围、定标系数、输出反射率均值都打印出来归档留档。这样下次处理新数据时如果结果跟历史统计差异较大能迅速定位是数据源问题还是处理流程问题。这个习惯救过我很多次也推荐你建立起来。如果你手头有国产高分、资源系列或其他非Landsat数据定标原理一样只是元数据解析部分需要按各载荷的格式调整。最笨但最有效的方法就是打开元数据文件一项项看不猜、不赌、不凭印象处理。按照本文的方法论走一遍你大概率能顺利产出第一批正确的定标产品。本文还有配套的精品资源点击获取
简介遥感图像的绝对辐射定标是将原始DN值转换为地表反射率或辐射亮度等物理量的关键预处理步骤。这一Python工程面向遥感初学者、地学研究人员及GIS开发人员演示如何利用GDAL库打开多波段影像、从JSON文件中读取增益与暗电流等定标参数并完成逐像元辐射校正同时给出大气校正与元数据整理的思路可结合实际任务灵活调整。压缩包共9个文件核心包括CrossRadiationCalibration.py可执行脚本和RadiometricCorrectionParameter.json参数文件另含工程配置文件与XML说明便于理解整个工程结构整体仅16KB轻量且易于阅读。目前已有896人学习下载适合希望快速掌握遥感辐射定标编程实现、理解传感器定标系数应用逻辑的入门者。借助该资源使用者能梳理从数据准备、参数解析、辐射值计算到结果保存的完整流程并基于示例代码扩展自己的处理脚本有效降低从理论到实践的试错成本。1. 概念拆解为什么绝对辐射定标是遥感定量化的第一道门槛拿到一景L1级影像很多人第一反应是直接拿DN值去算NDVI。我之前也这么干过直到一次做长时间序列的植被指数变化分析发现同一条样带的NDVI曲线在季度交界处出现莫名其妙的跳变。排查到最后问题出在影像之间的辐射基准不一致——说白了DN值只是传感器把接收到的能量量化成的一个整数它跟地物真实反射太阳辐射的能力之间还隔着一层“翻译”过程。这层翻译就叫绝对辐射定标。绝对辐射定标解决的问题是把影像记录的DN值转换成具有物理意义的辐亮度或者反射率。DN值的大小不仅受地物本身反射率影响还跟太阳高度角、传感器增益设置、大气路径辐射等因素有关。如果不做定标两景不同时间、不同角度、甚至不同传感器获取的影像根本没法在数值层面直接对比。这一点在定量遥感里是底线在目视解译里无所谓但在做反演、做光谱分析、做机器学习样本标注时踩不踩这一步结果能差出一个量级。用生活化的方式理解DN值像是一张照片上像素的明暗程度同样的白墙中午拍和傍晚拍亮度完全不一样辐亮度则是把相机换成照度计测的是这面墙实际反射了多少光能量。绝对辐射定标就是帮你把“照片的明暗”换算成“物理上的亮度”后面再做大气校正、反射率反演、地表温度反演才有共同的语言基础。本文面对的读者我默认是遥感或地学相关专业的学生、刚入行做遥感数据处理的技术人员以及需要在Python里批量处理遥感影像的开发者。我会从原理讲到代码实现再到坑点排查全程用Landsat系列和国产高分系列的数据做例子尽量把每一步讲透让你拿到自己的影像后能直接照着做。2. 核心原理与数据准备理解定标公式背后那些参数2.1 绝对辐射定标的标准公式与参数含义先给出最常见的绝对辐射定标公式$$L DN \times Gain Offset$$其中L表示大气顶层的辐亮度单位一般是 W/(m²·sr·μm)DN就是影像原始记录的灰度值Gain是增益系数Offset是偏置系数。这两个系数称为定标系数由传感器在发射前实验室定标确定并在在轨运行期间通过定标场、月球观测、交叉定标等方式不断更新。对Landsat 8/9 OLI传感器来说Gain和Offset不会直接给数字而是以元数据文件里的RADIANCE_MULT_BAND_x和RADIANCE_ADD_BAND_x字段形式存在。x是波段号从1到11包含海岸波段和卷云波段。对国产卫星如高分一号PMS、资源三号等定标系数通常放在影像附属的XML文件里字段名可能是Gain、Offset也可能是CalibrationCoefficient不同载荷差异很大处理前一定要先翻元数据。需要注意的是有些老一辈遥感从业者习惯把公式写成L (DN - Offset) / Gain这是因为不同数据产品的定标系数定义方式不同。用之前先确认你的元数据里到底给的是“乘以再加”还是“先减再除”搞反了结果会差几倍甚至变成负值。2.2 从辐亮度到大气顶层反射率为什么要多走一步辐亮度是物理量但它还依赖太阳照射条件。同一个地物夏天太阳高、辐照度强辐亮度就大冬天太阳低辐亮度就小。这给不同季节影像对比又添了一层麻烦。所以实际应用中一般会把辐亮度进一步转化为大气顶层反射率Top of Atmosphere Reflectance简称TOA反射率也叫表观反射率。计算公式为$$\rho \frac{\pi \cdot L \cdot d^2}{ESUN \cdot cos(\theta)}$$先解释每个变量的意义d是日地距离修正因子单位是天文单位AU用来修正地球绕太阳公转导致的日地距离变化取值通常在0.983到1.017之间ESUN是大气顶层太阳光谱辐照度单位是W/(m²·μm)这个值对同一传感器是固定的在定标文件里可以查到θ是太阳天顶角注意是90度减去元数据里的太阳高度角Sun Elevation别直接用高度角代入余弦。我在实际处理中常规流程是直接输出TOA反射率作为标准产品。原因很简单辐亮度受太阳高度角影响而同一区域的时序影像太阳高度角每天都在变用TOA反射率才能消除这个变量。当然TOA反射率还包含大气散射和吸收的影响真正的地表反射率需要额外做大气校正比如6S、MODTRAN、FLAASH等。但TOA反射率作为绝对辐射定标的输出产物已经能满足大部分相对对比分析需求。2.3 基于Python的工具链与安装环境做定标计算Python生态里最顺手的组合是rasterionumpypandas再加一个matplotlib做结果可视化验证。rasterio负责读写GeoTIFFnumpy处理数组运算pandas整理波段清单和元数据。不同机器的安装我建议直接走conda因为rasterio对GDAL的依赖比较敏感用pip装偶尔会碰到proj库版本冲突conda能把这个麻烦降到最低。安装命令非常简单conda create -n rs_calibration python3.10 conda activate rs_calibration conda install rasterio numpy pandas matplotlib -c conda-forge如果你机器上已经装好了Python环境不打算再建虚拟环境直接用pip安装也可以pip install rasterio numpy pandas matplotlib装完之后可以快速验证一下rasterio是否能正常读取数据。我建议你在正式动手前先跑通环境避免处理到一半才发现GDAL相关依赖有问题。import rasterio print(rasterio.__version__)只要能输出版本号环境就算就绪了。3. 实操流程用Python对Landsat影像做绝对辐射定标3.1 读取元数据与定标系数一般Landsat 8/9 L1级数据包的MTL.txt文件里就存了我们需要的定标参数数据包内结构类似这样LC08_L1TP_123032_20220615_20220615_01_T1/ |-- LC08_L1TP_123032_20220615_20220615_01_T1_MTL.txt |-- LC08_L1TP_123032_20220615_20220615_01_T1_B1.TIF |-- LC08_L1TP_123032_20220615_20220615_01_T1_B2.TIF |-- ...用Python读取MTL文件非常容易不需要额外库正则表达式就能搞定核心字段。我一般这样写import re def read_mtl(mtl_path): with open(mtl_path, r) as f: content f.read() mult_dict {} add_dict {} # 提取每个波段的增益和偏置 for band in range(1, 12): mult_pattern rfRADIANCE_MULT_BAND_{band}\s*\s*([\d.Ee-]) add_pattern rfRADIANCE_ADD_BAND_{band}\s*\s*([\d.Ee-]) mult_match re.search(mult_pattern, content) add_match re.search(add_pattern, content) if mult_match and add_match: mult_dict[band] float(mult_match.group(1)) add_dict[band] float(add_match.group(1)) # 提取太阳高度角和日地距离 sun_elevation float(re.search(rSUN_ELEVATION\s*\s*([\d.Ee-]), content).group(1)) earth_sun_distance float(re.search(rEARTH_SUN_DISTANCE\s*\s*([\d.Ee-]), content).group(1)) return mult_dict, add_dict, sun_elevation, earth_sun_distance这段代码的逻辑是把MTL里每个波段的两个关键系数抽出来同时把太阳高度角和日地距离也一并取到。日地距离这个字段在Landsat 8/9的MTL里叫EARTH_SUN_DISTANCE单位就是AU。如果你的数据是Landsat之前的旧版本如Landsat 5/7可能没有这个字段那需要用经验公式计算后面我会补充。3.2 单波段辐射定标的核心代码拿到定标系数后用rasterio把波段影像读成numpy数组再套公式转换成辐亮度和TOA反射率最后写回GeoTIFF。下面这段代码基本上可以直接复用import numpy as np import rasterio import math def apply_rad_calibration(band_path, output_path, gain, offset, sun_elevation, earth_sun_distance, esun, wavelength_band): with rasterio.open(band_path) as src: dn src.read(1).astype(np.float32) profile src.profile.copy() # 处理填充值Landsat的填充值为0 dn[dn 0] np.nan # 某些数据里填充值是65535或-9999根据实际情况调整 # dn[dn 65535] np.nan # 绝对辐射定标DN - 辐亮度 radiance dn * gain offset # 辐亮度 - TOA反射率 cos_theta math.cos(math.radians(90 - sun_elevation)) reflectance (math.pi * radiance * earth_sun_distance ** 2) / (esun[wavelength_band] * cos_theta) # 将无效值转回0或者填nan取决于后续需求 reflectance[np.isnan(reflectance)] 0 # 反射率理论范围0-1将异常值裁剪 reflectance np.clip(reflectance, 0.0, 1.0) profile.update(dtyperasterio.float32, count1, compresslzw, nodata0) with rasterio.open(output_path, w, **profile) as dst: dst.write(reflectance, 1) return reflectance有几个细节值得单独说明。第一读进来的DN数组必须转成float32再做运算不然整型数组计算会截断小数。第二填充值不处理的话会被算成一个异常大的反射率导出来的图全是麻点。Landsat L1级数据的填充值通常是0但你也可能面对其他处理级别的数据务必先检查直方图确认填充值是多少。第三clip到0-1区间是因为TOA反射率物理上不可能超过1出现超过1的情况通常意味着定标系数用错了或者大气条件极端。ESUN的值不同传感器不一样。Landsat 8 OLI的ESUN在USGS官方文档里有表格Landsat 5 TM、Landsat 7 ETM也各有自己的ESUN表。我习惯把这些值存成一个字典方便多波段循环时查表调用。这里给出Landsat 8 OLI的ESUN参考值单位W/(m²·μm)分别为波段1到7esun_oli { 1: 2067.0, 2: 1893.0, 3: 1603.0, 4: 972.6, 5: 245.0, 6: 79.72, 7: 33.99 }需要强调ESUN这类参数在学术圈有不同版本不同文献给的数值微有差别但对大多数应用场景来说差异极其微小不会对结果产生实质性影响。3.3 日地距离修正因子的两种获取方式Landsat 8/9的新版本MTL里会直接给EARTH_SUN_DISTANCE字段直接用就行。但如果面对Landsat 5/7或一些国产卫星数据元数据没给这个参数就得自己算。一个经典的经验公式是$$d 1 - 0.0167 \cdot \cos(2\pi \cdot (J-3) / 365)$$其中J是儒略日Julian Day也就是这一年在第几天。1月1日是第1天。用Python计算非常简洁from datetime import datetime def julian_day(date_str): dt datetime.strptime(date_str, %Y-%m-%d) return dt.timetuple().tm_yday def earth_sun_distance(jday): return 1 - 0.0167 * np.cos(2 * np.pi * (jday - 3) / 365)这个公式算出来的d误差在0.0001 AU级别对定量分析完全够用。需要小心的是有些国产卫星的元数据文件名里带的日期字符串格式可能是“20220615”而不是“2022-06-15”解析前先做格式转换。3.4 批量处理整个数据包多波段循环与输出管理真实项目中很少只处理单个波段。一个Landsat场景至少要用到可见光到短波红外的六七个波段国产高分系列也类似。如果每个波段都手写一遍上面那些代码不仅重复劳动还容易在复制粘贴时把增益系数搞混。我通常用外层循环一次性批量处理import os from glob import glob def batch_calibrate(data_dir, mtl_path, output_dir, band_list): mult_dict, add_dict, sun_elev, earth_dist read_mtl(mtl_path) esun { 1: 2067.0, 2: 1893.0, 3: 1603.0, 4: 972.6, 5: 245.0, 6: 79.72, 7: 33.99 } os.makedirs(output_dir, exist_okTrue) for band in band_list: band_files glob(os.path.join(data_dir, f*_B{band}.TIF)) if not band_files: print(f波段{band}文件未找到跳过) continue band_path band_files[0] output_path os.path.join(output_dir, fTOA_B{band}.tif) apply_rad_calibration( band_path, output_path, gainmult_dict[band], offsetadd_dict[band], sun_elevationsun_elev, earth_sun_distanceearth_dist, esunesun, wavelength_bandband ) print(f波段{band}定标完成 - {output_path})这段代码中band_list你可以只给需要的波段比如做植被分析就给4、5波段红波段和近红外这样能大幅减少磁盘占用和计算时间。处理完的结果是多个单波段GeoTIFF后续如果需要合成RGB或者计算NDVI直接用rasterio再合并或者做波段运算即可。这里提一句多波段定标不推荐一次性把所有波段读进内存再算尤其是一景Landst 8全波段有11个波段每个1.5GB左右16位整型全部读入内存很容易爆掉。逐个波段流式处理更稳健。4. 常见问题与实战排查技巧4.1 反射率全部为负或超过1这是定标最常翻车的地方。反射率为负大概率是DN值里包含了填充值且没剔除干净或者Offset大于DN乘以Gain的结果——本质是填了无效像元。反射率超过1常见于ESUN查表错误、太阳天顶角用错了角度制或者元数据里的太阳高度角没转换成天顶角就代入了余弦计算。我调试时习惯先画一个直方图看看反射率分布范围。正常的TOA反射率分布应该在0到0.5之间集中水体在0.05以下植被在0.3左右超过1就一定有bug。画图用matplotlibimport matplotlib.pyplot as plt plt.hist(reflectance[reflectance 0].ravel(), bins100, range(0, 1.2)) plt.xlabel(TOA Reflectance) plt.ylabel(Pixel Count) plt.show()如果直方图在0.999附近出现一个尖峰说明大量像元被clip到了1.0往往代表增益系数或太阳高度角的符号出了问题。4.2 批量处理中元数据字段缺失国产卫星数据元数据格式五花八门。有些高分系列的XML文件里增益系数命名可能是CalibrationCoefficient甚至同一个卫星的A/B相机命名还不一样。碰到这种情况我建议直接打开XML文件人工核对字段名再针对性写解析规则。不要赌字段名更不要在没确认字段结构的情况下批量跑几百景数据那样出了问题排查成本极高。另外有些L1级数据处理软件比如某些国产预处理系统会在文件名或附属文件中标注“已定标”或“表观反射率产品”这类数据拿到手就已经是TOA反射率不需要再做一次绝对辐射定标。判断方法很简单打开原文件头信息看单位如果单位是reflectance或者percent reflectance就直接跳过多算一步的环节如果单位还是DN或者count再执行定标。4.3 输出文件精度选择与磁盘空间规划我见过不少人定标完直接把结果存成uint16类型结果就是浮点型的反射率被四舍五入成0和1两个整数整个图像变成黑白二值图。绝对辐射定标的结果必须保存为浮点型建议float32这样能保证反射率的小数精度。还有一点容易被忽略L1级原始数据通常是16位整型单波段文件体积可能在150MB左右定标后转成float32文件体积直接翻倍到300MB左右。如果一景影像做8个波段的批量处理输出就是2.4GB。时间序列动辄几十上百景存储规划要在处理前就想清楚。压缩推荐用LZW无损压缩并且对浮点型遥感数据压缩率很高有时能压到原来的一半以下。输出路径组织方面我习惯这样管理output_dir/ |-- LC08_20220615/ | |-- TOA_B2.tif | |-- TOA_B3.tif | |-- TOA_B4.tif | |-- TOA_B5.tif |-- LC08_20220701/ | |-- TOA_B2.tif | |-- ...一个场景一个文件夹文件夹名带上日期和轨道号后面做时间序列管理会非常省事。4.4 定标结果叠加显示时的投影与范围不一致单波段定标结果如果直接用来做RGB合成要确保每个波段的空间范围和分辨率完全一致。Landsat L1级产品理论上各波段已经配准到一个网格了但实际中偶尔会遇到某景影像部分波段范围偏移半个像元。用rasterio打开两个波段对比一下transform和bounds就能发现。如果有细微偏移用rasterio的reproject或者warp功能重采样对齐统一到参考波段的网格上再做合成。在我自己的流程中定标完成后的第一步永远是先目检结果。用rasterio读入红、绿、蓝三个波段的TOA反射率用matplotlib的imshow拉伸显示一遍看看地物纹理是否自然、水体是否呈暗色、植被是否突出。这一步花不了半分钟能挡住98%的明显错误。别急着往下游流程走等发现问题再回头排查成本会高好几倍。5. 进阶经验从定标到应用的延伸5.1 定标后的NDVI计算示例定标完成后TOA反射率可以直接用于计算各种光谱指数。拿NDVI举例需要红波段Landsat 8对应的波段4和近红外波段波段5代码极短但很能说明问题import rasterio import numpy as np def calc_ndvi(red_path, nir_path, output_path): with rasterio.open(red_path) as src_red: red src_red.read(1).astype(np.float32) profile src_red.profile.copy() with rasterio.open(nir_path) as src_nir: nir src_nir.read(1).astype(np.float32) ndvi (nir - red) / (nir red 1e-10) ndvi np.clip(ndvi, -1, 1) profile.update(dtyperasterio.float32, count1, compresslzw, nodata-9999) with rasterio.open(output_path, w, **profile) as dst: dst.write(ndvi, 1)分母加一个1e-10是为了避免除以零。如果直接用没有定标过的DN值算NDVI比值虽然能约掉一部分增益影响但大气和太阳高度角的影响依然存在做时序对比时误差会非常明显。定标后再算指数不同时相的指数才具备可比性。5.2 Python环境与性能优化心得在处理大范围高分辨率影像时例如高分一号PMS全色波段分辨率为2米一景全色影像可能达到2GB以上。直接用numpy整幅读取有时会内存告急这时候可以用rasterio的窗口读取方式分块处理with rasterio.open(band_path) as src: for ji, window in src.block_windows(1): dn_block src.read(1, windowwindow) # 对block做定标计算 radiance_block dn_block * gain offset # 写回输出文件的对应窗口 dst.write(radiance_block, 1, windowwindow)这样能把内存峰值降下来一个数量级。此外如果每块的计算都是纯数组操作且波段之间没有依赖还可以用concurrent.futures做线程池并行多核CPU环境下提速显著。我个人在使用中发现rasterio本身在读写大文件时性能瓶颈往往在磁盘IO而不在计算。使用SSD存储数据、处理过程中避免同时读写同一块磁盘、及时关闭文件句柄这些看似不起眼的习惯对整体效率影响非常大。建议批处理时先在一个小范围测试场景上跑通全流程再扩展到全场景批处理能少踩很多无谓的坑。5.3 关于定标系数版本更新的提醒定标系数不是永远不变的。卫星传感器在轨运行会逐渐老化探测元件响应度下降导致同样的地物反射输入DN输出值逐渐漂移。所以传感器定标团队会定期更新定标系数比如Landsat 8每次数据处理版本更新定标系数都可能变化。因此建议你从USGS官网下载L1级数据时注意数据包的处理级别和生成时间用MTL里自带的最新系数不要拿着2015年某个教程里的固定系数套用到2024年的数据上。国产卫星数据同样如此各省中心或高校接收的高分数据定标系数通常每季度或每半年更新一次从数据接收方渠道获取最新版本最可靠。此外在做多期影像对比时尽量统一数据来源和处理级别。不要混用同一传感器的不同处理级别比如L1T和L2A混着用否则定标系数差异可能造成虚假的时序变化信号。这几个细节注意好你的定标结果才能经得起审稿人或者项目验收的推敲。6. 写在最后的几点实践体会做遥感数据处理这几年我最大的感受是绝对辐射定标这个步骤看似基础却决定了后面所有分析结果的可信度。很多人在数据集上花了大量时间调模型参数却不舍得花几分钟把影像定标做对最后模型精度上不去还以为是算法问题。如果你正被光谱指数异常、时序曲线跳变这类问题困扰不妨先回头检查一下定标环节。另一个想分享的小技巧是做任何批处理前先做一景完整流程并保存一个验证脚本。把输入影像的行列数、DN值范围、定标系数、输出反射率均值都打印出来归档留档。这样下次处理新数据时如果结果跟历史统计差异较大能迅速定位是数据源问题还是处理流程问题。这个习惯救过我很多次也推荐你建立起来。如果你手头有国产高分、资源系列或其他非Landsat数据定标原理一样只是元数据解析部分需要按各载荷的格式调整。最笨但最有效的方法就是打开元数据文件一项项看不猜、不赌、不凭印象处理。按照本文的方法论走一遍你大概率能顺利产出第一批正确的定标产品。本文还有配套的精品资源点击获取
