LGHAP: the Long-term Gap-free High-resolution Air Pollutant concentration dataset, derived via tensor-flow-based multimodal data fusion

Developing a big data analytics framework for generating the Long-term Gap-free High-resolution Air Pollutant concentration dataset (abbreviated as LGHAP) is of great significance for environmental management and Earth system science analysis. By synergistically integrating multimodal aerosol data acquired from diverse sources via a tensor-flow-based data fusion method, a gap-free aerosol optical depth (AOD) dataset with a daily 1 km resolution covering the period of 2000–2020 in China was generated. Specifically, data gaps in daily AOD imageries from the Moderate Resolution Imaging Spectroradiometer (MODIS) aboard Terra were reconstructed based on a set of AOD data tensors acquired from diverse satellites, numerical analysis, and in situ air quality measurements via integrative efforts of spatial pattern recognition for high-dimensional gridded image analysis and knowledge transfer in statistical data mining. To our knowledge, this is the first long-term gap-free high-resolution AOD dataset in China, from which spatially contiguous PM_2.5 and PM₁₀ concentrations were then estimated using an ensemble learning approach. Ground validation results indicate that the LGHAP AOD data are in good agreement with in situ AOD observations from the Aerosol Robotic Network (AERONET), with an R of 0.91 and RMSE equaling 0.21. Meanwhile, PM_2.5 and PM₁₀ estimations also agreed well with ground measurements, with R values of 0.95 and 0.94 and RMSEs of 12.03 and 19.56 µg m⁻³, respectively. The LGHAP provides a suite of long-term gap-free gridded maps with a high resolution to better examine aerosol changes in China over the past 2 decades, from which three major variation periods of haze pollution in China were revealed. Additionally, the proportion of the population exposed to unhealthy PM_2.5 increased from 50.60 % in 2000 to 63.81 % in 2014 across China, which was then reduced drastically to 34.03 % in 2020. Overall, the generated LGHAP dataset has great potential to trigger multidisciplinary applications in Earth observations, climate change, public health, ecosystem assessment, and environmental management. The daily resolution AOD, PM_2.5, and PM₁₀ datasets are publicly available at https://doi.org/10.5281/zenodo.5652257 (Bai et al., 2021a), https://doi.org/10.5281/zenodo.5652265 (Bai et al., 2021b), and https://doi.org/10.5281/zenodo.5652263 (Bai et al., 2021c), respectively. Monthly and annual datasets can be acquired from https://doi.org/10.5281/zenodo.5655797 (Bai et al., 2021d) and https://doi.org/10.5281/zenodo.5655807 (Bai et al., 2021e), respectively. Python, MATLAB, R, and IDL codes are also provided to help users read and visualize these data.

1 Introduction

Atmospheric aerosols not only impact regional climate by changing the Earth radiation budget but also significantly influence air quality at the ground level (Fuzzi et al., 2015; Gao et al., 2018; Shen et al., 2020; Sun et al., 2015; Yang et al., 2020; Zheng et al., 2020). Monitoring aerosol loading in the atmosphere is thus of great significance for climate change attribution and haze pollution assessment. Aerosol optical depth (AOD), an indicator of aerosol bulks distributed within a column of air from the Earth's surface to the top of the atmosphere, has been monitored for decades to map global aerosol loading in the atmosphere. Compared with sparsely and unevenly distributed ground-based aerosol monitoring stations (e.g., the Aerosol Robotic Network – AERONET), satellite instruments can map AOD with vaster spatial coverage at even sub-hourly sampling frequency (e.g., geostationary satellite). An overview of sensors, algorithms, and AOD datasets that are widely used in the community can be found in the literature such as in Sogacheva et al. (2020) and Wei et al. (2020).

Due to negative impacts of bright surfaces (e.g., snow cover) and clouds, as well as algorithmic restrictions, satellite AOD retrievals often suffer from extensive data gaps, significantly reducing the downstream application potential such as mapping particulate matter (PM) concentrations at the ground surface (e.g., Bai et al., 2019a; J. Wei et al., 2021​​​​​​​). Also, excessive data gaps in AOD imageries may result in large uncertainty when assessing aerosol impacts on weather and climate (Guo et al., 2017; Li et al., 2019; Zhao et al., 2020; Zheng et al., 2018). Over the years, versatile gap-filling methods have been developed (e.g., Bai et al., 2016, 2020b; Chang et al., 2015). Nonetheless, filling data gaps in satellite-based AOD retrievals is still challenging due to extraordinary nonrandom missing values and high aerosol dynamics in space and time.

Wei et al. (2020) provided a short review of methods that have been frequently applied to deal with data gaps in AOD products. In general, merging AOD data acquired from diverse instruments and/or platforms is the most popular approach to improve AOD spatial coverage (Sogacheva et al., 2020). Statistical methods such as linear regression (Bai et al., 2019a; Wang et al., 2019; Zhang et al., 2017), inversed variance weighting (Chen et al., 2018; Ma et al., 2016; Sogacheva et al., 2020), and maximum likelihood estimation (Xu et al., 2015) are often applied to account for systematic bias among different datasets. Data fusion methods such as Bayesian maximum entropy can be applied to blend AOD products with different resolutions (Tang et al., 2016; X. Wei et al., 2021). Another way is to reconstruct missing AOD values using either neighboring observations in space and time or external data sources such as AOD simulations from numerical models (Li et al., 2020; Xiao et al., 2017) and even meteorological factors (Bi et al., 2018).

Although there exist a variety of gap-filling methods, spatially gap-free AOD datasets are still rare, particularly high-resolution AOD datasets from satellites, significantly limiting downstream applications such as PM_x concentration mapping. In spite of versatile PM_2.5 concentration prediction models (e.g., Di et al., 2019; Fang et al., 2016; Hu et al., 2014; Li et al., 2016; Lin et al., 2016; Liu et al., 2009; B. Wang et al., 2021), to date, there are few publicly accessible PM_x concentration datasets that can be used to examine haze pollution variations regionally and globally. Several typical datasets, e.g., the one generated by Dalhousie University (van Donkelaar et al., 2010, 2016), CHAP (J. Wei et al., 2021), and TAP (Geng et al., 2021), have been widely applied to advance our understanding of aerosol impacts across China and the globe. However, these datasets more or less still suffer from drawbacks in terms of spatial and/or temporal resolution, spatial coverage, and data accuracy. To meet contemporary needs, Zhang et al. (2021) provided a more comprehensive review of the widely used PM_x concentration mapping approaches. With a thorough review of PM_2.5 concentration mapping techniques, an optimal full-coverage PM_2.5 modeling scheme was proposed, in which diverse aerosol datasets were fused toward a full-coverage AOD map based on a multimodal approach (Bai et al., 2022). In parallel with these efforts, some have attempted to improve AOD data coverage over space with high accuracy by merging AODs observed at adjacent times directly (Li et al., 2022).

With such prior knowledge, the current study developed a big data analytics framework for generating the Long-term Gap-free High-resolution Air Pollutant concentration dataset (abbreviated as LGHAP hereafter), aiming at providing gap-free AOD, PM_2.5, and PM₁₀ concentration data with a daily 1 km resolution in China for the period of 2000 to 2020. To achieve such a goal, multimodal aerosol data acquired from diverse sources including satellites, ground stations, and numerical models were synergistically integrated via high-order singular value decomposition (HOSVD) to form a tensor-flow-based data fusion framework in the current study. Full-coverage PM_2.5 and PM₁₀ concentration data were then estimated on the basis of the gap-filled AOD dataset. This 21-year-long gap-free high-resolution (daily and 1 km) aerosol dataset was then compared against ground-based AOD and PM_x observations to validate the data accuracy of each product, particularly their performance in spatial pattern recognition and temporal trend assessment. These advances enabled a better assessment of long-term variability in haze pollution in China as well as of the corresponding population exposure over the past 2 decades.

2 Data sources

Table 1 provides a brief summary of the multisource datasets used in this study to generate the LGHAP dataset. As shown, 6 satellite-based AOD products, 5 numerical simulations of AOD and aerosol components, 11 meteorological factors, and 6 datasets of ground-based AOD and air pollutant concentration measurements, as well as a set of land cover, topographic, and socioeconomic parameters, were employed. Descriptions of these datasets are given in the following subsections.

Table 1 Summary of the data sources used in this study to generate gap-free high-resolution AOD and PM_x concentration datasets. NDVI denotes the normalized difference vegetation index.

2.1 Gridded aerosol products

In many previous studies, coarse AOD and/or aerosol component simulations acquired from numerical models have oftentimes been used as the primary data source to help derive full-coverage AOD and/or PM_2.5 concentration maps (e.g., Park et al., 2020; Y. Wang et al., 2021). However, due to the lack of a high-accuracy near-real-time emission inventory, simulated AOD and/or aerosol components are often prone to large uncertainty, which could be inevitably introduced into the final PM_2.5 estimations if no observational data are applied for possible bias correction. In such a research context, here we used six satellite-based AOD products with a relatively long temporal coverage (>5 years) to help better reconstruct historical AOD variations over space and time, though geostationary satellites can provide AOD observations at even an hourly resolution. The reasons are twofold. On the one hand, the operational AOD product from the recent Chinese FY-4 satellite is still unavailable. On the other hand, the AOD product from Himawari-8 cannot provide observations in the northwest region of China.

The latest AOD product derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) on board Terra using the multiangle implementation of atmospheric correction (MAIAC) algorithm (Lyapustin et al., 2011, 2018), was hereby used as the baseline dataset for the generation of gap-free AOD maps. This AOD product has not only a finer spatial resolution (1 km) but also a comparable and perhaps even better accuracy, when comparing it with those derived from the Dark Target and Deep Blue algorithms (Goldberg et al., 2019; Lyapustin et al., 2018). In addition, AOD products derived from MODIS on board Aqua, the Multi-angle Imaging SpectroRadiometer (MISR) on board Terra, the Visible Infrared Imaging Radiometer Suite (VIIRS) on board Suomi NPP, the Advanced Along-Track Scanning Radiometer (AATSR) on board Envisat, and POLarization and Directionality of the Earth's Reflectances (POLDER) on board PARASOL were also employed. The ultimate goal was to reduce the bias level in the final full-coverage AOD product by providing as many observational AODs as possible. Accuracies of these AOD products have been extensively validated in previous studies, e.g., de Leeuw et al. (2018), Xiao et al. (2016), Wei et al. (2019)​​​​​​​, and Che et al. (2019), to name a few. A brief description of these satellite-based AOD products can be found in Sect. S1 in the Supplement​​​​​​​.

In addition to satellite-based AOD products, numerically simulated aerosol diagnostics from MERRA-2, including AOD and aerosol components such as black carbon, organic carbon, dust, and sulfate, were also applied to help reconstruct missing AOD information and to predict PM_2.5 and PM₁₀ concentrations at the ground level. The aerosol components were used here as a proxy for an emission inventory when predicting PM_x concentrations. Big data analytics procedures applied to these datasets will be described in Sect. 3.

2.2 In situ AOD and air quality measurements

AOD observations from AERONET were used as the ground truth to evaluate the data accuracy of the generated gap-free AOD product and also as the learning target to infer AOD from air pollutant concentration and atmospheric visibility. Considering few valid data were provided in the Level 2.0 dataset, here we used the Level 1.5 AOD data to guarantee adequate in situ AOD data coverage in space and time. To validate the gridded AOD products in this study, each in situ AOD observation was registered with the gridded mean AOD over a 50 × 50 km window.

Near-surface air pollutant concentrations including PM_2.5, PM₁₀, NO₂, and SO₂ that were sampled at state-controlled monitoring sites were also applied, not only to help establish machine-learned regression models for PM_x prediction (PM_2.5 and PM₁₀) but also to infer AOD over air quality monitoring sites given their dense distributions across China. The gauged air pollutant concentration data have been released online on an hourly basis by the China National Environmental Monitoring Centre since late 2013. For quality control, outliers were first detected and removed from each pollutant dataset by following the criteria used in our previous study (Bai et al., 2020a). The missing values were then reconstructed using the diurnal-cycle-constrained empirical orthogonal function (DCCEOF) method proposed in Bai et al. (2020b).

The 3 h resolution atmospheric visibility data acquired from 4052 weather stations, at which in situ air quality measurements were not available, were employed to help generate gap-free AOD maps before 2014. Previous studies have attempted to predict PM_2.5 concentration from atmospheric visibility data with good accuracies (Liu et al., 2017), indicative of great potential for estimating AOD. Specifically, visibility data were used as an important predictor for site-specific AOD prediction, and the resulting AOD predictions were then used as critical prior information for reconstructing AOD distributions over space, especially over those regions without satellite AOD observations. Given the availability of abundant air quality measurements and the fact that automatic visibility sensors have been widely used across China since 2014, atmospheric visibility data after 2014 were thereby excluded to guarantee data consistency (Li et al., 2018​​​​​​​). For quality control, the consistency of visibility data was examined using an outlier detection method; i.e., the annual mean should not exceed 3 times the standard deviation of data over a 5-year time window (Zhang et al., 2020). Those data with apparent jumps and drifts in visibility time series were excluded. Meanwhile, visibility data on rainstorms and foggy days were eliminated as well.

2.3 Auxiliary data

As shown in Table 1, 11 meteorological factors, including air temperature at the near surface, wind speed and direction, relative humidity, surface pressure, boundary layer height, total column water vapor, downward solar radiation, and instantaneous moisture flux, were used to help resolve nonlinear relationships between PM_x and AOD, as well as to downscale AOD from MERRA-2. These data were acquired from the fifth-generation ECMWF atmospheric reanalysis (ERA5), and the first three factors were extracted at the levels of not only the ground surface but also 850 and 500 hPa so as to indicate the vertical structure of the atmosphere. Additionally, population data from WorldPop, land cover from CLCD during 2000 to 2019 (Yang and Huang, 2021) and GlobeLand 30 in 2020 (Chen et al., 2014), and elevation data from the Global Digital Elevation Model (GDEM) version 2, as well as the monthly composited 1 km normalized difference vegetation index (NDVI) from MODIS were employed to resolve the socioeconomic and ecological contributions to haze pollution. Properties of these datasets can be found in Table 1, and datasets with a finer resolution were upscaled to 0.01^∘ via a cubic interpolation method.

Figure 1 Flowchart of the proposed big data analytics framework for generating the long-term gap-free high-resolution air pollutant concentration dataset (LGHAP), taking aerosol optical depth (AOD) and PM_x (PM_2.5 and PM₁₀) concentration in China as illustration. HOSVD is an abbreviation of high-order singular value decomposition. MET, LULC, DEM, and POP denote variables of the meteorology, land use/land cover, digit elevation model, and population, respectively.

3 Methodology

With the aim of the generation of LGHAP aerosol datasets to advance environment management and Earth system science analysis, here we developed a big data analytics framework via a seamless integration of tensor-flow-based multimodal data fusion with ensemble-learning-based PM_x concentration estimation. The proposed method transformed a set of data tensors of AOD and other related datasets such as air pollutant concentration and atmospheric visibility, which were acquired from diverse sensors or platforms via integrative efforts of spatial pattern recognition for high-dimensional gridded data analysis to achieve data fusion and multiresolution image analysis, as well as knowledge transfer in statistical data mining. The proposed method consists of three major procedures in general, including multisensory data homogenization, tensor-flow-based AOD reconstruction, and ensemble learning for PM_x concentration estimation. The analytical framework of the big data analytics is depicted in Fig. 1 and described in detail in the following subsections.

3.1 Multisensory data homogenization

Since a set of aerosol products with different types, resolutions, and accuracies were applied to support the reconstruction of gap-free AOD imageries, harmonizing cross-platform biases and scale differences between these diverse datasets is crucial to multisensory data integration. In this study, machine-learned regression models were established to harmonize these heterogeneous aerosol datasets. A baseline dataset was first selected to be used as the learning target, and other datasets were calibrated to the level of baseline dataset to make them comparable. Given the finer resolution and higher proportion of data coverage in space and time, the MAIAC AOD product from Terra (AOD_Terra) was selected as the baseline dataset. Consequently, six machine-learned regression models were established between AOD_Terra and each gridded AOD product (i.e., five satellite-based AOD products plus MERRA-2 AOD simulations) using the random forest method. Meteorological factors (MET), land cover types (LULC), topography (DEM), and population (POP) were used as covariates to help downscale these multimodal AOD products to have a resolution the same as AOD_Terra while accounting for cross-mission biases arising from temporal and algorithmic differences.

Considering data gaps are extensive in satellite AOD products, especially over regions with vast cloud cover, providing prior AOD information over such regions is thus of great value in support of the reconstruction of missing AOD values. As indicated in our recent studies, AOD can be accurately predicted from ground-measured air pollutant concentration, showing an accuracy even higher than some satellite AOD retrievals (Li et al., 2022; Bai et al., 2022). To support AOD reconstruction over regions with fewer or even without valid satellite AOD observations, we attempted to infer AOD over air quality monitoring sites from in situ air pollutant concentration measurements via a machine learning approach. Similarly, machine-learned regression models were established using random forest by taking AOD_Terra as the learning target while ground-measured air pollutant concentrations, meteorological factors, land cover, and terrain information were used conjunctively as predictors.

The transformation of ground-measured air pollutant concentration data to AOD allows for providing external observational AOD data to supplement satellite observations, especially over regions suffering from significant data gaps. Since air pollutant concentration data were not available before 2013, atmospheric visibility data sampled at dense weather stations were used as an alternative for site-based AOD prediction, by applying a similar prediction model to that described above and used for air pollutant concentration. Figure S1 show the ground-based validation results of AOD inferred from atmospheric visibility and air pollutant concentrations, indicative of a generally good accuracy of these inferred AOD values. All efforts led to aggregating a set of multimodal aerosol data with different properties for multisensory data fusion to achieve gap-free AOD mapping as the next step.

3.2 Tensor-flow-based AOD reconstruction

The core of generating full-coverage AOD imageries is to fill in data gaps in AOD_Terra. Previous studies have demonstrated that merging satellite AOD retrievals at adjacent time steps can help improve the observational AOD coverage at each single snapshot, while the involvement of numerical AOD simulations can help bridge AOD data gaps (Li et al., 2022; Bai et al., 2022). In this study, a tensor completion method was particularly designed and applied to fulfill the gap filling in AOD_Terra. Specifically, the incomplete AOD_Terra imageries were deemed the hard data (true AOD state), while other AOD datasets (e.g., the downscaled AOD datasets and site-specific AOD predictions inferred from air pollutant concentration and atmospheric visibility) were used as the soft data (complementary data) to help reconstruct AOD distribution in AOD_Terra via tensor-flow-based pattern recognition. Detailed procedures for gap filling are outlined as follows.

3.2.1 Initial AOD tensor construction

Due to extensive data gaps in satellite-based AOD retrievals, it is insufficient to reconstruct all missing AOD information in AOD_Terra for a given date by simply merging the harmonized satellite-based AOD data synchronously. To fulfill AOD gap filling, the tensor completion method was thus applied to synergistically integrate AOD acquired from diverse sources. Consequently, creating the data tensor of AOD is of critical importance. In this study, the data tensor of AOD was constructed by incorporating not only observational AOD from both satellites and those inferred from in situ air quality indicators on the same date but also historical AOD retrievals from MODIS instruments (AOD_Terra and AOD_Aqua) and part of data from the downscaled MERRA-2 AOD (denoted as AOD_M2 hereafter). The latter two were applied to provide knowledge of AOD distributions over space to guide the reconstruction of missing values in AOD_Terra.

For the screening of historical observations resembling AOD_Terra distribution on the given date to be reconstructed, AOD_M2 was used in concert with AOD_Terra and site-based AOD estimations to identify similar imageries. Toward this goal, site-specific AOD estimations and 5 % randomly selected downscaled AOD_M2 data were merged directly with valid AOD_Terra to form a new image on each date. Subsequently, correlations and biases were estimated between AOD_Terra on the given date to be reconstructed and each newly merged historical AOD_Terra image. To avoid the inclusion of imageries with distinct variation patterns, only those closely resembling AOD_Terra on the date to be reconstructed were finally retained in terms of their correlations and biases subject to thresholds of R>0.7 and RMSE <0.2. Once sufficient historical imageries were obtained, the data tensor of AOD was constructed by compiling the observed AOD imageries on the given date with historical imageries into a three-dimensional data array $\mathbf{A}\in {\mathbf{R}}^{N_{\mathrm{1}}\times N_{\mathrm{2}}\times N_{\mathrm{3}}}$ (composed of N₃ images with a size of N₁×N₂). Considering satellite AOD retrievals suffer from extensive data gaps, we injected data values of site-specific AOD estimations and 1 % randomly selected downscaled AOD_M2 data directly onto grids where AOD_Terra values were missing on each specific date as prior knowledge. This not only accelerates convergence speed during the reconstruction process but also avoids large reconstruction errors over regions with tremendous data gaps in satellite-observed AOD imageries.

Figure 2 Scatterplots between ground-observed and satellite-based AOD data in different regions of China. (a–e) Original Terra MODIS AOD, (f–j) reconstructed AOD, and (k–o) combined AOD between original and reconstructed data. BTH, YRD, SC, and WC refer to the regions of Beijing–Tianjin–Hebei, Yangtze River Delta, South China, and West China, respectively.

3.2.2 Gap filling via tensor completion

Previous studies have well demonstrated the good performance of matrix decomposition methods such as empirical orthogonal functions and singular value decomposition (SVD) for missing value imputation (Bai et al., 2020b; Beckers and Rixen, 2003). However, these methods can only work on a two-dimensional matrix mathematically, namely the matrix domain. To integrate spatial features of AOD revealed by datasets to generate a smooth AOD distribution with complete coverage, in this study, the HOSVD, a specific orthogonal Tucker decomposition, was applied. More detailed descriptions of HOSVD can be found in the literature such as in Sun et al. (2021), Tucker (1966), Kolda and Bader (2009), and Sidiropoulos et al. (2017).

In Algorithm 1, we provided a stepwise description of the algorithm used to fill data gaps in AOD_Terra by integrating AOD features recognized in different imageries as the data tensor of AOD via HOSVD. To initiate the tensor decomposition, grids with missing values in the original AOD tensor were first filled with the spatial average of valid AOD data in each individual image. Then, the AOD tensor was decomposed along each of the three dimensions, while the dominant features in each dimension determined by the corresponding rank values were applied to reconstruct the data tensor. By gradually increasing the rank values and iteratively updating the initial filled values, the tensor can be reconstructed to better delineate AOD distribution over space after several iterations.

To confirm the convergence, a small portion of observational AOD values were randomly held out in advance, and the reconstructed values over these grids in each iteration were compared with these held-out data till the difference between them was lower than 0.01 (a threshold to determine convergence, a.k.a ε₁ in Algorithm 1). Meanwhile, to make the computational burden manageable, the study region (China in this study) was divided into 40 subregions (refer to Fig. S2 for the spatial distribution of these subregions), and the tensor completion was then performed over each individual region. Finally, the reconstructed imageries were mosaicked to attain a national gap-free AOD map on each specific date. During this step, a smooth filter was applied to solve the boundary effect when mosaicking two adjacent maps. Specifically, the data value on each overlapped grid at the boundary (50 km on the edge of a subregion) was averaged via an inverse distance (the distance to the edge) weighting scheme. In the end, the mosaic AOD_Terra image was retained as the final gap-free AOD product.

3.3 PM_x concentration estimation

In this study, the widely used random forest method was applied to establish regression models for PM_2.5 and PM₁₀ concentration estimation. Ground-measured PM_2.5 (or PM₁₀) concentration data were used as the learning target, and gap-filled AOD, aerosol components (AER_comp), meteorological factors (MET), the digital elevation model (DEM), the NDVI, land cover information (LC), and population were used as regressors. The random forest (RF) regression model can be generally formulated as

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathrm{PM}}_x& =\mathrm{RF}(\mathrm{AOD},{\mathrm{AER}}_{\mathrm{comp}},\mathrm{MET},\mathrm{DEM},\mathrm{NDVI},\mathrm{POP},\\ & \mathrm{LC},\mathrm{month}),\end{array}\end{array}$$

where month is a categorical variable that was used to account for monthly varying relationships between AOD and PM_x. For validation, PM_2.5 and PM₁₀ measurements from 10 % of monitoring sites were randomly held out to evaluate the predictive performance of each regression model. During the training process, 500 regression trees were used in each RF model, and each tree was grown on a bootstrap sample. The learning dataset was randomly divided into two parts during the training process, with 80 % used as the training set and the remaining 20 % for testing. In order to guarantee a larger value of PM₁₀ than PM_2.5, PM_2.5 estimations from Eq. (1) were used as one predictor in addition to factors used to predict PM_2.5 when estimating PM₁₀ concentration. Such a model can also significantly improve the prediction accuracy of PM₁₀ given the prior PM_2.5 information.

3.4 Point–surface data fusion

Ground-measured PM_2.5 and PM₁₀ concentration data were further fused with their gridded estimations to enhance the data accuracy of PM_x data after 2014. Here, the well-known optimal interpolation (OI) method was applied to perform point–surface fusions between two different types of dataset. Please refer to Bai et al. (2022) and Li et al. (2022) for a more detailed description of the OI method used to fuse PM_x concentration data. In this study, a modified scheme was developed to select neighboring observations. To avoid an isotropic interpolation effect, here we only used 30 ground observations with land cover, terrain, and atmospheric conditions similar to those at the analyzed grid cell to estimate the innovation that should be assigned to the background value at the given grid. In other words, a similarity measure was first estimated between the analyzed grid cell and neighboring sites in terms of land cover, the DEM, and atmospheric conditions. The 30 observations with similar background fields were then used in the OI procedure to correct possible bias in gridded PM_x estimations. Such a treatment can help exclude those observations with different ambient background, e.g., one site not far from the given grid but separated by a high mountain, thereby avoiding the possible propagation of antiphase corrections to data over adjacent grids.

4 Results and discussion

4.1 Data accuracy of gap-free AOD in LGHAP

Table 2 summarizes the data accuracy of the gap-free AOD dataset generated in this study. For comparison, the data accuracy of each original AOD dataset was also assessed. Since in situ AOD measurements were not used as data input when reconstructing missing AOD information, the gap-free AOD can be directly compared with in situ AOD measurements from AERONET. As indicated, all these AOD datasets are in good agreement with in situ AOD measurements. Generally, AODs from MODIS on board Terra and Aqua have an almost identical data accuracy, which is also among the highest when comparing with other datasets (R=0.95, and RMSE = 0.14). AODs from AATSR show a comparable accuracy to that of MODIS but with a relatively low correlation with ground-based AOD measurements. AODs from MISR, POLDER, and VIIRS exhibit a similar bias level, with R varying from 0.80 to 0.92 and RMSE ranging from 0.22 to 0.29. In contrast, AOD_M2 data have the poorest accuracy among these eight gridded AOD datasets (R=0.77, and RMSE = 0.36), even though AOD data from AERONET and satellite observations like MODIS had been already assimilated. This indicates the presence of large biases in AOD_M2, and thus these AOD_M2 data cannot solely be used to delineate AOD distributions over space.

Table 2 Data accuracy of original and gap-free AOD datasets used and/or generated in this study. The expected error (EE) was defined as $\pm \mathrm{0.05}+\mathrm{0.15}$ × AOD_site. MAE denotes the mean absolute error and was calculated as the absolute value of bias.

Compared to the first seven gridded AOD datasets, the LGHAP AOD dataset has an accuracy slightly worse than the original MODIS AOD product but comparable to AODs from MISR, POLDER, and MERRA-2, with an R of 0.91 and RMSE equaling 0.21 compared to ground-based AOD observations. Nevertheless, the gap-filled AOD appeared to overestimate ground-based AOD observations, and this could be due to the involvement of AODs from VIIRS and POLDER as these two products significantly overestimated ground AOD observations, which can be indicated by the proportion of data pairs above the expected error (EE). On the other hand, significant underestimations in AOD_M2 were not introduced into the LGHAP AOD as the below-EE ratio was 32.97 % in the former and was only 12.27 % in the latter. These results indicate that the LGHAP AOD data are more likely to resemble AOD distributions revealed by satellite observations rather than AOD_M2, endorsing the advantages of involving multisensory satellite AOD observations to support missing AOD reconstruction. Figure 2 further compares the data accuracy of original AOD_Terra and the reconstructed data over different regions of China. It is indicative that the purely reconstructed data have an accuracy (R=0.88, and RMSE = 0.26) lower than the original AOD_Terra (R=0.95, and RMSE = 0.13) across China, especially in South China where the reconstructed data significantly underestimated the ground-based AOD observations. Possible reasons for this effect could be attributed to extensive data gaps in satellite AOD retrievals due to frequent and extensive cloud cover there (refer to Fig. S3 for the distribution of the mean data integrity of AOD_Terra during 2000–2020), and the scarce AOD observations significantly limit the learning capacity in space and temporal domain during the tensor completion process. In other words, limited observations in satellite imageries greatly reduced the learning performance from the sparse tensor. Nevertheless, the purely reconstructed data exhibit a bias level comparable to AOD retrievals from several satellite instruments, e.g., MISR, VIIRS, and POLDER. This demonstrates the good performance of the proposed tensor completion method in reconstructing missing AOD information. By combining the reconstructed data with original AOD_Terra, we obtained a 21-year-long gap-free high-resolution (daily and 1 km) AOD dataset with satisfying accuracy (R=0.91, and RMSE = 0.21).

In Fig. 3 we present a comparison of AOD time series between the LGHAP dataset and ground observations at three AERONET sites under different air pollution levels. As shown, the AOD time series from LGHAP are temporally continuous, whereas data gaps are common in AERONET observations. Generally, AODs from LGHAP are well reconstructed with respect to the temporal variations in aerosol loading at these three sites, with R ranging from 0.77 to 0.90 and RMSE varying between 0.11 and 0.21. For illustration, Fig. 4 compares the spatial distribution of original and gap-filled AOD on four dates with different AOD_Terra coverage over space. As shown, the missing AOD values were well reconstructed after gap filling, resembling a smooth and reasonable AOD distribution over space, even over regions with very limited prior AOD observations from Terra MODIS (e.g., Fig. 4d). As indicated in Fig. 4a and c, the high AOD loading was also properly reconstructed even though no prior information was provided by AOD_Terra. Since AERONET AOD observations were not used as a data input when generating the LGHAP AOD dataset, these independent validation results clearly demonstrated the high accuracy of the LGHAP AOD product as well as a good performance of the proposed AOD gap-filling approach.

Figure 3 Comparison of monthly AOD time series from LGHAP and AERONET at three different stations in China. Latitude and longitude information of each site is given in parentheses.

Figure 4 Spatial patterns of the reconstructed AOD under different baseline AOD coverage ratios. In each panel, the upper part presents the original AOD distribution from Terra MODIS and the gap-filled imagery is shown below. The zoomed-in views of the outlined regions are shown on the right. Dates are given in the format year/month/day.

Since the final gap-free AOD product was generated mainly by integrating a set of data tensors of gridded AOD with AOD estimations from in situ air quality measurements, the relative contribution of each product to the final gap-free dataset is worth investigating. In this study, a data-coverage-ratio-weighted nonlinear correlation coefficient was proposed to examine the relative contribution of each gridded product to the LGHAP AOD dataset. The nonlinear correlation coefficient was used to assess the mutual information between two variables (Sun et al., 2021; Wang et al., 2005), while the data coverage ratio was multiplied to indicate the overall contribution of one product to the final fused dataset (refer to Sect. S2 for the definition of this indicator). As shown in Fig. 5, the relative contribution of each gridded product varied with time and the input data sources. In the 2 early years (2000–2001), the AOD distribution in gap-free imageries was determined largely by AOD_Terra (81 %), whereas this ratio decreased to about 30 % when many other products were involved, especially AOD from Aqua and PARASOL. With the advent of VIIRS and the loss of PARASOL after 2012, the relative contribution changed drastically as AOD from MODIS and VIIRS played the dominant roles in reconstructing AOD distribution. Note the relative contribution of AOD_M2 remained lower than 10 %, indicative of the greater importance of satellite observations in generating the LGHAP AOD product.

Figure 5 Spatiotemporal variations in the relative contribution of each gridded AOD product to the generation of the LGHAP AOD dataset. The relative contribution was estimated as the data-coverage-ratio-weighted nonlinear correlation coefficient (please refer to Sect. S2 in the Supplement for the arithmetic theory to calculate this measure). The annual mean shown outside is the national averaged contribution in each individual year, while the regional mean shown on the map was averaged over the past 21 years in each subregion.

With respect to the temporally averaged contribution in each subregion, it shows that the relative contribution of each product also varied significantly across regions. Generally, AOD from MODIS aboard Terra and Aqua played the most important role (>60 %) in generating the LGHAP AOD product, except over the southwest part of the country (Tibetan Plateau) where AOD_M2 contributed most. This is largely associated with the fact that data gaps are abnormally high in satellite observations over this region because of the vast and long-lasting snow cover (refer to Fig. S3 for the data integrity distribution). Consequently, AOD_M2 would play an important role in reconstructing AOD distribution over such regions. Note that the relative contribution of AOD estimations from in situ air quality measurements was not accounted for in the current analysis because of incomparable spatial coverage of in situ data in contrast to gridded AOD products, and this does not imply the contribution of in situ AOD estimations was negligible. Overall, the results shown here clearly highlight the success of big data analytics in generating the LGHAP AOD dataset via integrative efforts from diverse data sources.

4.2 Data accuracy of PM_2.5 and PM₁₀ estimations

By taking advantage of the gap-filled AOD, daily 1 km resolution PM_2.5 and PM₁₀ concentration data in China were then estimated via an ensemble learning approach. Figure S4 shows the sample-based cross-validation accuracy of two prediction models. It shows that the original daily PM_2.5 prediction model had a sample-based cross-validation R² of 0.79 and RMSE of 20.04 µg m⁻³. This accuracy is comparable to that of our previous study (Bai et al., 2019a) but slightly worse than those reported in some recent studies (Table 3). In contrast, PM₁₀ had a much higher prediction accuracy, with an R² of 0.90 and RMSE of 21.06 µg m⁻³ for the daily product. This good performance should be attributed to the involvement of PM_2.5 estimations as a predictor in the PM₁₀ prediction model. Figure 6 shows the site-specific (held-out in advance) validation accuracy of daily, monthly, and annual mean PM_2.5 and PM₁₀ concentration in LGHAP. As shown, the site-specific validation results indicated that the final full-coverage (gap-free) daily PM_2.5 and PM₁₀ concentration data are in good agreement with ground-based measurements, with an R of 0.95 and RMSE of 12.03 µg m⁻³ for PM_2.5 and an R of 0.94 and RMSE of 19.56 µg m⁻³ for PM₁₀. Overall, PM_x data in LGHAP not only are spatially complete with a finer resolution but also have a comparable accuracy with previous studies.

Table 3 Comparison of the data quality of PM_2.5 from LGHAP with other related studies.

Figure 6 Scatterplots between observed and estimated PM_2.5 and PM₁₀ concentration. Panels (a)–(c) denote daily, monthly, and annual mean PM_2.5 validation results, respectively, while panels (d)–(f) are for PM₁₀ concentration. The ground measurements were acquired from 30 independent air quality monitoring sites that were randomly held out before the model training.

Figure 7 presents a 2-year-long comparison of PM_2.5 concentration time series from LGHAP and two other open-access datasets with PM_2.5 measurements sampled at four United States embassies in China. Since this ground-based dataset has been seldom noticed and used, it can be applied as an independent dataset to fairly evaluate the accuracy of these three machine-learned PM_2.5 estimations. As shown, all these three datasets well reconstructed temporal variations in PM_2.5 from 2019 to 2020. Temporally, LGHAP and TAP are continuous, while CHAP suffers from significant data gaps because no gap filling was applied when generating the dataset. Compared with the other two datasets, LGHAP PM_2.5 data had a better agreement with ground-based PM_2.5 measurements. This high accuracy could be partially due to the fusion of in situ PM_2.5 data measured at adjacent sites via the OI method. Figure S5 compares PM_2.5 time series from LGHAP with PM_2.5 measurements sampled at five United States embassies in China. It is indicative that historical PM_2.5 variations over these five cities were well reconstructed in LGHAP, even over years before 2014 during which PM_2.5 measurements from state-control monitoring sites were not available. Note PM_2.5 estimations appeared to significantly underestimate PM_2.5 concentration sampled at the embassy in Beijing before 2013. Considering the reconstructed AOD time series agreed well with AERONET AOD in Beijing (Fig. 3a) and the model performed well in predicting historical PM_2.5 in Shanghai during the synchronous time period (Fig. S5b), we are more willing to attribute this issue to significant PM_2.5 overestimations by the US Embassy during that period. Overall, these independent validation results collectively indicate a good accuracy of PM_2.5 in the LGHAP dataset.

Figure 7 Comparison of PM_2.5 concentration time series between LGHAP (red line) and two open datasets (blue: TAP; green: CHAP). Here, hourly PM_2.5 concentrations measured by four United States embassies in China from 2019 to 2020 (grey bar) were used as an independent PM_2.5 dataset to validate these three daily products. CHAP and TAP are two open-access datasets providing PM_2.5 concentration that were created by J. Wei et al. (2021) and Geng et al. (2021), respectively.

In Fig. 8 we compared the spatial distribution of PM_2.5 that was reconstructed by different datasets. Compared to LGHAP and TAP, PM_2.5 data from CHAP are not gap-free since the spatial coverage is determined by the AOD data coverage in the MAIAC product. Compared to TAP, LGHAP PM_2.5 data have a finer resolution (1 km versus 10 km), enabling us to examine PM_2.5 variations in space with more details. Overall, LGHAP performs better in reconstructing PM_2.5 spatial distributions than the other two datasets. Reasons could be attributed to the following two aspects. Firstly, in situ PM_2.5 measurements were fused with gridded PM_2.5 estimations using the OI method when generating the final PM_2.5 product in LGHAP. This can help correct modeling biases in original PM_2.5 estimations. Secondly, a set of satellite-based AOD retrievals were incorporated when generating the full-coverage AOD product, which greatly helps reduce large biases in numerical AOD simulations, yielding more accurate PM_2.5 estimations in turn. This also highlights the great advantages of using big data analytics methods to advance air pollution assessment.

Figure 8 Comparison of PM_2.5 distribution reconstructed by different PM_2.5 concentration datasets. From the left to right, it shows in situ PM_2.5 concentration measurements, CHAP, TAP, and LGHAP, respectively. Dates are given in the format year/month/day.

To illustrate the fine resolution of the LGHAP dataset, we compared the annual mean PM₁₀ concentration in 2019 with the proportion of impervious surface that was derived from 30 m resolution land cover data in eastern China. As shown in Fig. 9, the finer resolution of the LGHAP dataset enables us to easily recognize the “hot spot” regions with high PM₁₀ loading. By referring to the impervious surface distribution on the right, we found that these hot spots are mainly over cities and towns, indicative of the presence of pollution islands in urban regions. Owing to the involvement of such high-resolution datasets, the spatial details of PM_2.5 and PM₁₀ can be well recognized in LGHAP. The finer-spatial-resolution advantage of the LGHAP dataset can also be demonstrated by comparisons of the spatial distribution of annual mean PM_2.5 concentration that was revealed by four different datasets shown in Fig. S6.

Figure 9 Comparison of annual mean PM₁₀ concentration with the proportion of areas coved by impervious surface in eastern China.

4.3 Long-term trends of haze pollution in China from 2000 to 2020

The aerosol pollution trends in China can be better examined by taking advantage of the LGHAP dataset given its long-temporal-coverage, gap-free, and high-resolution superiority. Severe haze pollution such as PM_2.5 is oftentimes observed during the winter half year (September–February). In this study, we first calculated mean PM_2.5 concentration in China during the winter half year from 2000 to 2020. As shown in Fig. 10, severe haze pollution events were mainly observed in North China during the wintertime, especially over the adjacent region in Hebei–Shandong–Henan provinces. In addition, the Sichuan Basin and Fenwei Plain also suffered from severe haze pollution. Temporally, severe haze pollution events occurred mainly from late 2002 to early 2017 and were significantly reduced after 2017. Similar pattern can also be inferred from PM₁₀ concentration distributions shown in Fig. S7.

Figure 10 Spatial distribution of mean PM_2.5 concentration from LGHAP during the winter half year (September–February) from 2000 to 2020 in China.

Figure 11 shows the temporal variations in the proportion of land areas covered by PM_2.5 concentration exceeding 35 µg m⁻³ (the national ambient air quality standard for 24 h PM_2.5 concentration given in GB 3095-2012​​​​​​​). As shown in Fig. 11a, severe PM_2.5 pollution occurred mainly during the wintertime in China, as more than one-third of land areas (indicated by the blue lines) were exposed to unhealthy PM_2.5 pollutants. Meanwhile, an apparent inflection was observed in 2007, after which the number of episode days decreased drastically for more than one-third of the land area covered by PM_2.5 concentration exceeding 35 µg m⁻³. According to the proportion of land area covered with annual mean PM_2.5 concentration greater than 35 µg m⁻³, the variation in haze pollution in China can be generally divided into three different periods during the past 2 decades (Fig. 11b). As indicated, an increasing trend was observed from 2000 to 2007, during which land areas covered by PM_2.5 concentration greater than 35 µg m⁻³ increased to nearly 40 % at a pace of 1.04 % a⁻¹. The second period was from 2008 to 2013, during which the land area coverage ratio decreased at a rate of −0.21 % a⁻¹. The third period started from 2014, after which the land area covered with PM_2.5 concentration more than 35 µg m⁻³ decreased drastically, at a pace of −2.23 % a⁻¹.

Figure 11 Temporal variations in the proportion of land areas covered with PM_2.5 concentration exceeding 35 µg m⁻³ and PM_2.5 trends during three different periods. (a) Temporal variations in the land coverage ratio with daily PM_2.5 concentration exceeding 35 µg m⁻³ from 2000 to 2020. (b) Same as (a) but for annual mean PM_2.5 concentration. (c–e) PM_2.5 trends during the periods of 2000–2007, 2008–2013, and 2014–2020. The dotted regions imply trend estimations are statistically insignificant at the 95 % confidence interval.

Figure 11c–e present the linear trend of PM_2.5 concentration during these three specific periods, from which we observed that significant PM_2.5 variations occurred mainly over the eastern part of the country where two-thirds of the population resides. A near-ubiquitous PM_2.5 increasing trend was observed during 2000–2007, with significant increase (>1.0 µg m⁻³ a⁻¹) mainly observed in eastern China. During the second period, PM_2.5 concentration over most regions showed a small decreasing trend except in the Hebei–Shandong–Henan region where an increasing trend was still observed. An apparent decreasing trend was observed over most parts of the country after 2014, indicative of significant reductions in PM_2.5 loading across China. This trend distribution is in line with our previous finding that was derived using the annual mean PM_2.5 concentration dataset generated by Dalhousie University (Bai et al., 2019b). However, differences were still observed in terms of the regions where significant decreasing trends were present. The most significant decreasing trends were mainly observed in the Sichuan Basin and Pearl River Delta in the previous study. However, regions with drastic PM_2.5 decreases were found mainly in North China where severe haze pollution events were oftentimes reported. Similar variation patterns can also be inferred from PM₁₀ (Fig. S8) and AOD (Fig. S9). Overall, the LGHAP dataset provides us with a gridded perspective to better examine long-term variations in haze pollution in China during the past 2 decades.

Figure 12 Spatial distribution of population-weighted PM_2.5 concentration and the proportion of the population exposed to PM_2.5 concentration greater than 35 µg m⁻³. Annual and daily LGHAP PM_2.5 concentration data were used for the calculation of weighted PM_2.5 and the proportion of population exposure, respectively. The diamond and red line indicate the annual mean and median population proportion, respectively.

4.4 Population exposure to PM_2.5 pollution in China

By taking advantage of fine-resolution LGHAP PM_2.5 concentration and gridded population data, population exposure to PM_2.5 pollution across China over the past 2 decades was estimated. Figure 12 shows the spatial distribution of population-weighted PM_2.5 concentration and the proportion of population exposed to PM_2.5 concentration greater than 35  µg m⁻³. As shown, the spatial distribution of population-weighted PM_2.5 concentration resembles the spatial pattern of annual mean PM_2.5 concentration, with high values observed mainly in eastern and central China as well as northwest China. Nonetheless, PM_2.5 sources in these two areas could be different. In northwest China, natural emissions could be the dominant source since a very limited number of people reside there. In contrast, most of the population lives in eastern and central China with a highly developed economy, and anthropogenic emissions thus might play more important roles in PM_2.5 formation (Xin et al., 2015; Yang et al., 2011). In regard to the proportion of the population exposed to the ambient environment with PM_2.5 concentration greater than 35 µg m⁻³, we observed that the annual mean population ratio exposure to unhealthy PM_2.5 increased gradually from 50.60 % in 2000 to 65.72 % in 2007. During 2007–2014, the ratio varied with small changes (<5 %), whereas a drastic decline was observed after 2014, when the annual mean proportion of population exposed to unhealthy PM_2.5 was reduced from 63.81 % in 2014 to 34.03 % in 2020, even though the total population increased from 1.37 billion to 1.41 billion during the synchronous period. Nonetheless, more than one-third of the population was still exposed to unhealthy PM_2.5, highlighting the requirement for further emission reduction actions to manage haze pollution in China.

5 Data availability

The LGHAP dataset, consisting of gap-free AOD, PM_2.5, and PM₁₀ concentration with a daily 1 km resolution from 2000 to 2020, is publicly accessible. All data are provided in the NetCDF format, and data in each individual year are archived in a zip file. For AOD, the dataset has a disk storage size of nearly 27 GB in total, which is available at https://doi.org/10.5281/zenodo.5652257 (Bai et al., 2021a). PM_2.5 (38 GB) and PM₁₀ (48 GB) concentration data can be acquired from https://doi.org/10.5281/zenodo.5652265 (Bai et al., 2021b) and https://doi.org/10.5281/zenodo.5652263 (Bai et al., 2021c), respectively. Additionally, monthly and annual mean datasets are also provided, which are publicly available at https://doi.org/10.5281/zenodo.5655797 (Bai et al., 2021d) and https://doi.org/10.5281/zenodo.5655807 (Bai et al., 2021e), respectively. In addition to these datasets, Python, MATLAB, R, and IDL codes that can be used to read and visualize these data are provided as well.

6 Conclusion

In this study, a big data analytics method was developed for generating the LGHAP dataset to advance research in Earth system science and environment management. With integrative efforts to fuse AOD features extracted from a set of AOD data tensors and knowledge transfer in statistical data mining from diverse air quality indicators, the LGHAP aerosol dataset providing 21-year-long (2000–2020) gap-free AOD, PM_2.5, and PM₁₀ concentration data with a daily 1 km resolution in China was generated. Gap-filled AOD imageries were firstly generated by reconstructing AOD distribution in AOD_Terra via synergistically fusing AOD features recognized from diverse satellites and numerical models as well as in situ data through tensor completion. Compared to ground-based AOD measurements, the gap-filled AOD data exhibit a satisfying prediction accuracy and good performance in delineating AOD variations over space and time. To our knowledge, this is the first attempt to generate a long-term high-resolution AOD dataset with a gap-free nature in China.

PM_2.5 and PM₁₀ concentration data were then estimated using an ensemble learning approach by taking advantage of the generated gap-free AOD imageries. Ground validation results also indicate good accuracies of these two gridded products, showing a bias level comparable with many previous studies. Compared with other open-access daily PM_2.5 concentration datasets, the LGHAP PM_2.5 dataset performs well due to the vantage of having gap-free and fine resolution products. With this gap-free and high-resolution dataset, the long-term variation trend of haze pollution in China over the past 2 decades was examined, and apparent inflections were observed in 2007 and 2014, when PM_2.5 concentration was found to turn from an increasing path to decreasing in 2007 with a more drastic decline observed starting from 2014. Moreover, the LGHAP dataset provides us with a gridded perspective to assess 2-decade-long population exposure to PM_2.5 pollution in China. In spite of a drastic decline in population exposure, there is still more than one-third of the population exposed to unhealthy PM_2.5 pollutants, highlighting the requirement of long-lasting actions to continue PM_2.5-related emission reduction.

Overall, these three gridded LGHAP aerosol products provide a long-term perspective on aerosol changes over different regions of China, and users are encouraged to use the LGHAP dataset to assess aerosol impacts on public health, air quality, climate, and ecosystems. The dataset has been publicly released online and is freely accessible via the links provided in Sect. 5. A global-scale dataset is on track and will be released to the public soon.