Scale Effects of the Relationships between Urban Heat Islands and Impact Factors Based on a Geographically-Weighted Regression Model

Urban heat island (UHI) effect, the side effect of rapid urbanization, has become an obstacle to the further healthy development of the city. Understanding its relationships with impact factors is important to provide useful information for climate adaptation urban planning strategies. For this purpose, the geographically-weighted regression (GWR) approach is used to explore the scale effects in a mountainous city, namely the change laws and characteristics of the relationships between land surface temperature and impact factors at different spatial resolutions (30–960 m). The impact factors include the Soil-adjusted Vegetation Index (SAVI), the Index-based Built-up Index (IBI), and the Soil Brightness Index (NDSI), which indicate the coverage of the vegetation, built-up, and bare land, respectively. For reference, the ordinary least squares (OLS) model, a global regression technique, is also employed by using the same dependent variable and explanatory variables as in the GWR model. Results from the experiment exemplified by Chongqing showed that the GWR approach had a better prediction accuracy and a better ability to describe spatial non-stationarity than the OLS approach judged by the analysis of the local coefficient of determination (R2), Corrected Akaike Information Criterion (AICc), and F-test at small spatial resolution (< 240 m); however, when the spatial scale was increased to 480 m, this advantage has become relatively weak. This indicates that the GWR model becomes increasingly global, revealing the relationships with more generalized geographical patterns, and then spatial non-stationarity in the relationship will tend to be neglected with the increase of spatial resolution.


Introduction
Rapid urban expansion and population growth leads to a common phenomenon, urban heat island (UHI), where urban areas have higher air and surface temperatures than their rural surroundings [1,2].The urban heat island (UHI) effect has received unprecedented attention because of its negative influences on the healthy development of the city and the health of the urban population, such as the increase of the heat and pollution-related mortality, the decrease of the habitats' comfort, and the elevation of the mean and peak energy demands of buildings [3].
With the development of remote sensing technology, the remote sensing data is widely used to study the UHI for the ground surface because it provides a relatively cheap and rapid method of acquiring up-to-date information over a large geographical area and obtaining data from inaccessible regions, etc.In recent years, most studies on UHI have been focused on spatial distribution, impact factors, changing trends, and so on.Among the above, the relationships between the UHI and related factors become a hot subject that wins wide-ranging and deep research.Price initially put forward the triangle relationships between vegetation index, soil moisture, and radiation temperature aided by remote sensing images [4].Carlson took advantage of remote sensing data to research urbanization and its impacts on the tropical rainforest, as well as that of the regional climate influenced by urbanization [5].Streutker proposed a two-dimensional Gaussian model to estimate the relationships between UHI magnitude and impact factors [6,7].Weng derived the vegetation fraction by a spectral mixture model as an alternative indicator of vegetation abundance, and demonstrated that the land surface temperature (LST) had a slightly stronger negative correlation with the vegetation fraction than with the Normalized Difference Vegetation Index (NDVI) for all land cover types at the different spatial resolutions [8].Some scholars, such as Yuan and Weng, compared the Normalized Difference Vegetation Index (NDVI), vegetation fraction and the percent of impervious surface area as indicators of UHI effects by investigating their relationships to the UHIs quantitatively [9,10].Jusuf and Callejas probed and verified that land use types had rather more influence on the increase of land surface temperature [11,12].Marc used spatial analysis methods to assess the UHI skin temperature amplitude and its relationship to land development intensity, size, and eco-environment for the 38 most populated cities in the United States [13].Several indices as the NDVI, the Normalized Difference Built-up Index (NDBI), the Normalized Difference Water Index (NDWI), and the Normalized Difference Bareness Index (NDBaI) were utilized to indicate the coverage of the vegetation, build-up, water and bare land, respectively, and their impacts on the UHI (LST) were studied deeply [14][15][16].Ashtiani and Mirzaei used artificial neural network (ANN) technology to predict the indoor thermal conditions on an urban heat island based on the outdoor conditions recorded at the weather stations [17,18]; Lee created a neural network predictive model to predict the urban heat island intensity and obtained good results with coefficients of correlation ranging from 0.95 to 0.99 [19].
The above-mentioned studies on the quantitative analysis of the relationships between UHI and impact factors are mainly based on the global regression models, like the ordinary least squares (OLS) and the neural network predictive model.The global regression models are dependent upon the assumption that the relationships are spatially invariant in the whole study area.However, in fact, the relationships are often characterized by local changes and the global regression model may hide the important details in the spatial distribution [20,21].
Recently, as an extension of traditional standard global regression techniques, geographically-weighted regression (GWR) was developed to explore spatially-varying relationships [20,22].A few studies have tested and verified the efficacy of the GWR model for investigating spatially-varying relationships in some fields, such as climatology [23,24], urban poverty [25,26], environmental equity [27], land use/land cover [21,28], urban landscape fragmentation [29], groundwater quantity [30], and forests [31].The GWR model can effectively solve the problem of the spatial non-stationarity by embedding the spatial location information into the regression parameters to explore the relationships between the explanatory variables and the dependent variable within the range of certain space.
In order to overcome the defects of the global regression model, a limited number of scholars applied the GWR model to quantify the relationships between the UHIs and impact factors, such as environmental factors, climate factors, and land cover types, etc. Li and Buyantuyev first introduced the GWR model into the relationships between the LST and the impact factors at the same time, and the results indicated that the GWR model not only provided a better fit than the traditional ordinary least squares (OLS) model, but also provided local detailed information about the spatial variation of LST [32,33].Su utilized the GWR model to test the spatial non-stationarity of the relationships between land cover types and the LST in TaoYuan, Taiwan, and found that the urban heat island intensity estimated by the global and GWR models for TaoYuan would equate to 2.63 • C and 3.17 • C, respectively, and the urban heat island was underestimated by the global model [34].Tian used the GWR model to quantify the relationships of the LST and related environmental factors in the Heihe River catchment, China, and found that all of the GWR models had better simulation with smaller Akaike Information Criterion (AICc) and higher coefficient of determination (R 2 ), compared with the OLS method [35].
However, the limited studies on the impacts of the impact factors on the UHI are focused on comparing the relative advantages of the GWR model over other models and quantitative descriptions of the GWR model.The scale effects, namely the change laws and characteristics at different spatial scales (spatial resolutions), is rarely involved.Li had studied the spatial scale-dependent relationships of UHI and related factors mainly from the influence of the bandwidth, namely a parameter of the GWR model, but he failed to investigate the scale effects of the GWR model from spatial resolutions [32].To solve the above problems, the Soil-Adjusted Vegetation Index (SAVI) [36], the Index-based Built-up Index (IBI) and the Soil Brightness Index (NDSI) [37,38], which indicate the coverage of the vegetation, built-up, and bare land, respectively, are utilized as the explanatory variables to investigate the scale effects, namely the change laws and characteristics at different spatial scales (spatial resolutions), on the UHI in a mountainous city by using the GWR model.

Study Area
Chongqing municipality, situated in the Southwestern China, is one of the fastest developing cities in China in the past twenty years.Characterized by rugged hills, Chongqing belongs to a typical inland mountain city and is known as a "mountain city".It is also one of the three famous "stoves" in China, the average temperature of the urban area is 5 • C-8 • C higher than that in rural area, and the urban heat island effects of the study area are obvious because of the air circulation difficulty by the terrain barrier, and the high population and building density caused by rapid urbanization.We chose the core area of the main city and its surrounding expansion area as the target area, and the specific location is shown as Figure 1.However, the limited studies on the impacts of the impact factors on the UHI are focused on comparing the relative advantages of the GWR model over other models and quantitative descriptions of the GWR model.The scale effects, namely the change laws and characteristics at different spatial scales (spatial resolutions), is rarely involved.Li had studied the spatial scaledependent relationships of UHI and related factors mainly from the influence of the bandwidth, namely a parameter of the GWR model, but he failed to investigate the scale effects of the GWR model from spatial resolutions [32].To solve the above problems, the Soil-Adjusted Vegetation Index (SAVI) [36], the Index-based Built-up Index (IBI) and the Soil Brightness Index (NDSI) [37,38], which indicate the coverage of the vegetation, built-up, and bare land, respectively, are utilized as the explanatory variables to investigate the scale effects, namely the change laws and characteristics at different spatial scales (spatial resolutions), on the UHI in a mountainous city by using the GWR model.

Study Area
Chongqing municipality, situated in the Southwestern China, is one of the fastest developing cities in China in the past twenty years.Characterized by rugged hills, Chongqing belongs to a typical inland mountain city and is known as a "mountain city".It is also one of the three famous "stoves" in China, the average temperature of the urban area is 5-8 °C higher than that in rural area, and the urban heat island effects of the study area are obvious because of the air circulation difficulty by the terrain barrier, and the high population and building density caused by rapid urbanization.We chose the core area of the main city and its surrounding expansion area as the target area, and the specific location is shown as Figure 1.The study area ranges from 29°25′-29°42′N, 106°22′-106°39′E, and its area is about 845.88 km 2 .The maximum and minimum elevations are 685 m and 90 m, respectively, and the average elevation is 284 m over the study area.The Yangtze River and the Jialing River go through the area roughly from south to north and from west to east, respectively.Separated by the terrain and the rivers, naturally, the study area forms a multi-center structure.The study area ranges from 29 • 25 -29 • 42 N, 106 • 22 -106 • 39 E, and its area is about 845.88 km 2 .The maximum and minimum elevations are 685 m and 90 m, respectively, and the average elevation is 284 m over the study area.The Yangtze River and the Jialing River go through the area roughly from south to north and from west to east, respectively.Separated by the terrain and the rivers, naturally, the study area forms a multi-center structure.

Methods
Firstly, the SAVI, IBI, NDSI, and LST were inversed from the Landsat-5 Thematic Mapper images, and the relationships between UHI (LST) and impact factors (SAVI, IBI, and NDSI) were established by the GWR model after they were aggregated from 30 m to 60 m, 120 m, 240 m, 480 m, and 960 m.Then, the advantages of the GWR models over the OLS models were analyzed and compared.Lastly, the scale effects of the relationships between UHI and impact factors by using a GWR model were analyzed further.In this paper, all of the local R 2 , coefficients and residuals are interpolated by the ordinary kriging method in order to make the results more intuitive.

Image Pre-Processing
Landsat-5 Thematic Mapper (TM) images dated on 18 June 2008 were used in this study.The images were obtained from the website of the Institute of Remote Sensing and Digital Earth Chinese Academic of Sciences, which has corrected the radiometric and geometrical distortions of the images to a quality level of 2 and has re-sampled the thermal infrared band data to a spatial resolution of 30 m before delivery.The image pre-processing was mainly performed by ENVI 4.8 software (Exelis Visual Information Solutions, Boulder, CO, USA).Firstly, a total of 25 ground control points, such as road intersections, river turns, and intersections, were selected based on the high-resolution IKONOS image, and the raw images were further rectified and re-sampled with pixel sizes of 30 m by 30 m for all bands, including the thermal band.The residuals at those control points ranged from 0.25 to 1.00 pixels.Then, the surface reflectivities of six optical band images were inversed through radiometric calibration and atmospheric correction by the ENVI FLAASH tool.

Dependent Variable: LST
The Landsat-5 Thematic Mapper (TM) thermal infrared band (10.4 to 12.5 µm) data were used to inverse the LST.After the digital number (DN) value was transformed to the thermal infrared spectral radiance received by the sensor at the top of the atmosphere, namely the top-of-atmospheric (TOA) radiance [39], and then the TOA radiance was converted to surface-leaving radiance (L sur f ace ) by eliminating the impacts of the atmosphere as follows [9,40]: In Equation (1), L TOA represents the TOA radiance, ε is the surface emissivity, and the three atmospheric parameters, namely τ, L upper and L down , represent the atmospheric transmission, the upwelling radiance and the downwelling radiance, respectively.
The atmospheric parameters (ε, L upper , and L down ) were obtained based on the atmospheric correction tool developed by Barsi et al. and available at the NASA website [40,41].In addition, the surface emissivity (ε) was estimated by utilizing the classification result and the NDVI value [42].Lastly, the LST was inversed from the surface radiance by the Landsat specific estimate of the Planck curve, expressed in Equation (2) [39]: In Equation ( 2), LST is the surface temperature in Kelvin (K), for Landsat-5 TM, K1 = 607.76(W•m −2 •sr −1 •µm −1 ), and K2 = 1260.56K.

Explanatory Variables: SAVI, IBI, and NDSI
The indices SAVI, IBI, and NDSI are positively correlated to the vegetation coverage, built-up coverage, and bare land coverage in the urban area, respectively, and the three indices are widely used for the land use mapping [36][37][38].Therefore, the indices SAVI, IBI, and NDSI were utilized as indicators of the coverage of the vegetation, built-up, and bare land, respectively, to analyze their effects on the land surface temperature.
The Soil Adjusted Vegetation Index (SAVI) was put forward by Huete to eliminate the soil background and noise.Compared to the NDVI, the SAVI is considered to be more suitable for the low-vegetation coverage region, like the urban area, due to the introduction of a soil adjusted factor, and it can be expressed as follows [36]: In Equation ( 3), l is the soil adjusted factor, and its value is ranging from 0 to 1. Usually, when l value is equal to about 0.5, the soil background difference could be reduced, and the noise would be further eliminated as well; ρ N IR and ρ RED represent the surface reflectivity of near infrared band and red band, respectively.
The Index-based Built-up Index (IBI) is constructed with the indices, namely NDBI, the Modified Normalized Difference Water Index (MNDWI) and SAVI and can be described as follows (4) [37]: where ρ GREEN , ρ RED , ρ N IR , and ρ MIR represent the surface reflectivity of the green band, red band, near-infrared band, and mid-infrared band, respectively.Xu put forward the Soil Brightness Index (NDSI) to extract bare land, and the NDSI can be described as follows [38]: where ρ GREEN and ρ RED represent the surface reflectivity of the green band and red band, respectively.

Data aggregation
The LST, SAVI, NDSI, and IBI at different spatial scales are obtained by aggregation processing, and the expressions is as shown below: In Equation ( 6), X (k, l) represents the original pixel value, X (i, j) represents the aggregated pixel value, X represents LST, SAVI, NDSI, and IBI, respectively, S is the aggregated area, and M is the number of pixels of the aggregated area.

Geographically-Weighted Regression
The OLS model is a global regression model; its regression coefficients are constant, and the relationships established by the OLS model are also treated as spatially constant in the whole study area.The traditional global regression model can be expressed as follows: where y i represents the dependent variable, x k represents the explanatory variable, β 0 is the intercept, β k is the estimated coefficient, and ε stands for the random error.
The GWR model was fully described by Fotheringham et al., who extended the concept of the global regression model by adding regional parameters, and can be rewritten as Equation ( 8) [22]: where (u i , v i ) denotes the space location of the point i, β 0 (u i , v i ) represents the intercept at the point i, β k (u i , v i ) represents the local coefficients estimated for the explanatory variable x k at the point i, and ε i is the random error at point i with distribution N 0, σ 2 .
In Equation ( 8), β 0 (u i , v i ) and β k (u i , v i ) can be estimated by solving the following matrix equation [20,22]: In Equation (8), β (u i , v i ) stands for the unbiased estimates of the regression coefficient β, X and Y represent the vector for the explanatory and dependent variables respectively, and W (u i , v i ) is the weight matrix which is used to ensure that observations near to the specific point have a larger weight.The Gaussian kernel function was used to compute the weight matrix as follows [22]: where b is called as the bandwidth of the kernel function, and d ij represents the Euclidean distance between point i and j.According to the Gaussian curve, the W ij value is gradually decreasing as the distance d ij is increasing [22,43].The weight value will be set to zero if the distance is greater than the basil width of the kernel function.
In order to compare the performances of the OLS model and the GWR model, four evaluation indices, including the coefficient of determination (R 2 ), Moran's index (Moran s I), the approximate likelihood ratio test based on the F-test (F), and the corrected Akaike Information Criterion (AIC c ), were utilized in this research.
The coefficient of determination (R 2 ), a measure of goodness-of-fit, is calculated by a comparison between the estimated and the observed values, a higher R 2 means the model is more perfectible.The Moran index (Moran s I), ranging from −1 to 1, a larger absolute value implies a more significant spatial autocorrelation [31], was calculated to test the spatial autocorrelation.
The F-test was put forwarded by Fotheringham et al., the null hypothesis is that the GWR model has no improvement over the OLS model, and the F-value is described as follows [22,35]: where RSS O and DF O represent the residual sum of squares and the degree of freedoms of the OLS model respectively, while RSS G and DF G represent the residual sum of squares and the degrees of freedom of the GWR model.The corrected Akaike Information Criterion (AIC c ) is also an important statistic to compare the model performances of the OLS model and the GWR model, and can be described as follows [44]: where n refers to the number of sampling points, σ refers to the estimate of the standard deviation of the residuals, and tr (S) refers to the trace of the hat matrix.As a general rule, a model with the smallest AIC c value is the best one which is closest to reality [45].
In this paper, both GWR and OLS models were built-up in ArcGIS 10.2 software (Environmental Systems Research Institute (ESRI), Redlands, CA, United States), and all of the evaluation indices, including R 2 , Moran s I, F value and AIC c , were also computed in it.

The Derivation Results of SAVI, IBI, NDSI, LST and Classification Result of the Study Area
The derivation results of SAVI, IBI, NDSI, LST, and the classification result of the study area were derived from the images with a spatial resolution of 30 m, and are shown in Figure 2.

The Derivation Results of SAVI, IBI, NDSI, LST and Classification Result of the Study Area
The derivation results of SAVI, IBI, NDSI, LST, and the classification result of the study area were derived from the images with a spatial resolution of 30 m, and are shown in Figure 2. The classification result of the study area was obtained by using the SVM classifier and verified by visual interpretation based on the Google Earth image, the overall accuracy of the classification result is 96.2%.As we can see in Figure 2e, The Yangtze River and the Jialing River go through the area roughly from south to north and from west to east, respectively, and the urban area forms a multi-center structure separated by the terrain and the rivers, naturally.From Figure 2d-e, the urban heat island phenomenon could easily be found in the study area, namely the surface temperature of urban areas is much higher than the surface temperature of rural areas.The surface temperature of urban areas are mainly distributed in the range of 28-56 °C, the average temperature is 32.3 °C; however, the surface temperature of the rural areas are mainly distributed in range of 22-32 °C, the average temperature is 26.8 °C.
Obviously, as we can see in Figure 2a-d, there is a non-linear relationship between LST and SAVI and NDSI in water areas, the pixels in the water areas are not suitable for building a linear model.Therefore, only the pixels in the non-water areas will be selected as the training samples in this paper.

Data Sampling
In order to satisfy the requirement of the GWR module in ArcGIS 10.2 (Environmental Systems Research Institute (ESRI), Redlands, CA, United States), the variables had to be further converted to vector formats, the conversion of data formats involved following steps: firstly, the training samples were uniformly generated in the non-water areas; then, training samples were converted into shapefile format by using the GeoDA0.9.5-i (Luc Anselin, Champaign, IL, USA).The number of training samples at different spatial resolutions is shown in Table 1.The classification result of the study area was obtained by using the SVM classifier and verified by visual interpretation based on the Google Earth image, the overall accuracy of the classification result is 96.2%.As we can see in Figure 2e, The Yangtze River and the Jialing River go through the area roughly from south to north and from west to east, respectively, and the urban area forms a multi-center structure separated by the terrain and the rivers, naturally.From Figure 2d,e, the urban heat island phenomenon could easily be found in the study area, namely the surface temperature of urban areas is much higher than the surface temperature of rural areas.The surface temperature of urban areas are mainly distributed in the range of 28 • C-56 • C, the average temperature is 32.3 • C; however, the surface temperature of the rural areas are mainly distributed in range of 22 • C-32 • C, the average temperature is 26.8 • C.
Obviously, as we can see in Figure 2a-d, there is a non-linear relationship between LST and SAVI and NDSI in water areas, the pixels in the water areas are not suitable for building a linear model.Therefore, only the pixels in the non-water areas will be selected as the training samples in this paper.

Data Sampling
In order to satisfy the requirement of the GWR module in ArcGIS 10.2 (Environmental Systems Research Institute (ESRI), Redlands, CA, United States), the variables had to be further converted to vector formats, the conversion of data formats involved following steps: firstly, the training samples were uniformly generated in the non-water areas; then, training samples were converted into shapefile format by using the GeoDA0.9.5-i (Luc Anselin, Champaign, IL, USA).The number of training samples at different spatial resolutions is shown in Table 1.

The Analysis of the Relationships between the LST and the Impact Factors at a Single Scale
In this section, the GWR was used to investigate the spatial non-stationarity relationships between the LST and the impact factors at a single spatial resolution (30 m).The impact factors include the SAVI, IBI, NDSI, and the combination of those three indices.For reference, the corresponding OLS models with the same dependent variable and explanatory variables in the GWR models were also explored.

Comparisons between the OLS and GWR Models
The GWR and OLS models at a spatial resolution of 30 m between the LST and the impact factors were established in ArcGIS 10.2 software (Environmental Systems Research Institute (ESRI), Redlands, CA, United States), and their evaluation indices are shown in Table 2.It is evident from Table 2 that all GWR models had much smaller AICc values, higher adjusted R 2 values, and better F-tests when compared with the corresponding OLS models.In single-factor analysis, the IBI had the greatest impact on the LST with an AICc value of 30,968 and adjusted R 2 value of 0.58 for the OLS model, and with an AICc value of 28,800 and adjusted R 2 of 0.70 for the GWR model, followed by SAVI and NDSI.The multi-factors relationship model, established with the IBI, SAVI, and NDSI simultaneously, provided the best performance with an AICc value of 28,800 and adjusted R 2 value of 0.73 for the GWR model, and with an AICc value of 30,379 and adjusted R 2 value of 0.62 for the OLS model.Additionally, the F-tests of the significance of improvement suggested the GWR performed better in all cases.It can be judged by the above comprehensive comparisons that, as a local statistical model, the GWR model's prediction accuracy is better than OLS and other traditional global models, which also verified the results of Li and Tian's research [32,35].

Spatial Non-Stationarity among Relationships
The results of the local coefficients and the local R 2 provide an effective measure to analyze the spatial varying relationships between the LST and the indices, which are shown in Figures 3 and 4.

Spatial Non-Stationarity among Relationships
The results of the local coefficients and the local R 2 provide an effective measure to analyze the spatial varying relationships between the LST and the indices, which are shown in Figures 3 and 4. It can be found from Figure 3 that the local coefficients for all the GWR relationship models with the LST to the indices change with the spatial location.Figure 3a-c and (d-f) are the coefficients for the SAVI, IBI, and NDSI in multi-factor models, respectively.Moreover, the yellow area represents water.
It can be found from Figure 3 that the local coefficients for all the GWR relationship models with the LST to the indices change with the spatial location.Figure 3a-c indicate the space distribution of the slopes for the SAVI, IBI, and NDSI in the single-factor models, respectively, and we can see that the SAVI is negatively correlated to the LST with coefficient values from −19.5 to −1.38, while the IBI and NDSI are positively correlated to the LST with coefficient values from 8.48 to 36.36, and from 10.53 to 48.39, respectively.The magnitude of the absolute value of the model coefficient indicates the degree of the effect of the explanatory variables on the LST.By the coefficient absolute value, we can also find that the spatial distribution of the SAVI' effect on the LST is similar to that of the IBI' effect on the LST, namely the SAVI and IBI had higher influence on the LST in the west and southeast, approximately along the distribution of the Zhongliang mountain and Tongluo mountain, and had lower influence on the LST in the other regions; while the influence of NDSI on LST is small in the west, southeast, and also the middle area.Moreover, Figure 3d-f represent the coefficients for the SAVI, IBI, and NDSI in the multi-factors GWR model (using the SAVI, IBI, and NDSI as the explanatory variables simultaneously), respectively.From Figure 3d-f, it can be seen that the coefficients for the SAVI, IBI and NDSI are not only changing with the spatial location, but also have both negative and positive correlations to the LST, with coefficient values from −24.07 to 5.20, from −3.51 to 28.86 and from −41.86 to 27.29, respectively.can also find that the spatial distribution of the SAVI' effect on the LST is similar to that of the IBI' effect on the LST, namely the SAVI and IBI had higher influence on the LST in the west and southeast, approximately along the distribution of the Zhongliang mountain and Tongluo mountain, and had lower influence on the LST in the other regions; while the influence of NDSI on LST is small in the west, southeast, and also the middle area.Moreover, Figure 3d-f represent the coefficients for the SAVI, IBI, and NDSI in the multi-factors GWR model (using the SAVI, IBI, and NDSI as the explanatory variables simultaneously), respectively.From Figure 3d-f, it can be seen that the coefficients for the SAVI, IBI and NDSI are not only changing with the spatial location, but also have both negative and positive correlations to the LST, with coefficient values from −24.07 to 5.20, from −3.51 to 28.86 and from −41.86 to 27.29, respectively.Local R 2 , ranging from 0 to 1, were used to detect and investigate the ability of the GWR model to fit the observations, a local model with higher local R 2 value will have a better performance.The spatial distributions of the local R 2 in the GWR models are shown as Figure 4. Unlike the OLS model, the spatial patterns of the local R 2 in the GWR model represent a marked regional differentiation.In single-factor GWR models, the local R 2 values for the LST-SAVI model, LST-IBI model, and LST-NDSI model are in the ranges of 0.016-0.737,0.165-0.757,and 0.02-0.71,respectively.Moreover, the local R 2 value for the multi-factor GWR model is from 0.198 to 0.774.For all GWR models, the spatial distributions of the local R 2 are similar, characterized by higher values (>0.5) in the north, southeast, and southwest part.In conclusion, both the coefficients and the local R 2 for all the GWR models are changing with the spatial location, and this indicates that the GWR has a good ability to characterize the non-stationarity of the relationships between the LST and the indices.Local R 2 , ranging from 0 to 1, were used to detect and investigate the ability of the GWR model to fit the observations, a local model with higher local R 2 value will have a better performance.The spatial distributions of the local R 2 in the GWR models are shown as Figure 4. Unlike the OLS model, the spatial patterns of the local R 2 in the GWR model represent a marked regional differentiation.In single-factor GWR models, the local R 2 values for the LST-SAVI model, LST-IBI model, and LST-NDSI model are in the ranges of 0.016-0.737,0.165-0.757,and 0.02-0.71,respectively.Moreover, the local R 2 value for the multi-factor GWR model is from 0.198 to 0.774.For all GWR models, the spatial distributions of the local R 2 are similar, characterized by higher values (>0.5) in the north, southeast, and southwest part.In conclusion, both the coefficients and the local R 2 for all the GWR models are changing with the spatial location, and this indicates that the GWR has a good ability to characterize the non-stationarity of the relationships between the LST and the indices.

Scale Effects of the Relationships Based on the GWR model at Multiple Scales
This section mainly analyzed the spatial scale effects of the relationships between LST and multi-factors at multi-scales.Firstly, the SAVI, IBI, NDSI, and LST images were aggregated from the original spatial resolution (30 m) to different spatial resolutions (60 m, 120 m, 240 m, 480 m, and 960 m).Then, the multi-factor relationships between LST and SAVI, IBI, and NDSI were built up based on a GWR model at different spatial resolutions (spatial scales).The maps of the local R 2 and the residuals of the multi-factor GWR models at different spatial scales are shown in Figures 5 and 6, respectively.

Scale Effects of the Relationships Based on the GWR model at Multiple Scales
This section mainly analyzed the spatial scale effects of the relationships between LST and multifactors at multi-scales.Firstly, the SAVI, IBI, NDSI, and LST images were aggregated from the original spatial resolution (30 m) to different spatial resolutions (60 m, 120 m, 240 m, 480 m, and 960 m).Then, the multi-factor relationships between LST and SAVI, IBI, and NDSI were built up based on a GWR model at different spatial resolutions (spatial scales).The maps of the local R 2 and the residuals of the multi-factor GWR models at different spatial scales are shown in Figures 5 and 6, respectively.It can be seen from Figure 5 that the local R 2 value of the entire study area in general is increasing with the increase of the spatial scale (Figure 5a-f).The spatial distributions of the local R 2 in the GWR models at small spatial scales (30 m and 60 m) are basically similar (Figure 5a,b), namely, the relatively low local R 2 value is mainly located in some area with large topographic variation and the intersection area of the Jialing River and the Yangtze River, and the local R 2 of other region is relatively high.Among them, the local R 2 values of the intersection regions of the Jialing River and Yangtze River, where the center region of the city with relatively higher human activities is, are obviously lower than the local R 2 values in other regions.This may be due to the influence of the higher human activities, which have complicated impacts on the SUHI.Therefore, further consideration needs to be given to the radiation effects of human activities in future work.
On the other hand, we can also intuitively find that the local difference of the spatial distribution of the local R 2 is gradually reduced with the increase of the spatial scale, especially after the spatial scale is decreased to the 960 m; this change is much more obvious and the local R 2 values of the entire study area tends to be consistent.Further statistics from Figure 5 show that with the increase of spatial scale, the distribution range of local R 2 changes from 0.19-0.77(30 m) to 0.83-0.92(960 m), the mean local R 2 increases from 0.53 (30 m) to 0.86 (960 m), the proportion of the area with high local R 2 value (>0.65) to the total area decrease from 21.5% (30 m) to 96.2% (960 m).However, the standard deviation of local R 2 decreases from 0.22 (30 m) to 0.03 (960 m), which implies that the GWR results become increasingly global, revealing the relationships with more generalized geographical patterns, and then spatial non-stationarity in the relationship tends to be neglected with the increase of the spatial resolution.
Furthermore, it is apparent from Figure 6 that with the increase of spatial scale, the degree of the spatial autocorrelation of the local residuals of the GWR model is gradually enhanced.Among them, the residuals distribution of the multi-factors GWR model at 30 m is approximately a random distribution (Figure 6a), the residuals distribution of the GWR model at 240 m appears to be a certain clustering (Figure 6d), and at the spatial resolution of 960 m, the values of the residuals are getting are closer and closer (Figure 6f).
The statistical results of residuals for GWR models at different spatial scales are shown in Table 3.It can be found from the Table 3 that with the increase of the spatial scale, the variation range of residuals is reduced from −2.78-4.10(30 m) to −1.38-1.25 (960 m), and the standard deviation of residuals is decreased from 1.83 (30 m) to 0.52 (960 m).According to the change trends of the distribution range and standard deviation of the residuals, the GWR models yield less accurate simulation results which reveal more generalized geographical patterns with the increase of the spatial scale, which is consistent with the analysis results from Figure 5.At the same time, from Table 3, the Moran's I values of the residuals gradually increase from 0.19 (30 m) to 0.39 (960 m) with the increase of the spatial scale.Among them, the Moran's I values at 30 m, 60 m, and 120 m are very close, ranging from 0.19-0.21,and the Moran's I values at 240 m, 480 m, and 960 m range from 0.31 to 0.39.According to the analysis of Moran's I values, with the increase of the spatial scale, the residuals of the GWR model show more and more strong spatial autocorrelation and spatial clustering.Furthermore, several performance indicators, including AICc, adjusted R 2 , and Moran's I, of GWR and OLS models at different spatial scales are shown in Table 4.  4. It can be seen from Table 4 that with the increase of the spatial scale, the AICc values of the GWR and OLS models are continuously decreasing, and the adjusted R 2 are continuously increasing.This indicates that, characterized by AICc and adjusted R 2 , the fitting degree of the GWR and OLS models are both constantly improved with the increase of spatial scale.
At the same time, we also can find from Table 4 that, characterized by higher R 2 values and lower AICc values, the GWR models have better ability than the OLS models to explain the relationships between SUHI and the impact factors at small spatial scales (30 m-240 m), and when the spatial scale is increased to 480 m and 960 m, this advantage becomes relatively weak.The main reason is, probably, that the distance of the adjacent pixels, which has a positive correlation with the bandwidth of the GWR model, is increasing with the increase of the spatial resolution, and the GWR results become increasingly global, revealing the relationships with more generalized geographical patterns, and then spatial non-stationarity in the relationship will tend to be neglected.When the scale space is small, the GWR model can better reveal the SUHI' local variation details, but when the spatial scale becomes large, the performance difference between the GWR model and OLS model is relatively small because of the decreases of the SUHI' local difference details.Therefore, if the spatial resolution of remote sensing data is less than 240 m, GWR mode is recommended to be used in the monitoring and analysis of SUHI in the mountain city, and when the spatial resolution is greater than 480 m, both the GWR and OLS models are suitable for the studies of SUHI because there are few performance differences between them.The GWR model has a relatively high fitting degree at large spatial scale with the cost of the smoothing of the detailed differences of the neighbor region; this is not conducive to reflecting the real information of the research area.Thus, we need to pay more attention to weigh the balance between the fitting degree of the model and the level of details of the spatial geographic information reflected by the model.

Discussion
The results reported in this paper showed significant spatial non-stationarity in the relationships between the land surface temperature (LST) and three explanatory variables in the mountain city.Previous studies have reported that the regression relationships between land surface temperature and impact factors are characterized by significant spatial non-stationarity and scale-dependence.Li found that the regression relationships between land surface temperature and eight explanatory variables are characterized by significant spatial non-stationarity [32].Cui found that the variations in temperature and its correlation with NDVI in different climatic zones and land use types were characterized by significant spatial non-stationarity [46].Thus, the hypothesis that a single set of regression coefficients could capture space-varying relationships between variables in large-scale analysis is not reasonable [47].Some other scholars have made much the same discovery that significant spatial non-stationarity is widespread and routine in the relationships between the land surface temperature (LST) and its impact factors [33][34][35].In this paper, the results from both single-factor and multi-factor models at spatial resolution of 30 m were used to analyze the spatial non-stationarity in the relationships between LST and three indices.The results from the single-factor models showed that the SAVI was negatively correlated to the LST and the IBI and NDSI were positively correlated to the LST.Regardless of the positive or negative association of LST with SAVI, IBI, and NDSI, the associations were stronger for the west and southeast areas (Figure 3a-c), approximately along the distribution of the two mountains, where the land surface temperature distribution was more sensitive to two mountain areas with high vegetation coverage and high altitude (Figure 4a-c).Meanwhile, the results from the multi-factors models showed that the coefficients for the SAVI, IBI, and NDSI were not only changing with the spatial location, but also had both negative and positive correlations to the LST (Figure 3d-f).However, the associations were stronger for south areas of the urban under the interaction of SAVI, IBI, and NDSI, where the land surface temperature distribution started to be more sensitive to some urban areas (Figure 4d).This indicated that the performance of the GWR model could be improved by using the combination of SAVI, IBI, and NDSI together.This result corresponded to previous findings that a multi-factors GWR model may be more appropriate if the explanation variables show spatial stationarity by the test [22,48].
Meanwhile, we found that the correlations between LST and SAVI, IBI and NDSI were not significant in some regions, especially the intersection regions of the Jialing River and Yangtze River, where the center region of the city is, with a relatively higher human activity.Previous studies suggested that economic development, population growth, and industry expansion have been proved as the significant driving factors for changes in the urban thermal environment [49][50][51].Zhao suggested that the correlation between the temperature and the NDVI was not significant in the regions with higher human activity and wet conditions [31].
In recent years, the detection of scale-dependent phenomenon and modelling of ecological processes across different scales have grown significantly [32,52].That research studied the spatial scale-dependent relationships of UHI and impact factors mainly from the influence of the bandwidth, namely a parameter of the GWR model.Li pointed out that by changing the bandwidth size of the GWR model, a series of spatial patterns of regression parameters at different spatial scales can be obtained, which are helpful in determining the operational scale range for each explanatory variable [32].However, they failed to investigate the scale effects of the GWR model from the spatial resolutions.In our study, we investigated the spatial scale-dependent relationships of UHI and impact factors from another aspect, namely investigating the scale effects of the GWR model from spatial resolutions, and we found that there was a contradictory phenomenon that the overall fitting degree of GWR model was gradually improved, but the local difference of the spatial distribution of the local R 2 was gradually reduced and the residuals of the GWR model showed stronger spatial autocorrelation and spatial clustering with the increase of spatial resolution.This indicated that the performance of the GWR model was not improved with the increase of spatial resolution.Further study on the choice of the optimal spatial resolution for the GWR model must be the emphasis of the future works.
In our study, the results showed that the GWR model is superior to the OLS model with better model performance and lower spatial autocorrelation of residuals at small spatial scales (30 m-240 m).Several scholars suggested that the GWR model was advantage over other models.Su utilized the GWR model to test the relationships between land cover types and the LST, and found that the urban heat island intensity estimated by the global and GWR models for TaoYuan would equate to 2.63 • C and 3.17 • C, respectively, and the urban heat island was underestimated by the global model [34].Tian used a GWR model to quantify the relationships of the LST and related environmental factors and found that all GWR models had better simulation with smaller AICc value and higher R 2 value, compared with OLS method [35].
However, when the spatial scale is increased to 480 m and 960 m, this advantage becomes relatively weak.One reason for this may be that the distance of the adjacent pixels, which has a positive correlation with bandwidth of the GWR model, is increasing with the increase of the spatial resolution, and the GWR results become increasingly global, revealing the relationships with more generalized geographical patterns, and then spatial non-stationarity in the relationship will tend to be neglected.Therefore, if the spatial resolution of remote sensing data is less than 240 m, GWR mode is recommended to be used in the monitoring and analysis of SUHI in the mountain city, and when the spatial resolution is greater than 480 m, both GWR and OLS models are suitable for the studies of SUHI because there are few performance differences between them.
In this paper, we have explored and analyzed the scale effects in the relationships between LST and SAVI, IBI, and NDSI in the mountainous city, but the result of this research still has some limitations, like lacking of the analysis of the scale effects of the relationships in non-mountainous cities, and lacking of the consideration of the influences of the complex topography in the study area.Therefore, this paper only provides a limited reference opinion for helping researchers to choose a suitable model for the non-mountainous cities to provide useful information for climate adaptation urban planning strategies.
In future research, several works are needed to be further studied.Due to the limitations of the complexities of UHI impact factors and the data, the conclusions of this paper just reflect the relationships between UHI and several surface indices in a mountainous city.Therefore, the comprehensive effects of the underlying surface properties, terrains, human activities, and other impact factors on the UHI of both mountainous and non-mountainous city will be the emphasis of future works.

Conclusions
In this paper, the SAVI, IBI, and NDSI were extracted from Landsat-5 TM images, and the LST was inversed by a thermal radiation transfer model to analyze the distribution feature of SUHI.Then, the regression models of the LST and SAVI, IBI, and NDSI were built at different spatial scales using the geographically-weighted regression approach, and the spatial non-stationary and the spatial scale effects of the relationships between the LST and the indices at both single scale and multi-scales were discussed quantitatively.Chongqing was selected to experiment for a case study, and the results are as follows: 1.
Both single-factor and multi-factors GWR models have better prediction accuracies, characterized by much smaller AICc values, higher adjusted R 2 values, and better F-tests, when compared with the corresponding OLS models.At the same time, both the coefficients and the local R 2 of the GWR models are changing with the spatial location, and this indicates that the GWR has a good ability to characterize the non-stationarity of the relationships between the LST and the indices.

2.
With the increase of spatial scales, the overall fitting degree of the GWR model is gradually improved based on the distribution range, the mean value of local R 2 .However, the standard deviation of the local R 2 and residuals are gradually reduced from 0.22 (30 m) to 0.03 (960 m) and 1.83 (30 m) to 0.52 (960 m), respectively.Meanwhile, the Moran's I values of the residuals gradually increase from 0.19 (30 m) to 0.39 (960 m).This indicates the GWR model becomes increasingly global, revealing the relationships with more generalized geographical patterns, and then spatial non-stationarity in the relationship tends to be neglected with the increase of the spatial resolution.

3.
Characterized by higher R 2 value and lower AICc value, GWR models have better ability than OLS models to explain the relationships between SUHI and impact factors (the SAVI, IBI, and NDSI) at small spatial scales (30 m-240 m), and when the spatial scale is increased to 480 m and 960 m, this advantage has become relatively weak because the GWR model becomes increasingly global, revealing the relationships with more generalized geographical patterns, and then spatial non-stationarity in the relationship tends to be neglected with the increase of the spatial resolution.Therefore, if the spatial resolution of remote sensing data is less than 240 m, GWR mode is recommended to be used in the monitoring and analysis of SUHI in the mountain city, and when the spatial resolution is greater than 480 m, both GWR and OLS models are suitable for the researches of SUHI because there are few performance differences between them.

Figure 1 .
Figure 1.Location of the study area ((a) China, (b) Chongqing, and (c) the study area).

Figure 1 .
Figure 1.Location of the study area ((a) China, (b) Chongqing, and (c) the study area).

Figure 2 .
Figure 2. The derivation results of SAVI, IBI, NDSI, LST, and the classification result of the study area at the spatial resolution of 30 m ((a) SAVI; (b) NDSI; (c) IBI; (d) LST; and (e) the classification result of the study area).

Figure 2 .
Figure 2. The derivation results of SAVI, IBI, NDSI, LST, and the classification result of the study area at the spatial resolution of 30 m ((a) SAVI; (b) NDSI; (c) IBI; (d) LST; and (e) the classification result of the study area).

Figure 3 .
Figure 3.The spatial distribution of the local coefficients estimated by GWR models at a spatial resolution of 30 m. (a-c) are the slope coefficients for the SAVI, IBI, and NDSI in single-factor models, respectively;and (d-f) are the coefficients for the SAVI, IBI, and NDSI in multi-factor models, respectively.Moreover, the yellow area represents water.
indicate the space distribution of the slopes for the SAVI, IBI, and NDSI in the single-factor models, respectively, and we can see that the SAVI is negatively correlated to the LST with coefficient values from −19.5 to −1.38, while the IBI and NDSI are positively correlated to the LST with coefficient values from 8.48 to 36.36, and from 10.53 to 48.39, respectively.The magnitude of the absolute value of the model coefficient indicates the degree of the effect of the explanatory variables on the LST.By the coefficient absolute value, we

Figure 3 .
Figure 3.The spatial distribution of the local coefficients estimated by GWR models at a spatial resolution of 30 m. (a-c) are the slope coefficients for the SAVI, IBI, and NDSI in single-factor models, respectively;and (d-f) are the coefficients for the SAVI, IBI, and NDSI in multi-factor models, respectively.Moreover, the yellow area represents water.

Figure 4 .
Figure 4.The spatial distribution of the local R 2 in the GWR models at a spatial resolution of 30 m, (a) for the LST-SAVI model; (b) for the LST-IBI model; (c) for the LST-NDSI model; and (d) for the multifactor model (using the SAVI, IBI, and NDSI as the explanatory variables simultaneously).Moreover, the yellow area represents water.

Figure 4 .
Figure 4.The spatial distribution of the local R 2 in the GWR models at a spatial resolution of 30 m, (a) for the LST-SAVI model; (b) for the LST-IBI model; (c) for the LST-NDSI model; and (d) for the multi-factor model (using the SAVI, IBI, and NDSI as the explanatory variables simultaneously).Moreover, the yellow area represents water.

Figure 5 .
Figure 5.The spatial distribution of the local R 2 in the multi-factor GWR models (using the SAVI, IBI, and NDSI as the explanatory variables simultaneously) at different spatial resolutions, (a-f) represent the spatial resolutions of 30 m, 60 m, 120 m, 240 m, 480 m, and 960 m, respectively.Moreover, the yellow area represents water.

Figure 5 .
Figure 5.The spatial distribution of the local R 2 in the multi-factor GWR models (using the SAVI, IBI, and NDSI as the explanatory variables simultaneously) at different spatial resolutions, (a-f) represent the spatial resolutions of 30 m, 60 m, 120 m, 240 m, 480 m, and 960 m, respectively.Moreover, the yellow area represents water.

Figure 6 .
Figure 6.The spatial distribution of the local residuals in the multi-factor GWR models (using the SAVI, IBI, and NDSI as the explanatory variables simultaneously) at different spatial resolutions; (a)-(f) represent the spatial resolutions of 30 m, 60 m, 120 m, 240 m, 480 m, and 960 m, respectively.Moreover, the blue area represents water.

Figure 6 .
Figure 6.The spatial distribution of the local residuals in the multi-factor GWR models (using the SAVI, IBI, and NDSI as the explanatory variables simultaneously) at different spatial resolutions; (a)-(f) represent the spatial resolutions of 30 m, 60 m, 120 m, 240 m, 480 m, and 960 m, respectively.Moreover, the blue area represents water.

Table 1 .
The number of training samples at different spatial resolutions.

Table 2 .
The performance comparisons of the OLS and GWR models.

Table 3 .
The residuals statistical of the GWR models at different spatial scales.

Table 4 .
Performance indicators of GWR and OLS models at different spatial scales.