Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
An Integrated Methodology of Subjective Investigation for a Sustainable Indoor Built Environment. The Case Study of a University Campus in Italy
Next Article in Special Issue
Thermo-Hygrometric Variability on Waterfronts in Negative Radiation Balance: A Case Study of Balneário Camboriú/SC, Brazil
Previous Article in Journal
Short-Term Effects of “Polish Smog” on Cardiovascular Mortality in the Green Lungs of Poland: A Case-Crossover Study with 4,500,000 Person-Years (PL-PARTICLES Study)
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Spatiotemporal Characteristics of the Surface Urban Heat Island and Its Driving Factors Based on Local Climate Zones and Population in Beijing, China

1
College of Urban and Environmental Sciences, Peking University, Beijing 100871, China
2
Key Laboratory for Earth Surface Processes of The Ministry of Education, Peking University, Beijing 100871, China
*
Author to whom correspondence should be addressed.
Atmosphere 2021, 12(10), 1271; https://doi.org/10.3390/atmos12101271
Submission received: 23 August 2021 / Revised: 20 September 2021 / Accepted: 25 September 2021 / Published: 29 September 2021
(This article belongs to the Special Issue Cool Cities: Towards Sustainable and Healthy Urban Environments)

Abstract

:
The increasing degree of urbanization has continuously aggravated the surface urban heat island (sUHI) effect in China. To investigate the correlation between spatiotemporal changes of sUHI and urbanization in Beijing, land surface temperature in summer from 2000 to 2017 and the distribution of local climate zones (LCZs) in 2003, 2005, 2010, and 2017 was retrieved using remote sensing data and used to analyze the sUHI area and intensity change. The statistical method GeoDetector was utilized to investigate the explanatory ability of LCZs and population as the driving factors. The year of 2006 was identified as the main turning year for sUHI evolution. The variation the sUHI from 2000 showed first an increasing trend, and then a decreasing one. The sUHI pattern changed before and after 2009. Before 2009, the sUHI mainly increased in the suburbs, and then, the enhancement area moved to the central area. The sUHI intensity change under different LCZ conversion conditions showed that the LCZ conversion influences the sUHI intensity significantly. Based on population distribution data, we found that the relationship between population density and sUHI gets weaker with increasing population density. The result of GeoDetector indicated that the LCZ is the main factor influencing the sUHI, but population density is an important auxiliary factor. This research reveals the sUHI variation pattern in Beijing from 2000 and could help city managers plan thermally comfortable urban environments with a better understanding of the effect of urban spatial form and population density on sUHIs.

1. Introduction

Since the 1950s, the speed of global urbanization has increased significantly, accompanied by an increasing transformation of the natural landscape into human living areas [1]. By 2030, the global urban area is predicted to increase by 1.2 million km2. Most of this increase will be in Asia, especially in China and India, which will account for more than half of the new urban area [2,3]. Additionally, China’s urban population is forecast to double by 2030 compared with 2000 [4].
Urbanization leads to changes in land use and cover [5,6], mainly manifesting as the replacement of large natural surfaces with impermeable surfaces and buildings. In addition, with the constant large number of people moving to urban areas, building density and human activity intensity also increases [7,8]. Rapid urban expansion directly alters local climate and hydrological systems [9] due to a variety of urban environmental problems, among which the intensification of the urban heat island (UHI) has attracted wide attention. The UHI effect refers to a phenomenon where the temperature in urban areas is higher than the surrounding suburban areas [10,11]. Some studies have demonstrated that the UHI effect has been increasing with the development of urbanization [12,13,14,15].
The main causes of sUHIs are the heat storage of impermeable surface materials and the anthropogenic heat brought by the urbanization development [12,16,17,18]. These causes correspond to the two most significant processes of urbanization, namely land and population urbanization.
The effect of land urbanization on sUHIs is usually discussed based on land use and cover distribution maps. However, land urbanization changes the urban spatial form, and the coarse land surface classification commonly used for the urban area does not allow for an accurate analysis of the effects of urban spatial form on UHIs. In this context, a new urban land surface classification system called Local Climate Zones (LCZs), which can reflect urban spatial form, has been increasingly used in UHI research [19,20,21,22,23]. The LCZ was proposed by Stewart and Oke of the University of British Columbia [24,25]. It divides regional climate zones into several local climate zones according to different underlying surface types, which allows for the correlation between climate information and urban planning [26,27]. The LCZ classification system is consists of two main categories, built and land cover types. The two categories include 17 sub-categories (Table S1) based on building density, geometric form, and surface characteristics. The LCZ classification basis includes seven geometric and surface cover properties (sky view factor, aspect ratio, building surface fraction, impervious surface fraction, previous surface fraction, height of roughness elements, and terrain roughness class), based on which the urban spatial form is characterized. In population urbanization, one of the most remarkable characteristics is the change in urban population density [28]. Nevertheless, few studies discuss the influence of population density distribution on UHIs because of the lack of spatial population data at a fine scale.
Beijing is a megacity that represents the political, cultural, and commercial center of China. Since the reform and opening-up, Beijing has experienced a rapid urbanization process, which will continue for a long time. Due to this rapid urbanization, the UHI in Beijing has been increasingly intensified, which has created a variety of ecological and social problems that affect urban development and the residents’ lives [29,30,31]. Based on land surface temperature (LST) retrieved from remote sensing data, this research employs LCZ classification maps and spatial population density data to explore the effects of urbanization on the sUHI. The specific objectives of this study are to (1) reveal the spatiotemporal characteristics of the sUHI in Beijing from 2000 to 2017; (2) discuss the relationship between LCZs, population density, and the sUHI, and (3) estimate the effect of LCZ and population on the sUHI.

2. Materials and Methods

2.1. Study Area

Beijing (39°28′–41°05′ N, 115°25′–117°30′ E), the capital of China, is a megacity with a population of over 20 million people. It is located in the northern tip of the North China Plain, and it is surrounded by mountains on the west, north, and northwest borders. Beijing has a temperate continental monsoon climate with hot rainy summers and cold dry winters. With the rapid and intensive urbanization, a significant UHI effect has been noticed in Beijing [32,33,34].
The study area for this research is the region within the 6th Ring Road of Beijing (Figure 1), which is a primary urban area that has experienced a rapid urbanization process since the reform and opening-up. The area is 2269 km2, and more than half of the population in Beijing live there. This area has experienced significant urban sprawl and population growth [35] and is considered the Beijing central area.

2.2. Data

2.2.1. Landsat Data

Eighteen Landsat-5/7/8 images were downloaded from the United States Geological Survey (USGS) website (http://earthexplorer.usgs.gov/ (accessed on 29 October 2018)) at a spatial resolution of 30 m. All data were selected on a summer day with no cloud within the study area, from 2000 to 2017. The Landsat data were acquired at approximately 10:30 am (Beijing time) according to the transit time of the satellite. The dates of selected data are listed in Table 1.

2.2.2. High-Resolution Remote Sensing Data

Google Earth data with pixel size of 0.46 m × 0.46 m were used to provide the high-resolution remote sensing image required for LCZ classification (see Methods). The images were acquired from a business software named 91weituzhushou (Beijing Qianfanshijing Technology Co., Ltd., Beijing, China). Due to the limitation of data availability, 2003, 2005, 2010, and 2017 were selected as the observation years. The data were acquired in summer, on the dates listed for the Landsat data, as shown in Table 2.

2.2.3. Population Data

Population density is an important indicator to characterize domestic living conditions. However, traditional census data cannot reflect the population spatial distribution at a fine resolution due to the adopted statistical method, which is based on administrative divisions. In this research, we used the population data provided by the Global High Resolution Population Denominators Project (https://www.worldpop.org/ (accessed on 29 October 2018)), which includes annual global population distribution for individual countries from 2000 to 2020 at a high spatial resolution (approximately 100 m). Details of the population data used in this study are listed in Table 3.
The unit of population data was translated from population/pixel to population density (people/hm2) by Equation (1).
POP _ Density = ( POP _ Pix / ( 81 . 7 × 81 . 7 ) ) × 10000
where POP_Density is population density (people/hm2); POP_Pix is the pixel value of the original data (number of people in the pixel), and 81.7 × 81.7 is the area of each pixel (m2).

2.3. Methods

To assess the sUHI, the basic information, such as distribution of land surface temperature (LST) and LCZ was prepared by using remote sensing retrieval method and high-resolution remote sensing image classification. In addition, relative LST was calculated by LST. Based on the obtained distribution maps of LST, relative LST, and LCZ, the sUHI intensity was estimated in two ways and the trend analysis was applied using the two kinds of sUHI intensity. Further analysis explored the explanatory power of two driving factors, the LCZ and population, on the spatio-temporal pattern of the sUHI using a recently proposed statistical method called the GeoDetector.

2.3.1. Local Climate Zone (LCZ) Classification

The LCZ classification system was used to study the influence of urban spatial form on the sUHI. Most existing research on sUHIs using LCZs is based on a single LCZ classification map. Therefore, they cannot reflect the spatiotemporal changes of LCZ distribution and their influence on the sUHI from a dynamic perspective. Consequently, there are few studies at the multi-year scale, especially over long periods (>10 years). For the present study, four LCZ classification maps were retrieved and used to understand the relationship between urban space form and the sUHI for an 18-year period from 2000 to 2017.
The World Urban Database and Access Portal Tools (WUDAPT) provide an instruction for LCZ classification with open access data and open-source software [36]. The LCZ system consists of 17 standard types, which include 10 built (LCZ 1–10) and seven land cover (LCZ A–G) types (Table S1). The LCZ areas 6, 7, 9, and 10, and LCZ C were rare in the study area. Therefore, they were excluded from our LCZ classification, thus resulting in 12 kinds of LCZs based on the Landsat and Google Earth images. Considering the lack of sub-meter high-resolution remote sensing images before 2003 within the 6th Ring Road in Beijing and the study period (2000–2017), images from 2003, 2005, 2010, and 2017 were selected for the LCZ classification mapping. For that, the ENVI 5.3 (Harris Geospatial Solutions, Boulder, Colorado) and SAGA GIS [37] software was used. The specific steps are as follows:
  • Data preprocessing
    The data preprocessing included coordinate system transformation, image registration, resampling, and clipping. Landsat data with 30 m resolution were resampled to 60, 90, 120, 150, 180, 210, and 240 m to compare the classification accuracy under different resolutions.
  • Selection of training and verification samples
    The ENVI 5.3 was used instead of Google Earth, which is recommended by WUDAPT, to select the training and verification samples using a visual interpretation method based on sub-meter level Google Earth images. Approximately 20–30 training samples and a similar number of verification samples were selected for each LCZ type.
  • LCZ classification using random forest classification method
    Based on the training samples and Landsat images, the LCZ classification results were obtained using the Random Forest Classification tool in SAGA GIS.
  • Accuracy evaluation and reclassification
    The confusion matrix, Kappa coefficient, and overall accuracy were obtained using the Confusion Matrix (polygons/grid) tool of SAGA GIS to estimate the LCZ classification accuracy. When the Kappa coefficient and overall accuracy were both larger than 0.7, the LCZ classification results were acceptable. Otherwise, the training and verification samples were corrected and reclassified for a new LCZ classification. This cycle continued until the LCZ classification met the accuracy requirement.

2.3.2. Land Surface Temperature (LST) Retrieved by a Mono-Window Algorithm

LST was retrieved from remote sensing data, in order to calculate the sUHI intensity. Among various remote sensing retrieval methods, mono-window algorithms proposed by Qin et al., (2001) [38] have been widely used in LST retrieval using Landsat images and can produce highly accurate results. Although the Landsat-8 Thermal Infrared Sensor (TIRS) has two thermal infrared bands, which allow the use of a split-window algorithm, a mono-window algorithm was adopted to ensure LST consistency over time. The LST was calculated by Equation (2),
T s = 1 C a ( 1 C D ) + b ( 1 C D ) + C + D T sensor DT a
where Ts (K) is the LST (K); Tsensor is the sensor’s brightness temperature (K); Ta is the average atmospheric temperature (K); a and b are regression coefficients, which have different values for different sensors and temperature ranges (Table S2); and C and D are intermediate variables, which are calculated based on emissivity (ε) and atmospheric transmissivity (τ) by Equations (3) and (4). The detailed calculation procedure of Ta, ε, and τ is in supplementary materials (Table S6, Equation (S4), Tables S3–S5).
C = ε τ
D = ( 1 τ ) 1 + ( 1 ε ) τ

2.3.3. Relative LST (LSTR)

To produce comparable LST within the study period, LST was normalized and relative LST (LSTR) was calculated based on the LST distribution map. LSTR reflected the difference between LST at the location of each pixel and the mean LST of the study area, as shown in Equation (5) [39],
LST R = Δ LST / LST m = ( LST i LST m ) / LST m
where LSTi is the LST value at the location of pixel i, and LSTm is the mean LST within the study area. The value of LSTR ranges between −1 and 1, and reflects the sUHI intensity at the location of the pixel. Based on LSTR, the sUHI intensity can be divided into five grades (Table 4). According to the pixel number variation, the trend of the sUHI area can be identified. Moreover, the pixel number variation of different sUHI grades can reflect the sUHI intensity change.

2.3.4. sUHI Intensity Based on LCZ Classification

Stewart and Oke (2012) defined UHI intensity as the difference between the mean temperature of the built types in LCZ with the largest proportion of the built area and the mean temperature of LCZ D, which is seen as the main LCZ type in the suburbs [24]. Accordingly, the sUHI intensity based on LCZ classification can be calculated by Equation (6),
sUHII LCZ = LST LCZ   X LST LCZ   D
where sUHIILCZ is the sUHI intensity based on the LCZ classification; LSTLCZ X is the mean LST of LCZ X, which occupies the largest area proportion among built types in the urban region; and LSTLCZ D is the mean LST of LCZ D.
The statistic results of our LCZ classification maps show that LCZ B and D were the main LCZ types in the suburbs during the summers 2003 and 2005 and of 2010 and 2017, respectively. In addition, LCZ 4 and 5 were the main LCZ types in the built area in the summers of 2003 and 2017 and of 2005 and 2010, respectively. Therefore, the sUHIILCZ in this research can be estimated by the following corrected equation:
sUHII LCZ = LST LCZ   4 LST LCZ   B LST LCZ   5 LST LCZ   B LST LCZ   5 LST LCZ   D LST LCZ   4 LST LCZ   D
Equation (7) represents a common method based on usual definitions. However, different equations used to estimate sUHIILCZ in different years do not ensure the comparability of the sUHIILCZ time series. Therefore, we proposed an area-weighted method to estimate the sUHIILCZ using a uniform equation:
AW _ sUHII LCZ = LST LCZ   4 & 5 LST LCZ   B & D
where AW_sUHIILCZ is the area-weighted sUHIILCZ; LSTLCZ 4&5 is the area-weighted mean LST of LCZ 4 and 5; and LSTLCZ B&D is the area-weighted mean LST of LCZ B and D.
To obtain the sUHIILCZ of the study area in the summers of 2000–2017, the corresponding LCZ classification maps of each year are needed. However, limited by the availability of high-resolution images and considering that the degree of urban form change is not large in adjacent years, the maps were not produced for all years. The LCZ classification maps of 2003, 2005, 2010, and 2017 were obtained and used to represent the LCZ distribution of two to three adjacent years to match the corresponding LST distribution maps (Table S7).
In this study, the annual LST distribution maps and the four-period LCZ classification maps were used to make the curves of sUHI intensity. It is worth noting that although the sUHI intensity for each year was obtained in this study, the sUHI effects of adjacent years are not comparable because of the different dates of Landsat images used for each year, and the sUHI intensity of each period only reflects the sUHI situation in the summer of that year approximately. Therefore, in order to obtain the approximate trends and major turning points of sUHI intensity in summer, a 5-year moving average method was used to deal with the annual series of sUHI intensity in the study. For the sUHI effect in a given region, because it takes 5 to 10 years for the population and land use to change significantly, significant changes often take a relatively long time. Additionally, the reason why longer time intervals such as five-year intervals were not directly used to analyze the trend is that the image acquisition dates for 2000, 2005, 2010, and 2015 were also mismatched and using only these four years of data to analyze the trend could easily be overly influenced by data fluctuation due to the short series. Therefore, the study used a 5-year moving average approach to annual sUHI intensity to reduce the impact of data fluctuation on the trend.

2.3.5. Quantitative Estimation of Explanatory Ability of Driving Factors

To estimate the correlation between LCZ, population, and the sUHI, GeoDetector was used due to its advantages in the quantitative analysis of spatial stratification heterogeneity. GeoDetector, which was proposed by Wang et al., (2010) [40], can be used to address categorical variables such as LCZ. Moreover, population density can be discretized to fit the data demand. The basic principle of GeoDetector is that if there is significant spatial consistency between an independent and a dependent variable, there is a correlation between them. The geographical detector consists of four sub-detectors: factor, risk area, interaction, and ecological detectors. The factor and interaction detectors were the primary detectors used in this study.
The factor detector is used to calculate the degree of interpretation of different independent variables for the dependent variable. The power of interpretation is measured by the q value, and the equation is as follows [41]:
q = 1 h = 1 L N h σ h 2 N σ 2
where h = 1, 2, L is the layer of the independent or dependent variable; Nh and N are the number of samples in layer h and the total number of samples, respectively; σ h 2 and σ2 are the variance of value in layer h and in the entire region, respectively. The value range of q is [0,1]. The closer the value is to 1, the stronger the power of interpretation for the independent variable to the dependent variable, and the value itself represents the potential of interpretation (the independent variable can interpret 100 × q% of the dependent variable).
The interaction detector can detect the interaction between bivariate variables, namely X1 and X2. By comparing the q value of X1, X2, and X1∩X2, which indicates a new layer produced by overlaying X1 and X2, the interaction type can be identified (Table 5) [26]. As mentioned above, the q value can be used to measure the explanatory power of the independent variables on the dependent variable. The explanatory power of the independent variables X1, X2 and X1∩X2, the variable reflecting their interaction, on the dependent variable are expressed as q(X1), q(X2) and q(X1∩X2), respectively. By comparing q(X1∩X2) with q(X1) and q(X2), it can be judged whether the explanatory power of the interaction of X1 and X2 on the dependent variable is enhanced or diminished with respect to a single specific variable or both variables. Further, it can be identified whether the explanatory power of X1 and X2 on the dependent variable is independent, i.e., if there is a linear relationship between q(X1), q(X2), and q(X1∩X2). Specifically, when q(X1∩X2) is less than q(X1) or q(X2), it means the interaction of X1 and X2 weakens the explanatory power compared to that of X1 or X2 alone. If q(X1∩X2) is between q(X1) and q(X2), the interaction of X1 and X2 shows more explanatory power than considering one of the two variables, but it weakens when considering the other factor. When q(X1∩X2) is greater than both q(X1) and q(X2), the interaction of X1 and X2 enhances the explanatory power compared to considering either variable only, and the q(X1∩X2) can be further compared to the sum of q(X1) and q(X2). It represents that it is independent between the explanatory contribution of X1 and that of X2 when q(X1∩X2) is equal to the sum of q(X1) and q(X2). Moreover, when q(X1∩X2) is greater than the sum of q(X1) and q(X2), the explanatory power of the interaction of X1 and X2 not only includes that brought by X1 and X2 itself, but also an additional nonlinear part.
In this research, the dependent variable in the model was the LSTR, which represented sUHI intensity, and the independent variables were the LCZ types and population density. Five thousand random points were used to extract the values of these three variables, and four groups of data (extracted from 2003, 2005, 2010, and 2017) were constructed. To simplify the data structure, population density data of more than 400 people/hm2 were divided into one class.
In addition to the section data analysis, the influence of LCZ and population density changes on the sUHI intensity was analyzed by using a geographic detector. This analysis was performed to reflect the dynamic relationship between LCZ, population density, and sUHI intensity. The dependent variable was the Theil–Sen regression slope of the LSTR time series, and the independent variables were the LCZ change types and Theil–Sen regression slope of population density time series. According to the sUHI change trend presented in this study, three time periods, 2000–2009, 2009–2017, and 2000–2017, were selected as study intervals. The LCZ types of 2000 and 2009 were replaced by those of 2003 and 2010, respectively, and the Theil–Sen regression slope of population density was also discretized.

3. Results

3.1. LST Retrieved from Landsat Images

Eighteen LST distribution maps of the summers from 2000 to 2017 were obtained based on Landsat data (Figure 2). The unit of LST was converted from Kelvin (K) to centigrade (°C), and the statistic for each LST distribution map is listed in Table S8. Due to the cloud disturbance in 2001 and 2013, the LST values of these two years were replaced by the average LST of two adjacent years.

3.2. LCZ Classification Maps

The classification accuracy of four LCZ classification maps for 2003, 2005, 2010, and 2017 at different resolutions were compared, and the results show that the classification accuracy performed well at a 60 m spatial resolution for all LCZ classification maps (Figure 3). Therefore, 60 m was selected as the spatial resolution of LCZ classification maps in this research.
Four LCZ classification maps (Figure 4) were produced after repeated iterative selection of training and verification samples. The classification accuracy was estimated by a confusion matrix (Tables S9–S12). Figure 4 shows clear LCZ distribution patterns. The LCZ types of typical regions were identified relatively accurately. For instance, the dense trees (LCZ A) of the mountainous region, low plants (LCZ D) of the suburbs, and water bodies (LCZ G) were accurately identified. In addition, the LCZ change in some regions was also clearly reflected. For example, for the terminal T3 in the Beijing Capital International Airport, the four LCZ classification maps clearly reflect the unconstructed, under construction, and completed construction stages.

3.3. Spatiotemporal Variation of sUHI

The pixel numbers of heat island (LSTR ≥ 0) and cool island (LSTR < 0) were identified for 2000–2017 (Figure 5). Furthermore, the number of heat islands at different intensity grades was also identified (Figure 6). To reduce the volatility in raw data and better identify the variation trend, a data smoothing method of moving average and linear regression for 5 years was applied.
From 2000 to 2017, the variation trend of the pixel number of the heat island, namely the total sUHI area, in the 6th Ring Road of Beijing increased first and then decreased, and there was an obvious turning point in 2006. In the limited region, the sum of heat and cold island pixels was fixed, so their variation trends were opposite.
Based on the statistic results of pixel numbers of different sUHI intensities, the proportional area of weak sUHI increased from 20.1% in 2000 to 31.1% in 2017. In contrast, the proportional area of medium and strong sUHIs decreased from 20 and 5.9% to 16.6 and 3.4%, respectively, for that same period.
From 2000 to 2017, the area of weak sUHI increased first and then decreased slightly. The turning point of the trend occurred in 2009. The variation in the trend of medium sUHI was similar to that of strong sUHI. Both showed a decrease from 2000 to 2009, after which they remained stable. The number of extreme sUHI pixels was significantly small (less than 2500), and the variation trend of pixels fluctuated.
To reflect the spatiotemporal variation characteristics of the sUHI, the slope of LSTR time series at the location of each pixel in three different period was calculated by the Theil–Sen regression method (Figure 7). A positive slope indicates an increasing sUHI, and a negative slope indicates a decreasing sUHI.
From 2000 to 2017, the most significant area of sUHI intensity enhancement was outside the 5th Ring Road, and the sUHI intensity in the region of the 3rd Ring Road increased slightly. In this period, a large area presented a sUHI intensity enhancement trend between the 3rd and 5th Ring Roads.
The sUHI intensity variation pattern showed a difference between the periods of 2000–2009 and 2009–2017. In the first period, the sUHI intensity increased in the U-shaped area of the westward opening outside the 4th Ring Road, and decreased in most areas within the 4th Ring Road. In the second period, the sUHI intensity within the 4th Ring Road increased significantly, and decreased in a large area outside the 4th Ring Road. In addition, the sUHI intensity in the southwest of the study area also increased from 2009 to 2017.
Furthermore, according to the change trend of sUHI intensity for each pixel in the three study periods, the pixels were divided into four classes (Figure 8), namely continuous increase (↗), first increase and then decrease (↗↘), first decrease and then increase (↘↗), and continuous decrease (↘). The main type within the 3rd Ring Road was first decrease and then increase, and since 2009, the sUHI intensity in this area increased significantly, which demands the attention of city managers. In recent years, the areas with a rising trend of sUHI intensity were concentrated in the southwest, southeast, and northwest. The areas with ring-shaped surroundings outside the 4th Ring Road showed an increasing trend in the early period and a decreasing or continuous decreasing trend in recent years, which indicates that the sUHI intensity near the central urban area has recently been alleviated.

3.4. sUHI Intensity Estimation Based on LCZ Scheme

Using Equations (7) and (8), the sUHI intensity based on the common and area-weighted methods was estimated from 2000 to 2017, as shown in Figure 9. The sUHI intensity was in the range of 2–6 °C. The results of the area-weighted method indicate that the minimum and maximum sUHI intensities occurred in 2014 and 2006, respectively. Two obvious turning points at approximately 2006 and 2014 or 2015 were observed.

3.5. Effect of LCZ and Population on sUHI Intensity

3.5.1. Thermal Characteristics of LCZ

According to the four LCZ classification maps of 2003, 2005, 2010, and 2017, the LSTR of each LCZ-type areas were extracted for statistical analysis, and corresponding violin charts were obtained (Figure 10). The results show that in 2005, 2010, and 2017, the rank of sUHI intensity in the compact zones was LCZ 1 < LCZ 2 < LCZ 3. However, in 2003, it was LCZ 1 < LCZ 3 < LCZ 2. The sUHI intensity rank in open zones (including LCZ 8) was consistent, being LCZ 4 < LCZ 5 < LCZ 8 in all four years.
The vegetation areas all presented the cold island effect. In 2005, 2010, and 2017, the intensity rank was LCZ A < LCZ D < LCZ B, whereas in 2003, it was LCZ D < LCZ A < LCZ B. The LSTR of LCZ E was stable at approximately 0.1. In 2005, the data distribution of the LSTR of LCZ F presented an irregular bimodal pattern, whereas in the other three years, the value of LSTR was approximately 0, thereby being smaller than that of LCZ E. The median LSTR value of LCZ G (water) was less than -0.1, thus being smaller than those of LCZ B and D. Nevertheless, its numerical relationship with LCZ A is uncertain.

3.5.2. LSTR Change with LCZ Conversion

The multi-phase LCZ classification maps provided a dynamic perspective for studying the influence of LCZ conversion on sUHI. The LCZ conversions from 2003 to 2010, from 2010 to 2017, and from 2003 to 2017 were obtained using the LCZ classification maps of 2003, 2010, and 2017, respectively. Combined with the slope of the LSTR time series in 2000–2009, 2009–2017, and 2000–2017, the relationship between LCZ conversion and sUHI intensity change was estimated (Figure 11).
The LCZ conversion can be divided into four categories: (1) mutual conversion between built types (LCZ 1–8 to LCZ 1–8), (2) built types to natural cover (LCZ 1–8 to LCZ A–G), (3) natural cover to built types (LCZ A–G to LCZ 1–8), and (4) mutual conversion between natural covers (LCZ A–G to LCZ 1–8).
For the first category, there were similarities and differences between the sUHI intensity changes of 2000–2010 and 2010–2017. One common characteristic was that when other built types were converted into LCZ 4, or when LCZ 8 was converted into other built types, the sUHI intensity decreased. The sUHI intensity also decreased when LCZ 3 was converted to LCZ 1/2, and LCZ 2 translated to LCZ 1. When LCZ 4 was converted to LCZ 2/3, the sUHI intensity increased. The most prominent difference was that when other built types were converted to LCZ 5, the sUHI intensity increased in 2000–2017 but decrease in 2000–2010.
Compared with the first category of LCZ conversion, the other three categories of LCZ conversion presented a similar effect on the sUHI intensity change from 2000–2010 and from 2010–2017. The second category led to a decrease in sUHI intensity, and only when built type was converted to bare soil (LCZ E), was there a difference between the sUHI intensity change in the two periods. The sUHI intensity increased in the third category, except for some individual conversion types. In the fourth category, the sUHI intensity increased when vegetation was converted to other types, and the opposite also occurred (sUHI intensity decreased when other types were converted to vegetation).

3.5.3. Relationship between LSTR and Population Density

To explore the relationship between human activities and sUHI, the population density and corresponding LSTR values were extracted at 5000 location points, which were randomly selected within the study area (excluding water bodies) for 2003, 2005, 2010, and 2017 (Figure 12).
The scatter plots of LSTR and population density for the four years showed a similar data distribution characteristic. Overall, LSTR increased with increasing population density. The median curve of LSTR had two turning points with increasing population density. As the population density increased from 0, the LSTR increased sharply until the first turning point. Between the first and second turning points, LSTR increased with increasing population density at a lower speed compared to the previous period. After the second turning point, LSTR no longer changed significantly.
The population density values of the two turning points are shown in Figure 12. Based on the statistical results of the four years, the population density values of the first and second turning points were approximately 70 and 300 people/hm2, respectively.
To detect the regions where population density had a significant direct impact on sUHI intensity, the Spearman rank correlation coefficient between LSTR and population density time series was calculated from 2000 to 2017 at the position of each pixel. The spatial distribution of correlation coefficients and p-values are shown in Figure 13 (only the area of p < 0.05 is shown).
The areas where the change in population density had a significant impact on sUHI intensity change were mostly distributed outside the 4th Ring Road, especially the outer suburbs outside the 5th Ring Road. Combined with the population density distribution map (Figure S1), these areas were distributed within low and medium population density areas.

3.5.4. Estimation of the Effect of LCZ and Population Density on sUHI Intensity Based on GeoDetector

For factor detection, the p values corresponding to the q values of the driving factors of LCZ type and population density were both less than 0.001 (Table 6), which indicates that the results passed the significance test. The q values indicate that LCZ can explain approximately 60–75% of the sUHI spatial differentiation, whereas population density can explain approximately 10–30% (Table 6). The interaction detection results (Table 6) showed that when considering the influence of both LCZ and population density on sUHI intensity, the explanation of the interpretation reached more than 70% (except for the results in 2005), and their interactions were enhanced and bivariate.
As for the effects of LCZ conversion and population density on sUHI intensity, the results of factor and interaction detection are listed in Table 7. The p value for the q values of each factor is less than 0.01 for all years, thus indicating that the results passed the significance test. The q value of LCZ conversion types was between 0.3 and 0.4, whereas that of population density change was between 0.01 and 0.02. The q value of the interaction detection results was greater than the sum of the q values of the two factors, which demonstrates an enhanced and nonlinear interaction.

4. Discussion

4.1. Characteristics of sUHI Spatiotemporal Variation

The LSTR and other similar LST normalization methods are widely used to analyze the sUHI effect [42,43,44,45]. Based on the pixel number trend of sUHI intensity grades, the characteristics of sUHI areas and intensity development in the area within the 6th Ring Road area of Beijing in 2000–2017 were investigated. The results showed that there were different stages in the sUHI development. In order to get the actual physical value of sUHI intensity and verify the results based on LSTR, this study also estimated the sUHI intensity based on the LCZ scheme, which was calculated as 2–6 K.
In terms of the evolution of sUHI intensity grades based on LSTR, for the sUHI area, 2006 was the turning point. From 2000 to 2006, the sUHI area increased, but it decreased after 2006. However, the turning point of sUHI intensity was in 2009. The sUHI intensity increased from 2000 to 2009, after which it remained stable. Nevertheless, the turning point of 2006, which was different from the LSTR scheme, was observed in the time series of sUHI intensity based on the LCZ scheme. The decreasing trend of sUHI intensity from 2006 to 2009 was consistent with that calculated based on the LSTR method. However, there was still a declining trend of sUHI intensity based on the LCZ scheme from 2009 to 2015. The different results indicate that the change in sUHI intensity might get lost in the division of sUHI intensity grade and therefore cannot accurately reflect the change in sUHI intensity. The rising trend of sUHI intensity in 2000–2006 was contrary to that of the estimation based on the LSTR method. To understand this, the pixel numbers of the medium and strong sUHI in 2000–2006 were further counted. The overall trend of the sum of the two sUHI grades increased, which is consistent with the trend of the LCZ-based sUHI intensity variation for the same period. Therefore, it can be concluded that when the LSTR method was used to obtain the medium and strong sUHI sequences, part of the increasing signal was masked by the overall downward trend in the average process.
According to the turning points of sUHI area and intensity, the development of sUHI in Beijing can be divided into three stages. (1) 2000–2006: both area and intensity of sUHI increased; (2) 2006–2015: both area and intensity of sUHI decreased; (3) 2015–2017: sUHI intensity increased (which was uncertain because of the short period), but the area decreased. The overall trend of an initial increase in sUHI intensity followed by a decrease presented in our study is consistent with other research on sUHI intensity variation in Beijing. Liu et al., (2020) investigated the magnitude and spatial patterns of sUHI from 2003 to 2018 in a similar area of Beijing [46]. The turning year of 2009, before which the sUHI intensity increased and after which it decreased, was also found in their study. Meng et al., (2018) also found an upward trend of sUHI intensity from 2003 to 2008 [47]. Combining all the studies above with ours, the sUHI intensity in Beijing increased from 2000 and then presents a downward trend. The turning year may be around 2006–2009. Furthermore, the sUHI intensity increased from 2016 to 2018 in the study by Liu et al., (2020) [46], which is similar to our finding. Therefore, the sUHI intensity may present upward trend recently, after 2015 or 2016, which needs further research.
For the full research period (2000–2017), the spatiotemporal pattern of sUHI change shows that areas with increased sUHI intensity were mainly concentrated in a parabola belt area of the westward opening outside the 4th Ring Road, where significant urbanization occurred [48], especially in the southeast area. In addition, the area in the southwest corner and the Capital International Airport Area in the northeast corner were also areas with significantly increased sUHI, while the sUHI intensity in the urban area within the 4th Ring Road showed a downward trend, with a slight increase in some areas. Since 2000, the sUHI intensity in the central area of Beijing decreased, or increased insignificantly, but it increased in peripheral areas. The spatial distribution of the Theil–Sen regression slope of LSTR in different periods reflected the spatiotemporal heterogeneity of sUHI intensity change. There was a general trend for sUHI enhancement in space, which was presented as the sUHI intensity increase from north to south, from east to west, and from periphery to center. The sUHI pattern changed before and after 2009. Before 2009, the sUHI mainly increased in the suburbs, developing from the inner to outer suburbs, and it decreased or slightly increased in the central urban area. After 2009, the sUHI intensity enhancement area moved to the central urban area, and sUHI decreased in a large area outside the 4th Ring Road. Nevertheless, the spatial patterns of sUHI change identified by Liu et al., (2020) [46] are opposite from ours. The main reason may be the difference in the definition of sUHI intensity used in their study and ours, which is based on the LST difference between urban and rural area identified by the proportion of impervious surface, and LSTR grade, respectively. Moreover, a similar increasing trend is identified in the rural area, similar to the area detected in our study in the research of Li et al., (2020) [49]. The research investigates the spatial pattern variation of Beijing from 2000 to 2013 using the difference between the LST of each pixel and the averaged of background areas as the sUHI intensity, which is similar to our LSTR based sUHI intensity definition. Therefore, the different method of calculating sUHI may lead to different or even opposite results for sUHI spatiotemporal patterns.

4.2. Relationship between LCZ, Population Density, and sUHI

The thermal differential of LCZs has been widely recognized [50,51,52,53,54]. According to the analysis of the relationship between LCZ type and LSTR, there are significant differences in the sUHI intensity of different urban forms. The order of sUHI intensity of built types showed that the sUHI intensity was closely related to the height of buildings. In general, the lower the building height, the stronger the sUHI intensity, which is consistent with previous studies [54,55,56]. For the same height level, the sUHI intensity in open zones was significantly lower than those in compact zones. For different heights, the difference of sUHI intensities in open zones was larger, which demonstrates that the change in sUHI intensity according to height is more obvious in an open environment.
The LCZ classification maps for different years allowed for analysis of the relationship between the change in urban spatial form and sUHI change. Wang et al., (2019) found the warming trend in the Pearl River Delta from 2000 to 2015 is related to conversion between and within LCZ types [57]. Based on the sUHI intensity change caused by the LCZ conversion in 2000–2009, 2009–2017, and the overall performance in 2000–2017, specific LCZ changes would lead to a stable sUHI intensity change trend. If other built types were transformed into LCZ 4, or LCZ 3 was transformed into other built types (except LCZ 8), the sUHI intensity would be effectively reduced, and the sUHI intensity would significantly increase if other LCZ types were converted to LCZ 8. Research by Wang et al., (2019) also shows that the expansion of LCZ 8 lead to LST increase [57]. An instability of sUHI intensity change caused by the LCZ transformation was also observed, especially in the transformation between built types. For example, the transformation of LCZ 2 into other built types (except LCZ 4) led to the decrease in sUHI intensity in 2000–2009 and the increase in 2009–2017.
Population has close relationship with UHI [58,59]. Zhang et al., (2013) used multiple linear regression to investigate the relationship between population density, which is extracted from statistical yearbooks, and other factors, and UHI in Shanghai [60]. A positive effect of population density on UHI was found in this study. However, Du et al., (2016) found that there was no significant correlation between population density and UHI intensity in the Yangtze River Delta Urban Agglomeration based on population data from statistical yearbooks and the linear regression method [61]. The statistical population data and the simple regression method may have led to the inconsistent results. In our study, based on the relationship between LSTR and the population density distribution map, increasing population aggravated the sUHI effect, and the impact varied with different population density ranges and regions. The influence was most significant when the population density was at a low level and began to increase. In that stage, the sUHI intensity increased sharply, but the sUHI was still weak (LSTR < 0.1). With the increase in population density, the direct impact of human activities on sUHI tended to be small, and the sUHI intensity approached or reached medium sUHI. After reaching a certain threshold (approximately 300 people/hm2), the increase in population density no longer increased the sUHI intensity effectively. In terms of space, the relationship between sUHI intensity and human activities was the most significant in the suburbs outside the 4th Ring Road, with a high positive correlation, but it was not significant in the central urban area.
The change in population density affected the sUHI differently according to different stages and regions. The reason for this may be that during the increase in population density from extremely low to the first turning point, not only the artificial heat increased, but the underlying surface also changed dramatically as the vegetation was converted to construction land [35]. These factors induced the LST to increase sharply, thus leading to the sUHI effect. From the first to the second turning point, urbanization occurred continuously. With the increase in population density, anthropogenic heat emissions, such as heat emissions from air conditioning in summer and transportation, gradually increase, which aggravates the sUHI effect [62]. However, when the population density reached the second turning point, the load of energy consumption reached its peak, and the anthropogenic heat emissions do not increase as they are limited by the urban carrying capacity. In addition, the process of increasing population density is often accompanied by the process of land urbanization, with a consequent increase in the area of building surfaces and impervious surfaces, especially in the early urbanization process, i.e., the early stage of rapid population density increase. During this stage, a large amount of natural surface, which has a large specific heat capacity and can remove heat from the surface through evaporation or transpiration, is replaced by building surfaces and impervious surfaces. In contrast, the lower specific heat capacity and lack of evaporation and transpiration on building surfaces and impervious surfaces can exacerbate the sUHI. However, as land urbanization reaches a mature stage, the increase in population will no longer significantly increase building surfaces and impervious surfaces due to the limitations of plot ratio and land development, which means that the area and ratio of building surfaces and impervious surfaces will reach a relatively stable condition. Consequently, its enhancement of sUHI will also reach a relatively stable condition. That is, when the population density reaches a certain value, its effect of enhancing the sUHI effect by influencing the increase in building surfaces and impervious surfaces will not change significantly after the population density reaches a certain value. Therefore, the increase in population density no longer affects the sUHI intensity significantly.
Recently, GeoDetector has increasingly been used in the analysis of driving factors of UHI [63,64,65], indicating that GeoDetector is appropriate for driving factor analysis for UHI. The quantitative method of GeoDetector helps in revealing the interpretability of LCZ, population, and their interaction to sUHI. Our driving factors analysis based on GeoDetector showed that both LCZ type and population density can effectively explain the spatial differentiation of sUHI, but the explanatory ability of LCZ type was significantly higher than that of population density. The interaction effect of the two factors can significantly improve the explanatory ability. Compared to population density, LCZ, which characterizes the urban space form, was the main factor influencing sUHI. Population density is an important factor in human activities and has a significant impact on sUHI, but its ability to explain the spatial differentiation of sUHI independently was poor. The interaction effect of the two factors was presented as enhancement and bivariate, thus indicating that the combination of urban space form and human activities can enhance the ability to explain the spatial differentiation of sUHI.
The dynamic analysis results show that LCZ conversion can contribute 30–40% of the spatial structure of sUHI intensity change, and the change in population density had a weak ability to explain the spatial structure of sUHI intensity change. This result shows that, compared to population density, LCZ conversion was the key factor leading to the sUHI intensity change. The phenomenon of urban spatial form being the main influencing factor, rather than the change in population density, has also been found in the center area of Kuala Lumpur [66]. Nevertheless, the results of interaction detection showed that when considering both LCZ conversion and change in population density, the q values of the three time periods were significantly improved. The explanatory power was about 10% higher than that of LCZ variation only, and significantly higher than that of population density variation only. The results indicate that despite the weak ability of population density change to interpret sUHI change independently, population density was indispensable to explain the changes in the sUHI effect.

5. Conclusions

From 2000 to 2017, the change in sUHI area and intensity within the 6th Ring Road in Beijing has an obvious turning year of 2006 and spatiotemporal heterogeneity. The sUHI area and intensity increased from 2000 to 2006 and then decreased. The upward trend of sUHI intensity after 2015 needs further confirmed. The sUHI development mode changed in 2009. Before 2009, the sUHI mainly increased in the suburbs, and developed from the suburbs to the outer suburbs, and the sUHI intensity in the central urban area decreased or increased slightly. After 2009, the increased area of sUHI moved to the central urban area, and a large sUHI declining area appeared outside the 4th Ring Road. Moreover, the sUHI intensity was quantitatively estimated, and ranged from 2 to 6 K. LCZ had an explanatory power in relation to the sUHI effect stronger than population density. Nevertheless, the combination of these two factors improved the overall explanatory power. From a dynamic perspective, LCZ conversion was the main factor leading to the sUHI intensity change. The ability of population density to explain the changes in the spatial structure of sUHI was weak, but when considering the interaction of LCZ conversion and population density change, the ability of interpretation was significantly higher than that of a single factor. The results indicate that the distribution of LCZ, representing the urban spatial form, led to the main distribution pattern of sUHI and the conversion of LCZ influenced the sUHI intensity significantly. In addition, population density only had a significant relationship with sUHI in low and medium population density areas and its change alone cannot cause the sUHI pattern to change dramatically, but it intensified the influence of LCZ conversion on sUHI.
The spatiotemporal characteristics of sUHI and the effect of driving factors of LCZ and population density on sUHI in Beijing recognized in this research could help predicting future sUHI changes with ongoing urbanization and planning thermal comfort environment by city managers. However, there are still some aspects of this study that need to be improved. Firstly, a single thermal image for every summer may not reflect the thermal condition of summer accurately, although the moving averaged method was used. Secondly, although the relationship between LCZ, population and sUHI that has been discussed, the mechanism of how LCZ and population affect sUHI is not clear. Thirdly, the absence of local meteorological data allows the influence of lower atmospheric conditions on sUHI to be unaccounted for. Although the study selected data from clear days that could exclude the influence of precipitation, it is not possible to confirm the influence of, for example, dust content and humidity, which affect the radiation reaching the ground and the counter-radiation of the atmosphere. Other atmospheric conditions such as wind can also affect sUHI by influencing ventilation with different cooling effects on the land surface. These effects will introduce uncertainties in the accuracy and comparability of the sUHI calculated using LST in this paper. Finally, the time series used for trend analysis is also relatively short. If the time interval is expanded to 5 years where trends can be seen more clearly, only three intervals, 2000–2005, 2005–2010, and 2010–2015, can be used for analysis, and trend analysis is vulnerable to data fluctuation. Therefore, more relevant data should be added to acquire a more accurate representation of summer land surface thermal conditions, and a mechanism model-based method should be applied to explore the path of LCZ and population to influencing sUHI in future studies. In order to identify the influence of the lower atmospheric conditions on the sUHI effect, local meteorological observations in the study area should also be introduced in future studies. The distribution and intensity of the sUHI effect at different times calculated by the LST should be checked to ensure that they are obtained under similar lower atmospheric conditions in order to exclude interference, or at least to identify possible major deviations. In addition, longer time series data are needed in the future study to obtain a clearer trend of the sUHI.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/atmos12101271/s1, Table S1. Local climate zones scheme (Stewart and Oke, 2012); Table S2. Regression coefficients for TM\ETM+ and TIRS for different temperature ranges; Table S3. Emissivity of typical land cover in thermal infrared bands of different sensors; Table S4. Estimation of the atmospheric transmittance for TM/ETM+; Table S5. Estimation of the atmospheric transmittance for TIRS for the water vapor content range of 0.5~3 g∙cm−2; Table S6. Equations for estimating the average atmospheric temperature; Table S7. Year matching between distribution maps of LST and LCZ. Table S8. Lowest, highest and mean LST (°C) for each LST distribution map; Table S9. Confusion matrix of LCZ classification of 2003; Table S10. Confusion matrix of LCZ classification of 2005; Table S11. Confusion matrix of LCZ classification of 2010; Table S12. Confusion matrix of LCZ classification of 2017; Figure S1. Population density distribution maps of 2003, 2005, 2010, and 2017.

Author Contributions

Conceptualization, Y.Z. and S.L.; Data curation, Y.Z.; Formal analysis, L.L., J.S. and F.W.; Methodology, Y.Z. and S.L.; Software, D.L. and Z.L.; Validation, D.L., L.L., Z.L., J.S. and F.W.; Visualization, Y.Z.; Writing—original draft, Y.Z.; Writing—review & editing, Y.Z. and S.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Natural Science Foundation of China, grant number 41590843. The APC was funded by National Natural Science Foundation of China, grant number 41590843.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Grimm, N.B.; Grove, J.G.; Pickett, S.T.A.; Redman, C.L. Integrated Approaches to Long-Term Studies of Urban Ecological Systems: Urban ecological systems present multiple challenges to ecologists—pervasive human impact and extreme heterogeneity of cities, and the need to integrate social and ecological approaches, concepts, and theory. Bioscience 2000, 50, 571–584. [Google Scholar] [CrossRef] [Green Version]
  2. Angel, S.; Parent, J.; Civco, D.L.; Blei, A.; Potere, D. The dimensions of global urban expansion: Estimates and projections for all countries, 2000–2050. Prog. Plan. 2011, 75, 53–107. [Google Scholar] [CrossRef]
  3. Seto, K.C.; Guneralp, B.; Hutyra, L. Global Forecasts of Urban Expansion to 2030 and Direct Impacts on Biodiversity and Carbon Pools. Proc. Natl. Acad. Sci. USA 2012, 109, 16083–16088. [Google Scholar] [CrossRef] [Green Version]
  4. Cao, G.-Y.; Chen, G.; Pang, L.-H.; Zheng, X.-Y.; Nilsson, S. Urban growth in China: Past, prospect, and its impacts. Popul. Environ. 2012, 33, 137–160. [Google Scholar] [CrossRef]
  5. Kalnay, E.; Cai, M. Impact of urbanization and land-use change on climate. Nature 2003, 423, 528–531. [Google Scholar] [CrossRef]
  6. Zhou, L.M.; Dickinson, R.E.; Tian, Y.H.; Fang, J.Y.; Li, Q.X.; Kaufmann, R.K.; Tucker, C.J.; Myneni, R.B. Evidence for a significant urbanization effect on climate in China. Proc. Natl. Acad. Sci. USA 2004, 101, 9540–9544. [Google Scholar] [CrossRef] [Green Version]
  7. Koellner, T.; Scholz, R.W. Assessment of land use impacts on the natural environment—Part 1: An analytical framework for pure land occupation and land use change. Int. J. Life Cycle Assess. 2007, 12, 16–23. [Google Scholar] [CrossRef] [Green Version]
  8. Chen, Y.C.; Lin, T.P.; Lin, C.T. A simple approach for the development of urban climatic maps based on the urban characteristics in Tainan, Taiwan. Int. J. Biometeorol. 2017, 61, 1029–1041. [Google Scholar] [CrossRef] [PubMed]
  9. Grimm, N.B.; Faeth, S.H.; Golubiewski, N.E.; Redman, C.L.; Wu, J.; Bai, X.; Briggs, J.M. Global change and the ecology of cities. Science 2008, 319, 756–760. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Oke, T.R. The distinction between canopy and boundary-layer urban heat islands. Atmosphere 1976, 14, 268–277. [Google Scholar] [CrossRef] [Green Version]
  11. Oke, T.R. The Energetic Basis of the Urban Heat-Island. Q. J. R. Meteorol. Soc. 1982, 108, 1–24. [Google Scholar] [CrossRef]
  12. Ma, T.; Zhou, C.H.; Pei, T.; Haynie, S.; Fan, J.F. Quantitative estimation of urbanization dynamics using time series of DMSP/OLS nighttime light data: A comparative case study from China’s cities. Remote Sens. Environ. 2012, 124, 99–107. [Google Scholar] [CrossRef]
  13. Li, D.; Bou-Zeid, E. Synergistic Interactions between Urban Heat Islands and Heat Waves: The Impact in Cities Is Larger than the Sum of Its Parts. J. Appl. Meteorol. Clim. 2013, 52, 2051–2064. [Google Scholar] [CrossRef] [Green Version]
  14. Sun, Y.; Zhang, X.B.; Ren, G.Y.; Zwiers, F.W.; Hu, T. Contribution of urbanization to warming in China. Nat. Clim. Chang. 2016, 6, 706–709. [Google Scholar] [CrossRef]
  15. Chapman, S.; Watson, J.E.M.; Salazar, A.; Thatcher, M.; McAlpine, C.A. The impact of urbanization and climate change on urban temperatures: A systematic review. Landsc. Ecol. 2017, 32, 1921–1935. [Google Scholar] [CrossRef]
  16. Allen, L.; Lindberg, F.; Grimmond, C.S.B. Global to city scale urban anthropogenic heat flux: Model and variability. Int. J. Climatol. 2011, 31, 1990–2005. [Google Scholar] [CrossRef]
  17. Lin, T.P.; Chen, Y.C.; Matzarakis, A. Urban thermal stress climatic mapping: Combination of long-term climate data and thermal stress risk evaluation. Sustain. Cities Soc. 2017, 34, 12–21. [Google Scholar] [CrossRef]
  18. Lin, T.P.; Lin, F.Y.; Wu, P.R.; Hammerle, M.; Hofle, B.; Bechtold, S.; Hwang, R.L.; Chen, Y.C. Multiscale analysis and reduction measures of urban carbon dioxide budget based on building energy consumption. Energy Build. 2017, 153, 356–367. [Google Scholar] [CrossRef]
  19. Budhiraja, B.; Pathak, P.; Agrawal, G. Spatio-Temporal Variability of Urban Heat Islands in Local Climate Zones of Delhi NCR. Remote Sens. Technol. Appl. Urban Environ. II 2017, 10431. [Google Scholar] [CrossRef]
  20. Skarbit, N.; Stewart, I.D.; Unger, J.; Gal, T. Employing an urban meteorological network to monitor air temperature conditions in the ‘local climate zones’ of Szeged, Hungary. Int. J. Climatol. 2017, 37, 582–596. [Google Scholar] [CrossRef]
  21. Khamchiangta, D.; Dhakal, S. Physical and non-physical factors driving urban heat island: Case of Bangkok Metropolitan Administration, Thailand. J Environ. Manag. 2019, 248, 109285. [Google Scholar] [CrossRef]
  22. Yang, X.S.; Peng, L.L.H.; Jiang, Z.D.; Chen, Y.; Yao, L.Y.; He, Y.F.; Xu, T.J. Impact of urban heat island on energy demand in buildings: Local climate zones in Nanjing. Appl. Energy 2020, 260, 114279. [Google Scholar] [CrossRef]
  23. Zhou, X.L.; Okaze, T.; Ren, C.; Cai, M.; Ishida, Y.; Watanabe, H.; Mochida, A. Evaluation of urban heat islands using local climate zones and the influence of sea-land breeze. Sustain. Cities Soc. 2020, 55, 102060. [Google Scholar] [CrossRef]
  24. Stewart, I.D.; Oke, T.R. Local Climate Zones for Urban Temperature Studies. Bull. Am. Meteorol. Soc. 2012, 93, 1879–1900. [Google Scholar] [CrossRef]
  25. Stewart, I.D.; Oke, T.R.; Krayenhoff, E.S. Evaluation of the ‘local climate zone’ scheme using temperature observations and model simulations. Int. J. Climatol. 2014, 34, 1062–1080. [Google Scholar] [CrossRef]
  26. Wang, J.; Xu, C. Geodetector: Principle and prospective. Acta Geogr. Sin. 2017, 72, 116–134. [Google Scholar]
  27. Cai, M.; Ren, C.; Xu, Y.; Lau, K.K.-L.; Wang, R. Investigating the relationship between local climate zone and land surface temperature using an improved WUDAPT methodology–A case study of Yangtze River Delta, China. Urban Clim. 2018, 24, 485–502. [Google Scholar] [CrossRef]
  28. Mao, Q.; Long, Y.; Wu, K. Spatio-temporal changes of population density and exploration on urbanization pattern in China: 2000–2010. City Plan Rev. 2015, 39, 38–43. [Google Scholar]
  29. Liu, K.; Su, H.B.; Li, X.K.; Wang, W.M.; Yang, L.J.; Liang, H. Quantifying Spatial-Temporal Pattern of Urban Heat Island in Beijing: An Improved Assessment Using Land Surface Temperature (LST) Time Series Observations From LANDSAT, MODIS, and Chinese New Satellite GaoFen-1. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2016, 9, 2028–2042. [Google Scholar] [CrossRef]
  30. Peng, J.; Xie, P.; Liu, Y.X.; Ma, J. Urban thermal environment dynamics and associated landscape pattern factors: A case study in the Beijing metropolitan region. Remote Sens. Environ. 2016, 173, 145–155. [Google Scholar] [CrossRef]
  31. Yue, W.Z.; Liu, X.; Zhou, Y.Y.; Liu, Y. Impacts of urban configuration on urban heat island: An empirical study in China mega-cities. Sci. Total Environ. 2019, 671, 1036–1046. [Google Scholar] [CrossRef]
  32. Wang, J.; Wang, K.; Wang, P. Urban heat(or cool) island over Beijing from MODIS land surface temperature. J. Remote Sens. 2007, 11, 330–339. [Google Scholar]
  33. Zhou, J.; Chen, Y.H.; Li, J.; Weng, Q.H. A Volume Model for Urban Heat Island Based on Remote Sensing Imagery and Its Application: A Case Study in Beijing. AGU Fall Meet. Abstr. 2007, 2007, B21A-0029. [Google Scholar]
  34. Quan, J.; Chen, Y.; Zhan, W.; Wang, J.; Voogt, J.; Wang, M. Multi-temporal trajectory of the urban heat island centroid in Beijing, China based on a Gaussian volume model. Remote Sens. Environ. 2014, 149, 33–46. [Google Scholar] [CrossRef]
  35. Ma, M.; Lu, Z.; Sun, Y. Population growth, urban sprawl and landscape integrity of Beijing City. Int. J. Sustain. Dev. World Ecol. 2008, 15, 326–330. [Google Scholar] [CrossRef]
  36. Bechtel, B.; Alexander, P.J.; Bohner, J.; Ching, J.; Conrad, O.; Feddema, J.; Mills, G.; See, L.; Stewart, I. Mapping Local Climate Zones for a Worldwide Database of the Form and Function of Cities. ISPRS Int. J. Geo-Inf. 2015, 4, 199–219. [Google Scholar] [CrossRef] [Green Version]
  37. Conrad, O.; Bechtel, B.; Bock, M.; Dietrich, H.; Fischer, E.; Gerlitz, L.; Wehberg, J.; Wichmann, V.; Böhner, J. System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. Geosci. Model Dev. 2015, 8, 1991. [Google Scholar] [CrossRef] [Green Version]
  38. Qin, Z.; Karnieli, A.; Berliner, P. A mono-window algorithm for retrieving land surface temperature from Landsat TM data and its application to the Israel-Egypt border region. Int. J. Remote Sens. 2001, 22, 3719–3746. [Google Scholar] [CrossRef]
  39. Zhang, Y.; Yu, T.; Gu, X.-f.; Zhang, Y.-x.; Chen, L.-f. Land surface temperature retrieval from CBERS-02 IRMSS thermal infrared data and its applications in quantitative analysis of urban heat island effect. J. Remote Sens. 2006, 10, 789. [Google Scholar]
  40. Wang, J.F.; Li, X.H.; Christakos, G.; Liao, Y.L.; Zhang, T.; Gu, X.; Zheng, X.Y. Geographical Detectors-Based Health Risk Assessment and its Application in the Neural Tube Defects Study of the Heshun Region, China. Int. J. Geogr. Inf. Sci. 2010, 24, 107–127. [Google Scholar] [CrossRef]
  41. Wang, J.F.; Zhang, T.L.; Fu, B.J. A measure of spatial stratified heterogeneity. Ecol. Indic. 2016, 67, 250–256. [Google Scholar] [CrossRef]
  42. Liu, L.; Zhang, Y.Z. Urban Heat Island Analysis Using the Landsat TM Data and ASTER Data: A Case Study in Hong Kong. Remote Sens. 2011, 3, 1535–1552. [Google Scholar] [CrossRef] [Green Version]
  43. Singh, P.; Kikon, N.; Verma, P. Impact of land use change and urbanization on urban heat island in Lucknow city, Central India. A remote sensing based estimate. Sustain. Cities Soc. 2017, 32, 100–114. [Google Scholar] [CrossRef]
  44. Shirani-Bidabadi, N.; Nasrabadi, T.; Faryadi, S.; Larijani, A.; Roodposhti, M.S. Evaluating the spatial distribution and the intensity of urban heat island using remote sensing, case study of Isfahan city in Iran. Sustain. Cities Soc. 2019, 45, 686–692. [Google Scholar] [CrossRef]
  45. Sultana, S.; Satyanarayana, A.N.V. Assessment of urbanisation and urban heat island intensities using landsat imageries during 2000-2018 over a sub-tropical Indian City. Sustain. Cities Soc. 2020, 52, 101846. [Google Scholar] [CrossRef]
  46. Liu, X.; Zhou, Y.Y.; Yue, W.Z.; Li, X.C.; Liu, Y.; Lu, D.B. Spatiotemporal patterns of summer urban heat island in Beijing, China using an improved land surface temperature. J. Clean Prod. 2020, 257, 120529. [Google Scholar] [CrossRef]
  47. Meng, Q.Y.; Zhang, L.L.; Sun, Z.H.; Meng, F.; Wang, L.; Sun, Y.X. Characterizing spatial and temporal trends of surface urban heat island effect in an urban main built-up area: A 12-year case study in Beijing, China. Remote Sens. Environ. 2018, 204, 826–837. [Google Scholar] [CrossRef]
  48. Wu, W.J.; Zhao, S.Q.; Zhu, C.; Jiang, J.L. A comparative study of urban expansion in Beijing, Tianjin and Shijiazhuang over the past three decades. Landsc. Urban Plan 2015, 134, 93–106. [Google Scholar] [CrossRef]
  49. Li, L.; Zha, Y.; Zhang, J.H. Spatial and dynamic perspectives on surface urban heat island and their relationships with vegetation activity in Beijing, China, based on Moderate Resolution Imaging Spectroradiometer data. Int. J. Remote Sens. 2020, 41, 882–896. [Google Scholar] [CrossRef]
  50. Alexander, P.J.; Mills, G. Local Climate Classification and Dublin’s Urban Heat Island. Atmosphere 2014, 5, 755–774. [Google Scholar] [CrossRef] [Green Version]
  51. Lelovics, E.; Unger, J.; Gal, T.; Gal, C.V. Design of an urban monitoring network based on Local Climate Zone mapping and temperature pattern modelling. Clim. Res. 2014, 60, 51–62. [Google Scholar] [CrossRef] [Green Version]
  52. Leconte, F.; Bouyer, J.; Claverie, R.; Petrissans, M. Using Local Climate Zone scheme for UHI assessment: Evaluation of the method using mobile measurements. Build. Environ. 2015, 83, 39–49. [Google Scholar] [CrossRef]
  53. Geletic, J.; Lehnert, M.; Dobrovolny, P. Land Surface Temperature Differences within Local Climate Zones, Based on Two Central European Cities. Remote Sens. 2016, 8, 788. [Google Scholar] [CrossRef] [Green Version]
  54. Wang, C.Y.; Middel, A.; Myint, S.W.; Kaplan, S.; Brazel, A.J.; Lukasczyk, J. Assessing local climate zones in arid cities: The case of Phoenix, Arizona and Las Vegas, Nevada. ISPRS J. Photogramm. Remote Sens. 2018, 141, 59–71. [Google Scholar] [CrossRef]
  55. Koc, C.B.; Osmond, P.; Peters, A.; Irger, M. Understanding Land Surface Temperature Differences of Local Climate Zones Based on Airborne Remote Sensing Data. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2018, 11, 2724–2730. [Google Scholar] [CrossRef]
  56. Quan, J.L. Diurnal Land Surface Temperature Characteristics of Local Climate Zones: A Case Study in Beijing, China. In Proceedings of the 2019 IEEE International Geoscience and Remote Sensing Symposium (Igarss 2019), Yokohama, Japan, 28 July–2 August 2019; pp. 7443–7446. [Google Scholar]
  57. Wang, R.; Cai, M.; Ren, C.; Bechtel, B.; Xu, Y.; Ng, E. Detecting multi-temporal land cover change and land surface temperature in Pearl River Delta by adopting local climate zone. Urban Clim. 2019, 28, 100455. [Google Scholar] [CrossRef]
  58. Kotharkar, R.; Surawar, M. Land Use, Land Cover, and Population Density Impact on the Formation of Canopy Urban Heat Islands through Traverse Survey in the Nagpur Urban Area, India. J. Urban Plan Dev. 2016, 142, 04015003. [Google Scholar] [CrossRef]
  59. Manoli, G.; Fatichi, S.; Schlapfer, M.; Yu, K.L.; Crowther, T.W.; Meili, N.; Burlando, P.; Katul, G.G.; Bou-Zeid, E. Magnitude of urban heat islands largely explained by climate and population. Nature 2019, 573, 55–60. [Google Scholar] [CrossRef]
  60. Zhang, H.; Qi, Z.F.; Ye, X.Y.; Cai, Y.B.; Ma, W.C.; Chen, M.N. Analysis of land use/land cover change, population shift, and their effects on spatiotemporal patterns of urban heat islands in metropolitan Shanghai, China. Appl. Geogr. 2013, 44, 121–133. [Google Scholar] [CrossRef]
  61. Du, H.Y.; Wang, D.D.; Wang, Y.Y.; Zhao, X.L.; Qin, F.; Jiang, H.; Cai, Y.L. Influences of land cover types, meteorological conditions, anthropogenic heat and urban area on surface urban heat island in the Yangtze River Delta Urban Agglomeration. Sci. Total Environ. 2016, 571, 461–470. [Google Scholar] [CrossRef]
  62. He, C.; Zhou, L.G.; Yao, Y.R.; Ma, W.C.; Kinney, P.L. Estimating spatial effects of anthropogenic heat emissions upon the urban thermal environment in an urban agglomeration area in East China. Sustain. Cities Soc. 2020, 57, 102046. [Google Scholar] [CrossRef]
  63. Yin, J.; Wu, X.X.; Shen, M.G.; Zhang, X.L.; Zhu, C.H.; Xiang, H.X.; Shi, C.M.; Guo, Z.Y.; Li, C.L. Impact of urban greenspace spatial pattern on land surface temperature: A case study in Beijing metropolitan area, China. Landsc. Ecol. 2019, 34, 2949–2961. [Google Scholar] [CrossRef]
  64. Hu, D.; Meng, Q.Y.; Zhang, L.L.; Zhang, Y. Spatial quantitative analysis of the potential driving factors of land surface temperature in different “Centers” of polycentric cities: A case study in Tianjin, China. Sci. Total Environ. 2020, 706, 135244. [Google Scholar] [CrossRef] [PubMed]
  65. Li, T.; Cao, J.F.; Xu, M.X.; Wu, Q.Y.; Yao, L. The influence of urban spatial pattern on land surface temperature for different functional zones. Landsc. Ecol. Eng. 2020, 16, 249–262. [Google Scholar] [CrossRef]
  66. Elsayed, I.S. Effects of population density and land management on the intensity of urban heat islands: A case study on the city of Kuala Lumpur, Malaysia. In Application of Geographic Information Systems; IntechOpen: London, UK, 2012; pp. 267–283. [Google Scholar]
Figure 1. Study area: (a) location map; (b) Google Earth high-resolution image of summer 2017.
Figure 1. Study area: (a) location map; (b) Google Earth high-resolution image of summer 2017.
Atmosphere 12 01271 g001
Figure 2. LST distribution maps from 2000 to 2017.
Figure 2. LST distribution maps from 2000 to 2017.
Atmosphere 12 01271 g002
Figure 3. Accuracy of LCZ classification results at different resolutions.
Figure 3. Accuracy of LCZ classification results at different resolutions.
Atmosphere 12 01271 g003
Figure 4. LCZ classification maps (2003, 2005, 2010, and 2017).
Figure 4. LCZ classification maps (2003, 2005, 2010, and 2017).
Atmosphere 12 01271 g004
Figure 5. Trend of pixel numbers of (a) heat island and (b) cool island from 2000 to 2017.
Figure 5. Trend of pixel numbers of (a) heat island and (b) cool island from 2000 to 2017.
Atmosphere 12 01271 g005
Figure 6. Variation trend of pixel numbers from 2000 to 2017 for (a) weak, (b) medium, and (c) strong sUHI.
Figure 6. Variation trend of pixel numbers from 2000 to 2017 for (a) weak, (b) medium, and (c) strong sUHI.
Atmosphere 12 01271 g006
Figure 7. Spatial distribution of the Theil–Sen regression slope of LSTR.
Figure 7. Spatial distribution of the Theil–Sen regression slope of LSTR.
Atmosphere 12 01271 g007
Figure 8. Classification of variation characteristics of sUHI intensity (↗: continuous increase; ↗↘: first increase and then decrease; ↘↗: first decrease and then increase; ↘: continuous decrease).
Figure 8. Classification of variation characteristics of sUHI intensity (↗: continuous increase; ↗↘: first increase and then decrease; ↘↗: first decrease and then increase; ↘: continuous decrease).
Atmosphere 12 01271 g008
Figure 9. sUHI Intensity variation based on the LCZ scheme from 2000 to 2017.
Figure 9. sUHI Intensity variation based on the LCZ scheme from 2000 to 2017.
Atmosphere 12 01271 g009
Figure 10. The statistical results of distribution of LSTR of each LCZ type in 2003, 2005, 2010 and 2017.
Figure 10. The statistical results of distribution of LSTR of each LCZ type in 2003, 2005, 2010 and 2017.
Atmosphere 12 01271 g010
Figure 11. Median slope of LSTR series for each LCZ conversion type (the vertical coordinate represents the type of transition out, and the horizontal coordinate represents the type of transition in).
Figure 11. Median slope of LSTR series for each LCZ conversion type (the vertical coordinate represents the type of transition out, and the horizontal coordinate represents the type of transition in).
Atmosphere 12 01271 g011
Figure 12. Scatter plot and median curve of LSTR-population density (2003, 2005, 2010, and 2017).
Figure 12. Scatter plot and median curve of LSTR-population density (2003, 2005, 2010, and 2017).
Atmosphere 12 01271 g012
Figure 13. Spatial distribution of Spearman rank correlation coefficient and p value between LSTR and population density (p < 0.05).
Figure 13. Spatial distribution of Spearman rank correlation coefficient and p value between LSTR and population density (p < 0.05).
Atmosphere 12 01271 g013
Table 1. Observation dates and sensor platform of Landsat data.
Table 1. Observation dates and sensor platform of Landsat data.
YearObservation Date (MM/DD)Satellite
200008/20Landsat 7
200107/06Landsat 7
200208/26Landsat 7
200307/28Landsat 7
200407/06Landsat 5
200507/09Landsat 5
200608/21Landsat 7
200709/09Landsat 7
200808/02Landsat 5
200907/20Landsat 5
201008/08Landsat 5
201107/26Landsat 5
201208/21Landsat 7
201307/23Landsat 7
201409/04Landsat 8
201508/22Landsat 8
201608/08Landsat 8
201707/10Landsat 8
Table 2. Observation dates of the high-resolution images.
Table 2. Observation dates of the high-resolution images.
YearObservation Date (MM/DD)
200309/13
200508/18
201009/27
201708/23
Table 3. Attributes of Population Data.
Table 3. Attributes of Population Data.
Data attributesDescription
Period2000–2017
Temporal resolutionAnnual
RegionBeijing
Spatial resolution3 arc (approximately 81.7 m in the location of Beijing)
UnitPopulation numbers/pixel
Table 4. sUHI intensity grade classification.
Table 4. sUHI intensity grade classification.
LSTRsUHI Intensity Grade
<0Cool island
0–0.1Weak sUHI
0.1–0.2Moderate sUHI
0.2–0.4Strong sUHI
>0.4Extreme sUHI
Table 5. Interaction types between two variables.
Table 5. Interaction types between two variables.
CriterionInteraction Type
q(X1∩X2) < Min(q(X1), q(X2))Weak and nonlinear
Min(q(X1), q(X2)) < q(X1∩X2) < Max(q(X1), q(X2))Weak and univariate
q(X1∩X2) > Max(q(X1), q(X2))Enhanced and bivariate
q(X1∩X2) = q(X1) + q(X2)Independent
q(X1∩X2) > q(X1) + q(X2)Enhanced and nonlinear
Table 6. Factor and interaction detection results (dependent factor: sUHI intensity; independent factors: LCZ and population density).
Table 6. Factor and interaction detection results (dependent factor: sUHI intensity; independent factors: LCZ and population density).
Factor DetectorInteraction Detector
YearLCZ Population DensityLCZ ∩ Population Density
qpqpqInteraction Type
20030.6498990.0000.3022040.0000.700385Enhanced and bivariate
20050.6015970.0000.1168240.0000.632329Enhanced and bivariate
20100.6966600.0000.1724580.0000.733208Enhanced and bivariate
20170.7226660.0000.2147580.0000.761741Enhanced and bivariate
Table 7. Factor and interaction detection results (dependent factor: sUHI intensity change; independent factors: LCZ conversion types and population density change).
Table 7. Factor and interaction detection results (dependent factor: sUHI intensity change; independent factors: LCZ conversion types and population density change).
Risk DetectionInteraction Detection
PeriodLCZ Conversion TypesPopulation Density ChangeLCZ Conversion Types ∩ Population Density Change
qpqpqInteraction Type
2000–20090.2961890.0000.0208200.0000.378720Enhanced and nonlinear
2009–20170.3083530.0000.0056470.0000.421530Enhanced and nonlinear
2000–20170.3984480.0000.0177390.0000.485913Enhanced and nonlinear
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhang, Y.; Li, D.; Liu, L.; Liang, Z.; Shen, J.; Wei, F.; Li, S. Spatiotemporal Characteristics of the Surface Urban Heat Island and Its Driving Factors Based on Local Climate Zones and Population in Beijing, China. Atmosphere 2021, 12, 1271. https://doi.org/10.3390/atmos12101271

AMA Style

Zhang Y, Li D, Liu L, Liang Z, Shen J, Wei F, Li S. Spatiotemporal Characteristics of the Surface Urban Heat Island and Its Driving Factors Based on Local Climate Zones and Population in Beijing, China. Atmosphere. 2021; 12(10):1271. https://doi.org/10.3390/atmos12101271

Chicago/Turabian Style

Zhang, Yatong, Delong Li, Laibao Liu, Ze Liang, Jiashu Shen, Feili Wei, and Shuangcheng Li. 2021. "Spatiotemporal Characteristics of the Surface Urban Heat Island and Its Driving Factors Based on Local Climate Zones and Population in Beijing, China" Atmosphere 12, no. 10: 1271. https://doi.org/10.3390/atmos12101271

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop