Improving the Accuracy of Daily PM2.5 ... - ACS Publications

Apr 4, 2016 - Fusion of Ground-Level Measurements with Aerosol Optical Depth. Observations, a Case Study in North China. Baolei Lv,. †,‡. Yongtao ...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/est

Improving the Accuracy of Daily PM2.5 Distributions Derived from the Fusion of Ground-Level Measurements with Aerosol Optical Depth Observations, a Case Study in North China Baolei Lv,†,‡ Yongtao Hu,§ Howard H. Chang,⊥ Armistead G. Russell,*,§ and Yuqi Bai*,†,‡ †

The Ministry of Education Key Laboratory for Earth System Modeling, Center for Earth System Science, Tsinghua University, Beijing 100084, China ‡ Joint Center for Global Change Studies, Beijing 100875, China § School of Civil and Environmental Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332, United States ⊥ Department of Biostatistics and Bioinformatics, Emory University, Atlanta, Georgia 30322, United States S Supporting Information *

ABSTRACT: The accuracy in estimated fine particulate matter concentrations (PM2.5), obtained by fusing of stationbased measurements and satellite-based aerosol optical depth (AOD), is often reduced without accounting for the spatial and temporal variations in PM 2.5 and missing AOD observations. In this study, a city-specific linear regression model was first developed to fill in missing AOD data. A novel interpolation-based variable, PM2.5 spatial interpolator (PMSI2.5), was also introduced to account for the spatial dependence in PM2.5 across grid cells. A Bayesian hierarchical model was then developed to estimate spatiotemporal relationships between AOD and PM2.5. These methods were evaluated through a city-specific 10-fold cross-validation procedure in a case study in North China in 2014. The cross validation R2 was 0.61 when PMSI2.5 was included and 0.48 when PMSI2.5 was excluded. The gap-filled AOD values also effectively improved predicted PM2.5 concentrations with an R2 = 0.78. Daily ground-level PM2.5 concentration fields at a 12 km resolution were predicted with complete spatial and temporal coverage. This study also indicates that model prediction performance should be assessed by accounting for monitor clustering due to the potential misinterpretation of model accuracy in spatial prediction when validation monitors are randomly selected.



INTRODUCTION Fine particulate matter, a complex mixture of particles with aerodynamic diameters of 2.5 μm or less (PM2.5), has been a major component of the severe air pollution levels experienced in many cities in China in recent years. PM2.5 is the dominant pollutant of concern, especially during the winter.1−3 Since PM2.5 can efficiently penetrate into human lungs and bronchi,4,5 long-term and short-term exposures to PM2.5 can increase premature mortality and morbidity.5 One major limitation with ground-level measurements of PM2.5 is the sparse spatial coverage of monitoring networks. Alternatively, satellite-based aerosol optical depth (AOD) observations have been used to estimate continuous spatial and temporal patterns of ground PM2.5 concentrations.6−8 Optical properties of atmospheric aerosols are sensitive to the optical properties of the particles that make up particulate matters, including size and composition. A large fraction of the particles that make up PM2.5 particle mass scatter sunlight, leading to a strong association between the mass of PM2.5 and observed AOD.9,10 All else being equal, as the number (and mass) of particles increases, light scattering increases, leading to © XXXX American Chemical Society

an increase in AOD. While the association between PM2.5 and AOD is typically strong, the relationship can change due to changes in the particle size distribution, particle composition, mixing height, humidity, and other factors. Early studies used various linear regression modeling techniques to estimate PM2.5 concentrations using AOD as well as other spatial and spatial-temporal predictors, such as land cover, elevation, meteorological parameters, and indicators for holidays.10,11 More complex models were then developed, including the linear mixed effect (LME) model,12,13 geographically weighted regression (GWR) model,14,15 remote sensing formula,16 semisupervised learning approach based on multiple factors,17 and complex ensemble models.18 Atmospheric chemistry models were also used to determine the relationships between AOD and PM2.5.8,19 Among all these studies, treating the linear slopes and intercepts between AOD Received: December 10, 2015 Revised: March 23, 2016 Accepted: April 4, 2016

A

DOI: 10.1021/acs.est.5b05940 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology

Figure 1. Study area and the locations of PM2.5 monitoring stations as shown by points. The right panel depicts the average satellite-based AOD values in 2014.

μg/m 3 in 2014, about nine times the World Health Organization (WHO) recommended standard of 10 μg/m3. Long-term exposure to such heavily polluted air is a major concern for more than 250 million people in this area.32,33 Many atmospheric chemistry modeling applications in this region, such as the use of the Community Multiscale Air Quality model (CMAQ), have been run at a 12 × 12 km grid resolution on a Lambert Conformal Conic map projection. Therefore, the same grid was chosen in this study to enable the comparison of and fusion with atmospheric chemistry model output in future studies. The whole area was divided into 6561 grid cells (81 × 81). Ground-Level PM2.5 Data. Hourly ground-level PM2.5 concentrations were obtained from the China National Urban Air Quality Real-time Publishing Platform (http://113.108.142. 147:20035/emcpublish/). Calibration and quality control of the monitors are conducted by the China National Environmental Monitoring Center.34 Two hundred ninety-eight monitors from 53 cities in the study area are located in 169 grid cells, described in Table S1, Supporting Information (SI). To facilitate computation in the data fusion process, at each time point we averaged PM2.5 concentrations from monitors that are located in a same grid cell. It is worth noting that point measurements of PM2.5 concentrations were used to represent the area-average PM2.5 levels over a grid cell. This may cause potential bias due to a change of support. Specifically, the gridaverage concentrations may not fully reflect the spatial variation in PM2.5 as captured by multiple monitors within a grid cell. MODIS AOD Data. MODIS is a sensor on board two of the U.S. National Aeronautics and Space Administration (NASA)’s Earth Observation System (EOS) satellites: Terra and Aqua. The sun-synchronous satellites provide two column aerosol observations at approximately 10:30 a.m. (Aqua) and 1:30 p.m. (Terra) local solar times every day. The MODIS AOD values are constrained between −0.05 and 5.0 to avoid bias introduced during the retrieval procedure. A detailed description of the MODIS AOD retrieval procedure is discussed in Remer et al.35 and Levy et al.21 In this study, MODIS AOD Level 2 data in 2014 were obtained from the Level 1 and Atmosphere Archive and Distribution System (LAADS, https://ladsweb.nascom. nasa.gov/). The nominal resolution of this data set is 10 × 10 km at nadir. A nearest neighbor approach was utilized to regrid the data to the 12 × 12 km setting. MODIS AOD Missing Data Pattern. For each of the 169 grid cells in which at least one ground monitoring station is available, averaged PM2.5 concentrations were computed among

and PM2.5 at different locations and at different time points as spatially and temporally correlated random effects has been shown to be promising in a study conducted in the southeast United States.20 However, to our knowledge, such a statistical model structure has not been previously applied in China. AOD data derived from instruments aboard satellites, such as the moderate resolution imaging spectroradiometer (MODIS), are often missing due to clouds, high surface albedos21 (e.g., snow and ice), and potentially, very high PM2.5 pollutant levels.22,23 If the missing AOD values are related to high PM2.5 pollution levels, the average PM2.5 concentrations estimated using only the retrieved AOD fields can be underestimated. Several attempts have been made to solve this issue. Kloog et al. filled in missing AOD values by smoothing the mean PM2.5 levels across the study area.24,25 Van Donkelaar et al. used the sampling bias correction factor method.7,26 Their calculations were usually month or year specific,19,27 such that the bias could not be effectively corrected when calculating the mean in a small number of days. Just et al. predicted PM2.5 levels within the grid cells with missing AOD by using PM2.5 levels in the neighboring grid cells, based on season-specific spatial patterns.28 Ma et al. used ordinary Kriging to interpolate the available retrieved AOD values.14 However, their method could have significant uncertainties if only a limited number of AOD values were present.29 In this study, we applied the Bayesian model proposed by Chang et al.20 in North China to estimate the spatially and temporally varied coefficients in a linear regression setting. We also considered using spatially interpolated PM2.5 concentrations as a predictor and filled the missing AOD values with novel methods. The rest of the paper is organized as follows. The study area, data, statistical model, and missing AOD gapfilling method are introduced in the Materials and Methods. Model estimation, performance evaluation, and PM2.5 prediction results are presented in the Results and Discussion.



MATERIALS AND METHODS Study Area. Figure 1 shows the North China study area. This region includes five provincial administrative divisions: Shandong, Hebei, Liaoning, Beijing and Tianjin. Economic development in this region is highly dependent on heavy industries, such as iron and steel factories and cement production, that are typical sources for emissions of PM2.5 and its precursors.30 Increasing vehicular population and electricity demand also worsen the air quality in this region.31 The annual mean PM2.5 concentrations in this region was 93 B

DOI: 10.1021/acs.est.5b05940 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology

Figure 2. (a) Average PM2.5 concentrations; (b) daily SVVR within the study area.

weights for predicting PM2.5 levels at unobserved locations. The uncertainties (prediction standard errors) of the interpolated PM2.5 concentration field are not spatially uniform, and the closer to the monitors a grid cell is, the lower its uncertainty will be.29 To ensure that the uncertainties of the interpolated values were generally uniform, a 10-fold leave-10%-cities-out method was used to generate the interpolated PM2.5 (PMI2.5) value using the OK method. Specifically, we first randomly removed the monitors in 10% of 53 cities within the study area. PMI2.5 values at the removed sites were obtained from kriging the monitors from the remaining cities. This procedure was repeated 10 times to obtain interpolated PMI2.5 values at all monitors. To further control the uncertainties with the variable PMI2.5, we defined the final PM2.5 spatial interpolator (PMSI2.5) as follows

days when AOD observations were available, among days when AOD values were missing, and across the whole year of 2014 (plotted, respectively, in black, red, and blue points in Figure 2a). The AOD averages were plotted in an ascending order based on the mean PM2.5 concentrations in that grid cell. It is clear that mean PM2.5 concentrations were consistently higher on days when AOD observation were missing than when AOD data were present (Figure 2a). To quantitatively measure the spatial and temporal pattern of AOD missing values, SVVR (spatial valid value ratio, defined in SI) and TVVR (temporal valid value ratio, defined in SI) were used. Daily SVVR values, shown in Figure 2b, indicate that missing AOD was more significant in winter than in other seasons. The annual TVVR was generally lower among grid cells in the northwest grassland areas and urban areas (Figure S1). In winter, the TVVR was quite low (as low as 10%) in most grid cells, e.g., in heavily polluted south Hebei. Because the missing AOD values were associated with higher pollution levels (Figure 2a), negative bias in long-term PM2.5 estimates will be present if only the retrieved AOD values are used in data fusion. MODIS AOD Gap Filling. To overcome the problem of missing AOD data, we proposed a novel two-step method. First, a city- and season-specific formula was introduced, based on a linear assumption of the relationship between PM2.5 and AOD:10 AODei , j = as

PMi , j PMs , j

AODos , j + bs

PMSI 2.5 = PMI 2.5

PMI 2.5 × f (PMI 2.5) σt = PMI 2.5 σ(PMI 2.5) SE(PMI 2.5) (2)

and f (x) = 10% ×

1

(

x α

)

+1

(3)

where σ(PMI2.5) and σt refer to, respectively, the uncertainties of PMI2.5 and the target uncertainties in that grid cell. The function f is a continuous inverse proportional function to constrain the uncertainties of PMI2.5 to be a relatively small proportion of the original interpolated values (PMI2.5). The proportion begins with 10% and decreases with increasing PMI2.5. The parameter α determines rate of decrease. When PMI2.5 equals α, f is 5%. In this study, α was determined to be 700 based on an optimization between decreasing uncertainties and decreasing variance information on PMSI2.5 as α increases. Other Spatiotemporal Predictors. First, we obtained meteorological data from the NCEP (National Centers for Environmental Prediction) FNL (Final) Operational Global Analysis that are on 1° × 1° grids and are prepared operationally every 6 h. This product is from the Global Data Assimilation System (GDAS), which continuously collects observational data from the Global Telecommunications System (GTS) and other sources. In this study, we used the parameters of ground temperature (T), relatively humidity (RH), and planetary boundary height (PBL) at 6:00 UTC (14:00 Beijing local standard time). The second data source included land cover and elevation data. Land cover data were obtained from the MODIS land-cover classification of the

(1)

AODei,,j refers to the estimated AOD value on day i in grid cell j. PMi,j refers to the daily mean PM2.5 concentration observed on day i in grid cell j. AODos,j refers to the average retrieved AOD values in grid cell j and in season s, where s = 1 denotes warm season from Apr 16−Oct 15 and s = 2 denotes cold seasons as the remaining days in a year. PMs,,j refers to the average PM2.5 concentrations in grid cell j in season s. The linear relationship was fitted separately in 53 cities. The cityspecific models were then used to predict missing AOD values in grid cells containing PM2.5 monitors. Second, we combined the estimated and retrieved AOD values to interpolate missing AOD values in the other grid cells without PM2.5 monitors using ordinary Kriging (OK) with exponential covariance function and obtained a final AOD data set with complete spatial coverage for each day. Interpolated Ground PM2.5 Data. We proposed a novel variable based on the interpolation of PM2.5 observations to be used in combination with AOD. A spatial interpolation utilizes the spatial dependence among observations to obtain optimal C

DOI: 10.1021/acs.est.5b05940 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology

model formulation based on out-of-sample prediction. We designed a workflow (depicted as the Figure S2) to evaluate the three methods using both 10-fold and spatial cross-validation experiments. For the 10-fold cross-validation, the data were randomly divided into 10 groups. For a particular experiment, nine groups were used to fill the missing AOD values and to fit two AOD− PM2.5 models (eq 4), with and without the variable PMSI2.5. The remaining group served as the test data set where the predicted PM2.5 values were compared to the observed measurements. The PM2.5 concentrations in the test data set could be predicted by either the retrieved AOD values or the filled AOD values. The evaluation was done by withholding a different group each time for 10 times. We performed a 10-fold spatial cross-validation by leaving 10% of the cities out each time or 10% of the monitors out each time. We considered leaving all observations out by cities instead of by individual monitors to address a weakness that has been reported in the literature.14,27 Specifically, given the dense number of monitors in cities, there were typically one or more monitors near any monitor. Hence, removing individual monitors one-at-a-time in a cross-validation experiment may falsely show that the model can perform spatial predictions well. Results from removing monitors by a group of cities, however, provide confidence in the model’s ability to predict at locations far away from monitors. Model prediction performance was evaluated using the normalized mean absolute error (NME), root mean square error (RMSE), R2, and the linear regression coefficients (along with their standard errors). NME and RMSE are defined in the SI.

International Geosphere-Biosphere Program (IGBP) product in 2010 (http://lpdaac.usgs.gov). The data set was assigned to the base grid. We used two parameters: green fractions and urban fractions. Finally, the mean elevation from the USGS GTOPO30 elevation data set (https://lta.cr.usgs.gov/ GTOPO30) for each grid cell was obtained. Statistical Model Framework. A linear regression framework was assumed to model the relationships between the MODIS AOD values and ground PM2.5 measurements20 (eq 4) PM(s , t ) = α0(s , t ) × AOD + α1(s , t ) + ε(s , t )

(4)

where PM(s,t) denotes the PM2.5 concentrations measured at location s and time t. Location s is the index of a monitorcontained grid cell, which is georeferenced as in Figure 1, and time t is the index of a day within the study period. The location and time specific slope and intercept are, respectively, designated as α0(s,t) and α1(s,t). The residual errors ε(s,t) are assumed to be normally distributed with mean zero and variance σ2. Model coefficients have both a temporal and spatially dependent structure. The coefficients of the linear regression are separated into a second level α0(s , t ) = β0(s) + β0(t ) + λ0 Z0

(5)

and α1(s , t ) = β1(s) + β1(t ) + λ1Ζ1

(6)

where βi(s) and βi(t), respectively, refer to the underlying independent spatial and temporal processes that determine the spatiotemporal variations of the slope (for i = 0) and intercept (for i = 1). The parameter vectors λi represent fixed effects for meteorological, land cover variables: Z0 for the slopes (for i = 0) and Z1 for the intercepts (for i = 1). Specifically, eq 5 models the relationships between AOD and PM2.5, with potential interactions between AOD and other predictors, and eq 6 models the underlying processes that have direct effects on the PM2.5 levels. PMSI2.5 is a predictor variable included in Z1 to capture addition variation in PM2.5 not explained by AOD and the other predictors. In this study, we used a tapered conditional Kriging method to interpolate the spatial coefficients βi(s) in eq 2 for spatial dependence.20 For those days without PM2.5−AOD pairs, we used a linear autocorrelation function to model the temporal coefficients βi(t), which specifies that βi(t) is proportional to the weighted average between the neighboring days.20 Therefore, we could obtain, with complete coverage, the spatially and temporally dependent coefficient fields through interpolating the fitted AOD field. We used the Markov Chain Monte Carlo (MCMC) method for model fitting under the Bayesian framework. One advantage of employing a Bayesian hierarchical framework is that the model could account for uncertainties in parameter estimates and PM2.5 predictions. This facilitates the model outputs being incorporated in human exposure and health impact evaluations that considers exposure uncertainties.36,37 It is worth noting that the fitted model can only be used to predict AOD values within the temporal extent of the available data, as day-specific AOD−PM2.5 relationships beyond the time scope (i.e., temporal extrapolation) were not estimated in our model structure. Model Evaluation Workflow. We evaluated the effectiveness and performance of the AOD gap-filling method, the novel variable PMSI2.5 as an additional predictor, and the statistical



RESULTS AND DISCUSSION Model Parameters Estimation. The mean PM2.5 concentrations from monitors for the final data set was 58.8 μg/m3, which is 1.7 times that of the China Class I standard. On average, AOD and PMSI2.5 were significantly associated with their 95% posterior intervals not including zero (Table 1), with Table 1. Statistical Assessment of the Estimated Coefficients terms in λ1

eEstc

SEd

Low95e

Up95f

intercept AOD PMSI2.5 Green.faca Urban.facb elevation temperature PBL RH

23.44 21.66 1.16 0.21 0.64 −0.012 0.04 −0.0013 −0.033

18.80 3.93 0.016 0.31 1.82 0.004 0.062 0.00047 0.016

−13.19 15.34 1.12 −0.4 −2.97 −0.02 −0.084 −0.0022 −0.065

60.99 30.81 1.19 0.82 4.20 −0.0044 0.163 −0.00035 0.000065

a

Green.fac refers to the fraction of green land in each of the grid cells. Urban.fac refers to the fractions of urban area in each of the grid cells. c Est denotes the estimated coefficients. dSE denotes the standard errors of the estimated coefficients. eLow95 denotes the lower bound of 95% posterior interval. fUp95 denotes the upper bound of 95% posterior interval. b

a positive coefficient indicating a positive linear relationship with the PM2.5 levels. The mean uncertainty of the variable PMSI2.5 was less than 5 μg/m3. Hence, due to the presence of PMSI2.5, an average uncertainty of approximately 5 μg/m3 was incorporated into the estimations considering its coefficient being close to 1 (Table 1). Both PBL and terrain height were D

DOI: 10.1021/acs.est.5b05940 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology

captured by the spatiotemporal predictors, including AOD. When using the random leave-10%-monitors-out crossvalidation method, the R2 was 0.68 (Figure 3d), which is similar to the R2 of 0.69 reported by Zheng et al.27 using the leave-one-monitor-out validation method but is 0.07 larger than that using the leave-10%-cities-out cross validation method in this study. Hence, performing cross-validation at the individual monitor level could exaggerate the model’s prediction performance. The annual PM2.5 levels were calculated by averaging the estimated PM2.5 concentrations by retrieved AOD values. The pollution levels in the cities in Hebei province along the foot of the Taihang Mountains were much higher than those in other cities. The full model with PMSI2.5 better characterized this belt-shape pollution zone (Figure S4). Moreover, the full model with PMSI2.5 predicted the lower pollution levels in the northern part of study area than the partial model, with lower prediction standard errors (Figure S4b,d). The standard errors in the urban area were generally lower due to the greater density of AOD−PM2.5 pairs. However, an estimation bias was apparent in the predictions even by the full model. The annual means of the observed PM2.5 concentrations in the heavily polluted Hebei cities were approximately 120 μg/m3. However, the highest estimated PM2.5 concentrations (full model) were no larger than 100 μg/ m3 (Figure S4c). Spatially, the estimated PM2.5 concentrations were not smooth, especially in the winter (Figure S5d). Moreover, the full model using retrieved AOD values failed to estimate the increased pollution area in the southern Hebei in the autumn and winter (Figure S5). In all of these cases, the model predicted PM2.5 concentrations only when AOD values were available. The observed mean PM2.5 levels with available AOD values were lower than the annual mean PM2.5 levels (Figure 2a), which accounts for the underestimations. AOD Gap-Filling Evaluation. There were approximately 40000 missing AOD values in the grid cells with PM2.5 monitors. Our missing AOD imputation method was effective at improving PM2.5 predictions. In the model-fitting procedure, the model overestimated AOD when the retrieved AOD values were low (Figure 4a). Although the R2 was low (0.36), the constructed AOD and retrieved AOD showed a good linear relationship when the retrieved AOD was less than 1.5. By visually checking the interpolated full AOD data field, we found that the spatial variations of the AOD values were well captured. There were almost no abnormal spatial patterns in

negatively associated with PM2.5 levels. Regarding the daily specific random intercepts and slopes (β0(t) and β1(t), respectively, in eq 5 and eq 6), they are larger in winter because the PM2.5 concentrations were higher (Figure S3, SI). Model Comparisons and PMSI2.5 Evaluations. In model fitting, R2 improved by 0.07 (from 0.66 to 0.73) when the variable PMSI2.5 was considered in the model. In the 10-fold leave-10%-cities-out cross validation procedure, the R 2 improved from 0.48 (Figure 3b) to 0.61 (Figure 3c) by a

Figure 3. Model evaluations using the 10-fold leave-10%-cities-out cross-validation method for the (a) Meteo (Meteorological variables) + LU (land use variables) + PMSI2.5 model, (b) AOD + Meteo + LU model, and (c) full model. Model evaluation using the 10-fold leave10%-monitors-out cross-validation method for the full model (d). No AOD gap-filling was performed in these models.

large margin of 0.13. The RMSE also decreased by 3.46 μg/m3, from 26.99 μg/m3 to 23.53 μg/m3. In general, the model with PMSI2.5 predicted PM2.5 more accurately. The variable PMSI2.5 explained more variance of PM2.5 than AOD values by comparing the model’s performance without using AOD values and that without PMSI2.5 (Figure 3a). In summary, PMSI2.5 contains additional spatial information PM2.5 that is not

Figure 4. (a) Scatterplots of the AOD values fitted by eq 1 and retrieved AOD values; (b) prediction performance using the reconstructed AOD with full model; (c) annual average AOD values with gap-filling. E

DOI: 10.1021/acs.est.5b05940 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology

Figure 5. Annual and seasonal mean estimations of the PM2.5 concentrations (μg/m3) using the full model and reconstructed AOD data set.

accounts for spatial clustering. They predicted the mean PM2.5 concentrations in southern Hebei province at around 90 μg/m3 in 2013 and 130 μg/m3 in the winter. In contrast, for 2014, we predicted improved PM2.5 levels to be at around 120 μg/m3 overall and 180 μg/m3 in the winter, which were closer to observations. Zheng et al.27 evaluated their LME model in Beijing−Tianjin−Hebei area with an R2 = 0.77 and NME = 0.22 using the more optimiztic leave-one-monitor-out cross validations method. They did not give estimated seasonal variations, even though they predicted annual mean PM2.5 concentrations well by using a correction factor. Xie et al. predicted PM2.5 levels in Beijing, with an R2 = 0.75 in crossvalidation.13 Their performance being better than that in this study is partially owing to its smaller study area. Without considering missing AOD values, they significantly underestimated annual mean PM2.5 levels. In our study, the data set with complete spatial coverage could well characterize the evolution of serious pollution episodes on a regional scale, and a case pollution episode is presented in the SI (Figure S6). In summary, our revised model and reconstructed AOD data set accurately predicts ground PM2.5 concentrations in North China. First, we used a Bayesian hierarchical framework to model the daily grid cell-specific linear relationships between the AOD and PM2.5 concentrations. Second, we used the spatial dependence in all the observed PM2.5 levels by developing an interpolated PM2.5 variable PMSI2.5. The new variable PMSI2.5 can significantly improve the model’s prediction performance. The R2 was improved by 0.13, to 0.61, in the cross-validation study. Third, we reconstructed the AOD values, and thus, we obtained daily estimated PM2.5 pollution maps with complete spatial and temporal coverage, which was shown to be useful for capturing the evolution of PM2.5 pollution episodes that occurred in North China. This study also demonstrates the potential misinterpretation of model accuracy when a completely random, leave-one (or 10%) cross-validation method is used to evaluate model

the daily AOD. When used in predicting PM2.5, the full model performed well with the reconstructed data set, giving an R2 of 0.78, an NME of 0.27, and an RMSE 33.39 μg/m3 in the crossvalidation experiment (Figure 4b). It is interesting that the cross-validation performance with the constructed AOD was better than the model with retrieved AOD. In part, this is because the linear relationship between AOD and PM2.5 enabled our model to leverage additional extra spatial information from PM2.5 monitors. PM2.5 Concentration Estimation. The model using the “full” AOD data set predicted the seasonal and annual average of the PM2.5 concentrations well (Figure 5). The PM2.5 concentrations were much higher in winter compared to those using the retrieved AOD only. The mean PM 2.5 concentration in winter was approached 200 μg/m3. The heavy pollution mainly occurred in the middle and southern Hebei province, the western Shandong province, Tianjin, and Beijing. The PM2.5 levels were also high in the middle of Liaoning province around Shenyang, its capital city. This area is part of the Northeast plain, which is traditionally a zone of heavy industry. The PM2.5 concentrations were generally lower than 40 μg/m3 in northwest Hebei Province and in the Inner Mongolia portion within the domain, which is mainly covered by grassland and forests. It is worth noting that there is a small area of higher PM2.5 levels in the east of Inner Mongolia and near Liaoning Province. The area is the location of an industrial city, Chifeng, which is the second most polluted city in Inner Mongolia, with an annual mean PM2.5 level of nearly 50 μg/m3. The successful prediction of this pollution hotspot among clean surroundings further demonstrated the model’s effectiveness. The PM2.5 levels were also lower in Shandong and the Liaodong peninsula and in the middle south of Shandong Province. A previous study by Ma et al.14 used the GWR method with 2 R = 0.64 in the 10-fold cross validation method, while we obtained an R2 = 0.61 by our stricter evaluation process that F

DOI: 10.1021/acs.est.5b05940 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology

space and the ground: effectiveness of emission control. Environ. Sci. Technol. 2010, 44 (20), 7771−7776. (7) Van Donkelaar, A.; Martin, R. V.; Levy, R. C.; da Silva, A. M.; Krzyzanowski, M.; Chubarova, N. E.; Semutnikova, E.; Cohen, A. J. Satellite-based estimates of ground-level fine particulate matter during extreme events: A case study of the Moscow fires in 2010. Atmos. Environ. 2011, 45 (34), 6225−6232. (8) Van Donkelaar, A.; Martin, R. V.; Park, R. J., Estimating groundlevel PM2.5 using aerosol optical depth determined from satellite remote sensing. J. Geophys. Res. 2006, 111, (D21).10.1029/ 2005JD006996 (9) Kahn, R.; Banerjee, P.; McDonald, D.; Diner, D. J. Sensitivity of multiangle imaging to aerosol optical depth and to pure-particle size distribution and composition over ocean. J. Geophys. Res. 1998, 103 (D24), 32195−32213. (10) Liu, Y.; Sarnat, J. A.; Kilaru, V.; Jacob, D. J.; Koutrakis, P. Estimating ground-level PM2.5 in the eastern United States using satellite remote sensing. Environ. Sci. Technol. 2005, 39 (9), 3269− 3278. (11) Wang, J.; Christopher, S. A., Intercomparison between satellitederived aerosol optical thickness and PM2.5 mass: implications for air quality studies. Geophys. Res. Lett. 2003, 30, (21).10.1029/ 2003GL018174 (12) Lee, H.; Liu, Y.; Coull, B.; Schwartz, J.; Koutrakis, P. A novel calibration approach of MODIS AOD data to predict PM 2.5 concentrations. Atmos. Chem. Phys. 2011, 11 (15), 7991−8002. (13) Xie, Y.; Wang, Y.; Zhang, K.; Dong, W.; Lv, B.; Bai, Y. Daily estimation of ground-level PM2.5 concentrations over Beijing using 3 km resolution MODIS AOD. Environ. Sci. Technol. 2015, 49, 12280. (14) Ma, Z.; Hu, X.; Huang, L.; Bi, J.; Liu, Y. Estimating ground-level PM2.5 in China using satellite remote sensing. Environ. Sci. Technol. 2014, 48 (13), 7436−7444. (15) Song, W.; Jia, H.; Huang, J.; Zhang, Y. A satellite-based geographically weighted regression model for regional PM 2.5 estimation over the Pearl River Delta region in China. Remote Sensing of Environment 2014, 154, 1−7. (16) Lin, C.; Li, Y.; Yuan, Z.; Lau, A. K.; Li, C.; Fung, J. C. Using satellite remote sensing data to estimate the high-resolution distribution of ground-level PM2.5. Remote Sensing of Environment 2015, 156, 117−128. (17) Zheng, Y.; Liu, F.; Hsieh, H.-P. In U-Air: When urban air quality inference meets big data. Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining, 2013; ACM, 2013; pp 1436−1444. (18) Lary, D. J.; Faruque, F. S.; Malakar, N.; Moore, A.; Roscoe, B.; Adams, Z. L.; Eggelston, Y. Estimating the global abundance of ground level presence of particulate matter (PM2.5). Geospatial Health 2014, 8 (3), 611−630. (19) Geng, G.; Zhang, Q.; Martin, R. V.; van Donkelaar, A.; Huo, H.; Che, H.; Lin, J.; He, K. Estimating long-term PM2.5 concentrations in China using satellite-based aerosol optical depth and a chemical transport model. Remote Sensing of Environment 2015, 166, 262−270. (20) Chang, H. H.; Hu, X.; Liu, Y. Calibrating MODIS aerosol optical depth for predicting daily PM2.5 concentrations via statistical downscaling. J. Exposure Sci. Environ. Epidemiol. 2014, 24 (4), 398− 404. (21) Levy, R. C.; Remer, L. A.; Kleidman, R. G.; Mattoo, S.; Ichoku, C.; Kahn, R.; Eck, T. Global evaluation of the Collection 5 MODIS dark-target aerosol products over land. Atmos. Chem. Phys. 2010, 10 (21), 10399−10420. (22) Tao, M.; Chen, L.; Su, L.; Tao, J., Satellite observation of regional haze pollution over the North China Plain. J. Geophys. Res. 2012, 117, (D12).10.1029/2012JD017915 (23) Engel-Cox, J. A.; Holloman, C. H.; Coutant, B. W.; Hoff, R. M. Qualitative and quantitative evaluation of MODIS satellite sensor data for regional and urban scale air quality. Atmos. Environ. 2004, 38 (16), 2495−2509. (24) Kloog, I.; Chudnovsky, A. A.; Just, A. C.; Nordio, F.; Koutrakis, P.; Coull, B. A.; Lyapustin, A.; Wang, Y.; Schwartz, J. A new hybrid

prediction performance. In our case, using a completely random, exhaustive leave-10%-monitors-out procedure led to an R2, NME, and RMSE of 0.68, 0.26, and 21.40 μg/m3, respectively. Using a process where monitors were removed after being grouped by city led to similar performance statistics of 0.61, 0.28, and 23.53 μg/m3, respectively. In cases where the air quality monitors are not uniformly distributed (as is almost always the case because monitors are preferentially located in cities), the random removal of monitors will typically lead to having one or more other monitors near the one removed. Hence, the values to be predicted can often be captured by measurements from an adjacent monitor that was not removed. Thus, model prediction performance should be assessed by accounting for geographical monitor clustering.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.est.5b05940. Variables NME, RMSE, SVVR, and TVVR, Table S1, and Figures S1−S6 (PDF)



AUTHOR INFORMATION

Corresponding Authors

*Tel: +1-404-894-3079. E-mail: [email protected]. *Tel: +86-10-62795269. E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS Y.H. and A.G.R.’s work for this publication was made possible in part by funding from the USEPA under Grant Nos. RD834799, RD833866, and RD835217. H.C.’s work for the Bayesian model was supported by National Institutes of Health Grant No. R21ES023763. Its contents are solely the responsibility of the grantee and do not necessarily represent the official views of the supporting agencies. Further, the US Government does not endorse the purchase of any commercial products or services mentioned in the publication. This study was also supported by the State Environmental Protection Key Laboratory of Sources and Control of Air Pollution Complex (No. SCAPC201406) and by Tsinghua University (Nos. 20131089277 and 553302001).



REFERENCES

(1) Wang, L.; Wei, Z.; Yang, J.; Zhang, Y.; Zhang, F.; Su, J.; Meng, C.; Zhang, Q. The 2013 severe haze over southern Hebei, China: model evaluation, source apportionment, and policy implications. Atmos. Chem. Phys. 2014, 14, 3151−3173. (2) Zhang, R.; Jing, J.; Tao, J.; Hsu, S.-C.; Wang, G.; Cao, J.; Lee, C.; Zhu, L.; Chen, Z.; Zhao, Y. Chemical characterization and source apportionment of PM2.5 in Beijing: seasonal perspective. Atmos. Chem. Phys. 2013, 13 (14), 7053−7074. (3) He, K.; Yang, F.; Ma, Y.; Zhang, Q.; Yao, X.; Chan, C. K.; Cadle, S.; Chan, T.; Mulawa, P. The characteristics of PM2.5 in Beijing, China. Atmos. Environ. 2001, 35 (29), 4959−4970. (4) Nel, A. Air pollution-related illness: effects of particles. Science 2005, 308 (5723), 804−806. (5) Pope, C. A., III; Dockery, D. W. Health effects of fine particulate air pollution: lines that connect. J. Air Waste Manage. Assoc. 2006, 56 (6), 709−742. (6) Lin, J.; Nielsen, C. P.; Zhao, Y.; Lei, Y.; Liu, Y.; McElroy, M. B. Recent changes in particulate air pollution over China observed from G

DOI: 10.1021/acs.est.5b05940 Environ. Sci. Technol. XXXX, XXX, XXX−XXX

Article

Environmental Science & Technology spatio-temporal model for estimating daily multi-year PM 2.5 concentrations across northeastern USA using high resolution aerosol optical depth data. Atmos. Environ. 2014, 95, 581−590. (25) Kloog, I.; Koutrakis, P.; Coull, B. A.; Lee, H. J.; Schwartz, J. Assessing temporally and spatially resolved PM2.5 exposures for epidemiological studies using satellite aerosol optical depth measurements. Atmos. Environ. 2011, 45 (35), 6267−6275. (26) Van Donkelaar, A.; Martin, R. V.; Pasch, A. N.; Szykman, J. J.; Zhang, L.; Wang, Y. X.; Chen, D. Improving the accuracy of daily satellite-derived ground-level fine aerosol concentration estimates for North America. Environ. Sci. Technol. 2012, 46 (21), 11971−11978. (27) Zheng, Y.; Zhang, Q.; Liu, Y.; Geng, G.; He, K. Estimating ground-level PM2.5 concentrations over three megalopolises in China using satellite-derived aerosol optical depth measurements. Atmos. Environ. 2016, 124, 232−242. (28) Just, A. C.; Wright, R. O.; Schwartz, J.; Coull, B. A.; Baccarelli, A. A.; Tellez-Rojo, M. M.; Moody, E.; Wang, Y.; Lyapustin, A.; Kloog, I. Using high-resolution satellite aerosol optical depth to estimate daily PM2.5 geographical distribution in Mexico City. Environ. Sci. Technol. 2015, 49 (14), 8576−8584. (29) Oliver, M. A.; Webster, R. Kriging: a method of interpolation for geographical information systems. International Journal of Geographical Information System 1990, 4 (3), 313−332. (30) Zhang, Q.; Streets, D. G.; Carmichael, G. R.; He, K.; Huo, H.; Kannari, A.; Klimont, Z.; Park, I.; Reddy, S.; Fu, J. Asian emissions in 2006 for the NASA INTEX-B mission. Atmos. Chem. Phys. 2009, 9 (14), 5131−5153. (31) Lv, B.; Zhang, B.; Bai, Y. A systematic analysis of PM2.5 in Beijing and its sources from 2000 to 2012. Atmos. Environ. 2016, 124, 98−108. (32) Yang, L.; Cheng, S.; Wang, X.; Nie, W.; Xu, P.; Gao, X.; Yuan, C.; Wang, W. Source identification and health impact of PM2.5 in a heavily polluted urban atmosphere in China. Atmos. Environ. 2013, 75, 265−269. (33) Madaniyazi, L.; Nagashima, T.; Guo, Y.; Yu, W.; Tong, S. Projecting Fine Particulate Matter-related Mortality in East China. Environ. Sci. Technol. 2015, 49 (18), 11141−11150. (34) Jiang, J.; Zhou, W.; Cheng, Z.; Wang, S.; He, K.; Hao, J. Particulate matter distributions in China during a winter period with frequent pollution episodes (January 2013). Aerosol Air Qual. Res. 2015, 15 (2), 494. (35) Remer, L. A.; Kaufman, Y.; Tanré, D.; Mattoo, S.; Chu, D.; Martins, J. V.; Li, R.-R.; Ichoku, C.; Levy, R.; Kleidman, R. The MODIS aerosol algorithm, products, and validation. J. Atmos. Sci. 2005, 62 (4), 947−973. (36) Gryparis, A.; Paciorek, C. J.; Zeka, A.; Schwartz, J.; Coull, B. A. Measurement error caused by spatial misalignment in environmental epidemiology. Biostatistics 2009, 10 (2), 258−274. (37) Szpiro, A. A.; Paciorek, C. J.; Sheppard, L. Does more accurate exposure prediction necessarily improve health effect estimates? Epidemiology (Cambridge, Mass.) 2011, 22 (5), 680.

H

DOI: 10.1021/acs.est.5b05940 Environ. Sci. Technol. XXXX, XXX, XXX−XXX