我手里刚好有一景 Landsat 8 OLI 影像云覆盖率 32%。这种数据要是直接拿去反演地表温度或者做地物分类结果基本没法用。多光谱光学遥感最烦人的一点就在这里云层不但遮住了地物信号还会在阴影区域造成假信息所以预处理的第一步永远是回答一个问题——云到底在哪。我这些年处理 Landsat、Sentinel-2 影像时最常用的就是 Fmask全称 Function of mask是 Zhu 和 Woodcock 从 2010 年开始迭代的一套自动云检测算法。今天把单张多光谱图像上跑通 Fmask 的完整过程捋一遍包括算法原理、数据准备、命令行执行、参数调优和踩坑记录适合刚入行的遥感小白也适合想把手头积压影像批量清洗一遍的研究人员和工程师。文章不会涉及太深的理论推导也不会只给一句“去官网下载跑一下”而是尽量把我实操中反复验证过的细节都写出来。1. 单张影像云检测的逻辑内核1.1 云不是一门好做的“目标检测”很多新手一开始觉得云检测很简单像素很亮、很白直接定一个阈值不就完了。我在早期也这么试过最好的一张影像大概能做对 90% 的厚云但碰到薄云和亮地表就大面积翻车。原因很简单云不只是“亮”和“白”的物体它随高度、厚度、相态变化很大。高层卷云很薄很透明在可见光波段的反射还不如一片裸露的干盐滩而低层厚云虽然亮但和积雪、城市高反射屋顶的光谱分布又高度相似。只用单波段亮度或单一阈值来切要么把雪山全判成云要么把半透明的薄云全部漏掉。所以主流云检测算法基本都是多项测试的组合亮度、卷云波段、温度、空间纹理、太阳几何把多个维度的投票结果综合起来才敢下结论。1.2 Fmask 怎么把“云判定”拆成几个可计算的问题Fmask 的设计目标其实很直接利用多光谱图像里云在物理性质上的三条硬特征。第一云顶反射率高尤其在可见光到短波红外这一段云通常比下垫面更亮第二云顶温度低热红外波段的亮度温度能明显把云和大多数地表分开第三卷云由细小冰晶组成对 1.38 微米波长有强散射所以专门设计的卷云波段对高层薄云特别敏感而这一波段里大多数陆表信号基本被水汽吸收掉了。Fmask 把这三条硬特征拆成几组规则检测先用 Otsu 自动阈值法找出候选云像元再对候选区域做连通域聚合接着结合太阳天顶角、方位角和云影搜索来定位云下面的阴影最后再用光谱规则把雪、水体等容易混淆的地物剥离出去。整个过程不需要人工给阈值算法会在每景影像内自适应计算这也是它这些年还能被广泛使用的原因。1.3 单张影像与多时相 Fmask 的差别为什么单独聊Fmask 后来还有时序处理版本利用两景以上影像中像素值变化不大、始终明亮的区域更有可能是云的判断把单张误判的沙漠、雪原给压下去。但我们日常做项目时经常只有单景影像可用比如突发灾害后的应急处理、历史存档中某一景独份数据的研究这时没法等待理想时相来凑时间序列。单张 Fmask 的优势是部署简单、输入只有一个目录不依赖任何其他日期的数据而且对绝大多数中等分辨率应用来说云和厚卷云的检测精度已经够用了。缺点是云影检测会比多时相明显差一些这一点后面我会重点说因为它直接关系到最后整景影像的可用像元统计。对比项单张 Fmask多时相 Fmask输入数据量1 景2 景及以上依赖因素太阳角度、光谱特征光谱特征 不同日期的变化云影精度中等偏差山区易漏检明显更稳适用场合应急、历史存档、单景任务长期序列产品生成2. 一张合格的多光谱影像该如何准备2.1 产品级别选 L1TP别随手拿 L2 表面反射率Fmask 在很多实现里是按 Level-1 的原始 DN 设计并用 Otsu 自适应切阈值的。用 L2 表面反射率产品直接喂容易遇到两个问题一是 L2 产品把热红外、卷云波段删掉或做了不同尺度缩放Fmask 会找不到它要的输入二是 DN 和表面反射率的动态范围差异会让阈值计算跑偏。所以我要么从 USGS EarthExplorer 下载 L1TP 产品要么在 Google Earth Engine 里把原始的 Level-1 波段导出保存。L1TP 是经过辐射校正和地形几何校正的标准产品每个像素的 DN 值能对应到固定的辐射亮度这正好满足 Fmask 的算法假设。老的 L1GT 或者没做几何精校正的数据我一般不用因为生成的云掩膜如果和地表影像空间错位后续做验证时会把问题搅成一锅粥。2.2 波段组织一个目录装下所有必备波段很多人第一次跑 Fmask 失败90% 的原因不是参数而是输入目录乱。Fmask 对输入场景有约定它会在你指定的目录下搜索该场景的各个波段文件并且靠文件名前缀和元数据文件里的场景 ID 去匹配。如果我把文件名随便重命名成 Band2.tif、Band3.tif基本就不可能跑得起来。正确做法是保留传感器原始文件名。以 Landsat 8/9 为例一个标准输入目录大概长这样LC08_L1TP_128044_20231216_20231220_02_T1/ ├── LC08_L1TP_128044_20231216_20231220_02_T1_B2.TIF ├── LC08_L1TP_128044_20231216_20231220_02_T1_B3.TIF ├── LC08_L1TP_128044_20231216_20231220_02_T1_B4.TIF ├── LC08_L1TP_128044_20231216_20231220_02_T1_B5.TIF ├── LC08_L1TP_128044_20231216_20231220_02_T1_B6.TIF ├── LC08_L1TP_128044_20231216_20231220_02_T1_B7.TIF ├── LC08_L1TP_128044_20231216_20231220_02_T1_B9.TIF ├── LC08_L1TP_128044_20231216_20231220_02_T1_B10.TIF └── LC08_L1TP_128044_20231216_20231220_02_T1_MTL.txt这个目录里 B8 全色波段一般用不到B10、B11 两个热红波段具体使用哪个看工具版本。Sentinel-2 的输入目录组织类似但要特别注意它的 L1C 产品一片场景可能被切成多个网格必须先把单景拼接好再送进工具。2.3 MTL / 元数据文件为什么不能省Landsat 的 MTL 文本文件记录太阳高度角、太阳方位角、成像时间、增益偏移等关键信息。Fmask 在阴影检测时要计算云的投影方向靠的正是这些角度。缺少 MTL 会导致两个后果要么程序直接退出要么跑出了一个没有地理方位意义的云影结果。我踩过一次坑从某些第三方网站下载的影像只给了各波段栅格没有附带完整 MTL当时以为云掩膜已经正常输出了结果验证时发现阴影掩膜整片偏了几公里方向还不对。所以数据到位后的第一件事不是急着运行而是先检查目录里有没有完整的 MTL或者对应的 XML 元数据。2.4 快速目检运行前先看一眼数据在真正运行前我会把真彩色合成在 QGIS 里快速拉一遍主要看三点成像区域是否有大面积积雪是否有明显条带或传感器坏行影像边缘是否被云层大面积覆盖。这看似多余但能帮你建立对结果的心理预期——待会儿跑完 Fmask如果雪区被标成云你不能直接说是算法错了因为很多情况下积雪和云在光谱上确实难分。提前知道自己手里是什么难度的场景调参时才不会手足无措。3. 核心执行从命令行到输出掩膜3.1 环境准备与可执行程序获取Fmask 4.1 是目前比较常用的版本官方发布包里同时提供了 MATLAB 源码和编译好的可执行文件。两者的差别在于源码版适合你后续想改规则、看中间变量或者做二次开发可执行文件则适合只想快速拿到云掩膜的工程场景。在 Linux 服务器上我一般直接下编译包然后在 conda 环境里把对应版本的 MATLAB Runtime 装上。版本对不上时会报库缺失这一点要特别留意。比如 Fmask 4.1 如果要求某个特定版本的 MCR直接硬跑会看到类似libmwmfl_*.so not found的错误常规解决办法是把对应 Runtime 加入LD_LIBRARY_PATH。如果是在 Windows 上做小范围实验也可以直接用官方提供的 Windows 版本但大批量处理时我还是更推荐 Linux原因很简单脚本循环、日志管理、定时任务都方便很多而且很多遥感服务器本身就在 Linux 环境里。3.2 运行命令到底怎么写我自己用的时候命令行并不是理解为有很多神奇参数——实际它相当简单就是把影像目录作为第一个参数传入。假设解压后的工具在/opt/Fmask_landsat_4.1那一段最简单的批处理脚本大概如下chmod x /opt/Fmask_landsat_4.1/Fmask_landsat /opt/Fmask_landsat_4.1/Fmask_landsat /data/LC08_L1TP_128044_20231216_20231220_02_T1运行需要一点时间一景 8000×8000 左右的 Landsat 影像在普通 CPU 上一般一到三分钟。如果场景里云和云影很多算法迭代复杂时间会明显拉长。我的习惯是先挑一景小范围或低分辨率测试确认流程通了再批量做不然一上来就开全场景批处理出错了排查成本很高。3.3 结果文件怎么读运行完后输出目录里会多出几个文件。掩膜类文件一般是Fmask4_scene.tif像元值编码为0 表示清晰陆地1 表示云2 表示云影3 表示雪4 表示水。有的版本还会输出一个概率图比如cloud_probability_scene.tif表示每个像元被判定为云的概率值或置信度这个文件用于后处理非常方便。此外还有shadow_shift_x.txt和shadow_shift_y.txt记录的是云影在像素空间上的偏移量说白了就是算法估计的“云把影子投到哪个方向、多远”。这些偏移文件可以用来帮助后续做阴影补偿但单张影像下它只是一个场景级估计不是每个像元都有可靠值使用时别把它当作精确物理量。提示处理历史数据时如果影像覆盖区域很大但云量很碎我建议把掩膜文件的压缩方式再整理一遍否则后续叠加分析时读写压力会很大。比如用gdal_translate转成-co COMPRESSDEFLATE文件会小不少读起来也不慢。3.4 用 Python 快速统计云覆盖率拿到掩膜后我最常做的一件事是统计整景影像里各类像元的占比。用 numpy 读 GeoTIFF 很快from osgeo import gdal import numpy as np ds gdal.Open(Fmask4_LC08_L1TP_128044_20231216_20231220_02_T1.tif) band ds.GetRasterBand(1) arr band.ReadAsArray() label [clear, cloud, shadow, snow, water] counts [(arr i).sum() for i in range(5)] total arr.size for name, cnt in zip(label, counts): print(f{name}: {cnt / total * 100:.2f}%)如果图像很大可以用 gdal_calc 或分块读取但单景 8-bit 掩膜内存占用不大直接读也扛得住。3.5 批量处理的脚本建议批量清洗数据集时我建议写一个简单的 shell 循环同时把中间日志保存下来for dir in /data/landsat/L1TP/*/; do scene$(basename $dir) /opt/Fmask_landsat_4.1/Fmask_landsat $dir /data/logs/fmask_${scene}.log 21 done这个脚本看起来不起眼但省下了大量重复手工操作。除此之外我会在每次批量前先跑一景确认日志里没有报错再全部启动。这样才能保证半夜跑完的数据不会因为一个输入目录缺文件而整体失败。4. 实操中的坑单张影像云检测的常见问题4.1 亮地表被误判成云单张 Fmask 最容易出的问题就是把亮地表误判成云。雪地、盐碱地、城市连片屋顶在可见光波段都和白花花的高反射目标很像。Fmask 内部有专门的雪测试来区分雪和云但雪和云的混合像元、部分融化的雪、或者干燥沙地反射率极高的时候仍会被标成云。我在处理高原影像时就遇到过整片现代冰川被标注为云的情况把掩膜和 RGB 叠加后一眼就能看出问题地物轮廓还在只是像元类别被换了。对这种情况我的处理经验是不要急着改算法参数而是先把明显属于地形连续区域的云像元清洗掉。如果只做后续反演可以把这些误检统一归为“无效像元”虽然损失了一部分有效面积但至少不会把错误类别混进统计分析里。4.2 云影漏检和偏移是单张处理的硬伤云影检测是 Fmask 流程里最难的部分。它的核心思路是把云物体当做一个“投影源”按太阳角度计算出阴影应该在的位置然后在附近区域搜索低亮度匹配。单张影像没有第二期视角可以参考只能靠灰度匹配于是在地形破碎、云影边缘被植被覆盖等情况下漏检很常见。还有一类问题就是山区阴影受到了地形遮蔽算法估算的云高和实际投影位置对不上结果阴影掩膜会整个偏移。另一个让人困惑的现象是云影被检测成了水。云影区域反射率低、色调整体偏暗和水体在某些波段上灰度非常接近所以输出里常出现大片阴影类被归为水的区域。我通常在检查掩膜时会把第 2、4 类一起调出叠加分析不会只盯云类像元。4.3 薄云与卷云阈值策略怎么平衡还有一个高频问题就是对薄云和卷云的检测力度。Fmask 使用卷云波段做阈值如果阈值设得很严很多高海拔的薄雾、大气散射也会被卷入候选云导致云类像元面积虚高如果阈值设得过松又会把真正盖在地物上方的卷云漏掉掩膜上的云区出现大量空洞。这本身就是一对矛盾。我的习惯是把卷云检测更多地看作产品需求问题如果后续要做的是高精度地表反射率反演那么宁可让算法多标一些云也不要漏掉薄云因为漏掉的薄云会污染反射率而不是简单地遮挡地表如果只是做影像接边的云区剔除就可以把阈值放松避免大量正常像元被误杀。说到底阈值不是固定值要根据下游任务倒推。4.4 掩膜验证的快速办法验证云检测结果我不想只靠肉眼看一般会做两类检查。第一类是把掩膜在 QGIS 里调成半透明叠在真彩色影像上随手检查几个异质区域比如雪线、水体、城市边缘。第二类更客观人工均匀随机抽几百个点逐一对比掩膜类别和目视类别。可以用已有的标签或者哨兵云分数等产品做参照。如果后续要做论文建议算一个简单的混淆矩阵from sklearn.metrics import cohen_kappa_score, confusion_matrix # gt: 人工抽样得到的类别, pred: Fmask 同位置类别 gt np.array([0, 1, 1, 0, 2, 0, 1, 0, 0, 3]) pred np.array([0, 1, 2, 0, 2, 0, 1, 0, 1, 3]) print(confusion_matrix(gt, pred)) print(cohen_kappa_score(gt, pred))需要注意公平验证时要给不同地表类型做分层抽样不能光挑看起来顺眼的地方取点否则精度数字会很好看但没有任何说服力。5. 流程落地与扩展从单张掩膜到业务数据5.1 完整流程清单把前面内容汇总成一张流程清单我每次处理都会对照一遍下载 L1TP 产品确认 MTL 和所有波段文件完整。在 QGIS 里目检影像判断场景复杂度记录云量初判。运行 Fmask检查日志没有报错。读取掩膜统计云、阴影、雪、水的像元比例。叠加验证看误检是否集中在哪些区域。根据下游任务决定是否需要后处理比如把阴影类并到无效类。这套流程看起来简单但在我经手的多个项目中大量时间其实花在第 5、6 步而不是运行算法本身。5.2 掩膜后处理清理破碎边界的技巧Fmask 输出的原始掩膜往往有大量细小斑块一个云团边缘会呈锯齿状还会夹杂不少单像素孤岛。做业务时如果希望掩膜更干净可以用形态学开闭运算整理一下但要控制窗口大小窗口太大会把真正的薄云边缘蚕食掉。我用 gdal 加一个小的 Python 脚本或者直接用 GDAL 内置的 sieve 功能去掉面积太小的斑块gdal_sieve.py -pixels 20 -nomask Fmask4_..._T1.tif clean_fmask.tif如果只想生成一张“到底是清晰还是不清”的二值掩膜可以直接用 gdal_calc 把类别 0 置为有效其他全置为无效gdal_calc.py -A clean_fmask.tif --outfilevalid_mask.tif \ --calcA0 ? 1 : 0 --NoDataValue0这样后续反演时就能很方便地套 mask 了。5.3 与深度学习云检测组合使用这几年基于深度学习的云检测在常规厚云上的召回率确实漂亮比如 Sentinel-2 上经常用的 s2cloudless对厚云的识别非常稳。但它在薄卷云和云影判断上并不总是优于 Fmask而且很多训练模型是拿欧洲场景做的拿到干旱区或寒区有时会水土不服。我在实际项目中会做混合策略先用 Fmask 打底保证物理含义可解释的结果再用深度学习方法作为独立结果并行对比最后取交集或者按置信度投票。这样得到的掩膜通常比单个模型更稳也不容易被人质疑算法是否只是记忆了某种场景。5.4 批量统计让积累的影像自动出报表当手头有几十上百景影像时单张看显然不现实。我会把 Fmask 跑完后生成一个 CSV 统计表记录每景影像的云占比、有效像元占比、平均云概率等信息。哪怕不做很复杂的分析这个表也能让你快速筛选出可用影像。以前有个项目需要挑出连续三年的无云影像做土地覆盖变化要是没有这种表光靠人工翻缩略图得翻一天。配合前文的脚本把cloud_percent、valid_area_km2追加到同一表格执行完看一眼就知道哪些数据能进下一步。这就是 Fmask 这类自动云检测工具带给实际生产的最大价值它把场景清理从手工劳动变成了流水线环节。