接手这个任务前先说说我为什么会写这篇东西。前阵子临时接到一个活儿需要快速处理一景当天上午刚获取的多光谱影像把云盖住的地方都标出来。我习惯性打开处理管线结果发现一套标准流程上来就要求输入前后十几天的时序数据——手里只有单张影像的时候整套方案直接瘫痪。后来我把思路切回到Fmask算法才意识到一个问题很多云检测教程默认你手里有时间序列但实际工作中大量需求恰恰是单张影像也要能出云掩膜。这篇博文就围绕这个场景展开。我会把Fmask云层检测的原理、单张影像上的实现方式和踩过的坑一次说清楚适合正在做遥感影像预处理、需要从单景多光谱图像里提取云掩膜的同行也适合刚入门、想知道Fmask到底在算什么的同学。1. 为什么单张是云检测的分水岭1.1 云检测的三大流派里只有物理规则法能单打独斗业内做云检测的方法大致分三类。第一类是时序分析法思路很直观同一块地在一个月里会被卫星拍好几次植被、土壤、水面反射率虽然会变但连续两期影像的差异通常是渐变的云不一样它今天盖住这儿、明天飘到那儿前后期变化是跳变的。这个思路做长时间序列很稳但它的前提就是必须有至少两期、最好是五六期以上的影像。手里只有单张的时候这类方法连启动都做不到。第二类是机器学习法包括随机森林、深度学习语义分割那一套。这类方法效果很好但需要大量标注样本而且模型训练时的影像季节、太阳高度角、地表类型都会影响泛化。你拿一个在欧洲训练好的云检测模型去处理热带雨林上空的影像误检率分分钟教会你做人。单张场景确实可以用它但前提是你手里得有个训练好的、且匹配当前地理区域的模型。第三类就是Fmask这种物理规则法。Fmask的全称是Function of mask最早是Zhu和Woodcock在2012年提出的针对Landsat影像设计。它的核心思想不是统计上像不像云而是物理上是不是云——把云的亮度高、温度低、光谱平坦等物理特征写成一连串的判断规则再组合成概率出来。正因为判断依据是物理量和光谱量而不是时序变化量所以它天然可以在单张影像上运行。这也是我最终选择它的根本原因。1.2 哪些实际任务逼着你在单张影像上做云掩膜以前我总觉得单张做云检测属于特殊需求后来干活多了才发现这其实是常态。灾后快速评估就是典型场景台风过境、洪水退去之后你拿到的往往就是灾后第二天的单景影像要的是当天立刻出产掩膜不可能等半个月攒齐时序再开工。还有一种情况是历史存档影像的修复。很多早期归档的多光谱影像只有单景存底没有配套时序。你要把这些数据纳入统一处理流程第一步就得给每张影像单独做云掩膜。另外像无人机多光谱、国产高分辨率卫星这类任务经常因为重访周期长或者任务定制原因同一区域很难快速攒出多期影像。这种一锤子买卖的影像要想用起来单张云检测就是绕不开的前置处理。理解了这一点你就会明白为什么Fmask在学术界和工业界沉淀这么多年还有生命力因为它解决的是最底层的刚需。2. Fmask盯着多光谱影像的六个特征2.1 每一个特征背后的物理原因Fmask判断云靠的是六个关键特征。第一个是亮度特征。云层在可见光和近红外波段的反射率通常很高尤其厚云反射率能超过0.5。地面物体里只有雪和部分裸地、人工地表能达到这个量级。第二个是白度特征。云的反射率在蓝、绿、红三个可见光波段之间差异很小看起来是白的而裸地、水体、植被在可见光波段有明显的谱形差异。Fmask计算可见光波段之间的光谱差异差异小、亮度高的像素就被认为是白的。第三个是亮温特征。云顶高度高、温度低热红外波段的亮温通常明显低于多数地表。这个特征在区分云和亮裸地时特别有效——沙漠在可见光里很亮但它热红外亮温高一测就露馅了。第四个是卷云波段特征。Landsat 8和Sentinel-2都有1.37-1.38微米附近的卷云探测波段。这个波段很有意思因为大气水汽对这个波段吸收极强地面信号基本到不了传感器但高空卷云在大气层顶之上它的反射信号能穿过来。所以这个波段一亮基本就是高空薄云而且这种云在可见光波段肉眼常常看不出来。第五个是NDSI雪检测特征。云和雪在可见光波段都是高反射、白光谱最难区分。但雪在短波红外波段反射率骤降云在短波红外依然保持高反射。利用归一化差异雪指数NDSI能把大部分雪像素先挑出去避免后面误判成云。第六个是光谱变异性特征。这个和可见光白度有点相似但参与计算的波段范围更宽。云在多个波段上的反射率曲线平滑而地表因为植被、土壤、阴影混合光谱曲线起伏大。Fmask用这个特征进一步压低裸土和建筑区的误检率。2.2 从特征到云概率阈值不是拍脑袋有了这六个特征之后Fmask并不是简单地大于某个值就是云而是把每个特征都映射成一个云概率值最后把所有概率综合起来。比如亮度高但亮温也很高说明是热的亮地表亮度这一项给云概率加分亮温这一项却会强烈减分两个证据互相制约最后综合概率可能只有0.3不会被判成云。这就是为什么Fmask误检率比那种全波段反射率大于0.5就标记为云的粗暴方法低很多。每个特征的概率映射都需要设定参数。以亮温为例算法里通常会用一个阈值区间亮温低于某值例如273K左右时云概率接近1高于某值例如295K时云概率趋近0中间区间则用线性或高斯曲线过渡。这些参数不是凭空来的是在Landsat影像上大量统计云顶和地表的典型亮温分布后得到的。但要注意这些参数只是通用初始值不同季节、不同纬度、不同地表类型都会有偏差完全照搬不是不行而是你会把不少陆地像素错误地卷进云里。3. 单张Landsat影像的Fmask落地流程3.1 数据准备你需要什么样的Landsat产品先说数据源。Landsat 8的Level 2产品已经包含了地表反射率和地表温度这是做Fmask最省事的输入。单张影像在USGS或相关公开数据池都能找到关键是下载时认准包含这两个数据层的产品级别。以Landsat 8为例你需要拿到这些波段蓝光Band 2、绿光Band 3、红光Band 4、近红外Band 5、短波红外Band 6和Band 7、热红外Band 10或Band 11、卷云波段Band 9。每个波段都是单独的GeoTIFF文件分辨率上可见光波段是30米热红外是100米重采样到30米。处理的第一步就是把所有波段统一到同一分辨率、同一范围。我习惯用rasterio写一个小函数批量读取顺手把无效值通常用0或-9999标记转成NaNimport rasterio import numpy as np def read_band(path): with rasterio.open(path) as src: band src.read(1).astype(np.float32) band[band 0] np.nan return band # 假设文件名遵循Landsat C2 L2的命名规则 blue read_band(LC08_xxx_B2.TIF) green read_band(LC08_xxx_B3.TIF) red read_band(LC08_xxx_B4.TIF) nir read_band(LC08_xxx_B5.TIF) swir1 read_band(LC08_xxx_B6.TIF) swir2 read_band(LC08_xxx_B7.TIF) cirrus read_band(LC08_xxx_B9.TIF) bt read_band(LC08_xxx_ST_B10.TIF)注意一点Level 2产品的浮点反射率文件名字段里通常带有SR_地表温度则是ST_开头别下错了。3.2 方案A直接用官方Fmask工具包如果不想从零实现算法可以走现成工具这条更实际的路。Fmask官方提供了MATLAB版本社区也移植了C和Python版本其中Python版本能直接读取GeoTIFF文件输入。它的输入是几类关键特征文件TOA反射率Landsat 8上是六个波段、亮温、卷云波段、以及太阳高度角、太阳方位角、卫星天顶角、卫星方位角这些观测几何信息。这些几何信息通常存在MTL元数据文件里。工具包的调用逻辑是先写一个输入配置文件把所有波段的路径、元数据里的角度参数填进去然后跑核心程序最后输出一个云和云阴影掩膜文件。掩膜里每个像素有编号0是清晰地表1是云2是云阴影3是水体4是雪。这里要提醒一句官方工具对输入波段顺序、单位有严格约定反射率默认是[0, 1]区间而不是百分比亮温默认是开尔文单位。不少人跑出来的掩膜全图都是云八成就是单位没换算。3.3 方案BPython教学版自实现用现成工具够快但如果你想真正搞懂Fmask在算什么还是推荐自己实现一版简化流程。我自己写的教学版本分四步代码量不大但每一步都对应原算法的核心思想。第一步计算潜在的云候选像素def potential_cloud(blue, green, red, nir, swir1, bt): # 亮度条件近红外反射率较高 bright (nir 0.30) | (red 0.30) # 温度条件亮温偏低 cold bt 275.0 # 白度条件可见光波段差异小 vis_mean (blue green red) / 3.0 whiteness np.abs(blue - vis_mean) np.abs(green - vis_mean) np.abs(red - vis_mean) white whiteness 0.07 return bright cold white, bright, cold, white第二步用NDSI把雪剔除def exclude_snow(swir1, green): ndsi (green - swir1) / (green swir1 1e-6) # NDSI大于0.4通常视为雪不参与云判定 return ndsi 0.4第三步卷云波段检测光学上不可见的薄云def cirrus_cloud(cirrus): return cirrus 0.01第四步把各证据综合成云掩膜def simple_fmask(bands, bt): blue, green, red, nir, swir1, cirrus bands cand, bright, cold, white potential_cloud(blue, green, red, nir, swir1, bt) snow exclude_snow(swir1, green) cloud_score (bright.astype(np.float32) * 0.4 cold.astype(np.float32) * 0.3 white.astype(np.float32) * 0.3) cloud (cand ~snow) | (cirrus 0.01) cloud_mask cloud (cloud_score 0.5) return cloud_mask注意这个教学版本的阈值是我按经验拍的比如白度阈值0.07、卷云阈值0.01、综合分数0.5。原版Fmask里这些特征都有更精细的概率函数不是简单的布尔值相加。这套代码的定位是帮你打通物理特征→判断逻辑这条路真要搞生产还是用方案A的完整工具或者去读原论文、按照概率模型逐步复刻。4. 单张场景的四个坑我都替你踩过4.1 坑一雪和云长得很像NDSI救不了全部第一次用简化版跑单张Landsat影像我选了一个冬季场景结果山区的云掩膜上出现了一大片误检。查了半天发现是雪亮度高、白度高、亮温也低三条特征和云完美重叠。NDSI确实能区分大部分雪但它对阴影里的雪和正在融化的薄雪并不可靠这两类雪在短波红外波段的衰减不明显NDSI也不够高。后来我的处理方法是给NDSI阈值做双轨制对于亮度条件极强近红外反射率大于0.4的像素把NDSI阈值从0.4上调到0.45压缩雪的生存空间对于亮度和白度都一般的像素保持0.4不动避免误伤真正的云。这种细节上的差异处理比改一个全局阈值要有效得多。4.2 坑二没有时序约束后阴影匹配变得保守Fmask原来有个很讨巧的设计——云阴影的确认会借助时序影像里的暗像素来验证。有前后期影像时云阴影位置的像素在前后期通常都是暗的这个约束能过滤掉不少误判。单张影像没有这个信息云的阴影只能靠太阳高度角、太阳方位角和卫星观测几何来估算投射位置。这意味着两件事第一你会碰到很多云匹配上了但阴影匹配不上的像素按原版规则它们会被置为清晰地表但实际观察里那片区域明显发暗第二为了不让阴影漏检太夸张我会把阴影匹配的搜索范围放宽一些。放宽的代价是误检增加尤其是山地阴影和沟壑阴影它们的几何形态和云阴影非常接近。单张场景下没有时序兜底只能靠DEM辅助剔除——山地的坡度朝向和阴影方向做一次交叉验证能去掉大部分假阳性。4.3 坑三缺热红外波段时温度项失效怎么办不是所有多光谱影像都有热红外波段尤其是无人机的多光谱载荷和部分高分辨率商业卫星。缺失了温度特征云和亮裸地的区分直接少了一条关键证据。我做过一组对照实验同一个区域有BT参与判断时裸地误检率在5%左右去掉BT之后误检率直接飙到18%主要误检来自矿区、干河床、城市边缘的亮屋顶。解决方案是把近红外和短波红外的约束加严。亮裸地通常在短波红外反射率特别高而云在短波红外虽然也高但相对可见光来说增幅没那么夸张。把短波红外反射率的云概率权重调低同时把亮度阈值往上提一档可以找回一部分误检。但说实话没有温度项的情况下Fmask的精度天花板就是比完整版低一截这也是物理现实没办法完全绕过。4.4 坑四亮地表误检阈值往哪个方向调亮地表误检是单张场景最常见的抱怨包括白色屋顶、干燥盐碱地、新修的水泥路。它们的问题是既亮又白天然适合云的可见光判定标准。我第一次跑城市区域时写字楼玻璃幕墙的亮斑全被框成了云整幅掩膜根本没法用。后来我的调参经验是这类误检的热红外亮温其实明显高于云顶所以只要温度项还在把冷阈值从275K往下压到270K就能甩掉一大片。但要注意夏天高纬度地区的卷云顶温也可能在270K以上阈值压太狠会把真云漏掉。所以正确做法是先跑一版默认参数把掩膜叠加到RGB影像上目视检查误检像素的位置再反过来微调。我一直觉得遥感处理里的参数调整本质上就是一次在和你的影像数据对话的过程没有哪个值能一劳永逸。5. 掩膜做出来以后怎么评估怎么打磨5.1 有QA波段时的像素级评估Landsat的Level 2产品自带的QA_Pixel波段里包含了云和水体的参考标记虽然它本身也有误差但用来做快速评估足够了。我的做法是在整景影像上随机采样500个像素手工目视确定真实类别云还是非云然后和Fmask结果做一张混淆矩阵。这里有个容易忽略的细节采样时不能只在均匀区域随机撒点要刻意覆盖云边界、薄云区、山区阴影这些麻烦地带。如果你全在厚云中心采样精度当然好看但对实际应用没有意义。我一般按三类区域分层采样厚云区150个点、薄云和云边缘区200个点、无云区150个点这样算出来的召回率和精确率才有参考价值。5.2 没有真值时怎么靠抽样和边缘检查把关单张影像处理的尴尬之处在于很多时候你并没有QA参考波段。这种情况下我的经验是两条腿走路。第一条是人工抽样目检把掩膜叠加在真彩色影像上按5x5网格浏览每个格里随机检查5到10个判为云的像素确认是不是真的云。第二条是边缘一致性检查真正的云边缘应该和影像中云的目视边缘吻合如果掩膜边缘比目视云区边缘宽出好几个像素那多半是阈值偏松了如果掩膜边缘在云中心里缩成一团那是白度阈值太严。这两条经验虽然主观但实际操作里比看混淆矩阵更快。5.3 形态学后处理与掩膜应用的细节算法得到的云掩膜通常会有不少斑驳的孔洞尤其是薄云的中间区域偶尔夹杂着清晰像素。物理上一片云覆盖的区域像素应该连续成块所以我会做一步形态学膨胀再腐蚀的闭运算把小孔补上。反过来对于单独孤立的云像素点我会再做一个腐蚀把它们当成椒盐噪声清理掉。还有一个特别容易被新手忽略的坑云掩膜的边缘必须适度向外扩张。云的边缘存在半透明过渡区反射率混入了地表信号算法常常漏判。如果掩膜后续要用于地表反射率产品的生产不向外扩一点最后得到的反射率产品会在云边界处残留一圏异常低值的像素。我一般对云区向外扩张2到3个像素这个值根据影像空间分辨率调整对Landsat 30米影像来说3个像素不算夸张。6. 迁移到其他多光谱传感器时怎么调整6.1 Sentinel-2与Landsat的波段差异很多人手里不只有Landsat还有Sentinel-2它的多光谱波段更丰富分辨率还更高。但迁移到Sentinel-2时要清醒一点它没有热红外波段温度特征直接作废。好在Sentinel-2有个10米分辨率的卷云波段B101.375微米附近对高云的敏感度比Landsat 8的卷云波段更强。Fmask 3.2版本实际上已经加入了Sentinel-2的适配逻辑思路就是重新分配权重大幅降低温度项依赖把卷云波段和短波红外的权重提上来。实际操作中我用Sentinel-2跑下来发现厚云的检测精度比Landsat 8略低但薄卷云的检测能力反而更强。这说明没有一个传感器是十全十美的算法的强项刚好匹配传感器的物理特性才是关键。6.2 其他多光谱载荷的适配思路国内的高分系列、资源系列卫星以及无人机多光谱载荷处理思路也是同一个框架先确定传感器上有哪几个可用波段再看这些波段能触发Fmask的哪些物理特征最后重新调整阈值。我整理过一个快速对照表方便按传感器类型做初步判断传感器情况可用的Fmask特征需要重点调整的参数有热红外可见光短波红外亮度、白度、亮温、NDSI完整套件保持默认逐场景微调无热红外但有卷云波段亮度、白度、卷云、NDSI降低温度权重上调卷云阈值无热红外也无卷云波段亮度、白度、NDSI收紧可见光白度阈值接受精度下降只有RGB多光谱亮度、白度只能做初步云检测输出需人工复核最后说一个非常实用的体会任何适配工作第一步永远是把传感器的光谱响应函数调出来看一眼。Google能搜到的公开波段信息加上自己做的简单光谱模拟就能预先判断云和主要地物的可分离性。有时候你折腾两天调参效果依然很差本质原因是传感器波段设置本身就不支持而不是你的算法写错了。这个判断越早做越省时间。我自己的习惯是在项目目录里保留一份各传感器的波段-特征映射表每次接新数据的传感器时先查表再决定投入多少精力做参数调试。这样既不会迷信默认参数也不会浪费时间在不可行的方案上。