Land Surface Temperature Retrieval from Landsat 8 TIRS — Comparison between Radiative Transfer Equation-Based Method , Split Window Algorithm and Single Channel Method

Accurate inversion of land surface geo/biophysical variables from remote sensing data for earth observation applications is an essential and challenging topic for the global change research. Land surface temperature (LST) is one of the key parameters in the physics of earth surface processes from local to global scales. The importance of LST is being increasingly recognized and there is a strong interest in developing methodologies to measure LST from the space. Landsat 8 Thermal Infrared Sensor (TIRS) is the newest thermal infrared sensor for the Landsat project, providing two adjacent thermal bands, which has a great benefit for the LST inversion. In this paper, we compared three different approaches for LST inversion from TIRS, including the radiative transfer equation-based method, the split-window algorithm and the single channel method. Four selected energy balance monitoring sites from the Surface Radiation Budget Network (SURFRAD) were used for validation, combining with the MODIS 8 day emissivity product. For the investigated sites and scenes, results show that the LST inverted from the radiative transfer equation-based method using band 10 has the highest accuracy with RMSE lower than 1 K, while the SW algorithm has moderate accuracy and the SC method has the lowest accuracy. OPEN ACCESS Remote Sens. 2014, 6 9830


Introduction
Land surface temperature (LST) is a key parameter in the physics of the earth surface through the process of energy and water exchange with the atmosphere, which plays an important role in a wide variety of scientific studies, such as ecology, hydrology, and global change studies [1,2].Thermal infrared (TIR) remote sensing provides a unique method for obtaining LST information at the regional and global scales since most of the energy detected by the sensor in this spectral region is directly emitted by the land surface [3].Many efforts have been devoted to establish methods for retrieving the LST from remote sensing data, and significant progresses have been made over the past decade [4].
The Landsat project provides a particular opportunity for the LST retrieval [5][6][7], as it has a relatively long data record period, with the launch of Landsat 3 in 1978 [8].From the Multispectral Scanner (MSS) of Landsat 3 to the Thematic Mapper (TM) of Landsat 4 and 5, and following by the Enhanced Thematic Mapper Plus (ETM+) of Landsat 7, there was only one thermal infrared channel available [8,9] (Figure 1).Therefore, a single-channel (SC) algorithm was developed to derive LST from this band.Influential researches were conducted mainly by Jimenez-Munoz et al. [10][11][12] and Qin et al. [13].Accurate determination of the LST using the SC method requires high-quality atmospheric transmittance/radiance code to estimate the atmospheric features involved in the radiative transfer equation [14].Previous TM and ETM+ sensors (Figure 1) have only one thermal band, while the Landsat 8 TIRS has two spectrally adjacent thermal bands.That is suitable for the split-window (SW) algorithm.The SW algorithm uses two thermal bands typically located in the atmospheric window between 10 and 12 μm [15].The basis of the SW algorithm is that the radiance attenuation for atmospheric absorption is proportional to the radiance difference of simultaneous measurements at two different wavelengths, each subject to different amounts of atmospheric absorption [3].It is the most widely used algorithm for LST retrieval due to the simplicity and robustness.Many operational LST products have been generated using different SW algorithms from various sensors, including the Advanced Very High Resolution Radiometer (AVHRR) [16], Advanced Along-Track Scanning Radiometer (AATSR) [17], Moderate Resolution Imaging Spectroradiometer (MODIS) [18], Spinning Enhanced Visible and Infrared Imager (SEVIRI) [19] and Geostationary Operational Environmental Satellites (GOES) [20].A review of the SW algorithm and different published researches can be found in [14,21].
Most recent studies for LST retrieval from Landsat 8 TIRS are by Rozenstein et al. [22] and Jimenez-Munoz et al. [23].The former one derived a SW algorithm and studied the parameters' (including water vapor content and land surface emissivity) sensitivity of the algorithm [22].While the latter introduced the SC method and another general SW algorithm from forward-simulated atmospheric profile databases.In this paper, three different LST retrieval approaches were explored and compared for the TIRS, including the radiative equation-based method, the SW algorithm and the SC method.For the radiative transfer equation-based method, the atmospheric profile was extracted from the NCEP data set and used to simulate atmospheric transmittance, downwelling and upwelling radiance from the MODTRAN (MODerate resolution atmospheric TRANsmission) model.For the SW algorithm, we chose the algorithm developed by Mao et al. [24,25], which only depends on atmospheric transmittance and land surface emissivity.The theory was originally created to estimate LST from Landsat TIR images by a single mono-window algorithm [13], then adapted to a split-window algorithm for MODIS and ASTER images [24,26].In this study, the coefficients of the SW algorithm were re-parameterized, corresponding to the TIRS' spectral response curve.Meanwhile, the atmospheric transmittance was simulated from water vapor content with typical atmospheric profiles by the moderate resolution atmospheric transmission (MODTRAN) model.For the SC method, a novel algorithm from Jimenez-Munoz et al. [12] was used.The land surface emissivity was derived from a NDVI thresholds method [27].Four SURFRAD sites were selected to conduct the validation and methods' comparison.The aims of this study were: (1) testing the suitability of three presented methods for LST estimation from the Landsat 8 TIRS imagery, and (2) comparing and analyzing retrieval results of those methods.Therefore, it showed different procedures to the readers for retrieving LST from Landsat 8 TIRS data in order to contribute more employments of this sensor in the future.

Radiative Transfer Equation and Radiative Transfer Theory Based Method
A simplified radiative transfer equation can express the apparent radiance received by a sensor [25]: where B i (T i ) the radiance received by channel i of the sensor with brightness temperature T i , the detail of the TIRS' calibration can be found at USGS Landsat Project website.B i (T s ) the ground radiance τ i (θ) atmospheric transmittance for channel i when view zenith angle is θ.TIRS is treated as nadir viewing since the view angle is no more than 7.5° [28].
According to Plank's law, B i (T s ) can be expressed as: where T s is the land surface temperature, c is the light speed (c = 2.9979×10 8 m/s), h is the Planck constant (h = 6.6261×10 −34 J•s), k is the Boltzmann constant (k = 1.3806×10 −23 J/K), λ i is the effective band wavelength for band i, which is defined as: ( ) where f i (λ) is the spectral response function for corresponding band.λ 1,i and λ 2,i are the lower and upper boundary of f i (λ). where With the thermal radiance measured at sensor level, accompany with the atmospheric parameters obtained with radiosounding, which can be used to estimate I i ↓ , I i ↑ and τ i from the radiative transfer model, the LST can be retrieved according to Equation (4).However, in this study, we only used the spatial and temporal closest atmospheric profile from the original 1° × 1° NCEP reanalysis data to simulate I i ↓ , I i ↑ , and τ i from the MODTRAN model for TIRS, since radiosonde profiles are not usually available except in dedicated field campaigns or around meteorological radiosounding stations with launching times when satellites passes.The atmospheric profile used in this study was the NCEP final operational global analysis data.This product is from the Global Data Assimilation System (GDAS), which continuously collects observational data from the Global Telecommunications System (GTS) and other sources for numerous analyses.The data is on a longitude/latitude grid and generated globally every 6 h (0:00, 06:00, 12:00 and 18:00 UTC).The extracted atmospheric profiles have 26 mandatory levels from 1000 to 10 hPa [29,30].Other vertical atmospheric parameters include the geopotential height, the air temperature and the relative humidity.We extracted the NCEP profiles from the on-line atmospheric correction tools [31].A recent study proved that the NCEP has enough vertical resolution for the atmospheric profile which is effective for the τ simulation in the radiative transfer model [29].However, there is an argument for the accuracy of NCEP data when it is applied to large and complex terrain areas due to spatial and temporal resolution [32].Comparison between NCEP, AIRS, MOD07 L2 products and the radiosonde measurement for LST retrieval from the single window algorithm revealed that the interpolated NCEP data is fairly accurate for local simulation [33].Consequently, we used the NCEP data to simulate I i ↓ , I i ↑ , and τ i from the MODTRAN model.As for the land surface emissivity, we will discuss it in the following section.

Qin et al. figured that I i
↓ and I i ↑ can be expressed as [13,26]: where B i (T a ) is the effective mean atmospheric radiance with effective mean atmospheric temperature T a .B i (T a

↓
) is the effective dowelling mean atmospheric radiance with effective dowelling mean atmospheric temperature T a ↓ .It was also noted that, the error of estimated LST caused by the difference of T a and T a ↓ is fairly insignificant [26].Therefore, radiance received by two channels (band 10 and 11) of Landsat 8 TIRS can be rewritten into Equations ( 7) and ( 8) [24]: Based on Plank's law, the spectral radiance emitted by an object is a nonlinear function.Thus, we linearized Planck's radiance function through the application of Taylor's expansion, since the simulated radiance received by two bands of the TIRS with spectral response curves (Figure 1) at the range of −10 °C to 50 °C (Figure 2) is close to a linear function.
( )  The Taylor expansion of Planck's function for band 10 and 11 (i) is presented as Equation (10), where T j refers to the brightness temperature T 10 , T 11 for respective bands, T a and T s as the atmospheric temperature and the land surface temperature.The parameter L i is defined in Equation (11).As mentioned before, we simulated the temperature range from −10 °C to 20 °C and 20 °C to 50 °C for calculating the piecewise linear relationship between L i and T. It has been suggested by Rozenstein et al. [22] and Qin et al. [26] and proved to have better performance for the linearization for Planck's function.Results are shown in Table 1.The high coefficient of determination (r 2 ) and low root mean square error (RMSE) indicate that the linearization for L i from T is reasonable in corresponding temperature range.
(10) (11) Overall, the Equations ( 7) and ( 8) can be overwritten as follows: ) L 10 and L 11 , can be calculated from Table 1 within a specific brightness temperature range for band 10 and 11 respectively.T s is hereby calculated from Equations ( 12) and ( 13) with the elimination of ∂B 10 (T 10 )/∂T and ∂B 11 (T 10 )/∂T as following:   Therefore, the unknown variables needed to inverse T S in Equation ( 18) are τ 10 , τ 11 , which are atmospheric transmittances for band 10 and 11 of TIRS; and ε 10 , ε 11 , which are land surface emissivities for respective bands.
Because of many technical difficulties, the atmospheric transmittance is usually not available in situ when satellite passes.Generally, the most practical way to determine the atmospheric transmittance is through the simulation of local atmospheric conditions, especially with water vapor content [13,24,26].In this study, we simulated the relationship between the atmospheric transmittance (τ) and the water vapor content (w) with two typical atmospheric profiles: the 1976 US standard and the Mid-latitude Summer atmospheric profile from MODTRAN.The range of water vapor content is from 0.2 g/cm 2 to 6.0 g/cm 2 , with a 0.2 g/cm 2 interval.Qin et al. [26] showed that it is better to divide the water vapor content range into several sections and evaluate each of them separately in order to achieve a better accuracy, when this relation is applied for a large range of values.Thus, we divided the water vapor range into two ranges: 0.2-3.0g/cm 2 and 3.0-6.0g/cm 2 , and used the quadratic regression for each range to simulate the relationship between τ and w, showed in Table 2.The estimation equations listed in Table 2 have high r 2 and low RMSE, which indicates that the estimation of transmittance with water vapor content by these equations has high accuracy.

Single Channel Method
Although the Landsat 8 TIRS has two spectrally adjacent channels which are suitable for the split window methods, the single channel (SC) method is still applicable for those bands.In this study, we used a general SW method proposed by Jimenez-Munoz et al. [12,34] for LST retrieval.The equation for this general SW method is as [11]: All the symbols have the same meaning in former equations.For Ψ 1 , Ψ 2 , Ψ 3 , they are derived from the water vapor content [11,12]: The parameters η, ξ, χ, φ are related to the effective band wavelength.According to the relative spectral response curves of TIRS (Figure 1), effective band wavelengths are 10.896 μm and 12.006 μm for band 10 and 11 respectively.Jimenez-Munoz et al. [12] suggested that those parameters can be parameterized as a third degree regression with the wavelength, at the wavelength range of 8-12 μm.Therefore, those parameters which we used to estimate Ψ 1 , Ψ 2 , Ψ 3 from the water vapor content (w) are listed in Table 3. Table 3. Parameters for estimation Ψ 1 , Ψ 2 , Ψ 3 (unitless) from water vapor content (w) in SC method.

Band 10
Band 11 The SC method has been applied to Landsat 5 TM, Landsat 7 ETM+, MODIS ASTER and ENVISAT AATSR sensors' TIR bands for LST retrieval [10,12,33,35].It has reasonable accuracy and only needs two parameters for LST inversion: the land surface emissivity (ε) and the water vapor content (w), which is very applicable for handling single channel TIR data.In this study, we used the SC method for LST inversion with Landsat 8 TIRS imagery.

Land Surface Emissivity Estimation
As mentioned in Sections 2.1-2.3, the land surface emissivity (ε) is indispensable for LST inversion.The emissivity of land, unlike that of oceans, can differ significantly from unity and vary with vegetation, surface moisture, roughness, and viewing angles [36].Although a series of simultaneous LST and LSE retrieval methods have been proposed [14,[37][38][39] and proved to be more accurate for satellite-derived LST.In this study, we still used the LST retrieval methods with prior known LSE, since the simultaneous methods need specific requirements and sophisticated algorithms [14,38,39], while the prior known LSE methods appear to be more practical with reasonable accuracy [27,40] for LST retrieval from Landsat imagery [11,41].
Three major methods were proposed for LSE estimation before LST inversion: classification-based emissivity method (CBEM) [42,43], NDVI-based emissivity method (NBEM) [44][45][46] and day/night temperature-independent spectral-indices (TISI) based method [47][48][49].The CBEM obtains the LSE image from a classification image, in which an emissivity value for each class is assumed in advance.However, this is not very operative because we need a good knowledge of the study area and emissivity measurements on the surfaces representatives of different classes coincident with the satellites transiting time [11].Meanwhile, several requirements may limit the usage of the TISI method.First of all, approximate atmospheric corrections and concurrence of both MIR and TIR data are required [50].Then, the surfaces must be observed under similar observation conditions, e.g., viewing angle, during both day and night [51,52].Additionally, accurate image co-registration must be performed [51].Due to the orbit and revisiting cycle, the TISI method is not applicable for Landsat 8 TIRS.An alternative, operative procedure is the NBEM [11].Because of its simplicity, this method has already been applied to various sensors with access to VNIR data [27,46,50,53,54].In this study, a NDVI thresholds method was used for LSE estimation from Landsat 8 imagery as the following Equation: ( ) The emissivities of vegetation (ε v ) and soil (ε s ) were calculated from the MODIS UCSB (University of California, Santa Barbara) emissivity library, using the following equation [55]: where ε i is the emissivity for channel i. ε i (λ) is the spectral emissivity.Other symbols have the same meaning with previous Equations.Soil and vegetation emissivities for Landsat 8 TIRS are listed in Table 4.The vegetation fraction (P v ) is derived from NDVI.For Landsat 8 imagery, it is calculated from red and near infrared bands (band 4 and 5) from the Operational Land Imager (OLI):  (28) where NDVI min = 0.2, NDVI max = 0.5, ρ 5 and ρ 4 are land surface reflectance after the atmospheric correction.
C in Equation ( 25) is a term which takes the cavity effect into account due to the surface roughness (C = 0 for flat surfaces).Sobrino et al. [50] suggested that C i can be estimated as the following: F' is the geometrical factor ranging between 0 and 1, depending on the geometrical distribution of the surface [27,40], which is typical 0.55 [56].

Validation Sites and Landsat 8 Imagery Processing
Four Surface Radiation budget network (SURFRAD) sites operated by the National Oceanic and Atmospheric Administration (NOAA) were selected for algorithm validation, corresponding to 41 scenes of Landsat8 imagery (Table 5).SURFRAD is the first network to operate across the United States.It began in 1995 with four stations and expanded to six in 1998 [58].The primary objective is to support climate research with accurate, continuous, long-term measurements of the surface radiation budget over the United States [59].Primary measurements in each SURFRAD site are the downwelling and upwelling components of broadband solar and thermal infrared irradiance.Ancillary observations include direct, diffuse solar and photosynthetically active radiation, UVB, spectral solar, and meteorological parameters [59][60][61].All the data are downloaded, quality controlled, and processed into daily files that are distributed freely to the public in near real time by anonymous FTP and the WWW server.At four validation sites, high-quality in situ measurements of upwelling and downwelling long wave radiations are provided [59].In this study, SURFRAD observations were used for the evaluation of LST retrieval.The surface skin temperature T s can be estimated from the following equation [62]: The F ↑ and F ↓ are upwelling and dowelling thermal infrared (3-50 μm) irradiance at the time when the Landsat 8 transits the sites.The ground measurement is 1-min interval, and the exact imaging time was accessed from the metadata for each scene.In the formula, σ is the Stefan-Boltzmann constant , and ε b is the broadband emissivity, which is converted from the two spectral (ε 31 , ε 32 ) emissivities of the 8-day MODIS Land Surface Temperature and Emissivity product (MOD11A2).Emissivities were extracted from the nearest date corresponding to the imaging date in the product for each validation site.We applied 3 × 3 pixels quality mask (emissivity error flag: 00, indicating average emissivity error ≤0.01) with the center of site's location.Pixels with good quality were averaged as the site's emissivity.The conversion from narrow bands emissivities to broad band emissivity is as following [63]: The water vapor content which used in the SW algorithm and SC method is derived from site's meteorology observations.The saturated water vapor pressure (e w * ) is calculated from dry air (T) temperature and air pressure (P) [ The unit of e is hectopascal and the unit of w is g/cm 2 , the convert factor is 0.098.T, P and RH are meteorology observations from the SURFRAD site.Forty-one scenes of Landsat 8 imagery were downloaded from the USGS EarthExplorer Website.We ordered the reprocessed imagery after 3 February 2014 to ensure that the original data is accurate enough for the LST inversion, as it has been announced that the offsets are removed about 2.1 K from Band 10 and about 4.4 K from Band 11, relative to products processed prior to 3 February 2014.The conversion from digital number (DN) to radiance, reflectance and at-satellite brightness temperature followed the guideline at the USGS website.Calibration parameters were directly accessed from the metadata file.
All images were atmospherically corrected with the ACTOR module in PCI Geomatica 2013 SP3.A cloud mask was generated based on the cirrus band (Band 9) and quality assessment (QA) band to avoid the disturbance of cloud.The geometric correction was performed by the same software, using a 1″/3 national elevation dataset provided by USGS.

Results from the Radiative Transfer Equation Based Method
Ground based LST (LST _g ) are shown in Table 6.LST inverted from Landsat 8 TIRS band 10 and 11 (LST _RT_b10 and LST _RT_b11 ), using the radiative transfer equation-based method are shown in Tables 7  and 8.
The bias (difference between estimated LST and ground LST), SD (standard deviation of the bias) and RMSE (root mean square error for estimated LST and ground LST) are also shown for each site.For both band 10 and 11 inverted LST, there is a general positive bias at selected sites, except for the Goodwin Greek site.For band 10 and 11, biases are 0.06 K and 0.05 K with all samples, indicating a slight overestimation.For band 10, RMSE at four sites are 0.87 K, 1.01 K, 0.93 K and 0.57 K, while for band 11, RMSE are 1.17 K, 1.19 K, 1.12 K and 0.75 K.The accuracy for LST_RT_b10 is higher than LST_RT_b11.
As it has been announced that the calibration variability of band 10 is 0.12 W/m 2 /sr/μm (~0.8 K) and 0.2 W/m 2 /sr/μm (~1.75 K) for band 11 by the USGS Landsat website.Although we used the reprocessed Landsat 8 data in this study, there is still large calibration uncertainty associated with band 11, which leads to the error for LST estimation.Meanwhile, this is expected, since band 11 is more affected by the water vapor continuum absorption and thus more sensitive to errors in atmospheric profiles [33].Moreover, the accuracy of the radiative transfer equation-based method depends mostly on how well the atmospheric profiles are able to represent the actual atmosphere over the site when the sensor captures image.In this study, the NCEP reanalysis data was used as the atmospheric profile, since the lack of the real time radiosonde data for each scene of Landsat 8 image.Wan and Li [65] proposed a method to evaluate the quality of atmospheric profiles by analyzing the difference between the LSTs derived from instruments which have split window bands at 11 and 12μm, such as MODIS and AATSR.Coll et al. [33] compared different sources of atmospheric profiles for land surface temperature retrieval from single channel thermal infrared data, including radiosonde, NCEP, Aqua AIRS and MOD07L2, and suggested that NCEP reanalysis profiles provide an alternative to radiosonde data.In this study, we compared the LST estimated from TIRS band 10 and 11, the result is shown in Table 9.There is a general overestimation trend for LST _RT_b11 than LST _RT_b10 at Bondville and Sioux Falls, while underestimation at Fort Peck and Goodwin Greek.However, all the difference is less than 1 K.The RMSE are 0.53 K, 0.26 K, 0.34 K and 0.25 K for each site respectively, suggesting that the NCEP is fairly accurate enough for retrieving LST from Landsat 8 TIRS band 10 & 11 based on the radiative transfer equation.

Results from the SW Algorithm
Results for LST estimated from the Landsat 8 TIRS (LST _SW ), using the SW algorithm are show in Table 10.There is a general underestimation for four sites except for the Sioux Falls site, as the bias is 0.07 K for this site.The bias for all sites is −0.15K, suggesting an underestimation for all scenes.
Meanwhile, the standard deviation of the bias for the Goodwin Greek site is 1.18 K, higher than other sites.Landscape around Goodwin Creek is rural pasture land with several ponds surrounded, which was checked from the high resolution image in Google Earth.Variations in the landscape leads to the emissivity estimation errors, since we only considered soil and vegetation for ε in this study.It may also note that the point scale ground measurement is largely incompetent to the pixel area retrieval from satellite [66].
Figure 3 illustrates the comparison of LST_g and inverted LST (LST_ SW ).It shows that they have high correlation (R 2 = 0.989).However, there is an underestimate trend for high LST condition, which can also be found from Table 3.Since the hypothesis behind Mao et al. [24] method is that the difference of T a and T a ↓ is fairly insignificant, it should be tested for high water vapor content and high air temperature conditions.That may also introduce external error for LST inversion.RMSE for Bondville and Fort Pack sites are 0.72 K and 0.73 K, while 1.18 K and 1.15 K for Goodwin Greek and Sioux Falls sites for (Table 10), and 1.025 K for all sites (Figure 3), suggesting that the SW algorithm has the potential for accurate LST inversion from TIRS imagery.

Results from the SC Theory
LST estimated from the Landsat 8 TIRS band 10 and 11 (LST _SC_b10 and LST _SC_b11 ), using single channel method are showed in Tables 11 and 12. Overestimation happens for all sites using both bands (bias > 0).For all samples, biases are 0.44 K with band 10 and 0.73 K with band 11.RMSE for LST _SC_b11 are 1.78 K, 1.43 K, 1.71 K and 1.34 K for each site respectively, which are higher than RMSE of LST _SC_b10 .It has the similar situation with radiative transfer equation based method, indicating that band 11 has larger calibration uncertainty.Moreover, we used the function proposed by Jimenez-Munoz and Sobrino [12] to calculate the parameters used for the SC method.However, it has to be noted that the general spectral response function with a full width half maximum (FWHM) of 1 μm applied to simulate the function may introduce extra errors [12,34].Moreover, Jimenez-Munoz and Sobrino [34] figured out that the most important error source in SC method is due to atmospheric effects, which leads to an error on the LST between 0.2 K and 0.7 K, and the land surface emissivity uncertainty, which leads to an error on the LST between 0.2 K and 0.4 K.The most recent study by Jimenez-Munoz et al. [23] indicates that higher water vapor content will introduce more error into the LST estimation, but can be partly solved by computing the atmospheric functions directly.In our study, this has been consolidated by comparing the result for radiative transfer equation based method and SC method for band 10 of TIRS, since LST estimated from the former has higher accuracy than the latter.As for band 11, the LST estimation accuracy is not acceptable for both methods.It has to be noted that we cannot ensure the effect of errors in LSE on LST inversion in this study, because of the lack of ground LSE measurement.It must be pointed out that, the uncertainty inside the MODIS product we used for the derivation of ground LST will introduce external error for result verification, although many other studies have used it for LST products' validation [62,67,68].

Comparison of Three Methods
RMSE of LST estimated from the radiative transfer equation-based method using band 10 and 11 are 0.903 K and 1.153 K for all 41 samples, compared with 1.39 K and 1.67 K for the SC method.RMSE for SW algorithm is 1.025 K (Figure 3) for all samples.It shows that the LST retrieval from band 11 has more uncertainty than band 10 (Tables 7, 8, 11 and 12).The radiative transfer equation-based method using band 10 has the highest accuracy with RMSE less than 1 K, while the SW algorithm has moderate accuracy.However, the difference is not quite obvious (with RMSE difference of 0.122 K).Specifically, for Good Greek and Sioux Falls sites, the former (RMSE: 0.93 K and 0.57 K) outperforms the latter (RMSE: 1.25 K and 1.10 K) (Tables 7 and 10).For the Bondville site, the latter has lower RMSE of 0.73 K than the former of 1.17 K (Tables 7 and 10).For the Fort Peck site, the difference is not obvious.As mentioned in Section 2.5, the stray light issue introduces disturbance from nearby area to the scene of TIRS, leading to larger uncertainty for LST retrieval.The USGS Landsat website has recommended not using band 11 for SW algorithm.In this study, the highest accuracy for LST retrieval comes from the radiative transfer equation-based method, using band 10 along with NCEP atmospheric profile.The SW algorithm has medium accuracy compared with other two methods.While the SC method has the lowest accuracy.Results for LST retrieval from band 11 show larger uncertainty than band 10, which is consistent with the announced TIRS calibration issue.
Furthermore, we applied the one side analysis of variance (ANOVA) [69] for these methods.The box plot of the estimated LST is shown in Figure 4. F-value and P-value for the ANOVA are 0.09 and 0.9865 respectively, which indicates that the differences between three presented methods are not significant.It consolidates the comparison of RMSE.Moreover, we used t test [70] for results from the radiative transfer equation-based method using band 10 and the SW algorithm.The p-value is 0.0972 (larger than 0.05, less than 0.1), which suggests that results from these two methods are moderately different.
Meanwhile, it should be noted that we derived a more general SW algorithm in this study.The parameter L i was estimated by the piecewise linear fitting (−10-20 °C and 20-50 °C, Table 1), while the atmospheric transmittance was estimated by the piecewise quadratic regression (0.2-3.0 g/cm 2 and 3.0-6.0g/cm 2 under two typical atmospheric profiles, Table 2).These division may not fall into the prevailing range of the validation sites' condition.However, it does not reduce the validity of the algorithm.The errors for L i estimation are shown in Figure 5 and absolute errors are less than 0.1 K under most conditions.As for the atmospheric transmittance (τ), Qin et al. [26] revealed that the even the |T 10 -T 11 | is less than 0.5 K, the split window algorithm can still produce a very accurate LST estimation in spite of a big transmittance error, although Rozenstein et al. [22] suggested that the contribution of τ error to the LST estimation is complex and also depends on the emissivity in both channels.

Conclusions
We applied the radiative transfer equation based method, the split window algorithm and the single channel method to the Landsat 8 TIRS data.For the first method, the NECP data was used to simulate the parameters needed as inputs for MODTRAN model.For the SW algorithm, coefficients were adapted based on spectral response functions of two TIRS bands (band 10 and 11).Atmospheric transmittance was derived from the MODTRAN model, using standard atmospheric profile.For the SC method, parameters were derived from the regression of a general spectral function corresponding to simulated atmospheric absorption profiles.Land surface emissivity were estimated by the NDVI threshold method.Forty-one scenes of imagery were used for validation at four selected SURFRAD sites with high frequency irradiance measurement and MODIS LSE product.For the investigated sites and scenes, results show that LST retrieval from the radiative transfer equation-based method using band 10 has the highest accuracy with RMSE ≤1 K, while the SW algorithm has moderate accuracy and the SC method has the lowest accuracy with all scenes.For those methods using single band, LST estimated from band 10 has higher accuracy than band 11.Future work should focus on theoretical evaluations of the effect for input parameters' (LSE, atmospheric transmittance) errors on estimated LST.Since the limitation of validation sites and scenes, the range of ground temperature in this study is mostly within 10°C to 30°C and the surrounding areas are relative homogeneous.More land surface types and different temporal scenes should be tested to verify the investigated three methods for LST estimation from Landsat 8 TIRS imagery.
9), T is the brightness temperature, other parameters have the same meanings with previous equations.

Figure 2 .
Figure 2. Relationship between temperature and radiance for Thermal Infrared Sensor (TIRS) band 10 and 11.

Figure 3 .
Figure 3.Comparison between SW inverted LST and ground LST.

Figure 5 .
Figure 5. L i estimation error for piecewise linear fitting.

Table 1 .
Linear fitting coefficients for Parameter L i .

Table 2 .
Relationship between atmospheric transmittance (τ) and water vapor content (w) with 1976 US standard and Mid-latitude Summer atmospheric profile.

Table 4 .
Emissivities of soil and vegetation for Landsat 8 TIRS band 10 and 11.

Table 5 .
Four Surface Radiation Budget Network (SURFRAD) sites' information and corresponding imaging date.

Table 6 .
Land Surface Temperature (LST) (LST _g ) from ground measurement at four selected sites for 41 scenes.

Table 7 .
Inverted LST (LST _RT_b10 ) using radiative transfer equation-based method from TIRS band 10 and the difference (Δ) between LST _RT_b10 and LST _g for four selected sites.

Table 8 .
Inverted LST (LST _RT_b11 ) using radiative transfer equation-based method from TIRS band 11 and the difference (Δ) between LST _RT_b11 and LST _g for four selected sites.

Table 9 .
Difference between LST inverted from TIRS band 10 based on radiative transfer equation based method (LST _RT_b10 ) and band 11 (LST _RT_b11 ) for four selected sites.

Table 10 .
Inverted LST (LST _SW ) using a split window (SW) algorithm from TIRS and the difference (Δ) between LST _SW and LST _g for four selected sites.

Table 11 .
Inverted LST (LST _SC_b10 ) using SC theory from TIRS band 10 and the difference (Δ) between LST _SC_b10 and LST _g for four selected sites.

Table 12 .
Inverted LST (LST _SC_b11 ) using SC theory from TIRS band 11 and the difference (Δ) between LST _SC_b11 and LST _g for four selected sites.