Spatiotemporal Heterogeneity and the Key Influencing Factors of PM2.5 and PM10 in Heilongjiang, China from 2014 to 2018

Particulate matter (PM) degrades air quality and negatively impacts human health. The spatial–temporal heterogeneity of PM (PM2.5 and PM10) concentration in Heilongjiang Province during 2014–2018 and the key impacting factors were investigated based on principal component analysis-based ordinary least square regression (PCA-OLS), PCA-based geographically weighted regression (PCA-GWR), PCA-based temporally weighted regression (PCA-TWR), and PCA-based geographically and temporally weighted regression (PCA-GTWR). Results showed that six principal components represented the temperature, wind speed, air pressure, atmospheric pollution, humidity, and vegetation cover factor, respectively, contributing 87% of original variables. All the local models (PCA-GWR, PCA-TWR, and PCA-GTWR) were superior to the global model (PCA-OLS), and PCA-GTWR has the best performance. PM had greater temporal than spatial heterogeneity due to seasonal periodicity. Air pollutants (i.e., SO2, NO2, and CO) and pressure were promoted whereas temperature, wind speed, and vegetation cover inhibited the PM concentration. The downward trend of annual PM concentration is obvious, especially after 2017, and the hot spot gradually changed from southwestern to southeastern cities. This study laid the foundation for precise local government prevention and control by addressing both excessive effect factors (i.e., meteorological factors, air pollutants, vegetation cover) and spatial-temporal heterogeneity of PM.


Introduction
Particulate matter (PM) is dispersed throughout the atmosphere, decreasing visibility, and interfering with plant photosynthesis [1][2][3][4]. In addition, numerous studies have demonstrated that airborne particulate matter has a negative impact on human health, increasing the chance of children's breathing problems [5]. It has been established that PM 2.5 and PM 10 (atmospheric particles with aerodynamic diameters of 2.5 µm and 10 µm, respectively) are associated with an increased risk of human mortality [6]. According to the Global Burden of Disease 2010 comparative risk assessment (GBD), ambient air pollution from PM 2.5 was ranked as the sixth highest overall risk factor for worldwide premature death [7].
According to the Global Urban Air Quality Index (AQI) Report, with the rapid urbanization and industrialization over the past three decades, Chinese cities have suffered serious air pollution challenges [8]. In January 2013, smog blanketed more than 1.3 million km 2 of China, impacting approximately 850 million people [9]. In the same dimensions [40,41]. A PCA cannot only reduce the dimensionality but also retain more information. Its purpose is to construct a set of new orthogonal variables known as the principal components for replacing the original high-dimensional dataset [42][43][44][45].
This study integrates a principal component analysis (PCA) into the local models to efficiently investigate the temporal and spatial heterogeneity of PM (PM 2.5 and PM 10 ) and explains the key factors that influence it simultaneously. Specifically, the objectives include the following: (1) to conduct a PCA on the most relevant variables (including air pollutants, meteorological factors, and vegetation factors); (2) to establish a principal component analysis-based geographically weighted regression (PCA-GWR), a temporally weighted regression (PCA-TWR), and a geographically and temporally weighted regression (PCA-GTWR) and compare them to the corresponding global model, a principal component analysis-based ordinary least square regression (PCA-OLS); and (3) to explore the spatiotemporal characteristics of PM and the key influencing factors on them in Heilongjiang Province in recent years. This research provides a scientific foundation and technical support for a better understanding of the temporal and spatial heterogeneities of PM as well as the factors that drive PM, all of which are critical for future PM prevention and control.

Study Area
The present study was conducted in Heilongjiang Province, China, which is located in the northernmost region of China, within the longitudes of 121 • 11 E to 135 • 05 E and the latitudes of 43 • 26 N to 53 • 33 N (Figure 1). Heilongjiang has a mountainous topography in the northwest, north, and southeast and flat land in the northeast and southwest [4]. The climate of Heilongjiang Province is a cold temperate zone in continental monsoon climate. The major climatic features include warm and dry springs, hot and humid summers, dry autumns, and cold and lengthy winters, resulting in a protracted heating season. With a forest area of 21.47 million hectares and a forest coverage of 47.3%, Heilongjiang Province has a considerable forest area. The forest coverage rate of the different cities varies greatly, with some cities having greater forest coverage rates, such as Da Xing'an Mountain, Heihe, Yichun, and Mudanjiang, while others have lower forest coverage rates [46].

Air Pollutants Data
For this study, daily records of 6-criterion air pollutants from 2014 to 2018, including SO 2 , NO 2 , PM 10 , CO, O 3 , and PM 2.5 , were collected by the Environmental Monitoring Station of Heilongjiang Province based on the 57 environmental monitoring sites in Heilongjiang Province. The details of measuring air pollutants was described in Ambient air quality standards of P.R. China [47]. Apart from CO (mg/m 3 ), all of the pollutants were measured in µg/m 3 .

Meteorological Data
Sixteen meteorological variables were received from the Resource and Environmental Science and Data Center's (http://www.resdc.cn/, (accessed on 12 September 2022)) daily data collection of meteorological station observations during the same time period. These included the average pressure (hPa), the maximum pressure (hPa), the minimum pressure (hPa), the average temperature ( • C), the daily maximum temperature ( • C), the daily minimum temperature ( • C), the average relative humidity, the cumulative precipitation from 8 pm to the next-day 8 pm (daily cumulative precipitation, mm), the average wind speed (m/s), the maximum wind speed (m/s), the extreme wind speed (m/s), the sunshine hours (h), the average surface temperature ( • C), the daily maximum surface temperature ( • C), and the daily minimum surface temperature ( • C). Maximum wind speed is the maximum of average 10-min wind speed in a given time period. Extreme wind speed is the maximum instantaneous wind speed in a given period of time.

Air Pollutants Data
For this study, daily records of 6-criterion air pollutants from 2014 to 2018, including SO2, NO2, PM10, CO, O3, and PM2.5, were collected by the Environmental Monitoring Sta tion of Heilongjiang Province based on the 57 environmental monitoring sites in Hei longjiang Province. The details of measuring air pollutants was described in Ambient ai quality standards of P.R. China [47]. Apart from CO (mg/m 3 ), all of the pollutants were measured in μg/m 3 .

Meteorological Data
Sixteen meteorological variables were received from the Resource and Environmen tal Science and Data Center's (http://www.resdc.cn/, (accessed on 12 September 2022) daily data collection of meteorological station observations during the same time period These included the average pressure (hPa), the maximum pressure (hPa), the minimum pressure (hPa), the average temperature (°C), the daily maximum temperature (°C), the daily minimum temperature (°C), the average relative humidity, the cumulative precipi tation from 8 pm to the next-day 8 pm (daily cumulative precipitation, mm), the average wind speed (m/s), the maximum wind speed (m/s), the extreme wind speed (m/s), the sunshine hours (h), the average surface temperature (°C), the daily maximum surface tem perature (°C), and the daily minimum surface temperature (°C). Maximum wind speed i the maximum of average 10-min wind speed in a given time period. Extreme wind speed is the maximum instantaneous wind speed in a given period of time.

MODIS NDVI
The NASA Earth Observing System (EOS) Terra and Aqua satellites are equipped with the moderate resolution imaging spectroradiometer (MODIS) [48]. It includes two types of satellites (i.e., Terra and Aqua) and the Aqua satellite was used in this study.
The less cloudy MODIS L1B data of Heilongjiang Province with spatial resolutions of 1 km from 2014 to 2018 were selected and preprocessed (https://ladsweb.modaps. eosdis.nasa.gov/ (accessed on 12 September 2022)). First, radiometric calibration and a geometric correction were performed to correct and eliminate the "bow" effect and other deformations [49,50]. Then, the atmospheric correction of MODIS was performed and the atmospheric correction used the simplified dark pixel method [51]. Finally, the normalized difference vegetation index (NDVI), which captures the forest ecosystem's structural information, was extracted using the binary classification approach. The NDVI equation is as follows [52]: where IR and R is the pixel value in the infrared band and the red band, respectively. The NDVI data were integrated with the air pollutants and meteorological data of 57 sites (i.e., 20 variables per site, one or two NDVI data per site per month, and 5945 records in total). The statistical indicators of the explanatory factors and the dependent variables are shown in Table 1.

Methods
To reduce the dimensionality and multicollinearity of variables, a principal component analysis (PCA) was conducted on the most relevant variables (including air pollutants, meteorological factors, and vegetation factors). Based on PCA, the global model (PCA-OLS) and local models (PCA-GWR, PCA-GTWR, and PCA-TWR) were established.

Principal Component Analysis
Principal component analysis is a statistical analysis method that transforms multiple variables into a set of mutually orthogonal vectors (i.e., principal components). These principal components (PCs) retain most of the information of the original variables, which are usually expressed as the linear combination of the original variables. The relationship between principal component y i and original variable x i is as follows [42,53,54]: Among them, y i is the i th PC, (x 1 , x 2 , x 3 , . . . , x p ) are the original variables, a ij represents the linear correlation coefficient for the i th principal component and the j th original variable, which is also known as the loading of the variable on the common factor. The PCs are independent of one another [55], so that: Cov y i , y j = 0 (i = j, i, j = 1, 2, . . . , p) Among them, y i and y j are the two PCs, and Cov y i , y j is the covariance between them. The variance of PCs decreases in turn, and the first PC contains the most information. the cumulative variance contribution rate is as follows: Among them, λ k (k = 1, 2, . . . , p) is the variance of the k th PC, ∑ p i=1 λ is the sum of variances for p variables.
The cumulative variance contribution rate of the first m variables is as follows: Among them, ∑ m i=1 λ is the sum of variances for the first m variables. When b reaches a certain value (in this study, b is not less than 85%), the first m components are selected as PCs. To make the PCs interpretable, they are rotated using the maximum variance orthogonal rotation approach.
Prior to PCA, the Kaiser-Meyer-Olikin (KMO) and Bartlett sphericity tests was used to determine whether the commonality of the variables was high. If KMO > 0.5 and p < 0.01, the variables are suitable for PCA [54].

Global Model (PCA-OLS)
The principal component analysis-based ordinary least square regression (PCA-OLS) is a global model used to describe the linear relationship between the response and independent variables [56], that is, PCs in this study. Its formula is as follows: where Y i is the response variable (i.e., PM 2.5 , PM 10 in this study), X k (k = 1, 2, . . . , p − 1) represents the effective or selected PCs. p is the total number of parameters to be estimated and n is the number of samples (in this study, n = 5945). β k is the model parameter, ε i is the model error term, with an expected value of zero and a normal distribution. Since PCA-OLS is a global model, model parameters β k are estimated by using data of the entire study area and the model coefficient vector where X and Y are the vectors of PCs and PM, the superscript T denotes the transpose of a matrix.

Local Models PCA-GWR Model
To explore geographical nonstationarity, principal component analysis-based geographically weighted regression (PCA-GWR) is conducted in this study. PCA-GWR extends the PCA-OLS regression framework as follows [57]: where Y i is the response variable (i.e., PM 2.5 , PM 10 in this study), are the intercept and a set of p − 1 slope parameters at the i th observation. X ik (k = 1, 2,..., p − 1) are a set of p − 1 PC at the i th location, p is the total number of parameters to be estimated, ε i is the error term of i th observation.
β(u i , v i ) could be estimated by the local least square method as follows [58]: in which the off-diagonal components are zero and the diagonal elements represent the geographic weights at the observation i, and the matrix is as follows: Among them, w ij is spatial weight and it is determined by the spatial kernel function, also called a distance-decay function.
A spatial kernel function is often classified into two types: fixed kernels and adaptive kernels [34]. The adaptive bisquare function is utilized to produce geographic weights in this study and the formula is as follows [59]: where h i is a non-negative metric known as bandwidth, d ij is the distance between locations i and j and it is calculated as follows: As the distance between the two locations increases, the spatial effect between the two gradually decays and disappear beyond the bandwidth, that is, w ij = 0 [58]. The Cross-validation (CV) score and Akaike Information Criterion (AICc) criterion are standard metrics for determining optimal bandwidth [60,61].

PCA-GTWR and PCA-TWR Model
While temporal variability is considered, principal component analysis-based geographically and temporally weighted regression (PCA-GTWR) is conducted with a weight matrix among i and other points is calculated in spatial and temporal space. The PCA-GTWR model can be expressed as follows [34]: where Y i is the response variable (i.e., PM 2.5 , PM 10 in this study), are the intercept and the parameter at the i th observation with the spatio-temporal coordi- The estimation of β(u i , v i , t i ) is similar to that in PCA-GWR (Equation (8)).
Because the spatio-temporal dimension is taken into account, a spatial-temporal weight matrix is constructed based on a spatio-temporal distance. In this study, spatio-temporal distance is as follows: where λ and µ are scaling factors used to balance the different effects of space and time, In this study, an adaptive Gaussian distance-decay function was applied to construct the spatial-temporal weight matrix, and the formula is as follows: where h 2 ST is a parameter of spatio-temporal bandwidth, and h 2 S and h 2 T are spatial and temporal bandwidth. When geographical variation is neglected (λ = 0), Equation (15) is changed as follows: where d T ij is the temporal distance. Based on d T ij , principal component analysis-based temporally weighted regression (PCA-TWR) is conducted [4]. PCA-TWR and PCA-GWR are special cases of PCA-GTWR without considering either spatial or temporal variation.

Model Assessment
Akaike's information criterion (AICc), adjusted coefficient of determination (R 2 a ), root mean squared errors (RMSE), and mean absolute error (MAE) are used to assess the model performance. The formula of AICc is shown as follows: Among them, n is the number of samples,σ is the estimated standard deviation of the error term, and tr(S) is the trace of hat matrix S(i.e.,Ŷ = SY). In the GWR model, the hat matrix is as follows: R 2 a will not exaggerate the explained percentage compared to the coefficient of determination (R 2 ), and the formula is as follows: where n is the number of the sample, and p is the number of coefficients (including all predictor coefficients and the intercept). The root mean squared errors (RMSE) and mean absolute errors (MAE) are also used in this study, and the formula is as follows: where y i is the observed values of the response variable andŷ i is the predicted values of the response variable.

Principal Component Analysis
The Kaiser-Meyer-Olikin (KMO) and Bartlett sphericity tests were conducted prior to the PCA, and the results (KMO = 0.76 > 0.5 and p < 0.01) revealed that the original 20 variables were highly correlated and PCA was appropriate. As shown in Table 2, the first six PCs were chosen due to the cumulative variance contribution rate of 87%. According to the principal component loadings, PC1~PC6 could represent the temperature, wind speed, air pressure, atmospheric pollution (SO 2 , NO 2 , CO), humidity, and vegetation cover from the original dataset, respectively.

Globel Model (PCA-OLS)
The stepwise PCA-OLS model for PM, a global model, was applied as a baseline for model comparison in this study. According to Table 3, the PM was negatively correlated with temperature (PC1), wind speed (PC2), and vegetation (PC6). PC5 was excluded for PM 10 due to its nonsignificance. PM 2.5 was positively correlated with pressure (PC3), atmospheric pollutants (PC4), and humidity (PC5), whereas PM 10 was only positively correlated with pressure (PC3) and pollutants (PC4). The PCA-OLS models modestly fitted the data according to the R 2 a , i.e., 61% and 54% of the total variation of PM 2.5 and PM 10 can be explained by the PCA-OLS models.

Local Models
To characterize the spatial and temporal heterogeneities, local models (i.e., PCA-GWR, PCA-TWR, and PCA-GTWR model) were utilized. The adaptive bisquare function was used, and AICc was employed to select a number of neighbors. Tables 4 and 5 summarize the statistics of the parameter estimates of the PCA-GWR, PCA-TWR, and PCA-GTWR for PM 2.5 and PM 10 , respectively. The sign of the global model parameters matches the average value of the local model parameters for both PM 2.5 and PM 10 . It indicates that the local models follow the same trend as the corresponding global model: PM positively correlated with air pressure, other atmospheric pollutants, and humidity, while negatively correlated with temperature, wind, and vegetation cover. The number of neighbors of PCA-GWR is roughly twice that of the PCA-GTWR model, whereas that between PCA-GWR and PCA-GTWR is quite similar (i.e., around 600) for both PM 2.5 and PM 10 . After considering spatial and temporal variation, PCA-GTWR performed better than PCA-GWR and PCA-TWR for PM according to the goodness-of-fit (i.e., R 2 a and AICc). However, the local model of PM 2.5 is superior (i.e., higher R 2 a ) than the corresponding model of PM 10 .  Table 6 depicts the goodness-of-fit of the PCA-based local and global models. No matter if it was the PM 2.5 or PM 10 , PCA-GTWR performed the best (largest R 2 a , smallest AICc, RMSE and MAE), followed by PCA-GWR, PCA-TWR, and PCA-OLS. The PCA-TWR is much superior to the PCA-GWR, with a greater R 2 a , a lower AICc, RMSE, and MAE. This demonstrates that the PM's temporal nonstationarity is more noticeable than its spatial nonstationarity.

Spatial Characteristics of PCA-GTWR Parameter Estimates for PM
According to Figure 2, each PC factor had a difference influence in various locations for PM 2.5 . Temperature (Figure 2a) was the inhibitor of PM 2.5 in most of Heilongjiang province. The inhibition was strong and obvious in the eastern area but weak in the west (e.g., Da Xing'an Mountain, Heihe, Qiqihar, Daqing, and Harbin). Wind speed (Figure 2b) primarily suppressed PM 2.5 , and the total inhibitory impact was comparable to temperature. Air pressure (Figure 2c) mostly promoted PM 2.5 , with a relative weak influence in Daqing, Harbin, Hegang, Jiamusi, and Shuangyashan. Atmospheric pollutants (Figure 2d) positively correlated with PM 2.5 , with three strata of correlations. The largest impact was observed the west of Harbin and the southern part around Daqing and Shuihua, whereas the smallest impact was observed in Da Xing'an Mountain, Heihe, and Qiqiha'er, and other places falling somewhere in the middle. Humidity (Figure 2e) had the ability to both enhance and inhibit PM 2.5 levels. In the western cities, especially Harbin, Daqing, Suihua, and Qiqiha'er, humidity plays a role of promoting PM 2.5 ; whereas in the eastern cities (e.g., Yichun, Hegang, Jiamusi, Shuangyashan, Qitaihe, and Jixi), it had a substantial negative effect. Vegetation cover (Figure 2f) obviously restricted PM 2.5 concentration. The restriction trend decreased towards the east, with the smallest impact in Hegang, Jiamusi, and Shuangyashan.
According to Figure 2, each PC factor had a difference influence in various locations for PM2.5. Temperature (Figure 2a) was the inhibitor of PM2.5 in most of Heilongjiang province. The inhibition was strong and obvious in the eastern area but weak in the west (e.g., Da Xing'an Mountain, Heihe, Qiqihar, Daqing, and Harbin). Wind speed (Figure 2b) primarily suppressed PM2.5, and the total inhibitory impact was comparable to temperature. Air pressure (Figure 2c) mostly promoted PM2.5, with a relative weak influence in Daqing, Harbin, Hegang, Jiamusi, and Shuangyashan. Atmospheric pollutants (Figure 2d) positively correlated with PM2.5, with three strata of correlations. The largest impact was observed the west of Harbin and the southern part around Daqing and Shuihua, whereas the smallest impact was observed in Da Xing'an Mountain, Heihe, and Qiqiha'er, and other places falling somewhere in the middle. Humidity (Figure 2e) had the ability to both enhance and inhibit PM2.5 levels. In the western cities, especially Harbin, Daqing, Suihua, and Qiqiha'er, humidity plays a role of promoting PM2.5; whereas in the eastern cities (e.g., Yichun, Hegang, Jiamusi, Shuangyashan, Qitaihe, and Jixi), it had a substantial negative effect. Vegetation cover (Figure 2f) obviously restricted PM2.5 concentration. The restriction trend decreased towards the east, with the smallest impact in Hegang, Jiamusi, and Shuangyashan. Similar to PM2.5, different PC factors had various influences on PM10 in space ( Figure  3). However, the magnitude of parameter estimates for PM10 was much larger than that Similar to PM 2.5 , different PC factors had various influences on PM 10 in space ( Figure 3). However, the magnitude of parameter estimates for PM 10 was much larger than that for PM 2.5 . Temperature (Figure 3a) could both enhance and inhibit PM 10 . Temperature tended to promote PM 10 in Da Xing'an Mountain, Harbin, and Heihe. Whereas in other places, it played an inhibitory role. The inhibitory impact was stronger in the eastern area (i.e., Hegang, Jiamusi, Shuangyashan, Jixi) than the western area (e.g., Qiqiha'er, Daqing, Suihua, and Yichun). Wind speed (Figure 3b) had an inhibitory influence on PM 10 , and the inhibitory effect was fairly uniform across the region. Both pressure (Figure 3c) and atmospheric pollutants (Figure 3d) positively correlated with PM 10 . The positive correlation between PM 10 and pressure tended to decrease from south to north, while its correlation with atmospheric pollutants presented the opposite trend. Vegetation cover (Figure 3e) inhibited PM 10 , with strong suppression in the cities of Qiqihar, Daqing, Suihua, Harbin, and Mudanjiang, and the suppression became weaker towards the east.
Suihua, and Yichun). Wind speed (Figure 3b) had an inhibitory influence on PM10, and the inhibitory effect was fairly uniform across the region. Both pressure (Figure 3c) and atmospheric pollutants (Figure 3d) positively correlated with PM10. The positive correlation between PM10 and pressure tended to decrease from south to north, while its correlation with atmospheric pollutants presented the opposite trend. Vegetation cover (Figure 3e) inhibited PM10, with strong suppression in the cities of Qiqihar, Daqing, Suihua, Harbin, and Mudanjiang, and the suppression became weaker towards the east.  Figure 4 shows that influence of each PC on PM has temporal non-stationarity and periodicity, particularly PC3 (air pressure) and PC4 (atmospheric pollutants). For PM, the GTWR parameter estimates of PC3 and PC4 in the heating (i.e., winter) and non-heating (i.e., spring, summer, autumn seasons) periods differed and periodically changed during the course of a year: high in the heating season (usually from October to April next year) and low in the non-heating seasons. This is because the higher concentrations of atmospheric pressure and other air pollutants in winter are correlated with higher PM concentrations than those in summer.

Temporal Characteristics of PCA-GTWR Parameter Estimates for PM
The parameter estimates of PC1 (temperature) and PC2 (wind speed) has larger temporal periodicity for PM10 than PM2.5. The inhibitory effect of PC2 (wind speed) in autumn (September, October, and November) tends to be larger than that in other seasons. Due to  Figure 4 shows that influence of each PC on PM has temporal non-stationarity and periodicity, particularly PC3 (air pressure) and PC4 (atmospheric pollutants). For PM, the GTWR parameter estimates of PC3 and PC4 in the heating (i.e., winter) and non-heating (i.e., spring, summer, autumn seasons) periods differed and periodically changed during the course of a year: high in the heating season (usually from October to April next year) and low in the non-heating seasons. This is because the higher concentrations of atmospheric pressure and other air pollutants in winter are correlated with higher PM concentrations than those in summer.

Temporal Characteristics of PCA-GTWR Parameter Estimates for PM
The parameter estimates of PC1 (temperature) and PC2 (wind speed) has larger temporal periodicity for PM 10 than PM 2.5 . The inhibitory effect of PC2 (wind speed) in autumn (September, October, and November) tends to be larger than that in other seasons. Due to the climate characteristics of Heilongjiang Province (cold temperate, and a temperate continental monsoon climate), the high wind speed in autumn tends to increase the inhibitory effect on PM. The inhibitory effect PC6 (vegetation coverage) on PM was obvious in 2014 but decrease and became stable afterwards. The impact of PC5 (humidity) on PM 2.5 is not stable because PC5 is comprehensive factor that mainly controlled by daily average relative humility, daily cumulative humility and daily sun hours.
the climate characteristics of Heilongjiang Province (cold temperate, and a temperate continental monsoon climate), the high wind speed in autumn tends to increase the inhibitory effect on PM. The inhibitory effect PC6 (vegetation coverage) on PM was obvious in 2014 but decrease and became stable afterwards. The impact of PC5 (humidity) on PM2.5 is not stable because PC5 is comprehensive factor that mainly controlled by daily average relative humility, daily cumulative humility and daily sun hours.

The Key Influencing Factors of PM
Recently, many studies have investigated the connections between PM and a number of variables, including air pollutants, meteorological variables (such as temperature, wind speed, and relative humidity), the normalized difference vegetation index (NDVI) derived from satellite imagery, aerosol optical depth (AOD), human activities (such as transportation emission variables, the density of industrial plants, land use, GDP), and topographical factors and more [62][63][64][65][66][67]. Although it tends to improve PM predictions using more factors, it is difficult to remove multicollinearity among these factors and interpret the spatio-temporal relationship between PM and a number of factors. In this study, PCA was used to reduce the dimension of independent variables (i.e., from 20 original variables to 6 PCs) and remove the multicollinearity, and then OLS, GWR, TWR, and GTWR models were conducted based on the PCs to predict and explore the temporal and spatial heterogeneity of PM in Heilongjiang Province. Air pollutants (mainly SO2, NO2, and CO) represented by PC4 in this study have the largest influence on both PM2.5 and PM10 (see Tables  2 and 3). Temperature expressed by PC1 had the second greatest impact on PM2.5, followed

The Key Influencing Factors of PM
Recently, many studies have investigated the connections between PM and a number of variables, including air pollutants, meteorological variables (such as temperature, wind speed, and relative humidity), the normalized difference vegetation index (NDVI) derived from satellite imagery, aerosol optical depth (AOD), human activities (such as transportation emission variables, the density of industrial plants, land use, GDP), and topographical factors and more [62][63][64][65][66][67]. Although it tends to improve PM predictions using more factors, it is difficult to remove multicollinearity among these factors and interpret the spatio-temporal relationship between PM and a number of factors. In this study, PCA was used to reduce the dimension of independent variables (i.e., from 20 original variables to 6 PCs) and remove the multicollinearity, and then OLS, GWR, TWR, and GTWR models were conducted based on the PCs to predict and explore the temporal and spatial heterogeneity of PM in Heilongjiang Province. Air pollutants (mainly SO 2 , NO 2 , and CO) represented by PC4 in this study have the largest influence on both PM 2.5 and PM 10 (see Tables 2 and 3). Temperature expressed by PC1 had the second greatest impact on PM 2.5 , followed by air pressure (PC3), wind speed (PC2), vegetation cover (PC6), and humidity factor (PC5). Whereas the impacts of temperature, pressure, and vegetation cover on PM 10 are roughly equivalent, followed by wind speed.
Although PM 2.5 and PM 10 are closely correlated [68], this study still shows some differences in the correlation between PM 2.5 /PM 10 and the influenced factors. The main distinction is that, in contrast to PM 2.5 , the humidity factor (PC5) is not significant for PM 10 and is eliminated in the global model. This may be because the humidity factor in this study is a combined variable of daily average relative humidity (positive), cumulative precipitation (positive), and sun hours (negative). It is challenging to demonstrate a meaningful relationship between the mixed variable and the comparatively large particle (PM 10 ). The second factor contributing to the difference was wind speed. The dilution effect of wind speed on PM 2.5 is stronger than that of PM 10 due to the relatively light quality, which is consistent with the previous study [68].
In general, PM is positively related to other air pollutants (mainly SO 2 , NO 2 , and CO). It is similar to the finding that 5 criteria air pollutants (i.e., PM 10 , SO 2 , NO 2 , CO, O 3 ) had a positive impact on PM 2.5 , except O 3 , in the previous study [4]. The concentration of PM reduces as the temperature (PC1) rises. When the temperature is high, the atmospheric tropospheric motion becomes more intense, resulting in the upward transport of PM, and the high temperature encourages Brownian particle motion, which is more favorable to diffusion [26]. Not surprisingly, the wind speed (PC2) has a negative effect on PM, which is consistent with previous studies [69,70]. This is because that the higher the wind speed, the stronger the particle dispersion capacity and the lower the particle concentration [71]. Atmospheric pressure (PC3) plays a positive role in increasing PM centration. When an area is subjected to high pressure, air convection is reduced, allowing contaminants to accumulate more easily, and vice versa [26,72]. The ability of the forest's complex canopy structure to absorb particulate matter has long been recognized as a critical tool for controlling PM [73,74]. In this study, vegetation cover (PC6) also has a considerable inhibitory influence on PM. The absolute magnitude of the coefficient indicates that it has a significant impact on preventing PM increases. The explanation could be that vegetation intercepts and absorbs particulate matter via Brownian diffusion, interception, and gravity deposition [75,76].

Global and Local Models
In the constructed model, the problem of the PCA-OLS model is that it uses a global model to represent the relationship between dependent and independent variables, ignoring the spatial or temporal effect of the studied variables. This spatial or temporal effect (spatial/temporal autocorrelation and heterogeneity) may violate the assumption of independent observation or constant variance, resulting in the biases of estimates of standard errors and imprecision of coefficient estimates [77,78]. PCA-GWR was better than PCA-OLS by considering spatial heterogeneity (11% and 17% improvement of adjusted R 2 of PM 2.5 and PM 10 , respectively). PCA-TWR was also superior to PCA-OLS by considering temporal heterogeneity (23% and 28% improvement of adjusted R 2 of PM 2.5 and PM 10 , respectively). The temporal heterogeneity is more significant than the spatial heterogeneity due to the obvious seasonal variation of PM, which is consistent with the previous study [4]. The PCA-GTWR model that considers both temporal and spatial heterogeneity has the best performance (compared with PCA-OLS, the adjusted R 2 increases by 25% and 31% for of PM 2.5 and PM 10 , respectively). However, the improvements of PCA-GTWR from PCA-TWR are very limited, especially for PM 2.5 . It indicates that temporal information is more effective than spatial information for modeling PM based on PCA.

The Spatial and Temporal Distribution of PM Concentrations in Heilongjiang Province
Heilongjiang Province is located in northern China and is noted for its high latitude and harsh winters. Heilongjiang Province is one of China's largest provinces, with a north-south latitude range of approximately 10 • and an east-west longitude range of approximately 4 • . According to Figures 5 and 6, the PM content in the southern area is significantly higher whereas that in the northern region is rather low due to the large population and coal burning in the south. Furthermore, the PM concentrations remained high from 2014 to 2017, especially in 2017. From spatio-temporal perspective, the hot spot of PM 2.5 concentration gradually changed from southwestern (e.g., Harbin, Daqing) to southeastern cities (e.g., Harbin, Mudanjiang, Qitaihe, and Jixi). According to Figures 5f and 6f, it can be seen that the temporal distribution of PM in Heilongjiang Province has heterogeneity and periodicity (high in winter and low in summer). The PM concentration in Heilongjiang Province grows dramatically during the heating season (from October to April next year). Coal-fired heating not only consumes energy but also harms the environment. However, due to the Heilongjiang Province's policies (Air Pollution Prevention and Control Action Plan) [10], the downward trend of annual PM concentration is obvious, especially after 2017, which indicates the progress in the prevention and control of air pollution in Heilongjiang. In the future, Heilongjiang Province should implement a multi-energy integrated heating system that makes full use of renewable energy sources such as wind, water, and solar energy. Additionally, it is critical to consider advancing science and technology, discovering methods to store renewable energy, and even storing energy in the summer for use in winter heating. Finally, increasing the rate of urban greening, focusing on forest ecosystem protection, and enhancing the ecosystem's ecological benefits are all necessary.
cantly higher whereas that in the northern region is rather low due to the large population and coal burning in the south. Furthermore, the PM concentrations remained high from 2014 to 2017, especially in 2017. From spatio-temporal perspective, the hot spot of PM2.5 concentration gradually changed from southwestern (e.g., Harbin, Daqing) to southeastern cities (e.g., Harbin, Mudanjiang, Qitaihe, and Jixi). According to Figures 5f and 6f, it can be seen that the temporal distribution of PM in Heilongjiang Province has heterogeneity and periodicity (high in winter and low in summer). The PM concentration in Heilongjiang Province grows dramatically during the heating season (from October to April next year). Coal-fired heating not only consumes energy but also harms the environment. However, due to the Heilongjiang Province's policies (Air Pollution Prevention and Control Action Plan) [10], the downward trend of annual PM concentration is obvious, especially after 2017, which indicates the progress in the prevention and control of air pollution in Heilongjiang. In the future, Heilongjiang Province should implement a multi-energy integrated heating system that makes full use of renewable energy sources such as wind, water, and solar energy. Additionally, it is critical to consider advancing science and technology, discovering methods to store renewable energy, and even storing energy in the summer for use in winter heating. Finally, increasing the rate of urban greening, focusing on forest ecosystem protection, and enhancing the ecosystem's ecological benefits are all necessary.

Conclusions
In this study, the spatial-temporal heterogeneity of PM (PM2.5 and PM10) concentration in Heilongjiang Province during 2014-2018 and the key impacting factors were investigated based on PCA-based local and global models. Six PCs with a contribution rate of 87% were selected, and each PC's (PC1-PC6) meaning was obvious after rotation, representing the temperature, wind speed, air pressure, atmospheric pollution, humidity, and vegetation cover, respectively. According to the model assessment, all the PCA-based local models (PCA-GWR, PCA-TWR, and PCA-GTWR) were superior to the PCA-based global (PCA-OLS) model, in which PCA-GTWR performed the best. Air pollutants (mainly SO2, NO2, and CO) have the largest influence on both PM2.5 and PM10. The temperature has the second greatest impact on PM2.5, followed by air pressure, wind speed, vegetation cover, and humidity factor, whereas the impacts of temperature, pressure, and vegetation cover on PM10 are roughly equivalent, followed by wind speed. In general, PM is positively related to other air pollutants (mainly SO2, NO2, and CO) and air pressure, and negatively correlated to temperature, wind speed, and vegetation cover.
This study demonstrated that the temporal heterogeneity was more pronounced than the spatial heterogeneity of PM in Heilongjiang Province. This work provides temporal and spatial heterogeneity evidence of the relationship between PM and meteorological factors, air pollutants, and vegetation cover. Furthermore, this work addresses both excessive impact factors and temporal and spatial heterogeneity and provides a theoretical foundation for precise local government prevention and control.

Conclusions
In this study, the spatial-temporal heterogeneity of PM (PM 2.5 and PM 10 ) concentration in Heilongjiang Province during 2014-2018 and the key impacting factors were investigated based on PCA-based local and global models. Six PCs with a contribution rate of 87% were selected, and each PC's (PC1-PC6) meaning was obvious after rotation, representing the temperature, wind speed, air pressure, atmospheric pollution, humidity, and vegetation cover, respectively. According to the model assessment, all the PCA-based local models (PCA-GWR, PCA-TWR, and PCA-GTWR) were superior to the PCA-based global (PCA-OLS) model, in which PCA-GTWR performed the best. Air pollutants (mainly SO 2 , NO 2 , and CO) have the largest influence on both PM 2.5 and PM 10 . The temperature has the second greatest impact on PM 2.5 , followed by air pressure, wind speed, vegetation cover, and humidity factor, whereas the impacts of temperature, pressure, and vegetation cover on PM 10 are roughly equivalent, followed by wind speed. In general, PM is positively related to other air pollutants (mainly SO 2 , NO 2 , and CO) and air pressure, and negatively correlated to temperature, wind speed, and vegetation cover.
This study demonstrated that the temporal heterogeneity was more pronounced than the spatial heterogeneity of PM in Heilongjiang Province. This work provides temporal and spatial heterogeneity evidence of the relationship between PM and meteorological factors, air pollutants, and vegetation cover. Furthermore, this work addresses both excessive impact factors and temporal and spatial heterogeneity and provides a theoretical foundation for precise local government prevention and control.