ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

遥感图像预处理全链路:辐射定标、大气校正与几何校正实操指南

遥感图像预处理全链路:辐射定标、大气校正与几何校正实操指南

第一次拿到的原始影像,看着光谱数不少、分辨率也不低,可一打开就发现整幅图又暗又糊:地物灰蒙蒙一片,纵向还带着几道莫名其妙的条纹,城区边缘跟底图对不上,山区的山脊上阴阳脸严重到看不出地物边界。做遥感图像预处理,本质就是把这一景"能看但不好用"的原始数据,变成"能被模型和公式信任"的定量产品,中间要过的几道工序——辐射定标、大气校正、几何校正与正射纠正、降噪裁剪镶嵌、重采样——顺序不能乱、参数不能糊弄。这篇内容写给两类人看:刚接触遥感影像处理、还在纠结"到底要不要做大气校正"的新手,以及已经能把流程跑通、但结果总觉得哪里不对劲、需要一套可复现参数与核查清单的实操者。下面我按一条完整处理链路把每一步拆开,讲清楚每一步为什么要做、卡在哪个参数上、踩过哪些坑。

1. 先把预处理的范围划清楚:哪几步必须做,哪几步是过度加工

1.1 从原始观测值到可比较物理量的完整链路

影像是传感器在大气层顶接收到的辐射能量记录,存储时为了让数据量小、便于传输,通常会被量化成整型数值,也就是我们常说的观测值。这个数值本身没有物理单位,不同时间、不同传感器拍出来的同一块地,数值完全没法横向比较。所以预处理的第一个动作一定是辐射定标,把无量纲的观测值换算成大气层顶辐亮度或者大气层顶反射率;第二个动作是大气校正,剥离大气散射与吸收的影响,落到地表反射率;第三个动作是几何校正与正射纠正,把像元摆到正确的地理位置上;第四个动作是降噪、去云、裁剪、镶嵌、重采样这些"整形"操作。

这条链路的顺序有硬约束。辐射定标必须在大气校正之前,因为大气校正模型吃的是有物理单位的辐亮度;几何校正理论上和辐射过程独立,但如果你要用暗像元法估气溶胶,最好先把影像裁到研究区,避免大片水体拉偏统计;重采样一定要放在几何纠正那一步一次性做完,不要在裁剪时重采样一次、投影时又重采样一次,两次插值会明显钝化边缘,尤其是建筑和道路这类高频地物。

我在实际项目里见过不少新手把顺序搞反:先裁剪、再重投影、再定标、再大气校正,一路下来重采样了三四次,最后拿去算植被指数的结果比参考值低了一截,查了半天以为是大气校正参数设错,其实是插值把窄波段的光谱特征磨平了。这个顺序问题排查起来特别费时间,所以一开始就把链路定死,比后面调参划算得多。

1.2 按应用目标决定步骤取舍,避免无意义的加工

不是所有应用都要走全流程。做目视解译、做变化检测的初步筛查、做土地利用分类的样本采集,其实做到正射纠正加适度增强就够了,上大气校正反而可能因为参数不准引入新的系统性偏差。但只要你涉及跨时相比对、定量反演、光谱匹配、指数阈值标定,大气校正就是绕不过去的一环——同一块农田在两个季节的观测值差异,可能一半来自物候、一半来自大气路径长度不同,不做校正就等于把两种信号混在一起。

判断标准可以简化成一句话:下游环节是否直接比较数值。只要答案是"是",就必须走到地表反射率;答案是"只看纹理和形状",做到正射纠正即可。我曾经为了省时间跳过一次大气校正,直接拿大气层顶反射率算植被指数做长时序趋势,结果雨季和旱季的曲线呈现出明显的相反趋势,后来补做校正才恢复正常——这个教训让我后来在所有涉及多时相的项目里,无条件先做校正再分析。

还有一个容易被忽略的取舍:空间分辨率是否需要重采样统一。多源数据融合时,很多人习惯把所有数据重采样到最高分辨率,觉得"越精细越好"。实际恰恰相反,把 30 米的数据三次卷积到 10 米,不会凭空产生细节,只会让相邻像元高度相关,后续做纹理特征提取时会出现大量伪高频信息,分类精度反而下降。合理的做法是把高分辨率数据降采样到较低分辨率,或者统一到一个中间尺度。

2. 辐射定标:把无量纲数值翻译成有物理意义的辐亮度

2.1 定标公式与参数来源,一次算给你看

定标的核心就是一个线性关系:辐亮度等于增益乘以观测值再加上偏移量。增益和偏移量这两个数不是我们算出来的,而是随数据一起下发的元数据文件里给的,通常在影像自带的描述文件、头文件或者是产品说明里,一个波段一组。有些产品直接给辐亮度增益和偏移,有些给的是最大最小辐亮度再加上量化级数,那就要自己按比例换算。

举个具体例子。假设某波段量化是 12 位,也就是观测值取值范围 0 到 4095,元数据给出该波段最大辐亮度为 620,最小辐亮度为 -6.2。那么增益等于(620 减负 6.2)除以(4095 减 0),约等于 0.1528;偏移量就是最小辐亮度,即 -6.2。这样一个观测值为 2000 的像元,定标后的辐亮度大约是 299.4。数值算出来之后一定要跟元数据里给的辐亮度范围核对一下,如果算出几百上千这种明显超出量级的数,八成是增益的单位搞错了——有的产品增益按"每千观测值"给,有的是"每观测值",差了一千倍。

如果要进一步换算成大气层顶反射率,还要引入太阳天顶角和日地距离因子。公式是反射率等于圆周率乘以辐亮度乘以日地距离平方,再除以波段大气外太阳辐照度与太阳天顶角余弦的乘积。日地距离因子随日期变化,在一月初最小、七月初最大,变化幅度约 3.4%。这个修正听起来很小,但在做长时序趋势分析时,3% 的系统性差异足以淹没真实的年际变化信号,所以凡涉及跨月比较,这个因子必须带上。

注意:定标系数一定要按波段分别读取,不要用一个波段的系数套到所有波段上。我见过最典型的一次事故,就是把全色波段的增益套用到多光谱波段,结果近红外波段辐亮度整体偏高了将近一个数量级,后续算出来的水体指数全是负值。

2.2 定标环节最容易翻车的四个细节

第一个是整型溢出。原始观测值是整型存储,如果你在代码里直接做乘法并且输出仍然声明为整型,超过范围的部分会被截断。正确做法是先把数据转成 32 位浮点再做运算,输出也声明为浮点格式。这个坑在批处理里特别隐蔽,因为单景看不太出来,只有统计直方图时才会发现高值端被整齐地压平了。

第二个是无效值与填充值。原始产品里通常用 0 或者一个特殊值表示无数据区域、卫星边缘、云掩膜外区域。这些值如果参与定标运算,会变成负的辐亮度,甚至污染后续的大气校正统计。处理前必须把填充值单独识别出来,写成掩膜层,参与统计时排除掉。我习惯在定标脚本里直接把 0 和元数据声明的填充值替换成 NaN,后面所有统计一律用忽略 NaN 的函数。

第三个是波段间空间配准。有些传感器的多光谱波段不是严格同时同视场成像的,波段之间有几十分之一像元的偏移,这在计算归一化差值指数时会表现为边缘出现彩色描边。如果发现指数图的地物边界有红蓝色边,先检查是不是波段配准问题,不要急着怀疑大气校正。

第四个是元数据版本不一致。同一个传感器不同时期下发的产品,定标系数可能更新过,尤其是早期在轨阶段的系数往往在后来的产品说明里被修订。批处理时如果统一用一个写死的系数表,跨年份数据就会出现台阶状跳变。

3. 大气校正:最难、最贵、也最容易被做过头的一步

3.1 什么情况下必须做,什么情况下可以不做

大气对观测的影响分两类:散射和吸收。散射让短波段整体变亮、对比度下降,蓝色波段受影响最重;吸收则主要是水汽和若干气体在特定波段打出的"坑"。做大气校正的目标,就是把这两部分影响从观测信号里剔除掉,还原地表的反射特性。

必须做的场景:跨时相定量比较、植被生理参数反演、水体悬浮物与叶绿素估算、地表温度反演的配套处理、矿物光谱匹配。这些应用对绝对数值的准确性要求高,尤其是水体,信号本身就弱,大气散射的贡献在水体上空可能占观测信号的八成以上,不做校正等于在雾里看水。

可以不做或者缓做的场景:单纯的目视解译、训练样本的目视勾画、以纹理和形状为主的目标检测。这些场景里,大气散射相当于给整幅图加了一层平滑的背景,对边界和形状的影响有限。

还有一个中间选项:只做相对大气校正。也就是不做物理模型的绝对反演,而是用暗像元法或者直方图匹配,把多时相影像的相对关系调一致。这种做法在变化检测、大区域镶嵌时非常实用,成本低、稳定性好,虽然拿不到真实的地表反射率量级,但时序可比性够用。

3.2 常用方法横向对比与选型逻辑

方法类别代表做法需要的辅助数据精度水平适用场景主要风险
辐射传输模型逐像元调用大气辐射传输计算观测几何、气溶胶光学厚度、水汽、臭氧、DEM、光谱响应函数高,可反演真实反射率定量反演、多时相对比参数不准时误差被放大,速度慢
商业软件内置模块集成式大气校正模块观测几何、气溶胶模式、能见度或光学厚度较高常规项目、中等数据量黑箱参数,出错难定位
暗像元法假定影像中存在反射率极低的暗目标仅影像自身中低,仅相对校正变化检测、镶嵌匀色缺暗目标时假设失效
快速大气校正基于影像统计的经验模型仅影像自身低到中快速出图、预研波段间关系不稳定时偏差大
地形校正附加项结合坡向坡度修正辐照DEM提高山区精度山区、丘陵过校正导致阴坡虚假变亮

选型逻辑我一般按三条判断:下游是否要做绝对量反演,决定要不要上辐射传输模型;研究区有没有大面积暗目标(深水、浓密植被阴影),决定暗像元法是否可用;数据量有多大、算力是否够,决定要不要退而求其次。

一个特别实际的经验:别一上来就全套参数硬上。先用快速方法跑一小块,看看校正前后的地物光谱曲线是否平滑——健康的植被光谱在红边附近应该有明显的抬升趋势,如果校正后曲线出现锯齿状的异常跳跃,说明某个波段的参数出了问题,这时候再回去查水汽波段和波段响应函数。

3.3 参数准备清单与实操中的坑

辐射传输模型法需要准备的参数比很多人想象的要多,而且每一项都会影响结果:

  • 观测几何:太阳天顶角、太阳方位角、观测天顶角、观测方位角。前两个通常元数据里有,后两个在近天底观测的产品里接近 0,但侧摆角大的产品必须读实际值。用错观测天顶角,路径长度估计就会偏,短波段误差尤其大。
  • 气溶胶光学厚度:这是最关键也最难获取的参数,可以用影像自身反演(暗像元法、对比度法),也可以用地基观测或再分析资料。近些年常用的做法是结合地表反射率先验做反演,精度比固定光学厚度提升明显。
  • 水汽含量:对近红外和短波红外波段影响显著,可以用影像自身的水汽吸收波段反演,也可以从再分析资料读取。
  • 臭氧含量:影响小但影响是系统性的,通常取区域月均值即可。
  • DEM:用于计算海拔带来的气压修正,山区尤其重要。
  • 光谱响应函数:这是最容易被忽略的一项。辐射传输模型计算的是单色辐亮度,而传感器实际接收的是波段范围内的积分。不做响应函数加权,波段越宽误差越大。

注意:大气校正输出一定要检查有没有负值。真实的暗像元,比如深水或者浓密植被阴影,反射率可以接近零但不会为负。如果出现大片负值,通常意味着气溶胶光学厚度给低了或者观测几何读错了。处理办法不是简单截断成零,而是回头检查参数。

我在山区做过一景校正,第一次结果整片阴坡都是负值,阳坡却还算正常。查下来是 DEM 用了粗分辨率产品,山谷和山脊高程差被抹平,路径长度估算失真。换成较精细的 DEM 之后负值基本消失。这也是为什么我一直强调:地形起伏大的区域,DEM 的质量比模型选型更重要。

4. 几何校正与正射纠正:让像素落到真正的地面上

4.1 粗校正与精校正的差别,别把两者混为一谈

粗校正用的是卫星轨道参数和姿态数据,通过有理多项式系数这类模型把影像从传感器坐标系映射到地面坐标系。这个过程不需要人工干预,速度很快,对于平坦地区,精度通常能到几个像元的水平,日常做区域概览、做在线地图底图够用。

精校正是在粗校正基础上,用地面控制点去修正残余误差。控制点可以从已有正射影像、高精度地形图、控制点库中获取,一般要均匀分布在整幅影像上,边缘和四角都要有。控制点数量没有绝对标准,但至少要保证每个象限有足够的点,且点的残差要落在合理范围内——我的经验是整体均方根误差控制在三分之一到二分之一个像元以内,超过一个像元就要考虑剔除粗差点重新选。

地平线以上还有一层必须做的事:正射纠正。粗校正只补偿了传感器位置和姿态带来的误差,没有考虑地形起伏导致的像元位移。在山区,同样的坡度,侧视观测和天底观测的位移方向和量级完全不同,一座 500 米高的山在较大侧摆角下可以造成十几个像元的位移,这个量级不做纠正是没法做小区域分析的。

4.2 控制点选取与重采样方法的选择

控制点选在什么位置,直接决定校正质量。我的原则是:选在稳定、清晰、小尺寸的目标上。道路交叉口、水库堤坝拐角、桥梁端点、大型建筑的直角角点都比较理想;尽量避免选在植被边界、水陆交界、施工区域,这些地方本身就年年变,控制点坐标今天准确明天就不准了。

重采样方法的选择有明确约束:

场景推荐方法原因
连续型光谱数据(辐射定标后、反射率产品)三次卷积或双线性保留光谱平滑性,避免台阶感
分类结果、土地利用图、掩膜层最近邻不产生不存在的类别值
快速预览、粗略检查最近邻速度快,便于批处理
需要保留边缘清晰度的指数产品双线性兼顾平滑与边缘

这里有一个反直觉的点:分类结果绝对不能双线性重采样。双线性会在类别边界插出 0.5 这类不存在的值,如果你把分类码当成整型存储,那些插出来的值会被截断成某个真实类别,图上会出现莫名其妙的细碎斑块,肉眼很难判断是分类错误还是重采样错误。我见过有人因为这个原因反复调分类器参数,折腾了两天。

4.3 投影选择与像元对齐的细节

投影不要随便选。研究区跨带的时候,要么按中心经线自定义投影,要么统一到一个覆盖范围更大的投影里,否则拼接处会出现明显的错位。像元对齐也很关键:如果两景影像的网格原点差半个像元,镶嵌之后的接缝处会出现一列模糊的重叠像元,做时序比对时会形成一条虚假的变化带。GDAL 的 warp 工具里有对齐选项,批处理时把目标网格的参考影像固定下来,让所有景都对齐到同一网格,能省掉后面一大堆麻烦。

5. 降噪、去云、裁剪、镶嵌:让数据能直接进模型

5.1 条带与坏像元的处理思路

推扫式成像的传感器在轨时间长了,某些探测元会老化,表现为单列或者单行像元的响应异常,图像上呈现为纵向或横向的细条纹。去条带的核心思路是用邻列统计量重新标定异常列:对正常列计算均值和方差,再让异常列匹配这个统计特性。这个做法对轻微条带效果好,但如果连续多列异常,或者条带宽度不均匀,强行匹配会带来新的伪纹理。

坏像元(也常叫死像元、热噪声像元)的修复要用邻域插值,但插值范围不要太大。我的经验是只对孤立坏像元用 3×3 或 5×5 邻域的中值替代,超过一定面积的连续缺失区域干脆标记成无数据,交给后续的时序填补或者干脆剔除,不要试图用插值"造"出一片区域。

5.2 云检测与云掩膜的重要性被严重低估

云和云影在分类和反演里是最大的污染源。检测手段主要有三条路:用产品自带的质量波段、用专门的云检测算法、用简单的阈值组合。质量波段最省事,但不同版本的产品编码规则会变,用之前一定要看清楚每一位代表什么。阈值法要结合波段特性设计,比如用蓝色波段亮度判断厚云、用短波红外与近红外比值判断薄云、用云影在近红外波段的低值特征配合太阳几何推算影子位置。

掩膜做完不要急着删除,把掩膜层单独存一份,后面做统计的时候用得上。我有个项目因为没保存掩膜,后期复核时发现某片区域的异常高值其实是云边缘的半透明薄云,只能重新跑一遍整个流程。

5.3 裁剪、镶嵌与接缝匀色

镶嵌最容易出问题的地方是接缝。相邻两景如果成像时间差几个月,整体色调差异会非常明显,直接拼在一起像贴了两张不同曝光的照片。处理办法有三层:先做辐射归一化,把各景的统计分布调到一致;再做重叠区域的加权过渡(羽化),让接缝处平滑过渡;最后如果还有残留的色差,在重叠区做直方图匹配。

提示:匀色不要做过头。把每景都强行匹配到同一参考分布,会抹掉真实的季节差异。做长时序镶嵌时,宁可保留一定的色调差异,也不要为了"好看"而破坏数值的真实性——匀色改动的是数值,不是只有外观。

6. 把流程工程化:批处理骨架与质量核查清单

6.1 目录结构与元数据管理

数据量一大,靠记忆管理文件必然出错。我习惯按"原始—中间—成果"三级目录组织,中间产物按处理步骤分文件夹,文件名里带上关键参数标识(比如定标版本、大气校正方法代号、投影代号)。这样出了问题能快速定位是哪一步引入的偏差,也方便只重跑某一步而不用从头来。

元数据必须和成果一起存。每景处理完,把输入的定标系数、大气参数、控制点文件、软件版本、处理时间写成一个伴随文件。半年后回来复核,最怕的就是"这景数据当时用什么参数做的"完全想不起来。

6.2 一个可直接套用的批处理骨架

下面是我常用的 Python 骨架,用 GDAL 做读写,把定标和大区域拼接串起来,思路可以直接照搬。

import os import numpy as np from osgeo import gdal gdal.UseExceptions() def read_band(path, band_idx): ds = gdal.Open(path, gdal.GA_ReadOnly) band = ds.GetRasterBand(band_idx) arr = band.ReadAsArray().astype(np.float32) nodata = band.GetNoDataValue() if nodata is not None: arr[arr == nodata] = np.nan return arr, ds.GetGeoTransform(), ds.GetProjection() def radiometric_calibrate(dn, gain, offset): # 先转浮点再运算,避免整型截断 return dn * np.float32(gain) + np.float32(offset) def to_reflectance(lum, esun, sun_zenith_deg, doy): # 日地距离因子,简化近似公式 d = 1.0 - 0.0167 * np.cos(2 * np.pi * (doy - 3) / 365.25) cos_zen = np.cos(np.deg2rad(sun_zenith_deg)) return (np.pi * lum * d * d) / (esun * cos_zen) def clip_by_mask(src, mask_shp, dst): gdal.Warp(dst, src, cutlineDSName=mask_shp, cropToCutline=True, dstNodata=np.nan) def mosaic_and_resample(vrt_path, out_path, res, epsg): gdal.Warp(out_path, vrt_path, dstSRS=f"EPSG:{epsg}", xRes=res, yRes=res, resampleAlg="bilinear", # 分类产品改 near srcNodata=np.nan, dstNodata=np.nan)

这段代码里刻意保留了两个容易被省掉的操作:一是先转浮点再乘增益,二是重采样算法作为参数暴露出来。批处理脚本写成"算法写死、参数可配"的形式,比把参数全部写死要安全得多——你永远不知道下一批数据的分辨率和投影是什么。

6.3 处理完的成品必须过一遍核查

我把核查清单固定成下面这几项,每次交付前逐条过:

核查项检查方法合格标准
数值范围统计最小最大值与分位数反射率大体落在 0 到 1,少量高值正常,无大面积负值
无效值覆盖率统计 NaN 占比与研究区形状吻合,无整条缺失
波段间一致性抽查典型地物光谱曲线曲线形态合理,无锯齿跳跃
几何精度与参考底图叠加查同名点偏移小于 1 像元
辐射一致性重叠区直方图对比两景统计分布接近
坐标系与网格查看元数据与项目约定一致,无半像元错位

这份清单看起来啰嗦,但每一条我都真实遇到过问题。尤其是"波段间一致性"这一项,很多错误在单波段显示时完全看不出来,只有把典型地物的光谱曲线画出来才发现某个波段整体偏移。

7. 常见问题速查与踩坑记录

7.1 典型故障对照表

现象最可能原因排查方向处理办法
整幅图偏暗、地物灰蒙未做大气校正或校正过轻对比深水区反射率是否为负提高气溶胶光学厚度或换方法
短波波段出现大片负值校正过重或观测几何读错核对天顶角数值与元数据修正参数,重新计算
纵向条纹明显探测元响应差异逐列统计均值方差列方向统计匹配
指数图地物边界有彩边波段间配准偏差检查波段间偏移量先做波段配准再计算
与底图整体偏移只做粗校正未做正射查看地形起伏与侧摆角引入 DEM 与控制点
镶嵌接缝明显多景辐射差异未归一对比重叠区直方图辐射归一加羽化
分类图出现碎斑分类结果被双线性重采样检查重采样方法改最近邻重跑
时序曲线出现台阶定标系数版本不一致查各景元数据系数来源统一系数表重算

7.2 几条只有踩过才知道的经验

先跑小区域再跑全域。不要等整景跑完才发现参数有问题。我现在的习惯是固定裁一块包含水体、植被、建筑、裸土的小区域做试跑,四类地物的光谱曲线都合理了,再铺开跑全图。这个小区域其实就是我的"试纸",几年下来没换过。

参数改动一次只改一个。大气校正的失控感很强,很多人一看到结果不对就同时改气溶胶、改水汽、改观测几何,改完发现结果变好了,但不知道是哪一项起了作用。下次遇到同类数据还是不会调。一次只动一个参数,记录每次的结果,两三景之后就有一套属于自己的经验值了。

地形校正不要滥用。地形校正对阴坡的提升效果很直观,很容易让人上瘾。但过校正会让阴坡变得比真实值更亮,导致同一类植被在阳坡和阴坡的光谱差异被人为反转,分类时反而更容易错。判断是否过校正,可以看同一片均匀林地在不同坡向上的反射率是否被拉平到接近水平——如果阴坡比阳坡还亮,那就是校正过头了。

保留一份不做任何处理的原始副本。这条听起来废话,但我确实遇到过中间产品被覆盖、想回到原始状态重跑的情况。原始数据放在只读位置,任何处理都在副本上进行,这个习惯能救命。

记录处理时间与数据量的关系。当数据量从几十景涨到上千景,很多原本忽略的问题会被放大:内存不够导致分块读取逻辑出错、临时文件占满磁盘、并行处理时读写冲突。项目一开始就把分块尺寸、并行度、临时目录规划好,比后期返工省力得多。

7.3 关于精度与效率的取舍

最后想说一个技术之外但很实际的事:预处理的精度不是越高越好,而是要和下游需求匹配。做区域尺度的物候分析,用快速大气校正就够了;做小流域水质监测,那就要老老实实上辐射传输模型,还得把气溶胶参数反演做扎实。我见过花两周时间把校正精度从 5% 压到 2%,结果下游分类精度只提升了零点几个百分点的情况,投入产出完全不划算。

反过来,某些环节省不得。几何精度就是典型:它的误差会传导到每一个像元的空间位置上,一旦研究区的田块尺度接近像元尺度,半个像元的偏移就足以把一个田块的混合像元比例彻底改变,后面再怎么优化模型也补不回来。这类"基础误差"环节,宁可慢一点,也要做扎实。

我自己现在的做法是给每个项目定一份"预处理等级":概览级只做粗校正加增强,分析级做定标加正射,产品级再加完整大气校正和精度验证。等级定下来,流程就不会随意摇摆,交付质量也稳定。这套做法不是从哪本书上抄的,是被返工逼出来的——每次返工几乎都能追溯到某一步被"顺手简化"了。

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进