牛晓君¹, 唐家奎¹, 张自力², 崔林丽³, 幸炜鹏¹, 宋艺¹
1. 中国科学院大学资源与环境学院, 北京 100049;
2. 浙江省环境监测中心, 杭州 310015;
3. 上海市气象科学研究所, 上海 200030
2018年3月22日 收稿; 2018年5月22日 收修改稿
基金项目: 中国科学院战略性先导科技专项资助（XDA20050103）和国家自然科学基金（41571352）资助
通信作者: 唐家奎，E-mail:jktang@ucas.edu.cn

摘要: 雾霾灾害的发生具有时空动态变化特征。葵花8号卫星是日本气象厅发射的新一代静止气象卫星，具有观测范围广、观测频率高、空间分辨率高等特点，非常适合雾霾过程的动态监测。本文提出利用葵花卫星可见光波段数据及地表反射率比值反演气溶胶光学厚度的新算法，并在北京区域选择2017年5月17日和18日两个雾霾日进行反演实验。反演结果分别与国际AERONET 3个站点Beijing、Beijing-CAMS和Beijing-PKU站观测数据，葵花官方气溶胶产品，MODIS官方气溶胶产品对比分析。结果表明，本文提出算法的反演结果，在反演精度、有效反演空间范围、有效反演时间范围、空间分辨率等方面均为最优。最后通过对北京区域两个雾霾日小时级气溶胶AOD反演结果的时空分析，探讨雾霾过程的时空动态变化。从AOD反演结果时间序列图像中，清楚地看到雾霾污染中心的移动方向及加重或减轻的动态变化状况，说明应用葵花卫星数据反演气溶胶对雾霾过程的监测具有较大的应用潜力。
关键词: 葵花8号卫星气溶胶雾霾光学厚度
The method of aerosol retrieval using Himawari-8 satellite data and its application in monitoring haze process
NIU Xiaojun¹, TANG Jiakui¹, ZHANG Zili², CUI Linli³, XING Weipeng¹, SONG Yi¹
1. College of Resources and Environment, University of Chinese Academy of Sciences, Beijing 100049, China;
2. Zhejiang Province Environmental Monitoring Center, Hangzhou 310015, China;
3. Shanghai Institute of Meteorological Science, Shanghai 200030, China


Abstract: Haze disaster is characterized by spatio-temporally dynamic changes.Himawari-8 is a new generation of geostationary meteorological satellite launched by Japan Meteorological Agency, and it is very suitable for haze process dynamic monitoring due to its observational ability in wide range, high frequency, and high spatial resolution. A new algorithm was proposed to retrieve aerosol optical thickness using the ratio of Himawari-8 visible band surface reflectivity, and the retrieval evaluation experiments were carried out on May 17, 2017 and May 18, 2017, the two haze days, in Beijing area. The retrieval results showed that the proposed method is the best in retrieval accuracy, retrieval effective spatial and temporal ranges, and spatial resolution by comparing the results with the observation data of the three international AERONET sites of Beijing, Beijing-CAMS, and Beijing-PKU, the Himawari-8 official aerosol products, and the MODIS official aerosol products. Finally, based on the analysis of hourly-level AOD retrieval results on the two haze dates, temporal changes of the haze process in Beijing area were discussed. Based on the retrieval images of the time series, we clearly observe the moving direction of the pollution center and the increasingly or decreasingly dynamic changes of the haze, which shows that the proposed Himawari-8 algorithm has potential application for the dynamic monitoring of haze process.
Keywords: Himawari-8 satelliteaerosolhazeAOD
近年来，随着经济和城市化迅猛发展，中国雾霾污染日趋严重, 雾霾遥感监测有迫切的应用需求。传统的雾霾监测手段主要是利用稀疏的地基站点进行PM_2.5测量，仅能获得观测点处的空气污染状况，对于大范围的雾霾监测，卫星遥感具有较强优势^[1]。中国雾霾灾害往往具有季节性、区域性特点，具有较多可见光谱段的极轨卫星已经初步展示雾霾监测的能力，以TERRA和AQUA卫星为例，每天能实现1~2次全球覆盖观测，搭载的MODIS传感器获取的气溶胶光学厚度(aerosol optical depth, AOD)产品是目前常用雾霾监测的指标^[2]。然而, 雾霾的主要成分气溶胶的生命周期短、时空变化强，每日1~2次的极轨卫星观测无法满足对其性质、动态变化等特性观测的需求，尤其对于雾霾的发生、发展、扩散迁移、消散等过程的监测，这需要更高频次的观测。静止卫星具有高时间分辨率的特点，可以对地面进行定点凝视观测，观测频率可达小时级甚至分钟级，因此，静止卫星对于雾霾动态过程的监测具有明显优势^[3]。
目前，雾霾(PM_2.5)卫星遥感监测工作主要是基于搭载卫星上的气溶胶监测传感器开展的，在卫星遥感反演气溶胶的基础上，再进一步进行雾霾监测分析^[1]。因此，气溶胶光学厚度反演成为雾霾遥感监测的关键环节。国内外已有****开展了利用静止卫星监测气溶胶光学厚度及雾霾灾害的研究工作。毛节泰和刘莉^[4]开展用GMS5卫星资料反演湖面上空气溶胶的可行性研究，结果表明，通过选取合适的湖面反射率, 卫星反演的气溶胶光学厚度和地面光度计遥感的月均值相对误差不超过30%。张军华等^[5]利用GMS卫星遥感中国大陆25个湖面上空气溶胶光学厚度，结果指出, 利用GMS卫星反演湖面光学厚度可以提供大陆上空气溶胶的信息。高玲等^[6]利用日本静止卫星MTSAT数据采用最小卫星反射率反演地表反射率的方法反演气溶胶光学厚度，并与AERONET观测及MODIS气溶胶光学厚度产品进行比较，结果表明，MTSAT反演的气溶胶光学厚度产品可以反映大气气溶胶光学厚度的日变化信息。任通等^[3]探讨利用中国风云2C静止卫星可见光资料反演气溶胶光学厚度的数值方法，研究结果表明，在东亚地区利用风云2C可见光资料反演的AOD产品可以反映气溶胶的分布，但算法存在不同地区AOD值的高估和低估现象。陈健等^[7]应用韩国静止卫星GOCI数据，采用国际深蓝算法进行长江三角洲地区气溶胶光学厚度的遥感反演研究，发现基于静止卫星GOCI数据的气溶胶产品具有监测气溶胶日变化的能力。Moulin等^[8]利用欧空局静止卫星Meteosat进行气溶胶反演，并计算1983—1994年期间由撒哈拉沙漠输送到大西洋和地中海上空的沙尘量。Wang等^[9]利用静止卫星GOES8反演大西洋区域的沙尘气溶胶，并与地面观测试验进行比对，结果表明，高时间分辨率的静止卫星可以有效弥补极轨卫星全球气溶胶遥感的不足。Prados等^[10]评估NOAA/NESDIS利用静止卫星GEOS可见光通道反演的气溶胶/烟尘产品(GASP)，通过与地基AERONET和MODIS官方产品比对，表明反演结果在北美地区具有较好的精度。
日本气象厅于2014年7月发射新一代静止气象卫星Himawari-8，Himawari-8可实现时间分辨率为10min/次的高频次对地观测，引起广泛关注。目前已有利用该卫星传感器应用于天气预报、海面温度、灰霾检测等方面的相关研究^[11-13]。Himawari-8搭载的AHI传感器具有对气溶胶敏感的可见光通道。目前已有****利用Himawari-8数据进行气溶胶遥感和雾霾监测的尝试：Daisaku等^[14]利用Himawari-8双通道(可见光0.64μm和近红外0.86μm波段)数据进行气溶胶反演试验；Sekiyama等^[15]将Himawari-8反演的气溶胶光学厚度用于全球气溶胶预测模型Kalman滤波系统，并预测了东亚的一次沙尘暴。葛邦宇等^[16]利用Himawari-8的可见光和近红外波段进行了气溶胶反演，并验证了用暗像元法反演葵花8号卫星数据得到的结果精度及方法的可行性。Wang等^[17]利用葵花卫星官方气溶胶产品，采用线性混合效应模型对2015年7月至2017年5月京津冀地区的PM_2.5相关分析建模，结果表明葵花卫星能够有效用于PM_2.5监测。
本文提出一种利用葵花卫星反演气溶胶光学厚度的算法，并将其结果与AERONET、MODIS等观测的气溶胶光学厚度值进行对比，对反演结果进行误差分析。最后以气溶胶光学厚度为雾霾指标，对北京区域一次雾霾过程进行反演分析。
1 卫星及数据葵花8号卫星(Himawari 8)是由日本气象厅于2014年发射的气象卫星(http://www.eorc.jaxa.jp/ptree/userguide.html)，经过性能检测后，于2015年正式投入使用。葵花8号卫星的星下点位于东经140.7°的赤道上空，距地高度约为35800km，搭载有AHI成像仪，该传感器共有3个可见光通道和13个近红外及红外通道，最高分辨率可达0.5km。AHI全盘扫描周期为10min，如选定时间对特定区域进行扫描，可达2.5min观测周期。相对于极轨卫星，葵花8号静止卫星具有空间分辨率高、观测频率高、观测范围广的应用优势，目前葵花卫星数据已经应用于天气、海洋、灰霾、火点监测等。由于AHI传感器设置了对气溶胶较为敏感的可见光波段，因此在气溶胶反演方面潜力巨大，已有葵花卫星官方产品包括AOD(aerosol optical depth)、AMV(atmospheric motion vector)、CSR(clear sky radiance)、HCAI(high-resolution cloud analysis information)。本文将采用葵花卫星AHI传感器可见光红、蓝波段数据(分辨率分别为0.5和1km)进行气溶胶光学厚度反演，并结合AERONET站点观测数据、MODIS气溶胶产品以及葵花卫星官方气溶胶产品(空间分辨率为5km)进行算法精度、空间分布和空间覆盖率方面的比对验证。
本研究采用的地面监测验证数据来自于国际AERONET观测站。AERONET (AErosol RObotic NETwork)是一个国际性的气溶胶地基观测网站(https://aeronet.gsfc.nasa.gov)，由NASA和PHOTONS在全球范围共建。主要提供3个质量等级的数据：level 1.0为未进行云检测的数据，level 1.5为进行了云去除的数据，level 2.0为进行了云去除和质量检测的数据。本文选用的对比站点为北京地区的Beijing、Beijing-CAMS、Beijing-PKU这3个站点，位置分布见图 1。由于2017年尚未发布level 2.0数据，验证数据所以采用level 1.5作为验证数据。
Fig. 1

	Download: JPG larger image
图 1 北京地区AERONET站点Beijing、Beijing-CAMS、Beijing-PKU分布图 Fig. 1 Location map of AERONET sites of Beijing, Beijing-CAMS, and Beijing-PKU in Beijing area 图 1 北京地区AERONET站点Beijing、Beijing-CAMS、Beijing-PKU分布图 Fig. 1 Location map of AERONET sites of Beijing, Beijing-CAMS, and Beijing-PKU in Beijing area -->

本文同时采用目前流行的极轨卫星MODIS官方气溶胶产品进行比对。MODIS是搭载在Terra和Aqua两颗卫星上的中分辨率成像光谱仪，其数据可以从官网(https://ladsweb.modaps.eosdis.nasa.gov)免费获得。本文采用MOD04暗像元法和深蓝算法结合的气溶胶产品作为交叉验证数据，其空间分辨率为10km。
2 葵花卫星反演AOD算法郭强等^[18]提出利用MODIS可见光波段反演陆地气溶胶光学厚度的算法，该算法无需事先获得地表反射率先验知识或者如经典暗像元法通过近红外通道与可见光通道之间比值关系估算地表反射率，而是利用在较短时间内地物可见光波段反射率比值变化较小的特性进行反演。本文将该算法改进应用于葵花卫星数据，算法原理如下：
根据Vermote等^[19]研究结果，在地表朗伯体、大气水平均一的假设条件下，卫星平台观测到的大气顶反射率(表观反射率)可以表示为

$\begin{aligned} \rho_{\mathrm{TOA}}\left(\mu_{\mathrm{s}}, \mu_{\mathrm{s}}, \varphi\right) &=\\ \rho_{0}\left(\mu_{\mathrm{s}}, \mu_{\mathrm{v}}, \varphi\right) &+\frac{T\left(\mu_{\mathrm{s}}\right) T\left(\mu_{\mathrm{v}}\right) \rho_{\mathrm{s}}\left(\mu_{\mathrm{s}}, \mu_{\mathrm{v}}, \varphi\right)}{\left[1-\rho_{\mathrm{s}}\left(\mu_{\mathrm{s}}, \mu_{\mathrm{v}}, \varphi\right) S\right]} ,\end{aligned} $

(1)

式中：ρ₀是大气路径辐射项等效反射率，ρ_s是地表反射率，T(μ_s)、T(μ_v)分别是入射方向和观测方向的大气透过率，S为大气下界面的半球反射率。在公式中T(μ_s)、T(μ_v)总是同时出现，因此将T(μ_s)T(μ_v)作为一个参数考虑。未知数ρ₀, S, T(μ_s)T(μ_v)与AOD、气溶胶粒子单次散射率、相函数、不对称指数等有关。在实际反演过程中，首先针对试验区域确定一个气溶胶模型(因本文研究区为北京，反演中默认选择大陆型气溶胶模式)，然后用6S、MODTRAN等辐射传输模型模拟在不同观测几何情况下AOD与上述3个参数之间的对应关系，建立查找表。若已知地表反射率和表观反射率，便可通过查找表计算得到AOD。求解式(1)的困难之处就在于很难获得准确的地表反射率值。
根据已有的研究，较短时间内(排除降雨、降雪等改变地表反射率情况)，同一地物的不同观测几何条件，可见光波段地表反射率比值变化较小，可忽略不计^[20]，因此，针对同一地物两次观测，可以近似表达为

$\frac{r_{m}\left(\lambda_{\mathrm{b} 1}\right)}{r_{n}\left(\lambda_{\mathrm{b} 1}\right)}=\frac{r_{m}\left(\lambda_{\mathrm{b} 2}\right)}{r_{n}\left(\lambda_{\mathrm{b} 2}\right)}, $

(2)

式中：r为地表反射率，m、n分别表示两个不同的观测几何条件，b1、b2表示不同的波段。
根据式(2)可以导出

$\frac{r_{m}\left(\lambda_{\mathrm{b} 1}\right)}{r_{m}\left(\lambda_{\mathrm{b} 2}\right)}=\frac{r_{n}\left(\lambda_{\mathrm{b} 1}\right)}{r_{n}\left(\lambda_{\mathrm{b} 2}\right)}=k.$

(3)

通过式(3)可以得出，对于同一地物在两个观测几何条件下，分别对应的两个波段反射率的比值相同。经典暗像元法采用较低反照率地表(暗像元)存在可见光波段(0.47、0.66μm)与短波近红外波段(2.13μm)反射率比值经验关系(即r_0.47μm=0.25r_2.13μm、r_0.66μm=0.5r_2.13μm)。因此，暗像元法的波段比值关系是存在于假设地物为朗伯体及地物光谱反射特性不随时间变化的情况下，是本文式(3)的特例。
因此，对于葵花8号卫星AHI传感器，红、蓝波段地表反射率之间可得如下关系

$\rho_{\mathrm{R}}^{S}=k \rho_{\mathrm{B}}^{S}, $

(4)

式中：ρ_R^S，ρ_B^S分别代表红蓝波段地表反射率。对于不同地物，比值k不同。将红、蓝波段地表反射率比值k代入式(1)中，可得如下方程组：

$\left\{ {\begin{array}{*{20}{l}}{\rho _{\rm{R}}^{{\rm{TOA}}}\left( {{\mu _s}, {\mu _v}, \varphi } \right) = }\\{{\rho _0}\left( {{\mu _s}, {\mu _v}, \varphi } \right) + \frac{{T\left( {{\mu _s}} \right)T\left( {{\mu _v}} \right)\rho _{\rm{R}}^S\left( {{\mu _s}, {\mu _v}, \varphi } \right)}}{{\left[ {1 - \rho _{\rm{R}}^S\left( {{\mu _s}, {\mu _v}, \varphi } \right)S} \right]}}, }\\{\rho _{\rm{B}}^{{\rm{TOA}}}\left( {{\mu _s}, {\mu _v}, \varphi } \right) = }\\{{\rho _0}\left( {{\mu _s}, {\mu _v}, \varphi } \right) + \frac{{T\left( {{\mu _s}} \right)T\left( {{\mu _v}} \right)\rho _{\rm{B}}^S\left( {{\mu _s}, {\mu _v}, \varphi } \right)}}{{\left[ {1 - \rho _{\rm{B}}^s\left( {{\mu _s}, {\mu _v}, \varphi } \right)S} \right]}}, }\\{\rho _{\rm{R}}^s = k\rho _{\rm{B}}^s.}\end{array}} \right.$

(5)

式中：未知数为红波段地表反射率ρ_R^S、蓝波段地表反射率ρ_B^S和AOD。因此只要确定比值k就能求解气溶胶光学厚度AOD。
目前葵花卫星尚无可靠的官方地表反射率产品发布，无法直接得到地表反射率比值k。MODIS官方地表反射率产品(MOD09)是目前国际上较为流行认可采用的产品，已有****采用MOD09用于估算其他卫星的地表反射率^[21-24]。
本文参考已有研究成果^[23-24]，采用奇异值分解(SVD)方法实现应用MODIS产品(MOD09)估算葵花卫星数据地表反射率，进而计算获得葵花卫星可见光波段地表反射率比值。Sayer等^[23]利用ASTER、USGS地物光谱库，分别模拟不同地物在MODIS以及AATSR观测时的光谱反射率，利用基于奇异值分解(SVD)的方法，模拟建立MODIS和AATSR两种传感器的不同光谱响应之间的转换关系，从而实现利用MODIS的地表反射率模拟计算AATSR的地表反射率。郭强^[24]同样采用SVD方法应用MODIS地表反射率产品计算HJ-1 A/1B CCD可见光波段地表反射率产品，进而用于气溶胶反演。前人研究结果表明，SVD分解方法可以很好地减少不同传感器之间因为波段响应函数的差异而产生的地表反射率误差，能够用于不同传感器观测的地表反射率之间的转换。
奇异值分解(singular value decomposition, SVD)是一种统计的方法，可以用来提取一组数据的特征。假设存在一个m×n阶矩阵M，那么矩阵M可以表示为

$\boldsymbol{M}=\boldsymbol{U} \boldsymbol{\mathit{\Sigma}} \boldsymbol{V}^{*}.$

(6)

本文利用Herold等^[25]建立的光谱数据库、MODIS光谱响应函数、Himawri-8光谱响应函数三者共同构建矩阵M，选取其中光谱曲线125条，包括绿色植被、枯萎植被、裸土、水体以及屋顶、沥青道路、水泥道路、停车场、铁轨等共26种典型地物，用以代表城市及周边大部分区域的典型地物类型。矩阵每行代表一个波段(前3行为MODIS、后3行为Himawari-8)，每一列代表一种光谱，矩阵M为6×125阶矩阵。对矩阵M进行奇异值分解后得到矩阵U，从矩阵U中选取前n个奇异向量构成矩阵S：

$\boldsymbol{S}=\left[\boldsymbol{U}_{1}, \boldsymbol{U}_{2}, \cdots, \boldsymbol{U}_{n}\right]$

(7)

取矩阵S的前3行构成S的子矩阵S_M，后3行构成子矩阵S_H。对于矩阵M，第i个列向量的前3个元素代表第i种光谱的MODIS反射率，记作向量V_im。同理，后三个元素代表H8反射率，记作向量V_ih，对于向量V_im，构建一个列向量C_i，存在

$\left\{\begin{array}{l}{\boldsymbol{V}_{i \mathrm{m}}=\boldsymbol{S}_{\mathrm{M}} * \boldsymbol{C}_{i}} \\ {\boldsymbol{V}_{i \mathrm{h}}=\boldsymbol{S}_{\mathrm{H}} * \boldsymbol{C}_{i}}\end{array}\right..$

(8)

求解后可以得到模拟MODIS地表反射率V_im与模拟葵花卫星反射率V_ih的转换关系

${\mathit{\boldsymbol{V}}_{i{\rm{h}}}} = {\mathit{\boldsymbol{S}}_{\rm{H}}}*\left( {\mathit{\boldsymbol{S}}_{\rm{M}}^{ - 1}*{\mathit{\boldsymbol{V}}_{i{\rm{m}}}}} \right).$

(9)

那么将MOD09作为$ {\mathit{\boldsymbol{\tilde V}}_{i{\rm{m}}}}$代入式(9)，就可以对葵花卫星地表反射率$ {\mathit{\boldsymbol{\tilde V}}_{i{\rm{h}}}}$求解。模拟拟合的精度与选取的奇异向量的个数n有关。所用奇异向量的个数n越大，奇异向量越能反应矩阵M的特性，因而拟合的精度越高。我们选取可见光3个波段，因而n的最大值为3。
本文采用上述SVD方法，利用MOD09地表反射率产品估算H8可见光波段的地表反射率产品，进而计算获得红蓝波段比值k。实际反演求解式(5)时，采用基于6S辐射传输模型的价值函数迭代求解方法，即在给定气溶胶模式、红蓝波段比值k情况下，对地表反射率和气溶胶光学厚度分别以一定步长间隔取值，代入6S辐射传输模型进行计算，分别得到模拟红、蓝波段卫星表观反射率，建立模拟的卫星表观反射率和观测卫星表观反射率之间的绝对损失函数(loss function)如下

${\rm{L}}\left( {{\rho _{6{\rm{s}}}}, {\rho _{\rm{s}}}} \right) = \left| {\rho _{{\rm{s}}, {\rm{R}}}^{{\rm{TOA}}} - \rho _{6{\rm{s}}, {\rm{R}}}^{{\rm{TOA}}}} \right| + \left| {\rho _{{\rm{s}}, {\rm{B}}}^{{\rm{TOA}}} - \rho _{6{\rm{s}}, {\rm{B}}}^{{\rm{TOA}}}} \right|.$

(10)

式中：L为6 s模拟计算与卫星观测表观反射率之间的绝对损失函数值；ρ_{s, R}^TOA、ρ_{s, B}^TOA分别为卫星观测红、蓝波段表观反射率，ρ_{6 s, R}^TOA、ρ_{6 s, B}^TOA分别为6S模拟计算红、蓝波段表观反射率。通过设置迭代搜索误差收敛条件，求解最优解。如果出现难以满足搜索误差收敛条件的像元，反演计算时进行标记，在后期预处理时利用邻近有效像元进行插值计算，最后即可实现AOD和地表反射率的同时反演计算。
3 AOD反演验证3.1 AOD反演的精度验证为了验证本文提出算法的反演精度，以北京作为研究区域，选取2017年5月中旬的一个雾霾过程开展研究。首先从日本官方网站下载到葵花卫星L1级可见光波段数据，并进行辐射校正、几何投影和云检测等预处理工作，由于云检测不是本文研究重点，本文采取目视解译方法严格剔除云像元。
葵花卫星可以获得每个10min的高频监测数据，为验证算法精度，本文首先针对3个AERONET站点位置上空给出每10min的葵花卫星反演结果，并分别与AERONET观测、葵花卫星官方产品(L2，空间分辨率为5km)、MODIS官方产品(空间分辨率为10km)进行对比验证。空间匹配时，本文算法反演结果和MODIS官方产品、葵花8号官方产品均采用对应AERONET站点所在经纬度位置为中心的3×3网格平均值；时间匹配时，选取靠近卫星观测最近时刻的AERONET站点观测数据用于对比。结果对比如图 2。
Fig. 2

Download: JPG larger image

图 2 AOD反演结果(H8retrieval_AOD)与AERONET(Beijing_AOD、Beijing_CAMS_AOD、Beijing_PKU_AOD代表站点观测值)、葵花官方产品(H8ARP_AOD)及MODIS官方产品(MOD/MYD04)的对比验证 Fig. 2 Comparison of AOD retrieval results(H8retrieval_AOD)with AERONET observation data(Beijing_AOD, Beijing_CAMS_AOD, and Beijing_PKU_AOD), Himawari-8 official products(H8ARP_AOD), and MODIS official products(MOD/MYD04) 图 2 AOD反演结果(H8retrieval_AOD)与AERONET(Beijing_AOD、Beijing_CAMS_AOD、Beijing_PKU_AOD代表站点观测值)、葵花官方产品(H8ARP_AOD)及MODIS官方产品(MOD/MYD04)的对比验证 Fig. 2 Comparison of AOD retrieval results(H8retrieval_AOD)with AERONET observation data(Beijing_AOD, Beijing_CAMS_AOD, and Beijing_PKU_AOD), Himawari-8 official products(H8ARP_AOD), and MODIS official products(MOD/MYD04) -->

从图 2可以看到，总体上，利用本文算法反演得到的AOD值与AERONET站点数据波动趋势较为一致，并且值的大小也较为接近，在16时之后，反演结果与观测值相比误差较大，且均大于观测值。可能因为16时之后太阳天顶角和相对方位角均较大的缘故导致。
统计可得，与AERONET 3站点Beijing、Beijing_CAMS、Beijing_PKU观测结果相比，本文反演结果的绝对误差平均值分别为0.06、0.05、0.03，相对误差平均值分别为12%、10%、7%。葵花官方产品仅在10时至15时之间有值，绝对误差平均值分别为0.12、0.11、0.07，相对误差平均值分别为20%、18%、14%，总体上看，葵花官方产品与AERONET站点数据相比在上午偏差较大，且存在低估现象。TERRA和AQUA双星MODIS官方气溶胶产品上下午各有一个观测值，绝对误差平均值分别为0.27、0.25、0.33，相对误差平均值分别为48%、46%、67%，MODIS官方产品误差最大，且存在明显高估现象。比较而言，不难看出，本文提出的算法反演精度明显优于葵花官方产品和MODIS官方产品，且在可反演的日时间范围上长于葵花官方产品。
本文反演算法的误差来源可能包括以下几个方面；1)经验设置的气溶胶类型参数，本文采取的是大陆型气溶胶模式，而对于北京区域并不一定符合实际状况；2)地表反射率比值的估算精度，本文应用MOD09地表反射率产品模拟计算获得地表反射率比值，并非真实的葵花地表反射率产品计算的比值；3)目视解译方法进行的云检测可能遗漏薄云像元；4)优化求解、迭代求解时精度设置等。本文的反演案例中，AOD均小于1，对于更为严重的雾霾过程，如AOD1时，需要进一步研究，因6S模型在能见度低于5km时的辐射传输计算可靠性降低，理论上分析，本文算法对于大气能见度低于5km时的AOD反演误差会进一步增大，因此较适用于能见度大于5km时的雾霾状况。
上述卫星数据反演结果与地面观测站点对比验证的误差，除受反演结果本身精度影响外，还与对比验证时的时空匹配精度有关。空间上，相对于10km分辨率的MODIS产品和5km的葵花产品，本文反演的1km空间分辨率产品更容易与地面站点匹配，另外对比时采用的是反演结果的3×3格网均值，因此在对比时空间匹配更具优势。时间上，相比于葵花卫星10min一次观测频率，MODIS只有上下午2次过境，与地面观测的时间匹配更难吻合，时空匹配精度低将进一步增大验证误差。
3.2 AOD反演的空间分布验证为进一步分析验证本文提出算法反演的空间分布的准确性和优势，本文将葵花卫星的反演结果与NASA官网获得的TERRA星上MODIS气溶胶产品(MOD04)进行对比，为了对比时相上的一致性，尽量选择MODIS过境时刻接近的葵花卫星反演结果。对比结果见图 3。
Fig. 3

	Download: JPG larger image
(a)5月17日上午10时55分MOD04气溶胶产品图像；(b)5月17日上午11时葵花卫星反演结果图像; (c)5月18日上午11时44分MOD04气溶胶产品图像; (d)5月18日上午11时50分葵花卫星反演结果图像 图 3 葵花卫星反演结果图像与MOD04气溶胶产品图像对比 Fig. 3 Comparison between retrieval AOD image and MODIS official AOD products image 图 3 葵花卫星反演结果图像与MOD04气溶胶产品图像对比 Fig. 3 Comparison between retrieval AOD image and MODIS official AOD products image -->

图 3中，颜色由蓝到红表示气溶胶光学厚度从低到高。由两颗卫星的AOD反演结果对比可以看出，总体上气溶胶光学厚度的空间分布较为一致，整体呈现东南高，西北低的状况。在东南部城市区域，MOD04值明显比本文算法的反演结果要高，这与前述图 2中MOD04明显高于本文算法反演结果和AERONET站点观测值的对比结果是一致的。另外，由于本文的反演结果空间分辨率(1km)明显优于MODIS官方产品(10km)，可以看出葵花卫星反演结果较MODIS产品空间分布细节清晰。另外也可以看出，MOD04产品存在较多的无值空白区域，这是由于MODIS产品采用的暗像元法局限性导致，而葵花卫星反演算法不受亮地表限制，反演结果有效区域明显大于MOD04。因此，在时间分辨率、空间分辨率以及应用范围的角度，葵花卫星反演结果均优于MODIS官方应用产品。
进一步地，将本文算法反演结果与葵花卫星官方产品进行对比。葵花官方气溶胶产品反演也是采用暗像元法，产品目前仍处于实验阶段。目前官方提供的两种质量级别的气溶胶产品，即L2(10min级)和L3(1h级)两级数据，空间分辨率均为5km，其中L3级数据是L2级数据的小时平均值。由于算法受太阳天顶角限制较大，在中国大陆地区每日10时后才有反演产品，另外除云覆盖影响外，由于算法限制，葵花官方气溶胶产品在中国覆盖率极低，大部分时间只有中国东部沿海地区有数据，除云覆盖影响外，也许与暗像元法可应用的区域局限性有关。如图 4所示，两个级别产品在中国均有较大范围的空白无效区域，其中L3产品是L2数据经过后续筛选处理获得，精度质量控制高于L2级数据，但时间分辨率由10min级降为1h级，有效监测空间范围进一步缩小，无法满足中国大范围区域监测应用的需求。
Fig. 4

	Download: JPG larger image
图 4 2017年5月17日11时葵花卫星官方气溶胶产品(ARP)L2级(a)和L3级(b) Fig. 4 Himawari-8 satellite official aerosol products (ARP) of level 2 (a) and level 3 (b) at 11 a.m. on May 17 2017 图 4 2017年5月17日11时葵花卫星官方气溶胶产品(ARP)L2级(a)和L3级(b) Fig. 4 Himawari-8 satellite official aerosol products (ARP) of level 2 (a) and level 3 (b) at 11 a.m. on May 17 2017 -->

为了兼顾比较L2级(更大空间覆盖范围)和L3级(更好的精度质量控制)数据的优势，本文同时选取葵花官方产品每日11时至14时整点时刻的L2和L3级数据进行对比验证(见图 5)，反演结果图像时间分别为：11、12、13、14整点时刻。
Fig. 5

	Download: JPG larger image
图 5 葵花卫星反演结果图像与葵花卫星官方气溶胶产品图像对比 Fig. 5 Comparison between retrieval AOD image and Himawari-8 official AOD products image 图 5 葵花卫星反演结果图像与葵花卫星官方气溶胶产品图像对比 Fig. 5 Comparison between retrieval AOD image and Himawari-8 official AOD products image -->

总体上三者AOD的空间分布趋势较为一致，2个日期均在北京东南部出现AOD高值区。比较而言，在2个日期里，L3级产品图像均在西北部出现大范围空白无效区域，这是因为L3产品为L2数据后续筛选处理的结果，说明葵花官方产品在优化产品精度时牺牲了产品空间覆盖率；L2级产品图像上零星分布无效空白像元，在整幅图像上呈明显斑点状，但空间覆盖范围明显优于L3级数据。综合空间覆盖范围和空间分布细节两个方面，本文算法反演结果均为最优。可以看出，在算法区域适应性方面，本文算法优于葵花官方产品采用的暗像元法，同时，由于本文算法反演结果空间分辨率(1km)明显高于葵花官方产品(5km)，能够更好满足较高空间分辨率AOD产品的需要。
4 雾霾过程监测雾霾的主要成因是对流层气溶胶，常用表征指标是PM_2.5浓度或空气质量指数AQI，而大气气溶胶属性常用表征指标是气溶胶光学厚度AOD，其定义为从大气底面到大气顶层的垂直气柱内气溶胶消光系数的积分。由于PM_2.5浓度或AQI常常采用的是近地面测量数据，而卫星遥感AOD是整层大气气溶胶消光系数的积分，因此卫星遥感监测雾霾常常采用AOD作为间接大气污染程度指标或应用AOD估算近地面PM_2.5浓度或AQI，诸多研究表明气溶胶光学厚度与大气污染程度呈正相关性，可以定性反映地面大气污染状况^[25-26]。如前所述，本文对北京区域2017年5月17日、18日进行小时级气溶胶AOD反演，反演结果如图 6所示。
Fig. 6

	Download: JPG larger image
图 6 2017年5月17日和18日8时—17时小时级时间序列气溶胶AOD反演结果图像 Fig. 6 Hourly time series AOD retrieval images from 8 a.m. to 17 p.m. on May 17 and 18, 2017 图 6 2017年5月17日和18日8时—17时小时级时间序列气溶胶AOD反演结果图像 Fig. 6 Hourly time series AOD retrieval images from 8 a.m. to 17 p.m. on May 17 and 18, 2017 -->

本文获得35个北京市空气质量监测站点2017年5月17日8时至17时及5月18日8时至16时(17时数据缺失)的PM_2.5监测数据。依据站点的经纬度地理位置以及监测数据采集时间，分别与本文的AOD反演结果图像进行时空匹配，提取出站点位置处的PM_2.5与AOD的数据对，进行相关性分析，如图 7所示。
Fig. 7

	Download: JPG larger image
图 7 北京35个站点空气质量监测站PM2.5数据与AOD相关性分析 Fig. 7 Correlation analysis between PM2.5 data and AOD at 35 air quality monitoring stations in Beijing 图 7 北京35个站点空气质量监测站PM2.5数据与AOD相关性分析 Fig. 7 Correlation analysis between PM2.5 data and AOD at 35 air quality monitoring stations in Beijing -->

本文反演的AOD与地面监测的PM_2.5之间直接对比的皮尔森相关系数为0.502，满足统计学上99%置信度的要求，这说明PM_2.5与AOD之间存在一定程度的线性正相关，但是相关性不是很强。前人的研究结果表明如果进一步对反演的AOD进行湿度订正与垂直订正，有望进一步提高AOD和PM_2.5的相关性^[26-27]。即便如此，亦可进一步说明气溶胶光学厚度与大气污染程度呈一定的正相关性，可以用来定性反映地面大气污染状况，作为城市地面颗粒物浓度监测的一个有效补充手段。
为验证本文算法对雾霾过程监测的可行性，以AOD为雾霾程度监测指标，分析北京区域此次雾霾过程的时空动态变化。从这两日反演AOD图像中，可以清楚地看到，总体上北京区域的雾霾分布呈现东南污染重，西北污染相对轻的趋势。总体上看，5月18日雾霾污染相比于5月17日更为严重。经查询中国空气质量在线监测分析平台，历史数据记录显示5月17日北京区域为中度污染，空气质量指数AQI日均值为184，5月18日为重度污染，AQI日均值为201。说明本文反演的AOD的日变化趋势与官方大气环境监测数据结果基本一致。
在时间序列上，本文采用的葵花卫星数据反演结果具有明显的时间分辨率优势，从5月17日8时至5月18日17时，每隔一小时的反演结果显示AOD空间分布一直在变化，尤其18日相对17日而言，气溶胶空间分布变化更快，说明雾霾时空分布的高度动态变化性特征。时间序列上，可以看出，5月17日8时起北京城区开始轻度雾霾，到下午17时，城中心雾霾加重，且在城南邻近河北区域出现雾霾严重区域。5月18日上午8时起北京城区域再次出现轻度雾霾，且从上午11时起，雾霾开始明显加重，至下午13点达到雾霾最重，然后至下午17时，雾霾开始减弱。另外也可以看出，随着时间变化，雾霾污染中心具有逐渐向北京城东北方向扩散移动的趋势。
从这两日反演结果时间序列图像中，可以清楚地看到雾霾污染中心的移动方向以及加重或减轻的动态变化状况，证明雾霾动态变化强，具有输送、扩散等特点。在实际应用中，使用传统的极轨卫星数据如MODIS等获得的气溶胶数据已无法满足对雾霾动态高频次监测的要求。目前雾霾监测预报主要使用地面监测站点，空间上分布稀疏，难以进行大范围监测和预报。雾霾的高空间分辨率、高时间分辨率以及大范围同时监测的需求，使得静止气象卫星有望成为未来雾霾监测较好的选择。本文研究表明，利用葵花8号卫星高时空分辨率可以进行每1h甚至每10min级的反演，能够动态描绘出污染物的扩散运动变化轨迹，如能结合当时的气象状况，有望成为研究污染物的溯源和扩散路径追踪的有效手段。
5 结论本文提出一种利用葵花8号卫星可见光波段反演高时空分辨率气溶胶光学厚度的算法，并通过与国际AERONET 3个站点Beijing、Beijing-CAMS、Beijing-PKU站观测数据，葵花官方气溶胶产品和MODIS官方气溶胶产品对比分析，验证算法的可行性。结果表明，本文提出算法的反演结果，无论在反演精度、有效反演空间范围、有效反演时间范围、空间分辨率等方面均为最优。该算法不受地表反射类型的限制，可应用于亮地表地区，突破了暗像元法的局限性。相比于MODIS气溶胶数据和葵花官方气溶胶数据，该算法的反演结果具有高覆盖度，高精度、高频次、高分辨率的特点，应用本文算法可以获得每10min级、空间分辨率为1km的气溶胶光学厚度日动态变化反演结果，更有利于对雾霾动态变化的监测。最后通过对北京区域2017年5月17日、18日小时级气溶胶AOD反演结果的时空分析，探讨北京区域此次雾霾过程的时空动态变化，从这两日反演结果时间序列图像中，可以清楚地看到雾霾污染中心的移动方向以及加重或减轻的动态变化状况，说明本文提出的应用葵花卫星监测雾霾过程的方法具有较好应用潜力。存在的问题及未来工作包括：1)云像元剔除由现在的目视解译改成利用自动云检测算法实现动态云掩膜；2)目前仅依据经验采用大陆型气溶胶模式，未来可根据实际情况采用自定义气溶胶模式，进一步提高AOD反演精度；3)本文的反演案例中，AOD均小于1，对于更为严重的雾霾过程，如AOD1时，需要未来进一步研究，理论上分析，本文算法较适用于能见度大于5km时的雾霾状况；4)目前计算地表反射率比值方法是基于MOD09产品的奇异值分解(SVD)模拟计算，未来将探索比较其他地表反射率比值的计算方法，进一步提高精度；5)针对雾霾的监测，本文直接应用葵花卫星反演的AOD对雾霾进行定性分析，未做近地面AOD订正处理，未来将开展反演AOD的湿度订正和垂直订正等后处理，获得近地面AOD后进一步估算颗粒物PM_2.5、PM₁₀浓度等研究。
感谢NASA官方网站提供的MODIS数据产品、日本气象厅官方网站提供的葵花8号卫星数据及国际AERONET站北京Beijing、Beijing-CAMS、Beijing-PKU 3站的数据支持。
参考文献

[1]	王桥, 厉青, 高健, 等. PM2.5卫星遥感技术及其应用[M]. 北京: 中国环境出版社, 2017.
[2]	戴羊羊, 李成范, 周时强, 等. 基于遥感的上海地区雾霾监测研究[J]. 测绘工程, 2015(12): 29-32. Doi:10.3969/j.issn.1006-7949.2015.12.007
[3]	任通, 高玲, 李成才, 等. 利用风云2C静止卫星可见光资料反演气溶胶光学厚度[J]. 北京大学学报(自然科学版), 2011, 47(4): 636-646.
[4]	毛节泰, 刘莉. GMS5卫星遥感气溶胶光厚度的试验研究[J]. 气象学报, 2001, 59(3): 352-359. Doi:10.3321/j.issn:0577-6619.2001.03.010
[5]	张军华, 斯召俊, 毛节泰, 等. GMS卫星遥感中国地区气溶胶光学厚度[J]. 大气科学, 2003, 27(1): 23-35. Doi:10.3878/j.issn.1006-9895.2003.01.03
[6]	高玲, 任通, 李成才, 等. 利用静止卫星MTSAT反演大气气溶胶光学厚度[J]. 气象学报, 2012, 70(3): 598-608.
[7]	陈健, 周杰, 李雅雯. 基于静止卫星GOCI数据的陆地上空气溶胶光学厚度遥感反演[J]. 遥感技术与应用, 2017, 32(6): 1 040-1 047.
[8]	Moulin C, Dulac F, Lambert C E, et al. Long-term daily monitoring of Saharan dust load over ocean using Meteosat ISCCP-B2 data:2. Accuracy of the method and validation using Sun photometer measurements[J]. Journal of Geophysical Research Atmospheres, 1997, 102(D14): 16 959-16 969. Doi:10.1029/96JD02598
[9]	Wang J, Christopher S A, Reid J S, et al. GOES 8 retrieval of dust aerosol optical thickness over the Atlantic Ocean during PRIDE[J]. Journal of Geophysical Research Atmospheres, 2003, 108(D19): 1-53.
[10]	Prados A I, Kondragunta S, Ciren P, et al. GOES aerosol/smoke product (GASP) over North America:comparisons to AERONET and MODIS observations[J]. Journal of Geophysical Research Atmospheres, 2007, 112(D15): 1-15.
[11]	Zou X, Zhuge X, Weng F. Characterization of bias of advanced Himawari imager infrared observations from NWP background simulations using CRTM and RTTOV[J]. Journal of Atmospheric & Oceanic Technology, 2016, 33(12): 2 553-2 567.
[12]	Kurihara Y, Murakami H, Kachi M. Sea surface temperature from the new Japanese geostationary meteorological Himawari-8 satellite[J]. Geophysical Research Letters, 2016, 43(3): 1 234-1 240. Doi:10.1002/2015GL067159
[13]	Shang H Z, Chen L F, Letu H S, et al. Development of a daytime cloud and haze detection algorithm for Himawari-8 satellite measurements over central and eastern China[J]. Journal of Geophysical Research Atmospheres, 2017, 122(6): 3 528-3 543. Doi:10.1002/2016JD025659
[14]	Daisaku U. Aerosol Optical Depth product derived from Himawari-8 data for Asian dust monitoring[R]. Meteorological Satellite Center Technical Note, 2016, 61: 59-63.
[15]	Sekiyama T T, Yumimoto K, Tanaka T Y, et al. Data assimilation of Himawari-8 aerosol observations:Asian dust forecast in June 2015[J]. SOLA-Scientific Online Letters on the Atmosphere, 2016, 12: 86-90.
[16]	葛邦宇, 杨磊库, 陈兴峰, 等. 暗目标法的Himawari-8静止卫星数据气溶胶反演[J]. 遥感学报, 2018, 22(1): 38-50.
[17]	Wang W, Mao F Y, Du L, et al. Deriving hourly PM2.5 concentrations from Himawari-8 AODs over Beijing-Tianjin-Hebei in China[J]. Remote Sensing, 2017, 9(8): 858. Doi:10.3390/rs9080858
[18]	郭强, 唐家奎, 何文通, 等. 利用MODIS可见光波段反演陆地气溶胶光学厚度[J]. 地理与地理信息科学, 2015, 31(2): 38-43. Doi:10.3969/j.issn.1672-0504.2015.02.009
[19]	Vermote E F, Tanre D, Deuze J L, et al. 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. Doi:10.1109/36.581987
[20]	Flowerdew R J, Haigh J D. An approximation to improve accuracy in the derivation of surface reflectances from multi-look satellite radiometers[J]. Geophysical Research Letters, 2013, 22(13): 1 693-1 696.
[21]	吕春光, 田庆久, 王磊, 等. 利用Landsat-8 OLI反演大气气溶胶的可见光谱段地表反射率估算[J]. 遥感信息, 2015(1): 43-50. Doi:10.3969/j.issn.1000-3177.2015.01.008
[22]	麻盛芳.MODIS地表反射率产品支持的HJ-1CCD数据高精度大气纠正[D].青岛: 山东科技大学, 2013. http://d.wanfangdata.com.cn/Thesis/Y2434319
[23]	Sayer A M, Thomas G E, Grainger R G, et al. Use of MODIS-derived surface reflectance data in the ORAC-AATSR aerosol retrieval algorithm:impact of differences between sensor spectral response functions[J]. Remote Sensing of Environment, 2012, 116(2): 177-188.
[24]	郭强.京津冀地区气溶胶光学厚度及空气质量指数遥感反演研究[D].北京: 中国科学院大学, 2014. http://d.g.wanfangdata.com.cn/Thesis_Y2954934.aspx
[25]	Herold M, Gardner M E, Roberts D. Spectral resolution requirements for mapping urban areas[J]. IEEE Transactions on Geoscience & Remote Sensing, 2003, 41(9): 1 907-1 919.
[26]	李成才, 毛节泰, 刘启汉, 等. 利用MODIS光学厚度遥感产品研究北京及周边地区的大气污染[J]. 大气科学, 2003, 27(5): 869-880. Doi:10.3878/j.issn.1006-9895.2003.05.08
[27]	王静, 杨复沫, 王鼎益, 等. 北京市MODIS气溶胶光学厚度和PM2.5质量浓度的特征及其相关性[J]. 中国科学院研究生院学报, 2010, 27(1): 10-16.