The Image-Based Integrated Method for Determining and Mapping Aerosol Optical Thickness

Abstract

This chapter focuses on the development of an image-based integrated method for determining and mapping aerosol optical thickness (AOT). Using the radiative transfer (RT) equation, a methodology is developed using nonvariant targets with the revised darkest pixel method to improve the accuracy to calculate AOT over urban areas. This method can be applied to create thematic maps using GIS that present AOT values in a given area. The methodology is applied to Landsat TM/ETM+ satellite images of Limassol, Cyprus over a period of time to visualize and assess the AOT levels over the urban area. An accuracy assessment of the method shows a strong correlation between the in-situ AOT values from the sunphotometers and the AOT values derived from the image-based integrated method.

1. Introduction

Aerosol optical thickness (AOT) is defined as the degree to which aerosols prevent the transmission of light through absorption or scattering of light. AOT can be retrieved by using radiometers on board satellites [1–3], ground-based sunphotometers, and satellite sensors, such as the moderate resolution imaging spectroradiometer(MODIS), which provide AOT measurements on a daily and monthly basis [4, 5]. However, the spatial resolution of the MODIS AOT data products is 10 × 10 km, which does not permit the identification of specific AOT distribution over urban areas or complex terrains [5–7]. Due to the variability of aerosols, atmospheric aerosol monitoring is difficult. Significant efforts to improve aerosol characterizations have included using in-situ measurements, ground-based remote sensing, and satellite observations [1–3]. AOT derived from satellite images often requires further validation [8]. The accuracy of satellite-derived AOT is frequently assessed by comparing satellite-based AOT with AErosol RObotic NETwork (AERONET, a program established by NASA) or field-based sun photometers [9, 10].

Research indicates that the determined AOT from satellite image data can be used as a tool to assess air pollution as well as identify the sources of local emissions [2, 7, 11–24]. AOT values can be used as a way of measuring air quality; determining AOT in large-scale pollution areas provides a synoptic, cost-effective means to further assess the air quality in such areas [16, 18–24]. Indeed, the AOT values derived from the atmospheric path radiance can be utilized to assess and monitor air quality and atmospheric pollution [8]. The image-based integrated method presented in this chapter can accurately calculate AOT values retrieved from satellite imagery by using radiative transfer (RT) equations and GIS. The method can be used to visually display AOT levels using thematic maps in order to identify concentrations of AOT over an urban area [11]. The image-based method can also be used with archived satellite images, thereby providing detailed information regarding spatial aerosol concentration overtime [11].

AOT is directly related to the atmospheric aerosol load, which is the main variable describing the effects of aerosols on radiative transfer in the Earth’s atmosphere. According to Guanter et al. [25], modeling atmospheric constituents and surface reflectance involves modeling the radiative transfer across the atmosphere. The key parameter for assessing atmospheric pollution is the aerosol optical thickness, which is also the most important unknown of every atmospheric correction algorithm for solving the radiative transfer equation and removing atmospheric effects from remotely sensed satellite images [8, 26–29]. Several researchers [5, 18, 26, 28, 29–37] found that using the radiative transfer and atmospheric modeling in conjunction with field measurements of aerosol optical thickness can yield more accurate atmospheric corrections instead of using simple image-based techniques.

The image-based integrated method discussed in this chapter combines the RT equation, satellite imagery, the modified darkest pixel (DP) method of atmospheric correction, and GIS to derive AOT measures. An example including 11 Landsat satellite images with in-situ measurements over a specific period of time is used to assess the AOT values based on the 30 × 30 m spatial resolution of Landsat over the city of Limassol, Cyprus and create thematic maps to display the AOT levels.

2. Methodology

The innovation of this methodology is the retrieval of AOT by using an image-based algorithm based on the radiative transfer equation to derive AOT using the reflectance values from the darkest pixel method of atmospheric correction. In the image-based integrated method, the RT equation retrieves AOT values from the satellite images with as much precision as possible. Also, the integrated use of field spectroscopy, GIS, and remote sensing for AOT retrieval is an important contribution to the field of remote sensing.

The first step of the methodology is the preprocessing of the satellite images, including geometric and radiometric correction. Also, the satellite images need to be atmospherically corrected. The darkest pixel atmospheric correction method can be applied, as it provides reasonable atmospheric correction, especially in cloud-free skies [28, 38, 39]. If there is a known nonvariant target or overpass in-situ measurements available in the selected site where the true ground reflectance values of the darkest target are known, then the modified DP method can be applied to result in more accurate reflectance values to the corrected satellite image [40]. The images are then processed using the RT equations in order to retrieve the AOT levels. Specific parameters from the satellite images, sun photometer measurements, and the reflectance values from the atmospherically corrected images are incorporated. In this methodology, the AOT values are calculated at the 500 nm wavelength, so Landsat band 1 will be used to derive AOT levels through the algorithm. Also, GIS analysis is done to produce thematic maps showing the AOT levels in the area of interest. In this chapter, an accuracy assessment is performed in the example of Limassol, Cyprus, where the AOT values retrieved from the GIS analysis are compared with the in-situ sunphotometer measurements from the sunphotometers. Figure 1 shows the methodology for the AOT retrieval.

Figure 1.

Methodology for AOT retrieval.

2.1. Detailed description of the image-based integrated method for determining AOT

According to Hadjimitsis and Clayton [40], RT equations can be used to retrieve AOT values from the path radiance. Therefore, the algorithm used in the image-based integrated method is formulated to solve for AOT, which is included in the equation for path radiance (L_p) and Mie scattering [11]. To calculate AOT, Eq. (1) is used.

Path radiance is the sum of the atmospheric path radiance due to Rayleigh (L_pr) and Mie scattering (L_pa) [11, 29, 40–43]:

$L_p=L_{\text{pr}}+L_{\text{pa}}$E1

where L_p is the path radiance (integrated band radiance measured in W/m²sr), L_pr is the atmospheric path radiance due to Rayleigh scattering in W/m²sr [44], and L_pa is the atmospheric path radiance due to Mie scattering in W/m²sr [44].

2.1.1. Calculating path radiance

The path radiance “‘L_p’” measured in Wm⁻² μm⁻¹ sr⁻¹, is calculated according to Eq. (2) [42, 45]:

$L_p=L_{\text{ts}}-\frac{\tau \left(\mu \right)\uparrow {\rho }_{\text{tg}}{E}_G}{\pi },$E2

where L_ts is the measured at-satellite radiance, τ(μ) ↑ is the direct (upward) target sensor atmospheric transmittance, ρ_tg is the ground target reflectance, and E_G is the global irradiance reaching ground.

Eq. (3) can be used to calculate the at-satellite radiance (L_ts) for Landsat satellite images that contain gain and calibration offset values in the satellite image’s header file, as described in the Landsat-7 Handbook [46]:

$L_{\text{ts}}=\text{gain DN}+\text{offset}$E3

The “‘E_G” is calculated according to Eq.(4) [45]:

$E_G=E_{0\text{λ}}\mu \text{exp}\left[-\left(\frac{1}{2{\tau }_r}+\frac{1}{6{\tau }_a}\right)/\mu \right]$E4

where E_0λ is the solar irradiance at the top of the atmosphere in Wm⁻², found in the calibration file of the satellite image, μ is the cosine of the solar zenith angle (cosθ₀), found on the satellite image header file, τ_r is the Rayleigh optical thickness, and τ_a is the AOT (which is unknown).

The “‘τ(μ) ↑’” [42] is calculated in Eq.(5):

$t\left(\mu \right)\uparrow =\text{exp}\left[-\left({\tau }_r+{\tau }_a\right)/\mu \right]$E5

Although several equations to calculate τ_r have been developed, most radiative transfer equations used the equation provided by Moller [47], as indicated in Eq.(6), which is considered to be the most appropriate [34]:

${\tau }_r\left({\lambda }_c\right)=0.00879{\lambda }_c^{-4.09}$E6

where λ_C is the central wavelength for each Landsat band.

2.1.2. Calculating Rayleigh scattering (L_pr)

To calculate path radiance due to Rayleigh scattering(L_pr), Eq.(7) is used [44]:

$L_{\text{pr}}=\left\{\frac{E_{0\text{λ}}\text{cos}{\theta }_0{P}_r}{4\pi \left(\text{cos}{\theta }_0+\text{cos}{\theta }_{\nu }\right)}\right\}\left\{1-\text{exp}\left[-{\tau }_r\left(\frac{1}{\text{cos}{\theta }_0}+\frac{1}{\text{cos}{\theta }_{\nu }}\right)\right]\right\}t{0}_{3\uparrow }t{0}_{3\downarrow }$E7

where P_r is Rayleigh scattering phase function and θ_v is the sensor viewing angle, as found on the satellite image header file. to_3↑ refers to the transmittance factor due to ozone in the upward direction, and to_3↓ refers to the transmittance factor due to ozone in the downward direction, both of which are considered negligible and therefore have a value of 1 [44].

The Rayleigh scattering phase function (P_r) is calculated using Eq.(8) [44]:

$P_r=\frac{3}{4}\left[+{\text{cos}}^2\left(180-{\theta }_0\right)\right]$E8

2.1.3. Calculating Mie scattering

To calculate Mie scattering, Eq. (9) is used [47]:

$\begin{array}{l}{L}_{\text{pa}}=\left\{\frac{E_{0\text{λ}}\text{cos}{\theta }_0{\omega }_a{P}^{\text{'}}{}_a}{4\pi \text{cos}{\theta }_0+\text{cos}{\theta }_{\nu }}\right\}\left\{1-\text{exp}\left[-{\tau }_a\left(\left(\frac{1}{\text{cos}{\theta }_0}\right)+\left(\frac{1}{\text{cos}{\theta }_{\nu }}\right)\right)\right]\right\}\\ \text{exp}\left[-{\tau }_r\left(\left(\frac{1}{\text{cos}{\theta }_0}\right)+\left(\frac{1}{\text{cos}{\theta }_{\nu }}\right)\right)\right]\end{array}$E9

where w_a is the aerosol single scattering albedo as calculated by Waggoner et al. [48] and P′_a is the phase function, as calculated by Gilabert et al. [44], Forster [29], and Turner et al. [45], provided in graph form or downloaded from AERONET stations taken with the Cimel sunphotometer (http://aeronet.gsfc.nasa.gov/).

2.1.4. Steps for performing the image-based integrated method

In order to retrieve AOT using the above equation and solve for τ_a, the below steps should be followed:

Step 1: Download the required Landsat image from the USGS website (glovis.usgs.gov) by selecting the area of interest and the date and time of the satellite overpass.

Step 2: Retrieve the solar zenith angle (θ₀) at the overpass time, according to the header file on the satellite image.

Step 3: Retrieve the sensor viewing angle (θ_v) at the overpass time, according to the header file on the satellite image.

Step 4: Retrieve E₀ according to satellite image calibration file.

Step 5: Calculate μ, which is the cos(θ₀) using the variable derived from step 2.

Step 6: Retrieve λ_C, using satellite image calibration file.

Step 7: Calculate τ_r using Eq. (6), using the variable determined in step 6.

Step 8: Solve for E_G using Eq. (4), using the variables derived from steps 4, 5, and 7. The τ_a (AOT) will be unknown.

Step 9:Find t(μ)↑ using Eq. (5), using the variables derived from steps 5 and 7. The τ_a (AOT) will be unknown.

Step 10: Calculate L_ts from the satellite image header file using Eq. (3).

Step 11: Calculate ρ_tg from the surface reflectance values.

Step 12: Calculate L_p using Eq. (2), using variables from steps 8–11.

Step 13: Calculate P_r using Eq. (8), using variables from steps 2 and 5.

Step 14: Solve for L_pr using Eq. (7), using variables from steps 2, 3, 4, 7, and 13 and using 1 for to_3↑ and to_3↓.

Step 15: Enter w_a which is between 0.73 and 0.87 in urban areas and 0.89–0.95 in remote areas.

Step 16: Enter P′_a using phase function graph or from AERONET.

Step 17: Calculate L_pa using Eq. (9), using variables from steps 2, 3, 4, 7, 15, and 16.

Step 18: Solve for τ_a (AOT) using Eq. (1), using variables from steps 14 and 17.

2.2. Image-based integrated method for GIS analysis

In order to solve for AOT and conduct a GIS analysis, it may be necessary to simplify the equation due to the complexity of the RT equations and their logarithmic components. The equations for path radiance, Rayleigh scattering, and Mie scattering can be simplified and then combined into one equation, to be solved in an image-based software. After the AOT values are calculated for every pixel in the satellite image, the values can be exported into a GIS geospatial database, with the AOT value and coordinates. In order to display AOT on a GIS thematic map, the Kriging method can be used to interpolate the AOT values and then overlay the data into GIS vector maps or satellite images of the area of interest to display a spatial distribution of AOT. The maps can also be overlaid with vector data including roads, building lots, parks, green areas, stadiums, industrial areas, and municipal boundaries to identify AOT distribution over the area of interest. In-situ AOT values, if available within the area of interest, can be correlated with the AOT values from the GIS maps to perform an accuracy assessment for a specific location.

3. AOT retrieval using the image-based integrated method: the example of Limassol, Cyprus

The example described in this chapter includes AOT and surface reflectance field measurements in Limassol Cyprus during Landsat-5 TM and Landsat-7 ETM+ satellite overpass. Only 11 Landsat satellite images with low cloud cover are used. These images were chosen due to their availability and clarity. The Landsat-5 TM images used are for 24/6/2010, 10/7/2010, 27/8/2010, and 28/9/2010, and the Landsat-7 ETM+ images used are for 13/4/2010, 29/4/2010, 31/5/2010, 16/6/2010, 7/11/2010, 9/12/2010, and 2/5/2011. During each satellite overpass, AOT measurements were taken using the Microtops II sunphotometer at five different locations in Limassol and the stationery Cimel sunphotometer, which is part of the AERONET network, located in the center of the study area. These AOT values are later compared to the AOT values derived from the algorithm using the satellite images. At the same time, the SVC HR-1024 spectroradiometer was used to measure the reflectance value of a large nonvariant target dark asphalt surface located near the center of the satellite image.

Geometric and radiometric correction is performed on the 11 Landsat-5 and Landsat-7 satellite images, as well as atmospheric correction using the modified DP method. The image-based integrated method is applied to the Landsat-5 and Landsat-7 images taken over specific days of 2010–2011. Eqs.(1–9) are used to retrieve AOT from bands 1–4 of the Landsat-5 and Landsat-7 satellite images, following the steps outlined. As the in-situ AOT measurements are at the 500 nm wavelength, Landsat band 1 is used to derive AOT levels in this method.

In this example, the model is run in the ERDAS Imagine and MATLAB software to generate an image that is consists of AOT values. The model is then applied to every surface using the Landsat band 1 satellite image, since the AOT required are in the 500 nm wavelength. The resulting image features the AOT values for every pixel in the image, as featured in Figure 2. Areas with no data, due to cloud presence or negative AOT values, appear in white. All AOT values can be exported to ASCII files to create a georeferenced AOT dataset, which is then imported into the ArcGIS software.

Using the ArcGIS software, each pixel is converted to a point and each point is associated with the calculated AOT value. The white area of the image contains values of “no data” (Figure 3, left). To create a GIS thematic map showing the AOT distribution in the urban area of Limassol, the interpolation method is used in order to estimate the values from the “no data” regions. The ordinary Kriging interpolation tool is applied, where the unknown values of the “no data” area are interpolated using the weighted average of neighboring samples. Figure 3 (right) shows the results of the Kriging interpolation for the Limassol example, where the dots indicate the points with AOT values, while the colored sections are the resulting AOT values determined from Kriging interpolation.

Figure 2.

AOT values derived from algorithm.

Figure 3.

(Left) Points with AOT dataset. (Right) Interpolation with Kriging method.

Once the interpolation is completed, an AOT thematic map is created, which is classified according to AOT values. The AOT values are displayed in different colors according to the AOT concentration and a legend of the AOT values is created. This facilitates the ability to visualize the concentration and distribution of AOT values over an urban area. Figure 4 features a map of Limassol following Kriging interpolation, which provides synoptic views of the AOT distribution. As is evident, high AOT levels are present throughout the city, especially in the industrial estates and on busy streets, with lower AOT levels present in parks, stadiums, and outside the city.

Figure 4.

Thematic map of Limassol (09/12/2010) after Kriging interpolation, with AOT values.

In the example, the thematic map is overlaid with GIS vector data from the Lands and Surveys Department for Limassol to facilitate the identification of sources of AOT values within the area of interest (Figure 5). Each polygon in the thematic map is connected to a GIS database, where data such as plot number and area information can be identified. In this way, the GIS database, combined with the AOT data, can create thematic maps to illustrate trends in different areas.

Figure 5.

Vector overlay.

In the example, the image-based integrated method is used to produce GIS maps indicating the AOT distribution over Limassol. In order to visualize high and low AOT levels, thematic maps are generated using colors ranging from blue to green to yellow to red for the specified AOT range. A separate thematic map is created for each of the 11 Landsat-5 and Landsat-7 satellite images. Figure 6 features the GIS thematic map using Landsat-7 satellite imagery from 24/6/2010.

Figure 6.

GIS thematic map, 24 June 2010, with vector overlay.

The AOT values derived from the algorithm are compared with the AOT values measured with the Microtops and Cimel sun photometers. Figure 7 shows all the locations of the in-situ AOT measurements using sun photometers. Locations 1–5 are measured with the Microtops sunphotometer, while location 6 is measured with the Cimel sunphotometer. In Figure 7, location 7 refers to the Department of Labor Air Quality Monitoring Site, which was used only for reference. An accuracy assessment is done to compare the AOT from the GIS map with the in-situ AOT as measured from the Microtops and Cimel sunphotometers during satellite overpass. A linear regression is conducted of the AOT levels measured with the sunphotometers and the AOT values derived from the algorithm, with a coefficient of determination R2 of 0.977.

Table 1 compares the AOT measurements taken on site using a sunphotometer against the AOT values derived from the GIS maps using the proposed algorithm by location and date, as well as the correlation coefficient (r) between the on-site AOT measurements and algorithm-derived AOT values from GIS, for each location. In order to determine the accuracy of the GIS model, the root mean square deviation (RMSD) is also calculated between the on-site AOT measurements and algorithm-derived AOT values from GIS for each location. The results of the RMSD, correlation coefficients, and the coefficient of determination verify that the AOT derived from the GIS map correlate strongly with the on-site AOT values as measured by sunphotometers.

Figure 7.

Location of in-situ AOT measurements.

Table 1.

Comparison of on-site AOT and AOT from GIS algorithm (4/2010–5/2011).

4. Conclusions

In this chapter, the image-based integrated method is developed to retrieve AOT values from satellite images and display these values in an urban area, using the example of Limassol, Cyprus. The findings show that the proposed integrated image-based method is able to accurately obtain AOT measurements from satellite images. The accuracy assessment shows strong agreement between the in-situ AOT values and the AOT values retrieved through the method described in the chapter. The proposed methodology is able to retrieve AOT using satellite images and visually display AOT over urban areas by using GIS to produce thematic maps. Also, the significance of the image-based integrated method is that the environmental data is not needed, as compared with other methods. This methodology is therefore an alternative to more sophisticated and complex methods of deriving AOT values from satellite images, especially where atmospheric or meteorological data are not accessible. Also, the methodology can be used with archived satellite images where environmental data are not available, thereby providing detailed information regarding spatial aerosol concentration overtime. Also, thematic maps can be generated to illustrate the distribution of AOT overtime in a specific area of interest and enabling trends to be identified.

Acknowledgments

The author would like to thank the Cyprus University of Technology, Department of Civil Engineering and Geomatics, Remote Sensing Lab for use of the sunphotometers and AERONET data.