第2章遥感图像预处理

第二章：遥感图像预处理——辐射定标与大气校正

各位同学，大家好。今天我们来聊聊遥感图像预处理里最核心的一环——辐射定标与大气校正。说实话，很多初学者拿到影像就直接做分类，结果精度惨不忍睹。我当年也犯过这个错，后来才明白：预处理没做好，后面全是白忙活。

核心观点：遥感图像预处理，说白了就是把传感器记录的“数字”还原成地物真实的“物理量”。这一步做不好，分类精度就是空中楼阁。

2.1 辐射定标——把DN值变成有意义的物理量

辐射定标，就是把影像的DN值（Digital Number）转换成辐射亮度值或反射率。为什么要做？因为DN值只是传感器记录的相对数值，不同传感器、不同时间拍的同一地物，DN值可能差很多。你想想看，这怎么比较？

我个人习惯把辐射定标分成两步：

绝对定标：把DN值转成辐射亮度值（单位：W·m⁻²·sr⁻¹·μm⁻¹）
相对定标：把辐射亮度值转成表观反射率（TOA Reflectance）

代码实现其实很简单，以Landsat 8为例：

import numpy as np
from osgeo import gdal

def radiometric_calibration(dn_image, gain, bias):
    """
    辐射定标：DN -> 辐射亮度
    gain: 增益系数
    bias: 偏置系数
    """
    radiance = dn_image * gain + bias
    return radiance

# 读取影像
ds = gdal.Open('LC08_B4.tif')
dn = ds.ReadAsArray()

# Landsat 8 波段4的增益和偏置（示例值）
gain = 0.0003342
bias = 0.1

rad = radiometric_calibration(dn, gain, bias)
print(f'辐射定标完成，辐射亮度范围：{rad.min():.2f} - {rad.max():.2f}')

小技巧：我在项目中遇到过，有些影像的元数据文件（MTL.txt）里直接提供了增益和偏置。但要注意，不同版本的Landsat数据，系数可能不一样。我建议每次都从元数据里读，别偷懒写死。

2.2 大气校正——去掉大气的“滤镜”

大气校正，就是把表观反射率转成地表真实反射率。为什么要做？因为大气中的水汽、气溶胶、臭氧等会吸收和散射电磁波，让影像看起来“雾蒙蒙”的。说白了，就是去掉大气这层“滤镜”。

常用的方法有三种：6S、FLAASH、Sen2Cor。我一个个说。

2.2.1 6S模型——学术界的“老大哥”

6S（Second Simulation of the Satellite Signal in the Solar Spectrum）是法国里尔大学开发的。它基于辐射传输方程，理论上最精确。但缺点也很明显——参数太多，配置复杂。

我曾经在项目中用6S做高光谱影像的大气校正，光是输入参数就调了三天。嗯，这里要注意：6S需要输入气溶胶光学厚度、水汽含量、臭氧浓度等参数，这些数据通常要从MODIS或气象站获取。

# 6S的输入文件示例（简化版）
# 几何参数
solar_zenith = 30.0      # 太阳天顶角
sensor_zenith = 0.0      # 传感器天顶角
relative_azimuth = 0.0   # 相对方位角

# 大气参数
aerosol_model = 1        # 1=大陆型, 2=海洋型
aerosol_optical_depth = 0.2  # 气溶胶光学厚度
water_vapor = 2.5        # 水汽含量(g/cm²)
ozone = 0.3              # 臭氧浓度(cm·atm)

2.2.2 FLAASH——ENVI里的“明星工具”

FLAASH（Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes）是ENVI软件里的大气校正模块。它基于MODTRAN辐射传输模型，但做了加速优化。我个人觉得，FLAASH是平衡精度和易用性的最佳选择。

FLAASH的输入参数相对简单：

影像中心经纬度
传感器类型（Landsat、SPOT、WorldView等）
飞行日期和时间
气溶胶模型（城市、乡村、海洋等）
能见度（或气溶胶光学厚度）

避坑指南：我曾经在FLAASH里选了“城市”气溶胶模型，结果校正后的影像偏蓝。后来才发现，那个区域其实是郊区，应该选“乡村”模型。所以，气溶胶模型的选择一定要结合实地情况，别想当然。

2.2.3 Sen2Cor——Sentinel-2的“专属搭档”

Sen2Cor是欧空局为Sentinel-2影像开发的大气校正工具。它最大的优点是自动化程度高，几乎不需要手动调参。你只需要输入影像文件，它就能自动完成大气校正、云检测、地形校正等。

我建议，如果你主要用Sentinel-2数据，直接上Sen2Cor。省时省力，效果还稳定。

# Sen2Cor命令行调用示例
# 安装后，在终端运行：
# L2A_Process --resolution 10 input.SAFE

# Python调用（通过subprocess）
import subprocess

def run_sen2cor(input_safe, resolution=10):
    cmd = f'L2A_Process --resolution {resolution} {input_safe}'
    subprocess.run(cmd, shell=True)
    print(f'Sen2Cor处理完成：{input_safe}')

# 使用示例
run_sen2cor('S2A_MSIL1C_20230101.SAFE', resolution=10)

2.3 三种方法的对比

为了让你看得更清楚，我整理了一个对比表：

方法	精度	易用性	适用数据	我的推荐
6S	高	低（参数多）	多光谱、高光谱	学术研究，不赶时间
FLAASH	中高	中（需调参）	多光谱（Landsat、SPOT等）	日常项目，精度和效率平衡
Sen2Cor	中	高（全自动）	Sentinel-2	批量处理，省心省力

2.4 几何校正——让影像“对号入座”

几何校正，就是把影像的像素坐标映射到地理坐标。说白了，就是让影像里的每个点，都能在地球上找到准确的位置。

正射校正是几何校正的一种，它考虑了地形起伏的影响。我建议，在山区或丘陵地带，一定要做正射校正，否则分类结果会“飘”。

from osgeo import gdal, osr

def orthorectification(input_tif, output_tif, dem_file):
    """
    正射校正（使用GDAL的Warp函数）
    input_tif: 原始影像
    output_tif: 校正后影像
    dem_file: DEM文件
    """
    # 设置正射校正参数
    warp_options = gdal.WarpOptions(
        format='GTiff',
        dstSRS='EPSG:4326',  # 输出坐标系
        rpc=True,            # 使用RPC文件
        dem=dem_file,        # DEM文件
        outputType=gdal.GDT_Float32
    )
    
    # 执行正射校正
    gdal.Warp(output_tif, input_tif, options=warp_options)
    print(f'正射校正完成：{output_tif}')

# 使用示例
orthorectification('LC08_B4.tif', 'LC08_B4_ortho.tif', 'srtm_30m.tif')

经验之谈：我在项目中遇到过，有些影像的RPC文件缺失或损坏。这时候，可以用地面控制点（GCP）来做几何校正。但GCP的精度直接影响校正结果，我建议至少选20个均匀分布的控制点。

2.5 图像融合与镶嵌

图像融合，就是把高空间分辨率（但低光谱分辨率）的影像和低空间分辨率（但高光谱分辨率）的影像融合，得到“两全其美”的结果。常用的方法有：

Brovey变换：简单快速，但容易产生色彩失真
Gram-Schmidt：精度高，是ENVI的默认方法
PC光谱锐化：基于主成分分析，适合多光谱影像

图像镶嵌，就是把多幅影像拼接成一幅。这里要注意：

色调一致性：不同影像的色调可能不同，需要做匀色处理
接缝线处理：用羽化或加权平均消除接缝
重叠区域：一般取平均值或最大值

2.6 图像裁剪与掩膜

图像裁剪，就是按研究区范围切出需要的部分。掩膜，则是用矢量或栅格数据“遮住”不需要的区域。

我常用的方法是：先用矢量边界裁剪影像，再用掩膜去除云、水体等干扰区域。

def clip_image_by_shapefile(input_tif, shapefile, output_tif):
    """
    用矢量边界裁剪影像
    """
    from osgeo import gdal
    
    # 打开矢量文件
    ds = gdal.OpenEx(shapefile)
    layer = ds.GetLayer()
    
    # 裁剪
    warp_options = gdal.WarpOptions(
        cutlineDSName=shapefile,
        cropToCutline=True,
        dstNodata=0
    )
    
    gdal.Warp(output_tif, input_tif, options=warp_options)
    print(f'裁剪完成：{output_tif}')

# 使用示例
clip_image_by_shapefile('LC08_ortho.tif', 'study_area.shp', 'LC08_clipped.tif')

注意：裁剪前一定要检查矢量和影像的坐标系是否一致。我曾经因为坐标系不匹配，裁剪出来的影像全是黑的。嗯，这个坑我踩过，你别再踩了。

2.7 本章知识体系

为了让你对本章内容有个整体把握，我画了一张流程图：

这张图把整个预处理流程串起来了。你跟着这个流程走，基本不会出大问题。

好了，第二章的内容就到这里。预处理是遥感分类的基石，花再多时间都值得。下一章，我们会讲几何校正的细节和图像融合的实战技巧。到时候见。

📚 课程目录

第二章：遥感图像预处理——辐射定标与大气校正

2.1 辐射定标——把DN值变成有意义的物理量

2.2 大气校正——去掉大气的“滤镜”

2.2.1 6S模型——学术界的“老大哥”

2.2.2 FLAASH——ENVI里的“明星工具”

2.2.3 Sen2Cor——Sentinel-2的“专属搭档”

2.3 三种方法的对比

2.4 几何校正——让影像“对号入座”

2.5 图像融合与镶嵌

2.6 图像裁剪与掩膜

2.7 本章知识体系