2. Sentinel-2数据获取与预处理

铁染与羟基蚀变识别依赖短波红外(SWIR)与可见光—近红外波段的组合比值,因此数据源的选择和预处理质量直接决定后续蚀变信息提取的可信度。Sentinel-2作为当前免费开放、重访周期约5天、空间分辨率最高10 m的多光谱数据,是区域尺度蚀变填图的常用选择。本节围绕“拿到一份可直接用于蚀变分析的L2A产品”这一目标,讲清数据获取、波段与分辨率认知、产品级别选择以及元数据检查四个关键环节。

一、为什么优先选择L2A级产品

Sentinel-2产品按处理级别分为L0、L1C和L2A。L1C是经过几何精校正的大气顶(TOA)反射率产品,未做大气校正;L2A是经过Sen2Cor大气校正后的地表反射率(BOA)产品,并附带场景分类层(SCL)。蚀变识别所用的波段比值(如铁染常用的B4/B2、B11/B8A,羟基常用的B11/B12)本质上比较的是地表反射特征,若使用TOA反射率,大气散射与气溶胶的影响会随波段波长不同而变化,导致比值失真,尤其在高气溶胶或薄云条件下误差显著。因此,除非有特殊研究需要,蚀变识别应直接选用L2A产品,避免自行从L1C做大气校正带来的参数不确定性。

需要注意的是,L2A产品在2017年之后的全球覆盖已较为完整,但早期或部分区域可能缺失。若目标区仅有L1C,则需用Sen2Cor自行处理,此时应记录气溶胶反演模式、水汽估计方式等参数,以便复现。

二、波段与分辨率认知

Sentinel-2 MSI共有13个波段,空间分辨率分为10 m、20 m和60 m三档。蚀变识别最关心的是20 m的SWIR波段(B11约1610 nm、B12约2190 nm)以及10 m的可见光—近红外波段(B2蓝、B3绿、B4红、B8宽近红外)和20 m的窄近红外B8A。铁染信息主要落在B2、B4与B11、B12的差异上,羟基信息则集中在B11与B12的吸收特征。由于不同波段分辨率不一致,做波段运算前必须统一重采样到同一网格,通常将10 m波段重采样到20 m,或全部重采样到20 m作为工作分辨率,避免因分辨率混用造成像元错位。

下表给出蚀变识别常用波段及其典型用途,便于在下载和后续处理时快速对照。

波段中心波长(约)原生分辨率蚀变识别中的典型用途
B2 蓝490 nm10 m铁染比值分母,反映铁氧化物在蓝端的吸收
B4 红665 nm10 m铁染比值分子,与B2组合增强铁染信号
B8 宽近红外842 nm10 m植被抑制、亮度参考
B8A 窄近红外865 nm20 m与SWIR组合,减少水汽影响
B11 SWIR11610 nm20 m羟基吸收特征核心波段,铁染辅助
B12 SWIR22190 nm20 m羟基与碳酸盐识别,黏土矿物区分

三、从欧空局获取L2A数据

常用获取渠道有Copernicus Data Space Ecosystem(CDSE)浏览器、Copernicus Open Access Hub以及通过API批量下载。以CDSE为例,操作流程为:注册账号并登录,在搜索面板中设定感兴趣区(可上传矢量或手绘多边形)、时间范围、云量阈值(蚀变识别建议初筛云量低于10%~20%),产品级别选择L2A,提交搜索后按云量和日期排序,选择覆盖目标区且云量最低的景,加入购物车后选择下载。单景L2A产品通常为数百MB至1 GB量级的压缩包,内含各波段JPEG2000文件、SCL分类层、元数据XML以及质量指示文件。

若需批量获取,可使用官方提供的检索接口,例如通过Python的sentinelhub或直接调用OData接口按产品标识下载。示例命令仅示意检索思路,实际使用时需替换为本人账号与目标区参数:

# 示例:按产品名称检索并下载(参数为示意,需替换)
# 产品命名示例:S2B_MSIL2A_20240315T030529_N0510_R075_T50TMK_20240315T045210
# 其中 T50TMK 为UTM瓦片编号,日期与处理基线可从名称直接读取

四、元数据检查:下载后第一件事

下载完成后不要急于做波段运算,应先检查元数据,确认产品可用。重点检查项包括:产品级别是否为L2A、处理基线(Processing Baseline)版本、云量与场景分类统计、太阳天顶角与方位角、投影与瓦片编号、以及是否存在异常的质量指示。处理基线4.00之后,L2A产品引入了偏移量(offset)机制,部分波段的数字量化值需减去1000再除以10000才得到反射率,若忽略这一点,反射率会出现系统性偏差,比值结果随之失真。因此,读取元数据中的BOA_ADD_OFFSET与QUANTIFICATION_VALUE是必做步骤。

场景分类层SCL同样关键。SCL将像元分为植被、裸土、水体、云、云影、雪等类别,蚀变识别前应利用SCL掩膜掉云、云影、水体和雪,只保留裸土与稀疏植被区域参与运算,否则云边缘和阴影会制造大量假异常。太阳几何参数则用于后续地形校正或光照归一化,若目标区地形起伏较大,应记录并考虑引入DEM做地形校正。

检查顺序建议:先看产品级别与处理基线,再看偏移量参数,然后统计SCL各类别占比,最后核对太阳几何与投影信息。任何一项异常都应记录在案,必要时更换日期或瓦片。

五、预处理流程的整体关系

下图把本节涉及的获取、检查与预处理环节组织为一条数据流,重点表现各环节的输入输出与判断分支,帮助理解“什么样的产品才能进入蚀变运算”。

Sentinel-2 L2A 获取与预处理数据流 CDSE 检索 AOI + 时间 + 云量≤20% 产品级别选 L2A L2A 产品包 B2/B4/B8/B8A/B11/B12 SCL + 元数据 XML 元数据检查 基线 / 偏移量 / 太阳角 SCL 类别占比 产品可用? 级别/偏移/云量 更换日期或瓦片 重新检索下载 反射率换算 (DN − offset) / 10000 SCL 掩膜 剔除云 / 云影 / 水体 / 雪 统一网格重采样 10 m → 20 m 工作分辨率 进入蚀变比值运算 B4/B2、B11/B8A、B11/B12 否 是

图中“产品可用?”这一判断节点是整条链路的关口:级别不对、偏移量未处理或云量过高,都应回到检索环节更换数据,而不是带着问题继续运算。通过检查后,先做反射率换算,再用SCL掩膜,最后统一重采样到20 m工作网格,输出的多波段反射率立方体即可进入铁染与羟基的比值运算。

常见错误:直接读取JPEG2000的DN值做比值,未减去偏移量也未除以量化值;或忽略SCL掩膜,把云边缘当作蚀变异常。这两类错误在结果图上往往表现为大面积连片“异常”,与真实矿化蚀变的空间形态明显不符,应在预处理阶段就排除。

完成上述步骤后,建议保留一份预处理记录,写明产品标识、处理基线、偏移量、SCL各类别占比、重采样方法与工作分辨率。这份记录既是结果可复现的依据,也是后续对比不同时相数据时判断差异来源的基础。