The Spatial Patterns of Land Surface Temperature and Its Impact Factors : Spatial Non-Stationarity and Scale Effects Based on a Geographically-Weighted Regression Model

Understanding the spatial distribution of land surface temperature (LST) and its impact factors is crucial for mitigating urban heat island effect. However, few studies have quantitatively investigated the spatial non-stationarity and spatial scale effects of the relationships between LST and its impact factors at multi-scales. The main purposes of this study are as follows: (1) to estimate the spatial distributions of urban heat island (UHI) intensity by using hot spots analysis and (2) to explore the spatial non-stationarity and scale effects of the relationships between LST and related impact factors at multiple resolutions (30–1200 m) and to find appropriate scales for illuminating the relationships in a plain city. Based on the LST retrieved from Landsat 8 OLI/TIRS images, the Geographically-Weighted Regression (GWR) model is used to explore the scale effects of the relationships in Zhengzhou City between LST and six driving indicators: The Fractional Vegetation Cover (FVC), the Impervious Surface (IS), the Population Density (PD), the Fossil-fuel CO2 Emission data (FFCOE), the Shannon Diversity Index (SHDI) and the Perimeter-area Fractal Dimension (PAFRAC),which indicate the vegetation abundance, built-up, social-ecological variables and the diversity and shape complexity of land cover types. Our findings showed that the spatial patterns of LST show statistically significant hot spot zones in the center of the study area, partly extending to the western and southern industrial areas, indicating that the intensity of the urban heat island is significantly spatial clustering in Zhengzhou City. In addition, compared with the Ordinary Least Squares (OLS) model, the GWR model has a better ability to characterize spatial non-stationarity and analyze the relationships between the LST and its impact factors by considering the space-varying relationships of different variables, especially at the fine spatial scales (30–480 m). However, the strength of GWR model has become relatively weak with the increase of spatial scales (720–1200 m). This reveals that the GWR model is recommended to be applied in the analysis of UHI problems and related impact factors at scales finer than 480 m in the plain city. If the spatial scale is coarser than 720 m, both OLS and GWR models are suitable for illustrating the correct relationships between UHI effect and its influence factors in the plain city due to their undifferentiated performance. These findings can provide valuable information for urban planners and researchers to select appropriate models and spatial scales seeking to mitigate urban thermal environment effect.


Introduction
Rapid and unprecedented urbanization has occurred across China over the past several decades-its urbanization rate, which can be defined as the proportion of urban population in the total population, has increased from 17.92% to 56.10% in 1978-2015 (National Bureau of Statistics of China, http://data.stats.gov.cn).With the development of urbanization, the urban areas have expanded and urban population has risen dramatically.The natural land cover is gradually being replaced by artificial materials (e.g., cement, asphalt, metals, etc.), which has changed the local urban landscape and climatic conditions, affected ecosystem functions and services and caused many ecological and environmental problems [1,2].One of the most prominent phenomena is the increasingly significant urban heat island (UHI) effect.The UHI effect refers to the phenomenon that the urban areas have higher temperatures than their surrounding suburban areas [3][4][5], which has attracted wide concern because of its adverse effect on human living environment, air quality and resident health [6][7][8][9].
In general, the UHI described in the term of air temperatures by the direct in-situ measurement and surface temperature from thermal remote sensing data [4,10].With the development of remote sensing technology, satellite-based estimation of Land Surface Temperature (LST) and other surface temperature (e.g., lake surface temperatures, etc.) derived from thermal infrared remote sensing imagery [11,12] becomes effective in interpreting UHI effect because it provides a relatively rapid and low-cost data information on landscape scale.Currently, numerous studies have been conducted to deal with the spatial characteristics, changing trends and impact factors of UHI effect based on remote sensing data and geographic information systems (GIS) technology [13][14][15][16][17].Among the above, the relationships between LST and related indicators become a hot topic that attracts wide-ranging research because these related findings have great significance to formulate urban planning policies to relieve UHI effect.Traditionally, some important indicators such as Normalized Difference Vegetation Index (NDVI), Fractional Vegetation Cover (FVC) and Normalized Difference Built-up Index (NDBI) have been commonly applied to examine the spatial variations of LST [18,19].Yuan and Bauer used the indicators including the NDVI and percent Impervious Surface Area (%ISA) based on Landsat imagery to quantitatively investigate the relationships between LST, percent impervious surface area and the NDVI from four different seasons [20].Several works have also used other biophysical parameters like, Normalized Difference Water Index (NDWI), the Normalized Difference Bareness Index (NDBaI) and other land-cover types as complimentary metrics to research UHI effect deeply by using remote sensing data [21][22][23].In addition, recent studies have developed various landscape pattern metrics such as landscape shape index, landscape configuration, patch density and edge density to measure the influence of these indicators on LST [24][25][26].However, there are still few UHI effect studies that have considered other indicators related to human activities into their respective researches, such as air pollution, energy consumption and demographic shifts [27,28].In this study, the Fossil-fuel CO 2 Emission and Population Density were incorporated into our impact factors analyses, in addition to other land use/land cover (LULC) variables and landscape pattern metrics.
Pervious literature on the relationships between LST and its influence factors commonly adopted Pearson correlations and the global regression models, like single or multiple linear regression and principal component regression [29][30][31].Although conventional correlation analysis or multivariate regression relationships are relatively well-established, these statistical analyses of existing studies have generally been aspatial and model the relationships without considering the spatially dependency of LST in the whole study area [25,32].In fact, the observed geographical and ecological patterns and processes tend to be spatial variable and the relationships between LST and its impact factors are often characterized by local changes [33,34].Commonly, this feature is referred to as spatial non-stationary [35].However, most of the existing researches are based on multiple linear regression analysis to model the relationship hiding the important details in the spatial variation and resulting in a failure to capture the spatial dependence of the data [36].This may generate misleading parameter estimates and uncertain significance test results [33].Hence, the impacts of related influence factors on LST considering the effects of spatial non-stationary need to be further investigated.In order to overcome the limitations of the aforementioned problem, this study examines the relationships between LST and its impact factors considering geographically-weighted regression (GWR) model, which developed a local regression technique to explore spatial non-stationarity in the relationship between variables at spatial scales [37][38][39].A limited number of studies have employed the GWR model when implementing a regression analysis of the relationships between LST and its impact factors [40,41].The studies by Buyantuyev [2] and Tian [42,43] showed that the GWR model can not only effectively solve the problem of the spatial non-stationarity of the relationships between LST and its impact factors but also provide better significance test results than conventional regression analysis models.
In addition, the existing studies have mainly investigated the relationships between LST and its influence factors at a single scale.However, LST and these impact variables as well as their spatial distribution and interactions exist scale dependence [33], namely, the spatial patterns of geographical and ecological processes may vary across scales, resulting in the need to explore the scaling-up effects on the relationships between LST and its various impact factors [43].Therefore, a multiscale analysis is necessary to uncover which the effects of related factors on LST vary across different scales.The main objectives of this study are as follow: (1) to analyze the spatial feature of LST and its impact factors, like landscape pattern metrics and social-ecological variables; (2) to explore and compare GWR model as well as traditional regression analysis such as the Ordinary Least Squares (OLS) model used to examine the relationships between LST and its impact factors; (3) to investigate scale effects on the association between LST and related explanatory variables, namely the change characteristics across different spatial resolutions (from 30-120, 240, 480, 720, 960, 1080 and 1200 m).Investigating and understanding the LST and related influence factors are important steps toward mitigating UHI effect and promoting urban sustainability [44], the findings of this study will contribute to our understanding of how the related impact factors may affect local LST at different spatial scales.This research will enrich the previous literature on the utility of GWR model at different scales and help develop appropriate urban planning policies to alleviate UHI effect.

Study Area
Zhengzhou City, situated in the north-central of China (112 • 42 E-114 • 14 E, 34 • 16 N-34 • 58 N), is an important central city and transportation hub in China (Figure 1).It covers an area of about 7446.2 km 2 and belongs to a typical inland plain city.This city is part of the North Temperate Zone continental monsoon climate and the annual average temperature is approximately 14.4 • C. The area has an annual precipitation of 640.9 mm and a terrain trend of higher terrain in the southwest and lower terrain in the northeast.The Yellow River goes through the area from west to east and delineates the northern boundary of Zhengzhou City.As the capital of Henan Province, Zhengzhou has been experienced rapid population growth and urbanization in the past several decades.According to the population census (http://www.zzstjj.gov.cn/), the urban population of Zhengzhou City increased from 1.25 to 4.78 million and the urbanization rate increased from 32.4% in 1978 to 68.3% in 2014 (Figure 2).Such a dramatic urbanization process inevitably brings urban heat island effect of the whole city because of the increasing built-up density and loss urban green spaces.We chose the central urban of Zhengzhou City and its surrounding expansion region as study area and the location map is shown as Figure 1.

Image Pre-Processing and Land-Cover Classification
In this study, the Landsat 8 OLI/TIRS (path124/row36) images acquired on 27 May 2014 were the primary data for deriving LST and Land use and Land cover (LULC) types.The images downloaded from the USGS data web (http://www.earthexplorer.usgs.gov)and cloud-free images or images with minimal cloud cover (less than 5%) were considered.Prior to LULC classification mapping and LST retrieval, these satellite images were subjected to a set of pre-processing procedures.The pre-processing included radiometric calibration, atmospheric correction (darkobject subtraction) using the TerrSet software (https://clarklabs.org/terrset/)and geometrical distortions correction.Then, the images were further re-sampled with pixel sizes of 30 m by 30 m for all bands, including the thermal band.

Image Pre-Processing and Land-Cover Classification
In this study, the Landsat 8 OLI/TIRS (path124/row36) images acquired on 27 May 2014 were the primary data for deriving LST and Land use and Land cover (LULC) types.The images downloaded from the USGS data web (http://www.earthexplorer.usgs.gov)and cloud-free images or images with minimal cloud cover (less than 5%) were considered.Prior to LULC classification mapping and LST retrieval, these satellite images were subjected to a set of pre-processing procedures.The pre-processing included radiometric calibration, atmospheric correction (darkobject subtraction) using the TerrSet software (https://clarklabs.org/terrset/)and geometrical distortions correction.Then, the images were further re-sampled with pixel sizes of 30 m by 30 m for all bands, including the thermal band.

Image Pre-Processing and Land-Cover Classification
In this study, the Landsat 8 OLI/TIRS (path124/row36) images acquired on 27 May 2014 were the primary data for deriving LST and Land use and Land cover (LULC) types.The images downloaded from the USGS data web (http://www.earthexplorer.usgs.gov)and cloud-free images or images with minimal cloud cover (less than 5%) were considered.Prior to LULC classification mapping and LST retrieval, these satellite images were subjected to a set of pre-processing procedures.The pre-processing included radiometric calibration, atmospheric correction (dark-object subtraction) using the TerrSet software (https://clarklabs.org/terrset/)and geometrical distortions correction.Then, the images were further re-sampled with pixel sizes of 30 m by 30 m for all bands, including the thermal band.
To produce the land-cover classification maps of the study area, a supervised classification and man-computer interactive interpretation were utilized in ENVI 5.1.The final urban landscape was divided into six classifications: water, impervious surfaces (IS), farmland, bareland, grassland and forest.During this process, Google Earth images with high spatial resolution images and other auxiliary data were used as source of reference information to help identify the land use types.Overall, the total accuracy of the land use classification achieves over 86.6%, meeting the recommended value by Janssen et al. [45].Hence, these data were available for further research.

Dependent Variable: LST
The land surface temperature (LST) was retrieved from the thermal infrared bands of Landsat satellite data.After the digital number (DN) values were converted to the thermal infrared spectral radiance, which were applied to obtain top of atmosphere (TOA) brightness temperature values expressed in Kelvin [15,46].We estimated the emissivity of surface materials (ε) based on the normalized difference vegetation index value (NDVI), expressed in Equation ( 1) [47]: where denotes the soil emissivity and ε s denotes vegetation emissivity, respectively.P v is the vegetation proportion [48].The NDVI is derived using the surface reflectance of Landsat images by Equation ( 2).

NDVI = (
where NIR = Band 5 (for Landsat 8) and Red = Band 4 (for Landsat 8) [15,19,49].Lastly, the emissivity-corrected LST values were measured as follows (Equation ( 3)) [50,51]: where T B = satellite brightness temperature in Kelvin; λ = wavelength of emitted radiance (λ = 10.8 µm for Landsat OLI/TIRS Band 10 [33]) The fractional vegetation cover (FVC) is considered as an important biophysical parameter to determine the land surface properties and widely used for the land mapping in previous researches [52,53].FVC mainly depicts the vegetation abundance of ground surface.Compared with the NDVI, FVC can better manifest the correlation with LST values that are retrieved from remote sensing images.Therefore, the FVC was utilized as an impact factor of vegetation coverage to analyze the effect on the land surface temperature.It can be obtained according to the NDVI and expressed as [19]: where NDVI i is NDVI value of a pixel i; NDVI S can be approximately equal to the minimum value of NDVI; and NDVI V can be about equal to the maximum value of NDVI.So, Equation ( 5) can be expressed as: Extraction of Impervious Surface (IS) According to the study by Wang et al. [54], the Impervious Surface (IS) indicator can be obtained from FVC and expressed in Equation ( 6): Population Density (PD) Data The high population density indicates the rapid expansion of settlements and dense impervious surface areas in a region, which can generate a mass of anthropogenic heat.Therefore, population density is considered to be one of the main contributors to the formation of urban heat island [55].In this study, the finest resolution (100 m) patterns of population distribution for mainland China are used to reveal the relationship between LST and population density in Zhengzhou City, this population data has been derived from the World Pop Project website and widely used by researchers [56][57][58].

Fossil-Fuel CO 2 Emission Data (FFCOE)
The fossil-fuel CO 2 emission has increased rapidly in recent years due to urbanization and industrialization [59], which leads directly to air pollution, global warming and other environmental problems [60,61].Hence, fossil-fuel CO 2 emission has an important effect on UHI.In this study, the fine-scale (1 × 1 km) fossil-fuel CO 2 emission data derived from the ODIAC Fossil Fuel Emission Dataset were employed to examine the correlation with LST [62].

The Shannon Diversity Index (SHDI)
In order to consider the effect of ecological processes and landscape patterns on LST, landscape metrics are applied to quantify landscape shape and configuration on remote sensing images.In this study, two selected landscape metrics were used to analyze the relationships between landscape patterns derived from LULC and LST.The Shannon diversity index (SHDI) based on various land-cover types was calculated, which is a useful metric of landscape configuration to measure the diversity in the whole community ecology [63].The Equation ( 7) according to Shannon and Weaver [64] as follows: where p i denotes the proportion of land-cover types in the study area and then multiplied by the natural logarithm (ln p i ).

Perimeter-Area Fractal Dimension (PAFRAC)
In landscape ecology, the Perimeter-area fractal dimension (PAFRAC) provides an index of the overall shape complexity of the patches across a wide range of spatial scales [63].PAFRAC can be expressed as follows [54]: where a ij denotes the area of patch ij, p ij denotes the perimeter of patch ij, N denotes the number of patches in the landscape of patch type.

Geographically-Weighted Regression (GWR) Model
In order to highlight the good performance of the GWR model, the OLS model, that is one of the classic global regression method, is also performed in the study area.The OLS model can be expressed as follows [24]: where y i denotes the dependent variable, x k denotes the explanatory variable, β 0 denotes the intercept, β k denotes the estimated coefficient and ε represents the random error.
The GWR model can be considered as a type of regression analysis method with geographically varying parameters, which was given a fully description of its algorithm and the principle by Fotheringham et al. [65].A conventional GWR model is described by the Equation ( 10): where y i , x k,i and ε i represent dependent variable, independent variable and the Gaussian error at location i, respectively; β 0 (u i , v i ) represents the intercept at the point i, (u i ,v i ) is the x-y coordinate of the i th location; and β k (u i , v i ) represents the local coefficients which are varying conditionals on the location and can be estimated by the Gaussian model.Such modelling is likely to obtain higher performance than conventional regression models.Additionally, in order to properly compare the performance of the OLS and GWR models, the following tests were utilized: Adjusted R-Squared (R 2 ), which is a measure of goodness-of-fit and a higher Adjusted R 2 means the method has better performance.The corrected Akaike Information Criteria (AICc), Global Moran's I index (Moran's I) and the approximate likelihood ratio test based on the F-value, which is presented for the null hypothesis that the GWR model has no correction over the OLS model [23].

Hot Spot Analysis
In order to investigate the spatial distribution of urban heat island in the study area, the Getis-Ord Gi* method was used to measure the cold-hot spots of LST using ArcGIS 10.2 software.The Getis-Ord G i * method was expressed by the following Equation ( 11) [66]: where X i and X j denote the variate value of pixels at location i, j, W ij (d) denotes a binary spatial weight matrix.The standardized Getis-ord statistic was applied to assess the significance of the Getis-ord statistic [67].The standardized version of the Getis statistic is Equation ( 12): where E(G i

Hot and Cold Spots in Zhengzhou City
The hot-cold spots analysis was based on the mean land surface temperature (LST) derived from remote sensing data.The threshold distance value (90 m) for this test was statistically determined by applying the Incremental Spatial Autocorrelation tool in ArcGIS 10.2 software.According to the Getis-Ord Gi* model, the Z-Score values in the range of 2.580-3.289(the 99% confidence level) indicated a statistically significant result and this threshold of Z(G i * ) value can be used to determine the hot spot zones, while the Z-Score Z(G i * ) values in the range of −3.289-−2.580(the 99% confidence level) showed a statistically significant result and this threshold of Z(G i * ) value can be used to determine the cold spot zones (Figure 3).The results show statistically significant hot spot zones in the center of Zhengzhou City, which stretched partway into the western and southern industrial areas, such as High-tech industrial development zone and Zhengzhou Economic and Technical Development Zone.
The city periphery areas (in the North and Southwest directions) are defined as statistically significant cold spots.As expected, the areas located in near the Yellow River green space and vegetation coverage show up as cold spots, suggesting that green space has a positive effect on reducing urban thermal environment.Overall, it can be indicated that the intensity of urban heat island is significantly spatial clustering in Zhengzhou City.
where E(Gi * ) is expected value of Gi * , Var(Gi * ) is variance value of Gi * .If Z(Gi * ) > 0 and statistically significant, this indicates a spatial clustering of values and belongs to hot spot zones.If Z(Gi * ) < 0 and statistically significant, this indicates a spatial clustering of values and belongs to cold spot zones.

Hot and Cold Spots in Zhengzhou City
The hot-cold spots analysis was based on the mean land surface temperature (LST) derived from remote sensing data.The threshold distance value (90 m) for this test was statistically determined by applying the Incremental Spatial Autocorrelation tool in ArcGIS 10.2 software.According to the Getis-Ord Gi* model, the Z-Score values in the range of 2.580-3.289(the 99% confidence level) indicated a statistically significant result and this threshold of Z(Gi * ) value can be used to determine the hot spot zones, while the Z-Score Z(Gi * ) values in the range of −3.289-−2.580(the 99% confidence level) showed a statistically significant result and this threshold of Z(Gi * ) value can be used to determine the cold spot zones (Figure 3).The results show statistically significant hot spot zones in the center of Zhengzhou City, which stretched partway into the western and southern industrial areas, such as High-tech industrial development zone and Zhengzhou Economic and Technical Development Zone.The city periphery areas (in the North and Southwest directions) are defined as statistically significant cold spots.As expected, the areas located in near the Yellow River green space and vegetation coverage show up as cold spots, suggesting that green space has a positive effect on reducing urban thermal environment.Overall, it can be indicated that the intensity of urban heat island is significantly spatial clustering in Zhengzhou City.

The Relationships between the LST and Explanatory Variables at a Single Scale
In this section, the bivariate relationships between LST as the dependent variable and the explanatory variables FVC, IS, PD, Fossil-fuel CO2 Emissions (FFCOE), SHDI and PAFRAC are presented in Table 1.The GWR model was employed to explore the spatial non-stationarity relationships between LST and related impact factors at a 30 m spatial scale (mesh size of each pixel).The corresponding OLS models with the same variables were also used as a reference.The derivation

The Relationships between the LST and Explanatory Variables at a Single Scale
In this section, the bivariate relationships between LST as the dependent variable and the explanatory variables FVC, IS, PD, Fossil-fuel CO 2 Emissions (FFCOE), SHDI and PAFRAC are presented in Table 1.The GWR model was employed to explore the spatial non-stationarity relationships between LST and related impact factors at a 30 m spatial scale (mesh size of each pixel).The corresponding OLS models with the same variables were also used as a reference.The derivation results of FVC and land cover classification of the study area derived from the Landsat 8 OLI/TIRS were shown in Figure 4.As shown in Figure 4a, the spatial distribution of FVC in study area exists significant differences, that is to say, the lower values of FVC are mainly distributed in the city center, while the surrounding areas exhibit high values of FVC.For the LULC map (Figure 4b), the spatial pattern of impervious surface (IS) areas cover the highest proportion in the study area, corresponding to low values of FVC and higher LST.The forest and grassland have the smaller proportion of area, corresponding to high values of FVC and lower LST.As shown in Figure 4a, the spatial distribution of FVC in study area exists significant differences, that is to say, the lower values of FVC are mainly distributed in the city center, while the surrounding areas exhibit high values of FVC.For the LULC map (Figure 4b), the spatial pattern of impervious surface (IS) areas cover the highest proportion in the study area, corresponding to low values of FVC and higher LST.The forest and grassland have the smaller proportion of area, corresponding to high values of FVC and lower LST.

Contrastive Analysis between the GWR and OLS Models
According to the evaluation results shown in Table 1, all the GWR models have much better performance than OLS models in the six explanatory variables.When compared with the evaluation indices (AICc and Adjusted R 2 ), it is obvious that all GWR models had higher Adjusted R 2 values, smaller AICc values than the corresponding OLS models.For the single impact factor analysis, the IS variable had the most important influence on the LST with an AICc value of 6793.64 and Adjusted R 2 value of 0.92 for the GWR model and with an AICc value of 15,626.42 and Adjusted R 2 value of 0.52 for the OLS model, followed by FVC, PD, FFCOE, PAFRAC and SHDI.For the multiple impacting factors results, evaluated with the IS, FVC, PD, FFCOE, SHDI, PAFRAC simultaneously, manifested

Contrastive Analysis between the GWR and OLS Models
According to the evaluation results shown in Table 1, all the GWR models have much better performance than OLS models in the six explanatory variables.When compared with the evaluation indices (AICc and Adjusted R 2 ), it is obvious that all GWR models had higher Adjusted R 2 values, smaller AICc values than the corresponding OLS models.For the single impact factor analysis, the IS variable had the most important influence on the LST with an AICc value of 6793.64 and Adjusted R 2 value of 0.92 for the GWR model and with an AICc value of 15,626.42 and Adjusted R 2 value of 0.52 for the OLS model, followed by FVC, PD, FFCOE, PAFRAC and SHDI.For the multiple impacting factors results, evaluated with the IS, FVC, PD, FFCOE, SHDI, PAFRAC simultaneously, manifested the best performance for the GWR model, with an AICc value of 11,818.52 and Adjusted R 2 value of 0.82, while the OLS model provided results with an AICc value of 15,460.78 and Adjusted R 2 value of 0.59.For the statistical comparison between two models in Table 2, there is a significant improvement in LST and its impact factors investigation by employing the GWR method.The simulated residuals of the GWR model (1567.66) is less than that of the OLS model (8744.31).Therefore, it is believed that the computational accuracy of GWR model is better than OLS model.As two kinds of statistical models, the GWR model provides the best performance for analyzing the relationships between LST and its impact factors.This result is consistent with previous research conclusions [36,43].Figure 5 shows the results of local coefficients estimated by GWR model to explore the spatial varying relationship between the LST and the six impact factors.It can be found that all explanatory variables exhibited a statistically significant contribution on impacting LST in Zhengzhou City and the local coefficients for all the indices change with the different space locations.As expected, the FVC has a negative relation with LST with coefficient values from −6.4 to −0.63, while the IS, PD and FFCOE have positive influence on the LST with coefficient values from 0.53 to 7.54, from 0.04 to 2.56 and from 0.01 to 0.94, respectively.The local coefficients values of the PAFRAC and SHDI are from −23.57 to 42.11 and from −35.08 to 17.92, respectively and have both negative and positive correlations to the LST.This indicates the unstable relationships among LST, PAFRAC and SHDI variables.Unlike the OLS model, the spatial distribution of the six explanatory variables effect on LST shows spatial differentiation in GWR model.The FVC estimated coefficient has higher negative influence on LST in the northern part of the city and a lower but still negative, influence in the city center and southeast part, mainly built-up and industrial part of Zhengzhou City.The IS has a higher positive influence on the LST mainly in the city center or in the western and southern part of the city and a slightly lower but still positive, influence in the northern part of the city.The PD and FFCOE had higher values in the city center and lower values in city surroundings.The PAFRAC estimated coefficient is negative relation in most of the city area and contrastingly, in the southern surroundings of the city, in a slightly positive relation.The SHDI local estimate has a negative influence on the LST in most of the city area, whereas it is positive relation in the other regions.In general, it can be reflected that the GWR model has good ability to characterize spatial non-stationarity among relationships in study area.

Effects of the Relationships between the LST and Explanatory Variables at Different Resolutions
In this section, we mainly investigated the spatial scale effects of the relationships between the LST and the six impact factors at multi-scales.The LST, FVC, IS, PAFRAC, PD, SHDI, FFCOE images were aggregated from 30-120, 240, 480, 720, 960, 1080 and 1200 m spatial resolutions using ArcGIS 10.2 software (ESRI, California, USA).Additionally, the relationships between the six impact factors and LST were estimated based on GWR model at multi-scales.We eventually obtained the maps of the local R 2 and the residuals of all impact factors at different spatial scales (mesh size) shown in Figures 6 and 7, respectively.Local R 2 is applied to explore the ability of the GWR model to fit the observations, which have high local R 2 values will have good performances.As shown in Figure 6, the local R 2 values of the whole Zhengzhou City become higher with the increase of the spatial resolution (Figure 6a-h).The distribution range of local R 2 changes from 0.20-0.82(30 m) to 0.4-0.95(1200 m) and the mean local R 2 value increases from 0.65 (30 m) to 0.82 (1200 m).However, the standard deviation (SD) values of local R 2 decreases from 0.44 (30 m) to 0.11 (1200 m).In addition, we also found that the spatial differentiation of the spatial pattern of the local R 2 was gradually reduced from the 30 m to 1200 m spatial scale, especially after the spatial resolution was increased to 720 m.It is suggested that the relationships are becoming more generalized spatial distributions and spatial non-stationarity tends to be neglected with the increase of the spatial scales.This finding is also verified by the research of Luo and Peng [23].For the spatial pattern of local R 2 , the values of local R 2 are changing with the spatial location, characterized by relatively higher values in the city surroundings and obviously lower values in the city center.This is because relatively frequent human activities in the city center seriously influence on urban thermal environment.
Furthermore, as Figure 7a-h shows, we can find that the spatial patterns of local residuals change with the spatial location and this indicates that GWR model has a good ability to show the spatial difference of multi-factors, with values from −2.45 to 2.98.Among them, the spatial distribution of the residuals at 30 m appears to be about a random distribution and with the increase of spatial scales, the residual distribution is becoming a certain clustering.From the statistical results of residuals in Table 3, it can be found that the variation range of residuals values is reduced from −2.45-2.98 (30 m) to −1.82-1.98 (1200 m) and the standard deviation (SD) values of residuals decrease from 0.79 (30 m) to 0.51 (1200 m).This indicates that with the increase of the spatial resolution, the GWR model tends to neglect the spatial non-stationarity and become more generalized geographical patterns.This result is consistent with above-mentioned study in Figure 5 and the previous researches [23].According to the Table 3, we also found that the Moran's I of residuals gradually increased from 0.12 (30 m) to 0.31 (1200 m), it can be revealed that the spatial autocorrelation of residuals is gradually enhanced and the trends of spatial clustering are becoming more and more stronger.

Variations of the GWR and OLS Models at Different Scales
In this section, we mainly analyzed performance parameters, including AICc and adjusted R 2 , of GWR and OLS models at different spatial resolutions.As shown in Table 4, it can be found that the AICc values of both GWR and OLS models are gradually reducing with the increase of the spatial scale, the values are in ranges from 11,818.52 (30 m) to 1951.86 (1200 m) and from 15,460.78 (30 m) to 2432.02 (1200 m), respectively.On the contrary, the adjusted R 2 values of GWR and OLS models are continuously increasing along with the increase of the spatial resolutions, these values are in ranges from 0.82 (30 m) to 0.89 (1200 m) and from 0.59 (30 m) to 0.75 (1200 m).This finding reveals that with the variation of spatial scales, the accurate degree of the GWR and OLS models' evaluation would be changed and the analysis results of both models have good consistency, that is, the fitting degree of both models are gradually improved with the spatial scales increase from 30 m to 1200 m.
In addition, for the parameters comparison between the two models, the GWR model gives higher adjusted R 2 and lower AICc values than those from OLS model, especially at the small spatial scales (30-480 m).It can be revealed that the GWR model have better performance than the OLS model to explore the relationships between LST and its impact factors at fine spatial resolutions (30-480 m).However, as the spatial scales increased from 720 m to 1200 m, the advantage of the GWR model becomes relatively weak and the AICc values of both models are nearly equal to each other, implying that the fitting degrees of both models are relatively consistent at resolutions coarser than 720 m.This may be explained by the reason that the distance of the adjacent pixels has a positive relation to the bandwidth of the GWR model, when the spatial resolution is increased, the bandwidth is also accordingly increased and the parameter values of the GWR model become gradually global, implying that the spatial pattern is more and more generalized and the spatial non-stationarity tend to be neglected [22].Therefore, when we use the GWR model to study the relationships between LST and its impact factors, choosing the suitable spatial scale becomes very important in future researches.When the spatial scale is finer, the GWR model can better uncover the spatial non-stationarity and local variation of LST and its impact factors in detail.However, when the spatial scale is coarser, compared with the OLS model, the analytical advantage of the GWR model is gradually weak and the performance difference of both models is relatively decreasing due to more generalized geographical patterns.For Zhengzhou City or other cities on the plains, if the spatial resolution of related remote sensing data is finer than 480 m, the GWR model can have a better ability than OLS models to characterize spatial non-stationarity of LST and explain the relationships between LST and its impact factors.At any spatial scale coarser than 720 m, both the GWR and OLS models are suitable to study LST and explore its relationships with impact factors, while the GWR model has a better fitting degree than OLS model at the large spatial scale.

Implications of the Relationships between LST and Its Impact Factors
In this study, we explored the relationship between LST and FVC, IS, PD, FFCOE, SHDI and PAFRAC by using the GWR and OLS models in Zhengzhou City.Study results from the single-factor models indicated that these influence factors can affect LST significantly.The most intensive warming in Zhengzhou City occurs mainly in the central city stretching to the western and southern industrial areas and is strongly positively related to the IS variable.However, FVC is negatively correlated to the LST and the associations are stronger in the northern part of the city, approximately along the distribution of the Yellow River green space and other vegetation coverage areas, where the spatial distribution of LST was lower with high vegetation coverage.So, the LST can decrease with the increment of vegetation coverage rate.This result corresponded to the existing studies that vegetation coverage (such as forests, parks and grasslands) can well reduce urban thermal environment and increase urban cooling areas through evapotranspiration and shading effect [68][69][70].
It should be pointed out that population density (PD) and Fossil-fuel CO 2 emission (FFCOE) variables are positively correlated to the LST and the associations are stronger in the central center of the city, where rapid urbanization and economic development occurred.This implies that these social-ecological variables can affect the LST distribution to some extent.Study findings are consistent with the previous studies, which have proved that population growth, urban expansion and CO 2 emission enhance UHI effect [7].Higher the Shannon Diversity Index (SHDI) values, indicating the diversity of land cover types, are correlated with lower LST in the study area, this is because the different physical attributes of land cover types can have various effects on the urban thermal environment.Hence, increasing the landscape diversity of the study area is one of the effective ways to reduce the UHI effect.In addition, the Perimeter-area Fractal Dimension (PAFRAC) variable, showing the shape complexity of the landscape patches, have both negative and positive correlations to the LST.This indicates that higher PAFRAC values of vegetation coverage types, including the forest and grassland landscape patches, are correlated with lower urban thermal environment, while higher PAFRAC values of impervious surface areas can have a positive effect on higher LST.Therefore, according to the urban vegetation landscape playing an important role in reducing UHI effect, especially in central city areas, it is more important to suggest urban planners to increase the shape complexity of urban green space when they consider mitigating urban heating in Zhengzhou City.

Implications for Choosing GWR and OLS Models
Selecting a suitable method is important for investigating the relationships between LST and its influence factors.Due to the significant spatial heterogeneity of LST, the hypothesis that a series of evaluation coefficients cannot consider the space-varying relationships between variables in global regression models is unreasonable [70].Hence, using the OLS model without considering the spatial location of variables maybe result in unreliable parameter estimates.The Strengths of the GWR model to explore the relationships between LST and its impact factors are obvious.It has become a helpful method for exploring the spatial non-stationarity phenomena that occurs within the spatial patterns of LST and its impact factors.Some scholars have also the same findings that the significant spatial non-stationarity widely existed among the LST and its impact factors and explained their relationships by using the GWR model [35,71].Furthermore, compared with the OLS model, the GWR model has a better performance to characterize spatial non-stationarity and analyze the relationships by increasing the fitting degree of coefficient estimates.However, some researches have pointed out that a relatively high overall fitness is not conducive to better understanding the relationships [32] and neglects the detail information of variables at large spatial resolution with the cost of the smoothing of the existing geographic differences in pixels [23].This study also verifies the same results and considers selecting appropriate spatial scales to discuss this issue.

Limitations
In this study, we have provided a comprehensive analysis framework to explore the relationships between LST and its impact factors at different spatial scales in Zhengzhou City.However, this research has some limitations as follows: Firstly, the study analyzed the spatial patterns of LST and related impact factors by only using one daytime remote sensing image due to limited data accessibility but seasonal variations and climatic conditions which may to some extent affect the urban thermal environment were not considered.Therefore, multiperiod remote sensing data should be used in future study.In addition, the main attention of this research is to examine the relationships between LST and its impact factors, the comprehensive effects of human activities, urban size, traffic flow and other impact factors on the LST, are also worth analyzing.Besides, limited by accuracy of remote sensing imagery interpretation, we have not further differentiated the types of urban land use and land cover in analysis of influence factors.Hence, further studies need to be conducted considering extraction urban land-cover information from high resolution remote sensing images.

Implications of Suitable Spatial Scales in GWR and OLS Models
With the development of remote sensing methods and technologies, there has been increasing research attention on the effects of spatial scales or spatial resolutions in modelling the landscape ecology processes by using remote sensing data [21,24,28].Multiscale analysis in UHI effect studies become gradually important.Although some studies have been conducted to explore the relationships between LST and its impact factors at different scales in GWR model, they failed to find the appropriate spatial scales to illustrate the correct relationships though comparison between GWR and OLS models.In this study, it is worth noting that the GWR and OLS models have been employed simultaneously to obtain related parameter estimates at different scales, while the performance of the GWR model is better than the OLS model with the increase of spatial scales.The results indicated that the GWR model obtains higher adjusted R 2 and lower AICc values than those from OLS model, especially at the fine spatial scales (30-480 m).However, as the spatial scales increased from 720 m to 1200 m, the strength of the GWR model becomes relatively weak and the AICc values of both models are nearly equal to each other, implying that the fitting degrees of both models are relatively consistent at resolutions coarser than 720 m.Therefore, considering the two important factors for finding the appropriate scales: a high R 2 and a low AICc value [32], we can think that the GWR model should be advised to be firstly applied in the analysis of UHI problems and related impact factors at scales finer than 480 m in the plain city.If the spatial scale of remote sensing data is coarser than 720 m, both OLS and GWR models are suitable for illustrating the correct relationships between UHI effect and its influence factors in the plain city due to their undifferentiated performance.Some researchers have also found that the GWR model had better simulation with higher R 2 and lower AICc values than the OLS model at small spatial scales when they used both models to explore the relationships between LST and related impact factors [40,42].Luo et al. [23] used both models to test the relationships between LST and three factors (SAVI, IBI and NDSI) at different spatial scales and found that GWR model has better performance than the OLS model to explain the relationships in mountain city when the spatial scales are less than 240 m.However, when the spatial scales are greater than 480 m, the GWR and OLS models have few performance differences in illustrating the relationships.Therefore, based on the above analysis, we can reveal that selecting the optimal scales to examine the relationships between LST and its impact factors by using the GWR model is variant in different types of cities.This finding can provide a valuable reference for UHI effect studies to choose appropriate models at different spatial scales in plain city.

Conclusions
Based on remote sensing data, our study analyzed the spatial distribution of LST and examined the quantitative relationships between LST and its impact factors, including FVC, IS, PD, Population Density (PD), Fossil-fuel CO 2 Emission, the Shannon Diversity Index (SHDI) and the Perimeter-area Fractal Dimension (PAFRAC), at multiple scales by using GWR model.The spatial non-stationary and scale effects of LST-influence factors relationships were also discussed in Zhengzhou City and the conclusions can be drawn as follows: (1) This study indicates that the intensity of UHI is significantly spatial clustering in Zhengzhou City.Hot spot zones were clustered in urban center and the western and southern industrial areas.The Cold spots zones were located in the city periphery areas (in the North and Southwest directions).( 2) Study results from the single-factor models indicates that these influence factors can affect LST significantly.LST is strongly positively related to the IS variable.However, a negative relationship exists between LST and FVC.It should be pointed out that population density (PD) and Fossil-fuel CO 2 emission (FFCOE) variables are positively correlated to the LST, implying that these social-ecological variables can affect the magnitude of LST to some extent.In addition, the SHDI and PAFRAC variables, indicating the diversity and shape complexity of land cover types, respectively, have both negative and positive correlations to the LST in different areas.This reveals the unstable relationships between LST and PAFRAC and SHDI variables.(3) Overall, compared with the OLS model, the GWR model has a better ability to characterize spatial non-stationarity and analyze the relationships between the LST and its impact factors by considering the space-varying relationships of different variables, especially at the fine spatial scales (30-480 m).However, the strength of GWR model has become relatively weak with the increase of spatial scales (720-1200 m).Given the important principles of the higher R 2 and lower AICc values for finding the optimal scales, we can infer that the GWR model is recommended to be applied in the analysis of UHI problems and related impact factors at scales

Sustainability 2018 , 21 Figure 1 .
Figure 1.The geographical location of the Zhengzhou City in China.(The Landsat image applied in study area shown in a false color composite).

Figure 2 .
Figure 2. The urban population and urbanization rate of Zhengzhou City from 1978-2014.

Figure 1 .
Figure 1.The geographical location of the Zhengzhou City in China.(The Landsat image applied in study area shown in a false color composite).

Sustainability 2018 , 21 Figure 1 .
Figure 1.The geographical location of the Zhengzhou City in China.(The Landsat image applied in study area shown in a false color composite).

Figure 2 .
Figure 2. The urban population and urbanization rate of Zhengzhou City from 1978-2014.

Figure 2 .
Figure 2. The urban population and urbanization rate of Zhengzhou City from 1978-2014.

Figure 3 .
Figure 3.The spatial distribution of hot-cold spot analysis in Zhengzhou City.

Figure 3 .
Figure 3.The spatial distribution of hot-cold spot analysis in Zhengzhou City.

Sustainability 2018 ,
10, x FOR PEER REVIEW 9 of 21 results of FVC and land cover classification of the study area derived from the Landsat 8 OLI/TIRS were shown in Figure 4.

Sustainability 2018 , 21 Figure 5 .
Figure 5.The spatial distribution of the local coefficients for FVC, IS, PAFRAC, PD, SHDI, FFCOE (af) at a spatial resolution of 30 m.Figure 5.The spatial distribution of the local coefficients for FVC, IS, PAFRAC, PD, SHDI, FFCOE (a-f) at a spatial resolution of 30 m.

Figure 5 . 21 Figure 6 .
Figure 5.The spatial distribution of the local coefficients for FVC, IS, PAFRAC, PD, SHDI, FFCOE (af) at a spatial resolution of 30 m.Figure 5.The spatial distribution of the local coefficients for FVC, IS, PAFRAC, PD, SHDI, FFCOE (a-f) at a spatial resolution of 30 m.

Figure 6 . 21 Figure 7 .
Figure 6.The spatial distribution of the local R 2 values for the six impact factors at different spatial scales, (a-h) denote the spatial scales of 30 m, 120 m, 300 m, 480 m, 720 m, 960 m, 1080 m and 1200 m, respectively.

Figure 7 .
Figure 7.The spatial distribution of the residuals values for the six impact factors at different spatial scales.(a-h) denote the spatial scales of 30 m, 120 m, 300 m, 480 m, 720 m, 960 m, 1080 m and 1200 m, respectively.

Table 1 .
The results comparisons of the GWR and OLS models.

Table 1 .
The results comparisons of the GWR and OLS models.

Table 2 .
The Residuals comparison between the GWR and OLS models for the multi-factors.

Table 3 .
The descriptive statistics and spatial autocorrelation of residuals at different spatial scales (mesh size).

Table 4 .
Parameters comparison between GWR and OLS model at different scales.