Stand delineation based on laser scanning data and simulated annealing

The use of airborne laser scanning (LS) is increasing in forestry. Scanning can be conducted from manned aircrafts or unmanned aerial vehicles (UAV). The scanning data are often used to calculate various attributes for small raster cells. These attributes can be used to segment the forest into homogeneous areas, called segments, micro-stands, or, like in this study, stands. Delineation of stands from raster data is equal to finding the most suitable stand number for each raster cell, which is a combinatorial optimization problem. This study tested the performance of the simulated annealing (SA) metaheuristic in the delineation of stands from grids of UAV-LS attributes. The objective function included three criteria: within-stand variation of the LS attributes, stand area, and stand shape. The purpose was to create delineations that consisted of homogeneous stands with a low number of small stands and a regular and roundish stand shape. The results showed that SA is capable of producing stand delineations that meet these criteria. However, the method tended to produce delineations where the stands often consisted of disconnected parts and the stand borders were jagged. These problems were mitigated by using a mode filter on the grid of stand numbers and giving unique numbers for all disconnected parts of a stand. Three LS attributes were used in the delineation. These attributes described the canopy height, the height of the bottom of the canopy and the variation of echo intensity within 1-m2 raster cells. Besides, a texture variable that described the spatial variation of canopy height in the proximity of a 1-m2 raster cell was found to be a useful variable. Stand delineations where the average stand area was about one hectare explained more than 80% of the variation in canopy height.


Introduction
Airborne laser scanning (LS) is being used increasingly in forest inventories . The method records the locations and intensities of the return echoes of the laser pulses sent from airplanes or unmanned aerial vehicles (UAV). The original echo data are often preprocessed by calculating summary attributes for small raster cells, such as percentiles of the height distribution of the echoes (Hyyppä et al. 2008). The data can be further processed into a canopy height model, which shows the difference between the height of the canopy and the ground (Hyyppä et al. 2008).
The LS data can be used for segmentation or for predicting tree and stand attributes. The objective of segmentation is to demarcate forested areas (Wang et al. 2008;Eysn et al. 2012), segment the forest into stands or delineate individual trees (Koch et al. 2014). The so-called area-based interpretation of LS data (Nasset 2014) first divides the forest into spatial units and then estimates the biophysical properties of these units using for instance predictive models or nearestneighbor imputation (Maltamo and Packalen 2014). The spatial units may be existing or manually delineated forest stands, grid cells or automatically delineated segments. If the size of the segments is smaller than the typical size of a stand compartment, they may be called micro-stands (Maltamo and Packalen 2014). In this case, the micro-stands may be further aggregated into larger units during the later stages of the calculation and forest planning process. The target units should be large enough for implementing cuttings and other treatments. Spatial optimization is one of the methods that can be used to group the micro-stands into larger treatment units (Lockwood and More 1993;Lu and Eriksson 2000;Heinonen et al. 2007Heinonen et al. , 2018Pukkala et al. 2014).
Region merging may be the most common method of creating spatial segments from LS raster data (Mustonen et al. 2008;Wulder et al. 2008;Chen et al. 2011;Maltamo and Packalen 2014;Dechesne et al. 2017;Sanchez-Lopez et al. 2018;Pascual et al. 2019;Brauchler and Stoffels 2020;Zhao et al. 2020). The idea of the method is to merge adjacent spatial units that are similar in terms of variables that are used in the process. Baatz and Schäpe (2000) developed a segmentation software that merged small units into spatially continuous homogeneous regions. Mustonen et al. (2008) delineated stands using such a region merging procedure which took into account the shape of the segments. Koch et al. (2009) suggested a method where the grid cells are first classified based on forest type, roughness of canopy and tree height. Adjacent cells are then merged if they belong to the same class. Wu et al. (2013) developed a method that uses the mean shift algorithm to generate "raw stands." The map is then refined by using the spectral clustering algorithm. Pukkala (2021) tested self-organizing maps as a means to delineate forest stands from grids of biophysical forest attributes.
Cellular automata have also been used successfully for segmentation in several recent studies (Pukkala 2019a, b;Jia et al. 2019). Cellular automata are self-organizing algorithms based on the assumption that the interaction between cells decreases rapidly with increase in distance (Strange et al. 2001). Each cell takes one of a limited number of states, which in the case of stand delineation are stand numbers. Previous studies show that cellular automata can also be used for other purposes such as simulating the spread of pathogens and pests (Möykkynen and Pukkala 2014;Möykkynen et al. 2017), optimizing land use (Strange et al. 2001) and optimizing forest management (Heinonen and Pukkala 2007a, b).
The idea of segmenting rater data into forest into stands is to find the optimal segment number for each cell of the raster. This is a combinatorial optimization problem. Many methods have been developed for combinatorial optimization. A part of the methods, for instance, linear programming, are exact in the sense that they find the optimal solution of the problem with certainty. However, segmentation problems are nonlinear and often very large, due to the high number of cells, making it difficult to use linear programming.
Another category of methods developed for combinatorial optimization are iterative search algorithms called heuristics or metaheuristics. Examples of these algorithms are simulated annealing (Reeves 1993, Lockwood andMore 1993), genetic algorithm (Srinivas and Patnaik 1994), tabu search (Glover 1989;Boston and Bettinger 2002), great deluge , threshold accepting  and ant colony optimization (Blum 2005;Dorigo et al. 2006;Jin et al. 2016). All these methods could be used for segmentation. Simulated annealing has often been found to be the best or among the best methods when alternative heuristics have been compared in forest planning problems Pukkala and Kurttila 2005;Pukkala and Heinonen 2006;Jin et al. 2016).
Taking into account that metaheuristics are widely used in forestry and are easy to implement as computer software, there is surprisingly little research on their use in the delineation of spatial segments. Cook et al. (1996) used simulated annealing for segmenting synthetic aperture radar (SAR) images. Karasulu and Korukoglu (2011), among others, used simulated annealing to segment greyscale images. There may be no studies on the use of metaheuristics for creating spatial segments from UAV-LS or other laser scanning data.
The current study aimed to test the applicability of simulated annealing in the segmentation of LS data. The analyses used echo data collected by unmanned aerial vehicle from the Mengjiagang forest farm, located in the Heilongjiang province of China. Simulated annealing was used to create homogeneous and compact spatial segments from the UAV-LS data. The study also examined the influence of the objectives of the stand delineation on the outcome of segmentation. The purpose was to create segments that correspond to forest stands. Therefore, the term stand delineation is used in the rest of this article.

Study area
The study area is located at the Mengjiagang forest farm of Huanan county in the Heilongjiang province of China (Fig. 1). The coordinates of the center of the study area are 46°26′N, 130°43′E. The total forest area of the case study forest is 15,503 ha, of which 4438 ha are natural forests and the rest are plantations. The elevation of the case study forest ranges from 180 to 450 m above mean sea level.
The most common tree species of the plantation forests are Pinus koraiensis, Pinus sylvestris, Larix gmelinii and Picea asperata. The average canopy height is 11 m, and the tallest trees are 30 m high. The breast height diameter of the 1 3 trees is mostly 5-40 cm and the age of the planted stands ranges from 12 to 52 years. Most of the stands have been planted using a density of 3300 stems/ha. The stands have been thinned regularly so that a typical density of a nearmature stand is 500 stems/ha.

Laser scanning data
Laser scanning data were acquired for research purposes on August 12, 2019 using a RIEGEL scanner in an unmanned aerial vehicle DJI M600 Pro. The purpose of the data collecting was to test and develop methodologies for automated stand delineation, individual tree detection and tree-level forest management planning (Packalen et al. 2020). The total scanned area was 429 ha, consisting of two sub-areas (209 and 220 ha). The scanning of the 429 ha took about one hour.
The used RIEGL VUX-1UAV device is a survey-grade laser scanner with a precision of 5 mm. The system also included a GPS device (Trimble) and an inertial unit (APPLANIX). The flight speed was 10 m s −1 , and the flight altitude was 180 m above ground. The field of view of the scanner was 330°. The LiDAR sensor operated at a 380 kHz pulse repetition rate with a maximum scan speed of 200 scans/second. The system had a multiple target capability, which increased the number of return echoes. In the UAV-LS data of this study, the highest number of return echoes was five. The average pulse density was 136 points/m 2 , and the maximum density was 1792 points/m 2 ( Table 1). The standard deviation of pulse density was 91.57 points/m 2 .

Data preprocessing
A 1.5 km × 1 km rectangular subarea (150 ha) of the UAV-LS data was used to develop the method. Another 150-ha subarea was used for validation purposes. Tree species and canopy height varied within the region, with distinct stand borders in some places and more gradual changes in the other parts. Before stand delineation, the UAV-LS data were preprocessed as follows. First, noise points were removed. Second, the echoes were classified in non-ground and ground categories using cloth simulation filtering (Zhang et al. 2016). A digital elevation model (DEM) with a 1-m spatial resolution was created from a triangular irregular network (TIN) of the ground echoes. Non-ground echoes were normalized using DEM. Finally, a set of attributes were calculated for 1-m 2 raster cells using the classified and normalized echo data. These attributes included the percentiles (1%, 5%, 10%, 20%, 25%, 30%, …, 75%, …, 95%, 99%) of the height and intensity distributions of the echoes, percentiles of the cumulative heights and intensities as well as the variance, coefficient of variation, skewness, mean, maximum, minimum, median, average absolute deviation, standard deviation and kurtosis of the heights and intensities of the echoes. The LiDAR360 software (www. lidar 360. com) was used to calculate the attributes. Topographic features  such as elevation and slope were derived from the DEM. The total number of attributes calculated for the 1-m raster cells was 87.

Principal component analysis
Since the number of LS attributes available for the raster cells was high (87) it was impossible to select the ones most suitable for segmentation by trying all possible attribute combinations. It was unnecessary to use all attributes either since many of them correlated strongly with each other. Therefore, principal component analysis (PCA) was employed to seek a low number of attributes that did not correlate with each other and which would explain most of the variation in the data (Silva et al. 2016). The first three principal components explained 57.6%, 16.7% and 7.9%, respectively, of the total variation of the data (82.2% in total). The rest of the variation was explained by a high number of additional principal components, each of which explained only a small fraction of the total variation. Therefore, only the three first principal components were used for stand delineation. The original LS attribute that correlated most strongly with the principal component was selected to represent the component. The selected LS attributes are shown in Table 2.
The first principal component, PC1, described canopy height. The variable that was selected to represent PC1 was the 95% height percentile of the distribution of echo heights, abbreviated as HP_95%. It correlated closely with the 90% height percentile of the echo heights (correlation coefficient 0.994) and the heights of the 99% (0.996) and 95% percentiles (0.996) of the accumulated height of the echoes. Previous studies show that these LS metrics correlate closely with the dominant height of the stand. For example, Guerra-Hernández et al. (2021) could explain 89.14% of the variation in dominant height by the 99% percentile of the distribution of the UAV-LS echo heights when the laser pulse density was 0.5-1 pulses/m 2 . Higher pulse densities can be expected to strengthen the relationship since the likelihood of hitting the treetops increases.
The representative variable of the second principal component was the height of the 5% percentile of the accumulated height of the echoes (AH_5%). It correlated closely for instance with the height of the 10% (0.971), 20% (0.946) and 25% (0.937) percentile of the accumulated height of the echoes. The variance of the echo intensity (IntVar) was selected to represent the third principal component. IntVar correlated closely with the standard deviation of the intensities (0.976), average absolute deviation of intensities (0.969) and the difference between 75 and 25% percentiles of the intensities (0.854).
Based on Jia et al. (2019), the original grids of the three selected attributes in the 1-m resolution were low-pass filtered using a 5 × 5 window. At the same time, the original grids were resampled to a 5-m resolution. After resampling and low-pass filtering, the size of the original grid of 1500 × 1000 cells was reduced to 300 × 200 cells.
Using 5-m cells instead of the original 1-m cells has been found to improve the delineation result since it decreases the variation of LS attributes within distances shorter than the diameter of tree crowns (Jia et al. 2019). The use of lowpass filtering and a larger cell size also contributes to more straight stand boundaries. The use of the original 1-m cells might result in delineations where stand borders would follow the perimeters of the crowns of large individual trees (Jia et al. 2019), producing very rugged stand borders. Jia et al. (2019) concluded that, when stand delineation is based on a canopy height model, the optimal cell sizes are 5 to 9 m. As an alternative, the obtained stand boundaries could be smoothed afterward using polygon regularization of some other method (e.g., Xie et al. 2017).

The simulated annealing algorithm
In this study, the simulated annealing (SA) metaheuristic was used to select the optimal stand number for each cell of a grid (Reeves 1993). The SA algorithm developed for stand delineation works as follows: 1. Set the values of the parameters of the algorithm: starting temperature (T ini ), freezing temperature (T end ), cooling multiplier (C), number of iterations per temperature (I) and the area of initial stands (A) 2. Generate an initial solution, i.e., give a starting stand number for each cell 3. Calculate the objective function value for the initial solution (OF Current ) 4. 4. Set the current temperature (T) equal to the starting temperature (T ← T ini ) 5. 5. Repeat for I iterations: • Generate a candidate move: • Select a random cell from the grid • Do nothing if all adjacent cells have the same stand number as the random cell • Otherwise, make a move: select randomly one of the stand numbers of adjacent cells and give it to the cell • Calculate the objective function value for the move (OF Candidate ) • Accept the move if it improves the objective function value (OF Candidate > OF Current ) • Otherwise, accept the move with the following probability (p): (1) p = exp OF Candidate − OF Current T 6. Multiply temperature by C (T ← C × T) 7. Repeat steps 5 and 6 until the current temperature (T) falls below the freezing temperature The algorithm represents the standard SA algorithm in all other respects except that only those cells that are located at stand borders can participate in a move. The functioning of the algorithm is visualized in Fig. 2 where a randomly selected cell is marked with blue color. In the left-hand-side diagram, all adjacent cells have the same stand number as the blue cell, and no candidate move is produced. In the middle diagram of Fig. 2, there is one candidate move where the number of stand 2 is given to the blue cell. In the right-hand-side diagram, the candidate is produced by selecting randomly the number of stand 1, 2 or 3 for the blue cell.
The initial solution was generated by dividing the grid into square-shaped 2-ha initial stands. The number of final stands may differ from the number of initial stands since a stand can disappear during the SA run or become divided into two or more stands. Two hectares were selected after preliminary runs with 1-ha, 2-ha and 3-ha initial stands.

Objective function
The objective function that was used to evaluate the quality of the delineation depended on three criteria: The stand number of the blue cell is not changed Select stand number 1 or 2 for the blue cell Select stand number 1, 2, 3 or 4 for the blue cell  Fig. 2 The principle of producing candidate moves in the simulated annealing algorithm developed for stand delineation. The current stand number of the randomly selected blue cell is the same as the stand number of the colored large rectangle within which the cell is located. On the left, no candidate moves are produced because all adjacent cells have the same stand number as the blue cell. In the middle, the only candidate move is to assign the ID number of stand 2 to the blue cell. In the right, the candidate move consists of selecting randomly the ID number of stand 1, 2 or 3 to the blue cell where w i is the weight and p i is the sub-priority function of criterion i. The criteria were stand area (Area), withinstand variance of stand attributes (Variance) and shape of the stand (Shape). Each criterion had an associated sub-priority function, which scaled the range of the original criterion between zero and one, expressed whether the criterion was maximized or minimized, and indicated how the criterion variable contributed to the priority (Fig. 3). The sub-priority functions were subjective since the acceptable level of within-stand variation, irregularity of stand shape and allowed range of stand area are all subjective requirements set for the stand delineation (see Baatz andSchäpe 2000, Haywood andStone 2010). The sub-priority for stand area was where Area is the area of the stand (ha) and a 1 and a 2 are parameters. The effect of the within-stand variance of stand attributes was calculated from where RelVar is the weighted mean of the relative variances of the stand attributes considered in the delineation and b 1 and b 2 are parameters. The relative variance of an attribute (RV) was calculated by dividing the variance of the cell values by the mean value of the attribute within the stand. This study used three ALS attributes ( Table 2). The weighted mean of their relative variances was calculated from where v i is the weight of attribute i. The weights of the attributes were proportional to the percentage of the total variation explained by the corresponding principal component. (3) The percentages of 57.6 (HP_95%), 16.7 (AH_5%) and 7.9 (IntVar) shown in Table 2 were converted to weights v 1 = 0.7001, v 2 = 0.2032 and v 3 = 0.0961, respectively.
The priority function for stand shape was where n is the number of cells that belong to the stand, RelDist is the relative distance of the cell from stand centroid, and c 1 and c 2 are parameters. RelDist was calculated by dividing the distance by the radius of a circle that had the same area as the stand. Therefore, RelDist is the distance of the cell from the stand centroid, expressed as the number of radii of a circular shape. Equation 6 calculates the standlevel shape index as the mean of cell-level distance scores. Sub-priority function p 3 is maximized when the shape of the stand is circular. The objective function value (Eq. 2) and the three sub-priority functions (Eqs. 3, 4, and 6) were calculated separately for each stand. Since a move of the SA algorithm changed the stand number of a single cell, it affected the objective function values of two stands. Therefore, OF Current in Eq. 1 was equal to the average objective function value of the two stands before implementing the move and OF Candidate was the average objective function value of the same two stands after implementing the move. The move was considered an improvement when the average objective function value of the two affected stands improved. If the move decreased the average objective function value, it was accepted with a probability calculated from Eq. 1.

Stand delineation cases
The performance of the SA algorithm was tested with different weights of the criteria of the objective function (Eq. 2). A low within-stand variance of the ALS attributes runs of this study. The vertical line of the right-hand-side diagram indicates the radius of a circle that has the same area as the stand was assumed to be the main objective of the delineation. Therefore, a delineation where the within-stand variance was minimized without other objectives was referred to as the baseline delineation (Table 3). Then, the effect of increasing weight of stand area or stand shape was analyzed. The parameters of the SA algorithm were as follows: initial temperature (T ini ) 0.1, freezing temperature (T end ) 0.00001, cooling multiplier (C) 0.95, number of candidates per temperature (I) 50,000. Randomly selected cells that were not located at the stand border (Fig. 2, left) were also counted as candidates (they were included in the count of 50,000 candidates per temperature) although their stand number was never changed. The SA parameter values were based on recommendations given in previous studies (Pukkala and Heinonen 2006;Jin et al. 2016) and preliminary optimizations where three different values of each parameter were tested. The preliminary optimizations showed that the delineation result was rather insensitive to the parameters if the cooling multiplier was 0.9 or higher, the number of candidates per temperature was at least 10,000, and the freezing temperature was lower than 0.0001.

Evaluation of stand delineations
The stand delineation results were numerically evaluated by calculating, for each of the three LS attributes (HP_95%, AH_5%, IntVar), the proportion of total variance explained by the delineation (R 2 ). The formula for R 2 was where SSE is the variation not explained by the delineation and SST is the total variation of the attribute within the grid. SST and SSE were calculated from the following formulas: where N is the number of stands, n j is the number of cells in stand j, y ij is the value of the attribute in cell i of stand j, y is the overall mean of the attribute, and y j is the mean value of the attribute among the cells that belong to stand j. Also, the average variation, area and shape scores of stands were calculated for each delineation, using Eqs. 3, 4 and 6. The average score was calculated using stand area as a weight variable.
It was seen that the SA algorithm often resulted in "gradual" stand borders, and the stands did not always form continuous areas. Because this kind of SA results would be hardly acceptable for forestry practice the delineation results were fine-tuned in two ways. First, the stand numbers assigned to the cells in the SA run were mode-filtered by replacing the original stand number with the most common stand number within a 5 × 5 window centered on the cell. Figure 4 shows that mode filtering greatly reduced the (7) ruggedness of the stand borders and the intermingling of stand numbers. It removed most of the isolated cells or small groups of cells that were located within adjacent stands. The other fine-tuning step was renumbering where the non-connected parts of a stand were given unique stand numbers. Figure 4 shows that mode filtering slightly decreased the R 2 of the three LS attributes used in delineation. As a consequence, the variation score decreased since the within-stand variation of the LS attributes increased. However, the average area and shape scores improved. Renumbering increased the R 2 of all attributes, which is a logical consequence of dividing disconnected parts of stands into separate new stands. Renumbering decreased the area score but increased variation and shape scores. All the numerical results reported below were calculated for modefiltered and renumbered stand delineations.

Effect of delineation criteria
Increasing the weight of the area criterion of Eq. 2 slightly decreased the degree of explained variance of the UAV-LS attributes (Fig. 5, top). However, the effect was smaller than the effect of an increasing weight of the shape criterion (Fig. 5, middle) or the effect of a simultaneous increase in the weights of area and shape (Fig. 5, bottom). In all delineations, the R 2 was the highest for HP_95% and lowest for IntVar, which is most probably a consequence of the weights of the UAV-LS attributes in Eq. 5.
Since an increasing weight of stand area had no clear effect on R 2 , it is logical that it had no clear effect on the average variation score of stands either (Fig. 6, top). A seemingly illogical result was that increasing the weight of the stand area did not systematically increase the mean area score of stands. This outcome is explained by the fact that the scores were calculated after renumbering where the non-connected parts of stands were given unique stand numbers. Increasing the weight of the stand area had no major effect on the average shape score of stands.
The effect of an increasing weight of the shape criterion was strong (Fig. 6, middle). Both the area and shape scores improved when the weight of the shape criterion was increased. The mean area-weighted variation score of the stands decreased slightly. The mean area score of stands also increased greatly when the weights of the area and shape criteria were increased simultaneously (Fig. 6, bottom). The shape score also increased but slower than the area score.
Increasing the weight of the shape score increased the average stand area and decreased the proportion of stands smaller than 0.1 ha (Table 3). This happened also when the weight of the shape criterion was increased simultaneously with the weight of the area. Increasing only the weight of the area did not have a similar effect, most probably because  Table 3 for the criteria weights   Table 3 suggest that using stand shape as a criterion led to delineations where fewer stands consisted of several disconnected parts before the renumbering step. The strong effect of the shape criterion can also be seen from Fig. 7 where different stands are shown with different colors. The baseline run, which minimized within-stand variation of the ALS attributes without any other criteria, led to a delineation with many small stands. Besides, the shape of several stands was irregular. Increasing the weight of the shape criterion, either alone or simultaneously with the weight of the stand area, led to more roundish stands and a lower number of small and irregular stands.
Based on the results presented above, it was concluded that the shape criterion should have a weight higher than zero in the objective function (Eq. 2). Increasing the weight of the area criterion also improved the delineation, but only when the shape criterion had a weight higher than zero. For example, the weight combination of Multi 3 (Table 3), where the weight of both area and shape was 0.15, resulted in good delineations (Fig. 8). Clear improvements in stand shape and area were achieved without large losses in within-stand homogeneity. Figure 8 shows that the delineation closely followed the changes in the three UAV-LS attributes. There were a few narrow strips of cells where the attribute value differed from the average value of the stand. These places were mostly narrow corridors that had been planted later than or with a different species as the adjoining areas. Also, a few very small stands remained in the delineation.
The effect of the initial stand area used in simulated annealing is shown in Fig. 9 (top). Using 1-or 2-ha squares as initial stands led to better delineations than the use of 3-ha initial stands. The average area on the stands was the Fig. 6 Effect of increasing weight of stand area and stand shape criteria in Eq. 2 on the mean variation, area and shape scores of stands. The mean scores are calculated by using the stand area as a weight variable. See Table 3  largest when 2-ha initial stands were used. Three-hectare initial stands led to delineations where the stand consisted of several disconnected parts, which were divided into separate stands in the renumbering step. The delineation result and the average variation, area and shape scores were not very sensitive to the weights of the three ALS metrics used in Eq. 5 (Fig. 9, middle). However, when all weight was on HP_95%, the area score of the delineation was worse than in cases where the segmentation was based on all three ALS metrics.
The boundaries delineated by the simulating annealing method coincided with manually drawn boundaries but, Shape 5 M ulti 5 Fig. 7 Final stands after the simulated annealing run, mode filtering and renumbering in a few stand delineation cases (see Table 3 for the explanation of the cases) compared with manual delineation, simulated annealing produced more detailed delineations (Fig. 9, bottom). The simulated annealing method produced more stands and stand boundaries than obtained in manual stand delineation. Especially the delineation of the current forest management map (blue lines in the lowest map of Fig. 9), which was drawn on aerial photograph in 2016, contained fewer and larger stands than obtained in simulated annealing.

Delineation results for an independent test area
The method developed for stand delineation was applied to another area in the same forest where the stand borders were less distinct. The delineation was done using the same criteria weights as in Multi 3 (0.15 for both area and shape, see Table 3). The SA parameters and the sub-priority functions were also the same as before.
Compared to the case study area, the delineation of the independent test area included a lower number of larger stands (Fig. 10), due mainly to the fact that the independent test region consisted of large continuous areas of similar homogeneous forest. The stand shape was usually rather circular, and stand boundaries were drawn in places where HP_95% (≈ canopy height) changed instantly. However, there were a few places without a stand boundary where a human interpreter most probably would have drawn a stand border. These places are marked with red ellipses in Fig. 10. In these places, the canopy height (HP_95%) did not change much but the stand density and crown size of the trees changed.
As a potential mitigation of the problem, an additional texture variable was calculated for the cells (Fig. 10, bottom) and used as a fourth ALS attribute in Eq. 5. The texture variable was calculated as follows. First, the variance of HP_95% within a 5 × 5 window centered on the cell was calculated for each cell of the original HP_95% grid (1-m resolution). The variance described the spatial variation of HP_95% in the surroundings of the cell. Then, the variances calculated for the 1-m 2 cells were low-pass filtered and resampled using a 5 × 5 window in the same way as was done the other attributes. The weight of the additional texture variable was 0.25 and the weights of the other variables were multiplied by 0.75 to keep the total weight of the attributes equal to one: The results obtained with the additional texture variable removed the problem of missing stand boundaries in places where the texture of the canopy changed (Fig. 10, bottom). The stand marked with the red line is an example of a fine-textured stand delineated differently than in the case where Eq. 5 was used (Fig. 10, top). Figure 10 (bottom) shows that there were stand boundaries in places where the texture of the HP_95% attribute changed.

Discussion
The study showed that simulated annealing can be used to delineate stands from UAV-LS data. The use of the method is not restricted to LS data, and the LS variables can be different from those used in this study. It seems clear that (10) RelVar = 0.526RV HP_95% + 0.152RV AH_5% + 0.072RV IntVar + 0.25RV Texture Fig. 8 The Multi 3 delineation overlaid with HP_95% (top), AH_5% (middle) and IntVar (bottom). A lighter tone indicates a higher value of the attribute. The weights of the criteria of the objective function (Eq. 2) were: within-stand variation 0.7, stand area 0.15 and stand shape 0.15 other metaheuristics developed for combinatorial optimization such as tabu search (Reeves 1993) and genetic algorithm could also be used. Additional metaheuristics that should also be tested in stand delineation include great deluge , threshold accepting  and ant colony optimization (Dorigo et al. 2006;Jin et al. 2016).
The SA algorithm used in this study corresponds to the basic SA algorithm with the difference that the stand numbers of cells that were not located at stand boundaries were

Comparison with manual delineation
Blue: Current sub-compartment map Red: Forester on canopy height model Yellow: Simulated annealing Fig. 9 Effect of initial stand area (top) and the weights of the ALS metrics (middle) on the stand delineation result when the weights of the criteria of the objective function (Eq. 2) were: within-stand variation 0.7, stand area 0.15, and stand shape 0.15. The lowest map shows the sub-compartment borders of the current forest map and a manual delineation of a forester. The 95% percentile of the distribution of echo heights is shown in black and white never changed. The algorithm is easy to use and program, and SA pseudocodes can be found on the internet. Taking into account the wide use of SA and other metaheuristics and the fact that assigning optimal stand numbers for grid cells is a combinatorial problem, it is surprising that these methods have not been used more in stand delineation. It might be expected that they should be the first option, offering benchmark delineations for the other methods.
The current study showed that SA did not always produce smooth stand borders. Sometimes the transitions of the stand numbers of cells between two stands were gradual (Fig. 4, top), which may correspond to reality in many forest ecosystems. The restriction that stand numbers may change only in cells that have adjacent cells with a different stand number did not prevent the algorithm from splitting stands to disconnected parts, many of which were only one-cell-wide. However, the problem was easily mitigated by applying a mode filter on the stand number grid obtained from the SA. The mode filter rectified and smoothed the stand borders with only small losses in the homogeneity of the LS attributes within the stands. Another post-optimization step was renumbering where unique stand numbers were given to the non-connected parts of the stands. Additional enhancement was used in the independent test area, where a spatial texture variable was used as an additional UAV-LS attribute.
It would be straightforward to test several other possibilities to improve the SA-based delineation, for instance using the renumbered solution as the starting solution for a new SA run. The SA solution (non-filtered, not renumbered) could also be used as the starting solution for a cellular automaton (Jia et al. 2019), which would remove all one-cell-wide stand parts and lead to less rugged stand borders. Another option to further improve the method would be to use the common border between a cell and the adjoining stands as a fourth criterion in the objective function (Eq. 2), as was done in the cellular automaton applications proposed for stand delineation (Pukkala 2019a(Pukkala , 2019bJia et al. 2019). This might also lead to smoother stand boundaries.
Additional possibilities to improve the delineation would be to use additional layers obtained from sources other than laser scanning. For example, a species layer would most probably be available for most forests, at least for plantations (Jia et al. 2019). Digital aerial photographs may be used to calculate vegetation indices (for instance NDVI) and other variables for the raster cells (Xu et al. 2015). A layer separating the forest from non-forest (roads, agricultural land, etc.) could also be created from base maps and used to mask off cells that do not represent forestland (Wang et al. 2008;Eysn et al. 2012;Pukkala 2019b).
The use of laser scanning data in forest inventory and forest management planning is still a rather new field. However, the popularity of LS-based forest inventories is increasing rapidly , which means that alternative ways to process and utilize the LS echo data should be developed and tested. This study aimed at creating homogeneous spatial units from UAV-LS data, to be used in the subsequent stages of forest inventory and planning calculations. The target size of the spatial units was rather large because we aimed at spatial units that corresponded to traditional stand compartments, which are large enough for the implementation of management actions.
Another possibility would be to create smaller segments, which are more homogeneous and may offer advantages in the prediction of growing stock attributes for the forest (Pascual et al. 2019). In this case, the segments should be further aggregated into larger treatment blocks. This could be done during the later stages of the planning process, using optimization (Heinonen et al. 2018). It is also possible to delineate segments that correspond to the crowns of individual trees (Vauhkonen et al. 2009;Koch et al. 2014;Packalen et al. 2020). In this case, spatial optimization may be used Fig. 10 Overlay of the stand delineation (Multi 3) and HP_95% in an independent test area (top). The red ellipses show places where canopy texture changes within a stand. The lower map shows the delineation where a texture variable was used as an additional ALS attribute. In the lower map, the stand boundaries are overlaid with the variance of HP_95% within a 5 × 5 window. Light tone indicates a high variance of HP_95% to aggregate the individually detected trees into stands or harvest blocks (Packalen et al. 2020).
The selection of the UAV-LS attributes that would lead to the best delineation results was not systematically examined in this article. Since the number of the attributes is high, the evaluation of all their possible combinations is impossible. Combinatorial optimization might be used in future studies to automatize and fasten the process of selecting the best combination of attributes for stand delineation. This study used principal component analysis to find a few features that would explain most of the variation among all UAV-LS attributes. One attribute was selected to represent each principal component. Although it was not shown that the three attributes used in our study are the optimal combination, we are confident that no drastic improvements could be obtained by using other attribute combinations.
It seems evident that an LS attribute that describes canopy height should be used in stand delineation. The other attributes used in our study described the height of the crown base and variation of echo intensity within a 1-m 2 cell. Analyses conducted for a test area suggested that variables that describe the variation on a different scale could improve the delineation result. We used the variance of canopy height between 1-m 2 cells within a 5 × 5 window to describe the texture of the 95% percentile grid. It helped to distinguish adjacent stands with equal canopy height but different textures. (Technology Development Team for High-efficient Silviculture of Forest Resources).
Funding Open access funding provided by University of Eastern Finland (UEF) including Kuopio University Hospital.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.