A Geographically and Temporally Weighted Regression Model for Ground-Level PM 2 . 5 Estimation from Satellite-Derived 500 m Resolution AOD

Regional haze episodes have occurred frequently in eastern China over the past decades. As a critical indicator to evaluate air quality, the mass concentration of ambient fine particulate matters smaller than 2.5 μm in aerodynamic diameter (PM2.5) is involved in many studies. To overcome the limitations of ground measurements on PM2.5 concentration, which is featured in disperse representation and coarse coverage, many statistical models were developed to depict the relationship between ground-level PM2.5 and satellite-derived aerosol optical depth (AOD). However, the current satellite-derived AOD products and statistical models on PM2.5–AOD are insufficient to investigate PM2.5 characteristics at the urban scale, in that spatial resolution is crucial to identify the relationship between PM2.5 and anthropogenic activities. This paper presents a geographically and temporally weighted regression (GTWR) model to generate ground-level PM2.5 concentrations from satellite-derived 500 m AOD. The GTWR model incorporates the SARA (simplified high resolution MODIS aerosol retrieval algorithm) AOD product with meteorological variables, including planetary boundary layer height (PBLH), relative humidity (RH), wind speed (WS), and temperature (TEMP) extracted from WRF (weather research and forecasting) assimilation to depict the spatio-temporal dynamics in the PM2.5–AOD relationship. The estimated ground-level PM2.5 concentration has 500 m resolution at the MODIS satellite’s overpass moments twice a day, which can be used for air quality monitoring and haze tracking at the urban and regional scale. To test the performance of the GTWR model, a case study was carried out in a region covering the adjacent parts of Jiangsu, Shandong, Henan, and Anhui provinces in central China. A cross validation was done to evaluate the performance of the GTWR model. Compared with OLS, GWR, and TWR models, the GTWR model obtained the highest value of coefficient of determination (R2) and the lowest values of mean absolute difference (MAD), root mean square error (RMSE), and mean absolute percentage error (MAPE).


Introduction
Epidemiologic studies demonstrated that fine particulate matters smaller than 2.5 µm in aerodynamic diameter (PM 2.5 ), associated with increased cardiopulmonary and cardiovascular mortality, have a significant negative impact on human health [1][2][3][4][5].The measurement of pollutant concentrations including ground-level PM 2.5 is therefore performed regularly in many industrialized countries that have established air quality monitoring stations a short distance from residential or manufacturing districts [6][7][8][9].Traditional air quality estimations from ground-based stationary ambient variability comparison between satellite-derived AOD and PM 2.5 , Li et al. [36] found diverse, even opposite, seasonal cycles in the PM 2.5 -AOD relationship owing to the varied height of the atmospheric mixing layer.Such studies highlighted the importance of spatio-temporal variability in PM 2.5 -AOD relations, but failed to explicitly and simultaneously incorporate variable factors into models to infer ground-level PM 2.5 concentrations from satellite AOD estimations.Referring to the shortage of ground-level PM 2.5 concentration of high spatio-temporal revolution for air quality monitoring and population exposure studies, we aim to develop a flexible estimation model to calibrate the satellite-derived AOD at 500 m spatial revolution to hourly ground-level PM 2.5 concentrations.
In this article, a geographically and temporally weighted regression (GTWR) model based on a satellite-derived AOD is proposed to estimate hourly ground-level PM 2.5 concentration images with fine spatial resolution at the regional scale.First, a simplified high-resolution MODIS aerosol retrieval algorithm (SARA) is selected to produce satellite-derived AOD images with 500 m spatial resolution.Second, the weather research and forecasting (WRF) model is employed to extract the hourly meteorological field data.Third, we examine and compare the quality of PM 2.5 estimations from GTWR to ordinary least squares (OLS), GWR, and temporally weighted regression (TWR) methods by means of case study.The case study region, named JSHA linking region, is a famous fossil energy base that covers the adjacent parts of Jiangsu, Shandong, Henan, and Anhui provinces in central China.

Study Region
With rapid economic development and urbanization, the consumption of fossil fuels including coal and petroleum has increased dramatically in China, resulting in severe regional air pollution and frequent haze episodes.The JSHA linking region takes Xuzhou as its center, which extends about 150 km to the north, west, and south from Xuzhou, and about 50 km to the east from Xuzhou (see Figure 1).The study region includes a large number of coal mines and thermal power plants, and numerous cities with a population of hundred thousands to several million.This particular economic region suffers frequently from haze episodes, especially in the winter heating period, due to coal burning in heating stations and power plants.
Remote Sens. 2016, 8, 262 3 of 21 variability comparison between satellite-derived AOD and PM2.5, Li et al. [36] found diverse, even opposite, seasonal cycles in the PM2.5-AOD relationship owing to the varied height of the atmospheric mixing layer.Such studies highlighted the importance of spatio-temporal variability in PM2.5-AOD relations, but failed to explicitly and simultaneously incorporate variable factors into models to infer ground-level PM2.5 concentrations from satellite AOD estimations.Referring to the shortage of ground-level PM2.5 concentration of high spatio-temporal revolution for air quality monitoring and population exposure studies, we aim to develop a flexible estimation model to calibrate the satellite-derived AOD at 500 m spatial revolution to hourly ground-level PM2.5 concentrations.
In this article, a geographically and temporally weighted regression (GTWR) model based on a satellite-derived AOD is proposed to estimate hourly ground-level PM2.5 concentration images with fine spatial resolution at the regional scale.First, a simplified high-resolution MODIS aerosol retrieval algorithm (SARA) is selected to produce satellite-derived AOD images with 500 m spatial resolution.Second, the weather research and forecasting (WRF) model is employed to extract the hourly meteorological field data.Third, we examine and compare the quality of PM2.5 estimations from GTWR to ordinary least squares (OLS), GWR, and temporally weighted regression (TWR) methods by means of case study.The case study region, named JSHA linking region, is a famous fossil energy base that covers the adjacent parts of Jiangsu, Shandong, Henan, and Anhui provinces in central China.

Study Region
With rapid economic development and urbanization, the consumption of fossil fuels including coal and petroleum has increased dramatically in China, resulting in severe regional air pollution and frequent haze episodes.The JSHA linking region takes Xuzhou as its center, which extends about 150 km to the north, west, and south from Xuzhou, and about 50 km to the east from Xuzhou (see Figure 1).The study region includes a large number of coal mines and thermal power plants, and numerous cities with a population of hundred thousands to several million.This particular economic region suffers frequently from haze episodes, especially in the winter heating period, due to coal burning in heating stations and power plants.

Satellite-Derived AOD
With the evolution of satellite remote sensing sensors and retrieval algorithms, more and more satellite-derived AOD products have been employed to estimate the ground-level PM2.5 concentration.

Satellite-Derived AOD
With the evolution of satellite remote sensing sensors and retrieval algorithms, more and more satellite-derived AOD products have been employed to estimate the ground-level PM 2.5 concentration.The MODIS AOD product with high temporal resolution has been widely used, twice a day, for environmental monitoring in numerous studies.However, its coarse spatial resolution of 10 km ˆ10 km does not allow detailed regional-or urban-scale air quality monitoring.In this article, the SARA algorithm developed by Bilal et al. [37], which shows higher accuracy and better reliability in AOD retrieval than the MODIS 04 algorithm over the complex surfaces under low and high aerosol loadings with mixed aerosol types [38], is used to generate high-resolution MODIS AOD images.SARA is applied with the Radiative Transfer Model (RTM) to directly calculate MODIS AOD at 500 m resolution based on MODIS (Terra/AQUA) level 1 (MOD02HKM calibrated radiance and MOD03 geolocation data) and level 2 (MOD09 surface reflectance) products with aerosol properties including single scattering albedo (SSA) and asymmetry parameter obtained from the Aerosol Robotic Network (AERONET) station at China University of Mining and Technology (http://aeronet.gsfc.nasa.gov/).The temporal averaging on the AERONET product is done within ˘30 min of the MODIS local overpass moment.

Ground-Level PM 2.5 Measurements
For the sake of monitoring and supervising the air quality in China, the Ministry of Environmental Protection of the People Republic of China has established more than 1497 monitoring stations over 367 cities and collects the measurements of PM 10 , PM 2.5 , SO 2 , NO 2 , O 3 , and CO in real time.In this study, hourly averaged PM 2.5 concentrations when the satellite passes over the study region were downloaded from the national air quality publishing platform (http://106.37.208.233:20035/).The locations of Xuzhou-CUMT AERONET station and 37 air quality monitoring stations in the study region are shown in Figure 1.

Auxiliary Data
In order to improve the performance of the PM 2.5 -AOD model, auxiliary data including PBLH, relative humidity (RH), wind speed (WS), and temperature (TEMP) are used as meteorological field inputs for air quality modeling.However, it is difficult to extract all the parameters from current remote sensing sensors while the ground-observed meteorological data are space limited.The WRF model v3.4.1 [39][40][41] is therefore considered to extract the needed meteorological field information based on the initial and boundary conditions generated by WRF Processing System (WPS) using NCEP FNL Operational Model Global Tropospheric Analyses dataset of 1 ˝ˆ1 ˝resolution (http://rda.ucar.edu/dsszone/ds083.2/).The WRF model is a mesoscale numerical weather prediction system designed for both atmospheric research and operational forecast, and serves as a wide range of meteorological applications across scales from tens of meters to thousands of kilometers.The nested domain scheme with 3 km horizontal grid space of WRF output is adopted in this paper, and the temporal resolution of WRF outputs was set at 1-h intervals.The physical options used in WRF include the single-moment 6-class (WSM6) microphysics, the Yonsei University (YSU) PBL scheme, the Rapid Radiative Transfer Model (RRTM) longwave and Dudhia shortwave radiation schemes, and Noah land surface model.Further, the simulated meteorological fields of 3 km ˆ3 km resolution are interpolated to 500 m ˆ500 m resolution, the same as SARA AOD products, based on cubic spline interpolation algorithm in Matlab.

Descriptive Statistics
The data from the above three sources, co-located in both space and time, were selected to establish a dataset for modeling during the period of 1 November 2014 to 28 February 2015.That is to say, only the data that were spatially co-located at the selected 37 ground stations and temporally matched at the MODIS overpasses were included in the following analyses.Finally the number of the dataset is 1105, which is about 1/8 of the observations under ideal observation conditions.
Descriptive statistics of variables in the proposed model are summarized in Figure 2. The SARA AOD data used in the model have a mean value of 0.479 and a standard deviation of 0.228.
The corresponding hourly ground-level PM 2.5 concentration ranges from 7.0 µg/m 3 to 405.0 µg/m 3 , with a mean value of 84.35 µg/m 3 and a standard deviation of 59.73 µg/m 3 (Figure 2a,b).This shows that SARA AOD and ground-level PM 2.5 exhibit a similar frequency distribution on the lower bound of their value ranges.For the WRF assimilated meteorological fields, PBLH (592 m ˘261 m), RH (31% ˘8%), and WS (4.4 m/s ˘2 m/s) are found to be less variable but show a nearly normal distribution (Figure 2c-e).WRF TEMP varies dramatically within the range of ´2 ˝C to 17 ˝C in our dataset, which also exhibits a nearly normal distribution, with a mean value of 6 ˝C and a standard deviation of 4.2 ˝C (Figure 2f).
For further understanding the response of ground level PM 2.5 concentrations to meteorological condition changes, the LOESS (locally weighted polynomial regression) trend lines were added in the scatter-plots between ground level PM 2.5 and the predictors in SPSS (Figure 3).LOESS combines the simplicity of linear least squares regression with the flexibility of nonlinear regression and fits simple models to localized subsets of the data [42].It shows that SARA AOD and RH seems to be significantly associated with ground level PM 2.5 concentration in a nearly positive linear model (Figure 3a,c), while PBLH and WS fit with a nearly negative linear model (Figure 3b,d).Ground-level PM 2.5 concentration increases when the TEMP is below around 7 ˝C, but decreases as the temperature gets higher than 7 ˝C (Figure 3e).AOD and ground-level PM2.5 exhibit a similar frequency distribution on the lower bound of their value ranges.For the WRF assimilated meteorological fields, PBLH (592 m ± 261 m), RH (31% ± 8%), and WS (4.4 m/s ± 2 m/s) are found to be less variable but show a nearly normal distribution (Figure 2c-e).WRF TEMP varies dramatically within the range of −2 °C to 17 °C in our dataset, which also exhibits a nearly normal distribution, with a mean value of 6 °C and a standard deviation of 4.2 °C (Figure 2f).For further understanding the response of ground level PM2.5 concentrations to meteorological condition changes, the LOESS (locally weighted polynomial regression) trend lines were added in the scatter-plots between ground level PM2.5 and the predictors in SPSS (Figure 3).LOESS combines the simplicity of linear least squares regression with the flexibility of nonlinear regression and fits simple models to localized subsets of the data [42].It shows that SARA AOD and RH seems to be significantly associated with ground level PM2.5 concentration in a nearly positive linear model (Figure 3a,c), while PBLH and WS fit with a nearly negative linear model (Figure 3b,d).Ground-level PM2.5 concentration increases when the TEMP is below around 7 °C, but decreases as the temperature gets higher than 7 °C (Figure 3e).

Geographically and Temporally Weighted Regression (GTWR) Model
Compared with traditional linear or nonlinear regression methods, the GWR model, as one of the spatial varying coefficient regression models, can significantly improve the estimation accuracy of ground-level PM2.5, especially in a study region with a complex spatial context.Furthermore, Huang et al. [43] developed a GTWR model to deal with the spatial and temporal non-stationarity issues simultaneously through incorporating temporal effects into the GWR model.In this article, to effectively delineate the spatio-temporal quantitative relationship between satellite-derived AOD, coupled with meteorological data, and ground-level PM2.5 concentrations at the regional scale, a GTWR model was developed to improve the local coefficient of determination parameter (R 2 ) for satellite-derived AOD and ground PM2.5 measurements on hourly basis.
Supposing the sample is labeled as i , the GTWR model for the PM2.5-AOD relationship can be mathematically expressed as: ) where ( , , ) , , ) denote the slope of AOD, PBLH, RH, TEMP, and WS, respectively; and i ε stands for the random error.The meteorological variables are used as sensitive factors to improve the association between satellite-derived AOD and ground-level PM2.5 concentrations because the RH and TEMP can be employed for indicating the changes in particle compositions and the PBLH and WS influence the vertical profile of particulate matters.
In order to calibrate the parameters in the GTWR model, a local weighted least squares algorithm is usually employed, assuming that the closer the measurements to point i in space-time coordinate system, the greater the weight the measurements have in predicting ( , , ) . Thus, the estimation of parameters ( , , ) is expressed as:

Geographically and Temporally Weighted Regression (GTWR) Model
Compared with traditional linear or nonlinear regression methods, the GWR model, as one of the spatial varying coefficient regression models, can significantly improve the estimation accuracy of ground-level PM 2.5 , especially in a study region with a complex spatial context.Furthermore, Huang et al. [43] developed a GTWR model to deal with the spatial and temporal non-stationarity issues simultaneously through incorporating temporal effects into the GWR model.In this article, to effectively delineate the spatio-temporal quantitative relationship between satellite-derived AOD, coupled with meteorological data, and ground-level PM 2.5 concentrations at the regional scale, a GTWR model was developed to improve the local coefficient of determination parameter (R 2 ) for satellite-derived AOD and ground PM 2.5 measurements on hourly basis.
Supposing the sample is labeled as i, the GTWR model for the PM 2.5 -AOD relationship can be mathematically expressed as: where pu i , v i , t i q represents the given coordinates of the training sample i in space pu i , v i q at time t i ; and α 5 pu i , v i , t i q denote the slope of AOD, PBLH, RH, TEMP, and WS, respectively; and ε i stands for the random error.The meteorological variables are used as sensitive factors to improve the association between satellite-derived AOD and ground-level PM 2.5 concentrations because the RH and TEMP can be employed for indicating the changes in particle compositions and the PBLH and WS influence the vertical profile of particulate matters.
In order to calibrate the parameters in the GTWR model, a local weighted least squares algorithm is usually employed, assuming that the closer the measurements to point i in space-time coordinate system, the greater the weight the measurements have in predicting αpu i , v i , t i q.Thus, the estimation of parameters αpu i , v i , t i q is expressed as: where Wpu i , v i , t i q is a square matrix computing the geographical and temporal weighted values of training datasets for measurement i by the diagonal elements, and the parameter X is a vector representing the parameters AOD, PBLH, RH, TEMP, and WS.
Considering the different scale effects of location and time, an ellipsoidal coordinate system was used to measure the spatio-temporal distance between a regressive point and its surrounding measurement data.Through combining the temporal distance d T with the spatial distance d S , the spatio-temporal distance is defined as: where b denotes different kinds of operators.Moreover, Gaussian distance decay-based functions and Euclidean distance were chosen to construct the spatio-temporal weight matrix, and the "+" operator was adopted to assess the space-time distance d ST .As a result, the d ST and the diagonal element ω ij of weight matrix Wpu i , v i , t i q can be obtained as follows: where h is a non-negative parameter known as bandwidth, which produces a decay of influence with distance.
where ϕ S is the scale factor of spatial distance, while ϕ T stands for the scale factor of temporal distance. - where h ST , h T , and h S are the parameters of spatio-temporal, spatial, and temporal bandwidths, respectively.Therefore, if the bandwidth of spatial and temporal is decided, the weight matrix Wpu i , v i , t i q and αpu i , v i , t i q can be calculated.In practice, ph S , h T q can be obtained automatically with an optimization technique by cross-validation through minimizing Equation (7) in terms of R 2 statistics or the corrected Akaike information criterion (AIC) [44].
where the function ŷi phq denotes the predicted value y i from GTWR, and the sum of the squared error is written as in the function CVRSSphq.

Statistical Analysis
Some statistical indicators were employed to evaluate the estimation accuracy of GTWR model through comparing the fitted PM 2.5 concentration with ground PM 2.5 measurements, such as R 2 , root mean square error (RMSE), mean absolute difference (MAD), and mean absolute percentage error (MAPE).The R 2 of a statistical model is an indicator describing the discrepancy between ground-observed PM 2.5 and satellite-estimated ground-level PM 2.5 , where higher values indicate better fitting accuracy.The R 2 is defined as: where n is the number of ground measurement; PM obs 2.5 and PM obs 2.5 stand for the original and averaged concentration of ground-observed PM 2.5 , respectively; and PM sat 2.5 stands for the concentration of satellite-estimated ground-level PM 2.5 .
The RMSE is employed to measure the differences between ground-observed PM 2.5 and satellite estimated ground-lever PM 2.5 , which is sensitive to both systematic error and random error.The RMSE is calculated as: The MAD used to measure the magnitude of mean error magnitude is computed as: ˇˇPM obs 2.5 piq ´PM sat 2.5piq ˇˇ (10) The MAPE expressing estimation accuracy of a statistical model as a percentage is defined by:  4 presents a flowchart of the proposed estimating method.The Terra/Aqua MODIS level-1 and level-2 products with aerosol properties, including SSA and asymmetry parameter obtained from AERONET, are first employed to generate AOD images with 500 m spatial resolution at the satellite overpassing moment based on SARA.Then, meteorological field images including PBLH, RH, WS, and TEMP are extracted from WRF assimilation, and are geo-registered and interpolated to the same coordinate system and the same spatial resolution of SARA AOD by using a cubic spline interpolation algorithm.Finally, hourly PM 2.5 concentration maps at 500 m resolution were produced using SARA AOD products and meteorological field datasets with four models including GTWR, ordinary OLS, GWR, and the TWR model, respectively.Only the data spatially co-located at the selected 37 ground stations and temporally matched with MODIS overpassing moments were included in the following analyses.More details are given in Table 1.

SARA AOD and Quality Contrast
In this study, MODIS Terra/Aqua satellite products were used to generate Terra/Aqua SARA AOD using XuZhou-CUMT AERONET data.High AOD values, 0.60 and 0.71 at 0.55 μm, were observed on 11 January 2015 and 13 February 2015, respectively, during the MODIS Terra/Aqua overpass moments at XuZhou-CUMT AERONET station.Figures 5d and 6d show the image of SARA AOD with 500 m resolution over the studied area on the haze days.Furthermore, MODIS Dark Target (DT) and MODIS Deep Blue (DB) AOD products at 0.55 μm were also obtained to compare with the SARA AOD.For both Terra and Aqua satellites, MODIS DT algorithm (Figures 5b and 6b) has many missing pixels and is incapable of retrieving AOD fully during the haze events owing to its selection of dark pixels when the observed MOD09 surface reflectance product at 2.21 μm is larger than 0.30 over the study region.Although the MODIS DB algorithm (Figures 5c and 6c) has fewer missing pixels and is able to retrieve AOD during the haze events as compared to the MODIS DT algorithm, the retrieved AOD has coarser spatial resolution than SARA AOD.The contrast shows that both

SARA AOD and Quality Contrast
In this study, MODIS Terra/Aqua satellite products were used to generate Terra/Aqua SARA AOD using XuZhou-CUMT AERONET data.High AOD values, 0.60 and 0.71 at 0.55 µm, were observed on 11 January 2015 and 13 February 2015, respectively, during the MODIS Terra/Aqua overpass moments at XuZhou-CUMT AERONET station.Figures 5d and 6d show the image of SARA AOD with 500 m resolution over the studied area on the haze days.Furthermore, MODIS Dark Target (DT) and MODIS Deep Blue (DB) AOD products at 0.55 µm were also obtained to compare with the SARA AOD.For both Terra and Aqua satellites, MODIS DT algorithm (Figures 5b and 6b) has many missing pixels and is incapable of retrieving AOD fully during the haze events owing to its selection of dark pixels when the observed MOD09 surface reflectance product at 2.21 µm is larger than 0.30 over the study region.Although the MODIS DB algorithm (Figures 5c and 6c) has fewer missing pixels and is able to retrieve AOD during the haze events as compared to the MODIS DT algorithm, the retrieved AOD has coarser spatial resolution than SARA AOD.The contrast shows that both MODIS DT and DB algorithms obtained a relatively high estimation of regional AOD, while SARA AOD obtained approximate values referring to ground-based AERONET AOD.In other words, SARA AOD not only displays similar spatial distribution pattern with MODIS DT/DB AOD over the study region, but also provides better accuracy and more details during the haze events.Hence the SARA AOD of 500 m spatial resolution is more suitable for studies on environmental problems since many studies need population-related PM 2.5 exposure information of a polluted factory or residential quarter.
Remote Sens. 2016, 8, 262 10 of 21 MODIS DT and DB algorithms obtained a relatively high estimation of regional AOD, while SARA AOD obtained approximate values referring to ground-based AERONET AOD.In other words, SARA AOD not only displays similar spatial distribution pattern with MODIS DT/DB AOD over the study region, but also provides better accuracy and more details during the haze events.Hence the SARA AOD of 500 m spatial resolution is more suitable for studies on environmental problems since many studies need population-related PM2.5 exposure information of a polluted factory or residential quarter.

Comparison between Fitted and Ground-Observed PM2.5
Hourly PM2.5 concentration values and matched estimations were obtained at the 37 monitoring stations in which regional haze episodes have occurred frequently over the study region from 1 November 2014 to 28 February 2015.About two-thirds of the observed measurements (group-A) extracted from each city is used as training dataset (N = 700) for model fitting and the other one-thirds (group-B) is used as a validation dataset (N = 405) for cross validation.For example if the number of monitoring stations in a city is three or four, the observed data points from the one selected station are used as validation datasets and the observed data points from the remaining two or three stations will be used as training datasets.The OLS, GWR, TWR, and GTWR models were tested using the same datasets.Unlike the OLS method, which directly uses numerical values of variables without any spatial or temporal considerations (taking space coordinates and time as exogenous determinants), the GWR model extends the traditional regression framework for dealing with spatial non-stationarity when different relationships exist at different points in space, and TWR integrates the temporal information in the weighting matrices to capture the temporal heterogeneity when different relationships exist at the same points in time.The results of the model fitting and cross validation are reported in Tables 2 and 3.

Comparison between Fitted and Ground-Observed PM 2.5
Hourly PM 2.5 concentration values and matched estimations were obtained at the 37 monitoring stations in which regional haze episodes have occurred frequently over the study region from 1 November 2014 to 28 February 2015.About two-thirds of the observed measurements (group-A) extracted from each city is used as training dataset (N = 700) for model fitting and the other one-thirds (group-B) is used as a validation dataset (N = 405) for cross validation.For example if the number of monitoring stations in a city is three or four, the observed data points from the one selected station are used as validation datasets and the observed data points from the remaining two or three stations will be used as training datasets.The OLS, GWR, TWR, and GTWR models were tested using the same datasets.Unlike the OLS method, which directly uses numerical values of variables without any spatial or temporal considerations (taking space coordinates and time as exogenous determinants), the GWR model extends the traditional regression framework for dealing with spatial non-stationarity when different relationships exist at different points in space, and TWR integrates the temporal information in the weighting matrices to capture the temporal heterogeneity when different relationships exist at the same points in time.The results of the model fitting and cross validation are reported in Tables 2 and 3.The tables show that the R 2 of GTWR (0.96 and 0.87) is two times more than that of OLS (0.35 and 0.41), while the R 2 of TWR (0.63 and 0.68) is a little bigger than that of GWR (0.59 and 0.60) both in model fitting and cross validation processes.Although the GTWR model is a bit over-fitting, i.e., R 2 decreased 0.09 from model fitting to cross validation, the GTWR obtained the highest value of R 2 in the cross-validation process.Furthermore, the GTWR model shows RMSEs of 11.47 µg/m 3 and 21.77 µg/m 3 , and MAD of 6.91 µg/m 3 and 12.92 µg/m 3 for model fitting and cross validation, respectively, between the ground-observed PM 2.5 and the fitted PM 2.5 concentrations.It shows also the PM 2.5 data estimated from the GTWR model agree well with the ground measurements, indicating that the spatial and temporal variation of model parameters should be considered in the PM 2.5 -AOD relationship.
Additionally, approximately 10.8% and 23.2% of the GTWR PM 2.5 estimation errors are obtained from model fitting and cross validation, respectively, which are far lower than those of the other three models.This proves the ability of the GTWR model to estimate ground-level PM 2.5 concentrations from MODIS SARA-AOD.The comparisons indicate that GTWR is able to accurately capture the spatio-temporal variability of the PM 2.5 -AOD relationship, and the PM 2.5 estimations agree well with the ground measurements over the study region even if haze episodes occur frequently.
To further evaluate the performance of the GTWR model, all the PM 2.5 measurements from 37 monitoring stations were used as a training dataset (N = 1105).Figure 7 quantitatively illustrates the mean difference between the GTWR estimations and ground measurements at 37 monitoring stations.It shows that the GTWR estimations agree well with the ground measurements.The PM 2.5 discrepancies vary from 6.61 µg/m 3 (as the ground measurement is 38.35 µg/m 3 ) to ´11.72 µg/m 3 (as the ground measurement is 117.79 µg/m 3 ).Around 94.6% of the GTWR estimations correlate well with ground measurements at very low discrepancies (e.g., with the error less than 4.0 µg/m 3 ).Figures 8 and 9 give the scatter plots between measured and fitted PM2.5 concentrations with OLS, GWR, TWR, and GTWR model, respectively.All the slopes (0.35, 0.56, 0.61, and 0.94 for model fitting, and 0.38, 0.57, 0.63, and 0.89 for cross validation, respectively, for the OLS, GWR, TWR, and GTWR models) of the above four models are smaller than 1.This characteristic implies that most of the PM2.5 concentrations are systematically underestimated both in model fitting and cross validation.Figures 8a  and 9a show that few estimations of OLS (r = 0.60 and 0.65 for fitting and validation, respectively) lie close to the 1:1 line, and a wide scattering of points is seen due to underestimation and overestimation.Figures 8b and 9c show that the GWR estimations have comparable correlation coefficients (r = 0.77 and 0.80 for fitting and validation, respectively) with TWR estimations, and the r of both GWR and TWR are greater than OSL.This tells us that spatial and temporal variability should be considered in the linear regression model.Figures 8d and 9d show that the majority of the GTWR estimations (r = 0.98 and 0.93 for fitting and validation, respectively) are close to the 1:1 line and fall in the 95% confidence envelope.They also show that when the ground measurement is lower than 100 μg/m 3 , the estimations have a relatively good quality except for OLS; when the ground measurements become larger than 150 μg/m 3 , the performance of GTWR is far better than the other three models.In a word, the GTWR model displays the best performance both for model fitting and cross validation.Figures 8 and 9 give the scatter plots between measured and fitted PM 2.5 concentrations with OLS, GWR, TWR, and GTWR model, respectively.All the slopes (0.35, 0.56, 0.61, and 0.94 for model fitting, and 0.38, 0.57, 0.63, and 0.89 for cross validation, respectively, for the OLS, GWR, TWR, and GTWR models) of the above four models are smaller than 1.This characteristic implies that most of the PM 2.5 concentrations are systematically underestimated both in model fitting and cross validation.Figures 8a  and 9a show that few estimations of OLS (r = 0.60 and 0.65 for fitting and validation, respectively) lie close to the 1:1 line, and a wide scattering of points is seen due to underestimation and overestimation.Figures 8b and 9c show that the GWR estimations have comparable correlation coefficients (r = 0.77 and 0.80 for fitting and validation, respectively) with TWR estimations, and the r of both GWR and TWR are greater than OSL.This tells us that spatial and temporal variability should be considered in the linear regression model.Figures 8d and 9d show that the majority of the GTWR estimations (r = 0.98 and 0.93 for fitting and validation, respectively) are close to the 1:1 line and fall in the 95% confidence envelope.They also show that when the ground measurement is lower than 100 µg/m 3 , the estimations have a relatively good quality except for OLS; when the ground measurements become larger than 150 µg/m 3 , the performance of GTWR is far better than the other three models.In a word, the GTWR model displays the best performance both for model fitting and cross validation.Figures 8 and 9 give the scatter plots between measured and fitted PM2.5 concentrations with OLS, GWR, TWR, and GTWR model, respectively.All the slopes (0.35, 0.56, 0.61, and 0.94 for model fitting, and 0.38, 0.57, 0.63, and 0.89 for cross validation, respectively, for the OLS, GWR, TWR, and GTWR models) of the above four models are smaller than 1.This characteristic implies that most of the PM2.5 concentrations are systematically underestimated both in model fitting and cross validation.Figures 8a  and 9a show that few estimations of OLS (r = 0.60 and 0.65 for fitting and validation, respectively) lie close to the 1:1 line, and a wide scattering of points is seen due to underestimation and overestimation.Figures 8b and 9c show that the GWR estimations have comparable correlation coefficients (r = 0.77 and 0.80 for fitting and validation, respectively) with TWR estimations, and the r of both GWR and TWR are greater than OSL.This tells us that spatial and temporal variability should be considered in the linear regression model.Figures 8d and 9d show that the majority of the GTWR estimations (r = 0.98 and 0.93 for fitting and validation, respectively) are close to the 1:1 line and fall in the 95% confidence envelope.They also show that when the ground measurement is lower than 100 μg/m 3 , the estimations have a relatively good quality except for OLS; when the ground measurements become larger than 150 μg/m 3 , the performance of GTWR is far better than the other three models.In a word, the GTWR model displays the best performance both for model fitting and cross validation.10a from Aqua).Then, the pollution is extended to the south and east, and brings higher PM2.5 to ZaoZhuang, Xuzhou, and the northwest of SuZhou (Figure 10b-e for 9 February 2015 and 10 February 2015).Up to 11 February 2015, the extremely high pollution is centered at ZaoZhuang, JiNing, and the northeast part of HeZe, with PM2.5 concentrations beyond 150 μg/m 3 .This shows that the PM2.5 concentrations over the study region slowly accumulated during this haze process except for the dramatic increase in ZaoZhuang and JiNing, which indicates that both local and external pollutant sources had an impact on the PM2.5 concentration.
Furthermore, it shows that the AOD-estimated PM2.5 concentration from Terra at about 10:30 a.m. is higher than that from Aqua at about 1:30 p.m. in the same day, and these temporal characteristics agree with the ground measurements over the study region (Figure 10b,c on 9 February 2015, and Figure 10d,e on 10 February 2015).Moreover, more spatial details and higher estimation accuracy are obtained near the monitoring stations.That is to say, if there are more ground measurements with spatial homogeneous distribution over the study region, the estimation accuracy of the GTWR model can be improved further.The AOD-estimated PM2.5 images with GTWR model are also compared with spatially interpolated ground measurements (Figure 11) with the ordinary Kriging method.The ground measurements can be regarded as relatively independent and the ordinary Kriging method is well proved in the literature [45,46], being more suitable than the inverse distance weighted (IDW) and three splines (TS) methods embodied in ArcGIS, in that the ground measurements are obtained from sparse and irregularly distributed air quality monitoring stations.Although the spatio-temporal patterns of ordinary Kriging-interpolated PM2.5 measurements are roughly similar to those of GTWRestimated PM2.5 concentration, the former have poor accuracy and less spatial details in places far from the monitoring stations, which might be attributable to the following reasons: (i) The ordinary Kriging-interpolated PM2.5 concentrations cannot exceed the maximum and minimum values of the ground measurements (e.g., the southeast of ShangQiu on 10 February 2015 and the northwest of XuZhou on 11 February); (ii) the coarse spatial coverage and irregular distribution of ground stations may introduce certain estimation errors into the ordinary Kriging method.Furthermore, it shows that the AOD-estimated PM 2.5 concentration from Terra at about 10:30 a.m. is higher than that from Aqua at about 1:30 p.m. in the same day, and these temporal characteristics agree with the ground measurements over the study region (Figure 10b,c on 9 February 2015, and Figure 10d,e on 10 February 2015).Moreover, more spatial details and higher estimation accuracy are obtained near the monitoring stations.That is to say, if there are more ground measurements with spatial homogeneous distribution over the study region, the estimation accuracy of the GTWR model can be improved further.
The AOD-estimated PM 2.5 images with GTWR model are also compared with spatially interpolated ground measurements (Figure 11) with the ordinary Kriging method.The ground measurements can be regarded as relatively independent and the ordinary Kriging method is well proved in the literature [45,46], being more suitable than the inverse distance weighted (IDW) and three splines (TS) methods embodied in ArcGIS, in that the ground measurements are obtained from sparse and irregularly distributed air quality monitoring stations.Although the spatio-temporal patterns of ordinary Kriging-interpolated PM 2.5 measurements are roughly similar to those of GTWR-estimated PM 2.5 concentration, the former have poor accuracy and less spatial details in places far from the monitoring stations, which might be attributable to the following reasons: (i) The ordinary Kriging-interpolated PM 2.5 concentrations cannot exceed the maximum and minimum values of the ground measurements (e.g., the southeast of ShangQiu on 10 February 2015 and the northwest of XuZhou on 11 February); (ii) the coarse spatial coverage and irregular distribution of ground stations may introduce certain estimation errors into the ordinary Kriging method.In conclusion, GTWR is capable of estimating ground-level PM 2.5 concentrations with good accuracy in regions lacking ground measurements but with satellite AOD available, and the estimated PM 2.5 images can display the detailed spatial patterns and dynamic evolution process of regional PM 2.5 concentrations.

Discussions
The GTWR model was employed to depict the relationship between satellite-derived AOD and ground-level PM 2.5 concentrations during the period between November 2014 and February 2015.The proposed GTWR model and estimation method in this paper has several benefits over the conventional statistical models used in previous studies.
Firstly, we use high spatial resolution SARA AOD of 500 m resolution to estimate ground-level PM 2.5 concentrations.Compared to ground-observed AOD from AERONET, the SARA AOD has much higher estimation accuracy than that from MODIS 04 algorithm during fine-particle pollution events [47], which indicates that SARA AOD is more suitable for estimating ground-level PM 2.5 concentrations.On the other hand, SARA AOD has a considerably larger spatial coverage than MODIS DT/DB AOD over areas with high surface reflectance or undergoing a haze episode.Therefore, the SARA-retrieved 500 m AOD and GTWR-estimated ground-level PM 2.5 concentrations can be utilized to monitor the urban air pollution status at a finer spatio-temporal scale.
Secondly, we use four models to depict the spatio-temporal dynamics in the relationships between satellite-retrieved AOD and ground-level PM 2.5 estimations.From Tables 2 and 3 we know that the estimation accuracy of TWR is higher than GWR, implying that the temporal non-stationarity has a more significant impact on model fitting based on the training datasets than spatial non-stationarity in this case study, which may be attributable to the following reasons: (i) in order to monitor the daily dynamic of air quality, the observations from MODIS TERRA and AQUA sensors are input to the statistic models for estimating ground-level PM 2.5 concentrations; (ii) haze episodes occur frequently in winter, resulting in low AOD value but high PM 2.5 concentrations; (iii) the topography and land use in the study region are pretty much even, without dramatic differences.The GTWR model allows for hourly-to-hourly variability in the PM 2.5 -AOD relationship by incorporating temporal effects into the GWR model.The lowest RMSE values and the high correlation coefficients from GTWR-based estimation illustrate that the fitted and observed PM 2.5 have a strong relationship (Figures 8 and 9).
Thirdly, we use hourly assimilated meteorological field data including PBLH, RH, TEMP, and WS extracted from WRF assimilation as auxiliary data to improve the estimation accuracy.The good performance of the GTWR-based estimation method shows that ground-observed meteorological field data can be substituted for simulated meteorological field datasets if hourly ground meteorological observations are lacking.

Conclusions
In this study, a GTWR model was proposed to incorporate the geographical difference and temporal variation of the PM 2.5 -AOD relationship, and to estimate ground-level PM 2.5 concentration from satellite-retrieved AOD products.SARA AOD of 500 m resolution is produced from MODIS observations and input to the GTWR model together with assimilated meteorological data, including the PBLH, RH, WS, and TEMP from the WRF model, to estimate ground-level PM 2.5 concentration.Other models such as OLS, GWR, and TWR were applied for comparison with a common case study in central China.This showed that the GTWR model produces the highest R 2 and r values as well as the lowest RMSE, MAD, and MAPE values among the four models.GTWR provides a simple but effective method of estimating hourly ground-level PM 2.5 concentrations at 500 m spatial resolution, and can substantially increase the estimation accuracy.The estimations are valuable for urban-and regional-scale air quality monitoring and government decision-making, especially for rural regions lacking ground measurements.
However, some limitations/errors remain in GTWR estimation because (i) no effective method is applied to fill the gaps in pixels where AOD is not retrieved due to SARA being unable to estimate AOD in missing pixels caused by the presence of thick clouds; (ii) no finer correction is made with LUCC data when a pixel involves water surface; (iii) the pixels covered by snow are temporally suggested to be ignored and set to null values in the current model and the model will be improved to be accountable for this special situation; (iv) the ground PM 2.5 monitoring stations are limited and irregularly distributed, which may reduce the estimation accuracy in areas far from stations; (v) it is effective in cases where the near-surface aerosol is the dominant load in the vertical direction, but will not work in long-range transported aerosol episodes such as dust storms for its entrance to the middle troposphere, when the-satellite observed AOD becomes less relevant to surface PM 2.5 concentration; and (vi) the spatial resolution and accuracy of meteorological field data simulated in the WRF model are not good enough.Future strategies for improvement include: (i) ground-level PM 2.5 estimation from other high spatial resolution satellite sensors, such as Landsat and HJ-1A, for the study of microenvironments of human exposure (e.g., coal mining areas, power plants, business quarters, and residential areas); (ii) more factors such as land use, wind direction, and fire product being cooperated in the GTWR model; (iii) physical and chemical properties of aerosol particulates being explored for PM 2.5 estimation; and (iv) air quality models such as temporal gap-filling being coupled with the GTWR model to estimate ground-level PM 2.5 concentrations at nighttime and in pixels where the AOD retrieval is not valid for a given day.

Figure 1 .
Figure 1.Location of AERONET station and environmental monitoring stations in the JSHA study region.

Figure 1 .
Figure 1.Location of AERONET station and environmental monitoring stations in the JSHA study region.

Figure 4 .
Figure 4. Flowchart of the proposed GTWR method for estimating hourly PM2.5 at 500 m resolution.

Figure 4 .
Figure 4. Flowchart of the proposed GTWR method for estimating hourly PM 2.5 at 500 m resolution.

Figure 7 .
Figure 7.Comparison between GTWR estimations and ground measurements of PM 2.5 concentration.

3. 3 .
Figure 10 illustrate the spatio-temporal distribution of PM2.5 hourly concentrations estimated from MODIS SARA-AOD with GTWR model during a haze formation process from the afternoon of

3. 3 .
Figure 10 illustrate the spatio-temporal distribution of PM2.5 hourly concentrations estimated from MODIS SARA-AOD with GTWR model during a haze formation process from the afternoon of

Figure 11 .
Figure 11.Spatial distribution of ground PM2.5 measurements over the study region during a haze formation process from the afternoon of 8 February 2015 to the morning of 11 February 2015, referring to the local passing moment of satellites Terra and Aqua ((a) Aqua 8 February 2015; (b) Terra 9 February 2015; (c) Aqua 9 February 2015; (d) Terra 10 February 2015; (e) Aqua 10 February 2015; and (f) Terra 11 February 2015).In situ measurements at training stations (circles) and validation (triangles) stations are also shown in the figures.

Figure 11 .
Figure 11.Spatial distribution of ground PM 2.5 measurements over the study region during a haze formation process from the afternoon of 8 February 2015 to the morning of 11 February 2015, referring to the local passing moment of satellites Terra and Aqua ((a) Aqua 8 February 2015; (b) Terra 9 February 2015; (c) Aqua 9 February 2015; (d) Terra 10 February 2015; (e) Aqua 10 February 2015; and (f) Terra 11 February 2015).In situ measurements at training stations (circles) and validation (triangles) stations are also shown in the figures.

Table 1 .
A summary of data variables used in the GTWR model.

Table 1 .
A summary of data variables used in the GTWR model.