GEOSTATISTICAL SOLUTIONS FOR DOWNSCALING REMOTELY SENSED LAND SURFACE TEMPERATURE

Remotely sensed land surface temperature (LST) downscaling is an important issue in remote sensing. Geostatistical methods have shown their applicability in downscaling multi/hyperspectral images. In this paper, four geostatistical solutions, including regression kriging (RK), downscaling cokriging (DSCK), kriging with external drift (KED) and area-to-point regression kriging (ATPRK), are applied for downscaling remotely sensed LST. Their differences are analyzed theoretically and the performances are compared experimentally using a Landsat 7 ETM+ dataset. They are also compared to the classical TsHARP method. * Corresponding author


INTRODUCTION
Land surface temperature (LST) is a crucial variable that impacts on many aspects of the human and physical environments.Remotely sensed images provide good sources for estimating LST, which can be retrieved from the thermal infrared (TIR) bands (e.g., in Landsat 5 TM, Landsat 7 ETM+ or Landsat 8 TIR sensor (TIRS) data).However, the TIR band often has a coarser spatial resolution than the corresponding multispectral bands.For example, the spatial resolutions of TIR bands of Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 TIRS are 120 m, 60 m and 100 m, but the corresponding multispectral bands of the Landsat data are all at a finer spatial resolution of 30 m.In real applications, users always need to obtain more detailed spatial information of LST.Thus, it is great interest to downscale the LST to a finer spatial resolution for more reliable applications.
The availability of finer spatial resolution (30 m) Landsat multispectral bands provides excellent opportunity to be fused with the LST image to produce 30 m LST.A classical method based on TsHARP (Agam et al., 2007) was developed for this issue.However, this method produces block artifacts in predictions.In this paper, four geostatiscal solutions, including regression kriging (RK) (Mukherjee et al., 2015), downscaling cokriging (DSCK) (Atkinson et al., 2008;Pardo-Igúzquiza et al., 2006;Rodriguez-Galiano et al., 2012), kriging with external drift (KED) (Sales et al., 2013), and area-to-point regression kriging (ATPRK) (Wang et al., 2015(Wang et al., , 2016a,b),b), were applied for downscaling LST retrieved from Landsat 7 ETM+ TIR band, with the aid of the finer spatial resolution Normalized Difference Vegetation Index (NDVI) image derived from the corresponding multispectral bands.Note that RK and DSCK were already applied for downscaling LST in earlier studies (Mukherjee et al., 2015, Rodriguez-Galiano et al., 2012).It is the first time to adopt KED and ATPRK for the issue in this paper.

METHODS
In this paper, the LST image was derived from the TIR band based on the radiative transfer equation method (Rodriguez-Galiano et al., 2012).Let EG , where G is the spatial resolution (zoom) ratio between the LST and NDVI images) in the NDVI image.The notations F and C denote the fine and coarse pixels, respectively.The objective of downscaling is to predict target variables ()

F
Z x for all fine pixels.The principles of the four geostatistical solutions are described below.

ATPRK
The ATPRK method consists of two steps: regression modelling (for overall trend estimation) and area-to-point kriging (ATPK) (for downscaling residual remaining after removal of the trend).Suppose ˆ() FR Z  x and ˆ() x are predictions of the regression and ATPK parts, the ATPRK prediction is ˆˆ( ) ( ) ( ) (1) For regression prediction, it makes full use of the fine spatial resolution information in the NDVI image, and is calculated as a linear combination of the NDVI image where a and b are two coefficients.

KED
In KED, for each fine pixel, the prediction is calculated as a linear combination of the N surrounding coarse LST pixels The kriging weight i  is estimated by the kriging matrix below  x .The kriging weights need to be calculated for each fine pixel.

DSCK
In DSCK, the coarse LSTand fine NDVI are considered jointly.The prediction is a linear combination of the N neighboring coarse LST pixels and M neighboring fine NDVI pixels.
The kriging weights i  and j  are estimated by the kriging matrix in (11), where  where i  is the kriging weight for the ith neighboring coarse pixel, and the weights are calculated by the matrix in ( 13) In ( 13), ( , ) ij  xx and ( , ) j  xx are estimated from the coarse residuals directly, where the coarse or fine pixel is treated as a centroid.The RK prediction is a combination of the regression prediction in (2) and the kriging prediction in ( 12)

Differences between the geostatistical solutions
The centroid-based RK method treats each coarse or fine pixel as a centroid and ignores the spatial support of pixels, and thus, its implementation is simple.Different from RK, the other three methods (DSCK, KED and ATPRK) take into account explicitly the size of support and predict variables from areal supports to points via semivariogram deconvolution to parameterize kriging for prediction (see the semivariograms in (5), ( 9) and ( 11)).
The differences between DSCK, KED and ATPRK lie mainly in model complexity and computational cost.Specifically, DSCK is a one-stage method which considers the coarse LST and fine NDVI simultaneously, by including both auto-semivariograms and cross-semivariogram in the single kriging matrix.The kriging weights need to be calculated only once.However, both the auto-semivariogram ( , )  xx in (11) involve complex deconvolution and convolution, which would sometimes require manual intervention.
KED does not need the cross-semivariogram that relates coarse LST with fine NDVI in DSCK.Alternatively, it includes the fine NDVI in the kriging matrix.The extension of the matrix can result in unstable matrix (also the case for DSCK), especially when the covariate (i.e., fine NDVI) does not vary smoothly in space (Goovaerts, 1997).Furthermore, the KED kriging weights are calculated locally for each fine pixel, which greatly increases the computational cost, especially for large areas.
As a recently developed solution, ATPRK noticeably simplifies the model complexity and reduces the computational cost.The utility of the cross-semivariogram in DSCK (i.e., accounting for the cross-correlation between coarse LST and fine NDVI) is achieved by the simple regression model in ATPRK instead (see ( 2) and ( 3)).ATPRK requires only auto-semivariogram modelling in downscaling residuals.The size of the kriging matrix in ATPRK in ( 5) is much smaller than that in DSCK in (11).Thus, ATPRK is much easier to automate and more userfriendly than DSCK.Compared to KED, the kriging weights need to be calculated only once for the entire study area.This is because ATPRK explicitly separates trend estimation (i.e., regression modelling) from ATPK-based spatial prediction of residuals.This makes ATPRK a fast image downscaling approach and free of the risk of instability in the kriging weights calculation.
With respect to the existing TsHARP method, it is also a twostep approach.The first step also involves regression modelling (as in ATPRK and RK).In the second step of TsHARP, however, the coarse residuals are added back to fine regression predictions straightforwardly.The inconsistency in spatial resolutions between the regression prediction and residual can result in blocky artifacts in TsHARP predictions.

EXPERIMENTS
The four geostatistical methods and the TsHARP method were tested using a data covering a 15 km by 15 km region in Granada, Spain.The data were acquired on 20 July, 2002.For objective evaluation, the 60 m LST and 30 m NDVI images were degraded to 120 m and 60 m, respectively.The task was to fuse the 120 m LST with the 60 m NDVI images to produce the 60 m LST image.The available 60 m LST image was used to evaluate the accuracy.
Fig. 1 shows the results of the five methods, and Fig. 2 shows the results of a sub-area for convenience of visual comparison.
From the results, we can observe that there are noticeable blocky artifacts in the TsHARP result.Due to the unstable kriging matrix, there are also blocky artifacts in the KED result.Compared to RK, DSCK and ATPRK can reproduce more spatial variations (see the restoration of the red pixels).
The quantitative assessment is in Table 1, where five indices are used, including the root mean square error (RMSE), correlation coefficient (CC), universal image quality index (UIQI) (Wang and Bovik, 2012), relative global-dimensional synthesis error (ERGAS), and coherence.Coherence (quantified by the CC) is an index measuring the relation between the coarse LST and the coarse LST obtained by upscaling the LST downscaling predictions.As seen from Table 1, ATPRK can produce more accurate results than the other methods, and can perfectly preserve the original coarse LST data.Although the TsHARP prediction has perfect coherence with the coarse LST, the RMSE and ERGAS are larger, and the CC and UIQI are smaller than that of DSCK, KED and ATRPK.Failing to account for the size of supports, RK leads to the smallest coherence.

CONCLUSION
This paper studies the use of the four geostatistical solutions (RK, DSCK, KED and ATPRK) for downscaling LST and compares their performances based on a case study for the Landsat 7 ETM+ dataset.RK fails to account for the change of support in downscaling and produces less accurate result than the other three geostatistical solutions.Their accuracies are generally greater than the classical TsHARP method.Amongst all methods, ATPRK can produce the most accurate result, and can perfectly preserve the original coarse LST data.

FFFF
-to-coarse residual semivariogram between fine and coarse pixels centered at x and j x , and  is the Lagrange multiplier.Suppose s is the Euclidean distance between the centroids of any two pixels.The semivariograms in (5) are calculated by convolving the fine-to-fine residual semivariogram () s is derived by deconvolution of the coarse residual semivariogram of the coarse residual image () R x .Details on deconvolution can be found in Wang et al. (2016a).


xx have the same meanings as those in (5).() Ci Y x is the NDVI for coarse pixel centered at i


xx are coarse and fine auto-semivariograms estimated from the original coarse LST pixels and fine NDVI pixels, respectively.xx is the fine-to-coarse semivariogram between fine and coarse LST pixels.xx need to be estimated by deconvolution and convolution processes.

Table 1 .
Quantitative assessment on the methods for the entire study area