Retrieval of All-Weather 1 km Land Surface Temperature from Combined MODIS and AMSR2 Data over the Tibetan Plateau

: Land surface temperature (LST) is one of the most valuable variables for applications relating to hydrological processes, drought monitoring and climate change. LST from satellite data provides consistent estimates over large scales but is only available for cloud-free pixels, greatly limiting applications over frequently cloud-covered regions. With this study, we propose a method for estimating all-weather 1 km LST by combining passive microwave and thermal infrared data. The product is based on clear-sky LST retrieved from Moderate-resolution Imaging Spectroradiometer (MODIS) thermal infrared measurements complemented by LST estimated from the Advanced Microwave Scanning Radiometer Version 2 (AMSR2) brightness temperature to ﬁll gaps caused by clouds. Terrain, vegetation conditions, and AMSR2 multiband information were selected as the auxiliary variables. The random forest algorithm was used to establish the non-linear relationship between the auxiliary variables and LST over the Tibetan Plateau. To assess the error of this method, we performed a validation experiment using clear-sky MODIS LST and in situ measurements. The estimated all-weather LST approximated MODIS LST with an acceptable error, with a coefﬁcient of correlation (r) between 0.87 and 0.99 and a root mean square error (RMSE) between 2.24 K and 5.35 K during the day. At night-time, r was between 0.89 and 0.99 and the RMSE was between 1.02 K and 3.39 K. The error between the estimated LST and in situ LST was also found to be acceptable, with the RMSE for cloudy pixels between 5.15 K and 6.99 K. This method reveals a signiﬁcant potential to derive all-weather 1 km LST using AMSR2 and MODIS data at a regional and global scale, which will be explored in the future.


Introduction
Land surface temperature (LST) is an important environmental variable that controls land surface energy exchanges and water balance [1][2][3][4]. LST data sets are essential for a wide range of applications in urban heat islands [5][6][7][8], drought monitoring [9][10][11], climate change [12,13], hydrological processes [14][15][16] and crop yield estimation [17,18]. Traditionally, LST is calculated from ground-measured radiation, but this approach is timeconsuming, labor-intensive and does not provide consistent LST measurements over large areas. The development of remote sensing technologies made it possible to consistently estimate LST at regional to global scales at high temporal and spatial resolutions.
Thermal infrared (TIR) sensors are widely used to retrieve satellite-based LST. There are several TIR LST retrieval algorithms, such as the single-window algorithm, the splitwindow algorithm and the multi-channel algorithm [19][20][21][22]. TIR-derived LST has a relatively low error (approximately 0.3-2 K) and moderate spatial resolution (approximately 1 km). However, because the thermal infrared signal cannot penetrate clouds, LST products

Study Area
The Tibetan Plateau is located in Central Asia, between 26 • 00 N~39 • 47 N and 73 • 19 E~104 • 47 E. With an area of about 2.5 × 10 6 km 2 and an average altitude >4000 m, it is the largest and highest plateau in the world. The terrain is higher in the west and lower in the east. The annual average air temperature in the central region of the plateau is below 0 • C. The Tibetan Plateau belongs to the alpine climate zone and is highly sensitive to global climate change. LST is of great significance to the study of climate change in the plateau and the world. Because of the complex topography and the harsh natural environment, ground measurements of LST are very difficult to take here. At present, remote sensing is Remote Sens. 2021, 13,4574 3 of 20 the main way to obtain spatiotemporally continuous LST measurements here. In recent years, soil temperature and moisture observation networks have been established in the Tibetan Plateau that can be used to validate remote sensing LST data in this region. Figure 1 shows the study area and soil temperature sites.
Remote Sens. 2021, 13,4574 3 of 21 plateau and the world. Because of the complex topography and the harsh natural environment, ground measurements of LST are very difficult to take here. At present, remote sensing is the main way to obtain spatiotemporally continuous LST measurements here. In recent years, soil temperature and moisture observation networks have been established in the Tibetan Plateau that can be used to validate remote sensing LST data in this region. Figure 1 shows the study area and soil temperature sites.

Satellite Observations
AMSR2 is the second generation of the Advanced Microwave Scanning Radiometer carried on the Global Change Observation Mission-Water (GCOM-W). GCOM-W was successfully launched by the Japan Aerospace Exploration Agency (JAXA) on 18 May 2012. In July of the same year, JAXA began to release global microwave brightness temperature observation data. The orbit height of AMSR2 is 700 km and its scanning width is 1450 km. Its radiometer includes 14 horizontal/vertical polarization channels at 6.9, 7.3, 10.65, 18.7, 23.8, 36.5 and 89.0 GHz. AMSR2 is mainly used to monitor global water distribution and energy cycles [46,47]. AMSR2 L3 10 km daily brightness temperature can be downloaded from the JAXA website https://gportal.jaxa.jp (accessed on 10 September 2019).
Aqua is an earth observation satellite that was launched by the National Aeronautics and Space Administration (NASA) on 4 May 2002. MODIS is an important sensor on Aqua. The overpass times of Aqua are 1:30 am and 1:30 pm local solar time, which are the same as those of AMSR2. The MODIS/Aqua LST daily L3 global 1 km grid product  chooses the best available pixel value from all the acquisitions from the 16 day period. The criteria used are low clouds, low view angle and the highest NDVI value. Daily or eightday composite products were excluded because they contain spatial gaps. Shuttle Radar Topography Mission (SRTM) 90 m elevation data were downloaded from the geospatial data cloud platform of the computer network information center of the Chinese Academy of Sciences http://www.gscloud.cn (accessed on 10 September 2019). Table 1 shows the details of the datasets used in this paper. The first four datasets were used to produce all-weather 1 km LST. The last two datasets were used to find ground station whose soil temperature can reasonably represent the 1 km LST of the station pixel. Land use data V1.0 of the Tibetan Plateau were obtained by fusing and revising many datasets, including GlobCover of the European Space Agency, MCD12Q1 of NASA, UMD Land Cover, and IGBP DISCover of USGS. Monthly vegetation index data from China were derived from SPOT/VEGETATION NDVI data by the maximum value composite method, which has been recording the monthly vegetation index of China since 1998 at a spatial resolution of 30 m. Land use and monthly vegetation index datasets were downloaded from the national Tibetan Plateau scientific data center https://data.tpdc.ac.cn/zh-hans/ (accessed on 10 September 2019). In order to compare RF with other methods, we downloaded the GEOS-5 soil temperature (ST) dataset from https://portal.nccs.nasa.gov/datashare/gmao_ops/pub/fp/das/ (accessed on 10 September 2019). The GEOS-5 FP (forward processing) products are generated by the GMAO (Global Modeling and Assimilation Office), and provide forecasts and assimilation products for near-real-time production. The ST dataset in GEOS-5 FP is hourly, at a spatial resolution of 0.3125 • by 0.25 • . In this study, 0-5 cm soil temperature was selected. We choose this ST dataset as a reference because Ma et al. (2021) had used ground soil temperature observations of approximately 800 stations worldwide to fully assess six model-and satellite-based surface ST products from April 2015 to December 2017; their results showed that the GEOS-5 exhibits the smallest averaged unbiased root mean square difference of 1.84 K [48]. All datasets mentioned above are from 1 May 2015 to 30 April 2016.

In Situ Measurements
To evaluate the error of the proposed method, the soil temperature dataset was also downloaded from the national Tibetan Plateau scientific data center https://data.tpdc. ac.cn/zh-hans/ (accessed on 10 September 2019). The soil temperature and moisture observation dataset of the Tibetan Plateau (2008-2016) records soil information every hour. The multi-scale soil temperature and moisture observation network dataset of the central Tibetan Plateau (2010-2016) records soil information every 30 min. In this study, 0-5 cm soil temperature was used to validate the all-weather 1 km LST product. Although soil surface temperature and LST are slightly different, the in-situ soil surface temperature can reasonably represent the LST in the 1 km grid domain because of the following reasons: (1) the land cover types of the sites are consistent with the station pixels; (2) the maximum NDVI of the station pixels is no more than 0.18 in 2015 and 2016, indicating there is no vegetation or little vegetation in these pixels; (3) the data acquisition time of the station is close to the transit time of the satellite [42]. Table 2 shows the details of each in situ station used in this study and their locations are shown in Figure 1. In this study, vertical polarization (v-pol) brightness temperatures at 10.65, 18.7 and 36.5 GHz of AMSR2 were used. As heavy rainfall, ice or snow have a great influence on the brightness temperature, the corresponding pixels need to be filtered. The Scattering Index (SI; Equation (1)) was used to find the pixels affected by heavy rainfall [49]: where T b,19v , T b,23v and T b,89v denote brightness temperatures at v-pol at 18.7, 23.8 and 89.0 GHz. Pixels with SI > 10 were identified as those affected by heavy rainfall. Equation (2) was used to filter the pixels affected by ice and snow [50]: where T b,37v denotes brightness temperature at v-pol at 36.5 GHz.
To ensure that all datasets were consistent with each other, the MODIS products sinusoidal projection was resampled to AMSR2 geographical projection using the MODIS reprojection tool developed by NASA. To estimate the passive microwave LST under a high resolution, the brightness temperature data of AMSR2 were resampled from 10 km to 1 km using the cubic convolution method [38]. In addition, we tested the method by comparing the correlation coefficient (r) between brightness temperature (36.5 GHz) and MODIS LST before and after cubic interpolation ( Figure 2); r changed from 0.62 to 0.61, which showed that cubic interpolation is suitable in the case of brightness temperature. LST data (LST_Day_1km, LST_Night_1km) and quality control data (QC_Day, QC_Night) were extracted from MYD11A1. The pixels with quality control data of "cloud", "average emission error > 0.04" or "average LST error > 3 K" were regarded as cloud-affected pixels and were removed. The NDVI data (1km_16_days_NDVI) were extracted from MYD13A2. SRTM DEM was also resampled from 90 m to a MODIS resolution of 1 km.

RF-Based LST Estimation
The spatial and temporal dynamics of LST are very complex since they depend on a variety of land surface variables, such as land cover and topography. The v-pol brightness temperature has been shown to be more sensitive to LST than the h-pol brightness temperature [44]. In addition, the 37 GHz v-pol brightness temperature shows a good linear relationship with LST and is thus considered the best microwave band for LST retrieval. However, it is not ideal to retrieve LST using 37 GHz v-pol brightness temperature in sparse and open shrubs because the brightness temperature in 37 GHz is affected by vegetation and soil water [24]. To identify the best variables for the RF regression, we created different groups of variables using 10.65 GHz,18.7 GHz and 36.5 GHz v-pol brightness temperatures. Each group also included DEM and NDVI, and half of them also included the microwave polarization difference index (MPDI). For each group, we performed 10-fold cross-validation, and the error was evaluated using the root mean square error (RMSE; Table 3). The plot of feature importance shows that the 37 GHz v-pol brightness temperature is of the highest importance and plays a leading role in LST retrieval ( Figure 3). The second most important variable is DEM, which describes the influence of topography on the LST in the Tibetan Plateau. We found that when the number of estimators is 131 and the max depth is 39, the RF yields the smallest RMSE. Consequently, the 36.5 GHz v-pol brightness temperature was used to constrain the range of LST, and the 10.65 GHz and 18.7 GHz v-pol brightness temperature were used to correct the effects of vegetation and soil moisture. different groups of variables using 10.65 GHz,18.7 GHz and 36.5 GHz v-pol brightness temperatures. Each group also included DEM and NDVI, and half of them also included the microwave polarization difference index (MPDI). For each group, we performed 10fold cross-validation, and the error was evaluated using the root mean square error (RMSE; Table 3). The plot of feature importance shows that the 37 GHz v-pol brightness temperature is of the highest importance and plays a leading role in LST retrieval ( Figure  3). The second most important variable is DEM, which describes the influence of topography on the LST in the Tibetan Plateau. We found that when the number of estimators is 131 and the max depth is 39, the RF yields the smallest RMSE. Consequently, the 36.5 GHz v-pol brightness temperature was used to constrain the range of LST, and the 10.65 GHz and 18.7 GHz v-pol brightness temperature were used to correct the effects of vegetation and soil moisture.     RF has been widely used in classification and nonlinear regression problems in remote sensing [51,52]. RF is robust to outliers and to unbalanced data, ensuring efficient performance. The proposed RF method is based on the statistical relationship between LST and auxiliary variables and involves the following two phases ( Figure 4): (1) training and (2) prediction. During training, we ensured that the brightness temperature, MODIS LST, NDVI and DEM had the same spatial resolution and that pixels acquired under clear- RF has been widely used in classification and nonlinear regression problems in remote sensing [51,52]. RF is robust to outliers and to unbalanced data, ensuring efficient performance. The proposed RF method is based on the statistical relationship between LST and auxiliary variables and involves the following two phases ( Figure 4): (1) training and (2) prediction. During training, we ensured that the brightness temperature, MODIS LST, NDVI and DEM had the same spatial resolution and that pixels acquired under clear-sky conditions were used to train the retrieval model with both ascending and descending acquisition. In the prediction phase, the model built through training was applied to the resampled brightness temperature and other auxiliary data to predict the 1 km LST for cloudy pixels. RF has been widely used in classification and nonlinear regression problems i mote sensing [51,52]. RF is robust to outliers and to unbalanced data, ensuring effi performance. The proposed RF method is based on the statistical relationship betw LST and auxiliary variables and involves the following two phases ( Figure 4): (1) tra and (2) prediction. During training, we ensured that the brightness temperature, MO LST, NDVI and DEM had the same spatial resolution and that pixels acquired under c sky conditions were used to train the retrieval model with both ascending and descen acquisition. In the prediction phase, the model built through training was applied t resampled brightness temperature and other auxiliary data to predict the 1 km LS cloudy pixels.   At present, there are two ways to retrieve high-resolution LST from passive microwave brightness temperature, namely, inversion-disaggregation and disaggregation-inversion. The first way uses the brightness temperature to retrieve low-resolution LST, then downscales the LST to a high resolution. The second way disaggregates the brightness temperature to a MODIS resolution, then uses the resampled brightness temperature to retrieve high-resolution LST. It should be noted that the process of retrieving LST from AMSR2 introduces great uncertainty, which is inherited by the downscaled results. Therefore, we compared these two methods at a 10 km resolution. First, we trained the RF model at a 10 km resolution and estimated the 10 km LST. Then, the 1 km resolution LST estimated by the RF disaggregation-inversion method was aggregated to the resolution of AMSR2. Finally, we calculated r and RMSE between the estimated 10 km LST and MODIS LST. As shown in Figure 5, the second method had smaller errors; thus, disaggregation-inversion was chosen in the present study.

Random Forest
fore, we compared these two methods at a 10 km resolution. First, we trained the RF at a 10 km resolution and estimated the 10 km LST. Then, the 1 km resolution LS mated by the RF disaggregation-inversion method was aggregated to the resolu AMSR2. Finally, we calculated r and RMSE between the estimated 10 km LST and M LST. As shown in Figure 5, the second method had smaller errors; thus, disaggreg inversion was chosen in the present study.

LST Fusion
The 1 km LST estimated can be merged directly with MODIS LST. After the fu 1 km LST for cloudy pixels and MODIS LST, some areas still had missing data beca the AMSR2 scanning gaps. In this study, we first used a temporally neighboring in lation, and then used a spatially neighboring interpolation to supplement the gap. B the gap is particularly wide at low latitudes, we were unable to obtain an acceptabl if we used spatially neighboring interpolation directly. Temporally neighboring in lation assumed that the cloud coverage was the same on adjacent days, and used th of the day before and after to interpolate the LST of a given date. If the data w missing after the temporally neighboring interpolation, spatially neighboring inte tion was performed. The procedure for the estimation method of 1 km LST is sho Figure 6.

LST Fusion
The 1 km LST estimated can be merged directly with MODIS LST. After the fusion of 1 km LST for cloudy pixels and MODIS LST, some areas still had missing data because of the AMSR2 scanning gaps. In this study, we first used a temporally neighboring interpolation, and then used a spatially neighboring interpolation to supplement the gap. Because the gap is particularly wide at low latitudes, we were unable to obtain an acceptable error if we used spatially neighboring interpolation directly. Temporally neighboring interpolation assumed that the cloud coverage was the same on adjacent days, and used the LST of the day before and after to interpolate the LST of a given date. If the data was still missing after the temporally neighboring interpolation, spatially neighboring interpolation was performed. The procedure for the estimation method of 1 km LST is shown in Figure 6. mated by the RF disaggregation-inversion method was aggregated to the resolution of AMSR2. Finally, we calculated r and RMSE between the estimated 10 km LST and MODIS LST. As shown in Figure 5, the second method had smaller errors; thus, disaggregationinversion was chosen in the present study.

LST Fusion
The 1 km LST estimated can be merged directly with MODIS LST. After the fusion of 1 km LST for cloudy pixels and MODIS LST, some areas still had missing data because of the AMSR2 scanning gaps. In this study, we first used a temporally neighboring interpolation, and then used a spatially neighboring interpolation to supplement the gap. Because the gap is particularly wide at low latitudes, we were unable to obtain an acceptable error if we used spatially neighboring interpolation directly. Temporally neighboring interpolation assumed that the cloud coverage was the same on adjacent days, and used the LST of the day before and after to interpolate the LST of a given date. If the data was still missing after the temporally neighboring interpolation, spatially neighboring interpolation was performed. The procedure for the estimation method of 1 km LST is shown in Figure 6.

Still exist missing value?
Spatio-temporal interpolation All-weather LST 1km LST on cloudy pixels

MODIS LST on clear days
Merge Figure 6. Procedure for merging estimated LST on cloudy pixels with MODIS LST. Figure 6. Procedure for merging estimated LST on cloudy pixels with MODIS LST.
Inverse distance weighting (IDW), regular kriging (Kriging) and regular spline function (Spline) are three widely used spatial interpolation methods. To find better spatially neighboring interpolation methods, the errors of three spatial interpolation methods in different land use types of the Tibetan Plateau were compared. Firstly, we used land cover products as criteria to obtain MODIS LST data under every land cover type. Secondly, we assumed the same cloud coverage area for every land cover type, and used three spatial interpolation methods to estimate the LST value under the assumed cloud coverage area. Finally, we calculated the RMSE between the estimated LST and original LST (Figure 7). For urban, cropland, water, forest, bare land and grass, Spline had the smallest error. For snow/ice, Kriging had the smallest error. As Spline interpolation had, in general, smaller errors (0.10-2.73 K) under different land use/covers, it was selected.
products as criteria to obtain MODIS LST data under every land cover type. Secondly, we assumed the same cloud coverage area for every land cover type, and used three spatial interpolation methods to estimate the LST value under the assumed cloud coverage area. Finally, we calculated the RMSE between the estimated LST and original LST (Figure 7). For urban, cropland, water, forest, bare land and grass, Spline had the smallest error. For snow/ice, Kriging had the smallest error. As Spline interpolation had, in general, smaller errors (0.10-2.73 K) under different land use/covers, it was selected.

Algorithm Evaluation Strategy
Because the LST retrieved by AMSR2 has a large uncertainty, it is necessary to evaluate its error to better evaluate the RF algorithm used in this paper, using r and RMSE as performance metrics. The AMSR2 LST under clear-sky pixels was compared to MODIS LST and in situ soil temperature. The AMSR2 LST under cloudy sky pixels was only compared to soil temperature because there was no MODIS LST for cloudy sky pixels. In addition, we chose the ST dataset from the GEOS-5 FP products as the reference to evaluate the RF method in this paper using the methods of Holmes, Zhao, and Zeng [24,44,53]. Since the datasets are spatially inconsistent (0.3125° by 0.25° for GEOS-5 ST, 1 km for the RF method, 10 km for the other three methods), the estimated LST was resampled to 0.25° by 0.25° by using the nearest-neighbor interpolation method, which is a feasible sampling approach and is extensively adopted in existing validation studies [54][55][56][57].
Holmes et al., (2009) developed a simple linear relationship of 37 GHz v-pol brightness temperature to obtain LST. The deviation of this method is within 1 K for 70% of the vegetated pixels. Due to its simple expression, it has been widely used as a soil temperature retrieval algorithm and in corresponding research [58][59][60]. It can be expressed as:

Algorithm Evaluation Strategy
Because the LST retrieved by AMSR2 has a large uncertainty, it is necessary to evaluate its error to better evaluate the RF algorithm used in this paper, using r and RMSE as performance metrics. The AMSR2 LST under clear-sky pixels was compared to MODIS LST and in situ soil temperature. The AMSR2 LST under cloudy sky pixels was only compared to soil temperature because there was no MODIS LST for cloudy sky pixels. In addition, we chose the ST dataset from the GEOS-5 FP products as the reference to evaluate the RF method in this paper using the methods of Holmes, Zhao, and Zeng [24,44,53]. Since the datasets are spatially inconsistent (0.3125 • by 0.25 • for GEOS-5 ST, 1 km for the RF method, 10 km for the other three methods), the estimated LST was resampled to 0.25 • by 0.25 • by using the nearest-neighbor interpolation method, which is a feasible sampling approach and is extensively adopted in existing validation studies [54][55][56][57].
Holmes et al. (2009) developed a simple linear relationship of 37 GHz v-pol brightness temperature to obtain LST. The deviation of this method is within 1 K for 70% of the vegetated pixels. Due to its simple expression, it has been widely used as a soil temperature retrieval algorithm and in corresponding research [58][59][60]. It can be expressed as: where T b,37v is the brightness temperature of 36.5 GHz v-pol. Zeng et al. (2014) established the LST inversion model in the Tibetan Plateau using ascending and descending AMSR-E brightness temperature. It can be expressed with Equations (4) and (5):  (6). We found that there is 4-digit precision in coefficient values of Equation (6). To assess whether 4-digit precision is necessary, we have conducted an additional experiment; the experiment is based on the MODIS LST and AMSR2 brightness temperature. Firstly, we used 4-digit, 3-digit, 2-digit and 1-digit coefficients to estimate the 10 km LST. Then, we resampled the 1 km MODIS LST into 10 km. Finally, we calculated the correlation coefficients (R) and RMSE between the estimated LST and MODIS LST. The RMSE value of 3-digit coefficients increased by 0.02 K compared with the 4-digit coefficients, which clearly showed that 4-digit precision in coefficients is necessary. where T Bv denotes the brightness temperature at the corresponding frequency. MPDI can be expressed with Equation (7), where a is 100 and f denotes frequency: There are a lot of alternatives to random forest; to make sure that the choice of random forest is justified, we made a comparison between the classic MLR-based methodology and regression techniques in machine learning (neural networks, nearest neighbor, ensembles such as random forests). The evaluated algorithms included MLR (multiple linear regression), MLP (multilayer perceptron), CART (Classification and Regression Tree), KNR (K-nearest neighbor regression), GBRT (gradient boost regression tree) and RF (random forest). These are representative regression models; however, few studies have been conducted to compare their quality in predicting LST [38,61,62]. For each model, we did 10-fold cross-validation to find the optimized parameters. Then, we calculated the R 2 and RMSE of each model under the optimized parameters. LST had some missing data due to cloud contamination, and was thus spatio-temporally discontinuous. Compared with the original MODIS LST, the results are spatio-temporally continuous. Spatially, the LST on the second row continues the spatial variation trend of the original MODIS LST. The LST decreases from southeast to northwest, which is consistent with the terrain of the Tibetan Plateau, which is higher in the northwest and lower in the southeast. Temporally, the all-weather 1 km LST captures well the monthly variation of temperature in the Tibetan Plateau, with a higher LST in summer (from June to August).

Different Methods Comparison
Based on the descending orbit data of 2 June 2015, we calculated the probability density histogram of the difference between the estimated LST and GEOS-5 LST (Figure 9). The distribution of the difference is close to the normal distribution for each method. The expectations (µ) for Holmes, Zhao, and our method are less than zero, which indicates an underestimation. In contrast, the µ of the Zeng method is 5.28, which indicates an overestimation. In addition, the µ of the Holmes method (µ = −2.23) and the method in this paper (µ = −2.84) are closer to zero. In terms of standard deviation (σ), the method in this paper yields the smallest σ of 6.19, which shows that the errors are more centralized.
In conclusion, compared with other methods, the method in this paper yields better expectation and standard deviation, with a difference that is more centralized around zero. expectations ( ) for Holmes, Zhao, and our method are less than zero, which indicates an underestimation. In contrast, the of the Zeng method is 5.28, which indicates an overestimation. In addition, the of the Holmes method ( = −2.23) and the method in this paper ( = −2.84) are closer to zero. In terms of standard deviation ( ), the method in this paper yields the smallest of 6.19, which shows that the errors are more centralized. In conclusion, compared with other methods, the method in this paper yields better expectation and standard deviation, with a difference that is more centralized around zero.  We calculated the r and RMSE between the estimated LST/MODIS LST and GEOS-5 LST. As shown in Figure 10, the Holmes method yields the worst result (r = 0.66 and RMSE = 12.48 K). This was followed by the Zhao method, which had an r of 0.33 and an RMSE of 12.15 K. The Holmes method and the Zeng method have the same r (0.66); the Zeng method has a smaller RMSE (8.89 K), indicating that the Zeng method performs better-perhaps because it is an improvement of the Holmes method. As the Zeng method is based on the 0-5 cm in situ soil temperature data from the Tibetan Plateau, it may be more suitable for the Tibetan Plateau. In conclusion, the RMSE values of the Holmes and the Zeng methods are both large: 12.48 K and 8.89 K, respectively. The r and RMSE of the method in this paper are 0.5 and 6.81 K; the error is close to that of the MODIS LST, showing that the RF method predicts LST with the smallest error compared with other methods. This could be because the RF method in this paper used the MODIS LST for training and therefore the LST it predicted was more similar to the GEOS-5 LST; in contrast, other methods used ground measurements as the true values.  Table 4 is the comparison between the classic MLR-based methodology and regression techniques in machine learning. The results confirm that classic MLR is outperformed by machine learning techniques and concretely, our experiments suggest that random forest regression outperforms the rest of the techniques.

Models
Optimized Parameters R 2 RMSE, K Training Time  Table 4 is the comparison between the classic MLR-based methodology and regression techniques in machine learning. The results confirm that classic MLR is outperformed by machine learning techniques and concretely, our experiments suggest that random forest regression outperforms the rest of the techniques.

Evaluation Using MODIS LST
Since MODIS LST products have been evaluated to have an error of less than 1 K at most areas, the AMSR2 LST predicted by RF were compared with the MODIS LST for validation using r and RMSE as performance metrics. We selected the dates of 1 July 2015, 1 October 2015, 1 January 2016 and 1 April 2016 to represent the summer, autumn, winter, and spring of the Tibetan Plateau, respectively. Figure 11 shows the scatter plot of predicted LST and MODIS LST. At night-time, the r was between 0.92 and 0.97, and the RMSE was between 1.16 K and 2.4 K. At daytime, the r was between 0.9 and 0.93, and the RMSE was between 2.35 K and 4.67 K. This shows that the night-time error is smaller than that of the daytime. This may be because the spatial heterogeneity of LST at night is smaller than that during the day, so the relationship between the LST and its predictors (NDVI, DEM, AMSR2 BT) is simpler and the RF model we obtain is more accurate. The figure also shows that the prediction error is smaller in autumn, winter, and spring than that in summer. The reason may be that there is more rainfall in summer, and the microwave signal is interfered with by the water vapor in the atmosphere, so the radiation transmission process is more complex.
To understand the error for the entire year, we calculated the r and RMSE between the predicted LST and MODIS LST for every day and night from 1 May 2015 to 30 April 2016. The predicted LST at night has a higher correlation with the MODIS LST ( Figure 12). During the night-time, most r values are close to 0.94. During the daytime, most r values are close to 0.92. The predicted LST at night has a smaller RMSE with the MODIS LST ( Figure 13). During the nighttime, the RMSE is in the range between 1.02 K and 3.39 K. During the daytime, the RMSE is in the range between 2.24 K and 5.35 K. The r and RMSE values can be compared with the r and RMSE values of previous machine learning methods. The ANN method Shwetha used to predict high spatio-temporal resolution land surface temperatures under cloudy conditions shows r values from 0.56 to 0.90 and RMSE values from 1.70 K to 2.12 K during the night-time, and has r values from 0.78 to 0.96 and RMSE values from 1.86 K to 4.00 K during the daytime [38]. Overall, there is a significant r and a small RMSE between the predicted LST and the MODIS LST both during the day and at night, indicating the ability of RF in predicting the LST.

Evaluation Using In Situ Measurements
The LST estimated by AMSR2 brightness temperature and RF was validated with a 0-5 cm soil temperature. The soil temperature was measured by a high-precision temperature measuring instrument. According to the station screening rules mentioned in the data preprocessing, the in situ soil surface temperature reasonably represents the 1 km LST in this study because these stations have an acceptable spatial heterogeneity. The stations are PL01, PL03, PL05, PL11, PL12, AL02 and SQ16. Figure 14 shows the scatter plot between the estimated LST and the in situ soil surface temperature. The estimated LST is close to the in situ surface temperature, and the RMSE is from 5.60 K to 10.56 K.
The retrieved AMSR2 LST under clear and cloudy conditions were respectively compared with the in-situ soil temperature. There is a scale difference between the in situ soil temperature and the 1 km LST, so the error caused by scale mismatch should be measured. Because the MODIS LST products have been verified to have an error of less than 1 K, this paper calculated the RMSE between the MODIS LST and in situ soil temperature to represent the error caused by the scale mismatch. Figure 15 shows the scatter plot between the in situ soil temperature and MODIS LST, and the AMSR2 LST under clear and cloudy conditions. Compared to the MODIS LST, the LST retrieved by RF is closer to the in-situ soil temperature. The retrieved LST under cloudy skies is closer to the in situ soil temperature than that under clear skies.
The LST estimated by AMSR2 brightness temperature and RF was validated with a 0-5 cm soil temperature. The soil temperature was measured by a high-precision temperature measuring instrument. According to the station screening rules mentioned in the data preprocessing, the in situ soil surface temperature reasonably represents the 1 km LST in this study because these stations have an acceptable spatial heterogeneity. The stations are PL01, PL03, PL05, PL11, PL12, AL02 and SQ16. Figure 14 shows the scatter plot between the estimated LST and the in situ soil surface temperature. The estimated LST is close to the in situ surface temperature, and the RMSE is from 5.60 K to 10      The retrieved AMSR2 LST under clear and cloudy conditions were respectively compared with the in-situ soil temperature. There is a scale difference between the in situ soil temperature and the 1 km LST, so the error caused by scale mismatch should be measured.  The retrieved AMSR2 LST under clear and cloudy conditions were respectively compared with the in-situ soil temperature. There is a scale difference between the in situ soil temperature and the 1 km LST, so the error caused by scale mismatch should be measured. Because the MODIS LST products have been verified to have an error of less than 1 K, this paper calculated the RMSE between the MODIS LST and in situ soil temperature to represent the error caused by the scale mismatch. Figure 15 shows the scatter plot between the in situ soil temperature and MODIS LST, and the AMSR2 LST under clear and cloudy from 0.385 (PL05) to 0.908 (AL02) and RMSE ranged from 5.132 K (AL02) to 12.349 K (SQ16). For the AMSR2 LST under a cloudy sky, the r ranged from 0.324 (PL12) to 0.745 (SQ16) and the RMSE ranged from 5.152 K (PL01) to 6.995 K (AL02). Overall, Table 5 indicates that the estimated AMSR2 LST has a similar error as the MODIS LST. In addition, the RMSE of the AMSR2 LST under a cloudy sky is similar with that under a clear skythis shows that RF can yield the same satisfactory accuracy under a cloudy sky as under a clear sky. The RMSE of the AMSR2 LST under a cloudy sky (5.152-6.995 K) is comparable to the RMSE of the ANN method proposed by Shwetha (2.9 K to 6.2 K for a cloudy sky) [38].   The r and RMSE between the remote sensing LST and in situ soil temperature were calculated for each in situ station for both clear and cloudy sky pixels (Table 5). For MODIS LST under clear-sky pixels, r ranged from 0.324 (PL12) to 0.926 (AL02) and RMSE ranged from 5.398 K (AL02) to 13.143 K (SQ16). For the AMSR2 LST under a clear sky, r ranged from 0.385 (PL05) to 0.908 (AL02) and RMSE ranged from 5.132 K (AL02) to 12.349 K (SQ16). For the AMSR2 LST under a cloudy sky, the r ranged from 0.324 (PL12) to 0.745 (SQ16) and the RMSE ranged from 5.152 K (PL01) to 6.995 K (AL02). Overall, Table 5 indicates that the estimated AMSR2 LST has a similar error as the MODIS LST. In addition, the RMSE of the AMSR2 LST under a cloudy sky is similar with that under a clear sky-this shows that RF can yield the same satisfactory accuracy under a cloudy sky as under a clear sky. The RMSE of the AMSR2 LST under a cloudy sky (5.152-6.995 K) is comparable to the RMSE of the ANN method proposed by Shwetha (2.9 K to 6.2 K for a cloudy sky) [38].

Conclusions
In this paper, a machine learning method was proposed to estimate all-weather 1 km LST over the Tibetan Plateau. The method resampled the brightness temperature of AMSR2 into a 1 km resolution, then used RF regression to define a nonlinear relationship between MODIS LST and surface variables to estimate the 1 km LST under cloudy skies. Following LST fusion and interpolation, all-weather 1 km LST data was generated.
In situ measurements and MODIS LST were used to evaluate the algorithm. The results show that the LST predicted by RF is consistent with the MODIS LST and 0-5 cm soil temperature. Additionally, the method fills the gaps caused by cloud cover, which can lead to incomplete data over a large area. The gap-filled LST captures the spatiotemporal variation trend of MODIS LST. This method makes it possible to obtain all-weather 1 km LST datasets, which can support regional and even global climate change research. Nevertheless, there is still room for further improvement of the method regarding the following aspects: (1) an AMSR2-derived LST downscaling method should be developed; downscaling is an important research topic that was simplified for this paper. State-of-the-art downscaling approaches should be further explored; (2) a more accurate time interpolation method is needed, which may use meteorological or other satellite data to decrease errors; (3) although this study has focused on the Tibetan Plateau, we believe that this method could be used over other regions-but additional validation will be needed.