Multi-Scale Relationship between Land Surface Temperature and Landscape Pattern Based on Wavelet Coherence: The Case of Metropolitan Beijing, China

: The relationship between urban landscape pattern and land surface temperature (LST) is one of the core issues in urban thermal environment research. Although previous studies have shown a significant correlation between LST and landscape pattern, most were conducted at a single scale and rarely involve multi-scale effects of the landscape pattern. Wavelet coherence can relate the correlation between LST and landscape pattern to spatial scale and location, which is an effective multi-scale correlation method. In this paper, we applied wavelet coherence and Pearson correlation coefficient to analyze the multi-scale correlations between landscape pattern and LST, and analyzed the spatial pattern of the urban thermal environment during the urbanization of Beijing from 2004 to 2017 by distribution index of high-temperature center (HTC). The results indicated that the HTC of Beijing gradually expands from the main urban zone and urban function extended zone to the new urban development zone and far suburb zone, and develops from monocentric to polycentric spatial pattern. Land cover types, such as impervious surfaces and bare land, have a positive contribution to LST, while water and vegetation play a role in mitigating LST. The wavelet coherence and Pearson correlation coefficients showed that landscape composition and spatial configuration have significant effects on LST, but landscape composition has a greater effect on LST in Beijing metropolitan area. Landscape composition indexes (NDBI and NDVI) showed significant multi-scale characteristics with LST, especially at larger scales, which has a strong correlation on the whole transect. There was no significant correlation between the spatial configuration indexes (CONTAG, DIVISION, and LSI) and LST at smaller scales, only at larger scales near the urban area has a significant correlation. With the increase of the scale, Pearson correlation coefficient calculated by spatial rectangle sampling and wavelet coherence coefficient have the same trend, although it had some fluctuations in several locations. However, the wavelet coherence coefficient diagram was smoother and less affected by position and rectangle size, which more conducive to describe the correlation between landscape pattern index and LST at different scales and locations. In general, wavelet coherence provides a multi-scale method to analyze the relationship between landscape pattern and LST, helping to understand urban planning and land management to mitigate the factors affecting urban thermal environment.

purpose of this study is to quantify the multi-scale relationship between LST and landscape pattern in Beijing metropolitan areas. Specifically, the study aims were to (1) quantify the dynamic of LST spatial pattern in Beijing from 2004 to 2017 by distribution index; (2) identify the multi-scale correlation between landscape pattern (landscape composition and spatial configuration index) and LST using wavelet coherence; (3) calculate the Pearson correlation coefficient between landscape indexes and LST by spatial rectangle sampling method and compared with wavelet coherence; (4) summarize the performance of wavelet coherence for the multi-scale correlation analysis of landscape pattern and the factors affected by them.

Study Area
Beijing is situated at 115.7°-117.4° E and 39.4°-41.6° N, with an area of 16,410.54 km 2 . It includes 16 administrative areas, the northern and western regions are mainly mountainous areas, and the urban area is distributed in the southern plain area [38][39][40]. According to Beijing's urban master plan

Data and Preprocessing
Two remote sensing images obtained in late summer (September) with highly clear atmospheric conditions (cloud amount less than 5%) were used in this study (Table 1), and both the two remote sensing images were downloaded from the official website of USGS (https://earthexplorer.usgs.gov/). The two Landsat data sets in 2004 and 2017 are basically in accordance with the start year of Beijing Master Plan (2004-2020) and the start year of Beijing Master Plan (2016-2035) released by the government with detailed land-use types, which are very helpful for the validation of the land use classification from the two Landsat data sets. The climate of Beijing has the characteristics of the semihumid continental monsoon climate in the north temperature zone. The weather is cold and dry in winter, with short duration in spring and autumn. Peng et al. studied the spatial pattern of Beijing's urban thermal environment, and testified that urban heat island effects of summertime are more significant than that of other seasons, which can be considered as the reasonable season to study the urban thermal environment [9].

Sensor
Scene ID Acquisition Data Season Landsat-5 TM LT51230322004252BJC00 8 September 2004 Summer Landste-8 OLI/TIRS LC81230322017255LGN00 12 September 2017 Summer The land surface temperature was retrieved based on the single window algorithm after preprocessing including geometric correction, radiation correction, and atmospheric correction. Land cover maps were obtained by supervised classification and visual interpretation with a resolution of 30 m (Figure 2), which were classified into six classed including impervious surface (built-up land, roads, parking lots and other artificial impervious areas), urban green space (urban grass and park), water (rivers, lakes, and reservoirs), forest land (forest and shrubs), cultivated land and bare land (cradled and sand land). In order to improve the classification quality and solve misclassification error or salt-and-pepper noises generated by spectral confusion, Sieve Classes of ENVI was applied to remove the isolated points (the parameters were set as minimum threshold 2 pixels, window size 8 pixels), and visual interpretation of the classified maps and manual correction was conducted according to the existing two periods land use data (2004 and 2017) of Beijing Master Plan (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). To verify the classification accuracy, 34,596 points was randomly selected based on the proportion of each land, and the existing two periods land use data (2004 and 2017) of Beijing Master Plan (2004-2020) and historical images of Google Earth were referenced to make judgments and statistics on the classification accuracy of the sample points. The classification confusion matrix calculated from the samples and the classification accuracy of each land cover type are shown in Table 2, where the overall accuracy (OA) of the land cover maps were 83.84% (2004) and 84.87% (2017), and the kappa coefficients were 0.737 (2004) and 0.770 (2017), respectively.

Land Surface Temperature Retrieval
The LST was retrieved from thermal infrared band of remote sensing images based on the single window algorithm [9,35]. The steps are as follows: (1) converting the digital number values of the thermal infrared band to surface-leaving radiance; (2) converting the surface-leaving radiance to atsensor brightness temperature; (3) converting the at-sensor brightness temperature to LST [41,42].
First, the digital number (DN) of the thermal infrared band was converted to sensor radiance values ( ) using Equation (1) [43].
where is the minimum quantized calibration digital number, is the maximum quantized calibration digital number, and are the ratios of spectral irradiance to and , respectively. Then the sensor radiance values can be converted to surface-leaving radiance through an atmospheric correction process using Equation (2).
where is the surface leaving radiance, and are the upwelling radiance and downwelling radiance, respectively. and are the atmospheric transmission and emissivity of the surface, respectively. , , and were calculated using the web-based atmospheric correction tool [42]. can be obtained by Equation (3): where and are the emissivity of vegetation and urban surface, respectively. According to previous studies, the and are typically 0.99 and 0.92 [9,42]. includes the geographical distribution of the natural surface and internal reflection, and is the vegetation fraction. and can be obtained by the following Equations: where and are the minimum and maximum values of the normalized difference vegetation index, respectively. and are the near-infrared and infrared bands of remote sensing images, respectively. F is the shape factor (usually the mean value is 0.55).
Second, assuming that the Earth's surface is a black body, and converting the surface-leaving radiance to at-sensor brightness temperature ( ) using Equation (7) [44]. Finally, the at-sensor brightness temperature was converted to LST using Equation (8) [45].
Then, the calculated LST value was converted (Kelvin) to degrees Celsius (°C) by subtracting the constant 273.

High Temperature Center and Landscape Metrics
The meteorological conditions and acquisition time were not exactly the same for the two images with different temperature ranges of LST. The Jenks natural breaks classification method was applied to perform the cluster analysis on the data to minimize the intra-class differences and maximize the inter-class differences [9]. The LST was divided into four categories: low temperature (<23 °C), sublow temperature (23-28 °C), sub-high temperature (28-32 °C), high-temperature (>32 °C), according to the natural discontinuity points obtained by the natural fracture method of ArcGIS. The hightemperature level was defined as the high-temperature center (HTC). The spatial pattern of LST was studied by calculating the proportion of HTC and distribution index (DI) of different administrative regions and land cover types [9,46]. The equation for calculating the distribution index is as follows: where and were the total area of the HTC and space unit, respectively. is the total area, and is the area of the HTC in this spatial unit. = 1 indicates the distribution of the HTC in the entire research area. If < 1, it means that the HTC distribution in the space unit is lower than the entire research area, has a negative contribution to the LST, which contributes to cooling. If > 1, the HTC distribution in the space unit is higher than the entire research area, has a positive contribution to the LST, which contributes to warming.
To study the multi-scale relationship between urban landscape metrics and LST, the following landscape composition and configuration metrics were selected ( Table 2): Normalized Difference Built-up Index (NDBI) and Normalized Difference Vegetation Index (NDVI) were selected to describe the spatial distribution of urban impervious surface and vegetation, respectively; Shannon's diversity index (SHDI) was selected to describe landscape Diversity; Contagion Index (CONTAG) was selected to describe the degree of agglomeration or extension of the landscape; Landscape Division Index (DIVISION) was selected to describe the dispersion of the landscape; Landscape Shape Index (LSI) was selected to describe the spatial form of the landscape [26,47,48]. These landscape metrics are calculated at the landscape level of 1 km window size. Table 3 lists the calculation equations and descriptions of the selected landscape metrics. The contagion index refers to the non-random or aggregation degree of patch types in the landscape. Additionally, a larger metric value means larger aggregation degree, conversely lower.

−
A measure of the fragmentation of land covers. DIVISION deals with the degree to which the landscape is broken up into separate patches.
Landscape shape index refers to the deviation degree between the real shape of patch and the circle or square with the same area, indicating the complexity of landscape. , Short Wave Infrared; , Near Infrared; , Red Infrared; , proportion of the district i occupied by land cover i; , is the number of adjacencies between pixels of patch types i and k based on the double-count method; m is the number of patch types present in the landscape; , area of district i; , is the total length of edge for class i and patch j.

Continuous Wavelet Transform
Wavelet transform is derived from Fourier transform and is an improvement of the Windowed Fourier transform. It has good localization properties in the time and frequency domain [49]. The multi-scale characteristics of a signal can be obtained by stretching and shifting the mother wavelet, where Haar wavelet, Mexican Hat wavelet, Meyer wavelet, Daubechies wavelet, and Morlet wavelet are commonly used [37]. This paper used the Morlet wavelet, which is defined as follows: where and are dimensionless spatial parameters and frequency, respectively. To follow the convergence condition = 6 [50]. Continuous wavelet transform is defined as the convolution between the original function and the mother wavelet function, as follows: where and denote the scale parameter and translation parameter, respectively. The amplitude characteristics of the original data at different scales and locations can be obtained by changing and . Calculating the wavelet transform in the Fourier frequency domain can improve the calculation speed. The Fourier transform of the discrete data sequence is defined as: Equation (11) can be expressed as: The wavelet power spectrum is defined as | ( )| . Since the data length of the original transect is limited, it is necessary to consider the error caused by the edge effect. This can be solved by introducing the influence of the cone, padding the edges of data with zeros before the wavelet transform, and then eliminating them afterwards [33,51].

Wavelet Coherency Analysis
Similar to the traditional Pearson correlation method, wavelet coherence is used to analyze the correlation between two spatial series X and Y at a given spatial scale and location. Multi-scale wavelet coherence analysis can more reliably reveal the correlation between the two factors than the single-scale Pearson correlation coefficient [33,35]. Wavelet coherent analysis is defined as follows [33]: where ( ) = ( ) ( ) represents the cross wavelet transform. S is a smoothing operator and defined as: where indicates smoothing along the scale and indicates smoothing along the spatial location. The formula is defined as follows: where ∏ denotes rectangular function and 0.6 is the empirical value of the scale decorrelation length based on the Morlet wavelet [52]. In this paper, the Monte Carlo method was used to test the significance level of coherencies, and the wavelet coherence image showed the correlation and intensity between landscape metrics and LST at different scales and locations. Wavelet transform, coherence analysis, and Monte Carlo significance test were performed using MATLAB code [33].

Pearson Correlation Coefficient
The traditional Pearson correlation coefficient can be used to evaluate the linear relationship of the entire transect [28]: where x and y are represented LST and landscape metrics. r indicates the degree of correlation, the value is between −1 and 1, a higher value of r indicates a greater correlation. In order to make the Pearson correlation coefficient represent the correlation at different scales and positions, the spatial rectangle sampling method was used to calculate the correlation coefficient, and the rectangle size was the scale obtained by the wavelet transform. The correlation of the center position of the rectangle was represented by the Pearson coefficient calculated over the entire rectangle range, and then moved the rectangle along the cell size to calculate the correlation of the next location, finally, the correlation coefficient of the whole transect at different scales and locations was obtained. It was compared with wavelet coherence and the advantages and disadvantages of two multi-scale analysis methods was analyzed.

Urban Land Surface Temperature Dynamics
In 2004, the temperature range of Beijing LST was 19.41-40.98 °C, with an average temperature of 26.12 °C. In 2017 the temperature range of Beijing LST was 13.62-49.65 °C, with an average temperature of 27.88 °C (Figure 3). Figure 4 illustrates the spatial pattern of Beijing's surface temperature in 2004 and 2017, an obvious high-temperature area appeared in the main urban zone and urban function extended zone, which showed a monocentric spatial pattern. By 2017, there was still an obvious high-temperature area in the main urban zone and urban functional extended zone, and the HTCs had spread to the far suburb zone and new urban development zone, which was related to the process of urban expansion, and showed a significant polycentric spatial pattern.    12%). Spatially, the distribution of HTCs tended to decrease from the main urban zone to the suburbs. In 2017, the proportion of HTCs the main urban zone (84.7%) and urban function extended zone (49.94%) were relatively reduced, while the new urban development zone and far suburb zone increased to 20.03% and 3.07%, respectively. The variations of HTC indicated a significant outward spread tendency, which was consistent with the direction of urban (impervious surface) expansion (Figure 2), and the DI value gradually decreased from the main urban zone to the far suburb zone. In 2017, the DI value of the new urban development zone and the far suburb zone increased, while the main urban zone and the urban functional extended zone relatively reduced, which revealed the trend of the HTC gradually expanding to the surrounding zone of urban. Among the four functional zones, only the far suburbs have a DI value of less than 1, indicating that the HTC distribution is lower than the mean level of Beijing and has a cooling effect on urban LST.   The DI values of water and cultivated land were higher than that of forest land, mainly because water has a higher specific heat capacity and radiation relative to forest land, and is easily affected by the surrounding impervious surface; soil was exposed to the air after crop harvested, which increased the LST of cultivated land. Generally, an obvious cooling effect should exist in the urban green space, and a DI value less than 1, however, the actual value was 1.95 in 2004 and 1.17 in 2017, mainly because Beijing's green space patches are relatively small and scattered, when the high-temperature air of the impervious surface flows to the green space, after the green space absorbed the surrounding heat, the temperature gradually rises, and the DI value was higher than 1.

The Wavelet Coherencies between Landscape Composition and Land Surface Temperature
The correlations between NDBI and LST exhibited significant multi-scale characteristics ( Figure  5). In 2004, there was no significant correlation between NDBI and LST at smaller scales (less than 0.24 km). With the increase of scales, several regions with a strong correlation began to appear at the scale of 0.24-4.82 km (mainly concentrated in areas with large NDBI values) and varied with spatial location. The correlation was relatively low around 30 and 45 km of the transect, mainly due to the existence of large area urban green space and water in this zone. When the scale exceeds 4.82 km, NDBI was positively correlated with LST on the whole transect, which has a warming effect on LST. In 2017, at smaller scales, the correlation of NDBI and LST tended to increase compared with that in 2004. The correlation mainly distributed at 0-30 km and 45-75 km as the scale increased to 0.48-7.68 km. At the scales over 10.86 km, the whole transects showed a strong positive correlation between NDBI and LST. The scale of low correlation near urban areas increased in 2017, mainly because the increase of impervious surface makes the distribution of urban green space and water more fragmented, which increased the scale where NDBI and LST showed a significant positive correlation in the whole transect.
There showed a strong multi-scale correlation between NDVI and LST ( Figure 5). In 2004, NDVI and LST showed a weak correlation at small scales (less than 0.24 km). At the scale of 0.24-0.96 km, a strong correlation begins to appear at several locations. At the scale of 0.96-7.68 km, there was a strong negative correlation in the range of 0-30 km and 50-75 km of the transect, while there was no significant correlation between 30 and 50 km (urban area), probably because the area of vegetation patches in the urban area was small and scattered, which made the NDVI value relatively low and the impact on LST was not obvious. When the scale exceeds 10.68 km, NDVI was negatively related to LST significantly on the whole transect, and has a cooling effect on LST. In 2017, the correlation between NDVI and LST on the scale of 0.96-10.68 km showed the same distribution characteristics as that in 2004. The difference is that when the scale is less than 0.96 km, the whole transects showed a strong negative correlation except several locations. It shows that at larger scales, the vegetation in the two-stage transect data has similar spatial distribution characteristics. The large difference below 0.96 km was probably because the government has built some recreational facilities, green belts, and parks to reduce the impact of urban heat islands on the living environment (for example, in order to host the 2008 Beijing Olympics, the government built the Olympic Park and some green belts to increase urban green space). These urban green space patches are small and scattered in the city, and have a large impact on the land surface temperature at small scales.
The multi-scale correlation between SHDI and LST was weaker than the other two landscape composition indices (Figure 5). At scales less than 3.84 km, wavelet coherent diagrams in 2004 and 2017 have significant correlations only at several locations. At scales of 3.84-15.36 km, there was a significant negative correlation in the urban area. In 2017, the correlation between SHDI and LST showed similar distribution characteristics to that of 2004, although the correlation area and location have expanded. Affected by the patch size of land cover types, the diversity increases first and then decreases with the change of scales. At the scale of 3.84-15.36 km, land cover types in urban areas presented higher diversity, and different land surface radiation within different land use types caused the various temperature distributions in these lands. The unbalanced distribution of LST in space accelerated the flow of land surface air and reduced the LST, which caused a significant negative correlation between SHDI and LST. In the wavelet coherence diagram, the X-axis represents the spatial position; the Y-axis represents the scale obtained from the correlation period; the solid black line represents the influence cone (COI); the thick black line represents the 95% confidence level; the arrow direction indicates the type of correlation, the arrow to the left indicates a negative correlation, while the arrow to the right indicates a positive correlation. There was no significant correlation between landscape configuration indexes and LST at smaller scales, while at larger scales there was a strong correlation in the urban area ( Figure 6). In 2004, the wavelet coherent diagram between CONTAG index and LST showed a strong positive correlation in the urban area at scales of 3.84-15.36 km, mainly because the impervious surfaces were the dominant patches, with high density and connectivity, the CONTAG index was larger in the urban area, and the large area of the impervious surface reduced the diffusion speed of LST. In 2017, the regions with strong positive correlations at these scales expand to both sides, which was the same trend as urban impervious surface expansion. For the DIVISION index, in 2004, there was a strong negative correlation in the urban area at the scales of 7.68-15.36 km. In 2017, the region with strong negative correlation expanded and the scale range increased to 3.84-15.36 km, probably because the urban land cover can be divided into several patches, which increased the diversity, and contributes to the diffusion of LST. For the LSI index, in 2004, a strong negative correlation was identified in the urban area at scales of 3.84-15.36 km. In 2017, the region with strong positive correlation expanded, mainly due to the expansion of the impervious surface, and the degree of fragmentation of the landscape and the complexity of the land type increased, conducive to spreading LST. In general, the influence of landscape configuration metrics on LST at small and medium scales were more complicated, and the opposite correlation characteristics may occur due to location. With the increase of scales, the correlation was more stable and significant.  Figures 7 and 8 show the Pearson correlation coefficient and wavelet coherence coefficient between landscape metrics and LST at multi-scales. The Pearson coefficient was calculated by the spatial rectangle sampling method, the rectangle size was set to the scale of the Y-axis of the wavelet coherence diagram (eight main scales were selected as an example), and the correlation coefficients at different scales was obtained by spatial the rectangle on the whole transect. At small scales, the rectangle used in calculating the Pearson coefficient was smaller, so that the correlation coefficient at different positions had large fluctuations, but the wavelet coherence coefficient was relatively smooth due to the smoothing operator of scale and position. With the increase of the scale, the Pearson correlation coefficient and the wavelet coherence coefficient have the same trend, although they had some fluctuations in several locations. The Pearson correlation coefficient still has jagged fluctuations at larger scales, while the wavelet coherence coefficient was smoother and closer to a straight line, which better reflected the true correlation between landscape metrics and LST. In general, at large scales, the correlation trends obtained by the two methods were similar, but the Pearson correlation coefficient calculated by the spatial rectangle sampling method was easily affected by the rectangle size and position, while the wavelet coherence was less affected by these factors. Wavelet coherence can better reveal the correlation between LST and landscape metrics at different scales and locations.

Discussion
In the past few decades, the accelerated urbanization of Beijing has rapidly increased the impervious surfaces such as roads, squares, and built-up land, which has dramatically changed the urban thermal environment [53]. To mitigate these problems, the satellite city development plans were proposed by government and urban planners to form a polycentric spatial pattern [9], as reflected in the land cover and LST maps of 2004-2017. In 2004, the high-temperature center of Beijing was concentrated in the main urban zone and urban function extended zone. In 2017, the hightemperature center was gradually spread to the peripheral urban zone, which made the thermal environment of the main urban zone improve apparently. In addition, according to the Beijing Urban Master Plan 2004-2020, the government has built some parks and green belts around the urban area to increase the urban green space and ameliorate the LST in the urban area [39]. NDVI and LSI showed a strong negative correlation at small scales in 2017, which also revealed the phenomenon of increased urban green space. However, with the in-filling development pattern between urban polycentric, most lands between the satellite city and the main urban zone were used for commercial finance and residential, and only smaller and scattered land were used for urban green space. The  impervious surface was gradually connected, which may have an adverse effect on the urban thermal environment. The significant impact of urban local thermal environment on urban ecology made it an obligation to understand the changes in land cover and spatial pattern of LST [54]. Previous studies have shown that LST were affected by a variety of factors, such as urban structure, land cover types, geographic location, and population [13,55,56]. Land cover type is one of the main factors affecting LST, compared to areas with abundant land cover types (especially the land types which could adjust the LST, such as water and vegetation), the single and continuous impervious surface contributes more to the urban LST. The water and vegetation in the Beijing urban zone are smaller and fragmented, showed a negative correlation with LST at small and medium scales, and have a cooling effect. While at large scales, the impermeable surface is the dominant patch with a significant positive correlation with LST, which reduced the cooling effect of water and vegetation. The spatial configuration indexes of the main urban area have a strong correlation with LST at large scales. Urbanization has increased the fragmentation of the Beijing landscape, and the urban land cover types have become more diverse and complex, which may be helpful to mitigate the LST in urban areas [29,38]. In addition, LST images obtained at different times are also affected by atmospheric radiation, solar radiation, wind speed, humidity, etc. The LST data obtained in 2004 and 2017 may not have been stable, which will have a certain effect on the analysis of the spatial pattern distribution characteristics of LST [57]. In this paper, the two LST data were selected at the beginning and end of Beijing's urban master plan (2004-2020), which represented the process of urban development of Beijing cities during a completed master plan period. The two periods of data sets provide the vision to investigate the process of land use change and the correlation of the corresponding temperature gradient during the pursuit of the 2004-2017 master plan. In addition, this paper focuses on the multiscale correlation between landscape index (landscape composition and spatial configuration) and the land surface temperature, and the two periods of data for subsequent wavelet coherent analysis are also helpful for the scale variation between the two periods.
The Pearson correlation diagram obtained by the spatial rectangle sampling at different locations and scales showed a similar trend to wavelet coherence, and the undulated tendency was mainly caused by the calculation of rectangle size and pixel size. In contrast, the wavelet coherent diagram was relatively smooth and less affected by the above factors, which was consistent with the real situation of the impact of landscape metrics on LST, and proved the effectiveness and adaptability of wavelet coherence in the multi-scale analysis. Although wavelet coherence provides a new perspective for multi-scale identification of the correlation between LST and landscape pattern, several deficiencies still existed. For example, the window size selected when calculating various types of landscape indices at the landscape level, and the transect selected when analyzing the correlation between the landscape index and the LST, these factors will inevitably affect the analysis results. Moreover, the ecological interpretation of wavelet coherent diagram is complex, and the scale needs to be more accurately and finely divided. Considering the diurnal and seasonal variation of temperature, the sensitivity of different indexes to the LST, as well as the impact of more time-series of LST images on the analysis results, more LST images, more effective landscape indexes, and precise scales may be applied in the future to study the LST dynamic and the correlation between urban landscape pattern and LST.

Conclusions
This study applied wavelet coherence and Pearson correlation coefficient to analyze the multiscale correlations between landscape pattern and LST, and analyzed the spatial pattern of the urban thermal environment during the urbanization of the Beijing from 2004 to 2017 using the distribution index of high-temperature centers. According to the experimental results, we drew the following conclusions:  The HTC of Beijing has a tendency to spread outwards from 2004 to 2017. In 2004, it mainly concentrated in the main urban zone and urban function extended zone, forming a spatial pattern of monocenter distribution. In 2017, the HTC gradually expanded to the new urban development zone and far suburb zone, with a spatial pattern of polycentric distribution, and gradually connected to the main urban zone and urban function extended zone, which has an adverse effect on urban LST.  Land cover types have a significant impact on the spatial pattern of LST. The mean LST and HTC distribution of the impervious surface were the highest, while the forest land and water body were lower. Urban green space has a cooling effect on the LST, but due to its small area and scattered distribution, easily affected by the large-area impervious surface, causing the mean temperature and HTC distribution to be higher than forest land and water body.  Wavelet coherent analysis showed that the correlation between landscape pattern metrics (landscape composition and configuration) and LST has significant multi-scale effects. The landscape composition indices NDBI and NDVI showed a correlation with LST at all scales, especially at large scales, which showed a strong positive correlation and negative correlation, respectively. The space configuration indices CONTAG, DIVISION, and LSI have no significant correlation with LST at smaller scales. With the increase of scale, it showed a strong correlation around the urban area.  The wavelet coherence and Pearson correlation coefficients showed that the landscape composition and spatial configuration have significant effects on the LST, but the landscape composition has a greater impact on the LST in the Beijing metropolitan area.  Totally, compared with Pearson correlation coefficient calculated by spatial rectangle sampling, the wavelet coherence diagram is smoother and less affected by the location and rectangle size, which is more conducive to describing the correlation between landscape pattern index and LST at different scales and locations.

Conflicts of Interest:
The authors declare that they do not have any conflicts of interest.