The BTH (113.45°E–119.85°E and 36.03°N–42.62°N) is one of the most important administrative, commercial, and cultural center in northern China, including Beijing and Tianjin, and 11 prefecture-level cities of Hebei Province (Figure 1). The region is densely populated and is a secondary industry that used coal as the main energy source emits various air pollutants, which causing relatively severe haze (Zhao et al., 2019). In particular, in the inland plains, coupled with unfavorable topography, it makes it more difficult for pollutants to spread (Lv et al., 2017). According to the statistics from the “China Environmental Bulletin” (http://www.cnemc.cn/jcbg/zghjzkgb/) during 2013–2020, the BTH included seven, eight, seven, six, six, five, four, and one, respectively, among the top 10 cities with poor air quality in China. Although the air quality in this region has improved during the past few years, we should still pay close attention to PM_2.5 pollution. Therefore, it is essential to analyze the spatiotemporal distribution and the trend of [PM_2.5].

In addition to adding AOD to the model for [PM_2.5] prediction, it has been recognized that combining meteorological and LU information can significantly improve the model predictability (Hu et al., 2017; Wang G. et al., 2021). In this study, for the proposed two-stage statistical regression (LME + GWR) model, a main independent predictor (AOD) and eight auxiliary predictors [i.e., planetary boundary layer height (PBLH), 2-m air temperature (TEMP), 10-m wind speed (WS), relative humidity (RH, specific humidity calculated), surface pressure (PRS), precipitation (PRCP), forest coverage (FC), and urban coverage (UC)] were utilized through variables selection and multicollinearity diagnosis. The datasets covered the period from January 1, 2013, to December 31, 2020. The detail information about the datasets is shown in Table 1.

1) PM_2.5 data. The PM_2.5 hourly concentration of 80 monitoring stations in BTH was obtained from the National Urban Air Quality Real-time Release Platform. In the process of fitting the daily mean [PM_2.5], we eliminated the [PM_2.5] (i.e., <2 and >500 μg/m³) that was not within the monitoring range of the National Ambient Air Quality Standard (NAAQS) (GB 3095-2012) to ensure the validity of the PM_2.5 data.

Considering the stability of the predictive model, the collinearity of the independent variables should be diagnosed. In this study, the variance inflation factor (VIF) and tolerance value (TV) were selected to diagnose the collinearity of the selected variables. The VIF and TV of all independent variables participating in the model satisfied VIF < 10 and TV > 0.1 for each year (Table 2), indicating that there was no collinearity problem among the independent variables and could be considered for model fitting.

A two-stage statistical regression model consisting of LME model and GWR model was used to simulate the spatiotemporal variation of the PM_2.5–AOD relationship. The LME model in the first stage was applied to correct the time-varying relationship of PM_2.5–AOD. The quadratic term of AOD (AOD²) and the interaction between PBLH and AOD (PBLH × AOD) were added to the model to explain the nonlinear relationship between AOD and PM_2.5. The specific structure of the model is as follows: PM 2.5 s t = ( β 0 + θ 0 ) + ( β 1 + θ 1 ) AOD s t + ( β 2 + θ 2 ) AOD s t 2 + ( β 3 + θ 3 ) PBLH s t + ( β 4 + θ 4 ) WS s t + ( β 5 + θ 5 ) TEMP s t + ( β 6 + θ 6 ) RH s t + ( β 7 + θ 7 ) PRS s t + ( β 8 + θ 8 ) PRCP s t + β 9 × PBLH s t × AO D st + β 10 × F C s + β 11 × U C s + ε s t ( θ 1   ∼ 8 ) ∼ N [ ( 0,0,0 , ψ ) ] , ε s t ∼ N ( 0 , σ 2 ) (1) where PM_2.5st is the [PM_2.5] at station s on day t; AOD _st is the AOD of the grid cell in which the station s is positioned on day t; AOD² st is the quadratic term for AOD at station s on day t; PBLHst, WSst, TEMPst, RHst, PRSst, and PRCPst are the planetary boundary layer height, wind speed at 10-m height, temperature at 2-m height, relative humidity, surface pressure, and precipitation at station s on day t, respectively; PBLH _st × AOD _st is the interaction between PBLH and AOD at station s on day t; FC _s and UC _s are the FC value and UC value at station s, respectively; β ₀ and θ ₀ are the fixed and random intercepts, respectively; β ₁ and β ₂ are the fixed slopes of square polynomials for AOD; β ₃, β ₄, β ₅, β ₆, β ₇, β ₈, β ₁₀, and β ₁₁ are the fixed slopes of PBLH, WS, TEMP, RH, PRS, PRCP, FC, and UC; β ₉ is the fixed slope of the interaction between PBLH and AOD; θ ₁ and θ ₂ are the daily random slopes of square polynomials for AOD; and θ ₃–θ ₈ are the daily random slopes of each meteorological variables, respectively.

The GWR model of the second stage was used to correct the spatial heterogeneity between PM_2.5 and AOD. The specific method was to model the residuals of the LME model. This GWR model was fitted once a day to account for temporal variability. In addition, the model using adaptive bandwidth selection methods calculated by minimizing the corrected Akaike Information Criterion (AIC) value. The specific expression is as follows: PM 2.5 _ resi st = β 0 ( u s , v s ) + β 1 ( u s , v s ) AOD s t + ε s t   ε s t ∼ N ( 0 , σ 2 )   (2) where PM_2.5_resi _st is the residual value from the LME model at station s in day t; AOD _st is the AOD value at station s on day t; (u _s , v _s ) is the spatial coordinates of the monitoring station s; and β ₀ (u _s ,v _s ) and β ₁ (u _s ,v _s ) represent the regression intercept and regression slope at station s, respectively.

For model verification, a 10-fold cross-validation (CV) method was conducted to detect the degree of overfitting of the model. The entire model-fitting dataset was randomly split into 10 subsets, with each subset containing approximately 10% of the dataset. In each CV time, we selected one subset as the testing sample and used the remaining nine subsets to fit the model for prediction on the testing sample. This process was repeated 10 times to ensure that all the subsets were predicted. We fitted a linear regression was performed between the measured and predicted [PM_2.5], and the fitted R ² , slope, RMSPE, and relative prediction error (RPE) were evaluated the performance of the model. They represented by Eqs. 3 and 4, respectively. R M S P E = ∑ i = 1 n ( y m o d , i − y o b s , i ) 2 n   (3) R P E = R M S P E × 100 % y ¯   (4) where y _mod,i is the estimated PM_2.5 at site i; y _obs,i is the observed PM_2.5 at site i; n is the total number of data samples; and y ¯   is the average of the observed PM_2.5.

The estimation process of the daily [PM_2.5] by the LME + GWR model is shown in Supplementary Figure S1.

As shown in Table 3, the daily minimum and maximum [PM_2.5] in the BTH were ranged from 2 to 3 μg/m³ and from 371 to 499 μg/m³, respectively, which indicated that the pollution degree of different areas in the BTH had considerable differences. The annual average [PM_2.5] during the investigated period from 2013 to 2020 in the BTH were 91.27, 85.93, 72.89, 68.22, 61.02, 53.53, 45.75, and 40.97 μg/m³, respectively, indicating that the PM_2.5 pollution has been on a downward trend in the past 8 years. However, it still exceeded the limit (35 μg/m³) of the national secondary standard for ambient air quality (GB3095-2012). The average annual AOD ranged from 0.37 to 0.69 during the same period. The great difference between the mean FC and UC reflected that most of the monitoring sites in the study area were located inside or around the city. In addition, the ranges of the meteorological variables from 2013 to 2020 are also shown in Table 3.

During the period from 2013 to 2020, the monthly [PM_2.5] monitored in the BTH demonstrated that the median monthly [PM_2.5] presented a U-shaped oscillation for each year (Figure 3). Overall, the [PM_2.5] displayed significant monthly differences, following the change pattern of “high in winter, low in summer, falling in spring and rising in autumn”. In detail, [PM_2.5] displayed a downward trend from January to May, a general stability from June to September, and an upward trend from October to December. The reason for the highest monthly [PM_2.5] in December and January was the combined effect of coal-fired heating in winter and unfavorable meteorological conditions in the BTH, such as low air humidity and weak wind speed.

The comparison of LME + GWR model fitting (Figure 4A) and 10-fold CV results (Figure 4B) from 2013 to 2020 indicated that the model displayed excellent performance in capturing daily [PM_2.5]. For model fitting, the data distribution was concentrated toward the regression line. The R ² ranged from 0.89 to 0.97, indicating that the two-stage model could effectively explain 89%–97% of the ground-level [PM_2.5] variation. The slope ranged from 0.89 to 1.04, indicating that only a small prediction bias remained in the model. In addition, the fitting results also displayed that the RMSPE and RPE were 6.85–24.60 μg/m³ and 16.67%–26.94%, respectively. Compared with model fitting, the 10-fold CV results showed that the CV-R ² , CV-RMSPE, and CV-RPE ranged from 0.85 to 0.95, 7.87 μg/m³ to 29.90 μg/m³, and 19.19–32.72%, respectively. The CV-R ² decreased, and CV-RMSPE and CV-RPE increased, indicating that the model had a slight overfitting. In addition, Figure 4 shows that, when the measured [PM_2.5] exceeds 400 μg/m³, the model had a slight “high value underestimation” phenomenon.

During the study period, the model performed best in 2020. Under the regression of the observed and predicted [PM_2.5] in the study area, CV- R ² was the highest at 0.95, and the CV-RMSPE and CV-RPE were the lowest at 7.87 μg/m³ and 19.19%, respectively. This is mainly attributable to the government’s series of “Air Pollution Prevention and Control Action Plans” (APPCAP), and the [PM_2.5] has been declining year by year. About 93.72% of [PM_2.5] data samples in 2020 were less than 100 μg/m³. In contrast, the model performed the worst in 2013, with the lowest CV-R ² and the largest forecast uncertainty. The main reason was that more than 32.38% of the data samples were more than 100 μg/m³, and relatively discrete data samples increased the difficulty of model fitting. Overall, the LME + GWR model that we have established was robust. Using the LME + GWR model combined with the 1-km MAIAC AOD product could excellently predict the daily near-surface [PM_2.5] with CV-R ² > 0.84, CV-RMSPE < 30 μg/m³, and CV-RPE < 33% in the BTH.

During the period from 2013 to 2020, [PM_2.5] in the BTH showed obvious spatiotemporal variations on different scales of annual, seasonal, and monthly. Overall [PM_2.5] in the BTH revealed a downward trend during the investigated period. Exploring the reasons for the decrease in the [PM_2.5] was inseparable from the national policy control, such as coal-to-gas and energy-saving transformation (Pan et al., 2021). In detail, the APPCAP implemented between 2013 and 2017 has successfully reduced [PM_2.5] (Yue et al., 2020), and the sharp decline in 2017 was closely related to the termination year of the APPCAP in 2017. Moreover, the [PM_2.5] in December 2015 and December 2016 were significantly higher than that in other years. The reason was that El Niño in 2015 enhanced the winter air pollution in northern China (Chang et al., 2016). The high concentration in the winter of 2016 might be influenced by anthropogenic factors (Ding et al., 2021). In addition, the [PM_2.5] in the winter of 2018–2020 decreased slowly compared with 2017, and the light-pollution areas (such as Langfang and Tangshan) slightly expanded. Furthermore, combined with the contribution of the suspension of work and production during new coronavirus disease (COVID-19) (Xian et al., 2021), the [PM_2.5] in the winter of 2020 dropped to 51.42 μg/m³, which was the lowest [PM_2.5] in winter during the study period.

In the BTH, [PM_2.5] presented the significant seasonal variation characteristics of “high in winter, low in summer, and transition between spring and autumn”, which were consistent with previous studies (Wu et al., 2016; Guo et al., 2021; Lu et al., 2021). The high [PM_2.5] in winter was concentrated in cities such as Shijiazhuang, Xingtai, and Handan. In the study areas, pollutant emissions were mainly due to the coal-fired heating and unfavorable meteorological conditions (Lv et al., 2017). Relevant studies have pointed out that the increase of boundary layer height and higher water vapor content in summer are the main reasons for the low [PM_2.5] (Qu et al., 2016; Ding et al., 2021). Moreover, the elevated [PM_2.5] levels in autumn were likely caused by the large scale straw burning in the rural areas and coal burning for heating in November (Duan et al., 2004; Lv et al., 2017). In addition, the spatiotemporal variation trends on the monthly scale follow the characteristics of seasonal changes, with the most polluted months appearing in December and January.

In previous studies, the CV-R ² value range of the satellite-based ground [PM_2.5] estimation model for the BTH was 0.54–0.95 (Table 4). Among these, the [PM_2.5] estimation model based on MODIS MAIAC AOD (CV-R ² up to 0.82–0.95) has been found to perform better than other [PM_2.5] estimation models (a maximum CV-R ² of 0.83), owing to its superior spatial resolution. Under the same high spatial resolution of AOD, our model showed similar or even better performance than other machine learning models. The performance statistics of the LME + GWR model developed was also comparable with other studies conducted in the United States that used the MODIS MAIAC AOD data (CV-R ² up to 0.62–0.84) (Hu et al., 2014a, 2014b; Chudnovsky et al., 2014; Stowell et al., 2020). For model, the LME model cannot estimate the daily value of PM_2.5 at non-monitoring points, even if there are abundant data available. Models such as TEFR and STLME also have this shortcoming (Wu et al., 2016; Wang W. et al., 2021). In addition, machine learning methods that account for complex nonlinear relationships between different variables by adding hidden nodes and layers exhibited good performance in estimating [PM_2.5] (Ni et al., 2018; Stafoggia et al., 2019; Sun et al., 2019; Zhao et al., 2020; Ding et al., 2021). However, the addition of more hidden nodes and layers consumed a lot of time in estimating PM_2.5 load and produced different results for each training (Wang W. et al., 2021). Therefore, the LME + GWR model has certain advantages in terms of performance and stability. For AOD product, AHI AOD was commonly used to estimate hourly [PM_2.5] due to its high resolution. However, AHI cannot retrieve AOD at nighttime, and the quality is slightly inferior to MODIS AOD (Sun et al., 2019; Wang W. et al., 2021).

The model proposed in this study has many advantages. First, the high spatiotemporal resolution MAIAC Terra/Aqua fusion AOD data were employed in the model and achieved satisfactory performance. Second, the AOD quadratic term (AOD²) and the interaction term of AOD and PBLH (PBLH×AOD) were introduced into the first-stage LME model to describe the nonlinear effect of the model. Third, we adopted the GWR model as the second-stage model to improve the spatial difference of the PM_2.5–AOD. The bisquare kernel bandwidth function and adaptive bandwidth method were selected owing to the difference between the daily sample data. After CV, the degree of overfitting was very small (compared with R ² , CV-R ² was only reduced by 0.01–0.04).

However, the model still has some limitations. One limitation was the mismatch in spatial resolution between MODIS MAIAC AOD (0.01° × 0.01°) and meteorological parameters (0.1° × 0.1° and 0.0625° × 0.0625°). Although the bilinear interpolation method for meteorological factors has proved to have better performance than linear interpolation and nearest neighbor interpolation algorithms (Zhao et al., 2020), more meteorological products with high spatial resolution were still needed. Another limitation was that we only keep three data records in some days to bring into the model. Related studies have pointed out that the overfitting degree of the two-stage model incorporating GWR decreases with the increase in the number of matching data records per day (Hu et al., 2014a; Wu et al., 2016). Therefore, too few observations in some days would lead to the GWR model to overfitting. We will explore the optimal threshold that matches the minimum number of data records later. In addition, some studies have indicated that PM_2.5 monitoring stations mostly located in cities and suburbs, and the PM_2.5 estimation in mountainous and rural areas was relatively poor (Zeng et al., 2020; Ding et al., 2021). Our study area is in the BTH with characteristics of urban industrial conditions. In particular, Hebei Province that has many rural administrative units also has a large number of factories. Provincial monitoring sites with a larger coverage area should be added to future research to increase the regional representation of the sample.