A National-Scale 1-km Resolution PM2.5 Estimation Modelover Japan Using MAIAC AOD and a Two-Stage Random Forest Model

Satellite-based models for estimating concentrations of particulate matter with an aerodynamic diameter less than 2.5 μm (PM2.5) have seldom been developed in islands with complex topography over the monsoon area, where the transport of PM2.5 is influenced by both the synoptic-scale winds and local-scale circulations compared with the continental regions. We validated Multi-Angle Implementation of Atmospheric Correction (MAIAC) aerosol optical depth (AOD) with ground observations in Japan and developed a 1-km-resolution national-scale model between 2011 and 2016 to estimate daily PM2.5 concentrations. A two-stage random forest model integrating MAIAC AOD with meteorological variables and land use data was applied to develop the model. The first-stage random forest model was used to impute the missing AOD values. The second-stage random forest model was then utilised to estimate ground PM2.5 concentrations. Ten-fold cross-validation was performed to evaluate the model performance. There was good consistency between MAIAC AOD and ground truth in Japan (correlation coefficient = 0.82 and 74.62% of data falling within the expected error). For model training, the model showed a training coefficient of determination (R2) of 0.98 and a root mean square error (RMSE) of 1.22 μg/m3. For the 10-fold cross-validation, the cross-validation R2 and RMSE of the model were 0.86 and 3.02 μg/m3, respectively. A subsite validation was used to validate the model at the grids overlapping with the AERONET sites, and the model performance was excellent at these sites with a validation R2 (RMSE) of 0.94 (1.78 μg/m3). Additionally, the model performance increased as increased AOD coverage. The top-ten important predictors for estimating ground PM2.5 concentrations were day of the year, temperature, AOD, relative humidity, 10-m-height zonal wind, 10-m-height meridional wind, boundary layer height, precipitation, surface pressure, and population density. MAIAC AOD showed high retrieval accuracy in Japan. The performance of the satellite-based model was excellent, which showed that PM2.5 estimates derived from the model were reliable and accurate. These estimates can be used to assess both the short-term and long-term effects of PM2.5 on health outcomes in epidemiological studies.


Introduction
Particulate matter (PM) is a mixture of solid particles and liquid droplets suspended in the air that can be classified into coarse (PM with an aerodynamic diameter less than 10 µm; PM 10 ) and fine (PM with an aerodynamic diameter less than 2.5 µm; PM 2.5 ) fractions based on their size [1,2]. PM 2.5 is considered more toxic than PM 10 because of its smaller diameter, meaning it can deposit into the alveolar region of the lungs and induce a series of immune responses [2]. Many epidemiological studies reported that exposure to PM 2.5 is associated with adverse health outcomes, such as allergic respiratory diseases, asthma, cardiovascular within 24° to 46° N and 123° to 149° E (Figure 1). The remote offshore islands, i.e. Bonin Islands, were excluded from the study. To facilitate data integration, the study was divided into 379,643 grids with a resolution of 45 s in longitude and 30 s in lati based on the Japan third mesh data (approximately 1-km × 1-km) [32].
The study collected satellite-based AOD, meteorological data, and land use data data sources, coverage years, and spatial and temporal resolutions are summarised in ble S1). The details of these data are described in the following sections.

MAIAC AOD
The 1-km-resolution MAIAC algorithm-based level-gridded AOD data (MCD1 collection 6 product) between 1 January 2011 and 31 December 2016 were downloa from the Level 1 and Atmospheric Archive and Distribution System Distributed A Archive Center (LAADS DAAC; https://ladsweb.modaps.eosdis.nasa.gov/search/ cessed on 12 September 2021). MAIAC AOD at 550 nm (SDS name: Optical_Depth_ from Terra (overpass at around 10:30 a.m.) and Aqua (overpass at around 1:30 p.m.) the best quality (QA.CloudMask = Clear (code: 001) and QA.AdjacencyMask = C (code: 000)) were retrieved separately. A simple linear regression between MAIAC A from Terra and Aqua was developed by using complete pairs of AOD from these satellites. If the MAIAC AOD from Terra was missing for a specific day, the MAIAC A from Aqua was used to impute the missing value by simple linear regression, and versa [23]. Then, the arithmetic means of AOD from Terra and Aqua were calculate represent daily values.

AERONET AOD
The AErosol RObotic NETwork (AERONET) programme is a global ground-b remote sensing aerosol network established by NASA and PHOTONS (PHOtom pour le Traitement Opérationnel de Normalisation Satellitaire; Univ. of Lille 1, CNES CNRS-INSU) with several collaborators, including national agencies, institutes, and versities across the world [33]. The AOD from the AERONET can be used as 'gro truth' to validate satellite-based measurements. Total optical depths based on the A The study collected satellite-based AOD, meteorological data, and land use data (the data sources, coverage years, and spatial and temporal resolutions are summarised in Table S1). The details of these data are described in the following sections.

MAIAC AOD
The 1-km-resolution MAIAC algorithm-based level-gridded AOD data (MCD19A2 collection 6 product) between 1 January 2011 and 31 December 2016 were downloaded from the Level 1 and Atmospheric Archive and Distribution System Distributed Active Archive Center (LAADS DAAC; https://ladsweb.modaps.eosdis.nasa.gov/search/, accessed on 12 September 2021). MAIAC AOD at 550 nm (SDS name: Optical_Depth_055) from Terra (overpass at around 10:30 a.m.) and Aqua (overpass at around 1:30 p.m.) with the best quality (QA.CloudMask = Clear (code: 001) and QA.AdjacencyMask = Clear (code: 000)) were retrieved separately. A simple linear regression between MAIAC AOD from Terra and Aqua was developed by using complete pairs of AOD from these two satellites. If the MAIAC AOD from Terra was missing for a specific day, the MAIAC AOD from Aqua was used to impute the missing value by simple linear regression, and vice versa [23]. Then, the arithmetic means of AOD from Terra and Aqua were calculated to represent daily values.

AERONET AOD
The AErosol RObotic NETwork (AERONET) programme is a global ground-based remote sensing aerosol network established by NASA and PHOTONS (PHOtométrie pour le Traitement Opérationnel de Normalisation Satellitaire; Univ. of Lille 1, CNES and CNRS-INSU) with several collaborators, including national agencies, institutes, and universities across the world [33]. The AOD from the AERONET can be used as 'ground truth' to validate satellite-based measurements. Total optical depths based on the AOD levels (Level 2.0: cloud screened and quality assured) of 12 stations in Japan ( Figure S1) were downloaded from the Goddard Space Flight Center (https://aeronet.gsfc.nasa.gov/new_ Remote Sens. 2021, 13, 3657 4 of 18 web/index.html, accessed on 12 September 2021). The AERONET did not provide the AOD at 550 nm and therefore this was interpolated from the AOD at 500 and 675 nm [34].

Ground PM Measurements
Hourly PM 2.5 measurements from 2011 to 2016 were downloaded from the environmental numerical database of the National Institute for Environmental Studies [35]. There were 1,068 stations measuring PM 2.5 that were maintained by local government in 2016 (Figure 1). The monitoring stations used a β-ray attenuation method, tapered element oscillating microbalance, and a light-scattering method to measure PM 2.5 concentrations continuously and hourly. The air quality data were collected, summarised, and released by the National Institute for Environmental Studies [35], and are available online (https://www.nies.go.jp/igreen/index.html, accessed on 12 September 2021).

Meteorological Variables
The in situ measurements of daily temperature, relative humidity, precipitation, and surface pressure recorded by the automated meteorological data acquisition system (AMeDAS) during 2011-2016 were downloaded from the Japan Meteorological Agency. There are around 1300 rain gauges across Japan, and approximately 921 and 147 monitoring stations observing daily temperature and relative humidity across Japan, respectively [36]. Missing values of meteorological variables were excluded. The regression-kriging method was applied to estimate daily temperature and surface pressure in a 1-km-resolution grid lacking weather-monitoring stations [37]. Additionally, ordinary kriging was utilised to estimate daily values of relative humidity and precipitation in 1-km grids.
Boundary layer height, 10-m-height zonal wind (u10), and 10-m-height meridional wind (v10) were obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) Re-analysis Interim (ERA-Int). Data at 00.00 and 12.00 UTC, and at four steps (i.e., 3, 6, 9, and 12 h) after the observation time points, with a spatial resolution of 0.125 • longitude × 0.125 • latitude, were downloaded from the ECMWF online depository (https://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/, accessed on 12 September 2021). The inverse distance weighted (IDW) method was used to interpolate the boundary layer height, u10, and v10 data to the centroids of grids. The fixed search radius was set to 15 km for each centroid of fixed grids, and all the observed values of boundary layer height within the 15 km searching radius of a specific grid were averaged. The same IDW method was applied in our previous study [23].
Five-kilometre cloud-fraction data were accessed from the Terra and Aqua MODIS Collection 6 Level cloud product (MOD06_L2 and MYD06_L2), which consisted of cloud optical and physical parameters downloaded from the LAADS DAAC. The cloud-fraction data from Terra and Aqua were averaged as daily values. The IDW method with a fixed search radius of 7.5 km was used to interpolate the cloud fraction values into each grid.

Land Use Data
Urban and built-up area data with 50 m resolution during 2006-2011 (ID number of the category = #2) were retrieved from the High-Resolution Land Use and Land Cover (HRLULC) map published by the Japan Aerospace Exploration Agency (JAXA) Earth Observation Research Center (https://www.eorc.jaxa.jp/ALOS/en/lulc/lulc_index.htm, accessed on 12 September 2021). Industrial area data created in 2009 were accessed from the Ministry of Land, Infrastructure, Transport and Tourism (https://nlftp.mlit.go.jp/ ksj/gml/datalist/KsjTmplt-L05.html, accessed on 12 September 2021). Road network data were retrieved from the Global Map Japan v2.2 Vector data released in 2016 by the Geoinformation Authority of Japan (https://www.gsi.go.jp/kankyochiri/gm_japan_e. html, accessed on 12 September 2021). The distances from the centroid of the grid to the nearest primary road and highway were calculated, as well as the total length of roads within each grid. The normalised difference vegetation index (NDVI), a measurement of plant health, was extracted from the Terra MODIS MOD13Q1 product with a 16-day interval and 250 m spatial resolution in the sinusoidal coordinate system.

Population Counts
Annual population count data with 100 m resolution during 2011-2016 were downloaded from the WorldPop website (https://www.worldpop.org/geodata/listing?id=29, accessed on 12 September 2021). The annual average population count in each 1-km grid was calculated.

Elevation
Elevation data with 30 m spatial resolution was accessed from the Advance Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model v002 (ASTGTM.002) announced in 2011. The average elevation within each 1-km grid was calculated.

Statistical Methods
The Pearson's correlation coefficient (R), RMSE, mean bias error (MBE), percentage of data within the expected error (EE; ±(0.05 + 0.15AERONET AOD)), and slope of the linear regression were selected to evaluate the consistency between MAIAC AOD and AERONET AOD [38,39].
Random forest is an ensemble learning model, which has outperformed most traditional machine learning algorithms [40]. The random forest model is an extended machine learning algorithm of bagging (bootstrap aggregation) that not only builds a number of decision trees on bootstrapped training samples, but also randomly selects a subset of predictors (m try ) from the full set of predictors when building these trees [40,41]. For regression, the random forest model aggregates predictions from trees by averaging the predicted values from trees. In the present study, a two-stage random forest model was applied to develop the satellite-based PM 2.5 model. In the first stage, random forest was used to impute the missing values of MAIAC AOD on a weekly basis. Since satellite-based AOD values are frequently missing due to the presence of clouds and in the mountainous areas, longitude, latitude, temperature, relative humidity, cloud fraction, and elevation were used to impute missing AOD values based on the previous study and our early work [19,42]. The formula of the first-stage random forest model is as follows: where AOD ij is the AOD value in grid i at day j; Long i and Lat i are longitude and latitude in grid i; Temp ij is temperature in grid i at day j; Hum ij is relative humidity in grid i at day j; CF ij is the cloud fraction in grid i at day j; and Elevation i is elevation in grid i. The numbers of trees and m try were set to 1000 and 2 (one-third of the number of predictors), respectively. The weekly models were evaluated by calculating the out-of-bag (OOB) coefficient of determination (R 2 ) and root mean square error (RMSE). The second-stage random forest model was utilised to estimate ground PM 2.5 concentrations. The predictors in the second-stage model were selected based on the previous studies and our early works [19,23,42,43]. All the predictors have plausible physical mechanisms on PM 2.5 concentrations. Temperature, relative humidity, boundary layer height, surface pressure, wind speed and direction, and precipitation significantly influence the distribution and formation of PM 2.5 [44]. Industrial areas, road networks, and human population are sources of PM 2.5 . Accessibility of data is also a consideration for including predictors. The other gaseous pollutants, such as nitrogen oxides and sulphur dioxides, are precursors of PM 2.5 , and we did not include them in the model since they are not available in grids without air quality monitoring stations. The formula of the second-stage random forest is as follows: where PM 2.5 ij is PM 2.5 value in grid i at day j; AOD ij is AOD value in grid i at day j; Days j is the day of the year; Temp ij is temperature in grid i at day j; Hum ij is relative humidity in grid i at day j; BLH ij is boundary layer height in grid i at day j; SP ij is surface pressure in grid i at day j; U10 ij is zonal wind at a height of 10 m in grid i at day j; V10 ij is meridional wind at a height of 10 m in grid i at day j; Prec ij is precipitation in grid i at day j; NDVI i is the normalised difference vegetation index in grid i; Indus_area i is the area of industrial facilities in grid i; Urban_area i is the urban and built-up area in grid i; Elevation i is elevation in grid i; Road_len i is the length of roads in grid i; Distoprim i is the distance from the grid centroid to the nearest primary road in grid i; Distohigh i is the distance from the grid centroid to the nearest highway in grid i; and Pop i is the average population count in grid i. The descriptive statistics of predictors from 2011 to 2016 are presented in Table S2. Similar to the first-stage model, the number of trees was set to 1000, and m try was tuned from 1 to 17. m try was finally set to 15 for the second-stage model according to the highest OOB R 2 . Ten-fold cross-validation was performed to evaluate the performance of the second-stage model. During the cross-validation process, the dataset was randomly divided into ten nonoverlapping and equally sized partitions, and nine partitions and one partition were used as training and validating data, respectively. The cross-validation process was repeated ten times until all partitions had been used to validate the model performance. The crossvalidation results were represented by cross-validation R 2 , RMSE, percentage relative error (RE; (mean(Obsearvations) − mean(Predictions))/mean(Observations) × 100%), and mean absolute error (MAE). A tree-based ensemble model-eXtreme Gradient Boosting (XG-Boost) [45] was performed to estimate ground PM 2.5 concentrations and compared with the performance of random forest. The hyperparameters tuned to optimise the performance of XGBoost is shown in Table S3.
In addition to the ten-fold cross-validation, a subsite validation that used the grids overlapping with the AERONET sites as the validating data was conducted to check the model performance. Additionally, the ten-fold cross-validation results were stratified by seasons (spring: March-May; summer: June-August, autumn: September-November; winter: December-February) and quartiles of AOD coverage (below 25th percentile: <16.61%; 25th-50th percentile: 16.61-23.79%; 50th-75th percentile: 23.79-33.26%; and above 75th percentile: ≥33.26%) to evaluate the model performance.
The relative importance of predictors in random forest model was determined by calculating permutation importance scores. Additionally, accumulated local effect (ALE) plots were used to explore the effects of predictors on predicted PM 2.5 values.

Validation of MAIAC AOD with Ground Measurements over Japan
The overall mean coverage of MAIAC AOD in Japan was 15.48%. The highest coverage was in autumn (19.48%), followed by spring (17.94%), winter (14.05%), and summer (10.51%). The spatial distribution of MAIAC AOD coverage is shown in Figure 2. The mean coverage of AOD was higher in the Kanto, southern Chubu, and Kyushu regions; however, the mean coverage was extremely low (<4.79%) in the Hokkaido, Tohoku, and Okinawa regions ( Figure 2). S4. Overall, the 1-km-resolution MAIAC AOD was highly correlated with AERONET AOD (R = 0.82), 74.62% of data falling within the EE (Figure 3). The MAIAC AOD slightly underestimated AOD values in Japan (MBE = −0.006) ( Figure 3). The performance of MA-IAC AOD was highest in the Hokkaido University (R = 0.93, RMSE = 0.12, MBE = 0.033, and 65.15% of data within the EE), while the performance of MAIAC was lowest in the Osaka-North station (R = 0.61, RMSE = 0.108, MBE = 0.036, and 58.62% of samples within the EE) ( Figure 4).   The descriptive statistics for MAIAC AOD and AERONET AOD are shown in Table S4. Overall, the 1-km-resolution MAIAC AOD was highly correlated with AERONET AOD (R = 0.82), 74.62% of data falling within the EE ( The descriptive statistics for MAIAC AOD and AERONET AOD are shown in Table  S4. Overall, the 1-km-resolution MAIAC AOD was highly correlated with AERONET AOD (R = 0.82), 74.62% of data falling within the EE (Figure 3). The MAIAC AOD slightly underestimated AOD values in Japan (MBE = −0.006) (Figure 3). The performance of MA-IAC AOD was highest in the Hokkaido University (R = 0.93, RMSE = 0.12, MBE = 0.033, and 65.15% of data within the EE), while the performance of MAIAC was lowest in the Osaka-North station (R = 0.61, RMSE = 0.108, MBE = 0.036, and 58.62% of samples within the EE) ( Figure 4).

Model Performance of the First-Stage Model
The average OOB R 2 and RMSE for the first-stage random forest model were 0.97 (range, 0.93-1.00) and 0.016 (unitless; range, 0.009-0.038), respectively. After imputing missing values, the mean MAIAC AOD was 0.182 ± 0.145 (mean ± SD) overall, and was highest Remote Sens. 2021, 13, 3657 9 of 18 in spring (0.264 ± 0.182), followed by summer (0.206 ± 0.158), winter (0.134 ± 0.087), and autumn (0.128 ± 0.079) ( Table 2). The seasonal trend of MAIAC AOD was consistent with in situ PM 2.5 measurements.   Figure S2. The performance of random forest was better than XGBoost in this study. contrast, increased relative humidity, precipitation, and elevation were negatively correlated with PM2.5 concentrations ( Figure S8). Day of the year, temperature, v10 wind, u10 wind, boundary layer height, NDVI, distance to the nearest highway, distance to the nearest primary road, and urban and built-up area showed non-linear relationships with ground PM2.5 concentrations ( Figure S8).

PM2.5 Estimates
The overall average PM2.5 estimate was 12.5 ± 5.5 μg/m 3 . The PM2.5 estimates were higher in southern Japan and in several metropolitan areas, including Tokyo, Nagoya, and Osaka, than in the northern region ( Figure 6). The seasonal mean PM2.5 estimates were highest in spring (14.7 ± 6.0 μg/m 3 ), followed by summer (13.2 ± 5.5 μg/m 3 ), autumn (11.1 ± 4.6 μg/m 3 ), and winter (11.1 ± 4.8 μg/m 3 ) ( Table 3). The seasonal trends were consistent with those determined from in situ PM2.5 measurements. The seasonal spatial distribution of estimated PM2.5 is presented in Figure S9. The PM2.5 concentrations showed a decreased trend from 2011 to 2016 (Figure 7), which is consistent with the estimation in a previous study [46]. The spatial distributions of the residuals (differences between estimated PM 2.5 and in situ measurements) are shown in Figure S2. Generally, the random forest model tended to slightly overestimate PM 2.5 in the grids close to the shore (mean, 0.1 ± 0.5 µg/m 3 ; range, −1.8 to 1.8 µg/m 3 ). Nevertheless, the residuals on most grids (75%) were below 4% of the mean measurement of PM 2.5 ( Figure S2). Additionally, the observed and predicted PM 2.5 showed a highly consistent trend during 2011-2016 according to the time series plot ( Figure S3).
We used the grids overlapping with the AERONET sites to validate the second-stage random forest model, and the model performance was excellent at these sites with a validation R 2 (RMSE, MAE, and RE) of 0.94 (1.78 µg/m 3 , 1.40 µg/m 3 , and 1.55, respectively) ( Figure S5). Additionally, the training and cross-validation results stratified by season are shown in Figure S6. The model performances were similar across the four seasons, while the model performance in autumn was slightly worse than the other seasons (cross-validation R 2 , RMSE, MAE, and RE of 0.83, 2.83 µg/m 3 , 2.03 µg/m 3 , and −0.63, respectively) ( Figure S6). The model performance under different AOD coverages is presented in Figure S7. The model performance increased as AOD coverage increased, and the model performance at the grids with coverage rate ≥33.26% (cross-validation R 2 , RMSE, MAE, and RE of 0.87, 2.74 µg/m 3 , 1.93 µg/m 3 , and −0.45, respectively) was slightly higher than the first, second, and third quartiles ( Figure S7).
Increased AOD, surface pressure, population density, road length, and area of industrial facilities were positively associated with ground PM 2.5 concentrations ( Figure S8). By contrast, increased relative humidity, precipitation, and elevation were negatively correlated with PM 2.5 concentrations ( Figure S8). Day of the year, temperature, v10 wind, u10 wind, boundary layer height, NDVI, distance to the nearest highway, distance to the nearest primary road, and urban and built-up area showed non-linear relationships with ground PM 2.5 concentrations ( Figure S8).

PM 2.5 Estimates
The overall average PM 2.5 estimate was 12.5 ± 5.5 µg/m 3 . The PM 2.5 estimates were higher in southern Japan and in several metropolitan areas, including Tokyo, Nagoya, and Osaka, than in the northern region ( Figure 6). The seasonal mean PM 2.5 estimates were highest in spring (14.7 ± 6.0 µg/m 3 ), followed by summer (13.2 ± 5.5 µg/m 3 ), autumn (11.1 ± 4.6 µg/m 3 ), and winter (11.1 ± 4.8 µg/m 3 ) ( Table 3). The seasonal trends were consistent with those determined from in situ PM 2.5 measurements. The seasonal spatial distribution of estimated PM 2.5 is presented in Figure S9. The PM 2.5 concentrations showed a decreased trend from 2011 to 2016 (Figure 7), which is consistent with the estimation in a previous study [46]. land use data to estimate monthly or annual mean PM2.5 concentrations. This shows that satellite-based AOD can provide additional spatial and temporal information about daily PM2.5 concentrations and improve the model performance. Nearly all previous studies used ensemble models (i.e., random forest, gradient boosting, or XGBoost) to develop satellite-based PM2.5 models, and only a few used a deep learning model or support vector machine (as summarised in Table S5) [51]. A recent study in Italy built models between 2013 and 2015 using an ensemble generalised additive model, and combined the results derived from a linear mixed-effect model, random forest, and XGBoost. They found that the cross-validation R 2 (RMSE) was 0.79 (6.56 μg/m 3 ), 0.79 (5.29 μg/m 3 ), and 0.81 (6.34 μg/m 3 ) for 2013, 2014, and 2015, respectively [17]. Overall, our model performance is comparable to or better than those in these studies.
Our results showed that day of the year was the most important variable in the random forest model for estimating PM2.5 (relative importance of 19.92%). Early studies that developed satellite-based AOD-PM2.5 estimation models by using the linear mixed effect model needed to include day-specific intercepts and slopes to achieve high model performance because the relationship between AOD and PM2.5 may change day-to-day due to differences in mixing height, relative humidity, particle composition, and vertical profiles [13]. Recent studies that applied the machine learning algorithm also found that day of the year (or Julian day) is an important variable for modelling PM2.5 [22,52,53]. According to the ALE plot ( Figure S8), the relationship between day of the year and PM2.5 estimates displayed a bimodal pattern with two peaks between Julian day 0-100 (spring) and after Julian day 300 (winter). Day of the year can be a surrogate of seasonal effects. PM2.5 concentrations in Japan displayed an obvious seasonal variation. Both foreign and domestic sources may cause the seasonal variations of PM2.5 in Japan. Long-range transportation from Asia (yellow dust storm and anthropogenic aerosols) contributes to high PM2.5 levels in western Japan, particularly in spring and winter [9,54]. Photochemical smog events frequently occurred in summer can contribute to high concentrations of PM2.5 in summer [55]. Besides, increase in operation time of coal-burning power plants and heat appliance usage can cause raising PM2.5 concentrations during winter [55]. AOD was the third most important variable in our model (relative importance of 11.89%). There was a strong and positive relationship between AOD and PM2.5 according to the ALE plot ( Figure S8). However, the results from previous studies are contradictory. A study conducted on an island, Great Britain of the UK, found that satellite-based AOD 0.47 and 0.55 μm values retrieved from the MAIAC AOD product were not important

Discussion
In the present study, we validated MAIAC AOD with ground measurements of AOD and developed a 1-km-resolution satellite-based model for daily PM 2.5 concentrations during 2011-2016 by using a random forest model in Japan. The MAIAC AOD showed a high retrieval accuracy in Japan (R = 0.82, with the EE: 74.62%), which was comparable to those conducted in South Asia (including Pakistan, India, and Bangladesh) (R = 0.882, within the EE: 72.22% and R = 0.887, within the EE: 73.5 for Aqua and Terra MAIAC AOD, respectively) [47] and better than in Central Asia (R = 0.730, within the EE: 58.7% for spring; R = 0.709, within the EE: 66.7% for summer; R = 0.729, within the EE: 62.4% for autumn, and R = 0.744, with the EE: 67.9% for winter) [38].
The satellite-based model achieved excellent performance with a training and 10-fold cross-validation R 2 (RMSE, MAE, and RE) of 0.98 (1.22 µg/m 3 , 0.85 µg/m 3 , and −0.20) and 0.86 (3.02 µg/m 3 , 2.14 µg/m 3 , and −0.52), respectively. The top ten important predictors of the model were day of the year, temperature, MAIAC AOD, relative humidity, v10 wind, u10 wind, boundary layer height, precipitation, surface pressure, and population density. The overall average PM 2.5 estimate was 12.5 ± 5.5 µg/m 3 , and the estimates were higher in southern Japan and metropolitan areas. The estimates showed that PM 2.5 concentrations in Japan decreased from 2011 to 2016. Long-range transportation is a main contributor of high PM 2.5 in western Japan during winter and spring. China implemented an Air Pollution Prevention and Control Action Plan in 2013 with a series of stringent clean air actions during 2013-2017, which led to PM 2.5 rapidly decreasing across China [48]. In addition, the domestic emissions, especially those derived from road transport, decreased in Japan [46]. Hence, the decrease in both foreign and domestic emission sources may cause a long-term decreasing trend in Japan.
Several studies developed nationwide PM 2.5 estimation models in Japan. Araki and colleagues modelled the spatial distributions of annual mean PM 2.5 concentrations using LUR and the regression-kriging method [49]. They included AOD at 500 nm from JAXA Satellite Measurements for Environmental Studies (JASMES) as a predictor in the models but excluded AOD during the backward variable selection process. Their final model retained the built-up area ratio, agricultural area ratio, population, distance to the nearest highway, distance to the coastline, precipitation, temperature, wind speed, and longitude as predictors. The leave-one-out cross-validation R 2 (RMSE) was 0.71 (1.2 µg/m 3 ) and 0.81 (1.0 µg/m 3 ) for LUR and regression-kriging, respectively, after removing the spatial outlier [49]. In 2020, Araki and colleagues in Japan used an artificial neural network and included meteorological variables (precipitation, wind speed, temperature, and relative humidity) and land use data (built-up area ratio, green area ratio, population density, distance to coastline, and longitude) as predictors to estimate monthly PM 2.5 concentrations during 1987-2016. Their model achieved a spatial and temporal cross-validation R 2 (RMSE) of 0.76 (1.9 µg/m 3 ) and 0.73 (2.0 µg/m 3 ), respectively [46]. The same group also applied the random forest model with ordinary kriging to estimate monthly PM 2.5 concentrations in Japan during 2010-2015 [50]. They used 5-fold cross-validation to validate the model performance and obtained a cross-validation R 2 of 0.79 with RMSE of 2.1 µg/m 3 for all stations [50]. Their results showed that the top ten important predictors of PM 2.5 were suspended PM, ozone (O 3 ), longitude, month, sulphur dioxides (SO 2 ), temperature, relative humidity, latitude, nitrogen dioxides (NO 2 ), and precipitation. Only a few studies attempted to estimate daily PM 2.5 concentrations in Japan. Shimadera [9]. We integrated MAIAC AOD with meteorological variables and land use data to develop the daily PM 2.5 estimation model in this study. Our model achieved a 10-fold cross-validation R 2 of 0.86 with RMSE of 3.02 µg/m 3 . Our model performance is better than those in studies that merely relied on meteorological variables and land use data to estimate monthly or annual mean PM 2.5 concentrations. This shows that satellite-based AOD can provide additional spatial and temporal information about daily PM 2.5 concentrations and improve the model performance.
Nearly all previous studies used ensemble models (i.e., random forest, gradient boosting, or XGBoost) to develop satellite-based PM 2.5 models, and only a few used a deep learning model or support vector machine (as summarised in Table S5) [17]. Overall, our model performance is comparable to or better than those in these studies.
Our results showed that day of the year was the most important variable in the random forest model for estimating PM 2.5 (relative importance of 19.92%). Early studies that developed satellite-based AOD-PM 2.5 estimation models by using the linear mixed effect model needed to include day-specific intercepts and slopes to achieve high model performance because the relationship between AOD and PM 2.5 may change day-to-day due to differences in mixing height, relative humidity, particle composition, and vertical profiles [13]. Recent studies that applied the machine learning algorithm also found that day of the year (or Julian day) is an important variable for modelling PM 2.5 [22,52,53]. According to the ALE plot ( Figure S8), the relationship between day of the year and PM 2.5 estimates displayed a bimodal pattern with two peaks between Julian day 0-100 (spring) and after Julian day 300 (winter). Day of the year can be a surrogate of seasonal effects. PM 2.5 concentrations in Japan displayed an obvious seasonal variation. Both foreign and domestic sources may cause the seasonal variations of PM 2.5 in Japan. Long-range transportation from Asia (yellow dust storm and anthropogenic aerosols) contributes to high PM 2.5 levels in western Japan, particularly in spring and winter [9,54]. Photochemical smog events frequently occurred in summer can contribute to high concentrations of PM 2.5 in summer [55]. Besides, increase in operation time of coal-burning power plants and heat appliance usage can cause raising PM 2.5 concentrations during winter [55].
AOD was the third most important variable in our model (relative importance of 11.89%). There was a strong and positive relationship between AOD and PM 2.5 according to the ALE plot ( Figure S8). However, the results from previous studies are contradictory. A study conducted on an island, Great Britain of the UK, found that satellite-based AOD 0.47 and 0.55 µm values retrieved from the MAIAC AOD product were not important predictors in a random forest model for estimating ground PM 2.5 concentrations, and the authors concluded that the contribution of satellite-based AOD in the model was very limited [22]. A study in Tehran, Iran also indicated that satellite-based AOD with an extremely high missing rate of 94.09% (3 km Terra Dark Target Algorithm AOD) did not significantly contribute to PM 2.5 estimations [56]. It should be noted that aerosol compositions (i.e., nitrate/sulphate ratio), presence of black carbon, and an arid environment with higher surface reflectance may cause poor correlations between satellite-based AOD and ground PM 2.5 concentrations [57]. Nevertheless, our results demonstrated that satellite-based AOD can provide valuable information for estimating ground PM 2.5 concentrations in Japan.
Meteorological variables had stronger influences on PM 2.5 estimates than population density and land use data according to the variable importance ( Figure S8). A previous study also observed that forest area and population density have lower importance values than meteorological variables due to a lack of temporal variation [53].
Temperature was the second-most-important variable in our model (relative importance of 12.25%). We observed a U-shaped relationship between temperature and PM 2.5 concentrations ( Figure S8). Temperature both positively and negatively influences PM 2.5 .
High temperature can increase PM 2.5 concentrations by promoting photochemical reactions and accelerating formation of PM 2.5 precursors and other secondary pollutants [44,58]. In winter, increased temperature may induce formation of a temperature inversion layer due to low surface temperature and result in accumulation of PM 2.5 [44]. On the other hand, high temperature can enhance thermal activities, facilitate dispersion of PM 2.5 , and increase evaporation loss of PM 2.5 (loss of vapor, ammonium nitrate, and volatile and semi-volatile pollutants), causing decreases in PM 2.5 concentrations [44]. In addition, low temperature can reduce atmospheric convection and lead to increases in PM 2.5 concentrations [44]. We observed negative correlations of PM 2.5 concentrations with relative humidity and precipitation ( Figure S8). These results are unsurprising; PM 2.5 concentrations may decrease as relative humidity increases beyond a threshold due to dry deposition [58]. Additionally, precipitation can decrease PM 2.5 concentrations via wet deposition [58].
Boundary layer height was the seventh most important variable in our model (relative importance of 7.48%). Boundary layer height can affect the vertical distribution of PM; low boundary layer height leads to accumulation of PM 2.5 near the earth's surface [44]. Tsai and colleagues in Taiwan reported that AOD normalised by haze layer height detected by ground-based LiDAR (i.e., the sum of boundary layer height and scaling height) remarkably improves the correlation between AOD and PM 2.5 concentrations [59]. Several studies also included boundary layer height from reanalysed datasets as a predictor in satellite-based AOD-PM 2.5 models [30,60,61]. However, normalising AOD by boundary layer height or including boundary layer height as a predictor did not improve the performance of a linear mixed-effect model in our previous work [23]. It is possible that the linear mixed-effect model cannot capture the complex interactions between AOD and boundary layer height. Nonetheless, we applied a random forest model in this study, and boundary layer height was a critical variable in this model. This demonstrates that random forest is suitable for developing satellite-based AOD-PM 2.5 models, especially when predictors (AOD and meteorological variables in this study) are highly inter-correlated, and their relationships are not linear or monotonic.
This study has several strengths. First, we comprehensively validated the MAIAC AOD with the AERONET AOD across Japan. The results can provide a basic understanding of accuracy of MAIAC AOD in Japan. Second, this is the first study in Japan that used machine learning algorithms with 1 km MAIAC AOD to develop a high spatiotemporal resolution PM 2.5 estimation model that can provide daily estimates of PM 2.5 concentrations at a national scale. We the assessed model performances of two machine learning algorithms-random forest and XGBoost for estimating PM 2.5 concentrations and found that random forest had the best performance. Our cross-validation results showed that the random forest model can provide accurate estimates in areas without air quality monitoring stations from 2011 to 2016. The Japan Environment and Children's Study (JECS), a nationwide birth cohort study, recruited 103,099 participants from 15 regional centres during 2011-2014 [62], and the residential addresses of the participants are being geocoded. The daily PM 2.5 estimates from the satellite-based model can be integrated with the geocodes of the JECS participants to analyse the acute and chronic effects of PM 2.5 on children's health, and to assess the time windows when children are vulnerable. Third, this study provided ALE plots that can describe how predictors influence PM 2.5 estimates in the model. Although ALE is not a new technology, it has rarely been used in previous studies. Here, we demonstrate the importance of ALE can help to quantify the effect of a specific predictor and its direction on PM 2.5 estimates.
This study has some limitations. First, missing AOD values are a major limitation and challenge when building satellite-based models. AOD values can be missing due to the presence of clouds and high surface reflectance, and thus this problem is much worse for islands, which are frequently cloudy and rainy. In the present study, the average coverage rate of the 1-km MAIAC AOD in Japan was only 15.48%. In our previous study in Taiwan, the average coverage rate of the 10 km Dark Target algorithm AOD was 22.99% [23]. Schneider and colleagues reported that the coverage rate of the 1 km MAIAC AOD in Great Britain, UK was 6-13% during 2008-2018 [22]. The coverage rates of AOD in island environments are frequently low based on these studies. Several methods, such as multiple imputation, random forest, and interpolation, were proposed to impute missing AOD values [19,53,63], but they may introduce uncertainties into the model. In 2014, Japan launched a geostationary satellite, Himawari-8, which can continuously provide AOD values every 10 min. A future study could use AOD values derived from this geostationary satellite to overcome the problem of missing AOD values. Second, our land-use data for industrial areas were from 2009 and thus not up-to-date. New data for industrial areas or points of interest in Japan are hard to access, which may constrain the estimation ability of the model. Third, emission inventory data and other criteria pollutants (i.e., nitrogen dioxides, sulphur dioxides, and ozone) were not included in our model. Interactions between pollutants can enhance the production of secondary pollutants, causing increases in PM 2.5 concentrations. A future study could include these data to improve the model performance.

Conclusions
This study applied a random forest model integrating satellite-based AOD, meteorological variables and land use data to generate full-coverage and daily PM 2.5 estimates at 1-km spatial resolution across Japan during 2011-2016. The model performance was excellent, with a 10-fold cross-validation R 2 of 0.86 and an RMSE of 3.02 µg/m 3 , which showed that PM 2.5 estimates derived from the satellite-based model were reliable and accurate. According to the variable importance, the determinants of outdoor PM 2.5 in Japan were day of the year, temperature, satellite-based AOD, relative humidity, wind vectors, boundary layer height, precipitation, surface pressure, and population density. These estimates can be used to assess both the short-term and long-term effects of PM 2.5 on health outcomes in epidemiological studies.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/rs13183657/s1, Table S1. Summary of data sources, coverage years, and spatial and temporal resolutions; Table S2. Descriptive statistics of 16 predictors from 2011 to 2016 (mean ± standard deviation); Table S3. Hyperparameters tuned to optimise the performance of eXtreme Gardient Boosting (XGBoost); Table S4. Descriptive statistics for the Multi-Angle Implementation of Atmospheric Correction (MAIAC) aerosol optical depth (AOD) and AErosol RObotic NETwork (AERONET) AOD; Table S5. Comparison of model performances between this and other studies that applied machine learning algorithms to develop satellite-based PM 2.5 models; Figure S1. Locations of 12 AErosol RObotic NETwork (AERONET) stations in Japan (Hokkaido University, Niigata, Noto, Toyama, Chiba University, Osaka-North, Nara, Osaka, Kobe, Shirahama, Fukuoka and Fukue); Figure S2. Model training (left panel) and 10-fold cross-validation (right panel) of eXtreme Gradient Boosting (XGBoost) model for estimating PM 2.5 concentrations. The red dashed line is the one-one line; Figure S3. The spatial distribution of average residuals (differences between estimated PM 2.5 and in situ measurements) based on the random forest model; Figure S4 [64,65] are cited in the supplementary materials.