孙灏，刘伟汉，王艳梅，周伟，蔡创创

(中国矿业大学(北京) 地球科学与测绘工程学院，北京 100083)

0 引言

京津冀地区是世界上人口密度最大的地区之一。改革开放以来，京津冀及周边地区经济发展迅猛，然而，经济快速增长伴随着化学燃料的快速消耗，也导致了颗粒物浓度的逐年增加，频繁引起灰霾天气，影响着气候环境与人们的生活质量[1]。

灰霾的核心物质是空气中悬浮的灰尘颗粒，气象学上称为气溶胶颗粒[2]。气溶胶是悬浮于地球大气中，沉降速度小，半径在0.001～100 μm的液态及固态粒子，显著影响着气候环境。气溶胶光学厚度(aerosol optical depth，AOD)是定量描述气溶胶对电磁波衰减作用的物理量，是气溶胶最重要的参数之一。遥感技术能够弥补气溶胶地基观测站空间性和连续性的不足，实现大范围实时性的气溶胶监测。遥感器观测到的表观反射率是受大气和地表共同影响的结果，且不同的气溶胶类型对表观反射率的影响也不同。因此，气溶胶光学厚度的卫星遥感反演的2个核心问题分别是从卫星传感器接收到的信息中去除地表噪声以及确定合适的气溶胶模式。

传统的暗像元反演算法利用在晴朗无云的浓密植被等低地表反射率上空，气溶胶使反射率增大这一原理反演AOD，利用对气溶胶不敏感的短波红外通道表观反射率获得红蓝波段的地表反射率[3]。研究表明，暗像元算法仅在暗地表能获得较好的反演结果[4]，冬季华北平原大部分的植被消失，不能满足暗像元的假定，监测精度大大降低，难以满足业务化需求[5]。针对地表反射率较高的旱季和中高纬度地区冬季，Hsu等提出的基于地表反射率库的深蓝算法已得到成功应用[6-7]。2013年，改进的深蓝算法借助地表反射率数据与归一化植被指数NDVI，将反演范围从干旱半干旱亮表面扩展到了包括植被等暗地表的广泛陆地区域，并进一步讨论和改进了气溶胶模式与云雪检测算法[8]。MODIS在V5.2版大气气溶胶反演中加入了深蓝算法，并在C6版本提供了改进的深蓝算法产品以及暗像元与深蓝算法融合产品[9]，但与MODIS反射率数据(1 km)相比，分辨率仍较低，无法满足区域气溶胶的高精度监测需求[10]。不同地区的气溶胶类型、地表光谱特性各不相同，MODIS气溶胶产品且作为全球范围内的统一算法产品，对特定研究区域的针对性较差。本文利用融合互补思路，集成深蓝算法适应范围广和暗像元算法精度高的优势[11-12]，结合华北平原的自然地理特征，开展了一种暗像元算法与深蓝算法结合反演冬季京津冀地区AOD的方法，并进行了连续的AOD监测。依据监测结果，结合辅助资料，分析了京津冀地区冬季的气溶胶分布特征与其影响因素。

1 研究区域和数据

1.1 研究区概况

注：该图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2016)1666号的标准地图制作，底图无修改。图1 研究区域

京津冀是中国“首都圈”，总面积约21.8×104km2，常驻人口约1.1亿，其地理位置及气溶胶地基站点分布情况如图1所示。作为人口密度集中、工业发达的经济圈，京津冀地区地处华北平原，东临渤海，北有燕山山脉，西有太行山脉，这种复杂的地理环境作用形成的海陆风环流、山谷风环流对区域大气气溶胶的扩散有着较大影响，灰霾天气频繁发生。对京津冀气溶胶进行全方位、长时间的连续监测是有必要的。

1.2 数据来源

本研究采用的数据包括AQUA星MODIS L1B数据MYD021KM，该数据是对经过初步处理的一级产品进行定位和定标后所生成的二级产品，以HDF格式存储，空间分辨率1 km，时间分辨率1 d，对于太阳反射波段(对应MODIS波段1～19和波段26)提供反射率和辐亮度2种类数据产品；MYD021KM对应的地理定位数据MYD03；地表发射率库的构建选择8日合成的MODIS地表反射率产品MOD09A1。验证数据选择AERONET(aerosol robotic network)气溶胶监测网提供的北京、香河地基站AOD数据。此外，还使用了MODIS 10 km分辨率气溶胶产品MOD04_L2以及地理空间数据云提供下载的SRTMDEM 90 m分辨率原始高程数据。MODIS和AERONET数据分别通过NASA和AERONET官网下载。1月通常是京津冀地区冬季的最冷月份，将研究的时间跨度扩充到3年，总时间扩充到6个月，时间范围定为2016—2018年1—2月。

2 研究方法

2.1 反演原理

遥感反演AOD的基本原理是气溶胶粒子对电磁波散射的吸收作用使遥感器接收的辐射能的性质和强度发生变化，通过测量这种变化就可以反演气溶胶的特性[13]。

假设大气水平均一、地表朗伯体，遥感器观测到的表观反射率可表示为[5]：

(1)

式中：ρ0是大气路径辐射项等效反射率；μS=cosθS，μV=cosθV，θS和θV分别是太阳天顶角和观测天顶角；φ和ρS分别是相对方位角和地表双向反射率；S和T分别是为大气下界半球反射率和大气透过率，T(μS)T(μV)作为一个大气参数考虑。因此，地表反射率和大气气溶胶模式的确定是AOD反演的2个重点。

通过辐射传输模型计算S、ρ0、T(μS)T(μV)3个大气参数和AOD之间的对应关系，建立包含各参数的查找表，通过查找表获取AOD。

2.2 反演过程

1)构建查找表。查找表通过6S(second simulation of the satellite signal in the solar spectrum)辐射传输模型构建，使用IDL调用6S软件，设定7参数进行辐射传输计算，包括9个太阳天顶角，12个观测天顶角，16个相对方位角。设定0、0.025、0.5、1.0、1.5和1.95共6个在波长0.55 μm处的AOD以供内插。根据研究区域与监测时间段，大气模式选择中纬度冬季大气。针对气溶胶模式的选择问题，国际气象与大气物理协会(IAMAP)由气溶胶基本组分含量的不同定义了标准辐射大气(SRA)下的大陆型、城市型及海洋型3种基本气溶胶类型(表1)[14-15]。其中大陆型和城市型气溶胶分别表示了陆地自然条件下和严重人为排放条件下的气溶胶模型。京津冀地区的实际气溶胶类型可以理解为在自然状态的大陆型气溶胶背景下叠加了一定程度的人为排放，介于城市型和大陆型之间。在不考虑海洋性粒子的情况下使用正交试验法将沙尘粒子和水溶粒子占比各分为3个等级，共得到7组候选气溶胶模型(煤烟粒子占比介于1%到22%之间)，分别使用暗像元算法计算2017年1月区间内的气溶胶光学厚度，通过匹配AERONET北京站实测数据，得出了误差最小的京津冀地区冬季气溶胶类型作为本研究使用的气溶胶模型。其中，沙尘性、水溶性和煤烟性粒子分别占比44%、37%和19%。

表1 气溶胶模型及粒子组分

2)构建深蓝算法地表反射率库。针对在高地表反射率地区，红蓝波段和近红外波段地表反射率无法构成线性关系，通过构建地表反射率库实现地气解耦。本研究以MODIS地表反射率产品MOD09 8天合成产品为基础构建地表反射率库，按日期存储在数据文件中。

3)数据预处理。使用MODIS的L1B产品读取相关波段的表观反射率及相应的定标系数，使用地理定位数据读取经纬度、太阳天顶角、太阳方位角、卫星天顶角、卫星方位角和高程等信息。由定标系数将相应的数据转换成真实物理值。利用海陆掩码文件进行海陆分离，使用多波段阈值法去除云像元，进而进行AOD反演[16-17]。

4)气溶胶光学厚度反演。当2.1 μm波段的表观反射率小于0.4时，该像元可被认为是暗像元[5]。通过角度几何参数求得散射角，通过受气溶胶影响较小的中红外波段计算出NDVISWIR后[18]，由2.1 μm波段的反射率获得0.47 μm和0.66 μm波段的地表反射率：

(2)

式中：a、b是散射角与NDVISWIR的函数。

按照上式计算出红蓝波段的地表反射率；在构建的查找表中线性内插由地理定位数据读取的角度参数，得到不同波段和不同AOD的3个大气参数S、ρ0和T(μS)T(μV)；将反射率和大气参数代入式(1)得到表观反射率，再对真实的表观反射率进行插值，最终得到AOD。

深蓝算法原理与暗像元原理类似，根据过境时间提取对应的蓝波段地表反射率，利用查找表插值得到AOD。

2.3 合成处理

2种算法通过6S查找表均得到0.55 μm波长处的AOD。使用优选的融合方式对2种方法反演的AOD产品进行合成，以暗像元算法为主导，对暗像元算法未能反演的区域使用深蓝AOD进行填充，并使用以5×5模板窗口为基础的掩膜平滑处理。在窗口内以中心像素(i，j)为基准点，制作共3种形状的9个掩膜窗口，分别计算各窗口内的灰度值方差(图2)。该方法基于的原则是掩膜窗口方差越小，该窗口的像素值越均匀，中心像素属于该窗口的可能性越大。取方差值最小的掩膜窗口进行平均，用该均值代替中心像素点的灰度值。对于噪声点，该方法可以对其有效平滑，又能较好地保留面状地物的边界[19-20]。

图2 选择式掩膜窗口

方差的计算公式为：

(3)

均值的计算公式为：

(4)

式中：N为各掩膜对应的像素个数。

最终结果以Tiff格式输出，整体流程见图3。

图3 气溶胶反演流程

3 结果

本研究根据上文方法，对2016—2018年冬季1—2月的MODIS L1B数据对京津冀地区进行了连续监测。

3.1 精度验证

选择AERONET气溶胶监测网北京站和香河站的气溶胶观测数据对本文的反演结果进行验证[21]。除AERONET尚未提供下载的2018年香河站相关数据以外，期间共匹配194对有效数据。将AERONET得到的AOD转换到与反演一致的550 nm波段，以MODIS气溶胶产品的误差Δτ=±0.05±0.2τ[22]作为误差线，并与只使用暗像元算法的反演结果作进一步对比(图4、图5、表2)。

图4 暗像元-深蓝AOD与地基AOD对比图

图5 观测周期匹配散点变化对比图

表2 反演结果验证

由表2可见，暗像元-深蓝反演合成AOD与地基AOD表现出了较高好相关性，R2约为0.79，整体结果略高于地基AOD，当AOD>1时，监测精度走低，绝对误差和相关系数低于整体精度水平；AOD>1.2时，由图4可见，多数监测结果明显偏高而脱离了误差线范围。相比之下，暗像元算法的有效匹配数降低至89对，其中香河站仅有33对，且只集中在AOD较小区域(<1)，相对误差与相关系数明显偏低，说明了暗像元算法在冬季的应用存在较大局限性，合成AOD结果匹配点数量的明显增加表现了深蓝算法与暗像元算法在冬季的高互补性。

注：该图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2016)1666号的标准地图制作，底图无修改。图6 2018年2月26日AQUA/MODIS真彩影像

以云量较少的一次典型灰霾天气(2018年2月26日)为例对暗像元算法与深蓝算法监测效果进行进一步比较，从AQUA/MODIS真彩影像(图6)红线标记区域附近可明显看出京津冀中南部的大范围灰霾现象[23]。图7分别展示了该日暗像元、深蓝及合成AOD结果，可见深蓝算法对冬季植被稀疏地表尤其是灰霾区域AOD监测具有极大优势，而暗像元算法难以实现此类区域的有效监测。

注：该图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2016)1666号的标准地图制作，底图无修改。图7 使用不同算法的AOD反演结果(2018年2月26日)

3.2 与MODIS AOD产品对比验证

为了进一步验证本研究在监测京津冀地区AOD空间分布上的准确性和适用性，分别选择了2016—2018年3次较为典型的冬季灰霾天气(2016年1月10日、2017年2月13日、2018年2月26日)，利用MODIS MYD04_L2暗像元-深蓝AOD融合产品与本研究的暗像元-深蓝AOD反演结果进行对比验证。

由图8可见，本研究反演合成的AOD与MYD04_L2 AOD产品在空间分布上具有较高的一致性，AOD的高、低值区均表现出较为一致的分布趋势。利用ArcGIS生成随机点的方式，对2种AOD进行了统计分析，统计结果见表3。

注：该图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2016)1666号的标准地图制作，底图无修改。图8 AOD反演结果与MYD04_L2同期AOD对比图

针对2016年1月10日、2017年2月13日、2018年2月26日3次AOD监测，本研究监测结果在ArcGIS生成的随机点中，有效AOD点位占比分别达到69%、98%、80%，远好于MYD04_L2同期产品的22%、74%、62%；在有效监测匹配点的对比中，二者呈现出了较为显著的相关性，相关系数均达0.9以上。

表3 反演结果与MODIS产品对比验证

通过本研究与MODIS AOD产品在京津冀地区冬季的AOD监测结果对比分析中可得，本研究AOD监测结果的空间分布趋势一致性高，有效监测值的相关性强；监测盲区显著减少，适用性更高；且空间分辨率提高到了1 km，AOD空间分布的表现更为具体，整体监测效果优于MODIS AOD产品。

3.2 京津冀AOD的空间特征与影响因素分析

研究表明，轻度灰霾天气的AOD均值在0.8左右[24]。以0.8为阈值统计了京津冀地区2016—2018年1—2月高AOD天数分布情况(图9(a))。能够看出京津冀地区AOD的空间分布呈现较为明显的西北低、东南高，且分布图的均质性整体表现较好，能够将西北部的张家口、承德划为低AOD区，将京津唐及以南各市划为高AOD区，说明了京津冀地区的气溶胶分布与地势的显著负相关性。从图9(b)可见，西北部连续的高大山体形成了气溶胶扩散流通的天然屏障，对北方吹来的冷空气有着较强的阻挡和削弱作用，南部的山西、河南、山东省部分地区也向京津冀持续输送着大气污染物，在山谷风和城市热岛环流和海陆风环流的共同作用下，京津冀平原地区气溶胶的流通与扩散较为缓慢，输送过程较为复杂[25]。

注：该图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2016)1666号的标准地图制作，底图无修改。图9 京津冀AOD的空间特征分析图

燃煤取暖依然是人们过冬采暖的主要手段，采暖燃煤气溶胶也是冬季气溶胶的主要类型[26]。除北京、天津两大直辖市外，保定、石家庄、邯郸、唐山、沧州等东部南部城市都是河北省的人口大市，燃煤取暖气溶胶排放量相对于西北部人口稀疏的张家口、承德市明显更高。北京、天津、唐山等传统的工业城市和各类企业的聚集地也通常是AOD的偏高区域。表明了能源消耗、人类活动对气溶胶也分布也有着较大影响[27]。

4 结束语

本文发展了冬季气溶胶光学厚度遥感反演算法，该算法以暗像元算法和深蓝算法为基础，对近3年1—2月京津冀地区大气气溶胶进行了连续监测，通过与AERONET北京和香河地基站实测数据的对比验证，评价了算法的性能，并根据连续监测结果和辅助资料分析了冬季京津冀地区AOD的空间特征及其影响因素，主要结论如下：

本算法相较于传统的暗像元算法，有效观测数明显增加，且具备较为良好的监测精度，说明了深蓝算法与暗像元算的较强的互补性，是暗像元算法的有效补充，对2种算法结果进行融合处理在冬季是一种良好的监测手段。在与MYD04_L2气溶胶产品的对比中，二者对AOD的空间分布趋势及AOD的监测值均表现出了较高一致性，但本算法对冬季京津冀地区的监测结果拥有更大的有效范围，且空间分辨率提高到了1 km，监测效果更好。

当AOD较高时，本算法可能出现AOD高估现象而使监测结果超出Δτ=±0.05±0.2τ的误差要求，但相较于暗像元算法在冬季的监测表现[5]，本算法的监测效果在精度与适用性方面均有明显提高。

通过对连续监测结果的统计，分析了冬季京津冀地区AOD的空间分布情况。京津冀冬季AOD整体分布呈现较为明显的西北低、东南高，与地势呈显著负相关，较为复杂的地理条件使得华北平原气溶胶积累容易，扩散缓慢。此外，冬季的燃煤取暖等人为气溶胶排放活动也显著影响着AOD的分布。

致谢：MODIS数据产品由NASA提供，地基太阳光度计站点验证数据由AERONET提供，在此表示感谢！