1. Introduction
In the process of urbanization, large-scale urban land development has become one of the primary ways in which human activities reshape the natural environment. One of the environmental impacts of urbanization is the phenomenon known as the urban heat island (UHI) effect, where the temperature in urban areas is higher than in the surrounding regions [
1]. This effect is exacerbated by the rapid urbanization and widespread transition from natural landscapes to impervious surfaces, leading to increased absorption of solar radiation at the surface and reduced evapotranspiration from natural vegetation [
2]. Furthermore, with the exacerbation of global warming, rapid urban expansion, and the increasing occurrence of extreme weather conditions, the UHI effect has significantly affected air quality, vegetation phenology, as well as the health and comfort of residents [
3]. Consequently, mitigating the UHI effect has become a focal point of research in various related fields.
The UHI effect can be assessed through air temperature or land surface temperature (LST). While air temperature is primarily measured by meteorological stations, the limited and scattered nature of these stations provides only partial insights into urban temperature variations. LST results from the energy flux and interactions between the Earth and the atmosphere, playing a regulatory role in the lower atmospheric temperature of cities [
4]. Scholars have conducted diverse studies on LST. For instance, Liu et al. utilized a three-step method, blending MODIS LST products and Landsat data to generate surface temperature for the summer months from 2003 to 2018. Their research revealed that the streets with high UHI effects in Beijing were mainly concentrated in the city center, while the streets with low UHI effects were predominantly in the suburbs [
5]. Yang et al., using a single-window algorithm in conjunction with Landsat 8 TM10 band, performed LST inversion and found a close correlation between the surface temperature in Shanghai and the distribution of buildings [
6]. David et al. investigated the relationship between simulated surface UHI (SUHI), thermal field variations, and land use indices. Their study identified the conversion of natural and agricultural land to urban or bare land as the primary cause of increased LST and SUHI [
7]. With the continuous development of satellite remote sensing technology, various satellite images (e.g., Landsat, MODIS) can rapidly and accurately reflect differences in urban surface temperature and urban morphology, providing a foundational dataset for local climate studies [
8]. Consequently, LST is widely employed to explore the relationship between the UHI effect and urban morphological indicators, such as land use/land cover types and landscape patterns [
9,
10].
Urban areas represent complex dynamic systems characterized by both two-dimensional (2D) and three-dimensional (3D) spatial features [
11]. Previous research on factors influencing LST has predominantly focused on 2D urban spatial morphology, examining the impact of green spaces, water bodies, and impervious surfaces on LST [
12]. Some studies have found that the impact of urban spatial morphology on LST is not always linear, but rather influenced by geographical location and seasonality [
1,
9,
13]. This implies that there may be variations in LST response among different regions and seasons. Additionally, research indicates a positive correlation between building area, main road area, and the Normalized Difference Built-up Index (NDBI) with LST. Conversely, water body area and Normalized Difference Vegetation Index (NDVI) have been found to be negatively correlated with LST [
14,
15]. This suggests that water bodies and vegetation play significant roles in reducing LST. Although these studies have selected various indicators at different scales to explore the relationship between urban spatial morphology and LST, they have predominantly focused on 2D spatial morphology. In contrast, 3D morphological information is a crucial characteristic of urbanization, with buildings serving as key components of urban structure and major influencers of the UHI effect. This is because buildings can alter the reflection and absorption of solar radiation and the heat diffusion within urban areas [
16].
With the acceleration of urbanization, the quantity and height of buildings in urban areas continue to rise, prompting researchers to focus on the relationship between the 3D spatial morphology of cities and LST. Factors such as building height, volume, density, and shape coefficient have been widely applied, and scholars have found significant correlations between these factors and LST [
17,
18]. As remote sensing technology continues to advance, an increasing number of scholars are contemplating the impact of both 2D and 3D urban morphology on LST. Research has found that both types of indicators influence LST. However, there is some controversy regarding whether 2D indicators or 3D indicators play a more significant role in shaping LST in urban areas. Alavipanah et al. argue that 3D indicators are more important in shaping the LST of different urban structures compared to 2D indicators [
19]. However, Huang et al. suggest that the influence of 2D morphology on LST is superior to that of 3D morphology [
20]. Therefore, empirical results regarding the impact of 2D and 3D urban environments on LST require further investigation.
In addition, many of these studies are based on entire cities, grids, or street units. For example, Koko et al. explored the relationship between land use/land cover changes and LST in Abuja, Nigeria. They found that LST is negatively correlated with NDVI and positively correlated with NDBI, indicating that urban expansion and reduced vegetation lead to increased LST [
21]. Chen et al. investigated the relationship between LST and urban morphology in the core area of Xi’an at different grid scales. They identified average building height as a seasonally stable factor with a cooling effect [
22]. Zhang et al. explored the relationship between LST and urban morphology in Beijing and Shanghai using streets as spatial units. They found that building density is the most important factor influencing LST, with 3D building forms being the most significant [
1]. However, these approaches cannot fully characterize the basic spatial units of a city, warranting further exploration, based on finer spatial units, to investigate factors influencing the urban thermal environment.
In addition to land cover and surface geometry, urban morphology encompasses a variety of functional zones associated with various human activities. Urban functional zones (UFZs) are delineated based on different physical characteristics, as well as social and economic functions, and they often serve as fundamental units in urban planning [
23]. Therefore, studying the 2D/3D urban morphology based on spatial units of UFZs holds greater practical significance for understanding the UHI effect. However, in the current literature, scholars have extensively studied the factors influencing LST at large scales, but there is a limited investigation into the influencing factors of LST at the finer scale of UFZs. Moreover, the intricate relationship between 2D and 3D urban morphology and LST requires a comparative analysis that considers 3D building features. This study aims to utilize multi-source geospatial data to more effectively represent the physical characteristics, as well as the social and economic functions within the study area, facilitating a fine-scale examination of the 2D/3D urban morphology’s response to the UHI effect.
Due to rapid urbanization, Beijing, as a mega-city in China, experiences a significant UHI effect. This study seeks to leverage multiple sources of geospatial data to retrieve the LST in Beijing and create a map of UFZs. Using Pearson correlation analysis and GeoDetector, the study aims to investigate the relationship between 2D/3D urban morphology and the UHI effect in different UFZs. The analysis results of Beijing have significant implications for other cities. The specific objectives of this study are as follows: (1) to retrieve and analyze the distribution of LST in Beijing; (2) to investigate the distribution characteristics of UFZs in Beijing; (3) to explore the relationship between different types of UFZs and LST; and (4) to examine the relationship between the 2D and 3D urban morphology parameters of UFZs and LST.
3. Method
The workflow of this study is illustrated in
Figure 2, encompassing four components:
(1) First, the extraction of LST within the Fifth Ring Road of Beijing was conducted using Landsat-8 remote sensing image.
(2) Second, using OSM road network data, we partitioned the urban area into basic units and conducted supervised and unsupervised identification of UFZs based on multi-source data and POI attributes.
(3) Third, the computation of 2D and 3D urban morphological factors for UFZs was carried out using land cover data and building vectors.
(4) Finally, analysis of the relationship between 2D/3D urban morphology and LST in UFZs was conducted using Pearson correlation analysis and GeoDetector.
3.1. LST Inversion
LST can be retrieved through various remote sensing sensors, such as MODIS and ASTER sensors, which are supported by the Terra and Aqua platforms [
28], as well as the Landsat satellite series, supported by sensors like ETM, TM, and OLI. Landsat-8, with its higher spatial resolution and suitable spectral bands for the study area, was employed in this study for LST inversion [
29]. Following the approaches of Feng et al. [
30] and Huang et al. [
20], an atmospheric correction method was applied to estimate the atmospheric impact on LST radiation. To ensure robust results, this study selected summer images for inversion, given that the UHI effect is more pronounced during this season. Therefore, using the Google Earth Engine platform, the average of LST inversion results from 1 August 2021 to 15 August 2021 was calculated. The computation of LST in this study is described as follows:
In the equations,
represents the top-of-atmosphere (TOA) atmospheric radiance;
denotes the surface emissivity;
signifies the thermal radiance derived from the Planck law at temperature
;
is the transmissivity of the atmosphere in the thermal infrared band; and
and
are the atmospheric downwelling and upwelling radiance, respectively. The relationship between
and
can be expressed as follows:
Thus, the LST (
) can be inverted from the Planck formula as follows:
For Landsat-8, the constants are K1 = 774.89 W/(m∙sr∙) and K2 = 1321.08 K.
3.2. UFZ Mapping
For UFZ mapping, this study proposes a two-stage method that combines both unsupervised and supervised classification to leverage the strengths of various data sources. Prior to UFZ classification, we initiated topological processing of OSM road network data. On this basis, correction of interruptions and duplicate routes was performed, establishing a fundamental framework for the classification of independent units of UFZs.
In the first stage, a Sentinel-2 image was initially utilized to compute the NDVI and NDWI values. Thresholds were then set to distinguish green areas and water bodies, respectively. Subsequently, the total area proportion of green areas and water bodies within each UFZ was calculated. Based on these proportions, green space zones were determined, as they typically have a larger proportion of green areas and water bodies. Next, the average nighttime light index for each UFZ was computed, and a threshold was set to identify commercial zones, as these zones typically exhibit higher levels of luminosity at night. Commercial and green space zones were extracted in this step. Following this, a supervised classification approach was employed to identify the types of the remaining UFZs.
In the second stage, the remaining UFZs were classified using multi-source features and a random forest classifier. We utilized multiple data sources as features for each UFZ, including the average values of NDVI, NDWI, nighttime light index, and five normalized kernel density estimates of POI measures. However, the quantity of POIs varies across different types, with commercial POIs outnumbering others. This results in an imbalanced distribution of point numbers among different POI types [
31]. To address these issues, we calculated the normalized kernel density of POI within each unit to determine UFZs. The function is expressed by the following formulas:
where
represents the kernel density function;
is the kernel function; n is the number of known points; h is the bandwidth; and
is the distance from the grid center point to the known point. In Equation (5),
is the normalized value of POI kernel density, and
and
are the minimum and maximum values of POI kernel density, respectively.
The categories of certain UFZs were determined through manual visual interpretation, along with the green spaces and commercial zones identified in the first stage, serving as the training dataset for the random forest classifier to identify the remaining UFZs. In this study, the dataset is randomly divided into a training set (comprising 70% of the data) and a testing set (30% of the data). Random forest was trained on the training set and evaluated on the testing set. The random forest algorithm integrates predictions from multiple decision trees, effectively reducing the risk of overfitting [
32].
To evaluate the accuracy of the proposed method, this study employed overall accuracy (OA) and Kappa coefficient metrics [
33]. OA represents the ratio of correctly classified samples to the total number of samples, while Kappa is used to assess the consistency between classification results and actual categories. The calculation formula for OA is as follows:
where
represents the number of UFZ categories,
denotes the number of correctly classified samples in each UFZ category, and
represents the total number of UFZ samples.
The calculation of OA is given by:
where
represents the chance-corrected agreement,
denotes the number of UFZ categories,
signifies the actual number of samples in each UFZ category, and
represents the number of samples misclassified into each category.
3.3. 2D/3D Urban Morphology Factors
Patches are fundamental units composing landscape patterns, representing relatively homogeneous non-linear areas distinct from the surrounding background [
34]. This study explored the relationship between the patchiness of land cover types and LST concerning UFZs, utilizing landscape metrics.
As shown in
Table 3, 2D urban morphology encompasses surface cover landscape indices, while 3D urban morphology involves building forms computed through architectural calculations. Three widely used landscape indices were employed in this study to assess urban landscape patterns: percentage of landscape area (PLAND), patch density (PD), and Shannon’s diversity index (SHDI). These indicators, calculated using FRAGSTATS 4.2, describe urban morphology from three perspectives: area proportion, shape complexity and diversity, and spatial arrangement. They effectively capture the ecological environment formed by the interaction of natural and anthropogenic factors within a region. Both PLAND and PD were calculated based on each land cover element, while SHDI was computed across all categories.
In addition, to focus on 3D architectural structures, six indices were constructed using building vectors, including the area-weighted mean shape index (shape index), building density (density), shape coefficient, mean height, height variance, and sky view factor (SVF). These six indicators encompass aspects of shape, composition, and distribution (
Table 3). Buildings represent the primary aspects of 3D urban morphology, and the selected indices cover differences in building height, shape, density, and other horizontal and vertical dimensions. They comprehensively measure the 3D morphology within UFZs. The calculation of the SVF is conducted using the UMEP plugin in QGIS software 3.16.
Figure 3 illustrates the building heights and SVF in a sample area. Considering the fact that the climate factors and the terrain factors are relatively consistent in the study area, this study only discusses the influence of 2D/3D urban morphological parameters on daytime urban LST.
3.4. Data Analysis Methods
Using machine learning models may introduce the risk of overfitting, especially with limited sample sizes. In contrast, Pearson correlation analysis and GeoDetector are simpler and more suitable for our research purposes. Additionally, they can provide direct relationships and the extent of influence among indicators, making the results easy to interpret and understand. Therefore, this study will first employ Pearson correlation analysis to identify variables that may have significant effects. Subsequently, variables will be input into the GeoDetector model to analyze the relationship between urban morphology parameters and LST. In addition, these two methods have been successfully applied in the related field [
20,
35].
3.4.1. Pearson Correlation Analysis
Considering the potential correlation between 2D/3D urban morphology indicators and LST, this study employed Pearson correlation analysis to investigate the relationship between 2D/3D urban morphology indicators and LST. Pearson correlation analysis is a statistical method used to measure the strength and direction of the linear relationship between two variables. It is based on the concept of covariance, calculated by dividing the covariance of two variables by the product of their respective standard deviations, resulting in a correlation coefficient ranging from −1 to 1 [
36]. The Pearson correlation formula is as follows:
where
represents 2D/3D urban morphology indicators,
represents LST, and n is the sample size. The correlation coefficient ranges from −1 to 1: when r > 0, it indicates a positive correlation between LST and the factor; when r < 0, it indicates a negative correlation; when r = 0, it indicates no correlation between LST and the factor.
3.4.2. GeoDetector
GeoDetector represents a set of statistical methods used to explore spatial variations and reveal the driving forces behind them. Its core concept is based on the assumption that if a certain independent variable significantly influences a dependent variable, their spatial distributions should exhibit similarity [
37]. In this study, the differentiation and factor detection module of GeoDetector is employed to assess the explanatory power of the influencing factors on the spatial variation of LST in UFZs, as indicated by the formula:
where
represents the explanatory power of a certain influencing factor on the spatial variation of LST in UFZs;
h = 1,2,3,…;
denotes the layers of LST or influencing factors;
and
represent the number of units in layer
h and the total number of UFZs, respectively; σ
h and σ
2 denote the variance of feature values in layer
h and UFZs, respectively. The
q value ranges from [0, 1], where a higher value indicates a stronger explanatory power of each influencing factor on the spatial variation of LST in UFZs [
38].
The interaction detection module is employed to detect the types of interactions between the independent variables (
Xi), determining whether the combination of different factors affects the explanatory power of the dependent variable (
Y). This method first calculates the
q value for each influencing factor and then computes the
q value for each pair of factors after interaction, as presented in
Table 4. A higher
q value indicates a stronger joint explanatory power of the two factors [
37].