基于TM图像的南京市气溶胶光学厚度反演

程 晨，陈 健，李鑫慧

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

0 引言

大气气溶胶是大气与悬浮在其中的固体和液体微粒共同组成的多相体系，是指悬浮在大气中的具有一定稳定性、沉降速度小的，尺度范围在10-3～102μm之间的分子团、液态或固态粒子所组成的混合物[1]。描述气溶胶性质的基本参数有光学厚度、尺度分布、折射率、散射相函数。在这些参数中，气溶胶光学厚度(aerosol optical depth，AOD)是各个特征参数的总体效应，集中体现了气溶胶的基本状况。传统的地面实时观测仍是研究AOD的重要手段[2-4]，但由于观测站点数量有限，且分布比较离散，无法反映整个区域的气溶胶空间分布情况。遥感技术连续、宏观、动态、快速的特点，为大范围的气溶胶观测和研究提供了可能。

20世纪90年代，Kaufman等利用绿色植被在红光波段和蓝光波段的反射率与2.1 μm处的反射率存在线性关系的这一性质，通过提取暗像元对陆地上的气溶胶光学厚度进行了反演[5];随后Robert等人针对“暗像元”的缺陷进行了改进，研究出了更加成熟的MODIS v5.2气溶胶产品业务化算法[6]。在国内，李成才等通过分析MODIS光学厚度遥感产品讨论了北京及周边地区的大气污染情况[7]。随着城市大气污染的监测需求日益增大，空间分辨率为30 m的Landsat TM数据在气溶胶监测中得到广泛应用。闵祥军[8]、宋巍巍[9]等利用暗像元法与 DEM模型结合的方法分别反演了敦煌地区和广州市的AOD，得到了较好的效果，但该方法对DEM的空间分辨率要求较高;刘小平[10]等通过改进暗像元获得AOD信息，从而实现了TM图像的快速大气校正，但是对于暗像元的选取要求较为苛刻;王耀庭[11]等利用两幅TM图像对比的结果对北京市的气溶胶空间分布进行了研究，但其前提是需要找到一幅“无污染”图像作为基准，使该方法的使用受到了限制。

本文以南京市为研究区，以Landsat5 TM图像为数据源，通过6S大气辐射传输模型建立查找表，在此基础上建立暗像元AOD的多元回归模型，最后利用克里金插值得到南京市的气溶胶光学厚度的空间分布数据，并对南京市AOD的空间分布特征进行了分析。

1 研究区概况和数据源

1.1 研究区概况

以整个南京市地区为研究区。南京市位于N31°14′～32°37′，E118°22′～119°14′之间。图 1 为南京市的TM图像。南京地区属于亚热带季风气候，雨量充沛，四季分明。全市拥有常住人口810万，南京市是一座人类活动密集、偏重于工业生产的大城市，每天会排放出大量的废气，特别是机动车排放出的尾气，矿产采集和化工生产的过程中扬起的粉尘都会进入大气中，成为气溶胶的一部分，严重影响当地的空气质量，以致影响居民的正常生活。

图1 研究区TM图像Fig.1 TM image of the study area

1.2 遥感数据与预处理

本研究选用的是一幅2006年5月20日的Landsat5 TM图像，轨道号120-38，从国际科学数据服务平台上下载获得。当天天气晴朗，少云。首先，通过已知的控制点对图像进行几何纠正，重采样方式为最邻近插值法;然后按照南京市的行政区划裁切出研究区范围;最后进行辐射定标，将图像中每个像元对应的DN值转换为辐亮度，进行AOD反演，详细计算过程和主要参数参见文献[12]。

2 气溶胶光学厚度反演原理与方法

在地表为朗伯体、大气水平均一的假设条件下，大气顶部的表观反射率[5]可以表达为

式中:ρ*为表观反射率(卫星观测到的反射率);ρa为大气的路径辐射对表观反射率的贡献;μ，φ，μ0，φ0分别为太阳天顶角的余弦、方位角和观测天顶角的余弦、方位角;T(μ0)和T(μ)分别为上行辐射和下行辐射的透过率;s为大气的球面反照率;ρ为地面的真实反射率。在实际反演AOD时，往往利用辐射传输模型计算出 AOD 与 s，ρa，T(μ0)和 T(μ)等参数的对应关系，据此建立查找表，通过查找表获取气溶胶光学厚度。本文的技术路线如图2所示。

图2 反演气溶胶光学厚度流程Fig.2 Flow chart of AOD retrieveing

2.1 计算表观反射率和地表反射率

图像在经过预处理后，通过计算将暗像元的表观辐亮度转化为表观反射率，即

式中:ρ为表观反射率;L为表观辐亮度，W·m-2·μm-1sr-1;D为日地之间距离(天文单位);ESUN为大气层顶的平均太阳辐照度，W·m-2;θ为太阳天顶角。

根据Kaufman的研究，在洁净大气的条件下，绿色植被在可见光的红、蓝波段和2.1 μm红外波段的反射率在图像上有相似的分布情况，进而推测红、蓝波段和2.1 μm的中红外波段可能存在着一定的线性关系，而这种线性关系对于红外波段波长为2.08 ～2.35 μm 的 TM 图像同样适用[13]。所以，如果TM图像上植被覆盖区在可见光的红、蓝波段的地表反射率可以通过红外通道的反射率估算出来，那么图像上得到的红、蓝波段的表观反射率和真实地表反射率之间的差异就是气溶胶造成的影响。TM图像的红光和蓝光波段地表反射率由以下经验公式[13]获得，即

式中，Rred，Rblue分别为TM图像红、蓝波段的地表反射率;MIR为TM图像中红外波段的表观反射率。

2.2 确定暗像元

在红、蓝波段遥感图像上，密集的植被相对于其他地物背景较暗，所以被称为“暗像元”。本文采用阈值法，通过设置NDVI和MIR的阈值来提取暗像元[14]，当像元满足 NDVI＞0.2，MIR ＜0.1 时，该像元被认定为暗像元。暗像元提取结果见图3。由图上可以看出，暗像元基本呈现均匀分布。

图3 暗像元提取结果Fig.3 Results of dark pixel extraction

2.3 确定6S辐射传输模型及其参数

6S辐射传输模型是一个被广泛接受和使用的大气辐射传输模型，用来模拟地气系统中太阳辐射的传输过程。其输入参数主要包括太阳和卫星几何路径参数、大气模式、气溶胶模式、光谱条件和地面反射率类型。

TM图像的太阳和卫星几何路径参数及光谱条件由6S模型提供;大气模式选择为中纬度夏季大气(确定大气温、压、湿廓线、臭氧、水汽含量);下垫面假定为均一表面，无方向性反射。而根据MODIS AOD业务化算法[15]，本文研究区域被划分为中等吸收气溶胶类型(单次散射反照)。陈骁强计算AERONET太湖站点2005年9月—2009年10月月平均单次散射反照率的平均值后发现，太湖站的单次散射反照率平均值为0.907 5，与大陆型气溶胶中等吸收型气溶胶对应的单次散射反照率值最为接近[16]，所以本文在6S辐射传输模型中气溶胶类型选择为大陆型气溶胶。

2.4 构建查找表

本文借助6S辐射传输模型来模拟TM图像红、蓝波段在一定的气溶胶光学厚度和气溶胶模式的条件下，表观反射率随AOD和地表反射率变化而发生的变化，从而建立查找表。建立查找表时，由于其他参数已经确定，主要考虑AOD和表观反射率的变化。通过对TM1和TM3这2个波段的像元统计发现，超过99.9%的像元反射率分别集中于0.09～0.15和0.06 ～0.13 区间内，所以在查找表中 TM1和TM3波段的表观反射率分别设定为0.09～0.15和0.06～0.13(都是以0.01为步长)。由于晴空下一般AOD＜2.0，且6S模型在气溶胶光学厚度过大(能见度小于5.0 km)时会出现程序错误，所以TM1和TM3波段的气溶胶光学厚度设定为0.21～1.21(以0.1为步长)，对于超出范围的采用外推计算，建立起了利用TM数据反演AOD的查找表。

2.5 多元回归建模

在得到TM1和TM3波段关于气溶胶光学厚度、地表反射率、表观反射率3者关系的查找表后，利用SPSS软件分别对两组数据进行多元线性回归分析。TM1和TM3波段的回归方程分别为

式中:z为气溶胶光学厚度;x为表观反射率;y为地表反射率。

2.6 反演气溶胶光学厚度

通过查找表分别得到TM1和TM3波段在550 nm上的暗像元AOD值以后，求平均获得南京市AOD空间分布数据。在ArcGIS中进行插值，为了减少计算量，在插值前将研究区TM图像重采样成60 m分辨率，重采样方式为最邻近像元法。由于在不均匀离散点分布的插值中，克里金插值比其他插值方法具有明显的优越性[17]，所以本研究选用克里金插值法。得到的结果如图4所示。

图4 2006年5月20日南京市气溶胶光学厚度空间分布Fig.4 AOD spatial distribution in Nanjing City，May 20，2006

3 结果与分析

3.1 AOD反演结果合理性检验

赵英时等认为，气溶胶散射会对TM图像的可见光波段产生较大的影响[18]，通过对比大气校正前后南京市真彩色合成图像和AOD反演结果，可以定性地检验AOD反演结果的合理性。因此，对TM数据进行了FLAASH大气校正。在大气校正前后的图像上各选取茂密植被区和裸地区，获取其随波长变化而发生变化的光谱曲线，分析大气校正前后地物光谱的变化。如图5所示。

图5 大气校正前(左)后(右)植被和裸地光谱反射率曲线Fig.5 Spectral signature of bare soil，vegetation before(left)and after(right)atmospheric correction

从图5可以看到，在可见光波段植被和裸地的表观反射率在校正前后发生了明显变化，其变化符合实际地物反射率变化规律，这说明大气校正取得了较好的效果。

图6分别为大气校正前后的真彩色合成图像以及把AOD反演结果取代校正后图像红波段进行彩色合成后的图像。

图6 南京市大气校正前(左)、后(中)图像与AOD反演结果(右)的比较Fig.6 Comparison of atmospheric correction before(left)and after(middle)and AOD inversion results(right)of Nanjing City

对比图6(左)和图6(中)可以看出，大气校正后南京市城区大部分气溶胶被去除，图像对比明显增强，地物的细节更加清晰(图6(中));对比图6(右)和图6(左)则可以发现，AOD空间分布与图6(左)中的气溶胶(相对模糊的区域)的强弱分布基本一致。综上所述，本研究中南京市当日的AOD反演结果较合理。

3.2 AOD空间分布特征分析

从图4中可以发现，南京市AOD分布呈现北部相对较高，而中部和南部较低的基本空间分布特征，尤其在长江两岸地势平坦、人类生活频繁、制造业发达的地区甚至出现了超过1.0的高值。在市区范围，除了燕子矶、紫金山周围茂密的植被覆盖地区AOD较低以外，其他地区的AOD都在0.7～0.9之间，处于相对较高的水平，这很可能跟汽车尾气排放等较密集的人类活动有关。从市区往南是江宁区、禄口镇、溧水县等区域，可以发现这一带是整个南京市AOD整体最低的一个区域，主要原因是这些区县相对市区人口较少，植被茂密。但是有小部分区域也出现了超过0.6的相对高值，如溧水县爱景山地区。这是因为该地区有全国最大的锶矿在开采，采矿过程中必然会产生很多粉尘，当粉尘散布到大气中时就会对气溶胶光学厚度产生显著影响。根据上述描述可以作出初步判断，植被、城市建成区、地形是影响AOD分布的主要因素。这个认识与李成才[19]、段婧[20]等的研究结果一致。

为了深入研究植被、城市建成区、地形对AOD空间分布的影响，本文引入了3个辅助数据，即归一化植被指数(normalized difference vegetation index，NDVI)、归一化城镇指数(normalized difference building index，NDBI)[21]和数字高程模型 (digital elevation model，DEM)，将 AOD 以每 0.01 步长进行分割，然后分别计算每0.01步长中所有像元的NDVI，NDBI，DEM的平均值。计算结果如图7所示。这里需要说明的是，在提取暗像元时，由于提取条件的限制，水体区并没有暗像元提出，水体上空的AOD由周围区域暗像元AOD插值得来，所以为了不影响准确性，在统计的过程中剔除了水体。

图7 AOD与NDVI(左)、NDBI(中)、DEM(右)间的变化规则Fig.7 Scatter plot of AOD corresponding with the average of NDVI(left)，NDBI(middle)and DEM(right)

从图上可以发现，AOD 与 NDVI，NDBI，DEM 都存在着良好的相关性，随着AOD的逐渐升高，3者都呈现了非常明显的变化。

从图7(左)中可以看出，AOD与NDVI呈现明显的负相关关系。随着AOD的增加，NDVI值相应迅速降低，然后趋于平缓，最后出现小幅上升。AOD低值的区域主要对应了NDVI较高的茂密森林区域，这说明植被茂密的程度与气溶胶的空间分布存在紧密的联系。当AOD处于高值时，对应的NDVI值变化不大，甚至有小幅的上升，这有可能是因为在一些AOD高值区域，气溶胶颗粒的来源不再跟下垫面类型有紧密的关系，而是主要来自于周围区域的扩散，这也可能就是南京市城区一些植被覆盖很好的公园、绿地AOD偏高的原因。

从图7(中)中可以看出，AOD与NDBI呈现明显的正相关关系。随着AOD的增加，NDBI值迅速升高，在0.7～0.8处出现了震荡，在趋于平缓后又出现小幅下降。这说明城镇建筑的集中程度与气溶胶的空间分布也存在紧密的联系。在AOD为0.7～0.8时出现的震荡，经过观察发现六合区东南部AOD相对较高，而这部分区域下垫面主要为农田，NDBI处于相对低值。造成这样情况的原因可能是这里地势平坦，向西是处于AOD高值区的大厂区域，受到扩散的作用;而南部是长江，水体周围的区域会受到水体蒸发产生液态悬浮性气溶胶颗粒的影响。当AOD处于高值时，NDBI值变化不大，甚至有小幅的下降，这可能由于在南京城市规划中，化工厂、发电厂等能够产生大量大气污染物的区域并不是位于市区中心区域，而是位于大厂镇等城市边缘乡镇。可见，造成一座城市AOD最高的区域常常不是人口和商业活动最集中的区域。

从图7(右)中可以看出，随着地形的升高，AOD呈现负指数衰减，这与徐希孺[22]等的观点一致。图7中AOD的低值对应的主要是南京市的一些地势较高的区域，例如紫金山、老山等海拔超过200 m的区域，而AOD的高值区主要位于海拔50 m以下的区域，这主要是由于地球重力作用，气溶胶颗粒密度随高度呈指数衰减[22]。对于地形较平坦的南京来说，地形的微小变化无法对气溶胶的扩散产生影响，但边界层的湍流运动导致城区气溶胶颗粒向四周大范围的扩散，这可能就是江宁北部和六合南部地区虽然工业较少，依然很难看到蓝天的缘故。

综上所述，利用TM图像获取60 m空间分辨率的AOD空间分布信息，对于在城市尺度上AOD的定量化研究会有所帮助，并且对于类似南京市这样的空间污染相对严重且分布具备明显规律的城市非常有应用价值。

4 结论与展望

本文利用TM数据，对南京市2006年5月20日的气溶胶光学厚度进行了卫星遥感反演，得到了当天南京市气溶胶光学厚度的空间分布图。主要获得以下结论:

1)基于TM图像，通过暗像元法和辐射传输模型建立查找表反演气溶胶光学厚度的方法是可行的，达到了气溶胶光学厚度精细化反演的目的。在提取足够的浓密植被作为暗像元过程中，针对南京市的植被覆盖情况和TM图像的特点，NDVI＞0.2，MIR＜0.1的阈值设定是有效的。

2)南京市的AOD与植被覆盖呈现明显的负相关，与城市建筑密集程度呈现明显的正相关，与地形高低呈现明显的负指数相关。气溶胶主要聚集于植被覆盖较低、城市建筑密集、人类活动较频繁、地形较低的区域。高空间分辨率图像比一般分辨率图像能够提供的更多的气溶胶空间分布信息。

结合辐射传输模型进行AOD反演的暗像元法针对植被覆盖条件较好的中国东南地区是可行的，但是也正是需要浓密植被覆盖这一条件局限了这一方法在其他植被覆盖稀少地区的应用。如果能够在以后的研究中减少对于浓密植被像元的依赖，并且更好地解决水体上空暗像元不足的问题，那么城市AOD反演的精细化研究一定可以得到进一步的发展。

[1]章澄昌，周文贤.大气气溶胶教程[M].北京:气象出版社，1995:2.Zhang C C，Zhou W X.Atmospheric aerosol tutorial[M].Beijing:China Meteorological Press，1995:2.

[2]夏祥鳌，王普才，陈洪滨，等.中国北方地区春季气溶胶光学特性地基遥感研究[J].遥感学报，2005，9(4):429-437.Xia X A，Wang P C，Chen H B，et al.Ground-based remote sensing of aerosol optical properties over North China in spring[J].Journal of Remote Sensing，2005，9(4):429-437.

[3]赵柏林，王 强，毛节泰，等.光学遥感大气气溶胶和水汽的研究[J].中国科学:B 辑，1983，13(10):951-962.Zhao B L，Wang Q，Mao J T，et al.Study on remote sensing of aerosol optical depth and vapor[J].Scientia Sinica Chimica:Serier B，1983，13(10):951-962.

[4]毛节泰，王 强，赵柏林.大气透明度光谱和浑浊度的观测[J].气象学报，1983，41(3):322-331.Mao J T，Wang Q，Zhao B L.The observation of the atmospheric transparency spectrum and the turbidity[J].Acta Meteorologica Sinica，1983，41(4):322-331.

[5]Kaufman Y J，Tanré D，Remer L A，et al.Operational remote sensing of tropospheric aerosol over land from EOS moderate resolution imaging spectroradiometer[J].Journal of Geophysical Research，1997，102(D14):17051-17067.

[6]Levy R C，Remer L A，Matto S，et al.A new algorithm for retrieving aerosol properties over land from MODIS spectral reflectance[J].Journal of Geophysical Research，2006.

[7]李成才，毛节泰，刘启汉，等.利用MODIS光学厚度遥感产品研究北京及周边地区的大气污染[J].大气科学，2003，27(5):869-880.Li C C，Mao J T，Liu Q H，et al.Research on the air pollution in Beijing and its surroundings with MODIS AOD products[J].Chinese Journal of Atmospheric Sciences，2003，5(27):869-880.

[8]宋巍巍，管东生.利用TM影像反演广州市气溶胶光学厚度空间分布[J].环境科学学报，2008，28(8):1638-1645.Song W W，Guan D S.The distribution of aerosol optical depth retrieved by TM imagery over Guangzhou，China[J].Acta Scientiae Circumstantiae，28(8):1638-1645.

[9]闵祥军，朱永豪.基于Landsat-TM图像自身的反射率反演方法[J].遥感技术与应用，1997，12(3):2-9.Min X J，Zhu Y H.A method of reflectance retrieval based on image of Landsat-5 TM[J].Remote Sensing Technology and Application，1997，12(3):2-9.

[10]刘小平，邓孺孺，彭晓鹃.基于TM影像的快速大气校正方法[J].地理科学，2005，25(1):88-92.Liu X P，Deng R R，Peng X J.A fast atmospheric correction method based on TM imagery[J].Scientia Geographica Sinica，2005，25(1):88-92.

[11]王耀庭，王 桥，杨一鹏，等.利用Landsat/TM影像监测北京地区气溶胶的空间分布[J].地理与地理信息科学，2005，21(3):20-22.Wang Y T，Wang Q，Yang Y P，et al.Monitoring spatial distribution of aerosols in Beijing based on Landsat TM Data[J].Geography and Geo-Information Science，2005，21(3):20-22.

[12]Chander G，Markham B.Revised Landsat 5 TM radiometric calibration procedures and post calibration dynamic ranges[J].IEEE Transactions on Geoscience and Remote Sensing，2003，41(11):2674-2677.

[13]Kaufman Y J，Wald A E，Remer L A，et al.The MODIS 2.1- μm channel-correlation with visible reflectance for use in remote sensing of aerosol[J].IEEE Transactions on Geoscience and Remote Sensing，1997，35(5):1286-1298.

[14]Guo Y W，Tsay S C，Cahalan R F，et al.Path radiance technique for retrieving aerosol optical thickness over land[J].Journal of Geophysical Research，1999，104(D24):31321-31332.

[15]Remer L A，Tanré D，Kaufman Y J，et al.Algorithm for remote sensing of tropospheric aerosol from MODIS:Collection 5[R].National Aeronautics and Space Administration，2006.

[16]陈骁强.基于环境一号卫星陆地气溶胶光学厚度反演方法研究[D].南京:南京师范大学，2011.Chen X Q.Study on the retrieval of aerosol optical thickness over land based on HJ-1 data[D].Nanjing:Nanjing Normal University，2011.

[17]刘光孟，汪云甲，张海荣，等.空间分析中几种插值方法的比较研究[J].地理信息世界，2011，9(3):41-45.Liu G M，Wang Y J，Zhang H R，et al.Comparative study of several interpolation methods on spatial analysis[J].Geomatics World，2011，9(3):41-45.

[18]赵英时.遥感应用分析原理与方法[M].北京:科学出版社，2003:21-30.Zhao Y S.The principle and method of analysis of remote sensing application[M].Beijing:Science Press，2003:21-30.

[19]李成才，毛节泰，刘启汉.利用MODIS资料遥感香港地区高分辨率气溶胶光学厚度[J].大气科学，2005，29(3):335-342.Li C C，Mao J T，Liu Q H.Remote sensing of high spatial resolution aerosol optical depth with MODIS data over Hong Kong[J].Chinese Journal of Atmospheric Sciences，2005，29(3):335-342.

[20]段 婧，毛节泰.长江三角洲大气气溶胶光学厚度分布和变化趋势研究[J].环境科学学报，2007，27(4):537-543.Duan J，Mao J T.Study on the distribution and variation trends of atmospheric aerosol optical depth over the Yangtze River Delta[J].Acta Scientiae Circumstantiae，2007，27(4):537-543.

[21]查 勇，倪绍祥，杨 山.一种利用TM图像自动提取城镇用地信息的有效方法[J].遥感学报，2003，7(1):38-39.Cha Y，Ni S X，Yang S.An effective approach to automatically extract urban Land-use from TM imagery[J].Journal of Remote Sensing，2003，7(1):38-39.

[22]徐希孺.遥感物理[M].北京:北京大学出版社，2003:198-312.Xu X R.Physics of remote sensing[M].Beijing:Peking University Press，2003:198-312.