Laser radar point cloud individual tree segmentation method considering optimal form, medium and equipment

Abstract

The invention discloses a laser radar point cloud individual tree segmentation method considering optimal form, a medium and equipment, and the method constructs a frame of initial segmentation, segmentation correction and individual tree form optimization, achieves the individual tree segmentation of the optimal form through a multi-level optimization strategy, gradually improves the segmentation accuracy and form integrity, and improves the segmentation efficiency. And the high-precision tree parameter extraction requirement is met. Specifically, during segmentation correction, according to the single-tree segmentation correction method based on the point density histogram, a wrong segmentation target is detected through point cloud density distribution in the vertical direction, segmentation points are detected through point cloud density distribution in the horizontal direction, and the precision of single-tree segmentation is remarkably improved. During form optimization, a graph cut optimization method of high-order energy constraint is introduced, and global optimization of individual tree forms is realized. Based on individual tree form priori knowledge, a high-order energy function fusing geometric features and form rules is constructed, an individual tree boundary is globally optimized by using a graph cut optimization algorithm, and the problem of crown form error in a traditional method is effectively solved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of photogrammetry and remote sensing, and particularly relates to a laser radar point cloud single tree segmentation method considering optimal shape, a medium and equipment. BACKGROUND

[0002] As the basic unit of forest, the structure parameters of single tree are crucial for forest resource investigation and ecological parameter estimation. Obtaining single tree information plays a key role in forest resource monitoring, carbon storage measurement and biodiversity assessment. At present, the traditional method for obtaining forest structure parameters mainly relies on manual field measurement or aerial image interpretation. Manual field measurement is time-consuming and labor-intensive, and is easily affected by subjective factors. It is difficult to accurately segment low canopy and irregular trees, and cannot meet the needs of continuous, large-scale, accurate and efficient forest resource investigation.

[0003] Traditional single tree segmentation methods are based on laser radar point cloud attributes and rules, and mainly fall into the following three categories:

[0004] (1) CHM grid-based method: By generating a canopy height model (CHM), local maximum detection combined with region growing, watershed and other algorithms are used to segment the tree crown. Although this method is computationally efficient, it is easily affected by tree crown overlap in high canopy density areas and terrain undulations, resulting in boundary errors and loss of details. The limitations of CHM-based single tree segmentation mainly include: ① In high canopy density forest scenes, the large overlap of tree crowns leads to incorrect segmentation of tree crown boundaries; ② In areas with large terrain undulations, normalized CHM can cause changes in tree crown shape, and image smoothing strategies can easily lead to loss of tree crown details. Therefore, this method is prone to incorrect segmentation, and the tree crown shape of the segmentation result is prone to distortion.

[0005] (2) Point-based method: Directly operating on point cloud data, clustering (such as K-means, mean shift) or trunk detection strategies are used, which preserves complete stand information, but has large computational redundancy, poor robustness, and relies on parameter settings; the limitations of point-based methods are: the amount of data involved in the calculation is much larger than the amount of data contained in CHM grid data, causing low computational efficiency.

[0006] (3) Voxel-based method: voxelization of point cloud data is used for feature classification and segmentation, balancing computational efficiency and data volume, but is sensitive to voxel size, and is prone to missing small trees or missegmenting large trees. Voxelization of point cloud data can reduce data redundancy, but the data volume is still a challenge compared to grid data.

[0007] Traditional methods have generality, but are generally faced with the defects of over-segmentation, under-segmentation and scene adaptability limitations. In recent years, machine learning technology has become a new direction to break through the traditional bottleneck with the advantages of high precision, strong generalization and automation. The single tree segmentation method based on machine learning significantly improves the accuracy through data-driven strategy, which mainly includes the following two types:

[0008] (1) Traditional machine learning (such as random forest): The trunk detection is optimized through feature engineering (such as geometric confidence), and the accuracy reaches 88.5%, but it is difficult to distinguish geometrically similar objects (such as street lamp poles and tree trunks). Machine learning method depends on the setting of initial learning threshold, and the migration is not strong.

[0009] (2) Deep learning method: Direction vector prediction (such as PED-net): Combined with semantic segmentation and voxel region growing, it is suitable for complex scenes, but the generalization is weak; 3D sparse convolution: Using circular detection box to match tree shape, accuracy 82.9% (IoU>0.5), efficient calculation, but limited by the lack of airborne point cloud features; Multi-stage network: Fusion of centroid prediction, attention mechanism and dynamic clustering, reduces the dependence on labeling, but the module cooperation is insufficient; Universal model (such as 3D CNN): Support multi-platform data, strong cross-scene applicability, but low accuracy in broad-leaved forest multi-layer crown segmentation. SUMMARY

[0010] The purpose of the present application is to provide a laser radar point cloud single tree segmentation method considering optimal shape, medium and equipment, which aims to solve the above problems.

[0011] The present application is mainly realized by the following technical solutions:

[0012] A laser radar point cloud single tree segmentation method considering optimal shape, comprising the following steps:

[0013] Step S1: preliminary segmentation; collecting unmanned aerial vehicle point cloud and ground point cloud of the target area, and processing to obtain crown height model CHM; based on the crown height model CHM, extracting preliminary segmentation cluster;

[0014] Step S2: segmentation correction;

[0015] Step S21: using the characteristics of tree point cloud in vertical dimension to extract the number of trunks of a single segmentation cluster; detecting the error segmentation target through the vertical direction point cloud density distribution;

[0016] Step S22: segmentation correction according to error type;

[0017] For clusters containing multiple trunks, perform adhesion segmentation; detect the center point line of the trunk, and use the point cloud density distribution in the horizontal direction to detect the segmentation point; for clusters without trunk, perform point cloud fusion;

[0018] Step S3: Optimizing the single tree morphology using a high-order energy-constrained graph cut optimization algorithm;

[0019] Step S31: performing an octree region growing algorithm on the point cloud to construct a supervoxel and Markov field graph structure;

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

[0021] Step S33: Graph cut optimization is performed using the α-expansion algorithm to obtain the optimal single tree morphology; the α-expansion algorithm maps the image into a network graph, constructs an energy function about the label, and then minimizes the energy function; the energy function is:

[0022]

[0023] in is the data item energy function;

[0024] D i (l i ) is a data item describing the supervoxel v i Assign label l i the price;

[0025] V is the voxel;

[0026] is the smoothing energy function;

[0027] S ij (l i , l j ) is a smoothing term that describes the supervoxel v with adjacent relationships i and v j Consistency cost of assigning labels;

[0028] l i and l j are the supervoxels v of the adjacency relationship i and v j Assign labels;

[0029] λ smooth is the smoothing parameter;

[0030] e ij is the connected voxel v i and v j edge;

[0031] E is energy and;

[0032] is a high-order energy function;

[0033] is a weight in the weight list;

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

[0035] S is a symmetric constraint range.

[0036] To better realize the present application, further, in the step S33, the data item D i (l i ) is:

[0037]

[0038] wherein: d i is the plane Euclidean distance from the node v i to the seed point;

[0039] d t is the average radius of the tree crown in the scene;

[0040] σ1 and σ2 respectively control the decay rate of the two ends;

[0041] C is a continuity constant;

[0042] The smoothing term S ij (l i , l j ) is:

[0043]

[0044] wherein: w ij is the edge weight.

[0045] To better realize the present application, further, the step S33 comprises the following steps:

[0046] e1: initialize the label, assign an initial label to each node, and use the initial segmentation result of the step S1 to assign values to the nodes;

[0047] e2: build a graph model, and calculate the initial total energy E current using the energy function E(L);

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

[0049] e4: build an expansion graph, define the source node as the label α, and define the sink node as the set of all non-α labels; define the data item edge as the source edge and the sink edge, wherein the source edge is the edge connected to the source node of each node v i , and represents the cost of selecting α as the label of the node, and the sink edge is the edge connected to the sink node of each node vi The edge connected to the sink node represents the cost of the node selecting non-α as the label; the smoothing item edge retains the original definition;

[0050] e5: using the maximum flow algorithm to solve the maximum flow in the extended graph, corresponding to the minimum cut of the graph, and dividing the graph into the part connected to the source node α and the part connected to non-α;

[0051] e6: updating the label assignment according to the minimum cut result, if v i is connected to the source node, then the label l i of v new = α, otherwise, retaining the original label;

[0052] e7: calculating the updated energy E new , if E new < E current , then retaining the label assignment of this expansion, and letting E new = E current ;

[0053] e8: repeating the steps e3-e7, performing α-expansion for each label, and gradually approaching the global optimal solution, if the energy no longer decreases or the maximum iteration number is reached after several iterations, then terminating the algorithm and outputting the final label assignment.

[0054] In order to better realize the present application, further, the step S1 comprises the following steps:

[0055] Step S11: data acquisition, rasterizing the forest point cloud and the ground point to obtain a digital surface model (DSM) and a digital elevation model (DEM), and subtracting the elevations in the corresponding grids to obtain a canopy height model (CHM);

[0056] Step S12: extracting local extreme points of the canopy height model (CHM) as seed points for region growing;

[0057] Traversing the model CHM and marking the candidate seed points; then, calculating the canopy protrusion degree, and retaining the points with the canopy protrusion degree greater than a set threshold value, finally, sequentially performing spatial deduplication and edge rejection to screen the seed points;

[0058] Step S13: performing region growing to extract the preliminary segmentation clusters.

[0059] In order to better realize the present application, further, the step S21 comprises the following steps:

[0060] a1: high-layer slice construction, performing vertical slicing on the primary segmentation single-tree point cloud, and calculating the number D(z) of points in each layer;

[0061] a2: using the normalized difference method to identify the demarcation point elevation of the trunk-canopy, and setting the elevation density distribution sequence as:

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

[0063] Wherein: D(z n ) is the number of points of the nth layer;

[0064] z min is the lowest height considered;

[0065] The slice step Δz is set to 0.1m, and the first-order forward difference is calculated:

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

[0067] Based on the difference sequence reflecting the density change rate between adjacent elevation layers, a positive transition appears at the boundary between the trunk and the crown; and the trunk layer height z trunk is located, and the trunk layer point cloud subset P trunk = {p i |z i ≤z trunk} is extracted.

[0068] a3: the trunk is preliminarily extracted by Euclidean clustering;

[0069] The Euclidean distance between points is calculated, and neighborhood search is performed on each point, and clustering is performed on all points; and the clustering result is filtered based on height constraint and principal component analysis constraint; the principal component analysis constraint includes principal direction constraint and linear feature constraint.

[0070] In order to better realize the present application, further, the step S22 comprises the following steps:

[0071] b1: the cluster containing multiple trunks is segmented by adhesion; the trunk center point connecting line is detected, and the segmentation point is found;

[0072] First, the point cloud is rotated around the Z axis, so that the axis is parallel to the X axis, the interval is divided, and the number of points falling in each interval is counted to obtain an under-segmented histogram, and the histogram contains two maximum value points; then, the minimum value points falling in the interval of the two maximum value points are taken as the cutting points;

[0073] b2: the point cloud fusion is performed on the cluster without trunk;

[0074] For the over-segmented single tree t i without trunk, all adjacent single trees T i containing a trunk are found,

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

[0076] In order to better implement the present invention, further, step S31 includes the following steps:

[0077] c1: Construct a point cloud octree and use the point similarity h 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 Represents the angle between two voxel normal vectors;

[0081] h Intensity is the difference between the average intensity values ​​of the points within two voxels;

[0082] c2: perform region growing to generate supervoxels;

[0083] c3: construct the Markov field graph structure;

[0084] Each supervoxel is treated as a node and given a unique label l i ∈{1, 2, ..., N}, where N is the number of tree trunk seed points; the initial value of the label is the corrected segmentation result 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.

[0085] In order to better implement the present invention, further, step S32 includes the following steps:

[0086] d1: Take a supervoxel node v from the supervoxel node list i ;

[0087] d2: computing node v i About the seed point P j Center symmetrical position

[0088]

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

[0090] d t is a prior threshold value;

[0091] dis(v i , P j ) is the planar Euclidean distance between the seed point P j and the node v i ;

[0092] d3: perform a kd-tree neighborhood search on the super voxel with the symmetric point set as the center in three increasing radii, and construct the node v i about the symmetric point set of the seed point P j ;

[0093] d4: record all the symmetric point sets of the node v i about the seed point P j and the corresponding weight list ;

[0094] d5: repeat steps d1-d4 until all super voxel nodes generate their symmetric lists.

[0095] A computer readable storage medium, having stored thereon a computer program, which, when executed by a processor, implements the above-mentioned laser radar point cloud single tree segmentation method considering the optimal shape.

[0096] An electronic device, comprising a memory and a processor; the memory has stored thereon a computer program; the processor is configured to execute the computer program in the memory to implement the above-mentioned laser radar point cloud single tree segmentation method considering the optimal shape.

[0097] The beneficial effects of the present application are as follows:

[0098] (1) The present application uses the elevation density feature to extract the single tree trunk to determine the error segmentation, and introduces the high-order energy constraint graph cut optimization method, so that the number is accurate, and the single tree segmentation is complete. In view of the problems of low segmentation accuracy and poor single tree shape of the existing single tree segmentation method, the present application establishes a single tree segmentation process framework of "initial segmentation-segmentation correction-shape optimization", and realizes the single tree segmentation with the optimal shape through a multi-level optimization strategy, gradually improves the accuracy and shape integrity of the segmentation, and meets the demand of high-precision tree parameter extraction.

[0099] (2) To solve the problem of insufficient initial segmentation accuracy, a single tree segmentation correction method based on point density histogram is proposed according to the characteristics of ground point cloud. The vertical direction point cloud density distribution is used to detect the error segmentation target, and the horizontal direction point cloud density distribution is used to detect the segmentation point, which significantly improves the accuracy of single tree segmentation. To solve the problem that traditional segmentation methods cannot consider the shape of single tree, a high-order energy constraint graph cut optimization method is introduced to realize the global optimization of single tree shape. Based on the prior knowledge of single tree shape, a high-order energy function is constructed by combining geometric features and shape rules, and a graph cut optimization algorithm is used to globally optimize the single tree boundary, effectively solving the problem of incorrect tree crown shape in traditional methods.

[0100] (3) The present application can effectively optimize the segmentation accuracy and tree crown shape integrity of the traditional CHM region growing segmentation method. The present application relies on two key components: (1) point cloud elevation density distribution feature guided segmentation correction (2) high-order energy constraint tree crown shape graph cut optimization method. Unlike traditional single tree segmentation methods based on CHM grid, this method is more stable in complex forest environments, has higher single tree segmentation accuracy, and effectively solves the problem of tree crown shape distortion in traditional segmentation methods. BRIEF DESCRIPTION OF DRAWINGS

[0101] Figure 1 The flowchart of the single tree segmentation method of the present application considering the optimal shape of laser radar point cloud;

[0102] Figure 2 The vertical direction point density distribution histogram generated in Example 1;

[0103] Figure 3 The histogram of under-segmentation in Example 1;

[0104] Figure 4 Example 1 of single tree segmentation before and after optimization;

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

[0106] Figure 6 Example 3 of single tree segmentation before and after optimization;

[0107] Figure 7 Example 4 of single tree segmentation before and after optimization. DETAILED DESCRIPTION

[0108] Example 1:

[0109] A single tree segmentation method of laser radar point cloud considering the optimal shape, as shown in Figure 1 , specifically comprising the following steps:

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

[0111] UAV and ground-based point clouds are collected for the target area, corresponding to the upper canopy point cloud and the trunk and lower vegetation point cloud, and the two are registered and fused, and the digital elevation model DEM of the region is extracted as the ground point.

[0112] (1) Data collection and canopy height model CHM generation;

[0113] a. Data collection: UAV LiDAR: flight height 50m, scanning angle 70°, point density ≥50 points / ㎡, obtain upper canopy point cloud; ground-based LiDAR: scanner erection height 1.5m, 360° full station scanning, point density ≥200 points / ㎡, obtain trunk and lower vegetation point cloud.

[0114] UAV and ground-based point cloud registration; manual registration: feature matching: select at least 3 groups of same name feature points (such as prominent trunk bifurcation points, rock top points) in two types of point clouds; ICP fine registration: use iterative closest point algorithm (ICP), set maximum iteration number 100 times, registration error threshold ≤0.1m; fusion verification: calculate the coincidence degree of overlapping area point cloud, require root mean square error (RMSE) <0.15m.

[0115] Ground point extraction and DEM generation; use progressive triangulation encryption filter (PTD), set slope threshold 8°, window size 1m×1m; generate digital elevation model (DEM), grid resolution 0.25m, interpolation method uses inverse distance weighting (IDW).

[0116] b. DSM construction: rasterize non-ground point cloud (0.25m resolution), take the maximum elevation value of point cloud in each grid;

[0117] c. Canopy height model CHM construction:

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

[0119] where CHM, DSM, DEM are all grid images, the value of the pixel represents the point cloud elevation of the corresponding position of the current pixel, (i,j) represents the image coordinates, if negative value appears, set to 0. Gaussian filter: use 3x3 convolution kernel to smooth CHM (σ=0.5), retain the main peak characteristics of the canopy.

[0120] (2) With seed point extraction;

[0121] Use 3x3 sliding window to traverse CHM, mark the center point of the window as the candidate seed point under the following conditions:

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

[0123] Seed point screening: Canopy convexity calculation:

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

[0125] Keep points with C > 0.7 and remove pseudo-extremes in flat areas. Spatial deduplication: If the distance between multiple extreme points is less than 1 meter, only the highest point is retained. Edge culling: Remove points less than 2 meters from the CHM boundary to avoid edge effects.

[0126] (3) Perform region growing to extract preliminary segmentation clusters;

[0127] ① The local maximum point is used as the starting seed point for the first batch of growth, and the eight pixels around the seed point of each growth are defined as the neighborhood pixels of the seed point. Each growth only performs conditional judgment on the neighborhood pixels of the current seed point. The neighborhood pixels that meet the conditions are merged into the area of ​​the seed point and serve as the seed point for the next cycle. Only when the eight neighborhood pixels of the current seed point have all undergone conditional judgment will the next cycle be entered.

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

[0129] H max is the height threshold of the crown vertex in the current area. Therefore, the height relationship between the growth point and the crown vertex in the current area should satisfy:

[0130] H max *0.55<h i <H max , in the formula h i Indicates the height of the point to be grown.

[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 height relationship between the point to be grown and the current seed point should satisfy the following formula, where h i is the height of the growth point, h j is 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 vertex in the current area, using the external input threshold parameter d t to control its growth range.

[0134] ⑤ Competition rules for adjacent pixels. If overlapping areas appear during region growing, set a reasonable competition rule 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 crown vertices, and calculate the height differences H1 and H2 from the point to the two crown vertices. Calculate the ratio of these two parameters. If the ratio of the height difference H1 to the distance D1 from the point to be grown to the crown vertex of Tree1 is less than the ratio of the height difference H2 to the distance D2 from the point to the other crown vertex, then assign the point to Tree1; otherwise, assign it to another seed point.

[0135] Step S2: Cluster segmentation and correction: Analyze the elevation density distribution characteristics of a single cluster, determine the type of misidentification and correct it.

[0136] (1) Since the tree point cloud shows a significant hierarchical characteristic 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: vertically slice the initially segmented single-tree point cloud with a step size of 0.1m and calculate the number of points in each layer:

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

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

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

[0141] z is the height value, the vertical coordinate;

[0142] Δz is the slicing step length; Δz = 0.1m;

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

[0144] b. Mutation threshold detection uses the standardized difference method to identify the elevation of the trunk-crown boundary point. The core of this method is to identify the discontinuous transition of the point cloud density in the vertical direction. Assume that the elevation density distribution sequence is:

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

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

[0147] z min is the lowest height considered;

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

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

[0150] The difference sequence reflects the rate of change of density between adjacent elevation layers, and a significant positive jump will occur at the trunk and crown boundary, as shown by the red arrow in the position. Figure 2

[0151] Through the above steps, the trunk layer height z trunk can be accurately located, and the trunk layer point cloud subset P trunk = {p i |z i ≤z trunk} is extracted. In the next step, the trunk is preliminarily extracted by Euclidean clustering. Euclidean clustering is a point cloud data clustering method based on Euclidean distance measurement, and its core idea is to classify points with similar Euclidean distances in point cloud as the same class. For two points p i (x i , y i , z i ) and p j (x j , y j , z j ) in point cloud space, the Euclidean distance is defined as:

[0152]

[0153] This distance measurement reflects the actual geometric distance between points in Euclidean space, and then the neighborhood search is performed for each point p i and the neighborhood radius r, and the 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. The neighborhood search is performed by calculating the distance between points p i ​The 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 the present invention is 0.2m. By performing a neighborhood search on all points, all points can be marked as a cluster. The minimum distance between all points in the cluster is less than the neighborhood radius.

[0156] In order to prevent shrubs and weeds close to the ground from affecting the accuracy of tree trunk extraction, the clustering results need to be further filtered to remove shrubs and weeds that are mistakenly identified as tree trunks. The filtering conditions used here are as follows:

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

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

[0159] Set the height threshold H t , if H m ≥H t , then the candidate trunk is retained, and the height threshold H used in this invention is t =1.2m.

[0160] ② Principal component analysis constraints: These constraints are divided into principal direction constraints and linear feature constraints. Their purpose is to preserve the vertical and approximately straight trunk structure. A cluster is considered to meet the principal component constraints only if it satisfies both the principal direction constraints and the linear feature constraints.

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

[0162]

[0163] Among them: a 1， b1 and c1 are three directional components.

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

[0165] The linear feature can be calculated by the three eigenvalues ​​λ1, λ2, and λ3 obtained by principal component analysis, and the linear index is defined as: L m =(λ1-λ2)÷λ3; It is stipulated that when L m When ≥0.7, the cluster is considered to be approximately a straight line.

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

[0167] a. Segmentation of multiple tree trunks. Samples containing multiple tree trunks after initial segmentation are considered under-segmented and should be optimized. Lines connecting the center points of the tree trunks should be detected to find appropriate segmentation points. To facilitate histogram calculation, the point cloud should first 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 trunks 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:

[0172]

[0173] Where: R is the rotation matrix;

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

[0175] like Figure 3 As shown, after rotation, the X-axis is divided into intervals with a step size of Δs = 0.1m. The number of points falling within each interval is counted to obtain an under-segmented histogram, which should contain two maximum points (tree trunk points). A minimum point will appear between the two maximum points. The minimum point falling within the interval formed by the two maximum points 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. Merge point clouds without trunks. After the under-segmentation optimization is completed, the next step is to optimize the over-segmented trees that do not contain trunks. Find all trees adjacent to the tree t that contain a trunk and form a set T. Calculate the Euclidean distance between the geometric center of all trees in T and the geometric center of tree t. Then sort them into the optimization queue according to the distance and then determine whether to merge them. In theory, we need to find an optimal tree T for tree t in the set T. i , the relationship between the two is t belongs to T i , find T i To ensure the morphological similarity of the merged trees, the 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 center of the nearest trunk should not exceed twice the average crown radius of the scene.

[0177] Step S3: high-order energy-constrained graph cut optimization algorithm is used to optimize the single tree shape.

[0178] (1) The region growing algorithm of octree is used to the point cloud to construct the super voxel and the graph structure of Markov field.

[0179] a. Construct the point cloud octree, first randomly select several seed points P seed in the leaf nodes of the octree, and try to ensure that the seed points are evenly distributed in space. The similarity of points is used to describe the characteristics of the octree nodes, and the point similarity h is measured by three characteristics: spatial distance, intensity value and normal vector:

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

[0181] Where: h coord represents the Euclidean distance calculated according to the three-dimensional coordinates of the voxel center point, h Normal represents the angle between the normal vectors of the two voxels, and h Intensity is the difference between the average intensity values of the two voxel points.

[0182] After obtaining the three similarity characteristics, in order to balance the influence of the three characteristics on similarity, normalization processing is needed:

[0183]

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

[0185] b. Generate super voxels by region growing, and the steps of region growing are as follows:

[0186] ① Add all seed points P seed to the growing queue;

[0187] ② Take out a seed point P i from the growing queue as the current growing seed point each time, search all voxels in its neighborhood using the octree, and generate a search queue for the seed point;

[0188] ③ Traverse the neighborhood voxels p i in the search queue, if the search queue is empty, it is considered that the super voxel has reached the maximum voxel radius, remove it from the growing queue and return to step ②.

[0189] ④Calculate the seed point P i and its neighboring voxel p i If the similarity is less than the threshold, the voxel p is marked. i Belongs to seed point P i of supervoxels and add the p i In the search sequence from the adjacent voxels to the seed point, the adjacent voxel points and the seed point P are required to be i The distance is smaller than the search radius R of the supervoxel;

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

[0191] c. Construct a Markov field graph structure.

[0192] Graph node: Each supervoxel is a node and is given a unique label l i ∈{1, 2, ..., N} (N is the number of tree trunk seed points). The initial value of this label is the corrected segmentation result obtained in step 2. Edge connection: 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.

[0193] (2) Calculate the central symmetry point of each supervoxel node about all tree trunk points and generate a symmetry list for multi-level search. The steps to generate the node symmetry table are as follows:

[0194] ① Take a supervoxel node v from the supervoxel node list i .

[0195] ②Calculate the node v i About a seed point P j Center symmetrical position Usually no nodes fall in exact symmetrical positions, so Is a virtual location node, if the node v i With seed point P j If the plane Euclidean distance of the node exceeds a certain threshold, the node will be j The symmetric point of is recorded as 0, and the calculation method is:

[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, d t is the prior threshold, usually set to the average radius of the tree crown in the scene, dis(v i, P j ) is the seed point P j and node v i The planar Euclidean distance of .

[0198] ③ In symmetrical position As the center, perform kd-tree neighborhood search on the supervoxel within three increasing radii to construct node v i About the seed point P j The symmetric point set of The weights of the symmetrical points within the three increasing radii are different and are designed as follows:

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

[0200] The second layer can tolerate slight deformation or data noise, and the search radius of this area is r2 = 0.3d t If the region contains a node and the node label is the same as v i The labels are different, so the weight w2=2 is used, which means a medium energy penalty.

[0201] The third layer is an asymmetric area, which is used when there are no nodes in the first and second layers. It can tolerate some missing data or incomplete crown morphology. The search radius of this area is set to r3 = 0.5d t If the region contains a node and the node label is the same as v i The labels are different, and the weight w3=1 is taken to minimize the energy penalty.

[0202] ④After the search is completed, record the node v i For each seed point P j The set of all symmetric points and their corresponding weight lists

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

[0204] (3) A symmetric constrained graph cut optimization algorithm is introduced, and the α-expansion algorithm is used for graph cut optimization to obtain the optimal single tree morphology. The energy function is as follows:

[0205]

[0206] in is the data item energy function, is the smoothing energy function, is a high-order energy function.

[0207] a. Data item constraints. Data item D i (l i ) describes the supervoxel v i Assign label l i The smaller the cost, the more consistent it is with the label l i The higher the consistency, the dynamic distance decay model is proposed for this item, and its function expression is:

[0208]

[0209] Among them, d i For node v i Planar Euclidean distance to the seed point, d t is the average radius of the tree crown in the scene, σ1 and σ2 control the attenuation rates at both ends respectively. The present invention sets these two parameters 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 , is calculated as:

[0211]

[0212] b. Smoothness term constraint. Smoothness term S ij (l i , l j ) describes the supervoxel v with adjacency i and v j The consistency cost of assigning labels uses the traditional Potts model, and its calculation formula is:

[0213]

[0214] w ij is the edge weight, and the edge weight is calculated as:

[0215]

[0216] where d t is a priori parameter representing the average crown radius of trees in the scene.

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

[0218]

[0219] Weight The calculation is obtained by the symmetric list, where δ(li , l j ) for judging label consistency, l i and l j are labels of v i and v j , the value of label consistency δ(l i , l j ) is 0, and vice versa is 1. In order to avoid repeatedly calculating the symmetric energy of a node, in the process of one iteration, for a certain node, the symmetric energy penalty is calculated only in the symmetric point set of the seed point with the same label as the node, and the maximum energy w max in the multi-level symmetric energy is taken as the energy penalty of the iteration.

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

[0221] ①Initialization of label, each node is given an initial label, in the present application, the initial segmentation result obtained above is used to assign values to the nodes, which can effectively reduce the iteration number of the algorithm and save the running time of the algorithm.

[0222] ②Construction of graph model, using the energy function model above, the edge weight (smooth term) of each node and its adjacent node, the energy penalty of the node assigned to a certain label (data term) and the high-order energy penalty under the symmetry constraint are calculated to obtain the initial total energy E current .

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

[0224] ④Convert the multi-label problem into a binary classification problem, i.e. α class and non-α class. Construct an expansion graph, define the source node as label α and the sink node as the set of all non-α labels; define the data term edge as the source edge and the sink edge, the source edge is the edge connecting each node v i to the source node, representing the cost of the node selecting α as the label, the sink edge is the edge connecting each node v i to the sink node, representing the cost of the node selecting non-α as the label; the smooth term edge remains the original definition.

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

[0226] ⑥Update the label assignment according to the minimum cut result, if v i is connected to the source node, then its label l i = α, otherwise, keep the original label.

[0227] ⑦Calculate the updated energy Enew , if E new <E current , then keep the label assignment of this extension, let E new =E current .

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

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

[0230] like Figures 5-7 As shown in Figure 2, (a) is a schematic diagram of single tree segmentation before optimization, and (b) is a schematic diagram of single tree segmentation after optimization. Figure 5 As shown in the red box in Figure (a), the crown is divided in a "one-size-fits-all" manner, without considering the growth morphology between the crowns. In the red box in Figure (b), it is obvious that the crown's branch and leaf morphology is more in line with objective facts, and the branches and leaves of different trees intersect with each other. Figure 6 As shown in the red box in Figure (a) before optimization, the crown segmentation method is "one size fits all" and does not take into account the growth morphology between crowns. In the red box in Figure (b) after optimization, it is obvious that the branch and leaf morphology of the crown is more in line with objective facts, and the branches and leaves of different trees intersect with each other. Figure 7 As shown in the red box in Figure (a) before optimization, the crown shape is obviously unreasonable. The crown shape should be roughly symmetrical, as shown in the red box in Figure (b) after optimization.

[0231] The above description is only a preferred embodiment of the present invention and does not limit the present invention in any form. Any simple modification or equivalent change made to the above embodiment based on the technical essence of the present invention shall fall within the scope of protection of the present invention.

Claims

1. A method for single tree segmentation of a laser radar point cloud considering morphological optimization, characterized in that, Comprising the following steps: Step S1: initial segmentation; collecting unmanned aerial vehicle point cloud and ground point cloud of the target area, and obtaining canopy height model CHM by processing; based on the canopy height model CHM, a preliminary segmentation cluster is extracted; Step S2: segmentation correction; Step S21: the number of tree trunks of a single segmentation cluster is extracted by using the characteristics of tree point cloud in the vertical dimension; the error segmentation target is detected by the vertical direction point cloud density distribution; Step S22: segmentation correction is performed according to the error type; For the cluster containing multiple tree trunks, adhesion segmentation is performed; the tree trunk center point connecting line is detected, and the segmentation point is detected by using the point cloud density distribution in the horizontal direction; for the cluster without tree trunk, point cloud fusion is performed; Step S3: a graph cut optimization algorithm with high-order energy constraint is used for single tree shape optimization; Step S31: the region growing algorithm of octree is performed on the point cloud, the super voxel is constructed, and the graph structure of Markov field is constructed; Step S32: the center symmetry point of each super voxel node with respect to all tree trunk points is calculated, and a multi-level search symmetry list is generated; Step S33: the graph cut optimization is performed by using the alpha-expansion algorithm to obtain the optimal single tree shape; the alpha-expansion algorithm is to map the image into a network graph, construct an energy function about the label, and then minimize the energy function; the energy function is: wherein is the data item energy function; D i (l i ) is a data item, describing a super voxel v i assigning a label l i the cost of V is a voxel; is a smoothing term energy function; S ij (l i ,l j ) is a smoothing term, describing the consistency cost of assigning labels to the super-voxels v i and v j with adjacency relationship; l i and l j are adjacent relationship super-voxels v i and v j assign labels; λ smooth is a smoothing parameter; e ij For connecting voxels v i and v j edges; E is an energy sum; is a high order energy function; w is a weight in the list of weights; delta(i, l i , l j ) is a function for determining label consistency; S is a symmetry constraint range.

2. The method of claim 1, wherein, In said step S33, the data item D i (l i ) is: where: d i is the node v i the planar Euclidean distance to the seed point; d t R is the average radius of the tree crown in the scene; σ1 and σ2 control the decay rates of the two ends, respectively; C is a continuity constant; Smooth term S ij (l i , l j ) is: where: w ij is an edge weight.

3. The method of claim 1 or 2, wherein, The step S33 comprises the following steps: e1: initializing a label, assigning an initial label to each node, and using the initial segmentation result of step S1 to assign the value of the node; e2: construct graph model, calculate initial total energy E using energy function E(L) current ; e3: selecting a label alpha from the label set as an expansion target; e4: construct an extended graph, define the source node as label a, and the sink node as the set of all non-a labels; define the data item edges as the source edges and the sink edges, the source edges are the edges from each node v i the edges connected to the source node represent the cost of the node selecting a as the label, and the sink edges are the edges from each node v i the edges connected to the sink node represent the cost of the node selecting non-a as the label; the smoothing item edges remain the original definition; e5: using the maximum flow algorithm to solve the maximum flow in the expansion graph, corresponding to the minimum cut of the graph, and dividing the graph into the alpha part connected with the source node and the non-alpha part; e6: update label assignment according to min-cut result, if v i Connect source node, then its label l i = a, otherwise keep original label; e7: compute updated energy E new if E new <E current then keep the label assignment of this expansion, let E new = E current ; e8: repeating steps e3-e7, performing alpha-expansion for each label, and gradually approaching the global optimal solution; if the energy does not decrease or reaches the maximum iteration number after several iterations, the algorithm is terminated, and the final label assignment is output.

4. The method of claim 1, wherein, The step S1 comprises the following steps: Step S11: data acquisition, rasterizing the forest point cloud and the ground point cloud to obtain a digital surface model DSM and a digital elevation model DEM, and subtracting the elevations in the corresponding grid to obtain a canopy height model CHM; Step S12: extracting local extreme points of the canopy height model CHM as seed points for region growing; Traversing the model CHM and marking the candidate seed points; then, calculating the canopy protrusion degree, and retaining the points with the canopy protrusion degree greater than a set threshold; finally, performing spatial deduplication and edge elimination in turn to screen the seed points; Step S13: region growing is performed to extract the preliminary segmentation cluster.

5. The method of claim 4, wherein, The step S21 comprises the following steps: a1: high layer slice construction, performing vertical slicing on the initial segmentation single tree point cloud, and calculating the number D(z) of points in each layer; a2: using the normalized difference method to identify the demarcation point elevation of the tree trunk-canopy, and setting the elevation density distribution sequence as: D = {D(z1), D(z2),... D(z n )}, z i = z min +iΔz; where: D(z n ) is the number of points in the nth layer. z min the starting position of the trunk, the lowest point; The slice step length Δz is set to 0.1 m, and the first-order forward difference is calculated: ΔD k = D(z k+1 )- D(z k ), k = 1, 2,... n - 1; Based on the differential sequence reflecting the density change rate between adjacent elevation layers, a positive transition appears at the boundary between the trunk and the crown; and then the trunk layer height z is located trunk , and the trunk layer point cloud subset P trunk ={p i |z i ≤z trunk} is extracted a3: preliminary extraction of tree trunks by Euclidean clustering; Calculate the Euclidean distance between points and points, and perform neighborhood search for each point, cluster all points, and filter the clustering results based on height constraints and principal component analysis constraints, including principal direction constraints and linear feature constraints.

6. The method of claim 5, wherein the method is a laser-radar point cloud single tree segmentation method considering the optimal shape. The step S22 comprises the following steps: b1: adhesion segmentation of clusters containing multiple tree trunks; detecting the center point connection line of the tree trunk and finding the segmentation point; First, rotate the point cloud around the Z axis to make the axis parallel to the X axis, divide the interval and count the number of points falling in each interval to obtain an under-segmented histogram, and the histogram contains two maximum points; then, the minimum point falling in the interval of the two maximum points is taken as the cutting point; b2: point cloud fusion for clusters without tree trunks; For each tree T, find all connected and containing one stem single trees t i , find all connected and containing one stem single trees t i , Calculate single tree T i With a single wooden t i The Euclidean distance of the single tree T i An optimized queue is formed and then fused based on the projection overlap ratio constraint.

7. The method of claim 1, wherein, The step S31 comprises the following steps: c1: constructing an octree of point clouds, and using the similarity h of points to describe the characteristics of the octree nodes; h = h coord + h Normal + h Intensity ; wherein: h coord represents the Euclidean distance calculated from the three-dimensional coordinates of the voxel center point; h Normal represents the angle between two voxel normal vectors; h Intensity is the difference between the average intensity values of the two in-vivo points; c2: region growing to generate a super voxel; c3: constructing a Markov field graph structure; Each super voxel is regarded as a node, and is given a unique label l i ∈{1,2,...,N}N is the number of stem seed points; the initial value of the label is the corrected segmentation result in step S2; if the spatial distance between two super voxels is less than 2m, an edge is established, and the weight is determined by the feature difference.

8. The method of claim 7, wherein, The step S32 comprises the following steps: d1 : take a voxel node v from the list of voxel nodes i ; d2: computing node v i Regarding the seed point P j Centrosymmetric position where (P x , P y ) and (v x , v y , v z ) represent the three-dimensional coordinates of the seed point and the node, respectively; d t is a prior threshold value; dis(v i , P j ) is the planar Euclidean distance between the seed point P j and the node v i ; d3: in symmetric position Perform a kd-tree neighborhood search of the super-voxel centered at v i Regarding the seed point P j The symmetric point set of P d4: record the node v after the search is completed i about the seed point P j all the symmetry point sets and their corresponding weight lists d5: repeat steps d1-d4 until all super voxel nodes generate their symmetry lists.

9. A computer-readable storage medium having stored thereon a computer program, characterized in that, The program is executed by the processor to realize the method of claim 1-8.

10. An electronic device, comprising: A computer device comprising a memory and a processor; the memory has a computer program stored thereon; the processor is configured to execute the computer program in the memory to realize the method of claim 1-8. A computer device comprising a memory and a processor; the memory has a computer program stored thereon; the processor is configured to execute the computer program in the memory to realize the method of claim 1-8.

Images