1. Introduction
Xinjiang province in northwest China is part of the Central Asian dust storm area, which is one of the main sources of dust aerosols [
1]. In recent years, many Xinjiang cities, especially the capital, Urumqi, have suffered a severe deterioration in air quality with significant contributions from atmospheric particulates [
2]. Aerosols can significantly influence the ecosystem, climate and hydrological cycle by affecting radiative forcing [
3] and its relation with the air quality indicators proportionated for sustainable development [
4]. High aerosol pollution events have a wide-ranging impact on visibility [
5] and human health [
6]. Scientific data about the spatial and temporal dynamic of dust aerosol in Xinjiang are needed by the local government to facilitate development of policies to protect the ecosystem and diminish dust storms. Due to the lack of ground stations in this vast area, satellite retrieval is the practical way to provide the spatial and temporal distribution of aerosol optical depth (AOD) [
7,
8,
9].
Many satellite sensors have released AOD products, including the total ozone mapping spectrometer (TOMS) [
10], geostationary operational environmental satellite (GOES) [
11], ozone monitoring instrument (OMI) [
12], medium resolution imaging spectroradiometer (MERIS) [
13], advanced very high resolution radiometer (AVHRR) [
14], multi-angle imaging spectroradiometer (MISR) [
15], sea-viewing wide field-of-view sensor (SeaWiFS) [
16], moderate-resolution imaging spectroradiometer (MODIS) [
17], and visible infrared imaging radiometer suite (VIIRS) [
18]. Yet, the quality of AOD products over arid/semiarid areas, such as Xinjiang, is relatively low due to a large bias in the surface reflectance estimation as well as the aerosol model used in the retrieval algorithms. The operational MODIS AOD product over land is based on two algorithms: the dark target (DT) and deep blue (DB) algorithms [
19]. In the Xinjiang area, the MODIS AOD product is mostly retrieved with the DB algorithm which is applied over bright areas, where the surface reflectance is relatively high, and distinguishing atmospheric aerosol contributions from the satellite sensor energy is difficult.
In the DB algorithm, for arid and semiarid regions, the surface reflectance was determined based on a pre-calculated surface reflectance database, which was compiled based on the minimum reflectivity method at the resolution of 0.1° × 0.1° for each season using MODIS images [
20]. The derived surface reflectance database therefore depends on the scattering angle, normalized difference vegetation index (NDVI), and season. Most of the validation studies concluded that, in general, MODIS DB retrieved aerosol products were comparable to aerosol robotic network (AERONET) data, and an expected error (EE) envelope could be defined that contained approximately 50–70% of the matchups [
21]. Bilal and Nichol reported up to 75–80% of the DB retrievals within the EE [
22]. One of the major error sources for the DB algorithm is the difference between the surface reflectance corresponding to the images and that from the pre-calculated database, as a result of the anisotropic surface reflectance. Another crucial aspect in AOD retrieval is the aerosol model. In the DB algorithm, the microphysical and optical properties of aerosols are based on a cluster analysis of the global AERONET database through 2010 [
17]. However, the aerosol characteristics vary locally, and they cannot be accurately described using a global aerosol model, which increases the uncertainty in AOD retrieval [
23].
In this paper, a monthly minimum land surface reflectance (MLSR) database for Xinjiang area (band: blue; resolution: 500 m; time span: 2010–2014) was established using MODIS surface reflectance product (MOD09A1), based on the minimum synthesis technique (MST) method to reduce the impact of cloud contamination [
20,
24,
25]. A bidirectional reflectance distribution function (BRDF) correction model for angle normalization was adopted to compensate the effect of surface reflectance anisotropy. In addition, the CE-318 sunphotometer data are used to determine the aerosol characteristics of the study area.
2. Field Measurements and Data Used
The CE-318 sunphotometer is an automatic ground-based radiometer measuring both direct solar irradiance and diffuse sky radiance for almucantar and principal solar planes with a 1.2° field of view. Although the channel wavelength configuration depends on the instrument version, filters at 440, 500, 675, 870, 940 and 1020 nm wavelengths are always present [
26], with an uncertainty of 0.01–0.02 [
27]. The Microtops II sunphotometer is a portable, manually-operated instrument which measures AOD through direct solar irradiance measurements in five wavebands (380, 500, 675, 870, and 1020 nm) with an uncertainty of ~0.015–0.02 [
28,
29].
In order to obtain the characteristic parameters of atmospheric aerosol in the study area, the CE-318 and Microtops II sunphotometers were used at two ground-observed sites (Dahuangshan and Wucaiwan site). The Dahuangshan site is located in a coal mining area with sparse vegetation coverage and the Microtops II sunphotometer observation was carried out from 10 July to 24 August 2014. The Wucaiwan site is located in a remote rural area near the Gurbantunggut Desert and the CE-318 sunphotometer observation was carried out from 15 July to 21 August 2014. A map of the study area is shown in
Figure 1, and the detailed information of sites and instruments is presented in
Table 1.
The MODIS (Terra) data products, namely calibrated radiance product (MOD02HKM), geolocation product (MOD03), surface reflectance product (MOD09A1), BRDF/albedo product (MCD43A1), and the Collection 6 (C6) aerosol product (MOD04), were obtained from the Level-1 and Atmosphere Archive & Distribution System (LAADS) Distributed Active Archive Center (DAAC) at the Goddard Space Flight Center (GSFC) (
http://ladsweb.nascom.nasa.gov). The information of these satellite products and their applications is listed in
Table 2.
3. Methodology
The top-of-atmosphere (TOA) reflectance (
), measured by a satellite sensor, is a function of surface and atmosphere optical parameters as well as solar/view zenith and azimuth angles; it can be estimated using Equation (1):
where
is the solar zenith angle,
is the view zenith angle,
is the relative azimuth angle,
is the aerosol optical depth,
is the aerosol reflectance resulting from multiple scattering in the absence of molecules;
is the Rayleigh reflectance resulting from multiple scattering in the absence of aerosols;
is the surface reflectance;
and
are the transmissions of the atmosphere on the sun-surface path and the surface-sensor path, respectively; and
is the atmospheric backscattering ratio to account for multiple reflections between the surface and atmosphere.
Equation (1) indicates that the TOA reflectance comes from two parts: the pure atmospheric contribution (the first two terms in Equation (1)) and the combination of the atmosphere and land (the last term in Equation (1)). A key step of the satellite-retrieved AOD is to remove the surface contributions from the satellite image. An overview of our retrieval algorithm is represented as a dataflow diagram in
Figure 2. This algorithm cannot retrieve cloud-contaminated pixels; therefore, before beginning the retrieval processing, we screened areas for the presence of clouds, and used the universal dynamic threshold cloud detection algorithm for cloud detection [
30]. The surface reflectance for a given pixel was then determined from the pre-calculated MLSR and angle information database, which was built using MOD09A1 products. The aerosol model is an important parameter that affects the precision of AOD retrieval and must be considered carefully. In this study, it was established by the ground-based measurements data.
3.1. Construction of Surface Database
In order to remove the surface contributions from the satellite signal, a database of surface reflectance was built for the AOD retrieval. The MOD09A1 dataset provides surface reflectance at 500 m resolution in sinusoidal projection. Each MOD09A1 pixel has the best possible observation during an 8-day period as selected by high observation coverage, small view angle, absence of clouds or cloud shadow, and low aerosol loading. Validation has been carried out for the MOD09 which indicated that 50.52% of the observations in the blue band were within the MODIS theoretical uncertainty of ±(0.005 + 5% ×
), where
is the surface reflectance [
31]. In this study, MOD09A1 datasets for 5 years (2010–2014) and six tiles (H23-H25 and V04-V05) were used to build the surface reflectance database to support AOD retrieval.
Ideally, the MOD09A1 product should be the surface reflectance corrected for aerosol effect and cloud masked in the quality flag. However, detection of thin cloud or subpixel cloud are difficult, and there is also uncertainly in the aerosol effect correction [
31]. Following the approach in the paper by Sun et al. [
32], we adopted the criterion of minimum land surface reflectance (MLSR) for database synthesis because it is the most efficient way to avoid cloud contamination. In each 500 m pixel, the lowest surface reflectance value in the twenty images time series, i.e., four images per month for 5 successive years, was identified as the clearest observation. In order to correct the effect of surface bidirectional reflectance in AOD retrieval, the angle information in MOD09A1 was stored in the MLSR database along with the minimum reflectance.
Figure 3 shows the distribution of surface reflectance in July, as well as the solar/view zenith and relative azimuth angles.
3.2. BRDF Correction Surface Reflectance
In the aerosol retrieval algorithm, the accurate determination of the surface reflectance is one of the most crucial concerns. For the DB algorithm, the surface reflectance is prescribed by one of several methods, dependent on location, season, and land cover type, from a global surface reflectance database in visible bands. However, all methods do not consider the influence of imaging geometry. As the surface has anisotropic reflection characteristics and the geometric angles of the satellite sensors vary, the surface reflectance of the MLSR is different from that of images used to retrieve AOD. Thus, the pre-calculated MLSR database cannot be directly used for AOD retrieval. In this study, we adopted the kernel-driven BRDF model for the correction of the effects of anisotropic reflection of the surface. The correction model was proposed based on the kernel-driven BRDF model. The kernel-driven BRDF model is a semi-empirical models which is derived as a simplification of physically based BRDF models, with the merit of its linear form and small number of model parameters. It can be generally described by Equation (2) [
33]:
where
and
are the volumetric and geometric kernels, respectively, which are functions of illumination and viewing geometry, describing volume and geometric scattering from surface elements
,
are the weights for the volumetric and geometric kernels, respectively, and
is to the weight of the isotropic reflectance.
By simply normalizing the kernel weights with the isotropic weight, we define the anisotropy shape factors as
and
. Then, the directional surface reflectance under geometric angle
,
,
can be predicted from the direction surface reflectance in the MLSR database under geometric angle
,
,
, with Equation (3):
Studies show that the combination of Ross–Thick and Li-SparseR (RTLSR) kernels works well with the observed data [
34], and the
,
, and
parameters of this combination are provided in the MCD43A1 products. However, as inversion of the BRDF is difficult, the BRDF parameters in MCD43A1 are noisy and cannot be directly used. Vermote et al. showed that the anisotropy of the surface is related to vegetation status [
35]. So, in this study, we classify the study area into three categories, i.e., sparse vegetation, median vegetation and dense vegetation, according to the threshold of the Enhanced Vegetation Index (EVI). Statistical mean values of A
1 and A
2 were calculated for each category from pixels marked as best quality (QC = 0) in the MCD43A1 product. The mean values of A
1 and A
2 were then applied for all pixels in the category. The statistical mean values of A
1 and A
2 are shown in
Table 3.
3.3. Aerosol Parameter Determination
To reduce the computation requirement, a look-up table (LUT) was constructed using the latest version (Version 2.1 Vector Code) of the 6S (second simulation of the satellite signal in the solar spectrum) (6SV) radiative transfer model (RTM) [
36]. The 6SV code is a widely-used radiative transfer code that simulates the satellite signal accounting for elevated targets. One of the practical advantages of 6SV is that it provides standard atmosphere and aerosol models. The 6SV model was used to construct a LUT of atmosphere optical parameters according to the parameters values listed in
Table 4.
In this work, the CE-318 ground-based data were used to determine the aerosol microphysical and optical parameters, such as the Ångström exponent (AE), complex refractive index (RI), single scattering albedo (SSA), and asymmetric factor (g). The AE can be used to determine aerosol models; for example, the AE of desert aerosols is −1.0 to 0.5 [
37]. For the desert model, the value of the real part of RI is 1.5–1.6 and the value of the imaginary part is approximately 0.01.
Table 5 shows the values of the AE and RI in the study. The aerosol optical parameters, including the SSA and g at four wavelengths (i.e., 440, 675, 870, and 1020 nm;
Figure 4), were retrieved from sky radiance almucantar measurements and direct sun measurements by a CE-318 sunphotometer. In this study, the first step is to calculate the average values of SSA and g at 440 and 675 nm, and then the inverse distance weighting interpolation [
38] is adopted to obtain the values of SSA and g at 550 nm, which is 0.827 and 0.910, respectively.
3.4. Rayleigh Correction for Elevation Effect
Within the MODIS blue wavelength range, where Rayleigh scattering is relatively important compared to longer wavelengths, neglect of polarization in the radiative transfer code leads to significant errors in the calculated reflectances [
39]. At sea level, the Rayleigh optical thickness (ROT) at visible channels can be estimated with an empirical function of wavelength
[
40]:
At a height above sea level, the ROT should be modulated by atmospheric pressure or elevation [
41]:
where
z is the ground height above sea level in kilometers, and 8.5 is the exponential scale height of the atmosphere. In this study, the MOD03 product was used for providing the altitude
z for each pixel.
3.5. Error Indicators
The error statistics of the algorithm were verified by comparing the satellite-retrieved AOD with ground-based measurements from sunphotometers. We used the following four indicators to evaluate the error: the correlation coefficient (
R), mean absolute error (
MAE), root mean square error (
RMSE), and
EE.
R is an indicator of relative agreement between satellite-retrieved AOD and ground-observed AOD; the MAE is the most natural measure of mean error magnitude; the
RMSE is used to measure the systematic and random differences between these two AOD observations; and the
EE, representing the confidence envelopes of the retrieval algorithm, is used to evaluate the quality of a new algorithm relative to MODIS C6 AOD. Good matches of satellite-retrieved AOD are reported when the satellite-retrieved AOD falls within the envelope. The statistical indicators are defined as follows:
where
is the ground-based sunphotometer measurement of AOD,
is the satellite-retrieved AOD;
and
are average values of
and
.
5. Conclusions
In this study, a modified aerosol retrieval algorithm was proposed for retrieving AOD over the arid/semiarid region of northern Xinjiang from MODIS data at 500 m spatial resolution. The assumptions in this algorithm are as follows: the variation in surface reflectance is small over a month, and single scattering albedo (SSA) and asymmetry factors (g) are regionally constant for a few days. The MODIS surface reflectance products (MOD09A1) were used to determine the surface reflectance, and a LUT was constructed based on the 6SV RTM, which uses SSA and g values from ground-based measurements. The retrieved AOD values were validated by ground-based sunphotometer observations in two sites and compared with the MODIS DB AOD products. The results show that the new algorithm accurately retrieved aerosol AOD over the arid/semiarid region of northern Xinjiang, and the retrieved aerosol distribution contained more spatial details and variability than that of the DB AOD products. In this study, a BRDF correction model was applied to reduce the effect of surface reflectance anisotropy in AOD retrieval; thus, the retrieval ability is improved especially over the area with obvious surface bidirectional reflectance characteristics. However, according to the uncertainty analysis, the proposed algorithm has some limitations that should be improved: (1) Errors may occur by using constant values of SSA and g for the day of retrieval; (2) the MLSR database approach was sometimes found to be unsuccessful over snow surfaces, particularly when seasonal changes are significant, such as the snow melting in March and accumulation in November, respectively; (3) the BRDF correction only considered the mean values of anisotropy shape factors; thus, it ignored the variation of BRDF shape within the categories of land surface. These limitations will be explored in our future studies.