Investigating the Impacts of Urbanization on PM2.5 Pollution in the Yangtze River Delta of China: A Spatial Panel Data Approach

Urbanization is a key determinant of fine particulate matter (PM2.5) pollution variability. However, there is a limited understanding of different urbanization factors’ roles in PM2.5 pollution. Using satellite-derived PM2.5 data from 2002 to 2017, we investigated the spatiotemporal evolution and the spatial autocorrelation of PM2.5 pollution in the Yangtze River Delta (YRD) region. Afterwards, the impacts of three urbanization factors (population urbanization, land urbanization and economic urbanization) on PM2.5 pollution were estimated by a spatial Durbin panel data model (SDM). Obtained results showed that: (i) PM2.5 pollution was larger in the north than in the south of YRD; (ii) Lianyungang and Yancheng cities had significant increasing trends in PM2.5 pollution from 2002 to 2017; (iii) the regional median center of PM2.5 pollution was observed in the Nanjing city, with gradual shifting to the northwest during the 16-year period; (iv) PM2.5 pollution showed significant and positive spatial autocorrelation and spillover effect; (v) population urbanization contributed more to the increase in PM2.5 pollution than land urbanization, while economic urbanization had no significant impact. The present study highlights the impacts of three urbanization factors on PM2.5 pollution which represent valuable and relevant information for air pollution control and urban planning.


Introduction
Fine particulate matter (PM 2.5 ) mostly originates from anthropogenic activities such as construction, traffic, industrial production, cooking, waste incineration, biomass burning, etc. [1,2]. With small diameter and strong toxicological properties, PM 2.5 can easily enter into the human body and induce cardiovascular and respiratory diseases [3][4][5]. Moreover, high levels of PM 2.5 reduce atmospheric visibility and disturb the earth's radiative balance [6,7]. Recently, numerous urban areas in China have suffered frequent heavy PM 2.5 pollution [8]. For instance, two large-scale and long-lasting haze  Table S1). The administrative boundary layer with a scale of 1:250,000 was obtained from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences [51].

Data
Firstly, the overall spatial pattern and dynamics of PM2.5 pollution were analyzed at the pixellevel. Then, we examined the spatial dependence of PM2.5 pollution by spatial autocorrelation analysis and explored the impacts of urbanization on PM2.5 pollution by spatial regression analysis, both of which were conducted at the city-level. Accordingly, the pixel-level and city-level annual PM2.5 data were prepared for further analyses (more details in Section 2.2.1). Regarding the impacts of urbanization on PM2.5 pollution, we mainly focused on three common forms of urbanization: population urbanization, land urbanization and economic urbanization. These urbanization forms can be expressed by the ratio of urban population, the ratio of urban built-up area and per capita gross regional product, respectively [11]. For the purpose of reducing omitted variable bias, we also considered secondary industry, vegetation coverage, precipitation and wind speed in spatial regression analysis [41,52]. Table 1 lists the specific definitions of all variables, while data sources for the PM2.5 concentration and the seven explanatory variables are presented in Sections 2.2.1 and 2.2.2, respectively.  Table S1). The administrative boundary layer with a scale of 1:250,000 was obtained from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences [51].

Data
Firstly, the overall spatial pattern and dynamics of PM 2.5 pollution were analyzed at the pixel-level. Then, we examined the spatial dependence of PM 2.5 pollution by spatial autocorrelation analysis and explored the impacts of urbanization on PM 2.5 pollution by spatial regression analysis, both of which were conducted at the city-level. Accordingly, the pixel-level and city-level annual PM 2.5 data were prepared for further analyses (more details in Section 2.2.1). Regarding the impacts of urbanization on PM 2.5 pollution, we mainly focused on three common forms of urbanization: population urbanization, land urbanization and economic urbanization. These urbanization forms can be expressed by the ratio of urban population, the ratio of urban built-up area and per capita gross regional product, respectively [11]. For the purpose of reducing omitted variable bias, we also considered secondary industry, vegetation coverage, precipitation and wind speed in spatial regression analysis [41,52]. Table 1 lists the specific definitions of all variables, while data sources for the PM 2.5 concentration and the seven explanatory variables are presented in Sections 2.2.1 and 2.2.2, respectively. The annual PM 2.5 dataset at a grid resolution of about 1 km was obtained from the Dalhousie University Atmospheric Composition Analysis Group for the period 2002-2017 [53]. The PM 2.5 data were produced by integrating multiple types of information, including satellite retrievals, simulations, and ground measurements [54][55][56]. The dataset has relatively high accuracy (R 2 = 0.81) [54] and was widely used at different spatial scales [12,16,57]. For the analyses of the spatiotemporal evolution of PM 2.5 pollution at the pixel-level, the original gridded PM 2.5 data (a resolution of 1 km) was employed. For the analyses of spatial autocorrelation and spatial regression at the city-level, the mean PM 2.5 values for each city were calculated based on administrative boundaries.

Urbanization and Other Factors
Data for seven explanatory variables for the period 2002-2017 were gathered from multiple sources. Urbanization and industrial information that includes the ratio of urban population, the ratio of urban built-up area, per capita gross regional product and secondary industry, were collected from the China City Statistical Yearbook (2003-2018) [58]. Notably, the per capita gross regional product data were converted to constant 2002 prices (Chinese Yuan) using GDP deflator [59]. The annual normalized differential vegetation index (NDVI) data with a 1 km spatial resolution, reflecting vegetation growth and coverage, were provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences [51]. The annual precipitation and average wind speed from meteorological stations in and around the YRD were derived from the China Meteorological Data Service Center [60]. The original meteorological measurements were interpolated to 1 km grid using the inverse distance weighting method in order to generate continuous surfaces of precipitation and wind speed [40]. In addition, NDVI, precipitation and wind speed were averaged within administrative boundaries for each city.

Variables Selection
For the spatial regression analysis in Section 2.3.4, all variables, including PM, PU, LU, EU, SI, NDVI, Prec and Wind (See Table 1 for their full names and descriptions) were appropriated natural logarithmic forms to reduce the potential heteroscedasticity (Table S2). In addition, variance inflation factors (VIF) were calculated for the seven explanatory variables. Obtained VIF values are less than 10 (Table S3), thus suggesting that the explanatory variables exhibit no multi-collinearity [61,62].

Trend Analysis
The nonparametric Mann-Kendall test is powerful in identifying the significant trend in time series data [63][64][65][66] and the nonparametric Sen's slope estimator [67,68] provides a robust estimate of the slope. Compared with the parametric approaches, nonparametric approaches are distribution free, less affected by outliers, and hence more applicable for time series analysis [69,70]. In this study, the Mann-Kendall and Sen's slope tests were used to identify trends in PM 2.5 time series from 2002 to 2017. Firstly, for each pixel, we performed the Mann-Kendall test to detect whether the PM 2.5 time series has a significant trend. If the p-value is less than 0.05, the corresponding trend is considered significant [71]. Then, the slope values of the PM 2.5 data were calculated by performing Sen's slope test in order to indicate the direction and the magnitude of the temporal changes in PM 2.5 concentration. The two tests were performed by using the trend package in R software (version 3.2.2, R Foundation for Statistical Computing, Vienna, Austria).

Standard Deviation Ellipse Analysis
We employed the standard deviation ellipse (SDE) method to depict the spatial pattern changes of PM 2.5 pollution during 2002-2017. Three basic parameters of the SDE, including median center, long-short axis ratio and azimuth, were used to measure the dispersion and directional tendencies of elements in our study [72]. The median center is the point that has a minimum sum of distances to all the elements, thus identifying the location where the elements are most concentrated [73]. Accordingly, the change in the median center can reflect the overall moving trace of elements. The long and short half axes illustrate the primary and secondary distribution direction of the elements. The long-short axis ratio (i.e., the ratio of the lengths of the long axis to the short axis), an indicator describing the shape of an ellipse, can determine whether the elements are densely distributed in a certain direction or evenly distributed in the ellipse. The smaller the value, the less probability the elements are concentrated in a certain direction [74]. The azimuth, which ranges from 0 to 180 degrees, refers to the angle deviating clockwise from the north direction to the long axis direction, revealing the leading direction of elements [75]. The SDE analysis was performed by using the median center and directional distribution tools in ArcGIS software (version 10.2, ESRI Inc., Redlands, CA, USA).

Spatial Autocorrelation Analysis
Spatial autocorrelation means that the nearby observations of a variable tend to be positively or negatively correlated [76]. To examine the necessity of constructing spatial regression models, we assessed the spatial autocorrelation of PM 2.5 pollution in the YRD using global Moran's I [77]. The Moran's I ranges from −1 to 1: a) when 0 < Moran's I < 1, it denotes a positive spatial autocorrelation and a clustering pattern of PM 2.5 pollution, b) when −1 < Moran's I < 0, it denotes a negative spatial autocorrelation and a dispersion pattern of PM 2.5 pollution and c) when Moran's I = 0, it denotes a random spatial distribution. The p-value reflects the significance level and additional details about this method can be found in the literature [78,79]. In this study, the global Moran's I of PM 2.5 pollution for each year between 2002 and 2017 was estimated based on the inverse distance spatial weight matrix [80,81] by employing the spatial autocorrelation tool in ArcGIS software (version 10.2, ESRI Inc., Redlands, CA, USA).

Spatial Regression Analysis
To cope with the potential spatial dependence in PM 2.5 pollution among the cities, we considered three types of spatial regression models when analyzing the impacts of urbanization on PM 2.5 pollution [82]. The first model is the spatial lag model (SLM), which includes a spatially lagged dependent variable and hypothesizes that the PM 2.5 pollution in the city i is influenced by the PM 2.5 pollution in the neighboring cities [83]. In our case, the SLM model can be established as follows: where β 1 , β 2, . . . , β 7 represent the seven regression coefficients of the explanatory variables; ρ denotes the spatial autoregressive coefficient measuring the intensity of spatial dependence in PM 2.5 pollution among cities; W stands for the spatial weight matrix, which is the inverse distance spatial weight matrix in this case; µ i and ξ t correspond to spatial effect and time effect, respectively. ε it is a normally distributed random error terms. The second model is the spatial error model (SEM), which contains a lagged error term and hypothesizes that the errors in neighboring cities are likely to interact. It can be constructed as follows: where φ it represents the spatial autocorrelation error term; λ is the spatial autocorrelation coefficient of the error term. The third model is the spatial Durbin model (SDM), which incorporates both the spatially lagged dependent variable and explanatory variables, considering not only the spatial dependence of the dependent variable but also that of the explanatory variables [84,85]. It can be defined as follows: where θ 1 , θ 2 , . . . , θ 7 refer to the spatial lag coefficients of the explanatory variables. Several model tests help to select an appropriate spatial regression model. Following the standard approach proposed by Elhorst [86], the spatial regression analysis starts with a non-spatial model estimated by OLS. Afterwards, the Likelihood ratio (LR) joint significance test can be performed in order to determine whether the model should control for spatial fixed or time fixed effects [86]. In this study, the random effect was not considered because the analysis focuses on the population rather than the sample [87]. Then, the Lagrange multiplier (LM) tests and Robust LM tests can be used to investigate whether there exist spatial interaction effects in dependent variables or errors. If the result justifies the existence of the spatial interaction effect, the SDM model is recommended to be established [85]. In this case, to access whether the SDM can be degraded into SLM or SEM, the Wald and LR tests should be carried out. By contrast, if there are no spatial interaction effects, the spatial regression model should not be selected.
The above spatial models can be estimated by the maximum likelihood (ML) method. Notably, since the SLM and SDM involve spatial lags of the dependent variable, the coefficients of explanatory variables (i.e., β 1 , β 2, . . . , β 7 ) in the two models could not represent the marginal effect [83]. Therefore, direct and indirect effects coefficients should be further estimated. In this study, the direct effect refers to the impact on PM 2.5 pollution of the city i resulted from a change in an explanatory variable of the city i. The indirect (spillover) effect illustrates the cumulative impacts on PM 2.5 pollution caused by changes in an explanatory variable of other cities. The sum of the two effects is the total effect. More details about the estimation of the direct and indirect effects can be found in LeSage and Pace [85]. The spatial regression analysis was performed in the MATLAB software (version R2016a, MathWorks Inc., Natick, MA, USA).

Spatial Pattern of PM 2.5 Pollution
The distinct spatial variation in the 16-

Trends of PM2.5 Concentrations
The significance levels (i.e., p-values) and slopes of the trends in PM2.5 concentrations were calculated by methods mentioned in Section 2.3.1. Figure 3a shows the p-values of the upward (slope > 0) and downward (slope < 0) trends, while Figure 3b shows the slopes of the trends. The results showed that approximately 16.12% and 0.98% of the study area had significant (p < 0.05) upward and downward trends, respectively ( Figure 3a). The north and central YRD were characterized by upward trends with larger parts of Lianyungang (LYG) and Yancheng (YC) cities having the greatest concentration growth (slope > 1.2) at a high significance level (p < 0.01) (Figure 3b). Xuzhou (XZ), Huai'an (HA), Yangzhou (YZ) and Taizhou (TZ) cities also witnessed the significant increase in PM2.5 concentration with slope values exceeding 0.9 in larger parts of these cities. By contrast, downward trends were noticed in the south YRD with most slopes being not significant (p ≥ 0.05). Only a few parts of the Lishui (LS) and Quzhou (QZ) cities illustrated significant (p < 0.05) decrease (slope < 0) in PM2.5 concentration.

Trends of PM 2.5 Concentrations
The significance levels (i.e., p-values) and slopes of the trends in PM 2.5 concentrations were calculated by methods mentioned in Section 2.3.1. Figure 3a shows the p-values of the upward (slope > 0) and downward (slope < 0) trends, while Figure 3b shows the slopes of the trends. The results showed that approximately 16.12% and 0.98% of the study area had significant (p < 0.05) upward and downward trends, respectively (Figure 3a). The north and central YRD were characterized by upward trends with larger parts of Lianyungang (LYG) and Yancheng (YC) cities having the greatest concentration growth (slope > 1.2) at a high significance level (p < 0.01) (Figure 3b). Xuzhou (XZ), Huai'an (HA), Yangzhou (YZ) and Taizhou (TZ) cities also witnessed the significant increase in PM 2.5 concentration with slope values exceeding 0.9 in larger parts of these cities. By contrast, downward trends were noticed in the south YRD with most slopes being not significant (p ≥ 0.05). Only a few parts of the Lishui (LS) and Quzhou (QZ) cities illustrated significant (p < 0.05) decrease (slope < 0) in PM 2.5 concentration.  Table 2). The moving directions and the rates of the median center in the three five-year intervals could be computed according to its coordinates ( Table 2).  Table 2), reflecting that the directionality of the PM 2.5 pollution was weakened and that the dispersion of pollution decreased in the long axis direction (i.e., "southeast-northwest") and expanded in the short axis direction (i.e., "northeast-southwest"). Furthermore, the azimuth changed only slightly downward trends, respectively (Figure 3a). The north and central YRD were characterized by upward trends with larger parts of Lianyungang (LYG) and Yancheng (YC) cities having the greatest concentration growth (slope > 1.2) at a high significance level (p < 0.01) (Figure 3b). Xuzhou (XZ), Huai'an (HA), Yangzhou (YZ) and Taizhou (TZ) cities also witnessed the significant increase in PM2.5 concentration with slope values exceeding 0.9 in larger parts of these cities. By contrast, downward trends were noticed in the south YRD with most slopes being not significant (p ≥ 0.05). Only a few parts of the Lishui (LS) and Quzhou (QZ) cities illustrated significant (p < 0.05) decrease (slope < 0) in PM2.5 concentration.   Table 2 displays the basic parameters of the ellipses. The median center of the pollution was found in northwestern Nanjing city (Figure 4), moving from 118.63° E and 31.87° N in 2002 to 118.50° E and 32.01° N in 2017 ( Table 2). The moving directions and the rates of the median center in the three five-year intervals could be computed according to its coordinates ( Table 2).  Table 2), reflecting that the directionality of the PM2.5 pollution was weakened and that the dispersion of pollution decreased in the long axis direction (i.e., "southeastnorthwest") and expanded in the short axis direction (i.e., "northeast-southwest"). Furthermore, the azimuth changed only slightly from 143.58° to 143.97° (Table 2), which indicated that the leading direction of PM2.5 concentrations had no distinct change.          Table 4 reports the results of the spatial diagnostic tests for non-spatial regression models. The results of the LR joint significance tests provided evidence about the existence of the spatial fixed effects (1370.892, p < 0.01) and time fixed effects (722.186, p < 0.01). This indicated that the panel data model including both the spatial and time effects was a suitable choice. Regarding the type of spatial regression model, as demonstrated in the first four rows (Table 4), the LM tests and Robust LM tests were statistically significant (p < 0.05) for all the model specifications. This result implies that the spatial interaction effects in both dependent variables and error terms could not be neglected. Thus, by following the general-to-specific approach, we established the SDM model with spatial and time fixed effects and subsequently implemented the Wald and LR tests. As highlighted in Table 5, all statistics of Wald and LR tests were significant at 1% level, which verified the superiority of SDM over SLM and SEM. Accordingly, the SDM model under the spatial and time fixed effects was chosen to describe the data, and the subsequent analysis was based on it.

Estimation Results
The SDM model was fitted with the maximum likelihood (ML) estimator by considering the spatial and time fixed effects ( Table 6). The spatial autoregressive coefficient ρ is highly significant (p < 0.01) with a value of 0.867, thus indicating a strong link between the PM 2.5 pollution in neighboring cities. Due to the feedback effect, the spatial regression coefficients cannot truly illustrate the marginal effects of urbanization, secondary industry, NDVI, precipitation and wind speed. Hence, we further computed the direct, indirect and total effects of each variable on PM 2.5 pollution (Table 7). Four variables (i.e., lnPU, lnLU, lnSI and lnNDVI) showed significant total effects with the coefficients of 4.084, 2.182, 1.379 and −5.988, respectively. For the direct effect, lnPU had the largest significant positive coefficient (0.310), followed by lnLU (0.149) and lnSI (0.116), suggesting that every 1% increase in lnPU, lnLU and lnSI will directly increase lnPM by 0.310%, 0.149%, and 0.116% in the local city, respectively. In contrast, lnNDVI, lnPrec and lnWind had negative direct effects with the respective coefficients of −0.541, −0.031 and −0.283, thus implying that NDVI, precipitation and wind speed exerted direct negative impacts on PM 2.5 pollution. By considering the indirect effect, lnPU, lnLU and lnSI showed significant and positive indirect effects. Specifically, every 1% growth in these three variables in nearby cities can increase lnPM by 3.774%, 2.033% and 1.263% in the local city, respectively. Moreover, lnNDVI had a considerable high negative and significant indirect effect (−5.447). For the lnEU, which represented the economic urbanization, neither the direct effect nor the indirect effect was significant.

The Impacts of Urbanization on PM 2.5 Pollution
The present study highlighted the different impacts of the three forms of urbanization on PM 2.5 pollution for the YRD during the study period (Table 7). Complying with the previous studies [41,45,88], population urbanization played a dominant role in increasing PM 2.5 pollution and land urbanization was also a large contributor to PM 2.5 pollution. Regarding economic urbanization, some studies argued that it can lead to an increase in PM 2.5 concentration [12,47]. However, it showed a negative and insignificant impact on PM 2.5 levels in the YRD. Based on our results, the underlying mechanisms that different forms of urbanization impact PM 2.5 pollution are illustrated in Figure 5.
Atmosphere 2018, 9, x FOR PEER REVIEW 12 of 19 [41,45,88], population urbanization played a dominant role in increasing PM2.5 pollution and land urbanization was also a large contributor to PM2.5 pollution. Regarding economic urbanization, some studies argued that it can lead to an increase in PM2.5 concentration [12,47]. However, it showed a negative and insignificant impact on PM2.5 levels in the YRD. Based on our results, the underlying mechanisms that different forms of urbanization impact PM2.5 pollution are illustrated in Figure 5. Population urbanization describes the migration of the population from rural areas to urban areas, mainly impacting PM2.5 pollution in two ways. On one hand, the urban population agglomeration is usually accompanied by high demands for new housing, public infrastructures and private cars, which will inevitably raise energy consumption [12]. The growing energy consumption contributes to increased exhaust emissions, thereby enhancing air pollution. On the other hand, as the population increases and the people's living standards improve, large amounts of domestic waste are produced in the urban area. The released aerosol particles due to the improper garbage disposal may further increase PM2.5 concentration [10].
Land urbanization refers to the process of urban sprawl. The impact of land urbanization on PM2.5 pollution may be explained in three aspects. Firstly, the larger the built-up area, the longer the commuting distance for residents, which will increase the traffic flows [89]. As a consequence, emissions from vehicles will contribute substantially to PM2.5 pollution. Secondly, the majority of the natural land covers (e.g., forest land and grassland) are replaced by artificial surfaces during the land urbanization process. Compared with artificial surfaces, forest and grass can purify the air through obstructing and absorbing the particulate matter. Accordingly, land cover changes may weaken the self-purification ability of urban ecosystems [18]. Lastly, the high-rise and high-density buildings can hinder the dilution and dispersion of air pollutants [90].
Economic urbanization can be regarded as the process of economic growth. It influences PM2.5 pollution from opposite directions. For one thing, rapid economic growth can induce more production activities that leads to increased pollution emissions [45]. For another thing, economic urbanization is usually accompanied by the increased energy efficiency, improved industrial technologies and raised public environmental awareness, which can alleviate air pollution. Population urbanization describes the migration of the population from rural areas to urban areas, mainly impacting PM 2.5 pollution in two ways. On one hand, the urban population agglomeration is usually accompanied by high demands for new housing, public infrastructures and private cars, which will inevitably raise energy consumption [12]. The growing energy consumption contributes to increased exhaust emissions, thereby enhancing air pollution. On the other hand, as the population increases and the people's living standards improve, large amounts of domestic waste are produced in the urban area. The released aerosol particles due to the improper garbage disposal may further increase PM 2.5 concentration [10].
Land urbanization refers to the process of urban sprawl. The impact of land urbanization on PM 2.5 pollution may be explained in three aspects. Firstly, the larger the built-up area, the longer the commuting distance for residents, which will increase the traffic flows [89]. As a consequence, emissions from vehicles will contribute substantially to PM 2.5 pollution. Secondly, the majority of the natural land covers (e.g., forest land and grassland) are replaced by artificial surfaces during the land urbanization process. Compared with artificial surfaces, forest and grass can purify the air through obstructing and absorbing the particulate matter. Accordingly, land cover changes may weaken the self-purification ability of urban ecosystems [18]. Lastly, the high-rise and high-density buildings can hinder the dilution and dispersion of air pollutants [90].
Economic urbanization can be regarded as the process of economic growth. It influences PM 2.5 pollution from opposite directions. For one thing, rapid economic growth can induce more production activities that leads to increased pollution emissions [45]. For another thing, economic urbanization is usually accompanied by the increased energy efficiency, improved industrial technologies and raised public environmental awareness, which can alleviate air pollution. Consequently, the final statistically insignificant direct and indirect effect coefficients of economic urbanization in the YRD may be attributed to the dual effects of economic urbanization. Liu et al. [41] also found that per capita GDP had no significant influence on air quality.

The Impacts of other Factors on PM 2.5 Concentrations
In addition to urbanization factors, we have introduced the secondary industry, NDVI, precipitation and wind speed into our model to reduce the omitted variable bias. The results showed that the secondary industry enhanced PM 2.5 levels (Table 7), which is in accordance with a previous study [47]. That is because the secondary industry, especially heavy industry and real estate, consume large amounts of energy and generates high concentrations of atmospheric pollutants. Vegetation coverage (Table 7), represented by NDVI, is able to partially absorb and block particulate matter [6], which contributed to the decline in PM 2.5 concentrations. Additionally, precipitation and wind speed played important roles in mitigating the local PM 2.5 pollution (Table 7) as precipitation can remove the particles by wet deposition while wind facilitates the dispersion of particles [91].

Policy Implications
The findings from this empirical study contributed to the development of policy implications with the aim to alleviate PM 2.5 pollution.
Firstly, the spillover effect of PM 2.5 pollution implies that inter-regional cooperation is necessary for preventing, mitigating and controlling air pollution. Accordingly, an information sharing system should be established at the regional level with the aim to provide the information related to air quality monitoring, pollution emissions and emergency response to air pollution accidents [92]. Meanwhile, unified environmental laws and regulations should be formulated and implemented to restrict the transfer of pollution among cities [12].
Secondly, governments should emphasize on the popularization of non-fossil fuel energy and citizen's environmental awareness. As demonstrated by the results, the growth in urban population and built-up area contributed substantially to the energy consumption and air pollutant emissions, thereby increasing PM 2.5 concentration. Hence, it is crucial to reduce coal consumption and encourage the usage of clean energy [88]. In addition, enhancing green consumption, promoting the use of public transportation and advocating garbage classification can increase the public's awareness of air protection [11].
Finally, attention should be devoted to industrial structure and urban planning. It is urgent to close the highly polluting industry, transform the traditional industries with advanced technologies and to encourage the development of green industry [12]. In addition, the impacts of vegetation coverage, precipitation and wind speed on PM 2.5 concentration should not be overlooked. To reduce the risk of air pollution exposure, urban planning measures should be adopted that include designing urban ventilation corridors, allocating large scale green space in the downwind direction of pollution source and arranging suitable tree species in the street [93].

Limitation
This present study has a few limitations. Firstly, the used socioeconomic data employed were derived based on the administrative boundaries, which may, to a certain extent, impact the estimation results. To more accurately characterize urbanization, grid data (e.g., satellite-derived nighttime light data and population data) should be considered in the follow-up studies. Secondly, the analysis was conducted at the annual time scale with the emphasis on detecting the impacts of urbanization on PM 2.5 pollution. Accordingly, the influencing mechanisms of precipitation and wind speed on air pollution could not be revealed completely and it would be beneficial to perform the spatial regression analysis on shorter time scales (e.g., day, month and season).

Conclusions
The present study analyzed the spatiotemporal evolution of PM 2.5 pollution at the pixel-level in the YRD from 2002 to 2017. Then, the spatial autocorrelation of PM 2.5 pollution among the 41 cities in the study domain was examined. Finally, we investigated the impacts of different urbanization factors (population urbanization, land urbanization and economic urbanization) on PM 2.5 pollution at the city-level. The key findings from this study are as follows: • Distinct variation in the spatial pattern of the PM 2.5 concentrations was identified in the YRD. Approximately 83.57% of the YRD had mean PM 2.5 concentrations higher than 35 ug/m 3 (WHO Interim Target 1). The pollution was more severe in the north than in the south part of the YRD.

•
The significant and largest increasing trends in PM 2.5 pollution was noticed in Lianyungang and Yancheng cities. However, the regional median center of the PM 2.5 pollution was located in the Nanjing city and was gradually shifting to the northwest during the research period.

•
The positive spatial autocorrelation and spillover effect of PM 2.5 pollution were verified by the global Moran's I values and the SDM model, respectively, thus suggesting that the prevention and control of air pollution should rely on inter-regional cooperation.

•
Population urbanization, land urbanization, secondary industry and NDVI exerted significant impacts on PM 2.5 pollution. Population urbanization had the largest positive total effect, followed by land urbanization and secondary industry, while NDVI showed a negative total effect. No significant effect was detected for economic urbanization.
In general, we have noticed that the impact mechanisms of different forms of urbanization on PM 2.5 pollution varied. The investigation of these impact mechanisms can help policymakers and urban planners in developing measures to abate air pollution.

Conflicts of Interest:
The authors declare no conflict of interest.