我手里刚好有一景 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 COMPRESS=DEFLATE，文件会小不少，读起来也不慢。
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 2>&1 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 --outfile=valid_mask.tif \ --calc="A==0 ? 1 : 0" --NoDataValue=0这样后续反演时就能很方便地套 mask 了。
5.3 与深度学习云检测组合使用
这几年基于深度学习的云检测在常规厚云上的召回率确实漂亮，比如 Sentinel-2 上经常用的 s2cloudless，对厚云的识别非常稳。但它在薄卷云和云影判断上并不总是优于 Fmask，而且很多训练模型是拿欧洲场景做的，拿到干旱区或寒区有时会水土不服。我在实际项目中会做混合策略：先用 Fmask 打底，保证物理含义可解释的结果，再用深度学习方法作为独立结果并行对比，最后取交集或者按置信度投票。这样得到的掩膜通常比单个模型更稳，也不容易被人质疑算法是否只是记忆了某种场景。
5.4 批量统计：让积累的影像自动出报表
当手头有几十上百景影像时，单张看显然不现实。我会把 Fmask 跑完后生成一个 CSV 统计表，记录每景影像的云占比、有效像元占比、平均云概率等信息。哪怕不做很复杂的分析，这个表也能让你快速筛选出可用影像。以前有个项目需要挑出连续三年的无云影像做土地覆盖变化，要是没有这种表，光靠人工翻缩略图得翻一天。配合前文的脚本，把cloud_percent、valid_area_km2追加到同一表格，执行完看一眼就知道哪些数据能进下一步。这就是 Fmask 这类自动云检测工具带给实际生产的最大价值：它把场景清理从手工劳动变成了流水线环节。
