A method for extracting flow lines in a glacier based on terrain intervention

Abstract

The application relates to the technical field of geographic information processing, and discloses a glacier flow line extraction method based on terrain intervention. A digital elevation model and glacier polygon contour data of a target glacier are acquired and subjected to mask processing, a boundary grid node set is screened to determine a potential glacier source point set and a glacier terminal point; contour line distribution and an auxiliary gradient scalar are calculated, a high-value overlapping area is taken to generate a limited terrain intervention mask patch; a basic cost tracking grid is generated based on the shortest distance from a grid point of an elevation grid in the glacier boundary to the glacier boundary and relative elevation; the patch is used to replace a corresponding position of the tracking grid to generate a composite constraint cost grid; the potential glacier source point set is used as a starting point, the glacier terminal point is used as an end point, and an initial flow line network is extracted in the composite constraint cost grid through optimization; after being subjected to smoothing and dimension reduction and redundant branches are removed, glacier flow line data are output. The application can avoid flow line penetration and remove redundant branches, and the extraction precision is improved.

Description

Technical Field

[0001] This invention relates to the field of geographic information processing technology, specifically to a method for extracting streamlines in glaciers based on terrain intervention. Background Technology

[0002] Streamlines in glaciers are fundamental spatial data for characterizing glacier geometry, measuring glacier length, and studying glacier dynamic evolution. Currently, automated extraction of glacier streamlines typically relies on digital elevation models and glacier contour data, utilizing spatial analysis algorithms such as distance transformation, hydrological flow direction analysis, or cost surface tracking to generate streamline networks.

[0003] However, existing methods for streamline extraction in glaciers face significant limitations in complex terrain environments. Actual glaciers typically exhibit topographical features such as undulating ice surfaces, exposed rock walls, or steep ridges. Conventional streamline extraction algorithms often rely solely on boundary Euclidean distance or global elevation descent trends for identification, lacking mechanisms for identifying and addressing local topographical obstacles within the glacier. This approach leads to the algorithm calculating streamlines that easily traverse physically impassable terrain barriers or deviate significantly from the actual flow direction of the glacier during path optimization. Furthermore, limited by the resolution of digital elevation model data and interference from local micro-topography, existing methods tend to generate dense starting points in areas of extreme local elevations when searching for streamline sources. These excessive starting points generate numerous branch streamlines that do not conform to glacier flow characteristics during subsequent path tracing, resulting in a complex final streamline data structure, reduced accuracy of streamline extraction, and increased workload for subsequent manual correction. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a method for extracting streamlines in glaciers based on terrain intervention. This method solves the problems of existing methods lacking identification and intervention mechanisms for local terrain obstacles, being susceptible to micro-topographic interference, and causing extracted streamlines to easily cross terrain obstacles, deviate from the actual flow direction, and generate a large number of redundant branches.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a method for extracting streamlines in glaciers based on topographic intervention, comprising the following steps: Obtain digital elevation model data and polygonal outline data of the target glacier.

[0006] The above digital elevation model data and glacier polygon outline data are processed by masking to obtain the elevation raster inside the glacier boundary, the glacier boundary and the set of boundary grid nodes. The set of boundary grid nodes is filtered by combining local elevation extremes and dynamic spatial distance to determine the set of potential glacier source points and glacier terminus points.

[0007] Based on the digital elevation model data, the contour line distribution gradient scalar and auxiliary gradient scalar are calculated. The high-value overlap area between the contour line distribution gradient scalar and the auxiliary gradient scalar is used to generate a restricted terrain intervention mask patch.

[0008] The basic cost tracking grid is generated based on the shortest distance and relative elevation from the grid points to the glacier boundary using the internal elevation raster of the glacier boundary. The grid values ​​within the restricted terrain intervention mask patch are replaced with preset auxiliary replacement grid parameters whose lower limit is greater than the upper limit of the basic cost tracking grid value. The corresponding positions in the basic cost tracking grid are replaced to generate a composite constrained cost grid.

[0009] Using the set of potential glacier source points as the starting point and the glacier terminus points as the ending point, the initial streamline network is extracted through optimization and iteration in a composite constrained cost grid. The initial streamline network is then smoothed and dimensionality reduced, and redundant branches are removed to output streamline data in the glacier.

[0010] Furthermore, after acquiring the digital elevation model data and glacier polygon outline data of the target glacier, and before masking the digital elevation model data and glacier polygon outline data, the following steps are also included: Extract spatial reference parameters from digital elevation model data and glacier polygon contour data, and unify the coordinate reference system in the spatial reference parameters into the same projected coordinate system; Using the glacier polygon outline data in the projected coordinate system as the clipping boundary, the digital elevation model data is spatially clipped to generate an initial glacier elevation matrix. The two-dimensional projected area of ​​the target glacier is calculated based on the spatial geometric coordinates contained in the glacier polygon outline data. A mapping table containing the mapping relationship between area threshold intervals and target spatial resolution is pre-established; the two-dimensional projected area of ​​the target glacier is matched in the mapping table to obtain the target spatial resolution associated with the matched area threshold interval; When there is a difference between the original spatial resolution and the target spatial resolution of the initial glacier elevation matrix, the bilinear interpolation algorithm is called to resample the initial glacier elevation matrix to generate elevation model data with standard resolution, thereby achieving scale unification of glacier data of different scales.

[0011] Furthermore, in determining the source and end points, starting from the lowest global elevation point in the elevation raster inside the glacier boundary, the outer grid is traversed along the edge of the glacier polygon outline data in a single direction, and the discrete two-dimensional boundary pixels are transformed into a one-dimensional index sequence as a set of boundary grid nodes. Nodes in the boundary grid node set that are locally higher than their neighboring points within the search window and whose absolute elevation is higher than the median elevation value of the elevation raster inside the glacier boundary are extracted. Local elevation extremes are obtained to generate an initial source point set. The initial source point set is then filtered by combining dynamic spatial distance to form a potential glacier source point set. The process involves traversing the elevation grid within the glacier boundary to find the grid cell with the lowest global absolute elevation. This grid cell is then identified as the glacier terminus. The initial set of source points is then filtered using dynamic spatial distance to form a set of potential glacier source points, including: Obtain the glacier area of ​​the target glacier and establish a dynamic distance threshold that includes the mapping relationship between glacier area and spatial distance as the dynamic spatial distance; Calculate the actual spatial Euclidean distance between every pair of points within the initial source point set; When the actual spatial Euclidean distance between adjacent points is less than the dynamic distance threshold, the point with the higher elevation is retained and the point with the lower elevation is removed to form a set of potential glacier source points, thereby reducing the generation of redundant starting points in a small area.

[0012] Furthermore, when extracting the restricted terrain intervention mask patches, the digital elevation model data is called to divide the elevation grid inside the glacier boundary into multiple elevation zones for initial step-coding. Based on the digital elevation model data, the slope value of each grid cell in the encoded grid matrix is ​​calculated to generate the contour line distribution gradient scalar. Extract the median of the code obtained from the initial step coding process and re-perform the elevation zone division operation on both sides to generate a secondary coding grid. Calculate the slope value of each grid cell in the secondary coding grid based on the digital elevation model data to generate an auxiliary gradient scalar. Extract the high-value overlap region corresponding to the contour line distribution gradient scalar and auxiliary gradient scalar in spatial location, construct a binarized mask matrix and assign values ​​to generate the initial steep slope region range.

[0013] The topographical obstacles that hinder glacier flow are identified by the two gradient scalars mentioned above. Then, dilation and erosion operations are performed sequentially on the discretely distributed high-value overlapping areas within the initial steep slope region to generate independent connected patches. Count the number of grid cells contained in each independent connected patch; Independent connected patches are removed when the number of grid cells is lower than the preset area lower limit threshold. Independent connected patches with a number of grid cells not lower than the preset area lower limit threshold are retained as restricted terrain intervention mask patches.

[0014] Furthermore, when constructing the basic cost tracking grid, the resistance to glacier flow is reflected by the superposition of two cost parameters; The difference between the maximum Euclidean distance between all grid points in the elevation raster inside the glacier boundary and the shortest Euclidean distance from the grid point to the nearest glacier boundary is calculated. The ratio of this difference to the maximum Euclidean distance between all grid points in the elevation raster inside the glacier boundary is then calculated to obtain the first term, which reflects the trend of grid points converging towards the center of the glacier. The ratio of the difference between the absolute elevation value of the grid point and the global minimum elevation of the elevation grid inside the glacier boundary, to the difference between the global maximum elevation and the global minimum elevation of the elevation grid inside the glacier boundary, is used to obtain the second term of the mapped relative elevation, which reflects the trend of elevation decline. The cost value of each grid point is generated by superimposing the first and second items. The cost values ​​of each grid point are then aggregated to generate a basic cost tracking grid. The corresponding positions of the basic cost tracking grid are replaced with restricted terrain intervention mask patches whose lower limit is greater than the upper limit of the basic cost tracking grid, thereby increasing the passage cost in steep slope areas.

[0015] Furthermore, in the streamline optimization extraction stage, the effective cells in the composite constraint cost grid are abstracted into independent nodes in the grid, and a bidirectional topological connection between adjacent nodes in the independent nodes is established based on the eight-neighborhood connectivity rule to map into a directed weighted graph. The passage cost weight between adjacent nodes is set as the product of the central tendency measure of the cost value mapped by the adjacent nodes in the composite constrained cost grid and the spatial geometric distance coefficient between them; Each source point in the set of potential glacier source points is designated as an independent topological search starting point, and the glacier terminus is set as a global target point; The hydrological flow direction minimum cost algorithm is executed to obtain multiple initial candidate streamlines connected to the global target point, and these multiple initial candidate streamlines are merged into an initial streamline network.

[0016] Finally, point-to-point collision detection of spatial coordinates is performed on the node sequence of the initial candidate streamlines in the initial streamline network. When two initial candidate streamlines originating from different source points intersect in space and their subsequent node sequences completely overlap, the overlapping streamlines are reduced to a single streamline in the main stream of the glacier. Calculate the actual spatial length of non-overlapping tributaries, obtain the maximum length of the midstream line of the target glacier, multiply it by a preset pruning ratio coefficient as the pruning threshold, and when the total length of non-overlapping tributaries is lower than the pruning threshold, the non-overlapping tributaries are stripped from the network to remove redundant branches. The curve fitting algorithm is invoked to perform spatial geometric smoothing on the discrete node sequence, outputting streamline data in the glacier.

[0017] This invention provides a method for extracting streamlines in glaciers based on topographic intervention. It has the following beneficial effects: 1. This invention calculates the high-value overlap area of ​​the contour line distribution gradient scalar and the auxiliary gradient scalar to generate a restricted terrain intervention mask patch. The grid values ​​in this patch are replaced with preset auxiliary replacement grid parameters whose lower limit is greater than the upper limit of the base cost tracking grid value. This can effectively identify physical obstacles in the terrain that hinder glacier flow and guide the streamline optimization algorithm to avoid these obstacle areas by increasing the passage cost of the corresponding area. This avoids the extracted glacier streamlines crossing terrain obstacles and improves the accuracy of streamline extraction under complex terrain.

[0018] 2. This invention generates a basic cost tracking grid based on the shortest distance and relative elevation from the grid points of the elevation raster inside the glacier boundary to the glacier boundary; the shortest distance parameter causes streamlines to converge towards the center of the glacier interior, and the relative elevation parameter causes streamlines to extend in accordance with the trend of decreasing altitude; at the algorithm level, it reflects the physical constraints of glacier flow, so that the extracted streamlines in the glacier can be stably kept in the central region of the glacier body and avoid the streamlines from deviating.

[0019] 3. This invention combines dynamic spatial distance to filter the initial source point set, and after extracting the initial streamline network, performs point-to-point collision detection of the node sequence of the initial candidate streamlines and removes non-overlapping branches based on the pruning threshold; it can reduce redundant source points caused by minor local topographic undulations in the initial stage, and remove redundant branches that do not meet the length conditions after streamlines converge, thereby simplifying the final output streamline data in the glacier and removing redundant streamline branches. Attached Figure Description

[0020] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is a schematic diagram of the system architecture of the present invention. Detailed Implementation

[0021] The technical solutions in 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, and not all embodiments. 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.

[0022] Please see the appendix Figure 1 This invention provides a method for extracting streamlines in glaciers based on topographic intervention, comprising the following steps: S100, acquire digital elevation model data and glacier polygon outline data of the target glacier; S200 uses a source-sink point identification module to perform masking processing on digital elevation model data and glacier polygon contour data to obtain the elevation raster inside the glacier boundary and the set of boundary grid nodes. Based on the local elevation extreme value and dynamic spatial distance judgment conditions, the set of boundary grid nodes is filtered to determine the set of potential glacier source points and glacier terminus points of the target glacier. S300 uses a steep slope terrain extraction module to discretize and calibrate digital elevation model data according to a preset elevation interval, calculates the contour distribution gradient scalar of spatial grid cells, calculates the elevation median benchmark offset of calibration parameters to generate an auxiliary gradient scalar, takes the high value overlap area between the contour distribution gradient scalar and the auxiliary gradient scalar, and generates a restricted terrain intervention mask patch through morphological operations. S400 uses the cost grid construction module to calculate the shortest Euclidean distance and relative elevation mapping value from different grid points to the boundary of the elevation raster inside the glacier boundary, generating a basic cost tracking grid. Based on the distribution range of the restricted terrain intervention mask patch, an auxiliary replacement grid is obtained. The grid cells in the restricted terrain intervention mask patch area are replaced with auxiliary replacement grid parameters whose lower limit is greater than the upper limit of the basic cost tracking grid, generating a composite constrained cost grid. S500 uses the centerline extraction module to start with the set of potential glacier source points, performs hydrological flow direction minimum cost optimization iteration in the composite constrained cost grid, obtains the initial vector streamline network connected to the glacier terminus points, performs dynamic window smoothing operation and node dimension reduction intersection processing on the initial vector streamline network, establishes a pruning threshold based on the length of the main streamline to remove redundant branch streamlines, and outputs streamline data in the glacier.

[0023] Please see Figure 2 The following section elaborates on step S100, which involves acquiring the digital elevation model data and polygon outline data of the target glacier. Specifically, it can be divided into the following sub-steps: S110, for the target extraction area, acquire the original digital elevation model data and the glacier polygon outline data reflecting the spatial distribution of the glacier. In this embodiment, the original digital elevation model data is represented as a regular raster array containing absolute elevation values, and the glacier polygon outline data is represented as a set of closed vector polygons containing the spatial coordinates of glacier boundary nodes. The acquisition methods for the original digital elevation model data and the collection methods for the glacier boundary vector data can be implemented by downloading from the internet; the acquisition process is well-known in the field and will not be described in detail here.

[0024] S120, based on the aforementioned basic data, extract the spatial reference parameters of both and determine whether the coordinate reference systems are consistent. If the reference systems are inconsistent, a coordinate transformation algorithm is used to unify the original digital elevation model data and the glacier polygon outline data to the same equal-area or equal-angle projection coordinate system, eliminating distance calculation deviations introduced by projection distortion. Using the glacier polygon outline data after unifying the coordinate system as the clipping boundary, the original digital elevation model data is spatially clipped. During this mask extraction process, the algorithm assigns null labels to the useless background raster data outside the glacier polygon outline to achieve physical isolation at the computational logic level, thereby accurately extracting the effective elevation raster inside the glacier polygon outline and generating the initial glacier elevation matrix.

[0025] S130. Due to the varying resolution requirements for elevation data in streamline extraction from glaciers of different scales, the initial glacier elevation matrix needs adaptive scaling based on the glacier area. The two-dimensional projected area of ​​the target glacier is calculated based on the spatial geometric coordinates contained in the glacier polygon outline data. The specific solution for this two-dimensional projected area depends on the sum of the cross products of the polygon vertex coordinates, and its calculation logic follows the conventional polygon area calculation rules in surveying. Based on the calculated two-dimensional projected area, a grid resolution adjustment operation is performed on the initial glacier elevation matrix. A mapping table containing multiple area threshold intervals is pre-established, with each area threshold interval associated with a specific target spatial resolution parameter; as a specific implementation, this mapping table can be set as follows: when the two-dimensional projected area of ​​the target glacier is less than 1 km²... 2 At that time, the target spatial resolution was set to 10m to preserve high-frequency terrain details; when the area is between 1km², the spatial resolution is set to 10m. 2 Up to 10km 2 When the area is between 10 km², the target spatial resolution is set to 30 m; when the area is greater than 10 km², the resolution is set to 30 m. 2 At this time, the target spatial resolution is set to 50m or lower. This hierarchical strategy can effectively avoid the risk of memory overflow caused by the massive redundant grid data in large glaciers. The physical basis for establishing this mapping table is that, in order to balance the accuracy and computational efficiency of the streamline, a higher resolution is suitable for small-area glaciers to improve accuracy, while a lower resolution is suitable for large-scale glaciers to improve computational efficiency.

[0026] The 2D projected area of ​​the target glacier is matched to the corresponding area threshold range to obtain the target spatial resolution. When there is a difference between the original spatial resolution of the initial glacier elevation matrix and the target spatial resolution, a bilinear interpolation algorithm is used to resample the initial glacier elevation matrix. The bilinear interpolation algorithm assigns smoothly transitioned elevation values ​​to the resampled grid nodes by calculating the distance-weighted average of the elevation values ​​of the four adjacent original grid cells around the center point of the newly distributed grid. In this distance-weighted calculation stage, to avoid the mathematical singularity collapse problem where the distance denominator approaches zero due to the complete overlap of the spatial positions of the new grid coordinates generated by resampling with the original grid nodes, a very small perturbation constant is added to the distance denominator of the interpolation calculation to ensure the completeness of the floating-point division operation. The resampling operation generates standard-resolution elevation model data, which serves as the basic input data for subsequent terrain intervention calculations and grid cost tracking.

[0027] Please see Figure 2 The working logic of step S200, which uses the source-sink point identification module to determine the set of potential glacier source points and glacier terminus points of the target glacier, is described in detail and can be divided into the following sub-steps: S210, based on the general technical principle that extracting streamline networks requires clearly defining the starting point of topology tracing, the entity range of the target glacier needs to be pre-extracted from the continuous elevation model. In this embodiment, based on the previously acquired elevation model data and glacier polygon outline data, masking is performed to separate the effective shape. During the masking process, to avoid numerical confusion between the absolute elevation of low-altitude glacier areas in the real terrain and the background data, the algorithm assigns null labels or preset minimum constants to the grid cells outside the glacier polygon outline, thereby achieving spatial topological isolation from the external data-free areas at the mathematical calculation level, thus obtaining an internal elevation raster of the glacier boundary containing only absolute elevation information. Based on this elevation raster, all continuous glacier boundary lines, including the outlines of internal islands and glacier walls, are extracted. To construct a one-dimensional topology search sequence on the boundary, the algorithm starts from the global lowest elevation and sequentially sorts the grid points on each boundary line to generate a set of boundary grid nodes for subsequent extreme value searches.

[0028] The position sorting operation can use the Moore neighborhood tracing algorithm or the eight-connected boundary tracing method. Taking the lowest global elevation point as the initial anchor point, the outer grid is traversed sequentially in a single clockwise or counterclockwise direction along the edge of the glacier mask. This transforms the discrete two-dimensional boundary pixels into a one-dimensional index sequence with a strict topological order, ensuring the continuity of subsequent local extremum comparisons.

[0029] S220, for the extracted boundary grid node set, establish local elevation comparison rules and global altitude distribution constraints, and preliminarily screen potential glacier source points using a multi-dimensional cross-validation method. Establishing potential glacier source points requires meeting two core criteria. The first criterion requires that the elevation of the boundary grid point is locally higher than that of adjacent points distributed along the glacier boundary line. In this embodiment, the search window for local comparison is set to 15 adjacent points on each side of the target point. The determination of this threshold number of neighboring points is based on balancing the sensitivity and anti-interference ability of source identification. If the search length set by the threshold is too long, some potential glacier source points of smaller branches may not be identified due to local topographic undulations; conversely, if the search length set by the threshold is too short, it will lead to too many small extreme points representing topographic noise along the boundary. As another screening criterion, the second criterion requires that potential glacier source points are usually located in higher altitude areas. The physical purpose of introducing this criterion is to avoid incorrectly guiding the main centerline to low-lying tributary areas. Considering that the glacier equilibrium line height is often missing due to lack of long-term monitoring, this scheme uses the median elevation value of the elevation grid within the entire glacier boundary as a benchmark to replace the equilibrium line height. Based on the above two criteria, nodes that simultaneously satisfy the conditions of having the largest local neighborhood elevation and an absolute elevation strictly higher than the median elevation value are extracted from the boundary grid node set, forming a preliminary set of potential glacier source points.

[0030] S230. Given the offset of glacier boundaries in high-altitude areas due to misjudgments of snow and ice during glacier cataloging, relying solely on local extrema and elevation limits easily generates densely distributed redundant false sources in nearby areas. To eliminate unreasonable potential glacier sources caused by this local topographic disturbance, a dynamic spatial distance function based on glacier area is established to perform proximity rejection filtering on the preliminary set of potential glacier source locations. A dynamic distance threshold is set, requiring the spatial distance between adjacent potential glacier source locations to be greater than the threshold value. This dynamic distance threshold With glacier area (Unit: km) 2 The relationship between the two functions is functional, and its calculation formula is as follows: ; in, This represents the dynamic distance threshold used to filter adjacent potential glacier source locations; Indicates the area of ​​the glacier; This represents the area weighting coefficient, and its value is set to 2 × 10. -6 m -1 ; This represents the basic distance constant, which is set to 500m to ensure that there is a basic interval between adjacent sources of large glaciers. This represents the preset maximum distance limit. The above area weighting coefficient... Distance constant from the base The values ​​are based on empirical statistical values ​​derived from multi-period historical glacier inventory data of the target area. Among them, the following is set: The 500m setting is to ensure a basic physical distance between adjacent sources, even in small glaciers, preventing the dense accumulation of false sources caused by localized micro-topographical undulations of moraine; [This setting is incomplete and requires further context.] 2×10 -6 m -1 This follows the natural spatial scale law of glacier accumulation zone expansion with overall area, that is, the larger the area of ​​a valley glacier, the wider its source catchment basin, and the required source repulsion interval threshold should increase linearly accordingly.

[0031] Based on the aforementioned distance function, the initial set of potential glacier source points is traversed, and the actual spatial Euclidean distance between each pair of points is calculated. When the actual Euclidean distance between adjacent points is less than the currently calculated dynamic distance threshold... When a node filtering operation is triggered, for source point pairs determined to be too close, the absolute elevation values ​​of their respective grids are compared. The pair with the higher elevation is retained as the true potential glacier source, while the pair with the lower elevation is discarded. As a mandatory safety net mechanism to maintain the completeness of the algorithm, if the set of potential glacier source points is completely cleared after the above distance filtering process, the highest global elevation point is directly regarded as the potential glacier source point.

[0032] S240, after extracting the source based on the above multi-dimensional filtering, the final water catchment target of the topological network needs to be established. The unique glacier terminus of the target glacier must be identified. For determining the glacier terminus, a global search is performed based on the fundamental geophysical law that glacial material flows continuously to the lowest point driven by gravity. The elevation grid within the glacier boundary is traversed to find the grid cell with the lowest global absolute elevation, and this lowest elevation point is strictly defined as the glacier terminus. This point will be designated as the unique mandatory termination sink in the subsequent streamline optimization iterative network.

[0033] Please see Figure 2 The working logic of step S300, which uses the steep slope terrain extraction module to discretize and calibrate the digital elevation model data according to the preset elevation interval and extract the limiting terrain intervention mask patches, is described in detail. Specifically, it can be divided into the following sub-steps: S310, based on the geological principle that steep rock walls and ridgelines are common in high-altitude glacier regions, and that such terrain features significantly influence the flow direction of the centerline, it is necessary to delineate the boundary region by calculating the elevation gradient. In this embodiment, the steep slope terrain extraction module discretizes and calibrates the elevation raster inside the glacier boundary generated in the previous step according to a preset elevation interval. The algorithm divides the continuous elevation model into multiple elevation zones at 50-meter intervals. The algorithm encodes the elevation zones from low to high as 1, 2, 3, and n, completing the initial step-wise encoding process of the elevation data. Based on this encoded grid sequence, the slope value of each grid cell is calculated using a conventional spatial partial derivative algorithm, and its contour distribution gradient scalar is obtained. The calculation formula is as follows: ; in, Represents the gradient scalar of contour line distribution; express The change in elevation per unit distance along an axis (one of the horizontal directions); express The change in elevation per unit distance along the axial direction (vertical direction). This occurs during the difference operation performed on the raster matrix. and The denominator of the partial derivative of the direction is the fixed spatial resolution size of the grid cell. Its physical properties determine that the denominator is always a positive real constant, thus ensuring that the algorithm will not trigger a division-by-zero error during global traversal.

[0034] S320, considering that relying solely on fixed-starting absolute elevation zone division inevitably introduces truncation errors at contour line boundaries, some originally gently sloping grids are assigned different codes after discretization simply because their elevations cross the boundary, leading to false responses that are systematically misjudged as high gradients. To eliminate this numerical interference, the algorithm introduces a baseline offset for secondary verification of the calibration parameters. Specifically, after extracting the median of the original codes, the elevation zone division operation is re-performed at 50-meter intervals to both sides. Based on the re-coded grids after this offset, the above calculation formula is reapplied to recalculate the contour gradients for each grid, generating auxiliary gradient data for spatial cross-comparison.

[0035] S330: Extract the common high-value grid region corresponding to the spatial locations of the contour line distribution gradient scalar and the auxiliary gradient scalar. Traverse the contour line gradient data of the two sets of grids. Establish a decision rule based on dual logic: if a grid cell has a non-zero value in both datasets, it indicates that the terrain steepness of the region has passed the misalignment test, and a value of 1 is assigned to it in the binarized mask matrix; if a value of 0 exists in either dataset, a value of 0 is assigned. This yields the initial steep slope region range, presenting a discrete scattering state.

[0036] S340 utilizes mathematical morphology operations to reconstruct the connectivity of discrete initial steep slope land regions. Due to limitations of the raster partial difference algorithm, the initially extracted scattered points failed to form continuous physical barriers. In this embodiment, the impdilate and imerode functions in the mathematical morphology software are called to process the discretely distributed scattered points. The impdilate operation expands the boundary of the target points to bridge discrete points with similar physical distances; the subsequent imerode operation shrinks the expanded outer boundary and removes isolated noise points that could not be connected. After the above processing, the scattered points have formed continuous restricted terrain intervention mask patches containing ridgelines, exposed rock walls, and steep glacier areas; to further improve the purity of the mask and reduce the complexity of subsequent cost grid calculations, a patch area threshold filtering mechanism is introduced after morphological connectivity reconstruction. The number of grid cells contained in each independent connected patch is counted. If the total number of grid cells in a patch is lower than the preset area threshold (e.g., less than 10 pixel cells), it is judged as invalid terrain random noise and removed. The large patches with substantial physical blocking significance are retained as the final restricted terrain intervention mask patches.

[0037] Please see Figure 2 The working logic of step S400, which uses the cost grid construction module to calculate the basic cost tracking grid and apply steep slope penalty constraints to generate a composite constraint cost grid, is described in detail. Specifically, it can be divided into the following sub-steps: S410, based on the general technical principle that hydrological topology optimization relies on spatial walk costs, requires pre-converting the geometric shape and elevation distribution into computer-recognizable constraint weight parameters. In this embodiment, based on the continuous restricted terrain intervention mask patches extracted in the previous steps and the initial contour of the glacier polygon, boundary grid points for cost grid calculation are merged. For effective grid cells inside the glacier, the shortest Euclidean distance between each current grid point and the nearest glacier boundary grid point is calculated. This distance parameter geometrically reflects the degree to which the grid cell is far from the glacier boundary. Simultaneously, the absolute elevation values ​​of the corresponding grid cells are extracted to provide input for the subsequent establishment of the cost function.

[0038] According to glacier dynamics, the tracking path of the main streamline in S420 must simultaneously satisfy the dual geometric trends of convergence towards the glacier center and continuous discharge towards the lower elevation terminus. To achieve this goal, a cost grid construction module is used to calculate the flow path at any grid point within the glacier. Cost value The calculation formula is as follows: ; in, Indicates the cost value; Represents grid points Euclidean distance to the nearest glacier boundary; This represents the maximum Euclidean distance between all grid points inside the glacier; Represents grid points The absolute elevation value; Indicates the maximum elevation of the glacier; This represents the minimum elevation of the glacier. In the actual algorithm execution logic, since any physical glacier has an objective two-dimensional coverage area and a three-dimensional spatial drop, its maximum distance and overall elevation difference must be positive real numbers greater than zero. This physical property fundamentally ensures that the denominator of the division is never zero.

[0039] In the above formula, the first term forces water to flow towards the center of the glacier, with the cost mapping value approaching 0 the farther away from the boundary, and the cost value of the grid closer to the boundary approaching 1. The second term drives water to flow towards lower elevation areas, with the cost approaching 0 in lower elevation areas and approaching 1 in higher elevation areas. By superimposing these two terms, the basic cost tracking grid for the glacier region is generated.

[0040] S430. Considering that the previously generated mesh cannot identify steep terrain, it needs to be re-identified to ensure the smooth progress of the centerline tracking process. The algorithm temporarily removes the force of continuous terrain intervention mask patches. Without considering the obstruction of steep slope areas, the above formula is reused to calculate the cost of each mesh globally, generating a set of auxiliary replacement meshes.

[0041] S440, based on the absolute mapping relationship of spatial coordinates, uses the cost grid in the auxiliary replacement grid as a substitute for missing areas in non-steep slope regions. However, the local gradient in steep slope regions is often greater than the global gradient. If only direct replacement is performed, the centerline will be dominated by the elevation difference and collide with the boundary. As a preferred terrain avoidance protection mechanism, the algorithm imposes mandatory numerical constraints on the grid located in steep slope regions. The minimum value of the replacement cost grid must be greater than the maximum value of the cost tracking grid in non-steep regions. Only in this way can centerline tracking extend from the potential energy head to the glacier terminus. Using this composite constraint, topological connectivity is preserved, and the elevation of steep slope regions is normalized so that centerline tracking follows the natural elevation variation trend, forcing the centerline in steep slope regions to mostly extend towards the glacier terminus.

[0042] Please see Figure 2 The working logic of step S500, which uses the centerline extraction module to perform path optimization and extract the target glacier trunk network based on the composite constraint cost grid, is described in detail and can be divided into the following sub-steps: S510, based on the general technical principle that streamline tracing in spatial hydrological analysis relies on discrete cost networks, requires pre-converting continuous two-dimensional raster data into a computer graph theory model with topological connectivity attributes. In this embodiment, the centerline extraction module is used to map the previously constructed composite constraint cost grid into a directed weighted graph. The algorithm abstracts each effective cell in the cost grid as an independent node in the graph theory network and establishes bidirectional topological connections between adjacent nodes based on the eight-neighborhood connectivity rule. To expand the applicability of the technical solution and prevent non-core algorithms from being easily circumvented, this embodiment uses universal graph theory logic to define passage resistance. The passage cost weight between adjacent nodes is defined as the product of the central tendency measure of the mapped cost value of the two nodes in the cost grid and the spatial geometric distance coefficient between them. When the two nodes are in orthogonal directions, the value of this distance coefficient is equal to the spatial resolution size of the grid; when the two nodes are in diagonal directions, the value of the distance coefficient is equal to the spatial resolution size multiplied by... This coefficient effectively eliminates the diagonal movement deviation introduced by the fixed topology of the square grid.

[0043] S520, relying on the directed weighted graph network established above, performs graph theory-based network path optimization. It calls the set of potential glacier source points and the unique glacier terminus point generated by the preceding module. In the algorithm's iterative logic, each potential glacier source point is designated as an independent topological search starting point, and the glacier terminus is set as the global target point. The algorithm calls the lowest-cost algorithm with a global optimal search guarantee to traverse network nodes; in this embodiment, the shortest path algorithm can specifically adopt Dijkstra's algorithm or the A* heuristic search algorithm. When using the A* algorithm, the heuristic function can be set as the weighted sum of the straight-line Euclidean distance from the current node to the glacier terminus and the global minimum elevation difference. This significantly reduces the iteration breadth in useless directions while ensuring the global minimum cost path is found, improving the traversal efficiency of the graph network. During the calculation process, the algorithm continuously evaluates and updates the cumulative cost from the initial source to the current affected node and records the trajectory of the predecessor node corresponding to the minimum cost. Considering the significant numerical penalty imposed in the preceding steps, some extremely isolated pseudo-sources may not be able to connect to the target point. To ensure completeness, a connectivity anomaly breaker mechanism is introduced. If the search traversal process of a source exhausts the connected subgraph without reaching the glacier terminus, the potential source is directly determined as an invalid pseudo-starting point and removed from the network. This yields multiple initial candidate streamlines that extend to the same terminus.

[0044] S530, considering the physical macroscopic morphology of multi-level tributary confluence in natural glacial basins, multiple initial candidate streamlines inevitably share the same discharge channel after flowing into the main glacier. Therefore, redundancy elimination and trunk fusion of their topology are necessary. All generated initial candidate streamline node sequences are traversed, and point-to-point collision detection of spatial coordinates is performed. When two candidate streamlines originating from different sources are determined to intersect in a certain space and their subsequent node sequences completely overlap, topology fusion logic is triggered, reducing the overlapping streamlines to a unique glacier trunk centerline. A micro-branch pruning mechanism is introduced. The actual spatial length of each independent non-overlapping tributary is calculated. If the total length of an independent branch is strictly lower than a threshold established based on the maximum length of the target glacier, it is determined to be a false tributary and removed from the network; this threshold can be calculated by obtaining the longest median streamline length of the target glacier as a benchmark and multiplying it by a preset pruning ratio coefficient (e.g., 5% to 10%). Using a relative dynamic ratio rather than an absolute length value for pruning ensures that the same algorithm is compatible with the adaptive extraction requirements of glaciers ranging from micro-hanging glaciers to giant valley glaciers, preventing the erroneous deletion of effective short tributaries in large glaciers.

[0045] S540, continuous smoothing processing is performed on the discrete streamline network after the above-mentioned fusion and pruning. Since the initially extracted paths are limited by the fixed pixel resolution of the raster data, their spatial morphology typically presents a jagged, broken-line appearance. To restore the continuous physical properties of fluid motion, a curve fitting algorithm is used to perform spatial geometric smoothing on the discrete node sequence. This operation reduces the step-turning interference caused by gridding without changing the macroscopic orientation by defining local control points and performing high-order continuous differentiable interpolation between the broken-line nodes. After processing, the algorithm finally outputs a set of glacier trunk network topology vector data that conforms to the geoscientific continuous smoothing properties and has complete projection information.

Claims

1. A method for extracting streamlines in glaciers based on topographic intervention, characterized in that, Includes the following steps: Acquire digital elevation model data and polygon outline data of the target glacier; The digital elevation model data and the glacier polygon outline data are masked to obtain the elevation raster inside the glacier boundary, the glacier boundary and the set of boundary grid nodes. The set of boundary grid nodes is filtered by combining local elevation extremes and dynamic spatial distance to determine the set of potential glacier source points and glacier terminus points. The digital elevation model data is called to divide the elevation grid inside the glacier boundary into multiple elevation zones for initial step-coding. Based on the digital elevation model data, the slope value of each grid cell in the encoded grid matrix is ​​calculated to generate a contour line distribution gradient scalar. Extract the median of the code obtained from the initial step-coding process and re-perform the elevation zone division operation on both sides to generate a secondary coding grid. Calculate the slope value of each grid cell in the secondary coding grid based on the digital elevation model data to generate an auxiliary gradient scalar. Extract the high-value overlap region between the contour line distribution gradient scalar and the auxiliary gradient scalar in spatial location, construct a binary mask matrix and assign values ​​to generate the initial steep slope region range; Dilation and corrosion operations are sequentially performed on the discretely distributed high-value overlapping areas within the initial steep slope region to generate independent connected patches; Count the number of grid cells contained in each independent connected patch; When the number of grid cells in an independent connected patch is lower than a preset area lower limit threshold, it is removed, and independent connected patches with a number of grid cells not lower than the preset area lower limit threshold are retained as restricted terrain intervention mask patches. The first term is obtained by subtracting the shortest Euclidean distance from the grid point to the nearest glacier boundary from the maximum Euclidean distance between all grid points in the elevation grid inside the glacier boundary. The ratio of this difference to the maximum Euclidean distance between all grid points in the elevation grid inside the glacier boundary is then calculated. The second term of the mapped relative elevation is obtained by calculating the ratio of the difference between the absolute elevation value of the grid point and the global minimum elevation of the elevation grid inside the glacier boundary, and the difference between the global maximum elevation and the global minimum elevation of the elevation grid inside the glacier boundary. The cost value of each grid point is generated by superimposing the first and second items; the cost values ​​of each grid point are then aggregated to generate a basic cost tracking grid. The grid values ​​within the restricted terrain intervention mask patch are replaced with preset auxiliary alternative grid parameters whose lower limit is greater than the upper limit of the base cost tracking grid value, and the corresponding positions in the base cost tracking grid are replaced to generate a composite constrained cost grid. Using the set of potential glacier source points as the starting point and the glacier terminus points as the ending point, an initial streamline network is extracted through optimization and iteration in the composite constrained cost grid. The initial streamline network is then smoothed and dimensionality reduced, and redundant branches are removed to output streamline data in the glacier.

2. The method for extracting streamlines in glaciers based on topographic intervention according to claim 1, characterized in that, After acquiring the digital elevation model data and glacier polygon outline data of the target glacier, and before the masking process of the digital elevation model data and glacier polygon outline data, the method further includes: Extract spatial reference parameters from the digital elevation model data and the glacier polygon outline data, and unify the coordinate reference system in the spatial reference parameters into the same projected coordinate system; Using the glacier polygon outline data within the projected coordinate system as the clipping boundary, the digital elevation model data is spatially clipped to generate an initial glacier elevation matrix.

3. The method for extracting streamlines in glaciers based on topographic intervention according to claim 2, characterized in that, After generating the initial glacier elevation matrix, the process includes: The two-dimensional projected area of ​​the target glacier is calculated based on the spatial geometric coordinates contained in the glacier polygon outline data. A mapping table containing the mapping relationship between area threshold intervals and target spatial resolution is pre-established; The two-dimensional projected area of ​​the target glacier is matched in the mapping table to obtain the target spatial resolution associated with the matched area threshold interval; When there is a difference between the original spatial resolution of the initial glacier elevation matrix and the target spatial resolution, the bilinear interpolation algorithm is invoked to resample the initial glacier elevation matrix and generate elevation model data with standard resolution.

4. The method for extracting streamlines in glaciers based on topographic intervention according to claim 1, characterized in that, The process of combining local elevation extremes and dynamic spatial distances to filter the boundary grid node set and determine the potential glacier source locations and glacier terminus locations includes: Starting from the global lowest elevation point in the elevation raster inside the glacier boundary, traverse the outer grid along the edge of the glacier polygon outline data in a single direction, and convert the discrete two-dimensional boundary pixels into a one-dimensional index sequence as the boundary grid node set. Extract nodes from the boundary grid node set that are locally higher than their neighboring points within the search window and whose absolute elevation is higher than the median elevation value of the elevation raster inside the glacier boundary. Obtain the local elevation extreme values ​​to generate a preliminary set of potential glacier source points. Combine the dynamic spatial distance to filter the preliminary set of potential glacier source points to form the potential glacier source point set. The grid cell with the lowest global absolute elevation is found by traversing the elevation grid inside the glacier boundary, and the grid cell with the lowest global absolute elevation is determined as the glacier terminus.

5. The method for extracting streamlines in glaciers based on topographic intervention according to claim 1, characterized in that, Using the set of potential glacier source points as the starting point and the glacier terminus points as the ending point, an initial streamline network is extracted through iterative optimization within the composite constrained cost grid, including: The effective cells in the composite constrained cost grid are abstracted as independent nodes in the grid, and bidirectional topological connections between adjacent nodes in the independent nodes are mapped to a directed weighted graph based on the eight-neighbor connectivity rule. The passage cost weight between adjacent nodes is set as the product of the central tendency measure of the cost value mapped by the adjacent nodes in the composite constrained cost grid and the spatial geometric distance coefficient between them; Each source point in the set of potential glacier source points is designated as an independent topological search starting point, and the glacier terminus is set as a global target point; The hydrological flow direction minimum cost algorithm is executed to obtain multiple initial candidate streamlines connected to the global target point, and the multiple initial candidate streamlines are merged into the initial streamline network.

6. The method for extracting streamlines in glaciers based on topographic intervention according to claim 5, characterized in that, The initial streamline network is smoothed and its dimensions reduced, and redundant branches are removed. The resulting streamline data in the glacier is output, including: Perform point-to-point collision detection of spatial coordinates by traversing the node sequence of the initial candidate streamlines in the initial streamline network; When two initial candidate streamlines originating from different source points intersect in space and the subsequent node sequences completely overlap, the overlapping streamlines are reduced to a unique glacier trunk centerline. Calculate the actual spatial length of the non-overlapping tributaries, obtain the length of the longest midstream line generated in the network of the target glacier and multiply it by a preset pruning ratio coefficient as a pruning threshold. When the total length of the non-overlapping tributaries is lower than the pruning threshold, the non-overlapping tributaries are stripped from the network to remove the redundant branches. The curve fitting algorithm is invoked to perform spatial geometric smoothing on the discrete node sequence, outputting the streamline data in the glacier.

Images