2.1 Multisource input data

Observational PM_2.5 data. The MEE began laying out a PM_2.5 monitoring network in January 2013, expanding the scope from key regions, including the North China Plain (NCP), the Yangtze River Delta (YRD), the Pearl River Delta (PRD), the Sichuan Basin (SB), municipalities directly under the Central Government and provincial capitals to 113 key and model cities for environmental protection and, eventually, to all cities above prefecture level. This expanded the number of observation sites from the initial 520 to over 1600. Since then, PM_2.5 mass concentrations have been recorded continuously using the β-absorption methods or a micro-oscillating balance following a standard protocol (Huang et al., 2021). Hourly PM_2.5 data of all sites between 2013–2020 are collected from the China National Environmental Monitoring Center (CNEMC, http://www.cnemc.cn, last access: 10 July 2022). To produce high-quality PM_2.5 data, a series of quality controls were conducted, including integrity checking, duplicate rejection, and outlier handling. All sites with a proportion of valid PM_2.5 records exceeding 60 % were considered. For each site, identical data for three consecutive hours were excluded first, and PM_2.5 values over three standard deviations from 24 h and 3 d moving average were regarded as outliers and discarded then. Eventually, PM_2.5 data from 1485 sites remained for model development and application. In addition, pre-2013 PM_2.5 measurements in US embassies in Beijing and Shanghai are used for independent validation evaluations (http://www.stateair.net/web/historical, last access: 10 July 2022).

Visibility and conventional meteorological data. The CMA established a national meteorological observation network in the 1950s with the station numbers exceeding 2000 at the beginning and stabilizing at ∼2450 afterward. The observation network can continuously record meteorological data on visibility and conventional meteorological variables, including temperature, pressure, wind, and RH. In recent years, meteorological observations, including 6-hourly records between 1960–2020 and gradually increasing hourly records after 2013, have been collected from the National Meteorological Information Center (NMIC). Due to the inconsistency of visibility data in terms of observation methods, we conducted a series of data conversions to ensure continuous and consistent data. Visibility data recorded on a scale ranging from 0 to 9 between 1960–1979 were converted to numerical data based on probability density distributions. Specifically, the probability density distribution of the visibility for each of the 10 years before and after 1980 was calculated at first. The numerical visibility from 1980 to 1989 was graded into classes, with the median value of each class being the corresponding value for each station, and finally, the class observations were converted into numerical observations. From September 2013 to 2016, visibility measurements gradually shifted from 6-hourly manual observations to 1-hourly automatic observations site-by-site. In keeping with manual measurements, the automatic records, which are slighter lower than manual measurements, were calibrated by dividing 0.75 following the guideline from the CMA (CMA, 2014).

Emission inventories and elevation. Historical anthropogenic emissions from 1960–2012 are taken from Peking global emission inventories, developed using a bottom-up approach with spatial resolution at 0.1^∘ × 0.1^∘ and temporal resolution at 1-month intervals (http://inventory.pku.edu.cn, last access: 10 July 2022) (H. Chen et al., 2016; Huang et al., 2014, 2015; Wang et al., 2014). Current anthropogenic emissions during 2013–2020 are from the Multi-resolution Emission Inventory in China (MEIC, http://meicmodel.org, last access: 10 July 2022) (Zhang et al., 2009; Zheng et al., 2018, 2021). Six emission variables from these two inventories are used as inputs for model development, including PM_2.5, NO_x, SO₂, NH₃, BC, OC, and CO. Thirty meter elevation data are collected from the global digital elevation model (GDEM) version 2 (https://earthexplorer.usgs.gov, last access: 10 July 2022). Both emission and elevation data are interpolated from grids to sites to match existing PM_2.5 sites.

Auxiliary data. Monthly normalized difference vegetation index (NDVI) products are downloaded from the Level-1 and Atmosphere Archive & Distribution System Distributed Active Archive Center (LADDS DAAC, https://ladsweb.modaps.eosdis.nasa.gov, last access: 10 July 2022). Land cover classification data are taken from the National Geographic Information Resources Catalogue Service System (https://www.webmap.cn/mapDataAction.do?method=globalLandCover, last access: 10 July 2022). Population data are taken from the Gridded Population of the World version 4 (GPWv4, https://sedac.ciesin.columbia.edu/data/collection/gpw-v4, last access: 10 July 2022) and are calibrated based on the total population in China City Yearbooks. NDVI, land cover, and population data are also interpolated according to PM_2.5 sites and trained as inputs for model development. However, during the model training process, we found that these data had little or no improvement in the hindcast capability of the model, and the time span of these data is insufficient for long-term historical retrieval. Hence, these auxiliary data are not used in model building.

Figure 1A conceptual scheme for constructing long-term historical site-based and gridded PM_2.5 records based on long-term visibility, conventional meteorological observations, emissions, and elevation.

2.2 Spatiotemporal feature extraction

For each PM_2.5 site, we extract five variables as temporal inputs, including the year, month, day, hour, and day of year. The longitude and latitude variables are taken out as location inputs (Fig. 1b). Visibility, RH, and temperature from the nearest meteorological station of each PM_2.5 are used as basic meteorological inputs. The distance between these two sites was also added as a feature. In addition to the influence of the nearest meteorological station, PM_2.5 concentrations at a site are also affected by surrounding conditions. For example, transport of pollution due to air movement is the main cause of heavy pollution episodes in the early stage (Zhong et al., 2017, 2018). Hence, we need to consider spatial effects from surrounding meteorological stations. Our previous study developed a novel feature engineering approach, which incorporated surrounding impacts by extracting spatial features (Zhong et al., 2021). Specifically, the remaining 19 nearest stations were matched for each PM_2.5 site, except the nearest meteorological station. Five variables, including longitude, latitude, temperature, visibility, and RH, were selected from the 19 stations. Then, we calculated the maximum, the minimum, the average value, the skewness value, and the standard deviation for each of the five variables. These produced features, which take advantage of surrounding conditions, are also considered as inputs. After spatiotemporal feature extraction, a total of 71 features were used as inputs for model training. To reduce computation and training time with guaranteed accuracy, the top 40 features in order of importance during small-sample testing processes are used for the following model training and hindcasting. These features included visibility, temporal features, spatial features, emission features, and elevation.

2.3 Gridded input construction

In the previous construction of input features for PM_2.5 sites, we used location information, time information, meteorological information from 20 surrounding meteorological stations, emission information, and elevation. If we assume that each cell in the grid cells is a virtual PM_2.5 site, then it is possible to generate input features for each grid point. After the model is trained based on input features and PM_2.5 concentrations at real PM_2.5 sites, we can feed the gridded input data into the model in turn and, consequently, construct a gridded PM_2.5 network. Therefore, we define a grid area at 0.25 ^∘ × 0.25^∘ with longitude from 70 to 150^∘ E and latitude from 10 to 60^∘ N and select the grid points covering mainland China. For each grid point, we performed spatiotemporal feature extraction and generated the same 71 input features as those of real PM_2.5 sites.

2.4 Model description

LightGBM is one of the state-of-the-art gradient boosting frameworks with better accuracy, lower memory usage, faster training speed, and the capability of handling large-scale data (Ke et al., 2017). Our previous research used this machine learning model to predict PM_2.5 mass concentrations, which shows an unprecedented predictive capacity on hourly, daily, monthly, and annual timescales (Zhong et al., 2021). This study will continue to use this algorithm and previously tuned hyperparameters for model development (Zhong et al., 2021). For hindcasting historical PM_2.5 datasets prior to 2013, a LightGBM model is trained and validated based on PM_2.5 observations and feature inputs from 2013 to 2020. The hindcast capability is validated using cross-validation methods, which are standard methods for parameter tuning and model validation in machine learning. The training dataset is divided into several parts, one of them is used as test data and the remaining parts are used as training data in turn. Each result yields a corresponding evaluation value, which is then averaged to provide an estimate of the model's accuracy. This estimation is quantified by two metrics: the coefficient of determination (R²) and root mean square error (RMSE). The hindcast capability is also validated using PM_2.5 observations from the US embassies in Beijing and Shanghai, which have been observing PM_2.5 data since as early as 2008. After model training and validation, historical 6-hourly input data are inputted into this model to reconstruct a site-based PM_2.5 dataset at 6 h intervals from 1960 to 2020, and gridded input data are inputted into the model to reconstruct a 0.25^∘ × 0.25^∘, 6-hourly, gridded PM_2.5 dataset. The daily, monthly, yearly, and decadal average PM_2.5 concentrations for each site and each grid are also calculated based on the two datasets. Monthly average values were obtained with daily values no less than 20 d; otherwise, they will be missing. Year-average values were calculated with 12 valid month values and decadal-average values were calculated with 10 valid year-average values. The flowchart for reconstructing PM_2.5 datasets is shown in Fig. 1.

3.1 Evaluation of model hindcast performance

The hindcast performance of our model is evaluated using two CV methods, including 10-fold CV and by-year CV. The 10-fold CV partitions the original training datasets into 10 subsamples, one of which is retained as the validation data in turn for testing the model and the remaining 9 subsamples are used as training data. This method is the most common CV that can be compared with results in other studies. However, 10-fold CV often overestimates the model's ability to hindcast continuous historical data. Therefore, we also use by-year CV, during which one year of data is selected sequentially for testing and the remaining data are used for model training. This method is specifically designed to evaluate the hindcast capability of the model.

Table 1Model performance in primary predictors, temporal resolution, and hindcast capability compared with other national PM_2.5 datasets in China.

Download Print Version | Download XLSX

Table 1 compares our dataset and the available datasets in primary predictors, temporal resolution, and CV results (Ma et al., 2016; Fang et al., 2016; Liu et al., 2017; Xiao et al., 2018; Xue et al., 2019; Liang et al., 2020; Huang et al., 2021; Wei et al., 2021a; Van Donkelaar et al., 2021; Geng et al., 2021; Bai et al., 2022). AOD-based datasets are only available from around 2000 at the earliest, with temporal resolutions ranging from daily scale to monthly scale. In contrast, our visibility-based dataset spans 61 years from 1960 to 2020 at 6-hourly intervals, showing a clear advantage in terms of time span and resolution. The R² and RMSE values of our 10-fold CV results are 0.78 and 21.14 µg m⁻³ for 6-hourly estimations, respectively, which indicates our model is quite robust in estimating PM_2.5. Due to a reduction in data amount, the R² and RMSE values further improved to 0.85 and 16.11 µg m⁻³ for daily estimations and 0.92 and 7.90 µg m⁻³ for monthly estimations. This result is comparable to or even better than those of other available datasets whose 10-fold CV R² ranges from 0.61 to 0.80 on a daily scale and from 0.71 to 0.93 on a monthly scale. Our by-year CVs R² and RMSE values are 0.71 and 25.63 µg m⁻³ for 6-hourly estimations, which indicates our model is still robust in hindcast performance. The by-year CV R² values for daily and monthly estimations (0.78 and 0.83) are higher than those in other available datasets (0.41–0.62 and 0.80), which might be partly attributed to spatial feature extraction and the large volume of our training dataset. Zhong et al. (2021) has shown that extracting spatial features can result in a better hindcast performance by fully representing dimensional heterogeneity. Compared to hundreds of thousands to millions of training samples in AOD-based models, the training samples for the visibility-based model are over 100 million. An increase in the order of magnitude for training datasets will yield better results in machine learning.

Figure 2Density scatterplots of observed PM_2.5 and estimated PM_2.5 across China for by-year CV from 2013 to 2020. The time resolution for CV results is hourly and 6-hourly between 2013-2017 and hourly between 2017–2020. Colors are probability distribution densities.

Download

The refined by-year CV results for each year between 2013–2020 are shown in Fig. 2. The by-year CV R² lies between 0.58 and 0.79 with better hindcast performance after 2014. The potential reasons why the R² value in 2013 is slightly lower than those in other years are as follows. First, the PM_2.5 observation network was just established in 2013, during which dehumidification systems, processing procedures, and data quality control methods were incomplete and, therefore, the overall data quality cannot be guaranteed. With the improvement of the observation network after 2014, both the quality and quantity of observations increase significantly. This situation where data quality is relatively low initially but increases over time is also found in O₃ observations. Second, the CMA began to convert some of the manual visibility observations to automatic observations in 2013, during which there were also some irregular procedures in instrument equipment, observation steps, and data quality control. Lastly, although we have corrected the biases between manual and automated observations, some biases may still exist. However, the biases are further reduced as we integrate all manual visibility observations in 2013 into our training dataset.

The model's hindcast capability is further evaluated independently using pre-2013 PM_2.5 observations. For the PM_2.5 data currently available, only the US embassies in Beijing and Shanghai have at least one year's PM_2.5 observations. Therefore, PM_2.5 data from these two sites are applied as an independent evaluation dataset. Figure 3 shows our estimated PM_2.5 are in close agreement with in situ measurements in Beijing and Shanghai, where the overall R² between observations and estimations is 0.74 and 0.79, respectively. For each year between 2008–2012 in Beijing, the R² values fluctuated between 0.70 and 0.81, reflecting a stable and accurate by-year hindcast capability. As shown in Fig. 3c–h, the low values, high values, and temporal variations in PM_2.5 measurements are all well estimated. In particular, PM_2.5 measurements are lacking at the US Embassy in Beijing in early 2008 and around 2009, but our model can provide reasonable and continuous estimations to fill in the gaps. This ability can also be used to fill in missing PM_2.5 observations of MEE from 2013 onwards, building a complete PM_2.5 dataset. Overall, the independent validation results show that historical PM_2.5 data can be well reconstructed by our model.

Figure 3(a) Density scatterplots of observed PM_2.5 and estimated PM_2.5 between 2008–2012 at the US Embassy in Beijing; (b) density scatterplots of observed PM_2.5 and estimated PM_2.5 in 2012 at the US Embassy in Shanghai; (c–g) time series of observed PM_2.5 and estimated PM_2.5 for each year between 2008–2012 at the US Embassy in Beijing; and (h) time series of observed PM_2.5 and estimated PM_2.5 for each year in 2012 at the US Embassy in Shanghai.

Download

The model's ability to make PM_2.5 predictions at locations outside the scope of the training stations is evaluated by spatial CV. For spatial CV, all the monitoring stations are randomly divided into five subsets, and the model is trained using data from four subsets and tested on the data from the remaining subset each time. As shown in Fig. 4, the R² for spatial cross-validation in different groups is between 0.75 and 0.79, reflecting robust predictive power for PM_2.5 concentrations at sites outside the training sites. Our previous study also examined this predictive ability using PM_2.5 data from 23 untouched regional PM_2.5 stations (Zhong et al., 2021).

Figure 4Density scatterplots of observed PM_2.5 and estimated PM_2.5 for each group of spatial CV results.

Download

Figure 5Spatial distribution of annual average PM_2.5 mass concentration at 1485 stations from 1960 to 2020.

3.2 Spatiotemporal variations in the site-based PM_2.5 dataset during 1960–2020

Figure 5 shows the spatiotemporal variations in annual average site-based PM_2.5 between 1960–2020. The trend of PM_2.5 in China experiences three major stages, corresponding to a slow increase under low concentrations between 1960–1978, a continuous accumulation with high concentrations reached between 1979–2013, and a rapid decrease between 2014–2020. During the first stage, though PM_2.5 pollution occurred in parts of the NCP and the Guanzhong Plain (GZP), PM_2.5 concentrations remain low in the vast majority of areas. This is mainly because anthropogenic emissions of PM_2.5 precursors and primary PM_2.5 grow slowly at a low base, resulting in relatively low total emissions in different regions. However, PM_2.5 pollution still occurring in the NCP and GZP, even with relatively low emissions, indicates the low environmental capacity of these two regions. During the second stage, PM_2.5 reached an unprecedentedly high concentration after a continuous increase in nearly all regions in China. The heaviest PM_2.5 pollution occurred in the NCP and the GZP. The SB and the Northeast China Plain (NeCP) are the polluted regions with the next highest PM_2.5 pollution. Even the YRD and the PRD also experienced PM_2.5 pollution during this stage. This worsening of PM_2.5 pollution is closely associated with massive anthropogenic emissions from rapidly increasing living and industrial activities after reform and opening-up policies. From 1979 to 2013, primary PM_2.5, NO_x, SO₂, NH₃, BC, OC, and CO from the Peking emission inventory increased by 98 %, 457 %, 159 %, 117 %, 45 %, −22 %, and 243 %, respectively. Despite a slow reduction in SO₂ after 2006, the total anthropogenic emissions each year still increased and thereby caused high-level PM_2.5 pollution after 2006. The results indicate that air pollutants cannot be emitted without restraint even in regions with high atmospheric capacity. Otherwise, PM_2.5 pollution will inevitably occur. In addition to anthropogenic emissions, sand and dust storms, resulting in high PM_2.5 concentrations in western Xinjiang, worsened PM_2.5 pollution by trans-regional transport from the desert regions. During the last stage, PM_2.5 decreased nationwide with the mass concentrations in nearly all stations approximately or below 35 µg m⁻³ in 2020, even in the NCP and the GZP with limited environmental capacity. The substantial declines in PM_2.5 illustrate the effectiveness of implementing the toughest ever clean air policy in China. The spatiotemporal variations of PM_2.5 between 1960–2020 clearly show the long-term impact of economic development and energy consumption on our air quality and the effectiveness of recent years' unprecedented emission control policies.

The specific turning points in annual PM_2.5 concentrations for different regions were investigated additionally. Figure 6 shows the temporal variations in national average monthly and yearly PM_2.5 mass concentrations and regional average 6-hourly, monthly, and yearly PM_2.5 mass concentrations in “2+26” cities of the NCP, the YRD, the PRD, and the SB. The national average yearly PM_2.5 reached a peak of 67 µg m⁻³ in 2007, declined in 2008, and then remained steady until 2013. A sharp fall followed after 2014, with PM_2.5 concentrations decreasing from 63 µg m⁻³ in 2013 to 34 µg m⁻³ in 2020. The annual PM_2.5 concentrations in the “2+26” cities also experienced similar changes with a peak in 2007 and a reduction in 2008, which might be related to emission reductions for the Beijing Olympics in 2008. For the YRD, the maximum value of PM_2.5 mass concentration occurred in 2013 without a striking peak in 2007. For the PRD, the annual PM_2.5 concentrations increased steadily between 1960–1978, then rose more and more steeply in the following years, with a steep increase in 2003 and 2004 and peaked in 2004. A steady decrease with slight fluctuation occurred from 2005 to 2013, and then a sharp fall followed after 2014. This trend is different from that in the “2+26” cities and the YRD. For the SB, the turning point occurred in 2013, before which the annual PM_2.5 concentrations increased steadily and remained steady.

Figure 6(a) The spatial distribution of average PM_2.5 mass concentrations between 1960–2020 and the (b–f) time series of average PM_2.5 mass concentrations for all sites in China (b), “2+26” cities (c), Yangtze River Delta (d), Pearl River Delta (e) and Sichuan Basin (f), respectively.

3.3 Detailed spatial distributions from gridded PM_2.5 datasets

Figure 7 shows the annual spatial variations in 0.25^∘ × 0.25^∘ gridded PM_2.5 between 1960–2020. Compared to site-based distributions, gridded PM_2.5 can portray the spatiotemporal variations in a clearer and more detailed way. For example, the most widespread and heaviest PM_2.5 pollution in western Xinjiang occurred in 1979. This abnormal pollution corresponds to the historical construction of northern severe dust storms, which recorded the event with the largest affected areas in April 1979 (Zhou and Zhang, 2003). As exposed to nearly the most frequent air stagnation in winter due to terrain and meteorological conditions (Wang et al., 2018), the NCP is the region with PM_2.5 pollution first to appear and last to disappear, except for areas affected by dust storms (Fig. 7). For year-to-year comparisons, it can be clearly seen that PM_2.5 concentrations in the NCP decreased slightly from 2007 to 2008 and from 2012 to 2013, respectively, and decreased significantly in 2014 relative to 2013. The PM_2.5 reduction is insignificant from 2015 to 2016, but striking from 2016 to 2017. In 2020, the nationwide PM_2.5 concentrations are comparable to those in the 1960s and close to the lowest level ever recorded in almost 61 years.

Figure 7Gridded distribution of annual average PM_2.5 mass concentration from 1960 to 2020.

Figure 8 shows inter-decadal spatial variations in gridded PM_2.5 between 1961–2020. PM_2.5 concentrations maintained at low levels in most areas over the first decade and increased to a certain extent in the NCP and western Xinjiang over the second decade. In the following decades, PM_2.5 pollution has worsened significantly in several key regions, including the NCP, the GZP, and the SB. This worsening was maintained until the last decade, during which PM_2.5 pollution mitigates significantly in nearly all populous and polluted regions in eastern China.

Figure 8Gridded distribution of decadal average PM_2.5 mass concentration from 1960 to 2020.

Figure 9Distribution of reconstructed PM_2.5 by different PM_2.5 datasets in 2000, 2005, 2010, 2015, and 2020. From top to down are TAP, GEFPM, LGHAP, CHAP, and our dataset.

The multi-year trend of our gridded PM_2.5 dataset is also compared with those of publicly available datasets, including the TAP data (Geng et al., 2021), the GEFPM data (Van Donkelaar et al., 2021), the LGHAP data (Bai et al., 2022), and the CHAP data (Wei et al., 2021a), which have been interpolated to the same grid resolution. Figure 9 shows the spatial distributions of PM_2.5 from those datasets at 5-year intervals between 2000–2020. One consistent trend across all datasets was that nationwide PM_2.5 mass concentrations experienced an increase following a decrease from 2000 to 2020. However, the turning points are different for different datasets. From 2010 to 2015, PM_2.5 pollution alleviated for TAP, CHAP, and our data, but worsened for GEFPM and LGHAP. For the time (2015 and 2020) with ground observations available, all PM_2.5 data show similar spatial distributions with the most severe pollution in the NCP in 2015 and significant improvement in nationwide air pollution in 2020. For the years (2000, 2005, and 2010) when ground observations were unavailable, significant disparities in pollution levels and regional distribution emerged from different datasets. Specifically, the LGHAP data are significantly lower than other data, while the TAP data are higher than others in nearly all regions except western Xinjiang. In western Xinjiang, PM_2.5 concentrations from the GEFPM data are the highest among all the datasets. Due to a lack of ground PM_2.5 observations before 2000, it is challenging to determine which dataset has the least bias and more reasonable distributions. In the future, applying an ensemble average to multi-datasets might be an effective way to eliminate systematic bias.