Analyzing Thermal Environment Contribution and Driving Factors of Lst Heterogeneity Based on Urban Different Development Zones

1. Introduction

The rapid urbanization process has resulted in an increasing impervious surface in urban and its surrounding areas, further impacting the urban thermal environment [1,2]. It is obvious that the city centers have higher land surface temperatures (LST) and experience higher LST heterogeneity than suburbs surrounding in terms of rural-urban temperature differences, a phenomenon known as the urban heat island (UHI) effect [3,4,5,6]. Therefore, accurate analysis of the urban landscape pattern and thermal environment is important for mitigating the UHI and improving urban livability [7].

Impervious surface area (ISA) and vegetation cover are two crucial parameters in the analysis of LST and urban thermal patterns [6,8]. Urban expansion leads to vegetation areas reduction, impervious surfaces increase, and alterations in regional climate [9,10]. To simulate the climatology of UHI phenomena corresponding to changes in urban landscape patterns, it is essential to have data on the percentage of ISA and fractional vegetation cover (FVC) that capture the landscape evolution [11,12]. Consequently, Landsat remote sensing data with a spatial resolution of 30 m × 30 m are commonly employed to generate updated urban landscape pattern data. Due to the mixed pixel problem, traditional per-pixel classifiers are often inadequate for handling the intricate fine-scale urban landscape patterns. Soft classification techniques, such as spectral unmixing, are typically used to extract sub-pixel information for characterizing urban land cover patterns and thermal environments [8,13].

To mitigate the UHI effect and adapt urban areas to climate change, it is crucial to understand how changes in ISA and vegetation cover, as well as their composition, contribute to urban LST [14,15,16]. Previous studies have shown that the percentage of ISA and FVC provide complementary information to hard-classified land cover types for analyzing the urban thermal environment [17,18,19,20,21]. However, the effects and spatial distribution of urban ISA and vegetation on the thermal environment need to be further analyzed, particularly in areas with varying ISA and vegetation densities due to different stages of urban development [2,22,23]. Therefore, it is necessary to analyze the impact of variations in fractional covers of ISA and vegetation on the thermal environment, as well as the spatial heterogeneity or agglomeration of LST in different urban development areas [2,24,25].

In urban areas, continuous percent ISA information with a value from 0% to 100% can reveal central business districts and urban residential areas with varying densities and landscape patterns, rural developed centers and relatively undeveloped areas [26]. The value of the percent ISA increased from 0 to 100% can result in increasing LST and UHI; inversely, FVC with a value increased from 0% to 100% can decrease LST and mitigate UHI [12]. Although some studies have indicated that urban landscape pattern or land-use types significantly affect LST, few studies explored the different landscape patterns and thermal contributions of urban different development areas [1,27]. In order to effectively adapt to climate change, it is important to analyze the thermal contribution of different landscape patterns to urban functional zones and the entire urban area.

The spatial heterogeneity of LST is a common geographical phenomenon in urban expansion [2,24,28]. Analysis of the spatial distribution and spatial heterogeneity or agglomeration of the urban thermal environment in the urbanization process can help planners understand characteristics such as the direction and degree of urban expansion, and further promote the healthy development of urbanization through reasonable urban planning [29,30,31]. Geodetector has been applied to analyze the factors controlling the spatial patterns of various geographical phenomena [32]. At present, there are few studies to measure the spatial heterogeneity of LST resulting from urban expansion and its driving forces. In order to effectively guide the construction of the urban ecological environment and cope with urban climate change, it is necessary to analyze the connection between the heterogeneity characteristics of LST and the influencing factors, which are vital for the identification, prevention, and remediation of the negative impact of the urban thermal environment [11,33]. In this study, Geodetector was used to reveal the spatial influence of impact factors on LST through driving forces analysis in urban different development categories.

This study selected a study area to explore the spatial influence of urban landscape patterns on the thermal environment. A fully constrained linear spectral unmixing method was used to extract subpixel fractional cover data. Two Landsat images were selected to quantify the impact of urban landscape patterns on the urban thermal environment and the factors controlling spatial heterogeneity of urban LST. The objectives of this study were to (1) investigate the characteristics of urban landscape pattern and LST in urban different development areas; (2) explore the variations of the thermal environment by analyzing and comparing the heat effect contributions based on different development zones and different dates; (3) reveal the internal regularities affecting the spatial pattern of LST by investigating and comparing the driving factors on LST spatial distribution in urban different development areas. Results of the study can provide information for urban planners on how to mitigate the UHI effect by rational landscape design and urban management.

2. Material and Methods

2.1. Study Area and Data

Fuzhou, the capital city of Fujian province in southeastern China, is selected as the study area (Figure 1). Fuzhou has experienced rapid urban expansion due to a significant increase in population over the last 30 years. This urban expansion entailed a transformation of land cover types, shifting from natural surfaces to impervious surface areas. These changes in the urban landscape pattern have led to alterations in thermal properties, surface radiative properties, and ecosystems. The city is surrounded by mountains and is located in a subtropical humid climate zone. The vegetation in the region is predominantly evergreen, and the fraction of vegetation cover remains nearly constant throughout different seasons. Fuzhou experiences hot summers, while winters have a mild climate, and variations in climate during spring, autumn, and winter are minimal.

Bi-temporal, cloud-free Landsat images were utilized in this study, including a Landsat 7 ETM+ image (acquired on September 12, 2004) and a Landsat 8 Operational Land Imager (OLI) and thermal infrared sensor (TIRS) images (acquired on September 27, 2021). Additionally, bi-temporal, high spatial resolution images (acquired in 2004 and 2021, respectively) with a 4-meter spatial resolution obtained from the Google Earth website were used to assess the accuracy of percent ISA and FVC derived through spectral unmixing from Landsat imagery. The land surface reflectance of visible and near-infrared bands of Landsat data was derived using the FLAASH model in ENVI.

To effectively analyze and compare urban landscape and thermal patterns in different development areas, the study area was divided into three circular zones (Zone 1, Zone 2, and Zone 3) from the city center to peripheral urban areas, representing different urban development densities (as delineated by yellow lines in the right figure of Figure 1). Zone 1 comprises the urban core area with high development density, Zone 3 represents an area with low development density in the peripheral urban areas, and Zone 2 encompasses an area with medium development density situated between Zone 1 and Zone 3. Zone 2 was not considered part of the urban area in 2004 but was included in the urban extent in 2021 due to urban expansion. These three divided zones exhibited distinct development densities. The data generated from zonal statistics aids in understanding the thermal environment contributions of various urban development areas and the driving forces behind the spatial heterogeneities of LST more accurately.

2.2. Methods

The methodology for this study consists of four main steps. The first step includes LST retrieval and fractional covers retrieval by spectral unmixing. The second step includes LST classification and accuracy assessment of subpixel fractional covers. The third step explores thermal environment contribution based on different urban zones, this step includes calculation of the thermal environment contribution index (CI_i), the weighted LST unit index (WLUI) and the regional weighted LST unit index (RWLUI). The forth step analyzes driving forces of the LST heterogeneity at different zones by calculating q-values. A flowchart of the study is illustrated in Figure 2.

2.2.1. LST Retrieval and LST Classification

Bi-temporal LST data were retrieved from Landsat 7 ETM+ and Landsat 8 TIRS/OLI images. The radiative transfer equation was employed to derive LST from the Landsat data. Firstly, Top-of-Atmosphere (TOA) radiance was calculated from the bi-temporal Landsat images through radiometric calibration. Subsequently, TOA radiance values for visible/near-infrared bands were converted to surface reflectance using the atmospheric correction tool MODTRAN 4.1. Simultaneously, TOA radiance values for thermal infrared bands (band 6 of Landsat 7 ETM+ and Band 10 of Landsat 8 TIRS) were transformed to surface-leaving radiance (L_T) by applying the atmospheric correction tool MODTRAN 4.1 and inverting Planck’s law [34]:

$$\mathrm{L}\mathrm{T}=\frac{{\mathrm{L}}_{\mathsf{\lambda }}{-\mathrm{L}}_{\mathsf{\mu }}-\mathsf{\tau }\left(1-\mathsf{\epsilon }\right){\mathrm{L}}_{\mathrm{d}}}{\mathsf{\tau }\mathsf{\epsilon }}$$

(1)

where L_λ is the TOA radiance image for the thermal infrared band, Lμ, L_d and τ are the upwelling radiance, downwelling radiance and transmission respectively. The atmospheric correction tool MODTRAN 4.1 was used to estimate the values of L_μ, L_d and τ. ε is the surface emissivity calculated according to the methods by [35,36].

Finally, LST were calculated by the Landsat-specific estimate of the Planck curve [37]:

$$\mathrm{L}\mathrm{S}\mathrm{T}=\frac{{\mathrm{k}}_2}{\ln \left(\right({\mathrm{k}}_1/{\mathrm{L}}_{\mathrm{T}})+1)}$$

(2)

where LST is the temperature in Kelvin (K), k₁ and k₂ are the pre-launch calibration constants. For Landsat 7 ETM+, k₁ = 666.09 W /(m²srμm) and k₂ = 1282.71 K. For TIRS Band 10, k₁= 774.89 W/(m²srμm), k₂= 1321.08 K.

The resolution of the visible/near-infrared bands of Landsat 7 ETM+ image and Landsat 8 OLI image is 30 m while band 6 of Landsat 7 ETM+ image and Band 10 of Landsat 8 TIRS image) are collected at a resolution of 60 m and 100 m respectively, the thermal infrared data were resampled to 30 m in this study. Figure 3a and Figure 3b are LST retrieved from the ETM+ image acquired on September 12, 2004 and LST retrieved from the TIRS image acquired on September 27, 2021 respectively.

The LST of a same area may show some differences between different dates even in same season because the multi-temporal images were acquired in different atmospheric conditions. In order to reduce the influence of such LST fluctuations, the four level LST classes were created by adding or subtracting standard deviation (SD) from mean LST of the whole study area:

LST₁: LST ⩽ (LST_mean − 2 SD)

LST₂: (LST_mean − 2 SD) < LST ⩽ LST_mean

LST₃: LST_mean < LST ⩽ (LST_mean + 2 SD)

LST₄: LST > (LST_mean + 2 SD) (3)

The LST increase from LST₁ to LST₄ and the LST₄ area is defined as a high temperature area. The classification of four levels LST can be used to analyze the evolution of the thermal environment. Four LST classes are calculated as in Eq. (3) and shown in Figure 3c and 3d.

2.2.2. Percent ISA and FVC Derived by Spectral Unmixing of Fully Constrained Least Squares and Accuracy Assessment

Continuous percent ISA information can reveal areas of urban development with varying densities and patterns. As one of the soft classification methods, linear spectral unmixing by fully constrained least squares (FCLS) is frequently adopted to extract sub-pixel land cover types for remote sensing imagery, especially for Landsat imagery. In this study, FCLS was used to extract sub-pixel ISA and vegetation of Landsat imagery. The FCLS method assumes that the reflectance of each pixel at each band is a linear combination of the reflectance of all endmembers within the pixel [38]. The spectral reflectance in band i can be described as follows:

${\mathrm{R}}_{\mathrm{i}}$= $\sum _{\mathrm{k}=1}^{\mathrm{n}}{\mathrm{f}}_{\mathrm{k}}{\mathrm{r}}_{\mathrm{i}\mathrm{k}}+{\mathsf{\epsilon }}_{\mathrm{i}}$(4)

$\sum _{\mathrm{k}=1}^{\mathrm{n}}{\mathrm{f}}_{\mathrm{k}}=1{\mathrm{f}}_{\mathrm{k}}\ge 0 \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{k}=1,2,3,\dots ,\mathrm{n}$(5)

where n is the number of endmembers, f_k the fraction of endmember k within the pixel, r_ik the reflectance of endmember k in band i, and ε_i is the residual error for band i. Eq. 5 is the constraint conditions for spectral unmixing, the fractions of one pixel must sum to 1 and all fractions must be ≥ 0.

Firstly, the water body was masked out from the images of the study area, and two-dimensional spectral mixing spaces were used to extract endmembers of vegetation, impervious surface (high- and low-albedo impervious surfaces), and soil. Then, applying the FCLS technique, the fractional cover of each component was estimated by Eqs. 4 and 5. Finally, the high- and low-albedo ISA classes were aggregated to total percent ISA.

Figure 4 displays the fractional covers of ISA and vegetation derived by FCLS from the Landsat data for 2004 and 2021. Because soil areas are limited and primarily distributed alongside the river, and soil does not have a significant effect on the analysis of the thermal environment impact of percent ISA and FVC, fractional soil is not shown here. With the lowest values (0) in black and the highest values (100%) in white, the fractional images of Figure 4 illustrate urban land cover composition and their changing patterns. The areas with a higher percent ISA value are those with more developed land and high-density residential areas. ISA and vegetation cover vary noticeably from zone 1 to zone 3. Generally, ISA increased at the periphery of the city due to urban expansion (Figure 4a and Figure 4c) with FVC decreasing (Figure 4b and Figure 4d) from 2004 to 2021. However, Figure 4b and Figure 4d show that the percentage of fractional vegetation for some pixels in zone 1 was higher in 2021 than those in 2004. One reason for this phenomenon is the greater attention to ecological landscape design in some developed areas, which resulted in more vegetation cover interspersed in the urban core area after 2004.

As per our previous studies [2,18], the accuracy of fractional covers derived from ETM+/OLI imagery by FCLS was assessed by comparing the cumulative area of fractional cover for ISA and vegetation in each of the selected test areas with the reference data extracted from the high-resolution Google Earth images. Area accumulation was carried out by multiplying the percent ISA with a pixel area of 900 m² (30 m * 30 m). Accuracy was assessed using five sampled plots shown in Figure 1 to represent the overall accuracy of the entire study area. Table 1 displays the results of the comparison between the reference areas obtained from high-resolution images and the corresponding areas derived from the Landsat ETM+/OLI imagery by FCLS. The visual interpretation method was employed to extract land cover types from the high-resolution data. This approach was feasible because the acquisition dates of high-resolution images were close to the dates of the corresponding ETM+/OLI imagery, and land cover changes were minimal between the reference data and corresponding Landsat imagery.

Table 1 indicates good agreement for ISA and vegetation cover between five sampled sites derived from Landsat and high-resolution images. Comparing with the reference data obtained from high-resolution images, the area of ISA and vegetation derived by FCLS from Landsat images in five test sites shows only small differences. In Table 1, the difference of ISA derived from ETM+/OLI imagery and corresponding high-resolution images was ⩽8.2% for five sites in 2004 and ⩽7.16% for five sites in 2021, and the difference of vegetation was ⩽3.99% and ⩽1.64% in 2004 and 2021 respectively.

In general, the accuracy of the vegetation fractions in the two test sites was higher than for ISA in both dates because the vegetation spectra can be more easily distinguished from other land cover types than those of ISA. Table 1 showed that the accuracy of five test sites in both dates was all less than 8.5% different between fractional covers of ISA and vegetation derived from ETM+/OLI imagery and those from high-resolution images. It showed that fractional covers derived by FCLS were reliable. Considering the heterogeneity of the urban landscape pattern and spatial resolution differences between Landsat imagery and high spatial resolution reference images, Table 1 showed that fractional covers of ISA and vegetation derived by FCLS are a good representation of reality in the study area.

In some areas, land cover change is taking place and detectable at the sub-pixel scale but cannot be detected at the pixel scale by a crisp classification. Fractional cover derived by soft classification can be used to measure the magnitude of change and further to quantify the urban landscape and thermal environment patterns.

2.2.3. Calculation of Thermal Environment Contribution Based on Different Urban Zones

Urban development areas exhibit varying landscape patterns and LST, thus contributing different heat effects to the urban thermal environment [19]. It is necessary to further identify the influence of distinct urban zones on the urban thermal environment and analyze their variations across different dates. This analysis is invaluable for urban planners in determining how to mitigate the UHI effect through judicious landscape design and urban management.

In this study area, the entire region was divided into three concentric zones (Zone 1, Zone 2, and Zone 3) based on differing urban development densities, as shown in Figure 1. Differing from the thermal environment contribution index (CI_i) calculated by [39] for various land cover types, in this study, the CI_i, the weighted LST unit index (WLUI), and the regional weighted LST unit index (RWLUI) were employed to evaluate the heat contribution effects across different zones of the study area. The urban CI_i was calculated as follows.

${\mathrm{C}\mathrm{I}}_{\mathrm{i}}=\frac{{\mathrm{C}\mathrm{I}}_{\mathrm{i}}^{\mathrm{*}}}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{C}\mathrm{I}}_{\mathrm{i}}^{\mathrm{*}}}\times 100\mathrm{\%}$(6)

${\mathrm{C}\mathrm{I}}_{\mathrm{i}}^{\mathrm{*}}=\frac{\sum _{\mathrm{j}=1}^{{\mathrm{n}}_0}\left(\mathrm{L}\mathrm{S}{\mathrm{T}}_{\mathrm{i}\mathrm{j}}-{\mathrm{L}\mathrm{S}\mathrm{T}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}\right)\times {\mathrm{n}}_0}{{\mathrm{L}\mathrm{S}\mathrm{T}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}\times \mathrm{N}}\times 100\mathrm{\%}$(7)

where ${\mathrm{C}\mathrm{I}}_{\mathrm{i}}^{\mathrm{*}}$ is the initial thermal environment contribution for the ith zone, LST_ij is the LST of the jth pixel in the ith zone with values that are larger than the mean LST of the whole study area (LST_mean), n₀ is the number of pixels with LST values larger than the value of LST_mean in the ith zone, and N is the number of pixels of the whole study area.

The WLUI describes the percentage of pixels with LST values larger than the LST_mean for the ith zone, whereas the RWLUI describes the percentage of the LST of the ith zone with values larger than the LST_mean of the whole study area. The two indices were calculated using

$\mathrm{W}\mathrm{L}\mathrm{U}\mathrm{I}=\frac{{\mathrm{n}}_{\mathrm{i}}}{{\mathrm{n}}_{\mathrm{z}\mathrm{o}\mathrm{n}\mathrm{e}\mathrm{i}}}\times 100\mathrm{\%}$(8)

$\mathrm{R}\mathrm{W}\mathrm{L}\mathrm{U}\mathrm{I}=\frac{{\mathrm{n}}_{\mathrm{i}}}{\mathrm{N}}\times 100\mathrm{\%}$(9)

where n_{zone i} is the number of pixels for the ith zone, and i is 1, 2 or 3 in this study. N is also the number of pixels of the whole study area as in Eq. 7.

In order to accurately describe the impacts of different landscape patterns on urban thermal environment, it is necessary to calculate CI_i, WLUI and RWLUI as Eqs. 6, 8 and 9 to further analyze the thermal environment impact of urban different development zones.

2.2.4. Spatial Statistical Analysis with Geographical Detector (Geodetector)

Spatial heterogeneity can reflect the spatial distribution of geographical phenomena and explain its decisive factors and their links [40,41]. The analysis of spatial heterogeneity can reveal the evolution of geographical features. In this study, the analysis of spatial heterogeneity was used to reveal the LST heterogeneity of urban expansion and its driving mechanism.

Geodetector is a statistical tool used to quantitatively interpret geospatial heterogeneity and the driving forces of spatial differentiation. Through the constructed Geodetector model, the influence of different factors on geospatial phenomena and spatial differences can be calculated. This is measured by the q-statistic [40,42] as follows.

q = 1 − $\frac{{\sum }_{\mathrm{H}=1}^{\mathrm{L}}{\mathrm{N}}_{\mathrm{h}}{\mathsf{\sigma }}_{\mathrm{h}}^2}{\mathrm{N}{\mathsf{\sigma }}^2}$ = 1-$\frac{\mathrm{S}\mathrm{S}\mathrm{W}}{\mathrm{S}\mathrm{S}\mathrm{T}}$

SSW =${\sum }_{\mathrm{H}=1}^{\mathrm{L}}{\mathrm{N}}_{\mathrm{h}}{\mathsf{\sigma }}_{\mathrm{h}}^2$, SST = Nσ² (10)

where N_h and N are the number of samples in sub-region h (hth stratum) and the entire region, respectively; h = 1, 2, …, L is the stratification of impact factor X; σ² and σ_h² are the total variance of Y across the entire study area and the sub-region h, respectively. SSW and SST are the sum of the variance in each layer and the total variance of the region, respectively. The value of q ranges from 0 to 1, the larger the value, the higher explanatory power or influence of X for the spatial differentiation of Y [43].

There are four detectors included in Geodetector, namely the factor detector, interaction detector, risk detector, and ecological detector. The factor detector and interaction detector are applied in this study. The factor detector can explore and determine to what extent a factor X explains the spatial differentiation of Y. Factors that represent or control the spatial pattern of urban land use should play an important role in determining the spatial distribution of LST. In this study, we use the factor detector to detect the spatial heterogeneity of LST and identify the key driver factors behind the variation in LST. It is represented by q calculated as Eq. 10; the q value is used to measure the influence of each driving factor on LST. The greater the q value, the greater the influence of the factor on LST, and the more obvious is the influence on LST spatial differentiation. The interaction detector can distinguish and determine the interaction types between different factors X_i on Y, and evaluate whether the joint action of factors X₁ and X₂ (q(x₁∩ x₂)) will strengthen or weaken the explanatory power of dependent variable Y, or whether their influences on Y are independent from each other.

Geodetector can be applied to identify the driving force behind elements by detecting their spatially layered heterogeneity [44]. In this study, three different factors X_i are percent ISA, FVC, and elevation; the dependent variable Y is LST. Geodetector requires explanatory variables to be discretized; hence each factor is pre-discretized into 5 intervals by the natural breakpoint method before analysis.

3. Results and Discussion

3.1. Analysis of Spatial Pattern of LST at Different Urban Zones

In this study, the spatial pattern of LST and its variation were analyzed based on different zones between bi-temporal images. Figure 3c and Figure 3d show that high LST (LST₄) and sub-high LST (LST₃) areas (areas with > LST_mean) were mainly distributed in the central (zone 1) and some areas of zone 2 in 2004. Zone 1 and part of zone 2 were assembled with ISA of buildings and roads, resulting in high and sub-high LST areas on both dates. Meanwhile, the distribution of LST₄ and LST₃ areas in zone 2 and zone 3 in 2021 was notably larger than in 2004. In zone 2 and zone 3, the extent of spatial distribution of LST₄ and LST₃ areas significantly increased from 2004 to 2021, especially in zone 2 and the western part of zone 3. This indicates that urban expansion resulted in an increase in LST primarily in zone 2 and zone 3 between the two dates.

The medium (LST₂) and low (LST₁) LST areas were mainly distributed in zone 2 and zone 3 of peri-urban areas in both images, especially zone 3. This was limited by terrain and topography because LST₂ and LST₁ areas were characterized by forested and water-covered regions in mountainous areas. Urban expansion led to a transformation in land cover from natural to ISA, mainly occurring in the south and west plain areas of zone 2 and zone 3 from 2004 to 2021. The land cover types in the northern areas of zone 2 and zone 3 changed little from 2004 to 2021 due to mountainous constraints; consequently, the medium and low LST areas experienced minimal change in these regions during the period.

To further quantify the impact of high and sub-high LST areas on the thermal environment, we need to analyze the change in areal proportion of high and sub-high LST areas with urban expansion. Table 2 shows the areal proportions of the LST₄ and LST₃ categories in each zone and the total LST₄ and LST₃ area of the study area. As shown in Figure 3, LST₄ and LST₃ areas dominated zone 1, accounting for 94.99% in 2004. However, the areal proportions of LST₄ and LST₃ areas in zone 1 decreased slightly to 92.98% in 2021, likely due to increased attention to urban ecology, consequently, greater vegetation cover interspersed in some areas of zone 1 after 2004. For zone 2 and zone 3, the areal proportions of LST₄ and LST₃ areas increased from 56.11% and 21.08% to 62.03% and 32.49%, respectively, from 2004 to 2021. This indicates that the high and sub-high LST areas exceeded LST_mean in this study area resulting from urban expansion primarily increased in zone 2 and zone 3 during the period.

As for the areal proportions of high and sub-high LST areas in each zone relative to the total high and sub-high LST areas of the entire study area, the areal proportions of LST₄ and LST₃ areas in zone 1 decreased from 26.58% to 20.78% during the period. In zone 2, areal proportions also decreased from 34.97% to 31.04%. However, the areal proportions of LST₄ and LST₃ areas in zone 3 relative to those of the entire study area increased from 38.45% to 48.18% in the same period. This indicates that, with urban expansion, high and sub-high LST areas mainly increased in zone 3, followed by zone 2. Although the area of high and sub-high LST also increased in zone 2 as shown in Figure 3, the rate of increase was slower compared to that of zone 3 when considering the entire study area (Table 2).

Obviously, zone 1 had the highest LST in general, followed by zone 2, and then zone 3. The high and sub-high level LST areas were heat sources contributing to the UHI effect. The lowest LST areas are primarily vegetation areas, which act as heat sinks. Although zone 2 has a cooling thermal environment due to its vegetation area, its cooling effect is limited due to the small area of vegetation. In zone 3, there are numerous vegetation areas that serve as major heat sinks, playing the most significant role in cooling the thermal environment in the study area. Additionally, the lowest LST areas also occur in the water bodies of all three zones.

3.2. Analysis of Thermal Environment Contribution at Each Urban Zone

The quantitative analysis of the thermal environmental impacts of different urban landscape patterns becomes important in urban ecology because of urban thermal environment problems. In this study, zone 1 obviously shows a higher LST, followed by zone 2 and zone 3. However, when considering the heat effect of each zone on the thermal environment and their changes between different dates, it needs to be further analyzed. In order to further quantify the heat contribution in detail, Eq. 6 was used to calculate the CI_i and to further analyze the thermal environmental impact of different zones in this study (Figure 5).

Urban expansion with ISA clearly resulted in an increase in LST. As shown in Figure 5, the CI_i of zone 1 and zone 2 were 24.84% and 33.44% in 2004 respectively. With the urban expansion resulting in LST variation, the CI_i of zone 1 and zone 2 decreased to 10.61% and 25.91% in 2021 respectively. However, the CI_i of zone 3 increased from 41.72% in 2004 to 63.49% in 2021. Though ISA increased in zone 2 and zone 3 in the study area, and the area of ISA with >LST_mean also increased in zone 2 from 2004 to 2021 (Table 2), the increase rate of ISA with >LST_mean was higher in zone 3 than in zone 2 in the period. This resulted in the decrease of CI_i in zone 2. Figure 5 indicates that urban expansion, coupled with areas exhibiting >LST_mean, has led to thermal environment issues, predominantly escalating in zones 2 and 3 during this period, particularly in zone 3. Zone 3 and zone 2 exhibited relatively high CI_i values on the two dates compared to zone 1, with zone 3 being the primary source of the heat effect in the study area and zone 2 in second place. Although CI_i values of zone 1 obviously decreased, the area with >LST_mean decreased only slightly between the two dates (Table 2). The increase of ISA with >LST_mean was more pronounced in zone 3 and zone 2 than in zone 1 from 2004 to 2021.

Compared with zone 1, the CI_i values of zone 2 decreased relatively little between the two dates. From 2004 to 2021, the CI_i of zone 3 increased significantly and contributed 63.49% of the heat effect in the entire study area in 2021, whereas the influence of zone 3 on the thermal environment was not significant in 2004. The thermal environment impact of zone 3 increased significantly in 2021 due to urban expansion. The CI_i of zones 2 and 3 reached from 75.16% in 2004 to 89.40% in 2021, indicating that the heat effect contribution was determined and dominated by zones 2 and 3 in the study area.

The different urban development zones exhibited distinct landscape patterns and LST features. To further quantify the impact and pattern of the thermal environment, WLUI and RWLUI were calculated based on Eqs. (8) and (9) to assess the thermal contribution effects of different zones in the study area (Figure 6). WLUI represents the percentage of pixels with values exceeding the LST_mean for each zone. Figure 6 showed that WLUI of zone 1 was high, accounting for approximately 93.07% and 93.22% from 2004 to 2021, with little change between the two dates. WLUI of zone 1 indicated that areas with high and sub-high LST (>LST_mean area) dominated in zone 1, with the areal ratio reaching over 93% in the urban core zone with high development density. Meanwhile, the WLUI values for zone 2 and zone 3, respectively decreased in sequence from zone 1 to zone 3 on both dates. Comparing WLUI values on both dates, WLUI values in zone 2 and zone 3 increased from 60.91% and 26.41% in 2004 to 65.25% and 34.25% in 2021 respectively. This reflects a higher increased rate of increase in WLUI in zone 3 compared to zone 2, indicating an amplified thermal environment impact in both zones from 2004 to 2021. The WLUI values of zone 2 and zone 3 indicate that >LST_mean areas increased in areas of medium development density and low development density with expansion in the period.

The RWLUI is the percentage of the high and sub-high LST area of each zone to the whole study area (Figure 6). The RWLUI values of zone 1 were 9.18% in 2004 and 9.21% in 2021. From the view of the whole study area, the thermal environment impact of zone 1 was nearly not strengthened from 2004 to 2021 because urban expansion mainly existed in zone 2 and zone 3 between the two dates. The RWLUI values for zone 2 and zone 3 continued to increase from 13.16% and 18.11% in 2004 to 14.09% and 23.47% in 2021 respectively.

The increase in thermal environment contribution indicates that high and sub-high LST areas increased with the rise in high percent ISA resulting from urban expansion. Figure 5 and Figure 6 indicated that, though the areal proportion of high and sub-high LST areas was dominant in zone 1 and the lowest in zone 3, the contribution of zone 3 to regional heat in the study area was higher than that of zone 2 and zone 1. This is because zone 3 contributes the most high and sub-high LST area among the three zones, followed by zone 2 and zone 1. Our analysis revealed the characteristics of the thermal environment effect and their variation in different development density areas at two different time points. Through analyzing CI_i, WLUI, and RWLUI of each zone and their variations, we can quantify the thermal environmental contribution of different development density areas with urban expansion. Different from other analysis of thermal environmental contribution based on different land cover types, this analysis is significant for urban planning, UHI mitigation, and urban ecology.

3.3. Driving Force Detection of Spatial Heterogeneity of LST by Geodetector

The spatial distribution and aggregation of LST were different in different urban development areas. The spatial distribution pattern of LST is related to the degree of urban expansion and impacted by different land cover types, the combination and configuration of various land cover types, and other factors. Detection of driving forces behind urban LST heterogeneity is crucial for alleviating UHI and adapting to urban climatic change. The spatial heterogeneity of LST results from the spatial distribution of different land cover types and other factors [45]. In this study, percent ISA, fractional vegetation cover, and elevation were selected as independent variables (X) and were all categorized into five types based on the Jenks natural breaks classification method. LST was selected as the dependent variable (Y). The influences of X on the spatial differentiation of Y in different zones in 2004 and 2021 were calculated using Geodetector and further comparatively analyzed.

In this study, the factor differentiation detector was used to reveal the impact of each factor on the spatial distribution of LST. Table 3 shows the q-values of three single factors on the spatial distribution of LST, and all were significant (p < 0.01). In the study area, the degree of influence of the three driving factors on the spatial distribution of LST varied in different development densities of the urban area. The q-values of three influencing factors in zone 1 were relatively low (⩽0.309) on both dates, and variations were also small in zone 1 between the two dates (Table 3). This indicates that the urban landscape pattern changes and urban expansion nearly did not exist in the urban core zone 1 from 2004 to 2021. Among the influencing factors, the FVC had a greater impact on LST in zone 1 in 2004. From 2004 to 2021, the q-value of FVC decreased from 0.309 to 0.144 in zone 1. According to the principle of Geodetector, the larger the q-value, the greater the driving force of the factor on LST. Table 3 indicates that the three factors had a lesser degree of influence on the heterogeneity of LST, and the spatial distribution patterns of LST varied minimally in zone 1 between the two dates.

The three factors had a greater influence in zone 2 and zone 3 than in zone 1 on both dates. When comparing the q values of the three influencing factors between the two dates, the q values in zone 2 and zone 3 in 2021 were larger than those in 2004, except for the q value of FVC at zone 2 in 2004, which was slightly larger than that of zone 2 in 2021. The results from the factor detection analysis showed that the largest q values were 0.581 and 0.653 in zone 2 in 2014 and 2021 respectively, and 0.494 and 0.681 in zone 3 in 2004 and 2021 respectively. This indicated that FVC had the most significant influence on LST heterogeneity in zone 2 and zone 3 in 2004, while elevation had the greatest influence on LST heterogeneity in zone 2 and zone 3 in 2021. This is because the landscape pattern changes due to urban expansion were influenced by topography (elevation) in zone 3 in 2021, which had a certain impact on the q value of zone 3. The elevation was ranked next, with values of 0.558 and 0.432 in zone 2 and zone 3 in 2004 respectively, and FVC was ranked next to elevation, with values of 0.572 and 0.602 in zone 2 and zone 3 in 2021 respectively.

When comparing q values between zone 2 and zone 3, the q values of the three influencing factors were smaller in zone 3 than in zone 2 in 2004. However, the q values of the three influencing factors were larger in zone 3 than in zone 2 in 2021. This is related to the change in urban landscape patterns in zone 2 and zone 3 resulting from urban expansion. FVC and elevation had a higher influence degree than percent ISA in zone 2 and zone 3 in 2004. The explanatory power (q value * 100%) of FVC and elevation decreased from 58.1% and 55.8% in zone 2 to 49.4% and 43.2% in zone 3 respectively, in 2004, and increased from zone 2 to zone 3 in 2021. This indicates that urban expansion also resulted in the change of landscape patterns in some higher elevation areas of zone 3. FVC and elevation had a higher influence degree on spatial heterogeneity of LST in zone 2 and zone 3 in both dates. However, the explanatory power of percent ISA significantly increased from 34.9% to 52.7% in zone 2 and from 22.4% to 54.3% in zone 3 from 2004 to 2021. Though the q values of FVC and elevation in zone 2 and zone 3 were larger than that of percent ISA in 2004, the q values of the influencing factor percent ISA significantly increased from 2004 to 2021 with urban expansion.

According to the above analysis, the three factors have weak explanatory power for the spatial heterogeneity of LST in zone 1 compared to zone 2 and zone 3. This means that the change in landscape pattern resulting from urban expansion increased the influence degree of the three driving factors on LST significantly in zone 2 and zone 3 compared to zone 1. In different urban development areas, the land cover types and their spatial distribution are different, which influences the value and spatial distribution of LST. In zone 1 of the urban core area, a high percent ISA is distributed homogeneously, resulting in the homogeneity of spatial distribution of high LST. Therefore, the influence degree of the three factors on LST heterogeneity is low. However, in zone 3 and zone 2, the three factors have a higher influence degree on the spatial distribution of LST. This indicates that they are the main driving factors affecting the LST in zone 2 and zone 3, especially the elevation factor in zone 3 in 2021. However, no one single factor dominantly affects LST distribution in zone 2 and zone 3 of the study area. Through analyzing and comparing q values of different zones, we can understand the driving force of LST distribution in different development areas, which will help in constructing a good urban ecological environment.

The spatial distribution of urban landscape and thermal patterns is usually the result of multifactor interactions. The influence of the interaction between different factors may be stronger than that of a single factor when each factor interacts with other factors on LST. Therefore, it is necessary to analyze the influence of the interaction between any two factors on the spatial difference of LST. Through interactive calculation of different driving factors in Geodetector, we analyzed the interactive detection of various driving factors on LST distribution for both dates. The q-value is defined as q (X₁∩X₂) and the results are shown in Table 4. Comparing with the results of single factor detection in Table 3, Table 4 shows that the q value increases significantly when each driving factor interacts with other factors, and its influence on LST increases, indicating that LST is a geographical pattern driven by multi-factor interaction.

The interaction detector creates a new stratum by overlaying factors X1 and X2to further compare q(X₁∩X₂) with q(X₁) and q(X₂) [46]. The results of interaction detection can show whether the combined effects of the evaluation indicators X1and X2 enhance, bi-enhance, nonlinearly enhance, weaken, bi-weaken, or nonlinearly weaken the explanatory power of LST heterogeneity, or whether the combined effects on LST are independent of each other. In this study, the explanatory power of the interactions of different factors in Table 4 was larger than that of the single factor in Table 3, and all the interactive factors showed nonlinear enhancement in 2004 and 2021. This indicates that the interaction between any two factors of percent ISA, FVC, and elevation amplifies their impact on LST distribution in three zones of two dates. The q values of interaction detection were significantly larger than those of single factor detection in Table 3 on both dates. This means that interactions of different factors have higher explanatory power on the LST heterogeneity at each zone than those of a single factor. The results from interaction detection analysis showed that the q values of interactions of three factors at zone 1 were lower than those of zone 2 and zone 3 (Table 4). This trend is the same as that in Table 3. This is because the landscape patterns of zone 1 are significant different from those in zone 2 and zone 3. Though the q values of interactions of three factors increased at zone 1 in both dates, the q values of interactions of three factors also significantly increased at zone 2 and zone 3 in both dates when compared with those of the single factor in Table 3.

As shown in Table 4, the results from interaction detection analysis showed that FVC ∩ EL and percent ISA ∩ EL at zone 2 were 0.743 and 0.749 (74.3% and 74.9% explanatory power on the LST heterogeneity) in 2004 and 2021, respectively, and were the most significant influences on LST heterogeneity at zone 2. The interactions of percent ISA ∩ EL and percent ISA ∩ FVC were ranked next at zone 2 in 2004 and 2021, with values of 0.637 and 0.591 in 2004 and 2021, respectively. As for zone 3, FVC ∩ EL was the most significant influence on LST heterogeneity in both dates; the q value was 0.621 and 0.767 in 2004 and 2021, respectively. Percent ISA ∩ FVC with a value of 0.510 and percent ISA ∩ EL with a value of 0.757 were ranked next at zone 3 in 2004 and 2021, respectively. Based on the analysis of the interaction between two factors, we can determine which kind of interaction has the strongest influence on the spatial heterogeneity of LST. This can help guide the composition and configuration of the driving factors in different development areas of urban for UHI.

In the study area, it was found that the influencing factors were different in three zones. Synthetically evaluating the results of factor and interaction detection, it is worth noting that the landscape pattern of urban expansion in zone 3 was limited by geomorphologic characterization, and elevation leads to an imbalanced distribution of LST in zone 3. When urban expansion was intense in zone 3, it is clear that most of the new urban area was distributed around the mountain. This led to an imbalance in urban expansion and also had a certain impact on the q value of zone 3. In addition, as a single factor, FVC had a high influence on LST heterogeneity in zone 2 and zone 3 of the study area.

Linear correlation cannot analyze the consistency of spatial distribution between two geographic variables. Compared with linear correlation analysis, geodetector can reveal the causal relationship between independent and dependent variables [46]. In this study, the Geodetector can show the driving force of different factors on the heterogeneity of LST resulting from urban landscape expansion through the q value. Spatial heterogeneity of LST may come from an imbalance in urban expansion. This analysis can reveal this imbalance and further help to understand the influencing factors of LST in different development areas of urban. It can provide significant information for urban planning and coping with urban climatic change.

4. Conclusions

Urban expansion changes the urban landscape pattern and also the spatial distribution of LST, and may have an effect on human health. This study aimed at exploring the impact of urban expansion on thermal environment at different development density areas. Bi-temporal Landsat images were used to derive LST, percent ISA, and FVC, and CI_i, WLUI, and RWLUI were further calculated to analyze the thermal environment contribution based on different urban zones. Geodetector was also used to explore the driving factors controlling the spatial heterogeneities of LST. This study compared the CI_i, WLUI, RWLUI_, and q between different urban zones and different dates, by revealing the thermal contribution of urban different landscape patterns and the factors controlling spatial heterogeneities of LST in urban different development areas. Though the research is focused only on one city, the analytical method can provide insights for quantifying urban landscape patterns and thermal environment in other cities. This could help alleviate negative effects of urbanization and also provide a reference for urban landscape design in other cities with different development areas. Results of the study can help mitigate the UHI effect and adapt to climatic change in cities. In this study, four conclusions were as follows:

(1) ISA increase resulted from urban expansion mainly existed in the south and west areas of zone 2 and zone 3 from 2004 to 2021, and the high and sub-high LST areas also mainly increased in zone 2 and zone 3 in the period.

(2) The thermal environment contribution of zone 3 significantly increased from 2004 to 2021. The CI_i of zones 2 and 3 reached over 89% in 2021, and the thermal environment contribution was determined and dominated by zone 2 and zone 3 in the study area in 2021.

(3) In the study area, driving factors affecting the spatial differentiation of LST are regionally different, and the influence degree of three factors on LST is weaker in zone 1 than those in zone 2 and zone 3. Compared with a single factor, the interaction between different factors can increase the influence on LST.

(4) The combined analysis of thermal environment contribution and Geodetector provides a suitable method to reveal the impact of the urban landscape pattern on urban thermal environment and driving mechanisms on LST heterogeneity. The analytic results can provide a profound understanding of the mechanisms of the UHI and provide recommendations for mitigating UHI.

In this study, Landsat images in summer were used to reveal thermal contribution and heterogeneous analysis of LST; all four seasons should also be considered in analysis in the future. In addition, we believe that the spatial scale of environmental factors is important in affecting LST, and heat contribution and LST heterogeneity may be different at different spatial scales. Therefore, multi-scale heat contributions, environmental factors, and LST should be taken into account in further comparative analysis in the future.