As the basic physical quantity of energy balance in the geogas system, land Surface Temperature (LST) has considered to be one of the most important indicators for assessing of the urban thermal environment. The urban thermal environment obviously affect by seasonal variation. Long time series temperature data is often used to analyze daily, monthly, seasonal, and annual changes in urban heat islands (HUI) . Some indicators such as atmospheric circulation, cloudiness and aerosol are used to analyze the seasonal trends of UHI . Compared with atmospheric temperature, spatial heterogeneity of LST existed strongly on urban scale  . LST is more strongly influenced by Land-Use and Land-Cover Change (LULC) and human activities  . What’s more, the impact of LULC on UHI is affected by seasonal variation  . Clarification as to how changes in LULC contribute to observe variations in LST by seasonal changeable is necessary in order to reduce UHI impact and enable urban areas to adapt to climate change effectively.
Due to extreme urbanization, more and more nonurban land cover in urban space have been replaced by Impervious Surface Area (ISA)  . The large impervious surface in urban will result in significant impacts in cities, including reduced evapotranspiration, increased surface runoff, increased storage and transfer of sensible heat, reduced air and water quality  . However, the effects of urban development density on LST are not well understood. Therefore, this study focuses on the impact of urban development density within UHI in different region. According to urban overall planning of Fuzhou, the study area was divided into three region, internal city region, intercity region, and rural-urban continuum.
In this paper, Fuzhou city was selected as the study area. Because 30 meter or lower resolution data (e.g. Landsat-class) are more suitable to map the extent and change of urban land cover   , fully constrained least squares linear spectral mixture analysis (LSMA) was used to obtain the sub-pixel LULC coverage information and the high-resolution image are used to verify the extracted results by area comparison method. According to the density of urban development, the study area is divided into three regions. The percent of ISA in each region based on pixels was divided to generate several of ISA coverage categories to analyze the relationship between LST and land cover patterns. Seasonal UHI plaque area and LST boxplots were used to analyze the impact of each region on UHI. Base on the vegetation-impervious-soil (VIS) model, ternary triangular chart method was applied for LST for each season to analyze the impact of each season on UHI .
2. Study Area and Data
2.1. The Study Region
Since 2000, the population of downtown of Fuzhou has grown rapidly, from about 1.485 million in 2000 to about 3.5 million in 2016  . With the expansion of the city and the rapid increase of the impervious surface, local high LST areas appear in the city, and the UHI effect is serious. Therefore, choosing Fuzhou as a research area has typical significance. According to the “Urban overall planning of Fuzhou (2011-2020)” , minimum enclosing rectangle in the central city was selected as the research area. Combined with urban space structure planning in Fuzhou , the research area was divided into three region, internal city region, intercity region, and rural-urban continuum in Figure 1.
2.2. Data and Preprocessing
In order to reflect the seasonal characteristics of UHI, the mid-season date of Landsat-8 images were selected to represent the spring, summer, autumn and winter. Because Fuzhou belongs to the subtropical marine monsoon climate, it’s difficult to obtain seasonal images in the same year in Cloud-prone and Raining Area. Therefore, bi-temporal Landsat-8 images (path 89/row 84) acquired on April 17, 2014, July 27, 2016, October 23, 2013 and December 13, 2014. The Google Earth image (acquired on October 1, 2013, April 11, 2014, December 18, 2014 and July 27, 2016) with 1 m spatial resolution were used as ancillary data for the assessment of the accuracy of sub-pixel LULC coverage information.
All data used in this study were projected to the Universal Transverse Mercator (UTM) based on the high-resolution images of the same phase. The RMSE of the projection was <0.5 pixels. For Landsat images, Second Simulation of the Satellite Signal in the Solar Spectrum (6S) and MODTRAN tools were used for radiometric calibration and atmospheric correction. In order to eliminate the influence of waters in the study area on the spectral analysis of low reflectivity features, modified normalized difference water index (MNDWI) was calculated to extract water .
Figure 1. The Landsat 8 OLI images of the study area (red, band 5; green, band 4; blue, band 3). The internal city region, intercity region, and rural-urban continuum was represented as Site 1, Site 2 and Site 3.
The methodology for this study consists of four steps. Step 1 includes retrieval of LST and fractional covers retrieval using LSMA. Step 2 includes accuracy assessment of sub-pixel LULC covers by area comparison method. Step 3 divides the fractional values of ISA base on the density of urban development and analyzes the differences in the effects of each LULC on the LST of the four seasons with ternary triangular chart method and contribution rate. Step 4 according to the urban overall planning, divides the study area into three regions and analyze the seasonal spatial distribution of thermal environment in each annular region.
3.1. Spectral Unmixing and Accuracy Assessment
In this paper, the spectrum of pure pixels in the study area is unknown. Therefore, minimum noise fraction rotation (MNF) and pixel purity index (PPI) used to extracting the spectrum of pure pixels in the study area and regards and (fi the fraction of end-member i within the pixel)as constraint condition. The fully constrained least squares of LSMA used to extract fractional LULC values from OLI imagery (Formula 1). The advantage over unconstrained methods is that: 1) data segmentation and root mean square (RMS)calculated for each loop computing; 2) not calculate mask pixels with all band values of 0 .
where is spectral reflectance in band b; N is the number of end-members; fi the fraction of end member i within the pixel; the reflectance of end member b in band i; is the residual error for band i.
Accuracy assessment of sub-pixel fractional covers by site investigation is time-consuming and laborious. Using the high-resolution image to extract fractional covers information, not only can obtain the higher precision than the low-resolution image extraction results, but also can be used as the validation of accuracy of the lower-resolution image .Six test sites within the study area about 2 km2 were selected for the accuracy estimation of fractional covers. The object-oriented classification method and the artificial visual interpretation were used to extract the LULC information of the preprocessing high-resolution images (the acquisition time is close to the Landsat image).Calculate the area of the Lands at image using (Formula (2)) and Mean Error (ME) between the Landsat and Google Earth images(Formula (3)).
where is Landsat images’ fractional covers area; is Landsat image’s spatial resolution; n is test site’s total pixel numbers.
where is Google Earth images’ fractional covers area.
3.2. Landsat Surface Emissivity Calculation and LST Retrieval
The Landsat 8 TIRS includes two thermal infrared bands. Since the band 10 is located in the lower atmospheric absorption are aand band 11 has larger uncertainty data as USGS said, so the band 10 is more suitable for single-band’s LST Retrieval calculation  . Convert the digital numbers of thermal bands to top of atmosphere (TOA) radiance using:
where is the TOA radiance image of the thermal band; ML is the band-specific multiplicative rescaling factor and AL is the band-specific additive rescaling factor; QCAL is the pixel digital number for thermal band 10.
There are many LST retrieval algorithms for Landsat images. The mono-window algorithm improved by Jimenezmunoz and Sobrino which requires less atmospheric correction parameters and has higher precision is widely used in Landsat TIRS LST retrieval.  (Formula (5)).
where ε is surface emissivity; γ and δ are the estimated parameters from the linear formula of Planck’s law; , and are three functions gained from actual atmospheric parameters (Formulas (6)-(8)).
where is atmospheric transmission; and are up welling and down welling radiance.
3.3. Contributions of the Landscape Pattern to the Urban Thermal Environment
Due to the range of LST is inconsistent in different season, it is impossible to directly compare the inversion results of the four seasons’ LST. In this paper, contribution index (CI) was used to analyze the difference in the average LST between each fractional covers area in different urban development region and entire study area (Formula (9)). CI estimates the contribution of landscape pattern in each urban development region to the LST dynamic increase or decrease, and it is based on spatial and temporal dimensions to analyze the contribution of each fractional category to urban heat island effect.
where CIi is the contribution of the each fractional category to the regional LST; LSTdif・ is the average LST difference between the each fractional category region and the entire area; Si is the area of fractional category region; Ssum is the area of study area.
3.4. Normalized LST Values
As the season changes, it will affect the relative strength of the surface temperature. In this paper, LST normalized was used to unify the LST distribution range from 0 to 1 without changing the spatial distribution characteristics of the thermal environment (Formula (10)).
where is LST for a given pixel; and is the minimum and maximum LST value in an image.
4. Results and Discussion
4.1. Accuracy Analysis of the Sub-Pixel Fractional Covers
The fully constrained least squares of LSMA was used to unmix mixed pixel, and V-I-S model was used to define categories of end members . Wherein, the category with high reflectivity and low reflectivity are superimposed and summed to represent the coverage value of the ISA. Figure 2 shows the each end member fractions covers in the study.
According to the result, the spring average RMS is 0.0008, the summer average is RMS 0.0011, the autumn average RMS is 0.0016, and the winter average RMS is 0.0006, which is less than the required 0.02 accuracy requirement . To further verify the LULC extraction results, each fractional covers category was extracted from Google Earth images, which was used to assess the accuracy of the fractional covers estimate from Landsat. Six test sites were chosen in Figure 2(1a) for the accuracy analysis. ME (Formula (4)) was calculated to analyze the difference of each fractional cover extraction result from two types of images. Table 1 shows the results of the accuracy assessment of LSMA. The area of fractional covers in the six test sites showed only small differences, which explained that the number of end members in this study is appropriate, and the fractional covers extraction results are reliable.
Table 1. Statistical table of average error of land use type extraction results (%).
Figure 2. Fractional images by LSMA (a, b, c represent Fractional impervious surface, the soil and vegetation, 1, 2, 3, 4 represent spring, summer, autumn, winter, and the red sites in 1a represent test sites used for accuracy assessment).
4.2. Urban Thermal Intensity Analysis on Season Variations
4.2.1. Influence of Sub-Pixel LULC on Seasonal LST
The “Formulas (4)-(8)” were used to retrieve the LST of the four seasons in the study area. Figure 3 shows the results of seasonal LST retrieval.
An analysis with percent ISA, vegetation and soil in Figure 2 and LST in Figure 3 shows that the change of each LULC in the past three years is small. However, the spatial distribution of surface temperature in the four seasons is obviously different. According to the density of urban development, the fractional values were divided to generate several of coverage categories. Base on the coverage categories, Formula (9) was calculated to analyze the impact of LULC on LST in different season in Table 2.
Figure 3. The spatial distribution of surface temperature in the four seasons. (a) Spring; (b) Summer; (c) Autumn; (d) Winter.
Table 2. The seasonal contribution rate of land use types at all levels to the surface temperature (%).
In Table 2, the positive contribution rate indicates that the fractional covers in this interval promote the increase of urban thermal environment in the current season; the negative contribution rate indicates a reduction of urban thermal environment in promotion; the absolute value of the contribution rate indicates the influence of the coverage categories on the thermal environment in the current season.
It can be seen from Table 2, when the density of the ISA in the pixel is greater than 40%, ISA will promote the urban thermal environment, and with the increase of density, the effect is strengthened. The intensity of the four seasons is: spring > autumn > summer > winter; when the density of the vegetation is greater than 60%, it will reduce the urban thermal environment. The intensity of the four seasons is: summer > autumn > spring > winter; the effect of soil coverage on the thermal environment is greatly affected by the season, and the effects in spring and autumn are relatively close.
Since the contribution rate can only analyze the degree of influence of a single type of fractional covers on the LST. In this paper, according to the V-I-S model in Figure 4, the LULC in the pixel was defined by proportion of three types of fractional covers in mixed pixels. And an LST variable was added to the V-I-S model to analyze the impact of each fractional cover on urban thermal environment in the four seasons in Figure 5.
In Figure 5, the red map indicates that this pattern of LULC promotes the urban LST increase; the blue map represents the pattern of LULC reduce the seasonal LST. Combined with the definition of V-I-S model can be found, in the spring, high LST areas appear in industrial, bare land and desert, with no obvious low LST areas; High LST areas in summer occur in the CBD, light industrial, and covercrops, and low LST areas appear in low-density residential and forest; the high LST areas in autumn appear in the CBD, high-density residential and range land, and the low LST areas appear in cover crops and lawn areas; in winter, light industry and heavy industry are the main areas of high LST distribution, the distribution of low-density and medium-density residential areas, crops and lawns will inhibit the urban heat island effect.
Figure 4. Illustration map of V-I-S (vegetation-impervious surface -soil) model.
Figure 5. Ternary triangular charts of land use type coverage and surface temperature. (a) Spring; (b) Summer; (c) Autumn; (d) Winter.
4.2.2. The Influence of Urban Development on Seasonal Thermal Environment
According to the “Urban overall planning of Fuzhou (2011-2020)”, the study area was divided into three regions, internal city region, intercity region, and rural-urban continuum. Boxplot is a method of depicting numerical data sets by quartile graphics. In this paper, boxplot was used to analyze the distribution of seasonal urban thermal environment base three regions in Figure 6. It can be seen from Figure 6 that the LST extremes of the intercity region differ greatly, followed by the rural-urban continuum and the interior of the city is smaller. With the season change, the maximum temperature difference between the various surface areas is the largest in summer, followed by spring and autumn, and the smallest in winter. Through the analysis of the mean and median LST of each region, it can find that the distribution of intercity region’s LST is more dispersed, followed by the rural-urban continuum, and the distribution of internal city region’s LST is more concentrated. The concentration is from high to low in winter and autumn, spring, summer.
Due to seasonal variation, it is hard to directly analyze the difference of thermal environment in each season by using LST. Therefore, in this paper, formula 10was used to normalize the LST of the four seasons, and use the density segmentation technology to classify the images. According to the order from weak to strong of the UHI, The thermal environment of the study area is divided into seven levels: extremely low LST region (ETL), low LST region (LT), sub-low LST region (SLT), medium LST region (MT), sub-high LST region (SHT), high LST region (HT) and extremely high LST region (EHT). According to previous research, the three levels of sub-high LST, high LST and extremely high LST region are the main plaques that constitute UHI . To analyze the spatial variation characteristics of four seasonal thermal environments in different urban development areas, the proportion of each LST level’s plaque in three regions was calculated (Table 3).
Table 3. Area of each LST level in every developing region (%).
Figure 6. The boxplot of the four seasons surface temperature.
Table 3 shows that the proportion of the UHI plaques in four seasons was not much changed. However, in the case of extremely high LST region, the area of internal city region in summer is the highest, and the lowest in winter; the intercity region in summer is highest, the lowest in autumn; the rural-urban continuum in spring is highest, and the lowest in winter. Combined with Figure 2, it can be found that the main LULC in internal city region is ISA; Most of the intercity region are impervious, and combined with the soil and vegetation; rural-urban continuum is most covered by soil and vegetable. Therefore, Tab.3 can also be used for verification the impact of each fractional cover on urban thermal environment in the four seasons.
Seasonal variation will have a certain impact on urban thermal environment when urban LULC are similar. In this paper, the study was divided into three region, internal city region, intercity region, and rural-urban continuum, according to urban overall planning of Fuzhou. Statistics of the four-season LST’s boxplots of three regions, which was used to analyze the distribution of seasonal urban thermal environment base different development. Combined with the LST normalization results, the ratio of the seven LST levels areas in each developing region was counted to analyze the spatial variation characteristics of seasonal thermal environment in city. According to the fractional covers obtained by the LSMA, the research area was quantitatively partitioned into 10 levels. Base on the 10 regions, the contribution rate of each fractional cover seasonal LST is calculated. Combined with the ternary triangular charts, contribution rate was used to analysis of the combined effects of three types of land use on LST.
The above results suggest three major conclusions:
1) Seasonal contribution rate of each region can effectively analyze the impact of LULC on seasonal LST. It was found that when the density of the ISA in the pixel is greater than 40%, it will promote the urban thermal environment. The intensity of the four seasons is: spring > autumn > summer > winter; when the density of the vegetation is greater than 60%, it will reduce the urban thermal environment. The intensity of the four seasons is: summer > autumn > spring > winter; the effect of soil coverage on the thermal environment is greatly affected by the season, and the effects in spring and autumn are relatively close;
2) Ternary triangular charts can used to analysis of the combined effects of three types of land use on LST. It was found that in the spring, high LST areas appear in industrial, bare land; high LST areas in summer occur in the CBD, light industrial, and cover crops; the high LST areas in autumn appear in the CBD, high-density residential; in winter, the distribution of low-density and medium-density residential areas, crops and lawns will inhibit the urban heat island effect;
Based on the degree of development, the difference in the area of seasonal UHI plaques in different regions is not obvious. However, there is a difference in LST distribution: the LST extremes of the intercity region differ greatly, followed by the rural-urban continuum and the interior of the city is smaller. With the season change, the extreme LST in summer is largest, followed by spring and autumn, and the smallest in winter; the distribution of intercity region’s LST is more dispersed, followed by the rural-urban continuum, and the distribution of internal city region’s LST is more concentrated. The concentration is from high to low in winter and autumn, spring, summer.
Project supported by the Natural Science Foundation of Fujian Province, China (2018J01739) and Special Scientific Research Fund of Research Institutes Public Welfare Profession of Fujian Province, China (2019R1102).
 Yesilirmak, E. (2013). Temporal Changes of Warm-Season Pan Evaporation in a Semi-Arid Basin in Western Turkey. Stochastic Environmental Research & Risk Assessment, 27, 311-321. https://doi.org/10.1007/s00477-012-0605-x
 Wu, C.G., Lin, Y., Wang, Y., Gong, Y., Cui, H. and Xia, L. (2015) Seasonal Variations of Urban Land Surface Thermal Environment and its Relationship to Land Surface Characteristics. Journal of Harbin Institute of Technology, 47, 26-30.
 Liu, D., Yan, L.I. and Kong, F. (2013) Spatial Distribution of Land Surface Temperature in Central City Proper and the Cooling of Objects Effect: A Case Study of Nanjing. Remote Sensing for Land & Resources, 25, 117-122.
 Nsubuga, F.W.N., Botai, J.O., Olwoch, J.M., Rautenbach, C.J.D., Kalumba, A.M., Tsela, P., et al. (2015) Detecting Changes in Surface Water Area of Lake Kyoga Sub-Basin Using Remotely Sensed Imagery in a Changing Climate. Theoretical & Applied Climatology, 127, 1-11.
 Han, F., Zhang, B.P., Zhao, F., Guo, B. and Liang, T. (2018) Estimation of Mass Elevation Effect and Its Annual Variation Based on Modis and NecpData in the Tibetan Plateau. Journal of Mountain Science, 15, 1510-1519. https://doi.org/10.1007/s11629-018-4865-x
 Arnold Jr., C.L. and Gibbons, C.J. (1996) Impervious Surface Coverage: The Emergence of a Key Environmental Indicator. Journal of the American Planning Association, 62, 243-258. https://doi.org/10.1080/01944369608975688
 Wilson, J.S., Clay, M., Martin, E., Stuckey, D. and Vedder-Risch, K. (2003) Evaluating Envi-ronmental Influences of Zoning in Urban Ecosystems with Remote Sensing. Remote Sensing of Environment, 86, 303-321. https://doi.org/10.1016/S0034-4257(03)00084-1
 Pal, S. and Ziaul, S. (2016) Detection of Land Use and Land Cover Change and Land Surface Temperature in English Bazar Urban Centre. Egyptian Journal of Remote Sensing & Space Science, 20. https://doi.org/10.1016/j.ejrs.2016.11.003
 Chen, X., Sun, H., Yang, H. and Zhang, Z. (2010) Urban Spatial Growth Analysis from Satellite-Derived Imperviousness in an Oasis City. International Conference on Environmental Science and Information Application Technology, 1, 517-520.
 Frazier, A.E. and Wang, L. (2011) Characterizing Spatial Patterns of Invasive Species Using Sub-Pixel Classifications. Remote Sensing of Environment, 115, 1997-2007. https://doi.org/10.1016/j.rse.2011.04.002
 Ridd, M.K. (1995) Exploring a V-I-S (Vegetation-Impervious Surface-Soil) Model for Urban Ecosystem Analysis through Remote Sensing: Comparative Anatomy for Citiesa. International Journal of Remote Sensing, 16, 2165-2185. https://doi.org/10.1080/01431169508954549
 Myint, S.W., Brazel, A., Okin, G. and Buyantuyev, A. (2010) Combined Effects of Impervious Surface and Vegetation Cover on Air Temperature Variations in a Rapidly Expanding Desert City. Mapping Sciences & Remote Sensing, 47, 301-320.
 Jiménez-Muñoz, J.C., Sobrino, J.A., Skokovic, D., Mattar, C. and Cristóbal, J. (2014) Land Surface Temperature Retrieval Methods from Landsat-8 Thermal Infrared Sensor Data. IEEE Geoscience & Remote Sensing Letters, 11, 1840-1843. https://doi.org/10.1109/LGRS.2014.2312032
 Coll, C., Caselles, V., Valor, E. and Niclòs, R. (2012) Comparison between Different Sources of Atmospheric Profiles for Land Surface Temperature Retrieval from Single Channel Thermal Infrared Data. Remote Sensing of Environment, 117, 199-210. https://doi.org/10.1016/j.rse.2011.09.018
 Jimenezmunoz, J.C., Cristobal, J., Sobrino, J.A., Soria, G., Ninyerola, M., Pons, X., et al. (2008) Revision of the Single-Channel Algorithm for Land Surface Temperature Retrieval from Landsat Thermal-Infrared Data. IEEE Transactions on Geoscience & Remote Sensing, 47, 339-349. https://doi.org/10.1109/TGRS.2008.2007125