Non-debris coverage glacier remote sensing mapping method based on multi-index probability threshold

Abstract

The invention mainly relates to the technical field of remote sensing, and provides a chipping-coverage-free glacier remote sensing mapping method based on a multi-index probability threshold value in order to improve the image accuracy and precision of a glacier distribution diagram. Establishing a multi-temporal image cube data set, constructing a glacier region recognition model based on the normalized snow index, the normalized water body index and the surface temperature, training the glacier region recognition model based on the multi-temporal image cube data set to generate a glacier region recognition result, and generating a glacier contour according to the glacier region recognition result. The limitation of a traditional single-index fixed threshold value method in a complex environment is overcome through a multi-index collaborative judgment strategy, and the image accuracy and precision of the glacier distribution diagram are improved.

Description

Technical Field

[0001] The present invention mainly relates to the field of remote sensing technology, and in particular to a remote sensing mapping method for glaciers without debris cover based on multi-indicator probability thresholds. Background Art

[0002] Glaciers are an essential component of the global water cycle, and their changes have profound impacts on regional water resource security and ecological stability. Currently, glacier mapping techniques fall into three main categories: multispectral index analysis based on optical satellite imagery, such as the Normalized Difference Snow Index (NDSI) and the Normalized Difference Water Index (NDWI). These methods are widely used to identify glaciers without debris cover due to their ease of use; object-based image segmentation methods; and supervised classification algorithms based on machine learning.

[0003] However, the currently used glacier identification methods generally have the following technical defects: (1) Image quality issues: In the presence of image noise such as thin cloud cover, rough ice surface, and mountain shadows, the single index threshold method is not robust enough and has a high misjudgment rate; (2) Spatiotemporal scale limitations: Existing methods are mainly applied to small areas or short periods of time and lack the ability to systematically monitor glacier changes in large areas over long time series; (3) Incomplete indicator selection: Existing methods often ignore the key role of surface temperature in distinguishing glaciers from glacial lakes; (4) Fixed threshold setting: Traditional methods use fixed threshold judgments and do not consider the probability distribution characteristics of indicators in multi-temporal images. Summary of the Invention

[0004] The technical problem to be solved by the present invention is: a remote sensing mapping method for glaciers without debris cover based on multi-indicator probability thresholds, the purpose of which is to improve the accuracy and precision of glacier mapping.

[0005] The technical solution adopted by the present invention to solve the above technical problems is:

[0006] A method for remote sensing mapping of debris-free glaciers based on multi-index probability thresholds, the method comprising:

[0007] Glacier remote sensing images are collected based on temporal and spatial conditions, preprocessed, and a multi-temporal image cube dataset is established. The glacier remote sensing images in the multi-temporal image cube dataset represent stacked glacier remote sensing image data at different times in the space of the same longitude and latitude.

[0008] A glacier area identification model was constructed based on the normalized snow index, normalized water index, and surface temperature;

[0009] The glacier area recognition model is trained based on the multi-temporal image cube dataset to generate glacier area recognition results;

[0010] Based on the trained glacier area recognition model, glacier area recognition results under target time and space conditions are generated, and glacier contours are generated based on the glacier area recognition results.

[0011] Furthermore, the preprocessing includes: performing atmospheric correction and orthorectification on the collected glacier remote sensing images, and setting the spatial resolution.

[0012] Furthermore, collecting glacier remote sensing images based on temporal and spatial conditions includes:

[0013] Define the spatial and temporal conditions for acquiring remote sensing images;

[0014] Based on the defined spatial and temporal conditions, glacier remote sensing images that meet the temporal and spatial conditions are screened out from the Landsat series image collection.

[0015] Furthermore, building a multi-temporal image cube dataset includes:

[0016] Select pixels with cloud scores less than a set threshold from the filtered glacier remote sensing images that meet the temporal and spatial conditions;

[0017] The selected pixels are grouped according to the set sliding time window, and the spectral data of multiple bands of the corresponding pixels are extracted to construct a multi-temporal image cube dataset.

[0018] Furthermore, constructing a glacier area identification model includes: calculating the normalized snow index, normalized water index and surface temperature for each pixel point in the multi-temporal image cube dataset. If the probability that the calculated values of the normalized snow index, normalized water index and surface temperature of the corresponding pixel point meet the set indicator threshold conditions all meet the set probability threshold conditions, the pixel point will be marked as a glacier area.

[0019] Furthermore, the indicator threshold condition of the normalized snow index is greater than 0.4, the indicator threshold condition of the normalized water index is less than 0.4, and the indicator threshold condition of the surface temperature is lower than -1°C.

[0020] Furthermore, the probability threshold condition for the normalized snow index is greater than 65%, the probability threshold condition for the normalized water index is greater than 70%, and the probability threshold condition for the surface temperature is greater than 55%.

[0021] Furthermore, generating a glacier outline based on the glacier area identification results includes:

[0022] Raster to vector conversion: convert pixel-based glacier area identification results into polygonal vector data;

[0023] Small area filtering: remove glacier polygons with an area smaller than the set value;

[0024] Hole filling: fill the holes inside the glacier polygon whose area is not greater than the set value;

[0025] Morphological processing: smoothing glacier boundaries and reducing jagged edges;

[0026] Boundary optimization: Optimize the glacier boundary using slope and elevation.

[0027] Furthermore, the method also includes: generating multi-period glacier distribution maps and glacier change statistics based on the glacier contours generated corresponding to different time windows.

[0028] Furthermore, the method also includes: evaluating the quality of the generated glacier profile based on correctness, completeness and F1-score indicators.

[0029] The beneficial effects of the present invention are as follows: by establishing a glacier area identification model for collaborative identification of multiple indicators including normalized snow index, normalized water index and surface temperature, glacier areas are identified, overcoming the limitations of the traditional single-indicator fixed threshold method in complex environments, especially showing significant advantages under interference conditions such as image noise, thin cloud occlusion and mountain shadows, and being able to better generate glacier images without debris cover. The accuracy of shadow area identification and seasonal snow differentiation are significantly improved, and the misjudgment rate at the glacier-lake junction is also significantly reduced, which has important reference value for regional water resources management and climate adaptation strategy formulation. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 This is a flow chart of the method for remote sensing mapping of glaciers without debris cover based on multi-index probability thresholds according to the present invention;

[0031] Figure 2 Create a flow chart for a multi-temporal image cube dataset;

[0032] Figure 3 Flowchart of the methodology for building a glacier region identification model;

[0033] Figure 4 A flow chart of the methodology was generated for glacier contours;

[0034] Figure 5 This is the architecture of the multi-indicator probability threshold glacier remote sensing mapping system. DETAILED DESCRIPTION

[0035] like Figure 1 As shown, the method for remote sensing mapping of glaciers without debris cover based on multi-index probability thresholds of the present invention includes the following steps:

[0036] Step 1: Collect glacier remote sensing images based on temporal and spatial conditions, perform preprocessing, and establish a multi-temporal image cube dataset.

[0037] like Figure 2 As shown in the figure, temporal and spatial conditions are first defined. Glacier remote sensing images are then collected based on these conditions. These images are then preprocessed, including atmospheric correction, orthorectification, and spatial resolution. Glacier remote sensing images can be directly acquired using preprocessed Landsat imagery collected using multiple generations of Landsat sensors (TM, ETM+, and OLI). This reduces the amount of data required for the preprocessing process and overcomes the limited temporal coverage of a single sensor. From the collected glacier remote sensing images, images from the glacier melt season (June-August) are selected to maximize the spectral difference between glaciers and non-glaciers. Pixel-level cloud detection is performed: a pixel-level cloud fraction calculation method is used, a cloud fraction threshold is set, and pixels with cloud fractions below the threshold are screened based on a cloud mask function and quality control relationships, achieving fine-grained data quality control. A sliding time window is set, and the selected pixels are grouped according to the set time window. Based on the grouping results, multiple independent image cubes are created to form a multi-temporal image cube dataset. The three dimensions of the multi-temporal image cube dataset are: spatial dimension X - the longitude direction of the image; spatial dimension Y - the latitude direction of the image; and temporal dimension T - the image sequence of different periods. The multi-temporal image cube dataset covers the entire study area spatially and spans multiple years temporally, facilitating the analysis of the changing characteristics of each pixel position over time. Spectral data of multiple bands of pixels in the multi-temporal image cube dataset are extracted for subsequent indicator calculations. The spectral data of multiple bands includes green light band, near infrared band, shortwave infrared band, and thermal infrared band data.

[0038] The mathematical expression of the cloud mask function and quality control relationship is as follows:

[0039] ;

[0040] in, Represents pixels In time The cloud mask value of Indicates the cloud quality score of the pixel; Indicates the set cloud fraction threshold. 1 indicates that the cloud fraction value of the pixel is less than the set cloud fraction threshold, and 0 indicates that the cloud fraction value of the pixel is not less than the set cloud fraction threshold.

[0041] The mathematical definition of the image cube is as follows:

[0042] ;

[0043] Where: represents the image cube; Indicates the pixel position ,wavelength ,time The corresponding spectral values; Indicates the A time window.

[0044] In this example, the time condition is set to 1988-2022, the cloud score threshold is 50, and the sliding time window is selected to be 5 years. The Landsat image acquisition dates and quantities are shown in Table 1. Taking 1988-1992 as an example, which is the first group of time windows, i=1.

[0045] Table 1 Landsat image acquisition date and number

[0046] date Time span sensor Number of Landsat images 1990 1988~1992 TM 2004 1995 1993~1997 TM 2696 2000 1998~2002 TM, ETM+ 2223 2005 2003~2007 TM, ETM+ 2206 2010 2008~2012 TM, ETM+ 2911 2015 2013~2017 OLI 4568 2020 2018~2022 OLI 4530

[0047] Step 2: Construct a glacier area identification model based on the Normalized Difference Snow Index (NDSI), Normalized Difference Water Index (NDWI) and surface temperature.

[0048] NDSI can effectively identify snow and ice, NDWI can effectively distinguish between water and ice, and surface temperature can further confirm the glacier area. A glacier area identification model is constructed based on NDSI, NDWI and surface temperature indicators to achieve multi-dimensional glacier feature discrimination and effectively deal with data missing and outliers. Figure 3 As shown in the figure, the construction of the glacier area identification model includes: setting the NDSI, NDWI and surface temperature index thresholds through quantitative experiments, and calculating the probability that the NDSI, NDWI and surface temperature index meet the corresponding thresholds for each pixel point in the multi-temporal image cube dataset. If the probabilities of the three indicators meeting the thresholds all meet the set probability threshold conditions, the pixel point is marked as a glacier area.

[0049] The NDSI index calculation includes:

[0050] For TM and ETM+ sensors, , for OLI sensors; Where: and For the green light band, and It is the shortwave infrared band.

[0051] The NDWI index calculation includes:

[0052] For TM and ETM+ sensors, ;For OLI sensors, Where: and For the green light band, is the near-infrared band, It is the shortwave infrared band.

[0053] Surface temperature (ST) calculation includes:

[0054] Landsat thermal infrared bands (Band 6 for TM and ETM+, and TIRS for OLI) are used to calculate land surface temperature using the radiation transfer equation and a single-channel algorithm. The radiation transfer equation for land surface temperature calculation is as follows:

[0055] ;

[0056] Where: is the surface temperature (℃); is the radiation brightness received by the sensor; and is the calibration constant. For Landsat-8TIRS, , for Landsat-5 / 7TM / ETM+, ;

[0057] The radiance received by the sensor: ;

[0058] Where: is the radiance of the top of the atmosphere; For the upward radiation of the atmosphere; Radiation downward from the atmosphere; is the atmospheric transmittance; is the surface emissivity.

[0059] For each pixel, the probability that the pixel satisfies the following conditions is calculated, and its normalized mathematical expression is:

[0060] ;

[0061] ;

[0062] ;

[0063] Where: is the indicator function, a is the set NDSI index threshold, b is the set NDWI index threshold, d is the set surface temperature threshold, Indicates the probability that the NDSI index is greater than a, represents the probability that NDWI is less than b, represents the probability that the surface temperature is less than d; Represents a time window The number of valid images in the image.

[0064] Taking the probabilities of the three indicators into consideration, a comprehensive probability model for glacier identification, namely the glacier area identification model, is constructed:

[0065] ;

[0066] Where: This is the glacier region identification result, 1 represents the glacier area, and 0 represents the non-glacier area; 、 and For the setting 、 、 The specific value of the probability threshold can be determined by comparison with reference glacier data (such as RGI6.0) and repeated experiments.

[0067] In this embodiment, a=0.4, b=0.4, d=-1°C are set; 、 、 .

[0068] Step 3: Based on the multi-temporal image cube dataset established in step 1, the glacier area recognition model established in step 2 is trained, and the glacier area recognition results of the target time and space conditions are generated based on the trained glacier area recognition model.

[0069] Step 4: Process the glacier area identification results to generate glacier outlines.

[0070] like Figure 4 As shown in the figure, the processing of glacier area identification results includes: (1) raster to vector conversion: converting pixel-based glacier identification results into polygonal vector data; (2) small area filtering: removing glacier polygons with an area smaller than a set value to reduce misidentification due to spectral confusion or noise; (3) hole filling: filling holes in the glacier polygon with an area of no more than a set square kilometer. These small holes are usually caused by spectral anomalies or gaps in satellite image strips; (4) morphological processing: applying morphological processing such as opening and closing operations to further smooth the glacier boundary and reduce jagged edges; (5) boundary optimization: using auxiliary information such as slope and elevation to further optimize shadow areas such as glacier boundaries.

[0071] The mathematical expression of morphological opening and closing operations is as follows:

[0072] ;

[0073] ;

[0074] ;

[0075] Where, is the original glacier binary image; and They are open and close operations respectively; and As the structural element, a disk-shaped structural element with a radius of 1-3 pixels is usually selected. and are the corrosion and dilation operations respectively.

[0076] Step 5: To evaluate the glacier mapping results, the glacier mapping results were calculated based on the reference debris-free RGI6.0 dataset results. 、 and Indicators, the mathematical definitions of each evaluation indicator are as follows:

[0077] ;

[0078] ;

[0079] ;

[0080] Where: represents the correctly mapped glacier area based on the method of the present invention, i.e., the intersection of the glacier area identified by the method and the reference data; represents the total area of glaciers mapped by this method; is the total area of the reference RGI6.0 glacier without debris cover; represents the glacier area identified by this method; Represents the glacier area of the reference data.

[0081] The debris-free glacier remote sensing mapping method based on multi-index probability thresholds described in the present invention can be completed through the JavaScript or Python API of the GEE platform. The system architecture is as follows: Figure 5 As shown in the figure, (1) data layer: including Landsat series satellite image datasets, reference glacier datasets (such as RGI6.0), digital elevation models (DEMs), etc.; (2) algorithm layer: including image preprocessing module, multi-index calculation module, probability statistics module, glacier identification module and post-processing optimization module; (3) application layer: including glacier mapping results display, glacier change analysis and climate change response analysis, etc.

[0082] The mathematical expression of the system workflow can be summarized as: ;

[0083] Where: Indicates the Glacier mapping results for a time window; Represents the original Landsat image set; represents a time window; Indicates the region of interest; 、 、 They represent the pre-processing function, glacier identification function and post-processing optimization function respectively.

Claims

1. A method for remote sensing mapping of debris-free glaciers based on multi-index probability thresholds, characterized by: The method comprises: Collect glacier remote sensing images based on temporal and spatial conditions, perform preprocessing, and establish a multi-temporal image cube dataset; A glacier area identification model was constructed based on the normalized snow index, normalized water index, and surface temperature; The glacier area recognition model is trained based on the multi-temporal image cube dataset to generate glacier area recognition results; Based on the trained glacier area recognition model, glacier area recognition results under target time and space conditions are generated, and glacier contours are generated based on the glacier area recognition results.

2. The method for remote sensing mapping of glaciers without debris cover based on multi-index probability thresholds according to claim 1, characterized in that: The preprocessing includes: performing atmospheric correction and orthorectification on the collected glacier remote sensing images, and setting the spatial resolution.

3. The method for remote sensing mapping of glaciers without debris cover based on multi-index probability thresholds according to claim 1, characterized in that: Collection of glacier remote sensing images based on temporal and spatial conditions includes: Define the spatial and temporal conditions for acquiring remote sensing images; Based on the defined spatial and temporal conditions, glacier remote sensing images that meet the temporal and spatial conditions are screened out from the Landsat series image collection.

4. The method for remote sensing mapping of glaciers without debris cover based on multi-index probability thresholds according to claim 3, characterized in that: Building a multi-temporal image cube dataset includes: Select pixels with cloud scores less than a set threshold from the filtered glacier remote sensing images that meet the temporal and spatial conditions; The selected pixels are grouped according to the set sliding time window, and the spectral data of multiple bands of the corresponding pixels are extracted to construct a multi-temporal image cube dataset.

5. The method for remote sensing mapping of glaciers without debris cover based on multi-index probability thresholds according to claim 1, characterized in that: Constructing a glacier area identification model includes: calculating the normalized snow index, normalized water index and surface temperature for each pixel point in the multi-temporal image cube dataset. If the probability that the calculated values of the normalized snow index, normalized water index and surface temperature of the corresponding pixel point meet the set indicator threshold meets the set probability threshold condition, the pixel point is marked as a glacier area.

6. The method for remote sensing mapping of glaciers without debris cover based on multi-index probability thresholds according to claim 5, characterized in that: The threshold condition for the normalized snow index is greater than 0.4, the threshold condition for the normalized water index is less than 0.4, and the threshold condition for the surface temperature is lower than -1°C.

7. The method for remote sensing mapping of glaciers without debris cover based on multi-index probability thresholds according to claim 5, characterized in that: The probability threshold condition for the normalized snow index is greater than 65%, the probability threshold condition for the normalized water index is greater than 70%, and the probability threshold condition for the surface temperature is greater than 55%.

8. The method for remote sensing mapping of glaciers without debris cover based on multi-index probability thresholds according to any one of claims 1 to 7, characterized in that: The glacier outline is generated based on the glacier area identification results, including: Raster to vector conversion: convert pixel-based glacier area identification results into polygonal vector data; Small area filtering: remove glacier polygons with an area smaller than the set value; Hole filling: fill the holes inside the glacier polygon whose area is not greater than the set value; Morphological processing: smoothing glacier boundaries and reducing jagged edges; Boundary optimization: Optimize the glacier boundary using slope and elevation.

9. The method for remote sensing mapping of glaciers without debris cover based on multi-index probability thresholds according to claim 8, characterized in that: The method further comprises: generating a multi-period glacier distribution map and glacier change statistics according to the glacier profiles generated corresponding to different time windows.

10. The method for remote sensing mapping of glaciers without debris cover based on multi-index probability thresholds according to claim 8, characterized in that: The method further includes: evaluating the quality of the generated glacier profile based on correctness, completeness and F1-score indicators.

Images