4. 农业遥感影像加载与预处理

作物长势监测依赖的是“可比较、可计算”的影像数据，而不是原始传感器记录值。多光谱遥感影像从卫星或无人机下传到本地时，往往以分波段文件、整景大范围、原始数字量化值（DN值）的形式存在，直接拿来计算NDVI等植被指数，会因波段未对齐、大气与太阳高度差异、像元尺度不一致而得到错误结论。本节围绕一条完整的预处理链展开：加载多光谱影像 → 波段组合 → 辐射定标 → 裁剪研究区 → 重采样统一分辨率，最终得到一套空间范围一致、物理量明确、可直接参与长势分析的多波段栅格。

一、加载多光谱影像：先认清数据组织方式

农业遥感常见数据源包括 Sentinel-2、Landsat 8/9 以及无人机多光谱相机（如蓝、绿、红、红边、近红外五波段）。在 QGIS 中加载前，先要判断数据是“单波段分文件”还是“多波段合成文件”。Sentinel-2 L2A 产品通常每个波段一个 JP2 文件，并附带一个 10 m、20 m、60 m 分辨率说明；Landsat Collection 2 则是每个波段一个 TIF，外加一个元数据 MTL 文本；无人机数据常为单文件多波段 TIF 或分波段 TIF。

加载方式有两种：一是直接把文件拖入图层面板，QGIS 会按文件类型自动识别；二是通过“图层 → 添加图层 → 添加栅格图层”，在对话框中指定源文件。对于分波段数据，建议先建立一个工程，把同一景的所有波段放在同一目录下，命名规范如 T31TCJ_20240612_B04.tif，便于后续批量处理。加载后应检查三件事：波段数是否正确、像元分辨率是否与元数据一致、坐标系是否为 UTM 等投影坐标系。若显示为地理坐标系（经纬度），应在处理前统一重投影，否则后续面积统计和重采样会出现偏差。

二、波段组合：把分散波段拼成一张多光谱栅格

波段组合（Band Stack / Layer Stacking）的目的是把同一景、同一范围、同一分辨率的多个单波段文件合并为一个多波段栅格，使后续辐射定标、指数计算可以一次性作用于整幅影像。QGIS 中可使用“栅格 → 杂项 → 合并”工具，或在处理工具箱中搜索“Merge”，将多个单波段图层按顺序合并。更规范的做法是使用 GDAL 的 gdal_merge.py 或 gdalbuildvrt 先建立虚拟栅格再转实体文件，便于控制波段顺序。

波段顺序必须与后续计算约定一致。以 Sentinel-2 为例，常用顺序为 B02（蓝）、B03（绿）、B04（红）、B08（近红外），对应 NDVI 计算时红波段为第 3 波段、近红外为第 4 波段。若顺序错乱，NDVI 会出现负值或异常高值。组合完成后，在图层属性中查看波段数，并用“识别要素”工具在典型地物（如农田、水体、道路）上读取各波段值，确认数值量级合理——植被在近红外波段应明显高于红波段，这是判断波段是否错位最直接的依据。

三、辐射定标：把 DN 值变成有物理意义的反射率

原始影像记录的是数字量化值（DN），受传感器增益、太阳高度角、大气散射影响，不同日期、不同景之间不可直接比较。辐射定标分两步：第一步是辐射亮度定标，把 DN 转为辐亮度；第二步是大气校正，把辐亮度转为地表反射率。对于作物长势监测，推荐直接使用已做大气校正的 L2A 级产品；若只有 L1C 级数据，则需在 QGIS 中借助 SCP 插件或外部工具（如 Sen2Cor）完成。

以 Landsat Collection 2 Level-1 为例，定标公式为：

Lλ = M_L × Q_cal + A_L
ρ  = (M_ρ × Q_cal + A_ρ) / sin(θ_SE)

其中 M_L、A_L 为辐射定标乘系数与加系数，M_ρ、A_ρ 为反射率定标系数，θ_SE 为太阳高度角，均可从 MTL 元数据文件中读取。在 QGIS 中可用“栅格计算器”逐波段套用公式，也可用 GDAL 的 gdal_calc.py 批量处理。示例参数（合理假设）：某景 B04 的 M_ρ = 2.0000E-05、A_ρ = -0.100000、太阳高度角 55°，则 DN=8000 时反射率约为 (2e-5×8000-0.1)/sin55° ≈ 0.073，处于植被红波段合理区间。定标后应检查值域是否落在 0–1，若出现大量负值或超过 1，说明系数或太阳高度角用错。

四、裁剪与重采样：统一空间范围与像元尺度

整景影像范围远大于研究区，直接计算会浪费算力并引入无关地物。裁剪使用“栅格 → 提取 → 按掩膜图层裁剪栅格”，以研究区矢量边界为掩膜，勾选“将输入图层的范围与掩膜图层对齐”，输出仅保留农田范围。裁剪前应确保矢量与栅格坐标系一致，否则会出现偏移或空结果。

重采样用于解决波段间分辨率不一致的问题。Sentinel-2 的 B08 为 10 m，而 B11、B12 为 20 m，若要把短波红外参与计算，需统一到同一网格。在 QGIS 中使用“栅格 → 投影 → 重采样”，设置目标分辨率 10 m，重采样方法按用途选择：连续型反射率数据用双线性插值，分类数据（如土地覆盖）用最近邻，避免引入不存在的类别。重采样后应核对输出像元大小、行列数与范围，确保与红、近红外波段完全对齐，否则波段运算会报“尺寸不一致”。

五、结果检查与常见问题

完成上述步骤后，用“栅格计算器”计算 NDVI 做一次快速验证：(B08-B04)/(B08+B04)。健康作物NDVI通常在 0.6–0.9，裸土在 0.1–0.2，水体为负值。若整幅影像NDVI接近 0 或出现大面积负值，优先排查三个原因：波段顺序错位、定标系数未套用、重采样后波段未对齐。此外，裁剪后若边缘出现黑边，说明掩膜范围小于影像且未设置 NoData，应在裁剪参数中指定无效值并勾选“裁剪后范围与掩膜一致”。

经过这套流程，输出的是一个范围限定在研究区、像元统一为 10 m、数值为地表反射率的多波段栅格。它可以直接进入下一阶段的植被指数计算与长势分级，也便于在不同时相之间做差值比较——这正是作物长势监测对影像数据“可比、可算、可复现”的基本要求。