Spatial Pattern Impact of Impervious Surface Density on Urban Heat Island Effect: A Case Study in Xuzhou, China

: Impervious surfaces (IPS) are the major source of urban heat island effect (UHI), and the relationships between IPS and land surface temperature (LST) have been widely studied. However, the spatial impact of landscape patterns of patches with different IPS density (IPSD) on the thermal environment remains largely unexplored. Based on three Landsat 8 images of the Xuzhou built-up area obtained in May and the corresponding ground observations from 2014 to 2020, the IPSD and LST maps were inversed through a linear spectral mixture analysis and mono-window algorithm, respectively. The landscape patterns of the ﬁve IPSD levels were characterized by four landscape-level and ﬁve class-level metrics. Finally, the spatial correlation between all landscape metrics and LST were analyzed using bivariate Moran’s I. The results were as follows: (1) The ﬁndings revealed that for the landscape-level metrics, LST had signiﬁcant positive spatial correlations with Shannon’s diversity index (SHDI), Shannon’s evenness index (SHEI), and patch density (PD), while showing a signiﬁcant negative correlation with contagion index (CONTAG), indicating that increasing the types, even distribution degree, and density of patches, or decreasing the aggregation degree of the ﬁve IPSD levels will lead to the enhancement of the thermal environment. (2) Furthermore, the class-level metrics of each IPSD level, percentage of landscape (PLAND), largest patch index (LPI), landscape shape index (LSI), aggregation index (AI), and patch cohesion index (COHESION) showed signiﬁcant correlations and LST, which signiﬁed that the spatial characteristics of patch proportion, predominance degree, shape complexity, aggregation degree, and natural connectivity degree of each IPSD level are important factors affecting UHI. In addition, the spatial correlations between LST and class-level metrics were signiﬁcantly positive for IPSD levels 4 and 5 with an evidently higher Moran’s I value, indicating that landscape patterns of IPSD levels 4 and 5 were the key factors in UHI enhancement. Furthermore, the impact weights of each class-level metric of IPSD levels 4 and 5 on LST were also analyzed by applying the principal component analysis and the multivariate regression standardization coefﬁcient. These results reveal the importance and impact mechanism of the IPSD spatial patterns on UHI evolution, which may provide a valuable reference for future urban planning and climate management. This study also suggests that regional UHI can be mitigated by reducing the area proportion, natural connectivity, and shape complexity of high-density impervious surfaces.


Introduction
Rapid urbanization has promoted socioeconomic development and infrastructure upgrades [1][2][3][4], but it is also accompanied by several environmental problems, including air pollution [5], vegetation destruction [6], energy crisis [5], agricultural land occupation [7] and the urban heat island (UHI) effect [8]. UHI is a phenomenon wherein urban areas experience higher temperatures than the surrounding non-urbanized areas [9,10]. Rapid urbanization has subjected natural landscapes to great pressure, eventually transforming them into construction zones [11]. Most built-up areas covered by impervious surfaces (IPS) can dramatically change the surface radiation, thermal properties, and humidity in urban areas, resulting in the UHI phenomena [9,10,[12][13][14]. Elevated temperatures in urban areas caused by UHI can aggravate energy and water consumption [15,16], urban air pollution [17], and human health risks [18]. Therefore, the urban thermal environment has attracted research interest from the fields of urban ecology, environment, and climate [19], and effective methods to mitigate the risks and negative consequences of the UHI effect are urgently needed.
As a large number of natural surfaces are transformed into IPSs during urbanization, the areas and densities of various buildings and structures are increasing [20]. Therefore, the area, abundance, landscape pattern, and other spatial characteristics of IPS also vary [21]. With the development of remote sensing technology, accurate large-scale IPS extraction has been realized [22], and the spatial distribution of land surface temperature (LST) can also be retrieved using satellite thermal infrared data [23,24]. With the support of multisource remote sensing data, research on the relationship between the IPS and the UHI effect has become popular. A large number of studies have shown that the IPS has a significant correlation with the urban thermal environment, including linear and non-linear positive correlations, which signifies that the regional LST will rise with an increase in the IPS area [25]. Some case studies have also found that an increase in construction density is positively correlated with LST [26][27][28][29]. However, these studies mainly focused on the numerical relationship between IPS indices and LST [30], and rarely involved the perspective of landscape ecology.
To explore the spatial effect of IPS on the UHI effect, some studies have revealed the influence of the IPS landscape patterns on land surface temperature [21,27]. Landscape patterns generally refer to the spatial arrangement and combination of landscape elements with different sizes and shapes, including the type, number, spatial distribution, and configuration of landscape constituent units [31], which can be quantified by various landscape metrics based on appropriate mathematical algorithms. Landscape patterns have always been related to the UHI phenomenon, and each land cover type constituting the landscape has unique radiation, thermal, and moisture characteristics, thus affecting the regional thermal environment [9]. In recent decades, owing to the increase in available satellite data, the relationship between landscape patterns and the UHI effect has been widely discussed from the perspectives of composition and configuration [32][33][34][35]. Studies have discovered that landscape features such as the size, shape, and spatial distribution of IPS can significantly affect the scale and extent of UHI [21,31,36]. However, in most studies, all IPSs are analyzed as a single landscape category without detailed classification. Therefore, the impact of the spatial pattern of IPS patches with various densities on LST remains unclear.
In this study, Landsat 8 remote sensing data and meteorological data were selected to extract the IPS density (IPSD) and LST of the study area by pixel unmixing and thermal infrared temperature retrieval, respectively. The study area was then classified into five levels according to their IPSD values. Finally, the spatial effects of spatial patterns of various IPSD levels on the UHI effect were revealed using landscape pattern analysis and geospatial analytic approaches.

Study Area
The built-up area of Xuzhou City (Figure 1) was selected as the study area, which is located in the northwest of Jiangsu Province, China (between 116 • 22 -118 • 40 E and 33 • 43 -34 • 58 N). In this study, the built-up area refers to the highly urbanized areas dominated by artificial buildings and roads in the main urban area of Xuzhou (including Quanshan District, Gulou District, Yunlong District, Tongshan District, and Jiuli District). Approximately 90% of the area of Xuzhou is dominated by plains, with an average altitude of approximately 40 m. Xuzhou is also located in a typical warm temperate zone and semi monsoon region, with an annual average temperature of 14 • C and an annual average rainfall of 800-930 mm. By the end of 2021, the GDP of Xuzhou reached 127.32 billion USD, a permanent population of 9.0285 million, and an urbanization rate of 66.2%. From 1995 to 2019, Xuzhou experienced rapid urbanization, with the scale of the urban built-up area increasing from 59 km 2 to 274 km 2 [37]. Xuzhou is also an important transportation hub, coal-producing area, and power-base in eastern China. In addition, Xuzhou has a largescale manufacturing industry. Rapid urbanization and industrialization have intensified the expansion of the IPS area and enhanced the UHI effect in the built-up areas of Xuzhou City. As a typical medium-sized city in China, Xuzhou has a high level of urbanization, industrialization, and transportation. It is not only widely distributed by impervious surfaces, but also has a significant UHI effect [20]. Therefore, the study area is very suitable for exploring the influence of IPS spatial characteristics on the heat island effect.
Quanshan District, Gulou District, Yunlong District, Tongshan District, and Jiuli District). Approximately 90% of the area of Xuzhou is dominated by plains, with an average altitude of approximately 40 m. Xuzhou is also located in a typical warm temperate zone and semi monsoon region, with an annual average temperature of 14 °C and an annual average rainfall of 800-930 mm. By the end of 2021, the GDP of Xuzhou reached 127.32 billion USD, a permanent population of 9.0285 million, and an urbanization rate of 66.2%. From 1995 to 2019, Xuzhou experienced rapid urbanization, with the scale of the urban built-up area increasing from 59 km 2 to 274 km 2 [37]. Xuzhou is also an important transportation hub, coal-producing area, and power-base in eastern China. In addition, Xuzhou has a large-scale manufacturing industry. Rapid urbanization and industrialization have intensified the expansion of the IPS area and enhanced the UHI effect in the built-up areas of Xuzhou City. As a typical medium-sized city in China, Xuzhou has a high level of urbanization, industrialization, and transportation. It is not only widely distributed by impervious surfaces, but also has a significant UHI effect [20]. Therefore, the study area is very suitable for exploring the influence of IPS spatial characteristics on the heat island effect.

Data
Three scenes of Landsat 8 images from 2014 to 2020 were downloaded from the United States Geological Survey (USGS) data center (https://glovis.usgs.gov/, accessed on 1 September 2022). Their multispectral (bands 1-7) and thermal infrared (band 10) bands provide spatial source data for land surface component extraction and LST inversion, respectively. The time interval of the three images is 3 years. During these two periods, Xuzhou City was in the process of rapid urbanization, and the IPS has experienced obvious evolution, which is conducive to reflecting the effects of different IPS spatial patterns on LST. In addition, the three images are all in May (close to summer), which can not only show a significant UHI effect, but also avoid the impact of seasonal differences on the research results. Another satellite dataset used in pixel unmixing validation was the GF-1 high-resolution remote sensing image obtained in April 2020, whose spatial resolution reached 2 m with the fusion of panchromatic and multispectral bands. It should be noted that there was only a 19-day difference in the acquisition date of the GF-1 image and the corresponding Landsat 8 data in 2020. Therefore, the land cover change during this period can be ignored.
In addition, ground observations at the transit time of the Landsat 8 satellite, including near-surface temperature (Tair) and air relative humidity (RH), were also neces-

Data
Three scenes of Landsat 8 images from 2014 to 2020 were downloaded from the United States Geological Survey (USGS) data center (https://glovis.usgs.gov/, accessed on 1 September 2022). Their multispectral (bands 1-7) and thermal infrared (band 10) bands provide spatial source data for land surface component extraction and LST inversion, respectively. The time interval of the three images is 3 years. During these two periods, Xuzhou City was in the process of rapid urbanization, and the IPS has experienced obvious evolution, which is conducive to reflecting the effects of different IPS spatial patterns on LST. In addition, the three images are all in May (close to summer), which can not only show a significant UHI effect, but also avoid the impact of seasonal differences on the research results. Another satellite dataset used in pixel unmixing validation was the GF-1 high-resolution remote sensing image obtained in April 2020, whose spatial resolution reached 2 m with the fusion of panchromatic and multispectral bands. It should be noted that there was only a 19-day difference in the acquisition date of the GF-1 image and the corresponding Landsat 8 data in 2020. Therefore, the land cover change during this period can be ignored.
In addition, ground observations at the transit time of the Landsat 8 satellite, including near-surface temperature (T air ) and air relative humidity (RH), were also necessary for LST retrieval. The ground observations of May 1, 2014, and May 16, 2017, were collected from the Collaborative Observation Test Site of China University of Mining and Technology in the study area (data recording frequency of 30 min), and the meteorological data of May 17, 2020, were obtained from the NOAA National Centers for Environmental Information (NCEI) of USA (https://www.ngdc.noaa.gov/, accessed on 7 September 2022) (recording frequency of 1 h). All of the satellite data and ground observations are presented in Table 1.

. Endmember Fraction Extraction
Urban land cover can be expressed as a combination of three ecological elements (vegetation, impervious surface, and soil) at a certain ratio, using the vegetation-impervious surface-soil (V-I-S) model [36]. In medium-and low-resolution remote sensing images, a single pixel often contains more than two land cover components. The linear spectral mixture analysis (LSMA) model assumes that the reflectance of each mixed pixel in a certain band has a linear relationship with the reflectance of all components in the pixel, and the linear coefficient is the area proportion of each component [22]. Therefore, the least-squares method can be used to solve the component fractions of the mixed pixels. However, without any constraints, an endmember fraction extracted by the LSMA may be greater than 1 or less than 0, or the sum of all endmember fractions in a mixed pixel is not equal to 1. Then, the fully constrained linear spectral mixture analysis (FCLS) [38] was developed by considering the two constraints. On this basis, a normalized spectral mixture analysis model (NSMA) was proposed to further solve the negative impact of shadows on mixed pixel decomposition results [39]. Based on the diversity of impervious surfaces in the study area, the endmembers in a mixed pixel were divided into grass (including farmland), forest, soil, and high-and low-albedo impervious surfaces (water generally exists independently and can be masked). The LSMA algorithms are as follows.
where R b is the standardized reflectance of a mixed pixel for band b; R b is the original reflectance of the mixed pixel for band b; N is the endmember number of the mixed pixel; f i is the fraction of endmember i; R i,b is the standardized reflectance of endmember i for Sample selection of each endmember should be completed in ENVI 5.3 before LSMA, as shown in Figure 2. First, the eigenvalues of all bands were calculated using the minimum noise fraction rotation (MNF), which can be used to eliminate nonmajor bands. The pure pixels of all endmembers were extracted from the major bands of the MNF using the pure pixel index (PPI). Finally, the pure pixels were superposed onto the major bands of the MNF as a region of interest to generate a scatter plot in the ENVI n-D visualizer, and each clustered scatter cloud was marked as the sample of an endmember. After the endmember fractions were extracted using the NSMA model, the IPSD was calculated as the sum of the fractions of the high-and low-albedo IPS.

=1
where is the standardized reflectance of a mixed pixel for band b; Rb is the original reflectance of the mixed pixel for band b; N is the endmember number of the mixed pixel; fi is the fraction of endmember i; , is the standardized reflectance of endmember i for band b; eb is the model calculation residual for band b; and RMS is the fitting accuracy of the LSMA model.
Sample selection of each endmember should be completed in ENVI 5.3 before LSMA, as shown in Figure 2. First, the eigenvalues of all bands were calculated using the minimum noise fraction rotation (MNF), which can be used to eliminate nonmajor bands. The pure pixels of all endmembers were extracted from the major bands of the MNF using the pure pixel index (PPI). Finally, the pure pixels were superposed onto the major bands of the MNF as a region of interest to generate a scatter plot in the ENVI n-D visualizer, and each clustered scatter cloud was marked as the sample of an endmember. After the endmember fractions were extracted using the NSMA model, the IPSD was calculated as the sum of the fractions of the high-and low-albedo IPS.

Endmember Fraction Validation
The endmember fractions extracted through the LSMA should be validated using high-resolution maps of land cover types. As shown in Figure 1b, all the sample points were set at intervals of 900 m along the four directions. Each pixel containing a sample point was regarded as the center pixel to create a 90 × 90 m 2 sample area. The fractions of IPS (including high-and low-albedo) and vegetation (including grass and forest) were then extracted from each sample area. All sample areas were superimposed on the GF-1 image, and manual interpretation was used to extract the proportions of land cover types in each sample area, which can be regarded as the true values of the component fractions. Although we have only found one phase of high-resolution satellite images whose transit time was close to the acquisition date of the Landsat 8 data of 2020, the pixel unmixing of the three Landsat 8 images was completely consistent with the rules and methods of endmember extraction. Therefore, the validation result of 17 May 2020 can represent the accuracy of the LSMA in this study.

Endmember Fraction Validation
The endmember fractions extracted through the LSMA should be validated using high-resolution maps of land cover types. As shown in Figure 1b, all the sample points were set at intervals of 900 m along the four directions. Each pixel containing a sample point was regarded as the center pixel to create a 90 × 90 m 2 sample area. The fractions of IPS (including high-and low-albedo) and vegetation (including grass and forest) were then extracted from each sample area. All sample areas were superimposed on the GF-1 image, and manual interpretation was used to extract the proportions of land cover types in each sample area, which can be regarded as the true values of the component fractions. Although we have only found one phase of high-resolution satellite images whose transit time was close to the acquisition date of the Landsat 8 data of 2020, the pixel unmixing of the three Landsat 8 images was completely consistent with the rules and methods of endmember extraction. Therefore, the validation result of 17 May 2020 can represent the accuracy of the LSMA in this study.

Landscape Metrics Selection for Five ISPD Levels
According to the IPSD values, the study area was classified into five levels using the mean-standard deviation method, including extremely high (LV5), sub high (LV4), medium (LV3), sub low (LV2), and low (LV1). The five levels were regarded as patches with various IPSDs to explore the effects of their landscape patterns on the UHI effect.
The proportion, shape complexity, predominance, natural connectivity, and aggregation are the five typical factors that can represent the spatial characteristics of patches, and the landscape metric is an effective approach for quantifying these spatial features [27]. Therefore, five class-level metrics-percentage of landscape (PLAND), landscape shape index (LSI), largest patch index (LPI), patch cohesion index (COHESION), and aggregation index (AI)-were selected to represent these spatial characteristics of each IPSD level, respectively. In addition, four landscape-level metrics: Shannon's diversity index (SHDI), patch diversity (PD), Shannon's evenness index (SHEI), and competition (CONTAG), were Land 2022, 11, 2135 6 of 20 selected to represent the diversity, number, distribution, and aggregation of all IPSD levels in the same region. All of the selected landscape metrics are listed in Table 2. Table 2. The landscape metrics selected in this study.

Landscape Metrics Formulas
Class level

Scale Effect Analysis
In the assessment of landscape patterns and processes, most landscape metrics are highly dependent on the observation scale, particularly in remote sensing-based research, as the scale effect may affect the survey results of landscape structure [40]. The scale effect has two aspects: the grain-size effect and spatial-extent effect [41], which represent the impact of the variations of the size of the smallest spatial analysis unit and the analysis extent on landscape structure extraction, respectively. First, the best grain size should be analyzed. Each IPSD level map was resampled at 10 m interval from 30-200 m to produce 18 images with different spatial resolutions (grain-sizes). FRAGSTATS 4.2 was used to calculate the class-level metrics of these resized images, and the response trend of each metric to the change in grain-size was recorded, as shown in Figure 3a-e (taking the data of 17 May 2020 as an example). It can be seen that the inflection points of the response curves of most metrics occur at the position where the grain-size is 50 m. In addition, according to the land area information loss index (Figure 3f) calculated using Equations (5) and (6), the land area accuracy loss at a grain-size of 50 m is relatively low. Therefore, 50 m can be selected as the best grain-size for landscape pattern analysis of IPSD levels.
where L j is the relative value of the patch area loss of IPSD level j, A j is the patch area of IPSD level j at a converted grain size, A bj is the patch area of IPSD level j at the original grain size, n is the number of IPSD levels, and S j is the land area loss index. The moving-window method was applied to calculate landscape metrics [42,43]. The landscape metrics within the spatial-extent of the moving-window window were assigned to the center pixel of the window and then moved to the next extent, and the spatial distribution map of each landscape metric was produced. The spatial-extent effect can be analyzed by adjusting the moving window size. To avoid pixel segmentation, the tested moving window sizes (spatial-extent) should be set as multiples of the best grain-size (n × 50 m). Therefore, six moving-window sizes at 250 m intervals from 250 × 250 m 2 to 1500 × 1500 m 2 were tested. FRAGSTATS 4.2 was applied to calculate the spatial distribution maps of four landscape-level metrics for six tested moving-window sizes, respectively. The sample points shown in Figure 1b were used to extract these metrics, as shown in Figure 4. The fluctuation curves of the four landscape-level metrics began to stabilize when the moving-window size increased to 750 × 750 m 2 , which indicates that 750 × 750 m 2 can be set as the best spatial-extent for analyzing the landscape patterns of the IPSD levels.
ing Equations (5) and (6), the land area accuracy loss at a grain-size of 50 m is relatively low. Therefore, 50 m can be selected as the best grain-size for landscape pattern analysis of IPSD levels.
where Lj is the relative value of the patch area loss of IPSD level j, Aj is the patch area of IPSD level j at a converted grain size, Abj is the patch area of IPSD level j at the original grain size, n is the number of IPSD levels, and Sj is the land area loss index. The moving-window method was applied to calculate landscape metrics [42,43]. The landscape metrics within the spatial-extent of the moving-window window were assigned to the center pixel of the window and then moved to the next extent, and the spatial distribution map of each landscape metric was produced. The spatial-extent effect can be analyzed by adjusting the moving window size. To avoid pixel segmentation, the tested moving window sizes (spatial-extent) should be set as multiples of the best grain-size (n×50 m). Therefore, six moving-window sizes at 250 m intervals from 250×250 m 2 to 1500×1500 m 2 were tested. FRAGSTATS 4.2 was applied to calculate the spatial distribution maps of four landscape-level metrics for six tested moving-window sizes, respectively. The sample points shown in Figure 1b were used to extract these metrics, as shown in Figure 4. The fluctuation curves of the four landscape-level metrics began to stabilize when the moving-window size increased to 750 × 750 m 2 , which indicates that 750 × 750 m 2 can be set as the best spatial-extent for analyzing the landscape patterns of the IPSD levels.

Landscape Surface Temperature Retrieval
Landsat 8 has two thermal infrared bands, 10 and 11. Owing to the uncertainty of information in band 11 [23,24], band 10 was selected as the thermal radiation data. The

Landscape Surface Temperature Retrieval
Landsat 8 has two thermal infrared bands, 10 and 11. Owing to the uncertainty of information in band 11 [23,24], band 10 was selected as the thermal radiation data. The mono-window algorithm [44,45] based on the thermal radiation transfer equation was used for LST inversion. The algorithm requires three parameters: land surface emissivity(ε), atmospheric transmittance (τ), and effective average atmospheric temperature (T air_e ). ε is determined by the thermal physical characteristics of the land cover types and is primarily affected by the atmospheric water vapor content (w), and T air_e is related to atmospheric temperature and water vapor distribution. Previous research reported that the inversion accuracy of the mono-window algorithm was ±1.4 • C [45], and the equations were as follows: T 10 = 1321.08/ln(1 + 774.89/L 10 ) where a = −62.7182 and b = 0.4339 are the regression coefficients of the Planck blackbody radiation function within 0-70 • C; T 10 is the brightness temperature calculated by thermal infrared band 10; L 10 is the thermal radiation parameter calculated by the DN value of band 10; P v is the proportion of vegetation area to pixel, which can be calculated by the Normalized Difference Vegetation Index (NDVI); R v and R x (R m or R s ) represent the radiation ratios of pure vegetation pixels, building pixels, and soil pixels, respectively; and d ε is the interference value of the interaction between vegetation and soil on the land surface emissivity. According to the simulation results of the Modern Resolution Atmospheric Transmission (MODTRAN 4) program, the linear relationships between τ and w for Landsat 8 TIR band 10 are shown in Table 3 [45]. Although the atmospheric water content data of the study area are not directly obtained, they can be calculated according to the empirical formula fitted from the historical data of the China Meteorological Station. The necessary parameters included the average altitude (H E = 40 m), the latitude (ϕ = 34.24 • ), and the relative air humidity (RH) of the study area, as shown in Equations (12) and (13) [46]. In the mid-latitude summer, the linear relationship between T air_e and T air is shown in Equation (14) [44].

Bivariate Moran's I Analysis
Two types of bivariate Moran's I [47] were applied to explore the spatial correlation between landscape metrics of IPSD levels and LST, namely bivariate global Moran's I and bivariate local Moran's I. Bivariate global Moran's I indicated whether there is spatial correlation between landscape metrics and LST in the whole region and the degree of their spatial correlation, while bivariate local Moran's I shows the spatial correlations in different spatial units. The equations are as follows: where I and I i are the bivariate global and local spatial autocorrelation indices, respectively; n is the total number of spatial units; W ij is the spatial weight matrix established by K adjacency method; x i and y j are the observed values of independent variable and dependent variable in spatial unit i and j respectively; S 2 is the variance of all samples; z i and z j are the normalized values of the variances of the observations in spatial units i and j. The values of I and I i range from −1 to 1. A positive value indicates a positive spatial correlation between landscape metrics and LST, meaning that pixels with high landscape metric values may be surrounded by pixels with high LST values. By contrast, a negative value represents a negative spatial correlation, indicating that pixels with high landscape metric values may be surrounded by pixels with low LST values. The larger the absolute values of I and I i , the stronger the spatial correlation between landscape metrics and LST. The four quadrants generated based on the local indicators of spatial association (LISA) of I i represent four types of local spatial autocorrelation [48]. Quadrant I, called High-High (HH) type, represents the high landscape metrics value surrounded by high LST value; Quadrant II is High-Low (HL), indicating the high landscape metrics value surrounded by low LST value; Quadrant III is Low-Hight (LH), indicating a low landscape metrics value surrounded by a high LST value; and Quadrant IV is Low-Low (LL), indicating the low landscape metrics value surrounded by low LST values.

Inversion Results of IPSD and LST
The IPSD maps calculated by LSMA are shown in Figure 5a-c. Compared with the high-resolution land cover data extracted from the GF-1 image (Figure 6), the inversed fractions (taking the components of IPS and vegetation on 17 May 2020, as an example) had a significant and strong linear correlation with the true fractions (with r > 0.90 and R 2 > 0.82, p < 0.001), which signified that the accuracy of the inversed IPSD was reliable. The IPS area of the study area has gradually expanded, and high-density impervious patches have significantly increased, which is consistent with the rapid urbanization of Xuzhou City during this period. The LST images retrieved by the mono-window algorithm are shown in Figure 5d-f. From 2014 to 2020, the study area also showed an apparent UHI effect, and the LST of the areas covered by IPS was higher than that of the surrounding vegetation, soil, and other natural surfaces. IPSD and LST values were extracted by the sample points shown in Figure 1b. The quantitative relationship between IPSD and LST is simulated using linear, quantitative, cubic and empirical regression models, respectively. The results are shown in Table 4. It can be seen that the fitting results of the four regression models are significant at p < 0.001, and the cubic model has the highest fitting goodness (with 0.447 < R 2 < 0.535) for three phases. Figure 7a-c reveals that the LST of three phases generally shows an upward trend with the increase of IPSD, with Pearson's correlation coefficients of 0.63, 0.69, and 0.66, respectively. In addition, as shown in Figure 7d, the average LST of each IPSD level indicated that the average LST also increased with an increase in the IPSD level, which further proved the enhancement effect of IPSD on the thermal environment. Therefore, it is valuable to further explore the spatial relationship between landscape patterns of IPSD levels and LST.

Spatial Correlation between Landscape Metrics of IPSD Levels and LST
The bivariate global Moran's I values between landscape metrics of the five IPSD levels and LST are listed in Table 5 Significant spatial correlations were observed between all the selected landscape metrics and LST. However, these spatial correlations show great differences with variations in landscape metrics and IPSD levels. For landscape-level metrics, SHDI, SHEI, and PD showed significant positive spatial correlations with LST (p < 0.001), indicating that increasing the number of IPSD levels, patch even distribution degree, and patch density could enhance the regional thermal environment. In contrast, CONTAG showed a significant negative spatial correlation with LST (p < 0.001), which signified that increasing the aggregation degree of the five IPSD levels would alleviate the UHI effect.
For class-level metrics which represent the independent spatial characteristics of the patches of each IPSD level for an area, PLAND, LPI, LSI, AI, and COHESION showed significant correlations (p < 0.001) with LST, indicating that the spatial characteristics of each IPSD level, such as patch proportion, predominance degree, shape complexity, ag-

Spatial Correlation between Landscape Metrics of IPSD Levels and LST
The bivariate global Moran's I values between landscape metrics of the five IPSD levels and LST are listed in Table 5 Significant spatial correlations were observed between all the selected landscape metrics and LST. However, these spatial correlations show great differences with variations in landscape metrics and IPSD levels. For landscape-level metrics, SHDI, SHEI, and PD showed significant positive spatial correlations with LST (p < 0.001), indicating that increasing the number of IPSD levels, patch even distribution degree, and patch density could enhance the regional thermal environment. In contrast, CONTAG showed a significant negative spatial correlation with LST (p < 0.001), which signified that increasing the aggregation degree of the five IPSD levels would alleviate the UHI effect. For class-level metrics which represent the independent spatial characteristics of the patches of each IPSD level for an area, PLAND, LPI, LSI, AI, and COHESION showed significant correlations (p < 0.001) with LST, indicating that the spatial characteristics of each IPSD level, such as patch proportion, predominance degree, shape complexity, aggregation degree, and natural connectivity degree, have significant spatial impacts on the thermal environment. However, the direction and degree of these spatial impacts varied at different IPSD levels. First, the bivariate global Moran's I between LST and class-level metrics of IPSD LV1 or LV2 were negative. Second, the bivariate global Moran's I between LST and class-level metrics of IPSD LV3 to LV5 were positive, and the value of Moran's I increased gradually as the level of IPSD increased.
In addition, LST have an obviously stronger spatial correlation with the class-level metrics of IPSD LV4 and LV5 (0.2 < Moran's I < 0.6) than IPSD LV3 (most Moran's I < 0.1). Therefore, the class-level metrics of IPSD LV4 and LV5 are key factors that affect the spatial characteristics of the thermal environment and deserve special attention. Based on the degree of spatial correlation between the LST and class-level metrics of IPSD LV4 or LV5, these metrics can be sorted as PLAND (average Moran's I = 0.57), LSI (average Moran's I = 0.56), COHESION (average Moran's I = 0.50), LPI (average Moran's I = 0.46), and AI (average Moran's I = 0.30). However, the spatial correlation degree cannot fully represent the contribution ratio of each class-level metric of IPSD LV4 and LV5 to LST; therefore, further discussion on the impact weights of these metrics is necessary.
The bivariate local Moran's I between the class-level metrics of IPSD LV4 or LV5 and LST are represented as LISA maps (Figure 8), which show the spatial distributions of the correlation between the two variables. The spatial clustering mode between the class-level metrics (PLAND, LSI, COHESION, LPI, and AI) and LST were mostly HH and LL. The superposition results of LISA maps and land cover classification maps (Figure A1 in Appendix A) showed that most of the HH areas were distributed in urban built-up areas with IPS coverage, whereas LL areas were mainly distributed in urban forests and parks with vegetation coverage, as well as farmland in the suburbs.
14 of 21 Figure 8. LISA maps between landscape metrics of IPSD LV4 or LV5 and LST (HH refers to the spatial aggregation of high landscape metrics and high LST; LL refers to the spatial aggregation of low landscape metrics and low LST; HL refers to the spatial aggregation of high landscape metrics and low LST; and LH refers to the spatial aggregation of low landscape metrics and high LST).

The Impact Weights of Class-Level Metrics of IPSD LV4 and LV5 on LST
The spatial correlation analysis in the previous sections proved that class-level metrics of IPSD LV4 and LV5 have significant negative spatial effects on LST. However, the contribution of each metric to the spatial impact on LST still needs to be discussed further. In general, the standardized coefficients of multiple linear regression (MRA) can be

The Impact Weights of Class-Level Metrics of IPSD LV4 and LV5 on LST
The spatial correlation analysis in the previous sections proved that class-level metrics of IPSD LV4 and LV5 have significant negative spatial effects on LST. However, the contri- bution of each metric to the spatial impact on LST still needs to be discussed further. In general, the standardized coefficients of multiple linear regression (MRA) can be used to express the contribution weight of each independent variable (class-level metrics of IPSD LV4 and LV5) to the dependent variable (LST). However, the multicollinearity among these metrics was revealed during testing, which prevented the direct use of MRA. The principal component analysis (PCA) can recombine multiple original highly correlated independent variables into a group of unrelated comprehensive variables, which are called major component variables (F i ) [49]. As collinearity was eliminated among the principal component variables, MRA could be directly used to extract standardized regression coefficients (β i ) between F i and LST. Finally, the impact weights of each class-level metric of IPSD LV4 and LV5 on LST were calculated using the eigenvalue (η i ) of the major component variables, the PCA component matrix, and the standardized regression coefficients.
All independent and dependent variables were standardized before PCA. The total variance explained by PCA (Table 6) showed that the Kaiser-Meyer-Olkin measure of sampling adequacy (KMO) values were all above 0.73, indicating that the PCA method was applicable to the independent variables of the class-level metrics of IPSD LV4 and LV5. Two principal components (F 1 and F 2 ) were extracted from the independent variables of the three phases, and all cumulative proportions of variance were more than 78%, indicating that F 1 and F 2 could effectively represent all independent variables. The component matrix of PCA represents the correlation coefficient (θ 1 and θ 2 ) of all the original independent variables (V1 to V10) with each principal component variable (F 1 and F 2 ), which are shown in Table 7 (where θ i < 0.2 was considered to be too low and was rejected). As shown in Equations (18) and (19) [49], F 1 and F 2 can be calculated as the weighted sum of the original variables and their contribution coefficient α ij (calculated using Equation (17)). Table 8 presents the multivariate linear regression results between the dependent variable normalized LST and independent variables F 1 and F 2 . The significant and strong linear relationship between independent variables and dependent variables for each phase is evident ((with p < 0.001 and r > 0.59). Here, the impact weights (β i ) of F 1 and F 2 on the LST were obtained. Then, by combining β i with the contribution coefficient α ij of the original variables to F 1 and F 2 , the impact weights of the class-level metrics of IPSD LV4 and LV5 on LST were calculated using Equation (20). Figure 9 shows that the impact weights of the class-level metrics of IPSD LV5 are higher than those of LV4. For the LV5 patches, the impact weights of the metrics were ranked from high to low as PLAND, COHESION, LSI, LPI, and AI. In addition, the impact weights of the first four metrics were similar, indicating that the spatial characteristics of IPSD LV5, such as area proportion, natural connectivity, shape complexity, and predominance degree have similar regulatory effects on the thermal environment. For the LV4 patches, the impact weights of the class-level metrics on LST were ranked differently from those of LV5, that is, LSI, PLAND, COHESION, LPI, and AI from high to low.  the original variables to F1 and F2, the impact weights of the class-level metrics of IPSD LV4 and LV5 on LST were calculated using Equation (20). Figure 9 shows that the impact weights of the class-level metrics of IPSD LV5 are higher than those of LV4. For the LV5 patches, the impact weights of the metrics were ranked from high to low as PLAND, COHESION, LSI, LPI, and AI. In addition, the impact weights of the first four metrics were similar, indicating that the spatial characteristics of IPSD LV5, such as area proportion, natural connectivity, shape complexity, and predominance degree have similar regulatory effects on the thermal environment. For the LV4 patches, the impact weights of the class-level metrics on LST were ranked differently from those of LV5, that is, LSI, PLAND, COHESION, LPI, and AI from high to low.

Discussion
In recent studies on the relationship between the IPS spatial pattern and UHI, IPS is generally regarded as a single land cover type (such as a building or built-up area) in the calculation of landscape pattern metrics [8,49], without considering the impact caused

Discussion
In recent studies on the relationship between the IPS spatial pattern and UHI, IPS is generally regarded as a single land cover type (such as a building or built-up area) in the calculation of landscape pattern metrics [8,49], without considering the impact caused by the difference of IPS density in the patches. Since the positive correlation between IPS and UHI has been proved by many studies [21,27], our study explored an innovative method to incorporate IPS density into landscape pattern analysis. According to the component fraction value of IPS, the land surface of the study area is classified into five levels (or can be regarded as five patch types), which corresponds to the land cover types in general landscape pattern analysis. Although the density index used for simple classification cannot completely simulate the urban IPS with high heterogeneous, it is still helpful to reveal the impact of the spatial pattern of the patches with various IPS densities on UHI.
It should be noted that landscape-level metrics represent the comprehensive spatial characteristics of the patches of all spatial elements in an area. The patches of the five IPSD levels have different ratios of natural surface and impervious surface (for example, level 1 contains higher natural surface and lower IPS ratios, while level 5 has the opposite ratios). Therefore, landscape-level metrics can be regarded as comprehensive spatial pattern indicators of vegetation and IPS patches with various proportions. This further indicates that the relationship between landscape-level metrics and LST results from the interaction between the IPS and natural patches in the region, including the enhancement effect of high IPS density patches (LV4 and LV5) on thermal environment and the mitigation effect of low IPS density patches (LV1 and LV2 with higher vegetation or soil density) on UHI.
In previous studies, IPS is generally regarded as one of the land cover types, and its class-level metrics PLAND, LPI, LSI, COHESION, and AI have been proved to have significant positive correlations with LST [50,51]. This is consistent with our results that these class-level metrics of IPSD LV3 to LV5 (with higher IPS proportion) show positive bivariate global Moran's I with LST (p < 0.001).
From IPSD LV1 to LV2, the spatial patterns of the patches are negatively correlated with LST, which because these patches are mainly covered by natural surfaces (vegetation and soil) with low IPS proportions, resulting in that the mitigation effect of vegetation on the thermal environment is dominant. On the contrary, from IPSD LV3 to LV5, the spatial patterns of the patches began to be positively correlated with LST, which owing to the increase in IPSD is accompanied by the enhancement of surface thermal radiation and anthropogenic heat emissions, and the heating effect of IPS thermal completely offsets the cooling effect of natural surfaces. With the gradual increase of IPS density, Moran's I tends to increase. In addition, the land cover type in IPSD LV5 patches are mostly IPS. However, IPSD LV4 patches are still mixed with more natural surfaces. This is an important reason for the difference in the impact weight rankings between the class-level metrics of LV4 and LV5. Previous research found that the landscape metric PLAND of building or built-up patches has the strongest correlation with LST [8], and the correlation between LPI and LST is stronger than that between AI and LST [51], which are also the same as our results. According to our research, priority should be given to reducing the area proportion, natural connectivity degree and shape complexity degree of high-density IPS regions, which could be one of the effective approaches to alleviate the regional thermal environment.

Conclusions
Based on three Landsat 8 images of the study area from 2014 to 2020, the LSMA and mono-window algorithms were used to extract the IPSD and LST, respectively. Landscape and geospatial analyses were applied to explore the spatial correlation between the landscape metrics of various IPSD levels and LST. It was revealed that the four landscape-level metrics, SHDI, SHEI, and PD all showed significant positive spatial correlations with LST, while CONTAG showed a significant negative correlation, indicating that increasing the regional number of IPSD levels, patch even distribution degree, and patch density, or reducing the aggregation degree between the five IPSD levels will lead to the enhancement of the thermal environment. It should be noted that the relationships between landscape-level metrics and LST reflect the comprehensive spatial effects of all IPSD level patches, which means that natural patches (mainly covered vegetation and soil) also contribute to these spatial effects on LST. For the five class-level metrics of each IPSD level, PLAND, LPI, LSI, AI, and COHESION, all showed significant correlations with LST. This means that the spatial characteristics of each IPSD level, such as patch proportion, predominance degree, shape complexity, aggregation degree, and natural connectivity degree, were all regulatory factors for the UHI effect. In particular, the degree of positive spatial correlation of the class-level metrics of IPSD LV4 and LV5 with LST were the highest. The impact weight of each class-level metric of IPSD LV4 and LV5 on LST was also obtained using PCA and MRA standardization coefficients, which can be sorted as PLAND_LV5 > COHESION_LV5 > LSI_LV5 > LPI_LV5 > AI_LV5 > LSI_LV4 > PLAND_LV4 > COHESION_LV4 > LPI_LV4 > AI_LV4. This study has revealed the importance of and differences in the spatial patterns of patches with various IPSDs in UHI effect regulation. It is suggested that in urban planning and environmental management, approaches to reduce the area proportion, natural connectivity degree, and shape complexity of regional high-density IPS are preferred to alleviate the UHI effect.
Considering the size, geographical location and industrial structure of the study area in this paper, our conclusions may be limited with regard to the applicability to medium-sized cities dominated by industry in central China. For other city types, such as agriculturedominated cities, commerce-dominated cities, and super large cities, the impact of IPS spatial pattern on the UHI effect needs to be discussed further. In addition, owing to the fact that the transit time of Landsat 8 images are basically fixed, the inversion time of UHI effect in this paper is limited to noon (around 11:00 AM, Beijing time), which cannot further reveal whether the impact of the IPS's spatial pattern on UHI varies with time. Another issue worth detailed study is the use of unmanned aerial vehicles (UAV) equipped with multispectral and thermal infrared sensors to obtain urban LST data at different times in one day, and to explore the temporal characteristics of the impact of IPS spatial patterns on the UHI effect.