基于天宫二号宽波段成像仪数据的气溶胶光学厚度遥感反演研究

杨 越，陈 健，崔嘉文

(南京信息工程大学遥感与测绘工程学院，南京210044)

1 引言

气溶胶一般是指由固体或液体颗粒悬浮在大气介质中形成的多相体系，其粒子直径在10-3～102μm之间。作为大气重要组成成分，气溶胶不仅影响大气环境与人类健康，还在地气系统中发挥着重要作用，对气候变化也有着重要的影响[1]。利用遥感技术获取大范围的气溶胶空间分布具有重要的理论意义和应用价值。

目前，已有研究人员开展了气溶胶的遥感反演工作。曾晶等[2]基于MODIS数据反演了长株潭地区气溶胶光学厚度(Aerosol Optical Depth，AOD)；郑尧等[3]针对HJ-1卫星CCD数据，反演了太湖地区AOD。在基于静止卫星反演AOD方面，张玉环等[4]利用GOCI数据对北京及周边城市的气溶胶状况进行了监测；葛邦宇等[5]利用Himawari-8数据反演了京津冀地区AOD；Husar等[6]利用AVHRR传感器的单通道信息获得了海洋上空的气溶胶光学厚度；Kaufman等[7]针对MODIS数据提出了经典的暗像元方法，成为全球MODIS AOD产品的业务化算法。以上针对极轨卫星、静止卫星开展的研究已经相当丰富，但目前针对空间实验室的气溶胶遥感研究仍然非常缺乏。

天宫二号是我国自主研发的第二代空间实验室，其搭载的宽波段成像仪在世界上首次实现了可见近红外、短波红外与热红外多光谱大视场全推扫成像的组合集成功能，具有空间分辨率高、数据新颖与数据获取个性化的特点，是对地观测的重要遥感工具，也为气溶胶遥感反演、大气环境监测提供了一种新的观测方式。

为填补气溶胶的空间站遥感研究的空白，本文选取黄河三角洲地区为研究区，基于天宫二号宽波段成像仪获取的可见近红外波段遥感影像，结合MODIS地表反射率产品，利用深蓝算法反演研究区上空气溶胶光学厚度，然后利用Himawari静止卫星的L2级AOD产品对反演结果进行对比研究，以探讨天宫二号数据用于AOD反演的可行性。

2 研究区与研究数据

2.1 研究区概况

黄河三角洲(38°15′N，118°5′E)，是指黄河入海口携带泥沙在渤海凹陷处沉积形成的冲积平原，面积约5400 km2，大部分区域隶属于山东省东营市。黄河三角洲地处中纬度地区，背陆面海，气候同时受到大陆与海洋的影响，属于暖温带半湿润大陆性季风气候区，年平均气温11.7～12.6℃，四季分明。研究区范围如图1所示。

2.2 研究数据

2.2.1 天宫二号数据

图1 黄河三角洲示意图(天宫二号宽波段成像仪RGB4、10、12 波段假彩色合成)Fig.1 Schematic diagram of Yellow River Delta(RGB4/10/12 false color composition of Tiangong-2 wide-band imager)

天宫二号搭载的宽波段成像仪具有图谱合一的特点，在获取影像信息的同时也获取地物光谱信息。可见近红外波段可获取14个通道(表1)，空间分辨率100 m，幅宽300 km的多光谱遥感影像。数据从载人航天空间应用数据推广服务平台官网(http：//www.msadc.cn)免费下载。本文选用2017年3月18日天宫二号宽波段成像仪可见近红外波段数据进行气溶胶光学厚度反演，影像范围包括整个黄河三角洲地区。由于获取的天宫二号原始数据为增大某一倍数的辐亮度值，因此除以头文件中所给放大倍数得到真实辐亮度值。

表1 天宫二号宽波段成像仪波段设置及参数信息Table 1 Parameter information and band settings of Tiangong-2 wide-band imager

2.2.2 MODIS数据

中分辨率成像光谱仪MODIS目前主要搭载于两颗卫星上：Terra和Aqua，具有分布在0.4～14 μm电磁波谱范围内的36个光谱通道。MODIS数据分为0～5级的6类产品[8]，产品的分辨率包括4类：250 m、500 m、1000 m、以及 5600 m(0.05°)。本文使用MODIS 8天合成地表反射率产品MOD09A1来构建适用于天宫二号宽波段成像仪的地表反射率数据库，由于天宫二号遥感影像与MOD09A1影像空间分辨率不同，考虑到重采样图像质量问题，将二者重采样为空间分辨率200 m的影像。然后针对二者蓝光波段波长范围与光谱响应函数存在的差异交叉辐射校正，以修正光谱响应不同而造成的误差[9]。

2.2.3 Himawari-8 AOD产品

Himawari-8是日本于2014年发射的静止气象卫星，提供了2～3级的多种AOD产品[10]，本文选用空间分辨率为5 km的L2级AOD产品与利用天宫二号数据反演得到的气溶胶光学厚度结果进行比较，验证结果精度。

3 研究方法

3.1 AOD反演原理

在进行AOD反演时，假设地表为朗伯体，并且大气水平均一条件下卫星传感器所接收得到的表观反射率ρTOA可以表示为式(1)：

式中μS与μV分别代表太阳天顶角θS的余弦和观测天顶角θ的余弦，φ则为相对方位角，ρ0为大气程辐射项的等效反射率，R表示朗伯体的地表反射率，S为大气整层向下路径的半球反射率，T为大气透过率。 其中，ρ0、S和T(μS)T(μV) 这3个参数都携带了大气气溶胶的信息，R反映了地表反射特性[11]。从公式中可以看出表观反射率耦合了大气与地表的信息，在进行AOD反演之前要进行地气解耦，去除地表贡献。

目前，针对地气解耦主要有两种方法来确定地表反射率，即暗像元法和深蓝算法。由于暗像元算法主要适用于存在大量浓密植被的区域，而本研究中遥感影像的成像时间为3月份，并非植被长势最茂盛的季节；且该区域地表状况复杂，盐碱化程度较高，地表反射率较高，因此不宜采用暗像元算法。相关文献表明，深蓝算法更适用于较广泛的地表类型去除地表反射率的影响，且应用于AOD反演具有较好的结果[12-14]。鉴于研究区地表土地利用较复杂且地表反射率较高的特点，本文利用深蓝算法反演研究区AOD。在地表反射率已知的情况下，利用6S辐射传输模型可以得到大气3参数、ρ0、S和T(μS)T(μV)，并进一步构建查找表获取影像每个像元AOD值[15-16]。

3.2 AOD反演

相对于 AVHRR、MODIS及葵花8号、GOCI静止卫星数据等各种用于气溶胶反演的遥感数据，天宫二号宽波段成像仪影像具有更高的空间分辨率，可提供区域尺度精细的气溶胶分布信息。由于天宫二号数据集中没有包含地表反射率产品，因此本研究选用MOD09A1地表反射率产品来构建适用于天宫二号的地表反射率库，从而进行地气解耦。并计算太阳角度与卫星观测角度，建立基于6S辐射传输模型的查找表。反演技术流程如图2所示，具体流程如下：

1)计算表观反射率。基于天宫二号宽波段成像仪可见近红外波段卫星影像，利用ENVI 5.3波段运算(Band Math)计算表观反射率如式(2)：

式中L为影像中各像元的辐亮度值；D为成像时刻的平均日地距离订正因子；Esun为影像该波段的太阳辐照度；θS为太阳天顶角。

2)构建地表反射率库。在利用MOD09A1地表反射率产品构建适用于天宫二号宽波段成像仪时，两者蓝波段波长范围不同，根据官方发布信息，天宫二号 12波段范围为 480～500 nm，而MOD09A1第三波段范围为459～479 nm。由于两者蓝光波段存在光谱响应差异，故而进行波段转换以减小不同光谱响应的差异。利用ENVI 5.3中自带的5种波谱数据库，选取典型的23种地物光谱数据，利用公式(3)所示蓝波段地表反射率转换公式计算得到2种传感器蓝光波段地表反射率。

图2 AOD反演流程图Fig.2 Flow chart of AOD inversion

式中，λ1、λ2为波长积分上下限，(S)λ为波长为λ时的传感器光谱响应函数，R(λ)为波长λ处波谱库中地物反射率。

根据计算得到的两种传感器蓝波段地表反射率可以建立线性回归模型，得到交叉辐射订正公式(4)进行光谱转换，以得到适用于天宫二号宽波段成像仪蓝波段的地表反射率。

两传感器线性回归模型如图3所示。

图3 MODIS和TG-2数据蓝光波段地表反射率散点图Fig.3 Scatter diagram of surface reflectance in blue light band of MODIS and TG-2 data

3)角度数据准备。深蓝算法中所需的角度数据包括观测天顶角、观测方位角、太阳天顶角、太阳方位角4个角度数据。由于天宫二号采用推扫式成像，所以星下点像元的观测天顶角为0°，并以星下点为中心沿垂直于卫星飞行方向向两侧递增至最大观测天顶角约21°；对于观测方位角，前向观测方位角与后向观测方位角相差180°，且原始数据提供了前向观测方位角，以此在IDL环境下逐像元计算出观测天顶角、方位角。对于太阳方位角、天顶角的计算，可根据成像时间、各像元的位置计算得到。

4)重采样及研究区裁剪。对天宫二号宽波段成像仪影像、MOD09A1地表反射率产品，4个角度数据重采样为相同的200 m分辨率，并利用研究区矢量图进行研究区裁切。

5)6S模型构建查找表。本文研究的天宫二号数据在6S模型中没有预设参数，所以需要自定义观测几何参数、大气模式参数、气溶胶模式、传感器光谱特性与地表反射率等参数。参考陈健等[15]研究中太阳天顶角与卫星天顶角对不同AOD值的表观反射率较敏感，而太阳方位角和卫星方位角较迟钝的研究结果，以4°为间隔设置卫星天顶角，以1°为间隔设置太阳天顶角，以12°为步长设置相对方位角(相对方位角＝卫星方位角-太阳方位角)。由于黄河三角洲地区濒临渤海，受区域位置及风向影响，气溶胶类型以大陆型或海洋型气溶胶为主。研究区影像成像时刻风向数据表明[17-19]：成像时刻1500 m高度风向主要为北风，由海洋吹向内陆，因此研究区气溶胶类型选用海洋型气溶胶。下垫面特性设置为均一的朗伯体表面，地表反射率参数设置为基于MOD09A1构建的消除蓝波段光谱响应差异的天宫二号地表反射率库。

6)AOD反演。做完上述流程后，利用深蓝算法进行研究区AOD反演，得到研究区2017年3月18日当日气溶胶光学厚度反演结果。

4 AOD反演结果与精度验证

4.1 反演结果

根据以上流程反演得到如图4所示2017年3月18日黄河三角洲地区气溶胶光学厚度分布图，反演的AOD分布结果主要呈现东部较高西部较低，沿海较高内陆较低的分布情况，主要受到地理位置与土地利用类型的影响。黄河三角洲东部沿海地区油气、地热等自然资源丰富，多建设有大型油田、矿井、盐田等，作业期间存在颗粒物排放。借助临近海岸带的地理优势，该区还存在大量人工渔业养殖场，经济较为发达，人口密集，气溶胶排放增加，且近海地区土壤类型以滨海潮盐土为主，土壤盐碱化严重，不利于植被生长，植被覆盖率较低，不利于空气污染调节[20-21]。西部内陆土壤肥沃，农田、林地等较多，受国家政策调控工业用地少，植被覆盖较多，AOD值普遍较小[21]。沿海少部分地区由于存在大量湿地、滩涂及芦苇荡等未开发地区及生态保护区，生态环境较好，空气污染小，AOD值也相对较小。此外，沿海地区AOD较高，内陆AOD较低的情况也可能受到海洋季风的影响，内陆地区受季风影响较小，气溶胶扩散较慢，空气质量相对较高。

图5为2017年3月18日黄河三角洲地区的Himawari-8 AOD产品，其空间分辨率为5 km。由图可以看出，同一时相的Himawari-8 AOD产品的空间分布与基于天宫二号结果基本一致，说明利用天宫二号宽波段成像仪进行气溶胶的遥感反演是完全可行的。另外，通过图4与图5的对比可以发现，基于天宫二号反演的AOD结果具有较高的空间分辨率，更能体现AOD分布的细节情况，更有利于该区AOD的精细化监测。

图4 黄河三角洲2017年3月18日200 m分辨率AOD反演结果Fig.4 AOD inversion results of the Yellow River Delta at a 200 m resolution on March 18， 2017

4.2 精度验证与误差分析

4.2.1 精度验证

由于黄河三角洲地区缺少实测站点数据，本文利用其他卫星AOD产品数据进行了间接的精度验证。以往研究多与 NASA官方发布的MOD04产品进行对比，但由于其研究区内2017年3月18日的AOD产品缺值较多，质量较差，故本文利用研究区范围内2017年3月18日Himawari-8 AOD L2产品对反演结果进行了精度评价。

图5 2017年 3月18日 5 km分辨率 Himawari-8AOD产品Fig.5 Himawari-8 AOD products with 5 km resolution on March 18，2017

首先对Himawari-8 AOD L2产品进行投影转换与研究区矢量裁切等预处理，获取空间分辨率为5 km的AOD产品。为了与本文反演结果更好地进行精度验证，先将空间分辨率为200 m的AOD反演结果重采样为与Himawari-8 AOD产品一致的5 km分辨率，然后再将重采样后的天宫二号数据反演的AOD结果与其相比较，结果如图6所示：两者的线性拟合度为0.739，表明2种数据的离散程度较小，有较高的一致性。但从图中也可以看出，天宫二号数据反演的AOD值高于Himawari-8 AOD L2产品AOD值，造成这种误差的原因可能主要与表观反射率、地表反射率与气溶胶模式选取等因素有关。

图6 天宫二号数据反演结果与Himawari-8 AOD L2产品对比图Fig.6 Comparison between Tiangong-2 data inversion results and Himawari-8 AOD L2 products

4.2.2 误差分析

1)地表反射率影响。卫星所接收的表观反射率信号主要受地表反射与大气散射的影响，针对不同AOD取值，利用6S辐射传输模型分析了地表反射率对表观反射率的影响，结果如图7所示。从图中可以看到地表反射率在AOD值小于0.5时呈现较高的敏感性，表观反射率随气溶胶浓度的增大而显著增大；但当AOD值大于0.5时呈现较低的敏感性，表观反射率随AOD的变化程度较小。此外，随着地表反射率的增加，表观反射率对AOD的变化越来越不敏感。

图7 地表反射率与表观反射率对气溶胶的敏感性分析Fig.7 Sensitivity analysis of surface reflectivity and apparent reflectivity to aerosols

本文中利用MOD09A1地表反射率产品构建的天宫二号地表反射率值相对较大，结合图7并参考程晨等[22]的研究，当地表反射率偏大时，表观反射率对气溶胶变化敏感度下降，反演结果精度会受到影响。另外，本研究中由官网下载的用于AOD反演的影像是经过辐射定标的影像，利用该影像计算得到的表观反射率数值偏小，也是影响反演结果精度的一个原因。

2)朗伯体假设的影响。朗伯体假设对AOD的反演会有一定的影响，且在不同的波段影响也不同。一般来讲，与近红外波段相比，可见光波段的非朗伯特性则几乎可以忽略。此外，天宫二号宽波段成像仪的最大观测天顶角为21°，参考以往研究[23]，较小的观测天顶角下地物BRDF特征不明显。并且MODIS地表反射率产品也已经经过了二向反射率校正处理。综上所述，可认为本研究中朗伯体假设对AOD的精度影响并不明显。

3)气溶胶模型选取的影响。在利用6S辐射传输模型构建查找表时，气溶胶模型的选择对反演结果影响显著[15]。本研究利用6S模型模拟了城市型、大陆型与海洋型3种典型气溶胶类型下表观反射率随AOD的变化，结果如图8所示。由图可以看出，大陆型与海洋型气溶胶类型模式下随AOD的增加，表观反射率迅速增大，而城市型气溶胶模式下表观反射率变化较小。大陆型与海洋型气溶胶的变化趋势相对一致，但当AOD值大于1.1时海洋型模式下表观反射率随AOD增大的变化更显著。

图8 不同类型气溶胶的敏感性分析Fig.8 Sensitivity analysis of different types of aerosols

结合研究区地理位置、风向等因素，本文中气溶胶模式选择为海洋型气溶胶模式。为了验证选取海洋型气溶胶的合理性，本文又针对大陆型与海洋型气溶胶分别做了反演实验，结果如图9所示。可以看出海洋型气溶胶反演结果与AOD真实值非常接近，吻合度较高，而大陆型气溶胶反演结果略高于真实AOD值。

图9 不同类型气溶胶反演结果Fig.9 Inversion results of different types of aerosols

4)空间分辨率的影响。本文用于AOD反演的天宫二号宽波段成像仪数据的空间分辨率为200 m，而Himawari-8AOD空间分辨率为5 km，低分辨率在一定程度上会令AOD变化范围变小，如天宫二号200 m AOD产品的数值范围为0～2.136，而空间聚合为5 km后取值范围减小为0～1.753。

5)其他因素影响。由于Himawari-8AOD产品反映的是波长500 nm处AOD结果，而基于天宫二号的AOD反演结果为波长550 nm处的AOD结果。由于Himawari-8 AOD产品仅有一个波段，且官方提供的Angstrom指数和大气浑浊度指数还在验证中，因此无法将二者统一为同一波长处的AOD分布，故两者存在一定偏差，对验证结果也有一定影响。

5 结论

本文反演AOD结果与Himawari-8 AOD产品的空间分布基本一致，两者的拟合度达到0.739，说明利用天宫二号宽波段成像仪进行气溶胶的遥感反演是完全可行的。并且天宫二号反演的AOD结果具有较高的空间分辨率，更能体现AOD分布的细节情况，更有利于该区AOD的精细化监测。

天宫二号数据具有波段范围宽、空间分辨率高与时间分辨率灵活等特点，但由于其缺少地表反射率产品，在反演气溶胶光学厚度时利用其他数据构建地表反射率会对反演结果造成一定影响，未来若能利用其丰富的多波段信息反演地表反射率，将能更深层次挖掘天宫二号数据的应用潜能，实现数据推广并应用于我国海陆上空气溶胶监测甚至更多研究领域。

致谢：天宫二号宽波段成像仪数据产品下载于“载人航天空间应用数据推广服务平台”(www.msadc.cn)，感谢载人航天工程提供宽波段成像仪数据产品。感谢中科院空间应用工程与技术中心覃帮勇老师在数据使用及论文写作中给与的帮助与宝贵建议。