做遥感这行拿到一景新数据的第一反应往往不是打开看看而是先翻元数据确认级别——因为同样叫一景影像L1级产品和L2级产品之间差的可能就是三四个小时的预处理工作量。我见过不少刚入门的同事把原始DN值直接丢进分类器跑出来的结果图看着挺花哨一验证精度惨不忍睹也见过有人对已经做过大气校正的产品又跑了一遍校正最后反射率被压得面目全非。遥感图像预处理这件事步骤本身在手册上都列得清清楚楚真正难的是判断哪些步骤必须做、哪些可以跳过、顺序怎么排。这篇就把辐射定标、大气校正、几何校正、正射校正、影像融合、镶嵌裁剪、重采样与降噪这一整条链路拆开讲重点放在每一步的判定依据和实操里那些文档不会写的细节。不管你是刚接触遥感图像处理的学生还是需要把流程工程化的开发者都能从中找到可以直接照做的东西。1. 先别急着点运行预处理清单要从数据级别倒推1.1 一景原始数据的出厂状态遥感卫星下传的原始数据本质上是一堆带位置标签的数字矩阵。传感器在轨记录的是入瞳处的辐射能量经过量化后存成整数这就是我们常说的DN值Digital Number。这个数字本身没有任何物理单位1和100的区别只代表相对强弱你不能拿它去和别的东西比较也不能直接代入物理模型反演。与此同时卫星平台在飞、地球在转、镜头有畸变所以每个像元的坐标也不是严格对齐地面某个点的。再加上大气层这一层毛玻璃挡在传感器和地表之间太阳光穿过它两次能量被散射和吸收掉一部分落进镜头的信号已经掺了杂质。这就是原始数据的三重不干净数值没有物理意义、位置有偏差、信号被大气污染。预处理的所有步骤本质都是在逐个消除这三类问题。理解这一点很关键因为一旦你知道每一步在解决哪个问题就能判断眼前这景数据到底缺不缺这一步而不是机械地照着流程图从第一个框点到最后一个框。1.2 产品级别与预处理工作量的对应关系不同级别的产品出厂时已经完成了不同程度的加工。下面这张表是我自己整理的心智模型实际项目中对照着看能省下不少重复劳动。产品级别已完成处理你还需要做的典型适用场景L0仅解包、分幅全部流程需要自行做辐射定标的研究L1辐射校正级辐射定标、部分几何粗校正大气校正、正射校正定量反演、指数计算L2几何精校正级系统几何校正、投影变换正射校正山区尤其需要制图、变化检测L3正射级几何精校正DEM正射大气校正若未做直接叠加分析L4及以上包含大气校正的地表反射率基本可直接使用时序分析、快速应用需要说明的是各家的级别命名并不统一有的把大气校正放在L2A有的放在L3所以在实际动手前一定要读产品说明书别只看级别编号。我踩过一次坑以为某个L2产品已经做过大气校正结果打开一看反射率值域明显不对一查文档才发现它只是几何精校正级白白浪费了半天时间返工。1.3 三类典型任务的预处理需求差异预处理清单不是越长越好它应该由下游任务倒推。做土地覆盖分类核心诉求是地物之间的可分性几何精度要够、辐射一致性要好但绝对反射率的准确性要求没那么苛刻大气校正可以用简化方法甚至跳过同景影像内部相对关系不变即可。做植被生化参数反演或者水质监测那就必须走完整的辐射定标加大气校正因为你要拿反射率去拟合物理模型绝对值的偏差会直接传导到结果上。做底图制图或者变化检测几何精度和色彩一致性是第一位正射校正和镶嵌的色彩均衡比大气校正更值得投入精力。这里有个常见的误区把所有任务都按最高标准处理。结果是一景数据跑两小时大气校正最后只是用来做一张示意性分布图性价比极低。反过来该做定标的地方图省事跳过后面所有定量分析全部作废损失更大。判断标准其实很简单——问自己一句下游算法用到的是像元之间的相对关系还是像元的绝对物理值前者可以偷懒后者一步都不能少。2. 辐射定标把DN值翻译成有物理意义的辐亮度2.1 DN、辐亮度、表观反射率的三层关系辐射定标要完成的是从DN到辐亮度的转换。辐亮度的单位是 W/(m²·sr·μm)表示单位面积、单位立体角、单位波长范围内接收到的辐射能量这才是有物理量纲的东西。再往上走一步是表观反射率TOA reflectance它把太阳辐照度和太阳天顶角考虑进来做了一次归一化好处是不同时间、不同传感器、不同光照条件下的影像可以放在同一个尺度上比较。这三层的关系可以这么理解DN值像是没有刻度的温度计你能看出今天比昨天热但说不出具体多少度辐亮度是装上了刻度能读数了但不同温度计之间还可能存在系统偏差表观反射率则是统一了标准所有人的读数可以直接对比。大气校正是在表观反射率的基础上再往下走一层剥离大气影响得到地表反射率。2.2 定标系数从哪来、怎么用定标公式本身很简单通常写成线性形式辐亮度 增益 × DN 偏移。增益和偏移这两个系数不在影像文件里而是在随数据下发的元数据文件中一般是一个XML或者文本文件字段名类似 RADIANCE_MULT_BAND_x 和 RADIANCE_ADD_BAND_x每个波段一组。有些传感器给的是绝对定标系数有些给的是相对定标系数需要配合太阳参数换算读元数据的时候一定要看清楚字段含义别想当然。用Python处理的话思路是先用栅格库把影像读成数组再逐波段套公式。伪代码大致是这样import numpy as np # 假设 bands 是一个形状为 (波段数, 高, 宽) 的数组 # gain 和 offset 是从元数据中解析出来的列表 radiance np.zeros_like(bands, dtypenp.float32) for i in range(bands.shape[0]): radiance[i] bands[i].astype(np.float32) * gain[i] offset[i]看起来三行搞定但有两个细节必须注意。第一原始数据一般是整型uint16居多如果不先转成浮点就直接乘增益结果会被截断尤其在增益很小的时候整幅图会变成一片零。第二定标后的数据类型要选 float32 而不是 float64前者足够用且内存占用减半处理大幅面影像时这个差别非常明显。如果要继续算表观反射率还需要太阳天顶角和日地距离这两个参数。公式是反射率 π × 辐亮度 × 日地距离² / (太阳辐照度 × cos(太阳天顶角))。太阳天顶角从元数据里读日地距离可以根据成像日期用经验公式估算也可以直接从元数据取。这里最容易错的是角度单位——有的元数据给的是度有的是弧度有的是天顶角有的是高度角代入前务必确认差一个90度或者一个π/180结果会离谱到你自己都认不出来。2.3 定标后必须做的合理性检查定标完成不等于做对了。我养成的习惯是每次定标后做三项检查第一统计各波段的辐亮度值域和该传感器的典型范围对比如果整体偏大或偏小一个数量级八成是系数用错了波段或者单位换算漏了第二检查是否有负值理论上辐亮度不应为负出现负值通常意味着偏移量符号搞反了第三随机取几个像元用元数据里的系数手工算一遍和程序结果核对。还有一点容易被忽略定标后的无效值处理。原始数据里常有填充值比如0或者65535定标后这些值会变成奇怪的数字混进统计里导致后续的拉伸和增强效果失真。正确做法是在定标前就用掩膜标记出无效像元或者定标后把它们统一设为NaN。这个小动作在出图阶段体现得特别明显不做的话直方图拉伸会被极端值拉偏整幅图看起来灰蒙蒙的。3. 大气校正最贵的一步也是最容易被滥用的一步3.1 大气到底改了哪些值太阳辐射从太空到地表再从地表回到传感器要穿过两层大气。上行和下行过程中气溶胶散射和分子散射会把一部分辐射偏折到别的方向同时水汽、臭氧、二氧化碳等成分会在特定波段吸收能量。最终传感器收到的信号里除了真正来自目标像元的反射还混进了两部分杂质一部分是大气散射直接进入视场的程辐射另一部分是邻近地物散射进来的邻近效应。这三部分叠加的后果是短波段蓝、绿被抬得最厉害因为散射强度与波长的四次方成反比蓝光波长短散射强所以原始影像看起来总是偏蓝偏灰。做归一化植被指数的时候如果没做大气校正蓝波段虚高红光波段也被抬升算出来的指数值偏低而且不同日期影像之间的可比性很差。这就是为什么做时序分析时大气校正几乎是必选项。3.2 辐射传输模型法与暗像元法的取舍大气校正的方法大致分两派。一派是辐射传输模型法基于物理方程模拟光在大气中的传输过程输入大气参数气溶胶光学厚度、水汽含量、臭氧浓度等输出地表反射率。这类方法精度高能处理复杂情况但计算量大而且对大气参数的准确性很敏感。常用的大气辐射传输模型有6S、MODTRAN等很多商业和开源软件都把它们封装成了现成模块。另一派是影像基方法最典型的是暗像元法假设影像中存在反射率接近零的暗地物深水体、浓密植被阴影把它的值当作程辐射估计量从整幅图里减掉。两者的选择不是看谁更高级而是看你的数据条件和精度要求。下面这张对照表是我实际项目中的判断依据维度辐射传输模型法暗像元法输入要求需要大气参数或辅助数据只需影像自身计算耗时单景几十分钟到数小时秒级到分钟级精度高绝对反射率可信中等相对关系可信适用场景定量反演、多时相定量比较快速出图、单景分类失败风险大气参数偏差导致系统性误差影像中无暗地物时严重失真我用暗像元法翻过车处理一景冬季北方影像整幅图被雪覆盖根本找不到真正的暗像元算法把雪地当成暗目标减出来的反射率大面积出现负值。所以用这类方法之前先看看影像里有没有可靠的暗地物没有就老老实实上模型法。3.3 校正结果的验证方式大气校正做完怎么知道对不对最直接的办法是找同步实测数据对比但实际项目中很少有这样理想的条件。退而求其次的做法有几种一是看植被在红光和近红外波段的反射率是否落在合理区间健康植被红光反射率通常在0.03到0.10之间近红外在0.3以上如果算出来红光0.3、近红外0.1那肯定是哪里搞反了二是看水体在近红外波段的反射率是否接近零这是物理上的硬约束三是比较校正前后的直方图正常情况下去除程辐射后短波段整体会下移长波段变化较小。还有一个实用技巧如果研究区里有大面积清澈深水体校正后水体的近红外反射率应该很低且稳定可以作为内部检查点。另外校正后的反射率值域一般落在0到1之间个别地物如云、雪可能略超1如果出现大面积的负值或者大于2的值基本可以判定参数设置有问题。4. 几何校正与正射校正让每个像元回到正确的位置4.1 几何误差的三种来源影像位置对不上的原因有三类。第一类是系统性的来自传感器自身的扫描方式、地球自转、轨道高度变化这类误差有规律可循可以用轨道模型和传感器参数精确计算并消除叫系统几何校正。第二类是随机性的来自卫星姿态的微小抖动、平台振动没有固定规律只能靠地面控制点来拟合。第三类是地形引起的山区尤其严重因为传感器是斜视成像山越高投影差越大同一座山在影像上的位置和它在平面地图上的位置能差出好几个像元。前两类误差在几何精校正里解决第三类必须靠正射校正。很多新手把这两件事混为一谈以为做了几何校正就不用做正射校正结果在平原地区看着还行一进山区叠图就明显错位。4.2 RPC加DEM的正射校正流程现在主流的正射校正走的是有理多项式系数RPC路线。卫星厂商在数据里附带了RPC文件它描述的是一组多项式关系把地面三维坐标映射到影像像素坐标。校正的时候对输出影像的每个像元根据它的平面坐标和DEM提供的高程用RPC反算出对应的原始影像位置再做灰度重采样。所以正射校正的三样必备材料是影像本身、RPC参数、覆盖研究区的DEM。DEM的分辨率选择有讲究。经验上DEM网格不宜粗于影像地面分辨率的3倍否则地形起伏的刻画不够细致正射效果会打折扣。但也不是越精细越好高分辨率DEM数据量大、处理慢而且如果DEM本身含有噪声反而会引入新的位置抖动。我一般按影像分辨率的2到3倍来选够用且不至于拖慢流程。另外要留意的是RPC参数的精度。有些数据的RPC是星上粗略估算的直接用它做正射平面误差可能达到十几米。这种情况需要用少量控制点对RPC做区域网平差优化把系统偏差先修掉再做正射。这一步在很多教程里被省略了但实际工程中如果精度要求高是跳不过去的。4.3 控制点选取的实操细节需要人工选控制点时有几个细节直接决定成败。控制点要均匀分布不能全挤在影像一角最好按四角加中心的布局来布数量上每景15到30个为宜太少拟合不稳太多边际收益递减还增加人工成本。控制点要选位置明确、不随时间变化的目标道路交叉口、水库大坝拐角、大型建筑轮廓都是好选择避开农田边界、河岸线、植被边缘这些季节性变化明显的地方。选点时的定位精度也需要注意。如果你是从已有正射影像上取参考点要确保参考影像本身的精度高于待校正影像的目标精度否则就是用一把不准的尺子去量东西。还有一个实操技巧选完点后先看残差把残差明显偏大的点剔除再重新解算通常一轮下来中误差能降不少。残差大于两个像元的点除非有充分理由否则不要保留。4.4 山区与城市的表现差异同一套参数在平原和山区效果完全不同。平原地区地形平坦正射校正前后的差别可能只有一两个像元肉眼几乎看不出来山区则可能出现十几甚至几十个像元的位移山体阴影的方向和长度都会变化。城市区域的问题在于高层建筑正射校正只能处理地形起伏处理不了建筑物本身的高度引起的投影差所以高楼会出现倾斜现象。这个问题在正射阶段无解需要专门的真正射处理代价很高一般项目里很少做。我的建议是拿到数据先看研究区地形用DEM算一下高程标准差起伏小于50米的区域正射校正的收益有限时间紧可以简化处理起伏超过200米的正射校正必须认真做而且DEM质量要重点把关。5. 融合、镶嵌、裁剪把数据整理成能直接用的底图5.1 全色锐化算法的效果差别多光谱影像波段多但分辨率低全色影像分辨率高但只有一个波段全色锐化pansharpening就是把两者的优势合起来生成既有高分辨率又有多光谱信息的产品。算法层面大致分三类基于成分替换的比如主成分变换、强度-色调-饱和度变换、基于多分辨率分析的比如小波变换、拉普拉斯金字塔、以及基于模型优化的比如各种变分方法。实际选型时我关注的指标是光谱保真度和空间细节注入的平衡。成分替换类方法空间细节注入得猛视觉上很锐利但容易出现光谱扭曲植被区域的颜色可能偏色明显。多分辨率分析方法光谱保持得好一些但空间细节可能不够锐利。基于模型的方法综合表现好但计算量大。如果你的下游任务是分类或者指数计算光谱保真度优先如果只是做视觉底图那锐利一点更好看。融合有个前提容易被忽略参与融合的多光谱和全色数据必须是同一时刻、同一视角获取的而且要精确配准。配准误差超过一个全色像元融合结果就会出现明显的重影看着像照片没对上焦。融合前务必做一次配准检查别等融合完了才发现问题。5.2 镶嵌的接边与色彩均衡研究区跨越多景影像时就要做镶嵌。简单的做法是按地理范围拼起来但这样拼出来的图会出现明显的拼接缝因为不同景的成像时间、光照条件、大气状态都不一样色调差异肉眼可见。要出能看的图必须做色彩均衡。色彩均衡有几种思路全局的直方图匹配把各景的直方图都往参考景上靠局部的方法在重叠区域做渐变过渡还有基于辐射归一化的方法从物理层面消除差异。实际操作中全局直方图匹配简单有效但遇到地物构成差异大的相邻景比如一景是城市一景是山地强行匹配会把某一景拉得失真。这种情况更适合用重叠区统计的方法只在重叠带内做线性调整。接边线的设计也有讲究。理想情况下接边线应该沿着自然地物边界走比如河流、山脊、道路这样拼缝最不显眼。手工描接边线费时但效果比直线切割好太多。另外要确保接边线两侧没有云、没有阴影突变否则拼出来像打了补丁。5.3 裁剪、掩膜与无效值处理裁剪看着最简单实际也有坑。按矢量边界裁剪时要确认矢量和影像的坐标系一致坐标系不匹配导致的裁剪错位是新手最常见的错误之一。裁剪后边缘会出现不规则形状无效区域需要统一设为NoData而这个NoData值的选择要和后续处理软件兼容有些软件默认用0但0在某些波段可能是合法值会造成误判。掩膜处理还包括云掩膜、水体掩膜、阴影掩膜等。云掩膜一般用专门的云检测算法生成也可以结合短波红外和热红外波段做阈值判断。这里提醒一句掩膜不是加上就完事后续所有统计和建模都要把掩膜区域排除否则云区的异常值会污染整个分析结果。我见过有人做完云掩膜直接丢了最后统计出来的均值偏高查了半天才发现是云没排除干净。6. 重采样、降噪与云处理收尾工序里的隐性成本6.1 重采样核的选择只要涉及几何变换、投影转换或者分辨率调整就一定会用到重采样。常用的核有三种最近邻、双线性、三次卷积。最近邻是直接取最近像元的值不产生新值所以最适合分类结果、掩膜这类离散数据的重采样缺点是边缘会出现锯齿。双线性用周围四个像元做加权平均结果平滑但会损失高频细节。三次卷积用周围十六个像元做更复杂的加权锐度保留得好但可能产生过冲在亮度突变的地方出现不自然的亮边暗边。选择逻辑很清晰连续型数据反射率、温度、高程优先考虑双线性或三次卷积离散型数据分类标签、掩膜、质量标记必须用最近邻。这条规则看着简单但在实际流程里特别容易搞混尤其是在一个脚本里批量处理多种数据的时候。我一般会在代码里显式判断数据类型防止误用。还有一个细节重采样会改变像元的实际覆盖范围多次重采样会造成信息损失累积。能用一次重采样解决的不要分成几次做。6.2 噪声类型与对应处理遥感影像的噪声主要有几类。条带噪声常见于推扫式传感器表现为整列或者整行的亮度异常处理方法有矩匹配、直方图均衡、频域滤波等。椒盐噪声表现为孤立的亮点或暗点中值滤波效果最好。随机噪声则表现为整体的颗粒感可以用高斯滤波或者非局部均值去噪。去噪的核心矛盾是噪声抑制和细节保留。滤波窗口开得大噪声压得干净但边缘和小目标也被磨掉了。我的经验是宁可轻一点也不要把细节抹平因为很多地物特征恰恰体现在高频部分。另外去噪最好在辐射定标之后、大气校正之前做因为噪声会影响大气校正参数估计的稳定性但如果噪声很轻微也可以留到最后统一处理避免反复重采样。6.3 云检测的工程化处理云是光学遥感绕不开的问题。薄云会降低对比度厚云直接遮挡地表云影还会让地表反射率系统性偏低。云检测的方法很多简单阈值法用蓝波段和短波红外的亮度组合判断机器学习方法则用训练好的分类器逐像元判别。工程上更在意的是处理策略是直接掩掉云区还是做云去除填补直接掩膜简单可靠适合云量小的情况代价是数据缺失。云去除填补比如时序合成、多时相插值能补出完整影像但填补区域的可信度需要标注清楚。我在做长时序分析时的做法是先用云掩膜把污染像元标出来再在时间维度上做插值填补同时对填补区域加一个质量标记波段后续分析时根据需求决定是否采用。这样既保证了数据完整性又不至于把插值结果当成真实观测。7. 踩坑实录顺序、位深、元数据这三处最容易翻车7.1 步骤顺序错了会怎样预处理的顺序不是随意排列的它遵循一个内在逻辑先修正数值再修正位置最后做拼接和裁剪。辐射定标和大气校正必须在几何变换之前因为大气校正的模型参数和像元的光照几何相关一旦做了重采样像元的原始观测几何就模糊了。几何校正必须在镶嵌之前否则各景之间的配准无从谈起。裁剪放在最后因为裁剪会丢弃边缘信息如果先裁剪再镶嵌边界处会缺数据。我见过最典型的错误是先做大气校正再做辐射定标顺序完全颠倒结果自然是错的也见过先镶嵌后做大气校正的因为各景成像条件不同同一套大气参数套在整幅镶嵌图上效果很差正确的做法是逐景校正后再镶嵌。7.2 数据类型与位深的陷阱整型和浮点的转换时机是另一个高频出错点。原始数据为了省存储一般是16位整型一旦开始做定标和校正就必须转成32位浮点。如果中途为了省空间又转回整型精度损失是不可逆的。反过来如果一开始就全程用64位浮点内存占用翻倍大幅面影像很容易撑爆内存。还有运算符的隐式类型问题。两个整型数组相除在某些语言里会做整除结果全是0或者1这个坑在各种脚本语言里都出现过。凡是涉及除法的计算务必先显式转成浮点。7.3 批量处理与元数据丢失批量处理时最容易出问题的是元数据。栅格文件里除了像素值还有坐标系、投影参数、NoData值、波段波长、成像时间等信息很多手写的处理脚本只顾着算像素忘了把这些信息写回去输出结果打开一看坐标全是默认值叠加到地图上位置完全不对。解决办法是用成熟的栅格库读写它们会自动继承大部分元数据但投影和NoData这两项最好显式设置一遍别依赖默认行为。另一个批量处理的隐患是参数硬编码。一景数据的参数写死在脚本里能跑通换个传感器就全错。稳妥的做法是把增益偏移、太阳角度这些参数从每景的元数据文件里动态读取脚本只负责逻辑不负责数值。最后分享一点个人体会整套预处理流程里真正花时间的从来不是点几下按钮而是判断这一步到底需不需要做、参数该怎么定、结果对不对。我现在的习惯是每做完一步就存一个中间结果命名里带上处理内容和日期虽然占点硬盘但出了问题能快速定位到是哪一步引入的比从头重跑高效得多。遥感图像预处理没有万能模板把每一步背后的物理逻辑想明白比记住操作步骤有用得多。
做遥感这行拿到一景新数据的第一反应往往不是打开看看而是先翻元数据确认级别——因为同样叫一景影像L1级产品和L2级产品之间差的可能就是三四个小时的预处理工作量。我见过不少刚入门的同事把原始DN值直接丢进分类器跑出来的结果图看着挺花哨一验证精度惨不忍睹也见过有人对已经做过大气校正的产品又跑了一遍校正最后反射率被压得面目全非。遥感图像预处理这件事步骤本身在手册上都列得清清楚楚真正难的是判断哪些步骤必须做、哪些可以跳过、顺序怎么排。这篇就把辐射定标、大气校正、几何校正、正射校正、影像融合、镶嵌裁剪、重采样与降噪这一整条链路拆开讲重点放在每一步的判定依据和实操里那些文档不会写的细节。不管你是刚接触遥感图像处理的学生还是需要把流程工程化的开发者都能从中找到可以直接照做的东西。1. 先别急着点运行预处理清单要从数据级别倒推1.1 一景原始数据的出厂状态遥感卫星下传的原始数据本质上是一堆带位置标签的数字矩阵。传感器在轨记录的是入瞳处的辐射能量经过量化后存成整数这就是我们常说的DN值Digital Number。这个数字本身没有任何物理单位1和100的区别只代表相对强弱你不能拿它去和别的东西比较也不能直接代入物理模型反演。与此同时卫星平台在飞、地球在转、镜头有畸变所以每个像元的坐标也不是严格对齐地面某个点的。再加上大气层这一层毛玻璃挡在传感器和地表之间太阳光穿过它两次能量被散射和吸收掉一部分落进镜头的信号已经掺了杂质。这就是原始数据的三重不干净数值没有物理意义、位置有偏差、信号被大气污染。预处理的所有步骤本质都是在逐个消除这三类问题。理解这一点很关键因为一旦你知道每一步在解决哪个问题就能判断眼前这景数据到底缺不缺这一步而不是机械地照着流程图从第一个框点到最后一个框。1.2 产品级别与预处理工作量的对应关系不同级别的产品出厂时已经完成了不同程度的加工。下面这张表是我自己整理的心智模型实际项目中对照着看能省下不少重复劳动。产品级别已完成处理你还需要做的典型适用场景L0仅解包、分幅全部流程需要自行做辐射定标的研究L1辐射校正级辐射定标、部分几何粗校正大气校正、正射校正定量反演、指数计算L2几何精校正级系统几何校正、投影变换正射校正山区尤其需要制图、变化检测L3正射级几何精校正DEM正射大气校正若未做直接叠加分析L4及以上包含大气校正的地表反射率基本可直接使用时序分析、快速应用需要说明的是各家的级别命名并不统一有的把大气校正放在L2A有的放在L3所以在实际动手前一定要读产品说明书别只看级别编号。我踩过一次坑以为某个L2产品已经做过大气校正结果打开一看反射率值域明显不对一查文档才发现它只是几何精校正级白白浪费了半天时间返工。1.3 三类典型任务的预处理需求差异预处理清单不是越长越好它应该由下游任务倒推。做土地覆盖分类核心诉求是地物之间的可分性几何精度要够、辐射一致性要好但绝对反射率的准确性要求没那么苛刻大气校正可以用简化方法甚至跳过同景影像内部相对关系不变即可。做植被生化参数反演或者水质监测那就必须走完整的辐射定标加大气校正因为你要拿反射率去拟合物理模型绝对值的偏差会直接传导到结果上。做底图制图或者变化检测几何精度和色彩一致性是第一位正射校正和镶嵌的色彩均衡比大气校正更值得投入精力。这里有个常见的误区把所有任务都按最高标准处理。结果是一景数据跑两小时大气校正最后只是用来做一张示意性分布图性价比极低。反过来该做定标的地方图省事跳过后面所有定量分析全部作废损失更大。判断标准其实很简单——问自己一句下游算法用到的是像元之间的相对关系还是像元的绝对物理值前者可以偷懒后者一步都不能少。2. 辐射定标把DN值翻译成有物理意义的辐亮度2.1 DN、辐亮度、表观反射率的三层关系辐射定标要完成的是从DN到辐亮度的转换。辐亮度的单位是 W/(m²·sr·μm)表示单位面积、单位立体角、单位波长范围内接收到的辐射能量这才是有物理量纲的东西。再往上走一步是表观反射率TOA reflectance它把太阳辐照度和太阳天顶角考虑进来做了一次归一化好处是不同时间、不同传感器、不同光照条件下的影像可以放在同一个尺度上比较。这三层的关系可以这么理解DN值像是没有刻度的温度计你能看出今天比昨天热但说不出具体多少度辐亮度是装上了刻度能读数了但不同温度计之间还可能存在系统偏差表观反射率则是统一了标准所有人的读数可以直接对比。大气校正是在表观反射率的基础上再往下走一层剥离大气影响得到地表反射率。2.2 定标系数从哪来、怎么用定标公式本身很简单通常写成线性形式辐亮度 增益 × DN 偏移。增益和偏移这两个系数不在影像文件里而是在随数据下发的元数据文件中一般是一个XML或者文本文件字段名类似 RADIANCE_MULT_BAND_x 和 RADIANCE_ADD_BAND_x每个波段一组。有些传感器给的是绝对定标系数有些给的是相对定标系数需要配合太阳参数换算读元数据的时候一定要看清楚字段含义别想当然。用Python处理的话思路是先用栅格库把影像读成数组再逐波段套公式。伪代码大致是这样import numpy as np # 假设 bands 是一个形状为 (波段数, 高, 宽) 的数组 # gain 和 offset 是从元数据中解析出来的列表 radiance np.zeros_like(bands, dtypenp.float32) for i in range(bands.shape[0]): radiance[i] bands[i].astype(np.float32) * gain[i] offset[i]看起来三行搞定但有两个细节必须注意。第一原始数据一般是整型uint16居多如果不先转成浮点就直接乘增益结果会被截断尤其在增益很小的时候整幅图会变成一片零。第二定标后的数据类型要选 float32 而不是 float64前者足够用且内存占用减半处理大幅面影像时这个差别非常明显。如果要继续算表观反射率还需要太阳天顶角和日地距离这两个参数。公式是反射率 π × 辐亮度 × 日地距离² / (太阳辐照度 × cos(太阳天顶角))。太阳天顶角从元数据里读日地距离可以根据成像日期用经验公式估算也可以直接从元数据取。这里最容易错的是角度单位——有的元数据给的是度有的是弧度有的是天顶角有的是高度角代入前务必确认差一个90度或者一个π/180结果会离谱到你自己都认不出来。2.3 定标后必须做的合理性检查定标完成不等于做对了。我养成的习惯是每次定标后做三项检查第一统计各波段的辐亮度值域和该传感器的典型范围对比如果整体偏大或偏小一个数量级八成是系数用错了波段或者单位换算漏了第二检查是否有负值理论上辐亮度不应为负出现负值通常意味着偏移量符号搞反了第三随机取几个像元用元数据里的系数手工算一遍和程序结果核对。还有一点容易被忽略定标后的无效值处理。原始数据里常有填充值比如0或者65535定标后这些值会变成奇怪的数字混进统计里导致后续的拉伸和增强效果失真。正确做法是在定标前就用掩膜标记出无效像元或者定标后把它们统一设为NaN。这个小动作在出图阶段体现得特别明显不做的话直方图拉伸会被极端值拉偏整幅图看起来灰蒙蒙的。3. 大气校正最贵的一步也是最容易被滥用的一步3.1 大气到底改了哪些值太阳辐射从太空到地表再从地表回到传感器要穿过两层大气。上行和下行过程中气溶胶散射和分子散射会把一部分辐射偏折到别的方向同时水汽、臭氧、二氧化碳等成分会在特定波段吸收能量。最终传感器收到的信号里除了真正来自目标像元的反射还混进了两部分杂质一部分是大气散射直接进入视场的程辐射另一部分是邻近地物散射进来的邻近效应。这三部分叠加的后果是短波段蓝、绿被抬得最厉害因为散射强度与波长的四次方成反比蓝光波长短散射强所以原始影像看起来总是偏蓝偏灰。做归一化植被指数的时候如果没做大气校正蓝波段虚高红光波段也被抬升算出来的指数值偏低而且不同日期影像之间的可比性很差。这就是为什么做时序分析时大气校正几乎是必选项。3.2 辐射传输模型法与暗像元法的取舍大气校正的方法大致分两派。一派是辐射传输模型法基于物理方程模拟光在大气中的传输过程输入大气参数气溶胶光学厚度、水汽含量、臭氧浓度等输出地表反射率。这类方法精度高能处理复杂情况但计算量大而且对大气参数的准确性很敏感。常用的大气辐射传输模型有6S、MODTRAN等很多商业和开源软件都把它们封装成了现成模块。另一派是影像基方法最典型的是暗像元法假设影像中存在反射率接近零的暗地物深水体、浓密植被阴影把它的值当作程辐射估计量从整幅图里减掉。两者的选择不是看谁更高级而是看你的数据条件和精度要求。下面这张对照表是我实际项目中的判断依据维度辐射传输模型法暗像元法输入要求需要大气参数或辅助数据只需影像自身计算耗时单景几十分钟到数小时秒级到分钟级精度高绝对反射率可信中等相对关系可信适用场景定量反演、多时相定量比较快速出图、单景分类失败风险大气参数偏差导致系统性误差影像中无暗地物时严重失真我用暗像元法翻过车处理一景冬季北方影像整幅图被雪覆盖根本找不到真正的暗像元算法把雪地当成暗目标减出来的反射率大面积出现负值。所以用这类方法之前先看看影像里有没有可靠的暗地物没有就老老实实上模型法。3.3 校正结果的验证方式大气校正做完怎么知道对不对最直接的办法是找同步实测数据对比但实际项目中很少有这样理想的条件。退而求其次的做法有几种一是看植被在红光和近红外波段的反射率是否落在合理区间健康植被红光反射率通常在0.03到0.10之间近红外在0.3以上如果算出来红光0.3、近红外0.1那肯定是哪里搞反了二是看水体在近红外波段的反射率是否接近零这是物理上的硬约束三是比较校正前后的直方图正常情况下去除程辐射后短波段整体会下移长波段变化较小。还有一个实用技巧如果研究区里有大面积清澈深水体校正后水体的近红外反射率应该很低且稳定可以作为内部检查点。另外校正后的反射率值域一般落在0到1之间个别地物如云、雪可能略超1如果出现大面积的负值或者大于2的值基本可以判定参数设置有问题。4. 几何校正与正射校正让每个像元回到正确的位置4.1 几何误差的三种来源影像位置对不上的原因有三类。第一类是系统性的来自传感器自身的扫描方式、地球自转、轨道高度变化这类误差有规律可循可以用轨道模型和传感器参数精确计算并消除叫系统几何校正。第二类是随机性的来自卫星姿态的微小抖动、平台振动没有固定规律只能靠地面控制点来拟合。第三类是地形引起的山区尤其严重因为传感器是斜视成像山越高投影差越大同一座山在影像上的位置和它在平面地图上的位置能差出好几个像元。前两类误差在几何精校正里解决第三类必须靠正射校正。很多新手把这两件事混为一谈以为做了几何校正就不用做正射校正结果在平原地区看着还行一进山区叠图就明显错位。4.2 RPC加DEM的正射校正流程现在主流的正射校正走的是有理多项式系数RPC路线。卫星厂商在数据里附带了RPC文件它描述的是一组多项式关系把地面三维坐标映射到影像像素坐标。校正的时候对输出影像的每个像元根据它的平面坐标和DEM提供的高程用RPC反算出对应的原始影像位置再做灰度重采样。所以正射校正的三样必备材料是影像本身、RPC参数、覆盖研究区的DEM。DEM的分辨率选择有讲究。经验上DEM网格不宜粗于影像地面分辨率的3倍否则地形起伏的刻画不够细致正射效果会打折扣。但也不是越精细越好高分辨率DEM数据量大、处理慢而且如果DEM本身含有噪声反而会引入新的位置抖动。我一般按影像分辨率的2到3倍来选够用且不至于拖慢流程。另外要留意的是RPC参数的精度。有些数据的RPC是星上粗略估算的直接用它做正射平面误差可能达到十几米。这种情况需要用少量控制点对RPC做区域网平差优化把系统偏差先修掉再做正射。这一步在很多教程里被省略了但实际工程中如果精度要求高是跳不过去的。4.3 控制点选取的实操细节需要人工选控制点时有几个细节直接决定成败。控制点要均匀分布不能全挤在影像一角最好按四角加中心的布局来布数量上每景15到30个为宜太少拟合不稳太多边际收益递减还增加人工成本。控制点要选位置明确、不随时间变化的目标道路交叉口、水库大坝拐角、大型建筑轮廓都是好选择避开农田边界、河岸线、植被边缘这些季节性变化明显的地方。选点时的定位精度也需要注意。如果你是从已有正射影像上取参考点要确保参考影像本身的精度高于待校正影像的目标精度否则就是用一把不准的尺子去量东西。还有一个实操技巧选完点后先看残差把残差明显偏大的点剔除再重新解算通常一轮下来中误差能降不少。残差大于两个像元的点除非有充分理由否则不要保留。4.4 山区与城市的表现差异同一套参数在平原和山区效果完全不同。平原地区地形平坦正射校正前后的差别可能只有一两个像元肉眼几乎看不出来山区则可能出现十几甚至几十个像元的位移山体阴影的方向和长度都会变化。城市区域的问题在于高层建筑正射校正只能处理地形起伏处理不了建筑物本身的高度引起的投影差所以高楼会出现倾斜现象。这个问题在正射阶段无解需要专门的真正射处理代价很高一般项目里很少做。我的建议是拿到数据先看研究区地形用DEM算一下高程标准差起伏小于50米的区域正射校正的收益有限时间紧可以简化处理起伏超过200米的正射校正必须认真做而且DEM质量要重点把关。5. 融合、镶嵌、裁剪把数据整理成能直接用的底图5.1 全色锐化算法的效果差别多光谱影像波段多但分辨率低全色影像分辨率高但只有一个波段全色锐化pansharpening就是把两者的优势合起来生成既有高分辨率又有多光谱信息的产品。算法层面大致分三类基于成分替换的比如主成分变换、强度-色调-饱和度变换、基于多分辨率分析的比如小波变换、拉普拉斯金字塔、以及基于模型优化的比如各种变分方法。实际选型时我关注的指标是光谱保真度和空间细节注入的平衡。成分替换类方法空间细节注入得猛视觉上很锐利但容易出现光谱扭曲植被区域的颜色可能偏色明显。多分辨率分析方法光谱保持得好一些但空间细节可能不够锐利。基于模型的方法综合表现好但计算量大。如果你的下游任务是分类或者指数计算光谱保真度优先如果只是做视觉底图那锐利一点更好看。融合有个前提容易被忽略参与融合的多光谱和全色数据必须是同一时刻、同一视角获取的而且要精确配准。配准误差超过一个全色像元融合结果就会出现明显的重影看着像照片没对上焦。融合前务必做一次配准检查别等融合完了才发现问题。5.2 镶嵌的接边与色彩均衡研究区跨越多景影像时就要做镶嵌。简单的做法是按地理范围拼起来但这样拼出来的图会出现明显的拼接缝因为不同景的成像时间、光照条件、大气状态都不一样色调差异肉眼可见。要出能看的图必须做色彩均衡。色彩均衡有几种思路全局的直方图匹配把各景的直方图都往参考景上靠局部的方法在重叠区域做渐变过渡还有基于辐射归一化的方法从物理层面消除差异。实际操作中全局直方图匹配简单有效但遇到地物构成差异大的相邻景比如一景是城市一景是山地强行匹配会把某一景拉得失真。这种情况更适合用重叠区统计的方法只在重叠带内做线性调整。接边线的设计也有讲究。理想情况下接边线应该沿着自然地物边界走比如河流、山脊、道路这样拼缝最不显眼。手工描接边线费时但效果比直线切割好太多。另外要确保接边线两侧没有云、没有阴影突变否则拼出来像打了补丁。5.3 裁剪、掩膜与无效值处理裁剪看着最简单实际也有坑。按矢量边界裁剪时要确认矢量和影像的坐标系一致坐标系不匹配导致的裁剪错位是新手最常见的错误之一。裁剪后边缘会出现不规则形状无效区域需要统一设为NoData而这个NoData值的选择要和后续处理软件兼容有些软件默认用0但0在某些波段可能是合法值会造成误判。掩膜处理还包括云掩膜、水体掩膜、阴影掩膜等。云掩膜一般用专门的云检测算法生成也可以结合短波红外和热红外波段做阈值判断。这里提醒一句掩膜不是加上就完事后续所有统计和建模都要把掩膜区域排除否则云区的异常值会污染整个分析结果。我见过有人做完云掩膜直接丢了最后统计出来的均值偏高查了半天才发现是云没排除干净。6. 重采样、降噪与云处理收尾工序里的隐性成本6.1 重采样核的选择只要涉及几何变换、投影转换或者分辨率调整就一定会用到重采样。常用的核有三种最近邻、双线性、三次卷积。最近邻是直接取最近像元的值不产生新值所以最适合分类结果、掩膜这类离散数据的重采样缺点是边缘会出现锯齿。双线性用周围四个像元做加权平均结果平滑但会损失高频细节。三次卷积用周围十六个像元做更复杂的加权锐度保留得好但可能产生过冲在亮度突变的地方出现不自然的亮边暗边。选择逻辑很清晰连续型数据反射率、温度、高程优先考虑双线性或三次卷积离散型数据分类标签、掩膜、质量标记必须用最近邻。这条规则看着简单但在实际流程里特别容易搞混尤其是在一个脚本里批量处理多种数据的时候。我一般会在代码里显式判断数据类型防止误用。还有一个细节重采样会改变像元的实际覆盖范围多次重采样会造成信息损失累积。能用一次重采样解决的不要分成几次做。6.2 噪声类型与对应处理遥感影像的噪声主要有几类。条带噪声常见于推扫式传感器表现为整列或者整行的亮度异常处理方法有矩匹配、直方图均衡、频域滤波等。椒盐噪声表现为孤立的亮点或暗点中值滤波效果最好。随机噪声则表现为整体的颗粒感可以用高斯滤波或者非局部均值去噪。去噪的核心矛盾是噪声抑制和细节保留。滤波窗口开得大噪声压得干净但边缘和小目标也被磨掉了。我的经验是宁可轻一点也不要把细节抹平因为很多地物特征恰恰体现在高频部分。另外去噪最好在辐射定标之后、大气校正之前做因为噪声会影响大气校正参数估计的稳定性但如果噪声很轻微也可以留到最后统一处理避免反复重采样。6.3 云检测的工程化处理云是光学遥感绕不开的问题。薄云会降低对比度厚云直接遮挡地表云影还会让地表反射率系统性偏低。云检测的方法很多简单阈值法用蓝波段和短波红外的亮度组合判断机器学习方法则用训练好的分类器逐像元判别。工程上更在意的是处理策略是直接掩掉云区还是做云去除填补直接掩膜简单可靠适合云量小的情况代价是数据缺失。云去除填补比如时序合成、多时相插值能补出完整影像但填补区域的可信度需要标注清楚。我在做长时序分析时的做法是先用云掩膜把污染像元标出来再在时间维度上做插值填补同时对填补区域加一个质量标记波段后续分析时根据需求决定是否采用。这样既保证了数据完整性又不至于把插值结果当成真实观测。7. 踩坑实录顺序、位深、元数据这三处最容易翻车7.1 步骤顺序错了会怎样预处理的顺序不是随意排列的它遵循一个内在逻辑先修正数值再修正位置最后做拼接和裁剪。辐射定标和大气校正必须在几何变换之前因为大气校正的模型参数和像元的光照几何相关一旦做了重采样像元的原始观测几何就模糊了。几何校正必须在镶嵌之前否则各景之间的配准无从谈起。裁剪放在最后因为裁剪会丢弃边缘信息如果先裁剪再镶嵌边界处会缺数据。我见过最典型的错误是先做大气校正再做辐射定标顺序完全颠倒结果自然是错的也见过先镶嵌后做大气校正的因为各景成像条件不同同一套大气参数套在整幅镶嵌图上效果很差正确的做法是逐景校正后再镶嵌。7.2 数据类型与位深的陷阱整型和浮点的转换时机是另一个高频出错点。原始数据为了省存储一般是16位整型一旦开始做定标和校正就必须转成32位浮点。如果中途为了省空间又转回整型精度损失是不可逆的。反过来如果一开始就全程用64位浮点内存占用翻倍大幅面影像很容易撑爆内存。还有运算符的隐式类型问题。两个整型数组相除在某些语言里会做整除结果全是0或者1这个坑在各种脚本语言里都出现过。凡是涉及除法的计算务必先显式转成浮点。7.3 批量处理与元数据丢失批量处理时最容易出问题的是元数据。栅格文件里除了像素值还有坐标系、投影参数、NoData值、波段波长、成像时间等信息很多手写的处理脚本只顾着算像素忘了把这些信息写回去输出结果打开一看坐标全是默认值叠加到地图上位置完全不对。解决办法是用成熟的栅格库读写它们会自动继承大部分元数据但投影和NoData这两项最好显式设置一遍别依赖默认行为。另一个批量处理的隐患是参数硬编码。一景数据的参数写死在脚本里能跑通换个传感器就全错。稳妥的做法是把增益偏移、太阳角度这些参数从每景的元数据文件里动态读取脚本只负责逻辑不负责数值。最后分享一点个人体会整套预处理流程里真正花时间的从来不是点几下按钮而是判断这一步到底需不需要做、参数该怎么定、结果对不对。我现在的习惯是每做完一步就存一个中间结果命名里带上处理内容和日期虽然占点硬盘但出了问题能快速定位到是哪一步引入的比从头重跑高效得多。遥感图像预处理没有万能模板把每一步背后的物理逻辑想明白比记住操作步骤有用得多。
