Downscaling Thermal Infrared Radiance for Subpixel Land Surface Temperature Retrieval.

Land surface temperature (LST) retrieved from satellite thermal sensors often consists of mixed temperature components. Retrieving subpixel LST is therefore needed in various environmental and ecological studies. In this paper, we developed two methods for downscaling coarse resolution thermal infrared (TIR) radiance for the purpose of subpixel temperature retrieval. The first method was developed on the basis of a scale-invariant physical model on TIR radiance. The second method was based on a statistical relationship between TIR radiance and land cover fraction at high spatial resolution. The two methods were applied to downscale simulated 990-m ASTER TIR data to 90-m resolution. When validated against the original 90-m ASTER TIR data, the results revealed that both downscaling methods were successful in capturing the general patterns of the original data and resolving considerable spatial details. Further quantitative assessments indicated a strong agreement between the true values and the estimated values by both methods.

Introduction

Satellite thermal infrared (TIR) imagery is the primary source to retrieve land surface temperature (LST) for various ecological and environmental studies at regional and global scales [1, 2, 3]. Finer-scale studies such as the analysis of urban heat island [4] and stream temperatures [5, 6], however, require LST at a higher spatial resolution than that provided by the currently available satellite thermal sensors. Moreover, due to the heterogeneous nature of land surfaces, mixed pixels of multiple anisothermal objects exist in TIR imagery at a relatively coarse spatial resolution. As a result, LST products retrieved from satellite TIR imagery are often composed of a mixture of different temperature components [7]. Therefore, there is a need for enhancing the spatial resolution of the current LST products. In particular, as coarse spatial resolution thermal sensors such as MODerate resolution Imaging Spectro-radiometer (MODIS) and Advanced Very High-resolution Radiometer (AVHRR) provide LST products with very high temporal frequencies, improvements on the spatial resolutions of these frequently used LST products will significantly augment their potential uses in many applications requiring LST with higher spatial resolutions.

Two basic approaches have been developed to enhance the LST products retrieved from satellite TIR imagery. The first approach uses spectral mixture analysis (SMA) to decompose mixed TIR pixels into multiple isothermal components [8, 9, 10, 11]. By modeling TIR radiance as a mixture of multiple thermal endmembers, SMA produces a set of component temperatures and their corresponding fractions for each mixed pixels. However, the spatial details within the mixed pixels remain unresolved in the decomposed component temperatures. In contrast, the second approach uses downscaling methods to disaggregate each mixed TIR pixel into multiple finer resolution TIR pixels for subpixel LST retrieval [5]. Specifically, it involves two steps: 1) downscaling TIR radiance using higher spatial resolution ancillary data (e.g. land cover maps), and 2) retrieving subpixel surface temperature from downscaled TIR radiance. Compared to the first approach, the second approach can provide substantial spatial details within mixed pixels.

In this paper, we present a downscaling approach for enhancing LST products using satellite TIR imagery. Given the fact that the algorithms for retrieving LST from TIR radiance are well established [12, 13, 14, 15], our focus is on the development of methods for downscaling TIR radiance from satellite sensors. Liu et al. [7] grouped scaling approaches to satellite-derived data into two broad categories: 1) scaling based on physical models, and 2) scaling based on statistical methods. In this paper, we develop both physical and statistical downscaling methods for subpixel TIR radiance estimation.

Downscaling Methods

Physical downscaling

The basic idea of the physical downscaling method is to establish a functional relationship between TIR radiance and ancillary data, which satisfies two conditions: (1) the functional relationship is physically meaningful and holds across different scales; and (2) the ancillary data can be easily scaled. Specifically, we modeled the TIR radiance as a linear combination of multiple land cover components weighted by their corresponding fractions and modified by atmospheric effects.

Under constant atmospheric conditions, the at-sensor TIR radiance R (i) observed in a given spectral channel at pixel i (i = 1, …, n; n is the total number of pixels in the image) can be modeled by the following equation [6]: where τ is the transmissivity of the atmosphere; ε is the emissivity of component k (k = 1, …, K; K is the total number of land cover components); f (i) is the fraction of component k within pixel i; B is the blackbody radiance defined as the integration of Plank function over the spectral channel given T (i), which is the surface temperature of component k within pixel i; and R is atmospheric thermal path radiance. When downscaling is considered, pixel i is disaggregated into m higher spatial resolution pixels, denoted by i, where j = 1, …, m. The downscaled at-sensor TIR radiance R (i) in the given spectral channel at pixel i is modeled by the following equation: where all parameters are defined similarly as in Equation (1).

Spectral mixture analysis attempts to retrieve T (i) directly from multiple equations, each of which corresponds to one spectral channel defined by Equation (1). In contrast, our objective here is not to directly retrieve T (i) but to downscale R (i) to R (i), from which subpixel temperatures T(i) can be retrieved. As such, we estimate the functional relationship between R (i) and f (i) by fitting the model in Equation (1) to observed data at the coarse spatial resolution. To avoid estimating the unknown component temperaturesT (i), we make the assumption of isothermal pixel (i.e. component temperatures T (i) are the same within pixel i) as suggested in [9], which allows us to drop the subscript k inT (i). The isothermal assumption may introduce some bias but it is necessary to reduce the number of unknown variables. An adjustment to this bias will be made later. Consequently, the model used to estimate the functional relationship between R (i) and f (i) becomes where T(i) is the LST at the coarse spatial resolution. Note that τ can be absorbed by ε such that ε̄ = τε the effective emissivity of component k after considering atmospheric effects. Equation (3) can then be rewritten as The least squares estimates of parameters are then obtained using observations over all pixels (i = 1, …, n) in the image at the coarse spatial resolution.

Similarly, the isothermal assumption indicates that T (i) = T(i). Hence, Equation (2) can be written as

Using the functional relationship obtained at the coarse spatial resolution, the downscaled at-sensor radiance R (i) at the higher spatial resolution is estimated by where are the plug-in estimates of parameters in (4); the superscript “0” in denotes the initial estimate. To correct the bias introduced by the isothermal assumption, the initial estimate, , is modified to R̂(i) by adding an additive term which is proportional to :

Statistical downscaling

A statistical approach to downscaling TIR radiance directly estimates the subpixel TIR radiance based on higher spatial resolution ancillary data (f (i), k =1,…,K) without referring to the physical models at different scales. In this paper, we applied a statistical method originally developed for disaggregating zonal census counts [16, 17] to downscale TIR radiance. As detailed in Figure 1, this statistical method obtains the high spatial resolution TIR radiance estimate R̂(i) in two stages: initialization and iteration. At the initialization stage, TIR radiance at high resolution pixel i is set equal to the TIR radiance at its corresponding coarse resolution pixel i. At the l-th iteration, a linear regression model between the TIR radiance estimate obtained from the previous iteration (l-1) and the land cover fractions (f (i), k =1,…, K) is fitted to all the pixels. Note that an intercept term is not used in the regression model in order to avoid co-linearity as the land cover factions sum to one. The model estimate of the TIR radiance is obtained as , where β̂1, …, β̂ are estimates of regression coefficients. To make the average of TIR radiance estimates of m high spatial resolution pixels i equal to the value of their corresponding coarse resolution pixel i (i.e. R(i)), are updated to using the equation: The iteration proceeds until the difference in the regression coefficient of determination r2 between two continuous iterations falls below a threshold value (e.g. 0.001).

Figure 1.

Procedure of statistical downscaling method.

Validation and Evaluation

ASTER data

The TIR radiance used to test our downscaling methods is from Advanced Spaceborne Thermal Emission Reflection Radiometer (ASTER). ASTER is one of the major satellite sensors used for land surface temperature retrieval [18]. It has 14 multi-spectral channels including three visible and near infrared (VNIR) channels with 15-m spatial resolution, six short wavelength infrared channels with 30-m spatial resolution, and five TIR channels with 90-m spatial resolution. ASTER level-1B data acquired on April 25, 2004 were collected in our study area, Yokohama City in Japan. The study area covers most inner portion of the city with most representative land use/land cover (LULC) types from the central business district (CBD) to suburban to rural areas in Japan. A false color display of the 15-m ASTER VNIR imagery of the study area is shown in Figure 2.

Figure 2.

The false color display of the 15-m ASTER VNIR imagery in the study area.

As the five ASTER thermal channels are highly correlated with each other, we only selected the thermal infrared data in channel 13 (10.25–10.95μm) in this study. The ASTER Level-1B data set contains radiometrically calibrated digital numbers (DNs) quantized as 16-bit unsigned integers for the TIR bands. We converted the scaled DNs in channel 13 to at-sensor TIR radiances by the formula [19]: Radiance = (DN value − 1) × 0.005693. The at-sensor TIR radiance image of ASTER channel 13 was shown in Figure 3(a).

Figure 3.

(a) The original 90-m TIR radiances of ASTER channel 13, (b) The simulated 990-m TIR radiances of ASTER channel 13, (c) The estimated 90-m TIR radiances by physical downscaling, and (d) The estimated 90-m TIR radiances by statistical downscaling. Unit for the scale bar is w / (m2 sr μm).

To evaluate the proposed downscaling methods, multi-resolution TIR radiances in a particular spectral region acquired at the same time and location are needed for validation. We simulated multi-resolution thermal data by upscaling the ASTER channel 13 from the original 90-m resolution to 990-m resolution. This allows us to apply the proposed methods to downscale the 990-m thermal data back to 90-m resolution, which will then be validated against the original 90-m thermal data. Specifically, the TIR radiances were upscaled by averaging every 11 × 11 TIR pixels at the 90-m resolution to one pixel at the 990-m resolution [7]. The simulated 990-m ASTER TIR imagery was shown in Figure 3(b). In addition to the simulated TIR radiances, LST at the 990-m resolution is also needed in the physical model. Unlike the TIR radiances, the upscaling function for LST is much more complicated than just a simple average [7]. Therefore, we used the 1000-m MODIS LST products acquired at the same time as the ASTER data over the study area as a good approximate in the physical model.

The three 15-m ASTER VNIR channels shown in Figure 2 were used to generate high resolution ancillary data (i.e. land cover fractions). Specifically the VNIR data were initially clustered into fifteen classes and further grouped into seven general land cover types (water soil two urban types and three vegetation types). The 90-m land cover fractions were then estimated as the proportion of the 15-m pixels of each land cover type within each 90-m pixel. The land cover fractions at 90-m resolution were aggregated into 990-m land cover fractions by simple averaging every 11 × 11 pixels of the 90-m fractions.

Physical downscaling results

The physical model between ASTER at-sensor TIR radiances and land cover fractions at the resolution of 990-m was well fitted to the data as indicated by a very high coefficient of determination r2 = 0.9238. The model parameter estimates in Equation (4) are summarized in Table 1. All parameters are statistically significant with all p-values being nearly close to zeros. The additive component due to atmospheric effects was estimated as 2.56986 w / (m2 sr μm). Note that all the effective emissivities are relatively smaller compared to those published in literature (Sabins, 1986; Lillesand et al., 2004; and the spectral library at UC Santa Barbara, http://www.icess.ucsb.edu). This is due to the fact that the effective emissivities were scaled by absorbing the transmissivity of the atmosphere into the true emissivities.

Table 1.

The parameter estimates of the physical model.

Parameters	Estimate	Standard Error	P-value
R^A	2.56986	0.18414	<2e-16
ε¯^_1	0.49519	0.03397	<2e-16
ε¯^_2	0.70226	0.01831	<2e-16
ε¯^_3	0.71396	0.01829	<2e-16
ε¯^_4	0.74050	0.02203	<2e-16
ε¯^_5	0.69724	0.01876	<2e-16
ε¯^_6	0.63358	0.02272	<2e-16
ε¯^_7	0.64320	0.02016	<2e-16

When the model parameter estimates obtained at the resolution of 990-m were applied to land cover fractions at the spatial resolution of 90-m, the 90-m TIR radiances were obtained. After bias correction (see Equation (7)), the final model estimates of the 90-m TIR radiances were shown in Figure 3(c). A visual comparison revealed that (1) the general patterns of the 90-m TIR radiances in Figure 3(a) were well captured in Figure 3(c); and (2) considerable spatial details in Figure 3(c), which were originally unresolvable in Figure 3(b), appeared in Figure 3(c). However, some hotspots (pixels with high TIR radiances) depicted in red in Figure 3(a) were not captured in Figure 3(c). This is due to the fact that there is a lack of individual patches (hotspots) comprised by more than 11 × 11 pixels of the 90-m resolution in the original thermal image. Additionally, a small portion of the map in Figure 3(c) illustrated some smoothing and blocking effects, which might be a consequence of the assumption of isothermal pixel. Further comparison between the model estimates and the true values (i.e. the observed 90-m TIR radiances of ASTER channel 13) is illustrated by a scatter plot in Figure 4. The scatter plot showed that most of the points lied along a 1:1 line (dotted line) with only a small portion of points departing far from the 1:1 line, indicating a general good agreement between the physical model estimates and the true values. To evaluate the results quantitatively, a simple linear regression model was fitted between the estimated values and the true values using all the 90-m pixels (n = 75625). The results demonstrated that the estimated values matched very well with the true values as evidenced by (1) a very high coefficient of determination (r2 = 0.777), (2) the closeness of the fitted line (solid line) and the 1:1 line in the scatter plot, and (3) a small residual standard error (0.2831).

Figure 4.

Estimated 90-m TIR radiances of ASTER channel 13 by physical downscaling. Solid lines are the fitted models and dotted lines are 1:1 lines. Units for both axis are w / (m2 sr μm).

Statistical downscaling results

The results of the 90-m TIR radiances by statistical downscaling were shown in Figure 3(d). The estimated values well captured the general patterns of the true values in Figure 3(a) and resolved substantial spatial details unseen in Figure 3(b). Compared to Figure 3(c), Figure 3(d) revealed more spatial details and illustrated less smoothing and blocking effects. Similar to Figure 3(c), some hotspots depicted in red in Figure 3(a) were not captured in Figure 3(d). This could be explained by the smoothing effect of the regression model used by the statistical downscaling. The scatter plot in Figure 5 showed almost the same patterns as that in Figure 4 and indicated a good agreement between statistical estimates and the true values. The results of quantitative evaluations were similar to those obtained in section 3.2. Compared to the physical downscaling, the statistical downscaling achieved a slightly higher coefficient of determination (r2 = 0.794) and smaller residual standard error (0.2723). These results further demonstrated the good performance of the statistical downscaling method.

Figure 5.

Estimated 90-m TIR radiances of ASTER channel 13 by statistical downscaling. Solid lines are the fitted models and dotted lines are 1:1 lines. Units for both axis are w / (m2 sr μm).

Conclusions

Subpixel LST retrieval is often needed in various environmental and ecological studies because the LST products retrieved from current satellite thermal sensors have limited spatial resolutions for finer-scale studies and may contain mixed pixels of multiple anisothermal objects in heterogeneous areas. In this paper, we developed two methods for downscaling coarse resolution TIR radiances in preparing for subsequent subpixel temperature retrieval. The first method was developed on the basis of a scale-invariant physical model on TIR radiances. From the physical model, a functional relationship between TIR radiances and land cover fractions was estimated using data at a coarse spatial resolution. The downscaled TIR radiances were then estimated by applying this functional relationship to land cover fractions at a high spatial resolution. The second method was based on a statistical relationship between TIR radiances and land cover fractions at a high spatial resolution. In this statistical downscaling method, high spatial resolution TIR radiances were initialized by the coarse spatial resolution observations and then iteratively regressed to the high spatial resolution land cover fractions until no significant improvements between two continuous iterations were achieved.

The two downscaling methods were applied to simulated 990-m TIR radiances of ASTER channel 13. The estimated 90-m TIR radiances were then validated against the original 90-m TIR radiances. The visual comparison of the results revealed that both downscaling methods successfully captured the general patterns of the original data and resolved considerable spatial details. The quantitative assessments indicated a strong agreement between the true values and the estimated values generated by both methods. Future research could consider the use of spatial dependence in the downscaling methods and explore other ancillary data. In conclusion, the downscaling methods developed in this paper showed promising results for further subpixel land surface temperature retrieval.