An Improved Aerosol Optical Depth Retrieval Algorithm for Moderate to High Spatial Resolution Optical Remotely Sensed Imagery

Abstract

To extract quantitative land information accurately and monitor the air pollution at city scale from moderate to high spatial resolution (MHSR) with a resolution no coarser than 30 m, optical remotely sensed imagery and aerosol parameters, especially aerosol optical depth (AOD), are a necessary step. In this paper, we introduce a new algorithm that can effectively estimate the spatial distribution of atmospheric aerosols and retrieve surface reflectance from moderate to high spatial resolution imagery under general atmosphere and land surface conditions. This algorithm has been improved in the following three aspects: (i) it has been developed for most of the moderate to high spatial resolution remotely sensed imagery; (ii) it can be applied to all kinds of land surface types including bright surface; and (iii) it is completely automatic. This algorithm is therefore suitable for operational applications. The derived AOD in Beijing from Landsat Thematic Mapper (TM), Landsat Enhanced Thematic Mapper Plus (ETM+), and Huan Jing 1 (HJ-1/CCD) data is validated with AErosol Robotic NETwork (AERONET) ground measurements at Beijng and Xianghe stations.

1. Introduction

A large number of land resource survey satellites has been launched in recent years. These satellites are usually loaded with optical sensors with moderate to high spatial resolution (MHSR), with a resolution no coarser than 30 m, and provide an enormous amount of earth observing data. These data have been extensively used for agricultural applications, forest management, geological surveys, water resource evaluation, disaster appraisals, and many other applications [1]. Instead of TOA reflectance from the sensor, surface reflectance is required to make consistent and quantitative geophysical and biophysical products due to the variable effect of atmosphere, especially aerosols, on optical wavelength in space and time [2,3]. Since data with a higher spatial resolution usually have lower temporal frequency and most of the data are contaminated by aerosols and other factors, only a limited number of data can be used without atmospheric correction. As the quantitative utility of these satellite data become more demanding, surface reflectance becomes increasingly important [2]. Atmospheric correction includes two major steps [2]: parameter estimation and surface reflectance retrieval. Once all atmospheric parameters are known, the retrieval of surface reflectance is relatively straightforward; therefore, aerosol retrieval is the key step.

Aerosol plays an important role in the earth’s radiation budget (ERB) by scattering and absorption of solar radiation, consequently affecting plant photosynthesis levels [4]. To quantitatively estimate the influence of aerosols on ERB, many studies on retrieving AOD from moderate to low spatial resolution remotely sensed data like MODerate resolution Imaging Spectroradiometer (MODIS) have been conducted. These studies include the densely dark vegetation (DDV) method [5,6], deep blue [7], second-generation operational algorithm [8], synergy of Terra and Aqua MODIS data [9], minimum reflectance technique [10], simplified high resolution MODIS aerosol retrieval algorithm [11], and the algorithm based on long time series MODIS imagery [12,13]. However, the retrieved AOD from moderate to low spatial resolution data is more suitable for global or regional scale aerosol spatial distribution studies [14].

Since aerosols are mainly composed of particulate matter (PM), they and what they have carried with them are the main pollutants influencing human health [15]. These pollutants are mainly induced by human activities in urban areas, so many studies focus on aerosol retrieval at city scale from MHSR remotely sensed data. The most used method is the dark object (DO) method [16] including DDV [5], and this method can only be applied to partial areas of images with dense vegetation, which is a serious limitation to imagery acquired over the winter season or in arid areas. The contrast reduction method [17] is used for stable surfaces with multi-temporal imagery and this method uses the variation of aerosol over a stable surface to extract AOD. The assumption of invariant surface reflectance limits its general applications because surface reflectance changes in both space and time under general conditions [1]. The histogram matching method [18] assumes that the surface reflectance histograms of clear and hazy regions are the same. After identifying clear sectors, the histograms of hazy regions are shifted to match the histograms of their reference sectors (clear regions) [18]. The difficulty of this method lies in how to identify the clear regions, which usually requires manual intervention [1]. There are also some hybrid methods such as ImAero-Landsat [19], which incorporates DDV and contrast reduction methods. However, it still needs a clear image as a reference for the contrast reduction method, or it fails.

With more MHSR remotely sensed imagery freely available, the potential to produce continental (and even global) scale products at a higher resolution has been launched, such as the NASA Web-enabled Landsat Data (WELD) project [20], and the Finer Resolution Observation and Monitoring of Global Land Cover (FROM-GLC) project [21]. Consequently, systematic atmospheric correction of MHSR data at the continental or even global scale is challenging due to the large data volume and the requirement of consistency among different scenes of imagery. The Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm [22] and MODIS-based algorithms [23,24] are proposed for this purpose. The LEDAPS uses the DDV method to derive the initial AOD and the void areas are interpolated by averaging neighboring values within a 7-km window, which will greatly degrade the effects of atmospheric correction. The MODIS-based algorithm uses the MODIS AOD product with a 10-km resolution, which ignores the local variation of aerosols.

With increased availability of high spatial resolution satellite imagery, such as Quickbird, Ikonos, and GeoEye, with spatial resolutions ranging from a few meters to one meter, urban areas have become one of the research hotspots of remote sensing [25], which usually requires a spatial resolution finer than five meters to improve the identification of buildings, street networks, and sparse vegetation [26]. Consequently, as one of the key steps for remote sensing preprocessing, the atmosphere characterization retrieval, especially the estimation of aerosol properties of urban areas has become a prerequisite for urban remote sensing research. OSIS (Observation of Shadows for aerosols Inversion over 3D Scenes) [25] was developed to retrieve the characterization of aerosol optical properties directly from very high spatial resolution images over urban areas, and employed a technique based on shadow/sun transitions [27]. OSIS was also applied to multispectral PELICAN airborne data, where it showed a high accuracy of aerosol optical thickness retrieval, even in complex areas such as cities [28]. However, the use of OSIS was limited by three aspects: the need to manually select the pixel couples for aerosol retrieval, the need to use the vector elevation model, and the long computation time [25]. Thus, it is very difficult to use for operational purposes.

In order to automatically retrieve accurate AOD over general atmosphere and land surface conditions from MHSR remotely sensed imagery, we present a new AOD estimation algorithm in this study. In this method, we assume that the atmosphere varies continuously in space, which is reasonable in practice. Consequently, the spatial distribution of aerosol can be derived through spatial expansion from a portion of known AODs, which can be accurately retrieved using the DO method. This algorithm is then successfully applied to the major MHSR remotely sensed imagery to retrieve AODs and conduct atmospheric correction. The estimated AODs in Beijing from (Enhanced) Thematic Mapper (TM) Plus (ETM+) and HJ-1/CCD are validated by Aerosol Robotic Network (AERONET) measurements, respectively. The two validation cases indicate that this algorithm can retrieve AODs from MHSR imagery reasonably well.

2. Data and Methodology

2.1. MHSR Remotely Sensed Data

In this study, a new algorithm was developed for the major MHSR remotely sensed data. These data have a spatial resolution not lower than 30 m including Landsat-TM/ETM+, Operational Land Imager onboard Landsat8 (Landsat/OLI), HJ-1/CCD, Wide Field Viewer onboard Gao Fen 1 (GF-1/WFV), Multispectral Imager on board ZiYuan 3 (ZY-1/MSI), Multispectral Imager on board GeoEye-1 (GeoEye-1/MSI), CCD and Wide Field Imager onboard Chinese and Brazil Earth Resource Satellite 02B (CBERS-02B/CCD/WFI), and so on. With the development of this algorithm, more MHSR remotely sensed data will be added. The detailed information of these data is listed in

Table 1. The detailed information of the major MHSR remotely sensed data used in this study.

Sensor(s)	Platform(s)	Spatial Resolution	Swath	Spectra for AOD Retrieval	Nation
TM, ETM+, OLI	Landsat series	30 m	185 km	VIS *, NIR †, SWIR ‡	USA
CCD	HJ-1/A&B	30 m	700 km (2 cameras)	VIS, NIR	China
CCD HR	CBERS-02B	20 m 2.36 m	113 km 27 km	VIS, NIR	China
CCD WFV	GF-1	8 m 16 m	60 km 800 km (4 cameras)	VIS, NIR	China
MSI	ZY-3	5.8 m	52 km	VIS, NIR	China
MSI	GeoEye-1	1.65 m	15.2 km	VIS, NIR	USA
MSI	THEOS	15 m	90 km	VIS, NIR	Thailand

* VIS means visible bands; ^† NIR means near infrared band; ^‡ SWIR means shortwave infrared band.

. All the data have visible and near infrared bands and only the Landsat series data have shortwave infrared bands, which are usually used for the DDV method. The spatial resolution ranges from 30 m to 1.65 m and in this study they were divided into two categories: moderate high resolution (coarser than 10 m and finer than 30 m, MHR) and high resolution (not coarser than 5 m, HR).

2.2. The Procedure

This new algorithm is a combination of the DO method [16] and the histogram matching method introduced by Liang [2]. Instead of determining the clear and hazy regions of an image, the DO method was used to obtain the initial AOD and the initial clear regions. Thus, the difficulty in determining the clear and hazy regions were avoided. Since most of the initial AOD regions are located in mountainous areas, especially in winter seasons and at arid areas, the AODs from other regions need to be retrieved through the spatial expansion method. The spatial expansion method is based on the assumption that the atmosphere, including aerosols, varies continuously in space. The procedure is illustrated in Figure 1, and the major steps are discussed in detail below.

Since longer wavelengths are less contaminated by aerosols, an unsupervised classification using NIR and two SWIR bands for the Landsat series of data, or Red and NIR bands for the other data was conducted to obtain the classification map of the input image. In this study, the K-means method was used, but other classification methods such as IsoData could also be employed. The number of classes needed to be specified in the first place, and the number of classes usually depends on the complexity of the landscape. Liang et al. indicated that 20–50 classes produced reasonable results [2]; therefore, we used the maximum classes of 50 in this study. The average reflectance of each class under different atmospheric conditions (from clear to hazy) was used for histogram matching.

In order to identify the clear regions of an image, the DO method combined with spatial expansion was used to derive the initial clear regions instead of drawing the clear regions manually [2]. The DO method was slightly modified in this study to satisfy multiple MHSR remotely sensed data, and is detailed in Section 2.3. The spatial expansion process is employed to ensure that the clear regions are large enough to incorporate more classes, and the detailed information is described in Section 2.4.

After the clear and hazy regions were determined, the histograms of the average reflectance of each class between clear and hazy regions were matched to obtain the surface reflectance of hazy regions. Given the surface reflectance, the AOD of hazy regions can be determined by using the Look-Up Table (LUT) method, which is described in Section 2.5.

Since the initial clear regions were derived from the DO method, the area of initial clear regions may be very small so that they do not include many of the classes of the input image, especially in winter seasons and in arid regions. Therefore, an iterated procedure was fulfilled before the last step. The area ratio of 90% between the coverage of the retrieved AOD and the whole image was used as the threshold. Once the threshold was satisfied, the iterated procedure stopped and went to the last step.

Finally, a spatial low-pass filter with a 5 × 5 window size was used to fill the regions without AOD retrieval, as aerosols are usually distributed smoothly in space, and the regions with less than 10% of the total area were randomly distributed over the image, which hardly influenced the final result.

2.3. The Modified DO Method

The DO method is probably the most widely-used method for AOD retrieval and atmospheric correction and has a long history [16,29,30,31,32,33]. Densely dark vegetation (DDV) is the most popular method within the DO methods, where image scenes containing dense vegetation pixels can be identified by shortwave infrared bands (around 2.1 μm), and the reflectance of this band has strong correlations with the reflectance at blue and red bands, from which the AOD of the image can be subsequently derived. However, for MHSR remotely sensed imagery, DDV can only be applied to Landsat series imagery (including TM, ETM+, and OLI) as other data do not include the shortwave infrared band. Moreover, even for Landsat series imagery, DDV becomes invalid under several situations: (1) in the winter seasons of high latitudes, and (2) in arid regions. For high spatial resolution imagery with narrow swathes such as GeoEye-1/MSI and GF-1/CCD, the possibility of dense vegetation missing is very high. Thus, the DO method was slightly modified in this study.

First, besides the dense vegetation, more dark objects were found, including shadows induced by mountains, clouds, and buildings. Mountain and cloud shadows are usually used by most MHSR data. Building shadows are only used by HSR data such as GF-1/CCD and GeoEye-1/MSI, in which building shadows can be identified. Since the surface reflectance at the blue band in mountainous areas is usually very low, the TOA radiance is mainly the path radiance induced by atmosphere, which is the same as the ordinary DO method.

The surface reflectance of shadows induced by clouds and buildings is not necessarily low, so the portion reflected by the surface needs to be considered. For a flat and Lambertian surface under a horizontally homogeneous atmosphere, the top-of-atmosphere (TOA) radiance L can be expressed by the following formula [31]:

$$L\left({\mu }_v\right)=L_p\left({\mu }_v\right)+E\times T\left({\mu }_v\right)\frac{\rho }{\pi \left(1-s\rho \right)}$$

(1)

where $L_p$ is the path radiance; $\rho$ is the surface reflectance; s is the spherical albedo of the atmosphere; $E$ is the solar irradiance at the bottom of atmosphere. T is the total transmittance from the surface to the sensor, and ${\mu }_v$is the cosine value of viewing zenith angles.

The total irradiance is the sum of direct and diffused irradiance, which can be expressed as follows:

$$E=E_{dir}+E_{dif}$$

(2)

where $E_{dir}$ and $E_{dif}$ are the direct and diffused portions of irradiance respectively.

For shadow areas, the direct irradiance is 0, so the TOA radiance $L_s$ is expressed as:

$$L_s\left({\mu }_v\right)=L_p\left({\mu }_v\right)+\frac{\Omega }{2\pi }{E}_{dif}\times T\left({\mu }_v\right)\frac{\rho }{\pi \left(1-s\rho \right)}$$

(3)

where $\mathsf{\Omega }$ is the sky-viewing solid angle, which means that the sky can be seen within this solid angle and is shaded out of this solid angle, and is expressed in sr. If there is no shaded area, the sky-viewing solid angle is 2π. Since cloud shadows are not shaded by other objects, $\mathsf{\Omega }$ is set to 2π and Equation (3) can be simplified as Equation (4). The shadow induced by a building is always accompanied by the walls of the building, so $\mathsf{\Omega }$ is usually less than 2π. $\mathsf{\Omega }$ can be calculated using a 3D model and the calculating method is referred to in Reference [28]. However, 3D models of cities are not readily accessible everywhere and the calculation is computational expensive; consequently, it is difficult for use in operational purposes. In practice, the shadow area close to the sun/shadow boundary was chosen and the influence of $\mathsf{\Omega }$ was lowered; subsequently, Equation (3) was simplified as Equation (4).

$$L_s\left({\mu }_v\right)=L_p\left({\mu }_v\right)+E_{dif}\times T\left({\mu }_v\right)\frac{\rho }{\pi \left(1-s\rho \right)}$$

(4)

By moving $L_p$ to the left and using Equation (4) to divide Equation (1), Equation (5) is derived:

$$\frac{L_s\left({\mu }_v\right)-L_p\left({\mu }_v\right)}{L\left({\mu }_v\right)-L_p\left({\mu }_v\right)}=\frac{ E_{dif}}{E}$$

(5)

$E_{dif}/E$ is the ratio of diffused irradiance within total irradiance, which can be derived from the output of 6SV (Second Simulation of a Satellite Signal in the Solar Spectrum–Vector) [34]. If the aerosol model is determined, $L_p$ and $E_{dif}/E$ are the monotonically increasing functions of AOD; consequently, AOD is derived using Equation (5). Therefore, Equation (5) was used to retrieve the AOD of shadow areas with high surface reflectance.

In practice, shadow detection is a prerequisite step for the DO method. The major methods of shadow detection were surveyed and categorized into six classes in Reference [32]: histogram thresholding, invariant color models, object segmentation, geometrical methods, physics-based methods, unsupervised, and supervised machine learning methods. Among these methods, the histogram thresholding on RGB and NIR channels [33] performed the best with an average accuracy of 92.5%. Moreover, the histogram thresholding needs few user inputs; therefore, the histogram thresholding method was chosen for shadow detection in this paper.

After the shadow areas were determined, a 9 × 9 window was used to scan the whole image. Only if the shadow area and non-shadow area were both not less than 1/3 of the 9 × 9 window (27 pixels), the average TOA radiance of the shadow areas and non-shadow areas (the boundary pixels between shadow and non-shadow areas were not used for averaging as they are partially influenced by direct solar irradiance) were used as $L$ and $L_s$ in Equation (5). The AOD was subsequently derived from the 9 × 9 window.

Second, most of MHSR data (except for the Landsat series of data) do not have shortwave infrared bands (around 2.1 μm); DDV is consequently not effective. In order to use dense vegetation to extract AOD from MHSR data without shortwave infrared bands, the following criteria were used:

$${\rho }_b=\{\begin{array}{cc}\hfill 0.02,& NDVI\ge 0.8\hfill \\ \hfill 0.06-0.05*NDVI,& 0.6\le NDVI<0.8\hfill \end{array}$$

(6)

where ${\rho }_b$ is the surface reflectance of the blue band and NDVI is the normalized difference vegetation index, which is defined as follows:

$$NDVI=\frac{{\rho }_{nir}-{\rho }_r}{{\rho }_{nir}+{\rho }_r}$$

(7)

where ${\rho }_r$ and ${\rho }_{nir}$ are the reflectance of the red and near infrared bands. The reflectance here is the TOA reflectance which is affected by atmosphere, so NDVI $\ge 0.6$ was used to assure that the selected pixels are dense vegetation and the densest vegetation is defined by $NDVI\ge 0.8$. The criteria were derived from the statistics from the MODIS data of northern China in July 2012 and the detailed procedure is as follows:

(1)

Calculate the NDVI of MODIS data using Equation (7).

(2)

Based on the criteria defined in the DDV method, use the MODIS band 7 (2.1 µa) to calculate the surface reflectance of the blue band.

(3)

Choose the NDVI and surface reflectance with the best quality, which are out of cloud and aerosol contamination.

(4)

Utilize the regression analysis between NDVI and the blue band reflectance to obtain the criteria in Equation (6).

Based on the two modifications, the DO method can be applied to more MHSR data where more patches of AOD can be retrieved, which is greatly helpful to the following processes.

2.4. Spatial Extension

Since atmosphere, including aerosols, is usually distributed smoothly in space, we can therefore assume that aerosols vary continuously and gradually in space. Subsequently, the patches of AOD obtained from the DO method in Section 2.2 were used as the initial values, and a spatial expansion process was used to expand the AOD patches within a reasonable range, in which the AOD of expanding areas can be derived by extrapolating from the initial AOD patches with little accuracy loss. Based on the assumption above, if there are some patches with known AOD in an image, the AOD over areas very close to the known AOD patches can be extrapolated. In this study, we defined the pixels within a 25-pixel distance to the boundary of the known AOD patches as the expanded areas, whose AOD can be extrapolated from known AOD patches. After the AODs were extrapolated, the surface reflectance of the expanded areas was derived by atmospheric correction. Figure 2 is an illustration of the spatial expansion process.

The left image shows an example of northern China in the non-growing season, in which most of the agricultural areas and urban areas are not densely vegetated, so the shadows and the dense vegetation in the mountainous areas (the areas filled with blue lines on the central image in Figure 2) were used as dark objects and the initial patches of AOD were retrieved (the central image in Figure 2 and an example of the initial AOD is shown in Figure 3b). However, the initial patches were almost in mountainous areas and the classes were very limited, so the histogram matching process could not be applied. With a 25-pixel expanding area from the initial patches, a large number of pixels in the agricultural areas (the areas contoured with yellow lines on the right image of Figure 2) were included and the histogram matching process subsequently worked effectively and efficiently. This process could be done iteratively until the entire image was filled with retrieved AODs. Figure 3c,d show two intermediate results of spatial expansion.

Figure 3 shows an example of the procedure of AOD expansion. Figure 3a shows the original image with a true color composite. Figure 3b is the initial AOD retrieved using the DO method and the areas are mainly mountains. Figure 3c is the first expansion with 25 pixels and the yellow depicts the expanded areas. Figure 3d shows the classification map using the K-means method. Figure 3e is the newly retrieved AODs from the histogram matching method based on the first expanded areas and the classification map; the greys are the newly retrieved AODs. Figure 3f is the final AOD map after several iterations and spatial low-pass filtering. Figure 3g is the land-cover classes used for the histogram matching.

In addition, while expanding the AODs, the cubic interpolation function of Matlab (griddata3 function) was used instead of simple extrapolation. The griddata3 function in Matlab first found the discrete points based on the three features around the insert point by the Delaunay method. A triangle was then created by the three discrete points. Finally, the linear interpolation/extrapolation or the cubic equations interpolation were chosen to calculate the value of any insert point within the triangle. This method took the initial patches and the expanding areas as a whole and could increase the accuracy of the interpolation/extrapolation.

This assumption is reasonable for most situations because the components of atmosphere are always moving from high concentration to low concentration and are eventually expected to be homogeneous spatially. The assumption is invalid for several scenarios as follows:

(1)

The elevation changes suddenly and dramatically and the height difference exceeds the aerosols spreading height so that aerosols cannot move towards areas with higher elevation. This scenario only happens in mountainous areas where human activities are infrequent and aerosols are usually low and homogeneous. Thus, its influence is neglected.

(2)

Close to pollution sources, such as industrial emission and large forest fires. The areas with heavy aerosols in this scenario usually cover very small areas of an image, so it does not affect the AOD retrieval unless the land surface simulated from neighbor areas is used to retrieve the AOD for these areas.

(3)

Combined with the classification map from unsupervised classification in Section 2.1, the spatial expansion process can be conducted several times until 90% of the pixels are filled with valid AOD values.

2.5. Creation of the LUTs

In the procedure, the AOD retrieval and atmospheric correction are employed very often and both of the processes need to resolve the radiative transfer equation. Instead of resolving the atmosphere radiative transfer equation online, which is very time consuming, the Look-Up Table (LUT) method is employed. The LUT method calculates different atmospheric quantities offline and is organized in the form of tables. These quantities can be determined online by searching the tables in operational procedures. In this study, the 6SV [34] was used to make the LUT. The following parameters were considered while constructing the LUT:

(1)

Solar and viewing geometry. The geometrical relationships among sun, target, and sensor can be described using the solar zenith, solar azimuth, sensor zenith, and sensor azimuth angles. The difference of solar and sensor azimuth angles is the relative azimuth angle. The geometry settings in the LUT were as follows: seven sensor zenith angles with a step of five between 0–30°; 10 solar zenith angles between 0° and 50°; and 19 relative azimuth angles with a step of 10 in the range 0–180°.

(2)

Atmospheric model. Mid-latitude summer, mid-latitude winter, tropical, subarctic summer, and subarctic winter were used, which are determined using latitude only.

(3)

Aerosol model. Urban, desert and continental models were used for urban, desert and other land surfaces, respectively. The aerosol model was determined using a global land cover map from the Finer Resolution Observation and Monitoring of Global Land Cover (FROM-GLC) project [21].

(4)

AOD at 550 nm. The setting was 0.01, 0.05, 0.1, 0.2, 0.4, 0.8, 1.0, 1.5, and 2.0.

(5)

Altitude. The setting was 0, 100, 200, 500, 1000, 2000, 3000, 4000, and 5000 m. In this study, the altitude was retrieved from the ASTER digital elevation model [35].

(6)

Water vapor. The setting was 0.1, 1.1, 2.1, 3.1, 4.1, and 5.1 g/cm². The water vapor was derived directly from the MODIS near-infrared water vapor bands (typical accuracy 5–10%) [36].

(7)

Surface reflectance. Since this study was for AOD retrieval, the blue band was the only one used for atmospheric correction and AOD retrieval in the procedure; therefore, the maximum surface reflectance of 0.2 for the blue band was set. Lambertian/isotropic land surface was assumed and 11 settings of surface reflectance between 0.0 and 0.2 with an increment of 0.02 were used.

3. Results and Validation

3.1. AOD Retrieval and Atmospheric Correction for Multiple MHSR Remotely Sensed Data

In this study, the new method was applied to a large amount of the MHSR remotely sensed data listed in Table 1. In order to demonstrate the effectiveness of the new algorithm, two examples are given: (1) the retrieved AODs and atmospheric corrected images for different types of MHSR imagery are shown in Figure 4. The multiple MHSR remotely sensed imagery in Figure 4 includes GF-1/WFV (a), Landsat/TM (b), HJ-1/CCD (c), Landsat8/OLI (d), ZY-3/MSI (e), THEOS/MSI (f), and GEOEYE-1/MSI (g). The images in Figure 4 also demonstrate the AOD retrieval and atmospheric correction over different land surface types including mountain (a–f), desert or Gobi (e and f), cropland (g) in both summer (d) and winter (e), bare soil (c, e and f), and urban (a and c). In addition, the atmospheric conditions of images in Figure 4 are very different. Therefore, Figure 4 comprehensively shows that the new algorithm has the ability to retrieve the AOD and surface reflectance for most of the MHSR remotely sensed imagery under general atmospheric and land surface conditions. (2) The retrieved AODs and atmospheric corrected images for multi-temporal Landsat TM/ETM+ imagery around Beijing are shown in Figure 5. The spatial distribution of the retrieved AODs (the right column of Figure 5) for every image in Figure 5 is consistent with the TOA reflectance (the left column of Figure 4) and the image after atmospheric correction (the central column of Figure 5) shows even tone throughout the imagery. The AODs retrieved at all the dates were also consistent. Therefore, Figure 5 shows the consistency of the retrieved AODs and the atmospheric correction effect from the new algorithm over time and space. The image taken on 25 May 2003 is not fully corrected because the central part of the image is too hazy (the land surface information is almost completely covered by aerosols) to be corrected. For the image at the top of Figure 5, the upper left portion of the image is not well corrected after atmospheric correction because the AOD is too high to be corrected.

Based on the examples shown in Figure 4 and Figure 5 and the analysis from over 10,000 results of MHSR remotely sensed imagery, the following conclusions are made:

(1)

The new algorithm can be applied to the major MHSR remotely sensed imagery.

(2)

The new algorithm can be applied at general atmosphere and land surface conditions.

(3)

The spatial distribution of the retrieved AOD is visually consistent with the distribution of aerosol, which is obvious in TOA reflectance imagery.

(4)

The retrieved AODs for multi-temporal Landsat TM/ETM+ imagery are also consistent, which shows the consistency of the new algorithm over time.

(5)

The effect of atmospheric correction for the major MHSR remotely sensed imagery is visually obvious. The correction effect for spatial consistency is very good, such as in Figure 3a,c,g, whose uneven distribution of aerosol are mostly corrected.

3.2. Validation

Although the new algorithm can be applied to the major MHSR remotely sensed imagery, there is insufficient imagery except for Landsat TM/ETM+ and HJ-1/CCD. Thus, to evaluate the performance of this algorithm quantitatively, AOD values at 550 nm wavelength derived from Landsat TM/ETM+ and HJ-1/CCD imagery at 30 m resolution around the Beijng area were compared with ground measured AOD data from AERONET sites including Bejing (39.976890°N, 116.381370°E) and Xianghe (39.753600°N, 116.961500°E), which are marked as red dots on Figure 4.

For Landsat TM/ETM+ imagery, since the temporal frequency was low (once per 16 days), only 14 images from 2001 to 2009 were retrieved. A total of 28 AERONET measurements from Beijing and Xianghe AERONET sites were collected and compared with the retrieved AODs, which are shown in the first plot of Figure 6. Results showed that the retrieved AODs agreed very well with the AERONET measurements and the R² value was as high as 0.949. The AOD discrepancies between those measured and those retrieved varied between −0.02 and 0.34, and around 57% of the derived AOD values correlated with AERONET AODs at low discrepancies (less than 0.15). The root mean square error (RMSE) was 0.184.

Since HJ-1/CCD imagery had a higher temporal frequency (once per two days), 84 images are selected from 2011 to 2013. Although some were influenced by small clouds and cloud shadows, about 70 AERONET measurements from the Xianghe and Beijing AERONET sites, respectively, were collected and compared with the retrieved AODs, which are shown in the second and third plots of Figure 6, respectively. Results showed that the retrieved AODs also agreed well with the AERONET measurements and the R² values for Xianghe and Beijing were 0.816 and 0.807, respectively, which were lower than that of Landsat TM/ETM+. The AOD discrepancies between those measured and those retrieved varied between −0.14 and 0.31, and around 50% of the derived AOD values correlated with AERONET AODs at low discrepancies (less than 0.15). The RMSEs for Xianghe and Beijing were 0.179 and 0.207, respectively.

The validation above indicates that AODs derived from the new algorithm agreed well with AERONET ground measured AODs. The AODs retrieved from Landsat TM/ETM+ were a little better than those from HJ-1/CCD due to three reasons: (1) the image quality of Landsat TM/ETM+ was better than that of HJ-1/CCD, especially the stability of retrieved radiance [37]; (2) HJ-1/CCD was short of 2.1 μm band, which could have improved the accuracy of the initial AOD retrieval; (3) HJ-1/CCD had much wider coverage than Landsat TM/ETM+, so the larger view zenith angle could also have influenced the accuracy of the AOD retrieval.

4. Conclusions

In this paper, a new algorithm that can effectively estimate the spatial distribution of atmospheric aerosols from the major MHSR imagery under general atmosphere and land surface conditions was developed. This algorithm was improved in the following four aspects: (i) it was developed for most of the MHSR remotely sensed imagery; (ii) it can be applied to all kinds of land surface types including bright surface; (iii) it is completely automatic; and (iv) the spatial resolution is the same as the input imagery and it can better present the spatial distribution of the AODs. This algorithm is therefore suitable for operational applications. The derived AOD around the Beijing area was validated by a comparison with AERONET measurements at Beijng and Xianghe sites, which indicated that the AOD data derived from the new algorithm agreed well with AERONET ground measured AODs. The AODs retrieved from Landsat TM/ETM+ were usually a little better than those from the other data like HJ-1/CCD, which was short of shortwave infrared bands.

Since the AOD retrieval depended on the initial clear regions, and because the clear regions for some of images were very small, especially in the winter seasons of higher altitudes and in arid areas, the spatial expansion procedure needed to be done iteratively and this procedure was very time consuming; thus, the efficiency of the algorithm was low in this situation. Operationally, the efficiency of the algorithm with small clear regions needs to be improved.