Analysis of PM 2.5 Variations Based on Observed, Satellite-Derived, and Population-Weighted Concentrations

: Fine particulate matter (PM 2.5 ), which can cause adverse human health effects, has been proven as the ﬁrst air pollutant in China. In situ observations with ground-level monitoring and satellite-based concentrations have been used to analyze the variations in PM 2.5 . However, variation analyses based on these two kinds of measurement have mainly focused on the concentration itself and ignored the effects on the population. Therefore, this study not only investigated these two kinds of measurements, but also performed weighted population analyses to study the variations in PM 2.5 . Firstly, daily models of timely structure adaptive modeling (TSAM) were constructed to simulate satellite-derived PM 2.5 levels from January 2013 to December 2016. Secondly, population-weighted concentrations were calculated based on TSAM-derived PM 2.5 surfaces. Finally, observed, TSAM-derived, and population-weighted concentrations were used to analyze the variations in PM 2.5 . The results showed the different importance of various input parameters; AOD had the highest rank. Additionally, TSAM models demonstrated good performance, ﬁtting R ranging from 0.86 to 0.91, and validating R from 0.82 to 0.89. According to the air quality standard in China, TSAM-derived PM 2.5 showed that the increase in area lower than Level II was 29.03% and the increase in population was only 14.81%. This indicates that the air quality exhibited an overall improvement in spatial perspective, but some areas with high population density showed a relatively low improvement due to uneven distributions in China. The population-weighted PM 2.5 concentration could better represent the health threats of air pollutants compared with in situ observations.


Introduction
Fine particulate matter (PM 2.5 ) with micro particle sizes can cause serious health consequences, such as morbidity, mortality, and cardiovascular and respiratory diseases [1][2][3][4]. Except these hazards, it could also be the corrosion of materials [5]. With rapid economic development and urbanization in recent decades, China has been facing the growing problem of poor atmospheric environment [6,7]. Air pollution not only poses threats to public health but also affects daily production and life, such as atmospheric visibility reductions and the destruction of the urban landscape [7,8].
To improve air quality, national and local governments have promulgated a series of plans and policies since 2012. Among them are the 12th Five-Year Plan on Air Pollution Prevention and Control in Key Regions (abbreviation: Key Regions) and the most recent Chinese Ambient Air Quality Standard (CAAQS) released on 10 September 2013. At the same time, nationwide air quality monitoring networks have gradually been established since 2013. Increasing numbers of studies have been conducted to analyze the temporalspatial variations in PM 2.5 based on ground-level monitoring measurements [9][10][11].
However, site-based observations, which collect concentrations of air pollutants from fixed sites, merely reflect ambient PM 2.5 concentrations at limited distances around the ground monitoring stations and hardly reflect spatial variations in air pollutants [12]. Additionally, researchers have focused on spatial interpolation methods to reveal regional variations in PM2.5, e.g., Kriging and inverse distance weight interpolation [13,14]. Spatial interpolation methods can derive surfaces of PM 2.5 concentrations. However, the uncertainties in these methods could be high when varying topography, meteorology, and local emissions among the monitoring sites are not accounted for [14][15][16][17].
A previous study proved that remote sensing could provide more information about the Earth's climate parameters (e.g., aerosols, temperature, pressure, etc.) [18]. Previous studies have proved that aerosol optical depth (AOD) and PM 2.5 have similar particle sizes, which provide the theoretical basis of modeling between PM 2.5 and satellite AOD [19,20]. At the same time, AOD with wide spatial coverage and acceptable accuracy which is satellite-retrieved can provide planar data from the macro scale and supplement ground observations. It has been increasingly utilized to simulate the regional distribution of PM 2.5 [6,21,22]. Spatial-temporal statistical modeling is one of the most popular methods in PM 2.5 concentration predictions and includes simple linear regression [23,24], multiple linear regression (e.g., land use regression (LUR) [14,17], geographically weighted regression (GWR) [21,25], linear mixed effect (LME) [26,27], timely structure adaptive modeling (TSAM) [22]), and artificial intelligence models (e.g., ANN) [16,[28][29][30]. While the contributing strength of model predictors (i.e., variables) are allowed to change in these daily modeling process, all these models (no matter whether they are linear or non-linear) for individual days over a study period have a fixed structure. In fact, both the models' predictors and their contributing strengths could vary with time. Then, TSAM takes temporal variations and magnitudes of predictors into consideration in the modeling process and provides a novel solution for continuously mapping the distributions of PM 2.5 [22].
Existing research has mainly focused on the variations in in situ or satellite-derived concentrations of PM 2.5 and then analyzed its spatial and temporal characteristics. Population density in China is not evenly distributed, and Hu demonstrated that 93.77% of the population live south-east of this "Hu Line", which is a contrast line of population density proposed by Hu in 1935 [31]. Regions with high population density would therefore have higher population exposure risks with the same concentrations of pollutants. Additionally, when analyzing the effect of air pollution, we not only need to analyze concentration variations in air pollutants, but we also need to consider population exposure. Therefore, the main aims of this study are: (1) to simulate the satellite-based distribution of PM 2.5 from 2013 to 2016 with the TSAM modeling framework; (2) to analyze the spatial-temporal variations in observed and TSAM-derived PM 2.5 ; and (3) to compare the difference between observed and population-weighted PM 2.5 levels for Key Regions.  Hourly ground PM 2.5 concentrations were collected from the official database, which had 661 fixed sites in 2013 and 1436 by the end of 2016. These were measured with the tapered element oscillating microbalance method (TEOM). Invalid and error measurements caused by instrument calibration failure and power failure were reduced based on the Chinese National Ambient Air Quality Standards (CNAAQS, GB3095-2012). In this process, repeatedly reported records for several successive hours, significantly abnormal measurements (e.g., more than 1000 µg/m 3 or less than 0 µg/m 3 ), or fewer than 20 records of PM 2.5 concentrations in a day were deleted. As a result, about 1% of records were discarded as abnormal during the study period. Screened hourly PM 2.5 measurements were averaged into daily values for each monitoring station. To maintain a consistent spatial resolution of PM 2.5 with the predictor of AOD, these daily PM 2.5 values were recalculated based on the 10 km grid cells. Additionally, PM 2.5 concentrations were calculated for grid cells with screened hourly observations which matched the satellite overpass times. Blue_Combined, NASA, Washington, DC, USA) with the highest quality assurance confidence flags (QA = 3) were screened as AOD values on specific days. To improve the spatial coverage, Terra AOD (MOD04) and Aqua AOD (MYD04) were extracted and combined, following an existing research method [22,32,33].

PM 2.5 Emission Data
Due to the unavailability of the national emission inventory, land use, traffic networks, and demographic data were regarded as alternatives for indirectly characterizing the PM 2.5 emissions. Population raster data with a resolution of 1 km in 2010 were collected and annually updated concerning the county-level population growth rate from 2013 to 2016. The variables of land use included the percentage of built-up area, forest area, grass area, and water area. The traffic variable was the total road length. Additionally, population density values were treated as proxy variables of demographic data.

PM 2.5 Dispersion Conditions Data
It has been acknowledged that meteorology and topography are key factors that can directly influence the dispersion of PM 2.5 [21,34]. About 900,000 daily meteorological records, including wind speed (WS, m/s), precipitation (PE, mm), relative humidity (RH, %), temperature (TEMP, K), and pressure (PS, hPa), were obtained from 824 sites of the meteorological monitoring network. Abnormal values were discarded in this process based on the quality control code and standard outlier elimination. These meteorological factors were interpolated with the inverse distance-weighted method into a 10 km raster. Moreover, as an indicator of topography, digital elevation model (DEM) data from the Shuttle Radar Topography Mission (SRTM) were collected with a spatial resolution of 90 m.

Structure of the TSAM Model
TSAM, a model which has been proven to be effective in PM 2.5 concentration estimations, was employed for satellite mapping in this study; details of TSAM were presented by Fang et al. [22]. The general formulation of TSAM is expressed in Equation (1), which consists of the dependent variable of PM 2.5 and independent variables including AOD, emission-related predictors, and dispersion-condition-related predictors.
In this case, PM 2.5 indicates the daily PM 2.5 concentrations; Satellite means the daily satellite-retrieved AOD; Emissions refers to the proxy variables of land use, road length, and population; and Dispersion includes variables of meteorology and topography. The basic structure for the daily TSAM model can be written as in Equation (2).
where PM 2.5gd represents the daily ground PM 2.5 measurements at cell g on day d, α 0,gd is the location-specific intercept, α 1,gd~α13,gd are the location-specific slopes for the corresponding predictors, AOD gd is used as a required variable for a specific model which denotes the MODIS AOD value, the flag t in the predictor means that the predictor is not constantly considered in the final model, and whether the predictor is selected or not will be determined based on the daily model performance. Var with the index i from 2 to 6 rep-resents the meteorological conditions, including wind speed, precipitation, temperature, pressure, relative humidity: Var with an index of 7 denotes the topography condition with the value of elevation; Var with an index of 8 is the total road length within the buffer with a radius of 5 km; Var with an index, i, from 9 to 12 represents the percentages of built-up, forest, grass, and water areas within the 5 km buffer around the site, respectively; and Var with an index of 13 denotes the population. ξ gd is the error term. Time average maps of all these input variables (Figures S1-S12) are provided in the Supplementary Materials.

TSAM Model Fitting, Validation, and Prediction
Firstly, all spatial data used in this study were unified to the same projection coordinate system (Asia_Lambert_Conformal_Conic). Then, the daily PM 2.5 values were resampled to 10 km for consistency with the size of the AOD pixel. Based on the in situ sites for regulatory PM 2.5 observations, variables of land use, road length, and population density were extracted using the buffer and overlay analysis of ArcGIS 10.1. The buffer radius for land use and road length was set as 5 km. The monitoring daily meteorological measurements were also processed into 10 km resolution using the IDW method. Subsequently, these variables were finally extracted to daily sample files. Details of the variable extraction and calculations were presented in our previous study [22]. The stepwise regression method was used to select the predictors in the specific modeling day. The optimal model structure was determined in line with comparison results of indicators of model performance (e.g., R 2 , RMSE, etc.) among different model structure combinations and used to fit the TSAM model for specific days.
Tenfold cross-validation tests were carried out, which randomly selected 90% as the modeling set and 10% as the validating set from the total samples for each modeling day. We repeated this step ten times until all samples used as validation sets were used to test the model robustness and modeling performance; detailed information refers to existing work [35]. Indicators, including the RMSE, mean prediction error (MPE), relative prediction error (RPE, which is defined as the RMSE divided by the mean observed value of samples), relative mean prediction error (RMPE, which is defined as the MPE divided by the mean observed value of samples), and R 2 , were used to indicate the performance of the TSAM model. Finally, a 10 km fishnet was created, with 92,122 grid cells in total. PM 2.5 values for each fishnet cell were estimated based on the validated TSAM model and associated modeling variables.

Calculation of Population-Weighted PM 2.5 Concentration
In this study, we not only considered the variations in PM 2.5 observations and TSAMderived concentrations analyzed, but also evaluated the exposure risk to PM 2.5 . We calculated population-weighted concentrations as indicators of exposure based on TSAMderived PM 2.5 concentrations [36]; the equation is shown below: where PopCon represents the population-weighted concentration of the target region, Pop i represents the population in grid i, and Con i is the TSAM-derived PM 2.5 concentration in grid i.

Analysis of TSAM Model Structure
Based on the framework of the TSAM model [22], 972 daily models were successfully constructed from 1 January 2013 to 31 December 2016. Additionally, the remaining days could not be modeled owing to the lack of sufficient samples. The modeling results demonstrated that AOD was the key variable with the highest rank among all the input parameters. The meteorological variables (wind speed, relative humidity, and precipitation) and land-use variables (e.g., percentage of built-up, forest, and grass) also exhibited high ranks. Other input parameters were relatively less used.
The spatial distributions of the intercept, variable coefficient, and corresponding standard error of the daily fitting model (e.g., 16 May 2016) are shown in Figure 2, presenting significant spatial aggregation. Intercepts ranged from 9.02 to 171.38, with a mean value of 29.79. The coefficient of AOD ranged from −0.2021 to 0.0058, and wind speed ranged from −31.08 to 15.55. A positive value of the variable coefficient indicated that it played a positive effect on PM 2.5 , and vice versa. The standard error of the intercept ranged from 2.49 to 29.75. Standard errors of the AOD coefficient ranged from 0.0035 to 0.028, and the values of wind speed ranged from 0.7227 to 9.4729. The standard error of coefficients, which measured the reliability of each coefficient estimate, were relatively lower than values of coefficients. At the same time, the residual standard deviations of the TSAM model in the modeling day ranged from −2.66 to 4.84; 97.76% of the modeling samples had residual standard deviations between −2.5 and 2.5.       Figure 5 shows the spatial distributions of TSAM-derived annual means over mainland China from 2013 to 2016. The blank regions in Figure 5 mean that there were no available annual mean values of TSAM-derived PM 2.5 , which was due to the missing satellite AOD with the influences of cloud or surface conditions (e.g., ice, snow, etc.). In Figure 5, PM 2.5 exhibits significant spatial aggregation. Overall, heavily PM 2.5 -polluted areas were centrally located in the north plain, center eastern part, the Sichuan Basin, and the Tarim Basin. The Beijing-Tianjin-Hebei urban agglomeration (BTH) was one of the most heavily polluted regions, with annual PM 2.5 levels ranging from 45 µg/m 3 to 105 µg/m 3 . Meanwhile, the Shandong Peninsula urban agglomeration and the Chengdu-Chongqing metropolitan area also had non-ignorable levels of PM 2.5 pollution, with annual means around 45~95 µg/m 3 . The spatial distributions of PM 2.5 also demonstrated obvious heterogeneous characteristics. Specifically, there were significant differences in PM 2.5 between the areas north and south of the Yangtze River; generally, the former levels were higher than the latter. At the same time, there were also clear differences within the same administrative region. For Beijing, the values of PM 2.5 ranged from 40 to 70 µg/m 3 in 2016. In terms of the spatial distribution, areas with high concentrations (Dongcheng, Xicheng, Chaoyang District, etc.) were located in the centers of the cities, and areas with low concentrations were located towards the north-west, such as Yanqing, Miyun, and Huairou District. Additionally, according to differences in TSAM-derived PM 2.5 between 2013 and 2016, over 85% of areas exhibited declines greater than 10 µg/m 3 in the study area. The decline in PM 2.5 ranged from 10 to 30 µg/m 3 . The Chengdu-Chongqing metropolitan region also exhibited a strong decline in PM 2.5 levels, with average decreases of about 15 to 30 µg/m 3 .

PM 2.5 Variation Analysis Based on Percentage of Area and Population
Summaries of the percentage of area and population exceeding and under air quality Level II standard (annual mean PM 2.5 concentration lower than 35 µg/m 3 ) based on TSAMderived PM 2.5 are shown in Figure 6. As Figure 6a

Comparison between Observed and Population-Weighted PM 2.5 Values for Key Regions
The national population-weighted PM 2.5 values decreased from 71.99 µg/m 3 to 48.04 µg/m 3 , with a decrement of 33.27%. To compare the differences among observed and population-weighted PM 2.5 , Table 2

Discussion
This study estimated the spatial-temporal distributions of satellite-derived PM 2.  Table S1. Among 65 variable combinations, the minimum, lower quartile, average, upper quartile, and maximum values of the model with AOD as the only explanatory variable were 0. 23, 8.47, 8.48, 8.62, and 8.66 µg/m 3 , respectively. With the addition of meteorological factors in explanatory variables, not all of these factors could reduce the model uncertainty. The z-score of temperature with the spatial autocorrelation (Moran's I) test was 37.51 Standard Deviations (std), which indicated that it had a significant spatial autocorrelation, and then it was not included in the TSAM model. Among the meteorological factors, only wind speed was added in the list of explanatory variables, and the model uncertainty decreased. This means that the wind speed, which could change the dispersion condition, had a significant effect on PM 2.5 variations. The addition of the precipitation variable did not reduce the standard error, which could be due to its transient effects on PM 2.5 such as the hourly scale, but the TSAM model in this study was developed with the daily scale. With the quantitative comparison of the standard error of 65 variable combinations, we found the model performance did not necessarily get better as the variables increased. The screening of explanatory variables is an important process for PM 2.5 simulation. With the screening, the explanatory variables of the TSAM model in this study day were defined as AOD, wind speed, percentage of forest, and percentage of water.
Both TSAM-derived and observed PM 2.5 values exhibited a significant decline during the study period, which proved that implementing air pollution prevention and control policies had a significant effect on improving the air quality. The annual and monthly concentrations of TSAM-derived PM 2.5 were lower than those observed for each region. The reasons for this could be that in situ observations revealed PM 2.5 levels near the ground monitoring station, which were generally located in urban areas with high PM 2.5 pollutant emissions.
TSAM-derived PM 2.5 showed a strong signature of seasonal variations, with the highest concentrations in winter and the lowest concentrations in summer, and intermediate ranges in spring and autumn during the study period. Winter heating could be the main cause of the high PM 2.5 levels, especially in northern and western China. During seasons when more heating was needed, the pollutant emissions sharply increased, and the PM 2.5 concentration rose correspondingly. In addition to free heating policies, the Spring Festival, which lasts from January to February, could be another key cause of seriously increased pollution in winter. During the Spring Festival, people set off large quantities of fireworks in celebration. Thus, PM 2.5 concentrations sharply increase in January, which is the most polluted month of the year. Meanwhile, there were some areas, especially in northern China, which suffered serious PM 2.5 pollution in spring (from March to May) due to dust storms. The percentage of polluted area increased by 29.03%, whereas the percentage of the population living under the Level II standard only increased by 14.81% during the study period. The difference between the percentages of area and population is attributable to the spatial heterogeneity of population distribution, i.e., some regions had relatively low population density and their air quality had been improved from 2013 to 2016. Moreover, the spatial continuity of air pollutant concentrations may be another reason for this result. This means that the implementation of air pollution prevention and control measures could not only directly improve local air quality, but also affect surrounding air quality at the same time.
The population-weighted PM 2.5 concentrations of the 13 Key Regions showed declines of more than 10 µg/m 3 from 2013 to 2016, except for the central and northern regions of Shanxi. Correspondingly, the declining proportion in population-weighted PM 2.5 was more than 15%. The declines in in situ observations of PM 2.5 ranged from 3.94 to 45.19 µg/m 3 for the 13 Key Regions. Ground-level monitoring measurements could only represent air pollutant concentrations in limited ranges around the stations; however, population-weighted PM 2.5 not only considered the wider spatial range but also took the population into consideration, which could better represent the health threats of air pollutants compared with merely in situ observations. The population-weighted PM 2.5 with a spatial resolution of 10 km could offer support for the large-scale analysis of air pollutant variations, but is relatively coarse to reflect features of the local variation in PM 2.5 . Urban air quality has clear spatial differentiations among small intra-urban scenarios. The limitations of operational satellite product could be improved with the development of a satellite retrieval algorithm. Moreover, the temporal variation analyses of PM 2.5 were implemented under relatively long-term temporal scales, such as monthly, seasonal, and annual. Long-term temporal scales could be used to represent the global trends in PM 2.5 variation, as well as to capture the short-term variation. With the development of stationary satellites, the temporal scale of satellite products can be improved to hourly with continuous sequences, which would capture dynamic variation in air pollutant levels. Additionally, this would greatly enhance the significance of satellitebased modeling and estimations of PM 2.5 for real-time air pollution exposure assessments.
More importantly, satellite data coverage could cause misunderstandings in annual mean difference analyses. For example, the value for Central Liaoning exhibited an increasing trend from 2013 to 2014 due to the low coverage of AOD systems in the winter of 2013 and a relatively high coverage in the winter of 2014. Thus, new algorithms for remedying the coverage of satellite-based AOD values are also urgently required for contemporary PM 2.5 modeling and estimation.

Conclusions
This study developed a national daily TSAM model for PM 2.5 concentration prediction using variables of satellite AOD, meteorology, land-use, etc., from January 2013 to December 2016. The TSAM model demonstrated good performance by taking timely variations and magnitudes of predictors into consideration in the modeling process. The TSAM-derived PM 2.5 could represent a large spatial range and reported consecutive spatial variations, which complements site-based observations of PM 2.5 . Both in situ observations and TSAMderived concentrations proved that the air quality showed great improvements from 2013 to 2016, due to the implementation of air pollution prevention and control policies. According to the air quality Level II standard, the majority of the population was subjected to PM 2.5 concentrations above this standard. The population based on the satellite-derived concentrations represents a complementary indicator to analyze the concentrations of air pollutants.