Analysis of Long Time Series of Summer Surface Urban Heat Island under the Missing-Filled Satellite Data Scenario

Surface urban heat islands (SUHIs) are mostly an urban ecological issue. There is a growing demand for the quantification of the SUHI effect, and for its optimization to mitigate the increasing possible hazards caused by SUHI. Satellite-derived land surface temperature (LST) is an important indicator for quantifying SUHIs with frequent coverage. Current LST data with high spatiotemporal resolution is still lacking due to no single satellite sensor that can resolve the trade-off between spatial and temporal resolutions and this greatly limits its applications. To address this issue, we propose a multiscale geographically weighted regression (MGWR) coupling the comprehensive, flexible, spatiotemporal data fusion (CFSDAF) method to generate a high-spatiotemporal-resolution LST dataset. We then analyzed the SUHI intensity (SUHII) in Chengdu City, a typical cloudy and rainy city in China, from 2002 to 2022. Finally, we selected thirteen potential driving factors of SUHIs and analyzed the relation between these thirteen influential drivers and SUHIIs. Results show that: (1) an MGWR outperforms classic methods for downscaling LST, namely geographically weighted regression (GWR) and thermal image sharpening (TsHARP); (2) compared to classic spatiotemporal fusion methods, our method produces more accurate predicted LST images (R2, RMSE, AAD values were in the range of 0.8103 to 0.9476, 1.0601 to 1.4974, 0.8455 to 1.3380); (3) the average summer daytime SUHII increased form 2.08 °C (suburban area as 50% of the urban area) and 2.32 °C (suburban area as 100% of the urban area) in 2002 to 4.93 °C and 5.07 °C, respectively, in 2022 over Chengdu City; and (4) the anthropogenic activity drivers have a higher relative influence on SUHII than other drivers. Therefore, anthropogenic activity driving factors should be considered with CO2 emissions and land use changes for urban planning to mitigate the SUHI effect.


Introduction
The surface urban heat island (SUHI)-a phenomenon in which the land surface temperature (LST) tends to be higher in urban center zones than surrounding suburban surfaces, is usually measured using satellite thermal remote sensing data [1,2].The SUHI effect is one of the greatest concerns for its adverse impacts on air and water quality, energy consumption, and urban dwellers' health during heat wave events [3].The SUHI phenomenon has been observed worldwide, especially in developing countries such as China [4].In the summer of 2022, China faced its most severe heatwave in over six decades [5].The province of Sichuan, situated in western China experienced recordbreaking temperatures.This exacerbated the challenges faced by urban residents, including power outages, which were compounded by a widespread drought that severely affected both food and factory production across the province [6,7].Consequently, accurately quantifying the intensity of the SUHI effect (SUHII) and better understanding the driving factors has become imperative.These measures not only aid in assessing potential heatrelated risks but also contribute to future city management strategies, guiding governmental decision-making [8].
LSTs retrieved from satellite thermal infrared (TIR) bands are key indicators for quantifying SUHIs [9].Satellite remote sensing has supplied effective and unique methods for acquiring LST data with frequent coverage.However, adverse atmospheric conditions coupled with long revisit cycles have largely limited satellite-derived LST applications in urban thermal environments [10].Especially in the most rainy and cloudy cities in China, the large missing rate of satellite data is a common and serious problem.For a single satellite sensor, a tradeoff occurs between spatial and temporal resolution.Thus, there is a significant requirement to develop a method capable of integrating remotely sensed data from diverse sensors to produce fine spatiotemporal resolution LSTs for a better understanding of SUHI dynamics [11].In the last ten years, numerous methods for spatiotemporal fusion have been suggested to achieve high-resolution LSTs by combining the high spatial resolution and the high temporal frequency of diverse remote sensing data sources [12].The available spatiotemporal data fusion approaches that have been experimentally tested on LST products are mainly classified into four categories: weighted function-based [13], unmixing-based [14], learning-based [15], and hybrid methods [16].The spatiotemporal fusion technique for LST data offers the opportunity to further understand the SUHI phenomenon [17,18].In the present review, the spatial and temporal adaptive reflectance fusion model (STARFM) [14], the enhanced STARFM (ESTARFM) method [19], the bilateral filter [9], the spatiotemporal adaptive data fusion algorithm (SADFAT) [20], and the spatiotemporal integrated temperature fusion model (STITFM) [21] in the weight function-based category; the pixel-based multi-spatial resolution adaptive fusion modeling framework (pMSRAFM) [22] in the unmixing-based category; the sparse-representation-based spatiotemporal reflectance fusion model (SPS-FTM) [23] in the learning-based category; and the flexible spatiotemporal data fusion (FSDAF) method [24] in the hybrid category.Although great progress has been made, most of these spatiotemporal fusion methods are unable to accurately capture the spatial details of LSTs, predict abrupt events, and preserve the spatial continuity of LSTs within urban areas simultaneously [25].The above methods have their own advantages and limitations.For example, the weight function-based category has the most methods developed, owing to high computation efficiency, simple parameters, and strong robustness.However, this category has weaknesses in the strong temporal variability in LSTs that makes them more sensitive to model parameters [16], particularly, the size of the moving window, and it is not feasible in heterogeneous urban areas.FSDAF can simultaneously predict dense time series LST data owing to its ability to predict abrupt changes and gradual change events, but also easily causes spatial discontinuity in urban LST data.Shi et al. [26] proposed a comprehensive flexible spatiotemporal data fusion (CFSDAF) method based on FSDAF and generated a high-spatiotemporal-resolution LST image, which can preserve the spatial continuity and spatial details of LST in urban areas.At present, the application of LST in urban environment studies requires more heat-related information at the urban district level with high spatial resolution [27].However, high-resolution LSTs may also be derived from the Landsat series TIR channels (i.e., Landsat 5, 8, and 9) at about 100 m but remain far from meeting the needs for improving SUHI monitoring accuracy.
In this study, Chengdu, a typical cloudy and rainy city in southwestern China has been selected as the study area.Chengdu has very few satellite images that are available for use, due to the annual average number of 340 days that experience cloudy and rainy weather.The main purposes of this study were to (1) propose a multiscale geographically weighted regression (MGWR) coupling CFSDAF method to generate a 30 m spatial resolution and 8 days temporal resolution summer LST dataset from 2002 to 2022, and produce higher accuracy in urban areas compared to other traditional spatiotemporal fusion methods; Sensors 2023, 23, 9206 3 of 25 and (2) perform quantitative analyses to investigate the influence of multiple natureanthropogenic driving factors on the summer SUHII in Chengdu City.

Study Area
Chengdu, the capital city of Sichuan province and the sixth largest city in China (103 •   C on seven days [28].The long-term extreme heat phenomenon easily leads to air pollution and public health problems.This study focuses on an area of 60 km × 54 km in Chengdu, which covers the core urban area, a smaller suburban area with 50% of the urban area, and a larger suburban area with 100% of the urban area (Figure 1).One challenging problem for monitoring the summer SUHI effect in Chengdu is the large rate of missing satellite LST data owing to many cloudy and rainy days throughout the whole year.
olution and 8 days temporal resolution summer LST dataset from 2002 to 20 produce higher accuracy in urban areas compared to other traditional spatiot fusion methods; and (2) perform quantitative analyses to investigate the influ multiple nature-anthropogenic driving factors on the summer SUHII in Chengdu

Study Area
Chengdu, the capital city of Sichuan province and the sixth largest city i (103°57′ E-104°20′ E, 31°15′ N-31°41′ N), has experienced rapid urbanization in century.In 2022, Chengdu's gross domestic product (GDP) reached 2080 billion lars.The population of the city has exceeded 21.2 million people, 15.4 million living in urban areas.Rapid urbanization induces significant SUHI effects, espe summer, which could lead to extreme heat events.From 5 to 24 August 2022, C experienced a record months-long heatwave, which exceeded 40 °C on seven d The long-term extreme heat phenomenon easily leads to air pollution and publi problems.This study focuses on an area of 60 km × 54 km in Chengdu, which co core urban area, a smaller suburban area with 50% of the urban area, and a lar urban area with 100% of the urban area (Figure 1).One challenging problem fo toring the summer SUHI effect in Chengdu is the large rate of missing satellite L owing to many cloudy and rainy days throughout the whole year.

Data Description and Preprocessing
The proposed MGWR-CFSDAF method mainly needs, at least, a pair of h low-spatial-resolution LST data on the prior date and one set of low-spatial-re LST data on the predicted date.In this study, due to limitations of cloudy an weather, we can only select seventeen Landsat LST and HJ-1B LST data from 2002 as the high spatiotemporal LST data through blending with low-spatial-re MODIS LST data for predicting high-spatiotemporal-resolution LST data.As sh Figure 2, there are ten high-spatial-resolution LST images without clouds and se images have a cloud cover ranging from 0% to 10%.

Data Description and Preprocessing
The proposed MGWR-CFSDAF method mainly needs, at least, a pair of high and low-spatial-resolution LST data on the prior date and one set of low-spatial-resolution LST data on the predicted date.In this study, due to limitations of cloudy and rainy weather, we can only select seventeen Landsat LST and HJ-1B LST data from 2002 to 2022 as the high spatiotemporal LST data through blending with low-spatial-resolution MODIS LST data for predicting high-spatiotemporal-resolution LST data.As shown in Figure 2, there are ten high-spatial-resolution LST images without clouds and seven LST images have a cloud cover ranging from 0% to 10%.(1) Landsat 5/8/9 LST data.The Landsat thermal infrared (TIR) channels have a minimum 16-day revisit cycle and spatial resolution of about 100 m, as Landsat 5 collects TIR channel data at 120 m spatial resolution while Landsat 8/9 has two TIR bands at 100 m spatial resolution.Landsat images are available from the U.S. Geological Survey (http://earthexplorer.usgs.gov/,accessed on 20 January 2023).Radiometric calibration and atmospheric correction were performed.We retrieved Landsat LST data from Landsat 5 TIR band 6 and Landsat 8/9 TIRS band 10 using a generalized single-channel method.For details of the generalized single-channel method, please refer to Jimenez-Munoz and Sobrino [29].The details of the Landsat, HJ-1B, and MODIS used in this study are summarized in Table 1.(1) Landsat 5/8/9 LST data.The Landsat thermal infrared (TIR) channels have a minimum 16-day revisit cycle and spatial resolution of about 100 m, as Landsat 5 collects TIR channel data at 120 m spatial resolution while Landsat 8/9 has two TIR bands at 100 m spatial resolution.Landsat images are available from the U.S. Geological Survey (http://earthexplorer.usgs.gov/,accessed on 20 January 2023).Radiometric calibration and atmospheric correction were performed.We retrieved Landsat LST data from Landsat 5 TIR band 6 and Landsat 8/9 TIRS band 10 using a generalized single-channel method.For details of the generalized single-channel method, please refer to Jimenez-Munoz and Sobrino [29].The details of the Landsat, HJ-1B, and MODIS used in this study are summarized in Table 1.
(2) HJ-1B LST data.The HJ-1B images used in this study are level-2 output products and were obtained from the China Center for Resources Satellite Data and Application (https://data.cresda.cn/#/mapSearch/,accessed on 25 January 2023).The spatial resolution of the HJ-1B TIR band is 300 m with a 4-day revisit cycle [30].The HJ-1B data were geometrically corrected using calibrated Landsat 8 images within the study area.The error was controlled within 0.5 pixels to meet the geometry correction requirements.Then, the HJ-1B data were radiometrically calibrated using calibration coefficients [31] to convert the digital number (DN) values of the raw HJ-1B images into satellite radiance images [32].Finally, the ENVI-FLAASH module was used for atmospheric correction on each HJ-1B CCD image after radiometric calibration.In this research, LST data were retrieved from the thermal band IRS4 of the HJ-1B imagery using the single-channel algorithm.For a more comprehensive description of the single-channel algorithm, please refer to the work of Duan et al. [33].(3) MODIS LST data.One kilometer spatial resolution Daily Terra MODIS daytime LST (MOD11A1) data and 8-day Terra MODIS daytime LST (MOD11A2) data were obtained using the generalized split-window algorithm from the Geospatial Data Cloud (http: //www.gscloud.cn/,accessed on 28 January 2023).Numerous research findings indicate that the root mean square error (RMSE) of the MODIS LST data are within 2.0 K and exhibit high accuracy in major global cities [34].MODIS LST data were re-projected to the same coordinate system as Landsat and HJ-1B using MODIS Reprojection Tools (MRT).Finally, we utilized the quality control band within MOD11A1 and MOD11A2 to identify pixels affected by cloud contamination, with the purpose of excluding them from subsequent analysis.
(4) In situ LST.In situ hourly LST data were collected from Chengdu Meteorological Office's 7 weather stations distributed across Chengdu City in summer (April to September) from 2002 to 2022.In situ LSTs were collected based on SI-111 infrared radiometers with an accuracy of ±0.2 K.
(5) Potential driving factors of SUHII.To explore the potential driving factors of SUHII in Chengdu City, thirteen driving factors were selected in this study and divided into four types: the satellite precipitation product (PRE), wind speed (WS), relative humidity (RH), and white sky albedo (WSA) form the climate types; nighttime light index (NLI), perpendicular impervious surface index (PISI), and PM 2.5 form the anthropogenic activity types; population counts (POP), population density (PD), and gross domestic products (GDP) form the population shift types; and enhanced vegetation index (EVI), bare-soil index (BI), and normalized difference water index (NDWI) form the natural land surfaces types.Table 2 shows the potential driving factors based on available data selected in this study.The population data we collected are two products of the WorldPop dataset [56].They describe the residential population where they actually live.

Generating High-Spatiotemporal-Resolution LST for SUHI Monitoring
In this study, in order to generate the 30 m spatial resolution and 8-day temporal resolution summer LST dataset from 2002 to 2022, a multiscale geographically weighted regression MGWR coupling CFSDAF method was proposed.The implementation consists of testing the proposed method part and monitoring the summer daytime SUHI part (Figure 3).If the performance of the first part is better, we can conduct the next part.

Generating High-Spatiotemporal-Resolution LST for SUHI Monitoring
In this study, in order to generate the 30 m spatial resolution and 8-day temporal resolution summer LST dataset from 2002 to 2022, a multiscale geographically weighted regression MGWR coupling CFSDAF method was proposed.The implementation consists of testing the proposed method part and monitoring the summer daytime SUHI part (Figure 3).If the performance of the first part is better, we can conduct the next part.In the first part (testing the proposed method), both downscaled high-spatial-resolution LST data using the MGWR model and MOD11A1 data captured on 20 April 2013, 16 April 2015, 2 April 2018, and 21 April 2022 were used as the base time (t 1 ) for the proposed method.While the other MOD11A1 data captured on 21 May 2013, 10 July 2015, 5 June 2018, and 7 May 2022 were used as the prediction time (t2) input base data for predicting the high-spatial-resolution LST at t2.In the following part, we call LST directly from MODIS, HJ-1B, Landsat as "observed LST", and LST derived from the proposed method of other spatiotemporal fusion methods as "predicted LST".The coefficient of determination (R 2 ), the RMSE and the absolute average difference (AAD), were computed between the predicted LST images and the observed LST images to validate the accuracy of the predicted LST data.
In the second part (monitoring summer daytime SUHI from 2002 to 2022), we selected thirteen pairs of downscaled LST with 30 m spatial resolution and MOD11A1 data with 1000 m spatial resolution during the same period as the input base LST data at t1 (Table 1).Afterwards, MOD11A2 at t2 was used to fuse the predicted LST datasets with a temporal resolution of 8 days and 30 m spatial resolution.Finally, the predicted summer daytime LSTs were averaged and used to monitor the summer SUHI from 2002 to 2022 in Chengdu City.

Downscaling LST Using MGWR
Classical geographically weighted regression (GWR) as a downscaling method cannot capture the spatial non-stationary relationship between LSTs and environmental variables [70].Unlike GWR, MGWR can build a nonstationary relationship between LSTs and multiple environmental variables [71].MGWR is introduced to analyze the scale differences in normalized vegetation index (NDVI), digital elevation model (DEM), slope, and aspect on the spatial pattern of LSTs.The MGWR model was employed to downscale LST changes from low-resolution LST data to high-resolution LST data.The mathematical expression of the MGWR is as follows [72]: where Y i is the predicted values of dependent variable (LST in our case) and i = 1, 2,3, . .., n; and β bw0 is a the intercept at optimal bandwidth.X ij represents the jth predictor variable with spatially varying regression coefficient (β bwj ) over spatial locations (µ i , v i ).The error term in the model is represented by ε i .
(1) LST retrieval from Landsat 5, Landsat 8, Landsat 9 and HJ-1B were aggregated to the spatial resolution of 1000 m.NDVI, DEM, slope, and aspect were extracted at a 30 m spatial resolution based on the Landsat imagery, HJ-1B imagery, and other auxiliary data, whereas these environmental variables were aggregated to the spatial resolution 100 m, 120 m, 300 m, and 1000 m, respectively.
(2) MGWR was used to establish a nonstationary relationship between LST 1000 , and NDVI 1000 as well as DEM 1000 , Slope 1000 , and Aspect 1000 , which can be expressed as LST 1000 = f NDVI 1000 , DEM 1000 , Slope 1000 , Aspect 1000 (2) where LST 1000 is the LST estimated by the scale conversion function at 1000 m spatial resolution scale; NDVI 1000 , DEM 1000 , Slope 1000 , and Aspect 1000 are environmental variables at 1000 m spatial resolution; f (.) is the MGWR converts the auxiliary variables to simulate LST.
(3) Influenced by soil moisture and other physical parameters, it is difficult to fully reflect the spatial heterogeneity of LST, which is manifested as LST residual information at low-spatial-resolution scales: where ∆LST s is the LST transformation residual at 1000 m spatial resolution; LST s is the LST data estimated by the MGWR; and LST s is the LST at a 1000 m spatial resolution.

Implementation of the Proposed Method
Figure 3 presents a detailed producer of the proposed method.In this study, the CFSDAF method was used to fuse high-spatiotemporal-resolution LST images in the study area, in order to monitor summer SUHII, by combining the MOD11A1, MOD11A2, and downscaled MGWR LST images.The calculation process can be expressed as follows:

Implementation of the Proposed Method
Figure 3 presents a detailed producer of the proposed method.In this study, the CFSDAF method was used to fuse high-spatiotemporal-resolution LST images in the study area, in order to monitor summer SUHII, by combining the MOD11A1, MOD11A2, and downscaled MGWR LST images.The calculation process can be expressed as follows: where LST t 2 x ij , y ij represents the predicted high-resolution LST image at prediction data t 2 ; LST t 1 x ij , y ij represents the high-resolution LST data at base time t 1 ; k is the kth similar pixel; n is the number of similar pixels for central pixel in a single window; and ∆LST x ij , y ij is the prediction of the total change of the target pixel x ij , y ij between t 1 and t 2 .
In this study, CFSDAF mainly includes the following six steps: (1) adjust the differences between high-spatial-resolution LST and low-spatial-resolution LST and high-spatialresolution LST; (2) classify high-spatial-resolution LST after extracting the endmembers; (3) obtain the temporal increments by the linear equation of spatial unmixing process; (4) obtain the spatial increments by inverse distance weighting (IDW) interpolation; (5) integrate the spatial and temporal increments; and (6) obtain the LST prediction by the information of neighborhood.For more detailed steps of the CFSDAF model kindly refer to previous studies [26].

SUHII Analysis
In this study, SUHII is the LST difference between urban (LST urban ) and suburban areas (LST suburban ) using the fused high spatiotemporal resolution summer LST dataset from 2002 to 2022 over Chengdu City.The formula is as follows [8]:

Boosted Regression Tree Model
The application of the machine learning statistical model, boosted regression tree (BRT), is employed to investigate the influences of thirteen potential driving factors on SUHI.The BRT model exhibits strong learning capabilities and adaptability to diverse data formats, even when handling complex data, without necessitating the consideration of interactions or correlations among independent variables.Furthermore, it offers significant advantages in exploring interactions between complex factors and making forecasts.The BRT model has found successful applications in a wide range of fields, including urban expansion, ecological modeling, and environmental science [73][74][75].
In this paper, the gbm package with the statistical programming software R (version 3.3.2) was used to analyze the contribution of potential driving factors to SUHI.The dependent variables are the SUHII, and the independent variables are the thirteen driving factors.The BRT model is a supervised learning method; three parameters were specified after testing.In this study, the learning rate, bagging fraction, and decision tree complexity were 0.01, 0.5 and 5, respectively.In this study, this model extracted 50% of the data points for training, with 50% of the data used to fit thirteen driving factors and the first regression tree is SUHI.

Land Cover Classification
In order to monitor the SUHII in Chengdu, defining the urban and suburban areas was the first step.In this study, support vector machine (SVM), as one of the machine learning algorithms, was used for image classification [76].Cloud-free Landsat 5 images were acquired on 25  Using Google Earth, we randomly chose 1600 sample points, 400 for each type, for the accuracy assessment (Figure 1).We also compared the performance of three machine learning classifiers-SVM, artificial neural networks (ANN), and maximum likelihood classification (MLC).Table 3 illustrates the overall accuracy and kappa coefficients.The result showed that SVM performed better than ANN and MLC.
Due to Chengdu being the sixth-largest city in China, its administrative boundaries encompass not only urban areas but also extensive suburban regions, which do not align with the requirements of the SUHI study.Therefore, in this study, urban and suburban areas were separated according to four land cover classification types from landcover maps.An urban area is defined as a high-intensity and densely occupied areas near a built-up area.After the urban area is determined, a suburban area is defined as the buffer zone that includes a smaller suburban area (50% of the urban area) and a larger suburban area (100% of the urban area) around the urban area (Figure 5).

Testing the Proposed Method
MGWR is an extension of the generalized linear regression, with NDVI, DEM, slope, and aspect as the nominated environmental variable set that is highly LST-related.Firstly, high-spatial-resolution LST at the t1, such as Landsat 9 LST observed on 21 April 2022 (Figure 6) with 100 m resolution, was aggregated to 1000 m.The LST downscaling results from MGWR from 1000 m to 100 m is shown in Figure 6, where four kinds of subareas, notably subarea (Figure 6a) in vegetation, subarea (Figure 6b) in built-up area,

Testing the Proposed Method
MGWR is an extension of the generalized linear regression, with NDVI, DEM, slope, and aspect as the nominated environmental variable set that is highly LST-related.Firstly, high-spatial-resolution LST at the t1, such as Landsat 9 LST observed on 21 April 2022 (Figure 6) with 100 m resolution, was aggregated to 1000 m.The LST downscaling results from MGWR from 1000 m to 100 m is shown in Figure 6, where four kinds of subareas, notably subarea (Figure 6a) in vegetation, subarea (Figure 6b) in built-up area, subarea (Figure 6c) in bare soil area, and subarea (Figure 6d) in waterbody area, were used to show the LST downscaling performance using MGWR.Visually, MGWR can extract more spatial texture information from land surface temperature data, effectively revealing temperature distribution variations within similar land cover types.To assess the accuracy of the method of MGWR in downscaling LST, GWR and thermal image sharpening (TsHARP), which possess the advantage of LST downscaling, were also used in this study.RMSE and the mean error (ME) as the evaluation metrics were used to quantitatively evaluate the performance of the three downscaling methods, which is shown in Table 4.We can see that the MGWR possesses lower RMSE and ME compared to GWR and TsHARP; it shows MGWR produces higher LST downscaling accuracy than other methods from 2002 to 2022 in the study area.Therefore, the downscaled LST images and MOD11A1 at t1 could be used as the LST base data of CFSDAF for predicting the high-spatial-resolution LST data at t2.The next step is testing the performance of the proposed method in the first part (Figure 3). Figure 7a,f,k,p was downscaled using MGWR LST at t1 on 20 April 2013, 16 April 2015, 2 April 2018, and 21 April 2022, respectively.Figure 7b,g,i,q was MOD11A1 as the similar time at t1.To assess the accuracy of the method of MGWR in downscaling LST, GWR and thermal image sharpening (TsHARP), which possess the advantage of LST downscaling, were also used in this study.RMSE and the mean error (ME) as the evaluation metrics were used to quantitatively evaluate the performance of the three downscaling methods, which is shown in Table 4.We can see that the MGWR possesses lower RMSE and ME compared to GWR and TsHARP; it shows MGWR produces higher LST downscaling accuracy than other methods from 2002 to 2022 in the study area.Therefore, the downscaled LST images and MOD11A1 at t 1 could be used as the LST base data of CFSDAF for predicting the high-spatial-resolution LST data at t 2 .The next step is testing the performance of the proposed method in the first part (Figure 3). Figure 7a,f,k,p was downscaled using MGWR LST at t 1 on 20 April 2013, 16 April 2015, 2 April 2018, and 21 April 2022, respectively.Figure 7b,g,i,q was MOD11A1 as the similar time at t 1 .Figure 7c,h,m,r was the MOD11A1 data at t 2 on 21 May 2013, 10 July 2015, 5 June 2018, and 7 May 2022, respectively, for predicting the LST data at 100 m and 300 m spatial resolution on the same predicted date at t 2 (Figure 7d,i,n,s).The observed LST data at t 2 (Figure 7e,j,o,t) can be used to evaluate the predicted LST results as the similar time at t 2 .Figure 8 shows scatter plots of correlations between observed LST and predicted LST on 21 May 2013, 10 July 2015, 5 June 2018, and 7 May 2022.We can see some of the scatters deviate a lot from the fitted line owing to the predicted LST images being affected by weather conditions like thin clouds and fog.However, the accuracy assessment shows that the R 2 ranges from 0.8103 to 0.9476, RMSE from 1.0601 to 1.4974, and AAD from 0.8455 to 1.3380 on the same dates, which proves the proposed method has a better performance for predicting the high-spatiotemporal-resolution LST data.Figure 8 shows scatter plots of correlations between observed LST and predicted LST on 21 May 2013, 10 July 2015, 5 June 2018, and 7 May 2022.We can see some of the scatters deviate a lot from the fitted line owing to the predicted LST images being affected by weather conditions like thin clouds and fog.However, the accuracy assessment shows that the R 2 ranges from 0.8103 to 0.9476, RMSE from 1.0601 to 1.4974, and AAD from 0.8455 to 1.3380 on the same dates, which proves the proposed method has a better performance for predicting the high-spatiotemporal-resolution LST data.In addition, CFSDAF and FSDAF were used to evaluate the performance of the predicted LST results using the proposed method (Figure 9).The R 2 , RMSE, and AAD between the predicted LST and the observed LST on 21 May 2013, 10 July 2015, 5 June 2018, and 7 May 2022 show that the proposed method can be used to improve the fusion accuracy of high-spatial-resolution LST.The proposed method, with higher accuracy than CFSDAF and FSDAF on different dates, which shows the performance of the spatiotemporal fusion model, is more sensitive to the spatial resolution scale.In addition, CFSDAF and FSDAF were used to evaluate the performance of the predicted LST results using the proposed method (Figure 9).The R 2 , RMSE, and AAD between the predicted LST and the observed LST on 21 May 2013, 10 July 2015, 5 June 2018, and 7 May 2022 show that the proposed method can be used to improve the fusion accuracy of high-spatial-resolution LST.The proposed method, with higher accuracy than CFSDAF and FSDAF on different dates, which shows the performance of the spatiotemporal fusion model, is more sensitive to the spatial resolution scale.

Monitoring Summer SUHII from 2002 to 2022
In this study, the proposed method has a better performance for predicting the high-spatiotemporal-resolution LST data in the first part (testing the proposed method), which suggests that we can conduct the next step to predict high-spatiotemporal-resolution LST data using the proposed method with 30 m spatial resolution and 8-day temporal resolution in summer for monitoring summer SUHII in Chengdu City from 2002 to 2022 (Figure 10).Since there is no real satellite-derived LST data at 30 m spatial resolution, in situ LST data were used to validate the predicted LST results.In addition, the CFSDAF model, the FSDAF model, and the observed MOD11A2 at t2 were also used to evaluate the performance of the proposed method.Table 5 shows the proposed method can produce higher accuracy predictions of high-spatiotemporal-resolution LST data than other spatiotemporal fusion methods.

Monitoring Summer SUHII from 2002 to 2022
In this study, the proposed method has a better performance for predicting the highspatiotemporal-resolution LST data in the first part (testing the proposed method), which suggests that we can conduct the next step to predict high-spatiotemporal-resolution LST data using the proposed method with 30 m spatial resolution and 8-day temporal resolution in summer for monitoring summer SUHII in Chengdu City from 2002 to 2022 (Figure 10).Since there is no real satellite-derived LST data at 30 m spatial resolution, in situ LST data were used to validate the predicted LST results.In addition, the CFSDAF model, the FSDAF model, and the observed MOD11A2 at t 2 were also used to evaluate the performance of the proposed method.Table 5 shows the proposed method can produce higher accuracy predictions of high-spatiotemporal-resolution LST data than other spatiotemporal fusion methods.Both averaged summer LSTs from 2002 to 2022 were computed (Figure 11).As shown in Figure 11, due to the rapid urbanization, the spatial distribution of the LSTs showed an irregular distribution, and the high LST areas changed from urban to suburban areas.The high LST areas were mainly concentrated in the urban high-density blocks, where buildings and population were highly concentrated.The low LST areas were mainly concentrated in the mountainous areas of the suburban, such as the Long- Both averaged summer LSTs from 2002 to 2022 were computed (Figure 11).As shown in Figure 11, due to the rapid urbanization, the spatial distribution of the LSTs showed an irregular distribution, and the high LST areas changed from urban to suburban areas.The high LST areas were mainly concentrated in the urban high-density blocks, where buildings and population were highly concentrated.The low LST areas were mainly concentrated in the mountainous areas of the suburban, such as the Longquan Mountain Range in the southeast of the study area.The spatial resolution characteristics of the LSTs from 2002 to 2022 were similar, such as the high LST areas covered by the built-up areas.The low LST areas were mainly distributed in the vegetation-covered areas, such as farmland and mountain areas outside the built-up areas.We can see that the increase in the SUHI effect was roughly in the "southeast-northwest" direction as the urban built-up area expanded.The expansion of the SUHI effect corresponds to the urban spatial growth pattern.quan Mountain Range in the southeast of the study area.The spatial resolution characteristics of the LSTs from 2002 to 2022 were similar, such as the high LST areas covered by the built-up areas.The low LST areas were mainly distributed in the vegetation-covered areas, such as farmland and mountain areas outside the built-up areas.We can see that the increase in the SUHI effect was roughly in the "southeast-northwest" direction as the urban built-up area expanded.The expansion of the SUHI effect corresponds to the urban spatial growth pattern.Figure 12 shows the summer SUHII from 2002 to 2022.The significantly increasing trends of the summer SUHII in Chengdu used the averaged 30 m predicted LSTs.The highest SUHII for summer occurred in 2022 (5.07 °C from a larger suburban area and 4.93 °C from a smaller suburban area).Summer SUHII increased from 2.32 °C in 2002 to 5.07 °C for a larger suburban area and increased by 2.85 °C in the same period for a smaller suburban area.This indicates a large SUHII in the summer from 2002 to 2022 over the Chengdu City.It not only increases the risk of heatwave extreme events but also presents a big challenge for scientists to mitigate serious SUHI effects.

Relationship between SUHI and Potential Driving Factors
As described in Table 2, thirteen driving factors were selected to evaluate t fluence on summer SUHI.The driving factors were divided into climate driving anthropogenic activity driving factors, population shift driving factors, and natur surfaces driving factors.Figure 13 presents the results of a BRT analysis for Ch City.The relative influence of each factor is scaled as a percentage [77].Overall, erage the most important factors are PISI, EVI, and NLI, with about 26.9%, 17.4 12.5%, respectively.The other influences range from high to low are POP, PD WSA, BI, PM2.5, NDWI, RH, PRE, and WS, with 9.7%, 9.5%, 9.0%, 3.1%, 3.0%, 2.9% 1.8%, 1.4%, and 1.1% on average, respectively.For the population shift driving (Figure 14), the natural land surfaces driving factors (Figure 15), the climate drivi tors (Figure 16), and the anthropogenic activity driving factors (Figure 17), PD WSA, and PISI are the most influential factors with the influence of 37.6%, 50.1% and 59.2%, respectively.

Relationship between SUHI and Potential Driving Factors
As described in Table 2, thirteen driving factors were selected to evaluate their influence on summer SUHI.The driving factors were divided into climate driving factors, anthropogenic activity driving factors, population shift driving factors, and natural land surfaces driving factors.Figure 13 presents the results of a BRT analysis for Chengdu City.The relative influence of each factor is scaled as a percentage [77].Overall, on average the most important factors are PISI, EVI, and NLI, with about 26.9%, 17.4%, and 12.5%, respectively.The other influences range from high to low are POP, PD, GDP, WSA, BI, PM 2.5 , NDWI, RH, PRE, and WS, with 9.7%, 9.5%, 9.0%, 3.1%, 3.0%, 2.9%, 1.9%, 1.8%, 1.4%, and 1.1% on average, respectively.For the population shift driving factors (Figure 14), the natural land surfaces driving factors (Figure 15), the climate driving factors (Figure 16), and the anthropogenic activity driving factors (Figure 17), PD, EVI, WSA, and PISI are the most influential factors with the influence of 37.6%, 50.1%, 50.8, and 59.2%, respectively.

Relationship between SUHI and Potential Driving Factors
As described in Table 2, thirteen driving factors were selected to evaluate their influence on summer SUHI.The driving factors were divided into climate driving factors, anthropogenic activity driving factors, population shift driving factors, and natural land surfaces driving factors.Figure 13 presents the results of a BRT analysis for Chengdu City.The relative influence of each factor is scaled as a percentage [77].Overall, on average the most important factors are PISI, EVI, and NLI, with about 26.9%, 17.4%, and 12.5%, respectively.The other influences range from high to low are POP, PD, GDP, WSA, BI, PM2.5, NDWI, RH, PRE, and WS, with 9.7%, 9.5%, 9.0%, 3.1%, 3.0%, 2.9%, 1.9%, 1.8%, 1.4%, and 1.1% on average, respectively.For the population shift driving factors (Figure 14), the natural land surfaces driving factors (Figure 15), the climate driving factors (Figure 16), and the anthropogenic activity driving factors (Figure 17), PD, EVI, WSA, and PISI are the most influential factors with the influence of 37.6%, 50.1%, 50.8, and 59.2%, respectively.Overall, each one of the potential driving factors had a comparable influence on SUHI.EVI has been widely used to characterize vegetation coverage.Previous studies have shown that SUHII is negatively correlated with EVI across 419 global big cities [78].From 2002 to 2019, the relative influence of SUHI on EVI is gradually weakening, owing to human activities.The contribution of PISI and NLI was relatively high, indicating that the built-up area and economic development are the main causes of SUHI, while the influence of the climate factors is relatively low during the study period.The results show that the relative influence of SUHI on the climate factors may not be significant.Therefore, the intensification of human activities and economic activities is the main reason for the aggravation of the SUHI effect in Chengdu.As shown in Figure 18, from 2002 to 2019, PD, POP, and GDP increased in Chengdu City and were mainly concentrated in the urban areas.This was mainly due to the increasing centralization of the city, with various industrial zones expanding around the city center.Overall, each one of the potential driving factors had a comparable influence on SUHI.EVI has been widely used to characterize vegetation coverage.Previous studies have shown that SUHII is negatively correlated with EVI across 419 global big cities [78].From 2002 to 2019, the relative influence of SUHI on EVI is gradually weakening, owing to human activities.The contribution of PISI and NLI was relatively high, indicating that the built-up area and economic development are the main causes of SUHI, while the influence of the climate factors is relatively low during the study period.The results show that the relative influence of SUHI on the climate factors may not be significant.Therefore, the intensification of human activities and economic activities is the main reason for the aggravation of the SUHI effect in Chengdu.As shown in Figure 18, from 2002 to 2019, PD, POP, and GDP increased in Chengdu City and were mainly concentrated in the urban areas.This was mainly due to the increasing centralization of the city, with various industrial zones expanding around the city center.Overall, each one of the potential driving factors had a comparable influence on SUHI.EVI has been widely used to characterize vegetation coverage.Previous studies have shown that SUHII is negatively correlated with EVI across 419 global big cities [78].From 2002 to 2019, the relative influence of SUHI on EVI is gradually weakening, owing to human activities.The contribution of PISI and NLI was relatively high, indicating that the built-up area and economic development are the main causes of SUHI, while the influence of the climate factors is relatively low during the study period.The results show that the relative influence of SUHI on the climate factors may not be significant.Therefore, the intensification of human activities and economic activities is the main reason for the aggravation of the SUHI effect in Chengdu.As shown in Figure 18, from 2002 to 2019, PD, POP, and GDP increased in Chengdu City and were mainly concentrated in the urban areas.This was mainly due to the increasing centralization of the city, with various industrial zones expanding around the city center.Overall, each one of the potential driving factors had a comparable influence on SUHI.EVI has been widely used to characterize vegetation coverage.Previous studies have shown that SUHII is negatively correlated with EVI across 419 global big cities [78].From 2002 to 2019, the relative influence of SUHI on EVI is gradually weakening, owing to human activities.The contribution of PISI and NLI was relatively high, indicating that the built-up area and economic development are the main causes of SUHI, while the influence of the climate factors is relatively low during the study period.The results show that the relative influence of SUHI on the climate factors may not be significant.Therefore, the intensification of human activities and economic activities is the main reason for the aggravation of the SUHI effect in Chengdu.As shown in Figure 18, from 2002 to 2019, PD, POP, and GDP increased in Chengdu City and were mainly concentrated in the urban areas.This was mainly due to the increasing centralization of the city, with various industrial zones expanding around the city center.Overall, each one of the potential driving factors had a comparable influence on SUHI.EVI has been widely used to characterize vegetation coverage.Previous studies have shown that SUHII is negatively correlated with EVI across 419 global big cities [78].From 2002 to 2019, the relative influence of SUHI on EVI is gradually weakening, owing to human activities.The contribution of PISI and NLI was relatively high, indicating that the built-up area and economic development are the main causes of SUHI, while the influence of the climate factors is relatively low during the study period.The results show that the relative influence of SUHI on the climate factors may not be significant.Therefore, the intensification of human activities and economic activities is the main reason for the aggravation of the SUHI effect in Chengdu.As shown in Figure 18, from 2002 to 2019, PD, POP, and GDP increased in Chengdu City and were mainly concentrated in the urban areas.This was mainly due to the increasing centralization of the city, with various industrial zones expanding around the city center.

Discussion
An MGWR-CFSDAF spatiotemporal fusion method was proposed to generate high-spatiotemporal-resolution LST data from Landsat, HJ-1B, and MODIS.Although the proposed method can preserve spatial detail and generate high-resolution LST images with high accuracy in Chengdu City, there are also some limitations.Firstly, the performance of the proposed method greatly relies on the pairs of temporally close LST images, which only allows for the clear-sky conditions because the TIR data is difficult to obtain due to cloud cover [79][80][81].If we want to acquire all-weather LST data, more effective cloud removal methods should be adopted to mitigate the influence of clouds.Secondly, since the overpass time of the Landsat, HJ-1B, and MODIS are different, within half an hour, a time normalization method should be applied to correct for possible inconsistencies in the future [82].Thirdly, the spatial distribution of the LST is significantly influenced not only by variations in surface thermal properties but also by a pronounced terrain effect [83].In this study, the spatiotemporal fusion accuracy of the LST data is less affected by mountainous terrain since the study area primarily comprises flat plains.However, this also suggests that further research is needed to explore whether the research method is applicable to urban areas with significant topographic variations.In addition, the main cause of urban thermal environmental change is carbon dioxide (CO2) [84].The spatiotemporal distribution of CO2 emissions has been affected by land use/cover change (LUCC) [85].Deng et al. [86] found that the potential changes in the LST were caused by LUCC.Therefore, in order to mitigate the SUHI effect and meet China's target of carbon neutrality, future studies are needed to explore the relationships between urban expansion, land use changes, CO2 emissions, and the SUHI effect.

Conclusions
This paper took Chengdu, a typical cloudy and rainy city that easily satisfies the missing-filled satellite data scenario, as a case study for SUHII monitoring and performed quantitative analyses to investigate the influence of thirteen potential driving factors on the SUHII from 2002 to 2019.Firstly, high-spatiotemporal-resolution LST dataset with 30 m spatial resolution and 8-day temporal resolution were predicted by the proposed method using an MGWR coupling CFSDAF method.The performance of the method could generate high-accuracy summer LST datasets better than the CFSDAF, and FSDAF methods for the MOD11A2 dataset.Secondly, significantly increasing trends in the SUHII in Chengdu from 2.32 °C in 2002 to 5.07 °C in 2022 were observed for larger suburban areas and increased 2.85 °C during the same period for a smaller suburban area.Finally, PISI, EVI, and NLI are the three most influential factors on SUHI.The total contribution for the driving factors (PISI > EVI > NLI > POP > PD > GDP > WSA > BI > PM2.5 > NDWI > RH > PRE > WS) indicated that the summer SUHI in Chengdu is highly affected by the anthropogenic factor.So, we recommend that the anthropogenic activity driving factor should be considered with CO2 emissions and land use changes for urban planning to mitigate the SUHI effect.

Discussion
An MGWR-CFSDAF spatiotemporal fusion method was proposed to generate highspatiotemporal-resolution LST data from Landsat, HJ-1B, and MODIS.Although the proposed method can preserve spatial detail and generate high-resolution LST images with high accuracy in Chengdu City, there are also some limitations.Firstly, the performance of the proposed method greatly relies on the pairs of temporally close LST images, which only allows for the clear-sky conditions because the TIR data is difficult to obtain due to cloud cover [79][80][81].If we want to acquire all-weather LST data, more effective cloud removal methods should be adopted to mitigate the influence of clouds.Secondly, since the overpass time of the Landsat, HJ-1B, and MODIS are different, within half an hour, a time normalization method should be applied to correct for possible inconsistencies in the future [82].Thirdly, the spatial distribution of the LST is significantly influenced not only by variations in surface thermal properties but also by a pronounced terrain effect [83].In this study, the spatiotemporal fusion accuracy of the LST data is less affected by mountainous terrain since the study area primarily comprises flat plains.However, this also suggests that further research is needed to explore whether the research method is applicable to urban areas with significant topographic variations.In addition, the main cause of urban thermal environmental change is carbon dioxide (CO 2 ) [84].The spatiotemporal distribution of CO 2 emissions has been affected by land use/cover change (LUCC) [85].Deng et al. [86] found that the potential changes in the LST were caused by LUCC.Therefore, in order to mitigate the SUHI effect and meet China's target of carbon neutrality, future studies are needed to explore the relationships between urban expansion, land use changes, CO 2 emissions, and the SUHI effect.

Conclusions
This paper took Chengdu, a typical cloudy and rainy city that easily satisfies the missing-filled satellite data scenario, as a case study for SUHII monitoring and performed quantitative analyses to investigate the influence of thirteen potential driving factors on the SUHII from 2002 to 2019.Firstly, high-spatiotemporal-resolution LST dataset with 30 m spatial resolution and 8-day temporal resolution were predicted by the proposed method using an MGWR coupling CFSDAF method.The performance of the method could generate high-accuracy summer LST datasets better than the CFSDAF, and FSDAF methods for the MOD11A2 dataset.Secondly, significantly increasing trends in the SUHII in Chengdu from 2.32 • C in 2002 to 5.07 • C in 2022 were observed for larger suburban areas and increased 2.85 • C during the same period for a smaller suburban area.Finally, PISI, EVI, and NLI are the three most influential factors on SUHI.The total contribution for the driving factors (PISI > EVI > NLI > POP > PD > GDP > WSA > BI > PM 2.5 > NDWI > RH > PRE > WS) indicated that the summer SUHI in Chengdu is highly affected by the anthropogenic factor.So, we recommend that the anthropogenic activity driving factor should be considered with CO 2 emissions and land use changes for urban planning to mitigate the SUHI effect.
57 E-104 • 20 E, 31 • 15 N-31 • 41 N), has experienced rapid urbanization in the 21st century.In 2022, Chengdu's gross domestic product (GDP) reached 2080 billion US dollars.The population of the city has exceeded 21.2 million people, 15.4 million of them living in urban areas.Rapid urbanization induces significant SUHI effects, especially in summer, which could lead to extreme heat events.From 5 to 24 August 2022, Chengdu experienced a record months-long heatwave, which exceeded 40

Figure 1 .
Figure 1.Location of the study area.

Figure 1 .
Figure 1.Location of the study area.

Figure 2 .
Figure 2. The selected high-spatial-resolution LST data for generating and evaluating the predicted LST results.

Figure 2 .
Figure 2. The selected high-spatial-resolution LST data for generating and evaluating the predicted LST results.

Figure 3 .
Figure 3. Flowchart of testing the proposed method and monitoring daytime SUHI.In the first part (testing the proposed method), both downscaled high-spatial-resolution LST data using the MGWR model and MOD11A1 data captured on 20 April 2013, 16 April 2015, 2 April 2018, and 21 April 2022 were used as the base time (t1) for the proposed method.While the other MOD11A1 data captured on 21 May 2013, 10 July 2015, 5 June 2018, and 7 May 2022 were used as the prediction time (t2) input base data for predicting the high-spatial-resolution LST at t2.In the following part, we call LST directly from MODIS, HJ-1B, Landsat as "observed LST", and LST derived from the proposed method of other spatiotemporal fusion methods as "predicted LST".

Figure 3 .
Figure 3. Flowchart of testing the proposed method and monitoring daytime SUHI.

( 5 ) 4 ) 5 )
If the validation of the MGWR possess is good, then we will perform upscaled LST at 1000 m spatial resolution to 100 m, 120 m, and 300 m.MGWR was used to downscale observed LST from 100 m, 120 m, and 300 m to 30 m.The specific steps of MGWR-based LST downscaling method are shown in Figure 4 and are summarized as follows: Sensors 2023, 23, x FOR PEER REVIEW 10 of 26 LST 100 f NDVI , DEM , Slope , Aspect ∆ (LST 120 f NDVI , DEM , Slope , Aspect ∆ (LST 300 f NDVI , DEM , Slope , Aspect ∆ (6) where LST 100 , LST 120 , and LST 300 are the downscaled LST data at a spatial resolution of 100 m, 120 m, and 300 m, respectively.NDVI 100 , NDVI 120 , NDVI 300 , DEM 100 , DEM 120 , DEM 300 , Slope , Slope , Slope , Aspect 100 , Aspect 120 , and Aspect 300 are the 100 m, 120 m, and 300 m, respectively, environmental variables after spatial aggregation; and ∆ , ∆ , and ∆ are the 100 m, 120 m, and 300 m spatial resolution conversion residual after spatial interpolation.(5) If the validation of the MGWR possess is good, then we will perform upscaled LST at 1000 m spatial resolution to 100 m, 120 m, and 300 m.MGWR was used to downscale observed LST from 100 m, 120 m, and 300 m to 30 m.

Figure 4 .
Figure 4. Flowchart of testing LST downscaling procedure based on MGWR.

Figure 4 .
Figure 4. Flowchart of testing LST downscaling procedure based on MGWR.
June 2002 and 9 May 2006, HJ-1B images were acquired on 20 April 2009, 23 August 2011, 27 April 2012, 20 April 2013, 28 July 2014, 16 April 2015, and 16 May 2016, Landsat 8 images were acquired on 1 May 2017, 2 April 2018, and 11 August 2019, and a Landsat 9 image was acquired on 21 April 2022.The landcover maps from 2002 to 2009 were classified using SVM based on the above cloud-free data.Land cover types are mainly built-up areas, water bodies, vegetation, and bare soil, which are typical in urban areas.

Figure 5 .
Figure 5.The delineation of urban and suburban areas over Chengdu City from 2002 to 2022.

Figure 5 .
Figure 5.The delineation of urban and suburban areas over Chengdu City from 2002 to 2022.
Figure 7c,h,m,r was the MOD11A1 data at t2 on 21 May 2013, 10 July 2015, 5 June 2018, and 7 May 2022, respectively, for predicting the LST data at 100 m and 300 m

Figure 6 .
Figure 6.The 1000 m aggregated Landsat 9 LST on 21 April 2022 of (a) 100 m downscaled LST of vegetation, (b) 100 m downscaled LST of built-up area, (c) 100 m downscaled LST of brae soil, and (d) 100 m downscaled LST of water body.

Figure 9 .
Figure 9.Comparison of the predicted LSTs using the proposed method and CFSDAF, FSDAF.(a) the coefficient of determination (R 2 ); (b) the root mean square error (RMSE); (c) the absolute average difference (AAD).

Figure 9 .
Figure 9.Comparison of the predicted LSTs using the proposed method and CFSDAF, FSDAF.(a) the coefficient of determination (R 2 ); (b) the root mean square error (RMSE); (c) the absolute average difference (AAD).

Figure 10 .
Figure 10.Spatial distributions of 30-m predicted LST using the proposed method.

Figure 10 .
Figure 10.Spatial distributions of 30-m predicted LST using the proposed method.

Figure 11 .
Figure 11.Spatial distribution of the summer averaged 30 m predicted LST using the proposed method from 2002 to 2022.

Figure 11 .
Figure 11.Spatial distribution of the summer averaged 30 m predicted LST using the proposed method from 2002 to 2022.

Figure 12 Figure 12 .
Figure 12 shows the summer SUHII from 2002 to 2022.The significantly increasing trends of the summer SUHII in Chengdu used the averaged 30 m predicted LSTs.The highest SUHII for summer occurred in 2022 (5.07 • C from a larger suburban area and 4.93 • C from a smaller suburban area).Summer SUHII increased from 2.32 • C in 2002 to 5.07 • C for a larger suburban area and increased by 2.85 • C in the same period for a smaller suburban area.This indicates a large SUHII in the summer from 2002 to 2022 over the Chengdu City.It not only increases the risk of heatwave extreme events but also presents a big challenge for scientists to mitigate serious SUHI effects.

Figure 13 .
Figure 13.The relative influence of SUHI of the driving factors in Chengdu City from 2002

Figure 12 .
Figure 12.Temporal changes of SUHII in the study area from 2002 to 2022.

Sensors 2023 , 26 Figure 12 .
Figure 12.Temporal changes of SUHII in the study area from 2002 to 2022.

Figure 13 .
Figure 13.The relative influence of SUHI of the driving factors in Chengdu City from 2002 to 2019.Figure 13.The relative influence of SUHI of the driving factors in Chengdu City from 2002 to 2019.

Figure 13 .
Figure 13.The relative influence of SUHI of the driving factors in Chengdu City from 2002 to 2019.Figure 13.The relative influence of SUHI of the driving factors in Chengdu City from 2002 to 2019.

Figure 14 .
Figure 14.The relative influence of SUHI of the population shift driving factors in Chengdu City from 2002 to 2019.

Figure 15 .
Figure 15.The relative influence of SUHI of the natural land surfaces driving factors in Chengdu City from 2002 to 2019.

Figure 16 .
Figure 16.The relative influence of SUHI of the climate driving factors in Chengdu City from 2002 to 2019.

Figure 17 .
Figure 17.The relative influence of SUHI of the anthropogenic activity driving factors in Chengdu City from 2002 to 2019.

Figure 14 . 26 Figure 14 .
Figure 14.The relative influence of SUHI of the population shift driving factors in Chengdu City from 2002 to 2019.

Figure 15 .
Figure 15.The relative influence of SUHI of the natural land surfaces driving factors in Chengdu City from 2002 to 2019.

Figure 16 .
Figure 16.The relative influence of SUHI of the climate driving factors in Chengdu City from 2002 to 2019.

Figure 17 .
Figure 17.The relative influence of SUHI of the anthropogenic activity driving factors in Chengdu City from 2002 to 2019.

Figure 15 . 26 Figure 14 .
Figure 15.The relative influence of SUHI of the natural land surfaces driving factors in Chengdu City from 2002 to 2019.

Figure 15 .
Figure 15.The relative influence of SUHI of the natural land surfaces driving factors in Chengdu City from 2002 to 2019.

Figure 16 .
Figure 16.The relative influence of SUHI of the climate driving factors in Chengdu City from 2002 to 2019.

Figure 17 .
Figure 17.The relative influence of SUHI of the anthropogenic activity driving factors in Chengdu City from 2002 to 2019.

Figure 16 . 26 Figure 14 .
Figure 16.The relative influence of SUHI of the climate driving factors in Chengdu City from 2002 to 2019.

Figure 15 .
Figure 15.The relative influence of SUHI of the natural land surfaces driving factors in Chengdu City from 2002 to 2019.

Figure 16 .
Figure 16.The relative influence of SUHI of the climate driving factors in Chengdu City from 2002 to 2019.

Figure 17 .
Figure 17.The relative influence of SUHI of the anthropogenic activity driving factors in Chengdu City from 2002 to 2019.

Figure 17 .
Figure 17.The relative influence of SUHI of the anthropogenic activity driving factors in Chengdu City from 2002 to 2019.

Table 1 .
The characteristics of inputs Landsat, HJ-1B, MOD11A1 and MOD11A2 LST satellite LST data of the proposed method.

Table 1 .
The characteristics of inputs Landsat, HJ-1B, MOD11A1 and MOD11A2 LST satellite LST data of the proposed method.

Table 2 .
The potential driving factors selected in this study.
[52,53] a spectral index used to estimate impervious surfaces, such as roads, buildings, and pavements, from the blue band (ρ blue ) and near-infrared band (ρ nir ) of the MOD09A1 data.The formula for PISI can be expressed as follows[52,53]: PISI = 0.8192ρ blue − 0.5735ρ nir +0.075 https://zenodo.org/records/6398971,accessed on 19 February 2023 ChinaHighPM 2.5 is a high-quality dataset in the CHAP series, providing comprehensive, high-res, long-term ground-level air pollutant data for China.Generated using AI and various data sources, it captures spatiotemporal air pollution variations, offering valuable insights into China's air quality.
://www.gis5g.com/,accessed on 16 March 2023 GDP is collected from the Geographic Data Sharing Infrastructure, global resources data cloud, which indicates the economic status within the city.Herein, the GDP density was selected to quantify surface urban heat island.

Table 3 .
Classification accuracy of SVM, ANN and MLC from 2002 to 2022.

Table 4 .
Downscaling statistics for MGWR, GWR, and TsHARP method in this study.

Table 4 .
Downscaling statistics for MGWR, GWR, and TsHARP method in this study.

Table 5 .
Comparison of the predicted LSTs using the proposed method, classical FSDAF and MOD11A2, respectively, with in situ LSTs during the summer from 2002 to 2022.

Table 5 .
Comparison of the predicted LSTs using the proposed method, classical FSDAF and MOD11A2, respectively, with in situ LSTs during the summer from 2002 to 2022.