作品介绍
1 应用背景
大气气溶胶是大气中重要的成分之一，是悬浮于大气中的固体和液体微粒与它们的气体载体共同组成的多相体系，其尺度大约在10-3到102 μm之间。大气气溶胶的特性对空气质量具有良好的指示作用，气溶胶的研究对空气质量的监测，污染物的溯源具有重大意义[1]。此外，气溶胶在地球大气辐射收支平衡和全球变化种扮演重要角色，对气溶胶光学特性的定量描述是其它研究开展的基础。
大气气溶胶光学厚度（Aerosol Optical Depth, AOD）是气溶胶的一个重要光学特征参量[2]，它不仅是大气浊度表征的质量因子，也是大气校正中最不可或缺的参数。当前对AOD的观测主要有两种方式：最精确的是气溶胶观测站组网，如AERONET（Aerosol Robotic Network）和SONET（Sun-Sky Radiometer Observation Network）[3]。它们虽然可以提供空间某点上具有时间连续性的气溶胶特性信息，但在空间上不连续，仅在有站点布设的小范围支持，无法满足环境监测实时、动态的要求；其次是基于遥感的反演观测，算法如暗目标法（Dark Target）[4]、深蓝算法（Deep Blue）[5]、结构函数法（Structure Function）[6]和多角度大气校正算法（MAIAC）[7]。相比于站点数据，卫星遥感反演AOD具有时空连续的优点，可直接作为气候模型的输入和影像精确校正，具有广阔的应用前景。
暗目标法适用于浓密植被的第地表反射率地表，对较高地表反射率的城市地区精度较差[8]。深蓝算法是一种用于反演高地表反射率地区AOD的算法。这两种算法均被应用于Moderate-Resolution Imaging Spectroradiometer（MODIS）生产AOD产品，但受限于数据源，它们的分辨率都较低（>10km），难以反应局部区域的气溶胶分布规律和污染源监测。目前Sentinel-2 A/B星协同观测可实现约5天的重访周期，提供13个波段的数据（10m、20m到60m）已成为监测高分辨率AOD的重要数据源。然而高空间分辨率影像的地表反射估算仍是一个挑战。因此，我们想到结构函数法（对比法），它无需对地表反射率估算的AOD反演算法[6,9]，是一种认为时间序列上（短期内）的两幅影像地表反射率不变，受到不同AOD影响的影像大气透过率不同（其中一个AOD/透过率已知）[6]，定义一种结构函数来表征透过率的变化程度[10]，进而反演未知AOD的方法。本质上它通过对比两景影像反演AOD，结构函数法只是用结构函数所表征的透过率来量化它们的变化程度。我们认为，结构函数法仍然存在不足：（1）机械地对一个窗口进行结构函数的运算，很容易受到异常点的影像，这引起随机误差；（2）仅仅使用透过率反演AOD极其容易受到气溶胶模型参数的影响，进而带来系统误差。
尽管如此，我们深受它的启发，并基于它的核心（对比法），提出一种消除上述误差的一种新算法。并将这种算法应用于一个典型区域——北京，与AERONET和广泛使用的MODIS的MAIAC AOD产品展开综合评估。
2 应用目标
（1）在结构函数法假设条件下，提出一种适用于高分辨率影像的新算法；
（2）采用Sentinel-2第一波段（443nm），将算法应用于北京地区的高分辨率AOD生产；
（3）与地面站点AERONET和MODIS的高精度的MAIAC AOD进行综合评估验证；
（4）尝试用本作品算法生产的高空间分辨率AOD监测低分辨率难以观测到的污染事件。
3 主要技术流程
图1 技术流程
为方便理解，介绍流程之前，有必要强调一些概念：
（1）本作品算法/结构函数法反演AOD需要两景影像，一景晴朗日（AOD已知）、一景污染日（AOD待求），晴朗日污染日地表反射率不变；
（2）反演的最小单元是窗口，窗口内大气特性稳定（AOD、透过率和程辐射不变），随着窗口的缩小，反演AOD分辨率随之提高；
（3）本作品算法气溶胶模型来自分析的AERONET站点（2017~2021）的参数。
本研究拟提出一种高效、稳定且参数不敏感的改进结构函数法，用于反演高分辨率遥感影像（<1km）的AOD，核心代码均用IDL实现，作品大致流程如图1所示：
本研究技术流程如下：
（1）首先准备三个数据集（Sentinel-2、AERONET和MODIS MAIAC）和一个辐射传输模拟工具6SV（Second Simulation of a Satellite Signal in the Solar Spectrum）；
（2）遍历匹配的晴朗（AOD<=0.1）、污染日影像（与晴朗日相差45天以内）逐窗口反演AOD；
（3）将反演AOD与MAIAC产品逐像元评估，与AERONET像元-点对比评估。
3.1 数据准备
3.1.1Sentinel-2数据
Sentinel-2是欧洲空间局哥白尼计划下的一个地球观测任务，该任务主要对地球表面进行观测以提供相关遥测服务，例如森林监测、土地覆盖变化侦测、天然灾害管理。它由两颗极地轨道卫星组成：Sentinel-2A（2015年6月23日发射）和Sentinel-2B（2017年3月7日发射），每颗卫星都携带一个多光谱仪器（MSI），在290公里宽的轨道上获取从可见光到短波红外的13个光谱段。每颗Sentinel-2卫星的重访时间为10天，两颗卫星合并后，重访时间缩短为5天。对外开放的Level-1C数据是经过几何精校正的表观反射率数据，包括13个波段，其中第一波段专为气溶胶探测设计[14]。Level-2A是经过官方工具Sen2Cor校正后的地表反射率影像，此外该工具处理过程中会对地物进行分类并生成场景分类图（Scene Classification Map, SCL）。SCL提供每个像元所属云、云阴影和水体等地物的索引，如表1所示。
表1场景分类图（Scene Classification , SCL）
从欧空局官网（https://scihub.copernicus.eu/dhus/#/home）可下载Level-1C数据和Level-2A数据。Level-1C数据包含反演所需的第一波段（与北京市重叠面积最大的瓦片ID为T50TMK，如图2所示）；Level-2A数据包含SCL波段[14]，是去云的重要参考。但因SCL的水体分类效果差，又考虑到水体变化较小，故采用Pekel等发布的全球水体数据作为水体掩膜[15]，而仅采用SCL提供的有关云、云阴影和雪分类索引（3, 8, 9, 10, 11）作为掩膜。
AERONET是以法国CIMEL公司生产的CE-318型太阳光度计为主要观测仪器的地基气溶胶遥感网络[16]，可提供全球超过500个站点的气溶胶光学厚度、Ångström指数、复折射指数、单次散射反照率、粒子谱分布等气溶胶特征参数，站点基本都分布在全球的主要区域，实现了每天采集数据的功能，将数据上传到网络中由统一的算法进行数据处理。本作品所用AOD时间分辨率为15min，粒子谱分布和复折射指数时间分辨率为3h。选取上述Sentinel-2瓦片包含的北京地区所有站点（共5个）名称分别为：Beijing、Beijing-CAMS、Beijing_PKU、Beijing_RADI和XiangHe，
如图2。
图2研究区与北京市相对位置关系和AERONET站点分布情况
AERONET数据
3.1.2 MAIAC
MODIS（Moderate-Resolution Imaging Spectroradiometer）的MAIAC提供1km分辨率的全球AOD产品，经过全球广泛验证，其精度高于暗目标和深蓝算法[17]。MAIAC AOD数据在Google Earth Engine（GEE）平台镶嵌后下载。
3.2 数据分析
3.2.1 气溶胶模型参数确定
在6SV中，气溶胶模型可由多个对数正态分布刻画，即[18]：
（1）
其中
表示粒子峰值体积浓度，表示粒子平均半径，表示标准差，表示粒子体积浓度随半径分布（后称粒子分布）。本作品认为m=2，即北京地区气溶胶的粒子分布是由粗、细单模态叠加形成的双模态粒子分布。我们用公式（1）拟合太阳光度计（AERONET）记录的离散粒子分布（2017~2021），然后按照不同AOD区间统计描述对数正态分布的三个参数（、和），最终分析得到这些参数（粗细模态三参数+复折射指数实部虚部共8个）与AOD的线性关系。最终结果如表2。
表2粒子体积分布参数和负折射指数
注：表中τ表示AOD，若参数与AOD有明显线性关系时，描述为AOD的函数，否则用均值描述。
3.2.2 晴朗日和污染日确定
污染日和晴朗日的确定步骤为：（1）在GEE平台对云量初筛，选择云占比<20%的所有影像（2017-2021）；（2）根据影像观测时间和AERONET确定研究区所有站点AOD均值，小于0.1为晴朗日，大于0.1为污染日；（3）目视检查云是否覆盖站点所在区域（可以反演但无法参与评估），进一步舍弃一部分影像。最终得到200景可用影像。
在实际反演时，并非所有污染日都可以参与反演，而是满足相隔天数低于一定阈值（45天）的影像对参与。这是考虑地表反射率短期不变的假设。
3.3 结构函数法误差来源
传感器接收的信息主来来自地表辐射和大气辐射两部分。假设大气水平均一、地表朗伯均一，传感器接收的表观反射率可表达为[19]：
（2）
其中，
– 太阳天顶角；
– 卫星高度角；
– 相对方位角；
– 大气光学厚度；
– 大气程辐射贡献反射率；
– 大气下行透过率；
– 大气上行透过率；
– 地表反射率；
– 大气半球反照率；
最简单的结构函数
定义为窗口内相距不远的两像元差的平方，然后将这种运算在窗口内逐像元执行即：
（3）
其中d为像元距离，i，j为窗口内像元索引。
结构函数法本质上是对比，它认为对于两景受到不同AOD影响的卫星影像，其中一景AOD已知（晴朗日），另一景（污染日）可通过与晴朗日对比确定，只是这种对比用结构函数的形式量化。具体来讲，结构函数法假设在一定窗口范围内大气特性稳定（AOD、透过率、程辐射不变），然后用污染晴朗日结构函数比值表征相应透过率的比值，进而在已知晴朗日AOD的情况下计算污染日AOD。从结构函数定义上，如公式（2），对窗口内所有像元（无选择地）计算得到的结构函数极易受到异常值的影像，异常值可能由于地表像元的非朗伯特性或是窗口内大气特性不再恒定。这容易产生随机误差。除此之外，因为结构函数法仅仅用透过率的关系来反演AOD，即：
（3）
其中，
和分别为晴朗日和污染日，这里的T为总透过率=，其中晴朗日AOD已知（)可求）。只利用透过率来连接反演的污染日AOD，极其容易受到气溶胶模型偏差的影响，误差会直接传播到反演的AOD上，从而带来系统误差。
3.4 本作品算法原理及详细步骤
针对上述误差来源，在结构函数法的启发下，提出一种新的算法（本作品算法）它克服结构函数法的两个误差来源，但仍基于对比法的思想。下面将详细介绍本作品算法原理和对上述误差来源的处理。
对于某一晴朗日和污染日窗口，当晴朗日表观反射率
、AOD和气溶胶模型已知（见3.2.1），晴朗日地表反射率的近似ρ可由6SV计算；当污染日气溶胶模型已知（见3.2.1）、表观反射率已知，可假设若干AOD得到若干污染日地表反射率近似值。我们认为最接近晴朗日地表反射率近似值所对应的AOD即为反演污染日AOD，如图1所示。即：
（4）
其中，F为损失；i,j为行列索引，k为预设的污染日AOD索引；m为窗口大小；
即最优AOD对应最小损失。这种考虑校正后代价函数最小的算法，不仅考虑透过率也将程辐射引入作为反演AOD的约束条件，可以一定程度上消除气溶胶模型带来的系统偏差；对于随机误差，在计算代价函数时，并非所有像元都参与计算，而是去除离群值后的像元，这就减弱了异常值对反演AOD的影响。
虽然运用迭代法可求得
，但与结构函数法相比效率很低，提出一种对上述介绍的迭代法更快速近似求解的插值方法。首先将小窗口所有像素值用一个统计量概括，本作品使用均值，在计算均值之前需要提出三倍中误差外的离群值。值得强调的是Hau等的算法实际上也经历这一步，Hau算法将20m分辨率的蓝光波段数据重采样至100m分辨率以降低配准误差。然后用6SV构建污染日条件（时间角度等）下AOD=[0.01, 0.25, 0.50, 1.00, 2.00]的查找表，经与晴朗日参数运算可得到公式（5）中slope、intercept和对应AOD的查找表，将晴朗日的统计量带入可得到模拟污染日统计量和对应AOD的查找表，最后再用quadratic插值（IDL中interpol函数，添加quadratic参数）可直接得到污染日AOD。
3.4.1 窗口大小确定
图3 基于北京地区Sentinel-2波段1反演和MODIS MAIAC AOD空间分布图
（前三行分别采用m=5、10和20的窗口反演AOD分布图，最后一行是MODIS MAIAC产品AOD分布图；不同列代表不同UTC时间）
基于上述算法，反演不同窗口（300m、600m和1200m，图前三行）的反演结果在空间上变化极小，说明本作品算法的稳定性（精度对窗口不敏感），可根据实际需要反演所需精度。但随着窗口减小，往往难以满足窗口内存在朗伯像元，结果会出现噪点。综合来看m=10，600m兼备高分辨率和较少噪点，因此在后续评估中将窗口设置为10。
3.5 数据评估
反演AOD分辨率为600m，将MODIS MAIAC（分辨率为1km）重采样至600m然后与反演AOD对比，由于Sentinel-2与MODIS过境时间不同这里的对比尚不能完全说明反演数据精度。因此，提取反演AOD部分栅格（站点所在）与AERONET前后半小时记录的AOD对比，这更能说明反演数据的真实精度。
3.5.1 与AERONET验证
将AERONET AOD插值至550nm并与本作品算法反演结果对比站点值（AERONET AOD 550nm）和站点所在栅格值（反演AOD 550nm）。分站点评估（a-e）和所有站点（f）参与评估结果如图所示。分站点评估相关系数>0.9，总体精度R=0.931，说明本作品算法反演北京地区AOD与AERONET一致性高，拟合线接近y=x说明气溶胶模型参数合理。
图4 分站点（a-e）和所有站点（f）与反演结果对比散点图
图5 与AERONET AOD的AOD偏差分布图
（灰点和误差线分别代表AOD-AEROENT AOD的50%、5%和95%分位点值）
从AOD偏差分布图（图）说明本作品反演AOD与AERONET AOD偏差在±0.2之间，低AOD值的较大负偏差可能与晴朗日的选取有关，本作品算法高AOD值表现稳定，有轻微高估倾向，这可能与大气稳定窗口减小有关。
3.5.2 与MAIAC验证
图6 所有参与反演栅格与MAIAC AOD对比散点密度图
MODIS MAIAC AOD产品是目前精度较高且广泛使用的AOD数据集，与它的验证对无AERONET站点的栅格精度具有一定指示作用。将所有参与反演的污染日的AOD与MAIAC AOD逐栅格对比，并去除明显不可信点，绘制得到散点密度图如图所示。R=0.907说明本作品算法与MAIAC具有极高的一致性，拟合线斜率<1，截距接近0（0.013），说明反演AOD倾向于低估高值AOD，而这一结论与AERONET验证所得结论相反，推测可能是MODIS与Sentinel-2在北京地区过境时间相差约6h所致。同样绘制AOD偏差分布图（图），可以发现对于低于0.8的AOD，偏差较小（±0.3之间）；当AOD>0.8，随着AOD增加，低估程度越来越大，这也与前述散点密度分析结果一直。
综上，与MAIAC对比弥补了AERONET仅仅可用于单点验证的不足，说明本作品算法所得AOD整体可信度较高，适合用于大规模生产AOD产品。
图7 与MAIAC AOD的AOD偏差分布图
4 关键技术
本作品有三个关键技术：其一是分析结构函数法的误差来源并在本作品算法中尽可能规避；其二是本作品算法的提出；其三是气溶胶模型选择。
4.1 结构函数法误差来源
有关结构函数法详细介绍和证明可参考周春艳等[9]。结构函数法认为在气溶胶的影响下，两窗口（晴朗污染）表观反射率的变化程度可以反应AOD，这一变化可由窗口内两距离不远的像元差的平方确定，并且为了减小配准误差将这种计算扩展至窗口多对像元，整个运算整理成函数就是结构函数[6]。在一系列假设下，污染日和晴朗日结构函数的比值等于对应总透过率（是AOD和气溶胶模型，角度等的函数）的比值，最后确定气溶胶模型后可由查找表确定最终AOD。结构函数法最初提出被用于沙漠地区反演，但针对城市地表，存在较多局部人为制造的气溶胶，难以保证地表朗伯，容易受到异常点的影像从而引起AOD的随机误差。
实际上，因为结构函数法仅仅使用透过率的关系反演气溶胶（在计算结构函数时将程辐射消除），所以对气溶胶模型比较敏感，在匹配性不强的气溶胶模型下往往会带来较大的系统误差。
4.2 本作品算法
假设两窗口AOD和气溶胶模型已知，认为大气校正后的地表反射率的估计值（由于地表复杂这里得到的只是地表反射率的近似）是接近的，这种接近程度用一种代价函数来刻画，而假设中的AOD稍有变化，便会破坏这种近似。
基于此，实际反演过程中，只有污染日AOD未知，可设置一系列AOD的离散值，带入代价函数，代价函数最小的情况下对应的AOD即为所求。在这个过程中，我们不仅使用了透过率，还考虑了程辐射的约束，这降低了算法对气溶胶模型的敏感程度可以消除系统误差。此外，在计算代价函数时，并非所有点都参与计算，而是去除一定的离群值后的，这可以消除异常点，减小随机误差。
4.3 气溶胶模型
在6SV[12]大气辐射传输模型中，气溶胶模型由包括粒子分布（认为研究区是双模态（Log-Norm）分布，每个分布包含三个参数）、复折射指数确定。从AERONET的lev15级2017~2021年间的数据中拟合出这些参数随AOD变化的一次函数系数[13]。最终得到研究区的气溶胶模型。
参考文献
[1]郭爽. 南京北郊PM_(2.5)中典型有机物及生物质燃烧示踪物的定量分析及来源解析[D]. 南京信息工程大学, 2021.
[2]贾臣, 孙林, 陈允芳, 吕喆鹏, 韦晶, 于会泳, 田信鹏. 查找表方法确定气溶胶类型[J]. 遥感学报, 2017, 21(03): 386–395.
[3]Li Z Q, Xu H, Li K T, Li D H, Xie Y S, Li L, Zhang Y, Gu X F, Zhao W, Tian Q J, Deng R R, Su X L, Huang B, Qiao Y L, Cui W Y, Hu Y, Gong C L, Wang Y Q, Wang X F, Wang J P, Du W B, Pan Z Q, Li Z Z, Bu D. Comprehensive Study of Optical, Physical, Chemical, and Radiative Properties of Total Columnar Atmospheric Aerosols over China: An Overview of Sun–Sky Radiometer Observation Network (SONET) Measurements[J]. Bulletin of the American Meteorological Society, 2018, 99(4): 739–755.
[4]Levy R C, Remer L A, Mattoo S, Vermote E F, Kaufman Y J. Second-generation operational algorithm: Retrieval of aerosol properties over land from inversion of Moderate Resolution Imaging Spectroradiometer spectral reflectance[J]. Journal of Geophysical Research: Atmospheres, 2007, 112(D13).
[5]Hsu N C, Tsay S-C, King M D, Herman J R. Aerosol properties over bright-reflecting source regions[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(3): 557–569.
[6]Tanre D, Devaux C, Herman M, Deschamps P Y. Estimation of Saharan aerosol optical thickness from blurring effects in Thematic Mapper data[J]. Journal of Geophysical Research, 1988, 931: 15955–15964.
[7]Lyapustin A, Wang Y, Laszlo I, Kahn R, Korkin S, Remer L, Levy R, Reid J S. 2. Aerosol algorithm[J]. Journal of Geophysical Research: Atmospheres, 2011, 116(D3).
[8]Gupta P, Levy R C, Mattoo S, Remer L A, Munchak L A. A surface reflectance scheme for retrieving aerosol optical depth over urban surfaces in MODIS dark target retrieval algorithm[J]. Atmospheric Measurement Techniques, 2016, 9(7): 3293–3308.
[9]周春艳. 环境一号卫星北京地区高分辨率大气气溶胶光学厚度反演算法研究[D]. 北京: 中国科学院遥感应用研究所, 2009.
[10]朱琳, 孙林, 杨磊, 徐菲菲, 徐青山. 结构函数法反演气溶胶光学厚度中像元的间隔设置[J]. 遥感学报, 2016, 20(04): 528–539.
[11]Vermote E F, Tanre D, Deuze J L, Herman M, Morcette J-J. Second Simulation of the Satellite Signal in the Solar Spectrum, 6S: an overview[J]. IEEE Transactions on Geoscience and Remote Sensing, 1997, 35(3): 675–686.
[12]Wei X, Chang N-B, Bai K, Gao W. Satellite remote sensing of aerosol optical depth: advances, challenges, and perspectives[J]. Critical Reviews in Environmental Science and Technology, 2020, 50(16): 1640–1725.
[13]Levy R C, Remer L A, Dubovik O. Global aerosol optical properties and application to Moderate Resolution Imaging Spectroradiometer aerosol retrieval over land[J]. Journal of Geophysical Research: Atmospheres, 2007, 112(D13).
[14]Louis J, Debaecker V, Pflug B, Main-Knorn M, Bieniarz J, Mueller-Wilm U, Cadau E, Gascon F. SENTINEL-2 SEN2COR: L2A Processor for Users[A]. L. Ouwehand. Proceedings Living Planet Symposium 2016[C]. Prague, Czech Republic: Spacebooks Online, 2016, SP-740: 1–8.
[15]Pekel J-F, Cottam A, Gorelick N, Belward A S. High-resolution mapping of global surface water and its long-term changes[J]. Nature, Nature Publishing Group, 2016, 540(7633): 418–422.
[16]Holben B N, Eck T F, Slutsker I, Tanré D, Buis J P, Setzer A, Vermote E, Reagan J A, Kaufman Y J, Nakajima T, Lavenu F, Jankowiak I, Smirnov A. AERONET—A Federated Instrument Network and Data Archive for Aerosol Characterization[J]. Remote Sensing of Environment, 1998, 66(1): 1–16.
[17]Lyapustin A, Wang Y, Levy R, Remer L, Hsu C, Reid J. A comparison between MODIS Dark Target, Deep Blue and MAIAC Aerosol Algorithms over Land[M]. AGU Fall Meeting Abstracts, 2010.
[18]Tanre D, Deuze J L, Herman M, Santer R, Vermote E. Second Simulation of the Satellite Signal in the Solar Spectrum, 6S: an overview[A]. 2002.
[19]Tanre D, Herman M, Deschamps P Y, Leffe A de. Atmospheric modeling for space measurements of ground reflectances, including bidirectional properties[J]. Applied Optics, 1979, 18(21): 3587–3594.
。
