Using mixed integer programming and airborne laser scanning to generate forest management units

Adrián Pascual · Sándor F.Tóth

Abstract Airborne laser scanning (ALS) has been widely applied to estimate tree and forest attributes, but it can also drive the segmentation of forest areas.Clustering algorithms are the dominant technique in segmentation but spatial optimization using exact methods remains untested.This study presents a novel approach to segmentation based on mixed integer programming to create forest management units (FMUs).This investigation focuses on using raster information derived from ALS surveys.Two mainstream clustering algorithms were compared to the new MIP formula that simultaneously accounts for area and adjacency restrictions, FMUs size and homogeneity in terms of vegetation height.The optimal problem solution was found when using less than 150 cells, showing the problem formulation is solvable.The results for MIP were better than for the clustering algorithms; FMUs were more compact based on the intravariation of canopy height and the variability in size was lower.The MIP model allows the user to strictly control the size of FMUs, which is not possible in heuristic optimization and in the clustering algorithms tested.The definition of forest management units based on remote sensing data is an important operation and our study pioneers the use of MIP ALS-based optimal segmentation.

Keywords Segmentation · Precision forestry · Optimization · Mathematical programming · Airborne laser scanning

Abbreviations

ALS Airborne laser scanning

3D Three dimensional

CHM Canopy height model

DTM Digital terrain model

FMUs Forest management units

GEOBIA Geospatial object-based image analysis

GIS Geographic information systems

LiDAR Light detection and ranging

OF Objective function

MIP Mixed integer programming

Introduction

The integration of remote sensing (RS) data has expanded the frontiers of large-scale mapping projects (White et al.2016).Laser scanning, in particular, has become a popular tool due to the capability of Light Detection and Ranging (LiDAR) to precisely describe the three-dimensional (3D) structure of forests (Wehr and Lohr 1999; Hyyppä et al.2008).The statistics computed from LiDAR data are used to assess forest ecosystems (Næsset et al.2004; Wulder et al.2012).The outputs of LiDAR-assisted forest inventory usually include wall-to-wall estimates of tree height in raster format and canopy height models (CHM).These layers are useful to aggregate forest stands of similar characteristics (Mustonen et al.2008; Chen and Hay 2011; Pukkala 2020).The term forest management unit (FMU) refers to a homogeneous forest area of similar bio-physical characteristics that can be managed as a unit.Decision models such as harvestscheduling are formulated based on these as spatial decision units (McDill et al.2016).The definition of FMUs is very important in forest planning as they represent the calculation units on which to implement planning decisions covering long-term intervals (Pukkala 2019a).The boundaries of FMUs have been traditionally considered as fixed in time, rather than dynamic.The traditional approach of permanent FMU delineations might, in some cases, be replaced by methods that allow dynamic boundaries subject to evolving management objectives (Pascual 2018).The widespread use of LiDAR to collect input data can help with the automation of these new algorithms (Kangas et al.2018; Goodbody et al.2019).The tactical delineation of harvest units within the boundaries of more permanent FMUs would benefit, in particular, from dynamic boundary generation methods.Using ALS data can better meet those demand constraints for various forest products subject to operational constraints.

The use of high-resolution outputs from LiDAR-assisted forest inventory has induced the development of clustering algorithms to delineate homogeneous forest stands for the purpose of management (Lucieer and Stein 2005; Maxwell et al.2018; Stereńczak et al.2018).The underlying algorithmic process is also known as “segmentation” (Koch et al.2009; Räsänen et al.2013; Olofsson and Holmgren 2014).The goal of segmentation has typically been to maximize the homogeneity of forest attributes and RS statistics within the FMUs.Delineating FMUs with minimum canopy height variance is not a trivial exercise, and yet it is a very important goal of segmentation because each FMU must be manageable on the ground as a unit (Saremi et al.2014; Pascual 2018).With the integration of LiDAR-and imagery-assisted methods, researchers can boost the detection of spatial discontinuities (Pascual et al.2019; Ordway and Asner 2020).The mainstream method to conduct segmentation and generate FMUs is to use clustering algorithms under geospatial object-based image analysis (GEOBIA), which provides fast results in large-scale mapping applications (Blaschke 2010).The spatial, spectral and textural properties of the segments can be controlled under GEOBIA options, although the clustering does not allow the user to specifically control the number of segments to be produced (Blaschke 2010; Liu and Xia 2010).Some widely used algorithms, such as mean-shift segmentation (MSS), restrict the simultaneous control of the number of desired FMUs and their size (Michel et al.2015).As a result, segmentation is hardly a one-time operation, as some degree of after-math is frequently required to increase the size of segments to the targeted value for FMUs by merging or removing segments below a certain area threshold (Hay et al.2005; Koch et al.2009).Despite these limitations, GEOBIA solutions are widely implemented on geographical information systems (GIS) such as the outstanding commercial grade eCognition software (Trimble, Sunnyvale, CA, USA), for which open sourced alternative solutions exist (Gonçalves et al.2019).Clustering rules differ to the principles of spatial optimization methods oriented to find an optimal solution within the feasible region defined by spatial constraints (Rönnqvist et al.2015).The benefit of viewing segmentation as a spatial optimization problem, rather than as a series of weighted clustering rules, motivated this study.

The conceptualization of segmentation as a spatially explicit mathematical formulation has been recently addressed in forest optimization (Pukkala 2019a, b).The examples cited used heuristics to tackle complex combinatorial problems using many thousands of raster cells.Heuristics require parametrization to avoid underperforming, which is difficult to measure in the absence of a global incumbent solution acting as best bound (Bettinger and Boston 2017).However, neither the size nor the number of FMUs can be strictly controlled during the search (Pukkala 2019a, b).As a proxy, the user can balance the weights that rule the clustering to increase the average size of the resulting FMUs (Pascual 2018).Heuristics can rapidly process thousands of cells (Öhman and Eriksson 2010).One alternative to heuristics is to use methods based on exact solutions such as mixed integer programming (MIP).

The decision for using MIP over GEOBIA and heuristics might be appropriate at the operational planning level, for which the array of cells could comprise some hundreds of units (Mustonen et al.2008; Räsänen et al.2013).Rapid improvement in computing power has made exact optimization techniques such as MIP an increasingly attractive choice for segmentation tasks.This research evaluates the tradeoffs using MIP over solutions based on clustering rules or heuristic optimization as carried out in other research areas (Saǧlam et al.2006).This research presents a new MIP optimization model to create FMUs using airborne LiDAR data.The model accounts for the intra-height variation of raster cells within a given FMU while controlling the mean height difference between adjacent FMUs.The size of FMUs is also controlled, which is a major contribution when compared to clustering algorithms and heuristics.This is just one of many formulations that are possible with our proposed MIP approach.For example, one could maximize the difference in mean canopy heights between adjacent units or use other silvicultural attributes as the basis of stand homogeneity within individual units.This flexibility is one of the reasons why our proposed mathematical optimization framework for segmentation should appeal to the forestry remote sensing community.

This study expands the frontiers of MIP spatial optimization towards stand delineation.We start by formally introducing our proposed mixed integer programming model.For demonstration purposes, this model was formulated to generate a pre-defined number of management units; (1) whose areas are roughly the same (+/-20%); (2) they are as compact as possible in terms of perimeterarea ratios; and, (3) the canopy height variation within each unit is kept to a minimum.The formal, mathematical formulation of the MIP is followed by a case study using ALS data from Spain to illustrate how the model works relative to existing benchmarks such as K-means or meanshift segmentation algorithm (MSS) under GEOBIA.We found that the performance of MIP is superior to these benchmarks with regards to controlling the size and shape of management units, as well as keeping canopy height variations within the units below a predefined maximum.

Methods

Model formulation

The proposed segmentation model is cast as a mixed 0-1 integer program.The mathematical notation is introduced first to provide a convenient section of reference for all components of the model.This section is divided into four parts.Using set notation, the model parameters are first defined accompanied by an explanation as to how they were acquired and/or estimated for the purpose of this study.The definition of decision variables is next.These are the variables whose values our model seeks to optimize.These are followed by the accounting variables, created for convenience so that certain attributes of the forest that result from the various values that the decision variables take in the model can be controlled by the user.A third class of variables, termed auxiliary variables, is used to linearize cross-product terms that would render the model non-linear and otherwise computationally intractable.Although these linearization techniques are standard routines in integer programming, we explain below how they work in detail.The model itself is presented in the form of an objective function that minimizes the total perimeter of the FMUs to control their shape, along with a set of inequalities that represent the constraints the user wishes to impose on segmentation, such the number of FMUs to be generated.The mathematical domains of the variables are defined as the last part of the model formulation.Finally, a detailed explanation is provided of the objective function as well as the constraints.

Parameters

C=set of raster cells indexed byc; set C represents the grid cells of a canopy height model (CHM) derived from ALS surveys.For computational testing, we used 144 cells;

S=set of FMUs indexed bys; 5 FMUs will be generated to have an average unit size that is manageable on the ground;

E=set of adjacent cells indexed byij; ESRI ArcMap’s Polygon Neighbors tool (ESRI Inc.2017) was used to derive 506 pairwise adjacencies;

hc=height of cellccomputed as the 95th height percentile of ALS echoes ranging within cellc;

ac=area of cellc; 144 cells each 1,600 m2was used;

pc=perimeter of cellcwhich was 160 m (40 m edge);

dmax=maximum allowable deviation from the mean heigh in each FMU for all cells in the FMU; 4 m for dmaxwas used.

Decision variables

xcs=1 if raster cellcis assigned to FMUs, 0 otherwise.

Accounting variables

as=the area of FMUs computed as the sum of all raster cellscassigned to FMUs,

ps=the perimeter of FMUscomputed as the sum of perimeters of all cells assigned to FMUsminus the sum of common boundary lengths between adjacent cells times two,

dcs=the height deviation of cellc(assigned to FMUs) from the mean height of FMUs.With parameter dmax, one can set a limit on dcsand thus control the variability of canopy height within each FMU.

Auxiliary variables for linearization

Model formulation

The objective function (Eq.1) minimizes the perimeters all FMUs.The purpose of this function is to keep the shape of the FMUs compacted rather than elongated for operational convenience.Eq.2) calculates the area of FMUs based on the number, and area of cells assigned to it.It also stores the areas of the FMUs in variable as.The pair of inequalities (3) (Eq.3) constrains the areas of the FMUs to be within +/-20% of each other.In the case study that follows, the computational performance of the model was tested with different levels of area variations: +/-10%, 20%, 30% and 40%.Eq.4 calculates the total perimeter of all FMUs by taking the perimeter of all cells, each with a perimeter of 4 m, minus the total length of shared boundaries between adjacent cells that are assigned to the same FMU.The cross-product term, xcsxjsin Eq.4 needs to be linearized to keep the model computationally tractable.In order to do so, we substituted xcsxcjwith a binary indicator variablethat turns on if both cellscandjare assigned to FMUs, i.e., xcs=xjs=1 .Otherwise,=0 .This is done by introducing the following pair of inequalities (11):

Eq.5 calculates the average height of all cells assigned to each of the FMUs.The cross-product between continuous variableand the binary xcsis linearized using the following set of three inequalities (12) that enforce the substitution ωcs=.

TheMconstant was set to 50 m to act as an upper boundary on average canopy height.When xcs=1 , the last inequality reduces to, forces ωcs=≤ωcswhich, together with ωcs≤.If, on the other hand xcs=0 , the first inequality in Eq.12 forces ωcs=0 , which completes the desired substitution.

Inequalities Eq.6 compute the height difference dcsbetween hc(canopy height in cellc) and the average canopy height of FMUs.When xcs=0 , these inequalities reduce to 0 ≤dcs, which means that the height deviation of Cellcfrom the mean height of FMUsis excluded from accounting when Cellcis not assigned to FMUs.Otherwise, when xcs=1 , then≤dcs, which in turn forces |||hc--hc≤dcsand hc-|||≤dcs.This, in conjunction with Constraint Eq.7, ensures that height deviations within each FMU cannot exceed the FMU’s mean canopy height by more than dmax=4m .Constraint Eq.7 forces the maximum allowed height deviation to be less than or equal to threhold dmax, which we set to 4 m.Working in concert with the objective function (1), Constraint Eq.8 requires that all raster cellscin setCare assigned to one and only to one FMU.Lastly, constraints Eqs.9 and 10 declare the decison variables as binary and the accounting variables as positive real.

Pilot study area

A forest patch within a Mediterranean pine forest in central Spain was selected to illustrate the model (Fig.1).The main species of the mixed forest arePinus pinasterAit.andPinus sylvestrisL.Forest management in the area targets sustained yield production using area-based regulation methods with rotation periods between 100 and 120 years, and planning intervals of 10-20 years for tactical planning.The limits of the FMUs defined in the first forest plans have been preserved since manual delineation was carried out using aerial images and field observations.The use of novel and dynamic methods to upgrade and update FMUs are of the highest interests to forest managers.The selected 23.0 ha area showed variability in terms of the spatial distribution of trees, presence of gaps and the balance between young versus mature forest patches.For a broader overview of the study area, see the ground data used in Pascual et al.(2018), although in this study only airborne LiDAR data was used.The ALS surveys were collected in April 2010 using the Leica ALS60 II model.The study area was scanned from an altitude of 1200 m above ground using a Leica ALS60 II sensor.The scanning angle of 12° resulted in a 26-cm footprint.The filtering algorithm from Axelsson (2000) was used to classify ALS echoes for generating the digital terrain model (DTM) from which the canopy height model (CHM) is derived by subtracting DTM elevations from ALS echoes classified as first returns.The original resolution for both the DTM and the CHM was 1-m2cell size, which was later upscaled to 40 m for the CHM, corresponding to 144 raster cells to cover the training area.

Spatial adjacencies and problem dimension

Raster grid cells were converted into vectorial polygons to compute the ALS statistics corresponding to the echoes within the extension of each CHM raster cell (Fig.1).In ALS-based forest inventory, height-based statistics and density metrics are calculated to approximate forest attributes using estimation models (Næsset et al.2004; Magnussen et al.2010).For instance,H95is the 95th height percentile of the ALS echoes distributed within the extension of the grid cell.The close relationship betweenH95and dominant height has been widely documented (Næsset et al.2004; White et al.2016), and therefore we usedH95as the ALS statistic for the model coefficients.The formulation of the MIP optimization problem requires the relative position of all cell polygons and the adjacency relationships.The pairwise adjacencies were calculated considering two cells are adjacent when sharing an edge.

Fig.1 Overview of the laser coverage in the region for the 2010 campaign; the canopy height model (CHM, m) is mapped on a grayscale while the testing site is in color.The 40-m grid-cell map and the CHM histogram are also presented

The MIP problem was resolved four times, considering the area fluctuation between FMUs ranging from 10% to 40%.The problem instances were read and solved with the IBM ILOG CPLEX 12.9 (64 bit) on 12 threads on an i7-4790 Intel CPU (3.60 GHz) machine with 20 GB of RAM and a 64-bit Windows 10 Pro.The problem comprised a total of 3496 variables (1456 non-negative and 2040 binary) and 7141 constraints.The rules to accept optimization results were either to consider the incumbent solution to the problem (optimality gap of 0%) or the best bound achieved for a time-limit of two hours.The specified time limit is feasible for operational planning but substantially greater than what is required using GEOBIA-based clustering.The performance of MIP optimization was benchmarked against two mainstream methods relying on clustering.

Segmentation using clustering algorithms

The array of possibilities to conduct clustering includes the K-means algorithm (Jain 2010), versions of the region growing approach (Lucieer and Stein 2005), or the robust MSS (Michel et al.2015).Ease of implementation, simplicity, efficiency, and empirical success are still the main reasons for the popularity of the K-means, whose goal is to partition features to conform groups of similar feature properties.The math formulation can be found in Jain (2010): (1) the algorithm first identifies seed features used to grow each group, for which the user declares the preferred number; (2) the first seed is selected randomly but the selection of the remaining seeds is based on weighting rules that favor the selection of subsequent seeds farthest in space from the existing set of seed features.The user needs only to specify the number of clusters (five for this study).A different approach is to segment an area by controlling the textural, spatial and compactness of the resulting segments, as happens with MSS.

MSS is a non-parametric estimator framed on featurespace analysis.The gradient ascent method discovers local maximums that progressively increase to be used as nexus to create spatial clusters (Cheng 1995).When applied to imagery, feature spectral differences and their proximity drives the segmentation.Another parameter is usually necessary during segmentation: the feature minimum size, which is the minimum number of, e.g., raster cells, to compose an individual feature.MSS does not allow to explicitly control the number of FMUs as happens with the K-means algorithm because the number depends on: (1) the weights assigned to spectral and spatial components; (2) the threshold for feature minimum size; and, (3) the extension of the area to segment.The parameters of the MSS algorithm (i.e., the spectral and spatial attributes, plus the minimum size allowed for a segment) were iteratively tuned to produce five FMUs.The values for the spectral and spatial ranged from 1 to 20 as coded in the algorithm ArcMap implementation (ESRI Inc.2017).Equal weights (10) were assigned.High values for the spectral help to distinguish objects with similar spectral attributes, while high values for the spatial boost the resolution goodness of the segments.

Results

Mixed integer programming

The progress of the optimizations under the novel MIP formula were first assessed by attending to their completion.The results showed that solving the problem is feasible but the computation time was insufficient to obtain a gap in optimality below 15%.The results for the MIP problem were assessed by considering the increment in parametera; allowing area variation from 10% to 40% had no impact on the homogeneity of FMUs in terms of height but it did allow more variability in their size (Fig.2).At the same time, the maximum deviation between FMUs mean heights and the maximum height for a cell inside a given stand progressively decreased with increasing value ofa(Table 1).The spatial arrangement of cell assignments was visually displayed to assess clustering and to identify multi-part FMUs.The presence of several multipartitioned and dis-jointed FMU was observed in all scenarios tested (Fig.3).

Fig.2 Size of FMUs using a=10%

Fig.3 Spatial layout of cell assignments to the five FMUs by increasing the value of a: a 10%, b 20%, c 30%, and d 40%

Table 1 Statistics of FMUs created with the optimization problem; variance and standard deviation (Std), maximum deviation between maximum cell height and the mean height of the given FMU are presented

Benchmarking MIP solutions and clustering algorithms

The standard deviation of cell heights within stands showed the MIP as the most efficient approach to minimize height variations within FMU boundaries.The mean heights for the set of five FMUs were also compared among the three methods.The results showed the K-means as creating the most balanced mosaic of FMUs (Fig.4).High inter-FMU variability and low intra-FMU height variability are the objectives of segmentation.The results for the K-means and MSS showed great disparity on FMU mean size compared to the MIP showcase (Fig.5).The performance of both algorithms succeeded in producing single-part and compact FMUs.However, the spatial layout of the FMUs was uneven, especially for the K-means algorithm.

Fig.4 Standard deviation of cell heights within the set of five FMUs a and standard deviation of the five mean heights b; the results are for the K-means algorithm, the novel MIP problem and the meanshift segmentation (MSS)

Fig.5 Distribution of FMU size and segmentation map using the K-means algorithm (a and b), and the mean-shift segmentation algorithm (c and d)

Discussion

The goal of the new formula tackles an important step in forest management planning: the generation of forest management units on which to implement decisions.The boundaries of FMUs, sometimes called stands, are usually respected in long-term planning but when it comes to operational planning, there is a real need to provide precise units to manage forest assets by maximizing opportunities remote sensing brings.The formula developed acknowledges the intra-variation of canopy heights within the units, a major improvement compared to alternative options on heuristic optimization (Pukkala 2019a, b).The MIP formulation accounted for segmentation principles (Koch et al.2009; Räsänen et al.2013), intra-group homogeneity and inter-group heterogeneity expressed in above-ground height.It allows the user to simultaneously control the minimum and maximum size of FMUs while maximizing the goodness of the units.The clustering algorithms tested were performed poorer than the MIP formula at producing homogenous FMUs in terms of height, although clustering algorithms mostly ensured single-part units as outcomes.The goodness of clustering was substantially better for clustering algorithms compared toMIP optimization, which might perform better with more computing time.

A main outcome of this study is that the formulation works and the problem is solvable.We showed how to make use of linearization procedures to tackle relevant non-linear relationships between the equations needed in the model system.The comparison might even look better for the MIP, considering that the optimality gap between 15-20% was achieved with the computer resources used and in the time of two hours to complete the search.Arguably, the time limit was not enough to address the combinatorial problem for the more complex problems, but the aim was to compare all the MIP optimizations in the same conditions.With the increasing availability of super-computation resources, the optimality gap can decrease rapidly without delaying optimization runs beyond reasonable limits for operational forest planning.The results showed the impact, and consequences, of optimality, as FMUs after the MIP search were not spatiallycontinuous (dis-jointed).In simple terms, a given FMU could consist of more than one multi-partitioned units, if the optimization turns incomplete as shown.Once the MIP formulation is proved to be feasible, the forest manager must decide to devote more time to optimize and reduce the optimality gap, or to accept less efficient but faster solutions.

The risk of accepting solutions remote from the incumbent might not be a problem when for instance, both Model I and Model II formulation for the harvest-scheduling are followed (McDill et al.2016).There might be no penalty using better FMUs generated with the MIP rather than through heuristics or clustering algorithms, despite the units being multi-partitioned.The restrictions might come from forest legislation due to the definition of FMUs, which varies between countries (Rönnqvist et al.2015).It is noteworthy that clustering algorithms might not always avoid the multipart problem as observed in preliminary tests using MSS with a larger dataset of cells (Fig.6).

Fig.6 Spatial layout of FMUs using the mean-shift algorithm with 2304 cells a in the training site and b over a nearby forest patch of the same dimension

Neighboring settings are also an important consideration when framing MIP problems and spatial problems in general dealing with adjacencies.Are two cells sharing a corner as neighbors as when sharing an edge? In this investigation, adjacency was restricted by considering only edges.The choice could have promoted the compactness of the generated FMUs for the MIP more than for the MSS for instance.The development of the MIP formula was simplified using squared cells-the format of a raster.But if would be interesting to further explore the formulation by adding a set of constraints to make the transition boundary between stands as clear as possible.

Imagine one cellaassigned to FMUiand cellbassigned to FMUj.Ifaandbare adjacent but belong to different FMUs, it would be useful to maximize the difference in height between those cells so that the division line between the two FMUs is easy to detect and paves the way to better recognition in the field.The conceptualization of the improved problem formulation will be the core of an additional investigation using nano-segments as input units, together with raster cells from LiDAR data.The use of nano-segments is promising (Heinonen et al.2018), and thedeveloped formula could follow a first clustering to group raster cells into nano-segments of irregular shape (Pascual 2018).Moreover, the level of patchiness could be integrated in the problem formulation as suggested by Pukkala (2019b).In this study, the solutions were benchmarked based on the AP ratio, which was not, as such, the driver of the optimization process.Simple but meaningful statistics on patchiness could be tested in further development of the problem formulation (San-Miguel et al.2020).

The use of LiDAR to derive canopy height models of forest ecosystems adds value to the development of new methods to generate forest management units.The developed formulation based on the MIP only requires canopy height models (CHM) and adjacency relationships to generate a set of FMUs.The information from the CHM reflects the vertical and horizontal distribution of the forest only from the raw ALS point cloud.To segment on the basis of model predictions for forest attributes somehow bias the homogeneity by adding an error in model prediction (Pukkala 2020).Model error might be assumable for dominant heights but perhaps not when using model predictions of basal area of volume to generate FMUs (Pascual 2018).

The timing of this study is suitable in the light of recent advances in satellite missions.For instance, the NASA’s GEDI mission has released LiDAR data worldwide since 2019, including vegetation statistics of height at a spatial resolution of 25 m (Dubayah et al.2020; Duncanson et al.2020).Therefore, there is good momentum to advance further integration of laser data into operational forest planning (Pascual 2018; Stereńczak et al.2018; Packalen et al.2020).

Conclusion

The definition of forest management units is an important decision on which MIP (mixed integer programming) optimization can make a difference.The research showed how to use MIP optimization to perform the delineation of management units.The description of the problem and the linearization of non-linear constraints is a good baseline to be upscaled to larger problems.Once the MIP model has shown to be valid, the selection of MIP over clustering rules is a matter of balancing time, spatial scales, and the desired quality of the delineation.The MIP formulation can be improved when it accounts for edge effects and the integration of irregular-shape units as calculation units.

AcknowledgementsThe authors thank the University of Eastern Finland and the University of Washington for financial, technical, and scientific support during the first stages of the research.The study also benefited from the research exchange platform provided by the SuFo-Run project (Marie Sklodowska-Curie Grant Agreement No.691149).The first author (I.P.in the scope of Norma Transitória-DL57/2016/CP5151903067/CT4151900586) was supported by Fundação para a Ciência e a Tecnologia through the MODFIRE project-A multiple criteria approach to integrate wildfire behavior in forest management planning (PCIF/MOS/0217/2017).