Estimation of Surface and Near-Surface Air Temperatures in Arid Northwest China Using Landsat Satellite Images

Near-surface air (Ta) and land surface (Ts) temperatures are essential parameters for research in the fields of agriculture, hydrology, and ecological changes, which require accurate datasets with different temporal and spatial resolutions. However, the sparse spatial distribution of meteorological stations in Northwest China may not effectively provide high-precision Ta data. And it is not clear whether it is necessary to improve the accuracy of Ts which has the most influence on Ta. In response to this situation, the main objective of this study is to estimate Ta for Northwest China using multiple linear regression models (MLR) and random forest (RF) algorithms, based on Landsat 8 images and auxiliary data collected from 2014 to 2019. Ts, NDVI (Normalized Difference Vegetation Index), surface albedo, elevation, wind speed, and Julian day were variables to be selected, then used to estimate the daily average Ta after analysis and adjustment. Also, the Radiative Transfer Equation (RTE) method for calculating Ts would be corrected by NDVI (RTE-NDVI). The results show that: 1) The accuracy of the surface temperature (Ts) was improved by using RTE-NDVI; 2) Both MLR and RF models are suitable for estimating Ta in areas with few meteorological stations; 3) Analyzing the temporal and spatial distribution of errors, it is found that the MLR model performs well in spring and summer, and is lower in autumn, and the accuracy is higher in plain areas away from mountains than in mountainous areas and nearby areas. This study shows that through appropriate selection and combination of variables, the accuracy of estimating the pixel-scale Ta from satellite remote sensing data can be improved in the area that has less meteorological data.


INTRODUCTION
Near-surface air temperature (Ta), usually refers to the air temperature at 2 m above the ground, is an essential factor affecting ecology, agriculture, and urban areas (Raja Reddy et al., 1997;Krüger and Emmanuel, 2013;Shamir and Georgakakos, 2014), and is also the basis for climate change studies (Alkama and Cescatti, 2016;Bathiany et al., 2018). The traditional method of obtaining Ta mainly relies on the temperature sensor installed at the meteorological station, and the interpolation method is frequently used to extend to regional-scale applications (Mostovoy et al., 2006). If the meteorological stations are sparse and unevenly distributed, the accuracy of the interpolation method will be greatly restricted due to the influence of underlying surface heterogeneity and heat conduction unevenness to air temperature (Chen et al., 2015). The air temperature has become the primary driving variable of many land surface models (Nieto et al., 2011), so its spatial fidelity must be higher than that obtained by interpolation of point observation data, and even the most complex geostatistical techniques cannot meet the requirements (Prince et al., 1998).
Satellite data has the characteristics of continuous spatial coverage (Czajkowski et al., 1997), which can obtain largescale atmospheric information and invert surface parameters, including global surface temperature, vegetation index, elevation, and other information Yoo et al., 2018). Due to the complexity of atmospheric radiation and its low proportion in remote sensing signals, Ta cannot be directly reversed (Xu et al., 2012;Li and Zha, 2018). However, Ta can be obtained by establishing the regression relationship between Ta and remote sensing inversion and auxiliary parameters, such as using land surface temperature (Ts), Normalized Difference Vegetation Index (NDVI), wind speed, geographic location and elevation, among which Ts is the most important parameter (Vancutsem et al., 2010;Hachem et al., 2012;Song and Wu, 2018).
Ts as the direct driving force of long-wave radiation and turbulent heat flux exchange at the surface-atmosphere interface is one of the most significant parameters in the physical process of surface energy and water balance at regional and global scales (Anderson et al., 2008;Li et al., 2013;Orhan et al., 2014;Folland et al., 2018). At present, there are three main methods for using Landsat to retrieve Ts: atmospheric correction method (also known as Radiative Transfer Equation: RTE), Single Channel Algorithm (SCA), and Split Window Algorithm (SWA). The SWA does not require any atmospheric profile information at the time of collection, which needs to use two thermal infrared channels . However, the United States Geological Survey has pointed out that Thermal Infrared Sensor (TIRS) 11th band has data reception abnormalities and calibration instability problems, which mainly affects the accuracy of the split window algorithm applied to Landsat-8 TIRS data to retrieve Ts (Xu, 2015). SCA and RTE rely on atmospheric transmissivity and upwelling and downwelling atmospheric radiances (Jimenez-Munoz et al., 2009;Sekertekin and Bonafoni, 2020). RTE removes the error caused by the atmosphere's thermal radiation on the surface and converts the thermal radiation intensity to the corresponding Ts (Ma and Pu, 2020). When using different data sets, the performance of the RTE method to retrieve Ts is also different. The Ts calculation result for Landsat TM 5 data is better than other methods in the same period with RMSE is 1-3°C (Sobrino et al., 2004;Ndossi and Avdan, 2016;Windahl and Beurs, 2016); however, the RMSE calculated based on Landsat 8 TIRS 10th band was 1.5-5°C, and most of them are inferior to other algorithms (Ndossi and Avdan, 2016;Wang et al., 2016;Sekertekin and Bonafoni, 2020). Moreover, the overestimation shown in the TIRS band will increase as the proportion of vegetation decreases (Xu and Huang, 2016).
Several studies have used surface information to estimate air temperatures, such as the temperature-vegetation index (TVX), energy balance, statistics, and machine learning methods (Zakšek and Schroedter-Homscheidt, 2009;Benali et al., 2012). Nemani and Running (1989), Goward et al. (1994) proposed the TVX approach to estimate near surface air temperature with promising results. The method is based on the assumption that there is a strong negative correlation between Ts and vegetation index (Goward et al., 1994;Czajkowski et al., 1997). Assuming that the Ta for fully covered vegetation is close to Ts, the value of full coverage NDVI (NDVImax) can be used to obtain an approximate value of Ta (Stisen et al., 2007;Nieto et al., 2011;Zhu et al., 2013). However, this assumption does not apply to all seasons, soil moisture, and ecosystem types, so estimation of Ta by using TVX method is not feasible in areas or seasons without high vegetation cover (Vancutsem et al., 2010).
The energy balance method has a physical mechanism, so it has well portability and versatility (Hou et al., 2013;Shen et al., 2020). According to the energy balance equation, Ta is related to surface temperature, and it depends on various environmental factors such as solar radiation, cloud cover, wind speed, soil moisture, and surface type (Prince et al., 1998). A large number of required parameters cannot be completely retrieved by remote sensing (Mostovoy et al., 2006), so it is difficult to use remote sensing to perform Ta inversion in the area. The issue of unclosed surface energy balance also brings additional uncertainty to this method (Zhang et al., 2015).
Statistical methods need to analyze the relationship between Ta and Ts and other auxiliary data, and then build an estimation model based on specific correlation (Cresswell et al., 1999;Park, 2011), which including simple statistical models, multiple linear regression (MLR) models, geographically weighted regression (GWR) models, and machine learning methods (Vogt et al., 1997;Vancutsem et al., 2010;Shen et al., 2020). Studies have shown that the linear regression models are more accurate in calculating the average daily temperature with a root mean square error (RMSE) ranging between 1.29-3.60°C (Chen et al., 2015;Shi et al., 2016;Yang et al., 2017), which can produce good results in a specific space and time range, but require a large amount of data involved in the calculation and training of the algorithms (Stisen et al., 2007). Geographically weighted statistical and machine learning methods usually have higher accuracy (Moser et al., 2015;Wang et al., 2017. Geographically and temporally weighted regression (GTWR) is an extension of the general linear regression model, which embeds changes in location and time into the regression equation and estimates regression coefficients for spatio-temporal variation by performing local regressions that can solve for constant-coefficient limits (Bai et al., 2016;. Machine learning methods can handle non-linear and highly correlated predictors (James et al., 2013) and estimate the temperature in areas with complex and heterogeneous underlying surfaces, mainly including neural networks (Jang et al., 2004), M5 model trees (Emamifar et al., 2013), and random forests , support vector machine (Moser et al., 2015). Random forests (RF) are widely used and have been verified in various terrains. Ho et al. (2014) indicated that the RF algorithm is very useful for mapping the variability of urban internal temperature. Meyer et al. (2016) have pointed out that compared with linear regression, generalized augmented regression model (GBM), and cubic regression, the RF algorithm performs poorly in the extremely cold Antarctica. However, Noi et al. (2017) have shown reliable results in mountainous areas. Therefore, the conclusion of RF estimation Ta still needs to be discussed.
Most of the current researches for estimating Ta is based on the MODIS due to its time continuity advantage, but its resolution cannot meet the demand for farmland or even smaller scales. The spatial resolution of Landsat 8 imagery is higher than that of MODIS, and the estimated Ta data have a more intuitive correspondence with land use. Generally, if the temperature estimate based on remote sensing data is accurate, the accuracy should be between 1 and 2°C (Vázquez et al., 1997). Research suggests that site density is positively correlated with model accuracy. In other words, the denser the sites, the higher the accuracy of the model (Shen et al., 2020). If Ta is estimated jointly in Northwest China, where meteorological stations are scarce, and in areas where other stations are d ense, the former regions often do not obtain more accurate Ta data (Chen et al., 2015;Li and Zha, 2018). Therefore, it is very necessary to separately model and estimate the area where the data is lacking.
The main objective of this study was to propose a statistical method based on Landsat 8 and auxiliary data for accurately estimating Ta, especially in arid northwest China, where meteorological stations are scarce and unevenly distributed. The specific objectives of this study were to 1) use the RTE method to estimate Ts based on the Landsat8 data, and improve accuracy. 2) Select reasonable independent variables, participate in the modeling of MLR and RF, and compare the results of Ta estimation. 3) Evaluate the performance of the optimal Ta estimation model in time and space scale.

Study Region and Meteorological Station
The study area is the Shiyang River basin in the arid region of Northwest China, which is located in the Hexi Corridor of Gansu Province and the coordinate range is between 101°40′E-104°20′E and 36°30′N-39°30′ N ( Figure 1). It covers an area of 41,600 km 2 , and the range of altitude is 1,157m-5012 m. This area is of a continental temperate arid climate, with the characteristics of aridity where precipitation is 300 mm per year and average annual pan evaporation is 2,000 mm.  Meteorological stations in Northwest China are sparse, with only five stations in and near the study area (Minqin, Wuwei, Yongchang, Wushaoling, Jintai station), which can be collected from the China Meteorological Science Data Center (http://data. cma.cn). The observation station (National Field Scientific Observation and Research Station on Efficient Water Use of Oasis Agriculture in Wuwei of Gansu Province) has a meteorological station to collect Ta data, SI-111 thermal infrared radiometers to obtain Ts data (the coordinates in Table 1). Google Earth Pro was used to verify the location of the meteorological station and ensure that Ta data corresponds totally to the site location. Then obtain the available satellite images from 2014 to 2019 and calculate and extract the Ts data corresponding to the location of the meteorological station through the GEE platform.

Satellite Data and Processing
Google Earth Engine (GEE) is a geospatial processing platform based on cloud computing developed by Google, which promotes fast analysis by using Google's computing infrastructure and providing a convenient platform for applications based on linking to the cloud computing engine (Becker et al., 2021). Most of the algorithms built into the GEE cloud computing platform use pixel-by-pixel calculation functions, so no matter what the area or proportion of the calculation and analysis is required, as long as the research area has available data. The platform is suitable for scientific researchers with a background in non-professional programming and can quickly realize globalscale remote sensing data processing and mining. This research uses the online JavaScript API of the GEE platform (https:// earthengine.google.com/) to access and analyze the data sets used from the public catalog, without downloading images, only outputting the processing results, which improves computing efficiency.
This study uses the Landsat 8 Raw and SR (Surface Reflectance product) dataset in the GEE platform, screening the images with cloudiness not exceeding 30% from 2014 to 2019 and perform cloud removal processing. To cover the entire area of the Shiyang River basin using images of four Landsat-8 tiles 131/033,131/034 and 132/033,132/034 ( Figure 1). By using the GEE platform, the remote sensing images were spliced efficiently, clipped, and parameter inversion, also the meteorological data were interpolated.

Land Surface Temperature
At present, there are three main methods for remote sensing to retrieve land surface temperature: atmospheric Radiative Transfer Equation (RTE) method, single-channel algorithm, and splitwindow algorithm. This study uses the Landsat 8 SR data set in the GEE platform, to retrieve the Ts based on the RTE method.
The expression of the radiation transfer equation of the thermal infrared radiation value (L λ ) received by the satellite sensor is Windahl and Beurs, 2016): where L λ is the spectral radiance value at the top of the atmosphere at the band λ (W · m −2 · μm −1 · sr −1 ); ε is the surface specific emissivity; B (T s ) is the blackbody thermal emissivity brightness (W · m −2 · μm −1 · sr −1 ); τ is the atmospheric thermal infrared band transmittance; L↓ is the downward radiance of the atmosphere after reflection on the ground (W · m −2 · μm −1 · sr −1 ); L↑ is the upward radiance of the atmosphere (W · m −2 · μm −1 · sr −1 ). Knowledge of land surface emissivity (LSE) is necessary to apply the above methods to a Landsat image. Considering different situations, obtain the emissivity value from NDVI (Sobrino et al., 2004;Orhan and Yakar, 2016): where NDVI Veg 0.7 and NDVI Soil 0.05 (Ma and Pu, 2020). When using Landsat8 images, take the 10th band to provide the thermal infrared radiance value. The calculation formula of B(T s ) is: The calculation of the surface temperature (T s ) uses the Planck formula: where Ts is surface temperature (°C); K 1 , K 2 can be obtained from the header file of the remote sensing data. For Landsat8 TIRS Band10, K 1 774.89W · m −2 · μm −1 · sr −1 , K 2 1321.08K. It can be seen that the use of Radiative Transfer Equation method to retrieve the Ts needs to have the atmospheric profile parameters, which can be obtained by entering the shadowing time, latitude, and longitude in the website provided by NASA (http://atmcorr.gsfc.nasa.gov/).

Auxiliary Data and Processing Flow
The near-surface air temperature (Ta) has a good correlation with the surface temperature (Ts) (Benali et al., 2012;Ruiz-Álvarez et al., 2019). Ts is a physical quantity that reflects the degree of cold and heat on the surface of a ground object because it is affected by the characteristics of the underlying surface, such as vegetation coverage and dry and wet conditions. Ta is a physical quantity reflecting the degree of cold and hot air in the atmosphere. The atmosphere has strong fluidity, which is easily affected by the surrounding environment (Xu et al., 2012;Gholamnia et al., 2017). Therefore, when looking for the correlation between Ts and Ta, the influence of various factors such as ground characteristic and environment must be considered. For the estimation of Ta, Jang et al. (2004) showed that Julian day is a more significant parameter than altitude or the solar zenith angle. In addition, we have chosen variables that are often selected as predictors in air temperature modeling literature, such as surface albedo, elevation (DEM), normalized vegetation index (NDVI), and wind speed (Riddering and Queen, 2006;Cristóbal et al., 2008;Hou et al., 2013). The full name of the DEM (WGS84/EGM96) data source is "NASA SRTM Digital Elevation 30 m," which is used directly on the GEE platform and provided by NASA JPL (Farr et al., 2007).
The broadband albedo (α) of the ground surface is a critical variable for many scientific applications, which is the ratio of the total radiant flux reflected by the ground surface to the incident flux (Liang et al., 2003). The calculation formula of α is: Where α is the surface reflectance, and its value is between 0-1.0; α 2 , α 4 , α 5 , α 6 , and α 7 are the 2, 4, 5, 6, and 7 bands of Landsat 8 surface reflectance products.
As has been commented in the previous section, the LSE can be retrieved from NDVI values. The data can be used to construct NDVI according to the following equation: Where NIR and R are the reflection values of the near-infrared and infrared bands, respectively, which are the fifth and fourth bands of Landsat 8. This study mainly considers the relationship between Ta and Ts, NDVI, surface albedo, DEM, wind speed and Julian day. The flow chart of the models used is shown in Figure 2, and the process can be summarized as Data preparation, processing and prediction. Data collection is the basis for establishing a reliable model, but it is also necessary to filter the input variables and try to use as few variables as possible under the condition of high model accuracy. Variables used in this study are readily available, which are closely related to the changes of Ta. Previous studies on Ta estimation lack the verification of Ts and analysis of related results. Therefore, this research explored whether it is necessary to improve the precision of Ts to achieve better results in actual application. After verification and comparison, the best model is selected from the combined methods to estimate Ta.
The spatial and temporal matching among the variables is the main issue for the reliability of the regression implementation. Table 1 shows the sources and resolutions of all data sets. On the spatial scale, the resolutions of spatial variability independent variables such as Ts, Albedo, NDVI, elevation, wind speed (resampled after interpolation) are all FIGURE 2 | The flow chart for estimating daily mean air temperature (Ta). Input variables include surface temperature (Ts), normalized vegetation index (NDVI), surface albedo (Albedo), elevation (DEM), wind speed, and Julian day, and build MLR and RF models. Variables are selected by indicators such as variance inflation factor (VIF), residual sum of squares (RSS), coefficient of determination ( R 2 ), Mallows Cp (Cp) value, and Bayesian Information Criterion (BIC). Finally, compare the simulation results of the four combined models with the measured values, then obtaining the Ta distribution in the study area.
Frontiers in Environmental Science | www.frontiersin.org December 2021 | Volume 9 | Article 791336 30 m × 30 m, which is spatially consistent. On the time scale, Ts, Albedo, and NDVI are instantaneous data, and only wind speed is daily data.

Multiple Linear Regression
The regression model aims to establish the corresponding functional relationship between several independent input parameters and output targets (Giacomino et al., 2011;Agha and Alnahhal, 2012;Williams and Ojuri, 2021). MLR is a linear regression technique and very useful for the best relationship between predictor variables and several independent variables, which is different from a simple linear regression analysis (Akan et al., 2015). R is an excellent tool for statistical calculation and statistical mapping. It is free, open-source software that does not require any license and is simple to operate (Williams and Ojuri, 2021). Using the "lm" function in R software, an MLR model of Ta and multiple correlation factors are established Expanding the MLR equation into the more commonly used form is: where A is the regression target variable; b 0 ∼b n are undetermined coefficients; X 1 ∼X n are independent variables.

Random Forest
The random forest (RF) is a non-parametric machine learning algorithm, which is more flexible than classical statistical models (Genuer et al., 2017;Li et al., 2019). It hardly requires statistical assumptions and is more tolerant of missing values and outliers. Due to the Law of Large Numbers, RF does not overfit. Injecting the correct randomness makes them accurate classifiers and regressors. RF algorithm can automatically distinguish the importance of each variable, give out the dependence between variables, and easily give explanations in combination with professional knowledge (Breiman, 2001). The RF method has been widely used for classification and regression in remote sensing applications (Ke et al., 2016;Park et al., 2016Park et al., , 2018Richardson et al., 2017). RF has begun to be used for Ta estimation in recent years (Ho et al., 2014;Zhang et al., 2016;Yoo et al., 2018). Use the default model parameter settings of R and its contribution packages to develop and apply statistical models (R Core Development Team, 2008;Ho et al., 2014;Liaw and Wiener, 2015).

Variable Adjustment
The absence of complete collinearity between any two independent variables is one of the assumptions of multiple linear regression. Variance Inflation Factor (VIF) is an index used to judge whether there is collinearity. If there is no linear relationship between the independent variables, the VIF value is 1, and a deviation from 1 indicates a trend of collinearity. From the effect of the multicollinearity test, multicollinearity can be tolerated when VIF <10. The VIF value of a variable greater than 10 indicates that there may be estimation problems. If there is multicollinearity, we would simply delete the variable directly, or use a biased estimate for processing (Shabani and Norouzi, 2015;Williams and Ojuri, 2021). The Ta peaks around the 200th day of each year, which is nonlinear with the increase of Julian Day, so assuming that there is a quadratic function relationship. Linearize the Julian Day and adjust it to x7 (J − 200) 2 as the independent variable. The Ts and the Julian day are independent variables with a strong correlation. Therefore, adding the relationship of x8 x7 × Ts as an interaction term to match the model, assuming that the slope of Ts depends on the value of Julian day.
Selecting appropriate variables can not only avoid overfitting but also increase the explanatory degree of the model. The idea of the optimal subset selection method is to model all the variable combinations, then select the model with the best result. The advantage of this method is that all possible combinations are tested, and the final choice must be the best result. However, as the number of candidate variables increases, the amount of calculation will increase exponentially. Therefore, this method is only suitable for situations with few independent variables.
Residual Sum of Squares (RSS), adjusted coefficients of determination (adjusted R 2 ), Mallows's Cp (Cp) value, and Bayesian Information Criterion (BIC) value are used to evaluate model statistics (Cristóbal et al., 2008). The closer Adjusted R 2 is to 1, and the other indicators are smaller, the better the model fits. The optimal model can be determined by comparing the indicators of each variable.

Validation Data and Indicators
The verification of Ts used the infrared sensor (SI-111) observation data in the uniform and widespread farmland from 2014 to 2018. Ta and wind speed data records for 2014-2019 come from the daily data set of surface climate data, obtaining from the China Meteorological Data Service Centre (http://data.cma.cn/). Use Ts and other independent variables from 2014 to 2017 as a training set to build the MLR model and use data from 2018 to 2019 to verify the accuracy. The RF model randomly extracts 70% of all the data as a training set, using the remaining data to validate the resulting model. Such a verification method can explore whether the model constructed by the data from the past years is also applicable to the future years, which achieve the expansion of the time scale.
A set of statistical parameters were calculated to assess the accuracy of the predicted air temperature, including coefficients of determination (R 2 ), root mean square error (RMSE), and model efficiency (MEF). Values of R 2 , RMSE, MEF and can be estimated using the following equations: Frontiers in Environmental Science | www.frontiersin.org December 2021 | Volume 9 | Article 791336 where n is the number of records in validation data sets used in this study, y i is the estimated variable, and O i is the observed variable, y is the mean value of the predicted for all validation sites, O is the mean value of the observed variable for all validation sites.
R 2 always calculated for a significance level of 0.005, was used as a measure of correlation and proportion of observed variability accounted by the model (Benali et al., 2012). RMSE is used to quantify the error (Willmott and Matsuura, 2005). RMSE is an indicator that shows the mean and spatial variance and is used to measure the quadratic error at a single level, which is also particularly sensitive to outliers (Janssen and Heuberger, 1995). The value of MEF is in the range between −1.0 and 1.0. If the performance of the estimation method is poor, the value will be lower (Zheng et al., 2013;Yang et al., 2017;Wang and Lu, 2018). Because MEF integrates correlation and error measurement, it is a robust statistical indicator for model consistency evaluation and reflects the adjustment of the 1:1 line, so it is used to measure the predictive ability of the model (Nash and Sutcliffe, 1970). The above parameters compared the

Surface Temperature Calibration
SI-111 collected Ts data for the 5 years from 2014 to 2018. Based on Landsat8, using RTE to invert the Ts, 40 points corresponding to the position of the ground instrument were extracted. The RTE method estimates Ts from 2014 to 2017, compared with observation data and its R 2 0.77, RMSE 3.27°C. This study found that the RTE method overestimated the value of Ts. Especially when the surface vegetation coverage is low, the error will be more obvious (Figure 3). To resolve this uncertainty, a correction method using NDVI is proposed to improve the accuracy of the RTE method for inversion of Ts.
Correcting RTE with NDVI, the expression is: where NDVI is Normalized Difference Vegetation Index; Ts RTE−NDVI is the surface temperature calculated using the Radiation Transfer Equation corrected by NDVI (°C); Ts RTE is the surface temperature calculated using the original Radiation Transfer Equation (°C). The 31 measured data from 2014 to 2017 were used as the target to revise the RTE ( Figure 4A), and the 2018 data were used as verification ( Figure 4B). After correction using NDVI of  2014-2017, compared with the verified data, the accuracy of the model is R 2 0.83, RMSE 2.09°C. And after the verification of the data in 2018, it is confirmed that Eq. 11 has improved the accuracy of RTE (before: R 2 0.63, RMSE 4.80°C; After: R 2 0.76, RMSE 2.41°C). The fitting line is close to 1:1, indicating that the estimation accuracy of the Ts is significantly improved.

Variable Importance and Selection
Seven variables were tested for multicollinearity, and there was no variable with VIF greater than 10, indicating that the regression model can be established. As shown in Figure 5, RSS monotonically decreases with the increase of the number of independent variables, which cannot be directly judged. Adjusted R 2 and Cp values tended to be the largest when the number of variables was 6 and 7. The adjusted R 2 for Ta models ranged from 0.87 to 0.95, which eventually fixed at 0.95. Further addition of the independent variable did not improve the adjusted R 2 , indicating that 6 independent variables were considered optimal. The change of BIC with the number of variables also proves that the model with six variables has the highest accuracy.
The statistical analysis indicated that Ts, Albedo, NDVI, DEM, x7, and x8 were significant independent variables for estimating Ta. As shown in Figure 6, R 2 was increased with the increase of the independent variables, so optimal combination of independent variables was not proposed in this study. However, in the absence of the wind speed, the adjusted R 2 reached a maximum while Cp and BIC arrived at a minimum, indicating that wind speed does not improve the accuracy of the Ta model. Except that the wind speed does not affect the model construction, the addition of other variables can improve the fitting accuracy of the model. The date of Ta estimation is almost always sunny when the average wind speed is slightly different in space. Therefore, there is no significant correlation between wind speed and Ta variation. Moreover, the average wind speed data is interpolated from the data obtained from meteorological stations, which may be inaccurate at the regional scale. To sum up, it is sufficient to use remote sensing data, elevation, and Julian days as the independent variables for the model.
After determining the independent variables, the formula of the constructed multiple linear regression equation is as follows: (12) where Ta is the near-surface air temperature; Ts is the surface temperature; a is the surface albedo; NDVI is Normalized Difference Vegetation Index; H is the altitude;  The random forest algorithm cannot provide an estimation form similar to MLR because it cannot be parameterized. Two parameters are critical to the operation of the random forest model: the number of trees to grow (ntree) and the number of variables to be used at each node (mtry). These two parameters are selected based on the percentage of model interpretation and mean of squared residuals. The mtry value was tested from 1 to 8, and the ntree value was tested using the following 6 values: 100, 300, 500, 1,000, 1,500, and 2,000. Figure 7 shows the changes in Mean of squared residuals and Percent variance explained after running a random forest with different ntree and mtry values. The highest Percent variance explained and the lowest Mean of squared residuals were obtained when mtry 3 and ntree 1,500 or 2,000. Since the higher the value of ntree, the more cost and time for calculation, so the value of ntree is set to 1,500, and the value of mtry is set to 3 to run the random forest model.
The amount of Ta-related data from 2014 to 2019 is 397. Randomly select 70% of the total amount as the training set to determine the variables and parameters, and the remaining data as the testing set. Figure 8 is a scatter plot of the simulated and observed values of the two models constructed from the training set data. The R 2 and RMSE of the MLR model were 0.953 and 1.74°C, respectively. Most of the points of the RF are concentrated and close to the 1:1 line. The comparison between the results of the RF model and the MLR model shows that the random forest has higher estimation accuracy.

Model Validation
To understand the universality of the built model, the MLR model was verified using the 2018-2019 dataset; the RF algorithm was verified using the reserved test dataset. At the same time, it is FIGURE 9 | The Ts involved in the Ta estimation were calculated through the radiation transfer equation (RTE, (A)) and the NDVI correction (NDVI-RTE, (B)). The coefficient of determination (R 2 ) and root mean square error (RMSE°C) were used to evaluate the accuracy of estimation Ta on the validation set for the two models and the two Ts data.
FIGURE 10 | Relations between the temporal distribution of observed near-surface air temperature (Ta°C) and model performance, using root mean square error (RMSE°C) to represent the accuracy in the estimation of air temperature (Ta).
Frontiers in Environmental Science | www.frontiersin.org December 2021 | Volume 9 | Article 791336 explored whether the improvement of the estimation accuracy of Ts will affect the estimation accuracy of Ta. The prediction result of Ta estimated through the validation data set (Figure 9) cannot reach the accuracy of the initial training set. When Ta is estimated from the validation data set of the RF, the prediction could not achieve the accuracy of the initial training set, respectively. When using NDVI to correct RTE, the R 2 and RMSE of the RF model was 0.943 and 2.12°C, respectively, indicating that the estimation accuracy of using RTE-NDVI is reduced. In contrast, the simulation results of MLR on Ta are relatively stable. The two calculation methods of Ts have little effect on MLR. Therefore, the accuracy of Ts inversion will not affect the accuracy of the two methods for estimating Ta. The estimation results of all models show good performance, and the fitting line is close to 1:1.

Spatial and Temporal Performance
Through the analysis of the average daily Ta and error estimated by MLR for 2018-2019, it was found that the temperature and error on the time scale have spatial variability. As shown in Figure 10, Ta had a clear trend with Julian day and season, and the estimation accuracy had an apparent seasonal trend. Sometimes, the meteorological station is covered by clouds and cannot participate in the calculation or verification, especially at the high-altitude Wushaoling station. The amount of remote sensing image data in winter is insufficient. Although the error is small, it cannot be fully verified. There was not much difference in the accuracy between spring and summer. The RMSE values in spring and summer were both 1.66°C, which had no difference in the estimation accuracy. With the increase of temperature, the error value was relatively stable and had no apparent change trend. In autumn, the RMSE was 2.35°C. At this time, the study area is in the rainy season, and the estimate of Ta will be affected by the increase in rainfall.
In 2018 and 2019, the MEF value of each station was between 0.84 and 0.97, and the overall MEF was 0.946, indicating that MLR has good performance. Analyze the error distribution of Ta estimated by the MLR model at each site. Figure 9 shows that Ta is more accurate at locations far from the mountainous area.
The MEF of the Minqin station for 2 years was close to 1, and the degree of the fitting was relatively high. The MLR model behaved differently in different altitude ranges. As the altitude increases, the performance of the MLR model gradually weakens, which may be caused by the complexity of the terrain. The RMSE has the maximum value at the highest altitude site (Wu Shaolin) ( Figures 11C,D). The average RMSE of Minqin and Jintai Station are below 1.5°C. The accuracy of these two stations is better than that of other stations, possibly because the terrain and environmental conditions of the stations are not complicated. The locations of these two stations are far away from the Qilian Mountains, which are as high as 5,000 m above sea level. The estimation model had the highest relevance and accuracy at the Minqin station, which is located on the north side of the basin, and the surrounding terrain is flat and less affected by high mountains. Although Jintai Station is close to the east of Qilian Mountain, which is far from the mountain range, the estimation accuracy of Ta is also very high.
Because the data of the two adjacent columns of the Landsat satellite are of different periods, remote sensing images of the entire study area cannot be obtained on the same day, so the basin can only be divided into two parts for comparison. Chosen four images with fewer clouds, then calculated the difference between the spatial interpolation result and the MLR model estimate. Inverse Distance Weighted (IDW) assumes that the influence of variables on the surrounding area decreases as the distance to the sample location increases, which was used to interpolate Ta data from the meteorological station. The error of Ta distribution obtained by IDW in the area is related to the underlying surface. Using meteorological and remote sensing data from July 31, August 22, 2018, May 21, May 28, 2019, Ta spatial interpolation and estimation maps were obtained and compared. It iswidely believed that Ta is related to altitude. Temperature decreased with the increase of altitude in most cases. In high altitude areas, where no meteorological station provides temperature data, the Ta of the regions was determined by data from nearby stations. Therefore, Ta in high altitude areas was often overestimated by interpolation (the red area below Figures 12A,C). Simultaneously, there are also large desert areas in the study area (the blue area in Figure 12). Since the temperature in the desert during the day will be much higher than in other zones, the interpolation method can underestimate the temperature by up to 15°C. The Ta of farmland is overestimated during spatial interpolation (the part of the red area above d in Figure 12B), which is also related to the location of meteorological stations. Most of the station is established in cities or border areas. The vegetation coverage of farmland is very high during the growth period of crops, which will affect its Ta. Therefore, when using interpolation methods such as IDW, more meteorological stations are needed to ensure the accuracy of the interpolation results, and at the same time, should pay attention to the unreliability of interpolation methods in mountains or deserts.

Model Selection
In areas with vegetation cover, the temperature value retrieved from remote sensing data is often a mixed value of Ts and canopy temperature (Tc). The accurate estimation of Ts still needs to be analyzed according to the vegetation situation and terrain (Zhang and Li, 2018). Comparing the Ts calculated by RTE with the measured values, it was found that RTE overestimated the Ts, especially when the vegetation coverage was low, consistent with the conclusions published by USGS in 2015 (Barsi et al., 2014;Xu and Huang, 2016). However, due to the characteristics of fewer parameters and suitable for any thermal infrared band, the improved accuracy can also be widely used. Ts is considered to be the most important independent variable the many estimated Ta models (Park, 2011;Zhang et al., 2016;Song and Wu, 2018). In this study, the RMSE of Ts directly calculated by RTE is 3-5°C, and it is 2-3°C after correction. The correction equation is proposed based on the local estimation results, and its universality still needs to be verified. The error is related to the surface emissivity or atmospheric profile parameters Sekertekin and Bonafoni, 2020).
The variables of this study to be selected are closely related to Ta in many studies (Cristóbal et al., 2008;Gholamnia et al., 2017;Shen et al., 2020). Although the correlation between time parameters and Ta was not strong, the adjusted x7 had a significant influence on model accuracy indicating that the adjustment using the Julian Day was useful. Studies have suggested that wind speed is significant, especially when using energy balance methods, which require wind speed to calculate aerodynamic resistance (Hou et al., 2013;Zhang et al., 2015). Through indicators comparison, it is found that when using a statistical method to build the model, the average wind speed is not helpful to improve Ta estimation accuracy (Stisen et al., 2007).
The selected variables participate in MLR and RF to estimate Ta, using the available satellite dataset. The RMSE of the two methods were both lower than 2.0°C, which indicated that they were both suitable for northwest China. Previously, Ho et al. (2014) Using the RF algorithm to estimate the maximum Ta in Vancouver, Canada, and the RMSE was 2.3°C. The optimal MLR model established by Yang et al. (2017) estimates the average Ta and the RMSE is 3.6°C. Therefore, the methods of estimating Ta in this study had relatively high accuracy. However, due to the 16-days revisit cycle of the Landsat satellite and the low temporal resolution of the images, there are limitations in monitoring daily Ta (Yoo et al., 2018).
When performing regression, the random forest cannot make predictions that exceed the range of the training dataset, which may lead to overfitting when modeling some data with specific noise. For many statistical modelers, the random forest is like a black box (Zheng et al., 2019), almost unable to control the internal operation of the model, can only try between different parameters and random allocation. The simulation results of the random forest algorithm for the existing data are very good, but when it is used in the prediction and estimation, its accuracy may drop suddenly. It is difficult to improve it because it cannot control the internal operation of the model. For this study, multiple linear regression is still recommended. The simulation accuracy of the MLR model is similar to that of validation, which is conducive to the evaluation of subsequent prediction results. From the comparison results of the models, the MLR and RF models that without improving the accuracy of Ts performed better, and the actual operation was simpler. The reason for this phenomenon may be that the NDVI involved in the Ts correction is also one of the independent variables of the Ta estimation model, which makes the internal adjustment of the model during operation can ignore the errors of Ts.

Spatial and Temporal Uncertainties of Model
The seasonal variation of the average temperature and the distribution of the station will affect the estimation accuracy of Ta (Holden et al., 2011;Chen et al., 2015). From the time point of view, the RMSE in autumn reaches 2.24°C, which is 0.57°C higher than that in spring and summer. This is consistent with the results of Golkar et al. (2018), which believe that the estimation of Ta in spring and summer has higher accuracy, while there are more uncertainties in autumn. Therefore, the estimation model of daily average temperature is more suitable for spring and summer days. Yang et al. (2017) integrated various statistical indicators and believed that each model performed better in spring, and the estimation accuracy decreased due to the influence of rainfall and cloudy weather. The study area is located in the inland of Northwest China, and rainfall is mainly concentrated in autumn. Benali et al. (2012) believe that the cloud cover of remote sensing images is inversely proportional to the model performance, and higher cloud cover harms model performance. That is, because cloud cover could reduce the accuracy of the Ts, which affects the accuracy of estimated values of daily average temperature. Therefore, the model can be optimized to the season or month scale, and the process of climate and environmental impacts can be added to adapt to the impact of climate change and cloud cover. The spatial distribution map of RMSE shows that the model in mountainous and plateau areas generally performs worse than in plain areas. This conclusion is consistent with the research results of Yang et al. (2017). Due to the influence of terrain differences, mountainous terrain has a broader spatial variation range than flat terrain, so mountain temperature changes are more complicated. On relatively flat terrain and vegetation, horizontal uniformity leads to stable atmospheric conditions (Blandford et al., 2008;Lin et al., 2016). Ta raster data is usually limited by site coverage, especially in mountainous areas where the density and elevation distribution of meteorological observations vary greatly (Holden et al., 2011). Shen et al. (2020) believe that in the case of large-scale use of MLR and RF, the accuracy of the estimation results is poor due to the scarcity of sites in Northwest China that can participate in model training. When MLR is estimated in the Shiyang River Basin, the MEF of the highest altitude station (3000 m) is higher than 0.83, and the RMSE 2.10°C, the model also performs well.
The method proposed by the research can obtain highprecision Ta estimation results when using the Landsat 8 data set, and improve the spatial resolution. However, the actual situation of only six meteorological stations limits the accuracy of Ta inversion in this area. We hope that more reliable verification data can be obtained in future studies.

CONCLUSION
The purpose of this study is to estimate Ts and Ta in the Shiyang River Basin. Perform parameter inversion based on Landsat8 images, use RTE to calculate and correct Ts. Using Ts, NDVI, elevation, and other parameters as variables to construct MLR and RF models for estimating Ta. The results show that: For MLR and RF, after calibrating Ts and participating in the estimation of Ta, both methods have accurate estimation results; The accuracy of RF training results is better than MLR, but the test set results of the two models are not significantly different; Although the topography of the study area is complex and the land cover conditions are different, the MLR model applies to the study area. In addition, compared with the MLR model, the comparison found that the interpolation method will underestimate the temperature in areas with low vegetation coverage like deserts, while the opposite is in the mountains and farmland areas.
This study can be used to accurately understand the distribution characteristics and changing trends of Ts and Ta in the study area. Further research can optimize the model in time to a month to adapt the climate change; make up the deficiency of Landsat8 and MODIS data by fusing remote sensing data, then obtaining Ta data with high temporal and spatial resolution.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.