CN104279967B

CN104279967B - 基于高光谱图像的气溶胶光学厚度反演方法

Info

Publication number: CN104279967B
Application number: CN201410222395.XA
Authority: CN
Inventors: 李传荣; 王宁; 刘耀开; 马灵玲; 高彩霞; 钱永刚; 王新鸿; 朱博
Original assignee: Academy of Opto Electronics of CAS
Current assignee: Academy of Opto Electronics of CAS
Priority date: 2014-05-23
Filing date: 2014-05-23
Publication date: 2017-04-19
Anticipated expiration: 2034-05-23
Also published as: CN104279967A

Abstract

本发明提供了一种基于高光谱图像的气溶胶光学厚度反演方法。该气溶胶光学厚度反演方法通过选择高光谱连续波段中对气溶胶光学厚度AOD比较敏感的波段，将气溶胶光学厚度AOD作为未知量引入到反演过程中进行求解不需要从图像中选取暗目标就能够反演气溶胶光学厚度AOD，适用于不同地表类型覆盖区域的气溶胶光学厚度AOD反演。

Description

基于高光谱图像的气溶胶光学厚度反演方法

技术领域

本发明涉及遥感技术领域，尤其涉及一种基于高光谱图像的气溶胶光学厚度反演方法。

背景技术

大气气溶胶是大气与悬浮在其中的固体和液体微粒共同组成的多相体系。气溶胶的光学特性主要表现在能够吸收和散射太阳光，进而影响地气系统辐射平衡，也是影响大气能见度的主要因素。整层大气因气溶胶散射和吸收造成的总光学厚度称为气溶胶光学厚度(Aerosol Optical Depth，AOD)。AOD是描述气溶胶光学特性的重要参数，也是影响遥感影像质量的重要因素，特别在定量遥感领域中，如何精确获取AOD信息关系到大气校正以及地表真实物理特性提取的精度，被遥感领域研究者所广泛关注与重视。此外，气溶胶作为大气颗粒物主要组成部分，获取AOD信息也能够在一定程度上反映区域气溶胶浓度和含量信息，为环境和健康保护领域提供决策支持。

气溶胶粒子的光学特性，例如吸收和散射等，是随波长而变化的。同时，由于气溶胶粒子类型多样、物理化学组成复杂、时空变化性强，导致AOD在时间、空间、光谱上都存在较大的差异。AOD难以通过直接的仪器测量得到，一般需要通过相应的反演过程。目前获取AOD的方法可大致分为地基观测和遥感观测两大类型。地基观测在近些年得到了巨大的发展，特别随着世界各国对气溶胶观测的重视，目前全球范围内已经建立了数百个地基测量站点，构建了气溶胶自动观测网(AERONET)，加拿大太阳光度计网(AEROCAN)，法国气溶胶网络(PHOTONS)等为代表的网络。我国目前也构建了中国气溶胶观测网CAeroNet和太阳分光辐射观测网CSHNET。地基观测多利用太阳光度计等高精度仪器，通过多个通道对大气进行观测，根据气溶胶在不同波长上消光能力的差异反演出AOD，由于反演中不受地面复杂背景影响，AOD反演精度较高。然而采用地基的方式仅能实现点观测，不能获得AOD空间分布信息。

遥感观测是目前比较常用的AOD获取方法，其优势在于可快速获取区域或者全球范围的AOD状况。目前，基于遥感数据反演AOD的代表性方法主要有：浓密植被暗像元法(Dense Dark Vegetation，DDV)和对比度降低法(Contrast Reduction，CR)。

根据大气辐射传输理论，遥感图像中同时反映地物反射以及大气吸收、散射等光学特性。若要精确反演AOD，前提是对地物反射有着比较好的估计。DDV法利用水体、浓密植被等低反射率物体对应的暗像元进行AOD反演。这些低反射率物体的反射率容易估计且精度较高，并且信号中大气信息占据重要地位适合于AOD的反演。DDV法首先需要找到图像中具有一定面积的低反射率区域，利用2.13μm附近通道几乎不受大气影响的特点，从遥感图像获取2.13μm通道地表反射率，进而估算该区域红光、蓝光波段反射率，最后结合暗像元在红光、蓝光波段的入瞳辐亮度、地表反射率以及对大气和气溶胶模式的先验知识，在辐射传输模型或者快速查找表的支持下反演得到AOD。但实际应用中，由于该方法在暗像元寻找、红光蓝光反射率估计等方面假设条件过多且苛刻，方法的应用受到局限，同时，该方法对于高反射区域(如城市、沙漠)的AOD反演会带来很大误差，甚至失效。

CR方法假定同一地区在一段时间内地表反射率不变的情况下，用清洁日大气作为污染日大气的参考进而反演AOD(Kaufman et al.，1990)。CR方法虽然弥补了DDV方法不适用于高反射地表上空的气溶胶光学厚度，但使用该方法时必须有同一地区的多幅遥感图像，而且要求其中一天为清洁日并已知当时的气溶胶和地表反射率状态，应用范围比较窄。此外，基于CR的AOD反演同样需要已知地表反射率作为基础，这在这一定程度上限制了遥感反演AOD的精度和应用范围。

发明内容

(一)要解决的技术问题

鉴于上述技术问题，本发明提供了一种基于高光谱图像反演气溶胶光学厚度的方法，以提高气溶胶光学厚度的测量精度。

(二)技术方案

本发明基于高光谱图像的气溶胶光学厚度反演方法包括：步骤A，在红光波段0.65μm～0.76μm之间选择第一观测波段；在蓝光波段0.4μm～0.5μm之间选择第二观测波段；在近红外或者短波红外波段之间选择第三观测波段；步骤C，在实测的高光谱图像中，抽取三个观测波段的图像数据，并在图像范围内选取P个均匀区域，分别获取三个观测波段的实际表观辐亮度，其中，表示在第p个均匀区域像元第i个观测波段的实际表观辐亮度；i＝1、2、3；p＝1、2、……、P；步骤D，获取反演过程所需的初始参数，该初始参数包括：初始气溶胶光学厚度AOD、地表反射率ρ_ground、大气水汽含量WVC、太阳天顶角SZA、相对方位角RAA、传感器观测高度OBS_H、地表高程ELEV、气溶胶模式、大气模式参数，对于每一均匀区域，均具有相应的地表反射率ρ_ground，步骤E，利用当前气溶胶光学厚度AOD及地表反射率ρ_ground、大气水汽含量WVC、太阳天顶角SZA、相对方位角RAA、传感器观测高度OBS_H、地表高程ELEV、气溶胶模式、大气模式参数，在表观辐亮度查找表中查询或插值得到三个观测波段对应的输出参数，该输出参数包括：太阳到地表的直射大气透过率太阳到地表的散射大气透过率t_d(μ_s)、地表到传感器的直射大气透过率地表到传感器的散射大气透过率t_d(μ_v)、大气程辐射L_path和大气半球反照率S；步骤F，依照如下公式，由输出参数计算P个均匀区域内像元三个观测波段的模拟表观辐亮度，其中表示在第p个均匀区域第i个观测波段的模拟表观辐亮度：

其中，大气层顶太阳辐照度E_s由MODTRAN辐射传输模型模拟得到，ρ_ground为所求区域的地表反射率；步骤G，在所选择的三个观测波段上，利用实际观测的不同区域像元表观辐亮度和模型模拟得到的表观辐亮度建立代价函数：

其中，λ_i第i个波段对应的波长；步骤H，判断上述代价函数是否达到最小值，如果是，执行步骤J；否则执行步骤I；步骤I，调整气溶胶光学厚度AOD，执行步骤E；以及步骤J，选取代价函数最小值对应的气溶胶光学厚度AOD作为最终的气溶胶光学厚度AOD。

(三)有益效果

从上述技术方案可以看出，本发明基于高光谱图像的气溶胶光学厚度反演方法具有以下有益效果：

(1)不需要从图像中选取暗目标就能够反演气溶胶光学厚度AOD，适用于不同地表类型覆盖区域，尤其是城市、沙漠等高反射区域，的气溶胶光学厚度AOD反演；

(2)利用高光谱遥感图像自身的信息反演气溶胶光学厚度AOD，直接为后续大气校正提供气溶胶光学厚度AOD的空间分布信息，克服了以往地表特性高精度反演需要同步实测气溶胶光学厚度AOD所带来的困难和问题；

(3)基于MODTRAN辐射传输模型构建查找表，可以通过从查找表中查找离线模拟的数据，大大提高了算法计算效率；

(4)考虑了不同大气参数、传感器观测几何条件在空间分布上的差异，因素更加全面，反演获取的气溶胶光学厚度AOD精度更高。

附图说明

图1为本发明实施例基于高光谱图像的气溶胶光学厚度反演方法的流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。需要说明的是，在附图或说明书描述中，相似或相同的部分都使用相同的图号。附图中未绘示或描述的实现方式，为所属技术领域中普通技术人员所知的形式。另外，虽然本文可提供包含特定值的参数的示范，但应了解，参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应的值。

本发明通过选择高光谱连续波段中对气溶胶光学厚度AOD和地表反射比较敏感的波段，将气溶胶光学厚度AOD和目标地物的地表反射率一并作为未知量引入到反演过程中进行求解，从而获取气溶胶光学厚度AOD的厚度。

首先，将本发明基于高光谱图像反演气溶胶光学厚度方法的原理进行如下介绍：

假设陆地表面均一并且朗伯(即地表反射各个方向的能量均一)的情况下，依据大气辐射传输原理，传感器观测到的表观辐亮度L_toa可以表示为：

其中，E_s为大气层顶太阳辐照度，ρ_ground为地表反射率，μ_s为太阳天顶角余弦值，和t_d(μ_s)分别为太阳到地表的直射和散射透过率，和t_d(μ_v)为地表到传感器的直射和散射透过率，L_path为大气程辐射，S为大气半球反照率。

由表观辐射亮度L_toa可以计算相应的表观反射率ρ_opp，可表示为：

其中，为整层大气反射率。需要注意，以上各物理量量均是波长的函数，为了表示简便在此省略了波长项。

从公式(2)可以看出，表观反射率ρ_opp既是气溶胶光学厚度AOD的函数，又是地表反射率ρ_ground的函数，具体来讲：

当地表反射率ρ_ground<<1时，卫星观测的表观反射率主要取决于大气贡献项，即公式(2)式右侧第1项，它是气溶胶光学厚度AOD和观测几何参数的函数；

地表反射率ρ_ground越大，地表的贡献，即公式(2)式右侧第2项的贡献也越大。

因此，若地表为一个黑体，即地表反射率ρ_ground为0，传感器观测到的反射率就是大气反射率，只与大气自身的光学特性有关。因此，根据该区域的气溶胶和大气模式，即可根据大气反射率推断出气溶胶光学厚度AOD。但对于实际的遥感图像而言，往往很难选取到地表反射率近似为0的地物，即使对于水体、浓密植被等近似黑体的地物，其反射率也在0.1左右，而且对于城市、沙漠等高反射地区，如何精确获取地表反射率并提高气溶胶光学厚度AOD反演精度是目前的重要技术瓶颈。

基于上述原理介绍，以下对本发明实施例进行详细说明。

在本发明的一个示例性实施例中，提供了一种基于高光谱图像反演气溶胶光学厚度的方法。请参照图1，本实施例基于高光谱图像反演气溶胶光学厚度的方法包括：

步骤A，在红光波段0.65μm～0.76μm之间选择第一观测波段；在蓝光波段0.4μm～0.5μm之间选择第二观测波段；在近红外或者短波红外波段之间选择第三观测波段；

高光谱数据的优势在于其图谱合一的特点，它同时提供数百幅图像和数百个连续波段的光谱信息，这对同时从空间和光谱维研究地物目标起到了很大的作用。本发明为了很好地应用这一优势，需要从高光谱图像中挑选出对气溶胶光学厚度AOD和地表反射敏感的连续波段。

需要注意的是，在红光波段和蓝光波段进行选择时，需要避开大气吸收波段(如O₂、水汽、CO₂和O₃的吸收带附近的波段)。

步骤B，在上述三个观测波段内，利用MODTRAN大气辐射传输模型，由不同的初始参数计算得到输出参数，利用初始参数和输出参数构建出三个观测波段内对应不同初始参数条件下的表观辐亮度查找表(LUT)；其中：

初始参数包括：气溶胶光学厚度AOD、地表反射率ρ_ground、大气水汽含量WVC、太阳天顶角SZA、相对方位角RAA、传感器观测高度OBS_H、地表高程ELEV、气溶胶模式、大气模式参数；

输出参数包括：太阳到地表的直射大气透过率太阳到地表的散射大气透过率t_d(μ_s)、地表到传感器的直射大气透过率地表到传感器的散射大气透过率t_d(μ_v)、大气程辐射L_path和大气半球反照率S；

其中，MODTRAN(MODerate resolution atmospheric TRANsmission)为一种常用的大气辐射传输模型，能够模拟0.2到100μm范围内电磁辐射在大气中的传输情况，关于该MODTRAN大气辐射传输模型的具体内容，可参见MODTRAN用户手册(Berk等，USER’S MANUAL)。

上述的气溶胶模式是指城市型、乡村型、海洋型气溶胶模型。

上述的大气模式是指中纬度夏季大气模式、中纬度冬季大气模式以及标准大气模式。

上述的大气透过率是指电磁波通过大气中某个给定路径长度后的辐射能与入射辐射能之比。

所构建的查找表建立了不同初始参数及输出结果之间的关系，可以看作是离线状态下大气辐射传输模型的近似。构建查找表时不可能穷尽所有的参数组合情况，因此，为了确定各初始参数的模拟间隔，需要对各初始参数进行敏感性分析，以保证在采样间隔内的模拟的输出参数接近线性变化。表1所示为经过敏感性分析所确定的各初始参数的模拟采样间隔。

表1MODTRAN辐射传输模型模拟的初始参数采样间隔

步骤C，在实测的高光谱图像中，抽取步骤A所选择的三个观测波段图像数据，并在图像范围内选取不少于2个均匀区域，分别获取实际观测的三个观测波段表观辐亮度其中下标p表示不同区域内的像元，其中，p＝1、2、……、P，P为均匀区域的个数；

本发明中，P≥2，优选地，P为4或者5。

步骤D，获取反演过程所需的初始参数；

该初始参数：包括初始气溶胶光学厚度AOD₀、地表反射率ρ_ground、图像获取时刻大气水汽含量WVC、太阳天顶角SZA、相对方位角RAA、传感器观测高度OBS_H、地表高程ELEV、气溶胶模式、大气模式参数；

本方法同时对气溶胶光学厚度AOD和地表反射率进行迭代反演。因此，需要给定气溶胶光学厚度AOD和地表反射率的初始值，其中，气溶胶光学厚度AOD初始值为标准大气情况下的估计值(如0.1)；不同区域像元的地表反射率ρ_ground的初始值利用公式2计算得到的每个像元的表观反射率ρ_opp，其中，太阳天顶角余弦值μ_s由图像获取时刻的太阳天顶角计算得到，大气层顶太阳辐照度E_s利用MODTRAN自带的数据库得到(Berk等，USER’S MANUAL)。

反演过程中需要的其它参数，包括CWV、SZA、RAA、OBS_H、ELEV，根据高光谱图像实际获取的时刻、高度等辅助信息中获取。

这些输入的用途是为了从LUT中获得对应的大气透过率、程辐射和大气半球反照率信息，并由此模拟图像像元的表观辐亮度。

步骤E，利用当前气溶胶光学厚度AOD及地表反射率ρ_ground、大气水汽含量WVC、太阳天顶角SZA、相对方位角RAA、传感器观测高度OBS_H、地表高程ELEV、气溶胶模式、大气模式参数，从LUT中查询插值得到三个观测波段对应的6个大气参数，这6个参数是太阳到地表的直射大气透过率太阳到地表的散射大气透过率t_d(μ_s)、地表到传感器的直射大气透过率地表到传感器的散射大气透过率t_d(μ_v)、程辐射L_path和大气半球反照率(S)；

步骤F，依照如下公式，由当前气溶胶光学厚度AOD、大气层顶太阳辐照度(E_s)、不同区域像元的地表反射率ρ_ground、太阳天顶角余弦值μ_s、太阳到地表的直射大气透过率太阳到地表的散射大气透过率t_d(μ_s)、地表到传感器的直射大气透过率地表到传感器的散射大气透过率t_d(μ_v)、程辐射L_path和大气半球反照率S计算不同区域像元三个观测波段的模拟表观辐亮度

其中，大气层顶太阳辐照度E_s可直接由MODTRAN辐射传输模型模拟得到，其在本领域内已经为技术人员所熟知，此处不再详细解释。每个区域像元模拟计算的方式相同，在公式中省略对应的下标。

步骤G，在所选择的三个观测波段上，利用实际观测的不同区域像元表观辐亮度和模型模拟得到的表观辐亮度建立代价函数，p＝1、2、……、P；

对于高光谱图像中从高光谱图像中提取出对应的表观辐射度在同样的太阳天顶角、方位角、观测天顶角、方位角等观测条件下，构建模拟和实际观测的表观辐亮度二者之间的代价函数如下：

其中，i为高光谱数据对应的波段，λ_i第i个波段对应的波长。

步骤H，判断上述代价函数是否小于预设值，如果是，执行步骤J；否则执行步骤I；

本实施例中，对于公式4所示的代价函数而言，该预设值取1×10^-6，当取其他形式的代价函数时，该预设值可以根据需要适当调整。

步骤I，调整气溶胶光学厚度AOD，执行步骤D；

通过不断调整气溶胶光学厚度值重新插值得到模拟的表观服辐亮度，并利用POWELL优化算法对代价函数进行优化求解。利用POWELL优化算法求的修正值的过程对于本领域技术人员是已知的，此处不再赘述。

步骤J，当代价函数达到最小值时，结束优化过程，选取代价函数最小值对应的气溶胶光学厚度AOD作为最终的气溶胶光学厚度AOD。

至此，已经结合附图对本实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明基于高光谱图像的气溶胶光学厚度反演方法有了清楚的认识。

此外，上述对各元件的定义并不仅限于实施方式中提到的各种具体结构或形状，本领域的普通技术人员可对其进行简单地熟知地替换。

综上所述，本发明提供一种基于高光谱图像的气溶胶光学厚度反演方法。该气溶胶光学厚度反演方法通过选择高光谱连续波段中对气溶胶光学厚度AOD比较敏感的波段，将气溶胶光学厚度AOD一并作为未知量引入到反演过程中进行求解不需要从图像中选取暗目标就能够反演气溶胶光学厚度AOD，适用于不同地表类型覆盖区域的气溶胶光学厚度AOD反演。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims

1.一种基于高光谱图像的气溶胶光学厚度反演方法，其特征在于，包括：

步骤B，在三个观测波段内，利用MODTRAN大气辐射传输模型，由不同的初始参数计算得到输出参数，利用初始参数和输出参数构建出三个观测波段内对应不同初始参数条件下的表观辐亮度查找表；

步骤C，在实测的高光谱图像中，抽取三个观测波段的图像数据，并在图像范围内选取P个均匀区域，分别获取三个观测波段的实际表观辐亮度，其中，表示在第p个均匀区域像元第i个观测波段的实际表观辐亮度；i＝1、2、3；p＝1、2、……、P；

步骤D，获取反演过程所需的初始参数；

该初始参数包括：初始气溶胶光学厚度AOD、地表反射率ρ_ground、大气水汽含量WVC、太阳天顶角SZA、相对方位角RAA、传感器观测高度OBS_H、地表高程ELEV、气溶胶模式、大气模式参数，对于每一均匀区域，均具有相应的地表反射率ρ_ground；

步骤E，利用当前气溶胶光学厚度AOD及地表反射率ρ_ground、大气水汽含量WVC、太阳天顶角SZA、相对方位角RAA、传感器观测高度OBS_H、地表高程ELEV、气溶胶模式、大气模式参数，在表观辐亮度查找表LUT中查询或插值得到三个观测波段对应的输出参数；

该输出参数包括：太阳到地表的直射大气透过率太阳到地表的散射大气透过率t_d(μ_s)、地表到传感器的直射大气透过率地表到传感器的散射大气透过率t_d(μ_v)、大气程辐射L_path和大气半球反照率S；

步骤F，依照如下公式，由输出参数计算P个均匀区域内像元三个观测波段的模拟表观辐亮度，其中表示在第p个均匀区域第i个观测波段的模拟表观辐亮度：

其中，大气层顶太阳辐照度E_s由MODTRAN辐射传输模型模拟得到，ρ_ground为所求区域的地表反射率；

步骤G，在所选择的三个观测波段上，利用实际观测的不同区域像元表观辐亮度和模型模拟得到的表观辐亮度建立代价函数；

步骤I，调整气溶胶光学厚度AOD，重新执行步骤E；以及

步骤J，选取代价函数最小值对应的气溶胶光学厚度AOD作为最终的气溶胶光学厚度AOD。

2.根据权利要求1所述的气溶胶光学厚度反演方法，其特征在于，所述步骤I包括：

利用POWELL优化算法，对代价函数进行优化求解，计算出气溶胶光学厚度AOD的修正值，并将此修正值加入前一气溶胶光学厚度上，得到调整后的气溶胶光学厚度AOD。

3.根据权利要求1所述的气溶胶光学厚度反演方法，其特征在于，所述P≥2。

4.根据权利要求3所述的气溶胶光学厚度反演方法，其特征在于，所述P＝4或5。

5.根据权利要求1所述的气溶胶光学厚度反演方法，其特征在于，所述步骤G中，代价函数如下：

其中，λ_i第i个波段对应的波长。

6.根据权利要求5所述的气溶胶光学厚度反演方法，其特征在于，所述步骤H中预设值为1×10^-6。

7.根据权利要求1至6中任一项所述的气溶胶光学厚度反演方法，其特征在于，所述步骤A中，在选择第一观测波段时，避开O₂、水汽、CO₂和O₃的吸收带。

8.根据权利要求1至6中任一项所述的气溶胶光学厚度反演方法，其特征在于，所述初始气溶胶光学厚度AOD为标准大气情况下的估计值。

9.根据权利要求8所述的气溶胶光学厚度反演方法，其特征在于，所述初始气溶胶光学厚度AOD为0.1。