基于深度学习反演区域气溶胶光学厚度 
下载： 1262次
1 引言
大气气溶胶通过辐射强迫和化学扰动在地球辐射收支中发挥重要作用。大气气溶胶粒子通过对太阳辐射和地球长波辐射的散射、吸收,直接影响地气系统的辐射能量收支。其次,气溶胶可以作为凝结核影响云的辐射特性,增加云量,间接影响气候[1]。气溶胶作为一种主要的空气污染物,严重危害人类健康 [2-5]。与此同时,它影响卫星成像的质量,使图像模糊,从而影响地表信息的精确提取,是遥感数据大气校正的重要输入参数之一。因此,气溶胶信息的准确获取在气候变化、环境污染监测和卫星数据利用中具有重要作用。其中,气溶胶光学厚度(Aerosol Optical Depth, AOD)是卫星遥感监测的重要参数之一,是监测气溶胶的重要指标。目前,许多卫星在气溶胶影响环境和气候的研究中发挥了重要作用,例如Moderate Resolution Imaging Spectroradiometer (MODIS)、the Visible Infrared Imaging Radiometer Suite (VIIRS)和the Multiangle Imaging Spectro Radiometer (MISR)等。相较于这些卫星,Landsat8具有更高的空间分辨率,能够提供更多的局部细节,有利于区域性气溶胶研究。
现有多种基于物理模型、结构函数和机器学习等方法的气溶胶光学厚度反演算法。其中应用较为广泛的是基于物理模型方法的气溶胶光学厚度反演算法,其利用辐射传输模型模拟地表、大气和大气顶层的关系,从而构建查找表并从中检索气溶胶光学厚度,例如浓密植被算法[6]和深蓝算法[7]。气溶胶算法的关键就是如何精准确定地表反射率,浓密植被法通过统计学方法,发现红、蓝波段的地表反射率较低且与中红外波段的地表反射率存在明显的线性关系,进而确定了红、蓝波段的地表反射率;深蓝算法利用SeaWiFS地表反射率产品,通过最小值合成法,构建了地表反射率库以提供红蓝波段地表反射率,进而求解气溶胶光学厚度。这些算法被成功地应用于MODIS传感器中并用于生产全球气溶胶产品(MOD04)。MOD04与地基观测结果具有较高的相关性且适合大尺度应用,但MOD04的空间分辨率仅为3 km和10 km,不利于小尺度分析。同时这些算法是建立在一些参数(例如地表反射率)假设基础上的,这些假设会在某些区域产生一定的误差,从而影响气溶胶光学厚度反演的精度。例如,浓密植被法确定的地表反射率的统计关系只在暗黑像元处稳定,在城市高亮地表处有较大偏差。
研究者针对地表反射率确定问题进行了诸多改进,提出了多个高精度的气溶胶光学厚度反演算法,例如Multi Angle Implementation of Atmospheric Correction(MAIAC)算法[8]和Simplified Aerosol Retrieval Algorithm (SARA)[9]。MAIAC算法通过分析MODIS的长时间序列数据来确定地表反射率,此算法在反演AOD的同时也能反演细粒子分布,最终生成具有1 km分辨率的气溶胶产品,有效地提高了AOD的反演精度与分辨率。 SARA基于 MODIS地表反射率产品(MOD09GA),提供了地表反射率,通过提取分析研究区域对应的AERONET站点参数,改进了区域气溶胶模型。此算法成功应用于京津冀以及香港地区。但是,在SARA中,下垫面被假设为朗伯表面且气溶胶参数在同一区域中被假设为短时间保持不变,这使得SARA在不同地区的应用会出现偏差。
为了避免不同的参数假设,Li等[10]直接使用站点和卫星数据,利用自回归综合移动平均(Autoregressive Integrated Moving Average, ARIMA)法,进行了长期气溶胶趋势分析并预测了气溶胶光学厚度[11]。结果表明,利用时间序列数据进行建模并预测气溶胶光学厚度是可行的,该方法在区域范围内表现出良好的精度。此外,Radosavljevic等[12-13]利用神经网络直接学习卫星观测气溶胶光学厚度与地面观测气溶胶光学厚度之间的关系。相较于物理模型,深度学习方法能够获取更高的精度。Huttunen等[14]对比分析了查找表方法、非线性回归方法和四种机器学习方法的优劣,结果表明,大多数机器学习方法能够得到更高精度的气溶胶光学厚度反演结果。
卫星表观反射率与气溶胶光学厚度之间的关系是一个复杂的非线性物理问题,但机器学习是一个强大的学习工具,可以基于样本数据进行学习并实现对复杂方程式的精确求解,已被用于多类时间序列预测中。因此,本文利用深度学习算法,直接从Landsat8 陆地成像仪(Operational Land Imager, OLI)影像中估算气溶胶光学厚度。具体来说,本文利用概率学习模型即深度置信网络(Deep Belief Network, DBN) [15],建立了地表气溶胶光学厚度与卫星大气顶部(Top of Atmosphere, TOA)反射率、观测角度之间的关系。为了使样本更具代表性并且能够充分学习非线性关系,本文选取分布在全球且具有长时间序列的AERONET站点数据,根据时间和空间信息,将其与Landsat8表观反射率数据匹配,从而形成模型的训练样本。在样本训练中,Landsat8提供TOA反射率以表征地表大气耦合的影响,利用归一化差分植被指数(Normalized Difference Vegetation Index, NDVI)来区分不同地表类型,利用观测几何数据表征各种成像条件。本文利用DBN对卫星数据和地表测量的AOD之间的关系进行训练建模,从而检索特定区域的气溶胶光学厚度。对2019年所有站点数据进行了预测并验证了它们的精度,对研究区域的小尺度气溶胶光学厚度进行了更精细的验证。
2 数据及研究区域
2.1 研究区域
本文选取美国西北部的爱达荷州以及英格兰东南部的大伦敦作为研究区域。爱达荷州的地势起伏较大,森林茂密,植被丰富,森林覆盖率达到了33%,并且水源丰富,其境内的AERONET-DRAGON网络可以对AOD估算质量进行验证。大伦敦地区的下垫面以城市人造地表为主,植被覆盖率约为13%;其地势较低,以平原为主,地形平坦;同时,伦敦受北大西洋暖流和西风的影响,气候湿润且多雨雾。
2.2 地面气溶胶光学厚度的测量
AERONET(Aerosol Robotic Network )是一个全球地基气溶胶观测网络,可提供稳定且易于获取的气溶胶光学特性数据库。它具有较高的时间分辨率,每隔 15 min便提供一次观测结果,并在全球范围内设置有400个以上有效观测站点,常被用于气溶胶研究和卫星遥感气溶胶反演验证等。AERONET提供三个级别的AOD观测数据:Level 1.0(未筛选过的产品),Level 1.5(云剔除产品)以及 Level 2.0(云剔除和质量控制产品)。AERONET新的V3版本数据使用自动化的近实时质量控制算法,AOD计算结果更加精确(不确定度为0.01~0.02),并利用自动云筛选代替了先前的手动质量控制, Level 2.0数据更加稳定[16]。本文使用Level 2.0站点数据进行训练样本收集,在全球范围内挑选455个AERONET站点以尽可能涵盖更多的气溶胶类型。本文需要收集波长550 nm处的AOD以组成训练样本并验证模型估算精度,但AERONET站点缺少波长550 nm的数据,因此使用Ångström指数[17]的方法进行波段转化,具体公式为
式中:τa(λ)为波长为λ的气溶胶光学厚度;a、β分别为Ångström指数和大气浑浊度系数。分别使用受水汽影响较小且与550 nm 最接近的 500 nm 和 675 nm作为波长 λ1、λ2,求解两系数值a和β。
2.3 卫星数据
Landsat8搭载有先进的Operational Land Imager (OLI)与Thermal Infrared Sensor (TIRS)传感器,其中OLI的空间分辨率为30 m,卫星每16 d就可以实现一次全球覆盖,多光谱波段覆盖了可见光到短波红外波段,波段设置如
表 1. Landsat8 OLI数据
Table 1. Landsat8 OLI data
|
OLI不仅具有较高的空间分辨率,而且波段覆盖广,其在气溶胶监测方面得到广泛应用。Google Earth Engine (GEE)是一个基于云的行星规模地理空间分析平台,在云端就可以实现全球尺度遥感数据的处理。本文利用GEE,在全球范围内提取了2013—2019年的Landsat8表观反射率数据。Landsat8数据含有B1~B7波段的TOA反射率数据、波长从可见光到近红外的辐射光谱、相应的观测几何数据(太阳高度角和相对方位角)、NDVI数据和质量控制波段的数据。NDVI计算公式 [18]为
式中:ND为NDVI;B4、B5分别为Landsat8 卫星B4、B5波段的反射率。计算后的NDVI按月合成月数据集,将其代入AOD估算模型以表征不同的土地覆盖类型。另外,在质量评估波段(Quality Assessment Band, QA),通过受仪器或云层影响的信号标示像元,能够快速识别污染像元。本文基于Landsat8卫星2013—2018年的数据进行样本收集,并利用QA波段进行云像元的剔除,其中2018年部分数据用于生成独立验证集,2019年的数据用于AOD估算。
2.4 数据时空匹配
为了进行时间匹配,本文选择卫星过境前后30 min内的AERONET站点测量值与卫星数据进行匹配。若区间内出现多组站点的测量值,则本文采取线性插值得到当前时刻的AOD:
式中:τ为成像时刻的AOD;ta、tb为与成像时刻最为接近的两个时间点;τa、τb分别为对应时刻的AOD;t为Landsat8卫星成像时刻。在空间匹配方面,本文利用AERONET站点的经纬度信息,通过影像投影进行逐像元匹配。本文共收集20003条训练样本。
3 DBN的气溶胶光学厚度估算
Landsat8卫星B7波段的波长较长且受大气影响较小,B2和B4波段信息可以表征蓝、红通道上的地表与大气耦合反射率,因此地表反射率可以通过B2、B4、B7三个波段之间的统计经验关系进行分离[6]。Landsat8卫星新增加的B1波段相较于蓝波段对气溶胶更加敏感,更易受到气溶胶的影响。因此,可以从这些输入变量中得到大气对卫星TOA反射率信号的影响程度;剩余的波段能够增加更多的信息,有利于模型学习其中的关系。但是,AOD、卫星TOA反射率、观测角度以及NDVI之间的关系非常复杂。深度学习在拟合非线性和复杂关系方面具有很大潜力,可以更好地表征这种关系。
本文利用DBN拟合参数之间的关系。DBN是一种概率生成模型,它通过采用逐层训练的方式,解决了深层次神经网络的优化问题,为整个网络赋予了较好的初始权值,使得网络只要经过微调就可以实现最优解。其中,逐层训练的关键是深度波尔兹曼机(Restricted Boltzmann Machine, RBM),RBM包含了一个用于输入训练数据的显层 (v),还有一层作为特征检测器的隐含层 (h),层内无连接,层与层之间全连接,这大大减少了计算量和计算难度。RBM 的训练过程实际上就是求解最大概率分布,由于这个分布的决定性因素在于权值W ,因此训练 RBM 的目标就是寻找最佳的权值。在DBN中,前一层的RBM隐含层连接着后一层的可见层,多层堆叠的RBM组合后通过连接一层后向传播(Back-Propagation,BP)层[19]来进行微调输出,后向传播层不仅可以应用于预测也可以应用于分类,本文利用其估算AOD。
式中:N为神经元个数;n为特征个数;μ为输出变量个数。为了充分学习参数之间的关系并防止过拟合,本文使用了两层RBM层,隐含层的神经元个数设置为20。利用深度置信神经网络估算AOD主要分为三大步骤,流程图如
1) 将卫星TOA反射率、观测角度和NDVI输入DBN模型,使用RBM对收集的样本进行逐层的非监督预训练,非监督预训练能够提取出必要特征并将其传递给下一层RBM。
2)利用步骤1)的结果生成初始权重并计算得到 AOD;随后计算模型预估的AOD与AERONET实际站点数据间的均方根误差(MSE),将误差返回BP层并进行反复微调,直到模型达到满意的性能。
3) 模型在使用前需要进行交叉验证,精度满足后利用其估算没有地面监测站地区的AOD值,从而进一步得到Landsat8卫星的气溶胶光学厚度分布。
本文采用十折交叉验证 (10-fold cross-validation) 技术 [21]来检验模型的拟合和估算能力,采用基于样本的交叉验证[22]创建独立的验证数据集。在基于训练样本的交叉验证中,本文将训练数据分为10个子集,其中9个用来训练,剩余的1个用来验证;并且在2018年的站点数据中,挑选与本文研究区域具有类似气溶胶特性的站点数据以形成独立的验证数据集,进而进行交叉验证。本文使用相关系数、均方根误差以及平均绝对误差来评价验证精度。
4 实验结果与分析
4.1 AERONET站点验证
为了初步验证AOD模型的估算能力,首先利用2019年收集的全部数据进行550 nm处的AOD估算并利用AERONET站点进行精度验证。本文在500 m分辨率的尺度下验证AOD的估算精度,由于估算的是30 m分辨率下的AOD,需要进行像元合成。为了尽可能减少云对验证结果的影响,本文采用如下步骤进行像元合成:1)以站点为中心,取大小为17 pixel ×17 pixel的窗口;2)将所有像元按照B2波段(蓝波段)的数值从小到大进行排序;3)删除蓝波段反射率最高的像元(占比30%)和蓝波段最暗的像元(占比20%);4)剩余像元按照QA波段信息进行云判断,非云的像元取平均即为最终结果。本文使用卫星过境前后30 min的站点数据进行验证,并利用皮尔逊积矩相关系数(R)、均方根误差(RM)、平均绝对误差(MA)和MODIS期望误差(EE)来评价精度。计算公式分别为
式中:Aprediction为模型估计结果;AAERONET为站点观测结果;Cov(·,·)为协方差;Var(·)为方差。本文在全球范围内共匹配到233个站点的数据,共计1962条有效观测数据。
图 3. AOD预测值与AERONET站点观测值的密度散点图
Fig. 3. Scatter density map of AOD predication value and AERONET site observation value
为了评价算法的区域性估算精度,本文获取了实验区域内2019年共10个站点的数据并进行区域精度评价,其中大伦敦地区2个,爱达荷州8个。
表 2. 误差统计
Table 2. Error statistics
|
图 4. 模型预测结果与不同站点实测值的对比
Fig. 4. Comparison between model predication results and measurement results at different sites
4.2 气溶胶的空间分布
在所选研究区域内评估气溶胶的空间分布和变化情况。为了估算更高精度的气溶胶光学厚度,本文首先利用Sun等[23]提出的基于GlobeLand30数据库的云检测算法进行云水像元的筛除,同时使用归一化冰雪指数(NDSI,NS)进一步去除雪体。NDSI的计算公式为
式中:B3为绿波段反射率;B6为短波近红外波段的反射率。
5 结论
提出了一种基于DBN的Landsat8-OLI气溶胶反演算法。使用GEE在全球范围内收集训练样本,将2013—2018年的AERONET站点数据与Landsat8表观反射率数据进行匹配以形成训练样本,并利用2019年的影像数据进行全球站点的预测反演。结果表明,在全球范围内,所提算法可以预测出具有较高精度的AOD,均方根误差和平均绝对误差分别为0.11和0.07,约有69.5%的数据落在EE误差线内。在小尺度范围内进行了AOD的验证,结果表明,AOD预测结果和站点测量结果之间有着较高的相关系数(0.79),预测结果具有较低的均方根误差和平均绝对误差,分别为0.045和0.034。验证结果证明了所提算法具有较高的精度且可靠性好。最后对气溶胶的空间分布进行了分析,所提算法能够在浓密植被区域提供稳定、连续的气溶胶空间分布结果,并能够在城市高亮地表精准地识别薄雾。与传统方法相比,所提算法仅利用单相卫星遥感数据即可实现气溶胶光学厚度反演,无需复杂的物理建模,简化了反演过程。所提算法仍存在一些不足,由于Landsat8卫星的回访周期较长,部分地区站点数据的不充足会造成训练样本分布不均匀,东南亚及非洲沙漠等地区的气溶胶评估结果变得不可靠。这些地区气溶胶变化剧烈,在这些区域反演出具有高分辨率和高精度的气溶胶光学厚度是未来着重的研究方向。
[2] Gauderman W J, Avol E, Gilliland F, et al. The effect of air pollution on lung development from 10 to 18 years of age[J]. The New England Journal of Medicine, 2004, 351(11): 1057-1067.
[3] Künzli N, Kaiser R, Medina S, et al. Public-health impact of outdoor and traffic-related air pollution: a European assessment[J]. The Lancet, 2000, 356(9232): 795-801.
[6] Kaufman Y J, Tanré D, Remer L A, et al. Operational remote sensing of tropospheric aerosol over land from EOS moderate resolution imaging spectroradiometer[J]. Journal of Geophysical Research: Atmospheres, 1997, 102(D14): 17051-17067.
[8] Lyapustin A, Wang Y, Laszlo I, et al. Multiangle implementation of atmospheric correction (MAIAC): 2. Aerosol algorithm[J]. Journal of Geophysical Research: Atmospheres, 2011, 116(D3): D03211.
[9] Bilal M, Nichol J E, Bleiweiss M P, et al. A simplified high resolution MODIS aerosol retrieval algorithm (SARA) for use over mixed surfaces[J]. Remote Sensing of Environment, 2013, 136: 135-145.
[12] RadosavljevicV, VuceticS, Obradovic Z. Aerosol optical depth retrieval by neural networks ensemble with adaptive cost function[EB/OL].( 2017-09-28)[2020-05-10]. https:∥pdfs.semanticscholar.org/6e7a/00df19d6e2b3698073d1f96a12f16f3c5a89.pdf.
[15] Hinton G E, Osindero S, Teh Y W. A fast learning algorithm for deep belief nets[J]. Neural Computation, 2006, 18(7): 1527-1554.
[16] Giles D M, Sinyuk A, Sorokin M G, et al. Advancements in the Aerosol Robotic Network (AERONET) Version 3 database : automated near-real-time quality control algorithm with improved cloud screening for Sun photometer aerosol optical depth (AOD) measurements[J]. Atmospheric Measurement Techniques, 2019, 12(1): 169-209.
[17] Ångström A. The parameters of atmospheric turbidity[J]. Tellus, 1964, 16(1): 64-75.
[18] Chen J M, Cihlar J. Retrieving leaf area index of boreal conifer forests using Landsat TM images[J]. Remote Sensing of Environment, 1996, 55(2): 153-162.
[19] Rumelhart D E, Hinton G E, Williams R J. Learning representations by back-propagating errors[J]. Nature, 1986, 323(6088): 533-536.
[20] Fletcher D, Goss E. Forecasting with neural networks[J]. Information & Management, 1993, 24(3): 159-167.
[21] Rodriguez J D, Perez A, Lozano J A. Sensitivity analysis of k-fold cross validation in prediction error estimation[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2010, 32(3): 569-575.
Article Outline
梁天辰, 孙林, 王永吉. 基于深度学习反演区域气溶胶光学厚度[J]. 光学学报, 2021, 41(4): 0401002. Tianchen Liang, Lin Sun, Yongji Wang. Retrieval of Regional Aerosol Optical Depth Using Deep Learning[J]. Acta Optica Sinica, 2021, 41(4): 0401002.
PDF全文
