Estimating Particulate Matter Concentration over Arid Region Using Satellite Remote Sensing: A Case Study in Makkah, Saudi Arabia

The air quality indicator approximated by satellite measurements is known as an atmospheric particulate loading, which is evaluated in terms of the columnar optical thickness of aerosol scattering. The effect brought by particulate pollution has gained interest among researchers to study aerosol and particulate matter. In this study we presents the potentiality of retrieving concentrations of particulate matter with diameters less than ten micrometer (PM10) in the atmosphere using the Landsat 7 ETM+ slc-off satellite images over Makkah, Mina and Arafah. A multispectral algorithm is developed by assuming that surface condition of study area are lambertian and homogeneous. In situ PM10 measurements were collected using DustTrak aerosol monitor 8520 and their locations were determined by a handheld global positioning system (GPS). The multispectral algorithm model shows that PM10 high during Hajj season compared to other season. The retrieval dataset gives the accuracy > 0.8 of R coefficient value over Makkah, Mina and Arafah. These results provide confidence that the multispectral algorithm PM10 models can make accurate predictions of the concentrations of PM10.


Introduction
Air pollution is currently one of the major problems in developed countries as well as in developing countries.The five pollutants which are ozone (O 3 ), nitrogen oxides (NOx), carbon monoxide (CO), sulphur dioxide (SO 2 ) and particulate matter (PM) are referred to as criteria of Air Pollution Index (API) by the Department of Environment (DOE) Malaysia.Department of environment, Malaysia (DOE) stated that most of API is based on PM10 which is higher than other pollutants.Particulate matter is a general term used for aerosols, small liquid droplets, or solid particles that are found in our air typically in the size range of 0.01 to 100 micrometers.These are much larger than individual molecules.The United States Environmental Protection Agency (EPA) uses the abbreviations PM10 and PM2.5 to specify certain sizes of particulates.PM2.5 stands for particulate matter less than 2.5 microns (also called micrometers, μm).PM10 refers to particles greater than 2.5 microns up to 10 microns.In some areas, particulate matter (PM) can be very heavy because of high levels of industrial activity or natural environmental conditions from a variety of sources, such as vehicles, factories, construction sites, farming, unpaved roads, burning wood, and blowing sand and dust in desert environments.In these types of environments, PM 10 particles are small enough to be inhaled and accumulate in the respiratory system especially in regard to cardio-vascular illnesses and reduce visibility by their scattering and absorption of radiation (Husar et al., 1981;Ball & Robinson, 1982).
Air quality monitoring at urban and regional scales has traditionally been done using a network of ground monitoring stations combined with dispersion models that predict air quality between monitor locations.For this regard, the availability of satellite remote sensing can provide a synoptic picture of air quality in a regional air shed, including information about sources and source locations for isolated events.The radiative properties of each component such as the path radiances/atmospheric reflectance at different wavelengths and satellite viewing geometries are calculated and archived in a look-up table.During the retrieval process, the path radiances/atmospheric reflectance of each aerosol mixture is calculated, and compared with the in situ data concentration.These sensors also can potentially be used to monitor air quality in rural or remote regions with no ground-based monitoring network.
In this study, a multispectral algorithm model of PM10 was developed to estimate the distribution concentration of particulate matter less than 10 micron (PM10) using Landsat 7 ETM+ slc-off data calibrated with in situ concentration measurements.The estimated model was then applied using PCI Geomatica image processing software to generate air pollution PM10 maps for particular date.The accuracy analyses were carried out by comparing the estimated data generated using multispectral algorithm model with in situ data.The results show that the retrieval dataset gives the accuracy > 0.8 of R coefficient value over Makkah, Mina and Arafah.

Study area
Figure 1 shows the selected study area of Makkah, Mina and Arafah over Saudi Arabia.The Holy City of Makkah was an arid-urban area (Latitude 21°25'19" North Meridian 39°49'46") is at an elevation of 277 m above sea level, and approximately 80 km inland from the Red Sea.The elevations of Makkah AI Mukarramah are a group of mountains and black rocky masses which are granitic basement rocks (Al-Jeelani, 2009).Mountains are traversed by a group of valleys, such as the Ibrahim valley.The Kaabah's location is in this valley.
Makkah climate is different from other Saudi Arabian cities, retains its warm temperature in winter (November to March), which can range from 17 °C at midnight to 25 °C in the afternoon.During summer (April to October), temperatures are considered very hot and break the 40 °C mark in the afternoon dropping to 30 °C in the evening.Rain is very rare with an average of 10-33 mm usually falls in December and January; and the humidity ratio is about 45-53 %.Winds are north-eastern most of the year time.Some unusual events often happen during the year, such as dust storms in summer, coming from the Arabian Peninsula's deserts or from North Africa (Al-Jeelani, 2009).

Methodology
The methodology process generally were divided into four major parts: data acquisition, pre-processing, data processing and finally, accuracy and validation of results.All data preprocessing and processing steps were carried out using PCI Geomatica 10.2 software.

Satellite image
The acquisition dates of the Landsat ETM+ Scenes employed in the air quality change detection process within seasonal variation (Hajj season and non-Hajj season).All Landsat 7 ETM+ scenes were downloaded from the United States Geological Survey (USGS) as a Level 8 product.These imageries were acquired in NLAPS format 30x30 m pixels.The ETM+ scenes, WRS Path/Row 169/45 were captured by the Landsat 7 ETM+ satellite on 29 th December 2006, and 19 th January 2009 respectively.All Landsat 7 ETM+ scenes were selected based on the minimum percentage of cloud cover (<10%) and the availability of ground truth data prior to acquisition.
To reduce effect due to zenith angle and surface reflectivity effect, both imageries were selected among the same value for sun zenith angle, azimuth angle.Image acquisitions after March 2003 are corresponding to SLC-off images which have gap line anomaly caused by scan problem.The malfunction of the SLC mirror assembly resulted in the loss of approximately 22 % of the normal scene area (Storey et al., 2005).Note that the SLC failure has no impact on the radiometric performance with the valid pixels.All scenes were affected by the failure in SLC (SLC-off mode) that occurred after 2003, so parts of the data in the scenes are lost.The acquisitions of Landsat 7 ETM+ scenes analyzed in this study are listed in Table 1.All meteorological data (Table 2) used in this study were taken from Weather Underground webpage.Cogliani (2001) and Rodr´ıguez et al. (2009) also used Weather Underground in their research.

Ground truth data
The ground truth data were obtained from the field survey covers all type of ground surface and scattered point throughout Makkah, Mina and Arafah using using DustTrak aerosol monitor 8520.The DustTrak aerosol monitor 8520 utilizes a specially shaped inlet to inertially separate particulate matter less than 10 μm fraction (which is collected on glass-fiber or quartz fiber filter media) and greater than 10 μm fraction (which is discarded).Heal et al. (2000) concluded that the DustTrak aerosol monitor 8520 (TSI Inc.) demonstrated excellent functionality in terms of ease of portability and real time data acquisition.The GPS measurements were taken for determining location of PM10 and AOT data collected which allows for users to compile ground measurements air quality data sets directly from the field as part of 'ground truthing' (Cunningham, 1998).The ground truth data of PM10 were divided into two groups; half of the numbers were randomly selected for calibration of algorithm and the other half for accuracy analysis.

Geometric and distortion correction
The Landsat 7 ETM+ satellite image was rectified using the second order polynomial coordinates transformation to relate groud control points in the map to their equivalent row and column positions in the Landsat 7 ETM+ scences.Corrected images were projected to Universal Transverse Mercator Projection with UTM 37 Q D000 WGS 1984 Datum.The reference points used to resample the satellite images were taken from 14 ground control point (GCP) collected at study area.

Radiometric and Atmospheric Correction
Radiometric correction is applied by transforming the values of DN to radiance or reflectance values There are different levels of radiometric calibration.The first converts the sensor DN to at-sensor radiances and requires sensor calibration information (Mather, 2004).The second is the transformation of the at-sensor radiances to radiances at the earth's surface.Radiometric correction is applied by transforming the values of DN to radiance or reflectance values through the algorithm as follows given by Chander et al. (2009): The value of Grescale (c1) and Brescale (c0) for Landsat 7 ETM+ used in this study can be found in Appendix A. Also can be expressed as: Where: For relatively clear Landsat scenes, a reduction in between-scene variability can be achieved through a normalization for solar irradiance by converting spectral radiance, as calculated above, to planetary reflectance or albedo.This combined surface and atmospheric reflectance of the Earth also known as top of atmosphere reflectance (TOA) is computed with the following formula (Mather, 2004): Where: Atmospheric correction was carried out using ATCOR2 available with PCI Geomatica using algorithms developed by Dr. Rudolf Richter (Richter 1996a(Richter , 1996b(Richter , 1997(Richter , 2005;;Ritcher et al., 2009).It calculates correction for flat areas applying constant or varying atmosphere accounting for adjacency effect.Atmospheric corrections widely used in hyper spectral imagery to derive surface reflectance without atmospheric effects.
Automatic calculate haze and cloud' would be the first run of ATCOR.The output files containing the haze and cloud mask.This mask can be edited if haze is not correctly assigned, e.g.defining additional haze areas or deleting wrongly assigned haze areas (Wen & Yang, 2008).Then ATCOR could be run again with 'Load Haze and cloud from file' employing this edited mask, which might yield better results for the haze removal.
The atmospheric correction algorithm calculates the surface reflectance using the default scale factor 4, i.e. the percent reflectance range 0 -100% is multiplied with the factor 4 in the output file.So an output value of DN=200 corresponds to a surface reflectance of 200/4=50% (or 0.5 for 0-1 reflectance range) and the output is coded as 8bit/pixel.Therefore, the maximum output value is 255, representing a surface reflectance of 255/4=63.75%(Geomatica Focus User Guide Version 10.0, 2005).Larger values will be truncated at 255.
ATCOR2 is based on a database of atmospheric correction functions stored in look-up tables.The database consists of a broad range of elevation information setup, sensor information, atmospheric information and correction parameter as in Tables 3. The result of ATCOR2 is a ground or surface reflectance image in each spectral band with a relative error of approximately 10 % (Lehner et al., 2004).

Data Processing
After undergo radiometric correction, the reflectance measured from the satellite (reflectance at the top of atmospheric, TOA) was subtracted by the amount given by the surface reflectance to obtain the atmospheric reflectance.The atmospheric reflectance was then related to the PM10 using the regression algorithm analysis.PCI Geomatica EASI modeling was used to input the developed multispectral algorithm.PM10 maps were generated using proposed algorithm based on the highest R and lowest RMSE values.The final results were in color coded image of PM10.

PM10 multispectral algorithm model
The Mie scattering theory was applied to compute the aerosol phase function and spectral optical depth, based on size distribution, real and imaginary index (King et al., 1999;Fukushima et al., 2000). (4) In the single scattering approximation (Popp et al., 2004), the path radiance is proportional to the aerosol optical thickness, a  , the aerosol scattering phase function, (5) The method can be applied to satellite imagery for which it is a priori known that vegetation are present, accounting on its geographic location and occurrence season in which the image was taken.The Equation ( 5) is rearranged to become Equation ( 6) (Xia, 2006).By neglecting molecule scattering due to Rayleigh (Paronis & Hatzopoulos, 1997;Kaufman & Tanre, 1998), Equation ( 6) becomes ( 7).The Mie theory, therefore, may be used for describing most spherical particle scattering systems, including Rayleigh scattering (Hahn, 2009).
(7) a a s v a s v P ( , , ) ( , , ) 4 So the algorithm of AOT for single band or wavelength (λ) is simplified as: Equation ( 8) is rewrite into three band equation as Equation ( 9).
Where i R  is the atmospheric reflectance (i = 1, 2 and 3 corresponding to wavelength for satellite), and j a is the algorithm coefficient (j = 0, 1 and 2) are empirically determined.
The relation between PM and AOT is derived for a single homogeneous atmospheric layer containing spherical aerosol particles.The mass concentration at the surface is obtained after drying the sampled air is given by Koelemeijer et al. (2006). (10) Hence, it can be expected that the parameter PM correlates better with AOT directly.Using the information, obtained by the spectral AOT retrieval, a method has been developed to retrieve particulate matter concentrations.
Several studies showed that PM10 and AOT have linear relationship correlation (Glantz et al., 2007).Chu et al. (2003) and Sifakis et al. (2002) also found that correlation as high as 0.78 to 0.95 are retrieved between the AOT values and PM10 measurements.As in literature PM 10 and AOT are better RMSE in linear correlation than exponential.Fraser et al. (1984), Kaufman et al. (1990), Gasso andHegg (1998, 2003) also attempts to relate PM and AOT measurements by linear equation.By substitute AOT in term of PM10, Equation ( 9) into (11), and the algorithm for single band or wavelength (λ) of PM10 is simplified as.
Where i R  is the atmospheric reflectance (i = 1, 2 and 3 corresponding to wavelength for satellite), and j a is the algorithm coefficient ( j = 0, 1 and 2) are empirically determined.

Accuracy and Validation of Results
The accuracy assessment and validation of results obtained was performing with ground truth data.Accuracy and validation analysis of results were perform using the new algorithm with PM10 ground truth values retrieved using DustTrak aerosol monitor 8520.

Results and Discussions
Distributions of PM10 with respect to atmospheric reflectance for the red band (b3), green band (b2) and blue band (b1) show in Figure 2 with respect to 2 sets of data 29 th December 2006 and 19 th January 2009 respectively.Table 4 shows the comparison values of R and RMSE for various type of algorithm using regression analysis.Then, PM10 maps were generated based on the highest R values and the lowest RMSE values, where the highest value of correlation coefficient, R is 0.888 and the linear regression model as equation ( 12).
Where PM10 is equal to PM10 concentration in (µg/m³), R λ1 , R λ2 and R λ3 are equal to the atmospheric reflectance/path radiance in blue, green and red band from Landsat 7 ETM+ satellite images respectively.The increasing of atmospheric reflectance/path radiance due to particulate matter are linear with concentration of PM10 in the study area.It is believed that the relatively high RMSE was due to limited number of air pollution stations used.Future study will consider of using more air pollution stations as well as other value-added ancillary data in order gain better and reliable accuracy.
The colour-coded PM10 maps generated from Landsat 7 ETM+ satellite images using developed multispectral algorithm were shown in Figures 3 and 4 for 29 th December 2006 and 19 th January 2009 respectively.PM10 concentrations are indicated through the color of red for high PM10 value and blue for low PM10 value.4) is quite high due to construction activities at Jamrah at Mina.Physical geography and topography of the area corresponding to the surrounding area of rocky mountain known as valley area where enables a weak air flow contribute significantly to the high concentration of PM10.Mahmood and Abdul Aziz (2008) also found that PM10 daily concentrations in the atmosphere of Mina valley ranged between 191 -262 μg/m 3 during the presence of pilgrims in Mina compared to the European standard of 50μg/m 3 during Hajj Seasons of 1424 and 1425 H (2004 -2005).This variation could be attributed the metrological conditions, despite the little variation between the two successive years, and Hajj circumstances.
The validation of the results generated using multispectral algorithm model with insitu data of PM10 shows that combination of two days dataset gives the accuracy > 0.8 of R coefficient value as shows in Figures 5 and 6.A good correlation agreement between measured PM10 and estimated PM10 obtained in this study indicates that the accuracy of the proposed algorithm model is high and can estimate the concentration of PM10 with reasonable value.

Conclusions
This study indicates that the air pollution PM10 can be mapped using Landsat 7 ETM+ slc-off imagery using the developed multispectral algorithm of PM10 model.The RGB algorithm PM10 model was used to generate PM10 concentration which gives highest R and lowest RMSE value shown from regression analysis.It add evident the linear relationship between PM10 and path radiances/atmospheric reflectance Landsat 7 ETM+ red band (R λ3 ), green band (R λ2 ) and blue band (R λ1 ).Air pollution is a growing problem in recent years due to the development of megacities like Makkah, which caused environmental consequences of air pollution of the construction activities and the vehicle smoke.The validation of the results generated using multispectral algorithm of PM10 shows that all dataset gives the accuracy > 0.8 of R coefficient value as shows in Figures 5  and 6.It shows that the developed multispectral algorithm is worked excellently in determining PM10 value within arid region of Makkah.
Radiance at the sensor's aperture in W/m 2 /sr/µm Grescale = Rescaled gain (the data product "gain" contained in the Level 1 product header or ancillary data record) in W/m 2 /sr/µm/DN Brescale = Rescaled bias (the data product "offset" contained in the Level 1 product header or ancillary data record ) in W/m 2 /sr/µm QCAL = the quantized calibrated pixel value in DN LMIN λ = the spectral radiance that is scaled to QCALMIN in W/m 2 /sr/µm LMAX λ = the spectral radiance that is scaled to QCALMAX in W/m 2 /sr/µm QCALMIN = the minimum quantized calibrated pixel value (corresponding to LMIN λ ) in DN

Figure 2 .R
Figure 2. Graph of PM10 data versus atmospheric reflectance for the three bands, red band (R λ3 ), green band (R λ2 ) and blue band (R λ1 ) for 2 days respectively 29 th December 2006 (S1) and 19 th January 2009 (S2) Table 5 shows that almost 30% of the pixels in the images are represent black colour indicating of unfilled gap and cloud .This is due to the presence a lot of cloud and also gap caused by SLC-off image.PM10 pollutants are decreased after-Hajj period of year 2009 compared to corresponding results of year 2006.High values of PM10 are observed during Hajj season on 29 th December 2006 where the maximum value is 151 µg/m³ and minimum value is 64 µg/m³.In contrast, low values of PM10 on 19th January 2009 are due to consequence of non Hajj season which reducing activities around study area where the maximum value is 139 µg/m³ and minimum value is 40 µg/m³.However, the PM10 distribution over Makkah and Mina on 19 th January 2009 (as shown in Figures area