A High-Resolution Land Surface Temperature Downscaling Method Based on Geographically Weighted Neural Network Regression

: Spatial downscaling is an important approach to obtain high-resolution land surface temperature (LST) for thermal environment research. However, existing downscaling methods are unable to sufﬁciently address both spatial heterogeneity and complex nonlinearity, especially in high-resolution scenes (<120 m), and accordingly limit the representation of regional details and accuracy of temperature inversion. In this study, by integrating normalized difference vegetation index (NDVI), normalized difference building index (NDBI), digital elevation model (DEM), and slope data, a high-resolution surface temperature downscaling method based on geographically neural network weighted regression (GNNWR) was developed to effectively handle the problem of surface temperature downscaling. The results show that the proposed GNNWR model achieved superior downscaling accuracy (maximum R 2 of 0.974 and minimum RMSE of 0.896 ◦ C) compared to widely used methods in four test areas with large differences in topography, landforms, and seasons. We also achieved the best extracted and most detailed spatial textures. Our ﬁndings suggest that GNNWR is a practical method for surface temperature downscaling considering its high accuracy and model performance.


Introduction
High-resolution land surface temperature (LST) products play a crucial role in various types of environmental research. The existing applications mainly include urban heat island effect research [1], climate change monitoring [2], soil moisture estimation [3], surface flux estimation [4], forest fire monitoring [5], drought monitoring [6], surface carbon inversion [7], and geological resource exploration [8]. The acquisition of highquality high-resolution LST data is particularly important for the development of those research areas.
With the development of remote sensing technology, spaceborne thermal infrared sensors are able to obtain data globally [9]. They provide a convenient and effective way to acquire thermal infrared data by remote sensing satellite. However, due to the limitation of thermal infrared sensor performance, it is difficult to obtain high-resolution thermal infrared data for research. As one of the core topics of remote sensing, spatial downscaling can help researchers obtain more high-resolution data, which is widely used in various remote sensing research fields, such as precipitation [10,11], water storage [12][13][14][15][16][17], and LST [18,19]. At present, the downscaling thermal infrared remote sensing image has become the main data source for high-resolution LST and thermal environment research.
The methods of obtaining high-resolution LST based on satellite thermal infrared remote sensing data can be roughly divided into two categories: statistical model-based downscaling methods, and image fusion-based downscaling methods [20]. Common image fusion-based methods include the generalized Laplacian pyramid (GLP) method [21], the co-kriging interpolation method [22], and Bayesian methods [23]. Compared to image fusion-based downscaling methods, statistical model-based downscaling methods for LST are easier to operate, and their downscaling results are accurate within an acceptable range, making them the most widely used. Statistics-based methods are based on the assumption of "relationship scale invariance". The key to successful downscaling is to deal with the spatial non-stationarity and complex nonlinearity in the complex regression relationship of "the LST-surface factors", and common algorithm models include the DisTrad (disaggregation procedure for radiometric surface temperature) downscaling algorithm [24], the TsHARP (Thermal sHARPening) algorithm [25], the Geographically Weighted Regression (GWR) algorithm [26], and the Multi-scale Geographically Weighted Regression (MGWR) downscaling algorithm [27]. In addition, applications of machine learning algorithms into this area are being developed [28], such as the random forest (RF) algorithm [29]. Although those statistics-based methods are widely used, they each have disadvantages. Traditional statistical models struggle to handle non-linear problems, while machine learning models face challenges in handling local spatial non-stationarity [30], which limits their ability to construct expressions for complex natural phenomena with significant spatiotemporal variability of LST, so that the model accuracy and visual effects are not satisfactory. In high-resolution scenes, due to the increase in the number of pixels and the demand for spatial details, the regression relationship between "the LST-surface factors" will become more complex, and the accuracy of the regression results should also be higher. In addition, "the thermal mixture effect" in hot pixels also makes downscaling more difficult [31]. To obtain high accuracy and reasonable high-resolution surface temperature, it is increasingly urgent to develop new and operational models that can handle spatial non-stationarity and complex nonlinearity in LST regression statistics.
Du et al. [32] recently proposed a Geographically Neural Network Weighted Regression (GNNWR) model, which combines the advantages of the Ordinary Liner Regression (OLR) and Neural Network models. By using the learning ability of neural networks and the local spatial interpretation ability of spatial weight, this model can deal with the spatial heterogeneity and complex nonlinearity of regression relationships, which has better fitting accuracy and prediction performance than the GWR and neural network models, with the potential to solve complex spatial relationships in many fields [33][34][35]. Specifically, the GNNWR model uses a spatially weighted neural network (SWNN) to replace the kernel function of GWR, so that the GNNWR model shows excellent interpretation ability in spatial non-stationary modeling, and its application in the study of LST downscaling can provide a new solution for improving the precision of LST downscaling and strengthening the spatial details of downscaling data. Therefore, our goal is to construct a high-resolution LST downscaling model based on GNNWR and demonstrate its superiority. In the model, Normalized Differential Vegetation Index (NDVI), Normalized Difference build-up Index (NDBI), DEM, and Slope serve as covariates to achieve accurate fitting of spatial heterogeneity and nonlinear features in the regression relationship between LST and multiple factors so as to obtain high-precision downscaling results of regional LST. Moreover, the downscaling accuracy of TsHARP, RF, GWR, and MGWR algorithms was quantitatively compared, and the universality of the constructed model was evaluated in different spatiotemporal environments. Moreover, Landsat8 Ready Analysis Data (ARD) were added to compare and evaluate the visual effects of downscaling model results. In other words, this study aims to prove from various aspects that the proposed GNNWR model has the best availability and accuracy, with the best reasonability and practicability in LST downscaling.

Study Area
In this study, Kansas City (study area A) and Phoenix (study area B) were selected as the study areas, the locations of which are shown in Figure 1. Study area A is located in the western part of Missouri, at the confluence of the Missouri River and the Kansas River. It is a typical plain landform with high vegetation coverage, as the mainland is covered by grassland and sparse forest. The climate in this area is temperate continental humid with four distinct seasons, which can represent most temperate areas in the Midwest of the United States. Study area B is located in the middle of Arizona, with a predominant land cover of bare soil. This area features a subtropical desert climate with low rainfall. prove from various aspects that the proposed GNNWR model has the best availability and accuracy, with the best reasonability and practicability in LST downscaling.

Study Area
In this study, Kansas City (study area A) and Phoenix (study area B) were selected as the study areas, the locations of which are shown in Figure 1. Study area A is located in the western part of Missouri, at the confluence of the Missouri River and the Kansas River. It is a typical plain landform with high vegetation coverage, as the mainland is covered by grassland and sparse forest. The climate in this area is temperate continental humid with four distinct seasons, which can represent most temperate areas in the Midwest of the United States. Study area B is located in the middle of Arizona, with a predominant land cover of bare soil. This area features a subtropical desert climate with low rainfall. In addition, February was selected as the winter sample and July the summer sample. Details are shown in Tables 1 and 2. The study areas cover common landform and land cover types, as well as winter and summer seasons, with a large span of surface temperatures, which is conducive to the study of the complex spatiotemporal relationship of downscaling LST. In addition, February was selected as the winter sample and July the summer sample. Details are shown in Tables 1 and 2. The study areas cover common landform and land cover types, as well as winter and summer seasons, with a large span of surface temperatures, which is conducive to the study of the complex spatiotemporal relationship of downscaling LST.

Data
NDVI has been widely used in the downscaling process of LST [36]. DEM and Slope are the key factors reflecting the differentiation of spatial thermal patterns under altitude differences [37], which indicates that the relationship model with LST can be well constructed by integrating factors such as NDVI, DEM, and Slope. Based on the above three factors, this paper adds NDBI as the independent variable input of the model. NDBI is easily affected by such factors as sparse vegetation and bare soil and can fill the deficiency of the representation of impervious surfaces in NDVI. The main factor data sources in this study were Landsat8 red band (R), near-infrared band (NIR), short-wave infrared band (SWIR1), thermal infrared band 10 (TIRS1), and NASA DEM. The experimental data are shown in Table 3, which mainly includes the data name, variable name, spatial resolution, and data source. (1) NDVI is sensitive to vegetation cover, with a value range of [−1, 1], indicating the presence of vegetation cover and increasing with the coverage degree. In this study, it was calculated from the Landsat 8 red band (R) and near-infrared band (NIR). The calculation formula is as follows: (5) Land surface temperature. Obtained from Landsat 8 thermal infrared band 10-TIRS using the mono-window algorithm proposed by Qin [38].
where T represents LST ( • C), T s represents the brightness temperature (/K) of the sensor, T a represents the average atmospheric temperature (/K), and a and b are reference coefficients. When the LST is between 0 and 70 • C, the values of a and b are −67.355351 and 0.458606, respectively. C and D are intermediate variables, which are calculated using the following formula: where ε is the emissivity of the surface, which is obtained according to Qin's empirical formula, and τ is the total atmospheric transmittance from the ground to the sensor, provided by NASA.
Based on the above data, this study aimed to obtain spatial downscaling results of LST with a spatial resolution of 30 m. However, due to the different sources of the above data, problems such as inconsistent projection information and different spatial resolution still existed. To establish the relationship between the spatial resolution of LST and surface factors, all factors were uniformly converted into the same projection coordinate system (WGS_1984_UTM). A 100 m grid was used as the spatial resolution to carry out data matching and the nearest resampling, and the spatial connection method was used to integrate all data into the same grid unit. A group of argument variables with a specification grid as the unit was formed.

Land Surface Temperature Downscaling Model Based on GNNWR
According to the "invariant relationship scale" hypothesis, downscaling of LST generally involves establishing a functional relationship between LST and surface factors at low resolution and then applying this regression function to high-resolution scenarios, along with residual correction, to obtain the predicted results. In this study, the relationship is shown by the following equation: LST 100 = f (NDV I 100 , NDBI 100 , DEM 100 , Slope 100 ) + ∆T 100 (6) where f represents the functional relationship between LST and surface factors at low resolution, and ∆T 100 is the regression residual. Then, the surface factors at high resolution are input into the regression function, and LST 30 is the desired result.
GNNWR obtains a spatial downscaling model of LST by constructing a functional relationship between LST and physical parameters at low spatial resolutions and then applying it to high spatial resolutions. However, unlike previous methods, GNNWR does not need to interpolate the residuals for model calibration [39].
Based on the geographical weighting idea similar to GWR, the GNNWR model considers the spatial variability of the regression relationship as the fluctuation level of the "OLR regression relationship" due to spatial non-stationarity at different locations. In this Remote Sens. 2023, 15, 1740 6 of 18 study, four factors, namely, NDVI, NDBI, DEM, and slope, were used. Therefore, in the LST downscaling experiment of this study, the GNNWR model structure is defined as follows: where s i = (u i , v i ) represents the spatial coordinates of the i-th sample point; the regression coefficients of the OLR model, denoted by β = (β 0 , β 1 , . . . , β 4 ), signify the overall level of the regression relationship between LST and the covariates across the entire study region.
The estimated matrix of the OLR coefficients is expressed as follows: where n represents the number of samples. So, substituting the OLR estimateβ m into Equation (8), we have LST i : where W(s i ) is a (1 + k) × (1 + k) diagonal matrix representing the geographical weights. k is the number of independent variables and can be expressed as: To capture the intricate interplay between spatial distance and spatial weight with precision, GNNWR devises a spatially weighted neural network that expresses the weight kernel function. More specifically, the SWNN leverages the estimated spatial distance between the prediction point and all modeling sample points as its input layer and generates the spatial weight matrix as the output layer, while flexibly adapting to the modeling requirements by selecting an appropriate number of hidden layers. Ultimately, the spatial weight of the prediction point i can be computed through the following formula: where [d s i1 , d s i2 , . . . , d s in ] is the Euclidean spatial distance from the point to the set of samples. Therefore, the LST downscaling model based on GNNWR is defined as shown in Figure 2.

Model Design and Implementation
To improve the optimization efficiency and solution capability of the spatially weighted neural network in the GNNWR model, a neural network architecture and implementation strategy was designed, as shown in Figure 3.

Model Design and Implementation
To improve the optimization efficiency and solution capability of the spatiall weighted neural network in the GNNWR model, a neural network architecture and im plementation strategy was designed, as shown in Figure 3. The spatial weighted neural network uses fully connected networks between eac layer and applies the dropout technique proposed by Srivastava [40] to improve th model's generalization ability. In addition, each hidden layer uses the parameter initial zation method proposed by He et al. [41] and the rectified linear unit (ReLU) activatio function to improve the model's optimization efficiency. At the same time, to reduce th influence of internal covariate shifts, batch normalization technology is introduced int the hidden layer to further improve the calculation ability of the GNNWR model [42 Regarding the number of neurons in each layer, this study utilized a neural networ

Model Design and Implementation
To improve the optimization efficiency and solution capability of the spatially weighted neural network in the GNNWR model, a neural network architecture and implementation strategy was designed, as shown in Figure 3. The spatial weighted neural network uses fully connected networks between each layer and applies the dropout technique proposed by Srivastava [40] to improve the model's generalization ability. In addition, each hidden layer uses the parameter initialization method proposed by He et al. [41] and the rectified linear unit (ReLU) activation function to improve the model's optimization efficiency. At the same time, to reduce the influence of internal covariate shifts, batch normalization technology is introduced into the hidden layer to further improve the calculation ability of the GNNWR model [42]. Regarding the number of neurons in each layer, this study utilized a neural network The spatial weighted neural network uses fully connected networks between each layer and applies the dropout technique proposed by Srivastava [40] to improve the model's generalization ability. In addition, each hidden layer uses the parameter initialization method proposed by He et al. [41] and the rectified linear unit (ReLU) activation function to improve the model's optimization efficiency. At the same time, to reduce the influence of internal covariate shifts, batch normalization technology is introduced into the hidden layer to further improve the calculation ability of the GNNWR model [42]. Regarding the number of neurons in each layer, this study utilized a neural network architecture comprised of seven layers, including one input layer, five hidden layers, and one output layer. The specific parameters used in this architecture are detailed in Table 4.  Figure 4 illustrates the training and validation process of the GNNWR model in which the dataset is partitioned into three subsets in an 8:1:1 ratio, namely, the training set, validation set, and test set, through random selection. The random division of data should make the three data subsets cover the whole research area, thereby minimizing the estimation bias arising from sample selection differences and enhancing the robustness and reliability of the model. Notably, the OLR regression coefficient denoting the average relationship is computed from the training set, which is primarily used to optimize the network weights and biases of the SWNN during the model training process. Moreover, the validation set is employed to monitor the overfitting of the model after each training epoch, whereas the test set is employed to evaluate the predictive performance of the model, both of which are excluded from the model training process [43].  Figure 4 illustrates the training and validation process of the GNNWR model in which the dataset is partitioned into three subsets in an 8:1:1 ratio, namely, the training set, validation set, and test set, through random selection. The random division of data should make the three data subsets cover the whole research area, thereby minimizing the estimation bias arising from sample selection differences and enhancing the robustness and reliability of the model. Notably, the OLR regression coefficient denoting the average relationship is computed from the training set, which is primarily used to optimize the network weights and biases of the SWNN during the model training process. Moreover, the validation set is employed to monitor the overfitting of the model after each training epoch, whereas the test set is employed to evaluate the predictive performance of the model, both of which are excluded from the model training process [43].  In this study, MSE was adopted as the loss function of the GNNWR model, and the MSE of the verification set was used as the overfitting evaluation index [33]. If the overfitting metric exhibits a persistent increasing or plateauing trend and exceeds the predetermined EarlyStop threshold, the model is deemed overfitting and the training process is terminated.

Model Evaluation
This study compared several representative downscaling algorithms for LST, including RF, TsHARP, GWR, and MGWR. Except for TsHARP using single-factor NDVI as the regression kernel, all other algorithms use the same factors as GNNWR. This study applied five commonly used indicators, including the coefficient of determination (R 2 ), root mean square error (RMSE), peak signal-to-noise ratio (PSNR), mean absolute error (MAE), and mean absolute percentage error (MAPE), in evaluating and comparing the spatial downscaling model accuracy of LST.
where y i represents the observed value,ŷ i represents the predicted value, and y i represents the mean value of the observed values.

Comparison and Analysis of Model Performance
This study used five performance evaluation metrics, namely, R 2 , RMSE, MAPE, PSNR, and MAE, to quantitatively evaluate the regression accuracy and fitting effect of the model and to compare the performance of the five models. Table 5 illustrates the specific regression evaluation indicators of the five downscaling algorithm models.
Among the four study areas, the GNNWR model showed higher accuracy, followed by the GWR and MGWR models, but with slight variations. In Area A, GWR performed better than MGWR, while in Area B, MGWR was superior to GWR. The GNNWR model demonstrated improved regression accuracy compared to both the MGWR and GWR models. For instance, compared to the best-performing models in each study area, the R 2 of the GNNWR model respectively increased by 0.025, 0.029, 0.009, and 0.012, and the RMSE respectively decreased by 2.33%, 4.91%, 4.49%, and 9.31%. The R 2 was above 0.9 in all three study areas, and the RMSE was mostly around 1 • C, with values less than 1 • C in three areas. Among them, the best-performing Area B-(2) had the highest R 2 and the smallest RMSE, which were, respectively, 0.974 and 0.896. Additionally, the GNNWR model showed average reductions of 5.73% and 5.34% in MAPE and MAE, respectively, and an increase of 2.445% in PSNR. These three models were superior to the RF and TsHARP models in terms of evaluation metrics. Figure 5 demonstrates the regression evaluation indexes of the GNNWR model, the MGWR model, and the GWR model in detail.
As seen in Figure 5, the GNNWR model had good performance in three areas, namely, A-(1), B-(1), and B-(2), with good R 2 and error metric evaluations. However, all models' R 2 values in the A-(2) region were slightly lower, which may be related to the heterogeneity of the region itself. Similarly, the models in the A region performed worse than those in the B region, indicating that the interpretation of spatial non-stationarity in complex terrains requires the consideration of more factors than in simpler terrains. This also indirectly indicates that the models' ability requirements are higher in complex terrains. Overall, the GNNWR model had the best performance in all indicators, and its performance was less affected by perturbations in different study areas, indicating good stability and an excellent capturing ability for spatial non-stationarity and complex nonlinearity. Among the four study areas, the GNNWR model showed higher accuracy, followed by the GWR and MGWR models, but with slight variations. In Area A, GWR performed better than MGWR, while in Area B, MGWR was superior to GWR. The GNNWR model demonstrated improved regression accuracy compared to both the MGWR and GWR models. For instance, compared to the best-performing models in each study area, the R 2 of the GNNWR model respectively increased by 0.025, 0.029, 0.009, and 0.012, and the RMSE respectively decreased by 2.33%, 4.91%, 4.49%, and 9.31%. The R 2 was above 0.9 in all three study areas, and the RMSE was mostly around 1 °C, with values less than 1 °C in three areas. Among them, the best-performing Area B-(2) had the highest R 2 and the smallest RMSE, which were, respectively, 0.974 and 0.896. Additionally, the GNNWR model showed average reductions of 5.73% and 5.34% in MAPE and MAE, respectively, and an increase of 2.445% in PSNR. These three models were superior to the RF and TsHARP models in terms of evaluation metrics. Figure 5 demonstrates the regression evaluation indexes of the GNNWR model, the MGWR model, and the GWR model in detail.  In terms of computation time, the order of models from fastest to slowest was TsHARP < RF < GWR < GNNWR << MGWR. However, except for GNNWR, all other models still required residual transformation interpolation to complete the downscaling process. Therefore, the total downscaling time for GNNWR was much shorter than that for the other models.
The 30 m GNNWR model regression results were resampled to 100 m and compared with the inverted LST data to analyze the absolute errors in each region. The error pixels were classified into 5 levels, namely, 0-1 • C, 1-2 • C, 2-3 • C, 3-4 • C, and >4 • C, and the proportion of pixels in each error level was calculated. The results are shown in Figure 6.
complex terrains requires the consideration of more factors than in simpler terrains. This also indirectly indicates that the models' ability requirements are higher in complex terrains. Overall, the GNNWR model had the best performance in all indicators, and its performance was less affected by perturbations in different study areas, indicating good stability and an excellent capturing ability for spatial non-stationarity and complex nonlinearity.
In terms of computation time, the order of models from fastest to slowest was TsHARP < RF < GWR < GNNWR << MGWR. However, except for GNNWR, all other models still required residual transformation interpolation to complete the downscaling process. Therefore, the total downscaling time for GNNWR was much shorter than that for the other models.
The 30 m GNNWR model regression results were resampled to 100 m and compared with the inverted LST data to analyze the absolute errors in each region. The error pixels were classified into 5 levels, namely, 0-1 °C, 1-2 °C, 2-3 °C, 3-4 °C, and >4 °C, and the proportion of pixels in each error level was calculated. The results are shown in Figure 6.  Table 6, it can be seen that around 60% of the pixels in the four study areas had errors within the 0-1 °C range, indicating that GNNWR had good accuracy in downscaling the surface temperature in different topographic regions and seasons. Errors >2 °C accounted for 12.66%, 5.29%, 3.44%, and 14.54% of the total pixels in the four study areas and were mainly distributed in areas with intensive human activities, such as urban and riverside regions, as well as mountainous regions where the model algorithm's outliers  Table 6, it can be seen that around 60% of the pixels in the four study areas had errors within the 0-1 • C range, indicating that GNNWR had good accuracy in downscaling the surface temperature in different topographic regions and seasons. Errors >2 • C accounted for 12.66%, 5.29%, 3.44%, and 14.54% of the total pixels in the four study areas and were mainly distributed in areas with intensive human activities, such as urban and riverside regions, as well as mountainous regions where the model algorithm's outliers were mainly located on the mountain ridges and valleys due to differences in elevation. In addition, the overall proportion of absolute error pixels in study area A-(2) was relatively good. However, combined with Figure 5, it can be seen that this area had a relatively large MAPE, indicating a significant level of error. The main reason is that the explanatory power of the independent variables in study area A-(2) was insufficient, such as the NDVI values not showing significant differences in winter, and trees being not as lush as in summer. Therefore, it may be necessary to replace or add more reasonable regression covariates to the model.

Spatial Distribution of High-Resolution Land Surface Temperature
To better analyze the downscaling results, this study compared the visual effects of Landsat8 ARD surface temperature product data, GNNWR downscaling results, and other models that performed well in each region in the four selected study areas. The Landsat 8 ARD surface temperature product is a high-resolution dataset (30 m) generated by the USGS through the reanalysis of Landsat 8 thermal infrared data applying georeferencing, atmospheric correction, and resampling techniques. Figure 7 shows the overall comparison of surface temperature in the four main study areas.

Spatial Distribution of High-Resolution Land Surface Temperature
To better analyze the downscaling results, this study compared the visual effects of Landsat8 ARD surface temperature product data, GNNWR downscaling results, and other models that performed well in each region in the four selected study areas. The Landsat 8 ARD surface temperature product is a high-resolution dataset (30 m) generated by the USGS through the reanalysis of Landsat 8 thermal infrared data applying georeferencing, atmospheric correction, and resampling techniques. Figure 7 shows the overall comparison of surface temperature in the four main study areas.  It can be seen from Figure 7 that the GNNWR downscaling results were in line with the original information of most LST, with more detailed and richer spatial textures. This indicates that while the model method achieved good results, it also maintained the consistency of the thermal characteristics distribution before and after downscaling. Especially in the B-(1) research area with flat and single topography in summer and the B-(2) research area with high altitude difference and single topography in summer, the regression results of the GNNWR downscaling model had good and reasonable visual detail effects. Compared with the ARD LST product, the GNNWR downscaling results had better spatial details, and the ARD LST product was an obvious resampled product. Compared with the better models in each study area, the LST distribution of the GNNWR downscaling results was more reasonable. For example, in the B-(1) and B-(2) regions, GNNWR did not have errors as did the MGWR model, which presented a large-scale data overestimation. In addition, GNNWR's results showed that in study area A, high-temperature areas are mainly concentrated in urban areas and higher than surrounding areas, indicating obvious urban heat island effects, while in study area B, due to the characteristics of desert areas, low temperatures are mainly concentrated around urban or high-altitude areas, which, in these areas mentioned above, shows that the LST obtained from GNNWR was basically consistent with satellite-retrieved LST.

Availability and Advantages of the GNNWR Model
Considerable progress has been made in satellite data downscaling. However, the current methods based on statistical models cannot simultaneously deal with the spatial non-stationarity and complex nonlinearity of the LST-surface factors relationship, which makes it difficult to obtain reasonable and accurate downscaling results, especially in the region of complex topography. Therefore, a GNNWR-based LST downscaling method was developed in this study to obtain high-precision and reasonable LST data.
It can be seen from Table 5 that the GNNWR model has the best performance, followed by the GWR and MGWR model, which have their advantages in each region and are far superior to the TsHARP and RF models. In particular, the bandwidth of each factor in the results of the MGWR model in study area A is equal, indicating that MGWR degenerates into GWR in complex landforms, which is not fully applicable to downscaling applications in this case because of its instability. At the same time, GNNWR has higher index accuracy than GWR and MGWR. For example, RMSE, MAPE, and MAE of the GNNWR model decreased by 5.26%, 5.73%, 5.34%, respectively, and PSNR increased by 2.445% on average, and the maximum R 2 reached 0.974. All these indicate that GNNWR has superior performance in constructing the LST and surface factors relationship. This also shows that the SWNN weight matrix of GNNWR is superior to the kernel function of GWR in space non-stationarity expression. In addition, in the comparison of image vision, the GNNWR model has a more accurate description and better expression of spatial details in the expression of ground object features and is more reasonable locally, which further proves the applicability of GNNWR in LST downscaling.
Using the GNNWR model for downscaling, reasonable design of deep learning model parameters can achieve excellent results, and it is more operational and theoretical, not only expressing the ability to cope with complex spatial non-stationarity but also compatible with the powerful computing power of deep learning. Therefore, it can be used for further research and testing of multi-region, multi-factor, and multi-data sets for common surface temperature downscaling scenarios based on "relationship scale invariance" to improve the universality and accuracy of model methods.

Veracity and Reasonability of LST Mappings
To better display the experimental results and assess the accuracy and reasonability of LST mappings of GNNWR, Figure 8 shows the comparison details in each area.
Upon examining the details in each study area, the thermal spatial distribution of each model or product is basically consistent. However, it is evident that the GNNWR model exhibits more reasonable and detailed temperature distributions for most land surfaces compared to other high-performing downscaling models. The GNNWR model outperforms other high-performing downscaling models in terms of not producing unrealistic temperature distributions, smooth color patterns, overall numerical abnormalities, or failing to express the geographic features of roads, buildings, and dense forest areas in complex terrain, such as rivers in region A-(1) and ridges in region B- (2). Although the high-performing GWR and MGWR models have achieved downscaling effects, they sometimes exhibit indistinctive and overly smooth land features, lack of details, and large-scale temperature estimation deviations. Specifically, in A-(1) and B-(1) regions, the regression results of the GWR model and the MGWR model show obvious "overestimation" in the river part and the plain part, respectively, while the regression results of the MGWR model show obvious "underestimation" in the hillside part in B- (2). Moreover, in the A-(2) region, a large colorful glossy area appears. At the same time, these unreasonable distributions also illustrate GNNWR's excellent regression prediction ability and high estimation accuracy. Upon examining the details in each study area, the thermal spatial distribution of each model or product is basically consistent. However, it is evident that the GNNWR model exhibits more reasonable and detailed temperature distributions for most land surfaces compared to other high-performing downscaling models. The GNNWR model outperforms other high-performing downscaling models in terms of not producing unrealistic temperature distributions, smooth color patterns, overall numerical abnormalities, or failing to express the geographic features of roads, buildings, and dense forest areas in complex terrain, such as rivers in region A-(1) and ridges in region B- (2). Although the high-performing GWR and MGWR models have achieved downscaling effects, they sometimes exhibit indistinctive and overly smooth land features, lack of details, and largescale temperature estimation deviations. Specifically, in A-(1) and B-(1) regions, the regression results of the GWR model and the MGWR model show obvious "overestimation" in the river part and the plain part, respectively, while the regression results of the MGWR model show obvious "underestimation" in the hillside part in B- (2). Moreover, in the A-(2) region, a large colorful glossy area appears. At the same time, these unreasonable distributions also illustrate GNNWR's excellent regression prediction ability and high estimation accuracy. Figure 9 visualizes the absolute errors of the regression results of the GNNWR model and the model with better performance in each study area at a spatial scale of 100 m.  From Figure 9, it can be observed that the errors of each model are mainly concentrated in densely populated areas such as cities and riverbanks, as well as in mountainous areas such as ridges, and the error proportions are relatively high in the GWR and MGWR models, which are consistent with our previous findings. In addition, the PSNR index of GNNWR in Table 5 is mainly between 30 and 40 and has the highest level, indicating that the GNNWR model also performs better in terms of quantitative visual accuracy. In summary, compared to other models, GNNWR results are more reasonable and accurate. From Figure 9, it can be observed that the errors of each model are mainly concentrated in densely populated areas such as cities and riverbanks, as well as in mountainous areas such as ridges, and the error proportions are relatively high in the GWR and MGWR models, which are consistent with our previous findings. In addition, the PSNR index of GNNWR in Table 5 is mainly between 30 and 40 and has the highest level, indicating that the GNNWR model also performs better in terms of quantitative visual accuracy. In summary, compared to other models, GNNWR results are more reasonable and accurate.

Model Applicability of Different Factors and Sensors
It can be seen from Table 5 that in the winter complex landform, the effect of NDVI weakens, and the regression fitting accuracy of the model decreases, which indicates that the effectiveness and selection of regression factors may lead to changes in the model's ability to explain the real LST. After eliminating the interference of the autocorrelation among the factors, the model's expressive power can be enhanced by adding suitable surface factors for complex terrain. In areas with dense vegetation, NDVI has a greater impact on the results during the downscaling process. In urban areas, NDBI plays a greater role in reflecting the surface dryness index. After canceling the NDBI factor, this study also conducted GNNWR downscaling analysis on the study area A-(1), as shown in Figure 10. Comparing the downscaling temperature prediction results in Table 5, it can be seen that, in addition to the commonly used factors of elevation and slope, the model's accuracy using the combined NDBI + NDVI is higher than that using only NDVI. Similarly, for different complex terrains and seasons, using surface factors that can represent distinct features of the area, such as slope direction, the normalized difference water index (NDWI), and the normalized difference bare soil index (NDBSI), endows the model with better interpretability for specific regions. As for data sources, currently, data from different satellite sensors are commonly used to obtain high temporal and spatial resolution surface temperatures. The imaging time and imaging methods of different sensors differ, resulting in differences in the produced surface temperature. Therefore, when establishing a regression relationship at a lower resolution, it will inevitably be affected by different variables. For example, when downscaling using MODIS data, the LST values are derived from the Terra/MODIS sensor, while the NDVI values are obtained from the OLI sensor of Landsat 8. As the two sensors have significant differences in various aspects, such discrepancies will result in errors in the final downscaling products. However, from the perspective of the regression model alone, the accuracy of the GNNWR model is still higher than traditional downscaling methods and machine learning methods. The next step is to apply the model to downscaling between different sensors.

Conclusions
By integrating four land surface factors as the independent variables in different sea- As for data sources, currently, data from different satellite sensors are commonly used to obtain high temporal and spatial resolution surface temperatures. The imaging time and imaging methods of different sensors differ, resulting in differences in the produced surface temperature. Therefore, when establishing a regression relationship at a lower resolution, it will inevitably be affected by different variables. For example, when downscaling using MODIS data, the LST values are derived from the Terra/MODIS sensor, while the NDVI values are obtained from the OLI sensor of Landsat 8. As the two sensors have significant differences in various aspects, such discrepancies will result in errors in the final downscaling products. However, from the perspective of the regression model alone, the accuracy of the GNNWR model is still higher than traditional downscaling methods and machine learning methods. The next step is to apply the model to downscaling between different sensors.

Conclusions
By integrating four land surface factors as the independent variables in different seasons, elevations, and land cover types in the research area, this study proposes a highresolution (<120 m) LST downscaling method based on GNNWR, which can well represent spatial heterogeneity and complex nonlinearity. The modeling accuracy is compared with four commonly used downscaling methods, and the visual effect is compared with the Landsat 8 ARD surface temperature product.
The results show that based on the comparison of visual effects and five evaluation metrics, the prediction results of the GNNWR model showed significant improvement in all metrics, and the mapping accuracy was also the highest with best spatial details in the thermal pattern distribution, especially in the summer season and under a single land cover type (maximum R 2 of 0.974 and minimum RMSE of 0.896 • C). This suggests that the GNNWR model, which combines the advantages of traditional geographical models for spatial interpretation and machine learning models for non-linear fitting, performs well in high-resolution downscaling applications. Different seasons, land cover types, elevations, and the random mechanism of the training set have certain impacts on the GNNWR model construction accuracy, but its overall effect is better than that of the widely used models, which means that the GNNWR model has good stability. On the other hand, in the selection of land surface factors, NDVI is still the main factor that affects the spatial distribution of surface temperature inversion results. After excluding the interference of autocorrelation between each factor, different regression factors can be used according to various research areas, such as adding the slope direction factor in mountainous areas and adding NDWI in areas with more water bodies, which improves the model fitting degree. In the future, we will continue to explore the impact of specific land cover types on the model.
Overall, the GNNWR-based downscaling method proposed in this study can improve the downscaling accuracy and spatial texture information of the original low spatial resolution surface temperature while ensuring spatial consistency. For the relevant research fields involved, it can provide researchers with a more convenient, fast, and accurate way to obtain high-resolution surface temperature.