Point-to-Surface Upscaling Algorithms for Snow Depth Ground Observations

: To validate the accuracy of snow depth products retrieved from passive microwave remote sensing data with a high conﬁdence level, the veriﬁcation method based on points of ground observation is subject to great uncertainty, due to the scale effect. Thus, it is necessary to use a point-to-surface scale transformation method to obtain the relative ground truth at the remote sensing pixel scale. In this study, by using the snow depth ground observations at different observation scales, the upscaling methods are conducted based on simple average (SA), geostatistical, Bayes maximum entropy (BME), and random forest (RF) algorithms. In addition, the cross-validation of the leave-one-out method is employed to validate the upscaling results. The results show that the SA algorithm is seriously inadequate for estimating snow depth variation in space, and is only suitable for regions with relatively ﬂat terrain and small variation of snow depth. The BME algorithm can introduce prior knowledge and perform kernel smoothing on observed data, and the upscaling result is superior to geostatistical and RF algorithms, especially when the observed data is insufﬁcient, and outliers appear. The results of the study are expected to provide a reference for developing a point-to-surface upscaling method based on snow depth ground observations, and to further solve the uncertainties caused by scale effects in snow depth and other land surface parameter inversion and validation, by using remote sensing data.


Introduction
The scale effect is one of the important research features in geography and remote sensing science [1,2]. Due to the spatial heterogeneity of the elements of the Earth's surface, the observed phenomena and summarized laws at one scale may be effective, similar, or need to be corrected at another scale [3,4]. Snow is one of the major elements of the cryosphere, and the spatial heterogeneity of snow is particularly strong, especially in mountainous areas. The snow depth varies greatly in near space. There are large errors and uncertainties in the construction and verification of the snow depth inversion model from passive microwave remote sensing data based on point-based snow depth ground observations. The main reason is that the ground observation points do not match the spatial scale of remote sensing pixels, and the observation points have difficulty representing the ground truth of pixel scale. Therefore, improving the accuracy and authenticity of remote sensing quantitative inversion of snow parameters by developing an upscaling method based on ground observation points, is of great significance.
The terrain is the main factor affecting the spatial differentiation of snow depth, including slope, aspect, and elevation, which determine the spatial distribution of snow depth in mountainous areas [5]. Terrain fluctuation, windblown snow, and other factors also lead to large differences in the spatial distribution of snow depth, which directly affect the accuracy of snow depth inversion by passive microwave remote sensing technology [6].
Passive microwave is the main data source for snow depth inversion because of its strong penetration ability [7]. However, due to the coarse resolution of the passive microwave remote sensing data, the problem of mixed pixels caused by the high spatial heterogeneity of snow cover is the main reason for the great uncertainties of passive microwave snow depth inversion [8]. Due to the scale effect, it is usually difficult to obtain reasonable and objective evaluation results when directly using ground 'point' measurement data instead of 'surface' to verify the accuracy of remote sensing products in areas with a large spatial heterogeneity of snow depth. Although increasing the number of observation points can compensate for the uncertainty caused by heterogeneity to some extent [9], the problem of spatial scale mismatch still exists. Therefore, the ground observation 'point' data upscaling to 'surface' into the pixel scale is an operational way to improve the accuracy of snow depth inversion as well as the authenticity of verification.
The point-to-surface upscaling transformation method can be divided into the simple average method, empirical regression method, geostatistical method, Bayesian method, etc. [10]. Based on the geography of Greenland, the simple averaging method was used to upscale weather station data and compare it with Climate-SAF surface albedo products; it concluded that the result was consistent with the in-situ observation, but there were problems in areas where the terrain and environment changed dramatically [11]. The Kriging method, which is a geostatistical method, is the most widely used in point-tosurface upscaling [12][13][14][15]. The conclusion indicated that the co-kriging method is usually superior to the ordinary kriging method, and that the greater the influence of the covariates on the target parameters, the better the co-kriging results [15]. It has also been shown that the reduced major axis (RMA) is the most practical upscaling method for modeling regional net primary productivity by comparing RMA and Kriging methods [16]. In addition, the machine learning algorithm has also been widely used in snow parameter inversion [17,18], but the research on machine learning in scale transformation remains to be further expanded.
It is important to obtain the true values of snow depth ground observations at the pixel scale for the construction and verification of snow depth remote sensing inversion models. In this paper, we compare the point-to-surface upscaling methods based on simple averaging, geostatistics, Bayesian maximum entropy, and random forest with the observed snow depth plots from three major snow areas in China, Northern Xinjiang, the Tibetan Plateau, and the northeast. The study is expected to provide a reference for developing a point-to-surface upscaling scheme for snow depth ground observation into the pixel scale, and further, to solve the uncertainties in the construction and evaluation of the snow depth inversion model based on passive microwave remote sensing data.

Ground Observation Data of Snow Depth
The snow depth ground plots dataset is derived from the National Cryosphere Desert Data Center (http://www.ncdc.ac.cn, accessed on 15 February 2022), which includes 25 km 2 and 500 m 2 snow plot datasets in typical snow areas in China [19,20]. The dataset contains 31 snow plots and 1440 snow observation points in total. There are 15 plots of 25 km 2 and 16 plots of 500 m 2 ( Table 1). The observation parameters include snow depth, air temperature, snow pressure, snow density, altitude, latitude, and longitude. Each 25 km 2 plot contains 20-82 observation points, with an average of 43 observation points. Each 500 m 2 snow plot contains 13 to 66 observation points, with an average of 41 observation points.

Digital Elevation Model (DEM)
The Shuttle Radar Topography Mission (SRTM), launched by NASA on 11 February 2000, acquired approximately 9.8 trillion bytes of radar image data from the 'Endeavour' spacecraft carrying the SRTM system. After more than two years of data processing, SRTM 1 and SRTM 3, two categories of DEM, were produced, corresponding to 30 m and 90 m spatial resolution, respectively [21]. In this paper, 90 m resolution SRTM 3 DEM data are used as auxiliary data in the point-to-surface upscaling process.

Method
The core element of point-to-surface upscaling is to convert the data of ground observation points to the satellite pixel scale through the functional relationship of Formula (1), so that the scale between ground observation and remote sensing pixel can be matched.
where Z is the pixel scale value, that is, the upscaling result; Z(x i ) is the observation data; x i is the ith field observation point.

Simple Average
The simple average (SA) method is an upscaling method that directly averages the true values of all observed data for each snow plot within a pixel scale. This process is conducted in the R language using the 'base' package.
where Z is the true value of the pixel scale; Z(x i ) is the observed value of field observation point i; n is the total number of field observation points.
The SA method has high requirements for spatial heterogeneity and spatial distribution of observation points. The results of upscaling are representative in the surface with less vegetation, and are relatively homogeneous, while large errors usually occur in the areas with complex terrain.

Geostatistical Method
The geostatistical method is based on the theory of regionalized variables, and uses the variogram function as the basic tool, which can model natural phenomena that have a certain randomness and a certain structure in spatial distribution [22][23][24]. The kriging method uses the measured values for weighting to estimate the predicted values of the unmeasured locations, and the general formula of the kriging method can be expressed as: where Z(x 0 ) is the value of the predicted position, λ i is the weight at position i, and Z(x i ) is the measured value at position i.
In this study, Ordinary Kriging (OK) and Simple Kriging (SK) functions provided by ArcGIS (Version 10.5) software are used to upscale the ground observation data of snow depth. In the OK process, the 'Measure Error' is set at 100% to ensure the nugget constant is entirely composed of measurement error, and there is no random variation caused by the microstructure of the variables.

Bayesian Maximum Entropy
Bayesian Maximum Entropy (BME) was first proposed by Christakos in 1990. It is a method to predict the distribution of time and space based on the Bayesian framework, combined with the concept of entropy in information theory [25,26]. The basic process of BME mainly includes the prior stage, intermediate stage, and posterior stage. In the prior stage, the prior probability density function is calculated by the generalized knowledge base and the maximum entropy principle. The expression is: represents the prior probability density function of random variables x map based on the generalized knowledge base (KG); x map represents random variables; C is the constant that plays the role of normalization; N C is the number of conditions; µ α is the Lagrange multiplier; g α is used for the known function relation of random variables.
The posterior probability function equation of the point to be estimated is: where f K (X 0 ) is the posterior probability density function of the point to be estimated; X so f t is the interval from soft data, and the value range is [α, β].
If the soft data is given in the form of the probability density function, the posterior probability function equation of the estimated point is: where f K (X 0 ) is the posterior probability density function of the point to be estimated; f S X so f t is soft data in the form of the probability density function.
In the prior stage, the probability density function was calculated based on latitude, elevation, slope, and aspect, for the snow depth at each observation point. The soft data were then set to a normal distribution with "Gaussian" type in the intermediate stage.
Combining with the "Gaussian" and "Holecos" covariance, the "nest number" was set to 2 or 3 in the posterior stage. Finally, the theoretical model was automatically fitted literately until the optimal upscaling model was fitted. The process was conducted in the STAR-BME package in the QGIS (version 2.18.28) software.

Random Forest
Random Forest (RF) is a machine learning algorithm based on a decision tree and bagging (bootstrap aggregation) proposed by Breiman [27]. Among the trees, the decision tree is divided into classification decision trees and regression decision trees, namely CART (Classification and Regression Trees) [28]. The RF algorithm in this paper combined random feature selection and the bagging method with CART as a weak learner to upscale for each snow plot. The snow depth of the field observation point was selected as the dependent Remote Sens. 2022, 14, 4840 5 of 11 variable of the model, and the longitude, latitude, elevation, slope, aspect, and observation time of the point selected as the covariates. The RF model predictions can be expressed as: where n is the number of covariates and y i (x 0 ) is the covariate at position x 0 (i = 1, 2, . . . , n). The package 'randomForest' in R language was used for building RF upscaling model in each snow plot. The number of trees "ntree" was started from 1000 to 2500, and the number of variables 'mtry' were set from 1 to 6.

Validation
The Leave-One-Out Cross Validation (LOOCV) and Mean Absolute Error (MAE) are used as indicators to test the accuracy of the above upscaling methods. The root means square error (RMSE) of prediction results was calculated to evaluate the upscaling method.
and N is the number of observation points.
The MAE describes the error between the observed value and the predicted value: the formula is and N is the number of observation points.

Upscaling Results Comparison
The snow depth observed points within the 25 km 2 plots were upscaled to 500 m and 25 km resolution, respectively, and the 500 m 2 plots were upscaled to 500 m resolution. Taking a 500 m resolution upscaling result of a 25 km 2 snow plot as an example, the upscaling results of the four methods except the SA method are shown in Figure 1. The results of OK, SK, and BME are consistent in space, and seriously affected by the spatial distribution of observed snow depth points. The maximum snow depth area and the minimum snow area are near the corresponding maximum and minimum snow depth observation points, showing an obvious 'island' characteristic. The upscaling results of RF are consistent with the other three methods in the high and low snow depth areas. Due to the introduction of elevation, slope, aspect, and other covariates, the influence of the spatial distribution of observation points on snow depth is reduced, and the lack of observation points may also reflect the characteristics of snow depth change in RF. Remote Sens. 2022, 14, x FOR PEER REVIEW 6 of 11  Table 2 shows the MAE of different upscaling methods. The results show that the MAE of BME is the smallest, and the MAE after 25 km 2 and 500 m 2 upscaling is 0.32 cm and 0.97 cm, respectively, with an average of 0.65 cm. The SA method has the largest MAE, with an average MAE of 2.67 cm. OK is slightly better than SK, and the MAE of RF is slightly worse than the two Kriging methods. In addition, the observation scale (plot size) has little effect on the MAE of the point-to-surface upscaling, and the results of different observation scales cannot explain which observation scale can achieve the optimal upscaling results. Table 3 shows the RMSE of different upscaling methods. The statistical results also show that the BME upscaling results are the best, with an average RMSE of 2.95 cm, followed by the SK, OK, RF, and SA methods. The RMSE of the SA method is 3.54 cm, and the error is the largest. The RMSE statistical results also show that the observation scale will affect the accuracy of the upscaling results. The upscaling results of 500 m 2 plots are significantly better than those of 25 km 2 plots. Different upscaling methods exhibit a decrease in RMSE with decreasing observed plot size. Combined with the MAE statistical results, it is shown that the SA method has difficulty obtaining the ideal upscaling results in complex areas with high spatial heterogeneity. In the authenticity test of snow depth remote sensing products, even if more snow depth observations in a remote sensing pixel are obtained, the use of the SA method to obtain the average snow depth observations to directly verify the snow depth products produces greater uncertainty.   Table 2 shows the MAE of different upscaling methods. The results show that the MAE of BME is the smallest, and the MAE after 25 km 2 and 500 m 2 upscaling is 0.32 cm and 0.97 cm, respectively, with an average of 0.65 cm. The SA method has the largest MAE, with an average MAE of 2.67 cm. OK is slightly better than SK, and the MAE of RF is slightly worse than the two Kriging methods. In addition, the observation scale (plot size) has little effect on the MAE of the point-to-surface upscaling, and the results of different observation scales cannot explain which observation scale can achieve the optimal upscaling results. Table 3 shows the RMSE of different upscaling methods. The statistical results also show that the BME upscaling results are the best, with an average RMSE of 2.95 cm, followed by the SK, OK, RF, and SA methods. The RMSE of the SA method is 3.54 cm, and the error is the largest. The RMSE statistical results also show that the observation scale will affect the accuracy of the upscaling results. The upscaling results of 500 m 2 plots are significantly better than those of 25 km 2 plots. Different upscaling methods exhibit a decrease in RMSE with decreasing observed plot size. Combined with the MAE statistical results, it is shown that the SA method has difficulty obtaining the ideal upscaling results in complex areas with high spatial heterogeneity. In the authenticity test of snow depth remote sensing products, even if more snow depth observations in a remote sensing pixel are obtained, the use of the SA method to obtain the average snow depth observations to directly verify the snow depth products produces greater uncertainty. It is considered that the use of the average error (RSME) may be affected by extreme values, and cannot reflect the stability and overall level of the various upscaling methods. Therefore, based on the cross-validation results of different snow plots, the RMSE boxplot is drawn to verify further the results of different point-to-surface upscaling methods (Figure 2). The median in the box diagram represents the overall level of each upscaling method, and the box height (the difference between the 25th percentile and 75th percentile) can indicate whether the RMSE is discrete, which corresponds to the stability of the different methods. The black dots in the box line diagram indicate outliers. In the 500 m 2 snow plots, the overall level gap between the five upscaling methods is small, among which the BME is the best, and the median RMSE is 1.91 cm, followed by SK, OK, SA, and RF. The BME method is the most stable, with a box height of only 1.87 cm, followed by SA, OK, RF and SK. The BME method is the best upscaling method for 500 m 2 observation plots. SK has the smallest RMSE median (RMSE = 3.4 cm) and the best stability (2.83 cm) in the 25 km 2 snow plot, indicating that the SK method is the best for 25 km snow plots. The observation scale has a great influence on the overall level and stability of the upscaling results, and the results show that the overall level and stability of different upscaling results at the 500 m 2 observation scale are better than those at the 25 km 2 observation scale. The 500 m 2 snow plot has a smaller observation extent and a relatively lower spatial heterogeneity, and the observation points are more representative in the upscaling process; therefore, the upscaling results show a better overall level and stability compared with the 25 km 2 plots when the number of observation points is basically the same. It is considered that the use of the average error (RSME) may be affected by extreme values, and cannot reflect the stability and overall level of the various upscaling methods Therefore, based on the cross-validation results of different snow plots, the RMSE boxplo is drawn to verify further the results of different point-to-surface upscaling methods (Fig  ure 2). The median in the box diagram represents the overall level of each upscaling method, and the box height (the difference between the 25th percentile and 75th percen tile) can indicate whether the RMSE is discrete, which corresponds to the stability of the different methods. The black dots in the box line diagram indicate outliers. In the 500 m snow plots, the overall level gap between the five upscaling methods is small, among which the BME is the best, and the median RMSE is 1.91 cm, followed by SK, OK, SA, and RF. The BME method is the most stable, with a box height of only 1.87 cm, followed by SA, OK, RF and SK. The BME method is the best upscaling method for 500 m 2 observation plots. SK has the smallest RMSE median (RMSE = 3.4 cm) and the best stability (2.83 cm in the 25 km 2 snow plot, indicating that the SK method is the best for 25 km snow plots The observation scale has a great influence on the overall level and stability of the upscal ing results, and the results show that the overall level and stability of different upscaling results at the 500 m 2 observation scale are better than those at the 25 km 2 observation scale The 500 m 2 snow plot has a smaller observation extent and a relatively lower spatial het erogeneity, and the observation points are more representative in the upscaling process therefore, the upscaling results show a better overall level and stability compared with the 25 km 2 plots when the number of observation points is basically the same.  Figure 3 compares the RMSE of upscaling results in two different observation scales for each plot. The cross-validation results of 25 km 2 snow plots show that the RMSE dif ferences of different upscaling methods are small, but the differences in 10th and 12th plots are larger. This is because the No.10 and No.12 snow plots were observed twice, a different times for the same location, and the terrain in this area is complex and of large relief. The elevation span of the observation point in the snow area is between 800 and 1273 m. The large elevation difference has a great influence on the different upscaling algorithms. The cross-validation results of 500 m 2 snow plots show that there are large differences of RMSE in No.2, and the RMSE of the BME method is the smallest. The area where the plot is located is bare and the terrain is relatively flat and relatively uniform in space. However, there is an obvious maximum snow depth in the observation point of the plot, and the snow depth is much larger at than other observation points, which leads to the poor results of other upscaling methods. The BME method can smooth the data in the upscaling process, eliminate the influence of the deviation of the sampling data on the  Figure 3 compares the RMSE of upscaling results in two different observation scales for each plot. The cross-validation results of 25 km 2 snow plots show that the RMSE differences of different upscaling methods are small, but the differences in 10th and 12th plots are larger. This is because the No.10 and No.12 snow plots were observed twice, at different times for the same location, and the terrain in this area is complex and of larger relief. The elevation span of the observation point in the snow area is between 800 and 1273 m. The large elevation difference has a great influence on the different upscaling algorithms. The cross-validation results of 500 m 2 snow plots show that there are large differences of RMSE in No.2, and the RMSE of the BME method is the smallest. The area where the plot is located is bare and the terrain is relatively flat and relatively uniform in space. However, there is an obvious maximum snow depth in the observation point of the plot, and the snow depth is much larger at than other observation points, which leads to the poor results of other upscaling methods. The BME method can smooth the data in the upscaling process, eliminate the influence of the deviation of the sampling data on the upscaling results, and achieve the effect of detrending, which is the possible reason for BME performing better in upscaling.

Error Analysis
emote Sens. 2022, 14, x FOR PEER REVIEW 8 upscaling results, and achieve the effect of detrending, which is the possible reaso BME performing better in upscaling.

Sensitivity Analysis
The results indicate that terrain, snow depth, upscaling method, and other fa have a certain impact on the upscaling results, and it is particularly important to sele appropriate upscaling method in a specific environment. The profile curvature also affects the upscaling algorithm. When the profile c ture is equal to 0, it represents the flat terrain, the negative value represents the de and the positive value represents the ascent. When the terrain is flat, the RMSE of dif observation scales is small. With the ascent or descent of the terrain, the RMSE incr This trend is most pronounced in the 25 km 2 plot, but the RMSE of the depression su is greater than that of the uplift surface, because the depression surface can accum more snow and the snow depth is usually thicker than on the uplift surface.
The effects of snow depth, terrain, and observation scale on the different upsc methods were compared. We can conclude that when the snow depth is shallow

Sensitivity Analysis
The results indicate that terrain, snow depth, upscaling method, and other factors have a certain impact on the upscaling results, and it is particularly important to select the appropriate upscaling method in a specific environment. The profile curvature also affects the upscaling algorithm. When the profile curvature is equal to 0, it represents the flat terrain, the negative value represents the descent, and the positive value represents the ascent. When the terrain is flat, the RMSE of different observation scales is small. With the ascent or descent of the terrain, the RMSE increases. This trend is most pronounced in the 25 km 2 plot, but the RMSE of the depression surface is greater than that of the uplift surface, because the depression surface can accumulate more snow and the snow depth is usually thicker than on the uplift surface.
servation scales are small, regardless of the flat, rough surface. In addition, when the ob servation scale is large (25 km 2 ), the difference between the SA method and other method is clear, and the error is large. However, at a small observation scale (500 m 2 ), the erro between different upscaling methods is close. In summary, the Bayesian maximum en tropy method has the best adaptability to areas with large snow depth changes and com plex terrain, especially in the range of large observation scales.   servation scale is large (25 km 2 ), the difference between the SA method and other method is clear, and the error is large. However, at a small observation scale (500 m 2 ), the erro between different upscaling methods is close. In summary, the Bayesian maximum en tropy method has the best adaptability to areas with large snow depth changes and com plex terrain, especially in the range of large observation scales.   The effects of snow depth, terrain, and observation scale on the different upscaling methods were compared. We can conclude that when the snow depth is shallow, the RMSEs of different upscaling methods are small, but with the increase of snow depth the BME algorithm is the best, and the errors are small under different observation scales. Under different observation scales, the errors of different upscaling methods are small at low altitudes. However, with the increase in altitude, the upscaling result of BME is best in the 25 km 2 plots. In the 500 m 2 snow plots, the RMSE of RF is small at high altitudes. The influence of surface curvature on BME is the smallest, and the errors at different observation scales are small, regardless of the flat, rough surface. In addition, when the observation scale is large (25 km 2 ), the difference between the SA method and other methods is clear, and the error is large. However, at a small observation scale (500 m 2 ), the error between different upscaling methods is close. In summary, the Bayesian maximum entropy method has the best adaptability to areas with large snow depth changes and complex terrain, especially in the range of large observation scales.

Conclusions
To obtain the relative true values of ground snow depth observations at the remote sensing pixel scale, this paper uses five methods: SA, SK, OK, BME, and RF, to study the upscaling of ground snow depth observations, and evaluates the accuracy of the five upscaling methods using the leave-one-out cross-validation method, obtaining the following main conclusions:

1.
The SA method is easy to calculate and consumes less time. It can obtain a better upscaling result in the spatial distribution area of snow depth with homogenization, but it is not applicable in areas with large spatial heterogeneity or outliers of snow depth.

2.
Geostatistical and BME methods can obtain better upscaling results than the RF algorithm, and the overall accuracy is better than RF. However, the RF algorithm can produce more consistent results with the spatial distribution characteristics of snow depth in areas lacking observation points under the background of introducing elevation, slope, aspect, and other covariates.

3.
When the amount of observation data is small and outliers appear, different upscaling methods will be greatly affected. The BME algorithm can smooth the data in the process of upscaling, achieve the effect of data detrending, and obtain a relatively better upscaling result. 4.
The observation scale, number of points, and terrain have great influences on the different upscaling methods. Therefore, to obtain more accurate pixel-scale ground observation true values, when the authenticity test of snow depth remote sensing products is carried out, the homogeneous surface with relatively flat terrain should be selected as much as possible, and the number of observation points should be increased, which can effectively improve the accuracy of upscaling results.
So far, passive microwave remote sensing is still the main data source for acquiring snow depth information on a regional and global scale. The traditional development and validation of the passive microwave snow depth inversion model, based on site observation, is subject to great uncertainty, because large errors may occur through the mismatch of observation scales from ground and from space. The point-to-surface upscaling method can obtain ground observation data at pixel scale and solve the uncertainty caused by the scale effect. This study is expected to provide a reference for developing a point-to-surface snow depth upscaling scheme, and to further solve the uncertainty problem caused by scale effects in passive microwave snow depth inversion and validation.  Data Availability Statement: The snow depth ground observed quadrat in two scales, 500 m 2 and 25 km 2 , respectively, are available at http://www.ncdc.ac.cn (accessed on 15 February 2022).

Conflicts of Interest:
The authors declare no conflict of interest.