大气气溶胶光学厚度反演软件系统设计和实现

刘翔龙，杨再昕

（河海大学 计算机与信息学院，江苏 南京 211100）

大气气溶胶是由大气介质和混合于其中的固体和液体颗粒物等组成的体系，大气气溶胶主要通过直接或间接改变地气系统的辐射收支来影响气候和环境。由于大气气溶胶光学厚度AOD（Aerosol Optical Depth）是反映大气浑浊度和确定气溶胶气候效应的重要依据，所以对气溶胶光学厚度的研究具有很现实的意义。在气溶胶光学厚度反演这一研究课题上国内外专家和学者已经做了很多工作。当前的气溶胶光学厚度主要是以地基监测为主，但是地基监测方法会受到很多方面因素的影响，例如气溶胶的生命周期较短，浓度的空间变化很大等。这样只能监测整层大气柱的气溶胶的总量和单点的信息，同时由于各地区地理环境条件对仪器监测有所限制，所以只能在有限的区域进行，所设的网点分布不能全部覆盖和反映各地信息。卫星遥感反演方法与传统的地基探测方法不同，它具有覆盖面积广、信息获取方便等特点，能更有效地获取气溶胶的信息，摆脱了地基探测局限性的缺点，使人们能够实时掌握大面积范围内的气溶胶的变化提供了可能。目前，在使用卫星遥感资料进行气溶胶光学厚度的反演是借助不同软件的部分功能模块加以整合得以实现的，而专门用来实现气溶胶光学厚度反演的一体化的软件几乎没有，这样就使研究者需要对反演方法的全面掌握和对不同软件的各部分功能模块有所了解，增加了不必要的工作，使反演过程复杂不便[1-2]。

为了实现对反演软件的开发，根据目前已有反演的理论基础，主要提出和设计利用GDAL开源库对MODIS L1B遥感数据提取和进行相关处理，选取暗像元后利用其参数信息结合6S模型建立查找表。使用QT设计用户操作界面，实现反演可视化。本软件能够实现气溶胶光学厚度的反演的功能。

1 反演原理

一般来说，当陆地表面是均匀朗伯表面，大气垂直均匀变化时，卫星测量值可用表观反射率ρ*表示：

式中L是卫星测量辐亮度，Es是大气顶的太阳辐射通量，τ0是整层大气光学厚度，（μs，Φs）表示太阳入射光的方向，（μv，Φv）表示卫星观测方向，μs，Φs，μv，Φv分别为太阳入射光和卫星观测天顶角的余弦和方位角。如果不考虑气体吸收，卫星观测的表观反射率为：

式中θs为太阳天顶角，θv为观测天顶角，Φ为相对方位角，Φs分别为太阳到地面和地面到卫星的整层大气 T（θv）透过率，S 为大气的球面反照率，ρ为地表反射率，Fd（θs）为表观反射率为零时归一化总的向下辐射通量。当地表反射率很小时，卫星观测反射率主要取决于大气贡献项。当地表反射率很大时，地表的贡献成为主要贡献项。这样可以得出反演的基础：如果已知地表反射率，确定大气气溶胶模型就可以得到气溶胶光学厚度。反之，如果已知气溶胶光学厚度和相应的大气参数，也可以得到地表反射率。

该软件的设计主要依据暗像元方法DDV（Dense Dark Vegetation），它是由Kaufman和Sendra反演稠密植被上空气溶胶光学厚度建立的。这种算法主要依据浓密植被和湖泊的反射率很低（约0.01～0.02），因此可以将有着大面积森林或水域的地方可以作为暗像元，暗像元算法是基于表观反射率的大气贡献项，即利用卫星观测的路径辐射反演气溶胶光学厚度，是气溶胶遥感应用比较常用的算法。利用MODIS图像反演气溶胶光学厚度主要是利用了空间分辨率为250 m的1波段（620～670 nm）、500 m 的 3 波段（459～479 nm）、500 m 的7波段（2 105～2 155 nm）3个波段。通过大量实验证明：中红外通道（2 100 nm）通道表观反射率除了受尘粒气溶胶的影响外，几乎不受其他气溶胶的影响而且其值接近地表反射率。因此暗像元的确定主要是在几何校正后通过第7波段在MODIS图像中寻找表观反射率小于0.05的点作为暗像元。

Kaufman等人通过研究大量的资料，得到红（670 nm）、蓝（470 nm）和中红外通道（2 100 nm）的地表反射率存在如下关系：

利用中红外通道的表现反射率和地表反射率近似相等的前提，可以采用7波段的表观反射率代替地表反射率，代入式子（3）、（4）估算求出1、3波段的地表反射率。并且选取合理的气溶胶模型，通过卫星观测的表现反射率根据公式（2）计算可得出气溶胶光学厚度[3-4]。

通过暗像元反演流程图如图1所示。

图1 反演流程图Fig.1 Flow chart of the inversion experiment

2 软件设计

2.1 软件设计的问题

该气溶胶光学厚度反演软件是基于windows 7系统操作环境下，编程采用了C++语言，开发工具选用QT C++，调用GDAL开源库、HDF驱动进行读取、处理HDF格式文件。系统结构功能图如图2所示。

图2 系统结构功能图Fig.2 Structure diagram of the system’function

对MODIS遥感图像数据需要解决如下问题：

1）MODIS遥感影像数据采用HDF文件格式储存。利用开源GDAL库的API函数提取HDF文件的栅格数据及属性实现对MODIS数据的读取和写入。

2）对MODIS数据进行几何纠正，其中包括一般的遥感图像校正和MODIS数据特有的“Bowtie”（蝴蝶结效应纠正）[5]。同时能够实现对遥感图像的一般操作，如图像的放大，缩小以及漫游等功能。

3）软件系统集成嵌入6S模块，实现利用6S大气辐射传输模式，建立起表观反射率—地表反射率—气溶胶光学厚度的查找表。

2.2 软件设计的方法

GDAL（Geospatial Data Abstraction Library）是一个在 X/MIT许可协议下的开源栅格空间数据转换库。它主要可以对栅格数据格式进行相关数据处理，它使用抽象数据模型来分析数据格式，抽象的数据模型包括元数据、栅格波段、子数据集域图、图像结构域、仿射地理坐标转换等。MODIS遥感数据采用HDF格式进行存储，这种文件格式可以存储不同类型的图像和数码数据。HDF包含有6种主要数据类型：栅格图象，调色板，科学数据库，注释，Vdata和Vgroup。在GDAL中每个文件有一个数据集，HDF文件中包含多个数据集，可以将这些数据集提取到抽象模型中的子数据域，将属性数据提取到子数据域里的元数据中[6]。利用GDAL对HDF的数据处理可以概括为：首先打开文件，获得子数据名列表，根据列表打开指定的数据集。然后获取子数据集的属性信息。打开需要的波段，提取波段数据。将数据信息转换后关闭数据集并显示。

1）运行环境配置

MODIS卫星数据一个文件包含多个SUBDATASETSGDAL，每个SUBDATASETSGDAL又包含很多波段的数据。GDAL库的默认支持文件不包含HDF格式，因此需要向GDAL源文件中添加HDF的驱动重新编译。在GDAL文件中的nmake.opt中修改GDAL_HOME和#Uncomment the following and update to enable NCSA HDF Release 4 support，添加HDF的函数库路径。然后用命令工具进行编译，将编译好的GDAL中的bin文件夹添加到系统环境变量中。

2）图像读写和处理模块

在 GDAL读取栅数据格式文件时首先用GDALAllRegister（）注册HDF文件驱动。打开指定的HDF文件的子数据集并获取该子数据集里的波段数据指针。GDALRasterBand::RasterIO（）函数可以用来提取多波段信息，可以根据情况自动进行数据转换，以及对数据窗口放大和缩小等。需要指出的是，波段数据是以行为单位进行读取的。通过以上步骤可以对HDF文件指定的波段进行处理，完成处理工作后，可以通过GDALClose（）函数关闭打开的子数据集和HDF文件。GDAL读取HDF关键代码实现如下：

MODIS数据虽然已经过辐射纠正，但是MODIS L1B数据仍存在自身缺陷，即具有蝴蝶结效应。所以在使用MODIS数据之前，必需要先进行几何校正。在几何纠正功能设计上采用投影—插值的方法进行相关编程，这种方法实际上是在投影变换的同时进行了去除蝴蝶结效应。软件的读写模块主要实现打开、保存和退出，相应的功能按钮为“打开影像”，“保存”和“退出”。图像处理模块主要实现对图像进行几何纠正（包括蝴蝶结效应纠正和一般的几何纠正），相应的功能按钮为“几何纠正”。

3）图像操作模块

操作模块包括在对当前地图窗口的图形进行缩放，将图像置为全图模式，漫游实现手动调整窗口位置。前后视图实现打开上下视图。在鹰眼图中能够实现对图像的选择和操作。 相应的功能按钮为“放大”，“缩小”，“漫游”，“全图”，“前视图”和“后视图”。

4）6S 模块

6S大气辐射模式能准确模拟太阳到目标物，目标物到传感器路径上的大气影响，是发展比较成熟的大气修正模式。由于6S模式本身相对精确和复杂，并且由FORTRAN语言编写。所以在添加6S模式时并没有进行C++语言的重新编写，而是将FORTRAN语言的源程序编译运行后生成的debug文件夹下的main.exe用ShellExecute函数调用，将输入界面的值保存到输入文件inputfile.txt中。在运行6S模型后，调用debug文件夹下的main.exe后，生成输出文件output.txt。6S功能模块相应的功能按钮为“6S”。

3 功能实现

3.1 读取和处理遥感图像

打开影像文件，利用软件对数据中1、3、7波段进行几何校正，目的是去除上文所述的蝴蝶结效应和进行一般的几何纠正。第7波段几何校正前后图如图3，4所示。

图3 第7波段几何校正前Fig.3 Band7 before the geometric correction

图4 第7波段几何校正后Fig.4 7 band after the geometric correction

经过几何校正后，可在第7波段图中寻找可以用于暗像元法反演的点，该点要满足表观反射率小于0.05，并且该点附近的暗像元点的边关反射率变化不大。本论文所采用暗像元点为北纬19.056 461，东经126.009 181的点。在该点处获得相关参数。

3.2 调用6S大气辐射传输模型

在6S选择项中跟着提示操作在几何参数中依次输入太阳天顶角、卫星天顶角、太阳方位角、卫星方位角4个变量。大气模式包括7种，依据暗像元所在地区选择中纬度夏季大气模式。气溶胶模式包括8种，将气溶胶类型设置为海洋性气溶胶。目标物和传感器高度选择100 m。地表状况处理设为地物面为均匀朗伯表面，无方向效应[7]。

3.3 建立查算表

在6S模型利用上述的几何参数，结合了不同大气模式、气溶胶模式，根据550 nm气溶胶光学厚度值和暗像元1、3波段的地表反射率反演出该暗像元点的地表反射率。1、3波段查算表建立完成后，可根据该暗像元点的实际地表反射率、表观反射率分别对照1、3波段查算表找出最匹配的AOD值，两波段对应AOD的算数平均值可近似等于该暗像元点550 nm气溶胶光学厚度。第1、3波段表观反射率查算表的部分如表格1、2所示。

表1 第1波段部分查算表Tab.1 Part of the band1’s lookup-table

4 结 论

文中设计和开发了一种基于暗像元法的大气气溶胶光学厚度反演的一体化软件，软件设计主要采用GDAL库等对MODIS遥感数据进行提取和相关处理。通过实际应用表明软件运行稳定可靠，GDAL和HDF能够实现完美匹配和兼容，最终能够实现气溶胶光学厚度的反演目的。这种软件的开发能够有效的简化了使用者的操作过程，降低了对使用者的要求，提高了效率。

[1]毛节泰，李成才，张军华，等.MODIS卫星遥感北京地区气溶胶光学厚度及与地面光度计遥感的对比 [J].应用气象学报，2002，13（特刊）:127-135.MAO Jie-tai，LICheng-cai，ZHANG Jun-hua，etal.The comparison of remote sensing aerosol optical depth from MODIS data and ground sun-photometer observations[J].Journal of Applied Meterological Science，2002，13 （specialissue）:127-135.

表2 第3波段部分查算表Tab.2 Part of the band3’s lookup-table

[2]张岩岫，王志清，刘立武，等.大气散射对激光制导武器对抗态势构建影响研究[J].现代电子技术，2012（21）：38-40，44.ZHANG Yan-xiu，WANG Zhi-qing，LIU Li-wu，et al.Influence of atmospheric scattering on confrontation building in laser guided weapon[J].Modern Electronics Technique，2012（21）：38-40，44.

[3]李晓静，刘玉洁，邱红，等.利用MODIS资料反演北京及其周边地区气溶胶光学厚度方法研究[J].气象学报，2003，61（5）:580-591.LI Xiao-jing，LIU Yu-jie，QIU Hong，et al.Retrieval method for optical thickness of aerosds over Beijing and it’s vicinity by using the MODIS data[J].Acta Meteorological Sinica，2003，61（5）:580-591.

[4]Kaufman Y J，Holben B N.Will aerosol measurements from Terra and Aqua polar orbiting satellites represent the daily aerosol abundance and properties[J].Geophysics Research Letters，2000，27（23）:3861-3864.

[5]郭广猛.关于MODIS卫星数据的几何纠正方法 [J].遥感信息，2002（3）:26-28.GUO Guang-meng.Geometric calibration of MODIS data[J].Remote Sensing Information，2002（3）:26-28.

[6]张莉，曾致远.基于HDF4文件格式的MODIS 1B影像数据提取的研究与实现[J].国土资源感，2004（4）:27-32.ZHANG Li，ZENG Zhi-yuan.The study and implementation of extraction MODIS level 1B image data based on a HDF4 file[J].Remote Sensing for Land&Resources，2004（4）:27-32.

[7]Verraote E F，Tanre D，Denze J L，et al.Second simulation of the satellite signal in the solar spectrum,6S:anoverview[J].IEEE Trans Geosci Remote Sens，1997（35）:675-686.