remote-sensing-analysis

遥感分析

用途：将原始地球观测影像转化为有依据的分析产品。这里的失误往往不易察觉——把未经校正的 DN 值当作反射率、把云误计为地表覆盖变化、在错误的波段上计算指数——因此，这项技能会优先做好各项检查。

数据访问（优先使用 STAC）

通过 STAC API 搜索，而不是使用各个数据提供方的门户；这样工作流统一且可脚本化：

import pystac_client
import odc.stac

catalog = pystac_client.Client.open("https://earth-search.aws.element84.com/v1")
items = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=[27.0, 38.3, 27.4, 38.6],
    datetime="2025-05-01/2025-09-30",
    query={"eo:cloud_cover": {"lt": 20}},
).item_collection()
ds = odc.stac.load(items, bands=["red", "nir", "scl"], resolution=10, chunks={})

主要数据集：sentinel-2-l2a（10 m 光学影像，地表反射率）、landsat-c2-l2（30 m，始于 1982 年）、sentinel-1-grd（SAR，不受天气影响）。Microsoft Planetary Computer 镜像了大多数数据集（需要通过 planetary_computer 签名）。对于洲际／全球范围或跨数十年的影像堆栈，改用 google-earth-engine，而不是下载数据。记录数据集、项目 ID 和搜索参数，以确保结果可复现。

处理级别规范

级别	含义	可直接用于分析？
L1C / L1TP	大气层顶（TOA）	指数基本可用；跨日期比较有风险
L2A / L2SP	地表反射率（BOA）	可以——默认选择
GRD（SAR）	检测振幅	需要地形校正和斑点滤波

务必说明所用级别。合成影像或时间序列中，切勿混用 TOA 和 BOA 场景。Landsat Collection 2 L2 需要应用其比例因子（reflectance = DN * 0.0000275 - 0.2）。

Sentinel-2 基线不连续——上述级别检查无法发现

自 2022 年 1 月 25 日起应用的处理基线 04.00，为 L2A 数字值添加了一个常量 BOA_ADD_OFFSET（当前为 −1000），以便对负地表反射率进行编码。该日期前后的两个场景都属于 L2A：上述级别检查看不出问题，但它们的 DN 值相差 1000。对它们求差会产生系统性的反射率偏移，看起来像真实变化，并且会通过之后应用的所有掩膜、阈值和精度报告。

• 从每个产品的元数据中读取 BOA_ADD_OFFSET 和 QUANTIFICATION_VALUE，不要将 −1000 和 10000 写死；这两个值按波段设定，且基线过去曾有变化。
• 使用 reflectance = (DN + BOA_ADD_OFFSET) / QUANTIFICATION_VALUE 进行转换。
• 在清单中记录每个场景的处理基线，而不只是产品级别。仅说明两个场景均为 L2A 并不足够。
• **不要重复校正。**经过协调的数据集——例如 Earth Engine 的 COPERNICUS/S2_SR_HARMONIZED 以及若干商业镜像——已将基线更新后的数据移回 2022 年前的范围。再次应用偏移量会反转误差，而不是消除误差。
• 如果任一场景的基线未记录，这种比较就缺乏可靠依据。应如实说明，而不是假定其属于 2022 年前或之后的数据。

云和质量掩膜——先于其他所有操作

• Sentinel-2：使用 SCL 波段进行掩膜（剔除类别 3 云影、8-9 云、10 卷云、11 积雪——谨慎保留 4 植被、5 裸地、6 水体、7 未分类）。
• Landsat C2：解码 QA_PIXEL 位字段（云、云影、卷云位）。
• 报告每个场景掩膜后有效像元的百分比；有效比例低于约 60% 的场景通常应予以排除。
• 对于无缝产品，应构建覆盖一个季节的中值合成影像，而不是挑选单个场景。

光谱指数

使用地表反射率计算，避免除以零，并明确标注波段——不同传感器的波段编号不同（Sentinel-2 的 NIR 为 B8，Landsat 8/9 的 NIR 为 B5）：

import numpy as np
import xarray as xr

def normalized_diff(a: xr.DataArray, b: xr.DataArray) -> xr.DataArray:
    """(a - b) / (a + b) with zero-denominator protection."""
    return xr.where(a + b == 0, np.nan, (a - b) / (a + b))

ndvi = normalized_diff(ds.nir, ds.red)     # vegetation
ndwi = normalized_diff(ds.green, ds.nir)   # open water (McFeeters)
ndbi = normalized_diff(ds.swir16, ds.nir)  # built-up

解读注意事项：NDVI 阈值取决于场景和季节；未经直方图检查，切勿硬编码“NDVI > 0.3 = 植被”。水体会干扰 NDBI；阴影在 NDWI 中会呈现为水体——应交叉核对各项指数，并与真彩色影像对照。

分类工作流

1. 定义图例，类别之间互斥且可通过影像区分。
2. 在场景中空间分散地采集训练样本；将其记录为带版本号的矢量文件。
3. 特征：波段 + 指数 + 纹理（GLCM）+ 多时相统计量（如有多期数据）。对于深度学习路线，交由 geo-deep-learning 处理。
4. 使用空间独立的测试集进行验证（参见 ml-experiment-standards → references/spatial-cv-protocol.md），并 报告各类别的 F1/IoU 及混淆矩阵——仅报告总体准确率会掩盖稀有类别的分类失败。
5. 绘制错误分布图：误分类的空间分布图可以揭示全局指标掩盖的系统性问题（地形阴影、城市/裸地混淆）。

SAR 注意事项（Sentinel-1）

预处理：轨道文件 → 热噪声去除 → 定标（σ⁰）→ 地形校正（使用 DEM 进行 Range-Doppler 校正）→ 斑点滤波（Lee/Refined Lee）→ 转换为 dB。统计时使用 dB 值；VV/VH 比值是区分水体和植被的有效特征。SAR 可以穿透云层——洪水制图和连续监测应优先使用 SAR。

陷阱检查清单

• 比较不同日期的场景时，未进行一致的大气校正。
• 混用 2022-01-25 基线变更前后的 Sentinel-2 场景，却未应用 BOA_ADD_OFFSET；或者对已协调一致的数据集合再次应用该参数。
• 忽略 Sentinel-2 波段混合带来的 20 m→10 m 分辨率差异（B11/B12 的原生分辨率为 20 m）。
• 在未应用比例因子的整数 DN 上计算指数 → 得到毫无意义的取值范围。
• 对 SAR 数据以线性单位制作中位数合成（统计时应使用 dB）。
• 训练像元和测试像元来自同一田块/多边形 → 准确率发生数据泄漏。
• 重投影后忘记应用 nodata 掩膜（边缘变成零 → 产生虚假的土地覆盖）。

执行约定

• **工作流：**定义现象和尺度；选择传感器、产品级别和日期；统一定标、掩膜、CRS 和分辨率；派生特征；分析；进行空间验证；发布溯源信息。
• **决策规则：**影像准备和经典分析使用此技能；明确的时序差分使用变化检测；神经网络训练使用深度学习；归档级执行使用 Earth Engine。
• **验证协议：**检查掩膜和有效像元数量；确认各场景的比例因子、偏移量和处理基线；确认波段分辨率；叠加输出结果；使用空间独立的验证数据；绘制错误分布图；测试季节或传感器敏感性。
• **失效模式：**若存在云或阴影泄漏、处理级别不可比、处理基线未记录或不匹配、重采样伪影、标签泄漏、nodata 污染，或结论超出传感器分辨率所能支持的范围，则应拒绝结果。
• **交付物：**可用于分析的影像或特征、处理清单、掩膜、派生产品、验证指标和错误分布图、可复现代码及局限性。
• **来源时效性：**执行时查阅权威来源登记表，了解产品、定标和目录变更。