一种气溶胶光学厚度反演方法及装置

摘要

本申请提供了一种气溶胶光学厚度反演方法及装置，涉及卫星遥感影像技术领域，包括：获取FY4A AGRI传感器的遥感影像，从遥感影像中获取每个像元的观测几何参数和蓝光数据；基于每个像元的蓝光数据得到第一蓝光波段表观反射率；利用所述像元的观测几何参数和预先建立的查找表，通过理论模拟方法计算每个像元的多个模拟气溶胶光学厚度对应的第二蓝光波段表观反射率；利用所述像元的第一蓝光波段表观反射率和多个模拟气溶胶光学厚度对应的第二蓝光波段表观反射率，拟合得到每个像元的一次线性方程；基于每个像元的第一蓝光波段表观反射率和对应的一次线性方程，计算每个像元的气溶胶光学厚度。本申请提高了气溶胶光学厚度反演精度。

说明书

技术领域

本申请涉及卫星遥感影像技术领域，尤其是涉及一种气溶胶光学厚度反演方法及装置。

背景技术

近年来，大气污染问题已成为影响人们身体健康、日常生活及城市可持续发展的重要因素，精确遥感反演气溶胶光学厚度，全面监测大气污染变化对环境治理和城市规划具有重要需求价值和科学意义。对FY4A AGRI数据进行几何校正等预处理步骤，并利用MODIS蓝光波段数据对FY4A AGRI蓝光波段进行交叉辐射定标，获得较为准确的DN值后，并利用深蓝算法（Deep Blue,DB）进行气溶胶反演。

目前的气溶胶产品数据主要是通过MODIS数据或者利用如Landsat卫星数据进行反演得到的。FY4A作为我国第二代静止轨道气象卫星，是中国新一代静止气象卫星的首发星。其搭载的干涉式大气垂直探测仪与静止轨道扫描成像辐射计“联手”，在全球首次实现静止轨道上三维大气的立体监测。此外，“风云四号”能够提供针对陆面、水体、闪电、空间天气的持续监测数据，能清晰区分云的不同形态和高、中层水汽，也具备捕捉气溶胶和雪的能力，但目前未有利用AGRI传感器反演气溶胶光学厚度的技术方案。

发明内容

有鉴于此，本申请提供了一种气溶胶光学厚度反演方法及装置，以解决上述技术问题。

第一方面，本申请实施例提供了一种气溶胶光学厚度反演方法，包括：

获取FY4A AGRI传感器的遥感影像，从遥感影像中获取每个像元的观测几何参数和蓝光数据；

利用预先建立的蓝光数据和蓝光波段表观反射率的线性方程，基于每个像元的蓝光数据得到每个像元的第一蓝光波段表观反射率；

利用所述像元的观测几何参数和预先建立的查找表，通过理论模拟方法计算每个像元的多个模拟气溶胶光学厚度对应的第二蓝光波段表观反射率；

利用所述像元的第一蓝光波段表观反射率和多个模拟气溶胶光学厚度对应的第二蓝光波段表观反射率，拟合得到每个像元的一次线性方程，所述一次线性方程以蓝光波段表观反射率为自变量，气溶胶光学厚度为因变量；

基于每个像元的第一蓝光波段表观反射率和对应的一次线性方程，计算每个像元的气溶胶光学厚度。

进一步，蓝光数据和蓝光波段表观反射率的线性方程的拟合步骤包括：

获取时空匹配的FY4A AGRI传感器数据和MODIS产品数据；

获取FY4A AGRI传感器数据中的卫星天顶角对应的MODIS产品数据的卫星天顶角；

获取FY4A AGRI传感器数据中的相对方位角对应的MODIS产品数据的相对方位角；

获取MODIS产品数据的卫星天顶角和相对方位角对应的蓝光波段表观反射率；

的蓝光波段表观反射率，拟合出蓝光数据和蓝光波段表观反射率的线性方程的系数

其中，

进一步，所述查找表的数据项包括：半球反射率、大气透过率、大气辐射沿大气传输路径的等效反射率、太阳天顶角、卫星天顶角、相对方位角和气溶胶光学厚度；其中，太阳天顶角的9个值包括：0，10，20，30，40，50，60，70或80；卫星天顶角的9个值包括：0，10，20，30，40，50，60，70或80；相对方位角的19个值包括：0，10，20，30，40，50，60，70，80，90，100，110，120，130，140，150，160，170和180。

进一步，观测几何参数包括：太阳天顶角、卫星天顶角和相对方位角；利用所述像元的观测几何参数和预先建立的查找表，通过理论模拟方法计算每个像元的多个模拟气溶胶光学厚度对应的第二蓝光波段表观反射率；包括:

利用所述像元的观测几何参数和预先建立的地表反射率线性方程，得到所述像元的地表反射率；其中，地表反射率线性方程的自变量为MODIS传感器的地表反射率，因变量为FY4A AGRI传感器的地表反射率；

利用所述像元的观测几何参数，在预先建立的查找表中，查询得到多个模拟气溶胶光学厚度对应的半球反射率、大气透过率和大气辐射沿大气传输路径的等效反射率；

利用所述像元的地表反射率、半球反射率、大气透过率和大气辐射沿大气传输路径的等效反射率，计算得到多个模拟气溶胶光学厚度对应第二蓝光波段表观反射率。

进一步，地表反射率线性方程的拟合步骤包括：

分别计算FY4A AGRI传感器和MODIS传感器在蓝光波段观测的地表反射率

其中，

利用多个

其中，

进一步，利用所述像元的观测几何参数和预先建立的地表反射率线性方程，得到所述像元的地表反射率；包括：

获取所述像元的太阳天顶角对应的MODIS产品数据的太阳天顶角；

获取所述像元的卫星天顶角对应的MODIS产品数据的卫星天顶角；

获取所述像元的相对方位角对应的MODIS产品数据的相对方位角；

获取MODIS产品数据的太阳天顶角、卫星天顶角和相对方位角对应的MODIS产品数据的地表反射率；

将MODIS产品数据的地表反射率作为自变量代入地表反射率线性方程，得到所述像元的地表反射率。

进一步，利用所述像元的观测几何参数，在预先建立的查找表中，查询得到多个气溶胶光学厚度对应的半球反射率、大气透过率和大气辐射沿大气传输路径的等效反射率，包括：

对于查找表中的每个气溶胶光学厚度，执行下述步骤：

获取所述像元的太阳天顶角相邻的两个太阳天顶角；获取所述像元的卫星天顶角在查找表中相邻的两个卫星天顶角；获取所述像元的相对方位角在查找表中相邻的两个相对方位角；

将两个太阳天顶角、两个卫星天顶角和两个相对方位角进行组合，得到8组太阳天顶角、卫星天顶角和相对方位角的组合；

通过查找表得到8个组合对应的8个半球反射率，8个大气透过率和8个大气辐射沿大气传输路径的等效反射率；

计算8个半球反射率的平均值，作为所述像元的半球反射率；

计算8个大气透过率的平均值，作为所述像元的大气透过率；

计算8个大气辐射沿大气传输路径的等效反射率，作为所述像元的大气辐射沿大气传输路径的等效反射率；

由此得到多个气溶胶光学厚度对应的半球反射率、大气透过率和大气辐射沿大气传输路径的等效反射率。

进一步，利用所述像元的第一蓝光波段表观反射率和多个模拟气溶胶光学厚度对应的第二蓝光波段表观反射率，拟合得到每个像元的一次线性方程；包括：

计算所述像元的第一蓝光波段表观反射率和每个第二蓝光波段表观反射率的差的绝对值，获取计算结果中最小和次小对应的两个第二蓝光波段表观反射率和对应的气溶胶光学厚度；

利用每个像元的两个第二蓝光波段表观反射率和对应的气溶胶光学厚度，拟合出每个像元的一次线性方程的系数。

第二方面，本申请实施例提供了一种气溶胶光学厚度反演装置，包括：

获取单元，用于获取FY4A AGRI传感器的遥感影像，从遥感影像中获取每个像元的观测几何参数和蓝光数据；

第一计算单元，用于利用预先建立的蓝光数据和蓝光波段表观反射率的线性方程，基于每个像元的蓝光数据得到每个像元的第一蓝光波段表观反射率；

第二计算单元，用于利用所述像元的观测几何参数和预先建立的查找表，通过理论模拟方法计算每个像元的多个模拟气溶胶光学厚度对应的第二蓝光波段表观反射率；

拟合单元，用于利用所述像元的第一蓝光波段表观反射率和多个模拟气溶胶光学厚度对应的第二蓝光波段表观反射率，拟合得到每个像元的一次线性方程，所述一次线性方程以蓝光波段表观反射率为自变量，气溶胶光学厚度为因变量；

第三计算单元，用于基于每个像元的第一蓝光波段表观反射率和对应的一次线性方程，计算每个像元的气溶胶光学厚度。

第三方面，本申请实施例提供了一种电子设备，包括：存储器、处理器和存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如本申请实施例的气溶胶光学厚度反演方法。

第四方面，本申请实施例提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令被处理器执行时实现如本申请实施例的气溶胶光学厚度反演方法。

本申请提高了气溶胶光学厚度反演精度。

附图说明

为了更清楚地说明本申请具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本申请的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本申请实施例提供的气溶胶光学厚度反演方法的流程图；

图2为本申请实施例提供的气溶胶光学厚度反演装置的功能结构图；

图3为本申请实施例提供的电子设备的结构图。

具体实施方式

为使本申请实施例的目的、技术方案和优点更加清楚，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本申请一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本申请实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本申请的实施例的详细描述并非旨在限制要求保护的本申请的范围，而是仅仅表示本申请的选定实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

首先对本申请实施例的设计思想进行简单介绍。

随着遥感数据定量化的广泛应用，对于卫星传感器辐射定标精度的要求也越来越高。传感器本身的辐射响应会随着使用时间增加而衰减，受到太阳辐射等因素的影响，导致传感器接收到的数据是有误差的，因此需要准确检测遥感探测器在运行期间内的辐射响应变化，及时更新辐射定标参数，以保证遥感数据定量化研究的精度。采用MODIS数据对FY4AAGRI传感器蓝光波段进行交叉辐射定标，在不同传感器太阳天顶角、卫星天顶角与相对方位角、通道光谱响应差异校正基础上，通过理论模拟方法计算FY4A AGRI蓝光波段表观反射率，并将其与传感器输出DN值进行位置匹配，进而计算FY4A AGRI传感器蓝光波段辐射定标系数。对于传感器高精度定量应用具有十分重要的意义。

本申请使用FY4A AGRI传感器反演AOD，利用AGRI与MODIS蓝光波段进行交叉辐射定标，得到了更为精确的蓝光波段的定标系数，为AOD的反演提供了更精确的数据基础；使用适用于蓝光波段的反演AOD的深蓝算法，其反演的精度相较于其他原理方法更加精确。

在介绍了本申请实施例的应用场景和设计思想之后，下面对本申请实施例提供的技术方案进行说明。

如图1所示，本申请实施例提供了一种气溶胶光学厚度反演方法，包括如下步骤：

步骤101：获取FY4A AGRI传感器的遥感影像，从遥感影像中获取每个像元的观测几何参数和蓝光数据；

步骤102：利用预先建立的蓝光数据和蓝光波段表观反射率的线性方程，基于每个像元的蓝光数据得到每个像元的第一蓝光波段表观反射率；

为了得到蓝光数据和蓝光波段表观反射率的线性方程，需要利用历史数据对数据进行拟合，但是由于FY4A AGRI传感器的蓝光波段表观反射率无法直接获取到，而从MODIS产品数据中能够直接获取到蓝光波段表观反射率，若将MODIS产品数据的蓝光波段表观反射率用于FY4A AGRI传感器的蓝光波段表观反射率，需要将MODIS产品数据和FY4A AGRI数据进行准确配准。同时，需要充分考虑成像时刻卫星天顶角、相对方位角等角度的影响，以及两个通道之间光谱响应的差异。

由于卫星天顶角变化，探测路径会因靠近扫描线边缘而增加，使得大气衰减严重，图像上显示越暗；同时，同一位置的卫星天顶角在不同轨道、不同时相及不同卫星影像上均有明显变化。基于以上影响，需要进行不同卫星天顶角数据校正。校正模型如下式（1）：

式中，

相对方位角的校正主要采用回归分析方法进行处理。将MODIS相对方位角对应的数据转化为0度相对方位角对应的数据。然后再将其转化为FY4A AGRI卫星相对方位角所对应的表观反射率。

将MODIS表观反射率经角度差异校正、光谱响应差异校正后，得到理论上的FY4A卫星蓝光波段表观反射率，将其与同一位置的FY4A AGRI蓝光波段原始DN值进行线性相关分析，可得出两者之间线性关系，即求得FY4A AGRI的蓝光波段表观反射率。

具体的，蓝光数据和蓝光波段表观反射率的线性方程的拟合步骤包括：

获取时空匹配的FY4A AGRI传感器数据和MODIS产品数据；

获取FY4A AGRI传感器数据中的卫星天顶角对应的MODIS产品数据的卫星天顶角；

获取FY4A AGRI传感器数据中的相对方位角对应的MODIS产品数据的相对方位角；

获取MODIS产品数据的卫星天顶角和相对方位角对应的蓝光波段表观反射率；

利用FY4A AGRI传感器数据中的蓝光数据和对应的MODIS产品数据的蓝光波段表观反射率，拟合出蓝光数据和蓝光波段表观反射率的线性方程的系数

其中，

步骤103：利用所述像元的观测几何参数和预先建立的查找表，通过理论模拟方法计算每个像元的多个模拟气溶胶光学厚度对应的第二蓝光波段表观反射率；

本实施例的大气辐射传输模型如下：

目前多数遥感反演手段均是基于大气辐射传输方程来反演气溶胶光学厚度值。卫星传感器所接收到的大气上界表观反射率

式中，

在晴朗无云的天气条件下，不考虑气体吸收情况，假设大气处于一个水平、均一的状态，对于均匀朗伯面的下垫面，卫星接收到的总辐射亮度

式中，

根据上式（4），利用垂直方向上入射的太阳辐射率

式中，

其中，

本实施例的查找表的数据项包括：半球反射率、大气透过率、大气辐射沿大气传输路径的等效反射率、太阳天顶角、卫星天顶角、相对方位角和气溶胶光学厚度；其中，太阳天顶角的9个值包括：0，10，20，30，40，50，60，70或80；卫星天顶角的9个值包括：0，10，20，30，40，50，60，70或80；相对方位角的19个值包括：0，10，20，30，40，50，60，70，80，90，100，110，120，130，140，150，160，170和180。如表1所示：

表1

根据大气辐射传输方程（5）可知，若获取对应光学厚度下的路径辐射反射率

本实施例中，观测几何参数包括：太阳天顶角、卫星天顶角和相对方位角；该步骤具体包括：

利用所述像元的观测几何参数和预先建立的地表反射率线性方程，得到所述像元的地表反射率；其中，地表反射率线性方程的自变量为MODIS传感器的地表反射率，因变量为FY4A AGRI传感器的地表反射率；

利用所述像元的观测几何参数，在预先建立的查找表中，查询得到多个模拟气溶胶光学厚度对应的半球反射率、大气透过率和大气辐射沿大气传输路径的等效反射率；

利用所述像元的地表反射率、半球反射率、大气透过率和大气辐射沿大气传输路径的等效反射率，计算得到多个模拟气溶胶光学厚度对应第二蓝光波段表观反射率。

其中，地表反射率线性方程的拟合步骤包括：

用FY4A AGRI传感器数据的蓝光波段进行AOD反演，采用MODIS地表反射率产品MOD09A1数据构建地表反射率数据，拟合地表反射率线性公式，其用8天合成地表反射率产品目的在于消除云覆盖的影响和地表反射率瞬时变化的影响，并且遥感影像周期较短，短时间内同一位置的地表反射率几乎可视为不变。

由于FY4A AGRI传感器与MODIS在蓝光波段的波段响应函数差异，需要采用以下地表反射率线性公式进行通道光谱响应差异校正：

其中，

以上两个公式中，

利用公式（8）和（9）拟合出不同地物对应的FY4A AGRI、MODIS蓝光波段反射率后，利用公式（7）即可回归得到通道光谱响应差异的校正系数

本实施例中，利用所述像元的观测几何参数和预先建立的地表反射率线性方程，得到所述像元的地表反射率；包括：

获取所述像元的太阳天顶角对应的MODIS产品数据的太阳天顶角；

获取所述像元的卫星天顶角对应的MODIS产品数据的卫星天顶角；

获取所述像元的相对方位角对应的MODIS产品数据的相对方位角；

获取MODIS产品数据的太阳天顶角、卫星天顶角和相对方位角对应的MODIS产品数据的地表反射率；

将MODIS产品数据的地表反射率作为自变量代入地表反射率线性方程，得到所述像元的地表反射率。

本实施例中，利用所述像元的观测几何参数，在预先建立的查找表中，查询得到多个模拟气溶胶光学厚度对应的半球反射率、大气透过率和大气辐射沿大气传输路径的等效反射率；包括：

对于查找表中的每个气溶胶光学厚度，执行下述步骤：

获取所述像元的太阳天顶角相邻的两个太阳天顶角；获取所述像元的卫星天顶角在查找表中相邻的两个卫星天顶角；获取所述像元的相对方位角在查找表中相邻的两个相对方位角；

将两个太阳天顶角、两个卫星天顶角和两个相对方位角进行组合，得到8组太阳天顶角、卫星天顶角和相对方位角的组合；

通过查找表得到8个组合对应的8个半球反射率，8个大气透过率和8个大气辐射沿大气传输路径的等效反射率；

计算8个半球反射率的平均值，作为所述像元的半球反射率；

计算8个大气透过率的平均值，作为所述像元的大气透过率；

计算8个大气辐射沿大气传输路径的等效反射率，作为所述像元的大气辐射沿大气传输路径的等效反射率；

由此得到多个气溶胶光学厚度对应的半球反射率、大气透过率和大气辐射沿大气传输路径的等效反射率。

步骤104：利用所述像元的第一蓝光波段表观反射率和多个模拟气溶胶光学厚度对应的第二蓝光波段表观反射率，拟合得到每个像元的一次线性方程，所述一次线性方程以蓝光波段表观反射率为自变量，气溶胶光学厚度为因变量；

具体的，计算所述像元的第一蓝光波段表观反射率和每个第二蓝光波段表观反射率的差的绝对值，获取计算结果中最小和次小对应的两个第二蓝光波段表观反射率和对应的气溶胶光学厚度；利用每个像元的两个第二蓝光波段表观反射率和对应的气溶胶光学厚度，拟合出每个像元的一次线性方程的系数。

步骤105：基于每个像元的第一蓝光波段表观反射率和对应的一次线性方程，计算每个像元的气溶胶光学厚度。

将每个像元的第一蓝光波段表观反射率代入步骤104的一次线性方程，得到每个像元的气溶胶光学厚度。

基于上述实施例，本申请实施例提供了一种气溶胶光学厚度反演装置，参阅图2所示，本申请实施例提供的气溶胶光学厚度反演装置200至少包括：

获取单元201，用于获取FY4A AGRI传感器的遥感影像，从遥感影像中获取每个像元的观测几何参数和蓝光数据；

第一计算单元202，用于利用预先建立的蓝光数据和蓝光波段表观反射率的线性方程，基于每个像元的蓝光数据得到每个像元的第一蓝光波段表观反射率；

第二计算单元203，用于利用所述像元的观测几何参数和预先建立的查找表，通过理论模拟方法计算每个像元的多个模拟气溶胶光学厚度对应的第二蓝光波段表观反射率；

拟合单元204，用于利用所述像元的第一蓝光波段表观反射率和多个模拟气溶胶光学厚度对应的第二蓝光波段表观反射率，拟合得到每个像元的一次线性方程，所述一次线性方程以蓝光波段表观反射率为自变量，气溶胶光学厚度为因变量；

第三计算单元205，用于基于每个像元的第一蓝光波段表观反射率和对应的一次线性方程，计算每个像元的气溶胶光学厚度。

需要说明的是，本申请实施例提供的气溶胶光学厚度反演装置200解决技术问题的原理与本申请实施例提供的气溶胶光学厚度反演方法相似，因此，本申请实施例提供的气溶胶光学厚度反演装置200的实施可以参见本申请实施例提供的气溶胶光学厚度反演方法的实施，重复之处不再赘述。

如图3所示，本申请实施例提供的电子设备300至少包括：处理器301、存储器302和存储在存储器302上并可在处理器301上运行的计算机程序，处理器301执行计算机程序时实现本申请实施例提供的气溶胶光学厚度反演方法。

本申请实施例提供的电子设备300还可以包括连接不同组件（包括处理器301和存储器302）的总线303。其中，总线303表示几类总线结构中的一种或多种，包括存储器总线、外围总线、局域总线等。

存储器302可以包括易失性存储器形式的可读介质，例如随机存储器（RandomAccess Memory，RAM）3021和/或高速缓存存储器3022，还可以进一步包括只读存储器（ReadOnly Memory，ROM）3023。

存储器302还可以包括具有一组（至少一个）程序模块3025的程序工具3024，程序模块3025包括但不限于：操作子系统、一个或者多个应用程序、其它程序模块以及程序数据，这些示例中的每一个或某种组合中可能包括网络环境的实现。

电子设备300也可以与一个或多个外部设备304（例如键盘、遥控器等）通信，还可以与一个或者多个使得用户能与电子设备300交互的设备通信（例如手机、电脑等），和/或，与使得电子设备300与一个或多个其它电子设备300进行通信的任何设备（例如路由器、调制解调器等）通信。这种通信可以通过输入/输出（Input /Output，I/O）接口305进行。并且，电子设备300还可以通过网络适配器306与一个或者多个网络（例如局域网（Local AreaNetwork，LAN），广域网（Wide Area Network，WAN）和/或公共网络，例如因特网）通信。如图3所示，网络适配器306通过总线303与电子设备300的其它模块通信。应当理解，尽管图3中未示出，可以结合电子设备300使用其它硬件和/或软件模块，包括但不限于：微代码、设备驱动器、冗余处理器、外部磁盘驱动阵列、磁盘阵列（Redundant Arrays of IndependentDisks，RAID）子系统、磁带驱动器以及数据备份存储子系统等。

需要说明的是，图3所示的电子设备300仅仅是一个示例，不应对本申请实施例的功能和使用范围带来任何限制。

本申请实施例还提供了一种计算机可读存储介质，该计算机可读存储介质存储有计算机指令，该计算机指令被处理器执行时实现本申请实施例提供的气溶胶光学厚度反演方法。

此外，尽管在附图中以特定顺序描述了本申请方法的操作，但是，这并非要求或者暗示必须按照该特定顺序来执行这些操作，或是必须执行全部所示的操作才能实现期望的结果。附加地或备选地，可以省略某些步骤，将多个步骤合并为一个步骤执行，和/或将一个步骤分解为多个步骤执行。

尽管已描述了本申请的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本申请范围的所有变更和修改。

最后应说明的是：以上各实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述各实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的范围。