Estimating Regional PM_2.5 Concentrations in China Using a Global-Local Regression Model Considering Global Spatial Autocorrelation and Local Spatial Heterogeneity

Abstract

Linear regression models are commonly used for estimating ground PM_2.5 concentrations, but the global spatial autocorrelation and local spatial heterogeneity of PM_2.5 distribution are either ignored or only partially considered in commonly used models for estimating PM_2.5 concentrations. Therefore, taking both global spatial autocorrelation and local spatial heterogeneity into consideration, a global-local regression (GLR) model is proposed for estimating ground PM_2.5 concentrations in the Yangtze River Delta (YRD) and in the Beijing, Tianjin, Hebei (BTH) regions of China based on the aerosol optical depth data, meteorological data, remote sensing data, and pollution source data. Considering the global spatial autocorrelation, the GLR model extracts global factors by the eigenvector spatial filtering (ESF) method, and combines the fraction of them that passes further filtering with the geographically weighted regression (GWR) method to address the local spatial heterogeneity. Comprehensive results show that the GLR model outperforms the ordinary GWR and ESF models, and the GLR model has the best performance at the monthly, seasonal, and annual levels. The average adjusted R² of the monthly GLR model in the YRD region (the BTH region) is 0.620 (0.853), which is 8.0% and 7.4% (6.8% and 7.0%) higher than that of the monthly ESF and GWR models, respectively. The average cross-validation root mean square error of the monthly GLR model is 7.024 μg/m³ in the YRD region, and 9.499 μg/m³ in the BTH region, which is lower than that of the ESF and GWR models. The GLR model can effectively address the spatial autocorrelation and spatial heterogeneity, and overcome the shortcoming of the ordinary GWR model that overfocuses on local features and the disadvantage of the poor local performance of the ordinary ESF model. Overall, the GLR model with good spatial and temporal applicability is a promising method for estimating PM_2.5 concentrations.

1. Introduction

PM_2.5, which denotes particulate matter with an aerodynamic diameter ≤2.5 μm, poses negative effects on human health [1,2,3]. The Global Burden of Disease (GBD) has reported that 2.9 million people died worldwide in 2013 due to PM_2.5 [4]. Numerous studies show that PM_2.5 has been implicated as among the key risk factors for cancer [5,6,7]. The combination of PM_2.5 and water vapor also produces the haze phenomenon, causing reduced visibility and inconvenience for people to live and work [8]. In recent years, with rapid urbanization and industrialization, PM_2.5 pollution in China has become more and more severe [9,10,11], so the monitoring of PM_2.5 is especially important. However, although China has been building PM_2.5 monitoring networks since 2013, these sparse and unevenly distributed PM_2.5 monitoring stations are insufficient to support large-scale, accurate PM_2.5 monitoring.

To fill the gap in ground measurements, large-scale remote sensing products were applied to estimate ground-level PM_2.5 concentrations [12,13,14]. Aerosol optical depth (AOD) has been widely recognized to be closely correlated to PM_2.5 concentration [15,16,17]. AOD data and other auxiliary factors, such as pollution source data, meteorological data, and elevation data, are commonly used to estimate PM_2.5 concentrations [18,19,20]. Continuous ground-level PM_2.5 concentrations can be obtained by many methods using continuous remote sensing data. Commonly used methods include machine learning methods, chemical transport models, and deep learning methods [21,22,23].

Linear regression models are commonly used in estimating PM_2.5 concentrations, including land use regression (LUR), multiple linear regression (MLR), and least squares regression (OLS) [24,25,26,27]. Researchers have applied these models to estimate PM_2.5 and explore the linear relationship between the variables and PM_2.5. However, the distribution of PM_2.5 is spatially autocorrelated and spatially heterogeneous, and the above models do not perform well in PM_2.5 concentration estimation which may be due to their lack of consideration of spatial effects [28].

To address the spatial autocorrelation problem, the spatial error model (SEM), the spatial lag model (SLM), and the spatial Durbin model (SDM) [29,30,31] are typically applied. In particular, the eigenvector spatial filtering (ESF) method, proposed by Griffith et al. [32], is widely used to filter out spatial autocorrelation by calculating and filtering the eigenvectors of the spatial weight matrix, which are then introduced as independent variables together with other variables to construct the model. The ESF regression model was developed to estimate PM_2.5 concentrations by introducing the spatial eigenvectors, which were selected from the spatial weight matrix constructed by the PM_2.5 monitoring stations [33]. However, although ESF effectively solves the problem of spatial autocorrelation, the ESF model with constant regression coefficients has poor local fit and ignores the spatial heterogeneity of PM_2.5 distribution.

Considering the spatial heterogeneity in PM_2.5 distribution, geographically weighted regression (GWR) [34,35,36], Bayesian maximum entropy (BME) [20,37], and Bayesian spatially varying coefficient (SVC) [38] have been applied to solve this issue and improve the accuracy of PM_2.5 concentration estimation. The GWR model, a local model proposed by Brunsdon et al. [39], is among the most commonly used local models, and its regression coefficients vary with spatial location, greatly improving the accuracy of local fits. Some researchers have continuously improved the GWR model and used it for PM_2.5 concentration estimation. Examples of these models include the multiscale GWR (MGWR) [40], geographically and temporally WR (GTWR) [41,42], geographically neural network WR (GNNWR) [43], and the principal component analysis-GWR (PCA-GWR) [44] models. All these models effectively solve the local spatial heterogeneity problem of PM_2.5 regression modeling and further expand the time series to optimize the model and improve the estimation accuracy, but these models ignore global spatial autocorrelation and overfocus on local features, and suffer from local overfitting. The above regression models either consider spatial autocorrelation or spatial heterogeneity and do not take both into account simultaneously.

In this paper, this study aims to address both global spatial autocorrelation and local spatial heterogeneity, and a global-local regression (GLR) model is proposed for estimating ground PM_2.5 concentrations in the Yangtze River Delta (YRD) and in the Beijing, Tianjin, Hebei (BTH) regions. The GLR model can effectively address both the global spatial autocorrelation and local spatial heterogeneity, and it also overcomes the shortcoming of the ordinary GWR model that overfocuses on local features and the disadvantage of the poor local performance of the ordinary ESF model. In order to verify the spatial applicability of the GLR model and exclude the possibility of accuracy inflation caused by spatial factors, we selected two regions as the study area. In addition, to verify the temporal applicability of the GLR model, we not only constructed monthly models but also conducted experiments on seasonal and annual time scales.

2. Materials and Methods

2.1. Study Areas

The YRD region (Figure 1a) is a flat alluvial plain covering an area of over 350,000 square kilometers on the east coast of China. The YRD region consists of Shanghai, Anhui province, Jiangsu province, and Zhejiang province. It is an important economic and cultural center in China, with a total population of over 150 million, and is an important international gateway to the Asia-Pacific region. However, along with its rapid economic growth, the YRD region has suffered from severe environmental pollution, especially air pollution, and people’s health is under serious threat [45]. There are 234 national monitoring stations distributed throughout the region, but most of them are located at low altitudes.

The BTH region (Figure 1b), consisting of the two municipalities of Beijing and Tianjin, and Hebei Province, is the most dynamic and largest economic region in northern China. With a population of more than 110 million and an area of over 210,000 square kilometers, it is among the most economically promising regions in China, and also the most polluted region in China [46]. Despite the high level of air pollution, there are only 79 national monitoring stations here.

2.2. Data Integration

In addition to the ground PM_2.5 data, the datasets used in this study include AOD, meteorological data including planetary boundary layer height (PBLH), surface pressure (PS), surface temperature (TS), and relative humidity (RH), the normalized difference vegetation index (NDVI), elevation (DEM), and the kernel densities of roads (ROAD) and factories (FACT). All data were measured from 1 December 2015 to 30 November 2016 over the YRD and BTH regions.

Ground PM_2.5 data were obtained from the China Environmental Monitoring Center (http://106.37.208. 233:20035 (accessed on 21 December 2021)). The dataset included 24 h concentrations of all air pollutants from all monitoring stations in China. The monthly, seasonal, and annual average PM_2.5 data were obtained from hourly PM_2.5 concentrations.

The details of remote sensing data and meteorological data are shown in

Table 1. The data sources of remote sensing and meteorological data.

Factors		Source	Spatial Resolution	Temporal Resolution
AOD		MCD19A2 https://ladsweb.modaps.eosdis.nasa.gov/search/ (accessed on 10 March 2022)	1 km (BTH) 3 km (YRD)	Daily
DEM		Shuttle Radar Topography Mission (SRTM) DEM http://srtm.csi.cgiar.org1 (accessed on 9 March 2022)	90 m	/
Meteorological Data	TS	European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 (https://cds.climate.copernicus.eu) (accessed on 9 March 2022)	0.25° latitude × 0.25° longitude	Hourly
	PBLH
	PS
	RH
NDVI		MOD13A3 https://ladsweb.modaps.eosdis.nasa.gov/search/ (accessed on 10 March 2022)	1 km	Monthly

. The meteorological data (PS, PBLH, RH, TS) has a spatial resolution of 0.25° latitude × 0.25° longitude, which is coarser than that of the AOD data (1 km in the BTH region and 3 km in the YRD region). For consistency with the AOD data, they were interpolated to 1 km (3 km in the YRD region) by ordinary kriging [47,48]. Daily AOD data and hourly meteorological data were averaged to obtain monthly, seasonal, and annual data. Likewise, the DEM data were also resampled to the same spatial resolution as the AOD data.

Pollution source data, including road networks and factory locations, were derived from BaiduMap (https://map.baidu.com (accessed on 14 March 2022)) and OpenStreetMap (https://www.openstreetmap.org/ (accessed on 14 March 2022)), respectively. Point-based factory location data and line-based road network data can be converted into kernel density raster data in ArcGIS, and the buffer range and spatial resolution were set to 24 km [33] and 1 km (3 km in the YRD region), respectively.

Subsequently, all values of independent variables at the monitoring stations were extracted from the corresponding raster. Further, all records with null data values were removed. We obtained a total of 1257 and 3651 records at monthly, seasonal, and annual time scales in the BTH region and the YRD region, respectively. In each period, there were approximately 74 records in the BTH region, and approximately 215 records in the YRD region.

2.3. Global-Local Regression Modeling Method

Our proposed GLR modeling method for PM_2.5 concentration estimation considering global spatial autocorrelation and local spatial heterogeneity can be divided into three steps: (1) the extraction of global spatial factors; (2) the construction of the GLR model; and (3) model assessment, comparison, and validation.

2.3.1. The Extraction of Global Spatial Factors

Global spatial factors are extracted to eliminate the effect of global spatial autocorrelation of the PM_2.5 spatial distribution. The ESF method is used to extract the global spatial factors by calculating and filtering the spatial eigenvectors of the global spatial weight matrix. The introduction of spatial eigenvectors can effectively solve the effect of spatial autocorrelation [49].

The extraction of global spatial factors can be expressed as three steps: (1) the construction of a centralized global spatial weight matrix W₀; (2) the calculation of the spatial eigenvectors of W₀; and (3) the filtering of the spatial eigenvectors.

Well-known schemes for the spatial weight matrix typically include spatially contiguous neighbors; “tri-cube” distance decline function; ranked distances; all centroids within distance d; and Gaussian distance decline [50,51]; inverse distances raised to some power. The centralized global spatial weight matrix W₀ for the stations in this study is constructed by a Gaussian kernel function:

$${\mathrm{W}}_{i,j}=\exp \left(-{\left(\frac{{\mathrm{d}}_{i,j}}{\mathrm{r}}\right)}^2\right)$$

(1)

where r is the bandwidth; ${\mathrm{d}}_{i,j}$ is the Euclidean distance between station i and station j.

After obtaining the n-dimensional centered spatial weight matrix W₀, where n is the number of monitoring stations, eigenvectors of this matrix are calculated, and only eigenvectors with eigenvalues λ_i > 0 and λ_i/λ_max ≥ 0.25 are retained [52], where λ_i is the eigenvalue of the spatial eigenvectors to be filtered and λ_max is the largest eigenvalue.

The formula for the ESF model is as follows [53]:

$$\mathrm{Y}=\mathsf{\beta }\mathrm{X}+\mathsf{\alpha }\mathrm{E}+\mathsf{\epsilon }$$

(2)

where Y is the PM_2.5 column vector, X is the matrix of the dependent variables, $\mathsf{\beta }$ is the corresponding coefficients matrix, E is the matrix of filtered spatial eigenvectors, and $\mathsf{\alpha }$ is the corresponding coefficient matrix. E is filtered out from all other candidate eigenvectors by stepwise regression [54].

2.3.2. The Construction of the GLR Model

The local spatial heterogeneity problem can be solved by the local regression model with spatially variable coefficients. The GWR method is used to obtain the spatially variable coefficients needed to construct the GLR model. Figure 2 shows the workflow of this step.

The construction of the GLR model involves (1) the introduction of global spatial factors; (2) the construction of the local spatial weight matrix; (3) the construction and judgment of the GLR model; and (4) obtaining the final GLR model.

The set of spatial eigenvectors (global spatial factors) obtained in the previous step (Section 2.3.1) need to be further filtered. In step (1), a spatial eigenvector is introduced as independent variable from the spatial eigenvectors in turn, and in step (3), we can get the temporary models with the introduction of different spatial eigenvectors. The model with the minimum Akaike information criterion (AIC) value is set as the GLR model, and the corresponding eigenvector is selected. Then, the AIC value of the current GLR model is compared with that of the previous GLR model (the first GLR model compared with the GWR model), and if the AIC value decreases, the above steps are repeated to select one more eigenvector from the remaining unselected eigenvectors and update the GLR model until the AIC value no longer decreases or all eigenvectors are selected; if the AIC value does not decrease, then the currently selected eigenvectors are the final selected eigenvectors and the current GLR model is the final GLR model.

Since the weight function that converts the distance relationships of the stations in the influence range into weight values is very sensitive to the bandwidth (the local spatial influence range), the bandwidth is continuously adjusted to obtain the optimal bandwidth before constructing the local spatial weight matrix in step (2). The most commonly used cross-validation (CV) method is applied to determine the optimal bandwidth [55,56].

The formula for the GWR can be written as:

$$y_i={\beta }_0\left(u_i,v_i\right)+{\displaystyle \sum }_{k=1}^p{\beta }_k\left(u_i,v_i\right)x_{ki}+{\epsilon }_i$$

(3)

where $y_i$ is the value of dependent variable,$x_i$ is the independent variable, ${\beta }_0$ and ${\beta }_k$ are the estimated coefficients, ε is the random error, and $\left(u_i,v_i\right)$ is the coordinate of sample point i.

The selected global spatial factors (E₁, E₂, …, E_q (q is the number of selected spatial eigenvectors)), which are filtered from the set E which is calculated in Section 2.3.1 are introduced as independent variables to construct the GLR model, and the GLR model can be expressed as:

$$y_i={\beta }_0\left(u_i,v_i\right)+{\displaystyle \sum }_{k=1}^p{\beta }_{ik}\left(u_i,v_i\right)x_{ik}+{\displaystyle \sum }_{k=1}^q{\alpha }_{ki}\left(u_i,v_i\right)E_{ki}+{\epsilon }_i$$

(4)

where $y_i$ is the value of the dependent variable at the i-th sampling point; $\left(u_i,v_i\right)$ is the coordinate of the i-th sampling point; ${\beta }_0\left(u_i,v_i\right)$ is the intercept term at this point; ${\beta }_{ik}\left(u_i,v_i\right)x_{ik}$ is the k-th independent variable and its coefficient at this point; p is the number of independent variables; ${\alpha }_{ik}\left(u_i,v_i\right)E_{ik}$ is the value of the k-th spatial eigenvector further selected from E (the E in Formula (2)) at this point and its coefficient; q is the number of spatial eigenvectors further selected from E; ${\epsilon }_i$ is the random error at this point.

In this study, the GLR model for estimating regional PM_2.5 concentrations in the YRD and BTH regions can be expressed as:

$$P{M}_{{2.5}_i}={\beta }_{0_i}+{\beta }_{1_i}AO{D}_i+{\beta }_{2_i}DE{M}_i+{\beta }_{3_i}T{S}_i+{\beta }_{4_i}P{S}_i+{\beta }_{5_i}R{H}_i+{\beta }_{6_i}PBL{H}_i+{\beta }_{7_i}NDV{I}_i+{\beta }_{8_i}ROA{D}_i+{\beta }_{9_i}FAC{T}_i+{\displaystyle \sum }_{k=1}^q{\alpha }_{k_i}{E}_{k_i}+{\epsilon }_i$$

(5)

where ${\beta }_{0_i}$ is the intercept term at point i; ${\beta }_{k_i}$(k = 1, 2 ..., 9) are regression coefficients at this point; ${\alpha }_{k_i}{E}_{k_i}$ is the value of the k-th spatial eigenvector further selected from E (the E in Formula (2)) at this point and its coefficient; q is the number of spatial eigenvectors further selected from E; ${\epsilon }_i$ is the random error at this point.

2.3.3. Model Assessment, Comparison, and Validation

To assess the performance of models and the accuracy of their fit to the true PM_2.5 concentrations, the adjusted R² (Adj.R²), AIC, and root mean square error (RMSE) are used as evaluation metrics. Moran coefficient (MC) test can be applied for detecting the effects of spatial autocorrelation.

$$\mathrm{MC}=\frac{n}{s_0}\frac{{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{w}_{ij}{Z}_i{Z}_j}{{{\displaystyle \sum }}_{i=1}^n{Z}_i^2}$$

(6)

where Z_i is the deviation between the value of geographical object i and its mean; n is the number of geographical objects; w_ij is the spatial weight between geographical objects i and j; S₀ is the total spatial weights: $S_0={{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{i=1}^n{w}_{i,j}$.

The higher the MC of PM_2.5, the stronger the spatial aggregation of PM_2.5. Additionally, if the MC for residuals are statistically insignificant or close to zero, the residuals can be considered as random errors.

A 10-fold CV is applied to validate the model’s ability to predict PM_2.5 concentrations in areas without monitoring stations. The PM_2.5 dataset is randomly divided into 10 subsets: 9 of them are used in turn for training and one for validation. The average RMSE for 10 times CV results is applied as an evaluation metric for the predictive power of models.

3. Results

3.1. Descriptive Statistics for PM_2.5 and Its Variables

What is shown in

Table 2. Summary statistics of PM_2.5 concentrations and all variables in the YRD and BTH regions (Std. Dev: Standard Deviation).

Factors	YRD				BTH
	Max	Min	Mean	Std. Dev	Max	Min	Mean	Std. Dev
PM2.5 (μg/m3)	134.80	7.00	49.78	22.11	264.96	17.99	76.04	42.34
AOD (/)	2.25	0.06	0.62	0.30	1.77	0.06	0.50	0.26
NDVI (/)	0.88	0.16	0.39	0.13	0.85	0.08	0.27	0.14
TS (K)	327.72	271.13	292.41	13.55	301.01	253.80	283.42	11.99
RH (/)	0.89	0.48	0.70	0.08	0.84	0.29	0. 53	14.08
PS (hPa)	1029.50	955.83	1010.78	13.04	1029.88	850.92	968.99	65.88
PBLH (m)	896.68	135.24	399.65	143.39	1135.44	174.19	529.17	187.07
DEM (m)	256.00	1.00	24.57	31.33	816.00	0.00	108.46	191.56
ROAD (km/km2)	3.56	0.217	1.30	0.729	5.63	0.45	1.75	1.37
FACT (count/km2)	0.69	0.00	0.21	0.14	0.54	0.03	0.26	0.13

is the summary statistics of the aforementioned variables and PM_2.5 concentrations at the monitoring stations in the YRD and BTH regions during the period from December 2015 to November 2016. The Pearson correlation coefficients shown in Tables S1 and S2 indicated that, in the YRD and BTH regions, the independent variables that were positively and significantly correlated with PM_2.5 were TS, PS, AOD, ROAD, and FACT. Additionally, the NDVI, PBLH, and DEM were significantly and negatively correlated with PM_2.5, but RH was very complexly associated with PM_2.5.

The distribution of the annual PM_2.5 observations in the YRD region and BTH region were shown in Figure 3a,b, which ranged from 22.62 to 66.75 μg/m³ in the YRD region, and ranged from 28.06 to 112.48 μg/m³ in the BTH region. The annual average PM_2.5 concentrations at the monitoring stations in the two regions were 49.61 and 75.99 μg/m³, respectively.

According to Table S3, the MCs for annual PM_2.5 were 0.605 and 0.968 in the YRD and BTH regions, respectively, indicating that there was spatial clustering of PM_2.5 in the two regions throughout the year with a strong spatial autocorrelation.

3.2. Model Results and Validation

3.2.1. The Assessment and Comparison of Models

Table 3. Model performance of the ESF, GWR and, GLR models in both regions at monthly, seasonal, and annual time scales.

Time		ESF			GWR			GLR
		Adj.R2	AIC	RMSE	Adj.R2	AIC	RMSE	Adj.R2	AIC	RMSE
YRD	Monthly average	0.574	1447.911	7.101	0.578	1404.906	5.854	0.620	1378.899	5.401
	Winter	0.686	1495.584	7.463	0.698	1466.581	7.155	0.703	1463.351	7.086
	Spring	0.449	1431.397	6.356	0.548	1354.051	5.08	0.564	1345.9	4.975
	Summer	0.382	1299.773	4.637	0.445	1240.45	3.912	0.534	1189.853	3.309
	Autumn	0.652	1262.926	4.777	0.709	1189.827	3.947	0.734	1173.29	3.796
	Annual	0.687	1193.754	4.334	0.733	1174.324	3.798	0.748	1114.99	3.427
BTH	Monthly average	0.799	528.971	7.733	0.797	504.516	6.907	0.853	467.312	4.885
	Winter	0.897	584.789	10.309	0.902	564.62	9.585	0.938	520.819	6.533
	Spring	0.77	513.331	5.973	0.758	496.641	5.685	0.809	475.005	4.748
	Summer	0.711	482.871	5.231	0.745	450.135	4.389	0.786	431.469	3.705
	Autumn	0.925	501.961	5.871	0.912	492.028	6.029	0.96	416.32	3.174
	Annual	0.926	486.544	5.362	0.943	451.56	4.437	0.959	410.66	3.138

shows the performance metrics of the monthly average, seasonal and annual models (ESF, GWR, and GLR) in the YRD and BTH regions. In the YRD region, the average Adj.R² of the monthly GLR model was 0.620, which was 8.0% and 7.4% higher than that of the ordinary monthly ESF (0.574) and GWR (0.578) models, respectively, and the average RMSE of the GLR (ESF, GWR) model was 5.401 μg/m³ (7.101 and 5.854 μg/m³). The seasonal average Adj.R² of the GLR model was 0.634, which was 16.9% and 5.7% higher than that of the ESF (0.542) and GWR (0.600) models, respectively, and the average RMSEs of the ESF and GWR models were 5.808 and 5.023 μg/m³, respectively, which were both higher than that of 4.792 μg/m³ for the GLR model. The annual Adj.R² and RMSE of the GLR model were 0.748 and 3.427 μg/m³, both better than those of the ESF and GWR models. The AIC value of the GLR model was the lowest among the three models at all time scales.

In the BTH region, the results show that the monthly average Adj.R² of the GLR model was 0.853, which was 6.8% and 7.0% higher than that of the ESF (0.799) and GWR (0.797) models, respectively, and the average RMSEs of the ESF, GWR and GLR model were 7.733, 6.907 and 4.885 μg/m³, respectively. The seasonal average Adj.R² of the GLR model was 0.873, which was 5.8% and 5.3% higher compared to that of the ESF (0.826) and GWR (0.829) models, respectively, and the seasonal average RMSE of the GLR (ESF, GWR) model was 4.540 μg/m³ (6.846 and 6.251 μg/m³). The annual Adj.R² and RMSE of the GLR model were 0.959 and 3.128 μg/m³, both better than those of the ESF and GWR models. Likewise, the AIC value of the GLR model was also the lowest among the three models at all time scales.

Tables S4 and S5 show the performance metrics of each monthly model in the YRD and BTH regions. The Adj.R² of the GLR model ranged from 0.393 to 0.719 in the YRD region, and 0.668 to 0.964 in the BTH region. The RMSE of the GLR model ranged from 2.893 to 8.237 μg/m³ in the YRD region, and 2.863 to 7.089 μg/m³ in the BTH region. The performance of all models in both regions decreases in summer.

In summary, higher average Adj.R², lower average RMSE and AIC suggested that the GLR model was superior to the ordinary GWR model and the ordinary ESF model in both regions at monthly, seasonal, and annual scales. Additionally, the GLR model had good spatial and temporal applicability.

3.2.2. The MCs for Residuals

The MCs of GLR model residuals are shown in

Table 4. The MCs for residuals of the GLR model at monthly, seasonal, and annual time scales.

Time	MC		Time	MC
	YRD	BTH		YRD	BTH
15_Dec	0.283 **	−0.371 **	16_Sep	0.132 **	−0.217 *
16_Jan	/	/	16_Oct	0.122 **	/
16_Feb	/	−0.298 **	16_Nov	/	/
16_Mar	/	/	Winter	0.238 **	−0.239 *
16_Apr	/	/	Spring	/	/
16_May	0.268 **	/	Summer	/	/
16_Jun	/	/	Autumn	/	/
16_Jul	/	−0.244 *	Annual	/	−0.218 *
16_Aug	/	/

/ insignificance; * significance at α = 0.05 level; ** significance at α = 0.01 level.

. The MCs of the GLR models were low and insignificant for most months and seasons in both regions. Overall, the residuals of the GLR model could be considered almost as random errors, indicating that the GLR model effectively weakens or even eliminates the effect of spatial autocorrelation on regression modeling.

3.2.3. Cross-Validation

A 10-fold CV was applied to evaluate the prediction accuracy of models. The RMSEs of monthly average, seasonal and, annual models are shown in

Table 5. The 10-fold CV RMSE of models at monthly, seasonal, and annual time scales.

Time	YRD			BTH
	ESF	GWR	GLR	ESF	GWR	GLR
Monthly average	7.650	7.286	7.024	10.526	11.001	9.499
Winter	7.789	7.665	6.932	8.037	8.984	8.261
Spring	6.680	6.419	6.073	7.229	6.868	6.848
Summer	4.840	4.669	4.790	5.968	5.719	5.662
Autumn	5.188	4.861	4.731	11.584	12.201	10.955
Annual	4.701	4.865	4.599	5.711	5.485	5.302

, and Table S6 shows the RMSE for the specific monthly model. The monthly average CV RMSEs of the models (ESF, GWR, and GLR) in the YRD region were 7.650, 7.286 and 7.024 μg/m³, respectively. In the BTH region, the monthly average CV RMSEs of the ESF, GWR, and GLR models were 10.526, 11.001 and 9.499 μg/m³. Meanwhile, the seasonal average CV RMSEs of the GLR (ESF, GWR) model were 5.63 μg/m³ (6.124 and 5.904 μg/m³) in the YRD region and 7.932 μg/m³ (8.205 and 8.443 μg/m³) in the BTH region. Except in summer, the RMSE of the GLR model was the lowest among the three models. Overall, the monthly, seasonal, and annual GLR models with lower RMSE had higher prediction accuracy than the other two models.

3.3. PM_2.5 Distribution Maps and Spatiotemporal Characteristics

3.3.1. Continuous PM_2.5 Distribution Maps

The coefficients were calculated by the GLR model, and by ordinary kriging could we obtain the coefficients for the entire YRD and BTH regions. The filtered spatial eigenvectors could also be interpolated by ordinary kriging into a raster image covering the entire region. Subsequently, we obtained PM_2.5 concentration maps throughout the YRD and BTH regions in different periods. The estimated monthly PM_2.5 results for the YRD and BTH regions are shown in Figures S1 and S2, respectively. The seasonal and annual PM_2.5 results for the YRD and BTH regions were shown in Figure 4 and Figure 5. In addition, for more realistic results, we set the negative values to 0 μg/m³ and removed all outliers. The missing parts of the PM_2.5 distribution maps caused by the missing values of the AOD images were obtained by inverse distance-weighted interpolation.

3.3.2. Spatiotemporal Distribution Based on the PM_2.5 Distribution Maps

Figure 6 shows the monthly PM_2.5 concentrations in the two regions obtained from the PM_2.5 estimates distribution maps. The monthly PM_2.5 concentrations in the two regions varied in a U-shape, with higher PM_2.5 concentrations in the winter months and lower PM_2.5 concentrations in the summer months. The most and least polluted months in the YRD were December 2015 and July 2016, when the average PM_2.5 concentrations reached 71.95 and 25.34 μg/m³, respectively. In the BTH region, the air was most polluted in December 2015, when the average PM_2.5 concentrations reached 139.64 μg/m³, and the cleanest month was August 2016, with an average PM_2.5 concentrations of 37.79 μg/m³.

As is shown in

Table 6. Estimated seasonal and annual PM_2.5 concentrations in the YRD and BTH regions.

Time	YRD		BTH
	Max	Mean	Max	Mean
Winter	96.27	66.86	153.36	94.91
Spring	74.03	45.68	76.96	61.21
Summer	48.83	26.77	83.55	48.94
Autumn	60.82	36.08	110.34	75.68
Annual	62.08	44.85	114.74	80.21

, the annual average PM_2.5 concentrations in the YRD and BTH regions were 44.85 and 80.21 μg/m³, respectively, and the annual average PM_2.5 concentration in the BTH region was almost twice as high as that in the YRD region. The annual average PM_2.5 concentrations in the YRD region and the BTH region exceed the WHO interim target-1 (35 μg/m³) by 28% and 129%, respectively.

According to the seasonal and annual maps, PM_2.5 pollution in both regions is significantly worse in winter than in the rest of the year. PM_2.5 pollution levels also varied by province and municipality. The order of PM_2.5 pollution from high to low in the YRD region is Anhui Province > Zhejiang Province > Jiangsu Province > Shanghai. Additionally, in the BTH region, Hebei province > Tianjin ≈ Beijing. The annual maps show that the most polluted areas of PM_2.5 in the YRD region were mainly located in the central part of Hefei in Anhui province and the western part of Wuxi and the eastern part of Changzhou in Jiangsu Province, while the most polluted areas of PM_2.5 in the BTH region were mainly located in the western part of Hengshui, the central part of Langfang and the eastern part of Handan and Shijiazhuang in Hebei Province. The annual PM_2.5 concentrations in the YRD region was higher in the northern part and lower in the southern part, and the PM_2.5 concentrations gradually decreased from the northwest to the southeast, while the PM_2.5 concentrations in the southern part of the BTH region was higher than that in the northern part.

4. Discussion

4.1. Method Improvement and Accuracy Enhancement

In this study, we proposed a GLR model considering global spatial autocorrelation and local spatial heterogeneity for PM_2.5 concentration estimation. Comprehensive results show that the GLR model outperformed the ordinary ESF and GWR models at all time scales and in two regions. This proves that our proposed model is a feasible and reasonable improvement to the ESF and GWR models. Compared with other models used in the YRD [49,57] and BTH [44] region, the GLR model still has good performance, and the average Adj.R² of the GLR model in two regions is also comparable to that of previous studies, or even improved. Further, the distribution map of PM_2.5 obtained from the GLR model is also similar to the maps in the above studies. Overall, the GLR model is an excellent model for estimating regional PM_2.5 concentrations, with good spatial and temporal applicability.

The spatial distribution of PM_2.5 has global spatial autocorrelation and local spatial heterogeneity, which are not fully considered by the ESF and GWR models. The ESF model only solves the problem of spatial autocorrelation, and can only obtain coefficients that are optimal for the global, with low accuracy in the local; considering the spatial heterogeneity, the GWR model overfocuses on local features and suffers from local overfitting. However, considering that the spatial distribution of PM_2.5 is spatially influenced by global (spatial autocorrelation) and local (spatial heterogeneity) factors, and the local spatial influence can be seen partly as the local expression of the global spatial influence and partly as the influence of the local environment itself, the GLR model solves the problems of spatial autocorrelation and spatial heterogeneity by introducing the spatial eigenvectors (global factors) in local modeling, making the model more consistent with the spatial pattern of the spatial distribution of PM_2.5. Therefore, the GLR model has higher accuracy than the ESF and GWR models and better prediction ability in areas without stations.

The climatic conditions in summer in both regions are very complex, leading to more complex non-linear relationships between PM_2.5 and factors, which may also be responsible for the reduced and fluctuating performance of all models during summer time in both regions, and this result is consistent with previous researches [33]. In the BTH region, there are fewer PM_2.5 monitoring stations, only approximately one-third of those in the YRD region, and they are extremely unevenly distributed. Therefore, when the 10-fold CV was performed, the impact of a few stations not participating in the modeling had a greater impact on the CV results, resulting in higher CV RMSEs in the BTH region than in the YRD region. Additionally, the distribution of stations in the BTH region is more concentrated and the spatial autocorrelation is stronger, resulting in a higher MC for the residuals than in the YRD region.

4.2. Comparison of PM_2.5 in the YRD and BTH Regions

The monthly, seasonal, and annual model results show that the model accuracy is lower in the YRD region than in the BTH region, and there may be some reasons that can explain this. The spatial resolution of the AOD product used in the YRD region is 3 km, while the spatial resolution of the AOD product used in the BTH region is 1 km. Coarser products may cause errors between the data used for modeling and the true values, resulting in reduced model performance. Further, the YRD region is larger than the BTH region, and its spatial and climatic conditions are much more complex than those in the BTH region, which leads to a more complex non-linear relationship between PM_2.5 and various factors in the YRD region. Therefore, the linear model is not able to estimate PM_2.5 concentrations well. Additionally, PM_2.5 pollution in the BTH region is more serious than that in the YRD region, and the PM_2.5 distribution is more aggregated. The spatial distribution of PM_2.5 is more stable, making it easier to estimate the PM_2.5 concentrations.

The YRD and BTH regions have hot and rainy summers with moist air, which is easily absorbed by particles that then gather and settle. High temperature improves the air diffusion efficiency and accelerates the diffusion of PM_2.5, while more rainfall intensifies the flushing of PM_2.5 [58]. Moreover, the light, moist air in summer is suitable for vegetation growth, which accelerates the absorption of PM_2.5. Therefore, in both regions, the summer is the best season for air quality, while the cold and dry winter is a period of severe PM_2.5 pollution, and this finding is consistent with previous studies [59,60]. PM_2.5 pollution in the BTH region is much more severe than that in the YRD region. Located at a higher latitude, the BTH region has a lower annual temperature and less precipitation than the YRD region. In addition, the BTH region is more industrialized, with a higher number of factories (37,135) than the YRD region (23,299), and coal dust from winter heating in the BTH region worsens the air [61].

According to the Pearson correlation coefficient, RH is very complexly associated with PM_2.5 in both regions, which is consistent with previous studies [33,49]. Low RH may increase the PM_2.5 concentrations due to the hygroscopic growth, while high RH may contribute to the settlement of fine particles, decreasing the PM_2.5 concentration at ground level [60,62]. Comparing the DEM (Figure S3) and factory data (Figure S4) with the PM_2.5 distribution maps of the two regions, we find that PM_2.5 pollution is worse at lower altitudes than at higher altitudes. Atmospheric pressure and temperature in high-altitude areas are lower than those in low-altitude areas. However, strong winds at high altitudes accelerate PM_2.5 diffusion [63], and people are more likely to build factories and road networks on level areas (plains), resulting in higher PM_2.5 concentrations at low altitudes.

4.3. Limitations and Future Works

The period of the data applied in this experiment is only one year, from December 2015 to November 2016, making the validation of the proposed method relatively weak on the annual time scale. More and better quality data are needed, and there is a serious lack of AOD data, a better method for processing AOD data are needed. Additionally, PM_2.5 monitoring stations with any null data were deleted, but the omitted PM_2.5 observations affect the construction of the global and local spatial weight matrix, which has a large impact on the extraction of spatial eigenvectors, thus affecting the model’s accuracy [57].

More data will be collected in a follow-up study, and we will address the above problems in follow-up research. In addition, although deep learning methods are less interpretable than linear models, they have higher accuracy and can explore the nonlinear relationship between factors and PM_2.5 [64,65]. We will try to combine linear models with deep learning methods to retain the interpretability of the models and improve the accuracy in the future.

5. Conclusions

In this paper, a GLR model considering global spatial autocorrelation and local spatial heterogeneity has been proposed for estimating ground PM_2.5 concentrations in the BTH region and the YRD region of China. Comprehensive results show that the GLR model is superior to the ordinary GWR model and ordinary ESF model in terms of Adj.R², RMSE, AIC, and CV RMSE at monthly, seasonal, and annual scales. The GLR model is able to address the effects of spatial autocorrelation and spatial heterogeneity on modeling, and to some extent, it overcomes the shortcoming of the ordinary GWR model that overfocuses on local features and the disadvantage of the poor local performance of the ordinary ESF model. Further, lower CV RMSE indicates that the GLR model has higher prediction accuracy in areas without monitoring stations. According to the annual map, annual average PM_2.5 concentrations in the YRD region and the BTH region exceed the WHO interim target-1 (35 μg/m³) by 28% and 129%, respectively. Overall, the GLR model with good spatial and temporal applicability is a promising method for estimating ground PM_2.5 concentrations.

As discussed in the 4.3 “Limitations and future works” section, more and better quality data will be tried in the future, and the linear and non-linear relationships between PM_2.5 and factor will be taken into account in the next step.