Complex waterline identification and extraction method based on point cloud

Abstract

The invention relates to the field of surveying and mapping science and technology, in particular to a complex waterline identification and extraction method based on point cloud. The method comprises the following steps: acquiring airborne laser point cloud and ground laser point cloud, and carrying out coordinate unification and fusion to obtain an original point cloud set; setting a height threshold based on the elevation distribution characteristics of the original point cloud set, and filtering noisy points of a near-shore water area to obtain a preprocessed point cloud set; repeated random sampling is carried out in the preprocessed point cloud set, and a preset threshold value is utilized to fit an optimal plane as a projection reference plane; therefore, by introducing a neighborhood geometric analysis and circle center backstepping judgment mechanism based on the projection point cloud, robust recognition of the waterline points in the complex terrain is realized, and the problem that a traditional method is easy to cause misjudgment in vegetation interference, local mutation and multi-scale slope break areas is solved; and the continuity and reliability of waterline extraction in a complex surface environment are improved.

Description

Technical Field

[0001] This invention relates to the field of surveying and mapping science and technology, and in particular to a method for identifying and extracting complex waterline based on point clouds. Background Technology

[0002] Waterline mapping is a crucial element of topographic maps, reflecting the complex processes of river flow dynamics and sediment transport. It provides a scientific basis for sediment deposition analysis, river evolution analysis, and waterway improvement projects. Currently, waterline mapping mainly relies on manual, walking surveying methods such as total stations and GPS. In complex river environments, the point and line characteristics of the waterline are difficult to grasp, resulting in low work efficiency and low mapping accuracy. Furthermore, sandbars and shoals are difficult to access due to poor transportation and sparse siltation, posing high safety risks and requiring intensive work.

[0003] For many years, scholars both domestically and internationally have primarily focused on interpreting waterline boundaries using optical and microwave remote sensing, extracting water boundaries based on the spectral and textural features of images. Some researchers have utilized SAR and panchromatic remote sensing images from different phases of the same tidal cycle, employing translation-invariant discrete wavelet transform algorithms for image fusion and employing a novel fusion rule to extract tidal zone topography. Others have used threshold segmentation, neural network classification, and object-oriented methods to extract waterlines from multispectral data of artificial coastlines, sandy coastlines, silty coastlines, and mangrove coastlines. Some researchers have used the K-means algorithm to coarsely segment island regions in remote sensing images, determining the initial contours of islands based on a grayscale dilation algorithm, and optimizing waterlines using a level set algorithm. However, waterlines are dynamic and variable, exhibiting different spatial morphological characteristics when observed under different time and scale conditions. Furthermore, remote sensing images are significantly affected by illumination intensity, cloud cover, and vegetation cover, making it difficult to meet the requirements for high-timeliness and high-frequency waterline monitoring. Summary of the Invention

[0004] Therefore, it is necessary to provide a method for identifying and extracting complex water edges based on point clouds to solve at least one of the above-mentioned technical problems.

[0005] To achieve the above objectives, a method for identifying and extracting complex waterline edges based on point clouds is provided, the method comprising the following steps: Step S1: Acquire airborne laser point clouds and ground laser point clouds, and perform coordinate unification and fusion to obtain the original point cloud set; Step S2: Based on the elevation distribution characteristics of the original point cloud, set a height threshold to filter out noise points in the nearshore waters and obtain a preprocessed point cloud; Step S3: Repeatedly sample the preprocessed point cloud and use the preset threshold to fit the best plane as the projection reference plane. Step S4: Project the preprocessed point cloud onto the reference plane to obtain the projected point cloud; calculate the basic geometric point cloud data based on the projected point cloud, and use the basic geometric point cloud data to determine the geometric edge detection radius; Step S5: Construct a geometric edge detection model based on the geometric edge detection radius; use the geometric edge detection model to traverse the projection point cloud and extract the water edge line point set.

[0006] The beneficial effects of this invention are as follows: by analyzing the elevation distribution characteristics of point clouds, noise points in nearshore waters are filtered out; by spatial projection, the point cloud is projected onto a fitting plane, preserving local features of land features while reducing data redundancy and effectively reducing the difficulty of computer processing; adaptive parameters of the detection model are calculated based on multi-scale factors such as the average distance, average curvature, and scale factor of the point cloud, improving adaptability; and a geometric edge detection model is constructed to detect point cloud edge features, enabling automatic extraction of complex water edges. The spatial characteristics of point clouds are analyzed mechanistically, enabling effective direct extraction of water edges from complex terrain, with advantages such as high extraction accuracy and strong timeliness. Therefore, this invention, by introducing a neighborhood geometric analysis and center-inverse determination mechanism based on projected point clouds, achieves robust identification of water edge points in complex terrain, solving the problem of misjudgment easily generated by traditional methods in areas with vegetation interference, local abrupt changes, and multi-scale slope break regions, thus improving the continuity and reliability of water edge extraction in complex surface environments. Attached Figure Description

[0007] Figure 1 This is a flowchart illustrating the steps of a method for identifying and extracting complex waterline based on point clouds. Figure 2 This is a schematic diagram for calculating the coordinates of the circle's center. Figure 3 This is a schematic diagram of the waterline detection. Figure 4 A schematic diagram of the terrain point cloud; Figure 5 This is the result image of the waterline extraction. The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0008] The technical method of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without inventive effort are within the scope of protection of this invention.

[0009] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0010] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0011] To achieve the above objectives, please refer to Figures 1 to 5 A method for identifying and extracting complex waterline based on point clouds, the method comprising the following steps: Step S1: Acquire airborne laser point clouds and ground laser point clouds, and perform coordinate unification and fusion to obtain the original point cloud set; Step S2: Based on the elevation distribution characteristics of the original point cloud, set a height threshold to filter out noise points in the nearshore waters and obtain a preprocessed point cloud; Step S3: Repeatedly sample the preprocessed point cloud and use the preset threshold to fit the best plane as the projection reference plane. Step S4: Project the preprocessed point cloud onto the reference plane to obtain the projected point cloud; calculate the basic geometric point cloud data based on the projected point cloud, and use the basic geometric point cloud data to determine the geometric edge detection radius; Step S5: Construct a geometric edge detection model based on the geometric edge detection radius; use the geometric edge detection model to traverse the projection point cloud and extract the water edge line point set.

[0012] In this embodiment of the invention, reference Figure 1 The diagram shows a flowchart illustrating the steps of a point cloud-based method for identifying and extracting complex water edges according to the present invention. In this example, the point cloud-based method for identifying and extracting complex water edges includes the following steps: Step S1: Acquire airborne laser point clouds and ground laser point clouds, and perform coordinate unification and fusion to obtain the original point cloud set; Of particular importance, step S1 includes: Step S11: Acquire airborne laser point cloud data and record the original coordinate system identifier and attitude parameters of the airborne point cloud to form the initial airborne point cloud data; Step S12: Acquire ground laser point cloud data and record the installation height information and scanning direction markings of the ground scanning equipment to form initial ground point cloud data; Step S13: Based on the acquisition parameters of the initial airborne point cloud data and the initial ground point cloud data, the spatial coordinates of the two types of point clouds are uniformly transformed to generate unified coordinate point cloud data. Step S14: Perform overlapping region matching on the unified coordinate point cloud data based on the airborne scanning range and the ground scanning range, and calculate the fusion transformation parameters based on the spatial offset of the corresponding points in the overlapping region; Step S15: Apply the fusion transformation parameters to the unified coordinate point cloud data, perform spatial fusion processing on the airborne point cloud and the ground point cloud, and generate fused point cloud data; Step S16: Perform point density checks and outlier removal on the fused point cloud data, and write the fused points that meet the density consistency requirements into the original point cloud set.

[0013] Specifically, at the data acquisition level, the original coordinate system identifiers, GNSS / IMU attitude and timestamps, sensor intrinsic and extrinsic parameters and scanning deployment information of airborne and ground point clouds are saved to form initial airborne point cloud data and initial ground point cloud data. Based on the sensor pose and the station reference, the two types of initial data are transformed to a unified reference system to generate unified coordinate point cloud data (the transformation process includes rigid transformation through recorded attitude parameters and interpolation or synchronization alignment of the time-series point cloud according to the timestamp). By identifying the overlapping coverage area of ​​airborne and ground point clouds and extracting feature points in the overlapping area or using global nearest neighbor pairs, iterative registration or local minimum residual optimization methods are used to calculate fusion transformation parameters to reduce the relative deviation between the two source points. The obtained fusion transformation parameters are applied to unified coordinate point cloud data to generate fused point cloud data. After fusion, point density statistics (such as based on voxels or neighborhood counting) and outlier detection (such as statistical outlier removal or removal based on neighborhood distance thresholds) are performed on the point cloud. Points that meet the preset density consistency and spatial continuity rules are selected and written into the original point cloud set.

[0014] Step S2: Based on the elevation distribution characteristics of the original point cloud, set a height threshold to filter out noise points in the nearshore waters and obtain a preprocessed point cloud; Preferably, step S2 includes the following steps: Step S21: Divide the original point cloud into several overlapping grid cells, and statistically analyze the elevation distribution within the overlapping grid cells to generate a grid elevation statistical record; Step S22: Calculate the elevation quantile value of the grid cell based on the grid elevation statistics record, and set the height threshold data according to the elevation quantile value; Step S23: Perform near-shore noise filtering on the original point cloud using a preset height threshold to generate a near-shore noise identifier set; Step S24: Write out the set of points in the original point cloud that are not marked by the nearshore noise identifier set as the preprocessed point cloud set.

[0015] It should be noted that the sliding window method is used, the grid size is m×m, and the overlap ratio is 50%.

[0016] In this embodiment, the technical means at the data level is to perform local statistical analysis on the original point cloud by overlapping grid division, adaptively generate a height threshold based on the elevation quantile characteristics of each grid, and perform near-shore noise judgment and removal based on the threshold on the point cloud, thereby outputting a consistent preprocessed point cloud.

[0017] Specifically, the original point aggregator is used to generate overlapping grid cells covering the entire area at a predetermined base scale to ensure neighborhood continuity and boundary smoothness. Within each overlapping grid cell, elevation samples of the points are collected and descriptive statistics (such as median, upper and lower quartiles, selected quantile values, etc.) are calculated to form a grid elevation statistical record. The statistical record contains both location index and sample size information for subsequent reliability assessment. Based on this grid elevation statistical record, several representative quantile values ​​are extracted and mapped to local height thresholds. The threshold setting takes into account local terrain fluctuations and sample density. By combining quantile values ​​with sample size and grid overlap ratio, the influence of extreme noise on the threshold is limited. The original point cloud is judged point by point based on the height threshold. Points located outside the threshold and conforming to the spatial distribution characteristics of nearshore noise are marked as nearshore noise and a nearshore noise identifier set is formed. At the same time, points within the threshold and meeting the minimum neighborhood density requirement are retained. The set of points not marked by nearshore noise is written out according to the original spatial index as a preprocessed point cloud.

[0018] Step S3: Repeatedly sample the preprocessed point cloud and use the preset threshold to fit the best plane as the projection reference plane. Preferably, step S3 includes the following steps: Step S31: Samples are drawn from the preprocessed point cloud after a specified number of random samplings; a candidate plane parameter set is fitted based on the samples; Step S32: Calculate the fitting residuals of the points within the sample for each candidate plane parameter, and generate fitting residual data; Step S33: Based on the fitted residual data, candidate planes are selected according to a preset deviation threshold and a preset in-point ratio threshold to form a set of candidate planes; Step S34: Select the best plane from the candidate plane set according to the minimum residual criterion, and use it as the projection reference plane.

[0019] Specifically, firstly, several sample sets are extracted from the preprocessed point cloud according to a preset number of random samplings. Each sample set contains a certain number of neighborhood points and forms a sample point cluster. For each sample point cluster, a robust fitting method is used to estimate the plane parameters and record them as candidate plane parameter items. The candidate plane parameter items, together with the sample point indices, form a candidate plane parameter set. For each candidate plane parameter, the corresponding sample points are projected onto the plane and the projection residual sequence is calculated. The residual statistics are summarized to form the fitted residual data. At the same time, the proportion of points in the sample that meet the residual threshold is calculated as the in-point ratio index. The candidate plane parameter set is screened using the preset residual threshold and the in-point ratio threshold as dual screening conditions to produce a candidate plane set that meets the requirements of stability and coverage. Based on the combined criteria of minimizing the overall residual and optimizing the sample coverage, a single optimal plane is selected from the candidate plane set, and its parameters are recorded as the optimal plane parameter record. The optimal plane parameter record is used as the parameter input for the projection reference plane in the subsequent step S4.

[0020] Preferably, step S33 includes the following steps: Step S331: Based on the fitted residual data, the candidate plane parameters are initially screened according to the preset deviation threshold and the preset in-point ratio threshold, and the candidate planes that meet the initial screening conditions are stratified according to the scale interval to obtain a scaled candidate plane set. Step S332: In the set of scaled candidate planes, statistically analyze the corresponding interior point ratio distribution and residual distribution characteristics of the candidate planes in each scale layer to generate scaled fitting statistics. Step S333: Calculate the stability score for each candidate plane based on the scaled fitting statistics, and sort them from high to low stability scores to form a stability sequence; Step S334: Perform local maximum detection on the stability sequence to identify candidate planes for cross-scale stable positions, and write the identified candidate planes into the candidate plane set.

[0021] In this embodiment, a set of candidate planes that have been scale-verified is output by performing threshold screening on the residuals and the proportion of interior points of the candidate planes, performing hierarchical statistical fitting on the features according to the scale, calculating the stability score based on the statistical features, and identifying stable candidate planes through cross-scale local maxima detection.

[0022] Specifically, each candidate plane is initially screened based on its fitted residual sequence and in-point ratio threshold according to pre-set residual deviation threshold and in-point ratio threshold. Candidates with excessively large residuals or excessively low in-point ratios are eliminated. The candidate planes that pass the initial screening are divided into several scale intervals according to the spatial / scale attributes of the samples or fitting windows used, resulting in a scaled candidate plane set. Within each scale layer, the distribution characteristics of the in-point ratio (such as median, quantile, and sample size) and the distribution characteristics of the residuals (such as mean, variance, and extreme values) of the candidate planes that pass the initial screening are statistically analyzed. These descriptors are then summarized into scaled fitting statistics to characterize the stability and consistency of the fitting within that scale layer. Based on the above-mentioned scaled fitting statistics, a stability score is calculated for each candidate plane. The stability score is mainly composed of low residuals, high in-point ratio and scale consistency (which can be obtained through weighted combination descriptors). All candidate planes are arranged in descending order of stability score to form a stability sequence. The local maximum detection strategy using a sliding window is employed to identify planes that exhibit prominent stable peaks in adjacent scales or adjacent candidate planes. These planes, which have been determined by local maximums and meet the minimum significance condition, are marked as cross-scale stable positions and added to the candidate plane set. At the same time, the set is accompanied by corresponding scale labels, stability scores, and fitting statistics for subsequent optimal plane selection.

[0023] Preferably, step S331 includes the following: Based on the fitted residual data, the candidate plane parameters are initially screened according to a preset deviation threshold and a preset in-point ratio threshold to obtain a subset of candidate planes that meet the initial screening requirements. The candidate plane subset that meets the initial screening requirements is divided into scale intervals, and the neighborhood coverage is calculated based on the ratio of inliers of the candidate planes within each scale interval. The coverage is then refined into a coverage rate description. A composite stability metric is constructed based on the coverage description and the fitted residual data, and a stability score sequence is generated for each candidate plane. The stability score sequence is subjected to local maximum determination according to the sliding window size to extract candidate planes for cross-scale structural stability, and the extracted candidate planes are reorganized into a scaled candidate plane set according to the scale interval.

[0024] Specifically, the original candidate plane parameters are initially screened item by item using the fitted residual data and the preset in-point ratio threshold. The candidate indices that pass and fail, as well as their residual distribution and in-point ratio, are recorded to form an initial screening result set. The candidate planes that pass the initial screening are divided into several scale intervals according to their corresponding sample window size or fitted scale attribute. In each scale interval, the proportion of in-points along the sample neighborhood of each candidate plane is statistically analyzed and the neighborhood coverage (i.e., the proportion of fitted points in the neighborhood to available neighborhood points) is further calculated. The coverage rate description is generated based on the distribution of coverage in the spatial and scale dimensions. This description includes the median of the coverage rate, the distribution range, and the local hole index to characterize the coverage continuity. The coverage description is combined with the fitting residual statistics of the candidate planes, and a composite stability measure is constructed through weighting or normalization. The corresponding stability score is recorded for each candidate plane in the scale-varying sequence, forming a stability score sequence. The sliding window local maximum detection is applied to the stability score sequence to identify the positions that show significant peaks in adjacent scales or adjacent candidates. The candidate planes corresponding to the identified local peaks are determined to be cross-scale structurally stable units. Finally, these candidate planes that have been determined by local maximums and meet the conditions of minimum coverage and residual significance are reorganized according to their scale labels to generate a scaled candidate plane set containing scale information, stability scores and coverage descriptions.

[0025] Preferably, the formula for extracting the projection reference plane in step S3 includes: Samples are drawn from the preprocessed point cloud at a specified number of random sampling times, typically 3 points as the minimum sample. A fitting plane is then constructed, and the normal vector of the fitting plane is calculated as follows. ; ; In the formula, i, j, k respectively along x, y, z The unit vector in the positive direction of the axis; x, y, z The algebra of the equation; ， ， ), ( ， ， ), ( ， ， These are the spatial coordinates of the three sampled points; 、 、 、 These are the normal vector of the fitted plane and the coefficients of the three equations, respectively. Then, the distance from the point cloud to the fitted plane is calculated by iterating through the points. The formula for calculating the distance from a point to the plane is as follows: ; In the formula, The distance from the point to the plane; ， ， () represents the current spatial coordinates. if ≤2 Then, the point cloud is marked as an interior point of the fitted plane, where, Let the standard deviation be 0.1; Repeat the above steps, iterate and sample 1000 times, and select the fitting plane with the most marked points as the projection reference plane.

[0026] Specifically, three-point samples are randomly extracted multiple times from the preprocessed point cloud according to a specified number of samplings. The geometric relationship of the three points is used to construct candidate plane parameters (represented by the coefficient set of the plane expression). For each candidate plane, the corresponding sample points and surrounding test points are projected onto the plane and the projection residual is calculated to obtain the distance sequence from each point to the plane. The distance sequence is summarized into fitted residual data and the proportion of points that meet the residual threshold is calculated as the in-point ratio index. The candidate planes are initially screened using the preset deviation threshold and the in-point ratio threshold, and the candidates that pass the initial screening are retained. Within the set of candidate planes that pass the initial screening, they are sorted or scored according to the residual statistical characteristics and coverage ratio. Finally, the single best plane is selected as the projection reference plane by using the minimum residual or composite optimal criteria.

[0027] In one implementation of this embodiment, laser point cloud scanning is performed on a section of the Yangtze River Three Gorges Reservoir area, and three spatial points are randomly sampled: Point 1: (10.2, 15.3, 8.5); Point 2: (12.5, 16.8, 8.7); Point 3: (11.0, 17.2, 8.6); Expand the determinant to calculate the normal vector =i(1.5×0.1−0.2×1.9)−j(2.3×0.1−0.2×0.8)+k(2.3×1.9−1.5×0.8); therefore, =−0.23, =−0.07, =8.5 (obtained by substituting the coordinates of point 1).

[0028] For example, to determine whether a point P(11.5, 16.0, 9.2) is an interior point; Substitute to obtain It is 4.34 meters Because of deviation 0..1, =4.34 is greater than 0.2, therefore this point is not an interior point.

[0029] Step S4: Project the preprocessed point cloud onto the reference plane to obtain the projected point cloud; calculate the basic geometric point cloud data based on the projected point cloud, and use the basic geometric point cloud data to determine the geometric edge detection radius; Preferably, step S4 includes the following steps: Step S41: Project the preprocessed point cloud onto the reference plane according to the plane parameters using the projection reference plane to generate the projection point cloud; Step S42: Calculate the point cloud density in each local projection region based on the projected point cloud set to generate basic geometric point cloud data; Step S43: Calculate several candidate detection radii in each projected local area based on the basic geometric point cloud data, and select and determine the detection radius based on the local stability of the candidate detection radii.

[0030] In this embodiment, the preprocessed point cloud is geometrically projected onto a selected reference plane to generate a projected point cloud in planar coordinates. Then, the point cloud density and local geometric features are statistically analyzed in the projection domain to construct basic geometric point cloud data. Based on these statistical features, several candidate detection radii are generated, and the final detection radius is determined by a local stability metric.

[0031] Specifically, rigid projection is performed on each 3D point in the preprocessed point cloud based on the optimal plane parameters to obtain the projected point cloud in the reference plane coordinate system. During the projection process, the original elevation deviation and projection index of each point are retained for subsequent backtracking. It is important to note that the projected point cloud is divided into several overlapping or non-overlapping local projection regions. For each local region, the point cloud density (e.g., based on fixed radius neighborhood counting or voxel lattice counting), the directional index of local point distribution, and the distribution characteristics of the projection height difference are statistically analyzed. The combination of these statistical results constitutes the basic geometric point cloud data of the region. Next, candidate detection radii are generated for each local region using multiple scale candidate values ​​(the scales range from small to large, covering the transition from local microstructures to larger terrain). For each candidate detection radius, corresponding local stability indices are calculated—such as the variance of the number of neighborhood points at that radius, the scale consistency of edge response, and the spatial continuity measure of the detection results. Finally, based on these stability indices, the candidate detection radius that exhibits low variance, high response consistency, and maintains spatial continuity at the scale is selected as the final detection radius for that local region.

[0032] It is important to note that basic geometric point cloud data includes point density, local curvature, normal vector, and elevation change rate, etc.

[0033] Of particular importance, step S4, determining the detection radius, includes: Index projection point cloud ( , ), traverse and calculate to each point in the neighborhood. The average distance, and calculate the average distance of the entire point cloud. ,as follows, ; In the formula, The number of point clouds in the neighborhood. For the point cloud in the neighborhood, Find the centroid of the neighborhood by traversing through it. And based on neighboring points Calculate the covariance matrix ,as follows, ; Decomposition of covariance matrix We obtain the three eigenvalues ​​as follows. ; In the formula, , , , , , , , , For each element in the covariance matrix, find its eigenvalues. , , , and 0 ; Calculate the average curvature of the overall point cloud based on the eigenvalues. as follows, ； In the formula, The overall average curvature of the point cloud; n Projection point cluster ; Based on average distance Mean curvature Calculate the detection radius ,as follows, ; In the formula, For the detection radius, This is a scaling factor (usually set to 1~3). This is the distance weight (usually set to 0.75). This is the curvature weight (usually set to 0.25).

[0034] In this embodiment, the average distance from the point to the center point within the neighborhood is calculated to quantify the scale. The covariance matrix of the neighborhood points is constructed and the eigenvalues ​​are obtained to characterize the local curvature / shape distribution. Then, the average curvature index derived from the eigenvalues ​​is used to adjust the candidate detection radius based on the average distance, thereby obtaining a detection radius that reflects both the neighborhood scale and the local geometric curvature.

[0035] Specifically, within each detection neighborhood, the average distance to the central projection point is calculated by traversing the neighborhood points to quantify the typical radial scale; subsequently, the sample average of the point coordinates is extracted for the same neighborhood, and the covariance matrix (the matrix elements are the mean of the products of the coordinate deviations) is calculated based on the sample bias. This process is then used to solve for the three eigenvalues. The average curvature index is constructed using the first two smallest feature values ​​of each neighborhood under normalization and semi-factoring (this index numerically reflects the flatness and curvature of the neighborhood); finally, the average distance and average curvature are coupled with empirical scaling factors and weights to calculate the detection radius.

[0036] In one implementation of this embodiment, it is assumed that five projection points are selected from a certain region in the projection point cloud, and the three neighboring points of each point are as follows: Center point 1: (5.0, 10.0), Neighboring points: (5.2, 10.1), (4.8, 10.2), (5.1, 9.9) This demonstration only focuses on center point 1. The distances from neighboring points to the center point will then be calculated separately. It is 0.224. It is 0.283. It is 0.141; Then the average neighborhood distance of center point 1 is =0.216; Assume the average neighborhood distances of the other four points are 0.198, 0.235, 0.209, and 0.227, respectively. Overall average distance: =0.217.

[0037] Suppose a neighborhood contains 4 points (k=4): P1:(5.0,10.0,8.5); P2:(5.2,10.1,8.6); P3:(4.9,10.2,8.4); P4:(5.1,9.9,8.5); First, calculate the centroid: =5.05, =10.05, =8.5; The deviations from P1 to P4 are (-0.05, -0.05, 0.0), (0.15, 0.05, 0.1), (-0.15, 0.15, -0.1), and (0.05, -0.15, 0.0). The eigenvalues ​​of the covariance matrix elements are calculated as follows: , , The values ​​are 0.003, 0.008, and 0.019, respectively. Then the curvature of the three neighborhoods 、 、 The values ​​are 0.183, 0.213, and 0.174, respectively. Then the overall average curvature 0.190.

[0038] Using the results of the previous calculations: 0.310 meters.

[0039] Step S5: Construct a geometric edge detection model based on the geometric edge detection radius; use the geometric edge detection model to traverse the projection point cloud and extract the water edge line point set.

[0040] Preferably, step S5 includes the following steps: Step S51: Define the search neighborhood based on the detection radius and calculate the model adaptive parameters to construct a geometric edge detection model based on the detection radius; Step S52: The geometric edge detection model calculates the geometric offset data of the neighborhood points based on the spatial distribution of each projection point relative to the detection radius within the neighborhood; and infers the coordinates of the center point within the neighborhood based on the geometric offset data. Step S53: Use the geometric edge detection model to traverse the cloud of projection points, perform edge determination on each projection point based on the offset of the center point coordinates and the geometric discontinuity of the neighboring points, and write the projection points that meet the edge determination criteria into the waterline point set.

[0041] In this embodiment, the detection radius determined a priori is used as the basis of the neighborhood. The adaptive parameters and center coordinates of the neighborhood are obtained through local statistics and geometric calculation. Then, the edge attributes of the projection points are determined point by point based on the numerical index of the center offset and the geometric discontinuity of the neighborhood, and the water edge line point set is output.

[0042] Specifically, for each projection point in the projection point cloud, a closed neighborhood is constructed with its corresponding detection radius, and the set of projection points within the neighborhood is extracted. The planar coordinates, radial distance from the center point, and relative orientation information of each neighborhood point are recorded. Then, adaptive parameters of the model are calculated on the neighborhood data—including radial distance distribution description (mean, variance, quantile), directional distribution description (azimuth histogram feature), and local density statistics of the neighborhood points—to characterize the geometric structure of the neighborhood and serve as a normalization factor for subsequent judgments. Meanwhile, it is important to note that the radial and directional information of neighboring points is used to calculate the set of possible center point coordinates within the neighborhood through aggregation and weighting methods (i.e., calculating the consistency of distances between several candidate center locations and neighboring points and recording the aggregation strength), and the aggregation degree and stability of this set of center coordinates are evaluated to obtain a single or finite number of center estimates and their aggregation degree indices; after obtaining the center estimates and adaptive parameters, the center offset (the deviation of the projected point from the estimated center relative to the detection radius or the normalized value of the neighborhood statistical scale) is calculated for the projected point, and the geometric discontinuity index of the neighborhood is calculated at the same time (e.g., the local slope change of radial distance or the abrupt change amplitude of azimuth distribution). According to the established judgment rules, the projection points are edge-determined using a joint threshold standard of center offset and geometric discontinuity. Projection points that meet the edge-determining conditions are marked and written into the waterline point set.

[0043] Preferably, step S53 includes the following: For each projection point in the projection point cloud, extract neighboring points within the detection radius centered on the projection point and construct a neighboring point set; calculate the direction vector and distance offset information of each neighboring point relative to the projection point for the neighboring point set; Based on the direction vector and distance offset information, the coordinates of the possible center points of the circle corresponding to the neighboring points are inferred to form a set of center point coordinates; the aggregation degree of the center point coordinate set is evaluated to identify whether there are stable center points in the neighborhood. When the set of center point coordinates does not show a stable aggregation and the directional distribution among neighboring points is significantly discontinuous, the projected point is written into the waterline point set; conversely, when the center point coordinates are highly concentrated and the neighboring points are evenly distributed, the point is excluded as a non-edge point.

[0044] In this embodiment, the direction vector and radial offset of the neighboring points of each projection point are calculated, and the candidate circle center is deduced from these geometric quantities. The aggregation of these candidate circle centers and the continuity of the neighborhood direction distribution are evaluated. Thus, the edge points are identified and the water edge line point set is output based on the criteria of low circle center aggregation degree and discontinuous direction distribution.

[0045] Specifically, for each projection point in the projection point cloud, neighborhood points are first extracted within a predetermined detection radius and a neighborhood point set is formed. The direction vector (the direction information obtained by subtracting the coordinates of the center point from the coordinates of the neighborhood point) and the corresponding radial distance offset (the magnitude of the vector or the radial projection length) are calculated point by point in the neighborhood point set relative to the center point. The direction and distance offset are used as neighborhood geometric description vectors. Based on these descriptive vectors, a set of candidate center point coordinates is generated by fitting several pairs or small samples of neighboring points (e.g., circle fitting of two or three points or minimum variance circle center estimation), and the center point coordinates are summarized to obtain a set of center point coordinates; the aggregation index (including spatial variance of the center position, cluster density and weight ratio of the main cluster points, etc.) is calculated for this set, the directional distribution description of the directional vector set of neighboring points is calculated (continuity is quantified by angle histogram or directional vector variance), and the radial consistency is examined in combination with the radial offset distribution of neighboring points; The following judgments are also required: if the set of center point coordinates shows low aggregation (sparse clustering, large variance, or low proportion of main clusters) and the distribution of neighborhood directions shows obvious breaks or high variance (indicating abrupt changes in the radial direction within the neighborhood), then the center projection point is determined to have edge characteristics and is written into the waterline point set; conversely, if the center coordinates are highly concentrated and the distribution of neighborhood directions is approximately uniform, then it is determined to be a non-edge point and is excluded.

[0046] Preferably, step S5, the extraction of the waterline point set, includes: The geometric edge detection model is used to traverse the projection point cloud to obtain the neighborhood point set, where the neighborhood point set is used... express; ; Among them, based on the detection radius α, the index distance from the projected point cloud is less than 2. The set of all neighborhood points, where the neighborhood point set is used express; The coordinates of the nth neighboring point set; Select the first point cloud in the neighborhood Based on the point to be detected and the neighboring point cloud and detection radius The coordinates of the two center points O1 and O2 are calculated using the following formula. ; ; ; ; ; In the formula, ( , ) Let O1 be the coordinates of the neighborhood point cloud. ( , ) O2 ( , ) The coordinates of the center point; The distance from the neighborhood point cloud coordinates to the circle; Traversal calculation The distances from the remaining point cloud points to O1 and O2; if the distance from all points to O1 or O2 is greater than the radius. If the value is not found, it is marked as a waterline point, and the judgment of that point is terminated. if The distances from the remaining point clouds to O1 and O2 are not all greater than the radius. If a point is not on the waterline, it is marked as a non-waterline point; the remaining point clouds are then used as points to be detected, the center coordinates of the circles are iteratively calculated, and the points to be detected are re-evaluated until... After all tests are completed, the waterline point set is obtained.

[0047] For clarity, the process of calculating the center coordinates of a circle based on neighborhood geometric relationships in this invention can be found in [reference needed]. Figure 2 The diagram shown illustrates the calculation of the center coordinates of the circle.

[0048] It shows the positional relationship between the projection point and the neighboring points in the planar coordinate system, as well as the geometric structure when constructing candidate circles based on the detection radius. Figure 2 The neighborhood reference point is given in the middle. Distance between the current projection point The diagram schematically illustrates two candidate center points, O1 and O2, determined by these two points.

[0049] In this embodiment, a set of neighborhood points is constructed for each projection point, a candidate circle is determined by the first neighborhood point and the point to be tested (with the given detection radius as the circle radius), and the coordinates of two possible circle centers are calculated. Then, the distance distribution from other points in the neighborhood to the circle centers is examined based on these circle centers. The criteria of circle center aggregation / point inclusion relationship and direction distribution determine whether the projection point belongs to the water edge line, thereby extracting the edge line point set in a point-by-point geometric examination manner.

[0050] Specifically, firstly, taking the current projection point as the center and according to the detection radius... Index its neighborhood point set on the projection plane. ( Only those items with a distance of less than 2 from the center point are included. (The points are selected to ensure that a circle with radius α can be formed between the two points), and the planar coordinates of each neighboring point are preserved; then a reference point is taken from the neighborhood. (coordinate ( , ) ) and the point to be measured ( , ) Pair up and calculate the two points. ,like ≤ 2 Then there exists with Let the circle have radius and pass through the two points. Using the midpoint of the two points as a reference, and along the perpendicular direction of the midpoint, we can identify two candidate positions for the circle's center. The center coordinates can be expressed as "midpoint coordinates ± perpendicular offset," where the offset is... (This quantity is a geometric quantity that determines the center of a circle at two points with a given radius); based on this, two sets of center coordinates O1, O2 are obtained; subsequently, the neighborhood point set... For each remaining point in the neighborhood, calculate its Euclidean distance to O1 and O2, and use the comparison between distance and radius as the criterion: if the distance from all points in the neighborhood to O1 (or O2) is greater than 1, then the distance from O1 (or O2) to O2 is greater than 1. If a point is not contained within the circle (i.e., the neighboring points are not included in the circle) or exhibits no stable aggregation at either of the two circle centers, then the point to be measured is considered to be located in the geometrically fractured zone at the boundary and is recorded as a waterline point; conversely, if there exists a circle center such that most of the neighboring points are concentrated within the radius of that circle center (the circle centers aggregate and the neighborhood directional distribution is approximately continuous), then the point to be measured is excluded as a non-edge point; if the determination is unclear, then the next point in the neighborhood is used as a reference. Repeat the above calculations until all resources are used up. .

[0051] For example: Let the detection radius be... =1.0, the point to be measured ( , ) =(0.0,0.0), reference neighborhood point If the distance between the two points is (1.0, 0.0), then the distance between the two points is... =1; the midpoint between the two points is (0.5, 0.0), and the vertical offset is... ≈0.866025; Therefore, the two centers are located at (0.5, +0.866025) and (0.5, −0.866025) respectively.

[0052] Calculate the distance from another point in the neighborhood, for example (0.0, 1.5), to O1: ==−0.5, =0.633975, Distance ≈ 0.807 (less than) This indicates that the neighboring point falls within the radius centered at O1, thus the center of the circle shows convergence, and the point to be tested should be judged as "non-edge"; if the distance from the vast majority of points in the neighborhood to both centers is greater than 1, then the neighboring point is considered to be "non-edge". If the point to be measured is determined to be an edge point, it will be written into the waterline point set.

[0053] To more clearly demonstrate the edge point determination logic of this invention, please refer to the appendix. Figure 3 The diagram shown illustrates the detection of waterline points.

[0054] In the diagram, red dots represent points along the waterline, while blue dots represent points not along the waterline.

[0055] Figure 3 The diagram depicts the process of calculating the set of circle centers based on the neighborhood region formed by the detection radius, and the schematic structure for determining waterline points using geometric discontinuities. The diagram illustrates the "non-edge case with highly clustered circle centers" and the "edge case with discrete circle centers and abrupt changes in direction distribution" through different types of neighborhood distributions.

[0056] Finally, to demonstrate the experimental point cloud features processed by this invention, see [link to relevant documentation]. Figure 4 The topographic point cloud diagram shown is as follows: Figure 5 The image shown is a result of the waterline extraction.

[0057] Figure 4The distribution of typical point cloud data in the reservoir area of ​​the Three Gorges Project on the Yangtze River is shown. (a) Data I is elongated, with the terrain gradually rising towards the center, a smooth outline, and only a small amount of weeds in the area, resulting in a relatively simple overall structure. In contrast, (b) Data II is an elliptical area with a generally flat terrain, a relatively rough edge outline, dense vegetation, and a small ditch structure in the lower left area. The two sets of point clouds show significant differences in morphology and complexity, which can be used to verify the stable identification capability of the present invention for waterline in various complex environments.

[0058] Appendix Figure 5 The results of waterline extraction for the aforementioned point cloud region are presented. As can be seen in the figure, this invention can directly extract continuous waterline structures from point clouds with significant topographic relief, dense vegetation, or small gullies, even when faced with various complex types of waterlines (such as straight lines, arcs, sharp angles, and gullies).

[0059] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the application are intended to be included within the invention.

[0060] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A method for identifying and extracting complex waterline edges based on point clouds, characterized in that, Includes the following steps: Step S1: Acquire airborne laser point clouds and ground laser point clouds, and perform coordinate unification and fusion to obtain the original point cloud set; Step S2: Based on the elevation distribution characteristics of the original point cloud, set a height threshold to filter out noise points in the nearshore waters and obtain a preprocessed point cloud; Step S3: Repeatedly sample the preprocessed point cloud and use the preset threshold to fit the best plane as the projection reference plane. Step S4: Project the preprocessed point cloud onto the reference plane to obtain the projected point cloud; calculate the basic geometric point cloud data based on the projected point cloud, and use the basic geometric point cloud data to determine the geometric edge detection radius; Step S5: Construct a geometric edge detection model based on the geometric edge detection radius; use the geometric edge detection model to traverse the projection point cloud and extract the water edge line point set.

2. The method for identifying and extracting complex waterline based on point cloud as described in claim 1, characterized in that, Step S2 includes the following steps: Step S21: Divide the original point cloud into several overlapping grid cells, and statistically analyze the elevation distribution within the overlapping grid cells to generate a grid elevation statistical record; Step S22: Calculate the elevation quantile value of the grid cell based on the grid elevation statistics record, and set the height threshold data according to the elevation quantile value; Step S23: Perform near-shore noise filtering on the original point cloud using a preset height threshold to generate a near-shore noise identifier set; Step S24: Write out the set of points in the original point cloud that are not marked by the nearshore noise identifier set as the preprocessed point cloud set.

3. The method for identifying and extracting complex waterline based on point clouds according to claim 1, characterized in that, Step S3 includes the following steps: Step S31: Samples are drawn from the preprocessed point cloud after a specified number of random samplings; a candidate plane parameter set is fitted based on the samples; Step S32: Calculate the fitting residuals of the points within the sample for each candidate plane parameter, and generate fitting residual data; Step S33: Based on the fitted residual data, candidate planes are selected according to a preset deviation threshold and a preset in-point ratio threshold to form a set of candidate planes; Step S34: Select the best plane from the candidate plane set according to the minimum residual criterion, and use it as the projection reference plane.

4. The method for identifying and extracting complex waterline based on point clouds according to claim 3, characterized in that, Step S33 includes the following steps: Step S331: Based on the fitted residual data, the candidate plane parameters are initially screened according to the preset deviation threshold and the preset in-point ratio threshold, and the candidate planes that meet the initial screening conditions are stratified according to the scale interval to obtain a scaled candidate plane set. Step S332: In the set of scaled candidate planes, statistically analyze the corresponding interior point ratio distribution and residual distribution characteristics of the candidate planes in each scale layer to generate scaled fitting statistics. Step S333: Calculate the stability score for each candidate plane based on the scaled fitting statistics, and sort them from high to low stability scores to form a stability sequence; Step S334: Perform local maximum detection on the stability sequence to identify candidate planes for cross-scale stable positions, and write the identified candidate planes into the candidate plane set.

5. The method for identifying and extracting complex waterline based on point clouds according to claim 4, characterized in that, Step S331 includes the following: Based on the fitted residual data, the candidate plane parameters are initially screened according to a preset deviation threshold and a preset in-point ratio threshold to obtain a subset of candidate planes that meet the initial screening requirements. The candidate plane subset that meets the initial screening requirements is divided into scale intervals, and the neighborhood coverage is calculated based on the ratio of inliers of the candidate planes within each scale interval. The coverage is then refined into a coverage rate description. A composite stability metric is constructed based on the coverage description and the fitted residual data, and a stability score sequence is generated for each candidate plane. The stability score sequence is subjected to local maximum determination according to the sliding window size to extract candidate planes for cross-scale structural stability, and the extracted candidate planes are reorganized into a scaled candidate plane set according to the scale interval.

6. The method for identifying and extracting complex waterline based on point clouds according to claim 3, characterized in that, The formula for extracting the projection reference plane in step S3 includes: Samples are drawn from the preprocessed point cloud at a specified number of random sampling times, typically 3 points as the minimum sample. A fitting plane is then constructed, and the normal vector of the fitting plane is calculated as follows. ； ； In the formula, i, j, k respectively along x, y, z The unit vector in the positive direction of the axis; x, y, z The algebra of the equation; ， ， ), ( ， ， ), ( ， ， These are the spatial coordinates of the three sampled points; 、 、 、 These are the normal vector of the fitted plane and the coefficients of the three equations, respectively. Then, the distance from the point cloud to the fitted plane is calculated by iterating through the points. The formula for calculating the distance from a point to the plane is as follows: ； In the formula, The distance from the point to the plane; ， ， () represents the current spatial coordinates. if ≤2 Then, the point cloud is marked as an interior point of the fitted plane, where, The standard deviation is set to a preset threshold of 0.

1. Repeat the above steps, iterate and sample 1000 times, and select the fitting plane with the most marked points as the projection reference plane.

7. The method for identifying and extracting complex waterline based on point cloud as described in claim 1, characterized in that, Step S4 includes the following steps: Step S41: Project the preprocessed point cloud onto the reference plane according to the plane parameters using the projection reference plane to generate the projection point cloud; Step S42: Calculate the point cloud density in each local projection region based on the projected point cloud set to generate basic geometric point cloud data; Step S43: Calculate several candidate detection radii in each projected local area based on the basic geometric point cloud data, and select and determine the detection radius based on the local stability of the candidate detection radii.

8. The method for identifying and extracting complex waterline based on point cloud as described in claim 1, characterized in that, Step S5 includes the following steps: Step S51: Define the search neighborhood based on the detection radius and calculate the model adaptive parameters to construct a geometric edge detection model based on the detection radius; Step S52: The geometric edge detection model calculates the geometric offset data of the neighborhood points based on the spatial distribution of each projection point relative to the detection radius; and infers the coordinates of the center point in the neighborhood based on the geometric offset data. Step S53: Use the geometric edge detection model to traverse the cloud of projection points, perform edge determination on each projection point based on the offset of the center point coordinates and the geometric discontinuity of the neighboring points, and write the projection points that meet the edge determination criteria into the waterline point set.

9. The method for identifying and extracting complex waterline based on point clouds according to claim 8, characterized in that, Step S53 includes the following: For each projection point in the projection point cloud, extract neighboring points within the detection radius centered on that projection point and construct a neighboring point set; Calculate the direction vector and distance offset information of each neighboring point relative to the projection point for the neighboring point set; Based on the direction vector and distance offset information, the coordinates of the possible center points of the circle corresponding to the neighboring points are inferred to form a set of center point coordinates; the aggregation degree of the center point coordinate set is evaluated to identify whether there are stable center points in the neighborhood; When the set of center point coordinates does not show a stable aggregation and the directional distribution among neighboring points is significantly discontinuous, the projected point is written into the waterline point set; conversely, when the center point coordinates are highly concentrated and the neighboring points are evenly distributed, the point is excluded as a non-edge point.

10. The method for identifying and extracting complex waterline based on point cloud as described in claim 8, characterized in that, Step S5, the extraction of the waterline point set, includes: The geometric edge detection model is used to traverse the projection point cloud to obtain the neighborhood point set, where the neighborhood point set is used... express; ； Among them, based on the detection radius Indexing distance less than 2 from the projected point cloud The set of all neighborhood points, where the neighborhood point set is used express; The coordinates of the nth neighboring point set; Select the first point cloud in the neighborhood Based on the point to be detected and the neighboring point cloud and detection radius The coordinates of the two center points O1 and O2 are calculated using the following formula. ； ； ； ； ； In the formula, ( , ) represents the coordinates of the neighborhood point cloud, O1( , O2 , () represents the coordinates of the center point; The distance from the neighborhood point cloud coordinates to the circle; Traversal calculation The distances of the remaining point clouds to O1 and O2; if the distance of all points to O1 or O2 is greater than the radius α, then mark them as waterline points and terminate the judgment of that point; if The distances from the remaining point clouds to O1 and O2 are not all greater than the radius. If a point is not on the waterline, it is marked as a non-waterline point; the remaining point clouds are then used as points to be detected, the center coordinates of the circles are iteratively calculated, and the points to be detected are re-evaluated until... After all tests are completed, the waterline point set is obtained.

Images