Land surface temperature retrieval from Landsat 8 data and validation with geosensor network

Abstract

A method for the retrieval of land surface temperature (LST) from the two thermal bands of Landsat 8 data is proposed in this paper. The emissivities of vegetation, bare land, buildings, and water are estimated using different features of the wavelength ranges and spectral response functions. Based on the Planck function of the Thermal Infrared Sensor (TIRS) band 10 and band 11, the radiative transfer equation is rebuilt and the LST is obtained using the modified emissivity parameters. A sensitivity analysis for the LST retrieval is also conducted. The LST was retrieved from Landsat 8 data for the city of Zoucheng, Shandong Province, China, using the proposed algorithm, and the LST reference data were obtained at the same time from a geosensor network (GSN). A comparative analysis was conducted between the retrieved LST and the reference data from the GSN. The results showed that water had a higher LST error than the other land-cover types, of less than 1.2°C, and the LST errors for buildings and vegetation were less than 0.75°C. The difference between the retrieved LST and reference data was about 1°C on a clear day. These results confirm that the proposed algorithm is effective for the retrieval of LST from the Landsat 8 thermal bands, and a GSN is an effective way to validate and improve the performance of LST retrieval.

Keywords

Land surface temperature (LST) / split-window algorithm / emissivity / Landsat 8

Introduction

Land surface temperature (LST) plays an important role in the research into urban heat island effects, environmental monitoring, and agricultural analyses, and significant progress has been made in these applications with a series of remote sensing satellites being launched ( Lu and Weng, 2006; Weng and Lu, 2008; Patel et al., 2012; Krishna and Sharma, 2013). Landsat 8 is the follow-on satellite to Landsat 7. Landsat 7 and Landsat 5 have only one thermal infrared band, but the Landsat 8 Thermal Infrared Sensor (TIRS) contains two thermal infrared bands, with wavelength ranges of 10.60–11.19 μm for band 10 and 11.50–12.51 μm for band 11. Using these two thermal infrared bands, two types of methods have been developed to retrieve LST from Landsat data: the single infrared channel methods and the split-window methods.

During the past decades, a number of single infrared channel methods and split-window methods have been developed for estimating LST from space. The single-channel methods require a radiative transfer model and atmospheric profiles, which can be provided either by satellite sounding or a conventional radiosonde ( Price, 1983; Susskind et al., 1984; Wan and Dozier, 1996). The common single-channel methods were developed for Landsat thermal infrared images ( Qin et al., 2001b; Jiménez-Muñoz, 2003; Jiménez-Muñoz et al., 2009). A number of algorithms, including single-channel algorithms and mono-window algorithms, have been modified with other new types of images, such as HJ-1B and CBERS thermal infrared images ( Yong et al., 2006; Sibo et al., 2008; Zhou et al., 2011). The split-window methods correct the atmospheric effects based on differential absorption in the adjacent infrared bands ( Price, 1984; Harris and Mason 1992; Sobrino et al., 1993; Qin et al., 2001a; Mao et al., 2005a; Du et al., 2014; Rozenstein et al., 2014; Wan, 2014). However, these split-window algorithms require many input parameters, such as land surface emissivities, and they also require more than two channels. Moreover, the spectral response functions and the wavelengths of different sensors vary, and some crucial parameters such as land surface emissivities and atmospheric transmittance need to be recalculated. Thus, the old functions are not suitable for use with a new sensor.

Although the U.S. Geological Survey (USGS) stated that there was some uncertainty about band 11 and suggested that only band 10 should be used for LST retrieval (http://landsat.usgs.gov/calibration_notices.php), we believe that it is appropriate to attempt to design a split-window algorithm for Landsat 8 data and to estimate the LST errors under this condition. Many researchers have developed modified algorithms to retrieve LST from Landsat 8 data ( Du et al., 2014; Jimenez-Munoz et al., 2014; Rozenstein et al., 2014; Yang et al., 2014). These methods are mainly based on algorithms proposed for MODIS or AVHRR imagery ( Qin et al., 2001a; Jiménez-Muñoz, 2003; Wan, 2014), and have been widely used in different fields, with mean errors of less than 1.5 K for Jimenez-Munoz’s method ( Jimenez-Munoz et al., 2014) and 0.93 K for Rozenstein’s method (Rozenstein et al., 2014) when tested with simulated data only. However, these methods just modify the parameters of the original algorithms for the Landsat 8 data, and they have not been validated with ground-truth data because the true LST values were not collected and compared.

In this paper, we propose a novel linear equation between band 10 and band 11 based on Planck’s radiation function, and we estimate the LST from the radiative transfer equation. Since the thermal infrared bands of Landsat 8 have different spectral response functions to Landsat 7, the parameters used to retrieve the LST have been redefined, including the emissivities of the typical land surface estimation equations. A sensitivity analysis was performed, and the algorithm was validated using ground data measured by a geosensor network (GSN) over the city of Zoucheng, Shandong Province, China. Finally, the errors in the LST retrieval were analyzed and are discussed based on the synchronous GSN data.

LST algorithms

The modified algorithm

We use Qin’s algorithm as the basis of the improved algorithm ( Qin et al., 2001b). The improved algorithm is based on correction of the radiative transfer equation, which states that the sensor-observed radiance is always impacted by the atmospheric transmittance and ground emission. The sensor-observed radiance for the thermal infrared bands can be expressed as:

(1) B_i(T_i)={\tau }_i{\epsilon }_i{B}_i(T_s)+[\mathrm{1}+{\tau }_i(\mathrm{1}-{\epsilon }_i)](\mathrm{1}-{\tau }_i)B_i(T_a).

The above Eq. (1) for the observed radiance was deduced by Qin ( Qin et al., 2001b), where i represents the order of the thermal infrared bands of Landsat 8, which can be 10 or 11. T_s , T_i , and T_a represent the LST, the brightness temperature of the thermal infrared bands, and the effective mean atmospheric temperature in Kelvin, respectively. Accordingly, B_i(T_s) , B_i(T_i) , and B_i(T_a) are the radiance of T_s , T_i , and T_a . {\tau }_i is the atmospheric transmittance, and {\epsilon }_i is the land surface emissivity. As is shown in Eq. (1), it is impossible to obtain T_s without these three parameters, T_a , {\tau }_i , and {\epsilon }_i , and thus Eq. (1) is an abnormal state equation. In this situation, we usually use three additional equations to calculate T_s . In general, {\tau }_i can be simulated by an atmospheric correction model (such as MODTRAN or LOWTRAN) with the total water vapor content ( Sobrino et al., 1993). Meanwhile, {\epsilon }_i can be estimated by the NDVI parameter ( Van de Griend and Owe, 1993). An alternative equation, which describes the relationship between radiance and temperature, can be deduced by simplifying Planck’s radiation function ( Qin et al., 2001b; Jiménez-Muñoz, 2003; Mao et al., 2005b).

In this section, we describe a novel equation for the relationship between the two infrared bands, based on Planck’s radiation function. According to Planck’s law, the thermal radiation of a black body at a given wavelength and temperature can be expressed as:

(2) B_{\lambda }(T)={\displaystyle \frac{C_{\mathrm{1}}}{{\lambda }^{\mathrm{5}}(e^{C_{\mathrm{2}}/\lambda T}-\mathrm{1})}},

where B_{\lambda }(T) is the spectral radiance of the black body, which is generally measured in Wm^‒2·sr^‒1·μm^‒1; \lambda is the wavelength in micrometers; C_{\mathrm{1}} and C_{\mathrm{2}} are the Planck’s constants, with C_{\mathrm{1}}=\mathrm{1.19104356}\times {\mathrm{10}}^{-\mathrm{16}} Wm²·sr^‒1 and C_{\mathrm{2}} =14387.685 μm·K; and T is the temperature in Kelvin.

For the different wavelength ranges of the Landsat 8 TIRS bands, the effective wavelength and the combined Planck radiance can be calculated as:

\{\begin{array}{l}\begin{array}{ll}{\lambda }_{\mathrm{ef}\mathrm{f}}={{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}\lambda f_i(\lambda )d\lambda /{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{f}_i(\lambda )d\lambda & (\mathrm{3}a)\end{array}\\ \begin{array}{ll}{B}_i(\lambda )={{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{B}_i(\lambda )f_i(\lambda )d\lambda /{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{f}_i(\lambda )d\lambda & (\mathrm{3}b)\end{array}\end{array}

where f_i(\lambda ) is the spectral response function. From Eq. (3a), we obtain the values of 10.901μm for {\lambda }_{\mathrm{10}} and 12.011 μm for {\lambda }_{\mathrm{11}} . B_i(\lambda ) is the Planck radiance for band 10 or band 11. The Landsat 8 TIRS has two thermal bands, and the spectral response functions of Landsat 8 are as shown in Fig. 1(a).

According to Planck’s radiation function and Eq. (3b), the relationships between temperature and radiance for band 10 and band 11 are shown in Fig. 1(b). The difference between band 10 and band 11 increases with T, and the dotted line indicates that they are satisfied with an approximately linear relationship. Thus, we presume that the relationship between the radiance of Landsat 8 band 10 and band 11 can be described as:

(4) B_{\mathrm{11}}(T)=a{B}_{\mathrm{10}}(T)+b.

The parameters and the coefficient of determination R^{\mathrm{2}} for the regression equation (Eq. (4)) are shown in different temperature ranges in Table 1. Values of R^{\mathrm{2}} greater than 0.999 indicate that the regression equation has a more effective correlation than the other algorithms for retrieving LST (Rozenstein et al., 2014; Yang et al., 2014). Thus, we consider that Eq. (4) is reliable for the relationship between B_{\mathrm{10}}(T) and B_{\mathrm{11}}(T) . Using the different features of the wavelength ranges and spectral response functions, we have different effective wavelengths and parameters a and b , as shown in Table 1.

In order to obtain T_s from Eq. (1), we linearize B_{\mathrm{11}}(T_a) and B_{\mathrm{10}}(T_s) by Eq. (3). Eq. (1) for Landsat 8 band 10 and band 11 can then be expressed as:

(5) B_{\mathrm{10}}(T_{\mathrm{10}})=C_{\mathrm{10}}{B}_{\mathrm{10}}(T_s)+D_{\mathrm{10}}{B}_{\mathrm{10}}(T_a),

(6) B_{\mathrm{11}}(T_{\mathrm{11}})=C_{\mathrm{11}}[a{B}_{\mathrm{10}}(T_s)+b]+D_{\mathrm{11}}[a{B}_{\mathrm{10}}(T_a)+b].

For simplification, we define:

(7) \{\begin{array}{l}{C}_{\mathrm{10}}={\tau }_{\mathrm{10}}{\epsilon }_{\mathrm{10}}\\ D_{\mathrm{10}}=[{\tau }_{\mathrm{10}}(\mathrm{1}-{\epsilon }_{\mathrm{10}})+\mathrm{1}](\mathrm{1}-{\tau }_{\mathrm{10}})\\ C_{\mathrm{11}}={\tau }_{\mathrm{11}}{\epsilon }_{\mathrm{11}}\\ D_{\mathrm{11}}=[{\tau }_{\mathrm{11}}(\mathrm{1}-{\epsilon }_{\mathrm{11}})+\mathrm{1}](\mathrm{1}-{\tau }_{\mathrm{11}})\end{array}.

The effective mean atmospheric radiance B_{\mathrm{10}}(T_a) and Eq. (6) can be transformed as:

(8) B_{\mathrm{10}}(T_a)={\displaystyle \frac{B_{\mathrm{10}}(T)-C_{\mathrm{10}}{B}_{\mathrm{10}}(T_s)}{D_{\mathrm{10}}}},

(9) B_{\mathrm{11}}(T_{\mathrm{11}})=C_{\mathrm{11}}[a{B}_{\mathrm{10}}(T_s)+b]+D_{\mathrm{11}}[a{\displaystyle \frac{B_{\mathrm{10}}(T_{\mathrm{10}})-C_{\mathrm{10}}{B}_{\mathrm{10}}(T_s)}{D_{\mathrm{10}}}}+b].

By Eq. (9), we have:

(10) B_{\mathrm{10}}(T_s)=E{B}_{\mathrm{10}}(T_{\mathrm{10}})+F{B}_{\mathrm{11}}(T_{\mathrm{11}})+G.

For simplification, we define:

(11) \{\begin{array}{l}E={\displaystyle \frac{-{\displaystyle \frac{D_{\mathrm{11}}}{D_{\mathrm{10}}}}}{C_{\mathrm{11}}-C_{\mathrm{10}}{\displaystyle \frac{D_{\mathrm{11}}}{D_{\mathrm{10}}}}}}\\ F={\displaystyle \frac{\mathrm{1}}{a{C}_{\mathrm{11}}-a{C}_{\mathrm{10}}{\displaystyle \frac{D_{\mathrm{11}}}{D_{\mathrm{10}}}}}}\\ G={\displaystyle \frac{-(C_{\mathrm{11}}+D_{\mathrm{11}})b}{a{C}_{\mathrm{11}}-a{C}_{\mathrm{10}}{\displaystyle \frac{D_{\mathrm{11}}}{D_{\mathrm{10}}}}}}\end{array}.

After obtaining the ground radiance of band 10, we get a backward equation of Planck’s radiation equation ( Wukelic et al., 1989) for T_s . Therefore, with this relationship, we obtain:

(12) \{\begin{array}{l}{T}_s={\displaystyle \frac{K_{\mathrm{2}}}{\ln (\mathrm{1}+K_{\mathrm{1}}/B_{\mathrm{10}}(T_s))}}\\ K_{\mathrm{1}}=C_{\mathrm{1}}/{\lambda }_{\mathrm{10}}^{\mathrm{5}}\\ K_{\mathrm{2}}=C_{\mathrm{2}}/{\lambda }_{\mathrm{10}}\end{array},

where {\lambda }_{\mathrm{10}} is the effective wavelength of band 10. Eqs. (10)–(12) are the proposed algorithms for Landsat 8 TIRS images. With this model, we require two parameters: one is the atmospheric transmittance estimated by MODTRAN 4, and the other is the land surface emissivity.

Land surface emissivity retrieval

Land surface emissivity is a proportionality factor that scales black-body radiance to predict emitted radiance, and it has a significant impact on the retrieval of LST in that an emissivity error of±0.025 can yield an LST error of±2 K ( Schädlich et al., 2001). Different approaches have been applied to predict land surface emissivity from NDVI values ( Van de Griend and Owe, 1993; Becker and Li 1995; Valor and Caselles 1996; Sobrino et al., 2004; Qin et al., 2006). Van de Griend and Owe (1993) found a significant relationship between NDVI and thermal emissivity, and the correlation coefficient was 0.94 in the logarithmic transformation. Valor and Caselles (1996) focused on the 10.5–12.5 μm region of the wavelengths, and they went on to explain the experimental behavior observed by Van de Griend and Owe. They then proposed a theoretical model that relates the emissivity to the NDVI ( Valor and Caselles, 1996). Becker and Li (1995) presented an improved temperature-independent spectral index for deriving emissivities from AVHRR data. Sobrino et al.(2004) made an assumption that the land surface is soil when the NDVI value is below 0.2, and vegetation when it is over 0.5. Qin et al. (2006) viewed urban areas as being composed of building surfaces and vegetation, and they estimated the emissivities by using the land surface thermal radiation equation and the known emissivities of different land-cover types. Here, we adopt the NDVI threshold method which was proposed by Qin based on Landsat TM6, and we improve the method for use with Landsat 8 data.

The proposed method contains two models: one uses a certain NDVI threshold to distinguish between soil pixels and vegetation pixels; the other model distinguishes between building pixels and vegetation pixels. For these mixed pixels, the method uses the following simplified equations:

(13) \epsilon ={\epsilon }_v{P}_v{R}_v+{\epsilon }_s(\mathrm{1}-P_v)R_s+{\delta }_{\epsilon },

(14) \epsilon ={\epsilon }_v{P}_v{R}_v+{\epsilon }_m(\mathrm{1}-P_v)R_m+{\delta }_{\epsilon },

where {\epsilon }_v , {\epsilon }_m , and {\epsilon }_s are the emissivities of soil, buildings, and vegetation, P_v is the proportion of vegetation, and {\delta }_{\epsilon } is the error correction term of the model. R_v , R_s , and R_m are the temperature ratios of vegetation, soil, and buildings, which can be expressed as:

(15) \{\begin{array}{l}{R}_v=\mathrm{0.9332}+\mathrm{0.0585}{P}_v\\ R_s=\mathrm{0.9902}+\mathrm{0.1068}{P}_v\\ R_m=\mathrm{0.9886}+\mathrm{0.1287}{P}_v\end{array},

where {\delta }_{\epsilon } can be calculated as:

(16) {\delta }_{\epsilon }=\{\begin{array}{ll}0{.0038\mathit{P}}_v;& P_v\le \mathrm{0.5}\\ \mathrm{0.0038}(\mathrm{1}-P_v);& P_v>\mathrm{0.5}\end{array}.

The proportion of vegetation P_v can be obtained from the NDVI, according to:

(17) P_v={\left[{\displaystyle \frac{\mathrm{NDVI}-{\mathrm{NDVI}}_s}{{\mathrm{NDVI}}_v-{\mathrm{NDVI}}_s}}\right]}^{\mathrm{2}},

where the {\mathrm{NDVI}}_s and {\mathrm{NDVI}}_v values are the thresholds of soil pixels and pixels of full vegetation. {\mathrm{NDVI}}_s=\mathrm{0.2} and {\mathrm{NDVI}}_v=\mathrm{0.7} ( Sobrino et al., 2004).

The typical land surface emissivities are shown in Table 2. Landsat 8 has two thermal infrared bands, and the wavelength ranges and the spectral response functions of the sensors are quite different to the Landsat 5 sensor. Thus, we need to modify some parameters of Qin’s retrieval algorithm, based on the characteristics of Landsat 8 TIRS.

In Qin’s land surface emissivity model, the water is extracted, and the emissivity of the water is defined as a fixed value. Thus, we define the emissivities of vegetation, soil, buildings, and water. Since the different thermal infrared bands have different emissivities for the same land-surface types, the typical land surface emissivities for the Landsat 8 TIRS bands need to be recalculated.

The spectral data used here were obtained from the Johns Hopkins University (JHU) spectral library, and were provided by Jack Salisbury at JHU (http://asterweb.jpl.nasa.gov/speclib/). The spectral functions are the reflectances of the typical land features. We can transform these reflectance data into an emissivity file, based on Kirchhoff’s law ( Sasamori, 1999). The emissivities of distilled water, three different kinds of vegetation, and four kinds of buildings and soils are used here to separately represent the typical land surface emissivities of water, vegetation, buildings, and soil (Fig. 2).

As for the typical land surface emissivity of the thermal infrared bands, the formula ( Wan and Dozier, 1996) is:

(18) {\epsilon }_i={\displaystyle \frac{{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}[p_{\mathrm{1}}{\epsilon }_{\mathrm{1}i}(\lambda )+p_{\mathrm{2}}{\epsilon }_{\mathrm{2}i}(\lambda )]f_i(\lambda )d\lambda }{{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{f}_i(\lambda )d\lambda }}.

This formula is used to calculate the emissivities of the mixed pixels. {\epsilon }_{i\mathrm{1}}(\lambda ) and {\epsilon }_{i\mathrm{2}}(\lambda ) are the emissivities of the different land-surface types in the pixels. p_{\mathrm{1}} and p_{\mathrm{2}} are the proportions. If it is a pure pixel, the emissivity can be calculated as:

(19) {\epsilon }_i={\displaystyle \frac{{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{\epsilon }_i(\lambda )f_i(\lambda )d\lambda }{{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{f}_i(\lambda )d\lambda }}.

For Qin’s emissivity model, the typical land surface emissivities are estimated by averaged values. Since the typical land-surface spectral functions are quite different from those of Landsat 5, we need to recalculate the emissivities of vegetation, bare land, buildings, and water, based on the original emissivity of Landsat TM6. Eq. (19) is therefore expressed as:

(20) {\displaystyle \frac{{\epsilon }_i}{{\epsilon }_j}}={\displaystyle \frac{{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{\epsilon }_i(\lambda )f_i(\lambda )d\lambda }{{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{f}_i(\lambda )d\lambda }}/{\displaystyle \frac{{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{\epsilon }_j(\lambda )f_j(\lambda )d\lambda }{{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{f}_j(\lambda )d\lambda }},

(21) {\epsilon }_i=[{\displaystyle \frac{{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{\epsilon }_i(\lambda )f_i(\lambda )d\lambda }{{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{f}_i(\lambda )d\lambda }}/{\displaystyle \frac{{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{\epsilon }_j(\lambda )f_j(\lambda )d\lambda }{{{\displaystyle \int }}_{{\lambda }_{\min }}^{{\lambda }_{\max }}{f}_j(\lambda )d\lambda }}]\times {\epsilon }_j.

Subscript j represents the parameters of Landsat 5 band 6 in Eq. (20) and Eq. (21), where {\epsilon }_i(\lambda ) is the emissivity of water, vegetation, soils, and buildings, as shown in Fig. 2. {\epsilon }_j is the emissivity from Qin’s Landsat 5 band 6 model. According to Eq. (21), it indicates that the results of the emissivities are correlated with the variation trends of land surface emissivity in Fig. 2 and the values of the typical land surface emissivities from TM6. Table 2 shows the different emissivities of the typical land-surface types for Landsat 8 TIRS and Landsat 5 TM6. This indicates that the values of these typical land surface emissivities for Landsat 8 TIRS are quite close to the values of Landsat 5, and the emissivities are only slightly different between band 10 and band 11. Thus, errors in the LST will be generated if different values of land surface emissivities are not used in the retrieval models.

Estimation of atmospheric transmittance

When the land surface emissivity is known, the LST still cannot be calculated without the atmospheric transmittance {\tau }_i . The atmospheric transmittance estimated by LOWTRAN ( Kneizys et al., 1983; Kneizys et al., 1988) and MODTRAN ( Berk et al., 1987, 1999) have been widely used in the development of LST algorithms, and water vapor content has also been extensively used as the determinant in the estimation of atmospheric transmittance. A number of papers have listed regression equations for Landsat 8 TIRS bands (Rozenstein et al., 2014; Yang et al., 2014), and these approaches obtain the atmospheric transmittance from MODTRAN 4 software that simulates the atmospheric conditions.

Parameter adjustment for Mao’s algorithm

In order to make a comparison with the proposed algorithm, we modified Mao’s algorithm ( Mao et al., 2005b) so that it was suitable for retrieving LST from Landsat 8 data. This algorithm was developed based on Qin’s algorithm ( Qin et al., 2001b), and it simplifies the Planck function in the linear equation:

(22) B_i(T_i)=a_i+b_i{T}_i.

Parameter i in Eq. (22) represents the band number of Landsat TIRS. According to the effect of the wavelength of the Landsat 8 TIRS bands, parameters a_i and b_i in different temperature ranges were computed, and the results are shown in Table 3.

Using Eq. (23), the simplified radiative transfer equation in Eq. (1) can be expressed as:

(23) a_i+b_i{T}_i={\tau }_i{\epsilon }_i(a_i+b_i{T}_s)+[{\tau }_i(\mathrm{1}-{\epsilon }_i)+\mathrm{1}](\mathrm{1}-{\tau }_i)(a_i+b_i{T}_a).

The emissivity {\epsilon }_i and atmospheric transmittance {\tau }_i can be estimated by the methods proposed above. For convenience, the coefficients of Eq. (23) are simplified as follows:

\{\begin{array}{l}{A}_i=b_i{\tau }_i{\epsilon }_i\\ B_i=a_i+b_i{T}_i-a_i{\tau }_i{\epsilon }_i\\ C_i=(\mathrm{1}-{\tau }_i)[\mathrm{1}+(\mathrm{1}-{\epsilon }_i){\tau }_i]b_i\\ D_i=-(\mathrm{1}-{\tau }_i)[\mathrm{1}+(\mathrm{1}-{\epsilon }_i){\tau }_i]a_i\end{array}.

Finally, the LST can be obtained through:

(24) T_s={\displaystyle \frac{C_{\mathrm{11}}(B_{\mathrm{10}}+D_{\mathrm{10}})-C_{\mathrm{10}}(B_{\mathrm{11}}+D_{\mathrm{11}})}{A_{\mathrm{10}}{C}_{\mathrm{11}}-A_{\mathrm{11}}{C}_{\mathrm{10}}}}.

Parameter adjustment for Qin’s algorithm

Qin’s algorithm is a popular split-window algorithm for AVHRR data, and this algorithm can be modified for use with Landsat 8 data (Rozenstein et al., 2014). The general form is:

(25) T_s=A_{\mathrm{0}}+A_{\mathrm{1}}{T}_{\mathrm{10}}-A_{\mathrm{2}}{T}_{\mathrm{11}},

where T_s is the LST, and T_{\mathrm{10}} and T_{\mathrm{11}} are the brightness temperatures of TIRS bands 10 and 11, respectively. A_{\mathrm{0}} , A_{\mathrm{1}} , and A_{\mathrm{2}} are the coefficients determined by the atmospheric transmittance and emissivities in both TIRS bands, and the related parameters are defined as:

(26) \{\begin{array}{l}{A}_{\mathrm{0}}=E_{\mathrm{1}}{a}_{\mathrm{10}}+E_{\mathrm{2}}{a}_{\mathrm{11}}\\ A_{\mathrm{1}}=\mathrm{1}+A+E_{\mathrm{1}}{b}_{\mathrm{10}}\\ A_{\mathrm{2}}=A+E_{\mathrm{2}}{b}_{\mathrm{11}}\\ A=D_{\mathrm{10}}/E_{\mathrm{0}}\\ E_{\mathrm{2}}=D_{\mathrm{10}}(\mathrm{1}-C_{\mathrm{11}}-D_{\mathrm{11}})/E_{\mathrm{0}}\\ E_{\mathrm{1}}=D_{\mathrm{11}}(\mathrm{1}-C_{\mathrm{10}}-D_{\mathrm{10}})/E_{\mathrm{0}}\\ E_{\mathrm{0}}=D_{\mathrm{11}}{C}_{\mathrm{10}}-D_{\mathrm{10}}{C}_{\mathrm{11}}\\ C_i={\epsilon }_i{\tau }_i\\ D_i=(\mathrm{1}-{\tau }_i)[\mathrm{1}+(\mathrm{1}-{\epsilon }_i){\tau }_i]\end{array},

where {\epsilon }_i and {\tau }_i are the same as in Mao’s algorithm, and the regression coefficients we chose were based on the condition of a temperature range of 10–50°C. Thus, a_{\mathrm{10}}=-\mathrm{64.6081} , b_{\mathrm{10}}=\mathrm{0.4399} , a_{\mathrm{11}}=-\mathrm{69.0215} , and b_{\mathrm{11}}=\mathrm{0.4756} .

Sensitivity analysis

The modified split-window algorithm for Landsat 8 requires both the land surface emissivity and atmospheric water vapor to retrieve the LST. In order to analyze the impact of the possible estimation errors for these critical parameters on the errors of the inverted LSTs, a sensitivity analysis is necessary. The following equation ( Qin et al., 2001a) is used to express the possible LST estimation error:

(27) \delta T=|T_s(x+\delta x)-T_s(x)|,

where \delta T is the error of the LST, x is the variable on which the sensitivity analysis needs to be performed, and \delta x is the possible error of this variable. T_s(x+\delta x) and T_s(x) indicate the LST obtained with the split-window algorithm when different values of x+\delta x and x are considered.

Sensitivity analysis for land surface emissivity

To analyze the sensitivity of the land surface emissivity for LST, we made an assumption that the value of atmospheric water vapor was 1.5 g/cm² in mid-latitude summer, and the transmittances for the TIRS bands were estimated by MODTRAN 4.0. Furthermore, due to the brightness temperature in band 11 always being less than that in band 10, we assumed it was 1.6 K smaller than band 10 in practice. Figure 3 shows the errors in the LST as a result of the errors in the land surface emissivity of band 11, with the emissivity of band 10 increasing. As shown in Fig. 3, we set the values of the emissivity error as 0.001, 0.002, 0.003, 0.005, and 0.01. As the emissivity error increases, the LST errors sharply increase. The curves indicate that the LST errors decrease as the emissivity increases. In summary, the LST has an error of 0.9°C with a 0.01 emissivity error in band 11.

Figure 3 displays the relationship between emissivity and LST errors at a brightness temperature of 300 K in band 10. As the brightness temperature increases, the errors in LST for band 11 are shown in Fig. 4. Clearly, a higher emissivity error produces a higher LST error, and the curves clearly indicate that the LST error increases as the brightness temperature increases. In this case, an 0.01 emissivity error at 50°C results in an error of more than 1°C in LST.

Because of the limited band numbers in Landsat 8, we chose the NDVI threshold method to retrieve the emissivity of band 10 and band 11. Thus, the emissivity errors will occur in both band 10 and band 11. Figure 5 shows the errors in the LST due to the errors in the land surface emissivity for band 10 and band 11. The results indicate a similar trend to Fig. 3, but the LST errors are lower than the emissivity error in the single band in Fig. 3. This shows that the LST error reaches 0.6°C with a 0.01 emissivity error in both band 10 and band 11.

Compared with Fig. 4, the LST errors in Fig. 6 are caused by the emissivity errors of both band 10 and band 11. Clearly, the curve trends are quite similar to Fig. 4, while the values of the LST errors are lower than those in Fig. 4. A 0.01 emissivity error at 50°C causes a less than 0.8°C error in the LST, as shown in Fig. 6.

Sensitivity analysis for atmospheric transmittance

In order to analyze the sensitivity of the atmospheric transmittance for the LST, we made an assumption that the values of the land surface emissivities of band 10 and band 11 were 0.970, and the at-sensor radiance of band 10 and band 11 could be computed from Eq. (1), while the value of the brightness temperature was 300 K in band 10. Figure 7 shows the errors in the LST as a result of the errors in the atmospheric transmittance of band 11 in the four possible conditions. The effects of the atmospheric transmittance for the LST are quite small when the atmospheric transmittance of band 11 is under 0.8, while the LST error increases when it is greater than 0.9.

Since the atmospheric transmittances of the two bands are calculated by the same atmospheric water vapor value, the errors from atmospheric water vapor have effects on both band 10 and band 11. Thus, we analyzed the conditions when errors in the atmospheric transmittance occur in both band 10 and band 11, as shown in Fig. 8, which plots the LST estimation error against the transmittance of band 10 for the four possible atmospheric transmittance errors in the two bands. The curves indicate that the effects of the LST errors are quite small when the atmospheric transmittance is less than 0.85, while the effects are greatly increased when it is greater than 0.85. The LST errors are less than 0.4°C when the errors in atmospheric transmittance are less than 0.02. In this case, when the atmospheric transmittance in band 10 is greater than 0.94, and the error reaches 0.05, the errors in LST are greater than 0.8°C. Since the atmospheric transmittance error of 0.05 corresponds to an error in atmospheric water vapor of 0.8 g/cm², and an error in the atmospheric water vapor of 0.1 g/cm² only causes a 0.2°C error in the LST, the errors in the LST as a result of the atmospheric transmittance in band 10 and band 11 can be kept within 0.5°C.

Validation of the algorithm

In order to validate the effectiveness of the modified split-window algorithm, it was used to retrieve the LST from Landsat 8 images from the city of Zoucheng in Shandong Province, China.

Data preprocessing

Two Landsat 8 images were chosen. One image was obtained on July 24, 2013; the other image was obtained on September 26, 2013. Firstly, the land surface was classified into four types: vegetation, bare land, buildings, and water. Figure 9(a) is the false-color composite image of the city of Zoucheng with Landsat 8 data. The data were recorded on July 24, 2013, and the atmospheric water vapor was 1.8 g/cm², which was retrieved by Landsat 8 Operational Land Imager (OLI) data using the FLAASH module in ENVI software. The image was classified by support vector machine (SVM), as shown in Fig. 9(b), and the classification accuracy was 97.8%.

Based on Eqs. (13) and (14), the land surface emissivities of band 10 and band 11 are shown in Fig. 9(c) and Fig. 9(d). The distributions of the emissivity are quite similar to the different land-surface types. It was found that the emissivity of water was greater than that of the other land-surface types, and the emissivity of soil and buildings was lower in these images.

For the other validation experiment, the second image was taken on September 26, 2013, and the atmospheric water vapor was 1.2 g/cm². Since the farmland was not covered with crops, the vegetation regions were divided into two parts: grass regions and farmland. The false-color composite image and the classification image are respectively shown in Fig. 10(a) and Fig. 10(b). The accuracy of the classification was 98.1%.

The land surface emissivities of band 10 and band 11 are respectively shown in Fig. 10(c) and Fig. 10(d). The values of the emissivities were much lower, in contrast to the results from July 24, especially in the farmland regions.

Experimental results and analyses

In order to analyze the accuracy of the proposed method for retrieving LST, we utilized a GSN at point C in Fig. 9(a), which consisted of GSN sensors at the locations shown in Fig. 11. Figure 12 shows a part of the GSN used in this fieldwork. The GSN could work both day and night, so it was used to obtain the temperature values when Landsat 8 was passing by, and the accuracy of the sensors was±0.5°C.

After obtaining the land surface emissivity and atmospheric transmittance of band 10 and band 11, we retrieved the LST for the city of Zoucheng by the use of the proposed split-window algorithm, and the results of July 24, 2013, are shown in Fig. 13(a). Figure 13(b) and Figure 13(d) show the LST results retrieved by Mao’s algorithm and Qin’s algorithm, respectively. The difference image shown in Fig. 13(c) illustrates that the water and the building regions have relatively large differences in the two algorithm results, while Fig. 13(e) indicates that the LST from Qin’s algorithm is lower than that from the proposed algorithm, and the maximum difference in the LST reaches 1.3°C for the buildings. However, it is confirmed that the ranges of the LST are quite close, and the heat island phenomenon of Zoucheng is also quite obvious: the high-temperature regions are buildings and bare land areas, while the lakes and rivers have a lower temperature.

Figure 14 displays the change in the vegetation temperature on July 24, 2013, from the wireless sensor that obtained temperature data every minute. In view of the passing-by time of Landsat 8 being 10:50 am, there are 100 temperature values from 10:00 am to 11:40 am. At the same time, the water temperature values were obtained from point A and point B in Fig. 9(a), the exact positions of which were located by GPS. High-precision thermometers were used, and the accuracy of the thermometers was 0.2°C. The results are shown in Table 4.

The LST results for September 26, 2013, are shown in Fig. 15(a), where it can be seen that the distribution of the LST is different to Fig. 13(a), due to the high temperature on farmland and the low temperature on buildings. The high temperature on farmland was caused by the decaying crops after reaping, while the low temperature in building regions was due to the low air temperature in September. Compared to the LST retrieved by Mao’s algorithm in Fig. 15(b) and Qin’s algorithm in Fig.15(d), the changing trends of the different land-surface types in the results of the proposed algorithm are quite similar. However, the water and the building regions show greater differences in Fig. 15(c), where the difference in the water regions is about 0.9°C, and the difference in the building regions is about -0.7°C. The greatest difference in Fig. 15(e) is about 1°C in the building and farmland regions, and the results indicate that the LST of Qin’s algorithm is lower than the LST of the proposed algorithm in the regions with low emissivity.

Figure 16 displays the temperature values on September 26, 2013, from 10:00 am to 11:40 am, which were obtained from the GSN every four minutes. The results from the ground survey and LST image are also shown in Table 4.

Table 4 shows the results of the LST for July 24 and September 26 retrieved from the remote sensing images and as measured by the ground temperature sensors. The results indicate that the difference values between the ground-measured values and the retrieved results are about 1°C. The LST errors for water are much greater than for the other types of land-surface features, and the errors for buildings and vegetation are less than 0.8°C. These errors are within the generally accepted level of 1.5°C ( Sibo et al., 2008; Zhou et al., 2011). Furthermore, the largest difference value with the ground-measured value for the proposed split-window algorithm is 1.144°C, which is higher than the 1.078°C for Qin’s algorithm, but is less than the 1.282°C for Mao’s algorithm. However, the difference values of water and buildings based on the proposed split-window algorithm are -0.617°C and -0.702°C, which are greater than those for Mao’s algorithm. Meanwhile, the difference value of buildings for Qin’s algorithm is 0.901°C, which is significantly larger than for the other two algorithms. The reason for this is that the ground temperature values obtained from the GSN were not the true values, and the results of the LST from these two algorithms are different, especially in the water regions, which have the greatest errors of nearly 0.75°C in Fig. 13(c) and 0.9°C in Fig. 15(c). Comparing the LST in Fig. 13(a) and Fig. 15(a), it can be seen that the distributions of the LST retrieved from these two algorithms are very similar. Furthermore, the air temperature from the meteorological station was 25–34°C for July 24 and 12–25°C for September 26. The average values of the GSN data and LST data are quite reasonable for the air temperature values because the air temperature is usually lower than the land surface temperature in daytime. Thus, it can be concluded that the estimated LST by the use of the proposed split-window algorithm is stable and close to the true value.

Discussion and conclusions

Many factors affect the retrieval of LST from remotely sensed thermal infrared images, and they can be mainly categorized into two aspects: atmospheric effects and land surface effects. The land-surface effects are caused by the different properties of the land-surface types, and thus different land-surface types have different emissivities. Although retrieving emissivity by NDVI is an effective approach, the relationship between the emissivity of different kinds of land cover and NDVI is still a challenging topic. In addition, due to the features of the wavelength ranges and spectral functions, different sensors have different values for the same land-surface type, but most of the published papers do not discuss the variation in emissivities for different sensors. We modified some parameters in Qin’s model and recalculated the typical land surface emissivities of vegetation, bare land, buildings, and water for Landsat 8. Although the changes are quite small when compared to Landsat 5, both band 10 and band 11 are slightly different. Atmospheric effects are the radiance effects of the atmosphere, which consist of absorption, reflection, and scattering by the atmosphere. In general, we can divide the atmospheric effects into the upwelling atmospheric radiance and the downwelling atmospheric radiance, which can be estimated by atmospheric transmittance and a standard atmospheric model. However, the estimated results depend on whether the standard atmospheric model is close to the true atmospheric conditions, otherwise the errors from the atmospheric parts will increase.

The proposed split-window algorithm for the retrieval of LST was specifically designed for Landsat 8 data. Based on the relationship of the Planck function between band 10 and band 11, we devised a novel linear equation with a coefficient of determination of about 0.9999 in different temperature ranges. As an additional equation for calculating the LST in the radiative transfer equation, this novel linear equation shows an effective correlation, and the coefficient of determination values are larger than for other models simplifying Planck’s radiation function (Rozenstein et al., 2014; Yang et al., 2014). In order to analyze the impact of the possible estimation errors for these critical parameters on the errors of the LST, sensitivity analyses for land surface emissivity and atmospheric transmittance were undertaken. In the sensitivity analysis for land surface emissivity, the results indicated that the LST errors in one band were greater than those in two bands. Furthermore, a 0.01 emissivity error in band 11 caused a 0.9°C LST error, and, in this case, a brightness temperature of 50°C would result in an error of more than 1°C. However, when a 0.01 error was given for both band 10 and band 11, it produced a 0.6°C LST error and a less than 0.8°C error at a brightness temperature of 50°C. Thus, when the errors of land surface emissivity are below 0.01 and the brightness temperature is less than 50°C, it is confirmed that the LST errors due to land surface emissivity are within 1°C. In the sensitivity analysis for atmospheric transmittance, the LST errors were caused by errors in the regression model of less than 0.25°C, and the errors in the LST as a result of the atmospheric transmittance in band 10 and band 11 were within 0.5°C. Overall, the LST errors caused by land surface emissivity and atmospheric transmittance were within 1.5°C. The sensitivity analysis demonstrated that this algorithm is very stable, and the LST error is within the generally accepted level of 1.5°C.

The validation of the algorithm was undertaken by experiments with data from the city of Zoucheng, and the results were compared with the results obtained from Mao’s and Qin’s split-window algorithms. The results of the LST displayed an apparent heat island effect in the city of Zoucheng for July 24, and a high temperature of farmland due to decaying crops was detected for September 26. The retrieved temperature values by the proposed algorithm in the experiments for the city of Zoucheng were very close to the ground-measured values obtained by the GSN, and the LST error was about 1°C.

References

[1] Becker F, Li Z L (1995). Surface temperature and emissivity at various scales: definition, measurement and related problems. Remote Sens Rev, 12(3‒4): 225–253

[2] Berk A, Anderson G, Acharya P, Chetwynd J, Bernstein L, Shettle E, Matthew M, Adler-Golden S (1999). MODTRAN4 user’s manual. Air Force Research Laboratory, Space Vehicles Directorate, 1–5

[3] Berk A, Bernstein L S, Robertson D C (1987). MODTRAN: a moderate resolution model for LOWTRAN. In: DTIC Document

[4] Du C, Ren H, Qin Q, Meng J, Li J ( 2014). Split-Window algorithm for estimating land surface temperature from Landsat 8 TIRS data. In: Geoscience and Remote Sensing Symposium (IGARSS), 2014 IEEE International, 3578–3581

[5] Harris A, Mason I (1992). An extension to the split-window technique giving improved atmospheric correction and total water vapour. Int J Remote Sens, 13(5): 881–892

[6] Jiménez-Muñoz J C ( 2003). A generalized single-channel method for retrieving land surface temperature from remote sensing data. J Geophys Res, 108,

[7] Jiménez-Muñoz J C, Cristóbal J, Sobrino J A, Sòria G, Ninyerola M, Pons X ( 2009). Revision of the single-channel algorithm for land surface temperature retrieval from Landsat thermal-infrared data. Geoscience and Remote Sensing. IEEE Transactions on, 47: 339–349

[8] Jimenez-Munoz J C, Sobrino J A, Skokovic D, Mattar C, Cristobal J ( 2014). Land surface temperature retrieval methods from Landsat-8 thermal infrared sensor data. IEEE Geosci Remote Sens Lett, 11(10): 1840–1843

[9] Kneizys F, Shettle E, Gallery W, Chetwynd J Jr, Abreu L (1983). Atmospheric transmittance/radiance: computer code LOWTRAN 6. Supplement: Program listings. In: DTIC Document

[10] Kneizys F X, Shettle E, Abreu L, Chetwynd J, Anderson G (1988). Users guide to LOWTRAN 7. In: DTIC Document

[11] Krishna A P, Sharma A ( 2013). Snow cover and land surface temperature assessment of Gangotri basin in the Indian Himalayan Region (IHR) using MODIS satellite data for climate change inferences. In, SPIE Remote Sensing (pp. 889314-889314-889319): International Society for Optics and Photonics

[12] Lu D, Weng Q (2006). Spectral mixture analysis of ASTER images for examining the relationship between urban thermal features and biophysical descriptors in Indianapolis, Indiana, USA. Remote Sens Environ, 104: 157–167

[13] Mao K B, Qin Z H, Shi J C, Gong P ( 2005a). The research of Split-Window Algorithm on the MODIS. Geomatics and Information Science of Wuhan University, 30(8): 703–707

[14] Mao K, Qin Z, Shi J, Gong P ( 2005b). A practical split-window algorithm for retrieving land-surface temperature from MODIS data. Int J Remote Sens, 26(15): 3181–3204

[15] Patel N, Parida B, Venus V, Saha S, Dadhwal V ( 2012). Analysis of agricultural drought using vegetation temperature condition index (VTCI) from Terra/MODIS satellite data. Environ Monit Assess, 184(12): 7153–7163

[16] Price J C (1983). Estimating surface temperatures from satellite thermal infrared data—A simple formulation for the atmospheric effect. Remote Sens Environ, 13(4): 353–361

[17] Price J C (1984). Land surface temperature measurements from the split window channels of the NOAA 7 advanced very high resolution radiometer. Journal of Geophysical Research: Atmospheres (1984–2012), 89: 7231–7237

[18] Qin Z, Dall’Olmo G, Karnieli A, Berliner P ( 2001a). Derivation of split window algorithm and its sensitivity analysis for retrieving land surface temperature from NOAA‐advanced very high resolution radiometer data. Journal of Geophysical Research: Atmospheres (1984–2012), 106: 22655–22670

[19] Qin Z, Karnieli A, Berliner P ( 2001b). A mono-window algorithm for retrieving land surface temperature from Landsat TM data and its application to the Israel-Egypt border region. Int J Remote Sens, 22(18): 3719–3746

[20] Qin Z, Li W, Gao M, Zhang H( 2006). Estimation of land surface emissivity for Landsat TM6 and its application to Lingxian Region in north China. Proc SPIE 6366, Remote Sensing for Environmental Monitoring, GIS Applications, and Geology VI, 636618

[21] Rozenstein O, Qin Z, Derimian Y, Karnieli A (2014). Derivation of land surface temperature for Landsat-8 TIRS using a split window algorithm. Sensors (Basel Switzerland), 14: 5768–5780

[22] Sasamori N (1999). Calculation of emissivity by measuring reflectance and transmittance. Report of Tokyo Metropolitan Industrial Technology Research Institute,2:45–48

[23] Schädlich S, Göttsche F, Olesen F S ( 2001). Influence of land surface parameters and atmosphere on METEOSAT brightness temperatures and generation of land surface temperature maps by temporally and spatially interpolating atmospheric correction. Remote Sens Environ, 75(1): 39–46

[24] Sibo D, Guangjian Y, Yonggang Q ( 2008). Two single-channel algorithms for retrieving land surface temperature from the simulated HJ-1B data. Prog Nat Sci, 18: 1001–1008

[25] Sobrino J A, Jiménez-Muñoz J C, Paolini L ( 2004). Land surface temperature retrieval from LANDSAT TM 5. Remote Sens Environ, 90(4): 434–440

[26] Sobrino J A, Li Z L, Stoll M P (1993). Impact of the atmospheric transmittance and total water vapor content in the algorithms for estimating satellite sea surface temperatures. Geoscience and Remote Sensing. IEEE Transactions on, 31: 946–952

[27] Susskind J, Rosenfield J, Reuter D, Chahine M (1984). Remote sensing of weather and climate parameters from HIRS2/MSU on TIROS-N. Journal of Geophysical Research: Atmospheres (1984–2012), 89: 4677–4697

[28] Valor E, Caselles V (1996). Mapping land surface emissivity from NDVI: application to European, African, and South American areas. Remote Sens Environ, 57(3): 167–184

[29] Van de Griend A, Owe M (1993). On the relationship between thermal emissivity and the normalized difference vegetation index for natural surfaces. Int J Remote Sens, 14(6): 1119–1131

[30] Wan Z ( 2014). New refinements and validation of the collection-6 MODIS land-surface temperature/emissivity product. Remote Sens Environ, 140: 36–45

[31] Wan Z, Dozier J (1996). A generalized split-window algorithm for retrieving land-surface temperature from space. Geoscience and Remote Sensing. IEEE Transactions on, 34: 892–905

[32] Weng Q, Lu D ( 2008). A sub-pixel analysis of urbanization effect on land surface temperature and its interplay with impervious surface and vegetation coverage in Indianapolis, United States. Int J Appl Earth Obs Geoinf, 10(1): 68–83

[33] Wukelic G E, Gibbons D E, Martucci L M, Foote H P (1989). Radiometric calibration of Landsat Thematic Mapper thermal band. Remote Sens Environ, 28: 339–347

[34] Yang L, Cao Y, Zhu X, Zeng S, Yang G, He J, Yang X ( 2014). Land surface temperature retrieval for arid regions based on Landsat-8 TIRS data: a case study in Shihezi, Northwest China. Journal of Arid Land, 6(6): 704–716

[35] Yong Z, Yu T, Xingfa G, Yuxiang Z, Fuliang Z, Shanshan Y, Wenjun Z, Xiaowen L ( 2006). Land surface temperature retrieval from CBERS-02 IRM SS thermal infrared data and its applications in quantitative analysis of urban heat island effect. Journal of Remote Sensing, 10: 789–797

[36] Zhou J, Li J, Zhao X, Zhan W F, Guo J X ( 2011). A modified single-channel algorithm for land surface temperature retrieval from HJ-1 B satellite data. Journal of Infrared and Millimeter Waves, 30(1): 61–67

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg