Power line residue automatic filtering method and system for airborne point cloud DSM production

Abstract

The application relates to a power line residue automatic filtering method and system for an airborne point cloud DSM production, which comprises the following steps: according to a DSM expression category rule, extracting points participating in DSM expression from point cloud data and analyzing local neighborhood features to generate a DSM expression point set with geometric features; screening candidate points based on a geometric feature threshold value and performing spatial clustering to obtain candidate point clusters; constructing a power line corridor by using a tower position to perform spatial constraint and combining an arc sag curve fitting to screen out point clusters conforming to a power line shape; finally confirming power line residue points through DSM gridding and multi-dimensional anomaly verification and removing the power line residue points from the DSM expression point set; performing surface repair on a hollow region formed after the removal, and iteratively optimizing and screening parameters through quality evaluation feedback, and outputting final DSM results after power line residue filtering.

Description

Technical Field

[0001] This invention relates to the field of photogrammetry and remote sensing technology, specifically to an automatic filtering method and system for power line residues in airborne point cloud (DSM) production. Background Technology

[0002] Airborne LiDAR technology can rapidly acquire high-density, high-precision 3D point cloud data of the Earth's surface, and has become an important data source for the production of Digital Surface Models (DSMs). As the fundamental spatial data representing the elevation of the Earth's surface and the tops of various features, DSMs are widely used in aerial photogrammetry, 3D urban modeling, line-of-sight analysis, and natural resource monitoring. The production of airborne point cloud DSMs typically involves several steps, including importing the original point cloud, flight strip adjustment, noise removal, point cloud classification, DSM representation category selection, gridding or interpolation to generate the DSM, and quality control. Point cloud classification is a crucial step in distinguishing different targets such as ground points, buildings, vegetation, power lines, towers, and noise into corresponding categories. When generating the DSM, operators select stable feature categories such as buildings, vegetation, bridges, and fixed structures based on the classification results, while excluding suspended linear targets such as power lines.

[0003] However, existing classification methods often misclassify power line points in forest areas, urban areas, mountain corridors, and areas where power lines cross buildings or tree canopies. These misclassifications are due to factors such as slender and sparse distribution, floating distribution, local breaks, uneven point density, and mixing with tree canopies or building edges, leading to misclassification as vegetation, buildings, or other categories involved in the DSM representation. A small number of residual power line points, once incorporated into the DSM gridding process, will form narrow linear elevation anomalies in the final DSM output, affecting the reliability of DSM profile analysis, visibility analysis, 3D modeling surface quality, and overall output quality control. Existing power line extraction methods primarily focus on identifying power line targets from the raw point cloud, while existing DSM production methods mainly generate the DSM directly based on classification results. There is a lack of a mechanism to verify, identify, and remove residual power lines that have already entered the DSM representation point set. Summary of the Invention

[0004] Based on this, the purpose of the present invention is to provide a method and system for automatic power line residue filtering for airborne point cloud DSM production, which can automatically verify and filter power line residues in the point set participating in DSM expression before DSM generation.

[0005] The objective of this invention is achieved through the following solution:

[0006] In a first aspect, the present invention provides an automatic power line residue filtering method for airborne point cloud (DSM) production, comprising the following steps:

[0007] S1: Read the preprocessed and initially classified airborne LiDAR point cloud data and DSM expression category rules, extract the points participating in DSM expression from the airborne LiDAR point cloud data according to the DSM expression category rules, construct the DSM expression point set, perform local neighborhood feature analysis on each point in the DSM expression point set, calculate the geometric feature index of each point, and generate a DSM expression point set with local geometric features.

[0008] S2: Based on the geometric feature indicators of each point in the DSM expression point set, feature screening is performed, and points that simultaneously meet the preset geometric feature conditions are extracted as candidate points. Spatial neighborhood clustering is performed on the candidate points, and the candidate points that are continuously distributed in space are grouped into the same cluster to generate a candidate point cluster set.

[0009] S3: Read the tower location data of the survey area, construct the power line corridor space range between each tower based on the tower location data, project each candidate point cluster in the candidate point cluster set onto the power line corridor space range for classification, filter out the candidate point clusters that fall within the corridor range, and perform spatial morphology fitting on each point cluster in the candidate point clusters that fall within the corridor range along the corridor direction, calculate the morphological feature parameters of each point cluster, confirm the point clusters that conform to the power line spatial morphology based on the morphological feature parameters, and generate confirmed power line residual point clusters;

[0010] S4: Construct a DSM grid based on the DSM expression point set, project the confirmed electric field line residual point clusters onto the DSM grid, perform anomaly verification on the grid cells at the projection position, confirm the points corresponding to the grid cells that pass the anomaly verification as electric field line residual points, remove the electric field line residual points from the DSM expression point set, and generate an updated DSM expression point set and hole region.

[0011] S5: Based on the updated DSM expression point set and the void region, perform surface repair processing on the void region by neighborhood interpolation or surface fitting to generate the repaired DSM raster, and perform quality evaluation of the repaired DSM raster by residual anomaly detection. Update the threshold parameters for candidate point selection according to the quality evaluation results and re-execute S2 to S5 until the quality evaluation results meet the preset quality requirements, and output the final DSM result.

[0012] Secondly, this invention provides an automatic power line residue filtering system for airborne point cloud (DSM) production, which is configured with the following modules:

[0013] The DSM feature point construction module is used to read the preprocessed and initially classified airborne LiDAR point cloud data and DSM expression category rules, extract the points participating in DSM expression from the airborne LiDAR point cloud data according to the DSM expression category rules, construct the DSM expression point set, perform local neighborhood feature analysis on each point in the DSM expression point set, calculate the geometric feature index of each point, and generate a DSM expression point set with local geometric features.

[0014] The candidate point cluster extraction module is used to perform feature screening based on the geometric feature indicators of each point in the DSM expression point set, extract points that simultaneously meet the preset geometric feature conditions as candidate points, and perform spatial neighborhood clustering on the candidate points to group the spatially continuously distributed candidate points into the same cluster, generating a candidate point cluster set.

[0015] The power line cluster identification module is used to read the tower location data of the survey area, construct the power line corridor space range between each tower based on the tower location data, project each candidate point cluster in the candidate point cluster set onto the power line corridor space range for classification, filter out the candidate point clusters that fall within the corridor range, and perform spatial morphology fitting on each point cluster in the candidate point clusters that fall within the corridor range along the corridor direction, calculate the morphological feature parameters of each point cluster, confirm the point clusters that conform to the power line spatial morphology based on the morphological feature parameters, and generate confirmed power line residual point clusters;

[0016] The power line point removal module is used to construct a DSM grid based on the DSM expression point set, project the confirmed power line residual point clusters onto the DSM grid, perform anomaly verification on the grid cells at the projection position, confirm the points corresponding to the grid cells that pass the anomaly verification as power line residual points, remove the power line residual points from the DSM expression point set, and generate an updated DSM expression point set and hole regions.

[0017] The DSM iterative repair module is used to perform surface repair processing on the void regions based on the updated DSM expression point set and void regions by performing neighborhood interpolation or surface fitting, generating a repaired DSM raster, and performing quality evaluation on the repaired DSM raster by residual anomaly detection. Based on the quality evaluation results, the threshold parameters for candidate point selection are updated and S2 to S5 are re-executed until the quality evaluation results meet the preset quality requirements, and the final DSM result is output.

[0018] Thirdly, this application provides a computer device, including a memory and a processor. The memory stores a computer program, and when the processor executes the computer program, it implements any of the above-mentioned automatic power line residue filtering methods for airborne point cloud (DSM) production.

[0019] Fourthly, this application provides a computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements any of the above-mentioned automatic power line residue filtering methods for airborne point cloud (DSM) production.

[0020] In summary, the automatic power line residue filtering method for airborne point cloud (DSM) production provided in this application constructs a complete technical chain from DSM expression point set feature analysis, candidate point screening, corridor constraint and spatial morphology fitting, DSM raster anomaly verification to closed-loop quality evaluation. This method can automatically identify and filter out power line points mixed into the DSM expression point set due to classification residues before DSM generation. This method effectively distinguishes power line remnants from easily confused features such as tree canopies, branches, and building edges by performing local neighborhood feature analysis and multi-level screening on the DSM expression point set, combined with the unique corridor spatial constraints and sag morphological features of power lines for dual confirmation, thereby reducing the risk of misjudgment. By performing multi-dimensional anomaly verification at the DSM grid scale, it ensures that the remnant point identification results are consistent with the actual DSM surface impact, avoiding discrimination bias caused by relying solely on local point cloud features. By performing layered surface repair on the removed areas and combining it with quality evaluation to form a closed-loop feedback mechanism, it can adaptively adjust the screening parameters and iteratively optimize the processing effect, thereby ensuring the elevation continuity and morphological consistency between the repaired area and surrounding features. Ultimately, it achieves the technical effect of eliminating narrow-band linear elevation anomalies in DSM results, improving DSM surface quality, and enhancing the level of production automation.

[0021] To better understand and implement this invention, the following detailed description is provided in conjunction with the accompanying drawings. Attached Figure Description

[0022] Figure 1 A flowchart illustrating an automatic power line residue filtering method for airborne point cloud (DSM) production, provided as an embodiment of this application;

[0023] Figure 2 A constraint diagram of the candidate range for pole corridors and power lines provided in the embodiments of this application;

[0024] Figure 3 The sag model fitting and residual point confirmation diagram provided in the embodiments of this application;

[0025] Figure 4 This application provides a DSM narrowband linear elevation anomaly verification map for embodiments of the invention.

[0026] Figure 5 This is a schematic diagram of a power line residue automatic filtering system for airborne point cloud (DSM) production, provided as another embodiment of this application. Detailed Implementation

[0027] To facilitate understanding of the present invention, a more complete description will be given below with reference to the accompanying drawings. Preferred embodiments of the invention are shown in the drawings. However, the invention can be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided to provide a thorough and complete understanding of the disclosure of the invention.

[0028] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used herein in the description of the invention is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0029] In one embodiment, such as Figure 1 As shown, an automatic power line residue filtering method for airborne point cloud (DSM) production is provided. This embodiment illustrates the method using a terminal as an example. It is understood that this method can also be applied to a server, or to a system including both a terminal and a server, and implemented through interaction between the terminal and the server. In this embodiment, the method includes the following steps:

[0030] S1: Read the preprocessed and initially classified airborne LiDAR point cloud data and DSM expression category rules. Extract the points participating in DSM expression from the airborne LiDAR point cloud data according to the DSM expression category rules, construct the DSM expression point set, perform local neighborhood feature analysis on each point in the DSM expression point set, calculate the geometric feature index of each point, and generate a DSM expression point set with local geometric features.

[0031] Specifically, the system reads airborne LiDAR point cloud data after preprocessing, flight strip adjustment, and initial classification. Each point cloud point contains three-dimensional coordinates, classification code, flight strip number, echo count, scan angle, reflection intensity, and GPS time information. The system also reads the expression category rule set used in DSM production, which defines the categories participating in DSM expression. Based on the expression category rule set, the system extracts all points whose classification codes belong to the expression categories from the original point cloud data, constructing a DSM expression point set. For each point in the DSM expression point set, the system establishes a local neighborhood search structure, searching for the K nearest neighbors of each point with a preset search radius. A covariance matrix is ​​constructed based on the K nearest neighbors of each point. The covariance matrix is ​​then decomposed into eigenvalues ​​to obtain three eigenvalues ​​and corresponding eigenvectors. Geometric characteristic indices such as linearity, flatness, and divergence are calculated based on the eigenvalues. Finally, the system uses points with lower elevations within the neighborhood to fit the local reference surface elevation, calculating the suspension height of each point relative to the local reference surface. The system assigns the calculated geometric feature indices and suspension height to the corresponding points in the DSM expression point set, generating a DSM expression point set with local geometric features. This point set serves as the basic data pool for all subsequent filtering and processing operations.

[0032] S2: Based on the geometric feature indicators of each point in the DSM expression point set, feature screening is performed, and points that simultaneously meet the preset geometric feature conditions are extracted as candidate points. Spatial neighborhood clustering is then performed on the candidate points to group the spatially continuously distributed candidate points into the same cluster, generating a candidate point cluster set.

[0033] Specifically, the system performs point-by-point discrimination of the geometric feature indicators of each point in the DSM representation point set based on preset linearity threshold, hover height threshold, cluster width threshold, and cluster length threshold. Points with linearity greater than the linearity threshold, hover height greater than the hover height threshold, a local cluster width less than the width threshold, and a continuous length greater than the length threshold along the main direction are marked as candidate points, generating a candidate point set. The system uses a density-based spatial clustering algorithm to perform spatial neighborhood clustering on the candidate point set, setting neighborhood radius and minimum number of points parameters. The system traverses each point in the candidate point set, searching for other candidate points within the neighborhood radius centered on that point. Points meeting the minimum number of points are grouped into the same cluster, and this cluster continues to expand outward until no new points are added. Through the above clustering process, the system groups spatially continuous candidate points into the same cluster, generating a candidate point cluster set. Each candidate point cluster contains the three-dimensional coordinates, original classification code, and local geometric feature indicators of its member points.

[0034] S3: Read the tower location data of the survey area, construct the power line corridor space range between each tower based on the tower location data, project each candidate point cluster in the candidate point cluster set onto the power line corridor space range for classification, filter out the candidate point clusters that fall within the corridor range, and perform spatial morphology fitting on each point cluster in the candidate point clusters that fall within the corridor range along the corridor direction, calculate the morphological feature parameters of each point cluster, confirm the point clusters that conform to the power line spatial morphology based on the morphological feature parameters, and generate confirmed power line residual point clusters.

[0035] Specifically, the system reads the tower location data of the survey area, which includes the three-dimensional coordinates of each tower. Based on the tower order in the location data, the system sequentially connects adjacent towers to generate the centerline of each tower corridor segment. Using the centerline as a reference, the system expands the corridor search width to both sides by a preset range, constructing the spatial range of the power line corridor between the towers. The system projects each candidate point cluster from the candidate point cluster set onto the power line corridor spatial range and calculates the proportion of points within the corridor range out of the total number of points in the cluster as the corridor assignment score. The system marks candidate point clusters with a corridor assignment score greater than a preset threshold as candidate point clusters falling within the corridor range.

[0036] For each candidate point cluster falling within the corridor area, the system establishes a local coordinate system based on the centerline direction of the corridor, transforming the candidate point cluster from the global coordinate system to this local coordinate system. In the local coordinate system, the system performs spatial morphology fitting on the candidate point cluster along the corridor direction, automatically selecting a straight line model, parabolic model, or catenary model for fitting based on the span of the point cluster along the corridor direction. The system uses a random sampling consensus algorithm to obtain initial model parameters and employs a nonlinear least squares method to optimize and solve the model parameters. The system calculates the root mean square error of the fitting residuals and the lateral width of the point clusters as morphological feature parameters. The system identifies point clusters with a lateral width less than a width threshold, a fitting residual less than a preset adaptive threshold, and consistent sag opening directions as point clusters conforming to the spatial morphology of power lines, generating confirmed residual power line point clusters.

[0037] S4: Construct a DSM grid based on the DSM expression point set, project the confirmed electric field residual point clusters onto the DSM grid, perform anomaly verification on the grid cells at the projection position based on the preset anomaly verification conditions, confirm the points corresponding to the grid cells that pass the anomaly verification as electric field residual points, remove the electric field residual points from the DSM expression point set, and generate an updated DSM expression point set and hole region.

[0038] Specifically, the system establishes a DSM grid at a preset resolution, rasterizes all points in the DSM representation point set, and takes the maximum elevation value of all points within the grid for each grid cell to generate the DSM grid. The system projects the points in the confirmed power line remnant point cluster onto the DSM grid, determines the grid cell position corresponding to each point, and marks the grid cells containing the confirmed power line remnant points as the initial anomalous grid. The system calculates the projection trajectory of the theoretical sag curve on the DSM grid based on the fitting model parameters of the confirmed power line remnant point cluster, performs a union operation on the grid cells traversed by the theoretical trajectory and the initial anomalous grid, and generates a complete set of anomalous grids. The system performs grid morphology verification on each grid cell in the anomalous grid set, calculates the ratio of the principal axis length along the corridor direction to the vertical width of the anomalous grid set, and determines whether this ratio is greater than a preset morphology threshold.

[0039] The system extracts DSM profile lines at anomaly locations along the vertical corridor direction, calculates the first-order difference of the profile elevation, and determines whether there is a narrow-band elevation abrupt change at the anomaly location, characterized by an initial rise followed by a fall. The system utilizes an established local neighborhood search structure to map point cloud neighborhoods to grid neighborhoods, calculates the elevation standard deviation within the anomaly grid neighborhood, and determines whether the anomaly grid elevation is significantly higher than the neighborhood mean and whether the neighborhood surface elevation change is gradual. The system identifies points corresponding to grid cells that simultaneously pass grid morphology verification, profile elevation anomaly verification, and neighborhood support verification as residual power line points. The system removes these residual power line points from the DSM representation point set, generates an updated DSM representation point set, and marks the spatial region covered by the anomaly grid set as a void region.

[0040] S5: Based on the updated DSM expression point set and the void region, perform surface repair processing on the void region by neighborhood interpolation or surface fitting to generate the repaired DSM raster, and perform quality evaluation of the repaired DSM raster by residual anomaly detection. Update the threshold parameters for candidate point selection according to the quality evaluation results and re-execute S2 to S5 until the quality evaluation results meet the preset quality requirements, and output the final DSM result.

[0041] Specifically, the system uses the updated DSM representation point set as the data base and the void region as the spatial scope for repair operations, performing surface repair processing on each void grid. The system automatically selects a repair strategy based on the size and shape of the void region. When the void region diameter is less than a preset grid multiple, an inverse distance weighted interpolation strategy is used, employing only the elevation values ​​of the effective grids surrounding the void for weighted interpolation. When the void region extends linearly, a profile trend repair strategy is used, extracting the effective DSM elevation values ​​on both sides of the void along the power line corridor direction to form a profile line, and filling the void along the profile direction using linear interpolation or spline interpolation. When the void region exceeds a preset size, a neighborhood surface fitting strategy is used, fitting a quadratic surface equation using the effective grid points surrounding the void, solving for the surface coefficients using the least squares method, and substituting the void grid coordinates into the surface equation to calculate the repair elevation value. The system calculates the repair elevation value grid by grid according to the selected repair strategy, generating a complete repaired DSM grid.

[0042] Preferably, the system performs residual anomaly detection on the repaired DSM raster, using the same method as the aforementioned anomaly verification to count the number of linear anomaly grids still existing in the repaired DSM and calculate the linear anomaly residual rate. The system calculates the gradient difference between the repaired elevation and the elevation of the neighboring effective grids in the repaired area, and calculates the slope difference between the repaired grid and the neighboring grids in the repaired area. The system evaluates the repair quality based on three indicators: linear anomaly residual rate, gradient difference, and slope difference. When the linear anomaly residual rate exceeds a preset residual rate threshold, the system lowers the linearity threshold and the suspension height threshold, and re-executes the candidate point extraction and subsequent processing steps based on the original DSM expression point set. When the gradient difference or slope difference value of the repaired area exceeds a preset range, the system traces the judgment criteria of the removed points according to the removal log, recovers some of the removed points, and re-executes the surface repair and evaluation steps. The system repeats the above iterative process until the quality evaluation result meets the preset quality requirements. When the number of iterations reaches a preset upper limit and the quality requirements are still not met, the system marks the anomaly area as a manually reviewed patch and generates an anomaly cause analysis report. The system outputs the repaired DSM raster, the set of power line residual correction points, the set of points to be manually reviewed, the repair layer of the removed areas, the DSM quality evaluation report, and the processing log as the final results.

[0043] In summary, the automatic power line residue filtering method for airborne point cloud (DSM) production provided in this application constructs a complete technical chain from DSM expression point set feature analysis, candidate point screening, corridor constraint and spatial morphology fitting, DSM raster anomaly verification to closed-loop quality evaluation. This method can automatically identify and filter out power line points mixed into the DSM expression point set due to classification residues before DSM generation. This method effectively distinguishes power line remnants from easily confused features such as tree canopies, branches, and building edges by performing local neighborhood feature analysis and multi-level screening on the DSM expression point set, combined with the unique corridor spatial constraints and sag morphological features of power lines for dual confirmation, thereby reducing the risk of misjudgment. By performing multi-dimensional anomaly verification at the DSM grid scale, it ensures that the remnant point identification results are consistent with the actual DSM surface impact, avoiding discrimination bias caused by relying solely on local point cloud features. By performing layered surface repair on the removed areas and combining it with quality evaluation to form a closed-loop feedback mechanism, it can adaptively adjust the screening parameters and iteratively optimize the processing effect, thereby ensuring the elevation continuity and morphological consistency between the repaired area and surrounding features. Ultimately, it achieves the technical effect of eliminating narrow-band linear elevation anomalies in DSM results, improving DSM surface quality, and enhancing the level of production automation.

[0044] In one embodiment, S1 of the automatic power line residue filtering method for airborne point cloud (DSM) production provided by the present invention specifically includes the following steps:

[0045] S11: For each point cloud point in the preprocessed and initially classified airborne LiDAR point cloud data, a spatial neighborhood retrieval data structure is constructed using the spatial distribution of the three-dimensional coordinates of the airborne LiDAR point cloud data as an index. The spatial neighborhood retrieval data structure is used to find the set of neighboring points of each point cloud point within a preset radius, and the spatial neighborhood index relationship of each point is generated.

[0046] S12: Based on the spatial neighborhood index relationship, extract the three-dimensional coordinate matrix of the neighborhood point set of each point in the DSM expression point set. After centering the three-dimensional coordinate matrix, construct the covariance matrix. Obtain three eigenvalues ​​by performing eigenvalue decomposition on the covariance matrix. Calculate the linearity of each point based on the distribution ratio of the three eigenvalues ​​and generate the linear feature parameters of each point.

[0047] S13: Extract the high-order sequence of the neighborhood point set of each point based on the spatial neighborhood index relationship. After sorting the high-order sequence from low to high, cut off the low end part as the local reference surface candidate point. Perform plane fitting on the local reference surface candidate point to obtain the local reference surface elevation. Subtract the local reference surface elevation from the elevation value of each point to obtain the floating height of each point and generate the floating feature parameters of each point.

[0048] S14: Associate and store the linear and dangling feature parameters of each point into the attribute field of each point to generate a DSM expression point set with local geometric features.

[0049] Specifically, the system uses the 3D spatial coordinates of each point in the preprocessed and initially classified airborne LiDAR point cloud data as an index to construct a spatial neighborhood retrieval data structure. This data structure is organized in a k-dimensional tree format, recursively dividing the point cloud space into hyperrectangular regions, enabling subsequent neighborhood retrieval operations to be performed with logarithmic time complexity. The system traverses each point in the DSM representation point set, sets a spatial retrieval radius centered on that point, and uses the k-dimensional tree structure to find all neighboring points within that radius, recording the index identifiers of these neighboring points to form the spatial neighborhood index relationship for that point. When constructing the index relationship, the system calculates an adaptive retrieval radius based on the local density distribution characteristics of the point cloud; the calculation formula is as follows:

[0050]

[0051] in, The total number of points in the DSM expression point set. For the first The Euclidean distance between a point and its nearest neighbor. The arithmetic mean of the nearest neighbor distances for all points. This is a dimensionless adjustment factor, ranging from 0.5 to 1.5, used to control the sensitivity of the search radius to local density fluctuations. The formula's dimension is meters; the first term, the average nearest neighbor distance, is in meters, and the second term, the standard deviation, is also in meters. The sum of these two terms results in identical dimensions. The system uses this adaptive radius as the initial search radius for each point and fine-tunes it based on the actual number of points retrieved within the neighborhood. If the number of retrieved points is below a set lower limit, the radius is increased; conversely, it is decreased, ensuring that the size of each point's neighborhood set remains within a stable range, ultimately generating a stable spatial neighborhood index relationship for each point.

[0052] Based on spatial neighborhood indexing, the system extracts the 3D coordinate matrix of its neighborhood points for each point in the DSM representation point set. This matrix has three rows (number of neighborhood points) and three columns (corresponding to the x-axis, y-axis, and elevation axis). The system centers this coordinate matrix by calculating the arithmetic mean of each column and subtracting the mean from each element in that column, resulting in a zero-mean coordinate deviation matrix. The system then uses this deviation matrix and its transpose to perform matrix multiplication to construct a 3D covariance matrix. This matrix is ​​a symmetric positive semi-definite matrix, where the diagonal elements represent the variance along each coordinate axis, and the off-diagonal elements represent the covariance between the coordinate axes. The system performs eigenvalue decomposition on this covariance matrix, using the Jacobi iteration method to solve for the eigenvalues ​​and eigenvectors, obtaining three non-negative eigenvalues ​​arranged in descending order. The system calculates the linearity of each point based on the distribution ratio of the three eigenvalues. This linearity index uses a weighted anisotropy ratio, calculated using the following formula:

[0053]

[0054] in, , , The characteristic values ​​are arranged in descending order, and their units are all square meters. It is a dimensionless, extremely small positive number used to prevent the denominator from being zero. This is a dimensionless adjustment parameter, ranging from 0.2 to 0.8. This parameter controls the weighting of the ratio of the second to third eigenvalues ​​on linearity. In the formula, the numerator of the first fraction is the difference between the largest and second largest eigenvalues, with dimensions in square meters; the denominator is the sum of the three eigenvalues, also with dimensions in square meters. Therefore, the fraction as a whole is dimensionless. The ratio in the exponential term... Dimensionless, the output value of the exponential function is dimensionless. This linearity index reflects the combined effect of the degree of aggregation and anisotropy intensity of the point cloud along the principal direction, and expresses the geometric properties of slender linear structures more stably than a single eigenvalue ratio. The system establishes a correlation between the calculated linear feature parameters and each point.

[0055] Based on spatial neighborhood indexing, the system extracts the elevation sequence of its neighborhood points for each point in the DSM expression point set. The system arranges this elevation sequence in ascending order of value, and extracts a significant portion of the sequence from the beginning of the sequence. These points correspond to points with lower elevations in the neighborhood, reflecting the elevation distribution of the local surface or stable features. The system performs plane fitting on the extracted low-end point set, using the least squares method to solve for the intercept term and two slope coefficients of the fitted plane. The elevation value of the fitted plane at the plane coordinates of the point to be analyzed is used as the local reference surface elevation. During the fitting process, the system assigns adaptive weights to each low-end point. The weight is related to the order of the point's elevation in the low-end sequence; the earlier the point is in the sequence, the greater the weight, to enhance the contribution of the most stable surface to the reference surface. The system subtracts the local reference surface elevation value from the elevation value of the point to be analyzed to obtain the suspension height of that point, calculated using the following formula:

[0056]

[0057] in, Let be the elevation value of the point to be analyzed, in meters. The total number of low endpoints captured is dimensionless. The serial numbers of the lower endpoints, sorted by elevation, are dimensionless. For the first The estimated elevation of the reference surface after plane fitting of the lowest endpoints, in meters. The dimensionless attenuation coefficient, ranging from two to five, controls the attenuation rate of the weight as the ranking position decreases. In the formula, the numerator is the weighted sum of the estimated datum elevations of each low endpoint multiplied by the exponential weights, in meters; the denominator is the sum of the weights, also in meters. Dividing the two yields the weighted average datum elevation, in meters. Final suspension height. The unit of measurement is meters. This weighting strategy assigns higher weights to reference surface estimates that are closer to the lowest point on the Earth's surface, thereby reducing the interference of vegetation or building base points on the reference surface. The system establishes a correlation between the calculated suspension characteristic parameters and each point.

[0058] Furthermore, the system writes the linear and levitation feature parameters of each point into the attribute fields of the corresponding points in the DSM expression point set. The system allocates two independent storage areas in the attribute fields to store the linear and levitation feature parameter values ​​respectively, while retaining the original 3D coordinates, classification code, and flight strip number of each point. During the writing process, the system performs numerical verification on the linear and levitation feature parameters of each point, checking whether the parameter values ​​fall within a reasonable range. If the parameter values ​​exceed the preset range, the system marks the point as a feature anomaly and records its index, but does not truncate or correct the parameter values ​​to ensure that subsequent processing can obtain the original calculation results. The system associates and stores the linear and levitation feature parameters of each point with the spatial neighborhood index, forming a complete data record structure. The system aggregates all points with completed feature parameter associations to generate a DSM expression point set with local geometric features. This point set contains the spatial coordinates, original classification attributes, spatial neighborhood index links, linear and levitation feature parameters of each point, serving as input data for subsequent screening of suspected residual power line points. The system writes the point set into the memory cache and generates a metadata description file for the point set, recording the calculation version number of the feature parameters, the neighborhood search parameters, and the feature value decomposition parameters.

[0059] In one embodiment, S2 of the automatic power line residue filtering method for airborne point cloud (DSM) production provided by the present invention specifically includes the following steps:

[0060] S21: Based on the linear feature condition, the linear feature parameters of each point in the DSM expression point set are judged one by one, and the points with linear feature parameters greater than the preset linearity threshold are marked as initial points to generate a set of initial linear feature points.

[0061] Among them, the linear feature condition is used to indicate the degree of concentration of point cloud distribution along a certain direction in a local neighborhood. The larger the linear feature parameter, the more significant the extension trend of the point cloud along the corresponding direction. It is used to distinguish between slender linear targets and area or scattered ground features.

[0062] S22: Based on the suspension feature conditions, the suspension feature parameters of each point in the initial selection set of linear features are judged one by one. Points whose suspension feature parameters are less than or equal to the preset suspension height threshold are screened out from the initial selection set of linear features, and a candidate point set that satisfies the dual condition constraints is generated.

[0063] Among them, the suspension feature condition is used to indicate the degree of elevation of the point cloud relative to the ground surface or stable ground surface below it. The greater the suspension height, the more likely the point cloud is to be an overhead target. It is used to distinguish the power line remnants from points that are close to the ground or attached to the surface of buildings or vegetation.

[0064] S23: Based on the spatial neighborhood index relationship, perform neighborhood connectivity analysis on each point in the candidate point set, classify the candidate points that satisfy the preset spatial adjacency relationship into the same connected component, and calculate the spatial extension span along the main direction and the lateral spread width perpendicular to the main direction for each connected component, and generate the geometric morphological parameters of each connected component.

[0065] S24: Based on geometric morphology parameters, each connected component is screened out, and connected components with a spatial extension span less than a preset length threshold or a lateral spread width greater than a preset width threshold are removed. The candidate points in the remaining connected components are used as suspected residual candidate points of the power line, and a set of candidate point clusters is generated.

[0066] Specifically, the system reads the linear feature parameters of each point in the DSM representation point set and performs numerical discrimination on each parameter. The system compares the linear feature parameter value of each point with a preset linearity threshold. Points with linear feature parameter values ​​greater than the threshold are marked as initial selection points, while points with linear feature parameter values ​​not greater than the threshold are not marked at this stage. After the system completes the discrimination of all points, it extracts all marked initial selection points to form a linear feature initial selection point set. During the discrimination process, the system synchronously records the spatial coordinates and attribute information of each initial selection point to ensure that the point set possesses complete spatial distribution characteristics. The linear feature parameter reflects the degree of concentration of the point cloud distribution along a certain direction within a local neighborhood. The larger the parameter value, the more significant the extension trend of the point cloud along the corresponding direction. This attribute allows slender linear targets to form distinguishable numerical differences from area-distributed or scattered targets. During the discrimination process, the system uses an adaptive weight adjustment mechanism to enhance the linear feature parameters. The calculation formula is as follows:

[0067]

[0068] in, This is the original weighted linearity value obtained from the calculation; this value is dimensionless. The effective neighborhood number is dimensionless and is determined by the ratio of the spatial distribution variance of the points within the neighborhood to the global distribution variance. and The lower and upper limits for the number of neighboring points are both dimensionless. The dimensionless gain coefficient takes values ​​between 0.1 and 0.5. The dimensionless term is used for extremely small positive numbers to prevent the denominator from being zero. The strengthening term within square brackets in the formula is dimensionless, and the whole thing is multiplied by . It remains dimensionless. This enhancement strategy assigns a higher linear feature response to points with a neighborhood number close to the upper limit, so that the linear structure in the dense region obtains a more stable pass rate in subsequent screening.

[0069] The system reads the levitation feature parameters of each point in the initial linear feature selection set and performs numerical discrimination on each parameter. The system compares the levitation feature parameter value of each point with a preset levitation height threshold. Points with levitation feature parameter values ​​greater than the threshold are retained in the point set, while points with levitation feature parameter values ​​not greater than the threshold are removed from the initial linear feature selection set. After completing the discrimination and removal operations for all points, the system aggregates the retained points into a candidate point set that satisfies dual constraints. Each point in this candidate point set simultaneously satisfies both the linear feature condition and the levitation feature condition, meaning that each point exhibits a thin, elongated linear distribution within its local neighborhood and is elevated relative to the underlying reference surface. During the removal process, the system indexes and records the removed points to trace the reasons for their exclusion in subsequent iterations. The levitation feature parameter reflects the degree of elevation of the point cloud relative to its underlying stable surface. A larger parameter value indicates that the point cloud is more likely to be an overhead target. This attribute allows residual power line points to form a distinguishable numerical difference from points that are ground-level or attached to building surfaces or vegetation surfaces. The system performs neighborhood normalization on the suspension height parameter during the discrimination process, and the calculation formula is as follows:

[0070]

[0071] in, The original suspension height value at this point is measured in meters. The number of neighborhood points is dimensionless. For the neighboring region The unit of the suspension height of each point is meters. Let the standard deviation of the suspension height values ​​at each point within the neighborhood be in meters. The dimensionless adjustment coefficient ranges from 0.5 to 2.0. The denominator is the sum of the mean and weighted standard deviation of the neighborhood suspension height, with the dimension in meters. The dimension is meters, and the two are divided by the following. Dimensionless. This normalization strategy allows the determination of suspension height at a local neighborhood scale, avoiding the global threshold failure problem caused by differences in terrain undulations or vegetation height in different regions.

[0072] Preferably, the system performs neighborhood connectivity analysis on each point in the candidate point set based on spatial neighborhood index relationships. Starting with each point in the candidate point set, it searches for neighboring candidate points using the spatial neighborhood index relationship. If the spatial distance between two candidate points is less than a set adjacency threshold, the two points are determined to have a spatial adjacency relationship. The system uses a traversal method to pass adjacency relationships to all candidate points, grouping directly adjacent or indirectly adjacent through intermediate points into the same connected component. Each connected component represents a spatially continuously distributed cluster of candidate points, constituting an independent processing unit. The system calculates the geometric parameters of each connected component, including the spatial extension span along the principal direction and the lateral spread width perpendicular to the principal direction. The system first performs principal component analysis on the spatial coordinates of each point within the connected component, obtaining the principal eigenvector direction as the principal direction. Then, it calculates the difference between the maximum and minimum values ​​of the projected coordinates of each point along this principal direction as the spatial extension span, and calculates the bounding box width of the projected point set of each point on the plane perpendicular to the principal direction as the lateral spread width. The system performs local neighborhood correction on the lateral spread width of each connected component, and the calculation formula is as follows:

[0073]

[0074] in, The number of points within a connected component is dimensionless. For the first The lateral offset of a point in the plane perpendicular to the principal direction, expressed in meters, is the distance from the projection of that point onto the plane perpendicular to the principal direction to the centroid of the projection point. The dimensionless adjustment coefficient ranges from 0.1 to 0.3. The first term in the formula is the square root of the arithmetic mean of the squares of the lateral offsets at each point, with units of meters. The second term is the sum of the squares of the ranges of the lateral offsets. The product of these units, with units in square meters, and the square root of the whole, has units in meters. This modified formula integrates the central tendency and dispersion of the lateral width, making the morphological parameters robust to outliers.

[0075] Furthermore, the system filters and classifies connected components based on their geometric morphological parameters. The system compares the spatial span of each connected component with a preset length threshold and its lateral width with a preset width threshold. For connected components with a spatial span less than the length threshold, the system determines that the component lacks the continuous extension length required for electric field line remnants and removes it from the candidate set. For connected components with a lateral width greater than the width threshold, the system determines that the component lacks the narrow morphological features required for electric field line remnants and removes it from the candidate set. The system retains connected components that simultaneously satisfy both the spatial span greater than the length threshold and the lateral width less than the width threshold, designating all candidate points in these retained components as potential electric field line remnant candidate points. The system aggregates all retained connected components into a candidate point cluster set, with each connected component corresponding to a point cluster in the candidate point cluster set. Each point cluster contains the spatial coordinates, original classification code, linear feature parameters, and levitation feature parameters of its member points. During the filtering process, the system records the removal reason label for each removed connected component, marking it as either insufficient length or excessive width. After filtering all connected components, the system performs integrity checks on the remaining candidate point clusters, checking whether the member points of each point cluster have continuous spatial neighborhood connections. For isolated points caused by filtering, the system removes them from their respective point clusters and recalculates the point cluster boundaries.

[0076] In one embodiment, S3 of the automatic power line residue filtering method for airborne point cloud (DSM) production provided by the present invention specifically includes the following steps:

[0077] S31: Read the tower location data of the survey area, and based on the three-dimensional coordinates of each adjacent tower in the tower location data, take the spatial line between two adjacent towers as the corridor centerline, and expand the preset width range to both sides with the corridor centerline as the center to construct the power line corridor spatial range corresponding to each span, and generate a set of corridor spatial constraint areas.

[0078] S32: Project each member point of the candidate point cluster in the candidate point cluster set to the corridor spatial constraint region set for spatial classification. Calculate the proportion of the number of points of each candidate point cluster falling into any corridor spatial constraint region to the total number of candidate point clusters as the corridor classification score. Based on the corridor classification score, select candidate point clusters with scores greater than the preset corridor classification threshold and generate a candidate point cluster set in the corridor.

[0079] S33: Based on the direction of the corridor's central axis to which each candidate point cluster belongs, transform the member points of each candidate point cluster in the candidate point cluster set within the corridor from the original three-dimensional coordinate system to a local coordinate system of the corridor with the corridor direction as the vertical axis, the perpendicular corridor direction as the horizontal axis, and the elevation as the vertical axis, thereby generating a local coordinate point set for each candidate point cluster.

[0080] S34: In the local coordinate system of the corridor, the local coordinate point set of each candidate point cluster is fitted with an arc curve along the vertical axis. The initial parameters of the fitted model are obtained by using the random sampling consensus algorithm. Then, the optimal model parameters are solved by least squares optimization. Based on the optimal model parameters, the root mean square error of the fitting residual and the lateral width of each candidate point cluster are calculated. The candidate point clusters with a root mean square error of the fitting residual less than a preset residual threshold and a lateral width less than a preset width threshold are confirmed as point clusters that conform to the arc shape of the power line, and the confirmed power line residual point clusters are generated.

[0081] Specifically, the system reads the tower location data of the survey area, which includes the three-dimensional spatial coordinates of each tower in the global coordinate system and the topological connection relationship between the towers. Based on the tower location data, the system extracts tower pairs with adjacent topological relationships, using the spatial line connecting two adjacent towers as the corridor's central axis. The direction vector of this central axis is determined by the three-dimensional coordinate difference between the towers at both ends. Using the corridor's central axis as the center, the system expands to both sides of a set width in the horizontal plane along a direction perpendicular to the central axis, constructing the power line corridor spatial range corresponding to each span. When constructing the corridor spatial range, the system adopts a segmented variable width strategy along the corridor direction. A basic expansion width is used at the tower positions at both ends of the central axis, and an expansion width is used in the area near the midpoint of the central axis. This width variation corresponds to the spatial shape of the power line sag, allowing the corridor to have a wider spatial coverage in the area with the largest sag. The system aggregates the corridor spatial ranges corresponding to each span into a set of corridor spatial constraint regions. The corridor width function used by the system during the construction process is:

[0082]

[0083] in, The arc length measured along the central axis of the corridor from the starting tower is expressed in meters. The basic extension width is measured in meters. This is a dimensionless adjustment coefficient. The span length between adjacent towers is measured in meters. The dimension of gravitational acceleration is meters per square second. The horizontal stress on the conductor is expressed in Newtons per square meter. The weight per unit length of the conductor is measured in Newtons per meter. The square root term in the formula is measured in meters, and this term represents the theoretical sag depth of the electric field line at the midpoint of the span. The square sinusoidal term is... and The value at this location is zero and dimensionless. The value at this point is dimensionless, and the entire second term has a dimension of meters, which is the same as the base width. The sum is expressed in meters. This width function causes the corridor space to narrow at the towers and widen in the middle of the span, matching the actual spatial distribution of the power lines.

[0084] Preferably, such as Figure 2 As shown, this application provides a constrained diagram of the candidate range of a pole corridor and power lines, which shows the spatial structure of the power line corridor within a single span. The poles are located at both ends of the corridor, serving as the start and end points of the corridor's central axis. The corridor's central axis is marked as a dashed line between two poles, representing the theoretical projection path of the power lines on the horizontal plane. Figure 2 The corridor direction is further indicated by the direction vector from the starting tower to the ending tower, which is used to indicate the orientation of the vertical axis when the candidate point cluster coordinates are transformed. Figure 2 The candidate range points are marked on both sides of the corridor's central axis. These points represent the candidate point cluster members that fall within the corridor's spatial constraints, and their spatial distribution reflects the degree of deviation of the candidate points from the corridor's central axis. Figure 2 The search buffer shown is an extended area on both sides of the central axis of the corridor. This area covers the space that the power line may occupy due to factors such as sag and wind deflection. Only candidate points that fall within this buffer are eligible to enter the subsequent sag fitting process. Figure 2 The profile constraint is marked at a specific location on the corridor's central axis. This location corresponds to the elevation profile extracted along the corridor direction during subsequent steps of sag model fitting. The profile constraint limits the effective data range along the longitudinal axis during sag fitting. Figure 2 As shown in the spatial constraint relationship, the system transforms the spatial attribution judgment of candidate point clusters from a globally unconstrained state to a locally constrained state along the corridor direction. This ensures that the sag fitting operation in subsequent steps is only performed on candidate point clusters that meet the prior conditions of the power line engineering location in space, thereby improving the spatial specificity of power line residual point identification.

[0085] The system performs spatial location determination on each member point in the candidate point cluster set, determining whether the point is located within any corridor space range of the constructed corridor spatial constraint region set. The system counts the number of points in each candidate point cluster that fall within the corridor space range, calculates the proportion of this number to the total number of points in the candidate point cluster, and uses this proportion as the corridor affiliation score for that candidate point cluster. When calculating the affiliation score, the system assigns a contribution weight to each member point that varies with its spatial location. Points located near the corridor's central axis receive a higher contribution weight, while points located at the corridor's edge receive a lower contribution weight, reflecting the difference in confidence in the affiliation determination based on different spatial locations. The system compares the corridor affiliation score of each candidate point cluster with a preset corridor affiliation threshold, selecting candidate point clusters with scores greater than the threshold and aggregating these candidate point clusters into a corridor-wide candidate point cluster set. The system excludes candidate point clusters with scores not greater than the threshold from subsequent processing. The weighted affiliation scoring formula used by the system in the scoring process is as follows:

[0086]

[0087] in, For the first A cluster of candidate points, For a member point in this point cluster, The horizontal distance from this point to the central axis of the corridor is measured in meters. The distance bandwidth parameter is measured in meters. For indicator functions, when point When a point is within the corridor space, its value is dimensionless; when a point is outside the corridor space, its value is zero (dimensionless). In the formula, the numerator is the sum of the products of the Gaussian weights and the indicator values ​​of each member point; its dimension is determined by the dimensionless nature of the indicator function, making it dimensionless overall. The denominator is the Gaussian weights and the dimensionless value; dividing the two yields... Dimensionless. This scoring formula assigns higher weights to points closer to the corridor's central axis and lower weights to points closer to the corridor's boundaries. This means that even if some points fall near the corridor boundaries due to positioning errors, their contribution to the score is reduced accordingly, thus improving the robustness of the attribution determination.

[0088] Furthermore, based on the direction of the corridor's central axis to which each candidate point cluster belongs, the system performs coordinate transformation on each candidate point cluster in the corridor's candidate point cluster set. The system uses the direction of the corridor's central axis to which the point cluster belongs as the vertical axis of the local coordinate system, the direction perpendicular to the corridor's central axis in the horizontal plane as the horizontal axis of the local coordinate system, and the elevation axis of the global coordinate system as the vertical axis of the local coordinate system. The system uses the starting tower position of the corridor as the origin of the local coordinate system and constructs a rotation and translation transformation matrix from the global coordinate system to the local coordinate system. The system substitutes the global coordinates of each member point in the candidate point cluster into the transformation matrix to calculate its coordinate values ​​in the local coordinate system, including the vertical projection distance along the corridor's direction, the horizontal offset distance from the corridor's central axis, and the elevation value, generating the local coordinate point set of the candidate point cluster. During the coordinate transformation process, the system normalizes the projection distance of each member point along the corridor's vertical axis; the calculation formula is as follows:

[0089]

[0090] in, The spatial location vector of this member point has units of meters. The spatial position vector of the pole at the starting point of the corridor has units of meters. Let the direction vector of the corridor's central axis have dimensions in meters, the numerator be the dot product of the two vectors with dimensions in square meters, and the denominator be the square of the magnitude of the direction vector with dimensions in square meters. Divide the two and... The unit of measurement is meters. This normalized projected distance represents the relative position of member points along the corridor's central axis, eliminating the influence of different span lengths on the vertical axis scale. The system simultaneously calculates the horizontal axis offset of each member point. The calculation formula is: The ratio of the cross product magnitude to the direction vector magnitude to the central axis, in meters, and the elevation coordinates. The elevation value of the global coordinates is directly taken in meters.

[0091] The system fits a sag curve along the vertical axis to the local coordinate set of each candidate point cluster in the local coordinate system of the corridor. The system uses a random sampling consensus algorithm to obtain the initial parameters of the fitted model. This algorithm randomly selects the smallest sample subset from the candidate point clusters, solves for the initial parameter values ​​of the model based on this subset, calculates the distance from all member points to the initial model, counts the number of member points with a distance less than a set threshold, and selects the model parameters containing the most consensus points as the initial parameters. Using the initial parameters obtained by the random sampling consensus algorithm as the starting point for iteration, the system uses a least squares optimization method to finely solve for the model parameters. Through iterative adjustment of the model parameters, the objective function value is gradually reduced until convergence. The weighted objective function used by the system in the optimization process is:

[0092]

[0093] in, The number of points within the candidate point cluster is dimensionless. For point index, and The first The local ordinate and elevation values ​​of each point are measured in meters. For model parameter vectors, The output dimension of the sag model function is meters. For the first The weight of each point is dimensionless and is determined based on the density of its local neighborhood. The dimensionless dimension of the sag curvature constraint coefficient is meters. Let the center position of the candidate point cluster along the vertical axis be in meters. The first term of the objective function is the weighted average of the squared residuals of each member point, in square meters. The second term is the absolute value of the second derivative of the model at the center position multiplied by the constraint coefficient, in square meters. The sum of the two terms gives the objective function in square meters.

[0094] Preferably, the system employs the Levenberg-Marquardt algorithm to iteratively optimize the objective function. The iteration terminates when the parameter change is less than a set value or the maximum number of iterations is reached. The system calculates the root mean square error (RMSE) of the fitting residuals and the lateral width of each candidate point cluster based on the optimal model parameters. The RMSE is obtained by taking the square root of the first term of the objective function, with units of meters. The lateral width is determined by the difference between the maximum and minimum values ​​of the horizontal axis coordinates of the member points, also with units of meters. The system identifies candidate point clusters with RMSEs less than a preset residual threshold and lateral widths less than a preset width threshold as point clusters conforming to the sag shape of power lines, generating confirmed power line residual point clusters. The system labels candidate point clusters with different sag shapes as linear model matching, parabolic model matching, or catenary model matching, and records the model type and optimized model parameters.

[0095] Preferably, during the fitting process, the system calculates the distance from the point to the model, the continuous length of the point cluster, the lateral width, the piecewise residual, the sag direction, and the height range. If a candidate point cluster meets the conditions of small lateral width, continuity along the corridor direction, a vertical trend conforming to a sag or linear suspension shape, and a stable height difference with the tree canopy or roof surface below, it is confirmed as a power line remnant. If the point cluster is linear but close to the roof edge, tree top, or road facilities, and does not conform to the sag model, it is marked as a non-power line candidate or an object for manual verification. Figure 3 As shown in the figure, the sag model fitting and residual point confirmation diagram provided in this application demonstrate that the residual points of power lines exhibit a slender, continuous spatial morphology with a sag trend between towers. Through model fitting, the actual residual power lines can be distinguished from linear interferences such as tree branches and roof edges.

[0096] In one embodiment, the preset anomaly verification conditions include morphological anomaly verification conditions, profile elevation anomaly verification conditions, and neighborhood elevation support verification conditions. S4 of the automatic power line residue filtering method for airborne point cloud (DSM) production provided by this invention specifically includes the following steps:

[0097] S41: Rasterize the DSM representation point set with a preset DSM grid resolution, count the elevation values ​​of all points falling into each grid cell and take the maximum value as the elevation value of the corresponding grid cell to generate a DSM raster.

[0098] S42: Map the three-dimensional coordinates of each point in the confirmed power line residual point cluster to the grid cell index of the DSM grid, mark the mapped grid cells as anomalous grids to be verified, and generate a set of anomalous grids to be verified.

[0099] S43: Perform morphological anomaly verification on the set of abnormal grids to be verified: Calculate the aspect ratio of the maximum extension length along the corridor direction to the maximum width perpendicular to the corridor direction of the set of abnormal grids to be verified, filter out abnormal grids with aspect ratio values ​​less than the preset aspect ratio threshold, and generate an abnormal grid set that passes the morphological verification.

[0100] S44: Perform profile elevation anomaly verification on the set of anomalies that have passed the morphological verification. Extract the DSM elevation profile sequence at each anomaly location along the vertical corridor direction. Calculate the first difference value of the elevation profile sequence near the lowest point of the sag. If the first difference value shows a peak shape that changes from positive to negative and the peak width is less than the preset width threshold, it is determined to pass the profile verification and generate an anomaly grid set that has passed the profile verification.

[0101] S45: Perform neighborhood elevation support verification on the set of anomalous grids that have passed profile verification. Determine the point cloud neighborhood corresponding to each anomalous grid according to the spatial neighborhood index relationship. After mapping the point cloud neighborhood to the DSM grid neighborhood, calculate the elevation standard deviation of the neighborhood grids. Calculate the comprehensive anomaly verification index of each anomalous grid. If the comprehensive anomaly verification index of the anomalous grid is greater than the preset comprehensive judgment threshold, it is determined to have passed the neighborhood support verification. The points corresponding to the anomalous grids that have passed the morphological verification, profile elevation anomaly verification, and neighborhood elevation support verification are identified as power line remnants. Remove the power line remnants from the DSM expression point set to generate an updated DSM expression point set and void regions.

[0102] Specifically, the system rasterizes the DSM representation point set using a preset DSM grid resolution, which is consistent along both the horizontal and vertical axes of the horizontal plane, forming square grid cells. The system calculates the index value of the grid cell to which each point belongs based on its planar coordinates, mapping each point in the DSM representation point set to its corresponding grid cell. The system performs elevation value statistics on all points within each grid cell, extracting the maximum elevation value and assigning it to the corresponding grid cell as its elevation value. During the statistical process, the system performs a neighborhood consistency check on the elevation values ​​of each point. If the elevation value of a point significantly deviates from the median elevation of its neighboring points by more than a set multiple, the point is marked as a potential outlier. The system reduces the weight of the elevation value of this potential outlier before participating in the maximum value extraction, thus suppressing the interference of isolated noise points on the grid elevation values. The system uses the following formula to allocate the weights of each grid cell:

[0103]

[0104] in, For the first in the grid The elevation values ​​of each point are in meters. The median elevation of the nearest neighbor point in the neighborhood of this point, in meters. The standard deviation of the neighborhood elevation values ​​is expressed in meters. This is a dimensionless scale parameter. The molecular weight in the exponential part of the weighted term is in square meters, the denominator is in square meters, and the overall exponential part is dimensionless. When extracting the maximum elevation value of the grid, the system multiplies each point's elevation value by its corresponding weight before comparison, thus suppressing the contribution of points whose elevation values ​​significantly deviate from their neighborhood distribution to the grid value. After extracting the elevation values ​​of all grid cells, the system arranges the grid cells according to their planar index to generate a DSM grid, which covers the entire spatial range of the survey area.

[0105] The system confirms that the 3D coordinates of each point in the power line remnant cluster are mapped to the grid cell index of the DSM raster. Based on the planar coordinates of each point, the system calculates the row and column index of the grid cell to which that point belongs, according to the resolution and origin position of the DSM raster. The system performs coordinate-to-index conversion on each member point in each confirmed power line remnant cluster. When multiple points are mapped to the same grid cell, the system marks that grid cell once, avoiding duplicate marking. The system extracts all mapped grid cells to form a set of anomaly grids to be verified. During the mapping process, the system evaluates the mapping confidence level, assigning higher confidence to points located in the center region of the confirmed power line remnant cluster and lower confidence to points located on the cluster boundaries. The mapping confidence level is calculated using the following formula:

[0106]

[0107] in, To confirm that the shortest distance from the residual point of the power line to the boundary of its cluster is measured in meters, Let the span of this point cluster along the principal direction be in meters. The standard deviation of the spatial location of the point cluster members is expressed in meters. This is a dimensionless adjustment coefficient. In the formula... The ratio to the denominator is dimensionless, and the confidence value is also dimensionless, ranging from zero to one. For grid cells with a mapping confidence level below a set threshold, the system adds a confidence label to the set of abnormal grid cells to be verified, allowing subsequent verification steps to adjust the severity of the verification based on the confidence level.

[0108] Further, the system determines the corridor direction of each grid cell in the set of anomalous grid cells to be verified. This direction information is derived from the recorded corridor centerline direction vector. The system performs spatial distribution analysis on the anomalous grid cell set along its corridor direction, calculating the maximum extension length of the set in the corridor direction and the maximum width in the direction perpendicular to the corridor. The system uses the ratio of the maximum extension length to the maximum width as the aspect ratio feature value of the anomalous grid cell set. The system uses morphological opening operations to preprocess the anomalous grid cell set to eliminate the interference of isolated anomalous grid cells caused by point cloud sparsity on the morphological feature calculation. The system compares the aspect ratio feature value with a preset aspect ratio threshold, filtering out anomalous grid cells with aspect ratio values ​​less than the threshold from the set. The system performs weighted processing on the extension length during the aspect ratio calculation process, and the calculation formula is as follows:

[0109]

[0110] in, The maximum extension length of the anomalous grid set along the corridor direction is in meters. The actual number of grid cells in the abnormal grid set is dimensionless. The total number of grid cells within the bounding rectangle of the anomalous grid set is dimensionless, and the ratio between the two is also dimensionless. The dimensionless filling gain coefficient. The maximum width of the anomalous grid set perpendicular to the corridor direction is in meters. The values ​​within square brackets in the numerator are dimensionless. The product is multiplied to the dimension of meters, then divided by... The unit of measurement is meter. Dimensionless. This weighting method enables anomalous grid sets with continuous morphology and complete internal filling to achieve a higher aspect ratio, increasing the probability that dense electric field line remnant regions will pass morphological verification.

[0111] Furthermore, the system performs profile elevation anomaly verification on the set of anomalous grids that have passed morphological verification. The system determines the planar position of each anomalous grid within the DSM grid and extracts a DSM elevation profile sequence at each anomalous location along the vertical corridor direction. This sequence covers several grid cells on both sides of the anomalous location. The system obtains the elevation values ​​of each grid cell in this profile sequence, forming a one-dimensional profile elevation array. Based on the sag model parameters of the confirmed power line remnant cluster at the anomalous location, the system determines the theoretical position of the lowest sag point. Near this position, the system calculates the first-order difference value of the profile elevation sequence. The first-order difference is the ratio of the difference between adjacent grid elevation values ​​to the grid spacing. The system judges the distribution pattern of the first-order difference value along the spatial location. When the difference value sequence successively shows positive values ​​followed by negative values, forming an alternating positive and negative peak pattern, and the width of this peak pattern along the profile is less than a preset width threshold, the system determines that the anomalous location has passed the profile elevation anomaly verification. The system's quantitative calculation formula for the peak amplitude is:

[0112]

[0113] in, The first-order difference sequence at the anomalous location has units of meters. The maximum value in the difference sequence has the dimension of meters. Let the minimum value in the difference sequence be in meters, and the difference between the two values ​​be in meters. The standard deviation of the difference sequence is measured in meters. The average value of the difference sequences has a dimension of meters, and the ratio of the two is dimensionless. This is a dimensionless adjustment coefficient. The exponent term in the formula is dimensionless, while the molecular weight is in meters. Dividing the two yields... The dimension is meters. This quantitative indicator reflects the significance of the peak amplitude of the difference sequence relative to the level of local fluctuations; the larger the peak amplitude, the higher the quantitative value.

[0114] The system performs neighborhood elevation support verification on the set of anomalous grids verified by profile analysis. Based on the spatial neighborhood index relationship established in step S11, the system determines the point cloud neighborhood corresponding to each anomalous grid, maps the planar coordinates of each point in the point cloud neighborhood to the DSM grid, and obtains the corresponding set of grid neighborhood cells. The system calculates the standard deviation of the elevation values ​​of each grid cell within the neighborhood as a quantitative indicator of the degree of neighborhood elevation fluctuation. The system combines the aspect ratio obtained in step S43, the peak difference obtained in step S44, and the neighborhood elevation standard deviation obtained in this step to calculate the comprehensive anomaly verification index for each anomalous grid. The formula for calculating the comprehensive anomaly verification index is:

[0115]

[0116] in, For abnormal grid The comprehensive anomaly verification index is dimensionless and ranges from zero to one. This is the weighting coefficient for the morphological verification index. This coefficient is dimensionless and ranges from 0.2 to 0.4. . For abnormal grid The aspect ratio of , which is dimensionless. This is a preset maximum aspect ratio reference value, which is dimensionless and used to... Normalized to the interval between zero and one, it becomes dimensionless after normalization. This is the weighting coefficient for the profile verification index; this coefficient is dimensionless. For abnormal grid The peak value of the first-order difference of the profile elevation, which is measured in meters. This is a preset maximum reference value for the profile differential peak value, measured in meters, used to... Normalized to the interval between zero and one, it becomes dimensionless after normalization. The weighting coefficient for the neighborhood support verification index is dimensionless. For abnormal grid The standard deviation of the neighborhood grid elevation, which is measured in meters. This is the preset maximum reference value for the standard deviation of elevation in the neighborhood, measured in meters, used to... Normalized to the interval between zero and one, it becomes dimensionless after normalization.

[0117] The system compares the comprehensive anomaly verification index with a preset comprehensive judgment threshold. Anomaly grids with an index greater than the threshold are judged to have passed neighborhood support verification. The system identifies the points corresponding to anomaly grids that simultaneously pass morphological verification, profile elevation anomaly verification, and neighborhood elevation support verification as power line remnants. The system removes these power line remnants one by one from the DSM expression point set, generating an updated DSM expression point set. Based on the spatial location of the removed grids in the DSM raster, the system marks the corresponding areas as void areas. These void areas, along with the updated DSM expression point set, are then included in the subsequent surface repair process.

[0118] Preferably, in steps S43 to S45, during the process of the system sequentially performing morphological anomaly verification, profile elevation anomaly verification, and neighborhood elevation support verification on the anomaly grid set to be verified, the corresponding Figure 4 The provided DSM narrowband linear elevation anomaly verification map shows the spatial morphological characteristics of the DSM narrowband linear elevation anomaly verification. Figure 4 The image shows typical anomalies formed by residual power line points in a DSM grid. Linear spikes are marked at specific locations in the DSM profile. These spikes are local protrusions in narrow bands extending along the corridor direction, where the grid elevation value is significantly higher than the surrounding background value. Their spatial distribution width corresponds to the physical diameter of the power line and the spatial resolution of the point cloud. They are the direct objects identified by calculating the aspect ratio in the morphological anomaly verification step. Figure 4 The protrusion in the profile is presented in a cross-sectional view perpendicular to the corridor direction. In the profile elevation sequence, the protrusion is a peak formed by the first-order difference value turning from positive to negative. Its peak amplitude and peak width characteristics correspond to the differential peak quantification index and width threshold judgment basis in the profile elevation anomaly verification, respectively. Figure 3 The neighborhood support check is presented in the form of anomaly grid and its neighboring grids. This check quantitatively assesses the stability of the surface around the anomaly grid by calculating the standard deviation of the neighboring grid elevation. When the elevation of the anomaly grid is significantly higher than the neighborhood mean and the neighborhood standard deviation is in a small range, it indicates that the anomaly grid does not have the elevation transition characteristics to blend with the surrounding stable surface. Figure 4The abnormal dragging in the data indicates a spatial tailing effect caused by the participation of power line remnants in DSM gridding, which leads to an abnormal increase in the elevation values ​​of adjacent grids. This effect causes the linear anomaly region to extend continuously along the direction of the power line.

[0119] pass Figure 4 The spatial morphological characteristics of the narrow-band linear elevation anomaly verification in DSM are shown. The system decomposes the manifestation of power line remnants on the DSM result surface into four verification dimensions: morphological characteristics of linear peaks, differential characteristics of profile protrusions, statistical characteristics of neighborhood support, and extension characteristics of anomaly dragging. This enables the morphological anomaly verification in step S43, the profile elevation anomaly verification in step S44, and the neighborhood elevation support verification in step S45 to obtain quantitative judgment basis with clear geometric correspondence in the DSM grid space.

[0120] In one embodiment, S5 of the automatic power line residue filtering method for airborne point cloud (DSM) production provided by the present invention specifically includes the following steps:

[0121] S51: Based on the updated DSM expression point set and the cavity region, extract the effective grid cells and corresponding elevation values ​​within a preset range around the cavity region as repair reference data, and generate a cavity neighborhood reference dataset.

[0122] S52: Based on the geometric morphological characteristics of the cavity region, the cavity region is repaired. If the cavity region is distributed as isolated points and its diameter is less than the preset aperture threshold, the inverse distance weighted interpolation method is used to repair the cavity region. The weighted average is calculated using the inverse square of the distance from each effective grid cell in the cavity neighborhood reference dataset to the center of the cavity grid as the weight to generate the interpolated repair elevation value. If the cavity region is distributed as a linear extension, the effective elevation profiles on both sides are extracted along the cavity extension direction. The cubic spline interpolation method is used to interpolate the cavity grid between the profiles point by point to generate the spline repair elevation value. If the cavity region is distributed as a planar feature and its diameter is greater than or equal to the preset aperture threshold, the effective grid cells in the cavity neighborhood reference dataset are used to construct a quadratic surface equation and the surface coefficients are solved by the least squares method. The planar coordinates of the cavity grid are substituted into the quadratic surface equation to calculate the surface repair elevation value.

[0123] S53: Write the interpolated repair elevation values, spline repair elevation values, and surface repair elevation values ​​of each hole grid into the corresponding positions of the DSM raster, generate the repaired DSM raster, and record the original elevation values ​​and types of repair elevation values ​​of each repair grid to generate a repair area record layer.

[0124] S54: Perform residual anomaly detection processing on the repaired DSM raster based on the repair area record layer. Perform anomaly re-inspection on the repaired area and its neighborhood in the repaired DSM raster. Calculate the proportion of the number of grids with still abnormal features to the total number of DSM grids as the anomaly residual rate. If the anomaly residual rate is greater than the preset residual rate threshold, reduce the preset linearity threshold and preset suspension height threshold and re-execute S2 to S5. If the anomaly residual rate is less than or equal to the preset residual rate threshold and the gradient change of the elevation transition in the repaired area is less than the preset gradient threshold, the quality is judged to be qualified, and the final DSM result is output. The final DSM result is used to indicate the digital surface model result that meets the quality requirements after power line residue filtering.

[0125] Specifically, based on the updated DSM expression point set and the hole region, the system performs spatial extent analysis on the hole region to determine the planar location and coverage of each hole grid. The system sets a neighborhood search radius centered on each hole grid and retrieves valid grid cells belonging to the updated DSM expression point set within this radius. The system extracts the planar coordinates and corresponding elevation values ​​of these valid grid cells and compiles the results into a hole neighborhood reference dataset. When retrieving valid grid cells, the system adjusts the neighborhood search radius based on the local point cloud density distribution characteristics of the hole region, using a smaller search radius for holes in dense point cloud regions and a larger search radius for holes in sparse point cloud regions to ensure that each hole obtains a sufficient number of repair reference data points. When constructing the reference dataset, the system evaluates the reference contribution of each valid grid cell. The contribution is determined by combining the distance of the valid grid cell to the center of the hole grid and the elevation confidence of the grid cell in the local region; grid cells that are closer in distance and have higher confidence receive higher contribution labels. The system associates and stores the reference dataset of the cavity neighborhood with the spatial location of the cavity area, so that each cavity grid corresponds to a set of spatially adjacent effective reference grid cells and their elevation values. This association serves as the data basis for subsequent repair processing.

[0126] Furthermore, the system classifies void regions based on their geometric morphological characteristics and performs corresponding repair processes: First, the system calculates the geometric dimensions and morphological parameters of each void region, including the diameter of its circumscribed circle, its extension length in each direction, and its aspect ratio. Based on these parameters, void regions are classified into three types: isolated point distribution, linear extension distribution, and area distribution. For isolated point distribution void regions with a diameter smaller than a preset aperture threshold, the system uses inverse distance weighted interpolation for repair. The system uses the inverse square of the distance from each effective grid cell in the void's neighborhood reference dataset to the center of the void grid as the weight to perform a weighted average of the elevation values ​​of each effective grid cell, generating an interpolated repair elevation value. The calculation formula is as follows:

[0127]

[0128] in, The total number of valid grid cells in the hollow neighborhood reference dataset is dimensionless. For effective grid cell indexing, For the first The planar distance from an effective grid cell to the center of a void grid is measured in meters. The smallest positive quantity with dimensions in meters is used to prevent the denominator from being zero. For the first The elevation value of each effective grid cell is measured in meters. The arithmetic mean of the distances from all valid grid cells to the cavity center is expressed in meters. In the formula, the numerator of the first fraction is the weighted sum of the elevation values ​​of each valid grid cell multiplied by the inverse square of the distance, also in meters; the denominator is the sum of the inverse squares of the distances, which is dimensionless. The overall dimension of the fraction is meters. The second term within square brackets is the normalized value of the distance variance, which is dimensionless and is used to correct the elevation value after multiplying the whole term by the first fraction. The unit of measurement is meters. This correction term appropriately reduces the interpolation result when the reference points are widely distributed, avoiding excessive influence of distant reference points on the repair value. For linearly extending cavity areas, the system extracts effective elevation profiles on both sides along the cavity extension direction, and uses cubic spline interpolation to interpolate the cavity grid between the profiles point by point to generate spline repair elevation values. The cubic spline interpolation uses the two profile points as control nodes and constructs a piecewise cubic polynomial function, so that the interpolation curve has continuous first and second derivatives at each control node, thereby ensuring a smooth transition of the repair elevation along the profile direction. For void regions with a planar distribution and a diameter greater than or equal to a preset aperture threshold, the system constructs a quadratic surface equation using effective grid cells in the void neighborhood reference dataset. The equation is a quadratic polynomial in planar coordinates, containing six undetermined coefficients: a constant term, a first-order term, a second-order term, and a cross term. The system solves for these six coefficients using the least squares method and substitutes the planar coordinates of the void grid into the solved quadratic surface equation to calculate the surface repair elevation value.

[0129] The system, based on the morphological type of each grid cell in the cavity region and the corresponding repair method, writes the interpolated repair elevation values, spline repair elevation values, and surface repair elevation values ​​calculated in step S52 into the corresponding grid positions of the DSM raster. The system iterates through each grid cell within the cavity region, determining the appropriate repair elevation value based on the morphological type of the cavity region to which the grid cell belongs, and replacing the original elevation value of that grid cell in the DSM raster with the corrected elevation value. During the write operation, the system simultaneously records the processing information for each repair grid cell, including the row and column index of the grid cell, the original elevation value before repair, the elevation value after repair, and the repair method identifier. The repair method identifier corresponds to inverse distance weighted interpolation, cubic spline interpolation, and quadratic surface fitting, respectively. The system aggregates the processing information of all repair grid cells into a repair area recording layer, which is stored in raster form, and its spatial coverage is completely consistent with the cavity region. When generating the repair area record layer, the system calculates the repair confidence score for each repair grid. The formula for calculating the repair confidence score is as follows:

[0130]

[0131] in, The standard deviation of the effective grid cell elevation values ​​in the reference dataset for the hole neighborhood is expressed in meters. The dimension of the range of effective grid cell elevation values ​​is meters, and the ratio between the two is dimensionless. The total number of grid cells in the current void region is dimensionless. The total number of effective grid cells in the neighborhood reference dataset is in meters. In the first fraction of the formula, the numerator is dimensionless, and the denominator is composed of one and the sum of ratios, also dimensionless. The overall dimensionless value of the fraction ranges from 0.5 to 1. The term within square brackets is the dimensionless neighborhood sufficiency factor; this value approaches one when there are sufficient neighborhood reference points and approaches zero when reference points are scarce. The system associates the repair confidence label with the corresponding grid repair record and stores it in the repair area record layer.

[0132] Furthermore, the system performs residual anomaly detection processing on the repaired DSM raster based on the repair area record layer. The system identifies the repaired grid cells and their neighboring grid cells in the repair area record layer as the detection focus area, and re-executes the anomaly detection process in this area, including three verification operations: morphological anomaly verification, profile elevation anomaly verification, and neighboring elevation support verification. The system then counts the number of grid cells within this area that still exhibit anomalous characteristics. The system calculates the proportion of these anomalous grid cells to the total number of DSM grid cells as the anomaly residual rate. Simultaneously, it calculates the gradient change in elevation transition between each repaired grid cell and its neighboring valid grid cells within the repair area. The formula for calculating the gradient change is:

[0133]

[0134] in, To repair the dimensionless total number of grid cells within the region, To repair the grid index, and The partial derivatives of the repaired DSM raster along the horizontal and vertical axes are respectively expressed in meters per grid spacing. The square root of the sum of the squares of the two terms is also expressed in meters per grid spacing. The elevation values ​​for repairing the grid are in meters. and These represent the maximum and minimum effective grid elevations of the neighborhood, both in meters. The term within square brackets represents the ratio of the repaired elevation to the neighborhood elevation range; this value is dimensionless. Overall gradient change. The unit is meters per grid spacing. The system compares the abnormal residual rate with a preset residual rate threshold. When the abnormal residual rate exceeds this threshold, the system lowers the linearity threshold and the suspension height threshold, and re-executes candidate point extraction, corridor screening, sag fitting, DSM anomaly verification, and removal operations with the adjusted threshold parameters. When the abnormal residual rate is not greater than this threshold and the gradient change in elevation transition of the repair area is less than the preset gradient threshold, the system determines that the DSM quality is qualified. The system outputs the final DSM raster that meets the quality requirements. This final DSM result is used to indicate the digital surface model result that meets the quality requirements after power line residual filtering.

[0135] It should be understood that although the steps in the flowcharts of the embodiments described above are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the flowcharts of the embodiments described above may include multiple steps or multiple stages. These steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least some of the steps or stages of other steps.

[0136] Based on the same inventive concept, this application also provides an automatic power line residue filtering system for airborne point cloud DSM production, which implements the above-described automatic power line residue filtering method for airborne point cloud DSM production. The solution provided by this system is similar to the implementation described in the above method. Therefore, the specific limitations of the one or more embodiments of the automatic power line residue filtering system for airborne point cloud DSM production provided below can be found in the limitations of the automatic power line residue filtering method for airborne point cloud DSM production described above, and will not be repeated here.

[0137] Preferably, such as Figure 5 As shown, the present invention provides an automatic power line residue filtering system 600 for airborne point cloud (DSM) production, which is configured with the following modules:

[0138] The DSM feature point construction module 610 is used to read the preprocessed and initially classified airborne LiDAR point cloud data and DSM expression category rules, extract the points participating in DSM expression from the airborne LiDAR point cloud data according to the DSM expression category rules, construct the DSM expression point set, perform local neighborhood feature analysis on each point in the DSM expression point set, calculate the geometric feature index of each point, and generate a DSM expression point set with local geometric features.

[0139] The candidate point cluster extraction module 620 is used to perform feature screening based on the geometric feature indicators of each point in the DSM expression point set, extract points that simultaneously meet the preset geometric feature conditions as candidate points, and perform spatial neighborhood clustering on the candidate points to group the spatially continuously distributed candidate points into the same cluster, generating a candidate point cluster set.

[0140] The power line cluster identification module 630 is used to read the tower location data of the survey area, construct the power line corridor space range between each tower based on the tower location data, project each candidate point cluster in the candidate point cluster set onto the power line corridor space range for classification, filter out the candidate point clusters that fall within the corridor range, and perform spatial morphology fitting on each point cluster in the candidate point clusters that fall within the corridor range along the corridor direction, calculate the morphological feature parameters of each point cluster, confirm the point clusters that conform to the power line spatial morphology based on the morphological feature parameters, and generate confirmed power line residual point clusters;

[0141] The power line point removal module 640 is used to construct a DSM grid based on the DSM expression point set, project the confirmed power line residual point clusters onto the DSM grid, perform anomaly verification on the grid cells at the projection position, confirm the points corresponding to the grid cells that pass the anomaly verification as power line residual points, remove the power line residual points from the DSM expression point set, and generate an updated DSM expression point set and hole region.

[0142] The DSM iterative repair module 650 is used to perform surface repair processing on the void regions based on the updated DSM expression point set and void regions by performing neighborhood interpolation or surface fitting, generating a repaired DSM raster, and performing quality evaluation on the repaired DSM raster by residual anomaly detection. Based on the quality evaluation results, the threshold parameters for candidate point selection are updated and S2 to S5 are re-executed until the quality evaluation results meet the preset quality requirements, and the final DSM result is output.

[0143] In one embodiment, this application also provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the above-described automatic power line residue filtering method for airborne point cloud (DSM) production.

[0144] In one embodiment, this application also provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the above-described method for automatic filtering of power line residues in airborne point cloud (DSM) production.

[0145] In the description of this specification, references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of this application. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of those different embodiments or examples.

[0146] For the device embodiments, since they basically correspond to the method embodiments, the relevant parts can be referred to in the description of the method embodiments. The device embodiments described above are merely illustrative. The components described as separate parts may or may not be physically separate, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this disclosure according to actual needs. Those skilled in the art can understand and implement this without creative effort.

[0147] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any person skilled in the art can easily conceive of various variations or substitutions within the technical scope disclosed in this application, and these should all be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for automatic filtering of residual power lines in airborne point cloud (DSM) production, characterized in that, Includes the following steps: S1: Read the preprocessed and initially classified airborne LiDAR point cloud data and DSM expression category rules, extract the points participating in DSM expression from the airborne LiDAR point cloud data according to the DSM expression category rules, construct the DSM expression point set, perform local neighborhood feature analysis on each point in the DSM expression point set, calculate the geometric feature index of each point, and generate a DSM expression point set with local geometric features; S2: Based on the geometric feature indicators of each point in the DSM expression point set, feature screening is performed, and points that simultaneously meet the preset geometric feature conditions are extracted as candidate points. Spatial neighborhood clustering is performed on the candidate points, and the candidate points that are continuously distributed in space are grouped into the same cluster to generate a candidate point cluster set. S3: Read the tower location data of the survey area, construct the power line corridor space range between each tower based on the tower location data, project each candidate point cluster in the candidate point cluster set onto the power line corridor space range for classification, filter out the candidate point clusters that fall within the corridor range, and perform spatial morphology fitting on each point cluster in the candidate point clusters that fall within the corridor range along the corridor direction, calculate the morphological feature parameters of each point cluster, confirm the point clusters that conform to the power line spatial morphology based on the morphological feature parameters, and generate confirmed power line residual point clusters; S4: Construct a DSM grid based on the DSM expression point set, project the confirmed electric field line residual point cluster onto the DSM grid, perform anomaly verification on the grid cell at the projection position based on preset anomaly verification conditions, confirm the point corresponding to the grid cell that passes the anomaly verification as electric field line residual point, remove the electric field line residual point from the DSM expression point set, and generate an updated DSM expression point set and hole region; S5: Based on the updated DSM expression point set and the void region, perform surface repair processing on the void region using neighborhood interpolation or surface fitting to generate a repaired DSM raster. Then, perform quality evaluation on the repaired DSM raster using residual anomaly detection. Update the threshold parameters for candidate point selection based on the quality evaluation results and re-execute S2 to S5 until the quality evaluation results meet the preset quality requirements, and output the final DSM result.

2. The method according to claim 1, characterized in that, S1 includes: S11: For each point cloud point in the preprocessed and initially classified airborne LiDAR point cloud data, a spatial neighborhood retrieval data structure is constructed using the spatial distribution position of the three-dimensional coordinates of the airborne LiDAR point cloud data as an index. The spatial neighborhood retrieval data structure is used to find the set of neighboring points of each point cloud point within a preset radius, and the spatial neighborhood index relationship of each point is generated. S12: Based on the spatial neighborhood index relationship, extract the three-dimensional coordinate matrix of the neighborhood point set for each point in the DSM expression point set. After centering the three-dimensional coordinate matrix, construct the covariance matrix. Obtain three eigenvalues ​​by performing eigenvalue decomposition on the covariance matrix. Calculate the linearity of each point based on the distribution ratio of the three eigenvalues ​​to generate the linear feature parameters of each point. S13: Based on the spatial neighborhood index relationship, extract the high-order sequence of the neighborhood point set of each point, sort the high-order sequence from low to high, and cut off the low end part as the local reference surface candidate point. Perform plane fitting on the local reference surface candidate point to obtain the local reference surface elevation. Subtract the local reference surface elevation from the elevation value of each point to obtain the floating height of each point, and generate the floating feature parameters of each point. S14: Associate and store the linear feature parameters and the suspending feature parameters of each point into the attribute field of each point to generate a DSM expression point set with local geometric features.

3. The method according to claim 1, characterized in that, The geometric feature conditions include linear feature conditions and levitation feature conditions, and S2 includes: S21: Based on the linear feature conditions, the linear feature parameters of each point in the DSM expression point set are judged one by one, and points with linear feature parameters greater than the preset linearity threshold are marked as initial selection points to generate a set of initial linear feature points. The linear feature condition is used to indicate the degree of concentration of point cloud distribution along a certain direction in a local neighborhood. The larger the linear feature parameter, the more significant the extension trend of the point cloud along the corresponding direction. It is used to distinguish between slender linear targets and area-shaped or scattered ground features. S22: Based on the suspension feature conditions, the suspension feature parameters of each point in the initial selection point set of linear features are judged one by one, and points whose suspension feature parameters are less than or equal to the preset suspension height threshold are screened out from the initial selection point set of linear features, thereby generating a candidate point set that satisfies the dual condition constraints. The suspension feature condition is used to indicate the degree of elevation of the point cloud relative to the ground surface or stable ground surface below it. The greater the suspension height, the more likely the point cloud is to be an overhead target. It is used to distinguish between power line remnants and points that are close to the ground or attached to the surface of buildings or vegetation. S23: Based on the spatial neighborhood index relationship, perform neighborhood connectivity analysis on each point in the candidate point set, classify the candidate points that satisfy the preset spatial adjacency relationship into the same connected component, and calculate the spatial extension span along the main direction and the lateral spread width perpendicular to the main direction for each connected component to generate the geometric morphological parameters of each connected component. S24: Based on the geometric morphology parameters, each connected component is screened out, and connected components with a spatial extension span less than a preset length threshold or a lateral spread width greater than a preset width threshold are removed. The candidate points in the remaining connected components are used as suspected residual candidate points of the power line, and a candidate point cluster set is generated.

4. The method according to claim 1, characterized in that, S3 includes: S31: Read the tower location data of the test area, and based on the three-dimensional coordinates of each adjacent tower in the tower location data, take the spatial line between two adjacent towers as the corridor centerline, and expand the preset width range to both sides with the corridor centerline as the center to construct the power line corridor space range corresponding to each span, and generate a set of corridor space constraint areas. S32: Project each member point of the candidate point cluster in the candidate point cluster set to the corridor spatial constraint area set for spatial belonging determination. Calculate the proportion of the number of points of each candidate point cluster falling into any corridor spatial constraint area to the total number of candidate point clusters as the corridor belonging score. Based on the corridor belonging score, select candidate point clusters with scores greater than the preset corridor belonging threshold and generate a candidate point cluster set in the corridor. S33: Based on the direction of the central axis of the corridor to which each candidate point cluster belongs, the member points of each candidate point cluster in the candidate point cluster set within the corridor are transformed from the original three-dimensional coordinate system to a local coordinate system of the corridor with the direction of the corridor as the vertical axis, the direction perpendicular to the corridor as the horizontal axis, and the elevation as the vertical axis, thereby generating a local coordinate point set for each candidate point cluster. S34: In the local coordinate system of the corridor, the local coordinate point set of each candidate point cluster is fitted with an arc sag curve along the vertical axis. After obtaining the initial parameters of the fitted model using the random sampling consensus algorithm, the optimal model parameters are solved by least squares optimization. Based on the optimal model parameters, the root mean square error of the fitting residual and the lateral width of each candidate point cluster are calculated. The candidate point clusters whose root mean square error of the fitting residual is less than a preset residual threshold and whose lateral width is less than a preset width threshold are confirmed as point clusters that conform to the arc sag shape of the power line, and the confirmed power line residual point clusters are generated.

5. The method according to claim 1, characterized in that, The preset anomaly verification conditions include morphological anomaly verification conditions, profile elevation anomaly verification conditions, and neighborhood elevation support verification conditions. S4 includes: S41: Rasterize the DSM representation point set with a preset DSM grid resolution, count the elevation values ​​of all points falling into each grid cell and take the maximum value as the elevation value of the corresponding grid cell to generate a DSM raster. S42: Map the three-dimensional coordinates of each point in the confirmed power line residual point cluster to the grid cell index of the DSM grid, mark the mapped grid cells as abnormal grids to be verified, and generate a set of abnormal grids to be verified. S43: Perform morphological anomaly verification on the set of abnormal grids to be verified: calculate the aspect ratio of the maximum extension length along the corridor direction to the maximum width perpendicular to the corridor direction of the set of abnormal grids to be verified, filter out abnormal grids whose aspect ratio is less than a preset aspect ratio threshold, and generate an abnormal grid set that passes the morphological verification. S44: Perform profile elevation anomaly verification on the abnormal grid set that has passed the morphological verification. Extract the DSM elevation profile sequence at each anomaly location along the vertical corridor direction. Calculate the first-order difference value of the elevation profile sequence near the lowest point of the sag. If the first-order difference value shows a peak shape that changes from positive to negative and the peak width is less than a preset width threshold, it is determined that the profile verification has been passed, and an abnormal grid set that has passed the profile verification is generated. S45: Perform neighborhood elevation support verification on the set of anomalous grids that have passed the profile verification. Determine the point cloud neighborhood corresponding to each anomalous grid according to the spatial neighborhood index relationship. After mapping the point cloud neighborhood to the DSM grid neighborhood, calculate the elevation standard deviation of the neighborhood grids. Calculate the comprehensive anomaly verification index of each anomalous grid. If the comprehensive anomaly verification index of the anomalous grid is greater than the preset comprehensive judgment threshold, it is determined to have passed the neighborhood support verification. Identify the points corresponding to the anomalous grids that have passed the morphological verification, profile elevation anomaly verification, and neighborhood elevation support verification as power line remnant points. Remove the power line remnant points from the DSM expression point set to generate an updated DSM expression point set and void regions.

6. The method according to claim 5, characterized in that, The formula for calculating the comprehensive anomaly verification index is as follows: in, For abnormal grid The comprehensive anomaly verification index is dimensionless and ranges from 0 to 1. , is the weighting coefficient of the morphological verification index, dimensionless, with a value ranging from 0.2 to 0.4, and ; For abnormal grid The aspect ratio is dimensionless; This is a preset maximum aspect ratio reference value, dimensionless, used to... Normalize to the 0 to 1 interval; The weighting coefficients for the profile verification indicators are dimensionless. For abnormal grid The peak value of the first-order difference of the profile elevation, in meters (m). The preset maximum reference value for the profile differential peak value, in meters (m), is used to... Normalize to the 0 to 1 interval; The weighting coefficients for neighborhood support verification indicators are dimensionless. For abnormal grid The standard deviation of the neighborhood grid elevation, in meters (m). The preset maximum reference value for the standard deviation of elevation in the neighborhood, in meters (m), is used to... Normalize to the 0 to 1 interval.

7. The method according to any one of claims 1-6, characterized in that, S5 includes: S51: Based on the updated DSM expression point set and the cavity region, extract the effective grid cells and corresponding elevation values ​​within a preset range around the cavity region as repair reference data, and generate a cavity neighborhood reference dataset. S52: Based on the geometric morphological characteristics of the cavity region, the cavity region is repaired. If the cavity region is distributed as isolated points and its diameter is less than a preset aperture threshold, the cavity region is repaired using inverse distance weighted interpolation. The weighted average is calculated using the reciprocal square of the distance from each effective grid cell in the cavity neighborhood reference dataset to the center of the cavity grid, generating an interpolated repair elevation value. If the cavity region is distributed as a linear extension, effective elevation profiles on both sides are extracted along the cavity extension direction. The cubic spline interpolation method is used to interpolate the cavity grid between the profiles point by point, generating a spline repair elevation value. If the cavity region is distributed as a planar feature and its diameter is greater than or equal to the preset aperture threshold, a quadratic surface equation is constructed using the effective grid cells in the cavity neighborhood reference dataset, and the surface coefficients are solved using the least squares method. The planar coordinates of the cavity grid are substituted into the quadratic surface equation to calculate the surface repair elevation value. S53: Write the interpolated repair elevation value, the spline repair elevation value and the surface repair elevation value of each hole grid into the corresponding position of the DSM raster to generate the repaired DSM raster, and record the original elevation value and the type of repair elevation value of each repair grid to generate a repair area record layer. S54: Perform residual anomaly detection processing on the repaired DSM raster according to the repaired area recording layer, and perform anomaly re-inspection on the repaired area and its neighborhood in the repaired DSM raster. Calculate the proportion of the number of grids with abnormal features to the total number of DSM grids as the anomaly residual rate. If the anomaly residual rate is greater than the preset residual rate threshold, reduce the preset linearity threshold and the preset suspension height threshold, and re-execute S2 to S5. If the anomaly residual rate is less than or equal to the preset residual rate threshold and the gradient change of the elevation transition in the repaired area is less than the preset gradient threshold, the quality is determined to be qualified, and the final DSM result is output. The final DSM result is used to indicate the digital surface model result that meets the quality requirements after power line residue filtering.

8. An automatic power line residue filtering system for airborne point cloud (DSM) production, characterized in that, The system includes: The DSM feature point construction module is used to read the preprocessed and initially classified airborne LiDAR point cloud data and DSM expression category rules, extract the points participating in DSM expression from the airborne LiDAR point cloud data according to the DSM expression category rules, construct the DSM expression point set, perform local neighborhood feature analysis on each point in the DSM expression point set, calculate the geometric feature index of each point, and generate a DSM expression point set with local geometric features. The candidate point cluster extraction module is used to perform feature screening based on the geometric feature indicators of each point in the DSM expression point set, extract points that simultaneously meet the preset geometric feature conditions as candidate points, and perform spatial neighborhood clustering on the candidate points to group the spatially continuously distributed candidate points into the same cluster, generating a candidate point cluster set. The power line cluster identification module is used to read the tower location data of the survey area, construct the power line corridor space range between each tower based on the tower location data, project each candidate point cluster in the candidate point cluster set onto the power line corridor space range for belonging determination, filter out the candidate point clusters falling within the corridor range, and perform spatial morphology fitting on each point cluster in the candidate point clusters falling within the corridor range along the corridor direction, calculate the morphological feature parameters of each point cluster, confirm the point clusters that conform to the power line spatial morphology based on the morphological feature parameters, and generate confirmed power line residual point clusters; The power line point removal module is used to construct a DSM grid based on the DSM expression point set, project the confirmed power line residual point clusters onto the DSM grid, perform anomaly verification on the grid cells at the projection position, confirm the points corresponding to the grid cells that pass the anomaly verification as power line residual points, remove the power line residual points from the DSM expression point set, and generate an updated DSM expression point set and hole regions. The DSM iterative repair module is used to perform surface repair processing on the void regions based on the updated DSM expression point set and the void regions, using neighborhood interpolation or surface fitting to generate a repaired DSM raster. The module then performs a quality evaluation on the repaired DSM raster using residual anomaly detection. Based on the quality evaluation results, the threshold parameters for candidate point selection are updated, and steps S2 to S5 are re-executed until the quality evaluation results meet the preset quality requirements, and the final DSM result is output.

9. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the method of any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method of any one of claims 1 to 7.

Images