Moderate-resolution snow depth product retrieval from passive microwave brightness data over Xinjiang using machine learning approach

ABSTRACT

Passive microwave (PM) remote sensing have been extensively used for snow depth (SD) estimation. However, current SD products from traditional PM data fail to capture the differentiation in mountainous and complex terrains with coarse resolution. Therefore, this study incorporates factors such as geographical location, topographic features, and land cover, along with various machine learning algorithms including Gaussian process regression (GPR), support vector machine (SVM), random forest (RF), extreme gradient boosting (XGBoost), and light gradient boosting machine (LightGBM), to construct and optimize SD estimation using enhanced-resolution PM data. The results demonstrate the following: (1)With the auxiliary variables, the SD product from LightGBM-based models exhibits the highest accuracy. (2)The performance of SD products from the LightGBM-based model varies monthly and annually, with shallow snow cover being slightly overestimated (30 cm). (3)The reliable SD product indicates spatial distribution characteristics in Xinjiang, with regions demonstrating no significant improvement being larger than those with no significant degradation. The above results illustrate the remarkable advantages of machine learning in capturing SD distribution and its spatio-temporal variation bolstered by enhanced PM data and auxiliary data.

KEYWORDS:

• Snow depth estimation
• machine learning
• brightness temperature
• enhanced resolution passive microwave data

1. Introduction

Snowpack, a prominent component of the cryosphere, is widely distributed across various regions (Koch et al. 2019; Li et al. 2019). It has unique physical properties, such as high albedo, low thermal conductivity, high emissivity, and sizable latent heat of ablation (Zolles and Born 2021). Snow depth (SD) is essential in regulating the energy between the Earth's surface and atmosphere, hydrological circulation, and carbon cycle (Liu et al. 2021b; Viallon-Galinier, Hagenmuller, and Lafaysse 2020; Xiao et al. 2020). Therefore, it is crucial to estimate SD distribution for a better understanding of regional water resources and climate changes (Adib et al. 2021).

Since the 1970s, remote sensing data has been utilized for SD inversion on a global or regional scale (Liang et al. 2015). However, the nonlinear/linear relation between reflectance and SD is unavailable for optical imagery (Wan et al. 2022). In contrast, SD retrieved by active/passive microwave data are immune to clouds, fog, or other weather conditions (Luojus et al. 2021; Zhu et al. 2021). As for active microwaves, the backscatter with high resolution is a mixture of snow and ground objects, which poses challenges for accurate SD estimation (Awasthi et al. 2021; Dai et al. 2022; Patil, Mohanty, and Singh 2020; Tedesco and Miller 2007). Passive microwave (PM) penetrating a certain snowpack depth can retrieve SD through the volume scatter of snowpack under dry conditions (Adib et al. 2021; Awasthi and Varade 2021). PM sensors provide a wider coverage, with extensive daily observations accumulated over decades (Mashtayeva et al. 2016; Merkouriadi et al. 2021; Wang et al. 2021). However, the spatial resolution of SD products derived from PM data is often too coarse (10–25 km) to accurately study snow hydrology at the watershed scale (Hao et al. 2019; Yang et al. 2019). To overcome the spatial limitations of PM data, NASA has introduced a new type of PM brightness temperature data with a spatial resolution ranging from 3.125 to 6.25 km. While this data have been used to estimate SD distribution in Afghanistan and its surrounding regions (Bair et al. 2018), its potential applicability to other regions remains unexplored.

Regarding PM data-based SD retrieval, various models have been developed to suit the specific frequencies, polarization, and characteristics of different study areas (Awasthi and Varade 2021; Feng et al. 2021). In general, these SD retrieval models can be categorized into two main approaches. The first approach involves a statistical method that relies on analyzing the linear and nonlinear relationship between SD and the gradient of multi-frequency brightness temperature (Che et al. 2008, 2016; Wei et al. 2021). The other approach is based on the utilization of brightness temperature and snow characteristic parameters, employing a radiative transfer model to derive the SD (Gu et al. 2019; Kang, Tan, and Kim 2019; Yang et al. 2021). To address the issue, various algorithms based on linear-described brightness temperature gradients, such as Chang, Foster, and Che, have been developed for SD inversion on both hemispherical and global scales (Chen et al. 2020; Chen, Muthu, and Cb 2021; Dai et al. 2012; Yang et al. 2019; Yu et al. 2012; Zhang et al. 2017). However, the use of fixed empirical coefficients do not conform to spatial–temporal variations of snow characteristics, resulting in errors in SD estimation (Yang et al. 2019). Radiative transfer models, such as HUT and MEMLS, require the incorporation of changes in snow particle size and snow density over time to retrieve the SD. This is achieved by establishing a relationship between brightness temperature and SD (Dai et al. 2017; Wiesmann and Mätzler 1999; Yang et al. 2021). However, owing to the complex mechanisms involved, these algorithms necessitate the utilization of sophisticated auxiliary data for accurate simulations, thereby limiting their global applicability. Therefore, the nonlinear-described brightness temperature gradients is gaining traction in SD estimation for shedding many surface empirical parameters, better accuracy over linear methods, and fewer geographical limitations (Han et al. 2023; Kelly 2009; Saberi et al. 2021). Nonlinear algorithms like artificial neural network (ANN), support vector machine (SVM), and random forest (RF) have been extensively employed for SD estimation (Ntokas et al. 2021; Xiao et al. 2018; Yang et al. 2020). With an extensive range of training samples, these algorithms incorporate the topographic variables of SD distribution, thereby improving estimation accuracy to a certain extent. Recently, gradient-boosting decision tree (GBDT) algorithms, such as RF, light gradient-boosting machines (LightGBM), and eXtreme Gradient Boosting (XGBoost), have emerged as state-of-the-art methods in addressing non-linear relationships within environmental problems. As an advanced ensemble algorithm, these GBDT methods effectively tackle nonlinearity and overfitting of data, while alsoproviding feature selection and ranking capabilities. This makes them especially well-suited for analyzing high-dimensional datasets and achieving superior accuracy. Consequently, there is a pressing need to investigate the application of GBDT algorithms in SD estimation.

Located at the center of Eurasia, Xinjiang is far from the ocean. This region serves as a representative landscape of arid regions in Eurasia and is a crucial component of the Pan-Third Pole Region (Liu et al. 2018). Remarkably, Xinjiang boasts the largest seasonal snow cover in China, accounting for up to one-third of the country’s total snow coverage (Yang et al. 2020). In contrast to maritime snow in Europe, Xinjiang's seasonal snow is dry-cold (Liu et al. 2020; Liu et al. 2021a). Based on SSM/I, SMMR, and AMSR-E PM data, many studies have used linear and nonlinear models (RF, SVM) to retrieve the SD distribution in Xinjiang (Xiao et al. 2018; Yang et al. 2020). Generally, the accuracy of snow depth (SD) estimation is affected by various factors including complex topography, uneven distribution, and fast diurnal variation in snow cover. Insufficient ground observation further hampers the accuracy of SD estimation (Chen et al. 2020; Chen, Muthu, and Cb 2021; Dai et al. 2012; Yu et al. 2012; Zhang et al. 2017). Some SD products, such as GlobSnow, do not provide SD in mountainous areas (Chen et al. 2020; Chen, Muthu, and Cb 2021; Dai et al. 2012; Luojus et al. 2021; Venäläinen et al. 2021; Yu et al. 2012; Zhang et al. 2017). In addition, most Xinjiang SD products are actually derived from global products, primarily operating at a resolution of 25 km. This is due to the coarse spatial resolution of PM bands. Therefore, while reflecting the differentiation characteristics of SD distribution in the mountain-oasis-desert system (MODS) is challenging, moderate-resolution SD could reveal its non-stationary variations.

Given the above considerations, we have developed a comprehensive framework for estimating the spatial distribution of SD. This framework includes resolution-enhanced PM brightness temperature data and auxiliary inputs, such as geographical location, topographic variables, and land cover. To create the most accurate and high-resolution SD product (covering the period from 2011 to 2019, with a spatial resolution of 3.125 km and temporal resolution of 1 d), we have utilized advanced ensemble-based machine learning techniques like LightGBM and XGBoost, and compared them with commonly used algorithms like RF, SVM, and GPR. The Shapley Additive Explanations (SHAP) value illustrates how the models estimate by revealing the contribution of each covariate. The specific objectives of this study are as follows:

1. To retrieve the best SD product in the MODS of Xinjiang, China, according to the performance of machine learning methods (XGBoost, LightGBM, RF, GPR and SVM).

2. To elaborate on the importance of auxiliary variables (geographical location, topographic variables, cover types, etc.) for improving SD production from global and local perspectives.

3. To clarify the temporal and spatial characteristics of SD distribution with a high spatial resolution (3.125 km).

2. Study region and data

2.1. Study area

As the core area of the Silk Road, Xinjiang is located in the northwestern part of China (Elevation: −192−8357 m; Lat: 34−50°N and Lon: 73−97°E) with an area of about 1,660,000 km² (Figure 1). The alpine regions in Xinjiang include the Altai Mountains in the north, the Kunlun Mountains in the south, and the Tianshan Mountains in the center. The Tianshan Mountains divide Xinjiang into two parts: northern and southern Xinjiang, where the Junggar and Tarim Basins are located (Liu et al. 2018).

Figure 1. Geographical location of study area and land cover.

Xinjiang experiences a temperate continental climate that is distinguished by its characteristic precipitation patterns, unlike the monsoon regions where atmospheric circulation is affected by topography and landforms. The average annual precipitation is about 150 mm. With a long winter and short spring and autumn, Xinjiang has abundant sunshine in summer and considerable temperature differences between day and night. With long, cold winters, the region's snow is generally from November to March (Liu et al. 2018).

2.2. Data collection and preprocessing

2.2.1. The enhanced-resolution PM data

NASA's program, Making Earth System Data Records for Use in Research Environments (MEaSUREs), has provided a new version of PM brightness temperature data, named Calibrated Enhanced-Resolution PM Daily EASE-Grid 2.0 (Equal-Area Scalable Earth Grid) Brightness Temperature. These PM data are from Level-2 satellite records of multiple sensors from 1978 to 2021. This enhanced resolution was obtained from the National Snow and Ice Data Center. We used data from 2011 to 2019 (November to March) to explore the feasibility of PM data in SD estimation. In this study, both horizontal (H) and vertical (V) polarizations of the three channels (19, 37, and 91 GHz) and 22 GHz with vertical polarization provided by the special sensor microwave image sensor (SSMIS, F-18) were used. These datasets have two spatial resolutions (19–22 GHz, 6.25 km; 37–91 GHz, 3.125 km). In order to minimize the influence of wet snow, this study exclusively utilized images captured during the early morning (03:52 AM) (Xiao et al. 2021).

After downloading the images from the official website, it was necessary to standardize their resolution. To achieve this, a bilinear interpolation method was employed in the study to resample the spatial resolution data from 6.25 km to 3.125 km.

2.2.2. Ground measurements

Field observation is still the most reliable data for SD estimation. The long-time series SD data in the study area are from the Meteorological Information Center of China's Meteorological Administration (Figure 1). There are 89 sites in Xinjiang, China, covering the daily measured SD data from 2011 to 2019.

To supplement the survey data, the research team of Tianshan Snow and Avalanche (Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences) conducted three field observations of snow characteristics around the Tianshan Mountains from 1 January 2018, to 24 February 2018, and 18 December 2018, to 28 February 2019. The field investigation route covers Urumqi, Changji, Shihezi, Kuitun, Jinghe, Huocheng, Yining, Zhaosu, Gongliu, Nilek, Xinyuan, and Hejing, where moisture, density, particle size, and other snow characteristics of snowpacks are measured by snow fork according to the geographical units of plain, desert, and mountain (Figure 1). All parameters were measured in strict accordance with the standards of the International Association for Snow and Ice Science (IACS), resulting in the collection of a total of 160 data points.

Table 1 shows that the spatial distribution of SD varied significantly in different cover types. In addition, forests, built-up areas, farmland, and grassland are high-value areas for SD, while low-value areas are primarily distributed in bare land and shrubland.

Table 1. Descriptive statistic of SD distribution in various cover types.

Land cover	Sample size	Minimum (cm)	Maximum (cm)	Mean (cm)	Std (cm)	Kurtosis	Skewness	K-S Value	CV (%)
Grassland	8930	1	66	14.3	10.4	3.56	0.96	0.95	72.42
Bareland	7545	1	52	12.4	12.8	3.60	1.26	0.84	103.00
Shrubland	1906	1	52	10.9	8.4	4.21	1.06	0.90	77.07
Forest	1353	1	83	27.6	11.23	3.88	0.53	1.00	40.66
Farmland	25821	1	65	15.8	11.84	3.64	1.02	0.93	74.86
built-up areas	10448	1	58	17.4	10.70	3.06	0.66	0.96	61.37

Note: standard deviation (Std); Coefficient of variation (CV).

2.2.3. Land cover products

Land cover was used to identify the potential impact on SD retrieval. The European Space Agency Climate Change Initiative (ESA-CCI) (https://www.esa-landcover-cci.org/) provides land cover datasets from 1992 to 2015. According to the United Nations (UN) Food and Agriculture Organization (FAO) criterion, land cover types are divided into 22 categories and 36 subcategories. It is also the global land cover product with the most extended time series and high spatial resolution (300 m). It is very suitable for large-scale and long-time series land cover changes.

This study employed ESA Climate Change Initiative Land Cover (CCI-LC) data in 2015 and refined it with ArcGIS 10.8 software based on the boundaries of the study area. A simple majority method was adopted to reclassify the land use types. The study areas contain nine land cover types: grassland, forest, bare land, cropland, shrubland, wetlands, water bodies, glaciers, and built-up areas (Figure 1). In addition to water bodies, glaciers, and wetlands, the study aimed to establish SD retrieval models for the six land cover types mentioned above.

2.2.4. Vegetation fraction data

Previous studies have shown that the land cover fraction poses a potential impact on the accuracy of SD retrieved from PM data. Therefore, this study chose the MODIS Vegetation Continuous Fields product (MOD44B), which can be accessed from the website https://search.earthdata.nasa.gov. With a spatial resolution of 250 m, MOD44B describes the percentage of tree canopy, non-tree vegetation, and non-vegetation in pixels. To match the cloud-free snow products, the original MOD44B data were spliced and resampled to a resolution of 3.125 km.

2.2.5. Topographic data

The topographic parameter was utilized as auxiliary data for the SD estimation. The Shuttle Radar Topography Mission (SRTM3) was used as DEM data, with a spatial resolution of 90 m (https://srtm.csi.cgiar.org/srtmdata/). The extracted DEM-based topographic data include longitude, latitude, altitude, slope, and aspect, which were resampled to a spatial resolution of 3.125 km prior to further processing.

3. Methodology

3.1. Snow detection

Considering that the microwave scatter characteristics of snow are similar to those of frozen soil, cold desert, and precipitation, we used Grody's decision tree to identify the scattered signals of snow (Table 2).

Table 2. Conditions for removing other scattering signatures.

Steps	Conditions
Scattering signature	Tb19V−Tb37V > 0 K
Precipitation	Tb22V > 258 K or 258 K ≥ Tb22 V ≥ 254 K and Tb19V − Tb37V ≤ 2 K
Cold deserts	Tb19V − Tb19H ≥ 18 K and Tb19V − Tb37V ≤ 10 K
Frozen ground	Tb19V − Tb19H ≥ 8 K and Tb19V − Tb37V ≤ 2 K

3.2. Support vector machine

Snow vector machine (SVM) is an algorithm that integrates both computer science and statistical theory. It has the characteristics of substantial generalization ability and high prediction accuracy. It can also obtain effective results in the case of small samples.

In SVM, an appropriate kernel function is significant for modeling complex nonlinear relationships. Equation (1) is a radial basis function: (1) $k(x_i,x)=\exp (\sigma {\Vert x_i,x\Vert }^2)$(1) The parameter σ is ‘speed’ in the core.

3.3. Gaussian process regression

In recent years, as a new machine learning tool, Gaussian process regression (GPR) has been gradually introduced into the field of remote sensing (Sun, Wang, and Xu 2014). GPR is a Bayesian method, which predicts by a mean function m(x) and a kernel (covariance) function k(x_i, x_j) (Equation (2)): (2) $y_i=GPR(m(x_i),k(x_i,x_j))$(2) Common kernel functions are as follows:

Squared Exponential Kernel (Equation (3)): (3) $k(x_i,x_j)={\sigma }_f^2\exp \left[-{\displaystyle \frac{1}{2}{\displaystyle \frac{{(x_i-x_j)}^T(x_i-x_j)}{{\sigma }_l^2}}}\right]$(3) Exponential Kernel (Equation (4)): (4) $k(x_i,x_j)={\sigma }_f^2\exp \left[-{\displaystyle \frac{r}{{\sigma }_l}}\right]$(4)

Matern 5/2 Kernel (Equation (5)): (5) $k(x_i,x_j)={\sigma }_f^2\left(1+{\displaystyle \frac{\sqrt{5}r}{{\sigma }_l}+{\displaystyle \frac{5r^2}{3{\sigma }_l^2}}}\right)\exp \left(-{\displaystyle \frac{\sqrt{5}r}{{\sigma }_l}}\right)$(5) Rational Quadratic Kernel (Equation (6)): (6) $k(x_i,x_j)={\sigma }_f^2{\left(1+{\displaystyle \frac{r^2}{2\alpha {\sigma }_l^2}}\right)}^{-\alpha }$(6) Where (x_i, y_i) is the data in the set, x_i can be understood as the input matrix, and y_i is the output matrix. In addition, the kernel parameters were based on the signal standard deviation σ_f and the characteristic length scale σ_l. Both σ_f and σ_l need to be greater than 0. In addition, r is the Euclidean distance between x_i and x_j; α is a positive-valued scale-mixture parameter.

3.4. Random forest

Random forest (RF) regression is an ensemble process with strong robustness and insensitivity to overfitting noise. While running in parallel, a manifold decision tree was established to reduce the variance in an algorithm. In this matter, the training data are initially divided into several subsets by taking recurrent samples to be trained independently, designated as bootstraps. Afterward, the final output is determined by combining and averaging the results of all models, characterized as aggregation (Equation (7)): (7) $y(x)={\displaystyle \frac{1}{N}\sum _{n=1}^N{\hat{f}}^n(x)}$(7) Where y is the final output and ${\hat{f}}^n(x)$ is n-th bootstrapped training set. N is the number of separate training sets.

3.5. Extreme gradient boosting

Extreme gradient boosting (XGBoost) is an improved algorithm based on gradient-boosting decisions. It can efficiently build boosted trees with parallel operations and solve classification and regression problems. In addition, XGBoost avoids overfitting and optimizing computation. This study adopted the regression function of the algorithm. The algorithm's core is to optimize the parameters of the objective function and implement the machine learning algorithm in the gradient boosting framework. Through parallel tree boosting, XGBoost can quickly and accurately solve data science problems, such as gradient-boosting decision trees and machines (Shin, Son, and Cha 2022).

The objective function of XGBoost usually consists of two parts: training loss and regularization, expressed by Equation (8): (8) $Obj(\theta )=L(\theta )+\mathrm{\Omega }(\theta )$(8) Where L is the training loss function and Ω is the regularization term. The training loss measures the model's performance on the training data. The purpose of the regularization term is to control the complexity of the model, such as overfitting. In addition, the complexity of each tree is typically calculated as follows, Equation (9): (9) $\mathrm{\Omega }(f)=\gamma T+{\displaystyle \frac{1}{2}\lambda \sum _{j=1}^T{w}_j^2}$(9) Where T is the number of leaves and ω is the vector of scores on leaves.

Finally, the structure score of XGBoost is the objective function defined as Equation (10): (10) $Obj=\sum _{j=1}^T\left[G_j{w}_j+{\displaystyle \frac{1}{2}(H_j+\lambda )w_j^2}\right]+\lambda T$(10) Where ω_j is independent of each other. The form is quadratic and the best ω_j for a given structure q(x).

3.6. Light gradient boosting machine

The light gradient boosting machine (LightGBM) is a novel and robust boosting algorithm that combines the advantages of random forest and gradient learning algorithms (Wang, Chang, and Liu 2022). Numerous practices have proved that LightGBM can solve the time-consuming problem of GBDT processing large samples of high-dimensional data, making GBDT more suitable for real industrial systems. LightGBM performs four optimizations on the traditional GBDT ensemble algorithm (Li et al. 2023). The first is to use the histogram algorithm to determine the optimal segmentation point. Second, to address the problem of large data volume, the gradient-based one-side sampling (GOSS) method only uses data instances with noticeable gradient changes to estimate the information gain, effectively saving computational time and space overhead. Third, exclusive feature bundling (EFB) bundles multiple contradictory features into one feature to achieve dimensionality reduction. Fourth, a leaf-wise leaf growth strategy with depth restrictions is adopted. A more detailed description of LightGBM can be found in the References section (Lin et al. 2022; Wang, Chang, and Liu 2022).

3.7. Shapley additive explanations (SHAP)

In SHAP, the contribution of each feature (ψ_i) to the model output v(N) is allocated based on the marginal contribution. Shapely values are represented by Equation (11): (11) ${\varphi }_i={\sum }_{S\subseteq N}{\displaystyle \frac{|S|!(M-|S|-1)!}{M!}}[v(S\cup \{i\})-v(S)]$(11) Where f is the model, M is the number of input features, and N is the set of all input features. The quantity v(S∪{i}) −v(S) expresses, for each single prediction, the deviation of Shapley values from their mean: the contribution of the i-th variable (Tan, Gan, and Wu 2023).

3.8. Estimation approach: theoretical foundation for SD estimation

The SD retrieval model flow can be divided into three steps (Figure 2):

Figure 2. Schematic overview of the modeling workflow.

Step 1: Selection of independent variables. PM brightness temperature (i.e. Tb19V, Tb19H, Tb37V, Tb37H), SD observation, topographic variables (latitude, longitude, elevation, slope, aspect, etc.), land cover types and fractions (bare land, grassland, cropland, forest, built-up areas, the percentage of the tree canopy, non-tree vegetation, and non-vegetation) were used as input variables.

Step 2: Model development. According to input variables and SD observation, the SD retrieval model was constructed using SVM, GPR, RF, XGBoost, and LightGBM. The model exhibiting the highest level of accuracy was further optimized using a 50% discount cross-verification method for different cover areas. To fully tap the performance of each model, this study adopted the appropriate kernel function to establish complex nonlinear relationships between SD and variables for SVM- and GPR-based SD retrieved models. For the RF, XGBoost, and LightGBM-based SD retrieved models, a Bayesian optimization methodology was used to tune the related parameters. Based on the SD inversion model constructed above, the optimal SD inversion model was compared and selected according to the different cover types.

The model selection process followed the guidelines outlined below: For the GPR model, the selection of kernel functions for different land cover types was determined. Specifically, exponential GPR was chosen for grassland, farmland, forest land, and construction. The Matern 5/2 Kernel was selected for bare land, while quadratic rational GPR was used for shrubland. In the case of the SVM function, a radial basis funtion was utilized as the kernel function for grassland, bare land, shrubland, farmland, forest, and built-up areas. In addition, RF, XGBoost, and LightGBM-based SD retrieved models use the tree-structured Parzen estimator (TPE) for hyperparameter optimization, with 100 iterations using the Hyperopt 0.2.7 library of Python 3.7.11(Tables 3 and 4).

Table 3. Hyperparameters tuning for RF model.

Model	Land cover	max_depth	n_estimators	criterion	min_samples_leaf	min_samples_split	max_features
RF	Grassland	7	100	mse	1	2	auto
	Bareland	4	160	mse	1	2	auto
	Shrubland	5	200	mse	1	2	auto
	Forest	12	300	mse	1	2	auto
	Farmland	8	120	mse	1	2	auto
	Built-up areas	2	220	mse	1	2	auto

Table 4. Hyperparameters tuning for XGBoost and LightGBM model.

Model	Land cover	Learning rate	N estimator	Max depth	Colsample by_tree	min_child_weight	gamma	Sub samples
XGBoost	Grassland	0.301	100	6	1	1	0	1
	Bareland	0.028	100	8	0.373	1	0	0.67
	Shrubland	0.015	220	3	0.819	2	0	0.99
	Forest	0.546	280	8	0.784	7	0	0.81
	Farmland	0.134	330	8	0.576	2	0	0.53
	Built-up areas	0.211	260	5	0.774	8	0	0.57
LightGBM	Grassland	0.1	100	1	1.01
	Bareland	0.131	100	2	0.77
	Shrubland	0.133	180	3	0.83
	Forest	0.504	230	6	0.87
	Farmland	0.606	150	6	0.61
	Built-up areas	0.011	170	8	0.73

Step 3: Model validation. We synthesized the SD retrieval results obtained by the best models of various cover types to obtain the SD distribution in Xinjiang. On this foundation, we used R², root mean square error (RMSE), and bias to evaluate the performance of the SD product.

3.9. Performance evaluation

In this study, determination coefficients (R²), RMSE, and mean absolute error (MAE) were selected to measure the accuracy of the evaluation. The R² close to 1 indicated that the model predictions were near perfect. The model with a smaller RMSE and MAE presented better prediction and robustness. The details are shown in Equations (12–14): (12) $R^2={\displaystyle \frac{\sum _{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}{\sqrt{\sum _{i=1}^n{(x_i-\overline{x})}^2{(y_i-\overline{y})}^2}}}$(12) (13) $RMSE=\sqrt{{\displaystyle \frac{1}{n}\sum _{i=1}^n{(x_i-y_i)}^2}}$(13) (14) $MAE={\displaystyle \frac{1}{n}\sum _{i=1}^n|x_i-y_i|}$(14)

Where $\overline{x}$is the average of the predicted value x_i,$\overline{y}$is the average of the actual value y; and n is the number of observations.

4. Results

4.1. Selection of optimal models

Previous research has indicated that incorporating geographic location and topographic variables into the snow depth (SD) retrieval model can greatly enhance the accuracy of spatio-temporal SD estimation. In line with this, our study employed XGBoost, LightGBM, RF, SVM, and GPR models to construct SD retrieval models. These models were augmented with geographic location data (longitude, latitude), topographic variables (elevation, slope, aspect), and cover types (percent tree cover, percent non-tree cover, percent non-vegetated). The estimation accuracy of these five models was then comparatively evaluated.

Figure 3 (a1–a6) reveals the consistency between estimated and observed snow depth (SD) in tree-based regression models. Upon closer examination, while all five SD retrieval models demonstrate a strong correlation with observed SD in the six land cover types (R² > 0.48), both the SVM and GPR models exhibit higher RMSE and MAE compared to XGBoost, LightGBM, and RF. Notably, the results indicate that the SD LightGBM-based retrieval model outperforms the other four models in capturing the intricate relationship between SD and auxiliary variables within the six land cover types.

Figure 3. Comparison of the accuracy obtained with five SD retrieved models.

4.2. Evaluation of model performance

4.2.1. Performance of lightGBM modeling

This section describes the SD LightGBM retrieval model for various cover types. There are multiple hyperparameters in LightGBM. All the above models were fitted to the training set and verified with the verification set. R², MAE, and RMSE were selected as statistical indicators for accuracy evaluation.

PM bands, auxiliary variables, SD observations, and indicators from LightGBM retrieval models differed in six land cover types (Figure 4). Specifically, the SD retrieval model tends to underestimates snow depth and displays the highest accuracy on bare land areas. This can be attributed to the presence of a single-type cover and flat terrain within pixels, which facilitates more accurate estimation. The accuracy of the SD retrieval model was significantly undermined in forests compared to other cover types because the forest canopy weakens the microwave scattering signal of snow.

Figure 4. Scatter plots of the 5 cross-validation of training accuracy for the LightGBM retrieval model with auxiliary variable: (a1)Grassland; (a2)Bareland; (a3)Shrubland; (a4)Forest; (a5)Farmland; (a6)Built-up areas.

4.2.2. Accuracy evaluation under different land covers

Figure 5 shows the scatter density diagrams of the actual SD and estimates of the six cover types using the LightGBM model. The figure shows that LightGBM-based SD products had good accuracy in the SD estimation of each covered area; its accuracy was consistent with the performance of the training model. Specifically, LightGBM-based SD products are more accurate in grassland, bare land, and shrubland than in farmland and built-up areas. Regarding forest cover, there was a substantial deviation between the actual SD and the estimated SD.

Figure 5. Scatter density diagrams of the SD product from LightGBM retrieval model with auxiliary variables: (a1)Grassland; (a2)Bareland; (a3)Shrubland; (a4)Forest; (a5)Farmland; (a6)Built-up areas.

4.3. Accuracy of the reconstructed snow depth product

4.3.1. Accuracy evaluation under different land cover

We used the models mentioned above to reconstruct long-time series SD products (2011–2019) and the SD observation to evaluate product accuracy. The results show that indicator performance fluctuates in both six cover types and interannual variability (Figure 6). Specifically, the accuracy of the models was the lowest in 2013, indicating a high deviation. In contrast, the years 2011 and 2017 displayed higher accuracy levels with relatively lower deviations. Additionally, the estimation accuracy of snow depth (SD) in bare areas was notably higher (R²: 0.784–0.965; MAE: 1.23–6.91 cm; RMSE: 2.16–12.29 cm), while the accuracy in forested areas was the lowest (R²: 0.480–0.792; MAE: 3.27–6.03 cm; RMSE: 4.57–8.75 cm) throughout the years. In addition, it can be inferred from Table 5 that the estimated SD in the six covered areas was lower than the actual SD; estimates of depth that were deeper than the actual SD resulted in more uncertainties.

Figure 6. Time series of statistical indicators for the SD product from LightGBM retrieval model with auxiliary variables: (a1)R²; (a2)RMSE; (a3)MAE.

Table 5. Statistic indices of snow-depth product for different snow-depth ranges in each cover type.

Cover types	Snow depth range	MAE/cm	RMSE/cm	Bias/cm
Grassland	<=10	1.78	3.42	−0.88
	10∼20	2.75	3.59	1.08
	20∼30	4.18	5.64	3.00
	30∼40	7.53	10.54	8.25
	>40	21.36	26.21	20.10
Bareland	<=10	1.39	2.90	−0.43
	10∼20	4.02	5.13	2.18
	20∼30	5.30	8.04	3.77
	30∼40	14.15	17.79	10.92
	>40	10.54	18.14	11.28
Shrubland	<=10	1.79	2.63	−0.96
	10∼20	2.65	3.74	1.67
	20∼30	5.22	7.08	2.47
	30∼40	15.70	18.00	14.29
	>40	25.97	28.47	32.50
Forest	<=10	6.65	9.81	−2.73
	10∼20	3.94	5.29	−1.07
	20∼30	4.73	6.43	−0.36
	30∼40	4.08	5.70	1.11
	>40	7.63	10.79	4.74
Farmland	<=10	1.86	3.66	−0.91
	10∼20	3.55	4.70	0.53
	20∼30	4.80	6.67	3.37
	30∼40	7.76	11.13	6.97
	>40	13.29	19.09	12.95
Build-up areas	<=10	2.66	4.80	−1.80
	10∼20	3.37	4.36	−0.09
	20∼30	4.09	5.49	2.72
	30∼40	7.35	9.45	6.58
	>40	14.55	19.48	12.73

4.3.2. Overall accuracy

The accuracy of the SD product was evaluated consistent with SD observations. The descriptive statistics illustrate that the LightGBM-based SD product with auxiliary variables (R²= 0.782, RMSE = 5.63 cm, MAE = 3.05 cm) achieved high accuracy (Figure 7).

Figure 7. Scatter density diagrams of the SD product from LightGBM retrieval model with auxiliary variable.

To explore the error distribution of the SD products, the SD observation was divided into five ranges (0–10 cm, 10–20 cm, 20–30 cm, 30–40 cm, and above). In each range, the RMSE, MAE, and bias were calculated according to the observed and retrieved SD (Table 6). Table 6 shows that the product overestimated SD when the actual SD was less than 10 cm (Bias = −1.43 < 0, RMSE ≥ 3.73). At the same time, it underestimated SD when the actual SD was higher than 10 cm (Bias ≥ 0). For the actual SD within 10–40 cm, the accuracy was higher; however, the predictive value was higher than for the observation. Furthermore, when the actual SD was greater than 40 cm, all the Bias and RMSE were greater than 10 cm, indicating a considerable underestimation.

Table 6. Statistic indices of snow-depth product for different snow-depth ranges.

	Snow depth range	RMSE/cm	MAE/cm	Bias/cm
LightGBM SD model with topographic variable	<=10	3.73	1.92	−1.43
	10∼20	4.18	3.17	0.06
	20∼30	6.02	4.37	2.05
	30∼40	9.73	7.08	6.24
	>40	19.14	13.21	12.71

4.4. Spatial-temporal of snow depth

Figure 8 shows the spatial distribution of the average SD during winter as derived from the SD product. The map shows that the high SD region (SD ≥ 20 cm) in northern Xinjiang was distributed in the northern and western Tianshan Mountains (112 900.39 km²). The high SD region (SD ≥ 10 cm) in southern Xinjiang was concentrated in the high-altitude area of the southern Kunlun Mountains. Meanwhile, the low-value region of SD was concentrated in the hinterland of the Gurbantonggut Desert, Taklimakan Desert, and surrounding regions. In addition, the maximum SD in northern Xinjiang was higher than that in southern Xinjiang, and the SD in mountainous areas was higher than in the oasis and desert regions.

Figure 8. Spatial distribution of average SD in Xinjiang’ winter from 2011 to 2019.

The spatial pattern of SD trends in Xinjiang has been effectively characterized through the utilization of Theil-Sen median analysis and the Mann-Kendall trend test. Figure 9 shows the interannual SD variation trend of SD product. Notably, approximately 25.89% of Xinjiang exhibited no significant improvement in SD over the period from 2011 to 2019. Around 18.00% of the region experienced no significant degradation, while another 18.00% remained relatively stable during the same period. Specifically, SD distribution with no significant improvement in northern Xinjiang was concentrated primarily in the northern foot of the Tianshan Mountains and the central Junggar Basin. In contrast, SD in southern Xinjiang was concentrated primarily in the southern mountain region. The SD distribution with no significant degradation was mainly distributed over the rims of the Junggar Basin, Yili Valley basin, and southern Tianshan Mountains. Meanwhile, there were stable areas between the areas with no significant improvement or degradation.

Figure 9. Spatial distribution of the variation trend of the inter-annual SD in Xinjiang from 2011 to 2019.

5. Discussion

5.1. Comparisons with existing studies

PM SD retrieved algorithms, including semi-empirical and machine learning algorithms, have been successfully applied to product generation (Hu et al. 2021; Xu et al. 2022; Yue, Tao, and Liyun 2022). Describing the nonlinear relationship between snow depth and multiple predictive variables, the SD estimation models are more accurate and bolstered by a batch of machine learning approaches, including ANN, SVM, and RF (Xiao et al. 2018; Yang et al. 2020; Zaerpour, Adib, and Motamedi 2020). For instance, coupled with PM brightness temperatures, auxiliary variables, and observations, SVM-based models could accurately estimate SD in Xinjiang, China, and Eurasia (Xiao et al. 2018; Yang et al. 2022). Regarding Tibet, China, the Karen River basin in the southwest of Iran, and the Northern Hemisphere, ANN-based models estimated SD with high accuracy (Cao, Yang, and Zhu 2008; Zaerpour, Adib, and Motamedi 2020). RF-based model estimates were highly accurate in northern Sweden, northwestern China, and the Northern Hemisphere (Yang et al. 2020; Zhang et al. 2021). Based on the studies above, the uncertainty of SD ANN-based models depends on the number of hidden layers and nodes. Since RF is considered the black-box model, its accuracy is undermined in terms of low-dimensional data and small samples.

Based on the above representative SD estimation algorithms, this study extended upon tree-based machine learning approaches, namely Random Forest (RF), XGBoost, and LightGBM. These approaches take into account the distribution characteristics of variables and are effective in capturing the nonlinear relationships among them. Through comparative analysis with observed SD data, the LightGBM model demonstrated the highest accuracy across all land cover types, thus being chosen as the final SD estimation model. The LightGBM-based model consistently outperformed other models in terms of accuracy. However, it should be noted that the statistical indices for SD products (Bias, RMSE, and MAE) varied within different land cover types and time periods. Particularly, when the actual SD exceeded 40 cm, the LightGBM-based model tended to underestimate SD. These findings suggest that machine/deep learning-based SD estimation models are well-suited for regions lacking direct observations of snow characteristics. Furthermore, these models have significantly enhanced the accuracy of SD distribution to a certain extent.

5.2. Influence of predictor variables

The RF, XGBoost, and LightGBM algorithms are called ensemble learning. Their learning ability does not have to include the specific relationship between predictors and the response variables (Guo et al. 2023). After setting the corresponding parameters, a result with good accuracy can be obtained. In order to provide further insights into the LightGBM model, this study evaluated the importance of different variables using feature_importance indices. The significance of each variable was then ranked accordingly, as illustrated in Figure 10.

Figure 10. The mean absolute value of the SHAP values and the sum of SHAP value magnitudes over all samples for (a1, b1)Grassland; (a2, b2)Bareland; (a3, b3)Shrubland; (a4, b4)Forest; (a5,b5)Farmland; (a6, b6)Built-up areas of the Catboost-based SD retrieved models.

The aforementioned findings emphasize the importance of considering the intricate topographic characteristics of the desert-oasis-mountain system in Xinjiang, China. By incorporating topographic variables such as longitude, latitude, elevation, slope, aspect, as well as cover conditions including percent non-tree cover (PNT), percent non-vegetated (PNV), and percent tree cover (PTC), the performance of the snow depth (SD) estimation models is significantly improved. For instance, in areas characterized as grassland, bare land, farmland, and built-up areas where SD estimation exhibits high accuracy, latitude has been identified as a crucial factor influencing the accuracy of the SD estimate. On the other hand, the impact of land cover fractions on the accuracy of SD estimation is relatively weak. This is reflected in the fact that the maximum mean SHAP value of the land cover fraction is less than 1. In areas characterized by shrubs and forests, where the accuracy of SD estimates is relatively low, topographic factors and cover conditions exhibit varying degrees of impact on the accuracy of SD estimation. Specifically, height and PNT are the dominant influencing factors for shrub areas, while cover conditions primarily influence SD estimates in forested areas. In contrast, the influence of topographic variables on the accuracy of the snow depth (SD) estimate is relatively minor. This finding indicates that forests and shrubland tend to weaken microwave signals emanating from the snowpack, thereby decreasing the sensitivity of the Passive Microwave (PM) bands to SD. Overall, this study reaffirms the crucial role of topographic variables and cover conditions in accurately estimating SD in regions characterized by complex terrain. These factors are essential for adjusting and improving SD estimates (Wang et al. 2019; Wang et al. 2020; Wei et al. 2019; Wei et al. 2021).

In contrast to previous studies, this research takes the lead in using enhanced PM bands in SD products, thereby significantly improving the spatial resolution to 3.125 km. Figures 8 and 9 show that SD products can provide more detailed SD distribution and obtain more accurate SD and variation characteristics in areas with reliable heterogeneity. Moreover, the experimental results highlight that the sensitivity of polarization difference to SD varies across different land cover types. While the 91 GHz bands are capable of detecting shallow snow cover, they are not incorporated into the SD retrieval model due to limitations imposed by the distinctive spatial distribution of shallow snow as well as the restricted spatial resolution of PM bands (Zhang et al. 2017).

5.3. Limitations and future work

5.3.1. Snow characteristics

While the four snow depth (SD) retrieval models effectively estimate the spatial distribution of SD using predictive variables, they exhibit limitations in capturing variations in snow characteristics and lack robustness across different time scales, such as monthly or yearly variations. Hence, it becomes essential to develop a model that can simulate the evolution of snow characteristics, allowing for an improved understanding of its relationship with microwave brightness temperature (Chen et al. 2020; Kim et al. 2019; Wójcik et al. 2008).

5.3.2. Mixed pixels

Four SD retrieval models were tailored to various cover types. However, with spatial heterogeneity in each pixel of enhanced PM bands, this study used coarse pixels as pure ones, while they were actually mixed pixels, undermining the SD estimate (Markus, Powell, and Wang 2006). Therefore, it is necessary to weight the bright-temperature land cover in the scanning range to reduce the error caused by the underlying surface as much as possible (Collados-Lara et al. 2020; Revuelto et al. 2020; Zheng et al. 2019).

The four snow depth (SD) retrieval models were customized for different land cover types. However, since there is spatial heterogeneity within each pixel of the enhanced Passive Microwave (PM) bands, this study treated coarse pixels as pure ones, even though they actually contained a mixture of land cover types. This approach compromised the accuracy of the SD estimation (Markus, Powell, and Wang 2006). Therefore, it becomes crucial to assign weights to the bright-temperature land cover within the scanning range to minimize errors arising from the underlying surface (Collados-Lara et al. 2020; Revuelto et al. 2020; Zheng et al. 2019).

5.3.3. Snow depth observation

Meteorological sites are sparse and unevenly distributed in Xinjiang, resulting in insufficient ground observations (Zhang et al. 2017). In addition, the functional data representation of the SD was not considered for the lack of continuous-time observations (Wang et al. 2020). In the future, strengthening the fixed-site observation and fine snow survey in a large area will help improve the accuracy of the SD retrieval model.

6. Conclusion

This study analyzed the application of enhanced PM data in SD estimation at temporal and spatial scales. With SVM, GPR, RF, XGBoost, and the LightGBM-based machine learning framework to retrieve SD, the LightGBM model with auxiliary variables presents best accuracy and performance through verification and comparison.

LightGBM-based SD products can describe the distribution characteristics of snow depth in Xinjiang more accurately. However, the accuracy of snow products is influenced by many factors (geographic location, topographic features, and land cover), and the variables’ influence in each cover type is different.

The accuracy of snow depth (SD) products is influenced by several factors such as geographic location, topographic features, and land cover. Furthermore, the influence of variables varies across different land covers. Although the results have demonstrated that approaches based on LightGBM provide highly accurate functional estimation, there are still some limitations. For instance, the coarse resolution PM data is improved through a series of algorithms rather than direct measurement, and the moderate-resolution PM data from scale conversion cannot precisely depict the snow brightness temperature under complex terrain and varied cover conditions. Therefore, the incorporation of deep learning techniques to address uncertainties at spatial and temporal scales becomes crucial for overcoming these limitations and achieving breakthroughs in SD estimation.

Author contributions

All co-authors of this manuscript significantly contributed to all phases of the investigation. They contributed equally to the preparation, analysis, review and editing of this manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This study was supported by the National Natural Science Foundation of China (NSFC Grant 42371146, 42001061), the National Science and Technology Basic Resources Survey Program of China (2019FY100202), the Third Xinjiang Scientific Expedition Program (2022xjkk0602; 2021xjkk1400), the Youth Innovation Promotion Association, CAS (Y970000375) and the K.C.Wong Education Foundation (grant number: GJTD-2020-14).