从短波红外与红光波段反演华北地区气溶胶

张钰萌，陈辉，马鹏飞，厉青，王中挺

(1.河南理工大学 测绘与国土信息工程学院，河南 焦作 454000；2.生态环境部卫星环境应用中心，北京 100094)

0 引言

华北地区是指位于中国北部的区域，包括北京、天津、河北、山西和内蒙古中部。该地区面积约84×104km2，人口约1.68 亿，2015年GDP为10 万亿元。该区域聚集有大量钢铁、化工、水泥等重污染产业，大量排放着氮氧化物、粉尘等大气污染物，而普遍缺少治污措施的“散乱污”企业也在排放大气污染物，同时，华北地区机动车总数超过了2 600万，成为大气污染物的主要来源之一。研究表明[1]，较为严重的大气污染，使得我国北方的人均预期寿命将减少5.5 年，同时也大大提高了肺癌、心脏病等重大疾病的发病率。

利用遥感数据获取华北地区气溶胶的空间分布，不仅能获得气溶胶污染的空间分布状况，而且能为颗粒物指标的反演提供数据源[2]。针对陆地地表，由于卫星接收到的信号中地表反射较强，需要去除地表反射率贡献得到大气信息，从而反演得到气溶胶，目前主要的算法有：假定浓密植被在红蓝波段存在线性关系的暗目标法[3-5]、假定地表反射随时间变化较慢的深蓝算法[6-8]、利用气溶胶散射的角度变化的多角度算法[9-11]、采用偏振信号的多角度偏振法[12-15]。在暗目标法基础上，可以将红蓝等波段存的线性关系扩展到低植被覆盖的地表，从而实现多种地表覆盖类型的气溶胶反演[16-17]。

搭载于美国的Terra以及Aqua星的MODIS传感器已持续获取数据近20 年，拥有覆盖可见、近红外、短波红外、热红外等36个波段，目前该数据已广泛应用于气象、环境、农业等方面。据不完全统计，我国已在气象、环境、农业等多个领域建立有数十个能直接接收MODIS数据的直收站。

然而，在利用MODIS数据进行华北地区气溶胶反演时，还存在着以下两个问题：一是冬季时，我国北方地区植被稀疏，暗目标较少，采用暗目标法较为困难，且精度偏低[18]；二是MODIS传感器已超期服役，器件开始老化，直收站接收到的数据在蓝光波段存在着明显的条带。因此，本文在Wang等[17]提出的算法基础上，结合短波红外波段从MODIS数据监测华北地区秋陆地气溶胶。

1 基本原理与算法

假定大气层中的气溶胶、大气分子水平分层分布，传感器接收到的反射率信号ρTOA为[19]：

(1)

式中：S、ρ0和T是3个代表大气状况参数，去除大气分子影响即可反演气溶胶光学厚度(AOD)。

除水体与雪外，大部分地表在短波红外波段(2.12 μm)的反射率与红波段(0.66 μm)反射率存在着线性关系[20]：

rSWIR=a×rred+b

(2)

式中：a、b为经验参数，根据随着地物类型而变化。Kaufman等[21]认为，b是大气校正中对气溶胶影响校正不足或校正过多带来的，因此，本研究中将b统一设置为0。

结合式(1)，可以得到方程组：

(3)

解上述方程组反演得到AOD。具体为：使用辐射传输软件计算各个AOD条件下短波红外波段和红波段的大气参数S、ρ0和T；利用式(3)和式(4)，将MODIS观测到的短波红外波段表观反射率ρTOA，SWIR和红波段ρTOA，red代入，得到rSWIR和rred；将rred代入式(5)，计算得到rSWIR；计算利用式(3)和式(5)得到的rSWIR之间的差，以差为自变量，AOD为因变量，进行线性拟合，求取差为0时的AOD值，即为所求。

在实际反演中，为加快运算速度，本文预先使用6S vector1.0[22]设定AOD从0～4变化，进行辐射传输运算存储在查找表中。根据MODIS的过境时间和观测特征，设置太阳天顶角为0～60°，观测天顶角0～60°，相对方位角为0～180°。气溶胶类型的设定参照了地面观测的结果[23-25]。华北地区的气溶胶由3种基本气溶胶混合而成，以水溶型为主，沙尘型和煤烟型为辅，设定各气溶胶组分比例为：沙尘型15%、水溶型83%、黑炭型2%。

2 地表反射率分析

为分析华北地区各地表类型的短波红外波段的反射率与红波段反射率的情况，本文收集了2015年MODIS的地表反射率产品，该产品校正了气溶胶、大气分子等对卫星观测反射率的影响[26]。提取2015年京津冀区域的红光、近红外和短波红外波段数据共计46景。然后，计算每景数据的红光和短波红外波段的比率和归一化植被指数(NDVI)。为消除云、雪等因素的影响，去除了地表反射率大于0.2以及NDVI小于0的数据。

随机选取了2015年度所有时间的城市、农田、森林、草地共4种典型土地类型的反射率比率、NDVI、红光波段反射率，分析它们随时间变化的情况。

图1表明红光和短波红外地表反射率比率整体较为稳定，随着时间变化较小，城市最高，在0.8左右；农田次之，在0.7左右，7—9月有明显的低值出现；草地更低，为0.6左右，随时间变化最小；森林最小，在0.4左右。

图1 城市、农田、森林、草地反射率比率随时间变化图

从图2可以看出，各典型地物的NDVI有较强季节变化规律，夏季植物茂盛达到最高，冬季植被凋零降到最低。森林整体最高，5—10月最高，在0.6以上，其他时间在0.4以下；农田次之，7—10月最高，在0.6以上，4—6月在0.4左右，其他在0.2左右；草地和城市较为接近，6—10月较高，在0.2以上。

图2 2015年度的城市、农田、森林、草地NDVI

从图3可以看出，各典型地物在红光波段的反射率也体现了季节变化规律。草地整体最高，在0.148左右；农田次之，7—9月最低，3—5月也较低；城市更低，7—9月达到最低；森林最低，2—4月最高，8—10月最低。

图3 2015年度的城市、农田、森林、草地红波段反射率

表1为各典型地物反射率比值、NDVI、红波段反射率统计结果。本文用标准差占平均值的百分比来衡量不同指标随季节变化的离散程度。可以看出，反射率比率的离散程度最小，百分比在15%以下；NDVI离散程度最大，百分比都在30%以上。因此，可以将红光和短波红外波段反射率比率的年平均值作为经验系数实现陆地气溶胶的反演。

表1 各典型地物反射率比值、NDVI、红波段反射率统计表

因此，本文将2015年的红波段与短波红外波段地表反射率比率作年平均，得到华北地区的反射率比率分布图，如图4所示。可以看出，华北地区的反射率比率随地物类型有明显的不同，河北北部的燕山、山西与河北交界的太行山等以森林覆盖为主等山区较低，在0.6以下，占据华北平原大部分地区的农田区域则在0.6～0.8之间，部分城市密集区以及河流高于0.8，内蒙古的草原区域在0.6左右。

注：该图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2017)1268号的标准地图制作，底图无修改。图4 华北地区年均反射率比率分布图

根据图4，本文统计了红波段与短波红外波段地表反射率比值出现的频率，如图5所示。可以看出，主要分布在0.4～0.7之间，其中，值在0.5和0.6出现的频率最多，超过了40%。

图5 华北地区年均反射率比值分布图

3 结果与验证

本文将上述分析得到的红波段与短波红外波段地表反射率比率年均值作为经验系数代入公式(3)至公式(5)，从而完成气溶胶的反演。具体的反演过程采用IDL语言实现。

3.1 反演实验

本文选取了2016年9月—2017年8月期间过境京津冀地区的MODIS数据进行了算法反演实验。图6为2016年10月13日、2017年1月17日、2017年3月19日、2017年5月4日4天比较典型的气溶胶状况反演获得的反演结果。由于可见红外波段无法穿透云层，本算法在云覆盖区域失效；同时，在冰雪覆盖区域，由于地表反射率过高，本算法也无法获得有效结果。

注：该图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2017)1268号的标准地图制作，底图无修改。图6 本算法的AOD反演结果图

如图6所示，2016年10月13日，华北地区的空气状况整体较好，仅在北京东南部、天津和河北中部出现了较高的气溶胶，部分地区高于1。2017年1月17日，华北地区的空气状况东南部出现了污染，较高的气溶胶主要分布在河北南部、河南北部、山东西北部，在河北河南交界处出现大片AOD大于2的区域，在河北北部有部分地区没有反演结果，可能是该区域属于山区，被雪覆盖后没有融化，雪的反射率过高导致算法失效。2017年3月19日，华北地区大部地区都出现了较重的空气污染，整体AOD在1以上，河北中南部AOD较高，在2以上，河北南部和河南由于被云覆盖没有获得反演结果。2017年5月4日，华北地区大部分地区受到沙尘天气，整体上环境大气质量较差，在河北中部出现了污染聚集团，AOD在2以上，但在中间最重区域算法失效，存在部分空白区域，原因是高浓度气溶胶与云的反射特征较为相近，本文将该处误识别为云。

3.2 结果验证

本文收集了卫星过境时AERONET北京Radi站和香河站气溶胶产品以验证本算法精度。北京Radi站位于北京市区北部，周围有大量居民点和商业区，属于典型的城市区域，香河站位于北京东南的香河县，周围以农田为主。AERONET是对气溶胶特性进行地面观测的网络，提供全球的气溶胶光学厚度、谱分布和散射相函数[27]观测结果。AERONET一般提供440 nm、675 nm等波段的气溶胶光学厚度，本文采用Angstrom公式[28]将AOD转化为550 nm结果，与本文结果相比较。

去除云、雾、雨、雪等天气影响的数据，在2016年9月1日至2017年8月31日TERRA和AQUA卫星过境的365 d中，本算法在北京Radi站共178 d获得有效验证数据，本算法在香河站共191 d获得有效验证数据(图7、图8)。

图7 本文结果与AERONET气溶胶产品时间变化(香河站)

图8 本文结果与AERONET气溶胶产品时间变化(北京Radi站)

从图7和图8可以看出，本算法获得的结果与地面结果较好地贴合在一起，但在AOD较高时部分反演结果与地面观测的结果有较大差异。而时间序列的观测结果也部分体现了季节变化的特征，香河站5—8月AOD明显低于其他月份，北京站6—8月的AOD也低于其他季节，但中间存在若干高AOD日。

图9 本文反演AOD结果与AERONET散点图(香河站)

图10 本文反演AOD结果与AERONET散点图(北京Radi站)

表2 反演结果验证表

注：r1为误差线范围内反演结果的比例；r2为比地面结果大同时不在误差线范围内反演结果的比例；r3为比地面结果小同时不在误差线范围内反演结果的比例。

3.3 分季节验证

本文分春、夏、秋、冬4个季节开展了算法验证，结果见图11、图12和表3。从图中可以看出，本算法在各个季节获得的结果都较为理想，相关系数都在0.9以上，RMSE也基本在0.1左右，在算法整体低估的情况下，秋季低估的状况最为明显；其中，香河站春季效果最好，夏季效果最差；北京站秋季和冬季效果较好，春季效果最差。可能是本文所用的气溶胶模型，在春季时，跟香河站的情况较为一致，但低估了北京地区的沙尘含量。

图11 本文反演AOD与AERONET产品之间分季节比较散点图(香河站)

图12 本文反演AOD与AERONET产品之间分季节比较散点图(北京Radi站)

表3 反演结果分季节验证表

4 结束语

本文探讨了通过假定短波红外波段与红光波段的地表反射率之间比率不变来获取华北地区陆地上空的气溶胶，主要取得如下结论：①与NDVI、红波段地表反射率相比，短波红外波段与红波段的地表反射率的比值较为稳定，利用该比值能够去除地表影响反演得到陆地气溶胶。②利用本文算法能够较好地获得华北地区城市和乡村区域的气溶胶分布，地面验证表明70%的结果误差范围在±0.05±0.15τ之内，且与地面结果的相关系数在0.9以上。

但本算法还需要在2个方面做进一步改进：①气溶胶各组分的比例固定，虽然整体也取得了不错的反演效果，但随着地表类型、季节的变化，会对反演结果带来±10%左右的偏差(图10、图11)，在今后的反演中，还需借助紫外、偏振等多源卫星数据提取气溶胶组分信息以提高反演精度。②本文算法将每个像元的红波段与短波红外波段地表反射率比率固定为其年均值，但还存在着一定的变化，统计表明(表1)，其与年均值偏移量在15%左右，考虑到红波段的反射在0.1左右，反射率比值在0.5～0.6，则其对地表反射率的误差在0.01左右，会造成约0.1气溶胶光学厚度误差[3]。

致谢：本文的MOD09数据由美国国家航空航天局提供，AERONET北京Radi站数据由中科院遥地所提供，AERONET香河站数据由中科院大气所提供，6S Vector1.1程序由Vermote E等提供，在此表示感谢。