资源一号04星WFI数据气溶胶反演与验证——以北京市及北京周边地区为例

丁宇，汪小钦，王峰

(1. 数字中国研究院(福建)，福州 350108；2.卫星空间信息技术综合应用国家地方联合工程研究中心，福州 350108；3.福州大学 空间数据挖掘和信息共享教育部重点实验室，福州 350108)

0 引言

气溶胶是大气层的重要组成成分，对地气系统、辐射平衡、气候变化有着显著影响[1]，高浓度的气溶胶也会加大人体患病的几率[2]。传统的气溶胶观测手段主要以地面实时监测为主，难以满足大范围监测和时空分布的要求，而遥感监测手段弥补了这一不足，能够对大气污染进行大范围的动态监测[3]。

大气气溶胶遥感实质上是研究如何实现地气解耦，而实现地气解耦的前提是得到准确的地表反射率信息。为此，研究者提出了暗像元[4]、偏振[5]、双星协同[6]、深蓝[7-8]等算法计算地表反射率，开展气溶胶反演研究。其中，暗像元法和深蓝算法在国内外MODIS[9-10]、OLI[11-12]、HJ-1[13-14]、GF-1[15-16]、FY-3[17]、Himawar-8[18]等传感器上得到广泛应用。暗像元法只对浓密植被等暗地表效果良好，深蓝算法主要针对城市、沙漠等高亮地表。为了解决不能同时反演包含明暗像元的复杂地表上空气溶胶光学厚度(aerosol optical depth，AOD)的问题，美国国家航空天局(National Aeronautics and Space Administration，NASA)最新发布的MODIS气溶胶产品C6版本包括了融合2种算法的产品。文献[19-20]将其与单一算法的MODIS气溶胶产品对比，发现融合后的产品更具空间完整性。国内学者结合2种算法在局部区域开展气溶胶反演研究[21-23]，也得到了相似的结论。

资源一号04星(CBERS04)由我国与巴西联合研制，并于2014年12月7日在太原卫星发射中心发射成功。其搭载的宽视场成像仪WFI有B1(0.45～0.52 μm)、B2(0.52～0.59 μm)、B3(0.63～0.69 μm)和B4(0.77～0.89 μm)共4个波段，空间分辨率为73 m，幅宽为866 km，重访周期高达3 d[24]。CBERS04WFI的时间分辨率和幅宽均优于HJ-1 CCD、GF-1 WFV、OLI等传感器，空间分辨率又优于MODIS、FY-3 MERSI、Himawar-8 AHI等传感器，但是目前利用WFI数据的应用较少，还未见到应用CBERS04WFI数据开展AOD反演的公开报道。

WFI在时空分辨率和幅宽上的优势，使其在气溶胶估算上具有较大的应用潜力。本文结合暗像元法与深蓝算法，探讨基于动态查找表的WFI数据AOD反演方法研究，为区域大气变化监测提供较高时空间分辨率的数据支持，同时扩展国产卫星的应用领域，拓展应用广度。

1 研究区和数据

1.1 研究区

本文以北京市及周边地区(115°25′E～117°30′E，39°26′N～41°03′N)作为研究区。研究区位于华北平原西北部，研究区东南为平原，西北为山地，总体地势为西北高、东南低(图1)。研究区气候属暖温带半湿润半干旱季风气候，年降水量时空分布不均匀，季节上夏多冬少，东北部和西南部山前迎风坡地区为相对降水中心，西北部和北部深山区降水较少[25]。因此一年中北京6—8月份的大气污染物的稀释和扩散较快，大气环境质量相对最好；冬季受气象条件和采暖的影响，春季受西北部沙尘暴的影响，大气环境质量相对较差[26]。

图1 研究区海拔及AERONET站点分布

1.2 数据

研究所用的数据包括CBERS04WFI数据、MODIS地表反射率数据和气溶胶产品，以及AERONET地基站点观测数据。

从中国资源卫星应用中心网站下载的2016—2018年113景CBERS04星WFI数据，包括tif文件、观测几何txt文件，和提供定标参数等信息的xml文件。

2016—2018年的MODIS产品，包括MOD09A1数据和MODIS气溶胶产品。MOD09A1数据是8 d合成的地表反射率产品，尽量降低了反射率短期变化以及云覆盖的影响[27]。MODIS气溶胶产品包括MOD04_3K(暗像元法，空间分辨率3 km)和MOD04_L2(深蓝+暗像元，空间分辨率10 km)，为获得较高时空分辨率的气溶胶产品，将2种产品进行融合，融合原则是以MOD04_3K数据为主，用MOD04_L2数据补充缺失区域。

采用AERONET站点Level 2.0级数据验证反演结果。因AERONET数据没有与反演结果对应的550 nm处AOD数据，利用二次多项式插值得到550 nm处AOD值作为验证数据[28]，如式(1)所示。

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

(1)

式中：τλ表示λnm处的AOD值；a1、a2为相应项未知系数。

2 研究方法

2.1 基本原理

假设地面为均匀的朗伯面且大气条件均一，卫星传感器的表观反射率为ρTOA，如式(2)所示。

(2)

式中：ρO为大气的路径辐射项等效反射率；ρs表示地表二向反射率；φ为相对方位角；μs=cosθs，μv=cosθv，μs与μv为太阳天顶角和观测天顶角；S为大气下界的半球反射率；T(μs)、T(μv)分别为大气分子上下行散射透过率[29]。

通过式(2)可以明显看出，要想获得准确的大气路径辐射项等效反射率，需要计算得到准确的地表反射率。因此，本研究结合深蓝算法和暗像元法，对亮暗像元进行区分，分别计算其地表反射率值。其中，暗像元法是根据浓密植被等暗地物在红蓝波段与短波红外波段的地表反射率间的线性关系实现气溶胶反演。但是由于众多传感器不具有短波红外波段，无法使用该算法进行AOD反演。因此许多研究基于红蓝通道固定地表反射率比值法进行气溶胶反演，式(3)为NASA的MODIS 数据V5.2气溶胶算法的红蓝通道地表反射率关系式[30]。

y=0.49x+0.05

(3)

式中：y为红光波段反射率；x为蓝光波段反射率。因WFI传感器缺失短波红外波段信息，选择利用红蓝波段间线性关系开展WFI传感器气溶胶反演研究。为得到WFI传感器红蓝波段间线性关系，将WFI和MODIS红蓝波段的光谱响应函数与ENVI光谱库中220种植被光谱(重采样为1 nm)进行卷积，计算WFI和MODIS红蓝通道植被地表反射率并进行回归分析(图2)。结合图2中拟合方程式和式(3)，计算得到WFI传感器红蓝波段地表反射率关系模型，如式(4)所示。

图2 MODIS与WFI红蓝波段地表反射率散点图

ρ3=0.500 4ρ1+0.046 8

(4)

式中：ρ1和ρ3分别为WFI蓝光与红光地表反射率。

针对高亮的城市、沙漠等区域，本研究选择采用的深蓝算法是依据蓝光波段气溶胶对卫星观测信号有较强贡献而地表反射较弱的特点，借助地表反射率库先验知识，以实现地气解耦。本文选择使用MOD09A1数据来构建地表反射率库，但不同传感器对相同地物成像得到的地表反射率会有差异。研究表明，0.01的地表反射误差带来约0.1的气溶胶光学厚度误差[31]。因此，需要通过波段修正的方法将MOD09A1数据转换为WFI地表反射率数据。选用ENVI光谱库中各类地物光谱共426条，分别与MODIS和WFI蓝光光谱响应函数卷积，计算光谱等效值并进行回归分析(图3)。利用二者拟合方程式将MOD09A1数据转换为WFI地表反射率数据。如图3所示，光谱修正公式与直线y=x极其接近，因此转换得到的WFI地表反射率与经过质量验证的MOD09A1地表反射率产品只有极小的差异，但经过光谱修正得到的WFI地表反射率更符合WFI传感器光谱特性。

图3 MODIS与WFI蓝波段地表反射率拟合图

2.2 WFI数据气溶胶反演流程

遥感反演气溶胶光学厚度的主要过程是将相关先验知识，如观测几何、气溶胶模型、观测时间等，迭代循环输入6SV辐射传输模型，计算得到3个大气参数ρo、S和T。结合暗像元法或深蓝算法计算得出的地表反射率参与辐射传输计算，得到最终AOD值，技术路线如图4所示。

图4 CBERS04星WFI数据气溶胶反演流程图

1)辐射定标。利用原始数据中的xml文件自带的辐射定标系数、太阳高度角和成像时间等参数，结合大气层外太阳辐照度[32]，将影像DN值转换为表观反射率，如式(5)所示。

(5)

式中：ρ为表观反射率；Gain为增益值；Bias为偏置值；D为日地距离；ESUNb为大气层顶太阳辐照度；θ为太阳天顶角[33]。

2)云水检测。影像上的云、水像元会对反演结果造成误差，需要在实验前对云、水进行检测。将归一化水体指数NDWI大于0.3的像元判定为水体像元[34]，将红光波段大于0.25且归一化植被指数NDVI[35]<0的像元判断为云像元[23]。

3)亮暗像元判断。通过增强型植被指数EVI[36]区分亮暗像元。将EVI大于0.4的像元判定为暗像元[37]，暗像元区域采用暗像元法反演AOD；EVI小于等于0.4的亮像元区域则采用深蓝算法(式(6))。

(6)

式中：ρ1为蓝光波段值；ρ3为红光波段值；ρ4为近红外波段值。

4)动态查找表构建。传统查找表只针对单个传感器，其构建方法是事先设置好参数组合，通过大气辐射传输模型计算得到，可移植性不强。动态查找表是直接在反演过程中，根据每一景影像的几何参数生成只针对该影像的查找表。查找表构建成功后，以临时文件的形式保存在计算机内存中，查找效率远超以文本形式保存的传统查找表。动态查找表可以输入多个传感器的光谱响应函数，同时适用多个传感器，具有较强的可移植性。构建方法：①将影像太阳、卫星天顶角数据四舍五入为整数值(WFI传感器的观测几何数据中这2个角度覆盖范围最大)，按照角度递增顺序划分区域；②对2个角度文件划分区域取并集，得到最终划分区域；③逐个区域计算观测几何参数，将预设AOD值(550 nm处AOD步长设置范围为0～1.95之间，步长为0.05)，以及气溶胶模型输入6SV辐射传输模型计算得到大气参数数据。本文选择能够较准确描述研究区气溶胶微粒状况的大陆型气溶胶模型，作为6SV模型输入条件，该模型提供了所需的散射相函数、单次散射反照率和不对称因子等关键参数[38]。

3 结果与分析

按照上述算法与数据处理流程，反演得到研究区空间分辨率为500 m×500 m的AOD分布。

3.1 气溶胶反演方法对比与分析

为分析本文方法应用于反演复杂地表上空AOD的性能，分别与暗像元法和深蓝算法结果进行比较(图5)。从图中可以发现：暗像元法较好地反演了浓密植被覆盖的山区，但无法反演亮像元的城市等区域；深蓝算法能够有效反演亮地表上空AOD，但对于其他区域的反演结果偏高，其反演结果均值为0.223 1，而本文方法反演结果均值为0.177 7，说明深蓝算法反演结果整体偏高，这也与张胜敏等[39]的研究结论一致。究其原因是因为深蓝算法使用的MOD09A1地表反射率产品是8 d数据合成构建而成的，这一构建方法会低估地表反射率，进而在反演中高估大气反射率，最终造成反演结果偏高[27]。本文方法有效集成了暗像元法和深蓝算法的优势，反演结果较单一算法更适合用于反演包含明暗像元的复杂地表上空AOD。

图5 2018年9月21日的反演结果

3.2 与MODIS产品对比

为较为全面地验证反演精度，本研究使用融合后的MODIS气溶胶产品与反演结果进行空间分布的对比。如图6所示，反演得到的 AOD与MODIS气溶胶产品高、低值分布情况基本一致，但覆盖范围更全；MOD04_3K(暗像元法)和MOD04_L2(深蓝+暗像元)气溶胶产品，都会在生产过程中去除不合适的像元(如云、沙漠、雪/冰和内陆水)，同时也剔除一部分最暗和最亮像元[40]，因此，融合后的MODIS产品还是会存在较多无效值，在覆盖范围上不如WFI AOD；AOD值最高的是研究区东南部的通州、大兴、廊坊等地，中部区域次之，西北部房山、门头沟、昌平、延庆、怀柔、密云等区域AOD值最低，即研究区AOD分布总体上呈现自东南城市区域向西北山区递减的趋势。与图1对比发现，AOD分布情况与海拔存在明显的对应关系，人类活动多的低海拔平坦区域AOD较高，山区AOD较低。

图6 WFI反演AOD与MODIS AOD空间分布对比图

3.3 基于AERONET站点监测数据的精度验证

本文采用AERONET站点Level 2.0级数据验证反演结果。为客观评价研究模型反演结果的精度，本文采用了相关系数(r)、均方根误差(RMSE)、平均绝对误差(MAE)和预期误差区间(expected error，EE)这4个统计指标对其进行精度评定。其中，预期误差区间EE的计算公式如式(7)所示[41]。

EE=±(0.05+0.2·τα)

(7)

式中：τα为AERONET站点实测AOD值。

如图7所示，本研究共得220个有效数据对，反演结果与AERONET站点实测值相关性较高(r>0.94)，RMSE为0.170 3，MAE为0.118 3，67%的反演结果在误差区间内，29%的反演结果高估，4%的反演结果低估。

图7 反演结果与AERONET AOD散点对比图

高估现象主要发生在低值AOD处，即AERONET站点AOD小于0.4时。这一现象主要有以下2点原因造成：①地表反射率相差0.01，AOD会出现约0.1的偏差[31]。因此在实际AOD值较低时，地表反射率微弱的偏差都可能导致反演结果较大的偏差。②高估的原因可能与MODIS地表反射率数据构建方式有关。MOD09A1地表反射率数据的构建尽量考虑高观测覆盖、低视角、无云及云的阴影以及气溶胶浓度的影响，选择8 d中最可能、最稳定的数据构建而成[27]。而在大气质量较好时，通过这种方式构建出的地表反射率数据，可能会因为考虑太多因素导致反射率略微偏低，从而导致气溶胶反演结果偏高。

综上所述，利用本文方法反演气溶胶光学厚度具有较高稳定性，反演结果的空间分布和精度都符合要求。

4 结束语

本文基于CBERS04星WFI数据特点，结合暗像元法和深蓝算法，实现了基于动态查找表的AOD估算。本文利用WFI数据结合暗像元法和深蓝算法较好地反演了包含亮暗像元的复杂地表上空AOD，避免了单一算法的不足。本文反演结果与MODIS气溶胶产品和AERONET站点AOD进行了对比验证，证明本文研究方法能较好地反映气溶胶的空间分布。反演结果空间分辨率和空间覆盖都优于MODIS气溶胶产品，且与地面站点实测值相关性显著(r>0.94)，但在空气质量较好时反演结果容易高估，这与地表反射率数据的精度及其构建方式有关。