A Spatiotemporal Fusion Model of Land Surface Temperature Based on Pixel Long Time-Series Regression: Expanding Inputs for Efficient Generation of Robust Fused Results

Abstract

Spatiotemporal fusion technology effectively improves the spatial and temporal resolution of remote sensing data by fusing data from different sources. Based on the strong time-series correlation of pixels at different scales (average Pearson correlation coefficients > 0.95), a new long time-series spatiotemporal fusion model (LOTSFM) is proposed for land surface temperature data. The model is distinguished by the following attributes: it employs an extended input framework to sidestep selection biases and enhance result stability while also integrating Julian Day for estimating sensor difference term variations at each pixel location. From 2013 to 2022, 79 pairs of Landsat8/9 and MODIS images were collected as extended inputs. Multiple rounds of cross-validation were conducted in Beijing, Shanghai, and Guangzhou with an all-round performance assessment (APA), and the average root-mean-square error (RMSE) was 1.60 °C, 2.16 °C and 1.71 °C, respectively, which proved the regional versatility of LOTSFM. The validity of the sensor difference estimation based on Julian days was verified, and the RMSE accuracy significantly improved (p < 0.05). The accuracy and time consumption of five different fusion models were compared, which proved that LOTSFM has stable accuracy performance and a fast fusion process. Therefore, LOTSFM can provide higher spatiotemporal resolution (30 m) land surface temperature research data for the evolution of urban thermal environments and has great application potential in monitoring anthropogenic heat pollution and extreme thermal phenomena.

1. Introduction

The urban thermal environment is one of the most significant phenomena of human activity on Earth [1], affecting various aspects, including water and air quality, microclimate, energy consumption, and human health [2,3,4]. Land surface temperature (LST), as long-term, wide-area data easily accessible by satellites, is the most common quantitative indicator used in climate and environmental studies [5,6,7,8]. Spatial resolutions greater than 50 m [9], represented by the Landsat series of satellites, are considered suitable for detailed studies of the urban thermal environment because they reflect the thermal structure of internal urban objects, such as streets and factories [10,11,12,13]. The actual availability of LST data from Landsat satellites is limited in the time series because of the 16-day revisit cycle, cloud cover, and sensor failures [14,15,16]. Whereas moderate resolution imaging spectroradiometer (MODIS) data, for example, can meet the temporal resolution requirements of LST data for many thermal studies [17,18], although it is difficult to adapt the spatial resolution to the needs of increasingly fine-scale urban studies.

Spatiotemporal fusion received widespread attention for producing higher-resolution data, which existing single-sensors cannot provide. The spatial and temporal adaptive reflectance fusion model (STARFM) [19] was one of the first spatiotemporal fusion methods to become popular and spawned other weight function-based spatiotemporal fusion methods [20,21,22], such as the enhanced STARFM (ESTARFM) [23]. The division of spatiotemporal fusion methods based on differences in principles [24,25] include unmixing-based [26,27], weight function-based [19], Bayesian-based [28,29], learning-based [30,31,32], and hybrid methods [33,34]. Among them, hybrid methods synthesize the above different types of methods to achieve complementary advantages, and flexible spatiotemporal data fusion (FSDAF) [35] is representative and derives many improved models [36,37,38].

The spatiotemporal fusion models were initially applied for generating reflectance products. Based on the need for the fusion of LST data, classical reflectivity spatiotemporal fusion models such as STARFM, ESTARFM, and FSDAF were attempted subsequently [39,40]. Additionally, in response to the spatiotemporal distribution characteristics of LST data, several spatiotemporal fusion methods specifically for LST were designed and proposed, including spatiotemporal image fusion models based on bilateral filtering [41], spatiotemporal adaptive data fusion algorithms for temperature mapping (SADFAT) [42], spatiotemporally integrated temperature fusion models (STITFM) [43], and some other methods [44,45,46].

To minimize changes in land cover, the image closest to the predicted date was selected as the reference image [47,48,49]. However, LST is also characterized by significant seasonal and weather variability [40], and even temporally close LST images can exhibit large numerical differences [50,51]. In this case, the effect of different single-image pair (minimum) inputs on the fusion results can be significant [52]. Some spatiotemporal fusion models enhance the fusion performance by adding a reference image to the minimum input [23,53]. Furthermore, the extended input provides a new perspective for reducing the uncertainty of single-image pairs.

In contrast to minimum input, extended input involves feeding as much data as possible while ensuring data quality. Based on the synthesis of available data, an artificial selection process was avoided. As illustrated in Figure 1, the implementation framework for extended input can be categorized into the following: (a) Weighting the fusion results obtained from different image pairs. (b) Combination of multiple image pairs to derive an image pair that is closer to the prediction time for fusion. (c) Fusion based on the mapping relationships derived from multiple image pairs. Among these, the (a) framework is the earliest to emerge [19]. However, as the number of image pairs increases, the time consumption grows linearly. Therefore, even though methods like STARFM and spatial-temporal data fusion approach (STDFA) [26] incorporate such extended input designs, in practical usage, single image pairs are often preferred to mitigate the time overhead. Compared to the (a) framework, the (b) combines inputs to avoid repetitive fusion processes. Spatiotemporal models such as virtual image pair-based spatiotemporal fusion (VIPSTF) [54] and robust optimization-based fusion (ROBOT) [55] have adopted this strategy. The (c) framework enhances model performance by delving into the mapping relationships among multiple image pairs and is one of the important trends in the development of extended input fusion [56]. However, the mapping relationships among multiple image pairs have not been exhaustively explored yet.

In comparison to traditional spatiotemporal fusion using single-reference image pairs, we endeavored to enhance robustness and fusion accuracy by expanding the model’s input. Specifically, our approach was based on mapping relationships between coarse and fine pixels over a long time series (2013–2022) and thus referred to as the long time-series spatiotemporal fusion model (LOTSFM). We proved the reasonableness of LOTSFM in principle, and three cities were experimentally verified. The results showed that fusion based on pixel-long time-series regression was feasible and exhibited high fusion accuracy and speed. In addition, relative to previous studies that treated the sensor difference term as a constant, we considered its variation in the spatiotemporal fusion of LST and introduced the variable Julian day for modeling, which optimized fusion accuracy.

The key content arrangement for Section 2 is as follows: In Section 2.3, the materialized extended input (c) framework is shown. Section 2.4 and Section 2.5 show that the linear and non-linear mapping parts of the fusion were modeled separately.

Table 1. Abbreviations and definitions.

Acronym	Definition
AD	Average difference
AAD	Average absolute difference
APA	All-round performance assessment
GAN	Generative adversarial network
GAN-STFM	GAN-based spatiotemporal fusion model
LBP	Local binary patterns
LOOCV	Leave-one-out cross-validation
LST	Land surface temperature
LOTSFM	Long time-series spatiotemporal fusion model
LTRM	Long time-series regression model
MODIS	Moderate resolution imaging spectroradiometer
MPF	Missing pixel filling
OLS	Ordinary least squares
PCCs	Pearson correlation coefficients
RC	Residual compensation
RMSE	Root-mean-square error
SDE	Sensor difference term estimation
SF	Spatial filtering
STARFM	Spatial and temporal adaptive reflectance fusion model
Coarse image	Image with relatively low spatial resolution;
Coarse pixel	Pixel in the coarse image;
Fine image	Image with relatively high spatial resolution;
Fine pixel	Pixel in the fine image;
Image pair	Coarse and fine image of the same location on the same date;
$C\left(x_i,t\right)$	$\mathrm{Value}\mathrm{of}\mathrm{coarse}\mathrm{pixel}{x}_i$$\mathrm{at}t$
$F\left(x_{ij},t\right)$	$\mathrm{Value}\mathrm{of}\mathrm{fine}\mathrm{pixel}{x}_{ij}$$\mathrm{at}t$
$C_F\left(x_i,t\right)$	$\mathrm{Value}\mathrm{of}\mathrm{ideal}\mathrm{coarse}\mathrm{pixel}{x}_i$$\mathrm{at}t,$ defined as Equation (6)
$\xi (x_i,t)$	$\mathrm{Value}\mathrm{of}\mathrm{sensor}\mathrm{difference}\mathrm{at}\mathrm{the}\mathrm{position}\mathrm{of}\mathrm{coarse}\mathrm{pixel}{x}_i$$\mathrm{at}t$
$a,b$	Time-series pixel regression parameters
${\alpha }_{\xi }, {\beta }_{\xi }{, {\gamma }_{\xi }, \theta }_{\xi }$	Time-series fitting parameters for sensor difference terms
$r$, $r_{CF}$	Fit residuals

provides a summary of the abbreviations and definitions used in this study.

2. Materials and Methods

2.1. Fusion Data Types

A variety of temperature inversion algorithms can be used to generate LST data [57,58,59]. One straightforward way is to use processed LST products. LST product data from MODIS and Landsat after multiple experiments and iterative evaluations are highly accurate and reliable, providing standardized LST results compared with different inversion algorithms [60,61,62] and are therefore used in many LST fusion studies [42,43,63]. In this study, LST products, including Landsat8/9 Level-2, MOD11A1, and MYD11A1, were also used as inputs to the model after band math processing.

Landsat 8/9 OLI_TIRS Collection 2 Level-2 product data covering the study area were obtained from the United States Geological Survey (USGS, http://glovis.usgs.gov/, accessed on 6 February 2023). Their LST products are recommended by the National Aeronautics and Space Administration (NASA) to avoid atmospheric parameter calculations. For the Level-2 product, Equation (1) calculation is performed to obtain the LST in degrees Celsius (°C), which is the default temperature unit in this study, and band 10 is the relevant LST band.

$${LST}_L=0.00341802\ast band10+149-273.15$$

(1)

This Level-2 product of Landsat 8/9 corresponds to LST data at a spatial resolution of 30 m. MODIS data were sampled to this spatial resolution after nearest-neighbor interpolation.

Daily MOD11A1 and MYD11A1 data covering the study area were obtained from NASA (http://ladsweb.nascom.nasa.gov/, accessed on 8 February 2023). MODIS images corresponding to Landsat dates were collected, projected, cropped, and resampled. Band 1, LST_Day_1km, was calculated by Equation (2) to obtain the LST.

$${LST}_M=0.02\ast band1-273.15$$

(2)

The acquired MODIS LST data were combined with Landsat LST data to form image-pair data for each date.

2.2. Regions and Date Range

The cities of Beijing, Shanghai, and Guangzhou are more evenly distributed across different latitudes in China, with the warm temperate semi-humid continental monsoon, northern subtropical maritime monsoon, and southern subtropical maritime monsoon climates, all of which are central cities in China and are important regions for the study of urban thermal environments [64,65].

LOTSFM is expected to be used to generate spatiotemporal data for LST that can be used for urban thermal environment studies; therefore, these important cities were selected for fusion as representatives. The locations of the three experimental regions, Beijing, Shanghai, and Guangzhou, are shown in Figure 2. Among them, the Guangzhou City region was selected, including a portion of the adjacent Foshan City, to avoid the images being affected by MODIS strips. The three square regions are of the same size, 36 km on one side, approximately 1300 km², and include 1.44 million Landsat 8/9 fine pixels.

In this study, all the available image pairs for each experimental region from 2013 to 2022 were collected (with as many as possible) to satisfy the extended input and validation of the experiments using the following process. (1) Cloud-free Landsat images of the study areas were screened based on cloud cover and quality bands. (2) MODIS images corresponding to cloud-free Landsat images were collected to form pairs. (3) Image pairs with large areas of missing pixels were filtered out.

Among the three experimental regions, Beijing used completely missing pixel-free samples for fusion because of the abundance of such samples. The percentage of missing pixels in the Shanghai and Guangzhou MODIS samples was maintained below 10% as far as possible to ensure the quality of filled and fused data. The final data filtered for input into the LOTSFM model are summarized in

Table 2. Summary of experimental data dates by regions.

Beijing		Shanghai		Guangzhou
Image ID	Date	Image ID	Date	Image ID	Date
MOD1-L1	2014.09.04	MOD1-L1	2013.05.25	MYD1-L1	2013.11.29
MOD2-L2	2014.10.06	MOD2-L2	2013.07.12	MYD2-L2	2013.12.31
MOD3-L3	2015.04.16	MOD3-L3	2013.08.13	MYD3-L3	2014.01.16
MOD4-L4	2015.05.18	MOD4-L4	2013.11.17	MYD4-L4	2014.10.15
MOD5-L5	2017.05.07	MOD5-L5	2013.12.03	MYD5-L5	2014.11.16
MOD6-L6	2017.07.10	MOD6-L6	2014.11.04	MYD6-L6	2015.01.03
MOD7-L7	2017.09.12	MOD7-L7	2015.03.12	MYD7-L7	2015.01.19
MOD8-L8	2017.09.28	MOD8-L8	2015.08.03	MYD8-L8	2015.10.18
MOD9-L9	2017.10.30	MOD9-L9	2016.01.26	MYD9-L9	2016.02.07
MOD10-L10	2018.04.08	MOD10-L10	2017.02.13	MYD10-L10	2016.12.07
MOD11-L11	2018.10.01	MOD11-L11	2017.04.02	MYD11-L11	2017.10.23
MOD12-L12	2018.10.17	MOD12-L12	2017.08.24	MYD12-L12	2019.09.27
MOD13-L13	2019.05.29	MOD13-L13	2018.03.04	MYD13-L13	2019.10.29
MOD14-L14	2019.06.14	MOD14-L14	2018.05.23	MYD14-L14	2019.11.14
MOD15-L15	2019.09.02	MOD15-L15	2018.12.17	MYD15-L15	2020.02.18
MOD16-L16	2019.09.18	MOD16-L16	2019.01.18	MYD16-L16	2021.01.19
MOD17-L17	2020.04.13	MOD17-L17	2019.07.29	MYD17-L17	2021.02.04
MOD18-L18	2020.08.03	MOD18-L18	2019.12.04	MYD18-L18	2021.02.20
MOD19-L19	2020.09.20	MOD19-L19	2020.01.21	MYD19-L19	2021.12.05
MOD20-L20	2021.05.02	MOD20-L20	2020.02.22	MYD20-L20	2022.09.03
MOD21-L21	2021.06.03	MOD21-L21	2020.05.12	MYD21-L21	2022.10.21
MOD22-L22	2021.06.19	MOD22-L22	2020.08.16
MOD23-L23	2021.09.07	MOD23-L23	2021.04.29
MOD24-L24	2022.01.29	MOD24-L24	2022.01.02
MOD25-L25	2022.03.02	MOD25-L25	2022.02.27
MOD26-L26	2022.03.26	MOD26-L26	2022.03.15
MOD27-L27	2022.04.19	MOD27-L27	2022.03.23
MOD28-L28	2022.05.13	MOD28-L28	2022.04.08
MOD29-L29	2022.05.21	MOD29-L29	2022.09.07

, with 29, 29, and 21 groups in Beijing, Shanghai, and Guangzhou, respectively, for a total of 79 image pairs of data.

In Table 2, L, MOD and MYD are Landsat 8/9, MOD11A1, and MYD11A1, respectively. MOD-L/MYD-L are image pairs. Landsat 9 is generally regarded as a replicated version of Landsat 8, with essentially the same sensor conditions and observation times, so the two data are considered homologous in this study and are mixed to increase the length of the time series. Although MOD11A1 and MYD11A1 are both MODIS LST data, they are respectively generated by the Terra and Aqua satellites. We specifically utilized their daytime land surface temperature products for our study. However, compared to Landsat 8 and 9, there are differences in the timing of daytime overpasses for MOD and MYD. Therefore, we did not mix their data in our research but instead used them separately. A reference table of approximate satellite crossing times for each data in this study can be seen in the Supplementary Materials.

2.3. Fusion Model Structure

LOTSFM is structured with reference to Fit-FC, a spatiotemporal fusion model that explicitly decomposes the structure into regression model fitting (RM), spatial filtering (SF), and residual compensation (RC) [66]. Many spatiotemporal fusion models can identify components similar to those of Fit-FC [67,68,69]. Additionally, the structure of the LOTSFM was adjusted to address the problem of long-time-series fusion for LST. Sensor difference term estimation (SDE) replaces the RC to reduce errors due to significant spatial and temporal variations in the difference term. An additional missing pixel-filling (MPF) component was added to ensure that as many time series as possible were considered for fusion, as shown in Figure 3.

The fusion process of LOTSFM is divided into two stages: training and prediction. The training stage was performed only once to obtain a time-series generic mapping relationship between the coarse and fine image pairs. The prediction stage generated fine images on each prediction date based on the mapping relationship. This is similar to the neural network working process; however, LOTSFM is based on the traditional principle that the training stage constructs mapping relations mainly through regression, which is expected to be faster. LOTSFM is based on time-series generic parameters obtained in the training stage, which are used in the fusion prediction for each date without the need for repeated calculations, thereby improving the efficiency of fusion.

2.4. Long Time-Series Regression Model

For the LST model, it is assumed that there is a high correlation between neighboring fine pixels in the time series. This was based on two considerations. (1) For fine pixels that are sufficiently adjacent, their values should correlate according to the first law of geography [70]. (2) Although abrupt changes such as land cover may occur over time and distance, the LST can still be strongly correlated in the time series due to the heat transfer between adjacent objects and similar conditions of solar thermal radiation. Based on this assumption, it is reasonable to establish the equation between any fine pixels in the neighborhood using coarse pixels as the neighborhood range, which can be described as

$$F\left(x_{ik},t\right)=a\left(x_{ik}{,x}_{ij}\right)F\left(x_{ij},t\right)+b\left(x_{ik}{,x}_{ij}\right)+r\left(x_{ik}{,x}_{ij},t\right)$$

(3)

where, $F(x_{ik},t)$ and $F\left(x_{ij},t\right)$ are the LST values of any two fine pixels in the coarse pixel $x_i$ at date$t$, respectively. $a\left(x_{ik}{,x}_{ij}\right)$ and $b\left(x_{ik}{,x}_{ij}\right)$ are coefficients that do not vary with t, $r\left(x_{ik}{,x}_{ij},t\right)$, which is the residual.

According to the pixel-temperature mixing principle [71], the value of the coarse pixel $x_i$, $C(x_i,t)$, is equal to the average of the fine pixel values $F\left(x_{ik},t\right)$ in it plus the $\xi$ term, which can be expressed as:

$$C\left(x_i,t\right)=\frac{1}{m}\sum _{k=1}^{m}F\left(x_{ik},t\right)+\xi \left(x_i,t\right)$$

(4)

where $\xi$ is the bias caused by the observation of different sensors and is referred to as the sensor difference term. Considering that $\xi$ of LST may be significantly different across times and coarse pixel locations, it is set to $\xi (x_i,t)$. $m$ is the number of fine pixels within a coarse pixel. From Equations (1) and (2), the pixel relationship between the different scales can be expressed as

$$C\left(x_i,t\right)=F\left(x_{ij},t\right)·\frac{1}{m}\sum _{k=1}^{m}a\left(x_{ik}{,x}_{ij}\right)+\frac{1}{m}\sum _{k=1}^{m}b\left(x_{ik}{,x}_{ij}\right)+\frac{1}{m}\sum _{k=1}^{m}r\left(x_{ik}{,x}_{ij},t\right)+\xi \left(x_i,t\right).$$

(5)

Similar to Equation (3), the relationship between the pixels at different scales exhibits linearity. The difference is that the sensor difference term $\xi$ in Equation (5) causes the correlation between pixels at different scales to weaken compared to those of the fine pixels at the same scale. We expect to exclude the effect of the difference term $\xi$ to retain a stronger correlation. Therefore, we set:

$$C_F(x_i,t)=C(x_i,t)\mathrm{-}\xi (x_i,t)$$

(6)

where, $C_F\left(x_i,t\right)$ is the coarse pixel value for the ideal case in which no observed difference exists between the sensors.

Based on Equations (5) and (6), the relationship between the ideal coarse pixel $x_i$ and the contained fine pixel $x_{ij}$ at any date $t$ can be expressed as

$$F\left(x_{ij},t\right)={\alpha }_F\left(x_{ij}\right)C_F\left(x_i,t\right)+{\beta }_F\left(x_{ij}\right)+r_{CF}\left(x_{ij},t\right)$$

(7)

where, $\alpha \left(x_{ij}\right)$ and$\beta \left(x_{ij}\right)$ can be estimated using the ordinary least squares method (OLS) [72,73,74]. The residual $r\left(x_{ik}{,x}_{ij},t\right)$ between fine pixels, is assumed to be an unpredictable fraction caused by abrupt spatial and temporal variations in land cover that cannot be reduced by a more macroscopic-scale least-squares fit between coarse and fine pixels. Thus, transforming Equation (5),$r_{CF}\left(x_{ij},t\right)$ for the predicted residuals is expressed as:

$$r_{CF}\left(x_{ij},t\right)=-{\left[\sum _{k=1}^{m}a\left(x_{ik}{,x}_{ij}\right)\right]}^{-1}·\sum _{k=1}^{m}r\left(x_{ik}{,x}_{ij},t\right)$$

(8)

Because $a\left(x_{ik}{,x}_{ij}\right)$ floats around 1, $r_{CF}\left(x_{ij},t\right)$ is approximated as the regional average of $r\left(x_{ik}{,x}_{ij},t\right)$. Thus, the ideal coarse pixel value $C_F$ constructed would have a high correlation with $F$ close to that between the fine pixels and, on average, make the correlation more stable. Therefore, a well-regressed relationship between $C_F$ and $F$ is expected.

In contrast to previous temporally dependent spatiotemporal fusion principles [24], LOTSFM generalizes this dependence principle for different locations and long time series according to Equation (5). Specifically, previous temporal dependence principles generally exist between the reference and predicted dates at the same location, and fusion parameters are derived from linear sections of the temporal profile, such as Fit-FC and ESTARFM. The LOTSFM generalizes the linear form applicable to multiple dates (long time series) at different locations. In Section 3.1, the principle is validated based on the fact that pixels at different scales exhibit strong correlations over a long time series.

2.5. Sensor Difference Term Estimation

The sensor difference term is generally regarded as a constant in surface reflectance spatiotemporal fusion models, including classical models such as STARFM and FSDAF. However, the effect of the sensor difference term on the LST fusion problem was significant across time and location, possibly because of data differences. For example, the significant variation of $\xi$ with time and location for the Beijing region is shown in Figure 4. Therefore, the sensor difference term is defined as Equation (9), which varies with the time and location of coarse pixel acquisition, was considered in this study.

$C\left(x_i,t\right)$ and $J\left(t\right)$ are two variables available for estimating the difference term of the sensor. According to earlier studies [42,75], the introduction of temporal variables has the potential to improve spatiotemporal fusion. $J\left(t\right)$ is the transformation of the date represented by $t$ into a Julian day, with the year discarded to facilitate its introduction into the estimation of the sensor difference term. It can also be regarded as the day of the year (DOY). For example, if $t$ represents 1 August 2017, the corresponding Julian date is 213. We provide the following fit model that can be used for the sensor difference term estimation:

$$\widehat{\xi }\left(x_i,t\right)={\alpha }_{\xi }\left(x_i\right)C\left(x_i,t\right)+{\beta }_{\xi }\left(x_i\right)+{\gamma }_{\xi }\left(x_i\right)sin\left(\frac{2\pi }{365.25}J\left(t\right)+{\theta }_{\xi }\left(x_i\right)\right)$$

(9)

where, ${\alpha }_{\xi }$, ${\beta }_{\xi }$, ${\gamma }_{\xi }$ and ${\theta }_{\xi }$ are fitting parameters that vary for the location $x_i$ of the coarse pixel $C\left(x_i,t\right)$.

The form of the Equation (9) fit is explained as follows: (1) The larger the coarse pixel value $C$, the more likely it is that a larger $\xi$ appears, showing a certain linear correlation. (2) Jules Day $J\left(t\right)$as a date variable: The sinusoidal function is more often used to express periodicity [76]. Furthermore, since the period of the Julian day is fixed, the frequency of the sinusoidal function is known to be $2\pi /365.25$. In Figure 4, estimated values $\widehat{\xi }$ tend to be within the interquartile range (Q3–Q1), demonstrating alignment with the observed variations.

2.6. Missing Pixel Filling

The missing pixel problem limits the use of spatiotemporal fusion models and reduces the spatiotemporal continuity of fused data. Based on the assumptions considered in Section 2.4 for pixel correlation, there should also be a time-series correlation between neighboring coarse pixels. Therefore, a solution to the problem of a limited number of missing pixels occurring in coarse-resolution images on the predicted date can be provided in view of pixel correlation. It is assumed that the level of correlation between coarse pixels is sufficient to fulfill the need for a suitable window range. The correlation between the coarse pixels is described in Section 3.1. The equation for filling any coarse pixels within the search window can be established as follows:

$${\widehat{C}}_{MPF}\left(x_u,t_p\right)=\alpha \left(x_u\right)C\left(x_v,t_p\right)+\beta \left(x_u\right),$$

(10)

$$x_v=\underset{x\ne x_u}{\arg \max }corrcoef(c(x),c(x_u)).$$

(11)

In Equation (10), ${\widehat{C}}_{MPF}\left(x_u,t_p\right)$ is the regression estimate of the missing pixels $C\left(x_u,t_p\right)$, and the coefficients $\alpha \left(x_u\right)$ and $\beta \left(x_u\right)$ can be solved by OLS. Equation (11) indicates that $C\left(x_v,t\right)$, which is used to establish the regression relationship with $C\left(x_u,t\right)$, is the coarse pixel with the highest time-series correlation.

2.7. Spatial Filtering

The primary cause of blocky artifacts in spatiotemporal fusion is the inconsistency in spatial resolution [67,77]. The raster boundaries of coarse-resolution images are not real land cover boundaries, although they are fused into spatial textures at fine resolution, which causes artifacts [78,79]. Spatial filtering by calculating the spectral differences to find similar pixels is a classic solution to artifact problems [19]. However, the efficiency of the pixel-level search approach is relatively low [24,80], and pure pixels and classification principles may not be suitable for continuous LST data [63].

Inspired by the All-around Performance Assessment (APA) framework [81], a fused image is considered to contain spectral and spatial information. The coarse-resolution image on the predicted date provides overall spectral information, whereas the fine-resolution image on the reference date provides detailed spatial information that is not provided by the coarse image. Therefore, we consider spatial filtering of coarse images to weaken the coarse spatial resolution spatial information as much as possible so that more fine-resolution image spatial information is inherited to suppress the appearance of artifacts. Specifically, the coarse resolution image is sampled to fine resolution based on nearest neighbor interpolation. The blocky texture caused by the coarse raster is removed based on mean filtering and then fed to the fuser as a normal coarse image.

3. Experiments and Results

3.1. Time-Series Correlation Levels between Pixels

Correlation statistics for the study regions can corroborate the principal assumption [82] and help anticipate the fusion performance of the LOTSFM. Correlations among coarse pixels, ideal coarse pixels, and fine pixels were assumed to construct LOTSFM. Therefore, we verified the variations in these three types of correlations using the Euclidean distance from the center of the pixel. In each region, 20,000 samples were taken for each type of correlation, and the distances were discretized to form Figure 5.

The average correlation for all types in the regions was greater than 0.95, which is in line with the principal assumption of high correlations between adjacent pixels. The correlations between coarse pixels (CC) were largely stable with distance over the experimental region range; therefore, the search window for missing pixel filling was set to be region global. The correlation between the ideal coarse pixels and fine pixels (CF) and between fine pixels (FF) decreases with distance similarly, which conforms to the expectation of ideal coarse pixels in principle, i.e., to retain the stronger correlation by excluding the difference term $\xi$.

3.2. Regional Replicability Assessment of LOTSFM

The LOTSFM fusion accuracy experiments were conducted in Beijing, Shanghai, and Guangzhou, China. Due to the small dataset size of the fusion, leave-one-out cross-validation (LOOCV) was considered appropriate, enabling as many rounds of validation as possible. LOOCV avoids randomness in the division of training and validation samples so that the accuracy estimates for a given dataset are constant and objective [83]. LOOCV was performed by removing one image pair at a time for validation and feeding the remaining data into the training set.

Accuracy is assessed with reference to the all-round performance assessment (APA) [81], which provides a comprehensive and standard model assessment framework adopted by new spatiotemporal fusion models, such as variation-based spatiotemporal data fusion (VSDF) [69] and the comprehensive flexible spatiotemporal data fusion (CFSDAF) [63]. Specifically, the APA proposes using metrics including the average difference (AD), root-mean-square error (RMSE), Robert’s edge (EDGE), and local binary patterns (LBP) as evaluation criteria for assessing the performance of spatiotemporal fusion algorithms.

Based on the data in Table 2, Beijing, Shanghai, and Guangzhou conducted 29, 29, and 21 rounds of validation, respectively. The results of each round are shown in Tables S1–S3 in the Supplementary Materials, and the statistics are summarized in

Table 3. Summary statistics of LOOCV results for each region of LOTSFM (AD and RMSE unit: °C).

City	Statistic	AD	RMSE	EDGE	LBP
Beijing	Average value	0.006038	1.603804	−0.206396	−0.000472
	Standard deviation	1.090600	0.599686	0.044764	0.002536
Shanghai	Average value	0.052643	2.169908	−0.151468	−0.000344
	Standard deviation	2.104711	1.207505	0.033107	0.004028
Guangzhou	Average value	−0.129306	1.715259	−0.093829	−0.001072
	Standard deviation	1.687730	1.058368	0.020060	0.001607

.

In the three-region experiment, the lowest absolute average AD value is 0.006 °C in Beijing, and the highest is about 0.13 °C in Guangzhou. The small average AD results for the three regions suggest that LOTSFM provides an approximately unbiased estimate for overall long time-series samples. The minimum RMSE is in Beijing, 1.60 °C, and the maximum is in Shanghai, 2.17 °C.

3.3. SDE Component Validity Assessment of LOTSFM

To verify the validity of the introduction of the Julian day for the estimation of the sensor difference term, we also conducted LOOCV for each region of LOTSFM without the Julian day input. The results of each round are shown in Tables S4–S6, and the statistics are summarized in Table S7 in the Supplementary Materials.

To provide a more visual representation of the optimization in each round, Figure 6 shows the results based on paired box plots. The AD of each region showed a decrease in variance; the values were more concentrated at approximately 0. The RMSE exhibited a significant decrease in the average value (p < 0.05). The introduction of Julian days does not include additional spatial information; therefore, the EDGE and LBP indicators, which are related to spatial detail, did not show significant changes.

The three-city fusion images obtained from the cross-validation experiments are shown in Figure 7. Despite the significant difference in resolution between MODIS (1 km) and Landsat (30 m), the LOTSFM algorithm can produce fused images with spatial details similar to those of Landsat using MODIS data from the predicted data. Although Figure 7 displays only thumbnail images, the thermal environment textures within the city, such as rivers, are visible in the Landsat and LOTSFM fusion images. In contrast, MODIS provides coarse information on the distribution of urban temperatures. This highlights the critical need to obtain high-resolution temperature data through spatiotemporal fusion for investigating urban thermal environments.

3.4. Fusion Accuracy Comparison with Other Models

To assess the LOTSFM fusion capability, STARFM, ESTARFM, FSDAF, and GAN-STFM were selected for comparison experiments, all of which are classical and easily accessible from open sources. STARFM, the earliest and most widely used spatiotemporal fusion model, is often used as a standard model for comparison [41,44,84]. ESTARFM is an enhanced version of STARFM, which has been used in many LST fusion studies [40,85]. FSDAF, as a hybrid fusion model, combines the strengths of weight function-based and unmixing-based methods and has achieved relatively good results in some fusion comparison experiments [86,87]. GAN-STFM [56] is a deep-learning fusion method that uses a classical generative adversarial network (GAN). GAN-STFM has a fusion process similar to that of LOTSFM, which involves training and prediction and is therefore used for comparison experiments.

A comparative experiment based on the Beijing region was conducted using four randomly selected image pairs for independent validation. The dates selected were 7 May 2017, 8 April 2018, 2 September 2019, and 29 January 2022. Based on the principle of the nearest date input and the respective model input requirements, STARFM and FSDAF use the nearest one-date image pair, whereas ESTARFM uses the nearest two-date image pairs as the fusion input. The GAN-STFM model allows as many image pairs as possible to be used for training; therefore, LOTSFM uses the rest of the 25 dates image pair as input.

The experimental accuracy results are visualized using the APA diagrams in Figure 8, and the results of the accuracy parameters are listed in

Table 4. Accuracy assessment of five spatiotemporal fusion models based on APA. The optimal values are indicated by asterisks (*) in each column.

Image Date	2017/5/7				2018/4/8
Metric	AD	RMSE	EDGE	LBP	AD	RMSE	EDGE	LBP
STARFM	−0.1530	2.8298	−0.2935	0.0018	−2.9707	3.3658	−0.3530	−0.0166
ESTARFM	−0.1428	2.5451	−0.2182 *	−0.0022	−1.7855	2.5301	−0.2278	−0.0113
FSDAF	0.3246	2.8239	−0.2580	−0.0084	−2.8314	3.3243	−0.4500	−0.0376
GAN-STFM	1.0141	3.0614	−0.4539	−0.0040	−1.7510	2.2203	−0.4331	0.0060
LOTSFM	−0.0547 *	1.7850 *	−0.2580	−0.0017 *	−0.7164 *	1.5016 *	−0.2211 *	−0.0050 *
Image Date	2019/9/2				2022/1/29
Metric	AD	RMSE	EDGE	LBP	AD	RMSE	EDGE	LBP
STARFM	−2.5883	2.8231	−0.2486	−0.0039	1.9178	2.1979	−0.3586	0.0069
ESTARFM	−0.1020 *	0.9795 *	−0.0856 *	−0.0014 *	3.0526	3.6576	−0.1013 *	0.0005
FSDAF	−2.3037	2.5587	−0.1911	−0.0074	2.0815	2.3886	−0.2895	0.0004 *
GAN-STFM	−0.7339	1.3711	−0.2766	0.0051	0.5294	1.3492	−0.5593	−0.0020
LOTSFM	−0.6541	1.1527	−0.1162	0.0014 *	0.2087 *	1.3394 *	−0.2971	0.0040

. In the APA diagram, the closer to the center of the semicircle, the better the RMSE accuracy, and the closer to the central axis, the better the spatial detail, using different colors to indicate the AD. The APA delineated green areas of good and grey areas of fairness based on the input to visually evaluate the fusion performance. The fair area indicates that the fused image is better than the input image in one of the spectral or spatial accuracy aspects, and a good area is better in both aspects. In Figure 8, ESTARFM and LOTSFM are categorized as good in (b) and (c). The contemporaneous optimal accuracy values are listed in Table 4, with 10 optimal values for LOTSFM, 6 for ESTARFM, and 1 for FSDAF. ESTARFM and LOTSFM were outstanding in the experiments; however, LOTSFM was more robust. In (d), the ESTARFM error is large, and LOTSFM is still able to maintain stable accuracy.

3.5. Fusion Time Comparison with Other Models

The fusion time consumption records for each model are listed in

Table 5. Time consumption statistics for each model in the independent validation experiments. “-” indicates that the process is not undergone.

Image Date	2017.05.07		2018.04.08		2019.09.02		2022.01.29
Fusion Stage	Train	Predict	Train	Predict	Train	Predict	Train	Predict
STARFM	-	141 s	-	142 s	-	142 s	-	146 s
ESTARFM	-	1336 s	-	1455 s	-	1421 s	-	1958 s
FSDAF	-	873 s	-	885 s	-	918 s	-	915 s
GAN-STFM	50,220 s	0.70 s	-	0.71 s	-	0.70 s	-	0.70 s
LOTSFM	247 s	0.05 s	-	0.05 s	-	0.05 s	-	0.05 s

. The three models, STARFM, ESTARFM, and FSDAF, undergo a complete fusion process each time, and the process parameters are not time-series generic. For comparison purposes, the total time consumption was noted as the prediction process consumption. Contrastingly, the GAN-STFM and LOTSFM are divided into two stages: training and prediction. The GAN-STFM time-series generic parameters are mainly neuron weights, and the LOTSFM are regression parameters.

The experimental workstation was configured with an Intel(R) Xeon(R) Silver 4110 CPU, NVIDIA Quadro P4000 GPU, and Linux system environment. STARFM, ESTARFM, FSDAF, and LOTSFM use the CPU for fusion calculations. Additionally, based on the features of LOTSFM, in which each pixel can be processed independently, multiple processes can be used to improve the fusion efficiency in the training stage. Depending on CPU performance, we used 12 processes for the experiments. The GAN-STFM model was trained for 500 iterations on a GPU, as described in a previous study [56]. According to the results in Table 5, STARFM is the least time-consuming model, 141 s, if the fusion is only performed for a certain date. LOTSFM requires a training process and is second only to STARFM in terms of time consumption, which is 248 s. When conducting multi-date experiments, LOTSFM and GAN-STFM only performed the prediction process based on the first training parameters, consuming less than 1 s.

3.6. Spatial Details Comparison with Other Models

A representative LST fusion experimental zone in Beijing was selected to compare the spatial details of each model, as shown in Figure 9. Summer Palace was a royal garden of the Qing Dynasty. The lake area in the Summer Palace is divided into three small lakes by an embankment: Kunming Hu 1, Kunming Hu 2, and Tuancheng Hu. Each of the three water surfaces comprised an island in the center of the lake, which appeared as an approximately circular normal-temperature region. The complex landscape and distinct temperature distribution in this area were used as contrasts for spatial details.

In Figure 9, image noise was observed in the FSDAF image, which is consistent with that in earlier studies that used FSDAF for LST fusion [63]. The GAN-STFM images exhibited excessive smoothness, which may be due to LST data differences, despite the use of multiscale structural similarity (MS-SSIM) to ensure visual performance. LOTSFM is intuitively the closest to Landsat, particularly in terms of the distribution of water LST in the Kunming Hu 2.

The spatial distribution of hot and cold urban areas is of concern for guiding the mitigation of thermal risks [88,89]. Therefore, representative areas with cooling effects, such as water bodies and green spaces, are selected in Figure 10 to reflect the detailed performance of each model in hot and cold areas.

Figure 10 Area 1 provides a detailed spatial representation of the lake area near the Summer Palace. Among the models, the LOTSFM low-temperature area demonstrates a better fit to the lake boundary. Moreover, the Tuancheng Hu Lake in the upper left stands out due to its relatively smaller area, resulting in a weaker cooling effect. As a result, the surface temperature of this lake does not exhibit the same deep blue hue as the other lake areas. Figure 10 Area 2 depicts a slanted square-shaped water body with a small land patch at its center. However, the models other than LOTSFM do not display the surface temperature of the central land area. Figure 10 Area 3 showcases the Forbidden City, where LOTSFM performs clearly in representing the cold area within the square-shaped moat. Figure 10 Area 4 represents the Temple of Heaven, which is predominantly covered in greenery, with two bare circular building areas clearly depicted in the LST of LOTSFM.

3.7. Error Distribution of Different Input Models

Because images with close dates have minimal land cover changes, they are generally used as reference image pairs for input; this principle was used for the experiments in Section 3.5. However, large differences in the values of the images with the closest dates among those that can be collected may also exist. The green and grey areas in the APA diagram reflect the reference image pair information. The large difference between the grey and green areas in Figure 8d is due to the excessive difference between the values of the reference pair and the predicted image, even though the reference and predicted dates were closest to each other in the sample dataset.

The extended input is intended to allow the model to synthesize the input data by feeding sufficient pairs of images to avoid artificial selection. Furthermore, more data input increases the upper limit of what can be achieved through fusion. Both STARFM and FSDAF achieved similar results when using a single image-pair input, as shown in Figure 11. ESTARFM, GAN-STFM, and LOTSFM expand the inputs. However, the expansion of ESTARFM is limited. Although it achieves good results in Figure 11 on 2019.09.02, the result on 2022.01.29 is unsatisfactory, indicating that its robustness still needs improvement. GAN-STFM showed ideal absolute difference distribution results for the two fusions on 2019.09.02 and 2022.01.29 in Figure 11, while LOTSFM performed relatively ideal on all dates.

Having a sufficient sample input is advantageous for learning-based methods; however, the fused reference image pairs that can currently be prepared for deep learning methods are still small datasets. The black-box approach used for small datasets can lead to high uncertainty in the fusion results. For instance, in Figure 11, the GAN-STFM shows concentrated image prediction errors in the middle of the region on 2017.05.07, while the other models have fused image errors that are low in the middle and high around the region. Furthermore, a small block artifact is visible in the upper left corner of the GAN-STFM result for 2022.01.29, while the other date images show no significant artifacts.

4. Discussion

4.1. Usage and Data Flexibility

In the experiments, the 36 km by 36 km area used was able to cover the main regions of the megacity. LOTSFM performs fusion in coarse pixels such that the fusion of larger or tiny areas is also easily achieved while allowing parallel computing to be easily applied to increase the speed of fusion. Based on the training and prediction processes, this allows potentially repetitive parameter calculation processes to be combined, with each subsequent fusion requiring only short prediction processes. This reduces the average computational time for multi-date fusion, enabling batch spatiotemporal fusion. The high spatial and temporal continuous LST remote sensing data generated have the potential to study urban thermal environment dynamics. For example, in industrial areas, it can be used to judge the intensity of work and monitor industrial thermal anomalies and pollution that may occur on any date. MODIS data do not accurately capture factory buildings, and Landsat data may suffer from a lack of key data. Based on the fused data available in the batch, the fine monitoring of industrial heat at a daily frequency can be achieved.

Additionally, the LOTSFM can accommodate the LST from different sensor sources. This is because (1) in Section 2.1, the coarse-resolution data used MOD11A1 and MYD11A1 from two different satellites, Terra and Aqua, respectively. Therefore, it is reasonable to believe that LOTSFM supports the fusion of other satellite LST data, such as FY-3C, ASTER, and Sentinel3. (2) In Section 3.1, even for the longest distance tested (40–50 km), a high correlation between ideal coarse pixels and fine pixels was observed (average PCCs > 0.95). Although MODIS and Landsat have coarse and fine resolutions of 1 km and 30 m, respectively, fusion with higher resolution ratios is still possible with the support of high correlation.

4.2. Fusion Performance Predictability

Based on the correlation results between the coarse and fine pixels during training, LOTSFM can predict the fused performance at different locations to reduce uncertainties. Figure 12a shows the spatial distribution of ideal correlations between the coarse and fine pixels during the training process, with high correlations in the middle of the region and low correlations around it. This is consistent with the absolute difference distribution of LOTSFM in Figure 11, particularly the results for 7 May 2017. Furthermore, based on the assumption of LOTSFM, heat transfer helps maintain stable temperature correlations between ground objects by promoting a temperature balance towards homogeneous conditions. However, the temperature of objects with high specific heat capacity is less affected by heat transfer, inhibiting its enhancement effect on the temperature correlation. Figure 12b shows a 10 m resolution land cover image of Beijing in 2017 derived from the FROM-GLC10 dataset [90]. Based on the heat transfer hypothesis of the LOTSFM, it is possible to explain the reason why the impermeable water areas with a low specific heat capacity in Figure 12b are more likely to have a high correlation and low absolute difference than the water areas.

Based on the 29 image pairs prepared for this study in Beijing, Figure 13 provides a time-series analysis of the potential influences of land cover. Figure 13a shows the time-series variations of LST for each land cover class. Distinct temperature differences exist among different land cover classes. The temperature ranks from high to low as follows: impervious surface, grassland, cropland, forest, and water bodies. While the average absolute difference (AAD) for different land cover types in Figure 13b is not highly stable along the time series, a more distinct observation can be performed by examining the box plot in Figure 13c, revealing that water bodies exhibit significantly higher AAD values and greater dispersion compared to impervious surfaces. This observation is consistent with the spatial patterns shown in Figure 12, indicating that LOTSFM has different absolute differences in fusion for different land cover classes. Based on this characteristic, it can be anticipated that LOTSFM is well-suited for monitoring thermal dynamics in urban areas characterized by the aggregation of impervious surfaces.

Furthermore, in Figure 13a, the lowest LST is observed on 2022.01.29, and the fusion AAD still exhibits favorable performance. This indicates that even under extreme temperature conditions, LOTSFM is capable of achieving stable predictions based on expanded inputs. However, in Figure 13b, a relatively large AAD is observed on 2022.05.13, which cannot be explained by land cover or correlations. Sensor differences may be the dominant factor contributing to this observation.

It should be recognized that spatial and temporal distribution of absolute differences did not completely correspond to the correlation distribution or land cover, and it was necessary to focus on the deviations caused by sensor differences. Earlier studies considered sensor differences, including bandwidth and solar geometry factors, to be constant [19,35]. In reflectance fusion, such considerations may already be justified because reflectance is relatively constant throughout the day. However, LST varies significantly throughout the day, and sensor observations may not be taken simultaneously, even on the same day. Treating the sensor difference term as a constant can lead to large deviations, which the LOTSFM considers. This is a potential reason for the performance drop in classical reflectance models when temperature–data fusion is conducted.

4.3. Limitations

LOTSFM is built upon the correlation of pixel time series. Insufficient input samples or poor correlation conditions can limit the performance of LOTSFM. In this case, the minimum input model may be the better choice. Therefore, the integration of minimum input models to form hybrid complementarities is prospected. Additionally, some minimum input models sensible to changes in ground cover are ideal for integration to address the difficulties in capturing land-cover changes.

Missing pixel filling in LOTSFM only serves to make the image meet the input integrity requirements. It is not discussed and validated in as much detail as a dedicated cloud missing reconstruction algorithm. Moreover, as with most temporal reconstruction algorithms, cloud filling in LOTSFM can only provide assumed LST values based on cloud-free conditions. Recent studies, such as the research conducted by Gong et al. [91], have proposed a reconstruction method using a nonlocality-reinforced network (NRN), examined various combinations of multimodal datasets, and evaluated them, which provides a useful reference for enhancing the missing pixel filling component of LOTSFM. Moreover, microwave data (e.g., FY-3D MWRI) is expected as input for LOTSFM to reconstruct missing pixel values closer to the real cloud conditions in the future.

The expansion of inputs in LOTSFM remains limited. We believe that a more comprehensive form of input expansion should encompass not only an increase in time series but also the incorporation of data from multiple sources. Currently, most spatiotemporal fusion models predominantly focus on space-based satellite data. Aerial data holds significant potential for integration into spatiotemporal fusion methods, especially given the rapid development of aerial unmanned drone technology. The fusion of aerial and satellite data can lead to the generation of higher-resolution and more accurate data, which has already been used for monitoring soil, crops, and forests [92,93,94]. However, while the integration of aerial data into surface temperature fusion is not yet widespread, it does not detract from recognizing it as one of the promising research directions.

5. Conclusions

High spatiotemporal resolution LST data are crucial for urban thermal environment studies; however, they are not always readily available or generated. Spatiotemporal fusion presents a potential solution; however, practical production applications face several challenges, including difficulties in sample selection, handling cloud cover, and fusion speed. Based on the above issues, a spatiotemporal fusion model, LOTSFM, based on long time-series data was proposed. The main contents and model features are summarized as follows:

• LOTSFM is a model with extended inputs, requiring only the necessary number of input samples to avoid sample selection and improve the robustness of the results.
• LOTSFM consists of two stages, training and prediction, and employs multi-process parallel computing to rapidly generate spatiotemporal fusion data in batches.
• LOTSFM utilized Julian days to estimate the sensor difference term, which was experimentally shown to significantly improve numerical accuracy.

Based on Landsat 8/9 and MODIS LST data, 79 image pairs were collected from three cities. The cross-validation results show that the average RMSE values are 1.60 °C, 2.17 °C, and 1.72 °C for Beijing, Shanghai, and Guangzhou, respectively. This capability enables effective monitoring of dynamic changes in urban heat island intensity, industrial heat emissions, and residential heat exposure. Particularly in the current global climate instability context, LOTSFM holds great potential in offering valuable guidance for proactively addressing and mitigating the detrimental impacts associated with extreme heat.

Nevertheless, it is important to acknowledge that further enhancements are required for LOTSFM. For instance, the fusion accuracy in areas with low correlation may be relatively diminished. To address this limitation, future research could explore the integration of alternative spatiotemporal fusion models to establish a complementary relationship, thereby enhancing the fusion accuracy in such areas. Conclusively, more than just the fusion of data, the complementary fusion of models will potentially be an important trend in the development of spatiotemporal fusion.