SPATIAL INTERPOLATION OF AEROSOL OPTICAL DEPTH POLLUTION : COMPARISON OF METHODS FOR THE DEVELOPMENT OF AEROSOL DISTRIBUTION

Air pollution is a growing problem arising from domestic heating, high density of vehicle traffic, electricity production, and expanding commercial and industrial activities, all increasing in parallel with urban population. Monitoring and forecasting of air quality parameters are important due to health impact. One widely available metric of aerosol abundance is the aerosol optical depth (AOD). The AOD is the integrated light extinction coefficient over a vertical atmospheric column of unit cross section, which represents the extent to which the aerosols in that vertical profile prevent the transmission of light by absorption or scattering. Seasonal aerosol optical depth (AOD) values at 550 nm derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor onboard NASA’s Terra satellites, for the 10 years period of 2000 2010 were used to test 7 different spatial interpolation methods in the present study. The accuracy of estimations was assessed through visual analysis as well as independent validation based on basic statistics, such as root mean square error (RMSE) and correlation coefficient. Based on the RMSE and R values of predictions made using measured values from 2000 to 2010, Radial Basis Functions (RBFs) yielded the best results for spring, summer and winter and ordinary kriging yielded the best results for fall.


Introduction
Aerosols are generally consists of minute solid and liquid particles held in suspense in the atmosphere, with aerodynamic diameter, which ranges from 10 −3 μm to 10 μm.These particles have significant effects on the aerial visibility, the earth's radiation budget as well as cloud formation, precipitation and human wellbeing (Solomon, 2007, Kaufman et al., 2002, Pope III et al., 2002).One of the most important parameters in the study of aerosols is the aerosol optical depth (AOD).AOD is the extinguishing phenomenon of optical beam power as a result of the presence of aerosols.Thus, AOD is used to determine the aerosol column concentration.For fine-mode aerosols, the AOD generally decreases as the wavelength increases.
The MODIS, which is mounted on the Terra satellites measures radiances across 36 visible, near infrared, and infrared channels from 415 to 14 235 nm (King et al., 1999) utilizing 0.25, 0.50, or 1 km spatial resolutions.This equipment has a viewing swath of about 2300 km in width, fragmented into 5 min "granules" per ~2030 km of length.An in depth description of the over ocean algorithm implemented on the MODIS for retrieving AOD, can be found in Remer et al. (2005); that of the dark-target algorithm can be found in Kaufman et al. (1997).However, as reported by Levy et al. (2010), the new order of overland dark-target algorithm was first proposed by Levy et al. (Levy et al., 2007, Levy et al., 2010), which addresses the limitations peculiar to the older versions of MODIS (Remer et al., 2005, Levy et al., 2007).In retrieving aerosol products using MODIS dark target land AOD, relationship between surface reflectance at VIS and SWIR is of specific importance, which generally holds for the majority of vegetative areas, but generally fails to hold for arid, semiarid, urban and desert regions (Hsu et al., 2004).On the other hand, the MODIS Deep Blue algorithm is developed for retrievals involving brightly reflecting surfaces.Such retrievals involve the use of blue bands where there exist low enough surface reflectances enabling the retrievals (Hsu et al., 2006).
In regions where there exist no previous measurements, spatial interpolation is a suitable technique for estimating the concentration of a pollutant (Isaaks and Srivastava, 1989, Denby et al., 2005, United, 2004) by employing a point-by-point forecasting utilizing the stochastic approach (Pisoni et al., 2009).The deterministic and stochastic spatial approaches have also been applied in some studies.Some of these deterministic approaches include the nearestneighbor and polynomial interpolation techniques (United, 2004, Isaaks andSrivastava, 1989); while geostatistical models such as kriging (Beelen et al., 2009, Jourdan, 2009) and cokriging (Cressie, 1992, Isaaks andSrivastava, 1989) constitute stochastic approaches.In particular, researchers have used kriging to map prior air pollution data involving such compounds as NO2, PM10, O3, SO2, CO (Beelen et al., 2009).For example, the study by Mulholland et al. (1998) utilized universal kriging to interpolate 1-h and 8-h data of 10 stations in the USA.Similarly, the study by Rojas-Avellaneda (2007) carried out a comparative investigation of the peak period statistics of 16 stations in Mexico City using inverse distance weighting; while Sanchez et al. (Fernández de Castro et al., 2003) employed kriging to interpolate data from 8 locations in Guadalajara urban area in Mexico.In addition, Son et al. (Son et al., 2010) used interpolation methods on data generated from 8h ozone concentration from 13 locations in Ulsan, Korea; while Montero et al. (Montero et al., 2010) utilized ordinary kriging to attain ozone data from over 27 locations within Madrid.

MODIS products
The spatially and temporally collocated MODIS data pairs spanning the years 2000_2010 for the full record of MODIS/Terra are acquired through the Multi-Sensor Aerosol Products Sampling System (MAPSS, http://giovanni.gsfc.nasa.gov/mapss/)(Ichoku et al., 2002, Petrenko et al., 2012).MODIS Level 2 Collection 5.1 MOD04 aerosol data from 4 July 2002 through 10 January 2011 are used.Note that: (a) vegetated surfaces are not 'dark' at the 550-nm wavelength and, therefore, the AOD at this wavelength over land is derived from the retrieved AODs at the 470-nm and 660-nm channels (Levy et al., 2010) and (b) the MODIS Dark Ocean product provides two AOD datasets, one from the inversion using the bestfitting aerosol model and another from the average of inversions using several well-fitting models (ATBD-2006; found online at http://modis-atmos.gsfc.nasa.gov/MOD04_L2/index.html); the latter is used for this research.The quality of each MODIS AOD retrieval is represented by its associated quality flag ranging from 3 (high confidence) to 0 (low or no confidence) (Levy et al., 2010).The Land_And_Ocean AOD dataset is generated from a union of AODs retrieved respectively by the Dark Land and Dark Ocean algorithms.It is noted, however, that Collection 5.1 has two different variable names for Land_ And_Ocean AOD; one is the 'Image_Optical_Depth_Land_ And_Ocean' that has no QA involved in its production and another is 'Optical_Depth_Land_And_Ocean' that requires quality flags_0 over land and]0 over ocean (ATBD, 2006); the latter data variable is consequently used here.However, unlike the individual Land and Ocean AOD datasets, the combination product does not report QA flags.

Spatial Interpolation Methods
In this subsection, the interpolation methods employed in the current study will be briefly illustrated.The study by (Cressie, 1992, Goovaerts, 1997, Chiles and Delfiner, 2009, Webster and Oliver, 2007) provide detailed insight into geostatistical theories.In Spatial interpolation, regionalized values are estimated at un-sampled points using a weight of observed regionalized values.Mathematically, spatial interpolation is expressed as: (1) Where Zg denotes the interpolated value at point g, Zsi denotes the observed value at point i, ns denotes the total number of observed points, while λ = λi represents the weight contributing to the interpolation.Here, computing the weights λ, which is required in the interpolation is somewhat non trivial.The various means of computing the weights will be illustrated.

Comparison of interpolation methods
The common methods of comparison of the interpolation methods is the cross validation as well as validation with an independent data set.However, given the limited sample size of the current study, the cross validation was applied.Cross validation comprises of the consecutive removal of a data point and interpolating the value from the remaining observations.The predicted value is then compared with the measured value (Mueller et al., 2004).The accuracy of the results was tested by using the root mean square error (RMSE): (2) Where z(xi) denotes the observed value at location I, z*(xi) denotes the interpolated value at location i, while n denotes the sample size.The smaller the RMSE, the fewer the errors.
Since the cross validation is only effective in validating the prediction accuracy at a sample site but cannot depict the spatial difference of the interpolation techniques; the raster analysis function of ESRI ArcGIS was used in the current study in order to afford the comparison of the area and spatial differences of contaminated areas estimated by various interpolation techniques.

Data Analysis
Statistical results indicated that the aerosol optical depth was non normally distributed (Figure 2).Data sets were analyzed with different software packages.Maps were produced with GIS software ArcGIS and its extension of Spatial Analyst.The geostatistical analysis and the probability calculations were carried out with GS+ and geostatistic extension of ArcMap.Data characteristic of seasonal mean AODs can be described by statistics, such as mean, median, skewness, kurtosis and so on.Table 1 shows the data characteristic and that after data transformation.Obviously, it is more approximate Gaussian distribution after transformation.Mean is closed to median in log transformation, and skewness is closer to 0, therefore, log transformation is adopted.

Spring
The Figures 3 and Table 2 illustrate the interpolation surfaces analyzed using Seven different techniques.These techniques include simple kriging, Disjunctive kriging, Ordinary kriging, IDW, RBF, GPI and LPI.These techniques are mapped using suitable parameters.For example,  = 2 was ideally chosen for the IDW technique.Each of the seven schemes also underwent cross-validation, and upon a number of computations, R and RMSE were used for validating the interpolation methods.Table 2 illustrates that performance is optimal when the RBF has maximum R and minimum RMSE during spring.On the other hand, the polynomial interpolation does not aptly interpolate data since it is unable to pass through the majority of the specified points.Also, the level of accuracy registered by the LP is more superior to That of the GP, due to the fact that the latter takes the global trends into consideration.Nevertheless, it provides a much smoother surface compared to the former.

Summer
Table 2 illustrate the interpretation of the interpolation surfaces using the seven techniques mentioned.Here, these techniques were mapped using optimal parameters.For example,  = 2 was ideally chosen for the IDW technique.Each of the seven schemes also underwent cross-validation, and upon a number of computations, R and RMSE were used for validating the interpolation methods.Table 2 illustrates that performance is optimal when the RBF has maximum R and minimum RMSE during the summer.On the other hand, the polynomial interpolation does not aptly interpolate data since it is unable to pass through the majority of the specified points.Also, the level of accuracy registered by the LP is more superior to that of the GP, due to the fact that the latter takes the global trends into consideration.Nevertheless, it provides a much smoother surface compared to the former.

Fall
Table 2 illustrate the interpretation of the interpolation surfaces using the seven techniques mentioned.Here, these techniques were mapped using optimal parameters.For example,  = 2 was ideally chosen for the IDW technique.Each of the seven schemes also underwent cross-validation, and upon a number of computations, R and RMSE were used for validating the interpolation methods.Table 2 illustrates that performance is optimal when the Ordinary kriging has maximum R and minimum RMSE during fall.On the other hand, the polynomial interpolation does not aptly interpolate data since it is unable to pass through the majority of the specified points.Also, the level of accuracy registered by the LPI is more superior to that of the GPI, due to the fact that the latter takes the global trends into consideration.Nevertheless, it provides a much smoother surface compared to the former.

Winter
Table 2 illustrate the interpretation of the interpolation surfaces using the seven techniques mentioned.Here, these techniques were mapped using optimal parameters.For example,  = 2 was ideally chosen for the IDW technique.Each of the seven schemes also underwent cross-validation, and upon a number of computations, R and RMSE were used for validating the interpolation methods.Table 2 illustrates that performance is optimal when the RBF has maximum R and minimum RMSE during the winter.On the other hand, the polynomial interpolation does not aptly interpolate data since it is unable to pass through the majority of the specified points.Also, the level of accuracy registered by the LPI is more superior to that of the GPI, due to the fact that the latter takes the global trends into consideration.Nevertheless, it provides a much smoother surface compared to the former.

Conclusion
The seven interpolation methods of seasonal mean AOD of 10 years over the world, was compared.The accuracy of interpolation was determined by cross-validation.RMSE and R were chosen as validation criteria because of high sensitivity in this study.Visual results were given in Table 2. From the comparison of cross-validation criteria, it has been observed that the best interpolation scheme is RBF for spring, summer and winter seasons and ordinary kriging for fall season.

Fig. 1
Fig. 1 Map of AERONET stations used in this study.a) Map of the 236 AERONET stations during the spring season.b) Map of the 239 AERONET stations during the summer season.c) Map of the 239 AERONET stations during the fall season.d) Map of the 202 AERONET stations during the winter season.

Figure 2 .
Figure 2. Histograms of raw values of mean AOD

Table 1 :
Data characteristic of seasonal mean AOD data