Tree-Level Harvest Optimization for Structure-Based Forest Management Based on the Species Mingling Index

This novel research investigated the use of a heuristic process to inform tree-level harvest decisions guided by the need to maximize the interspersion of tree species across a forest. In the heuristic process, a species mingling value for each tree was computed using both (1) neighbors that were simply of a different species than the reference tree and (2) neighbors that were uniquely different species from both the reference tree and other neighbors of the reference tree. The tree-level species mingling value was averaged for the stand, which was then subject to a maximization process. Constraints included residual tree density levels and minimum tree volume harvest levels. In two case studies, results suggest that the species mingling index at the stand level can be significantly increased over randomly allocated harvest decisions using the heuristic process described. In the case studies, we illustrate how this type of process can inform management decisions by suggesting the distance between residual trees of similar species given the initial stand structure and the objectives and constraints. The work represents a unique tree-level optimization approach that one day may be of value as new technologies are developed to map the location of individual trees in a timely and efficient manner. OPEN ACCESS Forests 2015, 6 1122


Introduction
The structure of forests, both the location of trees and their characteristics (marks), is a significant factor in the processes that regulate growth and promote sustainability [1] and, thus, has an impact on future stand development and ecological stability [2].Forest structure can be described by the arrangement of trees across a landscape and their associated characteristics [3], and the manipulation of forest structure suggests that the growth and health of forests can be altered, for better or worse.A consideration of structural features of stands during operational decision-making processes can perhaps facilitate or enhance the biological diversity of the forest [4].Structure-based approaches to forest management can also possibly increase forest quality through changes in species composition and structural diversity, and optimization of these approaches can thus possibly accelerate stand development, rather than allow gradual natural succession to control a stand's destiny [5].These types of approaches can also be used to transition an even-aged plantation to a structurally-complex forest or be used to manage mixed-species uneven-aged forests, as an understanding of spatial forest structure has been noted as a key to the management of uneven-aged forests [1].
In terms of measuring the interspersion of trees across a landscape, the species mingling index can be used to describe spatial diversity or the dissimilarity of tree species given a reference tree and a set of neighboring trees [4,6,7].It can be viewed as a surrogate measure of biodiversity that accounts for basic tree diversity considerations [8].The size of a tree (diameter, height) does not necessarily matter in the computation of the species mingling index unless a minimum threshold is assumed.For example, in some inventories, trees that are very small (seedlings with no diameter recorded and less than 1 m in height) may be included along with larger, more mature trees.However, some analyses using the mingling index may only involve the larger trees.At the tree level, the species mingling value is inversely proportional to the relative tree density of each tree species.At the stand level, the more homogenous the collection of trees, the lower the species mingling index becomes.The species mingling index has been used to identify the effects of structure-based approaches to forest management on stand spatial structure [5].Natural development of some forests in the absence of major disturbances can lead to increases in the species mingling index, and thus, the structural development of unmanaged forests might be used as a reference for planning [4].While generally applied to live trees, the mingling concept has also been applied to the presence of snags (dead trees) in an assessment of the natural condition of forests [9] and would seem to complement other resource measurements regarding the ecological characteristics of forests [10] and to help assess forest biodiversity [11].
Methods for simulating tree marking regimes for uneven-aged forest management practices using basal area levels and the variation in basal area and harvested tree diameters as constraints have also been described [12].Tang et al. [13] illustrated how the spatial structure for partial harvests of Phyllostachys edulis (Haubrich) (moso bamboo) in China can be improved through the use of Monte Carlo simulation.Further, Tscheschel and Stoyan [14] illustrated a process using simulated annealing to reconstruct a point pattern, which was further explored in [15] for reconstructing landscape-level individual tree distributions from sample plot distributions.Finally, the species mingling index has been optimized for the development of spatially diverse hypothetical forests, whereby trees of different species can swap locations [1].Rather than reconstruct tree point pattern distributions from sample information, we concentrate an initial distribution of trees and the selection of those trees to harvest (remove) while maximizing the species mingling index.The only tree-level harvest optimization process located during our literature review was one where a heuristic was used to select the trees to harvest from within a 0.25-ha uneven-aged forest plot located in Spain, in order to optimize utility (based on the market value of harvested trees and the residual forest structure), while adhering to wood-flow and gap (opening) size constraints [16].
The main objective of this paper is to demonstrate a tree-level optimization approach for uneven-aged management of forests in the United States that: (1) involves maximizing the species mingling index subject to minimum harvest levels and minimum residual basal area levels that are informed by a habitat suitability model; and (2) involves maximizing the species mingling index subject to minimum harvest levels and minimum residual basal area levels that are informed by a higher-level stand-level optimum management regime, which itself accounts for changes in stand density and certain operational constraints (minimum level of extracted volume).Given a spatial representation of trees across a landscape, the approach we describe can be used to develop an operational plan optimized to enhance the diversity and structural complexity of a forest.However, as we note in the Materials and Methods and Discussion of the paper, the initial spatial distribution of trees and information concerning the first-period harvest schedule are necessary requirements.One aspect of interest that the paper does not specifically address is whether the spatial arrangement of the final locations of the residual trees form point patterns that are in fact independent patterns, through the use of random labelling or random reallocation tests (as used in [17]).The results of such tests would indicate whether the patterns produced deviate from a hypothesized pattern of complete spatial randomness.

Calculation of the Species Mingling Index
In forestry, the species mingling index represents the proportion of neighboring trees to a reference tree that are not the same tree species as the reference tree [18].Therefore, the fundamental unit of information in the calculation is the point-related characteristic (a tree location and its species, size and status), rather than the test location on the ground.Gao et al. [5] described how the species mingling index could be applied to individual trees to produce qualitative statements of localized areas.For each tree within a stand, the species mingling value is calculated as: where i = the tree around which the analysis occurs.j = a tree from the set (Ni) of n closest neighboring trees to tree i. Mi = the mingling value of tree i. n = the number of closest neighbors to the reference tree that is assumed in the computation.Ni = the set of n closest neighboring trees to tree i. vij = 1 if speciesi ≠ speciesj; 0 otherwise.
The result is a value between 0 and 1 that represents no (0) or significant (1) mingling of tree species around the reference tree (i).For each tree, there are only (n + 1) possible values.It has been suggested that n = 4 was best for developing forest structure decisions [1,5,19,20].This, for example, helps one to understand five tree diversity scenarios from no mixing to strongly mixed.However, it seems that any number of neighbors (n) can be used to develop a stand-level average species mingling index.While often used alone to help assess the diversity of forest species within a forest [11], the species mingling index has been widely used as a nearest neighbor statistic to describe the spatial structure of a forest [7,21].The basic premise of the index, to assess whether marks of the locations of pairs of trees are within a certain distance, also informs second-order characteristics, such as the mark mingling function, which can help one determine the extent of the spatial correlation between forest characteristics, such as tree species interspersion or structural diversity [17].
If a harvest decision is associated with the species mingling index computation through a tree-level planning process, the index is calculated as: where ti = a binary integer representing the tree continuing to be alive (1) or being harvested (0) after a planned harvest operation.The average species mingling index for an entire stand is of interest here, as it was also in [5].For an entire stand within a forest, an average species mingling index could be used to describe the spatial dispersion of the tree species located there.where MDn = the species mingling index for a stand using n neighbors.T = the number of trees in the stand.

Case Study Areas
In order to illustrate the method, we apply it to two 1-ha forests described by both real and hypothetical tree locations.The first is an uneven-aged eastern oak-hickory (Quercus spp.-Carya spp.) deciduous forest located in the southern United States (Georgia), and the second is an uneven-aged mixed-conifer forest that may likely be found in the intermountain west (eastern Oregon) of the United States.

Eastern United States Oak-Hickory Forest
The eastern United States oak-hickory forest inventory is real, and the locations of the trees across a 1-ha area (Figure 1) were measured using a compass direction and a distance from an established control point.A summary of the inventory of the deciduous stand (Table 1) suggests that there are 687 trees per ha, amounting to 22.6 m 2 per ha of basal area.Over half of these (381 trees per ha) are small eastern hophornbeam (Ostrya virginiana Mill.) with diameters smaller than 12 cm (Figure 2).Oaks (Quercus spp.) represent about 55.3% of the basal area; hickories (Carya spp.) represent about 11.6% of the basal area; and eastern hophornbeam represents only about 9.0% of the basal area.In this stand, there are 229 trees per ha greater than 13 cm in diameter at breast height (DBH), and the total merchantable volume is slightly more than 245 metric tons per ha.The trees are represented geographically as a distribution of marks (species and diameters) that are known.Distances between trees were represented by the computed Euclidean distance, based on the x and y coordinates representing tree centers.The volumes were developed using the methodology for southern U.S. deciduous forests described in [22] and the methodology for southern U.S. coniferous forests described in [23].In this forest, a harvest was assumed to occur immediately, with the goals of maximizing tree species mingling and further promoting the development of an uneven-aged stand of oaks and hickories that white-tailed deer (Odocoileus virginianus Zimmermann) favor.In a habitat model for white-tailed deer [24], it has been suggested that the retention of a greater number of oak tree species and the retention of a greater oak basal area level will improve habitat quality.However, a minimum harvest level is required to ensure the feasibility of the treatment, and the minimum level that we assume is 90 metric tons of wood per ha from trees greater than 15 cm in diameter.Further, the amount of oak basal area is currently about 13.7 m 2 per ha, which is less than the optimal level described in [24] for white-tailed deer habitat, yet some of the oak volume needs to be harvested to meet the minimum harvest level, therefore, we are assuming that the minimum basal area to be retained of oaks is 11.5 m 2 per ha from trees that are 25 cm or greater in diameter, which represents the range of basal area required to achieve a 0.5 habitat suitability index level for this variable [24].Given this scenario, the tree-level problem that is the subject of this manuscript is designed to maximize the post-harvest tree species mingling index given the harvest level suggested.

Western United States Mixed Conifer Forest
The western United States mixed conifer forest inventory is real, yet the locations of the trees across a 1-ha area are hypothetical and randomly assigned.A summary of the inventory of the western mixed-conifer stand (Table 2) suggests that there is a larger number of very small Abies grandis (Dougl.ex D. Don) Lindl., Abies lasiocarpa (Hook.)Nutt.and Pinus contorta (Dougl.ex.Loud.) trees.Many of these were recorded as having diameters smaller than 2.5 cm and heights lower than 0.3 m.The 432 trees per hectare that represented trees greater than 13 cm in breast height diameter were randomly allocated to a one-hectare plot (Figure 3) using an independent marking procedure similarly described in [25], yet the distribution of marks (species and diameters) are known a priori from the inventory of the stand.Distances between trees were represented by the computed Euclidean distance between trees, based on the x and y coordinates representing tree centers.The tree list that represented this stand was projected 100 years into the future using the stand-level optimization model described in [26] that utilizes the methodology contained in the Forest Vegetation Simulator [27] growth and yield system.The goal of the stand-level harvest problem was to maintain as closely as possible the stand density to between 35% and 55% of the Reineke stand density index throughout the 100-year projection period (time horizon), subject to maintaining a minimum residual (post-harvest) basal area of 18.4 m 2 per hectare and requiring a harvest entry to remove at least 21 m 2 per hectare.Of the projected harvested trees in the first decade of the plan, 56 were Abies grandis, 38 were Abies lasiocarpa and were 152 Picea engelmannii (Parry ex.Engelm.).Given this scenario, a tree-level problem was designed to maximize the species mingling index during the first decade (post-harvest) given the harvest level suggested.

Optimization Method
From a formulaic point of view, the tree-level operational plan for a single stand that we illustrate is designed as a non-linear optimization problem that maximizes the species mingling index for the stand, subject to constraints on the type of variable employed and accounting for rows that facilitate the summation of the species mingling index.

Maximize MDn
(3) Subject to: where hi = a variable representing the harvest of tree i (1 = yes).Tb,c = the set of pre-harvest trees that represent species group b, tree diameter class c.Hb,c = an assumed number of harvested trees that arise from species group b, tree diameter class c. bai = the basal area (m 2 per hectare) represented by tree i. BA = the post-harvest basal area target for the stand.voli = the volume (m 3 per hectare) represented by tree i. VOLUME = the harvest volume target for the stand.
Unfortunately, because MDn is computed and optimized in a post-harvest manner, the number of neighboring trees necessary to compute Mi is not known a priori, (i.e., the members of the set Ni may not include the n closest trees to each reference tree, if some of these are scheduled for harvest), which precludes the use of standard or non-linear mixed integer programming methods that require each possible relationship (combination of neighboring live trees) to be constructed prior to solving the model.Therefore, a heuristic method (threshold accepting) was developed to generate near-optimal solutions to the problem.
Threshold accepting was developed initially by Dueck and Scheuer [28] and has been applied to forest-level optimization problems in [29][30][31][32], among others.Threshold accepting is an iterative search process for developing near-optimal solutions to problems.Typically, an initial feasible solution is randomly generated, and depending on how the heuristic has been developed, subsequent changes to the solution are made to either ensure feasibility or penalize objective function values for infeasibilities that occur.In our case (Figure 4), feasibility with respect to the constraints is always maintained as a new solution is developed.The heuristic process begins with the random scheduling of harvest of set Hb,c.For the eastern United States oak-hickory forest case study, the process then involves selecting one tree for a change in its current status (e.g., scheduled to unscheduled or vice versa, representing a 1-opt move).The constraints are assessed, and the species mingling index is computed.Since the western United States mixed conifer forest case study needs to adhere to a specific schedule of harvest for certain tree size classes, the process involves selecting two trees of the same species and diameter class, one tree that has been scheduled for harvest and one tree that has not.Their status is then switched (representing a 2-opt move).As with the eastern United States oak-hickory forest case study, the constraints are assessed, and the species mingling index is computed.For the western United States mixed conifer forest case study, the number and type of trees in Hb,c is prescribed, yet for the eastern United States oak-hickory case study, the number and type of trees in Hb,c is not.Therefore, Equation ( 7) is ignored in the eastern United States oak-hickory case study.
Compared to simulated annealing, a threshold accepting process will accept every solution created along the Markov chain that is either superior or not much worse than the previous solution that was located.The assessment of "not much worse" involves comparing the potential objective function value to a lower bound (in the case of a maximization process) determined by either (1) the best solution's objective function value minus the current threshold or (2) the previous solution's objective function value minus the current threshold.In our case, the best solution's objective function value is used in this analysis.Simulating annealing uses a probability function to determine whether to formally accept a worse solution rather than the deterministic method employed by threshold accepting.As the heuristic progresses, the threshold value contracts, forcing the process to only accept changes to the solution that are very close to (in terms of objective function value) the previous feasible solution.Therefore, depending on the initial threshold value, the process may allow significant diversification of the search process in the early stages and, depending on the termination criteria, may allow significant intensification of the search process in the late stages.In our implementation of threshold accepting to this process, the initial threshold (in terms of MDn) varies from 0.001 to 0.005, a range that is based on several trials.To illustrate the process and to assess the sensitivity of this parameter, we allow either 10 or 25 feasible iterations (changes to tree status) at each threshold level.The rate of change of the threshold, based on the previous threshold value, is either 0.99, 0.995 or 0.9975, and the process terminates when the threshold value is 0.00001 or less.Finally, to ensure that the heuristic process does terminate, the threshold is changed if 100 consecutive non-feasible choices are considered.Finally, in the development of the heuristic, a process was enabled to assess Mi for each tree in a manner where the neighbors simply had to be of a different species (and thus, where vij = 1) than the reference tree.However, a separate process was also enabled to assess Mi for each tree in a manner where the neighbors had to be uniquely different species from both the reference tree and other live residual neighbors of the reference tree.Both of these representations of the species mingling index are accommodated in this work.This approach is slightly different from the approach suggested in [19], where the Mi for each tree is multiplied by a ratio of the number of unique tree species around the reference tree over the total potential number of unique tree species in the stand.
Initially, 3000 random harvest schedules were generated for each case study stand, given the direction provided by the stand-level guidance.For each combination of threshold accepting search parameters and assumptions of species mingling, 100 solutions were then generated.In sum, for each case study, 6000 feasible solutions were generated using threshold accepting to illustrate the usefulness of the process for tree-level harvest optimization based on the species mingling index (2 types of mingling × 5 initial thresholds × 2 iterations per threshold × 3 threshold change rates × 100 runs each).Since some of the sets of 100 solution values developed with threshold accepting were not normally distributed, a non-parametric method, the two-tailed Mann-Whitney U-test, was used to determine whether unpaired data of the sample sets had significantly different median values using a significance level of 0.05.

Eastern United States Oak-Hickory Forest
In every result from using the tree-level optimization process, the minimum basal area (11.5 m 2 per ha from trees that are 25 cm or greater in diameter for oak species) constraint seemed to be the most binding of the two main constraints.Some flexibility was noticed in the minimum harvest volume required, as many solutions were 2%-3% above the minimum (90 metric tons of wood per ha from trees greater than 15 cm in diameter).The best results when maximizing the species mingling index using non-unique oak species neighbors seemed to be 0.4809 (Table 3).A map of the harvest plan for this case is illustrated in Figure 5. Four sets of threshold accepting parameters seemed to produce the best results in the case where non-unique oak species neighbors were used to compute the species mingling index.In each case, 25 feasible iterations were used at each threshold level, and the rate of change of the threshold was 0.9975.However, the initial thresholds were 0.001, 0.002, 0.003 and 0.005.When an initial threshold of 0.001 was employed, the mean solution value (0.4723) was the highest among all of the results from using the potential sets of threshold accepting parameters.With these four parameter sets, the greatest number of solutions with a value of 0.4809 was located: 8, 11, 13 and 11 for initial thresholds of 0.001, 0.002, 0.003 and 0.005, respectively.In addition, the greatest number of solutions with a value that was within 1% of 0.4809 was located with these four parameter sets: 37%, 41%, 39% and 34% for initial thresholds of 0.001, 0.002, 0.003 and 0.005, respectively.These were among the parameter sets that allowed threshold accepting to search longer (25 iterations and a slow rate of change) at each threshold level.According to the results of the two-tailed Mann-Whitney U-tests, of these four parameter sets, in only one case were two sets of unpaired sample data significantly different (p < 0.05), when the initial threshold was 0.002 and when it was 0.005.From this and previous results, one might conclude that using an initial threshold of 0.005 would have produced weaker results than using an initial threshold of 0.002, yet in the latter case, the results were not significantly different than when the initial threshold was 0.001 or 0.003.Further, of these four parameter sets, with an initial threshold of 0.005, the mean solution value was lower and the coefficient of variation was higher than the results using the other three sets of parameters.Table 3. Results of the application of threshold accepting to a tree-level harvest schedule of an uneven-aged oak-hickory stand in Georgia (eastern United States), based on maximizing the species mingling index and using non-unique oak tree species as neighbors.The best results when maximizing the species mingling index using unique oak species neighbors seemed to be 0.3397 (Table 4), and this solution was much more elusive to obtain as compared to maximizing the species mingling index using non-unique oak species neighbors.In only three of 3000 attempts was the solution located.One set of threshold accepting parameters seemed to produce the best results in the case where unique oak species neighbors were used to compute the species mingling index.Here, 25 feasible iterations were used at each threshold level; and the rate of change of the threshold was 0.9975, and the initial threshold was 0.001.With this set of parameters, the mean solution value (0.3302) was the highest among all of the results from using the potential sets of threshold accepting parameters, and the coefficient of variation among solutions was the lowest.However, according to the results of the two-tailed Mann-Whitney U-tests, the set of results produced using this parameter set was not significantly different (p < 0.05) than when the initial threshold was 0.002, the rate of change was 0.9975 and 25 feasible iterations were used at each threshold level.Table 4. Results of the application of threshold accepting to a tree-level harvest schedule of an uneven-aged oak-hickory stand in Georgia (eastern United States), based on maximizing the species mingling index and using unique oak tree species as neighbors.Results obtained from the generation of near-optimal sets of harvest plans for the eastern oak-hickory forest also suggest that our attempts to optimize the arrangement of harvested trees led to a greater species mingling than the random selection of trees.In the case where non-unique oak species are used to maximize the mingling index, of 3000 feasible random solutions generated, the mean mingling index was 0.3468, the maximum (best) was 0.4179 and the coefficient of variation was 6.9%.When compared to the 3000 solutions generated from all 30 threshold accepting parameter sets, only 3% of the optimized solutions were lower in value than the best randomly generated solution.In the case where unique oak species are used to maximize the species mingling index, of 3000 feasible random solutions generated, the mean species mingling index was 0.2298, the maximum (best) was 0.2863 and the coefficient of variation was 8.1%.Again, when compared to the 3,000 solutions generated from all 30 threshold accepting parameter sets, only 3% of the optimized solutions were lower in value than the best randomly generated solution.
Figure 5. Map of the harvest tree locations (circled) assuming unique tree species in the species mingling index computations, for the uneven-aged oak-hickory stand in Georgia (eastern United States).The initial threshold was 0.004; the rate of change was 0.9975; and the iterations per threshold were 25.

Western United States Mixed Conifer Forest
Results obtained from the generation of near-optimal sets of harvest plans for the western United States mixed conifer forest suggest that our attempts to optimize the arrangement of harvested trees lead to a greater species mingling index than the random selections of trees.In fact, the lowest-valued threshold accepting solution, using either non-unique or unique tree species in the species mingling index computation, was better than the very best randomly-generated solution, regardless of the parameters assumed in the threshold accepting process.
The sensitivity of the results (MDn) to the parameters is interesting.Table 5 illustrates that generally, as the rate of change in threshold accepting increases from 0.99 to 0.9975, higher value results (best and mean) are produced.This implies that the slower the threshold value is reduced during the operation of the heuristic, the better the solution produced.In essence, the heuristic search process is allowed to explore the solution space more broadly (diversify the search).Further, as the number of iterations per threshold level increase, the search explores more broadly, as well, in the initial phases and more intensively in the latter phases of the process.The best initial threshold level for this particular problem seems to be around 0.003 to 0.005 for situations when both the neighbors of the reference tree simply had to be of a different species (Table 5) and 0.001 to 0.005 when the neighbors had to be of a uniquely different species (Table 6).Thus, the sensitivity of the results with respect to the initial threshold was weaker.An example of the harvest plan developed for uniquely different species using the very best set of parameters is illustrated in Figure 6.When non-unique tree species were used in the maximization process, initiating the search with a smaller threshold (0.001) did not seem to allow the process enough of an opportunity to locate higher quality solutions (i.e., the threshold contracted too soon).Initiating the search with a larger threshold (0.005) produced high-quality solutions, but perhaps allowed the process to diversify the search too much.Rather than locate very high-quality areas of the solution space and concentrate the effort there, the search was allowed to expand beyond these with returning.The best set of solutions when computing the species mingling index with non-unique neighbors was obtained when the initial threshold was either 0.002, 0.003 or 0.005, the rate of change was 0.9975 and the iterations per threshold were 25.However, none of the sets of solutions generated using an initial threshold of 0.001, 0.002, 0.003, 0.004 or 0.005 were significantly different (p < 0.005).Similar results were found when computing the species mingling index using unique neighbors.The best set of solutions was obtained when the initial threshold was 0.001, 0.003 or 0.004, the rate of change was 0.9975 and the iterations per threshold were 25 in this case.However, the set of solutions generated using an initial threshold of 0.001 was significantly different (p < 0.005) than the other four sets of data, where the only difference in search parameters was the initial threshold (i.e., 0.002 to 0.005).Further sets of solutions with initial thresholds of 0.002 to 0.005 were not significantly different (p > 0.005).
Table 5. Results of the application of threshold accepting to a tree-level harvest schedule of an uneven-aged mixed-conifer stand in eastern Oregon (western United States), based on maximizing the species mingling index and using non-unique tree species as neighbors.Since the locations of trees in the western mixed conifer forest were randomly located, we compared the average distances between residual trees of the same species to provide some insight into the types of general rules one might employ in the field.For the randomly generated solutions to the tree-level problem, the post-harvest average distance to the first residual tree of the same species was 10.3 m for Abies lasiocarpa, 6.5 m Picea engelmannii and 5.7 m for Abies grandis.For the species mingling index that was computed using non-unique neighbors and using the five sets of 100 solutions (initial threshold = 0.001 to 0.005, rate of change 0.9975, iterations per threshold 25), the post-harvest average distance to the first residual tree of the same species was 12.8 m for Abies lasiocarpa, 8.8 m Picea engelmannii and 6.4 m for Abies grandis.Each of these distances were significantly different (p < 0.05) than the solutions generated randomly.
For the species mingling index computed using unique neighbors, the post-harvest average distance to the first residual tree of the same species was 14.0 m for Abies lasiocarpa, 7.0 m Picea engelmannii and 5.9 m for Abies grandis.Each of these distances was again significantly different (p < 0.05) from the solutions generated randomly and was also significantly different from the solutions developed for the non-unique neighbors.Table 6.Results of the application of threshold accepting to a tree-level harvest schedule of an uneven-aged mixed-conifer stand in eastern Oregon (western United States), based on maximizing the species mingling index and using unique tree species as neighbors.To further illustrate this type of guidance applied to the case study stand, the locations (x, y) of the trees were randomized to form nine other unique, alternative descriptions of the locations of the trees in the forest.These were then subjected to the threshold accepting process (initial threshold = 0.003, rate of change 0.9975, iterations per threshold 25).The average distances to the first residual tree of the same species are illustrated in Table 7 and seem to be consistent for a stand of these characteristics regardless of where the trees are located.The solution results produced by threshold accepting suggest that the distance between residual Abies lasiocarpa and Picea engelmannii should be 1-2 m longer than the random results suggest in order to provide the best species mingling index among the tree species.

Discussion
In many areas of the world, tree species diversity within forests is an important aspect of forest management and conservation efforts [3].If the interspersion of tree species is considered an important proxy of diversity across a forested landscape, then some specific operational direction may be necessary for the field personnel charged with managing the forest.The outcome of studies such as this might help serve as a basis for tree marking decisions and be used as simulated references for marteloscope training systems, where field personnel practice marking a stand of trees for harvest and receive immediate feedback on their decisions to help improve their understanding of the management of forests [33].A tree-level optimization process would complement this type of training system, which would further require growth and yield models and estimates of economic and environmental outcomes from the trees selected for harvest.One advantage of the proposed process we designed is the ability to use stand-level guidance and the spatial structure of a forest to assess the tree species present and to schedule those to be harvested, so that the interspersion of tree species remains high.What remains to be acknowledged is that when individual trees are removed, the gaps created may become regenerated with any number of native or planted species.The location and density of these are generally impossible to predict with certainty.Therefore, beyond the immediate (first time period) computation of the species mingling index (with or in absence of harvest), future development of a forest can only be modeled using assumptions of regenerated tree species, density and location, which is a venture more appropriate for gap simulators to model.
In the general definitions of the species mingling that are described in the literature, vij generally does not take into account that all of the neighboring trees (j) to a reference tree could be of the same tree species.In cases where there are several types of potential tree species, Mi may be the same value for a reference tree regardless of whether the nearest neighbors are all of the same species or all of different species.Therefore, we developed a process to accommodate both cases and illustrate the obvious: that when neighbors need to be unique tree species, MDn is lower than when this is ignored.However, the correlation between the species mingling index and more direct measurements of biodiversity, such as the presence or abundance of other species of concern [1], is still an issue that is only implied through the use of the mingling index.Further, other measures of spatial patterns of trees [17] may be of interest in tree-level optimization efforts, such as the one described here.
To account for edges and to minimize the bias in the calculation of the species mingling index (as in [34]), all trees within a certain distance of the edge of the plot may not be included in the computation of MDn; therefore, the species mingling index may be based on the following:  (10) where Te = the set of trees in a stand that are inside the assumed edge buffer.
In an earlier phase of this work and to account for edges, trees within 15 m of the edge of the plot were not included in set Te.They were, however, used to represent neighbors of other trees within the set and were also available for harvest in the optimization process.We made this assumption because in the western mixed-conifer case, every tree had at least seven neighbors that were less than 15 m away from their location, and therefore, we considered the 15-m buffer to be conservative given that some of the trees within this buffer might be scheduled for harvest.However, this assumption was met with negative review and, therefore, was ignored in the current work.We further felt that the nearest neighbor methods of [15] were not applicable to this work, because in our case, some trees will be removed (harvested) and [15] no trees are removed.Therefore, the appropriate distance to the edge of the tract cannot be determined simply based on the number of neighbors of each tree, as it is unknown a priori how many of these neighbors will remain after harvest.While this (not using a buffer when assessing the mingling index) is a limitation of the current study, a buffer zone assumption has been shown to reduce estimation error [15].
With respect to deciduous forests in North America, oak-dominated ecosystems are an important resource within the broad landscape, and concerns regarding the sustainability of the resource, prompted in part by the exclusion of fire and the employment of disturbance processes that promote the development of other tree species, have been noted [35].Silviculture decisions associated with uneven-aged deciduous forest regeneration in the southern United States can be complex and can be influenced by the level of residual basal area desired, the type of residual structure desired and the composition of residual tree species desired [36].Several other factors can influence the composition of deciduous forests managed through uneven-aged methods, including species occurrence, site quality and tree selection practices [37].For example, the likelihood of coppice sprouting among deciduous trees may vary depending on tree size prior to cutting [38], thus further complicating the silvicultural goals of oak and hickory forests, such as the one studied here.Further, the case study stand contains a large number of generally non-merchantable trees (Ostrya virginiana) that may require removal to promote the regeneration of desirable tree species (e.g., oaks).Although one goal of management may be to promote higher quality deer habitat from the current mature forest, deer exclusion may also be necessary to control browsing on oak seedlings and saplings that are regenerating within canopy gaps [39].Finally, small animal consumers of acorns produced by oak trees can complicate the desired goals associated with using silvicultural treatments that promote oak regeneration, particularly when the timing of harvest entries attempts to coincide with the periodicity of acorn production from mature trees [40].
This work extends the recent postulation of Holopainen et al. [41], who suggested that tree-level data collected via laser scanning of forest resources might be used for simulation or optimization of forest management decisions.The tree-based structure of the algorithm we described also requires information that is similar to what is used in distance-dependent growth and yield models.LiDAR technology may be useful at some point in both identifying the location of trees [42] and estimating the species type based on crown shape and stem form, although the identification of tree species may require aerial photography, as well [43].Through this information, one may be able to imply the tree species.It has been suggested that this type of model is difficult to apply in forests with more than three types of tree species.However, others have suggested that at this time, these types of methods are widely recognized as being useful for providing ecological insight into the growth and health of forests, but are not necessarily adequate for inventories of large areas [44].In the future, processes might be developed and employed to acquire the type of spatial information necessary to represent stem maps and, thus, information needed to minimize the effect of edges (e.g., [4]) when calculating structural indices, such as the species mingling index.Machine-mounted scanners on sophisticated harvesting equipment may soon be able to collect data and produce tree maps [41] that might then be employed to direct machine operators to the trees of interest for harvesting or even to calculate localized species mingling values based on residual and harvested trees within reach or along temporary trails.Unfortunately, the envisioned process would require further research on operational issues related to logging planning, as the decisions regarding felling, bucking and yarding of trees would be made in near-real time.Further, one might foresee that in the future, people charged with marking residual trees for partial harvests might employ laser technology to ensure that the distance between those of similar species is within the recommended guidelines.In the absence of the availability of specific tree locations, simulations of hypothetical forests (e.g., numerous representations of random locations of trees) may help further refine guidance for tree-level decisions.
Finally, the problem described in Section 2.3 was formulated and solved with a heuristic method because the selection of harvested trees affects the set of neighbors that are used to compute the species mingling index.Given that the harvested trees are the main decision, formulation of this problem through mixed integer non-linear programming was deemed quite difficult.Therefore, threshold accepting was chosen due to its speed and reputation for producing high-quality results in comparison with other heuristics [29].Further, the results from the heuristic were generated from random initial solutions in order to develop statistically independent samples.While the final solutions of each run of the model can only be considered near-optimal, the quality of these were assessed in a self-validation manner that required an examination of the sensitivity of results to changes in the search parameters.These processes demonstrate the application of a well-understood heuristic to a new management problem and represent an adequate level of validation [45].

Conclusions
Tree-level operational harvest scheduling optimization, using the species mingling index as the objective and recognizing other operational constraints, can be accomplished mathematically using heuristic programming methods.We have shown that near-optimal solutions generated from the use of a threshold accepting heuristic can produce harvest plans that are significantly better with respect to tree species mingling than those produced at random.However, implementation of these plans would require sophisticated databases (tree locations, species, etc.) and methods for selecting the appropriate trees in a forest in a manner that is not random.As spatial, computational and other (e.g., harvesting) technologies progress, perhaps processes like the one described here can one day be successfully accomplished by field personnel.

Figure 1 .
Figure 1.Map of the tree locations for the uneven-aged oak-hickory stand in Georgia (eastern United States).

Figure 2 .
Figure 2. Diameter distribution for the uneven-aged oak-hickory stand in Georgia (eastern United States).

Figure 3 .
Figure 3. Map of the tree locations for the uneven-aged mixed-conifer stand in eastern Oregon (western United States).

Figure 4 .
Figure 4. Generalized process flow for the threshold accepting heuristic when applied to a tree-level harvest schedule based on maximizing the species mingling index.

Figure 6 .
Figure 6.Map of the harvest tree locations (circled) assuming unique tree species in the species mingling index computation, for the uneven-aged mixed-conifer stand in eastern Oregon (western United States).The initial threshold was 0.003; the rate of change was 0.9975; and the iterations per threshold were 25.

Table 1 .
Characteristics of the case study uneven-aged oak-hickory stand in Georgia (eastern United States).

Species Stems per Hectare Average DBH d (cm) Average Mass (Metric Tons per Tree) Mass per Hectare (Metric Tons)
a Quercus coccinea Muenchh., Quercus falcata Michx., Quercus nigra L., Quercus rubra L., Quercus stellata Wangenh; b Acer rubrum L., Cornus florida L., Fagus grandifolia Ehrh., Liquidambar styraciflua L., Liriodendron tulipifera L., Oxydendrum arboreum L., Prunus serotina Ehrh., Ulmus spp.; c Pinus echinata Mill., Juniperus virginiana L.; d Diameter at breast height; e Many of these trees were too small in diameter for the volume tables to report a volume; therefore, the estimate of mass per tree is low.

Table 2 .
Characteristics of the case study uneven-aged mixed-conifer stand in eastern Oregon (western United States).

Table 7 .
Results of the application of threshold accepting to a tree-level harvest schedule of an uneven-aged mixed-conifer stand in eastern Oregon (western United States), for the original and nine other alternative representations of the case study forest.The data arose from 100 runs of the model for each alternative.