Estimating Ground Level NO 2 Concentrations over Central-Eastern China Using a Satellite-Based Geographically and Temporally Weighted Regression Model

People in central-eastern China are suffering from severe air pollution of nitrogen oxides. Top-down approaches have been widely applied to estimate the ground concentrations of NO2 based on satellite data. In this paper, a one-year dataset of tropospheric NO2 columns from the Ozone Monitoring Instrument (OMI) together with ambient monitoring station measurements and meteorological data from May 2013 to April 2014, are used to estimate the ground level NO2. The mean values of OMI tropospheric NO2 columns show significant geographical and seasonal variation when the ambient monitoring stations record a certain range. Hence, a geographically and temporally weighted regression (GTWR) model is introduced to treat the spatio-temporal non-stationarities between tropospheric-columnar and ground level NO2. Cross-validations demonstrate that the GTWR model outperforms the ordinary least squares (OLS), the geographically weighted regression (GWR), and the temporally weighted regression (TWR), produces the highest R2 (0.60) and the lowest values of root mean square error mean (RMSE), absolute difference (MAD), and mean absolute percentage error (MAPE). Our method is better than or comparable to the chemistry transport model method. The satellite-estimated spatial distribution of ground NO2 shows a reasonable spatial pattern, with high annual mean values (>40 μg/m3), mainly over southern Hebei, northern Henan, central Shandong, and southern Shaanxi. The values of population-weight NO2 distinguish densely populated areas with high levels of human exposure from others.


Introduction
High ground level nitrogen oxides (NO x = NO + NO 2 ) are identified to be deleterious to human health, including decreased lung function and an increased risk of respiratory symptoms [1,2].In addition, NO x can also produce other photochemical pollutants like O 3 in photochemical reactions, and acts as a gaseous precursor of aerosols and acid rain.Thus, the NO x concentration has been included in multi-pollutant health indices [3] and its monitoring with complete spatial coverage is needed for exposure assessment.Since 1995, a series of satellites sensors, e.g., the Global Ozone Monitoring Experiment (GOME), the Scanning Imaging Absorption Spectrometer for Atmospheric Cartography (SCIAMACHY), and the Ozone Monitoring Instrument (OMI) have been successfully used to retrieve vertical NO 2 columns [4][5][6][7][8].A dramatic increase in tropospheric NO 2 columns was revealed by the GOME and SCIAMACHY observations over China [9][10][11][12], the world's largest developing country along with the fastest growing economy.
Given that the existing ambient monitoring stations are sparse and unevenly distributed, there is a growing interest in the top-down satellite approach to obtain timely map of the spatial variations of surface concentrations of NO 2 .A close relationship between ground level NO 2 concentrations and satellite-retrieved tropospheric NO 2 columns is expected based on two facts: (1) ground level NO 2 accounts for the majority of tropospheric NO 2 columns since human activities are their main source; and (2) the short lifetime of near-surface NO 2 results in little transport, both vertically and horizontally [13].Petritoli et al. [14] demonstrated a significant correlation between in situ NO 2 measurements and the GOME tropospheric NO 2 columns.Recently, satellite observations were combined with land use regression models to provide spatio-temporally resolved ambient NO 2 [15][16][17].In addition, an approach proposed by Lamsal et al. [18] that combines the vertical profiles of NO 2 generated by the chemical transport model and satellite tropospheric NO 2 columns, has been widely used to estimate ground level NO 2 concentrations [19,20].However, the emission inventories used for the model simulations are based on outdated statistical data about human activities.These model-based profiles may not capture the actual vertical distribution of NO 2 , especially where anthropogenic NO x emissions are undergoing rapid changes such as in China [21].Kim et al. [22] estimated the surface NO 2 volume mixing ratio by using multiple regression models with OMI data.
In this study, a geographically and temporally weighted regression (GTWR) model is introduced to estimate the ground level NO 2 concentrations by using the OMI tropospheric NO 2 columns over central-eastern China.The GTWR model is adapted from the geographically weighted regression (GWR) model [23][24][25][26] by taking into account spatio-temporal non-stationarity, which has been proven to effectively establish the relation between satellite-retrieved aerosol optical depth and fine particulate matter (PM 2.5 ) [27,28].Furthermore, population-weighted ground level NO 2 concentrations are calculated to evaluate population exposure levels in different regions.

Study Area
This study focuses on the central-eastern China with a geographic scope of 20 • N-45 • N and 105 • E-124 • E (major populated areas in China, see left panel in Figure 1).The study area covers 20 province-level administrative units in mainland China, including the regions of the North China Plain, Yangtze River Delta, and Pearl River Delta that are most polluted.715 ambient monitoring stations are located in this study area (see the right panel in Figure 1).

OMI Tropospheric NO2 Columns
OMI is a Dutch-Finnish nadir-viewing hyperspectral instrument onboard the Earth Observing System Aura satellite in a Sun-synchronous orbit with an equatorial crossing time of approximately 13:45 local time.It measures sunlight backscattered radiances from the Earth in three channels covering a wavelength range of 270 to 500 nm (UV-1: 270 to 310 nm; UV-2: 310 to 365 nm; and, visible: 365 to 500 nm) at a spectral resolution of 0.45 to 0.63 nm [29].OMI makes simultaneous measurements in a swath of width 2600 km, divided into 60 fields of view (FoVs).The FoVs vary in size from ~13 × 26 km near nadir to ~40 × 250 km at the outermost FoVs.The OMI measurements in the spectral range 402-465 nm are used to retrieve the NO2 columns.First, NO2 slant columns are determined from the OMI calibrated earthshine radiance spectra by using the differential optical absorption spectroscopy (DOAS) algorithm [30].Second, the slant columns are then converted into the vertical columns using air mass factors (AMFs) calculated from radiative transfer models.Finally, the stratospheric and tropospheric column amounts are derived separately under the assumption that the two quantities are largely independent [31].
Here, we used the Version 3 Aura OMI NO2 Standard Product (OMNO2) available from the NASA Goddard Earth Sciences Data and Information Services Center (http://disc.gsfc.nasa.gov/Aura/OMI/omno2_v003.shtml).The major improvements include: (1) an improved spectral fitting algorithm for retrieving slant column densities, including the use of monthly mean solar spectral irradiances; (2) improved Global Modeling Initiative model-based monthly a priori NO2 and temperature profiles [32].For further details, please refer to [33].The main error sources in determining tropospheric NO2 columns are associated with uncertainties in the surface albedo, aerosols, cloud interference, and the NO2 vertical profile [34][35][36][37].Overall, OMI retrievals tend to be lower in urban regions and higher in remote areas, but generally agree with other measurements within ±20% [38].

OMI Tropospheric NO 2 Columns
OMI is a Dutch-Finnish nadir-viewing hyperspectral instrument onboard the Earth Observing System Aura satellite in a Sun-synchronous orbit with an equatorial crossing time of approximately 13:45 local time.It measures sunlight backscattered radiances from the Earth in three channels covering a wavelength range of 270 to 500 nm (UV-1: 270 to 310 nm; UV-2: 310 to 365 nm; and, visible: 365 to 500 nm) at a spectral resolution of 0.45 to 0.63 nm [29].OMI makes simultaneous measurements in a swath of width 2600 km, divided into 60 fields of view (FoVs).The FoVs vary in size from ~13 × 26 km near nadir to ~40 × 250 km at the outermost FoVs.The OMI measurements in the spectral range 402-465 nm are used to retrieve the NO 2 columns.First, NO 2 slant columns are determined from the OMI calibrated earthshine radiance spectra by using the differential optical absorption spectroscopy (DOAS) algorithm [30].Second, the slant columns are then converted into the vertical columns using air mass factors (AMFs) calculated from radiative transfer models.Finally, the stratospheric and tropospheric column amounts are derived separately under the assumption that the two quantities are largely independent [31].
Here, we used the Version 3 Aura OMI NO 2 Standard Product (OMNO2) available from the NASA Goddard Earth Sciences Data and Information Services Center (http://disc.gsfc.nasa.gov/Aura/OMI/omno2_v003.shtml).The major improvements include: (1) an improved spectral fitting algorithm for retrieving slant column densities, including the use of monthly mean solar spectral irradiances; (2) improved Global Modeling Initiative model-based monthly a priori NO 2 and temperature profiles [32].For further details, please refer to [33].The main error sources in determining tropospheric NO 2 columns are associated with uncertainties in the surface albedo, aerosols, cloud interference, and the NO 2 vertical profile [34][35][36][37].Overall, OMI retrievals tend to be lower in urban regions and higher in remote areas, but generally agree with other measurements within ±20% [38].
The data were filtered using a number of criteria [39] to ensure retrieval quality including: (1) cloud radiance fraction <0.3, (2) surface albedo <0.3, (3) solar zenith angles <85 • , (4) 10 < cross-track positions < 50, and (5) root mean squared error of fit <0.0003.In addition, the cross track pixels affected by row anomaly (http://www.knmi.nl/omi/research/product/rowanomaly-background.php) were excluded, which was first noticed in the data in June 2007.Then, the NO 2 tropospheric column densities from the Level-2 OMNO2 Swath product were binned on to a 0.1 × 0.1 • grid by calculating the area-weighted averages at each grid cell.

Ambient Monitoring Station Data
The Ministry of Environmental Protection of Republic of China has built 1497 ambient monitoring stations over 367 cities in order to assess the air quality in China.Hourly mean concentrations of air pollutants including PM 2.5 , PM 10 , NO 2 , SO 2 , and O 3 are available since 2013 in the national air quality publishing platform (http://106.37.208.233:20035/).In this study, hourly mean ground-based NO 2 concentrations of 715 stations in central-eastern China from 1 May 2013 to 30 April 2014 (13:00-15:00 local time) were included.The locations of these stations are shown in Figure 1.

Meteorological Data
In order to improve the performance of our regression model, a number of meteorological parameters such as air temperature, relative humidity, planetary boundary layer height, wind speed, and air pressure from the Weather Research & Forecasting Model (WRF, version 3.4.1)were used.NCEP FNL Operational Model Global Tropospheric Analyses dataset of 1 × 1 • resolution (http://rda.ucar.edu/dsszone/ds083.2/)was adopted in the WRF model.The WRF model is a mesoscale numerical weather prediction system designed for both atmospheric research and operational forecast, and serves as a wide range of meteorological applications across scales from tens of meters to thousands of kilometers.The nested domain scheme with 30 km horizontal grid space of WRF output centered at 115 • E, 32.5 • N was adopted, and the temporal resolution of WRF outputs was set 1 h intervals.The number of altitude levels is 30 and the top-level pressure is 50 hPa.The physical options used in WRF include the single-moment 3-class (WSM3) microphysics, the Yonsei University (YSU) PBL scheme, the Rapid Radiative Transfer Model (RRTM) longwave and Dudhia shortwave radiation schemes, and Noah land surface model.Then, the hourly mean meteorological data from 13:00 to 15:00 local time with a spatial resolution of 30 km was interpolated to a 0.1 × 0.1 • grid same as the NO 2 tropospheric column products.

Population Data
Worldwide gridded population data are available at 5-year intervals from 1995 to 2020 from the NASA Socioeconomic Data and Applications Center (Gridded Population of the World, v4; http://sedac.ciesin.columbia.edu/).The population data in 2013 was obtained by linearly-interpolating the data in 2010 and 2015 using 0.1 × 0.1 • resolution.

GTWR Model
The GTWR model for the relationship of ground NO 2 concentrations and satellite tropospheric columns can be expressed as [40]: where (u i , v i , t i ) represents the given coordinates of the training sample i in location (u i , v i ) at time t i .NO 2_ground(i) is the ground level NO 2 concentration observed by the ambient monitoring station at (u i , v i , t i ).NO 2_Trop(i) is the OMI NO 2 column density, β 0 (u i , v i , t i ) indicates the intercept of the GTWR model, β 1 (u i , v i , t i ) is a coefficient describing the unique spatial and temporal relationship between NO 2_ground(i) and NO 2_Trop(i) .ε i is the random error.
We introduced a number of meteorological parameters to the GTWR, i.e., air temperature at 2 m above the ground (T), relative humidity (RH), wind speed at 10 m above the ground (WS), planetary boundary layer height (PBLH), dew point temperature at 2 m above the ground (T d ), and the ambient pressure near ground (P).Akaike's information criterion (AIC) [41] was used to judge whether the GTWR performance could be improved with the addition of each specific meteorological parameter.The AIC value for the GTWR model is expressed as: where σ is the maximum likelihood estimation of the standard deviation for random error ε i (i = 1, 2, . . ., n). S is the hat matrix of the dependent variable.tr(S) is the trace of matrix S. S and σ are calculated using Equations ( 17) and (18), respectively.The smaller AIC is, the better the model performance will be.As indicated in Table 1, the model performance improves substantially when the meteorological parameters of PBLH, RH, WS, T, and P are included.This is because that: (1) high temperature can increase photochemical reactions and hence reduce the lifetime of NO 2 ; (2) high relative humidity is related to low NO 2 concentration since it enhances the conversion rate of secondary aerosol from NO X ; (3) high PBLH is often related to low NO 2 concentration when it is supposed that NO 2 are well-mixed and confined within the planetary boundary layer; (4) high wind speed is favorable to pollutant dispersion that will result in the decrease of NO 2 concentration; and (5) high pressure increases atmospheric stability, leading to less atmospheric general circulation and thus more NO 2 .The GTWR can be modified as: , and β 6 (u i , v i , t i ) denote the slope of T, RH, PBLH, WS, and P, respectively.In the GTWR model, a local weighted least squares algorithm is employed to determine the parameters of β(u i , v i , t i ): where W(u 0 , v 0 , t 0 ) is a square matrix comprising the geographically and temporally weighted values of training datasets for measurement i by the diagonal elements.X and Y are, respectively, expressed as: Remote Sens. 2017, 9, 950 The temporal distance d t i0 and the spatial distance d s i0 are given by: By combining the temporal distance d t i0 and the spatial distance d s i0 , the spatio-temporal distance is defined as: where ⊗ denotes different kinds of operators.Here, the "+" operator is adopted, the d st i0 is hence computed by: where λ and µ stand for the scale factors of temporal and spatial distance, respectively.Furthermore, an ellipsoidal coordinate system is used to calculate the d st i0 : Gaussian distance decay-based functions and Euclidean distance are chosen to construct the spatio-temporal weight matrix W(u 0 , v 0 , t 0 ).The diagonal element w i (u 0 , v 0 , t 0 ) of the W(u 0 , v 0 , t 0 ) can be obtained by: where h ST , h T and h S are the parameters of spatio-temporal, spatial, and temporal bandwidths, respectively.Adaptive spatio-temporal bandwidths are adopted according to the density of sample points around the given point (u 0 , v 0 , t 0 ).When many sample points are closely distributed around the given point, the bandwidths are small.On the contrary, if there are not enough sample points near it, the bandwidths are larger when THE sample points are sparsely distributed.In practice, the bandwidths are determined with an optimization technique by cross-validation through minimizing Equation (14).
where the function ŷi (h ST ) denotes the predicted value from the GTWR which is built without sample i.
The ground level NO 2 at (u i , v i , t i ) is estimated by: where , WS (i) , P (i) , and where S is the hat matrix of Y and is calculated as: The maximum likelihood estimation of the standard deviation for rand error is calculated as: where RSS is the residual sum of squares between estimated ground level NO 2 concentrations and observed ones:

Population-Weighted NO 2
The population data are introduced to calculate the population-weighted NO 2 (PNO 2 ) for different province-level administrative units: where PNO j 2 is the population-weighted NO 2 for province j, NO j,k 2 and Population j,k are the NO 2 concentration and population data of pixel k in province j respectively.

Implementation Process and Statistical Indicators
To correlate the ground-based measurements with satellite data, the 715 ambient monitoring stations in the central-eastern China were merged into 509 stations by averaging all of the measurements within a grid of 0.1 × 0.1 • .For the 509 grid cells, the total numbers of satellites and ambient monitoring observations are 54,867 (Figure 2a) and 110,545 (Figure 2b), respectively.Combining the satellite and ground observations, there are 31,463 valid data pairs (Figure 2c).The spatial distribution of the numbers of filtered satellite observations in Figure 2a shows a north-south difference, which is likely due to a higher cloud fraction over southern China.These 509 stations with total 31,463 dataset were divided randomly into 10 groups.The model fitting and cross-validation process was repeated 10 times, for every time one group was used for the cross-validation and the rest were used to train the fitting model until all groups were entered into the cross-validation once, thereby creating out-of-sample predictions for all the stations [42].To be more specific, all of the 31,463 datasets were used both in the fitting and the cross-validation.Some statistical indicators were employed to quantitatively assess the model performances.They are the coefficient of determination (R 2 ), whose higher value indicating better fitting accuracy, the root mean square error (RSME), that is sensitive to both systematic and random errors, the mean absolute difference (MAD), that measures the mean error magnitude, and the mean absolute percentage error (MAPE), which characterizes the prediction accuracy of a statistical model.with total 31,463 dataset were divided randomly into 10 groups.The model fitting and cross-validation process was repeated 10 times, for every time one group was used for the cross-validation and the rest were used to train the fitting model until all groups were entered into the cross-validation once, thereby creating out-of-sample predictions for all the stations [42].To be more specific, all of the 31,463 datasets were used both in the fitting and the cross-validation.
Some statistical indicators were employed to quantitatively assess the model performances.They are the coefficient of determination (R 2 ), whose higher value indicating better fitting accuracy, the root mean square error (RSME), that is sensitive to both systematic and random errors, the mean absolute difference (MAD), that measures the mean error magnitude, and the mean absolute percentage error (MAPE), which characterizes the prediction accuracy of a statistical model.

Spatio-Temporal Non-Stationarities between Tropospheric-Columnar and Ground Level NO2
According to previous studies [43,44], the tropospheric NO2 profiles show a large spatial-temporal variation.It is necessary to assess the impact of the spatio-temporal non-stationarities on the satellite-estimated ground level NO2 concentrations.Lamsal et al. [38] showed that OMI retrievals are underestimated in urban regions and overestimated in remote areas about 20%.To isolate the influence of different land covers types, 293 pure urban grids and corresponding ambient stations were picked out from the total 509 grids and stations.As shown in Figure 3a, the mean values of tropospheric-columnar and corresponding ground level NO2 over three provinces in eastern China (see also Figure 1) i.e., Shandong, Zhejiang, and Hunan, are compared.The mean values of OMI tropospheric NO2 columns of the three provinces are different when the column data is composited with respect to the ground level NO2 mass concentrations from ambient monitoring stations.This is related to the spatial difference in tropospheric NO2 profiles due to different topographies and meteorological conditions.Moreover, the mean values of OMI tropospheric NO2 columns in summer (May to July 2013), autumn (August to October 2013), winter (November 2013 to January 2014), and spring (February to April 2014) are compared in Figure 3b.The relationship between the NO2 columns and ground level NO2 shows a significant seasonal variation.The NO2 columns in winter and autumn are higher than those in summer and spring when the values of ground level NO2 are at the same level.This seasonal difference is more notable when ground concentrations increase, which is likely because of the longer lifetime of NO2 in winter and autumn as compared to that in summer and spring.Consequently, it can exist for a longer time in the upper layer in the case of high ground emissions.The numbers of satellite observations used in Figure 3a,b are given in Tables 2 and 3, respectively.It should be pointed out

Spatio-Temporal Non-Stationarities between Tropospheric-Columnar and Ground Level NO 2
According to previous studies [43,44], the tropospheric NO 2 profiles show a large spatial-temporal variation.It is necessary to assess the impact of the spatio-temporal non-stationarities on the satellite-estimated ground level NO 2 concentrations.Lamsal et al. [38] showed that OMI retrievals are underestimated in urban regions and overestimated in remote areas about 20%.To isolate the influence of different land covers types, 293 pure urban grids and corresponding ambient stations were picked out from the total 509 grids and stations.As shown in Figure 3a, the mean values of tropospheric-columnar and corresponding ground level NO 2 over three provinces in eastern China (see also Figure 1) i.e., Shandong, Zhejiang, and Hunan, are compared.The mean values of OMI tropospheric NO 2 columns of the three provinces are different when the column data is composited with respect to the ground level NO 2 mass concentrations from ambient monitoring stations.This is related to the spatial difference in tropospheric NO 2 profiles due to different topographies and meteorological conditions.Moreover, the mean values of OMI tropospheric NO 2 columns in summer (May to July 2013), autumn (August to October 2013), winter (November 2013 to January 2014), and spring (February to April 2014) are compared in Figure 3b.The relationship between the NO 2 columns and ground level NO 2 shows a significant seasonal variation.The NO 2 columns in winter and autumn are higher than those in summer and spring when the values of ground level NO 2 are at the same level.This seasonal difference is more notable when ground concentrations increase, which is likely because of the longer lifetime of NO 2 in winter and autumn as compared to that in summer and spring.Consequently, it can exist for a longer time in the upper layer in the case of high ground emissions.The numbers of satellite observations used in Figure 3a,b are given in Tables 2 and 3, respectively.It should be pointed out that the numbers of satellite observations for high ground level NO 2 (>100 µg/m 3 ) are less than five in Hunan and in summer.

Comparison between Model Fitted and Ground-Observed NO2
The ordinary least squares (OLS), GWR, temporally weighted regression (TWR), and GTWR models were tested using the same datasets.As shown in Tables 4 and 5, the OLS performance reveals that the tropospheric NO2 columns are potentially useful for ground level NO2 with R 2 of 0.45 and 0.44 for fitting and validation, respectively.The TWR outperforms the GWR with significant increases of R 2 values from 0.55 and 0.49 to 0.61 and 0.55.This suggests that the temporal non-stationarity is more dominant than the spatial non-stationarity between the tropospheric NO2 columns and ground level NO2.Among the four models, the GTWR has the best performance in both model-fitting and cross-validation with the highest R 2 and lowest errors (RMSE, MAD, and MAPE).Nevertheless, the GTWR regression shows a slight over-fitting, i.e., the R 2 generated from the cross-validation is 0.09 smaller than that from the model-fitting.In addition, the scatter plots in Figure 4 shows the largest correlation slope and the smallest intercept for the GTWR model.It is worth noting that all of the regression line slopes for the four models are less than 1. Figure 5 is present to assess the impact of the numbers of valid observations on the GTWR performance.The R 2 over Hunan (Figure 5a) is smaller than those over Shandong (Figure 5b) and Zhejiang (Figure 5c), due to less observations.

Comparison between Model Fitted and Ground-Observed NO 2
The ordinary least squares (OLS), GWR, temporally weighted regression (TWR), and GTWR models were tested using the same datasets.As shown in Tables 4 and 5, the OLS performance reveals that the tropospheric NO 2 columns are potentially useful for ground level NO 2 with R 2 of 0.45 and 0.44 for fitting and validation, respectively.The TWR outperforms the GWR with significant increases of R 2 values from 0.55 and 0.49 to 0.61 and 0.55.This suggests that the temporal non-stationarity is more dominant than the spatial non-stationarity between the tropospheric NO 2 columns and ground level NO 2 .Among the four models, the GTWR has the best performance in both model-fitting and cross-validation with the highest R 2 and lowest errors (RMSE, MAD, and MAPE).Nevertheless, the GTWR regression shows a slight over-fitting, i.e., the R 2 generated from the cross-validation is 0.09 smaller than that from the model-fitting.In addition, the scatter plots in Figure 4 shows the largest correlation slope and the smallest intercept for the GTWR model.It is worth noting that all of the regression line slopes for the four models are less than 1. Figure 5 is present to assess the impact of the numbers of valid observations on the GTWR performance.The R 2 over Hunan (Figure 5a) is smaller than those over Shandong (Figure 5b) and Zhejiang (Figure 5c), due to less observations.There are three possible errors in the estimation of ground level NO2 concentrations using the satellite-based GTWR model.First, satellite data are collected over an area of hundreds of km 2 , while in-situ measurements are point observations.Second, the errors in the retrieval of OMI tropospheric NO2 columns are underestimated in urban regions and overestimated in remote areas by about −20% and 20%, respectively [38].Third, the uncertainty in the meteorological parameters can affect the vertical distribution of tropospheric NO2.Zhang et al. [45] validated the NCEP FNL data against meteorological station data over Henan, China during 2012, and they found the errors of air temperature and pressure are −3~2 K and −10~10 hPa, respectively.We introduced the expected random errors (Gaussian distribution) from the tropospheric NO2 columns, air temperature, and air pressure, to assess their impact on the performance of the GTWR model.As shown in Table 6 and Figure 6, our model uncertainties are relatively low with the expected uncertainties from the model parameters.In the GTWR model, the smaller the spatio-temporal distance between two samples is, the greater weight coefficients are given.As illustrated in Figure 7, the GTWR performs much better than the OLS for the samples whose distances to the ambient monitoring stations are within 100 km, whereas the GTWR performance is worse than the OLS (0.41 versus 0.44 for R 2 ) for the samples that are more than 100 km away from the ambient monitoring stations.In the regions like Anhui, Jiangxi, and Fujian, where the ambient monitoring stations are very sparse and unevenly distributed, the nearest samples are mostly more than 100 km away.Hence, we used adjustable bandwidths according to the sample distances rather than the fixed ones.As compared to Figure 7, the R 2 in Figure There are three possible errors in the estimation of ground level NO 2 concentrations using the satellite-based GTWR model.First, satellite data are collected over an area of hundreds of km 2 , while in-situ measurements are point observations.Second, the errors in the retrieval of OMI tropospheric NO 2 columns are underestimated in urban regions and overestimated in remote areas by about −20% and 20%, respectively [38].Third, the uncertainty in the meteorological parameters can affect the vertical distribution of tropospheric NO 2 .Zhang et al. [45] validated the NCEP FNL data against meteorological station data over Henan, China during 2012, and they found the errors of air temperature and pressure are −3~2 K and −10~10 hPa, respectively.We introduced the expected random errors (Gaussian distribution) from the tropospheric NO 2 columns, air temperature, and air pressure, to assess their impact on the performance of the GTWR model.As shown in Table 6 and Figure 6, our model uncertainties are relatively low with the expected uncertainties from the model parameters.There are three possible errors in the estimation of ground level NO2 concentrations using the satellite-based GTWR model.First, satellite data are collected over an area of hundreds of km 2 , while in-situ measurements are point observations.Second, the errors in the retrieval of OMI tropospheric NO2 columns are underestimated in urban regions and overestimated in remote areas by about −20% and 20%, respectively [38].Third, the uncertainty in the meteorological parameters can affect the vertical distribution of tropospheric NO2.Zhang et al. [45] validated the NCEP FNL data against meteorological station data over Henan, China during 2012, and they found the errors of air temperature and pressure are −3~2 K and −10~10 hPa, respectively.We introduced the expected random errors (Gaussian distribution) from the tropospheric NO2 columns, air temperature, and air pressure, to assess their impact on the performance of the GTWR model.As shown in Table 6 and Figure 6, our model uncertainties are relatively low with the expected uncertainties from the model parameters.In the GTWR model, the smaller the spatio-temporal distance between two samples is, the greater weight coefficients are given.As illustrated in Figure 7, the GTWR performs much better than the OLS for the samples whose distances to the ambient monitoring stations are within 100 km, whereas the GTWR performance is worse than the OLS (0.41 versus 0.44 for R 2 ) for the samples that are more than 100 km away from the ambient monitoring stations.In the regions like Anhui, Jiangxi, and Fujian, where the ambient monitoring stations are very sparse and unevenly distributed, the nearest samples are mostly more than 100 km away.Hence, we used adjustable bandwidths according to the sample distances rather than the fixed ones.As compared to Figure 7, the R 2 in Figure In the GTWR model, the smaller the spatio-temporal distance between two samples is, the greater weight coefficients are given.As illustrated in Figure 7, the GTWR performs much better than the OLS for the samples whose distances to the ambient monitoring stations are within 100 km, whereas the GTWR performance is worse than the OLS (0.41 versus 0.44 for R 2 ) for the samples that are more than 100 km away from the ambient monitoring stations.In the regions like Anhui, Jiangxi, and Fujian, where the ambient monitoring stations are very sparse and unevenly distributed, the nearest samples are mostly more than 100 km away.Hence, we used adjustable bandwidths according to the sample distances rather than the fixed ones.As compared to Figure 7, the R 2 in Figure 8 improves from 0.41 to 0.50 for the samples with larger distances (>100 km) after adjusting the bandwidth.8 improves from 0.41 to 0.50 for the samples with larger distances (>100 km) after adjusting the bandwidth.To compare the GTWR method with the chemistry transport model (CTM) approach, we generated the tropospheric NO2 profiles by using a WRF-Chem model with the monthly MIX Asian    To compare the GTWR method with the chemistry transport model (CTM) approach, we generated the tropospheric NO2 profiles by using a WRF-Chem model with the monthly MIX Asian To compare the GTWR method with the chemistry transport model (CTM) approach, we generated the tropospheric NO 2 profiles by using a WRF-Chem model with the monthly MIX Asian anthropogenic emission inventory [46].This emission inventory has a spatial resolution of 0.25 × 0.25 • and involves four emission categories, including industry, power, transport, and residential.The model has 20 vertical levels and the top level pressure is 200 hPa.The RADM2 chemical mechanism is used for the gas-phase chemical reaction calculations.The Modal Aerosol Dynamics Model for Europe-MADE/SORGAM is chosen for the aerosol scheme.Then, we estimated the ground level NO 2 concentrations in January 2014 over central-eastern China using the approach described by Lamsal et al. [18,19].As shown in Figure 9, the result of the GTWR fitted is much better than the WRF-Chem.Recently, Gu et al. [43] estimated the ground level NO 2 over China using the chemistry transport model approach with the Community Multi-scale Air Quality (CMAQ) model by considering the influence of China's high atmospheric pollution on obtaining the vertical distribution of tropospheric NO 2 profiles.They achieved a correlation coefficient (R) of 0.80 for January 2014, which is comparable to the coefficient of determination (R 2 ) of 0.60 obtained by the GTWR.
Remote Sens. 2017, 9, 950 13 of 19 anthropogenic emission inventory [46].This emission inventory has a spatial resolution of 0.25 × 0.25° and involves four emission categories, including industry, power, transport, and residential.The model has 20 vertical levels and the top level pressure is 200 hPa.The RADM2 chemical mechanism is used for the gas-phase chemical reaction calculations.The Modal Aerosol Dynamics Model for Europe-MADE/SORGAM is chosen for the aerosol scheme.Then, we estimated the ground level NO2 concentrations in January 2014 over central-eastern China using the approach described by Lamsal et al. [18,19].As shown in Figure 9, the result of the GTWR fitted is much better than the WRF-Chem.Recently, Gu et al. [43] estimated the ground level NO2 over China using the chemistry transport model approach with the Community Multi-scale Air Quality (CMAQ) model by considering the influence of China's high atmospheric pollution on obtaining the vertical distribution of tropospheric NO2 profiles.They achieved a correlation coefficient (R) of 0.80 for January 2014, which is comparable to the coefficient of determination (R 2 ) of 0.60 obtained by the GTWR.To further evaluate the performance of the GTWR model, the comparison of the annual mean of NO2 concentrations between the model-fitted and ground-observed data is given in Figure 10.Overall, the NO2 concentrations estimated by the GTWR model agree well with the ground-based measurements.More than 90% of the cross-validation stations possess low mean discrepancies of less than 10 μg/m 3 .

Spatial Distribution of GTWR Fitted Ground-Observed NO2
The spatial distributions of annual mean NO2 values are shown in Figure 11.The fitted ground-observed NO2 concentrations by GTWR in (a) have similar spatial patterns to the satellite tropospheric NO2 columns in (b).The concentrations are comparable to the interpolated in situ observations using the Kriging method in (c) over the region with high values.Importantly, in the areas without monitoring stations (e.g., southern Jiangxi and northern Fujian), Figure 11a provides more reasonable estimations that are overestimated in Figure 11c.In Figure 11a, high NO2 concentrations are clustered in the regions of North China Plain, Yangtze River Delta, and Pearl River Delta.Especially, the NO2 concentrations in southern Hebei, northern Henan, central Shandong, and southern Shaanxi exceeded the Level 2 standard of the Chinese National Ambient Air Quality Standard (40 μg/m 3 ).Figure 12 denotes dramatic seasonal changes in the spatial distribution of GTWR fitted ground level NO2.Unparalleled high values are found in winter, while the lowest values are found in summer.To further evaluate the performance of the GTWR model, the comparison of the annual mean of NO 2 concentrations between the model-fitted and ground-observed data is given in Figure 10.Overall, the NO 2 concentrations estimated by the GTWR model agree well with the ground-based measurements.More than 90% of the cross-validation stations possess low mean discrepancies of less than 10 µg/m 3 .

Spatial Distribution of GTWR Fitted Ground-Observed NO 2
The spatial distributions of annual mean NO 2 values are shown in Figure 11.The fitted ground-observed NO 2 concentrations by GTWR in (a) have similar spatial patterns to the satellite tropospheric NO 2 columns in (b).The concentrations are comparable to the interpolated in situ observations using the Kriging method in (c) over the region with high values.Importantly, in the areas without monitoring stations (e.g., southern Jiangxi and northern Fujian), Figure 11a provides more reasonable estimations that are overestimated in Figure 11c.In Figure 11a, high NO 2 concentrations are clustered in the regions of North China Plain, Yangtze River Delta, and Pearl River Delta.Especially, the NO 2 concentrations in southern Hebei, northern Henan, central Shandong, and southern Shaanxi exceeded the Level 2 standard of the Chinese National Ambient Air Quality Standard (40 µg/m 3 ).Figure 12 denotes dramatic seasonal changes in the spatial distribution of GTWR fitted ground level NO 2 .Unparalleled high values are found in winter, while the lowest values are found in summer.

Population-Weighted Ground Level NO 2 Concentrations
Given the NO 2 toxicity to human health, it is necessary to evaluate population exposure levels over different provinces.Traditionally, the province-level mean NO 2 concentration provided by the Chinese environmental protection agencies are the arithmetic means of all values in administered cities.Here, we calculated the annual mean population-weighted NO 2 (AMPNO 2 ) concentrations by using Equation (20).The annual mean NO 2 concentrations (AMNO 2 ) and AMPNO 2 of 17 provinces in central-eastern China are summarized in Table 7. AMPNO2 is higher than AMNO2 for all of the provinces, especially for densely populated provinces, e.g., Hebei, Beijing, and Guangdong.People from these 17 provinces except Anhui, Fujian, Jiangxi, and Hunan, are exposed to high-level NO 2 concentrations (>30 µg/m 3 ).Heibei, Tianjin, and Beijing suffer from the most serious NO 2 pollution, with more than 70% of people affected by high-level NO 2 .The satellite-estimated ground level NO 2 concentration is observed at afternoon (13:00-15:00) leading to underestimated annual mean values.

Conclusions
In this study, a satellite-based GTWR model has been applied to estimate ground level NO 2 concentrations over central-eastern China.OMI tropospheric NO 2 columns, together with ambient monitoring station measurements and meteorological data from May 2013 to April 2014 were considered.The results show that the GTWR model produces the highest cross-validation R 2 (0.60) and the lowest errors (RMSE, MAD, and MAPE), in comparison with other models, i.e., OLS, GWR, and TWR.The model performance is significantly correlated with the meteorological parameters that likely describe the NO 2 vertical profile shapes.Our method is better than or comparable to the CTM method.
The satellite-estimated spatial distribution of annual mean NO 2 shows a similar spatial pattern to the tropospheric NO 2 column and possesses similar value with the in situ observation.High annual mean NO 2 concentrations (>40 µg/m 3 ) are found in southern Hebei, northern Henan, central Shandong, and southern Shaanxi.Seasonal changes in the spatial distribution of ground level NO 2 are easily identifiable with unparalleled high values in winter and the lowest values in summer.The population-weighted NO 2 demonstrates that people who lived in densely populated areas are more likely to be exposed to high NO 2 pollution.
One of the major error sources in the estimation of ground level NO 2 concentrations using OMI data is the spatial gradient and the horizontal inhomogeneity between individual satellite pixels.In September 2017, the TROPOMI/S5P will be launched and measure tropospheric NO 2 columns with a higher spatial resolution (7 km × 7 km) [47], which enables an improved accuracy in the ground level NO 2 estimation.

Figure 1 .
Figure 1.Study area and locations of ambient monitoring stations.

Figure 1 .
Figure 1.Study area and locations of ambient monitoring stations.

Figure 2 .
Figure 2. Numbers of satellite observations (a), ambient monitoring observations (b), and valid satellite and ground observation pairs (c) for 593 grid cells.

Figure 2 .
Figure 2. Numbers of satellite observations (a), ambient monitoring observations (b), and valid satellite and ground observation pairs (c) for 593 grid cells.
Remote Sens. 2017, 9, 950 9 of 19 that the numbers of satellite observations for high ground level NO2 (>100 μg/m 3 ) are less than five in Hunan and in summer.

Figure 3 .
Figure 3. Mean values of tropospheric NO2 columns (×10 16 molec/cm 2 ) in different provinces (a) and different seasons (b) when the column data is composited with respect to the ground level NO2 mass concentrations observed at 293 pure urban ambient monitoring stations.Error bars stand for one standard deviation

Figure 3 .
Figure 3. Mean values of tropospheric NO 2 columns (×10 16 molec/cm 2 ) in different provinces (a) and different seasons (b) when the column data is composited with respect to the ground level NO 2 mass concentrations observed at 293 pure urban ambient monitoring stations.Error bars stand for one standard deviation.

Figure 4 .
Figure 4. Scatter plots between the observed NO2 and predicted NO2 concentrations using OLS (a), GWR (b), TWR (c), and GTWR (d) for cross validation over central-eastern China from May 2013 to April 2014.

Figure 4 .
Figure 4. Scatter plots between the observed NO 2 and predicted NO 2 concentrations using OLS (a), GWR (b), TWR (c), and GTWR (d) for cross validation over central-eastern China from May 2013 to April 2014.

Figure 5 .
Figure 5. Scatter plots between the observed NO2 and predicted NO2 concentrations for cross validation over Shandong (a), Zhejiang (b), and Hunan (c) from May 2013 to April 2014.

Figure 5 .
Figure 5. Scatter plots between the observed NO 2 and predicted NO 2 concentrations for cross validation over Shandong (a), Zhejiang (b), and Hunan (c) from May 2013 to April 2014.

Figure 5 .
Figure 5. Scatter plots between the observed NO2 and predicted NO2 concentrations for cross validation over Shandong (a), Zhejiang (b), and Hunan (c) from May 2013 to April 2014.

Figure 6 .
Figure 6.Scatter plots between the observed NO 2 and GTWR predicted NO2 concentrations with random errors from tropospheric NO 2 columns (20%) (a), air temperature (2 K) (b), and air pressure (10 hPa) (c) over central-eastern China from May 2013 to April 2014.

Figure 8 .
Figure 8. Scatter plots between the observed NO2 and GTWR predicted NO2 concentrations for samples with different distances ((a) <20 km; (b) 20-100 km; (c) >100 km) after adjusting the bandwidth for cross validation over central-eastern China from May 2013 to April 2014.

Figure 7 .
Figure 7. Scatter plots between the observed NO 2 and predicted NO 2 concentrations using GTWR ((a-c) and OLS (d-f) for samples with different distances (a,d) <20 km; (b,e) 20-100 km; (c,f) >100 km) for cross validation over central-eastern China from May 2013 to April 2014.

Figure 8 .
Figure 8. Scatter plots between the observed NO2 and GTWR predicted NO2 concentrations for samples with different distances ((a) <20 km; (b) 20-100 km; (c) >100 km) after adjusting the bandwidth for cross validation over central-eastern China from May 2013 to April 2014.

Figure 8 .
Figure 8. Scatter plots between the observed NO 2 and GTWR predicted NO 2 concentrations for samples with different distances ((a) <20 km; (b) 20-100 km; (c) >100 km) after adjusting the bandwidth for cross validation over central-eastern China from May 2013 to April 2014.

Figure 9 .
Figure 9. Scatter plots between the observed NO2 and predicted NO2 concentrations by GTWR (a) and WRF-Chem (b) for cross validation over central-eastern China in January 2014.

Figure 9 .
Figure 9. Scatter plots between the observed NO 2 and predicted NO 2 concentrations by GTWR (a) and WRF-Chem (b) for cross validation over central-eastern China in January 2014.

Figure 10 .
Figure 10.Annual means of observed NO2 and fitted NO2 by the GTWR model for all 509 stations.Figure 10.Annual means of observed NO 2 and fitted NO 2 by the GTWR model for all 509 stations.

Figure 10 .
Figure 10.Annual means of observed NO2 and fitted NO2 by the GTWR model for all 509 stations.Figure 10.Annual means of observed NO 2 and fitted NO 2 by the GTWR model for all 509 stations.

Figure 11 .
Figure 11.Spatial distribution of annual mean values of (a) ground level NO2 concentrations fitted by GTWR model, (b) NO2 tropospheric columns, and (c) ground level NO2 concentrations interpolated by Kriging method.Figure 11.Spatial distribution of annual mean values of (a) ground level NO 2 concentrations fitted by GTWR model, (b) NO 2 tropospheric columns, and (c) ground level NO 2 concentrations interpolated by Kriging method.

Figure 11 .
Figure 11.Spatial distribution of annual mean values of (a) ground level NO2 concentrations fitted by GTWR model, (b) NO2 tropospheric columns, and (c) ground level NO2 concentrations interpolated by Kriging method.Figure 11.Spatial distribution of annual mean values of (a) ground level NO 2 concentrations fitted by GTWR model, (b) NO 2 tropospheric columns, and (c) ground level NO 2 concentrations interpolated by Kriging method.

Figure 12 .Figure 12 .
Figure 12.Spatial distribution of seasonal mean values of GTWR fitted (upper) and Kriging interpolated (lower) ground level NO2 concentrations.Figure 12.Spatial distribution of seasonal mean values of GTWR fitted (upper) and Kriging interpolated (lower) ground level NO 2 concentrations.

Table 1 .
Akaike 's information criterion (AIC) values when satellite, planetary boundary layer height (PBLH), relative humidity (RH), wind speed at 10 m above the ground (WS), air temperature at 2 m above the ground (T), and ambient pressure near ground (P) data are included respectively in the geographically and temporally weighted regression (GTWR) model.

Table 2 .
Numbers of satellite observations used in Figure3a.

Table 3 .
Numbers of satellite observations used in Figure3b.

Table 2 .
Numbers of satellite observations used in Figure3a.

Table 3 .
Numbers of satellite observations used in Figure3b.

Table 6 .
Quantitative assessment of GTWR cross-validation with random errors from tropospheric NO 2 columns, air temperature, and air pressure.

Table 6 .
Quantitative assessment of GTWR cross-validation with random errors from tropospheric NO2 columns, air temperature, and air pressure.

Table 6 .
Quantitative assessment of GTWR cross-validation with random errors from tropospheric NO2 columns, air temperature, and air pressure.