Analyzing the interpretative ability of landscape pattern to explain thermal environmental effects in the Beijing-Tianjin-Hebei urban agglomeration

Classification of landscape

Previous research has shown that the relationship between a single landscape and surface temperature has been well studied (Zhao et al., 2018a). In real life, however, the surface temperature is the result of the mutual influence of various landscapes. Therefore, we explored the relationship between surface temperature and two kinds of integrated landscape pattern, as well as a single landscape.

Since the change of landscape pattern is related to human activities (Li et al., 2016), we divided the landscapes of the Beijing-Tianjin-Hebei urban agglomeration into six landscapes according to the degree of human influence on the landscape shaping process and primary classification results of land cover: (1) ForLS, mainly referring to the forest land in the land cover classification; (2) grassland landscape (GraLS), mainly referring to grassland in the land cover classification; (3) WetLS, mainly referring to rivers, lakes, and other water areas in the land cover classification; (4) farmland landscape (FarmLS), mainly referring to cultivated land in land cover classification; (5) ArtLS, mainly referring to the artificial surface in land cover classification; (6) bare land landscape, mainly referring to the unused land in the land cover classification. Ecological landscape refers to the types of land cover mainly providing ecosystem services, including forest, garden, grassland, water area, wetland, and other natural land, such as glaciers. Green landscape generally refers to the area with green plants, including park green space, production green space and protection green space, and so on where NDVI values tend to be high. Therefore, in this study, landscapes were integrated into two kinds of patterns, the ecological landscape pattern (EcoLSP, made up of three kinds of landscapes, ForLS, GraLS, and WetLS) and green space landscape pattern (GreenLSP), which was expressed using NDVI. In the study, we used all the single landscapes and integrated landscape patterns to explore the relationship between landscape pattern and surface temperature using the Pearson correlation analysis method and OLS regression method.

Division of study area according to terrain

In mountainous areas with complex environmental conditions, key topographic factors, such as altitude and slope, have a significant impact on surface temperature (Jiao et al., 2019; Duan et al., 2015). The terrain in the Beijing-Tianjin-Hebei urban agglomeration is complex, including hills that include the Taihang, Yanshan, and Jundu mountains and parts of the northern China plain, with significant differences in elevation between these different types of terrain.

To explore the influence of altitude on surface temperature, based on the grid of 1 × 1 km, we statistically analyzed the average surface temperature and average DEM values in different grids. We found a significant correlation between altitude and surface temperature using the Pearson correlation analysis method (the correlation coefficient was −0.68**). To reduce the influence of topographic factors on the relationship between landscape pattern and surface temperature, we divided the region in the Beijing-Tianjin-Hebei urban agglomeration into three types of terrain—plain, hills, and plateau (Fig. 3).

Spatial information expression of surface temperature

The Getis‐Ord G_i^* local hot-spot analysis method can be used to identify the spatial clustering regions with high value and low value of statistical significance. Therefore, in this part, we used this method to carry out clustering analysis and spatial heterogeneity expression of surface temperature. We combined the results of the hot-spot analysis with the landscape pattern and terrain to reveal landscape composition information for hot and cold spots in different types of terrain in this region. The formula of Getis‐Ord G_i^* (Feng et al., 2018) was as follows: (1) $$G_i^{\ast }={\displaystyle \frac{{\sum }_{j=1}^n{w}_{ij}{x}_j-\overline{x}{\sum }_{j=1}^n{w}_{ij}}{\sqrt{{\displaystyle \frac{{\sum }_{j=1}^n{x}_j^2}{n}-x^2}}\sqrt{{\displaystyle \frac{\left[n{\sum }_{j=1}^n{w}_{ij}^2-{\left({\sum }_{j=1}^n{w}_{ij}\right)}^2\right]}{n-1}}}},}$$where n represents the total number of factors in the study area; w_ij is the spatial weight coefficient between the factor i and factor j, reflecting the spatial relationship; x_j is the attribute value of the factor j; and x̄ is the average attribute value of n factors.

The method to determine the spatial weight coefficient is the key to determining the degree of influence between two factors. In this part, we selected the inverse distance weighting method to calculate the spatial weight coefficient. The traditional inverse distance weighting method was constructed primarily in the two-dimensional coordinate system. In hills, however, altitude has a significant impact on the Euclidean distance between two points. Therefore, we combined DEM data and the two-dimensional coordinate system to construct the three-dimensional coordinate system to calculate the distance and weight in the three-dimensional coordinate system. The formula of w_ij follows: (2) $$w_{ij}={\displaystyle \frac{{\left({\left(X_j-X_i\right)}^2+{\left(Y_j-Y_i\right)}^2+{\left(Z_j-Z_i\right)}^2\right)}^{-{\displaystyle \frac{p}{2}}}}{{\sum }_{a=1}^m{\left({\left(X_j-X_i\right)}^2+{\left(Y_j-Y_i\right)}^2+{\left(Z_j-Z_i\right)}^2\right)}^{-{\displaystyle \frac{p}{2}}}}\left(i\ne j\right),}$$where m represents the total number of factors in the local neighborhood of factor i, and it is 8 according to the best practice criterion; X_j, Y_j, and Z_j are the X, Y, and Z coordinate values of the factor j; X_i, Y_i, and Z_i are the X, Y, and Z coordinate values of the factor i; and p is the parameter, which is determined by the minimum value of the root mean square error, where it is 2. To avoid the situation in which the distance between two elements is 0, w_ij = 1 when i = j.

Z score, p values, and confidence level can reflect a statistically significant aggregation or discrete model of factors or the value. The z-score (Z(G_i^*)) is the statistic of G_i^*, which is a multiple of the standard deviation; p is the probability. The relationship among the three is shown in Table 2, and the relevant contents in the table refer to the help document of ArcGIS 10.2 software. The formula of Z(G_i^*) is as follows:

where E(G_i^*) is the expected value of the local exponent of Getis‐Ord G_i^*; and VAR(G_i^*) is the variance of the local exponent of Getis‐Ord G_i^*.

Exploration of the optimal analytical window granularity of the relationship between landscape pattern and surface temperature

The relationship between surface temperature and landscape pattern is correlated significantly with the analytical window granularity. If the analytical window granularity is too small, the calculation will be too large, the information of landscape pattern inside the window will be single, and the excessive analytical window granularity will affect the full embodiment of heterogeneity (Estoque, Murayama & Myint, 2017). Therefore, based on the resolution of surface temperature data (1 × 1 km), we selected eight kinds of analytical window granularity, including 1 × 1, 3 × 3, 5 × 5, 7 × 7, 9 × 9, 11 × 11, 13 × 13, and 15 × 15 km, and explored the optimal analytical window granularity of the relationship between landscape pattern and surface temperature in the Beijing-Tianjin-Hebei urban agglomeration.

Following the principle of grid statistical analysis, we divided the research area into multiple grids. On the basis of this principle, we then analyzed the average temperature, the PLAND of different landscapes and the average NDVI for each grid. Therefore, the phenomenon that the resolution of basic data was not identical has little effect on our research results.

We used the Pearson correlation analysis method to express the relationship between two factors by observing the correlation coefficient r (Lu et al., 2018; Miao, Ni & Borthwick, 2010). We determined the analytical window granularity by comparing r².

Interpretative ability of landscape pattern to explain surface temperature

On the basis of the optimal analytical window granularity, we combined OLS and GWR model to explore the interpretative ability and degree of influence of landscape pattern on surface temperature from different spatial levels.

The OLS model reflected the positive and negative relationship between the interpretative variables and explained the variable and the strength of the relationship by the regression coefficient (β). The R² of the results reflected the fitting effect of the evaluation model. The formula of the regression model and R² (Farahani et al., 2010) is as follows:

In formula (4), y_i represents the explained variable value of the sampling point i; x_i represents the interpretative variable value of the sampling point i; β₀ is the intercept, and β is the slope or regression coefficient of the interpretative variable; and ε_i is the residual. In formula (5), the means of x_i and y_i are the same as the means in formula (4), x̄ represents the average value of the interpretative variable, ȳ represents the average value of the explained variable, and n represents the total number of explained variables.

Although the OLS model can express the interpretative ability of landscape pattern to explain surface temperature as a whole, it cannot accurately reflect the spatial heterogeneity. On the basis of OLS model and GWR model, we incorporated spatial characteristics of the data and used subsamples from nearby observations to estimate each point in the study area. The regression results reflected the spatial heterogeneity of the interpretative ability and the degree of influence. The GWR formula (Guo, Ma & Zhang, 2008; Zhao et al., 2018a) is as follows: (6) $$y_i={\text{β}}_0(u_i,v_i)+\underset{p}{\overset{k=1}{{{\displaystyle \sum }}^{\text{}}}}{\text{β}}_k(u_i,v_i)x_{ik}+{\epsilon }_i,$$where y_i represents the value of the explained variable of the sampling point i, x_ik represents the value of the interpretative variable k at the sampling point i, (u_i, v_i) represents the coordinates of the sampling point i, β₀(u_i, v_i) represents the statistical regression constant of the sampling point i, β_k(u_i, v_i) represents the regression coefficient k of the sampling point i, its purpose is to estimate the adjacent space observations, p is the number of variables participating in the regression at a sampling point, and ε_i is the residual. We adopted the kernel function for the estimation function of β_k(u_i, v_i), and the following formula of β_k(u_i, v_i): (7) $${\text{β}}_k(u_i,v_i)={(X^{T}W(u_i,v_i)X)}^{-1}{X}^{T}W(u_i,v_i)y,$$

where X and X^T is the matrix and transpose of the factor values in the neighborhood of the observation point (u_i, v_i); and W(u_i, v_i) represents the distance weight matrix of the observation point (u_i, v_i), which is the observation function between the observation point and other points and is composed of factor W_ij in the distance weight matrix. The formula of W_ij is as follows: (8) $$W_{ij}=\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{d_{ij}^2}{h^2}}\right),$$

where h is the kernel bandwidth; and d_ij is the distance between the observation point i and the factors j in the neighborhood. In this study, we used Akaike information criterion (AIC) to determine the optimal bandwidth, and the weight outside the bandwidth was 0. AIC was a standard used to measure the goodness of statistical model fit, which was based on the concept of entropy and provided a standard to weigh the complexity of estimation model and the goodness of fitting data. In general, when AIC was minimum, bandwidth h was the best. The expression formula is as follows: (9) $$\mathrm{A}\mathrm{I}\mathrm{C}=2k-2\mathrm{l}\mathrm{n}\left(L\right),$$where k is the number of model parameters and represents the complexity of the model; and L is the likelihood function, which can reflect the degree of difference between different models.

Interpretative ability of landscape pattern to explain surface temperature in the study area

According to Table 4, we found that the interpretative ability of ForLS, FarmLS, ArtLS, GreenLSP, and EcoLSP to explain surface temperature was strong, whereas the interpretative ability of GraLS and WetLS to explain surface temperature was weaker. The area ratio of ForLS and EcoLSP increased by 10%, the surface temperature dropped by 0.69 and 0.63 °C, respectively; the NDVI value increased by 0.1 which could represent GreenLSP, the surface temperature dropped by 1.45 °C. The area ratio of ArtLS and FarmLS increased by 10%, and the surface temperature increased by 1.26 and 0.65 °C, respectively, which indicated that ArtLS and FarmLS was the main thermal landscape, moreover, the influence of ArtLS on surface temperature is twice as great as that of FarmLS.

Interpretative ability of landscape pattern to explain surface temperature in different terrains

The interpretative ability of ForLS, FarmLS, ArtLS, GraLS, and EcoLSP was stronger in hills (Table 5); the interpretative ability of the GreenLSP was stronger in hills and the plateau; and the interpretative ability of WetLS was stronger in the plain.

The area ratio of ForLS and EcoLSP increased by 10% and the surface temperature dropped by 0.64 and 0.63 °C in the hills and dropped by 0.45 and 0.41 °C in the plains and dropped by 0.99 and 0.18 °C in the plateau, respectively. The NDVI value increased by 0.1, the surface temperature dropped by 2.21 °C in the hills, whereas the surface temperature dropped by 0.10 °C in the plains. It was obvious that the cooling effect of the ForLS, EcoLSP, and GreenLSP in the hills was more significant.

Artificial landscape and FarmLS were the main thermal landscapes. The area ratio of the ArtLS increased by 10%, the surface temperature rose by 1.69 °C in the hills and increased by 0.32 and 0.17 °C, respectively in the plateau and plain region. The area ratio of FarmLS increased by 10%, the surface temperature rose by 0.67 °C in the hills and rose by 0.59 and 0.12 °C, respectively in the plateau and plains. The degree of influence of ArtLS and FarmLS was more significant in the hills.

Although the interpretative ability of GraLS and WetLS to explain surface temperature was not strong in the Beijing-Tianjin-Hebei urban agglomeration, the relationship among GraLS, WetLS and surface temperature was significantly positive in the hills and the plateau, whereas it was significantly negative in the plains.

Interpretative ability of landscape pattern to explain surface temperature using the GWR model

In the large-scale area of the Beijing-Tianjin-Hebei urban agglomeration, the area of GraLS and WetLS was small, and the relationship between GraLS, WetLS, and surface temperature was not obvious. Therefore, we studied only the interpretative ability of the ArtLS, FarmLS, GreenLSP, and EcoLSP to explain surface temperature and we divided the research results into five grades according to the quantile method (Fig. 7).

On the whole, EcoLSP and GreenLSP had a significant cooling effect. Increasing the area of ArtLS and FarmLS could promote the increase in surface temperature (Tables 6 and 7), which was consistent with the results given in Tables 4 and 5. Spatial heterogeneity was evident, however, in the interpretative ability and degree of influence of the four landscapes on surface temperature.

In hills, the class IV and V of EcoLSP R² values accounted for 76.13% and 62.03% of the terrain, respectively, and the class IV and V coefficients accounted for 61.32% and 78.85%, respectively (Table 7). This indicated that the interpretative ability and cooling ability of EcoLSP to explain surface temperature were greater in hills than in the plains and the plateau, and the interpretative ability was stronger at the border of the plains and hills, such as in the Yanshan and Taihang mountains (Fig. 7).

In hills, the class IV and V of GreenLSP R² values accounted for 75.25% and 83.82% of the terrain, respectively, and the class I and II coefficients accounted for 88.29% and 87.17% of the terrain, respectively (Table 7). The class V of GreenLSP R² was concentrated in the central hills, including the Yanshan Mountains. This coefficient showed a trend of gradual decrease from the southeast to the northwest in space (Fig. 7), which indicated that GreenLSP had a stronger interpretative power and a greater degree of influence on the surface temperature in the hills.

In hills, the class IV and V of FarmLS R² values accounted for 64.74% and 74.81% of the terrain, respectively, and the class IV and V coefficients accounted for 86.40% and 92.69% of the terrain, respectively (Table 7). These results showed that the interpretative ability and degree of influence of FarmLS on surface temperature in the hills was stronger. In hills, the class IV and V of ArtLS R² values accounted for 64.26% and 43.52% of the terrain, respectively. The percentage of the class V of ArtLS R² values in the plains was 56.20% and was focused on the surrounding areas of Beijing, Tianjin, and other cities, as well as on the junction of the hills and plains. In hills, the class IV and V coefficients accounted for 85.10% and 92.56% of the terrain, respectively (Table 7). This result showed that the interpretative ability of ArtLS to explain surface temperature was stronger in the hills and the plains, and the degree of influence of ArtLS on surface temperature was greater in the hills.