CN116297068A

CN116297068A - 基于地表反射率优化的气溶胶光学厚度反演方法及系统

Info

Publication number: CN116297068A
Application number: CN202310593490.XA
Authority: CN
Inventors: 康苒; 徐天河; 司璐璐; 尹会英; 胡松伯; 邓彩云; 陈婷婷
Original assignee: Shandong University
Current assignee: Shandong University
Priority date: 2023-05-25
Filing date: 2023-05-25
Publication date: 2023-06-23
Anticipated expiration: 2043-05-25
Also published as: CN116297068B

Abstract

本发明公开一种基于地表反射率优化的气溶胶光学厚度反演方法及系统，涉及气溶胶光学厚度数据识别技术领域，包括：确定气溶胶光学厚度最小值的取值范围，根据取值范围内的气溶胶遥感产品，确定采集时间范围；根据采集时间范围内的卫星遥感影像，确定最小表观反射率的卫星遥感影像；根据卫星过境时间获取实测地基遥感影像；对确定的卫星遥感影像选择蓝光波段、红光波段、近红外波段和短波红外波段在设定像元范围内的值，构建地表反射率库和表观反射率库；根据观测角度、地理位置、大气参数、植被指数、表观反射率和地表反射率，得到气溶胶光学厚度实测值，优化地表反射率库的同时，提高反演精度。

Description

基于地表反射率优化的气溶胶光学厚度反演方法及系统

技术领域

本发明涉及气溶胶光学厚度数据识别技术领域，特别是涉及一种基于地表反射率优化的气溶胶光学厚度反演方法及系统。

背景技术

气溶胶是地球-大气-海洋系统的重要成分，气溶胶对空气质量有很大影响。气溶胶光学厚度（Aerosol Optical Depth，AOD）描述了大气中气溶胶总柱对阳光的衰减，是气溶胶遥感最重要的光学性质之一，通常以此推算大气中的气溶胶含量，由此进行大气环境质量的监测分析。

表观反射率以及地表反射率是影响AOD测算的两个关键指标，在城市等地表类型比较复杂的地区准确估算AOD，需要在计算大气程辐射的时候把地表反射率和表观反射率准确地分离，否则，地表反射率的误差会被放大10倍传递给AOD反演结果中。目前，已有研究大多基于低空间分辨率地表反射率产品来解决地表反射率或表观反射率估计问题，由于云的影响、卫星传感器问题仍然面临着地表反射率估算精度不高、数据缺口等诸多挑战。

目前，遥感监测AOD的方法分为地基监测和卫星遥感监测；对于地基监测，存在AOD站点分布不均且局部站点个数少等问题，降低了算法的通用性，特别是无地基站点区域，无法满足AOD高时空分辨率反演精度要求。

对于卫星遥感监测，传统研究方法通常基于物理模型，即通过模拟大气层顶（Topof Atmosphere, TOA）反射率（亦称表观反射率）的辐射传输模型，利用查找表模拟辐亮度与卫星观测辐亮度进行比较来反演气溶胶特性。该方法费时、费力，且需要建立详细的查找表数据库。随着机器学习方法进入大气科学领域，卫星数据被应用大尺度的AOD反演中，AOD定量遥感反演数据的预测精度、分辨率和覆盖范围得到了很大提高。考虑到在小空间尺度上人口密度和道路密度较高的地区，该区域上空气溶胶时空分布存在差异，仍需要建立综合多因素的深度学习模型以获取较高的AOD时空分布结果。

综上所述，目前AOD的遥感反演仍存在以下问题：随着卫星遥感数据量增大，如何选取“最干净”的影像时间构建精确的地表反射率库，以及如何利用深度学习实现少站点、高精度的AOD精确估算。

发明内容

为了解决上述问题，本发明提出了一种基于地表反射率优化的气溶胶光学厚度反演方法及系统，构建了精确的地表反射率库，引入能够表征地表参数的重要波段及比值，实现地表反射率库优化的同时，提高AOD反演精度。

为了实现上述目的，本发明采用如下技术方案：

第一方面，本发明提供一种基于地表反射率优化的气溶胶光学厚度反演方法，包括：

根据区域内对应的气溶胶遥感产品，确定气溶胶光学厚度最小值的取值范围，根据处于取值范围内的气溶胶遥感产品，确定采集时间范围；

根据采集时间范围内的卫星遥感影像，确定具有最小表观反射率的卫星遥感影像；

对具有最小表观反射率的卫星遥感影像选择蓝光波段、红光波段、近红外波段和短波红外波段在设定像元范围内的值，以构建地表反射率库和表观反射率库；

获取历史卫星遥感影像及对应的实测地基遥感影像，并基于地表反射率库和表观反射率库确定对应的表观反射率和地表反射率，再结合卫星观测角度、地理位置、大气参数和植被指数，对构建的反演神经网络进行训练，以此得到待测区域的气溶胶光学厚度实测值。

作为可选择的实施方式，对具有最小表观反射率的卫星遥感影像，通过卫星自带的质量波段进行按位运算去除云或云阴，并对去云处理后的缺失数据进行填补。

作为可选择的实施方式，对具有最小表观反射率的卫星遥感影像选择蓝光波段、红光波段和近红外波段在20×20像元范围内的波段值，以及短波红外波段在10×10像元范围内的波段值，并将各个波段选择的波段值去除异常值后求平均。

作为可选择的实施方式，获取卫星过境时间的前后设定时间段内的实测地基遥感影像，并对实测地基遥感影像内的值求平均。

作为可选择的实施方式，所述地表反射率库包括红光波段、近红外波段、短波红外波段、红光波段与近红外波段的比值、红光波段与短波红外波段的比值、近红外波段与短波红外波段的比值。

作为可选择的实施方式，所述表观反射率包括蓝光波段、红光波段和近红外波段。

作为可选择的实施方式，所述大气参数包括臭氧和水汽，所述反演神经网络以观测角度、六种地表反射率、三种表观反射率、臭氧、水汽、地理位置和植被指数为输入，以实测地基遥感影像中气溶胶光学厚度实测值为输出，且反演神经网络包括一个含有14个节点的输入层，2个隐层，且每个隐层各有12个节点以及含有1个节点的输出层。

第二方面，本发明提供一种基于地表反射率优化的气溶胶光学厚度反演系统，包括：

采集时间确定模块，被配置为根据区域内对应的气溶胶遥感产品，确定气溶胶光学厚度最小值的取值范围，根据处于取值范围内的气溶胶遥感产品，确定采集时间范围；

最小表观反射率确定模块，被配置为根据采集时间范围内的卫星遥感影像，确定具有最小表观反射率的卫星遥感影像；

反射率库构建模块，被配置为对具有最小表观反射率的卫星遥感影像选择蓝光波段、红光波段、近红外波段和短波红外波段在设定像元范围内的值，以构建地表反射率库和表观反射率库；

反演模块，被配置为获取历史卫星遥感影像及对应的实测地基遥感影像，并基于地表反射率库和表观反射率库确定对应的表观反射率和地表反射率，再结合卫星观测角度、地理位置、大气参数和植被指数，对构建的反演神经网络进行训练，以此得到待测区域的气溶胶光学厚度实测值。

第三方面，本发明提供一种电子设备，包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令，所述计算机指令被处理器运行时，完成第一方面所述的方法。

第四方面，本发明提供一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时，完成第一方面所述的方法。

与现有技术相比，本发明的有益效果为：

本发明构建了精确的地表反射率库，以实现AOD的精确反演，具体通过确定气溶胶光学厚度最小值的取值范围，并根据处于取值范围内的气溶胶遥感产品，确定采集时间范围，可剔除受到大气影响较多的影像，以确定可用高分辨率的卫星遥感影像；再结合采集时间范围内的卫星遥感影像，确定具有最小表观反射率的高分辨率的卫星遥感影像，最后与站点实测的地基遥感影像进行精确时空匹配，使得AOD精确反演重要参数之一的地表反射率库更加精确。

本发明提出一种基于地表反射率优化的气溶胶光学厚度反演方法及系统，除了引入植被指数、大气参数等信息，还引入了能够表征地表参数的重要波段及比值，实现地表反射率库优化的同时，提高了AOD反演的精度。

本发明的气溶胶光学厚度反演方法对于少点站数据、省级尺度的AOD精确反演具有可复制性，能够有效解决地基站点分布不均、布设较少等问题，为AOD的高精度时空反演提供了重要参考手段。

本发明附加方面的优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

图1为本发明实施例1提供的基于地表反射率优化的气溶胶光学厚度反演方法流程图；

图2为本发明实施例1提供的地表反射率构建流程示意图；

图3为本发明实施例1提供的BP神经网络结构图；

图4（a）-图4（b）为本发明实施例1提供的训练集预测结果和测试集预测结果对比图。

具体实施方式

下面结合附图与实施例对本发明做进一步说明。

应该指出，以下详细说明都是示例性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

实施例1

本实施例提供一种基于地表反射率优化的气溶胶光学厚度反演方法，如图1所示，包括：

在本实施例中，基于相同时期或连续一段时间内的地表反射率不变，表观反射率最小的像元受到大气影响很弱的原理，采用MODIS AOD产品中的影像数据，以确定气溶胶光学厚度最小值的取值范围及确定采集时间范围；然后获取采集时间范围内的卫星遥感影像，从而确定具有最小表观反射率的卫星遥感影像；最后基于该卫星遥感影像构建地表反射率库。

下面结合图2对地表反射率库的构建进行详细说明。

为构建时空连续的地表反射率库，本实施例将时间窗口设置为30天，提取历史时间（如2018-2020年）区域内AOD站点对应的气溶胶遥感产品（本实施例采用MODIS影像数据生成的气溶胶遥感产品），由此可得到AOD取值分布范围，继而确定AOD最小值的取值范围，如AOD最小值集中在0.02-0.1之间，那么AOD值大于0.1的影像日期则代表该天受气溶胶污染影响大；

根据气溶胶光学厚度最小值的取值范围，筛选处于取值范围内的气溶胶遥感产品，并确定采集时间范围，在该时间范围内则表示该时期影像受大气影响较弱，更有利于获取真实的地表反射率；

确定采集时间范围后，获取该采集时间范围内的卫星遥感影像（本实施例采用Sentinel-2卫星遥感影像），确定最小表观反射率对应的影像时间，得到该时间窗口内“最干净”的卫星遥感影像；并对确定的卫星遥感影像，通过Sentinel-2卫星自带的质量波段QA60进行按位运算去除云（Bit 11）或者云阴（Bit 10）；去云处理后的缺失数据选用时间相近且云量小于20%的Landsat-8影像地表反射率数据进行填补，以实现无缝隙AOD空间反演。

在本实施例中，区域内地基站点采用的是SONET站点进行观测，由于Sentinel-2卫星观测与SONET站点观测的时空尺度不一致，SONET站点观测的数据是按照15min时间间隔获得点状连续观测，而遥感卫星是覆盖区域的面状瞬时观测；此外，Sentinl-2卫星的空间分辨率最高为10m，过于精细的空间分辨率容易使模型受到地形起伏和地表覆盖的影响，降低信噪比，给AOD反演结果带来误差；如果空间采样窗口过小或为中心像元时，在减少数据量的同时也会出现畸值，但是如果空间采样窗口太大，则无法真实地反映气溶胶状况。

因此，为了增强数据可信度，平衡空间异质性和数据细节，本实施例采用时空匹配方法，且时间匹配考虑到研究范围为区域尺度；

因此，对于空间匹配，以站点经纬度为中心，对处理后的具有最小表观反射率的卫星遥感影像选择蓝光波段（B2）、红光波段（B4）和近红外波段（B8）在20×20像元范围内的值，以及短波红外波段（B12）在10×10像元范围内的波段值，统一到200m×200m范围内采样，并将各个波段选择的波段值去除异常值后求平均，从而构建地表反射率库和表观反射率库；

已有研究表明，短波红外波段2.12um波段受大气的影响远小于较可见光和近红外波段受大气的影响，早期的地表反射率估算方法也将2.12um波段的表观反射率近似为地表反射率。因此，为了能从Sentinel-2卫星遥感影像中精确估计地表反射率，本实施例考虑引入短波红外波段（B12），由于在红、蓝波段的反射率和短波红外的2.12um波段间的反射率存在稳定关系，所以可以合成多天的晴空像元来构建地表反射率库；由此，将B4、B8、B12三个不同波段进行组合计算，生成 B4、B8、B12、

、/>

、/>

，以此构建地表反射率库；对于表观反射率库基于B2、B4和B12构建。

对于时间匹配，本实施例获取历史卫星遥感影像及对应的实测地基遥感影像，其中根据卫星过境时间确定区域内地基站点的采集时间范围，以此获取实测地基遥感影像；具体是以卫星过境时间的前后设定时间段（本实施例采用15 min）内的SONET站点数据求平均，以此是为了获取卫星遥感影像对应的真实AOD值，为了训练反演神经网络。

由此，基于上述两部分（即Sentinel-2卫星遥感影像中的数据和SONET站点的实测数据）实现了时空匹配。

在本实施例中，构建地表反射率库后，基于库中被优化的地表反射率，再结合观测角度、地理位置、大气参数、植被指数、表观反射率等，可对AOD实测值进行反演。

本实施例的反演神经网络基于BP神经网络构建，主要包括设置层之间的节点联接权重及偏置、隐含层数及神经元节点数，明确各层的传递函数类型、确定训练参数，训练参数包括观测角度、大气参数、地理位置、地表反射率、表观反射率、植被指数；

其中，观测角度为卫星的观测天顶角；大气参数为臭氧、水汽，用于大气分子吸收性校正；地理位置为经纬度，用于确定最小反演像元位置；地表反射率为地表反射率库中的信息，即B4、B8、B12、

、/>

、/>

；表观反射率基于大气敏感波段构建，即B2、B4、B12；对于植被指数的研究，通常采用红外波段和近红外波段计算，但这种计算方式受气溶胶影响很大，由于在红、蓝波段的反射率和短波红外的2.13 um波段间的反射率存在稳定关系，所以本实施例利用短波红外波段计算归一化植被指数筛选暗像元，以此得到植被指数

，/>

和/>

分别为近红外波段和短波红外波段的表观反射率。

本实施例从影响气溶胶反演精度的三个重要方面出发，包括大气参数、卫星参数及地表信息，重点刻画地表参数信息变量，引入精确地表信息的神经网络模型能够减少输入参数的同时，得到较高精度的区域尺度AOD反演结果。

在构建反演神经网络之前，为了提高结果合理性和解决多类别样本数据值域差别大的问题，需要对数据进行归一化处理，即把不同种类的数据统一到0-1的范围，消除数据量纲，以加快BP网络的运算速度，归一化公式为：

（1）

式中，

和/>

分别为归一化之后和归一化之前的变量值，/>

和/>

分别为变量/>

的最大值和最小值。

设BP神经网络的输入层、隐藏层和输出层分别有

、/>

、/>

个神经元，各层的输入输出关系为：

输入层：

，/>

（2）

隐含层：

，/>

（3）

输出层：

，/>

（4）

其中，

表示输入层/>

的第/>

个节点的输入变量，/>

为偏差，/>

为激活函数；/>

是线性传递函数；/>

表示第k个节点的输出。

在连接权值更新阶段，为了使均方误差

的值最小化，BP神经网络基于梯度下降策略，以目标的负梯度方向对参数进行调整，逐步修正输出层与隐藏层之间的连接权值和偏差，偏差的更新步骤与连接权值的更新步骤类似，具体为：

误差性能函数为：

（5）

输出层：

（6）

隐含层：

（7）

其中，

为真值，/>

和/>

分别为神经元的梯度项，/>

和/>

为学习速率；/>

表示输入层/>

的第i个节点的输入变量，/>

，n为输入变量的个数，本实施例设为14，分别代表观测角度、地表反射率（B4、B8、B12、/>

、/>

、/>

）、表观反射率（B2、B4、B12）、臭氧、水汽、/>

、经纬度；/>

表示第一个隐含层/>

中的第/>

个节点与输入层中的第/>

个节点之间的连接权值，/>

表示更新后的连接权值；/>

表示第二个隐含层中第/>

个节点和输入层中的第/>

个节点之间的连接权值，/>

表示更新后的连接权值；/>

表示输出层节点和第二个隐含层中的第/>

个节点之间的连接权值；/>

为激活函数的导函数；/>

是线性传递函数的导函数。

在本实施例中，隐含层的激活函数采用双曲正切s型激活函数tanh；将观测角度、地表反射率(B4、B8、B12、

、/>

、/>

)、表观反射率（B2、B4、B12）、臭氧、水、/>

、经纬度作为输入，将AOD作为输出，由此可以确定BP神经网络输入层有14个节点，输出层1个节点；通过反复测试，当最佳隐含层节点数为12时，收敛速度最快；同时相关研究表明，增加隐含层个数，可以提高仿真精度，因此，本实施例设计的BP神经网络有2个隐层，每个隐层各有12个节点，由此可得BP神经网络的结构为14-12-12-1，如图3所示，其中，/>

表示第一个隐含层中各节点的值，/>

表示第二隐含层中的各节点的值，/>

代表输出层节点的值。

在本实施例中，输入层到隐含层的激活函数采用relu函数，输出层的传递函数采用线性激活函数；学习率设置为0.01，迭代次数为1000；通过以上步骤训练样本数据库，获得各层节点间的联接权重和各节点的偏置。

在本实施例中，神经网络精度评价采用相对分析误差RPD，表示预测能力：

（8）

其中，STDEV为标准差，RMSE均方根误差的比值；当1.5<RPD<2.0时，表明神经网络只能进行粗略估测，当2.0≤RPD<2.5时，表明神经网络具有较好的预测能力，当2.5≤RPD<3.0时，神经网络具有很好的预测能力，当RPD≥3.0时，神经网络具有极好的预测能力。

以无间隙高分辨率大气污染物浓度数据集（the Long-term Gap-free High-resolution Air Pollutant concentration dataset, LGHAP），其时间分辨率为1天，空间分辨率为1 km，与地基AOD观测数据具有良好的一致性，将其作为反演结果缺失部分的弥补数据协同地基站点作为验证数据。

数据训练集为站点时空匹配数据，验证集样本为生成的随机抽样点，训练集与验证集比例为3:2，最终神经网络训练结果为RPD=3.23、RMSE=0.12，预测结果RPD=3.12、RMSE=0.16，如图4（a）-图4（b）所示，表明具有极好的预测能力。

实施例2

本实施例提供一种基于地表反射率优化的气溶胶光学厚度反演系统，包括：

此处需要说明的是，上述模块对应于实施例1中所述的步骤，上述模块与对应的步骤所实现的示例和应用场景相同，但不限于上述实施例1所公开的内容。需要说明的是，上述模块作为系统的一部分可以在诸如一组计算机可执行指令的计算机系统中执行。

在更多实施例中，还提供：

一种电子设备，包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令，所述计算机指令被处理器运行时，完成实施例1中所述的方法。为了简洁，在此不再赘述。

应理解，本实施例中，处理器可以是中央处理单元CPU，处理器还可以是其他通用处理器、数字信号处理器DSP、专用集成电路ASIC，现成可编程门阵列FPGA或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。

存储器可以包括只读存储器和随机存取存储器，并向处理器提供指令和数据、存储器的一部分还可以包括非易失性随机存储器。例如，存储器还可以存储设备类型的信息。

一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时，完成实施例1中所述的方法。

实施例1中的方法可以直接体现为硬件处理器执行完成，或者用处理器中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器、闪存、只读存储器、可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。该存储介质位于存储器，处理器读取存储器中的信息，结合其硬件完成上述方法的步骤。为避免重复，这里不再详细描述。

本领域普通技术人员可以意识到，结合本实施例描述的各示例的单元即算法步骤，能够以电子硬件或者计算机软件和电子硬件的结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims

1.基于地表反射率优化的气溶胶光学厚度反演方法，其特征在于，包括：

2.如权利要求1所述的基于地表反射率优化的气溶胶光学厚度反演方法，其特征在于，对具有最小表观反射率的卫星遥感影像，通过卫星自带的质量波段进行按位运算去除云或云阴，并对去云处理后的缺失数据进行填补。

3.如权利要求1所述的基于地表反射率优化的气溶胶光学厚度反演方法，其特征在于，对具有最小表观反射率的卫星遥感影像选择蓝光波段、红光波段和近红外波段在20×20像元范围内的波段值，以及短波红外波段在10×10像元范围内的波段值，并将各个波段选择的波段值去除异常值后求平均。

4.如权利要求1所述的基于地表反射率优化的气溶胶光学厚度反演方法，其特征在于，获取卫星过境时间的前后设定时间段内的实测地基遥感影像，并对实测地基遥感影像内的值求平均。

5.如权利要求1所述的基于地表反射率优化的气溶胶光学厚度反演方法，其特征在于，所述地表反射率库包括红光波段、近红外波段、短波红外波段、红光波段与近红外波段的比值、红光波段与短波红外波段的比值、近红外波段与短波红外波段的比值。

6.如权利要求1所述的基于地表反射率优化的气溶胶光学厚度反演方法，其特征在于，所述表观反射率包括蓝光波段、红光波段和近红外波段。

7.如权利要求5或6所述的基于地表反射率优化的气溶胶光学厚度反演方法，其特征在于，所述大气参数包括臭氧和水汽，所述反演神经网络以观测角度、六种地表反射率、三种表观反射率、臭氧、水汽、地理位置和植被指数为输入，以实测地基遥感影像中气溶胶光学厚度实测值为输出，且反演神经网络包括一个含有14个节点的输入层，2个隐层，且每个隐层各有12个节点以及含有1个节点的输出层。

8.基于地表反射率优化的气溶胶光学厚度反演方法，其特征在于，包括：

9.一种电子设备，其特征在于，包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令，所述计算机指令被处理器运行时，完成权利要求1-7任一项所述的方法。

10.一种计算机可读存储介质，其特征在于，用于存储计算机指令，所述计算机指令被处理器执行时，完成权利要求1-7任一项所述的方法。