尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

6S辐射传输模型实战:大气校正参数设置与Python代码实现

6S辐射传输模型实战:大气校正参数设置与Python代码实现 简介这份资源是6S模型的操作说明文档面向从事定量遥感、大气校正与卫星遥感数据分析的科研人员和工程师帮助理解大气传输过程对可见光与近红外观测的影响。6S模型由法国大气光学实验室在5S基础上改进而来可模拟平面观测、高层目标及非朗伯反射等复杂情形。文档围绕吸收效应、散射效应、内在大气反射率、方向效应、大气校正方案以及吸收与散射的交互作用展开并配有计算机代码说明、输入输出示例与子程序描述便于读者掌握模型原理与使用方式。资源包内含1个PDF文件大小约664KB结构完整、便于查阅。目前已有343人学习下载适合需要系统了解6S模型大气传输机制、开展遥感数据大气校正与地表参数反演的研究者参考使用。1. 从一次气溶胶反演偏差说起6S模型到底在算什么做遥感定量反演的人大多踩过同一个坑同一景影像用不同的大气校正参数跑出来的地表反射率在蓝光波段能差出百分之十几。排查到最后问题往往不在传感器定标而在辐射传输这一环——气溶胶光学厚度、水汽柱含量、观测几何这些量怎么进模型直接决定了大气程辐射和透过率算得准不准。6SSecond Simulation of the Satellite Signal in the Solar Spectrum辐射传输模型就是干这件事的给定太阳—地表—传感器这条路径上的大气状态算出大气顶信号里有多少是程辐射、多少是地表反射经大气衰减后的贡献。它属于逐次散射近似加SOS方法的辐射传输求解器覆盖0.25到4微米波段支持均一朗伯体、非均一朗伯体、BRDF地表能输出大气校正系数、球面反照率、偏振分量等。对做Landsat、Sentinel-2、MODIS大气校正的人来说6S是绕不开的参考实现很多业务化校正链比如早期LEDAPS、部分国产卫星地面处理系统内部都嵌了它的算法。这一篇不讲公式推导讲的是怎么把6S跑起来、参数怎么设、大气传输过程在代码里对应哪几个量以及结果怎么验证。2. 6S模型的大气传输过程拆解与参数映射2.1 太阳—大气—地表—传感器这条链路里发生了什么6S把大气顶反射率拆成三部分路径反射率程辐射、地表反射经大气两次透过后的贡献、以及地表与大气多次反射的耦合项。用公式粗写就是 ρ_TOA ρ_path T_down·T_up·ρ_surf/(1−s·ρ_surf)其中s是大气球面反照率。大气传输过程的核心就是算清楚ρ_path、T_down、T_up、s这四个量它们由气溶胶光学厚度、气溶胶类型、水汽、臭氧、观测几何共同决定。理解这一点很关键很多人以为6S只是输入参数输出反射率其实它输出的是一组大气校正系数业务系统拿这组系数去逐像元反算地表反射率。所以参数设错误差会以乘性因子的形式传导到整景影像。2.2 输入参数分组几何、大气、光谱、地表6S的输入按功能分四组理解分组比死记参数名有用分组关键参数典型取值/说明几何太阳天顶角、方位角、观测天顶角、方位角、月日角度单位度方位角以太阳为参考大气气溶胶光学厚度AOT、气溶胶模式、水汽、臭氧AOT 550nm模式选大陆/海洋/城市光谱光谱条件、波段号或自定义响应函数内置传感器或自定义波段地表地表类型、BRDF参数均一朗伯体最常用几何参数里最容易错的是方位角定义。6S要求的是相对方位角即观测方位角减太阳方位角很多人直接填绝对方位角结果程辐射算偏。气溶胶模式的选择影响散射相函数城市型气溶胶单次散射反照率低程辐射会比大陆型小蓝光波段尤其明显。2.3 用Python封装6S的最小可跑示例6S官方是Fortran程序输入靠一个文本文件或交互式问答。实际工程里一般用Py6S这个Python封装它把参数对象化避免手写输入文件出错。from Py6S import * # 初始化6S对象 s SixS() # 几何参数太阳天顶角、观测天顶角、相对方位角、月日 s.geometry Geometry.User() s.geometry.solar_z 35.0 # 太阳天顶角单位度 s.geometry.solar_a 120.0 # 太阳方位角 s.geometry.view_z 10.0 # 观测天顶角 s.geometry.view_a 145.0 # 观测方位角 s.geometry.month 7 s.geometry.day 15 # 大气参数气溶胶光学厚度、气溶胶模式、水汽、臭氧 s.aero_profile AeroProfile.PredefinedType(AeroProfile.Continental) s.atmos_profile AtmosProfile.UserWaterAndOzone(2.5, 0.3) # 水汽2.5g/cm2臭氧0.3atm-cm s.aot550 0.2 # 550nm气溶胶光学厚度 # 光谱条件以Landsat 8 OLI蓝光波段为例 s.wavelength Wavelength(0.48) # 单位微米 # 地表均一朗伯体反射率0.15 s.ground_reflectance GroundReflectance.HomogeneousLambertian(0.15) # 运行 s.run() # 取结果大气校正系数与反射率分量 print(路径反射率:, s.outputs.path_radiance) print(大气顶反射率:, s.outputs.apparent_reflectance) print(地表反射率:, s.outputs.ground_reflectance) print(下行透过率:, s.outputs.transmittance_down) print(上行透过率:, s.outputs.transmittance_up) print(球面反照率:, s.outputs.spherical_albedo)这段代码的逻辑是先固定几何再设大气状态最后指定光谱和地表s.run()触发Fortran内核计算。参数说明上aot550是550nm处的气溶胶光学厚度业务上一般来自MODIS或AERONETUserWaterAndOzone两个参数分别是水汽柱含量g/cm²和臭氧柱含量atm-cm。Wavelength(0.48)是单波长计算如果要模拟整个波段响应得用Wavelength配合传感器响应函数做积分。提示Py6S安装依赖Fortran编译器和6S源码Windows下建议用conda装py6sLinux下先编译6S再pip install Py6S否则run()会报找不到可执行文件。3. 用6S做大气校正的完整操作流程3.1 从影像元数据提取几何参数大气校正的第一步不是跑6S而是把影像的几何参数提出来。以Landsat 8为例元数据MTL文件里有SUN_ELEVATION和SUN_AZIMUTH需要转成天顶角import math sun_elev 55.3 # 来自MTL的SUN_ELEVATION sun_azim 120.0 # 来自MTL的SUN_AZIMUTH solar_z 90.0 - sun_elev # 天顶角 90 - 高度角 solar_a sun_azim # 方位角直接用 print(f太阳天顶角: {solar_z:.2f}, 太阳方位角: {solar_a:.2f})观测天顶角和方位角对Landsat这种近星下点传感器一般取0但Sentinel-2的视场角大边缘像元的观测天顶角能到10度以上必须逐像元处理。这一步的误差会直接进6S导致边缘和中心的大气校正结果不一致。3.2 气溶胶光学厚度的获取与插值AOT是6S里最敏感的参数。常见做法有三种用MODIS的MOD04产品插值到目标影像、用AERONET站点数据、或者用暗像元法从影像自身反演。工程上最稳的是MODIS插值import numpy as np from scipy.interpolate import griddata # MOD04的AOT点数据经纬度、AOT值 modis_lon np.array([116.1, 116.5, 116.9, 116.3]) modis_lat np.array([39.7, 39.9, 39.6, 40.1]) modis_aot np.array([0.18, 0.22, 0.25, 0.20]) # 目标影像中心经纬度 target_lon, target_lat 116.4, 39.8 # 反距离插值 aot griddata((modis_lon, modis_lat), modis_aot, (target_lon, target_lat), methodlinear) print(f插值AOT: {aot:.3f})逻辑说明MODIS AOT空间分辨率1km或10km目标影像可能30m插值到整景时如果AOT空间变化大建议分块插值而不是整景取一个均值。参数上methodlinear适合点分布均匀的情况点稀疏时用nearest更稳但会引入块状效应。3.3 逐波段运行6S并生成校正系数查找表业务化校正不会逐像元跑6S而是先生成一张查找表LUT把AOT、观测天顶角、地表反射率离散成网格每个格点跑一次6S存下校正系数再对影像逐像元查表插值。from Py6S import * import numpy as np aot_grid np.arange(0.05, 0.55, 0.05) # AOT离散 vza_grid np.arange(0, 15, 5) # 观测天顶角离散 lut {} for aot in aot_grid: for vza in vza_grid: s SixS() s.geometry Geometry.User() s.geometry.solar_z 35.0 s.geometry.view_z vza s.geometry.month 7 s.geometry.day 15 s.aot550 aot s.aero_profile AeroProfile.PredefinedType(AeroProfile.Continental) s.atmos_profile AtmosProfile.UserWaterAndOzone(2.5, 0.3) s.wavelength Wavelength(0.48) s.ground_reflectance GroundReflectance.HomogeneousLambertian(0.15) s.run() lut[(round(aot,2), vza)] { path: s.outputs.path_radiance, tdown: s.outputs.transmittance_down, tup: s.outputs.transmittance_up, s: s.outputs.spherical_albedo } print(LUT条目数:, len(lut))这段是LUT生成的核心。参数说明aot_grid步长0.05是精度和速度的折中AOT变化剧烈时缩到0.02vza_grid对Landsat可以只取0对宽视场传感器要加密。每个格点存四个量后续反算地表反射率时用ρ_surf (ρ_TOA − ρ_path)/(T_down·T_up s·(ρ_TOA − ρ_path))。注意LUT的维度不要盲目加AOT、VZA、水汽三维全离散格点数会指数增长。水汽对可见光影响小一般固定用影像过境时的再分析数据即可。3.4 反算地表反射率并检查异常值拿到LUT后对影像逐像元查表插值代入公式反算。反算完必须做异常值检查地表反射率出现负值或大于1说明程辐射扣多了或大气参数设错。def correct_toa(toa, path, tdown, tup, s): # 6S大气校正反算公式 rho_surf (toa - path) / (tdown * tup s * (toa - path)) return rho_surf # 示例某像元TOA反射率0.25查表得系数 rho correct_toa(0.25, 0.08, 0.85, 0.88, 0.12) print(f地表反射率: {rho:.4f}) # 异常值统计 if rho 0 or rho 1: print(异常检查AOT或几何参数)逻辑上path是程辐射tdown和tup是上下行透过率s是球面反照率。负值通常出现在蓝光波段且AOT设得过大时因为程辐射被高估。这时回查AOT来源或者检查气溶胶模式是否选错。4. 6S参数敏感性分析与常见报错排查4.1 哪些参数对结果影响最大不是所有参数都同等重要。做敏感性分析能帮你把精力放在关键量上参数变化范围对蓝光波段TOA的影响优先级AOT0.1→0.4程辐射增大约2倍高气溶胶模式大陆→城市程辐射差10%~20%高水汽1→3 g/cm²近红外影响明显中臭氧0.2→0.4 atm-cm绿光以上影响小低观测天顶角0→15度边缘像元差5%中AOT和气溶胶模式是绝对的重点。很多人AOT用对了但气溶胶模式默认用了大陆型而实际是城市污染气溶胶结果蓝光波段程辐射偏大校正后地表反射率偏低。4.2 运行6S时的典型报错与处理Py6S跑不起来八成是环境问题。常见报错和处理SixSException: Error running 6SFortran可执行文件路径不对检查PY6S_PATH环境变量或重新编译6S。ValueError: Wavelength must be between 0.25 and 4.0波长超范围6S只覆盖0.25~4微米紫外和热红外不能用。输出全为0或NaN几何参数里月日设成了0或者方位角填了负值6S对输入范围敏感。结果和预期差很多先查方位角是不是相对方位角再查AOT单位是不是550nm。# Linux下检查6S可执行文件是否就位 which sixs echo $PY6S_PATH # 手动跑一次6S看是否报错 sixs input.txt排查顺序建议从环境到参数先确认6S能独立运行再确认Py6S能调用最后才怀疑参数。参数问题里几何和AOT占九成。4.3 用AERONET数据验证6S输出验证是大气校正里最容易被跳过的一步。有AERONET站点的区域可以拿站点的AOT和实测地表反射率反推对比6S输出。# AERONET实测AOT与6S输入AOT对比 aeronet_aot 0.21 # 站点实测 sixs_aot 0.20 # 6S输入 diff abs(aeronet_aot - sixs_aot) / aeronet_aot print(fAOT相对偏差: {diff*100:.1f}%) # 偏差超过20%时考虑用AERONET值替换MODIS插值 if diff 0.2: print(建议用AERONET实测值重新跑6S)逻辑说明AERONET的AOT精度高于MODIS站点附近优先用实测值。偏差大说明MODIS插值或时空匹配有问题这时用AERONET值重跑LUT再对比校正后的地表反射率与地面实测形成闭环验证。5. 把6S嵌进批量处理链的进阶技巧单景跑通只是开始业务化要处理成百上千景。第一个技巧是LUT复用同一传感器、同一季节、同一区域几何和大气状态相近LUT可以跨景复用只对AOT做微调。把LUT存成HDF5或npz加载比重新跑6S快两个数量级。import numpy as np # 保存LUT np.savez(lut_landsat8_blue.npz, aot_gridaot_grid, vza_gridvza_grid, pathnp.array([lut[k][path] for k in lut]), tdownnp.array([lut[k][tdown] for k in lut])) # 加载复用 data np.load(lut_landsat8_blue.npz) print(复用LUTAOT网格:, data[aot_grid])第二个技巧是并行化。6S单次运行几十毫秒但LUT格点多时串行很慢用multiprocessing按AOT分块并行核数翻倍速度翻倍。第三个技巧是波段响应积分单波长计算和真实波段响应有差异宽波段传感器要用响应函数加权积分Py6S支持Wavelength配合PredefinedWavelengths做积分蓝光波段积分和单波长能差3%~5%。最后一个容易忽略的点6S输出的球面反照率s在多次散射强时不能忽略高反射地表如雪、云边缘必须带上s项否则反算的地表反射率会偏高。验证方法是拿已知反射率的地面目标如水泥地、水体对比偏差在5%以内说明参数链没问题。本文还有配套的精品资源点击获取
返回列表