A method for reconstructing snow-covered area in mountainous regions based on machine learning

Abstract

This invention discloses a machine learning-based method for reconstructing snow cover area in mountainous areas, comprising: acquiring and preprocessing optical remote sensing data and topographic data; the optical remote sensing data includes MOD09GA data and Landsat8OLI data; the preprocessing includes topographic correction and snow cover identification; using the snow cover area generated from the Landsat8OLI data using a snow cover identification algorithm as the ground truth label, and constructing a sample dataset based on the preprocessed MOD09GA data and topographic data; dividing the sample dataset into a training set and a validation set at a 4:1 ratio to train and validate an XGBoost model; the XGBoost model improves prediction performance by combining multiple decision trees; and using the trained XGBoost model to predict snow cover area on new input data. This method can generate daily 30-meter (m) snow cover area data, which has important research significance and value for mountainous snow cover monitoring, hydrological process simulation, and ecosystem evolution.

Description

Technical Field

[0001] This invention relates to the field of snow cover monitoring, and relates to, but is not limited to, a method for reconstructing snow cover area in mountainous areas based on machine learning. Background Technology

[0002] Snow cover is the most widespread component of the cryosphere, the second largest component, and one of the most dynamic elements on the Earth's surface. It is also the fastest-changing variable in the cryosphere with seasonal variations. Accurate monitoring of the spatiotemporal changes in snow cover is of great significance for research on climate change, water resource management, and disaster prevention.

[0003] Snow cover area is a crucial variable reflecting the spatiotemporal distribution of snow and is also important input data for some hydrological models. There are two main methods for obtaining snow cover area: traditional observation methods, using meteorological stations and snow gauge lines, but these are limited by station distribution and cannot accurately monitor snow distribution and changes in inaccessible mountainous areas. Satellite remote sensing, with its rapid, all-weather, and large spatial scale, has become the preferred method for large-scale snow cover parameter retrieval. It can acquire parameters such as snow cover area over large areas and inaccessible locations, greatly compensating for the shortcomings of traditional observation methods. Remote sensing for snow cover area can be divided into optical remote sensing and microwave remote sensing. However, the lower spatial resolution of microwave remote sensing limits its application in mountainous and other complex terrain areas. Therefore, most research on snow cover area acquisition uses satellite optical remote sensing data. Hyperspectral and high spatiotemporal resolution optical remote sensing sensors can effectively identify the distribution range and temporal trends of snow cover and have been widely used for snow cover area monitoring.

[0004] With the continuous progress and development of remote sensing technology, the algorithms and products for retrieving snow cover area have become mature. The main forms of snow cover area retrieved by remote sensing are binary snow cover area (SCE) and fractional snow cover (FSC). Binary snow cover area uses 0 or 1 to indicate whether the pixel is a snow-covered pixel, while fractional snow cover uses a value range of 0 to 1 to represent the proportion of snow cover within a pixel.

[0005] Snow cover exhibits strong spatiotemporal heterogeneity and is one of the fastest-changing large-scale features on Earth. Traditional monitoring methods are insufficient to effectively monitor large-scale spatiotemporal changes in snow cover and cannot provide sufficient snow cover extent data for watershed-scale hydrological models. With the development of remote sensing technology, an effective solution has been provided for rapidly and accurately obtaining large-scale snow cover. However, existing snow area products, including Moderate-Resolution Imaging Spectroradiometer (MODIS), still have inherent limitations. For example, the MOD10A1 daily product has a relatively coarse spatial resolution (500m) that cannot accurately depict the spatial details of snow cover in mountainous areas. The US Landsat series, which has a better spatial resolution, has a 16-day revisit period that is difficult to capture the temporal characteristics of snowmelt changes. Due to the limitations of sensor performance, high temporal and spatial resolution are often mutually exclusive. When faced with the highly heterogeneous distribution of mottled snow cover in mountainous areas and the modeling and simulation of runoff and snowmelt equivalent in mountainous areas, it is inadequate and cannot simultaneously meet the spatial needs of highly spatially heterogeneous snow cover distribution in mountainous areas and the temporal series continuity needs of hydrological simulation.

[0006] Therefore, developing a high-precision, high-spatiotemporal-resolution snow cover area downscaling inversion algorithm to produce daily 30m snow cover area data is of great research significance and value for snow cover monitoring in mountainous areas, hydrological process simulation, and ecosystem evolution. Summary of the Invention

[0007] In view of this, embodiments of the present invention provide a method for reconstructing snow cover area in mountainous areas based on machine learning, which at least overcomes the limitations of existing technologies in obtaining snow cover area data with high spatiotemporal resolution.

[0008] The specific technical solutions of this invention are as follows:

[0009] This invention provides a method for reconstructing snow cover area in mountainous areas based on machine learning, the method comprising:

[0010] Optical remote sensing data and terrain data are acquired and preprocessed. The optical remote sensing data includes MOD09GA data and Landsat8 OLI data. The preprocessing includes terrain correction and snow cover recognition. The snow area generated by the snow cover recognition algorithm in the Landsat8 OLI data is used as the ground truth label. Based on the preprocessed MOD09GA data and terrain data, a sample dataset is constructed. The sample dataset is divided into a training set and a validation set at a ratio of 4:1 to train and validate the XGBoost model. The XGBoost model improves prediction performance by combining multiple decision trees. The trained XGBoost model is used to predict the snow area on new input data.

[0011] In some embodiments, the MOD09GA data includes seven surface reflectance bands, as well as four bands related to the observation state: sensor zenith angle, sensor azimuth angle, solar zenith angle, and solar azimuth angle; the terrain data includes 30-meter resolution elevation images, slope data, and aspect data.

[0012] The acquisition and preprocessing of optical remote sensing data and topographic data includes: performing topographic correction on MOD09GA data and Landsat8 OLI data; and using a snow identification algorithm combining three remote sensing indices—Normalized Difference Snow Index (NDSI), Normalized Difference Vegetation Index (NDVI), and Normalized Difference Forest Snow Index (NDFSI)—to identify snow cover in the Landsat8 OLI surface reflectance image, thereby obtaining a snow area image.

[0013] In some embodiments, the following terrain correction model is used during the terrain correction process:

[0014] L H (λ)=L I (λ)-a(λ) * (IC-IC H );

[0015]

[0016] In the formula, Z is the solar zenith angle, S is the slope, and φ is the angle of inclination. z The azimuth of the sun. For slope aspect, IC is defined as a function of solar zenith angle, slope, solar azimuth angle, and slope aspect, with a value range from -1 to 1; L H (λ) represents the corrected pixel reflectance value in the λ-band, L I (λ) represents the pixel reflectance value of the λ band before correction, and a(λ) represents the radiance value L from the original image. I The relationship between (λ) and its corresponding IC L I (λ)=a(λ)*IC+b(λ) to obtain, IC H IC is the cosine of the solar zenith angle on the horizontal plane.

[0017] In some embodiments, the snow area generated by the snow recognition algorithm from Landsat8 OLI data is used as the ground truth label, and a sample dataset is constructed based on the preprocessed MOD09GA data and terrain data, including: cropping the MOD09GA image and elevation image based on the range of the Landsat8 OLI surface reflectance image at the same time on the same day.

[0018] The cropped image was resampled to the same 30-meter spatial resolution as the Landsat8 OLI image using nearest neighbor interpolation, and the number of rows and columns of all used images were matched. The 30-meter snow area image obtained by using the snow cover recognition algorithm from the Landsat8 OLI surface reflectance image at the same time on the same day was used as the ground truth label, and cloud pixels were removed using the QA band of Landsat8 OLI. All data unified to the 30-meter spatial resolution were stitched, segmented, filled with missing values, and normalized. The processed data was saved as an npy file.

[0019] In some embodiments, the method further includes: constructing the XGBoost model using the Jupyter Notebook software platform; training the input data separately according to the row and column numbers of Landsat8 OLI during model training; wherein, the positions of different row and column numbers are used to train and establish snow area reconstruction rules applicable to different regions, and the position of a single row and column number is used to set a preset number of effective pixels for training.

[0020] In some embodiments, the snow cover identification algorithm using a combination of three remote sensing indices—NDSI, NDVI, and NDFSI—to identify snow cover in Landsat 8 OLI surface reflectance images includes: in non-forested and forested areas, when NDSI is greater than 0.4 and NIR reflectance in the near-infrared band is greater than 0.11, it is identified as snow cover; in forested areas, when NDSI is between 0 and 0.4, and when NDFSI is greater than 0.4 and NDVI is less than 0.6, it is identified as snow cover.

[0021] In some embodiments, the method further includes: selecting 12 Landsat 8 OLI images that fully cover the study area as a test dataset; using MOD09GA and elevation data as input, reconstructing snow area data through a trained XGBoost model; generating snow area as real data using the surface reflectance of cloudless Landsat 8 OLI images; wherein the cloud cover in the cloudless images is less than a specific value; comparing the reconstructed snow area data with the real data to evaluate the performance of the trained XGBoost model.

[0022] The beneficial effects of the technical solutions provided in the embodiments of the present invention include at least the following:

[0023] In this embodiment of the invention, surface reflectance data and topographic data from MOD09GA and Landsat8 OLI are acquired. After data preprocessing, a sample dataset is created, and a machine learning model based on the XGBoost algorithm is trained to directly reconstruct a 30-meter (m) high-resolution snow cover area. This method uses high temporal and spatial resolution optical remote sensing data and topographic data. Combined with the XGBoost algorithm, it can better learn the complex patterns and relationships in these data, achieving more efficient and accurate prediction and reconstruction of snow cover area. Compared with traditional downscaling methods, the new method has higher accuracy and can achieve a daily 30m snow cover area downscaling effect, improving the ability to monitor snow cover area. Attached Figure Description

[0024] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort, wherein:

[0025] Figure 1 A flowchart illustrating the method for reconstructing snow cover area in mountainous areas based on machine learning, as provided in an embodiment of the present invention.

[0026] Figure 2 A comparison chart of the changes in NDSI calculated from the Landsat 8 OLI surface reflectance before and after terrain correction, provided in an embodiment of the present invention.

[0027] Figure 3 A schematic diagram illustrating the training process of the XGBoost model provided in an embodiment of the present invention;

[0028] Figure 4 A comparison image of the snow accumulation effect at row and column number 135033 provided in this embodiment of the invention, showing the effect of true color synthesis and XGBoost reconstruction. Detailed Implementation

[0029] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. The following embodiments are used to illustrate the present invention, but are not intended to limit the scope of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0030] In the following description, references are made to “some embodiments,” which describe a subset of all possible embodiments. However, it is understood that “some embodiments” may be the same subset or different subsets of all possible embodiments and may be combined with each other without conflict.

[0031] It should be noted that the terms "first, second, and third" used in the embodiments of the present invention are only used to distinguish similar objects and do not represent a specific ordering of objects. It is understood that "first, second, and third" can be interchanged in a specific order or sequence where permitted, so that the embodiments of the present invention described herein can be implemented in an order other than that illustrated or described herein.

[0032] It will be understood by those skilled in the art that, unless otherwise defined, all terms used herein (including technical and scientific terms) have the same meaning as commonly understood by one of ordinary skill in the art to which these embodiments of the invention pertain. It should also be understood that terms such as those defined in general dictionaries should be understood to have the same meaning as in the context of the prior art and should not be interpreted in an idealized or overly formal sense unless specifically defined as herein.

[0033] Machine learning can handle high-dimensional, multi-feature data and uncover implicit relationships between different features within complex data sets. This invention uses the superior Extreme Gradient Boosting (XGBoost) model to reconstruct 30m high-resolution snow cover area data.

[0034] Figure 1 A flowchart illustrating a method for reconstructing snow cover area in mountainous areas based on machine learning, as provided in an embodiment of the present invention, is shown below. Figure 1 As shown, the method includes at least the following steps:

[0035] Step S110: Acquire optical remote sensing data and terrain data and perform preprocessing; the optical remote sensing data includes MOD09GA data and Landsat8 OLI data; the preprocessing includes terrain correction and snow cover identification.

[0036] Here, MOD09GA data is a standard product of MODIS, providing reflectance information of the Earth's surface across multiple wavelengths. It is of great significance for studying the Earth's surface radiation characteristics, vegetation cover, land use change, and climate change. This data can be obtained through NASA's official website or other data providers.

[0037] The specific characteristics of MOD09GA data are as follows: 1) Band range: MOD09GA products include surface reflectance data for bands 1-7, covering the visible and near-infrared spectral range; 2) Spatial resolution: MOD09GA products offer two spatial resolutions: 500 meters and 1000 meters. Reflectance data for bands 1-7 is provided at 500-meter resolution, while additional geographic information (such as observation count and quality assessment level) is provided at 1000-meter resolution; 3) Temporal resolution: MOD09GA products are daily data, providing surface reflectance observations once a day; 4) Data format: MOD09GA data is typically stored in HDF4 (Hierarchical Data Format version 4), a file format used for storing and organizing large-scale data.

[0038] Landsat 8 OLI data, provided by NASA, includes nine bands with a spatial resolution of 30 meters, including a 15-meter panchromatic band. This data is available through the Geospatial Data Cloud website or the U.S. Geological Survey (USGS) website.

[0039] The terrain data used was NASA's SRTM V3 digital elevation model (DEM) data. Elevation data for the study area was extracted, and slope and aspect topographic factors were calculated based on the elevation data. Elevation (DEM), slope, and aspect data served as spatial feature auxiliary data for the algorithm and as important input parameters for terrain correction.

[0040] The preprocessing of the acquired data includes: performing terrain correction on MOD09GA and Landsat8 OLI data, and using terrain features such as DEM, slope, and aspect as spatial feature auxiliary data for the algorithm; and performing snow cover recognition processing on the Landsat8 OLI data. The purpose of terrain correction is to eliminate the influence of terrain effects on the algorithm. By performing geometric transformation on reflectance, the reflectance of pixels is converted to the same reference surface, eliminating the changes in radiance values ​​caused by terrain variations, thus allowing pixels to reflect the true spectral characteristics of ground features. In the snow cover recognition processing, appropriate snow cover recognition algorithms (such as thresholding and classification methods) are used to extract snow cover area from the Landsat8 OLI data.

[0041] In other methods, preprocessing operations also include: masking the MOD09GA band data based on cloud and water masks provided by the quality control QA band to remove cloud and water pixels from the surface reflectance, retain land pixel information under clear sky, and avoid the influence of clouds and water on the algorithm inversion.

[0042] Step S120: Using the snow area generated by the snow recognition algorithm from the Landsat8 OLI data as the ground truth label, construct a sample dataset based on the preprocessed MOD09GA data and terrain data.

[0043] Here, the snow area generated by the Landsat 8 OLI data using a snow area recognition algorithm is used as the ground truth input during model training, and cloud pixels are removed using the quality control QA band. Then, the terrain-corrected MOD09GA data, Landsat 8OLI snow area data, and terrain data (elevation, slope, aspect) are integrated into a sample dataset.

[0044] Step S130: The sample dataset is divided into a training set and a validation set in a 4:1 ratio to train and validate the XGBoost model; the XGBoost model improves prediction performance by combining multiple decision trees.

[0045] Here, XGBoost is an ensemble algorithm based on gradient boosting decision trees, capable of handling classification and regression problems with weak supervision. In snow cover area reconstruction, the XGBoost algorithm combines multiple weak classifiers into a single strong classifier to achieve better classification performance.

[0046] The XGBoost algorithm is as follows: Given a dataset D = {(x ...} containing n samples and m features. i ,y i )}(|D|=n,x i ∈R m ,y i ∈R n The sample prediction results can be expressed as:

[0047]

[0048] In the formula, For the predicted results, K is the number of trees, x i Representative sample, f k The tree structure corresponding to the predicted score, Let T represent the spatial function of the ensemble of K trees, where T is the number of nodes in the trees, s is the node label corresponding to the sample, and w is the number of nodes in the ensemble. s (x) represents the score at node q.

[0049] XGBoost performs feature splitting by continuously adding trees. Each added tree learns a new function, and each decision tree is trained based on the prediction results of the previous round to minimize the value of the objective function. This step-by-step iterative training method gives the model good generalization ability.

[0050] Step S140: Use the trained XGBoost model to predict the snow cover area on the new input data.

[0051] Here, a pre-trained XGBoost model is used to predict new input data, which includes the same feature data as the training sample data. Based on the prediction results, a 30m high-resolution snow cover area reconstruction map of the mountainous area is generated. The reconstructed 30m high-resolution snow cover area can be used in fields such as snow cover monitoring and climate change research. In other implementations, the pre-trained XGBoost model can also be validated using data from a validation set to evaluate the model's predictive performance.

[0052] This invention enables the generation of accurate snow cover area data using high-resolution Landsat 8 OLI data. Furthermore, by combining MOD09GA and DEM data, the XGBoost algorithm is employed to improve the accuracy and efficiency of snow cover area reconstruction. This machine learning-based snow cover area reconstruction method offers advantages such as high efficiency, accuracy, and scalability, and is of significant importance for fields such as snow cover monitoring and climate change research.

[0053] In some embodiments, the MOD09GA data includes seven surface reflectance bands, as well as four bands related to the observation state: sensor zenith angle, sensor azimuth angle, solar zenith angle, and solar azimuth angle; the terrain data includes 30-meter resolution elevation images, slope data, and aspect data.

[0054] The above step S110, "acquiring optical remote sensing data and topographic data and performing preprocessing", includes the following steps: performing topographic correction on MOD09GA data and Landsat8 OLI data; using a snow cover identification algorithm that combines three remote sensing indices, NDSI, NDVI and NDFSI, to identify snow cover on the Landsat8 OLI surface reflectance image and obtain a snow cover area image.

[0055] In mountainous areas, due to the undulating terrain and relatively complex topography, significant variations in slope aspect and gradient result in substantial differences in the irradiance received by pixels on different slope aspects. For example, shaded slopes have weaker irradiance and lower brightness values, leading to lower reflectivity values, while sunny slopes have stronger irradiance and higher brightness values, resulting in higher reflectivity values. This causes the same ground feature to exhibit different spectral information depending on its slope aspect, while different ground features may exhibit the same spectral information on different slope aspects. Therefore, inversion estimation methods that assume a flat surface will result in significant errors when applied to mountainous areas. To eliminate the influence of topographic effects, topographic correction has become an essential step in the preprocessing of remote sensing images in mountainous areas.

[0056] In some embodiments, the following terrain correction model is used during the terrain correction process:

[0057] L H (λ)=L I (λ)-a(λ) * (IC-IC H );

[0058]

[0059] In the formula, Z is the solar zenith angle, S is the slope, and φ is the angle of inclination. z The azimuth of the sun. For slope aspect, IC is defined as a function of solar zenith angle, slope, solar azimuth angle, and slope aspect, with a value range from -1 to 1; L H (λ) represents the corrected pixel reflectance value in the λ-band, L I (λ) represents the pixel reflectance value of the λ band before correction, and a(λ) represents the radiance value L from the original image. I The relationship between (λ) and its corresponding IC L I (λ)=a(λ)*IC+b(λ) to obtain, IC H IC is the cosine of the solar zenith angle on the horizontal plane.

[0060] In the above embodiments, a semi-empirical model based on an empirical rotation model is used to perform terrain correction on the study area imagery. This model introduces an empirical rotation model of a curved surface that does not assume a Lambertian body, transforming the reflectivity of pixels to the same reference surface and eliminating the variation in radiance values ​​caused by terrain changes. The model formula considers the relationship between pixel reflectivity values ​​before and after correction, as well as the influence of the solar incidence angle (IC) on reflectivity. The IC value ranges from -1 to 1 and is determined by the solar zenith angle, slope, solar azimuth angle, and aspect. The slope and aspect are calculated based on the DEM.

[0061] This terrain correction model can effectively eliminate the effects of terrain. Figure 2 This image shows a comparison of the changes in NDSI topographic correction calculated using Landsat 8OLI surface reflectance before and after correction. As shown in Figure 2, a cropped image from December 20, 2020, at row and column number 136033, shows a comparison of the mountainous undulating area before and after correction. Before topographic correction, from... Figure 2 As can be seen, the image is affected by the undulating mountain terrain. Areas on the shady slopes of mountains that should be covered in snow have lower pixel radiance, resulting in lower calculated NDSI values ​​that fail to reflect the true information about the terrain features. After terrain correction, the low radiance values ​​in the shaded areas of the shady slopes are improved; for example, the NDSI in the low-value areas of the northeastern part of the cropped area is higher than before. This demonstrates that terrain-corrected imagery can effectively reduce the mountain effect caused by topographic relief.

[0062] In some embodiments, step S120 above is further implemented through the following steps: based on the range of the Landsat 8 OLI surface reflectance image at the same time on the same day, the MOD09GA image and elevation image are cropped; the cropped image is resampled to the same 30-meter spatial resolution as the Landsat 8 OLI image using the nearest neighbor interpolation method, and the number of rows and columns of all used images are matched to be consistent; the 30-meter snow area image obtained by the snow recognition algorithm from the Landsat 8 OLI surface reflectance image at the same time on the same day is used as the ground truth label, and cloud pixels are removed using the QA band of Landsat 8 OLI; all data unified to the 30-meter spatial resolution are stitched, segmented, filled with missing values, and normalized, and the processed data is saved as an npy file.

[0063] Here, the cropping operation is used to crop feature data such as MODIS and DEM images to the same range as the Landsat8 OLI image. The resampling operation uses the nearest neighbor difference method to interpolate and resample all images to 30m resolution, ensuring that the number of rows and columns in all images matches consistently. The cloud pixel removal operation uses the Landsat8 OLI's quality control QA bands to remove cloud pixels. Then, the data unified to 30m resolution undergoes stitching, segmentation, and missing value imputation operations. Data normalization is performed to improve the efficiency and effectiveness of model training. The processed data is saved in .npy format for convenient use in subsequent model predictions.

[0064] When inputting the model, a 30m resolution snow area image was generated from Landsat 8OLI surface reflectance imagery using a snow cover recognition algorithm as the ground truth label. Each sample should contain multiple features and corresponding labels (i.e., the snow area or snow extent generated at the current time). These features include seven reflectance bands from the MOD09GA data, four bands related to the observation state (sensor zenith angle, sensor azimuth angle, solar zenith angle, and solar azimuth angle), and three terrain features (DEM, slope, and aspect) such as surface reflectance bands and terrain factors, totaling 14 feature data. The process of training the model using the constructed sample dataset is as follows: Figure 3 As shown, 14 feature data of each sample were input into the XGBoost model to reconstruct a 30m high-resolution snow cover area.

[0065] In some embodiments, the method further includes: during model training, training the input data separately according to the row and column numbers of Landsat8 OLI; wherein, the positions of different row and column numbers are used to train and establish snow area reconstruction rules applicable to different regions, and the positions of a single row and column number are used to set a preset number of effective pixels for training.

[0066] Here, this embodiment of the invention implements the XGBoost model on the Jupyter Notebook software platform. This is a snow cover area downscaling model. Preprocessed data is input into the model for training.

[0067] Since different row and column numbers may have different geographical characteristics, data from different Landsat 8OLI rows and columns were trained separately to build snow cover area reconstruction models for different geographical regions, improving the model's applicability and accuracy. Approximately 360 million valid pixels were used for training at a single row and column number, providing a large amount of training data for the model.

[0068] In some embodiments, the snow cover identification algorithm using a combination of three remote sensing indices—NDSI, NDVI, and NDFSI—to identify snow cover in Landsat 8 OLI surface reflectance images includes: in non-forested and forested areas, when NDSI is greater than 0.4 and NIR reflectance in the near-infrared band is greater than 0.11, it is identified as snow cover; in forested areas, when NDSI is between 0 and 0.4, and when NDFSI is greater than 0.4 and NDVI is less than 0.6, it is identified as snow cover.

[0069] Currently, the widely used SNOMAP snow cover identification algorithm mainly relies on the NDSI threshold to distinguish between snow-covered and non-snow-covered land features. However, research shows that the snow cover discrimination threshold varies under different vegetation cover conditions. This invention employs a multi-remote sensing index dynamic identification algorithm combining NDSI, NDVI, and NDFSI to more accurately identify snow cover on different land cover types, thus improving the accuracy of snow cover identification.

[0070] In some embodiments, the method further includes: selecting 12 Landsat 8 OLI images that fully cover the study area as a test dataset; using MOD09GA and elevation data as input, reconstructing snow area data through a trained XGBoost model; generating snow area as real data using the surface reflectance of cloudless Landsat 8 OLI images; wherein the cloud cover in the cloudless images is less than a specific value; comparing the reconstructed snow area data with the real data to evaluate the performance of the trained XGBoost model.

[0071] here, Figure 4 The comparison of the reconstructed snow cover area at row and column number 135033 is shown. Figure 4 As can be seen, the terrain at row and column number 135033 is complex, with many fine, patchy snow accumulations, and the snow distribution is quite varied. The snow area reconstructed by the XGBoost algorithm can better represent the snow distribution.

[0072] This invention effectively mitigates the terrain effects caused by mountainous terrain through terrain correction. Based on MODIS and Landsat8 OLI data, as well as terrain data, the XGBoost algorithm is used to directly reconstruct snow cover area with high spatiotemporal resolution. The XGBoost model is chosen as the core algorithm because, compared to traditional random forest algorithms, the superior performance of the Extreme Gradient Boosting (XGBoost) model significantly improves efficiency when processing large-scale data. The model's generalization ability is also significantly enhanced, enabling it to handle higher-dimensional data with more features and to uncover implicit relationships between different features within complex data. This is crucial for snow cover area reconstruction in terrain-complex regions.

[0073] The algorithm of this invention can effectively improve the accuracy of snow cover area reconstruction, achieving good reconstruction results with high accuracy of reconstructed snow cover pixels and a low probability of misclassified snow cover pixels. The reconstructed snow cover area can clearly depict the spatial details of snow cover distribution, displaying a more accurate snow cover distribution. Furthermore, it exhibits high stability across different indices such as OA, PA, and Kappa coefficient values. By combining machine learning with remote sensing data, high spatiotemporal resolution (30m) snow cover area reconstruction can be achieved daily in mountainous areas.

[0074] It should be understood that the phrase "one embodiment" or "an embodiment" throughout the specification means that a specific feature, structure, or characteristic related to the embodiment is included in at least one embodiment of the invention. Therefore, "in one embodiment" or "in an embodiment" appearing throughout the specification does not necessarily refer to the same embodiment. Furthermore, these specific features, structures, or characteristics can be combined in any suitable manner in one or more embodiments. It should be understood that in the various embodiments of the invention, the sequence numbers of the above-described processes do not imply a sequential order of execution; the execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the invention. The sequence numbers of the above-described embodiments of the invention are merely descriptive and do not represent the superiority or inferiority of the embodiments.

[0075] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0076] In the several embodiments provided by this invention, it should be understood that the disclosed methods can be implemented in other ways. The methods disclosed in the several method embodiments provided by this invention can be arbitrarily combined without conflict to obtain new method embodiments. The features disclosed in the several method embodiments provided by this invention can be arbitrarily combined without conflict to obtain new method embodiments.

[0077] The above description is merely an embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for reconstructing snow cover area in mountainous areas based on machine learning, characterized in that, include: Acquire optical remote sensing data and terrain data and perform preprocessing; The optical remote sensing data includes MOD09GA data with a spatial resolution of 500 meters and Landsat8 OLI data with a spatial resolution of 30 meters; the preprocessing includes terrain correction and snow cover identification; the terrain data includes 30-meter resolution elevation imagery, slope data, and aspect data. The snow area generated by the snow recognition algorithm on the preprocessed Landsat8 OLI data was used as the ground truth label. The preprocessed MOD09GA data and terrain data were resampled to the same 30-meter spatial resolution as the Landsat 8 OLI data to construct a sample dataset. The sample dataset is divided into a training set and a validation set in a 4:1 ratio for training and validation of the XGBoost model; the XGBoost model improves prediction performance by combining multiple decision trees. The trained XGBoost model is used to predict the snow cover area of ​​the new input data, and a 30-meter high-resolution snow cover area reconstruction map of the mountainous area is generated based on the prediction results. The new input data includes the same feature data as the training sample data. The feature data includes seven reflectance bands, four bands related to the observation state, and three terrain features from the MOD09GA data.

2. The method according to claim 1, characterized in that, The MOD09GA data includes 7 surface reflectance bands, as well as 4 bands related to the observation status: sensor zenith angle, sensor azimuth angle, solar zenith angle, and solar azimuth angle. The acquisition and preprocessing of optical remote sensing data and terrain data includes: Terrain correction was performed on MOD09GA and Landsat8 OLI data; A snow cover identification algorithm combining three remote sensing indices—NDSI, NDVI, and NDFSI—was used to identify snow cover in Landsat 8 OLI surface reflectance images, yielding snow cover area images.

3. The method according to claim 2, characterized in that, The following terrain correction model is used in the terrain correction process: ； ； In the formula, Z The zenith angle of the sun. S For slope, The azimuth of the sun. Slope direction, IC Defined as a function of solar zenith angle, slope, solar azimuth angle, and slope aspect, its value ranges from -1 to 1; For correction Pixel reflectance values ​​in the band, Before correction Pixel reflectance values ​​in the band, The original image radiance value Its corresponding IC Relationship Obtain For horizontal plane IC That is, the cosine of the solar zenith angle.

4. The method according to any one of claims 1 to 3, characterized in that, Using snow cover area generated from Landsat 8 OLI data using a snow cover recognition algorithm as the ground truth label, and based on preprocessed MOD09GA data and terrain data, a sample dataset was constructed, including: Based on the range of the Landsat 8 OLI surface reflectance image at the same time on the same day, the MOD09GA image and elevation image are cropped. The cropped image was resampled to the same 30-meter spatial resolution as the Landsat 8 OLI image using nearest neighbor interpolation, and the number of rows and columns of all images used were matched. The 30-meter snow area image obtained by using the snow recognition algorithm from the Landsat 8 OLI surface reflectance image at the same time on the same day is used as the ground truth label, and cloud pixels are removed using the QA band of Landsat 8 OLI. All data, unified to a spatial resolution of 30 meters, are spliced, segmented, missing value filled, and normalized. The processed data is then saved as an npy file.

5. The method according to any one of claims 1 to 3, characterized in that, The method further includes: The XGBoost model was constructed using the Jupyter Notebook software platform. During model training, the input data is trained separately according to the row and column numbers of Landsat8 OLI; the positions of different row and column numbers are used to train and establish snow area reconstruction rules applicable to different regions, while the positions of a single row and column number are used to set a preset number of effective pixels for training.

6. The method according to claim 2, characterized in that, The snow cover identification algorithm, which combines three remote sensing indices—Normalized Snow Index (NDSI), Normalized Vegetation Index (NDVI), and Normalized Forest Snow Index (NDFSI)—to identify snow cover in Landsat 8OLI surface reflectance images, includes: In both non-forested and forested areas, when the NDSI is greater than 0.4 and the NIR reflectance in the near-infrared band is greater than 0.11, it is identified as snow. In forested areas, snow is identified when NDSI is between 0 and 0.4, and when NDFSI is greater than 0.4 and NDVI is less than 0.

6.

7. The method according to any one of claims 1 to 3, characterized in that, The method further includes: Twelve scenes of Landsat 8 OLI data covering the entire study area were selected as the test dataset. Using MOD09GA and elevation data as input, snow cover area data is reconstructed through a trained XGBoost model. Snow cover area is generated using the Landsat 8 OLI surface reflectance of cloudless imagery as the actual data; the cloud cover in the cloudless imagery is less than a specific value. The reconstructed snow area data is compared with the real data to evaluate the performance of the trained XGBoost model.

Images