云检测与云掩膜工程实践

引言

云是光学遥感最大的数据损耗源。全球平均云覆盖率约 67%，在热带雨林与赤道带常年超过 80%，这意味着单景影像里真正可用的像元往往不足三成。如果云像元不剔除，它会同时污染分类结果、扭曲植被指数、制造虚假的变化检测信号，而且这类错误是系统性的，不会因为样本量增大而消失。

云检测的难点在于「云」这个概念本身没有清晰边界。薄卷云、碎积云、云边缘的半透明像元、以及雪与亮色建筑物的混淆，构成一条连续谱而非二值判定。算法要在漏检与误检之间取舍：漏检会把云当成地表，误检会把地表当云丢弃，两者在不同业务里的代价完全不同。农业监测宁可多丢也不愿把云算成作物，而干旱监测则不能把亮色沙地误判为云。

第二个难点是阴影。云阴影与云本体在几何上配对出现，但阴影区亮度极低，容易被误判为水体或地形阴影。掩膜只去云不去阴影，会导致大量暗像元进入统计；去阴影时又容易误伤真实水体。这个配对关系是工程实现中最需要小心处理的部分。

本文按「物理特征到工程实现」的顺序展开。先讲云与阴影的光谱特性，再拆解 Fmask 与 Sen2Cor 的阈值逻辑，然后覆盖机器学习方法、多时相检测、掩膜后处理与形态学操作，最后落到时序合成去云与精度评估。

理解波段与质量波段结构是前提，基础见 遥感影像基础与元数据 ；掩膜与辐射校正的顺序关系见 辐射定标与大气校正 ；去云后的时序数据是 变化检测工程实践 的输入。

目录

1. 云对遥感产品的影响
2. 云的物理特性与光谱特征
3. 阈值法与 Fmask 算法
4. Sen2Cor 与官方掩膜产品
5. 机器学习与云概率方法
6. 多时相与时间序列云检测
7. 掩膜后处理与形态学操作
8. 时序合成与去云重建
9. 精度评估与验证方法

1. 云对遥感产品的影响

云的影响不止是遮挡。厚云完全阻断地表信号，薄云与卷云则部分透射并叠加散射分量，使像元值既包含地表信息又包含云的信息，这类像元最难处理，因为无论标记为云还是地表都不完全正确。

对下游产品的具体影响：

量化一下量级：一景 30 米分辨率的 Landsat 影像约 3.6 亿个像元，若云覆盖 40% 则有 1.4 亿像元不可用。在做区域统计时，如果不做掩膜，云像元会把平均值拉偏，且偏差方向固定，不会被随机误差抵消。

工程上的第一原则是「先掩膜，再计算」。任何指数计算、分类、统计都应在掩膜之后进行，而不是先算再过滤，因为云像元在计算过程中可能产生异常值（如负的 NDVI、溢出的比值），污染中间结果。

1.1 云覆盖的时空分布规律

云覆盖不是随机的，它有明确的气候与季节规律，理解这些规律才能合理规划数据需求。

区域类型        年均云覆盖    有效观测占比    典型挑战
赤道雨林        80%~90%       10%~20%        常年多云，时序几乎不可用
热带季风        60%~75%       25%~40%        雨季集中，旱季窗口短
温带大陆        45%~60%       40%~55%        季节差异大，冬季雪云混淆
干旱沙漠        10%~25%       75%~90%        云少但亮色地表易误检
高纬地区        60%~75%       25%~40%        极夜与低太阳角，雪覆盖

工程含义很直接：在雨林区域，指望单传感器光学数据做月度监测是不现实的，必须依赖 SAR 或拉长合成周期。而在沙漠区域，问题不是数据量而是误检率，需要更保守的阈值与更强的后处理。

云覆盖还有日内规律。积云在午后发展旺盛，上午的过境时刻（如 Landsat 的 10:30 左右）云量通常低于午后。这也是多数光学卫星选择上午过境的原因之一。

2. 云的物理特性与光谱特征

云由水滴与冰晶组成，在可见光到短波红外波段表现为高反射、低吸收，且各波段反射率接近，因此视觉上呈白色。这一「光谱平坦」特性是云检测的基础判据之一。

关键光谱特征：

• 蓝光与红光高反射：厚云反射率常在 0.7 以上，而多数地物低于 0.4。
• 近红外与短波红外仍高：雪在短波红外（1.6 微米）反射率急剧下降，而云保持高值，这是区分云雪的核心判据。
• 热红外低温：云顶温度低，与地表形成明显温差，厚云可低 20 K 以上。
• 卷云的特殊性：薄卷云在可见光下几乎透明，但在 1.38 微米水汽吸收带呈强反射，这是卷云检测专用通道。

云 / 雪 / 亮色地物判别表（反射率近似值）
目标     蓝(0.48)  红(0.66)  近红外(0.86)  短波红外(1.6)  热红外
厚云      0.75      0.75      0.72          0.65           低温
雪        0.85      0.82      0.80          0.10           低温
沙漠      0.30      0.38      0.45          0.50           高温
城市      0.20      0.22      0.25          0.28           高温
水体      0.06      0.04      0.02          0.01           中温

短波红外是云雪分离的关键：雪在 1.6 微米吸收强烈，云不吸收。没有短波红外通道的传感器（如 Sentinel-2 的 10 米分辨率只有四个波段）云雪分离能力显著下降，这也是 Sentinel-2 的 20 米短波红外波段在云检测中不可替代的原因。

云阴影的成因是云的几何遮挡：太阳光被云阻挡，地面出现与云形状相似但沿太阳方位角偏移的暗区。偏移距离取决于云高与太阳高度角，云越高、太阳越低，偏移越远。这使得阴影检测必须与云检测配对，且需要估计云高。

3. 阈值法与 Fmask 算法

阈值法是最经典的云检测路线，核心是用多个波段的组合判据逐步筛选。Fmask（Function of mask）是应用最广的开源实现，最初为 Landsat 设计，现已扩展到 Sentinel-2。

Fmask 的处理流程分几步：

1. 潜在云层检测（基于亮度和温度）
   - 亮度概率：可见光波段反射率高于阈值
   - 温度概率：热红外亮温低于阈值
   两者组合得到潜在云

2. 云雪分离（用短波红外与归一化雪指数）
   - 短波红外反射率 > 0.03 判为云
   - 否则判为雪

3. 卷云检测（1.38 微米通道）

4. 云阴影检测（几何投影 + 暗像元判据）

5. 形态学膨胀（补偿边缘模糊与配准误差）

import numpy as np

def fmask_lite(blue, red, nir, swir1, bt, nir_swir_ratio=0.8):
    # 简化版潜在云判据，非完整 Fmask 实现
    brightness = (blue + red + nir) / 3.0
    potential_cloud = (brightness > 0.35) & (bt < 300.0) & (nir / (swir1 + 1e-6) > nir_swir_ratio)
    # 云雪分离：短波红外高者为云
    cloud = potential_cloud & (swir1 > 0.03)
    snow = potential_cloud & (swir1 <= 0.03)
    return cloud, snow

阈值法的参数需要按区域标定。同一个 0.35 的亮度阈值，在沙漠区域会大量误检，在热带雨林则可能漏检。生产系统应当按气候区维护多套参数，或用局部自适应阈值替代全局固定值。

Fmask 的 Python 版本（fmask 库）支持 Landsat 4/5/7/8/9 与 Sentinel-2：

python -m fmask --sensor S2 --source ./S2A_L1C.SAFE --output ./mask.tif
python -m fmask --sensor L8 --source LC08_L1TP_118038_20230512 --output ./mask.tif

4. Sen2Cor 与官方掩膜产品

多数业务场景不需要自己实现云检测，直接使用官方掩膜产品是更经济的选择。Sentinel-2 的 Sen2Cor 处理器生成 L2A 产品，其中包含场景分类图（Scene Classification Layer，SCL），每个像元标注为 11 类之一。

SCL 类别定义：

工程上通常把 3、8、9、10、11 一并掩膜，7 按业务容忍度决定。把 8 与 10 也掩掉会损失部分有效像元，但能避免边缘污染，这是保守策略的常见选择。

Landsat 的官方掩膜在 QA_PIXEL 波段里以位标志给出：

import numpy as np

def decode_qa_pixel(qa):
    # Landsat Collection 2 QA_PIXEL 位定义
    fill = (qa & (1 << 0)) != 0
    dilated_cloud = (qa & (1 << 1)) != 0
    cirrus = (qa & (1 << 2)) != 0
    cloud = (qa & (1 << 3)) != 0
    cloud_shadow = (qa & (1 << 4)) != 0
    snow = (qa & (1 << 5)) != 0
    clear = (qa & (1 << 6)) != 0
    water = (qa & (1 << 7)) != 0
    return {"fill": fill, "cloud": cloud, "shadow": cloud_shadow,
            "cirrus": cirrus, "snow": snow, "clear": clear, "water": water}

位运算的顺序容易写错，建议直接用 rio-cogeo 或 xarray 生态里的解码函数，并写单元测试固定住位定义。官方产品的优势是可复现、跨团队一致，劣势是算法版本更新会改变掩膜结果，需要记录版本号。

4.1 质量波段的解读纪律

官方质量波段有一个共同陷阱：位标志是独立的，可能同时置位。例如一个像元可以同时是「云」和「云阴影」，也可能是「云」和「雪」同时为真（算法未完全收敛）。解读时必须定义优先级。

def qa_priority_mask(flags):
    # 按优先级从高到低叠加，后写覆盖前写
    mask = np.zeros(flags["cloud"].shape, dtype=np.uint8)
    mask[flags["water"]] = 1          # 水体
    mask[flags["clear"]] = 2          # 晴空
    mask[flags["snow"]] = 3           # 雪
    mask[flags["shadow"]] = 4         # 云阴影
    mask[flags["cirrus"]] = 5         # 卷云
    mask[flags["cloud"]] = 6          # 云
    mask[flags["fill"]] = 7           # 无数据
    return mask

优先级原则是「最保守的判定优先」：无数据 > 云 > 卷云 > 阴影 > 雪 > 晴空。这样能保证即使标志冲突，也不会把可疑像元当作有效数据使用。

另一个纪律是记录掩膜的派生信息：用了哪个波段、哪套阈值、处理器版本、生成时间。这些信息在半年后排查精度问题时是唯一的线索。建议把掩膜与元数据一起存为带属性的 GeoTIFF，属性里写入处理链版本。

5. 机器学习与云概率方法

阈值法在复杂场景下的精度瓶颈催生了机器学习方法。思路是把云检测建模为逐像元分类或概率回归，用人工标注或已有掩膜作为训练标签。

主流方案对比：

MODIS 的 MOD35 云掩膜产品提供 1 公里分辨率的四档置信度（confident clear、probably clear、probably cloudy、confident cloudy），是很多研究的参考真值来源。它的思路是贝叶斯判别：为每个判据计算云与晴空的似然，组合成后验概率。

import numpy as np

def naive_cloud_probability(features, weights, thresholds):
    # 多判据加权投票的简化实现
    votes = np.zeros_like(features[0], dtype=np.float32)
    total_w = 0.0
    for feat, w, thr in zip(features, weights, thresholds):
        votes += w * (feat > thr).astype(np.float32)
        total_w += w
    return votes / total_w

深度学习方法（CloudNet、U-Net 变体）在 Sentinel-2 上报告的总体精度可达 95% 以上，但存在两个工程问题。其一是跨传感器泛化差：用 Sentinel-2 训练的模型直接用在 Landsat 上精度会掉十几个百分点，因为波段响应函数与分辨率不同。其二是标注数据稀缺：像素级云标注成本极高，多数公开数据集规模在几十到几百景之间，容易过拟合到特定区域。

实用建议是：区域固定、数据量大时训练专用模型；区域多变、数据量小时用官方掩膜或 Fmask 加后处理。混合策略也常见，用官方掩膜做主体，用轻量模型修补边缘。

6. 多时相与时间序列云检测

单时相检测的固有局限是无法区分「云」与「一直很亮的地物」。多时相方法利用时间维信息：如果某像元在前后时相都是暗的，而当前时相突然变亮，则大概率是云。

核心判据是时间序列上的突变检测。以 NDVI 为例，云会使 NDVI 骤降，云阴影会使 NDVI 骤升，两者都与相邻时相形成明显差异。

import numpy as np

def temporal_cloud_flag(series, window=5, k=2.5):
    # series: (t, y, x) 的时间序列，返回突变掩膜
    t = series.shape[0]
    flags = np.zeros_like(series, dtype=bool)
    for i in range(window, t - window):
        ref = np.concatenate([series[i - window:i], series[i + 1:i + 1 + window]])
        med = np.median(ref, axis=0)
        mad = np.median(np.abs(ref - med), axis=0) + 1e-6
        z = np.abs(series[i] - med) / (1.4826 * mad)
        flags[i] = z > k
    return flags

用中位数绝对偏差（MAD）而非标准差，是因为时间序列中云本身是离群值，标准差会被它拉大从而降低灵敏度。1.4826 是使 MAD 与标准差在正态分布下可比的一致性常数。

多时相方法依赖两个前提：一是影像已精确配准，否则像元位置错位会产生假突变；二是已做辐射校正，否则季节性的亮度变化会被误判。这两点在 坐标系统与投影变换 与辐射校正篇中已覆盖。

另一个思路是时序拟合，用谐波模型或双逻辑斯蒂模型拟合物候曲线，把偏离拟合曲线的点判为云。这类方法对季节性明显的区域效果好，对常年多云的热带区域则因有效观测太少而失效。

7. 掩膜后处理与形态学操作

原始掩膜几乎总是有噪声：云边缘的半透明像元被漏检，阴影检测在深色地表上产生零星误检，配准误差导致云与阴影错位。后处理的目标是把掩膜修得「干净且不过度膨胀」。

常用操作按顺序：

import numpy as np
from scipy import ndimage

def postprocess_mask(mask, close_px=3, open_px=3, dilate_px=2):
    m = mask.astype(np.uint8)
    # 闭运算填补云内部空洞
    m = ndimage.binary_closing(m, structure=np.ones((close_px, close_px))).astype(np.uint8)
    # 开运算去除孤立小斑
    m = ndimage.binary_opening(m, structure=np.ones((open_px, open_px))).astype(np.uint8)
    # 膨胀补偿边缘半透明像元与配准误差
    m = ndimage.binary_dilation(m, structure=np.ones((dilate_px, dilate_px))).astype(np.uint8)
    return m.astype(bool)

顺序有讲究。闭运算在前填补云内部因薄云被漏检的孔洞，开运算在中去除孤立噪声，膨胀在后覆盖边缘。如果先膨胀再开运算，噪声也会被放大后难以去除。

膨胀半径是关键参数。经验值是 2 到 3 个像元，对应 Sentinel-2 的 20 到 60 米。设小了边缘半透明像元仍会进入下游，设大了有效像元损失显著。可以用地面验证数据做敏感性分析，找出精度与数据保留率的平衡点。

阴影处理有两条路线。一是独立检测：用几何投影把云按太阳方位角与云高偏移，得到候选阴影区，再用亮度阈值确认。二是配对膨胀：把云掩膜按几何关系平移后再膨胀，简单但精度依赖云高估计。

def shadow_from_cloud(cloud_mask, cloud_height_m, sun_elev_deg, sun_azim_deg, res_m):
    offset = cloud_height_m / np.tan(np.radians(sun_elev_deg)) / res_m
    dx = int(round(offset * np.sin(np.radians(sun_azim_deg + 180))))
    dy = int(round(offset * np.cos(np.radians(sun_azim_deg + 180))))
    return np.roll(np.roll(cloud_mask, dy, axis=0), dx, axis=1)

np.roll 的环绕特性会在影像边界产生错误，实际实现应改用带边界填充的平移。云高未知时可用 2 到 4 公里作为默认值，或从热红外的云顶温度反推。

8. 时序合成与去云重建

掩膜的最终目的是获得可用的时间序列。当单景有效观测不足时，需要跨时相合成，用邻近时相填补被云遮挡的像元。

主流合成策略：

import numpy as np

def quality_mosaic(stack, quality, valid_threshold=0.5):
    # stack: (t, y, x)，quality: (t, y, x) 越大越好
    q = np.where(np.isnan(stack), -np.inf, quality)
    idx = np.argmax(q, axis=0)
    return np.take_along_axis(stack, idx[np.newaxis], axis=0)[0]

时间序列插值要格外小心。线性插值在长云期（连续 30 天以上）会产生与真实物候不符的直线，用这种数据做物候提取会把生长季起止日期算错。稳妥做法是限制插值跨度，超过阈值就标记为缺失而不是强行填补。

遥感领域常用的 Gap-filling 方法有基于谐波的 HANTS、基于时空加权的 NSPI、以及基于邻域相似的 GNSPI。这些方法的共同假设是「地物变化是平滑的」，对火灾、洪水、砍伐这类突变事件不成立，会把突变抹平成渐变。做 变化检测工程实践 时，绝不能使用填补后的数据，必须用原始观测。

SAR 数据不受云影响，是云区时序分析的重要补充。在光学数据长期缺失的区域，用 Sentinel-1 的 C 波段后向散射可以维持监测连续性，方法见 SAR 与 InSAR 处理 。

8.1 合成周期与观测频率的权衡

合成周期决定了每个像元的可用观测数量，进而决定合成结果的可靠性。周期太短则有效观测不足，合成值由个别观测主导；周期太长则时间分辨率下降，错过真实变化。

合成周期与观测数估算（Sentinel-2 A/B 双星，5 天重访）
区域        5 天窗口   10 天窗口   30 天窗口   建议周期
赤道雨林    0.5~1      1~2         3~6        30 天以上
温带大陆    2~3        4~6         12~18      10~15 天
干旱沙漠    4~5        8~10        25~30      5~10 天

经验规则是每个像元在合成窗口内至少有 3 个有效观测，中值合成才有意义。低于 3 个时应拉长周期或退回质量优先策略。

降尺度技巧：先在高频窗口做质量优先合成得到「最优观测」，再在低频窗口做统计合成，兼顾时间分辨率与抗噪能力。这在物候提取中比单一周期效果更好。

需要强调的是，合成产品的元数据里必须记录每个像元的有效观测数与合成方法。没有这两个信息，下游无法判断某个像元值的可信度，也无法复现结果。

9. 精度评估与验证方法

云掩膜的精度评估必须区分不同类别的代价。用总体精度会掩盖问题：如果云只占影像 10%，把所有像元判为晴空的朴素算法也能拿到 90% 的总体精度。

正确的指标组合：

import numpy as np

def mask_metrics(pred, truth):
    tp = int(np.sum(pred & truth))
    fp = int(np.sum(pred & ~truth))
    fn = int(np.sum(~pred & truth))
    tn = int(np.sum(~pred & ~truth))
    precision = tp / (tp + fp + 1e-9)
    recall = tp / (tp + fn + 1e-9)
    iou = tp / (tp + fp + fn + 1e-9)
    f1 = 2 * precision * recall / (precision + recall + 1e-9)
    return {"precision": precision, "recall": recall,
            "iou": iou, "f1": f1, "tp": tp, "fp": fp, "fn": fn, "tn": tn}

关键指标的解读：

• 召回率（Recall）：真实云中被检出的比例。低召回意味着云污染进入下游，是最危险的错误。
• 精确率（Precision）：判为云的像元中确实是云的比例。低精确意味着有效像元被浪费。
• IoU：交并比，对类别不平衡更鲁棒，是云检测论文最常报告的主指标。
• F1：精确率与召回率的调和平均，平衡两者的综合指标。

生产系统的验收标准通常是召回率优先。以 Sentinel-2 为例，Fmask 在厚云上的召回率可达 95% 以上，但在薄卷云上常低于 70%，这是当前技术的主要短板。

验证数据的来源有三类：人工目视解译（最可靠但成本高）、其他掩膜产品交叉比对（快速但有循环论证风险）、地面云观测站（空间尺度不匹配）。实践中常用组合方案，用少量人工标注做基准，用大量交叉比对做趋势监控。

验证还要按地物类型分层。云在雪地、沙漠、亮色城市上最容易误检，这些子集的精度必须单独报告。一个在总体上有 93% IoU 的掩膜，在雪地上的 IoU 可能只有 60%，如果业务涉及积雪区，这个数字才是有意义的。

权衡取舍

• 官方掩膜 vs 自实现：官方可复现、免维护，自实现可针对区域优化，建议默认官方、特殊需求才自研。
• 召回优先 vs 精确优先：多数业务应召回优先，宁多丢像元也不放云进下游，仅在数据稀缺区才放宽。
• 膨胀半径：小则边缘污染，大则有效像元损失，2 到 3 像元是常见折中，需按业务做敏感性分析。
• 阈值法 vs 机器学习：阈值法零训练、可解释，ML 精度高但需标注与跨传感器验证，按数据规模选择。
• 单时相 vs 多时相：多时相精度更高但依赖配准与校正质量，且需要足够的时序长度。
• 插值填补 vs 标记缺失：插值让时序连续但会抹平突变，突变敏感业务必须保留缺失标记。
• 光学 vs SAR 互补：SAR 不受云影响但物理含义不同，混用时需注意后向散射与反射率不可直接比较。

常见坑清单

• 用总体精度评估：现象是掩膜看起来精度很高但下游仍被污染，原因是云占比低导致类别不平衡，规避方法是报告召回率与 IoU。
• 只掩云不掩阴影：现象是时序中 NDVI 出现异常高峰，原因是阴影区被当作正常地表，规避方法是用几何投影配对生成阴影掩膜。
• 阴影误伤水体：现象是水体面积统计忽高忽低，原因是阴影与水体光谱接近，规避方法是加入近红外低反射与形状约束联合判定。
• 雪被误判为云：现象是冬季高纬影像有效像元极少，原因是未使用短波红外分离云雪，规避方法是引入 1.6 微米通道判据。
• 膨胀半径过大：现象是有效像元损失超过 20%，原因是统一用了大半径，规避方法是按业务精度需求做参数扫描。
• 先膨胀后开运算：现象是噪声被放大后无法去除，原因是形态学顺序错误，规避方法是闭、开、膨胀的标准顺序。
• 掩膜与影像配准不一致：现象是云边缘残留亮边，原因是掩膜与影像分辨率或网格不同，规避方法是掩膜重采样到影像网格且用最近邻。
• 用填补数据做变化检测：现象是火灾与砍伐被抹平，原因是插值假设变化平滑，规避方法是变化检测只用原始观测。
• 忽略掩膜版本：现象是同一景不同时间处理结果不一致，原因是官方算法版本更新，规避方法是记录处理器版本号。
• 薄卷云漏检：现象是影像整体偏亮且分类精度下降，原因是卷云在可见光下几乎透明，规避方法是引入 1.38 微米通道或薄云专用算法。

小结

云检测的本质是在漏检与误检之间做有意识的取舍，而这个取舍必须由业务决定而不是由算法默认。理解云与阴影的光谱配对关系、短波红外在云雪分离中的不可替代性、以及召回率优先的评估纪律，是避免踩坑的基础。

工程落地建议分三层：默认使用官方掩膜（Sentinel-2 的 SCL、Landsat 的 QA_PIXEL）作为主体，用形态学后处理修补边缘，用多时相突变检测补齐薄云漏检。只有在区域固定、数据量大、精度要求极高时才自研模型，且必须做跨传感器与跨区域的泛化验证。

下一步建议把掩膜接入 光谱指数计算 的流水线，观察掩膜前后 NDVI 时序的差异作为质量指标；也可以结合 影像分类工程实践 理解云污染如何在分类精度上暴露，最后在 瓦片服务与切片金字塔 中把掩膜作为独立图层发布，供前端交互式开关。