LEAST SQUARE APPROACH FOR ESTIMATING OF LAND SURFACE TEMPERATURE FROM LANDSAT-8 SATELLITE DATA USING RADIATIVE TRANSFER EQUATION

Land Surface Temperature (LST) is one of the significant variables measured by remotely sensed data, and it is applied in many environmental and Geoscience studies. The main aim of this study is to develop an algorithm to retrieve the LST from Landsat-8 satellite data using Radiative Transfer Equation (RTE). However, LST can be retrieved from RTE, but, since the RTE has two unknown parameters including LST and surface emissivity, estimating LST from RTE is an under the determined problem. In this study, in order to solve this problem, an approach is proposed an equation set includes two RTE based on Landsat-8 thermal bands (i.e.: band 10 and 11) and two additional equations based on the relation between the Normalized Difference Vegetation Index (NDVI) and emissivity of Landsat-8 thermal bands by using simulated data for Landsat-8 bands. The iterative least square approach was used for solving the equation set. The LST derived from proposed algorithm is evaluated by the simulated dataset, built up by MODTRAN. The result shows the Root Mean Squared Error (RMSE) is less than 1.18°K. Therefore; the proposed algorithm can be a suitable and robust method to retrieve the LST from Landsat-8 satellite data.


INTRODUCTION
Land Surface Temperature (LST) is a significant parameter in climatic, hydrologic, ecological, biogeochemical, and related studies (Wan 1999).LST is a major factor in global change studies, estimating radiation budgets in heat balance studies and as a control index for climate models (Rozenstein et al. 2014).Furthermore, accurate calculation of outgoing long wave radiation emitted from the earth's surface needs accurate LST and emissivity values (Liang 2005).Up to now, several techniques have been proposed to estimate LST from satellite passive sensors data.Typical of these methods is a single channel and split window.The single channel technique is based on the linearization of the radiative transfer equation.The split window technique is estimated LST by the use of brightness temperature of two TIR bands which are located in the atmospheric window between 10 and 12 μm.This technique is based on the correcting thermal observations distorted by atmospheric effects, utilizing the differential absorption in two adjacent TIR channels.In these techniques, Land Surface Emissivity (LSE) is assumed as known parameter.However, estimation of the accurate LSE from thermal radiance measurements can be complicated.Estimating LST and LSE from radiative transfer equation is an under determined problem because N thermal bands have N observation (thermal radiance) with N+1 unknown (N spectral emissivities and one LST).Up to now, many techniques are developed for estimating LST and LSE from different remote sensors.Most of these techniques are a review in two publications (Li, Tang, et al. 2013;Li, Wu, et al. 2013).These methods can be categorized into two groups: 1) estimating LST by assuming LSE as a known parameter 2) Estimating LSE and LST by solving an undetermined problem some of these techniques includes MODIS Day/Night LST algorithm which has been proposed by the MODIS Science Team based on assumptions that LSE ratios are the same or do not change significantly between two times (Wan 1999), MODIS CEBLST method developed by Momeni and Saradjian for MODIS data by adding two constraint equations based on relation between LSE and NDVI (Momeni and Saradjian 2008) and ASTER TES algorithm which has been proposed by the ASTER Science Team based on Relationship between the minimum LSE and spectral contrast (Gillespie et al. 1998).Most recent studies for LST retrieval from Landsat 8 TIRS are by (Rozenstein et al. 2014;Jiménez-Muñoz et al. 2014;Moghaddam, Akhoondzadeh, and Saradjian 2015a;Tan et al. 2017;Yu, Guo, and Wu 2014;Wang et al. 2015).In this study, we propose a method for LST retrieval from the Landsat-8 TIR Data based on the equation set of RTE for band 10 and 11 Landsat-8 and solving the underdetermined problem by adding new equations.Indeed, the four equations were made for this proposition.Two of them were developed based on the RTE for band 10 and 11, and two of them were developed based on the relation between NDVI and emissivity of two thermal bands of Landsat-8.The report issued by U.S. Geological Survey expressed a stray light caused calibration problem of TIRS bands.The Landsat-8 team is trying to puzzle out the issue, and planned reprocessing of Landsat-8 imagery.Thus, in this article, algorithms were tested by an extensive simulated data set

MATERIALS AND METHODS
Two simulated datasets were built up in this study, one of them (SD1) based on Plank's law another one (SD2) based on MODTRAN 4. These datasets (SD1 and SD2) includes brightness temperatures, surface temperature, emissivity and atmospheric parameters (atmospheric transmission, upwelling, and downwelling radiance) for the TIRS bands.Using equation 1 reflectance of Landsat-8 bands have been simulated for different material based on two separate spectral library data and relative spectral response of Landsat-8 reflective and thermal wavelengths.
Where R is the reflectance of simulated data in the range of (λmin, λmax), F is the spectral response function, and R is the real reflectance.
ASTER Spectral Library (ASL, http://speclib.jpl.nasa.gov)was used to build up the simulated dataset.ASL contains directional hemispherical reflectance of the different type's area.Table I provides information about number, types, and the name of the reference library of the simulated data which are built up.Table II shows these datasets characteristics.
The method used in this study includes two observed radiance equations for bands 10 and 11 based on the RTE and two additional equations based on the relation between the NDVI and LSE of TIRs band which is developed by (Moghaddam, Akhoondzadeh, and Saradjian 2015b).Figure 1 shows the flowchart of this algorithm.Equation 2 can be written in the matrix form equation 3: (3) Where L is observation matrix, and its components are L10, L11, f () and g ().X is an unknown matrix that includes LST and two emissivities.Since the equation 3 is nonlinear, we used simple linearized by using Taylor's first approximation series.Therefore equation 3 can be written as: (4) The Jacobian matrix can be written as: (5) (6) An iterative weighted least squares solution can solve X. Equation 7explains this solution: (7)

Standard atmospheric profiles
In this study, the Q matrix was estimated based on Planck's radiation law and Landsat-8 spectral responses of the bands 10 and 11 by using equation 8: Where B (T, λ) is the Planck's law, and ψ(λ) is the Landsat-8 spectral response.The λ1 and λ2 are the spectral range of Landsat-8 bands 10 or 11.
Therefore, the initial value for equation 7 are suggested below:

EXPERIMENTAL RESULTS
The Proposed algorithm was applied to the SD1 and SD2 test cases.The static parameters that were used for analysis of error include the mean error (bias), standard deviation error and RMS were derived.These values are provided in Table III.
Table III depicts that, the RMSEs of LST for simulated data that were built up by using MODTRAN are around 1.18°K and for SD1 that were created by using Plank's law, RMSE is equal to 0.99°K.Since the result shows that the mean temperature error is almost close to zero (i.e.: LSTs which are estimated by using the proposed algorithm, are close to the simulated LST), the proposed algorithm is numerically free of bias.The results show that the maximum of mean error, standard deviation error, and RMSE of SD2 are 0.32°K, 1.23°K, and 1.24°K respectively, while their minimums are 0.10°K, 1.09°K, and 1.12°K respectively.
The RMSEs of LSE for SD2, are around 0.0063 for band 10 and 0.0066 for band 11, and RMSEs of LSE based on the SD1 data set is equal 0.0065 for band 10 and 0.0069 for band 11.The land surface emissivity uncertainty, which leads to an error on the LST.

CONCLUSION
In this study, we have presented an algorithm to retrieve LST from Landsat-8 satellite data using equation sets based on the RTE and additional constraint equations.The least square approach was used to solve.The relation regression between the NDVI and LSE based on the spectral library was used to create the additional constraint equations.The algorithm presented was tested with simulated data sets which were built up based on different atmospheric conditions and for a variety of land cover such as bare soil, rock, vegetated, and so on.The main advantage of the proposed method is retrieving LST from RTE without any approximation and auxiliary data.Also, this method is solving the emissivity and LST simultaneously.
in the SD2, Landsat-8 TIRS band (i.e.: band10 and band11) brightness temperatures are obtained from the RTE by inversion of the Planck's law, surface temperature was chosen based on the temperature of the first layer of the atmospheric profiles (T0) as T0 − 5°K, T0, T0 + 5°K, T0 + 10°K, and T0 + 20°K (Jiménez-Muñoz and Sobrino 2008), the emissivity was extracted from spectral library that explained in table I and the atmospheric parameters had been simulated using the MODTRAN-4 based on the standard atmospheric profiles of MODTRAN (includes: tropical (TRO), mid-latitude summer (MLS), mid-latitude winter (MLW), subarctic summer (SAS), U.S standard (USS), sub-arctic winter (SAW)) and the spectral response functions of Landsat-8 TIRS.

Figure 1 -
Figure 1 -Flowchart of NDVI-based emissivity estimation for Landsat-8 TIRS bands Therefore, the equation set can be written as equation 2 for weighted least squares solution

Table I -
property of simulated data

Table II -
Simulated datasets characteristics

Table III .
Result of the validation using simulated data