Optimizing the Spatial Structure of Metasequoia Plantation Forest Based on UAV-LiDAR and Backpack-LiDAR

Abstract

Optimizing the spatial structure of forests is important for improving the quality of forest ecosystems. Light detection and ranging (LiDAR) could accurately extract forest spatial structural parameters, which has significant advantages in spatial optimization and resource monitoring. In this study, we used unmanned aerial vehicle LiDAR (UAV-LiDAR) and backpack-LiDAR to acquire point cloud data of Metasequoia plantation forests from different perspectives. Then the parameters, such as diameter at breast height and tree height, were extracted based on the point cloud data, while the accuracy was verified using ground-truth data. Finally, a single-tree-level thinning tool was developed to optimize the spatial structure of the stand based on multi-objective planning and the Monte Carlo algorithm. The results of the study showed that the accuracy of LiDAR-based extraction was (R² = 0.96, RMSE = 3.09 cm) for diameter at breast height, and the accuracy of R² and RMSE for tree height extraction were 0.85 and 0.92 m, respectively. Thinning improved stand objective function value Q by 25.40%, with the most significant improvement in competition index CI and openness K of 17.65% and 22.22%, respectively, compared to the pre-optimization period. The direct effects of each spatial structure parameter on the objective function values were ranked as follows: openness K (1.18) > aggregation index R (0.67) > competition index CI (0.42) > diameter at breast height size ratio U (0.06). Additionally, the indirect effects were ranked as follows: aggregation index R (0.86) > diameter at breast height size ratio U (0.48) > competition index CI (0.33). The study realized the optimization of stand spatial structure based on double LiDAR data, providing a new reference for forest management and structure optimization.

1. Introduction

Forest structure includes both spatial and non-spatial structures, both of which are attributes that must be measured and regulated in forest management planning [1,2]. The non-spatial structure of a forest is often used to describe the average state of forest attributes such as average stand age, average tree height, average densities, and so on [1]. Forest spatial structure refers to the spatial structure of species within the same forest community, i.e., the horizontal distribution pattern of trees and the spatial arrangement of their attributes [3], including the aggregation index [4], the diameter at breast height size ratio [5], the competition index [6], and the openness [7]. The main difference with non-spatial structure is that forest spatial structure is dependent on the spatial location of the stand, and it also determines to a large extent the stability of the forest, the possibilities for development, and the size of the management space [8], making forest spatial structure more important than non-spatial structure [9].

Improving the quality of forest ecosystems and thus their carbon sequestration capacity by optimizing the stand spatial structure is a feasible approach [10]. As an important measure to optimize stand regulation, thinning is a method to improve the growing environment and promote stand growth and regeneration by adjusting the number of standing trees in the stand [11,12,13]. However, how to develop a reasonable thinning plan to optimize stand spatial structure is a hot issue in forest management [14]. Traditional forest management schemes generally take the non-spatial structure adjustment of forest stands as their main means, i.e., to achieve the purpose of stand structure optimization by adjusting stand density, average age, and tree species composition [15]. With the proposal of structural forest management theory and new quantitative indicators for the spatial structure of forest stands, stand nurturing and thinning of forest stands for optimizing the spatial structure of forest stands have been gradually attracting attention. However, the spatial structure of forest stands is very complex, and each evaluation index of the spatial structure of forest stands not only involves different aspects but also is contradictory and interdependent with each other [16], so the optimization of the spatial structure of forest stands cannot be achieved by only setting a single objective. Based on the theory of multi-objective planning, by setting the conflicting evaluation indexes as different objective functions with corresponding constraints, the combined optimization of all objective functions could be achieved by solving the function. Subsequently, Dong et al. [17], Qiu et al. [18], Sun et al. [19], and Bettinger et al. [20] started based on different objectives and used heuristic algorithms to solve the problem, and all of them finally succeeded in solving the forestry planning problem, which also provides an effective method for forest management control based on the optimization of forest stand spatial structure.

The acquisition of forest parameters with spatial location information is a prerequisite for measuring and evaluating the spatial structure of forests. However, related studies around the world have shown that most forest parameters are obtained through manual field surveys, which have disadvantages such as low efficiency, high cost, and small scale [21,22]. It is therefore important for forest resource monitoring and management to minimize the time and labor costs of field surveys while achieving accurate estimates of forest parameters and their structures [23]. In recent years, the use of Light Detection and Ranging (LiDAR) technology has provided a new approach for accurate estimation of forest parameters and has also been widely applied [24,25]. LiDAR can generate accurate 3D point cloud data of a target feature by calculating the time difference between the emission of a laser pulse and its return to the receiver. Based on the point cloud data, we can reconstruct indoors the 3D structural features of individual trees to extract forest parameters and other information, or we can realize parameter estimation by extracting the characteristic variables of the point cloud data and building a regression model with the ground measurement data [26,27]. Depending on the installation platform and mode of operation, LiDAR can be further divided into two categories: airborne LiDAR and ground-based LiDAR [28]. The biggest advantage of unmanned aerial vehicle LiDAR (UAV-LiDAR), as a representative of airborne LiDAR, was that it can acquire high-density point cloud data in the study area at low cost, which was very useful for the extraction of single tree parameters. However, the point cloud density obtained by UAV-LiDAR usually decreases as the height of the stand decreases, i.e., the point cloud density is higher at the top of the canopy, but the number of point cloud distributions below the canopy cover is lower [29], which results in a high error in the estimation of parameters covered by the canopy. Unlike UAV-LiDAR, backpack-LiDAR uses a bottom-up data acquisition method, which is well suited for the estimation of parameters such as diameter at breast height and branch height, as it can record the under-canopy characteristics of trees and present them as a high-density point cloud [30]. However, due to the mutual occlusion of the canopy, it is difficult for the laser pulses to effectively penetrate the canopy and thus capture information from the treetops, which leads to a biased estimation of tree height by the backpack LiDAR [31]. In view of the shortcomings of single-platform LiDAR in parameter estimation, many scholars have tried to combine multi-platform LiDAR systems to overcome the above limitations. For example, Lu et al. [32] combined UAV and backpack LiDAR to estimate the above ground biomass of forests, and the results showed that the combination of dual radar can improve the estimation accuracy. Bazezew et al. [33] and Giannetti et al. [34] integrated data obtained from airborne laser scanners (ALS) and terrestrial laser scanners (TLS) and estimated forest structure parameters such as above ground biomass, crown diameters, and tree height, respectively. The results showed that the combined use of ALS and TLS data can effectively improve the accuracy of parameter estimation. Fekry et al. [29] fused ground-based and UAV-LiDAR data and estimated forest parameters, resulting in significant improvements in the estimation accuracy of tree height, tree volume, and crown volume over single LiDAR. These demonstrate that the combination of multiple platforms of the LiDAR system can effectively compensate for each other’s shortcomings, achieve complementary advantages, and improve the estimation accuracy of forest structure parameters.

As a fast-growing tree species, Metasequoia is widely planted in more than 50 countries and regions around the world and has a broad research base. At the same time, Metasequoia is an important species for the creation of planted carbon sink forests in China, and optimizing its stand spatial structure is also important for the management of carbon sink forests. Therefore, Metasequoia plantations were used as the target of this study. Firstly, UAV-LiDAR and backpack-LiDAR were used together to obtain accurate information on the diameter at breast height, tree height, and spatial location of the standing trees. Then, four spatial structure parameters were calculated: aggregation index, diameter at breast height size ratio, competition index, and openness. Finally, a dynamic multi-objective optimization model of the spatial structure of the forest was constructed and solved using the Monte Carlo algorithm to optimize the spatial structure of the forest stand. At the same time, a structural equation model was used to quantitatively analyze the contribution of each spatial structure parameter to the optimization of the spatial structure of the forest. The results of the study provided new ideas and methods for intelligent optimization and quality improvement of forest spatial structure based on UAV-LiDAR and backpack-LiDAR.

2. Materials and Methods

This study aimed to use UAV-LiDAR and backpack-LiDAR data to accurately estimate forest parameters and optimize the spatial structure of forest stands on this basis. As shown in Figure 1, it can be seen that the realization of this study was divided into a total of four phases. Additionally, the detailed steps of these four phases are explained in Section 2.3, Section 2.4, Section 2.5 and Section 2.6.

2.1. Study Area

As shown in Figure 2, the study area was in Qingshan Lake National Forest Park, Lin’an District, Hangzhou, Zhejiang Province (119°45′40″E, 30°15′26″N~119°45′46″E, 30°15′31″N). The area is a typical Metasequoia pure plantation, introduced from South America in the 1960s, with a stand density of about 1200 plants·hm⁻². The topography of the area is dominated by low-elevation hills, with a slight rise from southwest to northeast, averaging 30 m above sea level. The climate belongs to the subtropical monsoon climate, with warm and humid conditions and four distinct seasons. The average annual temperature is 15.8 °C, the extreme minimum temperature is −12 °C, and the extreme maximum temperature is 41 °C. The annual precipitation is approximately 1628.6 mm, mainly in the summer.

2.2. Data Acquisition

2.2.1. UAV-LiDAR Data

The UAV-LiDAR point cloud data was acquired on 13 December 2022, using a LiDAR system on board a DJI Matrice 600 Pro six-rotor UAV. The system integrates a Velodyne Puck LITE laser scanner (the sensor parameters are shown in

Table 1. LiDAR sensor parameters for multiple systems.

Parameters	UAV-LiDAR	Backpack-LiDAR
Sensors	Velodyne Puck LITE×1	Velodyne Puck VLP-16×2
Max. range	100 m	100 m
Ranging accuracy	±3 cm	±3 cm
Wavelength	903 nm	903 nm
Vertical FOV	± 15°	± 15°
Horizontal FOV	360°	360°
Scanning speed	300,000 pts/s	600,000 pts/s

), a GNSS receiver, and an IMU inertial navigation system, with the relevant equipment and components shown in Figure 3. The data was collected at an altitude of 60 m, a flight speed of 8 m/s, a course spacing of 25 m, a 50% collateral overlap of data sampling, and a final average point cloud density of 200 points/m². An example of the point cloud data is shown in Figure 3c.

2.2.2. Backpack-LiDAR Data

The backpack-LiDAR data was collected using the LiBackpack DGC50 LiDAR system on 14 December 2022. The system is equipped with a microcomputer, GNSS antenna, panoramic camera, hand-held touch tablet, and Velodyne Puck VLP-16 LiDAR sensors in both horizontal and vertical directions (the sensor parameters are shown in Table 1), and the associated equipment and components are shown in Figure 3. During the data collection process, the surveyor walked slowly through the forest area along a zigzag route with the backpack LiDAR equipment on his shoulder to ensure stability. An example of the point cloud data is shown in Figure 3d.

2.2.3. Field Inventory Data

In parallel with the UAV-LiDAR and backpack-LiDAR data acquisition, 113 Metasequoia plants in the study area were measured and located in this study. The diameter at 1.3 m was measured with a diameter-checking steel ruler, and the height of the tree was measured using a Blume–Leiss height gauge, while the position of the single tree was also accurately located using a Huace Smart RTK (better than 0.5 m accuracy) [22,35]. The statistical characteristics of the field-measured data are shown in

Table 2. Table of statistical characteristics of field-measured data.

Attributes	Maximum	Minimum	Range	Average	Standard Deviation
DBH/cm	58	14.3	43.7	29.2	8.9
H/m	29.5	18.2	11.3	24.9	2.2

.

2.3. LiDAR Data Pre-Processing

In this study, the same pre-processing operations were performed on the data collected by UAV-LiDAR and backpack-LiDAR. The raw point cloud data were first processed by denoising and filtering [36]. The ground and non-ground points were then classified using the improved progressive TIN densification (IPTD) algorithm proposed by Zhao et al. [37] to obtain ground and non-ground points, and the resulting ground points were further processed by Kriging interpolation to obtain a digital elevation model (DEM) with a spatial resolution of 0.5 m. Finally, the DEM was used to normalize the point cloud data in order to reduce the influence of the topography on the tree height estimation [38].

2.4. UAV-LiDAR and Backpack-LiDAR Matching

Although the backpack-LiDAR device is equipped with its own GNSS antenna, the satellite signal is affected by the tree canopy as it is under the tree canopy for most of the data acquisition process, resulting in inaccurate positioning obtained by backpack-LiDAR [39]. Therefore, the data acquired by backpack-LiDAR were considered relative coordinates, while the data acquired by UAV-LiDAR, which flies above the tree canopy with minimal interference from the outside, were considered absolute coordinates in this study. Then, the same name feature point matching method was used [40,41], i.e., finding the same points in the 3D point cloud for coordinate transformation, thus enabling backpack-LiDAR and UAV-LiDAR to be in the same coordinate system.

2.5. Forest Parameter Extraction and Accuracy Verification

2.5.1. The Diameter at Breast Height Extraction

In this study, we first used the density-based spatial clustering of applications with noise (DBSCAN) algorithm on normalized backpack-LiDAR data to achieve accurate segmentation of tree trunks. The DBSCAN algorithm could realize classless clustering of point cloud data by using the two parameters of minimum cluster points (MinPts) and neighborhood radius (Eps). Based on the characteristics of standing trees, this study set MinPts to 100 and Eps to 0.2 m. Meanwhile, to avoid the results caused by low-ground vegetation, only point cloud data above 1 m were used for single tree cluster segmentation [42]. After that, this study extracted the vertical point cloud slices of each standing tree at a specified height (labeled as SH) at 1.3 m from the ground, and finally, the least squares combined with the cylindrical fitting method was used to obtain the diameter at breast height and the spatial location of the standing trees (as shown in Figure 4a) [43].

2.5.2. Tree Height Extraction

In this study, the point cloud segmentation (PCS) algorithm [44], which was a single tree segmentation algorithm integrating the region-growing method and the thresholding method, was used for segmenting the single tree canopy on normalized UAV-LiDAR data. The algorithm first realized the segmentation of the first single wood by finding the canopy apex as the first seed point and then sequentially judging the distance between the surrounding points and this seed point and comparing it with the set threshold value [45]. Additionally, repeat the operation for the remaining point cloud data until all the standing trees were realized to be segmented. However, the traditional PCS algorithm makes it difficult to accurately find the canopy apex in the adhering and overlapping canopies, which leads to the poor accuracy of single tree segmentation. So this study referred to the study of Lu et al. [32], where the standing tree position information extracted from the backpack LiDAR was used as an assumed seed point to assist in the segmentation of a single tree (as shown in Figure 4b). Finally, based on the results of the single tree segmentation, we could get the tree heights of all the standing trees.

2.5.3. Accuracy Verification

In this study, the accuracy of the diameter at breast height and tree height extracted from LiDAR data was assessed by measuring the diameter at breast height and tree height on the ground, and the evaluation indicators included coefficient of determination (R²), root mean square error (RMSE), and relative root mean square error (rRMSE), with equations under (1)–(3):

$$R^2=1-\frac{\sum _{i=1}^n{(x_i-\widehat{x_i})}^2}{\sum _{i=1}^n{(x_i-{\overline{x}}_i)}^2},$$

(1)

$$RMSE=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{(x_i-{\widehat{x}}_i)}^2},$$

(2)

$$rRMSE=\frac{RMSE}{\overline{x}}\times 100\%,$$

(3)

where $x_i$ is the measured height or diameter at breast height of the ith standing tree; $\widehat{x_i}$ is the height or diameter at breast height value extracted from the LiDAR data of the ith standing tree; and $n$ is the number of trees in the standing tree.

2.6. Spatial Structure Indicators and Optimization Model Construction Techniques

2.6.1. Spatial Structural Unit of Forest Stands

The stand spatial structure unit is the basis for calculating the stand spatial structure parameters and evaluating the spatial characteristics of the stand. In this study, the widely used four-neighbor method was used to determine the spatial structural unit of the stand, i.e., a central tree was considered the object tree, and then the four standing trees closest to it horizontally were selected as competing trees to form a spatial structural unit [46,47]. At the same time, in order to avoid edge effects, the area 5 m from the boundary line of the sample plot was set as a buffer area, and trees within the buffer area only participated in the calculation as competing trees of the object trees, but not as object trees [48,49,50].

2.6.2. Indicators for Evaluating the Spatial Structure of Forest Stands

In this study, four spatial structure parameters were selected as evaluation indicators, including the aggregation index R to evaluate the distribution pattern of trees, the diameter at breast height size ratio U to reflect the degree of size differentiation of trees, the competition index CI to reflect the competition status of trees, and the openness K to reflect the growth space of standing trees. The formulas are defined as shown in (4–7):

(1)

Aggregation index (R)

$$R_i=\frac{r_i}{\frac{1}{2}\sqrt{\frac{F}{N}}},$$

(4)

where $r_i$ is the distance from the ith individual tree to its nearest neighbor (m), $F$ is the sample plot area (m²), and $N$ is the number of plants in the sample plot. The aggregation index is divided into three intervals [0, 1); (1, 2.1491); and 1, representing the aggregated distributions, uniform distributions, and random distributions, respectively;

(2)

The diameter to breast height size ratio (U)

$$U_i=\frac{1}{n}{\sum }_{j=1}^n{V}_{ij},$$

(5)

where $n$ is the number of adjacent tree plants of the target tree $i$, if the diameter at breast height of adjacent tree $j$ is larger than the diameter at breast height of reference tree $i$, $V_{ij}=1$, otherwise $V_{ij}=0$. The diameter at breast height size ratio U was judged on five levels: [0, 0.25); [0.25, 0.5); [0.5, 0.75); [0.75, 1); and 1, representing the dominance, sub-dominance, intermediate, sub-inferiority, and inferiority of the target tree within a spatial structural unit, respectively;

(3)

Hegyi Competition Index (CI)

$${CI}_i={\sum }_{j=1}^n\frac{d_j}{d_i·L_{ij}},$$

(6)

where $L_{ij}$ is the distance (m) between object tree $i$ and competitor tree $j$, $d_i$ is the diameter at breast height (cm) of object tree $i$, $d_j$ is the diameter at breast height (cm) of competitor tree $j$, and $n$ is the number of competing trees of object tree $i$. A higher competition index CI means more intense competition between stands, which is detrimental to the overall growth of the stand;

(4)

Openness (K)

$$K_i=\frac{1}{n}{\sum }_{j=1}^n\frac{Lij}{H_{ij}},$$

(7)

where $L_{ij}$ is the horizontal distance (m) between object tree $i$ and nearest tree $j$, $H_{ij}$ is the tree height (m) of nearest tree $j$ from $i$, and $n$ is the number of competing trees in object tree $i$. The openness K is within five intervals of (0, 0.2]; (0.2, 0.3]; (0.3, 0.4]; (0.4, 0.5]; and (0.5, ∞], indicating severely inadequate, inadequate, relatively adequate, largely adequate, and very adequate internal growth space in the stand, respectively.

2.6.3. Spatial Structure Dynamic Multi-Objective Optimization Model Construction

In the case of planted forests, a more ideal stand spatial structure means less competition within the stand, more space for growth, and a more even horizontal distribution pattern [18,51]. Therefore, this study used the multiplication and division method to multiply the parameters that were expected to become larger and divide the parameters that were expected to become smaller, and finally, the individual objective functions were combined to form a composite function Q for evaluating and optimizing the stand spatial structure. The multi-objective optimization model for the stand spatial structure is shown in the equation in (8).

$$Q=f\left\{\max \left(R\right)·\mathrm{m}\mathrm{a}\mathrm{x} \left(K\right)·\mathrm{m}\mathrm{i}\mathrm{n} \left(CI\right)·\mathrm{m}\mathrm{i}\mathrm{n} \left(U\right)\right\}=max\left\{\frac{\left(1+R\right)\times (1+K)}{(1+CI)\times (1+U)}\right\},$$

(8)

subject to:

$$\overline{R}>R,$$

(9)

$$\overline{K}>K,$$

(10)

$$\overline{CI}<CI,$$

(11)

$$\overline{U}<U,$$

(12)

$$\overline{M_T}=M_O,$$

(13)

$$\overline{N_T}<N_{\mathrm{O}}\times P,$$

(14)

where Q is the multidimensional stand spatial structure index, i.e., the objective function value, which reflects the degree of superiority or inferiority of the stand spatial structure. The higher the value of the objective function, the more desirable the spatial structure of the stand. R, K, CI, and U are the aggregation index, openness, Heygi competition index, and diameter at breast height size ratio of the stand in the initial state of the study area, respectively; $\overline{R}$, $\overline{K}$, $\overline{CI}$, and $\overline{U}$ are the aggregation index, openness, Heygi competition index, and diameter at breast height size ratio of the stand after thinning, respectively; $M_O$ and $N_{\mathrm{O}}$ represent the number of diameter classes and the total number of trees within the plots in their initial state, respectively; $\overline{M_T}$ represents the number of diameter classes in the stand after thinning; $\overline{N_T}$ represents the number of thinning trees; and P indicates the intensity of stand harvest. Because of the possibility of a zero value for the CI, 1 was added to each parameter in this study to avoid invalidating the fractional equation.

In addition, to elucidate the extent to which each spatial structural parameter contributes to the variation of the objective function value, this study also used structural equation mode (SEM) to quantitatively analyze the direct and indirect effects of aggregation index R, competition index CI, the diameter at breast height size ratio U, and openness K on the variation of the objective function value Q [52].

2.6.4. Monte Carlo Algorithm Solving

The Monte Carlo algorithm, as a kind of heuristic algorithm, is simple and easy to implement [53], so in this study, the Monte Carlo algorithm was used to solve the dynamic multi-objective optimization model of the spatial structure (Equation (8)) [54]. The core idea of the Monte Carlo algorithm was to randomly select a tree as timber to be felled and then recalculate the characteristic parameter values and diameter classes of the spatial structure of the felled stand and compare them with the characteristic parameter values and diameter classes of the spatial structure of the stand in its initial state. If the set constraints were met, the proposed tree was output as a feasible solution for thinning while the values of the spatial structural characteristics parameters and diameter classes of the stand in the initial state were updated, and then the next tree was randomly selected. If the constraints are not met, the program jumps straight to selecting the next tree to be felled, and the cycle repeats until the harvesting intensity or simulated number of thinnings is reached (the initial harvest intensity in this study was set at 30%, and the number of simulated thinnings was 30,000) [51].

3. Results

3.1. Diameter at Breast Height Extraction Based on Backpack-LiDAR

In order to accurately extract the diameter at breast height of Metasequoia in the study area, six-point cloud slices of different heights were set up for diameter at breast height fitting in this study, and the results are shown in Figure 5. In general, point cloud slices based on different heights could accurately extract the diameter at breast height of Metasequoia from backpack-LiDAR data, and their R² was 0.96. However, a comprehensive analysis and comparison of the results in Figure 5a–f showed that the diameter at breast height RMSE extracted from point cloud slices at a height of 35 cm was the smallest (3.09 cm), and the rRMSE was only 10.45%. Therefore, 35 cm was chosen as the optimal point cloud slice height for the subsequent extraction of the diameter of all standing trees in the study area.

3.2. Tree Height Extraction Based on UAV-LiDAR

The correlation between tree height extracted using different single tree segmentation methods based on UAV-LiDAR data and field-measured data is shown in Figure 6. In general, both methods were able to extract tree heights more accurately for the identified single tree, with a high accuracy R² of 0.83 for direct single wood segmentation based on UAV-LiDAR data and a high accuracy R² of 0.85 for single tree segmentation with the help of auxiliary seed points. Secondly, the accuracy of single tree segmentation was compared, with the number of stems obtained from the single tree segmentation and the number of stems matching the actual number of stems of 108 with the help of auxiliary seeds (as shown in Figure 6b). In contrast, direct single tree segmentation yielded only 55 stem plants (as shown in Figure 6a). This suggested that the method used in Figure 6b, whereby backpack-LiDAR-extracted stand locations were used as seed points, could indeed improve the accuracy of single tree segmentation of UAV-LiDAR data in high-density stands. Accurate single tree segmentation also contributed to the height extraction, which explained the slightly higher accuracy of the trees extracted by this method than those extracted by direct single tree segmentation.

3.3. Results of Stand Structure Optimization

Figure 7a shows the results of the dynamic multi-objective optimization of the spatial structure of Metasequoia stands in the study area obtained by solving the Monte Carlo algorithm, and Figure 7b shows the change process of four spatial structure indexes, namely the aggregation index R, the diameter at breast height size ratio U, the competition index CI, and the openness index K. Analysis of Figure 7 showed that the objective function value Q increased with the number of simulated thinnings, indicating that the spatial structure of the stand was continually being optimized and strengthened. In terms of spatial structure indicators, an increase in the aggregation index R indicates a further tendency for the stand to be evenly distributed, while an increase in the openness index K indicates a better widening of the growing space of the standing trees within the stand. In contrast, a reduction in the competition index CI reflected a reduction in stand competition within the stand, and a simultaneous reduction in the diameter at breast height size ratio U indicated a gradual dominance of target trees within a stand spatial structure unit. When the number of simulations reached 15,300, the objective function and the spatial structure parameters no longer changed and remained stable, which meant that the spatial structure of the stand had reached a relative optimum.

The degree of direct and indirect influence of each spatial structure parameter on the value of the objective function is shown in Figure 8. The analysis of the graphs showed that the direct effects of the aggregation index R, the competition index CI, the diameter at breast height size ratio U, and the openness K on the values of the objective function are not entirely consistent. In order of direct influence (absolute value of the throughput coefficient): openness K (1.18) > aggregation index R (0.67) > competition index CI (0.42) > diameter at breast height size ratio U (0.06), it could be seen that the parameter with the greatest direct influence on the value of the objective function is openness. Next, we analyzed the indirect effect of each spatial structure parameter on the objective function value. The indirect influence of the aggregation index R on the value of the objective function through the other spatial structure parameters was 0.86; the indirect influence of the competition index CI on the value of the objective function through the other spatial structure parameters was 0.33; and the indirect influence of the diameter at breast height size ratio U on the value of the objective function through the other spatial structure parameters was 0.48. This also suggested that the aggregation index has the greatest indirect effect on the value of the objective function.

Changes in the spatial structure of the stand before and after thinning are shown in Figure 9a and

Table 3. Changes in stand spatial structure before and after thinning.

Parameters	Before Thinning	After Thinning	Change Tendency	Magnitude of Change/%
R	1.43	1.55	Increased	+8.40
U	0.50	0.46	Decreased	−8.00
CI	1.77	1.47	Decreased	−17.65
K	0.09	0.11	Increased	+22.22
Q	0.63	0.79	Increased	+25.40

, where the objective function value Q improved from 0.63 to 0.79 after thinning, an overall increase of 25.40%. In terms of spatial structure parameters, both aggregation index R and openness K had been significantly improved, with openness K increasing from 0.09 before thinning to 0.11, an increase of 22.22%, while aggregation index R was slightly less effective than openness K, increasing from 1.43 before thinning to 1.55, an increase of only 8.40%. On the contrary, the competition index CI and the diameter at breast height size ratio U showed a decreasing trend, with the competition index CI decreasing from 1.77 to 1.47, a decrease of 17.65%, while the diameter at breast height size ratio U also showed a decreasing trend, but it only decreased from 0.50 to 0.46, a change of 8.00%, which was not obvious. In addition, although all spatial structure parameters had been optimized to varying degrees compared to the pre-thinning period, the study area as a whole still exhibited a severe lack of internal growth space, strong overall competition, and an overall subdominant state.

In terms of non-spatial structure, the characteristics of the variation in the distribution of the total number of diameter classes and diameter at breast height intervals in the stand are shown in Figure 9b. Similar to the pre-thinning period, the average diameter at breast height of the stands remained around 26 cm. The distribution of diameter at breast height showed that the majority of trees below 26 cm were of small to medium diameter, which was conducive to the growth and regeneration of the stand, but the distribution of stand diameter did not show an inverted J-shape, which reflected that the optimized stand still lacked sufficient vitality and stability.

The relationship between the objective function value and the number of harvesting trees is shown in Figure 10. From the graph, it was clear that as the number of harvesting trees increased, the value of the objective function continued to increase until the end of the simulated thinning, when a total of 331 trees were cut by the model, giving a cutting intensity of 23.71%.

4. Discussion

Tree height extraction based on airborne LiDAR had been widely used and validated [55,56], but most studies were predicated on accurate single tree segmentation results. Although the PCS algorithm achieved good segmentation accuracy in low-density coniferous forests, it usually failed to achieve satisfactory segmentation accuracy in forests with higher stand densities. The use of known stem positions as seed points could, in theory, effectively improve the PCS algorithm’s ability to detect a single tree, but in practice, there were still certain biases in the process. As for simultaneous localization and mapping, the simultaneous (SLAM) algorithm used by the backpack LiDAR accumulated errors during the data stitching process [57], i.e., the longer the device was used, the greater the bias in the positioning of the stems. The lack of accurate stand locations affected the accuracy of single tree segmentation, which explained why the tree heights extracted in this study did not reach a more accurate height. Although this problem with SLAM algorithms could often be improved by reducing operation time and area, this could also reduce operational efficiency and operational flexibility.

An important point in the extraction of the diameter at breast height using the adaptive cylindrical fitting method is the selection of a suitable point cloud slice height. In this study, six different point cloud slice heights were used to fit the diameter at breast height. The results showed that the 35 cm point cloud slice height had the highest accuracy of the diameter at breast height, but there was no significant difference in accuracy with other point cloud slice heights. This was not entirely consistent with the results of Xie et al. [58], which might be related to the growth environment and species characteristics of the standing trees. The subject of this study was a plantation that has been in a state of conservation and has no obvious branches below 2 m, and the stem cleanliness is extremely high, which might not lead to significantly different heights of different point cloud slices when extracting diameter at breast height. Secondly, although the accuracy of the extraction of diameter at breast height was at a high level in this study, it was overall underestimated, a result that was consistent with the findings of Su et al. [59]. This might be due to two parts: firstly, the location of the diameter at breast height measured with a breast diameter scale was not exactly at 1.3 m, resulting in an overestimation of the measured diameter at breast height value; secondly, it was due to an algorithmic issue with diameter at breast height extraction, as LiDAR could detect and reflect very fine information, such as grooves on the bark surface, so the estimation of diameter at breast height was usually not based on the outermost point cloud of the bark but on the stem, which had the highest density of point clouds within the bark [60]. In addition, the results of the diameter at breast height extraction showed a certain amount of overestimation, which was found to be mainly caused by incomplete point clouds around the stem by comparison with the point cloud data. This was consistent with the conclusion of Heo et al. [61] that missing point clouds around tree stems led to large LiDAR-extracted diameters at breast height values and also confirmed the findings of Olofsson et al. [62] that the diameter at breast height could be more accurately depicted and measured by using a multi-station TLS.

Forest spatial structure optimization aims to produce healthy and stable forests through stand restructuring, so establishing tree-level harvesting plans with clear spatial location information is important for both forest structure restructuring and management [17]. The results of this study showed that all the spatial structure parameters could be improved by reasonable thinning, and among the four selected spatial structure parameters, the competition index CI and the openness K had the most obvious effect. Qiu et al. [18] showed that selective harvesting could reduce competition for limited resources and increase the growing space of individual trees by reducing the number of standing trees in a stand, which could explain the above results. Secondly, the results of this study also showed that the aggregation index R of the stand increased further after optimization, indicating that the stand tended to be more uniformly distributed, whereas Dong et al. [1] showed in their study that the spatial distribution pattern of the stand tended to be randomly distributed after optimization. The main reason for this difference is that the two studies were based on different optimization principles and had different objectives. The optimization of the spatial structure of the stands in this study was based on plantations with consistent forest age, which was fundamentally different from the optimization principles established by Dong et al. for natural forests based on inconsistently aged trees, and therefore optimization results were not the same for both. The direct influence of the spatial structure parameters on the objective function values showed that openness had the greatest direct influence on the objective function values, which indicated that the optimization of the spatial structure of the stand could be achieved directly and effectively by increasing the growing space of the standing trees in high-density stands. In terms of the degree of indirect influence of the spatial structure parameters on the value of the objective function, the aggregation index had the greatest influence on the value of the objective function. This was mainly because the aggregation index reflects the spatial distribution pattern of the stand, and an adjustment of the aggregation index meant that the spatial distribution pattern of the stand changed, which in turn led to changes in all spatial structure parameters, which explained why the aggregation index has the widest influence on the other parameters and why it has the greatest influence. In addition, although the spatial structure of the stand had been optimized, there were still certain problems, partly due to the disadvantages of a single structure of the plantation at the early stage of planting and partly because the optimization of the spatial structure of the stand could not be solved in one go [63], and it was often necessary to establish multiple objective functions to carry out multiple optimizations to achieve an overall improvement in the spatial structure of the stand. However, even so, this study still had strong applicability in planted forests by establishing a stand selective thinning model and optimizing the spatial structure of forest stands from four different aspects.

Finally, this study took the perspective of optimizing the spatial structure of forest stands to achieve management of plantation forests, but there were still some limitations in practice. Firstly, this study combined UAV-LiDAR and backpack LiDAR to obtain parameters of standing trees in the study area. While this method overcame the time-consuming nature of manual work and provided more accurate forest parameters than single-source LiDAR estimates, it did not accurately identify diseased and dead trees, which might result in diseased trees being retained and healthy trees being harvested during the selection process for harvesting. In addition, only the spatial structure of the stand was selected as an optimization objective in this study, while some important non-spatial structural parameters, such as the spacing distribution of diameter at breast height, were neglected, which explained why the spacing distribution of diameter at breast height was not significantly improved throughout the study area.

5. Conclusions

This study combined both UAV-LiDAR and backpack-LiDAR data to extract forest structure parameters. The use of dual LiDAR data overcame the shortcomings of a single data source and could provide a more reliable guarantee for the accurate extraction of forest structure parameters. Secondly, based on intelligent and accurate extraction of forest parameters, the objective function value Q was improved by 25.40% by constructing a dynamic multi-objective optimization model of spatial structure, indicating that the method could be effectively applied to the management of plantation forests. Finally, from the perspective of structured forest management, the direct contribution of openness K to the objective function value Q reached 51.30%, and the indirect contribution of the aggregation index R to the objective function value Q reached 51.50%, which indicated that the openness K and aggregation index R are the objects that needed to be adjusted first in the process of stand spatial structure optimization.