CN109974665A - 一种针对缺少短波红外数据的气溶胶遥感反演方法及系统 - Google Patents

Abstract

本发明公开了一种针对缺少短波红外数据的气溶胶遥感反演方法及系统。在该方法中考虑到角度和地表类型对波段关系的影响，通过使用多角度数据集构建分角度格网分地物类型的波段回归系数查找表，反演得到不同地物类型表面的气溶胶光学厚度。针对气溶胶光学厚度反演结果对波段之间关系有较强的敏感性的问题，引入气溶胶校正指数(ACI)和NDVI季节性变化特点作为限制条件。本发明的核心技术为基于地物光谱波段直接的相关关系，结合不同太阳观测角度对波段关系的影响，建立气溶胶光学厚度查找表和波段相关关系查找表，通过筛选条件可以快速、准确的反演得到不同地物的气溶胶光学厚度反演结果。通过上述技术手段，提高运行效率和反演精度。

Description

一种针对缺少短波红外数据的气溶胶遥感反演方法及系统

技术领域

本发明属于遥感影像处理技术领域，具体涉及一种针对无短波红外波段数据的气溶胶遥感反演方法及系统，快速反演得到高分辨率的气溶胶光学厚度产品。

背景技术

气溶胶光学厚度是研究大气污染、气溶胶辐射、大气校正等问题的重要地球物理参数，通过卫星影像可以大范围地监测气溶胶变化，在空间尺度上解释气溶胶对地球和大气的影响，在卫星遥感影像中包含地表和大气两部分的信息，两者相互影响，大气中气溶胶的作用导致图像细节损失以及为地表真实信息的定量反演应用带来较大误差，同样因为地表信息的影响也为准确的反演气溶胶带来困难。获取准确的气溶胶光学厚度对地表真实信息的反演以及空气质量研究具有重要意义。

目前使用遥感影像反演气溶胶的算法主要有：暗目标法和深蓝算法，其中暗目标法应用最为广泛，利用在浓密植被覆盖地区可见光波段地表反射率与2.1微米短波红外波段的线性关系反演气溶胶光学厚度，并在MODIS影像中实现业务化运行，对于地表反射率较高的城市、沙漠、裸土覆盖等地区多使用深蓝算法反演气溶胶光学厚度。

然而，受制于数据源空间分辨率的影响，目前的一些气溶胶产品存在一个低空间分辨率的缺点(如MODIS，最高空间分辨率的气溶胶产品为3km)，而在地表结构复杂的城市等地区则需要更高空间分辨率的气溶胶光学厚度信息。然而，对于缺少2.1微米短波红外波段的卫星数据，常规的暗目标法无法使用，比如：环境一号卫星CCD传感器具有30m的高空间分辨率的优点，可以作为反演高空间分辨率气溶胶光学厚度的数据源，但是环境一号CCD数据只包含三个可见光波段和一个近红外波段。

与暗目标法相比深蓝算法需要其他地表反射率数据辅助，因而影响其通用性。暗目标方法中忽略了角度对波段之间相关关系的影响，环境一号CCD数据单台传感器幅宽为360km，观测天顶角最大能达到35度左右，所以太阳观测几何对气溶胶反演的影响不容忽视，另外不同地物之间的线性关系也有较大差异。

因此，针对大幅宽、多角度、高分辨率卫星影像存在的问题有：

①缺少短波红外波段。

②太阳、观测角度对波段相关关系的影响。

③地物类型对波段相关关系的影响。

④上述暗目标法与深蓝算法通用性差。

⑤使用波段相关关系反演误差较大。

发明内容

本发明的目的是为解决上述问题，而提供一种针对缺少短波红外数据的气溶胶遥感反演方法及系统。

一种针对缺少短波红外数据的气溶胶遥感反演方法，它包括：

S1：卫星影像以及入射/观测角度数据获取；

S2：气溶胶光学厚度查找表的建立；

S3：构建波段间随地物类型和角度信息变化的回归关系；

S4：使用在S1、S2、S3步骤中的卫星影像数据、角度影像数据、气溶胶光学厚度查找表、不同地物类型和角度信息的回归关系计算目标影像对应的地表反射率红蓝波段关系下的气溶胶光学厚度；

S5：通过S4步骤得到的气溶胶光学厚度进行卫星影像的大气校正，得到的地表反射率必须通过ACI指数和NDVI年际变化函数两个筛选条件的匹配才能得到准确的气溶胶光学厚度影像；

S1所述的数据获取包括：获取光学卫星的高分辨率影像，影像至少包含蓝、绿、红、近红外四个波段的大气层顶反射率；以及通过影像自带的经纬度信息以及观测时间，计算得到卫星影像对应的太阳天顶角、观测天顶角、太阳方位角、观测方位角影像；

对S1获取的数据进行数据预处理，所述的预处理包括：对卫星影像进行定标、裁剪、云检测、计算大气层顶反射率，对角度信息进行三角网插值；

所述的预处理为：

在环境卫星影像查找影像自带的XML文件中的定标参数(g,L0)通过下式将影像的DN值计算为辐亮度L，单位为W*m^(-2)*sr^(-1)*um^(-1)；

L＝DN/g+L₀ (1)

使用shp格式的矢量文件对定标影像进行裁剪，得到所需范围内的影像数据，由于云的影响，无法从影像中计算后续步骤，因此需要将影像中云覆盖的区域去除，根据环境卫星的波谱特征，采用特征阈值的方法识别云覆盖像元，特征阈值算法分为两步，第一步使用环境卫星中的第三红光波段B₃设定阈值C1，由：

B₃>C1 (2)

C1取0.2，在此基础上进行第二步定义NDVI：

NDVI＝(B₄-B₃)/(B₄+B₃)>C2 (3)

C2取-0.015；B₃和B₄分别为第三红光波段和第四近红外波段表观反射率值；将云检测后的辐亮度L数据通过下式转换为表观反射率：

上式中，ρ^*为表观反射率，L_λ为波段辐亮度，d为天文单位日地距离，θ_s为太阳天顶角，E_sunλ为波段太阳表观辐射率均值，单位为W*m^(-2)*sr^(-1)*um^(-1)；

S1中还需生成将环境卫星自带的角度文件转换为栅格数据，其中太阳天顶角和太阳方位角经过像素点的经纬度和观测数据计算，观测天顶角和观测方位角使用角度文件将离散像素点使用三角网插值得到角度影像信息；

S2中所述的气溶胶光学厚度查找表的建立，包括：通过6S辐射传输模型程序计算在不同太阳观测天顶角和相对方位角、不同气溶胶光学厚度下的地表反射率以及大气层顶反射率，进而计算得到气溶胶光学厚度查找表中的ρ_a、S、T三个参数，其中太阳天顶角范围为0～66度以6度为间隔，观测天顶角范围为0～39度以3度为间隔，相对方位角范围为0～180度以12度为间隔，气溶胶类型选择大陆型；

所述的ρ_a为大气程辐射，S为大气向下反射率，T为大气透射率；

所述的ρ_a、S、T三个参数计算方法如下：

假设地表朗伯，大气水平均一，卫星传感器接收到的辐射能量来自地面辐射和大气程辐射两部分贡献，其传感器接收的表观反射率ρ^*表示为：

式中θ_s和θ_v为太阳天顶角与观测天顶角，ρ_a为大气程辐射，ρ_t为地表反射率，S为大气向下反射率，T为大气透射率。

令T＝T(θ_s)T(θ_v)，上式可以简写为：

式中三个参数ρ_a、S、T可以通过假设存在三个地表真实反射率分别为0、0.5、0.8，以及确定的观测几何和大气状态参数情况下使用6S进行辐射传输模拟得到三个表观反射率ρ₁、ρ₂、ρ₃，建立三个方程组即可求得三个参数ρ_a、S、T的唯一解。

在S3步骤中基于POLDER多角度数据集统计得到分角度、分土地覆盖类型的地表反射率波段间回归系数，对应角度信息以及地物类型得到地表红蓝波段间的回归系数查找表；

所述的回归系数(查找表)，由下述方法计算：

基于USGS波谱库的统计分析，可以发现地物暗目标的可见光波段间具有较好的相关关系，可以建立蓝光与绿光波段之间的线性回归关系查找表，

B₂＝a(theta,cover types)B₁+b(theta,cover types) (8)

为了评价回归系数的在不同角度和地类的适用性，我们通过真实值与回归模型预测值的差异来分析，将上式改写为：

Diff＝B₂-(B₁·a+b) (9)

式中：B₁，B₂分别为HJ-1A/B CCD数据的蓝光和绿光波段，a，b分别为回归系数的斜率和截距；

在S4步骤中，将上述处理后的数据，以及生成的气溶胶光学厚度查找表和回归系数查找表带回式5、7、9，得到初始气溶胶光学厚度；

在S5步骤中，因为计算气溶胶需要较高的波段间回归参数拟合精度，所以，气溶胶光学厚度的反演需要极高的波段回归精度，但是在应用中我们使用的是所有数据拟合的一组回归参数，在气溶胶光学厚度反演计算中所需回归参数与拟合直线垂直距离越远反演误差越大，而垂直距离直接改变截距大小，所以在回归系数中加入了对截距b的偏移因子i(以0.005作为间隔，bi＝b+(i*0.005-0.1),i＝0～40)，计算所有偏移条件下的气溶胶光学厚度值。

由于在截距b中加入了偏移因子i，可以计算得到多个气溶胶光学厚度值，因此，本发明定义筛选条件作为筛选气溶胶光学厚度指标：气溶胶校正指数(Aerosol calibrationindex,ACI)，如下式所示，

式中：B₄,B₁为加入偏移因子i反演的气溶胶光学厚度进行大气校正后的近红外波段和蓝光波段

该式使用不同偏移因子计算的气溶胶光学厚度进行大气校正后的近红外和蓝光波段的地表反射率，可以得到不同偏移因子对应是ACI值；ACI呈现先增大后减小的趋势，而在ACI的极大值之前选用年NDVI变化函数作为筛选条件，进一步确定气溶胶光学厚度，NDVI年际变化曲线为：

y＝-7.43e-6·x²+0.00298·x+0.32 (11)

式中：x为时间，y为NDVI。随着气溶胶光学厚度的变化大气校正后的NDVI也随之改变，取与先验NDVI最接近一组NDVI计算值对应的气溶胶光学厚度。

本发明通过上述流程方法提供一种针对缺少短波红外数据的气溶胶遥感反演系统，主要包括四个模块：

①数据预处理模块，用于对卫星影像的下载定标，裁剪，云检测，计算大气层顶反射率，以及对应的太阳观测角度数据。

②波段间回归关系计算模块，在红蓝波段通过分地物类型，分角度格网进行波段间回归计算，得到回归系数查找表，一次计算得到一个文本文件，后续反演气溶胶光学厚度的过程中无需再次计算，直接查表即可，提高系统运行效率。

③反演气溶胶光学厚度模块，使用S2步骤的计算方法，得到不同角度下气溶胶光学厚度查找表，也是一个文本文件，使用方法同波段间回归关系计算模块，直接查表即可。

④气溶胶光学厚度筛选模块，使用S5步骤的计算方法，通过多次迭代得到气溶胶光学厚度的精确结果，直接输出为二进制影像文件，文件扩展名为.dat，可用ENVI，ARCGIS软件直接打开使用，具有同输入的影像数据相同的地图投影坐标系统。

利用上述根据本发明提供的一种针对缺少短波红外数据的气溶胶遥感反演系统，基于AERONET(Aerosol Eobotic NETwork)网站下载的北京地区四个站点数据进行验证，利用相关的分类原则，通过实验观测结果以及站点观测结果进行匹配验证，采用均方根误差RMSE(Root Mean Square Error)、平均绝对误差MAE(Mean Absolute Error)及NASA预期误差EE(Expected Error)四种统计指标，其中EE为反演结果在误差线(下式)内的数据占总数的百分比。

EE＝τ_i±(0.05+0.20τ_i)

式中，τ_i为AREONET站点实测值。

经过验证得到反演气溶胶光学厚度结果的精确程度量化指标，为气溶胶光学厚度的遥感定量估算提供必要的前提。

本发明提出的针对缺少短波红外数据的气溶胶遥感反演方法及系统可以解决以下几个方面的问题：

①可以处理缺少短波红外波段的数据，只需要蓝、绿、红、近红外波段的数据。

②去除入射、观测角度对反演结果的影响。

③不同的地物类型给出不同的反演参数，提高算法普适性与通用性。

④提高气溶胶光学厚度反演精度。

基于上述过程使用IDL计算机程序语言开发的针对缺少短波红外数据的气溶胶遥感反演系统，其详细计算过程在具体实施方式中指出，方法流程以及系统模块之间关系和反演的部分气溶胶结果示例在附图中给出。

本发明包含建立气溶胶光学厚度查找表和波段间回归关系查找表的文本文件，计算在所有地物覆盖以及所有太阳观测角度下不同气溶胶光学厚度反演参数，在系统运行中直接查找对应的反演参数，提高系统运行效率，另外，通过筛选条件限制得到更精确的气溶胶反演结果。

附图说明

图1为针对缺少短波红外数据的气溶胶遥感反演方法流程图；

图2为针对缺少短波红外数据的气溶胶遥感反演系统部署结构图；

图3为HJ-1B CCD2影像数据及气溶胶光学厚度结果图(影像文件名：HJ1B-CCD2-457-68-20140726-L20001182135，成像日期：20140726)；

图4为HJ-1A CCD1影像数据及气溶胶光学厚度结果图(影像文件名：HJ1A-CCD1-3-64-20121006-L20000860427，成像日期：20121006)；

具体实施方式：

下面结合附图对本发明作详细的描述：

本发明中使用的是环境一号卫星影像数据，全称为环境减灾小卫星星座(HJ-1A/B/C)包括HJ-1A、B和C星三颗卫星，是中国用于环境和灾害监测预报的卫星，单台卫星的重访周期为4天，HJ-1A和HJ-1B两台卫星组网后重访周期可达到2～4天；其中HJ-1A/B每颗卫星上分别搭载了2台多光谱CCD传感器，包含蓝、绿、红和近红外4个波段，空间分辨率为30米；实验数据下载自中国资源卫星应用中心陆地观测卫星数据服务平台；环境卫星数据只是本发明的实验数据，不限制本发明的使用范围，其他类型数据只要满足数据输入的条件即可使用；

地物暗目标在可见光波段存在着显著的相关关系，可以基于这种相关关系通过气溶胶查找表来估算气溶胶光学厚度。

在本发明中，首先基于USGS光谱库确定暗目标在可见光波段之间相关关系，统计选取相关性最高的光学波段建立线性回归关系；

但是这种波段的相关关系会随着观测角度和地类发生变化(如图1中S3)，因此需要基于BRDF数据集建立分角度、地类的回归关系，然后通过气溶胶查找表反演气溶胶光学厚度。

考虑到气溶胶光学厚度的反演结果对回归系数较为敏感，通过改变波段间回归系数截距b，划分不同的间隔，得到不同间隔下对应的气溶胶光学厚度结果(如图2中S5)，我们定义气溶胶指示指数(ACI)和NDVI年际变化函数作为筛选条件，通过多次迭代计算选择最精确的气溶胶光学厚度影像结果。

在本发明中，具体通过以下步骤进行实施：

首先，在图1S1中经过数据预处理对卫星影像进行定标，裁剪，云检测，以及计算影像的表观反射率；具体为在环境卫星影像查找影像自带的XML文件中的定标参数(g,L₀)通过下式将影像的DN值计算为辐亮度L，单位为W*m^(-2)*sr^(-1)*um^(-1)。

L＝DN/g+L₀ (1)

使用shp格式的矢量文件对定标影像进行裁剪，得到所需范围内的影像数据，由于云的影响，无法从影像中计算后续步骤，因此需要将影像中云覆盖的区域去除，根据环境卫星的波谱特征，采用特征阈值的方法识别云覆盖像元，特征阈值算法分为两步，第一步使用环境卫星中的第三红光波段B₃设定阈值C1，由：

B₃>C1 (2)

C1取0.2，在此基础上进行第二步定义NDVI：

NDVI＝(B₄-B₃)/(B₄+B₃)>C2 (3)

C2取-0.015；B₃和B₄分别为第三红光波段和第四近红外波段表观反射率值；将云检测后的辐亮度L数据通过下式转换为表观反射率：

上式中，ρ^*为表观反射率，L_λ为波段辐亮度，d为天文单位日地距离，θ_s为太阳天顶角，E_sunλ为波段太阳表观辐射率均值，单位为W*m^(-2)*sr^(-1)*um^(-1)。

另外，在图1S1中还需生成将环境卫星自带的角度文件转换为栅格数据，其中太阳天顶角和太阳方位角经过像素点的经纬度和观测数据计算，观测天顶角和观测方位角使用角度文件将离散像素点使用三角网插值得到角度影像信息。

其次，在图1S2步骤中，假设地表朗伯，大气水平均一，卫星传感器接收到的辐射能量来自地面辐射和大气程辐射两部分贡献，其传感器接收的表观反射率ρ^*表示为：

式中θ_s和θ_v为太阳天顶角与观测天顶角，ρ_a为大气程辐射，ρ_t为地表反射率，S为大气向下反射率,T为大气透射率。

令T＝T(θ_s)T(θ_v)，上式可以简写为：

式中三个参数ρ_a、S、T可以通过假设存在三个地表真实反射率分别为0、0.5、0.8，以及确定的观测几何和大气状态参数情况下使用6S进行辐射传输模拟得到三个表观反射率ρ₁、ρ₂、ρ₃，建立三个方程组(下式)即可求得三个参数ρ_a、S、T的唯一解。

环境一号卫星观测天顶角最大角度在35度左右，综合考虑计算精度和计算机计算效率，我们将太阳天顶角、观测天顶角、相对方位角分别以6度、3度、12度为间隔划分，使用6S重复以上过程模拟不同观测几何和大气状态情况下(气溶胶光学厚度以0.05为间隔划分)对应的ρ_a、S、T，建立气溶胶查找表，查找表生成后不再进行上述计算，在后续中直接查表计算后续步骤。

再次，在图1S3步骤中，基于USGS波谱库的统计分析，可以发现地物暗目标的可见光波段间具有较好的相关关系，可以建立蓝光与绿光波段之间的线性回归关系，

B₂＝a(theta,cover types)B₁+b(theta,cover types) (8)

为了评价回归系数的在不同角度和地类的适用性，我们通过真实值与回归模型预测值的差异来分析，将上式改写为：

Diff＝B₂-(B₁·a+b) (9)

式中：B₁，B₂分别为HJ-1A/B CCD数据的蓝光和绿光波段，a，b分别为回归系数的斜率和截距。

通过气溶胶查找表查找卫星像元对应的太阳天顶角、观测天顶角和相对方位角，可以确定在此观测几何下40种变化梯度的气溶胶浓度对应的ρ_a、S、T值；代入该像元的表观反射率计算得到40组各波段地表反射率，然后求解Diff＝0时一组地表反射率其对应的气溶胶光学厚度即为反演结果。

为了确定角度和地物对回归系数的影响，我们使用POLDER 3多角度观测数据集进行分析；POLDER产品的空间分辨率为(6×7km)，是目前能获得的最新多角度卫星遥感数据，有丰富的角度、光谱和极化信息，是多角度遥感的理想数据之一。

POLDER传感器的特点在于每一轨可获取多达16个角度的观测，最大观测角度达到60度，因此每月的合成数据集可以包含数百个观测，其中POLDER 3BRDF数据集包含经纬坐标，经过云去除、大气校正等处理得到所有轨道观测的地表反射率；POLDER BRDF数据文件均为ASCII文件，每个文件包含三行头文件，第1、2行头文件提供了像元经纬坐标、地物类别代码(GLC2000土地覆盖类型)、轨道号、有效观测次数(单次观测周期最多16个观测值)、以及该像元内地物的均匀性。

将PLODER BRDF数据集转换通过波段转换系数转换为HJ-1A/B CCD传感器下的BRDF数据集，之后分角度格网(太阳天顶角0～70°以2°为间隔，观测天顶角0～40°以2°为间隔，相对方位角以5°为间隔进行格网划分)分地物类别(GLC2000)对HJ-1A/B CCD传感器的蓝绿波段进行线性回归，建立回归系数查找表，回归系数查找表由太阳-观测角度、地物类型以及蓝绿波段之间的拟合斜率a和截距b构成，同气溶胶查找表生成后使用方式一样，在后续步骤中直接查表计算即可。

在图1S4步骤中，将上述处理后的数据，以及生成的气溶胶光学厚度查找表和回归系数查找表带回式5、7、9，得到初始气溶胶光学厚度。

在图1S5步骤中，因为计算气溶胶需要较高的波段间回归参数拟合精度，所以，气溶胶光学厚度的反演需要极高的波段回归精度，但是在应用中我们使用的是所有数据拟合的一组回归参数，在气溶胶光学厚度反演计算中所需回归参数与拟合直线垂直距离越远反演误差越大，而垂直距离直接改变截距大小，所以在回归系数中加入了对截距b的偏移因子i(以0.005作为间隔，bi＝b+(i*0.005-0.1),i＝0～40)，计算所有偏移条件下的气溶胶光学厚度值。

由于在截距b中加入了偏移因子i，可以计算得到多个气溶胶光学厚度值，因此，本发明定义筛选条件作为筛选气溶胶光学厚度指标：气溶胶校正指数(Aerosol calibrationindex,ACI)，如下式所示，

式中：B₄,B₁为加入偏移因子i反演的气溶胶光学厚度进行大气校正后的近红外波段和蓝光波段

该式使用不同偏移因子计算的气溶胶光学厚度进行大气校正后的近红外和蓝光波段的地表反射率，可以得到不同偏移因子对应是ACI值；ACI呈现先增大后减小的趋势，而在ACI的极大值之前选用年NDVI变化函数作为筛选条件，进一步确定气溶胶光学厚度，NDVI年际变化曲线为：

y＝-7.43e-6·x²+0.00298·x+0.32 (11)

式中：x为时间，y为NDVI；随着气溶胶光学厚度的变化大气校正后的NDVI也随之改变，取与先验NDVI最接近一组NDVI计算值对应的气溶胶光学厚度作为本发明的最终结果。

实施例

本发明涉及一种针对缺少短波红外数据的气溶胶遥感反演方法及系统，以下结合附图和实施例进一步详细说明本发明技术方案所涉及的各个细节；本实施例使用个人计算机(PC机)进行仿真实现，其软件基于64位Windows 10操作系统和IDL二次开发平台。

采用的实施例为在在中国资源卫星应用中心下载的环境卫星影像，根据技术方案(图1)编写IDL处理系统程序，系统程序结构如图2，调试运行程序得到气溶胶光学厚度结果。

编写反演气溶胶光学厚度程序：根据技术方案中的①卫星遥感图像数据预处理，包含定标、裁剪云检测、计算大气层顶反射率；②计算气溶胶光学厚度查找表；③建立波段间回归关系查找表；④反演气溶胶；⑤筛选气溶胶光学厚度影像，包含在四个模块中，如图2。

实施例中使用两幅影像文件名分别为：HJ1A-CCD1-3-64-20121006-L20000860427，HJ1B-CCD2-457-68-20140726-L20001182135；成像日期分别为2012年10月6日和2014年7月26日，卫星过境时间在当地时间10点左右，像素点行列数为16086×14627、15610×13978，WGS_1984_UTM_Zone_50N如图3、4左图，影像为4、3、2波段假彩色合成。

在影像自带压缩包中有与影像配套的角度文件，通过IDL编写的角度文件读取插值模块，生成与对应影像相同行列数相同坐标系的角度文件，包含：太阳天顶角、太阳方位角、观测天顶角、观测方位角。

对两幅影像进行定标，通过读取对应影像配套的XML文件中的定标系数，采用式1对整幅影像进行定标。

选择覆盖影像的部分矢量文件，使用IDL进行掩膜运算提取所需处理影像，如果不选，默认处理整幅影像。

由式2、3在IDL中使用波段运算函数进行影像的云检测。

之后通过式4计算得到影像表观反射率。

以上计算过程为数据预处理模块。

下面使用IDL编写运行气溶胶光学厚度模块：

在这一模块中需要通过IDL调用6S辐射传输程序，计算在大陆型气溶胶模式，大气模式选择中纬度夏季，通过循环计算在不同太阳天顶角、观测天顶角、相对方位角，以及不同的气溶胶光学厚度下使用式6、7得到对应不同角度不同气溶胶光学厚度对应的三个参数ρ_a、S、T，储存在txt文件中。

波段间回归关系模块：

读取PLODER多角度数据集，统计在不同地物类型中，使用核函数拟合地物的BRDF信息，通过式8、9在所有太阳观测角度下进行波段间回归，得到对应不同地物，不同太阳观测角度的回归参数a、b，储存在txt文件中。

反演以及筛选气溶胶光学厚度：

由式6、7反演气溶胶光学厚度，通过改变波段间回归关系的截距b值计算得到多幅对应的气溶胶光学厚度影像结果，采用式10、11进行对气溶胶光学厚度进行筛选，得到精确的气溶胶光学厚度结果，气溶胶光学厚度具有与卫星影像相同的行列数和投影坐标系；图3、4右图为使用本发明提出的方法经过运行得到对应左图卫星影像的精确气溶胶反演结果。

Claims (8)

1.一种针对缺少短波红外数据的气溶胶遥感反演方法，它包括：

S1：卫星影像以及入射/观测角度数据获取；

所述的数据获取包括：获取光学卫星的高分辨率影像，影像至少包含蓝、绿、红、近红外四个波段的大气层顶反射率；以及通过影像自带的经纬度信息以及观测时间，计算得到卫星影像对应的太阳天顶角、观测天顶角、太阳方位角、观测方位角影像；

对获取的数据进行数据预处理，所述的预处理包括：对卫星影像进行定标、裁剪、云检测、计算大气层顶反射率；生成将环境卫星自带的角度文件转换为栅格数据，其中太阳天顶角和太阳方位角经过像素点的经纬度和观测数据计算，观测天顶角和观测方位角使用角度文件将离散像素点使用三角网插值得到角度影像信息。

S2：气溶胶光学厚度查找表的建立；

S3：构建波段间随地物类型和角度信息变化的回归关系查找表；

S4：使用在S1、S2、S3步骤中的卫星影像数据、角度影像数据、气溶胶光学厚度查找表、不同地物类型和角度信息的回归关系计算目标影像对应的地表反射率红蓝波段关系下的气溶胶光学厚度；

S5：通过S4步骤得到的气溶胶光学厚度进行卫星影像的大气校正，得到的地表反射率必须通过ACI指数和NDVI年际变化函数两个筛选条件的匹配才能得到准确的气溶胶光学厚度影像。

2.根据权利要求1所述的一种针对缺少短波红外数据的气溶胶遥感反演方法，其特征在于：

所述的预处理为：

在环境卫星影像查找影像自带的XML文件中的定标参数(g,L0)通过下式将影像的DN值计算为辐亮度L，单位为W*m^(-2)*sr^(-1)*um^(-1)；

L＝DN/g+L₀ (1)

使用shp格式的矢量文件对定标影像进行裁剪，得到所需范围内的影像数据，由于云的影响，无法从影像中计算后续步骤，因此需要将影像中云覆盖的区域去除，根据环境卫星的波谱特征，采用特征阈值的方法识别云覆盖像元，特征阈值算法分为两步，第一步使用环境卫星中的第三红光波段B₃设定阈值C1，由：

B₃>C1 (2)

C1取0.2，在此基础上进行第二步定义NDVI：

NDVI＝(B₄-B₃)/(B₄+B₃)>C2 (3)

C2取-0.015；B₃和B₄分别为第三红光波段和第四近红外波段表观反射率值；将云检测后的辐亮度L数据通过下式转换为表观反射率：

上式中，ρ^*为表观反射率，L_λ为波段辐亮度，d为天文单位日地距离，θ_s为太阳天顶角，E_sunλ为波段太阳表观辐射率均值，单位为W*m^(-2)*sr^(-1)*um^(-1)。

3.根据权利要求2所述的一种针对缺少短波红外数据的气溶胶遥感反演方法，其特征在于：

S2中所述的气溶胶光学厚度查找表的建立，包括：通过6S辐射传输模型程序计算在不同太阳观测天顶角和相对方位角、不同气溶胶光学厚度下的地表反射率以及大气层顶反射率，进而计算得到气溶胶光学厚度查找表中的ρ_a、T三个参数，其中太阳天顶角范围为0～66度以6度为间隔，观测天顶角范围为0～39度以3度为间隔，相对方位角范围为0～180度以12度为间隔，气溶胶类型选择大陆型；

所述的ρ_a为大气程辐射，为大气向下反射率，T为大气透射率。

4.根据权利要求3所述的一种针对缺少短波红外数据的气溶胶遥感反演方法，其特征在于：

所述的ρ_a、T三个参数计算方法如下：

假设地表朗伯，大气水平均一，卫星传感器接收到的辐射能量来自地面辐射和大气程辐射两部分贡献，其传感器接收的表观反射率ρ^*表示为：

式中θ_s和θ_v为太阳天顶角与观测天顶角，ρ_a为大气程辐射，ρ_t为地表反射率，为大气向下反射率，T为大气透射率。

令T＝T(θ_s)T(θ_v)，上式可以简写为：

式中三个参数ρ_a、T可以通过假设存在三个地表真实反射率分别为0、0.5、0.8，以及确定的观测几何和大气状态参数情况下使用6S进行辐射传输模拟得到三个表观反射率ρ₁、ρ₂、ρ₃，建立三个方程组即可求得三个参数ρ_a、T的唯一解。

5.根据权利要求1、2、3或4所述的一种针对缺少短波红外数据的气溶胶遥感反演方法，其特征在于：

在S3步骤中基于POLDER多角度数据集统计得到分角度、分土地覆盖类型的地表反射率波段间回归系数，对应角度信息以及地物类型得到地表红蓝波段间的回归系数查找表；

所述的回归系数查找表，由下述方法计算：

基于USGS波谱库的统计分析，可以发现地物暗目标的可见光波段间具有较好的相关关系，可以建立蓝光与绿光波段之间的线性回归关系，

B₂＝a(theta,cover types)B₁+b(theta,cover types) (8)

为了评价回归系数的在不同角度和地类的适用性，我们通过真实值与回归模型预测值的差异来分析，将上式改写为：

Diff＝B₂-(B₁·a+b) (9)

式中：B₁，B₂分别为HJ-1A/B CCD数据的蓝光和绿光波段，a，b分别为回归系数的斜率和截距。

6.根据权利要求5所述的一种针对缺少短波红外数据的气溶胶遥感反演方法，其特征在于：

在S4步骤中，将s1预处理数据，以及生成的气溶胶光学厚度查找表和回归系数查找表带回式5、7、9，得到初始气溶胶光学厚度。

7.根据权利要求6所述的一种针对缺少短波红外数据的气溶胶遥感反演方法，其特征在于：

在S5步骤中，在回归系数中加入了对截距b的偏移因子i(以0.005作为间隔，bi＝b+(i*0.005-0.1),i＝0～40)，计算所有偏移条件下的气溶胶光学厚度值。

由于在截距b中加入了偏移因子i，可以计算得到多个气溶胶光学厚度值，因此，本发明定义筛选条件作为筛选气溶胶光学厚度指标：气溶胶校正指数ACI，如下式所示，

式中：B₄,B₁为加入偏移因子i反演的气溶胶光学厚度进行大气校正后的近红外波段和蓝光波段

该式使用不同偏移因子计算的气溶胶光学厚度进行大气校正后的近红外和蓝光波段的地表反射率，可以得到不同偏移因子对应是ACI值；ACI呈现先增大后减小的趋势，而在ACI的极大值之前选用年NDVI变化函数作为筛选条件，进一步确定气溶胶光学厚度，NDVI年际变化曲线为：

y＝-7.43e-6·x²+0.00298·x+0.32 (11)

式中：x为时间，y为NDVI。随着气溶胶光学厚度的变化大气校正后的NDVI也随之改变，取与先验NDVI最接近一组NDVI计算值对应的气溶胶光学厚度。

8.一种针对缺少短波红外数据的气溶胶遥感反演系统，主要包括四个模块：

①数据预处理模块，用于对卫星影像的下载定标，裁剪，云检测，计算大气层顶反射率，以及对应的太阳观测角度数据。

②波段间回归关系计算模块，在红蓝波段通过分地物类型，分角度格网进行波段间回归计算，得到回归系数查找表，一次计算得到一个文本文件，后续反演气溶胶光学厚度的过程中无需再次计算，直接查表即可，提高系统运行效率。

③反演气溶胶光学厚度模块，使用权利要求1所述的一种针对缺少短波红外数据的气溶胶遥感反演方法S2步骤的计算方法，得到不同角度下气溶胶光学厚度查找表，也是一个文本文件，使用方法同波段间回归关系计算模块，直接查表即可。

④气溶胶光学厚度筛选模块，使用权利要求1所述的一种针对缺少短波红外数据的气溶胶遥感反演方法S5步骤的计算方法，通过多次迭代得到气溶胶光学厚度的精确结果，直接输出为二进制影像文件，文件扩展名为.dat，可用ENVI，ARCGIS软件直接打开使用，具有同输入的影像数据相同的地图投影坐标系统。