6S大气校正原理、实现方法

6S大气校正新方法

本文实现目标：利用python py6s package实现更加便捷的6S大气校正

首先，简单叙述6s大气校正原理

6s模型里需要输入一系列的与模拟成像日期大气情况的参数，这些参数主要包括：几何参数、气溶胶光学厚度、水蒸汽、臭氧、高程等，同时，这些参数的输入都以查找表这一简单直观的方式获取。最后生成一个线性大气校正公式，逐一得到每波段每像素的6s大气校正后的地表反射率值。

1.几何参数igeom：
如果需要使用6s校正高分辨率卫星影像，一般设置igeom = 0，
同时需要用户自定义输入各必须的几何参数：

solar_z =             # 太阳天顶角 = 90° - 太阳高程角
solar_a =             # 太阳方位角
view_z =              # 卫星天顶角，如果影像为垂直90度观测，则view_z = 0，view_a = 0，但是大部分情况会有一些角度。
view_a =              # 卫星方位角

month =               # 影像成像月份
day =                 # 影像成像具体天

2.大气模式idatm：
（一共有7种内置的大气模式，从0-6，无大气吸收，热带，中纬度夏季，中纬度冬季，亚寒带夏季，亚寒带冬季以及美国标准62大气模式）
我们可以通过影像中心点的经纬度坐标+影像成像月份两个信息共同判断：

# 大气模式类型
if sLatitude > -15 and sLatitude <= 15:
   s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.Tropical)

if sLatitude > 15 and sLatitude <= 45:
   if s.geometry.month > 4 and s.geometry.month <= 9:
      s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeSummer)
   else:
     s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.MidlatitudeWinter)

if sLatitude > 45 and sLatitude <= 60:
   if s.geometry.month > 4 and s.geometry.month <= 9:
      s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticSummer)
  else:
     s.atmos_profile = AtmosProfile.PredefinedType(AtmosProfile.SubarcticWinter)

3.选择气溶胶模式iaer：
（一共有7种内置的气溶胶模式，从0-7，0-无气溶胶模式，1-大陆型，2-海洋型，3-城市型，5-沙漠背景的shettle模型，6-生物燃烧模型，7-平流层模型）

aero_profile = AeroProfile.PredefinedType(AeroProfile.Continental)
# ***这里注意下，是AeroProfile.PredefinedType而不是AtmosProfile.PredefinedType***

4.设置 550nm气溶胶光学厚度,根据日期从MODIS的月度气溶胶产品中获取。
6S模型中可以通过定义两种值的其中之一来确定，可见度v：visibility或者气溶胶光学厚度aot550；然而需要注意py6s中目前只能使用气溶胶光学厚度aot550参数来定义。
该值可以通过查询研究区最邻近站点对应日期的550nm对应的值；
查询网址：

https://aeronet.gsfc.nasa.gov/cgi-bin/type_piece_of_map_opera_v2_new

aot550 = aot

5.输入目标海拔的高度，可使用研究区四个角上的经纬度值确定一个矢量研究区范围，用该范围掩膜DEM值，得到该研究区的平均高程值：

# 通过研究去区的范围去求DEM高度。
pointUL = dict()
pointDR = dict()
pointUL["lat"] = point1lat
pointUL["lon"] = point1lon
pointDR["lat"] = point4lat
pointDR["lon"] = point2lon
meanDEM = (MeanDEM(pointUL, pointDR)) * 0.001    # MeanDEM（）为自定义计算高程的函数

6.传感器Xpp高度的输入，xpp= -1000 表明传感器为星载传感器，xpp=0 表明传感器在地面上-100;由此可以看出，6S模型也可以对无人机等传感器获得的影像进行大气校正。

# 研究区海拔、卫星传感器轨道高度
s.altitudes = Altitudes()
s.altitudes.set_target_custom_altitude(meanDEM)
s.altitudes.set_sensor_satellite_level()        # 选择卫星传感器

7.光谱条件iwave的输入，其中6S内置了很多传感器的光谱条件，如MODIS，TM，AVHRR，SPOT，NOAA等。然而高分辨率商业卫星影像都需要自定义光谱条件。其中光谱响应函数的计算是以2.5nm为间隔, 一般官方提供的光谱响应函数文件的光谱间隔为10nm, 需要自己插值制作为2.5nm间隔。当然还有一个便捷的途径，从ENVI软件里自带的许多光谱响应函数文件.sli中得到，注意ENVI中的光谱间隔为1nm,因此反而需要将间隔增大2.5倍！！！
经过本渣渣艰难的摸索，光谱响应函数是影像你的6S大气校正是否正确的最关键因素，需要特别小心，一种检验方法，可以与ENVI中的Flaash大气校正结果相互印证。
还需要注意的是，有的朋友每个波段最后得到的参数完全一样，那代表你的结果一定有问题，可以先从你的光谱响应矩阵入手检查！！

待校正的影像数据类型、存储方式的要求、单位、缩放尺度

1.存储方式的要求：
从6s大气校正的过程来看，影像的存储格式并没有要求，可以是TIL、XML、DAT、HDF等栅格影像的存储方式；相反重要的是获取存储6s所需参数的头文件信息。

2.输入的数据类型：
首先如果你手里拿到的是一景高分辨率卫星影像的原始数据格式–DN值，你需要进行辐射定标；
一个非常重要的直击灵魂的问题：使用6S模型进行大气校正的输入影像是辐射定标后的表观辐亮度值还是表观反射率值？？
补充答案：需要的是没有经过任何尺度放大或者缩小的表观辐亮度，因此，根据自己的原始文件进行判断。如果手里的文件是DN值，根据官方提供的Gain和Bias值定标；如果手里的文件是TOA数据，则需要根据TOA与TOR的转换公式计算得到！！

2.1计算表观辐亮度值（）的公式：

Radiance=gain*DN+offset    

每个波段对应的gain和offset值，都可以通过影像的头文件中对应名查询到。

2.2计算表观反射率（Top of Atmosphere Reflectance）的公式：

ρ=π*L*d2/（ESUN*cos（θ））
rad = toa/ ρ

其中ρ为表观反射率，L为表观辐亮度，d为日地距离，ESUN为太阳平均辐射强度，θ为太阳高度角。

…还有些没写完，最近忙着做实验，空了补充

❤❤如果这篇文章对你有所帮助，希望大家不吝啬给我一个鼓励的赞~~❤❤