辐射定标与大气校正

引言

遥感影像里的像素值本身没有物理意义。传感器记录的 DN 值只是探测器响应的相对强度，受增益、积分时间、探测器老化、太阳高度角、大气散射与吸收层层影响。要让不同时间、不同传感器、不同地点的影像可比，必须先把 DN 还原成地表反射率这个物理量，这条链路就是辐射定标与大气校正。

工程上的难点有三处。第一是参数可得性：定标系数由厂商在元数据里给出，但格式各异、单位各异，历史存档数据的系数可能经过多次修订，用错版本会让整幅影像系统性偏亮或偏暗。第二是大气参数的不确定性：气溶胶光学厚度是最大的误差源，实测站网稀疏，多数场景只能靠影像自身反演或气候态默认值，误差常在 10% 到 30%。第三是链条的完整性：反射率转换只是其中一环，太阳几何、地形、邻近效应、BRDF 各占一份，漏掉任何一环都会在特定场景下暴露。

本文按「物理量到实现」的顺序展开。先厘清辐射传输链条上的各个量，再讲 DN 到辐亮度、辐亮度到反射率的两步换算，然后覆盖大气散射吸收机制与 6S、MODTRAN 模型的实际配置，最后落到暗像元与 QUAC 这类无参数场景的经验方法、地形校正与质量评估。

本文假定读者已经理解影像波段与元数据结构，相关基础见 遥感影像基础与元数据 。校正后的反射率主要用于计算植被指数，方法见 光谱指数计算 ；校正质量直接影响时序分析结论，与 变化检测工程实践 中的阈值设定强相关。

目录

1. 辐射传输链条与物理量
2. DN 到辐亮度的定标
3. 表观反射率与太阳几何
4. 定标系数与产品级别
5. 大气散射与吸收机制
6. 6S 与 MODTRAN 模型实操
7. 暗像元与 QUAC 经验校正
8. 地形校正与邻近效应
9. 交叉验证与质量评估

1. 辐射传输链条与物理量

从太阳到传感器，能量经过的环节可以拆成一条清晰的链条。太阳辐射以辐照度 E（W/m²）到达大气顶，穿过大气时部分被散射与吸收，到达地表后按地物反射率 ρ 反射，反射能量再穿过大气上行到达传感器，被记录为辐亮度 L（W/(m²·sr·μm)）。

链条上的关键物理量：

两个容易混淆的概念：辐照度描述「来了多少能量」，辐亮度描述「往某个方向走了多少能量」。传感器测量的是后者，因为探测器有明确的瞬时视场。地表反射率是半球反射率，与观测方向相关，严格的表达需要 BRDF 描述，但多数产品近似为朗伯体。

完整的大气校正目标是把 ρ_TOA 反演为 ρ_surf，公式可概括为：

ρ_TOA = ρ_path + T_down · T_up · ρ_surf / (1 - s · ρ_surf)

其中 ρ_path 是大气程辐射，T_down 与 T_up 是下行与上行透过率，s 是球面反照率。这四个参数都由辐射传输模型给出，工程上不可能解析求解，必须借助查找表。

2. DN 到辐亮度的定标

第一步是把 DN 转成辐亮度，公式为线性关系：

L = gain × DN + offset

多数传感器的元数据直接给出增益与偏移。Landsat 8/9 OLI 在 MTL 文件里以 RADIANCE_MULT_BAND_x 与 RADIANCE_ADD_BAND_x 命名，单位是 W/(m²·sr·μm)/DN。Sentinel-2 的 L1C 产品则不同，元数据给的是 QUANTIFICATION_VALUE（约 10000）与 RADIANCE_ADD_BAND，需要用 DN 除以量化值再乘修正系数。

import numpy as np

def dn_to_radiance(dn, mult, add):
    dn = dn.astype(np.float32)
    return dn * mult + add

mult = 0.012246     # RADIANCE_MULT_BAND_4，Landsat 8 OLI Band 4 红光
add = -61.22947     # RADIANCE_ADD_BAND_4
rad = dn_to_radiance(band4, mult, add)

需要注意几个陷阱。其一，DN 的饱和值不能参与计算，Landsat 8 的 L1 产品饱和像元 DN 为 65535，直接换算会得到异常高的辐亮度，必须先掩膜。其二，某些产品已做过部分定标，元数据里标注为 reflectance 而非 radiance，重复换算会得到量级错误的结果。其三，探测器阵列的各像元响应不一致，严格的定标还包含相对辐射定标（去条带），这在高分辨率商业卫星上尤其重要。

gdalinfo -json LC08_L1TP_118038_20230512_MTL.json | jq '.["metadata"][]'
grep -i "RADIANCE_MULT" LC08_L1TP_118038_20230512_MTL.txt   # 查看定标系数

2.1 绝对定标与相对定标

绝对辐射定标解决「DN 对应多少物理辐亮度」，相对辐射定标解决「同一波段内各探元响应是否一致」。推扫式传感器的探元多达数千个，任何一个探元的响应偏离均值都会在影像上形成一条贯穿的条带。

条带检测与去除的常见做法：对整景计算各列（沿飞行方向）的均值曲线，用中值滤波提取平滑趋势，残差超过阈值即判定为条带。

import numpy as np
from scipy.signal import medfilt

def destripe(band, kernel=101, threshold=0.02):
    col_mean = band.mean(axis=0)
    trend = medfilt(col_mean, kernel)
    gain = trend / col_mean          # 每列增益
    gain = np.where(np.abs(gain - 1) > threshold, gain, 1.0)
    return band * gain[np.newaxis, :]

相对定标通常在传感器出厂时完成并内置于产品中，但探测器老化会导致响应漂移，长期存档数据需要按年度更新的相对定标系数重处理。这也是为什么同一景数据在不同时期下载，数值可能略有差异。

3. 表观反射率与太阳几何

辐亮度到表观反射率的换算引入了太阳几何与日地距离修正：

ρ_TOA = (π · L · d²) / (ESUN · cos θ_s)

d 是日地距离天文单位修正因子，ESUN 是波段内的太阳平均光谱辐照度（W/(m²·μm)），θ_s 是太阳天顶角。d² 项来自椭圆轨道：地球近日点（1 月初）时 d≈0.9833，远日点（7 月初）时 d≈1.0167，辐照度差异约 6.7%。忽略这一项会在南北半球同一季节的影像间引入系统性偏差。

import math

def earth_sun_distance(doy):
    return 1 - 0.01673 * math.cos(math.radians(0.9856 * (doy - 4)))

def toa_reflectance(rad, esun, sun_elev_deg, doy):
    d = earth_sun_distance(doy)
    cos_theta = math.sin(math.radians(sun_elev_deg))
    return math.pi * rad * d * d / (esun * cos_theta)

太阳高度角在元数据里常以 SUN_ELEVATION 给出，注意它是高度角不是天顶角，cos θ_s = sin(高度角)。这两个角度互换是新手最常犯的错误，在低太阳高度角时误差被放大。此外，影像范围内太阳高度角并不恒定，一景 185 公里的 Landsat 影像跨约 1.7 度纬度，太阳高度角差异可达 1 度以上，高精度处理需要逐像元计算。

各传感器的 ESUN 值不同，Landsat 8 OLI Band 4 约为 1895.33 W/(m²·μm)，Sentinel-2A MSI Band 4 约为 1665.48。这些常数随传感器版本变化，必须从官方文档获取而不能凭记忆。

3.1 热红外的亮温换算

热红外波段不能转反射率，而是转亮温。Landsat 8 的 TIRS 波段用普朗克反函数：

import math

K1 = 774.8853      # Band 10 定标常数 K1，单位 W/(m²·sr·μm)
K2 = 1321.0789     # Band 10 定标常数 K2，单位 K

def radiance_to_bt(rad):
    return K2 / math.log(K1 / rad + 1.0)

亮温是传感器视角的等效黑体温度，要得到真实地表温度还需做大气校正与比辐射率校正。比辐射率依赖地物类型，植被约 0.98、裸土约 0.95、水体约 0.99，用 NDVI 阈值法可以近似估计。忽略比辐射率校正会带来 1 到 3 K 的系统偏差。

注意 Landsat Collection 2 之后 Band 10 的定标常数经过修订，早期 Collection 1 的常数会产生约 0.3 K 的差异。跨版本拼接温度时序数据前必须统一到同一常数集。

4. 定标系数与产品级别

不同产品级别的处理程度差异极大，选错级别会让后续工作白做。

Landsat Collection 2 的 L1 产品提供 TOA 反射率与亮温，L2 提供地表反射率与地表温度。Sentinel-2 的 L1C 是 TOA 反射率，L2A 是地表反射率，欧洲航天局用 Sen2Cor 处理器生成。

产品级别的另一处差异是重采样与投影。L1C 常保留原始几何，L2A 已重采样到 10/20/60 米统一网格。混用不同级别的数据做时序分析，会因几何配准误差掩盖真实的辐射变化。

L2A_Process --resolution 10 S2A_MSIL1C_20230512T024551_N0509_R132_T51SUP_20230512T050213.SAFE   # 用 Sen2Cor 生成 L2A

gpt Sen2Cor.xml -Ssource=./S2A_L1C.SAFE   # 或调用 ESA 哨兵工具箱命令行

工程上要建立「产品级别台账」：每份入库数据记录级别、处理链版本、生成时间。同一分析任务内不允许混用级别，这是比算法选择更重要的纪律。

5. 大气散射与吸收机制

大气对辐射的影响分散射与吸收两类。散射不改变总能量但改变方向，吸收把能量转换成热。二者在不同波段的表现差异极大。

散射有三种机制：

• 瑞利散射：分子尺度散射，强度与波长四次方成反比，蓝光波段最强。它造成天空呈蓝色，是程辐射的主要来源，在可见光短波端可达表观反射率的 30% 以上。
• 米氏散射：气溶胶尺度散射，强度与波长关系较弱，受气溶胶类型与浓度控制。城市与沙尘区域差异巨大。
• 非选择性散射：云滴等大粒子散射，各波段近似相同，这也是云在影像上呈白色的原因。

吸收主要由水汽、臭氧、二氧化碳、甲烷造成，形成明确的吸收带。水汽在 940 纳米与 1130 纳米附近有强吸收带，臭氧在 550 纳米附近有查普伊斯吸收带，影响所有可见光波段约 1% 到 3%。

主要吸收带与传感器影响
O3    500~700 nm    可见光整体轻微衰减
H2O   940 nm        水汽反演通道，Sentinel-2 B9
H2O   1130 nm       植被水分胁迫
CO2   2010 nm       短波红外，Landsat B7 边缘

工程含义是波段选择：植被指数多选红光与近红外，因为这两个窗口的大气透过率高；水体指数用短波红外，因为水汽吸收增强了水陆差异。而大气校正的关键在于气溶胶光学厚度（AOD）的估计，它决定了可见光波段的程辐射量级。

5.1 气溶胶光学厚度的获取途径

AOD 是大气校正中最大的单一误差源，获取方式按可靠性排序：

import numpy as np

def dark_pixel_aod(red, blue, min_red_ref=0.01):
    # 暗像元法反演 AOD：找红光反射率接近零的像元
    dark_mask = red < min_red_ref
    if dark_mask.sum() < 100:
        return 0.2                       # 无足够暗像元，退回默认值
    return float(np.clip(blue[dark_mask].mean() * 5.0, 0.02, 1.0))

选择原则是「能用实测就用实测，能用卫星产品就用产品，实在没有才用默认值」。在流水线中应当把 AOD 来源与取值一并记录，作为后续质量评估的依据。同一景影像用不同 AOD 校正，红光波段的地表反射率差异可达 0.03 以上，足以改变植被指数的分类结论。

6. 6S 与 MODTRAN 模型实操

6S（Second Simulation of the Satellite Signal in the Solar Spectrum）是遥感领域最常用的辐射传输模型，覆盖 0.25 到 4 微米，考虑瑞利散射、气溶胶散射、多次散射与气体吸收，速度远快于 MODTRAN，适合批量生成查找表。MODTRAN 精度更高、波段更宽（可到热红外），但计算成本高、商业授权受限。

Py6S 是 6S 的 Python 封装，参数配置是核心工作：

from Py6S import SixS, AeroModel, AtmosProfile, Geometry

s = SixS()
s.geometry = Geometry.User()
s.geometry.solar_z = 35.0        # 太阳天顶角
s.geometry.solar_a = 140.0       # 太阳方位角
s.geometry.view_z = 0.0          # 观测天顶角（星下点）
s.geometry.view_a = 0.0

s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer)
s.aero_profile = AeroModel.Urban
s.aero_profile.aot550 = 0.25     # 550 nm 气溶胶光学厚度

s.wavelength = Wavelength(0.665, 0.680)   # 红光波段范围
s.run()
print(s.outputs.apparent_reflectance)
print(s.outputs.transmittance_total_scattering.down)

关键参数与取值建议：

批量处理的标准做法是预生成多维查找表，维度为太阳天顶角、观测天顶角、相对方位角、AOD，然后用插值快速反演。这样可以避免逐像元跑 6S，把单景处理时间从小时级压到分钟级。

import numpy as np

def interpolate_lut(lut, sza, vza, raa, aod):
    from scipy.interpolate import RegularGridInterpolator
    interp = RegularGridInterpolator(
        (sza_axis, vza_axis, raa_axis, aod_axis), lut, bounds_error=False,
        fill_value=None)
    return interp((sza, vza, raa, aod))

6.1 查找表的维度与规模设计

查找表是工程化的核心。维度太少插值误差大，维度太多内存爆炸。经验取值如下：

维度与步长建议
solar_z     0~70°   步长 5°   → 15 档
view_z      0~40°   步长 5°   →  9 档
rel_azimuth 0~180°  步长 15°  → 13 档
aot550      0~1.0   步长 0.1  → 11 档
合计 15 × 9 × 13 × 11 = 19305 条记录

近红外波段气溶胶影响弱，可以把 AOD 维度压到 3 档，内存立减。整表以 float32 存储约 1.5 MB，可常驻内存。生成时间取决于模型速度，6S 单条约 10 毫秒，全表约 3 分钟，一次性生成后可复用整个传感器生命周期。

插值要注意边界：太阳天顶角超过 70 度时 6S 的平面平行大气假设失效，球面大气修正不可忽略，此时应限制在 70 度内并对外插结果做标记。极区影像的太阳天顶角常年偏大，这类场景建议直接使用官方 L2A 产品。

7. 暗像元与 QUAC 经验校正

当缺乏气溶胶实测数据时，用影像自身估计大气参数的经验方法更实用。

暗像元法（DOS，Dark Object Subtraction）假设影像中存在反射率接近零的暗目标（深水体、浓密植被阴影），其表观反射率即等于程辐射。取影像各波段的最小值或 1% 分位数作为程辐射估计，从全图减去。

import numpy as np

def dos_correction(toa, percentile=1):
    dark = np.percentile(toa, percentile)
    return np.clip(toa - dark, 0, None)

DOS 的优点是零参数、快，缺点是假设过于简化：它忽略了透过率与下行辐照度的变化，在气溶胶浓度高的区域严重低估校正量。改进版 DOS1 加入透过率估计，DOS2 再考虑地表反射率与程辐射的关系。

QUAC（Quick Atmospheric Correction）是 ENVI 里的经验算法，不依赖元数据与实测参数。它利用影像内部不同地物的光谱统计关系，反演出一条增益与偏移曲线，把影像拉到「合理」的反射率范围。

envipy quac --input=toa.dat --output=quac.dat --sensor=unknown   # ENVI 命令行调用 QUAC

python -c "import spectral; print(spectral.__version__)"   # 用 spectral 库的 Python 实现

QUAC 的适用边界很明确：多波段、含多种地物、气溶胶空间分布均匀的场景效果好；单波段、地物单一、或气溶胶剧烈变化的场景不可靠。它输出的是相对反射率，跨影像不可比，适合单景分类而不适合时序分析。要在时序分析中用反射率，必须走物理模型校正或直接使用官方的 L2A 产品。

8. 地形校正与邻近效应

山区影像的辐射畸变来自地形起伏：向阳坡接收的辐照度远大于背阴坡，同一种植被在两侧坡面呈现完全不同的亮度，直接做分类会把阴影误判为水体或裸土。

地形校正的经典方法是 C 校正与 Minnaert 校正。C 校正假设地表朗伯，用太阳入射角余弦做归一化：

ρ_corrected = ρ_observed · (cos θ_s + c) / (cos i + c)

其中 i 是太阳入射角与坡面法线的夹角，c 是经验常数（通常取坡面余弦的均值）。Minnaert 校正引入非朗伯指数 k，更适合粗糙表面。两者都需要坡度与坡向数据，由 DEM 计算。

gdaldem slope dem.tif slope.tif -compute_edges -alg Horn
gdaldem aspect dem.tif aspect.tif -compute_edges
gdaldem hillshade dem.tif hillshade.tif -z 1.5 -az 315 -alt 45

邻近效应是另一类容易被忽视的畸变：目标像元接收到的辐射中包含周围地物散射进来的能量，在雾霾或高气溶胶条件下，邻近像元贡献可达 10% 到 20%。6S 提供邻近效应参数，但多数业务处理忽略它，代价是在水体与陆地边界、云边缘出现模糊的光晕。

地形校正的副作用同样需要注意：过校正会在阴影区产生异常高值，尤其当 c 取值不当时。稳妥做法是对校正前后的影像做同一地物的光谱统计，确认阴影与向阳坡的同类地物光谱趋于一致再放大规模处理。校正后的影像更适合做 云检测与云掩膜 ，因为阴影误判会显著减少。

8.1 BRDF 与多角度效应

地表并非朗伯体，反射率随太阳入射角与观测角变化，这种方向性用 BRDF 描述。同一地物在不同观测几何下反射率差异可达 20% 到 40%，这是宽幅传感器（如 MODIS 的 2330 公里幅宽）边缘与中心亮度不一致的根本原因。

业务上常用核驱动模型把 BRDF 归一化到固定几何：

ρ(θ_s, θ_v, φ) = f_iso + f_geo · K_geo + f_vol · K_vol

三项分别对应各向同性散射、几何光学散射（阴影效应）与体散射（冠层内部多次散射），f 系数由多角度观测拟合得到。MODIS 的 MCD43 产品即基于此模型，提供天底反射率（NBAR）。

对单角度传感器（Landsat、Sentinel-2），无法直接拟合 BRDF，只能用 C 因子法或固定核系数近似。实际工程中更常见的做法是接受这一误差，但在时序分析时用同一季节、相近太阳角的数据，避免 BRDF 效应被误读为地物变化。这也是为什么严格的时序分析要求数据在同一天文条件下采集。

9. 交叉验证与质量评估

辐射定标与大气校正没有绝对真值，只能通过交叉验证控制误差。常用三条路线：

第一条是与官方 L2A 产品对比。取同一景数据的自处理结果与 Sentinel-2 L2A 或 Landsat L2 逐像元比较，计算偏差（bias）、均方根误差（RMSE）与相关系数。城市与裸地等均质区域偏差应控制在 0.01 到 0.02 以内，植被区因 BRDF 效应可放宽到 0.03。

第二条是与地面实测对比。用 ASD 光谱仪等设备在过境时刻同步测量地表反射率，比较同点位像元值。这是最可靠的方法，但受云、时间同步、尺度不匹配限制，实测点数量通常很少。

第三条是内部一致性检查。同一传感器相邻日期的影像，在稳定地物（沙漠、深水体）上的反射率应接近；同一景内水体的近红外反射率应接近零，如果显著大于零说明校正不足或存在云污染。

import numpy as np

def evaluate(est, ref, mask=None):
    if mask is not None:
        est, ref = est[mask], ref[mask]
    diff = est - ref
    return {
        "bias": float(np.mean(diff)),
        "rmse": float(np.sqrt(np.mean(diff ** 2))),
        "r": float(np.corrcoef(est, ref)[0, 1]),
        "n": int(est.size),
    }

质量评估还要关注波段间关系。植被在红光与近红外的关系、水体在近红外与短波红外的关系，都是物理约束。如果校正后这些关系被破坏，说明某一步引入了错误。

权衡取舍

• 物理模型 vs 经验方法：6S 精度高但需气溶胶参数与计算成本，QUAC 零参数但只有相对反射率，按是否有实测数据决定。
• AOD 来源：AERONET 站网精度最高但空间稀疏，影像反演覆盖广但有反演误差，气候态默认值最省事但误差最大。
• 官方 L2A vs 自处理：官方产品一致性好、可复现，自处理灵活可加地形校正，建议默认用官方、特殊需求才自处理。
• 地形校正：山区必须做，但过校正会制造假高值，c 参数需按区域标定。
• 邻近效应：雾霾场景收益明显，晴朗场景可忽略，按业务精度要求决定是否实现。
• 查找表 vs 逐像元：查找表快百倍但需插值，逐像元精度略高但工程上不划算。
• 时间成本：全链路物理校正单景分钟级，经验法秒级，按数据量级选择。

常见坑清单

• 太阳高度角当天顶角：现象是反射率整体偏高数倍，原因是用错角度定义，规避方法是确认 cos θ = sin(高度角)。
• 忘记日地距离修正：现象是冬夏影像亮度系统性差异，原因是忽略 d² 项，规避方法是按儒略日计算修正因子。
• 定标系数版本错：现象是整景影像偏亮或偏暗，原因是用了他批次或他传感器的系数，规避方法是从元数据文件直接读取而非硬编码。
• 饱和像元未掩膜：现象是云顶或高反射目标出现异常高值，原因是 DN 饱和值参与换算，规避方法是先掩膜再计算。
• 对 L2A 重复校正：现象是反射率被压得过低，原因是产品已含地表反射率，规避方法是检查产品级别台账。
• QUAC 用于时序分析：现象是不同日期影像反射率不可比，原因是 QUAC 输出相对反射率，规避方法是时序分析只用物理模型或官方 L2A。
• 波段范围与 ESUN 不匹配：现象是反射率量级偏差，原因是 ESUN 对应波段不同，规避方法是按传感器版本查官方文档。
• 忽略气溶胶类型：现象是城市与清洁区校正不一致，原因是统一用了大陆型气溶胶，规避方法是按区域调整气溶胶模型。
• 地形校正参数未标定：现象是阴影区出现异常高反射率，原因是 c 取值过大，规避方法是用同区域样本做标定。
• 未做质量评估：现象是错误在时序分析阶段才暴露，原因是跳过交叉验证，规避方法是把 bias 与 RMSE 纳入流水线门禁。

小结

辐射定标与大气校正的本质是把传感器记录的相对值还原成有物理含义的地表反射率。链路分两步：DN 到辐亮度依赖厂商系数，辐亮度到反射率依赖太阳几何与大气参数。前者的风险在参数版本，后者的风险在气溶胶估计。

方法选择取决于数据可得性。有实测 AOD 与计算资源就走 6S 或 MODTRAN，没有就走 DOS 或 QUAC 这类经验方法，但必须清楚经验方法输出的相对反射率不适合时序分析。地形校正与邻近效应是山区与雾霾场景的必备补充，也是最能体现工程细节的两环。

下一步建议把校正结果接入 光谱指数计算 ，观察 NDVI 在时序上的稳定性作为校正质量的间接指标；也可以结合 GDAL 栅格数据读写 把整条校正链封装成可并行、可复现的处理管道。