The Effect of Socioeconomic Factors on Spatiotemporal Patterns of PM2.5 Concentration in Beijing–Tianjin–Hebei Region and Surrounding Areas

The study investigated the spatiotemporal evolution of PM2.5 concentration in the Beijing–Tianjin–Hebei region and surrounding areas during 2015–2017, and then analyzed its socioeconomic determinants. First, an estimation model considering spatiotemporal heterogeneous relationships was developed to accurately estimate the spatial distribution of PM2.5 concentration. Additionally, socioeconomic determinants of PM2.5 concentration were analyzed using a spatial panel Dubin model, which aimed to improve the robustness of the model estimation. The results demonstrated that: (1) The proposed model significantly increased the estimation accuracy of PM2.5 concentration. The mean absolute error and root-mean-square error were 9.21 μg/m3 and 13.10 μg/m3, respectively. (2) PM2.5 concentration in the study area exhibited significant spatiotemporal changes. Although the PM2.5 concentration has declined year by year, it still exceeded national environmental air quality standards. (3) The per capita GDP, urbanization rate and number of industrial enterprises above the designated size were the key factors affecting the spatiotemporal distribution of PM2.5 concentration. This study provided scientific references for comprehensive PM2.5 pollution control in the study area.


Introduction
Atmospheric pollution significantly influences human health, climatic environment, and sustainable urban development [1][2][3][4]. According to a World Health Organization (WHO) report in 2014, atmospheric pollution causes more than seven million deaths worldwide each year [5]. A recent study based on the Global Exposure Mortality Model estimated that 8.9 million people died globally in 2015 [6]. With China's rapid economic, industrial, and urban development, atmospheric pollution has become an increasing problem. East China, especially the Beijing-Tianjin-Hebei region, has witnessed frequent occurrences of severe haze since 2013. The Beijing-Tianjin-Hebei region is responsible for the majority of North China's economic development, and coordinated development in this region is one of three national strategies in China [7]. Therefore, comprehensive atmospheric pollution governance in the region has attracted wide attention from China's government. PM 2.5 (atmospheric particulate energy consumption, industrial structure, population, and environmental background. Ignoring these differences is very likely to cause biased or suspicious conclusions.
Given that the national observation network was not completed until the end of 2014, this study focused on the spatiotemporal evolution of PM 2.5 concentration in 28 cities of the Beijing-Tianjin-Hebei region and surrounding areas during 2015-2017, and then identified its socioeconomic determinants. Firstly, an estimation model for high spatial-resolution PM 2.5 estimation was created based on the reconstructed AOD missing gaps and the spatiotemporal heterogeneous relationship between PM 2.5 and auxiliary variables, which disclosed the spatial distribution of PM 2.5 concentration in the study area. Secondly, based on the analysis of the spatiotemporal evolution of PM 2.5 concentration, the socioeconomic factors that influence local PM 2.5 concentration were investigated by a spatial panel model. The research study's conclusions provide scientific references for local atmospheric pollution control.

Study Area
The study area includes 28 cities of the Beijing-Tianjin-Hebei region and surrounding areas, covering an area of~275,000 km 2 . As shown in Figure 1, the terrain is high in the west and low in the east, with elevation ranging from sea level to >2000 m. The region experiences four distinct seasons, with hot and rainy summers due to the East Asian monsoon, and cold and dry winters due to subtropical high-pressure systems. The Beijing-Tianjin-Hebei region is the most economically developed area in northern China and is the area with the most PM 2.5 pollution. and environmental background. Ignoring these differences is very likely to cause biased or suspicious conclusions. Given that the national observation network was not completed until the end of 2014, this study focused on the spatiotemporal evolution of PM2.5 concentration in 28 cities of the Beijing-Tianjin-Hebei region and surrounding areas during 2015-2017, and then identified its socioeconomic determinants. Firstly, an estimation model for high spatial-resolution PM2.5 estimation was created based on the reconstructed AOD missing gaps and the spatiotemporal heterogeneous relationship between PM2.5 and auxiliary variables, which disclosed the spatial distribution of PM2.5 concentration in the study area. Secondly, based on the analysis of the spatiotemporal evolution of PM2.5 concentration, the socioeconomic factors that influence local PM2.5 concentration were investigated by a spatial panel model. The research study's conclusions provide scientific references for local atmospheric pollution control.

Study Area
The study area includes 28 cities of the Beijing-Tianjin-Hebei region and surrounding areas, covering an area of ~275,000 km 2 . As shown in Figure 1, the terrain is high in the west and low in the east, with elevation ranging from sea level to >2000 m. The region experiences four distinct seasons, with hot and rainy summers due to the East Asian monsoon, and cold and dry winters due to subtropical high-pressure systems. The Beijing-Tianjin-Hebei region is the most economically developed area in northern China and is the area with the most PM2.5 pollution.

Data
Ground PM2.5 observations were from the public platform of the China National Environmental Monitoring Center (http://www.cnemc.cn). Before release, these data had been calibrated and quality controlled to meet the national environmental air quality standards of China (GB3095-2012). The present study used daily average PM2.5 concentration (DAPC) data from January 2015 to December 2017 at 256 observation stations in the study area and surrounding areas. Data of monthly average

Data
Ground PM2.5 observations were from the public platform of the China National Environmental Monitoring Center (http://www.cnemc.cn). Before release, these data had been calibrated and quality controlled to meet the national environmental air quality standards of China (GB3095-2012). The present study used daily average PM 2.5 concentration (DAPC) data from January 2015 to December 2017 at 256 observation stations in the study area and surrounding areas. Data of monthly average PM 2.5 concentrations (MAPC) were generated by averaging the DAPC. Spatial distribution of observation stations is shown in Figure 1. Among 256 observation stations, 83 were used for constructing the PM 2.5 estimation model, and the remaining were used to verify the accuracy of the estimation results.
The latest C6 version of daily AOD (DAOD) was used to construct PM 2.5 concentration estimation model. The C6 version DAOD had higher spatial resolutions (3 km) compared with the previous C5 version. Data of DAOD of Aqua (DAODA) and DAOD of Terra (DAODT) from January 2015 to December 2017 were collected from the website https://ladsweb.modaps.eosdis.nasa.gov/.
According to practical situations in the study area, air temperature (AT), wind speed (WS) at 10 m, surface pressure (SP), and boundary layer height (BLH) were chosen to assist the estimated spatial distribution of PM 2.5 concentration (Table 1). These data came from the ERA-Interim reanalysis data (http://apps.ecmwf.int/datasets/) of the European Centre for Medium-Range Weather Forecast, with a spatial resolution of 0.125 • . Socioeconomic data were collected from the China Statistics Yearbook, China City Statistics Yearbook and statistical yearbooks of provinces and regions in the study area (http://data.cnki.net/ Yearbook). Six factors were chosen: person density (PD), per capita gross regional product (PGRP), urbanization rate (UR), the proportion of secondary industry in GDP (PSIGDP), industrial smoke (dust) emissions (ISDE), and the number of industrial enterprises above designated size (NIEDS) ( Table 1). All PGRP data were transformed uniformly to a constant price in 2015. Additionally, logarithmic transformations were performed on all socioeconomic data to eliminate their heteroscedasticity.

PM 2.5 Estimation
Due to the impact of clouds, rain, and other weather conditions, there are a lot of gaps (no data region) in DAOD data. Therefore, it was necessary to fill the missing data gaps. DAOD data were constructed based on the complementarity between DAODT and DAODA on spatial coverage and significant correlation. The data were filled as follows: DAODT and DAODA of a month were used to establish the relationship: where m refers to month. a T,m , b T,m , a A,m and b A,m are regression coefficients between DAODT and DAODA. ε T and ε A are error terms. Next, the missing data gaps of DAODT and DAODA were reconstructed based on the acquired relationships in Equations (1) and (2), generating reconstructed DAODT (RDAODT) and reconstructed DAODA (RDAODA). For example, if DAODT has a value V test at location L test , but DAODA does not, we can use V test and Equation (2) to estimate the value of DAODA at location L test . In this way, DAOD was estimated as (RDAODT + RDAODA)/2 because Aqua and Terra measure AOD in the morning and afternoon, respectively. Finally, the monthly DAOD averages were calculated, which were used to generate the monthly average AOD data (MAOD).
Although AOD is an indicator of PM 2.5 concentration, PM 2.5 concentration is also significantly influenced by air temperature, precipitation, and other climatic factors [34,35]. This study used MAOD, AT, WS, BLH, and SP as the auxiliary variables to estimate the spatial distribution of PM 2.5 concentration. We assumed that the relationship between PM 2.5 concentration and auxiliary variables changes with time and space and the following model was constructed: whereMAPC m is the estimated average PM 2.5 concentration during the month m. u refers to spatial position. α m is intercept, and β m,0 and β m,i are coefficients of MAOD m and other auxiliary variables AUX m,i . n is a variable with a value range of < =4. When n = 0, no climatic factor is chosen. When n = 4, the AT, WS, BLH, and SP were all used to construct the model. With regards to the temporal heterogeneous relationship between PM 2.5 and auxiliary variables, the auxiliary variables of the model were chosen based on the following criteria with consideration to temporal changes of the relationship between PM 2.5 and auxiliary variables: (1) the chosen auxiliary variables were significantly correlated with PM 2.5 ; (2) the chosen auxiliary variables improved the interpretation of the model to PM 2.5 variation. In the views of the spatial heterogeneous relationship between PM 2.5 and auxiliary variables, a local regression method, geographically weighted regression (GWR) [36], was applied to assess and describe the relationship.
For each month, observations of the training stations were used to construct the PM 2.5 estimation model, and observations from the validation stations were used to validate the accuracy of the estimated result ( Figure 1). Some statistical indexes, including correlation coefficient, mean absolute error (MAE), and root-mean-square error (RMSE), were chosen to evaluate the effectiveness of the proposed model.

Effect of Economic and Social Factors on PM 2.5 Concentration
Socioeconomic data from the statistical yearbook were based on city-scale annual statistics. Hence, the derived MAPC data should be processed accordingly, which generated the city-scale annual average PM 2.5 concentration (AAPC). The logarithms of AAPC were calculated to ensure consistency with the pre-processing of Socioeconomic data.
The PM 2.5 distribution presented strong trans-regional characteristics and inevitably affected nearby regions. Global Morans' I analysis [37] was used to measure the spatial correlation of PM 2.5 concentrations. In addition, local Morans' I analysis [37] was applied to describe the spatial heterogeneity of PM 2.5 concentrations in different geographical units.
The effects of socioeconomic factors on PM 2.5 concentration in the urban scale were analyzed by a spatial panel model [38,39]: where i refers to city and t is the year. y it is the explained variable, which is equal to AAPC of city i in year t. lnX it is the explanatory variable which refers to socioeconomic factors, and β is the corresponding coefficients. µ i is the spatial effect and η t is the time-period effect. w i,j refers to elements in the spatial weight matrix W. ρ is the spatial autoregression coefficient of dependent variables. γ denotes the spatial autocorrelation vector of explanatory variables. λ is the spatial autocorrelation coefficient of the error term.

Construction of the Estimation Model
Auxiliary variables were chosen monthly according to the selection criteria in Section 3.1 ( Table 2). The estimation models of all months involved MAOD, which again confirms that AOD was a good indicator of PM 2.5 concentration. The number of other chosen auxiliary variables changes with time, indicating that although PM 2.5 concentration was greatly affected by the climatic conditions, there was a significant temporal change in sensitivity. This validates the justifiability of the proposed assumption. Many studies [6,40] have reported that precipitation affects PM 2.5 concentration, and the effect is more significant in the time dimension or in the large spatial range. However, in this study, we constructed an estimation model for each month, which reduces the effect in the time dimension. Next, unlike other statistical methods, GWR is a local spatial regression method-only the data within the local range participates in the model construction, thereby weakening the effect in the large spatial range. As a result, precipitation was excluded in this study. This is consistent with other studies [41][42][43] that build PM 2.5 estimation models based on GWR.
The proposed spatiotemporal heterogeneous model (SHM) was compared with the uniform relationship model (UM) based on multiple linear regression. The construction accuracy of the estimation model throughout the study period is shown in Figure 2. UM showed a relatively lower goodness of fit and high temporal fluctuation, with minimum and maximum values of R 2 being 0.17 and 0.67, respectively. Comparatively, SHM increased interpretation to changes in PM 2.5 concentration. Its average of R 2 (0.77) was significantly higher than that of UM (0.45), while its average of RMSE (8.87 µg/m 3 ) was considerably smaller than that of UM (13.81 µg/m 3 ). In addition, SHM demonstrated better stability. All these indicate that it is necessary to consider the spatial heterogeneity of the relationships between PM 2.5 and auxiliary variables.
The proposed spatiotemporal heterogeneous model (SHM) was compared with the uniform relationship model (UM) based on multiple linear regression. The construction accuracy of the estimation model throughout the study period is shown in Figure 2. UM showed a relatively lower goodness of fit and high temporal fluctuation, with minimum and maximum values of R 2 being 0.17 and 0.67, respectively. Comparatively, SHM increased interpretation to changes in PM2.5 concentration. Its average of R 2 (0.77) was significantly higher than that of UM (0.45), while its average of RMSE (8.87 g/m 3 ) was considerably smaller than that of UM (13.81 g/m 3 ). In addition, SHM demonstrated better stability. All these indicate that it is necessary to consider the spatial heterogeneity of the relationships between PM2.5 and auxiliary variables.

Accuracy Validation and Estimation Results
The validation of estimated MAPC in the study area during 2015-2017 is shown in Figure 3. R, MAE and RMSE of UM were 0.89, 11.25 g/m 3 , and 15.55 g/m 3 , respectively. In contrast, SHM significantly increased the estimation accuracy of MAPC, increasing the correlation coefficient by 3% and decreasing MAE and RMSE by as much as 18% and 16%, respectively. As expected, the proposed model achieved a higher estimation accuracy of AAPC than UM. The correlation coefficient of the

Accuracy Validation and Estimation Results
The validation of estimated MAPC in the study area during 2015-2017 is shown in Figure 3. R, MAE and RMSE of UM were 0.89, 11.25 µg/m 3 , and 15.55 µg/m 3 , respectively. In contrast, SHM significantly increased the estimation accuracy of MAPC, increasing the correlation coefficient by 3% and decreasing MAE and RMSE by as much as 18% and 16%, respectively. As expected, the proposed model achieved a higher estimation accuracy of AAPC than UM. The correlation coefficient of the proposed model increased by 5%, while the MAE and RMSE decreased by as much as 19% and 17%, respectively ( Figure 4). Huang et al. [44] Figure 5. There is a significant spatiotemporal variation in PM2.5 concentration. In terms of spatial variation, high MAPC in the southeast region of the study area from January-April was observed, whereas PM2.5 concentration in the northwest region towards the central area was relatively higher from May-December. With regards to temporal variation, PM2.5 concentrations in January, February, November, and December were significantly higher than those in other months. The reason behind this may be related to indoor heating and climatic conditions. Some studies have reported that the burning of biomass and fossil energy for heating in winter generated huge PM2.5 emissions [47,48]. Dust storms, which frequently occur in North China in late winter and spring, is another major contribution that aggravates PM2.5 concentrations [49]. The lowest average MAPC (40.33 g/m 3 ) occurred in August, and the highest (138.49 g/m 3 ) was in December. The AAPC in the study area during 2015-2017 is shown in Figure 6. Generally, AAPC decreased from the central areas of Shijiazhuang, Baoding, Hengshui, and Xingtai to surrounding areas, accompanied with obvious concentration characteristics. The average APPC decreased from 77.3 g/m 3 in 2015 to 64.85 g/m 3 in 2017.     Figure 5. There is a significant spatiotemporal variation in PM2.5 concentration. In terms of spatial variation, high MAPC in the southeast region of the study area from January-April was observed, whereas PM2.5 concentration in the northwest region towards the central area was relatively higher from May-December. With regards to temporal variation, PM2.5 concentrations in January, February, November, and December were significantly higher than those in other months. The reason behind this may be related to indoor heating and climatic conditions. Some studies have reported that the burning of biomass and fossil energy for heating in winter generated huge PM2.5 emissions [47,48]. Dust storms, which frequently occur in North China in late winter and spring, is another major contribution that aggravates PM2.5 concentrations [49].   Figure 5. There is a significant spatiotemporal variation in PM 2.5 concentration. In terms of spatial variation, high MAPC in the southeast region of the study area from January-April was observed, whereas PM 2.5 concentration in the northwest region towards the central area was relatively higher from May-December. With regards to temporal variation, PM 2.5 concentrations in January, February, November, and December were significantly higher than those in other months. The reason behind this may be related to indoor heating and climatic conditions. Some studies have reported that the burning of biomass and fossil energy for heating in winter generated huge PM 2.5 emissions [47,48]. Dust storms, which frequently occur in North China in late winter and spring, is another major contribution that aggravates PM 2.5 concentrations [49].

Spatiotemporal Analysis of City-Scale PM 2.5 Concentration
Spatiotemporal variations of city-scale AAPC in the study area are shown in Figure 7. The AAPC in Taiyuan, Yangquan, Changzhi, and Jincheng changed slightly, but AAPC in the other cities decreased year by year, especially in Hengshui, Liaocheng, Dezhou, and Jinan. This might be related to the relatively high PM 2.5 pollution present in these cities. Nevertheless, AAPC in all cities were still higher than the national environmental air quality standard of 35 µg/m 3 (GB3095-2012), and much higher than the health standard recommended by the WHO of 10 µg/m 3 . This demonstrated that more work is required to control PM 2.5 concentrations in the study area.

Spatiotemporal Analysis of City-Scale PM2.5 Concentration
Spatiotemporal variations of city-scale AAPC in the study area are shown in Figure 7. The AAPC in Taiyuan, Yangquan, Changzhi, and Jincheng changed slightly, but AAPC in the other cities decreased year by year, especially in Hengshui, Liaocheng, Dezhou, and Jinan. This might be related to the relatively high PM2.5 pollution present in these cities. Nevertheless, AAPC in all cities were still higher than the national environmental air quality standard of 35 g/m 3 (GB3095-2012), and much higher than the health standard recommended by the WHO of 10 g/m 3 . This demonstrated that more work is required to control PM2.5 concentrations in the study area.   To further identify the local aggregation pattern of city-scale AAPC, the local Moran's I analysis was performed (Figure 9). The AAPC was dominated by high-high (HH) and low-low (LL) aggregation types; however, high-low (HL) or low-high (LH) types were not found. This indicated that AAPC in the study area had a local spatial positive correlation. HH-type regions, also known as the high-value aggregation region of AAPC, stably locate in the study area center, whereas the

Spatiotemporal Analysis of City-Scale PM2.5 Concentration
Spatiotemporal variations of city-scale AAPC in the study area are shown in Figure 7. The AAPC in Taiyuan, Yangquan, Changzhi, and Jincheng changed slightly, but AAPC in the other cities decreased year by year, especially in Hengshui, Liaocheng, Dezhou, and Jinan. This might be related to the relatively high PM2.5 pollution present in these cities. Nevertheless, AAPC in all cities were still higher than the national environmental air quality standard of 35 g/m 3 (GB3095-2012), and much higher than the health standard recommended by the WHO of 10 g/m 3 . This demonstrated that more work is required to control PM2.5 concentrations in the study area.   To further identify the local aggregation pattern of city-scale AAPC, the local Moran's I analysis was performed (Figure 9). The AAPC was dominated by high-high (HH) and low-low (LL) aggregation types; however, high-low (HL) or low-high (LH) types were not found. This indicated that AAPC in the study area had a local spatial positive correlation. HH-type regions, also known as the high-value aggregation region of AAPC, stably locate in the study area center, whereas the To further identify the local aggregation pattern of city-scale AAPC, the local Moran's I analysis was performed (Figure 9). The AAPC was dominated by high-high (HH) and low-low (LL) aggregation types; however, high-low (HL) or low-high (LH) types were not found. This indicated that AAPC in the study area had a local spatial positive correlation. HH-type regions, also known as the high-value aggregation region of AAPC, stably locate in the study area center, whereas the distribution of LL-type regions is unstable. The LL-type regions were concentrated in the west of the study area during 2015-2016, but in the northeast region of the study area in 2017.
distribution of LL-type regions is unstable. The LL-type regions were concentrated in the west of the study area during 2015-2016, but in the northeast region of the study area in 2017.

Effects of Socioeconomic Factors on PM2.5 Concentration
A series of tests were required to select the optimal analysis model. Firstly, to determine whether the spatial panel model should be applied, the Lagrange Multiplier (LM), and robust LM (RLM) tests were applied to assess the spatial correlation of the errors of the classic panel model. The LM error, LM lag, and RLM error under each condition (unfixed effects, spatial fixed effects, time-period fixed effects, and spatial and time-period fixed effects) all passed the 5% significance test (Table 3), indicating a significant spatial correlation of the errors in the classic panel model. Hence, the spatial panel model should be applied. The likelihood ratio (LR) tests of spatial fixed effects and time-period fixed effects both exceeded the 1% significance level, which proved the superiority of spatial and time-period fixed effects to the spatial fixed effects or time-period fixed effects. Subsequently, a spatial panel Durbin model with spatial and time-period fixed effects was constructed to test whether it could be simplified into SPLM or SPEM ( Table 4). The Hausman test was significant at the 1% level, indicating that the model with random effects was rejected; the Wald and LR tests were significant at the 1% level. Therefore, the spatial panel Dubin model could not be simplified.  The spatial panel Dubin model with spatial and time-period fixed effects was chosen to analyze how economic and social factors influence PM2.5 concentration. As shown in Table 5, UR has a significantly negative effect on PM2.5 concentration, and increasing the UR of a city by 1% could decrease PM2.5 concentration in the city by 0.78%. This is because, with the increasing demands for

Effects of Socioeconomic Factors on PM 2.5 Concentration
A series of tests were required to select the optimal analysis model. Firstly, to determine whether the spatial panel model should be applied, the Lagrange Multiplier (LM), and robust LM (RLM) tests were applied to assess the spatial correlation of the errors of the classic panel model. The LM error, LM lag, and RLM error under each condition (unfixed effects, spatial fixed effects, time-period fixed effects, and spatial and time-period fixed effects) all passed the 5% significance test (Table 3), indicating a significant spatial correlation of the errors in the classic panel model. Hence, the spatial panel model should be applied. The likelihood ratio (LR) tests of spatial fixed effects and time-period fixed effects both exceeded the 1% significance level, which proved the superiority of spatial and time-period fixed effects to the spatial fixed effects or time-period fixed effects. Subsequently, a spatial panel Durbin model with spatial and time-period fixed effects was constructed to test whether it could be simplified into SPLM or SPEM ( Table 4). The Hausman test was significant at the 1% level, indicating that the model with random effects was rejected; the Wald and LR tests were significant at the 1% level. Therefore, the spatial panel Dubin model could not be simplified.  The spatial panel Dubin model with spatial and time-period fixed effects was chosen to analyze how economic and social factors influence PM2.5 concentration. As shown in Table 5, UR has a significantly negative effect on PM 2.5 concentration, and increasing the UR of a city by 1% could decrease PM 2.5 concentration in the city by 0.78%. This is because, with the increasing demands for people to live in a heathy environment, China has implemented stricter environmental regulations than before. Considering the spatial interaction of UR, the PM 2.5 concentration in a city could decrease by 2.23% if the UR in the surrounding cities is increased by 1%. Our findings are different from studies before 2010 [50], which suggested that UR has a positive effect on PM 2.5 concentration. The reason behind this is that after China made carbon emission reduction commitments at the 2009 Copenhagen Climate Conference, the government, enterprises, and society have been vigorously promoting ecological sustainability, which effectively curbed the aggravation of air pollution, and reduced PM 2.5 concentration. Note: ***, ** and * represent significance at the 1%, 5% and 10% levels, respectively.
The coefficients of lnPGRP and (lnPGRP) 2 were significantly positive and negative, respectively, indicating the existence of an inverted U-shaped environmental Kuznets curve (EKC) of PM 2.5 concentration in the study area. With the increase of the per capita income level, the study area experienced a process of pollution first and then treatment, resulting in the PM 2.5 concentration first increasing and then decreasing, which is consistent with some existing research [51][52][53]. Although most cities in the study area except for Beijing, Tianjin, and Qinhuangdao are still in the stage of industrialization [54] and have a secondary industry-dominated structure, some cities shut down many small-sized, high-pollution, and high-energy-consumption enterprises in order to meet environmental quality requirements, and reached the peak of pollution ahead of schedule at the cost of low and medium economic growth. This is further demonstrated by the insignificant coefficient of PSIGDP.
Nevertheless, NIEDS has a significant positive effect on PM 2.5 concentrations. PM 2.5 concentrations in cities may increase by 0.15% and adjacent areas by 0.67% when NIEDS is increased by 1%. This is related to the fact that many NIEDS are resource-and energy-consuming enterprises. For example, the Hebei Iron and Steel Group is an ultra-large iron and steel group that ranks as the first in China and the second in the world in terms of crude steel output. Such enterprises provide important support for local employment and economic development. Given the path-dependence of their development, it is difficult to carry out energy saving and emission reduction measures immediately and thoroughly [55]. These enterprises will generate a large amount of dust pollution during production activities. Accordingly, the PM 2.5 concentration of a city is increased 0.043% when the ISDE of adjacent cities is increased by 1%.
In addition to the above socioeconomic factors, the PM 2.5 concentration of a city was also influenced by those in surrounding cities. The coefficient of the spatial lag term of PM 2.5 concentration was 0.6337 and passed the 1% significance test, which was mainly due to the transmission and diffusion of PM 2.5 . Wang et al. [56] also reported that the contribution rate of foreign sources to PM 2.5 concentrations in the Beijing-Tianjin-Hebei region was 23.4%.
To further identify the influence of different social factors, we calculated the direct, indirect, and total effects of socioeconomic factors on PM 2.5 concentration (Table 6). Among seven factors, the total effect of lnPGRP, (lnPGRP) 2 , lnUR, lnISDE, and lnNIEDS passed the significance test, indicating that these five factors influenced the spatiotemporal distribution of PM 2.5 concentration in the study area. The order of the degree of influence for these five factors was: lnPGRP > lnUR > lnNIEDS > (lnPGRP) 2 > lnISDE. The other factors may not significantly influence PM 2.5 concentration. To be specific, lnPGRP and lnUR are primary factors that influenced PM 2.5 concentration, and the lnPGRP, lnUR, and lnNIEDS have spillover effects on the PM 2.5 concentration in surrounding cities. Note: ***, ** and * represent significance at the 1%, 5%, and 10% levels, respectively.

Conclusions
This study investigated the effects of socioeconomic factors on the spatiotemporal distribution of PM 2.5 in Beijing-Tianjin-Hebei and surrounding areas during 2015-2017. First, an estimation model considering spatiotemporal heterogeneous relationships was developed to depict the spatiotemporal pattern of PM 2.5 concentration in the study area. Then, on the basis of analyzing the spatiotemporal evolution of PM 2.5 concentration, a spatial panel Dubin model was applied to determine how socioeconomic factors affect PM 2.5 concentration. Major conclusions of this research include: