Measuring the Urban Land Surface Temperature Variations in Zhengzhou City Using the Landsat-Like Data

Satellite-based remote sensing technologies are utilized extensively to investigate urban thermal environments under rapid urban expansion. Current MODIS data is, however, unable to adequately represent the spatially detailed information because of its relatively coarser spatial resolution, while Landsat data can’t explore temporally the refined analysis due to the low temporal resolution. In order to resolve this situation, we used MODIS and Landsat data to generate “Landsat-like” data by using the flexible spatiotemporal data fusion method (FSDAF), and then studied spatiotemporal variation of land surface temperature (LST) and its driving factors. The results showed that 1) The estimated “Landsat-like” data have high precision; 2) By comparing 2013 and 2016 datasets, LST increases ranging from 1.8°C to 4°C were measurable in areas where the impervious surface area (ISA) increased, while LST decreases ranging from -3.52°C to -0.70°C were detected in areas where ISA decreased; 3) LST has a strongly negative relationship with the Normalized Difference Vegetation Index (NDVI), and a strongly positive relationship with Normalized Difference Built Index (NDBI) in summer; and 4) LST is well correlated with Building density (BD), in a complex conic mode, and LST may increase by 0.460°C to 0.786°C when BD increases by 0.1. Our findings can provide information useful for mitigating undesirable thermal conditions and for long-term urban thermal environmental management.


Introduction
Worldwide urban expansions are blowing at various spatiotemporal scales under the rapid population and economy increase [1][2][3][4].The World Urbanization Prospects 2014 Revision report stated that 54%of the world's population resided in urban regions in 2014 and was projected to grow to 66% in 2050 [5].Inevitably, urban expansion caused land use and land cover change and brought about many environment or ecological system problems [6].Many natural vegetation and farmland were replaced by man-made impervious surface area leading to dramatic change of land surface characteristic, local urban climate and urban thermal environment [7][8][9].
With the development of satellite-based thermal remote sensing technologies, multi-parameter remotely sensed datasets were widely used to characterize land surface temperatures (LST) and to investigate thermal variations across urban environments [10][11][12][13].Because of the high-frequency of its multiple daily data-acquisition passes, the Moderate Resolution Imaging Spectroradiometer (MODIS) with 1000 meters resolution was most popular in studies that required a long series of scans or a high temporal resolution time scale [14][15][16][17].Landsat data with a higher spatial resolution of 60m or 100m in the thermal band is, however, more suitable for detailed investigation of LST, impervious surface cover (ISC), and land use changes, although the Landsat satellite passes less frequently [18][19][20][21][22].Both MODIS and Landsat data have their strengths and limitations.The low spatial resolution of MODIS scans can only be used to validate large-scale research.It is insufficient for producing detailed descriptions of LST variation, or for analyzing thermal driving forces [23].Similarly, the temporally sparse resolution of Landsat data cannot produce the short-interval datasets necessary to investigate LST dynamics [24].Recently, data fusion methods have been used to deal with these spatiotemporal resolution problems by combining different remote sensing datasets [25][26][27].For example, by utilizing the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM), Shen et al. [28] obtained "Landsat-like" LST datasets from 1988 to 2013, and analyzed the mode, evolutionary character, and thermal mechanics of the Wuhan city urban heat island.Huang et al. [29] also proposed a spatiotemporal fusion model based on a bilateral filter method (STBFM), and testified to its higher precision compared to the STARFM method.Wu P, et al. [30] put forward the spatiotemporal integrated temperature fusion model (STITFM), which can be used to integrate data from multiple sensors with flexibility.Besides the above three methods (STARFM, STBFM, and STITFM), a Flexible Spatiotemporal Data Fusion method (FSDAF) was proposed by Zhu et al. [31], which had a high calculation efficiency and minimum data requirements, although this final method is no longer used extensively for investigating LST.
Zhengzhou is a provincial city in the Central plains of China.It has the highest urban population density in mainland of China, is a key transportation junction, and a major economic center in the middle of China.It is currently undergoing a Great-leap-forward development [32].Moreover, the urban built-up area has increased rapidly, by 65.83% over the last decade.There is, therefore, an urgent requirement to quantitatively evaluate urban thermal environmental change, and explore the principal driving mechanisms.In this study, we used the FSDAF method to derive "Landsat-like" LST data, then used it to analyze the spatiotemporal variations in urban expansion and LST, and finally, we identified some land surface factors that effected the LST, including NDVI, NDBI, and building density.
Three Landsat Enhanced Thematic Mapper Plus (ETM+) images, free of cloud, taken on 28th June 2013, 22nd January 2014, and 20th June 2016, and one Landsat 8 Operational Land Imager/ Thermal Infrared Sensor(OLI-TRIS) image, with a thin cloud cover, taken on 22nd January 2017, were selected for this study.All four images (path 124/row 36) were downloaded from the USGS Earth Explorer (https://earthexplorer.usgs.gov)website.In addition, Daily MODIS LST data (MODIS 11A1) was obtained from the NASA Earthdata website (https://ladsweb.modaps.eosdis.nasa.gov).Finally, terrestrial meteorological data (including air pressure, temperature, and relative humidity, etc.), used for estimating land surface temperatures, was obtained from the Chinese National Meteorological Information Center (http://data.cma.cn).2013 vector map datasets for the study area, including the Zhengzhou administrative division map with properties of building footprints described by groundlevel boundaries, were obtained from the Geographical Information Monitoring Cloud Platform (http://www.dsac.cn).

Data Preprocessing and Image Classification
The preprocessing of ETM+/OLI-TRIS images included: Geometric correction, radiometric calibration, and atmospheric correction.All procedures were performed using the ENVI5.1 software package.First, the thermal infrared band was resampled to 30m.Then, geometric correction was carried out so that the root mean square error was guaranteed to be within 0.5 for the cell range.Next, after completing radiometric calibration by converting the image digital number (DN) values into the radiance values, atmospheric correction was performed, based on the meta information derived from images using the FLASSH model.Furthermore, for solving the problem of missing pixels in Landsat ETM + SLC-off images taken since 2003, the nearest-neighborhood interpolation method was used to predict missing pixel values, in order to reduce the errors in subsequent data processing.
After pre-processing, land use/land coverage was sorted into four types, based on a support vector machine (SVM) supervised classification method.This process incorporated an evaluation of selected training samples in which the types were designated according to a visual interpretation of their primary use or state.The categories defined were Impervious Surface Area (ISA, including building land, tarmac/cement road, and tarmac/cement ground), Vegetation, Water, and Bare Soil.For assessing the accuracy of the classification of the results, the User Accuracy and Overall Accuracy were evaluated based on the revalidation of samples.The overall accuracy of the classified images in the four study groups were 89.25%, 88.33%, 91.27%, and 94.69%, respectively, which met the accuracy requirements of this study.NDVI is a measure of the amount of vegetation at the surface [33,34],and NDBI has always been used to extract built-up areas [35,36] in the study area.They could be calculated with the equations (1) and ( 2) Where: NIR R is the reflectance in the near infrared band, which is corresponding to Landsat ETM+ band4 (0.775-0.900μm) and Landsat OLI-TIRS band5 (0.845-0.885μm Where:   S B T represents black-body radiance; L  is the spectral radiance at the sensor's aperture in W•m -2 •sr -1 •μm -1 ; S T stands for real land surface temperature (K); is the atmospheric transmittance, as well as L  and L  is the upward and downward atmospheric thermal radiance, respectively.The three parameters were estimated using the Atmospheric Correction Parameter Calculator provided by NASA Web site in conjunction with image meta information and meteorological data;  is the land surface spectral emissivity ,which is related to NDVI and was calculated based on NDVI threshold method in this study [37].
When NDVI < 0.05, the LULC is mainly bare soil and water, and the  of bare soil is 0.973, but for water,  = 0.995 [38]; When NDVI>0.7,it's regarded as full of vegetation coverage area, and  =0.99;When 0.05<NDVI<0.7,LULC is a mixed object type, and  is calculated with equation (4):.0.004 0.986 Where v P is the vegetation proportion calculated by equation ( 5): Where NDVI V and NDVI S represents the NDVI values of vegetation and soil, respectively.In this study, due to the lack of detailed surface spectral information, 0 DVI 5 N 0. V  and .NDVI 0 7 S  were defined to simplify the calculation process.

The Flexible Spatiotemporal Data Fusion Method
By combining the advantages of various spatiotemporal data fusion models, the Flexible Spatiotemporal Data Fusion Method (FSDAF) [31] ensures that minimum input requirement of data sources are met, while capturing the features of both gradual and abrupt changes, in order to produce fine-resolution images that accurately predict the LST distributions.Because of the strengths of this approach, the FSDAF method was chosen and applied to obtain high spatial-temporal resolution LST images in this study.Before using the FSDAF method to produce predicted fine LST images, the MODIS LST images (1000m resolution) were resampled to 30m resolution (the same as Landsat LST images), based on the nearest-neighbor algorithm method, and the input reference image was cropped to show the same regions.The implementation process was divided into 6 steps: 1) Classify Landsat LST images at time T1; 2) estimate the temporal changes taking place for each class of coarseresolution MODIS LST images, from T1 to T2; 3) predict the fine-resolution LST images at T2, based on the predicted temporal changes, and calculate pixel residuals of MODIS LST images; 4) use the thin plate spline interpolation function to predict the high-spatial-resolution LST images, based on the MODIS LST images at T2; 5) allocate the residuals to the predicted high-spatial-resolution LST images with the TPS interpolation function; 6) generate final fine resolution "Landsat-like" LST images at T2 based on the weights of pixels in the moving windows, which are assigned by nearestneighborhood information.The calculation process is as shown in ( 7)-( 13): represents the DN values of fine spatial resolution pixels;  , ,  ( ) x y b  is the residual which is allocated to the high spatial resolution pixel j from MODIS LST pixel i 。 The calculation of weight is carried out with reference to previous studies [26].And the TPS function that guides the residual distribution is conducted as shown in equations ( 9)-( 13): , , ) ( , , ) ( , ,  ( ) Where m is the number of sub-pixels in MODIS LST images; , , ( ) After several experiments in software using the ENVI IDL8.5 software package, parameter settings were defined and used in the study, includes as follows: 30 Pixels × 30 Pixels window size, 4 classes of LST image classification, and a similar pixel search threshold of 30.

Accuracy Assessment
As can be seen from Figure 2, the three types of LST images show consistent spatial distribution trends in the study area.The LST in the center of all the images is significantly higher in summer (top row) than the LST of the surrounding areas, but lower in winter (bottom row).From the overall data distribution maps, the rough image texture from the MODIS LST data products is clearly visible due to the low spatial resolution (1000m).Contours of partial vegetation coverage areas and water areas, as well as their corresponding low-temperature characteristics, are not clearly visible in the coarseresolution MODIS LST images.The 30m resolution "Landsat-like" LST images predicted by the FSDAF method, however, have the same level of fine texture details as the ETM + LST images obtained at the same time.In order to verify the accuracy level of predicted high-spatial resolution LST data, ETM+LST and MODIS LST images for 28th June 2013 (representing summer) and 22nd January 2014 (representing winter) were compared.A linear correlation comparison between predicted LST and observed LST was performed by the way of generating random points, as shown in predicted and ETM+LST data was, therefore, only 2.46°C and 1.73°C, respectively, and the RMSE was less than 2.3°C between predicted and MODIS LST data.There were high correlation coefficients (R²) between predicted and observed LST data, from 0.670 to 0.703.In addition, the regression coefficients between predicted "Landsat-like" LST and MODIS LST of 0.96 and 1.01 are closer to 1 than those between "Landsat-like LST" and ETM+ LST, of 0.68 and 0.81, which is due to the fact that the LST data is predicted based on the features of the simultaneous MODIS LST images.On the whole, the small errors and the high accuracy level of the predicted LST data fully meet the needs of this study.Figure 4 shows high-spatial-resolution LST images for each month of 2013, including ETM+/OLI-TRIS LST images, and four 30m resolution predicted "Landsat like" LST images (* represents the predicted LST images).The "Landsat-like" LST images are generated based on monthly MODIS LST data (the maximum synthesis of daytime MODIS LST products) using the FSDAF method.As shown in Figure 4, the high spatial-temporal resolution monthly LST datasets could provide sufficient data support for research concerned with the spatial-temporal distributions of LSTs and changes in the regional thermal environment over the long-term and with high-precision.

Spatiotemporal variations of LST
Figure 5 shows the changes in land use/land cover (LULC) ratios over the study time-series period.The total study area is 1028.35km 2 , and by using the classified results for four dates, it was found that the percentage for the total area classed as ISA changed from 37.91%, to 40.23%, to 47.12%, and finally to 48.88%, from June 28th, 2013 to January 22nd, 2017.The overall percentage increase due to ISA expansion was 10.97%, and the corresponding area increased by 112.84 km 2 .The Vegetation and Bare Soil coverage areas showed a significant seasonal variation due to the climatic characteristics of the study area: In summer, the area of Vegetation was larger than the area of Bare Soil; but in winter, the opposite was true, because the two classification groups are inversely related to some extent.The total area of Vegetation and Bare Soil, however, was also continuously decreasing, which is due to the transformation some of these areas into other LULC classification groups, most notably into ISA in some regions during urban expansion.The Water area is, in fact, largely unchanged over the study period.(2014/01/22 and 2017/01/21), the LST of the central urban area is significantly lower than that of the surrounding suburbs, which demonstrates a significant relationship; that the heat islands of the summer are cold islands in winter.This is consistent with the results of Zhou and Yao et al. [16,17].Furthermore, with regard to the distribution of LST ranges in the central urban area, it is clear that the LST for most regions in summer of 2016 was higher than in 2013, while for winter the differences were not significant.It can also be seen that the relatively higher LST regions in summer are concentrated in the central urban area, while the lower LST regions are mainly distributed in the southwest part of the study area, where urban rivers, lakes, and the northern Yellow River are located.On the other hand, in winter, the higher-LST characteristics are more prominent in the surrounding suburbs, while the lower-LST zones mostly appear in the central urban area, as well as in the vicinity of rivers and lakes.Comparing the spatial distribution changes in LST between 2013 and 2016 also reveals that while the higher-LST regions in summer increased significantly overall, they showed a trend of expansion towards the southwest, while the LST also rose slightly in the west and southeast of the study area.There were, however, few changes or trends in the spatial distribution of LSTs in the winter months.Figure 7 Shows statistics for the mean values of LST, NDVI, and NDBI, respectively, in the five different municipal zones.From 2013 to 2016, the mean LST in EQ district increased by 2.25°C in summer, and this was undoubtedly the most prominent LST change, compared to JS where the temperature decreased by 0.05°C, GC where it decreased by 0.44°C, ZY where it decreased by 0.68°C, and HJ where it decreased by 0.48°C.In winter, the mean LST for the five zones all showed a decreasing trend over time, especially in GC and EQ.At the same time, from Fig 7 (b) and (c), it can be seen that the mean values of NDVI and NDBI for each municipal zone also changed, and the distribution trends for these two indices and LST indicate an obvious seasonal correlation.In summer, the changes in mean LST and mean NDBI are greatest for GC, and smallest in EQ in 2013 and HJ in 2016.However, it is also clear that the maximum difference in NDBI is in zone EQ in summer, and the minimum is in zone JS in summer.In winter, however, there is little correlation between the spatial distribution of mean LST and mean NDBI.summer, mean LST values of ISA and Bare Soil basically remain at the maximum, while the mean LST in Vegetation areas is relatively lower.In winter, however, the mean LST remains highest in the Bare Soil areas, but is slightly lower in the ISA areas than in the Vegetation regions.
In order to further explore the difference of mean LST corresponding to different LULC types, randomly sampled points (n=6000) were used to make statistics regarding the LST distribution for each type of LULC, as are shown in Figure 8.For the different LULC types over the four dates, it can be seen that the LST medians for the sampled points have exactly the same distribution patterns as the mean LST values in Table 2, and the order is as follows: Bare Soil>ISA>Vegetation>Water (2013/06/28), ISA>Bare Soil>Vegetation>Water (2016/06/20), and Bare Soil>Vegetation>ISA>Water (2013/06/28 and 2017/01/21).Based on previous studies, one of the principle factors contributing to this phenomenon is that different LULC types have different specific heat capacities [9,39].Building materials of building land covered by ISA, such as industrial lands, residential areas and commercial regions, have lower specific heat capacities, so that the process of heat absorption and heat dissipation across these areas is more rapid.The relative LST, therefore, becomes much higher, whereas the larger specific heat capacities of the Vegetation and Water areas result in the smaller mean LST.  Figure 9(a) shows the distribution map of impervious surface area on June 28th, 2013 and June 20th, 2016, generated using the NDBI threshold method.As for the overall spatial distribution, areas where the ISA increased were concentrated around the borders of the various municipal zones, especially in the southwest of JS, the northwest of GC, the northeast of EQ, and the southeast of ZY.This shows a large-scale intensive growth trend and presents a pattern of clusters of expansion located around the main urban center.The areas where ISA decreased are mainly located in discrete regions of JS, GC, EQ, and ZY. Figure 9(b) shows the spatial distribution of LST difference between the two periods.This demonstrates that the areas in which the LST increased are mostly concentrated in the area where the JS, GC, EQ, and ZY zones meet, while the areas where the LST reduces occurs mainly in the northwest of JS, the southwest of GC, and eastern ZY.Comparison of Figure 9(a) and Figure 9(b) shows that there is a strong correlation between the ISA regional variation and the LST distribution differences.Most areas where ISA increased show an increasing trend in LST, especially in the main urban center where the borders of JS, GC, EQ, and ZY intersect, where the phenomenon is most significant, while areas where ISA decreased are mostly characterized by a reduction in LST.In the central and southern parts of the ZY area, the decrease of ISA is matched by the increase in LST.
Figure 10 quantitatively shows the mean LST difference for the three different ISA variation groups (more, same, less), across the five LULC types, for the summer of 2013 and 2016.Clearly, the mean LST difference for the area where the ISA increased in each zone remained the largest, ranging from 1.8°C to 4°C, meaning that the land surface in that zone experienced significant increase of temperature.The mean LST difference in the area where the ISA showed No Change was slightly smaller, ranging from -0.15°C to 1.98°C, and except for zone EQ, the LST was effectively unchanged in the other four zones where the LST difference was close to 0°C.The mean LST difference in the areas where the ISA decreased was between -3.52°C and -0.70°C, constituting the largest negative change, and in all regions the temperature change was below -2°C, apart from in zone EQ, showing a significant negative trend in the LST.From this data, it can be concluded that the LST significantly rises when other types of LULC are converted to ISA, and that it is relatively stable in areas where the ISA shows No Change, but in regions where the ISA decreased, the LST showed a significant negative trend.This is consistent with the results of previous studies [40,41], further indicating that ISA regional changes have a large impact on the spatial distribution of LST.

Driving factors of LST
NDVI is widely used to describe land surface vegetation coverage [33,42], while NDBI is a remote sensing index used for extracting data about the ISA and Built-up regions.These can quantitatively describe the regional distribution characteristics of these land use types, and the close correlation between NDVI-LST and NDBI-LST has also been discussed extensively [43].from this part of the study due to its special characteristics [44].Figure 11 shows the correlation between LST and NDVI and NDBI.Plotting NDVI against LST, for the summers of 2013 and 2016, shows that LST has a negative linear relationship with NDVI (the correlation coefficients R 2 are 0.425 and 0.549 respectively, p<0.001.);whereas there is a less pronounced positive linear relationship in winter.LST always relates strongly positively, however, with NDBI, although seasonal difference still exists, as the correlation in summer (R 2 =0.601 and 0.609) is stronger than in winter (R 2 =0.308 and 0.214).Furthermore, LST related more closely and robustly to NDBI than to NDVI, further illustrating that ISA changes during urban expansion have a more significant impact than vegetation coverage changes, on LST.These results are consistent with those obtained by Ma et al [45].The most prominent growth of ISA in areas undergoing rapid urbanization is due to a sharp increase in land areas used for building.Building density (BD) is the most effective parameter reflecting the expansion of building land, and it undoubtedly impacts LST [46].Based on the building vector map data for the study area in 2013, several square buffer zones with a spacing of 0.5 km were selected around the main urban center.LULC classification of images from 28th June 2013 (start date) to 22nd January 2017 (end date) allowed the ISA changes for each buffer zone to be statistically analyzed.When the side length of the buffer zone was 19 km, the area difference of ISA was stable to within 4 km 2 , which could be regarded as no change in building land coverage in this region.As shown in Figure 12, by using the method for calculating building density described in Guo et al. [47], a complex relationship between BD and LST was obtained.A nonlinear quadratic polynomial correlation model, like a conic mode, was developed for the corresponding four periods, and R 2 was 0.867, 0.782, 0.931, and 0.843, respectively (p < 0.001).From the distribution of sample points shown in Figure 12 (a), (b) and (d), BD-Mean LST shows a negative correlation first and then a positive correlation, while the BD values at the inflection point range from 0.2 to 0.26.Clearly, it can be concluded that the LST in areas with a BD of less than 0.26 is lower than that in medium and highdensity areas.BD correlates negatively with LST in this area, while BD and Mean LST have a strongly positive correlation when BD is greater than 0.26.Notably, on 20th June 2016, BD always remained positively related to LST, as shown in Figure 12(c).With the information shown in Table 2 and Figure 8, it is clear that this correlation pattern was due to differences in LST distribution among various LULC types on different dates.For example, the mean LST of Bare Soil (low building density) areas was highest during three dates (all except 20th June 2016), showing that in low density areas mean LST was most deeply influenced by Bare Soil, which was the main coverage type for these fishing grids.Bare Soil areas were gradually replaced by building land/ISA as the building density increased, leading to a reduction in overall LST in low density areas.In the medium and high-density areas, however, the larger proportion of building land areas played a decisive role in the changes of LST.however, the mean LST of the ISA area reached its maximum, so the BD-Mean LST always showed a positive trend, regardless of the building density of the area.The increase of building densities also resulted in an increase of population density, in the frequency of human activities, and in energy consumption, all of which contributed significantly to the changes in the thermal environment and the rise of LST.As Figure 12 shows, a positive linear correlation between BD and Mean LST was created in the medium-density and high-density regions where the BD was above 0.26, and it can be clearly seen that in those areas the LST increased by 0.691°C to 0.786°C for each 0.1 increase of BD in summer, which is more severe than that in winter when it was 0.460°C to 0.570°C.
Figure 12.The scatter plot of correlation between Mean LST and Building density for summer and winter in 2013 and 2016.A 120m*120m fishing grid (n=8971) was generated in the selected area, and the corresponding building area for each grid was evaluated and used as the building density (BD).
The mean LST was also calculated.Then the building density interval [0-1] was equally divided into copies, and the mean LST, at increments of 0.01 from the minimum and maximum BD, was calculated.Finally, the correlation between building density and LST was calculated, and the relationship was plotted.
In order to comprehensively and quantitatively analyze the driving force of each parameter on the LST, this study extracted the mean values of NDVI and NDBI from the fishing grid generated using BD statistics (n=8971).A multivariate nonlinear regression model was then constructed with three parameters and LST, as shown in Table 3.It appears that the correlation coefficients between multiple factors and LST, ranging from 0.791 to 0.931, are significantly greater than those between any single factor and LST (NDVI-LST, NDBI-LST and BD-LST).There are also still seasonal differences, in that the correlation in summer is stronger than that in winter.Furthermore, it seems that that the spatial distributions of LST in areas of urban expansion is closely related to human activities, which effect local climate characteristics and produce the changes of LULC.In conclusion, LST changes across the urban thermal environment results from the combined effects of multiple driving factors.

Discussion
The study of changes in the thermal environment due to urban expansion was based on Landsat LST images for many years, during which time researchers often used single day LST data to represent annual data, because of the satellites' long return time cycle, and frequent cloud contamination of images [42,48,49].This made it difficult to accurately describe the temporal variation characteristics of local thermal environments.The FSDAF method cited in this study can, however, be used to predict missing monthly high-spatial-resolution LST images, by using MODIS LST data products.Compared with the results of previous studies, the RMSE between predicted LST and the observed LST (Landsat LST images and MODIS LST products are considered to be the observed LST data) can be kept within 2.5°C [29,30] and the predicted LST has a high precision level.
In addition, the flexibility and performance of the FSDAF method is better than other spatio-temporal fusion models due to its higher efficiency and lower input requirement.This study analyzes the correlation patterns between LULC and LST under urban expansion from a holistic perspective.Based on the response of ISA transition regions to LST changes, the correlation characteristics of spatial distribution variations in ISA and LST were analyzed qualitatively and quantitatively.In addition, the impact of various factors on LST changes were quantitatively analyzed based on random sampling of points combined with calculations of building densities and NDVI and NDBI surface indices.A multivariate nonlinear regression model relating the three parameters and LST was also developed.Following the conclusions of most research projects in this subject area [21,43,50], it was seen that in summer months NDVI-LST shows a negative correlation, while NDBI has a quite strongly positive correlation with LST.However, it is different to the persistent positive correlation trend proposed by Guo et al. [47], as the current study indicated that there is a negative correlation first in the low-density areas where BD is less than 0.26, and then a positive correlation in the medium and high-density areas, between BD and Mean LST.In addition, the characteristics of summer and winter variations were studied and compared, and this showed that whether a single-factor or multi-factor analysis was carried out, there was always a seasonal difference, although the three driving factors correlated strongly in summer but related weakly in winter, to LST.
The 30m-resolution "Landsat-like" LST images produced by the FSDAF method are highly consistent with MODIS LST products, with respect to both their spatial distribution trends and data fitting, but analysis of the linear correlation between predicted LST and the ETM + LST data clearly shows that in areas with low LST, the predicted LST is higher than the observed LST.In areas with high LST the predicted value is lower than the observed LST, and these differences seem to be strongly related to the LULC distribution system for sampling points in the different LST ranges.Furthermore, the current study indicates that the BD-Mean LST relationship has a completely opposite correlation trend before and after the lowest point in the conic mode, which differs from the results obtained by Guo et al. [47], but also further indicates that the mechanism of LST change in building land areas is complicated and influenced by many factors, and factors such as energy consumption and human activities should be taken into consideration.
Changes in the proportion of LULC types in the fishing grid, before and after the lowest point on the trend line, and the spatial location of the inflection point, led to the development of a nonlinear correlation model relating BD and Mean LST, something that will be elucidated in a follow-up study.One notable issue is that the LST datasets used in this research were all collected during daytime, whereas in previous studies incorporated differences such as the spatial-temporal distribution variation of LST in Zhengzhou City between day and night [39,51,52].It is clear that using only daytime LST data is not a basis for comprehensively exploring changes in local thermal environments.The FSDAF method could, nevertheless, produce predictions of high-spatialresolution LST images taken at night, based on the corresponding characteristics taken from MODIS LST data products, and this methodology could be further developed to advance the study of spatialtemporal LST distribution in Zhengzhou city.

Conclusions
In this study, a complete dataset allowing the qualitative and quantitative spatial analysis of LST distribution was produced by combining Landsat LST images with 30m-resolution "Landsat-like" LST datasets predicted using the FSDAF method.The LST distribution was investigated on the basis of different LULC types, and the driving forces and driving patterns behind LST changes were evaluated by calculating surface indices NDVI/NDBI and building density.In this way, the spatialtemporal distributions of LST during the expansion of Zhengzhou city was analyzed.The results show that the high-spatial-resolution "Landsat-like" LST datasets predicted by the FSDAF method have a high accuracy level, and the overall trend is consistent with the Landsat LST images and MODIS LST products for the same time periods.This is useful for solving the problem of the lack of high-spatial-resolution LST datasets caused by infrequent satellite passes and cloud contamination of high-resolution images.With regards to the spatial distribution of LST in the study area, LST was higher in the main urban area than in the surrounding suburbs in the summer, showing an obvious characteristic of heat islands, in June 2013 and 2016.On the other hand, the LST appeared much higher in the surrounding suburbs due to the effects of the cold island in winter.Meanwhile, due to urban expansion, the areas with relatively high LST in summer became larger, especially in the EQ zone, which is closely related to the rapid expansion of ISA in those regions and the decrease of vegetation coverage.Compared to other LULC types, the seasonal variation of LST distribution was particularly prominent in ISA regions, which means that the LST is the highest in summer but lowest in winter in these areas.Notably, between the summer of 2013 and 2016, there was a strong correspondence between the regional changes in ISA coverage and the distribution of LST differences in the five municipal zones.LST increases ranging from 1.8°C to 4°C were measurable in the areas where ISA increased, while the LST presented a negative trend, ranging from -3.52°C to -0.70°C, in the areas where the ISA decreased.From the perspective of a quantitative analysis, there were also differences in the driving patterns and driving forces due to the different driving factors that produced LST changes.The NDVI-LST plot showed significant seasonal differences, for example in summer of 2013 and 2016 this revealed a negative linear correlation, with R 2 being 0.425 and 0.549 respectively, whereas there was almost no correlation in winter.On the other hand, the NDBI-LST plot always showed a strongly positive correlation trend in summer, with an R 2 value up to 0.609, but a less positive relationships in winter, where is reached a minimum of 0.214.The impact of NDBI changes on LST changes was also significantly stronger than the impact of NDVI changes.Unlike the situation with the former two indices mentioned, the relationship between building density and LST showed an extremely strong nonlinear quadratic polynomial correlation, with a negative correlation in low density regions but a positive correlation mode in medium and high-density regions.It can be seen, for example, that for an increase of 0.1 in BD when BD is greater than 0.26, the LST increases by 0.460°C to 0.786°C.Multivariate nonlinear regression analysis including the three parameters correlated more strongly with LST changes than the relationship between any single factor and LST, which further indicated that the spatial-temporal distribution of LST was the result of multiple factors.
In order to further investigate changes in the thermal environment during urban expansion, "Landsat-like" high-resolution LST data images will be predicted for nighttime in future research, and the FSDAF method will be used to compare differences between daytime and nighttime LST spatial-temporal distributions.The influences of driving factors will be also be comparatively analyzed under those conditions.Meanwhile, it is also possible to carry out more long-term exploratory research, on monthly, quarterly, and inter-annual time scales, to analyze the development of the thermal environment over different time scales.The climatic characteristics of the study area during the different time periods, such as air temperature and humidity, should be taken into consideration when investigating the impact of changes in the thermal environment.It is clear that this study, and future studies of this type, can provide practical information and guidance to assist local thermal environmental management, and to advise decision-making bodies about good land use practices to be applied during urban expansion.
of spatial resolution from T1 to T2; k  is the weight of the kth similar pixel; of class a of high spatial resolution LST data on band B from T1 to T2;


Landsat LST pixels and the fine-resolution LST pixels predicted based on the temporal changes; are high-spatial-resolution LST pixel values at T2 based on temporal changes and optimized TPS interpolation function parameters respectively, and b is the difference between two types of LST pixels; is the TPS function corresponding to band B.

Fig 3 .
Most of the random points on the predicted "Landsat-like" LST images for the two periods were evenly distributed beside the line Y=X.The root-mean-square error (RMSE) between the Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 11 September 2018 doi:10.20944/preprints201809.0192.v1

Figure 4 .
Figure 4.The distribution of monthly LSTs in 2013.* represents the predicted "Landsat-like" LST image produced using the FSDAF method.

Figure 6
Figure6shows the high-spatial-resolution LST images for the four study dates.From the seasonal overall trends, it can be determined that in summer (2013/06/28 and 2016/06/20) the relatively higher LST regions are mainly concentrated in the center of study area, while in winter

Figure 6 .
Figure 6.LST images for summer and winter in 2013 and 2016 (unit: °C).

Table 2
Shows statistics relating to the mean LST for each LULC type, corresponding to the four dates.With regard to the overall trends, the lowest LST is always located in a Water area, but the other three types of LULC have a more prominent seasonal difference in mean values for LST.In Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 11 September 2018 doi:10.20944/preprints201809.0192.v1

Figure 8 .
Figure 8. Box graphs of the relationship between LULC and LST (unit: °C).

Figure 9 .
Figure 9. Spatial correlation between ISA and LST regional changes in summer, between 2013 and 2016.(a) shows the ISA regional change, which is generated using the NDBI threshold method, and (b) shows the spatial distribution of LST difference.

Figure 10 .
Figure 10.Mean LST difference in the ISA change area over different zones in summer, between 2013 and 2016.
Figure 7 shows the relationships between LST and NDVI and NDBI for the different municipal zones.Based on these two surface indices, the correlations between NDVI/NDBI and LST are quantitatively analyzed in order to describe the temperature driving effects of LULC changes (mainly focusing on Vegetation and ISA land use types) on LST in urban expansion areas.Water type has been excluded Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 11 September 2018 doi:10.20944/preprints201809.0192.v1

Table 1 .
Accuracy assessment of Land Use/Land Cover (LULC) classification.

Table 2 .
Statics of the land surface temperature (°C) by different land use/land cover types.

Table 3 .
The multivariate linear regression models of the LST, BD, NDVI and NDBI