The Changes of Heat Contribution Index in Urban Thermal Environment: A Case Study in Fuzhou

: With the acceleration of global warming and urbanization, the problem of the thermal environment in urban areas has become increasingly prominent. In this paper, Fuzhou was selected to quantify the impact of land use cover change (LUCC) on land surface temperature (LST). The results showed that from 1993 to 2016, the land use/cover types of the study area changed greatly, especially the change of construction land, which led to an obvious change in the spatial pattern of LST. From 1993 to 2016, the spatial and temporal distribution of LST contributions in Fuzhou was uneven. The central urban area had a positive contribution to the rise of LST, while Minqing and Yongtai had a negative contribution. From the perspective of different land use/land cover types, forest or grass land, cultivated land, and water all made a negative contribution to the increase of surface temperature, while construction land made a positive contribution. Outcomes provided by the multi-distance spatial cluster analysis (Ripley’s K function) showed that there was a scale effect in the concentration and dispersion of LST; from 1993 to 2016, the concentration range of LST in the study area gradually expanded and the degree of concentration increased.


Introduction
In October 2018, the Intergovernmental Panel on Climate Change (IPCC) proposed a special report "Global Warming of 1.5 • C", which stated that in comparison to the time before industrialization, the temperature of today's world was increased significantly (0.8~1.2 • C) because of human activities, and a 1.5 • C increase in global temperature may be achieved in 2030 [1]. Rapid urbanization is taking place in the whole world, particularly in China. In January 2017, the State Council of China released the National Land Planning Outline, which predicted that China's urban area would increase to approximately 1,167,000 km 2 by 2030, compared to 89,000 km 2 in 2015 [2]. The rapid development of urbanization in China has caused serious problems, including the artificial transformation of underlying surfaces, greenhouse gas emissions, urban heat island (UHI) effects, air and water pollution, etc., which have attracted national and global concerns [3][4][5]. Typical examples of problems caused by the artificial transformation of underlying surfaces are its negative impacts on the ecological environment, the health of residents, the air quality, and the sustainability of urban economic development [6]. The deterioration of the thermal environment of cities has become a characteristic of climate change of the worldwide urbanized areas. It seriously hinders the process of urbanization and the sustainable development of the urban ecological environment [7]. Therefore, it is necessary to study the characteristics and laws of the temporal and spatial evolution of the urban thermal environment and to explore the process and mechanism of thermal environment change. These issues are very important and have a real influence on the regional ecological environment protection and sustainable development of cities [8]. Figure 1 shows the location of Fuzhou. This study was conducted in Fuzhou (118 • 22 ~119 • 57 E, 25 • 20 ~26 • 38 N). This region covers an area of 11,236 km 2 , which consists of five urban districts (Gulou, Taijiang, Cangshan, Jin'an, and Mawei) and seven county-level municipalities (Changle, Minhou, Minqing, Yongtai, Luoyuan, Lianjiang, and Fuqing). This region has a northern subtropical monsoon climate, with an average annual temperature of approximately 18~21 • C and average annual precipitation of approximately 900~2100 mm. The study area tilts from west to east, and the Minjiang River runs through the urban area and flows into the sea [21]. In 2016, the population of the study area was 6,870,600, and the GDP per capita was CNY 82,251 [22]. This study was conducted in Fuzhou (118°22′~119°57′ E, 25°20′~26°38′ N). This region covers an area of 11,236 km 2 , which consists of five urban districts (Gulou, Taijiang, Cangshan, Jin'an, and Mawei) and seven county-level municipalities (Changle, Minhou, Minqing, Yongtai, Luoyuan, Lianjiang, and Fuqing). This region has a northern subtropical monsoon climate, with an average annual temperature of approximately 18~21 °C and average annual precipitation of approximately 900~2100 mm. The study area tilts from west to east, and the Minjiang River runs through the urban area and flows into the sea [21]. In 2016, the population of the study area was 6,870,600, and the GDP per capita was CNY 82,251 [22].

Data Source and Pre-Processing
This study used cloud-free and geometrically corrected Landsat imagery from the official website (https://earthexplorer.usgs.gov/) (accessed on 25 August 2021) of the United States Geological Survey (USGS). To match the resolution of the multispectral band, USGS had resampled the resolution of the thermal infrared band to 30 m before the data product announcement (see Table A1). The pre-processing of images mainly included geometric correction, radiometric calibration, atmospheric correction, band fusion, and image cutting.

Classification of Different Land Use/Land Cover Types
Based on the standard of Land Use Status Classification (GB/T 21010-2017), 400 training samples were selected in the study area by referring to the remote sensing monitoring database of land use status in China issued by Resource and Environment Science and Data Center, Chinese Academy of Sciences (https://www.resdc.cn/) (accessed on 25 August 2021). The maximum likelihood method was used to classify the land use/land cover (LULC) types in the study area. Figure 2 shows that the LULC types include six components: water, forest or grass land, cultivated land, construction land, wetland, and bare land. By using ERDAS tools, 300 sampling points were generated in the study area by utilizing a stratified random distribution method. The overall accuracy of the classification results was greater than 85%, and the Kappa coefficient was greater than or equal to 0.75, which met the accuracy requirements (see Table A2).

Data Source and Pre-Processing
This study used cloud-free and geometrically corrected Landsat imagery from the official website (https://earthexplorer.usgs.gov/) (accessed on 22 August 2021) of the United States Geological Survey (USGS). To match the resolution of the multispectral band, USGS had resampled the resolution of the thermal infrared band to 30 m before the data product announcement (see Table A1). The pre-processing of images mainly included geometric correction, radiometric calibration, atmospheric correction, band fusion, and image cutting.

Classification of Different Land Use/Land Cover Types
Based on the standard of Land Use Status Classification (GB/T 21010-2017), 400 training samples were selected in the study area by referring to the remote sensing monitoring database of land use status in China issued by Resource and Environment Science and Data Center, Chinese Academy of Sciences (https://www.resdc.cn/) (accessed on 22 August 2021). The maximum likelihood method was used to classify the land use/land cover (LULC) types in the study area. Figure 2 shows that the LULC types include six components: water, forest or grass land, cultivated land, construction land, wetland, and bare land. By using ERDAS tools, 300 sampling points were generated in the study area by utilizing a stratified random distribution method. The overall accuracy of the classification results was greater than 85%, and the Kappa coefficient was greater than or equal to 0.75, which met the accuracy requirements (see Table A2).

Retrieval of Land Surface Temperature
The retrieval of LST of Landsat 5, Landsat 7 and Landsat 8 thermal infrared bands was carried out by using the radiative transfer equation (RTE). By making use of this method, the influence of the atmosphere on the surface thermal radiation can be estimated, then the influence of this part from the total thermal radiation observed by satellite sensors can be subtracted to obtain the surface thermal radiation intensity. Ultimately, the surface thermal radiation intensity can be converted into the corresponding surface tem-

Retrieval of Land Surface Temperature
The retrieval of LST of Landsat 5, Landsat 7 and Landsat 8 thermal infrared bands was carried out by using the radiative transfer equation (RTE). By making use of this method, Sustainability 2021, 13, 9638 4 of 18 the influence of the atmosphere on the surface thermal radiation can be estimated, then the influence of this part from the total thermal radiation observed by satellite sensors can be subtracted to obtain the surface thermal radiation intensity. Ultimately, the surface thermal radiation intensity can be converted into the corresponding surface temperature [9,[23][24][25]. The equation is given as follows: where ε is the surface emissivity; B(T s ) is the ground radiance; T s is the LST of pixels(K); L ↓ is the atmospheric downward radiance; τ is atmospheric transmittance; and L ↑ is the atmospheric upwelling radiance. The normalized vegetation index threshold method proposed by Sobrino et al. [23] was used to estimate the surface specific emissivity as follows: where ε v is the vegetation emissivity, which has a value of 0.99; ε u is the urban surface emissivity, and the value is taken as 0.92; P v is the vegetation coverage; dε includes the geometric distribution of natural surfaces and the effects of internal reflections; F is the shape factor, and the value is equal to 0.55; and the values of NDVI max and NDVI min are 0.5 and 0.2, respectively. The equation of T s was built based on Planck's law, B(T s ) can be expressed as Equation (5). Then, T s can be calculated using Equation (6) as follows: where T s is the LST of pixels (K); for TM, K 1 = 607.76 W/(m 2 ·sr·µm), K 2 = 1260.56 K; for ETM+, K 1 = 666.09 W/(m 2 ·sr·µm), K 2 = 1321.08 K.
where LST is the actual LST of pixel ( • C). Standardized treatment of LST is conducted by: (8) where T N is the standard value of surface temperature, its interval is (0, 1); LST is the surface temperature; LST max and LST min are the maximum and minimum of the surface temperatures. Figure 3 shows that T N is divided into five classes: low-temperature region, sub-low temperature region, medium-temperature region, sub-high temperature region, and high-temperature region, based on the natural discontinuous point classification method.
Twenty LST sampling points (see Figure A1) were chosen to compare the inversion LST to the collected field measured LST, respectively, and the error was within an acceptable range of about 1 (°C) [27][28][29], indicating that the inversion LST data had certain reliability.  To achieve the interannual comparability of images in different years and months, relative land surface temperature (RLST) was used to compare the contribution of the urban thermal environment in different years [26]. The formula for RLST is shown as follows: where RLST i is the relative LST of the ith year ( • C); LST i is the LST of the ith year ( • C); and LST i is the average LST of the ith year ( • C). Figure 4 shows that RLST is divided into 6 classes: RLST < 0, 0 ≤ RLS < 2, 2 ≤ RLST < 4, 4 ≤ RLST < 6, 6 ≤ RLST < 8, and RLST ≥ 8 [26].

Calculation of Contribution Index
The spatial utilization pattern of land use has large spatial heterogeneity, and the influence mechanism on the thermal environment also shows spatial difference [30]. The warming or cooling extent of per LULC type, taking into account the proportion of the total area it occupies, is quantified using the Contribution Index (CI). The CI is used to link spatial structure, as well as long-term changes in LULC, to LST intensities. The index can be calculated by using the following formula [17,30,31]: where i is ith kind of land use/cover types; LST i is the average land surface temperature Twenty LST sampling points (see Figure A1) were chosen to compare the inversion LST to the collected field measured LST, respectively, and the error was within an acceptable range of about 1 ( • C) [27][28][29], indicating that the inversion LST data had certain reliability.

Calculation of Contribution Index
The spatial utilization pattern of land use has large spatial heterogeneity, and the influence mechanism on the thermal environment also shows spatial difference [30]. The warming or cooling extent of per LULC type, taking into account the proportion of the total area it occupies, is quantified using the Contribution Index (CI). The CI is used to link spatial structure, as well as long-term changes in LULC, to LST intensities. The index can be calculated by using the following formula [17,30,31]: where i is ith kind of land use/cover types; LST i is the average land surface temperature of ith land use/cover types; LST is the average land surface temperature of the study area; S i is the area of ith land use/cover types; and S is the total area of the study area. When CI ≥ 0, the ith land use/cover type contributes positively to the LST rise; when CI < 0, the ith land use/cover type contributes negatively to LST elevation. Similarly, the CI values of the thermal environment of different regions were also calculated by taking the difference between the average LST of different regions and the average LST of the study area and the proportion of the region to the total area of the study area.

Analysis of Spatial Pattern
In this paper, the average nearest neighbor analysis (Nearest Neighbor Analysis) and multi-distance spatial clustering analysis (Ripley's K function) was used to analyze the spatial pattern distribution of the thermal environment. The method has been widely used in the study of landscape patterns and spatial distribution of the population of different species [32,33]. The nearest proximity distance is an index of the degree of proximity of point elements in geographical space, which can accurately and objectively determine the spatial distribution of point patterns [34]. The average nearest neighbor analysis is a spatial analysis method based on spatial distance judgment point spatial agglomeration. There are great differences in the pattern distribution information obtained from different spatial scales [35]. Therefore, Ripley's K function is introduced to analyze the multi-scale spatial distribution pattern of the thermal environment [36]. Ripley's K function is a kind of global statistical index, which is used to analyze the spatial distribution characteristics of point objects in the study area at any scale. According to the set radius distance, the relation function is generated statistically [37].

Average Nearest Neighbor Analysis
The geographic information system software is used to transform the RLST surface and point data, calculate the average maximum proximity ratio (Nearest Neighbor Ratio, NNR), and analyze the spatial pattern characteristics of the urban thermal environment characterized by LST. The equations are given as follows: where D O is the average observation distance between RLST; D E is the average expected distance between RLST under random spatial distribution; n is the RLST number of elements; d i is the distance between each RLST point element and its nearest point element; and A is the area of the study area. When NNR = 1, the RLST is completely random; when NNR < 1, the RLST is clustered; and when NNR >1, the RLST is dispersed.

Multi-Distance Spatial Clustering Analysis
To maintain the stability of the variance, Ripley's L(d) after the square was proposed by Besag to replace Ripley's K(d) [37,38]: where A is the area of study area; n is the RLST number of elements; and d is the size of spatial scale. When L(d) > 0, the RLST point elements are clustered on the corresponding scaled; when L(d) < 0, the RLST point elements are dispersed on the corresponding scale; when L(d) = 0, the RLST point elements are completely randomly distributed on the corresponding scale.

Characteristics of LULC Conversion
Based on Figure 2 and Tables 1 and 2, from 1993 to 2016, the LUCC in the study area changed significantly. The study area mainly consists of forest or grass land, cultivated land, and construction land. From 1993 to 2016, the area of construction land increased from 416.83 km 2 to 1647.86 km 2 , with a net increase of 1231.04 km 2 and a total change rate of 295.33%. Construction land transformed mainly from cultivated land and forest or grass land, with a conversion area of 820.01 km 2 . Most of the new construction land expanded radially along the Minjiang River and the urban transportation network. The area of cultivated land area decreased from 2400.29 km 2 in 1993 to 1539.42 km 2 in 2016, and the total change rate decreased by 35.87%. It is clearly shown that a large amount of the cultivated land has been occupied by construction land, with an area of 820.01 km 2 . The change of the area of forest or grass land showed a fluctuating trend, but the total area decreased by 297.57 km 2 in 24 years, and the total change rate decreased by 3.69%. Notably, the area of wetland has decreased from 35.85 km 2 in 1993 to 28.43 km 2 in 2016. Although the absolute amount of change in total area was not significant (7.42 km 2 ), the decreasing proportion was high (20.71%).

Temporal and Spatial Variation of Land Surface Temperature
According to Figure 3 and Table 3, it can be found that the medium-high temperature areas (medium-temperature areas, sub-high temperature areas, and high-temperature areas) are mainly distributed in the urban centers and designated towns. These areas generally have intensive buildings, concentrated industries and commerce, and intense population activities. From 1993 to 2016, along the Minjiang River, the area of the mediumhigh temperature region increased from 235.84 km 2 in 1993 to 3310.19 km 2 in 2016. The net increase was 3074.35 km 2 . It should be noticed that from 1993 to 2016, the area of the medium temperature region increased by 2681.2 km 2 , with a high growth rate of 1148%. Low and sub-low temperature areas are mainly distributed in water areas and mountainous areas with high forest coverage. The proportion of low and sub-low temperature areas in the study area showed a decreasing trend, which was 97.9%, 85.3%, 75.8%, and 70.5% for the four periods, respectively. Combined with Table 4, the changes between the different T N grades in Fuzhou can be further analyzed. From 1993 to 2016, the thermal environment of Fuzhou is mainly reflected in the reduction of the area of the low-temperature zone and sub-low-temperature zone, and the increase of the area of medium-high temperature zone (medium-temperature zone, subhigh-temperature zone, and high-temperature zone). The area of the low-temperature zone decreased significantly, from 2558.28 km 2 (1993) to 24.48 km 2 (2016), with a reduced area of 2534.32 km 2 , mainly being transformed into a sub-low temperature zone (2318.79 km 2 ). The area of the medium-temperature zone increased significantly, from 233.5 km 2 (1993) to 2914.70 km 2 (2016), with an increased area of 2681.20 km 2 , mainly being transformed from the sub-low temperature zone (2590.82 km 2 ). It can be seen that, from 1993 to 2016, the thermal environment of Fuzhou changed dramatically, mainly from a low-temperature zone to a high-temperature zone. According to Figure 4 and Table 5, the RLST of Fuzhou has changed significantly from 1993 to 2006. With the acceleration of urbanization and more intense human activities, the natural surface was gradually replaced by the artificial surface, and the impervious area increased. The area showed a trend of rapid growth with the year, from 75.73 km 2 in 1993 to a net increase of 389.70 km 2 , to 359.96 km 2 in 2016. The high-temperature area of RLST ≥8 was mainly concentrated in Gulou and Taijiang in the central urban and appeared  It is worth noticing that the direction of expansion of medium-high temperature areas and the RLST ≥8 area was approximately the same as the expansion direction of the new construction land. Therefore, the rapid expansion of urban construction land is one of the important factors which promotes the rise of LST.

Temporal and Spatial Differences in Thermal Environment Contribution
It can be observed from Figure 5 that all the CI values of Fuzhou's urban areas were higher than 0 from 1993 to 2016. This indicated that the average LST in this area was higher than the one in the study area. Consequently, Fuzhou's urban areas made a positive contribution to the thermal environment in the study area and were the main cause of LST increase (CI reached the maximum value of 0.148 in 2016). The CI values of Minqing and Yongtai were both fewer than 0, which indicated that they made a negative contribution to the thermal environment in the study area and were the main source that slowed down the increase of LST. The CI value of Luoyuan, Changle, and Lianjiang showed a fluctuating trend from a negative to a positive value, and the CI value kept increasing. The maximum CI values of Minhou and Fuqing were reached in 2000 and 2008, which were 0.044 and 0.08, respectively. study area. From 2008 to 2016, the CI value of wetland changed from negative to positive, so its contribution to the increase of LST also changed from negative to positive. During this period, the urbanization process was accelerated, and the wetland around the city were greatly developed, which reduced the area of the wetland and destroyed its functions. The following is a typical example: in 2016, the Puxiazhou Wetland was completely converted to construction lands (Fuzhou Strait International Convention and Exhibition Center was completed in 2010). Since the distribution of bare land was scattered and only accounted for a small proportion of the study area, their contributions to LST in the past 24 years were fluctuating. Figure 7 shows the correlation relationship between the change of LULC type and LST. The data in 1993 and 2016 were taken as examples to conduct a regression analysis for the area proportion of land use/cover type conversion and its corresponding LST. Among them, the LULC types with a CI value of fewer than 0 were defined as the cold sources in the study area. The cold source includes forest or grass land, water, cultivated land, and wetland. The results showed that the expansion of construction land and the loss of cold sources were positively correlated with LST from 1993 to 2016. In other words, the large-scale expansion of construction lands and the loss of forest or grass land, water, cultivated land, and wetland will lead to the rapid rise of LST.  In 2000, the area of Minhou University Town was under construction, and the urban area expanded from the central city to Minhou in an "enclave" mode. In 2013, Fuzhou New District (Mawei, Cangshan, Changle, and Fuqing) constructed a spatial structure of "one core (Mawei New City), two districts (South Wing Development Area, North Wing Development Area), and three axes (Fuping Comprehensive Development Axis, South Wing Regional Development Axis, and North Wing Regional Development Axis)". The radiation of the university town and the tilt of the policy has promoted the development of the economy and urban construction of Minhou, Fuzhou, Changle, and Fuqing. In addition, they also promoted the rise of LST in the study area.
Taking the LULC type of Fuzhou as the sampling unit and combining it with Figure 2, the contribution of different LULC types on the urban thermal environment from 1993 to 2016 was quantified and is shows in Figure 6. From Figure 6, the CI values of forest or grass land, cultivated land, and water in the four study periods were all fewer than 0, indicating that these three LULC types had a negative contribution to the increase of LST. Among them, forests/grass land contributed the most. the large-scale expansion of construction lands and the loss of forest or grass land, water, cultivated land, and wetland will lead to the rapid rise of LST.   In 2000, the CI of forest or grass land reached the minimum value of −3.090. The high-coverage forest or grass land absorbed and blocked solar radiation and the water vapor generated from the transpiration of plants, which promoted the exchange of energy and matter between the surface and the near-surface atmosphere. As a result, a low ground temperature can be maintained. The CI value of construction land continued increasing from 0.068 in 1993 to 0.598 in 2016, which was the main contributor to the increase of LST. The construction land expanded rapidly in the past 24 years, and the underlying surface properties have changed, which lead to a significant increase in LST and an intensification of the thermal environmental effect. The CI scores of wetland in 1993, 2000, 2008, and 2016 were −0.018, −0.016, −0.003, and 0.001, respectively. Due to the small proportion of wetland in the total study area, the absolute value of CI was small. From 1993 to 2000, the CI value of wetland did not change significantly, but all of them were fewer than 0. During this period, the coverage of vegetation and shrubbery was high within the wetland in the study area. From 2008 to 2016, the CI value of wetland changed from negative to positive, so its contribution to the increase of LST also changed from negative to positive. During this period, the urbanization process was accelerated, and the wetland around the city were greatly developed, which reduced the area of the wetland and destroyed its functions. The following is a typical example: in 2016, the Puxiazhou Wetland was completely converted to construction lands (Fuzhou Strait International Convention and Exhibition Center was completed in 2010). Since the distribution of bare land was scattered and only accounted for a small proportion of the study area, their contributions to LST in the past 24 years were fluctuating. Figure 7 shows the correlation relationship between the change of LULC type and LST. The data in 1993 and 2016 were taken as examples to conduct a regression analysis for the area proportion of land use/cover type conversion and its corresponding LST. Among them, the LULC types with a CI value of fewer than 0 were defined as the cold sources in the study area. The cold source includes forest or grass land, water, cultivated land, and wetland. The results showed that the expansion of construction land and the loss of cold sources were positively correlated with LST from 1993 to 2016. In other words, the large-scale expansion of construction lands and the loss of forest or grass land, water, cultivated land, and wetland will lead to the rapid rise of LST.

Spatial Agglomeration Effect of Land Surface Temperature
Based on Figure 4, the regions with RLST > 2 were extracted by using the GIS tool. The regions with RLST greater than 2 °C were defined as the regional heat islands [39,40], and Table 6 shows the average NNR. From Table 6, the NNR values from 1993 to 2016 were all fewer than 1 (p < 0.01) and decreased from 0.84 in 1993 to 0.74 in 2016, indicating that the LST in the study area was not randomly distributed in space. Instead, it showed a strong spatial agglomeration effect. Due to the rapid economic development of Fuzhou, the urban area sprawled like a "big cake". The built-up areas expanded extensively, and the LST agglomeration intensity became more significant. By analyzing Figure 8 and Table  7, it can be seen that the NNR values in four periods all decreased as the RLST grade increased. This result showed that the agglomeration effect of high temperature became more significant. When 2 < RLST ≤ 4, 6 < RLST ≤ 8, and RLST > 8, there was a positive correlation between the RLST area and NNR from 1993 to 2016. When 4 < RLST ≤ 6, both the RLST area and the NNR values from 1993 to 2016 showed similarities. In 1993 and 2000, when RLST > 8, the degree of spatial agglomeration of the relative ground temperature was the highest. During the period from 2008 to 2016, when 6 < RLST ≤ 8, the degree of spatial agglomeration of the relative ground temperature was the highest. Therefore, the high-temperature concentration region of RLST > 6 should be considered more in the future.
As it can be observed in Figure 9, by using five kinds of spatial resolutions and utilizing the multi-distance spatial clustering analysis function of the GIS spatial statistics tool, the L(d) function analysis was applied on LST. The spatial scale (distance increment) was set as 5 km, and the confidence interval was taken as 90%. The spatial scale values when LST appears random distribution in space [L(d) = 0] are defined as the critical values [36]. In other words, the critical values are located at the differentiation point where the agglomeration and dispersion appear. According to Figure 9 under the spatial resolution of 300 m × 300 m, when the spatial scale value was fewer than the critical value of 56 km in 2000 and 2008, L(d) was larger than 0, which indicated that the LST presented an agglomeration trend. On the contrary, when the spatial scale was larger than the critical value of 56 km and L(d) < 0, the LST showed a dispersion trend. When the spatial scale value was fewer than the critical value of 55 km in 1993 and 65 km in 2016, the surface temperature showed a trend of agglomeration. Contrarily, when the spatial scale was larger than the critical value of 55 km and 65 km, respectively, it presented a dispersion trend. From 1993 to 2016, the peak values of LST appeared at the spatial scale of 25, 20, 20, and 25 km, respectively. At this time, L(d) reached the maximum value of 6.92, 4.86, 6.09, and 8.86, respectively, indicating a relatively high concentration of LST at this scale. Overall, the spatial scale increased from 55 km to 65 km from 1993 to 2016, and the L(d) value increased from 6.92 to 8.86. These outcomes showed that the spatial agglomeration range

Spatial Agglomeration Effect of Land Surface Temperature
Based on Figure 4, the regions with RLST > 2 were extracted by using the GIS tool. The regions with RLST greater than 2 • C were defined as the regional heat islands [39,40], and Table 6 shows the average NNR. From Table 6, the NNR values from 1993 to 2016 were all fewer than 1 (p < 0.01) and decreased from 0.84 in 1993 to 0.74 in 2016, indicating that the LST in the study area was not randomly distributed in space. Instead, it showed a strong spatial agglomeration effect. Due to the rapid economic development of Fuzhou, the urban area sprawled like a "big cake". The built-up areas expanded extensively, and the LST agglomeration intensity became more significant. By analyzing Figure 8 and Table 7, it can be seen that the NNR values in four periods all decreased as the RLST grade increased. This result showed that the agglomeration effect of high temperature became more significant. When 2 < RLST ≤ 4, 6 < RLST ≤ 8, and RLST > 8, there was a positive correlation between the RLST area and NNR from 1993 to 2016. When 4 < RLST ≤ 6, both the RLST area and the NNR values from 1993 to 2016 showed similarities. In 1993 and 2000, when RLST > 8, the degree of spatial agglomeration of the relative ground temperature was the highest. During the period from 2008 to 2016, when 6 < RLST ≤ 8, the degree of spatial agglomeration of the relative ground temperature was the highest. Therefore, the high-temperature concentration region of RLST > 6 should be considered more in the future. space expansion [41]. It can be seen from Table 8 that the spatial scale has gradually decreased as the spatial resolution increases. This phenomenon illustrates that the spatial agglomeration effect of LST will be weakened because of the increase in spatial resolution. The results of multi-distance spatial clustering analysis showed that the LST in the study area was not distributed randomly but agglomerated on a certain spatial scale. Additionally, with the rapid development of urbanization, the LST will develop in a more agglomerated trend.     As it can be observed in Figure 9, by using five kinds of spatial resolutions and utilizing the multi-distance spatial clustering analysis function of the GIS spatial statistics tool, the L(d) function analysis was applied on LST. The spatial scale (distance increment) was set as 5 km, and the confidence interval was taken as 90%. The spatial scale values when LST appears random distribution in space [L(d) = 0] are defined as the critical values [36]. In other words, the critical values are located at the differentiation point where the agglomeration and dispersion appear. According to Figure 9 under the spatial resolution of 300 m × 300 m, when the spatial scale value was fewer than the critical value of 56 km in 2000 and 2008, L(d) was larger than 0, which indicated that the LST presented an agglomeration trend. On the contrary, when the spatial scale was larger than the critical value of 56 km and L(d) < 0, the LST showed a dispersion trend. When the spatial scale value was fewer than the critical value of 55 km in 1993 and 65 km in 2016, the surface temperature showed a trend of agglomeration. Contrarily, when the spatial scale was larger than the critical value of 55 km and 65 km, respectively, it presented a dispersion trend. From 1993 to 2016, the peak values of LST appeared at the spatial scale of 25, 20, 20, and 25 km, respectively. At this time, L(d) reached the maximum value of 6.92, 4.86, 6.09, and 8.86, respectively, indicating a relatively high concentration of LST at this scale. Overall, the spatial scale increased from 55 km to 65 km from 1993 to 2016, and the L(d) value increased from 6.92 to 8.86. These outcomes showed that the spatial agglomeration range of LST in the study area expanded, and the agglomeration degree also increased. These mainly resulted from the diffusion of heat islands from urban centers to suburbs [26]. In recent years, following the guidance of "eastward expansion, southward expansion, and westward expansion", the urban built-up area of Fuzhou has been constantly spreading to the periphery, and the east, west, and south have become the main directions of urban space expansion [41]. It can be seen from Table 8 that the spatial scale has gradually decreased as the spatial resolution increases. This phenomenon illustrates that the spatial agglomeration effect of LST will be weakened because of the increase in spatial resolution. The results of multi-distance spatial clustering analysis showed that the LST in the study area was not distributed randomly but agglomerated on a certain spatial scale. Additionally, with the rapid development of urbanization, the LST will develop in a more agglomerated trend. recent years, following the guidance of "eastward expansion, southward expansion, and westward expansion", the urban built-up area of Fuzhou has been constantly spreading to the periphery, and the east, west, and south have become the main directions of urban space expansion [41]. It can be seen from Table 8 that the spatial scale has gradually decreased as the spatial resolution increases. This phenomenon illustrates that the spatial agglomeration effect of LST will be weakened because of the increase in spatial resolution. The results of multi-distance spatial clustering analysis showed that the LST in the study area was not distributed randomly but agglomerated on a certain spatial scale. Additionally, with the rapid development of urbanization, the LST will develop in a more agglomerated trend.

Discussion
Urbanization significantly influenced trends of the LST. The increase of city size and the proportion of per LULC both contributed to urban warming [42][43][44]. Rapid and continuous urbanization has greatly changed the LULU of Fuzhou [45]. The population increased by almost 2 million in Fuzhou between 1993 and 2016, while the population density also increased. In addition to the increasing population sizes, construction land areas were also expanding in Fuzhou. The expansion of construction land in Fuzhou was mainly concentrated in the central urban area and the eastern coastal areas. The western region expanded to the surrounding areas with the county and city center as the origin. The degeneration of forest or grass land was mainly concentrated in the south of Minqing, the northwest of Fuqing, the coastal area of Luoyuan, and a small part of Minhou. However, forest or grass land still has the largest area in the study area.
As expected, LST responded strongly to spatiotemporal dynamics of LULC in Fuzhou. High temperatures were observed in construction land and its extent increased with time as the city was expanding. The influence of construction land (CI > 0) explained why medium-and high-temperature zones (medium-temperature zone, sub-high-temperature zone, and high-temperature zone) and the area with RLST ≥ 8 were recorded in the central urban area and established town areas of various counties. The shape of medium and high-temperature zones in Fuzhou closely mimics that of the construction land, indicating their strong warming influence. This concurs with Connors J. P. et al. [46], who observed that surface urban heat island expands with the expansion in construction land. This results in large amounts of stagnant heat and high temperatures, especially in closely packed and high-rise buildings [17]. Forest or grass land has been indicated to be a strong mitigation measure (CI < 0) against the rise of LST in Fuzhou. In Fuzhou, low-temperature zones (low-temperature zone and sub-low temperature zone) and areas with RLST < 8 remained regions that had high coverage of forest or grass land. This concurs with Cai et al. [16] who, in the Fuzhou municipality, showed that the LST reduction effect of forest or grassland increased with the percentage of the total area covered. Although water and farmland also had a cooling effect (CI < 0), their contribution had remained minimal over the years due to the low proportion that they occupied in Fuzhou.
Above all, the contribution index is scientific and convenient in quantifying the temperature effect caused by the evolution of underlying surfaces at the urban scale. The research clarified the influence process of the evolution of urban underlying surfaces on the LST and reveals the internal relationship. Meanwhile, the results can provide help to effectively cope with the UHI effect, promote the sustainable and healthy development of cities, and provide a reference for urban planning and management.
This study used Landsat remote sensing image as basic data, which has the advantages of good timing, stable data, and easy access. It has been widely used in the field of heat environment research. However, the accuracy of extracting LULC will be limited because of its spatial resolution of 30 m. With the continuous progress of remote sensing technology, more and more high-resolution images have gradually emerged. Considering higher resolution images to improve the accuracy of LULC classification and provide reference for internal details, it provides the basic conditions for the study of the surface.
The factors affecting the surface temperature are complex. In addition to the urban underlying surface, it also includes wind speed, precipitation, artificial heat discharge, landscape configuration, building height, urban structure and socioeconomic factors, which will be one of the important research directions in the future research.

Conclusions
Based on remote sensing images from 1993 to 2016, this study quantitatively investigated the temporal and spatial pattern changes and characteristics of the thermal environment in Fuzhou during the process of urbanization. The study content includes supervised classification, LST inversion, CI calculation, and multi-scale spatial model. The key points are concluded as follows: (1) From 1993 to 2016, the LULC type in the study area changed significantly. The area of construction land increased by 1231.04 km 2 , with a change rate of 295.33%. Cultivated land had been occupied by construction land. The reduction in the absolute number of wetland was small, but the decrease in proportion was high. Construction land expansion and cold source loss were both positively correlated with LST. (2) Based on the urban areas and LULC types of Fuzhou, the CI was calculated. The results showed that the urban areas of Fuzhou made a positive contribution to the rise of LST, while Minqing and Yongtai was a negative contributor from 1993 to 2016. Overall, the temporal and spatial distribution of LST contributions in Fuzhou counties and urban areas were found to be uneven. Forest or grass land, wetland, and water had cooling effects, while construction land had warming effects. Wetland's contribution to the LST increase changed from negative to positive. (3) The medium and high-temperature regions spread gradually along the axis of the Minjiang River, and the low temperature and sub-low temperature areas decreased significantly. From 1993 to 2016, the spatial agglomeration range and degree of LST were continuously expanded, and the larger spatial scale showed an increasing agglomeration trend.
The following three suggestions are given: (1) Controlling the expansion of construction land and allocating urban green lands reasonably. The expansion of construction lands has resulted in the disappearance of a large number of urban green lands. What is worse, both the expansion of construction lands and the loss of urban green lands will lead to the increase of LST. Therefore, it is necessary to control the expansion of construction lands, increase the coverage of green vegetation reasonably and rationally, and improve the rate of urban green lands.
(2) Optimizing the urban spatial layout. In the process of urban planning and construction, road traffic systems should be properly managed, and the density and height of buildings in construction lands should be controlled. In addition, more space should be left for corridors used for ventilation in the urban areas should be set aside, and the air flow rate should be accelerated. (3) Implementing the wetland ecological red line strictly. The laws and regulations of wetland protection should be improved, and the wetland protection zones should be established. People's awareness of wetland protection should be enhanced so that the stress of urbanization on wetlands can be effectively alleviated. Consequently, the sustainable utilization and development of wetland resources can be achieved.
Fuzhou is considered an important area in China because it is the China (Fujian) Pilot Free Trade Zone, the Pan-Pearl River Delta Region, and the "prior to carry and try" Economic Zone on the West Coast of the National Straits. Hence, it is urgent to carry out research on the impact of the environmental thermal effect in Fuzhou and formulate regional sustainable development countermeasures to cope with the impact of urbanization. The research results can promote the healthy development of urban areas, provide scientific decision-making references, beneficial experience, and inspiration for China. By conducting more research and obtaining more reliable experimental outcomes, more advanced urbanization development strategies can be formulated, and the harmonious development of cities and natural habitats can be promoted.

Acknowledgments:
We are sincerely grateful to U.S. Geological Survey for Landsat image and data processing, to Chinese Academy of Sciences for reference data. We greatly appreciate numerous volunteers for field survey. Furthermore, we feel indebted to three anonymous reviewers for their valuable comments and suggestions on adequately advancing the quality of this manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.   Figure A1. Location of LST sampling points in the study area.