📚 课程目录第3章亮温辅助的云顶温度掩膜

3. 亮温辅助的云顶温度掩膜

为什么可见光阈值会把雪和沙地判成云

厚云在可见光波段反射率高，雪面和亮沙地同样高。仅用蓝光反射率阈值时，这三类目标在特征空间里几乎重叠，阈值一放松就误检，一收紧就漏掉薄高云。热红外亮温提供了另一维判据：云顶越高，云顶温度越低。中高云云顶亮温常明显低于同季节地表亮温，而雪面和沙地虽然反射率高，其表面温度接近甚至高于周围地表。把“高反射率”和“低亮温”两个条件联合，就能把高云从亮地表目标中分离出来。

亮温数据的准备与单位统一

热红外波段通常以DN值存储，需要先转为辐亮度再反演亮温。以Landsat 8/9 第10波段（TIRS Band 10）为例，采用官方定标系数：

# 示例参数（Landsat 8/9 TIRS Band 10，具体以元数据为准）
K1 = 774.8853      # W/(m^2·sr·μm)
K2 = 1321.0789     # K
RADIANCE_MULT = 3.3420e-4
RADIANCE_ADD  = 0.10000

def dn_to_brightness_temp(dn):
    """DN -> 辐亮度 -> 亮温(K)"""
    L = RADIANCE_MULT * dn + RADIANCE_ADD          # 辐亮度
    T = K2 / math.log(K1 / L + 1.0)                # 亮温，单位K
    return T

def k_to_celsius(t_k):
    return t_k - 273.15

关键点：亮温必须与反射率影像严格配准、同网格重采样，否则云边缘会出现“反射率高但亮温不低”的错位像素，导致掩膜边缘破碎。重采样建议用最近邻，避免亮温被插值平滑后削弱云顶低温特征。

双条件联合判据的构建

设蓝光反射率为 ρ_blue，亮温为 T_b（单位K）。联合判据写成：

cloud = (rho_blue > R_TH) & (T_b < T_TH)

其中 R_TH 沿用厚云初筛的反射率阈值（示例取 0.25），T_TH 为亮温阈值。T_TH 的确定不能拍脑袋，应基于本景亮温直方图：先统计非云亮地表（雪、沙地、裸土）的亮温分布，取其低分位作为下限参考，再结合已知云区亮温确定。示例场景中，夏季地表亮温约 300–315 K，高云云顶约 255–275 K，取 T_TH = 285 K 可形成明显间隔。

若影像跨越较大纬度或地形高差，单一全局 T_TH 可能不够。可先按高程或纬度分带，再对每带单独定阈值；本例先用全局阈值建立可复现基线。

案例实现：从两幅栅格到一张云掩膜

下面给出可直接运行的完整实现，输入为已配准的蓝光反射率数组和热红外DN数组，输出为0/1云掩膜，并统计被亮温条件剔除的“亮地表误检”像素数，用于验证联合判据是否真正解决了雪/沙地误判。

import numpy as np
import math

# ---------- 亮温反演 ----------
K1, K2 = 774.8853, 1321.0789
RADIANCE_MULT, RADIANCE_ADD = 3.3420e-4, 0.10000

def dn_to_bt(dn):
    L = RADIANCE_MULT * dn + RADIANCE_ADD
    return K2 / np.log(K1 / L + 1.0)

# ---------- 联合判据 ----------
def cloud_mask_bt(rho_blue, dn_tir, r_th=0.25, t_th=285.0):
    bt = dn_to_bt(dn_tir)                       # 亮温(K)
    refl_cloud = rho_blue > r_th                # 可见光高反射
    cold_cloud = bt < t_th                      # 云顶低温
    mask = (refl_cloud & cold_cloud).astype(np.uint8)
    # 仅靠反射率会判为云、但亮温不低的像素 = 疑似亮地表误检
    false_positive_fixed = int(np.sum(refl_cloud & ~cold_cloud))
    return mask, bt, false_positive_fixed

# ---------- 示例数据（教学假设值，非实测） ----------
# 构造 6x6 小场景：高云、雪面、沙地、普通地表
rho_blue = np.array([
    [0.30, 0.32, 0.28, 0.10, 0.12, 0.11],
    [0.31, 0.29, 0.33, 0.11, 0.10, 0.13],
    [0.27, 0.30, 0.31, 0.12, 0.11, 0.10],
    [0.60, 0.62, 0.58, 0.55, 0.57, 0.56],   # 雪/沙地：高反射
    [0.61, 0.59, 0.63, 0.54, 0.56, 0.58],
    [0.58, 0.60, 0.62, 0.57, 0.55, 0.59],
], dtype=float)

# 对应亮温DN（示例）：高云区DN低，雪/沙地DN高
dn_tir = np.array([
    [28000, 28500, 28200, 40000, 40500, 39800],
    [28300, 28100, 28600, 40200, 39900, 40100],
    [28400, 28200, 28000, 40300, 40000, 39700],
    [30000, 30200, 29800, 39500, 39800, 39600],  # 雪/沙地亮温高
    [30100, 29900, 30300, 39700, 39600, 39900],
    [29800, 30000, 30200, 39800, 39700, 39500],
], dtype=float)

mask, bt, fp_fixed = cloud_mask_bt(rho_blue, dn_tir)
print("亮温范围(K):", round(bt.min(),1), "~", round(bt.max(),1))
print("云掩膜:\n", mask)
print("被亮温条件剔除的亮地表误检像素数:", fp_fixed)

运行后应观察到：左上 3×3 高云区亮温约 270–275 K，满足 T_b < 285 K，被保留为云；右下 3×3 雪/沙地区亮温约 300–305 K，虽 ρ_blue > 0.25 但亮温不低，被联合判据剔除。fp_fixed 统计的正是这些“仅靠反射率会误判”的像素，是验证本例目标是否达成的直接指标。

阈值敏感性检查与结果判读

亮温阈值不是固定常数，应做一次小范围扫描，观察云像元数与误检剔除数的变化，选择两者都稳定的区间：

for t in [275, 280, 285, 290, 295]:
    m, _, fp = cloud_mask_bt(rho_blue, dn_tir, r_th=0.25, t_th=t)
    print(f"T_TH={t}K  云像元={int(m.sum())}  剔除误检={fp}")

判读原则：T_TH 过低会漏掉较暖的低云；过高会把冷地表（如高海拔雪面在冬季）重新纳入云。若扫描发现云像元数在某一区间内基本不变，说明该区间内云与地表亮温分离良好，取区间中值作为工作阈值最稳妥。

结果验证与常见问题

云边缘破碎：多为亮温与反射率未严格配准，或重采样用了双线性导致云顶低温被平滑。改用最近邻重采样并检查配准残差。
冬季高海拔雪面被误判为云：雪面亮温可能低于 285 K。此时应结合高程分带降低 T_TH，或引入短波红外反射率辅助区分雪与云。
低云漏检：低云云顶温度接近地表，亮温条件失效。本例定位是高云与亮地表分离，低云需另用其他判据，不应通过抬高 T_TH 强行纳入。
验证方式：抽取云区、雪/沙地区、普通地表各若干样本点，核对掩膜值；同时记录 fp_fixed，确认联合判据确实剔除了仅靠反射率会误判的像素。