Full-coverage 1 km daily ambient PM_2.5 and O₃ concentrations of China in 2005–2017 based on a multi-variable random forest model

The health risks of fine particulate matter (PM_2.5) and ambient ozone (O₃) have been widely recognized in recent years. An accurate estimate of PM_2.5 and O₃ exposures is important for supporting health risk analysis and environmental policy-making. The aim of our study was to construct random forest models with high-performance and estimate daily average PM_2.5 concentration and O₃ daily maximum of 8 h average concentration (O₃-8 hmax) of China in 2005–2017 at a spatial resolution of 1 km × 1 km. The model variables included meteorological variables, satellite data, chemical transport model output, geographic variables and socioeconomic variables. Random forest model based on 10-fold cross-validation was established, and spatial and temporal validations were performed to evaluate the model performance. According to our sample-based division method, the daily, monthly and yearly estimations of PM_2.5 from test datasets gave average model-fitting R² values of 0.85, 0.88 and 0.90, respectively; these R² values were 0.77, 0.77 and 0.69 for O₃-8 hmax, respectively. The meteorological variables and their lagged values can significantly affect both PM_2.5 and O₃-8 hmax estimations. During 2005–2017, PM_2.5 concentration exhibited an overall downward trend, while ambient O₃ concentration experienced an upward trend. Whilst the spatial patterns of PM_2.5 and O₃-8 hmax barely changed between 2005 and 2017, the temporal trend had spatial characteristics. The dataset is accessible to the public at https://doi.org/10.5281/zenodo.4009308 (Ma et al., 2021a), and the shared dataset of Chinese Environmental Public Health Tracking (CEPHT, 2022) is available at https://cepht.niehs.cn:8282/developSDS3.html.

1 Introduction

Air pollution is becoming a main concern of modern society due to various health risks. According to the latest Global Burden of Disease (GBD) report, air pollution has caused approximately 6.67 million deaths (95 % uncertainty intervals: 5.90–7.49 million) and ranked fourth on the global list of death-related risk factors in 2019 (Health Effects Institute, 2020; Murray et al., 2020). Ambient fine particulate matter (PM_2.5) and ambient ozone (O₃) have been identified and proven to be related to many health outcomes. China is known to be one of the countries with the most serious air pollution in the world. Strict pollution control measures (including the Air Pollution Prevention and Control Action Plan and Three-Year Action Plan To Fight Air Pollution) were enacted by the Chinese government to control and reduce air pollution since 2013. The implementation of these measures has resulted in a markable drop of emissions and PM_2.5 concentration. However, the occasional pollution events, as well as the short development history of air quality monitoring network, have brought many difficulties to accurately capture the temporal and spatial patterns of PM_2.5 and O₃ concentrations. Therefore, it is difficult to develop a complete decision-making basis for handling air pollution. In addition, there are gaps in epidemiological studies linking air pollutants to health outcomes, due to the lack of accurate measurements of PM_2.5 and ambient O₃ concentrations. To this end, an accurate estimate of PM_2.5 and O₃ exposures is essential to support health risk analysis and environmental policy-making.

Suitable model variables and advanced estimation method are important to achieve accurate modeling. Basically, PM_2.5 is jointly affected by both natural conditions and human activities over space and time, e.g., aerosol optical depth (AOD), meteorological conditions, geographic factors and human-related features (Wei et al., 2021), while O₃ is a secondary pollutant, which is produced by a series of complex photochemical reactions on the basis of precursors including nitrogen oxides (NO_x) and volatile organic compounds (VOCs) under the action of high temperature and strong radiation. These complex characteristic puts forward higher requirements on the ability of the modeling method to handle multi-variable and capture the non-linear relationships between variables and air pollutants. Many models have been developed to estimate the spatiotemporal distribution of PM_2.5 and O₃ concentrations in China. Machine-learning approaches (e.g., random forest (RF), extreme gradient boosting and deep belief network models) can mine useful information from a large amount of input data and explore the nonlinear relationship as well as bring a better performance in modeling work (Chen et al., 2018, 2019; Di et al., 2017; Li et al., 2017; Wei et al., 2019; Zhan et al., 2018). However, most of these estimation datasets cannot balance long time series and high spatiotemporal resolution. Besides, there is no long-term estimation dataset for both PM_2.5 and O₃ concentrations with high temporal and spatial resolution for supporting epidemiological research. Therefore, by incorporating multi-source data into random forest models, this study makes an attempt to estimate the high-resolution (1 km × 1 km) ambient PM_2.5 and O₃ concentrations of China in 2005–2017.

2 Method

The model variables of this study include meteorological variables, geographical variables, socioeconomic variables, satellite data and chemical transport model output in 2013–2017. Daily average PM_2.5 and O₃ daily maximum of 8 h average concentration (O₃-8 hmax) monitoring data of 1479 sites in 2013–2017 were obtained (Figs. 1; S1 and S2 in the Supplement). A 1 km × 1 km standard grid is created across the country (35.55^∘ N to 43.12^∘ N, and 112.95^∘ E to 120.35^∘ E) with a total of 9495025 grid cells. The coordinate system of the grid is WGS-84, and the projection of the grid is the Albers conical equal-area projection. We construct high-performance random forest models (temporal resolution: daily; spatial resolution: 1 km × 1 km) and estimate the grid daily average PM_2.5 concentration and O₃-8 hmax concentration of China in 2005–2017.

Figure 1 Station distribution in China and average ground monitoring concentration based on the available data of PM_2.5 (a) and O₃-8 hmax (b) from 2013 to 2017.

2.1 Dataset

The model variables used in this study mainly include Aqua AOD for PM_2.5 modeling, GEOS-Chem chemical transport model output for O₃ modeling and some variables shared by PM_2.5 and O₃: 13 meteorological variables (includes boundary layer height, surface pressure, 2 m dew point temperature, evaporation, albedo, low cloud cover, medium cloud cover, high cloud cover, total precipitation, 10 m U wind component, 10 m V wind component, 2 m surface temperature and surface solar radiation downwards) and their lag 1 and lag 2, geographic and socioeconomic variables, such as the digital elevation model (DEM), normalized difference vegetation index (NDVI), population, gross domestic product (GDP), road network and dummy variables (includes season, month, and province). A more detailed description of the model variables is given in Table S1 in the Supplement. The processing method has been described in detail in our earlier studies (Ma et al., 2021b; Zhao et al., 2019). Briefly, most of the model variables are processed into 1 km × 1 km resolution based on the standard grid using interpolation methods (such as inverse distance weighted and bilinear algorithm) in ArcGIS 10.2 and Python 2.7. AOD is processed by ENVI 5.3+IDL and extracted into standard grid using ArcPy, and then the inverse-distance-weighted interpolation is carried out to obtain the 1 km × 1 km resolution data. For the long-term variables, the corresponding monthly and annual level value is assigned to each day. Subsequent modeling work was carried out based on the dataset covering monitoring data and all variables.

2.2 Random forest model

Random forest is an ensemble machine learning method consisting of many individual decision trees growing from bagged data, and its prediction is the result of the voting process (Breiman, 2001). The RF algorithm primarily integrates learning principles, trains several individual learners and finally forms a strong learner through a certain combination strategy; through multiple rounds of training, multiple prediction results are obtained, and the final results are obtained after average aggregation.

The random forest models are established using the 10-fold cross-validation method. First, this method randomly divides the modeling dataset into 10 parts, then nine of them are used for modeling, and the remaining one is used for estimation and is compared with observations. The verification is repeated until every part is predicted. In this way, the modeling and verification of estimation are repeated 10 times in total, and the average values of the 10 runs are taken as the final result, i.e., the CV-R². The formulae of the models are as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\mathrm{PM}}_{\mathrm{2.5}i,j}& =f({\mathrm{METE}}_{i,j},\mathrm{lag}\mathrm{1}{\mathrm{METE}}_{ij},\mathrm{lag}\mathrm{2}{\mathrm{METE}}_{i,j},\\ & {\mathrm{AOD}}_{i,j},{\mathrm{LD}}_j,{\mathrm{ROAD}}_j,{\mathrm{NDVI}}_j,{\mathrm{ELE}}_j,{\mathrm{GDP}}_j,\\ & {\mathrm{POP}}_j,{\mathrm{SEASON}}_i,{\mathrm{MON}}_i,{\mathrm{PRO}}_j),\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{\mathrm{O}}_{\mathrm{3}}\text{-}\mathrm{8}{\mathrm{hmax}}_{i,j}& =f({\mathrm{METE}}_{i,j},\mathrm{lag}\mathrm{1}{\mathrm{METE}}_{ij},\\ & \mathrm{lag}\mathrm{2}{\mathrm{METE}}_{i,j},{\mathrm{GEOS}}_{i,j},{\mathrm{LD}}_j,{\mathrm{ROAD}}_j,\\ & {\mathrm{NDVI}}_j,{\mathrm{ELE}}_j,{\mathrm{GDP}}_j,{\mathrm{POP}}_j,\\ & {\mathrm{SEASON}}_i,{\mathrm{MON}}_i,{\mathrm{PRO}}_j),\end{array}\end{array}$$

where PM_2.5 i,j and O₃-8 hmax_i,j are the PM_2.5 and O₃-8 hmax concentrations on day i in grid cell j; METE_i,j represents 13 meteorological variables on day i in grid cell j, and lag 1 METE_i,j and lag 2 METE_i,j represent the corresponding 1 d lag and 2 d lag values, respectively; GEOS_i,j and AOD_i,j are the GEOS-Chem model output and AOD value on day i in grid cell j; LD_j, ROAD_j, NDVI_j, ELE_j, GDP_j and POP_j are the land use coverage, length of a variety of roads, NDVI, elevation, GDP and population in grid cell j, respectively; and SEASON_i, MON_i and PRO_j are the season and month of day i and the province of grid cell j, respectively.

In general, the random forest parameters that need to be adjusted include n_estimators (number of decision trees) and the max_depth (maximum depth of the trees). Unlike the previous methods of manually adjusting parameters, the parameters of random forest were optimized using GridSearchCV, which can realize a cross-validated grid search over a parameter grid. After GridSearchCV, we set max_depth as 36 and n_estimators as 200 for PM_2.5 modeling. For O₃-8 hmax modeling, we set max_depth as 54 and n_estimators as 200.

2.3 Validation method

To comprehensively verify the model performance, we construct the main models using sample-based division method. Models using spatial-based and temporal-based division method are further constructed to test the model performance on a spatial and temporal scale.

The dataset was randomly divided into the training set (90 % of the records) and test set (10 % of the records) by using the sample-based division method. We construct the main model using the training set with a 10-fold cross-validation. Since the data in the test set are not used in the main model, true model performance can be verified. The coefficient of determination (R²) of the main model on the test set (test-R²) and the verification indicators of model uncertainty, the root mean square error (RMSE) and mean absolute error (MAE) are calculated for the PM_2.5 and O₃-8 hmax model, respectively. The monthly and yearly test-R² values are also calculated.

For the spatial verification, 90 % of the monitoring stations are randomly selected. The monitoring data of these stations are used as the training set, and the monitoring data of remaining stations are used as the testing set. For the temporal verification, 90 % of dates in 2013–2017 are randomly selected, and the data in theses dates are used as the training set; the data in the remaining 10 % of the dates in 2013–2017 are used as the testing set. After that, the test-R², RMSE and MAE are calculated.

2.4 Estimation of daily PM_2.5 and ambient O₃ of China from 2005 to 2017

Based on the final models of PM_2.5 and O₃-8 hmax, we estimate the gridded daily average PM_2.5 concentration and O₃-8 hmax concentration of China in 2005–2017. The spatial pattern and temporal trend of PM_2.5 and O₃-8 hmax concentrations are analyzed and compared with other modeling products.

The modeling and estimations are performed in Python 2.7.13 using the scikit-learn 0.20.3 and GridSearchCV packages. The workflow of this study is displayed in Fig. 2.

Figure 2 The workflow of modeling process in the study.

3 Results and discussion

A total of 981 744 monitoring data records were used in the final model-fitting dataset. The mean ± standard deviation of PM_2.5 and ambient O₃ concentrations in 2013–2017 were 59.60 ± 45.85 and 86.72 ± 47.73 µg m⁻³, respectively. The results of descriptive analysis for variables included in PM_2.5 and O₃-8 hmax model are shown in Table S2.

3.1 Model fitting and validation

The cross-validation results indicate that the estimated PM_2.5 and O₃-8 hmax concentrations matched reasonably with the observed PM_2.5 and O₃-8 hmax concentrations, with high fitted test-R² values. According to our sample-based division method, the test-R² values of the estimated daily, monthly and yearly PM_2.5 concentrations were 0.85, 0.88 and 0.90, respectively (Fig. 3). Likewise, the test-R² values of the estimated daily, monthly and yearly O₃-8 hmax concentrations were 0.77, 0.77 and 0.69, respectively (Fig. 4). The RMSE and MAE for PM_2.5 in daily level were 17.72 and 9.37 µg m⁻³; for O₃-8 hmax, the values were 23.10 and 15.43 µg m⁻³. The model performance is comparable to previous studies (Di et al., 2017; Li and Cheng, 2021; Liu et al., 2020; Wei et al., 2021, 2020, 2019). At the provincial/city level, the model performance of PM_2.5 estimations of Shanghai, Beijing, Hubei, Hebei and Sichuan ranked in the top five with relatively high test-R² (≥ 0.90), while those of Tibet, Qinghai, Gansu, Anhui and Yunnan were less accurate with relatively low test-R² values (< 0.70). The model performance of O₃-8 hmax estimations of Beijing, Chongqing, Shanghai, Tianjin and Henan ranked in the top five with relatively high test-R² values (≥ 0.83), while those of Gansu, Anhui, Heilongjiang, Guizhou and Tibet were poorer with relatively low test-R² values (< 0.62) (Table S3).

Figure 3 The density plot of PM_2.5 model. From left to right are the different temporal scales: daily, monthly and yearly; from top to bottom are the different validation methods: sample-based, spatial-based and temporal-based.

Figure 4 The density plot of O₃-8 hmax model. From left to right are the different temporal scales: daily, monthly and yearly; from top to bottom are the different validation methods: sample-based, spatial-based and temporal-based.

The spatial and temporal test-R² of our models explained the uncertainty to some content (Figs. 3 and 4). The spatial test-R² values for daily, monthly and yearly PM_2.5 estimation were 0.83, 0.87 and 0.85, respectively, while those of daily, monthly and yearly O₃-8 hmax estimations were 0.74, 0.77 and 0.68, respectively. The relatively high spatial test-R² demonstrates the reasonable performance of our models in areas without monitoring stations. The temporal test-R² values of daily, monthly and yearly PM_2.5 estimations were 0.49, 0.65 and 0.76, respectively, while those of daily, monthly and yearly O₃-8 hmax estimations were 0.58, 0.63 and 0.56, respectively. These results indicate the uncertainty of our models when modeling data in the historical period, although the performance is among the best compared with previous studies. The simulation accuracy is a universal issue in the present studies of air pollutant concentrations in the historical period without monitoring data. Further efforts are need to improve the model performance of historical estimations.

3.2 Feature importance

The feature importance of the variables in our random forest models is presented in Tables S4-1 and S4-2. Similar to previous studies (Chen et al., 2018; Zhan et al., 2018), the meteorological factors and their lagged values can significantly affect both PM_2.5 and O₃-8 hmax modeling. Moreover, the specific features for PM_2.5 and O₃, AOD and GEOS-Chem output, also demonstrated high importance in modeling work.

For PM_2.5 modeling work (Table S4-1), the meteorological variables (boundary layer height, evaporation, 2 m dew point temperature) and their lagged effect were among the top 10 important factors, totaling 33.6 % in modeling work. The lagged effects greatly contributed to PM_2.5 modeling. For example, the lag 1 boundary layer height ranked first (17.2 %) in our study, which is similar to previous studies (Zhao et al., 2019). The interpolated AOD (5.6 %), DEM (4.9 %) and season (3.7 %) also demonstrated high importance, which showed crucial effects of satellite data, terrain distribution characteristics in the study area and study period on PM_2.5 modeling. The relative contributions of land use, NDVI, population density, road length and GDP are negligible (the importance scores less than 1 %). Unlike DEM, these factors are subjected to the influence of socioeconomic status in study area. In the future study, the integration of these factors with a higher temporal resolution might change its contribution to the estimation.

The feature importance of ambient O₃ is consistent with its formation and dissipation mechanism: surface solar radiation downwards and its lagged effect according for 39.2 % in modeling work (Table S4-2). Other meteorological factors (2 m temperature, boundary layer height, 10 m V wind component and low cloud cover) according for totaling 9.54 % importance scores. Our analysis also suggests the high importance of GEOS-Chem model (7.2 %), altitude (1.9 %), and dummy factors including year (2.2 %) and province (1.6 %) in O₃ modeling. By contrast, the relative contributions of land use, NDVI and road length are negligible (the importance scores less than 1 %). The high importance rank of population and GDP might be attributed to the relatively high sensitivity of O₃ to anthropogenic emission sources (compared to PM_2.5).

3.3 The spatial characteristics and temporal trend of PM_2.5 and ambient O₃ of China from 2005 to 2017

During 2005–2017, PM_2.5 showed an overall downward trend, while ambient O₃ showed an upward trend in recent years (Figs. 5, S3–S6). Relative to 2005, PM_2.5 concentration has increased by 2.60 µg m⁻³ in 2013. Nevertheless, after the implementation of the Air Pollution Prevention and Control Action Plan, a strict pollution control measure, PM_2.5 concentration has declined by 11.041 µg m⁻³ in 2017 (relative to 2013). This has resulted in a downward trend of PM_2.5 concentration in 2005–2017: PM_2.5 concentration in 2017 has decreased by 8.44 µg m⁻³ relative to 2005 (Figs. 5 and S3). In key pollution areas, with the implementation of various air pollution prevention and control policies, PM_2.5 levels in the Beijing–Tianjin–Hebei region have dropped the most, but the overall concentration levels are still higher than those in the Yangtze River Delta and Pearl River Delta (Fig. S4). For O₃-8 hmax, the upward trend barely changed. Relative to 2005, O₃-8 hmax concentrations in 2013 and 2017 have increased by 0.39 and 7.83 µg m⁻³, respectively. The upward trend during 2005–2017 was mostly due to the significant changes between 2013 and 2017: relative to 2013, the O₃-8 hmax concentration has increased by 7.44 µg m⁻³ in 2017 (Figs. 5 and S5). The Beijing–Tianjin–Hebei region has shown an obvious upward trend since 2013, while the Pearl River Delta region change trend is not obvious (Fig. S6). During the strict pollution control period, VOC emissions were not effectively controlled, which could be one of the main reasons. Therefore, integrated management of VOCs and NO_x in key industries and areas is important.

Figure 5 The temporal trend of PM_2.5 and O₃-8 hmax concentration in China from 2005–2017. The black dots represent the monthly average PM_2.5 and O₃-8 hmax concentration from 2005 to 2017, the blue color band represents the range of the monthly average PM_2.5 and O₃-8 hmax concentration plus or minus the RMSE value from 2013–2017 (period with monitoring data), and the green color band represents the range of the monthly average PM_2.5 and O₃-8 hmax concentration plus or minus the MAE value from 2013–2017 years.

The seasonal distributions of PM_2.5 and O₃-8 hmax concentrations were obvious during 2005–2017 (Figs. S7 and S8). The lowest seasonal PM_2.5 concentration occurred in summer, with an average concentration of 33.6 ± 11.39 µg m⁻³, and the highest seasonal PM_2.5 concentration occurred in winter, with an average concentration of 57.4 ± 21.76 µg m⁻³. In winter, temperature inversion occurs frequently, and the thickness of the mixed layer is low, which is not conducive to the diffusion of pollutants, which leads to the accumulation of PM_2.5 near the ground (Sun et al., 2014). In contrast, the lowest seasonal O₃-8 hmax concentration was in winter, with an average concentration of 72.65 ± 6.28 µg m⁻³; the highest seasonal O₃-8 hmax concentration was in summer, with an average concentration of 97.44 ± 13.58 µg m⁻³. Temperatures and solar radiation conditions in summer increase the incidence of severe O₃ pollution events, which is consistent with its formation and dissipation mechanism.

The PM_2.5 concentrations in the Beijing–Tianjin–Hebei, Chengdu–Chongqing and Xinjiang regions are higher than other regions, followed by the central China. The PM_2.5 concentrations in the southwestern regions (Yunnan and Tibet) and western part of Sichuan Province are the lowest, followed by the inner-north regions and the south and southeastern regions (Figs. 6, S3 and S4; Table S5). The O₃-8 hmax concentrations in the Bohai Rim, Yangtze River Delta, Pearl River Delta, and other economically developed regions like southern Xinjiang, Inner Mongolia, and northeastern Gansu are relatively high (Figs. 6, S5 and 6; Table S5). This spatial pattern barely changed during 2005–2017 (Figs. S3 and S5), but the temporal trend showed spatial characteristics (Figs. 6; S4 and S6). For PM_2.5 concentration, the key pollution areas were severely polluted during 2005–2013. The air pollution control measures of these regions were strict during 2013–2017; thus, the decline was obvious, especially for the Beijing–Tianjin–Hebei region. For O₃-8 hmax concentration, the growth rate was not obvious (except for the eastern part of Hubei Province) during 2005–2013. However, after 2013, there was a clear upward trend across the country, especially in the northern China.

Figure 6 Estimated annual mean and difference of PM_2.5 and O₃-8 hmax concentration in China during 2005 to 2017. The first row shows maps of PM_2.5-related indicators, and the second row shows maps of O₃-8 hmax related indicators. From left to right are the average concentrations during 2005–2017, the difference between 2017 and 2005, the difference between 2013 and 2005, and the difference between 2017 and 2013.

3.4 Evaluation of the PM_2.5 and O₃ concentration products with comparison with other products

Our estimation datasets include the PM_2.5 and O₃-8 hmax concentration data of China in 2005–2017 with a spatial resolution of 1 km × 1 km resolution. With high spatial and temporal resolutions, our validation results are comparable with other modeling work (see Table S6). Considering the future application in epidemiological research, our estimation datasets would be useful: for acute effects studies, the high spatial resolution would effectively reduce exposure errors; for chronic effects studies, long-term exposure data are essential for the development of cohort studies.

Nevertheless, our estimation datasets also contain some limitations. First, we did not use emission data in our model limited by coarse resolution. However, the newly published high-resolution emission inventory of China (http://meicmodel.org/, last access: 18 February 2022) may be utilized in future estimation studies to improve accuracy. Second, our modeling still has spatial and temporal uncertainties. In areas where monitoring sites are sparsely distributed, such as western China, it may be difficult to accurately capture the association between air pollution concentrations and variables. The model validation of the historical period is also limited. Third, the interpolation process of model features inevitably introduces systematic errors. Therefore, more high-quality and high-resolution basic data would be needed in the future.

4 Data availability

The estimated PM_2.5 and O₃ data are freely accessible at https://doi.org/10.5281/zenodo.4009308 (Ma et al., 2021a), and the shared dataset of Chinese Environmental Public Health Tracking (CEPHT, 2022) is available at https://cepht.niehs.cn:8282/developSDS3.html.

5 Conclusions

We constructed random forest models for simulating of daily average PM_2.5 and O₃-8 hmax concentrations of China during 2005–2017, with a feature list and comparable model performance. The estimation dataset would be useful for supporting both long-term and short-term epidemiological studies. The model can be further used for simulating daily concentrations of longer time period. The key findings are summarized as follows. First, the RF model proved its superiority in our study and can be further used in the future estimation of air pollutant concentration. Second, meteorological data are the most sensitive to PM_2.5 and O₃ modeling. For PM_2.5 modeling work, boundary layer height, evaporation, 2 m dew point temperature and its lagged effects showed the highest sensitivity. For O₃ modeling work, surface solar radiation downwards and its lagged effect were the most sensitive. Third, PM_2.5 concentration has trended downward in China, and the key polluted areas during 2005–2013 were effectively controlled during 2013–2017. O₃ concentration has trended upward in China, especially in the northern China during 2013–2017.