A New Tree-Level Multi-Objective Forest Harvest Model (MO-PSO): Integrating Neighborhood Indices and PSO Algorithm to Improve the Optimization Effect of Spatial Structure

Abstract

Accurate, efficient, impersonal harvesting models play a very important role in optimizing stand spatial structural and guiding forest harvest practices. However, existing studies mainly focus on the single-objective optimization and evaluation of forest at the stand- or landscape-level, lacking considerations of tree-level neighborhood interactions. Therefore, the study explored the combination of the PSO algorithm and neighborhood indices to construct a tree-level multi-objective forest harvest model (MO-PSO) covering multi-dimensional spatial characteristics of stands. Taking five natural secondary forest plots and thirty simulated plots as the study area, the MO-PSO was used to simulate and evaluate the process of thinning operations. The results showed that the MO-PSO model was superior to the basic PSO model (PSO) and random thinning model Monte Carlo-based (RD-TH), DBH dominance (DOMI), uniform angle (ANGL), and species mingling (MING) were better than those before thinning. The multi-dimensional stand spatial structure index (L-index) increased by 1.0%~11.3%, indicating that the forest planning model (MO-PSO) could significantly improve the spatial distribution pattern, increase the tree species mixing, and reduce the degree of stand competition. In addition, under the four thinning intensities of 0% (T1), 15% (T2), 30% (T3), and 45% (T4), L-index increased and T2 was the optimal thinning intensity from the perspective of stand spatial structure overall optimization. The study explored the effect of thinning on forest spatial structure by constructing a multi-objective harvesting model, which can help to make reasonable and scientific forest management decisions under the concept of multi-objective forest management.

1. Introduction

Forest planning mainly consists of three stages: assessment, planning, and monitoring [1,2]. It is a complex task highly dependent on mathematical planning and information technology [3], which provides a decision-making framework for forest managers: when, where, how, and how much to harvest [4]. Modern forest planning involves the coordination of multi-function relationships. The trade-offs between forest ecosystem services supply and forest management costs may especially complicate the decision-making process [5,6]. Therefore, combinatorial optimization is the key to optimizing stand structure, maintaining species diversity, improving forest productivity, and increasing habitat heterogeneity [7,8,9].

Characterization of stand spatial structure is an essential part of sustainable forest management [10,11,12]. Among them, neighborhood indexes that consider the location of individual trees and interactions between neighboring trees have received widespread attention. For instance, DBH (diameter at breast height) dominance, species mingling, and uniform angle are distance-dependent, neighborhood-based indices. DBH dominance is widely used to describe tree competition, Species mingling is used to evaluate the degree of tree species mixing, and uniform angle represents tree spatial distribution pattern [11]. Some studies evaluated forest spatial structure using neighborhood indices and developed corresponding strategies to increase the diversity of tree species, optimize stand spatial distribution pattern, and reduce the degree of competition, thus enhancing forest quality [11]. Unfortunately, few studies use the multi-dimensions neighborhood indices to develop forest planning strategies.

Thinning, as a silvicultural treatment in forest planning, can promote forest growth, improve stand structure, maintain biodiversity, and prevent pests and fires [12]. Trees should be removed, and strategies should be taken in tree-level forest planning. These are the key elements, as the thinning process is limited by multiple factors. Previous studies focused on constructing thinning models using stem quality attributes [13], tree species diversity [11], adjacencies between trees [14], and economic criteria [15,16] as the criteria for selecting harvested trees, which have certain limitations, including the ignoration of the potential impact of the spatial distribution characteristics [17,18]. Multi-dimension optimization objectives often conflict with each other, showing nonlinearity and discontinuity characteristics. Therefore, multi-objective optimization algorithms are needed in tree-level forest planning.

Recently, heuristics are increasingly being used to solve forest planning problems [19,20]. For instance, Fotakis et al. [21] proposed the constrained NSGA-II, which is used for multi-objective forest planning problems, including maximization of timber volume, minimization of division levels, and even the flows and minimum timber yield. Bettinger and Tang [11] combined the threshold accepting heuristic and the species mingling index to achieve tree-level operational harvest scheduling optimization. Compared with other algorithms, the particle swarm optimization (PSO) algorithm has significant advantages in solving large-scale mathematical optimization problems by using memory and feedback mechanisms: (1) fast calculation speed and strong global search ability; (2) faster convergence than evolutionary algorithms and genetic algorithms; (3) simple parameter setting, strong local search ability, and low possibility of falling into local optimum [22,23]. Multi-objective forest planning is also a similar optimization process: calculate the changes of comprehensive indices before- and after-harvesting under the constraints of multiple management objectives and in a given solution space (e.g., trees in the plot) to determine the optimal harvesting intensity.

However, the application of PSO in tree-level forest planning has not been explored in depth. Therefore, from the perspective of tree-level multi-objective forest planning, it is a novel concept and method to combine particle swarm optimization (PSO) with multi-dimensional stand structure indices to construct a harvest model. It contributes to the expansion and extension of sustainable forest management. Based on five natural secondary forest plots in South Dongting Lake and thirty simulated plots, the following studies were conducted: (1) develop a tree-level multi-objective forest planning model by integrating PSO algorithm and neighborhood indexes (DOMI, ANGL, and MING). (2) explore the response of species diversity, competition, and forest spatial distribution pattern to forest thinning, as well as the optimal thinning intensity.

2. Materials and Methods

2.1. Data

2.1.1. Natural Plots

The natural plots were located in the Dongting Lake region in Hunan province, China. Dongting Lake is the second largest freshwater lake in China, located in the middle of the Yangtze River (28°39′ N–30°14′ N, 111°42′ E–113°93′ E). The zonal vegetation in the area was the subtropical evergreen broad-leaved forest in the middle subtropics. Cyclobalanopsis glauca and Schima superb were dominant tree species. The frost-free period was 263~277 d. The mean annual rainfall was 1400~1500 mm. The mean annual temperature was 16.6~17.0 °C, the lowest temperature in January was 4.1~4.5 °C, and the highest temperature in July was 28.7~29.2 °C.

We prepared five natural plots, which were located in Lutou Forest Farm (LT), Daweishan Nature Reserve (DW), Longhushan Forest Farm (LH), Dashanchong Forest Farm (DS), and Wuyunjie Nature Reserve (WT). The area of each plot was 20 m × 30 m. The tree height, DBH (≥5 cm), tree position, tree crown width, and tree species of all surveyed trees were recorded. The tree species were recorded and provided by forest farm staff. The area of each plot was 20 m × 30 m. The tree height, DBH (≥5 cm), tree position, tree crown width, and tree species of all surveyed trees were recorded. The DBH was measured using a ruler. Tree height was measured using a laser altimeter. The crown width was measured with a tape measure and the spatial position (x,y) of each tree was recorded using RTK. The spatial distributions of five natural plots and the statistical characteristics of trees are illustrated in Figure 1 and

Table 1. Basic statistical characteristics of the five studied plots.

Site	DW	LH	DS	LT	WT
Altitude (m)	1300	68	217	335	713
Slope (°)	47	35	27	15	38
Aspect	E	W	ES	EN	S
Canopy density	0.75	0.65	0.70	0.80	0.60
Mean DBH (cm)	12.7	10.3	14.7	11.3	9.4
Mean Height (m)	11.3	9.1	13.8	9.7	7.1
Mean Crown (m)	2.4	2.7	2.0	3.1	2.5
Number of species	8	7	10	6	8

.

2.1.2. Simulated Plots

The study sample plots included two parts: simulated plot and natural plot. The spatial distribution pattern of simulated plots was divided into three types: uniform distribution (A), random distribution (B), and aggregated distribution (C). Each type contained ten plots, for a total of 30 plots (Figure 2). In our study, Python (version 3.6.0) was used to generate 30 plots for each type, numbered 1–30, with 30–70 trees per plot, and a plot area was 20 m × 30 m, the dominant species were Cyclobalanopsis glauca and Schima superb. The standard deviation of tree size (DBH, tree height, crown) in each group increased as the number of groups increased: the minimum standard deviation of tree size of group 1 was 0.014, and for group 30, the standard deviation of tree size was at 0.613, which was the largest.

2.2. Neighborhood Indices

In order to accurately demonstrate the relationship of neighboring trees, such as the degree of tree species mixing, tree competition, and tree spatial distribution pattern, we compiled the widely used and representative stand spatial structure indices (

Table 2. Stand spatial structure index.

Index	Calculation Formula	Variable Definition
Uniform angle (ANGL)	$W_i=\frac{1}{n}{\displaystyle {\displaystyle \sum }_{j=1}^n}{Z}_{ij}$	$W_i$ is the uniform angle of central tree i, n is the number of neighboring trees, $Z_{ij}$ is the variable of uniform angle. When the angle between central tree i and neighboring tree j is less than the standard angle, $Z_{ij}$ = 1, otherwise, $Z_{ij}$ = 0.
DBH Dominance (DOMI)	$U_i=\frac{1}{n}{\displaystyle {\displaystyle \sum }_{j=1}^n}{k}_{ij}$	$U_i$ is the neighborhood comparison of central tree i, $k_{ij}$is the value variable of neighborhood comparison. When the DBH of neighboring tree j is smaller than that of central tree i, $k_{ij}$ = 0, otherwise, $k_{ij}$ = 1.
Species mingling (MING)	$M_i=\frac{1}{n}{\displaystyle {\displaystyle \sum }_{j=1}^n}{v}_{ij}$	$M_i$ is the species mingling of central tree i, $v_{ij}$ is the value variable of the mingling degree. When the central tree i and neighboring tree j are the same trees, $v_{ij}$ = 0, otherwise $v_{ij}$ = 1.

). Hui et al. proposed the “Structure-based Forest Management” theory. Among them, the division of forest spatial structure units by “1 + 4” had been recognized by the international public and widely used. According to the recommendation and the “1 + 4” criterion, the number of neighboring trees n selected in this paper was 4 in our study [24,25]. The index distribution of five natural plots was shown in Figure 3.

2.3. Construction of Dynamic Multi-Objective Optimization Model of Forest Spatial Structure

The multi-objective optimization of forest spatial structure is a dynamic optimization problem. Each goal is independent and affects the others. It is necessary to consider the tree species diversity, the inter- and intra-specific competition, and the tree distribution to make the overall forest structure healthy and stable. Therefore, no matter which forest management mode is chosen, stand spatial structure should be evaluated comprehensively instead of depending on a single structure index. Considering that tree species mixing, competition, and distribution pattern of trees are the main influencing factors of forest spatial structure, three objective functions of MING, DOMI, and ANGL are selected as the constraints of the evaluation method. In addition, the stand spatial structure optimization strategy first adjusts the horizontal distribution pattern of stands, and judges the horizontal distribution pattern of stands according to the average value of uniform angle index. It is generally believed that if the stand is not seriously disturbed, the horizontal distribution pattern should be random distribution after long-term forest development and succession [26,27,28]. Therefore, when making forest harvested strategy, we should first analyze the uniform angle index of the stand, and give priority to adjusting the horizontal distribution pattern of the stand from non-random distribution to random distribution.

The purpose of constructing the multi-objective optimization model is to optimize the stand structure, so the output is the target trees that affect the stand structure. The trees with small MING, large DOMI, and large ANGL significantly affect the overall stand structure [29,30], which needs to be removed to optimize the stand structure. Our study was based on utility theory [31,32], with the constraints of maximizing MING, minimizing DOMI and ANGL, and then constructing a multi-dimensional stand spatial structure index (L) to evaluate the overall stand spatial structure after optimization. The L is constructed as follows:

$$\mathrm{maximize}{M}_i=\frac{1}{n}{{\displaystyle \sum }}_{j=1}^n{v}_{ij},$$

(1)

$$\mathrm{minimize}{U}_i=\frac{1}{n}{{\displaystyle \sum }}_{j=1}^n{k}_{ij},$$

(2)

$$\mathrm{minimize}{W}_i=\frac{1}{n}{{\displaystyle \sum }}_{j=1}^n{Z}_{ij},$$

(3)

$$L_i=f\left(M_i,W_i,U_i\right),$$

(4)

$$L=\frac{1}{N}·\frac{{\omega }_m}{{\omega }_w·{\omega }_u}{\sum }_{i=1}^N \frac{(1+M_i)·(1+{\sigma }_{M_i})}{(1+W_i)·(1+{\sigma }_{W_i})·(1+U_i)·(1+{\sigma }_{U_i})},$$

(5)

where $M_i$, $W_i$ and $U_i$ are the species mingling, uniform angle, DBH dominance of the target tree $i$, respectively; $N$ is the total number of trees in the stand; $n$ is the number of the neighbors of the target tree $i$; ${\sigma }_{M_i}$, ${\sigma }_{W_i}$ and ${\sigma }_{U_i}$ are the standard deviation of species mingling, uniform angle, DBH dominance of the target tree $i$, respectively; ${\omega }_m$, ${\omega }_w$, and ${\omega }_u$ are the weights corresponding to the three optimization constraint objectives, respectively; L represents the overall level of stand spatial structure. The smaller L, the simpler the composition of tree species, the greater the competition among trees in certain areas, and the spatial distribution of stands is closer to non-random distribution.

The more ideal the stand spatial structure, the higher the degree of tree species mixing, the lower the degree of competition, and the more random the spatial distribution of trees. In addition, the thinning intensity is also an essential part of the forest harvest scheduling model. Therefore, quantitative evaluation of the impact of different harvesting intensities on the stand spatial structure characteristics plays an important role in the early stages of forest development and succession. In this study, according to the Technical Regulations on the Investigation and Design of Forest Cutting Areas in Hunan Province, the cutting intensity should be considered according to the comprehensive factors such as the forest management objectives, site conditions, and stand conditions (http://www.hunan.gov.cn, accessed on 2 December 2022). The tending cutting intensity was less than 20%~25%; the release cutting intensity was less than 40%, and the accretion cutting intensity was less than 30%. The cutting intensity varied in different regions. The plots in our study were in Hunan Province. In order to meet the actual thinning intensity and situation, we set four different levels of cutting intensity, namely 0% (T1), 15% (T2), 30% (T3), and 45% (T4), to explore the impact of harvesting intensity on the stand spatial structure.

2.4. Construction of Tree-Level Multi-Objective Forest Harvest Model (MO-PSO)

2.4.1. Particle Swarm Optimization

Particle swarm optimization (PSO) is an efficient and powerful population-based stochastic search technique for solving global optimization problems. PSO algorithm abstracts each solution of the problem to be solved as a particle in the solution space. The solution process is actually an iterative optimization process in the solution space by a particle swarm composed of $\mathrm{m}$ particles. Particles search for new positions by constantly updating their positions and speeds and determining whether the searched positions are optimal by using the fitness function.

$$v_{id}=\omega v_{id}+c_1{r}_1(p_{id}-x_{id})+c_2{r}_2\left(p_{gd}-x_{id}\right),$$

(6)

$$x_{id}=x_{id}+v_{id}+min\sqrt{{\left(x_i-x^{\prime }\right)}^2+{\left(y_i-y^{\prime }\right)}^2},$$

(7)

where $X_i=\left(x_{i1},x_{i2},\cdots ,x_{iD}\right)$ represents the position of the i^th particle in a D-dimensional search space, $V_i=\left(v_{i1},v_{i2},\cdots ,v_{iD}\right)$ is the velocity of i^th particle, $P_{\mathit{besti}}=\left(p_{i1},\dots ,p_{iD}\right)$ is i^th particle’s personal best position, $P_{\mathit{gbest}}=\left(p_{g1},\dots ,p_{gD}\right)$ is global best position. Acceleration constants $c_1$, $c_2$ and inertia weight $\omega$ are predefined by the user, and $r_1$ and $r_2$ are the uniformly generated random numbers in the range (0, 1).

When the particle position is updated according to Equation (7), it is possible that the particle position is not within the definition domain. Such particles are called transboundary-particles. When the transboundary-particles appear, we multiply the speed in Equation (7) by the adjustment coefficient $\beta$ to pull the transboundary-particles back into the definition domain. The new location formula is shown in Equation (8). If the particles remain outside the boundary after adjustment, reduce the $\beta$ value according to Equation (9) until the particles return to the decision space.

$$x_k^i=x_k^i+\beta v_k^i+min\sqrt{{\left(x_i-x^{\prime }\right)}^2+{\left(y_i-y^{\prime }\right)}^2},$$

(8)

$$\beta =1/\left({\delta }^2+1\right),$$

(9)

where $\delta$ is the number of adjustments.

2.4.2. The Method of MO-PSO Model Construction

In this study, each tree in the forest was regarded as a solution in the solution space, the forest space corresponded to the target solution space of the PSO algorithm, and the spatial coordinates of trees in the forest corresponded to the position of particles, and the Function was taken as the fitness function, thus transforming the multi-objective planning problem of stand spatial structure into an iterative optimization process of the particle swarm in the solution space.

For an ordinary particle swarm, N particles are randomly initialized, and then these particles are divided into $M$ sub-swarm by K-means algorithm. In each sub-swarm, the particle with the shortest distance from other particles in the particle swarm is selected as the central particle of the particle swarm. The central particle is calculated by the following formula:

$$x_k^{\mathit{Center}}=arg\left(min{{\displaystyle \sum }}_{j=1}^{\frac{N}{M}}\parallel x_k^i-x_k^j\parallel \right), i\in \left[1,\frac{N}{M}\right],$$

(10)

$$dis\left(x_k^i,x_k^j\right)=\parallel x_k^i-x_k^j\parallel =\sqrt{\left({{\displaystyle \sum }}_{l=1}^n{\left(x_k^{il}-x_k^{jl}\right)}^2\right)/n},$$

(11)

in which $x_k^{\mathit{Center}}$ represents the central position of each particle swarm, dis is the distance between i^th particle and j^th particle, $x_k^i$ and $x_k^j$, respectively, represented as i^th particle and j^th particle in the particle swarm, and $n$ is the dimension of the decision variable $x$. If the distance between the $x_k^i$ and other particles in the same swarm is the smallest, then the central particle of the $k$^th swarm $x_k^{\mathit{Center}}=x_k^i$.

After determining the center particle of the $M$ sub-swarm, the distance between the sub-swarm $k$ and other sub-swarms is calculated, respectively, and the dynamic particle group is constructed according to the results obtained. The upper limit and lower limit of the distance between particle swarms are $D_{max}$ and $D_{min}$, respectively, and the particle swarm with the smallest distance from the particle swarm $k$ is a particle swarm a. If dis$\left(x_k,x_a\right)>D_{max}$, a new particle swarm $x_{M+1}$ is generated. If the minimum distance between a particle swarm and another particle swarm is less than $D_{min}$, a particle swarm is deleted. The calculation formula of the newly generated particle swarm is:

$$x_{M+1}^{il}=\left(x_k^{il}+x_a^{il}\right)/2+c_1{\left(-1\right)}^{\mathit{round}\left(0.5+c_2\right)}\left|x_k^{il}-x_a^{il}\right|/2,$$

(12)

where $x_{M+1}^{il}$ represents the l dimension component of the i^th particle in the particle swarm $M+1$, $c_1,c_2$ is the random number in $\left[0,1\right]$, and round $\left(\cdot \right)$ represents the rounding function.

The particle swarm searches for trees that meet the optimization objectives in the stand space using a parallel mechanism, and the final output of the algorithm is the target tree that needs to be regulated or logged in the stand. The construction process of the multi-objective optimization harvesting model (MO-PSO) integrating neighborhood indices and PSO algorithm is shown in Figure 4 and Figure 5. The simulation of the tree-level harvesting process was performed by c++ and Python.

Step 1: Import spatial position information and individual tree attribute data of all trees $Q$, such as DBH, tree height, and crown width—import data.

Step 2: Randomly generate an initial population $P$ with the size of N. Use the $\mathrm{K}$-means algorithm to divide the randomly generated particles into $M$ sub-swarm, and randomly generate the initial position and velocity of each particle—initial population.

Step 3: Calculate the Euclidean distance ${\left\{D_i\right\}}_{i=1}^Q$ from each particle to the tree, adjust the position of the particle to make it fall on the nearest tree, and calculate the individual fitness value according to the initial particle—Calculate distance.

Step 4: Initialize the external solution set $p_k^i$ of $i$^th particle in particle swarm $k$ to be empty, the external set $L_k$ of particle swarm $k$ to be empty, and the external set $G$ of the global optimal solution of all particles to be empty—initialize the external solution set.

Step 5: Store the non-dominated solution in the external set and update the optimal position—update particle position.

(1) The current position of each particle and the current position of all particles in each particle swarm are stored in the external set of each particle and the external set of the swarm, respectively, according to the non-dominated solution rules, and the optimal position $pbes{t}_k^i\left(t\right)$ of i^th particle at the time $t$ and the optimal position $lbes{t}_k\left(t\right)$ of the particle swarm $k$ are calculated.

(2) The current positions of all particles are stored in the external set $G$ of the global particle swarm according to non-dominated rules, and then the optimal particle of the population is found from the external set $G$ of the global optimal solutions of all particles, and its position is the global optimal position $gbest\left(t\right)$.

Step 6: Update the speed and position of particles according to Equation (6) and Equation (7), adjust the position of particles, and calculate the current fitness value. If the new position exceeds the boundary of the definition domain, the out-of-boundary particles shall be treated according to Equation (8) and Equation (9)—calculate the current fitness value.

Step 7: Calculate the distance between particle swarm, insert or delete particle swarm according to the method described in Section 2.4.2—insert or delete particle swarm.

Step 8: For the inserted new particle swarm, initialize its particles and generate particle external set and group external set—update new particle swarm.

Step 9: Judge the change of L-index. If the L-index does not change, stop searching and output the target tree. Otherwise, return to Step 6 and continue to search for the optimal solution—judge the value of L-index.

3. Results

3.1. Model Performance

The simulated harvesting results in

Table 3. The model performance for the MO-PSO model, PSO model, and the RD-TH model.

Distribution Pattern	MO-PSO			PSO			RD-TH
	AITE	ALIV	$\overline{RIP}$/%	AITE	ALIV	$\overline{RIP}/\%$	RUNN	ALIV	$\overline{RIP}/\%$
Uniform	26.9	0.697	31.94	34.2	0.664	30.34	100,000	0.631	27.01
Random	37.4	0.702	29.51	49.6	0.693	27.72	100,000	0.676	24.36
Aggregated	43.9	0.604	24.71	54.3	0.597	20.33	100,000	0.554	19.77

RIP: the relative increased proportion; $\overline{\mathrm{RIP}}$: average of relatively increased proportion; AITE: average iterations; ALIV: average of L-index; RD-TH: random thinning model Monte Carlo-based; RUNN: random cutting times.

show that the MO-PSO model was superior to the basic PSO model (PSO) and random thinning model Monte Carlo-based (RD-TH). For the uniform distribution, random distribution, and aggregated distribution stand, under the MO-PSO model, the average of the relative increased proportion ($\overline{\mathrm{RIP}}$) were as large as 31.94%, 29.51%, and 24.71%, respectively, which were slightly larger than PSO model. However, they were far larger than RD-TH model. At the same time, the average iterations of MO-PSO (AITE) were about 26~43 times, which was far less than RD-TH (10,000 times simulating). As far as MO-PSO and PSO models, when the MO-PSO model integrated the multi-swarm objectives strategy, the ALIV and AITE were also evidently better than PSO. Obviously, our proposed approach has advantages in multi-dimensional stand spatial structure optimization.

In scheduling individual tree harvests in our study, compared with the conventional single-objective optimization model, the MO-PSO model performed better in improving the stand multi-dimensional spatial structure before- and after-thinning. As shown in Figure 6 and

Table 4. The statistical characteristics (mean (std)) for the five plots with alternative thinning intensity.

Plot Code	Intensity	DBH/cm	MING	RIP/%	DOMI	RIP/%	ANGL	RIP/%	L-Index	RIP/%
DW	0%	12.710 (3.143)	0.473 (0.223)		0.590 (0.351)		0.560 (0.259)		0.612 (0.225)
	15%	13.107 (2.967)	0.491 (0.201)	3.81	0.425 (0.332)	−27.97	0.484 (0.246)	−13.57	0.671 (0.203)	9.60
	30%	13.412 (2.762)	0.502 (0.198)	2.24	0.336 (0.303)	−20.94	0.442 (0.246)	−8.68	0.719 (0.213)	7.19
	45%	13.430 (2.784)	0.515 (0.190)	2.59	0.303 (0.278)	−9.82	0.407 (0.224)	−7.92	0.749 (0.197)	4.15
LH	0%	10.360 (3.853)	0.572 (0.210)		0.485 (0.350)		0.550 (0.275)		0.594 (0.198)
	15%	11.797 (3.521)	0.589 (0.199)	2.97	0.431 (0.338)	−11.13	0.514 (0.265)	−6.55	0.661 (0.176)	11.32
	30%	11.903 (3.326)	0.603 (0.192)	2.38	0.342 (0.314)	−20.65	0.474 (0.257)	−7.78	0.708 (0.188)	7.07
	45%	12.001 (3.117)	0.612 (0.189)	1.49	0.290 (0.283)	−15.20	0.447 (0.245)	−5.70	0.732 (0.172)	3.33
DS	0%	6.710 (2.976)	0.514 (0.223)		0.499 (0.357)		0.529 (0.250)		0.607 (0.215)
	15%	7.286 (2.853)	0.539 (0.192)	4.86	0.447 (0.346)	−10.42	0.498 (0.240)	−5.86	0.672 (0.197)	10.63
	30%	7.074 (2.793)	0.548 (0.189)	1.67	0.350 (0.315)	−21.70	0.470 (0.246)	−5.62	0.708 (0.201)	5.44
	45%	7.295 (2.707)	0.564 (0.167)	2.92	0.311 (0.292)	−11.14	0.452 (0.242)	−3.83	0.738 (0.193)	4.23
LT	0%	11.320 (3.194)	0.604 (0.231)		0.493 (0.349)		0.531 (0.265)		0.605 (0.208)
	15%	12.036 (3.013)	0.630 (0.203)	4.30	0.445 (0.339)	−9.74	0.493 (0.252)	−7.16	0.672 (0.184)	11.02
	30%	12.213 (2.936)	0.641 (0.204)	1.75	0.357 (0.316)	−19.78	0.449 (0.249)	−8.92	0.715 (0.194)	6.45
	45%	12.501 (2.903)	0.657 (0.189)	2.50	0.312 (0.291)	−12.61	0.431 (0.246)	−4.01	0.722 (0.177)	1.01
WT	0%	8.430 (3.723)	0.511 (0.231)		0.506 (0.356)		0.504 (0.274)		0.615 (0.225)
	15%	8.989 (3.237)	0.533 (0.206)	4.35	0.456 (0.349)	−9.88	0.467 (0.262)	−7.34	0.681 (0.208)	10.75
	30%	9.035 (2.898)	0.540 (0.197)	1.16	0.370 (0.319)	−18.86	0.423 (0.262)	−9.42	0.721 (0.211)	5.85
	45%	9.127 (2.836)	0.549 (0.194)	1.83	0.312 (0.292)	−15.68	0.409 (0.261)	−3.31	0.746 (0.205)	3.49

RIP: the relative increased proportion; DBH: diameter at breast; W-index: uniform angle; U-index: DBH dominance; M-index: species mingling.

, with the number of algorithm iterations increasing (Figure 6), the L index of five plots increased significantly, ranging from 1.0% to 11.3%. Under different thinning intensities, the average of MING increased by about 2.77% and showed an upward trend in all plots after thinning, and the variation in DS was the highest (about 4.86%). ANGL and DOMI decreased compared with those before thinning, and all five plots showed the same trend.

In addition, we also conducted simulated thinning on 30 simulated plots with uniform distribution (A), random distribution (B), and aggregated distribution (C). As shown in Figure 6, in the process of algorithm solution, the L index of A, B, and C all increased. The magnitude of fitness value (L index) change in A and B stands was slightly larger than that in C. The optimization time of A and B was shorter. The result showed that the optimized stand spatial structure was significantly better than before, and the PSO algorithm was feasible in solving the multi-objective optimization problem of stand spatial structure.

3.2. Thinning Intensities

As shown in Figure 7 and Table 4, the MING increased significantly with the thinning intensity increases for all the natural plots. However, the DOMI and ANGL values for all plots between two successive thinning intensities decreased significantly with the increase in thinning intensity. For the relatively increased proportion (RIP) of MING, DOMI, and ANGL, three similar inverse U-shaped tendencies were observed with the gradual increases in thinning intensity, L-index represents a multi-dimensional comprehensive stand spatial structure index, covering MING, DOMI, and ANG. L-index increased and the maximum values of the RIP of L-index usually appeared under the intensity of T2, indicating that the overall stand spatial structure was affected to varying degrees under different thinning intensities, and removing 15% of the trees (T2) was the optimal thinning intensity from the perspective of stand spatial structure overall optimization in our study.

3.3. Competition, Structure and Spatial Distribution Pattern

To explore the effects of MING, DOMI, and ANGL on the selection of harvested trees, we analyzed the distribution of MING, DOMI, and ANGL of harvested trees. As shown in Figure 8, the spatial structure indices distribution of harvested trees is different among the five plots. Of all five levels of the MING index distribution ($M_i=0,0.25,0.5,0.75,1$), trees with lower values of MING index were more likely to be harvested ($M_i=0$ or $M_i=0.25$, that is, 1 or 0 of the nearest four neighboring trees are the same as the target tree), and the mean value of relative frequency was 33.7%. However, trees with higher or moderate MING index ($M_i=0.75$ or $M_i=1$) were relatively fewer. Similarly, trees with higher DBH dominance ($U_i=0.75$ or $U_i=1$) or the trees showed aggregated distribution ($W_i=0.75$ or $W_i=1$) were more likely to be harvested, and the mean values of the relative frequency of DOMI and ANGL were 47.3% and 35.6%, respectively. The results showed that the individual trees with higher DBH dominance, lower species mingling, and horizontal distribution aggregation were more likely to be harvested, indicating that these trees could have a negative influence on stand structure, tree species composition, or the intensification of tree competition.

In addition, changes in DOMI were the largest, indicating that thinning may significantly impact the competition among neighboring trees, followed by distribution pattern and tree species composition (Figure 9 and Table 4). Tree competition, tree diversity, and horizontal spatial distribution of trees should be fully considered in the forest spatial structure optimization.

4. Discussion

In the mixed forest, the stand spatial distribution pattern, tree species mixing. and competition are important criteria in the harvested tree selection, and they are key endogenous constraints to tree growth. For instance, same- and different-species of trees through asymmetrical competition for light or symmetrical competition for water and nutrients affect the stand spatial structure and its succession process [33,34]. However, the traditional indices describing competition, stand distribution pattern. and tree species mixing, such as total basal area, tree density, Weibull distribution. and Gini coefficient, are all distance-independent, and cannot accurately quantify the impact of neighborhood interaction on the stand spatial structure [35,36]. Therefore, in this study, three distance-independent neighborhood indices, DBH dominance (DOMI), uniform angle (ANGL), and species mingling (MING), were selected to describe the degree of stand competition, spatial distribution of forest, and tree species diversity, respectively. DOMI measures the proportion of n (n = 4 in our study) nearest neighbors larger than the target tree, ANGL is used to describe the angular uniformity of the n nearest neighbors around the target tree, and MING is used to measure the proportion of n (n = 4) nearest neighbors not belonging to the same species as the target tree [10]. In order to explore the impact of thinning on forest spatial distribution pattern, mixing and competition, this study attempted to combine the PSO algorithm with multi-dimensional stand spatial structure indices (MING, DOMI, and ANGL) to construct a tree-level multi-objective forest planning model (MO-PSO), which is used to evaluate and simulate the impact of different thinning intensities on different stands. The results indicated that the MING, DOMI, and ANGL of different stands indices were better than those before thinning under alternative thinning intensities, and the multi-dimensional stand spatial structure (L-index) had been improved to varying degrees. Therefore, the model can provide guidance and help for multi-objects forest planning of natural secondary mixed forests.

In addition, the results indicated that DOMI, ANGL, and MING were better criteria to optimize tree selection. Obviously, the competition between target trees and neighboring trees and the degree of forest aggregation are important factors affecting forest growth. Trees in high aggregation and fierce competition environments have smaller tree size characteristics (DBH, height, crown) than those in free-growing (or less competitive) environments [30]. In our study, trees with larger size neighboring trees (DOMI > 0.5), highly aggregated distribution (ANGL > 0.5), and lower tree species mixing degree (MING < 0.5) were more likely to be harvested. The simulation thinning analysis results also showed that thinning had a greater impact on the degree of stand competition. Under different thinning intensities, with the improvement of the overall stand structure (L-index increases), the individual trees were more evenly distributed. Our results were consistent with some previous studies [11,12,37,38].

In previous studies, when the complexity of forest spatial management planning problems increased, showing nonlinear or discontinuous characteristics, the multiple optimization objectives conflicted. Researchers proposed to combine the heuristic algorithm with a single structure index for a solution [11,39], but the processing capacity and application scope of this technology is seriously limited. Therefore, in our study, the PSO algorithm was used to construct a multi-objective forest planning model, which quickly and fairly sought approximate optimal solutions and excellent local search ability. In addition, the initialization particle swarm generated randomly by the PSO algorithm avoids the subjectivity of harvesting tree selection and makes the solution more objective. The results indicated that the model had a significant effect on the optimization of multiple objectives, such as improving the degree of tree species mixing and reducing the aggregation degree of tree horizontal distribution and competition pressure, which are consistent with previous studies [12,40,41]. However, forest planning is affected by various factors such as topography, climate, and soil properties. Due to the limitations of research conditions and data accessibility, this study only incorporated the stand spatial structure characteristics into the model. At the same time, the logging gaps will produce naturally regenerated tree species over time [11], thus the model is unsuitable for monitoring and analyzing the dynamics in forest ecosystems. Therefore, further perfecting the forest planning model to improve the applicability and evaluation accuracy of the model under dynamic conditions will be the core and difficulty of the subsequent research.

Thinning reduces the number of trees, decreases the competition for resources in limited supply (e.g., light, nutrients, and water) among trees and expands the living space of individual trees [42]. However, the optimal thinning intensity varied greatly under different optimization objectives. In this study, the multi-dimensional stand spatial structure index (L-index) increased significantly under all four thinning intensities, and the value of some forest structural attributes showed certain fluctuations. Compared to the unharvested, L-index increased by 23.7% or more, but the relatively increased proportion (RIP) increased and then decreased significantly with increasing thinning intensity. When 15% of the trees were removed from the studied plots, the RIP values all reached the maximum, indicating that removing 30% of the trees might be the most effective thinning intensity to improve the multi-dimensional stand spatial structure. However, other studies have shown that the optimal thinning intensity to increase species diversity is about 60%, sustainable wood production is about 25~40%, promote stand growth and yield are 45% and 38%~52%, and improve soil physicochemical properties and microbial community are about 30%–45% [43,44,45,46,47]. Compared with other intensities, the optimal thinning intensity of stand structure optimization in this study was much less than other intensities. Therefore, whether removing 30% of the trees is the optimal thinning intensity, its impact on various aspects of the forest ecosystem and its differences, and whether it complies with relevant local forest management laws and regulations need further research.

In addition, some studies reported that advanced remote sensing technologies such as airborne laser scanning (ALS), ground laser scanning (TLS), and digital aerial photogrammetry (DAP) have been widely used in forest inventory and management. Integrating ALS and DAP data (data fusion) will help habitat modeling or drawing the forest attribute map [48], which might be of great help in calculating DOMI, ANGL, and MING. Combined with the MO-PSO model in this study, real-time tree-level forest planning decision-making will be possible. Although there are still many challenges in integrating ALS and DAP data to quantify three-dimensional spatial structure and attributes of stands [49], taking into account the synergy among MO-PSO, ALS, and DAP are crucial for making tree-level decisions, which will help to maintain or increase species diversity, bird habitat composition diversity, and habitat suitability.

5. Conclusions

Under the guidelines of the structure-considered forest management strategy, we presented an accurate, efficient, and comprehensive thinning model for the selection of harvesting trees by integrating the PSO algorithm and stand spatial structure indices with species diversity, competition, and spatial distribution pattern. The simulated cutting results showed that the MO-PSO model was superior to the basic PSO model (PSO) and random thinning model Monte Carlo-based (RD-TH). The analysis of the effects of MING, DOMI, and ANGL index on the selection of harvested trees showed that the individual trees with higher DBH dominance, lower species mingling, and horizontal distribution aggregation were more likely to be harvested. Thinning especially had a significant influence on tree competition from neighboring trees. A comprehensive analysis of four thinning intensities on the overall forest structure showed that removing 15% of the trees is the optimal thinning intensity in our study. The multi-dimensional stand spatial structure index increased by approximately 1.0% ~ 11.3% in the five natural plots. However, forest management and planning is a continuous, multi-factorial, and cyclic optimization process. The stand structure characteristics and research scope covered in our study were limited. Therefore, further research will consider more factors, such as the environment, temperature, and terrain.