Spatial Downscaling of Nighttime Land Surface Temperature Based on Geographically Neural Network Weighted Regression Kriging

Abstract

Land surface temperature (LST) has a wide application in Earth Science-related fields, and spatial downscaling is an important method to retrieve high-resolution LST data. However, existing LST downscaling methods have difficulties in simultaneously constructing and expressing spatial non-stationarity, spatial autocorrelation, and complex non-linearity during the LST downscaling process, which limits the performance of the models. Moreover, there is a lack of research on high-resolution nighttime land surface temperature (NLST) reconstruction based on spatial downscaling, which does not meet the data needs for urban-scale nighttime urban heat island (UHI) studies. Therefore, this study combined Geographically Neural Network Weighted Regression (GNNWR) with Area-to-Point Kriging interpolation (ATPK) to propose a Geographically Neural Network Weighted Regression Kriging (GNNWRK) model for NLST downscaling. To verify the model’s generality and robustness, this study selected four study areas with different landform and climate type for NLST spatial downscaling experiments. The GNNWRK was compared with four benchmark downscaling methods, including TsHARP, Random Forest, Geographically Weighted Regression, and GNNWR. The results show that compared to these four benchmark methods, the GNNWRK method has higher accuracy in NLST downscaling, with a maximum Pearson’s Correlation Coefficient (Pcc) of 0.930 and a minimum Root Mean Square Error (RMSE) of 0.886 K. Moreover, the validation based on MODIS NLST data and ground-measured NLST data also indicates that the GNNWRK model can obtain more accurate, high-resolution NLST with richer and more detailed texture. This enhances the potential of NLST in studying the effects of urban nighttime heat islands at a finer scale.

1. Introduction

Land surface temperature (LST) is a crucial indicator and key parameter for characterizing the dynamic balance of surface energy, describing the process of energy exchange at the earth’s surface, and studying climate change [1]. It is widely used in climate change monitoring [2], vegetation dynamics monitoring [3,4], water resources management [5,6], soil moisture estimation [7], and surface carbon retrieval [8], among other applications [9,10].

Nighttime land surface temperatures (NLST) are generally lower than their daytime counterparts, primarily due to the absence of solar radiation, which allows the land surface to lose heat through longwave radiation [11]. This cooling effect is influenced by land surface material, with impervious surfaces retaining more heat compared to vegetated or natural surfaces, resulting in different spatial patterns of LSTs during the daytime and nighttime [12,13,14]. In addition, the mechanisms which drive LST are also allegedly different between day and night [15,16,17]. Given that NLST has a major effect on the urban heat island and human health [18,19,20], the lack of high-resolution NLST limits our understanding.

While ground measurement stations can provide relatively accurate land surface temperature, their limited number and sparse spatial distribution cannot offer continuous LST data over large areas. Spaceborne thermal infrared sensors can provide continuous high-resolution thermal infrared data globally, offering a more operable and effective way to retrieve LST data [21,22].

At present, there is a spatiotemporal resolution trade-off of the satellite-derived LST data [23]. To avoid increasing satellite costs, spatial downscaling of spaceborne thermal infrared data has been extensively employed to address the trade-off issue between spatial and temporal resolution. Spatial downscaling technology can make predictions in higher spatial resolution than the original image, i.e., extrapolating information from a large spatial range to a smaller one, thereby obtaining more information-rich remote sensing images [24]. Therefore, spatial downscaling of spaceborne thermal infrared remote sensing imagery can yield more detailed and continuous LST data.

The methods for downscaling LST can be classified into two groups, namely image fusion-based methods and kernel-based methods [25,26]. Kernel-based downscaling methods are based on the assumption of “relationship scale invariance”, which posits that the regression relationship between LST and various kernels (auxiliary factors) is scale-invariant [27]. The kernel-based methods typically build a regression relationship of low-resolution LST and auxiliary factors at first, subsequently introducing auxiliary factors of high-resolution into the regression. After residual correction, a high-resolution LST is generated [28]. Therefore, the accuracy of modeling the regression of between “LST and auxiliary factors” determines the LST downscaling method’s performance.

Classical downscaling methods were often focused on the linear statistical relationship between the vegetation index and LST [29,30]. Later, machine learning methods such as random forests were introduced to construct the nonlinear relationships between LST and auxiliary factors [28,31]. The “spatial invariance” assumption of global regression models has limitations in geographical studies [32]. Geographically Weighted Regression (GWR) has found extensive applications in urban heat island research and LST downscaling, as it can measure the spatial non-stationarity of LST [33]. Due to the insufficient ability of classical GWR to express the local characteristics of spatial relationships, the Multi-scale Geographically Weighted Regression (MGWR) model has been employed to construct multi-scale spatial non-stationary relationships. Methods based on MGWR have shown superior downscaling performance compared to those GWR-based methods [34].

Du et al. [35] introduced a Geographically Neural Network Weighted Regression (GNNWR) model, which employs Spatially Weighted Neural Networks (SWNN) as a novel approach, diverging from the traditional kernel function utilized in GWR. This enhances the interpretability of the GNNWR model in the context of spatial non-stationarity modeling [36,37,38]. Although it has been successfully applied to daytime LST downscaling (Liang et al. [39]), the inherent spatial autocorrelation of the data has not been sufficiently considered. Spatial non-stationarity and spatial autocorrelation are the two core effects that need to be considered in the context of modelling spatial data [40], neglecting spatial autocorrelation limits the accuracy of spatial modelling. Area-to-point Kriging (ATPK) is a geostatistical interpolation and downscaling method based on spatial autocorrelation [41]. Using ATPK to interpolate and reconstruct high-resolution regression residuals has the potential to enhance the accuracy of downscaling methods. Currently, it is extensively used in various domains for downscaling, including LST and soil moisture [42,43,44,45].

Therefore, this study combined GNNWR with ATPK to develop a method called GNNWRK that takes spatial non-stationarity, spatial autocorrelation, and nonlinearity into account, downscaling the NLST from a spatial resolution of 1000 m to 100 m. In light of data relevance and accessibility, we conducted the downscaling study using Landsat data, ASTER NLST data, and OSM road network data, and we quantitatively compared our results with four benchmark downscaling methods including TsSHARP, RF, GWR, and GNNWR. We used the ASTER NLST data product AST_08 and nighttime ground measured LST data to assess the accuracy of the GNNWRK method. Moreover, the generalizability of GNNWRK method was evaluated in different spatiotemporal environments.

2. Materials

2.1. Study Area

To emphasize the generalizability of the downscaling model in complex scenarios, this paper incorporated considerations like varying landcover and climate types, and we selected Tucson (referred to as study area A) and Indianapolis (referred to as study area B) as the study areas. Study area A is the second-largest city in Arizona, bordered by mountains on all sides, with the urban area situated among dry riverbeds and alluvial plains. The region has a subtropical desert climate, with a dry environment and abundant sunlight. Study Area B is located in the central region of Indiana and is characterized by plains. The surface is predominantly covered with farmland and sparse forests, with a high level of vegetation cover. The climate is a humid continental climate, distinguished by four distinct seasons. It serves as a typical representative of the temperate Great Plains region in the Midwest of the United States. The geographical locations of study areas are shown in the Figure 1.

2.2. Data

The study employed Landsat-8/9 imagery and ASTER NLST to assess the effectiveness and superiority of the GNNWRK method. Landsat imagery at a preprocessing level of L1 was used to extract the auxiliary factors for downscaling, and ASTER NLST served to simulate MODIS NLST and to verify the accuracy of the downscaled results. Moreover, the auxiliary factor of road network kernel density was derived by computing the road network data obtained from Open Street Map (https://www.openstreetmap.org, accessed on 1 December 2023). In addition, to further evaluate the performance of GNNWRK in downscaling real MODIS NLST, we also used MOD11A1 NLST and in situ LST measured at the ground stations. All the satellite data used are summarized in

Table 1. Characteristics of Landsat, MODIS, and ASTER satellite imageries used in this study.

Study Area	Code	Landsat8/9	MODIS	ASTER
Tucson, Arizona, USA	A-(1)	13 October 2023	16 October 2023	16 October 2019
				16 October 2023
	A-(2)	06 May 2023	08 May 2023	26 May 2022
				08 May 2023
	A-(3)	18 December 2018	09 December 2018	17 December 2015
				09 December 2018
	A-(4)	27 May 2019	25 May 2019	13 March 2016
				25 May 2019
	A-(5)	25 June 2018	02 July 2018	13 June 2017
				02 July 2018
Indianapolis, Indiana, USA	B-(1)	15 September 2022	15 September 2022	08 September 2022
				15 September 2022
	B-(2)	03 June 2022	04 June 2022	17 June 2021
				04 June 2022

, with comprehensive descriptions provided in the subsequent sections.

2.2.1. Spaceborne Thermal Infrared Sensors and NLST Data

The MOD11A1 LST product was used as the low-resolution NLST for downscaling in this study. The MOD11A1 dataset is a daily land surface temperature/emissivity (LST/E) data product derived from the MODIS sensors on Terra satellite, which includes both daytime and nighttime LSTs, with a spatial resolution of 1000 m [46]. The MOD11A1 LST product used the split-window algorithm to retrieve LST from thermal infrared sensors [47]. After validation with the ground measured LST obtained by field and weather stations, the RMSE of MODIS LSTs at most stations was less than 2 K. This study used MOD11A1 LST obtained from MODIS LST Collection 6.1.

The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) is a spaceborne sensor mounted on the Terra satellite. This study utilizes ASTER NLST as high-resolution NLST data. ASTER was capable of capturing LST data through its thermal infrared bands, providing a fine spatial resolution of 90 m. The AST_08 surface kinetic temperature product was generated based on ASTER data through the Temperature/Emissivity Separation (TES) algorithm and maintained the same spatial resolution of 90 m. ASTER LST has been used for the validation analysis of MODIS LST, and cross-satellite comparison studies have found that they have a low error range (0.2–1.5 K) [48].

2.2.2. Auxiliary Factors

Synthesizing the existing research, the retrieval of spectral indexes based on Landsat multispectral high-resolution imagery was an important means of obtaining auxiliary factors in the process of downscaling LST [25,49,50]. Among these, the Normalized Difference Vegetation Index (NDVI) has been widely used in studies of LST downscaling [29]. Humidity (WET) and the Normalized Multi-band Drought Index (NMDI) have also been used to characterize complex land cover conditions [25]. Impervious surfaces help enhance the accuracy of LST downscaling in urban areas, therefore, the use of the Normalized Difference Built-up Index (NDBI) and road network data provided by Open Street Map supplements the limitations of NDVI in characterizing impervious surfaces [51].

The NDVI represents vegetation cover which is closely related to land surface thermal properties. The index is derived from the Red(R) and Near-Infrared (NIR) bands of Landsat using the following formula:

$$NDVI=\frac{NIR-R}{NIR+R}$$

(1)

The NDBI is influenced by the coverage of impervious surfaces and is capable of reflecting information regarding human-made impervious surfaces such as urban built-up areas. This index is derived from the NIR and Short Wavelength Infrared (SWIR1) bands of Landsat.

$$NDBI=\frac{SWIR1-NIR}{SWIR1+NIR}$$

(2)

WET is one of the three feature components derived from tasseled cap transformation, and it is closely associated with the moisture content of vegetation and soil. In this study, WET is calculated using data from Landsat 8/9.

$$\begin{gathered} WET=0.1509\times BLUE+0.1973\times GREEN+0.3279\times RED+0.3406\times NIR \\ -0.7112\times SWIR1-0.4572\times SWIR2 \end{gathered}$$

(3)

The soil moisture content and its spatial distribution plays a crucial role in influencing LST, which can be expressed by NMDI. In this study, NMDI is calculated using data from Landsat 8/9.

$$NMDI=\frac{\left(NIR-\left(SWIR1-SWIR2\right)\right)}{\left(NIR+\left(SWIR1-SWIR2\right)\right)}$$

(4)

Auxiliary factors that directly reflect the thermal characteristics of the land surface can provide useful information about land cover from a thermal perspective. Dong et al. [26] used the daytime LST derived from Landsat at a previous time(t1) as the auxiliary factors for the downscaling model at the target time(t2), developing a hybrid framework with downscaling performance superior to classic kernel-based regression models and spatiotemporal fusion models. Considering that surface factors derived from daytime satellite might not accurately represent the land surface thermal properties at night, Yoo et al. [48] used ASTER LST data as input for the local linear forest model, downscaling the MODIS NLST data from an original spatial resolution of 1000 m to a finer resolution of 250 m. Similarly, we introduced the ASTER NLST from a previous time (t1) as an auxiliary factor LST_t1 to establish the “LST-auxiliary factor” regression relationship for the target time (t2).

2.2.3. In-Site LST Data

This study used ground temperature data measured at sites provided by the U.S. Climate Reference Network (USCRN) for comparative analysis with the downscaled NLSTs. The USCRN is an extensive data collection resource that provides comprehensive meteorological and climatic information from 114 sites. These sites are developed, deployed, monitored, and maintained by the National Oceanic and Atmospheric Administration (NOAA). This includes data on air temperature, relative humidity, and surface radiative temperature on an hourly basis [52]. USCRN sites use Apogee Instruments IRTS-P infrared temperature sensors to monitor the ground temperature of thermal infrared radiation, and its data has been widely used in the ground validation of spaceborne thermal infrared remote sensing data [21,53,54].

Study area A has a USCRN site located on a hill west of Tucson city, with the coordinates of 111.17°W, 32.24°N. The data were obtained through the National Centers for Environmental Information (NCEI).

3. Methods

3.1. GNNWR

Classical GWR uses a kernel function to calculate spatial weight matrices, but its structure is relatively simple and struggles to express the complex relationships in complicated geographic processes. Du et al. [35] introduced an artificial neural network model to express the complex non-linear effects inherent in spatial non-stationary relationships, proposing the GNNWR model. This model employed artificial neural networks to calculate spatial weight matrices, significantly enhancing the model’s ability to express spatial non-stationary relationships. The GNNWR is defined as:

$$y_i=w_0\left(u_i,v_i\right)\times {\beta }_0\left(u_i,v_i\right)+\sum _{k=1}^n{w}_k\left(u_i,v_i\right)\times {\beta }_k{x}_{ik}{+\epsilon }_i\hspace{1em}i=1,2,\cdots ,n$$

(5)

The above expression is derived from the classical GWR model, the regression constant term ${\beta }_0\left(u_i,v_i\right)$ is transformed into $w_0\left(u_i,v_i\right)\times {\beta }_0\left(u_i,v_i\right)$, while $w_k\left(u_i,v_i\right)\times {\beta }_k{x}_{ik}$ is equivalent to the regression coefficient ${\beta }_k\left(u_i,v_i\right)$ in the GWR. Therefore, $w_0\left(u_i,v_i\right)$ represents the spatial non-stationary weight of ${\beta }_0$, and the non-stationary weight of the regression coefficient ${\beta }_k$ is represented as $w_k\left(u_i,v_i\right)$. By substituting the estimated value ${\beta }_k$ of ${\widehat{\beta }}_k$ obtained from the OLR model into the above formula, the mathematical expression for the estimated value ${\widehat{y}}_i$ of $y_i$ is:

$${\widehat{y}}_i=\sum _{k=0}^p{\widehat{\beta }}_k\left(u_i,v_i\right)x_{ik}=\sum _{k=0}^p{w}_k\left(u_i,v_i\right)\times {\widehat{\beta }}_k\left(OLR\right)x_{ik}$$

(6)

If converted into matrix form, it is as follows:

$${\widehat{y}}_i=x_i^T\cdot \widehat{\beta }\left(u_i,v_i\right)=x_i^T{W\left(u_i,v_i\right)\left(X^{T}X\right)}^{-1}{X}^{T}y$$

(7)

$W(u_i,v_i)$ is the spatial non-stationary weights matrix at sample point i, and its mathematical expression is as follows:

$$W\left(u_i,v_i\right)=\left[\begin{array}{cccc}{w}_0\left(u_i,v_i\right)& 0& \cdots & 0\\ 0& w_1\left(u_i,v_i\right)& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & w_p\left(u_i,v_i\right)\end{array}\right]$$

(8)

For one point $i$, the spatial weight should be jointly determined by the spatial distance of that point to all known sample points. Therefore, the overall non-linear relationship between the target point $i$ and the positions of n sample points will affect the solving precision of the spatial weight matrix. GNNWR has designed a type of SWNN to execute the neural network expression of the weight kernel function; its mathematical model is defined as follows:

$$W\left(u_i,v_i\right)=SWNN\left({\left[d_{i1}^s,d_{i2}^s,\cdots ,d_{in}^s\right]}^T\right)$$

(9)

$\left[d_{i1}^s,d_{i2}^s,\cdots ,d_{in}^s\right]$ represents the spatial proximity characterization of target point i, manifested as the spatial distance vector from target point i to the other sample points. The output generated by the SWNN is the spatial weight matrix $W\left(u_i,v_i\right)$ at point i. By substituting the weight matrix into the equation of GWR, the estimated values of the dependent variable y for each sample point, denoted as $\widehat{y}$ can be expressed as follows:

$$\widehat{y}=\left[\begin{array}{c}{\widehat{y}}_1\\ {\widehat{y}}_2\\ ⋮\\ {\widehat{y}}_n\end{array}\right]=\left[\begin{array}{c}{x}_1^T{SWNN\left({\left[d_{11}^s,d_{12}^s,\cdots ,d_{1n}^s\right]}^T\right)\left(X^{T}X\right)}^{-1}{X}^T\\ x_2^T{SWNN\left({\left[d_{21}^s,d_{22}^s,\cdots ,d_{2n}^s\right]}^T\right)\left(X^{T}X\right)}^{-1}{X}^T\\ ⋮\\ x_n^T{SWNN\left({\left[d_{n1}^s,d_{n2}^s,\cdots ,d_{nn}^s\right]}^T\right)\left(X^{T}X\right)}^{-1}{X}^T\end{array}\right]y=S_{GNNWR}y$$

(10)

3.2. ATPK for Downscaling the Regression Residuals

The use of spatial interpolation methods to reconstruct high-resolution regression residuals can improve the accuracy of the downscaled model. ATPK is a geostatistical methodology that enables predictions on finer spatial supports than those of the original data [41]. Unlike traditional Kriging interpolation, ATPK considers the size of the spatial support, spatial autocorrelation and point spread function (PSF), with each observation point as the centroid [55].

Based on n low-resolution pixels, an estimation of high-resolution values at point x can be made through ATPK. The mathematical expression can be defined as follows:

$$\widehat{\epsilon }\left(x_0\right)=\sum _{i=1}^n{\lambda }_i\epsilon \left(u_i\right)$$

(11)

In this case, ${\lambda }_i$ represents the weights of each low-resolution pixel, $\epsilon \left(u_i\right)$ stands for the low-resolution pixel values. ATPK estimates the weights through optimal unbiased estimation, where the optimal estimation condition is the minimum error variance, and the unbiased estimation condition is that the sum of the weights equals 1. The kriging system mathematical formula for ATPK is as follows:

$$\left\{\begin{array}{l}{\displaystyle \sum _{j=1}^k}{\lambda }_{i}Cov\left(u_i,u_j\right)+{\mu }_s=Cov\left(u_i,x\right)\\ {\displaystyle \sum _{j=1}^k}{\lambda }_j=1\end{array}\right.$$

(12)

In this formula, $Cov\left(u_i,u_j\right)$ denotes that the covariance between low-resolution pixels $u_i$ and $u_j$, and $Cov\left(u_i,x\right)$ is the covariance between the low-resolution pixel $u_i$ and the target point x. ${\mu }_s$ is the Lagrange multiplier introduced at point x, which aims to constrain the weights. Assuming that the low-resolution pixel value is the mean value of the high-resolution pixel values within the spatial range of that pixel, then the PSF is as follows:

$$h_V\left(x\right)=\left\{\begin{array}{l}\frac{1}{S_v},x\in V\left(x\right)\\ 0,otherwise\end{array}\right.$$

(13)

In which, $S_v$ represents the size of pixel V, and $V\left(x\right)$ denotes the spatial support of pixel V centered around x. Therefore, the area-to-area covariance and area-to-point covariance can be calculated with the definition of PSF:

$$Cov\left(u_i,u_j\right)=\frac{1}{N\left(u_i\right)N\left(u_j\right)}\sum _{k=1}^{N\left(u_i\right)}\sum _{l=1}^{N\left(u_j\right)}Cov\left(v_k,v_l\right),v_k\in u_i,v_l\in u_j$$

(14)

$$Cov\left(u_i,x\right)=\frac{1}{N\left(u_i\right)}\sum _{k=1}^{N\left(u_i\right)}Cov\left(v_k,x\right),v_k\in u_i$$

(15)

In the above formula, $N\left(u_i\right)$ is the number of high-resolution pixels within the spatial range of the original low-resolution pixel $u_i$,$Cov\left(v_k,v_l\right)$ is the covariance between any high-resolution pixel $v_k$ within the low-resolution pixel $u_i$ and any high-resolution pixel $v_l$ within the low-resolution pixel $u_j$; $Cov\left(u_i,x\right)$ represents the covariance between any high-resolution pixel corresponding to the target point x and any high-resolution pixel $v_k$ within the low-resolution pixel $u_i$.

The kriging system of ATPK in matrix form is defined as follows:

$$\left[\begin{array}{cccc}Cov\left(u_1,u_1\right)& \cdots & Cov\left(u_1,u_k\right)& 1\\ ⋮& \ddots & ⋮& ⋮\\ Cov\left(u_k,u_1\right)& \cdots & Cov\left(u_k,u_k\right)& 1\\ 1& \cdots & 1& 0\end{array}\right]\left[\begin{array}{c}{\lambda }_1\\ ⋮\\ {\lambda }_k\\ 0\end{array}\right]=\left[\begin{array}{c}Cov\left(u_1,x\right)\\ ⋮\\ Cov\left(u_k,x\right)\\ 1\end{array}\right]$$

(16)

The key to calculating the kriging weights in ATPK lies in constructing the point-to-point covariance model. The covariance function is a similarity function, which decreases as the distance increases; the semi-variogram function is a dissimilarity function, the variance of the difference in the semi-variogram function will increase as the distance increases, and their sum is constant. Thus, the estimation of the covariance matrix can be converted into the estimation of the semi-variogram function. The point-to-point covariance model at high resolution can be calculated by the deconvolution of the semi-variogram function. An empirical method for deconvolution calculation has been proposed by Wang et al. [55].

3.3. GNNWRK

This article considered the spatial autocorrelation, spatial non-stationarity, and complex non-linearity in spatial relationship modeling and proposed a GNNWRK model suitable for spatial downscaling. The model consisted of the following two parts: regression relationship modeling based on the GNNWR, and residual downscaling based on the ATPK. Its formula can be expressed as:

$${\widehat{Z}}_{GNNWRK}\left(x_0\right)={\widehat{Z}}_{GNNWR}\left(x_0\right)+{\widehat{r}}_{ATPK}\left(x_0\right)$$

(17)

In the formula, ${\widehat{Z}}_{GNNWR}\left(x_0\right)$ is the estimated value at $x_0$ based on the GNNWR, and ${\widehat{r}}_{ATPK}\left(x_0\right)$ is the estimated value of the residual at $x_0$ obtained through ATPK.

In this study, the following six auxiliary factors were used: NDVI, NDBI, WET, NMDI, OSM road, and LST_t1. Hence, the GNNWR model can be expressed as:

$$\begin{gathered} LST=w_0\left(u_i,v_i\right)\times {\beta }_0\left(u_i,v_i\right)+w_1\left(u_i,v_i\right)\times {\beta }_1\times NDVI\left(u_i,v_i\right) \\ \hspace{1em}\hspace{1em}\hspace{1em}+w_2\left(u_i,v_i\right){\beta }_{2}NDBI\left(u_i,v_i\right)+w_3\left(u_i,v_i\right){\beta }_{3}WET\left(u_i,v_i\right) \\ \hspace{1em}\hspace{1em}\hspace{1em}+w_4\left(u_i,v_i\right){\beta }_{4}NMDI\left(u_i,v_i\right)+w_5\left(u_i,v_i\right){\beta }_{5}road\left(u_i,v_i\right) \\ +w_6\left(u_i,v_i\right){\beta }_{6}LST\_t1\left(u_i,v_i\right)+{\epsilon }_i\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \end{gathered}$$

(18)

In the equation, $\left(u_i,v_i\right)$ are the spatial coordinates of target point $\mathrm{i}$, $\beta =({\beta }_1,{\beta }_2,\cdots ,{\beta }_7)$ are the regression coefficients of each factor in the OLR model, and $w=(w_1,w_2,\cdots ,w_7)$ are the spatial non-stationary weights of each factor at point $\mathrm{i}$, solved by SWNN. ${\epsilon }_i$ is the regression residual of the GNNWR model at point i. By dividing the model into the regression term and the residual term, the low-resolution LST (coarse scale LST) can be expressed as:

$${LST}_c=f_{GNNWR}\left(x_c\right)+{\epsilon }_c$$

(19)

where $f_{GNNWR}\left(x_c\right)$ represents the low-resolution regression trend term, and ${\epsilon }_c$ is the low-resolution residual term. Using ATPK for interpolating the residuals, we obtained the high-resolution residuals. Applying the regression relationship $f_{GNNWR}$ to a high spatial resolution, the expression for the high-resolution LST (fine scale LST) ${LST}_f$ is:

$${LST}_f=f_{GNNWR}\left(x_f\right)+f_{ATPK}\left({\epsilon }_f\right)$$

(20)

The GNNWRK model for downscaling NLST is defined as shown in Figure 2.

3.4. GNNWRK Modelling and Implementation

The GNNWRK model replaces the kernel function of classic GWR with SWNN. The structure of SWNN in this study and its implementation strategy are shown in Figure 3.

The input to the SWNN consisted of spatial distance vectors that represented the distances between the target point and all available sample points, and the output layer was the spatial weight matrix of each independent variable [39]. SWNN used Dropout technology [56] to enhance the model’s generalization ability. To improve the efficiency of model iterative optimization, the activation functions of each hidden layer are rectified linear units (ReLU), and the the He parameter initialization method [57] was used. The SWNN adopted the Early Stopping strategy to detect and avoid the overfitting problem of neural network training. In addition, to accelerate the convergence speed of the model and to alleviate the “ Gradient Vanishing “ problem of deep neural networks, Batch Normalization (BN) [58] is also introduced. The parameters settings of the SWNN are shown in

Table 2. Parameter settings of SWNN.

Input Layer	Hidden Layer1	Hidden Layer2	Hidden Layer3	Hidden Layer4	Hidden Layer5	Output Layer
flexible	2048	1024	512	256	128	7
Maximum epoch		Learning rate	dropout	Batch maximum		optimizer
1024		0.0008	0.01	256		Adamax

, where the input shape of the input layer is determined by the size of the study area.

Figure 4 shows the GNNWRK’s implementation process. To enhance the robustness of the model, the dataset was randomly divided into three subsets with ratios of 8:1:1, respectively. The training sets serve the primary purpose of updating and refining the model parameters. Then, the validation set was used as a diagnostic tool to identify any instances of overfitting within the model. The test set was employed to assess the model’s efficacy in fitting data as well as its capacity for generalization. Mean Square Error (MSE) quantified the average disparity between the fitted values and the actual values, which is extensively employed as a metric to assess the accuracy and fitting capabilities of regression models. The study used MSE as the loss function and the MSE of the validation set was utilized as an indicator of overfitting. If the overfitting indicator stopped decreasing or showed an upward trend, and it would have exceeded the pre-set early stopping patient threshold, then it can be considered that the model has begun to overfit, and the training stops.

3.5. Validation and Evaluation

3.5.1. Method Validation

When validating the LST downscaling methods, it was important to consider the systematic errors associated with different thermal infrared sensors. The “upscaling–downscaling” validation strategy has frequently been adopted in previous research endeavors [26,31,49,59]. This strategy is able to validate downscaling results without in-site measurements of LST or alternative high-resolution LSTs thus avoiding the impact of systematic errors among multiple sources of LST on method validation. The “upscaling–downscaling” strategy uses high-resolution LST to simulate the low-resolution LST. In this study, the 100 m ASTER NLST was aggregated to 1000 m as a simulated MODIS NLST, then the simulated MODIS NLST was downscaled to 100 m and validated with the original 100 m ASTER NLST.

3.5.2. Method Evaluation

The study employed four common metrics as evaluation indices for method accuracy, including the following: Pearson’s correlation coefficient (Pcc), root mean square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE). In terms of method comparison, random forest (RF), TsHARP, GWR, and GNNWR were selected. Except for TsHARP, which employed the NDVI as the sole factor in its regression, the other methods employed identical auxiliary factors as those used in GNNWRK. The TsHARP, RF, and GNNWR methods conducted residual correction with coarse model residual, while GWR interpolated the residual with ordinary kriging before residual correction [28,29,33,39].

$$Pcc\left(y,\widehat{y}\right)=\frac{cov\left(y,\widehat{y}\right)}{{\sigma }_y{\sigma }_{\widehat{y}}}$$

(21)

$$RMSE=\sqrt{\frac{{\sum }_{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2}{n}}$$

(22)

$$MAE=\frac{{\sum }_{i=1}^n\left|y_i-\widehat{y_i}\right|}{n}$$

(23)

$$MAPE=\frac{1}{n}\sum _{i=1}^n\left|\frac{y_i-\widehat{y_i}}{y_i}\right|\times 100\%$$

(24)

In the above formula, $y_i$ represents the observed value, ${\widehat{y}}_i$ represents the predicted value, and ${\overline{y}}_i$ represents the mean of the observed values.

4. Results

4.1. LST Downscaling Evaluation

This study validated the downscaling results based on ASTER nighttime LST data, using four performance evaluation indices, Pcc, RMSE, MAPE, MAE, to quantitatively assess the model’s effectiveness. Moreover, it compared the performance of five downscaling methods.

As shown in

Table 3. Quantitative Results of different downscaling models.

Area Code	Model	Metrics
		Pcc	RMSE	MAE	MAPE
A-(1)	GNNWRK	0.922	1.707	1.503	0.506
	GNNWR	0.909	1.809	1.591	0.535
	GWR	0.885	1.987	1.749	0.589
	RF	0.830	2.125	1.798	0.604
	TsHARP	0.870	1.912	1.619	0.544
A-(2)	GNNWRK	0.930	0.886	0.652	0.223
	GNNWR	0.921	0.938	0.691	0.236
	GWR	0.895	1.084	0.799	0.273
	RF	0.843	1.301	0.977	0.334
	TsHARP	0.781	1.605	1.148	0.392
B-(1)	GNNWRK	0.808	0.998	0.738	0.253
	GNNWR	0.786	1.051	0.777	0.267
	GWR	0.748	1.150	0.835	0.287
	RF	0.749	1.115	0.831	0.285
	TsHARP	0.613	1.523	1.121	0.385
B-(2)	GNNWRK	0.763	1.098	0.843	0.291
	GNNWR	0.746	1.134	0.871	0.301
	GWR	0.723	1.180	0.897	0.310
	RF	0.699	1.219	0.932	0.322
	TsHARP	0.497	1.775	1.318	0.455

, the GNNWRK method performed with the highest accuracy in all four study areas, followed by the GNNWR method, while the GWR method had lower accuracy than GNNWR. Compared to the GNNWR, the GNNWRK’s Pcc in the four study areas increased by 0.013, 0.009, 0.022, 0.017, respectively; and RMSE decreased by 5.64%, 5.54%, 5.31%, 3.17%, respectively. The Pcc of the GNNWRK method in three areas exceeded 0.8, with areas A-(1) and A-(2) exceeding 0.9, and the RMSE was below 2 K in all areas, with areas A-(2) and B-(1) below 1 K. In addition, compared to the GNNWR, the GNNWRK method’s MAE and MAPE decreased on average by 4.85% and 4.87%, respectively, with the lowest MAE reaching 0.652 in all areas. Overall, reconstructing high-resolution residuals based on ATPK helped to improve method accuracy. Figure 5 provides a detailed presentation of the accuracy evaluation indices of the downscaling methods in the four study areas.

As shown in Figure 5, the GNNWRK performed well in all four study areas, with all Pcc values greater than 0.7, as well as lower error indices. Comparatively speaking, methods based on geographically weighted theory showed significant improvement over random forest in areas A-(1) and A-(2), while in area B-(1), the performance of GWR was slightly inferior to that of the random forest method. Also, the accuracy of GNNWRK in area B was lower than in area A. Reviewing the characteristics of the study areas, it is known that the landform in area A was more complex than in B, suggesting that the spatial heterogeneity of the area itself has a significant impact on the performance of models that consider spatial non-stationarity. Furthermore, indices including RMSE, MAE, and MAPE reveal that the GNNWRK method achieved the lowest error indices compared to GNNWR and GWR, indicating that its downscaling results provided LST estimates that were closest to the ASTER LST.

Overall, the GNNWRK downscaling method outperformed the benchmark methods, and its performance was less disturbed by varying landform conditions and climate types, demonstrating superior generality and resistance to disruption.

4.2. Spatial Distribution of High-Resolution NLST

Figure 6 shows the overall comparison of NLST spatial distribution across the four study areas. To further analyze the ability of GNNWRK to the mapping of the high-resolution NLST, this study used the ASTER LST product AST_08 as a benchmark, comparing the visual effects of the GNNWRK’s downscaling results with those of other methods.

The overall comparison of NLST spatial distribution suggests that the downscaled NLSTs of the GNNWRK method were closest to the NLST retrieved by the ASTER thermal sensor, with its spatial pattern and NLST differences aligning most closely with the AST_08 product. This suggests that the GNNWRK has better constructed and expressed the spatial relationship between NLST and auxiliary factors. Compared to other downscaling methods, the spatial texture of the GNNWRK downscaling results was more accurate and detailed, clearly depicting the boundaries between the high-temperature and low-temperature regions. Among the benchmark methods, the GNNWR’s details were similar to those of the GNNWRK, but it was less accurate in terms of NLST spatial distribution; the GWR showed some areas being overly smoothed, while overestimating the NLST of roads and buildings, leading to incorrect spatial texture details; the RF exhibited NLST differences that were too subtle, along with more fragmented spatial textures; whereas the TsHARP was significantly affected by “boxy artifacts”, indicating that using regression models with stronger interpretive capabilities and spatial interpolation of residuals can significantly improve the visual effects in the downscaling process over larger spatial resolution spans.

In summary, the GNNWRK downscaling results better captured the spatial distribution characteristics of NLST. For instance, in study area A-(2), the GNNWRK model effectively highlights the high-temperature phenomenon of large impervious surfaces such as parking lots and airports, whereas the GNNWR does not reflect these high-temperature zones as well. In study area B-(1), both the GNNWR and GWR models overstated the central UHIs, while the RF underestimates the extent of these high-temperature regions. This indicates that the GNNWRK accurately reflected the spatial autocorrelation of LST and could capture clusters of high and low temperatures, which was beneficial for identifying urban heat islands and urban cool islands.

The GNNWRK downscaling results show that in study area A, the high LSTs were predominantly cluster within the urban area and extended to the adjacent mountains and dunes, with clear temperature differences between bare soil surfaces and farmland or water bodies, reflecting the thermophysical properties of desert regions. In study area B, the LSTs of the urban confines were noticeably higher than those of the villages and farmlands, demonstrating the urban heat island effect. Unlike area A, the water bodies in area B belong to high-temperature regions, with temperatures even exceeding some urban central areas, highlighting the impact of climate and landcover conditions on the spatial pattern of NLST.

5. Discussion

5.1. The Analysis of the Model Superiority in the Detail Description of the Local Region

In order to visually compare the superiority of the GNNWRK more intuitively, we made detailed comparisons of the downscaled NLST of each study area. Figure 7 shows that the results of each method are quite similar to the spatial distribution pattern of NLST of AST_08, and the GNNWRK has the most reasonable detail spatial textures among all the downscaling methods. GNNWRK performed well in expressing impervious surfaces, correctly constructing the LSTs of the roads and building areas in A-(2) and B-(1). There were significant thermal property differences among the mountains, ridges, and valleys in A-(1) and A-(2), and the GNNWRK downscaled result reflected this difference well. While the GNNWR method showed commendable downscaling capabilities, it was somewhat limited in expressing the spatial distribution of NLSTs. It may exhibit significant deviations in NLST estimates, excessive smoothing, or difficulty in delineating the contours of land surface features. Specifically, the GNNWR underestimated the high-temperature areas present in the mountainous regions of study area A and the water bodies of study area B, while the boundary between impervious surfaces and vegetation in study area B appeared overly smooth. This highlights the necessity of considering spatial autocorrelation within the GNNWR. Urban areas have the urban heat island effect at night, and this phenomenon was captured in the GNNWRK downscaled results in study areas A-(2) and B-(1). Compared with A-(2), the urban heat island effect in area B-(1) was weaker, and the high NLST area did not show an obvious road network structure. This may be due to the fact that the climate in area B-(1) was more humid, which made the impervious surfaces cool down faster at night, and the higher vegetation coverage also inhibited the nighttime urban heat island effect.

Through comparative study on the characteristics of the spatial distribution of NLST, it was found that the GNNWRK has higher fitting accuracy for NLST, and its spatial texture was richer and more detailed compared to other methods. It better explained the spatial autocorrelation and spatial heterogeneity of NLST. Compared with GNNWR, GNNWRK less often exhibited the phenomenon of excessive smoothing and more reasonably captured areas of high and low LSTs. In summary, GNNWRK has the capability to map high-precision NLST in heterogeneous areas, and the GNNWRK downscaled result was more accurate and reasonable compared to the benchmark methods.

5.2. Validation of Downscaling Real MODIS NLST

With the simulated NLST generated by ASTER NLST, the study adopted an “upscaling–downscaling” strategy to circumvent the influence of systematic errors across different thermal infrared sensors on the validation of results. It must be acknowledged that the performance of the downscaling method on simulated NLST may differ from the actual MODIS NLST. Therefore, it was necessary to validate the method’s ability to downscale real MODIS NLST. Specifically, this study used the GNNWRK to downscale the MOD11A1 NLST product to 100 m and used ASTER NLST to validate the downscaling results.

The results shown in Figure 8 indicate that GNNWRK also performed well in downscaling real MODIS NLSTs. Specifically, the downscaling NLST obtained by GNNWRK had a Pcc ranging from 0.748 to 0.925 with ASTER NLST, and their RMSE is less than 1.755 K. In addition, Figure 9 shows that the downscaled NLST well retains the spatial distribution characteristics of MOD11A1 NLST, while also achieving rich spatial details and textures. Notably, GNNWRK still performed well in complex urban areas (A-(2) and B-(1)), proving its potential in obtaining high-resolution NLST in urban regions.

To analyze the differences between the GNNWRK downscaled NLST and ASTER NLST, we calculated the absolute errors between the downscaled NLST and ASTER NLST, as shown in Figure 10. It was found that the errors in study area A-(1) were significantly larger than those in other areas, illustrating the RMSE that was notably greater than that of the other areas. The mountain ridges and valleys of study area A showed larger errors, indicating that the elevation and terrain had a certain impact on the accuracy of downscaling method. Additionally, more errors occurred around the edges of the water body and land in study area B, suggesting that the current auxiliary factors did not adequately account for the contours of water bodies, and there was a need to consider introducing more appropriate auxiliary factors.

Overall, these results demonstrate that GNNWRK has a satisfactory performance in downscaling real MODIS NLST across diverse and complex areas.

5.3. Validation with In-Site NLST

To further evaluate the accuracy of the GNNWRK method in estimating high-resolution NLST, this study compared the NLST ground measured at the USCRN climate observation stations in study area A with the 100 m resolution NLST downscaled by GNNWRK and included the bilinear interpolated MODIS NLST as resampled MODIS for comparison. The USCRN station network provided hourly LST monitoring data, where the hourly average LST was the corrected average of LST measured within an hour; thus, this study selected monitoring data covering the satellite overpass time as the ground validation benchmark. For ease of comparison, the satellite overpass times mentioned are in UTC time.

Table 4. In-site validation of GNNWRK on MOD11A1 LSTs.

Area Code	Time	In-Site NLST/K	GNNWRK NLST/K	Error/K	Resampled MODIS/K	Error/K
A-(1)	02 July 2018 05:30:00	299.15	299.17	0.02	299.42	0.27
A-(2)	16 October 2023 05:01:35	298.55	299.19	0.64	299.27	0.72
A-(3)	08 May 2023 05:05:57	295.65	295.55	0.1	295.17	0.48
A-(4)	09 December 2018 05:28:47	281.15	282.48	1.33	282.69	1.54
A-(5)	25 May 2019 05:34:56	292.25	287.51	4.74	282.90	9.35
mean				1.366		2.472

reveals that the absolute value of error between the GNNWRK downscaled NLSTs and in-site NLSTs was at a mean of 1.366 K, which is a substantial decrease compared to the 2.472 K from the MODIS resampled NLSTs. This suggests that the GNNWRK is capable of retrieving more accurate NLST data, thereby proving the necessity of NLST downscaling. Comparing the ground measured NLSTs with the downscaled and resampled results shows that the downscaled results are consistent with the original MODIS data, meaning they are often simultaneously lower or higher than the ground measurements. Moreover, the downscaled results are closer to the ground measurements, indicating that by establishing a regression relationship of “LST—auxiliary factors”, the downscaling method captures the thermal physical properties of the land surface, making the estimated NLST closer to the actual NLST. Overall, verification based on ground measurements demonstrates that GNNWRK possesses the capability to map high-precision NLST.

5.4. Spatial Texture and Influencing Factors of NLST

To further analyze the spatial texture of NLST and explore the main factors of NLST distribution, we compared the spatial characteristics of NLST with auxiliary factor NDVI, NDBI, WET, and NMDI.

It can be found from Figure 11 that the NLST had obvious spatial differentiation patterns in all four study areas. In study area A-(1), the NLST at the bottom of the valley was significantly lower than the slopes on both sides, which showed similar pattern with WET. The high-temperature area in study area A-(2) was mainly composed of mountain slopes and urban area, with the characteristics of high in the northwest and low in the southeast. At the same time, the NDBI in this area showed its limitations in distinguishing impervious surfaces from bare soil. Study area B had a higher vegetation coverage compared to study area A, and it could be seen that there was a negative correlation between NDVI and NLST in study area B-(1), reflecting the cooling effect of vegetation.

GeoDetector is a statistical method that can quantitatively measure the contribution rates of influencing factors to the dependent variable with its q-value [60,61,62]. To investigate the influencing factors that led to the spatial distribution of NLSTs, we used the NDVI, NDBI, WET, and NMDI as input into the GeoDetector model. The factor detector results showed the contribution rates (%) of the influencing factors to the spatial heterogeneity of NLSTs.

As shown in Figure 12, the ability of influencing factors to explain the spatial heterogeneity of NLSTs varies significantly across different study areas. In study area A-(1) and A-(2), WET had the largest q-value, indicating that the moisture content of the soil and vegetation was the dominant factor in the spatial distribution of NLSTs. NDBI and NDVI also contributed to the NLSTs, while NMDI showed a small influence with q-value less than 5%. Compared to study area A, study area B-(1) had a higher vegetation coverage and a wetter climate, resulting in a significant increase in NDVI’s q-value. In addition, it was found that NDBI plays a more important role in NLST spatial heterogeneity in areas which have more urban build-up area and impervious surfaces.

We also used the interaction detector to identify the interactive effects of influencing factors on NLST distribution and explore the explanatory ability of associated factors [63]. In this study, two types of interactions were identified. The first was a nonlinear enhancement interaction, where the q-values exceeded the combined sum of the individual factors, indicating a nonlinear enhancement of explanatory ability. Another type was a bivariate enhancement interaction, characterized by q-values that were higher than those of each single factor but lower than their combined sum. It was found from Figure 13 that the interactive effects are dominated by nonlinear enhancement, indicating that when selecting a NLST downscaling model, its ability to fit non-linearity should be taken into consideration. Overall, the interactive effects enhanced the explanatory ability of influencing factors, suggesting that the spatial distribution of NLSTs is the result of a combination of multiple factors.

6. Conclusions

This study has taken the spatial non-stationarity, spatial autocorrelation, and complex non-linearity inherent in spatial relationships into account, proposed the GNNWRK model, which combines GNNWR and ATPK, and applied it to downscaling NLST. This method utilized GNNWR to establish the spatial non-stationarity relationship between LST and land surface auxiliary factors that represent landcover conditions and thermal characteristics, as well as used ATPK to interpolate the regression residuals to high resolution.

The GNNWRK method was evaluated in two regions with different geographical and climatic features and compared with four benchmark methods, including the traditional LST downscaling method TsHARP, the classic machine learning model RF, the traditional geographical model GWR, and GNNWR. Experimental results indicate that GNNWRK has higher accuracy and stability compared to the benchmark methods. First, GNNWRK showed improvement in all metrics compared to benchmark methods, for example, Pcc increased by up to 0.022, and RMSE decreased by up to 5.64%. Secondly, although the accuracy of GNNWRK varied across different study areas, its overall performance remained superior to benchmark methods, which means that GNNWRK has good stability and generality. Additionally, validation based on in-site NLST data confirmed that GNNWRK could be used to obtain high-precision NLST.

Overall, the proposed GNNWRK method helped to overcome the spatiotemporal resolution trade-off of satellite TIR sensors and has provided technical support for obtaining high spatiotemporal resolution NLST.