Land-surface temperature retrieval at high spatial and temporal resolutions based on multi-sensor fusion

Abstract

Land-surface temperature (LST) is of great significance for the estimation of radiation and energy budgets associated with land-surface processes. However, the available satellite LST products have either low spatial resolution or low temporal resolution, which constrains their potential applications. This paper proposes a spatiotemporal fusion method for retrieving LST at high spatial and temporal resolutions. One important characteristic of the proposed method is the consideration of the sensor observation differences between different land-cover types. The other main contribution is that the spatial correlations between different pixels are effectively considered by the use of a variation-based model. The method was tested and assessed quantitatively using the different sensors of Landsat TM/ETM+, moderate resolution imaging spectroradiometer and the geostationary operational environmental satellite imager. The validation results indicate that the proposed multisensor fusion method is accurate to about 2.5 K.

Keywords:

• land-surface temperature
• image fusion
• remote sensing
• resolution

1. Introduction

Land-surface temperature (LST) retrieved from remote-sensing thermal infrared observations has proved to be vital in applications such as determining the land-surface energy exchange with the atmosphere (Sellers et al. 1997; Xu, Liang, and Liu 2011), the impact of urbanization on global climate and environment change (Jin, Dickinson, and Zhang 2005; Yu et al. 2009), atmospheric stability and local wind systems (Tarpley 1979), land-surface characteristics (Li and Avissar 1994), and evapotranspiration (Anderson et al. 2012). However, the available satellite LST products have either low spatial resolution or low temporal resolution, which constrains their potential applications. For instance, the Landsat satellites, which retrieve LST at high spatial resolutions (~100 m), can only revisit the same location on earth at intervals of half to one month. The moderate resolution imaging spectroradiometer (MODIS) can make repeated LST observations in one day or a half day, but with a 1 km spatial resolution (Wan and Li 1997; Wan et al. 2004). Higher temporal sampling of LST is achievable with the geostationary operational environmental satellite (GOES), but with a coarser spatial resolution (about 4 km) (Sun and Pinker 2005; Inamdar et al. 2008). As a result, it is extremely difficult to acquire satellite images with high temporal and spatial resolutions due to trade-offs among these resolutions (Zhan et al. 2013). One possible cost-effective solution is to explore the data fusion methods that can blend the two types of LST products from the different sensors to generate synthetic LST products with high resolutions in both space and time.

A large number of data fusion methods have been developed (for a review see Ehlers 1991; Hall and Llinas 1997; Pohl and Van Genderen 1998; Luo, Chih-Chen, and Kuo Lan 2002). Classical remote-sensing image fusion techniques include panchromatic/multispectral fusion (Wang et al. 2005; Thomas et al. 2008), multispectral/hyperspectral fusion (Eismann and Hardie 2005), and multitemporal fusion (Shen et al. 2009). These fusion techniques are aimed at extracting information and enhancing visual effects; however, they are not effective in simultaneously enhancing temporal resolution and spatial resolution, and cannot provide high-spatial-resolution observations with a dense time series.

Only in recent years have we witnessed the emergence of spatiotemporal fusion models for predicting an unknown high-resolution image. For example, Inamdar et al. (2008) proposed an Normalized Difference Vegetation Index (NDVI)-based fusion method by blending observations from MODIS and GOES for half-hourly LST values, at a 1 km spatial resolution. However, the method needs auxiliary NDVI data, and LST with a 1 km spatial resolution limits the sensors’ ability to quantify land-surface processes in heterogeneous landscapes. Gao et al. (2006) proposed a spatial and temporal adaptive reflectance fusion model (STARFM) that does not require any auxiliary data. This approach was first used to blend Landsat and MODIS data to predict daily surface reflectance at Landsat spatial resolution and MODIS temporal frequency (Gao et al. 2006). Since then, the STARFM framework has also been applied to the generation of gross primary productivity (Singh 2011), the analysis of dryland forest phenology (Walker et al. 2012), the examination of virus dissemination (Liu and Weng 2012), and the estimation of daily evapotranspiration (Anderson et al. 2011).

In the study presented herein, we extend the STARFM framework (Gao et al. 2006) and present a spatiotemporal fusion method for remote-sensing LST by blending observations from multiple sensors (Landsat TM/ETM+, MODIS, and the GOES Imager). One of the significant improvements of this method is the consideration of sensor observation differences based on within-class regression, with M-estimation (Huber 1964) to improve the fusion accuracy. Meanwhile, the spatial correlations between different pixels are considered using a variation-based prior constraint model. Based on the proposed method, our objective was to obtain daily and/or hourly high-spatial-resolution LST products by blending observations from Landsat with MODIS and GOES.

2. Methodology

2.1. STARFM

The premise of STARFM is that the preprocessed and corrected reflectances of different sensors are consistent and comparable (Gao et al. 2006, Masek et al. 2006). STARFM predicts reflectances based upon a spatially weighted difference computed between the Landsat and the MODIS scenes at the acquisition date and one or more MODIS scenes at the prediction date (Gao et al. 2006). The prediction algorithm in STARFM is given by the following:

(1)

where F and C denote the fine-resolution and coarse-resolution reflectances, respectively, (x_i,y_j) denotes the pixel location, t_k is the acquisition date, k∊(1,∞), t₂ is the prediction date, w is the size of the moving window, and (x_w/2,y_w/2) is the central pixel.

The spatial weighting function W_ijk determines how much each similar pixel contributes to the central pixels. Ensuring that the sum of all the weights is 1, W_ijk can be normalized as follows:

(2)

where E_ijk may be related to the spectral difference, temporal difference or geometric distance in the moving window, and can be expressed as follows (Gao et al. 2006):

(3)

where A is a scale factor, equal to 10,000 (Gao et al. 2006), S_ijk represents the spectral difference, T_ijk represents the temporal difference, and D_ijk represents the geometric distance. Equations (2) and (3) show that one pixel may be assigned a higher weight W_ijk if the pixel has smaller values of S_ijk, T_ijk, and D_ijk.

2.2. The preliminary calculation function

Although the observation data can be preprocessed before the fusion process, systematic observation differences may still exist between the different sensors. The observation differences may come from the differences in the spectral ranges and/or calibration algorithms. To minimize the effects of the observation differences on the fusion results, this study proposes to fit regression curves for different land-cover types when measuring the radiance difference (similar to the spectral difference when fusing reflectance as STARFM).

To consider the observation differences, a simple linear regression relationship between the coarse-resolution image and the corresponding fine-resolution image at time t can be described as follows:

(4)

where Fand C (and subsequently) denote the fine-resolution LST and coarse-resolution LST, a is the gain parameter of the linear regression function between the coarse- and fine-resolution pixels, and b is the bias. Thus, the radiance difference can be calculated as follows:

(5)

However, the linear regression coefficients (a and b) between the two sensors may be different for different land-cover types. The radiance difference, Equation (5), can be rewritten as follows:

(6)

where the subscript c represents the land-cover types which can be obtained by classifying the fine-resolution image.

For convenience, the land-cover map can be determined using an unsupervised classification algorithm. The linear regression coefficients for each land-cover type can then be solved by various optimization estimation methods. The simplest method is the least-squares approach; however, the least-squares method is sensitive to the existence of outliers and the distribution assumption of noise. When the noise distribution is not normal or outliers exist in the observations, the estimation may be unsatisfactory (Wu and Tam 2001). To improve the accuracy of the estimated parameters, this paper uses the M-estimation method with the Huber cost function (Huber 1964) to solve the coefficients.

In Equation (1), STARFM utilizes a weighting function to incorporate the spatial correlation information from the neighboring coarse- and fine-resolution image pixels. The spatial correlations used by STARFM are very useful for reflectance data, but may introduce a degraded effect for LST. Therefore, the ‘correct’ neighboring pixels and a more suitable window size need to be selected in the proposed method. Meanwhile, the STARFM algorithm allows for the input of either one or two base image pairs. For the predicted image at t ₂, one base image pair may be at t₁ or t₃, while the two base image pairs may be at t₁ and t₃. According to the different input base image pairs, the preliminary calculation functions may be expressed as the following two cases:

Case 1: A single image pair as the t₁ or t₃ input, and we eliminate the temporal difference within the spatial weighting function (Hilker et al. 2009). Therefore, the spatial weighting function E_ij is only related to the radiance difference S_ijc and the geometric distance D_ij (Hilker et al. 2009), and can be expressed as follows:

(7)

Here, D_ij can be easily obtained by the method of Gao et al. (2006), and S_ijc denotes the radiance difference considering the observation difference. W_ijc can be normalized as follows:

(8)

Equations (7) and (8) mean that a smaller distance and a smaller radiance difference from the similar pixel to the central pixel should produce a higher weight. The fine-resolution LST can be calculated using the following:

(9)

Case 2: Two image pairs as the t₁ and t₃ inputs, and we adopt a temporal weight (Zhu et al. 2010). Either of the fine-resolution LSTs at t₁ or t₃ can be used as the base date to predict the fine-resolution LST of the prediction date t₂, and they are marked as F(x_w/2, y_w/2, t_1–2) and F(x_w/2, y_w/2, t_3–2), respectively. A consideration when obtaining LST at t₂ is to use a temporal weighted combination of the two prediction results (Zhu et al. 2010). The temporal weight can be expressed as follows (Zhu et al. 2010):

(10)

Equation (10) means that with a smaller change of resampled coarse-resolution LST within the window w between the time t_k (k=1 or 3) and t₂, a larger temporal weight should be assigned. The predicted fine-resolution LST at the prediction time t₂ is then calculated as follows:

(11)

2.3. Variation-based model

By far the majority of the research into similar fusion problems has been based on filter-based methods, in which the spatial correlations between different pixels are not effectively considered. In this study, the spatial correlations are considered by using a variation-based method. Compared with the filter-based fusion methods, we can solve all the unknowns as a whole and consider the correlations between them. As Ludusan and Lavialle (2012) noted, the variation-based fusion methods are more robust, and have proved to be more suitable for practical image-fusion applications. We present the following variation-based fusion model:

(12)

where the first term is the data fidelity term, acts as the regularization spatial constraints term, λ can be called the regularization parameter, is the final result at the prediction time t₂, and is a function of the fine-resolution LST at t₁ and coarse-resolution LSTs at t₁ and t₂. The spatial constraints can be expressed as follows:

(13)

where ϕ₁ is the number of pixels of the whole image and p is the pixel of the moving window. ρ(·) is a constraints prior, and when it is a quadratic potential function, as in Equation (14), the corresponding prior is regarded as Gauss–Markov (Shen and Zhang 2009):

(14)

As for d_p(Z_i), the second-order differences are computed between similar pixels and the center pixel within the moving window in the image:

(15)

Complete details of the computing processes can be found in Shen and Zhang (2009).

It is noted that the initial value should, theoretically, be updated with the spatial weighting function. In this paper, however, it is fixed due to the computational load consideration. It is therefore similar to implementing a variational denoising process on the initial results solved by the method in Section 2.2.

3. Data

The primary streams of data used in this study comprise the parameters retrieved from cloud-free MODIS, Landsat, and GOES images. These are described below.

3.1. MODIS 1 km LST

MODIS is an EOS instrument on both the Terra and Aqua satellites that serves as the keystone for global studies of the atmosphere, land, and ocean processes, scanning ±55° from the nadir in 36 bands ranging from visible to thermal infrared regions (TIR) (Wan and Li 1997). The bands in the thermal atmospheric window are designed for the remote sensing of surface temperature and emissivity. The LST product was proposed by the MODIS Science Team as a daily daytime and nighttime 1 km global land product, and includes derivative products at lower temporal frequencies and spatial resolutions. The accuracy of the MODIS daily LST is better than 1°C in the range from 10 to 50°C, as validated from in situ measurement data collected in field campaigns between 2000 and 2002 (Wan et al. 2004).

We use the MODIS LST/Emissivity, Daily, Level 3, Global 1 km SIN grid product known as MOD11_A1, which is available from the International Scientific Data Service Platform (http://datamirror.csdb.cn) and The Next Generation Earth Science Discovery Tool (http://reverb.echo.nasa.gov/reverb/). The 1 km LST/Emissivity data are produced by a generalized split-window algorithm that uses radiances from the MODIS TIR window channels 31 (10.8–11.3 µm) and 32 (11.8–12.3 µm) (Wan and Dozier 1996).

3.2. Landsat TM/ETM + LST

Landsat images offer the longest continuous global record of the earth's surface and are a unique resource for global change research and applications in agriculture, cartography, geology, forestry, surveillance, education, and national security (Güçlüer, Bayram, and Maktav 2010). Landsat offers significant advantages over MODIS in the spatial resolution. The Landsat TIR channel has about a 100 m spatial resolution (Landsat TM = 120 m, Landsat ETM + = 60 m) and is available from The Next Generation Earth Science Discovery Tool (http://reverb.echo.nasa.gov/reverb/). This is an important spatial resolution because it is coarse enough for global coverage, yet detailed enough to characterize human-scale processes such as the dynamics of urban growth (Masek, Lindsay, and Goward 2000) and urban heat islands (Kumak, Bhaskar, and Padmakumari 2012). We retrieve the Landsat LST using a generalized single-channel method (Jimenez-Munoz and Sobrino 2003) from the Landsat TIR window channel 6 (10.4–12.5 µm).

3.3. GOES imager LST

GOES provides approximately half-hourly temporal resolution earth observation data, which constitute an important element of the land-surface processes. The GOES 10 and GOES 13 Imager instruments consist of five spectral channels (one visible band and four infrared bands) (Sun and Pinker 2003; Sun, Pinker, and Basara 2004), which measure the radiant and reflected solar energy from sampled areas of the earth. GOES Imager data are available in a web-based data archive and distribution system known as NOAA's Comprehensive Large Array-data Stewardship System (http://www.class.noaa.gov/nsaa/products). For the present study, the infrared image data at a nadir resolution of about 4 km stored as GOES Variable Format (GVAR) counts packaged in 10-bit words in Network Common Data Format are used. The 10-bit (0–1023) GVAR count value is converted to brightness temperatures for the corresponding channels, utilizing the calibration coefficients (Weinreb, Johnson, and Han 2011) for the GOES Imager. In this study, the LST data at an hourly or half-hourly temporal resolution from GOES have been retrieved according to the corresponding algorithms (Sun and Pinker 2003; Sun, Pinker, and Basara 2004).

3.4. Ground truth data

The network of surface radiation measurement sites (SURFARD) was established in 1993 through the support of the NOAA Office of Global Programs. Its primary objective is to support climate research with accurate, continuous, long-term measurements pertaining to the surface radiation budget over the United States (Wang and Liang 2009). The Sioux Falls station (SXF), which covers crop data from South Dakota (96.62°W, 43.73°N), and the Desert Rock station (DRA), which covers sand data from Nevada (116.02°W, 36.62°N), are selected in this study. The data change from a 3-min to 1-min resolution after 1 January 2009, and can be downloaded from http://www.esrl.noaa.gov/gmd/dv/data/. The site provides continuous radiometric instrumentation to measure upwelling and downwelling thermal infrared irradiance, which are used to retrieve the ground truth of the LST, according to Inamdar et al. (2008). As in the previous studies (Sun and Pinker 2003; Sun, Pinker, and Basara 2004; Inamdar et al. 2008), the ground truth of the LST is used to directly evaluate the fused results from the GOES and TM LST products.

4. Experimental results

In this section, we test the performance of the proposed method by fusing multiple LST products from the different sensors. The MODIS LST and GOES Imager LST products are reprojected to the UTM projection and resampled to the same size as the Landsat TM/ETM + LST products using a cubic interpolation approach (with available georegistration information). To evaluate the quality of the fusion results quantitatively and visually, several representative metrics are employed in this study. The average absolute difference (AAD) and the root mean squared error (RMSE) between the predicted synthetic LST and the actual LST are shown to directly reflect the deviation of the predicted LST. The scatter plots of the predicted against actual images are also shown, and provide an intuitive comparison between the estimated and actual LST. In addition, the ground truths from the SXF data are also used to verify the effectiveness of the proposed method.

4.1. Fusing of MODIS and ETM + LST products

In the first series of experiments, the proposed method is tested on sub-images from Beijing and Hubei, China. Figure 1 shows the MODIS (upper row) and Landsat-7 ETM+(second and third rows) LSTs from 11 November 2001, 29 December 2001, and 15 February 2002, respectively, from left to right. Using STARFM and the proposed method, we reconstruct the Landsat-7 ETM + LST from 29 December 2001, given the two LST pairs from 11 November 2001 (Figure 1a and d), and 15 February 2002 (Figure 1c and f) and its MODIS counterpart (Figure 1b). The actual ETM + LST image (Figure 1e) can be used to evaluate the predicted result. The images of the third row in Figure 1 are the colored LSTs, corresponding to the second row, which show the changes in LST between the different months. Figure 1j and l are the predictions using the one input pair of 11 November 2001, and 15 February 2002, respectively, while Figure 1k is the prediction using the two input pairs from both November and February.

Figure 1. The actual and predicted LSTs. (a), (b), and (c) are the actual MODIS LST on 11 November 2001, 29 December 2001, and 15 February 2002, respectively. (d), (e), and (f) are the actual ETM + LST on 11 November 2001, 29 December 2001, and 15 February 2002, respectively. (g), (h), and (i) are the actual ETM + LST, but colored corresponding to (d), (e), and (f), respectively. (j), (k), and (l) are the predictions using one input date pair of November 11 as the base data, two input date pairs of November 11 and February 15 as the base data, and one input date pair of February 15 as the base data, respectively.

The effect of the two fine-resolution images at t₁ and t₃ for the center-predicted image at t₂ may be different. It is possible that incorrect information may be brought in when a disturbance event results in land-cover type changes (e.g. fires and other human activities). Gao et al. (2006) used both one input pair and two input pairs to predict the fine-resolution image, and did obtain reliable predicted results. However, that study did not compare the predictions. Similarly, some applications based on a single pair of inputs (Hilker et al. 2009; Singh 2011; Walker et al. 2012; Watts et al. 2011) or two pairs of inputs (Anderson et al. 2011; Liu and Weng 2012) with the STARFM method have shown that the results can meet the application demands, but also without a comparison. In another study (Zhu et al. 2010), the authors discovered that STARFM works better for reflectance with one input date pair rather than two input pairs, because of the similar phenology between the one input date pair and the predicted data. The image predicted by STARFM using two input pairs resulted in an unrealistic image due to the large differences between the two input pairs (Zhu et al. 2010). Whether the input pairs at t₁ and/or t₃ may bring about a distinguishing effect when fusing LST has not been researched to date. This study quantitatively compares the predictions of LST by the use of both one input date pair and two input pairs.

Here, we set the data of 11 November 2001, 29 December 2001, and 15 February 2002, as t₁, t₂, and t₃, respectively. Although there are no distinct visual differences between Figure 1j–l, some detailed regions in Figure 2 are selected for the convenience of visual judgment. Figure 2 shows the comparisons between the actual and predicted LST. Figure 2a–c are the predicted LST at t₂ by STARFM using the data at t₁, the data at t_3, and the data at t₁ and t₃, respectively, from left to right. Figure 2d–f are the predicted LSTs at t₂ by the proposed method using the same data as the upper row, and (g) is the actual LST at t₂. It can be seen that the predictions of the proposed method are better than those of STARFM, when compared with the actual LST. Meanwhile, the predictions using the two input image pairs at t₁ and t₃ are better than the predictions using one input image pair at t₃, but are inferior to the predictions using one input image pair at t₁. The results indicate that the input pair at t₃ may cause larger uncertainties for estimating LST changes at the 100 m scale from 1 km information.

Figure 2. Comparisons between actual and predicted LSTs. (a), (b), and (c) are predicted LSTs at 29 December 2001, by STARFM, using the data from 11 November 2001; the data from 11 November 2001, and 15 February 2002; and the data from 15 February 2002, respectively, from left to right. (d), (e), and (f) are predicted LSTs on 29 December 2001, by the proposed method using the same data as the upper row. (g) is the actual LST on 29 December 2001.

Figure 3 shows the scatter plots and histogram between the actual observed and predicted ETM + LST. Taking the predictions of the two input image pairs as an example, Figure 3a and b are the correlations between the actual and predicted LSTs of STARFM and the proposed method, respectively. The data of Figure 3b are closer to the 1-1 line than those of Figure 3a, indicating that the proposed method improves the prediction of changes between the December and November/February images. The histogram of the percentage probability distribution versus LST difference for STARFM and the proposed method is shown in Figure 3c. Compared to STARFM, the proposed method improves the mean LST difference from 0.35 to 0.24 K.

Figure 3. Comparisons between the actual LST and the predicted LST using scatter plots by STARFM (a), scatter plots by the proposed method (b), and a histogram of the predicted minus observed LST for each method (c), respectively.

To quantitatively compare the proposed method with STARFM for the different input pairs, the quantitative values of AAD and RMSE are listed in Table 1. The AAD and RMSE values of the proposed method are lower than those of STARFM, indicating that the proposed method can reconstruct the ETM + LST more precisely than STARFM. The prediction using one input date pair from 11 November 2001 obtains higher evaluation scores than the others, which again agrees with Figure 2.

Table 1. AAD and RMSE between the predictions from different input data and actual LST of the first series of experiments.

The base data as input		11/11/2001	02/15/2002	11/11/2001 and 02/15/2002
AAD (K)	STARFM	1.1826	1.5901	1.2496
	Proposed	1.1642	1.5257	1.1880
RMSE (K)	STARFM	1.4961	2.0014	1.6175
	Proposed	1.4714	1.9494	1.4986

To further test the performance of the proposed method, a second series of experiments with a 48×48 km subset from Hubei province, China, are conducted. Figure 4 shows the MODIS (upper row) and Landsat-7 ETM+(second row and bottom row) LSTs from July 9, October 13, and November 30, 2002, respectively, from left to right. Using STARFM and the proposed method, we reconstruct the Landsat-7 ETM + LST from 13 October 2002, based on one input image pair from July 9 or November 30 and two input image pairs from July 9 and November 30. Figure 5 shows the STARFM (Figure 5a–c) and proposed method (Figure 5d–f) colored predictions in comparison with the actual Landsat ETM + LST (Figure 5g). From Figure 5, we can see that better predictions can be obtained with only one pair from November 30 as the base data. The quantitative evaluation results shown in Table 2 demonstrate that the predictions using one pair of data from November 30 are better than the others. The highest AAD and RMSE values appear with the predictions from the use of one pair of data from July 9, which suggests that disturbance events may have happened between July 9 and October 13. Furthermore, the predictions of the proposed method also have lower AAD and RMSE values than those of STARFM. Similarly, the scatter plots in Figure 6 show that the predicted LST using one input pair from November 30 (Figure 6b) more closely matches the actual LST (1:1 line) than the results from the use of one input pair from July 9 (Figure 6a) and the results from the use of two input pairs on July 9 and November 30 (Figure 6c).

Figure 4. The actual Landsat and MODIS LST. (a), (b), and (c) are the MODIS LSTs on July 9, October 13, and November 30, 2002, respectively. (d), (e), and (f) are the Landsat TM LSTs on July 9, October 13, and November 30, 2002, respectively.

Figure 5. Comparisons between the actual and predicted LSTs (colored). (a), (b), and (c) are the predicted LSTs on October 13 by STARFM using the data from July 9, the data from November 30, and the data from July 9 and November 30, respectively, from left to right. (d), (e), and (f) are the predicted LSTs on October 13 by the proposed method using the same data as the upper row. (g) is the actual LST on October 13 (the change is in the range of 290–310 K).

Figure 6. Scatter plots of predicted LSTs using the proposed method against actual LST for the different input data pairs. (a) is the scatter plot of the predicted LST using the one input data pair from July 9 against the actual LST. (b) is the scatter plot of the predicted LST using the one input data pair from November 30 against the actual LST. (c) is the scatter plot of the predicted LST using the two input data pairs from July 9 and November 30 against the actual LST.

Table 2. AAD and RMSE between the predictions from different input data and actual LST of the second series of experiments.

The base data as input		07/09	11/30	07/09 and 11/30
AAD (K)	STARFM	3.4262	0.7511	1.6869
	Proposed	3.2759	0.7295	1.5771
RMSE (K)	STARFM	3.6018	0.9523	2.0714
	Proposed	3.3158	0.9159	1.8099

4.2. Fusing of the GOES and Landsat TM LST products

In the third series of experiments, the proposed method is tested on sub-images from South Dakota, USA. From the downloaded header file of the Landsat TM data, we can see that the time of the center scan was at 17:14:23 UTC time. The GOES 13 data was observed between 10:00 and 22:00 UTC time on 4 September 2010. All the TM and GOES LSTs are listed in Figure 7. The TM (Figure 7a) and the GOES LSTs observed at around 17:00 UTC time (Figure 7i) are the input base data, and the other GOES LST (Figure 7b–n) is the input data of the prediction time. The predictions from 10:00 to 22:00 UTC time are shown in Figure 8. Although we cannot give a quantitative comparison, on account of the lack of corresponding ground truth, the predictions can show the interdiurnal temperature variation, as shown in Figure 8b–m. In particular, we have no need to obtain the prediction at around 17:00 UTC time because the actual TM LST (Figure 8a) is observed at a similar time. Figure 9 shows the average LST of Figure 8 between 10:00 and 22:00 UTC time.

Figure 7. The actual Landsat LST and GOES LST. (a) is the Landsat LST observed on 4 September 2010; and (b)–(n) are the GOES LSTs observed between 10:00 and 22:00 UTC time during 4 September 2010.

Figure 8. Comparisons between the actual TM and predicted LST (colored). (a) is the actual Landsat TM LST (colored), and (b)–(n) are the predictions (colored) at Landsat spatial resolution (‘TM’) between 10:00 and 22:00 UTC time on 4 September 2010, respectively.

Figure 9. The average LST of Figure 8 between 10:00 and 22:00 UTC time on 4 September 2010.

To validate the performance of the proposed method, the fourth set of experiments utilizes the ground truth from the SXF and DRA data. The TM LST and the GOES 13 LSTs (images containing the SXF site) from 00:00 to 04:00 UTC time and 15:00 to 23:00 on 18 June 2010 are listed in Figure 10 (the GOES observations between 05:00 and 14:00 UTC time are contaminated by clouds or otherwise reduced visibility). The TM LST (Figure 10a) and the GOES LST observed at around 17:00 UTC time (Figure 10i) are the base data, and the other 13 GOES LSTs (Figure 10b–n) are the input data of the prediction time. Figure 11b–n are the predictions from 00:00 to 04:00 UTC time and 15:00 to 23:00 UTC time. Similarly, we have no need to obtain the prediction at 17:00 UTC time because the actual TM LST (Figure 11a) is observed at this time.

Figure 10. The actual Landsat LST and GOES LST. (a) is the Landsat LST observed on 18 June 2010, and (b)–(o) are the GOES LSTs observed between 00:00 and 04:00 UTC time and 15:00 to 23:00 UTC time on 18 June 2010, respectively.

Figure 11. Comparisons between the actual TM and predicted LST (colored). (a) is the actual Landsat LST at 17:00 UTC time (colored), and (b)–(n) are the predictions (colored) at Landsat TM spatial resolution (‘TM’) from 00:00 to 04:00 UTC time and 15:00 to 23:00 UTC time on 18 June 2010, respectively.

Evaluations of the LSTs from the SXF site and the predictions from 18 June 2010 are shown in Figure 12. The time units shown on the abscissa (in UTC) in Figure 12a comprise a single day (24 h). The dotted curve in red in Figure 12a represents the LST measurements at the SXF site, while the predicted LSTs are denoted by a circle in blue. Note that the solid circle at 17:00 UTC time represents the observed TM LST based on the work of Jimenez-Munoz and Sobrino (2003). Some values are interrupted where the GOES observations are lacking due to clouds or otherwise reduced visibility. The LSTs of the predictions show a downward trend between 00:00 and 04:00 UTC time and then show an increasing trend followed by a downward trend between 15:00 and 23:00 UTC time. The trends of the available predictions are in agreement with the ground measurement. The scatter plot of the predicted LST versus in situ LST is shown in Figure 12b with 0.95 for R² and 2.14 K for RMSE.

Figure 12. Evaluation of estimates of LST based on proposed method with the SXF site observations. (a) is the predicted LST (circles in blue) compared directly with the SXF observations (dotted line in red) on 18 June 2010, the solid circle at 17:00 UTC time represents the actual TM LST, the data gaps in GOES from 05:00 to 14:00 UTC time, due to clouds or otherwise reduced visibility, are not predicted, (b) the scatter plot of the predicted LST versus in situ LST.

Similar experiments are also undertaken for a different surface type using the data from the DRA site. The retrieved ETM + LST and the GOES 10 LSTs (images containing the DRA site) observed on 20 August 2002 are used. There are 45 GOES LSTs with 30-min resolutions over that day (the GOES observations between 08:30 and 09:30 UTC times are absent). The ETM + LST and the GOES LST observed at around 18:00 UTC time are the base data, and the other GOES LSTs are the input data of the prediction time. The 46 inputs and the 44 predictions are not listed here; however, the evaluations between the LSTs from the DRA site and the predictions are shown in Figure 13. The time units shown on the abscissa (in UTC) in Figure 13a comprise a single day (24 h). The dotted curve in red in Figure 13a represents the LST measurements at the DRA site, while the predicted LSTs are denoted by a circle in blue. The scatter plot of the predicted LST versus the DRA observations is shown in Figure 13b with 0.98 for R² and 1.87 K for RMSE. From Figures 12 and 13, we can see that the proposed method can obtain hourly and even half-hourly high-spatial-resolution LST, although the absolute mean errors of the predictions are about 2.5 K. The sources of errors may include the following:

1. The data preprocessing contains registration errors, and cloud contamination may also have a significant negative impact on the validation process.

2. The retrieval methods of the different sensors (see Section 3). For the same time and area, the LSTs retrieved from different sensors may be different. The retrieval methods themselves may also have certain errors.

3. The land surface is typically heterogeneous over satellite pixel areas, while in situ LSTs are usually collected over significantly smaller and more homogeneous areas (Yu et al. 2012). The scale difference may be a major source of error (Wang and Liang 2009).

Figure 13. Evaluation of estimates of LST based on proposed method with the DRA site observations. (a) is the predicted LST (circles in blue) compared directly with the DRA observations (dotted line in red) on 20 August 2002, the solid circle at 18:00 UTC time represents the actual ETM + LST, the data gaps in GOES from 08:30 to 10:00 UTC time, due to the GOES observations being absent, are not predicted, (b) the scatter plot of the predicted LST versus in situ LST.

5. Conclusion

This paper presents an extended spatiotemporal fusion method for remote-sensing LST by blending observations from multiple sensors (Landsat, MODIS, and GOES). The objective was to obtain daily or hourly LST at a Landsat spatial scale. Through consideration of the observation differences between the sensors and the spatial correlations between pixels, the proposed method can improve the current state of temporal surface temperature prediction. We also compare the predictions based on different input LST pairs, and the comparison reveals that the input LST pairs having large differences with the predicted LST result in incorrect contributory information, the same as with the predicted reflectance. Although the predictions are improved somewhat when the initial value is fixed, further improvement will be realizable when it is updated iteratively. In addition, the hourly LSTs at the Landsat spatial scale can be obtained with our method only if the Landsat data are available at the correct revisit cycle. This limiting factor could be solved by simultaneously fusing GOES, MODIS, and Landsat data.

Acknowledgements

This work was supported by the Major State Basic Research Development Program (973 Program) under Grant 2011CB707103, National High Technology Research and Development Program (863 Program) under Grant 2013AA12A301, National Natural Science Foundation of China under Grant 41271376, the Hubei Natural Science Foundation under Grant 2011CDA096, and the Fundamental Research Funds for the Central Universities under Grant 2012205020205. Many thanks to the anonymous reviewers. We would like to thank F. Gao, F. Li, and D. L. Sun for their help.