★ FEATURED ARTICLE 前阵子一个朋友从数据平台下载了两景高分影像打开软件一看山体、河流、村庄清清楚楚他说这图这么清晰直接拿来分类应该没问题。我劝他先做预处理他不听结果出了分类结果以后同一片裸地和耕地反复互串山体阴影被识别成了水体他回头来问我到底是分类器的问题还是数据的问题。我说多半是你跳过了高分遥感数据预处理这一步。所谓高分遥感数据预处理指的是在正式分析之前对来自卫星或航空平台的高分辨率遥感影像做辐射定标、大气校正、几何校正、正射校正等系列标准化操作把传感器记录的数字量化值转换成能够真实反映地表辐射特征、空间位置准确可靠的影像产品。它解决的核心问题是原始数据“看着能用”和“数据真的能用”之间的巨大差距。这篇文章适合刚接触遥感的学生、地理信息从业者也适合所有准备拿高分影像做分类、反演、变化检测的人。我会把流程拆开讲清楚并且把我这几年实际操作中踩过的坑一并写出来。1. 为什么高分数据不能跳过预处理直接使用1.1 传感器记录的是灰度值不是地面真实信息很多人第一次接触高分影像会有一个错觉影像存储的都是像元值看起来是“颜色”于是认为这些像元值可以直接反映地面属性。实际上高分影像原始数据里存的是一串DN值Digital Number也就是传感器接收到的辐射能量经过量化以后的离散数值。对多数光学高分卫星来说DN值通常以16位整型存储像元取值范围是0到65535。这个DN值本身没有物理含义它受多种因素叠加影响大气分子和气溶胶对太阳辐射的吸收与散射、太阳高度角随纬度和季节的变化、传感器自身增益和偏移设置、地形起伏引起的阴阳坡差异等等。换句话说同一块地物在晴天和薄云天气下拍出来的DN值不同在上午十点和下午两点拍的也不同传感器侧摆大一点小一点更是会变。如果拿着这样的DN值直接去做分类、做植被指数、做定量反演相当于把大气条件、观测几何都当成地物信息一起算了进去。我在实际项目里见过一个很典型的例子同一研究区、相隔半个月的两景影像NDVI直方图整体偏移了0.1以上导致变化检测结果出现大面积伪变化。后来重新做了辐射定标和大气校正两景影像的NDVI分布才基本对齐。这就是预处理的必要性最直观的证明。1.2 跳过预处理的连锁反应分类、变化检测和参数反演全遭殃如果只是做目视解译原始DN值或许还能用毕竟人眼对相对亮度和纹理比较敏感。但一旦进入定量分析问题就接二连三地出现。先拿植被指数来说。NDVI的定义是近红外波段反射率减去红光波段反射率再除以两者之和。大气对红光和近红外的影响程度不一样红光波段受气溶胶散射影响更强所以不做大气校正时NDVI通常会被低估植被区域看起来“不绿”。反之如果传感器增益设置得不合理DN值饱和又会让NDVI在某些高亮区域出现虚假的极值。再说分类。监督分类算法本质上是根据训练样本的像元值分布来划分特征空间。原始DN值没有经过定标不同时相、不同景的影像统计特征不一致训练样本在一景影像上的特征到另一景影像上就漂移了。我实际测试过用同一套训练样本去分类两景同地区不同日期的未预处理影像耕地与裸地的分类边会整体偏移好几十米这根本不是分类器参数能调回来的问题。变化检测更是重灾区。不同时相的影像太阳高度角不同、大气条件不同不做辐射归一化直接做差值或比值分析得到的变化区域里混入了大量“大气变化”和“观测角度变化”的伪信号。所以行业里有个共识预处理花的时间不是浪费它决定了后续分析的下限。1.3 不同高分载荷的预处理起点不一样高分遥感是个大筐常见的载荷类型包括全色影像、多光谱影像、高光谱影像、立体测绘像对甚至还有SAR雷达影像。它们的预处理侧重点差异非常大。我先用一张表说明后面再重点展开光学高分影像的流程。数据类型典型问题预处理重点全色影像几何畸变、地形投影差几何校正、正射校正多光谱影像辐射失真、大气散射辐射定标、大气校正、几何校正高光谱影像大气吸收影响严重、条带噪声大气校正、坏线修复、光谱重建立体像对高程定位偏差区域网平差、DEM生成SAR影像斑点噪声、地形辐射畸变多视滤波、地形辐射校正这篇文章主要讨论光学高分影像的预处理这也是绝大多数土地分类、目标识别、植被监测项目的起点。SAR有另一套逻辑以后有机会单独写。2. 高分数据预处理全流程拆解定标、校正与融合的先后顺序2.1 辐射定标从DN值到反射率的物理换算辐射定标是整个预处理流程的第一步目标是把DN值换算成有物理意义的辐亮度或反射率。换算公式不复杂L Gain × DN Bias其中Gain和Bias是定标系数通常可以从影像自带的元数据XML文件、头文件或者数据分发方的文档里查到。不同传感器的定标系数命名可能不一样有的叫RadianceGain/RadianceBias有的叫ScaleFactor/Offset需要仔细分辨。换算成辐亮度以后很多分析任务还需要进一步把它转换成地表反射率或表观反射率。如果不引入大气参数至少可以计算大气层顶的反射率TOA Reflectance公式为ρ π × L × d² / (ESUN × cosθ)这里的d是日地距离天文单位ESUN是太阳光谱辐照度平均值θ是太阳天顶角。d可以通过儒略日查表或近似公式计算ESUN在传感器定标文档里一般有对应波段的值太阳天顶角则在元数据里直接给出。我在处理某高分多光谱影像时元数据XML里除了定标系数还有太阳高度角、卫星侧摆角这些参数后面都会用到。实操建议是第一步永远是把元数据完整读一遍不要想当然地套用上一个项目的参数。我见过有的数据源在不同批次更新了定标系数旧文件里的数值已经过时直接用会导致反射率整体偏差百分之几对定量分析影响很大。2.2 大气校正给影像摘掉那层“滤镜”辐射定标做完影像还是笼罩在大气散射和吸收造成的“滤镜”里。大气校正是为了反演真实的地表反射率。主流方法分两类物理模型和经验模型。物理模型的代表是6S模型和FLAASH模块。它们通过大气辐射传输方程模拟太阳辐射在大气中的吸收、散射过程需要输入成像时间、中心经纬度、传感器高度、大气模式热带、中纬度夏季/冬季、亚极地等、气溶胶模型等参数。FLAASH还可以利用影像里有水汽吸收波段的特性自动反演水汽含量精度比较高适合定量反演类应用。经验模型的代表是DOS暗像元减法和QUAC。DOS假设影像里存在暗像元如清洁水体、浓密阴影这些像元的地表反射率接近于零那么传感器在这些像元接收到的信号主要来自大气路径辐射把它减去就能近似消除大气影响。QUAC更省事不需要输入复杂的成像参数直接从影像统计特征估算适合快速预处理和参数不全的情况。我自己的经验是如果是做NDVI、叶面积指数反演、水体参数反演这类定量任务尽量用物理模型如果只是做地物分类或者目标识别且来不及精细处理QUAC已经很够用。需要强调的是水体区域的暗像元选择要特别小心如果影像里水体已经发生藻华或泥沙含量高它不是真正的暗像元把它当作零反射率去减会让整景影像都偏暗。城区由于多次散射严重物理模型的误差也偏大处理城市高分影像时不要对反射率绝对值期望过高。2.3 几何校正与正射校正解决“位置不准”的问题几何问题在高分影像上是绕不开的。卫星在轨成像时传感器姿态角变化、地形起伏、地球曲率都会造成像点位置偏移。尤其是山区高差带来的投影差快则几十米慢则上百米不处理的话影像叠加到矢量边界上会明显错位。高分影像现在普遍带有RPC有理多项式系数文件。正射校正的做法是把RPC文件、精细DEM和影像一起输入到校正工具中通过有理多项式模型纠正地形引起的像点位移输出指定投影坐标系通常选UTM或国家坐标系的正射影像。操作上要注意几个细节DEM分辨率最好不低于影像分辨率的10倍。比如2米分辨率影像用30米分辨率的DEM就够了但如果到了亚米级30米DEM会造成地形边缘的残余误差建议换更高精度的DEM。输出像元大小要和原始影像一致不要随意放大比如0.8米影像如果重采样成2米细节就丢了。如果影像没有RPC文件那就只能用地面控制点做多项式纠正。控制点要选在道路交叉口、田埂角、稳定水体边界这类特征明显的位置均匀分布至少20到30个高阶多项式虽然能降低残差但容易过度扭曲一般用二阶多项式已经是上限。精度评估用RMSE均方根误差比较直观检查点在影像上的预测位置与实测位置之间的偏差控制在1个像元以内是比较理想的。2.4 影像融合让一幅图既有光谱信息又有细节很多高分卫星同时携带全色传感器和多光谱传感器全色影像分辨率高比如0.5米但只有单波段灰度信息多光谱影像有丰富颜色蓝绿红近红外等但分辨率较低比如2米。影像融合就是把两者优势合并输出一幅同时具备高空间分辨率和多光谱信息的影像。常用的融合方法有Brovey变换、PCA变换、Gram-Schmidt变换、NNDiffuse等。高分数据我做得最多的是Gram-Schmidt和NNDiffuse。Gram-Schmidt对光谱保真度控制得比较好适合做真彩色合成和地物解译NNDiffuse对空间细节的保持更锐利适合做精细目标识别。融合前要确认一个关键顺序先做辐射定标、大气校正、几何校正再做融合。如果先融合再校正全色和多光谱的配准误差会被融合过程放大后续校正的效果也会打折。另外融合后要检查波段顺序。很多专业软件在处理时能自动识别波段波长但如果你把融合结果导出成普通GeoTIFF原来看似正常的RGB通道可能会变成BGR顺序直接拿到WebGIS上显示会红蓝对调。全色与多光谱的分辨率比值通常在4:1到10:1之间。比例越大融合越依赖算法猜测空间细节伪纹理现象越明显。如果比值超过10:1我通常建议直接用多光谱结果做分析别强行融合。3. 预处理工具链选型从图形界面到可复用的脚本方案3.1 三个主流方案的特点对比遥感预处理工具的选择很大程度上取决于数据量多少、是否要重复执行、以及使用者的编程基础。我常用的方案有三类各自特点对比如下方案优点缺点适合场景专业遥感桌面软件功能完善、参数界面成熟、交互直观授权成本高、批处理能力弱单景或少量的精细处理、新手学习开源GIS桌面软件免费、插件生态丰富、社区教程多部分流程需要手工点选、自动化不够彻底中小规模处理、没有预算限制的项目Python/命令行脚本完全可复现、批处理能力强、便于集成到业务线上手有门槛、需要自己理解数据结构大批量数据、需要定期重复的流程如果你只是处理三五景影像用专业遥感桌面软件点一遍半小时也就完成了完全没必要写脚本。但如果项目是几十上百景数据或者每个月都要处理新增影像脚本化几乎是唯一选项。我现在的做法是GUI用于单景探索性处理和参数试算脚本用于正式批处理。3.2 Python脚本化高置信、可复现的预处理骨架Python生态里处理栅格数据最常用的库是rasterio和GDAL底层绑定库。下面我给出一个辐射定标加重投影加裁剪的脚本骨架它可以直接套用到大部分光学高分影像上。import numpy as np import rasterio as rio from rasterio.warp import reproject, calculate_default_transform from rasterio.crs import CRS def rad_to_toa(src_path, dst_path, gain, bias, sun_zenith, d_au): 把DN值转成大气层顶反射率简化版不含ESUN按波段修正 cos_sza np.cos(np.deg2rad(sun_zenith)) with rio.open(src_path) as src: profile src.profile.copy() profile.update(dtypefloat32, driverGTiff) with rio.open(dst_path, w, **profile) as dst: for band in range(1, src.count 1): dn src.read(band).astype(float32) rad dn * gain[band - 1] bias[band - 1] ref (np.pi * rad * d_au ** 2) / (cos_sza) dst.write(ref, band) def resample_to_crs(src_path, dst_path, dst_crsEPSG:32650, resolutionNone): 把影像重投影到目标坐标系 with rio.open(src_path) as src: transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds, resolutionresolution ) kwargs src.profile.copy() kwargs.update({ crs: dst_crs, transform: transform, width: width, height: height, driver: GTiff }) with rio.open(dst_path, w, **kwargs) as dst: for band in range(1, src.count 1): reproject( sourcerio.band(src, band), destinationrio.band(dst, band), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingrio.enums.Resampling.bilinear )这段代码里有两个地方值得说明。辐射定标里我把ESUN修正省略了严格做TOA反射率时不同波段的ESUN不一样需要从查找表读取并逐个波段计算。这里简化是为了把流程讲清楚实际项目中务必把ESUN带上。重投影我用了双线性重采样适合连续表面数据如果做完分类再重投影就应该用最近邻避免改变类别值。还有一点很关键重投影不等于正射校正。如果影像有明显的山区地形畸变光做重投影是纠正不了投影差的必须先利用RPC和DEM做正射校正。rasterio本身支持RPC校正但需要配合高程数据代码会比上面复杂不少。3.3 批处理前把“元数据纪律”建立起来预处理脚本写好后最怕的是不同批次影像的参数不一样批量跑出来一堆废数据。我建议项目启动前先建立一个数据清单CSV至少包含以下字段影像文件名、采集日期、传感器名称、中心经纬度、投影坐标系EPSG代码、太阳高度角、定标增益偏移文件名、是否存在RPC文件、是否存在云掩膜文件。批处理时脚本读取这个CSV来驱动遇到缺失字段直接报警而不是猜默认值。这样做一次虽然繁琐但能避免后面成百上千景数据里混入坐标错误、参数错误、日期错误的隐患。我在一个多时相项目中就因为忽略了一景旧批次影像的定标参数差异导致最终合成产品出现一条明显的“接缝”后来花了大量时间重新排查才找到问题。4. 五类高频翻车现场我的排查过程与处理办法4.1 坐标漂移投影信息错认的常见来源现象是影像叠加到矢量边界上偏离了。新入行的同事第一反应是“数据源精度不行”。但实际上大部分坐标漂移是投影信息识别错误导致的。有一次处理某高分影像文件扩展名是.tif没有附带投影说明同事直接按WGS84经纬度坐标系加载结果影像和研究区矢量边界偏了大概几百米。我打开元数据XML看到里面其实写明了UTM Zone 50N而我们的矢量边界在CGCS2000坐标系下两者基准面和投影方式都不一样当然对不上。排查链路是这样的先用GDAL命令行工具查看影像的完整信息。gdalinfo 影像文件名.tif重点看Coordinate Reference System字段、Origin和Pixel Size字段。如果Origin坐标数值范围和你预期差很远大概率是投影识别错了。再用专业软件把影像和已知参考图层叠加多选几个特征点对比基本就能定位问题。处理办法很简单重新设定或转换投影到统一坐标系即可但要在数据入口就把这道关把住否则后续分析的误差都白算了。4.2 位深与波段顺序看着正常算出来全错高分影像常见16位整型存储但分为有符号和无符号两种。数据说明文件里写的是“16-bit signed”如果你用无符号整型方式读取负值会被解释成很大的正数影像上就会出现大量异常亮点。这个问题从“影像显示”层面很难看出来直方图却会拖出一条长长的尾巴。我排查的步骤是先查元数据里的Pixel Type字段再用工具读一遍直方图和min/max统计值。如果最小值是-1000这种说明是有符号数据没按有符号读如果最小值是0且最大值接近65535但影像上很多像元过曝要警惕是否真的包含了负值。波段顺序是另一个高频坑。不少高分数据的波段存储顺序是BGR而不是RGB因为传感器原始波段按波长排列时蓝、绿、红、近红外都有存储时可能把蓝放在第一个。如果你直接在通用GIS软件里把它当作红绿蓝显示颜色会红蓝互换整个影像看起来非常诡异。解决办法是在读取时按BandName或波长信息排序或者在处理脚本里显式指定波段映射。4.3 大气校正后影像发灰或偏色有次做大气校正参数填的是中纬度夏季大气模式和城市气溶胶模型跑完之后整景影像灰蒙蒙的暗像元反而比校正前更亮。排查过程我按三个方向来第一检查暗像元选择。我用的影像里有一大片清洁水体按理是完美的暗像元但那一景水体边缘有薄云残留自动暗像元统计把它也算了进去导致大气路径辐射被高估。改用掩膜掉云区后的水体统计结果立刻正常。第二检查气溶胶参数。当时输入的气溶胶光学厚度AOT偏高模型过度补偿了散射地物反射率被压暗。换成从影像自身反演的AOT后反射率恢复到了合理区间。第三如果以上都没问题那就只是显示问题。大气校正后的反射率值是0到1的小数直接拉伸显示会显得很暗或者发灰。我通常会把反射率乘以10000转成整型存储再做2%线性拉伸保证可视化效果。注意这只是显示拉伸不影响数值。4.4 融合结果出现“彩色噪声”与伪纹理高分影像融合后最常遇到的问题是颜色失真和像棉花一样发花的纹理。这不是融合算法本身不行多半是输入数据没对齐。排查链路我建议这样走先检查全色和多光谱影像之间的配准偏移。在GIS软件里把两个波段叠加用半透明方式查看道路、屋顶边缘是否完全重合。哪怕只差半个像元融合后边缘就会出现双重轮廓放大看特别明显。解决办法是用配准工具先把两者的几何偏差修掉再做融合。再检查融合时的波段选择。Gram-Schmidt融合需要指定一个低分辨率多光谱波段作为参考波段通常选近红外波段效果比较稳因为它携带的光谱信息最丰富跟全色波段的亮度分布更接近。如果选了波段1蓝波段融合结果往往颜色偏冷、偏暗。出现这种问题时重新选参考波段再跑一遍就够了。最后检查直方图匹配。有些融合工具支持先做直方图匹配让全色波段和多光谱亮度范围一致能有效减少融合后颜色断层。对于NNDiffuse这类方法输入数据的辐射定标精度直接影响融合质量所以前面定标那一步不要偷懒。4.5 大范围影像内存爆掉分块处理的正确姿势高分影像动辄几个GB16位多光谱四波段10000×10000像元就是800MB如果做浮点运算再翻几倍一次性读进内存很容易把机器拖垮。我遇到过一次连着跑几景影像时程序直接被杀掉排查后发现是脚本用了src.read()一次性读取全图。解决办法是分块处理rasterio支持按窗口读取和写入。大致思路是import rasterio as rio window_size 1024 with rio.open(src_path) as src: profile src.profile.copy() profile.update(compresslzw) with rio.open(dst_path, w, **profile) as dst: for row in range(0, src.height, window_size): for col in range(0, src.width, window_size): window rio.windows.Window(col, row, min(window_size, src.width - col), min(window_size, src.height - row)) data src.read(windowwindow) # 在这里做各类算法处理 dst.write(data, windowwindow)分块大小根据内存量调整一般1024×1024像元是比较稳妥的选择。另外输出GeoTIFF时开启LZW或DEFLATE压缩能显著减少磁盘占用。高分影像有的原始存储根本没压缩一景8GB的影像开压缩后能压到3GB以下存储压力小很多。5. 预处理完成不等于处理完成质量检查与数据交付5.1 七项质量检查清单预处理跑完之后千万别急着拿去训练模型或者做分类先花十五分钟做一遍质量检查。我总结了一份七项清单每次都按这个来检查项检查方法验收标准辐射值范围统计每波段min/max/mean/std反射率0到1或0到10000之间无明显截断直方图形状画直方图连续光滑无断层、无长尾异常影像时相信息对比元数据日期与结果文件日期完全一致地理范围在GIS里叠加行政边界或矢量影像范围覆盖目标区域无明显缺口像元大小查看Pixel Size与设计分辨率一致边界错位叠加参考底图或高精度矢量特征点偏差小于1个像元波段合理性读取典型地物光谱曲线水体低反射、植被近红外高反射有次检查时发现同一批次数据里有几景影像的经纬度范围与邻近影像差了0.1度最后定位到是某次下载时把数据导错了因为提前做了检查才没把错误数据带进后续流程。5.2 一个快速验证流程除了检查清单我习惯再做一个快速验证整体花不了几分钟但对结果信心提升非常大。先统计所有波段的均值。对于多光谱影像近红外波段均值应该明显高于可见光波段如果近红外比红光还低要么是波段顺序错了要么是处理过程把波段搞混了。再选三个典型地物点建立光谱剖面分别是植被区、清洁水体和裸土。植被在近红外波段反射率要显著抬升水体从红光到近红外一路下降裸土各波段相对平缓。如果这三条曲线符合物理规律说明辐射处理没问题。最后做一个简单的叠加目视检查在GIS软件里把正射后的影像叠在上一期历史影像或高精度矢量上沿着一条公路断面画一条线观察地物边界是否对齐。这一步能直观看出几何精度是否达标。5.3 输出格式与数据组织建议预处理结果建议统一输出为GeoTIFF格式坐标参考系统统一并在文件名中带上处理版本号。如果后续要拼接多景影像用GDAL的VRT虚拟格式先建索引不要一开始就拼一个大文件能节省大量I/O时间。数据目录组织我通常采用这样的结构原始数据目录只读权限保留元数据和原始文件中间产品目录存放定标、校正后的分步结果方便回溯问题出在哪一步最终产品目录存放经过质量检查的正射融合影像和QA报告。QA报告里记录每景影像的处理参数、检查结果、发现的问题及处理方式。这个习惯帮过大忙。有次三个月后要复现某批数据的处理过程靠着当时的QA报告很快确认了所有参数省去了重新摸索的时间。预处理是遥感分析里最不“性感”但又最不能省的一环。最后分享三点个人体会第一拿到新格式的影像数据先花半小时读元数据比后面返工节省的时间多得多第二不要迷信某套处理流程和参数放之四海而皆准每个数据源、每个地区、每个季节都可能需要微调第三脚本允许的话加上参数记录至少每一景处理完留一个log文件。我刚开始也经历过拿原始图直接分析、结果一团糟的阶段后来养成了预处理后检查的习惯项目的返工率明显低了很多。如果你还没吃到教训希望在读到这篇文章的此刻就开始重视它如果已经踩过坑我想你会对我的总结会心一笑。遥感数据处理这条路慢就是快。
阅读完成 · 觉得有帮助？
