Analysis of the Urban Heat Island Effect in Shijiazhuang , China Using Satellite and Airborne Data

The urban heat island (UHI) effect resulting from rapid urbanization generally has a negative impact on urban residents. Shijiazhuang, the capital of Hebei Province in China, was selected to assess surface thermal patterns and its correlation with Land Cover Types (LCTs). This study was conducted using Landsat TM images on the mesoscale level and airborne hyperspectral thermal images on the microscale level. Land surface temperature (LST) was retrieved from four scenes of Landsat TM data in the summer days to analyze the thermal spatial patterns and intensity of surface UHI (SUHI). Surface thermal characteristics were further examined by relating LST to percentage of imperious surface area (ISA%) and four remote sensing indices (RSIs), the Normalized Difference Vegetation Index (NDVI), Universal Pattern Decomposition method (VIUPD), Normalized Difference Built-up Index (NDBI) and Biophysical Composition Index (BCI). On the other hand, fives scenes of airborne TASI (Thermal Airborne Spectrographic Imager sensor) images were utilized to describe more detailed urban thermal characteristics of the downtown of Shijiazhuang city. Our results show that an obvious surface heat island effect OPEN ACCESS Remote Sens. 2015, 7 4805 existed in the study area during summer days, with a SUHI intensity of 2–4 °C. The analyses reveal that ISA% can provide an additional metric for the study of SUHI, yet its association with LST is not straightforward and this should a focus in future work. It was also found that two physically based indices, VIUPD and BCI, have the potential to account for the variation in urban LST. The results concerning on TASI indicate that diversity of impervious surfaces (rooftops, concrete, and mixed asphalt) contribute most to the SUHI, among all of the land cover features. Moreover, the effect of impervious surfaces on LST is complicated, and the composition and arrangement of land cover features may play an important role in determining the magnitude and intensity of SUHI. Overall, the analysis of urban thermal signatures at two spatial scales complement each other and the use of airborne imagery data with higher spatial resolution is helpful in revealing more details for understanding urban thermal environments.


Introduction
LCTs have played a significant role in environmental change, and accurate and real-time information of LCTs is critical to environmental monitoring and management [1,2].The interaction between LCTs and climatic variability is so complex that advanced models are required to describe this process, especially for urban-rural areas.Since the change from nature areas to impervious areas can result in significant regional climate changes, urban thermal environments have been identified in some studies from the perspective of the scientific community for meteorology and remote sensing [3][4][5].The differences in surface temperature between urban areas and surrounding rural areas, a phenomenon termed as urban heat island (UHI) effect and due to urbanization, have been widely observed [6,7].UHI considerably increases the energy demands and degrades the outdoor air quality especially in summertime, as accompanying the negative impact on the local mseteorology conditions.
The utilization of remote sensing in the assessment of surface thermal properties as well as surface UHI has been carried out by various researchers.These studies can be conducted using satellite sensors such as Moderate Resolution Imaging Spectroradiometer (MODIS) and Advanced Very High Resolution Radiometer (AVHRR), which have been found to be feasible to describe the coarse thermal spatial distributions, yet, insufficient to investigate the accurate relationships between LST and surface characteristics [8][9][10][11].On the other hand, satellite sensors of medium spatial resolution, such as Landsat Thematic Mapper (TM)/Enhanced Thematic Mapper Plus (ETM+) thermal infrared data and ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer), have been extensively utilized for studies of SUHI at the local scale [12][13][14].For example, [10] used ETM images to assess the thermal environment of major cities in Greece; Chen, et al. [15] examined the changes in LCTs and the associated intensity and spatial pattern of the SUHI effect in a rapidly changing area of the Pearl River Delta region in China.
Although satellite data are useful in analyzing SUHI effect at a coarse or medium scale, they cannot lend themselves to well establish a better understanding of what specific LCTs contribute to the SUHI effect.Analysis of thermal energy characteristics for discrete LCTs of the urban landscape generally requires being carried out at a very fine spatial scale, in order to adequately resolve these land cover features and their attendant thermal energy regimes.Compared with space-borne data mostly of medium spatial resolution, high spatial resolution thermal data obtained from aircraft can provide more detailed information and be used as an alternative approach to assess the SUHI effect carefully and profoundly [16,17].Gluch, et al. [18] and Pu, et al. [19] for example used airborne thermal data to examine the thermal response of individual land covers in urban areas.According to this context, high spatial resolution thermal data is necessary to investigate and reveal to what extent typical LCTs within the city contribute to increases in urban thermal emissions or vice versa.
Impervious surface area (e.g., asphalt, building and paved roads), featured by impeding the surface water infiltration, is often used as a major indicator of the urbanization and environmental quality [20,21].Percentage imperviousness is reported to be one significant factor to increase heat emissions in urban areas, further justifying its practical application in assessing the urban thermal environment, since it can capture more spatial information on complex urban landscapes [22][23][24].Although the positive correlation between ISA and urban LST has been known for a while, their roles and interactions across complex urban environments remain elusive.Therefore, a thorough understanding of the effect of the biophysical parameter of ISA% on LST needs to be further investigated, in order to effectively adapt to the urban land cover changes and climate changes.It was reported that thermal responses of urban landscapes vary with different biophysical properties especially of vegetation covers and building [4].RSIs, which are largely determined by land cover types, can be used to characterize surface thermal patterns as well as SUHI [25,26], although these indices have both strengths and weaknesses.RSIs such as Normalized Difference Vegetation Index (NDVI), Normalized Difference Built-up Index (NDBI), and the Normalized Difference Water Index (NDWI) have been employed to establish their relationships with LST using regression analysis [15,27,28].However, most of these relationships are subject to variation due to the predominantly bare ground surfaces or undestroyed vegetation covers [24,29].Simple remote sensing indices alone may not be sufficient to quantitatively determine SUHI, and therefore more physically based indices need to be investigated to formulate the correlations between surface parameters and LST.
To date, studies on the quantitative relationship between SUHI and typical LCTs based on multi-scale thermal remote sensing observations have not been publicly reported.In this study, a combination of aerial thermal data and spaceborne thermal data was used to investigate how LCTs influences the spatial pattern and intensity of the SUHI with a case study in Shijiazhuang, China.Analyses focusing on the mesoscale level with Landsat TM images and on the microscale level with airborne hyperspectral thermal data, respectively, were conducted to comprehensively discover urban thermal patterns in terms of LCTs, RSIs, ISA% and LST.Thermal sensors of the two different spatial resolutions were used to improve urban thermal study.The specific objectives of this research are: (1) analyzing the spatial patterns as well as the intensity of SUHI using four Landsat ETM+ images acquired in the summer of 2006, 2007, 2009 and 2010; (2) quantifying the correlations between LST and some indices (NDVI, VIUPD, NDBI, BCI, and ISA%) over the study region; and (3) discovering thermal values over individual land cover and capturing the spatial distribution of urban thermal patterns, using the fine spatial resolution surface temperature derived from airborne thermal images (TASI) collected for the downtown of Shijiazhuang.

Study Site and Datasets
Shijiazhuang is a city located in North China (Lat.38°03′N, Lon.114°29′E).Situated in the low latitudes of Eurasia and having a semiarid monsoon climate, its annual average precipitation is only about 500 mm/year, mostly occurring in July and August.Shijiazhuang is the largest city of the Hebei Province with a population of around 5 million.In this research, the city of Shijiazhuang and its suburban area were selected as the study area (Figure 1).The dominant land cover in this region is dense evergreen forest, agriculture, bare soil and diverse urban materials.Urban materials mainly consist of roof materials, concrete and asphalt surfaces with varying age, and bare soil and water pools either in the undeveloped areas or park.Other vegetation covers include small patches of urban lawns and wetlands.Overall, the high diversity of urban materials and a mixture of various vegetation types make it an ideal area to analyze SUHI effects.
To quantitatively analyze the urban thermal patterns and explore the urban heat island intensity in the city of Shijiazhuang, four cloud-free TM images (5 September 2006;23 August 2007;12 August 2009; and 15 August 2010) from Landsat 5 were selected.Two scenes (Path/Row: 124/33, 124/34) were encompassed in the study region for each of the acquisition date.Landsat images provide six visible and near infrared bands (1-5 and 7) with a spatial resolution of 30 m, and a thermal infrared band (band 6) with a spatial resolution of 120 m.Five high-resolution airborne images collected by a Thermal Airborne Spectrographic Imager (TASI) sensor were also acquired, as part of the National High Technology Research and Development Program of China (863 Program), to capture the urban thermal spatial distribution at the microscale level.The basic specifications of TASI are given in Table 1 [30,31].These TASI images were acquired at an altitude of 0.5 or 1 km from 25 July 2010 to 15 August 2010 including three time phases (morning, noon, and evening).The coverage of these airborne thermal images is centered at the downtown of Shijiazhuang, extending 10 km in north-south and 5 km in east-west.In addition, a WordView-2 (WV2) image was acquired in September 2010.WV2 provides a higher spatial (2 m MS and 0.5 m Panchromatic) resolution.A short-wave infrared image (0.9-2.5 um, 101 bands) collected by a Short-Wave-Infrared Airborne Spectrographic Imager (SASI) sensor on 15 August 2010, was used as ancillary data to identify specific landscape features.All data used in this study are listed in Table 2.

Methodology
As previously mentioned, the current study was conducted using Landsat TM images on the mesoscale level and airborne hyperspectral thermal images on the microscale level.The overall processing scheme is displayed in Figure 2.

Image Pre-Processing
The multi-temporal Landsat TM images were first geo-referenced to the Universal Transverse Mercator (UTM) coordinate system based on the high-resolution Google Earth TM images [32] and WV imagery, and then re-sampled using the nearest neighbor algorithm with a pixel size of 30 m × 30 m for all bands, including the thermal band.The root-mean-square error (RMSE) of the rectification was less than one pixel in this study.Atmospheric correction was applied to the visible and near infrared bands using the Fast Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) module in the ENVI 4.8 software.A mosaic of two Landsat TM images acquired during each day was compiled.The mosaic images were finally clipped to the study area.Assessing the accuracy of surface reflectance retrievals is significant for the calculation of various remote sensing indexes.In this study, we validated the obtained surface reflectance via comparisons with limited ground based reflectance observations.During the Landsat-5 overpass on 15 August 2010, we collected 15 ground measurements of surface reflectance spectra of dominant cover types using the ASD (Analytical Spectral Devices) spectroradiometer (covering the spectrum of 0.35-2.5 um).More detailed information of this ground measurements can be found in our previous study in [31].By comparing the reflectance retrieval from TM with the ASD measurements aggregated using the Landsat 5 TM spectral response functions, it was demonstrated that the reflectance of all six bands agrees well.The absolute errors in terms of root mean squares error are 0.044 for the infrared bands and 0.018 for the visible bands.Given the uncertainty of the ASD measurements, aerosol scattering model, and slight surface heterogeneity over the test sites, the accuracy of the surface reflectance retrieval is reasonable.Based on the assumption that the radiometric response function is linear, TASI images were calibrated using coefficients retrieved by measuring two known standard black-body sources.Atmospheric correction on the basis of the radiative transfer calculation was performed to transform the radiance into surface radiance.The atmospheric upwelling and downwelling radiance and the atmospheric transmittance required for the atmospheric correction were obtained using the MODTRAN radiative transfer code and local atmospheric profiles in the study area.The atmospheric profiles of air temperature and humidity were retrieved from near real time NCEP reanalysis data.In addition, the in situ air temperature measurements obtained from a local weather station and the in situ atmospheric column water vapor content measured by a CE318 automatic sun tracking photometer were used to correct the profiles of air temperature and humidity.Other parameters were set as following: ozone (mid-latitude summer) and carbon dioxide (from MODTRAN simulations).The retrieved spectra of welling path radiance and downwelling sky radiance were both integrated to TASI channels on the basis of the corresponding filter function.More detailed information about the data pre-processing can be found in our previous study [33].

Retrieval of Impervious Surface Area from Landsat TM Images
Many methods have been developed in the literature for estimating ISA% based on remotely sensed imagery [34].Spectral mixture analysis, an image-processing method that determines the fraction of fundamental components, provides a suitable approach to extract ISA% [35].A physical image-processing method integrating normalized spectral mixture analysis (NSMA) [36] with multiple endmember spectral mixture analysis (MESMA) [37] was used to estimate ISA from Landsat images in this study [31].This method involves endmember selection, unmixing analysis, and calculation of fraction images.
First, the pixel purity index (PPI) approach was combined with an N-dimensional visualizer to identify the major candidate imagery endmember, referring to the vegetation-impervious surface-soil (V-I-S) model [38].Endmember spectral containing high albedo, low albedo, bare soil, and vegetation were selected from the image itself, resulting in a total of 8, 8, 7, 7 endmembers for each image of 2006, 2007, 2009, and 2010, respectively.Specially, water area was masked out and not considered in this process.Second, NSMA was applied to all endmember combinations (in total 47, 47, 35, 35 combinations for each image).Of all the endmember combinations for each pixel, the one that brought about all non-negative endmember fractions and the lowest root mean squared error was considered the optimal combination.ISA fraction was finally calculated by summing up the low albedo and high albedo fractions.
To assure the quality of the retrieved impervious surface area, the estimated result was evaluated by comparing the Landsat-based ISA% with that obtained from high spatial resolution imagery.A total of 150 sample plots with 3 × 3 pixels (i.e., 90 m × 90 m) were identified on the impervious fraction image of August 2010, while the corresponding actual ISA% was digitized from the WV2 image.
In order to investigate the inherent relationship between imperviousness density distribution and urban thermal pattern, the urban landscape was divided into different classes based on the imperviousness density.According to a rule-based analysis, pixels were classified as urban when the percent ISA were greater than 30%, whereas those pixels of less than 30% were considered non-urban.Moreover, three ISA threshold values were used to define different levels of the urban development: 30%-50% for low density; 50%-75% for medium density and >75% for high-density.

Derivation of Remote Sensing Indices from Landsat TM Images
Four remote sensing indices were utilized in this the application to demonstrate the significance in discovering the urban thermal patterns.
NDVI [39] is widely used to estimate the density of vegetation and can be expressed as: where ρ4 and ρ3 represent the reflectance value of near-infrared band (band 4) and red band (band 3), respectively, both from the Landsat images.NDBI [40] is an index that can describe the information of built-up areas.It can be noted as: ) where b5 and b4 represent the digital value of mid-infrared band (band 5) and near-infrared band (band 4), respectively, for the Landsat images.
VIUPD is a vegetation index that is computed based on the universal pattern decomposition approach [41].The index can reflect some biophysical factors, such as vegetation concentrations and the degree of terrestrial vegetation vigor, more sensitively than the NDVI and EVI.More importantly, VIUPD is sensitive to bare soil and dead vegetation coverage.Unlike the traditional broadband vegetation indices usually computed with near-infrared and red bands, VIUPD is constructed on the basis of the all observed wavelengths.It can be computed as a combination of four pattern decomposition coefficients (the patterns of vegetation, soil, water, and supplementary) and can be derived as: where Cv, C5, Cw and C4 represent the respective standard pattern reflectance of vegetation, soil, water, and supplementary yellow-leaf; a is the coefficient according to the standard soil pattern.Information of biophysical composition information is essential in delineating urban ecological morphology for environmental modeling and management.Traditional methods to obtain biophysical composition information were either dependent on various land cover classification approach or too complicated to carry out.BCI [25] is a biophysical index proposed for effectively extracting urban biophysical compositions.By means of Tasseled Cap (TC) transformation and V-I-S triangle model, the formula of BCI can be given as: where H is the representative spectral response of high albedo; L is the representative spectral response of low albedo; and V is the representative spectral response of vegetation.These three factors can be retrieved with the following formula: min max min 11 11 min max min 22 22 min max min 33 33 where TCi (i = 1, 2, and 3) are the three normalized TC components, representing high albedo, vegetation, and low albedo, respectively; TCmax and TCmin are the maximum and minimum values of the ith TC component.

Retrieval of Surface Temperature from Landsat TM Images
The derivation of LST from Landsat TM data mainly involves three steps: radiometric calibrations, atmospheric and surface emissivity corrections, and characterization of spatial variability.The digital number of geometrically corrected Landsat TM band 6 (10.44-12.42um) were first converted to spectral radiance using Equation ( 8), and then converted to at-sensor brightness temperature with Equation ( 9): 0.0056322 0.1238 where Tb is the brightness temperature in Kelvin, L is spectral radiance; K1 and K2 are the calibration constants in m•W•cm −2 •sr•μm −1 (K1 = 60.776,K2 = 1260.5).Since the brightness temperature obtained above is according to a black body, corrections for spectral emissivity were required to retrieve the surface temperature [42].The formula is as follows: where λ is the wavelength of emitted radiance (λ = 11.5 μm), σ is the Boltzmann constant (1.38 × 10 −23 J/K), and h is the Planck's constant (6.626 × 10 −34 Js), C is the velocity of light (2.998 × 10 8 m/s).Accurate land surface emissivity (LSE) is essential in retrieving LST and various approaches have been utilized to acquire LSE [43,44].An easy-to-apply method, called NDVI Thresholds Method (NDVI THM ), was applied to obtain emissivity values.This method is a modified semi-empirical approach and more detailed description is available at [45].

Retrieval of Surface Temperature from TASI Images
Temperature Emissivity Separation algorithm (TES) reported by [46] was employed to obtain surface temperature for TASI images.This method is accomplished by three modules: Normalized Emissivity Method (NEM), Ratio (RATIO), and Maximum-Minimum Difference (MMD).In TES, NEM is used to estimate spectral contrast for all TASI bands with an initial LSE (typically 0.98).The NEM derived emissivity is rationed to relative spectral in the RATIO module and variations of emissivity spectral MMD can be characterized through a logarithmic function of the minimum emissivity value εmin.Yang, et al. [30] established a semi-empirical relationship among εmin and MMD for TASI data, and the result of statistical regression was expressed as:  [15], five distinctive surface types (i.e., vegetation, bare soil, water bodies, built-up areas, and large rooftops) were identified.Large rooftops were regarded as a single land-use type since a diversity of manufactory, markets, factories, and planned residential areas were covered with huge metal material.The supervised signature extraction with the support vector machine (SVM) algorithm was employed to extract main landscapes.Training the SVM with a Gaussian radial basis function requires two critical parameters: C is a regularization parameter that controls the trade-off between maximizing the margin and minimizing the training error, while γ describes the kernel width [47].These two input parameters were tuned using a 10-fold cross validation method with the help of library LIBSVM [48], and the C and γ was set to 300 and 3. Accuracy assessment of classification is important; thus, the statistical measures involving overall accuracy (OA) and Kappa statistic (KC) were calculated for different image data using a variety of test samples collected from WV2 and Google Earth TM , as shown in Table 3.In order to further examine the spatial patterns of the thermal emission in the urban at the microscale level, a high-resolution classification map was produced using the WV2 image and the SASI image.In this context, a method on the basis of SVM classifier in combination with an object-oriented approach was used to identify different urban categories (Table 4).A watershed approach integrating the texture gradient with a marker selection was implemented to perform image segmentation, the result of which was fused with the SVM classification map [49,50].The input parameters required for running SVM were fixed using a 10-fold cross validation method through library LIBSVM, and the slack trade-off C and bandwidth γ for radial basis kernel function was set as 8 and 1.More specifically, the SASI image acquired on 15 August 2010 was used to finely discriminate concrete from mixed asphalt based on distinctive spectrum features in infrared.Accuracy assessment was conducted using the test samples collected from field survey, WV2 and Google Earth TM for the classification result, and the OA and KC was 88.37% and 0.902, respectively.Additionally, the classification result was further checked visually and refined manually, in order to have the most reliable product for the next step.The results of LST map and associated descriptive statistics are shown in Figure 4 and Table 5, respectively.In 2006, the maximum mean LST was recorded for metal rooftop, followed by built-up areas and bare soil.Compared with other land cover features, the vegetation had a lower average LST of 23.5 °C.For the other three years, similarly to the results of 2006, the vegetation all exhibited a lowest average LST, while impervious surface and bare soil both had an average LST over 33 °C.This is attributed to the thermal capacities of vegetation which allows it to easily release heat stored through canopy transpiration especially in hot weather.Especially, influenced by microclimate factors and the surrounding ecosystem, bare soil distributed in the North China Plain usually has low thermal inertia due to the low moisture content, which means that less energy is stored effectively and it is more readily dissipated as sensible heat, resulting in a similar thermal response behavior as that of the urban impervious surface.It should be noted the temperatures of water were higher than vegetation for all four years.This phenomenon can be attributed to various reasons.First and foremost, the type of cropland in our study area was irrigable land, which emits less thermal energy than most other farmed areas.Second, the small area of water bodies and wetlands, especially in sub-urban regions of Shijiazhuang usually surrounded by large area of bare soils, could easily be identified as mixed pixels and classified as water at the 30 m resolution of TM.Third, the vegetation covers in our study area included a considerable fraction of forest, which has much lower LST than water due to the fact that the dense vegetation reduced the amount of heat stored through transpiration.Both of these factors mean the cropland tends to behave slightly like water bodies.It should be mentioned that [51,52] reported this phenomenon in another city (Beijing), which is located also in North China and close to Shijiazhuang.In addition, a large area of rooftop surfaces was converted from fallow land or vegetation areas and showed the maximum average LST and standard deviation (SD).This clearly indicated the process of urbanization, in which much of the land cover changed from pervious to impervious surfaces, giving rise to an increase in LST.Interestingly, the SD value was relatively large for water compared to vegetation, implying that the surface temperature of water actually had a higher variation.This may be caused by the fact that water bodies in the study area included diverse compositions such as rivers, lakes, reservoirs and ponds, the LST value of which all varied greatly.The LST statistical description of each imperviousness category is shown in Table 6.It was observed that the high-density urban area (percent ISA categories > 75% ISA) had the highest mean LST values in all the years.The mean LST values of relatively developed areas (percent ISA > 30%) were about 2-4 °C higher than those of rural areas (percent ISA < 30%).Comparison between the four cases reveals that the larger the percentage of ISA, the higher the mean temperature in urban areas.SUHI intensity was used to further interpret the retrieved LST, since it can minimize the difference in meteorological variation and atmospheric conditions resulting from different observing times of the Landsat overpasses [8].SUHI intensity is a useful indicator for assessing urban heat islands, and is computed as the mean temperature difference between the area considered as urban and that excluding urban zone.The estimated result of SUHI intensity is computed according to the administrative zoning map of Shijiazhuang provided by Shijiazhuang land resources bureau, and listed in Table 7.It was observed that the SUHI intensities in summer were about 2-4 °C in Shijiazhuang at approximately 10:50 am, revealing that the average surface temperature of the urban area is generally higher than that of the surrounding suburban areas for the summer season.This phenomenon is more obvious in many other megacities in China, and previous studies have shown that higher SUHI intensity was found in Beijing, Shanghai, and Shenzhen [12,15,53].

Relationships among LST and RSIs, ISA%
A scatterplot was used to evaluate the accuracy of percentage of impervious surface area retrieved by MESMA.As illustrated in Figure 5a, the regression analysis exhibits a relationship with an R 2 of 0.80 and a slope of 0.70.Our study showed an overall RMSE of 8.80%.RMSE of 9.49% and 7.83% were achieved for less-developed areas (with ISA less than 30) and the developed area (with ISA greater than or equal to 30), respectively.Moreover, compared with the results of some prior studies based on the spectral mixture analysis, such as [54] (R 2 = 0.70, 200 samples), [55] (R 2 = 0.73, 50 samples) and [56] (R 2 = 0.78, 64 samples), the result reported by our study was acceptable and reliable.A total of 500 randomly generated sample points obtained from the 2010 imagery were employed to investigate the relationship between LST and ISA%. Figure 5b reveals a relatively weak positive correlation between LST and ISA%, with the determination coefficient of only 0.42.Several studies including those conducted by [57,58] reported that ISA% is positively correlated with the surface temperature.However, it should be mentioned that [24,27] reported that the coefficients of determination of LST to ISA% were fairly low (R 2 was less than 0.4).In our study, the relationship demonstrated is not obvious, indicating there might be more complicated associations between them.This probably can be partly attributed to that the LST of the impervious materials is complex and susceptible to the arrangement and configuration of the buildings and paved road [59].Nonetheless, most of the higher LST values (average, maximum and minimum) were associated with higher percentage ISA for the selected samples, asserting the role of impervious surfaces in shaping the UHI effect.To further explore this elusive contribution, a zonal analysis was utilized to evaluate the average LST at each 1% increment of ISA%.A good linear relationship between the average LST and ISA% statistically determines by the analysis and shown in Figure 5c, confirming ISA% can describe the variations in LST.
In this study, four RSIs including NDVI, VIUPD, NDBI, and BCI were employed to establish relationships with LST. Figure 6 shows the distribution of the NDVI, VIUPD, NDBI, and BCI indices from the Landsat TM image on 15 August 2010.In addition to a sample with 500 randomly generated points available for scatter plot analysis, the corresponding average temperature values at 0.01 intervals were also calculated for each index.The coefficients of determination R 2 between average LST and various remote sensing indexes were computed.Figure 7a,b shows that NDVI was roughly correlated to LST.For the region where density of vegetation was lower, this correlation was weak with some slight fluctuations.However, when the vegetation coverage was above a certain level (e.g., NDVI ≥ 0.10), an obvious decreasing trend in temperature was observed.Concerning VIUPD (shown in Figure 7c,d), when the density of vegetation cover was lower, the correlation between LST and VIUPD was positive with the coefficients of determination (R 2 ) greater than 0.98.Conversely, for the region where vegetation density was high, the correlation between LST and VIUPD was negative with R 2 around 0.99.The value 0.52 for VIUPD could be considered a threshold to split dominated land cover.When VIUPD was larger than 0.52, the dominated land type was vegetation covers and urban impervious surfaces; while when VIUPD was less than 0.52, the dominated land type is understory vegetation and the mixture of vegetation covers, bare soil and a small account of water pools.It should be noted that since the range of 0.4 and 0.5 was dominated by understory vegetation and the mixture of vegetation, bare soil and a small account of water bodies, it is reasonable that LST increased in this range.The regression analysis between VIUPD, NDVI and LST presents the relationships quantitatively and clearly in Equations ( 13) and (14).Results show that LST is sensitive to vegetation covers in terms of a negative correlation; yet, the retrieved surface temperature is more sensitive to the changes of VIUPD, especially for the sparse vegetation covers since the impact of soil background and understory vegetation is taken into consideration and eliminated by VIUPD.
NDBI was used to quantitatively analyze the correlation between built-up area and LST.The relationship among them reveals a statistical correlation (Figure 8a,b).The lower NDBI corresponded to a low LST, whereas the higher NDBI coverage was associated with high-density built-up areas that lead to higher LST.However, researchers previously found that bare soil and impervious surfaces experienced a wider variation in surface temperature than vegetation cover [60,61], and our analysis above is also consistent with previous findings, and demonstrated that bare soil has a similar thermal response to urban impervious surfaces.To remove the influence of bare soil in the low-density built-up class, a quantitative analysis was carried out using the BCI, which is effective in separating impervious surfaces and bare soil.Figure 8c,d shows that positive correlations among the LST and BCI was more significant than that of NDBI and LST and larger variations response to surface temperature was exhibited in Figure 8d.For the region occupied by vegetation and bare soil (NDBI < −0.3 and BCI < 0.0), there is a relatively large range in remote sensing values for BCI compared to NDBI to describe LST.The regression analysis between NDBI, BCI and LST presents the relationships quantitatively and clearly in Equations ( 15) and ( 16).

Analysis of Urban Thermal Patterns at the Microscale Level
An evaluation of surface temperature retrieval is important for SUHI applications.In this study, 15 field measurements, collected almost simultaneously with the TASI overpasses, were used to evaluate the accuracy of temperature and emissivity derived from TASI.More specifically, four predominant land covers including construction material (five samples), soil (two samples), vegetation (four samples) and water body (two samples) collected on August 15, 2010 were used to assess the performance of TASI.Results [33] showed a satisfactory agreement with ground data in the same program of Shijiazhuang, with root mean square error value of 0.046 for broad-emissivity and 2.2 °C for temperature.This suggests the airborne based LST can be acceptable and reliable for further analysis.
The finer LCTs map obtained from high spatial resolution data is shown in Figure 9. Surface temperature of much finer spatial resolution retrieved from TASI data were shown in Figure 10.The microscale study in this research mainly focused on profoundly examining the thermal response of urban land covers using multi-temporal TASI imagery.Mean temperature and SD of different land cover types was calculated for the TASI image of different time phases, and results are shown in Table 8.In addition, the frequency distributions of LST over various LCTs (Figure 11) were utilized to comprehensively investigate the thermal behavior of diversified urban materials.Concerning the time phase of early morning, the temperature differences among various urban covers were smaller compared to other periods of time (noon and evening).Most land cover had a relatively similar thermal response at this time, and their surface temperature differences were minimal.Water was observed to show the highest temperature with an average value of 23.24 °C.For the noon and evening, as can be clearly seen, the most important surface components contributing heat to urban areas were rooftops, asphalt, and concrete which have lower thermal inertia.Contrary to the morning, water and vegetation area had lower temperatures during the noon.When it came to the evening, the temperature of impervious surface and bare soil was still higher than other cover types due to significantly different thermal bulk properties.For the two most commonly used impervious materials (concrete and mixed asphalt), both of them had a higher variability and wider extent in temperature likely due to the heterogeneity of imperviousness surfaces.This heterogeneous characteristic is determined by many factors such as thermal properties and anthropogenic activity [62,63].Interestingly, concrete was warmer than asphalt during the noon and evening, but appeared cooler in the morning.This can be attributed to the thermal conductivity and heat capacity of concrete which allow it to absorb more radiation during the day.In addition, a description of surface temperature reveals that metal rooftops comprise both the highest maximum and the lowest minimum temperature; with a maximum mean LST of 43.99 °C and a maximum SD of LST of 5.09 °C.Contrarily, vegetated areas were cooler by about 3 °C than impervious areas during the noon by average.Moreover, a surface temperature difference for regular vegetation (mainly grassland and cropland) and tree canopy showed that the former was slightly cooler than the latter, with a temperature difference of about 0.5-1 °C not considering the systemic errors from TES algorithm.Slight temperature deviation between the regular vegetation and trees is probably due to many factors, such as the water content of the vegetation and the mixed energy emission from the ground beneath and vegetation covers.The relatively high water content in the grassland and crop is perhaps the major reason for this phenomenon in our study.We found that a considerable fraction of (a) ( the short grassland was interspersed with urban wetland (park and pools), which would contribute to the lower LST.In addition, cropland in this region basically was irrigable land, thus making the LST of cropland tending to be slightly lower than the trees in urban areas.
In order to carefully discover the surface temperature responses in the local domain, the two subset images covering the same region acquired on the noon and morning of 7 August 2010 were selected and subtracted.Experiments found that the maximum value of averaged temperature difference was recorded for rooftop, followed by concrete and asphalt, with a respective value of 13.82 °C, 11.52 °C, and 11.11 °C.As expected, water body has the minima of 5.58 °C.Concerning the percentage pixels occupied by temperature differences with respect to three impervious surface types (rooftop, concrete and asphalt), the STD values was 3.96 K, 2.73 K and 4.92 K.Moreover, it was observed that about 80% of the pixels fall under the temperature difference level of 12.6 K, 10.3 K and 9.8 K, respectively, for these three impervious surface types.This reveals that the surface temperature differences vary significantly in the microscale domain and the impervious materials are more responsible for urban thermal variations.Chudnovsky, et al. [64] emphasized that the land cover variability, which is due to the complexity of the urban composition, manifests itself in the development of the UHI phenomenon.Furthermore, since the temperature difference of the rooftop, concrete and asphalt is obvious, our experiments suggest that this discrepancy information can be treated as auxiliary data to classify fine urban material.References such as [65][66][67] even used thermal properties of surface temperature or surface emissivity to classify urban regions.Complementing these studies, we believe that the temperature difference obtained from high-resolution sensors can facilitate in discriminating urban land cover.
Based on the urban thermal spatial patterns, two profile transects along the north-south direction (in Figure 1) across different TASI images were constructed to examine the distribution of UHI in the city of Shijiazhuang (Figure 12).The two profiles were drawn in a degraded spatial resolution of 30 m to include most of the urban cover types in the study area.For the first temperature profile, numerous peaks and valleys occurred on the right-hand side, especially for noon.This is associated with the fact that built-up areas were mostly localized in the southern part (right-hand side) of the study area and the suburban areas were mainly in the northern (left-hand side).Results clearly show the larger temperature fluctuations existed in the built-up land such as commercial and industrial areas, whereas the lower fluctuations were observed in the left-hand area, which was mainly covered by cropland and low dense built-up.Similarly, the highest peaks were recorded in the central urban area for the second temperature profile.Furthermore, it was noted that as the two profiles moved from suburb areas to urban covers, the increase of surface temperatures continued as they entered into more impervious surfaces.This demonstrates that the urban area is warmer than the surrounding suburb areas, and thereby further confirming the variation of urban thermal patterns existed in the study area.In addition, the thermal variation in the early morning was much smaller than that during noon and the evening.This heat variation was found not only among the city center and surrounding areas but also between highly-dense built-up and low-density built-up areas.

Comparison of the TASI and LANDSAT TM LST
In our study, the airborne imagery was acquired at a high-spatial resolution and at a local time different from the overpass time of the Landsat.The spatial and temporal disparity for the two sensors makes it hard to make a comparison between TASI LST with LANDSAT LST.To provide an assessment of the retrieved LST of two used sensors, we carried out the comparisons of LST from the sensors in the urban-suburb region when both sensors had observations available on 15 August 2010.The overpass time was about 13:30 (Local time) for TASI and about 10:50 (Local time) for Landsat TM on that day.Before the comparison, a simple spatial resampling was applied to both the Landsat and TASI pixels for resampling them to 30 m which is the resolution the visible bands of TM.
A scatter plot of the TASI and Landsat LST is shown in Figure 13a.The determination correlation (R 2 ) is 0.62 with a bias of 3.71 K, indicating that TASI observed a higher LST compared to Landsat TM for most of the pixels.The Root Mean Square Error (RMSE) is 5.31 K.Most of the relatively larger differences with a value above 5 K are found in the urban region (south in Figure 10).In addition, the pixel scale comparisons were categorized and examined according to the four dominant land-cover types, including impervious surfaces, vegetation, waterbody and bare soil.For the waterbody, LST from TASI has a relatively smaller bias, which is 1.47 K relative to the Landsat LST and the corresponding RMSE is 2.36 K.The vegetation has a positive bias of 1.84 K and an RMSE of 2.79 K between the TASI and TM LST.Most of the pixels with a large deviation between the TASI LST and the Landsat LST are found in the impervious surface areas and bare soil types, and the bias and RMSE is 4.03 K and 5.62 K for impervious surfaces, and 3.71 K and 5.31 K for bare soil.This implies that a larger surface temperature rise occurred on the impervious surfaces from the overpass time of Landsat (about 11:50 in local time) to the overpass time of TASI (about 13:30 in local time) on 15 August 2010.
To understand the spatial variability, the LST profile of TASI LST and LANDSAT LST was examined on the two line transects in Figure 1 (see in Figure 13b).TASI observed a higher LST than Landsat for most of the pixels on the line mainly due to the time lag of about two hours.At the same time, TASI showed a much higher variability of LST especially for the intensive built up area.Figure 13b confirms that more detailed information of the surface temperature can be captured by the TASI sensor.The LST from the TASI sensor was assessed in a prior study [33].An accuracy of 2.

Discussion
Through the inspection of thermal values of LANDSAT TM (Table 5) and TASI (Table 8), we found a basically similar thermal response.A decrease of thermal values was recorded with impervious and soils, and followed by vegetation and water body.Our results at both levels indicate the composition of urban cover significantly affects the magnitude of LST.The land cover feature that contributes most to the magnitude of LST is the impervious surface.Within the study area, the mean LST of impervious surfaces was about 2-5 °C higher than that of the vegetation covers.It is particularly worth noting that the integration of high-precision community classification schemes and fine thermal values is effective in characterizing thermal pattern variability in urban areas.The TASI retrieved LST over each land cover category was basically similar to that of Landsat images, further supporting the finding that diversity of impervious surfaces (rooftops, concrete, and mixed asphalt) contributes most to the urban heat island.More specifically, the high percentage of impervious surfaces in the CBD creates a high value in LST at noon, revealing the influence of various dark impervious surfaces (i.e., concrete and asphalt) on energy irradiance.The thermal response in residential areas, which contains a blend of impervious surfaces and vegetated covers, is a function of the density of the houses, the amount of vegetation and the age of the neighborhood.
Through the inspections of the correlations of RSIs, ISA% and LST, we found that the increase of vegetation cover would decrease LST; whereas the increase of impervious surfaces (buildings and paved roads) would increase LST as well as facilitating the urban heat island phenomena.This is mainly due to the fact that changes in percentage of vegetation covers and impervious surface area would lead to modifications of the land surface characteristics that greatly affect urban heat energy budget and balance, such as albedo and evapotranspiration [68].On the other hand, although urban thermal environments can be depicted with an urban ecology landscape metric in terms of the ISA% using remote sensing, the effect of impervious surfaces on LST is complicated, and revealing the linkage between them is not straightforward.Since urban LST varies with the combined factors of building arrangements and shape such as largest patch index and patch density of paved surfaces [59,69], it is surprising to see that LST cannot be explained by ISA% with a high determination coefficient in Figure 5.Despite all this, the analysis of relationships between LST and percentage ISA reveals that surface temperature variations could be better accounted for by percentage ISA than by the commonly used NDVI based on the zonal analysis.This research replenished previous achievements, which suggested that percent imperviousness is suitable for LST studies.According to this context, our study suggests that ISA% retrieved from remote sensing provides an additional valuable metric in the study of UHI effects and attention should be focused on it further research work.
Given that RSIs are commonly used in the assessment of land type coverage, relationships between four RSIs and LST were established by using regression analysis in the present research.Results show that with different densities of vegetation covers and impervious surfaces, there are different patterns of relationships between LST and RSIs.When the density of fallow land and impervious surface was predominant, and in the correlations between LST and NDVI, NDBI remained unstable because these two indices are intended to highlight only one land cover feature, yet confusions among other land cover features (particular impervious surfaces and bare soil) could not be effectively addressed.In this study, two physically based remote sensing indices involving VIUPD and BCI were proposed to account for the variation of surface temperature.According to the conceptual V-I-S framework, urban land cover types (excluding water) can be regarded as a combination of several basic biophysical components.On this basis, VIUPD and BCI were proposed to quantify biophysical compositions for an urban area.These two indices both well consider the role of lower density land covers, and therefore have greater potential in capturing the LST trend.A full validation of the VIUPD and BCI was not conducted in our study; however, the results from the comparison of them with commonly used index were promising.The comparisons with NDVI and NDBI demonstrated the strength of these two indexes in the context of urban thermal modelling.However admittedly, RSIs obtained from medium resolution images is not good enough to describe urban thermal environmental due to its limitations in capturing the fine surface thermal characteristics.Therefore, the RSIs' approach needs further study to determine whether it a reliable tool to analyze SUHI effect.
Regarding the two temperature profiles obtained from TASI images (shown in Figure 12), both profiles exhibited similar properties during the three time phases.Surface temperatures increased significantly with the increase of solar radiation.These temperature profiles exhibited higher variability, and larger temperature fluctuations in the built-up land such us commercial and industrial areas, whereas the lower fluctuations were observed in cropland and low dense built-up.It should be noted that these heat variations were found not only among the city center and surrounding areas, but were also seen between dense built-up and low density built-up areas with greater fluctuation.This indicates that the distribution and arrangement of urban land covers plays a significant role in the SUHI effect, and moreover, it reveals that the difference in building density may lead to different thermal patterns.It is closely related to the fact that the modification of the variability or the complexity of buildings can change solar heat absorbed by the impervious surfaces because of the varied exposure surfaces.
By explicitly investigating the spatial pattern of urban thermal emissions on two different spatial scales, the approach used in this study expands and strengthens our understanding of the roles and interactions of LCTs on LST in urban landscapes.Study suggests that careful consideration is required for land use planning, particularly in urban-suburb areas where a widespread change in LCTs can significantly impact thermal responses.Regrettably, based on the above context, it was found that, over that five year period since 2006, a significant increase of built-up area and a decrease in vegetation covers (cropland, forest, and grassland) had taken place in the city of Shijiazhuang resulting from urban sprawl.The built-up area increased from 638.815 km 2 in 2006 to 681.259 km 2 in 2010, while vegetation decreased from 1660.359 km 2 to 1629.085 km 2 in the same periods.In terms of the areal extent of the diversity categories based on percentage ISA, overall the areal extent of percentage ISA > 30% was 616.252 km 2 in 2006, and it increased to approximately 688.293 km 2 in 2010.A sharp increase also occurred in the categories of >75% ISA.These LCTs changes inevitably lead to a higher frequency of UHI effect [28,52].Just as the current study being reported, the Landsat thermal images have revealed obvious heat island effect in the Shijiazhuang city during summer days since 2006, with SUHI intensities of 2-4 °C.Therefore, urban planners attempting to mitigate the negative impact of urban sprawl on urban heat island should give serious consideration to both the area percentage of various land-cover types and their spatial distribution as well.

Conclusions
In this paper, the surface thermal characteristics of Shijiazhuang, China, was studied at two different spatial scales.Multi-scale remote sensing images mainly from Landsat TM and TASI were applied to analyze the spatial variations of LCTs as well as the SUHI effect.To assess the influences of land cover types on LST, four scenes of Landsat images acquired in summer days were used to retrieve LCTs, ISA% and LST, as well as quantifying the relationships among them.Meanwhile, detailed urban thermal characteristics of the central Shijiazhuang city were examined using multi-temporal airborne thermal images.
Results show that satellite data (LANDSAT TM) are very useful for investigation of the SUHI effect at a mesoscale level; however, they do not lend themselves to gaining a better understanding of which specific surfaces contribute to the urban heat island effect and to what degree.Nevertheless, high resolution thermal image data can offer detailed characterization for the complex urban landscapes, and hence achieve a more realistic and through assessment of the urban heat island.
Overall, several conclusions were drawn from our research: (1) SUHI existed in the Shijiazhuang in summer days, and apparent discrepancies in LST were found among urban land cover types with a difference at noon as high as 7 °C for satellite images and as high as 10 °C for airborne images; (2) quantitative analysis between LST and various indices revealed that a large variation of LST even existed for the same land cover type (i.e., impervious surface, vegetation); however, these LST variation can be described by the commonly used indexes (with R 2 of 0.74 and 0.65 for NDVI and NDBI, respectively).The other two indexes VIUPD and BCI which take into account biophysical

Figure 1 .
Figure 1.Study area in Shijiazhuang, China.(a) WordView-2 RGB color composite of the study area.Line 1 (Left) and 2 (Right) are two north-south direction profile transects.(b) Location of the study area and Landsat TM imagery.

4. 1 .Figure 3 .
Figure 3a,b show the two examples of the classified map in 2006 and 2010.The percentage of impervious surface areas retrieved from the improved spectral mixture analysis in 2006 and 2010 is mapped and shown in Figure 3c,d.

Figure 9 .
Figure 9. LCTs map derived from WV2 image and SASI image using SVM and watershed segmentation method.

Figure 10 .
Figure 10.LST derived from TASI imagery on (a) the morning of 7 August 2010; (b) the noon of 7 August 2010; (c) the noon of 15 August 2010; (d) the evening of 25 July 2010; and (e) the evening of 27 July 2010.(Note that the white regions are those could not pass the data quality control.)

Figure 11 .
Figure 11.Histogram of LST over all land covers on (a) the morning of 7 August 2010; (b) the noon of 7 August 2010; (c) the noon of 15 August 2010; (d) the evening of 25 July 2010; (e) the evening of 27 July 2010.

Figure 13 .
Figure 13.Comparison of the LST from TASI and Landsat TM on 15 August 2010 (a) Scatter plot of the LST; (b) the LST profile on two line transects (Lines marked in Figure 1).

Table 2 .
Characteristics of the data used in this study (the overpass time is in local time).
Retrieval of Land Cover Patterns from Landsat and TASI Images Image classification was performed to map LCTs patterns for Landsat TM images acquired on 5 September 2006, 23 August 2007, 12 August 2009 and 15 August 2010.Referring to the previous researches of LCTs classification regarding to UHI

Table 3 .
Accuracy assessment of classification of Landsat TM images of 2006-2010.

Table 5 .
Descriptive statistics of LST (in °C ) derived from Landsat TM images over different land cover types.

Table 6 .
Descriptive statistics of LST (in °C) derived from Landsat TM images over categories of percent impervious surface area.

Table 8 .
Descriptive statistics of LST (in °C ) derived from TASI LST images over different land cover types in 2010 (the overpass time is in local time).