Impacts of Spatial Configuration of Land Surface Features on Land Surface Temperature across Urban Agglomerations, China

Abstract

Booming urbanization triggers a significant modification of surface landscape configuration and hence complex urban climates. Considerable concerns exist regarding impacts of impervious surface area (ISA) and/or urban green space (UGS) on land surface temperature (LST). However, a knowledge gap still exists concerning the influence of urban landscape components and related spatial configuration on LST. To date, case studies have usually focused on individual cities, while few reports have addressed the impacts of urban surface components and relevant spatial configurations on LST within cities of different sizes, at different latitudes, and with different climatic backgrounds. Considering case studies from different latitudes and various climatic backgrounds can assist in obtaining comprehensive viewpoints about impacts of urban surface features on LST in both space and time. In this paper we analyzed data from three urban agglomerations, Beijing–Tianjin–Hebei (BTH), the Yangtze River Delta (YRD) and the Pearl River Delta (PRD), over the period 2000–2015. These three regions are densely populated with the most developed socio-economy across China, and are also dominated by booming urbanization. Based on Landsat remotely sensed data, we included the spatial pattern of surface components and related configuration into our analysis, quantifying impacts of spatial configuration of surface components on LST in both space and time. We found generally rising LST over all cities, which can be attributed to continuous urban expansion-induced decreased UGS. Generally, LST over ISA was 0.96–7.96 °C higher than that over UGS. We investigated the impacts of spatial pattern of land surface components on LST and found that the joint effect of the composition and spatial configuration of land surface components had the most significant impact on LST. Specifically, ISA and UGS had higher impact on LST than the impact of geometry of the ISA and UGS on LST. In the future, continuous expansion of ISA and continuous shrinking of UGS will drive the rising tendency of LST. Moreover, a larger rising tendency of LST will be observed in larger sized cities than smaller sized cities.

1. Introduction

Recent decades have witnessed a continual movement of the global population to cities [1]. According to statistics, 54% of the global population live in cities, and the number of megacities increased from 10 in 1990 to 28 in 2014, with a population of 453 million, accounting for 12% of the world’s total population [2]. Rapid urbanization results in the continual transformation of natural landscapes such as vegetation and water bodies, to impervious surface area and other artificial surfaces [3,4], which has modified the thermal exchange between the land surface and the atmosphere in urbanized areas and the thermodynamic properties of the urbanized underlying surface [5], resulting in higher air and land surface temperatures in the urban area, than in the surrounding rural areas. This phenomenon is also known as the urban heat island (UHI) effect [6,7]. UHI is usually divided into surface UHI (SUHI) measured by the LST, and atmospheric UHI measured by air surface temperature [7]. In this study, we mainly examine the relationship between land surface components and land surface temperature, therefore we focus on SUHI. With the development of satellite thermal infrared technology, remote sensing images such as Landsat, MODIS, etc., have been widely used to analyze the radiation energy from land surface components (LSC) such as vegetation, water bodies, unreclaimed wasteland and other natural surfaces, and roads, building roofs and other man-made surfaces [4]. Therefore, direct linkages between LST and the spatial pattern of land cover components exist [8]. Investigation of spatiotemporal connections between LST and LSC can improve understanding of the impacts of LSC or land surface properties on the UHI and can help to provide theoretical grounds for urban planning [9,10,11].

More and more research has appeared addressing the relationships between LST, impervious surface area (ISA) and urban green space (UGS) (e.g., [9,12,13,14]). The spatial pattern of ISA and UGS includes land surface compositions and related spatial configuration [15]. Previous studies have mainly focused on the relationship between specific LSC and related proportions and LST [12,13,16], such as positive relation between ISA, built-up areas and LST [17,18,19], and negative relation between LST and natural land surfaces such as forests, wetlands and other water bodies [13,20,21]. Landscape configuration can delineate the shape, size, location, and other spatial characteristics of a specific LSC, which is significant for quantifying the impacts of the spatial pattern of the LSC on LST [4].

However, we have limited knowledge about the coupled relation between LST and components, and related spatial configurations of the land surface properties [1,4]. Therefore, we cannot obtain a complete understanding of the spatial heterogeneity of the UHI and relevant impacts from LSC on LST. To fill this knowledge gap, our study combined the LSC of the ISA and UGS, and the spatial compositions and landscape configuration of the LSC, quantifying impacts of the spatial pattern of different LSC on the spatiotemporal heterogeneity of the LST over different urban sizes and levels.

Furthermore, in the backdrop of fast urban expansion, prediction of future urban expansion and land cover changes will help to optimize the allocation of land resources and scientific urban planning of the city, boosting inter-regional coordination and socioeconomic sustainable development [22]. In addition, the LST simulation for the future period, based on the spatial configuration of the LSC components, is significant for further understanding of future changes of the UHI and for mitigation of the UHI [23], providing theoretical grounds for scientific urban planning. Previous studies mainly focused on the relationship between LSC and LST for past periods, or the influences of land surface components on UHI. Few studies are available addressing future evolution of LSC and LST changes, and the relation between these two variables. Techniques for modeling urban expansion and land use simulation include cellular automata (CA) model [24], CLUE-S series model [25], InVEST model [26], CA-Markov model [27], etc. However, the use of only one specific modeling method cannot comprehensively describe the changes and transformations of different LSC. Liu et al. [28] proposed that the FLUS (future land-use simulation) model can integrate cellular automata (CA), and the neural network (ANN) algorithm model can effectively handle the conversion of different LSC, improving the accuracy of the modeling results. Regarding the prediction of the LST, the methods available are Markov chain, cellular automata (CA) and logistic regression [29]. Meanwhile, prediction of the LST can also be performed using a linear model developed by the LSC and the LST [30]. Due to the complexity of the land surface structure [31] and the seasonal variation of the LSC [13], the relationship between LSC and LST should be nonlinear. As such, we adopted the ANN algorithm in the simulation of future LST.

From 2002 to 2011, China’s urbanization rate has grown at an average annual rate of 1.35 percentage points, and the urban population has increased by an average of 20.96 million people per year. The study period selected, 2000–2015, is a typical representative period of the development stage of urbanization. In addition, there were many cities involved in this study, increasing the difficulty of collecting remote sensing images that met the study conditions. The development of urban agglomerations is the most remarkable feature of urbanization in China [32]. The literature shows studies relating the urban heat island effect and the urban component on urban agglomeration in China. Liu et al. [33] focused on 1288 urban clusters across China to study the relationship between SUHI and urban size, and their results showed that SUHI could effectively be mitigated by moderately dispersed, polycentric, and decentralized urban size. Su et al. [34] examined the effect of urban form on LST across varies cities, and the results indicated that compactness and urban shape significantly impacted LST in small sized, medium sized cities, and megacities. Our article selects three major urban agglomerations as the study area (Figure 1), namely, Beijing–Tianjin–Hebei (BTH), the Yangtze River Delta (YRD), and the Pearl River Delta (PRD), which, together, contribute 40% of China’s GDP, and 22.7% of China’s total population [35]. These three major urban agglomerations are densely populated, economically developed, and the UHI effect is prominent [4,36], having far-reaching impacts on the ecological environment in their vicinity. In the analysis of land surface temperature on specific cities, different authors used various indicators such as the land use land cover. He et al. [37] divided seven types of land use to study the impacts of environmental temperatures on the relationship between LST and land cover, and the results showed that the relation varied greatly with background temperature. Therefore, in a different climate background, understanding the impacts of the spatial pattern of the LSC of cities of different sizes and levels within these three major urban agglomerations on the temporal and spatial heterogeneity of the LST, can help provide scientific and theoretical grounds for mitigation of the UHI effect and for improvement of the urban ecological environment.

The principle objectives of this study are: (1) clarify the main LSC that affect the changes in the LST over cities of different sizes and levels, from 2000 to 2015, and characterize the temporal variations of the LSC; (2) identify impacts of spatial configurations of the LSC on LST, and quantify to what degree the components and spatial configuration of the land surface features on LST with integrated consideration of the impacts of the LSC and related spatial pattern on LST, changes; and (3) characterize the future changes of LST and LSC over cities of different sizes and levels. The structure of this study is as follows. Section 2 presents further details of the study region and the data. Section 3 details the methods used in this study. Section 4 presents the results and discussion, and Section 5 displays the major findings and conclusions of this current study.

2. The Study Region and the Data

2.1. The Study Region

The three major urban agglomerations, BTH, YRD and PRD (Figure 1), are highly populated and highly socio-economically developed regions of China. They are located along the coastal regions, from north China to southeast China, with different climate types and different ecological environments. Specifically, BTH is dominated by the temperate monsoon climate; YRD is characterized by the subtropical monsoon climate; and PRD is featured by the maritime monsoon climate [38]. Due to the rapid development of these three urban agglomerations, there are few cities with a resident population of less than one million. Therefore, based on the 2019 resident population data and the national classification standard, the cities Beijing, Tianjin, and Langfang within the BTH urban agglomeration were chosen as the supercity (with a resident population of more than 10 million), megacity (with a resident population of less than 10 million but over 5 million) and large city (with a resident population of more than one million but less than five million), respectively (Figure 1). The supercity, megacity and large city chosen within YRD are respectively Shanghai, Ningbo, and Nanjing (Figure 1). The supercity, megacity and large city chosen within PRD are respectively Guangzhou, Dongguan and Zhongshan (Figure 1). The above-mentioned cities are taken as case studies in this current study.

2.2. The Data

The remote sensing images used in this study included Landsat 5 TM and Landsat 8 OLI/TIRS (

Table A1. Landsat 5 TM and Landsat 8 OLI/TIRS images information of the study area.

Landsat Satellites	Landsat Data Identification	Data Acquisition Time
2000 Landsat 5 TM	LT51230322000145BJC00	24 May 2000
	LT51230332000145BJC00	24 May 2000
	LT51220322000170BJC00	18 June 2000
	LT51220332000170BJC00	18 June 2000
	LT51180382000158BJC02	06 June 2000
	LT51180392000158BJC02	06 June 2000
	LT51180402001080BJC00	21 March 2001
	LT51200372000140BJC00	19 May 2000
	LT51200382000204BJC00	22 July 2000
	LT51220432001060BJC00	01 March 2001
	LT51220442001060BJC00	01 March 2001
	LT51220452001060BJC00	01 March 2001
2007 Landsat 5 TM	LT51230322007148IKR00	28 May 2007
	LT51230323007148IKR00	28 May 2007
	LT51220322009258IKR00	15 September
	LT51220332009242IKR00	30 August 2009
	LT51180382007209BJC00	28 July 2007
	LT51180392007209BJC00	28 July 2007
	LT51180402007209BJC00	28 July 2007
	LT51200372006140BJC01	20 May 2006
	LT51200382007207IKR00	26 July 2007
	LT51220432008208BKT00	26 July 2008
	LT51220442008208BKT00	26 July 2008
	LT51220452009290BJC00	17 October 2009
2015 Landsat 8 OLI/TIRS	LC81230322015106LGN00	16 April 2015
	LC81230332015106LGN01	16 April 2015
	LC81220322015227LGN01	15 August 2015
	LC81220332015275LGN00	02 October 2015
	LC81180382015215 LGN00	03 August 2015
	LC81180392015215LGN00	03 August 2015
	LC81180402015215LGN00	03 August 2015
	LC81200372016344LGN01	09 December 2016
	LC81200382016344LGN01	09 December 2016
	LC81220432015291LGN00	18 October 2015
	LC81220442015291LGN01	18 October 2015
	LC81220452015003LGN00	03 January 2015

). These images covered the study regions effectively (http://earthexplorer.usgs.gov/ accessed on 9 January 2021). The spatial resolution of these remote sensing images was 30 m × 30 m. Images for each city in this study were selected from the same season, mainly the summer season, for the comparison of the spatial pattern of LST over time. In addition, this allowed the analysis of the remote sensing data to be less influenced by cloud coverage and free of data quality problems due to seasonal differences in vegetation growth. The above-mentioned remote sensing image data were used to extract the land surface components of each city and the LST changes.

3. Methods

3.1. Land Surface Temperature (LST) Retrieval

We used the radiative transfer equation method to retrieve LST during 2000–2015 for each city considered in this study, based on the thermal infrared band of Landsat 5 TM images (band 6, wavelength of 10.40–12.50 μm, with a resolution of 120 m) and Landsat 8 OLI/TIRS images (band 10, wavelength of 10.60–11.19 μm; band 11, wavelength of 11.5–12.51 μm, with a resolution of 100 m). Before we retrieved LST, we performed preprocessing of the Landsat remote sensing images, including radiometric calibration (both multispectral bands and thermal infrared bands), and FLAASH atmospheric correction based on ENVI software [39,40,41].

The retrieval procedure was as follows. First, we converted the DN value of the thermal infrared bands into top-of-atmospheric radiance. Then, based on the top-of-atmospheric radiance, we calculated the surface-leaving radiance [12,23]. In the third step, assuming earth is a black body, we converted the surface-leaving radiance value to at-satellite brightness temperature [4,41].

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

(1)

where ${\mathrm{L}}_{\lambda }$ is the top-of-atmospheric radiance, ${\mathrm{L}}_{\mathrm{u}}$ is the upwelling radiance values, L_d is the downwelling radiance values, τ is the atmospheric transmission. Among them, τ, ${\mathrm{L}}_{\mathrm{u}}$ and L_d were assessed using the Atmospheric Correction Parameter Calculator online tool (http://atmcorr.gsfc.nasa.gov). The tool uses the National Centers for Environment Prediction (NCEP) to model the atmospheric global profile at the specific date, time, and location as the input data [42]. The method was easy to obtain synchronously and simplified the retrieval process of the algorithm. In addition, when using synchronous atmospheric parameter data, the accuracy RMSE of radiative transfer equation method was 0.6 K, which showed high accuracy [43]. Before this, we carried out atmospheric correction on the remote sensing images, which greatly eliminated the influence of atmospheric molecular factors such as water vapor, carbon dioxide, oxygen and aerosol scattering on the reflection of land surface objects in the atmosphere, and by obtaining land surface reflectance information accurately, improved the accuracy of land surface temperature retrieval.

ε is the surface emissivity [44,45]. According to Qin et al. [46], from the perspective of the scale of satellite pixels, remote sensing images can be roughly regarded as composed of three types: water surface, building (including urban and rural, mainly consisting of roads and various buildings) and natural surface (mainly referring to natural land surface, woodland and farmland, etc.). Therefore, it was necessary to estimate the surface emissivity according to the proportion of main types of land surface, and the surface emissivity ε was divided into three categories: water, building and natural surface (Equation (2)).

$${\left\{\begin{array}{c}{\mathsf{\epsilon }}_{water}=0.995 \left(\mathrm{NDVI}\le 0\right)\\ {\mathsf{\epsilon }}_{building}=0.9589+0.086\times {\mathrm{P}}_{\mathrm{V}}-0.0671\times {\mathrm{P}}_{\mathrm{V}}{}^2 (0<NDVI<0.7)\\ {\mathsf{\epsilon }}_{natural}=0.9625+0.0614\times {\mathrm{P}}_{\mathrm{V}}-0.0461\times {\mathrm{P}}_{\mathrm{V}}{}^2 \left(\mathrm{NDVI}\ge 0.7\right)\end{array}\right.}_{}$$

(2)

where ε_water, ε_building and ε_natural denote the water emissivity, building emissivity and natural surface emissivity. P_V denotes the vegetation coverage by Equation (4) [47].

$${\mathrm{P}}_{\mathrm{V}}={(\frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }})}^2$$

(3)

where NDVI is the Normalized Difference Vegetation Index.

$$\mathrm{NDVI}=\left(\frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red}}\right)$$

(4)

$${\mathrm{T}}_{\mathrm{B}}=\frac{{\mathrm{K}}_2}{\ln \left(\frac{{\mathrm{K}}_1}{{\mathrm{L}}_{\mathrm{T}}}+1\right)}$$

(5)

where L_T denotes the surface-leaving radiance; K₁ and K₂ are the calibration parameters of the thermal infrared band, respectively. For Landsat 5 TM, K₁ = 607.76 Wm² sr⁻¹ µm⁻¹, K₂ = 1260.56 K; for Landsat 8 OLI/TIRS, K₁ = 774.89 W m² sr⁻¹ µm⁻¹, K₂ = 1321.08 K.

The final step was to calculate the surface temperature based on the brightness temperature of the thermal infrared band and the land surface emissivity [48]. Moreover, we converted the Kelvin unit of the LST to a Celsius unit.

$$\mathrm{LST}=\frac{{\mathrm{T}}_{\mathrm{B}}}{1+\left(\mathsf{\lambda }\times \frac{{\mathrm{T}}_{\mathrm{B}}}{\mathsf{\rho }}\right)\ln \mathsf{\epsilon }}$$

(6)

where T_B is the brightness temperature in the thermal infrared band, λ is the center wavelength of the thermal infrared band, and the center wavelength of the Landsat 5 TM thermal infrared band (band 6) is 11.5 µm. The Landsat 8 OLI/TIRS thermal infrared band (band 10) has a center wavelength of 10.9 µm, ρ = h × c/σ (1.438 × 10⁻² mK).

3.2. Extraction of the Land Surface Components (LSC)

This study classified the LSC of each city based on the spectral index of the remote sensing images, and four types were identified, namely, impervious surface area, urban green space, water body, and others. First, we calculated the Modified Normalized Difference Water Index (MNDWI) (Equation (7)) [49]. Based on the Google Earth images, we determined the segmentation threshold and extracted and masked the water body.

$$\mathrm{MNDWI}=\frac{\mathrm{Green}-\mathrm{SWIR}1}{\mathrm{Green}+\mathrm{SWIR}1}$$

(7)

where Green, SWIR1 denote the second and the fifth band of the Landsat 5 TM and the third and the sixth band of the Landsat 8 OLI/TIRS. The biophysical composition index (BCI) was used to extract the ISA of each city [50]. This method was based on the concept of the V–I–S (Vegetation–Impervious–Soil) triangle model proposed by Ridd [51], which assumed that the urban land surface included three parts, i.e., vegetation, ISA, and soil. It is an environmental index that can distinguish urban features [7]. Before calculation of the BCI index, the tasseled cap transformation (TC transformation) was performed. The above-mentioned image after the water body mask analysis was used. The specific calculation is as follows:

$$\mathrm{BCI}=\frac{\frac{\mathrm{H}+\mathrm{L}}{2}-\mathrm{V}}{\frac{\mathrm{H}+\mathrm{L}}{2}+\mathrm{V}}$$

(8)

where H is the normalized brightness component (TC1) with high reflectivity, L is the normalized humidity component (TC3) with low reflectivity, and V is the normalized greenness component (TC2) of the vegetation. H, L, and V can be calculated as:

$$\mathrm{H}=\frac{\mathrm{TC}1-\mathrm{TC}{1}_{\min }}{\mathrm{TC}{1}_{\max }-\mathrm{TC}{1}_{\min }}$$

(9)

$$\mathrm{V}=\frac{\mathrm{TC}2-\mathrm{TC}{2}_{\min }}{\mathrm{TC}{2}_{\max }-\mathrm{TC}{2}_{\min }}$$

(10)

$$\mathrm{L}=\frac{\mathrm{TC}3-\mathrm{TC}{3}_{\min }}{\mathrm{TC}{3}_{\max }-\mathrm{TC}{3}_{\min }}$$

(11)

where TC1, TC2, and TC3 are the first three TC components obtained after the TC transformation of the remote sensing data. TC_min and TC_max are the minimum and maximum values of each TC component respectively.

In the BCI, the ISA value is positive and relatively high, the value of the bare soil is close to 0, and the value of the vegetation is low or negative, which can effectively be distinguished from other land cover types [50]. The BCI can effectively reflect the biophysical composition of the urbanized environment. Compared with the NDBI (normalized difference built-up index), the BCI can better distinguish between soil and ISA with high albedo [7]. The BCI has advantage over other indices in the analysis of ISA. After calculation of the BCI of each city based on Equations (8)–(11), the segmentation threshold was determined, and the ISA was extracted.

The final step was to extract the green space of each city from the remote sensing images based on the NDVI. In this step, we masked the water bodies and ISA from the remote sensing images, and the segmentation threshold was determined by comparing the original remote sensing image and the Google Earth images. The green vegetation area was identified as a green space. The remaining parts that were not classified as ISA, green spaces and water bodies, were classified as others.

In summary, four surface components were extracted for each city: ISA, green spaces, water bodies, and others. Comparison was performed between the classification results of the components of the urbanized land surface from 2000 to 2015, and historical images of the same period in Google Earth. To test the accuracy of the extraction results, we selected 200 random sample points for each city, and the overall classification accuracy of each city reached as high as 80%.

3.3. Spatial Pattern of the Land Surface Components

In order to compare impacts of the composition of various land surface components and relevant spatial pattern on LST, we characterized the spatial modes of the LSC using kernel density estimation [19] and landscape metrics [52]. Kernel density estimation is mainly used to calculate the density of point elements or line elements and their neighboring points. This method takes the characteristics of the spatial proximity of the feature elements into account, and can calculate the contribution of the surrounding points to match the proximity effect of the LST caused by the surrounding pixels [19,53]. The kernel density estimation can be written as follows:

$$\mathrm{D}\left({\mathrm{x}}_0\right)=\frac{1}{\mathrm{n}}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\mathrm{K}\left(\frac{{\mathrm{x}}_0-\mathrm{x}\left(\mathrm{i}\right)}{\mathrm{r}}\right)$$

(12)

where K() is the kernel density function; r is the search radius which quantifies the distance from the estimated point element to the sample point x_o; $\mathrm{x}\left(\mathrm{i}\right)$ represents the neighbor within the circular neighbor areas; n is the number of the point elements within the search radius. The kernel density estimation results are constrained within the range of 0–100% through standardization (Equation (13)), and the results of the spatial density of the LSC can be computed as:

$$\mathrm{ISAD}\left(\mathrm{UGSD}\right)=\frac{\mathrm{D}\left(\mathrm{x}\right)-\min \left(\mathrm{D}\left(\mathrm{x}\right)\right)}{\max \left(\mathrm{D}\left(\mathrm{x}\right)\right)-\min \left(\mathrm{D}\left(\mathrm{x}\right)\right)}\times 100\%$$

(13)

The ISAD and UGSD show the spatial density of the ISA and UGS pixels neighboring to the ISA and UGS within a certain search radius. In this current study, based on three principles, i.e., being important in practice and theory [21,54], being easy to calculate and explain [1,21], and being minimum in redundancy [21,55,56], we adopted six landscape configuration metrics that reflected the spatial area, size, shape, and aggregation degree of the ISA and the UGS, i.e., the mean patch shape index (SHAPE_MN), the mean patch size (AREA_MN), the area-weighted fractal dimension index (FRAC_AM), the largest patch index (LPI), the landscape division index (LDI), and the aggregation index (AI). The above-mentioned landscape metrics are described in

Table 1. The landscape metrics used to delineate the shape configuration of the land surface components (LSC) analyzed in this study.

Landscape Metrics (Abbreviation)	Description	Unit (Value Range)
Mean patch shape index (SHAPE_MN)	The value of a given patch type divided by the total number of patches	None
Mean patch size (AREA_MN)	The average area of a given patch type within the study unit	Hectare
Area-weighted fractal dimension index (FRAC_AM)	The measure of the spatial shape complexity of a certain type of patch	None
Largest patch index (LPI)	The proportion of the largest patch of a given patch type divided by the total landscape area	Percent (0 < LPI ≤ 100)
Landscape division index (LDI)	The difference between the maximum value of the diversity index and the calculated value	None
Aggregation index (AI)	The number of similar adjacencies of the corresponding type divided by the maximum value when the type is maximally clustered into one patch	Percent (0 ≤ AI ≤ 100)

and calculated using Fragstats 4.2 software [57].

3.4. Pearson Correlation Analysis

The Pearson correlation coefficient is a measure of the linear correlation between two different variables [58]. This study used the Pearson correlation coefficient to quantify the relationship between the distribution density of the LSC of each city and the landscape metrics and the LST from 2000 to 2015.

3.5. Spatial Regression Model

The ordinary least squares (OLS) multiple linear regression model and spatial regression model were used to study the influence of the distribution patterns of the LSC (such as ISA and UGS) on LST changes [21,59]. The OLS multiple linear regression model assumed that the error terms are independent. The Moran I test was used to evaluate whether there was spatial autocorrelation between the error terms of the OLS model. The Moran index (Moran’s I) of the Moran I test ranges from −1 to 1. When Moran’s I is greater than 0, it means that the data have a positive spatial correlation. The larger the value, the more obvious the spatial correlation is; when Moran’s I is less than 0, it indicates that the data present a negative spatial correlation, and the smaller the value, the greater the spatial difference is; when Moran’s I is 0, the spatial pattern is random. If there was a spatial autocorrelation between the error terms of the OLS multiple linear regression model (P < 0.01), the spatial regression model fused with the spatial autocorrelation and the OLS multiple linear regression model would be introduced to jointly analyze impacts of the spatial distribution pattern of the LSC on the LST [21,60]. The spatial regression model includes the Spatial Lag Model (SLM) whose response variable y is spatial autocorrelation and the Spatial Error Model (SEM) whose error term is spatial autocorrelation. The spatial lag model is presented as:

$$\mathrm{y}=\mathsf{\rho }{\mathrm{W}}_1\mathrm{y}+\mathsf{\beta }\mathrm{x}+\mathsf{ϵ}$$

(14)

where y is the response variable; x is the explanatory variable; W₁ is the spatial weight matrix reflecting the spatial trend of the response variable y; β represents the spatial regression coefficient of the explanatory variable x; ϵ is the error distribution; ρ is the spatial lag coefficient and its value ranges between 0 and 1. The closer ρ is to 1, the more similar the value of the response variable y is to its neighboring areas. The spatial error model assumes that a spatial effect cannot be fully explained by the explanatory variable x in the error term. The spatial error model is as follows:

$$\mathrm{y}={\mathsf{\lambda }\mathrm{W}}_2\mathsf{\mu }+\mathsf{\beta }\mathrm{x}+\mathsf{ϵ}$$

(15)

where y is the response variable; x is the explanatory variable; W₂ is the spatial weight matrix reflecting the spatial trend of the residuals; μ represents the error term of the spatial change; β represents the spatial regression coefficient of the explanatory variable x; ϵ is the error distribution, and λ is the spatial error coefficient whose value ranges between 0 and 1. The closer that the value λ is to 1, the more similar the value of the response variable y is in adjacent areas.

The determination of the appropriate spatial regression model for use in this study was by application of the following criteria [61]: whether or not the residuals were independent, whether the Lagrange Multiplier (LM) and the Robust Lagrange Multiplier (R-LM) were statistically significant, R² values, and AIC (the Akaike’s Information Criterion). Based on above-mentioned criteria, the spatial error model was found to be the right model for this current study. In particular, the spatial error model with the maximum likelihood method was selected for this research. The above-mentioned regression analysis was performed using Geoda software and R Spdep package.

3.6. Variance Partitioning

Variance partitioning is used to quantify the contribution rate of the spatial pattern of the LSC to the LST changes. The variance of the response variable is decomposed into independent or joint explanatory variables (or variable groups) to explain different parts of the impacts of LSC on LST changes [62,63,64]. In this study, based on the spatial regression model, variance partitioning was used to quantify the relative importance of the impact of the compositions and spatial distribution of the ISA and UGS of each city on the LST changes [21,58]. The fractional contribution of the ISA and UGS to the LST was subdivided into four parts: (1) the unique effect of ISAD (or UGSD); (2) the unique effect of ISA (or UGSD) spatial configuration metrics; (3) the joint effects of ISAD (or UGSD) and ISA (or UGS) spatial configuration metrics; (4) unexplainable parts. The above-mentioned variance partitioning analysis was performed using Geoda software and R Spdep package.

3.7. Modeling of Future Changes in the LSC

In this study, we used the FLUS model to simulate changes in the LSC of each city considered in this study, and predicted the future LSC changes of each city using the Markov chain model [28]. First, the multi-layer feedforward artificial neural network algorithm (BP-ANN) was used to obtain the suitability probability (the probability that one region is suitable for urbanization) based on land use changes and other various driving factors such as terrain, traffic, location, and policies. Meanwhile, based on historical surface composition data, the Markov chain model was used to quantify the magnitude of LSC change during a specific year [61,65]. This model quantified the amount of change in LSC by comparing LSC of two cities during two specific periods. In this analysis, a transition probability matrix was obtained to quantify the probability that one land surface component category transformed to other land surface component categories. Then, the LSC data of each city in 2007 was taken as input (initial year), and the LSC data in 2015 was taken as the variable to verify the modeling accuracy of the FLUS model. The changes of the LSC in 2023 was predicted. The above-mentioned entire analysis was performed using GeoSOS-FLUS software (http://www.geosimulation.cn/flus.html) accessed on 9 January 2021.

3.8. Prediction of the LST

In this study, the multi-layer feed forward back propagation neural network (MFPNN) method of BP algorithm (BP neural network) was used to model and predict the future LST based on the historical LST records [66]. The BP neural network consists of an input layer, an output layer and one or more hidden layers. According to the complexity of the modeling object, we chose an appropriate network structure to realize the mapping of any nonlinear function from the input layer to the output layer. When using the BP neural network to predict the future LST based on the historical LST records, the LST data from 2000 to 2007 were used as input, and the LST data in 2015 were used as output. By setting the initial weight, learning rate, decay rate, maximum number of iterations and other parameters, we constructed the training network and tested the network until the network performance analysis and prediction accuracy met the requirements. Then we used the trained network to predict the LST of each city in 2023. The above calculation process was carried out by Matlab software. The analysis procedure is shown as Figure 2.

4. Results and Discussion

4.1. Spatial Pattern of the LST and LSC

4.1.1. Changes in LST during 2000–2015

The spatiotemporal changes in the LST from 2000 to 2015 based on Landsat 5 TM and Landsat 8 OLI/TIRS remote sensing images were demonstrated in Figure 3. In general, the LST over the cities considered in this study was increasing. In the BTH urban agglomeration, from 2000 to 2015, the LST of supercity Beijing was generally higher than that of megacity Tianjin, and large city Langfang. When compared with the LST during 2000–2007, the maximum and average LST over Beijing increased significantly in 2015, being 56.7 °C and 32.5 °C, respectively. The maximum LST was 6.5 °C and 14.6 °C higher than in 2000 and 2007, respectively. The average LST was 7.2 °C higher than in 2000, and 8.8 °C higher than in 2007. The spatial patterns of the LST over Beijing during the past three years were similar, i.e., the high-temperature area was observed mainly in the main urban area of southeast Beijing and the LST decreased radially toward the surrounding suburbs. The maximum LST experienced no significant change in Tianjin from 2000 to 2015, being around 42 °C. The average LST showed an increasing trend year by year, being respectively 27.3 °C, 28.1 °C and 29.5 °C. In addition, the high-temperature area of Tianjin expanded year by year from the central city to the surrounding suburbs. The LST of Langfang was concentrated in some small regions in the central and southern parts of Langfang. The LST in 2015 over Langfang was higher than during 2000–2007, and the maximum LST in 2015 was 5.7 °C and 7 °C higher, respectively, than during 2000–2007. The average LST in 2015 was 29.2 °C, which was 2.6 °C and 5.4 °C higher than in 2000 and 2007, respectively. In the YRD urban agglomeration, from 2000 to 2007, the LST in Shanghai (the supercity) was generally higher than Ningbo (the megacity) and Nanjing (the large city). The high LST was found mainly in central Shanghai and gradually spread to the surrounding areas. The LST of Ningbo showed an increasing trend from 2000 to 2015 with the highest LST of 36.7 ℃, 44.8 ℃ and 49.7 °C, and the average LST of 24.6 °C, 30.4 °C and 32.2 °C. The high LST was found in the eastern parts of Ningbo, and the regions with high LST continued to expand outward. The highest LST in Nanjing was observed in 2015 with a maximum LST of 42.1 °C, which was 0.3 °C and 0.5 °C higher than in 2000 and 2007, respectively. The areas with highest LST in Nanjing were concentrated along the north and south banks of the Yangtze River. In the PRD urban agglomeration, the LST over Guangzhou was generally higher than that in Dongguan and Zhongshan. From 2000 to 2015, the LST in Guangzhou increased year by year. The highest LST regimes were 36.1 °C, 40.2 °C and 44.1 °C, and the average LSTs were 20.8 °C, 22.9 °C and 26.4 °C. The regional difference between high and low temperatures was remarkable. The highest LST in Dongguan was observed in 2015, with the highest LST of 40.7 °C and the average LST of 28.3 °C. From 2000 to 2015, the highest LST in the Zhongshan area was generally around 34–35 °C. The regions with high LST were concentrated in the northwest and central Zhongshan, and gradually extended to the periphery.

According to the distribution of LST in different scales of the cities, our study found there is a certain relationship between city scale and LST. This result is consistent with the research of Su et al. [67]. In addition to driving LST change, internal factors such as natural elements (climatic zone, vegetation and waterbodies) and urban morphological factors such as city scale, are the external driving forces affecting LST. The scale of the city not only directly affects the urban thermal environment by changing the physical properties of the underlying surface, but also indirectly affects the LST through changes in urban ventilation, traffic demand, energy consumption, and contact with surrounding areas. The expansion of cities of different scales and levels has transformed a large number of natural surfaces into impervious surfaces. This leads to a decrease in latent heat flux and an increase in sensible heat flux, which in turn causes an increase in LST [68].

4.1.2. Changes in Spatial Pattern of the LSC during 2000–2015

The LSC from 2000 to 2015 were extracted for each city considered in this study, including ISA, UGS, water body and others. Classification of the LSC of cities of different sizes and levels across the BTH, YRD and PRD urban agglomerations was performed. The conversion of the LSC was analyzed using the land use transition matrix. The proportions of the LSC and the average LST for each LSC are shown in Figure 4, Figure 5 and Figure 6. From 2000 to 2015, the percentage of ISA in Beijing increased from 10.9% to 24.7%, and expanded ISA occurred mainly in central Beijing and continued to expand from central Beijing to the surrounding areas. Due to the expansion and encroachment of ISA, the percentage of UGS decreased annually from 80.3% to 69.4%. The highest LST was observed mainly over the ISA. From 2000 to 2015, the average LST across the ISA in Beijing was 31.6 °C, 26.9 °C and 37.6 °C, respectively. The average LST across the UGS was 24 °C, 22.8 °C and 30.7 °C, respectively. The difference in LST over ISA and UGS was about 4–8 °C, indicating that UGS could effectively alleviate the UHI effect. The ISA over Tianjin gradually expanded from the central city to the periphery. In 2015, the percentage of the ISA accounted for 28%. The percentage of UGS decreased, accounting for 40.1%. The average LST was 25.5 °C. This conversion process of the land surface components or modification of the land surface features caused significantly rising LST. Meanwhile, from 2000 to 2015, the ISA in Langfang also showed an increase from 12.2% to 18.7%, and a decrease in the proportion of the UGS from 87.6% to 57.8%. The average LST over the water bodies and UGS was generally lower than that over the ISA-dominated regions, showing a significant mitigation effect by the natural land surface, such as water bodies and vegetation, on the UHI effect.

Shanghai, the supercity located in the YRD, experienced a significant encroachment of ISA from 2000 to 2015 with increased percentage of ISA from 22.7% to 39.8% and decreased UGS from 66.6% to 39.8%. The difference in LST between ISA and UGS from 2000 to 2015 was 2.4 °C, 2.5 °C and 1.6 °C, respectively. From 2000 to 2015, the ISA of Ningbo was observed mainly in the northern and eastern parts of Ningbo, gradually expanding to the periphery. The percentage of ISA increased from 8.4% to 18.5%, and the percentage of UGS decreased from 87.3% to 62%. The average LST from 2000 to 2015 over the ISA was 27.2 °C, 33.8 °C, and 35.1 °C, showing an increasing trend year by year. Spatial distribution of ISA change was highly similar to that of the high LST. The ISA of Nanjing expanded along the north and south banks of the Yangtze River. By 2015, the ISA of Nanjing had increased to 25.1%, and the UGS accounted for 56.2%.

The ISA of Guangzhou, the supercity located in the PRD, continued to expand from the central parts to the north and southeast parts of Guangzhou from 2000 to 2015, and the percentage of ISA increased to 21.2%. The percentage of UGS decreased from 80% to 65.8%. The difference of the average LST between ISA and UGS showed an amplifying trend, being respectively 2.1 °C, 2.6 °C and 3.3 °C. The spatial pattern of ISA in Dongguan changed significantly from 2000 to 2015. Taking the main urban area as the center, the ISA encroached along the main traffic lines. The ISA increased from 28.3% in 2000 to 44.9% in 2015. The UGS in 2015 reduced to 29.7%. In addition, the difference in the average LST between ISA and UGS increased annually from 2000 to 2015, being, respectively, 1.6 °C, 3 °C and 3.7 °C. The ISA of the northern and central-southern areas of Zhongshan continued to expand to the periphery. From 2000 to 2015, ISA accounted for 16.3%, 23.3%, and 31%, respectively, and UGS accounted for 50.6%, 43%, and 35.1%, respectively. The LST over the ISA was 1.78 °C, 2.15 °C and 1.75 °C higher, than over the UGS. Additionally, the average LST corresponding to the UGS was similar to the LST of the water body.

In summary, from 2000 to 2015, the proportion of ISA in cities of different scales and levels expanded, and the proportion of UGS decreased. The LST corresponding to ISA was generally higher than that of UGS. The average LST corresponding to water bodies and UGS was generally lower than the average LST corresponding to ISA, which fully demonstrated the mitigation effect of natural surfaces such as water bodies and vegetation on the SUHI.

This paper mainly studied the relationship between ISA, UGS in the LSC and LST, therefore, the BCI index and NDVI index were mainly used to extract ISA and UGS. This classification method is consistent with Estoque et al. [45]. Overall, the LST of ISA in all cities was generally more than 1 °C higher than UGS. Other scholars have also studied the LST difference between ISA and UGS. Estoque et al. found the mean LST of impervious surface was about 3 °C higher than that of green space in Bangkok, Jakarta and Manila. Bokaie et al. [69] found the difference between the mean LST of ISA was more than 6 °C above UGS. Among the LSC, ISA and UGS are important factors that affect urban LST and thermal environment. Increasing urban greening can help alleviate the urban heat island effect.

4.2. Impacts of Spatial Pattern of LSC on LST

4.2.1. Relation between LSC and LST

The relationship between ISA density and UGS density and LST across three major urban agglomerations with different urban sizes and levels is shown in Figure 7, Figure 8 and Figure 9, respectively. We found significantly positive relation between LST and ISA density and significantly negative relation between UGS density and LST (P < 0.001). These findings corresponded closely with previous research results [12,70,71] and indicated that ISA could enhance UHI effects, and UGS could help to alleviate UHI effects. Rational vegetation planning in urban planning can effectively alleviate the urban thermal environment [70]. It can be seen from Figure 7 that ISA had the strongest impacts on LST change in Beijing with slope values of 0.135 and 0.116 in 2000, and 2015, respectively, which were higher than those in Tianjin and Langfang. Similarly, the impacts of UGS on LST in Beijing were also the strongest, with slope values of −0.105, −0.105, and −0.103, respectively, from 2000 to 2015, being stronger than those of Tianjin and Langfang. These observations indicated that Beijing, a city with larger urban sizes and levels than Tianjin and Langfang, was subject to stronger UHI effects than Tianjin and Langfang. By comparison, the slope values of relation between UGS and LST were generally smaller than between ISA and LST, implying higher impacts of ISA on LST, than of UGS on LST [45]. Figure 8 depicts the relationship between ISA density, UGS density and LST in Shanghai, Ningbo and Nanjing. The Pearson correlation coefficients (R) and slope values of the regressive relation between ISA density and LST were generally higher than that of UGS and LST. These results also indicated more fractional contribution of the impermeable land surface to LST changes than other driving factors. The impact of impervious land surfaces and UGS on LST changes varied in cities of different sizes and levels in different years. The impact of ISA on LST in Ningbo was generally stronger (the slope values were 0.059, 0.079, and 0.087, from 2000 to 2015) than Shanghai and Nanjing.

The relationships between ISA density, UGS density and LST in Guangzhou, Dongguan and Zhongshan from 2000 to 2015 is shown in Figure 9. Stronger impacts of the ISA on LST changes were observed in Guangzhou than in Dongguan and Zhongshan. The slope values of relation between ISA and LST in Guangzhou were, respectively, 0.048, 0.059 and 0.072 during 2000–2015, and were followed by Dongguan and Zhongshan, respectively. Meanwhile, the impact of the UGS on LST change in Guangzhou was generally stronger than in Dongguan and Zhongshan. The slope values of the relation between UGS and LST were respectively −0.031, −0.047 and −0.055, from 2000 to 2015, respectively. The correlation coefficients and slope values of the relation between ISA and LST were higher than those between UGS and LST.

In summary, the slope of UGS and LST under the three major urban agglomerations was generally lower than the slope of ISA and LST. Compared with UGS, ISA had a greater impact on LST. The results of this influence were consistent between 2000 and 2015. By comparing all cities, we found that in cities with larger urban sizes and levels, the impacts of ISA and UGS on LST were usually higher than smaller cities in size and scale.

4.2.2. Relation between Spatial Configuration of LSC and LST

In this study, six landscape metrics were adopted, i.e., SHAPE_MN, AERA_MN, FRAC_AM, LPI, LDI, and AI. These landscape metrics represented the shape, size, fractal dimension, proportion of the largest patch area, dominance, and aggregation of the LSC. They were selected to fully depict impacts of spatial configuration of LSC on LST changes.

Table 2. Analysis results of the landscape configuration metrics of the ISA within different urban sizes and levels considered in this study. Meanings of these abbreviations of the variables can be referred to in Table 1.

City	Year	SHAPE_MN	AERA_MN	FRAC_AM	LPI	LDI	AI
Beijing	2000	1.27	2.96	1.29	36.02	0.87	78.69
	2007	1.26	4.34	1.33	51.12	0.74	83.49
	2015	1.23	5.40	1.39	69.56	0.52	84.64
Tianjin	2000	1.27	3.37	1.31	43.58	0.81	80.81
	2007	1.27	2.63	1.31	42.46	0.82	78.43
	2015	1.23	5.25	1.37	63.08	0.60	86.71
Langfang	2000	1.24	2.67	1.22	6.27	0.99	79.43
	2007	1.24	2.39	1.22	7.49	0.99	78.31
	2015	1.25	2.81	1.24	7.74	0.99	78.60
Shanghai	2000	1.23	2.46	1.35	59.89	0.64	77.81
	2007	1.19	3.13	1.36	46.93	0.76	82.02
	2015	1.25	4.32	1.39	46.45	0.74	81.32
Ningbo	2000	1.23	1.98	1.26	14.30	0.96	77.28
	2007	1.25	2.59	1.29	18.80	0.94	78.26
	2015	1.27	2.90	1.33	25.34	0.93	75.69
Nanjing	2000	1.18	0.86	1.25	27.92	0.91	66.37
	2007	1.29	1.77	1.30	33.04	0.88	71.61
	2015	1.27	2.04	1.33	26.16	0.91	70.85
Guangzhou	2000	1.24	1.63	1.29	29.91	0.91	71.69
	2007	1.27	2.28	1.33	46.59	0.78	74.90
	2015	1.25	2.72	1.37	58.83	0.65	76.68
Dongguan	2000	1.24	3.40	1.35	23.35	0.89	78.39
	2007	1.24	4.73	1.40	50.63	0.70	80.70
	2015	1.21	6.63	1.44	90.06	0.19	84.91
Zhongshan	2000	1.23	1.25	1.28	20.75	0.94	68.61
	2007	1.25	1.78	1.29	21.16	0.93	72.34
	2015	1.26	2.61	1.35	51.14	0.73	76.26

and

Table 3. The analysis results of the landscape configuration metrics of the urban green space for different urban sizes and levels considered in this study.

City	Year	SHAPE_MN	AERA_MN	FRAC_AM	LPI	LDI	AI
Beijing	2000	1.19	40.04	1.40	96.12	0.08	96.25
	2007	1.21	33.91	1.39	91.88	0.16	95.81
	2015	1.26	19.89	1.37	80.30	0.35	94.19
Tianjin	2000	1.32	13.34	1.30	9.10	0.97	89.85
	2007	1.24	19.77	1.37	39.95	0.80	92.49
	2015	1.28	12.12	1.28	8.46	0.97	91.01
Langfang	2000	1.19	32.41	1.39	31.79	0.80	92.53
	2007	1.14	46.92	1.40	81.51	0.31	95.25
	2015	1.27	11.32	1.33	11.73	0.96	88.67
Shanghai	2000	1.25	18.92	1.35	17.68	0.92	91.48
	2007	1.29	8.76	1.33	20.45	0.94	86.30
	2015	1.28	4.65	1.29	7.99	0.98	82.68
Ningbo	2000	1.18	89.40	1.35	79.43	0.33	97.41
	2007	1.22	42.97	1.37	76.76	0.38	95.90
	2015	1.26	13.13	1.34	59.64	0.61	93.26
Nanjing	2000	1.30	14.59	1.39	20.68	0.89	87.07
	2007	1.28	30.23	1.36	37.57	0.80	93.23
	2015	1.28	11.24	1.37	17.47	0.90	87.56
Guangzhou	2000	1.30	13.53	1.36	62.50	0.60	91.30
	2007	1.27	24.60	1.37	78.37	0.38	93.93
	2015	1.29	14.34	1.35	69.68	0.51	92.66
Dongguan	2000	1.28	6.87	1.25	17.41	0.85	88.14
	2007	1.30	5.10	1.22	13.71	0.96	86.37
	2015	1.27	3.64	1.22	16.01	0.95	85.08
Zhongshan	2000	1.27	8.68	1.28	40.66	0.82	89.89
	2007	1.32	5.16	1.27	38.74	0.84	85.50
	2015	1.28	3.91	1.24	36.7	0.86	85.22

display the average values of the urban ISA and UGS landscape configuration metrics respectively for each city. Figure 10 shows the correlation between the urban ISA and UGS landscape configuration metrics and the LST in cities with different sizes from 2000 to 2015. Generally speaking, the LST and patch shape configuration index of the ISA were subject to significant positive correlation [1]. LDI reflects the degree of patch separation. Except for LDI, the LST and all other landscape metrics were subject to significant negative correlation [45,60].

Specifically, the three cities in the BTH urban agglomeration witnessed increased ISA from 2000 to 2015, and the patches tended to be complex and clustered. The landscape configuration metrics of the ISA in Beijing were generally higher than Tianjin and Langfang, indicating that the ISA patch area in Beijing was larger, more complex, and more concentrated when compared with Tianjing and Langfang. This was also the reason why the LST over the ISA in Beijing was generally higher than that of Tianjin and Langfang. From 2000 to 2015, contrary to the ISA, the UGS landscape configuration metrics decreased year by year, the area gradually decreased, and the degree of aggregation decreased. Beijing was dominated by a relatively higher degree of spatial concentration, dominance, and shape complexity of the UGS. It can be seen from Figure 10 that the highest correlation can be observed between LST, LPI and AREA_MN, indicating a stronger UHI given the larger ISA patch area, a more concentrated ISA and higher patch of the ISA [45]. Among the six UGS landscape configuration metrics, a higher correlation stood between LST and LPI, AREA_MN and AI. The larger and uninterrupted UGS patches had a more obvious cooling effect on the LST and hence a more obvious cold island effect [60,72]. In addition, the correlation between the LST and the ISA landscape configuration metrics in each year was generally higher than that between LST and the UGS landscape configuration metrics, indicating that the ISA had a higher impact on the LST.

In YRD, the period from 2000 to 2007 witnessed increased ISA landscape configuration metrics over the cities within the PRD with increased ISA patch area and increased complexity and concentration of ISA. However, the UGS landscape configuration metrics showed a downward trend in 2015, indicating that the area of UGS decreased, and the degree of concentration and complexity of UGS also decreased. By comparison, the concentration and complexity of the ISA in Shanghai was higher than in Ningbo, and the concentration and complexity of the ISA in Nanjing was the lowest. However, Ningbo had the largest UGS patch area, and the highest complexity and concentration degree of UGS, followed by Shanghai and Nanjing. The spatial configuration of the ISA and the UGS significantly affected the LST. Among all the indicators that reflected the landscape configuration, the LST was generally highly correlated with the LPI, AERA_MN, and AI of the ISA and the UGS. The increase in the complexity and concentration of the ISA shape meant it could absorb more solar radiation, leading to an increase in the LST. At the same time, the more complex the shape of the ISA and the greater the difference in ISA shape, the greater the benefit in terms of energy exchange between the interior of the city and the vegetation, resulting in cooling effects [73]. This phenomenon clearly indicates the complex thermal environment characteristics of the urbanized regions [1].

In the PRD urban agglomeration, the period from 2000 to 2015 witnessed a continually enhanced level of the ISA landscapes in Guangzhou, Dongguan, and Zhongshan, and this process was accompanied by increased ISA patch area, a higher degree of ISA shape complexity and an incremented aggregation of ISA. Specifically, we observed the highest ISA landscape configuration metrics in Dongguan, followed by Guangzhou and Zhongshan. We observed a decreased shape complexity of the UGS patches and a continually lower aggregation degree of the UGS patches. Meanwhile, Guangzhou had the highest green space landscape configuration metrics, followed by Dongguan and Zhongshan. Furthermore, we identified a higher correlation between LST and LPI, AREA_MN, AI of the ISA and the UGS, indicating that the larger and more complex patches of the LSC usually had stronger impacts on LST. In addition, the marginal characteristics of the LSC patches also had a certain degree of influence on the LST changes. For example, the increased edge of the LSC patches could strengthen the energy flow in the UGS and its surroundings, resulting in a decrease in the LST [60].

In summary, the LST and the ISA landscape configuration metrics were in significantly positive correlation. The LST and the UGS landscape configuration metrics were in significantly negative correlation [74,75]. Meanwhile, correlation between LST and ISA landscape configuration metrics was stronger than between LST and the UGS landscape configuration metrics, indicating that the ISA had a stronger impact on the LST changes when compared with UGS landscape configuration metrics. With the increase in the ISA landscape configuration metrics and the decrease in the UGS landscape configuration metrics, the warming effects of the ISA and the cooling effects of the UGS were affected by many factors such as city size, vegetation type and climatic conditions [76]. In this study, the larger the size of the city, the higher the average value of the ISA landscape configuration metrics, and hence the higher the LST. This conclusion that the shape configuration of LSC in different urban scales is related to LST is more consistent with the existing research, such as Zhou et al. [77]. This may be due to the greater complexity of the forms of megacities and supercities, leading to a more significant impact on the LST [68].

It may also be attributed to the higher degree of intensification of urban land development, which leads to an increase in anthropogenic heat emissions and an increase in LST in cities with high levels of scales [78]. In addition, we observed a high correlation between LST and LPI, AREA_MN, AI of ISA and UGS, implying that a concentrated larger LSC patch area can drive higher LST. In this sense, scientific urban planning considering optimal design and planning of UGS and ISA will greatly alleviate UHI, and it is particularly important for large size cities such as Beijing, Shanghai and Guangzhou in this study. In general, in addition to controlling urban expansion, proper dispersion of LSC and multi-center distribution can effectively alleviate SUHI.

4.2.3. The Contribution of the Composition and Configuration of the LSC to LST

This study also considers the combined effects of the composition and configuration of the LSC to explore the influence of the distribution pattern of the LSC on LST. The OLS multiple linear regression model and the spatial error model were both used to establish the relation between LST and the composition and configuration of ISA (or UGS). The results are shown in Appendix A

Table A2. The modeling results by the OLS model for the LST, the ISAD and the ISA landscape configuration metrics within different cities considered in this study from 2000 to 2015. * indicates p < 0.05; ** indicates p < 0.01; *** indicates p < 0.001.

City	Year	ISAD	SHAPE-MN	AREA-MN	FRAC-AM	LPI	LDI	AI	R2	Moran’I	AIC
Beijing	2000	0.674 ***	0.047 ***	0.064 ***	−0.158 ***	−0.059 **	0.296 ***	0.014 ***	0.592	0.724 ***	541094
	2007	0.541 ***	−0.058 ***	0.256 ***	−0.212 ***	−0.126 ***	0.304 ***	−0.026 ***	0.325	0.644 ***	682229
	2015	0.788 ***	−0.095 ***	0.077 ***	0.707 ***	−0.349 ***	−0.413 ***	0.03 ***	0.759	0.658 ***	393801
Tianjin	2000	0.777 ***	−0.172 ***	0.929 ***	−0.482 ***	−0.713 ***	0.501 ***	−0.021 ***	0.469	0.524 ***	382749
	2007	0.049 ***	−0.203 ***	0.016	0.125 ***	0.75 ***	0.024	−0.054 ***	0.456	0.491 ***	379852
	2015	0.874 ***	−0.083 ***	0.373 ***	0.025	−0.444 ***	0.089 ***	0.001	0.603	0.68 ***	329310
Langfang	2000	0.769 ***	−0.294 ***	1.138 ***	−0.578 ***	−0.777 ***	0.642 ***	−0.094 ***	0.551	0.579 ***	226556
	2007	0.593 ***	−0.403 ***	1.44 ***	−1.275 ***	−0.689 ***	1.296 ***	−0.161 ***	0.366	0.56 ***	264969
	2015	0.686 ***	−0.222 ***	0.996 ***	−0.891 ***	−0.453 ***	0.835 ***	−0.101 ***	0.497	0.658 ***	239078
Shanghai	2000	0.823 ***	−0.015	0.34 ***	−0.93 ***	0.128 ***	0.782 ***	−0.042 ***	0.607	0.437 ***	192958
	2007	0.909 ***	−0.093 ***	0.791 ***	−0.884 ***	−0.216 ***	0.798 ***	−0.057 ***	0.69	0.391 ***	168755
	2015	0.814 ***	−0.098 ***	0.684 ***	−0.326 ***	−0.362 ***	0.355 ***	−0.039 ***	0.61	0.47 ***	191986
Ningbo	2000	0.605 ***	0.015	0.132 ***	−0.164 ***	−0.169 ***	0.307 ***	0.018 **	0.405	0.645 ***	318874
	2007	0.816 ***	0.037 **	0.067 *	−0.068	−0.13 ***	0.154 ***	0.011 *	0.667	0.549 ***	235489
	2015	0.789 ***	0.014	0.149 ***	−0.174 ***	−0.105 ***	0.239 ***	0.009**	0.684	0.617 ***	229122
Nanjing	2000	0.722 ***	−0.107 ***	0.421 ***	−0.828 ***	−0.107 *	0.764 ***	−0.046 ***	0.471	0.588 ***	215496
	2007	0.48 ***	−0.24 ***	0.469 ***	0.184 *	−0.469 ***	0.0189	−0.009	0.194	0.863 ***	264010
	2015	0.014 ***	−0.109 ***	−0.371 ***	0.0403	0.911 ***	0.053	−0.052 ***	0.218	0.612 ***	257677
Guangzhou	2000	0.481 ***	0.074 ***	−0.028	−0.438 ***	0.11**	0.449 ***	0.003	0.296	0.599 ***	284862
	2007	0.592 ***	−0.011	0.145 ***	−0.328 ***	−0.021	0.433 ***	0.026 ***	0.427	0.708 ***	265101
	2015	0.774 ***	0.001	0.05*	−0.048	−0.036	0.183 ***	0.004	0.696	0.6 ***	189168
Dongguan	2000	0.472 ***	0.051 *	−0.092	−0.358 ***	0.29 ***	0.388 ***	−0.011	0.314	0.515 ***	92150.1
	2007	0.677 ***	0.018	−0.099 **	−0.018	0.169 ***	0.153 ***	−0.002	0.601	0.543 ***	73249.3
	2015	0.799 ***	−0.021	−0.001	0.121 ***	−0.007	−0.02	−0.000001	0.741	0.439 ***	55702.8
Zhongshan	2000	0.692 ***	−0.192 ***	0.551 ***	−0.854 ***	−0.095	0.716 ***	−0.083 ***	0.393	0.606 ***	59497.9
	2007	0.719 ***	−0.061 *	0.413 ***	−1.179 ***	0.171 **	0.987 ***	−0.092 ***	0.499	0.575 ***	55615.7
	2015	0.759 ***	−0.125 ***	0.351 ***	−0.412 ***	−0.025	0.447 ***	−0.064 ***	0.529	0.585 ***	52061.1

,

Table A3. The modeling results of the LST, the UGSD and the UGSD landscape configuration metrics within different cities considered in this study from 2000 to 2015. * indicates p < 0.05; ** indicates p < 0.01; *** indicates p < 0.001.

City	Year	UGSD	SHAPE-MN	AREA-MN	FRAC-AM	LPI	LDI	AI	R2	Moran’I	AIC
Beijing	2000	−0.67 ***	0.018 ***	0.141 ***	−0.329 ***	0.381 ***	0.567 ***	−0.067 ***	0.623	0.717 ***	519350
	2007	−0.017 ***	−0.072 ***	0.299 ***	−0.215 ***	−0.501 ***	0.286 ***	−0.036 ***	0.297	0.629 ***	693643
	2015	−0.784 ***	−0.033 ***	0.241 ***	−0.024 ***	−0.222 ***	0.029 *	0.042 ***	0.707	0.703 ***	448522
Tianjin	2000	−0.586 ***	0.175 ***	−0.867 ***	0.159 ***	0.562 ***	−0.305 ***	0.055 ***	0.338	0.578 ***	421253
	2007	−0.643 ***	0.011 *	−0.5 ***	0.376 ***	−0.044	−0.541 ***	−0.004	0.421	0.563 ***	390533
	2015	−0.695 ***	0.092 ***	−0.521 ***	0.087 ***	0.363 ***	−0.22 ***	−0.0007	0.514	0.697 ***	364348
Langfang	2000	−0.72 ***	0.039 ***	−0.699 ***	0.289 ***	−0.061 ***	−0.83 ***	−0.004 ***	0.532	0.568 ***	231241
	2007	−0.749 ***	−0.082 ***	0.203 ***	0.303 ***	−0.393 ***	−0.275 ***	0.031 ***	0.582	0.497 ***	199182
	2015	−0.75 ***	−0.033 ***	−0.425 ***	0.479 ***	−0.241 ***	−0.734 ***	0.007	0.584	0.53 ***	218060
Shanghai	2000	−0.745 ***	0.078 ***	−0.644 ***	0.28 ***	0.122 ***	−0.478 ***	0.03 ***	0.562	0.478 ***	203939
	2007	0.909 ***	−0.093 ***	0.791 ***	−0.884 ***	−0.216 ***	0.798 ***	−0.057 ***	0.69	0.391 ***	168755
	2015	−0.767 ***	−0.073 ***	−0.232 ***	0.688 ***	−0.191 ***	−0.633 ***	−0.016 **	0.523	0.527 ***	212543
Ningbo	2000	−0.444 ***	0.026 ***	0.109 ***	0.002	−0.394 ***	−0.073*	−0.045 ***	0.412	0.625 ***	317177
	2007	−0.715 ***	0.015 ***	−0.005	−0.092 ***	0.241 ***	0.204 ***	−0.065 ***	0.54	0.649 ***	279376
	2015	−0.809 ***	0.004	−0.072 *	−0.355 ***	0.563 ***	0.355 ***	−0.063 ***	0.605	0.66 ***	259549
Nanjing	2000	−0.5 ***	−0.009	−0.199 ***	0.014	0.079 *	−0.253 ***	−0.067 ***	0.304	0.651 ***	242417
	2007	0.002	0.027 ***	−0.037	0.17 ***	−0.746 ***	−0.445 ***	0.022 **	0.123	0.849 ***	272507
	2015	−0.303 ***	0.013	−0.115 ***	0.332 ***	−0.517 ***	−0.45 ***	0.042 ***	0.178	0.659 ***	262691
Guangzhou	2000	−0.437 ***	0.033 ***	0.245 ***	−0.676 ***	0.535 ***	0.824 ***	−0.094 ***	0.349	0.568 ***	275906
	2007	−0.542 ***	0.039 ***	−0.032	−0.168 ***	0.366 ***	0.374 ***	−0.048 ***	0.381	0.73 ***	274000
	2015	−0.728 ***	0.037 ***	−0.02	−0.313 ***	0.436 ***	0.445 ***	−0.046 ***	0.653	0.632 ***	204422
Dongguan	2000	−0.453 ***	−0.0002	0.259 ***	−0.494 ***	0.237 ***	0.589 ***	−0.079 ***	0.311	0.514 ***	92337.4
	2007	−0.687 ***	−0.103 ***	0.417 ***	−0.767 ***	0.374 ***	0.865 ***	−0.115 ***	0.567	0.564 ***	76395.8
	2015	−0.786 ***	−0.044*	0.243 ***	−0.354 ***	0.035	0.293 ***	−0.044 ***	0.651	0.524 ***	66765.5
Zhongshan	2000	−0.501 ***	0.013	−0.976 ***	0.607 ***	0.291 **	−0.73 ***	0.007	0.24	0.657 ***	65224
	2007	−0.605 ***	−0.135 ***	−0.392 ***	0.221 **	0.439 ***	−0.16 **	−0.072 ***	0.246	0.678 ***	66221
	2015	−0.54 ***	−0.043*	−0.35 ***	−0.416 ***	0.754 ***	0.366 ***	−0.072 ***	0.302	0.669 ***	61937.3

,

Table A4. The results by the SEM for the LST, the ISAD and the ISA landscape configuration metrics in different cities considered in this study from 2000 to 2015. * indicates p < 0.05; ** indicates p < 0.01; *** indicates p < 0.001.

City	Year	ISAD	SHAPE-MN	AREA-MN	FRAC-AM	LPI	LDI	AI	R2	AIC
Beijing	2000	0.708 ***	0.017 *	0.007	−0.144 ***	0.023 *	0.126 ***	0.02 ***	0.883	258028
	2007	0.688 ***	0.011	0.051 **	−0.177 ***	−0.029	0.14 ***	0.023 ***	0.737	473592
	2015	0.863 ***	−0.068 ***	0.096 ***	0.088 ***	−0.095 ***	−0.025 *	0.018 ***	0.906	183419
Tianjin	2000	0.953 ***	−0.122 ***	0.964 ***	−1.043 ***	−0.601 ***	0.899 ***	−0.069 ***	0.707	302773
	2007	0.008 ***	−0.321*	0.555 ***	−0.608 ***	0.522 ***	0.666 ***	−0.101 ***	0.682	309919
	2015	1.074 ***	−0.052 ***	0.482 ***	−0.336 ***	−0.446 ***	0.296 ***	−0.026 ***	0.864	181627
Langfang	2000	1.066 ***	−0.281 ***	1.241 ***	−1.65 ***	−0.716 ***	1.521 ***	−0.167 ***	0.784	163139
	2007	0.876 ***	−0.316 ***	1.339 ***	−1.95 ***	−0.588 ***	1.789 ***	−0.183 ***	0.682	205517
	2015	0.951 **	−0.192 ***	1.08 ***	−1.554 ***	−0.495 ***	1.368 ***	−0.144 ***	0.805	155246
Shanghai	2000	0.889 ***	0.014	0.374 ***	−0.664 ***	−0.129 ***	0.538 ***	−0.028 ***	0.746	159735
	2007	1.011 ***	−0.075 ***	0.723 ***	−0.814 ***	−0.278 ***	0.717 ***	−0.059 ***	0.78	142922
	2015	0.952 ***	−0.059 ***	0.589 ***	−0.476 ***	−0.271 ***	0.462 ***	−0.041 **	0.759	154920
Ningbo	2000	0.682 ***	0.08 ***	−0.061 **	−0.035	−0.037	0.016	0.029 ***	0.768	215893
	2007	0.888 ***	0.036**	−0.037	−0.074	−0.031	0.073 *	0.026 ***	0.827	166308
	2015	0.886 ***	0.016	0.007	−0.127 ***	0.003	0.113 ***	0.019 ***	0.863	139198
Nanjing	2000	0.873 ***	0.012	0.271 ***	−0.651 ***	0.079 *	−0.253 ***	−0.067 ***	0.746	158792
	2007	0.695 ***	0.028*	0.157 ***	−0.257 ***	−0.095 ***	0.227**	−0.009	0.877	102932
	2015	0.003 *	−0.103 ***	0.114 ***	−0.258 ***	0.326 ***	0.275 ***	−0.034 ***	0.685	187073
Guangzhou	2000	0.465 ***	0.059 ***	−0.026	−0.146 *	0.043	0.112	0.011*	0.696	209781
	2007	0.612 ***	0.024 *	0.03	−0.14 **	−0.008	0.118 **	0.004	0.823	154786
	2015	0.795 ***	0.032 ***	−0.019	−0.082**	0.012	0.071 **	0.013 ***	0.87	112694
Dongguan	2000	0.537 ***	0.063**	−0.075 *	−0.12	0.103 *	0.093	0.007	0.631	74375.7
	2007	0.73 ***	0.015	−0.095 ***	0.028	0.085 *	0.009	0.007	0.798	53227.4
	2015	0.865 ***	0.004	−0.041	−0.0007	−0.003	0.007	−0.044 ***	0.837	42713.4
Zhongshan	2000	0.249 ***	−0.077 **	0.619 ***	−0.795 ***	−0.269 **	0.648 ***	−0.051 ***	0.719	44037.4
	2007	0.817 *	−0.05 ***	0.448 ***	−0.772 **	−0.132 **	0.64 ***	−0.066 ***	0.762	40533.4
	2015	0.803 ***	−0.036 *	0.347 ***	−0.41 ***	−0.139 **	0.352 ***	−0.043 ***	0.788	36459.4

and

Table A5. The modeling results by the SEM of LST, UGSD and the UGS landscape configuration metrics in different cities considered in this study from 2000 to 2015. * indicates p < 0.05; ** indicates p < 0.01; *** indicates p < 0.001.

City	Year	UGSD	SHAPE-MN	AREA-MN	FRAC-AM	LPI	LDI	AI	R2	AIC
Beijing	2000	−0.7 ***	0.022 ***	−0.029 **	−0.017 ***	0.009	0.005	−0.021 ***	0.882	256111
	2007	−0.003**	0.085 ***	−0.052 ***	−0.015	−0.421 ***	−0.078 ***	−0.006	0.722	490252
	2015	−0.841 ***	−0.01 ***	0.011	0.047 ***	−0.063 ***	−0.075 ***	−0.02 ***	0.902	202098
Tianjin	2000	−0.712 ***	0.081 ***	−0.558 ***	0.481 ***	0.157 ***	−0.445 ***	0.045 ***	0.683	321000
	2007	−0.675 ***	0.014 **	−0.39 ***	0.275 ***	0.032 ***	−0.386 ***	0.016 ***	0.682	309919
	2015	−0.771 ***	0.023 ***	−0.323 ***	0.345 ***	0.114 ***	−0.299 ***	0.014 ***	0.846	205080
Langfang	2000	−0.891 ***	0.016 ***	−0.552 ***	0.335 ***	−0.044 *	−0.739 ***	0.023 ***	0.768	170234
	2007	−0.741 ***	0.102 ***	−0.76 ***	0.281 ***	0.043	−0.857 ***	0.0345 ***	0.678	206887
	2015	−0.852 **	−0.027 ***	−0.46 ***	0.439 ***	0.009	−0.581 ***	0.021 ***	0.784	162180
Shanghai	2000	−0.785 ***	0.04 ***	−0.398 ***	0.191 ***	0.191 ***	−0.3 ***	0.024 ***	0.737	164506
	2007	−0.774 ***	−0.043 ***	0.096 ***	0.05*	−0.072 **	−0.037	−0.022 ***	0.757	156928
	2015	−0.83 ***	−0.036 ***	−0.293 ***	0.463 ***	0.069 **	−0.406 ***	0.012*	0.74	164919
Ningbo	2000	−0.496 ***	0.033 ***	−0.028	0.041 ***	−0.273 ***	−0.139 ***	0.004	0.757	221564
	2007	−0.707 ***	0.013 **	−0.018	−0.016 *	0.041	0.031 *	−0.023 ***	0.813	181500
	2015	−0.853 ***	0.019 ***	−0.068 ***	−0.025	0.104 ***	−0.001	−0.015 ***	0.85	154105
Nanjing	2000	−0.621 ***	0.007	−0.176 ***	0.143 ***	0.013	−0.232 ***	−0.009	0.731	167733
	2007	0.003 *	0.04 ***	−0.148 ***	0.091 ***	−0.216 ***	−0.213 ***	0.011 **	0.863	114682
	2015	−0.43 ***	0.009	−0.109 ***	0.103 ***	−0.069 **	−0.164 ***	0.003	0.694	184050
Guangzhou	2000	−0.427 ***	0.032 ***	0.071 **	−0.175 ***	0.021	0.208 ***	−0.013*	0.693	209151
	2007	−0.498 ***	0.031 ***	−0.064 ***	0.009	−0.033	−0.047 **	0.002	0.819	157696
	2015	−0.678 ***	0.034 ***	−0.076 ***	−0.031 **	0.034	0.013	0.002	0.861	121155
Dongguan	2000	−0.479 ***	0.029 ***	0.106 *	−0.101 **	−0.083	0.106 *	−0.015	0.626	74827.3
	2007	−0.712 ***	−0.014	0.165 ***	−0.183 ***	−0.003	0.199 ***	−0.033 ***	0.789	54989.2
	2015	−0.783 ***	−0.019	0.099*	−0.172 ***	0.021	0.142 **	−0.015 *	0.818	48173.7
Zhongshan	2000	−0.619 ***	0.005	−0.595 ***	0.392	0.254 ***	−0.438 ***	0.015	0.723	44648.6
	2007	−0.611 ***	−0.011	−0.362 ***	0.298 **	0.157 **	−0.321 ***	0.01	0.742	43877.6
	2015	−0.41 ***	0.032	−0.288 ***	−0.004	0.264 ***	−0.007	0.005	0.763	40265.9

. In general, the ISAD and UGSD in each city were the main predictors of LST changes (Appendix A Table A2 and Table A3), which were positively and negatively correlated with LST, respectively. In terms of standard coefficients, the composition of ISA and UGS had more significant impact on LST than variables representing landscape configuration. Among the various ISA landscape configuration metrics, the positive correlation between the AREA_MN, which measures the complexity of the ISA, and the LDI, which measures the degree of ISA landscape separation, and LST, was relatively obvious, but the composition of the ISA had a stronger influence on LST than the landscape morphology. Among the various UGS landscape configuration metrics, the AREA_MN showed a more significant negative correlation with LST. The contribution of the UGS composition to the change of LST was generally higher than the UGS spatial form, while in some cities, the spatial form of UGS had a slightly significant impact on LST. However, it can be seen from Appendix A Table A2 and Table A3 that the Moran’I of each city from 2000 to 2015 was significant, which means that there was spatial autocorrelation in the error term of the OLS model. Therefore, this study also used SLM and SEM to continue to explore the influence of the distribution pattern of the LSC on LST. The residual error, Lagrange multiplier, robust Lagrange multiplier, R² values, and AIC criterion, were used to compare the two models, and finally select the SEM for this research.

Appendix ATable A4 and Table A5 show the SEM results of the impact of ISA and UGS spatial composition and configuration on LST, respectively. Overall, the results of the SEM were similar to the OLS. However, the R² values of the SEM were generally higher, and the AIC was significantly lower than the OLS. Moreover, the standard coefficients of the ISAD and UGSD prediction factors were relatively high, and the standard coefficients of each landscape configuration metrics were relatively low, which also illustrated the necessity of such studies to consider the effects of spatial autoregression and the superiority of the SEM compared with the OLS [21].

Variance partitioning was used to quantify the contribution rate of the LSC to LST, and the results are shown in Figure 11 and Figure 12. In general, the joint effect of the composition and configuration of the LSC had the highest contribution rate to the LST. However, in different cities and different years, the unique effects of the LSC composition and the unique effects of the LSC configuration had different impacts on LST.

In the BTH urban agglomeration, the unique effect of ISA composition in the supercity of Beijing gradually increased and had stronger impact on LST than the ISA landscape configuration, being 1.1%, 1.5% and 1.6% higher than ISA landscape configuration during 2000–2015. In 2000–2007, the impact of UGS landscape configuration on LST was slightly stronger than that of UGS composition. In the spatial error model, a lower standard coefficient of the UGS density reflected higher AREA_MN, FRAC_AM which reflected the complexity of the UGS shape, and higher LPI metrics of the concentration of UGS patches, indicating enhanced impacts of landscape configuration on LST given the increasingly concentrated and complicated shape of the UGS [69]. In Tianjin and Langfang, the composition of the ISA had a higher impact on the LST than the shape configuration of ISA. Meanwhile, the composition of the UGS had a higher impact on the LST than the shape configuration of the UGS. However, different from Beijing, the period of 2000 to 2015 witnessed weakening effects of the spatial pattern and landscape configuration of the ISA and UGS on LST over time. In addition, in the BTH urban agglomeration, the composition and shape configuration of the ISA and the UGS on the LST in Langfang were generally higher than that of Tianjin and Beijing, being 1–2% higher than Tianjing and 3–4% higher than Beijing.

In the YRD urban agglomeration, during the study period from 2000 to 2015, the ISA composition made a 2.49%, 3.27% and 3.4% higher contribution rate to the LST changes than the configuration of ISA in Shanghai. Additionally, the composition and configuration of the ISA had increasing impacts on LST change. Similarly, from 2000 to 2015, the composition of Shanghai’s UGS had 1.8%, 2.3%, and 2.1% higher impact on the LST than UGS configuration. Furthermore, during the period of 2000 to 2015, the composition of the ISA in Ningbo had 1.6%, 3.2% and 2.7% higher impact on the LST than the shape of the ISA. During 2000 to 2015, the impact of UGS composition on LST showed an upward trend, especially in 2007, and the impact of the UGS composition was the strongest. During 2000 to 2015, the impact of the ISA composition on the LST in Nanjing showed a decreasing trend. The impact of the ISA composition on the LST during 2000 and 2007 was 3.4% and 1.4% higher than that of the ISA configuration, respectively. However, in 2015, the impact of the ISA configuration was slightly higher than that of the ISA composition. The standard coefficient of the ISA density was lower than that of the ISA landscape configuration. Similarly, the impact of UGS composition on the LST also showed a decreasing trend during the period from 2000 to 2015.

The joint effects of the LSC and landscape configuration of Guangzhou, the supercity in the PRD urban agglomeration, on LST changes were remarkable. The composition of ISA and UGS had a slightly higher impact on LST than the shape configuration of the ISA and UGS. From 2000 to 2015, the effect of the ISA composition on the LST of Dongguan was stronger than that of the ISA configuration. The effect of the ISA composition increased over time. The effect of the UGS composition on the LST from 2000 to 2015 was also stronger than the shape configuration of the UGS. The effect of the UGS composition on the LST increased, particularly during 2007. The impact of ISA composition on the LST from 2000 to 2015 in Zhongshan was stronger than that of ISA configuration on LST. However, the contribution rate of ISA to LST showed a decreasing trend with time, indicating that with the expansion of the impervious surface area, the patches became more complex and concentrated, and the influence of ISA configuration on the LST increased to a certain extent. From 2000 to 2015, the contribution rate of the UGS composition to the LST was slightly higher than that of the UGS shape configuration, the influence of the UGS composition on the LST increased annually, and the impact of the UGS shape configuration on the LST decreased annually. All these results indicated that the decreased proportion of the UGS resulted in sporadic and fragmented distribution of the UGS, and the cooling effect of the UGS on the LST was significantly weakened.

The aforementioned results indicated that the joint effects of the composition and shape configuration of the LSC were the strongest on the LST. There was a difference in magnitude of the effects of the shape configuration of the LSC on LST [21]. When compared with the unique effect of the composition and the unique effect of configuration of the LSC on the LST, the unique effect of the LSC of each city was generally higher than the unique effect of the shape configuration of the LSC on LST [69]. The composition and configuration of the ISA had the potential to drive the increase in the LST, and the composition and configuration of the UGS had a cooling effect on the LST. However, in some specific cities during some specific periods, the configuration of the LSC was slightly higher than that of the LSC, which may be affected by the climate background. In different climatic conditions such as temperature, solar radiation, air pressure, etc., the proportion of LSC and spatial morphology of different cities have different impacts on the LST [79]. In addition, it was also affected by the difference in the morphology of the LSC of the year [80].

In summary, the results of this study show that LST is significantly positively correlated with ISAD and negatively correlated with UGSD. These conclusions are consistent with other research (i.e., Myint et al. [9]; Estoque et al. [45]). From the perspective of landscape configuration, the results of all the cities show that the correlation between LST and ISA landscape configuration is higher. This indicates that ISA has a stronger influence on LST than UGS. Moreover, the ISA patches of cities with higher scale levels are relatively larger, more complex, and more concentrated, which is also the reason for the higher LST. Our results are consistent with the research of Estoque et al. [45], Li et al. [60], and Zhang et al. [81]. In terms of the LSC of UGS, the UGS connected in patches with relatively large area is more helpful to reduce the LST. From the perspective of joint effect of ISAD (or UGSD) and ISA (or UGS) spatial configuration metrics, the results of all cities show that the joint effect of ISAD and ISA spatial configuration metrics have a stronger influence on LST than UGS has on LST. This joint effect is stronger than the unique effect, which is consistent with the results of Zhou et al. [21]. Moreover, results from the variation partitioning of all cities (Figure 11 and Figure 12) also indicate that LSC density plays a more important role than that of shape configuration of LSC.

4.3. Modeling of Future LST and LSC

4.3.1. Modeling of the Future LSC

This study used the spatial pattern of the LSC in each city in 2007 as the initial input. Additionally, we used 14 variables, i.e., distance from the city center, distance from the county center, distance from the town center, distance from the commercial center, distance from the subway, distance from the highway, distance from the provincial road, distance from the county road, distance from the village road, DEM, aspect, slope, GDP, population density, and spatial distribution of water bodies, and the neural network model (BP-ANN model) in modeling the suitability probability of the development of various LSC. In addition, based on the Markov chain model (Markov Model), we predicted the LSC. The LSC data in 2015 were used to test the validity of prediction of the spatial pattern in the LSC in 2023 (Figure 13). The kappa coefficient of the FLUS model simulation results of each city was generally higher than 0.75, and the overall accuracy was higher than 80%.

From the prediction results, under natural conditions, the spatial patterns of the LSC of three cities in the BTH urban agglomeration are different. In 2023, the ISA of Beijing will expand mainly in the southeast direction, encroaching on green space and resulting in a decrease in the proportion of the UGS to 63.8%, with no significant change in the proportion of water and other areas. The proportion of the ISA area in Tianjin will increase to 35.4%, expanding to the periphery mainly by central and coastal ports. This is consistent with the development model of Tianjin’s “double-center and one-axis”, that is, the central city and Binhai New Area are used as the “double centers”, and the traffic axis between the central city and Binhai New Area is the “one axis” for further development [82]. The proportion of the UGS will reduce to 47.4%. It is necessary to pay more attention to the protection of the UGS and the ecological environment during the urbanization processes due to the continuous development of construction land. In addition to expanding outwards along the original urban center, the ISA of Langfang will contain many scattered ISA areas, which will account for 26.3%, while the proportion of the UGS will decrease to 71.8%. Compared with 2000, the growing rate of the percentage of the ISA in Beijing, Tianjin and Langfang in 2023 will be 20.2%, 21.7% and 14.1%, respectively.

Under natural conditions, the percentage of the ISA in Shanghai, the supercity in the YRD urban agglomeration, will increase sharply from 22.7% in 2000 to 44.5% in 2023, and will expand substantially from the central city to the surrounding area. The ISA of Ningbo will expand to the periphery with the northern and central regions as the center, and gradually will connect to form a belt. From 2000 to 2023, the percentage of the ISA will increase from 8.4% to 26.8%, and the proportion of the UGS will decrease from 87.3% to 67.2%. The ISA of Nanjing not only will expand along the edge of the existing ISA, but will also fill the interior of the city. The ISA will grow rapidly and the percentage of ISA will increase from 11.5% in 2000 to 29.7% in 2023.

Under natural conditions, the expansion of the ISA of Guangzhou, the supercity in the PRD, will be centered on the northwest, central and south parts of the PRD. In comparison, the expansion of ISA to the south is more obvious than other directions. The ISA expansion to the north due to topographical factors is insignificant. From 2000 to 2023, the proportion of the ISA will increase from 12.2% to 25.7%, while the UGS will decrease from 80% to 65.3%. By 2023, the ISA will further expand to the periphery, accounting for 51.9%. The expansion of the ISA will encroach on the UGS and water bodies. From 2000 to 2023, the proportion of the UGS will decrease from 45.3 to 25.6%, and the percentage of water bodies will decrease from 8% to 6.8%. The ISA of Zhongshan will expand to the surroundings centered on the north and the middle, and there will be scattered flaky ISA in the southern area of Zhongshan. From 2000 to 2023, the proportion of ISA will increase from 16.3% to 35.4%, the UGS will reduce from 50.6% to 35%, and the proportion of water bodies will also decrease significantly, from 14.8% to 13.3%.

The simulation of the LSC in 2023 can play a certain guiding role for future urban development policy. For the above cities, compared with 2000, the ISA will increase sharply in 2023 at the expense of UGS and water bodies. Therefore, in the process of rapid urbanization, it is necessary to coordinate the proportion of various types of land in the city, protect the ecological environment, and coordinate development. In addition, by comparing the results of all cities, we found that the ISA of each city increased significantly from 2000 to 2015, but the increase in the percentage of ISA of cities of different sizes and levels differed from 2015 to 2023. Generally speaking, the above cities with higher scale levels show a weaker growth trend of the ISA, while the ISA of cities with relatively low scale levels increased more significantly than cities with larger sizes or scales.

In summary, the prediction results of each city indicate that the LST will continue to rise. In the future, while maintaining the rapid development of economic construction, it is necessary to take further measures to alleviate the increasingly serious SUHI. Based on the existing studies [83], it is possible to alleviate the SUHI with various measures such as rational planning of urban layout, reducing anthropogenic heat emissions, enhancing urban greening, and protecting urban water bodies.

4.3.2. Modeling of Future LST Changes

In this study, the BP ANN algorithm was used to simulate the LST in 2023 under the scenario of natural surface composition expansion. The future LST simulation results can provide a theoretical basis for the evolution of the future UHI effect and alleviation of the UHI effect. The LST data in 2015 was used as the verification data to verify the accuracy of the prediction results. The Pearson correlation coefficients between the predicted LST in 2015 and the observed LST were statistically significant, and the RMSE values were statistically low, indicating the validity of the modeling performance of the BP ANN algorithm used in this study. The prediction accuracy was used to simulate the LST in 2023. Among them, the Pearson correlation coefficients of the observed and the predicted LST in Beijing, Tianjin and Langfang in 2015 were 0.88, 0.77, and 0.74, respectively, and the RMSE values were 1.95, 1.74, and 1.5, respectively; the Pearson correlation coefficients of the observed and the predicted LST in Shanghai, Ningbo and Nanjing in 2015 were 0.86, 0.78 and 0.72, and the RMSE values were 0.89, 0.88 and 1.17, respectively. The Pearson correlation coefficients of the observed and predicted LST in Guangzhou, Dongguan and Zhongshan in 2015 were 0.84, 0.88 and 0.76, respectively, and the RMSE values were 1.22, 1.19 and 1.16, respectively.

Under the scenario of natural LSC expansion, the simulation results of the LST in cities with different urban sizes in the BTH, YRD, and PRD urban agglomerations in 2023 are shown in Figure 14, Figure 15 and Figure 16, respectively. The regions with high predicted LST in Beijing (the supercity in BTH) were highly similar to those dominated by ISA, and the regions with low LST were consistent with those dominated by UGS and water bodies. From 2015 to 2023, with the expansion of the ISA, the LST will also be increasing. The average LST will increase from 32.5 °C to 35.3 °C, and the average LST over the ISA will increase from 37.6 °C to 40.96 °C. The temperature difference over the UGS will increase from 6.83 °C to 8.42 °C. There are mainly two centers in the regions with high LST in Tianjin, including the central city and the coastal port area, which are gradually connected in a flaky shape. From 2015 to 2023, the average LST will rise slightly from 29.5 °C to 30.6 °C, and the average LST over the ISA will not change significantly. However, the LST over UGS will increase from 28.3 °C to 30.1 °C. The difference in LST between the ISA and the UGS is 1.59 °C. The regions with high LST in Langfang follow a fragmented spatial distribution pattern. The average LST will increase from 29.2 °C to 32.8 °C from 2015 to 2023. In 2015, the difference in LST between the non-commitment land surface and the UGS was 3.62 °C. However, no remarkable difference, only 0.53 °C, in LST over ISA and UGS, can be identified. In contrast, the higher the scale of the BTH urban agglomeration, the more obvious the temperature rise of the ISA, and the greater the difference in the LST between ISA and the UGS.

In Shanghai, the supercity in the YRD, regions with high LST expand outward from the main urban area. The LST over the ISA and UGS in 2015 and 2023 is 1.6 °C and 1.05 °C respectively. The regions with high LST are found mainly in the north and central Ningbo and are similar to those dominated by ISA. The average LST in 2023 will be 34.4 °C, which is an increase of 2.2 °C when compared with 2015. The average LST over ISA will increase by 1.3 °C, and the average LST over the UGS will increase by 2 °C. The areas with high LST in Nanjing are still concentrated on the north and south banks of the Yangtze River, and are expanding to the south. Similar to 2015, the difference in LST between ISA and UGS in 2023 will be 0.71 °C.

The spatial distribution of the LST in Guangzhou, the supercity in the PRD, shows a decreasing trend from southwest to northeast. Compared with 2015, the average LST in 2023 over the ISA and UGS will increase significantly, respectively, being about 3.2 °C, 2.95 °C and 3.28 °C. The regions with high LST are widespread and low LST will be observed mainly over the water bodies. In 2015–2023, the average LST will increase from 28.3 °C to 31.3 °C, and the LST over the ISA will increase from 30.1 °C to 32.2 °C. The LST over the UGS will increase from 26.4 °C to 29.6 °C. In Zhongshan, regions with high LST were concentrated mainly in the northern and central regions, expanding to the south, consistent with the spatial distribution of ISA. The low LST was found mainly in regions dominated by UGS and water bodies. From 2015 to 2023, the difference in LST between ISA and UGS will be 1.75 °C and 1.62 °C, respectively.

In summary, by comparing the LST in 2015 and the predicted LST in 2023, the higher the scale of the city, the more significant the difference in LST between ISA and UGS. Comparing the spatial pattern of the LST and LSC, the expansion of the ISA density remarkably affects the UHI effect, and the cooling effect of water bodies and UGS is more apparent in cities of large urban sizes. In the process of urban planning, an increase of the proportion of water bodies and UGS will have a very significant effect on alleviation of the UHI effect.

5. Conclusions

Based on Landsat 5 TM and Landsat 8 OLI/TIRS remote sensing images, spatial regression models, variance partitioning and spatial statistical methods, we analyzed the impacts of the composition and shape configuration of the LSC on LST in cities with different urban sizes. This study explored impacts of the LSC on LST of different sizes and levels of cities across three typical urban agglomerations with different climate backgrounds across China. Additionally, artificial neural network algorithms and other methods were used to simulate the future spatial pattern of LSC and LST, providing a scientific ground for mitigating UHI effects and sustainable urban planning. Major interesting and important findings and conclusions were obtained as follows:

(1) From 2000 to 2015, the LST of each city generally showed an increasing trend over time. Additionally, we identified continuous expansion of ISA across the cities considered in this study. With the encroachment and expansion of ISA, the proportion of UGS decreased year by year.

(2) Comparison of LST over LSC indicated a highly similar spatial pattern of high LST across cities with different urban sizes. The LST over ISA was generally higher than that across UGS. The difference in LST between ISA and UGS in different cities in different years was different, the lowest was 0.96 °C and the highest was 7.96 °C. The LST of the UGS and the water body was similar, indicating that these two factors significantly reduce the LST.

(3) The LST and the ISA density were in significantly positive correlation, while LST and the UGS density were in significantly negative correlation. Among LSC, the influence of ISA on LST was greatest. Additionally, we also found that in larger sized cities, ISA had greater impacts on LST, than in smaller sized cities. This may be due to the fact that cities with higher scales and levels have a faster urbanization process and a larger proportion of ISA, so the impact on LST is more significant. In addition, based on the landscape metrics reflecting the shape of the LSC, we analyzed the effect of the configuration of the LSC on the LST. We hold the opinion that the positive correlation between the LST and the ISA landscape configuration metrics was stronger than the negative correlation between LST and the UGS landscape configuration metrics. Therefore, ISA had stronger impact on LST than UGS. Meanwhile, the ISA landscape configuration metrics generally increased, and the UGS landscape configuration metrics decreased over time. The higher the scale of the city, the higher the average ISA landscape configuration metrics; the corresponding LST was also relatively high. It showed that the larger size and scale of the city was closely related to the high degree of influence of the LSC on the LST.

(4) We investigated the influence of the spatial pattern of the LSC on the LST and found that whether ISA or UGS, the joint effect of the composition and configuration of ISA and/or UGS had the most significant impact on the LST. The joint effect differed in magnitude from the unique effect of the composition or configuration of the LSC on the LST. Generally speaking, the unique effect of the ISA and the UGS had a stronger influence on the LST than the unique effect of its configuration.

(5) We simulated the future LSC and LST in both space and time in each city considered in this study. We found that the average LST in the future will also show an upward trend, the proportion of the ISA will continue to expand, and the proportion of the UGS will continue to shrink. The higher the size of the city, the weaker the growing percentage of ISA, while the lower the size of the city, the stronger the growing percentage of ISA. By comparing the average LST over the LSC in the future, it is interesting to find that the larger the size of a city, the more significant is the difference in the LST between ISA and UGS.