基于地表反射率关系库的Sentinel-3A AOD反演——以中原城市群为例

朱 世 鑫，潘 竟 虎

(西北师范大学地理与环境科学学院，甘肃 兰州 730070)

0 引言

大气气溶胶是悬浮在空气中的分散颗粒，包括初级气溶胶(以颗粒的形式直接排放到大气中)和次级气溶胶(由大气中的主要污染物转化而来)，直径为0.001～100 μm[1]。大气气溶胶的浓度变化直接影响人体健康与空气质量，并通过辐射强迫的直接或间接效应影响地气收支平衡和气候变化[2，3]，因此，对气溶胶的空间分布和变化进行监测非常重要。气溶胶光学厚度(Aerosol Optical Depth,AOD)是描述气溶胶特性的重要参数之一，可表示大气的浑浊度[4]。传统气溶胶监测多为地基监测，然而由于站点布设有限，难以获取大范围的气溶胶空间分布信息。近年来，利用具有高时效性和大尺度观测能力的卫星遥感技术可获取空间分布连续的气溶胶数据，在宏观环境和污染分布监测上具有较大潜力[5]。

AOD反演的关键是确定地表反射率信息，以便于实现地气解耦[6]，反演算法包括暗目标算法[7]、深蓝算法[8,9]、结构函数法[10]、偏振算法[11]、多角度算法[12]等。其中，结构函数法需找到“清洁日”影像作为基准且对几何校正要求较高，偏振算法仅能应用于反演细粒子气溶胶，多角度算法需要特定的传感器支持，3种算法的业务化应用较为困难。近年来暗目标算法不断得到改进与应用[13,14]，但其仍具有明显局限性，不仅需要短波红外波段的支持，且不适用于亮像元地区(城市、沙漠等区域)。深蓝算法可以反演出亮像元地区的AOD，但其精度低于暗目标算法，而且需要外部的地表反射率产品，目前常使用MOD09A1数据，空间分辨率为500 m；由于单波段的反演需要具体的地表反射率栅格值才能实现地气解耦，在反演过程中，如果表观反射率数据的空间分辨率小于500 m，通常会将其重采样为500 m[15-17]，最终得到500 m分辨率的反演结果，因此，深蓝算法反演结果受外部地表反射率数据空间分辨率限制。如何能在不损失原有数据空间分辨率的前提下，同时反演出暗、亮像元地区气溶胶光学厚度，是当前气溶胶光学厚度反演研究的一个重点。

Sentinel-3A卫星搭载的OLCI(Ocean Land Colour Instrument)传感器有21个通道，波段范围为400～1 020 nm，辐射分辨率和信噪比较高，AOD反演潜力很好，且其空间分辨率为300 m，高于常使用的MODIS传感器，但目前对于Sentinel-3A OLCI传感器的气溶胶反演研究还较少。Mei等[18]利用XBAER算法进行OLCI数据的AOD反演，但其构建地表反射率数据库的过程限制了算法应用。王峰等[19]采用红蓝通道比值算法，根据NASA发布的MODIS红蓝通道地表反射率关系式，经过光谱转换后获取适用于Sentinel-3A OLCI传感器的固定关系式，对中国台湾暗像元地区(浓密植被区)进行AOD反演的结果较好，但该反演结果不包括亮像元地区，且根据NASA得到的固定关系式没有考虑季节性和区域性差异，最终影响了反演精度。

本文基于现有理论，参考深蓝算法的地表反射率库获取方法，并根据红蓝通道比值算法的思想，以中原城市群作为研究区，拟合不同NDVI等级下的地表反射率关系式，建立逐月的地表反射率关系库，利用拟合得到的系数和截距反演AOD，以期消除外部地表反射率数据对反演结果空间分辨率的限制，获取同时包括暗、亮像元的高空间分辨率AOD反演结果，并利用SONET(Sun-Sky Radiometer Observation Network)地基观测数据检验AOD反演结果，用MODIS AOD产品与反演结果进行空间分布对比分析，以此检验反演结果的精度和可靠性。

1 研究区和数据

1.1 研究区概况

根据国家发改委2016年12月发布的《中原城市群发展规划》，中原城市群是中国东西向经济发展的过渡地带，涵盖河南省18个地级市以及河北、山西、安徽、山东等省份的部分地级市。该地区资源压力大，产业结构以第二产业为主，近年来发展速度较快，人们在日常生活和密集性生产中向城市大气中排放的气溶胶粒子日益增多，大气环境污染问题日渐严峻，尤其在污染物不易扩散的秋冬季更为严重[20，21]，将其作为研究区具有代表性。研究区内可用的SONET地基观测站点有焦作站和嵩山站(图1)。

图1 中原城市群范围和SONET 站点位置Fig.1 Scope and SONET site location of Central Plains Urban Agglomeration

1.2 数据与预处理

1.2.1 AOD反演所需数据

(1)MOD09A1数据。该数据是经大气校正后的地表反射率数据，选取8天内质量最佳的数据组合，包含MODIS传感器的Band1-Band7波段，用于构建地表反射率关系库。由于OLCI传感器与MODIS传感器之间存在光谱响应差异，因此需要利用实测光谱数据对MOD09A1数据进行修正，而光谱库在一定程度上可以解决缺少实测数据的问题[22]。将ENVI光谱库中的不同地物光谱重采样至1 nm间隔，与MODIS和OLCI的红、蓝波段的光谱响应函数分别进行卷积计算，得到不同地物在MODIS或OLCI的红、蓝波段下的真实反射率(式(1))；对不同地物真实反射率进行散点拟合，由式(2)模拟计算MODIS红(蓝)波段与OLCI红(蓝)波段的光谱转换统计模型，获取Sentinel-3A的地表反射率数据。本文采用最小值合成法处理每月的地表反射率数据，获得每月地表反射率最小值，以去除云、阴影等影响，并将结果作为当月的地表反射率真实值。

(1)

式中：Γ(λi)表示波长为λi时的传感器响应函数；R(λi)为λi处的地物光谱值；Δλ采用0.0001微分值。

ROLCI=ARMODIS+B

(2)

式中：ROLCI和RMODIS分别表示由式(1)计算得到的OLCI和MODIS的真实地表反射率；A、B为对应的拟合参数。

(2)Sentinel-3A数据。该数据来源于欧空局(https://scihub.copernicus.eu)，本文选取研究区晴空像元相对较多的L1B级影像数据进行AOD反演，时间范围为2019年3月至2020年2月；除各波段的辐射亮度数据外，辅助数据包括太阳(卫星)天顶角、太阳(卫星)方位角数据。利用SNAP软件对原始辐射亮度数据和几何数据进行预处理，得到投影转换后的表观反射率数据和角度数据，并使用编程语言进行数据裁剪及波段合成。

(3)DEM数据。该数据来源于地理空间数据云(http://www.gscloud.cn)的SRTMDEMUTM 90M数据集，用于对AOD反演结果进行高程校正。为便于与Sentinel-3A数据匹配，将90 m的DEM数据重采样为300 m。

1.2.2 验证数据

(1)SONET观测数据。SONET是由中国科学院联合中国多所高等院校和科研院所等单位在中国典型区域建立的观测网，相比于集中分布于京津冀、台湾等地区的AERONET(Aerosol Robotic Network)观测数据，SONET站点在中国境内分布更均匀，覆盖全国多个省域，目前共有23个站点。此外，SONET具有与AERONET相当的数据采集能力，且不确定性小于AERONET[23]。因此，本文选取研究区内的焦作站和嵩山站的Level 1.5级观测数据(http://www.sonet.ac.cn)，用于验证AOD反演结果的精度。由于反演结果是550 nm处的AOD，但站点观测数据没有550 nm处的AOD值，不便于直接对比，所以需要使用已有数据计算出550 nm处的AOD值作为验证数据。相较于Angstrom波长指数法，通过二次多项式法计算得到的AOD插值数据结果更可靠[24]。该方法使用3个已知通道的AOD拟合未知参数，然后计算出550 nm处的AOD值(式(3))，本文选取675 nm、500 nm、440 nm 3个通道进行计算。

ln(τλ)=a0+a1ln(λ)+a2(ln(λ))2

(3)

式中：τλ为λnm处AOD值；a0、a1、a2为对应的拟合参数。

(2)MOD04_3K数据。本文选用MOD04_3K AOD作为反演结果空间分布的检验数据，其与Sentinel-3A卫星过境时间相近，是NASA发布的Level 2级第6版气溶胶数据，空间分辨率3 km，在浓密植被区反演精度较高。

2 研究方法

2.1 基本原理

假设大气分子的分布与变化均一，且地表是均匀朗伯面，则卫星接收到表观反射率(ρTOA)的计算公式为：

ρTOA(us,uv,φ)=

(4)

式中：us、uv分别为太阳天顶角和卫星天顶角的余弦值；φ为相对方位角；Tg为总气体分子透过率；ρ0为大气路径程辐射项；T(us)、T(uv)分别为大气上行、下行透过率；ρs为地表反射率；S为大气下界的半球反射率。将式(4)变换可得地表反射率的求解方程：

ρs(us,uv,φ)=

(5)

由式(4)、式(5)可知，表观反射率中包含大气和地表的信息，因此，AOD反演的前提是实现地气解耦，去除地表信息，然后从大气信息中得到AOD的反演结果。具体流程(图2)为：利用辐射传输模型获取预设气溶胶类型条件下的参数查找表，结合表观反射率与角度文件，迭代输入不同AOD值得到相应的大气参数，并根据式(5)计算出对应的红蓝波段地表反射率；然后与已知的真实地表反射率关系式进行对比，当二者最接近时，当前AOD迭代值即为AOD反演值。

图2 研究流程Fig.2 Flow chart of the study

2.2 查找表构建

计算地表反射率所用的参数一般通过大气辐射传输模型获取，目前常使用6SV(Second Simulation of the Satellite Signal in the Solar Spectrum vector code)模型。但6SV 2.1模型中并没有OLCI传感器的光谱响应函数，需要利用FORTRAN语言对模型源码进行修改，添加OLCI传感器的光谱响应支持。在构建查找表时使用添加的波段，可减少由光谱响应函数不同造成的误差[25]。

若直接利用大气辐射传输模型计算，速度缓慢，为提高反演效率，通常采用查找表的方式进行AOD反演。本文在查找表的设置中，将太阳天顶角和卫星方位角的步长设置为6°，相对方位角的步长设置为12°，550 nm处AOD值的步长设置为0.05，高程设置为0，大气类型设置为中纬度冬季和中纬度夏季；此外，相较于其他气溶胶模式，大陆型气溶胶模式反演结果精度较高[26]，因此，本文将气溶胶模式设置为大陆型气溶胶。最终得到550 nm处AOD值和3种几何角度在不同组合条件下的大气参数查找表，参数包括：总气体分子透过率(Tg)、大气路径程辐射项(ρ0)、总大气透过率(T(us)×T(uv))、大气下界的半球反射率(S)，对于查找表中不存在的几何角度所对应参数，通过线性插值的方式获取。

2.3 云水掩膜

云像元的存在会对大气参数的精度造成影响，因此云掩膜不充分会导致AOD反演结果出现较大误差。常用的云掩膜方法需要借助热红外波段，从云层与地表存在温差的角度剔除云像元。由于OLCI传感器没有热红外波段，同时考虑到根据表观反射率及其空间变化特征设置合适的阈值，可有效地实现云检测[27]，本文采用如下针对OLCI传感器的快速云检测方法(图3)：在厚云区域，OA3通道中的表观反射率ρTOA3较高，利用此性质，以0.35作为阈值识别厚云像元[19]；根据薄云区域的非均一性，当ρTOA3的3×3像元标准差大于0.002时，或当ρTOA6的3×3像元标准差大于0.05时，判别为云像元；OA13通道中的表观反射率ρTOA13随云顶高度的增加而增加，而OA12通道中的表观反射率ρTOA12对云顶高度的依赖性非常弱[28]，因此二者在云像元的反射率比值大于非云像元的反射率比值，设定比值大于0.315时为云像元；此外，对初步的云识别结果边界进行5×5像元扩充处理，以去除云边缘和云阴影的不确定性。从图4可以看出，由于缺少热红外波段的支持，卷积云区域出现轻微漏检，但云掩膜算法整体运行稳定，能较准确地识别出云像元区域。

注:STD3×3表示3×3像元标准差，ρTOAi表示OAi通道中的表观反射率。

利用OA8和OA17通道的表观反射率计算NDVI，将NDVI小于0的区域判定为水体像元进行水掩膜。

2.4 地表反射率关系库构建

本文基于红蓝通道比值算法的思想，通过拟合不同NDVI下的红蓝波段地表反射率线性关系，建立地表反射率关系库，利用固定关系式，根据系数和截距进行反演，避免了计算的地表反射率与真实地表反射率之间的比较，从而实现高空间分辨率的气溶胶光学厚度反演。

首先，将OLCI表观反射率重采样到500 m，以便与地表反射率数据匹配。考虑到浓密植被区与非浓密植被区的差异，通过表观反射率计算NDVI，根据NDVI划分等级分别进行拟合：浓密植被区的红蓝波段线性关系较稳定，参考王峰等[19]的研究，将NDVI大于0.45的像元判断为浓密植被区；非浓密植被区由于植被稀疏程度不同，地表反射率关系较为复杂，因此，每0.05的NDVI设置一个级别，拟合不同NDVI等级下的红蓝波段地表反射率关系。最后，将拟合出的地表反射率关系式应用到原始的Sentinel-3A数据中，反演出300 m分辨率的AOD。拟合过程如下：假设短时间内地表反射率不发生变化或变化较小，以最小值合成法合成的地表反射率数据作为当月的真实地表反射率，使用每月的多景待反演影像计算NDVI，去除云、水像元后取NDVI均值，以保证有足够的像元参与地表反射率关系的拟合；然后读取不同NDVI下对应位置的红蓝波段地表反射率像元值，以NDVI等级分类并分别进行拟合，得到一阶线性方程式，建立地表反射率关系库。考虑到拟合关系的季节性差异，该地表反射率关系库是逐月建立的，表1为2020年2月红蓝波段地表反射率关系拟合结果。

图4 云识别效果与影像对比Fig.4 Comparison between cloud recognition effect and remote sensing image

表1 2020年2月地表反射率关系库Table 1 Relationship database of land surface reflectivity in February 2020

2.5 AOD反演

利用经过云、水掩膜处理后的300 m表观反射率数据计算出NDVI，根据地表反射率关系库确定该NDVI下使用的拟合关系式；然后以0.001作为初始AOD进行计算，步长β设置为0.05，根据预处理后的表观反射率和几何数据，代入式(5)迭代计算不同AOD条件下的红、蓝波段地表反射率ρ8和ρ4。由于表观反射率会随着AOD的增加而增加，且波长更短的波段对气溶胶增加更敏感[29]，因此，随着AOD的增加，(ρ4-b)/ρ8的结果会不断减小，当(ρ4-b)/ρ8小于系数a时停止迭代，并根据迭代前后的数据插值获取初步反演结果(图5)。

图5 AOD反演算法流程Fig.5 Flow chart of AOD retrieval algorithm

2.6 高程校正

单个查找表无法兼顾研究区内所有高程信息，建立多个高程条件下的查找表会影响计算效率，因此，本文以高程为0设置查找表，对AOD初步反演结果进行高程校正[30]。在目标高度为Z(单位为km，标准大气高度为8.5 km)时，大气分子的光学厚度可表示为：

τλ(z=Z)=τλ(z=0)exp(-Z/8.5)

(6)

2.7 精度验证

本文使用SONET站点数据对反演所得OLCI AOD进行精度验证，评价指标包括决定系数(R2)、均方根误差(RMSE)、期望误差(EE)[31]，计算公式如下：

(7)

EE=±(0.05+0.2×τ0)

(8)

式中：τ为反演AOD值；τ0为站点AOD值；n为验证样本数。

3 结果与验证

3.1 空间分布一致性检验

利用上述算法获得研究区内300 m空间分辨率的OLCI AOD数据，并与MOD04_3K AOD数据相比较，选择的影像类型包括：1)研究区AOD值区间较大，出现区域性雾霾天气(图6a和图6d)。OLCI反演结果与MOD04_3K AOD均可有效识别出研究区中部的焦作市、新乡市和郑州市与北部的邢台市、聊城市、邯郸市和安阳市的区域性雾霾天气，但OLCI AOD的空间分辨率更高，可清晰地识别出中部地区雾霾天气的细节信息；在城市群北部的邢台市，本文估算结果的有效像元比MODIS产品更多，MODIS产品在邢台市AOD高值区出现很多缺失像元，可能是由于重雾霾天气下的大气散射作用导致该区域反射率偏大，被误判为亮像元，因此MODIS在反演AOD时剔除了这些像元。2)研究区AOD值区间较小，空气质量较好(图6b和图6e)。OLCI反演结果与MOD04_3K AOD仍保持一致，AOD值较低的区域主要分布于西部秦岭等地区和西北部太行山区，在东部平原地区能明显看出城市区域AOD值较高,表明人类活动对气溶胶分布具有较显著影响。3)MOD04_3K AOD产品在晴空像元区域出现大量数据缺失(图6c和图6f)。北方地区冬季地表植被覆盖相对较少、地表反射率高，因此基于暗目标算法的MOD04_3K AOD产品在研究区内数据缺失较多，而本文反演的OLCI AOD仍可获取大范围的AOD分布，从而反映出研究区的雾霾分布特征。

图6 OLCI AOD与MOD04 AOD空间分布对比Fig.6 Comparison of spatial distribution between OLCI AOD and MOD04 AOD

综上可知，OLCI AOD与MOD04_3K AOD之间高值、低值区域相互对应，具有较好的空间分布一致性，但本文基于地表反射率关系库反演的AOD的空间分辨率更高，展现的细节信息更多，并且空间覆盖有效范围也更大。

3.2 SONET站点验证

首先，为使验证结果更合理，对过境时间前后的两个观测结果进行线性插值，以此作为验证数据；然后进行验证数据与反演结果的相关分析，并计算出精度验证评价指标，评估本文反演结果的精度(图7)。从图7可以看出，本文反演结果与SONET观测结果之间相关性较好，决定系数R2达到0.9099，但部分反演结果存在高估现象，RMSE低至0.1909，偏差整体较小，有67.44%的OLCI AOD在期望误差(EE)区间内，进一步说明本研究的反演结果整体反演精度和可靠性较高。

图7 OLCI AOD与SONET AOD精度验证Fig.7 Precision verification of OLCI AOD based on SONET AOD

根据图8a可以发现嵩山站点的反演结果较好，OLCI AOD与SONET的观测值和变化趋势较为接近，但当AOD接近0.1时，反演结果出现明显低估，这可能是由查找表的构建和插值方法所致[22]。焦作站点的反演结果(图8b)与观测值相比整体偏高，尤其是3月的反演结果高估现象较为明显,主要原因是该月份与焦作站点同等级NDVI的像元数量较少，影响了关系库的拟合精度，进而影响到反演结果，此外，地表反射率的数据质量也会影响反演结果。

图8 OLCI AOD与SONET AOD时序对比Fig.8 Comparison of time series between OLCI AOD and SONET AOD

3.3 误差分析

3.3.1 拟合误差影响 本文AOD反演是基于红蓝通道比值算法中固定系数与截距的思想。系数和截距是影响反演结果的重要因素，有必要进一步分析其差异性导致的误差，因此本文基于6SV模型定量模拟二者对AOD反演造成的误差。红蓝波段地表反射率比值多处于1.65～2.25之间[32]，因此定量模拟中假设比值为2.0，探究不同反射率地表下系数和截距所造成的误差，参数设定与分析结果如图9所示，其中，OZA，SZA，RAA，ρ8分别表示卫星天顶角、太阳天顶角、相对方位角和红光波段地表反射率，大气类型为中纬度夏季。

图9 拟合误差影响分析Fig.9 Analysis of the influence of fitting error

由图9a可知，0.01的截距误差大约会造成0.1的反演误差，不同地表反射率下的误差较接近。由图9b可知，随着地表反射率增大，系数误差引起的反演误差也增大，且偏差较大：当ρ8小于0.15时，0.1的系数误差约造成0.1的反演误差；当ρ8大于0.15时，0.1的系数误差至少造成0.2的反演误差。因此，随着地表反射率增大，拟合误差对反演结果的影响会更大。

3.3.2 误差来源分析 拟合误差是反演结果最重要的影响因素，此外，以下因素也在一定程度上影响反演结果。1)数据质量与数量：在使用地表反射率数据计算地表反射率关系式时，采用最小值合成法去除云、阴影的影响，但该方法不一定能完全剔除受影响像元，从而导致拟合误差，进而影响反演结果；此外，部分月份内某些NDVI等级所对应的像元数量不足，也会导致拟合质量不佳，从而造成反演误差。2)气溶胶类型：本文选择了普适性更强的大陆型气溶胶，未考虑气溶胶组分所引起的变化。随着AOD值的增加，气溶胶类型的影响会越来越重要[33]。3)地表反射率变化：本文假设地表反射率在短期内无变化或变化较小，然而在自然或人为影响下(如冰雪消融、农作物收割等),地表反射率会发生较大改变，造成反演误差[2]。4)太阳位置和传感器位置的差异：本文虽然已对MODIS地表反射率进行光谱转换以消除传感器不同产生的误差，但太阳位置、传感器位置等差异也会对反演结果产生一定影响。

4 结论

本文考虑区域性、下垫面不同及季节性变化，假设短期内地表反射率不变或变化较小，通过逐月拟合不同下垫面条件下的红蓝波段地表反射率关系式，建立地表反射率关系库，实现了高空间分辨率的AOD反演，得到以下结论：1)红蓝通道比值算法依据固定关系式进行反演，建立关系库可以避免与具体地表反射率栅格值的对比，不受外部地表反射率数据的限制，从而获取更高空间分辨率的反演结果。2)与SONET站点数据的验证结果显示，本文算法不仅能够获取更高空间分辨率的反演结果，而且反演质量也较高，决定系数R2达到0.9099，RMSE低至0.1909，67.44%的OLCI AOD落在预期误差区间(EE)内。3)系数和截距是影响反演结果的重要因素，0.1的系数误差会造成至少0.1的反演误差，0.01的截距误差大约会造成0.1的反演误差。

本文借助真实地表反射率数据建立了逐月的地表反射率关系库，由于小区域范围内并不能保证每一个NDVI等级内都有足够的像元进行关系拟合，从而导致部分月份反演误差相对较大，未来可以考虑引入历史同期数据，保证关系库建立所需的数据数量及质量达到要求；此外，还需进一步研究太阳位置、传感器位置等不同导致的反演结果误差。