出版日期: 2018-01-25
点击次数：
下载次数：
DOI: 10.11834/jrs.20187033
2018 | Volumn22 | Number 1
上一篇 | 下一篇

暗目标法的Himawari-8静止卫星数据气溶胶反演

葛邦宇^1,2 , 杨磊库¹ , 陈兴峰² , 李正强² , 梅笑冬³ , 刘李⁴

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

2. 中国科学院遥感与数字地球研究所国家环境保护卫星遥感重点实验室，北京 100101

3. 中国科学院大学资源与环境学院，北京 100049

4. 中国资源卫星应用中心，北京 100094

收稿日期: 2017-02-17; 优先数字出版日期: 2018-01-01

基金项目: 国家自然科学基金(编号：41401403，41501399，41601392和41401424)；河南省科技攻关项目(编号：162102310089)；河南理工大学博士基金(编号：B2014-018，B2016-14)

第一作者简介: 葛邦宇(1989—　)，男，硕士研究生，研究方向为气溶胶卫星遥感反演。E-mail：geby2015@163.com

通讯作者简介: 杨磊库(1980—　)，男，副教授，研究方向为气溶胶卫星遥感反演。E-mail：leiku.yang@foxmail.com

中图分类号: P407.4

文献标识码: A

摘要

Himawari-8 (H8) 是由日本气象厅发射的新一代静止气象卫星，可实现10 min/次的高频次对地观测，搭载的AHI (Advanced Himawari Imager) 传感器设置有与MODIS暗目标气溶胶反演算法所需的类似波段。本文参考暗目标算法构建了针对该卫星传感器的陆地气溶胶反演算法：首先，通过基于地基站点观测数据的精确大气校正，统计得到短波红外与可见光波段的地表反射率比值关系，将此作为先验知识用于地—气解耦时的反射率估计；然后，初步假设大陆型气溶胶类型，利用辐射传输模型建立查找表；最后，通过模拟与卫星观测的表观反射率误差最小实现气溶胶光学厚度反演解算。选取2016年5月覆盖京津冀地区的观测数据进行测试，将反演结果与对应时间的MODIS气溶胶光学厚度产品进行对比验证，空间分布趋势一致、相关性较高，相关系数R达到0.852；通过与地基观测网AERONET站点实测数据对比验证，所有站点的相关系数R²均大于0.88，精度较高。利用反演的高时间分辨率产品，分析了京津冀地区的大气空间分布和日变化情况，结果表明：采用暗目标法对H8静止卫星陆地气溶胶光学厚度反演具有一定的潜力和可行性，能反映气溶胶的高时间变化信息，有望成为大气环境污染变化监测新的重要手段。

关键词

葵花8 (Himawari-8), 暗目标法, 地表反射率, 气溶胶光学厚度, 京津冀

Study on aerosol optical depth retrieval over land from Himawari-8 data based on dark target method

GE Bangyu^1,2 , YANG Leiku¹ , CHEN Xingfeng² , LI Zhengqiang² , MEI Xiaodong³ , LIU Li⁴

1.School of Surveying and Land Information Engineering, Henan Polytechnic University, Jiaozuo 454003, China

2.State Environmental Protection Key Laboratory of Satellite Remote Sensing, Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing 100101, China

3.College of Resources and Environment, University of Chinese Academy of Sciences, Beijing 100049, China

4.China Centre for Resources Satellite Data and Application, Beijing 100094, China

Abstract

Himawari-8 (H8), as a new generation of geostationary meteorological satellites that observes full-disk images (images of the Earth as seen from the satellite) per 10 min, was launched by Japan Meteorological Agency to investigate aerosol characteristics. The Advanced Himawari Imager (AHI) onboard Himawari-8 has similar spectral bands with Moderate Resolution Imaging Spectroradiometer (MODIS). This study applies the Dark Target (DT) method for Aerosol Optical Depth (AOD) retrieval from AHI data. The atmospheric effect is established for the AHI data over AErosol RObotic NETwork sites. The ratio of surface reflectance between the shortwave infrared and visible bands is then obtained. This ratio serves as a priori knowledge for the surface reflectance estimation in the atmosphere-surface coupling model. Assuming that the aerosol type is continental, we build a look-up table through the radiative transfer model. With the retrieval algorithm, the AOD is determined in the case of a minimum difference between simulated apparent reflectance and satellite observations. The algorithm was used to retrieve AOD over the Beijing-Tianjing-Hebei area of China in May 2016. H8 AOD products were compared with MODIS products, and the results revealed a good spatial coincidence, with the correlation coefficient R being 0.852. The H8 AOD products were validated with AERONET observations, and they showed good linear relations, with the correlation coefficient R² being better than 0.88. The high temporal resolution products were used to analyze aerosol spatial distribution and diurnal variation in the Beijing-Tianjin-Hebei region. Results show that H8 AOD retrieval based on the DT method has certain feasibility and potential. AOD products can express the high temporal variation of aerosol and are thus potentially useful in atmospheric environmental pollution monitoring.

Key words

Himawari-8, dark target method, surface reflectance, aerosol optical depth, Beijing-Tianjin-Hebei area

1 引　言

大气气溶胶对地气系统、辐射平衡、大气环境有显著影响。根据IPCC第5次报告(IPCC，2014)，气溶胶引起的辐射强迫不确定性达100%，被认为是影响气候模型不确定性的最大因素之一(Hansen 等，2011)；同时高浓度的气溶胶造成的雾霾污染危害人类健康。因此对气溶胶性质、时空分布、变化趋势的监测尤为重要，其中卫星遥感监测目前被认为是一种重要手段。

气溶胶的卫星遥感反演研究始于19世纪70年代中期，根据传感器特性不同发展出不同的反演方法，如AVHRR(Advanced Very High Resolution Radiometer)的单通道法(Ignatov和Stowe，2000)、POLDER(Polarization and directionality of the Earth’s reflectances)的偏振算法(Herman 等，1997)、MISR(Multi-angle Imaging SpectroRadiometer)的多角度算法(Liu 等，2004)、MODIS(Moderate Resolution Imaging Spectroradiometer)的暗目标法(Kaufman 等，1997；Chu 等，2003；Levy 等，2007, 2013)。MODIS气溶胶产品应用最为广泛，海洋上空产品精度可达到±(0.03+5%)，陆地上空可达到±(0.05+15%)(Levy 等，2013)。上述传感器均搭载在极轨卫星上，每天1—2次实现全球覆盖观测，然而气溶胶的生命周期短、时空变化性强，为研究其动态变化、传输特性等，极轨卫星传感器时间分辨率已不能满足要求。

静止卫星数据具有高时间分辨率的特点，可用来研究气溶胶光学厚度AOD(Aerosol Optical Depth)的日变化，甚至小时变化(张军华等，2003)。毛节泰等人(2001)、张军华等人(2003)利用日本的静止气象卫星GMS5可见光通道反演气溶胶光学厚度；高玲等人(2012)利用日本静止卫星MTSAT采用最小卫星反射率反演地表反射率的方法反演AOD；国际上利用静止卫星进行AOD的反演也有大量成果，如Moulin等人(1997)利用Meteosat静止气象卫星可见光通道反演了1983年—1994年期间由撒哈拉沙漠输送到大西洋和地中海上空的沙尘量；Wang等人(2003)、Prados等人(2007)利用GOES数据反演气溶胶；Guerrieri等人(2007)利用MSG反演地中海区域的气溶胶。由于上述静止卫星传感器缺少对气溶胶敏感的蓝光波段设置，AOD反演精度不高，多处于试验阶段，且很少有AOD产品发布。

日本气象厅于2014年7月发射了新一代静止气象卫星H8(Himawari-8)(Bessho 等，2016)，可实现全区域拍摄时间分辨率为10 min/次的高频次对地观测，引起了广泛关注。目前针对该卫星传感器逐步开展了如天气预报、海面温度、灰霾检测等相关研究(Kurihara 等，2016；Zou 等，2016；Shang 等，2017)。H8搭载的AHI(Advanced Himawari Imager)传感器波段设置包含对气溶胶敏感的蓝光通道，在陆地气溶胶反演方面具有很大的潜力，然而目前针对该卫星传感器的算法多处于实验研究阶段。Uesawa(2016)(http://www.data.jam.go.jp/mscweb/technotes/msctechrep61-6:pdf[2016-12-20])利用H8数据采用双通道算法(可见光0.64 μm和近红外0.86 μm波段)进行了气溶胶反演试验，该算法基于沙尘反演算法，未利用对气溶胶敏感的蓝光波段，且文中没有明确地表波段关系。Yumimoto等人(2016)将H8反演的气溶胶光学厚度用于全球气溶胶预测模型Kalman滤波系统，并预测东亚的一次沙尘暴(Sekiyama 等，2016)。AHI传感器设置了短波红外2.3 μm通道，该波段受大气影响较小，与MODIS波段设置类似，使得采用类似MODIS暗目标算法进行H8气溶胶反演成为可能。本文通过结合地基气溶胶观测数据对站点上空的H8数据进行精确大气校正，初步统计出H8的地表反射率波段关系用于气溶胶反演时的地表估计，并对配套的云、水、冰雪掩膜、查找表制备方法等进行了改进，从而构建了H8陆地AOD反演算法。

2 数据与方法

2.1 数据

2.1.1 Himawari-8数据

H8静止卫星搭载的AHI传感器具有16个波段，第3波段星下点空间分辨率0.5 km，1、2、4波段空间分辨率1 km，5—16波段空间分辨率2 km，如表1。反演方法利用1—6波段，其中云、水、雪检测基于第1(0.46 μm)、2(0.51 μm)、4(0.86 μm)、5(1.6 μm)波段；AOD反演主要选择与MODIS陆地暗目标方法相似的1(0.46 μm)、3(0.64 μm)、6(2.3 μm)波段，已有研究将AHI与NPP VIIRS数据进行交叉对比，发现上述波段二者吻合度较高，表明该传感器定标精度较高(Yu和Wu，2016)。H8数据由日本宇宙航空研究开发机构JAXA(Japan Aerospace Exploration Agency)发布，用户通过网站注册后可以免费获取数据(http://www.eorc.jaxa.jp/ptree/index.html[2017-02-17])。

表 1 AHI波段参数设置
Table 1 The band settings of AHI

下载CSV

波段	中心波长/μm	空间分辨率/km
1	0.46	1
2	0.51	1
3	0.64	0.5
4	0.86	1
5	1.6	2
6	2.3	2
7	3.9	2
8	6.2	2
9	7.0	2
10	7.3	2
11	8.6	2
12	9.6	2
13	10.4	2
14	11.2	2
15	12.3	2
16	13.3	2

2.1.2 MODIS数据

MODIS是搭载在极轨卫星Terra、Auqa上的中分辨率成像光谱仪。目前国际上公开发布的AOD产品中，MODIS的AOD产品具有较高的验证精度，得到广泛认可，常用于检验其他传感器气溶胶产品或新算法。Terra赤道平均过境时间是地方时10:30左右，Aqua赤道平均过境时间是地方时13:30左右。本文将每天两次过境的MODIS气溶胶产品，通过时间匹配用于H8 AOD反演结果的交叉对比分析。

2.1.3 AERONET数据

利用地基站点气溶胶观测数据与卫星AOD反演结果进行对比验证被认为是目前最为精确的一种检验方式。美国国家宇航局(NASA)组建了气溶胶自动监测网络AERONET(AErosol RObotic NETwork)，利用法国CIMEL公司的CE-318型太阳辐射计，采用统一的反演方法，其发布的AOD产品不确定性为0.01—0.02(Eck 等，1999)，为全球气溶胶特性研究提供了大量资料。

本文选择京津冀地区AERONET长期观测站点：Beijing (39.977°N，116.381°E)、Beijing_CAMS (39.933°N，116.317°E)、Beijing_RADI (40.005°N，116.379°E)、XiangHe (39.754°N，116.962°E)、SONET_Xingtai (37.182°N，114.360°E)，其中SONET_Xingtai站点是中国太阳—天空辐射计观测网SONET(Sun-sky radiometer Oberservation NETwork)(Li 等，2015)与AERONET的共享站点。观测数据优先采用经人工筛选质量级别高的Level 2.0产品，没有Level 2.0产品的站点选用Level 1.5产品。

2.2 反演原理

卫星对地观测获得的表观信息耦合了地表和大气两部分信号，为获取气溶胶信息，首先进行地气解耦合。关键技术在于：(1) 地表反射率的准确估计；(2) 气溶胶模式的准确确定(陈良富等，2011)。MODIS暗目标法的基本原理是在没有大气影响的情况下，植被、部分裸地等暗目标地物反射率较低，且在短波红外波段与可见光红、蓝波段具有较好的线性相关性，可利用短波红外波段实现地表参数的估计，通过气溶胶模型假设，利用辐射传输方程建立查找表LUT(Look Up Table)，实现AOD反演。

假设大气水平均一，传感器观测的地表目标为均匀朗伯体，在大气窗口波段不考虑气体吸收的情况下，卫星接收到的大气层顶表观反射率ρ_TOA表达为(Kaufman 等，1997)

${\rho _{{\rm{TOA}}}}\left({{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \phi } \right) \!=\! {\rho _0}\left({{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \phi } \right) + \frac{{T\left({{\mu _{\rm{s}}}} \right) \cdot T\left({{\mu _{\rm{v}}}} \right){\rho _{\rm{s}}}\left({{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \phi } \right)}}{{\left[ {1 \!-\! {\rho _{\rm{s}}}\left({{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \phi } \right)S} \right]}}$

(1)

式中，μ_s=cosθ_s，μ_v=cosθ_v，θ_s与θ_v分别为太阳天顶角与观测天顶角； $ \phi $ 为相对方位角；ρ₀是大气的路径辐射项，是大气中分子和气溶胶散射的反射率；S为大气下界的半球反射率；T为大气透过率；ρ_s为地表反射率，当ρ_s很小时，ρ_TOA主要取决于大气贡献项，当ρ_s很大时，地表成为主要贡献。因此，ρ_TOA既是气溶胶光学厚度的函数，又是下垫面反射率的函数，如果已知下垫面反射率，则卫星观测的表观反射率扣除下垫面贡献即可得到气溶胶AOD值。

2.3 反演方法

暗目标法反演气溶胶光学厚度的原理，是通过统计波段间的经验关系作为先验知识，利用受大气影响小的短波红外波段(2.3 μm)估计红(0.64 μm)、蓝(0.47 μm)波段的地表反射率，结合假设的气溶胶类型利用辐射传输模型建立查找表LUT，利用辐射传输模型式(1)模拟卫星观测的表观反射率TOA(Top Of Atmosphere)，将模拟值与卫星真实观测的表观反射率进行匹配，按照一定的价值函数约束反演AOD。

H8反演流程如图1所示，主要分为3部分：(1) 数据的准备与输入；(2) 数据预处理：对H8原始数据进行格式投影转换和裁剪，对裁剪区域进行云、水、冰雪像元的识别，计算研究区域太阳天顶角、方位角、观测天顶角、方位角、水陆掩膜；(3) AOD反演：在预处理图像中选择暗像元，利用短波红外求出红、蓝波段在暗像元区域的地表反射率，根据6S辐射传输模型建立LUT，找出在相同几何条件以及地表反射率的情况下，使红、蓝波段表观反射率和查找表中模拟表观反射率之间平方差和最小的0.55μm AOD作为反演结果。

图 1 H8卫星数据反AOD演流程图

Fig. 1 Flow chart of H8 satellite data retrieval AOD

短波红外、红、蓝波段之间的反射率比值关系是实现地气解耦合的关键因素之一，也是本文研究要解决的关键问题。

2.3.1 利用6S辐射传输模型建立查找表

式(1)概括性地描述了卫星观测的表观反射率、地物地表反射率、大气气溶胶的参数之间的关系，然而实际的地气系统包含了更多地气特性参数，一个单独参数的反演往往无法通过辐射传输方程求得解析解。在气溶胶遥感反演中，利用建立LUT反演气溶胶获得最优解是一种通用且高效的方法。利用6S辐射传输模型(Kotchenova 等，2008)，根据H8的波段设置建立查找表，输入参数包括：太阳天顶角分别为0°，6°，12°，24°，36°，48°，54°，60°，66°；观测天顶角范围0°—72°，间隔6°；相对方位角范围0°—180°，间隔12°；作为算法测试，气溶胶模式选择大陆型气溶胶(王中挺等，2015)；大气模式根据观测时间和经纬度确定；550 nm气溶胶光学厚度分别为0，0.25，0.5，1，2，3，5；地表反射率分别为0，0.01，0.05，0.08，0.1，0.12，0.15，0.18，0.2，0.23，0.25；H8红、蓝波段光谱响应函数。根据上述参数模拟卫星接收到的表观反射率，提取每次辐射传输计算的相关地气参数建立查找表(王新强等，2003)。

2.3.2 云、水、冰雪掩膜

在云、水、冰雪覆盖区域，无法应用暗目标法反演AOD。根据H8数据云、水、冰雪像元的反射率统计特征，获取合适的阈值检测和去除云、水、雪像元(Remer 等，2005)，具体公式如下：

(1)云像元检测：当满足式(2)或(3)，即0.51 μm波段3×3窗口的表观反射率 $ {\rho _{0.51}} $ 标准差大于0.006，或蓝波段的表观反射率 $ {\rho _{0.47}} $ 大于0.25，则判定该窗口内所有像元为云像元。

${\sigma _{0.51}} = \sqrt {\frac{1}{9}\mathop \sum \limits_{i = 1}^9 {{\left( {{\rho _{0.51}} - {\mu _{0.51}}} \right)}^2}} > 0.006$

(2)

${\rho _{0.47}} > 0.25$

(3)

(2)陆地水体检测：当满足式(4)，即植被指数NDVI小于0时，则判定为水体像元。

${\rm{NDVI}} = \displaystyle\frac{{{\rho _{0.86}} - {\rho _{0.64}}}}{{{\rho _{0.86}} + {\rho _{0.64}}}} < 0$

(4)

(3)冰雪检测：当满足式(5)，即冰雪指数NDSI(ρ_0.51和ρ_1.6)大于0.13时，则判定为冰雪像元。

${\rm{NDSI}} = \displaystyle\frac{{{\rho _{0.51}} - {\rho _{1.6}}}}{{{\rho _{0.51}} + {\rho _{1.6}}}} > 0.13$

(5)

2.3.3 地表波段关系

AOD反演的关键是实现地气解耦合，实现地气解耦合的关键之一在于地表反射率的准确估计，地表反射率0.01的误差能导致AOD反演误差达到0.1(Kaufman 等，1997)。由于波段设置、成像方式以及定标的差异，已有的MODIS波段比值关系不能直接用于H8。因此，建立针对H8数据的地表波段关系，通过H8短波红外波段表观反射率准确估计红、蓝波段的地表反射率，是本文利用暗目标法反演AOD的关键问题。

在MODIS暗目标气溶胶反演算法C5版本之前，暗目标区域仅指短波红外2.1 μm波段反射率小于0.1的浓密植被区，且各波段地表反射率服从如下关系

$\left\{ {\begin{array}{*{20}{c}}{\rho _{\rm{R}}^{\rm{s}} = 1/2\rho _{{\rm{SWIR}}}^*}\\[7pt]{\rho _{\rm{B}}^{\rm{s}} = 1/4\rho _{{\rm{SWIR}}}^*}\end{array}} \right.$

(6)

式中， ${\rho _{\rm{R}}^{\rm{s}}}$ 、 ${\rho _{\rm{B}}^{\rm{s}}}$ 分别表示红、蓝波段的地表反射率， ${\rho _{{\rm{SWIR}}}^*}$ 为短波红外波段的表观反射率，即暗目标红蓝波段的地表反射率由短波红外波段的表观反射率按比值关系获得。C5版本算法后，暗目标被扩展到反射率小于0.25的地表，即部分反射率较低的城市和裸地。同时算法改进了地表波段关系，并考虑了地表类型(短波红外植被指数)和方向性(散射角)的影响，红波段的地表反射率通过短波红外的表观反射率估计得到，蓝波段的地表反射率通过红波段估计的地表反射率得到，关系如下

$\left\{ {\begin{array}{*{20}{l}}{\rho _{\rm{R}}^{\rm{s}} = f\left( {\rho _{{\rm{SWIR}}}^*} \right)}\\[7pt]{\rho _{\rm{B}}^{\rm{s}} = g\left( {\rho _{\rm{R}}^{\rm{s}}} \right)}\end{array}} \right.$

(7)

式中， ${\rho _{\rm{R}}^{\rm{s}}}$ 、 ${\rho _{\rm{B}}^{\rm{s}}}$ 分别表示红、蓝波段估计的地表反射率， $\,{\rho _{{\rm{SWIR}}}^*}$ 为短波红外波段的表观反射率，地表波段关系函数f、g请见参考文献(Levy 等，2007)。H8与MODIS波段设置相似，然而由于二者光谱响应函数、成像方式和定标的差异，无法直接使用MODIS暗目标算法中的地表波段关系。

为此，本文结合地基站点的大气观测信息对其上空的H8观测数据进行大气校正，获得较准确的地表波谱反射率，进而统计分析构建适用于H8的地表波段关系。针对京津冀地区，选取干洁大气且分布均匀的图像，以AERONET的Beijing站点、Xianghe站点、Sonet_Xingtai站点为中心进行裁剪得到子图像，对子图像进行大气校正，得到地表反射率，最后对所有子图像的地表反射率统计3个波段之间的关系。具体流程如下：

(1) 暗像元筛选：对经过云、水、冰雪掩膜处理的子图像，根据短波红外波段的表观反射率 $\;{\rho _{{\rm{2}}{\rm{.3}}}^*}$ 按照 $0.01 \leqslant \rho _{2.3}^{\rm{*}} \leqslant 0.25$ 条件筛选暗像元。

(2) 大气校正：输入子图像的几何角度、获取时间、海拔，根据所处地理位置和实际的天气状况选择大气模式，气溶胶模式选择大陆型气溶胶模式，气溶胶值由子图像中心AEROENT站点观测的AOD转换得到，通过6S辐射传输计算大气校正参数，进而得到子图像的地表反射率(张婷媛等，2009)。

(3) 地表关系统计：根据AERONET观测数据，筛选了不同季节清洁大气条件下8100个像元进行大气校正，计算了红(0.64 μm)和蓝(0.47 μm)、短波红外(2.3 μm)波段地表反射率，并对其进行统计回归分析(图2)，得到地表波段关系如下

$\left\{ {\begin{array}{*{20}{l}}{\rho _{0.64}^{\rm{s}} = 0.66\rho _{2.3}^{\rm{s}}}\\[7pt]{\rho _{0.47}^{\rm{s}} = 0.49\rho _{0.64}^{\rm{s}} - 0.005}\end{array}} \right.$

(8)

式中， ${\rho _{0.47}^{\rm{s}}}$ 、 ${\rho _{0.64}^{\rm{s}}}$ 、 ${\rho _{2.3}^{\rm{s}}}$ 分别表示红、蓝、短波红外大气校正后的地表反射率。式(8)可用于H8 AOD反演时可见光波段的地表反射率估计。然而由于中国区AERONET站点较少，地表类型缺乏多样化，H8数据时间积累不充足，目前较难统计出类似MODIS暗目标算法考虑植被指数NDVI和散射角的地表波段关系。

2.3.4 AOD反演

借鉴MODIS暗目标算法，剔除云、水体、冰雪像元，依据短波红外2.3 μm波段的反射率 $0.01 \leqslant \rho _{2.3}^{\rm{*}} \leqslant 0.25$ 的条件，在5×5窗口中选出暗像元，去除最亮的50%和最暗的20%像元后取平均值用于反演计算(Levy 等，2007)，聚合处理后的数据空间分辨率为10 km。

按照式(8)，通过短波红外波段(2.3 μm)估计出红、蓝波段的地表反射率。根据太阳几何、卫星观测几何，以及红蓝波段地表和表观反射率，利用查找表进行查找，如式(9)，当查找出的0.55 μm AOD使蓝、红波段的实际观测的表观反射率( $\;\rho _{0.47}^{\rm{*}}$ 和 $\;\rho _{0.64}^{\rm{*}}$ )和理论的表观反射率( $\;\rho _{0.47}^{\rm{l}}$ 和 $\;\rho _{0.64}^{\rm{l}}$ )平方差和Cost最小时，对应的AOD值即为反演结果。

${\rm{Cost}} = {\left( {\rho _{0.47}^* - \rho _{0.47}^{\rm{l}}} \right)^2} + {\left( {\rho _{0.64}^* - \rho _{0.64}^{\rm{l}}} \right)^2}$

(9)

图 2 H8大气校正地表反射率关系图

Fig. 2 The correlation of H8 atmospheric correction surface reflectance

3 结果与分析

本文以京津冀为研究区域(经纬度范围为113°E—120°E，36°N—43°N)，为确保研究区域植被茂密、云量较少，选取5月份世界时(UTC)2:00—8:00逐小时数据，该时间段中国区域太阳天顶角较小，光照条件好，从而保证足够的反演结果和验证点。另外，选取UTC2:30和UTC5:30两个时刻，用于与MODIS对比验证。

本文试验区选取京津冀地区，经纬度范围为113°E—120°E，36°N—43°N。考虑到植被的茂密度和云量的影响，依据H8的快视图浏览经验地选择5月份数据进行测试，以确保更多的反演结果和验证点。按照UTC选取2:00—8:00的京津冀地区的每小时观测数据，这个时间段对应北京时间10:00—16:00，整个京津冀区域均具有较好的光照条件，可用于日内时间变化分析。另外，选取UTC2:30和UTC5:30每天两个时刻，用于与MODIS对比验证。利用前述方法对测试数据进行反演，最终AOD反演结果的空间分辨率为10 km。

3.1 MODIS结果对比

将H8暗目标方法与MODIS暗目标方法反演的相同空间分辨率(10 km)的0.55 μm AOD产品进行对比分析，图3分别为2016年5月21日UTC 2:30(图3(a))和5:30(图3(b))京津冀地区的H8假彩色图像。图4(a)(c)和(b)(d)分别为2016年5月21日UTC 2:30和5:30的H8和MODIS反演的AOD。可以看出：H8和MODIS反演的AOD结果在空间分布上一致，本文反演结果空间覆盖的有效范围更大，但存在细微差别，两者差异可能来自于卫星观测角度不一致和云检测方法不同。图5为2016年5月21日UTC 5:30两颗卫星相同地点的AOD结果对比，从回归分析结果看出：反演结果总体与MODIS相关性较好，相关系数R达到0.852，表明算法具有可行性和一定潜力；存在一些异常点，可能由于其在云边缘，二者云检测差异造成的；从回归方程来看，斜率小于1，截距大于0，表明目前的算法相比MODIS产品，在AOD低值区偏高、高值区偏低，可能与地表估计和气溶胶类型假设相对较粗有关。

图 3 2016年5月21日UTC 2:30和UTC 5:30京津冀地区H8假彩色图

Fig. 3 H8 false color map over Beijing-Tianjin-Heibei region in UTC 2:30 and UTC 5:30 May 21, 2016

图 4 2016年5月21日UTC2:30和5:30 H8反演的AOD结果与MODIS反演的AOD结果对比

Fig. 4 Comparison between the AOD from H8 and MODIS (UTC 2:30 and 5:30) in May 21, 2016

图 5 2016年5月21日UTC 5:30 H8 AOD反演结果与MODIS AOD数据对比

Fig. 5 AOD product of H8 versus MODIS (UTC 5:30) in May 21, 2016

3.2 地基实测验证

H8卫星反演为0.55 μm波段处的AOD，而AERONET没有直接观测的0.55 μm的AOD产品，一般根据0.44 μm、0.5 μm、0.67 μm采用二次多项式拟合的方法计算得到0.55 μm的AOD(Eck 等，1999)。卫星反演AOD验证多采用时空匹配验证方法(Levy 等，2013)：(1) 时间匹配上，将卫星观测时间前后半小时地基观测数据的平均值作为该时刻测值；(2) 空间匹配上，以站点为中心，将空间半径25 km范围内H8反演的AOD像元进行平均，作为该站点处AOD反演结果。

按照上述验证方法将AOD反演结果与地面站点观测进行时空匹配，验证结果见图6，图6(a)为2016年5月H8反演AOD和5个站点验证的结果，可以看出，AOD反演结果和所有站点验证的相关性较高，相关系数R²达到0.882，标准偏差为0.09，但反演结果相对于AERONET观测具有在AOD低值区偏高、高值区偏低的趋势，这与图5和MODIS的反演对比分析结果相似。图6(a)—(f)分别为2016年5月H8反演AOD和AERONET的Beijing、Beijing_CAMS、Beijing_RADI、SONET_Xingtai、Xianghe站点观测的验证结果，5个站点总体相关性都很高，特别是SONET_Xingtai、Xianghe两个非城市站点，卫星与地基AOD结果一致性较好；而Beijing、Beijing_CAMS、Beijing_RADI等3个城市站点的验证结果中卫星AOD偏高，从而表明图6(a)卫星反演值的总体偏高可能主要来自于城市站点的贡献。暗目标气溶胶反演方法在乡村植被覆盖地表类型反演结果优于人工目标集中的城市较亮地表区域，与目前MODIS的研究结果存在的问题相似(Gupta 等，2016)。

3.3 京津冀高时间分辨率监测应用

京津冀地区位于环渤海区域的西部、华北平原北端，其西侧、北侧背靠太行山和燕山山脉，因污染物人为排放较高和被山脉围拢的特殊地形，是国内大气污染最为严重的区域之一(王跃思等，2014)。

图7为京津冀地区2016年5月UTC5:30 H8(图7(a))和MODIS(图7(b))月平均AOD，从中可以看出H8和MODIS反演的AOD空间分布基本一致，整个地区的AOD分布呈现出北低南高的特征。不仅与京津冀地区特殊的地形有关，还与工业化水平、人口密度有很大关系。图8为5月UTC5:30 H8和MODIS的月平均AOD对比，两者整体相关性相对较好，相关系数R为0.711；部分异常点可能由于云检测、月平均聚合方法的不同引起的；AOD低值区偏高、高值区偏低的区域可能继承于本文算法本身对单景的反演结果。图9为5月份UTC 2:00—8:00(北京时间10:00—16：00)月平均AOD的日变化，从整体来看，整个研究区AOD北低南高；从一天内的时间变化来看，UTC 2:00(北京时间10点)至5:00(北京时间13点)，整个区域AOD值较高，随时间呈升高趋势，可能与交通早高峰和人为排放有关(贺克斌等，2011)；下午UTC 5:00至8:00，整个区域AOD逐渐降低，可能与交通负荷逐渐减少有关。总的来说，应用H8反演高时间分辨率AOD在气溶胶变化研究中具有很大的潜力，能为区域大气环境污染时间变化特性监测提供有力的补充手段。

图 6 H8 2016年5月京津冀地区0.55 μm的AOD卫星产品与AERONET站点观测值对比

Fig. 6 H8 AOD product versus AERONET AOD product in May 2016 for stations of all stations, Beijing, Beijing_RADI, Beijing_CAMS, SONET_Xingtai and Xianghe

图 7 2016年5月UTC 5:30的京津冀地区H8和MODIS的月平均AOD对比

Fig. 7 Comparison between the averaged AOD from H8 and MODIS in 5:30(UTC)May 2016

图 8 2016年5月UTC 5:30的京津冀地区H8和MODIS的月平均AOD对比

Fig. 8 H8 averaged AOD versus MODIS averaged AOD in 5:30(UTC)May, 2016

图 9 2016年5月京津冀地区平均AOD日变化

Fig. 9 Diurnal images of the mean AOD over Beijing-Tianjing-Heibei in May 2016

4 结　论

本文通过对地基气溶胶站点上空的H8数据进行大气校正，初步统计出地表波段关系用于地表反射率估计，并对云、水、冰雪掩膜等进行了改进，构建了基于暗目标算法的H8陆地AOD反演算法。针对京津冀地区，通过对2016年5月UTC 2:00—8:00的H8卫星观测数据进行反演，将反演结果与MODIS气溶胶产品、AERONET地基观测资料进行对比验证，结果具有较好的相关性和验证精度，表明利用暗目标算法对新一代静止卫星陆地气溶胶光学厚度反演具有一定的潜力。

通过结果对比分析，发现可对算法做进一步改进，提高AOD反演的精度：(1) 地表波段关系的进一步优化：目前由于中国区AERONET站点较少且地表类型缺少多样化，H8数据时间积累不充足，本算法仅采用了固定的线性地表波段关系，后续研究通过大量的数据积累可统计出类似MODIS的关于植被指数NDVI和散射角的地表波段关系，从而实现地表反射率的精确估计；(2) 气溶胶模型的确定：目前算法仅采用了大陆型气溶胶，初步表明算法的可行性及应用潜力，下一步可考虑动态气溶胶模型来提高反演精度；(3) 云检测方法的改进：由于H8缺少1.38 μm波段，对卷云污染的检测能力不足，可能是导致本文算法反演结果与MODIS对比分析出现异常点的原因，下一步将考虑引入热红外波段帮助提高卷云的掩膜精度。(4) 亮地表的反演：从反演结果来看，目前本文采用的暗目标法还不能反演胡焕庸线以西的干旱与半干旱的裸地、沙漠等亮目标上空的AOD，而H8数据又缺少深蓝算法需要的0.412 μm波段，进一步研究可以考虑通过合成晴朗天的地表反射率作为先验知识，提高地表估计精度，从而实现亮目标区域的气溶胶反演。

总体来说，本文算法研究结果表明新一代静止卫星用于气溶胶反演的潜力，将为中国同类型的FY-4(已发射，在轨测试中)等卫星的气溶胶反演提供方法参考。该类型静止卫星在大气环境监测中的应用，将为灰霾时间演变、空间输送等信息监测提供有效手段，同时对灰霾预报提供更加有效的数据支撑。

志谢感谢日本宇宙航空研究开发机构提供的H8卫星数据、美国国家宇航局(NASA)哥达德空间飞行中心(GSFC)提供的MODIS气溶胶产品。使用了AERONET的Beijing、Beijing_CAMS、Beijing_RADI、Xianghe、SONET_Xingtai站点的观测数据，作者在此对上述机构和维护AERONET站点的首席科学家表示衷心的感谢！

参考文献(References)

Bessho K, Date K, Hayashi M, Ikeda A, Imai T, Inoue H, Kumagai Y, Miyakawa T, Murata H, Ohno T, Okuyama A, Oyama R, Sasaki Y, Shimazu Y, Shimoji K, Sumida Y, Suzuki M, Taniguchi H, Tsuchiyama H, Uesawa D, Yokota H and Yoshida R. 2016. An introduction to Himawari-8/9—Japan’s new-generation geostationary meteorological satellites. Journal of the Meteorological Society of Japan Serise II, 94 (2): 151–183. [DOI: 10.2151/jmsj.2016-009]

Chen L F, Li S S, Tao J H and Wang Z T. 2011. Research and Application of Aerosol Remote Sensing Quantitative Inversion. Beijing: Science Press (陈良富, 李莘莘, 陶金花, 王中挺. 2011. 气溶胶遥感定量反演研究与应用. 北京: 科学出版社)

Chu D A, Kaufman Y J, Zibordi G, Chern J D, Mao J T, Li C C and Holben B N. 2003. Global monitoring of air pollution over land from the Earth Observing System-Terra Moderate Resolution Imaging Spectroradiometer (MODIS). Journal of Geophysical Research Atmosphere, 108 (D21): 4661 [DOI: 10.1029/2002JD003179]

Eck T F, Holben B N, Reid J S, Dubovik O, Smirnov A, O’Neill N T, Slutsker I and Kinne S. 1999. Wavelength dependence of the optical depth of biomass burning, urban, and desert dust aerosols. Journal of Geophysical Research: Atmospheres, 104 (D24): 31333–31349. [DOI: 10.1029/1999JD900923]

Gao L, Ren T, Li C C, Yang D W, Shi G M and Mao J T. 2012. A retrieval of the atmospheric aerosol optical depth from MTSAT. Acta Meteorologica Sinica, 70 (3): 598–608. [DOI: 10.11676/qxxb2012.049] ( 高玲, 任通, 李成才, 杨东伟, 石光明, 毛节泰. 2012. 利用静止卫星MTSAT反演大气气溶胶光学厚度. 气象学报, 70 (3): 598–608. [DOI: 10.11676/qxxb2012.049] )

Guerrieri L, Corradini S, Pugnaghi S and Santangelo R. 2007. An aerosol optical thickness retrieval algorithm for Meteosat Second Generation (MSG) data over land: applications to the Mediterranean area//Proceedings of the SPIE Volume 6745, Remote Sensing of Clouds and the Atmosphere XII. Florence, Italy: SPIE: 67450D [DOI: 10.1117/12.736648]

Gupta P, Levy R C, Mattoo S, Remer L A and Munchak L A. 2016. A surface reflectance scheme for retrieving aerosol optical depth over urban surfaces in MODIS Dark Target retrieval algorithm. Atmospheric Measurement Techniques, 9 (7): 3293–3308. [DOI: 10.5194/amt-9-3293-2016]

Hansen J, Sato M, Kharecha P and von Schuckmann K. 2011. Earth’s energy imbalance and implications. Atmospheric Chemistry and Physics Discussions, 11 (9): 27031–27105. [DOI: 10.5194/acpd-11-27031-2011]

He K B, Yang F M, Duan F K and Ma Y L. 2011. Atmospheric Particulate Matter and Regional Complex Pollution. Beijing: Science Press (贺克斌, 杨复沫, 段凤魁, 马永亮. 2011. 大气颗粒物与区域复合污染. 北京: 科学出版社)

Herman M, Deuzé J L, Devaux C, Goloub P, Bréon F M and Tanré D. 1997. Remote sensing of aerosols over land surfaces including polarization measurements and application to POLDER measurements. Journal of Geophysical Research: Atmospheres, 102 (D14): 17039–17049. [DOI: 10.1029/96JD02109]

Ignatov A and Stowe L. 2000. Physical basis, premises, and self-consistency checks of aerosol retrievals from TRMM VIRS. Journal of Applied Meteorology, 39 (12): 2259–2277. [DOI: 10.1175/1520-0450(2001)040<2259:PBPASC>2.0.CO;2]

IPCC. 2014. Climate Change 2013: The Physical Science Basis. Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge, United Kingdom and New York: Cambridge University Press

Kaufman Y J, Tanré D, Remer L A, Vermote E F, Chu A and Holben B N. 1997. Operational remote sensing of tropospheric aerosol over land from EOS moderate resolution imaging spectroradiometer. Journal of Geophysical Research: Atmospheres, 102 (D14): 17051–17067. [DOI: 10.1029/96JD03988]

Kotchenova S Y, Vermote E F, Levy R and Lyapustin A. 2008. Radiative transfer codes for atmospheric correction and aerosol retrieval: intercomparison study. Applied Optics, 47 (13): 2215–2226. [DOI: 10.1364/AO.47.002215]

Kurihara Y, Murakami H and Kachi M. 2016. Sea surface temperature from the new Japanese geostationary meteorological Himawari-8 satellite. Geophysical Research Letters, 43 (3): 1234–1240. [DOI: 10.1002/2015GL067159]

Levy R C, Mattoo S, Munchak L A, Remer L A, Sayer A M and Hsu N C. 2013. The Collection 6 MODIS aerosol products over land and ocean. Atmospheric Measurement Techniques Discussions, 6 (11): 159–259. [DOI: 10.5194/amtd-6-159-2013]

Levy R C, Remer L A, Mattoo S, Vermote E F and Kaufman Y J. 2007. Second-generation operational algorithm: Retrieval of aerosol properties over land from inversion of Moderate Resolution Imaging Spectroradiometer spectral reflectance. Journal of Geophysical Research Atmospheres, 112 : D13211 [DOI: 10.1029/2006JD007811]

Li K T, Li Z Q, Li D H, Li W, Blarel L, Goloub P, Benjamin T, Xu H, Xie Y S, Hou W Z, Li L and Chen X F. 2015. Transfer method to calibrate the normalized radiance for a CE318 Sun/sky radiometer. Chinese Optics Letters, 13 (4): 041001 [DOI: 10.3788/COL201513.041001]

Liu Y, Sarnat J A, Coull B A, Koutrakis P and Jacob D J. 2004. Validation of Multiangle Imaging Spectroradiometer (MISR) aerosol optical thickness measurements using Aerosol Robotic Network (AERONET) observations over the contiguous United States. Journal of Geophysical Research: Atmospheres, 109 (D6): D06205 [DOI: 10.1029/2003JD003981]

Mao J T, Liu L and Zhang J H. 2001. GMS5 remote sensing of aerosol optical thickness over chaohu lake. Acta Meteorologica Sinica, 59 (3): 352–359. [DOI: 10.11676/qxxb2001.037] ( 毛节泰, 刘莉, 张军华. 2001. GMS5卫星遥感气溶胶光学厚度的试验研究. 气象学报, 59 (3): 352–359. [DOI: 10.11676/qxxb2001.037] )

Moulin C F, Guillard F, Dulac and Lambert C E. 1997. Long-term daily monitoring of Saharan dust load over ocean using Meteosat ISCCP-B2 data: 1. Methodology and preliminary results for 1983–1994 in the Mediterranean. Journal of Geophysical Research: Atmospheres, 102 (D14): 16947–16958. [DOI: 10.1029/96JD02620]

Prados A I, Kondragunta S, Ciren P and Knapp K R. 2007. GOES Aerosol/Smoke Product (GASP) over North America: comparisons to AERONET and MODIS observations. Journal of Geophysical Research: Atmospheres, 112 (D15): D15201 [DOI: 10.1029/2006JD007968]

Remer L A, Kaufman Y J, Tanré D, Mattoo S, Chu D A, Martins J V, Li R R, Ichoku C, Levy R C, Kleidman R G, Eck T F, Vermote E and Holben B N. 2005. The MODIS aerosol algorithm, products, and validation. Journal of the Atmospheric Sciences, 62 (4): 947–973. [DOI: 10.1175/JAS3385.1]

Sekiyama T T, Yumimoto K, Tanaka T Y, Nagao T, Kikuchi M and Murakami H. 2016. Data assimilation of Himawari-8 aerosol observations: Asian dust forecast in June 2015. Scientific Online Letters on the Atmosphere Sola, 12 : 86–90. [DOI: 10.2151/sola.2016-020]

Shang H Z, Chen L F, Letu H, Zhao M, Li S S and Bao S H. 2017. Development of a daytime cloud and haze detection algorithm for Himawari-8 satellite measurements over central and eastern China. Journal of Geophysical Research: Atmospheres, 122 (6): 3528–3543. [DOI: 10.1002/2016JD025659]

Uesawa D. 2016. Aerosol Optical Depth product derived from Himawari-8 data for Asian dust monitoring[DB/OL].

Wang J, Christopher S A, Reid J S, Maring H, Savoie D, Holben B N, Livingston J M, Russell P B and Yang S K. 2003. GOES 8 retrieval of dust aerosol optical thickness over the Atlantic Ocean during PRIDE. Journal of Geophysical Research: Atmospheres, 108 (D19): 8595 [DOI: 10.1029/2002JD002494]

Wang X Q, Yang S Z, Zhu Y H and Yi W N. 2003. Aerosol optical thickness retrieval over land from MODIS data based on the inversion of the 6S model. Chinese Journal of Quantum Electronics, 20 (5): 629–634. [DOI: 10.3969/j.issn.1007-5461.2003.05.025] ( 王新强, 杨世植, 朱永豪, 易维宁. 2003. 基于6S模型从MODIS图像反演陆地上空大气气溶胶光学厚度. 量子电子学报, 20 (5): 629–634. [DOI: 10.3969/j.issn.1007-5461.2003.05.025] )

Wang Y S, Zhang J K, Wang L L, Hu B, Tang G Q, Liu Z R, Sun Y and Ji D S. 2014. Researching significance, status and expectation of haze in Beijing-Tianjin-Hebei Region. Advances in Earth Science, 29 (3): 388–396. [DOI: 10.11867/j.issn.1001-8166.2014.03.0388] ( 王跃思, 张军科, 王莉莉, 胡波, 唐贵谦, 刘子锐, 孙扬, 吉东生. 2014. 京津冀区域大气霾污染研究意义、现状及展望. 地球科学进展, 29 (3): 388–396. [DOI: 10.11867/j.issn.1001-8166.2014.03.0388] )

Wang Z T, Xin J Y, Jia S L, Li Q, Chen L F and Zhao S H. 2015. Retrieval of AOD from GF-1 16 m camera via DDV algorithm. Journal of Remote Sensing, 19 (3): 530–538. [DOI: 10.11834/jrs.20154110] ( 王中挺, 辛金元, 贾松林, 厉青, 陈良富, 赵少华. 2015. 利用暗目标法从高分一号卫星16m相机数据反演气溶胶光学厚度. 遥感学报, 19 (3): 530–538. [DOI: 10.11834/jrs.20154110] )

Yu F F and Wu X Q. 2016. Radiometric inter-calibration between Himawari-8 AHI and S-NPP VIIRS for the solar reflective bands. Remote Sensing, 8 (3): 165 [DOI: 10.3390/rs8030165]

Yumimoto K, Nagao T M, Kikuchi M, Sekiyama T T, Murakami H, Tanaka T Y, Ogi A, Irie H, Khatri P, Okumura H, Arai K, Morino I, Uchino O and Maki T. 2016. Aerosol data assimilation using data from Himawari-8, a next-generation geostationary meteorological satellite. Geophysical Research Letters, 43 (11): 5886–5894. [DOI: 10.1002/2016GL069298]

Zhang J H, Si Z J, Mao J T and Wang M H. 2003. Remote sensing aerosol optical depth over China with GMS-5 satellite. Chinese Journal of Atmospheric Sciences, 27 (1): 23–35. [DOI: 10.3878/j.issn.1006-9895.2003.01.03] ( 张军华, 斯召俊, 毛节泰, 王美华. 2003. GMS卫星遥感中国地区气溶胶光学厚度. 大气科学, 27 (1): 23–35. [DOI: 10.3878/j.issn.1006-9895.2003.01.03] )

Zhang T Y, Lin W P, Chen J Z and Zong W. 2009. Comparative study on atmospheric correction methods of FLAASH and 6S model based on Spot 5. Journal of Optoelectronics·Laser, 20 (11): 1471–1473. [DOI: 10.16136/j.joel.2009.11.032] ( 张婷媛, 林文鹏, 陈家治, 宗玮. 2009. 基于FLAASH和6S模型的Spot 5大气校正比较研究. 光电子·激光, 20 (11): 1471–1473. [DOI: 10.16136/j.joel.2009.11.032] )

Zou X, Zhuge X and Weng F. 2016. Characterization of bias of advanced himawari imager infrared observations from NWP background simulations using CRTM and RTTOV. Journal of Atmospheric and Oceanic Technology, 33 (12): 2553–2567. [DOI: 10.1175/JTECH-D-16-0105.1]

摘要

关键词