A Practical Split-Window Algorithm for Estimating Land Surface Temperature from Landsat 8 Data

This paper developed a practical split-window (SW) algorithm to estimate land surface temperature (LST) from Thermal Infrared Sensor (TIRS) aboard Landsat 8. The coefficients of the SW algorithm were determined based on atmospheric water vapor sub-ranges, which were obtained through a modified split-window covariance–variance ratio method. The channel emissivities were acquired from newly released global land cover products at 30 m and from a fraction of the vegetation cover calculated from visible and near-infrared images aboard Landsat 8. Simulation results showed that the new algorithm can obtain LST with an accuracy of better than 1.0 K. The model consistency to the noise of the brightness temperature, emissivity and water vapor was conducted, which indicated the robustness of the new algorithm in LST retrieval. Furthermore, based on comparisons, the new algorithm performed better than the existing algorithms in retrieving LST from TIRS data. Finally, the SW algorithm was proven to be reliable through application in different regions. To further confirm the credibility of the SW algorithm, the LST will be validated in the future.


Introduction
Land surface temperature (LST) is a key parameter in the physics of land-surface processes regionally and globally; LST has been used in soil moisture estimation [1] and in climatic, hydrological, ecological and biogeochemical studies.Currently, LST can only be obtained over large spatial and temporal scales through remote sensing data, which have attracted much attention in the last three decades [2,3].Consequently, many LST retrieval methods have been proposed from remotely sensed data, particularly multi-channel thermal infrared (TIR) data, and these methods can be roughly grouped into three categories: the single-channel algorithm, multi-channel methods (e.g., the split-window algorithm [4][5][6] and the temperature and emissivity separation method [7]) and multi-time methods (e.g., the temperature-independent spectral indices method [8], two-temperature method [9,10] and the physical day and night algorithm [11]).Some algorithms, such as the iterative spectrally smooth temperature-emissivity separation [12] and the alpha-driven emissivity method [13], have also been proposed to retrieve both LST and emissivity from hyperspectral TIR data.Amongst these methods, the split-window algorithm is the most commonly used, given that this algorithm removes the atmospheric effect and obtains the LST from the linear or nonlinear combination of the brightness temperatures of two adjacent channels centered at 11 and 12 µm.
Landsat 8, also previously called Landsat Date Continuity Mission, extends the remarkable 40-year Landsat record and has enhanced capabilities, including new spectral bands in visible and thermal infrared wavelengths, an improved sensor for signal-to-noise performance and associated improvements in radiometric resolution, and an improved duty cycle to collect significantly more images daily [14].The Thermal Infrared Sensor (TIRS) instrument is one of the major payloads aboard this satellite which can observe the land surface by using the split-window thermal infrared channels (CH10: 10.6 µm to 11.2 µm; CH11: 11.5 µm to 12.5 µm) at a resolution of 100 m.Compared with the TIRS predecessors, namely, the Thematic Mapper (TM) and the Enhanced Thematic Mapper Plus (ETM+), which only have one thermal infrared channel, the TIRS instrument possesses two advantages.First, the TIRS has two thermal infrared channels in the atmospheric window that provide a new LST retrieval opportunity using the widely used split-window algorithm rather than the single-channel method.Second, as shown in Figure 1, the spectral filters of TIRS two bands present narrower bandwidth than that of the thermal band onboard TM and ETM+, the two thermal infrared channels have narrower bandwidths in the TIRS, which can capture finer land surface information, as shown in Figure 1 [15][16][17].Consequently, the LST retrieval algorithm improves our understanding of the new sensor (TIRS) especially on response characteristics to land surface energy, and this algorithm will provide a valuable reference for future studies on the potential capability in terms of LST measurement.
This study aims to develop an operational split-window algorithm (SW) for LST estimation using two thermal infrared channels (TIRS 10: 10.60 µm to 11.19 µm;and TIRS 11: 11.50 µm to 12.51 µm) of the TIRS images.This paper is organized as follows: Section 2 describes the principle associated with LST retrieval through the SW algorithm, and this section also presents the algorithm development for the TIRS data.Section 3 shows the SW algorithm coefficients and the ways to obtain land surface emissivity and atmospheric water vapor.Section 4 illustrates the conducted model sensitivity and consistency using several key parameters, such as brightness temperature noise, land surface emissivity (LSE) and water vapor.Section 5 focuses on the application of the SW algorithm in TIRS images for LST product estimations.Section 6 presents the conclusion.

Split-Window Algorithm Principle
Based on radiative transfer theory, for a cloud-free atmosphere under thermodynamic equilibrium, the channel radiance Bi(Ti) measured on top of the atmosphere (TOA) in a thermal infrared channel of the sensor aboard the satellite is provided with a significant approximation as follows [16,18]: where τi is the effective transmittance of the atmosphere in channel i, and εi is the channel effective surface emissivity.Bi is the Planck function, and Bi (Ts) is the measured radiance if the surface is a black body with a surface temperature Ts (K), _ ↑ and _ ↓ are the upward and downward atmospheric thermal radiance.The first term on the right side of Equation (1) represents the surface emission attenuated by the atmosphere.The second term represents the downward atmospheric thermal radiance reflected by the surface and which reaches the sensor.The third term represents the upward atmospheric emission towards the sensor.From Equation (1), we can deduce that LST retrieval requires knowledge of surface emissivity and atmospheric information.
The split-window algorithm removes the atmospheric effect through differential atmospheric absorption in the two adjacent thermal infrared channels centered at about 11 and 12 µm, and the linear or nonlinear combination of the brightness temperatures is finally applied for LST estimation.Given that this algorithm does not require accurate information about the atmospheric profiles during satellite acquisition, such algorithms have been widely used in LST retrieval from several sensors.A new refinement of the generalized split-window algorithm proposed by Wan (2014) [19] is added with a quadratic term of the difference amongst the brightness temperatures (Ti, Tj) of the adjacent thermal infrared channels, which can be expressed as where Ti and Tj are the TOA brightness temperatures measured in channels i (~11.0 μm) and j (~12.0 µm), respectively; ε is the average emissivity of the two channels (i.e., ε = 0.5 [εi + εj]), whilst Δε is the channel emissivity difference (i.e., Δε = εi − εj); bk (k = 0,1,...7) are the algorithm coefficients derived in the following simulated dataset.

Algorithm Development for Landsat 8
Given the unavailability of a database of in situ LST measurements that coincide with the Landsat 8 overpass, the coefficients bk in Equation (2) are obtained through numerical simulation with different atmospheric and surface conditions.
In our simulation, Thermodynamic Initial Guess Retrieval (TIGR) atmospheric profiles are considered.The TIGR database is constructed by the Laboratoire de Meteorologie Dynamique, and TIGR represents a worldwide set of atmospheric situations (2311 radio soundings) from polar to tropical atmospheres with a column water vapor of 0.1 g/cm 2 to 8 g/cm 2 [20,21].In one of the levels, the profiles with a relative humidity greater than 90% are discarded as under a cloudy condition, which results in a total of 946 atmospheric situations under clear skies with water vapor ranging from 0.06 g/cm 2 to 6.3 g/cm 2 .According to these profiles, the MODTRAN 5.2 atmospheric transmittance/radiance code is used to calculate the channel atmospheric parameters ( , ) in Equation (1) of each atmospheric profile with a spectral integration of the response filters of the two TIRS channels.
For the surface conditions, the LST in the simulation are designed with 7 levels according to the bottom atmospheric temperature T0 of each atmospheric profile, that is, the LST range from T0 − 10 K to T0 + 20 K in a 5 K step.In addition, a total of 53 emissivity spectra in 3 µm to 14 µm are selected from the American Advanced Spaceborne Thermal Emission Reflection (ASTER) emissivity database [22], including 5 water types, 8 man-made target types, 4 vegetation types, 5 rock types, 30 soil types and 1 mineral type.Notably, by contrast to other studies [23,24], man-made targets are considered in this study because the LST products from TIRS are used in urban environmental studies, given the advantages of these products in terms of the finer resolution compared with the Moderate-resolution Imaging Spectroradiometer (MODIS) and Advanced Very High Resolution Radiometer (AVHRR) products, and also the free charge compared with the ASTER product with a resolution of 90 m.
Finally, combined with the atmospheric parameters ( , ), LST and emissivity, the channel radiance at the TOA is determined according to Equation (1), and the brightness temperatures Ti and Tj are obtained from the inverse of Planck's law in the two channels.Given that the field of view (FOV) of the TIRS is about 15 degrees and almost observes the land surface at a nadir direction, meanwhile, the angular effect of atmospheric data, land surface emissivity, and LST as reported in our previous study is also not remarkable [25,26] in FOV = 15 degrees, it is reasonable to ignore the angular variation in the development of SW algorithm.Thus, a total of 350,966 different groups of Ti and Tj, LST, and ε and ∆ε are obtained (946 atmospheres × 53 emissivity × 7 surface temperatures).Thus, the coefficients b0-b7 in Equation (2) can be determined through a statistical regression method.

Algorithm Coefficients
Given that the radiation in the thermal infrared wavelength was attenuated by the atmospheric column water vapor (CWV), we calculated the coefficients b0-b7 in Equation ( 2) independently on the CWV to improve the LST retrieval accuracy.Thus, in our SW algorithm, the CWV was divided into 5 sub-ranges and an overlap of 0.5 g/cm 2 was considered between 2 adjacent sub-ranges, which resulted in [0.0, 2.5], [2.0, 3.5], [3.0, 4.5], [4.0, 5.5] and [5.0, 6.3] g/cm 2 .The CWV was retrieved from a modified split-window covariance and variance ratio method, as stated in Section 3.3.However, given the somewhat unsuccessful CWV retrieval, a group of coefficients for the entire CWV range needed to be calculated to ensure the spatial continuity of the LST product.Table 1 displays the coefficients in different CWV sub-ranges and the root-mean-square error (RMSE) of the temperature that are estimated based on the simulation data.The table shows a significant variation in the coefficients with the CWV sub-ranges, particularly in the heavy CWV loading.The RMSE was smaller than the variation from 0.34 K to 0.93 K, which increased as the CWV increased.To investigate the error details from the algorithm, we further presented the histograms of the temperature difference between the actual Ts in the simulated dataset and the Ts estimated through the SW algorithm using the coefficients (see Table 1) of the 5 CWV sub-ranges.In those figures, the temperature difference evidently fell in the ranges of [−1.0, 1.0] K, which covered about 97.92%, 94.42%, 91.67%, 86.77%, 75.99% and 92.11% of the total cases in Figure 2a-f, respectively.The maximum temperature error was about −3.09K (see Figure 2e), which corresponds to the heavy CWV content.In the case of the entire water vapor range, the RMSE was about 0.87 K.

Determination of LSEs
The classification-based emissivity method (CBEM) [23,[27][28][29] that estimated LSEs from emissivity look-up tables (LUT) according to conventional land cover classification information was applied in this study to obtain LSEs for the SW algorithm.The CBEM is the simplest method in terms of processing, and this method can provide accurate LSEs for LST retrieval as long as the land surfaces are accurately classified and each class has familiar LSEs.The Finer Resolution Observation and Monitoring of Global Land Cover (FROM-GLC), which was the first 30 m resolution global land cover map generated from numerous Landsat TM and ETM+ data [30], is newly released to the public and can be downloaded freely.The FROM-GLC contains 10 types of land covers at level-1, namely, cropland, forest, grassland, shrubland, wetland, waterbody, tundra, impervious, barren land and snow-ice.Table 2 lists the resultant classification scheme with a level-2 hierarchy, which involves the cover-type end-component only.Compared with other popular global land covers, such as the International Geosphere-Biosphere Program (IGBP) [31,32], the FROM-GLC considers the impervious urban class.The high resolution of the FROM-GLC, which combines 100 m resolution TIRS data can significantly improve the accuracy of LSEs.According to the classification scheme in Table 2, three BRDF kernel (bi-directional distribution function) models, i.e., geometrical, volumetric and specular models [24,27], were used to calculate the scene emissivity for each land cover.Considering the variation of the biophysical and structural characteristics in the various land covers, different combinations of these BRDF kernel models were applied for the various land covers.According to the FROM-GLC, we classified 20 land cover types (level-2) into the volumetric BRDF kernel model, and the other 8 types were classified into the specular BRDF kernel model.Given that sparse vegetation was not in the FROM-GLC, we did not consider the geometrical BRDF kernel model.The details are shown in Table 2.Moreover, in this table, the parameter bF = −ln(1 − A), where A is the projected fractional area of the volumetric portion (i.e., vegetation cover).
The volumetric model can be applied in the 50% to 100% coverage range for trees, crops and shrubs, and over the entire range for grasses.σ is the surface roughness and is equal to 0.2.
The vegetation and ground emissivity spectra were selected from the MODIS University of California, Santa Barbara (UCSB) Emissivity Library [33] on the composition of the level-1 and level-2 products in the FROM-GLC, and the spectra were entered in the BRDF kernel models to obtain the scene emissivity.Although the emissivity was reported to vary with the view zenith angle, we found that the angular variation of emissivity was insignificant in the FOV (i.e., 15 degrees) of the TIRS instrument and can be disregarded without leading to obvious errors to the LST.Moreover, we used Fraction of vegetation cover (FVC) to obtain land surface emissivity (LSE) of land cover with temporal variation.It is well known that LSE is a key parameter used in land surface dynamics and varies with the composition of land cover, especially for the natural surface.Assumption that pixel LSE is a linear combination of vegetation and bare soil emissivities of the different land covers, the FVC is estimated from the NDVI, calculated from the red and near-infrared reflectance of Operational Land Imager, another payload on Landsat8, by using the method proposed by Carlson (1997) and Sobrino (2001) [34,35].Finally, an emissivity LUT based on the land cover and FVC was established, and the average band emissivities for the TIRS bands 10 and 11 over different land covers are listed in Table 3.During the retrieval process, the land covers of the pixel from the FROM-GLC and FVC that were calculated from the Normalized Difference Vegetation Index (NDVI) was used to determine the emissivity LUT to obtain the pixel channel emissivities.

Determination of Atmospheric CWV
Inspired by several previous studies [36][37][38][39], a modified split-window covariance and variance ratio (MSWCVR) method was developed to retrieve CWV from the TIRS data.With a vital assumption that the atmosphere is unchanged over the neighboring pixels, the MSWCVR method relates the atmospheric CWV to the ratio of the upward transmittances in two thermal infrared bands, whereas the transmittance ratio can be calculated based on the TOA brightness temperatures of the two bands.Considering N adjacent pixels, the CWV in the MSWCVR method is estimated as In Equation (3a), c0, c1 and c2 are the coefficients obtained from the simulated data; τ is the band effective atmospheric transmittance; N is the number of adjacent pixels (always excluding water and cloud pixels) in a spatial window size n (i.e., N = n × n); Ti,k and Tj,k are the respective brightness temperatures (K) of bands i and j at the TOA level for the kth pixel; and and are the mean or median brightness temperatures of the N pixels for the two bands.Using the aforementioned 946 cloud-free TIGR atmospheric profiles, we first used the new high accurate atmospheric radiative transfer model MODTRAN 5.2 to simulate the band effective atmospheric transmittance, and then we obtained the coefficients through regression, which resulted in c0 = −9.674,c1 = 0.653 and c2 = 9.087.The model analysis indicated that this method will obtain a CWV RMSE of about 0.5 g/cm 2 .The details about the CWV retrieval can be found in [40].

Sensitivity Analysis to Instrument Noises
The noise-equivalent-change-in-temperature (NEΔT) of the TIRS instrument was designed as 0.80 K at 240 K, 0.4 K at 300 K and 0.27 K at 360 K for channel 10, and 0.71 K at 240 K, 0.4 K at 300 K and 0.29 K at 360 K for channel 11 [41].However, this NEΔT was much greater than that of the similar thermal infrared bands on orbit in the MODIS and AHVRR [42,43] and caused significant uncertainty in the land surface temperature retrieval.Fortunately, Ren recently estimated the radiometric noise for the two bands of the TIRS by using images over uniform ground areas, where the actual NEΔT of the two bands was found to be about 0.1 K, which was better than the design specification [44].
To investigate the effect of noise on the LST error, we added a Gaussian noise with a zero mean and a standard deviation equal to 2 × NEΔT (= 0.1, 0.2 and 0.4 K) to the TOA brightness temperatures Ti and Tj in Equation (2).For CWV ∈ [0, 2.5] considering most parts of the 946 atmospheric profiles, we selected the coefficients in this CWV range (see Table 1) as examples to check the variation of the RMSEs affected by the given NEΔT.Thus, the LST error was 0.355 K for NEΔT = 0.1 K, 0.397 K for NEΔT = 0.2 K and 0.531 K for NEΔT = 0.4 K. Compared with the LST error (i.e., 0.34 K) for the no instrument noise, NEΔT = 0.1 K contributed only 4.5% of the error in the retrieved LST, whereas NEΔT = 0.2 K contributed about 16.8% and NEΔT = 0.4 K contributed up to 56.3%.

Sensitivity Analysis to LSEs
According to Equation (2), and ∆ determined the sensitivity from the LSE uncertainties, which can be expressed as 2 5 2 2 To determine the algorithm sensitivities in the different CWV sub-ranges, we consider the coefficients in all the CWV sub-ranges.By using the same regression method mentioned in Section 2.2 and by combining with Equation (4), we can obtain the variation of α and β.Table 4 shows the range, mean and standard derivation of α and β for all the sub-ranges.Based on Table 4, the absolute value of α and β notably increased as the CWV in the atmosphere increased, which indicated better convergence of the sensitivities of LST to and ∆ in a wet atmospheric condition compared with those in a dry atmospheric condition.According to Equation ( 2), the combined uncertainty from and ∆ contributed to the LST error δLST, which can be written as Table 5 shows that the LST error was caused by 1% of the LSE uncertainties for the different CWV ranges.From this table, the LST error was determined as [0.73, 1.52] K with a mean of 1.02 K and a standard deviation of 0.10 K for the CWV sub-range of [0.0, 2.5] g/cm 2 .The error generally decreased as the CWV increased, and the minimum was obtained for the highest atmospheric condition consequently.In addition, the range of the LST error in the wet atmospheric condition was narrower than that in the dry atmospheric condition.This finding indicates that in wet atmospheric conditions, the dominant factor that affects the LST error is the uncertainty from the atmospheric information, rather than from that of the surface conditions.In this case, the LST retrieval was insensitive to the error included in the LSEs.

Sensitivity Analysis to the Atmospheric CWV
As shown in Figure 2, the CWV is important in improving the accuracy of LST retrieval.However, we did not directly use the CWV in Equation ( 2) to estimate the LST, but we applied this parameter to determine the algorithm coefficients (see Table 1).As mentioned in Section 3.3, the CWV retrieval had an error of about 0.5 g/cm 2 theoretically.Thus, misclassifying CWV from the correct sub-range to the incorrect sub-range was possible, which can result in wrong coefficients for the SW algorithm.Table 6 lists the LST errors caused by using the wrong coefficients in the adjacent CWV sub-ranges.Given that the CWV error varies from −0.5 g/cm 2 to 0.5 g/cm 2 , we considered only the cases of misclassifying CWV in the adjacent CWV sub-ranges.Thus, for a true CWV value of 2.6 g/cm 2 , the retrieved water vapor value may be ranged in [0.0, 2.5], [2.0, 3.5] or [3.0, 4.5] g/cm 2 .In Table 6, the diagonal direction shows the LST errors with the correct coefficients in the corresponding CWV sub-ranges, which were similar to the results in Figure 2. The other values are the results from the incorrect coefficients, which are denoted as δLSTinc for convenience in the following discussion.The δLSTinc was less than 1.5 K with a CWV of less than 4.0 g/cm 2 , but δLSTinc increased dramatically with a CWV of more than 4.0 g/cm 2 .However, for a CWV in the sub-ranges of [2.5, 3.0], [3.5, 4.0], [4.5, 5.0] or [5.5, 6.0] g/cm 2 , the CWV retrieval value using the MSWCVR method may fall in the CWV adjacent sub-ranges, which consequently decreases the LST accuracy.For instance, a CWV is 4.8 g/cm 2 and belongs to a sub-range of [4.5, 5.0] g/cm 2 .If this CWV is misclassified into [3.5, 4.0], the δLSTinc will be 1.29 K, whereas if this CWV is misclassified into [5.5, 6.0], the δLSTinc will reach 2.45 K. To reduce the influence of the CWV error on the LST, for a CWV within the overlap of two adjacent CWV sub-ranges, we first use the coefficients from the two adjacent CWV sub-ranges to calculate the two initial temperatures and then use the average of the initial temperatures as the pixel LST.For example, the LST pixel with a CWV of 2.1 g/cm 2 is estimated by using the coefficients of [0.0, 2.5] and [2.0, 3.5].This process initially reduces the δLSTinc and improves the spatial continuity of the LST product.

Comparison amongst Different Split-Window Algorithms
To date, two other SW algorithms were proposed for LST retrieval from TIRS data.Jiménez-Muñoz et al. (2014) [45] presented an SW algorithm based on the structure suggested by Sobrino et al. (1996) [46], and Rozenstein et al. (2014) [47] utilized a first-order Taylor-series linearization of the radiative transfer equation and addressed a general form based on the work of Qin et al. (2001).The formats of the current SW algorithms are expressed as follows: Jiménez-Muñoz et al.: In the above equations, dk (k = 0, 1…6) and ek (k = 1, 2, 3, 4) are the algorithm coefficients; w is the CWV; ε and ∆a are the average emissivity and emissivity difference of two adjacent thermal channels, respectively, which are similar to Equation (2); and fk (k = 0 and 1) is related to the influence of the atmospheric transmittance and emissivity, i.e., fk = f(εi,εj,τi,τj).Note that the algorithm (Equation (6a)) proposed by Jiménez-Muñoz et al. added CWV directly to estimate LST.Rozenstein et al. used CWV to estimate the atmospheric transmittance (τi, τj) and optimize retrieval accuracy explicitly.Therefore, if the atmospheric CWV is unknown or cannot be obtained successfully, neither of the two algorithms in Equations (6a) and (6b) will work.By contrast, although our algorithm also needs CWV to determine the coefficients, this algorithm still works for unknown CWVs because the coefficients are obtained regardless of the CWV, as shown in Table 1.We first obtained the coefficients of Equations (6a) and (6b) using the same simulated dataset in our algorithm development, and then we analyzed the difference between the SW algorithms.Table 7 presents the LST errors caused by the different algorithms for various CWV sub-ranges.From this table, we know that the errors of all the algorithms were close to one another for a CWV of less than 3.5 g/cm 2 , which is much less than 1.0 K.However, under wet atmospheric conditions, the LST error increased quickly as the CWV increased, especially for the algorithm of Rozenstein et al., Jiménez-Muñoz's algorithm had similar results with our algorithm.Moreover, we added an uncertainty of ±0.5 g/cm 2 to the CWV, and we found that the LST error was about 0.8 K for the CWV of −of t tha 2 and 1.1 K for the CWV of +0.5 g/cm 2 for the algorithm of Jiménez-Muñoz et al.; the results were somewhat better than those of our algorithm, as shown in Table 7, probably because of the direct usage of CWV in the algorithm of Jiménez-Muñoz et al. (Equaiton (6a)) to reduce the influence of the CWV on the LST retrieval accuracy.As stated above, all three algorithms relied on the CWV input and were impractical without this parameter.To deal with this potential problem, our algorithm also obtained the coefficients in Equation (2) for the entire CWV range, and this group of coefficients only resulted in an LST error of about 0.87 K for all the CWV, about 0.46 K for the CWV sub-range of [0.0, 2.5] g/cm 2 and about 1.11 K for a CWV of less than 3.5 g/cm 2 .This CWV range contains most cases of atmospheric moisture in polar, mid-latitude and tropical profiles.Therefore, compared with the other two algorithms, our algorithm is more practical and can even obtain LST with high accuracy for cases even without known CWV information.

Application of the Split-Window Algorithm
Based on the above discussions, the input parameters to drive the new SW algorithm includes the brightness temperatures (Ti and Tj) of the two adjacent bands of the TIRS, FROM-GLC land cover products and emissivity lookup table, which are a fraction of the FVC that can be estimated from the red and near-infrared reflectance of the Operational Land Imager (OLI).These parameters can be easily obtained, which makes our algorithm generally operational in practice.Figure 3 shows the main flowchart of the LST retrieval LST from the TIRS data, where the clouds in the images are eliminated by using the band quality files along with the OLI and TIRS data, which are consequently removed before the LST retrieval.The output can contain LST products and emissivity images in the two channels and in the CWV product.Based on Figure 4b, the temperatures of the urban impervious surface are higher than those around the vegetated surface and water bodies, which consequently lead to the urban heat island effect.In Figure 4e, the temperatures in the deserts are much higher than those in the oasis areas and water bodies, and the temperature difference in this area can reach up to 20 or 30 K. In this region, the LST is highly correlated to the surface types and soil moistures, and then becomes a significant indicator of the environmental conditions.Compared with the histograms shown in Figure 4c and f, the spatial temperature difference of Figure 4d is more significant than that of Figure 4a, which is probably because the atmospheric moisture in northwest China is generally low from the lack of surface evapotranspiration and input water vapor from other places; the LST is mainly determined by the surface conditions, whereas the regions of Figure 4a has a relatively higher water vapor, and the spatial temperature difference is smoothened by the atmospheric effect.Note that, the validation of the LST estimated from the Landsat 8 TIRS data was not conducted in this paper because it is discussed in another separate paper by using in situ measured temperature synchronously with Landsat 8 overpass in several uniform surface and also with some other products, such as MODIS and ASTER LST products.

Discussions
As stated above, this paper proposed a practical algorithm to obtain LST from Landsat 8, but for its good performance several attentions should be paid.First, as pointed by some studies [44,49], data gap existed in some of current TIRS images because of the uncertainly of pixel-to-pixel radiometric calibration.As a result, the used water vapor will show spatial discontinuity, and consequently, the LST products will show some obvious blocks in their images.Therefore, it is very necessary and important to investigate the pixel-to-pixel radiometric calibration the TIRS image.Besides, the spatial scales of land cover product (30 m), water vapor (1 km), TIRS pixel (100 m), and FVC data (30 m) can reduce the LST retrieval accuracy if spatial mis-registration exists between different input auxiliary data.
As mentioned in Section 2.2, the used 946 TIGR atmospheric profiles is composed of 812 profiles for CWV ∈ [0.0, 2.5] g/cm 2 , 100 profiles for CWV ∈ [2.0, 3.5] g/cm 2 , 76 profiles for CWV ∈ [3.0, 4.5] g/cm 2 , 36 profiles for CWV ∈ [4.0, 5.5] g/cm 2 and 15 profiles for CWV ∈ [5.0, 6.3] g/cm 2 .The dry atmosphere takes great part in the used atmospheric profiles, but our algorithm developed several groups of coefficients depending on the CWV sub-ranges, so the current algorithm will not cause significant uncertainty to the final LST retrieval with known CWV, but it will lead some error to the LST result for the cases without input CWV.To reduce this effect, we are trying to find more representative profiles to optimize the current algorithm in the coming future work.
In addition, newly released from U.S. Geological Survey (USGS), discrepancies have been noted between calibrated Landsat 8 TIRS Bands 10 and 11 data and these calibration uncertainty results about more than 2 K errors compared with surface-water temperature field campaign.Although, the radiometric and geometric refinements have been implemented in thermal infrared band offsets soon and improved the data accuracy to −2.1 ± 0.8 K for band 10 and −4.4 ± 1.75 K for band 11 at a 300 K brightness temperature respectively, this accuracy still cannot meet the requirements of algorithm application [48,49].Meanwhile, according to the reprocessing campaign results, this uncertainty is scene-dependent and probably related to out-of-field response in the TIRS instrument, and gives us an obstacle in algorithm validation, for we cannot analyze whether errors are caused by TIRS sensor or the algorithm.

Conclusions
Land surface temperature, as one of the key variables, describes surface states and processes critical in studies of climate, hydrology, ecology and human health.This study proposed a practical split-window algorithm for LST retrieval from the Landsat 8 TIRS data using two thermal infrared channels: CH10 (10.6 µm to 11.2 µm) and CH11 (11.5 µm to 12.5 µm).The new SW algorithm estimated the LST by using a nonlinear combination of brightness temperatures, and the algorithm coefficients were designed to be dependent on water vapor, which was obtained from the TIRS data itself using the MSWCVR method.According to the principle of the CBEM, channel emissivities were obtained from the FROM-GLC products at 30 m, and FVC was calculated from the red and near-infrared images aboard the Landsat 8.
The simulation results showed that the new algorithm can obtain LST with an accuracy of better than 1.0 K.A model analysis with respect to the noise that included in the measured brightness temperature, emissivity and water vapor was conducted, and the results showed that the new algorithm was robust in LST retrieval.
Our algorithm was compared with two other existing algorithms, and the results showed that the accuracy of our algorithm was close to that of Jiménez-Muñoz and was better than that of Rozenstein.However, because we used water vapor to determine the coefficients of the SW algorithm rather than directly applying this parameter in the algorithm equations as other algorithms do, and also because we obtained the algorithm coefficients for all the water vapor ranges, the new SW algorithm in this paper was more practical than the two existing algorithms, especially in conditions where the water vapor is unknown.The new SW algorithm was applied to two different regions, and the estimated LST presented higher value in urban and desert areas than the vegetated surface and waterbody.The results further proved the credibility and reliability of our proposed SW algorithm.
Considering the strip problem, ghost signal caused by stray light and a time-varying absolute calibration error for TIRS, the validation exercise is still a tough problem and will be discussed in another separate paper by using ground measured surface temperature at different uniform areas and with other surface temperature products, such as MODIS [6] and ASTER [7] products.

Figure 1 .
Figure 1.Spectral response functions of the thermal bands of the different sensors on board the Landsat platforms.

Figure 4 .
Figure 4.The LST (Land Surface Temperature) retrieved from TIRS (Thermal Infrared Sensor) data at urban area of Beijing city and Gobi in northwest China.(a) and (d) are the false color images; (b) and (e) are the retrieved LST; while (c) and (f) present LST histograms.

Table 1 .
The coefficients bk (k = 0,1,…7) in different atmospheric column water vapor (CWV) sub-ranges and the root-mean-square error (RMSE) of the temperature estimated based on the simulation data.

Table 2 .
Land cover schemes in FROM-GLC (Finer Resolution Observation and Monitoring of Global Land Cover) at levels 1 and 2 and their components.

Table 3 .
Average emissivity for two TIRS (Thermal Infrared Sensor) channels at different land covers of FROM-GLC (Finer Resolution Observation and Monitoring of Global Land Cover)

Table 5 .
LST (Land Surface Temperature) errors caused by 1% uncertainties in LSEs (Land Surface Emissivities) for different ranges of CWV (Column Water Vapor).

Table 6 .
LST (Land Surface Temperature) errors caused by using wrong coefficients in adjacent CWV (Column Water Vapor) sub-ranges.

Table 7 .
LST (Land Surface Temperature) error for different split-window (SW) algorithms at different sub-ranges of CWV (Column Water Vapor).The second, third and fourth columns correspond to the RMSEs (Root-Mean-Square Errors) of Jiménez-Muñoz et al., Rozenstein et al., and our algorithm, while the last column is RMSE from our algorithm using the coefficients derived from all water vapor range [0, 6.3] g/cm 2 when no CWV information can be obtained.