1. Introduction

2. Materials and Methods

2.1. Materials

2.1.1. Study Area

The Northeast China Hupao National Park is located in the northeastern region of China, along the border of Jilin Province and Heilongjiang Province, at 42°31′06″N-44°14′49″N, 129°5′0″E-131°18′48″E. The park is situated in the southern region of the Changbai Mountain Range’s Laoyeling branch, with the terrain mainly consisting of middle and low mountains, canyons and hills, reaching an elevation of 1477.4 meters at Laoyeling. The park’s soil is predominantly dark brown soil and marshy soil. The park is located at the heart of the temperate coniferous and broad-leaved mixed forest ecosystem in Asia, with a continental humid monsoon climate, boasting an extremely rich variety of temperate forest plant species. The forest coverage rate is 93.32%, with the main vegetation type being temperate coniferous and broad-leaved mixed forests.

Figure 1. Images and point cloud data of eight sample plots in Northeast China Hupao National Park. The red circles represent the plots areas.

Figure 1. Images and point cloud data of eight sample plots in Northeast China Hupao National Park. The red circles represent the plots areas.

2.1.2. Sample Plots Data

In the study area, eight circular sample plots with a diameter of 30 meters are established, and various forest parameters for each tree are measured within the plots. These parameters include tree height, north-south canopy width, east-west canopy width, height under branches (HUB), DBH, tree species, tree location and three indirectly measured parameters: stock volume, above-ground biomass (AGB), and below-ground biomass (BGB). For each sample plot, the average value or the cumulative value of the measured parameters are recorded. As these are natural forests, the trees vary in height and canopies overlap. The terrain, tree species distribution, and tree growth conditions are all different for the eight sample plots, which meets the requirements for sample plots verification and updating tests.

Additionally, the manually obtained information from the sample plots is confirmed in the outfield, and the confirmed data is considered accurate, so it can be used to compare with the estimated results for accuracy.

Table 1. Parameters for the eight sample plots.

Table 1. Parameters for the eight sample plots.

Plot ID		1	2	3	4	5	6	7	8
Number of Trees		50	42	45	41	29	41	31	39
Dominant Tree Species		pine	linden	linden	linden	linden	poplar	oak	linden
Average Tree Height (m)		14.5	13.5	16.9	8.2	14.3	16.4	13.8	14.3
Average DBH (cm)		18.6	17.8	16.9	10.5	20.7	20.8	20.1	18.3
Average HUB (m)		4.8	3.6	3.2	2.3	4.0	5.1	3.6	4.6
Average Canopy Width (m)	East-West	2.72	2.49	2.61	2.27	2.60	2.83	2.82	2.27
	North-South	2.72	2.10	2.31	2.16	1.82	2.97	2.20	1.98
Cumulative Stock Volume (m3)		14.4	9.7	8.1	3.5	9.5	14.6	9.8	9.1
Cumulative AGB (t/ha)		8390.7	7743.4	5879.0	2239.1	7158.5	8565.2	9342.1	7717.1

2.1.3. LiDAR Data

The LiDAR data is acquired using an airborne LiDAR scanning system, combined with IMU/DGPS-assisted surveying technology. Specifically, a Cessna 208b aircraft is selected to carry the riegl-vq-1560i LiDAR payload (with seamless integration of an Inertial Measurement Unit (IMU) and global navigation satellite system (GNSS)). The Cessna 208b aircraft also jointly observes data with ground Cooperating Operating Reference Stations (CORS) and artificial base stations. The LiDAR data acquisition flight altitude is set to an average height above ground level (AGL) of 1000 meters, with the laser recorder capturing four echoes for each laser pulse. With a flight line overlap of 20%, the average laser point density is 20 points per square meter (ppm²). In accordance with the project’s task specifications, the reported horizontal accuracy is between 15 and 25 centimeters (cm), while the vertical accuracy is approximately 15 cm.

To conduct the experiment, the point cloud data is clipped using vectors of the sample plots boundaries, with an additional buffer zone of 30m radius reserved. The resulting point cloud data within the boundaries is used for the experiment.

Figure 2. Point cloud data of the eight sample plots after clipping.

Figure 2. Point cloud data of the eight sample plots after clipping.

2.2. Methods

Figure 3 illustrates the workflow of this experiment, which consists of two parts: data preparation and sample plot verification. The data preparation includes the pre-processing of the point cloud data and the segmentation of trees in the plots. The sample plot verification involves plots matching, tree species identification and extraction of forestry parameters.

2.2.1. Point Cloud Pre-processing

Pre-processing includes filtering and normalization of point cloud data, aiming to obtain a point cloud that eliminates terrain effects and reflects the real horizontal and vertical structure of the forest.

Point cloud filtering is the process of separating ground points and non-ground points (mainly tree points), and the obtained ground points can be used to generate a digital elevation model (DEM). Point cloud filtering is implemented using the LiDARForest software, which uses the progressive triangulated irregular network (TIN) densification algorithm (PTD) proposed by Axelsson [43]. The main idea of this algorithm is to select some low points as seed points to construct an initial TIN, and then gradually judge other points. If a point is a ground point, it is added to the TIN, otherwise it is filtered out. The point is classified as a ground point when its angle and vertical distance to the endpoints on the triangle plane are below the set thresholds. Repeat the above process until there are no more points that meet the threshold conditions. The experimental results show that point cloud filtering based on the PTD algorithm has fewer iterations and is suitable for most terrain conditions.

The elevation values of LiDAR point cloud are a combination of ground elevation and object elevation, which makes it difficult to reflect the actual vegetation height. Therefore, before conducting relevant research, it is common to normalize the point cloud. Point cloud normalization refers to the process of subtracting the elevation values of DEM from point cloud, which is generated by interpolating the ground points obtained from filtering. The normalized point cloud CNTP eliminates the influence of terrain and reflects the true horizontal and vertical distribution of trees. The highest point of the segmented individual tree can be used as the tree height.

Figure 4. Point cloud pre-processing procedure. (a) Original point cloud; (b) Filtered point cloud; (c) Normalized point cloud. Ground points are shown in brown and tree points in green in (b) and (c).

Figure 4. Point cloud pre-processing procedure. (a) Original point cloud; (b) Filtered point cloud; (c) Normalized point cloud. Ground points are shown in brown and tree points in green in (b) and (c).

2.2.2. Tree Segmentation

To implement the tree segmentation, the positions of the trees within the plots need to be determined. Based on the structure of the trees, the position of the highest point within individual tree range can be selected as the center of the tree. Therefore, the local highest points in C_NTP are selected as the initial treetop points.

The C_NTP is partitioned into grids, with the grid size, d, depending on the average distance between points. The elevation value Z of the highest tree point within each grid is assigned to the grid. If there are no tree points within the grid, it is assigned a value of 0. The treetop points are selected from the C_NTP according to the following criteria: the grid containing the treetop point should have the maximum value within its neighborhood, and it should not be adjacent to any unit that does not contain any tree points. To prevent the loss of short treetop points, the grid size d is generally set to an acceptable small value, which results in some short treetop points actually being redundant parts of adjacent taller trees. To address this, a distance threshold D is defined, and the shorter treetop points are removed when the distance between two treetop points is less than D. Finally, the collection of all the selected treetop points forms the initial treetop point set C_TOP.

Figure 5. Workflow for obtaining initial treetop points.

Figure 5. Workflow for obtaining initial treetop points.

After obtaining the initial treetop points, many studies use a clustering method to segment the canopies. This study focuses on the vertical profile structure of initial points and utilizes the rotational profiles to search for the edges of the individual tree canopy, thereby performing the tree segmentation and refining the positions of the trees [44].

As shown in Figure 6, for a treetop point T₀, the longitudinal profile at a certain angle may have the following cases: (a) the profile only contains the canopy of T₀; (b) there are N treetop points T_i (1 ≤ i ≤ N) in addition to T₀ within the profile, which is the general situation. The profile is divided into M sub-segments S_j (1 ≤ j ≤ M) of equal length and Z_j is the maximum Z value of the tree points in the sub-segment S_j, the corresponding tree point is P_j. The edge point P_EDGE and the intersection point P_CROSS are searched in P_j. For situation (a), where only the canopy of T₀ exists in the profile, only P_EDGE needs to be searched. In this case, the edge point is located at both ends of the profile, as shown in Figure 6. For situation (b), in addition to the edge point, the intersection point P_CROSS between the canopy of T₀ and other canopies needs to be found. According to the overlapping situation of tree canopies, the intersection point is regarded as a valley point in the profile. Therefore, Equation 1 for obtaining P_EDGE and P_CROSS are derived.

$$P_j=\left\{\begin{array}{c}{P}_{EDGE} , Z_j\in min\left(Z\right) \\ P_{CROSS} , Z_j<Z_{j-1},Z_j<Z_{j+1}\end{array}\right.$$

(1)

where P_j: the vertex of the subsection S_j of the profile, (1 ≤ j ≤ m-1)

Z_j: Z value of P_j

Figure 6. Sub-segments and the location of P_EDGE and P_CROSS. (a) Single canopy; (b) Multiple canopies.

Figure 6. Sub-segments and the location of P_EDGE and P_CROSS. (a) Single canopy; (b) Multiple canopies.

The profile varies in different directions. To accurately identify individual trees, the edges can be extracted by analyzing a series of rotational profiles, which are divided at intervals according to a certain angle. The steps are as follows:

(1)

Arrange C_TOP in descending order of elevation to obtain the sequence of treetop points T_i;

(2)

Get a profile at the given angle for treetop T₁. Let Z_j be the maximum Z value of the tree points in sub-segment S_j, and P_j is the corresponding tree point. Find P_EDGE and P_CROSS in P_j (Figure 6 and Equation 1) to obtain the canopy range under this profile;

(3)

Rotate the profile around T₁ to another angle to obtain the new P_EDGE and P_CROSS. Repeat the rotation until the rotation angle reaches 180°, then all possible edges are obtained, as well as the rough range Rt of the tree;

(4)

Iterate over the remaining treetop points, removing the point from C_TOP if it lies within Rt. All high tree points within Rt can be grouped into the tree’s canopy;

(5)

Select the next treetop in descending order of elevation and repeat (2) - (4) until all treetop points have been segmented;

(6)

High points within a linear distance less than D_MAX from the tree’s axis are grouped into the corresponding canopy set. The tree’s axis is the line from the canopy’s center of gravity to the ground, and this method can exclude some isolated high points.

Figure 7. Workflow of rotational profile algorithm.

Figure 7. Workflow of rotational profile algorithm.

Figure 8. Treetop points location after segmentation. (a) Original point cloud; (b) Segmented point cloud; (a) A partial profile of (b). The treetop points are shown in rectangles or circles in (b) and (c).

Figure 8. Treetop points location after segmentation. (a) Original point cloud; (b) Segmented point cloud; (a) A partial profile of (b). The treetop points are shown in rectangles or circles in (b) and (c).

2.2.3. Plots Matching Based on Dominant Trees

The spatial accuracy of measured data is limited by the measurement environment, while point clouds are more accurate (see Table 2 for accuracy). The measured data and point cloud data do not match in spatial position, so plot matching is required to register the measured data to the point cloud data, achieving spatial consistency between the two data sets before further operations can be carried out.

The concept of matching degree R is proposed for plots matching: Δh is the difference between the measured tree height and the point cloud height. When the absolute value of Δh is less than the threshold D_T, the position of the sample tree is considered correct. The matching degree of the plot is defined as the ratio of the number of trees that satisfy Δh < D_T to the total number. When the matching degree of the plot is greater than 90%, the spatial position of the site can be considered accurate.

$$R=\frac{n}{N}\times 100\%$$

(2)

where R: the matching degree

N: the number of trees in the plot

n: The number of trees in the plot that satisfy Δh < D_T

When the matching degree does not meet the requirement, it is necessary to perform plots matching based on dominant trees. Dominant trees in forestry refer to the largest trees in terms of diameter, tree height, tree canopy height, occupying the largest space, receiving the most sunlight and experiencing almost no compression. Based on the positive correlation between tree height, diameter and tree canopy height, we define the top 10% of tree heights (or 5% when there are many trees) in the plot as dominant trees. Using the position information of dominant trees, a translation and rotation coordinate transformation matrix is established between the plot and the LiDAR point cloud data. The plot is rotated and translated so that the position of the dominant trees in the plot matches that of the segmented point cloud, and plot attributes are updated when the matching degree exceeds 90%.

In fact, the positions of the dominant trees in measured data and the segmented point cloud cannot achieve complete consistency. This is because of the presence of inclines and tree canopies, which make the treetop position not fully represent the tree center. Therefore, deviation value D (Equation 3) is proposed to quantitatively reflect the deviation between the positions of the dominant trees in the measured data and the segmented point cloud. It is required that the matching degree R is greater than 90%, and the deviation value D is also smaller than the set threshold. At this point, the matching of the two data sets is considered complete, and can be used for subsequent parameter extraction.

$$D=\frac{{\sum }_{i=1}^n∆d_i}{n}$$

(3)

where D: deviation value

n: the number of the dominant trees in the plot

Δd_i: deviation of the position of the ith dominant tree in measured data and point cloud. (1 ≤ i ≤ n)

Figure 9. Plots matching based on dominant trees. (a) Locations of the measured trees and extracted trees (from segmented point cloud) before matching; (b) Locations after matching.

Figure 9. Plots matching based on dominant trees. (a) Locations of the measured trees and extracted trees (from segmented point cloud) before matching; (b) Locations after matching.

2.2.4. Tree Species Identification and Verification

Traditional field surveys are currently the most reliable way to obtain tree species information. However, errors in tree species identification and recording are inevitable. Nonetheless, the majority of tree species information obtained through field surveys is correct. Therefore, this paper utilizes correct tree species information obtained from field surveys to automatically identify and verify tree species information in the sample plots.

Currently, many scholars use the combination of canopy point cloud and spectral information to achieve tree species identification [45, 46]. However, the spectral information of a certain tree species changes with season, time, and leaf age. Therefore, in addition to deformation caused by occasional lightning strikes or rock invasion, the geometric shape of the tree canopy is the more reliable information for tree species identification. Different tree species often have different tree canopy shapes, as shown in Figure 10. In this paper, we sampled the canopy width of typical species according to height to represent the geometric shape of the canopy profile. The obtained width set is used as a sample to train the DBN classifier, and the specific details are as follows:

(1)

Extraction of key points based on parallel lines

For a given set of segmented tree points of a single tree, draw a vertical line L from the treetop to the ground. Choose several profiles with a certain width D_W centered on the treetop and use the selection criterion that they should not be connected to other tree canopies. The size of D_W is related to the average point distance dis of the point cloud, generally set as D_W = 2dis. On a profile, draw parallel lines along L with a certain spacing Δd to obtain the positions of endpoints within the parallel lines space. After all parallel lines are scanned all endpoints can be defined as key points of the canopy under the profile.

Figure 11. Different profiles of the segmented point cloud from a single tree.

Figure 11. Different profiles of the segmented point cloud from a single tree.

Figure 12. Extraction of key points based on parallel lines. The yellow points are the key points.

Figure 12. Extraction of key points based on parallel lines. The yellow points are the key points.

(2)

Extraction of geometric morphological features of canopy profiles

Connect key points on the same parallel lines, record the position l_i and length gi of obtained intersection (1 ≤ i ≤ n, the number of the lines is n+1). The intersection positions and lengths obtained from parallel lines can summarize the shape of the tree canopy, but the number of intersection lines for each tree may not be the same due to the influence of tree height. Due to the limitation of point cloud density, the length accuracy is reduced when the same number of intersection lines are forced to be intercepted for each tree. To facilitate the subsequent training of the DBN, the cubic spline interpolation method is used to interpolate the number of intersection lines for each tree to be consistent. The cubic spline interpolation method is highly accurate and can avoid the Runge phenomenon, which causes violent oscillations in the interpolation function curve during polynomial interpolation.

S(l) is set as the interpolation function. Each small interval of S(l) is a cubic function, as shown in Equation 4 and Equation 5:

$$S\left(l\right)=\left\{s_i\left(l\right),l\in \left[l_{i-1},l_i\right],i=1,\cdots ,n\right\}$$

(4)

$$s_i\left(l\right)=a_i{l}^3+b_i{l}^2+c_{i}l+d_i , \left(\mathrm{i}=1,\cdots ,\mathrm{n}\right)$$

(5)

The function S(l) should satisfy the value of g_i at node l_i, as shown in Equation 6:

$$S\left(l_i\right)=g_i , \left(i=1,\cdots ,n\right)$$

(6)

At all nodes (except the first node ${\mathrm{l}}_0$ and last node${\mathrm{l}}_{\mathrm{n}}$), the first-order derivative and second-order derivative are continuous, ensuring the same slope and degree of curvature at the nodes, as shown in Equation 7:

$$\left\{\begin{array}{c}{s’}_i\left(l_i\right)={s’}_{i+1}\left(l_i\right)\\ {s’’}_i\left(l_i\right)={s’’}_{i+1}\left(l_i\right)\end{array}\right. \left(\mathrm{i}=1,\cdots ,\mathrm{n}-1\right)$$

(7)

The natural spline condition: the second derivative of the first node and last node is 0, allowing the slope of the endpoints to remain balanced at a certain position and minimizing the oscillation of the interpolation function curve. Equation 8:

$$S’’\left(l_0\right)=S’’\left(l_n\right)=0$$

(8)

Based on Equations 6, 7, and 8, the coefficients a_i, b_i, c_i and d_i in Equation 5 can be solved to determine the interpolation function S(l). Then, by reselecting appropriate positions l’_i and substituting them into S(l), the corresponding intersection lengths g’_i can be obtained. This process yields a set of lengths with equal dimensions for each profile, which can be used as a feature vector to describe the geometric characteristics of the profile.

(3)

Learning and updating of tree species by DBN

Deep belief network (DBN) [47] is chosen for learning and updating of tree species. Compared to most other deep learning models, DBN is more thoroughly trained on data, achieving better convergence in small-scale samples and making it more applicable in forest environments with limited sample sizes.

The DBN is formed by stacking multiple Restricted Boltzmann Machines (RBMs), whose structure is shown in Figure 13. The RBM consists of two layers: a hidden layer and a visible layer. When stacking RBMs to form a deep learning network, the output layer (hidden layer) of the previous RBM serves as the input layer (visible layer) of the next RBM, and so on. Finally, an output layer is added to form the basic structure of the DBN.

RBM is a fully connected layer with no connections within each layer, meaning that neurons in the same layer are independent of each other. Set $v$ as the network input; the neuron likelihood functions of the visible and hidden layer can be expressed as follows:

$$P\left(\mathit{h}|\mathit{v}\right)=\prod _{j=1}^{N}P\left(h_j|\mathit{v}\right)$$

(9)

$$P\left(\mathit{v}|\mathit{h}\right)=\prod _{i=1}^{M}P\left(v_j|\mathit{h}\right)$$

(10)

where $\mathit{v}={\left[v_1,v_2,\dots ,v_M\right]}^{\mathrm{T}}$ , is the visible layer input vector

$\mathit{h}={\left[h_1,h_2,\dots ,h_N\right]}^{\mathrm{T}}$ , is the hidden layer input vector

M and N are the number of neurons in the visible and hidden layers, respectively

The energy function between the visible layer and the hidden layer is defined as follow:

$$E\left(\mathit{v},\mathit{h};\theta \right)=-\sum _{i,j=1}^{M,N}{\omega }_{i,j}{v}_i{h}_j-\sum _{i=1}^M{b}_i{v}_i-\sum _{j=1}^N{c}_j{h}_j$$

(11)

where ${\omega }_{i,j}$: weight between ith visible layer neuron and jth hidden layer neuron

$b_i$: bias of ith visible layer neuron

$c_j$: bias of jth hidden layer neuron

$\theta =\left\{\mathit{W},\mathit{b},\mathit{c}\right\}$, are the network parameters, obtained by training

The activation functions of visible layer and hidden layer neurons are:

$$P\left(h_j|\mathit{v}\right)=Sigmoid\left(c_j+\sum _i{\omega }_{i,j}{v}_i\right)$$

(12)

$$P\left(v_i|\mathit{h}\right)=Sigmoid\left(b_i+\sum _j{\omega }_{i,j}{h}_j\right)$$

(13)

$$Sigmoid\left(x\right)=\frac{1}{1+e^{-x}}$$

(14)

The DBN is formed by stacking multiple RBMs, and Figure 14 shows the DBN structure. During the training process, the previous layer’s RBM is fully trained before the current layer’s RBM is trained, repeating until the final layer where a Softmax classifier is added to form a complete DBN structure.

In this study, the set of interpolated cubic spline intersection lengths g’_i are used as input, while the corresponding tree species are used as the target output to train the DBN to achieve the profile-based tree species identification. According to the experiment, this study uses a DBN with three layers of RBMs for the identification, which shows the fastest convergence speed and the highest accuracy among all models tested.

Figure 15. Tree species information extraction and checking process.

Figure 15. Tree species information extraction and checking process.

2.2.5. Forestry Parameters Extraction

The segmented point cloud reflects the structure of individual tree canopies, and some forestry parameters are directly reflected in it. For example, the elevation of the treetop point can be used as the tree height. East-west canopy width can be obtained based on the maximum distance in the east-west direction of the canopy point set of an individual tree and the north-south canopy width can be obtained similarly. Some parameters may not be directly obtained, but they can be indirectly calculated through estimating equations using known tree species, tree height, canopy width and other information extracted from the segmented point cloud as parameters. The methods for obtaining forestry parameters and estimating equation factors are shown in Table 3.

DBH refers to the diameter of a tree trunk at a height of 130cm above the ground, and can be estimated by inputting the tree height into the model. Forest biomass is the energy foundation and material source of the forest ecosystem; it is divided into AGB and BGB, usually expressed in terms of dry matter or energy accumulation per unit area or unit time. AGB can be estimated from tree height and DBH. The stock volume is the total volume of living matter of a tree, expressed in cubic meters, and it can be estimated from the DBH. The estimation equations for forest parameters are closely related to tree species, and the models for different tree species are not universal. After tree species identification and updates the tree species information can be considered accurate and applied to the estimation of forestry parameters.

The estimation factor of the DBH estimation model is the tree height (H, in m), and the output is the DBH (D, in cm). The specific equations and application conditions are shown in Table 4.

The factors of the AGB models are tree height (H, in m) and DBH (D, in cm), and the output is AGB (M, in t/ha). The specific equations are shown in Table 5.

The same type of equation is used for the calculation of the stock volume of different tree species.

$$V=a·D^b·{(d+\frac{e}{D+k})}^c$$

(15)

where V: stock volume, in m³

D: DBH, in cm

The rest of the parameters are shown in the Table 6.

3. Results

3.2. Tree Segmentation

By setting the grid size d, the initial position of trees can be obtained. Using rotating profiles with an angular step of θ, segmented canopies can be generated while optimizing the tree positions. The initial positions are the basis for subsequent treetop points acquisition and final segmentation. Therefore, in order to obtain more accurate initial positions, it is necessary to set a more appropriate grid size. The setting of d is based on the size of the tree canopy in the sample plots, which is concentrated in the range of 2-4 meters in diameter. The value of d should be close to the radius of the tree canopy to ensure that the positions of the trees can be identified to the greatest extent possible. In the experiment, d is set to 1m, 1.5m and 2m to segment the eight sample plots.

Figure 16. The number of trees segmented from each sample plot.

Figure 16. The number of trees segmented from each sample plot.

Due to the complexity of natural forest conditions, overlapping canopies can cause short trees to be mistaken as parts of tall trees, resulting in multiple trees being recognized as one tree. Additionally, some trees are divided into multiple trees because of their large canopy or complex shape. These discrepancies cause under-segmentation and over-segmentation in the segmentation results. To evaluate the correct segmentation, under-segmentation and over-segmentation of each plot under different grid sizes d, recall, precision and F-score [48] are used as indicators to verify and evaluate the accuracy of individual tree segmentation (Equations 16, 17 & 18). The results are shown in Table 8.

$$r=\frac{TP}{TP+FN}$$

(16)

$$p=\frac{TP}{TP+FP}$$

(17)

$$F=\frac{2\left(r\times p\right)}{r+p}$$

(18)

where r, p& F: recall, precision and F-score

TP: number of correctly segmented trees

FN: number of under-segmented trees

FP: number of over-segmented trees

Table 8. The TP, FN and FP of different grid sizes d.

Table 8. The TP, FN and FP of different grid sizes d.

Plot ID	Real Number	d= 1 m			d= 1.5 m			d= 2 m
		TP	FN	FP	TP	FN	FP	TP	FN	FP
1	50	49	1	14	47	3	2	42	8	1
2	42	41	1	9	40	2	0	34	8	0
3	45	45	0	11	44	1	0	39	6	0
4	41	41	0	7	40	1	1	37	4	1
5	29	29	0	6	29	0	2	25	4	1
6	41	41	0	8	40	1	0	36	5	0
7	31	31	0	5	31	0	1	29	2	0
8	39	39	0	5	37	2	0	35	4	0
All	318	316	2	65	308	10	6	277	41	3

Among the three selected grid sizes, when d=1.5m, the total number of under-segmented and over-segmented trees is the lowest. The number of correctly segmented trees is similar to the actual number of trees in each plot, with a total of 308 correctly segmented trees, which was close to the total actual number of 318 trees. The recall, precision and F-score are calculated in Table 9, and the F-scores for all eight plots exceed 0.95 when d=1.5m, with the highest overall F-score of 0.976 in three grid sizes, indicating that the rotational profile algorithm achieves good segmentation results. Therefore, the part of the correctly segmented trees with d=1.5m was selected for subsequent species identification and parameters extraction.

Figure 17. Segmented point cloud of each plot.

Figure 17. Segmented point cloud of each plot.

3.3. Plots Matching and Tree Species Identification

The main purpose of plot matching is to solve the problem of positional differences between measured data and point cloud. To this end, this paper proposes the concept of matching degree and uses dominant trees for plot matching. However, due to the complexity of the actual situation, it is difficult to achieve complete consistency between the two datasets. Therefore, the deviation of the position of each tree is statistically collected and the deviation values are calculated, as shown in Table 10. Plot matching is considered complete when the matching degrees of the eight plots are greater than 90% and the deviation values are less than the threshold of 8m.

Based on the plots matching, accurate tree positions can be obtained. In forest surveys, the trees species information in sample plots is mostly accurate. This paper utilizes species information from sample plots to identify and extract individual erroneous tree species information, thereby providing correct species information for obtaining accurate forestry parameters in subsequent forestry surveys.

This article successfully converts the classification problem of 3D point clouds into 1D feature vector classification. It also effectively preserves the detailed characteristic information, fully utilizing spatial structural information to summarize the geometric morphology of each tree, thus achieving tree species identification, classification and updating. The main tree species in the plots include pine, oak, birch, elm, poplar, maple and other species. Several typical profiles are selected as training samples to obtain the morphological structure and parameters of each species. Table 10 shows that when the sample size of a certain tree species is insufficient, more data can be obtained by selecting profiles in multiple directions.

The trained model is applied to achieve classification. To evaluate the classification performance, the correct rate is defined as the ratio of the number of successfully classified trees to the number of misclassified trees. The method proposed is used to update the tree species information of the trees in eight plots, and the results are shown in Table 11. Including "Others", which is a collection of trees with extremely low numbers, the overall correct rate reaches 90.9%, providing a reliable tree species classification method and model for forest surveys.

The details are listed in the form of a confusion matrix in Table 13, which can more accurately indicate the error correction information for each tree species.

Table 12 and Table 13 show that most tree species can be corrected with 100% accuracy, and the few misclassified tree species are due to different degrees of damage to the tree canopy or inaccurate profile structural information caused by being close to adjacent trees during segmentation. Overall, using LiDAR data to correct tree species information is a reliable method with high accuracy, suitable for classifying most normally growing tree species.

Table 12. Accuracy of the tree identification.

Table 12. Accuracy of the tree identification.

Tree Species	Correct Classification	Type Ⅰ Error 1	Type Ⅱ Error 2	Corrected Tree Number	Correct Rate	Overall Correct Rate
Pine	39	1	2	3	100%	90.9%
Oak	25	0	2	2	100%
Birch	14	0	1	1	100%
Elm	17	2	1	3	100%
Linden	107	7	0	6	85.7%
Poplar	18	1	1	2	100%
Maple	45	0	3	2	66.7%
Others	32	0	1	1	100%

¹ Type I Error: record a tree belonging to the correct tree species as another species; ² Type II Error: record a tree belonging to another species as the current tree species.

Table 13. Confusion matrix of tree species misinformation correction.

Table 13. Confusion matrix of tree species misinformation correction.

	Collected	Pine	Oak	Birch	Elm	Linden	Polar	Maple	Others	CR 1
True
Pine		39	0	0	1	0	0	0	0	100%
Oak		0	25	0	0	0	0	0	0	100%
Birch		0	0	14	0	0	0	0	0	100%
Elm		1	0	0	17	0	0	1	0	100%
Linden		1	2	1	0	107	1	1	1	87.5%
Poplar		0	0	0	0	0	18	1	0	100%
Maple		0	0	0	0	0	0	45	0	66.7%
Others		0	0	0	0	0	0	0	32	100%

¹ CR: Correct Rate.

3.4. Forestry Parameters Extraction

After tree species identification and updating, information is taken from the major tree species’ segmented point cloud, and is utilized to estimate forestry parameters. Linear regressions are performed between the measured forestry parameters and the estimated parameters.

Figure 18. Linear regression plots of six parameters. Red lines are linear regression lines. (a) Tree height; (b) East-west canopy width; (c) North-south canopy width; (d) DBH; (e) AGB; (f) Stock volume.

Figure 18. Linear regression plots of six parameters. Red lines are linear regression lines. (a) Tree height; (b) East-west canopy width; (c) North-south canopy width; (d) DBH; (e) AGB; (f) Stock volume.

To quantitatively describe the estimation results of each parameter, the coefficient of determination (R²), the mean absolute error (MAE) and the root mean square error (RMSE) are chosen to evaluate the estimation accuracy (Table 14).

$$R^2=1-\frac{{\sum }_{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(18)

$$MAE=\frac{1}{n}\sum _{i=1}^n\left|y_i-\widehat{y_i}\right|$$

(19)

$$RMSE=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2}$$

(20)

where $y_i$: the measured value. (1 ≤ i ≤ n, n is the number of measured values)

$\widehat{y_i}$: the estimated value

$\overline{y}$: average of measured values

Table 14. Estimation accuracy of forestry parameters.

Table 14. Estimation accuracy of forestry parameters.

Parameters	Tree Height	Canopy Width		DBH	AGB	Stock Volume
		East-West	North-South
R2	0.893	0.757	0.694	0.840	0.896	0.891
RMSE	1.793	0.719	0.638	3.735	70.914	0.090
MAE	1.465	0.575	0.511	3.015	48.875	0.066

The estimation results of four forestry parameters, namely, tree height, DBH, AGB, and stock volume, show strong consistency with the measured values, with R² exceeding 0.8. Tree height shows higher consistency in the range of 20-30 m, and in the range of 5-20 m, the estimated tree heights are slightly higher than the measured heights in general, which was similarly reflected in the estimation of DBH. The accuracy of canopy width estimation is relatively lower (R²=0.757 for east-west and R²=0.694 for north-south), and the estimated values are smaller than the measured values in general. This also confirms that the use of LiDAR data is expected to become a powerful tool for future forestry surveys, significantly reducing human involvement and achieving rapid and accurate acquisition of large-scale forest parameters.

4. Discussion

5. Conclusions

References

1. Eckert, S. Improved Forest Biomass and Carbon Estimations Using Texture Measures from WorldView-2 Satellite Data. Remote Sens. 2012, 4, 810-829. [CrossRef]
2. He, X.; Ren, C.; Chen, L.; Wang, Z.; Zheng, H. The Progress of Forest Ecosystems Monitoring with Remote Sensing Techniques. Scientia Geographica Sinica 2018, 38, 97-1011. [CrossRef]
3. MALHI, Y.; BALDOCCHI, D. D.; JARVIS, P. G. The carbon balance of tropical, temperate and boreal forests. Plant, Cell and Environment 1999, 22, 715-740. [CrossRef]
4. Dong, X.; Li J.; Chen, H.; Zhao L.; Zhang, L.; Xing S. Extraction of individual tree information based on remote sensing images from an Unmanned Aerial Vehicle. Journal of Remote sensing 2019, 23, 1269-1280. [CrossRef]
5. Liu, F.; Tan, C.; Zhang, G.; Liu, J. Estimation of Forest Parameter and Biomass for Individual Pine Trees Using Airborne LiDAR. Transactions of the Chinese Society for Agricultural Machinery 2013, 44, 219-224,242. [CrossRef]
6. Ouma, Y. O. Optimization of Second-Order Grey-Level Texture in High-Resolution Imagery for Statistical Estimation of Above-Ground Biomass. Journal of Environmental Informatics 2006, 8, 70-85. [CrossRef]
7. Marshall, M.; Thenkabail, P. Advantage of hyperspectral EO-1 Hyperion over multispectral IKONOS, GeoEye-1, WorldView-2, Landsat ETM+, and MODIS vegetation indices in crop biomass estimation. ISPRS Journal of Photogrammetry and Remote Sensing 2015, 108, 205-218. [CrossRef]
8. Mohammadi, J.; Shataee, S.; Babanezhad, M. Estimation of forest stand volume, tree density and biodiversity using Landsat ETM+Data, comparison of linear and regression tree analyses. Procedia Environmental Sciences 2011, 7, 299-304. [CrossRef]
9. Franklin, S. E.; Hall, R. J.; Smith, L.; Gerylo, G. R. Discrimination of conifer height, age and crown closure classes using Landsat-5 TM imagery in the Canadian Northwest Territories. International Journal of Remote Sensing 2003, 24, 1823-1834. [CrossRef]
10. Brown, S.; Pearson, T.; Slaymaker, D.; Ambagis, S.; Moore, N.; Novelo, D.; Sabido, W. CREATING A VIRTUAL TROPICAL FOREST FROM THREE-DIMENSIONAL AERIAL IMAGERY TO ESTIMATE CARBON STOCKS. Ecological Application 2005, 15, 1083–1095. [CrossRef]
11. Cloude, S. R.; Papathanassiou, K. P. Forest Vertical Structure Estimation using Coherence Tomography. IGARSS 2008 - 2008 IEEE International Geoscience and Remote Sensing Symposium, Boston, USA, 2008. [CrossRef]
12. Blomberg, E.; Ferro-Famil, L.; Soja, M. J.; Ulander, L. M. H.; Tebaldini, S. Forest Biomass Retrieval From L-Band SAR Using Tomographic Ground Backscatter Removal. IEEE Geoscience and Remote Sensing Letters 2018, 15, 1030–1034. [CrossRef]
13. Matasci, G.; Hermosilla, T.; Wulder, M. A.; White, J. C.; Coops, N. C.; Hobart, G. W.; Zald, H. S. J. Large-area mapping of Canadian boreal forest cover, height, biomass and other structural attributes using Landsat composites and lidar plots. Remote Sensing of Environment 2018, 209, 90–106. [CrossRef]
14. Luckman, A.; Baker, J.; Honzák, M.; Lucas, R. Tropical Forest Biomass Density Estimation Using JERS-1 SAR: Seasonal Variation, Confidence Limits, and Application to Image Mosaics. Remote Sensing of Environment 1998, 63, 126–139. [CrossRef]
15. Means, J. Use of Large-Footprint Scanning Airborne Lidar To Estimate Forest Stand Characteristics in the Western Cascades of Oregon. Remote Sensing of Environment 1999, 67, 298–308. [CrossRef]
16. Hyyppä, J.; Yu, X.; Hyyppä, H.; Vastaranta, M.; Holopainen, M.; Kukko, A.; Kaartinen, H.; Jaakkola, A.; Vaaja, M.; Koskinen, J.; Alho, P. Advances in Forest Inventory Using Airborne Laser Scanning. Remote Sens. 2012, 4, 1190-1207. [CrossRef]
17. Kaartinen, H.; Hyyppä, J.; Yu, X.; Vastaranta, M.; Hyyppä, H.; Kukko, A.; Holopainen, M.; Heipke, C.; Hirschmugl, M.; Morsdorf, F.; Næsset, E.; Pitkänen, J.; Popescu, S.; Solberg, S.; Wolf, B.M.; Wu, J.-C. An International Comparison of Individual Tree Detection and Extraction Using Airborne Laser Scanning. Remote Sens. 2012, 4, 950-974. [CrossRef]
18. Kampa, K.; Slatton, K.C. An adaptive multiscale filter for segmenting vegetation in ALSM data. IGARSS 2004. 2004 IEEE International Geoscience and Remote Sensing Symposium, Anchorage, USA, 2004. [CrossRef]
19. Lee, H.; Slatton, K. C.; Roth, B. E.; Cropper, W. P. Adaptive clustering of airborne LiDAR data to segment individual tree crowns in managed pine forests. International Journal of Remote Sensing 2010, 31, 117–139. [CrossRef]
20. Mei, C.; Durrieu, S. Tree crown delineation from digital elevation models and high resolution imagery. proceedings of the isprs working group part. 2004.
21. Chen, Q.; Baldocchi, D.; Gong, P.; Kelly, M. Isolating Individual Trees in a Savanna Woodland Using Small Footprint Lidar Data. Photogrammetric Engineering and Remote Sensing 2006, 72, 923–932. [CrossRef]
22. Tang, F.; Zhang, X.; Liu, J. Segmentation of tree crown model with complex structure from airborne LIDAR data. Geoinformatics 2007: Remotely Sensed Data and Information. 2007. [CrossRef]
23. Leckie, D.; Gougeon, F.; Hill, D.; Quinn, R.; Armstrong, L.; Shreenan, R. Combined high-density lidar and multispectral imagery for individual tree crown analysis. Canadian Journal of Remote Sensing 2003, 29, 633–649. [CrossRef]
24. Reitberger, J.; Schnörr, Cl.; Krzystek, P.; Stilla, U. 3D segmentation of single trees exploiting full waveform LIDAR data. ISPRS Journal of Photogrammetry and Remote Sensing 2009, 64, 561–574. [CrossRef]
25. Lu, X.; Guo, Q.; Li, W.; Flanagan, J. A bottom-up approach to segment individual deciduous trees using leaf-off lidar point cloud data. ISPRS Journal of Photogrammetry and Remote Sensing 2014, 94, 1–12. [CrossRef]
26. Wang, Y.; Weinacker, H.; Koch, B. A Lidar Point Cloud Based Procedure for Vertical Canopy Structure Analysis And 3D Single Tree Modelling in Forest. Sensors 2008, 8, 3938–3951. [CrossRef]
27. Morsdorf, F.; Meier, E.H.; Allgöwer, B.; Nüesch, D. CLUSTERING IN AIRBORNE LASER SCANNING RAW DATA FOR SEGMENTATION OF SINGLE TREES. 2003.
28. Li, W.; Guo, Q.; Jakubowski, M. K.; Kelly, M. A New Method for Segmenting Individual Trees from the Lidar Point Cloud. Photogrammetric Engineering & Remote Sensing 2012, 78, 75–84. [CrossRef]
29. Holmgren, J.; Persson, Å. Identifying species of individual trees using airborne laser scanner. Remote Sensing of Environment 2004, 90, 415–423. [CrossRef]
30. Othmani, A.; Lew Yan Voon, L. F. C.; Stolz, C.; Piboule, A. Single tree species classification from Terrestrial Laser Scanning data for forest inventory. Pattern Recognition Letters 2013, 34, 2144–2150. [CrossRef]
31. Lin, Y.; Hyyppä, J. A comprehensive but efficient framework of proposing and validating feature parameters from airborne LiDAR data for tree species classification. International Journal of Applied Earth Observation and Geoinformation 2016, 46, 45–55. [CrossRef]
32. Kim, S.; McGaughey, R. J.; Andersen, H.-E.; Schreuder, G. Tree species differentiation using intensity data derived from leaf-on and leaf-off airborne laser scanner data. Remote Sensing of Environment 2009, 113, 1575–1586. [CrossRef]
33. Ørka, H. O.; Næsset, E.; Bollandsås, O. M. Classifying species of individual trees by intensity and structure features derived from airborne laser scanner data. Remote Sensing of Environment 2009, 113, 1163–1174. [CrossRef]
34. Kashani, A.; Olsen, M.; Parrish, C.; Wilson, N. A Review of LIDAR Radiometric Processing: From Ad Hoc Intensity Correction to Rigorous Radiometric Calibration. Sensors 2015, 15, 28099–28128. [CrossRef]
35. Solodukhin V. I.; Zukov A. J.; Mazugin I. N. Laser aerial profiling of a forest. Lew NIILKh Leningrad Lesnoe Khozyaistvo 1977, 10, 53-58.
36. Shrestha, R.; Wynne, R.H. Estimating Biophysical Parameters of Individual Trees in an Urban Environment Using Small Footprint Discrete-Return Imaging Lidar. Remote Sens. 2012, 4, 484-508. [CrossRef]
37. Bortolot, Z. J.; Wynne, R. H. Estimating forest biomass using small footprint LiDAR data: An individual tree-based approach that incorporates training data. ISPRS Journal of Photogrammetry and Remote Sensing 2005, 59, 342–360. [CrossRef]
38. Wang, D.; Xin, X.; Shao, Q.; Brolly, M.; Zhu, Z.; Chen, J. Modeling Aboveground Biomass in Hulunber Grassland Ecosystem by Using Unmanned Aerial Vehicle Discrete Lidar. Sensors 2017, 17, 180. [CrossRef]
39. García, M.; Riaño, D.; Chuvieco, E.; Danson, F. M. Estimating biomass carbon stocks for a Mediterranean forest in central Spain using LiDAR height and intensity data. Remote Sensing of Environment 2010, 114, 816–830. [CrossRef]
40. DONOGHUE, D.; WATT, P.; COX, N.; WILSON, J. Remote sensing of species mixtures in conifer plantations using LiDAR height and intensity data. Remote Sensing of Environment 2007, 110, 509–522. [CrossRef]
41. Jin, J.; Yue, C.; Li, C.; Gu, L.; Luo, H.; Zhu, B. Estimation on Forest Volume Based on ALS Data and Dummy Variable Technology. FOREST RESOURCES WANAGEMENT 2021, 0, 77-85. [CrossRef]
42. PANG, Y.; LI, Z.-Y. Inversion of biomass components of the temperate forest using airborne Lidar technology in Xiaoxing’an Mountains, Northeastern of China. Chinese Journal of Plant Ecology 2013, 36, 1095–1105. [CrossRef]
43. Axelsson, P. Processing of laser scanner data—algorithms and applications. ISPRS Journal of Photogrammetry and Remote Sensing 1999, 54, 138–147. [CrossRef]
44. Qian, C.; Yao, C.; Ma, H.; Xu, J.; Wang, J. Tree Species Classification Using Airborne LiDAR Data Based on Individual Tree Segmentation and Shape Fitting. Remote Sens. 2023, 15, 406. [CrossRef]
45. Goodenough, D. G.; Chen, H.; Dyk, A.; Richardson, A.; Hobart, G. Data Fusion Study Between Polarimetric SAR, Hyperspectral and Lidar Data for Forest Information. IGARSS 2008 - 2008 IEEE International Geoscience and Remote Sensing Symposium. Boston, USA, 2008, . [CrossRef]
46. Alonzo, M.; Bookhagen, B.; Roberts, D. A. Urban tree species mapping using hyperspectral and lidar data fusion. Remote Sensing of Environment 2014, 148, 70–83. [CrossRef]
47. Hinton, G. E.; Osindero, S.; Teh, Y.-W. A Fast Learning Algorithm for Deep Belief Nets. Neural Computation 2006, 18, 1527–1554. [CrossRef]
48. Goutte, C.; Gaussier, E. A Probabilistic Interpretation of Precision, Recall and F-Score, with Implication for Evaluation. Advances in Information Retrieval 2005, 345–359. [CrossRef]