A laser radar point cloud single tree segmentation method considering optimal morphology, medium and equipment

Abstract

This invention discloses a method, medium, and device for single-tree segmentation from lidar point clouds that considers optimal morphology. The invention constructs a framework for initial segmentation, segmentation correction, and single-tree morphology optimization. Through a multi-level optimization strategy, it achieves morphologically optimal single-tree segmentation, progressively improving segmentation accuracy and morphological integrity to meet the requirements of high-precision tree parameter extraction. Specifically, during segmentation correction, a single-tree segmentation correction method based on point density histograms detects erroneous segmentation targets by using the vertical point cloud density distribution and detects segmentation points using the horizontal point cloud density distribution, significantly improving the accuracy of single-tree segmentation. During morphology optimization, a graph cut optimization method with high-order energy constraints is introduced to achieve global optimization of single-tree morphology. Based on prior knowledge of single-tree morphology, a high-order energy function integrating geometric features and morphological rules is constructed, and a graph cut optimization algorithm is used to globally optimize the single-tree boundary, effectively solving the problem of crown morphology errors in traditional methods.

Description

Technical Field

[0001] This invention belongs to the field of photogrammetry and remote sensing technology, specifically relating to a method, medium, and equipment for single-tree segmentation of lidar point clouds that takes into account optimal morphology. Background Technology

[0002] As the basic building block of a forest, the structural parameters of individual trees are crucial for forest resource surveys and ecological parameter estimations. Obtaining information on individual trees plays a key role in forest resource monitoring, carbon storage measurement, and biodiversity assessment. Currently, traditional methods for obtaining forest structural parameters mainly rely on manual field measurements or aerial image interpretation. Manual field measurements are time-consuming, labor-intensive, and susceptible to subjective factors. They also struggle to accurately distinguish between low-lying canopies and irregularly shaped trees, failing to meet the needs of continuous, large-scale, accurate, and efficient forest resource surveys.

[0003] Traditional single-tree segmentation methods are based on LiDAR point cloud attributes and rules, and can be mainly divided into the following three categories:

[0004] (1) CHM-based method: By generating a canopy height model (CHM), the tree canopy is segmented using algorithms such as local maximum detection combined with region growing and watershed. Although computationally efficient, it is easily affected by canopy overlap and terrain undulation in high-density forest areas, leading to boundary errors and loss of details. The limitations of CHM-based single-tree segmentation are mainly reflected in: ① In high-density forest scenes, the large overlap area of ​​the tree canopy leads to incorrect segmentation of the tree canopy boundary; ② In areas with large terrain undulations, normalized CHM will cause changes in the tree canopy morphology, while image smoothing strategies are prone to loss of tree canopy details. Therefore, this method is prone to incorrect segmentation, and the tree canopy morphology of the segmentation result is prone to distortion.

[0005] (2) Point-based methods: Directly manipulate point cloud data and use clustering (such as K-means, mean shift) or tree trunk detection strategies. Although complete forest stand information is preserved, the computational redundancy is large, the robustness is poor, and it depends on parameter settings. The limitation of point-based methods is that the amount of data required for calculation is much greater than the amount of data contained in CHM raster data, resulting in low computational efficiency.

[0006] (3) Voxel-based methods: Point clouds are voxelized and segmented by feature classification, which balances computational efficiency and data volume. However, they are sensitive to voxel size and are prone to missing small trees or missegmenting large trees. Although voxelizing point clouds can reduce data redundancy, the data volume is still a challenge compared to raster data.

[0007] While traditional methods possess generalization capabilities, they generally suffer from oversegmentation, undersegmentation, and limitations in scene adaptability. In recent years, machine learning technology, with its advantages of high accuracy, strong generalization, and automation, has become a new direction for overcoming traditional bottlenecks. Machine learning-based single-tree segmentation methods significantly improve accuracy through data-driven strategies and can be mainly divided into the following two categories:

[0008] (1) Traditional machine learning (such as random forest): Tree trunk detection is optimized through feature engineering (such as geometric confidence), achieving an accuracy of 88.5%, but it is difficult to distinguish geometrically similar objects (such as street lamp poles and tree trunks). Machine learning methods rely on the setting of the initial learning threshold and have poor transferability.

[0009] (2) Deep learning methods: Direction vector prediction (e.g., PED-net): Combines semantic segmentation and voxel region growth, adapts to complex scenes, but has weak generalization; 3D sparse convolution: Uses circular detection boxes to match tree structure, with an accuracy of 82.9% (IoU>0.5), and is computationally efficient, but is limited by the lack of airborne point cloud features; Multi-stage network: Integrates centroid prediction, attention mechanism and dynamic clustering, reduces label dependence, but has insufficient module collaboration; General model (e.g., 3D CNN): Supports multi-platform data, has strong cross-scene applicability, but has low accuracy in multi-layer canopy segmentation of broadleaf forests. Summary of the Invention

[0010] The purpose of this invention is to provide a method, medium, and device for single-tree segmentation of lidar point clouds that takes into account optimal morphology, in order to solve the above-mentioned problems.

[0011] This invention is mainly achieved through the following technical solutions:

[0012] A method for single-tree segmentation of lidar point clouds that takes into account optimal morphology includes the following steps:

[0013] Step S1: Initial segmentation; Collect UAV point cloud and ground-based point cloud data of the target area and process them to obtain the canopy height model (CHM); Based on the canopy height model (CHM), extract the initial segmentation clusters;

[0014] Step S2: Segmentation correction;

[0015] Step S21: Extract the number of tree trunks in a single segmentation cluster using the characteristics of tree point clouds in the vertical dimension; detect incorrectly segmented targets by the point cloud density distribution in the vertical direction;

[0016] Step S22: Perform segmentation and correction based on the error type;

[0017] For clusters containing multiple trunks, adhesion segmentation is performed; the connection between the center points of the trunks is detected, and the segmentation points are detected using the horizontal point cloud density distribution; for clusters without trunks, point cloud fusion is performed.

[0018] Step S3: Use a high-order energy-constrained graph cut optimization algorithm to optimize the morphology of a single tree;

[0019] Step S31: Perform an octree region growing algorithm on the point cloud to construct a graph structure of supervoxels and Markov fields;

[0020] Step S32: Calculate the central symmetric point of each supervoxel node with respect to all trunk points, and generate a list of symmetries for multi-level search;

[0021] Step S33: Use the α-expansion algorithm to perform graph cut optimization to obtain the optimal single-tree morphology; the α-expansion algorithm maps the image into a network graph, constructs an energy function with respect to the labels, and then minimizes the energy function; the energy function is:

[0022]

[0023] in Energy function for data items;

[0024] D i (l i ) represents a data item describing the supervoxel v i Assigning label l i The cost;

[0025] V represents a voxel;

[0026] The energy function for the smoothing term;

[0027] S ij (l i , l j ) represents the smoothing term, describing the adjacency of supervoxels v. i and v j The consistency cost of assigning labels;

[0028] l i and l j These are the supervoxels with adjacent relationships, v i and v j Assign tags;

[0029] λ smooth For smoothing parameters;

[0030] e ij For connecting voxels v i and v j The edge;

[0031] E represents energy;

[0032] It is a higher-order energy function;

[0033] These are the weights in the weight list;

[0034] δ(l i , l j () is used to determine label consistency;

[0035] S is the range of symmetric constraints.

[0036] To better implement the present invention, further, in step S33, data item D i (l i )for:

[0037]

[0038] Where: d i For node v i The Euclidean distance to the seed point;

[0039] d t The average radius of the tree canopy in the scene;

[0040] σ1 and σ2 control the attenuation rates at both ends, respectively;

[0041] C is a continuity constant;

[0042] Smoothing term S ij (l i , l j )for:

[0043]

[0044] Where: w ij The edge weight is denoted as .

[0045] To better implement the present invention, step S33 further includes the following steps:

[0046] e1: Initialize labels, assign an initial label to each node, and assign values ​​to the nodes using the initial segmentation results from step S1;

[0047] e2: Construct a graphical model and calculate the initial total energy E using the energy function E(L). current ;

[0048] e3: Select a label α from the label set as the expansion target;

[0049] e4: Construct an extended graph, defining the source node as label α and the sink node as the set of all non-α labels; define data item edges as source edges and sink edges, with the source edge being each node v. i An edge connecting to the source node represents the cost at which that node chooses α as its label, and the sink edge is for each node v.i An edge connected to a sink node represents the cost at which that node chooses a non-α label; edges for smoothing terms retain their original definitions.

[0050] e5: Use the maximum flow algorithm to solve for the maximum flow in the extended graph, which corresponds to the minimum cut of the graph, and divide the graph into the α part connected to the source node and the non-α part;

[0051] e6: Update the label assignment based on the minimum cut result, if v i Connecting the source node, its label l i =α, otherwise retain the original label;

[0052] e7: Calculate the updated energy E new If E new <E current If so, the label assignment from this expansion will be retained, allowing E to... new =E current ;

[0053] e8: Repeat steps e3 to e7, perform α-expansion for each label, and gradually approach the global optimum. If the energy no longer decreases or the maximum number of iterations is reached after several iterations, terminate the algorithm and output the final label assignment.

[0054] To better implement the present invention, step S1 further includes the following steps:

[0055] Step S11: Data acquisition. The forest point cloud and ground points are rasterized to obtain the digital surface model (DSM) and digital elevation model (DEM) respectively. The elevations within the corresponding raster are subtracted to obtain the canopy height model (CHM).

[0056] Step S12: Extract local extrema of the canopy height model CHM as seed points for region growth;

[0057] Traverse the CHM model and mark candidate seed points; then, calculate the canopy convexity and retain points with canopy convexity greater than a set threshold; finally, perform spatial deduplication and edge removal in sequence to filter and obtain seed points.

[0058] Step S13: Perform region growing to extract preliminary segmented clusters.

[0059] To better implement the present invention, step S21 further includes the following steps:

[0060] a1: High-level slice construction, vertical slices are made on the initially segmented single-tree point cloud, and the number of points D(z) in each layer is calculated;

[0061] a2: The standardized difference method is used to identify the elevation of the trunk-crown boundary point. The elevation density distribution sequence is given as:

[0062] D = {D(z1), D(z2), ..., D(z...}} n )},z i =z min +iΔz;

[0063] Where: D(z) n ) represents the number of points in the nth layer;

[0064] z min The lowest height considered;

[0065] With the slice step size Δz set to 0.1m, calculate the first-order forward difference:

[0066] ΔD k =D(z) k+1 )-D(z k k = 1, 2, ..., n-1;

[0067] Based on the differential sequence reflecting the density change rate between adjacent elevation layers, a positive transition occurs at the boundary between the trunk and canopy; thus, the trunk layer height z can be located. trunk Extract the point cloud subset P of the trunk layer. trunk ={p i |z i ≤z trunk};

[0068] a3: Preliminary extraction of tree trunks using Euclidean clustering;

[0069] The Euclidean distance between points is calculated, and a neighborhood search is performed for each point to cluster all points. The clustering results are then filtered based on height constraints and principal component analysis constraints, which include principal direction constraints and linear feature constraints.

[0070] To better implement the present invention, step S22 further includes the following steps:

[0071] b1: Segment clusters containing multiple trunks by connecting them together; detect the connection between the center points of the trunks and find the segmentation point;

[0072] First, rotate the point cloud around the Z-axis so that the axis is parallel to the X-axis. Divide the data into intervals and count the number of points falling within each interval to obtain an under-divided histogram containing two maxima. Then, use the minima falling within the interval of the two maxima as the cut points.

[0073] b2: Perform point cloud fusion on clusters without tree trunks;

[0074] For over-segmented single trees without trunks i Find all adjacent single tree Ts that contain a single tree trunk. i ,

[0075] Calculate the T of a single tree i With single wood t i The Euclidean distance, and based on the Euclidean distance, the single-wood T i An optimized queue is formed, and then fusion is performed based on the projection overlap ratio constraint.

[0076] To better realize the present invention, step S31 further includes the following steps:

[0077] c1: Construct a point cloud octree and use the similarity h of the points to describe the characteristics of the octree nodes;

[0078] h = h Coord +h Normal +h Intensity ;

[0079] Where: h coord Represents the Euclidean distance calculated based on the three-dimensional coordinates of the voxel center point;

[0080] h Normal This represents the angle between the normal vectors of two voxels;

[0081] h Intensity The difference between the average intensity values ​​of two points within a voxel;

[0082] c2: Performs regional growth to generate supravoxels;

[0083] c3: Construct the Markov field diagram structure;

[0084] Each hypervoxel is treated as a node and assigned a unique label. i ∈{1, 2, ..., N}, where N is the number of tree trunk seed points; the initial value of this label is the segmentation result corrected in step S2; if the spatial distance between two hypervoxels is less than 2m, then an edge connection is established, and the weight is determined by the feature difference degree.

[0085] To better realize the present invention, step S32 further includes the following steps:

[0086] d1: Retrieve a supervoxel node v from the list of supervoxel nodes. i ;

[0087] d2: Compute node v i Regarding seed point P j Central symmetric position

[0088]

[0089] Among them (P) x P y ) and (v x vy v z () represents the three-dimensional coordinates of the seed point and the node, respectively;

[0090] d t This is the prior threshold;

[0091] dis(v i P j (P is the seed point) j and node v i The plane Euclidean distance;

[0092] d3: in symmetrical position Centered on a node v, a kd-tree neighborhood search is performed on the supervoxels within three increasing radii to construct node v. i Regarding seed point P j symmetric point set

[0093] d4: Record node v after the search is complete. i Regarding seed point P j The set of all symmetrical points and its corresponding weight list

[0094] d5: Repeat steps d1 to d4 until all supervoxel nodes generate their symmetric lists.

[0095] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the above-described method for morphologically optimal single-tree segmentation of lidar point clouds.

[0096] An electronic device includes a memory and a processor; the memory stores a computer program; the processor is configured to execute the computer program in the memory to implement the above-described method for morphologically optimal single-tree segmentation of lidar point clouds.

[0097] The beneficial effects of this invention are as follows:

[0098] (1) This invention uses elevation density features to extract individual tree trunks to identify incorrect segmentation, and introduces a graph cut optimization method with higher-order energy constraints to achieve accurate and morphologically complete individual tree segmentation. Addressing the problems of low segmentation accuracy and poor individual tree morphology in existing individual tree segmentation methods, this invention establishes a "preliminary segmentation - segmentation correction - morphological optimization" individual tree segmentation process framework. This invention achieves optimal morphological segmentation of individual trees through a multi-level optimization strategy, gradually improving the accuracy and morphological integrity of the segmentation, and meeting the requirements for high-precision tree parameter extraction.

[0099] (2) To address the problem of insufficient initial segmentation accuracy, and considering the characteristics of the ground point cloud, a single-tree segmentation correction method based on point density histograms is proposed. This method detects erroneous segmentation targets by using the vertical point cloud density distribution and identifies segmentation points using the horizontal point cloud density distribution, significantly improving the accuracy of single-tree segmentation. To address the difficulty of traditional segmentation methods in considering the morphology of individual trees, a graph cut optimization method with higher-order energy constraints is introduced to achieve global optimization of the single-tree morphology. Based on prior knowledge of the single-tree morphology, a higher-order energy function integrating geometric features and morphological rules is constructed. The graph cut optimization algorithm is then used to globally optimize the single-tree boundary, effectively solving the problem of crown morphology errors in traditional methods.

[0100] (3) This invention can effectively optimize the segmentation accuracy and crown morphology integrity of the traditional CHM region growth segmentation method. This invention relies on two key components: (1) segmentation correction guided by point cloud elevation density distribution characteristics and (2) a crown morphology map-cutting optimization method with higher-order energy constraints. Unlike the traditional single-tree segmentation method based on CHM grids, this method is more stable in complex forest environments, has higher single-tree segmentation accuracy, and effectively solves the problem of crown morphology distortion in traditional segmentation methods. Attached Figure Description

[0101] Figure 1 The flowchart of the LiDAR point cloud single-tree segmentation method that takes into account the optimal morphology of the present invention is shown below.

[0102] Figure 2 This is the vertical point density distribution histogram generated in Example 1;

[0103] Figure 3 This is the histogram of undersegmentation in Example 1;

[0104] Figure 4 Example 1: Single-tree segmentation before and after optimization;

[0105] Figure 5 Example 2 of single-wood segmentation before and after optimization;

[0106] Figure 6 Example 3 shows the single-wood segmentation before and after optimization;

[0107] Figure 7 Example 4 shows the single-wood segmentation before and after optimization. Detailed Implementation

[0108] Example 1:

[0109] A method for single-tree segmentation of lidar point clouds that takes into account optimal morphology, such as Figure 1 As shown, the specific steps include:

[0110] Step S1: Data acquisition and initial segmentation.

[0111] The drone point cloud and ground-based point cloud of the target area are collected, and the corresponding point cloud of the upper canopy and the tree trunk and lower vegetation are obtained. The two are then registered and fused, and the digital elevation model (DEM) of the area is extracted as the ground points.

[0112] (1) Data acquisition and generation of the canopy height model (CHM);

[0113] Data acquisition: UAV LiDAR: flight altitude 50m, scanning angle 70°, point density ≥50 points / ㎡, to acquire point cloud of the upper canopy; Ground-based LiDAR: scanner erection height 1.5m, 360° full-station scanning, point density ≥200 points / ㎡, to acquire point cloud of tree trunk and lower vegetation.

[0114] UAV and ground-based point cloud registration; manual registration: feature matching: select at least 3 sets of corresponding feature points from the two types of point clouds (such as significant tree branch points, rock vertices); ICP fine registration: use the Iterative Closest Point (ICP) algorithm, set the maximum number of iterations to 100, and the registration error threshold ≤ 0.1m; fusion verification: calculate the overlap of point clouds in the overlapping area, requiring the root mean square error (RMSE) < 0.15m.

[0115] Ground point extraction and DEM generation; using progressive triangulation densification filtering (PTD), setting a slope threshold of 8°, and a window size of 1m×1m; generating a digital elevation model (DEM) with a raster resolution of 0.25m, and using inverse distance weighting (IDW) as the interpolation method.

[0116] b. DSM construction: Rasterize the non-ground point cloud (0.25m resolution) and take the maximum elevation value of the point cloud within each grid cell;

[0117] c. Construction of the Canopy Height Model (CHM):

[0118] CHM(i,j)=DSM(i,j)-DEM(i,j);

[0119] CHM, DSM, and DEM are all raster images. The pixel value represents the point cloud elevation at the current pixel location, and (i,j) represents the image coordinates. Negative values ​​are set to 0. Gaussian filtering: CHM is smoothed using a 3×3 convolution kernel (σ=0.5) to preserve the main peak features of the canopy.

[0120] (2) Extraction of seed points;

[0121] Using a 3×3 sliding window to traverse the CHM, the condition for marking the center point of the window as a candidate seed point is as follows:

[0122] CHM(x,y)>CHM(x+i,y+j)(i,j∈{-1,0,1}, (i,j)≠(0,0));

[0123] Seed point selection: Canopy elevation calculation:

[0124] Where N is the number of grids within a radius of 1 meter;

[0125] Retain points where C > 0.7 and remove spurious extrema from flat regions. Spatial deduplication: If the distance between multiple extrema points is less than 1 meter, retain only the highest point. Edge removal: Remove points less than 2 meters from the CHM boundary to avoid edge effects.

[0126] (3) Perform region growing to extract preliminary cluster segments;

[0127] ①Use the local maximum point as the starting seed point for the first batch of growth. Define the eight pixels around the seed point for each growth as the neighboring pixels of the seed point. Each growth only performs condition judgment on the neighboring pixels of the current seed point. The neighboring pixels that meet the conditions are merged into the region of the seed point and used as the seed point for the next cycle. Only when all eight neighboring pixels of the current seed point have been conditionally judged will it enter the next cycle.

[0128] ② Calculate the height difference between the point to be grown and the current crown vertex. Only points whose height is less than the height of the current seed point and greater than a set height threshold can be merged into the region of the current seed point. Generally, the tree vertex is the highest point of a tree, and the height of crown points should be less than the tree vertex. Under natural growth conditions, points at the edge of the crown will not fall near the ground but will be greater than a certain threshold. This invention sets this threshold as H. max *0.55,

[0129] H max This represents the height threshold of the current region's tree canopy vertices. Therefore, the height relationship between the point to be grown and the current region's tree canopy vertices should satisfy:

[0130] H max *0.55<h i <H max In the formula, h i This indicates the height of the growth point.

[0131] ③ Calculate the height difference between the point to be grown and the current seed point, and set the adjacent height difference threshold to 0.8. Therefore, the relationship between the height of the point to be grown and the height of the current seed point should satisfy the following formula, where h i h represents the height of the growth point. j This represents the height of the current seed point.

[0132] h j *0.8<h i <h j ;

[0133] ④ Determine the Euclidean distance between the growth point and the crown vertices within the current region, using the externally input threshold parameter d. t To control its growth range.

[0134] ⑤ Adjacent Pixel Competition Rule. If overlapping areas occur during region growth, a reasonable competition rule is set to determine the ownership of the point to be grown. The rule used is to calculate the planar Euclidean distances D1 and D2 from the point to be grown to the two competing tree crown vertices, and calculate the height differences H1 and H2 from the point to the two tree crown vertices. The ratio of these two parameters is calculated. If the ratio of the height difference H1 from the point to be grown to the tree crown vertex of Tree1 to the distance D1 is less than the ratio of the height difference H2 from the point to the other tree crown vertex to the distance D2, then the point is assigned to Tree1; otherwise, it is assigned to the other seed point.

[0135] Step S2: Cluster segmentation correction; analyze the elevation density distribution characteristics of individual clusters, determine the type of misidentification, and correct it.

[0136] (1) Since tree point clouds exhibit significant hierarchical characteristics in the vertical dimension, this characteristic is used to extract the number of tree trunks in a single segmentation cluster.

[0137] a. High-level slicing construction: Vertical slices are made on the initially segmented single-tree point cloud with a step size of 0.1m, and the number of points in each layer is calculated:

[0138] D(z)=∑δ(z i ∈[z i , z i +Δz));

[0139] Where: z i This is the height value, i.e., the z-axis coordinate value;

[0140] δ is a function of the number of points within the calculation range;

[0141] z represents the height value, the ordinate;

[0142] Δz is the slice step size; Δz = 0.1m;

[0143] like Figure 2 As shown, a histogram of point density distribution in the vertical direction is generated.

[0144] b. Mutation threshold detection: The standardized difference method is used to identify the elevation of the trunk-crown boundary point. Its core lies in identifying discontinuous transitions in point cloud density in the vertical direction. Let the elevation density distribution sequence be:

[0145] D = {D(z1), D(z2), ..., D(z...}} n )},z i =z min +iΔz;

[0146] Where: D(z) n ) represents the number of points in the nth layer;

[0147] z min The lowest height considered;

[0148] The slice step size Δz is set to 0.1m. The first-order forward difference is calculated:

[0149] ΔD k =D(z) k+1 )-D(z k k = 1, 2, ..., n-1;

[0150] This differential sequence reflects the rate of density change between adjacent elevation layers, and a significant positive transition will occur at the boundary between the trunk and the canopy, such as... Figure 2 The location indicated by the red arrow in the middle.

[0151] The above steps can be used to accurately determine the height z of the tree trunk layer. trunk Extract the point cloud subset P of the trunk layer. trunk ={p i |z i ≤z trunk The next step is to initially extract the tree trunk using Euclidean clustering. Euclidean clustering is a point cloud data clustering method based on the Euclidean distance metric. Its core idea is to group points in the point cloud that are close in Euclidean distance into the same class. For two points p in the point cloud space... i (x i y i , z i ) and p j (x j y j , z j Its Euclidean distance is defined as:

[0152]

[0153] This distance metric reflects the actual geometric distance between points in Euclidean space. Then, a neighborhood search is performed for each point, given a point p. i And its neighborhood radius r, its neighborhood is defined as:

[0154] N(p i ,r)={p j ∈P|d(p i p j )≤r};

[0155] Where: P represents the point cloud data to be clustered. Neighborhood search calculates the value of point p. iThe Euclidean distance to all other points is used to filter out points whose distance is less than the neighborhood radius r. The neighborhood radius r used in this invention is 0.2m. By performing a neighborhood search on all points, all points can be marked as a certain cluster. The minimum distance between all points within a cluster is less than the neighborhood radius.

[0156] To avoid the influence of ground-level shrubs and weeds on the accuracy of tree trunk extraction, the clustering results need to be further filtered to remove shrubs and weeds that are misidentified as tree trunks. The filtering criteria used here are as follows:

[0157] ① Height constraint, the goal is to exclude low vegetation and non-tree features. For each cluster C m Calculate its elevation span:

[0158] H m =max(z) m )-min(z m );

[0159] Set height threshold H t If H m ≥H t If the candidate tree trunk is selected, the height threshold H used in this invention will be retained. t =1.2m.

[0160] ② Principal component analysis constraints, which consist of principal direction constraints and linear characteristic constraints, aim to preserve the trunk structure that is perpendicular to the ground and approximately straight. Only when a cluster satisfies both the principal direction constraints and the linear characteristic constraints is it considered to satisfy the principal component constraints.

[0161] For clustering C m Perform principal component analysis to obtain its first principal component vector v1 = (a1, b1, c1). Define the angle between the principal direction and the z-axis:

[0162]

[0163] Where: a 1， b1 and c1 are the three directional components.

[0164] When C m When the angle is ≤15°, the cluster is considered to be approximately perpendicular to the ground.

[0165] Linearity can be calculated using the three eigenvalues ​​λ1, λ2, and λ3 obtained from principal component analysis, and the linearity exponent is defined as: L m = (λ1-λ2)÷λ3; it is stipulated that when L m When the value is ≥0.7, the cluster is considered to be approximately a straight line.

[0166] (2) Perform segmentation correction based on error type; segment clusters containing multiple trunks and merge clusters without trunks.

[0167] a. Multi-trunk adhesion segmentation. For samples containing multiple trunks after initial segmentation, this is considered undersegmentation, and optimization is performed to address this error. This requires detecting the lines connecting the center points of the trunks and finding suitable segmentation points. To facilitate histogram calculation, the point cloud needs to be rotated around the Z-axis so that the axis is parallel to the X-axis.

[0168] For example, assuming the centers of the tree trunk are c1(x1,y1) and c2(x2,y2), calculate the slopes of c1 and c2:

[0169]

[0170] The rotation angle α can be obtained from the slope: α = arctan(k);

[0171] The rotation transformation is as follows:

[0172]

[0173] Where: R is the rotation matrix;

[0174] p′ i For point p i Points after rotation transformation;

[0175] like Figure 3 As shown, after rotation, the points are divided into intervals along the X-axis with a step size Δs = 0.1m. The number of points falling within each interval is counted to obtain an under-segmented histogram, which should contain two maxima (tree trunk points). Minimum points will appear between the two maxima. The minimum point falling within the interval formed by the two maxima is used as the cutting point, and the dividing line is perpendicular to the line connecting the tree trunks. After segmentation, the point cloud is inversely transformed back to the original coordinate system.

[0176] b. Merging of Tree-less Point Clouds. After undersegmentation optimization, the next step is to optimize oversegmented trees that do not contain trunks. Find all trees adjacent to tree t that contain a trunk and form a set T. Calculate the Euclidean distance between the geometric centers of all trees in T and the geometric center of tree t, and sort them according to the distance to enter the optimization queue. Subsequently, determine whether to merge them. Theoretically, the goal is to find an optimal tree T for tree t within set T. i The relationship between the two is that t belongs to T. i Find T i The iteration terminates and exits. To ensure the similarity of tree morphology after merging, a projection overlap ratio is defined to constrain the merging process. It is stipulated that the ratio of the major axis to the minor axis of the merged single tree should not exceed 2, and the distance between the geometric center of t and the nearest trunk center should not exceed twice the average crown radius of the scene.

[0177] Step S3: Perform single-tree morphology optimization using a high-order energy-constrained graph cut optimization algorithm.

[0178] (1) Use the octree region growing algorithm on the point cloud to construct the graph structure of supervoxels and Markov fields.

[0179] a. Construct an octree for the point cloud. First, randomly select several seed points P from the leaf nodes of the octree. seed To ensure the seed points are distributed as evenly as possible in space, the characteristics of octree nodes are described using point similarity. Point similarity h measures spatial distance, intensity value, and normal vector:

[0180] h = h Coord +h Normal +h Intensity ;

[0181] Where: h coord h represents the Euclidean distance calculated based on the three-dimensional coordinates of the voxel center. Normal h represents the angle between the normal vectors of two voxels. Intensity This is the difference in average intensity values ​​between two voxel points.

[0182] After obtaining the three similarity features, in order to balance the influence of the three features on the similarity, they need to be normalized:

[0183]

[0184] Where (x1, y1, z1) and (x2, y2, z2) represent the center coordinates of two voxels, r is the resolution of the octree, θ(n1, n2) represents the angle between the normal vectors of the two voxels, and I1 and I2 represent the mean values ​​of the point intensities of the two voxels. max This represents the maximum intensity value of the point cloud data.

[0185] b. Perform region growth to generate supravoxels. The steps of region growth are as follows:

[0186] ① Put all seed points P seed Join the growth queue;

[0187] ② Each time, a seed point P is taken from the growth queue. i As the current seed point, use octree growth to search all voxels in its neighborhood to generate the search queue for the seed point.

[0188] ③ Traverse the neighborhood voxels p in its search queue i If the search queue is empty, it is assumed that the supervoxel has reached the maximum voxel radius, and it is removed from the growth queue and the process returns to step ②.

[0189] ④ Calculate the seed point P i Its neighboring voxel p i If the similarity h is less than a threshold, then the voxel p is labeled. i Belongs to seed point P i The supervoxel, and added with p i In the search sequence from adjacent voxels to the seed point, it is required that the adjacent voxel points are related to the seed point P. i The distance is less than the search radius R of the hypervoxel;

[0190] ⑤ Repeat steps ② to ④ until all voxels are assigned to supervoxels or the queue is empty.

[0191] c. Construct the Markov field diagram structure.

[0192] Graph nodes: Each hypervoxel is treated as a node and assigned a unique label. i ∈{1, 2, ..., N} (N is the number of seed points in the trunk), the initial value of this label is the corrected segmentation result obtained in the second step. Edge linking: if the spatial distance between two supervoxels is less than 2m, an edge link is established, and the weight is determined by the feature difference.

[0193] (2) Calculate the central symmetry point of each supervoxel node with respect to all trunk nodes, and generate a list of symmetries for multi-level search. The steps for generating the node symmetry table are as follows:

[0194] ①Retrieve a supervoxel node v from the list of supervoxel nodes i .

[0195] ② Calculate the node v i Regarding a certain seed point P j Central symmetric position Typically, no node will fall on a precisely symmetrical position, therefore It is a virtual location node, if node v i With seed point P j If the planar Euclidean distance exceeds a certain threshold range, then the node with respect to P... j The symmetrical point is denoted as 0. The calculation method is as follows:

[0196]

[0197] Among them (P) x P y ) and (v x v y v z ) represent the three-dimensional coordinates of the seed point and the node, respectively, and d t As a priori threshold, it is usually set to the average radius of the tree canopy in the scene, dis(v iP j (P is the seed point) j and node v i The plane Euclidean distance.

[0198] ③ In symmetrical position Centered on a node v, a kd-tree neighborhood search is performed on the supervoxels within three increasing radii to construct node v. i Regarding seed point P j symmetric point set The weights of the symmetrical points within the three layers of increasing radii are different, and the design is as follows:

[0199] The first layer is a strictly symmetrical region with a search radius of r1 = 0.1d. t If the region contains nodes, and the node labels are the same as v i For different labels, a weight w1=3 is taken, resulting in the most severe energy penalty;

[0200] The second layer can tolerate slight deformation or data noise; the search radius for this region is r2 = 0.3d. t If the region contains nodes, and the node labels are the same as v i The labels are different, so a weight w2=2 is applied, resulting in a moderate energy penalty.

[0201] The third layer is an asymmetric region used when neither the first nor the second layer has any nodes. It can tolerate some missing data or incomplete canopy shapes. The search radius for this region is set to r3 = 0.5d. t If the region contains nodes, and the node labels are the same as v i With different labels, a weight w3=1 is taken, resulting in the minimum energy penalty.

[0202] ④ Record the node v after the search is complete. i For each seed point P j The set of all symmetrical points and its corresponding weight list

[0203] ⑤ Repeat steps ① to ④ until all supervoxel nodes generate their symmetric lists.

[0204] (3) A graph cut optimization algorithm with symmetric constraints is introduced. The α-expansion algorithm is used for graph cut optimization to obtain the optimal single-tree shape. The energy function is shown below:

[0205]

[0206] in Let the energy function be the data item. The energy function for the smoothing term. It is a higher-order energy function.

[0207] a. Data item constraints. Data item D i (l i Description of supervoxel v i Assigning label l i The lower the cost, the better it usually indicates with respect to the label l. i The higher the consistency, the better. For this item, a dynamic distance decay model is proposed, whose functional expression is:

[0208]

[0209] Where, d i For node v i The planar Euclidean distance to the seed point, d t Let σ1 and σ2 be the average crown radius in the scene, and let σ1 and σ2 control the decay rate at both ends, respectively. In this invention, these two parameters are set to σ1 = 1.5 and σ2 = 2, which can control the Euclidean distance from the node to the seed point to increase slowly when it is less than 1.5 times the average crown diameter and increase rapidly when it is greater than 1.5 times the crown diameter.

[0210] C is a continuity constant, the purpose of which is to control the piecewise function in d. i =σ1d t The continuity at a given point is calculated using the following formula:

[0211]

[0212] b. Smoothing term constraint. Smoothing term S ij (l i , l j Describes adjacency relationships of supervoxels v. i and v j The consistency cost of label assignment, calculated using the traditional Potts model, is as follows:

[0213]

[0214] w ij Here are the edge weights, and the method for calculating edge weights is as follows:

[0215]

[0216] Where d t is a priori parameter, representing the average canopy radius of trees in the scene.

[0217] c. Symmetry constraints. The energy design for symmetry is as follows:

[0218]

[0219] Weight The calculation is obtained from a symmetric table, where δ(l)i , l j Used to determine label consistency, l i and l j v i and v j The label, the label is consistent δ(l) i , l j The value of ) is 0, and the value of ) is 1 otherwise. To avoid repeatedly calculating the symmetry energy for a node, in one iteration, for a given node, the symmetry energy penalty is only calculated in the set of symmetry points of the seed point with the same label, and the maximum energy w is taken from the multi-level symmetry energy. max This serves as an energy penalty for the current iteration.

[0220] d. Algorithm flow for α-expansion graph cut optimization. The specific algorithm is as follows:

[0221] ① Initialize labels: Assign an initial label to each node. In this invention, the initial segmentation result obtained above is used to assign values ​​to the nodes, which can effectively reduce the number of algorithm iterations and save algorithm running time.

[0222] ② Construct a graph model. Using the energy function model described above, calculate the edge weights (smoothing terms) of each node and its neighboring nodes, the energy penalty assigned to a node for a certain label (data term), and the higher-order energy penalties under symmetry constraints to obtain the initial total energy E. current .

[0223] ③ Select a label α from the label set as the expansion target.

[0224] ④ Transform the multi-label problem into a binary classification problem, namely, α class and non-α class. Construct an extended graph, defining the source node as the label α and the sink node as the set of all non-α labels; define data item edges as source edges and sink edges, with the source edge being each node v. i An edge connecting to the source node represents the cost at which that node chooses α as its label, and the sink edge is for each node v. i An edge connected to a sink node represents the cost of that node choosing a non-α label; smoothing edges retain their original definition.

[0225] ⑤ Use the maximum flow algorithm to solve for the maximum flow in the extended graph, which corresponds to the minimum cut of the graph. This cut divides the graph into two parts: one part is connected to the source node α, and the other part is not α.

[0226] ⑥ Update the label assignment based on the minimum cut result, if v i Connecting the source node, its label l i =α, otherwise retain the original label.

[0227] ⑦ Calculate the updated energy Enew If E new <E current If so, the label assignment from this expansion will be retained, allowing E to... new =E current .

[0228] ⑧ Repeat steps ③ to ⑦ to perform α-expansion for each label, gradually approaching the global optimum. If the energy no longer decreases or the maximum number of iterations is reached after several iterations, terminate the algorithm and output the final label assignment.

[0229] like Figure 4 As shown, (a) is a schematic diagram of tile segmentation before optimization, and (b) is a schematic diagram of tile segmentation after optimization. The segmentation effect before optimization is distorted, and the various shapes are deformed, while the segmentation effect after optimization is more realistic and effective.

[0230] like Figures 5-7 As shown, (a) is a schematic diagram of single-tree segmentation before optimization, and (b) is a schematic diagram of single-tree segmentation after optimization. Specifically, as... Figure 5 As shown, in the red box of Figure (a), the tree crown is divided in a "one-cut" manner, without considering the growth patterns between the crowns. In the red box of Figure (b), it is clear that the branch and leaf morphology of the tree crown is more in line with objective reality, with branches and leaves of different trees intersecting each other. Figure 6 As shown, in the red box of the unoptimized figure (a), the tree crown is divided in a "one-cut" manner, without considering the growth patterns between the crowns. In the red box of the optimized figure (b), it is clear that the branch and leaf morphology of the tree crown is more in line with objective reality, with branches and leaves of different trees intersecting each other. Figure 7 As shown in the red box of Figure (a) before optimization, the tree crown shape is obviously unreasonable. The tree crown shape should roughly present a symmetrical structure, as shown in the red box of Figure (b) after optimization.

[0231] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Any simple modifications or equivalent changes made to the above embodiments based on the technical essence of the present invention shall fall within the protection scope of the present invention.

Claims

1. A method for single-tree segmentation of lidar point clouds that takes into account optimal morphology, characterized in that, Includes the following steps: Step S1: Initial segmentation; Collect UAV point cloud and ground-based point cloud data of the target area and process them to obtain the canopy height model (CHM); Based on the canopy height model (CHM), extract the initial segmentation clusters; Step S2: Segmentation correction; Step S21: Extract the number of tree trunks in a single segmentation cluster using the characteristics of tree point clouds in the vertical dimension; detect incorrectly segmented targets by analyzing the point cloud density distribution in the vertical direction. Step S22: Perform segmentation and correction based on the error type; For clusters containing multiple trunks, adhesion segmentation is performed; the connection between the center points of the trunks is detected, and the segmentation points are detected using the point cloud density distribution in the horizontal direction; for clusters without trunks, point cloud fusion is performed. Step S3: Use a high-order energy-constrained graph cut optimization algorithm to optimize the morphology of a single tree; Step S31: Perform an octree region growing algorithm on the point cloud to construct a graph structure of supervoxels and Markov fields; Step S32: Calculate the central symmetric point of each supervoxel node with respect to all trunk points, and generate a list of symmetries for multi-level search; Step S33: Use the α-expansion algorithm to perform graph cut optimization to obtain the optimal single-tree morphology; the α-expansion algorithm maps the image into a network graph, constructs an energy function with respect to the label, and then minimizes the energy function; the energy function is: ； in: Energy function for data items; D i (l i ) is a data item, describing a super voxel v i assigning a label l i the cost; V represents a voxel; The energy function for the smoothing term; For smoothing terms, describe the adjacency of supervoxels v. i and v j The consistency cost of assigning labels; l i and l j are respectively adjacent relationship of super voxel v i and v j assigning labels; λ smooth is a smoothing parameter; e ij For connecting voxels v i and v j edges; E represents energy; It is a higher-order energy function; These are the weights in the weight list; Used to determine label consistency; S is the range of symmetric constraints; Step S33 includes the following steps: e1: Initialize labels, assign an initial label to each node, and assign values ​​to the nodes using the initial segmentation results from step S1; e2: construct a graph model, calculate the initial total energy E using the energy function E(L) current ; e3: Select a label α from the label set as the expansion target; e4: Construct an extended graph, defining the source node as label α and the sink node as the set of all non-α labels; define data item edges as source edges and sink edges, with the source edge being each node v. i An edge connecting to the source node represents the cost at which that node chooses α as its label, and the sink edge is for each node v. i An edge connected to a sink node represents the cost at which that node chooses a non-α label; edges for smoothing terms retain their original definitions. e5: Use the maximum flow algorithm to solve for the maximum flow in the extended graph, which corresponds to the minimum cut of the graph, and divide the graph into the α part connected to the source node and the non-α part; e6: Update the label assignment based on the minimum cut result, if v i Connecting the source node, its label l i =α, otherwise retain the original tag; e7: Calculate the updated energy E new If E new <E current If so, the label assignment from this expansion will be retained, allowing E to... new =E current ; e8: Repeat steps e3 to e7, perform α-expansion for each label, and gradually approach the global optimum. If the energy no longer decreases or the maximum number of iterations is reached after several iterations, terminate the algorithm and output the final label assignment.

2. The method for single-tree segmentation of lidar point clouds considering optimal morphology as described in claim 1, characterized in that, In step S33, data item D i (l) i )for: ； Where: d i For node v i The Euclidean distance to the seed point; d t The average radius of the tree canopy in the scene; σ1 and σ2 control the attenuation rates at both ends, respectively; C is a continuity constant; Smoothing terms for: ； Where: w ij The edge weight is denoted as .

3. The method for single-tree segmentation of lidar point clouds considering optimal morphology as described in claim 1, characterized in that, Step S1 includes the following steps: Step S11: Data acquisition. The forest point cloud and ground points are rasterized to obtain the digital surface model (DSM) and digital elevation model (DEM) respectively. The elevations within the corresponding raster are subtracted to obtain the canopy height model (CHM). Step S12: Extract local extrema of the canopy height model CHM as seed points for region growth; Traverse the CHM model and mark candidate seed points; then, calculate the canopy convexity and retain points with canopy convexity greater than a set threshold; finally, perform spatial deduplication and edge removal in sequence to filter and obtain seed points. Step S13: Perform region growing to extract preliminary segmented clusters.

4. The method for single-tree segmentation of lidar point clouds considering optimal morphology as described in claim 3, characterized in that, Step S21 includes the following steps: a1: High-level slice construction, vertical slices are made on the initially segmented single-tree point cloud, and the number of points D(z) in each layer is calculated. a2: The standardized difference method is used to identify the elevation of the trunk-crown boundary point. The elevation density distribution sequence is given as: ， ； Where: D(z) n () represents the number of points in the nth layer; z min This is the starting point of the tree trunk, at its lowest point; With the slice step size Δz set to 0.1m, calculate the first-order forward difference: ； Based on the differential sequence reflecting the density change rate between adjacent elevation layers, a positive transition occurs at the boundary between the trunk and the canopy; thus, the height of the trunk layer can be located. Extracting the point cloud subset of the tree trunk layer ; a3: Preliminary extraction of tree trunks using Euclidean clustering; The Euclidean distance between points is calculated, and a neighborhood search is performed for each point to cluster all points. The clustering results are then filtered based on height constraints and principal component analysis constraints, which include principal direction constraints and linear feature constraints.

5. The method for single-tree segmentation of lidar point clouds considering optimal morphology as described in claim 4, characterized in that, Step S22 includes the following steps: b1: Segment clusters containing multiple trunks by connecting them together; detect the connection between the center points of the trunks and find the segmentation point; First, rotate the point cloud around the Z-axis so that the axis is parallel to the X-axis. Divide the data into intervals and count the number of points falling within each interval to obtain an under-divided histogram containing two maxima. Then, use the minima falling within the interval of the two maxima as the cut points. b2: Perform point cloud fusion on clusters without tree trunks; For over-segmented single trees without trunks i Find all adjacent single tree Ts that contain a single tree trunk. i , Calculate the T of a single tree i With single wood t i The Euclidean distance, and based on the Euclidean distance, the single-wood T i An optimized queue is formed, and then fusion is performed based on the projection overlap ratio constraint.

6. The method for single-tree segmentation of lidar point clouds considering optimal morphology as described in claim 1, characterized in that, Step S31 includes the following steps: c1: Construct a point cloud octree and use the similarity h of the points to describe the characteristics of the octree nodes; ； Where: h coord Represents the Euclidean distance calculated based on the three-dimensional coordinates of the voxel center point; h Normal This represents the angle between the normal vectors of two voxels; h Intensity The difference between the average intensity values ​​of two points within a voxel; c2: Performs regional growth to generate supravoxels; c3: Construct the Markov field diagram structure; Each hypervoxel is treated as a node and assigned a unique label. N is the number of tree trunk seed points; the initial value of the label is the segmentation result after correction in step S2; if the spatial distance between two supervoxels is less than 2m, an edge connection is established, and the weight is determined by the feature difference degree.

7. A method for single-tree segmentation of lidar point clouds considering optimal morphology as described in claim 6, characterized in that, Step S32 includes the following steps: d1: Retrieve a supervoxel node v from the list of supervoxel nodes. i ; d2: Compute node v i Regarding seed point P j Central symmetric position : ； Among them (P) x P y ) and (v x v y v z () represents the three-dimensional coordinates of the seed point and the node, respectively; d t This is the prior threshold; dis(v) i P j (P is the seed point) j and node v i The plane Euclidean distance; d3: in symmetrical position Centered on a node v, a kd-tree neighborhood search is performed on the supervoxels within three increasing radii to construct node v. i Regarding seed point P j symmetric point set ; d4: Record node v after the search is complete. i Regarding seed point P j The set of all symmetrical points and its corresponding weight list ; d5: Repeat steps d1 to d4 until all supervoxel nodes generate their symmetric lists.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by the processor, the program implements a single-tree segmentation method for lidar point clouds that takes into account morphological optimization, as described in any one of claims 1-7.

9. An electronic device, characterized in that, It includes a memory and a processor; the memory stores a computer program; the processor is used to execute the computer program in the memory to implement the morphologically optimal lidar point cloud single-tree segmentation method according to any one of claims 1-7.

Images