Spatio–Temporal Relationship and Evolvement of Socioeconomic Factors and PM2.5 in China During 1998–2016

A comprehensive understanding of the relationships between PM2.5 concentration and socioeconomic factors provides new insight into environmental management decision-making for sustainable development. In order to identify the contributions of socioeconomic development to PM2.5, their spatial interaction and temporal variation of long time series are analyzed in this paper. Unary linear regression method, Spearman’s rank and bivariate Moran’s I methods were used to investigate spatio–temporal variations and relationships of socioeconomic factors and PM2.5 concentration in 31 provinces of China during the period of 1998–2016. Spatial spillover effect of PM2.5 concentration and the impact of socioeconomic factors on PM2.5 concentration were analyzed by spatial lag model. Results demonstrated that PM2.5 concentration in most provinces of China increased rapidly along with the increase of socioeconomic factors, while PM2.5 presented a slow growth trend in Southwest China and a descending trend in Northwest China along with the increase of socioeconomic factors. Long time series analysis revealed the relationships between PM2.5 concentration and four socioeconomic factors. PM2.5 concentration was significantly positive spatial correlated with GDP per capita, industrial added value and private car ownership, while urban population density appeared a negative spatial correlation since 2006. GDP per capita and industrial added values were the most important factors to increase PM2.5, followed by private car ownership and urban population density. The findings of the study revealed spatial spillover effects of PM2.5 between different provinces, and can provide a theoretical basis for sustainable development and environmental protection.


Introduction
The rapid development of China's economy in recent decades has caused serious environmental pollution, among which atmospheric pollution is particularly serious [1,2]. The frequent haze weather across the country has seriously affected the urban environment [3,4] and the physical and mental health of residents [5][6][7]. The main pollutant forming haze weather is fine particulate matter with a diameter of less than 2.5 µm (PM 2.5 ). PM 2.5 can reduce visibility. It is harmful to people's life, especially in health effects [8,9], so PM 2.5 pollution has become a research hotspot. The previous studies mainly involved two aspects, namely micro aspects and macro aspects. The micro aspects focus mainly on chemical components [10][11][12] and physical and mental health effects of PM 2.5 [13][14][15][16], etc. The macro aspects focus mainly on the influencing factors [17,18], spatiotemporal variations and The socioeconomic data were obtained from the National Bureau of Statistics of the People's Republic of China (http://data.stats.gov.cn/). In this study, four major socioeconomic factors were collected in 31 provinces of Mainland China during 1998-2016, namely GDP per capita (GDPP), industrial added values (IVA), urban population density (UPD) and private car ownership (PCO). In order to eliminate the influence of dimension, z-scores were used to standard all factors and variables.

Methods
Five methods were used in this paper, namely unary linear regression model, Spearman's rank correlation analysis, univariate spatial autocorrelation, bivariate spatial correlation analysis and spatial regression analysis. The main methods were analyzed in detail below.

Unary Linear Regression Model
In order to analyze the temporal trend, the slope of socioeconomic factors and PM 2.5 concentration were calculated by using unary linear regression model. The slope is expressed as: where slope is the trend gradient, Y t denotes the variable (PM 2.5 concentration or GDPP or IVA or UPD or PCO) in the t-th year, T is the study period of 1998-2016. A positive (negative) slope means that the variable increases (decreases) over the years. The greater the absolute value of the slope, the faster the increase or decrease of speed.

The Univariate Spatial Autocorrelation Analysis
Moran's I, as the most commonly used indicator of global spatial autocorrelation, was initially suggested by Moran [42]. In essence, it represents the cross product statistics of a variable and its spatial lag. The degree to which the feature values of a position are similar or different from those of its spatial neighbors is measured by spatial autocorrelation. The global spatial association of PM 2.5 concentration across China was explored by global Moran's I in this paper. To explore the local spatial association (spatial clustering or spatial dispersion) in adjacent provinces, we chose a local indicator of spatial association (LISA) [43] as the analysis method. The global Moran's I and local Moran's I are calculated by: where I stands for global Moran's I for the whole study region, I i is the Moran's I for province i, x i donates PM 2.5 concentration at province i, x j donates PM 2.5 concentration at all the other provinces (where j = i). Also, x is the mean PM 2.5 concentration of 31 provinces in China, n represents the total number of provinces. σ is the standard deviation of the PM 2.5 concentration of 31 provinces. w ij is the spatial weight matrix, representing province i is adjacent to province j, neighboring provinces were 1 and non-adjacent provinces were 0. The values of I or I i ranged from −1 to 1. A positive (negative) I or I i value indicates positive (negative) spatial autocorrelation in the provinces. Positive autocorrelation indicates that provinces with similar PM 2.5 concentration are closely distributed in space, whereas negative spatial autocorrelation indicates that PM 2.5 concentration of neighboring provinces are dissimilar. A zero I or I i value indicates a random spatial pattern. The size of the absolute value of I or I i can reflect the strength of the spatial correlation.

The Bivariate Spatial Correlation Analysis
The spatial correlation between PM 2.5 and socioeconomic factors were examined by global bivariate Moran's I and local bivariate Moran's I. Global bivariate Moran's I reflects the global spatial associations between PM 2.5 and another variable (GDPP, IVA, UPD or PCO) across the whole region, whereas local bivariate Moran's I explores the local spatial correlations within different provinces [44][45][46]. Global bivariate Moran's I and local bivariate Moran's I are given by: where I xy is the global bivariate Moran's I, and I i xy is the local bivariate Moran's I in province i. n is the total number of provinces, and w ij is the queen contiguity weight matrix. Z x i is the standardized z-scores of PM 2.5 concentration in the i-th province, Z y j is the standardized z-scores of socioeconomic factors (GDPP, IVA, UPD or PCO) in the j-th province. The values of I xy or I i xy is in the range [−1,1]. The values of I xy or I i xy greater than 0, less than 0, equal to 0 indicate positive spatial correlation, negative spatial correlation, or no correlation between PM 2.5 concentration and socioeconomic factors, respectively. The size of the absolute value of I xy or I i xy can reflect the strength of the spatial correlation.

The Spatial Regression Model
Spatial lag model (SLM) and spatial error model (SEM) were based on the ordinary least squares (OLS) [47,48]. SLM can be used to explore whether PM 2.5 concentration diffuses in one province, whereas SEM can be used to interpret the dependence of spatial error [44]. The SLM and SEM can be defined as follows: where Y it denotes PM 2.5 concentration in province i in the t-th year, α represents a constant term. β 1 , β 2 , β 3 , and β 4 are the parameters to reveal the correlations between PM 2.5 and GDPP, IVA, UPD, and PCO, respectively. wY it is a spatial lag-dependent variable vector, it reflects the endogenous interaction effects among Y it , ρ is a spatial regression coefficient that denotes the spatial dependence of the sample observations. w µ reflects the interaction effects among the disturbance term of different provinces. The spatial autoregressive coefficient λ denotes the spatial dependence of the residuals; ε is the random error term, µ represents the spatially autoregressive error terms. In order to determine whether SLM or SEM is more suitable for the simulation of PM 2.5 , a Lagrange multiplier (LM) test and robust Lagrange multiplier (RLM) test should be estimated by the OLS. Anselin et al. proposed the criterion that if SLM-LM and SEM-LM are not significant, the OLS model was selected as the final model. If SLM-LM is significant and SEM-LM is not significant, SLM will be selected, and vice versa for SEM; if both SLM-LM and SEM-LM are insignificant, SLM-RLM is significant but SEM-RLM is not significant, SLM will be selected; if both SLM-LM and SEM-LM are insignificant, SEM-RLM is significant but SLM-RLM is not significant, SEM will be selected [46].
The univariate spatial autocorrelation analysis, the bivariate spatial correlation analysis and the spatial regression analysis were conducted in GeoDa software (GeoDa Press LLC, Chicago, IL, USA), and we chose queen contiguity weight matrix in GeoDa software.

The Spatial Distribution of Socioeconomic Factors and PM 2.5 in China
From Figure 1, 31 provinces in China were classified and mapped according to the values of PM 2.5 and socioeconomic factors. From Figure 1a, it could be found that PM 2.5 concentration only in Tibet met the WHO Air Quality Guideline (AQG) level (10 µg/m 3 ) in 1998. PM 2.5 concentration in most provinces were observed between 10 µg/m 3 and 35 µg/m 3 . Furthermore, PM 2.5 concentration in some provinces, such as Tianjin, Anhui, Shandong, Gansu, Ningxia, and Xinjiang, were found to be higher than 35 µg/m 3 . From Figure 1b, in 2016, obvious changes mainly occurred in some provinces of China. For example, PM 2.5 concentration increased obviously (>35 µg/m 3 ) in Liaoning, Beijing Hebei, Jiangsu, Shanghai and Henan; however, PM 2.5 concentrations in Gansu and Ningxia were found to have fallen below 35 µg/m 3 . The distributions of socioeconomic factors were similar to the distribution of PM 2.5 both in 1998 and 2016 generally. From Figure 1c,d, provinces with GDPP below 10,000 yuan accounted for more than 80% in 1998. Obviously, GDPP in all provinces was higher than 10,000 yuan in 2016, some of which had a GDPP of more than 100,000 yuan, such as Shanghai, Beijing and Tianjin.  Figure 1g). In 2016, the UPD in most provinces was higher than 2000 person/per square kilometer, and some provinces even exceeded 3000 person/per square kilometer (Figure 1h). From Figure 1i,j, we could find that PCO in Mainland China was less than 1 million in 1998. PCO in 31 provinces was more than 1 million, except Tibet and Qinghai in 2016. spatial regression analysis were conducted in GeoDa software (GeoDa Press LLC, Chicago, IL, USA), and we chose queen contiguity weight matrix in GeoDa software.

The Spatial Distribution of Socioeconomic Factors and PM2.5 in China
From Figure 1, 31 provinces in China were classified and mapped according to the values of PM2.5 and socioeconomic factors. From Figure 1a, it could be found that PM2.5 concentration only in Tibet met the WHO Air Quality Guideline (AQG) level (10μg/m 3 ) in 1998. PM2.5 concentration in most provinces were observed between 10μg/m 3 and 35μg/m 3 . Furthermore, PM2.5 concentration in some provinces, such as Tianjin, Anhui, Shandong, Gansu, Ningxia, and Xinjiang, were found to be higher than 35 μg/m 3 . From Figure 1b, in 2016, obvious changes mainly occurred in some provinces of China. For example, PM2.5 concentration increased obviously (>35 μg/m 3 ) in Liaoning, Beijing Hebei, Jiangsu, Shanghai and Henan; however, PM2.5 concentrations in Gansu and Ningxia were found to have fallen below 35μg/m 3 . The distributions of socioeconomic factors were similar to the distribution of PM2.5 both in 1998 and 2016 generally. From Figure 1c and Figure 1d, provinces with GDPP below 10,000 yuan accounted for more than 80% in 1998. Obviously, GDPP in all provinces was higher than 10,000 yuan in 2016, some of which had a GDPP of more than 100,000 yuan, such as Shanghai, Beijing Figure 1g). In 2016, the UPD in most provinces was higher than 2000 person/per square kilometer, and some provinces even exceeded 3000 person/per square kilometer (Figure 1h). From Figure 1I,j , we could find that PCO in Mainland China was less than 1 million in 1998. PCO in 31 provinces was more than 1 million, except Tibet and Qinghai in 2016.    Table 1, and Table 2 showed the temporal variation trend (the fitted slope) of 1998-2016. In 1998, PM2.5, GDPP, IVA, UPD and PCO in Mainland China were 23.97μg/m 3 , 6860 yuan/person, 3413.49 billion, 459 person/per square kilometer and 4.24 million private cars, respectively; and reached to 29.68μg/m 3 , 53,935 yuan/person, 24787.78 billion, 2408 person/per square kilometer and 163.3 million private cars in 2016, respectively. The fitted slope of PM2.5, GDPP, IVA, UPD and PCO in Mainland China were 0.138, 0.173, 0.174, 0.164 and 0.165, respectively. These indicated that GDPP, IVA, UPD, PCO and PM2.5 generally increased from 1998 to 2016, but the increased trend of IVA, GDPP, PCO and UPD was faster than the increased trend of PM2.5. Figure 2 showed their temporal variations intuitively. It could be found that GDPP, IVA, UPD and PCO in Mainland China showed an increase trend gradually in 1998-2016. The PM2.5 concentration also increased generally but began to fluctuate sharply from 2010. It indicated that the increasing trend of PM2.5 concentration was similar to that of GDPP, IVA, UPD and PCO in 1998-2016, and this increasing trend was significant, especially before 2006.   Table 1, and Table 2 Figure 2 showed their temporal variations intuitively. It could be found that GDPP, IVA, UPD and PCO in Mainland China showed an increase trend gradually in 1998-2016. The PM 2.5 concentration also increased generally but began to fluctuate sharply from 2010. It indicated that the increasing trend of PM 2.5 concentration was similar to that of GDPP, IVA, UPD and PCO in 1998-2016, and this increasing trend was significant, especially before 2006.

The Temporal Variation of Socioeconomic Factors and PM 2.5 in the Seven Geographical Subareas
The seven regions of China are showed in Figure 3. From East China. To better analyze the variation of PM 2.5 (GDPP, IVA, UPD or PCO) in seven sub-regions, this slopegraph in Figure 4 can be used to show the increases/decreases between just two fixed points (1998 and 2016) for different factors. Most importantly, slopegraph focused on the overall macro change between two periods points, not changes in each year or intervening period. Slopegraph is a great visualization method for focusing on that aspect of the macro change.

The Temporal Variation of Socioeconomic Factors and PM2.5 in the Seven Geographical Subareas
The seven regions of China are showed in Figure 3. From Table 1, in 1998, the GDPP and UPD in East China and North China, the IVA in East China, and the PCO in North China and Central China were far higher than other geographical subareas. The GDPP, IVA, UPD, and PCO in Northwest China were relatively lower. However, Northwest China had the highest PM2.5 From Table 2 and Figure 4, it could be found that PM 2.5 concentration in subareas except Northwest China and Southwest China presented an obviously increasing trend, and PM 2.5 concentration in Southwest China increased slowly over the years. However, PM 2.5 concentration in Northwest China presented a descending trend. Some literatures suggested that sand and dust was the major cause of affecting PM 2.5 concentration in Northwest China [38]. The possible reason on the minus slope (−0.015) for PM 2.5 in Northwest China may be an increase in vegetation coverage [49], and the decrease of dust events in Northern China in recent decades. The reduction of the wind speed in the northern hemisphere was the main reason for the decrease of dust event incidence [50]. The socioeconomic factors in the seven geographical subareas all presented an increasing trend, likely leading to the increase of PM 2.5 concentration between 1998 and 2016. changes in each year or intervening period. Slopegraph is a great visualization method for focusing on that aspect of the macro change.
From Table 2 and Figure 4, it could be found that PM2.5 concentration in subareas except Northwest China and Southwest China presented an obviously increasing trend, and PM2.5 concentration in Southwest China increased slowly over the years. However, PM2.5 concentration in Northwest China presented a descending trend. Some literatures suggested that sand and dust was the major cause of affecting PM2.5 concentration in Northwest China [38]. The possible reason on the minus slope (−0.015) for PM2.5 in Northwest China may be an increase in vegetation coverage [49], and the decrease of dust events in Northern China in recent decades. The reduction of the wind speed in the northern hemisphere was the main reason for the decrease of dust event incidence [50]. The socioeconomic factors in the seven geographical subareas all presented an increasing trend, likely leading to the increase of PM2.5 concentration between 1998 and 2016.
.  changes in each year or intervening period. Slopegraph is a great visualization method for focusing on that aspect of the macro change.
From Table 2 and Figure 4, it could be found that PM2.5 concentration in subareas except Northwest China and Southwest China presented an obviously increasing trend, and PM2.5 concentration in Southwest China increased slowly over the years. However, PM2.5 concentration in Northwest China presented a descending trend. Some literatures suggested that sand and dust was the major cause of affecting PM2.5 concentration in Northwest China [38]. The possible reason on the minus slope (−0.015) for PM2.5 in Northwest China may be an increase in vegetation coverage [49], and the decrease of dust events in Northern China in recent decades. The reduction of the wind speed in the northern hemisphere was the main reason for the decrease of dust event incidence [50]. The socioeconomic factors in the seven geographical subareas all presented an increasing trend, likely leading to the increase of PM2.5 concentration between 1998 and 2016. .

The Spatial Distribution of Temporal Trends for Socioeconomic Factors and PM2.5 in Different Provinces
The slope values of different provinces were mapped in Figure 5. From Figure 5a, it could be found that PM2.5 in most provinces of China increased rapidly. The provinces with a slower growth in PM2.5 were mainly distributed in Inner Mongolia, Sichuan, Chongqing, Guizhou and Yunnan. On the contrary, PM2.5 of Gansu, Ningxia, and Shaanxi presented showed a downward trend. From Figures 5b-e, we can see that the fitted slopes of GDPP, IVA and PCO in different provinces were all more than 0.155, indicating that the increased trends of GDPP, IVA and PCO were rapid in provinces in 1998-2016. UPD increased rapidly in most provinces except Beijing, Ningxia, Jiangsu and Hainan, and the fitted slope of UPD only in Beijing was negative, this may be because of Beijing's population control policies.  The slope values of different provinces were mapped in Figure 5. From Figure 5a, it could be found that PM 2.5 in most provinces of China increased rapidly. The provinces with a slower growth in PM 2.5 were mainly distributed in Inner Mongolia, Sichuan, Chongqing, Guizhou and Yunnan. On the contrary, PM 2.5 of Gansu, Ningxia, and Shaanxi presented showed a downward trend. From Figure 5b-e, we can see that the fitted slopes of GDPP, IVA and PCO in different provinces were all more than 0.155, indicating that the increased trends of GDPP, IVA and PCO were rapid in provinces in 1998-2016. UPD increased rapidly in most provinces except Beijing, Ningxia, Jiangsu and Hainan, and the fitted slope of UPD only in Beijing was negative, this may be because of Beijing's population control policies.

The Spatial Distribution of Temporal Trends for Socioeconomic Factors and PM2.5 in Different Provinces
The slope values of different provinces were mapped in Figure 5. From Figure 5a, it could be found that PM2.5 in most provinces of China increased rapidly. The provinces with a slower growth in PM2.5 were mainly distributed in Inner Mongolia, Sichuan, Chongqing, Guizhou and Yunnan. On the contrary, PM2.5 of Gansu, Ningxia, and Shaanxi presented showed a downward trend. From Figures 5b-e, we can see that the fitted slopes of GDPP, IVA and PCO in different provinces were all more than 0.155, indicating that the increased trends of GDPP, IVA and PCO were rapid in provinces in 1998-2016. UPD increased rapidly in most provinces except Beijing, Ningxia, Jiangsu and Hainan, and the fitted slope of UPD only in Beijing was negative, this may be because of Beijing's population control policies.

The Correlation between Socioeconomic Factors and PM 2.5 in Mainland China
The Spearman's rank correlation coefficients (r-GDPP, r-IVA, r-UPD and r-PCO) between PM 2.5 concentration and GDPP, IVA, UPD, and PCO in Mainland China in 1998-2016 were shown in Figure 6. Most of the correlation coefficients in Figure 6 were positive, indicating that PM 2.5 was positively correlated with socioeconomic factors. From Figure 6a The Spearman's rank correlation coefficients (r-GDPP, r-IVA, r-UPD and r-PCO) between PM2.5 concentration and GDPP, IVA, UPD, and PCO in Mainland China in 1998-2016 were shown in Figure 6. Most of the correlation coefficients in Figure 6 were positive, indicating that PM2.5 was positively correlated with socioeconomic factors. From Figures 6a,b, Figure 6c,d, most of the correlation coefficients were positive except for a few years. All the p-values were higher than 0.05, indicating that PM2.5 had a positively correlation with UPD and PCO in most years, but the correlations were not significant during the research period. The weak correlation between PM2.5 and PCO increased obviously before 2003. However, the correlation coefficients between PM2.5 and UPD presented a downward trend since 2001, indicating that the impact of UPD on PM2.5 was getting weaker and weaker. In the Figure 7, dark green represents a significant negative correlation, light green is a negative correlation, dark yellow means a significant positive correlation, and pale yellow represents a positive correlation. From    This was consistent with anthropogenic effects on the dust loading in East China was far higher than near desert source regions in Northwest China [50]. There were other factors which determines the PM 2.5 trend in Northwest China, PM 2.5 in Northwest China was mainly affected by sand and dust [38].
Previous studies have shown a positive correlation between air temperature and PM 2.5 concentration in summer in Northwest China [51]. affected by sand and dust [38]. Previous studies have shown a positive correlation between air temperature and PM2.5 concentration in summer in Northwest China [51].

The Spatial Statistical Relationship between Socioeconomic Factors and PM2.5
3.4.1. Global Spatial Autocorrelation of PM2.5 From Figure 9, the global Moran's I values of PM2.5 were positive at the 95% confidence level and increased over time, but fluctuated around 0.5 since 2003, indicating that PM2.5 exibited significantly positive spatial autocorrelation and spatial homogeneous, and spatial autocorrelation of PM2.5 in 31 provinces of China strengthened gradually. In other words, PM2.5 at one province tended to be similar to those of their neighboring provinces, the spatial spillover effect had been increasing in different provinces.

Global Spatial Autocorrelation of PM 2.5
From Figure 9, the global Moran's I values of PM 2.5 were positive at the 95% confidence level and increased over time, but fluctuated around 0.5 since 2003, indicating that PM 2.5 exibited significantly positive spatial autocorrelation and spatial homogeneous, and spatial autocorrelation of PM 2.5 in 31 provinces of China strengthened gradually. In other words, PM 2.5 at one province tended to be similar to those of their neighboring provinces, the spatial spillover effect had been increasing in different provinces. affected by sand and dust [38]. Previous studies have shown a positive correlation between air temperature and PM2.5 concentration in summer in Northwest China [51].

The Spatial Statistical Relationship between Socioeconomic Factors and PM2.5
3.4.1. Global Spatial Autocorrelation of PM2.5 From Figure 9, the global Moran's I values of PM2.5 were positive at the 95% confidence level and increased over time, but fluctuated around 0.5 since 2003, indicating that PM2.5 exibited significantly positive spatial autocorrelation and spatial homogeneous, and spatial autocorrelation of PM2.5 in 31 provinces of China strengthened gradually. In other words, PM2.5 at one province tended to be similar to those of their neighboring provinces, the spatial spillover effect had been increasing in different provinces. In order to identify the provinces with significant spatial correlation and type of spatial clusters of PM 2.5 , LISA was calculated and mapped in Figure 10. High-high (HH) clusters means that PM 2.5 concentration of one province and its neighbors were higher than the annual average PM 2.5 concentration in Mainland China. While, low-low (LL) clusters refers to the provinces with low PM 2.5 concentration being surrounded by neighbors with low PM 2.5 concentration, whose value is lower than the annual average values. High-low (HL) outliers means that high PM 2.5 concentration had low PM 2.5 concentration in the neighboring provinces and vice versa for the low-high (LH) outliers. The HH and LL clusters can reflect the similar PM 2.5 concentration clustering, indicating spatial autocorrelation of PM 2.5 is positive; spatial dispersion of PM 2.5 concentration is reflected in the HL and LH outliers, it indicates that PM 2.5 concentrations have a negative spatial autocorrelation. From Figure 10 In order to identify the provinces with significant spatial correlation and type of spatial clusters of PM2.5, LISA was calculated and mapped in Figure 10. High-high (HH) clusters means that PM2.5 concentration of one province and its neighbors were higher than the annual average PM2.5 concentration in Mainland China. While, low-low (LL) clusters refers to the provinces with low PM2.5 concentration being surrounded by neighbors with low PM2.5 concentration, whose value is lower than the annual average values. High-low (HL) outliers means that high PM2.5 concentration had low PM2.5 concentration in the neighboring provinces and vice versa for the low-high (LH) outliers. The HH and LL clusters can reflect the similar PM2.5 concentration clustering, indicating spatial autocorrelation of PM2.5 is positive; spatial dispersion of PM2.5 concentration is reflected in the HL and LH outliers, it indicates that PM2.5 concentrations have a negative spatial autocorrelation. From Figure 10  In this paper, global bivariate Moran's I was used to determine if PM2.5 in one province were spatial correlated with socioeconomic factors of its neighbors across the study region. From Figures

Spatial Correlations between PM 2.5 and Socioeconomic Factors
In this paper, global bivariate Moran's I was used to determine if PM 2.5 in one province were spatial correlated with socioeconomic factors of its neighbors across the study region. From Figure 11a,b,d, the global bivariate Moran's I values presented positive growth trends at a 95% confidence level in 1998-2016. These indicated that spatial correlations between PM 2.5 at a province and GDPP, IVA and PCO of its adjacent provinces were positive and significant, and the positive spatial correlations increased during the study period. From Figure 11c As shown in Figure 12, the bivariate local Moran's I values for PM 2.5 and socioeconomic factors were calculated. HH clusters means that the provinces with high PM 2.5 concentration clustered the neighboring provinces with high values of GDPP, IVA, UPD and PCO, and their values were higher than their annual average values. LL clusters means that the provinces with low PM 2.5 concentration were near predominantly the provinces with low values of GDPP, IVA, UPD and PCO, and their values were lower than their annual average values. HL outliers occur where the neighbors of the provinces with high PM 2.5 concentration have low GDPP, IVA, UPD and PCO. LH outliers mean that there were low values of PM 2.5 concentration in one province, and there were high values of GDPP, IVA, UPD and PCO in the adjacent provinces. We used data from 2016 as an example to analyze the local spatial correlations between PM 2.5 concentration and socioeconomic factors. In 2016, from Figure 12c, Hubei was the only province where appeared HH clusters of PM 2.5 concentration and UPD. From Figure 12a,b,d, the LL clusters of PM 2.5 concentration and GDPP (IVA or PCO) were mostly covered in some provinces of Southwest China. Shanghai appeared a HH clusters of PM 2.5 concentration and GDPP. Shandong, Jiangsu, Shanghai and Anhui were the provinces that had a HH clusters of PM 2.5 concentration and IVA. High PM 2.5 concentration and high PCO clustered in Henan, Shandong, Jiangsu, Shanghai and Anhui. The place where was the HL outliers of PM 2.5 concentration and GDPP (IVA or PCO) was Xinjiang. The LH outliers of PM 2.5 concentration and IVA were covered in Fujian and Jiangxi. Fujian was the province that appeared a LH outliers of PM 2.5 concentration and PCO. In this paper, global bivariate Moran's I was used to determine if PM2.5 in one province were spatial correlated with socioeconomic factors of its neighbors across the study region. From Figures

Regression Results of the Spatial Regression Model
The analysis results of spatial autocorrelation confirmed the existence of spatial dependence of PM2.5, so spatial regression models were used to further confirm the spatial dependence of PM2.5 concentration. First, the estimated results of OLS were calculated in Table 3. LM test and RLM test were performed for the residuals of OLS regression. The values of SLM-LM and SEM-LM were significant (p < 0.05) except for 1999 and 2016. The values of SLM-RLM were significant (p < 0.1) in most years, while SEM-RLM was not significant in 1998-2016. Therefore, the SLM model was adopted. The results of spatial lag model regression in 1998-2016 were shown in Table 4

Regression Results of the Spatial Regression Model
The analysis results of spatial autocorrelation confirmed the existence of spatial dependence of PM 2.5 , so spatial regression models were used to further confirm the spatial dependence of PM 2.5 concentration. First, the estimated results of OLS were calculated in Table 3. LM test and RLM test were performed for the residuals of OLS regression. The values of SLM-LM and SEM-LM were significant (p < 0.05) except for 1999 and 2016. The values of SLM-RLM were significant (p < 0.1) in most years, while SEM-RLM was not significant in 1998-2016. Therefore, the SLM model was adopted. The results of spatial lag model regression in 1998-2016 were shown in Table 4  Notes: Log-L denotes Log likelihood; *, **, ***represent coefficients are significant at the 10%, 5%, 1% levels, respectively. Notes: Log-L denotes log likelihood; *, **, ***represent coefficients significant at the 10%, 5%, 1% levels, respectively.

Spatial Distribution and Temporal Variation of PM 2.5 and Socioeconomic Factors
Based on four socioeconomic factors dataset and PM 2.5 concentration dataset, this study examined the spatial distribution and relationships between socioeconomic factors and PM 2. Previous studies have also shown that high PM 2.5 concentration were mainly distributed in economically developed areas [38]. From Figure 2, the temporal variations showed that the overall increase trend of PM 2.5 is the same as that of GDPP, IVA, UPD and PCO during 1998-2016, but PM 2.5 exhibited a downward trend from 2006 to 2012. The reason for this phenomenon may be the implementation of sustainable development policies of energy conservation, pollutant reduction and green development proposed in the eleventh five-year plan [52]. The external cause could be meteorological factors. For example, the nitrate and secondary organic aerosols formation was greatly facilitated by high humidity [53]. Wind speed is conducive to the diffusion of PM 2.5 [54]. The chemical reaction rate of PM 2.5 precursor pollutants accelerates with the increase of temperature and solar radiation. [55,56]. From Tables 1 and 2 [19,50]. From Figure 5, a downward trend of PM 2.5 presented in Gansu, Ningxia, and Shaanxi is mainly attributed to implementation of clean air policies in recent years. UPD only in Beijing showed a downward trend, which is in line with the requirements of the Beijing-Tianjin-Hebei coordination to "strictly control the increase, dredge the stock, dredge the combination" of Beijing's population size. In addition to the policy factor of Beijing-Tianjin-Hebei cooperation, the negative growth of Beijing's permanent population is also related to the overall trend of population returns in labor-exporting provinces.

The Relationships between PM 2.5 and Socioeconomic Factors
From Figure 7, the Spearman's rank correlation analysis indicated that four socioeconomic factors produced an increase of PM 2.5 in most provinces of China. However, GDPP, IVA and PCO appeared a negative correlation with PM 2.5 in Shaanxi, Gansu and Ningxia. Gansu, Ningxia and Beijing were the negative correlation between PM 2.5 and UPD. From Figure 8, the socioeconomic factors had strong impact on PM 2.5 concentration in North China, Northeast China and East China, but in contrast less affected PM 2.5 concentration in Northwest China. It's worth noting that meteorological factors and urban fugitive dust also contributed to PM 2.5 concentration [57,58]. Soil and desert dust was the major cause of high Fe and K contents in urban fugitive dust in Northern China, and PM 2.5 was more affected by soil dust in northern China than in southern China [59]. Furthermore, coal combustion produced fugitive dust which increased PM 2.5 concentration. This influence was especially strong in Northern part of China [60]. In addition, it was reported that desert dust and soil dust often affected Northwest China. So, sand and dust played an important role in influencing PM 2.5 concentration in Northwest China [38,61,62]. GDPP and IVA appeared significantly positive correlations with PM 2.5 in most years in Figure 6a,b. While the correlation between PM 2.5 and UPD (or PCO) was all insignificant in Figure 6c,d. Furthermore, when spatial factors were considered in Figure 11a,b,d, GDPP (IVA or PCO) imposed a positive externality on PM 2.5 ; that is, the increase of GDPP (IVA or PCO) in one province may cause the increase of PM 2.5 in the neighboring provinces. The reason is that the pollution particles, generated by the activities of residents, emissions from factories and private cars, may be passed from one province to the surrounding provinces through atmospheric movements such as wind speed, wind direction temperature. High temperature and wind speed can promote the convection of air. This can create better conditions for the dilution and dispersion of particulate matter. Notably, the Spearman's rank correlation analysis and bivariate spatial correlation analysis gave a consistent conclusion for the downward trend on the UPD's impact on PM 2.5 concentration in Figures 6c and 11c. It indicated that the impact of UPD on PM 2.5 was getting smaller and smaller. This may be because of the population control policy. The population size of some provinces has been gradually controlled since the population control policy was implemented. The impact of UPD on PM 2.5 may be closely related to population size [63]. Figure 12 showed the local bivariate cluster maps for PM 2.5 and socioeconomic factors in 31 provinces of China, in 2016. For the provinces of the HH clusters, the development of the social economy in their adjacent provinces had positive radiation effect on these provinces. Their economic development in the local provinces have also brought about a number of pollution sources that have indirectly increased PM 2.5 . Some provinces with slow economic growth in Northwestern and Southwestern China had fewer pollution sources, which was easy to form LL clusters. Furthermore, the HL type provinces were mainly distributed in Xinjiang in Figures 10 and 12. As a region of severe sandstorm and abundant coal resources, the exploitation and utilization of these coal resources have produced many pollutants in Xinjiang and destroyed the ecological balance of atmospheric environment. The complex topography of Xinjiang is also not conducive to the diffusion of atmospheric pollutants. In addition, Qinghai and Tibet with underdeveloped industry have less pollution resources and lower PM 2.5 concentration. Rich vegetation in Sichuan and Yunnan can effectively reduce PM 2.5 concentration. The combination of these factors formed the obvious HL outliers around Xinjiang.

The Spatial Spillover Effect of PM 2.5
The spatial spillover effect means that the changes of PM 2.5 concentration in one province can impact on PM 2.5 concentration of other provinces. In this paper, spatial spillover effect of PM 2.5 concentration in adjacent provinces can be reflected by the global Moran's I. From Figure 9, there was a positive increasing trend of the global Moran's I values of PM 2.5 concentration during the study period. It indicated that the spatial correlation of PM 2.5 gradually became stronger over time. The spatial autoregressive coefficient (W*PM25) were all significant (p < 0.01) in 1998-2016 (column 2 of Table 4). These meant that the spatial spillover effect is becoming more and more significant. From Figures 10  and 12, the HH clusters of PM 2.5 concentration (HH clusters of PM 2.5 and socioeconomic factors) were mainly distributed in some provinces of economically developed area (i.e., North China, East China). At the same time, provinces of economically backward areas (i.e., Southwest China) appeared to have LL clusters. These indicated that the spatial spillover effect of North China, East China and Southwest China were higher than other regions. All the three regions have strong PM 2.5 pollution homogeneity. In other words, there were spatial spillover effects in different provinces, but were particularly severe in North China, East China and Southwest China. However, we note that the spatial spillover effects on PM 2.5 pollution for all regions are non-negligible. So local governments should consider the policies of adjacent provinces and coordination with adjacent provinces is indispensable. To further analyze the correlation between PM 2.5 and aerosol emission density, we used dataset of GDP per area and IVA per area to calculate the Spearman's rank correlation coefficients. Due to the absence of data in 2016 and Tibet, the time sequence of the experiment was from 1998 to 2015, and Tibet was excluded. We have tested the correlations from a spatial perspective. The experiment was designed for the effects of GDPP, GDP per area, IVA and IVA per area on PM 2.5 from two different scale, including regional scale and provincial scale. On regional scale, the correlation between PM 2.5 and GDP per area (IVA per area, GDPP, or IVA) were significant in North China, Northeast China, East China, Central China and South China in Figure 13. On provincial scale, Table 5 showed that most provinces of the other five geographical regions except Southwest China and Northwest China presented a significant correlation between PM 2.5 and GDP per area (IVA per area, GDPP, or IVA). Although there were some slight differences in the correlations values and p values under two different scales, the overall trend was consistent. These indicated that the increase of GDP and industry has a strong positive impact on PM 2.5 , especially in North China, Northeast China, East China, Central China and South China. However, the influence was not strong in Northwest China. The reason may be that PM 2.5 concentration in Northwest China is more affected by sandstorms. These further validated the idea in this article. That is, human activities contribute to PM 2.5 concentration, but are not the only factor.    Notes: *, ** represent coefficients are significant at the 5%, 1% levels, respectively.

Conclusions
This paper estimated spatial distribution, temporal variations and relationships of socioeconomic factors and PM 2.5 in 31 provinces of China using a unary linear regression model, Spearman's rank correlation analysis method, univariate spatial autocorrelation analysis method, bivariate spatial correlation analysis method and the spatial regression analysis during the period of 1998-2016. Results demonstrated that PM 2.5 generally increased with the increase of socioeconomic factors from 1998 to 2016, but there were different temporal variations trend and relationships in different provinces and regions. Socioeconomic factors and PM 2.5 concentration in most provinces in East China, South China, Central China, North China, and Northeast China had rapid growth trend, and socioeconomic factors were significantly correlated with PM 2.5 concentration. Although the growth trend of socioeconomic factors in Southwest China and Northwest China were also fast, PM 2.5 presented a slowly growth trend in Southwest China and a descending trend in Northwest China, and socioeconomic factors were weakly correlated with PM 2.5 concentration. Urban population density was not an important influencing factor in affecting PM 2.5 concentration. GDP per capita and industrial added values in the local and adjacent provinces were the key influencing factors for the increase of PM 2.5 concentration. Private car ownership also contributed to PM 2.5 concentration. PM 2.5 in neighboring provinces were also an important factor to increase the local PM 2.5 concentration. The results of the research can provide effective guidelines for urban sustainable development and further protect the environment of cities.
Author Contributions: J.L. conceived and designed the study; Y.Y. analyzed the data and wrote the paper. J.L., Y.Y., G.Z and Q.Y. contribute to the editing and reviewing of the paper.