Optimizing Forest Spatial Structure with Neighborhood-Based Indices: Four Case Studies from Northeast China

: The ﬁne-scale spatial patterns of trees and their interactions are of paramount importance for controlling the structure and function of forest ecosystems; however, few management techniques can be employed to adjust the structural characteristics of uneven-aged mixed forests. This research provides an accurate, e ﬃ cient, and impersonal comprehensive thinning index ( P -index) for selecting candidate harvesting trees; the index was proposed by weighting the commonly used quantitative indices with respect to stand ﬁne-scale structures, competition status, tree vigor, and tree stability. The applications of the proposed P -index in evaluating and simulating the process of thinning operations were examined using four 1-ha mapped plots with di ﬀ erent forest types, namely, natural secondary forest, natural pine-broadleaved mixed forest, natural larch-birch mixed forest, and natural oak forest, which were widely distributed across the Heilongjiang Province in Northeast China. The results indicated that the proposed P -index could e ﬀ ectively a ﬀ ect the structural di ﬀ erentiations between di ﬀ erent forest types and alternative thinning intensities. The marginal beneﬁts of alternative thinning intensities on the integrated forest structure indicated that removing 10% of the trees from the plots might be the optimal thinning intensity from the perspective of optimizing stand structure, in which the P -index values could be increased by approximately 5%–11% for the four tested plots. The main conclusion from This paper was that the proposed P -index could be used as a quantitative tool to manage uneven-aged mixed forests. mixed larch-Betula forest, the Korean pine-broadleaved forest, the natural secondary forest, and the natural oak forest, respectively.


Introduction
Forest structure, including spatial and nonspatial characteristics, is one attribute contributing to the sustainable management of mixed uneven-aged forests. The nonspatial forest structure usually describes the average states of forest attributes, such as the mean diameter at breast height (DBH), mean tree height (HT), and stand density, which have been widely used in forest resource surveys, management, and monitoring [1]. Spatial forest structures exist at a variety of spatial scales and are of paramount importance for determining habitat and species diversity, which can be roughly divided into four levels: alpha, beta, gamma, and delta [2,3]. In forest ecosystems, the gamma and delta levels usually operate on landscape scales or at large forest areas, reflecting the fragmentation degrees of habitat and landscape; in contrast, the beta level refers to the variation between forest stands, reflecting the substitutability degree of different species under habitat gradients (e.g., elevation and succession stage), and the alpha level usually operates within forest stands, mainly reflecting the species richness, relative abundance, and uniformity. Many studies have been carried out at landscape levels and to evaluate stand spatial structure characteristics according to the structure-based forest management (SBFM) strategy [7,25]; however, the applicability of these comprehensive indices to optimizing forest spatial structure has not been widely evaluated. The study by Li et al. [26] is an example that has utilized a comprehensive index for selecting candidate trees for harvest according to the bivariate distributions of the three spatial parameters (i.e., U-, M-, and D-indices). However, the nonspatial indices of individual trees were all severely absent within the above studies. In fact, stand spatial and nonspatial attributes are almost equally important in controlling the timber values and the function of forest ecosystem services. Therefore, integrating spatial and nonspatial stand indicators together into the process of optimizing forest spatial structure may be more practical from the perspective of forest management.
The purpose of our study was to put forward a comprehensive thinning index based on tree neighborhood-based spatial relationships and nonspatial individual growth status. The comprehensive index was then employed to implement parametric thinning simulations for four different forest types across the province of Heilongjiang in Northeast China. Our hypotheses were as follows: (1) the proposed comprehensive thinning index can reflect the structural differentiations between different forest types, (2) the SBFM strategy can improve the complexity and diversity of stand spatial structures at fine scales, and (3) there is an optimal thinning intensity that can effectively improve the stand spatial structure characteristics.

Case Study Areas
Our test stands included four different forest types across the entire Heilongjiang Province in Northeast China (Figure 1), which are the major forest types in Northeast China [27]. The first stand is a mixed larch-birch (Larix gmelinii (Ruprecht) Kuzeneva-Betula platyphylla Suk.) deciduous forest (Plot a) located in the Cuigang Forest Farm (123 • 20 -124 • 21 E, 52 • 16 -52 • 47 N), which is part of the Daxing'an Mountains [28]. The terrain in This area is of mainly low-to medium-elevation mountains that trend significantly from northeast to southwest, with many valleys and overlapping peaks. The elevation of This region varies between 180 and 1530 m a.s.l., and the average elevation is approximately 573 m. The climate belongs to a cold-temperate continental monsoon climate with typical short, warm, humid summers and long, cool, dry winters. The average annual temperature is −2.8 • C, with a maximum temperature of 40.6 • C and a minimum temperature of −52.3 • C. The average annual precipitation ranges between 450 and 500 mm, and precipitation occurs primarily in the summer. The periods of snow cover in forests are as long as five months, and the thickness of snow is close to 50 cm. The soils are classified as dark brown coniferous forest soils that are formed under the combined influences of heat and moisture in mixed forests; however, meadow and swamp soils are also common in regions with low elevations. Natural L. gmelinii forests are historically the primary top zonal vegetation; however, most have disappeared, since excessive cutting was carried out between 1950 and 2000. Large trees were frequently removed from the forests for structural timbers, leaving a mass of small and unhealthy trees. During the long period of natural succession after the cutting disturbance, some native broadleaf species have naturally regenerated in This forest community, gradually forming a mixed and uneven-aged larch-birch forest. The forest community diversities are somewhat low, mainly due to the cold temperature, and L. gmelinii, B. platyphylla, and Populus davidiana Nakai are the most common species in these mixed forests.
The second plot is a mixed Korean pine (Pinus koraiensis Siebold et Zuccarini)-broadleaved forest (Plot b) that was once the climax community of the Xiaoxing'an Mountains. The plot is located in the Danqinghe Experimental Forest Farm (129 • 11 -129 • 26 E, 46 • 32 -46 • 39 N). The reserve is a low-altitude area (190-1028 m), located approximately 45 km north of the Songhua River [16]. This area has a temperate continental monsoon climate, with a mean annual rainfall of 600 mm and a mean annual evaporation of 1250 mm. The mean annual temperature is 2 • C, and the minimum temperature is about −35 • C in the winter (January), when frozen snowfall accumulates to a depth of 20-30 cm. Dark brown forest soils also predominate across This region. The vegetation primarily consists of more than 20 tree species, classified as Xiaoxing'an Mountain flora, among which B. platyphylla and the "three great hardwoods" in Northeast China (e.g., Fraxinus mandshurica Rupr., Phellodendron amurense Rupr., and Juglans mandshurica Maxim.) predominate in valleys lower than 300 m a.s.l., while Quercus mongolica Fischer ex Ledebour and B. platyphylla gradually become more common at higher elevations (300-700 m). Abies fabri (Mast.) Craib and Picea asperata Nakai usually appear at altitudes above 700 m. The specific study site was located in the 27th compartment, which belongs to the core area of This forest farm. The study site experienced several degrees of selective cutting before 1997, followed by complete natural restoration. To optimize the structure and promote the growth of the forest, a slight under-tending treatment (e.g., removing less than 15% of smaller trees and only those with less economic value) was implemented in 2014.
Forests 2020, 11, x FOR PEER REVIEW 4 of 19 mean annual evaporation of 1250 mm. The mean annual temperature is 2 °C, and the minimum temperature is about −35 °C in the winter (January), when frozen snowfall accumulates to a depth of 20-30 cm. Dark brown forest soils also predominate across this region. The vegetation primarily consists of more than 20 tree species, classified as Xiaoxing'an Mountain flora, among which B. platyphylla and the "three great hardwoods" in Northeast China (e.g., Fraxinus mandshurica Rupr., Phellodendron amurense Rupr., and Juglans mandshurica Maxim.) predominate in valleys lower than 300 m a.s.l., while Quercus mongolica Fischer ex Ledebour and B. platyphylla gradually become more common at higher elevations (300-700 m). Abies fabri (Mast.) Craib and Picea asperata Nakai usually appear at altitudes above 700 m. The specific study site was located in the 27th compartment, which belongs to the core area of this forest farm. The study site experienced several degrees of selective cutting before 1997, followed by complete natural restoration. To optimize the structure and promote the growth of the forest, a slight under-tending treatment (e.g., removing less than 15% of smaller trees and only those with less economic value) was implemented in 2014. The remaining two sites are the natural secondary forest (Plot c) and natural oak (Q. mongolica) forest (Plot d), located in Maoershan Forest Farm (127°18′-127°41′E, 45°2′-45°18′N), which covers an area of 26,496 ha, with a south-north length of approximately 30 km and an east-west width of 20 km [29]. It is located in a low-altitude area of the Zhangguangcai Mountains, approximately 80 km from Harbin city, the capital of Heilongjiang Province. The climate is a temperate continental monsoon climate with a mean annual temperature of 3.1 °C, a mean annual precipitation of 723 mm, a mean annual evaporation of 1094 mm, and a mean annual cumulative temperature (i.e., ≥ 10 °C) of 2283 °C. The frost-free period is approximately 120-140 days, and the relative humidity is 75%. The zonal soil is dark-brown forest soil, while azonic soils such as albic, meadow, and pear soils are common in seasonal or perennial catchment areas. The forestlands are mainly covered by natural regenerated secondary forests where the primary Korean pine-broadleaved forests have been severely damaged. Since ownership was transferred from the state to Northeast Forestry University as a trial forest farm in 1952, a large portion was subjected to complete restoration, gradually forming highly mixed natural secondary forests with more than 20 tree species. The specific study The remaining two sites are the natural secondary forest (Plot c) and natural oak (Q. mongolica) forest (Plot d), located in Maoershan Forest Farm (127 • 18 -127 • 41 E, 45 • 2 -45 • 18 N), which covers an area of 26,496 ha, with a south-north length of approximately 30 km and an east-west width of 20 km [29]. It is located in a low-altitude area of the Zhangguangcai Mountains, approximately 80 km from Harbin city, the capital of Heilongjiang Province. The climate is a temperate continental monsoon climate with a mean annual temperature of 3.1 • C, a mean annual precipitation of 723 mm, a mean annual evaporation of 1094 mm, and a mean annual cumulative temperature (i.e., ≥ 10 • C) of 2283 • C. The frost-free period is approximately 120-140 days, and the relative humidity is 75%. The zonal soil is dark-brown forest soil, while azonic soils such as albic, meadow, and pear soils are common in seasonal or perennial catchment areas. The forestlands are mainly covered by natural regenerated secondary forests where the primary Korean pine-broadleaved forests have been severely damaged. Since ownership was transferred from the state to Northeast Forestry University as a trial forest farm in 1952, a large portion was subjected to complete restoration, gradually forming highly mixed natural secondary forests with more than 20 tree species. The specific study site was located in the 36th compartment for the natural secondary forest, with a slope that was shady (north-facing), downward, and gentle (13 • ), while the 16th compartment for the natural oak forest had a slope that was sunny (south-facing), medium, and steep (23 • ). To optimize species composition and promote the growth of the natural secondary forest, approximately 500 P. koraiensis trees per hectare were replanted within the forest gaps in 2004, and This was followed by complete restoration. The tree species richness is relatively higher, and F. mandshurica, P. davidiana, Ulmus pumila Linn., P. amurense, P. koraiensis, and B. platyphylla are dominant. For the natural oak forest, the study site has experienced little anthropogenic disturbance (e.g., cutting and replanting) over the last 60 years; however, the structural complexity and species richness are both lower than those of the natural secondary forest, mainly due to the significant differences between their terrains.

Plot Measurements
From June to September in 2016 and 2017, four 100 m × 100 m permanent plots were selected and estimated from the three forest farms described. To minimize the closure error of the plot area (i.e., keep it < 1/400) and the measurement errors of the positions of trees, the entire plot was divided into 100 subplots (10 m × 10 m). In each plot, the precise positions of any tree with a diameter at breast height (DBH; 1.3 m) ≥ 5 cm were recorded in three dimensions. Other tree characteristics, such as species, DBH, HT, crown width, height of living branches, and health status (e.g., survival and the presence of diseases or insects), were also recorded. The spatial distributions and basic statistical characteristics of all surveyed trees for the four plots are illustrated in Table 1.

Data Analysis
Reasonable forest structures, including reasonable spatial and non-spatial structural characteristics, are usually regarded as the basis for realizing sustainable forest management. Therefore, a set of ecological-and economic-oriented indicators were carefully selected to develop a composite harvesting index, which was further employed to optimize stand structure. It should be noted that all of the selected indices in This study are relative indicators (i.e., ratio indicators), which help to prevent obvious selection priority tendencies (e.g., the preferential harvesting of large or small trees and the priority harvesting of broad-leaved trees or conifers) for harvesting trees of the proposed composite thinning index.

Nonspatial Structural Parameters
The stand nonspatial structural characteristics were evaluated using two indices: one was related to tree vigor, and the other was related to tree stability. The tree crown is the main site for a series of important physiological activities of forests and trees, such as photosynthesis and respiration; thus, the size, structure, and distribution of the canopy can comprehensively reflect the growth vigor and development trend of trees. Previous studies have indicated that the crown size usually increases significantly with tree diameter, presenting a typical sigmoid tendency [29,30]. To eliminate the preference of tree vigor indices, the ratio of crown width with respect to DBH was employed to evaluate tree vigor: where V i represents the ratio of crown width with respect to the DBH of tree i; CW i represents the mean crown width of tree i and was calculated as the mean value of crown width in the east (CW i e ), west (CW i w ), south (CW i s ), and north (CW i n ) directions; and DBH i represents the diameter at breast height of tree i. This index presents a significant exponential decay trend with increasing DBH and usually varies from 0 to 0.5 for the main tree species in Northeast China [31]. Trees with higher vigor have index values near 0.5, while trees with lower vigor have index values smaller than 0.1. The stability of trees against snow, ice, and wind, which is known as the slenderness coefficient, was calculated by dividing the total height by the DBH of the tree [32]: where S i is the height to diameter ratio of tree i, which usually has a range from 0 to 1; HT i is the total height of tree i; and DBH i is the diameter at breast height. The relationship between S and DBH was also a significant exponential decay function. In general, higher values of S indicate a higher position of the centre of gravity of trees with longer crown lengths, but a lower stability than trees with smaller S values [33,34]. Additionally, the S-index has significant effects on the mechanical properties of wood [35], i.e., trees with a higher S usually have a smaller maximum bending moment than do trees with a smaller S value.

Spatial Structural Parameters
To quantitatively determine the potential of harvesting trees from the perspective of optimizing forest spatial structural characteristics, four fine-scale spatial structural parameters were selected and employed, including the mingling index (M) [3,7,13], dominance index (D) [15], uniform angle index (U) [3,7,12], and Heygi competition index (H) [36]. In This paper, a spatial structural unit was defined as the combination of an arbitrary reference tree and its four nearest neighbor trees, as implemented in Programming [3], Li et al. [16], and Bettinger and Tang [37]. To eliminate the edge effects, a buffer area width of 5 m was employed; thus, the trees located in the core area were treated as reference trees, and the corresponding four parameters were calculated, while the other trees located in the buffer area were treated only as neighboring trees [16]. The M refers to the proportion of different species between the reference tree and its four nearest neighbors, reflecting the local-scale species diversity [13]. The D characterizes the diametric (or tree height) differentiation between a reference tree and its four nearest neighbors, which was put forward by Zhao et al. [15]. The U refers to the horizontal distribution pattern of the reference tree and its four nearest neighbors [12]. If the U values fall in the interval of 0.475-0.517, it implies a random distribution pattern [38]. Otherwise, the distribution pattern can be regular (U < 0.475) or clumped (U > 0.517). The above described parameters (i.e., M, D, and U) all include five possible values: 0.00, 0.25, 0.50, 0.75, and 1.00, which represent five different ecological grades for different aspects in the forest (Table 2). Obviously, values close to 0 indicate a lower level of species mingling, an absolutely inferior tree growth status, and a very regular distribution pattern of trees in the horizontal dimension, while high values document high species diversity, predominant growth status, and clumped patterns. The competition index (H), which might be regarded as the modified version of Hyegi's index [36], also refers to the competitive status, based on the properties (i.e., diameter and distance), between the reference tree and its four nearest neighbors. This index is a continuous variable, and the values may theoretically range from 0 to +∞. The formulations of the four spatial structural indices are stated as: where M i , D i , U i , and H i represent the mingling degree, diametric differentiation, spatial pattern, and competition pressure for reference tree i, respectively; n is the number of neighbor trees for any reference tree; and v ij , k ij , and z ij are discreteness variables and are set as different values according to the following criteria: (1) v ij = 1, if reference tree i and its neighbor tree j are of different tree species, otherwise, v ij = 0; (2) k ij =1, if neighbor tree j is smaller than the reference tree i, otherwise, v ij = 0; and (3) z ij = 1, if the angle α of two neighbor trees (e.g., the angle of neighbor trees j and j-1) is smaller than the expected standard angle α 0 (α 0 = 72 • ) [39], or otherwise z ij = 0. sp i and sp j represent the tree species of reference tree i and its neighbor tree j, respectively; dbh i and dbh j represent the diameter at breast height of reference tree i and its neighbor tree j, respectively; and L ij represents the Euclidean distance between reference tree i and its neighbor tree j. Details about the essential characteristics and classification criteria of these parameters are illustrated in Table 2.

Comprehensive Thinning Index
The quantitative spatial and non-spatial indices described above all have clear biological significance and can be rapidly assessed in the field, making the selection of candidate trees possible by evaluating the relationship between each reference tree and its four nearest neighbors. Since the magnitudes and measuring units of these indices are all markedly different, the comprehensive thinning index was formulated in utility theoretic as follows [40]: where P i is the comprehensive thinning index of tree i, N is the number of selected quantitative indices of tree growth status and spatial structural characteristics, w k and f k () are the relative importance (or weight) and sub-utility function of the k-th index, respectively, and x ik is the quantity of the k-th index for tree i. Obviously, smaller values of P for any reference tree indicate a higher probability of the tree being harvested. As emphasized here, the relative importance of objectives (w k () in Equation (7)) should be considered carefully as a part of the comprehensive harvesting index design; however, little guidance can be drawn from the literature about how to assign the weights of each index in SBFM. Therefore, two different weighting coefficient schemes were employed for testing purposes. The first scheme assumed that all objectives had equal importance, i.e., that all w k values in Equation (7) were 0.1667 (precisely 1/6). However, the weights of all considered indices within the second scheme were assumed to be 0.10 for the dominance index; 0.15 for tree vigor, tree stability, and uniform angle index; 0.20 for the mingling index; and 0.25 for competition pressure. In addition, the intensity of selective cutting is also an inherent part of any forest harvest scheduling model; thus, the effects of four different levels of selective cutting intensity (i.e., removing 10%, 20%, 30%, or 40% of the trees) on the stand spatial structural characteristics were quantitatively evaluated. The intensities of 20% and 30% reflect the present situation of forest management practices in Northeast China [41,42], while the intensities of 10% and 40% are somewhat under-or over-estimated, respectively, compared to the actual ranges.
The sub-utility functions (f k () in Equation (7)) for tree vigor (V), mingling index (M), and dominance index (D) increased linearly, so that a quantity (x ik in Equation (7)) equal to zero gave a sub-utility of zero, while the maximum possible quantity of the goal variable produced a sub-utility of one ( Figure 2a). For tree stability (S) and competition pressure (H), a linear decreasing function was employed, in which the maximum possible quantity of each index corresponded to a sub-utility of 1.0, and the minimum possible quantity produced a sub-utility of zero ( Figure 2b). The sub-utility of the uniform angle index (U) increased linearly from zero to one when the quantity increased from zero to 0.5 and then decreased to zero when the quantity further increased to one (Figure 2c). As emphasized here, the minimum, mean, and maximum values were calculated for each plot separately. China [41,42], while the intensities of 10% and 40% are somewhat under-or over-estimated, respectively, compared to the actual ranges. The sub-utility functions (fk() in Equation (7)) for tree vigor (V), mingling index (M), and dominance index (D) increased linearly, so that a quantity (xik in Equation (7)) equal to zero gave a sub-utility of zero, while the maximum possible quantity of the goal variable produced a sub-utility of one (Figure 2a). For tree stability (S) and competition pressure (H), a linear decreasing function was employed, in which the maximum possible quantity of each index corresponded to a sub-utility of 1.0, and the minimum possible quantity produced a sub-utility of zero ( Figure 2b). The sub-utility of the uniform angle index (U) increased linearly from zero to one when the quantity increased from zero to 0.5 and then decreased to zero when the quantity further increased to one (Figure 2c). As emphasized here, the minimum, mean, and maximum values were calculated for each plot separately.  (7)) for tree vigor, mingling index, and dominance index (a); for stability and competition pressure (b); and for uniform angle index (c).

Results
The statistical characteristics of different spatial parameters for the four forest types with alternative selective thinning intensities (i.e., removing 0%, 10%, 20%, 30%, and 40% of the trees) are  (7)) for tree vigor, mingling index, and dominance index (a); for stability and competition pressure (b); and for uniform angle index (c).

Results
The statistical characteristics of different spatial parameters for the four forest types with alternative selective thinning intensities (i.e., removing 0%, 10%, 20%, 30%, and 40% of the trees) are illustrated in Tables 3 and 4  The reason for This phenomenon might be that some larger broad-leaved trees have been protected from harvesting during recent decades. The comprehensive thinning index (i.e., P-index) is also a logical negative indicator; thus, the higher the values of the P-index are, the lower the urgency with which the forest stand must be harvested in consideration of optimizing the spatial structure. The management urgency of Plot d was significantly higher than that of the other three plots, regardless of which weighting strategies were employed.
With respect to the equal weights for all considered indicators (Table 3), the non-spatial variables of DBH and HT and the spatial indices of D and M all increased significantly for the four plots with an increase in selective thinning intensity. This result also corresponded to an increase of approximately 0.0243 cm for DBH, 0.0131 m for HT, 0.0004 for the D-index, and 0.0035 for the M-index when the thinning intensity was increased by 1% in Plot a. However, the spatial indices of U and H both decreased significantly with increases in the selective thinning intensity for all the plots, implying that a decrease of approximately 0.0003 for the U-index and 0.0246 for the H-index were observed in Plot a when the thinning intensity was increased by 1%. For the spatial indices of S and V, two typical inverse U-shaped tendencies were observed with the gradual increases in selective thinning intensity, in which the maximum values of S and V usually appeared with the intensity of removing 10% of the trees. As emphasized here, the values of increase or decrease for all the considered indices might be different in the other three plots; however, the variation tendencies were all consistent with those of Plot a. The mean values of the P-index logically increased with the selective thinning intensity for all four plots, which corresponded to an increase of approximately 0.0015 in Plot a, 0.0018 in Plot b, 0.0019 in Plot c, and 0.0022 in Plot d when the respective thinning intensities were increased by 1%. The relative increased proportion (RIP) of the P-index between any two successive selective thinning intensities decreased significantly with the increase in selective thinning intensity, indicating that removing 10% of the trees might be the optimal intensity from the perspective of optimizing the forest stand spatial structure. Note: 1) The number of trees within the core area for each plot. 2) The relative increased proportion of the comprehensive harvesting index (i.e., P-index in Equation (7)) between any two successive selective thinning intensities. Table 4. The statistical characteristics (Mean (Std)) of different optimized parameters for the four forest types with alternative selective thinning intensity, in which unequal weights for each index were employed.

Plot Intensity Number/Trees 1) DBH/cm HT/m D-index U-index M-index S-index V-index H-index P-index RIP /% 2)
Plot a  1) The number of trees within the core area for each plot; 2) The relative increased proportion of comprehensive harvesting index (i.e., P-index in Equation (7)) between any two successive selective thinning intensities.
Compared with the equally weighted strategy, the quantity values of each indicator might be different with the use of unequal weights for the four plots with alternative selective thinning intensity. However, similar variation tendencies were observed for most of the considered indicators, namely, the variables of DBH and M both increased, and the variables of U and H both decreased with increasing selective thinning intensity, while an inverse U-shaped tendency of the variable S was observed for all of the plots. Taking Plot a as an example, the above results also corresponded to an increase of approximately 0.0183 cm in the DBH and 0.0046 in the M-index, and a decrease of approximately 0.0001 in the U-index and 0.0226 in the H-index, when the thinning intensity was increased by 1%. The variation tendencies of the other three indictors (i.e., HT, D-, and V-indices) were somewhat complicated among different plots; however, the values of these indictors with the implementation of alternative selective thinning were all significantly larger than with no harvest. The mean values of the P-index also increased with alternative selective thinning intensity, which corresponded to an increase of approximately 0.0017 in Plot a, 0.0021 in Plot b, 0.0022 in Plot c, and 0.0025 in Plot d. The RIP values further indicated that the optimal intensity of selective thinning from the perspective of optimizing the forest stand spatial structure was 10%. More importantly, it should be noted that the values of the increases in the P-index and RIP were both higher than those when equal weights were employed.
The frequency distribution of assigned harvest trees with respect to the entire core area for different tree species within the four tested plots is illustrated in Figure 3. For natural mixed larch-Betula forests, the harvest possibilities of L. gmelinii were significantly higher than those of other tree species, and were as much as 6.07% greater than those for lower harvesting intensity (i.e., removing 10% of the trees) and 28.83% higher than those for higher harvesting intensity (i.e., removing 40% of the trees) when equal weights were employed. Similarly, significantly higher harvest possibilities of P. koraiensis, namely, 6.88% of that for lower harvesting intensity and 16.09% of that for higher harvesting intensity, were observed for natural mixed Korean pine-broadleaved forest. With respect to natural oak forests, Q. mongolica was the predominant harvest tree species, accounting for approximately 6.88% of that for the lower harvesting intensity and 16.09% of that for the higher harvesting intensity. However, the distribution pattern of harvest tree species of the natural secondary forest was more complicated than that in the other three plots, in which the harvest possibility of F. mandschurica varied from 2.85% to 13.41%, that of P. davidiana varied from 3.17% to 8.13%, that of P. koraiensis varied from 0.21% to 8.66%, and that of U. pumila varied from 1.58% to 4.86%, with an increase in harvesting intensity. As emphasized here, the numerical value of the frequency distribution of the assigned harvest tree species might be different when unequal weights were employed; however, the final verdicts were all consistent with the equal weights. The spatial distribution pattern of the assigned harvest trees for the four studied plots when the harvesting intensity was 10% and unequal weights were employed, in which the improvements in the spatial structural characteristics are evident and clearly discernible with the naked eye, is illustrated in Figure 4.  1) The number of trees within the core area for each plot; 2) The relative increased proportion of comprehensive harvesting index (i.e., P-index in Equation (7)) between any two successive selective thinning intensities.
Compared with the equally weighted strategy, the quantity values of each indicator might be different with the use of unequal weights for the four plots with alternative selective thinning intensity. However, similar variation tendencies were observed for most of the considered indicators, namely, the variables of DBH and M both increased, and the variables of U and H both decreased with increasing selective thinning intensity, while an inverse U-shaped tendency of the variable S was observed for all of the plots. Taking Plot a as an example, the above results also corresponded to an increase of approximately 0.0183 cm in the DBH and 0.0046 in the M-index, and a decrease of approximately 0.0001 in the U-index and 0.0226 in the H-index, when the thinning intensity was increased by 1%. The variation tendencies of the other three indictors (i.e., HT, D-, and V-indices) were somewhat complicated among different plots; however, the values of these indictors with the implementation of alternative selective thinning were all significantly larger than with no harvest. The mean values of the P-index also increased with alternative selective thinning intensity, which corresponded to an increase of approximately 0.0017 in Plot a, 0.0021 in Plot b, 0.0022 in Plot c, and 0.0025 in Plot d. The RIP values further indicated that the optimal intensity of selective thinning from the perspective of optimizing the forest stand spatial structure was 10%. More importantly, it should be noted that the values of the increases in the P-index and RIP were both higher than those when equal weights were employed.   Figure 3. The quantity distribution of harvested trees with respect to the entire core area for different tree species within the four tested plots, in which Plots a-d represent the mixed larch-Betula forest, the mixed Korean pine-broadleaved forest, the natural secondary forest, and the natural oak forest, respectively; symbols W1 and W2 represent the equally and unequally weighted scenarios, respectively; and C1-C4 represent the different harvesting intensities (i.e., removing 10%, 20%, 30%, or 40% of the trees).
The frequency distribution of assigned harvest trees with respect to the entire core area for different tree species within the four tested plots is illustrated in Figure 3. For natural mixed larch-Betula forests, the harvest possibilities of L. gmelinii were significantly higher than those of other tree species, and were as much as 6.07% greater than those for lower harvesting intensity (i.e., removing 10% of the trees) and 28.83% higher than those for higher harvesting intensity (i.e., removing 40% of the trees) when equal weights were employed. Similarly, significantly higher harvest possibilities of P. koraiensis, namely, 6.88% of that for lower harvesting intensity and 16.09% of that for higher harvesting intensity, were observed for natural mixed Korean pine-broadleaved forest. With respect to natural oak forests, Q. mongolica was the predominant harvest tree species, accounting for approximately 6.88% of that for the lower harvesting intensity and 16.09% of that for the higher harvesting intensity. However, the distribution pattern of harvest tree species of the natural secondary forest was more complicated than that in the other three plots, in which the harvest possibility of F. mandschurica varied from 2.85% to 13.41%, that of P. davidiana varied from 3.17% to 8.13%, that of P. koraiensis varied from 0.21% to 8.66%, and that of U. pumila varied from 1.58% to 4.86%, with an increase in harvesting intensity. As emphasized here, the numerical value of the frequency distribution of the assigned harvest tree species might be different when unequal weights were employed; however, the final verdicts were all consistent with the equal weights. The spatial distribution pattern of the assigned harvest trees for the four studied plots when the harvesting intensity was 10% and unequal weights were employed, in which the improvements in the spatial structural characteristics are evident and clearly discernible with the naked eye, is illustrated in Figure 4. . The quantity distribution of harvested trees with respect to the entire core area for different tree species within the four tested plots, in which Plots a-d represent the mixed larch-Betula forest, the mixed Korean pine-broadleaved forest, the natural secondary forest, and the natural oak forest, respectively; symbols W1 and W2 represent the equally and unequally weighted scenarios, respectively; and C1-C4 represent the different harvesting intensities (i.e., removing 10%, 20%, 30%, or 40% of the trees).

Discussion
The spatial patterns of trees and their interactions are of paramount importance in controlling the competition status, seedling growth, and survival and crown formation of trees, and knowledge of these patterns contributes to improving our understanding of the history, functions, and future potential of a particular forest ecosystem. It is almost impossible to maintain or generate a complex and reasonable stand structure promoting stand productivity, biodiversity, and stability without considering the detailed spatial characteristics of the forest. Therefore, how to utilize thinning techniques to improve the rationality and diversity of stand fine-scale stand structure has become the main concern of forest managers. The analysis presented in This paper successfully constructed a comprehensive thinning index by weighting the widely used spatial and nonspatial structure indices. The index was then employed to evaluate and simulate the effects of different thinning intensities on four different forest types across the province of Heilongjiang in Northeast China. The obtained results indicated that the proposed comprehensive thinning index can effectively affect the structural differentiations between different forest types and alternative thinning intensities and therefore can be used as a quantitative tool to manage uneven-aged mixed forests.
Our first hypothesis, that the proposed comprehensive thinning index can reflect the structural differentiations between different forest types and alternative thinning intensities, appears to be perfectly supported by our results. Since the four studied plots were extracted from different climatic regions (Figure 1), the differentiations between their spatial and nonspatial structural characteristics were obvious. With almost similar stand basal areas (i.e., 20 m 2 /ha), their stem numbers and the sizes of individuals differed significantly (Table 1), in which larch-birch mixed forests had approximately twice as many trees as the others. The differences among these plots were also reflected clearly by the number of tree species, which is another traditional measure of forest structure. For the two more novel nonspatial indices (i.e., the Sand V-indices), Plot b obtained the smallest stability value and the largest vigor value, simultaneously, when compared with the other stands, indicating the unique stand vigor and stability characteristics of pine-broadleaved mixed forests, which have suffered some anthropogenic disturbances during recent decades.
Spatial attributes provided further information about the plots. The differentiations of tree neighborhood-based spatial parameters (i.e., U-, M-, and D-indices) were mainly reflected by the mixed species compositions, in which the tree size differentiations were all very uniform, and the horizontal patterns of the trees all exhibited clumped distributions, while the mingling values differed significantly from each other. Plots b and c were both stands with strongly mixed compositions; however, Plots a and d had relatively lower mingling values, both belonging to the moderate mixed composition (Tables 3  and 4). The significant positive relationships between species richness and mingling values implied that the complexity of the neighborhood-based stand structure of the four tested plots may highly depend on the traditional measures of biodiversity, which is consistent with the results from previous studies on the comparisons of stand structure with different species compositions [26,43]. In addition, since a slight decreasing tendency of the values of the U-index with increasing stand density was observed, the number of trees per hectare might be another traditional measure of stand structure that can contribute to improving the complexity of fine-scale stand structure; however, its contribution was far less than that of the M-index. The phenomena were also detected in natural Scots pine (P. sylvestris) [44] and artificial larch (L. principis) [25] forests, both indicating that a single tree species with a large number of individuals will significantly reduce the species mingling degree and slightly increase the uniform angle degree. With respect to the distance-dependent measure, the Heygi competition index indicated that the trees grown in Plot c were suffering significantly more competition pressures than the others, mainly due to the existence of some large trees (42 trees per hectare of DBH≥30 cm) and coppice trees that shared the same (or similar) coordinates with each other (Figure 4). The differences in integral stand structure between different plots can be distinguished clearly using the proposed comprehensive thinning index. Additionally, because of the longer natural succession process (approximately 60 years) and proper human disturbance (i.e., planted P. koraiensis under canopy) in Plot c, the P-index values were significantly larger than those of Plot b, which was converted from the original pine-broadleaved mixed forests through repeated intensive selective cuttings.
Our second hypothesis that the SBFM strategy can improve the complexity and diversity of stand spatial structure at fine scales was also perfectly supported by our results. Unlike the thinning priority index provided by Gadow et al. [23], Pastorella and Paletto [7], Li et al. [26], and Ye et al. [25], which considered only part (or all) of the stand structure characteristics based on tree neighborhood relationships, the comprehensive thinning index proposed in This paper included four critical aspects with respect to stand fine-scale structures, competition status, tree vigor, and stability. The simulated results indicated that the thinning treatment based on the SBFM strategy significantly affected the spatial structure of the four tested plots to different degrees. With increasing thinning intensities, D and M both increased; however, U and H both decreased, while S and V fluctuated with intensity. These results indicated that thinning could increase the tree dominance and mixing degree and decrease the aggregation extent of the horizontal distribution and competition pressure of reference trees, which was consistent with the results from some successful practices [16,23,26]. However, the detailed response mechanisms among tree vigor, stability, and thinning intensities were more complex when the SBFM strategy was implemented, and therefore, these require further research. The differences in the P-index between different thinning intensities for most of the simulation cases were also statistically significant (Tables 3 and 4). As emphasized here, the benefits of using the neighborhood-based parameters are that the multifaceted relations among trees are clearly defined and can be investigated very quickly and precisely. Although the distance-dependent measure of the Heygi competition index requires detailed tree coordinates in space, which usually require very time-and cost-consuming data collection methods, it can be estimated using the universal and stable power functions between the DBH and H-index of individuals, a method that has been verified for a wide range of stand characteristics [45,46]. Therefore, the application of our proposed comprehensive thinning index can not only improve the speed and accuracy of field surveys, but also reduce the subjectivity of the choice of trees to be harvested in uneven-aged mixed forests.
The final hypothesis, that there would be an optimal thinning intensity that could effectively improve the stand spatial structure characteristics, was not supported. The integral forest structure characteristics (i.e., P-index) all increased significantly, although some attributes may have been sacrificed. The P-index values increased as much as 10% or more when compared with those with no harvest, while the rates of these increases (i.e., RIP) decayed gradually with increasing thinning intensities in all four tested plots. The maximum RIP values were all observed when 10% of the trees were removed from the studied plots, indicating that removing 10% of the trees might be the most effective thinning intensity in terms of improving the integral stand structure characteristics, especially for the areas that have suffered repeated intensive selective cutting and in which the valuable trees have almost been depleted. The intensity is also coincidently satisfied with forest management policies and regulations from the State Forestry Bureau of China [47], in which harvesting for timber purposes from natural forests has been strictly prohibited, instead of allowing low-intensity tending operations to improve the quality of forest ecosystems. The assigned harvest within This intensity mainly focused on the trees with low health, lower vigor (e.g., small crown), and high competition status, while only a handful of trees were harvested for the purpose of adjusting the neighborhood-based structure ( Figure 4). However, the optimal intensity from the perspective of stand structure optimization was far less than that for the purpose of promoting stand growth and yield (45% in [48]; 38%-52% in [49]), species diversity (44% in [50]), carbon sequestration and stocks (45% in [51]), and soil physical and chemical properties (30% in [52]). Therefore, the effects of removing 10% of the trees on various aspects of forest ecosystems, and the differences compared with other intensities, require more research.
The proposed comprehensive thinning index could be considered as accurate, efficient, and impersonal for selecting thinning trees in practice; however, there may be some technical reasons for the hinderance of direct application. The first, and perhaps the most important, was that two subjective weights (i.e., equal and unequal weights) were employed in the six portions of the objective function for all four plots, without considering the details of stand characteristics, management objectives, and the willingness of other stakeholders. Therefore, when one goal is weighted differently from the others, the results may also be biased. The second important technical reason was that the information regarding the vertical dimension for tree distribution was not also considered, which is also ecologically important [53,54]; thus, it would also be interesting to explore whether the inclusion of vertical structure such as the stand layer index (L) [53] or other similar indices could improve the simulation of optimizing stand structure. The last, but not least, reason is that the time and money consumed by mapping the locations of the trees usually increased geometrically with increasing sampling size.

Conclusions
Following the guidance of the SBFM strategy, the case study presented here proposed an accurate, efficient, and impersonal comprehensive thinning index for the selection of candidate harvesting trees by weighting the commonly used quantitative indices with respect to stand fine-scale structure, competition status, tree vigor, and tree stability. The applications of This index in evaluating and simulating the process of thinning operations for the four widely distributed forest types across the province of Heilongjiang in Northeast China indicated that the proposed comprehensive thinning index could effectively affect the structural differentiations between different forest types and alternative thinning intensities. The ranks of the diversity and complexity of integral forest structure characteristics were as natural secondary forest (0.600 and 0.638), pine-broadleaved mixed forest (0.588 and 0.623), larch-birch mixed forest (0.566 and 0.597), and natural oak forest (0.552 and 0.585) when equal and unequal weights were implemented, respectively. The marginal benefits of alternative thinning intensities on the integral forest structure showed that removing 10% of the trees from the stands was the optimal thinning intensity for optimizing the integral stand structure characteristics, in which the comprehensive thinning index could be improved by approximately 5%-11% in the four test plots. However, forest management is a long-term project in nature, and one optimization cannot solve all issues; thus, the processes of the optimization, monitoring, and evaluation of forest ecosystems should be implemented cyclically.