杨磊库，胡秀清，王涵，何兴伟，刘培，徐娜，杨忠东，张鹏

1.河南理工大学测绘与国土信息工程学院,焦作454003;

2.国家卫星气象中心,北京100081

1 引言

气溶胶是地球大气系统的重要组成部分（Lenoble 等，2013），通过对太阳辐射的吸收和散射作用引起的直接辐射强迫效应、以及与云的相互作用引起的间接辐射强迫作用，影响地球辐射收支及气候，在辐射强迫估计中是不确定最大的变量之一（IPCC，2013）。气溶胶严重时还会引起环境污染，造成空气质量下降，从而威胁人类健康，近年来得到广泛关注（Li等，2017）。

卫星遥感从全球尺度提供气溶胶的持续观测信息，可以提高对气溶胶特性的认识能力（Kokhanovsky，2013；Guo 等，2016），被国际上认可为气溶胶监测最有潜力的手段之一（Lenoble等，2013）。近年来，用于气溶胶观测的卫星平台在不断发展，具有全球观测的极轨卫星，如美国国家航空航天局NASA（National Aeronautics and Space Administration）的TERRA、AQUA、美国国家海洋和大气管理局NOAA（National Oceanic and Atmospheric Administration）的JPSS（Joint Polar Satellite System）、欧空局ESA（European Space Agency）的Envisat、Sentinel以及中国的风云三号（FY-3）卫星；具有区域高时间频次观测的新一代静止卫星，如日本的Himawari-8；美国的Goes-R；中国的风云四号（FY-4）（卢乃锰和谷松岩，2016；Zhang等，2019）。搭载的传感器不断改进，从MODIS（MODerateresolution Imaging Spectro-radiometer）、VIIRS（Visible infrared Imaging Radiometer）、MERSI（Medium Resolution Spectral Imager）等多光谱观测，发展到法国的POLDER（POLarization and Directionality of the Earth’s Reflectances）和中国高分五号搭载的DPC（Directional Polarimetric Camera）传感器基于多角度偏振观测对气溶胶特性的反演能力（Li 等，2018；Dubovik 等，2019）。也因此涌现了大量的气溶胶产品和相应的算法，如NASA MODIS、欧空局Aerosol_cci 计划的气溶胶产品（Kokhanovsky 和De Leeuw，2009；Popp等，2016）。

与海洋相比，陆地表面反射率高且气溶胶类型复杂，其上空气溶胶反演仍是当前研究的热点和难点，主要取决于地—气解耦的问题，归结为两点：地表的准确估计和气溶胶类型的合理假设（Kokhanovsky 和De Leeuw，2009）。多角度卫星观测通过观测量增加可以在一定程度上消除地表的影响和帮助气溶胶类型的合理确定，从而提高陆地气溶胶的反演精度，如针对MISR（Multi-Angle Imaging SpectroRadiometer）传感器算法（Garay等，2020）；Dubovik 等（2011）针对POLDER 传感器的多角度偏振观测开发的GRASP（Generalized Retrieval of Aerosol and Surface Properties）算法不仅气溶胶光学厚度AOD（Aerosol Optical Depth）反演精度较高，且Ångström 指数AE（Ångström Exponent）和细粒子比FMF（Fine-Model Fraction）、单散射反射率SSA（Single-scattering albedo）等气溶胶特性参数也具有较好的反演验证精度（Tan 等，2019；Wei 等，2020）；尽管多角度传感器对气溶胶反演具有很大优势，单角度观测的传感器仍是目前气溶胶观测的一个重要手段，如MODIS 气溶胶产品验证精度较高、稳定可靠且每天全球覆盖1 次（TERRA 和AQUA 双星组合可以实现2 次覆盖），被广泛使用。接替MODIS 传感器的VIIRS 传感器搭载于Suomi NPP 和NOAA-20（也称作JPSS-1）上陆续发射升空。

单角度观测传感器由于观测量有限，陆地产品主要为气溶胶光学厚度AOD，一般需要借助于先验知识对反演模型中的地表反射率参数进行估计。典型的算法如MODIS 的暗目标算法DT（Dark Target）和深蓝算法DB（Deep Blue）（Hsu等，2013；Levy 等，2013）。DT 算法一般假设波长更长的2.12 μm 波段受气溶胶的影响较小甚至可以忽略，然后通过统计的地表波段比值关系作为先验知识来估算红光0.65 μm和蓝光0.47 μm波段的反射率，进而实现气溶胶的反演。DT 算法首次应用在MODIS 上实现业务化产品（Kaufman 等，1997），这也是卫星遥感对气溶胶探测的一个里程碑。MODIS 气溶胶产品及DT 算法已经历几代的发展，从最初的C4 版（Collection 4）算法产品（Kaufman等，1997），发展到第二代的C5 版（Collection 5）算法产品，随后的C6 版（Collection 6）（Levy 等，2013），以及最新发布C6.1 版算法产品。基于DT算法的MODIS 气溶胶产品具有较高的全球验证精度（Levy 等，2010；Munchak，2014；Sawyer 等，2020），并被移植到许多类似的传感器如VIIRS、新一代静止卫星Himawari-8（H8）上搭载的AHI（Advanced Himawari Imager）传感器（Jackson 等，2013；Levy 等，2015；苏城林等，2015；葛邦宇等，2018；Sawyer 等，2020），这些传感器均具有对气溶胶信号敏感的蓝光0.47 μm 和红光0.65 μm波段，以及与MODIS 2.12 μm相似的短波红外波段用于地表估计。DB 算法主要针对沙漠、干旱/半干旱的亮地表而设计，利用深蓝波段（或近紫外）地表反射率比较低、而气溶胶信号贡献较大的特性，通过挑选晴朗天像元并剔除瑞利散射，经多天合成的地表反射率库作为先验知识估算地表反射率，该算法在C5.1 版中首次被纳入MODIS 业务气溶胶产品。经过多年的数据积累，Hsu 等（2013）在C6 版DB 算法中合成了随散射角变化及考虑BRDF 特性新的地表反射率先验知识库。DB算法气溶胶产品总体具有较好的验证精度（Sayer等，2013），但对多年数据积累并建立地表反射率库依赖性强，如在一些地区依然存在严重低估的现象（Tao 等，2017），为此C6.1 版算法重新修订了这些地区地表反射率库。DB 算法在亮地表上空的气溶胶产品目前已得到基本认可，算法移植到VIIRS 上（Hsu 等，2019），并取得了较好的验证效果（Sayer 等，2019）。其通过多天合成地表反射率库作为先验知识进行地表估计的方法也被广泛采用，尤其适合缺少类似2.12 μm短波红外波段的传感器，如针对欧空局MERIS（MEdium Resolution Imaging Spectrometer）和OLCI（Ocean and Land Colour Instrument）传感器的XBAER（Xtensible Bremen AErosol Retrieval）算法（Mei等，2018），韩国GOCI（Geostationary Ocean Color Imager）静止卫星传感器的YAER（Yonsei Aerosol Retrieval）算法（Choi 等，2018）。另外需要提及的是，MODIS 地表验证小组的Lyapustin 等（2011）利用上/下午卫星16 d 的MODIS 单角度观测数据构建了基于多时相多角度的气溶胶反演MAIAC（Multi-Angle Implementation of Atmospheric Correction）算法，以联合求解地表和大气信息，目前最新C6 版产品也具有较好的全球验证精度（Lyapustin等，2018）。

风云三号（FY-3）系列卫星也属于极轨卫星，其上搭载的MERSI 系列传感器与MODIS 在波段设计上有许多相似之处，也可用于气溶胶反演（Yang 等，2019），将对长时间序列地球气溶胶卫星观测提供重要补充（Levy 等，2015；Sawyer 等，2020）。然而，目前还缺少精度可靠、稳定且全球适用的业务化AOD产品。最新发射（2017-11）的风云三号第4 颗星FY-3D，其上搭载了升级版MERSI 传感器（MERSI-II），在各方面性能上有明显的提高（胡秀清等，2018），挖掘其气溶胶反演定量能力非常必要和迫切。由于MERSI-II/FY-3D目前积累的数据较少，无法构建可靠的地表先验知识利用DB 算法进行反演，本研究初步借鉴DT算法进行MERSI-II/FY-3D 气溶胶反演能力测试。然而，MERSI 的波段设置在光谱响应和空间分辨率方面与MODIS 又存在一些差异。为此，本文在地表波段关系先验知识建立、像元滤除与聚合等方面进行了针对性的改进和优化，从而构建了适合MERSI-II的气溶胶反演算法。

2 数据介绍及辐射量转换

2.1 MERSI-II/FY-3D数据介绍

FY-3 是中国建设的第二代极轨系列气象卫星，具有全球观测的能力，搭载了包括紫外、可见光至红外及微波探测的多台遥感仪器（董超华等，2010；张鹏等，2012；杨忠东等，2013）。2008年5月首颗FY-3A 星上天，FY-3B、FY-3C 星分别于2010年11月和2013年3月陆续升空（杨忠东等，2017），最近又发射了FY-3D 业务星（2017-11），后续计划还将发射4颗卫星，包括上午星、下午星、晨昏轨道卫星、降水卫星，可实现全球每天多频次覆盖观测（唐世浩等，2016）。FY-3D 属于下午轨道卫星，过赤道升交点时间约为13：30，轨道高度平均836 km（Xu等，2018；Yang等，2019）。

经过几代人多年的努力探索，FY-3 系列卫星及其载荷也在不断成熟和稳定。MERSI 是FY-3 系列卫星的主要光学载荷之一，属于可见光红外波段扫描成像仪，可用于云、气溶胶、水汽、陆表特性、海洋水色等参数反演，为全球生态环境、灾害监测和气候评估提供重要观测数据（张鹏等，2012；杨忠东等，2017；Yang等，2019）。MERSI系列传感器总体设计扫描观测天顶角范围±67.8°，扫描幅宽约3000 km，可实现全球每天一次无缝观测，而MODIS幅宽为2330 km，在赤道附近还不能保证无缝观测。已发射的A、B、C 共3 颗卫星上搭载了第一代MERSI 传感器，最新发射的D 星（以及随后的卫星）属于业务星，其上搭载了升级版的第二代MERSI-II 传感器，在光谱响应和成像质量等性能上有明显的提高（Xu 等，2018；胡秀清等，2018）。MERSI-II在0.412—12.0 μm波段范围内增加到25 个波段。MERSI 系列传感器与MODIS 设计相似，具有用于气溶胶反演的主要波段，如引言中所述的0.47 μm、0.65 μm，2.1 μm共3 个关键波段。表1 给出了DT 算法用于气溶胶反演的MODIS 和MERSI-II 对应波段和相应的空间分辨率，以及各个波段在本文气溶胶反演算法中的用途。可以看到，除过前述的3 个关键波段之外，MERSI-II 还增加了与MODIS 相似的1.38 μm波段用于卷云滤除，0.87 μm、1.64 μm 近红外以及11.0 μm、12.0 μm 热红外波段可用于冰雪、内陆水体掩膜等。另外，MERSI 系列传感器在可见光蓝、绿、红3个真彩色通道，以及近红外波段和两个热红外的空间分辨率为250 m，其余波段均为1 km。相比之下，MODIS 仅在红光和近红外波段的空间分辨率为250 m；而MODIS在短波近红外波段为500 m，MERSI 为1km。MERSI 的这种设计将有助于更高空间分辨率的气溶胶反演和反演结果信噪比的提高。

由于可见光、近红外与短波近红外在空间分辨率上的差异影响，本文算法测试采用MERSI-II聚合后的1 km 分辨率数据集，这一点会影响后续10 km×10 km 像元聚合规则与MODIS DT 算法有一定的差别。表1 同时也可以看出，MERSI-II 的波段设置和MODIS 又存在一些差异，如：（1）相同通道光谱响应函数差异带来的中心波长差异；（2）MODIS 设置有1.24 μm 短波红外通道，而MERSI 设置在1.03 μm 处。这两点均会影响后续反演算法中查找表的重新构建及地表估计先验知识的建立。另外，MERSI-II 1 km 几何定位数据提供了经纬度、太阳与观测天顶角和方位角、高程和水陆掩膜等信息，也作为反演算法输入。

表1 FY-3D MERSI-II波段设置与MODIS对比及在本文算法中的用途说明Table 1 Spectral configuration of MERSI-II/FY-3D and MODIS and usage in the aerosol algorithm

2.2 辐射量转换

国家卫星气象中心负责MERSI 数据的总体分发，用户可以进入其网站进行免费获取。每一轨道观测数据被切割为5 min/景，数据存储格式为HDF5 格式。在进行定量应用前，需要结合定标系数进行DN（Digital Number）值到辐射值或表观反射率TOA（Top-of-Atmosphere）的转换。数据文件中数据集EV_1KM_RefSB 和EV_250_Aggr.1KM_RefSB 分别为1 km 反射通道（波段编号5—19）和250 m 反射通道（波段编号1—4）升尺度聚合为1 km后的数据DN值。式（1）给出了将DN值经过辐射定标计算转换为TOA 的具体步骤：首先利用定标系数将DN 值转换为反射率量（百分比反射率乘以100）；最后借助日地距离、太阳天顶角信息对反射率量进行校正，得到最终的TOA（Xu 等，2018；徐娜，2019）：

式中，Cal0、Cal1、Cal2分别为数据集VIS_Cal_Coeff中相应通道的定标系数（分别对应第1，2，3列）。这里需要说明的是第一代MERSI具有非线性响应，定标系数为二次多项式拟合，MERSI-II 性能得到提升后传感器响应为线性，二次项系数Cal2一般为0；DES为日地距离，取自数据文件的属性数据集；μ表示太阳天顶角余弦，天顶角信息可以从几何定位数据中获得。

数据集EV_1KM_Emissive 和EV_250_Aggr.1KM_Emissive分别为1 km发射通道（波段编号20—23）和250 m 聚合为1 km 的发射通道（波段编号24—25）的放大后辐亮度值RAD0，因此红外通道不需要使用定标系数进行辐射亮度计算。对于本文算法采用11 μm、12 μm 的热红外通道，需要计算亮度温度作为冰雪掩膜的判断条件。式（2）给出了辐射值到亮温计算的具体计算步骤（徐娜，2019）：首先对原始亮度值进行复原，得到辐亮度RAD，单位mW/（m2·cm-1·sr）；然后基于中心波长以及通道辐亮度，通过普朗克Plank 公式计算到等效黑体亮温Te；进一步利用通道亮温修正系数，将Te转换为通道黑体亮温Tbb。

式中，a和b为复原公式的斜率和截距，可从EV_1KM_Emissive 和EV_250_Aggr.1KM_Emissive 数据集的属性中得到，CW为MERSI-II 等效中心波长（表1，单位μm），A和B为修整系数，均可以从L1文件里面的属性数据集中获得；Te为等效黑体亮温；Tbb为修正后的亮温。

3 反演方法

3.1 原理

对于无云且大气水平均一的地—气系统，朗伯地表（各向同性反射）上空卫星观测的TOA 方程可以表达为（Kaufman等，1997）：

式中，μ=cos（θ），θs、θv分别表示太阳入射、卫星观测天顶角，φ表示太阳—卫星观测之间的相对方位角；ρ*(μs，μv，φ)表示卫星观测的表观反射率TOA，ρ0(μs，μv，φ)表示大气反射率（或称程辐射），ρ代表地表反射率；T（μ）代表大气直射透过率e-τ/μ和漫射透过率td之和，其中T（μs）和T（μv）分别表示下行和上行总的大气透过率；S表示大气向下的半球反照率；τ代表大气光学厚度。注：式（3）不考虑气体吸收的影响，大气散射包括大气分子和气溶胶的共同作用；第1项表示大气自身的程辐射贡献，第2 项表示地表与大气多次散射相互作用后的贡献；另外，ρ0(μs，μv，φ)，T（μ），S分别是关于气溶光学厚度AOD 和光学特性参数的函数。

式（3）作为地—气辐射传输前向模型，被广泛的用于目前的气溶胶卫星遥感反演算法中。要从式（3）反演气溶胶光学厚度AOD，还必须估计出地表反射率和气溶胶特性（类型）。如前所述，通常这两个参量通过先验知识来解决，对于单角度多光谱卫星观测，地表反射率参量一般可通过统计的波段关系（如暗目标DT 方法）或多天晴朗日合成（如深蓝DB 算法）来得到。对于气溶胶特性，通常首先会利用辐射传输代码将事先假定好的几种气溶胶类型计算成相应的查找表LUT（Look-Up Table）；对于某个反演来讲，在入射/观测几何（天顶角和方位角）给定和地表先验知识得到的情况下，通过前向模型（地—气耦合辐射传输模型）计算不同气溶胶类型、不同光学厚度下的TOA，并与观测的TOA 进行比较，将匹配误差最小的一组对应的气溶胶类型和光学厚度作为最终的反演结果。

图1为本文算法的流程图，首先需要校正依附在表观反射率中的气体（如水汽、臭氧、二氧化碳以及其他气体）吸收的影响。这是因为利用前向模型模拟计算的表观反射率却不包含气体吸收影响，不进行校正将会导致反演时卫星观测表观反射率与前向模拟表观反射率的不正确匹配。另外，从图1 可以看出，反演之前还需要去除云/水体/冰雪等不适合反演的像元；然后在10 km×10 km范围内，将气体吸收订正后的多波段表观反射率数据，进行暗像元判断与挑选，并结合云/冰雪/水体掩膜结果，进行像元的筛选与聚合，得到用于最终反演的表观反射率数据；气溶胶类型与查找表LUT的建立方法保持和DT 算法一致；地表反射率估计是本文算法进行气溶胶反演的关键，仿照DT 算法建立了适合MERSI-II 的地表估计比值模型；算法中前向模拟计算TOA 采用针对粗/细粒子比加权方式，然后通过与卫星观测TOA 比较最小进行模型的求解，最后输出气溶胶光学厚度AOD、Ångström 指数、细粒子比等反演结果，并依据反演数据输入和过程情况对反演结果加入质量标识。下面将给出各个环节详细的介绍和处理过程。

图1 本文MERSI-II气溶胶反演算法流程图Fig.1 The flowchart of MERSI-II aerosol retrieval algorithm

3.2 气体吸收订正

在气溶胶反演时，对气体（如水汽、臭氧、二氧化碳以及其他气体）吸收的订正是必不可少的一步，使得反演时卫星观测表观反射率与前向模拟表观反射率的可匹配。对于气体吸收订正，MODIS DT 算法一般事先通过大量辐射传输模拟，然后建立气体含量与气体光学厚度的关系模型，进而计算气体的吸收透过率来校正依附卫星观测表观反射率中气体吸收的影响。一般情况下，水汽和臭氧的含量被认为随时空变化较大，可从1°×1°分辨率的NCEP 再分析数据中获得；而二氧化碳和其他气体含量则采用全球平均数据代替。然后利用气体含量计算气体总的光谱吸收透过率，将光谱反射率数据除以光谱吸收透过率，实现气体吸收订正（Levy等，2013；Patadia等，2018）：

3.3 云/水体/冰雪掩膜

辐射传输模型针对的是晴朗无云情况，因此反演前需要将包含云的像元滤除。另外，在图像中存在一些“湿”像元包括冰雪、沼泽和内陆水体表现出极差的波段比值关系，不满足DT 算法的反演要求，在气溶胶反演时也需要经过各种检测后掩膜掉。

（1）云掩膜。参考MODIS DT 算法，主要依靠可见光蓝光0.47 μm 波段和近红外1.38 μm 波段的表观反射率和空间变异（像元3×3 邻域内所有9 个像元的反射率的标准差）进行云检测和掩膜（Levy等，2007）。其中，1.38 μm 波段主要用于卷云探测，对地表不敏感。如式（5）所示，当这两个波段的表观反射率或标准差大于某个阈值的4个条件中任意一个满足时，像元被判定为云：

上述云掩膜方法易使一些气溶胶高值区被当成云滤除掉，为此，C6版DT算法引入了加权平均标准差（Levy等，2013，2015），表达式为

式中，表示3×3 邻域内9 个像元的平均值，n=9（3×3个像元）。

依据这个参数，对式（5）第2 个判定条件进行了改进，当加权平均标准差且时，像元被判定为云。本文针对MERSI-II采用了这种最新的判断方法。

（2）内陆水体掩膜。MODIS DT 算法内陆水体的掩膜采用0.65 μm 红光波段和0.87 μm 近红外波段表观反射率计算的植被指数NDVI（Normalized Difference Vegetation Index），当NDVI＜0.1 时，判断为水体然后掩膜掉，不进行反演（Levy 等，2013）。值得注意的是，本研究发现这个判定条件会使得雾霾区被当成内陆水体而被掩滤除掉。这是因为基于表观反射率的NDVI 受大气条件影响而变化，当气溶胶变大时，NDVI 会变小。严重雾霾时，NDVI会小于阈值0.1。为此，本文借助短波近红外2.13 μm 波段受气溶胶影响较小、且容易被水体吸收而反射率小的特性，加入该波段的表观反射率作为限定条件（式（7）），从而实现雾霾区反演。这也是本文一个重要发现和算法改进。

（3）冰雪掩膜。MODIS DT 算法冰雪掩膜通过对0.87 μm近红外和1.24 μm短波红外波段计算的雪指数NDSI（Normalized Difference Snow Index）并设定阈值，同时结合热红外波段（11.0 μm或12.0 μm）计算的亮温值作为判断条件（Levy 等，2013）。然而，MERSI波段设置缺少1.24 μm波段，经过反复测试，可以用短波红外1.64 μm 波段替代。如式（8）所示，本文采用0.87 μm 和1.64 μm 波段的雪指数NDSI 并设定阈值，并结合热红外10.8 μm 波段计算的亮温值作为限定条件，完成针对MERSI 数据的冰雪检测与掩膜。

3.4 像元筛选与聚合

最终的气溶胶反演针对的是聚合后的10 km×10 km 数据，因此需要对这个范围内的所有像元进行筛选和平均处理，主要包括暗像元的挑选和不适合像元的去除两个方面。暗像元挑选条件为2.12 μm 反射率小于0.25，与MODIS DT 算法保持一致（Levy 等，2013，2015）；然后将云/冰雪/水体标识的像元去除，剩下的像元按照反射率大小进行排序，去掉部分较暗像元（20%），去掉部分较亮像元（50%），将剩余像元（30%）进行平均处理，以去除粗差干扰，增加输入数据信噪比。这里需要说明的是MODIS 采用的是500 m×500 m数据，在10 km×10 km 范围内有20×20 个像元，而本文采用升尺度聚合后的MERSI-II 1 km 数据，在10 km×10 km 范围内有10×10 个像元，因此信噪比相比要低一些，对反演精度会有一定的负影响。

3.5 地表估计模型

地表反射率的估计是气溶胶反演至关重要的一步，是影响气溶胶反演精度的主要因素之一。DT 算法假设短波红外2.1 μm 波段受气溶胶影响较小，然后利用统计的波段比值关系作为先验知识，实现0.65 μm 红光波段和0.47 μm 蓝光波段的地表反射率的估计（Kaufman 等，1997），地表估计线性模型如下式所示：

早期的C4版DT算法仅限浓密植被地表上空气溶胶反演，因此算法命名为浓密植被法DDV（Dense Dark Vegetation），算法中采用0.5 这一固定的地表波段比值关系。Levy 等（2007）开发了C5版（Collection 5）算法，将暗地表的条件放宽到反射率（2.1 μm 波段）小于0.25 的城市、裸土等地表，进而扩展了算法反演的覆盖范围。并借助短波红外植被指数NDVISWIR与散射角两个参数实现考虑地表类型和方向反射更为精确的地表估计，该地表波段关系模型延续至目前的C6 版算法。最新发布的C6.1 版算法又专门针对城市地表上空气溶胶反演进行了改进，主要利用地表覆盖分类图计算出城市建成区占比来确定地表波段比值关系，以解决城市地表上空比周围地类上空气溶胶反演明显高估的现象（Gupta等，2016）。

然而由于波段设置差异，MODIS 的地表估计比值关系模型不能直接应用于MERSI。Levy 等（2015）将DT算法移植到VIIRS传感器上也采用了新的地表估计模型。为此，本研究借助AERONET（AErosol RObotic NETwork）地基站网观测数据，在全球范围内裁剪站点周围一定区域内（3 km×3 km）的MERSI-II L1b 多光谱数据和几何定位数据，并加入NCEP 的水汽和臭氧等辅助数据用于气体吸收订正。以此数据集为基础，结合对应的地基精确气溶胶观测，通过对大气状况（如AOD＜0.2 的清洁天气，大气校正误差会更小）和气溶胶类型（如AE＞1.6 的细粒子）的判定、以及云/水体/冰雪掩膜等处理去掉不符合数据，然后利用辐射传输模型进行精确的大气校正，得到地表反射率数据集，最后建立适合MERSI-II 的波段比值关系地表估计模型。通过对大气校正后的地表反射率进行分析，给出了MERSI-II 地表反射率关于植被指数和散射角参数估计模型，如下式（10）—式（13）所示。对比MODIS DT 算法波段比值关系（Levy等，2007），可以看出斜率a和截距b整体具有较大差别。此外，本研究发现MODIS 0.65 μm/2.13 μm波段比值关系随短波红外植被指数NDVISWIR变化的特性，在MERSI-II上并不明显，如式（10）所示，斜率和截距为固定值；而0.47 μm/0.65 μm波段比值关系随植被指数NDVISWIR变化更加明显，如式（11）所示。对于散射角亦是如此，且本研究发现0.47 μm/0.65 μm 波段比值关系随太阳天顶角的变化更加明显，对0.47 μm/0.65 μm 波段斜率在式（11）的基础上进一步依据太阳天顶角的大小进行修正，得到式（13）。

式中，θs为太阳观测天顶角；NDVISWIR为1.03 μm与2.13 μm 波段的植被指数，这里需要注意的是MODIS的NDVISWIR是利用1.24 μm与2.13 μm短波红外通道计算得到的，而MERSI波段设置没有1.24 μm波段，这里采用1.03 μm 波段替代，这也是针对MERSI 进行需要重新建立地表估计模型的原因之一。这里需要说明的是，由于目前MERSI-II 积累的数据还较少，本研究没有像MODIS C6.1 版本算法那样对城市地表关系进行专门处理，本文算法与MODIS C6 版算法更为相似，这可能给反演精度带来一定的影响。

3.6 气溶胶类型与查找表

气溶胶类型假设是影响气溶胶反演的主要因素之一，本文采用MODIS DT算法的5种气溶胶类型：3种球形细粒子为主的气溶胶，包括非吸收型（nonabsorbing）、中度吸收型（moderately absorbing）和强吸收型（strongly absorbing）；1 种非球形粗粒子（spheroid or dust）为主的气溶胶；1 种大陆型气溶胶（Continental）。5种气溶胶类型双正态分布下的粒子谱分布及复折射指数（Complex refractive indice）等参数，请参考文献（Levy 等，2013）。注意，细粒子气溶胶类型按照季节类型被指定到全球不同区域（Levy 等，2013），算法反演时每一个地点（或像元）只能提取出一种细粒子气溶胶类型，然后与给定的一种粗粒子气溶胶类型，按照细粒子比混合出新的气溶胶类型。

气溶胶类型光学特性参数计算时，细粒子气溶胶类型被假设为球形粒子，并利用MIEV 米散射代码计算光学特性参数；对于粗粒子为主的非球形气溶胶，米散射不再适用，利用改进的T-matrix代码来计算光学特性参数（Dubovik 等，2002）；最后将计算得到的粗、细粒子散射相函数和消光系数等光学特性参数作为RT3 的输入。在查找表构建时，选择7个离散的气溶胶光学厚度作为输入（索引τ0.55=0.0、0.25、0.5、1.0、2.0、3.0 和5.0），15 个离散的太阳和观测天顶角（θs，θv=0°、6°、12°、24°、36°、48°、54°、60°、66°、72°、78°和86°），以及16 个相对方位角（φ=0°—180°，以12°的间隔递增）。然后，用RT3 辐射传输代码计算出上述所有情况下的式（1）大气程辐射、透过率等辐射传输参数，进而可建立查找表LUT（Look-Up Table）。由于MERSI-II 与MODIS 传感器光谱响应函数差异引起中心波长的差异（表1），本研究针对MERSI-II 传感器，按照上述方法，建立了0.47 μm、0.55 μm、0.65 μm 可见光和2.13 μm 短波红外4个波段的查找表。

3.7 反演求解

反演求解一般是在卫星时刻太阳入射与观测角度，不同气溶胶类型与光学厚度下，通过式（1）前向模型模拟计算出表观反射率TOA，并与卫星观测的TOA 进行比较，将匹配误差最小的一组对应的气溶胶类型（本文算法指的是细粒子比）和光学厚度作为最终的反演结果。由于卫星观测TOA 还受气体吸收的影响，因此反演前，应对卫星观测的TOA 进行气体吸收订正处理（见3.2节），和模拟计算量对应起来。

反演处理时，对于每一个像元，首先根据经纬度和时间信息从细粒子分布图上提取出细粒子气溶胶类型（Levy 等，2013），然后结合查找表，利用式（1）分别模拟计算出细粒子、粗粒子类型下3 个波段（0.47 μm、0.65 μm 和2.13 μm）的表观反射率和：

如果将卫星传感器观测到的总的表观反射率看作是由细粒子为主和粗粒子为主的气溶胶共同贡献的加权（权重η），则光谱表观反射率可表达为（Levy等，2007，2013）：

通过离散化后的不同细粒子比η（0—1，间隔0.1），前向模拟计算得到一组表观反射率，并与0.47 μm蓝光、0.65 μm红光和2.13 μm短波红外3个波段观测的表观反射率进行比较，寻找最接近的模拟值即为反演方程的解，最终反演结果为0.55 μm波段的气溶胶光学厚度τ0.55。反演结果输出时，参考MODIS DT 算法，依据10 km×10 km 范围内污染像元数、参与反演的像元数、反演求解的匹配误差等对MERSI-II 气溶胶反演结果设置了质量标识QAF（Quality Assurance Flag），值大小0—3分别表示了反演结果质量等级（Levy 等，2013）。其中，QAF=3 为质量最高的反演结果，一般被推荐给对产品精度要求较高的应用；QAF=2，1，0 分别代表反演质量较好、一般、较差3个等级，对产品精度要求不高的应用可以依据情况使用。

4 结果分析

4.1 单景反演结果

图2—4 分别给出了3 个单景反演结果示例，并选用MODIS 气溶胶产品进行反演结果对比。MODIS 分别搭载在上午星Terra（过赤道平均时10：30）和下午星Aqua（过赤道平均时13：30）上，MODIS 气溶胶产品分别被命名为MOD04 和MYD04。FY-3D 过境时间与下午星Aqua 最为接近，因此选用最新C6.1 版的MYD04 气溶胶产品与本文算法MERSI-II 的结果进行对比，以及用于随后的结果验证对比。这里，挑选MERSI-II 和MODIS 成像时刻轨道相近、观测范围重叠度大的数据，将单景AOD 反演结果进行对比，两者空间分辨率均为10 km×10 km。图2—4 分别展示了3 个典型区域结果示例：（1）案例1为2018-01-03非洲西部区域的一次观测，MODIS 成像时间为UTC 时12：45，MERSI-II 为UTC 时13：00，两者成像时间比较接近，仅相差15 min，从图2（a）的RGB彩色图可以看出，该区域有厚气溶胶存在；（2）案例2为2018-10-03 中国东部区域一次观测，MODIS 与MERSI-II 成像时间相同，均为UTC 时5：25，从图3（a）的RGB彩色图可以看出，该区域有轻度气溶胶分布；（3）案例3为2019-02-11中国东部区域一次观测，MODIS与MERSI-II成像时间非常接近，分别为UTC 时5：55 和5：50，从图4（a）的RGB 彩色图可以看出，京津冀地区有雾霾过程影响。

图2 MERSI-II气溶胶反演结果及与MODIS对比（案例1：非洲西部，2018-01-03）Fig.2 Comparison between MERSI-II and MODIS for the result of aerosol retrieval at 0.55 μm，Case one：Western Africa，Jan 3，2018

从3 个案例的MERSI-II 与MODIS 单景反演结果对比来看，由于MERSI 的扫描成像幅宽较大，因此单景反演左右扫描方向范围要大于MODIS，见图2—图4 的MODIS 结果在MERSI-II 结果图中的紫色虚线框范围。对比案例1（图2（c）（d））和案例2（图3（c）（d））的MERSI-II与MODIS反演结果，无论对于案例1 厚气溶胶或案例2 低气溶胶情况下，两者AOD 空间分布均比较相似，值域大小一致性也较好。将二者的反演结果几何校正到同一投影系统下、配准并裁剪相同区域，并生成对比散点图，如图2（b）和图3（b）所示，AOD反演值相关性非常高，相关系数R均达到0.9以上。从图4 案例3 的反演结果对比来看，MODIS 气溶胶产品将雾霾比较严重的区域丢失，在大部分雾霾区域几乎没有反演结果（图4（b）），而本文算法可实现MERSI-II 对该区域的反演（图4（c））。同时，本文结果没有引入更多的噪声点，说明我们找到了雾霾漏反演的根本原因，且关于内陆水体掩膜的改进方案是可行的。

图3 MERSI-II气溶胶反演结果及与MODIS对比（案例2：中国东部，2018-10-03）Fig.3 Comparison between MERSI-II and MODIS for the result of aerosol retrieval at 0.55 μm（Case two：Eastern China，Oct 3，2018）

图4 MERSI-II气溶胶反演结果及与MODIS对比（案例3：中国东部雾霾，2019-02-11）Fig.4 Comparison between MERSI-II and MODIS for the result of aerosol retrieval at 0.55 μm（Case three：Haze over Eastern China，Feb 11，2019）

4.2 基于地基站网验证分析

本文选用2018年10月、2019年2月、2019年5月的3个月MERSI-II全球白天观测数据进行反演测试实验，平均约155 景/d，共计1.4 万景数据。对MERSI-II 反演结果的真实性检验采用全球气溶胶自动观测站网AERONET 产品，该站网仪器以CE318 为主。目前，该站网产品已经升级到第3 版（Version 3），并完全取代早期的版本Version 2 和Version 1（Sinyuk等，2020）。在每个版本中又分为不同级别的产品：Level 1.0为初级产品，没有经过任何筛选，质量等级较低；Level 1.5为更高一级产品，经过了自动云去除和质量控制，质量等级较高；Level 2.0为加入人工筛选且经过严格质量控制产品，质量等级最高，本文选用该产品对反演结果进行验证。其中，Level 1.0和Level 1.5计算机自动处理产品，更新速度快，平均1 次/周，数据的现势性较好；而Level 2.0 加入人工判断，数据具有一定的滞后性，这也是本文选用更早时期的MERSI-II 观测数据进行反演测试的因素之一。最终全球有约有450 个站点参与了对MODIS 和MERSI-II气溶胶反演结果的验证与对比分析。

在进行卫星气溶胶结果验证时，传感器波段设置通常会与地基CE318 观测数据波段不对应，导致无法直接进行验证。这时，需要借助地基多波段AOD 观测值插值得到对应卫星结果相应波段处的AOD。Eck 等（1999）研究发现AOD 和波长二者的对数之间具有较好的二次多项式关系，如式（16）所示。可利用多波段（至少3 个）AOD观测，求解出多项式的系数，进而将卫星传感器相应波段代入可得到，对应波段的AOD 结果（Levy 等，2010）。一般验证时，多展示0.55 μm 波段处的AOD验证结果。

式中，λ表示地基观测相应的波长，τλ为对应波长的AOD 值；a0、a1、a2为二次多项式的系数。为了得到更稳定的拟合二次项系数的解，本文采用了地基观测0.44 μm、0.675 μm、0.87 μm、1.02 μm四波段参与插值运算，主要由于AERONET 站网内各种型号的CE318仪器均具有这4个波段观测。

为了对比分析，本文同时也对与MERSI-II 测试数据对应的3 个月的全球MODIS 气溶胶产品进行验证。如前所述，两种数据反演结果均为10 km×10 km，且数据组织方式相似，均切割为5 min/景。因此，本研究采用相同的验证方法。首先在时空上进行匹配，时间上：提取卫星过境时间前后30 min 内2 个以上的地基观测数据，并取平均值；空间上：以地面站点所在位置为圆心，25 km 为半径，对落在范围内的3 个以上的反演结果取平均（Levy 等，2015；Sawyer 等，2020）。然后将对应的结果生成散点图，并计算常用的用统计量（相关系数，均方根误差、平均偏差）来衡量反演精度，式（17）—（19）给出了具体计算公式。

（1）相关系数：

（2）均方根误差：

（3）平均偏差：

式中，AODs、AODg为匹配的卫星和地基气溶胶结果，分别为卫星和地基AERONET 所有匹配点的平均值。

另外，进一步计算得到散点落入期望误差EE（Expected Error）范围之内、之上、之下的百分比，对应图5的%within EE、%above EE、%below EE，也可作为精度衡量参数。用落在EE 之上、之下的点数比例对比，可简单地判断反演结果偏高或偏低，用落入EE 之内的点数比例来判断总体精度高低。本文取EE=±（0.05+0.15τ），通常落入点数比例大于等于2/3，认为达到期望误差要求（Levy等，2010，2013）。

为了更好地评定MERSI-II气溶胶反演的精度和潜力，参考MODIS气溶胶产品验证分析经验（Levy等，2010，2013，2015；Sawyer等，2020），本文仅对两种数据质量标识等级最高QA=3 的反演结果进行验证及对比分析。验证结果散点图如图5所示，图5（a）为对MODIS下午星Aqua气溶胶产品0.55 μm波段的验证结果，图5（b）为MERSI-II 0.55 μm波段验证结果。对比图5（a）和（b）可以看出，MERSIII的匹配点数为5008，高于MODIS 匹配点数4167（约高出20%左右），这主要与MERSI成像幅宽大于MODIS的优势、以及与本文算法对内陆水体掩膜的改进有关，另外MODIS更严格的质量控制也对此会有一定的影响。MERSI-II 反演结果整体验证精度较好，相关系数R 高达0.866；均方根误差、平均偏差也较小，说明与地基观测比较接近；落入期望误差范围EE=±（0.05+0.15τ）点数比例为65.14%，基本达到2/3的要求，说明本文MERSI-II的气溶胶反演结果已接近国际同类型气溶胶产品精度要求。同时，对比MODIS 验证结果，MERSI-II 验证点的整体离散程度要更大一些，特别是AOD低值区域有一些反演过高的噪声点，这可能与本研究对基于NDVI 的内陆水体掩膜改进造成更多污染像元的处理还不够精细有关；另外，MODIS验证结果的相关系数、均方根误差、平均偏差以及落入EE 的比例等因子均优于本文MERSI-II 验证结果，这可能与本研究采用的数据分辨率、像元筛选与聚合、地表估计以及质量控制有关，将在第5节进行更细致地讨论。

图5 MERSI-II气溶胶光学厚度0.55 μm波段地基验证结果及与MODIS对比Fig.5 Comparison between MERSI-II and MODIS for AOD validation against AERONET at 0.55 μm

4.3 月平均结果对比

对上述3 个月MERSI-II 全球反演测试结果，将每一景结果（10 km×10 km，约为0.1°×0.1°）经过几何校正、拼接等处理，在0.1°×0.1°的基础上，根据10×10 范围内有5 个以上的有效反演像元的规则，取平均处理聚合成1°×1°的日合成产品（Sawyer等，2020）。进一步当一个月时间范围内有3 个（或3 d）以上的有效值时取平均，生成分辨率1°×1°的月平均结果（如图6（b）所示）。同时图6（a）给出了MODIS 对应月份的月平均结果。图6（a）、（b）、（c）（审图号：GS（2021）5922号）分别表示MODIS月平均、MERSI-II月平均和结果对比散点图；第1、2、3 行图分别表示2018年10月、2019年2月、2019年5月共3 个月份。对比图6（a）和（b）系列图的3 个月份MERSI-II 与MODIS月平均结果，从AOD 值域高低和空间分布来看，整体均具有较好地一致性。进一步将二者月平均结果生成对比散点图，如图6（c）所示，3个月份结果的线性相关性均非常高，相关系数R 均达到0.93 左右。再一次印证本文算法MERSI-II 气溶胶反演的定量能力和潜力，以及结果已接近国际同类气溶胶产品。

图6 MERSI-II气溶胶光学厚度月平均结果及与MODIS对比（审图号：GS（2021）5922号）Fig.6 Comparison between MERSI-II and MODIS for monthly mean AOD at 0.55 μm

但同时对比结果也存在一定的差异：（1）对比图6（a）和（b）系列图，MERSI-II月平均结果高值区整体高于MODIS，反映在结果对比散点图中低值区更接近1∶1线（图6（c）系列图中黑线），高值区偏离1∶1线倾向上方，这可能由于本文算法改进内陆水判别方法，实现更多气溶胶高值区（如雾霾）的反演，从而引起MERSI-II高值区高于MODIS；（2）MERSI-II月平均结果在北纬高纬（尤其是北美洲）气溶胶分布边界处（对比图6（a）和（b）系列图），有明显的高值边界，这可能与本文算法冰雪掩膜不够严格有一定的关系，还需要进一步完善；（3）在赤道附近的沙漠地区，出现个别离散高值点，本文算法改进内陆水体判别方法实现雾霾天气反演的同时引入了一些噪声点，因此后续工作还需要对该方法做进一步完善。

5 结论

气候与空气质量监测对卫星遥感气溶胶定量产品的需求非常迫切，然而中国仍然没有摆脱依靠国外资料的现状，如被广泛使用的MODIS 气溶胶产品。风云三号卫星搭载的MERSI 与MODIS 属于同类型传感器，具备气溶胶反演的相应通道，且具有全球每日覆盖的稳定业务观测能力，但目前MERSI 还没有精度可靠、稳定的、全球适用的气溶胶业务化产品。进一步提高风云气象卫星的定量化应用水平以及国际影响力和竞争力，成为当前和今后的一个重要发展目标。本文面向这一需求，借鉴MODIS DT算法思想，针对FY-3D星搭载的MERSI-II 传感器数据，在传感器波段设置差异、地表估计模型、云/冰雪/水体检测、像元筛选与聚合等方面进行专门改进与优化，发展了针对MERSI-II 全球适用的陆地气溶胶反演算法。为测试MERSI-II 传感器气溶胶反演的定量能力，并与MODIS 气溶胶产品具有可对比性，本文算法这里在各方面处理尽可能与DT算法保持一致。

将本文算法MERSI-II 气溶胶反演结果进行地基验证、并与MODIS气溶胶产品单景、月平均进行对比，结果表明：（1）MERSI-II 单景反演结果与MODIS 具有较好的一致性，相关系数R 达到0.9 以上；（2）本文算法改进基于NDVI 的内陆水体判别方法后，可实现雾霾高值区有效反演，而MODIS在该区域漏反演；（3）MERSI 具有幅宽大于MODIS的优势，以及加上本算法对内陆水体判别方法的改进，使得MERSI-II 可反演有效像元数整体要高于MODIS，表现在单景反演结果范围和验证散点图匹配点数具有多于MODIS结果；（4）利用AERONET地基站网观测数据对MERSI-II全球3个月的反演结果进行验证，整体具有较高的验证精度，落入期望误差EE=±（0.05+0.15τ）点数比例为65.14%，基本达到2/3 的比例要求（Levy 等，2010；Sawyer 等，2020）；（5）通过对3个月份MERSI-II与MODIS 的月平均结果对比，空间分布一致性较好，结果的线性相关性非常高，相关系数R 均达到0.93 左右。总体说明本文算法气溶胶反演结果已接近国际同类型产品的精度要求，MERSI 具有较好的定量能力和应用潜力，将对全球长时间序列气溶胶观测是一个重要补充。

同时，与MODIS 验证结果相比，本文算法MERSI-II 反演验证结果各方面精度指标整体略低于MODIS。可能存在以下几个方面的原因：（1）本研究采用的是在轨测试初步得到的定标系数，随着仪器趋于稳定和数据积累增加，期待优化后的定标系数对反演结果精度提高有所帮助；（2）MODIS气溶胶产品反演时采用500 m 数据，通过20×20 像元聚合到10 km×10 km（Levy 等，2013），而本文MERSI-II 反演测试采用1 km 数据，通过10×10 像元聚合到10 km×10 km，因此信噪比低于MODIS算法，验证散点图离散度要大于MODIS 验证结果，后续做产品时，采用原始250 m 分辨率有望提高反演精度；（3）本文算法改进内陆水判别方法，可实现雾霾等高值区气溶胶反演，但同时这种改进可能会引入一些NDVI 较小、而不适合反演的点（如北纬高纬度的雪融地区，赤道附近的沙漠区），这些区域地表关系模型适用性变差，故MERSI-II月平均结果中出现一些离散高值点，表现在验证结果散点图的低值区也存在一些离散高值点，这一方面在后续工作中还需要进一步加入限定条件进行优化；（4）如前所述，MODIS C6.1版算法针对城市地表设置了专门的比值估计模型（Gupta 等，2016），由于目前MERSI-II观测积累较少，本文算暂未对该方面进行特别处理，这也是影响精度差异的主要原因之一，下一步在更多观测数据积累的基础上，进一步优化地表估计模型；（5）从4.3 节验证散点图的对比来看，产品的质量控制也会是产生验证结果差异的因素之一，MODIS 气溶胶产品在质量控制方面更为精细和严格，验证结果散点图中噪声点更少，下一步需要在这些方面进行更精细的改进和优化。

志 谢本研究使用了AERONET 提供的全球地基站点数据用于结果验证，使用了NASA LAADS DAAC 提供的MODIS 气溶胶产品用于结果对比分析。在此表示衷心的感谢！