武汉地区环境卫星气溶胶光学厚度反演

刘宇昕，王斌，曹诗颂，张峰，王兆徽

(1.国家卫星海洋应用中心，北京 100081；2.自然资源部 空间海洋遥感与应用重点实验室，北京 100081；3.河海大学 海洋学院，南京 210098；4.中国气象局武汉暴雨研究所，武汉 430074；5.北京建筑大学 测绘与城市空间信息学院，北京 102616；6.北京集思睿景科技有限公司，北京 102401)

0 引言

近百年来，伴随着社会工业化的发展以及城市化进程的加速，多种大气污染源的过度排放，使城市空气状况受到极大影响，已经远远超过城市生态承载力。大气污染对城市人群健康以及建筑等各方面产生的影响，已成为我国城市和区域可持续发展需要面对的一个难题。同时，区域性的大气污染也已成为众多学者关注的重点，诸如京津冀、珠三角、长三角等城市群的大气污染越来越受到人们重视，并已开展了大量相关研究。此外，在中部地区如武汉及周边区域，工业化和城市化带来的大气污染也是环境管理部门所面临的棘手问题。在诸多大气污染物之中，气溶胶与大气环境状况密切相关，气溶胶颗粒的增加可引发雾霾天气的产生，严重影响人们的身体健康和正常生活。

卫星遥感可以迅速、准确、动态地监测较大区域的大气状况及变化，用以及时研究相关应对方案和规划，降低污染给人们所带来的影响。气溶胶光学厚度(aerosol optical depth，AOD) 作为衡量城市空气质量的重要指示因子，对大气污染情况的判定以及定量遥感研究都具有很大意义。作为气溶胶观测的一种有效方法，气溶胶光学厚度遥感定量反演研究至今已有40多年的历史，多年来针对不同种类的传感器，开发出了多种气溶胶光学厚度反演算法。如：基于改进的甚高分辨率辐射计(advanced very high resolution radiometer，AVHRR)的单通道和双通道气溶胶光学厚度反演算法；基于总臭氧测绘光谱仪(total ozone mapping spectrometer，TOMS)的气溶胶指数和近紫外光通道算法；基于宽视场海洋观测传感器(sea-viewing wide field-of-view sensor，SeaWiFS)、中分辨率成像光谱辐射计(moderate resolution imaging spectroradiometer，MODIS)、中分辨率成像光谱仪(medium-resolution imaging spectrometer，MERIS)的多通道算法；基于偏振探测器(polarization and directionality of the Earth’s reflectance，POLDER)，以及暗目标法、深蓝算法、结构函数法、MAIAC(multi-angle implementation of atmospheric correction)算法[1]、SARA(simplified aerosol retrieval algorithm)算法[2]等气溶胶光学厚度反演算法。结构函数法是目前普遍采用的一种气溶胶反演算法，其在陆地卫星专题测图仪(landsat thematic mapper，Landsat TM)、MODIS、AVHRR以及环境(HJ)卫星数据中都得到了广泛的应用。

使用结构函数法反演气溶胶光学厚度至今已有30多年的历史。在国外研究方面，Tanré等[3]率先提出了结构函数算法用于干旱地区的气溶胶滚光学厚度反演，并使用该方法成功反演出Saharan沙漠地区的气溶胶光学厚度；Holben等[4]使用结构函数法实现了AVHRR等数据的气溶胶光学厚度反演；Levy等[5]使用MODIS数据反演气溶胶光学厚度，并对结构函数法进行了相应的改进。国内研究方面，Liu等[6]使用结构函数法实现了地球观测卫星系统(systeme probatoire d’observation de la terre，SPOT)数据的气溶胶光学厚度反演并判定大气污染状况；李晓静等[7]利用MODIS数据运用结构函数法，成功反演得到北京和周边地区的气溶胶光学厚度；黄艇[8]以兰州地区为例，结合MODIS数据，使用结构函数法，反演得到了该地区的气溶胶光学厚度，并与多种方法的结果进行比较；孙林[9]使用改进的结构函数法，并结合双向反射分布函数(bidirectional reflectance distribution function，BRDF)模型，针对MODIS数据成功反演得到北京地区的气溶胶光学厚度，提高了反演精度；周春艳等[10]针对高分辨率HJ卫星数据使用结构函数法获取了北京地区的气溶胶光学厚度；朱琳等[11-12]研究了使用已有的地表反射率产品为待反演图像提供清晰图像的方法，以推动结构函数法的业务化应用，并研究获取了最佳的像元间隔设置。武汉地处中部，地形包括低山、丘陵、平原等，土地利用类型多样且利用程度较高，湖泊塘堰众多，工业发达且市内较多重工业，大气污染较严重，但武汉这一较具代表性区域还缺乏气溶胶光学厚度反演相关研究。相关研究已经证明，对于城市等亮地表地区，结构函数法在反演TM、MODIS等不同分辨率载荷的AOD方面都有较好的结果。本文选择武汉市作为研究区域，使用结构函数法反演环境一号(HJ-1)A/B卫星的气溶胶光学厚度，并对结果进行了验证。

1 结构函数法反演模型

由于卫星传感器接收到的辐射能量来源于大气散射以及地表反射，当地表反射很少时，卫星观测到的信号主要来自于大气部分；当地表反射增大时，地表反射又成为主要部分。通过卫星遥感数据反演气溶胶光学厚度方法的实质，就是将来自于大气和地表反射的辐射分离开来，也就是地气解耦。

在地表反射率较高的区域，地表反射在卫星接收到的能量中占主要部分，大气气溶胶的散射较弱，所以大气气溶胶信息会因使用地表反射率模型带来的误差而丢失。不同地点地表覆盖类型的差异,通过遥感反射率在空间上的变化率反映出来，就是此区域的结构函数。结构函数法通过大气透过率得到气溶胶光学厚度，能够较好解决地表反射率模型误差带来的问题。

1.1 结构函数法原理

使用结构函数法反演一组遥感数据，其中包括了一张晴天无云状态下气溶胶光学厚度非常小的影像，通过计算其像元间地表反射率变化率后作为参考图像，代表该地区地表反射率的分布规律。对于该组中其他影像，假设该地表空间分布规律在此组遥感数据对应时间内保持不变，则该时间段内遥感数据表观反射率在空间分布的变化是由气溶胶所引起，依据其分别与参考影像的差别，可以对气溶胶光学厚度进行计算。结构函数法由于较少受到地面反射率的限制，因此在干旱和半干旱地区的应用较为广泛，原理如下[3]：假设大气水平均匀，陆地表面为朗伯体，信号用表观反射率ρ*表示，定义如式(1)所示。

(1)

式中：θs、θv、φv是观测几何参数中的太阳天顶角、观测天顶角和方位角；ρa(θs，θv，φv)是大气反射率；T(θs)是大气下行透过率；td(θv)为环境辐射大气上行透过率；S是大气向下半球反照率；ρc(M)是目标物M的表面反射率；<ρc(M)>是平均地表反射率；τ是大气光学厚度；μv为cosθv。

(2)

(3)

在上述假设条件下，Δρi，j是一个常数，表观反射率是光学厚度的一个函数，给出一个地表结构函数见式(4)。

(4)

式中：N为区域总像元数，通过式(3)和式(4)，可得式(5)。

(5)

式中：M*和M分别表示卫星观测数据和真实地表情况的结构函数。当t1和t2的地表反射率没有发生变化，M2(d，t1)=M2(d，t2)，根据公式(5)可得式(6)。

(6)

根据经验，ln[T(θs)exp(-τ/μv)]与气溶胶光学厚度成线性关系[3]。因此，当已知一个时刻的气溶胶光学厚度时，即可以计算其他时刻的气溶胶光学厚度。

1.2 结构函数模型

在结构函数法反演气溶胶光学厚度的过程中，非常重要的一项内容就是确定结构函数模型，包括结构函数公式、窗口大小以及距离值的确定。其反演原理是以朗伯体假设为出发点，使得地表反射特性不受观测几何或地形的影响；然而实际中大部分地表是非朗伯体，意味着气溶胶光学厚度的反演结果受地形、观测几何等影响。公式(4)定义的结构函数只考虑了一个方向，不能很好地对地表结构进行建模，可以通过增加(ρi，j-ρi+d，j)2和(ρi，j-ρi+d，j+d)2项，将表面结构模型从一维扩展到多维，以减小模型误差，包含了行、列以及对角线方向，如公式(7)所示[6]。

(7)

式中：m、n确定结构函数窗口范围；d确定计算结构函数值时两像元之间的距离。很明显，不同距离d下的气溶胶光学厚度计算结果会有差别，计算过程中需考虑不同的结构函数窗口，并选择不同的距离值进行匹配。

2 武汉地区气溶胶光学厚度反演实例

2.1 数据源

本文使用的HJ-1A/B卫星CCD(charge-coupled device)数据，来源于环境减灾卫星，是我国专门用于环境与灾害监测的小卫星星座，星座中包含2颗光学卫星(HJ-1 A星和HJ-1 B星)，这2颗光学卫星上分别装有2台宽覆盖多光谱可见光相机(简称CCD相机)。另外，HJ-1 A星上装有一台超光谱成像仪，HJ-1 B星上装有一台红外扫描仪。环境一号系列卫星的宽覆盖多光谱相机、红外相机、超光谱成像仪的重访观测周期分别为48、96和96 h，基本具备了在全国范围内开展环境和灾害大范围、全天候观测的能力。

2.2 数据预处理流程

结构函数法反演HJ-1 A/B卫星遥感影像的预处理主要有辐射定标、计算表观反射率、参考影像大气校正、几何校正、波段合成以及生成6S查找表等。

1)辐射定标是指建立遥感传感器的数字量化输出值(digital number，DN)与其所对应视场中辐射亮度值之间的定量关系，HJ卫星数据的定标参数来自中国资源卫星应用中心[13]。首先利用绝对定标系数，将CCD影像DN值转换为辐亮度影像，然后将辐亮度转化为表观反射率。

2)大气校正部分借助ENVI(the environment for visualizing images)软件提供的快速大气校正工具(quick atmospheric correction，QUAC)，通过从遥感图像上采集不同地物的波谱特征，依靠得到的经验值实现大气校正。目前QUAC支持的多光谱和高光谱波谱范围是400～2 500 nm，输入数据可以是辐射亮度值、表观反射率、无单位的raw数据，数据储存格式和类型没有特殊要求，但必须提供多光谱和高光谱传感器数据每个波段的中心波长信息。QUAC可以利用很少的参数快速地对多光谱和高光谱数据进行处理，自动化程度较高，已被很多专家学者采用[14]。

3)几何校正部分仅针对结构函数法，原因是结构函数法是只针对一个时间序列内的影像进行气溶胶光学厚度反演的算法，这就要求该时间序列内的所有影像必须在地理位置上可以很好地重合，尽量保证不同影像的几何偏移在1个像元之内。几何校正部分使用ERDAS 软件的自动几何校正模块实现，具体流程是：首先选取一景几何校正后效果较好的遥感影像作为标准，并对其他影像进行校正。操作过程中，首先要手动选取3个或3个以上地理位置比较明确的控制点，借助自动校正模块，选取其余的控制点，并将误差超过1个像元的控制点删除，以此为准，对该时间序列内的所有影像进行几何校正，误差控制在0.5个像素以内。

4)波段合成的目的是通过遥感影像不同波段进行组合，从而突出与气溶胶光学厚度反演相关的信息。对于结构函数法，选择HJ卫星的第1(蓝光波段)、第3(红光波段)、第4(近红外波段)波段进行合成，然后通过影像裁剪，获得武汉地区的多波段合成结果。

5)查找表法。6S模式是“第二代太阳短波辐射的卫星信号模拟(second simulation of satellite signal in the solar spectrum，6S)”的简称，由法国大气光学实验室和美国马里兰大学地理系Veromote E，在太阳光谱波段卫星信号模拟程序(simulation of satellite signal in the solar spectrum，5S)模型的基础上开发的辐射传输模型，目前己成为世界上发展最完善的大气辐射校正模型之一[15]。该模式主要用来模拟星载或机载遥测仪器在0.25～4.0 μm光谱上无云情况下卫星传感器理论上应接收到的辐射值，模式的光谱积分步长为2.5 nm[8]。6S大气辐射传输模式主要包括观测几何参数、大气参数、气溶胶参数、传感器的光谱特性以及地表反射率。根据以下6S参数设置生成查找表：9个太阳天顶角为0°、6°、12°、24°、35.2°、48°、54°、60°和66°；12个观测天顶角为0°、6°、12°、18°、24°、30°、36°、42°、48°、54°、60°、66°；16个相对方位角为0°、12°、24°、36°、48°、60°、72°、84°、96°、108°、120°、132°、144°、156°、168°、180°；大气模式为中纬度夏季和中纬度冬季；大气气溶胶模式为大陆型、城市型以及自定义类型；气溶胶光学厚度预设值分别为0、0.25、0.5、1.0、1.5和1.95。

2.3 武汉地区气溶胶光学厚度反演

在遥感数据预处理结果的基础上，反演过程中需要使用到HJ卫星CCD影像第1、3、4波段的表观反射率、遥感影像元数据以及生成的查找表，通过编程实现结构函数法反演气溶胶光学厚度。通过数据查询，选取2012年9月1日至2012年10月2日中8 d的数据作为基础数据集，具体的数据状况如表1所示。

表1 反演数据列表

以上数据集在时间序列覆盖整个9月份，数据基本无几何畸变，且云覆盖量很少。通过在湖北省环境保护厅网站查询得到这段时间内武汉市的大气污染状况，如图1所示，以此为基础，获取相应时间的空气污染指数(air pollution index，API)。根据空气污染指数，选择2012年9月15日为“清洁日”，标准是选取一段时间序列内空气质量次优的一天。即空气污染指数第二小的日期作为“清洁日”。其原理是“清洁日”要求大气洁净，大气污染程度较轻，同时为避免观测数据的随机性和偶然性带来的影响，因此放弃API最小的日期，选择API第二小的日期作为“清洁日”。

图1 2012年9月1日至2012年10月2日武汉大气污染状况

以2012年9月20日为例，依次对原始影像进行辐射定标、表观反射率计算、波段融合、几何校正和裁剪等操作，对于“清洁日”的影像，裁剪完成后需要进行大气校正得到参考影像。实验结果表明，武汉地区由于水体面积覆盖范围较大，河湖纵横，水体信号对于气溶胶光学厚度的反演造成了很大的影响，导致长江、东湖、南湖等水体面积较大区域的气溶胶光学厚度明显偏高。并且由于结构函数法反演气溶胶光学厚度是基于窗口进行反演，这样水陆边界区域就存在大量的混合像元，这些地区的气溶胶光学厚度值也偏高，因此需要寻求合理的方案解决水体对于气溶胶光学厚度反演过程的影响。

剔除水体影响的具体方法是使用归一化植被指数(normalized difference vegetation index，NDVI)进行水体和陆地的判别，并对影像进行掩膜，从武汉地区的表观反射率影像中提取陆地部分。经过处理得到9月15日和9月20日的武汉市地区陆地区域的表观反射率影像，如图2所示。

图2 2012年9月15日与2012年9月20日武汉地区陆地区域的表观反射率合成影像

以图2所示的2景影像为基础，首先要根据表观反射率,计算得到2 d的结构函数值，经过比较分析，结构函数窗口选择10、距离值选择5，为最佳的参数组合，计算得到2天的结构函数值影像。接下来，对查找表进行插值，结构函数法查找表插值计算过程中需要使用到查找表中的太阳天顶角、观测天顶角、大气下行辐射透过率3项参数，通过插值运算结合结构函数值影像进行计算，最终得到气溶胶光学厚度。

计算得到的2012年9月20日的气溶胶光学厚度如图3所示，水体部分已经全部掩膜，对应于图中的白色部分，陆地部分的气溶胶光学厚度值介于0～2.0之间。

图3 2012年9月20日气溶胶光学厚度

3 验证

对于以上反演结果，本文通过2种方式进行验证：一是依据全自动太阳光度计CE318实测数据对反演结果进行精度验证；二是依据湖北省环保厅提供的武汉地区大气污染状况，从时间上的分布趋势对反演结果进行验证。

1)实测数据验证。2012年9月1日、2012年9月5日、2012年9月18日、2012年9月19日、2012年9月20日、2012年9月28日、2012年10月2日的气溶胶光学厚度如图4所示，选取的实测数据与影像进行匹配的时间窗口小于1 h，验证情况如表2所示。

经过以上验证可知，本次反演结果的绝对误差均在0.15以内，相对误差均在15%以内，表明反演

图4 气溶胶光学厚度反演结果

表2 气溶胶光学厚度反演精度验证

结果具有较高精度。同时，从空间分布特征可以看出，气溶胶光学厚度较大的区域主要包括市区以及长江沿岸，并且武汉三镇中汉口和汉阳地区较高，这一趋势与武汉市周边地区的实际状况基本一致。

2)时间分布趋势验证。针对反演结果，对武汉市地区的气溶胶光学厚度的平均值进行排序，得出气溶胶光学厚度从大到小的分布情况为：2012年9月19日>2012年9月20日>2012年10月2日>2012年9月18日>2012年9月1日>2012年9月28日>2012年9月5日。结合湖北省环保厅查询到的大气污染状况数据，可得到2012年9月1日、2012年9月5日、2012年9月18日、2012年9月19日、2012年9月20日、2012年9月28日、2012年10月2日的API分别为：88、66、83、119、105、71、85。

以上述数据为基础，得到AOD与API的分布趋势如图5所示。其中，横轴代表时间，纵轴(左)代表API，纵轴(右)代表AOD。可以看出，气溶胶光学厚度的反演结果与大气污染状况的时间分布趋势基本一致，从时间序列分布状况反映了反演结果的正确性。

图5 气溶胶光学厚度时间序列分布趋势

3)与MODIS气溶胶光学厚度产品验证结果对比。根据国家卫星气象中心运行的MODIS气溶胶遥感产品质量检验分析结果，1 km分辨率MODIS L1B AOD产品的绝对误差平均值为0.041 1，均方根误差为0.27[16]。本文反演结果的绝对误差平均值为0.051 1，精度略低于MODIS产品，均方根误差为0.057 3，准确度高于MODIS产品。综合考虑，反演结果与1 km分辨率MODIS气溶胶产品具有相当的质量评价，但较高的时空分辨率决定其具有更高的细致性评价。

4 结束语

本文以武汉及周边为研究区域，选取适用于干旱、半干旱以及城市等亮地表区域的结构函数法，进行气溶胶光学厚度反演，通过对结构函数法的理论基础和模型进行研究，确定了结构函数公式、反演窗口和距离值。然后编程实现对武汉地区2012年9月至10月的8景HJ卫星CCD影像的气溶胶光学厚度反演。并结合CE318实测数据、湖北省环境保护厅网站提供的武汉市大气污染状况及分布趋势数据、MODIS气溶胶产品质量检验分析结果，从多方面对反演结果进行了综合验证。HJ-1 A和HJ-1 B双星联合观测情况下，全球覆盖的时间分辨率是2 d，空间分辨率30 m，与MODIS和TM等载荷相比，具有较高的时空分辨率，适用于结构复杂、气溶胶厚度变化较快的城市地区业务化监测，是一个较有意义的尝试。且验证结果表明，产品具有较高精度，与实际特征相一致，使用结构函数法进行HJ-1 A/B卫星的CCD影像气溶胶光学厚度反演在武汉地区具有较高的可靠性，可应用于该地区的气溶胶分布和变化趋势分析，以及大范围的区域大气环境监测等。