Spatiotemporal Evolution of the Urban Thermal Environment Effect and Its Influencing Factors: A Case Study ofBeijing, China

Rapid urbanization has led to significant changes in land surface temperature (LST), which in turn affect the urban thermal environment effect and the health of residents. Exploring the causes of the urban thermal environment effect will provide guidance for promoting sustainable urban development. The spatiotemporal evolution of the urban thermal environment effect within the sixth ring road of Beijing was analyzed by inversion of remote sensing data to obtain the LST in 2004, 2009, 2014, and 2019. In addition, based on multivariate spatial data, we applied the standard deviation ellipse (SDE), spatial principal component analysis (PCA), and other methods to analyze and identify the relationships between the urban thermal environment effect and its influencing factors. The results show that from 2004 to 2019, the spatial distribution of urban development and LST within the sixth ring road of Beijing were closely related, the heat island area showed a small increasing trend, and differences in the thermal environment effect between different administrative regions in different periods were obvious. The main factors affecting the urban thermal environment effect were urban construction intensity, vegetation and water bodies, socioeconomic activities, and geomorphology. It is noteworthy that human factors had a greater impact than natural factors. Among them, the positive effect of the normalized difference impervious surface index (NDBBI) and the negative effect of the fractional vegetation cover (FVC) were the most prominent. This study provides theoretical support for mitigating the urban thermal environment effect and promoting sustainable urban development.


Introduction
The process of urbanization has a profound impact on the ecological environment of a region. Urbanization is first reflected spatially in the replacement of natural landscapes by urban landscapes and in changes to the overall state of the regional ecological environment by affecting the material cycle and energy flow of the ecosystem [1]. The urban thermal environment effect is an important element that can characterize the energy exchange between anthropogenic activities and natural systems and is a manifestation of the impact of urbanization on a regional climate [2]. Moreover, it offers an important basis for understanding how ecosystems in urban areas respond to landscape evolution. First proposed by Howard in 1818 [2], the urban thermal environment effect is specifically reflected in the significantly higher temperatures in cities than in the outer suburbs, which are usually measured by the urban heat island (UHI) effect. Since it was first proposed, it has received extensive academic attention [3][4][5].
Land surface temperature (LST) is considered to be an important parameter that directly controls the urban thermal environment effect. LST reflects all surface-atmosphere interactions and energy fluxes between the atmosphere and the ground and plays an important role in studying thermal environment effects [6,7]. Since the implementation of The Reform and Opening-up Policy, China's accelerated rate of urbanization in recent decades has led to an increase in the LST in urban areas, with direct impacts on urban air quality [8,9], local climate [10,11], energy utilization [12], thermal comfort [13][14][15], and biological community [16]. Furthermore, a high-temperature environment can pose a threat to human health from both the physiological and psychological perspectives. High temperatures are not only a direct cause of heat stroke and dehydration but also affect sleep quality and induce irritability and mental health disorders in individuals, thereby increasing the risk of depression and other psychological disorders [17][18][19]. These effects degrade the quality of the thermal environment, which in turn slows urban development. The establishment of ecological environmental protection and sustainable, environment-friendly living practices for humans continues to garner increasing attention [20,21], and alleviating the thermal environment effect and designing sustainable urban living environments have become urgent urban ecological environment goals [22,23]. Therefore, clarifying the correlation between the urban thermal environment effect on urban human settlements and various influencing factors will serve to advance the current understanding regarding mutual feedback between the landscape pattern and ecological processes in urban areas. These advances will not only benefit the sustainable development of urban environments but also assist in improving the quality of human life [24,25].
Numerous studies have been conducted on the spatiotemporal evolution [26,27], formation causes [28][29][30], morphology and structure [31], environmental impacts [32,33], mechanisms [34,35], and simulations [36,37] of the urban thermal environmental effect. The factors influencing thermal environment effects have been identified as [38] (1) land cover, including land use cover change [39,40], the normalized difference vegetation index (NDVI) [41,42], normalized difference built-up index [43,44], and normalized difference impervious surface index (NDBBI) [45,46]; (2) the characterization of urban development using socioeconomic factors, including the gross domestic product (GDP) [47], population density [32,48], and industrial production activities [49]; and (3) natural factors, mainly including elevation and slope. However, it is difficult to use a single geographical element to analyze the urban thermal environment effect comprehensively. Spatial principal component analysis (PCA) can eliminate correlation effects between the evaluation indices. Moreover, numerous rating indices can be simplified into a smaller number of integrated indices while retaining the vast majority of information. In the integrated evaluation function, the weight of each principal component is its contribution rate, which can reflect the proportion of the amount of information of the principal component containing the original data to the total amount of information so that the determination of the weights is objective and reasonable [50]. This method is conducive to the comprehensive analysis of the relationship between various influencing factors and the urban thermal environment effect. In addition, many studies have emphasized the importance of landscape composition and configuration in mitigating the urban thermal environment effect [24,51,52]. However, the lack of comparative studies of different administrative regions within cities at different time periods makes it difficult to make specific recommendations for the mitigation of thermal environment effects in different regions.
Beijing is one of the most populous metropolitan areas in the world. Intensive human activities and urban development within the sixth ring road have led to a significant increase in LST, posing a major threat to human comfort and the thermal environment. Therefore, in this study, the spatiotemporal evolution of the urban thermal environment effect within the sixth ring road of Beijing was analyzed by inversion of remote sensing data to obtain the LST in 2004, 2009, 2014, and 2019. In addition, based on multivariate spatial data (Landsat 4-5 Thematic Mapper (TM), Landsat 8 Operational Land Imager Thermal Infrared (OLI TIRS), digital elevation model (DEM), road network, point of interest (POI), and National Polar-Orbiting Operational Environmental Satellite System Preparatory Project/Visible Infrared Imaging Radiometer Suite (NPP/VIIRS) nighttime light data), we obtained various natural and human factors closely related to the urban thermal environment effect and applied the standard deviation ellipse (SDE), spatial PCA, and other methods to identify and analyze the relationships between the urban thermal environment effect and its influencing factors. The results of this study provide theoretical support for mitigating the urban thermal environment effect and promoting sustainable urban development.

Study Area
Beijing is located in the northern part of the North China Plain and is the capital of the People's Republic of China ( Figure 1). The city has a temperate continental monsoon climate with an annual average temperature of 12.3 • C and an annual average precipitation of 600 mm. In the past 20 years, the average temperature and precipitation in Beijing have fluctuated, showing an overall upward trend ( Figure 2). According to the Master Plan of Beijing City (2016-2035), the total administrative area of Beijing comprises 16,410 km 2 , and the permanent population of Beijing is limited to 23 million. The sixth ring road represents the boundary of the most concentrated area for various activities and the population of Beijing. The area within the sixth ring road is a key location for the construction of urban ecological infrastructure and offers a reasonable area for conducting research on the urban thermal environment effect.
cific recommendations for the mitigation of thermal environment effects in different regions.
Beijing is one of the most populous metropolitan areas in the world. Intensive human activities and urban development within the sixth ring road have led to a significant increase in LST, posing a major threat to human comfort and the thermal environment. Therefore, in this study, the spatiotemporal evolution of the urban thermal environment effect within the sixth ring road of Beijing was analyzed by inversion of remote sensing data to obtain the LST in 2004, 2009, 2014, and 2019. In addition, based on multivariate spatial data (Landsat 4-5 Thematic Mapper (TM), Landsat 8 Operational Land Imager Thermal Infrared (OLI TIRS), digital elevation model (DEM), road network, point of interest (POI), and National Polar-Orbiting Operational Environmental Satellite System Preparatory Project/Visible Infrared Imaging Radiometer Suite (NPP/VIIRS) nighttime light data), we obtained various natural and human factors closely related to the urban thermal environment effect and applied the standard deviation ellipse (SDE), spatial PCA, and other methods to identify and analyze the relationships between the urban thermal environment effect and its influencing factors. The results of this study provide theoretical support for mitigating the urban thermal environment effect and promoting sustainable urban development.

Study Area
Beijing is located in the northern part of the North China Plain and is the capital of the People's Republic of China ( Figure 1). The city has a temperate continental monsoon climate with an annual average temperature of 12.3 °C and an annual average precipitation of 600 mm. In the past 20 years, the average temperature and precipitation in Beijing have fluctuated, showing an overall upward trend ( Figure 2). According to the Master Plan of Beijing City (2016-2035), the total administrative area of Beijing comprises 16,410 km 2 , and the permanent population of Beijing is limited to 23 million. The sixth ring road represents the boundary of the most concentrated area for various activities and the population of Beijing. The area within the sixth ring road is a key location for the construction of urban ecological infrastructure and offers a reasonable area for conducting research on the urban thermal environment effect.

Data Sources
In this study, Landsat 5 TM and Landsat 8 OLI multispectral and thermal infrared remote sensing images were obtained from the United States Geological Survey (USGS) data center (https://glovis.usgs.gov/, accessed on 3 October 2021). The time of the collected remote sensing image data was hot summers with clear weather, and the cloud cover in the study area was less than 5%. The remote sensing image data were preprocessed by ENVI 5.3 software using radiometric calibration, atmospheric correction, and area of interest clipping [45]. Then, the images were classified into seven land cover categories, namely cropland, forest, shrub, grassland, water, impervious land, and barren land ( Figure 3). Additional data included POI data, DEM data, various road network data, NPP/VIIRS nighttime light data, and a vector map of the administrative boundaries within the sixth ring road (Table 1).

Data Sources
In this study, Landsat 5 TM and Landsat 8 OLI multispectral and thermal infrared remote sensing images were obtained from the United States Geological Survey (USGS) data center (https://glovis.usgs.gov/, accessed on 3 October 2021). The time of the collected remote sensing image data was hot summers with clear weather, and the cloud cover in the study area was less than 5%. The remote sensing image data were preprocessed by ENVI 5.3 software using radiometric calibration, atmospheric correction, and area of interest clipping [45]. Then, the images were classified into seven land cover categories, namely cropland, forest, shrub, grassland, water, impervious land, and barren land ( Figure 3). Additional data included POI data, DEM data, various road network data, NPP/VIIRS nighttime light data, and a vector map of the administrative boundaries within the sixth ring road (Table 1).

Methods
This study used remote sensing data obtained within the sixth ring road of Beijing from 2004 to 2019 to extract the urban surface temperatures in 2004, 2009, 2014, and 2019 and then to analyze the spatial and temporal evolution patterns of the urban thermal environment effect. Based on previous studies, various natural and human factors close-

Methods
This study used remote sensing data obtained within the sixth ring road of Beijing from 2004 to 2019 to extract the urban surface temperatures in 2004, 2009, 2014, and 2019 and then to analyze the spatial and temporal evolution patterns of the urban thermal environment effect. Based on previous studies, various natural and human factors closely related to the urban thermal environment were obtained from multivariate spatial data (Landsat 4-5 TM, Landsat 8 OLI TIRS, DEM, road network, POI, and NPP/VIIRS nighttime light data), including the normalized difference impervious surface index (NDISI), fractional vegetation cover (FVC), modified normalized difference water index (MNDWI), normalized difference between bare land and building index (NDBBI), surface albedo, NPP/VIIRS nighttime light data (NTL), elevation, slope, road network density index, and POI. We applied the SDE, spatial PCA, and other methods to analyze and identify the relationships between the urban thermal environment effect and its influencing factors. A flowchart depicting the methodology employed in this study is provided in Figure 4.

Retrieval of Surface Temperature
LST is considered to be one of the important parameters that directly controls the urban thermal environment effect. LST reflects all surface-atmosphere interactions and energy fluxes between the atmosphere and the ground and plays an important role in studying thermal environment effects. In this study, the atmospheric correction method was used to extract the LST using Equations (1) and (2) [53] as follows: where B(T s ) denotes the blackbody thermal radiation brightness, T s denotes the surface temperature, L down and L up represent the downward and upward radiation brightness of the atmosphere, respectively, and K 1 and K 2 are constants. The L up , L down , K 1 , and K 2 values for the Landsat 5 TM image were 0.38 W/(m 2 ·sr·µm), 3.65 W/(m 2 ·sr·µm), 607.76, and 1260.56, respectively. The L up , L down , K 1 , and K 2 values for the Landsat 8 OLI image were 1.25 W/(m 2 ·sr·µm), 2.37 W/(m 2 ·sr·µm), 774.89, and 1321.08, respectively. Using the standard deviation method [54], the study area was divided into seven thermal class zones based on LST: extremely high temperature, high temperature, relatively high temperature, medium temperature, relatively low temperature, low temperature, and extremely low temperature (Table 2). Since the LST values in the extremely high temperature, high temperature, relatively high temperature, and medium temperature zones of the study area were generally higher than the average LST, these areas were considered urban heat island zones [55].

Retrieval of Surface Temperature
LST is considered to be one of the important parameters that directly controls the urban thermal environment effect. LST reflects all surface-atmosphere interactions and energy fluxes between the atmosphere and the ground and plays an important role in studying thermal environment effects. In this study, the atmospheric correction method was used to extract the LST using Equations (1) and (2) [53] as follows: where B(Ts) denotes the blackbody thermal radiation brightness, Ts denotes the surface temperature, Ldown and Lup represent the downward and upward radiation brightness of Note: u represents the average surface temperature in the study area, and std represents the standard deviation of the surface temperature.

Acquisition of Surface Information
(1) The NDISI is the normalized difference impervious surface index. The impervious surface is the main component of the underlying surface of a city; it can prevent water from infiltrating the ground and is calculated as follows [56]: where MNDWI is the modified normalized difference water index and NIR, SWIR1, and TIRS1 correspond to bands 5, 6, and 10 in the Landsat 8 remote sensing images, respectively.
(2) The FVC is the proportion of the vertical projection area of ground vegetation to the total calculated area, calculated as where NDVI is the normalized differential vegetation index, NDVI min and NDVI max are the normalized vegetation index values of bare soil and vegetation and take the values of 0.05 and 0.70, respectively, and Red and NIR correspond to bands 4 and 5, respectively, in the Landsat 8 remote sensing images.
(3) The MNDWI is the normalized difference processing carried out for image bands containing water body information, calculated as where Green and SWIR1 correspond to bands 3 and 6, respectively, in the Landsat 8 remote sensing images.
(4) The NDBBI is the normalized difference processing of remote sensing images containing information such as bare land and building land after removing information on water body and vegetation, and the associated formula is [57] where Green, NIR, and SWIR2 correspond to bands 3, 5, and 7, respectively, in the Landsat 8 remote sensing images. (5) The surface albedo is the ratio of incident solar radiation energy that can be reflected by the surface, calculated as [58] where a 2 , a 4 , a 5 , a 6 , and a 7 correspond to the gray values of bands 2, 4, 5, 6, and 7, respectively, in the Landsat 8 remote sensing images. (6) The NPP/VIIRS nighttime light data were obtained from The National Oceanic and Atmospheric Administration of the United States (http://www.ngdc.noaa.gov/, accessed on 3 October 2021) with a spatial resolution of approximately 500 m.
(7) The elevation was the vertical distance from the ground point to the height starting surface. (8) The slope was the amplitude of the surface fluctuation in a certain area. Both the elevation and slope were obtained from the DEM data.
(9) The social and economic activity index was used to evaluate the impact of regional social and economic activities on an urban thermal environment. The POI and road network data of the research area were used to represent the social and economic activity index [55].
The POI refers to the point data in the Internet electronic map, which basically contains four attributes: name, address, coordinates, and category. It mainly involves the types of restaurants, entertainment, accommodation, medical care, education, and sports that can characterize socioeconomic activities. The road network includes highways, national roads, provincial roads, county roads, urban roads, and other roads at all levels. The density analysis tool in the ArcGIS 10.8 software was used to process the POI and road network data to extract the socioeconomic activity index.

Spatial PCA
Spatial PCA uses a geographic information system (GIS) to assign each spatial variable to a matrix and distributes the influence degree of relevant spatial variables on dependent variables to the corresponding principal component factors through PCA [59]. In this study, SPSS 24.0 and ArcGIS were used for the correlation analysis and spatial PCA of various natural and human factors.

SDE
The SDE was first proposed by Lefever [60] to reveal the spatial distribution characteristics of geographical elements and has since been widely used in the field of spatial statistics. With the rapid development of GIS technology, the SDE method based on geographic information has become a conventional spatial statistics module tool [61]. In this study, the SDE method was used to identify the migration trajectory of the heat island.

Spatial Distribution of the Urban Thermal Environment
The heat island within the sixth ring road of Beijing from 2004 to 2019 presented a polycentric irregular agglomeration distribution ( Figure 5). The four temperature zones (extremely high, high, relatively high, and medium temperature) representing heat island areas were concentrated in built-up urban areas with dense buildings and populations. The highest temperature occurred in urban development zones. The cold island areas were primarily distributed in rivers, inner lakes, and high mountains, with the lowest temperatures detected at the junction of the Mentougou and Haidian districts. The alpine area of the study area comprised low and extremely low temperature zones, with an obvious "cold island" effect. According to a comparison of four time points, the heat island distribution from 2004 to 2019 was relatively the same as that of the built-up area.
In 2004, the heat island area was relatively small and primarily distributed in the central and southwestern portions of the research area. The highest and lowest surface temperatures reached 41.02 • C and 21.22 • C, respectively ( Figure 6, Table 3). The heat island areas were concentrated in the Chaoyang, Fengtai, Haidian, and Daxing districts, with a mean surface temperature of 32.898 • C. The extremely high temperature areas were concentrated in southern Chaoyang, eastern Fengtai, central Daxing, and eastern Fangshan, accounting for 24.5%, 24.43%, 16.57%, and 15.63%, respectively, of the total extremely high temperature area. The cold island areas were concentrated in Haidian, Chaoyang, Changping, and Tongzhou. The extremely low temperature areas were concentrated in eastern Mentougou and in central Haidian and Changping, accounting for 43.67%, 31.24%, and 22.78%, respectively, of the total extremely low temperature area.
In 2009, the heat island area expanded compared with that in 2004, and the mean surface temperature of the heat island area was 27.959 • C. The highest surface temperature within the study area reached 51.96 • C, while the lowest reached 15.76 • C. After 2004, the speed of urban development within the research area accelerated, resulting in significant development within the Chaoyang, Fengtai, Haidian, and Daxing districts. Consequently, the concentrated heat island areas were relatively maintained within these districts, with the proportion of the heat island area not significantly changing in each district, while the spatial location of the heat island shifted with an extension direction that corresponded to that of urban construction (i.e., extending from southwest to northeast). Meanwhile, the heat island areas of Changping, Shunyi, and Tongzhou expanded, resulting in their even distribution throughout the research area. The extremely high temperature agglomeration areas were primarily distributed in central Daxing and eastern Fangshan, accounting for 23.60% and 14.36%, respectively, of the total extremely high temperature area. The change in the spatial location of the cold island was insignificant and was concentrated in Chaoyang, Haidian, Changping, and Tongzhou. The cold island proportion in Chaoyang increased, that in Changping and Tongzhou decreased, and that in Haidian remained stable. The extremely low temperature agglomeration areas were concentrated in central and western Haidian, eastern Mentougou, and central Changping, accounting for 53.12%, 34.60%, and 7.22%, respectively, of the total extremely low temperature area.   Table 3). The heat island areas were concentrated in the Chaoyang, Fengtai, Haidian, and Daxing districts, with a mean surface temperature of 32.898 °C. The extremely high temperature areas were concentrated in southern Chaoyang, eastern Fengtai, central Daxing, and eastern Fangshan, accounting for 24.5%, 24.43%, 16.57%, and 15.63%, respectively, of the total extremely high temperature area. The cold island areas were concentrated in Haidian, Chaoyang, Changping, and Tongzhou. The extremely low temperature areas were concentrated in eastern Mentougou and in central Haidian and Changping, accounting for 43.67%, 31.24%, and 22.78%, respectively, of the total extremely low temperature area.      The Beijing Municipal Government Report 2010 stated that Beijing's environmental problems required urgent attention. It introduced strict control policies such as total population control and traffic restrictions, as well as measures to protect water sources and strengthen planning constraints. The impacts are seen in the 2014 data, with the heat island area slightly decreasing, and the mean surface temperature of the heat island area was 35.734 • C. In 2014, the highest and lowest surface temperatures of the research area reached 51.96 • C and 25.65 • C, respectively. The heat island areas continued to be primarily concentrated in the Chaoyang, Fengtai, Haidian, and Daxing districts. The heat island proportion increased slightly in Chaoyang, decreased slightly in Daxing, and remained stable in Haidian and Fengtai. The extremely high temperature agglomeration areas were primarily distributed in southern Chaoyang and central Fengtai, accounting for 35.96% and 13.58%, respectively, of the total extremely high temperature area. The cold island areas were concentrated in Haidian, Chaoyang, Tongzhou, Changping, and Daxing. The cold island proportion in Haidian, Chaoyang, and Changping decreased slightly, whereas that in Changping and Daxing significantly increased. The extremely low temperature agglomeration areas were concentrated in central Haidian, Fengtai, and Tongzhou, accounting for 32.69%, 25.01%, and 15.79%, respectively, of the total extremely low temperature area.
Owing to opposition between urban development and local ecological protection policies, the heat island area in 2019 expanded compared with that in 2014, with a mean surface temperature of 31.186 • C. In 2019, the highest surface temperature of the study area reached 49.34 • C, while the lowest reached 20.98 • C. The heat island areas were still largely concentrated in the Chaoyang, Fengtai, Haidian, and Daxing districts. The heat island proportion in Chaoyang decreased slightly while remaining relatively stable in the other districts. The extremely high temperature agglomeration areas were primarily distributed in central Fengtai, southern Chaoyang, eastern Daxing, and southern Tongzhou, accounting for 19.04%, 18.17%, 15.01%, and 11.11%, respectively, of the total extremely high temperature area. The cold island areas were concentrated in Chaoyang, Haidian, Changping, and Tongzhou. The cold island proportions increased slightly in Chaoyang and Changping, whereas they decreased slightly in Tongzhou, Daxing, and Haidian. The extremely low temperature agglomeration areas were concentrated in eastern Mentougou and western Haidian, accounting for 57.11% and 39.81%, respectively, of the total extremely low temperature area. Figure 7 shows that the spatial development axis of the heat island within the sixth ring road of Beijing from 2004 to 2019 remained in a northeast-southwest direction. However, the azimuth angle continued to increase, the flatness decreased, the directionality of the SDE decreased, and the dispersion increased ( Table 4). The heat island center was located at the junction of the Xicheng and Dongcheng districts next to Beihai Park. Owing to the 2008 Beijing Olympic Games and the frequent use of Beijing Capital International Airport as an international transportation hub (Figure 1), urban construction in the northeast of the study area (Chaoyang and Shunyi) developed rapidly, resulting in a 1.21-km northeastern shift of the heat island center from 2004 to 2009. Influenced by urban construction in the northwest region of the study area (Changping and Haidian), the center of the heat island shifted 1.03 km to the northwest from 2009 to 2014. Finally, the center of the heat island shifted 0.68 km to the southeast from 2014 to 2019, primarily owing to development of the Daxing and Tongzhou districts in the southeast region of the study area. Overall, the center of the heat island shifted 1.77 km to the northeast from 2004 to 2019. In summary, the change in the spatial pattern of the thermal environment was correlated with a change in urban construction intensity. The thermal environment effect expanded along a northeastsouthwest direction, and the difference in the thermal environment effect continued to decrease along the southeast-northwest and northeast-southwest directions. owing to development of the Daxing and Tongzhou districts in the southeast region of the study area. Overall, the center of the heat island shifted 1.77 km to the northeast from 2004 to 2019. In summary, the change in the spatial pattern of the thermal environment was correlated with a change in urban construction intensity. The thermal environment effect expanded along a northeast-southwest direction, and the difference in the thermal environment effect continued to decrease along the southeast-northwest and northeast-southwest directions.

Temporal Evolution of the Urban Thermal Environment
From Table 5, the heat island area exhibited a small and variable change over time. The total area of the heat island was 1523.  Over the study period, the mean heat island areas of Dongcheng and Xicheng exceeded 90% of the area of each district, with a significant heat island effect. The heat island area showed a trend of first decreasing and then continuously increasing, with a minor decreasing trend overall. From 2004 to 2019, the medium temperature area proportion decreased, whereas that of the high temperature area increased. Minor changes occurred in other temperature zones ( Table 6). In Fengtai and Daxing, the mean heat island area was >80% of the area of each district over the study period. From 2004 to 2019, the heat island area showed a trend of first continuously decreasing and then slightly increasing, with an overall significantly decreasing trend. From 2004 to 2019, the proportion of the medium temperature area to the total study area increased, and that of the relatively high and high temperature areas to the total study area decreased. Meanwhile, that of the relatively low temperature zone increased, and the other temperature zones exhibited minimal changes.  The mean heat island areas in Fangshan, Chaoyang, and Shijingshan accounted for >70% of the area of each district over the study period, indicating a significant heat island effect. The heat island area in Fangshan first increased, then significantly decreased, and was followed by an increase, with an overall significantly decreasing trend. From 2004 to 2019, the proportion of the medium temperature area increased, that of the relatively high and high temperature zones decreased, and that of the relatively low temperature zone increased. The other temperature zones exhibited minor changes. In Chaoyang and Shijingshan, the heat island area first decreased, then increased, and finally decreased again, with a small decrease overall. From 2004 to 2019, the proportion of the relatively high temperature area to the total study area decreased, that of the low temperature area in the Shijingshan district increased, and that of the other temperature zones exhibited minor changes. The mean heat island areas of Haidian, Tongzhou, Changping, and Shunyi all exceeded 50% over the study period, with a significant increasing trend overall. The heat island areas of Haidian, Tongzhou, and Shunyi first increased, then decreased, and finally increased again, whereas that of Changping showed a continuously increasing trend. The proportion of the relatively high and medium temperature areas to the total study area increased, that of the relatively low temperature zone decreased, and those of the other temperature zones all exhibited minor changes. The proportion of the mean heat island area of Mentougou to the total area of the district was <30% over the study period. Furthermore, the heat island area showed a trend of first continuously increasing and then decreasing, with an overall decreasing trend. The proportion of relatively high and medium temperature zones decreased, those of the low and extremely low temperature zones increased, and those of the other temperature zones all exhibited minor changes.

Correlation Analysis of Influencing Factors
In this study, 600 sample points were randomly selected in the study area using the sampling tool in the Data Management Tools module of ArcGIS 10.8 (Figure 8), and various factors were extracted in order to obtain the Pearson correlation coefficients of the LST and each index (Tables 7 and 8). Among them, the FVC, elevation, slope, and MNDWI were significantly negatively correlated with the LST, indicating that vegetation, terrain, and water can mitigate the UHI effect. The NDIBBI, road network density, NTL, POI density, NDISI, and land surface albedo were significantly positively correlated with the LST. As such, these indices reflect the effects of human activities and urban construction intensity on the urban thermal environment. Table 7. Index of urban thermal environmental impact factors.

First-Level Indicators Second-Level Indicators Third-Level Indicators
Natural factors

Spatial PCA
Four main factors were identified as affecting the urban thermal environment pattern in the study area (Tables 9 and 10): urban construction intensity, vegetation and water distribution, socioeconomic activities, and topography (index values of 0.525, 0.277, 1.245, and 0.572, respectively).
The first principal component reflects the influence of heat generated by urban construction on the thermal environment, including the building index, and the heat island spatial distribution was consistent with the urban built-up area. The second principal component reflects the impact of vegetation and water bodies on the thermal environment, including the water index and FVC. The third principal component reflects the impact of social and economic activities on the urban thermal environment, including the POI and road network density. The fourth principal component reflects the impact of the slope and elevation on the thermal environment. The topography of Beijing is high in the west and low in the east, and urban buildings are concentrated in low-altitude and flat areas, which is not conducive to thermal diffusion and significantly impacts the urban thermal environment.

Relationship between Principal Component Simulation and LST
Linear regression fitting analysis was conducted to test the fitting effect of the principal component factors on the urban thermal environment's spatial pattern. From Table 11, the regression equation for the LST of the study area, as well as the natural and abovementioned human factors, fit well. The F value of the regression equation was 288.84, and the significance coefficients of the first four principal components were all <1%. Therefore, these four principal components are all important factors affecting the urban thermal environment.
where X * 1 -X * 10 is the standardized variable value. For every 1 unit of change in the 10 influencing factors (NDISI, NDBBI, MNDWI, slope, DEM, FVC, albedo, NTL, POI, and road network density index), the LST of the research area changed by 0.068, 1.056, −0.128, 0.073, −0.205, −0.910, 0.371, 0.383, 0.28, and 0.34, respectively, according to Equation (14). Moreover, human factors intensified the urban thermal environment effect, among which NDIBBI had the most prominent positive effect, with a contribution of 1.056, and FVC had the most significant negative effect, with a contribution of 0.910. The NTL value, surface albedo, road network density, and POI density had significant positive effects. The elevation and MNDWI had significant negative effects, and the NDISI and slope had minor effects.
A simultaneous change in each of the above influencing factors will result in temperature changes of 2.498 and −1.316 • C and an overall temperature increase of 1.182 • C (i.e., human factors have a greater impact than natural factors on exacerbating the urban thermal environment effect).

Result Analysis
The results show that the spatial distributions of urban development and LST within the sixth ring road of Beijing from 2004 to 2019 were closely related, with high temperatures mainly occurring in the well-developed central part of the city and low temperatures mainly occurring in the less-developed urban fringe. In 2004, the heat islands were concentrated in the central and southwestern regions of the study area. In 2009, urban construction accelerated, and the northeastern region of the city was greatly developed, resulting in the heat island area expanding from southwest to northeast. However, by 2014, the heat island area had decreased, owing to the introduction of relevant ecological protection policies. In 2019, the heat island area had expanded, owing to the opposition between urban development and ecological protection policies. Some studies have pointed out that from 1994 to 2013, the urban scale of central Beijing became larger and construction land increased [62]. The surface temperature of nearly 60% of central Beijing increased by more than 4 • C [63]. Overall, the heat island area of the study area showed an increasing trend from 2004 to 2019. The thermal environment effect extended along a northeast-southwest direction, and differences along the southeast-northwest and northeast-southwest directions decreased.
From 2004 to 2019, the heat island area showed a small increasing trend, while spatial variation was obvious. The heat island area in the northern study area increased, while in the southern study area, it decreased. The proportion of the heat island area to the total area of each administrative region varied greatly between different periods, with a high proportion in the administrative regions located in the middle of the study area and a low proportion in the administrative regions located at the edge of the study area. A previous study pointed out that from 2003 to 2017, the high temperature area increased in Beijing, with the centroid concentrated in the center of the city. In contrast, the centroid of the low temperature area diffused to the surrounding areas [64], indicating that thermal environment effects in different regions of the city are significantly different.
The main factors affecting the urban thermal environment effect are urban construction intensity, vegetation and water bodies, socioeconomic activities, and geomorphology, in which the NDBBI had the most prominent positive impact (a contribution of 1.056 • C) and FVC had the most prominent negative impact (a contribution of −0.910 • C). Some scholars have pointed out that the spatial distribution of surface temperature is closely related to the NDVI [65], while the FVC and building density have the greatest impact on the surface temperature [66], which is consistent with the results of this study. Additionally, the NTL value, surface albedo, road network density, and POI density had prominent positive effects, elevation and MNDWI had prominent negative effects, and NDISI and slope had minor effects. Meanwhile, when all the influencing factors simultaneously changed by 1 unit, the temperature changed by 2.498 or −1.316 • C, with an overall temperature increase of 1.182 • C; that is, human factors had a greater influence than natural factors on aggravating the urban thermal environment effect.

Recommendations for Sustainable Development
The present study provides new insight into the causes of the urban thermal environment effect and provides guidance for sustainable urban planning. The results show that the spatial distribution of urban development and the LST within the sixth ring road of Beijing from 2004 to 2019 were closely related, and differences in the thermal environment effects between different administrative regions and different periods were obvious. The main factors affecting the pattern of the thermal environment were urban construction intensity, vegetation and water bodies, socioeconomic activities, and geomorphology. Therefore, when planning for different areas within cities, specific analyses should be conducted in conjunction with these factors in order determine corresponding strategies [67]. For example, in areas where industries are concentrated, optimizing neighborhood space and reducing waste emissions could facilitate a reduction in LST. However, in areas with dense buildings and few green spaces, changing the landscape composition by increasing urban green space and limiting impervious surface expansion is not an effective way to mitigate the increase in urban surface temperature. More attention should be paid to various costeffective means [68], such as road cooling [69][70][71], urban green roofs [72,73], the cooling of building materials [74,75], increasing vertical greening and vegetation density [76,77], and demolishing dilapidated buildings to increase air circulation. Stricter greening policies should be implemented in new expansion areas which, combined with the optimal park cooling range and scale, would help mitigate thermal environment effects [78]. Meanwhile, the results of this study indicated that increasing water system areas would effectively mitigate thermal environment effects. Water systems act as important cooling factors, and their effect on mitigating LST can be higher than that of vegetation. Thus, strengthening the management and optimization of water systems can play a positive role in mitigating thermal environment effects [79].

Limitations
This study involved a comprehensive analysis of the urban thermal environment effect and its influencing factors in Beijing; however, certain shortcomings remain. Owing to the complexity of surface temperature changes, limited availability of basic research data and other practical challenges, including the potentially inappropriate temporal resolution of the cross-section data used, were unavoidable. In this study, remote sensing image data of a few sunny days were selected. Although the meteorological conditions of the days on which the sample data were collected were similar and relatively good, the data from these 4 days were not fully representative of the 15-year period of surface temperature changes in the city. Therefore, the magnitude of the temperature sample data should be increased in future research to ensure model accuracy. Finally, the spatial resolution of the remote sensing data must be improved by using higher-resolution remote sensing data in the future.

Conclusions
The spatiotemporal evolution of urban thermal environment effects in Beijing (defined as the area within the sixth ring road) was analyzed by inversion of remote sensing data to obtain the LST in 2004, 2009, 2014, and 2019. In addition, based on multivariate spatial data (Landsat 4-5 TM, Landsat 8 OLI TIRS, DEM, road network, POI, and NPP/VIIRS nighttime light data), we obtained various natural and human factors closely related to the urban thermal environment effect and applied the standard deviation ellipse, spatial PCA, and other methods to analyze and identify the relationship between the thermal environment effect of urban human settlements and their influencing factors. The following conclusions can be summarized.
The spatial distribution of urban development and the LST within the sixth ring road of Beijing from 2004 to 2019 was closely related. High temperatures mainly occurred in the well-developed central part of the city, and low temperatures mainly occurred in the less-developed urban fringe. Thermal environment effects extended along a northeastsouthwest direction, and the differences between the southeast-northwest and northeastsouthwest directions continuously decreased.
From 2004 to 2019, the heat island areas showed a small increasing trend with obvious spatial variation. The heat island area in the northern part of the study area increased, and that in the southern part decreased. The proportion of the heat island areas to the total area of each administrative region varied greatly between different periods within the city, with high proportions in administrative regions located in the middle of the study area and low proportions in administrative regions located at the edge of the study area. The differences in the thermal environment effects between different administrative regions in different periods were obvious.
The main factors affecting the urban thermal environment effect were urban construction intensity, vegetation and water bodies, socioeconomic activities, and geomorphology, and the human factors had a greater impact than natural factors. Among them, the positive effect of the NDBBI was the most prominent, with a contribution of 1.056 • C. The negative effect of the FVC was the most prominent, with a contribution of −0.910 • C. The overall urban temperature increased by 1.182 • C under the influence of various influencing factors. Therefore, when planning for different areas within the city, specific analyses should be conducted in relation to the main factors affecting the thermal environment effect in order to determine the appropriate strategy. This study provides theoretical support for mitigating the urban thermal environment effect and promoting sustainable urban development.

Data Availability Statement:
The data presented in this study are available upon reasonable request from the corresponding author.