Print

出版日期: 2016-07-25
点击次数：
下载次数：
DOI: 10.11834/jrs.20164287
2016 | Volumn20 | Number 4





                              上一篇|





下一篇
大气遥感专栏
结构函数法反演气溶胶光学厚度中像元的间隔设置
expand article info 朱琳1,2 , 孙林3 , 杨磊4 , 徐菲菲5,6 , 徐青山1
1. 中国科学院安徽光学精密机械研究所 中国科学院大气成分与光学重点实验室, 安徽 合肥 230031
2. 中国科学技术大学, 安徽 合肥 230027
3. 山东科技大学测绘科学与工程学院, 山东 青岛 266590
4. 国家海洋局第一海洋研究所, 山东 青岛 266061
5. 中国科学院上海技术物理研究所, 上海 200083
6. 中国科学院大学, 北京 100049
收稿日期: 2015-04-13; 修改日期: 2015-12-21; 优先数字出版日期: 2016-07-25
基金项目: 国家自然科学基金（编号：41171270）；中央级公益性科研院所基本科研业务费专项资金资助项目（编号：2014623）
第一作者简介: 朱琳(1989-),女,博士研究生,主要研究方向大气探测与遥感、海洋及陆地气溶胶卫星反演。E-mail:zhulin0706@hotmail.com
中图分类号: TP79
文献标识码: A
文章编号: 1007-4619(2016)04-0528-12

摘要

结构函数法气溶胶光学厚度反演精度受像元间隔设置影响很大，且并非所有像元都能获得较好的反演结果，因而研究像元间隔的设置能够提高反演精度，研究反演误差小的像元能够提高算法效率。为了获得最佳的像元间隔设置，本文以胶州湾地区为例，利用250 m和500 m两种分辨率数据计算了不同像元间隔时的结构函数值，分别利用单一像元间隔法、均值法、坡高法以及线性区域均值法获得待反演像元最终的结构函数值反演550 nm处的气溶胶光学厚度，并依据CE318观测数据进行精度验证，通过分析点对点反演结果和光学厚度的空间分布，确定反演误差小、受分辨率影响小的像元间隔设置。实验发现线性区域均值法在一定程度上提高了反演精度和稳定性。此外，通过对反演结果可接受像元的地表反射率结构函数值的统计和分析，发现500 m分辨率时可接受像元比例优于250 m，当地表反射率结构函数值大于0.02时反演结果较好，而这些像元往往分布在山麓、山涧、海岸线、河流、城乡结合部等地理要素的突然改变的地区。

关键词

结构函数法 , 地表反射率 , 表观反射率 , 像元间隔 , 气溶胶光学厚度

Pixel distance settings in aerosol optical depth retrieval through the structure function method
expand article info ZHU Lin1,2 , SUN Lin3 , YANG Lei4 , XU Feifei5,6 , XU Qingshan1
1.Key Laboratory of Atmospheric Composition and Optical Radiation, Anhui Institute of Optics and Fine Mechanics, Chinese Academy of Sciences, Hefei 230031, China
2.University of Science and Technology of China, Hefei 230027, China
3.Shandong University of Science and Technology, Qingdao 266590, China
4.The First Institute of Oceanography, SOA, Qingdao 266061, China
5.Shanghai Institute of Technical Physics of the Chinese Academy of Science, Shanghai 200083, China
6.University of Chinese Academy of Science, Beijing 100049, China

Abstract

Retrieval of Aerosol Optical Depth(AOD)over land through satellite remote sensing is difficult. Traditional methods to retrieve aerosol over land involve the use of the Dense Dark Vegetation(DDV) algorithm. However, this algorithm requires the support of an infrared band and dark pixels. Compared with DDV, the structure function method can be utilized to retrieve AOD over bright areas with high surface reflectance. However, the precision of retrieving AOD through the structure function method can be significantly influenced by pixel distance. Not all pixels can be retrieved with high precision; this condition results in a large amount of unnecessary computations. As a result, determining the pixel distance settings and the pixels with a small retrieval error is important in improving the retrieval precision and efficiency of the structure function method. To acquire the best pixel distance settings, we calculated the structure function values with different pixel distances and obtained the final structure function value of the target pixel through single-pixel distance, average, slopeheight, and average linear area methods to retrieve AOD at 550 nm over Kiaochow Bay. The first band of nine-scene MODIS L1B data with resolutions of 250 m and 500 m was used. Surface reflectance was the minimum value of the first band of MOD09Q1 in 2012. By validating the AOD data measured by CE318 and by analyzing the retrieval result and spatial distribution of AOD, we obtained the pixel distance settings with high precision and minimal effects of resolution. The experiment indicated that the average linear area method has a small retrieval error whether the resolution is 250 m or 500 m. According to statistics and by analyzing the structure function value of the surface reflectance of pixels with an acceptable error (the absolute value of the absolute error is 0.1 for measured AOD less than 0.6 and 0.2 for measured AOD greater than 0.6), we found that the ratio of pixels with are solution of 500 m is obviously larger than that with 250 m resolution. Pixels with a surface reflectance structure function value larger than 0.02 have a good result and areaways distributed over areas with a sudden change in geographical elements, such as piedmonts, mountain streams, coasts, rivers, and rural-urban continua. In this study, we utilized different methods and historical surface reflectance data with two resolutions to calculate the final structure function value. The retrieval precision and distribution of AOD were analyzed to acquire the best pixel distance settings, and its features were identified through the statistics of the structure function values with high retrieval precision. The average linear area method is stable, and its result has an insignificant relationship with data resolution. The retrieved AOD with high precision is concentrated on pixels with a surface reflectance structure function value larger than 0.02. The efficiency of this method increases when the measured AOD is greater than 0.5.

Key words

structure function method , surface reflectance , apparent reflectance , pixel distance , AOD(Aerosol Optical Depth)

1 引 言

结构函数法利用同一地区的多幅图像反演气溶胶光学厚度(AOD)，从中选取地表分布最清晰的一幅作为参考，假设这些图像对应的地表反射率空间分布并没有因为时间、人为等因素发生变化，那么从这一地区其他图像的表观反射率空间分布中去除地表反射率的贡献，就可以确定大气气溶胶的贡献，从而得到气溶胶光学厚度(Tanré等，1988)。

结构函数法通过不同图像的结构函数值与大气透过率之间的关系确定气溶胶含量，能够获取荒漠、城市等亮地表上空的气溶胶光学厚度，反演得到的光学厚度不受气溶胶散射相函数的影响(Tanré等，1988)。地表反射率和表观反射率的结构函数值的比值可以理解为大气透过率，结构函数值的准确度直接影响到大气透过率的取值，更关乎查找表插值得到的气溶胶光学厚度。然而，结构函数模型的实质是计算一定窗口内一定像元间隔下反射率差异值的均方根，因此，影响结构函数值准确度的主要因素就是像元间隔的设置。Tanré等人(1988)在提出结构函数法时，指出规则的结构函数值是先随像元间隔d的增加而变大，而后保持平稳，反演利用了结构函数值从小到大的坡高确定待反演像元的结构函数值。Liu等人(2002a)；Holben等人(1992)在利用SPOT HRV数据反演台湾地区的气溶胶光学厚度用于监测空气质量时，发现受地形、观测几何等的影响，结构函数值随d的变化并不是规则的，因此提出了3个方向对比的结构函数模型降低d的影响，并利用像元间隔从1到10的结构函数值的均值法反演气溶胶光学厚度。Liu等人(2002a)通过像元间隔从1到10的频数统计，发现像元间隔超过5时，结构函数值随像元间隔的变化存在负相关的现象，剔除了像元间隔超过5的结构函数值，然后再利用均值法反演气溶胶光学厚度，反演结果与不剔除相比，相对误差最大可以减小45%，表明像元间隔对反演结果影响很大。孙林等人(2006)针对城市地区复杂多变的地表改进了结构函数模型，提出了点与区域对比的模型，降低了像元间隔对结构函数值的影响，在使用MODIS数据反演北京城区的光学厚度时，考虑了城市BRDF特性，取得了较高的气溶胶反演精度。已有研究表明，结构函数法需要从一系列图像中选取最清晰的一幅进行大气校正获取地表反射率信息，清晰图像的选取直接影响了结构函数值的计算结果(Holben等，1992；Liu等，2002a；孙林，2006)。由于云、雾、雨等天气的影响，清晰图像在短时间内是无法获取的，然而，对历史地表反射率产品进行最小值合成可以获得受天气等因素影响小的最接近实测值的地表反射率数据。据此，徐菲菲等提出使用历史地表反射率数据代替参考图像支持结构函数法，采用点与区域对比的结构函数模型反演北京地区的气溶胶光学厚度，并与利用实时地表反射率数据反演的气溶胶光学厚度比较，取得了较好的反演结果(徐菲菲，2012)。尽管有相当多的研究改进了结构函数法，但算法在实际反演中仍然不稳定，只有一部分像元的反演结果可以接受，分析和研究这些像元结构函数值的特点可以提高算法反演效率。然而，结构函数值的计算关键在于像元间隔的设置。因此，获取最佳的像元间隔设置对结构函数法的实际应用意义重大。

为了获取最佳的像元间隔设置，得出反演结果可接受像元的地表反射率结构函数值特点，本文以胶州湾地区为例，利用250 m和500 m两种分辨率数据计算了不同像元间隔d从1到10时的结构函数值，根据结构函数值随d的变化关系，利用坡高法、均值法、单一像元间隔法和线性区域均值法得到最终的结构函数值反演气溶胶光学厚度。通过分析不同方法得到的气溶胶光学厚度的分布直方图、点对点反演误差和空间分布特征，确定反演精度较高、受分辨率影响小的像元间隔设置，统计反演结果可接受和不可接受像元的地表反射率结构函数值，根据结果可接受与不可接受像元的数量差异确定最佳的数据分辨率和相应的地表反射率结构函数值特点。本文对三方向对比的结构函数模型像元间隔的研究，得出了适合胶州湾地区的像元间隔设置和适于反演的地表反射率结构函数值特点，提高了结构函数法的反演精度与稳定性。

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

假设地表是朗伯体、大气水平分布均匀，传感器接收到的表观反射率${\rho _{_{_{{\rm{TOA}}}}}}({\mu _{\rm{s}}}，{\mu _{\rm{v}}}，\phi)$可表示为(Tanré等，1988；Liou，2002a)：

$\begin{align} & {{\rho }_{\text{TOA}}}({{\mu }_{\text{s}}}\text{, }{{\mu }_{\text{v}}}\text{,}\phi )={{\rho }_{\text{0}}}({{\mu }_{\text{s}}}\text{,}{{\mu }_{\text{v}}}\text{,}\phi )+ \\ & \frac{T({{\mu }_{\text{s}}})\exp (-\tau /{{\mu }_{\text{v}}}){{\rho }_{\text{s}}}({{\mu }_{\text{s}}}\text{,}{{\mu }_{\text{v}}}\text{,}\phi )}{[1-{{\rho }_{\text{s}}}({{\mu }_{\text{s}}}\text{,}{{\mu }_{\text{v}}}\text{,}\phi )S]} \\ \end{align}$ (1)

式中，${\rho _{_{_{\rm{0}}}}}({\mu _{\rm{s}}}，{\mu _{\rm{v}}}，\phi)$为大气程辐射等效反射率；${\rho _{\rm{s}}}({\mu _{\rm{s}}}，{\mu _{\rm{v}}}，\phi)$为地表二向反射率，当地表为朗伯体时等于地表反射率；${\mu _{\rm{s}}} = \cos {\theta _{\rm{s}}}$，${\mu _{\rm{v}}} = \cos {\theta _{\rm{v}}}$，${\theta _{\rm{s}}}$与${\theta _{\rm{v}}}$分别为太阳天顶角和观测天顶角；${\phi}$为相对方位角；S为大气半球反射率；T(μs)为大气下行透过率；exp(-τ/μv)为上行透过率。

由式(1)可知，要从${\rho _0}({\mu _{\rm{s}}}，{\mu _{\rm{v}}}，\phi)$中提取气溶胶信息，需要确定气溶胶模式和地表反射噪声。当地表反射率较低时，传感器接收到的辐射值主要受大气及气溶胶影响，由于水的反射率很小，地表反射噪声易去除，因此早期气溶胶反演主要集中在海洋等水体上空(Tanré等，1997)。对于水体分布有限的大陆地区，气溶胶遥感主要采用浓密植被法，根据浓密植被等暗像元在近红外波段处的表观反射率确定红蓝波段的地表反射率，利用红蓝波段的反射率与气溶胶光学厚度的对应关系确定光学厚度(Kaufman等，1997)。随着地表反射率的增加，气溶胶对传感器接收到的辐射值的影响逐渐变小，当地表反射率增加到一定程度，辐射值将不随气溶胶光学厚度变化。研究表明，在反射率较高的地区，地表反射贡献难去除，而浓密植被法由于缺少暗像元不能实现高反射率地区的气溶胶反演(Durkee等，1986；Fraser和Kaufman，1985)。为实现高反射率地区的气溶胶反演，Tanré等人(1988)从存在“模糊效应”的一组图像中反演了干旱地区的气溶胶光学厚度，该方法假设地表反射率空间分布相对不变，那么图像的“模糊效应”主要是由气溶胶的变化引起的，从而确定气溶胶光学厚度，由于反射率的空间分布是用结构函数值描述的，所以称这种基于“模糊效应”的反演方法为结构函数法。

假设大气条件不变，根据式(1)可知图像上相邻两像元(x，y)和(x，y+1)的表观反射率差异$\Delta \rho _{x，y}^*({\tau _{\rm{a}}}，{\mu _{\rm{s}}}，{\mu _{\rm{v}}})$与相应的地表反射率差异△ρx，y满足：

$\Delta \rho _{x\text{,}y}^{*}({{\tau }_{\text{a}}}\text{,}{{\mu }_{\text{s}}}\text{,}{{\mu }_{\text{v}}})\approx \Delta {{\rho }_{x\text{,}y}}\frac{T({{\mu }_{\text{s}}})\exp (-\tau \text{ }/{{\mu }_{\text{v}}})}{1-\bar{\rho }S}$ (2)

式中，$\bar \rho$表示相邻两像元的地表反射率均值。

从式(2)可知，当地表反射率不变时，可通过表观反射率的变化确定光学厚度。由于结构函数法需要清晰图像地表反射率的空间分布作为参考，因此对参考图像计算经大气校正所得地表反射率的结构函数值${M_{{\rho}}}(d，{t_1})$，并将待反演图像表观反射率的结构函数值记为${M_{{\rho }}}^*(d，{t_2})$，忽略地表与大气之间的多次散射(Tanré等，1981)则有：

${{M}_{\rho }}^{*}(d\text{,}{{t}_{2}})={{M}_{\rho }}(d\text{,}{{t}_{2}})T({{\mu }_{\text{s}}})\exp (-\tau /{{\mu }_{\text{v}}})$ (3)

而结构函数法假设地表反射率保持不变，所以${M_{{\rho }}}(d，{t_2})= {M_{{\rho }}}(d，{t_1})$，式(3)可转化为

${{M}_{\rho }}^{*}(d\text{,}{{t}_{2}})={{M}_{\rho }}(d\text{,}{{t}_{1}})T({{\mu }_{\text{s}}})\exp (-\tau /{{\mu }_{\text{v}}})$ (4)

由式(4)可知，结构函数法获取气溶胶光学厚度关键在于地表反射率和表观反射率结构函数值的计算。目前，用于计算结构函数的模型主要有单方向对比、3个方向对比、点与区域对比3种。单方向对比的结构函数模型是与同一行间隔一定的像元对比(Holben等，1992)，3个方向对比的结构函数模型将像元与其同一行、同一列及对角方向一定间隔的像元对比(Liu等，2002a)，孙林等针对城市地区复杂的地表分布定义了点与间隔一定的像元集合对比的结构函数模型(孙林，2006)。

$\ln [{{M}_{\rho }}^{*}(d\text{,}{{t}_{2}})]\cong \ln [{{M}_{\rho }}(d\text{,}{{t}_{1}})]+\ln [T({{\mu }_{\text{s}}})]-\tau /{{\mu }_{\text{v}}}$ (5)

气溶胶反演是根据式(5)实现的，由于结构函数值受模型参数的影响很大(孙林等，2006)，已有的结构函数法气溶胶反演建立在大量的结构函数值的计算上，过程相对繁琐，而单一的参数设置带来的反演误差很大。为简化结构函数值的反演，降低参数设置对反演结果的影响，确定最佳的参数组合，本文以胶州湾地区为例，使用了多种方法获取的结构函数值进行结构函数法气溶胶光学厚度反演。

3 实验数据及流程

胶州湾地区地处崂山山脉以西，大珠山小珠山以东，包含青岛市主城四区、青岛经济技术开发区、城阳、胶州南部、胶南东部等人居区域，围绕胶州湾从沿岸向内陆辐射的团带状经济地带。但从地形地貌上讲，胶州湾是海滨丘陵地区，地势东高西低，南北两侧隆起。虽然胶州湾地区植被覆盖度较高，也有不少水体，十分适宜用暗目标法、浓密植被法反演气溶胶光学厚度，但是这两种算法至少需要两个波段的支持，然而结构函数法只需要使用一个波段。此外，结构函数法是在对比算法的基础上发展的，像元间反射率差异越大，反演效果越好，而植被和水体的存在会提高反射率间的差异，因此，本文以胶州湾地区为研究对象，研究适用于胶州湾地区的最佳结构函数模型参数组合。实验采用历史地表反射率产品进行结构函数法反演，利用2012年全年的地表反射率产品(MOD09Q1)的红光波段使用最小值合成得到的如图 1所示的地表反射率数据用于结构函数法气溶胶光学厚度反演，待反演图像是MODIS L1B数据(MOD02)经过重投影、辐射定标后得到的(Remer等，2002；Chu等，2002；Levy等，2004)，同时为了研究分辨率对结构函数值的影响，本文采用250 m、500 m两种分辨率的数据。


图 1 地表反射率图像及CE318观测位置
Fig. 1 The image of surface reflectance and the position of CE318

综合考虑实验区域复杂的地形和3种结构函数计算模型的耗时问题，采用式(6)所示的3个方向对比的模型用于结构函数值的计算(Liu等，2002a，2002b)：

$ \begin{split} {M^2}(d)=\frac{1}{{3(n - d)(m - d)}}\sum\limits_{i = 1}^{n - d} {\sum\limits_{j = 1}^{m - d} {[{{({\rho _{i，j}} - {\rho _{i，j + d}})}^2} + } } \\ {({\rho _{i，j}} - {\rho _{i + d，j}})^2} + {({\rho _{i，j}} - {\rho _{i + d，j + d}})^2}] \end{split} $ (6)

式中，模型参数m=n=15，d=1，2，3，4，5，6，7，8，9，10。

由式(5)可知，得到地表结构函数值和表观函数值后，就可以根据大气透过率与光学厚度之间的关系得出光学厚度值。为得到大气透过率与气溶胶光学厚度的对应关系，本文采用6S模型构建不同太阳天顶角、观测天顶角、气溶胶光学厚度和$\ln [T({{\mu }_{\text{s}}})]-\tau /{{\mu }_{\text{v}}}$的查找表(Vermote等，1997；Li等，2014)。利用表观反射率与地表反射率的结构函数比值与$\ln [T({{\mu }_{\text{s}}})]-\tau \text{ }/{{\mu }_{\text{v}}}$进行线性插值得到550 nm处的气溶胶光学厚度，从而实现结构函数法气溶胶光学厚度的反演。其中，建立查找表时气溶胶光学厚度设置如下：0.00001、0.05、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9、1.0、1.2、1.4、1.6、1.8、2.0、2.5。

为验证结构函数法的反演结果，本文利用安置在山东科技大学测绘科学与工程学院(120°12′E，36°0′N)的CE318太阳光度计的同步观测数据进行误差分析。CE318是由CIMEL公司生产的可以自动追踪太阳测量可见光和近红外的不同波段、不同方向、不同时间的辐射量度，进而推算大气气溶胶光学厚度(Holben等，1998；Li等，2008；Huang等，2008)，实验所用验证数据如表 1所示。需要指出的是，实验依据6S查找表插值获得的气溶胶光学厚度是550 nm处的，而CE318的观测波段并不包含550 nm，因此要利用式(7)所示关系式计算出550 nm处的观测值(Khoshsima等，2014；Simha等，2013；Beegum等，2009；Balakrishnaiah等，2011)。

$ {\tau _a}(\lambda)= \beta {\lambda ^{ - \alpha }} $ (7)

式中，${\tau _a}(\lambda)$是波长$\lambda $处的气溶胶光学厚度，α为Angstrom指数，描述气溶胶光学厚度与波长之间的关系，β为大气浑浊度系数。

表 1 与MOD02成像时间临近的CE318观测值
Table 1 The AOD values at the imaging time of MOD02 measured by CE318 sunphotometer

下载CSV 
日期 440 nm 500 nm 670 nm 870 nm 550 nm
2015-01-01 0.183 0.154 0.102 0.079 0.139
2015-01-02 0.526 0.443 0.296 0.215 0.392
2015-01-17 0.416 0.369 0.262 0.187 0.320
2015-02-14 0.39 0.348 0.242 0.174 0.299
2015-02-17 0.994 0.933 0.731 0.604 0.844
2015-02-23 0.742 0.717 0.657 0.633 0.704
2015-03-01 0.705 0.66 0.55 0.488 0.625
2015-03-06 0.863 0.716 0.46 0.327 0.628
2015-03-11 0.659 0.556 0.367 0.261 0.487

为了合理选取验证数据，实验对2015年1月1日至2015年3月11日CE318的所有观测数据进行计算和统计，得到了如图 2所示的550 nm气溶胶光学厚度统计图，并在图 2中标注了0.2、0.4、1.0 3条分界线。


图 2 550 nm处的气溶胶光学厚度统计图
Fig. 2 Statistical graph of AOD at 550 nm

由图 2可知，气溶胶光学厚度小于0.2出现的次数很少，0.2到0.4之间的相对较多，而且光学厚度越大出现的次数越多，大多数的取值都在0.4到1.0之间，也有部分观测结果大于1，因此，验证数据的选取也遵循这个规律，小于0.2的有1天，介于0.2到0.4的有3天，0.4到1之间有5天。

4 实验结果与分析

为了得到适宜胶州湾地区的像元间隔设置，本文以地面CE318观测位置处的像元为例，探究胶州湾地区结构函数值随像元间隔d的变化，计算了d从1到10时地表反射率与待反演图像表观反射率的结构函数值，结果如图 3所示。


图 3 CE318所处像元的结构函数值随d值的变化
Fig. 3 The variation of structure function values due to the distance d at the pixel of CE318

由图 3可知，250 m分辨率数据得到的结构函数值随像元间隔的变化显著大于500 m分辨率数据，这是由于250 m分辨率数据能够捕捉到较为细小的反射率变化。但是，两种分辨率数据对应的地表反射率和表观反射率的结构函数值随d的变化都不是规则的，因此，无法选取理想的斜坡高度进行反演，而结构函数值的变化使得d=1到d=10的均值失去了代表性。然而，当d=1到d=4时，大部分数据的结构函数值都在一条直线上。为了降低像元间隔对反演结果的影响，确定适合胶州湾地区的最佳像元间隔设置，本文分别使用(1) 单一像元间隔法：m=n=10，d=5的结构函数值；(2) 线性区域均值法：m=n=15，d=1到d=4的结构函数值的均值；(3) 均值法：m=n=15，d=1到d=10的结构函数值的均值；(4) 坡高法：m=n=15，d=4与d=1的结构函数值之差四种不同的方法进行结构函数法反演，4种方法用两种分辨率数据反演得到的550 nm气溶胶光学厚度如图 4、图 5所示。


图 4 500 m分辨率数据下4种方法所得550 nm气溶胶光学厚度的分布直方图
Fig. 4 Distribution histogram of AOD at 550 nm retrieved by the four methods with data resolution of 500 m

图 4和图 5中两条蓝线之间的区域是根据观测值确定的反演结果可接受范围，对于观测值小于0.6取±0.1的绝对误差，对于观测值大于0.6的取±20%的相对误差。从图 4可知，除法1外，其他3种方法的反演结果分布具有较高的一致性，且分布直方图都存在一个较明显的峰值。此外，除前四天的反演结果外，其余五天的反演结果中法2的可接受像元明显多于其他3种方法。结合图 4、图 5可知，250 m分辨率下4种方法反演结果分布的差异相对500 m要小，直方图上并不存在明显的峰值，法2的反演优势只在第五天才可以看出。为了深入分析4种方法的反演特点，实验对两种分辨率得到的结果进行了精度验证，由于CE318的观测局限在特定的经纬度，而遥感数据是面元，因此，实验利用CE318所在经纬度提取了最邻近的像元，反演结果如表 2所示。


图 5 250 m分辨率数据下4种方法所得气溶胶光学厚度的分布直方图
Fig. 5 Distribution histogram of AOD retrieved by the four methods with data resolution of 250 m

表 2 CE318所在像元的反演结果
Table 2 The retrieval results of different methods at the pixel of CE318

下载CSV 
日期 观测值 分辨率/m 方法1 方法2 方法3 方法4
2015-01-01 0.139 500 0.7598 0.5937 0.6091 0.6042
250 0.8223 0.5664 0.6388 0.6119
2015-01-02 0.392 500 0.6491 0.6380 0.6419 0.6435
250 0.5159 0.5553 0.5735 0.5836
2015-01-17 0.320 500 0.9299 0.6913 0.7268 0.6908
250 0.9821 0.7364 0.8131 0.8167
2015-02-14 0.299 500 0.8871 0.5796 0.6184 0.5801
250 1.0625 0.6916 0.7914 0.7684
2015-02-17 0.844 500 0.8245 0.7862 0.7395 0.8258
250 0.9762 0.8494 0.8864 0.9469
2015-02-23 0.704 500 1.1386 0.8293 0.8728 0.8346
250 0.9991 0.8639 0.9260 0.9339
2015-03-01 0.625 500 0.8262 0.7123 0.7117 0.7111
250 0.6553 0.8635 0.8515 0.8675
2015-03-06 0.628 500 0.7649 0.6026 0.6281 0.6134
250 0.6792 0.6432 0.6810 0.7180
2015-03-11 0.487 500 0.8530 0.5705 0.6231 0.5590
250 0.7583 0.6490 0.7129 0.6391

由表 2可知，4种方法在前4天的反演误差都很大，除了方法2外，其他4种方法的反演结果受分辨率的影响都较大。对于最后一天，250 m的反演误差明显很大，而500 m的反演误差相对较小，但方法1和方法3的反演误差明显较大。对于剩余4天的反演结果，除了方法1的反演误差时大时小外，其他3种方法的反演结果基本可以接受。综合来说，方法1的反演结果不稳定，而且受分辨率的影响很大，其他3种方法的反演结果各有千秋，但方法2的反演结果受分辨率的影响较方法3和方法4要小。尽管方法2较其他3种方法具有一定的稳定性，但利用两种分辨率数据得到的前四天的反演结果误差都很大。然而，图 2和表 2可知，方法2虽然不能反演小于0.2的气溶胶光学厚度，对0.2到0.5之间的气溶胶光学厚度反演结果也不太理想，但其反演结果比其他4种方法都要小，且反演结果都大于0.5，当光学厚度超过0.5时，反演结果在误差范围内，且250 m分辨率数据的反演结果相对500 m更接近观测值，这表明方法2的结构函数值确定方法可以用于反演胶州湾地区的气溶胶光学厚度。

为了进一步研究4种方法的特点，本文选取了2015年2月14日和2015年3月11日两种分辨率数据的反演结果，对比了4种方法所得气溶胶光学厚度的空间分布，如图 6所示，绿色区域则表示图 4、图 5中标注的反演结果可接受范围。

由图 6可知，在使用相同的色阶进行密度分割时，250 m分辨率的反演结果空间分布较500 m分辨率的离散。然而，无论哪一种分辨率，方法1与方法4的气溶胶分布较为离散，方法2与方法3则相对集中。将方法2与方法3反演结果的空间分布和表观反射率图像比较，可知方法2的分布较方法3更符合表观反射率的分布规律，而250 m分辨率数据更为细致。因此，综合考虑四种方法的时耗、点对点反演结果及反演结果的空间分布，认为在胶州湾地区采用方法2的像元间隔设置能够提高算法的反演精度与稳定性。图 6结合图 1所示的地表反射率图像可知，反演误差较小的像元集中在山麓、山涧、海岸线、河流、城乡结合部等明暗差异明显的地区。考虑待反演图像的表观反射率是不可控制的，无法通过表观反射率研究算法对结构函数值的要求，需要通过地表反射率确定适于反演的地表反射率结构函数值，图 7所示为两种分辨率下的地表反射率结构函数值分布图。


图 6 不同分辨率及不同方法得到的气溶胶光学厚度分布图
Fig. 6 The distribution of AOD get by different methods with different resolutions


图 7 不同分辨率下地表反射率结构函数值分布图
Fig. 7 The distribution of surface reflectance structure function values with different resolutions

由图 7可知，250 m分辨率对应的结构函数值总体比500 m的小，且分布相对离散，结合图 6可知，当地表反射率结构函数值介于0.02到0.03时，反演结果相对较好。由于图 6仅是两天的反演结果，为了得到更准确的结构函数值范围，根据线性区域均值法在两种分辨率下九个待反演图像的反演结果，实验分别统计了反演结果可接受和不可接受像元对应的地表反射率结构函数值，如图 8、图 9所示。


图 8 500 m分辨率数据所得地表反射率结构函数值统计结果
Fig. 8 Statistics of surface reflectance structure function values with 500 m resolution data

由图 8、图 9可知，除了前四天的反演误差较大，导致可接受的像元数量较少，结构函数值并不存在显著的特点，其他情况下，反演结果可以接受的像元地表反射率结构函数值大致成正态分布。其中，500 m分辨率时，地表反射率结构函数值取值在0.01到0.055之间，主要分布在0.02到0.04之间。250 m分辨率时，地表反射率结构函数值取值在0.007到0.05之间，且以0.01到0.03之间的像元数居多。此外，当观测值大于0.5时，250 m分辨率下反演结果不可接受的像元明显比500 m分辨率要多。因此，500 m分辨率数据更适于反演气溶胶光学厚度，而利用线性区域均值法反演气溶胶光学厚度对于地表反射率结构函数值大于0.02的像元具有较强的反演能力，而这些像元往往集中在地形、人为景观、环境等发生突然改变的地区，这些地区由于受地理要素等影响导致地表反射率变化显著，得到的结构函数值较大，因此反演得到的气溶胶光学厚度精度相对较高，更适于结构函数法的实际应用，并且结构函数值越大反演结果越好。然而，当地表反射率结构函数值小于0.01时，也有部分像元可以反演，但是这部分像元在图像中的比例很小。


图 9 250 m分辨率数据所得地表反射率结构函数值统计结果
Fig. 9 Statistics of surface reflectance structure function values with 250 m resolution data

不难发现，500 m分辨率下反演误差可接受的像元较误差过大的小，导致这个结果的原因可能是没有考虑像元实际的太阳天顶角和观测天顶角等角度信息，而是对整个图像设定了相同值。此外，气溶胶光学厚度的空间分布并不均匀，以观测点处的观测值衡量其他像元的反演结果，本身就会带来误差，在以后的工作中将重点分析研究两种误差下像元比重差异显著的原因。

5 结 论

本文以胶州湾地区为例，基于两种分辨率的历史地表反射率产品使用多种方法确定待反演像元最终的结构函数值，分析了不同方法的反演精度及光学厚度空间分布，统计了高精度像元的地表反射率结构函数值。实验发现不同的计算方法得到的气溶胶光学厚度均存在一定的误差，单一像元间隔设置的误差显著且不稳定，坡高法和均值法的反演结果受分辨率的影响大，而线性区域均值法的反演结果稳定且受分辨率影响较小，得到的气溶胶光学厚度空间分布相对集中，简化了算法流程，在一定程度上提高了算法的稳定性与精度。

通过对线性区域均值法在两种分辨率下所得反演结果可接受与不可接受像元的地表反射率结构函数值的统计与分析，发现反演精度与地表反射率结构函数值的空间分布密切相关，在山区、丘陵等地形变化明显的地区，湖泊、河流等反射率较低的水体附近以及人为景观变化明显的城市边缘等受地理要素、人为等因素影响导致地表反射率变化明显的地区，结构函数值相对较大，反演结果较为稳定，可以预先设置地表反射率结构函数值的阈值作为能否进行结构函数法反演的判断标准，超过阈值则进行反演。其中，500 m分辨率数据下可接受与不可接受像元在数量上的差异较250小许多，在利用阈值反演时能够提高效率。对于胶州湾地区而言，这个阈值可以设为0.02，大于0.02时就可以进行反演，一定程度上简化了反演流程，提高了结构函数法气溶胶光学厚度反演的效率，推动了算法的实际应用。

文中是针对胶州湾地区提出的参数组合与地表反射率结构函数值阈值，对于其他沿海及丘陵地区是否适用并未验证，在以后的工作中，将重点分析误差可接受与误差过大的像元地表反射率结构函数值的差异，并将结构函数法用于大区域的气溶胶反演。

参考文献(References)

  • Balakrishnaiah G, Raghavendra KumarK, Suresh Kumar Reddy B, Rama GopalK, Reddy R R, Reddy L S S, NazeerAhammed Y, Narasimhulu K, Krishna Moorthy K, Suresh Babu S.2011.Analysis of optical properties of atmospheric aerosols inferred from spectral AODs and Ångström wavelength exponent. Atmospheric Environment, 45 (6) : 1275–1285 . [DOI:10.1016/j.atmosenv.2010.12.002]
  • Beegum S N, Moorthy K K, Babu S S.2009.Aerosol microphysics over a tropical coastal station inferred from the spectral dependence of Angstrom wavelength exponent and inversion of spectral aerosol optical depths. Journal of Atmospheric and Solar-Terrestrial Physics, 71 (17/18) : 1846–1857 . [DOI:10.1016/j.jastp.2009.07.004]
  • Chu D A, Kaufman Y J, Ichoku C, Remer L A, Tanré D, Holben B N.2002.Validation of MODIS aerosol optical depth retrieval over land. Geophysical Research Letters, 29 (12) : MOD2-1–MOD2-4 . [DOI:10.1029/2001GL013205]
  • Durkee P A, Jensen D R, Hindman E E, VonderH T H.1986.The relationship between marine aerosol particles and satellite-detected radiance. Journal of Geophysical Research:Atmospheres, 91 (D3) : 4063–4072 . [DOI:10.1029/JD091iD03p04063]
  • Fraser R S, Kaufman Y J.1985.The relative importance of aerosol scattering and absorption in remote sensing. IEEE Transactions onGeoscience and Remote Sensing, GE-23 (5) : 625–633 . [DOI:10.1109/TGRS.1985.289380]
  • Holben B, Vermote E, Kaufman Y J, Tanré D, Kalb V.1992.Aerosol retrieval over land from AVHRR data-application for atmospheric correction. IEEE Transactions onGeoscience and Remote Sensing, 30 (2) : 212–222 . [DOI:10.1109/36.134072]
  • Holben B N, Eck T F, Slutsker I, Tanré D, Buis J P, Setzer A, Vermote E, Reagan J A, Kaufman Y J, Nakajima T, Lavenu F, Jankowiak I, Smirnov A.1998.AERONET-A federated instrument network and data archive for aerosol characterization. Remote Sensing of Environment, 66 (1) : 1–16 . [DOI:10.1016/S0034-4257(98)00031-5]
  • Huang J P, Zhang W, Zuo J Q, Bi J R, Shi J S, Wang X, Chang Z L, Huang Z W, Yang S, Zhang B D, Wang G Y, Feng G H, Yuan J Y, Zhang L, Zuo H C, Wang S G, Fu C B, Chou J F.2008.An overview of the semi-arid climate and environment research observatory over the Loess Plateau. Advances in Atmospheric Sciences, 25 (6) : 906–921 . [DOI:10.1007/s00376-008-0906-7]
  • Kaufman Y J, Tanré D, Remer L A, Vermote E F, Chu A, Holben B N.1997.Operational remote sensing of tropospheric aerosol over land from EOS moderate resolution imaging spectroradiometer. Journal of Geophysical Research:Atmospheres, 102 (D14) : 17051–17067 . [DOI:10.1029/96JD03988]
  • Khoshsima M, Ahmadi-Givi F, Bidokhti A A, Sabetghadam S.2014.Impact of meteorological parameters on relation between aerosol optical indices and air pollution in a sub-urban area. Journal of Aerosol Science, 68 : 46–57 . [DOI:10.1016/j.jaerosci.2013.10.008]
  • Levy R C, Remer L A, Kaufman Y J.2004.Effects of neglecting polarization on the MODIS aerosol retrieval over land. IEEE Transactions onGeoscience and Remote Sensing, 42 (11) : 2576–2583 . [DOI:10.1109/TGRS.2004.837336]
  • Li C, Xue Y, von Hoyningen-Huene W, Zhang J H, Pan P Y.2014.Quantitative evaluation of uncertainties in satellite retrieval of dust-like aerosols induced by spherical assumption. Journal of Quantitative Spectroscopy and Radiative Transfer, 136 : 45–57 . [DOI:10.1016/j.jqsrt.2014.01.003]
  • Li Z Q, Blarel L, Podvin T, Goloub P, Buis J P, Morel J P.2008.Transferring the calibration of direct solar irradiance to diffuse-sky radiance measurements for CIMEL sun-sky radiometers. Applied Optics, 47 (10) : 1368–1377 . [DOI:10.1364/AO.47.001368]
  • Liu K N.2002.An Introduction to Atmospheric Radiation. Vol. 84.2nd ed. San Diego:Academic Press Liu G R, Chen A J, Lin T H, Kuo T H. 2002a. Applying SPOT data to estimate the aerosol optical depth and air quality. Environmental Modelling and Software, 17 (1) : 3–9 . [DOI:10.1016/S1364-8152(01)00047-0]
  • Liu G R, Lin T H, Kuo T H.2002b.Estimation of aerosol optical depth by applying the optimal distance number to NOAA AVHRR data. Remote Sensing of Environment, 81 (2/3) : 247–252 . [DOI:10.1016/S0034-4257(02)00002-0]
  • Remer L A, Tanré D, Kaufman Y J, Ichoku C, Mattoo S, Levy R, Chu D A, Holben B, Dubovik O, Smirnov A, Martins J V, Li R R, Ahmad Z.2002.Validation of MODIS aerosol retrieval over ocean. Geophysical Research Letters, 29 (12) : MOD3-1–MOD3-4 . [DOI:10.1029/2001GL013204]
  • Simha C P, Devara P C S, Saha S K.2013.Aerosol pollution and its impact on regional climate during Holi festival inferred from ground-based and satellite remote sensing observations. Natural Hazards, 69 (1) : 889–903 . [DOI:10.1007/s11069-013-0743-6]
  • Sun L.2006.Remote Sensing of Aerosols over Urban Areas. Beijing:Institute of Remote Sensing Applications Chinese Academy of Sciences
  • Sun L, Liu Q H, Liu Q andChen L F.2006.Aerosol optical thickness retrieving over bright surface:progress and prospect. Progress in Geography, 25 (3) : 70–78 .
  • ( 孙林, 柳钦火, 刘强, 陈良富. 2006. 高反射率地区气溶胶光学厚度遥感反演:现状及展望. 地理科学进展, 25 (3) : 70–78. [DOI:10.11820/dlkxjz.2006.03.009] )
  • Tanré D, Kaufman Y J, Herman M, Mattoo S.1997.Remote sensing of aerosol properties over oceans using the MODIS/EOS spectral radiances. Journal of Geophysical Research:Atmospheres, 102 (D14) : 16971–16988 . [DOI:10.1029/96JD03437]
  • Tanré D, Herman M, Deschamps P Y.1981.Influence of the background contribution upon space measurements of ground reflectance. Applied Optics, 20 (20) : 3676–3684 . [DOI:10.1364/AO.20.003676]
  • Tanré D, Deschamps P Y, Devaux C, Herman M.1988.Estimation of Saharan aerosol optical thickness from blurring effects in thematic mapper data. Journal of Geophysical Research, 93 (D12) : 15955–15964 . [DOI:10.1029/JD093iD12p15955]
  • Vermote E F, Tanré D,Deuzé J L, Herman M and Morcrette J J. 1997.6S user guide Version 2. Lille:Laboratoired'Optique-Atmosphérique
  • Xu F F. 2012. AOD Retrieval Applying Structure Function Method Supported by Surface Reflectance Products.Qingdao:Shandong University of Science and Technology