数据预处理：辐射定标与大气校正、几何校正与配准、云掩膜与质量控制

各位同学，咱们今天聊点实在的。遥感数据拿到手，看着挺漂亮，但直接扔进模型里？那基本等于自杀。我刚开始做遥感那会儿，就吃过这个亏——拿了一景Landsat 8，没做预处理就硬算NDVI，结果出来的趋势图跟心电图似的，完全没法解释。后来才明白，数据预处理这步，省不得。

说白了，预处理就是给原始数据“洗澡”。把传感器误差、大气干扰、几何变形这些脏东西洗掉，才能看到地物的真实面貌。今天咱们就掰开揉碎了讲：辐射定标、大气校正、几何校正、云掩膜，还有质量控制。

1. 辐射定标：从DN值到物理量

辐射定标是预处理的第一步，也是最基础的一步。你想想看，传感器记录的是DN值（Digital Number），一个无量纲的整数。但我们要做时间序列分析，需要的是物理意义明确的辐射亮度或者反射率。DN值会随着传感器增益、偏移量变化，不同时相的数据如果不统一，趋势分析就是扯淡。

我个人习惯用绝对定标系数来做。Landsat数据在头文件里会给出增益（Gain）和偏移（Bias），公式很简单：

# 辐射定标：DN → 辐射亮度
import numpy as np

def radiance_calibration(dn, gain, bias):
    """
    DN值转辐射亮度
    gain: 增益 (头文件中获取)
    bias: 偏移 (头文件中获取)
    """
    radiance = dn * gain + bias
    return radiance

# 举个例子，Landsat 8 OLI波段4
dn_array = np.array([[120, 135], [98, 145]])
gain = 0.0003342  # 单位：W/(m²·sr·μm)
bias = -0.1
rad = radiance_calibration(dn_array, gain, bias)
print(rad)

定标完得到辐射亮度，再转表观反射率（TOA Reflectance）。这一步要考虑太阳天顶角和日地距离。嗯，这里要注意：Landsat 8之后的传感器，头文件里直接给了反射率定标系数，省事不少。

2. 大气校正：把雾霾去掉

辐射定标只是第一步。大气中的水汽、气溶胶、臭氧会散射和吸收电磁波，导致地表反射率被“污染”。我做过一个对比实验：同一景影像，不做大气校正的NDVI最大值0.6，做完之后0.85。差了将近30%！你说这能直接分析趋势吗？

常用的方法有几种：

• 6S模型：物理模型，精度高，但需要大气参数（气溶胶光学厚度、水汽含量等）。我一般在MODIS数据里找这些参数。
• FLAASH：ENVI里集成的模块，基于MODTRAN，适合高光谱和多光谱数据。
• DOS（暗目标法）：不需要外部参数，假设影像中存在暗目标（如清洁水体），其反射率接近0。简单粗暴，适合快速处理。

我个人在时间序列分析中，更推荐6S模型或者FLAASH。为什么？因为时间序列需要一致性，DOS方法对不同时相的校正结果可能不稳定。我曾经用DOS处理过一组30年的Landsat数据，结果发现1980年代的影像校正后偏暗，2000年后的偏亮——这就是方法本身引入的误差。

3. 几何校正与配准：让像素对齐

时间序列分析的核心是“同一个位置，不同时间”。如果像素没对齐，分析出来的趋势就是空间误差的叠加。Landsat Level-1数据已经做过系统几何校正，但精度一般在30米左右。对于时间序列分析，这还不够。

我建议做亚像素级配准。具体来说：

1. 选一个基准影像（比如第一景，或者质量最好的那一景）
2. 用自动匹配算法（如相位相关、SIFT特征点）找到控制点
3. 用多项式模型（一般2阶就够了）做重采样
4. 检查RMSE（均方根误差），控制在0.5个像素以内

# 用GDAL做几何配准的伪代码
from osgeo import gdal

def image_registration(reference_path, target_path, output_path):
    """
    影像配准：将target影像配准到reference影像
    """
    ref_ds = gdal.Open(reference_path)
    target_ds = gdal.Open(target_path)
    
    # 获取参考影像的投影和地理变换
    ref_proj = ref_ds.GetProjection()
    ref_geotrans = ref_ds.GetGeoTransform()
    
    # 创建输出文件
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(output_path, 
                           ref_ds.RasterXSize, ref_ds.RasterYSize,
                           target_ds.RasterCount, gdal.GDT_Float32)
    out_ds.SetProjection(ref_proj)
    out_ds.SetGeoTransform(ref_geotrans)
    
    # 重采样（推荐使用三次卷积或兰索斯）
    gdal.ReprojectImage(target_ds, out_ds, 
                        target_ds.GetProjection(), ref_proj,
                        gdal.GRA_Lanczos)
    return output_path

4. 云掩膜与质量控制：剔除坏数据

遥感影像最怕什么？云！云覆盖的区域，地表信息完全丢失。时间序列分析中，如果某个时相有云，那个像素的值就是无效的。不处理的话，趋势会被严重扭曲。

Landsat数据自带的QA波段（Quality Assessment）提供了云、云阴影、雪等掩膜信息。我一般这样处理：

• 读取QA波段，提取云和云阴影的比特位
• 生成二值掩膜（1表示好像素，0表示坏像素）
• 对坏像素做插值或直接剔除

# Landsat 8 QA波段云掩膜提取
def cloud_mask_from_qa(qa_band):
    """
    从QA波段提取云掩膜
    QA波段比特位定义：
    Bit 3: 云
    Bit 4: 云阴影
    Bit 5: 雪
    """
    # 云掩膜：Bit 3 = 1 表示云
    cloud_mask = (qa_band & (1 << 3)) != 0
    # 云阴影：Bit 4 = 1
    shadow_mask = (qa_band & (1 << 4)) != 0
    # 合并
    bad_mask = cloud_mask | shadow_mask
    return ~bad_mask  # 返回好像素掩膜

有了掩膜之后，时间序列中缺失的像素怎么办？我个人习惯用时间序列插值。比如用前后时相的中值、或者用Savitzky-Golay滤波来填补。但要注意：如果连续多个时相都有云，插值结果可能不可靠。这时候我宁可舍弃这段数据，也不强行填补。

5. 质量控制：数据预处理的自检清单

预处理做完之后，别急着跑模型。先做一轮质量控制。我给自己定了个检查清单：

检查项	方法	合格标准
辐射定标一致性	对比不同时相的同质区域（如沙漠、水体）	反射率差异 < 5%
大气校正效果	检查植被光谱曲线是否出现“红边”特征	红边位置在700-750nm
几何配准精度	计算控制点RMSE	RMSE < 0.5像素
云掩膜效果	目视检查云覆盖区域是否被正确标记	误检率 < 5%
时间序列完整性	统计每个像素的有效观测次数	有效观测 > 总时相的70%

做完这些检查，数据才算真正“干净”了。我记得有一次，一个学生拿着预处理完的数据来找我，说趋势分析结果很奇怪。我一查，发现他大气校正的参数设置错了，把温带地区设成了热带大气模式。改过来之后，结果就正常了。所以，预处理这步，真的不能马虎。

好了，数据预处理就讲到这里。下一章咱们开始聊时间序列重建和趋势分析方法。到时候见。

---

公众号:蓝海资料掘金营，微信deep3321