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

: Satellite-based remote sensing technologies are utilized extensively to investigate urban thermal environment under rapid urban expansion. Current Moderate Resolution Imaging Spectroradiometer (MODIS) data are, however, unable to adequately represent the spatially detailed information because of its relatively coarser spatial resolution, while Landsat data cannot explore the temporally continued analysis due to the lower temporal resolution. Combining MODIS and Landsat data, “Landsat-like” data were generated by using the Flexible Spatiotemporal Data Fusion method (FSDAF) to measure land surface temperature (LST) variations, and Landsat-like data including Normalized Di ﬀ erence Vegetation Index (NDVI) and Normalized Di ﬀ erence Built Index (NDBI) were generated to analyze LST dynamic driving forces. Results show that (1) the estimated “Landsat-like” data are capable of measuring the LST variations; (2) with the urban expansion from 2013 to 2016, LST increases ranging from 1.80 ◦ C to 3.92 ◦ C were detected 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 signiﬁcant negative correlation with the NDVI and a strong positive correlation with NDBI in summer. Our ﬁndings can provide information useful for mitigating undesirable thermal conditions and for long-term urban thermal environmental management.


Introduction
Worldwide urban expansions are developing at various spatiotemporal scales under the rapid population and economy growth [1][2][3][4]. The World Urbanization Prospects 2014 Revision report stated that 54% of the world's population lived in urban regions in 2014 and the number would be 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]. Especially, the coupling of rapid urban expansion and climate warming has increased the heat stress in many mega cities [7]. A great deal of natural vegetation and farmland was replaced by man-made impervious surface area, leading to dramatic change of land surface temperature (LST) and local climate [8][9][10].
LST is an indicator showing the biosphere-land-ocean-atmosphere interactions and surface energy balance, and it can reflect the thermal response to urban canopy, building height, surface coverage, anthropogenic heat, and energy consumption [11][12][13].
Traditionally, near-surface air temperatures as a proxy of in-situ LST can be accurately measured by using in-situ thermometers but are limited in station-based observation [14,15]. With the development of earth observation technologies, a range of remotely sensed datasets have been used to characterize LST and to investigate thermal variations across urban environments [16][17][18]. Satellite thermal infrared (TIR) LST, theoretically based on the radiative transfer equation, has been widely studied focusing on single-channel, split-window, and multi-channel algorithms for solving emissivity and atmospheric effects [19][20][21]. For example, Landsat and Moderate Resolution Imaging Spectroradiometer (MODIS) data are intensively used on different spatiotemporal scale LST studies [22][23][24][25][26][27][28][29][30][31] (Table1). Landsat data with a higher spatial resolution of 60 m or 100 m in the thermal band are more suitable for detailed investigation of LST, however, Landsat data has lower temporal resolution (16 days) [32][33][34][35][36]. MODIS data, with higher temporal resolution (twice a day), are most popular in continuous studies satisfying a high temporal resolution time scale [37,38]. The low spatial resolution of MODIS scans can only be used to coarse-scale research. It is insufficient for producing detailed descriptions of LST variation or for analyzing thermal driving forces [39]. Similarly, the temporally sparse resolution of Landsat data cannot produce the short-interval datasets necessary to investigate LST dynamics [40].
Both MODIS and Landsat data have their strengths and limitations. To obtain simultaneously the higher spatiotemporal resolution data, data fusion methods have been used to deal with these spatiotemporal resolution problems by generating the fusion LST data of Landsat and MODIS, satisfying the need of long-term and fine-scale regional thermal environment research [41][42][43][44]. For example, by utilizing the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM), Shen et al. [45] obtained "Landsat-like" LST datasets from 1988 to 2013 and analyzed the thermal mechanism of the Wuhan city urban heat island. Huang et al. [46] proposed a spatiotemporal fusion model based on a bilateral filter method (STBFM) and testified its higher precision compared with the STARFM method. Weng et al. [47] developed a new fusion algorithm to predict daily LST. Wu et al. [48] put forward the spatiotemporal integrated temperature fusion model (STITFM) to integrate data from multiple sensors with flexibility. Besides the above three methods (STARFM, STBFM, and STITFM), a Flexible Spatiotemporal Data Fusion method (FSDAF) [49], which has a strong advantage in generating high spatiotemporal resolution data due to high calculation efficiency and minimum data requirements, was adopted by Zhang et al. [50] to generate monthly time series LST data. Although there is different research focusing on the new data fusion algorithms to generate higher spatial and temporal resolution land surface factors, there are still few documents to indicate simultaneously the LST dynamics and the driving forces by combining the Landsat-like data Normalized Difference Vegetation Index (NDVI) and Normalized Difference Vegetation Index (NDBI). Zhengzhou is a provincial city in the Central plains of China. It has the highest urban population density in mainland of China and is a major economic center in the middle of China. It is currently undergoing a great-leap-forward development [51]. 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. The purpose of this study includes: (i) measuring urban land surface temperature dynamics under fast urban expansion by using Landsat-like LST, and (ii) analyzing LST driving forces combining Landsat-like NDVI and NDBI and other factors.

Materials
The study area shown in Figure 1 is the Zhengzhou, located from 113°26′ to 113°51′ E, and 34°35′ to 34°57′ N. It comprises five municipal zones: Erqi (EQ), Zhongyuan (ZY), Huiji (HJ), Jinshui (JS), and Guancheng (GC). The built-up area of this conurbation expanded significantly, from 328.01 km 2 in 2011 to 443.04 km 2 in 2016 according to the dataset of National Bureau of Statistics of the People's Republic of China.
Data used in this paper mainly include two parts (Table S1)

Data Preprocessing and Image Classification
The preprocessing of ETM+/OLI-TRIS images (Level 1) includes radiometric calibration and atmospheric correction by using the ENVI5.1 software package. First, the thermal infrared band was corrected with "Thermal Atmospheric Correction" tool in ENVI5.1, and then it was resampled to 30 m by using the bilinear method. Next, for optical bands, atmospheric correction was performed by using the FLASSH function. Finally, the digital number (DN) of the thermal band was validated into the radiance, while the DN of the optical band was validated into reflectance. Furthermore, for solving the problem of missing pixels in Landsat ETM + SLC-off images taken since 2003, the nearestneighborhood interpolation method was used to predict missing pixel values. In addition, a 120 m*120 m fishnet (n=8971) was generated to calculate the mean impervious surface area (ISA) ratio. The interval of ISA ratio was equally divided into 100 segments from zero to one, and corresponding mean LST in each interval was calculated.
The support vector machine (SVM) supervised classification method was used for image classification. As a machine learning method based on the statistical learning theory, it has a strong precision level [52][53][54]. The land use includes four types, impervious surface area (ISA), vegetation, water, and bare soil. Furthermore, we used the random sampling method to select more than 100 samples to assess the accuracy of the classification results.

Calculation of NDVI, NDBI, and LST
NDVI (Normalized Difference Vegetation Index) and NDBI (Normalized Difference Built-Up Index), as two typical land surface factors, can be calculated as follows: where 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).
Landsat LST calculation is based on the radiative transfer equation: where ( ) represents the surface radiance; is the spectral radiance at the sensor's aperture in W·m -2 ·sr -1 ·μm -1 ; stands for real land surface temperature (K); is the atmospheric transmittance; and ↑ and ↓ are the upward and downward atmospheric thermal radiance, respectively. The three parameters were estimated using image meta information and meteorological data as well as the atmospheric correction parameter calculator [55,56] provided by NASA Website; is the land surface spectral emissivity, which is critical for accurate calculation of LST.
In this study, was calculated based on NDVI threshold method to differentiate types of land use/land cover (LULC). When NDVI < 0, the LULC was mainly water, but when 0≤NDVI<0.15, LULC was dominated by urban impervious areas and bare soil [57]. When NDVI≥0.727, it was regarded as full of the vegetation coverage area [57]. However, when 0.15≤NDVI<0.727, LULC became a mixed object type [58]. Therefore, was calculated with Equation (4): Finally, (℃) corresponding to the real LST was obtained with the Planck formula: The values of and are associated with the effective wavelength of the thermal band of each sensor [59]. For the ETM+ image, is 666.09 W/(m 2 *μm*sr) and is 1282.71 K, but for the OLI-TIRS image, is 774.89 W/(m 2 *μm*sr) and is 1321.08 K.

The Flexible Spatiotemporal Data Fusion Method
By combining the advantages of various spatiotemporal data fusion models, the Flexible Spatiotemporal Data Fusion Method (FSDAF) [49] ensures the minimum input requirement of data sources while capturing the features of both gradual and abrupt changes. The MODIS LST data (1000 m resolution) were resampled to 30 m resolution based on the nearest-neighbor algorithm method. Additionally, all input images needed to be cropped into the same area with the predetermined number of rows and columns before performing the spatiotemporal fusion method. The input data included two parts: 1) a pair of Landsat LST and MODIS LST data at T1, which were used to estimate the spatial differences between fine-resolution pixels and coarse-resolution pixels; 2) one MODIS LST at T2 , which was used as the spatial feature reference and applied to calculate the time differences between T1 and T2 for predicting fine-resolution LST at T2. The implementation process was divided into six steps: 1) classify Landsat LST at time T1; 2) estimate the temporal changes taking place for each class of coarse-resolution MODIS LST from T1 to T2; 3) predict the fine-resolution LST at T2 based on predicted temporal changes and calculate pixel residuals of MODIS LST; 4) use the thin plate spline (TPS) interpolation function to predict the high-spatial-resolution LST based on the MODIS LST at T2; 5) allocate the residuals to predicted high-spatial-resolution LST with the TPS interpolation function; and 6) generate final fine resolution "Landsat-like" LST at T2 based on the weights of pixels in the moving windows, which are assigned by nearest-neighborhood information. The calculation process is as shown in (6)- (12): where , , represents the DN values of fine spatial resolution pixels; , , is Landsat pixel values;∆ ( , , ) is the change of spatial resolution from T1 to T2; is the weight; ∆ ( , ) represents the change of class a of high spatial resolution data on band b (represents LST, NDVI, and NDBI data) from T1 to T2; , , is the residual, which is allocated to the high spatial resolution pixel from MODIS LST pixel 。 The calculation of weight was carried out with reference to previous studies [41], and the TPS function that guides the residual distribution was conducted as shown in Equations (8)-(12): , , , = ( , , ) , , = , , − , , , , = ( , ) where m is the number of sub-pixels in MODIS LST ; ( , , ) represents the residual values between Landsat LST pixels and the fine-resolution LST pixels predicted based on the temporal changes; ∆ , , represents the pixel value changes of band b in MODIS LST from T1 to T2; , , and , , are high-spatial-resolution LST pixel values at T2 based on temporal changes and optimized TPS interpolation function parameters, respectively, and ( , , ) is the difference between two types of LST pixels; , , represents the weights of guiding residual allocation; , , is the normalized weight of , , ; is homogeneous coefficient; and ( , ) is the TPS function corresponding to band b.
To better understand the procedures of how to generate the Landsat-like LST by using FSDAF, a workflow is shown in Figure 2. By using the ENVI IDL8.5 software package, the FSDAF model could be operated after many parameter settings were defined, including 30 pixels × 30 pixels window size, four classes of LST classification, and a similar pixel search threshold of 30. Table 2 shows the input and the output of the FSDAF method for predicting fine-resolution LST data. A pair of inputs (one MODIS LST and one Landsat LST) in T1 and another single input [MODIS LST data (as a reference of spatial feature)] in T2 were used in the FSDAF method to generate the output "Landsat-like" LST data in T2. To assure better fusion effects, the pair of input data should generally be cloud-free and close to the acquisition date of the single input data.
[ ]  As for generating monthly-series "Landsat-like" LST data, there were two steps: 1) a pair of daily LST data (Landsat LST and MODIS LST <MOD11A1>) and 8-day MODIS LST (MOD11A2) were selected as inputs in the FSDAF model, and the output-8-day "Landsat-like" LST-could be predicted; and 2) four predicted 8-day "Landsat-like" LSTs in a month were integrated with maximum synthetic method to obtain monthly "Landsat-like" LST. If there was no Landsat image available in one month, the Landsat image of the previous or the next month was chosen, as well as the corresponding MOD11A1 in that date. If there were two Landsat images available in one month, four predicted 8-day "Landsat-like" LST were calculated according to the first and the second half month, respectively. In the study, for further discussing the FSDAF model significance on a serial analysis while data missing or cloud contaminating, the integrated monthly LST dynamics in a whole year were performed.

Land-like LST Accuracy Assessment
The accuracy of Landsat-like LST was assessed by comparing it with MODIS LST and Landsat LST in three periods by using random points (n=3000). On 28 June 2013 and 22 January 2014, Landsatlike LST data were generated for verification with Landsat LST on the same date ( Figure 3). Because of the cloud coverage on 22 January 2017, by using MODIS LST on 21 January 2017, the predicated Landsat-like LST was generated and assessed ( Figure 4). The root-mean-square errors (RMSEs) between the Landsat-like LST and the ETM+LST data were 2.46 °C and 1.73 °C, and the RMSEs were slightly smaller between Landsat-like LST and MODIS LST data, corresponding to 2.30 °C and 1.30 °C, respectively. In addition, the regression coefficients between Landsat-like LST and MODIS LST (0.96 and 1.01) were closer to one than those between Landsat-like LST and ETM+ LST (0.68 and 0.81) because the generation of Landsat-like LST was based on the MODIS LST. Figure 3 shows the predicted Landsat-like LST was highly similar to MODIS LST with a 0.976 regression coefficient (R 2 =0.753). Spatially, all 30 m resolution Landsat-like LST had approximately the same "cold-hot" texture details distribution with Landsat LST and MODIS LST.  As shown in Table 3, for various LULC types, the RMSEs between the Landsat-like LST and the Landsat LST data were also mostly higher than those between the Landsat-like LST and the MODIS LST data except for the RMSE of water on 28 June 2013. The RMSEs in summer were higher than those in winter. The applicability of the FSDAF method in generating high spatial-temporalresolution LST is greatly affected by different LULC effects [48,60].

LST Variations under Urban Expansion
Urban expansions and land use/land cover (LULC) changes over the period are illustrated in Figure 5 and  (Table  S2).  In summer, the relatively higher LST (more than 50 °C) regions were mainly concentrated in the center of the study area, where the highest LST on 20 June 2016 was up to 61.89 °C, while the air temperature was 30.4 °C that day. In winter, the LST of the central urban area was significantly lower than that of the surrounding suburbs (Figure 7), and the lowest LST was -3.92 °C on 22 January 2014, while the air temperature was 2.5 °C. It is clear that the LST of most regions in summer of 2016 was higher than that of 2013, while for winter, the differences were not significant. Although the heterogeneity of the land surface causes the varieties of LST spatially, LSTs are mainly dominated by the distribution of LULC. In summer, the lower LST regions were mainly distributed in the southwest and the north covered by water bodies and grass, while in winter, the higher LST characteristics were more prominent in the surrounding suburbs, while the lower LST zones mostly appeared in the central urban area. Comparing the spatial distribution changes in LST between 2013 and 2016, we found the higher LST regions in summer increased significantly, showing a trend of expansion towards the southwest, while the LST also increased slightly in the west and the southeast of the study area.  (Figure 8b). The mean values of NDVI and NDBI also changed (Figure 8c, d).
In the EQ region with a big park and green land, the mean NDVI was the highest, while the NDBI was the smallest in summer. The NDBI in winter seasons was larger than that in summer because of the withering vegetation in grassland and few crops covering some farmlands. By comparing the changes in two periods, both the maximum increase of NDBI and the maximum decrease of NDVI in summer were identified in EQ region. The spatial distributions of LST and three indices (including ISA ratio, mean NDVI, mean NDBI) indicate an obvious consistent or opposite trend change for each municipal zone in summer.   Urban expansion includes not only the total area increase but also internal urban land increase. Figure 10 shows spatially the ISA changes and the corresponding LST differences between two summers of 2013 and 2016. The ISA increase regions existed in various municipal zones, especially in the southwest of JS, the northwest of GC, the northeast of EQ, and the southeast of ZY (Figure 10a). The LST increase regions mostly occurred in the area where the LST decreased ( Figure 10b).
As shown in Figure 11(a), the mean LST differences of two periods were calculated for parts where ISA changed or not in each zone. For ISA increase type, LST increased with the range 1.9 to 3.5 °C, showing a "warming effect" of the building land. LST decreased with the range -3.8 to -0.5 °C where ISA decreased. The mean LST difference presented a small increase and ranged from -0.15 °C to 1.98 °C where ISA had no change. By comparing two term changes in EQ zone, the largest LST increased 3.5 °C for the parts where ISA increased. LST increased where other types of LULC were converted to ISA; on the contrary, LST decreased where ISA regions were converted to non-ISA regions. As shown in Figure 11(b) and (c), including ISA increased region and ISA decreased region, a strong nonlinear quadratic polynomial correlation between ISA ratio change and LST change existed (R 2 s were 0.698 and 0.652, respectively, p < 0.001). As for the ISA increased areas, LST showed 0.27 °C growth as ISA ratio increased by 0.1. However, for ISA decreased regions, LST showed two different trends. LST increased by 0.25 °C as ISA ratio increased by 0.1 in the range from -1 to -0.3, while LST decreased by 0.22 °C as ISA ratio increased by 0.1 in the range from -0.3 to 0.

Driving Factors on LST
NDVI and NDBI are two important factors reflecting the LULC change. Figure 12 shows correlations of NDVI-LST and NDBI-LST to reveal the driving forces (mainly focusing on vegetation and ISA) on LST. LST had a negative linear relationship with NDVI (R 2 = 0.425 and R 2 = 0.549, p<0.001), whereas there was a less obvious linear relationship in winter. NDBI had a positive relationship with LST that was relatively stronger in summer (R 2 = 0.601 and 0.609) and weaker in winter (R 2 = 0.308 and 0.214). The NDBI had a closer relationship with LST than NDVI, further illustrating that the building land increase had a more significant impact on LST than vegetation coverage during urban expansion.

Integrated Monthly LST Dynamics Based on Landsat-like Data
For generating monthly "Landsat-like" LST using MODIS 11A2 products and Landsat LST, four 8-day "Landsat-like" LST data in each month were generated to integrate the monthly LST with maximum value composite method. For example, the first two predicated LST (Figure 13a1  As shown in Figure 14, with the available Landsat images of each month (Table 4), the integrated monthly "Landsat-like" LST of 2013 were generated. The continuous high resolution monthly LST data can provide sufficient information to illustrate the serial LST dynamics in the regional thermal environment over the long-term and with high precision. The monthly LST values of different percentiles in 2013 were calculated, and the distribution trend is shown as Figure 15. From January to December, distinct seasonal variations showed the highest LST (36.45 °C) in summer and the lowest LST (9.29 °C) in winter.  The relationship of NDVI-LST and NDBI-LST is illustrated in Figure S1. The correlation coefficient R 2 s of NDVI-LST and NDBI-LST in summer were above 0.4 (June, July, and August of 2013), while R 2 was rather weak in winter ( Figure 16). From the monthly series of correlation analysis, the positive correlation of NDBI-LST was more stable and stronger than the negative relationship of NDVI-LST. At the same time, the surface parameter NDVI had the largest driving force on LST. The seasonal difference in the correlation of NDVI-LST was significant, especially among summer and the other seasons, while the correlation of NDBI-LST only presented the lowest in winter and remained stably strong during the year of 2013. As shown in Figure 17, air temperature exhibited a similar monthly trend to LST (R 2 =0.878). In August, the highest monthly air temperature was 30 °C, while the highest monthly LST was up to an average 37.44 °C for the whole city. However, the highest pixel LST was more than 50 °C, which was almost 13 °C higher than the air temperature in August. A moderate negative correlation (R 2 =0.696) was detected between particulate matter PM2.5 and LST. Generally, PM2.5 can reduce the amount of surface solar irradiances reaching the ground surface and cause a cooling effect on the surface. In contrast, there were no strong results showing effects of wind speed and rain on LST variations.

Discussion
Spatiotemporal fusion models can make up the limitations of remotely sensed data missing or cloud contamination in the related thermal environment researches, and they provide a new method for previous researchers who had to use single daily LST data representing continuous data, because of the satellites' long return time cycle [50,61]. Especially, the flexibility and the performance of the FSDAF method has higher efficiency and lower input requirement. Compared with the results of previous studies [48], the RMSE between predicted LST and the observed LST in this study could be kept within 2.50 °C, showing a high precision level. However, the result of Landsat-like LST is dependent on the date of input Landsat image acquisition and is especially influenced by using the replacement data because of data missing or cloud contamination. Figure 18 shows the closer the input Landsat image acquisition date to the predicted date is, the higher the accuracy of the result is. By selecting Landsat LST of 28 June 2013 as the reference, "Landsat-like" LST predicted 4 June 2013, which was the closest to the output predicted date showing the smallest differences from actual LST [ Figure 18 Figure 18(c)]. Also, the selection date of the paired input data also affects the predicted Landsat LST. The input coarse MODIS LSTs of various dates have differences in the satellite viewing angles and sun geometry, among others, which also cause some deviations between actual LST and predicted LST. As to the serial LST dynamics analysis, 8-day Landsat-like LSTs were predicated and then were combined to generate the integrated monthly Landsat-like LST by using the maximum value composite method. However, due to the cloud coverage influence, uncertainties and new methods should be solved while using the FSDAF method for large scale LST analysis [62,63]. The urban area of Zhengzhou city expanded by 28.94% from 2013 to 2016. This is similar to previous studies that revealed a continuous and rapid expansion in the process of Zhengzhou urbanization for many years [64,65]. The urban expansion affected LST with great variations. In summer, the center of the urban area had a relatively higher LST than the surrounding suburbs, while it presented a completely opposite LST distribution in winter. This phenomenon is widespread in most major cities in north and central-north China. In Zhou's research, spatial distributions of surface Urban Heat Island Intensity (SUHII) of Zhengzhou ranged from 2 °C to 3 °C in summer and from -0.5 °C to 0 °C in winter [32]. Yao et al. also proposed that the SUHII of Zhengzhou ranged from 2 °C to 4 °C in summer and from -2 °C to 0 °C in winter [40]. However, more attention should be paid to explain uncertainties about the urban cooling effect, especially by combining the results of the numerical simulation and the remote sensing retrieval method. As for the driving forces on LST, similarly to other research [66], in summer months, NDVI-LST showed a negative correlation, while NDBI had a quite strongly positive correlation with LST. However, due to the seasonal fluctuations of vegetation, the relationship of NDVI-LST changed dramatically and showed weaker correlations in winter, as reported in other research [67,68]. Different LULC types had different specific heat capacities, and vegetation was transferred to built-up land, which increased LST [69]. This is consistent with the results of previous studies [70], indicating that ISA changes have a large impact on the spatial distribution of LST.
In addition, except for ISA, NDVI, and NDBI, a range of climatic factors including air temperature, wind speed, and rain can also affect the LST variations [23]. Research has disclosed climatic multiple effects on LST and pointed out the potential strategies for heat stress alleviation and urban planning [71,72]. For example, the strong correlation between air temperature and LST could be used to perform estimation of spatially distributed near surface air temperature with constructing models [73,74]. However, according to Figure 17a, we can see the biggest difference between LST and air temperature was 15 °C in March. The reason for this may have been the fast vegetation growth inducing the LST increase. In the whole year, LST fluctuated in some months while there was only one peak in the air temperature curve. That shows the complex thermal variations mechanics of land surface in respect to radiance and emission. Although we pointed out the negative impacts of PM2.5 on LST, there are still great uncertainties. For example, Cao et al. pointed out urban haze/aerosol pollution was also a great contributor that intensifies the high surface temperature agglomeration effect on the center of the urban area at night [75]. The wind effects may differ according to different location and topography. The change of the urban boundary layer wind field reduces the LST of the urban center, especially in coastal cities, and good ventilation can better improve the effects of the urban thermal environment [76]. With long-term and multi-site data, it could be disclosed that there is a negative correlation presented between LST and precipitation [73]. In addition to these significant factors that have a great impact on the LST changes, there are still some possible contributors, e.g., incoming surface radiation, duration and intensity of sunlight, elevation, etc. [77]. LST variations are affected by the contribution from multi factors, and comprehensive analysis is helpful for fully understanding the relevant driving forces and driving modes [78].

Conclusions
In this study, to overcome the shortage of coarser spatial resolution of MODIS data and lower temporal resolution of Landsat data, "Landsat-like" datasets were generated using the FSDAF method, including LST, NDVI, and NDBI. Under the urban expansion, LST variations and related driving factors were analyzed.
The results show that high spatiotemporal-resolution "Landsat-like" LST has a high accuracy level, and the overall trend is consistent with Landsat LST and MODIS LST for the same date. From July 2013 to July 2016, the urban area expanded from 389.85 km 2 to 484.56 km 2 . LST increases ranging from 1.80 °C to 3.92 °C were detected in areas where the impervious surface area increased, while LST decreases ranging from -3.52 °C to -0.70 °C were detected in areas where ISA decreased. The relationship between NDVI and LST showed significant seasonal differences; in summer, it revealed a negative linear relationship (R 2 = 0.425 and R 2 = 0.549, p<0.001). NDBI presented a strongly positive correlation with LST in summer with R 2 more than 0.6. Monthly Landsat-like LSTs were generated from four predicated 8-day "Landsat-like" LST in each month to indicate comprehensive dynamics. Also, monthly Landsat-like NDVI and NDBI, air temperature, PM2.5, and other factors were analyzed to disclose the driving forces. In future research, temporal and spatial fusion efficiency and precision can be improved by using multi-source data or deep-learning algorithms. Also, urbanization effects on thermal environment should be analyzed by coupling different dimensional data including land surface factors, energy consumption, and urban morphology. This study can provide practical information and guidance to assist local thermal environmental management and to advise decision-makers about good land use practices to be applied during urban expansion.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1: Remote sensing datasets, Table S1: The monthly series about the correlations of LST with NDVI and NDBI in 2013, Table S2: Accuracy assessment of Land Use/Land Cover (LULC) classification.