Land Surface Temperature Retrieval from Passive Microwave Satellite Observations: State-of-the-Art and Future Directions

: Land surface temperature (LST) is an important variable in the physics of land–surface processes controlling the heat and water ﬂuxes over the interface between the Earth’s surface and the atmosphere. Space-borne remote sensing provides the only feasible way for acquiring high-precision LST at temporal and spatial domain over the entire globe. Passive microwave (PMW) satellite observations have the capability to penetrate through clouds and can provide data under both clear and cloud conditions. Nonetheless, compared with thermal infrared data, PMW data su ﬀ er from lower spatial resolution and LST retrieval accuracy. Various methods for estimating LST from PMW satellite observations were proposed in the past few decades. This paper provides an extensive overview of these methods. We ﬁrst present the theoretical basis for retrieving LST from PMW observations and then review the existing LST retrieval methods. These methods are mainly categorized into four types, i.e., empirical methods, semi-empirical methods, physically-based methods, and neural network methods. Advantages, limitations, and assumptions associated with each method are discussed. Prospects for future development to improve the performance of LST retrieval methods from PMW satellite observations are also recommended.


Introduction
Land surface temperature (LST) is an important parameter in the interaction between the Earth's surface and the atmosphere [1][2][3]. It is an essential variable of land surface energy budget and a crucial parameter for studying land surface physical processes at local and global scales [4][5][6][7][8]. LST is extensively applied in various studies, including hydrology, climatology, ecology, meteorology, agronomy, environmental monitoring, and public health [9][10][11][12][13][14][15][16]. However, due to the strong temporal and spatial heterogeneity of land surface parameters and land surface-atmospheric energy exchange, LST varies rapidly with time and space, making it difficult to characterize the temporal variation and spatial distribution of LST using traditional methods [17]. Owing to the progress of space-borne remote sensing technique, satellite data provides an exclusive feasible way for estimating global-scale LST with satisfying temporal and spatial measurements other than point values obtained from traditional methods.
The acquisition of LST using remote sensing technology is mainly on the basis of thermal infrared (TIR). Space-borne TIR measurement is directly related to LST by way of radiative transfer equation (RTE). LST retrieval from TIR measurements has drawn a lot of attention since 1970s, and a variety of algorithms were put forward, e.g., single channel algorithm [18], split-window algorithm [19,20], temperature and emissivity separation algorithm [21], and day/night algorithm [22]. TIR-derived LST has a decent retrieval accuracy (1-2 K) and spatial resolution (e.g., 1 km for Moderate Resolution Imaging Spectroradiometer (MODIS) data). However, TIR radiation is easily influenced by atmospheric water vapor, wind speed, cloud, and rain-fall. This disadvantage shall be considered since cloud-covered surface occupies 60% of the surface of the Earth [23].
In comparison with TIR remote sensing, passive microwave (PMW) remote sensing is capable to penetrate cloud, rain, snow, vegetation, and other arid surface features, which can cover the shortage of TIR data with its ability to acquire all-weather surface information. Consequently, there have been a large series of researches on LST retrieval algorithms from PMW satellite observations. PMW-based LST retrieval methods can be mainly categorized into four types, i.e., empirical methods, semi-empirical methods, physically-based methods, and neural network methods [24][25][26][27]. Thus, summaries and comparisons of different PMW-based LST retrieval methods are urgently required and indispensable for us to better understand the mechanism of microwave radiative transfer process.
This paper provides an overview of various categories of algorithms that were established to estimate LST from PMW satellite observations. This paper is organized with six sections. In Section 2, we provide the theoretical basis for retrieving LST from PMW satellite observations. In Section 3, we briefly discuss the underlying obstacles in LST retrieval from PMW satellite observations. In Section 4, we gave a brief recall of concrete methods for retrieving LST from PMW observations from empirical methods to physically-based methods. For each method, we shall detail the main theory and assumptions involved in the model development, and highlight its advantages, drawbacks, and potential. In Section 5, the validation of PMW-derived LST is given in two aspects, including (i) direct comparison and (ii) inter-comparison. Finally, future development and perspectives are provided in Section 6.

Theoretical Background
All objects with an absolute temperature above 0 K emit radiation in the form of electromagnetic waves. Natural surfaces emit thermal radiation mainly in the far infrared region; however, it extends all through the electromagnetic spectrum into microwave region.
Planck's law is a function to describe the spectrum of emitted radiation from a blackbody with wavelength as independent variable. In the microwave region, the function is usually written with frequency as independent variable and can be expressed as where B(f,T) is the spectral radiance (W m −2 Hz −1 sr −1 ) at frequency f (Hz) of a blackbody at temperature T (K), c is the speed of light in vacuum (2.9979 × 10 8 m s −1 ), h is the Planck constant (6.6262 × 10 −34 J s), and k is the Boltzmann constant (1.38 × 10 −23 J K −1 ). In the range of low frequencies (i.e., long wavelengths), the Rayleigh-Jeans approximation (i.e., hf /kT << 1) of Planck's law is valid. Equation (1) is approximated as Therefore, the Planck function is approximated by the Rayleigh-Jeans' law in the microwave region. The radiation measured by a satellite-borne microwave sensor can be expressed as linear equation of temperature instead of radiance. Because most natural objects do not behave as a blackbody, the idea of emissivity is introduced to relate the radiance of an object at temperature T to that of a blackbody at the same temperature. The spectral radiance of a non-blackbody can be calculated by the product of spectral emissivity multiplying Planck's law as shown in Equation (2).
PMW remote sensing measures thermal radiation from the land surface in the centimeter wave band and is largely determined by the physical temperature and the emissivity of a radiating body. The emitted radiation in the microwave region (λ = 1-1000 mm) is low as compared with infrared radiation (λ = 8-14 µm). Different from infrared emission, which represents thermal radiation from the very shallow top of the surface, PMW emission represents thermal radiation from the integral of a shallow upper layer of the surface.

Bare Soil Surfaces
For non-scattering plane-parallel atmosphere, brightness temperature observed from the top of the atmosphere (TOA) over bare soil surfaces consists of three radiation components: (1) surface-emitted radiation attenuated by the atmosphere, (2) surface-reflected and then atmosphere-attenuated atmospheric downwelling radiation and cosmic background radiation, and (3) atmospheric upwelling radiation. The PMW RTE over bare soil surfaces can be written as where T b , T s , and T c are the brightness temperature at the TOA, the LST, and the cosmic background brightness temperature (approximately 2.7 K), respectively; Γ a = exp(-τ a ) is the atmospheric transmittance; τ a is the optical depth of atmosphere; ε s is the soil surface emissivity; and T ad and T au are the downwelling and upwelling atmospheric brightness temperatures, respectively. T ad and T au can be approximated by Equation (4) for PMW sensor with an Earth incidence angle at 53 • or 55 • : where T ae is the atmospheric effective temperature, which is highly relying on the vertical distribution of atmospheric parameters, e.g., atmospheric temperature, relative humidity, and liquid water [28]. Atmospheric optical depth τ a along atmospheric slant path is written as [29]: where τ o is the oxygen optical depth at nadir, a w and q w are the extinction coefficient and at nadir optical depth of atmospheric water vapor, respectively, a l and q l are the extinction coefficient and at nadir optical depth of cloud liquid water, respectively, and θ is the viewing zenith angle, i.e., the Earth incident angle of a PMW sensor.

Vegetated Surfaces
For vegetated surfaces, vegetation is expressed as a single-scattering layer lying above the rough soil. The impact of vegetation on the microwave radiation observed above the canopy can be attributed in two ways. On the one hand, the vegetation scatters and absorbs the radiation emitted from the soil. On the other hand, the vegetation itself emits radiation. These two effects tend to have a counteract effect between each other [30]. As illustrated in Figure 1, brightness temperature at the TOA over vegetated surfaces is approximated as the sum of five radiation components: (1) soil-emitted radiation attenuated by the vegetation and the atmosphere; (2) vegetation upwelling radiation attenuated by the atmosphere; (3) vegetation downwelling radiation reflected by the surface and attenuated by the vegetation and the atmosphere; (4) atmospheric downwelling radiation and cosmic background radiation attenuated by the vegetation, reflected by the surface, and attenuated by the atmosphere; and (5) atmospheric upwelling radiation. The PMW RTE over vegetated surfaces is written as [31]: where T v is the vegetation surface temperature, ω v is the vegetation single scattering albedo, Γ v is the vegetation transmittance and equals to exp(-τ v ) where τ v is the vegetation optical depth, which is can be expressed as a function of vegetation water content using Equation (7) [32]: where b v is a parameter that depends on vegetation types, and w v is the vegetation water content.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 28 and attenuated by the atmosphere; and (5) atmospheric upwelling radiation. The PMW RTE over vegetated surfaces is written as [31]: where Tv is the vegetation surface temperature, ωv is the vegetation single scattering albedo, Γv is the vegetation transmittance and equals to exp(-τv) where τv is the vegetation optical depth, which is can be expressed as a function of vegetation water content using Equation (7) [32]: where bv is a parameter that depends on vegetation types, and wv is the vegetation water content.

PMW Satellite Sensors and Specifications
PMW satellite observations were available since 1978 after the launch of Nimbus-7. Combining the continuous PMW observations shows a great potential in the generation of global-scale long term LST record. Table 1 lists the performance characteristics of PMW satellite sensors. A brief overview of PMW satellite sensors is described in the following sections. (2) vegetation upwelling radiation attenuated by the atmosphere; (3) vegetation downwelling radiation reflected by the surface and attenuated by the vegetation and the atmosphere; (4) atmospheric downwelling radiation and cosmic background radiation attenuated by the vegetation, reflected by the surface, and attenuated by the atmosphere; and (5) atmospheric upwelling radiation.

PMW Satellite Sensors and Specifications
PMW satellite observations were available since 1978 after the launch of Nimbus-7. Combining the continuous PMW observations shows a great potential in the generation of global-scale long term LST record. Table 1 lists the performance characteristics of PMW satellite sensors. A brief overview of PMW satellite sensors is described in the following sections.  Instrument  SMMR  SSM/I  TMI  AMSR-E  AMSR2  MWRI   Satellite  Nimbus-7  DMSP F8, F10, F11,  F13, F14, F15  TRMM  Aqua  GCOM-W1  FY-3A, 3B, 3C, [33,34]. It provides dual-polarized microwave brightness temperature measurements at five channels every other day, i.e., 6.63, 10.69, 18.0, 21.0, and 37.0 GHz. The instrument scans the Earth with a scanning width of 780 km and an incidence angle of 50.2 • [35]. The SMMR daily EASE-Grid brightness temperature is available from the National Snow and Ice Data Center (NSIDC) website (https://nsidc.org/data/NSIDC-0071).

Special Sensor Microwave/Imager (SSM/I)
The Special Sensor Microwave/Imager (SSM/I) was equipped on a series of the Defense Meteorological Satellite Program (DMSP) platforms designated F8, F10, F11, F13, F14, and F15 [36,37]. The series of DMSP satellites are in a near circular, sun synchronous, polar orbit. The first satellite F8 was launched on 18 June 1987, while the last one F15 was launched on 12 December 1999. Equator crossing times vary among different satellites. The instrument measures dual-polarized microwave brightness temperatures at 19.35, 22.235, 37.0, and 85.5 GHz, and vertical-polarized brightness temperatures at 22.235 GHz [38]. The swath width of SSM/I is 1394 km and the incidence angle is 53.1 • . The daily EASE-Grid brightness temperatures observed by SSM/I is downloadable from the NSIDC website (https://nsidc.org/data/NSIDC-0032).

Tropical Rainfall Measuring Mission (TRMM)-Microwave Imager (TMI)
The Tropical Rainfall Measuring Mission (TRMM)-Microwave Imager (TMI) was carried on the TRMM satellite launched on 27 November 1997, jointly by the National Aeronautics and Space Administration (NASA) and the Japan Aerospace Exploration Agency (JAXA) [39,40] and ceased on April 8, 2015. The orbit of TRMM satellite is near-equatorial, other than sun-synchronous. The spatial coverage is between approximately 38 • north and south latitude. The TMI microwave radiometer operates at five frequencies: 10.65, 19.35, 37.0, and 85.5 GHz in horizontal and vertical polarizations and 21.3 GHz in vertical polarization [41]. The swath width of TMI is 759 km and the incidence angle is 52.8 • [42]. The TMI brightness temperature is approachable from the EARTHDATA website (https://disc2.gesdisc.eosdis.nasa.gov/data/).

Advanced Microwave Scanning Radiometer for Earth Observing System (AMSR-E)
The Advanced Microwave Scanning Radiometer for Earth Observing System (AMSR-E) was carried on the Aqua satellite operated by NASA [43]. The Aqua satellite was sent into space on 2 May 2002, and ceased operation on 4 December 2011. It runs in a near-circular sun-synchronous orbit with an orbit altitude of 705 km. The AMSR-E instrument contains twelve channels including 6.925, 10.65, 18.7, 23.8, 36.5, and 89.0 GHz in both horizontal and vertical polarizations [44]. It scans the Earth's surface with a swath width of 1445 km and an incidence angle of 55 • [45]. The AMSR-E daily EASE-Grid brightness temperatures can be downloaded from the NSIDC website (https://nsidc.org/data/NSIDC-0301).

Advanced Microwave Scanning Radiometer 2 (AMSR2)
The Advanced Microwave Scanning Radiometer 2 (AMSR2) was launched by JAXA on 18 May 2012 with the Global Change Observation Mission-Water "SHIZUKU" (GCOM-W1) satellite [46]. The orbit of GCOM-W1 satellite is near-circular sun-synchronous, with an altitude of 700 km. The AMSR2 instrument is the descendent of the AMSR-E onboard the Aqua satellite. It measures microwave brightness temperatures at 6.925, 7.3, 10.65, 18.7, 23.8, 36.5, and 89.0 GHz in both horizontal and vertical polarizations [47]. An extra C-band (7.3 GHz) channel was added to reduce the impact of the Radio Frequency Interference (RFI). The swath width of AMSR2 scans is 1450 km and the incidence angle of AMSR2 is 55 • [48]. The AMSR2 brightness temperatures can be approached from the G-Portal address through ftp protocol (https://gportal.jaxa.jp/gpr/).

Microwave Radiation Imager (MWRI)
The Microwave Radiation Imager (MWRI) is a microwave radiometer onboard the FengYun-3 (FY-3) satellite series developed by the China National Space Administration (CNSA) [49]. The first satellite FY-3A was launched on 27 May 2008. The FY-3 series satellites are in a near-circular sun-synchronous orbit at 836 km above the sea level. The channel settings of MWRI instrument include ten channels, i.e., 10.65, 18.7, 23.8, 36.5, and 89.0 GHz in both horizontal and vertical polarizations [50]. The swath width of MWRI of 1400 km and incidence angle is 53 • [51]. The MWRI brightness temperatures are published by the Chinese National Satellite Meteorological Center (CNSMC) and are valid on its website for registered users (http://satellite.nsmc.org.cn/portalsite/default.aspx).

Difficulties and Problems for LST Retrieval from PMW Satellite Observations
LST retrieval from PMW satellite observations is still challenging. There are several difficulties and obstacles regarding to the LST retrieval as follows: (1) LST retrieval from PMW observations is mathematically underdetermined and unsolvable. Based on the PMW RTE, if brightness temperature at the TOA measured in N channels (i.e., N equations), there will be N + 1 unknown (one LST and N emissivities), even if other quantities are assumed to be known. (2) It is hard to achieve the decoupling of LST, surface emissivities, and downwelling atmospheric brightness temperature from the measured brightness temperature. The downwelling atmospheric brightness temperature and LST are coupled together through emissivities. On the one hand, the non-unity emissivity reduces the radiation emitted from a natural surface. On the other hand, surface emissivity increases the reflection of downwelling atmospheric brightness temperature back to the atmosphere. These two effects compensate each other, making the decoupling of LST and downwelling atmospheric brightness temperature more difficult. (3) Atmospheric correction is difficult to implement [52]. Although PMW observations are much more immune to clouds and water vapor than TIR observations, the atmospheric impact on radiance obtained by a PMW radiometer at the TOA cannot be fully neglected, especially under cloudy conditions. The atmospheric impact is mainly caused by the absorption and emission of three atmospheric components: water vapor, liquid water existing in clouds and rain, and oxygen molecule. Thus, an accurate prior knowledge of the vertical distribution of highly-variable atmospheric parameters, i.e., water vapor, liquid water, and temperature are strongly in demand. (4) Compared with surface emissivities in the TIR region, emissivities in microwave region have strong variation correlating to surface properties, including soil moisture, surface roughness, and vegetation cover [53,54]. Based on the PMW RTE, the uncertainty of emissivity leads to large error of retrieved LST when the variation of emissivity is large. Since the microwave emissivity occurs at vertical polarization is relatively smaller than that at horizontal polarization, vertical-polarized channels are of higher priority in retrieving LST. Furthermore, it is difficult to quantitatively estimate the perturbation of emissivity caused by the effects of soil moisture, vegetation cover, and surface roughness. (5) The LST retrieved from PMW satellite observations represents subsurface temperature whereas LST retrieved from TIR satellite observations normally represents skin temperature [23]. The depth of the retrieved LST from PMW satellite observations changes with wavelength and other surface parameters (e.g., soil moisture and soil texture). (6) Validation of LST derived from PMW observations at the satellite pixel level is also a challenging problem. Validation is ambiguous due to the strong temporal variation and spatial heterogeneity of LST. It is difficult to conduct in situ LST measurements and obtain LST data representing the PMW pixel level. The direct comparison between in situ temperature measurements and PMW-retrieved LST observations is also difficult to realize.

Methods for LST Retrieval from PMW Satellite Observations
Over the past several decades, a series of methods were put forward to retrieve LST from PMW satellite observations based on different assumptions and approximations of atmospheric parameters and surface emissivity. These methods can be roughly sorted as four categories: empirical methods, semi-empirical methods, physically-based methods, and neural network methods. Although all of these methods can be used to retrieve LST from PMW satellite observations, these methods have their unique advantages and drawbacks, and their application range alters significantly. In most methods, LST is explicitly expressed as a function of brightness temperature at the TOA, surface emissivity, and/or atmospheric parameters. In addition to these parameters, one important procedure is to determine regression coefficients in the LST retrieval methods. Three types of temperature dataset are commonly used to determine regression coefficients: (1) ground-based air temperature/LST, (2) data simulated by radiative transfer model, and (3) satellite-derived TIR LST.

Empirical Methods
The empirical methods use the site-specific relationship between ground-based LST measurements and PMW satellite observations at single or multiple channels to retrieve LST. The advantage of the empirical methods is the simplicity and convenience. The main weakness is that empirical methods need determination of site-specific parameters. This weakness has more or less limited the applications of these methods at global scale under various surface and atmospheric conditions.

Single-Channel Regression Method
The single-channel regression method uses the empirical relationship between LST and passive microwave brightness temperature at a single channel to retrieve LST. The Ka band (37 GHz) is thought as the most practicable channel in microwave region for LST retrieval because Ka band is less sensitive to land surface characteristics and has a relatively higher atmospheric transmittance than other bands [24]. Compared with the horizontally polarized channel, the vertically polarized channel is more suitable for LST retrieval because vertical-polarized emissivity is relatively less affected by soil moisture variations at the earth incidence angle of 50 • -55 • . In previous studies, several researchers found that the relationship between LST and the vertically polarized brightness temperature at 37 GHz can be regressed a linear function [24,[55][56][57][58][59][60]. Thus, the single-channel regression method for LST retrieval can be described as where T b,37V is the brightness temperature at 37 GHz vertically polarized channel, and a s and b s are the regression coefficients.
To improve the performance of single-channel regression method, some studies [55,60] established the regression relationships between LST and the brightness temperature at 37 GHz vertically polarized channel for daytime/nighttime (descending/ascending) data, respectively. These studies indicate that both of the accuracies of the separated regression equations for daytime and nighttime data are relatively higher than those of the regression equation for all data.

Multi-Channel Regression Method
The multi-channel regression method builds an empirical equation between LST and PMW brightness temperatures at multiple channels to estimate LST [61]. In addition to the PMW brightness temperature, the combinations of brightness temperatures at different channels and/or index variables derived from these channels (e.g., microwave polarization difference index) are also selected as potential predictors in the regression relationship [62][63][64][65]. The multi-channel regression method is written as 10.7, 19, 22, 37, 89 GHz (9) where T b,xiH and T b,xiV are the brightness temperature at x i (x i = 6.9, 10.7, 19, 22, 37, and/or 89) GHz horizontally and vertically polarized channels, and a m , b m , c m , d m , and e m are the regression coefficients.
To further upgrade the LST retrieval result of the multi-channel regression method, some studies [66,67] built the regression relationship for each land cover type. Furthermore, the stepwise regression method was used to remove redundant regression variables [67,68].
In this method, the combinations of brightness temperatures and index variables are always used to indicate the impact of atmosphere and vegetation on LST retrieval. Thus, compared to the single-channel regression method, this method considers some influence factors on LST retrieval process and naturally improves the accuracy of retrieved LST. However, the multi-channel regression method is still an empirical-based method, which is lack of physical meaning. It cannot be applied at global scale.
The advantage of the multi-channel regression method is that usage of brightness temperature from multiple channels avoids more complicated physical or mathematical derivations. The disadvantage of multi-channel method is that it cannot explain how surface conditions function in the LST retrieval procedure and the improvement of accuracy basically relies on repeated trial and error.

Semi-Empirical Methods
The semi-empirical methods introduce some empirical constraints of emissivity to make the LST solvable from radiative transfer model. This can partially restrain the errors of global scale algorithms when used in local regions because its calculation is on the basis of the radiative transfer model.
The disadvantage of semi-empirical methods is that the empirical constraint of surface emissivity varies with the land surface conditions, which limits its application. In addition, semi-empirical methods may improperly introduce or ignore some environmental variables from radiative transfer model in the retrieval process, which can bring errors in the retrieved LST.

Single-Channel Emissivity Method
Single-channel emissivity method estimates surface emissivity at single channel using auxiliary data, then uses the emissivity to retrieve LST from radiative transfer model. Gao et al. [25] estimated surface emissivity at 19 GHz using the polarization ratio (PR) of brightness temperature. The relationships between emissivity and PR for forested and non-forested areas are as follows: where ε f,19H and ε nf,19H are the horizontal emissivities over the forested and non-forested areas, respectively, a 1 , b 1 , c 1 , a 2 , b 2 , and c 2 are the coefficients to be regressed, and PR 6 and PR 19 are the PR values at 6.9 and 19 GHz, respectively. The PR is calculated by horizontal brightness temperatures dividing vertical brightness temperatures (PR = T bH /T bV ).
Taking into account that low-frequency microwave observations are barely interfered by atmospheric conditions, the coefficients in Equation (10) were calculated from collocated microwave brightness temperature and TIR LST during a dry season. Gao et al. [25] assumed that the regression relationships established during a dry season can be applied to all seasons.
Once the emissivity is available, LST can be retrieved from microwave brightness temperature: This method is a practical method to retrieve LST under all-weather conditions over tropical forest areas [69]. Except for land cover type image, no other auxiliary data are required. However, the assumption that the regression relationships established during a dry season (under clear-sky conditions) is applicable for all seasons (under all-weather conditions) should be further verified.

Dual-Polarization Emissivity Method
Dual-polarization emissivity method builds the relationship between horizontally polarized and vertically polarized emissivity, which can make the LST solvable from radiative transfer model. A simple method for LST retrieval from PMW brightness temperatures over snow-and ice-free surfaces was developed by Fily et al. [70]. This method is based on an empirical linear relationship between vertical and horizontal polarized emissivities at 37 GHz: where ε 37V and ε 37H are the vertical and horizontal emissivities at 37 GHz, and m and n are the regression coefficients. Assuming that atmospheric parameters are known, LST can be deduced from the PMW RTE for two polarizations at 37 GHz by means of the empirical linear regression. The polarized emissivity relationship method is written as To facilitate the applications of this method to retrieve LST from long-term PMW brightness temperatures in high-latitude areas when there is no available prior knowledge on atmospheric status, a simple atmospheric correction for PMW satellite observations at 37 GHz was performed by Fily et al. [70] on the assumption that the atmosphere can be described as no clouds and a mean constant integrated water vapor content of 1.5 g cm −2 . Atmospheric parameters were then given as a constant [71,72].
This method is less insensitive to cloud in the atmosphere. However, for global application, better performance would be realized when using actual atmospheric parameters derived from ancillary data in the dual-polarization emissivity method. The limitation of this method is that it is easily affected by the uncertainties of the regression coefficients m and n, especially for dense vegetation surface when the coefficients is close to 0, which cause large error for the retrieved LST. This method malfunctions when applied at global scale because the assumption of constant atmospheric conditions and empirical coefficients got at local scale are not universal globally. Therefore, André et al. [73] estimated pixel-based coefficients m and n using least square regression during a certain time period assuming that the polarized emissivity relationship is pixel dependent (i.e., variable with space) and temporally invariable (i.e., m and n are assumed as constant during the given period).

Multi-Channel Emissivity Method
Multi-channel emissivity method builds the relationship of emissivity among multiple channels; thus, it makes the LST solvable from radiative transfer model. A three-channel emissivity method was proposed by Basist et al. [74]. This study assumes that surface emissivity can be expressed as the difference between the nominal emissivity over dry surfaces and an emissivity adjustment item: where ε 19V is the surface emissivity at 19 GHz channel, ε 0,19V is the nominal emissivity over dry surfaces at 19 GHz channel, and ∆ε 19V is the change of emissivity at 19 GHz channel due to surface wetness. The studies of Basist et al. [74] showed that surface emissivity changes with frequency and surface wetness ( Figure 2). As surface wetness decreases, surface emissivity increases and the gradient of emissivity between low and high frequencies decreases. Because the variation of emissivity along with frequency over wet surfaces is nonlinear, the gradient of emissivity between two or more frequencies is used to describe the change of emissivity due to surface wetness: where ε 37V and ε 89V are the surface emissivities at 37 and 89 GHz channels, respectively, and α 1 and β 1 are the regression coefficients.
where ε37V and ε89V are the surface emissivities at 37 and 89 GHz channels, respectively, and α1 and β1 are the regression coefficients.
. Neglecting atmospheric effects on microwave radiation, the microwave brightness temperature estimation can be written as Combining Equations (14)- (16), the vertically polarized brightness temperature at 19 GHz channel can be expressed as LST can then be deduced from Equation (17): Neglecting atmospheric effects on microwave radiation, the microwave brightness temperature estimation can be written as Combining Equations (14)- (16), the vertically polarized brightness temperature at 19 GHz channel can be expressed as LST can then be deduced from Equation (17): where T b,37V = ε 37V T s and T b,89V = ε 89V T s are the vertically polarized brightness temperature at 37 and 89 GHz channels, respectively. In this method, brightness temperature at 19 GHz vertically polarized channel is used as the primary channel to retrieve LST, and the combinations of brightness temperatures at various vertically polarized channels are used to further compensate the impact of surface emissivity on LST retrieval due to the variation in surface wetness.
The performance of this method was further analyzed and evaluated by Williams et al. [75] and Basist et al. [76]. This method has three assumptions: (1) dry surfaces have a nominal emissivity of 0.95 at 19 GHz, (2) atmospheric effects can be neglected, and (3) the change of emissivity due to surface wetness can be described as a function of the slope of emissivity between two or more frequencies.
The advantage of this method is its applicability to satellite PMW observations over snow-and ice-free surfaces at global scale without other auxiliary data once the regression coefficients α 1 and β 1 are available. The limitation of this method is that the emissivity of dry surfaces in fact is not a constant of 0.95, and it changes with soil type and texture.
To improve the performance of the three-channel emissivity method, Han et al. [77] expressed the emissivity adjustment term as the combination of the emissivity at different channels: where ε 0,89V is the nominal emissivity over dry surfaces at 89 GHz.
According to the simulation analysis, the first item in the right-hand side of Equation (19) can be quantified as a quadratic function of the difference between emissivities at 37 and 89 GHz (i.e., ε 89V -ε 37V ). Moreover, the last item in the right-hand side of Equation (19) can be quantified as a constant. Consequently, surface emissivity at 19 GHz is expressed as where α 2 , β 2 , and γ 2 are the regression coefficients. Similarly, assuming that atmospheric effects are negligible and combining Equations (16) and (20), LST can be described as As mentioned previously, because the emissivity of dry surfaces is related to soil texture and soil type, it is hard to accurately determine the emissivity of dry surfaces. In addition, the variation in the denominator item in the right-hand side of Equation (21) (i.e., ε 0,19V -γ 2 ) is small. Therefore, replacing it using its average value is practicable. To facilitate the application of this method, Equation (21) can be rewritten as [77]: where A-E are the regression coefficients. LST is retrieved based on the initial value of T s estimated from brightness temperature at 37 GHz vertically polarized channels, and the Newton iterative method is used to solve Equation (22).

Physically-Based Methods
The physically-based methods have strong theoretical basis and clear physical meaning. However, the complexity of the RTE determines that it is difficult to be applied in extensive application. Simplification and assumption of the RTE in actual retrieval restrain the applicability of the physically-based methods.

Dual-Channel Physical Method
Dual-channel physical method introduces the relationship of emissivity/transmittance between adjacent channels into radiative transfer model, which makes LST solvable. On the basis of the PMW RTE, the TOA brightness temperature can be written as follows by substituting Equation (4) into Equation (3): where ∆T = T s -T ae is the difference between LST and atmospheric effective temperature. According to the simulation analysis, Weng and Grody [26] pointed out that the second term in the right side of Equation (23) is insensitive to emissivity and can be parameterized as a linear function of atmospheric water vapor content (WVC) for high emissivity surfaces (ε > 0.9). Equation (23) can be rewritten as where γ is the frequency-dependent regression coefficient.
Combining the brightness temperatures at 19 and 22 GHz vertically polarized channels, the following relationship can be deduced from Equation (24): where T b,22V is the brightness temperature at 22 GHz vertically polarized channel, ε 22V is the surface emissivity at 22 GHz, τ 19V and τ 22V are the atmospheric transmittance at 19 and 22 GHz, respectively, γ 19V and γ 22V are the regression coefficients at 19 and 22 GHz, respectively, and W is the atmospheric WVC. Weng and Grody [26] assumed that (1) surface emissivities at 19 and 22 GHz are nearly identical, and (2) the transmittance of clouds (as well as oxygen) at 19 and 22 GHz is nearly identical. The difference between brightness temperature measurements at the two adjacent frequencies is mainly caused by the different absorption effect of atmospheric water vapor. Therefore, the ratio of atmospheric transmittance at 19 and 22 GHz in the right-hand side of Equation (25) can be quantified as a function of WVC in the form of exponent, i.e., τ 19V /τ 22V = e cW , where c is the coefficient to be regressed. On the basis of the two assumptions, Weng and Grody [26] proposed a dual-channel physical method for retrieving LST from brightness temperature measurements at 19 and 22 GHz. This method can be written as Taking into account that atmospheric WVC is difficult to retrieve from PMW satellite observations over land, Weng and Grody [26] analyzed the WVC and T s in field measurements and observed an empirical relationship between them to estimate atmospheric WVC in Equation (26). Based on the initial value (or the first guess) of T s estimated from the brightness temperature at 89 GHz vertically polarized channel, the Newton iterative method was used to solve Equation (26) for retrieving LST.
The dual-channel physical method is developed on the basis of the PMW RTE and has a clear physical meaning. This physically-based method can be applied to retrieve LST from PMW satellite observations over snow-and ice-free surfaces at global scale. The limitation of the method is that large errors were found under dry atmosphere conditions with low WVC.
According to the simulation analysis, Huang et al. [78] found that the parameterization scheme is still valid for surfaces with medium emissivity (ε = 0.8-0.9). Furthermore, surface emissivity at 19 GHz does not equal with it at 22 GHz but an obvious linear relationship between them can be observed. Therefore, one additional regression coefficient a between 1-ε 19V and 1-ε 22V in the right side of Equation (25) was used to improve the performance of the dual-channel physical method. The improved method can be written as where a is the regression coefficient between 1-ε 19V and 1-ε 22V .
To make the improved dual-channel physical method more robust, Huang et al. [78] estimated WVC from atmospheric reanalysis data (i.e., NCEP data) rather than the empirical relationship between WVC and T s used in Weng and Grody [26].

Multi-Temporal Physical Method
Multi-temporal physical method introduces the temporal information to the RTM, which makes the LST solvable from RTM. This method assumes that the surface emissivity is invariant over a time period. Pulliainen et al. [62] proposed a multi-temporal physical method to retrieve LST from PMW satellite observations based on the hypothesis that surface emissivity over boreal forest areas in specific location is close to stable (i.e., invariable in one period of time) over snow-free surfaces. This method is implemented in two stages.
In the first stage, the maximum likelihood algorithm was used to estimate LST for each measurement and mean surface emissivity from multi-temporal measurements. It is assumed that the impact of atmosphere can be characterized by the mean atmospheric conditions. Furthermore, the duration of multi-temporal measurements should last sufficiently long to ensure that the atmospheric conditions can be expressed as an average value. The inversion problem can be expressed as a constrained least-square minimization between the simulated and observed brightness temperatures for multiple channels and multi-temporal measurements: with respect to ε i and T s,j (28) where N is the number of microwave channels, M is the number of multi-temporal measurements, T obs b,ij is the observed brightness temperature for ith channel and jth measurement, T sim b,ij is the simulated brightness temperature for ith channel and jth measurement, which is a function of surface emissivity (ε i ) at ith channel, T s,j at jth measurement, and a constant (χ) that characterizes the mean atmospheric conditions.
Once the mean emissivity for a given location were determined in the first stage, these values were further used to retrieve LST for any measurements in the second stage by assuming that the estimated emissivity are invariant with time over the time period of the multi-temporal measurements. In this case, the atmospheric effects are not necessary to be characterized by the mean atmospheric conditions. The atmospheric parameter χ, with a given average value and standard deviation, can be expressed using a Gaussian random variable or a free parameter, rather than a constant value as in the first stage [79]. The maximum "a posterior" estimation method was used to solve the constrained least-squares minimization between the simulated and observed brightness temperatures for multiple channels and one measurement: where σ i is the standard deviation of brightness temperature at channel i, λ is the standard deviation of χ, andχ is the mean value of χ. Furthermore, Xiang and Smith [80] assumed that surface emissivity is invariant in two observations during a day, e.g., satellite descending and ascending overpass, to simultaneously retrieve LST and emissivity from PMW brightness temperature measurements. The inversion problem can then be expressed as a non-linear least-squares minimization between the simulated and observed brightness temperatures for two times at multiple channels: with the convergence test defined as where R k is the kth residual between the simulated and observed brightness temperatures, k is the iteration step, N is the number of microwave channels, T obs b,ij is the observed brightness temperature at ith channel for jth time, T sim b,ij is the simulated brightness temperature at ith channel for jth time, which is a function of LST at jth time (T s,j ) and emissivities at ith channel (ε i ).
The Levenberg-Marquardt algorithm was used to perform the least-squares minimization. The constraints in the retrieval method were set as 280-340 K for LST and 0.7-1.0 for surface emissivity, respectively. The convergence is typically achieved no more than two iterations for most pixels [80]. However, the convergence process is very sensitive to the value used for the convergence test.
The gridded temperature-moisture profiles obtained from the reanalysis data, e.g., European Centre for Medium-Range Weather Forecasts (ECMWF) and National Centers for Environmental Prediction (NCEP), were used to account for atmospheric effects in the radiative transfer calculations. If significant discrepancies between the gridded temperature-moisture profiles and actual atmospheric conditions exist, this method may not find solutions for LST and emissivity. Therefore, it is necessary to remove cloud-contaminated pixels with significant amounts of cloud liquid water prior to the application of the multi-temporal physical method [80].
On the assumption that surface emissivity for a given location is invariant with time, the multi-temporal physical method proposed a novel idea to retrieve LST and emissivity from multi-temporal observations. Since the snow cover considerably changes surface emissivity, this method can only be applied under snow-free conditions.

Multi-Channel Physical Method
Multi-channel physical method simultaneously retrieves M geophysical parameters (e.g., LST and soil moisture) from RTE using PMW satellite observations at N channels (e.g., 10.65 GHz, 18.7 GHz, and 36.5 GHz), where M should be less than N in a robust retrieval. An iterative least-squares minimization method is always used to solve this problem. The minimization method is based on the weighted-sum of squared differences (χ 2 ) between the simulated and observed brightness temperatures: (32) where N is the number of microwave channels, T obs b,i is the observed brightness temperature at channel i, T sim b,i is the simulated brightness temperature at channel i, x is the set of geophysical parameters (e.g., LST and soil moisture), and σ i is the standard deviation of measurement noise at channel i.
The Levenberg-Marquardt algorithm was employed to calculate the set of geophysical parameters x that makes the value of χ 2 minimal [68]. At each iteration, the algorithm starts with prior values x 0 of the geophysical parameters and regulates these values until the functions converge to the minimum of χ 2 .
Taking into account that the polarization difference brightness temperature (∆T b = T bV − T bH ) exhibits higher sensitivity to various land surface characteristics, Guha and Lakshmi [81] used the weighted-sum of squared differences (χ 2 ) between ∆T obs b,i and ∆T sim b,i to retrieve the geophysical parameters from PMW observations. The advantage of this method is that it can retrieve several land surface parameters simultaneously. The limitation is that some atmospheric and vegetation variables may be mistakenly assumed or neglected in the retrieval process to make the unknowns solvable, which leads to large error for the retrieval.

Neural Network Methods
The neural network methods are computationally efficient for the estimation of geophysical parameters [82]. It has satisfactory performance in solving nonlinear problems by exploring the inherent statistical correlations among the retrieved parameters.
As one of the most widely used neural networks, multi-layer perceptron (MLP) consists of the input layer, the output layer, and the hidden layer. Each layer includes one or more neurons directionally linked with the neurons from the previous and the next layer. A key step of constructing MLP is the determination of the topological structure, i.e., the dimension of the input vectors, the number of hidden layers, the number of neurons in each hidden layer, and the dimension of the output vectors. The number of inputs and outputs in the MLP can be fixed entirely by the problem itself. However, determining the number of hidden layers and the number of neurons in each hidden layer appears to be more ambiguous and is difficult to determine through theoretical analysis. Consequently, these parameters must be treated differently from one case to another depending on the complexity of the problem at hand. Figure 3 shows an example of a 3-layer perceptron with three inputs, two outputs, and the hidden layer containing five neurons. A neural network inversion method with first guess information was proposed to simultaneously retrieve LST, atmospheric water vapor, cloud liquid water, and emissivities over snow-and ice-free surfaces from PMW satellite observations [27]. To constrain the inversion process, this method optimizes the use of first guess information. The neural transfer function is written as where y  is the retrieved geophysical parameters (e.g., LST and emissivities), gW is the neural network g with parameters W, y b is the first guess for the geophysical parameters, x o = RTM(y) + η are the simulated observations using a radiative transfer model with y as the initial input, and η is the noise in observation. Further analysis of the theoretical errors in the neural network inversion method can be found in Aires [83] and Aires et al. [84,85]. Prigent et al. [86] evaluated the LST retrieved by neural network inversion method relative to surface 2-m air temperature by means of a combined analysis of the differences between LST and surface air temperature as a function of solar flux, soil characteristics, and cloudiness. Catherinot et al. [87] further evaluated the LST retrieved by neural network inversion method by comparing retrieval results with in situ measurements in various environmental conditions within one-year period. The neural network inversion method was also used to estimate LST over snow-covered surfaces from combined PMW and TIR measurements [88].
The advantages of the neural network methods are that they are relatively easy to implement and can approximate any function regardless of its linearity. Moreover, the neural network methods are attractive for complex problems. The disadvantages of these methods are often abused in cases where simpler solutions such as linear regression would be best. The neural network methods are black boxes, which are incapable to specify the theoretical grounding and reasoning process. Furthermore, the neural network methods require a large number of cases to train the neural network. The methods cannot function well when there are insufficient training data.
In addition to the traditional neural network methods, other machine learning methods, e.g., random forest and convolutional neural network, have attracted more and more attention in recent years. These more advanced methods provide an alternative and feasible way to retrieve LST from

Input layer
Hidden layer output layer Figure 3. Topological structure of a 3-layer perceptron with three inputs, two outputs, and the hidden layer containing five neurons.
A neural network inversion method with first guess information was proposed to simultaneously retrieve LST, atmospheric water vapor, cloud liquid water, and emissivities over snow-and ice-free surfaces from PMW satellite observations [27]. To constrain the inversion process, this method optimizes the use of first guess information. The neural transfer function is written aŝ whereŷ is the retrieved geophysical parameters (e.g., LST and emissivities), g W is the neural network g with parameters W, y b is the first guess for the geophysical parameters, x o = RTM(y) + η are the simulated observations using a radiative transfer model with y as the initial input, and η is the noise in observation. Further analysis of the theoretical errors in the neural network inversion method can be found in Aires [83] and Aires et al. [84,85]. Prigent et al. [86] evaluated the LST retrieved by neural network inversion method relative to surface 2-m air temperature by means of a combined analysis of the differences between LST and surface air temperature as a function of solar flux, soil characteristics, and cloudiness. Catherinot et al. [87] further evaluated the LST retrieved by neural network inversion method by comparing retrieval results with in situ measurements in various environmental conditions within one-year period. The neural network inversion method was also used to estimate LST over snow-covered surfaces from combined PMW and TIR measurements [88].
The advantages of the neural network methods are that they are relatively easy to implement and can approximate any function regardless of its linearity. Moreover, the neural network methods are attractive for complex problems. The disadvantages of these methods are often abused in cases where simpler solutions such as linear regression would be best. The neural network methods are black boxes, which are incapable to specify the theoretical grounding and reasoning process. Furthermore, the neural network methods require a large number of cases to train the neural network. The methods cannot function well when there are insufficient training data.
In addition to the traditional neural network methods, other machine learning methods, e.g., random forest and convolutional neural network, have attracted more and more attention in recent years. These more advanced methods provide an alternative and feasible way to retrieve LST from PMW satellite observations. Mao et al. [89] used a deep dynamic learning neural network method to retrieve LST from AMSR2 observations. Tan et al. [90] constructed a deep learning convolutional neural network to retrieve LST from AMSR2 data by the combination of twelve vertical and horizontal polarization channels of 7.3, 10. 65, 18.7, 23.8, 36.5, and 89 GHz. Yoo et al. [91] estimated all-weather LST from AMSR2 brightness temperature and other auxiliary parameters using a random forest model.

Validation of PMW-Derived LST
Validation is a crucial procedure to evaluate the accuracy of PMW-derived LST products. It provides a reliable message in regard to the quality of LST products to potential LST users, as well as feedback information to the developers of LST retrieval methods for further improvement. As previously mentioned, validation of PMW-derived LST is still challenging due to the strong temporal variation and spatial heterogeneity of LST. Two methods were widely used for validation of PMW-derived LST, i.e., direct comparison method and inter-comparison method. The direct comparison method directly compares PMW-derived LST with ground-based temperature measurements at the satellite overpass time [82,87,92,93]. The inter-comparison method is a cross-validation procedure in which a well validated LST product was used as reference LST and compared with the PMW-derived LST [70,82,94,95]. Each method has its advantages and disadvantages. Note that both these two methods cannot be considered as "actual" validation of PMW-derived LST due to the difficulty in obtaining actual values at the microwave satellite pixel scale. We prefer to refer to this validation as comparison.

Direct Comparison Method
The direct comparison method is widely used for the comparison of PMW-derived LST and in situ LST measurements [24,58,67,82,87,92,96]. Unlike the inter-comparison method, the direct comparison method can be performed in both clear-and cloudy-sky cases. The disadvantage of this method is the scale mismatch issue, i.e., PMW-derived LST is spatially integrated measurements at pixel scale, whereas in situ LST is point measurements. In theory, the direct comparison of PMW-derived LST versus in situ measurements should be performed over large and homogeneous surfaces. However, it is difficult to select large and homogeneous surfaces at the pixel scale of approximately 25 km × 25 km all over the world.
The PMW-derived LST and in situ LST measurements should be collocated in space and time before comparison. For a meaningful comparison of PMW-derived LST and in situ LST measurements, sites located in heterogeneous surfaces should be discarded, especially the sites located in coastal regions because emissivity of land surfaces is much higher than that of water surfaces. Similarly, the iceand snow-covered pixels should be removed. Taking into account the large spatial heterogeneity during daytime, the comparison of PMW-derived LST and in situ LST measurements are often restricted to nighttime to take advantage of surface thermal stability. In addition to in situ LST, previous studies also used surface air temperature measured at meteorological stations as reference data to compare PMW-derived LST because surface air temperature is strongly related to LST [25,26,70,71]. Table 2 summarizes the accuracies of different methods for retrieving LST from PMW satellite observations.

Inter-Comparison Method
TIR-derived LST is widely used as reference data to evaluate PMW-derive LST [70,82,94]. There are three advantages for this method. First, the comparison of TIR-and PMW-derived LST is conducted at the satellite pixel scale instead of the comparison of in situ point measurements and PMW-derived LST with a large footprint. Second, this comparison can be performed over various land cover types at global scale [94]. Third, the quality control in the TIR-derived LST product can be used to remove water bodies and cloudy pixels. This is because water bodies with very low surface emissivity have significant influences on the PMW LST retrieval. However, the comparison of TIR-and PMW-derived LST can only be performed in clear-sky cases because TIR remote sensing cannot measure LST in overcast areas.
The comparison of TIR-and PMW-derived LST is very suitable for cases where both TIR and PMW sensors are on-board the same satellite [94], e.g., both MODIS and AMSR-E sensors are on-board the Aqua satellite and both Visible and Infrared Radiometer (VIRR) and MWRI sensors are on-board the Chinese FY-3 satellite. Some studies also conduct this comparison for TIR and PMW sensors on board different satellites by taking into account the acquisition time difference, e.g., the comparison of AATSR-and SSM/I-derived LST with a maximum time difference of 1 h [82] and of Advanced Very High Resolution Radiometer (AVHRR) and SSM/I-derived LST with a maximum time difference of 30 min [70]. This time difference could lead to large discrepancies between TIR-and PMW-derived LST.
To perform the inter-comparison, TIR-derived LST with higher spatial resolution should be aggregated to match the pixel scale of PMW-derived LST. The fraction of TIR clear-sky pixels within a PMW pixel is determined by means of TIR cloud masking product. Taking into account the balance between comparison accuracy and the number of available pixels, different fraction thresholds of TIR clear-sky pixels within a PMW pixel are used to determine whether the PMW pixel is considered as a clear-sky pixel, e.g., 10% [73], 50% [82], and 98% [97]. To minimize the effects of clouds, Ermida et al. [94] considers a PMW pixel as a clear-sky pixel only when all TIR clear-sky pixels within the PMW pixel have best quality. According to the fraction threshold, TIR-derived clear-sky LST within a PMW pixel is averaged to compare with PMW-derived LST [58]. In addition to TIR-derived LST products, the LST products from the International Satellite Cloud Climatology Project (ISCCP [97]) and the gridded reanalysis data (e.g., ECMWF and NCEP) are also used as reference data to inter-compare PMW-derived LST [71,73,81,82,87].

Physical Meaning of PMW-Derived LST
The physical meaning of LST should be clearly defined with absolute certainty. However, the physical significance of PMW-derived LST is not fully clarified, especially when the surface is heterogeneous and non-isothermal. Different from skin temperature retrieved from TIR measurements, PMW-derived LST can be considered as radiometric temperature from a certain depth, normally ranging from a few millimeters to several centimeters underneath the surface owe to the penetration effects of microwave radiation. However, effective radiative depth of microwave radiation is difficult to determine due to the lack of thorough understanding of microwave emission. Radiative depth of different microwave channels varies with each other, making the definition of PMW-derived LST even more complicated. Over dense forest covered surfaces, PMW-derived LST refers to forest canopy temperature, which is different from the temperature over bare soils. Therefore, more efforts should be made to further refinement of physical meaning of PMW-derived LST and the conversion between different LST because directly using PMW-derived LST in climatological or ecological models without conversion can cause unexpected uncertainties in relevant applications.

Establishment of Microwave Spectral Emissivity Library
Surface emissivity is critical for LST retrieval from PMW observations. To date, several models for microwave emissivity simulation were developed, e.g., the community radiative transfer model, (CRTM) [98] and advanced integral equation model (AIEM) [99]. These models can simulate spectral emissivity at different polarizations and frequencies with preset surface roughness parameters and volumetric soil moisture. Moreover, satellite-derived microwave emissivity products are also publicly available, e.g., the AMSR-E/Aqua Monthly Global Microwave Land Surface Emissivity product [100]. Both the simulated emissivity and the satellite-derived emissivity are very useful to characterize the spatio-temporal variations of microwave emissivity [101][102][103][104][105]. However, these emissivity products with large modeling or retrieval error cannot represent the spectral variations of actual microwave emissivities. In fact, compared with surface emissivities in TIR region, microwave emissivities are highly variable when surface properties such as, vegetation cover, soil moisture, and surface roughness vary. To develop new physically-based methods for retrieving LST from PMW observations, it is necessary to have a better understanding of the spectral variations of actual microwave emissivities over various land cover types. However, there is a lack of microwave spectral emissivity libraries similar to those established in TIR region, e.g., the ASTER Spectral Library [106] and the MODIS UCSB Emissivity Library. Therefore, it is urgently needed to establish microwave spectral emissivity library of both laboratory and field measurements.

Development of New Physically-Based Methods for Retrieving LST from PMW Observations
Existing methods for retrieving LST from PMW observations are reviewed in Section 3. Among these methods, physically-based LST retrieval methods outperform the other methods in large-scale applications thanks to their clear physical meaning and general applicability. However, existing physically-based methods over-simplified the PMW RTE, leading to relatively low accuracy of LST retrieval.
Although PMW remote sensing has the capability of penetrating clouds, the cloud effects cannot be completely neglected, especially at high-frequency channels where microwave radiation is more sensitive to atmospheric effects [107]. The studies of Han et al. [52] showed that atmospheric effects at 36.5 and 89.0 GHz vary from approximately 3 K to 22 K over surfaces with different emissivities. The neglect or over-simplification of atmospheric effects could degrade the accuracy of PMW-derived LST. Therefore, it is necessary to develop a reasonable parameterization scheme for atmospheric effects.
In addition to atmospheric effects, the correction of surface emissivity effects is also crucial for LST retrieval from PMW observations. Some methods for PMW LST retrieval assume that surface emissivities have certain relationship at different frequencies and polarizations. Previous studies built the relationships between emissivities at different frequencies and polarizations using satellite-derived or model-based emissivities, but these relationships are based on empirical-statistical methods, which are lack of physical meaning [70,78,108,109]. It is necessary to investigate the intrinsic relationship between emissivities at different frequencies and polarizations using laboratory and field measurements. To promote the performance of LST retrieval from PMW observations, new physically-based methods should be developed by comprehensively accounting for the reasonable parameterization of atmospheric effects and the intrinsic relationship between emissivities at different frequencies and polarizations.

Simultaneous Retrieval of LST and Soil Moisture from PMW Observations
One of the main reasons for relatively low accuracy of LST retrieval from PMW observations is the large uncertainty in surface emissivity determination. Compared with surface emissivity in the TIR region, microwave emissivity varies strongly with soil moisture and shows large spatial, temporal, and spectral variability. Taking into account the advantages of PMW remote sensing with multi-frequency and dual-polarization characteristics, one of the most promising ways to improve the performance of LST retrieval methods is the simultaneous estimation of LST and soil moisture from PMW observations. Njoku and Li [110] made some attempts to retrieve LST, soil moisture, and vegetation water content simultaneously using PMW observations at 6-18 GHz. With the evolution of PMW remote sensing technique, more available frequencies are of great significance to make simultaneous retrieval of multi-parameter possible.
Two aspects should be paid more attention to the simultaneous retrieval of LST and soil moisture from PMW observations. First, the PMW RTE should be solvable, i.e., the number of unknowns should be no higher than the number of equations. It requires reducing the number of unknowns by combining some parameters into an effective parameter. For example, Pan et al. [111] and Zeng et al. [112] combined vegetation transmittance and surface roughness into one unknown because the two parameters appear as a product in the PMW RTE. Second, the solution of the PMW RTE should be stable. It requires choosing appropriate channels to perform the simultaneous retrieval according to the sensitivity of brightness temperature of the channels to the unknowns and the correlation between channels. Furthermore, the optimal initial guesses method, e.g., ANN, can be incorporated into the simultaneous retrieval of LST and soil moisture to make the solution of PMW RTE fast convergent.

Merging TIR-and PMW-Derived LST
LST can be retrieved from TIR and PMW satellite observations. TIR-derived LST has the advantages of higher retrieval accuracy (generally 1-2 K) and finer spatial resolution (e.g., 1 km for MODIS data), but it is excluded from cloudy-sky applications because of its incapability to penetrate clouds. Moreover, TIR-derived LST is surface skin temperature, which corresponds to the radiation emitted from a depth of the order of the penetration depth, i.e., of the order of the wavelength [113]. PMW satellite observations are applicable to retrieve LST in all-weather cases, but it has the limitations of low spatial resolution (e.g., 25 km for AMSR-E data) and low LST retrieval accuracy (generally 3-5 K). Furthermore, PMW-derived LST is subsurface temperature, which represents an equivalent soil temperature from specific depth below the surface to the land surface [114]. TIR-and PMW-derived LST is complementary. Therefore, combining the advantages of these two LST is promising for producing an all-weather LST product at a high spatial resolution [23,[115][116][117][118][119][120][121].
Two main concerns should be taken into account in the merging of TIR-and PMW-derived LST. First, PMW-derived subsurface temperature should be converted to surface skin temperature, which corresponds to LST retrieved from TIR data [122][123][124]. To better combine the TIR-and PMW-derived LST, it is indispensable to establish a physically-based model to convert subsurface temperature to surface skin temperature in terms of soil heat conduction equation [95,125]. Second, PMW-derived LST should be downscaled to the scale of TIR data for the spatial match between TIR-and PMW-derived LST. Up to now, there are no robust methods for PMW-derived LST downscaling.

Validation of PMW-Derived LST
Validation is the procedure of independently quantifying the uncertainty of satellite-derived products. Without validation, any method, model, and product cannot be confidently used. The biggest issue of PMW-derived LST validation is the spatial representativeness of ground-based LST measurements at satellite pixel scale. Although ground-based LST is regarded as the most reliable way for PMW-derived LST validation, it suffers from the difficulty in representing LST at the satellite pixel scale, which is normally tens of kilometers for PMW observations, especially in the cases when the surface is heterogeneous or fluctuant. To better validate PMW-derived LST, multiple-site measurements with an optimal spatial sampling scheme over relatively homogeneous areas, e.g., deserts and grasslands, are necessary. Furthermore, ground-based LST measurements are critical for validation of PMW-derived LST under cloudy conditions taking into account that clouds lead to relatively homogeneous LST over extended areas [77,78,126].