Improving retrievals of regional fine particulate matter concentrations from Moderate Resolution Imaging Spectroradiometer (MODIS) and Ozone Monitoring Instrument (OMI) multisatellite observations

Abstract

A combination of multiplatform satellite observations and statistical data analysis are used to improve the correlation between estimates of PM_2.5 (particulate mass with aerodynamic diameter less that 2.5 µm) retrieved from satellite observations and ground-level measured PM_2.5. Accurate measurements of PM_2.5 can be used to assess the impact of air pollution levels on human health and the environment and to validate air pollution models. The area under study is California's San Joaquin Valley (SJV) that has a history of poor particulate air quality. Attempts to use simple linear regressions to estimate PM_2.5 from satellite-derived aerosol optical depth (AOD) have not yielded good results. The period of study for this project was from October 2004 to July 2008 for six sites in the SJV. A simple linear regression between surface-measured PM_2.5 and satellite-observed AOD (from MODIS [Moderate Resolution Imaging Spectroradiometer]) yields a correlation coefficient of about 0.17 in this region. The correlation coefficient between the measured PM_2.5 and that retrieved combining satellite observations in a generalized additive model (GAM) resulted in an improved correlation coefficient of 0.77. The model used combinations of MODIS AOD, OMI (Ozone Monitoring Instrument) AOD, NO₂ concentration, and a seasonal variable as parameters. Particularly noteworthy is the fact that the PM_2.5 retrieved using the GAM captures many of the PM_2.5 exceedances that were not seen in the simple linear regression model.

Implications:

Particulate Mass (PM) in the air is a concern because of its effect on climate and human health. PM concentrations retrieved from satellite observations of aerosol optical depth can provide broad regional coverage that is not attained by surface sites. The techniques developed in this paper have resulted in greatly improved correlations between PM retrieved from satellite observations and PM from surface measurements in areas where the correlation is typically low. These improved retrievals can be used to fill in the gaps between surface sites and validate air quality models that are used for air quality forecasts and epidemiological studies.

Supplemental Materials: Supplemental materials are available for this paper. Go to the publisher's online edition of the Journal of the Air & Waste Management Association for information on the effects of grid size and MODIS data quality flags on the GAM results. The coefficients for the GAMs used in this study are also listed.

Introduction

Particulate mass (PM) in the air is a concern in many regions of the United States and in the world because of its effect on climate and human health (Solomon et al., 2007; Samet et al., 2000; Pope et al., 2009; Krewski et al., 2009). The Intergovernmental Panel on Climate Change has identified PM as one area in which more research is needed to reduce uncertainty in estimates of climate change (Solomon et al., 2007). Also, epidemiology studies have shown that life expectancy is positively correlated with a decrease in the long-term exposure to fine-PM concentration (Pope et al., 2009). Elevated concentrations of PM_2.5 (PM with aerodynamic diameter less that 2.5 μm) in California are associated with about 18,000 premature deaths each year. (Tran et al., 2009). Environmental agencies around the world monitor PM concentrations at surface sites. In the United States, there are approximately 1600 U.S. Environmental Protection Agency (EPA)-approved PM_2.5 monitoring sites; however, these measurements are made at a single point. Point measurements are not necessarily representative of regional variability and satellite observations can help “fill in the gaps” between these measurement stations. In recent years, application of satellite observation to surface air quality has advanced considerably. Researchers have been working for years to use satellite observations to improve air quality forecasts and to provide information to assist regulators and policy makers (Wang and Christopher, 2003; Chu et al., 2003; Al-Saadi et al., 2005).

Satellite observations provide a total-column aerosol optical depth (AOD) product, which is a measure of the extinction by aerosol total column and it can be related to PM (Koelemeijer et al., 2006). Correlations between AOD with PM using single-variate, least-squares linear regression (a simple linear regression) show considerable variation, with correlation coefficients (r ²) between 0 and 0.85 depending on the season and location, or region, of comparison. Hoff and Christopher (2009) review and discuss over 30 papers that have explored the relationship between surface PM and satellite AOD measurements. In some cases, notably in the western United States, there is little or no correlation (Al-Saadi et al., 2005; Engle-Cox et al., 2004; Kahn et al., 2005; Gupta et al., 2006). In general, satellite retrievals of AOD are most accurate over areas of dark vegetation that is more typical in the eastern United States (Levy et. al. 2010).

Many schemes have been applied to improving satellite-surface correlations. Satellite AOD measurements from MODIS (Moderate Resolution Imaging Spectroradiometer), MISR (Multiangle Imaging SpectroRadiometer), POLDER (POLarization and Directionality of Earth's Reflectance), and GOES (Geostationary Operational Environmental Satellite) have been the subject of these comparisons. Gupta et al. (2006) investigated linear correlations between MODIS AOD and PM_2.5 at 25 urban areas. They found that correlations varied between 0.11 and 0.85. Kacenelenbogen et al. (2006) used simple linear regressions to compare POLDER AOD and surface PM_2.5 measurements over France and obtained a correlation coefficient of 0.55 for the general area. Hoff and Christopher (2009) suggest that the poorer correlations in the western United States were due to a wider variation in aerosol types (more nitrate than sulfate and more biomass burning smoke than over the eastern part of the United States), higher surface reflectivity making passive satellite AOD retrieval difficult, and more elevated plumes in the AOD signatures. Elevated aerosol plumes, typically associated with long-range transport, are integrated into the column AOD retrieval but are not reflected in the surface PM measurements. Satellite AOD measurements have been combined with other parameters, such as the vertical structure of aerosol as calculated from a chemical transport model (van Donkelaar et al., 2010; Liu et al., 2009b), or estimates of boundary layer height and relative humidity from models (Liu et al., 2005; Kacenelenbogan et al., 2006). Liu et al. (2009a) used fractional component AOD from MISR and aerosol vertical distributions from the GEOS-Chem transport model to constrain estimates of surface PM_2.5. (GEOS—Chem is a global 3-D chemical transport model. See http://acmg.seas.harvard.edu/geos/index.html.) Liu et al. (2009a) used generalized linear regression models (GLMs) and model-generated meteorological fields to compare the ability of MISR and MODIS AOD to predict surface PM_2.5 in the St. Louis region during 2003. The correlation coefficient between retrieved and measured PM using a GLM was r ² = 0.62 and r ² = 0.51 for MISR and MODIS, respectively.

A number of different statistical techniques have been used in an effort to improve the correlation between AOD and PM. Gupta and Christopher (2008) used multiple linear regression (MLR) and neural network techniques to include planetary boundary layer (PBL) height, location, temperature, and relative humidity (RH) in an analysis in the southeastern United States and achieved an improved r ² of 0.7 using MLR and 0.83 using the neural network technique compared with r ² of 0.41 for a simple linear regression model. Pelletier et al. (2007) applied generalized additive model (GAM) techniques to retrieve PM from the combination of AERONET (AErosol RObotic NETwork) AOD and National Center for Environmental Prediction (NCEP) meteorological data and compared with surface PM collected at two sites in northwestern Europe. Their model increase r ² to 0.76 compared with an r ² of 0.27 derived from a simple linear regression. Vidot et al. (2007) used the Pelletier et al. GAM method to relate surface PM with SeaWiFS (SeaWide Field-of-view Sensor), AOD, and NCEP meteorological parameters, achieving an r ² of 0.61 compared with 0.48 for the linear case. Chudnovsky et al. (2012) have applied a mixed-effects model approach to AOD retrievals from the Geostationary Operational Environmental Satellite (GOES) to predict PM_2.5 concentrations over New England. These studies motivated our use of GAM methods to improve estimates of PM only using satellite observations.

This paper demonstrates that combining observations from multiple satellite using advanced statistical techniques can dramatically improve the estimates of PM. These results are obtained using the remote satellite observations only, without surface measurements or input from meteorological models. Specifically, our technique uses GAMs to improve retrievals of PM from a combination of A-Train satellite observations. These improved retrievals can be used to fill in the gaps between surface measurements and validate air quality models that are used for air quality forecasts and epidemiological studies.

Methodology

Study area

Our study concentrated on California's San Joaquin Valley (SJV) because of its importance as an agricultural area and its persistent problems with poor air quality. The SJV is home to over 3 million people and encompasses nearly 64,000 km². Improved large-scale information on PM will be useful in validating the air quality models used in pollution and epidemiology studies in this area. The SJV is located southeast of San Francisco in the center of the state between the Coastal Mountain Range to the west and the Sierra Nevada Range to the east (see Figure 1). This topography leads to northwesterly winds during the spring, summer, and fall and a thick well-mixed boundary layer. During the winter, the winds abate and cold temperatures form an inversion layer keeping the boundary layer low and inhibiting surface mixing. The highest PM_2.5 concentrations are typically found during the winter and fall. The winter and fall are also characterized by high concentrations of NH₄NO₃ (Battye et al., 2003). Fresno, located in the center of the SJV, has a population of more than 500,000 and experiences frequent hospitalizations for asthma (Watson, 2000). From 1991 to 1996, PM_2.5 annual averages ranged from 18 to 24 μg/m³, with the highest 24-hr averages ranging from 56 to 93 μg/m³. The 24-hr averages exceeded the National Ambient Air Quality Standards (NAAQS) for PM_2.5, which are 15.0 μg/m³ for the annual average and 35 μg/m³ for the 24-hr average.

Figure 1. Map of the San Joaquin Air Pollution Control District in California showing the ground stations used in this study.

Elevated concentrations of PM_2.5 can be associated with about 18,000 (ranging from 5600 to 32,000) premature deaths in California each year based on the 2004–2006 air quality data (Tran et al., 2009). Ground-based measurements are sparse in California and do not provide enough coverage for air quality monitoring. Further, there is little or no correlation between MODIS AOD and surface PM measurements in this region using single-variate, least-squares linear regressions (Engle-Cox et al., 2004). Previous attempts to improve this correlation by taking into account surface reflectance included in the MODIS Deep Blue retrieval for bright surfaces (Hsu et. al., 2006) versus the standard MODIS dark-target retrieval models have met with little success (Ballard et al., 2008; Justice et al., 2009).

Data used

Surface data.

 The PM_2.5 surface measurements for the years 2004–2008 came from sites maintained and operated by the California Air Resources Board (CARB) (Table 1, Figure 1). Air quality indicators such as chemical compounds and PM concentrations, ozone levels, and meteorological conditions are monitored at these sites (Bucsela et al., 2008). The sites used for this study are listed in Table 1 (also see Figure 1).

Table 1. List of sites used in this study and their latitude and longitude

Site	Latitude	Longitude
Bakersfield	35.36˚N	−119.06˚W
Fresno	36.78˚N	−119.77˚W
Modesto	36.64˚N	−120.99˚W
Stockton	36.95˚N	−121.27˚W
Tracy	36.68˚N	−121.44˚W
Turlock	36.48˚N	−120.84˚W

In this study, we use PM_2.5 obtained using both the Federal Reference Method (FRM) gravimetric-filter-based samplers for daily PM and the beta attenuation monitor (BAM) for the hourly PM. The FRM samplers are the only samplers currently approved by the EPA to provide data for determining the attainment status for an area. The BAM continuous mass monitors are used to collect data for understanding the diurnal and episodic behavior of fine particles, transport assessment, and for use by health scientists investigating exposure patterns. BAM measurements are a confirmed U.S. Federal Equivalence Measurement relative to the 24-hr FRM standard (Zhu et al., 2007). The BAM determines the deposited mass by the attenuation of high-energy electrons through the sample filter. Accuracy is estimated to be ±8% of indication for hourly readings and ±2% for daily readings when compared with the FRM measurements (Ecotech, 2006).

Hourly PM_2.5 and daily PM₁₀ averages were obtained using the organization's online “Real-Time Query Tool”: http://www.arb.ca.gov/aqmis2/aqdselect.php All available relative humidity data, collected since 2002, were gathered using CARB's “Meteorological Data Query Tool”: http://www.arb.ca.gov/aqmis2/metselect.php

Satellite data.

 The satellite data used in this study come from the Ozone Monitoring Instrument (OMI) onboard National Aeronautics and Space Administration's (NASA's) EOS-Aura (Earth Observing System-Aura) satellite and the MODIS onboard NASA's EOS-Aqua (Levelt et al., 2006). Aqua and Aura were launched into a sun-synchronous orbit on May 4, 2002, and July 15, 2004, respectively. Aqua observes the SJV between 1:00 and 2:00 p.m., Pacific Standard Time. There is an approximate 8-min time lapse between the Aqua and Aura overpasses. MODIS AOD was used in this study because of its good spatial resolution, effective clouds mask algorithms, and its near-daily global coverage.

MODIS acquires daily global measurements with 36 spectral bands (from 0.41 to 14 μm) at three different spatial resolutions (250 m, 500 m, and 1 km). AOD is reported at 10 × 10-km² resolution. MODIS high spatial and temporal resolutions are well suited for evaluating air quality on regional and global scales (Levy et al., 2007a, 2007b). The MODIS data products used in this study were the Level II Aerosol Product, Corrected Optical Depth Land, Surface Reflectance Land, Quality Assurance Land, Deep Blue Aerosol Optical Depth Land, Deep Blue Surface Reflectance Land, and Optical Depth Land and Ocean. Both the standard MODIS dark-target and the Deep Blue algorithm require a surface reflectance value to retrieve AOD values. The standard MODIS AOD retrieval algorithm uses the dark-target approach over land. The land surface reflectance parameterization is based on an empirical relationship between surface reflectances measured in the visible channels and at 2.13 µm. The relationship is a function of geometry and a mid-infrared (IR) vegetation index measuring the greenness of the vegetation (Levy et al., 2007a). This algorithm is accurate for determining AOD above highly vegetated surfaces, but performs poorly over highly reflective surfaces such as cities and deserts. The SJV surface reflectance is a patchwork of bright and dark targets. We studied the effects of both parameters on estimating PM_2.5 in this paper. The Deep Blue algorithm has proven to be successful at deriving AOD values over bright surfaces such as deserts (Hsu et al., 2006). Deep Blue uses a set of lookup tables based on a polarized radiative transfer model, allowing it to simulate the radiance for “a range of solar and viewing geometries at the top of the atmosphere” (Hsu et al., 2004). The MODIS dark-target AOD product lies within an expected error envelope of ±(0.05 + 15%) at 550 nm (Levy et al., 2010).

The OMI is an imaging spectrometer that measures solar light backscattered by the Earth's atmosphere and surface (Bucsela et al., 2006). The instrument consists of two spectrometers, one measuring the ultraviolet (UV) spectral range from 270 to 365 nm in two subranges (UV1: 270–314 nm, resolution: 0.42 nm, sampling: 0.32 nm; UV2: 306–380 nm, resolution: 0.45 nm, sampling: 0.15 nm), the other measuring the UV-visible spectrum from 350 to 500 nm (resolution: 0.63 nm, sampling: 0.21 nm). OMI uses a charge-coupled diode (CCD) array with one dimension resolving the spectral features and the other dimension allowing a 114˚ field of view, providing a 2600-km viewing swath transverse to the orbit track. Its nadir spatial resolution ranges from 13 × 24 to 24 × 48 km², depending on the instrument's operating mode and observing geometery. For this study, we used OMI aerosol optical depth, aerosol absorption optical depth (AAOD), and NO₂ data product version 3. Good agreement is generally seen between the OMI NO₂ measurements and ground-based measurements, with tropospheric columns underestimated by 15–30% (Celarier et al., 2008).

For this study, time-coordinated data sets from October 2004 to July 2008 were used. Satellites provide AOD usually only once a day and correlations are typically better with hourly rather than daily PM_2.5 surface measurements (Gupta and Christopher, 2008; Ballard et al., 2008; Justice et al., 2009). Gupta and Christopher (2008) found that linear correlations increased from 0.52 to 0.62 when hourly data were used instead of daily. In this study, the surface PM_2.5 measurement was taken at the nearest hour to the satellite overpass.

Relating PM and AOD

For cases where there are no significant aerosol plumes above the boundary layer, PM can be related to AOD (Koelemeijer et al., 2006) through the relationship

1

where H is the well-mixed boundary layer over which the AOD is measured, f(RH) is the ratio of ambient and dry extinction coefficients, Q _ext,dry is the Mie extinction efficiency (a function of particle size, composition, and wavelength), ρ is the particle density, and r _eff is the particle effective radius (the ratio of the third and second moments of the size distribution). Equation 1 shows that the relationship between PM and AOD is dependent on the vertical distribution of particles in the atmosphere, their size and composition, and the efficiency with which they interact with light.

Ground-based measurements relate to the aerosol PM that is resident in the PBL. PBL height is a function of the local meteorology. For example in the SJV, winters are characterized by small PBL, resulting in high concentrations of PM near the surface. In the summer, the PBL is larger, resulting in decreased concentrations of PM near the surface due to dilution. The presence of f(RH) in eq 1 shows that the surface RH also affects AOD because the uptake of water vapor onto particles increases their size and extinction efficiency. In this study, these meteorological effects are represented with a “seasonality” parameter, θ, that is modeled as a simple function of “day of year,” Admittedly, θ will not take into account episodes where large amounts of particles reside above the PBL, affecting AOD but not PM. Incorporating these occurrences in the present algorithm would likely improve results; however, information on the vertical distribution of particles in the atmosphere are not routinely available in the SJV from ground-based lidars, and data from the space-based Calipso lidar is too infrequent to be useful for this study. Finally, particle- and gas-phase constituents are interrelated through the processes of nucleation, condensation, and evaporation; therefore, a relationship between the satellite gas-phase product and the PM response can reasonably be expected. In this study, column concentrations of NO₂ from the OMI satellite were found to be significant in the retrieval of PM.

Statistical methods

In a simple linear regression model, the relationship between a response y and parameters x ₁, …, x _p is given by

2

where the B_j are constants and ϵ is the residual error. In the simplest case, when PM is compared with AOD only, eq 2 becomes

3

Linear regressions can be extended to include multiple parameters. In many cases, however, the relationships between the parameters and the response are not linear. Additive models are used to construct a regression function from smooth additive functions of the determining parameters, the f_j (x_j ). In an additive model, we replace the linear relationship ( eq 1) between the parameters and response with a functional relationship so that an additive model takes the form

4

where the f_j are simple, smooth arbitrary functions replacing the constant coefficients B_j . Replacing the B_j with f_j is necessary to capture any nonlinear dependency of PM_2.5 with other parameters. Some parameters may have a large effect on the response in only a particular part of their range. For example, a variable may be dominated by random error at very low values or have other undesirable effects at very high values, such as MODIS Deep Blue AOD as shown in Figure 4a. For these reasons, it is appropriate to make the simple extension of linear to nonlinear regressions. GAMs are such an extension, retaining additive terms (Wood, 2006, 2007).

We explored many combinations of parameters looking for the best fit of “retrieved” to measured PM. Both the standard MODIS AOD data product and MODIS AOD derived using the Deep Blue algorithm (DB_AOD) were used because linear regressions of measured PM_2.5 with DB_AOD gave marginally better correlations than using the standard MODIS product (Justice et al., 2009). Preliminary studies showed a good correlation between PM_2.5 and surface NO_x (NO_x = NO + NO₂) measurements indicating that particulate pollution from transportation contributed to PM_2.5 concentration. Because the goal of this study was to use only remote observations, OMI NO₂ was used as a proxy for surface NO_x and this proxy gave promising results. OMI aerosol optical depth (OMI_AOD) and aerosol absorption optical transmittance (OMI_AAOD) were included because they are derived from observations at shorter wavelengths and may be more sensitive to aerosols at higher altitudes. Finally, a seasonal parameter (Θ) was used to account for season effects on relative humidity and PBL height. With this full set of parameters, eq 4 becomes

5

The requirement that all of these parameters were available for a given data set resulted in a small sample size. This is discussed in section “Effect of parameter selection on sample size.” Models for use in the GAM were constructed using subsets of these parameters to explore the tradeoff between improvements in retrieved PM and sample size. These models are shown in Table 2. Equation 5 corresponds to model M. The notation, s(x_j , k = n) designates the use of a spline function with n knots to model the effect of parameter x_j . The spline coefficients for all these models are listed in Supplementary Materials.

Table 2. GAMs used in this study

Symbol	Formula
M	PM25_Hourly ∼ s(AOD_47, k = 4) + s(DB_AOD_47, k = 4)
	+ s(OMI_AOD, k = 4) + s(OMI_AAOD, k = 4)
	+ s(OMI_NO2, k = 3) + s(θ, bs = “cc”, k = 5) + 1
M–St	PM25_Hourly ∼ s(DB_AOD_47, k = 4) + s(OMI_AOD, k = 4)
	+ s(OMI_AAOD, k = 4) + s(OMI_NO2, k = 3)
	+ s(θ, bs = “cc”, k = 5) + 1
M–DB	PM25_Hourly ∼ s(AOD_47, k = 4) + s(OMI_AOD, k = 4)
	+ s(OMI_AAOD, k = 4) + s(OMI_NO2, k = 3)
	+ s(θ, bs = “cc”, k = 5) + 1
M–AA	PM25_Hourly ∼ s(AOD_47, k = 4) + s(DB_AOD_47, k = 4)
	+ s(OMI_AOD, k = 4) + s(OMI_NO2, k = 3)
	+ s(θ, bs = “cc”, k = 5) + 1
M–AA,–St	PM25_Hourly ∼ s(DB_AOD_47, k = 4) + s(OMI_AOD, k = 4)
	+ s(OMI_NO2, k = 3) + s(θ, bs = “cc”, k = 5) + 1
M–AA,–DB	PM25_Hourly ∼ s(AOD_47, k = 4) + s(OMI_AOD, k = 4)
	+ s(OMI_NO2, k = 3) + s(θ, bs = “cc”, k = 5) + 1
Notes: PM25_Hourly is the hourly surface measurement of PM2.5. ∼ is the R code symbol denoting that the response on the left is related to the parameters on the right; see eq 4. AOD_47 is the MODIS AOD retrieved using the standard algorithm. DB_AOD_47 is the MODIS AOD retrieved using the Deep Blue algorithm. OMI_AOD is the OMI AOD (or aerosol optical thickness). OMI_AAOD is the OMI aerosol absorption optical thickness. OMI_NO2 is the OMI NO2 measurement. Θ is the seasonal parameter. s(xj , k = n) designates the use of a spline function with n knots to model the effect of parameter xj . s(θ, bs = “cc”, k = 5) designates a cyclical spline function with 5 knots to model the seasonal parameter.

The implementation used for the GAMs in this study was contained in the mgcv package of the R-project code (R Development Core Team, 2008; Wood, 2006, 2011). For all models, 80% of the data were randomly selected as a training data set for the GAM. The resultant model was tested on the remaining 20% of the data. This procedure was followed five times and the results averaged to ensure that the results were not affected by the random selection of the points. Two standard checks were employed to ensure that results were not influenced by some fortuitous aspects of the data set (see Wood, 2008; Faraway, 2002). All forms given for the spline terms were checked using generalized cross-validation (the GCV of routine mgcv), which is a mathematically precise and informative method of selective data denial. Also, all terms, simple and spline-fit, were tested using visual analysis of partial residuals, which reveals any idiosyncrasies of the samples.

Results and Discussion

Linear regressions between parameters

When a particular response, in this case PM, is a function of multiple parameters, it is instructive to explore the relationship between different parameters. This information can guide the selection of parameters to include in a multiple regression and help explain the results of the analysis. Excellent correlations exist between AOD retrieved by the same method or instrument at different wavelengths. For example, the correlation between MODIS AOD at 0.47 and 0.55 μm using the standard retrieval (AOD_47 and AOD_55) is 0.97. The correlation for the Deep Blue retrieval (DB_AOD_47 and D_AOD_55) was also 0.97. The high correlation between AOD at different wavelengths may be due to the constraints that the MODIS retrieval puts on the wavelength dependence. MODIS AOD at 0.55 μm is derived from the measurements at 0.47 and 0.66 μm; however, the close correlation between the 0.47 and 0.55 μm values and the retrieval constraints indicate that no useful additional information is to be had by considering multiple wavelengths. This practice is typical of other satellite retrievals and indicates that no new information will be obtained from including both wavelengths from the same retrieval. Correlations for the linear regressions between the parameters of interest in this study are shown in Table 3. It is interesting that the AOD retrieved using different methods does not correlate very well: the correlation between AOD_47 and DB_AOD_47 was 0.35; the correlation between AOD_47 and OMI AOD was 0.15; and the correlation between DB_AOD_47 and OMI AOD was 0.007. The correlation between OMI AOD and AAOD was 0.55. These low correlations are surprising. They suggest that the different retrieval algorithms and inputs play an important role in determining AOD. They also suggest that different retrievals or observations may provide different information about the aerosol optical depth. This may be because the measurements are made at different wavelengths (OMI observations are made at shorter wavelengths than MODIS) or different assumptions are made for aerosol models or surface reflectance are used in the retrieval algorithm. The other surprising fact from the data in Table 3 is the generally low level of correlation between the other parameters. Low correlations between PM and satellite observations in this region have been observed and commented on before. It is possible that observations from different satellites may represent different geographical areas or dependence on their respective grid and pixel coverage.

Table 3. Correlation coefficient (r ²) matrix for variables used in regressions

Variable	PM2.5_Hourly	AOD_47	DB_AOD_47	OMI_AOD	OMI_AAOD	OMI_NO2	Day of Year
PM2.5_Hourly	1.00	0.12	0.16	0.004	−0.06	0.05	0.01
AOD_47		1.00	0.35	0.15	0.05	0.07	0.01
DB_AOD_47			1.00	0.007	0.01	0.01	0.01
OMI_AOD				1.00	0.55	−0.02	−0.03
OMI_AAOD					1.00	−0.01	0.02
OMI_NO2						1.00	0.04
Day of year							1.00

The use of a mixed-effects model is suggested because this data set contains several different kinds of AOD. A mixed-effects model approach is a generalization of the standard linear model that incorporates both fixed and random effects, and can account for repeated data sets. The authors explored a mixed-effects modeling approach, using the GAMM (generalized additive mixed model) function within the R code mgcv package (Wood, 2011); however, no particular advantage was found over the GAM. To understand why, one needs to look at the details of the “repeated measures” that were used in the data set. The standard MODIS AOD was obtained using a retrieval algorithm that got surface reflectance using the standard dark-target model. The MODIS Deep Blue AOD was obtained using an algorithm that got surface reflectance using the Hsu et al. (2006) method. The OMI AOD was obtained using radiances from a different instrument and a retrieval algorithm different from either algorithm used by MODIS. One can think of these different AODs as different interpretations of the same physical phenomenon. This is borne out by the poor correlation coefficients between the AODs shown in Table 3.

Correlations for Fresno site

As an example of how the GAM can work for an individual site, we present here an analysis for the Fresno site. Figure 2a shows a scatter plot of PM_2.5 derived from a linear regression of MODIS Deep Blue AOD (DB_AOD_47) and ground-based measured PM_2.5. The ordinate is labeled as “modeled PM_2.5” and the abscissa, “measured PM_2.5.” The correlation coefficient, r ², is 0.12. Note the lack of agreement especially for large values of PM_2.5. These are the most important points to match because they are associated with the exceedance of EPA criterion pollutant standards, which is 35 μg/m³ for a 24-hr average PM_2.5. Figure 2b shows a scatter plot of PM_2.5 modeled using a GAM versus the measured PM_2.5 for the same data. The correlation coefficient has improved to r ² = 0.72 for PM_2.5 measurements; therefore, for the Fresno site, the use of the GAM showed an improvement in the retrieved PM_2.5 over the simple linear regression.

Figure 2. Scatter plots of retrieved PM_2.5 versus measured PM_2.5 for the Fresno site. (a) PM_2.5 retrieved from a linear regression to MODIS AOD versus measured PM_2.5 and (b) PM_2.5 retrieved from GAM M_–AA–DB versus measured PM_2.5. The linear regression line is shown as a solid line and the 1:1 line is shown as a dashed line for reference.

Multisite retrieval

Data from all of our SJV sites are combined into one data set. A site-specific parameter was sued to differentiate between sites to account for topological and population differences, but it had a negligible contribution to the retrieval and it was not included in the parameter sets considered here.

An examination of the data shows that there are many more data points for low PM_2.5 events than high PM_2.5 events. This causes the GAM to put more weight on the lower PM events; however, it is the higher PM_2.5 events (the exceedances) that are of greater importance to the community, regulators, and epidemiologists. For the purposes of this study, we chose to weight the response using a weighting factor related to the value of PM_2.5 itself. In the general retrieval process, PM_2.5 is not known a priori so weighting was accomplished using an iterative process. First the GAMs were run without weighting. These results can be seen in Table 4 in the row titled “r ² for GAM” and range between 0.55 and 0.59 for the various models used. This is a great improvement over the simple linear regression shown in the row titled “r ² for linear,” which range from 0.08 to 0.18. The resultant retrieved PM values were used to establish a weighting and the GAMs were run again. These data are shown in Table 4 in the row titled “r ² for weighted GAM” and show a marked improvement over the unweighted GAM. Values of r ² for the 6 models varied from 0.74 to 0.77. The model that includes all of the parameters, M__AA__DB, gave the best performance.

Table 4. Correlation coefficients, number of valid points, and parameter significance using 3 × 3 MODIS grids for several GAMs compared with a linear regression and with the GAMs weighted as described in the “Multisite retrieval section.

Parameter	M	M–St	M–DB	M–AA	M–AA,–St	M–AA,–DB
r ^2 for GAM	0.59	0.57	0.57	0.58	0.56	0.55
r ^2 for weighted GAM	0.77	0.76	0.75	0.76	0.74	0.74
r ^2 for linear	0.08	0.16	0.18	0.10	0.17	0.17
No. of Points	549	871	1306	852	1298	2070
St AOD	****	X	****	****	X	****
DB AOD	****	****	X	****	****	X
OMI AOD	**	**	**	****	****	****
AAOD	****	*	****	X	X	X
NO2	*	****	****	****	****	****
θ	****	****	****	****	****	****
Notes: Significance codes used were significance as Bayesian P values; Nychka (1988); and R mgcv.pdf package documentation February 2012, http://cran.r-project.org/web/packages/mgcv/mgcv.pdf, also available in R package history. The P value categories do indicate low sampling variability of our estimates, even when parameter influence is small, reflecting the size of the samples (Schervish, 1996). X = not included in GAM formula; ****P < 0.001; ***P < 0.01; **P < 0.1; *P < 1.

Taylor diagrams (Taylor et al., 2001; Taylor, 2008) are used to visually compare the performance of statistical models, taking into account the standard deviation and correlation coeffiecent. A Taylor diagram is a polar plot with the standard deviation as the radial axis and the correlation coefficient as the angular axis. Figure 3 is an example of a Taylor diagram for the six models described in Table 2. The point of the abscissa is the standard deviation of the reference measurement, in this case the measured PM_2.5. The curved dashed lines represent the distance from the reference point and indicate the root mean square difference between the model and reference measurement. Figure 3 shows that all six GAMs used in this study lie very close to each other, demonstrating that there is not much difference in the performance of these models. We chose to concentrate on model M_{– AA, – DB} because it includes the most sample points.

Figure 3. Taylor diagram comparing the performance of the six GAMs of Table 2.

Figure 4 shows scatter plots of modeled/measured PM_2.5 concentrations for three different cases. Figure 4a is a scatter plot of PM_2.5 retrieved from a linear regression to MODIS AOD versus measured PM_2.5. Figure 4b is a scatter plot of PM_2.5 retrieved from GAM M_–AA,–DB versus measured PM_2.5, and Figure 4c represents the same data plotted in weighted space. When using a weighting function, it is appropriate to look at the correlation as a function of the weighted variables. A standard correlation between two variables, x and y, is

Figure 4. Scatter plots comparing retrieved/measured PM_2.5 for data points valid with the M_–AA–DB model. (a) PM_2.5 retrieved from a linear regression to MODIS AOD versus measured PM_2.5; (b) PM_2.5 retrieved from GAM M_–AA–DB versus measured PM_2.5; and (c) PM_2.5 retrieved from the weighted GAM, M_–AA–DB, versus measured PM_2.5 weighted space. The linear regression line for the retrieved/measured PM_2.5 is shown as a solid line.

6

where μ is the mean value and σ is the standard deviation of the parameter x or y. The weighting function, w_n , used was based on the retrieve values

7

to insure that the weighting function was normalized. The weighted correlation becomes

8

Figure 4c shows the scatter plot for the weighted observed and retrieved PM_2.5. It is clear that the weighted results matched a larger number of the higher values of PM_2.5.

Parameter sensitivities

From Table 4 it is clear that OMI NO₂ was a significant contributor to predicting PM_2.5. This is probably because OMI NO₂ is a measure of how much the SJV is affected by NO_x-generating polluters, in this case, vehicle traffic. NO₂ is formed in the environment from primary emissions of oxides of nitrogen. Although there are natural sources of NO_x (e.g., forest fires), the combustion of (fossil) fuels has been, and remains, the major contributor in the SJV. Traffic pollution and power plants are two of the biggest sources of NO₂ pollution in the United States. Sulfate aerosol make up only about 5% of the total aerosol in the Fresno area, whereas nitrate makes up about 25% and black carbon (BC) and organic carbon (OC) make up about 18% and 21%, respectively (McMurry et al., 2004). Nitrate is often a major component, and substantial conversion from NO₂ to nitrate may occur within 30–60 km, or one to two OMI footprint regions (Valin et al., 2011). Locally measured PM_2.5 and local nitrogen oxide measurements correlate geographically, responding to strong local urban (traffic) sources. In turn, local nitrogen oxide measurements correlate with satellite-derived NO₂ column concentrations.

Figure 5 shows the relationships between the response (PM_2.5) and the parameters (e.g., AOD, NO₂). Figure 5a shows the contribution that the MODIS Deep Blue AOD (DB_AOD_47) parameter has on the response, PM_2.5. This contribution is the coefficient, f_j in eq 4, which corresponds to the parameter MODIS Deep Blue AOD. The contribution is monotonically increasing as expected. The comb on the x-axis shows where the value of DB_AOD_47 data points lie. There is a high concentration of points in the middle range from about 1.75 to 2.75. At the very low and high values of AOD, the contribution to PM is no longer linear. In the low-AOD range, the contribution to PM_2.5 is smaller, probably reflecting the sensitivity of the instrument. In the high-AOD range, PM_2.5 is more sensitive to DB_AOD_47.

Figure 5. Plot of smoothed splines used to account for the contribution of PM_2.5 due to (a) MODIS Deep Blue AOD (DB_AOD_47) and (b) seasonality parameter, θ. These sensitivity plots give a partial contribution to the PM_2.5 in μg m⁻³, which adds to other contributions.

Meteorological conditions, such as vertical mixing in the PBL, RH, long- and short-range transport, and temperature and sunlight effects on the formation of secondary organic aerosol strongly affect PM_2.5 concentrations. Previous studies have shown that there is a seasonal variation of PM_2.5 in the SJV. The largest values of PM_2.5 at the surface are measured in winter when wind speeds are low and cold temperatures form an inversion layer keeping the boundary layer low and inhibiting surface mixing (Watson et al., 2000, and references cited within; Justice et al., 2009). Others have considered the affects of pollution sources on the air quality in the region (e.g., Chow et al., 2006), and found that emissions from motor vehicles are essentially constant during the year. Because PM_2.5 is a dry measurement and AOD is not, the seasonal variation in RH will also be reflected in the seasonal functional dependence.

Our premise was to use only remote observations to retrieve surface PM_2.5 values and no reliable observations of PBL or RH could be obtained from remote sources. Attempts to use PBL and RH from assimilations models were fruitless because these proved to be inaccurate unless meteorology was included in the assimilation. The usefulness of our retrievals for validating model predictions of PM_2.5 would be compromised if the same meteorological model was used in estimating PBL and RH. Thus, introducing a seasonality parameter was the logical way to account for the strong seasonal variation in PM_2.5 in the SJV and the day-of-year parameter was the best choice for this study. It is not uncommon to take seasonal variation into account in retrieving PM_2.5 from remote measurements. For example, Liu et al. (2009a) discriminated between spring and nonspring seasons in their work using GAMs and model-generated meteorological fields.

One limitation with using day-of-year parameter is that the effects of long- and short-range transport at high altitude are not taken into account. Investigating the relationship between sunphotometer-derived AOD to PM at Ispra, Italy, Barnaba et al. (2010) found that maximum AOD occurred in spring/summer and maximum PM values occurring in winter. They were able to use coincident lidar measurements to demonstrate that the variation in aerosol vertical distribution was important in this discrepancy. Justice et al. (2009) also found that surface measurements of PM_2.5 in SJV were greatest in the winter, whereas MODIS AOD values were lowest in the winter. However, because there are no lidar sites in the SJV, we were unable to explore the affect of the vertical resolution of aerosol in this study. The lack of information on the vertical structure of aerosol in the SJV remains a serious limitation to the study of aerosol in the SJV and for improved retrievals of PM_2.5 from satellite observations.

Figure 5b shows that the seasonal parameter (day-of-year) has a strong effect on retrieved PM_2.5 and Table 4 shows that it is statistically significant. Figure 5b shows that the effect of the seasonal parameter on PM_2.5 is stronger during the winter season (days 0–50 and 300–356) and weaker during the rest of the year (days 100–250). This result is not surprising given the seasonal dependence of PM_2.5 in the SJV. Most studies take the seasonal dependence into account to some degree to have success in retrieving PM (Gupta and Christopher, 2008; Liu et al., 2009a).

Effect of parameters selection on sample size

The number of “valid” data points for a particular model is dependent upon the parameters in the model. In order for a point to be “valid” for a particular model, all parameters in the model must have a valid value. The number of valid data points that can be used with a particular model is dependent upon the number and kind of parameters in the model. For example, in model M (see Table 2), where all six parameters are used, the number of valid points is 549, whereas the number of valid points increases to 817 in model M__St when only the MODIS standard AOD is dropped. Table 5 compares the results of the six models listed in Table 2 for the same points as in the most restrictive model, M, with the largest number of parameters. Predictably, the r ² for the various models are the same or higher than for cases where all available points are used.

Table 5. Effect of adding observation type on the number of data points with valid data

Observation type		No. of Days
A	Days in series	8244
B	Total MODIS days	4455
C	Total valid 3 × 3 MODIS days	4340
D	Days in row C without PM25_Hourly	1299
E	Days in row C without MODIS _AOD	97
F	Days in row C without MODIS_DB_AOD	1993
G	Days in row C without OMI_NO2	230
H	Days in row C without OMI_AOD	531
I	Days in row C without OMI_AAOD	1955

The number of sample points decreases dramatically when satellite data are included. Hidy et al. (2009) note that for MODIS data in Fresno from 2002 to 2008, cloud-free days were available 43%, 61%, 94%, and 81% of the time in winter, spring, summer, and fall, respectively. Table 5 shows how the number of sample points decreases as more parameters are added. The data set is limited to days that have OMI data, from October 1, 2004, to July 6, 2008. This gives a total of 8244 possible samples for six sites. Of these possible samples, MODIS data existed for 4455 points. The criterion for assigning a valid 3 × 3 pixel was that >50% of the pixels had to have valid data. This restriction is more severe than is typically used for these types of studies and (slightly) reduced the number of points to 4340. The reduction in points resulting from the merge with satellite data is largely due to cloudy pixels. Because OMI has a larger footprint than MODIS, the number of “cloudy” pixels is larger.

The effect of the different GAMs on the number of valid points can also be seen in Table 4. The GAM with the most parameters gives the best results; however, the number of points is small. By not including two of the most restrictive parameters, OMI_AAOD and DB_AOD, the number of points using GAM M_–AA–DB has increased nearly 4 times to 2070 and the r ² has decreased to 0.74 for the case of the 3 × 3 MODIS grid. As discussed, adding more parameters significantly reduced the sample size. In this case the sample size was reduced from 2070, when only four parameters were considered, to 549 when all six parameters were considered. Thus, the GAM M_–AA–DB gave almost the same performance as the six-parameter model, M (see Figure 3), and considered nearly four times the valid data points. The significance of the parameters is shown in Table 4 by significance codes. It can be seen in Table 4 that when the parameters used in a GAM are changed, the functions (e.g., f_j in eq 4) and significance of the new set of parameters also change. This would not occur if all the predicting variables were orthogonal to each other, but they rarely are (Draper and Smith, 1998; Wood, 2006).

Retrieved and observed PM trends

To further illustrate the utility of the multisatellite GAM retrieval, the measured PM_2.5 is compared with the retrieved values and plotted versus the day of the study period in Figure 6. Figure 6a shows the measured PM_2.5 as black dots and the retrieved PM_2.5 as gray dots. The retrieved PM_2.5 were obtained using the results of a linear regression of MODIS AOD to PM_2.5. The black and gray lines are loess fits to the measured and retrieved points, respectively. A straight dot-dashed line marks the EPA 24-hr criterion for reference. It is seen that the retrieved PM_2.5 using the linear regression model does not capture the trend of the measured values and never does reach the USEPA criterion.

Figure 6. Trends of retrieved and measured PM_2.5 in SJV. The gray and black points are the retrieved and measured PM_2.5, respectively. The gray and black lines are loess fits to the retrieved and measured PM_2.5, respectively. (a) Linear model; (b) GAM. The EPA 24-hr criterion is shown as a dot-dashed line.

In Figure 6b, the black dots are once again the measured PM_2.5 and the gray dots represent the retrieved values of PM_2.5 from the GAM M_–AA–DB. The black and gray lines are loess fits to the measured and retrieved points, respectively, and the dot-dashed line marks the EPA 24-hr criterion. It is seen that the retrieved PM_2.5 from the GAM more closely matches the measured values. Seasonal trends can be seen in the retrieved GAM values that match those of the measured values. In particular, many exceedances are captured. Of the 2070 data points in the set, the surface measurements recorded 149 exceedances. The linear fit recorded 6 and the GAM fit recorded 167. The GAM correctly identified 69% of the exceedances and 14% of the retrieved exceedances were false positives. Although this is a big improvement, it underscores the fact that there is a long way to go before surface measurements will be replaced by satellite retrievals for regulatory purposes; however, we contend that the real near-term utility of satellite-derived PM will be in filling in the gaps of surface measurements to improved validation of air quality models and for epidemiology studies.

Conclusion

This paper shows that using GAMs with multiplatform satellite observations can improve the correlation between observed and retrieved PM_2.5 in California's SJV. PM_2.5 retrieved using a simple linear regression with MODIS AOD yields correlation coefficients between 0.08 and 0.18 for this region. The GAM implementation used was from the mgcv routine of R, a commonly available, free statistical program. The parameters used are MODIS AOD (standard dark-target algorithm, MODIS Deep Blue AOD, OMI AOD, OMI AAOD, OMI NO₂ concentration, and a seasonal parameter. Using the GAM techniques and combinations of these parameters, the correlation coefficients between retrieved and measured PM_2.5 ranged between 0.74 and 0.77. The best retrieval used all of these parameters, r ² = 0.77; however, the number of data points was low, 549. The model that included MODIS standard dark-target AOD, OMI AOD, OMI NO₂, and a seasonal parameter gave a correlation coefficient of 0.74 but increased the number of points to 2070. Particularly noteworthy is the fact that the PM_2.5 retrieved using the GAM captures many of the PM_2.5 exceedances that were not captured with the simple linear regression. It was also found that the MODIS grid of 3 × 3 yielded better results than either a 1 × 1 or a 5 × 5 grid. This grid size compared well with the OMI pixel size.

Further improvements to any new technique are always needed. The GAMs and input data can be refined by using other combinations of parameters and including more measurement sites, especially rural sites. Certainly, the inclusion of other available data, especially information on the vertical distribution of particles, when available, may further improve these results. The question of the general applicability of this technique also needs to be addressed in the future. We intend to apply this technique to other areas globally that have similar topography to the SJV and that have demonstrated poor correlations between AOD and PM_2.5. We expect that the combinations of parameters used in the SJV will be useful in these areas. For example, the meteorology, surface conditions, and pollution sources particular to the SJV can also be found in other areas of the planet, such as the Po Valley in Italy, China's Red Basin, and the Indo-Gangetic Plain. In other areas, it is likely that the parameters will have to be changed to improve correlation of retrieved and measured PM_2.5; however, the methodology presented here should be applicable. We feel that the techniques demonstrated in this study can be used to greatly enhance the utility of satellite-retrieved PM_2.5 for air quality purposes. The GAM allows one to see the relationships between the parameters and the response that may lead to improved understanding of the physical phenomena involved.

Improved retrieved PM_2.5 from satellite data will have profound benefits. They will be valuable in validating emission inventories used in global climate models and for validation of air quality models. They will be especially useful in sparsely populated regions where no data are available and areas where residents are exposed to significant PM_2.5 concentrations but where measurements are not available. Such data can be used to track pollution transport, suggest the placement of future surface monitoring stations, and support epidemiological studies that seek to identify the sources of the most toxic air pollutants and the susceptible populations. These techniques will prove useful in improving satellite observations of particulate as well as gas-phase air quality, not only in the western United States, but globally.

Acknowledgment

MODIS data were obtained from the Goddard Earth Sciences (GES) Data and Information Services Center (DISC) and OMI data were obtained from the OMI International Science Team for the OMI satellite data. The authors wish to acknowledge the MODIS and OMI teams for the use of their data products. Drs. Lorraine Remer and Richard Kleidman gave the authors valuable insight on the use of MODIS aerosol data products and on the techniques used in this study. Dr. Omar Torres gave the authors valuable insight on the use of the OMI aerosol data products. The comments of Ms. Meytar Sorek-Hamer and Dr. Raymond Hoff were very valuable.

The authors acknowledge the personnel of CARB for the effort required to collect the surface data used in this study. Much of the surface measurements used in this study were compiled by two teams of students working at Ames under the NASA DEVELOP program under the leadership of Dr. J. Skiles and Ms. C. Schmidt during the summers of 2009 and 2010. Dr. A. Strawa was the mentor for these groups.