Improving Satellite-Driven PM 2.5 Models with VIIRS Nighttime Light Data in the Beijing–Tianjin–Hebei Region, China

: Previous studies have estimated ground-level concentrations of particulate matter 2.5 (PM 2.5 ) using satellite-derived aerosol optical depth (AOD) in conjunction with meteorological and land use variables. However, the impacts of urbanization on air pollution for predicting PM 2.5 are seldom considered. Nighttime light (NTL) data, acquired with the Visible Infrared Imaging Radiometer Suite (VIIRS) aboard the Suomi National Polar-orbiting Partnership (S-NPP) satellite, could be useful for predictions because they have been shown to be good indicators of the urbanization and human activity that can affect PM 2.5 concentrations. This study investigated the potential of incorporating VIIRS NTL data in statistical models for PM 2.5 concentration predictions. We developed a mixed-effects model to derive daily estimations of surface PM 2.5 levels in the Beijing–Tianjin–Hebei region using 3 km resolution satellite AOD and VIIRS NTL data. The results showed the addition of NTL information could improve the performance of the PM 2.5 prediction model. The NTL data revealed additional details for predication results in areas with low PM 2.5 concentrations and greater apparent seasonal variation due to the seasonal variability of human activity. Comparison showed prediction accuracy was improved more substantially for the model using NTL directly than for the model using the vegetation-adjusted NTL urban index that included NTL. Our ﬁndings indicate that VIIRS NTL data have potential for predicting PM 2.5 and that they could constitute a useful supplemental data source for estimating ground-level PM 2.5 distributions.


Introduction
In recent decades, rapid economic development in China has caused severe environmental pollution, particularly air pollution [1]. As a major air pollutant in many regions of China, fine particulate matter (aerodynamic diameter < 2.5 µm; PM 2.5 ) is associated with various adverse health outcomes including cardiovascular and respiratory diseases [2,3] and consequently, it has received considerable attention [4]. Accurate assessment of air pollution levels and efficacious solutions to public health concerns require measurement of PM 2.5 concentrations at high temporal and spatial resolutions [5,6].
Measurements of PM 2.5 concentrations can be acquired by ground-based monitoring networks and satellite remote sensing. Although PM 2.5 measurements obtained from ground-based monitoring sites are accurate, such networks are typically sparse, which means the spatial coverage of routine measurements is limited [5,7,8]. Because of the large spatial coverage and reliability of repeated measurements, satellite remote sensing provides a potentially cost-effective method for predicting PM 2.5 concentrations where ground-monitoring stations are unavailable or sparsely distributed [7,9,10]. Satellite-derived aerosol optical depth (AOD), which reflects the integrated quantity of particles in This study developed a mixed-effects model to estimate ground-level PM2.5 concentrations using the MODIS 3 km resolution AOD products, together with meteorological and VIIRS NLT data. The objectives of this study were to investigate how NTL data could be incorporated into a statistical model for predicting ground-level PM2.5 concentrations, and to assess whether this parameter could improve PM2.5 prediction accuracy significantly. The remainder of this paper is organized as follows. Section 2 describes the study area and the data used. The method to establish the relationship between AOD and PM2.5 concentration using NTL data is presented in Section 3. The results of the performance of the model are given in Section 4, while discussion of the results and the findings is presented in Section 5. Finally, the derived conclusions are stated in Section 6.

Study Area
The study area was the Beijing-Tianjin-Hebei (BTH) region, which is situated in northern China ( Figure 1). It includes Beijing, Tianjin, and eight neighboring cities in China's Hebei Province. This area covers the Greater Beijing Economic Hub and its surrounding area, and it is also the core region of the Bohai Bay Trade Hub, having long been the most important political, economic, and cultural center of northern China. Following rapid development, the BTH region is also one of the most polluted city clusters in China [20].

Ground PM2.5 Measurements
Hourly PM2.5 concentrations, observed at 79 monitoring stations in the BTH region from 2016, were obtained from the website of the China Environmental Monitoring Center (http://106. 37.208.233:20035/). We also collected additional data from the website of the Beijing

Ground PM 2.5 Measurements
Hourly PM 2.5 concentrations, observed at 79 monitoring stations in the BTH region from 2016, were obtained from the website of the China Environmental Monitoring Center (http://106. 37.208.233: 20035/). We also collected additional data from the website of the Beijing Municipal Environmental Monitoring Center (http://zx.bjmemc.com.cn/). Overall, data from 114 monitoring sites were included in the present study ( Figure 1). The ground concentrations measured between 10:00 and 14:00 local time (LT) at each site were averaged to derive daily PM 2.5 concentrations that matched the overpass time of the MODIS Terra and Aqua satellites. The study period comprised the entire year of 2016, spanning 366 days.

MODIS 3 km AOD Products and Calibration
We obtained MODIS Collection 6 Level 2 aerosol products with 3 km spatial resolution for 2016 from NASA Level 1 and the Atmosphere Archive and Distribution System (accessed at https: //ladsweb.nascom.nasa.gov/data/search.html). The retrieval algorithm of the higher resolution AOD product is similar to that of the 10 km AOD product, but it averages 6 × 6 pixels in a single retrieval box rather than 20 × 20 pixels after cloud screening and other surface mask processes [16,31]. The MODIS 3 km AOD products have been validated by ground-based data from six Aerosol Robotic Network (AERONET) sites in China and the spatiotemporal variations of 3 km AOD showed good agreement with the AERONET [32].
In this study, Aqua (overpass at~13:30 LT) and Terra (overpass at~10:30 LT) MODIS AOD data at 550 nm were first combined and averaged to improve the spatial coverage. In grid cells that had both Terra and Aqua AODs, the averaged value represented the mean of the AOD distribution from 10:00 to 14:00 LT. However, there were many days when only one of the two retrievals was available. Therefore, when combining the two MODIS products, the averaged AOD of a grid cell was biased toward either the morning or afternoon condition. To overcome this bias, Lee et al. [10] used the average Terra AOD/Aqua AOD ratio to estimate the missing AOD values. In this study, to exploit the measurements of both satellites fully, we fitted a simple linear regression to define the relationship between the daily mean AOD values of Terra and Aqua [6]. The relationship between the Aqua and Terra AODs for different seasons was found similar and therefore, we built the regression model using annual data. The derived regression equation was as follows: where τ is the AOD. A value of R 2 = 0.77 was obtained for both regression models. The missing AOD value was estimated using the above regression equation and then averaged. Consequently, each AOD grid cell contained a mean value that better represented the average conditions from 10:00 to 14:00 LT.

NDVI Data
We considered the Normalized Difference Vegetation Index (NDVI) in the model development process because a previous study has reported that forest cover can affect the AOD-PM 2.5 relationship [5]. The Terra MODIS NDVI data used in this study constituted the MOD13A2 products, which are provided every 16 days at 1 km spatial resolution. MOD13A2 data are computed from atmospherically corrected bidirectional surface reflectance [33]. The 16-day mean values of NDVI were used as covariates in the PM 2.5 model.

VIIRS NTL and Urban Index Based on NTL
In this study, we examined the impact of urbanization on the relationship between PM 2.5 and AOD; the urbanization indicator was expressed by the VIIRS NTL. The VIIRS DNB monthly composites, as a new generation of NTL image products, have been produced and released by the National Oceanic and Atmospheric Administration's National Centers for Environmental Information [34,35]. As a low NTL detecting band, the DNB can detect NTL within the radiation range of 3.00 nW · cm −2 · sr −1 to 0.02 W · cm −2 · sr −1 , with a span from 75 • N to 65 • S, 3000 km swath, and a 12-h revisiting period. The Remote Sens. 2017, 9, 908 5 of 14 spatial resolution is 15 arcseconds (742 m), with a spectral range of 0.5-0.9 µm [36,37]. Monthly VIIRS NTL values were used in this study.
As the primary purpose of this study was to investigate how/whether NTL data could be used to improve prediction accuracy of PM 2.5 concentrations, an urban index using VIIRS NTL was also used in the model to compare model performances with different NTL expressions. We selected the Vegetation-Adjusted NTL Urban Index (VANUI) as another urbanization indicator. Following Zhang et al., the VANUI was calculated in this study based on both the VIIRS NTL data and the NDVI data [38]. The equation used for calculating the VANUI is as follows: where VANU I (i,j) is the VANUI value for the i-th pixel in the j-th evaluation area, NDV I (i,j) is the NDVI value for the i-th pixel, and V I IRS nor (i,j) is the normalized DN value of NTL for the i-th pixel. V I IRS nor (i,j) is calculated using the following equation: W cm sr     , with a span from 75°N to 65°S, 3000 km swath, and a 12-h revisiting period. The spatial resolution is 15 arcseconds (742 m), with a spectral range of 0.5-0.9 μm [36,37]. Monthly VIIRS NTL values were used in this study.
As the primary purpose of this study was to investigate how/whether NTL data could be used to improve prediction accuracy of PM2.5 concentrations, an urban index using VIIRS NTL was also used in the model to compare model performances with different NTL expressions. We selected the Vegetation-Adjusted NTL Urban Index (VANUI) as another urbanization indicator. Following Zhang et al., the VANUI was calculated in this study based on both the VIIRS NTL data and the NDVI data [38]. The equation used for calculating the VANUI is as follows: (1 )

Meteorological Data
The meteorological data were obtained from the Goddard Earth Observing System Data Assimilation System Forward Processing (GEOS-5 FP). These data included planetary boundary layer height (PBLH, m), relative humility (RH, %), wind speed (WS, m/s), and surface pressure (PS, hPa). The GEOS-5 FP is the latest version of the GEOS-5 meteorological data, with spatial resolution of 0.25 • (latitude) × 0.3125 • (longitude) in a nested grid that covers China, which has been available since April 2012 [39]. GEOS-5 FP data from 2016 were used in this study for model development and hourly GEOS-5 FP data from 10:00 to 14:00 LT were averaged as daytime meteorological parameters to match the satellite overpass times.
All data were re-projected using the Albers equal area conic projection and resampled to spatial resolution of 3 km. Then, all the data were clipped following administrative boundary data of the study area. Finally, all variables in 2016 were matched by grid cell and date for model fitting.

Model Development
Previous studies have shown that time-varying parameters such as RH, PM 2.5 vertical and diurnal concentration profiles, and PM 2.5 optical properties can influence the PM 2.5 -AOD relationship [10,17]. Therefore, we developed a linear mixed-effects (LME) model [10] to account for daily variations of the AOD-PM 2.5 relationship. Considering the regional spatial scale of this study, we hypothesize that these time-varying parameters exhibit little spatial variability and consequently, the PM 2.5 -AOD relationship has minimal spatial variation on any given day. Daily variations of the relationship between additional variables and PM 2.5 were not considered in our model because they led to over-fitting of the models. In addition, whether a predictor should be incorporated needs to be tested in the model. We tested other meteorological parameters (e.g., air temperature) and land use types not included in the final models, and no significant improvement was found in the regression performance.
We first fitted three LME models to compare the performances of models with and without NTL predictors: the LME model only with meteorological data (referred to as the Meteor_LME, Equation (4)), LME model without VIIRS NTL data (referred to as the NDVI_LME, Equation (5)), and LME model with VIIRS NTL data (referred to as the NTL_LME, Equation (6)). In addition, we also fitted the LME model using VANUI (referred to as the VANUI_LME, Equation (7)) instead of NTL and NDVI to establish the form of urbanization indicator using NTL that benefitted model performance the most. These models can be expressed as follows: where PM s,t (µg/m 3 ) is the averaged observed PM 2.5 concentration at grid cell s on day t; b 0 and b 0,t are the fixed and day-specific random intercepts, respectively, b i (i = 1, 2 . . . . . . ) are the fixed slopes for the predictors, and b 1,t is the random slope for AOD at each day; AOD s, t refers to the MODIS AOD value (unitless) at grid cell s on day t; PBLH s,t (m), RH s,t (%), WS s,t (m/s), and PS s,t (hPa) are the meteorological parameters at grid cell s on day t (definitions in Section 2.4); NTL nor s,t (unitless) is the normalized value of VIIRS NTL data at grid cell s on day t, NDVI s,t (unitless) is the MODIS NDVI value at grid cell s on day t, and VANUI s,t (unitless) is the index value obtained by combining the NTL and NDVI values (definitions in Section 2.3.3) at grid cell s on day t. In these statistical models, the AOD fixed effect represents the average effect of AOD on PM 2.5 for all studied days. The AOD random effects explain the daily variability in the PM 2.5 -AOD relationship. These equations were applied to the entire fitting dataset to generate the fixed effects intercept and slopes for all the days and the random effects intercept and slopes for each individual day.

Model Validation
To validate the performances of the four models, commonly used statistical indicators between fitted values from the models and observations were calculated. We calculated the coefficient of determination (R 2 ) to indicate the degree of agreement of fitting between predicted and observed PM 2.5 concentrations. The mean prediction error (MPE) and the root mean square prediction error (RMSE) were adopted to evaluate the prediction accuracy of each of the four models. The MPE and RMSE are defined as and whereŷ i is the predicted PM 2.5 concentrations, y i is the observed values, and n is the total number of observations. The goodness of fit of each model was checked using the Akaike information criterion (AIC) and the Bayes information criterion (BIC) [40]. Smaller values of the AIC and the BIC indicate better fitting of the model to the data. In addition, a 10-fold cross-validation (CV) method was applied to assess potential model over-fitting [41]. The dataset was first split into 10 folds with approximately 10% of the total data points in each fold. In each round of the CV, one fold was used as the test data, the remaining nine folds were used to fit the model, and predictions were made on the held-out test fold. This process was repeated 10 times until each fold had been tested. Furthermore, the same statistics (e.g., R 2 , MPE, and RMSE) were calculated for the cross-validation results.

Descriptive Statistics
The descriptive statistics of the dependent and predictor variables used in the model fitting are presented in Table 1. There were 3674 matched data records used for the model fitting. The overall mean PM 2.5 concentration is 55.20 µg/m 3 , which is higher than China's National Ambient Air Quality Standard (NAAQS) of 35.00 µg/m 3 for annual mean PM 2.5 , and the standard deviation (SD) is 42.07 µg/m 3 . The mean value of MODIS AOD is 0.71, although some sites have high AOD values of >2.00. The normalized NTL for 2016 has a mean value of 0.22 and the mean NDVI value is 0.42, which means there is relatively high vegetation cover in the BTH region. The RH has a relatively small mean value of 33%, which implies a relatively dry atmospheric environment. The PBLH (1636.88 ± 623.41 m) is more variable in comparison with other meteorological variables.  Table 2 summarizes the fixed effect and random estimates of the four LME models with different combinations of predictors. Daily AOD-PM 2.5 relationships of 203 days were calculated using the LME models. The fixed terms of the intercept and slopes for the Meteor_LME and VANUI_LME models are different to these for the NDVI_LME and NTL_LME models, while the standard errors (SEs) of the fixed terms for the four models are similar. The SEs of the daily specific intercepts and slopes for the four models are about 23.00 and 26.00, respectively. The fixed effects and random effects results show all four LME models can reflect the temporal variations of AOD-PM 2.5 relationships, irrespective of whether the predictor variables are different. Table 2 also shows that incorporating the VIIRS NTL data into the LME model produced the smallest values of the AIC and the BIC, indicating better model fitting.  Figure 3 shows the model fitting and CV results for the different models. For the four model fittings, the overall values of R 2 between the predicted and observed PM 2.5 concentrations are all >0.75, with RMSEs of around 20 µg/m 3 and MPEs of around 14.5. The CV R 2 values of the four models are all >0.70, which reflects a decrease of 0.04 compared with the model fitting results. The CV RMSE values increase by about 1 µg/m 3 for the four models, which suggests the models are not substantially over-fitted. Comparison of the CV statistics of the Meteor_LME model and the NDVI_LME model, by adding the NDVI predictor to the model, shows the overall prediction accuracy was improved. The CV RMSPE decreased by 0.06 µg/m 3 , representing a 0.29% reduction. The NTL_LME fitting and CV R 2 values are 0.7512 and 0.7092 (Figure 3c), respectively, which indicate this model performed best of all the models. Comparison of the other three models (Figure 3a,b,d) with the NTL_LME model (Figure 3c) shows the NTL model marginally increased the R 2 values. The CV results indicate all models have a trend to overestimate the PM 2.5 concentrations in relation to relatively low PM 2.5 observations and to underestimate the PM 2.5 values in relation to relatively high PM 2.5 observations. The Meteor_LME model produced the worst model fitting and it yielded greater bias than the other models.

Model Fitting and Validation
Both the NTL_LME and VANUI_LME models included NTL data but their results are different. The NTL_LME results improved the performance of the Meteor_LME model, but the VANUI_LME results had little effect on the Meteor_LME model. This result suggests that the form of the urbanization indicator using NTL could affect the results of the models, and that NTL as a predictor variable added to the model directly performed better than using the form of the urban index as a predictor variable.
Both the NTL_LME and VANUI_LME models included NTL data but their results are different. The NTL_LME results improved the performance of the Meteor_LME model, but the VANUI_LME results had little effect on the Meteor_LME model. This result suggests that the form of the urbanization indicator using NTL could affect the results of the models, and that NTL as a predictor variable added to the model directly performed better than using the form of the urban index as a predictor variable.   Figure 4 shows the annual mean PM 2.5 concentrations predicted by the different models. The results show that the patterns of PM 2.5 estimates from the four models are largely similar. Spatially, there is a southeast-to-northwest gradient in the predicted PM 2.5 distribution. High concentrations are found in Beijing and the south of Hebei Province, while low concentrations appear in the north of Hebei Province. From the details of Beijing, highly urbanized and populated areas such as Haidian, Shijingshan, and Fengtai exhibit higher PM 2.5 levels compared with rural areas such as Huairou.

Annual and Seasonal Mean PM 2.5 Concentrations
Furthermore, higher PM 2.5 levels are predicted over southern regions and lower values predicted over northern and western suburban districts.
Comparison of the PM 2.5 concentrations estimated by the different models shows those models using NTL, either integrated as an urban index (VANUI_LME) or directly as a variable (NTL_LME), generated greater detail of the predicted surface than those not using NTL data (i.e., Meteor_LME and NDVI_LME), especially in areas with low PM 2.5 values (e.g., Chengde and Zhangjiakou). In this region, PM 2.5 estimates generated from the models with NTL were higher than those estimated from the models without NTL. In addition, the NTL_LME and NDVI_LME models provided higher predictions than the Meteor_LME and VANUI_LME models in Shijiazhuang and Baoding. Figure 4 shows the annual mean PM2.5 concentrations predicted by the different models. The results show that the patterns of PM2.5 estimates from the four models are largely similar. Spatially, there is a southeast-to-northwest gradient in the predicted PM2.5 distribution. High concentrations are found in Beijing and the south of Hebei Province, while low concentrations appear in the north of Hebei Province. From the details of Beijing, highly urbanized and populated areas such as Haidian, Shijingshan, and Fengtai exhibit higher PM2.5 levels compared with rural areas such as Huairou. Furthermore, higher PM2.5 levels are predicted over southern regions and lower values predicted over northern and western suburban districts.
Comparison of the PM2.5 concentrations estimated by the different models shows those models using NTL, either integrated as an urban index (VANUI_LME) or directly as a variable (NTL_LME), generated greater detail of the predicted surface than those not using NTL data (i.e., Meteor_LME and NDVI_LME), especially in areas with low PM2.5 values (e.g., Chengde and Zhangjiakou). In this region, PM2.5 estimates generated from the models with NTL were higher than those estimated from the models without NTL. In addition, the NTL_LME and NDVI_LME models provided higher predictions than the Meteor_LME and VANUI_LME models in Shijiazhuang and Baoding. Because of limited satellite data in winter (January, February, and December), the data were grouped into two seasons: warm (15 April-14 October) and cold (15 October-14 April). Figure 5 shows the two seasonal mean distributions of PM2.5 concentration. The PM2.5 pattern of the BTH region exhibits distinct seasonal variational characteristics, with the highest pollution occurring in the cold season and the lowest pollution occurring in the warm season (except for Beijing). The predicted PM2.5 concentrations within the warm season range from 8 to 190 μg/m 3 , while the predicted PM2.5 concentrations within the cold season range from 2 to 300 μg/m 3 . Such seasonal variation can be attributed to the mixed contributions from meteorological conditions and local emission sources.
Comparison of Figure 5 shows the mean PM2.5 estimates generated from models with NTL have higher values in the cold season and lower values in the warm season than models without NTL. This indicates that models with NTL exhibit more significant seasonal variation than models without NTL, especially in northern areas of the BTH region. For example, the average PM2.5 concentration within the cold season predicted by the NTL_LME model is 54.65 μg/m 3 , which is higher than the PM2.5 concentration predicted by the NDVI_LME model (53.09 μg/m 3 ). Conversely, Because of limited satellite data in winter (January, February, and December), the data were grouped into two seasons: warm (15 April-14 October) and cold . Figure 5 shows the two seasonal mean distributions of PM 2.5 concentration. The PM 2.5 pattern of the BTH region exhibits distinct seasonal variational characteristics, with the highest pollution occurring in the cold season and the lowest pollution occurring in the warm season (except for Beijing). The predicted PM 2.5 concentrations within the warm season range from 8 to 190 µg/m 3 , while the predicted PM 2.5 concentrations within the cold season range from 2 to 300 µg/m 3 . Such seasonal variation can be attributed to the mixed contributions from meteorological conditions and local emission sources.
Comparison of Figure 5 shows the mean PM 2.5 estimates generated from models with NTL have higher values in the cold season and lower values in the warm season than models without NTL. This indicates that models with NTL exhibit more significant seasonal variation than models without NTL, especially in northern areas of the BTH region. For example, the average PM 2.5 concentration within the cold season predicted by the NTL_LME model is 54.65 µg/m 3 , which is higher than the PM 2.5 concentration predicted by the NDVI_LME model (53.09 µg/m 3 ). Conversely, the average PM 2.5 concentration within the warm season predicted by the NTL_LME model is 47.95 µg/m 3 , which is lower than the PM 2.5 concentration predicted by the NDVI_LME model (48.44 µg/m 3 ). In addition, the differences between PM 2.5 concentrations estimated by models with and without NTL are larger in the cold season than in the warm season. The level of coal and biomass combustion associated with heating increases in winter, which represents human activities. NTL is used as an indicator of urban areas, human activities, and population density. Consequently, the contribution of NTL to predictions of PM 2.5 concentrations is correspondingly higher in winter.
Remote Sens. 2017, 9,908 11 of 14 the average PM2.5 concentration within the warm season predicted by the NTL_LME model is 47.95 μg/m 3 , which is lower than the PM2.5 concentration predicted by the NDVI_LME model (48.44 μg/m 3 ). In addition, the differences between PM2.5 concentrations estimated by models with and without NTL are larger in the cold season than in the warm season. The level of coal and biomass combustion associated with heating increases in winter, which represents human activities. NTL is used as an indicator of urban areas, human activities, and population density. Consequently, the contribution of NTL to predictions of PM2.5 concentrations is correspondingly higher in winter.

Discussion
This paper evaluated the performance of VIIRS NTL data when used as a predictor for PM2.5 concentration estimations. In brief, it was proven that VIIRS NTL data have considerable potential to provide additional predictive power for estimations of PM2.5 concentrations.
Our results showed that the incorporation of NTL data as a predictor could improve PM2.5 prediction accuracy, and that NTL data might partially account for PM2.5 variability in the models, which suggests that elevated PM2.5 levels cannot be explained fully by AOD. The experimental results also showed that models with NTL exhibited more significant seasonal variation than models

Discussion
This paper evaluated the performance of VIIRS NTL data when used as a predictor for PM 2.5 concentration estimations. In brief, it was proven that VIIRS NTL data have considerable potential to provide additional predictive power for estimations of PM 2.5 concentrations.
Our results showed that the incorporation of NTL data as a predictor could improve PM 2.5 prediction accuracy, and that NTL data might partially account for PM 2.5 variability in the models, which suggests that elevated PM 2.5 levels cannot be explained fully by AOD. The experimental results also showed that models with NTL exhibited more significant seasonal variation than models without NTL. Because NTL can reflect human activity, which varies seasonally, the impact of human activity on the AOD-PM 2.5 relationship can account for the seasonal variability exhibited by those models with NTL. Furthermore, the differences between the PM 2.5 concentrations estimated by models with and without NTL were larger in the cold season than in the warm season. In winter, the main contribution to elevated PM 2.5 concentrations is from coal and biomass combustion for heating. This information might be not captured by MODIS AOD, whereas NTL as an indicator of human activity can reflect the effect of heating to some extent. Thus, NTL data might introduce into the models additional PM 2.5 emission sources missed by MODIS AOD, which could lead to predictions of PM 2.5 that are more accurate. As the impact of human activity on PM 2.5 concentrations is more apparent in the cold season than in the warm season, the differences between models with and without NTL are larger in the cold season.
It should be noted that the form in which NTL data are incorporated within the models could affect the applicability of NTL data in air pollution monitoring studies. We compared the performances of models with NTL data and with VANUI, which is an urban index using NTL and NDVI. The model based on NTL and the NDVI variable appeared to provide better estimation performance than the model based solely on VANUI. This might be because the VANUI focuses on inner-urban variations more than on urbanization characteristics [38]. Further research is needed to explore the potential data fusion approaches for better use of NTL data in statistical models for PM 2.5 concentration predictions.
It should also be mentioned that NTL data can be acquired easily and are free to download, while it is difficult to obtain other urban indicator data such as accurate spatial distributions of population. In addition, VIIRS can provide data from 2013 to the present, and the DMSP/OLS NTL has a long historical archive. Our statistical model for predicting PM 2.5 concentrations could be further used to produce long-term time series of serious PM 2.5 exposure, which could then be applied in long-term epidemiological studies. Furthermore, we used high spatial resolution MODIS 3 km AOD data to estimate high spatial resolution daily average PM 2.5 concentrations. Higher resolution data could help reveal additional local details of the relationship and thus, the acquisition of such data is desirable.
Although the developed LME models with NTL produced promising results, some limitations remain. We combined Terra MODIS data and Aqua MODIS data to improve the spatial coverage of the MODIS 3 km AOD products, but there were still a number of missing data. Further studies could introduce additional estimators to fill the gaps in areas where AOD is not retrieved. Furthermore, the PM 2.5 -AOD relationship varies in both space and time, but the LME model used in our study incorporates only the temporal variation in the PM 2.5 -AOD association. We will develop new statistical models with NTL that can include both the spatial and the temporal variabilities in the PM 2.5 -AOD relationship.

Conclusions
In this paper, we assessed the applicability of VIIRS/DNB NTL data in estimating surface PM 2.5 concentrations using an LME model in the BTH region of China in 2016. The results indicated that NTL data were an effective predictor of PM 2.5 concentrations. Of the models examined, the model incorporating NTL provided better performance and illustrated more detail in the predicted surface, especially in areas with low PM 2.5 values. Our analysis also showed that the introduction of NTL into PM 2.5 prediction models increased the seasonal variability. The results of this study provide insight into the possibility of enhancing predictions of surface PM 2.5 concentrations by incorporating NTL data. The adoption of our method could help in investigations of the impact of urbanization on PM 2.5 concentrations.