A Deep Learning-Based Method for Separating Branches and Leaves from Forest Tree Laser Point Clouds

Abstract

This invention discloses a method for separating branches and leaves from forest tree laser point clouds based on deep learning. The method includes acquiring forest tree laser point cloud data; dividing the de-noising laser point cloud data into aboveground points and ground points; performing branch and leaf separation on the aboveground points using machine learning algorithms and manual correction; performing voxelization on the aboveground points after branch and leaf separation; using the point cloud data with branch and leaf separation and voxelization as a training sample dataset; augmenting the training sample dataset using data augmentation methods; training a deep learning network using the augmented training sample dataset to obtain a trained deep learning network model; and achieving branch and leaf separation within voxels using the trained deep learning network model. The deep learning network model constructed in this invention can more effectively extract global and local feature information from forest tree point clouds, achieving branch and leaf separation with high accuracy and good stability.

Description

Technical Field

[0001] This invention belongs to the field of forestry technology, specifically relating to a method for separating branches and leaves from laser point clouds of trees based on deep learning. Background Technology

[0002] Accurate reconstruction of tree skeleton spatial models plays a crucial role in the research of forest information resource surveys and the inversion of tree phenotypic structural features. In recent years, methods for tree model reconstruction can be broadly divided into two categories: image-based and laser point cloud-based tree reconstruction methods. These two methods differ in data acquisition but share certain similarities in their methods of extracting tree spatial features.

[0003] Image-based tree trunk reconstruction methods typically employ binocular or multi-view vision and depth camera techniques to capture tree texture and depth information, combining curvature constraints to reconstruct the tree trunk. Alternatively, multiple images from a multi-view vision system can be used to reconstruct tree trunks through interactive editing, and leaf models can be reconstructed through leaf segmentation. While image-based trunk reconstruction methods are simple and easy to implement, the extracted tree trunks are in a two-dimensional plane, making it difficult to intuitively reflect the three-dimensional morphology of the tree. Furthermore, the extraction results are susceptible to interference from various factors such as diverse environmental backgrounds, complex tree trunk spatial topology, stereo matching errors, and viewpoint occlusion.

[0004] With the advancement of laser measurement technology, the digital modeling of trees has gradually shifted from methods based on two-dimensional images to generation from three-dimensional laser point clouds. Laser scanning boasts high precision and density, and can be combined with theoretical algorithms from computer graphics and machine vision to process and analyze three-dimensional point clouds of trees. For example, machine learning methods can be used to classify branches and leaves in forest point clouds; graph theory can be used to store the collected point cloud data as an octree structure; the PROSAC algorithm, which combines point normals, can be used to model and extract the tree trunks; or space-filling triangular patches can be used to construct each leaf to build a tree canopy model, etc.

[0005] With the continuous breakthroughs in the application of next-generation artificial intelligence technologies and the practical application of laser scanning technology, deep learning technology has begun to emerge in 3D point cloud processing, giving rise to several well-known deep network models. For example, the PoineNet network proposed in 2016 designed a multilayer perceptron based on symmetric functions and the T-Net architecture based on spatial transformation invariance, solving the problems of disorder and partial rotation invariance in input point clouds. Three deep learning networks for laser point cloud processing emerged in 2017: VoxelNet, which divides 3D point clouds at the voxel scale and performs classification, detection, and position regression; PointNet++, which borrows the idea of ​​multilayer receptive fields from convolutional neural networks (CNNs) and combines spherical neighborhoods to extract local spatial features of point clouds; and PointCNN, which designs χ-transform and continuous weighted function convolution operations to achieve spatial target point cloud convolution, local feature extraction, and classification. In addition, deep network modules DPAM based on graph theory for point cloud sampling, grouping, and pooling, and 3DContextNet, which combines a hierarchical data structure (Kd-tree), have also been proposed.

[0006] Although deep learning has seen frequent innovations in industrial applications in recent years, its application in forestry remains limited. Some works involve using convolutional neural networks to identify or directly analyze spatial forest point cloud data projected onto a two-dimensional plane. Examples include constructing convolutional neural networks and deep belief networks to process digital surface models (DSMs) of the University of Kentucky Robinson Forest and lateral projection maps of trees in Xiamen City generated from laser point clouds for classification of coniferous and broadleaf tree species; or employing deep network models such as Superpoint graph, PointNet, U-Net, Kd-Net, and Faster R-CNN to perform tasks at the voxel scale, such as branch and leaf separation of European beech, single-tree crown segmentation of subtropical forests in Chizhou, Anhui, identification of some tree species in Finland, Canada, and parameter extraction of rubber forests in Hainan.

[0007] Despite some progress made in laser point cloud analysis using deep learning, the following problems still exist: 1) The branches and trunks in the tree canopy have diverse topological structures and forms, and are subject to mutual occlusion, which places high demands on the robustness and universality of the network. How to design a more accurate deep network for branch and leaf point cloud segmentation is still under research; 2) Deep learning can map the two-dimensional data to be analyzed into a high-dimensional space through convolution operations to augment and compress the feature dimensions. However, given the disorder and rotational changes of three-dimensional point clouds, it is difficult to set up convolution operations in a multi-scale space and extract the effective global and local features of the point cloud. Summary of the Invention

[0008] The technical problem to be solved by the present invention is to provide a method for separating branches and leaves of forest laser point clouds based on deep learning, which addresses the shortcomings of the prior art. This method for separating branches and leaves of forest laser point clouds constructs a deep learning network for classifying branch and leaf point clouds with multiple feature encoding layers and feature decoding layers. It can more effectively extract global and local feature information of forest point clouds. The point cloud branch and leaf separation is achieved through this deep learning network model with high accuracy and good stability.

[0009] To achieve the above-mentioned technical objectives, the technical solution adopted by the present invention is as follows:

[0010] A deep learning-based method for separating branches and leaves from laser point clouds of trees includes:

[0011] Step 1: Acquire forest tree laser point cloud data;

[0012] Step 2: De-noise the collected laser point cloud data and separate the de-noiseed laser point cloud data into ground points and ground points;

[0013] Step 3: Use machine learning algorithms and manual correction to perform branch and leaf separation on the ground points, and then perform voxelization on the ground points that have undergone branch and leaf separation.

[0014] Step 4: Use the point cloud data that has already undergone branch and leaf separation and voxelization as the training sample dataset, and use the data augmentation method to augment the training sample dataset to obtain a new training sample dataset.

[0015] Step 5: Construct a deep learning network and train it using a new training sample dataset to obtain a trained deep learning network model.

[0016] Step 6: Collect laser point cloud data of the forest to be tested, perform noise reduction processing on the collected laser point cloud data, and divide the noise-reduced laser point cloud data into ground points and ground points; perform voxel subdivision on the ground points, and input the point cloud data in a single voxel into the trained deep learning network model to achieve branch and leaf separation of the point cloud data in the voxel.

[0017] As a further improved technical solution of the present invention, the deep learning network includes a feature encoding layer and a feature decoding layer; the feature encoding layer is used to extract point cloud feature information using a downsampling method, and the feature decoding layer is used to transmit point cloud feature information using an upsampling method;

[0018] The feature encoding layer includes a sampling module, a first grouping module, and a first PointConv module; the feature decoding layer includes a feature interpolation module, a second grouping module, and a second PointConv module.

[0019] As a further improvement of the present invention, there are four feature encoding layers and four feature decoding layers.

[0020] As a further improvement to the present invention, the calculation process of the sampling module in the j-th feature coding layer is as follows:

[0021] Input size is (d+C) j )×N j Given a point cloud matrix, N points are selected from the point cloud matrix using the farthest point sampling method. j+1 Each sub-sampling point yields a size of (d+C) j+1 )×N j+1 The point cloud matrix;

[0022] Where: (d+C j )×N j It refers to having d-dimensional coordinates and C j N-dimensional features j Point, (d+C) j+1 )×N j+1 It refers to having d-dimensional coordinates and C j+1 N of high-dimensional features j+1 Each sampling point.

[0023] As a further improvement of the present invention, the calculation process of the first grouping module in the j-th feature coding layer is as follows:

[0024] The d×N output of the input sampling module j+1 The coordinates of the sub-sampling points and the initial input (d+C) j )×N j The feature point set is classified using the nearest neighbor rule method, and N is extracted from the sampling module. j+1 For each of the sampling points, find the K nearest neighbors in the neighborhood of each point to form a group, resulting in d×K×N. j+1 The set of coordinate indices of nearest neighbors relative to each sampling center point, and the result of (d+C) based on the coordinate indices. j )×K×N j+1 The feature set of nearest neighbors, where K is N j+1 The number of nearest neighbors in the neighborhood of each sampling center point.

[0025] As a further improvement of the present invention, the input of the first PointConv module in the j-th feature encoding layer is three parts of data; wherein the first part of data is: N j+1 The features of K nearest neighbors around a center point, with a size of (d+C) j )×K×N j+1 The second part of the data is: N j+1The local coordinates of K nearest neighbor points relative to the center point, with a size of d×K×N. j+1 The third part of the data is: N j+1 The density of the K nearest neighbors of a center point is 1×K×N. j+1 The first PointConv module outputs a size of (d+C). j+1 )×N j+1 A new local area.

[0026] As a further improvement of the present invention, the calculation process of the first PointConv module in the j-th feature coding layer is as follows:

[0027] Step (1): Calculate the center point using the kernel density estimation method. The kernel density is estimated for each input point among its K nearest neighbors. Then, the inverse density transform of the K nearest neighbors is calculated using the kernel density estimate. The result of the inverse density transform is input into a multilayer perceptron (MLP2) for nonlinear transformation of the activation function, yielding inverse density coefficients of size 1×K. The inverse density coefficients are then subjected to C0... in The tiling process is repeated 1 time, and the output size is C. in The inverse density tensor of ×K;

[0028] Step (2): Input the local coordinates of the K nearest neighbors relative to the center point into the multilayer perceptron MLP1, and output the size (C). in ×C out A weight tensor W of )×K, where C in For the number of input features, C out To output the number of features;

[0029] Step (3): Set the output of step (1) to C. in The inverse density tensor of ×K and the feature F of the K nearest neighbors in Perform a dot product operation, and the output size is C. in ×K matrix Again Perform C out This tiling operation, repeated 10 times, outputs a size of (C). in ×C out The intermediate tensor of )×K;

[0030] Step (4): Set the output of step (2) to (C) in ×C out The weight tensor W of )×K and the size of the output of step (3) are (C in ×C out Perform a dot product operation on the intermediate tensor of )×K, and the output size is (C) in ×C out A matrix of size K

[0031] Step (5): The matrix obtained in step (4) The center point is obtained by summing the first and last dimensions. The output feature F out ;

[0032] The center point For N j+1 One of the center points.

[0033] As a further improvement of the present invention, the calculation process of the feature interpolation module in the j-th feature decoding layer is as follows:

[0034] N obtained from the sampling module in the (j-1)th feature coding layer j Each point in the sampling points Centered on the j-th feature encoding layer, the sampling module obtains N. j+1 Among the sampling points, select the distance N j Each point in the sampling points The three nearest neighbors Three neighboring points were obtained. In N j+1 The indices of the sampling points, where t = 1, 2, 3;

[0035] Three neighboring points were calculated using the reciprocal of the distance. Relative to N j The corresponding center point among the sampling points weight value The calculation formula is:

[0036]

[0037] in Point With point The Euclidean distance between them, where μ represents the weighted power exponent;

[0038] The initial input (d+C′) of the j-th feature decoding layer is obtained by using the indices of these three neighboring points. j )×N j+1 The feature points obtained from the feature point set are the features (d+C′) corresponding to these three neighboring points. j )×3, then combined with the corresponding weight values The summation and averaging of the features of these three neighboring points yields the following result. The corresponding (d+C′) j Interpolation characteristics of 1×1 Interpolation features The calculation formula is:

[0039]

[0040] Because the (j-1)th feature encoding layer has N j There are 10 sampling points, so we get (d+C′) j )×N j Interpolation characteristics;

[0041] We get (d+C′) j )×N j After determining the interpolation features, the skip connection method is used to connect (d+C′). j )×N j The interpolation features and the output of the (j-1)th feature coding layer (d+C) j )×N j The features are concatenated to obtain the final output of the feature interpolation module: (2d+C) j +C′ j )×N j ;

[0042] Where (d+C′) j )×N j This indicates that there are d-dimensional coordinates and C′. j N-dimensional features j One point.

[0043] As a further improvement of the present invention, the calculation process of the second grouping module in the j-th feature decoding layer is as follows:

[0044] The final output of the input feature interpolation module (2d+C) j +C′ j )×N j The d×N obtained by the sampling module in the (j-1)th feature coding layer j The coordinates of the sub-sampling points are classified using the nearest neighbor rule method in N. j For each of the sampling points, find the K' nearest neighbors in the neighborhood of each point to form a group, resulting in d×K′×N. j The set of coordinate indices of nearest neighbors relative to each sampling center point, and the result of (2d+C) based on the coordinate indices. j +C′ j )×K′×N j The feature set of nearest neighbors, where K′ is N j The number of nearest neighbors in the neighborhood of each sampling center point.

[0045] As a further improvement of the present invention, the input of the second PointConv module in the j-th feature decoding layer is three parts of data; wherein the first part of data is: N j The features of K′ nearest neighbors around a center point are given, with a size of (2d+C).j +C j ′)×K′×N j The second part of the data is: N j The local coordinates of K′ nearest neighbors relative to the center point, with a size of d×K′×N, are given by the center point. j The third part of the data is: N j The density of K′ nearest neighbors around a center point is 1×K′×N. j The output size of the second PointConv module is (d+C′). j-1 )×N j New local area.

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

[0047] This invention uses voxelized segmentation of forest laser point cloud data with branch and leaf labels as the training set; it constructs a deep learning network for branch and leaf point cloud classification with multiple feature encoding and feature decoding layers. This deep learning network model includes an improved PointConv module that uses kernel density estimation to calculate local point cloud features and a point cloud feature interpolation module that uses multi-scale feature fusion. It can more effectively extract global and local feature information of forest point clouds. The point cloud branch and leaf separation is achieved through this deep learning network model with high accuracy and good stability, without the need for manual parameter adjustment and result intervention correction. Attached Figure Description

[0048] Figure 1 (a) is a schematic diagram showing the location of Danzhou City, Hainan Province, the study area.

[0049] Figure 1 (b) is a schematic diagram of three rubber plantations in a rubber tree plantation in Danzhou City, Hainan Island.

[0050] Figure 2 (a) is a schematic diagram of the branch and leaf segmentation and voxel subdivision results obtained by using machine learning methods and manual correction for some sample plots of rubber tree variety PR107.

[0051] Figure 2 (b) is a schematic diagram of the branch and leaf segmentation and voxel subdivision results obtained by using machine learning methods and manual correction for a portion of the sample plots of the rubber tree variety CATAS 7-20-59.

[0052] Figure 2 (c) is a schematic diagram of the branch and leaf segmentation and voxel subdivision results obtained by using machine learning methods and manual correction for some sample plots of rubber tree variety CATAS 8-79.

[0053] Figure 3 This is a diagram of a deep learning network architecture for branching and separating point clouds.

[0054] Figure 4 (a) is the center point p i A schematic diagram illustrating the calculation of local features in conjunction with its nearest neighbors.

[0055] Figure 4 (b) is the improved PointConv module framework diagram.

[0056] Figure 5 (a) is a schematic diagram of the classification accuracy curve of a deep learning network.

[0057] Figure 5 (b) is a schematic diagram of the loss function curve of a deep learning network.

[0058] Figure 6 (a1) is a schematic diagram of laser point cloud data and branch and leaf separation results of the test plot of rubber tree variety PR107.

[0059] Figure 6 (a2) shows the laser point cloud data and branch-leaf separation results of the test plot of rubber tree variety CATAS 7-20-59.

[0060] Figure 6 The middle (a3) ​​shows the laser point cloud data and branch and leaf separation results of the test plot of rubber tree variety CATAS 8-79.

[0061] Figure 6 The middle (b1) diagram shows the fitting results of the cylindrical branches of the test plot of rubber tree variety PR107.

[0062] Figure 6 Figure (b2) shows the fitting results of the cylindrical branches of the test plot of rubber tree variety CATAS 7-20-59.

[0063] Figure 6 Figure (b3) shows the fitting results of the cylindrical branches of the test plot of rubber tree variety CATAS 8-79.

[0064] Figure 6 The middle (c1) diagram shows the classification results of the main branches and primary branches of the rubber tree variety PR107 in the test plot.

[0065] Figure 6 The middle (c2) diagram shows the classification results of the main branches and primary branches of the rubber tree variety CATAS 7-20-59 in the test plot.

[0066] Figure 6 The middle (c3) diagram shows the classification results of the main branches and primary branches of the rubber tree variety CATAS 8-79 in the test plot.

[0067] Figure 6The middle (d1) plot shows the leaf clustering results and leaf cluster convex hull display of the test plot of rubber tree variety PR107 based on spatial watershed and first-order branching.

[0068] Figure 6 The middle (d2) plot shows the leaf clustering results and leaf cluster convex hull display of the test plot of rubber tree variety CATAS 7-20-59 based on spatial watershed and first-order branching.

[0069] Figure 6 The middle (d3) plot shows the leaf clustering results and convex hull display of leaf clusters based on spatial watershed and first-order branching in the test plot of rubber tree variety CATAS 8-79.

[0070] Figure 7 The figure in the middle (a) shows the comparison between the measured and calculated values ​​of the diameter at breast height (DBH) of the trees in the three rubber tree plots.

[0071] Figure 7 Figure (b) shows the comparison between the calculated angles between the main branches and first-order branches of the trees in the three rubber tree plots and the measured values.

[0072] Figure 7 The middle (c) plot shows the correlation between the width of the first-order branches of the trees and the corresponding leaf cluster volume in the three rubber tree plots. Detailed Implementation

[0073] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings:

[0074] This embodiment proposes a tree skeleton reconstruction and forest parameter inversion method that integrates deep learning and computer graphics. First, voxelized segmentation of rubber forest laser point cloud data with branch and leaf labels is used as the training set. Next, a deep learning network for branch and leaf point cloud classification with multiple feature encoding and decoding layers is constructed. This network model includes an improved PointConv module that uses kernel density estimation to calculate local point cloud features and a point cloud interpolation module that uses multi-scale feature fusion, which can more effectively extract global and local feature information from the forest point cloud. Third, for the segmented tree branch and trunk point cloud, computer graphics methods are used to establish the connectivity of the skeleton's center points, and cylinders are used to perform segmented branch and trunk fitting, thereby reconstructing a three-dimensional model of the entire forest skeleton from the branch and trunk point cloud. Finally, compared with manually measured data, each tree in the sample plot is finely characterized at the leaf cluster scale, and the topological structure and phenotypic parameters of different rubber tree varieties are analyzed.

[0075] 1.1 Study Area and Data Acquisition:

[0076] The study area is located in Danzhou City, Hainan Province, specifically in the experimental field of the Rubber Research Institute of the Chinese Academy of Tropical Agricultural Sciences (19°32′47.89″N, 109°28′29.33″E). Figure 1 As shown, Hainan Island, China's second-largest rubber production base, has a tropical monsoon climate, with a rainy season from May to October and a dry season from November to April. The average annual rainfall is 1815 mm, and the average annual temperature is approximately 23.2℃. This climate is conducive to the development of the rubber industry, and the planting area of ​​rubber trees in the region is constantly increasing. After decades of introduction and experimentation, rubber tree varieties such as RRIM660, PR107, PB86, Haiken No. 2, CATAS 7-20-59, CATAS 8-79, and Reken 126 have gradually developed into excellent rubber tree varieties suitable for large-scale planting in Hainan, characterized by resistance to various natural disturbances, strong regeneration ability, and high yield. Among them, PR107 has a moderate growth rate, high dry rubber content, resistance to irritation, and good wind resistance, making it an excellent high-yielding rubber tree variety with wind resistance; CATAS 7-20-59, a hybrid offspring of PR107 and RRIM600, is fast-growing, high-yielding, and relatively late-maturing, making it a new rubber tree variety with good wind resistance and currently widely promoted; CATAS 8-79 has strong regeneration ability, early maturity, high yield, and stable yield, making it a high-yielding line with excellent comprehensive traits. Therefore, as... Figure 1 As shown, rubber trees of variety PR107 (plot 1), variety CATAS 7-20-59 (plot 2), and variety CATAS 8-79 (plot 3) were selected as typical tree species for the experiment.

[0077] LiDAR data was acquired on October 15, 2019, using a Velodyne HDL-32E high-definition LiDAR sensor that can operate in backpack mode. This sensor has 32 pairs of lasers / detectors and was used to measure rubber tree maps with the following settings: +10.67° to -30.67° vertical field of view (FOV) with an angular resolution of 1.33°, 360° horizontal field of view with an angular resolution of 0.16°, a frame rate of 10 Hz, and a measurement range of 70 m. The Velodyne HDL-32E scanning system was carried by an experimenter, and the scanner was set to "continuous shooting mode," acquiring data at a speed of 10 revolutions per second. The experimenter carried the Velodyne laser scanning system and walked within three rubber tree plots according to a predefined measurement route. The survey route was programmed as a predefined "back-and-forth parallel" route (e.g., ...). Figure 1 (Dashed line) Two parallel mapping routes were spaced 6m apart to cover as many trees as possible across the three study plots. Due to the complex terrain of the rubber plantations and the heavy scanning equipment, the researchers moved along the survey routes at a speed of 0.5m / s. The Velodyne lidar system, combined with Simultaneous Localization and Mapping (SLAM) technology, rapidly completed registration for each scan and generated high-density point clouds for each target rubber tree. The average resolution of the lidar data obtained from the three rubber tree plots was approximately 2cm.

[0078] Figure 1 To provide an overview of the study area, Figure 1 Image (a) shows the location of the study area, Danzhou City, Hainan Province. Figure 1 Image (b) shows three rubber plantations within a rubber plantation in Danzhou City, Hainan Island. Figure 1 (b) shows a remote sensing image of three rubber plantations on Google Maps. Different rectangles mark the edges of the different rubber plantations, and dashed lines indicate the detection routes of the airborne lidar.

[0079] 1.2 Training Samples and Test Samples:

[0080] After acquiring point cloud data of trees using a laser scanner, Gaussian filtering was used to denoise the collected point cloud, and ground point filtering (CSF) was applied to separate the denoised point cloud into aboveground and ground points. Then, traditional machine learning algorithms combined with manual visual inspection and result correction were used to perform branch and leaf separation on the aboveground points in the experimental area. Three subsets were created from three rubber plantation plots as training samples for the subsequent deep network. Each subset consisted of an area of ​​approximately 200 × 100 m representing the corresponding rubber plantation plot. Figure 2 The left side of (a), (b), and (c) shows the branch and leaf separation results of partial subsets of data from three rubber plantation plots. Since the different rubber tree varieties (PR107, CATAS7-20-59, and CATAS 8-79) on the three plots have the same east-west (3m) and north-south (7m) planting spacing, and there are few understory shrubs, the same voxel size (3m long, 7m wide, 7m high) was set for the three training plots to voxelize the ground tree point clouds, ensuring that each tree is located at the center of each voxel. After partitioning, all laser point clouds from the three training sets were assigned to corresponding voxels, ensuring that the point cloud of a single tree could be partitioned by two consecutive voxels above and below. Next, a thinning process was used, and the point cloud within each voxel was randomly sampled into 4096 points. Figure 2 The right side of (a), (b), and (c) shows the extracted point clouds of branches and leaves of a single tree, bounded by upper and lower voxels. After the above operations, a total of 2743 rubber trees (approximately 900 trees per variety) with separated branches and leaves were collected as the training set for deep learning.

[0081] The test sample consisted of three subsets of plots representing three rubber tree varieties, ensuring no overlap with the training sample. Each test plot was a 35×35m area. Plot data were acquired using a backpack LiDAR system at the same scanning resolution. Specific parameters of the rubber trees throughout the study area are shown in Table 1, where crown width and crown volume were calculated using a single-tree separation algorithm.

[0082] Table 1 shows the parameters of rubber trees in the study plot:

[0083]

[0084] A large and sufficient training sample is fundamental for optimizing deep learning network parameters. To avoid overfitting and improve the network's generalization ability, it is necessary to collect as many training sample sets as possible to participate in the training process of the deep learning network to optimize neuron weights. In this study, data augmentation techniques are used to obtain more training samples. Without collecting and processing new tree point cloud data, based on the existing tree point cloud within voxels, a method is proposed to move each point in the voxel along a random vector with a small offset. The strategy, namely Gaussian noise jitter, is shown in equations (1) and (2).

[0085]

[0086] Where a and b are two random numbers between 0 and 1, λ is a power amplification factor, and r is the generated perturbation. Then, the perturbation at each point... Substituting into equation (2), the original tree point cloud is changed from p i Change to p′ i In equation (2) μ noise and σ noise Let be the mean and variance of the Gaussian noise.

[0087]

[0088] Alternatively, by rotating the point cloud within the voxel along the Z-axis (i.e., the central axis of the tree crown) by a small angle (e.g., 1–2°) according to the spatial rotation equation, and by adjusting the equation parameters, a smaller standard deviation (range 1–5 cm) compared to the original tree point cloud can be ensured after the transformation. In summary, approximately 1100 samples were added to each rubber tree variety, increasing the total number of training samples to 6000.

[0089] Figure 2 The image shows the branch and leaf segmentation and voxel subdivision results of the point cloud data from three rubber plantation training plots obtained using a machine learning and manual correction method. The right side displays magnified views of the point clouds of branches and leaves of a single tree, bounded by voxels. The categorized leaf point clouds are shown in gray, and the categorized branch point clouds are shown in black. These data will be used as the training sample set to input into the deep learning network for parameter acquisition. Figure 2 (a) is a schematic diagram of the branch and leaf segmentation and voxel subdivision results obtained by using machine learning methods and manual correction for some sample plots of rubber tree variety PR107. Figure 2 (b) is a schematic diagram of the branch and leaf segmentation and voxel subdivision results obtained by using machine learning methods and manual correction for a portion of the sample plots of the rubber tree variety CATAS 7-20-59. Figure 2(c) is a schematic diagram of the branch and leaf segmentation and voxel subdivision results obtained by using machine learning methods and manual correction for some sample plots of rubber tree variety CATAS 8-79.

[0090] 1.3 Tree branch and leaf separation based on deep learning network:

[0091] The deep learning network constructed in this study was trained using point clouds with branch and leaf labels within a single voxel as input, and the weights of the neurons in the deep learning network were optimized. The specific network architecture is as follows: Figure 3 As shown, the network here mainly consists of four feature encoding layers and four feature decoding layers.

[0092] 1.3.1 Feature Encoding Layer in Deep Learning Networks:

[0093] Deep learning networks have four feature encoding layers. Figure 3 The layers in the network, from left to right, are designated as Layer 1, Layer 2, Layer 3, and Layer 4, respectively, and are denoted as Feature Coding Layer 1, Feature Coding Layer 2, Feature Coding Layer 3, and Feature Coding Layer 4. The feature coding layers use a downsampling method to extract point cloud features. Each feature coding layer includes a sampling module, a first grouping module (referred to as the grouping module), and a first PointConv module (i.e., the optimized PointConv module, which can be simply referred to as the PointConv module). The j-th feature coding layer starts with a (d+C)... j )×N j The point cloud matrix is ​​used as input, that is, the input has d-dimensional coordinates and C... j N-dimensional features j There are 1 point. The output is a (d+C) j+1 )×N j+1 The matrix, i.e., the output has d-dimensional coordinates and C j+1 Summary of high-dimensional features and local context N j+1 A two-dimensional matrix composed of sub-sampling points.

[0094] The sampling module of the feature encoding layer: inputs the point set {p1, p2, ..., p} within the current voxel. n A subset is selected from it using the farthest sampling method. This process extracts some of the more important center points, thereby reducing the number of point cloud samples and preserving the original spatial structure of the point cloud to a great extent. The input to the sampling module is a set of data of size 3×N. j The coordinates of the point are given, and the output is of size 3×N. j+1 The coordinates of the sub-sampling points.

[0095] The first grouping module of the feature encoding layer: uses the nearest neighbor rule classification method, extracting N from the sampling module. j+1 For each of the given points, find the K nearest neighbors within its neighborhood to form a group. Combine these neighbors with the center N. j+1 Feature calculation is performed on each point. The input to the grouping module is the d×N output of the sampling module. j+1 The coordinates of the sub-sampling points and the initial input size are (d+C) j )×N j Given a feature point set, output two results: d×K×N. j+1 The set of coordinate indices of nearest neighbors relative to each sampling center point and (d+C) j )×K×N j+1 The feature set of nearest neighbors, where K is the number of nearest neighbors in the neighborhood of the sampling center point, is set to K=32 in the coding layer.

[0096] The first PointConv module of the feature encoding layer, namely the optimized PointConv module, solves the convolution difficulty caused by the disordered spatial arrangement of point clouds: it uses the density function to optimize the N extracted by the first grouping module. j+1 The features of K nearest neighbors around the center point are weighted and convolved to obtain the j-th feature encoding layer N. j+1 The local correlation features of each center point are used to input the feature values ​​of this layer into the (j+1)th feature encoding layer. The input of the first PointConv module consists of three parts of data, namely N j+1 The grouped feature of the K nearest neighbors of each center point is of size (d+C). j )×K×N j+1 N j+1 The local coordinates of K nearest neighbors relative to the center point, grouped_xyz, are of size d×K×N. j+1 N j+1 The density of the K nearest neighbors of a center point is 1×K×N. j+1 The output data size after passing through the first PointConv module is (d+C) j+1 )×N j+1 New local area.

[0097] For example, Figure 3The first feature encoding layer (first dashed box from the left in the first row of the figure) is shown. The input is a point cloud matrix of (3)×4096, where (3) is the sum of the coordinate dimension d=3 and the feature dimension C=0 of the point cloud, and 4096 is the total number of points in the point cloud. This input first enters the sampling module, which uses the farthest point sampling method to find N2=2048 points and obtains the coordinates of the 3×2048 sampling center points (sampled_xyz). The coordinates of the sampling center points and the feature point set of (3)×4096 of the original input then enter the first grouping module. The first grouping module finds the nearest K=32 nearest neighbors in the neighborhood of each of the N2 center points extracted by the sampling module to form a group, and obtains the coordinate index set of the nearest neighbors relative to each sampling center point, 3×32×2048 (grouped_xyz). Based on the coordinate index, the feature set of the nearest neighbors (3)×32×2048 (grouped_feature) is obtained. The coordinate dimension in the first feature decoding layer is the number of features. Meanwhile, the coordinate index set (grouped_xyz) of the nearest neighbor points obtained by the first grouping module (3)×32×2048 is used to obtain the density of the nearest neighbor points (1×32×2048) through the kernel density estimation algorithm. This density of the nearest neighbor points, together with the coordinate index set (grouped_xyz) of the nearest neighbor points and the feature set (grouped_feature) of the nearest neighbor points, is used as the input of the first PointConv module. According to the multilayer perceptron (MLP) and other operations in the first PointConv module, the final output of the feature encoding layer is obtained, that is, the output feature set of 2048 points of the first feature decoding layer, which is (64)×2048.

[0098] The second, third, and fourth feature encoding layers in the deep learning network constructed in this study operate on similar principles. Through the encoding layers, center points within voxels are continuously extracted, reducing their number from 4096 to 72, while simultaneously increasing the information contained in the center points (the feature dimension of the center points increases from 3 to 512). After four feature encoding layers, high-dimensional feature extraction of the point cloud is achieved.

[0099] 1.3.2 Feature Decoding Layer in Deep Learning Networks:

[0100] In the deep learning network constructed in this paper, the number of feature decoding layers is the same as the number of feature encoding layers, which is also 4 layers. Figure 3The network layers, arranged from right to left in the order of network operation, are designated as Layer 4, Layer 3, Layer 2, and Layer 1, respectively, and are denoted as Feature Decoding Layer 4, Feature Decoding Layer 3, Feature Decoding Layer 2, and Feature Decoding Layer 1. The feature decoding layers utilize upsampling to transmit point cloud feature information. Each feature decoding layer includes a feature interpolation module, a second grouping module (referred to as the grouping module), and a second PointConv module (i.e., the optimized PointConv module, which can be simply referred to as the PointConv module). At the j-th feature decoding layer, the point cloud features augmented after feature decoding are processed from (d+C′). j )×N j+1 Reduced to (d+C′) j-1 )×N j The number of point clouds gradually increases after decoding, with N... j+1 and N j It is the size of the point set of the input and output of the j-th layer of feature decoding.

[0101] Feature interpolation module: N obtained by combining the sampling modules in the (j-1)th feature coding layer j sampling points The N obtained by the sampling module in the j-th feature encoding layer j+1 sampling points Feature association is performed, and then inverse distance weighted interpolation (see Equation 4) is used to achieve feature propagation. For example, N obtained after sampling from the (j-1)th feature encoding layer... j Taking each of the sampling points as the center, N is obtained after sampling in the j-th feature coding layer. j+1 Among the sampling points, a distance of N is selected. j Each point inside The three nearest neighbors Three neighboring points were obtained. In N j+1 The index of each sampling point is used to calculate the three neighboring points based on the reciprocal of the distance. Relative to N j The corresponding point among the points weight value The calculation is shown in equation (3), where the L2 norm ||||2 represents the Euclidean distance between the two points, and μ is the weighted power exponent, with a default value of 2.

[0102]

[0103] Next, using the indices of these three points, the initial input (d+C′) of the j-th feature decoding layer is... j )×N j+1 The feature points obtained from the feature point set are the features (d+C′) corresponding to these three neighboring points. j )×3, then combined with the corresponding weight values The characteristics of these three points are summed and averaged to obtain... The corresponding (d+C′) j Interpolation characteristics of 1×1 Because the (j-1)th feature encoding layer has N j There are sampling points, so a total of (d+C′) are obtained. j )×N j interpolation features ( Figure 3 (The small box with the subscript "Interpolation Result" is shown in the middle), where each interpolation feature... The specific calculation is as follows:

[0104]

[0105] We get (d+C′) j )×N j After interpolating the features, skip links are used to combine the interpolated features with the output features (d+C) of the (j-1)th feature encoding layer. j )×N j Perform a connection operation on (lj-1_points) to obtain the final output of the feature interpolation module (2d+C). j +C j ′)×N j ( Figure 3 (The small box with a gray background is shown below).

[0106] The second grouping module of the feature decoding layer is similar to the first grouping module in the feature encoding layer. The input is the final result of the feature interpolation module (2d+C). j +C j ′)×N j The d×N obtained by the sampling module in the (j-1)th feature coding layer j The coordinates of the sampling center point (lj-1_xyz). Using these N... j For each of the given points, select K′ nearest neighbors within its own neighborhood, for a total of (2d+C). j +C′ j ) features, the expanded results are (2d+C j +C′ j )×K′×N j and d×K′×N j This is the output of the second grouping module. Here, K′ is the number of nearest neighbors within the center point's neighborhood, which is set to K′ = 16 in the decoding layer.

[0107] The second PointConv module in the feature decoding layer is similar to the first PointConv module in the feature encoding layer. The input to the second PointConv module consists of three parts of data, namely N...j The grouped feature of K′ nearest neighbors around a center point is of size (2d+C). j +C′ j )×K′×N j N j The local coordinates of K′ nearest neighbors relative to the center point are grouped_xyz, with a size of d×K′×N. j N j The density of K′ nearest neighbors around a center point is 1×K′×N. j The output data size after passing through the second PointConv module is (d+C′). j-1 )×N j New local area.

[0108] For example, such as Figure 3As shown in the fourth feature decoding layer, the input is a feature point set of (512)×72. Taking each point in the 3×128 sampling points (l3_xyz) obtained from the third feature encoding layer as the center, three points that are closest to it and do not include itself are selected from the 3×72 sampling points (l4_xyz) obtained from the fourth feature encoding layer. The indices of these three points in the 72 sampling points (l4_xyz) are obtained, and their respective weight values ​​are calculated according to the inverse distance. Then, the corresponding features (512)×3 of these three points are found in the (512)×72 feature point set (l4_input) initially input to the fourth feature decoding layer through the indices of these three points. The features of these three points are then summed and averaged in combination with the weight values ​​calculated by the inverse distance to obtain the average feature of (512)×1. Since there are 128 sampling points in the third feature encoding layer, the interpolated features of (512)×128 are obtained. The interpolated feature is connected to the (256)×128 feature point set (l3_points) output by the third feature coding layer through skip links, resulting in a (768)×128 feature point set. This feature point set, along with the 3×128 sampling point coordinates (l3_xyz) obtained from the third feature coding layer, is then fed into the second grouping module. The second grouping module searches for the nearest K′ = 16 nearest neighbors in the neighborhood of each of the 128 sampling points to form a group, resulting in a 3×16×128 set of coordinate indices of the nearest neighbors relative to each sampling point (grouped_xyz). Based on the coordinate indices, the feature set (grouped_feature) of the (768)×16×128 nearest neighbors is obtained. Meanwhile, the 3×16×128 nearest neighbor coordinate index set (grouped_xyz) obtained by the second grouping module is used to obtain a point cloud density of 1×16×128 through a density estimation algorithm. This point cloud density, together with the 3×16×128 point coordinate index set (grouped_xyz) and the point feature set (768)×16×128 (grouped_feature), is used as the input of the second PointConv module. According to the multilayer perceptron (MLP) and other operations in the second PointConv module, the final output of the feature decoding layer is obtained, that is, the output result of the 128 points of the 4th feature decoding layer (512)×128. Figure 3 (As shown in feature decoding layer 4 on the right).

[0109] The second, third, and fourth feature decoding layers in deep networks work in a similar way. Thus, these four feature decoding layers propagate features from the subsampled point cloud to a denser point cloud. That is, as the number of layers increases, the number of center points increases (from 72 to 4096), but the feature information contained in each center point decreases (the feature dimension of the center point decreases from 512 to 128). After passing through these four feature decoding layers, the features of all input points are obtained.

[0110] Finally, the fully connected layer compresses the 128-dimensional point cloud features within the voxel to 2 dimensions through a 1×1 convolution operation, which is the confidence level of belonging to the branch and leaf categories. Based on the confidence level, the tree branches and leaves are classified point by point.

[0111] 1.3.3 Improved PointConv modules (i.e., the first PointConv module and the second PointConv module):

[0112] The calculation principle of the first PointConv module and the second PointConv module is the same.

[0113] Traditional convolution operates on two-dimensional image data, which is typically represented as a dense, regularly arranged grid where the relative positions of each pixel are always ordered. In a convolutional neural network, each filter is confined to a small local region, and the convolution operation is performed by sliding the window sequentially. In contrast, three-dimensional point cloud data is a set of unordered three-dimensional points in space, a point cloud set {p}. i |i=1,2,...,n}, where each point contains a set of position vectors (x,y,z) or other information such as color, surface tangents, surface tangent vectors, etc. Unlike images, point clouds have a more flexible representation. Therefore, the order and relative position of a point cloud are different in different local regions. Traditional discrete convolutional filters for images cannot be directly applied to point clouds. Therefore, this paper proposes an improved PointConv network to extract local and global point cloud features.

[0114] Figure 4 The display shows the sampling point p i The operation is performed on a local region consisting of K nearest neighbors of a central point, and the results of the computation with point-by-point input are shown. Figure 4 A framework diagram for the improved PointConv module. Figure 4 (a) represents a local region of the center point and its nearest neighbors, where p i Let p1, p2, p3, and p4 be the center points, and p1, p2, p3, and p4 be the points of interest. i The nearest neighbors of a point are f1, f2, f3, and f4, which are the features of that point. Figure 4 (b) indicates that at point p iThis refers to the process of performing PointConv operations on a local region centered on p. The input features originate from p... i The feature set of the K nearest neighbors of the center is given, and the output is p. i The corresponding feature information.

[0115] The improved PointConv module's input is divided into three parts: 1) grouped_feature is p i The features of the K nearest neighbors, of size C in ×K, C in 1) Input feature number. 2) grouped_xyz represents the K nearest neighbors relative to the center point p. i The local coordinates of point P are 3×K; 3) density is the density of the K nearest neighbors, with a size of 1×K. The output of the PointConv module is the center point P. i Features, size C out ×1, C out This is the number of output features.

[0116] The optimized PointConv module performs the following operations:

[0117] Step 1: Calculate the inverse density: Use the kernel density estimation (KDE) method to estimate the kernel density of each input point in the K nearest neighbor grouped_xyz (i.e., the density of the K nearest neighbors, see Equation 5). Then calculate the inverse density transformation of the K nearest neighbors (see Equation 8). Input the results into MLP2 for nonlinear transformation of the activation function to obtain inverse density coefficients of size 1×K. Finally, perform C... in The next copy tiling The output size is C in The inverse density tensor of ×K.

[0118] Kernel density estimation (KDE) formula:

[0119]

[0120] in, Let h be the density function, and h be the window width. This indicates that for the k-th point p k Select a neighboring point p m Total number The points are obtained from the input set, and the selected points are... For kernel functions, the kernel function used Mainly:

[0121]

[0122] Where B(a,b)=Γ(a)Γ(b) / Γ(a+b), and the gamma function is... When τ = 0, it is a uniform kernel; when τ = 1, it is an Epanechnikov kernel; when τ = 2, it is a Biweight kernel. After repeated experiments, the Biweight kernel was adopted here.

[0123] The optimal window width h is calculated using the thumb rule, and the formula is as follows:

[0124]

[0125] in It is the standard deviation of the distance between the input sample point and its neighboring points. The dimension of the kernel density estimation object is 1 here; the inverse density transformation S is:

[0126]

[0127] Step 2: Calculate the weights: Input the local coordinates of the K nearest neighbors (grouped_xyz) into MLP1, and the output size is W = (C in ×C out The weight tensor is 1×K. The 1×1 convolution acts like a multilayer perceptron, with a ReLU nonlinear layer performing a nonlinear transformation after the 1×1 convolution.

[0128] Step 3: Convert the output of Step 1 into C in The inverse density tensor of ×K and the feature F of the K nearest neighbors in That is, grouped_feature performs a dot product operation (multiplying corresponding elements of two matrices of the same size), and the output size is... The matrix, then... Perform C out The next copy tiling The output size is (C in ×C out The intermediate tensor of )×K.

[0129] Step 4: Convert the output of step 2 (C) in ×C out The weight tensor of )×K and the output of the third step (C) in ×C out Perform a dot product operation on the intermediate tensor of )×K, and the output size is The matrix is ​​given by equation (9), where · represents the dot product.

[0130]

[0131] Step 5: Analyze the output of step 4. Summing the first and last dimensions yields the center point P. i The output feature F out =C out ×1.

[0132] 1.3.4 Loss Function:

[0133] Deep learning networks use the negative log likelihood loss function as the network's loss function, as shown in equations (10) and (11). During training, the loss function is defined as follows:

[0134]

[0135] In formula (10), l = {1, 2} represents the two categories of branches and leaves. i,l This represents an index related to the classification category. If the calculated category corresponds to the i-th point p in the voxel... i If the current category l is the same, then the indicator is assigned a value of 1; otherwise, it is assigned a value of 0. Represent each point p i The probability of the current category is obtained by inversion through a deep learning network. This represents p in the voxel obtained through the Softmax function. i The confidence level of a point belonging to the l-th category.

[0136] 1.4 Tree skeleton reconstruction:

[0137] Based on the branch and trunk point cloud data after branch and leaf classification, the previously proposed method was used to reconstruct the branch and trunk skeleton of rubber trees in three experimental forest plots. The steps are as follows:

[0138] First, the tree trunk height is layered, and the center point of each layer is determined. Specifically, based on the vertical height of a single tree, the tree trunk point cloud is divided into several layers from bottom to top at certain height intervals. The layering of each tree depends on the tree's height and the set height intervals; the number of layers... Δh was uniformly set to 0.3m. After obtaining the branch point cloud for each layer, cluster centers were extracted from the branch point cloud for each layer using the DBSCAN density clustering algorithm. The distance threshold for the DBSCAN algorithm was set to 0.10m, based on the average point cloud resolution of 0.02m for the trees in the three sample plots and multiple experiments. `layer` represents the first layer, used to mark the branch distribution at different levels, thereby obtaining a preliminary outline of the branch connection points of the entire tree.

[0139] Secondly, the extracted center points are classified into root nodes, branch nodes, and edge nodes. Starting from the root node, the tree skeleton is transformed into a chain of connections from the root node to each layer's center point, ending at the edge node, based on the shortest distance to the node above it. The connection chain between each center point belonging to an adjacent layer is represented by a cylinder with spatial directionality and an adaptive radius. The fitting process is performed, where r is the radius of the cylinder, and the central axis of the cylinder is determined based on the corresponding branch segment's scanned point cloud P. b It is obtained through fitting. The specific steps are as follows: in Represents the branches that make up this section of the trunk A point cloud, Let U be the center coordinates of all point clouds in this branch segment, i.e., the average value of the point cloud coordinates. SVD represents singular value decomposition, and the decomposed U is... The matrix, Σ is The matrix V is a 3×3 matrix. Taking the first column of V, V(1) = [V(1,1),V(1,2),V(1,3)], we get the vector of the fitted line. Combined with the center point of the branch point cloud for that segment... The equation of the fitted spatial line L1 is constructed, and its point-normal form is:

[0140]

[0141] Next, calculate all point clouds of this branch segment according to formula (13). The shortest distance d to the fitted line L1 i Among them, q 1, q2 represents any two points on the fitted line L1, and ||2 represents the L2 norm.

[0142] The radius of the cylinder is then set as the average distance from all point clouds of that branch segment to the fitted straight line, i.e.

[0143] Next, the main branches and primary branches of the tree are classified. Specifically, the tree skeleton is divided into main branches and many primary branches. The main branches are identified by finding the connecting chains on the trunk starting from the root node, based on the shortest Dijkstra distance algorithm and the minimum change rule of the tree's growth angle. Primary branches are chains formed by searching from the branch nodes on the main branch chain to the edge node at the other end.

[0144] Finally, based on plant physiological theory, a leaf cluster is defined as a collection of leaves attached to the same main branch or primary branch. Using a spatial three-dimensional watershed segmentation method, the extracted main branches and primary branches are used as cluster centers to segment the leaf clusters. After leaf cluster segmentation, the volume of leaf clusters on a single tree can be calculated using convex hull calculations. Therefore, the three-dimensional reconstruction of a tree can be simplified into a model of different branch skeletons and several leaf cluster combinations.

[0145] 2.1 Training and testing results of deep learning networks:

[0146] Training and testing of the deep learning network were both conducted on a Windows 10 64-bit server equipped with an Intel(R) Core(TM) i7-10750H CPU @ 2.80GHz processor (Intel Inc., Santa Clara, CA, USA) and 16GB of RAM. Since deep learning requires a large amount of training data, an NVIDIA GTX 1650Ti GPU (NVIDIA Inc., Santa Clara, CA, USA) was used instead of a CPU to reduce training time and improve computational power. In the deep learning network constructed in this study, the learning rate was 0.0001, the batch size was 16, and the epochs were 200. Training accuracy and training loss are as follows: Figure 5 As shown, the total training time is approximately 48 hours.

[0147] Figure 5 The training accuracy and loss function curves of the deep learning network in this paper are shown below for different number of epochs. Figure 5 (a) is a schematic diagram of the classification accuracy curve of a deep learning network. Figure 5 (b) is a schematic diagram of the loss function curve of a deep learning network.

[0148] As the learning process progresses, the training accuracy of the training samples (point clouds in each voxel) shows an upward trend, while the training loss shows a downward trend, indicating that the deep learning network in this study is a global optimization process. The training accuracy and training loss show significant increases and decreases respectively in the first 25 epochs, indicating that the weights of neurons in the deep network quickly meet the requirements of the branch and leaf classification task. During training, the neural network may encounter some complex samples in a batch, such as partial voxels containing multiple trees, data gaps caused by occlusion, or blurred leaf and branch shapes. These weaken the model's learning efficiency and lead to drastic fluctuations in the regression loss function value. However, the overall upward and downward trends of the curves indicate that the deep network in this paper has good convergence during training. After 200 epochs, the accuracy and loss of the training samples converge to 97.33% and 0.05 respectively, indicating that the deep learning network constructed in this study has the ability to classify point clouds into branches and leaves.

[0149] Table 2 shows the comparison results between deep learning methods and traditional machine vision algorithms, illustrating the performance of the algorithms in terms of classification accuracy, IOU, and classification time. IOU is defined as follows:

[0150]

[0151] In the above formula and Each represents a point cloud p i The corresponding real tags Category labels identified by the algorithm == means that if they are equal, it returns "True"; if they are not equal (≠), it is recorded as "False". Equation (14) gives the ratio of the number of correctly classified points to the number of incorrectly classified points in all point clouds.

[0152] Test samples were then fed into a pre-trained deep learning network with weights to perform branch and leaf classification. The classification results of different methods are shown in Table 2 and... Figure 6As shown in (a1), (a2), and (a3), traditional machine learning is based on feature description and detection methods. The extraction of target features relies on manual design and acquisition. The entire algorithm lacks forward and backward propagation of features and weights, resulting in limited learning and expressive capabilities, achieving only a 75.24% segmentation accuracy. For forest areas with irregular canopy shapes and complex internal structures, the high density of branches and leaves and overlapping occlusion between canopies can lead to decreased segmentation accuracy. In contrast, deep learning-based 3D point cloud recognition typically employs methods such as multi-view projection, voxel partitioning, and direct input of raw point clouds. It can automatically extract effective features from 3D point clouds using pooling, convolution, activation, feature compression, and augmentation. Compared to traditional machine learning methods, deep learning methods directly use disordered 3D point clouds as input, aggregating more feature information through concatenation of local and global features, and achieving point cloud segmentation based on category confidence. However, PointNet, through multilayer perceptrons and max pooling operations, can only retain globally optimal feature information, losing local point cloud features contained in local groups, thus affecting segmentation accuracy, which is 82.75%. The PointConv method used in this study considers both global feature acquisition and weighted convolution and transformation of local point cloud features through kernel density functions, and combines an interpolation module to achieve multi-scale feature fusion of the point cloud, enhancing the characterization of semantic features of tree point clouds at different scales, thus achieving higher branch and leaf segmentation accuracy (90.32%). Table 2 lists the comparison results of the three methods for branch and leaf segmentation of tree point clouds. The table shows that although the deep learning method takes longer to segment, it undoubtedly achieves better quantitative results in terms of both IoU value and overall segmentation accuracy. Therefore, the deep learning framework in this study, combined with the local features of tree point clouds, further improves the accuracy and performance of segmentation.

[0153] Table 2 shows the point cloud branch and leaf classification detection results using different methods:

[0154] method Overlap (branches / leaves) / % Overall accuracy / % Classification time (training / testing) / h Traditional machine learning methods 60.52 / 54.28 75.24 25.45 / 0.04 PointNet 70.16 / 66.53 82.75 43.94 / 0.06 This research method 78.24 / 72.74 90.32 48.14 / 0.06

[0155] 2.2 Results of rubber tree framework reconstruction:

[0156] Based on the skeleton data extracted by deep learning, and according to the tree skeleton reconstruction method in 1.4, tree skeleton reconstruction was carried out on three different types of rubber tree test plots. Figure 6 The results of branch and leaf separation and tree skeleton reconstruction were shown in three rubber tree plots. Figure 6 (a1) is a schematic diagram of laser point cloud data and branch and leaf separation results of the test plot of rubber tree variety PR107.

[0157] Figure 6(a2) shows the laser point cloud data and branch-leaf separation results of the test plot of rubber tree variety CATAS 7-20-59. Figure 6 The middle (a3) ​​shows the laser point cloud data and branch and leaf separation results of the test plot of rubber tree variety CATAS 8-79. Figure 6 The middle (b1) diagram shows the fitting results of the cylindrical branches of the test plot of rubber tree variety PR107. Figure 6 Figure (b2) shows the fitting results of the cylindrical branches of the test plot of rubber tree variety CATAS 7-20-59. Figure 6 Figure (b3) shows the fitting results of the cylindrical branches of the test plot of rubber tree variety CATAS 8-79. Figure 6 The middle (c1) diagram shows the classification results of the main branches and primary branches of the rubber tree variety PR107 in the test plot. Figure 6 The middle (c2) diagram shows the classification results of the main branches and primary branches of the rubber tree variety CATAS 7-20-59 in the test plot. Figure 6 The middle (c3) diagram shows the classification results of the main branches and primary branches of the rubber tree variety CATAS 8-79 in the test plot. Figure 6 The middle (d1) plot shows the leaf clustering results and leaf cluster convex hull display of the test plot of rubber tree variety PR107 based on spatial watershed and first-order branching. Figure 6 The middle (d2) plot shows the leaf clustering results and leaf cluster convex hull display of the test plot of rubber tree variety CATAS 7-20-59 based on spatial watershed and first-order branching. Figure 6 The middle (d3) plot shows the leaf clustering results and convex hull display of leaf clusters based on spatial watershed and first-order branching in the test plot of rubber tree variety CATAS 8-79.

[0158] Figure 6 Figures (b1), (b2), and (b3) show the cylinder fitting results for the three sample plots. During the fitting process, the branch point cloud was first divided into multiple layers from bottom to top at certain height intervals (0.3m). The center points of each layer were extracted through clustering, and connected chains were established based on the shortest Dijkstra distance. Then, cylinders with directionality and adaptive radius were used to perform segmented fitting of the branches. After fitting, each rubber tree branch was composed of approximately 180 cylinders on average, and each cylinder was labeled with a different color.

[0159] Secondly, based on the connection chain between the growth angle and the center point, the tree skeleton is divided into main branches and many first-order branches, such as... Figure 6As shown in (c1), (c2), and (c3), the main branches were extracted, and the various first-order branches extending from the main branches are represented by different shades of color. It can be seen that the PR107 rubber tree exhibits greater variation in branch growth angle, with branches spreading outwards, resulting in a radiating tree structure. Compared to the PR107 rubber tree, the CATAS 7-20-59 and CATAS 8-79 rubber trees show less variation in branch growth angle, with branches extending mostly upwards, resulting in a vase-shaped tree structure.

[0160] Subsequently, based on plant physiological theory, leaf clusters were defined as sets of leaves growing on the same main stem or first-order branch. Using a spatial three-dimensional watershed segmentation method, the leaf clusters of each tree were segmented using the extracted point clouds corresponding to the main stem and first-order branches as cluster centers. The segmented leaf clusters were then represented by different colors. The segmentation results are shown below. Figure 6 As shown in (d1), (d2), and (d3), among the three rubber tree plots, the PR107 variety has more primary branches, resulting in approximately 3 leaf clusters per tree; while the CATAS 7-20-59 and CATAS 8-79 varieties have fewer primary branches, resulting in an average of about 2.5 leaf clusters per tree. Finally, the convex hull algorithm was used to mark the leaf clusters of each tree and calculate their volume.

[0161] 2.3 Analysis of Rubber Tree Parameter Inversion Results:

[0162] Table 3 presents the actual measured values ​​and calculated values ​​of growth parameters for trees in three rubber tree test plots (PR107, CATAS 7-20-59, and CATAS 8-79), including diameter at breast height (DBH), diameter of primary branches, and the angle between the main branch and the primary branch. Additionally, the values ​​are compared using three comparative indicators: R0... 2 The effectiveness of the method in this study was analyzed by quantifying the root mean square error (RMSE) and the relative root mean square error (rRMSE). Figure 7 Specific parameter comparison results are provided.

[0163] Table 3 compares the stand parameters obtained by the method used in this study with the field measurements:

[0164]

[0165] Note: (F) represents the actual measured value, and (O) represents the calculated value using the methods of this study.

[0166] Figure 7 (a) shows the comparison between the measured and calculated diameter at breast height (DBH) values ​​of rubber trees in the three rubber tree plots obtained using the methods of this study. The comparison results between the PR107 and CATAS 7-20-59 rubber trees and the field observations are as follows: (R...) 2=0.95, RMSE=3.01cm, rRMSE=12.22%) and (R 2 =0.94, RMSE=2.78cm, rRMSE=12.53%. Compared with the first two plots, the algorithm achieved better results in extracting tree branch parameters of CATAS 8-79 rubber tree (R = 0.94, RMSE = 2.78cm, rRMSE = 12.53%). 2 =0.97, RMSE=2.63cm, rRMSE=10.55%. This is mainly because the point cloud data of PR107 and CATAS 7-20-59 rubber trees were collected during the leafy stage. The mutual occlusion between different organs of the tree led to data loss, and the wind disturbance during the scanning process caused the tree to shake, generating noise, which in turn affected the final parameter inversion results. In contrast, the CATAS 8-79 rubber tree is not resistant to cold damage, and the data was collected in the leafless state. Therefore, the obtained branch and trunk data were more complete. In addition, the branches and trunks of the CATAS 8-79 rubber tree are thicker, and the impact of tree vibration caused by wind disturbance is smaller, resulting in higher quality point clouds, less noise, and higher accuracy of parameter inversion.

[0167] Figure 7 Figure (b) shows the comparison between the calculated angles between the main branches and first-order branches of the trees in three rubber tree plots and the measured values. It can be seen that good accuracy was obtained for different varieties of rubber trees. Among them, the PR107 variety (R... 2 =0.94, RMSE=4.94°, rRMSE=11.23%) and CATAS 8-79 variety (R 2 =0.92, RMSE=3.48°, rRMSE=10.86%) correlation coefficient relative to CATAS 7-20-59 variety (R 2 The slightly higher values ​​(rRMSE = 0.91, RMSE = 3.81°, rRMSE = 11.45%) can be explained by the fact that the PR107 rubber tree variety has clearer branch characteristics and a divergent skeletal structure, with a larger angle between the primary branches and the main branches, making it easier to obtain accurate results in field measurements and algorithmic angle calculations. In contrast, the sparse foliage in the CATAS 8-79 rubber tree forest results in less shading, making the branch angle estimation more accurate.

[0168] Figure 7Table (c) shows the correlation between the width of primary branches and the corresponding leaf cluster volume in the three rubber tree plots. Among them, the PR107 variety of rubber tree, due to its larger branching angle, provides more space for leaf growth, resulting in a larger crown volume and leaf cluster volume corresponding to the primary branches. The CATAS 7-20-59 variety of rubber tree has a smaller branching angle range, forming a vase-shaped crown, thus its leaf clusters are of medium size and evenly distributed. The CATAS 8-79 variety of rubber tree, due to frost damage and leaf drop, has smaller leaf cluster volume corresponding to the primary branches. Figure 7 The distribution is relatively concentrated in the middle. Figure 7 (c) shows that the width of the primary branch and the corresponding leaf cluster volume generally satisfy a positive correlation distribution, that is, the thicker branch supports the larger leaf cluster. This is consistent with the principle that plant branches are responsible for nutrient transport and the gravity bearing of organs such as leaves, flowers and fruits. The resulting local deviation can be explained by local changes caused by human pruning or natural environmental disturbances.

[0169] 3. Conclusion:

[0170] Accurate reconstruction of tree skeleton models plays a crucial data support role in analyzing tree phenotypic structure, phototropism, and the influence of forest competition factors. First, this study proposes an improved PointConv deep learning algorithm for backpack-based mobile LiDAR data from rubber plantations. By performing voxel subdivision of the point cloud and designing a four-layer encoding and decoding layer, along with an improved PointConv module and interpolation feature propagation module, the algorithm automatically acquires the global optimum and local neighborhood kernel density features of the forest point cloud at different scales, aiming to accurately classify branch and leaf point clouds. Analysis of test data from three different rubber plantation varieties shows that the classification accuracy of the deep learning method in this study reaches 90.32%. Second, for the extracted branch point cloud, computer graphics algorithms are used to automatically extract forestry parameters such as main branch diameter at breast height (DBH), first-order branch diameter, branch angle, and leaf cluster volume. Finally, by comparing with measured values ​​from the sample plots, the correlation coefficient R between the calculated and actual values ​​of the main branch and first-order branch diameters of different rubber tree varieties is shown. 2 All values ​​were above 0.94, and the RMSE was below 3.01 cm; the correlation R between the estimated and actual values ​​of the angle between the main branch and the first-order branch was [value missing]. 2The RMSE is above 0.91 and below 4.94°, and the leaf clusters of rubber trees show an approximately positive correlation with the diameter of the corresponding first-order branches. In ecological studies of forest trees, the apparent differences between different species of the same family and genus, the spatial allometric growth patterns they exhibit, and the productive driving forces generated by their interaction with the natural environment can all be found in the apparent morphological structure of the tree. The deep learning network and computer graphics algorithm proposed in this paper can automatically reconstruct the skeletal spatial model of trees from forest tree scanning data, and is applicable to parameter inversion of different varieties of rubber trees. This represents one of the current applications of artificial intelligence technology in the forestry field.

[0171] The scope of protection of this invention includes, but is not limited to, the above embodiments. The scope of protection of this invention is defined by the claims. Any substitutions, modifications, or improvements to this technology that are easily conceived by those skilled in the art fall within the scope of protection of this invention.

Claims

1. A method for separating branches and leaves from a laser point cloud of a forest based on a deep learning method, characterized in that: The method comprises the following steps: Step 1: acquiring forest laser point cloud data; Step 2: performing de-noising processing on the collected laser point cloud data, and dividing the de-noised laser point cloud data into above-ground points and ground points; Step 3: performing branch and leaf separation operation on the above-ground points by using a machine learning algorithm and manual correction method, and voxelizing and sectioning the above-ground points that have undergone the branch and leaf separation operation; Step 4: taking the point cloud data that have undergone the branch and leaf separation operation and voxelization and sectioning as a training sample data set, and expanding the training sample data set by using a data augmentation method to obtain a new training sample data set; Step 5: constructing a deep learning network, training the deep learning network by using the new training sample data set, and obtaining a trained deep learning network model; Step 6: collecting laser point cloud data of a forest to be measured, performing de-noising processing on the collected laser point cloud data, dividing the de-noised laser point cloud data into above-ground points and ground points, voxelizing and sectioning the above-ground points, and inputting the point cloud data in a single voxel into the trained deep learning network model to realize branch and leaf separation of the point cloud data in the voxel. The deep learning network comprises a feature encoding layer and a feature decoding layer; The feature encoding layer is configured to extract point cloud feature information by using a down-sampling method, and the feature decoding layer is configured to transmit point cloud feature information by using an up-sampling method; The feature encoding layer comprises a sampling module, a first grouping module and a first PointConv module; The feature decoding layer comprises a feature interpolation module, a second grouping module and a second PointConv module; There are four feature encoding layers in total, and there are four feature decoding layers in total; The calculation process of the sampling module in the jth feature encoding layer is as follows: Input a point cloud matrix with size (d+C j )×N j , select N j+1 sub-sampling points from the point cloud matrix by the farthest point sampling method to obtain a point cloud matrix with size (d+C j+1 )×N j+1 ; Where: (d+C j )×N j It refers to having d-dimensional coordinates and C j N-dimensional features j Point, (d+C) j+1 )×N j+1 It refers to having d-dimensional coordinates and C j+1 N of high-dimensional features j+1 Individual sampling points; The calculation process of the first grouping module in the jth feature encoding layer is as follows: The input sampling module outputs d x N j+1 sub-sampling points of the coordinate and the initial input (d+C j ) x N j feature point set, using the nearest neighbor rule classification method, in the neighborhood of each point in the N j+1 sampling points extracted by the sampling module, find the nearest K neighbors to form a small group, get d x K x N j+1 neighbor point coordinate index set relative to each sampling center point, and get (d+C j ) x K x N j+1 neighbor point feature set according to the coordinate index, wherein K is the number of neighbor points in the neighborhood of each point in the N j+1 sampling center points. The input of the first PointConv module in the jth feature encoding layer is three parts of data; wherein the first part of data is: the features of K neighboring points near the N j+1 center points, with a size of (d+C j )×K×N j+1 ; the second part of data is: the local coordinates of K neighboring points near the N j+1 center points relative to the center points, with a size of d×K×N j+1 ; and the third part of data is: the density of K neighboring points near the N j+1 center points, with a size of 1×K×N j+1 ; the first PointConv module outputs a new local region with a size of (d+C j+1 )×N j+1 ; The calculation process of the first PointConv module in the jth feature encoding layer is as follows: Step (1), calculating the center point by using the kernel density estimation method The kernel density estimation of each input point in the K nearest neighbor points nearby the center point, and then the inverse density transformation of the K nearest neighbor points is calculated by the kernel density estimation. The result of the inverse density transformation is input into the multi-layer perception MLP2 to perform the nonlinear transformation of the activation function, and the inverse density coefficient with the size of 1×K is obtained. The inverse density coefficient is copied and tiled for C in times, and the inverse density tensor with the size of C in ×K is output. Step (2), input the local coordinates of K nearest neighbors relative to the center point into a multi-layer perceptron MLP1, output a weight tensor W of size (C in × C out ) × K, where C in is the number of input features and C out is the number of output features; Step (3): Set the output of step (1) to C. in The inverse density tensor of ×K and the feature F of the K nearest neighbors in Perform a dot product operation, and the output size is C. in ×K matrix Again Perform C out This tiling operation, repeated 10 times, outputs a size of (C). in ×C out The intermediate tensor of )×K; Step (4), performing a dot product operation on the weight tensor W of size (C in × C out ) x K output from step (2) and the intermediate tensor of size (C in × C out ) x K output from step (3), and outputting a matrix of size (C in × C out ) x K Step (5), summing the first and last dimensions of the matrix obtained in step (4) to obtain the output feature F of the center point out ;​ wherein center point is N j+1 a point in the center point.

2. The forest laser point cloud branch and leaf separation method based on the deep learning method according to claim 1, wherein: The calculation process of the feature interpolation module in the jth feature decoding layer is as follows: N j sample points obtained by the sampling module in the (j-1)th feature coding layer as the center, select three neighboring points closest to each of the N j+1 sample points obtained by the sampling module in the jth feature coding layer j sample points obtained by the sampling module in the jth feature coding layer sample points obtained by the sampling module in the jth feature coding layer sample points obtained by the sampling module in the jth feature coding layer sample points obtained by the sampling module in the jth feature coding layer j+1 sample points obtained by the sampling module in the jth feature coding layer The three neighborhood points are calculated according to the reciprocal of the distance The weight value of the corresponding center point in the N j sampling points The calculation formula is: ​ wherein represents the Euclidean distance between points and μ represents the weighted power exponent; The three neighboring points are indexed by the index of the initial input (d+C')xN j )×N j+1 The three neighboring points in the feature point set are obtained (d+C')x3, and the corresponding weight value is combined The features of the three neighboring points are summed and averaged to obtain The corresponding (d+C')x1 interpolation feature j )×1 interpolation feature The calculation formula of the interpolation feature is Because the j-1th feature encoding layer has N j sampling points in total, the interpolated feature of (d+C′ j )×N j is obtained. obtained (d+C′ j )×N j )×N j )×N j )×N j )×N j )×N j +C′ j )×N j ; Where (d+C′) j )×N j This indicates that there are d-dimensional coordinates and C′. j N-dimensional features j One point.

3. The forest laser point cloud branch and leaf separation method based on the deep learning method according to claim 2, wherein: The calculation process of the second grouping module in the jth feature decoding layer is as follows: the final output of the input feature interpolation module (2d+C j +C′ j )×N j and the coordinates of the d×N j sub-sampled points obtained by the sampling module in the j-1th feature encoding layer, using the nearest neighbor rule classification method, find the nearest K' neighbors in the neighborhood of each point in the N j sampled points to form a small group, obtain the d×K'×N j coordinate index set of the nearest neighbor points relative to each sampled center point, and obtain the feature set of the (2d+C j +C′ j )×K'×N j nearest neighbor points according to the coordinate index, where K' is the number of nearest neighbor points in the neighborhood of each point in the N j sampled center points.

4. The forest laser point cloud branch and leaf separation method based on the deep learning method according to claim 3, wherein: The input of the second PointConv module in the jth feature decoding layer is three parts of data; wherein the first part of data is: the features of K' neighboring points near the N j center points, with a size of (2d+C j +C' j )×K'×N j ; the second part of data is: the local coordinates of K' neighboring points near the N j center points relative to the center points, with a size of d×K'×N j ; and the third part of data is: the density of K' neighboring points near the N j center points, with a size of 1×K'×N j ; the second PointConv module outputs a new local region with a size of (d+C' j-1 )×N j .

Images