Single wood volume nondestructive measurement method based on foundation laser radar

Abstract

The invention provides a single wood volume nondestructive measurement method based on a foundation laser radar, belongs to the field of single wood volume, and focuses on feature extraction of trunk structures in different forms and calculation of the trunk volume. According to the method, masson pine is taken as a research object, 52 individual trees are taken as samples, true values of different height diameters are acquired, the volume of the individual trees is calculated, and meanwhile, a high-density three-dimensional laser point cloud is scanned by using a foundation laser radar system. The method comprises the following steps: firstly, sampling original point cloud data, constraining growth of seed points by using normal vectors and angles, and accurately extracting trunk point clouds; simulating interception of a real disc through point cloud slices, and calculating a single-layer point cloud area by using an improved convex hull algorithm; in the shielded region, constructing a fitting circle by using the local point cloud to calculate the area of the layer; and finally, the single wood volume is obtained through the thought of an infinitesimal method. The correlation R2 between the calculated wood volume and the real calculation reaches 0.995, the RMSE is 0.017 m < 3 >, and the nRMSE is 2.38%, which indicates that the tree trunk can be accurately extracted and the single wood volume can be calculated under the non-contact condition.

Description

Technical Field

[0001] The present invention relates to the field of tree trunk algorithms, and in particular to a non-destructive measurement method for the volume of a single tree based on a ground-based laser radar. Background Art

[0002] With the growing demand for refined management of forest resources, accurate and efficient extraction of individual tree trunk information has become a key link in forestry digitization, especially in areas dominated by economic forests. The accurate extraction of trunk parameters (such as diameter at breast height and volume) is of great significance for wood quality assessment, biomass estimation, felling strategy formulation, and forest health monitoring. Traditionally, the acquisition of trunk parameters mainly relies on manual measurement, which is a highly labor-intensive and inefficient process. Some key structural attributes, such as the change in trunk diameter along the height, trunk curvature, and cross-sectional profile, often need to be measured section by section or obtained using destructive sampling methods, which is not only time-consuming and labor-intensive, but also causes irreversible damage to the sample site structure. Under conditions of tall trees or dense canopies, problems such as partial obstruction of trunks and poor accessibility also limit the feasibility and accuracy of measurements.

[0003] LiDAR technology is widely used in tree parameter extraction because it can use echoes to obtain high-precision point clouds, thereby completing the three-dimensional reconstruction of the target object. Ground-based LiDAR systems (TLS) are increasingly used in forest trunk modeling and parameter estimation due to their high spatial resolution. However, the complex understory environment, the irregular spatial distribution of trees, the interlaced occlusion of stems, and the noise problem in point cloud data all pose significant challenges to the automatic extraction of individual tree trunks. At the same time, there is a lack of more sophisticated trunk surface reconstruction methods, and the use of cylindrical fitting to represent the growth of tree trunks cannot meet the requirements for accurate calculation of the volume of individual trees.

[0004] The present invention provides a non-destructive measurement method for single-tree volume based on ground-based laser radar, focusing on extracting structural features of tree trunks of different shapes and calculating trunk volume. The method uses Pinus massoniana as the research object and 52 single trees as samples. The real values ​​of diameters at different heights are collected and the volume of single trees is calculated. At the same time, a high-density three-dimensional laser point cloud is scanned using a ground-based laser radar system. First, the original point cloud data is sampled and processed, and the growth of seed points is constrained by normal vectors and angles to accurately extract the trunk point cloud. Then, the point cloud is sliced ​​to simulate the interception of a real disk, and the area of ​​a single layer of point cloud is calculated using an improved convex hull algorithm. In the occluded area, the local point cloud is used to construct a fitted circle to calculate the area of ​​that layer. Finally, the volume of single trees is obtained using the idea of ​​the differential element method. Summary of the Invention

[0005] This patent provides a nondestructive method for measuring the volume of a single tree using ground-based LiDAR. Based on a two-dimensional convex hull algorithm, this patent proposes a dual-constraint convex hull algorithm. This algorithm optimizes the convex hull boundary to ensure a close fit between the boundary and the tree trunk point cloud. A distance threshold is also added to improve computational efficiency while maintaining accuracy, thereby enhancing the model's accurate description of the tree trunk shape.

[0006] The technical solution of the present invention provides: a non-destructive measurement method of a single wood volume based on a ground-based laser radar, comprising the following steps:

[0007] S101: Data Collection

[0008] The main tree species studied was Pinus massoniana. The diameter at breast height, tree height, diameters at different relative positions, height below branches, and crown width of 153 Pinus massoniana trees were collected and measured. The bark volume of each tree was calculated based on the diameters at different relative positions. A FARO Focus3D terrestrial 3D laser scanner was used to scan the Pinus massoniana to obtain high-precision 3D point clouds, and multi-site clouds were stitched together using a public target sphere.

[0009] S102: Data Preprocessing

[0010] The average density of the data point cloud collected by the ground-based lidar is as high as 5000p / m 2 In order to avoid the reduction of computational efficiency caused by the large amount of data from ground-based lidar, the voxel downsampling method is used to process the original ground-based lidar. This method calculates the centroid of each voxel unit and replaces all the points in the unit with the centroid coordinates to complete the point cloud downsampling under clearing conditions while ensuring that the distribution characteristics of the point cloud remain unchanged.

[0011] Based on the distribution of points, this patent sets the voxel side length r to 0.1m;

[0012] Through morphological comparison and screening, 52 typical trees with different morphologies were finally selected as experimental subjects;

[0013] S103: Trunk point cloud extraction

[0014] In order to achieve non-destructive measurement of the volume of a single tree, it is still necessary to accurately identify the point cloud of the tree trunk from the scattered and disordered laser point cloud. This patent proposes a tree trunk extraction algorithm based on normal vector constraints. The specific steps of the algorithm are as follows:

[0015] Step 1: Initialize all point sets Points, seed point set Seed Points, adjacent point set AdjacentPoints, trunk point set TrunkPoints, and temporary point set Temp Points;

[0016] Step 2: Read the pre-processed point cloud data of a single tree and establish a normal vector for the point cloud data; for each scan point p in the point cloud, search for its K nearest neighboring points, and then calculate the local plane P of these points in the least squares sense. This plane can be expressed as:

[0017]

[0018] Where, is the normal vector of plane P, and d is the distance from P to the origin of the coordinate system;

[0019] It can be considered that the normal vector of the plane fitted by the K nearest neighbor points is the normal vector of the current scanning point; the normal vector of the plane P can be obtained by principal component analysis (PCA). It is known from the operation that P passes through the centroid P0 of its K neighbor points, and the normal vector satisfy First, perform eigenvalue decomposition on the covariance matrix M in formula (2) to obtain the eigenvalues ​​of M. The eigenvector corresponding to the minimum eigenvalue of M is the normal vector of P;

[0020]

[0021] The point cloud normal vector represents the direction of the plane formed by its neighboring points. This patent uses a 1m search radius and 10 nearest neighbor points to calculate the point cloud normal vector after establishing a kd-tree for the entire point cloud.

[0022] Step 3: Based on the establishment of the normal vector, initialize the Points point set, assign each point a normal vector attribute value and set the category of all points to unclassified;

[0023] Step 4: When When Z n Can be and The cross product of Z n Indicates the verticality of the point cloud normal vector. When Z n When it approaches 0, it means that the normal vector is almost perpendicular to the Z axis; therefore, select Z n The area with a value greater than 0.05 and a point cloud height greater than 0.05m and less than 0.2m is used as the initial seed point for subsequent growth, and the points that meet the conditions are added to SeedPoints;

[0024] Step 5: To further improve algorithm efficiency and reduce the amount of data traversed by the loop, the algorithm adopts a segmented growth method and adds the seed points of the previous segment as the ignition seed points of the new segment. In this patent, a 0.3m interval and a 0.3m layer thickness are used to layer the single tree point cloud, and the layered point cloud is added to the TempPoints collection. The pseudo code for point cloud slicing and layering is as follows:

[0025]

[0026] Step 6: For each layer of temporary points, loop and add each seed point, and calculate the Euclidean distance of the seed point to the 10 nearest neighboring points, and determine whether the angle between each point and the normal vector of the seed point is less than the angle threshold A. T ; Change the category of the points that meet the conditions to trunk points, and at the same time change the points that meet the conditions to seed points Z n The difference is less than Z T The point is recorded as the new seed point;

[0027]

[0028]

[0029] Step 7: Loop through each slice point cloud and gradually add the points that meet the conditions to the trunk point cloud set to generate the final trunk point cloud;

[0030] S104: Double Constrained Convex Hull Algorithm

[0031] Based on the Graham convex hull algorithm, a double-constrained convex hull algorithm is proposed. This algorithm mainly improves the Graham convex hull algorithm. The tree trunk point cloud is sliced ​​in the same way as the slicing algorithm in S103 and input into the optimization algorithm one by one to accurately extract the boundary points and calculate the single-layer area. The algorithm mainly uses the angle and CHL convex hull side length threshold control to gradually approximate the actual point cloud distribution, making the convex hull more closely fit the tree trunk surface. Then, the discretized Green's formula is used to calculate the convex hull area of ​​each layer.

[0032] Sort the generated convex hull in counterclockwise direction. At this time, for the S vertices of the convex set (x1,y1),(x2,y2),...,(x n ,y n ), the area calculation formula is (Formula (3)):

[0033]

[0034] Where, Area is the convex hull area calculated by the optimized convex hull algorithm, (x i ,y i ) The coordinates of the convex hull vertices arranged in counterclockwise order, where n is the number of convex hull points; the initial point is usually repeatedly added to the end of the point set to close the polygon;

[0035] S105: Hough circle fitting algorithm

[0036] Due to the scanning method of ground-based lidar data, areas above an average height of 15 meters are difficult to capture, which easily leads to missing point cloud data and incomplete tree trunk point clouds. The double-constrained convex hull algorithm has difficulty completing the point cloud when data is missing, and the calculated single-layer area is too small. Hough circle fitting is used to fit the shape of the missing points at higher locations and calculate the area to compensate for the inability to accurately calculate the upper layer area. The specific principles of the Hough circle fitting algorithm are as follows:

[0037] Traverse each edge point (x i ,y i ): Calculate the possible center of the circle for all possible radii (a k ,b k ), and corresponding to the accumulator A[a k ][b k ][r k ] plus 1; search for local maxima exceeding the threshold in the accumulator, and these maxima are the parameters (a*, b*, r*) of the detected circle.

[0038] S106: Calculation of volume of single timber

[0039] The volume of a single tree is calculated using the idea of ​​the infinitesimal method. The trunk of each slice layer is assumed to be a number of frustum-like cones. The frustum calculation formula is used to calculate the volume of the trunk of each layer (Formula (4)):

[0040]

[0041] Because the improved convex hull algorithm DB is more sensitive to missing point cloud data, circle fitting is used to calculate the area when the point cloud data is incomplete. To further realize the automatic calculation of material volume, the algorithm designs an area ratio threshold ARatio to automatically determine whether the area ratio of the DB algorithm and the circle fitting is greater than ARatio. If it is greater than the threshold, it indicates that the DB algorithm is greatly affected by the missing data, and the area of ​​this layer is calculated using the result of circle fitting.

[0042] To further correct the area results of each layer and prevent the occurrence of larger or smaller values ​​of fitting, the algorithm uses a standard deviation threshold to assist in constraining the area of ​​each slice layer. The reliability of each layer area is verified by comparing the area of ​​the current layer with the area of ​​the previous layers. A decrease in area usually indicates that the normal stem is gradually becoming thinner, while a larger mutation indicates that there is an incorrect fitting. In the case of an area mutation, the algorithm will discard the data of the current layer to prevent incorrect calculation (Formula (5)):

[0043] Diff>K*S (5)

[0044] where Diff={a i -a j |a i ,aj ∈P,i≠j}, represents the difference between the area of ​​the current layer and the average area of ​​the previous 10 layers; K represents the control multiple of the threshold; S represents the standard deviation of Diff;

[0045] Although the ground-based lidar system can collect high-density point clouds, its ability to collect information about objects above 15 meters is weak, and it cannot capture relatively complete point cloud information, resulting in the inability to accurately identify tree trunks at high altitudes. The algorithm uses cone fitting to calculate the volume of tree trunks at high altitudes (Formula (6)):

[0046]

[0047] Where S n represents the area of ​​the highest point cloud detected; H represents the height of a single tree; slice represents the number of slice layers; Δh represents the slice layer height;

[0048] The volume of a single log is the sum of the volume of the logs in the following layers and the volume of the fitted cone in the top layer, that is,

[0049] V t =V i +V j (7).

[0050] Beneficial effects of the present invention:

[0051] This application is based on TLS (Terrain LiDAR Systems) data and proposes a non-destructive measurement method for single wood volume based on ground-based laser radar, focusing on the extraction of trunk structural features of different forms and the calculation of trunk volume. This method mainly takes Pinus massoniana as the research object and 52 single trees as samples. It destructively collects the true values ​​of diameters at different heights and calculates the volume of single wood. At the same time, it uses a ground-based laser radar system to scan high-density three-dimensional laser point clouds. First, the original point cloud data is sampled and processed, and the growth of seed points is constrained by normal vectors and angles to accurately extract the trunk point cloud; then, the point cloud is sliced ​​to simulate the interception of a real disk, and the improved convex hull algorithm is used to calculate the area of ​​a single layer of point cloud; in the occluded area, the local point cloud is used to construct a fitting circle to calculate the area of ​​the layer; finally, the volume of single wood is obtained through the idea of ​​the differential element method. The correlation R between the volume calculated by the present invention and the actual calculated volume 2 Up to 0.995, RMSE is 0.017m 3 , nRMSE is 2.38%, indicating that the invention can accurately extract tree trunks and calculate the volume of single wood without contact, and the results are reliable. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] To more clearly illustrate the technical solutions of the present disclosure, the following briefly introduces the drawings required for use in some embodiments of the present disclosure. Obviously, the drawings described below are only drawings of some embodiments of the present disclosure, and those skilled in the art can also derive other drawings based on these drawings. Furthermore, the drawings described below are schematic diagrams and are not intended to limit the actual dimensions of the products, actual processes of the methods, actual timing of signals, and the like involved in the embodiments of the present disclosure.

[0053] Figure 1 Schematic diagram of the proposed algorithm flow; (a) data preprocessing; (b) trunk extraction; (c) single-layer slice diameter and area fitting; (d) single wood volume calculation.

[0054] Figure 2 The figure shows the result of accurately identifying the tree trunk point cloud using this algorithm. The red points represent the tree trunk.

[0055] Figure 3 To improve the accuracy of volume calculation using this algorithm, the black dashed line represents the 1:1 line. DETAILED DESCRIPTION

[0056] In order to enable those skilled in the art to better understand the technical solutions in this application, the following will provide a clear and complete description of the technical solutions in the embodiments of this application in conjunction with the drawings in the embodiments of this application. Obviously, the embodiments described are only part of the embodiments of this application, not all of the embodiments. Based on the embodiments in this application, all other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of this application.

[0057] The following will be combined with the drawings in the embodiments of this application to clearly and completely describe the technical solutions in the embodiments of this application. Obviously, the embodiments described are only part of the embodiments of this application, not all of the embodiments. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without making creative efforts are within the scope of protection of this application.

[0058] The present invention is described in further detail below by way of examples.

[0059] In one embodiment, a method for non-destructive measurement of a single wood volume based on a ground-based laser radar comprises the following steps:

[0060] Data collection

[0061] The main tree species studied was Pinus massoniana. The diameter at breast height, tree height, diameters at different relative positions, height below branches, and crown width of 153 Pinus massoniana trees were collected and measured. The bark volume of each tree was calculated based on the diameters at different relative positions. A FARO Focus3D terrestrial 3D laser scanner was used to scan the Pinus massoniana to obtain a high-precision 3D point cloud. Multi-site cloud splicing was completed based on a public target sphere. The average point cloud density was 4410p / m 2 .

[0062] Data preprocessing

[0063] The average density of the data point cloud collected by the ground-based lidar is as high as 5000p / m 2 In order to avoid the reduction of computational efficiency caused by the large amount of data of ground-based lidar, the original ground-based lidar is processed by voxel downsampling. This method calculates the center of gravity of each voxel unit and replaces all the points in the unit with the coordinates of the center of gravity to complete the point cloud downsampling under clearing to ensure that the distribution characteristics of the point cloud remain unchanged.

[0064] In this patent, based on the distribution of points, the voxel side length r is set to 0.1m.

[0065] Through morphological comparison and screening, 52 typical trees with different morphologies were finally selected as experimental subjects.

[0066] Trunk point cloud extraction

[0067] In order to achieve non-destructive measurement of the volume of a single tree, it is still necessary to accurately identify the point cloud of the tree trunk from the scattered and disordered laser point cloud. This patent proposes a tree trunk extraction algorithm based on normal vector constraints. The specific steps of the algorithm are as follows:

[0068] Step 1: Initialize all point sets Points, seed point set Seed Points, adjacent point set AdjacentPoints, trunk point set TrunkPoints, and temporary point set Temp Points;

[0069] Step 2: Read the pre-processed point cloud data of a single tree and establish a normal vector for the point cloud data. For each scan point p in the point cloud, search for its K nearest neighboring points, and then calculate the local plane P of these points in the least squares sense. This plane can be expressed as:

[0070]

[0071] Where, is the normal vector of plane P, and d is the distance from P to the origin of the coordinate system;

[0072] It can be considered that the normal vector of the plane fitted by the K nearest neighbor points is the normal vector of the current scanning point; the normal vector of the plane P can be obtained by principal component analysis (PCA). It is known from the operation that P passes through the centroid P0 of its K neighbor points, and the normal vector satisfy First, perform eigenvalue decomposition on the covariance matrix M in formula (2) to obtain the eigenvalues ​​of M. The eigenvector corresponding to the minimum eigenvalue of M is the normal vector of P.

[0073]

[0074] The normal vector of a point cloud represents the direction of the plane formed by its neighboring points. This patent establishes a kd-tree for all point clouds and uses a 1m search radius and 10 nearest neighbor points to calculate the point cloud normal vector.

[0075] Step 3: Based on the establishment of the normal vector, initialize the Points point set, assign each point a normal vector attribute value and set the category of all points to unclassified;

[0076] Step 4: When When Z n Can be and The cross product of Z n Indicates the verticality of the point cloud normal vector. When Z n When it approaches 0, it means that the normal vector is almost perpendicular to the Z axis. Therefore, select Z n The area with a value greater than 0.05 and a point cloud height greater than 0.05m and less than 0.2m is used as the initial seed point for subsequent growth, and the points that meet the conditions are added to SeedPoints.

[0077] Step 5: To further improve the efficiency of the algorithm and reduce the amount of data traversed by the loop, the algorithm adopts a segmented growth method and adds the seed points of the previous segment as the ignition seed points of the new segment. This patent uses a 0.3m interval and a 0.3m layer thickness to layer the single tree point cloud and add the layered point cloud to the TempPoints collection. The pseudo code for point cloud slicing and layering is as follows:

[0078]

[0079]

[0080] Step 6: For each layer of temporary points, loop and add each seed point, and calculate the Euclidean distance of the seed point to the 10 nearest neighboring points, and determine whether the angle between each point and the normal vector of the seed point is less than the angle threshold A. T ; Change the category of the points that meet the conditions to trunk points, and at the same time change the points that meet the conditions to seed points Z n The difference is less than Z TThe point is recorded as the new seed point.

[0081]

[0082] Step 7: Loop through each slice point cloud and gradually add the points that meet the conditions to the trunk point cloud set to generate the final trunk point cloud.

[0083] Double-constrained convex hull algorithm

[0084] Based on the Graham convex hull algorithm, a dual-constrained convex hull algorithm is proposed, which primarily improves the Graham convex hull algorithm. Using the same slicing method as in S103, the tree trunk point cloud is sliced ​​and fed into the optimization algorithm one by one, accurately extracting boundary points and calculating the area of ​​each layer. This algorithm primarily uses angle and CHL convex hull edge length thresholds to gradually approximate the actual point cloud distribution, ensuring that the convex hull more closely matches the tree trunk surface. The discretized Green's formula is then used to calculate the area of ​​each layer's convex hull.

[0085] Sort the generated convex hull in counterclockwise direction. At this time, for the S vertices of the convex set (x1,y1),(x2,y2),...,(x n ,y n ), the area calculation formula is (Formula (3)):

[0086]

[0087] Where, Area is the convex hull area calculated by the optimized convex hull algorithm, (x i ,y i ) The coordinates of the convex hull vertices are arranged in counterclockwise order, and n is the number of convex hull points; the initial point is usually repeatedly added to the end of the point set to close the polygon.

[0088] Hough circle fitting algorithm

[0089] Due to the scanning method of ground-based LiDAR data, areas above an average height of 15m are difficult to capture, which easily leads to missing point cloud data and incomplete tree trunk point clouds. The double-constrained convex hull algorithm is difficult to complete the point cloud when data is missing, and the calculated single-layer area is too small. This patent uses Hough circle fitting to fit the shape of the missing point cloud at a higher position and calculate the area, thus compensating for the inability to accurately calculate the upper layer area. The specific principles are as follows:

[0090] Traverse each edge point (x i ,y i ): Calculate the possible center of the circle for all possible radii (a k ,b k ), and corresponding to the accumulator A[a k ][bk ][r k ] plus 1; search for local maxima exceeding the threshold in the accumulator, and these maxima are the parameters (a*, b*, r*) of the detected circle.

[0091] Calculation of volume of single timber

[0092] The volume of a single tree is calculated using the idea of ​​the infinitesimal method. The trunk of each slice layer is assumed to be a number of frustum-like cones. The frustum calculation formula is used to calculate the volume of the trunk of each layer (Formula (4)):

[0093]

[0094] Because the improved convex hull algorithm DB is more sensitive to missing point cloud data, circle fitting is used to calculate the area when the point cloud data is incomplete. To further realize the automated calculation of material volume, the algorithm designs an area ratio threshold ARatio to automatically determine whether the area ratio of the DB algorithm and the circle fitting is greater than ARatio. When it is greater than the threshold, it indicates that the DB algorithm is greatly affected by the missing data, and the area of ​​this layer uses the result of circle fitting.

[0095] To further correct the area results of each layer and prevent the occurrence of larger or smaller values ​​from being fitted, the algorithm uses a standard deviation threshold to assist in constraining the area of ​​each slice layer. The reliability of each layer's area is verified by comparing the area of ​​the current layer with the areas of the previous layers. A decrease in area usually indicates that a normal stem is gradually becoming thinner, while a larger mutation indicates that there is an incorrect fitting. In the case of a mutation in area, the algorithm will discard the data of the current layer to prevent incorrect calculation (Formula (5)):

[0096] Diff>K*S (5)

[0097] where Diff={a i -a j |a i ,a j ∈P,i≠j}, represents the difference between the area of ​​the current layer and the average area of ​​the previous 10 layers; K represents the control multiple of the threshold; S represents the standard deviation of Diff.

[0098] Although the ground-based lidar system can collect high-density point clouds, its ability to collect information about objects above 15 meters is weak, and it cannot capture relatively complete point cloud information, resulting in the inability to accurately identify tree trunks at high altitudes. The algorithm uses cone fitting to calculate the volume of tree trunks at high altitudes (Formula (6)):

[0099]

[0100] Where S nIt represents the area of ​​the highest point cloud detected; H represents the height of a single tree; slice represents the number of slice layers; Δh represents the slice layer height.

[0101] The volume of a single log is the sum of the volume of the logs in the following layers and the volume of the fitted cone in the top layer, that is,

[0102] V t =V i +V j (7)

[0103] Algorithm performance evaluation

[0104] Evaluation of volume calculation accuracy

[0105] The volume of 52 individual trees surveyed in the field was used as the actual value to compare with the volume calculated by the algorithm proposed in this invention. 2 ) and the root mean square error (RMSE) were used to evaluate the results (Formula (8-10)).

[0106]

[0107] Among them, Y i represents the true value of the field survey; represents the timber volume calculated by the algorithm proposed in the present invention; n represents the number of trees; and Y max Y min The mean, maximum, and minimum values ​​of Y i .

[0108] Volume calculation results

[0109] The extracted volume is compared with the true value ( Figure 3 ), R 2 is 0.995, and the RMSE is 0.017m 3 , with an nRMSE of 2.38%. This evaluation shows that the error between the volume calculated by the proposed algorithm and the actual destructive measurement of individual timber volume is minimal, demonstrating that the algorithm has extremely high accuracy in calculating individual timber volume. This low nRMSE value fully demonstrates the algorithm's high stability and reliability, and the calculated timber volume meets the practical needs of high-precision forestry resource surveys and monitoring.

[0110] The above describes the basic principles and main features of the present invention and the advantages of the present invention. It is obvious to those skilled in the art that the present invention is not limited to the details of the above exemplary embodiments, and that the present invention can be implemented in other specific forms without departing from the spirit or essential characteristics of the present invention. Therefore, from all points of view, the embodiments should be regarded as exemplary and non-restrictive. The scope of the present invention is defined by the appended claims rather than the above description, and it is intended that all changes that fall within the meaning and range of equivalents of the claims are included in the present invention. Any reference signs in the claims should not be construed as limiting the claim to which they relate.

[0111] In addition, it should be understood that although this specification is described in terms of implementation methods, not every implementation method contains only one independent technical solution. This narrative method of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment can also be appropriately combined to form other implementation methods that can be understood by those skilled in the art.

Claims

1. A non-destructive measurement method for single wood volume based on ground-based laser radar, characterized by: The following steps are involved: S101: Data Collection The main tree species studied was Pinus massoniana. The diameter at breast height, tree height, diameters at different relative positions, height below branches, and crown width of 153 Pinus massoniana trees were collected and measured. The bark volume of each tree was calculated based on the diameters at different relative positions. A FARO Focus 3D terrestrial 3D laser scanner was used to scan the Pinus massoniana to obtain a high-precision 3D point cloud. The cloud was then spliced ​​from multiple sites using a public target sphere. S102: Data Preprocessing The voxel downsampling method processes the original ground-based lidar. This method calculates the center of gravity of each voxel unit and replaces all the points in the unit with the center of gravity coordinates to complete the point cloud downsampling in a clear manner while ensuring that the point cloud distribution characteristics remain unchanged. Based on the distribution of points, this patent sets the voxel side length r to 0.1m; Through morphological comparison and screening, 52 typical trees with different morphologies were finally selected as experimental subjects; S103: Trunk point cloud extraction A tree trunk extraction algorithm based on normal vector constraints is proposed. The specific steps of the algorithm are as follows: Step 1: Initialize all point sets Points, seed point set Seed Points, adjacent point set AdjacentPoints, trunk point set TrunkPoints, and temporary point set Temp Points; Step 2: Read the pre-processed point cloud data of a single tree and establish a normal vector for the point cloud data; for each scan point p in the point cloud, search for its K nearest neighboring points, and then calculate the local plane P of these points in the least squares sense. This plane can be expressed as: Where, is the normal vector of plane P, and d is the distance from P to the origin of the coordinate system; It can be considered that the normal vector of the plane fitted by the K nearest neighbor points is the normal vector of the current scanning point; the normal vector of the plane P can be obtained by principal component analysis (PCA). It is known from the operation that P passes through the centroid P0 of its K neighbor points, and the normal vector satisfy First, perform eigenvalue decomposition on the covariance matrix M in formula (2) to obtain the eigenvalues ​​of M. The eigenvector corresponding to the minimum eigenvalue of M is the normal vector of P; The point cloud normal vector represents the direction of the plane formed by its neighboring points. After building a kd-tree for all point clouds, a 1m search radius and 10 nearest neighbor points are used to calculate the point cloud normal vector. Step 3: Based on the establishment of the normal vector, initialize the Points point set, assign each point a normal vector attribute value and set the category of all points to unclassified; Step 4: When When Z n Can be and The cross product of Z n Indicates the verticality of the point cloud normal vector. When Z n When it approaches 0, it means that the normal vector is almost perpendicular to the Z axis; therefore, select Z n The area with a value greater than 0.05 and a point cloud height greater than 0.05m and less than 0.2m is used as the initial seed point for subsequent growth, and the points that meet the conditions are added to SeedPoints; Step 5: The algorithm uses a segmented growth method and adds the seed points of the previous segment as the ignition seed points of the new segment; the single tree point cloud is layered with a 0.3m interval and a 0.3m layer thickness, and the layered point cloud is added to the Temp Points collection. The pseudo code for point cloud slicing and layering is as follows: Step 6: For each layer of temporary points, loop and add each seed point, and calculate the Euclidean distance of the seed point to the 10 nearest neighboring points, and determine whether the angle between each point and the normal vector of the seed point is less than the angle threshold A. T ; Change the category of the points that meet the conditions to trunk points, and at the same time change the points that meet the conditions to seed points Z n The difference is less than Z T The point is recorded as the new seed point; Step 7: Loop through each slice point cloud and gradually add the points that meet the conditions to the trunk point cloud set to generate the final trunk point cloud; S104: Double Constrained Convex Hull Algorithm The double-constrained convex hull algorithm is an improvement on the Graham convex hull algorithm. Using the same method as the slicing algorithm in S103, the tree trunk point cloud is sliced ​​and fed into the optimization algorithm one by one to accurately extract boundary points and calculate the area of ​​a single layer. This algorithm primarily uses angles and CHL convex hull side length thresholds to gradually approximate the convex hull boundary established by Graham to the actual distribution of the point cloud, making the convex hull more closely fit the tree trunk surface. The discretized Green's formula is then used to calculate the convex hull area of ​​each layer. Sort the generated convex hull in counterclockwise direction. At this time, for the S vertices of the convex set (x1,y1),(x2,y2),...,(x n ,y n ), the area calculation formula is (Formula (3)): Where, Area is the convex hull area calculated by the optimized convex hull algorithm, (x i ,y i ) The coordinates of the convex hull vertices arranged in counterclockwise order, where n is the number of convex hull points; the initial point is usually repeatedly added to the end of the point set to close the polygon; S105: Hough circle fitting algorithm Due to the scanning method of ground-based lidar data, areas above an average height of 15 meters are difficult to capture, which easily leads to missing point cloud data and incomplete tree trunk point clouds. The double-constrained convex hull algorithm has difficulty completing the point cloud when data is missing, and the calculated single-layer area is too small. Hough circle fitting is used to fit the shape of the missing points at higher locations and calculate the area to compensate for the inability to accurately calculate the upper layer area. The specific principles of the Hough circle fitting algorithm are as follows: Traverse each edge point (x i ,y i ): Calculate the possible center of the circle for all possible radii (a k ,b k ), and corresponding to the accumulator A[a k ][b k ][r k ] plus 1; search for local maxima exceeding the threshold in the accumulator, these maxima are the parameters (a*, b*, r*) of the detected circle; S106: Calculation of volume of single timber The volume of a single tree is calculated using the idea of ​​the infinitesimal method. The trunk of each slice layer is assumed to be a number of frustum-like cones. The frustum calculation formula is used to calculate the volume of the trunk of each layer (Formula (4)): Because the improved convex hull algorithm DB is more sensitive to missing point cloud data, circle fitting is used to calculate the area when the point cloud data is incomplete. To further realize the automatic calculation of material volume, the algorithm designs an area ratio threshold ARatio to automatically determine whether the area ratio of the DB algorithm and the circle fitting is greater than ARatio. If it is greater than the threshold, it indicates that the DB algorithm is greatly affected by the missing data, and the area of ​​this layer is calculated using the result of circle fitting. To further correct the area results of each layer and prevent the occurrence of larger or smaller values ​​of fitting, the algorithm uses a standard deviation threshold to assist in constraining the area of ​​each slice layer. The reliability of each layer area is verified by comparing the area of ​​the current layer with the area of ​​the previous layers. A decrease in area usually indicates that the normal stem is gradually becoming thinner, while a larger mutation indicates that there is an incorrect fitting. In the case of an area mutation, the algorithm will discard the data of the current layer to prevent incorrect calculation (Formula (5)): Diff>K*S (5) where Diff={a i -a j |a i ,a j ∈P,i≠j}, represents the difference between the area of ​​the current layer and the average area of ​​the previous 10 layers; K represents the control multiple of the threshold; S represents the standard deviation of Diff; Although the ground-based lidar system can collect high-density point clouds, its ability to collect information about objects above 15 meters is weak, and it cannot capture relatively complete point cloud information, resulting in the inability to accurately identify tree trunks at high altitudes. The algorithm uses cone fitting to calculate the volume of tree trunks at high altitudes (Formula (6)): Where S n represents the area of ​​the highest point cloud detected; H represents the height of a single tree; slice represents the number of slice layers; Δh represents the slice layer height; The volume of a single log is the sum of the volume of the logs in the following layers and the volume of the fitted cone in the top layer, that is, V t =V i +V j (7)

Images