A Practical Split-Window Algorithm for Retrieving Land Surface Temperature from Landsat-8 Data and a Case Study of an Urban Area in China

This paper proposes a practical split-window algorithm (SWA) for retrieving land surface temperature (LST) from Landsat-8 Thermal Infrared Sensor (TIRS) data. This SWA has a universal applicability and a set of parameters that can be applied when retrieving LSTs year-round. The atmospheric transmittance and the land surface emissivity (LSE), the essential SWA input parameters, of the Landsat-8 TIRS data are determined in this paper. We also analysed the error sensitivity of these SWA input parameters. The accuracy evaluation of the proposed SWA in this paper was conducted using the software MODTRAN 4.0. The root mean square error (RMSE) of the simulated LST using the mid-latitude summer atmospheric profile is 0.51 K, improving on the result of 0.93 K from Rozenstein (2014). Among the 90 simulated data points, the maximum absolute error is 0.99 °C, and the minimum absolute error is 0.02 °C. Under the Tropical model and 1976 US standard atmospheric conditions, the RMSE of the LST errors are 0.70 K and 0.63 K, respectively. The accuracy results indicate that the SWA provides an LST retrieval method that features not only high accuracy but also a certain universality. Additionally, OPEN ACCESS Remote Sens. 2015, 7 4372 the SWA was applied to retrieve the LST of an urban area using two Landsat-8 images. The SWA presented in this paper should promote the application of Landsat-8 data in the study of environmental evolution.


Introduction
Land surface temperature (LST) plays a significant role in environmental evolution, including heat energy circulation and the water cycle.Retrieving LST from remote sensing data is an important capability for researching the spatial evolution of the regulation of the heat environment at local and global scales [1][2][3][4].Studies on the derivation of LST using satellites have been conducted primarily using NOAA AVHRR data [5,6] at a regional scale.Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) thermal infrared (TIR) data have been used in local-scale studies of LST distributions [7][8][9].Data from medium-resolution thermal infrared sensors, such as MODIS and AVHRR with daily temporal resolution, and high-resolution thermal infrared sensors, such as the Landsat series that pass once every 16 days, have been supplied to the public in recent decades [10].Many studies on the retrieval of LSTs from satellite TIR data have been performed since the 1990s.The methods used can be broadly classified into three categories: single-channel methods, multi-channel methods, and multi-angle methods [11].The single-channel method, also called the mono-window algorithm, utilizes the radiance measured by a sensor in a single channel to retrieve LST using the general radiance transfer equation [12][13][14].The multi-channel method, also called the split-window algorithm (SWA), uses the different absorptions of two TIR channels, linearizing or nonlinearizing the radiance transfer equation with respect to the temperature or wavelength.Many efforts have been made to extend the SWA because these algorithms assume that the land surface emissivity (LSE) values in both TIR channels are known [15][16][17][18][19][20][21][22][23][24][25][26][27].The multi-angle method is in accordance with different atmospheric absorptions because of the differential path-lengths associated with observing the same object from various viewing angles [21,24,[28][29][30][31][32].The launch of the Landsat-8 satellite on February 11, 2013 adds to the remote sensing data of the Landsat family, TM and ETM+.The Landsat-8 satellite carries the Operational Land Imager (OLI) with nine bands and Thermal Infrared Sensor (TIRS) with band10 and band11.Using an SWA, the TIRS band10 and band11 provide the atmospheric rectification for the thermal infrared data [33,34].Jimenez-Munoz and Sobrino et al. [35], Rozenstein et al. [10] and Du et al. [36] proposed three different SWAs for LST retrieval from Landsat-8 TIR data.The objective of this paper is to provide another SWA for retrieving LSTs from Landsat-8 data.We attempt to improve the convenience and accuracy of the LST retrieval of Landsat-8 TIR data and provide an alternative.Although the idea of the SWA presented in this paper is based on the practical split-window algorithm proposed by Mao et al. [27], we did not merely transplant Mao's principle; we also studied the input parameter determination and algorithm structures.We conclude that the SWA idea proposed in this study can be not only applied to Landsat-8 TIR data but can also be used with other kinds of remotely sensed thermal infrared data, such as ASTER and MODIS.
In this paper, we have developed a practical SWA for Landsat-8 TIRS data from another perspective and subjected it to a sensitivity analysis and accuracy assessment.Additionally, a specific calculated case for an urban area in China is provided.Specifically, this paper will (1) describe the SWA model and derive the LST calculation; (2) determine the atmospheric transmittance and the emissivity of the ground for the Landsat-8 TIRS; (3) analyse the sensitivity of the essential input parameters and assess the accuracy via simulated data using MODTRAN 4.0; and (4) retrieve the LST of an urban area in China.

Theoretical Derivation of the SWA
The theoretical basis of the LST retrieval algorithm is Planck's law, which is based on the thermal radiance of the ground and the heat transfer from the ground through the atmosphere to the remote sensor [27].According to Planck's law, the thermal emittance from an object can be expressed as Planck's radiance function [37], as follows in Equation (1): where Bn(T) denotes the object's spectral radiance, T its absolute temperature, k the Boltzmann constant, h the Planck constant, and c the speed of light.The units of Bn(T) are Wm −2 sr −1 µm −1 .The general radiance transfer equation [38] for the remote sensing of LST can be formulated as follows in Equation ( 2): where Ts is the LST, Ti is the brightness temperature in channel i, εi is the ground emissivity and τi is the atmospheric transmittance in band i. Bi(Ts) is the ground radiance, and Ii ↓ and Ii ↑ are the downward and upward path radiances, respectively.Ii ↓ and Ii ↑ can be expressed by Equations ( 3) and (4), respectively [27]: (1 ) ( ) (1 ) ( ) where a T  is the average temperature of the downward radiance of the atmosphere, and Ta is the average temperature of the upward radiance of atmosphere.Every term of the radiance transfer equation includes the Planck function.According to Qin's [39] analysis, few differences result from using Ta instead of a T  .Therefore, Equation (2) can be written as follows: which can be consolidated as follows: In Equation ( 6), the value of Bi(Ti) can be calculated by substituting the brightness temperature i T into the Planck function.The terms Bi(Ts) and Bi(Ta) need to be simplified.In this paper, considering the scope of the application and the desired retrieval accuracy, we applied different simplified fitting methods for Bi(Ts) and Bi(Ta).Because Bi(Ts) has a greater impact on the retrieval precision and it contains inversion parameters, we decided to nonlinearly fit it with a quadratic function, whereas Bi(Ta), which has less impact on the retrieval accuracy, can be linearly fit.Then, Equation ( 6) can be written as follows: If the equation is condensed by setting 2 3 (1 )(1 ( 1) ) (1 )(1 ( 1) ) ( ) Then, Equation ( 7) can be expressed as follows: (8) For Landsat-8 TIRS band10 and band11, Equation (8) generates the following equation set:

A T B T C T D A T B T C T D
or, upon consolidation, When this equation is solved for Ts, the result can be expressed as which can be generalized to: where , , ,    and  are parameters computed by 10 10 11 , , ,    and 11  , respectively.

Determination of the Fitting Parameters for ( )
This study simplified Bi(Ts) via a quadratic fitting method and Bi(Ta) via a linear simplification approach.To strengthen the universality of the model, the temperature range is set from 180 K-363 K (−90 °C-90 °C) for Landsat-8 thermal infrared bands 10 and 11.The selection of the temperature range is based on a full consideration of the existence of extreme surface temperature regions.The extreme minimum temperature was reported to have reached −71 °C in Oymyakon village in Russia, and the surface temperature reached as high as 82.3 °C in several desert areas in the summer.Based on the TIR spectral response function of Landsat-8, the effective wavelengths [40] and wavenumbers of these TIRS are calculated to be 10.9034 μm and 917.1417608 cm −1 for the minimum temperature, respectively, and 12.0028 μm and 833.1387464 cm −1 for the maximum temperature, respectively.
Then, the thermal radiation values can be obtained using Equation (1).The corresponding Planck curve is shown in Figure 1.The quadratic fitting function is shown in Table 1, the linear regression formula is listed in Table 2, and the fitting parameters Bi(Ts) and Bi(Ta) are listed in Table 3.

Determination of Atmospheric Transmittance
The atmospheric transmittance is an important factor in the split-window algorithm.According to previous studies, the atmospheric transmittance strongly depends on the atmospheric water vapour content, especially for thermal bands [27].This paper studied the relationship between the Landsat-8 TIRS' atmospheric transmittance estimates and the water vapour in the atmosphere via MORTRAN 4.0 simulations that use the Mid-Latitude Summer atmosphere profile.The simulation data are listed in Table 4. Cubic polynomial fits are adopted for the relationship between the atmospheric water vapour content and transmittance, and the fitting expressions are shown in Table 5.The agreement is very satisfactory.Because the errors in the water vapour content are 0.20 g•cm −2 [41] and 0.10 g•cm −2 , the corresponding atmospheric transmittance evaluation errors are less than 0.031 and 0.016, respectively.The water vapour content w was estimated using the same daytime MODIS data.Although the times of MODIS and Landsat-8 cannot be strictly the same, the Terra satellite carrying MODIS transits the study area approximately 43 minutes earlier than Landsat-8 ((approximately 02:35 versus approximately 03:08).The study area is ringed by mountains on three sides, and we have therefore assumed that the water vapour content w is constant during this time.Kaufmann, Y. J. et al. [42] provided a method to compute the water vapour content in the atmosphere using the reflectance of band2 and band19 of MODIS; thus, we can obtain the water vapour content of the atmosphere using Equation (13).
where w is the water vapour content in g•cm −2 , 0.020 , and 19 r and 2 r denote the reflectance of MODIS band19 and band2, respectively.

Determination of the Emissivity of the Ground
In consideration of the high signal-to-noise ratio (SNR), this study employed the NDVI threshold method, which considers the differing impacts of distinct land cover types (e.g., water, vegetation, bare soil and impervious surfaces), proposed by Sobrino et al. [43] to determine the land surface emissivity (LSE).
We have also referred to the ASTER spectral database [44] to determine the emissivity values of water, vegetation and non-vegetation for Landsat-8 TIRS band10 and band11 (Table 6).Water pixels with negative NDVI values (NDVI < 0) were set to LSE values of 0.991 and 0.986 for band10 and band11, respectively.Fully vegetated pixels with a NDVI value greater than 0.50 were set to LSE values of 0.984 and 0.980.Non-water pixels with NDVI less than 0.20 are considered to be non-vegetated areas, and the LSE values of these pixels were set to 0.964 and 0.970.Finally, pixels with NDVI values greater than 0.20 and less than 0.50 are considered to be mixed pixels incorporating various degrees of vegetation and non-vegetation.The LSE values of these pixels were calculated using Equations ( 14) and ( 15) [45][46][47]. ( where εimix is the emissivity of the mixed pixels for band10 and band11, εiv is the emissivity of the fully vegetated pixels for band10 and band11, εin is the emissivity of the non-water and non-vegetated pixels for band10 and band11, and F = 0.55 is a shape factor considering the geometrical distribution [48].The parameter pv is the scaled NDVI value, in which NDVImin and NDVImax are the NDVI values for non-vegetated and fully vegetated land covers, respectively.Please note that, in recent years, colourful steel plates mainly composed of electro-galvanized steels have been installed in architectural fields, with blue-coloured plates being especially common.As a result, regions with blue roof plates in our study area are relatively prevalent.Because the corresponding NDVI of these blue colour steels is relatively large, generally above 0.20 and sometimes even as high as 0.53, we separated these areas from the 0.20-0.50range of NDVI to produce amplitude ratio values.According to the ASTER spectral database, the Roofing Materials--Galvanized Steel Metal category features wavelengths of 10.9009 μm and 12.0108 μm and corresponding reflectances of 0.041 and 0.038.Therefore, we calculated the Galvanized Steel Metal emissivity values in band10 and band11 to be 0.959 and 0.962, respectively.

Steadiness Analysis of the SWA
To assess the steadiness of the SWA proposed in the paper, we analysed the impact of misestimated SWA input parameters on the LST derivation error.The atmospheric transmittance calculated through the atmosphere's water vapour content and the LSE play a major role in the SWA; thus, the influences of these two input parameters were analysed.To enable comparison with the existing SWA proposed by Rozenstein et al. (2014), the parameter values were kept almost consistent with that of Rozenstein et al. (2014).

Steadiness Analysis of Atmospheric Transmittance
The atmospheric transmittance is estimated from the water vapour content of the atmosphere.The impact of misestimating the input parameter of the atmosphere transmittance can be assessed by the water vapour content in the column, which is taken to be an important input parameter for the SWA.Because the LST error is related to the brightness temperatures of band10 and band11, we set T10 = 283.0-333.0K, with an interval of 10.0 K, for a T10-T11 difference of −3-3 K to calculate the LST estimation error, assuming e10 = 0.967 and e11 = 0.971 and undervaluing the atmospheric water vapour content by 0.1 and 0.2 g•cm −2 , as was done by Rozenstein et al. (2014) The water vapour content range was 1.0-4.0g•cm −2 .
The LST estimation error is large when the atmospheric water vapour content is 4.0 g•cm −2 and T10 is 3 K lower than T11.In addition, for constant values of the other parameters, higher T10 values correlate with higher LST errors.The specific error variation is shown in Figure 2A,B.In the case of underestimating the water vapour content by 0.1 g•cm −2 and 0.2 g•cm −2 , the corresponding maximum LST errors are 0.56 K and 1.11 K, respectively, and the RMSEs of the LST error are 0.30 K and 0.59 K, respectively.The contribution of the calculated atmospheric transmittance error to the retrieved LST is not only related to the water vapour content but also correlates with the LSEs of band10 and band11.The ratio of the LSEs of band10 and band11 is invariant in the analysis process.When the water vapour content column is undervalued by 0.1 and 0.2 g•cm −2 , the corresponding maximum LST errors are 0.21 K and 0.41 K, respectively, and the RMSEs of the LST error are 0.10 K and 0.19 K, respectively, as illustrated in Figure 3A,B.

Stability Analysis of LSE
Due to the possibility of having separate LSE errors in both bands, the impact on the LST estimation error of simultaneous errors in the LSE values of both band10 and band11 is discussed in the following.The LSE estimation error also generates an LST calculation error when the atmospheric transmittance has been precisely estimated.In general, the LST calculation error has a positive correlation with the T10 and T10-T11 values and has a negative correlation with LSE, as illustrated in Figure 4A,B.To be specific, the maximum LST calculation error, 0.44 K, occurs when T10 is 330 K, the T10-T11 value is 3 K and the LSE is 0.900 due to undervaluing by 0.005.In the case of underestimating the LSE by 0.001, the maximum LST calculation error is 0.09 K.In Figure 4A,B, the atmospheric water vapour content was kept constant at 1.5 g•cm −2 .The sensitivity of the LSE to the LST estimation error was analysed as the atmospheric water vapour content varied from 1.0-4.0g•cm −2 , as seen in Figure 5A,B.The LST calculation error has a positive correlation with the T10 value and has a negative correlation with LSE and water vapour content.The maximum LST calculation error, 0.45 K, occurs when T10 is 333 K, the water vapour content is 1.0 g•cm −2 , and LSE is 0.900 due to undervaluing by 0.005.
Generally speaking, the SWA presented in this paper is stable with respect to the impact of incorrectly estimating the LSE and the atmospheric transmittance, i.e., the atmospheric water vapour content.Furthermore, the LST evaluation error is negatively correlated with the LSE estimate.

Evaluation of the Accuracy of the Proposed SWA
MODTRAN 4.0 was used to simulate a series of thermal radiation processes from the ground to the sensor based on the Mid-Latitude Summer atmospheric profile with specified LST, LSE and water vapour content values to evaluate the accuracy of the SWA.First, the simulated thermal radiance values were converted to brightness temperatures corresponding to the Landsat-8 TIRS, band10 and band11, as Ti input parameters for the SWA.Then, according to our SWA and the given input parameters, the LST was calculated.Finally, the LST evaluation errors were obtained by calculating the difference between the given and calculated LSTs.The 90 combination scenes are listed in Table 6.The water vapour content values were 1.0, 2.0, or 3.0 g•cm −2 , the range of the assumed LST values was 10-60 °C, with an interval of 10.0 °C, and the range of the LSE values was 0.940-0.980,with an interval of 0.010.The accuracy assessment results, shown in Table 7, are satisfactory, and the RMSE of the LST estimation errors is 0.51 °C, outperforming the result of 0.93 °C from Rozenstein et al. (2014).To assess the general applicability of this SWA to different atmospheric conditions, we also evaluated the accuracy of this SWA under a tropical model and 1976 US standard atmospheric conditions.The assessment results are listed in Tables 8 and 9.The RMSEs of the LST errors under the tropical model and the 1976 US standard atmospheric conditions are 0.70 °C and 0.63 °C, respectively.

Study Area
The city of Taiyuan-Yuci, China, was chosen as the urban study area.This urban area is located in north-central China, between 37°24'-38°04'N and 112°09'-113°06'E.Taiyuan-Yuci covers a land area of 3336 km 2 and has approximately 4.41 million inhabitants (China, the data of the Sixth Population Census, 2010).The region features a warm continental monsoon climate.

Data Used
Two cloud-free Landsat-8 images, acquired at UTC 03:08:38 on June 27, 2013 and UTC 03:06:26 on June 30, 2014, were rectified to a Universal Transverse Mercator (UTM) projection with WGS-84 datum and UTM zone 49.Two daytime MODIS-Terra images acquired at 02:35 and 04:10 on June 27, 2013 and June 30, 2014, respectively, are also used in this paper to calculate the water vapour content of the atmosphere.Because the MODIS and Landsat-8 acquisition times on the same day were not rigidly consistent, we used the average value of the two temporal MODIS water vapour content estimates to calculate the water vapour content of the atmosphere when the Landsat-8 satellite passed over the study area.

LST Retrieval
Using the SWA proposed in this paper, the LST of the urban area was retrieved via the workflow shown in Figure 6.Three sets of input parameters, 10   can be confirmed with the NDVI images and Equations ( 14) and ( 15), and the emissivity values are listed in Table 6.The three input parameters have now all been confirmed and entered into the SWA formula, and the LST images have been produced.Figure 7 shows the LST estimated from the Landsat-8 TIR data.The LST on June 27, 2013, ranges from 296.6 K-323.8K, with a mean value of 304.5 K and a standard deviation of 2.953 K.The LST on June 30, 2014, ranges from 292.4 K-322.9K, with a mean value of 303.7 K and a standard deviation of 3.490 K.

Conclusions
A practical SWA for Landsat-8 TIRS data was proposed in this paper.The analysis of sensitivity shows that the SWA has good stability.In five sensitivity analysis scenarios using the Mid-Latitude Summer atmospheric profile, the minimum LST error is −0.839K, and the maximum LST error is 0.786 K.In general, the LST error is very sensitive when the water vapour content and T10-T11 change together, as is seen in Figure 2B, especially when the water vapour content is 0.5-1.0gcm −2 .Based on the validation results of the simulation data generated by MODTRAN 4.0, the SWA performs very well.The verification result is satisfactory in that the maximum and minimum LST errors are 0.993 K and 0.016 K, respectively, and the LST error RMSE is 0.51 K under Mid-Latitude Summer atmospheric conditions.The LST error RMSEs under Tropical and 1976 US standard atmospheric conditions are 0.70 K and 0.63 K, respectively.The assessment results demonstrate that the SWA is not only steady and accurate but also universal for different atmospheric conditions.
In this paper, we offer a strategy for evaluating the atmospheric water vapour content using MODIS data and for determining the relationship between the atmospheric transmissivity and water vapour content from Landsat-8 TIRS data using simulation data from MODTRAN 4.An exploratory attempt was performed on a corresponding LSE estimation based on Landsat-8 TIRS data, presented by the LSE of three classic earth surfaces and one new earth surface (galvanized coloured roofs), which further improved the SWA's integrality and practicality.Additionally, a complete application was performed by retrieving the LST of a study area.The proposed SWA in this paper will contribute to the completion of the inversion algorithm system of LST and will promote the application of remote sensing technology in the research of the thermal environment of the Earth.

Figure 2 .
Figure 2. (A) The LST evaluated error (K) in the case of undervaluing the atmospheric water vapour content by 0.1 g•cm −2 for a water vapour content in the range of 1.0-4.0g•cm −2 and an interval of 1.0 g•cm −2 over six gradations of T10 = 283.0-333.0K and a T10-T11 value range of −3.0-3.0K.e10 = 0.967, and e11 = 0.971.(B) The LST evaluation error (K) for undervaluing the atmospheric water vapour content by 0.2 g•cm −2 under the same conditions as in (A).

Figure 3 .
Figure 3.The LST evaluation error (K) in the case of undervaluing the atmospheric water vapour content by 0.1 g•cm −2 (A) and 0.2 g•cm −2 (B) for water vapour contents in the range of 1.0-4.0g•cm −2 at an interval of 1.0 g•cm −2 over six gradations of T10 = 283.0-333.0K and a T10-T11 value of 1 K.The LSE variable is independent.

Figure 4 .Figure 5 .
Figure 4.The illustration of the LST calculation error (K) in the case of undervaluing LSE by 0.001 (A) and 0.005 (B) over four gradations of T10 = 270.0-330.0K with an interval of 20 K.The atmospheric water vapour content was constant at 1.5 g•cm −2 , and the T10-T11 values were 1.0, 2.0, and 3.0 K.The LSE variable is independent.

Figure 7 .
Figure 7.The LST Spatial distribution of the study urban area retrieved from Landsat-8: (A) on June 27, 2013, and (B) on June 30, 2014.

Table 1 .
The quadratic fitting results for the Landsat-8 TIRS.

Table 2 .
The linear regression results for the Landsat-8 TIRS.

Table 3 .
The fitting parameters of Bi(Ts) and Bi(Ta) for the Landsat-8 TIRS.

Table 4 .
The water vapour of the atmosphere and the corresponding atmospheric transmittance as simulated by MORTRAN 4.0.

Table 5 .
The cubic fit regression functions of the atmospheric transmittance and water vapour content range of 0.5-3.0g•cm −2 .

Table 6 .
The emissivity values of water, vegetation and non-vegetation for Landsat-8 TIRS band10, and band11.

Table 7 .
Estimation errors of LST for different simulated combinations of specified LST, LSE, and water vapour content values based on the Mid-Latitude Summer atmospheric profile.

Table 8 .
Estimation errors of LST for different simulated combinations of specified LST, LSE, and water vapour content values based on the Tropical model atmospheric profile, RMSE = 0.70 °C.

Table 9 .
Estimation errors of LST for different simulated combinations of specified LST, LSE, and water vapour content values based on the 1976 US Standard atmospheric profile, RMSE = 0.63 °C.