RS and GIS Supported Urban LULC and UHI Change Simulation and Assessment

Key Laboratory of Spatio-temporal Information and Ecological Restoration of Mines (MNR), Henan Polytechnic University, Jiaozuo, Henan 454003, China School of Surveying and Mapping Land Information Engineering, Henan Polytechnic University, Jiaozuo, 454003 Henan, China College of Surveying and Geo-informatics, Tongji University, Shanghai 200092, China Collaborative Innovation Center of Aerospace Remote Sensing Information Processing and Application of Hebei Province, China


Introduction
Land use and land cover changes (LULCCs), as one of the most significant processes related to earth ecological environment problems and social progress issues [1][2][3][4][5], related to global and regional changes [6,7], have largely affected earth biochemical cycle [8,9], sustainable use of resources [10], biodiversity [11], and urban planning and policymaking [12]. As a specific type of LULCCs, urban sprawl which plays an important role in urban intelligent growth and modernization is a sign of development and progress of human civilization and urbanization [13]. By prediction, more than 60% of the human population will live in cities by 2030 [14] which results in increasing replacement of natural landscape by the human-made landscape and will cause the temperature in the urban area to be higher than that in the suburban surrounding or rural area [15], which is also known as UHI. UHI phenomenon will not only accelerate urban environmental temperature and air pollution but also significantly increase energy consumption, increase urban temperature, and reduce quality of life. UHI is realized as a significant factor leading to global warming, heat-related mortality, and nonforecast climate change. A comprehensive study of the influencing factors of the UHI effect is critical for formulating reasonable urban planning policies and mitigating the effects of UHI [16][17][18].
The urban sustainability development which can comprehensively consider the ecological environment and human environment at the natural, economic, and social levels in the process of urban growth, and reflect the integrity of human settlement environment and health status of the ecosystem as a whole, is regarded as the basis for regional urban environmental system governance and prevention policy formulation [19,20]. The urban expansion was found to have a significant impact on local temperatures, in Chapman et al.'s review article which found that in some cases by up to 5 degrees [16]. UHI is closely associated with urban structure and will further increase by urban sprawl [21]. In the past decades, most researchers examined LULCCs and UHI in isolation, with few considering their combined effect.
Furthermore, existing satellite-based studies have typically evaluated and assessed the status of UHI at any given time but have limitations for studying the dynamic progress of LULCCs and UHI. Most research focused on evaluating the current or past status of LULCCs or UHI, while there has been little attempt to simulate or predict future change even though this information is crucial to inform effective sustainable urban development policy. Change simulation can provide valuable information for future prediction, as well as can indicate anthropogenic impact and identify degradation and deforestation which is useful for urban development planning. Various dynamic prediction models, including empirical-statistical models, optimization models, agent-based models, and Cellular Automata-Markov (CA-Markov) models, have been used for LULCC simulation [22], while they are seldom used for UHI change trend prediction.
In this research, we proposed a strategy to assess and analyze the impact of LULCCs on UHIs, as well as to simulate, predict, and explore the relationship and interaction between LULCCs and thermal environment trend in the future. For this purpose, we employed multitemporal remotely sensed data captured in 1986, 1996, 2006, and 2016, GIS spatial analysis methods, and CA-Markov trend simulation model. We test the method on Zhengzhou city, one of the fastest-growing metropolitan cities. Outputs are expected to contribute to urban planning, urban security assessment, and sustainable development in urban environments.
The rest part of the paper is organized as follows: a brief introduction about the research area, preliminary work, and preprocessing of datasets is given in Section 2. The proposed research framework based on remote sensing (RS) is drawn and applied to the research area in Section 3. The results are shown in Section 4; analysis and discussion are given in Section 5. Finally, the conclusions are drawn in Section 6.

Study Area and Datasets
The study site is located in Zhengzhou city (Figure 1), the capital of Henan Province in the central part of the P.R. China, with a total area of 7507 km 2 as well as a population of 9 878 000 inhabitants. Zhengzhou is one of the National Central Cities in China and serves as the political, economic, technological, and educational center of the province, as well as a major transportation hub in China (http://en.wikipedia .org/wiki/Zhengzhou). The annual average temperature of the city is 14.5°C, and the general terrain trend is tilt from southwest to northeast. The study area is undergoing the accelerating of Chinese agglomeration, economic development, and urban expansion.
The primary dataset in this study is multitemporal Landsat remotely sensed data, with data covering the period from 1986 to 2016, which were acquired from USGS (https:// earthexplorer.usgs.gov/). As a crucial instrument aboard Landsat satellites to collect imagery, the spatial resolution of TM/ETM+ is 30 m for the VIR-NIR band and 60 m/120 m for the TIR band. Landsat data have been widely used for LULCCs or UHI monitoring, while few of the research attempted to track and combine the long-term dynamic of LULCCs and UHI for environment security assessment. The second major dataset used is GlobelLand30, which we utilized as a reference source (http://www.globallandcover.com). The vector data, demographic data, and simulated data in the future are also collected for future stats analysis. A more detailed description about datasets, including captured date, sensor types, and resolution, can be found in Table 1. The data availability statement: all datasets used in this paper including the primary remotely sensed data, boundary vector data, and processed results can be obtained from hyperlink: https://pan.baidu .com/s/1d9GwkEpwryWYgHU_78qV5g, with password: c8rp. Researchers who are interested in this topic can download the data from the above hyperlink or contact the corresponding author to obtain source data to conduct secondary analysis.
The selected remotely sensed data were preprocessed to overcome geometric and atmospheric conditions, distortion, and errors through atmospheric and radiometric correction. To reduce the geometric distortion and radiometric difference, the selected remotely sensed data were preprocessed through geometric and radiometric correction. Specifically, following the previous works in [23], the radiometric calibration procedure is finished using the commercial software (ENVI). In addition, the FLAASH module of ENVI software was used to complete the atmospheric correction.

Methodology
The methodology employed in this research includes three main stages: (1) LULC and UST maps depicted from 1986 to 2026 with the help of SVM, MWA algorithms, and CA-Markov model; (2) analyzing and discussing of spatial distribution and temporal change of LULC and UST with respect to expansion intensity, buffer zone analysis, human settlement transfer, etc.; (3) quantitative and qualitative evaluation relationship between LULC and UST, more details can be found in Figure 2 3.1. Land Use and Land Cover Classification. Land use and land cover (LULC) maps were retrieved from four-time nodes remotely sensed dataset over the metropolitan city of Zhengzhou. In this image classification processing, one of the most state-of-art machine learning algorithm, support vector machine (SVM), was selected to classify research area into five different categories (agricultural land, vegetation, waterbody, bare land, and building land) according to the Classification Criteria for Land Use Status/GB-T21010-2015 and GlobeLand30 standard products [24]. The process for LULC classification using the SVM classifier includes training samples selection, SVM kernel determination, feature vectors inputs, and SVM classifier application. About 100 training samples of each class are selected by an expert in the field of RS; the most robust radial basis function (RBF) was selected as the SVM classifier: where γ is the bias term in kernel function for polynomial and sigmoid kernels, and Landsat spectral bands except thermal infrared are selected as inputs. Classification accuracy will affect subsequent change analysis; in this research, overall accuracy (OA) and kappa coefficient (KC) were chosen as the evaluation criterion. OA and KC are a nonparametric test which can reflect the consistency of labelled value and predicted value [12]. The mathematical model of OA and KC can be expressed as where q is the number of classes, n represents the total number of considered pixel, n ii is the diagonal element of the confusion matrix, n i+ represents the marginal sum of the rows in the confusion matrix, and n +i represents the marginal sum of the columns in the confusion matrix.
3.2. Retrieval of Urban Surface Temperature. By literature review, there are three algorithms, including single-channel algorithm, monowindow algorithms, and radiative transfer equation algorithm, which are widely used for urban surface temperature retrieval from a single-band Landsat thermal band [25]. In this research, Qin et al.'s monowindow algorithm (MWA) was selected for retrieval urban surface temperature from Landsat TM/ETM+ image. The mathematical model of MWA can be expressed as where UHI is the urban surface temperature; T sensor is the brightness temperature, in a normal case, constant a 6 = − 67:355351 and constant b 6 = 0:458606; ε is the ground surface emissivity; τ 6 is the atmospheric transmittance; T a is the effective mean atmospheric temperature; and T 0 is the near-surface air temperature. While like many other urban heat island analysis studies [15,26], we chose a method to obtain TOA spectral radiance and focus on the dynamic evolution of urban surface temperature in spatial and temporal. The standardized method used in this research is as follows: where N i is normalization temperature value of the ith pixel position; UHI i means the urban heat temperature value of the ith pixel position; UHI max and UHI min represent the maximum and minimum value of UHI before normalization.
Considering the actual situation of the research area, the statistical mean-standard deviation method is selected to divide the thermal field into different levels. According to the statistics of temperature, the average and standard deviation of UHI are used as the demarcation point of the division interval. The mean values of average and standard deviation of UHI maps obtained in 1986,1996,2006, and 2016, where T = 0:512, T s = 0:162 were taken as the demarcation point of the thermal interval. The divided interval results are shown in Table 2. 3.3. CA-Markov Model. The CA model [27] has a capacity for spatial and temporal change simulation which can be defined as where f is the local transition rule of the cell, S is a set of cellular states, N is the cellular field, and t + 1 and t represent the start and end time. The Markov model [28] that has been widely used for trend simulation at various scenarios can be expressed as   Journal of Sensors where S t+1 and S t are statuses at time t + 1 and t, respectively, P ij is the transition probability matrix in a state. The CA-Markov model is a combination of the CA model and the Markov model, both of which are dynamic models with discrete states [29]. And when combined with the CA model, the CA-Markov model can overcome the limitation that the Markov model failed to catch the spatial distribution in the future and can be used to simulate spatial-temporal changes.
In this research, LULC and UST maps in 2006 and 2016 were selected as the study years to calculate the transition area matrix using the Markov model. And a standard 5 by 5 contiguity filter was selected for the CA model which means the condition of the future pixel is not only decided by information from the previous state but also considered by corresponding surrounding pixels.

Urban Expansion
Intensity. The rate of urban expansion (RUE) and the intensity of urban expansion (IUE) are the most common methods for describing urban expansion. The RUE describes the annual average area change of the built-up area during the research period, while IUE refers to the proportion of urban land use expansion area of a space unit in the total urban area during the research period. The RUE and IUE only study the quantitative change of urban built-up area at the beginning and end of a certain period, while ignoring the dynamic process of urban growth. In this paper, based on urban growth intensity, the concept of urban expansion intensity index (UEII) was proposed [30], which describes the degree of differentiation of urban expansion in different directions and denotes the growth of the urban areas of a spatial unit as a percentage of the total area of the land unit in the study period [31] and can be used to study the urban growth process more reasonably and accurately. The mathematical model of UEII is as where UEII is the urban expansion intensity index of the spatial unit during periods between t and t + n; UA t+n and UA t are the urban area in the spatial unit at time t and t + n, respectively; T A is the total land area; n is the research period.
The UEII can be divided into different levels according to relevant literature [31] and shown in Table 3.

Urban Gravity Center
Movement. The gravity center model reflects the direction of movement and distance to the center of gravity over time [32,33]. The city's center of gravity reflects the geometric equilibrium of urban space to a certain extent, in which spatial position will be constantly moving during the growth of the urban built-up area. The movement direction of the gravity center reflects which direction of urban growth. To make further analysis of the city's spatial expansion, the urban's gravity center is introduced in this research. The mathematic model of the gravity center can be described as where X t and Y t are the horizontal and vertical coordinates of the gravity center of the t th year, a ti means the area of i th polygon of the t th year, and x i and y i mean the horizontal and vertical gravity center coordinates of the i th polygon. The mathematical equation of gravity center distance can be expressed as follows: where D is the distance of the urban gravity center shifted from the beginning to ending time point; X i and Y i represent the horizontal and vertical urban center gravity coordinates, respectively, in the initial moment of research, while X j and Y j represent corresponding coordinates in the ending moment.

LULC Maps and Urban Expansion
Intensity. There are a total of 15538, 13829, 15991, and 13937 pixels that were selected for training in remotely sensed data captured on 12 August 1986, 15 April 1996, 9 April 2006, and 23 April 2016, respectively. The number of selected test pixels (ground truth points) is 5643 points. The classification maps from 1986 to 2016 are shown in Figure 3; overall accuracy and kappa coefficient evaluation results are shown in Table 4. On the basis of classification maps of 2006 and 2016, the LULC map of 2026 predicted with the CA-Markov model is shown in Figure 4.
According to the above-mentioned urban expansion intensity model, the statistical results of building area classification and simulation in 1986,1996,2006,2016, and 2026 (shown in Figures 3 and 4), and urban growth intensity graded criteria, the level of urban growth intensity is obtained (Table 5).

Levels
Criteria Interval Urban Gravity Center Changes. Based on the gravity center model, the gravity center and changing footprint of Zhengzhou's urban built-up area from 1986 to 2026 were calculated. The gravity centers and the trends of gravity center development are shown in Table 6.
The results showed that the direction and distance of the city's center of gravity movement are quite different in a different period. The center of gravity movement of each period during the study period has the following characteristics: (1) The city center of Zhengzhou moved 9.65 km to the northwest from 1986 to 1996. The dominant driving forces was the rapid development of urban construction in Xingyang, a country-level city now becomes as a district of the metropolis of Zhengzhou, from 1986 to 1996. The development of Jinshui District of Zhengzhou was earlier and faster than other urban areas, which led to the gravity center movement to Jinshui and Xingyang district. (2) The city center moved 15.69 km to the southeast from 1996 to 2006, mainly because with the acceleration of urbanization, the Zhengzhou station was expanded; the Xinzheng international airport, the city around expressway, and some large-scale construction project were gradually put into service, which expanded the surface area of man-made land type in the southeast direction and causing the shift of gravity obviously. (3) The city center moved 7.87 km to the north from 2006 to 2016, mainly due to the large-scale construction of Zhengzhou High-tech Industrial Zone and the operation of Zhengzhou East High-Speed Rail Station operated in 2012, which resulted in a rapid increase in the man-made surface in the north. Meanwhile, the development of old district Erqi city of Zhengzhou has been saturated. Therefore, the development speed of Erqi District was lower than that of Jinshui District during 2006 to 2016, which made the gravity center move to the north in the past 10 years. (4) According to the prediction result of LULC using the CA-Markov model, the gravity center will move 1.14 km westward from 2016 to 2026. With the continuous urbanization process, the urbanization of Zhengzhou has basically become saturated and the expansion speed has slowed down.  Table 7).
According to the above results, in general, the strongest urban growth intensity was in the NEE and SWW directions from 1986 to 2026, which is the dominant direction of urban growth in the study area. Combined with the geographical location of the administrative area of the study area, the strongest urban growth intensity is concentrated in Jinshui District, Erqi District, and Xinmi, mainly because the urbanization    Temporal dimension studies showed that the urban growth intensity in NEE and NWW orientation was relatively strong from 1986 to 1996, which is mainly due to the rapid development of Erqi District and Xingyang District. Urban growth was mainly in the direction of NEE, SWW, and SSW from 1996 to 2006, and during this period, the urban growth pattern was consistent, and the growth of other cities was relatively slow, which made the city's center of gravity move to the southeast. From 2006 to 2016, the urban growth intensity in the NEE, NWW, and SWW directions was relatively strong. During this period, there was vigorously built activities in transportation facilities and high-tech industries in Zhengzhou city, which further strengthened the land development efforts and increased the intensity of urban growth. In the meantime, the land use development status of prefecturelevel cities such as Gongyi and Xinmi was growing rapidly because of an unsaturated state. According to the prediction of LULC results, there will be a similar trend of urban growth between the period of 2016 to 2026 and the period of 2006 to 2016. The urbanization construction of Zhengzhou city has reached a new stage, and urban growth will be relatively stable and mature at that time.

UHI Results and Statistical
Features. The results of UHI from 1986 to 2016 are retrieved using WMA, and the UHI degree of 2026 is simulated using the CA-Markov model. UHI distribution maps and UHI degrees from 1986 to 2026 are shown in Figure 6 and Table 8.
The spatial distribution maps of UHI ( Figure 6) in the study area showed that the low temperature and sub-low temperature zone were mainly located in the Yellow River basin in the northeast of the study area. And the high-temperature zone is mostly located in the urban built-up zone, north of the central part. Based on the statistical results of Figure 7 and Table 8, the rank of different temperature areas has characteristics as that the dominant areas by proportion are medium temperature and sub-high temperature level, followed by high temperature, sub-low temperature, and low-temperature areas.

Spatial-Temporal Changes of Urban Thermal Environment.
According to the statistical results of the UHI degree, the dynamic curve of the UHI level is obtained (Figure 8). The changes and proportions of UHI in the study area between 1986 to 1996, 1996 to 2006, and 2006 to 2016 were calculated (shown in Tables 9-12). The methodologies used for LULCC analysis were adopted to analyze spatial-temporal changes of the urban thermal environment (shown in Figures 7-9).
The dynamic curve of the UHI trend result (Figure 8) showed that the area of the low-temperature zone and sublow temperature zone is shrunken during the period from 1986 to 2026. The dominated UHI zones during this period are medium, sub-high, and high-temperature zone. The area of the low-temperature zone has a downward trend, and the area of the sub-low temperature zone experiences a slight increase and then a continuous decline. The change trend of the sub-high temperature zone is opposite to that of the sub-low temperature zone. The middle-temperature zone and the high-temperature zone are both in a wave-like rising trend, while the trend of medium temperature zone is opposite to that of the high-temperature zone.
Combined with the distribution map of the urban thermal environment in the study area for the past 40 years, a long-term UHI comparison was carried out (Figure 9). Between 1986 and 1996, the thermal environment change process of Zhengzhou city was complicated, with a large change range, fast change dynamic. The area of sub-low temperature and medium temperature zone increased; the area of low-temperature, sub-high temperature, and hightemperature zone decreased. The dynamic of the medium temperature zone changed the most, that is, more than 10%. Between 1996 and 2006, the thermal environment of Zhengzhou has a tendency to change to a high-temperature area. The area of medium and low temperature decreased, the area of sub-high temperature increased, and the hightemperature area changed the most, that is, more than 85%. While the change range of other degrade is relatively small, in particular, the area of the low-temperature zone is almost unchanged. From 2006 to 2016, the trend of the thermal environment was similar to that between 1986 and 1996. The area occupied by extreme temperature level decreased; the area occupied by intermediate temperature level increased. And the urban thermal environment had a tendency to concentrate to medium and high-temperature levels. According to the prediction result, between 2016 and 2026, the overall thermal environment of the city has a small change and is in a state of dynamic equilibrium. It is predicted that by 2026, urban surface temperature will further evolve into hightemperature areas, and the problem of the UHI effect will become more prominent.      Figure 10). Profiles of LST in a different period showed that (1) there are significant jagged jumps in the EW and NS profiles of LST in each period. The urban central area and the suburbs are characterized by "upward convexity" and "depressed." This is mainly due to the complex structure of the underlying surface of the urban surface, and the change of the underlying surface characteristics in a small area makes the surface temperature abrupt; (2) in this study, LST in the central area of the city was generally higher, and there was obvious "bumps" with uneven "peak" and "low valley" morphological features. In the suburbs, the LST was low, while the "jumping" phenomenon with a large jump and fast frequency of change is more significant than that of the city center. The cause of that related to the underlying surface of the central area is mainly composed of steel, cement, and masonry, and the structure is relatively simple, while the underlying of suburb mainly consists of green land, water area, soil, cement, and masonry, and the structure is relatively complicated; (3) the comparison of various characteristics of EW and NS LST profiles showed that the EW direction changes rapidly and complex than that of NS direction. This phenomenon showed that the underlying surface of EW section of Zhengzhou has a more diverse structure type than NS direction, and changes are more complicated; (4) the perspective of the long-term changes in the distribution pattern of LST, EW, and NS profiles showed that with the continuous development of urbanization construction, the overall average temperature of the city was gradually increasing, while the frequency and jump range of urban surface temperature profile had a decreasing trend. This is mainly because the urbanization construction of Zhengzhou city is gradually improved, the characteristic structure of the underlying surface of the city tends to be stable, and the temperature field structure of the urban thermal environment tends to be simplified.
Buffer analysis is an important spatial analysis method used to determine the proximity of research elements in GIS. In this research, the buffer analysis method was adopted to study the urban thermal environment. With the support of this methodology, the relationship between LST and the location of the city center from 1986 to 2026 was discussed, the spatial characteristics of LST distribution within a certain distance were described, and the pattern of thermal environment space of Zhengzhou city was analyzed. Four circular       10 Journal of Sensors buffers with a distance of 5 km with the Erqi Memorial Tower taken as the center of the circle was drawn ( Figure 11). The LST grade vector data of 1986, 1996, 2006, 2016, and 2026 and buffer vector data were superimposed and analyzed to obtain LST distribution results of five different distances from the center of the circle (Figure 12). The proportion area of different grades in the different buffer zone (Figure 12) indicates that during the 40 years' urbanization in Zhengzhou, the area of low-temperature zone gradually decreased and the area of intermediate grade temperature zone experienced a process of increasing first and then decreasing; the area of higher temperature zone increased continuously. Therefore, the pattern of urban thermal environment in Zhengzhou city showed a trend of agglomeration to an intermediate level firstly and then change to the extremely high temperature.
By longitudinally comparing the temperature levels in the buffer zone, the results showed that in the 0 km to 10 km buffer zone, the temperature of each grade changes greatly, the area of low-temperature zone gradually decreases, and the area of high-temperature zone continues to increase rapidly. The area ratio of high temperature is more than 80% in 2016, and this number will close to 90% in 2026 by prediction. In the buffer zone of 10 km to 20 km, the area of lowtemperature area is generally declining, the area of medium temperature zone is stable, and the area of high-temperature zone is changed from the least to the most. It is predicted that by 2026, the area of the lower temperature zone will be close to 10 km 2 . In the 20 km to 30 km buffer zone, the area of low temperature and sub-low temperature continues to decline, the area of medium temperature zone and subhigh temperature zone is in a wave-like state, and the area of high-temperature zone increases slightly. In the 30 km to 40 km buffer zone, the medium and sub-high temperature zones are the main temperature grades. The proportion of high-temperature zone increased at the beginning and then stabled at 20%; the proportion of low-temperature zone is small and in a reduced state. The above phenomenon is mainly due to the continuous acceleration of the urbanization process in Zhengzhou city, the continuous expansion

11
Journal of Sensors of the human-made area around Erqi Tower, and the increasing surface area of artificial land, which lead to continuous temperature rising in the high buffer zone, and the lowtemperature zone reduced gradually.
The horizontal comparison (in the 10 km, 20 km, 30 km, and 40 km buffer zone) of the temperature levels showed that the change trend of low-temperature zone increased and then decreased, which displayed as a "n" shape change trend. The area of sub-low temperature zone increased gradually during the period of 2006 to 2026. The overall change range of medium temperature zone during the period from 1986 to 1996 is small, while it is increasing gradually during the period of 2006 to 2026. The sub-high temperature zone occupied a larger proportion in the 10 km buffer zone than

Relationship between LULC and UHI.
In this subsection, the relationship between LULC and UHI is analyzed using three quantitative indices such as normalized difference vegetation index (NDVI), normalized difference building index (NDBI), and normalized difference waterbody index (NDWI) which indicate relative important LULC types vegetation, built-up area, and waterbody. NDVI is a standardized way to measure vegetation, which quantifies vegetation by measuring the difference between red (R) and near-infrared (NIR) spectral reflectance value from remotely sensed data. NDVI can be calculated in its formula NDVI = ð NIR − RÞ/ðNIR + RÞ. NDBI is another solution for easily calculating of the built-up area because it is simple, rapid, and accurate in urban area mapping. NDBI can be calculated using the formula as NDBI = ðSWIR − NIRÞ / ðSWIR + NIRÞ, where SWIR and NIR represent the spectral reflectance value of shortwave infrared and NIR band. The NDWI, which can be calculated using green and NIR spectral reflectance value using its formula NDWI = ðGREEN − NIRÞ /ðGREEN + NIRÞ, is most appropriate for waterbody mapping. In this research, there are 30 sets of NDVI, NDBI, and MNDWI index and corresponding UST data were randomly The regression equations (15), (16), and (17) showed that there was a significant negative correlation between NDVI, NDWI, and LST; the correlation coefficients are 0.835 and 0.868, respectively. Waterbody has a more obvious effect on relieving LST than vegetation. There was a significant positive correlation between NDBI and LST, with a correlation coefficient of 0.821.

Analysis and Discussion
The urban expansion intensity level (Table 5) from 1986 to 2026 indicates that the urban growth intensity maintains an increase during this period in general. There was a medium-speed growth stage from 1986 to 1996 and a rapid growth from 1996 to 2006. And according to the predicted results of LULC maps, there will be a slight slow down to some extent in the next 10 years, while the urban growth is still at a medium speed.
What is the relationship between LULC type and UHI? In order to discover this relationship and to effectively and   (Tables 13-17). The spatial-temporal change of thermal environment from 1986 to 2026 showed that the proportion of LULC types in a temperature grade changed significantly. For example, the built-up area in the high-temperature zone has rapidly increased from less than 10% in 1986 to more than 45% in 2016. In the low-temperature zone, the proportion of water has been close to 65% from the beginning of the study period and has continued to drop to less than 30% in 2016, and it is predicted that the proportion will further decrease significantly by 2026. Priyankara et al.'s research [34] demonstrates similar findings that mean LST has a strong significant positive relationship with a fraction of impervious surface and persistent impervious surface, while a strong negative relationship with a fraction of forest surface and new added impervious surface. Priyankara et al. suggested that more vegetation areas are recommended in both horizontal and vertical directions to reduce the UHI effect.
The change of proportion of the same LULC type in different temperature grades from 1986 to 2026 showed that the changes of the built-up area and bare land area are similar, that is, in the early stage of research, these two land cover types are mainly concentrated in medium temperature zone. With the expansion of the urban area, the proportion of the high-temperature zone increased. When combining the spatial distribution of built-up and bare land in the study area, it fully demonstrates the significant positive impact and contribution of these two LULC types on urban thermal environment. The distribution of vegetation is mainly in a medium level of UHI grades, which indicates that vegetation plays an important role in balancing surface temperature. Waterbody is distributed in a low-temperature zone, and proportion was gradually decreased in the low-, medium-, and high-temperature zone. With the acceleration of urbanization process, the proportion of waterbody in a lowtemperature zone will be less than 8% in 2026 by prediction, but it still has an obvious advantage over other land cover types showing that the waterbody is indispensable for reducing UHI effect and maintaining the balance and stability of urban thermal environment. LULC types of distribution on all temperature grades were discussed here, while the contribution of other factors might be differentiated. Ranagalage et al.'s research [35] revealed that mean LST was positively correlated with the increase of fraction ratio of building area and forest area and with the decrease of fraction ratio of agricultural and forest area. Building density is a crucial element in increasing LST.

Conclusions
In this research, comprehensive research of LULCCs on urban heat environment assessment was performed using the RS and GIS spatial technique. The spatial and temporal changes of urbanization of Zhengzhou city from 1986 to 2026 were analyzed, and conclusions can be drawn as follows.
Land use and land cover changes of 40 years in Zhengzhou city have been studied and analyzed. In the past forty     A combination analysis of natural and thermal environmental has been yielded. During the period from 1986 to 2026, LST in the study area is distributed mainly in medium, sub-high, and high-temperature zone. The thermal environment change process in Zhengzhou is relatively complicated, and the dynamics of spatial-temporal change were dramatical. In the early stage of the study period, the temperature grade trends to medium zone, the middle trends to hightemperature zone, and the later stage has a tendency to change to the medium zone. The temperature changes in the east-west direction were faster than that in the northsouth direction. There is a significant correlation between vegetation, water and urban surface temperature, and LST, and a positive correlation between the built-up area and LST.
The apparent drawback of the paper is that our study is limited to daytime UHI due to the limit of the Landsat dataset. Therefore, future works tend to integrate the advantages of Landsat and nighttime light datasets and then extend our study to nighttime.