Prediction of Urban Thermal Environment Based on Multi-Dimensional Nature and Urban Form Factors

: The urban thermal environment is affected by multiple urban form and natural environment factors; research on the accurate prediction of the urban thermal environment, considering the interaction among different urban environmental factors, is still lacking. The development of a machine learning model provides a good means of solving complex problems. This study aims to clarify the relationship between urban environmental variables and the urban thermal environment through high-precision machine learning models as well as provide scenarios of future urban thermal environment developments. We deﬁned an urban thermal environment index (UTEI), considering twelve urban form and natural indicators sourced from the remote sensing data of 150 cities in the Jing-Jin-Ji region from 2000 to 2015. We achieved accurate predictions of UTEI through training a gradient-boosted regression trees model. By unpacking the model, we found that the contribution rate of elevation (ELEV) was the highest. Among all the urban form indicators, the elongation index (ELONG), urban population (POP), nighttime light intensity (NLI), urban area size (AREA), and urban shape index (SHAPE) also had high contributions. We set up ﬁve scenarios to simulate the possible impact of different urban form factors on the overall urban thermal environment quality in the region. Under extremely deteriorated patterns that do not control urban expansion and vegetation reduction, the average UTEI could be as high as 0.55–0.76 ◦ C in summer and 0.24–0.29 ◦ C in winter, yet in the extremely optimized situation, UTEI decreased by 0.69 ◦ C in summer and 0.56 ◦ C in winter. Results showed that better urban form improves the quality of urban environments and can provide important insights for urban planners to mitigate urban heat island problems. elevation (ELEV) data were retrieved from the Shuttle Radar Topography Mission (SRTM; 30 m resolution) (https://earthdata.nasa.gov/, accessed on 30 March 2019). Wind speed (WIN), the precipitation rate (PRE), air relative humidity (RHU), and solar radiation (RADI) were sourced from the China Meteorological Forcing Dataset (CMFD) [36]. The CMFD is a high spatial–temporal resolution gridded near-surface meteorological dataset developed speciﬁcally for studies of land surface processes in China, with a temporal resolution of three hours and a spatial resolution of 0.1 ◦ . rate (PRE),


Introduction
Cities are growing at an alarming rate. As the United Nations predicts, nearly 66% of the world's population will live in cities by 2050 [1]. Therefore, the improvement of urban thermal environment quality, which has an important impact on human health and energy consumption, is more and more significant in order to promote sustainable development and improve human health and well-being [2,3]. The changes in surface heat flux, moisture flux, and air circulation caused by thermal environment change have a negative impact on water resources, air quality, and human health [2,4,5]. How to effectively evaluate the negative impact of urban development on the thermal environment in the current century has become a key issue on the recent global agenda, according to The 2030 Agenda for Sustainable Development [6] and The New Urban Agenda [7].
The deterioration of the urban thermal environment, which is often measured by the magnitude of the urban heat island effect and the rise of land surface temperature, is affected by complicated urbanization processes and natural environment factors [8][9][10][11][12][13]. In recent years, more and more studies have shown that urban form and natural factors have multi-path effects on urban thermal environment quality, and their relationship has been proved to be nonlinear in previous studies [14]. For example, some studies have shown that the effect of city size on the intensity of the urban heat island effect shows a marginal decreasing effect [15,16]. Urban layout and population density have a strong interaction with urban air pollution or other environmental problems [17], and their contributions to the urban thermal environment might demonstrate unique change patterns along the urban development gradient [14].
Although some scholars have done research in the related field, some have only focused on a few indicators that influence the urban thermal environment, and they lack consideration of other dimensions of urban environmental factors [18,19]. Otherwise, many studies have adopted the ordinary least squares (OLS) model and other traditional regression models that cannot accurately detect the nonlinear influence of urban form on the thermal environment and the complex relationship between urban form and natural factors [20][21][22].
Recently, machine learning models have been successfully applied for scenario prediction and pattern recognition in various real-world situations [23][24][25][26] since they are useful in identifying the underlying relationship between different sources of data [27,28]. Antonio et al. [29] found that the modeling of environmental noise in urban environments was a nonlinear problem and proved that machine learning regression methods widely outperformed multiple linear regression models. There are also other cases of using machine learning models to solve urban environmental problems. Alavipanah et al. [30] used a boosted regression tree (BRT) to explore the influence of 2D and 3D building indicators and vegetation coverage on the thermal environment. Zhang et al. [31] predicted future changes in land use/land cover and land surface temperature based on an artificial neural network (ANN) model. Sun et al. [32] found that urban density exhibited significant positive effects on LST in Ningbo city based on random forest regression (RF). Machine learning models can better deal with identifying the complex relationship between urban natural and social environmental elements and can help to accurately predict the changes in the urban thermal environment under different urban development scenarios. Nonetheless, using a machine learning model to explore the influencing factors of the urban thermal environment is still in the early stages and comparatively lacking.
To fill the gaps, the objective of this study is to analyze the multivariate correlation between urban indicators and the urban thermal environment through a high-precision prediction model and simulate the possible development changes under typical urbanization scenarios. We used urban form indicators (UFIs), including urban area size (AREA), nighttime light intensity (NLI), urban population (POP), the urban shape index (SHAPE), the elongation index (ELONG), the urban aggregation index (AI), and urban vegetation cover (UVC), combining natural factors such as urban altitude (ELEV), wind speed (WIN), air relative humidity (RHU), precipitation (PREC), and solar radiation intensity (RADI) in the Jing-Jin-Ji region from 2000 to 2015. Six linear regression and nonlinear machine learning models were constructed for model comparison. We finally verified the optimal prediction performance of a gradient-boosted regression trees model and calculated the relative importance of each factor. According to historical trends, we built five scenarios and simulated the possible impact of different urban form indicators on the overall urban thermal environment in the region by calculating the PDF value. The innovation of this study is to provide deeper insights on how to alleviate the UHI effect and improve the urban thermal environment by optimizing the local urban form.

Study Area
The study area is located in north China and is commonly known as the Jing-Jin-Ji region ( Figure 1). It covers an area of 218,000 km 2 and is one of the fastest developing urban agglomerations in China. During the study period of 2000 to 2015, the Jing-Jin-Ji region experienced rapid population growth and urbanization processes, causing severe The study area is located in north China and is commonly known as the Jing-Jin-Ji region ( Figure 1). It covers an area of 218,000 km 2 and is one of the fastest developing urban agglomerations in China. During the study period of 2000 to 2015, the Jing-Jin-Ji region experienced rapid population growth and urbanization processes, causing severe environmental deterioration and urban thermal environment problems [33]. The region mainly belongs to the temperate monsoon climate zone [34]; most areas have four distinct seasons, which are characterized by hot and rainy summers, dry and rainy springs, sunny autumns with sand storms, and cold winters. The region is comprised of a large number of cities, and we extracted 150 urban patches based on land use products in 2015.

Urban Thermal Environment Index
We defined the urban thermal environment index (UTEI) to characterize the thermal environment of the study area. It was calculated by the average of the land surface temperature (LST). LST has the advantage of being spatially explicit, while the commonly used meteorological station data cannot accurately reflect the temperature difference between each city region [35]. Therefore, we used remotely sensed land surface temperature data from the MOD11A2 dataset (http://www.gscloud.cn, accessed on 15 March 2019) with 1 km resolution as the source of land surface temperature. The MOD11A2 Version 6 product provides an average 8-day per-pixel LST, so we calculated monthly and yearly UTEI using Equation 1 for the years 2000, 2005, 2010, and 2015: where xm,n stands for a single pixel value in the LST grid, and m×n calculates the number of pixels in the region. The mean LSTs for both summer and winter were summarized by overlaying the regulatory observation grid unit layers and are specified as the dependent variables for further analysis.

Urban Thermal Environment Index
We defined the urban thermal environment index (UTEI) to characterize the thermal environment of the study area. It was calculated by the average of the land surface temperature (LST). LST has the advantage of being spatially explicit, while the commonly used meteorological station data cannot accurately reflect the temperature difference between each city region [35]. Therefore, we used remotely sensed land surface temperature data from the MOD11A2 dataset (http://www.gscloud.cn, accessed on 15 March 2019) with 1 km resolution as the source of land surface temperature. The MOD11A2 Version 6 product provides an average 8-day per-pixel LST, so we calculated monthly and yearly UTEI using Equation where x m,n stands for a single pixel value in the LST grid, and m × n calculates the number of pixels in the region. The mean LSTs for both summer and winter were summarized by overlaying the regulatory observation grid unit layers and are specified as the dependent variables for further analysis.

Natural Indicators
We selected five commonly used natural indicators [14], including elevation (ELEV), the precipitation rate (PRE), solar radiation (RADI), air relative humidity (RHU), and wind speed (WIN). Illustration of the indicators' major impact pathways on UTEI is represented in Figure 2, and the physical definitions of these natural indicators are listed in Table 1. In order to achieve better spatiotemporal resolution and accuracy to reflect the urban natural environment, we also used remotely sensed land surface data as data sources. The elevation (ELEV) data were retrieved from the Shuttle Radar Topography Mission (SRTM; 30 m resolution) (https://earthdata.nasa.gov/, accessed on 30 March 2019). Wind speed (WIN), the precipitation rate (PRE), air relative humidity (RHU), and solar radiation (RADI) were sourced from the China Meteorological Forcing Dataset (CMFD) [36]. The CMFD is a high spatial-temporal resolution gridded near-surface meteorological dataset developed specifically for studies of land surface processes in China, with a temporal resolution of three hours and a spatial resolution of 0.1 • . the precipitation rate (PRE), solar radiation (RADI), air relative humidity (RHU), and wind speed (WIN). Illustration of the indicators' major impact pathways on UTEI is represented in Figure 2, and the physical definitions of these natural indicators are listed in Table 1. In order to achieve better spatiotemporal resolution and accuracy to reflect the urban natural environment, we also used remotely sensed land surface data as data sources. The elevation (ELEV) data were retrieved from the Shuttle Radar Topography Mission (SRTM; 30 m resolution) (https://earthdata.nasa.gov/, accessed on 30 March 2019). Wind speed (WIN), the precipitation rate (PRE), air relative humidity (RHU), and solar radiation (RADI) were sourced from the China Meteorological Forcing Dataset (CMFD) [36]. The CMFD is a high spatial-temporal resolution gridded near-surface meteorological dataset developed specifically for studies of land surface processes in China, with a temporal resolution of three hours and a spatial resolution of 0.1°.  Solar radiation RADI W*m −2 3-hourly mean (from −1.5 hr to +1.5 hr) surface downward shortwave radiation.

Urban Form Indicators
Urban form is a multi-dimensional concept containing urban structures and geometric forms as well as urban development factors. We used three widely applied landscape indexes to represent the urban geometric form, including the urban shape index (SHAPE), the urban aggregation index (AI), and the elongation index (ELONG). These indicators were calculated in fragstats 4.2, and the calculation methods are shown in Table 2.

Urban Form Indicators
Urban form is a multi-dimensional concept containing urban structures and geometric forms as well as urban development factors. We used three widely applied landscape indexes to represent the urban geometric form, including the urban shape index (SHAPE), the urban aggregation index (AI), and the elongation index (ELONG). These indicators were calculated in fragstats 4.2, and the calculation methods are shown in Table 2.
We used night light data to represent the level of urban development. The nighttime light intensity (NLI) was from the Defense Meteorological Satellite Program Operational Linescan System (DMSP-OLS), with an approximately 1 km resolution. We used the ridge line method to correct the nighttime light data time series to obtain consistent time series data and then obtained the final urban nighttime night index (NLI) [37]. Due to the lack of night light data in 2015, we used the data from 2013 instead. Urban area (AREA) indicates urban size and was obtained from the remote sensing monitoring data of land use status in China. It was downloaded from the Resource and Environment Data Cloud Platform (http://www.resdc.cn/DataList.aspx, accessed on 1 May 2020), with 30 m resolution. Urban population density was used to measure urban density, which was obtained from the Center for International Earth Science Information Network at Columbia University (http://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-density-rev10, accessed on 12 March 2020), with a resolution of 1 km. The population density data and area data were log-transformed to make their distributions more normal and to increase the predictive ability of some linear models for comparison [38]. We also achieved the monthly normalized differential vegetation index (NDVI) from MODIS datasets (MOD13Q1) of July and December, with 250 m spatial resolution to represent urban vegetation cover (UVC). Table 2. Computational method and description of geometric urban form indicators.

Name Parameter Description
Shape Index SHAPE = 1 when the patch is square and increases without limit as the patch shape becomes more irregular.
Aggregation Index (AI) AI = g ii max→g ii (100) g ii = number of like adjacencies (joins) between pixels of patch type (class) i based on the single-count method. max-> g ii = maximum number of like adjacencies (joins) between pixels of patch type (class) i (see below) based on the single-count method.
AI increases as the focal patch type is increasingly aggregated and equals 100 when the patch type is maximally aggregated into a single, compact patch.

Methods
Based on the natural and urban form indicators obtained in Part 2.2, our research methods are mainly divided into three steps. The first step is to screen and establish the optimal prediction model. In this process, we divided the sample dataset set into a training set, a validation set, and a test set for supervised training and improved the generalization performance of the model by regularization methods. Then, we obtained the best model by comparing the results of cross-validation. The second step is to analyze the relative importance of urban form and natural indicators based on our interpretable machine learning model. The third step is to simulate the possible impacts of urban form changes in the next fifteen years by setting up different scenarios through the statistical data of historical changes. An overview of all processes is shown in Figure 3. e 2022, 13, x FOR PEER REVIEW 6 of 14

Linear Regression Models
Linear regression is the simplest and most widely used statistical technique for predictive modeling. We implemented three widely used linear regression models: the lasso regression model, the Bayesian ridge regression model, and the elastic net regression model. The lasso regression model is a regression analysis method that performs both variable selection and regularization in order to enhance the prediction accuracy and interpretability of the statistical model it produces [39]. The Bayesian regression model is an approach to linear regression in which statistical analysis is undertaken within the context of Bayesian inference [40]. When the regression model has errors that have a normal distribution, and if a particular form of the prior distribution is assumed, explicit results are available for the posterior probability distributions of the model's parameters. The elastic net regression model is a regularized regression method that linearly combines both penalties L1 and L2 of the lasso and ridge regression methods [41]. It is useful when there are multiple correlated features.

Nonlinear Regression Models
We also used three nonlinear models for supervised training: the random forests (RF) model, the Adaboost model, and gradient-boosted regression trees (GBRT) model.
The random forest regression model is an ensemble learning method for classification, regression, and other tasks that operate by constructing a multitude of decision trees at training time and outputting the class that is the mode of the classes (classification) or the mean prediction (regression) of the individual trees [42][43][44]. In prediction tasks, each decision tree is constructed using the best split for each node among a subset of predictors randomly chosen at that node. In the end, a simple majority vote is taken for the prediction.
Adaboost (short for adaptive boosting) is a popular boosting algorithm introduced to fit a sequence of weak learners on repeatedly modified versions of the data [45]. It is sensitive to noisy data and outliers. In some problems, it can be less susceptible to the overfitting problem than other learning algorithms.
The gradient-boosted regression trees (GBRT) model is a generalization of boosting to arbitrary differentiable loss functions and is an accurate and effective off-the-shelf procedure that can be used for regression problems. This model has been used in a variety of areas, with advantages that include the ability to find nonlinear transformations, the abil-

Linear Regression Models
Linear regression is the simplest and most widely used statistical technique for predictive modeling. We implemented three widely used linear regression models: the lasso regression model, the Bayesian ridge regression model, and the elastic net regression model. The lasso regression model is a regression analysis method that performs both variable selection and regularization in order to enhance the prediction accuracy and interpretability of the statistical model it produces [39]. The Bayesian regression model is an approach to linear regression in which statistical analysis is undertaken within the context of Bayesian inference [40]. When the regression model has errors that have a normal distribution, and if a particular form of the prior distribution is assumed, explicit results are available for the posterior probability distributions of the model's parameters. The elastic net regression model is a regularized regression method that linearly combines both penalties L1 and L2 of the lasso and ridge regression methods [41]. It is useful when there are multiple correlated features.

Nonlinear Regression Models
We also used three nonlinear models for supervised training: the random forests (RF) model, the Adaboost model, and gradient-boosted regression trees (GBRT) model.
The random forest regression model is an ensemble learning method for classification, regression, and other tasks that operate by constructing a multitude of decision trees at training time and outputting the class that is the mode of the classes (classification) or the mean prediction (regression) of the individual trees [42][43][44]. In prediction tasks, each decision tree is constructed using the best split for each node among a subset of predictors randomly chosen at that node. In the end, a simple majority vote is taken for the prediction.
Adaboost (short for adaptive boosting) is a popular boosting algorithm introduced to fit a sequence of weak learners on repeatedly modified versions of the data [45]. It is sensitive to noisy data and outliers. In some problems, it can be less susceptible to the overfitting problem than other learning algorithms.
The gradient-boosted regression trees (GBRT) model is a generalization of boosting to arbitrary differentiable loss functions and is an accurate and effective off-the-shelf procedure that can be used for regression problems. This model has been used in a variety of areas, with advantages that include the ability to find nonlinear transformations, the ability to handle skewed variables without requiring transformations, computational robustness, and high scalability [23,46,47]. This method progressively fits the base learner to the current pseudo-residual by the least squares method in every iteration. Initially, the average Y value is used as an assumption for predicting all of the observations, which is similar to fitting a linear regression model.
We first adopted a grid search for two parameters (number of estimators and maxleaf nodes) to make the GBRT model have the least mean square error. Then, in order to effectively overcome the overfitting problem of the model, we adopted two regularization methods, shrinkage and subsampling. The shrinkage strategy was proposed by Friedman to introduce the learning rate to scale the step length of the gradient descent procedure [48]. Empirical evidence suggests that small values of the learning rate favor better test errors [49]. The subsampling strategy allowed us to combine gradient boosting with bootstrap averaging (bagging). For each iteration, the base classifier was trained on the basis of a fraction subsample of the available training data. Root mean squared error (RMSE), mean absolute error (MAE), coefficient of determination (R 2 ), and explained variance (EV) were used to evaluate and compare the nonlinear regression models.

Calculation Method of Feature Importance
We further calculated the contribution of each feature to the prediction results based on the trained GBTR model. In the model, since the decrease in impurity from splitting them to create a normalized estimate of the predictive power of that feature was combined with the fraction of samples a feature contributes to, we calculated the contribution of each variable to the reduction of impurity, and the common contribution of several variables could be obtained in the scikit-learn framework. The equation is given as: where p(t) is the proportion of samples reaching t and where j t denotes the identifier of the variable used for splitting node t. Here, we used the Gini index as the impurity function, which is known as the Gini importance [50].

Set-Up of Simulation Scenarios
We adopted a method to determine scenario settings by statistical analysis of historical changes, which is widely used in spatial-temporal predictive simulation [51]. We built five scenarios of UFI projection to 2030 based on the UFI in 2015 and the historical urban UFI change from 2000 to 2015. We first calculated the UTEI change every five years using the dataset of the 150 cities and fit the probability density function (PDF) of the fiveyear changes. Then, we counted the data distribution of five-year changes in all urban morphological indicators, including AREA, NLI, POP, NDVI, SHAPE, AI, and ELONG. These UFIs can reflect the changes in urban characteristics brought about by urbanization, and they are also aspects of urban planning that need to be focused on. We defined the indicators that increase the intensity of the heat island effect as positive indicators and defined other indicators as negative indicators.
From the PDF of each UFI change, we drew the extremely low (10%), low (25%), medium (50%), high (75%), and extremely high (90%) annual UFI change. We defined 90% of the change in positive indicators and 10% of the change in negative indicators as extreme deterioration scenarios (EDs). In this changing situation, this means that the urban form in the next five years will advance in a direction that is not conducive to the improvement of thermal environment quality. Similarly, we took the first 75% of the change of positive indicators and the first 25% of the change of negative indicators as subdeterioration scenarios (SDs), 50% of the change of positive indicators and 50% of the change of negative indicators as keeping the existing trend scenarios (KTs), 25% of the change of positive indicators and 75% of the change of negative indicators as slight optimization scenarios (SOs), and 10% of the change of positive indicators and 90% of the change of negative indicators as an extreme optimization scenario (EO). Note that the optimization and deterioration we mentioned here are only for the urban thermal environment.

Model Comparison and Evaluation Results
Optimization of the GBRT model parameters was obtained by grid searching; the number of estimators was best set to 944, while max depth was best set to 40 (the result is shown in Figure 4a). The effect of regularization strategies to prevent overfitting in the test dataset is shown in Figure 4b. As depicted, when shrinkage was not used, the yellow curve decreased the fastest, which indicated that better deviation can be achieved in fewer iterations, but with the increase of iterations, the deviation will increase to an extent. By contrast, the best result was that only shrinkage was used, of which the deviation of the test set was significantly lower than that of other methods. Therefore, in the final prediction model, we set the learning rate of the model to 0.1.
Atmosphere 2022, 13, x FOR PEER REVIEW 8 of 14 negative indicators as an extreme optimization scenario (EO). Note that the optimization and deterioration we mentioned here are only for the urban thermal environment.

Model Comparison and Evaluation Results
Optimization of the GBRT model parameters was obtained by grid searching; the number of estimators was best set to 944, while max depth was best set to 40 (the result is shown in Figure 4a). The effect of regularization strategies to prevent overfitting in the test dataset is shown in Figure 4b. As depicted, when shrinkage was not used, the yellow curve decreased the fastest, which indicated that better deviation can be achieved in fewer iterations, but with the increase of iterations, the deviation will increase to an extent. By contrast, the best result was that only shrinkage was used, of which the deviation of the test set was significantly lower than that of other methods. Therefore, in the final prediction model, we set the learning rate of the model to 0.1. All six models were firstly trained with default parameter data, and the results of the five-fold cross-validation are shown in Table 3, including the resultant R 2 , EV, MAE, and RMSE for the validating and testing processes. All the urban form indicators and natural indicators were considered in the model estimation process.
All the nonlinear regression models had higher R 2 and lower RMSE values than the three linear regression models in both summer and winter. Regarding the values of RMSE, all the models had slightly higher performances in winter than in summer, which seemed to relate to the much larger dynamic range of summer UTEI than winter UTEI [32].
The GBRT model explained 89% (RMSE = 0.77 °C, MAE = 0.61) of the UTEI in summer and captured 89% (RMSE = 1.11 °C, MAE = 0.83) of the UTEI in winter. Among all the six models, the prediction performance of GBRT was the most stable in both summer and winter. However, the results of other methods in cross-validation results were not so good, with the R 2 being low or fluctuating sharply. All six models were firstly trained with default parameter data, and the results of the five-fold cross-validation are shown in Table 3, including the resultant R 2 , EV, MAE, and RMSE for the validating and testing processes. All the urban form indicators and natural indicators were considered in the model estimation process.
All the nonlinear regression models had higher R 2 and lower RMSE values than the three linear regression models in both summer and winter. Regarding the values of RMSE, all the models had slightly higher performances in winter than in summer, which seemed to relate to the much larger dynamic range of summer UTEI than winter UTEI [32].
The GBRT model explained 89% (RMSE = 0.77 • C, MAE = 0.61) of the UTEI in summer and captured 89% (RMSE = 1.11 • C, MAE = 0.83) of the UTEI in winter. Among all the six models, the prediction performance of GBRT was the most stable in both summer and winter. However, the results of other methods in cross-validation results were not so good, with the R 2 being low or fluctuating sharply.

Relative Importance of Urban Form and Natural Factors
The contribution of all factors to predictive ability was obtained by factor importance analysis, as shown in Figure 5. The results showed that the order of contribution of all factors was similar in summer and winter, only slightly different. The contribution rate of Atmosphere 2022, 13, 1493 9 of 14 ELEV was the highest in both seasons (about 21% in both summer and winter). Relative humidity and solar radiation also had high contribution rates (about 10%). In addition to natural factors, urban form factors also played an important role. Among them, POP and ELONG had the highest contribution rates; their contribution rates reached 7-10%. The contribution ratio of other UFIs was about one order of magnitude with climate factors, mostly within 10%. Among them, AI, SHAPE, and AREA contributed less but also reached 3-6%, which was close to that of UVC. Furthermore, by comparing the contribution rates of each season in a year, we found that natural factors contributed more in summer while urban form factors contributed more in winter.

Relative Importance of Urban Form and Natural Factors
The contribution of all factors to predictive ability was obtained by factor importance analysis, as shown in Figure 5. The results showed that the order of contribution of all factors was similar in summer and winter, only slightly different. The contribution rate of ELEV was the highest in both seasons (about 21% in both summer and winter). Relative humidity and solar radiation also had high contribution rates (about 10%). In addition to natural factors, urban form factors also played an important role. Among them, POP and ELONG had the highest contribution rates; their contribution rates reached 7-10%. The contribution ratio of other UFIs was about one order of magnitude with climate factors, mostly within 10%. Among them, AI, SHAPE, and AREA contributed less but also reached 3-6%, which was close to that of UVC. Furthermore, by comparing the contribution rates of each season in a year, we found that natural factors contributed more in summer while urban form factors contributed more in winter.

Scenario Simulation Results
Based on previous studies [15,19,[52][53], we took AREA, NLI, POP, SHAPE, and AI as positive indicators because their increases played enhanced roles in the UTEI. We used UVC and ELONG as negative indicators. According to our scenario setting method, we noted the change of each indicator over the course of fifteen years in the future under each scenario (Table 4). Taking 2015 as the baseline, we further calculated the presumptive value of each urban form index in 2030. In all scenarios, we set the value of urban form factors beyond the reasonable range of the urban form index as the critical value. Based

Scenario Simulation Results
Based on previous studies [15,19,52,53], we took AREA, NLI, POP, SHAPE, and AI as positive indicators because their increases played enhanced roles in the UTEI. We used UVC and ELONG as negative indicators. According to our scenario setting method, we noted the change of each indicator over the course of fifteen years in the future under each scenario (Table 4). Taking 2015 as the baseline, we further calculated the presumptive value of each urban form index in 2030. In all scenarios, we set the value of urban form factors beyond the reasonable range of the urban form index as the critical value. Based on the estimation results of urban form changes under five scenarios, we forecasted the changes of UTEI in 2030 under five scenarios, as shown in Figure 6 (specific data in Table 4 is available in Table A1 in Appendix A).  The results show that under the KT scenario, the average UTEI of Jing-Jin-Ji cities will continue to increase. The average summer and winter UTEIs increase to 0.43 and 0.06 °C, respectively, by 2030. This means that if the government does not take measures to restrain the development of the UTEI, the problem of the urban thermal environment will become more and more serious in the future, especially in summer, when the heat pressure is higher. Under the two scenarios of optimizing the urban form, the average level of UTEI will decrease. In the SO scenario, the UTEI declined slightly, ranging from 0.05 to 0.13 °C, while in the EO scenario, the heat island effect was significantly alleviated, with decreases in the UTEI reaching 0.69 °C in summer and 0.56 °C in winter, respectively. In deteriorating situations, the UTEI enhancements were very significant, especially in summer. The average UTEI enhancements were as high as 0.55-0.76 °C in summer. In winter, the average level of the UTEI rose relatively little, but it also reached 0.24-0.29 °C. Our findings indicate that the UTEI will not be aggravated only under the SO and EO scenarios. Nevertheless, in order to improve the thermal environment of the city, it is better to develop the EO scenario.
Specifically, in the context of rapid urban development, compared with natural factors, it is more operable and meaningful to optimize urban form factors to alleviate urban thermal environment problems. At present, many studies have proved that urban form is an effective tool to alleviate the negative impact of urbanization on the thermal environment [54][55][56][57]. In spite of this, as the changes in urban structures are varied in cities and regions, the effectiveness of urban forms in mitigating heat-related risks also varies. The results show that under the KT scenario, the average UTEI of Jing-Jin-Ji cities will continue to increase. The average summer and winter UTEIs increase to 0.43 and 0.06 • C, respectively, by 2030. This means that if the government does not take measures to restrain the development of the UTEI, the problem of the urban thermal environment will become more and more serious in the future, especially in summer, when the heat pressure is higher. Under the two scenarios of optimizing the urban form, the average level of UTEI will decrease. In the SO scenario, the UTEI declined slightly, ranging from 0.05 to 0.13 • C, while in the EO scenario, the heat island effect was significantly alleviated, with decreases in the UTEI reaching 0.69 • C in summer and 0.56 • C in winter, respectively. In deteriorating situations, the UTEI enhancements were very significant, especially in summer. The average UTEI enhancements were as high as 0.55-0.76 • C in summer. In winter, the average level of the UTEI rose relatively little, but it also reached 0.24-0.29 • C. Our findings indicate that the UTEI will not be aggravated only under the SO and EO scenarios. Nevertheless, in order to improve the thermal environment of the city, it is better to develop the EO scenario.
Specifically, in the context of rapid urban development, compared with natural factors, it is more operable and meaningful to optimize urban form factors to alleviate urban thermal environment problems. At present, many studies have proved that urban form is an effective tool to alleviate the negative impact of urbanization on the thermal environment [54][55][56][57]. In spite of this, as the changes in urban structures are varied in cities and regions, the effectiveness of urban forms in mitigating heat-related risks also varies. Therefore, city managers need to combine regional characteristics, establish links between the thermal environment and various social and natural factors, clarify the direction and intensity of each factor's effect on the thermal environment, and exercise targeted macrocontrol on urban building expansion directions and surface vegetation coverage in the process of urban planning. Attention should be paid to changes in urban expansion trends, urban vegetation cover, and other natural-human factors that would influence the UTEI in the process of urban development. Changes that are conducive to alleviating the UTEI problem should be encouraged, and trends that promote the UTEI problem should be adjusted and improved.

Conclusions and Limitations
This study explores the influences of urban form and natural indicators on the urban thermal environment in the Jing-Jin-Ji urban agglomeration. Through the comparison of six linear and nonlinear models, we selected the best prediction model and calculated relative feature importance. We also simulated the possible urban variables and urban thermal environment changes under typical urbanization scenarios.
Results showed that nonlinear models (Adaboost, RF, and GBRT) performed better in predicting the UTEI than linear models (Lasso, BN, and BR models), among which the GBRT model had the highest prediction accuracy (the value of its ten-fold cross-validation reached 0.89). Feature relative importance results showed that natural factors and urban form factors had certain influence on the UTEI. Relative humidity and solar radiation have high contribution rates (about 10%). ELONG, POP, NLI, AREA, and SHAPE have high contributions to the generation of the UTEI among all the urban form indicators. We also confirmed the seasonal differences in the contribution of the urban environmental factors to the UTEI. Natural factors contribute more in summer, while urban form factors contribute more in winter.
The five scenarios' analyses showed that under the extremely deteriorated patterns, the average UTEI could be as high as 0.55-0.76 • C in summer and 0.24-0.29 • C in winter. On the other hand, in an extremely optimized situation, UTEI decreased by 0.69 • C in summer and 0.56 • C in winter. Overall, it indicated that cities need to develop under the two optimization scenarios in order to improve the urban thermal environment. City managers should pay more attention to changes in urban expansion trends, urban vegetation cover, and other natural-human factors that would influence the UTEI in the process of urban development. Changes that are conducive to alleviating the UTEI problem should be encouraged, and trends that promote the UTEI problem should be adjusted and improved. The research results also have important reference significance for future macro planning at the urban scale. Our research has successfully built high-precision machine learning models to clarify the influencing factors of the UTEI and set scenarios for future urban environment development. As there is a growing need for practical solutions to decreasing the urban thermal environment in the process of urbanization, it is a good attempt to apply machine learning models to solve the problem.
Several limitations of this paper can be summarized as follows: Firstly, this paper mainly studied the cities of the Jing-Jin-Ji region in China, most of which are small-and medium-sized cities. For bigger cities, the forecasting effect needs further confirmation. Secondly, the factors of the urban form we used are still limited. It is difficult to obtain more information on urban forms on an urban scale for a large number of cities. However, the multi-center structure of a city, the three-dimensional characteristics of a city, and so on, can provide more information about the city's social and economic activities and physical characteristics, which can be used as the direction of future research. Thirdly, this study considered the interaction between urban form factors and selected natural factors but did not elaborate on the mediating or moderating effects of the urban form and how much the mediating or moderating effects account for the change of the UTEI. We will further explore it in the future. Another limitation of this study was that some of our scenario settings may be difficult to achieve for all cities. In the future, the scenario setting method based on historical change statistics and expert experience could be combined with our model to explore the potential impact of urban form change on urban thermal environments.