《遥感地学分析》课程教学资源（文献资料）Satellitebased mapping ofdaily highresolution ground PM 2.5 in China viaspacetime regression modeling.pdf_大学文库

Contents lists available at ScienceDirect Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse Satellite-based mapping of daily high-resolution ground PM2.5 in China via space-time regression modeling Qingqing Hea,b , Bo Huanga,b,c,⁎,1 aDepartment of Geography and Resource Management, The Chinese University of Hong Kong, Hong Kong b Big Data Decision Analytics (BDDA) Research Centre, The Chinese University of Hong Kong, Hong Kong c Institute of Space and Earth Information Science, The Chinese University of Hong Kong, Hong Kong ARTICLE INFO Keywords: Aerosol optical depth Geographically and temporally weighted regression (GTWR) Interior point algorithm (IPA) PM2.5 China ABSTRACT The use of satellite-retrieved aerosol optical depth (AOD) data and statistical modeling provides a promising approach to estimating PM2.5 concentrations for a large region. However, few studies have conducted high spatial resolution assessments of ground-level PM2.5 for China at the national scale, due to the limitations of high-resolution AOD products. To generate high-resolution PM2.5 for the entirety of mainland China, a combined 3 km AOD dataset was produced by blending the newly released 3 km-resolution Moderate Resolution Imaging Spectroradiometer (MODIS) Dark Target AOD data with the 10 km-resolution MODIS Deep Blue AOD data. Using this dataset, surface PM2.5 measurements, and ancillary information, a space-time regression model that is an improved geographically and temporally weighted regression (GTWR) with an interior point algorithm (IPA)- based efficient mechanism for selecting optimal parameter values, was developed to estimate a large set of daily PM2.5 concentrations. Comparisons with the popular spatiotemporal models (daily geographically weighted regression and two-stage model) indicated that the proposed GTWR model, with an R2 of 0.80 in cross-validation (CV), performs notably better than the two other models, which have an R2 in CV of 0.71 and 0.72, respectively. The use of the combined 3-km high resolution AOD data was found not only to present detailed particle gradients, but also to enhance model performance (CV R2 is only 0.32 for the non-combined AOD data). Furthermore, the GTWR's ability to predict PM2.5 for days without PM2.5-AOD paired samples and to generate historical PM2.5 estimates was demonstrated. As a result, fine-scale PM2.5 maps over China were generated, and several PM2.5 hotspots were identified. Therefore, it becomes possible to generate daily high-resolution PM2.5 estimates over a large area using GTWR, due to its synergy of spatial and temporal dimensions in the data and its ability to extend the temporal coverage of PM2.5 observations. 1. Introduction Fine particulate matter, a complex of solid and liquid particles suspended in the air with aerodynamic diameters of 2.5 μm or less (PM2.5), has been associated with adverse effects on public health by many epidemiological studies (Kampa and Castanas, 2008; Madrigano et al., 2013; Pope and Dockery, 2006). The rapid urbanization and industrialization of China has led to PM2.5 becoming a dominant factor in air pollution, especially over urban areas, and thus an unprecedented issue of public concern (Bi et al., 2014; Che et al., 2014; Guo et al., 2009). Given the advent of severe PM2.5 levels in China, it is urgent to assess the health effects of PM2.5 exposure across China, using PM2.5 data with a high spatial resolution. However, currently, such an assessment is seriously hindered by the limited measurements available from the sparse network of surface monitoring stations. Because satellite remote sensing has the capacity to provide data with large (even global) spatial coverage, satellite-derived aerosol optical depth (AOD) information is an alternative method of inferring ground-level PM2.5 concentrations (Engel-Cox et al., 2004; Wang and Christopher, 2003). To better investigate PM2.5 exposure for air pollution assessment and epidemiological research, satellite-derived PM2.5 concentrations at high spatial resolution are needed. High-resolution studies in this research field have been increasingly carried out in North America using 1-km Multiangle Implementation of Atmospheric Correction (MAIAC) AOD (Hu et al., 2014; Kloog et al., 2014). However, equivalent studies for China are scarce, due to the lack of high-resolution satellite-derived AOD products. Studies on a national scale have only presented PM2.5 variations with resolutions ranging from 10 to https://doi.org/10.1016/j.rse.2017.12.018 Received 13 March 2017; Received in revised form 2 December 2017; Accepted 15 December 2017 ⁎ Corresponding author at: Department of Geography and Resource Management, The Chinese University of Hong Kong, Hong Kong. 1 Co-first author. E-mail address: bohuang@cuhk.edu.hk (B. Huang). Remote Sensing of Environment 206 (2018) 72–83 0034-4257/ © 2017 Elsevier Inc. All rights reserved. T

50 km, employing widely used AOD products, e.g., the Moderate Resolution Imaging Spectroradiometer (MODIS) Level 2 10 km AOD product and Multi-angle Imaging SpectroRadiometer (MISR) Level 2 17.6 km AOD product (Fang et al., 2016; Ma et al., 2014; Ma et al., 2016). In early 2014, a 3-km aerosol product was included in the new MODIS Collection 6 (C6). It is the first daily global aerosol dataset with a relatively high spatial resolution (Remer et al., 2013). The fine aerosol gradients can benefit air quality applications, i.e., improve the performance of statistical models and provide details of spatial variation for fine particulate matter (Xie et al., 2015). However, the MODIS 3 km AOD data were retrieved using the Dark Target (DT) algorithm that derives AOD over dark surfaces, which may cause a large number of missing values in the daily AOD image, especially for a large area with a complex land surface (He et al., 2017). In other words, it is impossible to generate a daily PM2.5 prediction at high spatial resolution for the entirety of China using only the MODIS 3 km AOD product as the primary predictor (Li et al., 2017; You et al., 2016). Thus, improving the daily coverage of the 3 km MODIS AOD for ground PM2.5 estimation over a large area with multiple types of land cover is essential. Various models have been developed to elucidate the quantitative association between satellite-based AOD and surface PM2.5, ranging from a simple linear regression model (Engel-Cox et al., 2004; Wang and Christopher, 2003) to advanced statistical models such as the linear mixed effect (LME) model (Lee et al., 2011; Xie et al., 2015), land use regression model (Kloog et al., 2011; Kloog et al., 2012), and geographically weighted regression (GWR) model (Hu et al., 2013; Ma et al., 2014; Song et al., 2014; Zou et al., 2016). To date, two main spatial and temporal regression methods have been used to establish the PM2.5-AOD relationship. One is the daily GWR model, which addresses the spatial variability between PM2.5 and AOD for individual days and achieves reasonable results (Ma et al., 2014; Song et al., 2014). However, this model ignores the fact that the PM2.5-AOD relationship could vary with time and may be dependent on previous days. Thus, daily GWR cannot make use of temporal autocorrelation existing in the data and fails to build a model for days with fewer PM2.5- AOD samples. The other method is a two-stage model that can account for both spatial and temporal non-stationarities (Hu et al., 2014; Wu et al., 2016). However, this model calibrates the spatial and temporal PM2.5-AOD relationship separately, i.e., using the LME model to account for its temporal variations in the first stage, and then the GWR model its spatial variations in the second stage. Consequently, this model does not simultaneously treat spatiotemporal variations of the PM2.5-AOD relationship; it is thus not in a better position to predict PM2.5 for days with insufficient or no PM2.5 observations. As such, it becomes difficult to use the two methods to generate high-resolution PM2.5 with high accuracy for a large area when there is a limited number of PM2.5-AOD paired samples. In recent years, the geographically and temporally weighted regression (GTWR) model proposed by Huang et al. (2010), which integrates spatial and temporal distance, has gained increasing popularity in environmental studies (Bai et al., 2016; Chu et al., 2015; Guo et al., 2017). Unlike the two spatiotemporal models already discussed, this space-time regression model can simultaneously incorporate temporal information from the estimation day or from previous days into the spatial variability through a spatial-temporal weighting mechanism. Nevertheless, the predictive power of GTWR for elucidating the relationship between PM2.5 and AOD, especially for days without samples, has yet to be explored. Optimization of parameter values for GTWR is also required to reduce the computational cost for a large dataset. As a corollary, there are still considerable challenges confronting a space-time model, such as GTWR, in exploring the PM2.5- AOD relationship, especially over an extensive area with a large but insufficient sample dataset. The objective of this study is to overcome such challenges by predicting ground PM2.5 at high spatial resolution on a daily basis for the entirety of mainland China, in the setting of a space-time regression model using the MODIS 3 km AOD as the major predictor, and meteorological and land use data as ancillary variables. A customized approach to filling in missing AOD values was devised to improve the availability of daily MODIS 3 km AOD. Subsequently, the GTWR model improved with the parameter optimization approach was developed to generate ground PM2.5 concentrations at high resolution, and its inferring power also explored comprehensively. A daily GWR model and a two-stage (LME + GWR) model were also implemented as benchmarks to examine the extent to which GTWR's simultaneous spatial and temporal weighting established a reliable association between surface PM2.5 and satellite AOD. 2. Data collection and preprocessing 2.1. Ground-level PM2.5 observations Daily average PM2.5 observations from 1 January 2015 to 31 December 2015 were obtained from the China Environment Monitoring Center (http://106.37.208.233:20035/). Ground PM2.5 mass concentrations of mainland China are measured by the tapered-element oscillating microbalance or beta-attenuation method, with calibration and quality control according to national standards GB3095-2012 (China's National Ambient Air Quality Standard, or CNAAQS) (China, 2012) and HJ 618-2011 (Determination of Atmospheric Articles PM10 and PM2.5 in Ambient Air by Gravimetric Method, available at http:// english.sepa.gov.cn/Resources/standards/Air_Environment/). At present, the observations come from 1456 monitoring stations in 329 cities (Fig. 1). Most are located in southeastern China, whereas the northwestern areas have little coverage. 2.2. MODIS AOD data 2.2.1. MODIS 3 km DT AOD products MODIS is a sensor that has been on board NASA Terra since 1999 and Aqua since 2002, providing two daily observations of columnar aerosol properties at approximately 10:30 am (Terra) and 1:30 pm (Aqua) local time. The DT algorithm was originally developed to derive aerosol properties over dark surfaces (e.g., dense vegetated areas) at a 10 km spatial resolution. To satisfy an identified need for monitoring fine-resolution pollution, aerosol datasets with a 3 km resolution (MxD04_3K, x is O for Terra and Y for Aqua) were introduced in the most recently released MODIS C6 in addition to the standard 10 km Level 2 aerosol products (MxD04_L2). The 3 km DT aerosol retrievals use a similar protocol to the 10 km DT aerosol products, but organize 6 × 6 pixels into a retrieval box for inversion, rather than 20 × 20, after all undesirable ones, i.e., cloud and snow/ice pixels, have been screened out and discarded (Remer et al., 2013). MODIS 3 km AOD data at 550 nm (scientific dataset name: Optical_Depth_Land_And_Ocean) in 2015 (downloaded from https://ladsweb.nascom.nasa.gov/) were used to model the PM2.5-AOD relationship. Data obtained before 2015 were used for calibration. 2.2.2. MODIS 10 km Deep Blue (DB) AOD products The MODIS DB algorithm aims to generate AOD retrievals over bright land (e.g., deserts) to fill the gaps left by the DT algorithm (Hsu et al., 2013). Unlike the DT retrieval process, the DB algorithm first produces aerosol properties at a 1 km resolution and then averages the individual retrievals into a 10 × 10 km grid. As with the 3 km AOD data, we calibrated DB aerosol retrievals for 2015 using DB AOD products obtained before 2015, and filled in some gaps where the 3 km DT AOD values were missing with data obtained in 2015. 2.2.3. MODIS AOD combination In MODIS Collection 6, the 3 km AOD data are retrieved using the DT algorithm, and the daily availability of 3 km AOD in the desert-like Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 73

area of northwestern China is very low (even no retrievals at all for many days). The MODIS DB algorithm was originally developed to obtain aerosol loading for bright targets, but it only provides the AOD product at a 10 km resolution. To address this issue, a four-step customized approach of combining the MODIS 3 km DT and 10 km DB AOD products was developed to improve the daily coverage of the MODIS 3 km DT AOD, as follows: (1) Re-gridding the 10 km DB AOD To match the grid cell size of the 3 km DT AOD, the 10 km DB AOD were resampled to a 3-km spatial resolution grid using the cubic convolution resampling technique in ENVI/IDL5.1. (2) Calibrating the MODIS DT/DB AOD According to the global validation of MODIS aerosol products and Fig. S2 (Supplementary information [SI]), there are different retrieval accuracies for DT and DB AODs and thus systematic bias between the DT and DB AOD values (Remer et al., 2013; Sayer et al., 2014). Therefore, to optimize the daily AOD values and fill the gaps between DT and DB AOD retrievals, we calibrated AOD retrievals for 2015 using ground-level aerosol observations from the Aerosol Robotic Network (AERONET) from 2001 to 2014 (Ma et al., 2014). Comparison against AERONET has been widely used to validate satellite-based AOD (Sayer et al., 2014; Tao et al., 2015). A linear analysis was conducted based on the long-term data before 2015 to define the relationship between AERONET and satellite AOD. According to previous studies (Li et al., 2009; Remer et al., 2013), in which the relationship between AERONET and MODIS AODs showed strong seasonality, we fitted linear regressions for each season. These relationships were also separately established for the Aqua 3 km DT, Terra 3 km DT, Aqua resampled-DB, and Terra resampled-DB AOD datasets. These fitted parameters were then used to calibrate AOD values for the four datasets during 2015 (SI, S1). (3) Filling in missing values for Aqua and Terra A simple linear regression was performed for the matched grid cells between Aqua-3 km and Terra-3 km AOD for each day to predict the missing values (i.e., estimation of the Aqua-3 km AOD using the available Terra-3 km value, and vice versa). The same protocol was applied to the resampled-DB AODs to predict the missing AOD values for Aqua and Terra (SI, S2). (4) Combining the four datasets To improve the daily AOD coverage, Aqua/Terra 3 km DT AOD and Aqua/Terra resampled-DB AOD data were merged together for spatiotemporal modeling. We used the same method as with the 3 km and resampled-DB AODs to average the Aqua and Terra values. Because the 3 km DT has higher spatial resolution than the 10 km DB to present finer aerosol gradients, we first averaged the Aqua- and Terra-3 km DT AOD values (both retrievals and predicted values) for each day and then applied the averaged DB data resampled to a 3-km grid to fill pixels when the 3 km AOD values were missing. That is, the averaged DB-resampled AOD values were applied as supplementary material to fill the gap when the 3 km AOD values for both satellites were missing. After the combination, the number of daily AOD values for each pixel increased significantly (see Section 5.2 for more details). When compared to the AERONET AOD observations from 2015, our combined AOD (fraction of retrievals falling within the expected error bounds (F) = 56.25% and RMSE = 0.20 for Aqua calibrated and predicted AOD, and F = 55.07% and RMSE = 0.19 for Terra) shows a Fig. 1. Geographical distribution of PM2.5 monitoring stations included in the study. Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 74

better accuracy than the MODIS operational 3-km DT AOD (F = 46.25%, RMSE = 0.38 for Aqua and F = 49.28%, RMSE = 0.36 for Terra), indicating the effectiveness of our proposed four-step approach. The detailed validation process is presented in S2 (SI). 2.3. Auxiliary data A number of meteorological and land-cover-related variables were used in this study. The meteorological variables were ground-based relative humidity (RH), temperature (T), wind speed (WS), pressure (P), and reanalysis-produced planetary boundary layer (PBL) data; the landcover-related variables are terrain (DEM) and Normalized Difference Vegetation Index (NDVI). Details about these variables are given in S3 (SI). In addition, the 1 km grid population dataset of 2010 (Fu et al., 2010) was used to calculate the population exposure in Section 4.4.1. 2.4. Data integration As the aim of this study was to predict daily average PM2.5 concentrations at 3 × 3 km spatial resolution, a 3-km grid was created (1,073,351 grid cells in total) and exactly matched to the combined AOD pixels. The PM2.5 observations and other variables were reprocessed with a consistent spatial size and temporal interval for integration, because these data vary in spatial and temporal resolution. PM2.5 measurements from multiple stations located within a 3-km grid were averaged. The meteorological station data for each day were interpolated to 3 km spatial resolution using an inverse distance weighted (IDW) technique, using the Spatial Analyst Tools in ArcGIS 10.2. Corresponding to the daily average PM2.5 measurements, the daily mean values obtained from averaging 3-hour PBL data were calculated. The two land-use datasets were separately aggregated to the 3 km grid cells. The final result was 60,810 matched samples incorporating ground PM2.5 and all independent variables, which were distributed over 325 days. 3. Development and validation Over the large area of mainland China, there is considerable spatial variability in the relationship between AOD and PM2.5. In addition, temporal variation significantly affects this relationship, and PM2.5 levels on previous days may affect the concentrations on the day of measurement to a certain degree. Therefore, we developed a space-time regression model, the GTWR model, to simultaneously address spatial heterogeneity and temporal influence on the PM2.5-AOD relationship. To confirm the better performance of GTWR over other relevant models in previous studies, a daily GWR model using the same variables and matched samples as the GTWR model and a two-stage (LME + GWR) model (Hu et al., 2014) were also implemented for detailed comparisons. 3.1. GTWR model development 3.1.1. General structure of the GTWR model Unlike the widely used GWR model, which only takes spatial variation into account when estimating an empirical relationship (Fotheringham et al., 2002), GTWR captures spatiotemporal heterogeneity based on a weighting matrix referencing both spatial and temporal dimensions. In this study, a GTWR model was fitted using the following structure: = + ×+ × + ×+ ×+ × + ×+ × + PM β μ υ t β μ υ t AOD β μ υ t NDVI β μ υ t RH β μ υ t WS β μ υ t T β μ υ t PBL β μ υ t DEM ε ( , ,) ( , ,) ( , ,) ( , ,) ( , ,) ( , ,) ( , ,) ( , ,) i i i i i ii i i ii i i ii i i ii i i ii i i ii i i ii i i 2.5 01 2 34 5 6 7 (1) where PM2.5i is the daily surface PM2.5 concentration of the sample i at location (μi,υi) on day ti; β0 denotes the intercepts at specific location (μi, υi) and day ti; and β1–β7 are the location-time-specific slopes for combined AOD, aggregated NDVI, interpolated RH, WS and T, averaged PBL, and aggregated DEM, respectively. The location (μi, υi) represents the central coordinates of a grid cell in which sample i is located. ti is the day of the year, and only the estimation for that day and previous days (ti ≤ testimation) are used for modeling. εi is the error term for sample i. To estimate the intercept of β0(μi,υi,ti) and the slopes of βi(μi,υi,ti) for each variable, the weight matrix W is introduced to account for the importance of sample j to the estimated parameters of sample i in terms of space and time (Huang et al., 2010). In this study, the widely used Gaussian distance decay-based functions and the combination method of spatial and temporal distance proposed in Huang et al. (2010) were applied to measure the spatiotemporal weight W. The details are presented in S4 (SI). 3.1.2. Exploratory analysis of independent variables Because of the relatively weak relationship between AOD and gridcollocated PM2.5 across China (correlation value of 0.43), meteorological and land-use data were incorporated into the GTWR model as covariates to improve the spatial consistency of the PM2.5-AOD relationships (SI, S3). However, if all parameters mentioned in Section 2 were used as independent variables, the result would exist a collinearity problem. Variance inflation factors (VIF) were calculated as the collinearity diagnostic for the independent variables (combined AOD, NDVI, RH, WS, T, P, PBL and DEM). The results demonstrated that the VIF values of pressure and DEM were very high (VIF > 10), which indicates a strong collinearity problem if all eight independent variables are included in the PM2.5 prediction model (SI, S5). To address this issue, we excluded pressure and kept the seven other variables to fit the final GTWR model (the equation is given in Eq. (1)). 3.1.3. Parameter selection and two GTWR models The spatiotemporal bandwidth hST and scale factor ρ are the two key parameters in GTWR modeling (Huang et al., 2010). In the original GTWR (Huang et al., 2010), the popular method of plotting the sum of squared error against the parameters hST and ρ through cross-validation (Fotheringham et al., 2002), i.e., calculating the scores of CVRSS(hST; ρ) for hST and ρ (Eq. (S2) in SI), was used to select the two optimal values. However, this method is inappropriate for GTWR modeling over a large area, because it requires the use of a nested-loop, which can be easily implemented but consumes a large amount of computation time in search of the two optimal key parameter values (hST and ρ) for our large dataset. Selecting the best hST and ρ is to efficiently find the minimal values for a constrained nonlinear multivariable problem (Huang et al., 2010). Thus, an optimization approach based on the interior point algorithm (IPA) (Byrd et al., 1999; Waltz et al., 2006) was devised in this study. Aided by a self-concordant barrier approach, IPA can expeditiously handle a nonlinear problem by approximately solving barrier sub-problems using a sequential quadratic programming iteration with trust region techniques (Byrd et al., 1999). Compared with other optimization methods such as the active-set algorithm (Nocedal and Wright, 2006), IPA can efficiently handle large, sparse problems, which is essential to a large dataset of PM2.5-AOD pairs. Using the IPA approach to determine the two optimal parameter values in this study, the computational time consumed for our sample dataset was reduced by approximately 60%, from about 2.5 days for the previous nested-loop method to about 1 day (SI, S4). In addition, to reduce the model overfitting problem, a fixed bandwidth regime rather than adaptive regime was applied to derive (SI, S4). With the GTWR model for the entire country using a general bandwidth of 120 km·day, we, however, could not predict robust PM2.5 concentrations extensively in northwestern China for some days due to the limited number of matched samples (SI, S6). To stabilize the model performance in areas where few matched samples were available, two Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 75

separate GTWR models were fitted for eastern (E model) and western (W model) China with the same model structure as Eq. (1). This approach decreased the effect of the lower number of samples for western China on model performance without degrading overall accuracy (SI, S6). Therefore, fixed bandwidths of 100 km ∙ day for eastern China and 290 km ∙ day for the west were applied in the final model. Western China was defined as the provinces of Xinjiang, Qinghai, Tibet, and the western parts of Inner Mongolia, Gansu, and Sichuan; eastern China contained the remaining areas (SI, Fig. S4). A scale factor ρ of 0.2 with the best model performance was obtained. A 10-fold cross-validation (CV) technique was used to examine model overfitting (Rodriguez et al., 2010). To evaluate model performance, several statistical indicators—R square (R2 ), root mean square error (RMSE), and mean prediction error (MPE)—were calculated by comparing differences between the estimated PM2.5 concentrations and surface PM2.5 measurements. The Pearson correlation coefficient was used in this study to quantify the correlation (r). 3.2. Daily GWR model and two-stage model To evaluate the performance of GTWR based on its integrated spatiotemporal weighting scheme, we also fitted daily GWR (D-GWR) models and two-stage models with the structure as in S7 (SI). The DGWR models separately calibrate the PM2.5-AOD relationship for individual days, and could thus be regarded as a special case of GTWR if the number of estimation days equals one. The two-stage model (Hu et al., 2014) performs an LME analysis in the first stage and a GWR analysis in the second stage. We first attempted two separate D-GWR (East and West) models using the same partitions as in the Seg-GTWR model, and the results show that the West D-GWR model failed to perform stably due to the insufficiency of samples in western China. For the two-stage model, the results of the West model show a large overfitting problem, with R2 decreasing from 0.82 in model fitting to 0.66 in model validation. Therefore, we only developed East D-GWR and two-stage models for comparison. 4. Results 4.1. Descriptive statistics Table 1 shows the descriptive statistics of the grid-averaged PM2.5 concentrations in 60,810 grid cells and grid-collocated independent variables distributed over the 325 days used in fitting the models. The mean PM2.5 concentration averaged for the entire study region was 60 μg/m3 with a standard deviation (Std. Dev) of 41 μg/m3 , and the mean combined AOD was 0.50 with a Std. Dev of 0.35. The mean PM2.5 values did not differ greatly between the two areas, but the mean aerosol loading was much higher in eastern versus western China (0.50 vs. 0.34). The values of the dependent and independent variables presented different seasonal variability. The maximum mean PM2.5 occurred in winter (86 μg/m3 ), and the minimum was reached in summer (44 μg/m3 ). The mean combined AOD values did not differ much among the four seasons, with the highest in summer (0.56) and the lowest in autumn (0.43). The detailed statistics for each season are listed in Table S7 (SI). The geographical distribution of the annual and seasonal mean PM2.5 measurements and combined AOD is shown in Fig. S5 (SI). Spatial ground PM2.5 exhibits a descending gradient from northeast to southwest due to differences in anthropogenic activities, land use, topography, and meteorological conditions. Generally, high levels of PM2.5 concentrations and AOD values occur in densely populated areas with thriving industries and low elevations (e.g., the North China Plain, Sichuan Basin, and Central China) and in the arid and semi-arid areas of northwestern China (e.g., the Taklimakan Desert). 4.2. Performance evaluation of the GTWR model Fig. 2 shows the scatter plots of GTWR model fitting and CV results for eastern China, western China, and the entirety of mainland China (the Seg model represents the E + W model, and the non-Seg model is the GTWR model using all matched samples without segmentation). For the model fitting, the respective overall R2 and RMSE values between the estimated and observed PM2.5 are 0.85 and 15.78 μg/m3 based on the results of the E + W model. For model validation, CV R2 and RMSE are 0.80 and 18.58 μg/m3 , which is a slight improvement compared to the R2 and RMSE of 0.76 and 19.77 μg/m3 of the non-Seg model. The CV results suggest that the Seg model is slightly overfitted but the nonSeg model is significantly overfitted, i.e., CV R2 decreases by 0.05 and 0.11 and RMSE increases by 2.80 and 5.21 μg/m3 for the two models, respectively. It is likely that the model itself could be improved by dividing mainland China into two parts, because the Seg model had comparable performance in model fitting, higher values in the CV processes, and much less overfitting than the non-Seg model. Therefore, we chose to fit two separate models to predict PM2.5 concentrations over China. The seasonal results of the Seg-GTWR model in fitting and CV processes based on the E and W models are shown in Fig. 3. It is clear that the model validation performance shows a similar seasonal pattern to that of model fitting. The values of R2 reach their maximum in autumn with 0.85 for model fitting and 0.80 for CV, followed by winter with respective values of 0.84 and 0.77. Correspondingly, a lower R2 for Table 1 Descriptive statistics of the grid-collocated data for all variables used in the GTWR model. Region Unit PM2.5 AOD NDVI RH WS T PBL DEM μg/m3 – – % m/s °C m m Whole (N = 60,810) Mean 60 0.50 0.35 57.01 2.22 15.32 650.99 370.49 Std. Dev 41 0.35 0.15 15.06 0.96 8.82 302.16 516.87 Min 1 0.00 0.00 9.39 0.15 −13.90 20.39 −13.47 Max 673 3.89 0.90 97.00 13.64 34.12 2137.74 4553.23 Median 50 0.43 0.34 59.01 2.03 16.94 602.01 97.49 East (N = 58,707) Mean 60 0.50 0.35 57.58 2.23 15.30 646.02 346.66 Std. Dev 40 0.35 0.15 14.85 0.96 8.84 296.68 492.17 Min 1 0.00 0.00 9.45 0.15 −13.90 20.39 4.09 Max 595 3.89 0.90 97.00 13.64 33.69 2137.74 3323.93 Median 51 0.43 0.34 59.65 2.03 16.95 598.48 87.58 West (N = 2103) Mean 57 0.34 0.29 40.93 1.97 16.03 789.72 1035.73 Std. Dev 62 0.32 0.11 11.80 0.70 8.22 403.83 710.44 Min 2 0.00 0.08 9.39 0.54 −8.44 75.84 −13.47 Max 673 2.94 0.67 76.89 6.66 34.12 2095.89 4553.23 Median 37 0.25 0.29 39.60 1.87 16.69 761.45 934.51 Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 76

model fitting and validation are observed in warm seasons (R2 and CV R2 are respectively 0.83 and 0.75 in spring and 0.80 and 0.71 in summer). However, contrary to R2 , summer is the season with the minimum error, with values of 12.21 and 15.61 μg/m3 for model fitting and validation and 7.54 and 10.18 μg/m3 for the RMSE and MPE, respectively. Winter is the season with the highest RMSE values of 20.36 and 24.82 μg/m3 (for model fitting and validation) and MPE values of 13.83 and 16.83 μg/m3 . The much higher RMSE and MPE values in autumn and winter are mainly attributable to the inclusion of more high PM2.5 values (> 100 μg/m3 ) during the cold seasons. However, model performance does not differ significantly over the four seasons (i.e., 95% of the daily local R2 values are > 0.75. The annual mean local R2 values for each grid vary from 0.33 to 0.99, and more than about 91% of the values are > 0.75. These results indicate that the performance of the daily local models is relatively high. In particular, the spatial distribution of the mean local R2 shown in Fig. 4 reveals that the values of local R2 tend to be higher in the densely populated areas of southeastern China with higher PM2.5 concentrations than in the less populated areas of western China with lower PM2.5 concentrations, which is partly attributable to the uneven distribution of monitoring stations. These spatial and temporal variations in the mean local R2 reveal a relatively stable performance for the fitted Seg-GTWR model, with some spatiotemporal variability for daily local regressions, and demonstrate its ability to resolve the large spatiotemporal heterogeneities of the PM2.5- AOD relationship. 4.3. Performance comparison between the models To comprehensively understand the difference between the GTWR, D-GWR, and two-stage models, D-GWR and two-stage models based on the matched sample dataset from eastern China were implemented. It was found that only 317 daily GWR models were fitted and models could not be constructed for 8 days due to the limited PM2.5-AOD pairs (daily samples < 6) on these days, because D-GWR cannot capitalize on temporal dependency in the data and thus make use of samples before the estimation days. However, the two-stage and GTWR models could generate PM2.5 estimations for all 325 days by taking advantage of both spatial and temporal dependency in the data. For this reason, the two-stage and GTWR models are much more robust in building the PM2.5-AOD relationship, especially for days without sufficient matched samples. It should be noted that the remaining 40 days could not be estimated by the two-stage and GTWR models mainly because of the unavailability of AOD retrievals due to, e.g., cloud contamination (see Section 2.2). y=0.81x+11.36 R2 =0.85 RMSE=15.28 MPE=10.16 N=58707 y=0.81x+11.69 R2 =0.80 RMSE=18.00 MPE=12.03 N=58707 (a) (e) y=0.81x+11.88 R2 =0.80 RMSE=18.58 MPE=12.27 N=60810 y=0.77x+12.71 R2 =0.82 RMSE=26.07 MPE=15.09 N=2103 y=0.81x+11.50 R2 =0.85 RMSE=15.78 MPE=10.33 N=60810 y=0.77x+14.31 R2 =0.72 RMSE=30.73 MPE=18.79 N=2103 y=0.83x+9.86 R2 =0.87 RMSE=14.56 MPE=9.53 N=60810 y=0.79x+13.02 R2 =0.76 RMSE=19.77 MPE=12.61 N=60810 (b) (c) (d) (f) (g) (h) Counts of points Fig. 2. Scatter plots for model fitting and cross validation. (a)–(d) are the model fitting results of the E, W, Seg-GTWR (E + W), and non-Seg GTWR models, respectively. (e)–(h) are the CV results of the E, W, Seg-GTWR, and non-Seg GTWR models, respectively. The blue line is the linear regression of the scattered dots; the dotted line is the 1-1 line as a reference. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 77

To compare the overall performance of the GTWR, D-GWR, and two-stage models in detail, the statistics of the results from the three East models are listed in Table 2. The D-GWR model performs the worst, with a significant overfitting problem and the lowest CV R2 of the three models. For the two-stage model with consideration of day-today variations in the PM2.5-AOD relationship, the first-stage LME model with an R2 of 0.63 shows a notable improvement from the estimates generated by a multivariate linear regression model (R2 of 0.34 in S3 (SI)). Adding GWR into the modeling process to account for spatial variability in the PM2.5-AOD relationship improves the performance of the two-stage model (R2 of 0.75 and CV R2 of 0.72). However, despite a combined modeling process, the two-stage model does not capture the spatiotemporal heterogeneities simultaneously within the PM2.5-AOD relationship as does the GTWR model. Through the synergy of spatial and temporal dimensions within the data, GTWR can calibrate the spatiotemporally heterogeneous PM2.5-AOD relationship on a local basis. As a result, the GTWR model, with a higher R2 of 0.85 and CV R2 of 0.80, outperforms the two-stage (LME + GWR) model, although the difference between the model fitting and validation is slightly greater than that generated by the two-stage model. These findings suggest that, despite a slight overfitting problem, the overall estimating performance of GTWR is significantly improved when spatial and temporal non-stationarities are simultaneously considered. 4.4. Predicted PM2.5 concentrations over China 4.4.1. General spatiotemporal variation of PM2.5 in China Surface PM2.5 maps for mainland China with a 3-km grid were derived using the Seg-GTWR model. The overall annual mean of predicted PM2.5 concentrations for the entirety of China is 44 μg/m3 , which is y=0.78x+12.01 R2 =0.83 RMSE=13.32 MPE=9.22 N=19426 y=0.75x+10.88 R2 =0.80 RMSE=12.21 MPE=7.54 N=11298 y=0.78x+12.65 R2 =0.75 RMSE=15.72 MPE=10.77 N=19426 y=0.72x+12.39 R2 =0.71 RMSE=15.61 MPE=10.18 N=11298 y=0.78x+12.82 R2 =0.85 RMSE=16.46 MPE=10.80 N=17187 y=0.80x+17.53 R2 =0.84 RMSE=20.36 MPE=13.83 N=12899 y=0.79x+12.26 R2 =0.80 RMSE=17.88 MPE=11.90 N=17187 y=0.79x+18.54 R2 =0.77 RMSE=24.82 MPE=16.83 N=12899 (e) (f) (g) (h) (a) (b) (c) (d) Counts of points Fig. 3. Scatter plots of Seg-GTWR's model fitting and cross validation. (a)–(d) are the model fitting results for spring, summer, autumn, and winter, respectively. (e)–(h) are the CV results for spring, summer, autumn, and winter. (a) (b) Local R2 Fig. 4. (a) Spatial distribution of the annual mean and (b) temporal variation of spatially averaged local R2 of the Seg-GTWR model for each grid with PM2.5 monitors. Note that we use a circle rather than a grid cell shape in the left panel because the 3 km grid located in the figure cannot display the variation of color. Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 78

slightly less than the annual mean value of observed PM2.5 (60 μg/m3 ) because the annual mean of PM2.5 estimates is calculated based on urban, rural, and depopulated zones, but the latter only represents the average level of urban areas. The difference between the annual predicted and measured PM2.5 may also result from the difference in the days with available AOD and PM2.5 measurements (due to cloud contamination, the AOD images on many days do not have enough retrievals for some areas). Fig. 5 shows the predicted maps of seasonal and annual mean PM2.5 concentrations. For comparison, the groundmeasured PM2.5 concentrations for areas in which the gridded AOD is available are also shown in Fig. 5, and the spatial and temporal patterns of both satellite-derived and ground-observed PM2.5 concentrations are in accordance. Fig. S6 (SI) also shows good agreement between the satellite-derived and measured PM2.5 for all days, because the differences for ~80% of monitoring stations fall within ± 10 μg/m3 . These findings confirm that GTWR is a promising modeling method for exploring population exposure to PM2.5. Spatially, the general distribution of annual PM2.5 concentrations demonstrates that high PM2.5 usually occurred in the economically and industrially developed areas of eastern and southern China, whereas low PM2.5 occurred in the rural and less-developed areas of western and northern China, corresponding well with the general distribution of the Chinese population (i.e., dense in southeastern China and sparse in northwestern China) (Hu, 1935) and inversely with the variation of terrain in China (i.e., low in the southeast and high in the northwest; see Fig. 1). Temporally, there were still remarkable seasonal differences in concentrations, with the most polluted season during 2015 being winter (seasonal mean value = 67 μg/m3 ) and the cleanest being summer (35 μg/m3 ), corresponding to the seasonal variation revealed by Ma et al. (2014) and Zhang and Cao (2015). Winter heating and unfavorable weather for fine particle dispersion contribute substantially to the wintertime maximum (Zhang and Cao, 2015). It is clear that over 70% of China, corresponding to > 1.25 billion or 92% of the Chinese population, comprises areas with risky levels of PM2.5 pollution according to the CNAAQS Level 2 standard (an annual mean value of PM2.5 concentration exceeding 35 μg/m3 ). Over 99% of China suffers from particle pollution exceeding the CNAAQS Level 1 standard of 15 μg/m3 . 4.4.2. PM2.5 hotspots Particularly heavily polluted areas are clustered in several centers (Fig. 6). The North China Plain, including Beijing, Tianjin, Hebei, Table 2 Comparisons of the model performance of D-GWR, two-stage and GTWR. Model Model fitting Model validation Days can be estimated R2 RMSE MPE R2 RMSE MPE GTWR 0.85 15.28 10.16 0.80 18.00 12.03 317 Fit line: y = 0.81x + 11.36 Fit line: y = 0.81x + 11.69 D-GWR 0.83 16.63 11.12 0.71 21.21 13.50 325 Fit line: y = 0.80x + 11.94 Fit line: y = 0.79x + 12.92 Two-stage Stage 1 0.63 24.22 16.91 0.61 24.72 17.35 325 Fit line: y = 0.62x + 22.62 Fit line: y = 0.61x + 23.16 Stage 2 0.75 19.87 13.91 0.72 21.01 14.59 Fit line: y = 0.74x + 15.97 Fit line: y = 0.73x + 16.67 (e) Autumn (f) Autumn (g) Winter (h) Winter (i) Annual (j) Annual (a) Spring (b) Spring (c) Summer (d) Summer Fig. 5. Seasonal and annual mean predicted PM2.5 and surface observed PM2.5 concentrations where gridded AOD is available. Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 79

Henan, Anhui, Shandong, and Jiangsu provinces, is one of the most severely particle-polluted areas in China. The annual mean PM2.5 estimate is 62 μg/m3 , and mean values for various subareas generally range from 20 to 104 μg/m3 . The heavy pollution is mainly attributable to activities related to rapid economic development and accelerated urbanization and is aggravated by the unfavorable topography of the Taihang mountains located in the west of Hebei province, which prevent fine particle dispersion (Fu et al., 2014; Quan et al., 2011; Tao et al., 2012). The Taklimakan Desert, located in the southern Xinjiang Autonomous Region, is another heavily polluted region. The annual mean PM2.5 level of this desert (the largest in China) is generally 62 μg/ m3 , with a maximum of 88 μg/m3 and minimum of 38 μg/m3 , due to the large amount of airborne dust and sand originating from the desert surface. Central China, including Hubei and Hunan provinces, is another hotspot with high pollution, with an annual mean predicted PM2.5 of 55 μg/m3 , which is related to low elevation, abundant bodies of water, thriving industries, and intense human activities (e.g., fossil fuel combustion) (Guo et al., 2012; Wang et al., 2015). The Sichuan Basin, which contains Chongqing and eastern Sichuan province, often suffers from episodes of heavy air pollution with an annual mean value of 49 μg/m3 , which are related to stagnant air circulation and high mountains at its margins (Li et al., 2015). Influenced by dust storms from the Taklimakan Desert (Han et al., 2009), the northern Tibetan Plateau has annual mean values of 25–55 μg/m3 , much higher than those in the southern part. The areas with low PM2.5 concentrations are located mainly in less densely populated areas, including northern Xinjiang, western Sichuan, southern Tibet, eastern Inner Mongolia, Hainan, and Fujian, with PM2.5 values generally lower than 35 μg/m3 . According to Ma et al. (2014), the general spatial pattern of PM2.5 pollution in China is similar, except that Central China is becoming a PM2.5 hotspot, and the northern Tibetan Plateau has a relatively high level of PM2.5 compared with 2013 levels. 5. Discussion 5.1. Predictive power of the GTWR model Comparisons between the GTWR, D-GWR, and two-stage models show that GTWR outperforms the latter two models in expressing the PM2.5-AOD relationship for days with PM2.5-AOD paired samples. Given that GTWR involves simultaneous consideration of spatiotemporal effects within the data, it is reasonable to assume that this model is capable of predicting PM2.5 concentrations for some days without samples. Therefore, we carried out two separate 10-fold CVs that randomly dropped 10% of monitoring stations (space CV) and days (time CV) for East China to examine the predictive power of GTWR. Compared to the model fitting results in Fig. 2(a), the CV results in Fig. 7(a–b) suggest that some overfitting exists when the CV technique is changed. However, because the distributions of PM2.5-AOD paired samples differ from day to day, the space and time CV not only examine overfitting in space and time separately, but also reveal the predictive power of GTWR modeling, especially for time CV. Despite the presence of overfitting in time CV, GTWR can explain 58% of variation in daily ground PM2.5 for days without samples (Fig. 7(b)), which indicates that the GTWR model has a better predictive ability than other statistical models (i.e., generating relatively reliable PM2.5 predictions for days without PM2.5 observations). Furthermore, we collected the required data in 2014 and processed them using the methods presented in Section 2 to examine GTWR's ability to predict historical PM2.5 concentrations. The relationships built in 2015 were used to predict PM2.5 in 2014 with the same assumption as in Ma et al. (2016) that the daily PM2.5-AOD relationship is PM2.5 (μg/m3 ) (a) (b) (c) Fig. 6. Zoom-in map of the predicted annual mean PM2.5 for pollution hotspots: (a) the North China Plain, (b) Sichuan Basin, and (c) Central China. Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 80

consistent for the same day of the preceding year. The results in Fig. 7(c) show an improved accuracy at higher spatial resolution as compared to the study of Ma et al. (2016) at a coarse scale (with CV R2 of 0.41 for the daily mean of historical PM2.5 predictions in China). These findings reveal that GTWR can extend the temporal coverage of PM2.5 observations and holds great potential in predicting historical PM2.5 concentrations. In addition, the model performance of GTWR is affected by the integration method of spatial and temporal distances. Fig. S7 (SI) shows that the CV R2 values vary with the scale factor ρ. Therefore, in addition to the widely examined parameter of bandwidth in GWR, the scale factor related to the metric system of spatiotemporal distance is also a key parameter for GTWR. 5.2. Effect of combined AOD The coverage of our combined AOD data is improved (Fig. S8, SI). A comparison of Fig. S8(a–c) (SI) shows that merging 3 km DT AOD by averaging Aqua and Terra AOD values (including the retrievals and predicted values) during the whole year (Fig. S8(c), SI) increased data availability temporally by 10%–50% over the original 3-km Aqua/Terra data (Fig. S8(a–b), SI), and the improvements were mainly located over eastern areas. Correspondingly, the matched PM2.5-AOD samples increased from 2917 (for the AOD dataset that only fused Aqua and Terra aerosol retrievals) to 31,118 (for the merged 3 km AOD dataset that fused the Aqua and Terra retrieved and predicted values). Subsequently, daily coverage of the combined AOD further improved after the combination step (Step 4 in Section 2.2.3) of consolidating the merged 3 km DT AOD with the resampled-DB AOD (Fig. S8(d), SI). Comparison with Fig. S8(c) (SI) shows that the sampling frequency for densely populated eastern and southern China increased by 10%–80% (mean value ~ 35%). The data availability for each pixel over the arid and semi-arid areas of northwestern China showed the most improvement, from null value to ~150 valid values during the whole year; the Tibetan plateau also showed significant improvement (20%–300%). Lastly, the samples for GTWR modeling increased to 60,810. According to our statistics, the model using the dataset of 2917 samples showed poor performance (R2 = 0.59, CV R2 = 0.32), indicating that it is impossible to estimate PM2.5 for the entirety of China using the original MODIS AOD retrievals without coverage improvement. However, the model performance (R2 = 0.77, CV R2 = 0.60) showed an improvement using the dataset of 31,118 samples compared with the results of the model using the dataset of 2917 samples, but still performed much worse than the GTWR model using the combined AOD dataset. These findings indicate that it is necessary to improve the data availability of satellite-derived AOD, because model performance can be substantially improved with greater inclusion of PM2.5-AOD pairs for GTWR modeling. As a result, the combined AOD not only presents better spatial details of PM2.5 (compare Fig. S9(a) and (b)) and makes it possible to generate high-resolution PM2.5 for the entirety of China, but also enhances the model performance of GTWR. 5.3. Comparisons with previous studies In previous studies, satellite-based ground PM2.5 estimates for the entirety of China could only reveal coarse gradients, with grid sizes of > 10 km (Fang et al., 2016; Li et al., 2017; Ma et al., 2014; Ma et al., 2016; You et al., 2016). Using the newly released MODIS 3-km AOD data, Li et al. (2017) and You et al. (2016) generated ground PM2.5 predictions at a coarse resolution (≥10 km). Our study has shown how the combined AOD dataset could improve spatial coverage for each day and maintain the spatial details of the MODIS 3 km AOD, directly benefiting the generation of high-resolution PM2.5 predictions on a daily basis. By considering the estimating power of modeling daily PM2.5-AOD relationships, our Seg-GTWR model accounts for 80% of the daily variations in the 10-fold model validation process for high-resolution PM2.5 modeling, surpassing the performance of previous models, such as the daily GWR models developed for coarse-resolution modeling over the entirety of China by Ma et al. (2014) with a CV R2 of 0.64 and You et al. (2016) with a CV R2 of 0.79, and the mixed effects model (CV R2 = 0.79) (Xie et al., 2015) and two-stage spatiotemporal model (CV R2 = 0.72) (Wu et al., 2016) developed for fine-resolution PM2.5 mapping over small regions. In recent years, Bai et al. (2016) and Guo et al. (2017) used the GTWR model to estimate PM2.5 concentrations, but restricted their study areas to small regions and did not tap the possibility of optimizing GTWR parameters or examine its estimating/ predictive power in elucidating the PM2.5-AOD relationship for days without samples or in historical years. In addition, the GTWR model developed in the present study also achieved a performance comparable to those achieved in the U.S. using MODIS 1 km AOD retrieved with the MAIAC algorithm (CV R2 = 0.62–0.84) (Chudnovsky et al., 2014 and Hu et al., 2014). 6. Conclusion Despite the increasing number of studies of satellite-derived PM2.5 for China, there is still little detailed information on the spatial pattern of fine particles due to the limitations of the spatial resolution of satellite AOD. In using combined high-resolution MODIS AOD data, this y=0.61x+25.13 R2 =0.58 RMSE=28.24 MPE=19.31 N=58707 y=0.79x+12.38 R2 =0.75 RMSE=20.73 MPE=13.31 N=58707 Counts of points (a) (b) y=0.54x+25.84 R2 =0.47 RMSE=37.57 MPE=24.51 N=25601 (c) Fig. 7. Scatter plots of model CV results for eastern China: (a) 10-fold CV results with 10% of monitoring sites dropped; (b) 10-fold CV results with 10% of days dropped; (c) evaluation of historical PM2.5 estimations in 2014 on a daily basis. Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 81