《遥感地学分析》课程教学资源（文献资料）Satellitebased mapping ofdaily highresolution ground PM 2.5 in China viaspacetime regression modeling

Remote Sensing of Environment 206 (2018)72-83 Contents lists available at ScienceDirect Remote Sensing of Environment ELSEVIER journal homepage:www.elsevier.com/locate/rse Satellite-based mapping of daily high-resolution ground PM2s in China via space-time regression modeling Qingqing He,Bo Huang . ARTICLE INFO ABSTRACT The of s bility to extend the 1.Introductio studies CK 2008:M dustrialization of Chir has led to PMas be ing a domi 014Ch 01:Guo et al Co e he th effects of PM-ex osure across Chin using PM and Re nt,The Chinese University of Hong Kong.Hong Kong ps/,0w￥10.10162017.12018 000ver Inc.gt reserved December 2017;Accepted 15 December 2017

Contents lists available at ScienceDirect Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse Satellite-based mapping of daily high-resolution ground PM2.5 in China via space-time regression modeling Qingqing Hea,b , Bo Huanga,b,c,⁎,1 aDepartment of Geography and Resource Management, The Chinese University of Hong Kong, Hong Kong b Big Data Decision Analytics (BDDA) Research Centre, The Chinese University of Hong Kong, Hong Kong c Institute of Space and Earth Information Science, The Chinese University of Hong Kong, Hong Kong ARTICLE INFO Keywords: Aerosol optical depth Geographically and temporally weighted regression (GTWR) Interior point algorithm (IPA) PM2.5 China ABSTRACT The use of satellite-retrieved aerosol optical depth (AOD) data and statistical modeling provides a promising approach to estimating PM2.5 concentrations for a large region. However, few studies have conducted high spatial resolution assessments of ground-level PM2.5 for China at the national scale, due to the limitations of high-resolution AOD products. To generate high-resolution PM2.5 for the entirety of mainland China, a combined 3 km AOD dataset was produced by blending the newly released 3 km-resolution Moderate Resolution Imaging Spectroradiometer (MODIS) Dark Target AOD data with the 10 km-resolution MODIS Deep Blue AOD data. Using this dataset, surface PM2.5 measurements, and ancillary information, a space-time regression model that is an improved geographically and temporally weighted regression (GTWR) with an interior point algorithm (IPA)- based efficient mechanism for selecting optimal parameter values, was developed to estimate a large set of daily PM2.5 concentrations. Comparisons with the popular spatiotemporal models (daily geographically weighted regression and two-stage model) indicated that the proposed GTWR model, with an R2 of 0.80 in cross-validation (CV), performs notably better than the two other models, which have an R2 in CV of 0.71 and 0.72, respectively. The use of the combined 3-km high resolution AOD data was found not only to present detailed particle gra￾dients, but also to enhance model performance (CV R2 is only 0.32 for the non-combined AOD data). Furthermore, the GTWR's ability to predict PM2.5 for days without PM2.5-AOD paired samples and to generate historical PM2.5 estimates was demonstrated. As a result, fine-scale PM2.5 maps over China were generated, and several PM2.5 hotspots were identified. Therefore, it becomes possible to generate daily high-resolution PM2.5 estimates over a large area using GTWR, due to its synergy of spatial and temporal dimensions in the data and its ability to extend the temporal coverage of PM2.5 observations. 1. Introduction Fine particulate matter, a complex of solid and liquid particles suspended in the air with aerodynamic diameters of 2.5 μm or less (PM2.5), has been associated with adverse effects on public health by many epidemiological studies (Kampa and Castanas, 2008; Madrigano et al., 2013; Pope and Dockery, 2006). The rapid urbanization and in￾dustrialization of China has led to PM2.5 becoming a dominant factor in air pollution, especially over urban areas, and thus an unprecedented issue of public concern (Bi et al., 2014; Che et al., 2014; Guo et al., 2009). Given the advent of severe PM2.5 levels in China, it is urgent to assess the health effects of PM2.5 exposure across China, using PM2.5 data with a high spatial resolution. However, currently, such an as￾sessment is seriously hindered by the limited measurements available from the sparse network of surface monitoring stations. Because satellite remote sensing has the capacity to provide data with large (even global) spatial coverage, satellite-derived aerosol op￾tical depth (AOD) information is an alternative method of inferring ground-level PM2.5 concentrations (Engel-Cox et al., 2004; Wang and Christopher, 2003). To better investigate PM2.5 exposure for air pollu￾tion assessment and epidemiological research, satellite-derived PM2.5 concentrations at high spatial resolution are needed. High-resolution studies in this research field have been increasingly carried out in North America using 1-km Multiangle Implementation of Atmospheric Cor￾rection (MAIAC) AOD (Hu et al., 2014; Kloog et al., 2014). However, equivalent studies for China are scarce, due to the lack of high-resolu￾tion satellite-derived AOD products. Studies on a national scale have only presented PM2.5 variations with resolutions ranging from 10 to https://doi.org/10.1016/j.rse.2017.12.018 Received 13 March 2017; Received in revised form 2 December 2017; Accepted 15 December 2017 ⁎ Corresponding author at: Department of Geography and Resource Management, The Chinese University of Hong Kong, Hong Kong. 1 Co-first author. E-mail address: bohuang@cuhk.edu.hk (B. Huang). Remote Sensing of Environment 206 (2018) 72–83 0034-4257/ © 2017 Elsevier Inc. All rights reserved. T

Q.He,B.Huang ote Sensing of cnt206(2018)72-83 50 km,employing widely used AOD products,e.g.,the Moderate R entire of mainland China,in the setting of a sp sio (MODIS)Level 2 el and Mult d m mprove d t relatively high spa ion(Ren 1.2013).The fine red compreher daily GWR nd and the exte ned a rel ween sun derives AOD o dark which may 2.Data collection and p reproces land s other ds.it is 2.1.Ground-level PMs observations rety of Chi MODIS 3 AOD 015 verage of the 3km MODIS AOD forg r2015 d f n the china ep (h .th nd chin d by th and sur ting 05.201 15) as N tAir Quality St r CNAAQS)(C 011:K PM2s in Ambie Air Method,availab le at ng et al. 2014:Zou et ).To date two the ee146emgesa329c he PM AOD rela ship.One is the daily GWR ad have little cov age 2.2.MODIS AOD date ih ti 2.2.1.MODIS 3 km DT A ral MODI is a se bee non board NASA Terra since 1999 od is 89103 mTera)and1:30p this model calib the and k surfa es (ede ted areasat ng the ts with am 04.31 x is O for introd MOD 10 n a 0p ol prod ol retrieval ds to g into ther 20×2 mber of PM -AOD paired sar ples ed out and di arded (Re er et aL 2013).MODIS 3km AOD da eogra nm d_ocean) ed b AOD 7).Unlike the ly discussed thi 12 m the estin fron days into the MODIS DB algo s to s rate AOD retrievals ove s,the er of GTWR for e lucidatin the re 2013.0 ke the T retrieval proce the DB algorithm firs et to be exr eter values for WR for a large sbredPBerosol is for 2015 usimg D AOD values were missing with data obtained. samp 2.2.3.MODIS AOD combin of this dicting spatial

50 km, employing widely used AOD products, e.g., the Moderate Re￾solution Imaging Spectroradiometer (MODIS) Level 2 10 km AOD product and Multi-angle Imaging SpectroRadiometer (MISR) Level 2 17.6 km AOD product (Fang et al., 2016; Ma et al., 2014; Ma et al., 2016). In early 2014, a 3-km aerosol product was included in the new MODIS Collection 6 (C6). It is the first daily global aerosol dataset with a relatively high spatial resolution (Remer et al., 2013). The fine aerosol gradients can benefit air quality applications, i.e., improve the perfor￾mance of statistical models and provide details of spatial variation for fine particulate matter (Xie et al., 2015). However, the MODIS 3 km AOD data were retrieved using the Dark Target (DT) algorithm that derives AOD over dark surfaces, which may cause a large number of missing values in the daily AOD image, especially for a large area with a complex land surface (He et al., 2017). In other words, it is impossible to generate a daily PM2.5 prediction at high spatial resolution for the entirety of China using only the MODIS 3 km AOD product as the pri￾mary predictor (Li et al., 2017; You et al., 2016). Thus, improving the daily coverage of the 3 km MODIS AOD for ground PM2.5 estimation over a large area with multiple types of land cover is essential. Various models have been developed to elucidate the quantitative association between satellite-based AOD and surface PM2.5, ranging from a simple linear regression model (Engel-Cox et al., 2004; Wang and Christopher, 2003) to advanced statistical models such as the linear mixed effect (LME) model (Lee et al., 2011; Xie et al., 2015), land use regression model (Kloog et al., 2011; Kloog et al., 2012), and geo￾graphically weighted regression (GWR) model (Hu et al., 2013; Ma et al., 2014; Song et al., 2014; Zou et al., 2016). To date, two main spatial and temporal regression methods have been used to establish the PM2.5-AOD relationship. One is the daily GWR model, which ad￾dresses the spatial variability between PM2.5 and AOD for individual days and achieves reasonable results (Ma et al., 2014; Song et al., 2014). However, this model ignores the fact that the PM2.5-AOD re￾lationship could vary with time and may be dependent on previous days. Thus, daily GWR cannot make use of temporal autocorrelation existing in the data and fails to build a model for days with fewer PM2.5- AOD samples. The other method is a two-stage model that can account for both spatial and temporal non-stationarities (Hu et al., 2014; Wu et al., 2016). However, this model calibrates the spatial and temporal PM2.5-AOD relationship separately, i.e., using the LME model to ac￾count for its temporal variations in the first stage, and then the GWR model its spatial variations in the second stage. Consequently, this model does not simultaneously treat spatiotemporal variations of the PM2.5-AOD relationship; it is thus not in a better position to predict PM2.5 for days with insufficient or no PM2.5 observations. As such, it becomes difficult to use the two methods to generate high-resolution PM2.5 with high accuracy for a large area when there is a limited number of PM2.5-AOD paired samples. In recent years, the geographically and temporally weighted re￾gression (GTWR) model proposed by Huang et al. (2010), which in￾tegrates spatial and temporal distance, has gained increasing popularity in environmental studies (Bai et al., 2016; Chu et al., 2015; Guo et al., 2017). Unlike the two spatiotemporal models already discussed, this space-time regression model can simultaneously incorporate temporal information from the estimation day or from previous days into the spatial variability through a spatial-temporal weighting mechanism. Nevertheless, the predictive power of GTWR for elucidating the re￾lationship between PM2.5 and AOD, especially for days without sam￾ples, has yet to be explored. Optimization of parameter values for GTWR is also required to reduce the computational cost for a large dataset. As a corollary, there are still considerable challenges con￾fronting a space-time model, such as GTWR, in exploring the PM2.5- AOD relationship, especially over an extensive area with a large but insufficient sample dataset. The objective of this study is to overcome such challenges by pre￾dicting ground PM2.5 at high spatial resolution on a daily basis for the entirety of mainland China, in the setting of a space-time regression model using the MODIS 3 km AOD as the major predictor, and me￾teorological and land use data as ancillary variables. A customized approach to filling in missing AOD values was devised to improve the availability of daily MODIS 3 km AOD. Subsequently, the GTWR model improved with the parameter optimization approach was developed to generate ground PM2.5 concentrations at high resolution, and its in￾ferring power also explored comprehensively. A daily GWR model and a two-stage (LME + GWR) model were also implemented as benchmarks to examine the extent to which GTWR's simultaneous spatial and tem￾poral weighting established a reliable association between surface PM2.5 and satellite AOD. 2. Data collection and preprocessing 2.1. Ground-level PM2.5 observations Daily average PM2.5 observations from 1 January 2015 to 31 December 2015 were obtained from the China Environment Monitoring Center (http://106.37.208.233:20035/). Ground PM2.5 mass con￾centrations of mainland China are measured by the tapered-element oscillating microbalance or beta-attenuation method, with calibration and quality control according to national standards GB3095-2012 (China's National Ambient Air Quality Standard, or CNAAQS) (China, 2012) and HJ 618-2011 (Determination of Atmospheric Articles PM10 and PM2.5 in Ambient Air by Gravimetric Method, available at http:// english.sepa.gov.cn/Resources/standards/Air_Environment/). At pre￾sent, the observations come from 1456 monitoring stations in 329 cities (Fig. 1). Most are located in southeastern China, whereas the north￾western areas have little coverage. 2.2. MODIS AOD data 2.2.1. MODIS 3 km DT AOD products MODIS is a sensor that has been on board NASA Terra since 1999 and Aqua since 2002, providing two daily observations of columnar aerosol properties at approximately 10:30 am (Terra) and 1:30 pm (Aqua) local time. The DT algorithm was originally developed to derive aerosol properties over dark surfaces (e.g., dense vegetated areas) at a 10 km spatial resolution. To satisfy an identified need for monitoring fine-resolution pollution, aerosol datasets with a 3 km resolution (MxD04_3K, x is O for Terra and Y for Aqua) were introduced in the most recently released MODIS C6 in addition to the standard 10 km Level 2 aerosol products (MxD04_L2). The 3 km DT aerosol retrievals use a similar protocol to the 10 km DT aerosol products, but organize 6 × 6 pixels into a retrieval box for inversion, rather than 20 × 20, after all undesirable ones, i.e., cloud and snow/ice pixels, have been screened out and discarded (Remer et al., 2013). MODIS 3 km AOD data at 550 nm (scientific dataset name: Optical_Depth_Land_And_Ocean) in 2015 (downloaded from https://ladsweb.nascom.nasa.gov/) were used to model the PM2.5-AOD relationship. Data obtained before 2015 were used for calibration. 2.2.2. MODIS 10 km Deep Blue (DB) AOD products The MODIS DB algorithm aims to generate AOD retrievals over bright land (e.g., deserts) to fill the gaps left by the DT algorithm (Hsu et al., 2013). Unlike the DT retrieval process, the DB algorithm first produces aerosol properties at a 1 km resolution and then averages the individual retrievals into a 10 × 10 km grid. As with the 3 km AOD data, we calibrated DB aerosol retrievals for 2015 using DB AOD pro￾ducts obtained before 2015, and filled in some gaps where the 3 km DT AOD values were missing with data obtained in 2015. 2.2.3. MODIS AOD combination In MODIS Collection 6, the 3 km AOD data are retrieved using the DT algorithm, and the daily availability of 3 km AOD in the desert-like Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 73

Q.He.B.Huang sing of Emnt cmt206(2018)72-83 90*E 10E N PM2.5 (ug/m Elevation (m Fig.1.Ge is ow (even n Th for the four data product at a 10k gesolhution afotreD8 (3) ng values forAqua and 3km and Terra-3 o p the km AOD for -3 km va sa)The e p ze of the 3km DT AOD,the 10km DB AOD for Ao and Terra (SL S2). 2 Terra apcd-DB Aop da ed the Fig.S (Supplem R kmandsampld-D s to average e the Aqua and Terre values (R we first av the gaps b n DT and DB AOD retriev we calibra dicted values)fo applied th ere missine.That is the averaged DB. npled aod values erial to fill t lidate mthe 3 kn d on the ch pixel increased sig es (Li et al 009 et a lity.we fitted linear egres for the Aqu3 m DT.Terrakm DT.Aqua rempled-DB.and

area of northwestern China is very low (even no retrievals at all for many days). The MODIS DB algorithm was originally developed to obtain aerosol loading for bright targets, but it only provides the AOD product at a 10 km resolution. To address this issue, a four-step cus￾tomized approach of combining the MODIS 3 km DT and 10 km DB AOD products was developed to improve the daily coverage of the MODIS 3 km DT AOD, as follows: (1) Re-gridding the 10 km DB AOD To match the grid cell size of the 3 km DT AOD, the 10 km DB AOD were resampled to a 3-km spatial resolution grid using the cubic convolution resampling technique in ENVI/IDL5.1. (2) Calibrating the MODIS DT/DB AOD According to the global validation of MODIS aerosol products and Fig. S2 (Supplementary information [SI]), there are different re￾trieval accuracies for DT and DB AODs and thus systematic bias between the DT and DB AOD values (Remer et al., 2013; Sayer et al., 2014). Therefore, to optimize the daily AOD values and fill the gaps between DT and DB AOD retrievals, we calibrated AOD retrievals for 2015 using ground-level aerosol observations from the Aerosol Robotic Network (AERONET) from 2001 to 2014 (Ma et al., 2014). Comparison against AERONET has been widely used to va￾lidate satellite-based AOD (Sayer et al., 2014; Tao et al., 2015). A linear analysis was conducted based on the long-term data before 2015 to define the relationship between AERONET and satellite AOD. According to previous studies (Li et al., 2009; Remer et al., 2013), in which the relationship between AERONET and MODIS AODs showed strong seasonality, we fitted linear regressions for each season. These relationships were also separately established for the Aqua 3 km DT, Terra 3 km DT, Aqua resampled-DB, and Terra resampled-DB AOD datasets. These fitted parameters were then used to calibrate AOD values for the four datasets during 2015 (SI, S1). (3) Filling in missing values for Aqua and Terra A simple linear regression was performed for the matched grid cells between Aqua-3 km and Terra-3 km AOD for each day to predict the missing values (i.e., estimation of the Aqua-3 km AOD using the available Terra-3 km value, and vice versa). The same protocol was applied to the resampled-DB AODs to predict the missing AOD values for Aqua and Terra (SI, S2). (4) Combining the four datasets To improve the daily AOD coverage, Aqua/Terra 3 km DT AOD and Aqua/Terra resampled-DB AOD data were merged together for spatiotemporal modeling. We used the same method as with the 3 km and resampled-DB AODs to average the Aqua and Terra va￾lues. Because the 3 km DT has higher spatial resolution than the 10 km DB to present finer aerosol gradients, we first averaged the Aqua- and Terra-3 km DT AOD values (both retrievals and pre￾dicted values) for each day and then applied the averaged DB data resampled to a 3-km grid to fill pixels when the 3 km AOD values were missing. That is, the averaged DB-resampled AOD values were applied as supplementary material to fill the gap when the 3 km AOD values for both satellites were missing. After the combination, the number of daily AOD values for each pixel increased sig￾nificantly (see Section 5.2 for more details). When compared to the AERONET AOD observations from 2015, our combined AOD (fraction of retrievals falling within the expected error bounds (F) = 56.25% and RMSE = 0.20 for Aqua calibrated and pre￾dicted AOD, and F = 55.07% and RMSE = 0.19 for Terra) shows a Fig. 1. Geographical distribution of PM2.5 monitoring stations included in the study. Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 74

Q.H Hng Rmate5arg时ntommct26(2018y72- 。aeaR-A and day and are the location-t and and T,average the day of th d cell in whic 2.3.Auxiliary data oestimate the intercept of )and the sopes of ofs lei in term variables tion 44 1 2.4.Data integration andgnd e PM2 0wr3sd AOD PM2. obse other varis bles w problen Variance inflation 季 ological station data for each day we WS,T,P,PBL DEM).The resul ariabl ndkept the variables tofit the final GTWR me (the equation is given in Eq.() 3.1.3 eter selec 3.Dev and validation error agains ting dati this relati me e.The mple a large am com a and p is sons 3.1.GTWR model development edal an hich when arge tational time method to about 1day ().In addition,to reduce the mode ntting probler dth regime rather than adaptiv PMg=民。，.+月4，）×A0D+月4，》×NDM (SL.S4 +Bu,h,)×H+B,u.h.t)×Ws+民(4，.)×T ntry dth of 120 (,i)x PBL+B()x DEM+ 1) ern C ina for s me days due to

better accuracy than the MODIS operational 3-km DT AOD (F = 46.25%, RMSE = 0.38 for Aqua and F = 49.28%, RMSE = 0.36 for Terra), indicating the effectiveness of our proposed four-step ap￾proach. The detailed validation process is presented in S2 (SI). 2.3. Auxiliary data A number of meteorological and land-cover-related variables were used in this study. The meteorological variables were ground-based relative humidity (RH), temperature (T), wind speed (WS), pressure (P), and reanalysis-produced planetary boundary layer (PBL) data; the land￾cover-related variables are terrain (DEM) and Normalized Difference Vegetation Index (NDVI). Details about these variables are given in S3 (SI). In addition, the 1 km grid population dataset of 2010 (Fu et al., 2010) was used to calculate the population exposure in Section 4.4.1. 2.4. Data integration As the aim of this study was to predict daily average PM2.5 con￾centrations at 3 × 3 km spatial resolution, a 3-km grid was created (1,073,351 grid cells in total) and exactly matched to the combined AOD pixels. The PM2.5 observations and other variables were re￾processed with a consistent spatial size and temporal interval for in￾tegration, because these data vary in spatial and temporal resolution. PM2.5 measurements from multiple stations located within a 3-km grid were averaged. The meteorological station data for each day were in￾terpolated to 3 km spatial resolution using an inverse distance weighted (IDW) technique, using the Spatial Analyst Tools in ArcGIS 10.2. Corresponding to the daily average PM2.5 measurements, the daily mean values obtained from averaging 3-hour PBL data were calculated. The two land-use datasets were separately aggregated to the 3 km grid cells. The final result was 60,810 matched samples incorporating ground PM2.5 and all independent variables, which were distributed over 325 days. 3. Development and validation Over the large area of mainland China, there is considerable spatial variability in the relationship between AOD and PM2.5. In addition, temporal variation significantly affects this relationship, and PM2.5 le￾vels on previous days may affect the concentrations on the day of measurement to a certain degree. Therefore, we developed a space-time regression model, the GTWR model, to simultaneously address spatial heterogeneity and temporal influence on the PM2.5-AOD relationship. To confirm the better performance of GTWR over other relevant models in previous studies, a daily GWR model using the same variables and matched samples as the GTWR model and a two-stage (LME + GWR) model (Hu et al., 2014) were also implemented for detailed compar￾isons. 3.1. GTWR model development 3.1.1. General structure of the GTWR model Unlike the widely used GWR model, which only takes spatial var￾iation into account when estimating an empirical relationship (Fotheringham et al., 2002), GTWR captures spatiotemporal hetero￾geneity based on a weighting matrix referencing both spatial and temporal dimensions. In this study, a GTWR model was fitted using the following structure: = + ×+ × + ×+ ×+ × + ×+ × + PM β μ υ t β μ υ t AOD β μ υ t NDVI β μ υ t RH β μ υ t WS β μ υ t T β μ υ t PBL β μ υ t DEM ε ( , ,) ( , ,) ( , ,) ( , ,) ( , ,) ( , ,) ( , ,) ( , ,) i i i i i ii i i ii i i ii i i ii i i ii i i ii i i ii i i 2.5 01 2 34 5 6 7 (1) where PM2.5i is the daily surface PM2.5 concentration of the sample i at location (μi,υi) on day ti; β0 denotes the intercepts at specific location (μi, υi) and day ti; and β1–β7 are the location-time-specific slopes for combined AOD, aggregated NDVI, interpolated RH, WS and T, averaged PBL, and aggregated DEM, respectively. The location (μi, υi) represents the central coordinates of a grid cell in which sample i is located. ti is the day of the year, and only the estimation for that day and previous days (ti ≤ testimation) are used for modeling. εi is the error term for sample i. To estimate the intercept of β0(μi,υi,ti) and the slopes of βi(μi,υi,ti) for each variable, the weight matrix W is introduced to account for the importance of sample j to the estimated parameters of sample i in terms of space and time (Huang et al., 2010). In this study, the widely used Gaussian distance decay-based functions and the combination method of spatial and temporal distance proposed in Huang et al. (2010) were applied to measure the spatiotemporal weight W. The details are pre￾sented in S4 (SI). 3.1.2. Exploratory analysis of independent variables Because of the relatively weak relationship between AOD and grid￾collocated PM2.5 across China (correlation value of 0.43), meteor￾ological and land-use data were incorporated into the GTWR model as covariates to improve the spatial consistency of the PM2.5-AOD re￾lationships (SI, S3). However, if all parameters mentioned in Section 2 were used as independent variables, the result would exist a collinearity problem. Variance inflation factors (VIF) were calculated as the colli￾nearity diagnostic for the independent variables (combined AOD, NDVI, RH, WS, T, P, PBL and DEM). The results demonstrated that the VIF values of pressure and DEM were very high (VIF > 10), which in￾dicates a strong collinearity problem if all eight independent variables are included in the PM2.5 prediction model (SI, S5). To address this issue, we excluded pressure and kept the seven other variables to fit the final GTWR model (the equation is given in Eq. (1)). 3.1.3. Parameter selection and two GTWR models The spatiotemporal bandwidth hST and scale factor ρ are the two key parameters in GTWR modeling (Huang et al., 2010). In the original GTWR (Huang et al., 2010), the popular method of plotting the sum of squared error against the parameters hST and ρ through cross-validation (Fotheringham et al., 2002), i.e., calculating the scores of CVRSS(hST; ρ) for hST and ρ (Eq. (S2) in SI), was used to select the two optimal values. However, this method is inappropriate for GTWR modeling over a large area, because it requires the use of a nested-loop, which can be easily implemented but consumes a large amount of computation time in search of the two optimal key parameter values (hST and ρ) for our large dataset. Selecting the best hST and ρ is to efficiently find the minimal values for a constrained nonlinear multivariable problem (Huang et al., 2010). Thus, an optimization approach based on the interior point al￾gorithm (IPA) (Byrd et al., 1999; Waltz et al., 2006) was devised in this study. Aided by a self-concordant barrier approach, IPA can expeditiously handle a nonlinear problem by approximately solving barrier sub-pro￾blems using a sequential quadratic programming iteration with trust region techniques (Byrd et al., 1999). Compared with other optimiza￾tion methods such as the active-set algorithm (Nocedal and Wright, 2006), IPA can efficiently handle large, sparse problems, which is es￾sential to a large dataset of PM2.5-AOD pairs. Using the IPA approach to determine the two optimal parameter values in this study, the compu￾tational time consumed for our sample dataset was reduced by ap￾proximately 60%, from about 2.5 days for the previous nested-loop method to about 1 day (SI, S4). In addition, to reduce the model overfitting problem, a fixed bandwidth regime rather than adaptive regime was applied to derive (SI, S4). With the GTWR model for the entire country using a general bandwidth of 120 km·day, we, however, could not predict robust PM2.5 concentrations extensively in northwestern China for some days due to the limited number of matched samples (SI, S6). To stabilize the model performance in areas where few matched samples were available, two Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 75

Q.He,B.Hang Remote Sensing of Ernir cnr206(2018)72-83 e gtwr models were fitted for bined AOD was 0.50 with a Std Dey of 0.35.The p D with e model structure as Eq(1).This ap reas,but the m acy (SI for the west were a licd in the final model.West r(/m),and the nimum was reached in summe with the highestin summer(and the ained the )A scale factor p of 0.2 with )The 10-fold l statistical indi eR.root mean al ground PM2s exhibits a de heast to and g dif timated pM aphy.and meteoro ditions.Ge t wa nd semi-arid areas 3.2.Daily GWR model and twe mode 4.2.Perfor ance evaluation of the GTWR mode thus be W m g al tive overall Rand RMSE values bet sand West)mode orm stably due to the .1and RM comparison d China into arts bec e the 4.Results 4.1.Descriptive statistics Table i shows the de stics of the that ode the entire s um with.5 form ed by Unit E ug/m Fa(N58,707 West (N =2103) 5535 76

separate GTWR models were fitted for eastern (E model) and western (W model) China with the same model structure as Eq. (1). This ap￾proach decreased the effect of the lower number of samples for western China on model performance without degrading overall accuracy (SI, S6). Therefore, fixed bandwidths of 100 km ∙ day for eastern China and 290 km ∙ day for the west were applied in the final model. Western China was defined as the provinces of Xinjiang, Qinghai, Tibet, and the western parts of Inner Mongolia, Gansu, and Sichuan; eastern China contained the remaining areas (SI, Fig. S4). A scale factor ρ of 0.2 with the best model performance was obtained. A 10-fold cross-validation (CV) technique was used to examine model overfitting (Rodriguez et al., 2010). To evaluate model perfor￾mance, several statistical indicators—R square (R2 ), root mean square error (RMSE), and mean prediction error (MPE)—were calculated by comparing differences between the estimated PM2.5 concentrations and surface PM2.5 measurements. The Pearson correlation coefficient was used in this study to quantify the correlation (r). 3.2. Daily GWR model and two-stage model To evaluate the performance of GTWR based on its integrated spa￾tiotemporal weighting scheme, we also fitted daily GWR (D-GWR) models and two-stage models with the structure as in S7 (SI). The D￾GWR models separately calibrate the PM2.5-AOD relationship for in￾dividual days, and could thus be regarded as a special case of GTWR if the number of estimation days equals one. The two-stage model (Hu et al., 2014) performs an LME analysis in the first stage and a GWR analysis in the second stage. We first attempted two separate D-GWR (East and West) models using the same partitions as in the Seg-GTWR model, and the results show that the West D-GWR model failed to perform stably due to the insufficiency of samples in western China. For the two-stage model, the results of the West model show a large overfitting problem, with R2 decreasing from 0.82 in model fitting to 0.66 in model validation. Therefore, we only developed East D-GWR and two-stage models for comparison. 4. Results 4.1. Descriptive statistics Table 1 shows the descriptive statistics of the grid-averaged PM2.5 concentrations in 60,810 grid cells and grid-collocated independent variables distributed over the 325 days used in fitting the models. The mean PM2.5 concentration averaged for the entire study region was 60 μg/m3 with a standard deviation (Std. Dev) of 41 μg/m3 , and the mean combined AOD was 0.50 with a Std. Dev of 0.35. The mean PM2.5 values did not differ greatly between the two areas, but the mean aerosol loading was much higher in eastern versus western China (0.50 vs. 0.34). The values of the dependent and independent variables pre￾sented different seasonal variability. The maximum mean PM2.5 oc￾curred in winter (86 μg/m3 ), and the minimum was reached in summer (44 μg/m3 ). The mean combined AOD values did not differ much among the four seasons, with the highest in summer (0.56) and the lowest in autumn (0.43). The detailed statistics for each season are listed in Table S7 (SI). The geographical distribution of the annual and seasonal mean PM2.5 measurements and combined AOD is shown in Fig. S5 (SI). Spatial ground PM2.5 exhibits a descending gradient from northeast to southwest due to differences in anthropogenic activities, land use, to￾pography, and meteorological conditions. Generally, high levels of PM2.5 concentrations and AOD values occur in densely populated areas with thriving industries and low elevations (e.g., the North China Plain, Sichuan Basin, and Central China) and in the arid and semi-arid areas of northwestern China (e.g., the Taklimakan Desert). 4.2. Performance evaluation of the GTWR model Fig. 2 shows the scatter plots of GTWR model fitting and CV results for eastern China, western China, and the entirety of mainland China (the Seg model represents the E + W model, and the non-Seg model is the GTWR model using all matched samples without segmentation). For the model fitting, the respective overall R2 and RMSE values between the estimated and observed PM2.5 are 0.85 and 15.78 μg/m3 based on the results of the E + W model. For model validation, CV R2 and RMSE are 0.80 and 18.58 μg/m3 , which is a slight improvement compared to the R2 and RMSE of 0.76 and 19.77 μg/m3 of the non-Seg model. The CV results suggest that the Seg model is slightly overfitted but the non￾Seg model is significantly overfitted, i.e., CV R2 decreases by 0.05 and 0.11 and RMSE increases by 2.80 and 5.21 μg/m3 for the two models, respectively. It is likely that the model itself could be improved by di￾viding mainland China into two parts, because the Seg model had comparable performance in model fitting, higher values in the CV processes, and much less overfitting than the non-Seg model. Therefore, we chose to fit two separate models to predict PM2.5 concentrations over China. The seasonal results of the Seg-GTWR model in fitting and CV processes based on the E and W models are shown in Fig. 3. It is clear that the model validation performance shows a similar seasonal pattern to that of model fitting. The values of R2 reach their maximum in au￾tumn with 0.85 for model fitting and 0.80 for CV, followed by winter with respective values of 0.84 and 0.77. Correspondingly, a lower R2 for Table 1 Descriptive statistics of the grid-collocated data for all variables used in the GTWR model. Region Unit PM2.5 AOD NDVI RH WS T PBL DEM μg/m3 – – % m/s °C m m Whole (N = 60,810) Mean 60 0.50 0.35 57.01 2.22 15.32 650.99 370.49 Std. Dev 41 0.35 0.15 15.06 0.96 8.82 302.16 516.87 Min 1 0.00 0.00 9.39 0.15 −13.90 20.39 −13.47 Max 673 3.89 0.90 97.00 13.64 34.12 2137.74 4553.23 Median 50 0.43 0.34 59.01 2.03 16.94 602.01 97.49 East (N = 58,707) Mean 60 0.50 0.35 57.58 2.23 15.30 646.02 346.66 Std. Dev 40 0.35 0.15 14.85 0.96 8.84 296.68 492.17 Min 1 0.00 0.00 9.45 0.15 −13.90 20.39 4.09 Max 595 3.89 0.90 97.00 13.64 33.69 2137.74 3323.93 Median 51 0.43 0.34 59.65 2.03 16.95 598.48 87.58 West (N = 2103) Mean 57 0.34 0.29 40.93 1.97 16.03 789.72 1035.73 Std. Dev 62 0.32 0.11 11.80 0.70 8.22 403.83 710.44 Min 2 0.00 0.08 9.39 0.54 −8.44 75.84 −13.47 Max 673 2.94 0.67 76.89 6.66 34.12 2095.89 4553.23 Median 37 0.25 0.29 39.60 1.87 16.69 761.45 934.51 Q. He, B. Huang Remote Sensing of Environment 206 (2018) 72–83 76