Automatic ridge line identification method and system based on hydrological analysis and density clustering
By combining hydrological analysis and density clustering methods with inverse distance weighted interpolation, priority flood algorithm, D8 method and DBSCAN algorithm, the ridgeline extraction process is optimized, which solves the problem of insufficient accuracy and adaptability of existing methods in complex terrain, and realizes efficient and accurate ridgeline identification.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-12-25
- Publication Date
- 2026-04-07
AI Technical Summary
Existing automatic ridgeline extraction methods are insufficient in terms of accuracy and adaptability. In particular, they are prone to misidentifying hillside points and cliffs in complex terrain. The results contain many noisy points, messy ridgeline branches, and low efficiency of depression filling algorithms, resulting in low accuracy and weak adaptability of ridgeline identification.
A method based on hydrological analysis and density clustering was adopted. Data preprocessing was performed using inverse distance weighted interpolation, depression filling was combined with priority flood algorithm, the D8 method was used to calculate the water flow direction, the DBSCAN algorithm was used for spatial clustering, and continuous processing was performed by combining skeletonization algorithm. Noise points were removed and ridge points with 50% elevation quantiles were retained, and ridge line extraction was optimized.
It significantly improves the adaptability and geometric accuracy of ridgeline extraction, increases processing efficiency, provides more reliable terrain data support, and provides a scientific basis for wind power site selection and wind resource assessment.
Smart Images
Figure CN121811085A_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The application belongs to the technical field of new energy planning, and relates to a mountain ridge line automatic identification method and system. BACKGROUND
[0002] As the core feature of the terrain, the morphology and trend of the mountain ridge directly dominate the regional airflow movement law. The wind turbines of the mountain wind farm are usually arranged along the mountain ridge line in the trend, which can maximize the capture of stable wind energy and achieve land composite utilization. As the terrain skeleton, the mountain ridge line directly determines the lifting effect of the mountain on the airflow and the wind speed amplification, and is the core basis for wind farm site selection. Therefore, the extraction accuracy of the mountain ridge line directly affects the reliability of the core work such as wind resource assessment, wind turbine layout optimization and power generation capacity estimation.
[0003] However, the traditional manual extraction of the mountain ridge is low in efficiency and strong in subjectivity, and cannot meet the rapid analysis needs of large areas and complex terrains. In the existing automatic extraction method of the mountain ridge line, hydrological analysis shows significant advantages due to its unique principle. It describes the constraint and guidance of the terrain to the wind field based on the surface water flow convergence relationship, and can more directly reflect the essential characteristics of the wind field affected by the terrain hydrodynamics in time and space. However, the current mountain ridge line extraction method based on hydrological analysis still has some problems, such as the need for manual setting of the accumulation threshold of the confluence, poor adaptability to different terrains, easy misidentification of mountain ridge points, such as misjudgment of mountain slope points, cliffs and local protruding points as mountain ridge points, no subsequent denoising and continuous processing, many noise points in the result, and branch chaos in the mountain ridge line, which easily occurs in complex mountainous scenes, low efficiency of the depression filling algorithm, slow processing speed, etc., resulting in low accuracy of the mountain ridge line identification and weak self-adaptability. SUMMARY
[0004] To solve the problem of low accuracy of the mountain ridge line identification and weak self-adaptability in the mountain ridge line extraction method based on hydrological analysis in the background art, the application provides a mountain ridge line automatic identification method and system based on hydrological analysis and density clustering, which is not only self-adaptive to different terrains, but also high in accuracy.
[0005] The method of the application comprises: obtaining digital elevation model data of a target area and performing data preprocessing by using an inverse distance weighted interpolation method to obtain preprocessed digital elevation model data; based on the preprocessed digital elevation model data, calculating to obtain positive terrain data of the target area; based on the preprocessed digital elevation model data, performing depression filling processing by using a priority flood algorithm to obtain digital elevation model data without depressions; based on the digital elevation model data without depressions, calculating to obtain water flow direction data of the target area by using a D8 method; Based on the water flow direction data of the target region, the flow accumulation of the target region is calculated; According to the positive terrain data of the target region and the flow accumulation of the target region, the region with zero flow accumulation and positive terrain is identified as a candidate ridge point, and a candidate ridge point set is obtained; The DBSCAN algorithm is used to perform spatial clustering on the candidate ridge point set, and the spatial distance is used to divide the clustering cluster and remove the noise points, and a clustering cluster set without noise points is obtained; The elevation screening is performed on each clustering cluster in the clustering cluster set without noise points, and the ridge points higher than 50% quantile of the cluster elevation are retained as the final ridge points; Based on the final ridge points, a skeleton algorithm is used for continuous processing to obtain a ridge line.
[0006] Further, the acquisition process of the preprocessed digital elevation model data includes: The digital elevation model data of the target region is obtained, denoted as H=[h(x,y)], where h(x,y) is the elevation value at coordinates (x,y), and the digital elevation model data of the target region is an N×N grid, N is a positive integer; the missing grid is filled by using the inverse distance weighted interpolation method, and the grid unit of the interpolation radius and the weight coefficient Determined according to actual needs, the formula of the inverse distance weighted interpolation method is: , Wherein, n represents the number of effective grids around the missing grid, represents the elevation value of the i-th effective grid, represents the distance between the i-th effective grid and the missing grid.
[0007] Further, the calculation process of the positive terrain data of the target region includes: The average elevation of each grid is calculated by using an M×M window, M is a positive integer, and M , Wherein, , Excluding the center grid itself; Then, the positive terrain judgment is performed, if Higher than the average elevation around The terrain mask data =1, all grid points are traversed in turn, and the positive terrain data of the target region .
[0008] Further, the acquisition process of the digital elevation model data without depression includes: The preprocessed digital elevation model data is filled with a priority flood algorithm, and points with lower elevations are added to the queue first. Each time, the point with the lowest elevation is taken out of the queue. If the elevation of a neighbor point is higher than the current point, the elevation of the neighbor point is lowered to the current point, and the neighbor point is added to the queue. The current point is marked as processed, and the queue is iterated until the queue is empty and there are no unprocessed points. Finally, the digital elevation model data without depressions is obtained .
[0009] Further, the calculation process of the water flow direction data of the target region comprises: The water flow direction of the digital elevation model data without depressions is calculated based on the D8 method. The water flow of each grid cell points to the grid with the lowest elevation in the 8 neighbors, and the 8 directions are represented by integers 1-8 (for example, 1 represents north, 2 represents northeast, and so on. The water flow direction matrix The calculation formula is as follows: , wherein, is the grid offset of the 8 directions; All grids are iterated in turn to obtain the water flow direction data of the target region .
[0010] Further, the calculation process of the flow accumulation of the target region comprises: According to the water flow direction data of the target region The flow accumulation of each grid cell is calculated , which is the total number of upstream water flow paths passing through the grid. The calculation formula is as follows: , , wherein, represents the neighbor point set in the direction of all water flow directions.
[0011] Further, the acquisition process of the cluster set of the noise points to be removed comprises: The grid cell with flow accumulation =0 and positive terrain =1 is taken as the candidate ridge point set, and the coordinates of the candidate ridge points are standardized; The DBSCAN algorithm is used, and the neighborhood parameters and the minimum number of samples of each cluster are set; Each ridge point is iterated and judged whether it is a core point. If the number of ridge points within the radius of a certain ridge point is greater than If the ridge point is within the radius of a core point, it is regarded as a boundary point. If the ridge point is within the radius of a core point, it is regarded as a boundary point. If the ridge point is within the radius of a core point, it is regarded as a boundary point. If the ridge point is within the radius of a core point, it is regarded as a boundary point. If the ridge point is within the radius of a core point, it is regarded as a boundary point. If the ridge point is within the radius of a core point, it is regarded as a boundary point. If the ridge point is within the radius of a core point, it is regarded as a boundary point. If the ridge point is within the radius of a core point, it is regarded as a boundary point. If the ridge point is within the radius of a core point, it is regarded as a boundary point. If the ridge point is within the radius of a core point, it is regarded as a boundary point. If the ridge point is within the radius of a core point, it is regarded as a boundary point. If the ridge point is within the radius of a core point, it is regarded as a boundary point.
[0012] Further, the process of obtaining the final ridge point includes: For each cluster, perform elevation filtering, and retain the ridge points higher than the 50th percentile of the cluster elevation as the final ridge point, the 50th percentile of the elevation of the i-th cluster is the elevation value that makes the cumulative probability equal to 50%, and the cumulative probability of the i-th cluster The specific calculation formula is as follows: , wherein, is the quantile, if then is the 50th percentile, represents the number of points belonging to the cluster represents the number of points in the cluster whose elevation value is less than or equal to h.
[0013] Further, the process of obtaining the ridge line includes: An initial ridge line segment of single grid width is extracted by using a skeletonization algorithm, then gaps are filled by using a closing operation, burrs are removed by using an opening operation, and finally segment connection is performed, that is, if the distance between the end points of adjacent ridge segments is less than or equal to 2 grids and the included angle of the strike is less than or equal to 25°, the connection is formed into a continuous ridge line, and the specific steps are as follows: First step, extract the center line of the original ridge area, and judge whether the edge grid of the ridge area meets the deletion condition, the deletion condition is as follows: taking each grid as the center, take a 3x3 grid range, and count the number of grids (denoted as N) belonging to the ridge area, if N 6, it is explained that the removable condition is met; the number of grids directly adjacent to each grid and belonging to the ridge area is recorded as C, if C=1 or C=2, it is deleted; the sequential traversal p2, p3,…, p9, p2 in the 8-neighborhood of each grid p1 is calculated, the number of transitions from non-ridge (0) to ridge (1) is recorded as A (p1), which is exactly equal to 1, then it is deleted; then the edge grid that meets the condition is deleted step by step, the center skeleton of the ridge area is reserved, and the line formed by all the remaining ridge area grids is single-grid width; In the second step, after the skeleton operation of the first step, small gaps of 1-2 grids appearing in the ridge line are filled by a closing operation, a preset structural element including a 3*3 size rectangle or a circular template is used for an inflation operation on the image, the single-grid line is expanded outward, the small gap of 1-2 grids is filled, and the lines on both sides of the gap are connected; In the third step, after the skeleton operation of the second step, for the small bifurcation of 1-2 grids appearing in the ridge line, a preset structural element including a 3*3 size rectangle or a circular template is used for an erosion operation on the image after the inflation operation of the second step, the expanded line is retracted to the original width, the single-grid characteristic is restored, and the filled continuous part is reserved; then a preset structural element including a 3*3 size rectangle or a circular template is used for an erosion operation, and small burrs with a length of not more than 2 grids are deleted; finally, a preset structural element including a 3*3 size rectangle or a circular template is used for an inflation operation, so that the main stem of the ridge line restores the single-grid width; In the fourth step, for the gap of >2 grids appearing in the ridge line, the endpoints of all ridge line segments are extracted, wherein there is only one adjacent grid in the single-grid line, the grid distance and the segment direction angle between each pair of endpoints are calculated one by one, the grid distance is the Euclidean distance, and the segment direction angle is the extension direction angle of the segment at the endpoint, if the distance between two points is not more than 2 grids and the direction angle is not more than 25°, the two segments are connected as a continuous line segment.
[0014] Based on the above method, the application provides a ridge line automatic identification system based on hydrological analysis and density clustering, which comprises a digital elevation model data preprocessing module, a positive terrain data acquisition module, a digital elevation model data acquisition module without depression, a flow direction data acquisition module, a flow accumulation quantity acquisition module, a candidate ridge point set acquisition module, a cluster set acquisition module for removing noise points, a final ridge point acquisition module, and a ridge line acquisition module.
[0015] The digital elevation model data preprocessing module acquires digital elevation model data of a target region and performs data preprocessing by using an inverse distance weighted interpolation method to obtain preprocessed digital elevation model data.
[0016] The positive terrain data acquisition module obtains the positive terrain data of the target region by calculation based on the preprocessed digital elevation model data.
[0017] The digital elevation model data without depression acquisition module obtains the digital elevation model data without depression by filling depression processing using a priority flood algorithm based on the preprocessed digital elevation model data.
[0018] The flow direction data acquisition module obtains the flow direction data of the target region by calculation using a D8 method based on the digital elevation model data without depression.
[0019] The flow accumulation quantity acquisition module calculates the flow accumulation quantity of the target region based on the flow direction data of the target region.
[0020] The candidate ridge point set acquisition module identifies a region with zero flow accumulation quantity and positive terrain as a candidate ridge point based on the positive terrain data of the target region and the flow accumulation quantity of the target region, and obtains a candidate ridge point set.
[0021] The cluster set with removed noise points acquisition module performs spatial clustering on the candidate ridge point set using a DBSCAN algorithm, divides cluster clusters according to spatial distance, removes noise points, and obtains a cluster set with removed noise points.
[0022] The final ridge point acquisition module performs elevation screening on each cluster in the cluster set with removed noise points, and retains a ridge point higher than 50% quantile of the cluster elevation as a final ridge point.
[0023] The ridge line acquisition module obtains a ridge line by performing continuous processing using a skeletonization algorithm based on the final ridge point.
[0024] Compared with the prior art, the present application has the following beneficial effects: (1) Adaptive terrain processing capability is improved: the present application uses "positive terrain recognition + zero flow accumulation" as the initial criterion for candidate ridge points, replacing the traditional method of relying on manual setting of the flow accumulation threshold, so that the preliminary identification of ridge points can be automatically completed according to the hydrological and morphological characteristics of the terrain itself, significantly improving the adaptive ability of the method to different terrains.
[0025] (2) Ridges line denoising and continuity optimization: In view of the problems of many noise points and branch disorder in the results of the existing method, the DBSCAN spatial clustering algorithm is introduced to cluster and analyze the candidate ridge points, so that the spatial isolated noise points such as misrecognized hillside points and local convex points are effectively removed; then, the points higher than the 50% quantile of the elevation in the cluster are reserved for purification, and the skeletonization algorithm is combined for continuity processing, so that the ridge line with good connectivity, stable geometric shape and clear main stem is finally obtained, and the defect of the traditional result of branch and broken is overcome; (3) Improvement of depression processing efficiency: The priority flood algorithm is used for depression filling processing of digital elevation model data, compared with the traditional iterative filling algorithm, the filling order is optimized by the priority queue, the calculation efficiency is significantly improved when processing large-scale digital elevation model data, and a high-quality depression-free digital elevation model basis is provided for subsequent hydrological analysis, so that the operation speed of the whole process is ensured; (4) Enhancement of comprehensive accuracy and reliability: by fusing hydrological analysis (reflecting the nature of water dynamics of air flow guided by terrain), spatial clustering, elevation screening, morphological processing and other multi-step optimization, the present application constructs a complete technical scheme of "preprocessing-positive terrain identification-concentration analysis-cluster denoising-elevation screening-continuity optimization", which systematically improves the geometric accuracy and topological rationality of ridge line extraction, provides a more reliable terrain skeleton data basis for wind resource assessment and wind turbine site selection and layout, and finally enhances the scientificity of wind power project planning and power generation estimation.
[0026] In summary, the present application automatically selects ridge points by positive terrain and zero convergence accumulation, effectively removes noise and short branches by combining DBSCAN clustering and elevation quantile purification, and efficiently fills depression and optimizes continuity by using the priority flood algorithm and the skeletonization algorithm, which significantly improves the adaptability, geometric accuracy and processing efficiency of ridge line extraction, and provides more reliable terrain data support for wind power site selection. BRIEF DESCRIPTION OF DRAWINGS
[0027] Figure 1 The method flowchart of the present application.
[0028] Figure 2 The system architecture diagram of the present application.
[0029] Figure 3 The original DEM graph of the embodiment.
[0030] Figure 4 The cumulative quantity comparison graph of DEM before and after depression filling of the embodiment.
[0031] Figure 5 The ridge line identification graph of the embodiment. DETAILED DESCRIPTION
[0032] In order to make the technical problems, technical solutions and beneficial effects of the present application clearer, the present application will be further described in detail below in conjunction with the drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and not to limit the present application.
[0033] The first aspect of the present application is based on a method for automatically identifying ridge lines of a mountain range based on hydrological analysis and density clustering, and a flowchart is shown as Figure 1 The specific steps are described as follows.
[0034] Firstly, the digital elevation model data of the target area is obtained and data preprocessing is performed by using the inverse distance weighted interpolation method to obtain the preprocessed digital elevation model data.
[0035] Digital elevation model (DEM) is a technology for digital simulation of ground terrain based on limited terrain elevation data. With the maturity and popularity of this technology, automatic extraction of terrain features has become a core requirement in the fields of geographic information science, environmental science and engineering construction.
[0036] Inverse distance weighted interpolation method: based on the principle that "the closer the known points to the unknown point, the greater the influence", the missing elevation value is filled by calculating the distance weight of the missing grid and the surrounding effective grid.
[0037] Specifically, the digital elevation model data of the target area is obtained, denoted as H=[h(x,y)], where h(x,y) is the elevation value at coordinates (x,y), and the digital elevation model data of the target area is an N×N grid, N is a positive integer; the missing grid is filled by using the inverse distance weighted interpolation method, the grid unit of the interpolation radius and the weight coefficient Determined according to actual requirements, the formula of the inverse distance weighted interpolation method is: , Wherein, n represents the number of effective grids around the missing grid, represents the elevation value of the i-th effective grid, represents the distance between the i-th effective grid and the missing grid.
[0038] Then, based on the preprocessed digital elevation model data, the positive terrain data of the target area is calculated.
[0039] Positive terrain refers to the terrain that is higher than the adjacent area, such as plateau and mountainous area.
[0040] Specifically, the average elevation of each grid is calculated by using a M×M window, M is a positive integer, and M , Wherein, , This indicates that the center grid cell itself has been excluded; Then perform a positive terrain assessment; if... Higher than the surrounding average elevation Then the terrain mask data =1, sequentially traverse all grid points to obtain the positive terrain data of the target area. .
[0041] Meanwhile, based on the preprocessed digital elevation model data, a priority flood algorithm is used to fill depressions and obtain digital elevation model data without depressions.
[0042] A depression is a closed low-lying area formed by a grid point and its connected neighborhood in a DEM. Its elevation is lower than the height of all surrounding accessible exit grids, preventing water from naturally flowing out of the boundary.
[0043] Priority-Flood algorithm: simulates the process of water starting from the edge or lowest point of the terrain and naturally flooding the entire area. By using a priority queue, it can accurately locate the points that need to be processed, and reduce the surrounding points so that they can naturally connect to the edge. The algorithm can process multiple depressions at once without repeated iterations.
[0044] Specifically, a priority flooding algorithm is used to fill depressions in the preprocessed digital elevation model data. Points with lower elevations are added to a queue first. Each time, the point with the lowest current elevation is retrieved from the queue. If a neighboring point has a higher elevation, its elevation is lowered to the current point's elevation, and the neighboring point is added to the queue. The current point is marked as processed. This process is repeated until the queue is empty and there are no unprocessed points, ultimately resulting in digital elevation model data without depressions. .
[0045] Next, based on the digital elevation model data without depressions, the D8 method was used to calculate the water flow direction data of the target area.
[0046] Specifically, the water flow direction in the D8 method is calculated based on the digital elevation model data without depressions. The water flow in each grid cell points to the grid cell with the lowest elevation in its eight neighborhoods. The eight directions are represented by integers from 1 to 8 (e.g., 1 represents due north, 2 represents northeast, and so on). The water flow direction matrix is then used. The calculation formula is as follows: , in, This refers to the grid offset in 8 directions; By iterating through all the grids, the water flow direction data for the target area can be obtained. .
[0047] Next, based on the water flow direction data of the target region, the flow accumulation amount of the target region is calculated.
[0048] Specifically, according to the water flow direction data of the target region The flow accumulation amount of each grid cell is calculated , that is, the total number of upstream water flow paths flowing through the grid, and the calculation formula is as follows: , , Among them, represents the set of neighborhood points whose water flow directions point to .
[0049] Next, according to the positive terrain data of the target region and the flow accumulation amount of the target region, the region with zero flow accumulation amount and positive terrain is identified as a candidate ridge point, and a set of candidate ridge points is obtained.
[0050] Next, the DBSCAN algorithm is used to perform spatial clustering on the set of candidate ridge points, and the clustering clusters and noise points are removed according to the spatial distance, and a set of clustering clusters with noise points removed is obtained.
[0051] Specifically, the grid cell with flow accumulation amount =0 and positive terrain =1 is taken as the set of candidate ridge points, and the coordinates of the candidate ridge points are standardized; The DBSCAN algorithm is used, and the neighborhood parameters and the minimum number of samples of each cluster are set; Each ridge point is iteratively judged to be a core point, if the number of ridge points within the radius of a certain ridge point is greater than , then the ridge point is considered as a core point; Each ridge point is iteratively judged to be a noise point, if the number of ridge points within the radius of a certain ridge point is less than , but the ridge point is within the radius of a certain core point, then the ridge point is considered as a boundary point; if the number of ridge points within the radius of a certain ridge point is less than , and the ridge point is not within the radius of a certain core point, then the ridge point is a noise point; The classification result of each ridge point is output, that is, the obtained clustering cluster set , where M is the number of clusters, that is, the number of core points, represents all the ridge points belonging to the cluster where the first core point is located, that is, the boundary point set of the first core point; The cluster set with noise points removed is obtained through the above steps.
[0052] Next, elevation screening is performed on each cluster in the cluster set with noise points removed, and ridge points higher than the 50th percentile of the cluster elevation are retained as final ridge points.
[0053] Specifically, elevation screening is performed on each cluster, and ridge points higher than the 50th percentile of the cluster elevation are retained as final ridge points, and the 50th percentile of the elevation of the i-th cluster is , which refers to the elevation value that makes the cumulative probability equal to 50%, and the cumulative probability of the i-th cluster is The specific calculation formula is as follows: , wherein, is the percentile, if , then is the 50th percentile, represents the number of points belonging to the cluster , and represents the number of points in the cluster whose elevation values are less than or equal to h.
[0054] Finally, based on the final ridge points, a skeletonization algorithm is used for continuous processing to obtain a ridge line.
[0055] Specifically, a skeletonization algorithm is used to extract single-grid-width initial ridge line segments, then gap filling is performed using a closing operation, burr removal is performed using an opening operation, and finally segment connection is performed, that is, if the distance between the end points of adjacent ridge segments is ≤2 grids and the strike angle is ≤25°, then connection is performed to form a continuous ridge line, and the specific steps are as follows: First step, the center line of the original ridge area is extracted, and it is judged whether the edge grid of the ridge area meets the deletion condition, and the deletion condition is as follows: taking each grid as the center, a 3x3 grid range is taken, and the number of grids belonging to the ridge area (denoted as N) is counted, if N 6, it is indicated that the deletion condition is met; the number of grids directly adjacent to each grid and belonging to the ridge area is denoted as C, if C=1 or C=2, then deletion is performed; the sequential traversal p2, p3,…, p9, p2 in the 8-neighborhood of each grid p1 is calculated, and the number of transitions from non-ridge (0) to ridge (1) is denoted as A(p1), which is exactly equal to 1, then deletion is performed; then the edge grids meeting the condition are gradually deleted, and the center skeleton of the ridge area is retained, until all the ridge area grids remaining form a line with a single grid width; Second step, after the skeletonization operation of the first step, fill in the small gap of 1-2 grids of the ridge line by closing operation, use the preset structural element including 3*3 size rectangle or circular template to perform inflation operation on the image, let the single grid line expand outward, fill in the small gap of 1-2 grids, connect the lines on both sides of the gap; Third step, after the skeletonization operation of the second step, for the small bifurcation of 1-2 grids of the ridge line, use the preset structural element including 3*3 size rectangle or circular template to perform erosion operation on the image after the inflation operation of the second step, let the expanded line shrink to the original width, restore the single grid characteristic, and at the same time keep the continuous part after filling; then use the preset structural element including 3*3 size rectangle or circular template to perform erosion operation, delete the small burr with a length of not more than 2 grids; finally, use the preset structural element including 3*3 size rectangle or circular template to perform inflation operation, so that the main stem of the ridge line restores the single grid width; Fourth step, for the gap of more than 2 grids of the ridge line, extract the endpoints of all ridge line segments, wherein there is only one adjacent grid point in the single grid line, calculate the grid distance and segment direction angle between each pair of endpoints one by one, the grid distance is the Euclidean distance, and the segment direction angle is the extension direction angle of the segment where the endpoint is located, if the distance between two points is not more than 2 grids and the direction angle is not more than 25°, the two segments are connected as a continuous line segment.
[0056] It should be pointed out that skeletonization means iteratively deleting image boundary grids to finally obtain a skeleton structure with a single grid width, while keeping the topological features of the original image. Morphological operation is a shape transformation operation based on grid, closing operation (inflation + erosion) is used to fill in the gap, and opening operation (erosion + inflation) is used to remove burr. Inflation: let the baseline in the grid "expand", which can fill in the small cavity in the baseline and connect the adjacent broken segments; Erosion: let the ridge line in the grid "shrink", which can remove small noise points, burrs, or disconnect the fine connection in the foreground.
[0057] The second aspect of the present application provides a ridge line automatic identification system based on hydrological analysis and density clustering, and the architecture diagram is as shown in Figure 2 The system is composed of a digital elevation model data preprocessing module, a positive terrain data acquisition module, a digital elevation model data acquisition module without depression, a flow direction data acquisition module, a flow accumulation quantity acquisition module, a candidate ridge point set acquisition module, a cluster set acquisition module for removing noise points, a final ridge point acquisition module, and a ridge line acquisition module.
[0058] The preprocessing module of digital elevation model data acquires digital elevation model data of a target region and performs data preprocessing on the digital elevation model data by using an inverse distance weighted interpolation method to obtain preprocessed digital elevation model data.
[0059] The positive terrain data acquisition module acquires positive terrain data of the target region by calculation based on the preprocessed digital elevation model data.
[0060] The digital elevation model data without depression acquisition module acquires digital elevation model data without depression by performing depression filling processing on the preprocessed digital elevation model data by using a priority flood algorithm.
[0061] The flow direction data acquisition module acquires flow direction data of the target region by calculation based on the digital elevation model data without depression by using a D8 method.
[0062] The flow accumulation quantity acquisition module calculates flow accumulation quantity of the target region based on the flow direction data of the target region.
[0063] The candidate ridge point set acquisition module identifies a region with zero flow accumulation quantity and positive terrain as a candidate ridge point based on the positive terrain data of the target region and the flow accumulation quantity of the target region to obtain a candidate ridge point set.
[0064] The cluster set with removed noise points acquisition module performs spatial clustering on the candidate ridge point set by using a DBSCAN algorithm, divides cluster sets according to spatial distance, and removes noise points to obtain a cluster set with removed noise points.
[0065] The final ridge point acquisition module performs elevation screening on each cluster in the cluster set with removed noise points, retains a ridge point higher than 50% quantile of the elevation of the cluster as a final ridge point.
[0066] The ridge line acquisition module performs continuous processing on the final ridge point by using a skeletonization algorithm to obtain a ridge line.
[0067] The specific implementation of each module in the system is consistent with the content described in the above method, and will not be repeated here.
[0068] The third aspect of the present application provides an electronic device, comprising a memory and a processor, which are connected in communication with each other, the memory stores computer instructions, and the processor implements the above-mentioned hydrological analysis and density clustering-based ridge line automatic identification method by executing the computer instructions.
[0069] The fourth aspect of the present application provides a computer readable storage medium, the computer readable storage medium stores a computer program, the computer program is executed by a processor to realize the automatic ridge line identification method based on hydrological analysis and density clustering as described above.
[0070] Embodiments Taking DEM data of a mountainous area (longitude range 121.25°E-121.75°E, latitude range 45.5°N-45.9°N) as the research object, the DEM data has a resolution of 30 m, and the original DEM graph is as shown in Figure 3 The ridge line is the highest connecting line along the uplift of the mountain, that is, the line formed by the continuous top edge of the mountain peak and the hill. This line is usually represented as a series of contour lines formed by connected mountain peak points or highlands on the map. The automatic ridge line identification method of the present application is as follows: Step 1, DEM data acquisition and preprocessing: first, acquire the DEM data of the target area, denoted as , wherein, is the elevation value (unit: m) at the coordinate , and the DEM data size in this embodiment is 600*600 grids; then, the inverse distance weighted interpolation method (IDW) is used to fill in the missing grids, the interpolation radius is set to 5 grid units, and the weight coefficient is 1 / d, and the inverse distance weighted interpolation method formula is: , wherein n represents the number of valid grids around the missing grid, represents the elevation value of the i-th valid grid, represents the distance between the i-th valid grid and the missing grid.
[0071] All missing grid points are traversed in turn to obtain the filled elevation data .
[0072] Step 2, positive terrain identification: first, calculate the average elevation around, mainly using a 50*50 window to calculate the average elevation of each grid, and the specific formula is: , wherein, , represents the center grid itself.
[0073] Then, the positive terrain is judged. If is higher than the average elevation around , the terrain mask data =1, and all grid points are traversed in turn to obtain ; Step 3, filling the pit process: in order to ensure the accuracy of subsequent water flow direction calculation, the priority flood algorithm is used to fill the pit of DEM data, which preferentially adds the point with lower elevation to the queue, takes out the point with the lowest elevation from the queue every time, if the elevation of the neighbor point is higher than that of the current point, the elevation of the neighbor point is reduced to the elevation of the current point, and the neighbor point is added to the queue, the current point is marked as processed, and the queue is traversed in turn until the queue is empty and there is no unprocessed point, and finally the DEM data without pit is obtained .
[0074] The comparison chart of cumulative amount before and after DEM filling is shown in Figure 4 , Step 4, water flow direction calculation: based on D8 method, the water flow direction of the DEM without pit is calculated, the water flow of each grid cell points to the grid with the lowest elevation in the 8 neighbors, which is represented by an integer from 1 to 8 (such as 1 represents north, 2 represents northeast, and so on), and the water flow direction matrix The calculation formula is as follows: Among them, is the grid offset of 8 directions.
[0075] All grids are traversed in turn to obtain the final water flow direction data .
[0076] Step 5, flow accumulation calculation: according to the water flow direction matrix The flow accumulation of each grid cell is calculated , that is, the total number of upstream water flow paths through the grid, and the calculation formula is as follows: , , Among them, represents the neighbor point set of all water flow directions .
[0077] Step 6, preliminary extraction of candidate ridge points: the grid cell with flow accumulation =0 and positive terrain =1 is identified as the candidate ridge point, and the candidate ridge point set is obtained.
[0078] Step 7, spatial clustering analysis: DBSCAN algorithm is used to perform spatial clustering on the candidate ridge points, and the clustering clusters are divided according to the spatial distance, and the noise points are removed; DBSCAN algorithm is a clustering method based on density scanning, in simple terms, DBSCAN algorithm is to start from a sample, find the neighbors, and continuously include the nearby points in the same category, and the specific steps are as follows: The first step is to standardize the coordinates of the candidate ridge points; Step 2: Set neighborhood parameters Minimum number of samples per cluster ; The third step is to iterate through and determine whether each ridge point is a core point. If a certain ridge point is a core point... The ridge point within the radius is greater than If so, then the ridge point is considered the core point; Step 4: Traverse and determine whether each ridge point is a noise point. If a certain ridge point... The ridge point within the radius is smaller than However, the ridge point is located at a certain core point. If a ridge point is within the radius, then that ridge point is considered a boundary point; if a certain ridge point The ridge point within the radius is smaller than Moreover, the ridge point is not located at a core point. If the ridge point is within the radius, then that ridge point is a noise point; Step 5: Output the classification results for each ridge point, i.e., the resulting set of clusters. Where M is the number of clusters (number of core points). This indicates that it belongs to the first core point. The set of all ridge points within the radius, i.e., the boundary points of the first core point.
[0079] The above steps yield the final set of clusters with noise points removed.
[0080] Step 8: Elevation Screening and Purification: Perform elevation screening for each cluster, retaining ridge points with elevations higher than the 50th percentile of the cluster's elevation as the final ridge points. The 50th percentile of the elevation for the i-th cluster is then used. The elevation value that makes the cumulative probability equal to 50% is the cumulative probability of the i-th cluster. The specific calculation formula is as follows: , in, Quantiles, if ,So It is the 50th percentile. Indicates belonging to a cluster The number of points, Indicates in cluster The number of points whose elevation value is less than or equal to h.
[0081] Step 9, continuous processing: the main skeleton (Zhang-Suen algorithm) algorithm is used to extract single grid width initial ridge line segment, and then the gap is filled by using the closing operation, the burr is removed by using the opening operation, and finally the segment connection is performed, that is, if the distance between the end points of adjacent ridge segments is less than or equal to 2 grids and the strike angle is less than or equal to 25°, the continuous ridge line is formed, and the specific principle and steps are as follows: First step, the original ridge is a large bright area (such as a yellow highland) on the map, and we need to extract the center line of this area, not the entire area. Single grid width is the basis for subsequent processing, and different line thicknesses will cause errors in subsequent gap filling and burr removal operations. The specific method is to repeatedly check the edge grid of the ridge area to determine whether it meets the deletable condition. The deletable condition is as follows: ① Take 3×3 grid range centered on each grid, count the number of grids belonging to the ridge area (denoted as N), if N 6, it means that the deletable condition is met; ② The number of grids directly adjacent to each grid and belonging to the ridge area is denoted as C, if C=1 or C=2, then delete; ③ Calculate the number of transitions from non-ridge (0) to ridge (1) when sequentially traversing p2, p3, …, p9, p2 in the 8-neighborhood of each grid p1, denoted as A(p1), which is exactly equal to 1, then delete.
[0082] Then, the edge grids that meet the condition are deleted step by step, and the center skeleton of the ridge area is retained, until the line formed by all the remaining ridge area grids is single grid width; Second step, after the skeletonization operation of the first step, the ridge line may have 1-2 grid small gaps, which is because the original terrain has small depressions, causing the line not to be connected during skeletonization. The closing operation can fill such small gaps. The specific method is to use a pre-set structure element (such as a 3×3 size rectangle or circular template) to perform dilation operation on the image, that is, to expand the single grid line outward to fill the 1-2 grid small gap, and connect the lines on both sides of the gap; Third step, after the skeleton operation of the first step, the ridge line may have 1-2 small branches (i.e. burrs) of grids, which is because the original terrain has a small protrusion, which is mistaken as a part of the ridge. These burrs will affect the subsequent segment connection and must be removed. The specific method is: first, use the same structural element to perform the erosion operation on the expanded image, so that the expanded line shrinks to the original width, restores the single grid characteristic, and at the same time preserves the filled continuous part; then use the same structural element to perform the erosion operation to delete the small burrs with a length of not more than 2 grids (the burrs are thin, and after erosion, they will disappear directly); finally, use the same structural element to perform the expansion operation to restore the single grid width of the ridge line trunk, and ensure that the main structure is not damaged; Fourth step, the closing operation can only fill in 1-2 small gaps of grids, and if the gap is larger (such as 3 grids), it cannot be filled; but these scattered segments may originally belong to the same ridge line, only with a slightly larger depression in the middle, so further connection is needed. The specific method is to extract the endpoints of all ridge line segments (only one adjacent grid point exists in a single grid line), and calculate the grid distance (Euclidean distance) and segment direction angle (the extension direction angle of the segment where the endpoint is located) between each pair of endpoints. If the distance between two points is not more than 2 grids and the direction angle is not more than 25°, then the two segments are connected as a continuous line segment.
[0083] Step 10, output the final result: the ridge line recognition map is as shown in Figure 5 The black (0 value) is a non-ridge point, and the white (1000 value) is a ridge point.
[0084] Those skilled in the art should understand that the embodiments of the present application can be provided as a method, a system, or a computer program product. Therefore, the present application can adopt a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer usable program code. The solutions in the embodiments of the present application can be implemented in various computer languages, such as object-oriented programming languages Java, C++, Python, and interpreted scripting language JavaScript.
[0085] The computer program instructions can also be loaded onto a computer or other programmable data processing electronic device to cause a series of operational steps to be performed on the computer or other programmable data processing electronic device to produce a computer-implemented process such that the instructions which execute on the computer or other programmable data processing electronic device provide steps for implementing the functions specified in the flowchart Figure 1 one or more flows and / or blocks Figure 1 one or more flows and / or blocks
[0086] The computer program instructions can also be loaded onto a computer or other programmable data processing electronic device to cause a series of operational steps to be performed on the computer or other programmable data processing electronic device to produce a computer-implemented process such that the instructions which execute on the computer or other programmable data processing electronic device provide steps for implementing the functions specified in the flowchart Figure 1 one or more flows and / or blocks Figure 1 one or more flows and / or blocks
[0087] The computer program instructions can also be loaded onto a computer or other programmable data processing electronic device to cause a series of operational steps to be performed on the computer or other programmable data processing electronic device to produce a computer-implemented process such that the instructions which execute on the computer or other programmable data processing electronic device provide steps for implementing the functions specified in the flowchart Figure 1 one or more flows and / or blocks Figure 1 one or more flows and / or blocks
[0088] While the preferred embodiments of the application have been described, additional variations and modifications can be employed by those skilled in the art. Therefore, the appended claims intend to cover all such modifications and variations as fall within the true spirit and scope of the application.
[0089] Obviously, numerous modifications and variations of the present application are possible in light of the above teachings. It is therefore to be understood that within the scope of the appended claims and their equivalents, the application can be practiced otherwise than as specifically described.
Claims
1. An automatic ridgeline identification method based on hydrological analysis and density clustering, characterized in that, include: The digital elevation model data of the target area is acquired and the data is preprocessed using the inverse distance weighted interpolation method to obtain the preprocessed digital elevation model data. Based on the preprocessed digital elevation model data, the positive terrain data of the target area is calculated. Based on the preprocessed digital elevation model data, a priority flood algorithm is used to fill depressions and obtain digital elevation model data without depressions. Based on the digital elevation model data without depressions, the D8 method was used to calculate the water flow direction data of the target area. Calculate the cumulative flow in the target area based on the water flow direction data of the target area; Based on the positive topographic data and the cumulative flow of the target area, areas with zero cumulative flow and positive topography are identified as candidate ridge points, thus obtaining a set of candidate ridge points; The DBSCAN algorithm is used to perform spatial clustering on the candidate ridge point set. The clusters are divided according to spatial distance and noise points are removed to obtain a set of clusters with noise points removed. For each cluster in the set of clusters where noise points have been removed, an elevation screening is performed, and ridge points that are higher than the 50th percentile of the cluster's elevation are retained as the final ridge points; Based on the final ridge point, a skeletonization algorithm is used for continuous processing to obtain the ridgeline.
2. The automatic ridgeline identification method based on hydrological analysis and density clustering according to claim 1, characterized in that: The process of acquiring the preprocessed digital elevation model data includes: Obtain the digital elevation model (DEM) data for the target area, denoted as H=[h(x,y)], where h(x,y) is the elevation value at coordinates (x,y). The DEM data for the target area consists of N×N grid cells, where N is a positive integer. Missing grid cells are filled using inverse distance weighted interpolation, with the interpolation radius defined by the grid cell and the weighting coefficient. The formula for inverse distance weighted interpolation is determined based on actual needs: ， Where n represents the number of valid rasters surrounding the missing raster. This represents the elevation value of the i-th valid raster. This represents the distance between the i-th valid grid cell and the missing grid cell.
3. The automatic ridgeline identification method based on hydrological analysis and density clustering according to claim 2, characterized in that: The calculation process for the positive terrain data of the target area includes: The average elevation of each grid cell is calculated using an M×M window, where M is a positive integer and M < N. The formula is as follows: ， in, , This indicates that the center grid cell itself has been excluded; Then perform a positive terrain assessment; if... Higher than the surrounding average elevation Then the terrain mask data =1, sequentially traverse all grid points to obtain the positive terrain data of the target area. .
4. The automatic ridgeline identification method based on hydrological analysis and density clustering according to claim 3, characterized in that: The process of acquiring the depression-free digital elevation model data includes: A priority flooding algorithm is used to fill depressions in the preprocessed digital elevation model (DEM) data. Points with lower elevations are added to a queue first. Each time, the point with the lowest elevation is retrieved from the queue. If a neighboring point has a higher elevation, its elevation is lowered to the current point's elevation, and the neighboring point is added to the queue. The current point is marked as processed. This process is repeated until the queue is empty and there are no unprocessed points, ultimately resulting in DEM data without depressions. .
5. The automatic ridgeline identification method based on hydrological analysis and density clustering according to claim 4, characterized in that: The calculation process for the water flow direction data of the target area includes: The D8 method is used to calculate the water flow direction in a depression-free digital elevation model (DEM). The water flow in each grid cell points to the grid cell with the lowest elevation in its eight neighborhoods. Eight directions are represented by integers from 1 to 8 (e.g., 1 represents true north, 2 represents northeast, and so on). The water flow direction matrix... The calculation formula is as follows: ， in, This refers to the grid offset in 8 directions; By iterating through all the grids, the water flow direction data for the target area can be obtained. .
6. The automatic ridgeline identification method based on hydrological analysis and density clustering according to claim 5, characterized in that: The calculation process for the cumulative traffic volume in the target area includes: Based on the water flow direction data of the target area Calculate the cumulative flow for each grid cell. This refers to the total number of upstream water flow paths passing through this grid, calculated using the following formula: ， ， in, Indicates the direction of all water flow. The set of neighborhood points.
7. The automatic ridgeline identification method based on hydrological analysis and density clustering according to claim 6, characterized in that: The process of obtaining the cluster set of noise-removed points includes: Accumulated traffic =0 and positive terrain The grid cells with a value of 1 are used as the set of candidate ridge points, and the coordinates of the candidate ridge points are standardized. The DBSCAN algorithm is used to set the neighborhood parameters. and the minimum number of samples for each cluster ; Iterate through the ridge points and determine if each ridge point is a core point. If a certain ridge point is a core point... The ridge point within the radius is greater than If so, then the ridge point is considered the core point; Iterate through each ridge point and determine if it is a noise point. If a certain ridge point is a noise point... The ridge point within the radius is smaller than However, the ridge point is located at a certain core point. If a ridge point is within the radius, then that ridge point is considered a boundary point; if a certain ridge point The ridge point within the radius is smaller than Moreover, the ridge point is not located at a core point. If the ridge point is within the radius, then that ridge point is a noise point; Output the classification result for each ridge point, i.e., the resulting set of clusters. Where M is the number of clusters, i.e., the number of core points. This represents the set of all ridge points belonging to the cluster containing the first core point, i.e., the set of boundary points of the first core point; The above steps yield a set of clusters with noise points removed.
8. The automatic ridgeline identification method based on hydrological analysis and density clustering according to claim 7, characterized in that: The process of obtaining the final ridge point includes: For each cluster, perform elevation screening, retaining ridge points higher than the 50th percentile of the cluster's elevation as the final ridge points. The 50th percentile of the elevation of the i-th cluster is then used. The elevation value that makes the cumulative probability equal to 50% is the cumulative probability of the i-th cluster. The specific calculation formula is as follows: ， in, Quantiles, if ,So It is the 50th percentile. Indicates belonging to a cluster The number of points, Indicates in cluster The number of points whose elevation value is less than or equal to h.
9. The automatic ridgeline identification method based on hydrological analysis and density clustering according to claim 8, characterized in that: The process of obtaining the ridgeline includes: The skeletonization algorithm is used to extract initial ridgeline segments with a single grid width. Then, the closing operation is used to fill the gaps, the opening operation is used to remove burrs, and finally, the segments are connected. That is, if the distance between the endpoints of adjacent ridgeline segments is ≤2 grids and the directional angle is ≤25°, they are connected to form a continuous ridgeline. The specific steps are as follows: Step 1: Extract the centerline of the original ridge area and determine whether the edge grid cells of the ridge area meet the deletion criteria. The deletion criteria are as follows: Taking each grid cell as the center, take a 3×3 grid area and count the number of grid cells belonging to the ridge area (denoted as N). If... N 6. This indicates that the conditions for deletion are met; the number of grids that are directly adjacent to each grid and belong to the ridge area is denoted as C. If C=1 or C=2, then delete; calculate the number of times non-ridge (0) to ridge (1) is transformed in the 8-neighborhood of each grid p1 when sequentially traversing p2, p3, …, p9, p2, denoted as A(p1). If this number is exactly equal to 1, then delete; then gradually delete the edge grids that meet the conditions, retain the central skeleton of the ridge area, until all the lines formed by the remaining ridge area grids are of single grid width; The second step, after the skeletonization operation in the first step, fills the small gap of 1-2 grids in the ridge line by closing the operation. Then, use the preset structural elements, including 3×3 rectangular or circular templates, to perform a dilation operation on the image, so that the single grid lines expand outward and fill the small gap of 1-2 grids, so that the lines on both sides of the gap are connected. In the third step, after the skeletonization operation in the second step, for the small forks of 1-2 grids that appear on the ridgeline, use a preset structural element including a 3×3 rectangular or circular template to perform an erosion operation on the image after the dilation operation in the second step, so that the expanded lines shrink back to their original width, restoring the single-grid characteristics, while retaining the continuous part after filling; then use a preset structural element including a 3×3 rectangular or circular template to perform an erosion operation to delete small burrs with a length of no more than 2 grids; finally, use a preset structural element including a 3×3 rectangular or circular template to perform a dilation operation, so that the main trunk of the ridgeline restores the single-grid width. Step 4: For ridgelines with gaps of more than 2 grids, extract the endpoints of all ridgeline segments. In a single grid line, there is only one point in an adjacent grid. Calculate the grid distance and segment orientation angle between each pair of endpoints. The grid distance is the Euclidean distance, and the segment orientation angle is the angle of the extension direction of the segment where the endpoint is located. If the distance between two points does not exceed 2 grids and the orientation angle does not exceed 25°, then connect the two segments into a continuous line segment.
10. An automatic ridgeline identification system based on hydrological analysis and density clustering that implements the method of any one of claims 1-9, characterized in that: It includes modules for preprocessing digital elevation model data, acquiring positive terrain data, acquiring digital elevation model data without depressions, acquiring water flow direction data, acquiring cumulative flow, acquiring candidate ridge point sets, acquiring cluster sets of noise-removed points, acquiring final ridge points, and acquiring ridge lines. The digital elevation model data preprocessing module acquires the digital elevation model data of the target area and performs data preprocessing using the inverse distance weighted interpolation method to obtain the preprocessed digital elevation model data. The positive terrain data acquisition module calculates the positive terrain data of the target area based on the preprocessed digital elevation model data. The module for acquiring digital elevation model data without depressions uses a priority flood algorithm to fill depressions based on the preprocessed digital elevation model data to acquire digital elevation model data without depressions. The water flow direction data acquisition module calculates the water flow direction data of the target area based on the digital elevation model data without depressions using the D8 method. The flow accumulation acquisition module calculates the flow accumulation in the target area based on the water flow direction data of the target area. The module for obtaining the candidate ridge point set identifies areas with zero cumulative flow and positive terrain as candidate ridge points based on the positive terrain data of the target area and the cumulative flow of the target area, thereby obtaining the candidate ridge point set. The module for obtaining the cluster set of noise-removed points uses the DBSCAN algorithm to perform spatial clustering on the candidate ridge point set, divides the clusters according to spatial distance and removes noise points to obtain the cluster set of noise-removed points. The final ridge point acquisition module performs elevation screening on each cluster in the cluster set after removing noise points, and retains ridge points that are higher than the 50th percentile of the cluster elevation as the final ridge points. The ridgeline acquisition module uses a skeletonization algorithm to perform continuous processing based on the final ridge point to obtain the ridgeline.