Ski field accumulated snow depth inversion method based on optimal feature screening constraint stacking neural network and SAR data

Abstract

The invention discloses a northeast ski field accumulated snow depth inversion method based on optimal feature screening and a constraint stacking neural network, and belongs to the technical field of satellite remote sensing image processing and accumulated snow parameter inversion. According to the method, sentry No.1 SAR data is taken as a core data source, multi-source auxiliary data such as weather station snow depth data, atmosphere reanalysis data, geographic information data, vegetation feature data and time data are fused, and an accumulated snow depth inversion initial feature set is constructed; according to the method, an optimal feature combination highly related to the snow depth is screened out by introducing a Spearman correlation coefficient-based feature screening method (FSSCC), and an integrated learning algorithm cNN-Stacking is constructed on this basis, so that adaptive weighted fusion of inversion results of various basic models is realized, and the snow depth inversion precision and stability are improved. An empirical model is adopted in the ski field ski track area, space fusion is conducted on the empirical model and the inversion result of the natural snow area, and accumulated snow depth distribution covering the whole ski field area with the spatial resolution being 20 m is obtained.

Description

Technical Field

[0001] This invention belongs to the technical field of satellite remote sensing image processing and application. Background Technology

[0002] Snow cover is an important component of cryosphere remote sensing, and its spatial distribution and variation characteristics are of great significance for climate change research, water resource assessment, and disaster monitoring. Snow depth is one of the fundamental physical parameters characterizing snow cover status. Accurately obtaining snow depth information helps in the rational assessment and utilization of snow resources, provides data support for the monitoring and early warning of natural disasters such as snowstorms, and is of great significance for ecosystem protection and biodiversity maintenance.

[0003] Snow cover monitoring methods based on optical remote sensing sensors typically offer high spatial resolution, enabling a relatively detailed description of the spatial distribution characteristics of snow cover. However, these methods are highly dependent on cloud cover, precipitation, and illumination conditions, and are prone to data gaps or discontinuous observations under complex weather conditions. In contrast, passive microwave remote sensing technology is unaffected by cloud cover and illumination conditions, enabling continuous acquisition of surface observation data under all-weather conditions, and has been used for snow depth and snow water equivalent estimation. However, limited by the sensor's operating mechanism, its spatial resolution is relatively low, making it difficult to characterize the fine variations in snow depth at local scales. In recent years, the development of active microwave remote sensing technology, especially synthetic aperture radar (SAR), has provided new technical means for high spatial resolution snow depth estimation. Furthermore, an increasing number of scholars have begun to explore the application of SAR data in snow depth estimation. Shi used multi-band, multi-polarization SAR data to analyze the relationship between radar backscattering characteristics and snow physical parameters, thus achieving the estimation of snow water equivalent or snow depth (see Shi J, Dozier J. Estimation of snow water equivalence using SIR-C / X-SAR. I. Inferring snowdensity and subsurface properties[J]. IEEE Transactions on geoscience and remote sensing, 2000, 38(6): 2465-2474.). Leinss et al. found that the differences and ratios of backscattering coefficients of SAR under different frequency bands and polarization modes are related to snow depth, and improved the semi-empirical model accordingly (see Leinss S, Parrella G, Hajnsek I. Snow height determination by polarimetric phasedifferences in X-band SAR data[J]. IEEE Journal of Selected Topics in AppliedEarth Observations and Remote Sensing, 2014, 7(9): 3794-3810.). However, due to the complex structure of snow layers, the relationship between snow depth and radar signals is usually significant and nonlinear. Therefore, the methods based on physical or semi-empirical models still suffer from limitations in accuracy and generalization ability under complex surface conditions. To overcome these shortcomings, some studies have introduced machine learning methods to model the mapping relationship between SAR data and snow parameters.Zhu et al. developed a snow depth (SD) inversion model for farmland and forest areas using SAR data and employing XGBoost, Random Forest (RF), and Support Vector Regression (SVR) algorithms (see Zhu X, Gu L, Li X, et al. Snow depthestimation based on parameter combinations selection and machine learning algorithm using C-band SAR data in northeast China[J]. IEEE Geoscience and Remote Sensing Letters, 2021, 19: 1-5.). Wei et al. integrated the physical processes of forest-snow radiative transfer into a machine learning framework, significantly improving the accuracy and spatiotemporal generalization ability of snow depth estimation in forested regions of the Northern Hemisphere (see Wei Y, Li X, Gu L, et al. Enhancing snow depth estimation in forested regions of the northern hemisphere: A physically-constrained machine learning approach with spatiotemporal dynamics[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2026, 231: 576-594.). In the field of deep learning, Yu et al. proposed a forest snow depth inversion algorithm called FF-DNN, which integrates optimized feature selection and deep neural networks (DNN), improving inversion performance while reducing uncertainty (see Yu G, Gu L, Li X, et al. Forest snow depth estimation based on optimized features and DNN network using C-Band SAR data[J]. IEEE Geoscience and Remote Sensing Letters, 2024, 21: 1-5.).

[0004] This study proposes a snow depth inversion method based on a constrained stacked neural network with optimal feature selection. This method utilizes Sentinel-1 SAR data and various auxiliary data, fully leveraging the advantages of multiple machine learning algorithms to obtain high-precision snow depth inversion results for ski resorts while maintaining a spatial resolution of 20m. Summary of the Invention

[0005] This invention selects Sentinel-1 SAR data, snow depth data from meteorological stations, and other auxiliary data to propose a snow depth inversion algorithm for ski resorts based on feature selection (FS) and constrained neural network stacking (cNN-Stacking).

[0006] The specific steps of the technical solution adopted in this invention are as follows:

[0007] Step 1: Preprocess the data required for snow depth inversion, match all the processed data according to time and space to obtain an initial dataset for snow depth inversion containing all features; the data includes Sentinel-1 SAR data, meteorological station snow depth data, atmospheric reanalysis data, geographic information data, vegetation feature data, time data, and snow cover data.

[0008] Step 2: Using the above data as initial features, construct a cNN-Stacking algorithm model to perform snow depth inversion;

[0009] Step 3: Use the Spearman correlation coefficient-based feature selection method to select all initial features and obtain the feature combination with the highest correlation to snow depth; and use this feature combination to perform snow depth inversion using the cNN-Stacking algorithm model constructed in Step 2.

[0010] Step 4: Invert snow depth in the ski resort area; By applying linear weighting to the snow cover, construct an empirical estimation model for the snow depth in the ski slope area. The calculation formula is as follows:

[0011]

[0012] The estimated results are spatially fused with the snow depth inversion results in non-ski slope areas to obtain complete snow depth inversion results covering the entire ski resort area;

[0013] The construction steps of the cNN-Stacking algorithm model in step two are as follows:

[0014] (a) Match all the processed data from step one according to time and space to obtain a complete training dataset containing all features;

[0015] (b) The set of sample indices for the complete training dataset Randomly divide the data into F non-overlapping subsets, denoted as... ,in This represents the sample set corresponding to the index of the o-th fold cross-validation;

[0016] Based on the above division, the feature set X and the corresponding subsets of the response variable Y in each trade-off are constructed for subsequent cross-validation training; where the response variable Y is the measured snow depth at the snow accumulation station.

[0017] (c) Initialize a space of size A matrix P, where k is the number of base models, is used to store the inversion results of each base model on the samples during cross-validation; the base models are selected from the K-nearest neighbor algorithm model, the random forest algorithm model, the lightweight gradient boosting machine algorithm model, and the gradient boosting decision tree model.

[0018] For each fold Each base model is trained sequentially, using a method other than the current fold during training. Other samples besides those used in training; after training, the model is then used to train the current 0th fold sample. Perform the inversion and store the output results in the corresponding positions in matrix P according to the sample index and model index to complete the construction of matrix P;

[0019] (c) Retrain all base models using the complete training dataset to obtain the final set of base models used in the inversion stage. ;

[0020] Subsequently, matrix P is input into the neural network model, and its output is a lower triangular matrix. In this model, the depth is 5 layers, with 100 nodes in each layer. The activation function is the ELU function, and nonnegativity and normalization constraints are introduced for the fusion weights of each base model. Apply constraints to make it satisfy Under constraints, the weight vectors of each base model It is calculated using the following formula:

[0021] (6)

[0022] in Let be a k-dimensional vector with all elements being 1;

[0023] (d) For any input sample cNN-Stacking adaptively calculates the corresponding weight vector based on its feature information, and then sums the inversion results of each base model according to the weights to finally obtain the snow depth inversion result:

[0024] (7)

[0025] Preferably, the selection and preprocessing methods for Sentinel-1 SAR data, meteorological station snow depth data, atmospheric reanalysis data, geographic information data, vegetation characteristic data, temporal data, and snow cover data in step one are as follows:

[0026] (1) Sentinel-1 SAR data was preprocessed using IW GRD level data with a spatial resolution of 20m. The processing included: orbit correction, thermal noise removal, radiometric calibration, speckle filtering, terrain correction and decibel conversion.

[0027] (2) Snow depth data for meteorological stations are from the China Meteorological Science Data Service Center;

[0028] (3) The atmospheric reanalysis data were calculated using the fifth generation of global climate atmospheric reanalysis results calculated from observational and model data from the European Centre for Medium-Range Weather Forecasts. The spatial resolution was approximately 11 km and the temporal resolution was hourly. The atmospheric reanalysis data preprocessing involved cropping and resampling according to the selected study area. The spatial resolution of the processed data was 20 m.

[0029] (4) Geographic information data selection: SRTMDEM 90M data, providing elevation data with a spatial resolution of 90 m; Geographic information data preprocessing: obtaining slope and aspect data through ENVI software, and then performing cropping, splicing and resampling according to the selected study area, with the spatial resolution of the processed data being 20 m.

[0030] (5) Vegetation feature data includes forest cover data, provided by the annual product of MOD44B VCF, with a spatial resolution of 250m. The vegetation feature data preprocessing was performed using ModisTool software to obtain tree cover, non-tree cover and non-vegetation (bare ground) cover data, and cropping, splicing and resampling were performed according to the selected study area. The spatial resolution of the processed data was 20m.

[0031] (6) The time data comes from meteorological station observation data. The period from November to March of the following year is a snow cover cycle. When the snow depth is greater than 0 for the first five consecutive days, the first day of the continuous process is defined as the starting day of the snow cover days and the snow cover days are recorded as 1. Subsequently, within the snow cover cycle, the snow cover days increase day by day in the order of the date until the end of the snow cover cycle.

[0032] (7) The snow cover data comes from Sentinel-2 optical data. Based on Sentinel-2 L2A level products, the optical images are processed by combining Normalized Difference Snow Index (NDSI) and Normalized Difference Forest Snow Index (NDFSI) to finally generate a snow cover data product with a spatial resolution of 20 m.

[0033] Preferably, the characteristic data obtained from the preprocessed Sentinel-1 SAR data includes: radar local incidence angle LIA, product local incidence angle PLIA, ellipsoidal incidence angle IAFE, and backscattering coefficients Sigma for different polarization modes. VH Gamma VH Beta VH Sigma VV Gamma VV Beta VV and the corresponding decibel (dB) Sigma VH dB Gamma VH dB Beta VH dB Sigma VV dB Gamma VV dB Beta VV Based on this, the following are added: polarization ratio dB Sigma PR, dB Gamma PR, dB Beta PR; polarization difference dB Sigma PD, dB Gamma PD, dB Beta PD; polarization sum dB Sigma PS, dB Gamma PS, dB Beta PS; normalized polarization index dB Sigma NDPI, dB Gamma NDPI, dB Beta NDPI; and radar vegetation index RVI2, dB RVI2. The calculation formula is as follows:

[0034] (1)

[0035] (2)

[0036] (3)

[0037] (4)

[0038] (5)

[0039] Among them BC VV BC VH BC represents the backscattering coefficients of VV polarization and VH polarization in the Sigma, Gamma, or Beta bands, respectively. PR BCPD BC PS BC NDPI These represent the polarization ratio, polarization difference, polarization sum, and normalized polarization index for the Sigma, Gamma, or Beta bands, respectively. RVI2 represents the radar vegetation index. , These represent the backscattering coefficients of VV polarization and VH polarization in the Sigma band, respectively;

[0040] The characteristic data of snow depth data from meteorological stations include: latitude, longitude, snow depth, and date;

[0041] The characteristic data of atmospheric reanalysis data include: 2m temperature, snow density, and snow albedo;

[0042] Geographic information data features include elevation, longitude, latitude, slope, and aspect.

[0043] Vegetation characteristic data include: tree cover, non-tree cover, and non-vegetation (bare ground) cover;

[0044] The features of time data include: year, month, day, date, and number of days with snow cover;

[0045] The characteristic data of snow cover data is the spatial distribution information of snow cover in the ski slope area of ​​the ski resort.

[0046] Preferably, in step two, batch normalization is introduced into each layer of the neural network model, and a 50% dropout rate is applied to the hidden layers; the Adam optimizer is used for training, and an early stopping strategy is adopted: if the model performance does not improve after 10 consecutive epochs, training is stopped.

[0047] Preferably, in step three, the feature selection method based on Spearman correlation coefficient is used to select all initial features to obtain the feature combination with the highest correlation to snow depth. The specific operation is as follows:

[0048] (1) The initial dataset for snow depth inversion prepared in step one is used to standardize all input data for each feature dimension of the dataset using the formula below;

[0049] (8)

[0050] Where x is the input data, and μ and σ are the mean and variance of each feature data, respectively;

[0051] (2) Feature selection based on Spearman correlation coefficient: For the initial feature dataset of snow depth inversion, the Spearman correlation coefficient between each feature and snow depth is calculated, and the features are ranked according to the magnitude of the absolute value of the correlation coefficient to evaluate the degree of association between different features and snow depth, and to provide a basis for subsequent feature selection; the formula for calculating the Spearman correlation coefficient is as follows:

[0052] (9)

[0053] in and These are the rankings of features x and y, respectively. and This represents the average ranking.

[0054] Root mean square error, mean absolute error, and coefficient of determination were selected as model performance evaluation indicators. The number of features in the snow depth inversion model was gradually increased according to the Spearman correlation coefficient between the features and the snow depth, and then input into the cNN-Stacking algorithm model for training. The changes in the model performance evaluation indicators were observed. When the model evaluation indicators no longer improved or tended to stabilize as the number of features increased, the feature addition process was stopped, and the current feature subset was determined as the optimal feature combination for snow depth inversion.

[0055] Further optimized, the optimal feature combination obtained in step three after feature filtering includes the radar local incidence angle (LIA) from Sentinel-1 data and the backscattering coefficients (Sigma) for different polarization modes. VH Gamma VH Beta VH and the corresponding decibel (dB) Sigma VH dB Gamma VH dB Beta VH dB Sigma PS, dB Gamma PS; 2m temperature, snow density, and snow albedo from ERA5 atmospheric reanalysis data; elevation, longitude, slope, and aspect from SRTMDEM 90M data; tree cover, non-tree cover, and non-vegetation cover from MOD44B VCF data; month, date, and number of days with snow cover from meteorological station observations.

[0056] The beneficial effects of this invention are:

[0057] The purpose of this invention is to achieve high-precision snow depth inversion in ski resort areas, enabling effective monitoring of snow depth at the remote sensing level. Combining Sentinel-1 SAR data, historical snow depth observation data from meteorological stations, and multi-source auxiliary data, a snow depth inversion method FS based on a constrained stacked neural network with optimal feature selection is proposed. SCC The cNN-Stacking approach, by introducing a feature selection strategy based on Spearman correlation coefficient, filters initial features and then uses the cNN-Stacking algorithm to adaptively weight and fuse the inversion results of multiple base models, thereby improving the accuracy and stability of snow depth inversion. Simultaneously, addressing the issue of snow property differences caused by artificial snowmaking in ski resort areas, an empirical estimation model is constructed by combining snow cover information and spatially fused with inversion results from non-ski resort areas. This establishes a snow depth inversion method suitable for ski resorts, obtaining spatially continuous snow depth data with a spatial resolution of 20m. Compared with existing snow depth products, this method has better applicability and completeness in ski resort areas, more accurately characterizing the spatial distribution characteristics of snow cover. This is beneficial for ski resort operation and management, rational utilization of snow resources, and research related to regional hydrology and climate change, providing reliable data support for related applications. Attached Figure Description

[0058] Figure 1 shows the research area, meteorological station, and ski resort location selected for this invention.

[0059] Figure 2 is based on FS. SCC Flowchart of the snow depth inversion algorithm of the cNN-Stacking algorithm.

[0060] Figure 3 is based on FS. SCC The snow depth inversion result of the -cNN-Stacking algorithm is shown in the figure at Wanfeng Ski Resort in Tonghua City, Jilin Province.

[0061] Figure 4 is based on FS SCC The snow depth inversion result of the -cNN-Stacking algorithm is shown in the image of the snow depth inversion at Lianhuashan Ski Resort in Changchun City, Jilin Province. Detailed Implementation

[0062] The technical solution of the present invention will be further explained and described below with reference to specific embodiments.

[0063] The specific steps of the embodiment are as follows:

[0064] Step 1: Data Preprocessing. The remote sensing and meteorological data required for snow cover inversion are preprocessed. The data includes Sentinel-1 data from November to March of each year from 2019 to 2023, January 11, 2025, and February 16, 2025, meteorological station data, ERA5 atmospheric reanalysis data, SRTMDEM 90M data, MOD44B VCF data, and temporal data.

[0065] (1) Sentinel-1 SAR data preprocessing: Sentinel-1 IW GRD level data from 46 training and test sites and 11 verification sites in Northeast China from November to March of the following year from 2019 to 2023, as well as two Sentinel-1 IW GRD level data from the Wanfeng Ski Resort area on January 11, 2025 and the Lianhuashan Ski Resort area on February 16, 2025, were preprocessed with a spatial resolution of 20m.

[0066] The preprocessing process includes: orbit correction, thermal noise removal, radiometric calibration, speckle filtering, terrain correction, and decibel conversion. The preprocessed feature data includes: radar local incidence angle (LIA) within the study area, product local incidence angle (PLIA), ellipsoidal incidence angle (IAFE), and backscattering coefficients (Sigma) for different polarizations. VH Gamma VH Beta VH Sigma VV Gamma VV Beta VV and the corresponding decibel-scaled backscattering coefficient dB Sigma VH dB Gamma VH dB Beta VH dB Sigma VV dB Gamma VV dB Beta VV Based on this, the following parameters are added: polarization ratio dB SigmaPR, dB GammaPR, dB BetaPR; polarization difference dB SigmaPD, dB GammaPD, dB BetaPD; polarization sum dB SigmaPS, dB GammaPS, dB BetaPS; normalized polarization index dB SigmaNDPI, dB GammaNDPI, dB BetaNDPI; and radar vegetation index RVI2, dB RVI2. The calculation formulas are as follows:

[0067] (1)

[0068] (2)

[0069] (3)

[0070] (4)

[0071] (5)

[0072] Among them BC VV BC VH BC represents the backscattering coefficients of VV and VH polarizations in the Sigma, Gamma, or Beta bands, respectively. PR BC PD BC PS BC NDPI These represent the polarization ratio, polarization difference, polarization sum, and normalized polarization index for the Sigma, Gamma, or Beta bands, respectively. RVI2 represents the radar vegetation index. , These represent the backscattering coefficients of VV polarization and VH polarization in the Sigma band, respectively;

[0073] (2) ERA5 Atmospheric Reanalysis Data Preprocessing: ERA5 atmospheric reanalysis data from 46 training and testing stations and 11 validation stations in Northeast China from November to March of each year from 2019 to 2023 were selected. ERA5 atmospheric reanalysis data from the Wanfeng Ski Resort area on January 11, 2025, and the Lianhuashan Ski Resort area on February 16, 2025, were also selected. The above data were spatially cropped according to the study area, and a unified spatial resolution of 20m was achieved through resampling to obtain meteorological characteristic data for snow depth inversion. The characteristic variables included 2m temperature, snow density, and snow albedo.

[0074] (3) SRTM DEM 90 m data preprocessing: SRTM DEM 90 m elevation data corresponding to 46 training and testing sites and 11 verification sites in Northeast China from November to March of each year from 2019 to 2023 were selected. SRTM DEM 90 m data for the Wanfeng Ski Resort area on January 11, 2025 and the Lianhuashan Ski Resort area on February 16, 2025 were also selected. The selected elevation data were cropped, stitched, and resampled. Slope and aspect information were calculated based on the elevation data to finally obtain topographic feature data for snow depth inversion, including elevation, longitude, latitude, slope, and aspect.

[0075] (4) MOD44B VCF data preprocessing: MOD44B VCF annual vegetation cover data corresponding to 46 training and testing sites and 11 validation sites in Northeast China from November to March of the following year from 2019 to 2023 were selected. MOD44B VCF data of Wanfeng Ski Resort area on January 11, 2025 and Lianhuashan Ski Resort area on February 16, 2025 were also selected. The above data were cropped, spliced ​​and resampled to extract vegetation features such as tree cover, non-tree cover and non-vegetation (bare ground) cover, which were used to characterize the spatial distribution of vegetation in the study area.

[0076] (5) Time data preprocessing: Snow cover days data from 46 training and testing sites and 11 verification sites in Northeast China from November to March of the following year from 2019 to 2023 were selected, sorted and preprocessed; at the same time, combined with the first snow date information of Wanfeng Ski Resort and Lianhuashan Ski Resort in the winter of 2024, the snow cover days corresponding to Wanfeng Ski Resort area on January 11, 2025 and Lianhuashan Ski Resort area on February 16, 2025 were calculated and used as the time feature input for snow depth inversion.

[0077] (6) Snow cover data preprocessing: Sentinel 2 L2A level products of Wanfeng Ski Resort area on January 11, 2025 and Lianhuashan Ski Resort area on February 16, 2025 were selected. The optical images were processed by combining Normalized Difference Snow Index (NDSI) and Normalized Difference Forest Snow Index (NDFSI) to obtain the spatial distribution information of snow cover in the ski resort ski slope area. Finally, a snow cover data product with a spatial resolution of 20 m was generated.

[0078] All processed data are matched according to time and space to obtain a complete training dataset containing all features.

[0079] Step 2: Complete the construction of the cNN-Stacking algorithm. The implementation steps of the algorithm are as follows:

[0080] (a) The set of sample indices for the complete training dataset Randomly divide the data into F non-overlapping subsets, denoted as... ,in Let X represent the sample set corresponding to the index of the o-th fold cross-validation. Based on the above partitioning, construct the corresponding subsets of feature set X and response variable Y in each fold for subsequent cross-validation training.

[0081] (b) Initialize a space of size Matrix P and basic model , where k is the number of base models, used to store the inversion results of each base model on the samples during the cross-validation process.

[0082] In this embodiment, the K-Nearest Neighbors (KNN), Random Forest (RF), Lightweight Gradient Boosting Machine (LightGBM), and Gradient Boosting Decision Tree (GBDT) models are selected as the base models. Therefore, k=4 in this embodiment.

[0083] For each fold Each base model is trained sequentially, using a method other than the current fold during training. Other samples besides those used in training; after training, the model is then used to train the current 0th fold sample. Perform the inversion and store the output in the corresponding position in matrix P.

[0084] Through the cross-validation process described above, we can obtain the inversion results generated by each base model for each sample without it participating in its own training, thereby effectively avoiding the problem of information leakage.

[0085] (c) After constructing matrix P, retrain all base models using the complete training dataset to obtain the final set of base models used in the inversion stage. .

[0086] Subsequently, using matrix P as input, a neural network model is trained, and its output is a lower triangular matrix. In the neural network model, the number of hidden layer nodes is fixed at 100, the network depth is set to 5 layers, the activation function is ELU function, batch normalization is introduced in each layer to accelerate convergence, and a 50% dropout rate is applied in the hidden layer to suppress overfitting.

[0087] Within this framework, nonnegativity and normalization constraints are introduced for the fusion weights of each base model. Apply constraints to make it satisfy Under constraints, the weight vectors of each base model It is calculated using the following formula:

[0088] (6)

[0089] in It is a k-dimensional vector with all elements being 1.

[0090] (d) For any input sample cNN-Stacking adaptively calculates the corresponding weight vector based on its feature information, and then sums the inversion results of each base model according to the weights to finally obtain the snow depth inversion result:

[0091] (7)

[0092] Step 3: Feature Selection based on Spearman correlation coefficient (FS) SCC The method filters all initial features to obtain the feature combination with the highest correlation to snow depth, and completes the FS (Feature Search). SCC The construction of the cNN-Stacking algorithm is shown in the following steps:

[0093] (1) Input data: The Sentinel-1 data, meteorological station snow depth data, atmospheric reanalysis data, geographic information data, vegetation characteristic data and time data of 46 training and test stations in Northeast China from November to March of the following year from 2019 to 2023 mentioned in step one are matched to create the initial dataset for snow depth inversion.

[0094] The initial dataset includes features such as radar local incidence angle (LIA), product local incidence angle (PLIA), ellipsoidal incidence angle (IAFE), and backscattering coefficients (Sigma) for different polarization modes, all derived from Sentinel-1 data. VH Gamma VH Beta VH Sigma VV Gamma VV Beta VV and the corresponding decibel (dB) Sigma VH dB Gamma VH dB Beta VH dBSigma VV dB Gamma VV dB Beta VVBased on this, the following data are added: polarization ratio dB Sigma PR, dB Gamma PR, dB BBeta PR; polarization difference dB Sigma PD, dB Gamma PD, dB Beta PD; polarization sum dB Sigma PS, dB Gamma PS, dB Beta PS; normalized polarization index dB Sigma NDPI, dB Gamma NDPI, dB Beta NDPI; radar vegetation index RVI2, dB RVI2; 2m temperature, snow density, and snow albedo from ERA5 atmospheric reanalysis data; elevation, longitude, latitude, slope, and aspect from SRTMDEM90M data; tree cover, non-tree cover, and non-vegetation (bare ground) cover from MOD44B VCF data; and year, month, day, date, and number of days with snow cover from meteorological station observation data.

[0095] (2) Standardize all input data for each feature dimension using the formula below.

[0096] (8)

[0097] Where x represents the input data, and μ and σ are the mean and variance of each feature data, respectively.

[0098] (3) Feature selection based on Spearman correlation coefficient: For the initial feature dataset of snow depth inversion, the Spearman correlation coefficient between each feature and snow depth is calculated, and the features are ranked according to the magnitude of the absolute value of the correlation coefficient to evaluate the degree of association between different features and snow depth, and to provide a basis for subsequent feature selection. Spearman's rank correlation coefficient (SCC) is a nonparametric statistical index used to measure the monotonic correlation between two variables. This method weakens the influence of outliers and scale differences on the correlation calculation results by converting the original variable values ​​into their corresponding ranks. Specifically, the Spearman correlation coefficient is calculated by ranking the two variables in the rank variable space and calculating their Pearson correlation coefficient to quantify the degree of correlation between them. Its value ranges from −1 to 1, where the larger the absolute value of the correlation coefficient, the stronger the correlation between the corresponding feature and snow depth; a positive correlation coefficient indicates a positive correlation, and a negative correlation indicates a negative correlation. The calculation formula is as follows:

[0099] (9)

[0100] in and These are the rankings of features x and y, respectively. and This represents the average ranking.

[0101] Root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (R²) were selected as performance evaluation metrics for the model. The number of features in the snow depth inversion model was progressively increased according to the Spearman correlation coefficient between features and snow depth, and then fed into the cNN-Stacking algorithm for training. The changes in RMSE, MAE, and R² were observed. When the model evaluation metrics no longer improved or tended to stabilize as the number of features increased, the feature addition process was stopped, and the current feature subset was determined as the optimal feature combination for snow depth inversion.

[0102] After feature selection, the optimal feature combination obtained includes the radar local incidence angle (LIA) from Sentinel-1 data and the backscattering coefficients (Sigma) for different polarization modes. VH Gamma VH Beta VH and the corresponding decibel (dB) Sigma VH dB Gamma VH dB Beta VH dB Sigma PS, dB Gamma PS; 2m temperature, snow density, and snow albedo from ERA5 atmospheric reanalysis data; elevation, longitude, slope, and aspect from SRTMDEM 90M data; tree cover, non-tree cover, and non-vegetation (bare ground) cover from MOD44B VCF data; month, date, and number of days with snow cover from meteorological station observations.

[0103] (4) Complete FS SCC - Construction of the cNN-Stacking algorithm. A cNN-Stacking snow depth inversion model for Northeast China was constructed using data from 46 training and testing sites in Northeast China from November to March of each year between 2019 and 2023, after feature selection. The constructed FS... SCC The optimal features for the cNN-Stacking algorithm are 22. The cNN-Stacking algorithm selects four regression models—KNN, RF, LightGBM, and GBDT—as its base models. In terms of network structure, the number of hidden layer nodes is fixed at 100, the network depth is set to 5 layers, and the Adam optimizer is used for training. When the validation loss no longer improves, the learning rate is automatically adjusted to improve training stability. The ELU function is used as the activation function. To accelerate convergence, batch normalization is introduced in each layer, and a 50% dropout rate is applied in the hidden layers to suppress overfitting. An early stopping strategy is adopted during model training: if the model performance does not improve after 10 consecutive epochs, training is stopped.

[0104] Step 4: Snow depth inversion in the ski resort area.

[0105] Field investigations of snow cover in the ski resort area revealed that the physical properties (such as density, crystal form, and layering structure) of the snow on the ski slopes differed significantly from the surrounding natural snow due to long-term human intervention such as artificial snowmaking, snow grooming, and repeated compaction. Compared to natural snow, artificial snow typically has higher density and a more stable snow layer structure, leading to significant changes in its characteristics. Since the snow depth samples in the training dataset used in this study were all derived from natural snow observations, the inversion model built based on this data could not accurately characterize the remote sensing response features of artificial snow in the ski slope area, resulting in considerable uncertainty in the snow depth inversion of the ski slope area. Therefore, the snow depth on the ski slopes and the snow depth outside the ski slopes were processed separately.

[0106] (1) Snow depth outside the snow track: According to steps one to three, FS is used. SCC -cNN-Stacking algorithm for snow depth inversion;

[0107] (2) Snow depth on the ski slope: An experience-based compensation estimation method is introduced. Based on multiple field surveys and existing engineering experience, the snow depth in the ski slope area of ​​the ski resort is usually maintained between approximately 50 cm and 60 cm. Considering the spatial differences in snow cover within the ski slope area, the snow cover data product mentioned in step one is selected as auxiliary information, with a value range of 0 to 1, to characterize the relative integrity of the snow distribution within the ski slope area. By linearly weighting the snow cover, an empirical estimation model for the snow depth in the ski slope area of ​​the ski resort is constructed, and its calculation formula is as follows:

[0108] (10)

[0109] The formula uses 50cm as the baseline value for snow depth in the ski area, and adjusts the snow depth within the range of 50-60cm according to the changes in snow cover, thereby reflecting the spatial differences in snow distribution in the ski area to a certain extent.

[0110] Using the FS introduced above SCC The cNN-Stacking algorithm was used to invert snow depth in the Wanfeng Ski Resort area on January 11, 2025, and the Lianhuashan Ski Resort area on February 16, 2025. The specific steps are as follows:

[0111] (1) Data Matching and Dataset Construction: The preprocessed multi-source remote sensing data of Wanfeng Ski Resort area on January 11, 2025 and Lianhuashan Ski Resort area on February 16, 2025 were matched and saved in a format that can be processed by the snow depth inversion algorithm. FS SCC The 22 features selected by the cNN-Stacking algorithm include the radar local incidence angle (LIA) from Sentinel-1 data and the backscattering coefficients (Sigma) for different polarization modes. VH Gamma VH Beta VH and the corresponding decibel (dB) Sigma VH dB Gamma VH dB Beta VH dB Sigma PS, dB Gamma PS; 2m temperature, snow density, and snow albedo from ERA5 atmospheric reanalysis data; elevation, longitude, slope, and aspect from SRTMDEM 90M data; tree cover, non-tree cover, and non-vegetation (bare ground) cover from MOD44B VCF data; month, date, and number of days with snow cover from meteorological station observations.

[0112] (2) Snow depth inversion and result fusion: The processed feature data is input into the trained cNN-Stacking model to perform snow depth inversion on the natural snow cover of the ski resort area, obtaining the spatial distribution results of snow depth in non-ski slope areas. Subsequently, combined with the snow cover products of Wanfeng Ski Resort area on January 11, 2025 and Lianhuashan Ski Resort area on February 16, 2025, the snow depth of the ski slope area is empirically estimated, and the estimation results are spatially fused with the snow depth inversion results of non-ski slope areas to obtain a complete snow depth inversion result covering the entire ski resort area. Figures 3-4 ).

[0113] Experimental results:

[0114] Based on Sentinel-1 SAR data, meteorological station data, and multi-source auxiliary data, a snow depth inversion model for ski resort areas in Northeast China was constructed, and FS was used. SCC The cNN-Stacking algorithm was used to validate the model's performance. Experimental results show that on the test dataset, the RMSE of this inversion model is 3.61 cm, the MAE is 2.70 cm, and the R² is 0.79, indicating that the model has a good fitting ability for snow depth.

[0115] A validation dataset was constructed using observational data from 11 independent validation sites to evaluate the model's generalization ability. The validation results show that the RMSE of the inversion model on the validation dataset is 4.41 cm, and the MAE is 3.37 cm, indicating that the model maintains relatively stable inversion accuracy even on data not used for training. Furthermore, using Lianhuashan Ski Resort on February 16, 2025, as a typical case study, the inverted snow depth results were compared with the measured snow depth at the same time. The RMSE was 6.32 cm, and the MAE was 5.32 cm, further validating the feasibility of this method in real-world ski resort applications.

Claims

1. A method for inverting snow depth in ski resorts based on a constrained stacked neural network with optimal feature selection and SAR data, characterized in that, The method includes the following steps: Step 1: Preprocess the data required for snow depth inversion, match all the processed data according to time and space to obtain an initial dataset for snow depth inversion containing all features; the data includes Sentinel-1 SAR data, meteorological station snow depth data, atmospheric reanalysis data, geographic information data, vegetation feature data, time data, and snow cover data. Step 2: Using the above data as initial features, construct a cNN-Stacking algorithm model to perform snow depth inversion; Step 3: Use the Spearman correlation coefficient-based feature selection method to select all initial features and obtain the feature combination with the highest correlation to snow depth; and use this feature combination to perform snow depth inversion using the cNN-Stacking algorithm model constructed in Step 2. Step 4: Invert snow depth in the ski resort area; By applying linear weighting to the snow cover, construct an empirical estimation model for the snow depth in the ski slope area. The calculation formula is as follows: The estimated results are spatially fused with the snow depth inversion results in non-ski slope areas to obtain complete snow depth inversion results covering the entire ski resort area; The construction steps of the cNN-Stacking algorithm model described in step two are as follows: (a) Match all the processed data from step one according to time and space to obtain a complete training dataset containing all features; (b) The set of sample indices for the complete training dataset Randomly divide the data into F non-overlapping subsets, denoted as... ,in This represents the sample set corresponding to the index of the o-th fold cross-validation; Based on the above division, the feature set X and the corresponding subsets of the response variable Y in each trade-off are constructed for subsequent cross-validation training; where the response variable Y is the measured snow depth at the snow accumulation station. (c) Initialize a space of size A matrix P, where k is the number of base models, is used to store the inversion results of each base model on the samples during cross-validation; the base models are selected from the K-nearest neighbor algorithm model, the random forest algorithm model, the lightweight gradient boosting machine algorithm model, and the gradient boosting decision tree model. For each fold Each base model is trained sequentially, using a method other than the current fold during training. Other samples besides those used in training; after training, the model is then used to train the current 0th fold sample. Perform the inversion and store the output results in the corresponding positions in matrix P according to the sample index and model index to complete the construction of matrix P; (c) Retrain all base models using the complete training dataset to obtain the final set of base models used in the inversion stage. ; Subsequently, matrix P is input into the neural network model, and its output is a lower triangular matrix. In this model, the depth is 5 layers, with 100 nodes in each layer. The activation function is the ELU function, and nonnegativity and normalization constraints are introduced for the fusion weights of each base model. Apply constraints to make it satisfy Under constraints, the weight vectors of each base model It is calculated using the following formula: (6) in Let be a k-dimensional vector with all elements being 1; (d) For any input sample cNN-Stacking adaptively calculates the corresponding weight vector based on its feature information, and then sums the inversion results of each base model according to the weights to finally obtain the snow depth inversion result: (7)。 2. The ski resort snow depth inversion method based on optimal feature selection, constrained stacked neural network, and SAR data as described in claim 1, is characterized in that, The selection and preprocessing methods for Sentinel-1 SAR data, meteorological station snow depth data, atmospheric reanalysis data, geographic information data, vegetation characteristic data, temporal data, and snow cover data in Step 1 are as follows: (1) Sentinel-1 SAR data was preprocessed using IW GRD level data with a spatial resolution of 20m. The processing included: orbit correction, thermal noise removal, radiometric calibration, speckle filtering, terrain correction and decibel conversion. (2) Snow depth data for meteorological stations are from the China Meteorological Science Data Service Center; (3) The atmospheric reanalysis data were calculated using the fifth generation of global climate atmospheric reanalysis results calculated from observational and model data from the European Centre for Medium-Range Weather Forecasts. The spatial resolution was approximately 11 km and the temporal resolution was hourly. The atmospheric reanalysis data preprocessing involved cropping and resampling according to the selected study area. The spatial resolution of the processed data was 20 m. (4) Geographic information data selection: SRTMDEM 90M data, providing elevation data with a spatial resolution of 90 m; Geographic information data preprocessing: obtaining slope and aspect data through ENVI software, and then performing cropping, splicing and resampling according to the selected study area, with the spatial resolution of the processed data being 20 m. (5) Vegetation feature data includes forest cover data, provided by the annual product of MOD44B VCF, with a spatial resolution of 250m. The vegetation feature data preprocessing was performed using ModisTool software to obtain tree cover, non-tree cover and non-vegetation (bare ground) cover data, and cropping, splicing and resampling were performed according to the selected study area. The spatial resolution of the processed data was 20m. (6) The time data comes from meteorological station observation data. The period from November to March of the following year is a snow cover cycle. When the snow depth is greater than 0 for the first five consecutive days, the first day of the continuous process is defined as the starting day of the snow cover days and the snow cover days are recorded as 1. Subsequently, within the snow cover cycle, the snow cover days increase day by day in the order of the date until the end of the snow cover cycle. (7) The snow cover data comes from Sentinel-2 optical data. Based on Sentinel-2 L2A level products, the optical images are processed by combining Normalized Difference Snow Index (NDSI) and Normalized Difference Forest Snow Index (NDFSI) to finally generate a snow cover data product with a spatial resolution of 20 m.

3. The ski resort snow depth inversion method based on optimal feature selection, constrained stacked neural network, and SAR data as described in claim 2, is characterized in that... The preprocessed feature data from Sentinel-1 SAR data includes: radar local incidence angle (LIA), product local incidence angle (PLIA), ellipsoidal incidence angle (IAFE), and backscattering coefficients (Sigma) for different polarizations. VH Gamma VH Beta VH Sigma VV Gamma VV Beta VV and the corresponding decibel (dB) Sigma VH dB Gamma VH dB Beta VH dB Sigma VV dBGamma VV dB Beta VV Based on this, the following are added: polarization ratio dB Sigma PR, dB Gamma PR, dB Beta PR; polarization difference dB Sigma PD, dB Gamma PD, dB Beta PD; polarization sum dB Sigma PS, dB Gamma PS, dB Beta PS; normalized difference vegetation index dB Sigma NDPI, dB Gamma NDPI, dB Beta NDPI; and radar vegetation index RVI2, dB RVI2. The calculation formula is as follows: (1) (2) (3) (4) (5) Among them BC VV BC VH BC represents the backscattering coefficients of VV polarization and VH polarization in the Sigma, Gamma, or Beta bands, respectively. PR BC PD BC PS BC NDPI These represent the polarization ratio, polarization difference, polarization sum, and normalized polarization index for the Sigma, Gamma, or Beta bands, respectively. RVI2 represents the radar vegetation index. , These represent the backscattering coefficients of VV polarization and VH polarization in the Sigma band, respectively; The characteristic data of snow depth data from meteorological stations include: latitude, longitude, snow depth, and date; The characteristic data of atmospheric reanalysis data include: 2m temperature, snow density, and snow albedo; Geographic information data features include elevation, longitude, latitude, slope, and aspect. Vegetation characteristic data include: tree cover, non-tree cover, and non-vegetation (bare ground) cover; The features of time data include: year, month, day, date, and number of days with snow cover; The characteristic data of snow cover data is the spatial distribution information of snow cover in the ski slope area of ​​the ski resort.

4. The ski resort snow depth inversion method based on optimal feature selection, constrained stacked neural network, and SAR data as described in claim 1, is characterized in that... In step two, batch normalization was introduced into each layer of the neural network model, and a 50% dropout rate was applied to the hidden layers. The Adam optimizer was used for training, and an early stopping strategy was adopted: if the model performance did not improve after 10 consecutive epochs, training was stopped.

5. The ski resort snow depth inversion method based on optimal feature selection, constrained stacked neural network, and SAR data as described in claim 1, is characterized in that... Step 3 uses a feature selection method based on Spearman correlation coefficient to filter all initial features and obtain the feature combination with the highest correlation to snow depth. The specific operation is as follows: (1) The initial dataset for snow depth inversion prepared in step one is used to standardize all input data for each feature dimension of the dataset using the formula below; (8) Where x is the input data, and μ and σ are the mean and variance of each feature data, respectively; (2) Feature selection based on Spearman correlation coefficient: For the initial feature dataset of snow depth inversion, calculate the Spearman correlation coefficient between each feature and snow depth, and sort the features according to the magnitude of the absolute value of the correlation coefficient to evaluate the degree of association between different features and snow depth, and provide a basis for subsequent feature selection; Root mean square error, mean absolute error, and coefficient of determination were selected as model performance evaluation indicators. The number of features in the snow depth inversion model was gradually increased according to the Spearman correlation coefficient between the features and the snow depth, and then input into the cNN-Stacking algorithm model for training. The changes in the model performance evaluation indicators were observed. When the model evaluation indicators no longer improved or tended to stabilize as the number of features increased, the feature addition process was stopped, and the current feature subset was determined as the optimal feature combination for snow depth inversion.

6. The ski resort snow depth inversion method based on optimal feature selection, constrained stacked neural network, and SAR data as described in claim 5, is characterized in that... In step three, after feature filtering, the optimal feature combination obtained includes the radar local incidence angle (LIA) from Sentinel-1 data and the backscattering coefficients (Sigma) for different polarization modes. VH Gamma VH Beta VH and the corresponding decibel (dB) Sigma VH dB Gamma VH dB Beta VH dB Sigma PS, dB Gamma PS; 2m temperature, snow density, and snow albedo from ERA5 atmospheric reanalysis data; elevation, longitude, slope, and aspect from SRTMDEM 90M data; tree cover, non-tree cover, and non-vegetation cover from MOD44B VCF data; month, date, and number of days with snow cover from meteorological station observations.

Images