A sub-pixel snow filling method based on medium-resolution remote sensing data

Abstract

The present application relates to the technical field of remote sensing information extraction and snow monitoring, and particularly relates to a sub-pixel snow mapping method based on medium-resolution remote sensing data, comprising the following steps: obtaining medium-resolution remote sensing images, a digital elevation model, and ground meteorological station snow observation data, and preprocessing the remote sensing images to obtain surface reflectivity; identifying cloud-covered pixels, filling in the cloud-free dataset by temporal and spatial interpolation and ground object spectral similarity matching; selecting snow, vegetation, and water as spectral endmembers and extracting features; constructing a linear mixing model, and using the least square method to retrieve the snow proportion in the pixel to obtain preliminary data; and generating a high temporal and spatial accuracy snow distribution map through meteorological station data verification, digital elevation model terrain correction, and model parameter optimization, which solves the problems of mixed pixels, cloud coverage, terrain errors, and the lack of accuracy verification.

Description

Technical Field

[0001] This invention relates to the field of remote sensing information extraction and snow cover monitoring technology, specifically to a sub-pixel snow cover mapping method based on medium-resolution remote sensing data. Background Technology

[0002] Snow cover, as one of the core reserves of global freshwater resources, is also a key factor in regulating surface energy exchange, influencing regional climate evolution, and the hydrological cycle. It holds irreplaceable strategic significance in areas such as water resource allocation, climate simulation, snow disaster early warning, and ice and snow industry planning. Traditional snow cover monitoring relies on ground-based stations, which, while providing precise data at single points (such as snow depth and snow quality), are limited by station density, covering only localized areas and failing to achieve large-scale, spatially continuous snow cover distribution monitoring. This is particularly challenging for monitoring the dynamic snow cover in complex terrain areas such as mountains and plateaus. With the development of remote sensing technology, its advantages of "large-area, periodic, and real-time" monitoring have become the core means to overcome this limitation. Among these, medium-resolution remote sensing data (such as MODIS and AVHRR series), with their high temporal resolution of daily or 8-day intervals, can promptly capture short-term dynamic changes in snow melting and accumulation, becoming the mainstream data source for large-scale snow cover monitoring.

[0003] However, existing medium-resolution remote sensing snow cover monitoring technologies still face multiple technical bottlenecks: First, medium-resolution data (spatial resolution 500m-1km) is prone to the "mixed pixel" problem—a single pixel often contains multiple land features such as snow, vegetation, and water bodies simultaneously. Traditional binary classification algorithms (which only determine "snow-covered pixels" or "non-snow-covered pixels") cannot accurately invert the actual proportion of snow within a pixel, resulting in a significant limitation on the accuracy of snow cover monitoring. Second, medium-resolution remote sensing data is significantly affected by cloud cover. Snow-covered areas (mostly high-latitude or high-altitude regions) experience frequent cloud cover in winter. Existing cloud removal methods mostly use single temporal interpolation or spatial replacement, resulting in poor spectral consistency between the filled area and surrounding land features, making it difficult to form a complete and accurate dataset of cloudless snow reflectance. Third, the undulations of complex terrain (such as mountains and hills) can lead to distortion of surface reflectance. Although there are relevant terrain correction techniques, such as the one authorized by CN116862798B, these techniques cannot completely eliminate the distortion. The proposed method for correcting topographic shading effects in optical satellite remote sensing images constructs a model by introducing shadow intensity factor (SIF), vegetation index factor (VIF), and band adjustment factor (BAF) to correct reflectance distortion caused by topographic shading. However, this method only addresses topographic shading correction for general optical satellite images and does not consider the specific characteristics of snow cover monitoring. It fails to deeply integrate topographic correction with core processes such as snow endmember extraction and mixed pixel decomposition, thus failing to address the accuracy issues of snow cover inversion under complex terrain. Fourth, existing sub-pixel snow cover mapping techniques lack a closed loop of "ground observation - model verification - parameter optimization," relying heavily on theoretical model assumptions. They do not verify the inversion results using snow observation data from ground meteorological stations (such as snow depth and snow condition records) or dynamically adjust model parameters based on verification errors, making it difficult to guarantee inversion accuracy. Fifth, some methods lack scientific rigor in endmember selection, failing to clearly define the spectral characteristics of endmembers such as snow cover, vegetation, and water bodies based on the target area's land cover features, and failing to verify endmember purity, further amplifying the errors in mixed pixel decomposition. In summary, existing technologies cannot simultaneously solve problems such as mixed pixel interference, cloud cover effects, complex terrain errors, and lack of accuracy verification in medium-resolution remote sensing data. They are unable to generate sub-pixel snow distribution maps with both high spatiotemporal resolution and high accuracy, and cannot meet the demand for accurate snow monitoring data in scenarios such as regional water resource management and ice and snow industry planning. Summary of the Invention

[0004] The purpose of this invention is to provide a sub-pixel snow cover mapping method based on medium-resolution remote sensing data, so as to solve the problems of mixed pixel interference, cloud cover influence, complex terrain error and lack of accuracy verification in medium-resolution remote sensing data.

[0005] To achieve the above objectives, the following technical solution is adopted.

[0006] A sub-pixel snow cover mapping method based on medium-resolution remote sensing data includes the following steps:

[0007] Step S1: Acquire medium-resolution remote sensing image data, digital elevation model data, and snow cover observation data from ground meteorological stations for the target area, and perform radiometric calibration and atmospheric correction preprocessing on the medium-resolution remote sensing image data to obtain surface reflectance data.

[0008] Step S2: Use cloud detection algorithm to identify cloud-covered pixels in surface reflectance data, and use spatiotemporal interpolation and ground object spectral similarity matching to fill in the areas where cloud-covered pixels are located, generating a complete surface reflectance dataset without cloud interference.

[0009] Step S3: Based on the land cover features of the target area, select snow cover, vegetation and water bodies as spectral endmembers, and extract the spectral curves of pure pixels from the complete surface reflectance dataset without cloud interference to determine the spectral features of each spectral endmember.

[0010] Step S4: Construct a linear mixing model. The reflectance of each pixel in the complete surface reflectance dataset without cloud interference is represented as a combination of the reflectance of each spectral endmember and its proportion. The least squares method is used to solve the linear mixing model to inversely determine the proportion of snow endmembers in each pixel, thus obtaining preliminary sub-pixel snow cover data.

[0011] Step S5: Use snow observation data from ground meteorological stations to verify the accuracy of the preliminary sub-pixel snow cover data, and combine it with digital elevation model data to perform terrain correction to correct reflectivity errors caused by terrain undulations. Improve the snow cover inversion accuracy by optimizing the parameters of the linear hybrid model, and finally generate a verified and corrected high spatiotemporal accuracy sub-pixel snow distribution map.

[0012] Optionally, in step S1, the medium-resolution remote sensing image data used is surface reflectance product data with a daily or 8-day temporal resolution and a spatial resolution of 500 meters or 1000 meters. The digital elevation model data is used to provide elevation, slope and aspect information of the target area. The snow cover observation data from the ground meteorological station includes snow depth and snow condition records, which are used for model verification and accuracy evaluation.

[0013] Optionally, the preprocessing steps for medium-resolution remote sensing image data may also include performing topographic correction on the radiometrically calibrated and atmospherically corrected surface reflectance data using digital elevation model data, in order to eliminate the influence of topographic shadows on surface reflectance.

[0014] Optionally, in the step of filling in the area where cloud-covered pixels are located, the spatiotemporal interpolation method is based on the historical reflectance data of the target pixels in the cloudless phase to perform time series interpolation or extrapolation, while the land cover spectral similarity matching method is to find cloudless pixels with similar terrain environment and land cover type to cloud-covered pixels in the spatial neighborhood and replace them with their spectral values. By combining these two methods to seamlessly fill in the cloud-covered area, a complete land reflectance dataset without cloud interference is generated.

[0015] Optionally, the step of selecting snow cover, vegetation, and water bodies as spectral endmembers is determined based on the main land cover type of the target area. The step of extracting the spectral curve of pure pixels from the complete surface reflectance dataset without cloud interference is to select a large, continuous, and spectrally uniform area on the image as the endmember sample selection area, and calculate the average reflectance of all pixels in the selected area as the representative spectral curve of the spectral endmember.

[0016] Optionally, after extracting the spectral curves of pure pixels, a spectral angle mapping algorithm is used to calculate the similarity between the spectrum of each pixel in the target area and the spectral curves of each extracted spectral endmember. This is used to verify the spectral purity of the selected endmember samples and to ensure that the endmember spectra used in the linear mixture model can accurately represent the main land cover types in the target area.

[0017] Optionally, in the step of solving the linear mixing model using the least squares method to invert the proportion of snow end-members, the linear mixing model assumes that the reflectance of each pixel is a linear combination of the reflectance of its constituent snow end-members, vegetation end-members, and water end-members, weighted by the proportion of their area within the pixel. By solving this system of linear equations using the least squares method, the proportion of each end-member is obtained that minimizes the error between the model's calculated value and the actual reflectance value of the pixel. The proportion of the snow end-member is the preliminary sub-pixel snow cover data.

[0018] Optionally, the step of using snow cover observation data from ground meteorological stations to verify accuracy is to compare the snow cover observation records at the meteorological station location with the sub-pixel snow cover data at the same location point retrieved by the model, calculate the coefficient of determination R² to quantify the retrieval accuracy, and adjust the parameters of the linear mixture model or the endmember spectral library according to the verification results until the coefficient of determination R² reaches or exceeds the preset threshold.

[0019] Optionally, the step of combining digital elevation model data to perform terrain correction to correct reflectivity errors caused by terrain undulation is to consider the influence of slope and aspect on the intensity of solar radiation received by the surface, establish a correction model using terrain parameters calculated from digital elevation model data, and correct the reflectivity values ​​of pixels significantly affected by terrain in the complete surface reflectivity dataset without cloud interference, thereby reducing the interference of terrain shadows on the accuracy of mixed pixel decomposition and snow cover inversion.

[0020] Optionally, after generating a verified and calibrated high spatiotemporal precision sub-pixel snow distribution map, the sub-pixel snow distribution map is also output as a raster format map product for analyzing the snow area change trend of the target area at different time scales, or for overlay analysis with a specific geographic boundary range to evaluate the snow resource distribution characteristics within that range and its supporting role in water resource management, ice and snow industry planning, or disaster prevention and mitigation decision-making.

[0021] Compared with the prior art, the present invention has the following beneficial effects:

[0022] This application constructs a complete technical process of "data preprocessing - cloud interference removal - endmember extraction and mixed pixel decomposition - accuracy verification and terrain correction - high spatiotemporal accuracy mapping," which accurately solves the problems of mixed pixel interference, cloud cover effects, complex terrain errors, and lack of accuracy verification in medium-resolution remote sensing data. Step S1 performs radiometric calibration and atmospheric correction preprocessing on medium-resolution remote sensing image data to obtain high-quality surface reflectance data, laying the data foundation for subsequent snow cover information extraction. Step S2 innovatively combines spatiotemporal interpolation and ground object spectral similarity matching to process cloud-covered pixels, which significantly improves the integrity and spectral consistency of cloudless data compared to existing single cloud removal methods, avoiding the loss of snow cover information caused by cloud interference. Step S3... Based on the land cover characteristics of the target area, snow cover, vegetation, and water bodies are selected as spectral endmembers, and pure pixel spectral curves are extracted to ensure the representativeness of the endmembers and reduce the decomposition error of mixed pixels from the source. Step S4 uses a linear mixture model and least squares method to invert the preliminary sub-pixel snow cover rate, specifically overcoming the limitations of mixed pixels in medium-resolution data and solving the problem of insufficient accuracy of traditional binary classification. Step S5 combines snow cover observation data from ground meteorological stations for accuracy verification, and uses digital elevation model data for terrain correction. Furthermore, by optimizing model parameters, a closed loop of "data-model-verification-optimization" is formed, effectively correcting the reflectivity error caused by terrain undulations. Finally, a high spatiotemporal accuracy sub-pixel snow cover distribution map is generated, perfectly meeting the triple requirements of "high temporal resolution + high spatial accuracy + high inversion reliability" for regional snow cover monitoring in the background technology.

[0023] Clearly define the temporal and spatial resolution of the medium-resolution remote sensing image data, as well as the specific uses of the digital elevation model data and ground meteorological station data, to ensure the data source is targeted and provide data assurance for mapping accuracy; add DEM topographic correction in step S1 preprocessing to eliminate the influence of topographic shadows on surface reflectivity in advance, avoiding interference from initial data errors in subsequent processes; refine spatiotemporal interpolation (interpolation of historical data from cloudless phases / The operation logic of extrapolation and spectral matching of ground features (spatial neighborhood similar pixel replacement) improves the operability and data consistency of cloud filling; it clarifies the basis for endmember selection and the method of pure pixel extraction, and verifies the purity of endmembers through spectral angle mapping algorithm to ensure endmember accuracy and fundamentally reduce decomposition errors; it refines the solution logic of the linear mixture model to make the snow endmember ratio inversion clearer and more reliable; it quantifies the accuracy verification standard to form measurable accuracy control indicators to ensure the credibility of the inversion results; it constructs a terrain correction model based on slope and aspect to specifically reduce the interference of complex terrain on snow inversion; and it outputs the mapping results as raster products and connects them with applications such as water resource management and ice and snow industry planning to realize the implementation of technical achievements in real scenarios, greatly improve practical value, and fully meet the goal of multi-scenario needs of snow monitoring services in the background technology. Attached Figure Description

[0024] Figure 1 This is a schematic flowchart illustrating the steps of an embodiment of a sub-pixel snow cover mapping method based on medium-resolution remote sensing data according to the present invention. Detailed Implementation

[0025] The present invention will now be described in detail with reference to the accompanying drawings and embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other.

[0026] The following detailed description is exemplary and intended to provide further detailed explanation of the invention. Unless otherwise specified, all technical terms used in this invention have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains. The terminology used in this invention is for the purpose of describing particular embodiments only and is not intended to limit the scope of exemplary embodiments according to the invention.

[0027] like Figure 1 As shown, a sub-pixel snow cover mapping method based on medium-resolution remote sensing data is as follows:

[0028] First, the target area and data source are determined. The target area must meet the requirements of diverse terrain and snow cover characteristics. Such areas include open plains that can provide uniform snow cover scenarios, as well as undulating mountains and hills that can verify the algorithm's adaptability to complex terrains. Furthermore, the long duration and significant dynamic changes of winter snow cover meet the requirements of sub-pixel snow cover mapping for the research object. Medium-resolution remote sensing image data should preferably be the MOD09GA surface reflectance product. This product can be obtained through publicly available remote sensing data sharing platforms. When acquiring the data, image data from the target area for the most recent three complete snow cover periods should be selected. Within each snow cover period, at least 8-10 cloudless or low-cloud phases should be selected to ensure a complete capture of the dynamic changes in snow cover from accumulation, stabilization, to melting. The data download format can be HDF-EOS, with the corresponding metadata file downloaded simultaneously. The metadata includes information such as solar zenith angle and solar azimuth angle, which will be used in subsequent atmospheric and topographic correction stages.

[0029] Digital elevation model (DEM) data with a resolution of at least 30m should be selected. This data accurately reflects the elevation, slope, and aspect characteristics of the target area. It is essential to ensure that the data completely covers the target area. If a small number of data gaps exist, they can be filled using the neighborhood fill function in the spatial analysis toolset of the GIS software. A 3×3 pixel neighborhood range can be selected during filling to ensure that the filled data still accurately reflects the local terrain undulations. The acquired DEM data needs to undergo coordinate system transformation to ensure consistency with the medium-resolution remote sensing image data. Both should use the WGS84 coordinate system and UTM projection to avoid positional deviations during subsequent spatial matching.

[0030] Snow cover observation data from ground-based meteorological stations can be obtained by applying to the meteorological data sharing platform of the target area. Select 15-20 ground-based meteorological stations evenly distributed within and around the target area, ensuring that the service radius of each station does not exceed 50 km to guarantee the spatial representativeness of the observation data. The observation data must include daily records of snow depth, snow type, and observation time during the snow cover period. The observation time should be as consistent as possible with the imaging time of the medium-resolution remote sensing image, with a time difference not exceeding 3 days, to minimize the interference of time differences on subsequent accuracy verification results.

[0031] After data acquisition, the medium-resolution remote sensing image data undergoes preprocessing. Preprocessing includes radiometric calibration, atmospheric correction, and topographic correction, all of which can be performed using ENVI series software. During radiometric calibration, the built-in radiometric calibration coefficients can be accessed through the radiometric correction tools in the software's menu bar to convert the image's DN values ​​to apparent reflectance. During calibration, all reflectance bands requiring processing must be selected, ensuring coverage of key bands such as blue, green, red, and near-infrared light. These bands are crucial for distinguishing snow cover from other ground features. After calibration, the data units are standardized to dimensionless, with values ​​ranging from 0 to 1.

[0032] Atmospheric correction can be performed using the FLAASH atmospheric correction module. Correction parameter settings must be combined with the climatic characteristics of the target area during the snow cover period. The atmospheric model can be a mid-latitude winter model, suitable for mid-latitude winter atmospheric conditions. The aerosol type can be a continental type, suitable for non-coastal land areas. Visibility parameters should prioritize the average visibility observed by ground meteorological stations in the target area during the same period. If specific observation values ​​are unavailable, the implementation report should indicate the use of software default values, and subsequent sensitivity analysis should be conducted to verify the impact of these default values ​​on the correction results. Parameters such as sensor height, imaging time, solar zenith angle, and solar azimuth angle are extracted from the metadata file of the remote sensing image and accurately entered to ensure the accuracy of the correction process. After atmospheric correction, the image data is converted from apparent reflectance to true surface reflectance, eliminating the influence of atmospheric scattering, absorption, and other interference factors.

[0033] The terrain correction process requires the integration of digital elevation model (DEM) data. Slope and aspect information for the target area can be calculated using GIS software's 3D analysis toolset. During the calculation, the output coordinate system must be consistent with the medium-resolution remote sensing image data to ensure spatial matching between terrain parameters and image data. Subsequently, based on slope and aspect data and solar parameters in the metadata, terrain-shadowed areas in the image are identified. The reflectance of pixels affected by shadows is adjusted to eliminate reflectance distortion caused by terrain undulations, providing high-quality surface reflectance data for subsequent cloud detection and endmember extraction.

[0034] Next, cloud cover pixel identification and filling processing is performed. Cloud detection employs a comprehensive algorithm based on reflectance thresholds and vegetation indices, combined with the spectral characteristics of medium-resolution remote sensing imagery. In the processing software, the blue and near-infrared bands of the surface reflectance data are first extracted. When the blue band reflectance is high and the near-infrared band reflectance is low, it is initially identified as a cloud-covered pixel. This is because clouds have a strong reflective ability for blue light and a strong absorption ability for near-infrared light. Simultaneously, the normalized vegetation index (NVI) is calculated. When the NVI is low and the blue band reflectance is high, it is also identified as a cloud-covered pixel. This is because cloud-covered areas typically lack vegetation cover, resulting in a low NVI value. To reduce the false positive rate, a cloud mask product can be introduced as an auxiliary reference. This product can be acquired synchronously with the medium-resolution remote sensing image data and contains cloud cover identification information for each pixel. After spatially registering the cloud mask product with the surface reflectance data, the cloud areas identified by the two are compared. For pixels that the algorithm identifies as clouds but the mask product does not, visual confirmation is performed using the software zoom-in function combined with the surrounding ground features. Similarly, for pixels that the algorithm identifies as non-clouds but the mask product does, visual verification is also performed. Finally, an accurate cloud cover mask map is generated, providing a clear target area for subsequent cloud filling.

[0035] The filling of cloud-covered areas employs a combination of spatiotemporal interpolation and ground object spectral similarity matching. For spatiotemporal interpolation, a time-series surface reflectance database of the target area must first be constructed. This involves collecting surface reflectance data from cloudless phases within the same orbit and sensor as the current cloud-covered image over the past three snow cover periods. These cloudless phase data must undergo a preprocessing procedure identical to that of the current image to avoid reflectance deviations caused by preprocessing differences. For cloud-covered pixels in the current image, if the current phase falls between two cloudless phases, linear interpolation can be used to calculate the reflectance of that pixel. If the current phase is at the beginning or end of the time series, linear extrapolation can be used to calculate the reflectance. The extrapolation time span should not exceed 15 days to control extrapolation errors within a reasonable range.

[0036] The spectral similarity matching method is used to correct the filling error of spatiotemporal interpolation. The spatial neighborhood is determined based on the topographic features of the target area. In mountainous areas with high heterogeneity, the neighborhood radius can be 5 pixels; in plains areas with high homogeneity, the neighborhood radius can be 10 pixels. Within this neighborhood, non-cloud pixels are screened. The topographic environment (slope, aspect differences) and land cover type (determined by normalized vegetation index) of these candidate pixels are compared with those of cloud-covered pixels. Candidate pixels with similar topographic environment and consistent land cover type to cloud-covered pixels are selected. The reflectance values ​​of each band of these candidate pixels are extracted, and the arithmetic mean method is used to calculate the average reflectance of each band as the final filling reflectance for cloud-covered pixels. If no candidate pixels meeting the criteria are found within the neighborhood, the neighborhood radius can be expanded to 15 pixels for rescreening. If no candidate pixels meeting the criteria are found after expansion, the preliminary filling reflectance obtained by spatiotemporal interpolation is used as the final filling reflectance. By combining the two methods, seamless filling of cloud-covered areas is achieved, generating a complete surface reflectance dataset free from cloud interference.

[0037] Subsequently, spectral endmembers were selected, representative spectral curves were extracted, and purity was verified. Based on the main land cover types in the target area, snow cover, vegetation, and water bodies were selected as spectral endmembers. These three land cover types account for the largest area in the target area, and their spectral characteristics are significantly different, which can effectively support the decomposition of mixed pixels. The selection of endmember sample areas can be completed using the region of interest tool in remote sensing processing software. Each endmember can select 3-5 sample areas at different locations to avoid the random errors of a single selection area. For the snow end-member sample selection, choose a large, continuous snow-covered area within the target region that is free of vegetation cover and water bodies. Examples include continuous snow-covered areas on open plains in winter or pure snow-covered areas on gentle mountain slopes. The selected area should be at least 50 pixels, and the near-infrared reflectance of the pixels within the selected area must match the spectral characteristics of pure snow. For the vegetation end-member sample selection, choose natural grasslands or farmland outside of snow cover within the target region. The selected area should be at least 50 pixels, and the normalized difference vegetation index (NDI) of the pixels within the selected area must match the characteristics of well-grown vegetation. For the water end-member sample selection, choose still or slow-flowing water bodies such as lakes and rivers within the target region. The selected area should be at least 30 pixels, and the near-infrared reflectance of the pixels within the selected area must match the spectral characteristics of pure water. After the sample selection is completed, a visual inspection using the software's zoom-in function is required to ensure that the selected area is free of cloud cover, terrain shadows, and other features.

[0038] When extracting the representative spectral curve of the endmember, for all sample selection areas of each endmember, extract the reflectance data of each band of each pixel in the selection area, and use the arithmetic mean method to calculate the average reflectance of each band to obtain the representative spectral curve of the endmember. After the calculation is completed, the representative spectral curves of the three endmembers can be exported as a text format file for easy use in subsequent linear mixing models.

[0039] The spectral purity of endmember samples was verified using a spectral angle mapping algorithm. This algorithm determines similarity by calculating the angle between the target pixel's spectrum and the representative spectral curve of the endmember; the smaller the angle, the higher the similarity and the better the endmember purity. The verification process can be performed in remote sensing software using a spectral angle matching tool. A complete surface reflectance dataset free from cloud interference and representative spectral curves of three endmembers are imported. A spectral angle threshold is set, which can be determined through multiple experiments, selecting the value that effectively distinguishes the spectra of the three types of land cover with the lowest misclassification rate. After running the spectral angle mapping algorithm, the proportion of pixels with spectral angles less than the threshold within each endmember sample selection area is counted out of the total number of pixels in the selection area. If the proportion is greater than 90%, the spectral purity of that endmember is considered to meet the standard; if the proportion is less than 90%, the sample selection area for that endmember is reselected, and the process of selecting the sample area and extracting the representative spectral curve is repeated until the spectral purity of all endmembers meets the standard.

[0040] Next, a linear mixture model is constructed, and the snow endmember proportion is solved using the least squares method. The core assumption of the linear mixture model is that the reflectance of each pixel in the cloud-free, complete surface reflectance dataset is a linear combination of the reflectances of the three endmembers (snow, vegetation, and water) weighted by their respective area proportions within the pixel. The model can be solved using a script in programming software. First, the cloud-free, complete surface reflectance dataset and representative spectral curves of the three endmembers are imported. The reflectance of each band for each pixel is extracted to form a pixel reflectance vector. This vector, along with the endmember spectral matrix, is substituted into the linear mixture model to construct an overdetermined system of equations. The least squares method is used to solve this system of equations, aiming to minimize the sum of squares between the model-calculated reflectance and the actual reflectance of each pixel in each band. Constraints are introduced during the solution process: the sum of the area proportions of the three endmembers must be 1, and each proportion value must be non-negative. This can be achieved using a constraint solution function in the programming software, which supports linear equality constraints and variable boundary constraints. The snow endmember ratio of each pixel is obtained by solving the problem, which is the preliminary sub-pixel snow cover data. This data is saved as a commonly used remote sensing data format file. The value of each pixel in the file is the preliminary snow cover, which ranges from 0 to 1, where 0 represents no snow and 1 represents a pure snow pixel.

[0041] After the preliminary sub-pixel snow cover data is generated, accuracy verification is required using snow cover observation data from ground meteorological stations. Topographic correction is then performed in conjunction with digital elevation model (DEM) data, and model parameters are optimized to improve inversion accuracy. For accuracy verification, the ground meteorological station observation data is first organized, and the latitude and longitude coordinates of each station are extracted. Using the "extract values ​​to points" tool in GIS software, the snow end-member ratio of the corresponding coordinate pixels in the preliminary sub-pixel snow cover data is obtained. These ratio values ​​are used as model inversion values ​​and recorded as an array (containing n values, where n is the number of meteorological stations participating in the verification, and each value corresponds to the snow cover inversion result for one meteorological station location). Simultaneously, the snow cover observation records for each meteorological station during the same period are converted into binary snow presence / absence data: a snow depth greater than 0 is recorded as 1 (indicating snow cover at that location), and a snow depth equal to 0 is recorded as 0 (indicating no snow cover at that location), forming a ground observation value array (also containing n values, corresponding one-to-one with the values ​​in the inversion value array to ensure that the inversion value of each meteorological station matches the observation value).

[0042] After data processing, linear regression analysis can be performed in programming software to calculate the coefficient of determination R². The specific steps are as follows: First, use the linear fitting function in the programming software to perform a first-order linear fit between the inverted value array (independent variable, representing the snow cover inverted by the model) and the observed value array (dependent variable, representing the snow cover state observed on the ground), obtaining the slope and intercept of the fitted line, and then generating the fitted value array corresponding to each inverted value; subsequently, calculate the average of the ground observed values, i.e., sum all the values ​​in the observed value array and divide by n; next, calculate the total sum of squares, ... The regression sum of squares and residual sum of squares are calculated as follows: the total sum of squares reflects the overall variability of the ground observations, calculated by summing the squares of the differences between each observation and the mean; the regression sum of squares reflects the explanatory power of the fitted values ​​for the variability of the observations, calculated by summing the squares of the differences between each fitted value and the mean; the residual sum of squares reflects the residual variability between the fitted values ​​and the observed values, calculated by summing the squares of the differences between each observed value and its corresponding fitted value, and satisfies the relationship that the total sum of squares equals the sum of the regression sum of squares and the residual sum of squares. Finally, the coefficient of determination R² is calculated as the ratio of the regression sum of squares to the total sum of squares. This value ranges from 0 to 1, with a value closer to 1 indicating better consistency between the model's inverted values ​​and the ground observations, and higher inversion accuracy. Alternatively, the correlation coefficient between the inverted value array and the observed value array can be calculated using the correlation coefficient calculation function in programming software, and then the correlation coefficient can be squared to obtain R². The results from both calculation methods are consistent and can be mutually verified to ensure calculation accuracy.

[0043] During the calculation process, outlier checks must first be performed on the array of inverted and observed values. If the inverted value of a certain weather station differs significantly from the observed value (for example, the observed value is 1 but the inverted value is 0, or the observed value is 0 but the inverted value is 1), it is necessary to check whether the latitude and longitude coordinates of the weather station are accurate and whether there is cloud residue or terrain shadow interference in the corresponding pixels. If it is confirmed that the difference is caused by data anomalies, the abnormal data must be removed and R² must be recalculated to avoid interference from outliers on the accuracy evaluation results. If the calculated R² is not lower than 0.85, the initial snow cover data is considered accurate, and no adjustment of model parameters is needed. If R² is lower than 0.85, the cause of the error needs to be analyzed and optimized accordingly: If the error mainly comes from the snow endmember spectrum (e.g., the inversion values ​​from multiple snow observation locations are generally low), the sample selection area for the snow endmember should be reselected to ensure that there is no vegetation or water mixed in within the selected area, and the endmember extraction and spectral verification process should be repeated; if the error mainly comes from model constraints (e.g., the snow inversion values ​​in water-covered areas are abnormal), the variable weights in the solution function of the programming software can be adjusted, or a reasonable range limit for the proportion of snow endmembers can be added to the objective function to reduce the interference of non-snow endmembers on the snow cover calculation. After optimization, the mixed pixel decomposition should be run again, and R² should be calculated again until R² reaches above 0.85, completing the accuracy verification.

[0044] Topographic correction is used to correct snow cover errors caused by topographic relief. It is based on slope, aspect, and solar parameters in metadata calculated from digital elevation model (DEM) data. Preliminary sub-pixel snow cover data and registered DEM data are loaded into GIS software. The solar incidence angle for each pixel is calculated using 3D analysis tools. Based on the solar incidence angle, the target area is divided into strongly shaded, weakly shaded, and unshaded areas. Different correction coefficients are used for different areas: the snow cover in strongly shaded areas is multiplied by a correction coefficient, the snow cover in weakly shaded areas is multiplied by another correction coefficient, and the snow cover in unshaded areas remains unchanged. The correction coefficients are determined by comparing historical snow cover data of similar features in unshaded areas to ensure consistency in snow cover characteristics between shaded and unshaded areas after correction. After correction, accuracy is verified again. If R² is not lower than 0.85, the results are output; if R² is lower than 0.85, the correction coefficients are further optimized until the accuracy meets the standard.

[0045] Finally, a verified and calibrated high-spatiotemporal-precision sub-pixel snow distribution map is generated and output as a raster product for subsequent application analysis. Visualization of the snow distribution map can be done in GIS software. The calibrated snow cover data is loaded, and a color mapping table is set: snow cover 0-0.2 is set to light white (representing light snow); 0.2-0.5 to white (representing moderate snow); 0.5-0.8 to grayish-white (representing heavy snow); and 0.8-1 to dark grayish-white (representing extremely heavy snow). Basic geographic information such as administrative boundaries of the target area, major rivers and roads, and the locations of ground meteorological stations are added to improve the readability of the distribution map.

[0046] The visualized snow distribution map is exported as a GeoTIFF raster product. During the export process, ensure the file includes geographic information such as coordinate system, cell size, top-left corner coordinates, and rotation angle. A world file can also be generated for easy integration into other GIS software. The exported raster product undergoes quality checks, including spatial positioning accuracy (verified by comparing ground control point coordinates with image coordinates), data integrity (no data holes), and spectral consistency (reasonable variations in snow cover for the same feature across different temporal images), ensuring the product quality meets application requirements.

[0047] The application analysis of this raster product needs to be conducted in conjunction with the actual needs of the target area. Regarding the analysis of snow cover area change trends, sub-pixel snow cover distribution maps for the target area over the past three snow seasons are loaded. Using county-level administrative divisions as statistical units, the monthly snow cover area for each county is calculated (the number of pixels with a snow cover rate greater than 0 multiplied by the area of ​​a single pixel), generating a monthly snow cover area statistical table. Snow cover area change curves can be plotted using office software or programming software to analyze the seasonal variation patterns of snow cover, providing data support for regional snow cover resource assessment. Regarding the maintenance of snow and ice venues, for ski resorts within the target area, snow cover rate data within the venue area is extracted using the cropping tool of GIS software. The average, maximum, and minimum snow cover rates within the venue are calculated. Combined with on-site observed snow thickness data, a "snow cover rate-thickness" correlation is established to analyze the snow thickness distribution and duration within the venue, assess the sustainable utilization cycle of the venue's snow cover resources, and provide a basis for venue maintenance personnel to formulate artificial snowmaking or snow melting and clearing plans. In water resource management, the study analyzes the changing trends of snow cover in major river basins within the target area, combines this with runoff monitoring data from hydrological stations within the basins to analyze the contribution of snowmelt to river replenishment, and predicts water supply and demand in different seasons. This provides decision-making reference for regional water resource management departments in formulating inter-regional water transfer plans. In disaster prevention and mitigation, the study compares snow cover distribution maps before and after heavy snowfall to identify areas with a sudden increase in snow cover. Combined with digital elevation model data, it identifies snow-covered areas with steep slopes (these areas are prone to avalanches) and issues timely avalanche warnings. Simultaneously, it marks townships with high snow cover and high population density to provide key area guidance for snow disaster relief and improve the efficiency of disaster prevention and mitigation response.

[0048] As is known from common technical knowledge, this invention can be implemented through other embodiments that do not depart from its spirit or essential characteristics. Therefore, the disclosed embodiments described above are merely illustrative in all respects and are not the only ones. All modifications within the scope of this invention or its equivalents are included in this invention.

Claims

1. A sub-pixel snow cover mapping method based on medium-resolution remote sensing data, characterized in that, Includes the following steps: Step S1: Acquire medium-resolution remote sensing image data, digital elevation model data, and snow cover observation data from ground meteorological stations for the target area, and perform radiometric calibration and atmospheric correction preprocessing on the medium-resolution remote sensing image data to obtain surface reflectance data. Step S2: Use cloud detection algorithm to identify cloud-covered pixels in surface reflectance data, and use spatiotemporal interpolation and ground object spectral similarity matching to fill in the areas where cloud-covered pixels are located, generating a complete surface reflectance dataset without cloud interference. Step S3: Based on the land cover features of the target area, select snow cover, vegetation and water bodies as spectral endmembers, and extract the spectral curves of pure pixels from the complete surface reflectance dataset without cloud interference to determine the spectral features of each spectral endmember. Step S4: Construct a linear mixing model. The reflectance of each pixel in the complete surface reflectance dataset without cloud interference is represented as a combination of the reflectance of each spectral endmember and its proportion. The least squares method is used to solve the linear mixing model to inversely determine the proportion of snow endmembers in each pixel, thus obtaining preliminary sub-pixel snow cover data. Step S5: Use snow observation data from ground meteorological stations to verify the accuracy of the preliminary sub-pixel snow cover data, and combine it with digital elevation model data to perform terrain correction to correct reflectivity errors caused by terrain undulations. Improve the snow cover inversion accuracy by optimizing the parameters of the linear hybrid model, and finally generate a verified and corrected high spatiotemporal accuracy sub-pixel snow distribution map.

2. The sub-pixel snow cover mapping method based on medium-resolution remote sensing data according to claim 1, characterized in that, In step S1, the medium-resolution remote sensing image data used is surface reflectance product data with daily or 8-day temporal resolution and spatial resolution of 500 meters or 1000 meters. Digital elevation model data is used to provide elevation, slope and aspect information of the target area. Snow cover observation data from ground meteorological stations, including snow depth and snow condition records, are used for model verification and accuracy evaluation.

3. The sub-pixel snow cover mapping method based on medium-resolution remote sensing data according to claim 1, characterized in that, The preprocessing steps for medium-resolution remote sensing image data also include using digital elevation model data to perform topographic correction on the surface reflectance data after radiometric calibration and atmospheric correction, in order to eliminate the influence of topographic shadows on surface reflectance.

4. The sub-pixel snow cover mapping method based on medium-resolution remote sensing data according to claim 1, characterized in that, In the process of filling in the area where cloud-covered pixels are located, the spatiotemporal interpolation method is based on the historical reflectance data of the target pixels during cloudless periods to perform time series interpolation or extrapolation. The land cover spectral similarity matching method is to find cloudless pixels with similar terrain environment and land cover type to cloud-covered pixels in the spatial neighborhood and replace them with their spectral values. By combining these two methods to seamlessly fill in the cloud-covered area, a complete land reflectance dataset without cloud interference is generated.

5. The sub-pixel snow cover mapping method based on medium-resolution remote sensing data according to claim 1, characterized in that, The steps of selecting snow cover, vegetation, and water bodies as spectral endmembers are determined based on the main land cover type of the target area. The step of extracting the spectral curve of pure pixels from the complete surface reflectance dataset without cloud interference is to select a large continuous and spectrally uniform area on the image as the endmember sample selection area, and calculate the average reflectance of all pixels in the selected area as the representative spectral curve of the spectral endmember.

6. The sub-pixel snow cover mapping method based on medium-resolution remote sensing data according to claim 5, characterized in that, After extracting the spectral curves of pure pixels, a spectral angle mapping algorithm is used to calculate the similarity between the spectrum of each pixel in the target area and the spectral curves of each extracted spectral endmember. This is used to verify the spectral purity of the selected endmember samples and to ensure that the endmember spectra used in the linear mixture model can accurately represent the main land cover types in the target area.

7. The sub-pixel snow cover mapping method based on medium-resolution remote sensing data according to claim 1, characterized in that, In the step of solving the linear mixing model using the least squares method to invert the proportion of snow end-members, the linear mixing model assumes that the reflectance of each pixel is a linear combination of the reflectance of its constituent snow end-members, vegetation end-members, and water end-members, weighted by the proportion of their area within the pixel. By solving this system of linear equations using the least squares method, the proportion of each end-member is obtained that minimizes the error between the model's calculated value and the actual reflectance value of the pixel. The proportion of the snow end-member is the preliminary sub-pixel snow cover data.

8. The sub-pixel snow cover mapping method based on medium-resolution remote sensing data according to claim 1, characterized in that, The steps for verifying the accuracy of snow cover observation data from ground meteorological stations are to compare the snow cover observation records at the meteorological station location with the sub-pixel snow cover data at the same location point retrieved by the model, calculate the coefficient of determination R² to quantify the retrieval accuracy, and adjust the parameters of the linear mixture model or the endmember spectral library according to the verification results until the coefficient of determination R² reaches or exceeds the preset threshold.

9. The sub-pixel snow cover mapping method based on medium-resolution remote sensing data according to claim 1, characterized in that, The steps for terrain correction to correct reflectance errors caused by terrain undulations by combining digital elevation model data are as follows: consider the influence of slope and aspect on the intensity of solar radiation received by the surface, establish a correction model using terrain parameters calculated from digital elevation model data, and correct the reflectance values ​​of pixels significantly affected by terrain in the complete surface reflectance dataset without cloud interference, thereby reducing the interference of terrain shadows on the accuracy of mixed pixel decomposition and snow cover inversion.

10. The sub-pixel snow cover mapping method based on medium-resolution remote sensing data according to claim 1, characterized in that, After generating a verified and calibrated high-spatiotemporal-accuracy sub-pixel snow distribution map, the sub-pixel snow distribution map is also output as a raster map product. This is used to analyze the trend of snow area change in the target area at different time scales, or to perform overlay analysis with a specific geographic boundary range to evaluate the distribution characteristics of snow resources within that range and its supporting role in water resource management, ice and snow industry planning, or disaster prevention and mitigation decision-making.

Images