An optimized semi-empirical physical approach for satellite-based PM_2.5 retrieval: embedding machine learning to simulate complex physical parameters

Satellite remote sensing of PM_2.5 (fine particulate matter) mass concentration has become one of the most popular atmospheric research aspects, resulting in the development of different models. Among them, the semi-empirical physical approach constructs the transformation relationship between the aerosol optical depth (AOD) and PM_2.5 based on the optical properties of particles, which has strong physical significance. Also, it performs the PM_2.5 retrieval independently of the ground stations. However, due to the complex physical relationship, the physical parameters in the semi-empirical approach are difficult to calculate accurately, resulting in relatively limited accuracy. To achieve the optimization effect, this study proposes a method of embedding machine learning into a semi-physical empirical model (RF-PMRS). Specifically, based on the theory of the physical PM_2.5 remote sensing (PMRS) approach, the complex parameter (VE_f, a columnar volume-to-extinction ratio of fine particles) is simulated by the random forest (RF) model. Also, a fine-mode fraction product with higher quality is applied to make up for the insufficient coverage of satellite products. Experiments in North China (35^∘–45^∘N, 110^∘–120^∘E) show that the surface PM_2.5 concentration derived by RF-PMRS has an average annual value of 57.92 µg m⁻³ vs. the ground value of 60.23 µg m⁻³. Compared with the original method, RMSE decreases by 39.95 µg m⁻³, and the relative deviation is reduced by 44.87 %. Moreover, validation at two Aerosol Robotic Network (AERONET) sites presents a time series change closer to the true values, with an R of about 0.80. This study is also a preliminary attempt to combine model-driven and data-driven models, laying the foundation for further atmospheric research on optimization methods.

1 Introduction

Epidemiological studies have indicated that PM_2.5 (fine particulate matter with an aerodynamic equivalent diameter no greater than 2.5 µm) can adversely affect human health, such as increasing the risk of diabetes and respiratory diseases (Bowe et al., 2018; Pope III et al., 2002; Xu et al., 2013), and accurate surface PM_2.5 concentration is the basis of air pollution health-related research. Satellite remote sensing has the advantages of high resolution and global coverage (Ma et al., 2014; Wu et al., 2020; He et al., 2022), including variables strongly associated with PM_2.5 such as aerosol optical depth (AOD). Therefore, it has become a mainstream method for fine-particle estimation (Zhang et al., 2021).

There are three main satellite-based ways of retrieving PM_2.5.

1. Chemical transport models-based method.

   This method calculates a scaling factor η between AOD and PM_2.5 simulated by atmospheric chemical transport models (CTMs) (Lyu et al., 2022; Xiao et al., 2022) and then transfers the proportional relationship to satellite AOD data when calculating surface PM_2.5 concentration (Geng et al., 2015; Van Donkelaar et al., 2006). However, the assumption of a constant factor between simulated and observed values has large spatiotemporal limitations.

2. Univariate/multivariate regression.

   This kind of data-driven method establishes a statistical model between AOD, auxiliary variables, and ground PM_2.5 observations. Machine learning is a common tool for such regression methods due to its powerful nonlinear fitting ability between multiple variables (Irrgang et al., 2021), but the regression algorithms in machine learning are affected by the distribution and density of ground stations (Gupta and Christopher, 2009; Li et al., 2017).

3. Semi-empirical physical approach.

   Taking the physical theory as the basis, surface PM_2.5 is derived through an empirical formula constructed from AOD and some PM-related key parameters, including an important empirical parameter related to the optical properties (S). The process steps are explicit and independent of ground station observations. Meanwhile, this approach has stronger physical interpretability than the previous two methods with a large space for optimization.

Due to the complexity of the physical parameters, many studies have optimized the semi-empirical physical approach. Based on 355 nm band radar observations, Raut and Chazette (2009) introduced a specific extinction cross-section to simplify the expression of S, and PM_2.5 concentration was estimated. Kokhanovsky et al. (2009) constructed a particle-effective radius model, which can obtain the particle concentrations throughout the atmospheric column. Furthermore, Zhang and Li (2015) proposed the physical PM_2.5 remote sensing (PMRS) method. It replaced S by defining a volume-to-extinction ratio of fine particles (VE_f) and used a quadratic polynomial of fine-mode fraction (FMF) to simulate VE_f, showing certain advantages (Li et al., 2016; Zhang et al., 2020).

However, the above semi-physical empirical models have some shortcomings. Firstly, the satellite data used in the models are blocked by clouds and fog in some areas; thus high-coverage and high-precision products need to be excavated and applied. Secondly, there are still large uncertainties in estimating physical parameters (such as a simple polynomial fit to S in the PMRS method), and their expressions need to be improved. To date, machine learning (ML) has developed rapidly (He et al., 2021). It can detect complex nonlinear relationships of multiple data and model their interaction (Yuan et al., 2020; Lee et al., 2022). This provides an idea for improving the accuracy of physical parameter acquisition so as to estimate high-precision PM_2.5 through semi-physical empirical models.

According to this idea, our study proposes an optimized semi-empirical physical model (RF-PMRS) based on the PMRS theory, which attempts to explore the possibility of combining physical models and ML. To be specific, we creatively embed ML (the random forest model) into the PMRS method to simulate the physical parameter (i.e., VE_f) derived from FMF and related variables, thus optimizing the previous polynomial expression. Moreover, to further improve the PM_2.5 retrieval accuracy, the physical–deep learning FMF (Phy-DL FMF) dataset generated by a hybrid retrieval algorithm of ML and physical mechanisms is introduced. Ultimately, we comprehensively validate the performance of the PM_2.5 obtained by our optimized approach.

The remainder of our article is as follows. Section 2 describes the experimental datasets. Section 3 illustrates the specific derivation process of the proposed method. Section 4 analyzes the evaluation results. Some supporting experiments are discussed in Sect. 5. The final part provides the conclusion.

2 Data

2.1 AERONET data

The Aerosol Robotic Network (AERONET) is a federation of ground-based sun–sky radiometer networks, providing worldwide remote sensing aerosol data for more than 25 years (Holben et al., 1998). The current revision of the dataset is Version 3 (Giles et al., 2017). Due to its high quality, the data from AERONET have been regarded as theoretical true values to evaluate satellite-based products in related studies (Chen et al., 2020; Gao et al., 2016; Wang et al., 2019). AOD, FMF, and volume size distribution products with Level 2.0 (quality-assured) are applied to calculate the true values of the physical parameters and then to implement our modeling purpose (not involved in PM_2.5 calculations). A total of nine AERONET sites corresponding to four typical aerosol types participate in the training. Table 1 shows the specific information.

Table 1 Data information on nine AERONET sites classified by aerosol types. Location indicates the latitude and longitude, where a negative number means a southern latitude and a western longitude. Two sites in bold font participate in the PM_2.5 validation experiment. GSFC: Goddard Space Flight Center.

2.2 MODIS AOD

MCD19A2, the Moderate Resolution Imaging Spectroradiometer (MODIS) Collection 6 Level 2 gridded (L2G) land AOD product (Lyapustin and Wang, 2015), is selected in this study. It is derived from the Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm, which can improve the accuracy of cloud detection and aerosol retrieval (Lyapustin et al., 2011). Moreover, this new advanced algorithm jointly combines MODIS Terra and Aqua into a single sensor (Lyapustin et al., 2014). The product is produced daily with a 1 km resolution, including aerosol parameters such as 470 and 550 nm AOD, quality assurance (QA), and uncertainty factors.

The processing of MCD19A2 data (in Hierarchical Data Format, HDF) is mainly divided into five steps: AOD–QA band extraction, best-quality AOD selection, Terra–Aqua data synthesis, missing information reconstruction, and mosaic. Finally, the daily AOD distribution in GeoTiff format is obtained.

2.3 Phy-DL FMF dataset

The original global land FMF products have poor data integrity and low accuracy. To enhance their reliability, Yan et al. (2022) have released a satellite-based dataset called Phy-DL FMF, which integrates physical and deep learning methods. Specifically, it selects the FMF data obtained by a physical method (i.e., lookup-table-based spectral deconvolution algorithm, LUT-SDA) as the optimization target (Yan et al., 2017). Then it combines the Phy-based FMF into a deep learning model along with multiple auxiliary data such as satellite observations for the final Phy-DL results. Note that the process is trained with AERONET data as the ground truth. The product has a spatial resolution of 1^∘ and covers 2001 to 2020 (daily scale). In the comparison experiment against the ground FMF, Phy-DL FMF shows a higher accuracy (R=0.78, RMSE =0.100) than MODIS FMF (R=0.37, RMSE =0.282) (Yan et al., 2022).

2.4 Meteorological data

The meteorological data are obtained from the ERA5 dataset, including the values of planetary boundary layer height (PBLH) and relative humidity (RH). As the fifth-generation reanalysis product released by the European Center for Medium-Range Weather Forecasts (ECMWF), ERA5 provides atmospheric data at 0.25^∘ every hour based on the data assimilation principle (Hersbach et al., 2018). It should be noted that RH is not archived directly in ERA5 and thus should be calculated by 2 m temperature T and dew point temperature T_d (refer to https://confluence.ecmwf.int/display/CKB/ERA-Interim:+documentation#ERAInterim:documentation-Computationofnear-surfacehumidityandsnowcover, last access: 20 July 2023):

$$\begin{array}{}\text{(1)}& \mathrm{RH}=\mathrm{100}\times {\displaystyle \frac{e_{\mathrm{s}}\left(T_{\mathrm{d}}\right)}{e_{\mathrm{s}}\left(T\right)}},\end{array}$$

where e_s(t) represents the saturation vapor pressure related to a temperature t in degrees Celsius (Simmons et al., 1999) of

$$\begin{array}{}\text{(2)}& e_{\mathrm{s}}\left(t\right)=\mathrm{6.112}\times \exp \left({\displaystyle \frac{\mathrm{17.67}\times t}{t+\mathrm{243.5}}}\right).\end{array}$$

2.5 Ground PM_2.5 measurements

The North China (NC) region is chosen as the main experimental validation area for the final PM_2.5 calculations. The near-surface hourly PM_2.5 values are obtained from the China National Environmental Monitoring Centre (CNEMC). Nowadays, over 1600 ground-based monitors are working continuously and a total of 232 stations (in 2017) participate in this work. Figure 1 displays the site distributions of the NC region.

Figure 1 The location of PM_2.5 ground monitoring stations in the NC region (35–45^∘, 110–120^∘). The red points represent the PM_2.5 stations.

3 Methods

Based on the basic physical properties of atmospheric aerosols, the semi-physical empirical approach starts from the integration of PM mass concentration and AOD. Then it combines several key factors related to PM_2.5 to derive the in situ PM_2.5 concentration through multiple remote sensing variables (Koelemeijer et al., 2006). The overall empirical relationship can be represented as

$$\begin{array}{}\text{(3)}& {\mathrm{PM}}_{\mathrm{2.5}}=\mathrm{AOD}{\displaystyle \frac{\mathit{\rho }}{H\cdot f\left(\mathrm{RH}\right)}}S,\end{array}$$

where ρ denotes the particle density and H denotes the atmospheric boundary layer height. f(RH) represents the hygroscopic growth factor related to relative humidity (RH). S is an optical characteristic parameter that should be simulated.

3.1 PMRS method

3.1.1 The expression of VE_f

To illustrate S more precisely, PMRS defines the columnar volume-to-extinction ratio of fine particles (i.e., VE_f), which can be regarded as the basis of our optimization method. So Eq. (3) is transformed into

$$\begin{array}{}\text{(4)}& {\mathrm{PM}}_{\mathrm{2.5}}=\mathrm{AOD}{\displaystyle \frac{\mathit{\rho }}{H\cdot f\left(\mathrm{RH}\right)}}{\mathrm{VE}}_{\mathrm{f}}.\end{array}$$

Related to particle size, aerosol extinction, and other properties, VE_f can be expressed as

$$\begin{array}{}\text{(5)}& {\mathrm{VE}}_{\mathrm{f}}={\displaystyle \frac{V_{\mathrm{f},\mathrm{column}}}{{\mathrm{AOD}}_{\mathrm{f}}}},{\mathrm{AOD}}_{\mathrm{f}}=\mathrm{AOD}\cdot \mathrm{FMF},\end{array}$$

where AOD_f is the fine-particle AOD and FMF is the fine-mode fraction. V_f,column can be expressed by the vertical integral of particle volume size distributions (PVSDs) within a certain aerodynamic diameter range of

$$\begin{array}{}\text{(6)}& V_{\mathrm{f},\mathrm{column}}=\underset{\mathrm{0}}{\overset{D_{\mathrm{p},\mathrm{c}}}{\int }}V\left(D_{\mathrm{p}}\right)\mathrm{d}{D}_{\mathrm{p}},\end{array}$$

where D_p,c represents the cutting diameter, the empirical value of 2.0 µm is chosen based on previous literature (Hand and Kreidenweis, 2002; Hänel and Thudium, 1977), and V(D_p) represents the PVSD corresponding to the geometric equivalent diameter (D_p).

3.1.2 Specific process and limitations

The PMRS method is developed from Eq. (4). Based on satellite AOD, the near-surface PM_2.5 can be obtained through multistep transformation. Figure 2a shows its specific process. Each arrow refers to a step, respectively, size cutting (output: AOD_f), volume visualization (output: V_f,column), bottom isolation (output: V_f, fine-particle volume near the ground), particle drying (output: V_f,dry, dry V_f), and PM_2.5 weighting. The overall expression is as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\mathrm{PM}}_{\mathrm{2.5}}=\mathrm{AOD}{\displaystyle \frac{\mathrm{FMF}\cdot {\mathrm{VE}}_{\mathrm{f}}\cdot {\mathit{\rho }}_{\mathrm{f},\mathrm{dry}}}{\mathrm{PBLH}\cdot f_{\mathrm{0}}\left(\mathrm{RH}\right)}},\text{(8)}& {\displaystyle }{f}_{\mathrm{0}}\left(\mathrm{RH}\right)={\left(\mathrm{1}-{\displaystyle \frac{\mathrm{RH}}{\mathrm{100}}}\right)}^{-\mathrm{1}},\end{array}$$

where FMF denotes the fine-mode fraction, ρ_f,dry denotes the dry mass density of PM_2.5, and PBLH represents the planet boundary layer height. f₀(RH) represents the approximation of f(RH) in Eq. (4), as expressed in Eq. (8). Considering the aerosol types in different regions, PMRS fits VE_f to a quadratic polynomial relation of FMF (Zhang and Li, 2015):

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\mathrm{VE}}_{\mathrm{f}}& =\mathrm{0.2887}{\mathrm{FMF}}^{\mathrm{2}}-\mathrm{0.4663}\mathrm{FMF}+\mathrm{0.356}\\ & (\mathrm{0.1}\le \mathrm{FMF}\le \mathrm{1.0}).\end{array}\end{array}$$

Figure 2 Surface PM_2.5 estimation flow of RF-PMRS. (a) The five steps of the PMRS method. Gray boxes are the intermediate outputs, blue boxes are the input data, and orange boxes denote the variables to be optimized. (b) The specific optimization of RF-PMRS: FMF dataset replacement and VE_f simulation by the RF model.

PMRS has strong physical significance; the calculation steps are well-defined and site-independent. Zhang and Li (2015) tested the performance of PMRS on 15 stations, and the validation results had an uncertainty of 34 %. Compared with the ground value of the city of Jinhua in China, a 31.3 % relative error was generated in Li et al. (2016). Moreover, Zhang et al. (2020) applied it to the PM_2.5 change analysis and prediction experiments in China over 20 years. However, there may be a more complex nonlinear relationship between VE_f and FMF, not just a simple quadratic formula. Since VE_f is related to the aerosol type, adding other spatiotemporal variables may optimize the fitting process. Additionally, high-quality FMF data are the basic guarantee for the estimated PM_2.5 quality. In a word, to further improve the physical method, a better nonlinear model between VE_f and related variables from reliable datasets needs to be explored.

Figure 3 Specific steps for simulating VE_f based on ML in our RF-PMRS method. The map used in step 1 is from NASA Visible Earth (https://visibleearth.nasa.gov/images/57752/blue-marble-land-surface-shallow-water-and-shaded-topography, last access: 15 June 2023). The red points in step 1 represent the distribution of the nine AERONET sites, and the two yellow quadrangles in the zoom-in view highlight the Beijing (BJ) and Beijing-CAMS (BC) sites.

3.2 Optimization method: RF-PMRS

Therefore, to overcome the above disadvantages, an optimized method called RF-PMRS is proposed. Figure 2b shows the process of our method, while optimizations for FMF and VE_f are described separately below.

1. FMF dataset selection.

   We introduce the Phy-DL FMF dataset into the PMRS method to improve the accuracy of size-cutting results. In terms of performance, it exhibits higher accuracy and wider space–time coverage than satellite products (Yan, 2021). See the Data section for details.

2. VE_f simulation based on ML.

   The main idea is to establish an ML model between the VE_f truth obtained from multiple AERONET sites and related variables, thus improving the subsequent VE_f simulation accuracy (Fig. 3).

   • Step 1: VE_f calculation.

     The VE_f true values are calculated concerning Eqs. (5)–(6). Due to the spatiotemporal variability in different aerosol types, we calculate the VE_f values at nine AERONET stations around the world (Table 1) to train a universal model. The first step in Fig. 3 shows their distribution characteristics. Among them, Beijing and Beijing-CAMS sites are highlighted since they participate in the subsequent point validation experiment.

   • Step 2: VE_f-related variable selection.

     According to the theory, FMF is selected as the most important modeling variable. Previous studies have also shown that the FMF–VE_f relationship has a good single-value correspondence, which is not affected by AOD. Compared with AOD_f and V_f,column, FMF is a better indicator for estimation (Zhang and Li, 2015). In addition, considering the spatiotemporal heterogeneity of VE_f, the latitude (lat), longitude (long), and data time (month, day) of each site are added to the training.

   • Step 3: RF model establishment.

     From step 2, VE_f can be expressed as

     $$\begin{array}{}\text{(10)}& {\mathrm{VE}}_{\mathrm{f}}=f(\mathrm{FMF},\mathrm{lat},\mathrm{long},\mathrm{month},\mathrm{day}).\end{array}$$

     We optimize the VE_f expression based on the random forest (RF) algorithm. RF is made up of multiple decision trees that can build high-accuracy models based on fewer variables (Ho, 1995; Yang et al., 2020). This ensemble ML method randomly samples the training dataset to form multiple subsets, and random combinations of features are selected in node splitting (Belgiu and Drăguţ, 2016). The specific process is to (1) generate training subsets, (2) build an optimal model, and (3) calculate the result (Fig. 3 shows its flowchart). Note that the station FMF values (S−FMF) from AERONET sites are used when training.

   • Step 4: accuracy validation.

     The VE_f estimation is also based on Eq. (10), where f is the optimal relationship after RF parameter adjustment, and Phy-DL FMF is applied to realize the extension of model results from the point to the surface. A 10-fold cross-validation (CV) (Rodriguez et al., 2009) and isolated validation (IV) are used to evaluate model performance (for details of the validation methods, see Appendix A1).

3. PM_2.5 value estimation and evaluation

   Then, we calculate PM_2.5 according to the corresponding process (Eq. 7). The variables (in Sect. 2.2 to 2.4) are spatially matched to ground sites at their respective resolutions. Based on UTC, the PM_2.5 validation is conducted on a daily scale in 2017. Because of the effective quantity of the AERONET public dataset and MODIS data, we choose 2017 as the representative year. Note that we select the measured empirical value of ρ_f,dry (i.e., 1.5 g cm⁻³) for the NC region from Gao et al. (2007).

   The statistical indicators used in the evaluation include the correlation coefficient (R), mean bias (MB), relative mean bias (RMB), root mean square error (RMSE), and mean absolute error (MAE). In addition, the relative predictive error (RPE) is added to validate the accuracy of the RF-based VE_f model. See Appendix A2 for specific information on these indicators.

4 Experiment results

Three main experiments are conducted to verify the proposed RF-PMRS method, and the specific information is shown in Table 2.

Table 2 A brief information summary of the experiments conducted in our study.

Table 3 Performance statistics of the RF model for training VE_f. N represents the number of data, and VE_f has no unit.

4.1 RF model performance for training VE_f

The simulation model of VE_f is trained based on the data in Table 1. Specifically, the 10-fold CV result is used to determine the optimal combination of parameters for the model (see Appendix A3 for the adjustment of the model parameters). Considering that the completeness of the training data will optimize the generalization performance of the model, the experiment fine-tunes the model based on all the original datasets (the training period of Table 1) under the optimal parameters, and then the final RF model is constructed. This is also the most common method for ML model construction. Next, the IV experiment provides independent time validation of the final model.

Table 3 shows the CV and IV results to, respectively, demonstrate the internal and external accuracy of the final RF model. It can be seen that RF can capture the complex relationship between VE_f and related variables well. R is as high as 0.974 (0.975), RMSE and MAE are both small, and RPE is around 30 %, which suggests the desired estimation accuracy. Overall, the CV results represent the great performance of the RF model for extracting information, i.e., the relationship of multisource data to VE_f. In the meantime, the statistical results in the CV and IV experiments are similar, indicating that the RF model has no obvious overfitting phenomenon.

4.2 Accuracy evaluation of PMRS and RF-PMRS at AERONET stations

The purpose of RF-PMRS is to construct an optimal model from the obtained point matching data pairs and generalize it to the space–time continuous surface data for VE_f derivation. In the subsequent experiments in Sect. 4.2 and 4.3, the VE_f values are obtained by introducing the Phy-DL FMF dataset (surface data) to the final RF model. At the same time, the Phy-DL FMF data are also applied to the PM_2.5 calculation process (FMF variable in Eq. 7) for a wide range of PM_2.5 concentration.

Then, the experiment compares PM_2.5 results of PMRS and RF-PMRS at the Beijing (BJ) and Beijing-CAMS (BC) AERONET sites in 2017. Here, RF-PMRS simulates VE_f based on RF and replaces the polynomial of the PMRS method. Note that the results of the two sites are compared with their respective nearest ground PM_2.5 stations (distances of 3.64 and 3.91 km, respectively, in line with the representative range of ground stations in previous studies; Shi et al., 2018). Figure 4 displays the time series of PM_2.5 values for different models at two sites. The blue line fits the red line better than the gray one, confirming that the PM_2.5 results of RF-PMRS are closer to the true values. Within the range of the black circles at positions 1 and 2, the variation in RF-PMRS results has better consistency with the ground truth, while the PMRS results show dislocation and excessive growth. The overall performance of the RF-PMRS estimations can signify the effectiveness of our proposed method framework. As observed in the red boxes at positions 3 and 4, both models have a certain degree of deviation, which is found to be consistent with the time regularity of the AOD high values. Meanwhile, Fig. B1 (in Appendix B) plots the bias time series between PMRS and RF-PMRS and in situ values. As can be seen, the bias of the optimization method (RF-PMRS) is stably distributed around zero, which greatly reduces the numerical uncertainty. It is worth noting that our method has mitigated the apparent overestimation of the original model (PMRS) well in the case of above-normal aerosol loadings. Furthermore, the average PM_2.5 values from ground stations, PMRS, and RF-PMRS are compared. As for the two sites, the RF-PMRS results are satisfactory. As depicted in Fig. 5, the RF-PMRS and station mean values are close, with a difference of 4.82 (BJ) and 2.73 µg m⁻³ (BC), suggesting a good estimation. Nevertheless, the PMRS results have deviations greater than 40 µg m⁻³, and overestimation exists at both sites. It can be inferred that, in our proposed method, the optimization of VE_f can greatly improve the PM_2.5 estimation accuracy.

Figure 4 Three PM_2.5 time series at the Beijing (BJ) and Beijing-CAMS (BC) sites under their respective DOYs in 2017. Here, DOY (valid) means the day of the year with valid AOD, FMF, and other PM_2.5-related data. Gray, blue, and red lines represent PM_2.5 values of PMRS, RF-PMRS, and the stations (STA), respectively. The red boxes and black circles select a specific period for analysis.

Figure 5 Annual average PM_2.5 values from the stations (left), RF-PMRS (middle), and PMRS (right) at the BJ and BC sites.

To visually compare the optimization effect, Fig. 6 plots the PM_2.5 bias distribution patterns for two methods. From the boxplot, the average PM_2.5 bias of RF-PMRS is close to zero (less than 5 µg m⁻³), which is greatly lower than that of PMRS. Moreover, PMRS PM_2.5 has a larger deviation range, which manifests in two aspects. One is the maximum bias; specifically, it has exceeded 100 µg m⁻³ at the BC site. The other is the overall distribution of the data bias; the BJ site ones are mostly distributed below 0, indicating an obvious overestimation. As for RF-PMRS, the above circumstances are not obviously reflected in it. In addition, as can be seen from the indicators, the RMSE and MAE of RF-PMRS PM_2.5 decrease by about half in comparison with PMRS. The experiment has confirmed that the RF-PMRS PM_2.5 values have a strong linear relationship with the ground truth at both sites, with R around 0.8 (0.82 at BJ and 0.78 at BC). Such a large optimization effect is attributed to the VE_f expression replacement to the fitted RF model.

Figure 6 Boxplots of (a) RF-PMRS and (b) PMRS PM_2.5 bias at the BJ and BC sites. The upper (lower) black line of each box represents the largest (smallest) value, the upper (lower) blue border represents the upper (lower) quartile, and the red line denotes the median. The yellow, orange, and gray symbols are the MB, RMSE, and MAE of the corresponding PM_2.5 concentration.

4.3 Generalization performance of RF-PMRS

Then, we estimate PM_2.5 based on PMRS and RF-PMRS within North China in 2017 (Fig. 1 exhibits the distribution pattern of the validation stations). Table 4 shows the accuracy statistics. It can be seen that RF-PMRS greatly reduces the bias (about 44.87 %), with an MB of about 2.31 µg m⁻³. Similar to the results at the sites, the RF-PMRS method can derive PM_2.5 concentration with practically no overestimation (underestimation). Although there is not much difference in the R values of the two models (R of RF-PMRS is only improved by 0.01), RMSE and MAE decrease by about 39.96 and 18.86 µg m⁻³, respectively. As a result, the optimized method deserves to be considered excellent.

Table 4 Validation results of PMRS and RF-PMRS PM_2.5 in North China.

Meanwhile, the PM_2.5 scatterplots are presented below. As depicted in Fig. 7, there are sufficient estimated samples (28 305) in the NC region, which guarantees the credibility of our validation results. In general, the RF-PMRS PM_2.5 values are distributed around the 1:1 reference line evenly, with a slightly higher R of 0.70 compared to that of the original method. The slope of the linear fitting relationship reaches 0.82, which indicates that the proposed method greatly reduces the overestimation of PMRS with a linear slope of 1.46. Although the overall performance of the RF-PMRS estimations maintains an excellent level, defects do remain. To be specific, in areas with high PM_2.5 concentration (especially greater than 150 µg m⁻³), RF-PMRS results exist with a slight underestimation. It may be caused by the relatively small number of high-value PM_2.5 points (only 1319 out of 28 305), which makes it difficult to adequately reflect the fitting effect of the method.

Figure 7 Validation scatterplots of PM_2.5 results from (a) PMRS and (b) RF-PMRS. Dashed red lines are 1:1 reference lines, and solid blue lines stand for the linear fits. The right legends show the point densities (frequency) represented by different colors.

As for RF-PMRS, the deviation is reduced to a large extent, so the maps of the probability density function based on the bias of PMRS and RF-PMRS are further drawn. Figure 8 visualizes the probability densities within different bias ranges. In terms of distribution characteristics, the overall bias of RF-PMRS from the value of 0 (solid black line) is small. About the curve shape, it is high and narrow, manifesting in the fact that the bias has a lower standard deviation (SD) and is more prone to appear around the mean. However, PMRS shows a more discrete distribution pattern, and there are many outliers outside the range of greater than 600 µg m⁻³. Simultaneously, as can be concluded from the three boxes, within the bias range of ±20 and ±40 µg m⁻³, the data numbers of RF-PMRS results increase by 8.32 % and 12.81 %, respectively. Outside the range of ±100 µg m⁻³, the number decreases by 9.10 %. Therefore, as far as the accuracy is concerned, RF-PMRS results have lower bias and better stability.

In addition to the above general performance comparison in Sect. 4.3, Fig. 9 presents the annual average RMSE spatial distribution of PMRS and RF-PMRS PM_2.5 at NC stations. The two methods show a large deviation in the middle and southeast, and the RMSE map of PMRS has more red points. However, RF-PMRS can weaken this phenomenon very well since its RMSE representative colors are generally light. In particular, the proportion of dark-red sites (RMSE greater than 60 µg m⁻³) decreases from 65.44 % (PMRS) to 4.15 % (RF-PMRS). In the areas where the ground stations are clustered, the deviation also reduces significantly.

Figure 8 Probability density functions of PMRS (yellow) and RF-PMRS (green) PM_2.5 bias. The dotted red, blue, and gray lines indicate the bias boundaries of ±20, ±40, and ±100 µg m⁻³, respectively. µ and σ represent the mean value and standard deviation of each data point.

Figure 9 RMSE of the yearly average PM_2.5 concentration values between different models and ground stations (a PMRS PM_2.5, b RF-PMRS PM_2.5). Note that the top red color of the RMSE legend indicates RMSE values equal to or greater than 60 µg m⁻³.

In a word, the above analysis demonstrates that compared with the simple quadratic polynomial relationship (Eq. 9), the established RF model in RF-PMRS can more accurately capture the relationship between VE_f and multiple variables, thereby improving the PM_2.5 estimation accuracy.

5 Discussion

5.1 Accuracy comparison of PMRS using MODIS FMF and Phy-DL FMF

To confirm the superiority of the Phy-DL FMF data adopted in our method framework, the experiment takes the BJ and BC sites as examples (in 2017) and then compares the PM_2.5 accuracy and the number of effective days calculated by PMRS based on different FMF values. Table 5 presents the overall day-level results. Here, “DOY” means the day of the year and “valid” means that all variables related to the PM_2.5 calculation are valid. As can be seen, after the FMF replacement, the number of valid DOYs grows (an increase of 113 d), which illustrates that the number of effective PM_2.5 concentrations has gone up by about 5 times. Moreover, the accuracy has been significantly enhanced, with R having increased by about 0.30 and RMSE and MAE having decreased by 26.14 % and 16.47 %, accordingly. On the whole, Phy-DL FMF contributes to the improvement in PMRS results, signifying the first step in optimization of the proposed RF-PMRS method is effective.

Table 5 Validation results of the PMRS method using different FMF data. The valid DOY refers to the number of days for which the AOD, FMF, and other data are not missing when calculating PM_2.5. Note that since the valid days of the two schemes are different, the MB and RMB are not compared.

5.2 Performance compared with other ML models

Different machine learning models are suitable for diverse research data, and decision tree (DT) models can better fit experiments with fewer variables, such as this study. For comparison, except for RF, the extremely randomized tree (ERT) (Geurts et al., 2006) and gradient boosting decision tree (GBDT) (Friedman, 2001) models have also been established. The results of training VE_f based on the above three DT models are presented in Tables 6 and 7. By contrast, RF performs best in CV and IV experiments, as indicated by the multiple accuracy indicators. Although the ERT and GBDT models are comparable to RF in some indicators, there exists a certain degree of overfitting in the above two models, which is manifested in the fact that their IV results are clearly worse than their respective CV ones. Thus, the RF model is applied to our study.

Table 6 Cross-validation results for comparison of the decision tree models for training VE_f. N represents the number of data, and VE_f has no unit.

Table 7 Isolated-validation results in comparison of the decision tree models for training VE_f. The indicators are the same as those in Table 6.

5.3 Feature importance of the embedded RF model

Additionally, the feature importance of RF is calculated to evaluate the contribution of model predictors to VE_f simulation. Figure B2 (in Appendix B) shows the results by normalization (taking 100 as the total). Without a doubt, FMF accounts for the largest proportion, about 76.4 %, which is consistent with the analysis when selecting the VE_f-related variables (see Sect. 3.2). The contribution of spatiotemporal variables is about $\mathrm{1}/\mathrm{3}$ of FMF, which indirectly affirms the credibility of RF feature learning. Also, it provides a basis for further uncertainty optimization of VE_f and PM_2.5 accuracy.

5.4 Advantages and disadvantages

5.4.1 Advantages of the RF-PMRS method

From the perspective of model parameter optimization, this paper embeds RF to replace the subprocess parameter of the semi-empirical physical model. As a result, the proposed method, RF-PMRS, reduces the uncertainty in the complex physical parameter (i.e., VE_f) based on the estimation steps of strong physical significance and realizes the coupling of machine learning and the model mechanism. The proposed method does not rely on the PM_2.5 values of ground stations and is not affected by the station density and distribution mode, which can estimate the PM_2.5 concentration independently.

Meanwhile, as for the method, we construct the VE_f model based on RF using high-precision point data and extend it to surface data for PM_2.5 estimations. The experimental results demonstrate the overall performance of the model (Sect. 4.1) and its applicability in North China (Sect. 4.2 to 4.3), showing that the method has certain universality from the point scale to the surface scale.

1. The overall performance of the model is high. We use the ground data of nine AERONET sites around the world to train the RF model and simulate the VE_f values; the site distribution is relatively uniform, and the number of training data is sufficient. Table 1 shows a total of 6463 data matching pairs in the training period, which is enough to establish a credible RF model. Table 3 results show that in IV experiments, the accuracy of the model is good and can be generalized in different periods. For VE_f, the model shows both high internal accuracy (CV) and external accuracy (IV), so it can be generalized in regions with different aerosol types.

2. In the subsequent PM_2.5 estimation, the model displays high applicability in North China. From the perspective of model construction, the four aerosol types are the classification basis of the training data, and comprehensive modeling can improve the generalization performance. Also, the addition of spatiotemporal variables can increase the model applicability in North China. On the other hand, the number of stations used in an area does not determine the regional accuracy of the established model, which can be derived from our results. Compared with the PM_2.5 ground measurements in the NC region, the relative deviation of the RF-PMRS PM_2.5 is only 2.31 µg m⁻³, which confirms that RF can represent the relationships within North China.

5.4.2 Limitations on the scope of the validation region

However, there are still some shortcomings, mainly manifested in the scope of the validation region. Due to limited experimental data, we only conduct experiments in North China (the main aerosol type is urban–industrial). The main reasons are as follows.

1. Insufficient ρ_f,dry value.

   As the empirical value in the semi-physical empirical model, the ρ_f,dry value is often obtained by field measurements and induction. The insufficient ρ_f,dry values hinder the derivation of PM_2.5 in other regions, and more research results are needed.

2. Disclosure limits on global PM_2.5 ground measurements.

   Accurate and sufficient in situ PM_2.5 values allow for the verification of estimated PM_2.5 results.

3. Fewer public AERONET sites.

   Therefore, only the BJ and BC sites in North China are used for representative point-scale validation.

5.4.3 Data differences and uncertainty analysis

In the RF-PMRS method, the VE_f model constructed by high-precision site data is generalized to surface data for validation, and the data types involved are as follows.

1. AERONET AOD vs. MODIS AOD.

   Two types of AOD are used for different experimental steps, among which AERONET AOD is applied to calculate the true values of VE_f for establishing the RF simulation model. The RF model construction is a step of PM_2.5 estimation (as the VE_f variable in Eq. 7). MODIS AOD is satellite AOD data, the most commonly used remote sensing data for large-scale retrieval of PM_2.5. It is an important variable for PM_2.5 estimation in RF-PMRS (as the AOD variable in Eq. 7). Thus, there is no error in the PM_2.5 calculation caused by AOD category replacement.

   As for uncertainty, AERONET AOD provides truth values for calculating VE_f, which theoretically has negligible uncertainty, and the simulation accuracy of VE_f represents its influence on estimating PM_2.5 to a certain extent. It is generally considered that MODIS AOD has guaranteed quality and sufficient accuracy to be used directly.

2. S-FMF vs. Phy-DL FMF.

   S-FMF is obtained directly from the AERONET monitoring sites and is one of the variables of the RF model (as the FMF variable in Eq. 10). In the point-to-surface extension, Phy-DL FMF is introduced into the RF model to replace S-FMF, and the 2017 VE_f values are obtained. The basis of the above replacement is that the accuracy of Phy-DL FMF is relatively consistent with that of S-FMF (Yan et al., 2022). Moreover, Phy-DL FMF data are applied to the PM_2.5 estimation steps (as the FMF variable in Eq. 7) for a wider range of validation experiments. The results show that the PM_2.5 concentration estimated by RF-PMRS has high accuracy, proving the credibility of Phy-DL FMF.

3. FMF uncertainty.

   Different surface data sources may affect the PM_2.5 results, introducing some uncertainty. Section 5.1 compares the PM_2.5 accuracy using two FMF data in 2017. The data missing time for MODIS FMF and Phy-DL FMF in North China are different, which can be found in the statistics on their respective available days (referred to as valid DOYs). There are far more valid days based on Phy-DL FMF than MODIS FMF (143 and 31 d), demonstrating the superiority of Phy-DL FMF. Although the specific validation time of two FMF varies, the overall accuracy of the PM_2.5 estimation (which can be regarded as the average accuracy over the year) shows that the Phy-DL FMF increases R to 0.68 (MODIS FMF: 0.38) with low uncertainty.

4. ρ_f,dry uncertainty.

   As introduced earlier, the ρ_f,dry value is often obtained by field measurements. In our study, we select 1.5 g cm⁻³ as the ρ_f,dry value for North China. There are certain variations in the empirical values of different regions, and there will be errors (uncertainty) between the values in Beijing and other places in the NC region. However, our experimental area is not large, and we use 1.5 g cm⁻³ to represent ρ_f,dry of the whole region, which has been applied in previous articles (Zhang and Li, 2015; Li et al., 2016).

5. Uncertainty between variable resolutions.

   In most experiments, the lowest resolution of all data will be taken as the unified resolution when obtaining data values. The different data may lose some spatial details during the upsampling–downsampling process, which brings uncertainty to the estimation results. In the RF-PMRS method, there is no such uncertainty problem. We set 1^∘ as the unified spatial unit and take the longitude and latitude of each cell's center as the reference longitude and latitude. The variables in the Data section are spatially matched to ground sites at their respective resolutions, and the space–time matching method has been described in the Methods section. So, all kinds of data uncertainties only exist in their instrument measurement or statistical release.

   Overall, RF-PMRS shows excellent estimation performance in North China, and the accuracy of surface PM_2.5 estimation based on remote sensing data is guaranteed. Next, with the improvement in related experimental data, we will verify our proposed method in a broader range and continuously optimize it from all aspects.

6 Conclusions

Among various satellite remote sensing methods for PM_2.5 retrieval, the semi-empirical physical approach has strong physical significance and clear calculation steps and derives the PM_2.5 mass concentration independently of in situ observations. However, the parameters of optical properties are difficult to express, requiring them to be optimized. Hence, the study proposes a method (RF-PMRS) that embeds machine learning in a physical model to obtain surface PM_2.5: (1) based on the PMRS method, select the Phy-DL FMF product with a combined mechanism, and (2) use the RF model to fit the parameter VE_f, rather than a simple quadratic polynomial. In the point-to-surface validation, RF-PMRS shows great optimized performance. Experiments at two AERONET sites show that R reaches up to 0.8. In North China, RMSE decreases by 39.95 µg m⁻³ with a 44.87 % reduction in relative deviation. In the future, we will further explore the combination of an atmospheric mechanism and machine learning and then research the PM_2.5 retrieval methods with physical meaning and higher accuracy.

Appendix A: Supplementary description

A1 The 10-fold cross-validation and isolated validation

The sample-based 10-fold cross-validation method is applied to tune the model parameters and test the internal accuracy of our model. The original dataset is randomly divided into 10 parts, 9 of which are used as the training set for model fitting, with the remaining 1 used for prediction; then the cross-validation process is repeated for 10 rounds until each data point has been used as the test set.

At the same time, when verifying the RF-based VE_f model, the dataset in the period that did not participate in the training in Table 1 is used for isolated validation.

A2 Statistical indicators

$$\begin{array}{c}{\displaystyle }R={\displaystyle \frac{\sum _{i=\mathrm{1}}^m\left(y_i-\stackrel{\mathrm{‾}}{y}\right)\sum _{i=\mathrm{1}}^m\left(f_i-\stackrel{\mathrm{‾}}{f}\right)}{\sqrt{\sum _{i=\mathrm{1}}^m{\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}^{\mathrm{2}}}\sqrt{\sum _{i=\mathrm{1}}^m{\left(f_i-\stackrel{\mathrm{‾}}{f}\right)}^{\mathrm{2}}}}},\\ {\displaystyle }\mathrm{MB}=\stackrel{\mathrm{‾}}{y}-\stackrel{\mathrm{‾}}{f},\\ {\displaystyle }\mathrm{RMB}=\mathrm{abs}\left({\displaystyle \frac{\stackrel{\mathrm{‾}}{y}-\stackrel{\mathrm{‾}}{f}}{\stackrel{\mathrm{‾}}{y}}}\right),\\ {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{m}}\sum _{i=\mathrm{1}}^m{\left(y_i-f_i\right)}^{\mathrm{2}}},\\ {\displaystyle }\mathrm{MAE}={\displaystyle \frac{\mathrm{1}}{m}}\sum _{i=\mathrm{1}}^m\left|y_i-f_i\right|,\\ {\displaystyle }\mathrm{RPE}={\displaystyle \frac{\sqrt{\frac{\mathrm{1}}{m}\sum _{i=\mathrm{1}}^m{\left(y_i-f_i\right)}^{\mathrm{2}}}}{\stackrel{\mathrm{‾}}{y}}},\end{array}$$

where m is the total number of observations; i is the number of measurements; y_i is the ith observation; f_i is the corresponding estimation result; and $\stackrel{\mathrm{‾}}{y}$ and $\stackrel{\mathrm{‾}}{f}$ are the averages of all observations and estimates, respectively.

A3 Parameter adjustments of the RF model

The four parameters of RF are adjusted; i.e., the correlation coefficient r changes with the (a) number of trees, (b) maximum depth, (c) maximum number of features when splitting, and (d) minimum number of split samples. Experiments show that the maximum depth varies greatly in a small range. To prevent overfitting, the four parameters of RF are adjusted to 60, 10, 2, and 8. It can ensure high accuracy while improving training efficiency.

Appendix B: Figures

Figure B1 The time series of PMRS and RF-PMRS PM_2.5 bias at the Beijing and Beijing-CAMS sites under their respective DOYs in 2017. The orange line represents the bias between the PM_2.5 values of PMRS and the stations, while the blue one indicates the PM_2.5 difference between RF-PMRS and the stations.

Figure B2 The predictor importance results (normalized) of the RF model for training VE_f.