辐射校正学习

辐射校正（Radiometric Correction）是遥感图像处理中的关键步骤，旨在消除传感器、大气、地形等因素对辐射值的影响，使图像数据更真实反映地物辐射特性。主要分为绝对辐射校正和相对辐射校正两类。

绝对辐射校正

将图像的数字值（DN）转换为物理辐射量（如辐射亮度或反射率），需结合传感器参数和大气模型。

典型方法：

1. 辐射定标（Radiometric Calibration）
   通过传感器提供的增益（Gain）和偏置（Bias）参数转换DN值：
   $$ L_{\lambda} = Gain \times DN + Bias $$
   其中，$L_{\lambda}$为辐射亮度（W·m⁻²·sr⁻¹·μm⁻¹）

2. 大气校正（Atmospheric Correction）
   使用模型（如6S、FLAASH、MODTRAN）去除大气散射、吸收影响，获取地表反射率。

相对辐射校正

对同一区域的多幅图像或多时相图像进行辐射归一化，消除非地物变化引起的差异。

常用方法：

1. 直方图匹配（Histogram Matching）
   将目标图像的直方图调整到参考图像的分布。
2. 伪不变特征法（PIFs）
   选择稳定的地物点（如道路、裸土）作为基准，进行线性或非线性拟合。

地形辐射校正

针对山区影像，消除地形阴影和太阳入射角差异的影响。

经典模型：

• C模型（Cosine Correction）
  $$ L_{H} = \frac{L}{cos(\theta)} $$
  $\theta$为太阳入射角，$L_{H}$为水平面辐射亮度

• Minnaert校正
  适用于非朗伯体表面，引入经验参数$k$：
  $$ L_{H} = L \cdot \left(\frac{cos \theta_{e}}{cos \theta_{i}}\right)^k $$
  $\theta_{e}$为太阳天顶角，$\theta_{i}$为坡度角

  import numpy as np
  from osgeo import gdal
  
  def radiance_calibration(input_path, output_path, gain, bias):
      """
      辐射定标：将DN值转换为辐射亮度值
      :param input_path: 输入影像路径
      :param output_path: 输出影像路径
      :param gain: 增益系数
      :param bias: 偏移系数
      """
      dataset = gdal.Open(input_path)
      band = dataset.GetRasterBand(1)
      dn = band.ReadAsArray()
      
      # 辐射定标公式: L = gain * DN + bias
      radiance = gain * dn.astype(np.float32) + bias
      
      driver = gdal.GetDriverByName('GTiff')
      out_dataset = driver.Create(output_path, 
                                dataset.RasterXSize, 
                                dataset.RasterYSize, 
                                1, 
                                gdal.GDT_Float32)
      out_dataset.GetRasterBand(1).WriteArray(radiance)
      out_dataset.SetGeoTransform(dataset.GetGeoTransform())
      out_dataset.SetProjection(dataset.GetProjection())
      out_dataset.FlushCache()
      dataset = None
      out_dataset = None
  
  def atmospheric_correction(input_path, output_path, Lp, tau, Ed):
      """
      大气校正：将辐射亮度值转换为地表反射率
      :param input_path: 输入辐射亮度影像路径
      :param output_path: 输出反射率影像路径
      :param Lp: 路径辐射
      :param tau: 大气透射率
      :param Ed: 下行辐照度
      """
      dataset = gdal.Open(input_path)
      band = dataset.GetRasterBand(1)
      L = band.ReadAsArray()
      
      # 大气校正公式: ρ = π * (L - Lp) / (tau * Ed)
      reflectance = np.pi * (L - Lp) / (tau * Ed)
      reflectance = np.clip(reflectance, 0, 1)  # 限制反射率在0-1之间
      
      driver = gdal.GetDriverByName('GTiff')
      out_dataset = driver.Create(output_path, 
                                dataset.RasterXSize, 
                                dataset.RasterYSize, 
                                1, 
                                gdal.GDT_Float32)
      out_dataset.GetRasterBand(1).WriteArray(reflectance)
      out_dataset.SetGeoTransform(dataset.GetGeoTransform())
      out_dataset.SetProjection(dataset.GetProjection())
      out_dataset.FlushCache()
      dataset = None
      out_dataset = None