LAND SURFACE TEMPERATURE ESTIMATION AT THE GLOBAL SCALE USING SATELLITE OBSERVATIONS

Aim of this paper is to offer a brief overview of satellite based LST estimation fundamentals, discuss the main design principles of the corresponding algorithms, present some examples of the global scale LST and land surface emissivity products from the literature that were derived using satellite observations, and finally highlight different challenges in the described measurement process.

temperatures are determined by the Wien's displacement law, see Section 2.1), or in the: b) mid-infrared (MIR) range (3-8 µm), which is more sensitive to higher temperatures (693 ÷ 89 °C), like the ones resulting from the forest fires. Sensor design is also influenced by the atmospheric opacity, i.e. the spectral range of electromagnetic waves that pass through the atmosphere without significant absorption under the clear sky conditions, since gases in the atmosphere have different absorptions across the spectrum. Generally speaking, radiation transmissivity is always imperfect, and these effects also need to be taken into account in LST estimation, for more details see Section 3.1. Topography of the terrain also can have some influence on the daytime measurements, since the slopes that are oriented towards the Sun are warmer than the slopes occluded by the shadows. The main motivation for LST estimation comes from its application in different energy-balance models, climate or geological studies. However, temperature estimation problem is considered as underdetermined and usually there is a need to incorporate additional measurements or information beside radiance measured by the main observation instrument in order to better constrain the set of possible solutions and improve the overall accuracy of the land surface temperature recovery. Approaches for satellite based LST retrieval can be categorized as procedures that estimate LST in a deterministic way, [1]; or on the other hand approaches that try to improve their estimate through regression based on some other quantities that are easier to measure accurately, e.g. vegetation or spectral indices in the visible or near infrared (VNIR) part of the spectrum, [2]. The rest of the paper is organized as follows: in Section 2 we introduce the physical quantities that are used in characterization of LST and that are measured by the satellite instruments. A relationship between a temperature of an idealized surface and its spectral radiance distributions is discussed in details. Section 3 starts with an overview of atmospheric influences and presents the main LST estimation approaches, and corresponding algorithms. In Section 4 are presented the main characteristics of the previously published and carefully validated Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global emissivity dataset product, [4], and its predecessor [5], which is followed by an example of singlechannel LST estimation using satellite image acquired by ASTER, [3]. Finally, in Section 5 we conclude the main exposition and point out future research directions.

PHYSICAL LAWS OF RADIATION AND SPECTRAL EMISSIVITY
Physical quantity that is related to surface temperature is the radiant flux or power R P [W], which corresponds to the total radiant energy per unit time that is emitted by some surface 0 A . However, measurement of the radiated energy generally depends on the relative position of the sensor that is observing some emitting surface. Therefore, a more sophisticated quantity that is usually used to describe surface radiation at a particular point in space is radiance R . It corresponds to the observed radiant power that is normalized per unit solid angle dΩ , and per unit projected area proj A ∂ that is defined by the observer's position, Fig. 1. Viewing direction angle θ corresponds to the offset between radiating surface normal and the radius vector defining observer's relative position. Radiance is expressed in units [W sr −1 m −2 ], and since in general it represents a directional quantity (depending on θ and azimuth angle φ ), it is usually denoted as R Ω , and by definition equal to: Note that the maximal solid angle Ω [sr] of 4π steradians corresponds to the area of a sphere with radius r that is divided by 2 r , i.e. the unit solid angle in spherical coordinate system is given by: On the other hand, any half sphere defines a solid angle of 2π sr. Cosine emission law, or more widely known as Lambert's emission law, states that the surface that can be considered as an ideally diffuse emitter (radiator) has a radiance that is independent of viewing direction θ , i.e. its radiance can be considered as isotropic. More exactly, radiant intensity: the Lambertian surface is proportional to the projection of emitting surface normal onto viewing direction θ , as illustrated by the term 0 cos I θ in Fig. 1, where 0 I represents its nominal radiant intensity in the vertical direction. Since radiance at some point, determined by the viewing direction θ , represents the corresponding radiant intensity in the direction θ that is normalized by the area of the observing surface, it follows that the radiance of the Lambertian surface at some point on the sphere is given by: which confirms that in such case Lambertian surface has the same radiance when viewed from any angle θ .
Radiant intensity of the Lambertian surface that is integrated over a half of the sphere gives the total emitted power of the radiant surface in the observer's hemisphere, and can be computed as: Expression (4) shows that the multiplicative or scaling factor between the radiant intensity of the Lambertian surface and the total power emitted into solid angle of 2π sr is equal to π . Similar relation also holds in the case of radiance 0 R in (3), which could be integrated over the same solid angle of 2π sr and in such case would give a radiant power or flux received by some hemisphere per unit area, which can be denoted by the symbol R However, since radiation is usually also spectrally varying property of the radiant surface, the quantity that is of more interest for accurate radiation measurements than the previously discussed radiance R Ω is the spectral radiance, When talking about effects related to reflection of radiation, a common quantity that is used in describing such process is the irradiance, I, or more exactly a radiant power received by a surface per unit area. Similar to radiance it also appears in the more sophisticated form of spectral irradiance, as well as in the directional form.

Planck's law, its approximation and the role of the temperature
In order to characterize radiation of some surface, its characteristics are compared to the radiation of the ideal emitting material, which is capable of absorbing all incident energy, but more importantly that is also capable of emitting all thermal energy possible. Such material model is also known as the blackbody and its radiation behaviour at different frequencies is described by one of the most important physical laws, Planck's law. It states that the spectral radiance of a blackbody at absolute temperature T, and frequency ν , is given by: where c is the speed of light, h is the Planck constant, and B k represents the Boltzmann constant, Fig. 2. Note that the spectral radiance in Fig. 2 is expressed in function of the wavelength, instead of frequency.
Surface of the blackbody in the thermodynamic equilibrium is considered to be Lambertian, which means that the spectral radiance is regarded as independent or uniform with respect to observer's angle of view. By integrating expression in (5) over all frequencies and over all solid angles corresponding to the half of the sphere (the other half is considered as not observable), total radiant power emitted per unit area of the black body is determined by the expression that is also known as the Štefan-Boltzmann law:  , which corresponds to the broadband measurement of surface radiation, by the measurement instrument that is not frequency selective, but aggregates emitted radiation that is coming from all directions and over all wavelengths. Such measurement would provide a way to measure the temperature of the blackbody by inversion of B R .
Distribution of the blackbody's spectral radiance that is defined by (5) depends on the temperature, as shown by the curves in Fig. 2-3. Radiance increase at some specific wavelength, which is a consequence of the rise in the blackbody's temperature, has an exponential character, as illustrated in Fig. 3. Also, the wavelength corresponding to the maximum of each spectral radiance distribution curve changes inversely with the change in temperature, i.e. with the rise of temperature, peaks of the curves move towards the shorter wavelengths, corresponding to higher frequencies (energies), Fig. 2. Described phenomenon is controlled by the Wien's displacement law and in a simple way indicates that the maximum, as well as the overall shape of the spectral radiance distribution, changes with temperature. The law states that for the blackbody, governed by (5), the maximum of the spectral radiance distribution is given by simple relation: where the displacement constant b is approximately equal to: . Such approximation is considered as more precise at the shorter wavelengths, as illustrated in Fig. 2.

Spectral emissivity
Effectiveness of some material in emitting energy as thermal radiation at some specific wavelength is described by its spectral emissivity, which represents a ratio of the material's surface spectral radiance and the spectral radiance of the blackbody that would be at the same temperature as the material and would emit radiation at the same wavelength: Spectral radiances in (9) are expressed with the arguments in the form of wavelength instead of frequency, similarly like in the Fig. 2, but more importantly, (9) represents the total hemispherical emissivity ratio of all spectral radiances observed over hemisphere H, i.e. the ratio of directional spectral radiances of the material and the corresponding spectral radiances of the blackbody at the same temperature and on the same wavelength, but that are integrated (averaged) over all possible directions of some half sphere.
On the other hand, it is also possible to define a more specific quantity, i.e. directional spectral emissivity, which in comparison to (9), besides wavelength, also depends on the specific direction of observation: and enables more detailed modelling. Emissivity of the surface cannot exceed 1, and it is fundamentally related to the surface's absorption capability. Namely, Kirchoff's law of thermal radiation states that for the body in thermodynamic equilibrium, meaning that a single temperature determines the current state of the system, like in the case of the blackbody, surface emissivity is equal to its absorptivity (input radiation is equal to output radiation). If it is assumed that the body is opaque and there is no transmissivity, based on the law of conservation of energy, it follows that the spectral reflectivity of the surface is given by the following relation:

TEMPERATURE INVERSION PROBLEM
Measured TOA radiance represents a superposition of emitted radiation originating from different types of sources. The most significant source for the purpose of LST estimation is the radiation emitted by the land surface, which in the relatively simplified model can be described by the spectral radiance of the blackbody at the given temperature T that is additionally attenuated by the coefficient describing real thermal radiation effectiveness of the given material in comparison to the idealized blackbody material with the emissivity equal to one, i.e. it can be described by the product of the material's spectral emissivity defined in (9) or (10) and the spectral radiance ( , ) B T λ defined in (5). Additional sources that can be indicated as having relatively significant contribution to the measured radiance value are factors that are taking into account influence of the atmosphere through which the originally emitted radiation propagates on its way to the TOA sensor.

Influence of the atmosphere on LST estimation and its modelling
These influences can be grouped into two main categories.
Namely, quantity denoted as the path radiance, L ↑ , is usually employed to describe the backscattering radiation that never reaches the Earth's surface, and that is reflected back to the TOA sensor by the particles and molecules in the atmosphere. It is caused by the Sun's radiation that is going through the upper layers of the atmosphere on its way to the Earth's surface. The second category of atmospheric influences is described by the quantity denoted as the sky irradiance, I ↓ , and accounts for the radiation effects originating from the radiation that is diffusely emitted by the atmosphere, that reaches the Earth's surface, and then finally reflects from the surface back to the sensor at TOA. Since the sky irradiance is diffusely reaching the Earth's surface from all possible directions in the hemisphere, it represents averaged or integrated down-welling radiance L ↓ of the atmosphere, i.e. it is equal to the L ↓ that is integrated over all directions corresponding to the hemisphere above the Earth's surface.
Since radiance L ↓ is considered as diffuse (same average It should be mentioned that in the most general consideration, all introduced radiance, irradiance, and material properties are considered as spectrally and directionally selective, however subscript Ω denoting directional dependence is usually omitted, in order to ease the notation or more importantly to simplify the general model by considering the corresponding quantities as undirected. However, such simplification is not always possible. Having in mind that emitted and reflected radiation from the surface are attenuated on their path through the atmosphere, which is described by the spectrally selective transmissivity of the atmosphere, λ τ , it follows that the observed TOA radiance ( , ) L T λ can be described as: LST term, T, appears in (12) as the argument of ( , ) B T λ , given by (5), and could be computed directly if estimates of other quantities would be available. Formally, there are five unknowns, out of which four are related to specific spectral properties of the material and atmospheric conditions. Therefore, LST estimation is considered as underdetermined problem, and there are various approaches that offer different solutions. Since satellite observations are usually made by multispectral instruments that are capable of providing radiance measurements at different wavelengths, it is possible to have a system of multiple equations of the form presented in (12), however, due to spectral sensitivity of described quantities there is always a problem of having an underdetermined system with four unknowns per each additional equation. Nevertheless, there exist several operational LST products that overcome this challenges and provide satisfactory LST estimates at the global scale using satellite observations. In the following few lines we will briefly discuss some of the challenges and approaches in constraining the solution set of (12). In addition, it should be also mentioned that the LST estimate obtained from (12) actually represents an average land surface temperature, which is spatially averaged over the spatial region represented by the single image pixel in the measurement model, but which in reality represents a heterogeneous surface with spatially varying properties. Similarly, it should be clear that the LST estimate depends on the characteristics of the land cover, acquisition time (day or night), shadows, spatially varying angle of view under which instrument observes different pixels in the same scene, as well as terrain topography. In order to determine characteristics of the atmosphere and compensate for the effects produced by the non-ideal transmissivity, presence of the additional path radiance, and the sky irradiance, often practice is to utilize independently observed atmospheric profiles of humidity, pressure and atmosphere temperature. This estimation procedure is usually performed by some suitably developed atmospheric radiative transfer model, like the MODTRAN, [6], which was e.g. utilized in one of the first global scale LST products, [7]. Such models usually also require additional information about the slope of the Earth's surface (the corresponding image pixel) and its height above the sea level, in order to determine internal parameters of the radiative transfer model that depend on the relative position of the Sun, as well as to incorporate diffuse Sun irradiance components and reflections from the ground in the case that the surface of the slope is not zero. Therefore, in order to achieve higher accuracy, estimation of these unknown quantities also relies on the quality of the utilized digital elevation model. All additional observations and measurements beside the main TOA radiance are usually performed by some other instruments that are either: a) on board of the same platform as the main TIR sensor, and which often operate in the shortwave infrared (SWIR) or VNIR spectral range; or b) on board of some other, specially designed probes (like the atmospheric balloons) or platforms that are acquiring observations from the same area in the relatively small temporal window around the radiance measurements performed by the main TIR instrument. If emissivity of the Earth's surface is somehow known or reliably estimated in advance, there are some approaches that utilize such special circumstances and perform estimation of required atmospheric radiances from (12) using only thermal infrared measurements, by exploiting differential absorption in adjacent TIR bands. Such methods are also known as the "split-window" algorithms, [8][9], and rely on knowledge base of emissivities of natural terrestrial materials, or global and regional emissivity products like the ones described in [3][4] and [10].

Effective wavelength of the measuring instrument
Usual assumption in LST estimation is that the corresponding LWIR and MIR sensors are designed as narrowband observation instruments, however sometimes this is only a partially true. Therefore, in order to achieve higher accuracy, there is a practice to utilize additional information preserved in the spectral response (sensitivity) of the instrument, which has been measured in advance under controlled, calibration conditions. A simplified procedure usually consists of determining effective wavelength of the instrument (filter), after which it can be considered as a narrowband instrument with an effective wavelength eff λ . An overview of different methods that are utilized in LST estimation was presented in [11]. It can be said that the general idea of the "effective value" is to find the weighted average of some quantity over the spectrum range of interest. By such approach original filter could be replaced by the idealized filter with simplified rectangular shape of width W that is idealistically performing measurements only at a single wavelength eff λ , but which aims to preserve the total power received by the real sensor at different wavelengths.
The right hand side of (13) will be equal to the desired approximation: where it follows that the effective wavelength is equal to: which can be interpreted as the expected value of λ over the distribution defined by the filter's spectral response function. In general, such "averaging" approach using T λ is favorable if we strive towards the higher accuracy, or if the sensor cannot be considered as narrowband.

Basic approaches in LST estimation
Generally speaking, computation of temperature through inversion of measured radiance can be performed directly from (5). In order to simplify the notation, introduce the following constants: As previously discussed, assuming that the instrument can be considered as narrowband, and by adopting simplifying assumptions that the measured ( , ) L T L λ  corresponds exactly to the spectral radiance ( , ) B T λ λ of the blackbody, which includes negligence of all atmospheric effects, as well as the real material emissivity resulting from the land cover characteristics, LST is obtained as: In the case when atmospheric parameters L ↓ , L ↑ , λ τ , from (12) can be estimated for some effective wavelength λ , as well as viewing direction angle θ , LST estimate can be improved by using the more general temperature inversion expression: Usual approach to significantly improve simple LST estimate in (16) is to incorporate some estimate of surface's spectral emissivity λ ε , i.e. use: . This can be achieved by using some global emissivity dataset like [4], or by estimating λ ε using some adequately designed method, e.g. algorithms based on measurements of normalized-difference vegetation index (NDVI) or fractional vegetation cover (FVC), [2], [13]. In the case that there is no vegetation cover, the emissivity of the bare land corresponds to soil emissivity and can be estimated from the empirical relationship of emissivity with the measurements in the visible red channel, [12].
Assuming that the emissivity estimate is known, and that atmospheric parameters are estimated by some other instrument or method, presented approaches can be considered as the "single-channel" algorithms, since only one TIR channel is needed to estimate LST, like e.g. the single-channel algorithm proposed in [13] that is utilizing information about water vapour content in the atmosphere. However, depending on the type of the surface, targeted temperature range, or atmospheric conditions, additional multispectral measurements over other wavelengths could still be required.
Another type of algorithms for LST retrieval are those relying on some combination of measurements in two or more TIR image channels. In the case of only two acquisition channels, like e.g. in [14], or [8], such approaches are named as "split-window" algorithms, [12]. Their main idea is to improve the final LST estimate s T by expressing it as a linear combination of individual LST estimates, brightness temperatures (16), that are obtained independently from each of the observation channels.
This means that measured radiances at specific wavelengths are used separately to compute each of the individual temperature estimates 1 T and 2 T , however question arises what is with other quantities that appear in the general measurement model (12). Parameters that still will be needed per each channel are the atmospheric transmissivity λ τ , and the surface emissivity λ ε , however the path radiance L ↑ and the sky irradiance L ↓ in (12)  T . This assu-mption is essential, since it enables one to replace L ↑ and L ↓ in (12) by the equations of type (15), putting the emphasis of estimation process on approximation of different radiance functions (15)  T and a T that appear in the system of two equations of type (12). Only requirement is to somehow remove derivatives that appear in all of the described first order expansions, and exactly this is achieved by the precomputed coefficients that enable linear fitting of radiance values against the temperatures in some predefined range of interest. Fitting coefficients are precomputed in advance by independent simulations corresponding to each of the observation channels.
Another, alternative and simpler approach to LST estimation that is also based on two measurement channels, is the one that is utilizing approximation of the Planck's law described in Section 2.1. It is based on the variations in shape of the radiance distribution curves that are described by (15) and illustrated in Fig. 1 It should be mentioned that in (18) it was assumed that: Since there is a logarithm of the ratio in (19), a usual practice in LST estimation based on this approach is to choose pairs of TIR wavelengths on larger mutual distance. As compared to the most of the presented design approaches, some recently proposed methods, like [15], try to overcome the requirement for accurate information about the properties of the atmosphere and the surface emissivity by directly modelling dependence between measured radiance and unknown LST using nonlinear model with several artificial parameters that have the role to aggregate original unknown parameters that were involved in the measurement process. Idea is to find enough similar, spatially adjacent pixels, where the corresponding similarity is determined based on some other multispectral observations independent of TIR measurements, and then numerically solve the corresponding system of nonlinear equations corresponding to several selected TIR measurements. In such a way, underdetermined nature of the original LST estimation problem hopefully can be resolved.

EXAMPLES OF SATELLITE BASED LST MEASUREMENTS
Operational LST products that are available at the global scale are of significant importance for different studies, as well as the golden standard in the current practice of LST measurements and development of new methods. We will focus our attention on products based on two instruments on board of NASA's Terra satellite: Moderate Resolution Imaging Spectroradiometer (MODIS) and ASTER, [16]. ASTER has five TIR bands with relatively high spatial resolution of 90 m, but with significantly longer revisit time than the coarse resolution MODIS with 1 km spatial resolution and revisit time of 1~2 days. ASTER based LST products are based on a specially designed temperature and emissivity separation (TES) algorithm, [17], which first estimates emissivity, by using multiple observations from five bands and regression of laboratory emissivity values, then iteratively removes the influence of reflected sky irradiance, and finally estimates LST. Based on such measurements that were performed between 2000-2008, recently was published the ASTER global emissivity database (GED), [3], which is considered as the most detailed emissivity map of the Earth, [10], providing average spectral emissivity at 100 m for each of the ASTER's five TIR bands, averaged over designated time frame, as well as providing monthly emissivity estimates at 5 km resolution for the time period between 2000-2015. This global scale database was an extension of NAALSED, [4], which was published in 2008. However, since 2008 there is an anomaly with ASTER SWIR sensors, reducing its cloud detection capabilities. In addition, there were also some research efforts in developing "split-window" algorithms for pairs of ASTER TIR bands, [18]. On the other hand, MODIS has much richer spectral capabilities, but at much coarser resolution. In order to produce emissivity estimates at 5 km resolution MODIS exploits day/night LST estimation algorithm, [7], which reduces the uncertainty about surface spectral emissivity by introducing day and night radiance measurements, but at the cost of having two LST temperatures instead of one, requiring at least seven radiance measurements at each acquisition time. The same specification, [7], also defines a "split-window" LST algorithm, which is used in another, alternative LST product provided by MODIS at 1 km resolution. It should be mentioned that MODIS instrument is also present on Aqua satellite, and together with another one on Terra provides LST products at the global scale on a daily basis. Performance of these satellite LST estimates has been also a subject of ground based studies, like [19].
In order to illustrate a simple single channel LST estimation using satellite observations acquired by ASTER TIR_Band12, [3], an image of computed LST in a bare land mining region of eastern Serbia is presented in Fig. 4.

CONCLUSIONS AND FUTURE WORK
Although LST estimation is a well understood problem that has been solved by various techniques in the past, it still represents a challenging task that requires detailed planning and validation of adopted measuring procedures. This is especially true in the current setting that is characterized by the increased availability of satellite based observations that are utilizing different types of sensors and have different mission objectives. We have summarized some basic concepts related to LST estimation using satellite observations and pointed out some of the main estimation approaches in this broad research domain. Global scale products, like the mentioned ASTER GED surface emissivity database, or daily MODIS LST estimates, are particularly important for further application of LST estimation in other analyses and experimental studies that can benefit from such information and rely on higher accuracy LST estimates. However, new generation of sensors and their novel applications will require refinement of such global scale products and development of new estimation techniques that will adequately address some of the challenges that were discussed in referenced works and this paper.