Airborne and SLAM point cloud automatic registration method for forest sample plot

Abstract

The invention relates to an airborne and SLAM point cloud automatic registration method for a forest sample plot, and belongs to the field of laser point cloud data processing. The method comprises the following steps: cooperatively acquiring point cloud data of a forest land area through an airborne laser radar and SLAM equipment, and performing feature extraction on SLAM point cloud by adopting a random forest model to obtain crown edge contour points; segmenting the airborne point cloud by using a multi-level watershed algorithm and extracting crown edge contour points; using an FPFH algorithm to match position features of computer-borne data crown edge contour points and SLAM data crown edge contour points; and a transformation matrix is calculated through an RANSAC method, and the source point cloud is transformed to the reference of the target point cloud, so that complete spatial information description of the forest land scene is realized. According to the method, the edge contour points of the tree point cloud can be automatically extracted, the airborne point cloud and the SLAM point cloud data are registered, and the forest point cloud vertical structure information is fused, so that the forest structure information is finely inverted.

Description

Technical Field

[0001] The invention relates to an automatic registration method of an airborne and SLAM point cloud of a forest plot, and belongs to the technical field of laser point cloud data processing. Background Art

[0002] Trees are an indispensable material resource for human survival and social development. The precise inversion of forest structure information is of great significance for exploring the growth status of forests and exploring the spatial distribution patterns of forests. LiDAR can obtain detailed forest spatial structure information by emitting laser energy and receiving returned information. For forest point clouds, it is difficult to completely scan the vertical structure information of the forest using a single airborne radar or ground radar. Airborne radar uses a top-down scanning method to accurately obtain complete forest middle and upper structural information, but due to the occlusion of branches and leaves inside the canopy, the vegetation information in the middle and lower layers is incomplete; SLAM equipment uses a bottom-up scanning method to obtain complete forest understory detailed information, but the upper crown information is easily blocked by the middle and lower crowns, and the point cloud data obtained has the problem of incomplete tree vertex clouds.

[0003] In order to solve the above defects, it is very necessary to combine the two. At present, there is an urgent need for an automatic registration method of airborne and SLAM point clouds for forest plots. Summary of the invention

[0004] In view of the shortcomings of the prior art, the present invention provides an automatic registration method of airborne and SLAM point clouds for forest plots, which solves the problem that the existing single point cloud acquisition platform cannot fully invert the vertical structure information of the forest.

[0005] The present invention adopts the following technical solution:

[0006] An automatic registration method of airborne and SLAM point clouds of forest plots, comprising the following steps:

[0007] S1: Use airborne laser radar equipment to collect the tree tops in the forest area and obtain airborne point clouds; use SLAM (Simultaneous Localization and Mapping) equipment to collect the tree bottoms in the forest area and obtain SLAM point clouds;

[0008] S2: Use the random forest model to classify and extract the tree point cloud from the SLAM point cloud, and extract the edge contour points of the crown point cloud;

[0009] S3: Use a multi-level watershed algorithm to segment the airborne point cloud into single tree point clouds and extract edge contour points of the crown point cloud;

[0010] S4: The edge contour points of the tree crown point cloud obtained by the airborne point cloud are used as source feature points, and the edge contour points of the tree crown point cloud obtained by the SLAM point cloud are used as target feature points. According to the coordinates of the source feature points and the target feature points, the FPFH algorithm is used to calculate the relative position relationship of each feature point relative to other points in the neighborhood range, and the position relationship of all relative neighborhood points is weightedly added according to the distance to obtain the position feature of each feature point;

[0011] S5: The edge contour points of the tree crown point cloud obtained from the airborne point cloud are used as source feature points, and the edge contour points of the tree crown point cloud obtained from the SLAM point cloud are used as target feature points. According to the positional features of the source feature points and the target feature points, the size-constrained RANSAC algorithm is used to preferentially match feature points with large tree crown sizes and the same positional features as points of the same name. The point cloud transformation matrix is ​​iteratively calculated based on the points of the same name, and the coordinates of each point in the source feature point are left-multiplied by the point cloud transformation matrix to obtain new coordinate information, thus completing the two-phase registration.

[0012] Preferably, step S2 includes the following sub-steps:

[0013] S21: First, filter and denoise the SLAM point cloud data, and use the Cloth Simulation Filtering (CSF) algorithm to separate the ground points and non-ground points of the SLAM point cloud data to complete the preprocessing; delete the ground points and construct an octree structure for the non-ground points to establish the spatial relationship of the point cloud;

[0014] S22: Use the random forest model to perform feature learning on the preprocessed SLAM point cloud data, and construct a multidimensional feature vector including geometric features and spatial features for the arc-shaped geometric features of the tree trunk point cloud. Use it as input data for the random forest model and train it to identify the tree trunk point cloud.

[0015] S23: Based on the trunk point cloud recognition results, the complete crown point cloud is searched and clustered to obtain the complete single tree point cloud; for the single tree point cloud, it is projected onto a two-dimensional plane, and the contour points are extracted using the alpha-shape algorithm to obtain the edge contour points of the crown point cloud.

[0016] In the process of extracting the edge contour points of the crown point cloud from the SLAM point cloud, the octree structure is first constructed for the point cloud data to establish the spatial relationship of the point cloud; then, the parameters of the random forest are set, and the random forest model is trained using the eigenvector matrix; finally, the crown point cloud is searched by neighborhood clustering based on the trunk point cloud data predicted by the model, the crown point cloud is projected onto a two-dimensional plane, and the alpha-shape algorithm (existing algorithm) is used to extract the contour points to obtain the edge contour points of the crown.

[0017] Preferably, step S3 includes the following sub-steps:

[0018] S31: De-noise the airborne point cloud data, perform ground filtering through cloth simulation to obtain ground points PC_2 and non-ground points PC_3; generate a digital surface model DSM and a digital terrain model DEM through interpolation, and obtain a canopy height surface model CHM by subtracting the two;

[0019] S32: normalize the canopy height surface model CHM to a grayscale space of 0-255, perform the first watershed segmentation using the foreground of the entire layer as a marker, and perform the second watershed segmentation using the upper layer of the segmented independent area as a marker;

[0020] S33: Taking the local highest point as the center and the median of the boundary detection as the crown radius constraint, the under-segmented area is segmented for the third time to obtain the edge contour points of the crown point cloud.

[0021] Multi-level watershed refers to the use of multi-level segmentation and marker control methods to effectively reduce the over-segmentation phenomenon that is prone to occur in traditional watershed methods. The first watershed uses the foreground of the entire layer as a marker to preliminarily divide different vegetation areas; the second watershed uses the upper layer of the independent area after segmentation as a marker to further subdivide the vegetation area; the third segmentation uses the local highest point as the center and combines the crown radius constraint to further refine the under-segmented area. Through a multi-level, gradually refined segmentation strategy, single trees can be accurately segmented. The single tree point cloud is clustered and projected onto a two-dimensional plane, and the alpha-shape algorithm is used to extract contour points to obtain the edge contour points of the crown.

[0022] Preferably, step S32 is specifically:

[0023] The canopy height surface model CHM is normalized to a grayscale space of 0-255 to generate a grayscale image srcA, srcA is binarized to obtain a binary image srcB, with pixel values ​​of 0 and 1, and the binary image srcB is morphologically eroded to remove small spots and noise in the binary image srcB for full-layer foreground marker extraction of the first watershed;

[0024] Cluster all foreground pixels with a value of 1 in the binary image srcB into multiple pixel clusters, and provide a mark index for each pixel cluster. The above information is recorded through the mark map marks;

[0025] The watershed segmentation method based on the marker map is used for segmentation and the segmentation result is optimized by the grayscale map. The final segmentation result is recorded in marks. Based on the marks result, the area Area of ​​each region is calculated according to the formula by counting the number of regional pixels gridsum. The maximum values ​​gradmax and gradmin of the regional pixels are counted to determine whether the pixel is located in the upper layer of the region. When the area Area of ​​the region is greater than the threshold NormalCrown, the upper layer pixels are recorded in the binary map srcTreeC with a value of 1. When the area Area of ​​the region is less than the threshold NormalCrown but greater than the threshold MinCrown, the pixels of the entire layer are recorded in the binary map srcTreeC with a value of 1, and the pixel values ​​of the remaining positions are 0. Binarization, erosion, clustering and the second watershed segmentation based on markers are performed again, and the segmentation results are recorded in marks2.

[0026] Step S33 is specifically as follows:

[0027] Taking the tree bottom points extracted from the SLAM point cloud as the position constraint, the under-segmented area in marks2 is segmented for the third time, and the canopy diameters searched along four directions from the highest point are counted. The four directions refer to the east, south, west, and north (up, down, left, and right of the node) directions along the horizontal plane orthogonal to the highest point. 1 / 2 of the median is taken as the canopy radius constraint threshold radius (the threshold radius refers to arranging the searched canopy diameters from large to small, and selecting half of the median as the threshold radius). The pixel values ​​in the area whose distance from the highest point is greater than the threshold radius are set to -1; the point cloud in each segmented area recorded in marks2 is assigned a color as a separate tree, the single tree segmentation of the original point cloud is completed, and the edge contour points of the crown point cloud are obtained.

[0028] Preferably, step S4 includes the following sub-steps:

[0029] S41: For the edge contour points of the tree crown point cloud obtained by the airborne point cloud and the SLAM point cloud, let P s Point is the point where the position feature needs to be calculated. t Point P s A point in the neighborhood of is defined as a local coordinate system (u, v, w) based on the normal vector of the two points and the vector of the line connecting the two points. The calculation formula of the three axes of the local coordinate system is as follows:

[0030] u=n s

[0031]

[0032] w=u×v

[0033] where n s For point Pt The normal vector, P t -P s For point P s To P t The vector of t -P s ‖ 2 Indicates P t and P s The Euclidean distance between them is based on the local coordinate system, and the normal vector n is represented by a set of angles. s and n t The deviation between them is calculated as follows:

[0034] α=v·n t

[0035]

[0036] θ=arctan(w·n t ,u·n t )

[0037] Among them, α is P t The normal vector n t The angle between the coordinate axis v and the coordinate axis; φ is P s The normal vector n s The angle between the two points and the line; θ is the normal vector n t In u P t The angle between the projection of the w plane and the u direction;

[0038] Therefore, the surface geometric features of the object can be described by the three parameters (α, φ, θ), and the calculated surface geometric features can be put into the histogram and divided into 11 intervals to form a simplified point histogram (Simplified Point Feature Histogram) descriptor between the source point and a neighboring point, that is, the SPFH descriptor; the SPFH descriptors of all neighboring points within the neighborhood of the source point are calculated to obtain the feature descriptor of each neighboring point relative to the source point, and the formula is:

[0039]

[0040] S42: Take the weighted average of the SPFH descriptors of the i-th source point, and use w as the weight coefficient of its j-th neighbor point ij Represented by, the FPFH descriptor formula of the i-th point is:

[0041]

[0042] Among them, k is the number of points in the neighborhood of the i-th point;

[0043] The FPFH descriptors of the crown edge contour points of the extracted airborne point cloud and SLAM point cloud are calculated respectively to obtain the position features of each point, thereby providing the point position features required for the matching of the same-name point pairs for the RANSAC algorithm registration in step S5.

[0044] Preferably, step S5 includes the following sub-steps:

[0045] S51: Calculate the minimum enclosing rectangle for each tree crown point cloud edge contour point separately, based on the predefined crown area threshold (the crown diameter of adult large and medium trees is about 6 to 12 meters, and the crown diameter of adult small trees is about 3 to 5 meters. The crown area should actually be the area of ​​the crown enclosing rectangle, so the crown area threshold can be taken as 25m 2 ) classify them into large crown sets and small crown sets;

[0046] S52: Estimate the point cloud transformation matrix T through the size-constrained RANSAC algorithm.

[0047] Preferably, in step S51, the process of calculating the minimum enclosing rectangle is:

[0048] A. Sort the edge contour points of the crown point cloud to obtain ordered edge contour points;

[0049] B. Based on the ordered edge contour points, for each edge contour point, select the next neighboring point adjacent to it to form a straight line, and use the straight line and the straight line perpendicular to it to determine the direction of the minimum enclosing rectangle, and calculate a minimum enclosing rectangle;

[0050] C. Calculate the vertical distance from the edge contour point to the four sides of the minimum bounding rectangle. If the vertical distance to the four sides is greater than the set threshold (the threshold is 12.5% ​​of the shortest side of the rectangle), it is a non-overlapping point, otherwise it is an overlapping point; the minimum bounding rectangle determined by any adjacent points is selected from the bounding rectangle with the largest overlap, and the ideal bounding rectangle can be obtained as the final minimum bounding rectangle;

[0051] D. Let the source feature point be P, the target feature point be Q, and each group of point clouds is classified according to the size of the crown area to form a set P = {P L ,P S} and Q = {Q L ,Q S}, where P L , Q L represents the large-scale crown set, P S , Q S Represents the set of edge points of small tree crowns.

[0052] Preferably, the specific steps of step S52 are as follows:

[0053] a. Find P through the nearest neighbor matching of SPFH descriptor L The three random sampling points in Q L The corresponding points in Q L The corresponding points in are taken as the initial three pairs of points with the same name;

[0054] b. According to P L and Q L The three pairs of points with the same name {(p i ,q i )}, substituting its coordinates into the calculation can make P L Transform the same-name point in Q L In the middle, click p i After rotation and translation, it is transformed into point q i , through p i With q i The coordinates of the points can be used to calculate the rigid transformation relationship between the points, that is, the rotation matrix R and the translation vector t:

[0055] q i =Rp i +t

[0056] Through the rotation matrix R and the translation vector t, the matrix T of the assumed change of the same name point can be obtained:

[0057]

[0058] Among them, p i represents the coordinates of the i-th point in the source set of points with the same name, q i represents the coordinates of the i-th point in the target set of points with the same name, R is a rotation matrix used to represent the rotation part in the rigid transformation; t is a translation vector used to represent the translation part in the rigid transformation;

[0059] c. L The coordinates of each point in the left multiply the assumed change matrix T to obtain the transformed P L The source feature point coordinates;

[0060] d. Through the Euclidean distance threshold (the Euclidean distance threshold can be 1m), at the converted source feature point P L and target feature point Q L Between, calculate Q L All points are the center of the circle, and the Euclidean distance threshold is the radius. The transformed P L In Q L The points within the range are regarded as inliers. If the number of inliers is less than 95% of the number of source feature points, return to step a and exclude the combination of these three points. Otherwise, the matrix is ​​considered to be a transformation matrix.

[0061] e. If the number of inliers after multiple iterations is always less than 95% of the number of source feature points, then P is introduced. S , Q S Supplement the points in and repeat steps ad;

[0062] f. Calculate the distance between the pairs of points with the same name. If the distance reaches the minimum value so far, it will be used as the final transformation matrix T.

[0063] The present invention adopts the size-constrained RANSAC algorithm, randomly and preferentially selects three points in the set with larger crown sizes from the source feature points, matches the corresponding three points with similar crown size information and the same position features in the target point cloud, calculates the transformation matrix, and transforms the source point cloud to the target point cloud benchmark through the matrix. The number of source point cloud points within the target point cloud distance threshold is counted, and it is continuously iterated until the number of points no longer increases, and then the points in the set with smaller crown sizes are introduced to participate in the matching of feature points, and the iteration is continued until the number of source point cloud points within the target point cloud distance threshold no longer increases.

[0064] For any details not provided in the present invention, please refer to the prior art.

[0065] The beneficial effects of the present invention are:

[0066] 1. The present invention can efficiently and accurately automatically extract the crown edge contour points of the airborne and SLAM point clouds, greatly reducing the burden of manual operation and improving the efficiency and accuracy of data processing. At the same time, the automatically extracted crown edge contour point data provides more reliable and accurate basic data for subsequent tree point cloud registration.

[0067] 2. The present invention can realize the accurate registration of point cloud data acquired by airborne radar and SLAM equipment in the same area, solve the problem of integrating vertical structural information between point cloud data from different sources, and enable users to comprehensively utilize multiple data sources for more comprehensive analysis and research. The registered point cloud data improves the integrity and accuracy of spatial information. BRIEF DESCRIPTION OF THE DRAWINGS

[0068] The drawings in the specification, which constitute a part of the present application, are used to provide further understanding of the present application. The illustrative embodiments of the present application and their descriptions are used to explain the present application and do not constitute improper limitations on the present application.

[0069] Figure 1 It is a flow chart of the method for automatic registration of airborne and SLAM point clouds of forest plots of the present invention;

[0070] Figure 2 A schematic diagram of local coordinates established by the FPFH algorithm of a certain embodiment of the present invention;

[0071] Figure 3A schematic diagram of the distribution of source points and neighborhood points in an embodiment of the present invention;

[0072] Figure 4 Schematic diagram of the minimum outer bounding rectangle of an embodiment of the present invention. DETAILED DESCRIPTION

[0073] In order to enable those skilled in the art to better understand the technical solutions in this specification, the technical solutions in the embodiments of the present invention are clearly and completely described below in conjunction with the drawings in the implementation of this specification, but are not limited to this. Anything not fully described in the present invention shall be based on the conventional technology in the art.

[0074] Example 1

[0075] An automatic registration method for airborne and SLAM point clouds of forest plots, such as Figure 1 As shown, the following steps are included:

[0076] S1: Use airborne laser radar equipment to collect the tree tops in the forest area to obtain airborne point clouds; use SLAM equipment to collect the tree bottoms in the forest area to obtain SLAM point clouds;

[0077] S2: Use the random forest model to classify and extract the tree point cloud from the SLAM point cloud, and extract the edge contour points of the crown point cloud;

[0078] S3: Use a multi-level watershed algorithm to segment the airborne point cloud into single tree point clouds and extract edge contour points of the crown point cloud;

[0079] S4: In the Gaussian plane coordinate system, 120 degrees east longitude is the central meridian. The intersection of the central meridian and the equator is taken as the coordinate origin, the projection of the central meridian is taken as the ordinate axis X, the X-axis is set to the north, the projection of the equator is taken as the abscissa axis Y, and the Y-axis is set to the east to form a coordinate system. Since the edge contour points extracted by the airborne and SLAM are both in this coordinate system and are two-dimensional points, the edge contour points of the crown point cloud obtained by the airborne point cloud are used as source feature points, and the edge contour points of the crown point cloud obtained by the SLAM point cloud are used as target feature points. According to the coordinates of the source feature points and the target feature points, the FPFH algorithm is used to calculate the relative position relationship of each feature point relative to other points in the neighborhood, and the position relationships of all relative neighborhood points are weightedly added according to the distance to obtain the position features of each feature point; its position features are given by Figure 2 It can be seen that they are: α, φ, θ. Therefore, each crown edge contour point of the airborne point cloud and the SLAM point cloud has its own unique position feature (α, φ, θ);

[0080] S5: The edge contour points of the crown point cloud obtained from the airborne point cloud are used as source feature points, and the edge contour points of the crown point cloud obtained from the SLAM point cloud are used as target feature points. According to the positional features of the source feature points and the target feature points, the size-constrained RANSAC algorithm (i.e., random sampling consistency algorithm) is used to preferentially match feature points with large crown sizes and the same positional features as points of the same name. The point cloud transformation matrix is ​​iteratively calculated based on the points of the same name, and the coordinates of each point in the source feature point are left-multiplied by the point cloud transformation matrix to obtain new coordinate information, thus completing the two-phase registration.

[0081] Example 2

[0082] An automatic registration method of airborne and SLAM point clouds for a forest plot is as described in Example 1, except that step S2 includes the following sub-steps:

[0083] S21: First, filter and denoise the SLAM point cloud data, and use the Cloth Simulation Filtering (CSF) algorithm to separate the ground points and non-ground points of the SLAM point cloud data to complete the preprocessing; delete the ground points and construct an octree structure for the non-ground points to establish the spatial relationship of the point cloud;

[0084] Octree structure process:

[0085] Determine the maximum and minimum values ​​on the three coordinate axes XYZ, build a minimum external bounding box containing all point cloud data as the root node, and then recursively divide the root node into eight sub-cubes of the same size, each of which becomes a child of the root node. For each child node, check whether it contains non-ground points. If so, continue to recursively split until the preset depth is reached, where the depth of the octree is calculated according to the formula:

[0086] D = ceil[log 2(L max -L min ) / r min )]

[0087] Where: L max is the maximum side length of the point cloud data; L min is the minimum side length; r min It is the minimum resolution of leaf nodes; ceil means rounding up.

[0088] S22: Use the random forest model to perform feature learning on the preprocessed SLAM point cloud data, and construct a multidimensional feature vector including geometric features and spatial features for the arc-shaped geometric features of the tree trunk point cloud. Use it as input data for the random forest model and train it to identify the tree trunk point cloud.

[0089] Specific:

[0090] The point cloud data in the node is sliced ​​longitudinally at equal intervals along the Z axis, and the RANSAC (Random Sample Consensus) algorithm is used to search for arcs on the sliced ​​projection point set {Pi (i = 1, 2, 3...)}. For the arc point set obtained by the search, feature parameters are used to describe its arc features, including geometric features: arc radius R i , fitting arc residual γ i , fitting interior point ratio I i etc.; Spatial characteristics: verticality of the center of the circle V i , the number of neighboring nodes C i And the root mean square error e i Etc. Construct a multidimensional feature vector F1 = [R i , I i , γ i , V i , e i , C i ] is used as input data for training the random forest model to realize the recognition and extraction of tree trunk point cloud.

[0091] S23: Based on the trunk point cloud recognition results, a complete crown point cloud is searched, and clustered to obtain a complete single tree point cloud. Specifically, the trunk point cloud is clustered using the Euclidean clustering method, and a column cluster area is constructed along the vertical direction of the pole (i.e., the height direction) based on the centroid of the cluster body. These areas are regarded as preliminary estimates of the trunk and its crown. The radius of the cluster area is set according to the diameter of the trunk; using the spatial proximity principle (the spatial proximity principle is based on Euclidean distance judgment, find the tree top node, and use the node as the initial calculation of the Euclidean distance between the surrounding eight direction nodes and the node. If it is less than 2 times the node side length, it is considered to be adjacent to the node; if there is a neighboring node in a certain direction, continue to search in the eight neighborhood directions of the neighboring node until there are no neighboring nodes), search for a complete crown point cloud in the eight neighborhood directions of the top of the trunk, project the crown point cloud to a two-dimensional plane, use the alpha-shape algorithm to extract edge contour points from the projection point set, and cluster the trunk point cloud and the crown point cloud into a single tree.

[0092] In the process of extracting the edge contour points of the crown point cloud from the SLAM point cloud, the octree structure is first constructed for the point cloud data to establish the spatial relationship of the point cloud; then, the parameters of the random forest are set, and the random forest model is trained using the eigenvector matrix; finally, the crown point cloud is searched by neighborhood clustering based on the trunk point cloud data predicted by the model, the crown point cloud is projected onto a two-dimensional plane, and the alpha-shape algorithm (existing algorithm) is used to extract the contour points to obtain the edge contour points of the crown.

[0093] Example 3

[0094] An automatic registration method of airborne and SLAM point clouds for a forest plot is as described in Example 1, except that step S3 includes the following sub-steps:

[0095] S31: cluster analysis is used to denoise the airborne point cloud data to remove noise points, such as flying points in the air and underground points, and ground filtering is performed through cloth simulation to obtain ground points PC_2 and non-ground points PC_3 (cloth simulation regards point cloud data as points on the cloth, and simulates the drooping process of the cloth under gravity by setting the hardness of the cloth. According to the simulation results, the ground point PC_2 at the bottom of the cloth after simulation and the non-ground point PC_3 above the cloth are identified);

[0096] Project the ground point PC_2 into a two-dimensional regular grid (a two-dimensional grid is a process of dividing a two-dimensional space into a series of rectangular grids of the same size, each grid being a grid unit. The process of projecting the ground point PC_2 into a two-dimensional regular grid is as follows: set the grid step size to 0.3, convert the Z coordinate of each ground point, project it from three-dimensional space to two-dimensional space, and determine the grid unit to which each point belongs based on the projected coordinates). Record the lowest point of the grid and interpolate the space grid to obtain a digital elevation model DEM (Digital Elevation Model). Similarly, project the non-ground point PC_3 into a two-dimensional regular grid, record the highest point of the grid to obtain a digital surface model DSM (Digital Surface Model), and make a difference between the two to obtain a canopy height surface model CHM (Canopy Height Model);

[0097] S32: normalize the canopy height surface model CHM to a grayscale space of 0-255, perform the first watershed segmentation using the foreground of the entire layer as a marker, and perform the second watershed segmentation using the upper layer of the segmented independent area as a marker;

[0098] S33: Taking the local highest point as the center and the median of the boundary detection as the crown radius constraint, the under-segmented area is segmented for the third time to obtain the edge contour points of the crown point cloud.

[0099] Multi-level watershed refers to the use of multi-level segmentation and marker control methods to effectively reduce the over-segmentation phenomenon that is prone to occur in traditional watershed methods. The first watershed uses the foreground of the entire layer as a marker to preliminarily divide different vegetation areas; the second watershed uses the upper layer of the independent area after segmentation as a marker to further subdivide the vegetation area; the third segmentation uses the local highest point as the center and combines the crown radius constraint to further refine the under-segmented area. Through a multi-level, gradually refined segmentation strategy, single trees can be accurately segmented. The single tree point cloud is clustered and projected onto a two-dimensional plane, and the alpha-shape algorithm is used to extract contour points to obtain the edge contour points of the crown.

[0100] Step S32 is specifically as follows:

[0101] The canopy height surface model CHM is normalized to a grayscale space of 0-255 to generate a grayscale image srcA. The normalization process is:

[0102] First, calculate the minimum value minVal and the maximum value maxVal in the canopy height surface model CHM, and then apply the linear transformation formula to each element in CHM: norVal = (val-minVal) * (255 / (maxVal-minVal)), where: norVal is the normalized value after transformation, ranging from 0 to 255; val is the current element value in CHM; minVal is the minimum value in the CHM matrix; maxVal is the maximum value in the CHM matrix. After the linear transformation, the minimum value minVal in CHM will be mapped to 0, the maximum value maxVal will be mapped to 255, and the other values ​​will be mapped to the corresponding positions between 0-255 in a linear proportion, and the processed CHM matrix will be converted into a grayscale image;

[0103] Binarize srcA to obtain a binary image srcB with pixel values ​​of 0 and 1, and perform morphological corrosion on the binary image srcB to remove small spots and noise in the binary image srcB for full-layer foreground marker extraction of the first watershed;

[0104] Cluster all foreground pixels with a value of 1 in the binary image srcB into multiple pixel clusters, and provide a mark index for each pixel cluster. The above information is recorded through the mark map marks;

[0105] The watershed segmentation method based on the marker map is used for segmentation and the grayscale map is used to optimize the segmentation results. The final segmentation results are recorded in marks. Based on the marks results, the area of ​​each region is calculated by counting the number of pixels in the region according to the formula gridsum. The maximum values ​​of the pixels in the region, gradmax and gradmin, are used to determine whether the pixel is located in the upper layer of the region.

[0106] The maximum value gradmax and the minimum value gradmin may not directly correspond to 1 and 0 in the binarized image. They refer to the upper and lower limits of the pixel values ​​retained or identified in the segmented area after a series of processing (such as morphological erosion, etc.). They are used to determine whether the pixel is located in the upper layer of the region. Specifically, when performing watershed segmentation of the upper foreground mark, by comparing the pixel value with these statistical values ​​(especially gradmax), combined with parameters such as the percentage of the top layer of the crown, it is determined which pixels belong to the upper layer of the crown, so as to perform segmentation more accurately.

[0107] When the area of ​​the region Area is greater than the threshold NormalCrown, the upper layer pixels are recorded in the binary map srcTreeC with a value of 1. When the area of ​​the region Area is less than the threshold NormalCrown but greater than the threshold MinCrown, the full layer pixels are recorded in the binary map srcTreeC with a value of 1, and the pixel values ​​of the remaining positions are 0. Binarization, erosion, clustering and the second watershed segmentation based on markers are performed again, and the segmentation results are recorded in marks2. In this embodiment, the threshold NormalCrown is set to 25, and the threshold MinCrown is set to 4.

[0108] Area=gridsum×Grid_Scale×Grid_Scale

[0109] Where gridsum is the number of pixels in the region; Grid_Scale is the length of the grid side.

[0110] Step S33 is specifically as follows:

[0111] Taking the tree bottom points extracted from the SLAM point cloud as the position constraint, the under-segmented area in marks2 is segmented for the third time, and the canopy diameters searched along four directions from the highest point are counted. The four directions refer to the east, south, west, and north (up, down, left, and right of the node) directions along the horizontal plane orthogonal to the highest point. 1 / 2 of the median is taken as the canopy radius constraint threshold radius (the threshold radius refers to arranging the searched canopy diameters from large to small, and selecting half of the median as the threshold radius). The pixel values ​​in the area whose distance from the highest point is greater than the threshold radius are set to -1; the point cloud in each segmented area recorded in marks2 is assigned a color as a separate tree, the single tree segmentation of the original point cloud is completed, and the edge contour points of the crown point cloud are obtained.

[0112] Example 4

[0113] An automatic registration method for airborne and SLAM point clouds of a forest plot is as described in Example 3, except that step S4 includes the following sub-steps:

[0114] S41: For the edge contour points of the tree crown point cloud obtained by the airborne point cloud and the SLAM point cloud, let P s Point is the point where the position feature needs to be calculated. t Point P s A point in the neighborhood of is defined as a local coordinate system (u, v, w) based on the normal vector of the two points and the vector of the line connecting the two points. The calculation formula of the three axes of the local coordinate system is as follows:

[0115] u=n s

[0116]

[0117] w=u×v

[0118] where n s For point P t The normal vector, P t -P s For point P s To P t The vector of ||P t -P s || 2 Indicates P t and P s The Euclidean distance between them is based on the local coordinate system, and the normal vector n is represented by a set of angles. s and n t The deviation between them is calculated as follows:

[0119] α=v·n t

[0120]

[0121] θ=arctan(w·n t ,u·n t )

[0122] Among them, α is P t The normal vector n t The angle between the coordinate axis v and the coordinate axis; φ is P s The normal vector n s The angle between the two points and the line; θ is the normal vector n t In u P t The angle between the projection of the w plane and the u direction;

[0123] Therefore, the surface geometric features of the object can be described by the three parameters (α, φ, θ), and the calculated surface geometric features can be put into the histogram and divided into 11 intervals to form a simplified point histogram (Simplified Point Feature Histogram) descriptor between the source point and a neighboring point, that is, the SPFH descriptor; the SPFH descriptors of all neighboring points within the neighborhood of the source point are calculated to obtain the feature descriptor of each neighboring point relative to the source point, and the formula is:

[0124]

[0125] S42: Take the weighted average of the SPFH descriptors of the i-th source point, and use w as the weight coefficient of its j-th neighbor point ij Represented by, the FPFH descriptor formula of the i-th point is:

[0126]

[0127] Among them, k is the number of points in the neighborhood of the i-th point;

[0128] like Figure 3 , red dot p q is any point in the feature point set, the yellow point p ki (1≤i≤5) is p q Points in the neighborhood, green point p ki (6≤i≤11) is p q Points outside the neighborhood.

[0129] The FPFH descriptors of the crown edge contour points of the extracted airborne point cloud and SLAM point cloud are calculated respectively to obtain the position features of each point, thereby providing the point position features required for the matching of the same-name point pairs for the RANSAC algorithm registration in step S5.

[0130] Example 5

[0131] An automatic registration method for airborne and SLAM point clouds of a forest plot is as described in Example 4, except that step S5 includes the following sub-steps:

[0132] S51: Calculate the minimum enclosing rectangle for each tree crown point cloud edge contour point separately, based on the predefined crown area threshold (the crown diameter of adult large and medium trees is about 6 to 12 meters, and the crown diameter of adult small trees is about 3 to 5 meters. The crown area should actually be the area of ​​the crown enclosing rectangle, so the crown area threshold can be taken as 25m 2 ) classify them into large crown sets and small crown sets;

[0133] S52: Estimate the point cloud transformation matrix T through the size-constrained RANSAC algorithm.

[0134] In step S51, the process of calculating the minimum enclosing rectangle is as follows:

[0135] A. Sort the edge contour points of the crown point cloud to obtain ordered edge contour points;

[0136] B. Based on the ordered edge contour points, for each edge contour point, select the next neighboring point adjacent to it to form a straight line, and use the straight line and the straight line perpendicular to it to determine the direction of the minimum enclosing rectangle, and calculate a minimum enclosing rectangle;

[0137] C. Calculate the vertical distance from the edge contour point to the four sides of the minimum bounding rectangle. If the vertical distance to the four sides is greater than the set threshold (the threshold is 12.5% ​​of the shortest side of the rectangle), it is a non-overlapping point, otherwise it is an overlapping point; the minimum bounding rectangle determined by any adjacent points is selected from the bounding rectangle with the largest overlap, and the ideal bounding rectangle can be obtained as the final minimum bounding rectangle;

[0138] D. Let the source feature point be P, the target feature point be Q, and each group of point clouds is classified according to the size of the crown area to form a set P = {P L ,P S} and Q = {Q L ,Q S}, where P L , Q L represents the large-scale crown set, P S , Q S Represents the set of edge points of small tree crowns.

[0139] The specific steps of step S52 are as follows:

[0140] a. Find P through the nearest neighbor matching of SPFH descriptor L The three random sampling points in Q L The corresponding points in Q L The corresponding points in are taken as the initial three pairs of points with the same name;

[0141] b. According to P L and Q L The three pairs of points with the same name {(p i ,q i )}, substituting its coordinates into the calculation can make P L Transform the same-name point in Q L In the middle, click p i After rotation and translation, it is transformed into point q i , through p i With q i The coordinates of the points can be used to calculate the rigid transformation relationship between the points, that is, the rotation matrix R and the translation vector t:

[0142] q i =Rp i +t

[0143] Through the rotation matrix R and the translation vector t, the matrix T of the assumed change of the same name point can be obtained:

[0144]

[0145] Among them, p i represents the coordinates of the i-th point in the source set of points with the same name, q i represents the coordinates of the i-th point in the target set of points with the same name, R is a rotation matrix used to represent the rotation part in the rigid transformation; t is a translation vector used to represent the translation part in the rigid transformation;

[0146] c. L The coordinates of each point in the left multiply the assumed change matrix T to obtain the transformed P LThe source feature point coordinates;

[0147] d. Through the Euclidean distance threshold (the Euclidean distance threshold can be 1m), at the converted source feature point P L and target feature point Q L Between, calculate Q L All points are the center of the circle, and the Euclidean distance threshold is the radius. The transformed P L In Q L The points within the range are regarded as inliers. If the number of inliers is less than 95% of the number of source feature points, return to step a and exclude the combination of these three points. Otherwise, the matrix is ​​considered to be a transformation matrix.

[0148] e. If the number of inliers after multiple iterations is always less than 95% of the number of source feature points, then P is introduced. S , Q S Supplement the points in and repeat steps ad;

[0149] f. Calculate the distance between the pairs of points with the same name. If the distance reaches the minimum value so far, it will be used as the final transformation matrix T.

[0150] The present invention adopts the size-constrained RANSAC algorithm, randomly and preferentially selects three points in the set with larger crown sizes from the source feature points, matches the corresponding three points with similar crown size information and the same position features in the target point cloud, calculates the transformation matrix, and transforms the source point cloud to the target point cloud benchmark through the matrix. The number of source point cloud points within the target point cloud distance threshold is counted, and it is continuously iterated until the number of points no longer increases, and then the points in the set with smaller crown sizes are introduced to participate in the matching of feature points, and the iteration is continued until the number of source point cloud points within the target point cloud distance threshold no longer increases.

[0151] The above is a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A method for automatic registration of airborne and SLAM point clouds of forest plots, characterized in that: The steps include: S1: Use airborne laser radar equipment to collect the tree tops in the forest area to obtain airborne point clouds; use SLAM equipment to collect the tree bottoms in the forest area to obtain SLAM point clouds; S2: Use the random forest model to classify and extract the tree point cloud from the SLAM point cloud, and extract the edge contour points of the crown point cloud; S3: Use a multi-level watershed algorithm to segment the airborne point cloud into single tree point clouds and extract edge contour points of the crown point cloud; S4: The edge contour points of the tree crown point cloud obtained by the airborne point cloud are used as source feature points, and the edge contour points of the tree crown point cloud obtained by the SLAM point cloud are used as target feature points. According to the coordinates of the source feature points and the target feature points, the FPFH algorithm is used to calculate the relative position relationship of each feature point relative to other points in the neighborhood range, and the position relationship of all relative neighborhood points is weightedly added according to the distance to obtain the position feature of each feature point; S5: The edge contour points of the crown point cloud obtained from the airborne point cloud are used as source feature points, and the edge contour points of the crown point cloud obtained from the SLAM point cloud are used as target feature points. According to the positional features of the source feature points and the target feature points, the size-constrained RANSAC algorithm is used to preferentially match feature points with large crown sizes and the same positional features as points of the same name. The point cloud transformation matrix is ​​iteratively calculated based on the points of the same name, and the coordinates of each point in the source feature point are left-multiplied by the point cloud transformation matrix to obtain new coordinate information to complete the registration.

2. The method for automatic registration of airborne and SLAM point clouds of forest plots according to claim 1, characterized in that: Step S2 includes the following sub-steps: S21: First, filter and denoise the SLAM point cloud data, and use the cloth simulation filtering algorithm to separate the ground points and non-ground points of the SLAM point cloud data to complete the preprocessing; delete the ground points and construct an octree structure for the non-ground points to establish the point cloud spatial relationship; S22: Use the random forest model to learn the features of the pre-processed SLAM point cloud data, and construct a multi-dimensional feature vector based on the arc-shaped geometric features of the tree trunk point cloud, which is used as the input data of the random forest model and trained to identify the tree trunk point cloud; S23: Based on the trunk point cloud recognition results, the complete crown point cloud is searched and clustered to obtain the complete single tree point cloud; for the single tree point cloud, it is projected onto a two-dimensional plane, and the contour points are extracted using the alpha-shape algorithm to obtain the edge contour points of the crown point cloud.

3. The method for automatic registration of airborne and SLAM point clouds of forest plots according to claim 2, characterized in that: Step S3 includes the following sub-steps: S31: De-noise the airborne point cloud data, perform ground filtering through cloth simulation to obtain ground points PC_2 and non-ground points PC_3; generate a digital surface model DSM and a digital terrain model DEM through interpolation, and obtain a canopy height surface model CHM by subtracting the two; S32: normalize the canopy height surface model CHM to a grayscale space of 0-255, perform the first watershed segmentation using the foreground of the entire layer as a marker, and perform the second watershed segmentation using the upper layer of the segmented independent area as a marker; S33: Taking the local highest point as the center and the median of the boundary detection as the crown radius constraint, the under-segmented area is segmented for the third time to obtain the edge contour points of the crown point cloud.

4. The method for automatic registration of airborne and SLAM point clouds of forest plots according to claim 3, characterized in that: Step S32 is specifically as follows: The canopy height surface model CHM is normalized to a grayscale space of 0-255 to generate a grayscale image srcA, srcA is binarized to obtain a binary image srcB, with pixel values ​​of 0 and 1, and the binary image srcB is morphologically eroded to remove noise for full-layer foreground marker extraction in the first watershed; Cluster all foreground pixels with a value of 1 in the binary image srcB into multiple pixel clusters, and provide a mark index for each pixel cluster. The above information is recorded through the mark map marks; Use the watershed segmentation method based on the marker image to perform segmentation and optimize the segmentation result through the grayscale image. The final segmentation result is recorded in marks. Based on the marks result, the area of ​​each region is calculated by counting the number of regional pixels gridsum according to the formula, and the maximum values ​​of regional pixels gradmax and gradmin are counted to determine whether the pixel is located in the upper layer of the region; when the area of ​​the region Area is greater than the threshold NormalCrown, the upper layer pixels are recorded in the binary map srcTreeC with a value of 1; when the area of ​​the region Area is less than the threshold NormalCrown but greater than the threshold MinCrown, the pixels of the entire layer are recorded in the binary map srcTreeC with a value of 1, and the pixel values ​​of the remaining positions are 0; binarization, erosion, clustering and the second watershed segmentation based on markers are performed again, and the segmentation results are recorded in marks2; Step S33 is specifically as follows: Taking the tree bottom points extracted from the SLAM point cloud as the position constraint, the under-segmented area in marks2 is segmented for the third time. The canopy diameters searched along the four directions of the highest point are counted, and 1 / 2 of the median is taken as the canopy radius constraint threshold radius. The pixel values ​​in the area whose distance from the highest point is greater than the threshold radius are set to -1; the point cloud in each segmented area recorded in marks2 is assigned a color as a separate tree, completing the single tree segmentation of the original point cloud and obtaining the edge contour points of the crown point cloud.

5. The method for automatic registration of airborne and SLAM point clouds of forest plots according to claim 3, characterized in that: Step S4 includes the following sub-steps: S41: For the edge contour points of the tree crown point cloud obtained by the airborne point cloud and the SLAM point cloud, let P s Point is the point where the position feature needs to be calculated. t Point P s A point in the neighborhood of is defined as a local coordinate system (u, v, w) based on the normal vector of the two points and the vector of the line connecting the two points. The calculation formula of the three axes of the local coordinate system is as follows: u=n s w=u×v where n s Point P t The normal vector, P t -P s Point P s To P t The vector of ||P t -P s ||2 means P t and P s The Euclidean distance between them is based on the local coordinate system, and the normal vector n is represented by a set of angles. s and n t The deviation between them is calculated as follows: α=v·n t θ=arctan(w·n t ,u·n t ) Among them, α is P t The normal vector n t The angle between the coordinate axis v and P s The normal vector n s The angle between the two points and the line; θ is the normal vector n t In u P t The angle between the projection of the w plane and the u direction; Therefore, the three parameters (α, φ, θ) are used to describe the surface geometric features of the object. The calculated surface geometric features are put into the histogram and divided into 11 intervals to form a simple point histogram descriptor between the source point and a neighboring point, namely the SPFH descriptor; the SPFH descriptors of all neighboring points within the neighborhood of the source point are calculated to obtain the feature descriptor of each neighboring point relative to the source point, and the formula is: S42: Take the weighted average of the SPFH descriptors of the i-th source point, and use w as the weight coefficient of its j-th neighbor point ij Represented by, the FPFH descriptor formula of the i-th point is: Among them, k is the number of points in the neighborhood of the i-th point; The FPFH descriptors of the tree crown edge contour points of the extracted airborne point cloud and SLAM point cloud are calculated respectively to obtain the position features of each point.

6. The method for automatic registration of airborne and SLAM point clouds of forest plots according to claim 5, characterized in that: Step S5 includes the following sub-steps: S51: Calculate the minimum enclosing rectangle for the edge contour points of the crown point cloud of each tree separately, and classify them into a large crown set and a small crown set based on a predefined crown area threshold; S52: Estimate the point cloud transformation matrix T through the size-constrained RANSAC algorithm.

7. The method for automatic registration of airborne and SLAM point clouds of forest plots according to claim 6, characterized in that: In step S51, the process of calculating the minimum enclosing rectangle is as follows: A. Sort the edge contour points of the crown point cloud to obtain ordered edge contour points; B. Based on the ordered edge contour points, for each edge contour point, select the next neighboring point adjacent to it to form a straight line, and use the straight line and the straight line perpendicular to it to determine the direction of the minimum enclosing rectangle, and calculate a minimum enclosing rectangle; C. Calculate the vertical distance from the edge contour point to the four sides of the minimum bounding rectangle. If the vertical distance to the four sides is greater than the set threshold, it is a non-overlapping point, otherwise it is an overlapping point. The minimum bounding rectangle determined by any adjacent points is selected from the bounding rectangle with the largest overlap, and the ideal bounding rectangle can be obtained as the final minimum bounding rectangle. D. Let the source feature point be P, the target feature point be Q, and each group of point clouds is classified according to the size of the crown area to form a set P = {P L ,P S } and Q = {Q L ,Q S }, where P L , Q L represents the large-scale crown set, P S , Q S Represents the set of edge points of small tree crowns.

8. The method for automatic registration of airborne and SLAM point clouds of forest plots according to claim 7, characterized in that: The specific steps of step S52 are as follows: a. Find P through the nearest neighbor matching of SPFH descriptor L The three random sampling points in Q L The corresponding points in Q L The corresponding points in are taken as the initial three pairs of points with the same name; b. According to P L and Q L The three pairs of points with the same name {(p i ,q i )}, substituting its coordinates into the calculation can make P L Transform the same-name point in Q L In the middle, click p i After rotation and translation, it is transformed into point q i , through p i With q i The coordinates of the points can be used to calculate the rigid transformation relationship between the points, that is, the rotation matrix R and the translation vector t: q i =Rp i +t Through the rotation matrix R and the translation vector t, the matrix T of the assumed change of the same name point can be obtained: Among them, p i represents the coordinates of the i-th point in the source set of points with the same name, q i represents the coordinates of the i-th point in the target set of points with the same name, R is a rotation matrix used to represent the rotation part in the rigid transformation; t is a translation vector used to represent the translation part in the rigid transformation; c. L The coordinates of each point in the left multiply the assumed change matrix T to obtain the transformed P L The source feature point coordinates; d. Through the Euclidean distance threshold, the converted source feature point P L and target feature point Q L Between, calculate Q L All points are the center of the circle, and the Euclidean distance threshold is the radius. The transformed P L In Q L The points within the range are regarded as inliers. If the number of inliers is less than 95% of the number of source feature points, return to step a and exclude the combination of these three points. Otherwise, the matrix is ​​considered to be a transformation matrix. e. If the number of inliers after multiple iterations is always less than 95% of the number of source feature points, then P is introduced. S , Q S Supplement the points in and repeat steps ad; f. Calculate the distance between the pairs of points with the same name. If the distance reaches the minimum value so far, it will be used as the final transformation matrix T.

Images