Next Article in Journal
Ecological Strategy Spectra for Communities of Different Successional Stages in the Tropical Lowland Rainforest of Hainan Island
Previous Article in Journal
Analysis of the NAC Gene Family in Salix and the Identification of SpsNAC005 Gene Contributing to Salt and Drought Tolerance
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Optimal Survey Design for Forest Carbon Monitoring in Remote Regions Using Multi-Objective Mathematical Programming

by
Sándor F. Tóth
1,
Kiva L. Oken
2,†,
Christine C. Stawitz
2,‡ and
Hans-Erik Andersen
3,*
1
School of Environmental & Forest Sciences, University of Washington, Seattle, WA 98195, USA
2
Quantitative Ecology and Resource Management, University of Washington, Seattle, WA 98195, USA
3
Pacific Northwest Research Station, USDA Forest Service, Seattle, WA 98195, USA
*
Author to whom correspondence should be addressed.
Present Address: Fishery Resource Analysis & Monitoring Division, Northwest Fisheries Science Center, National Marine Fisheries Service, NOAA, Seattle, WA 98112, USA.
Present Address: NOAA Fisheries Office of Science and Technology, Seattle, WA 98115, USA.
Forests 2022, 13(7), 972; https://doi.org/10.3390/f13070972
Submission received: 27 April 2022 / Revised: 3 June 2022 / Accepted: 7 June 2022 / Published: 22 June 2022
(This article belongs to the Section Forest Inventory, Modeling and Remote Sensing)

Abstract

:
Cost-effective monitoring of forest carbon resources is critical to the development of national policies and enforcement of international agreements aimed at reducing carbon emissions and mitigating the impacts of climate change. While carbon monitoring systems are often based on national forest inventories (NFI) utilizing a large sample of field plots, in remote regions the lack of transportation infrastructure often requires heavier reliance on remote sensing technologies, such as airborne lidar. The challenge motivating our research is that the efficacy of estimating carbon with lidar varies across the various carbon pools within forest ecosystems. Lidar measurements are typically highly correlated with aboveground tree carbon but are less strongly correlated with other carbon pools, such as down woody materials (DWM) and soil. Field measurements are essential to both (1) estimate soil and DWM carbon directly and (2) develop regression models to estimate tree carbon indirectly using lidar. With limited budgets and time, however, decision makers must find an optimal way to combine field measurements with lidar to minimize standard errors in carbon estimates for the various pools. We introduce a multi-objective binary programming formulation that quantifies the tradeoffs behind the competing objectives of minimizing standard errors for tree carbon, DWM carbon, and soil carbon. Using NFI and airborne lidar data from a remote boreal forest region of interior Alaska, we demonstrate the operational feasibility of the method and suggest that it is generalizable to other carbon sampling projects because of its generic mathematical structure.

1. Introduction

Forests sequester carbon from the atmosphere as part of photosynthesis, and carbon is stored within several ecosystem pools: aboveground, belowground, deadwood, litter, and soil. In the face of climate change that is triggered by increasing concentrations of atmospheric carbon, the role of forests in storing and emitting carbon is critical in the global carbon cycle. While changes in climate are projected to be more pronounced in higher latitudes [1], the carbon cycling process in boreal forests is poorly understood, even though these ecosystems are estimated to contain more carbon than any other terrestrial ecosystem on earth [2]. Forest biomass and soil carbon down to a depth of 1 m in boreal forests was estimated to be 559 gigatons by Watson et al. [3]. In interior Alaska, an estimated 1.3 billion tons of carbon are sequestered within live trees alone [4]. Soils are estimated to contain even more carbon, mostly in permafrost layers that have been frozen for many years and remained mainly undisturbed by human activities [5]. However, forests in Alaska are subject to various disturbance processes, such as fire, insect activity, and permafrost degradation, that can significantly affect carbon storage and emissions within this boreal region. Climate forecasts predict that air temperatures will rise more rapidly at higher latitudes [6] and rising annual temperatures across Alaska during the past 60 years confirm that this warming is already happening.
Due to the complex nature of the interactions between a changing climate, carbon cycle, fire regimes, and human disturbance, it is difficult to predict how carbon stocks in Alaskan boreal forests will change in the future. Monitoring these forests for changes in carbon stocks is crucial in understanding how, when, and at what pace climate and human changes will affect this ecosystem. In turn, these ecosystem changes will likely affect rural Alaskans and have a feedback effect on the velocity of climate change itself.
The Forest Inventory and Analysis (FIA) program, which is the national forest inventory of the United States conducted by the USDA Forest Service, aims to monitor the status and trends of all public and private forests in the nation via annual assessments [7]. The program collects data to estimate the weight of forest carbon pools, among other forest attributes, according to a nationally standardized field sampling design and nationally consistent plot and subplot sizes [7]. However, as of 2020, Alaska is the only state in the United States where the annual FIA system has not been fully implemented. Although it is estimated that Alaska contains 17% of all forest area in the United States [8], most of its forest lands have not been inventoried in a systematic manner that is consistent with the national FIA program [4,9]. At the same time, the lack of infrastructure in this region makes the application of the standardized national field sampling design cost-prohibitive. Remote sensing techniques, such as airborne lidar (light detection and ranging) and digital imaging, may allow for more cost-effective estimation of carbon across large areas. However, lidar techniques can only be used to accurately estimate some, but not all, of the carbon pools that comprise total forest carbon. For example, aboveground tree carbon can be estimated with lidar quite well (the correlation coefficient for tree carbon based on observed pilot data in our study was 0.91), but field surveys are more appropriate for soil carbon pools. In other words, estimation of total forest carbon across the multiple pools must entail a combination of lidar methods, such as lidar, aerial or satellite imaging, and ground sampling techniques. Because different sampling designs have varying projected standard errors and costs, finding the optimal combination of ground vs. lidar protocols subject to budgetary, logistical, and quality constraints is a non-trivial task. In many remote regions, logistical issues might be the binding constraint on sampling. Internationally, developing nations with significant forest inventories but limited budgets and infrastructures wishing to adopt REDD+ (reducing emissions from deforestation and forest degradation, plus the sustainable management of forests and the conservation and enhancement of forest carbon stocks) regimes in exchange for financial compensation will also have to face the same constraints. The REDD+ program requires national sampling, reporting and certification protocols for forest carbon stocks and emissions that are transparent and verifiable [10]. While the requirements of REDD+ programs are often a primary driver for the implementation of forest carbon monitoring systems, there are other national-scale reporting activities that depend upon forest monitoring, including the FAO global forest resources assessment [11] and nationally determined contributions (NDC) submissions to reduce national emissions and adapt to the impacts of climate change as required by the Paris Agreement [12].
The application of mathematical programming can both (1) inform the design of new carbon sampling strategies in remote areas with no monitoring systems (i.e., REDD+), and (2) provide a quantitative basis for evaluating tradeoffs between sampling alternatives in existing inventory programs (i.e., the FIA program in the United States). However, while previous studies have evaluated sampling alternatives for forest inventory through a straightforward comparison of standard errors for various combinations of lidar and field plots [13], to date very little work has been done to document the tradeoffs behind the survey design choices and the associated costs using more advanced mathematical optimization techniques. One notable exception was the work of Köhl et al. [14], who compared the projected percent standard errors of five sampling design alternatives in two different cost scenarios for remotely sensed data acquisitions in Puerto Rico. While the authors’ tradeoff analysis did provide future users with a rough strategic guide based on the expected precision of the five design alternatives, it failed to take actual logistical and budgetary constraints into account. For example, available budgets and skilled field labor might change, and some field plots might simply be inaccessible, even in developed countries such as the United States (Alaska). These problems could be even more pronounced in developing nations aspiring to be part of REDD+.
In this study, we present a multi-objective mathematical programming model that can help analysts identify the set of so-called Pareto optimal combinations of field and lidar sampling strategies for boreal forest carbon inventory. We seek those alternative combinations of field and lidar plot intensities that would lead to the smallest possible expected standard error in the estimated mean carbon tonnage for the different pools. A field-lidar plot intensity combination is “Pareto optimal” (i.e., efficient or non-dominated) if no other intensity combination would yield a smaller expected standard error for one carbon pool without increasing the standard errors for any one of the other pools. Multi-objective mathematical programs are sets of functions that represent the objectives of the analysis, such as minimizing standard errors, and inequalities that represent the constraints, such as budgetary restrictions, within which the objectives can be pursued. Since lidar and field sampling methods have different associated costs and are projected to lead to different standard errors for the different carbon pools, our objective is to minimize standard errors. The constraints are finite budgets of both money and time to be spent on ground plots, as well as time to be spent on processing remotely sensed data. Using data collected from an FIA pilot project in the Tanana Valley of interior Alaska (USA), we demonstrate how multi-objective mathematical programming can help identify the best compromise sampling intensity combinations.

2. Materials and Methods

For simplicity of illustration, we combined the five carbon pools, mentioned previously, into three that can be measured directly: tree carbon, down woody material (DWM) carbon, and soil carbon. The first pool, “tree”, included both live and dead standing trees; DWM included both fine (<3 inch (7.62 cm) diameter) and coarse (≥3 inch (7.62 cm) diameter) down woody materials; and “soil” included litter. We chose not to include belowground biomass in our analysis because it is typically predicted based on aboveground tree biomass. Our case study site for the estimation of the three carbon pools listed above was a 13.7 million ha forest area in the Tanana valley in Alaska (Figure 1). Forty percent of this area has been previously inventoried [4]. The region is one of the most sparsely populated in the United States and has very few roads and only one small city. As such, measuring nominal 5642 ground plots per standard FIA program specifications is not economically viable (standard FIA plot intensity is 1 plot per 6000 acres (2428 ha)). Thus, there is interest in using airborne lidar technologies to monitor the area in conjunction with a limited number of field plots.
The rest of this Materials and Methods section is structured as follows. We start by providing a mathematical description of the proposed multi-objective linear programming model. We then parameterize the model for our case study using parameter estimates derived exclusively from observed data (provided by a pilot study conducted within the area of interest). The primary objective of this exercise is to give the reader a sense of how this model works and what benefits it can offer. Third, we use a computer algorithm, implemented in R, to iterate through all possible solutions to the problem and identify those that are Pareto optimal. Last, we perform a sensitivity analysis to determine how changes in model parameter values can affect the solution set. The goal of the sensitivity analyses is to demonstrate the utility of the proposed model in situations where parameter estimates based on observed data are not available.

2.1. The Multi-Objective Linear Programming Model

To formulate this problem as a linear program, we first enumerated the feasible combinations of field and lidar plot intensities, where intensity refers to the number of plots of each type to sample in the region. We use the term “lidar plot” to describe a location where only remotely sensed metrics (i.e., no field measurements) are available. In our model, we use a combination of lidar and field plots in a so-called “double sampling” inventory design, where remotely sensed metrics (lidar-based canopy heights, etc.) are collected on a large sample of relatively inexpensive “lidar plots” within a region of interest, and then expensive field measurements (i.e., direct measurements of forest carbon) are collected on a smaller subsample of these plots (“field plots”) [15]. For field plots, the maximum intensity was set to be approximately the maximum number of field plots that could fit within the given budget: 250 plots. For lidar plots, the maximum intensity was set to 5480, as requested by the decision makers. Next, the minimum intensity was set at a small number relative to the maximum intensity; we chose 10 for field plots and 180 for lidar plots. Since the objective is to minimize standard error, the model formulation will seek to maximize the number of plots subject to the constraints. Thus, the solutions the model formulation returns should not be sensitive to the chosen minimum intensity values. The decision variables, denoted by x i j   , are therefore defined as binary, with one decision variable for every combination of field plot and lidar plot intensity that we considered, as represented by subscripts i and j, respectively. If the intensities evaluated included every integer value between the minimum and maximum intensity values, there would be more than one million decision variables. Since our goal was to illustrate the model, we chose to only examine intensities at certain intervals (10 for field plots, 100 for lidar plots) in order to reduce the number of decision variables to a more feasible number. In other words, we considered 10, 20, 30, 40 … and up to 220 field plots and 180, 280, 380 … and up to 5480 lidar plots. This resulted in 25 choices for field plots ((250 − 10)/10 + 1 = 25) and 54 for lidar plots ((5480 − 180)/100 + 1 = 54)—a total of 1350 field vs. remote sampling combinations (25 × 54 = 1350).
The linear programming formulation is as follows:
min ( i , j ε i j t x i j ,     i , j ε i j w x i j ,     i , j ε i j s x i j )
subject to:
i , j c i j x i j B
i , j f i x i j F  
i , j l j x i j L  
i , j x i j = 1  
x i j { 0 , 1 }             i , j  
In the above series of equations, i and j denote the discrete intensity options for the number of field and lidar plots to sample, respectively. Objective function coefficients ε i j t , ε i j w ,   and ε i j s are the calculated standard errors for tree carbon, down woody material carbon, and soil carbon estimates (in tC/ha) at field plot intensity i and lidar intensity j. Constraint coefficient c i j is the cost of measuring i field plots and j lidar plots (in USD), and B is the total budget in USD. Coefficient f i is the time needed to sample i field plots (in hours), and F is the total number of hours that can be spent in the field in two seasons, assuming three field crews who work 8 h per day. Coefficient l j is the time needed to collect data from j lidar plots (in hours), and L is the total number of hours allocated to this purpose during the two-year time period. Equation (1) is a composite objective function that seeks to minimize the standard error in estimates of tree carbon, down woody material carbon, and soil carbon, respectively. Inequality (2) is the budget constraint that forces the sum of the total field plot measurement cost for i field plots and the cost of lidar analyses for j lidar plots to be below budget B. Inequality (3) ensures that the field plot measurement time fits within the time allocated for fieldwork during the two-year planning horizon (F). Inequality (4) ensures that the time needed to collect remotely sensed data does not exceed what is available in a two-year time window (L). Inequality (5) allows the model to choose only one combination of field and lidar plot densities. Equation (6) defines the decision variables x i j in the formulation as binary.

2.2. Values of Parameters and Constants

The standard error coefficients in objective function (1) were calculated as a function of field plot intensity i and lidar plot intensity j, based on the following formula developed for double sampling (Cochran, 1977) [15]:
ε i j k = ν k ( 1 ρ k 2 ) i + ν k ρ k 2 j ν k Z  
Parameter ν k is carbon variance for pool k, where k = t, w, and s for tree carbon, down woody material carbon, and soil carbon variances, respectively; ρ k is the lidar correlation coefficient for pool k, where k = t, w, and s for tree carbon, DWM carbon, and soil carbon correlations; and Z is the maximum number of lidar plots that can fit into the 13.7 million hectare area (assumed to be one per ha). We calculated the carbon variances and lidar correlations based on a preliminary analysis of tree carbon, down woody material carbon, and soil carbon using a pilot field and lidar dataset in the Tanana valley of interior Alaska [16]. Our estimates for tree carbon, down woody material carbon, and soil carbon variances were 845.8, 32.4, and 2575.6 Mg/ha, respectively.
The correlation between lidar metrics and field plot-based estimations for each carbon pool is also a component of the standard error, ε i j k , which denotes the standard error for carbon pool k using field plot intensity i and lidar plot intensity j. In this study, we assumed that lidar technology would be used for remotely sensed sampling. The lidar-based metric used for the lidar correlations was mean plot-level tree segment height. An individual tree segmentation algorithm was run on a 1 m lidar canopy height model, and the maximum height within each segment (above 1.37 m or breast height) was extracted as the segment height [17]. The plot-level lidar metric was then calculated as the mean lidar segment height for all segments within the 0.06744 ha footprint of the FIA plot (this corresponds to the combined area of all four subplots on an FIA plot). Although there are a large number of possible lidar metrics, we chose this metric, representing the mean tree height within the plot, since it is easy to interpret and has been shown to be highly correlated with aboveground biomass (AGB) in boreal forests [18]. For the purpose of this study, tree carbon was based on all trees with a breast-height diameter of 2.54 cm. Likewise, although there are many other possible approaches to estimating soil carbon and DWM carbon with satellite and/or airborne lidar, restricting our analysis to a single lidar-based estimation provides a straightforward yet realistic basis for comparing the relative costs and trade-offs of lidar and field sampling in an optimization framework. Furthermore, our intention in this paper is to demonstrate the application of a multi-objective optimization approach that could be parameterized with different standard error values, correlations, cost levels, and logistical constraints, depending on the specific inventory context within which it is applied. Based on our pilot observation data, our estimates for tree carbon, down woody material carbon, and soil carbon correlations were 0.91, 0.23, and 0.05 Mg/ha, respectively.
The rest of the model parameters were provided by USFS Forest Inventory and Analysis staff in Alaska. The total budget, B, USD 2 M. The amount of time to measure i field plots, f i was calculated as i × 8 h, and the amount of time to analyze j lidar plots, l j , was calculated as j × 2 h. Costs were assumed to be USD 9285 for each field plot and USD 100 for each lidar plot. The total time for fieldwork, F, and lidar analysis, L, were calculated by assuming a two-year time period with 45 field days in each year and 50 weeks of analysis time per year, assuming 40 h workweeks as the standard.

2.3. Finding Pareto Optimal Sampling Intensity Combinations

We solved the multi-objective program (1)–(6) for Pareto optimality by first enumerating all feasible solutions with respect to the budget and time constraints. Then, we checked which subset of these solutions were non-dominated or Pareto optimal. This was done by iterating through each feasible solution and testing if any of the other feasible solutions had smaller standard errors with respect to one of the three carbon pools without having a larger standard error with respect to any of the other pools. We were able to do this due to the relatively small number of discrete plot intensity combinations that we decided to analyze for illustration purposes. If a greater number of plot intensity combinations are to be evaluated, specialized multi-objective solution algorithms could be used, such as those discussed in Tóth et al. [19] or in Tóth and McDill [20].

2.4. Sensitivity Analyses

As briefly discussed in Section 2, estimates of variance for forest carbon pools are not always available in remote regions. In order to test our model’s robustness in the face of this uncertainty, we re-solved the model with perturbed variance estimates of 300, 400, 500, 600, 700, 800, 900, and 1000 Mg/ha for tree carbon (in addition to our best estimate of 845.8 Mg/ha), with 10, 20, 30, 40, 50, 100, 200, 300, 400, and 500 Mg/ha for down woody material carbon (in addition to our best estimate of 32.4 Mg/ha), and with 1500, 2000, 2500, 3000, 3500, and 4000 Mg/ha for soil carbon (in addition to our best estimate of 2575.6 Mg/ha). We also examined the effect of changing the carbon correlations between remote and field estimates of soil carbon to 0, 0.1, 0.2, 0.4, 0.6, and 0.9 (in addition to our best estimate of 0.05), to 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.88, and 0.9 for down woody material carbon (in addition to our best estimate of 0.23), and to 0.4, 0.5, 0.6, 0.7, 0.88, and 0.95 for tree carbon (in addition to our best estimate of 0.91). Last, we also re-solved the model for more optimistic budget scenarios of USD 2.25 M, USD 2.5 M, and USD 2.75 M (in addition to our baseline budget of USD 2.0 M).

3. Results

We found four Pareto optimal sampling intensity combinations (Table 1). These were 210 field plots and 480 lidar plots, 200 field plots and 1380 lidar plots, 190 field plots and 2280 lidar plots, and 180 field plots and 3280 lidar plots. At the baseline budget of USD 2 M, having more than 210 or fewer than 160 field plots is not optimal in any of the variance and correlation scenarios that we considered, even though there were feasible solutions at greater and smaller sampling intensities.
We also found that only the budget constraint affected feasibility; the time constraints were not binding in any of the variance and correlation scenarios. When we increased the budget, the number of Pareto optimal sampling intensity combinations first increased from four at USD 2 M to five at USD 2.25 M; then it went back down to four at USD 2.5 M and two at USD 2.75 M (Table 2).
As expected, the standard errors decreased with increasing budgets from a range of 1.01–1.47 tC/ha for tree carbon at USD 2 M to 0.9–0.99 tC/ha at USD 2.5 and 0.86 tC/ha at USD 2.75. For DWM carbon, the expected standard errors decreased from 0.39–0.41 tC/ha at USD 2 M to 0.35–0.37 tC/ha at USD 2.5 M and 0.35–0.36 tC/ha at USD 2.75 M. Last, for soil carbon, we found that the standard errors went down from 3.5–3.78 tC/ha at USD 2 M to 3.21–3.42 tC/ha at USD 2.5 M and 3.21–3.27 at USD 2.75 M. Please note that standard error ranges are provided because different Pareto optima at a given budget resulted in a range of different associated standard errors (Table 1 and Table 2).
Changing the variances of the three carbon pools did not affect which sampling schemes were Pareto optimal. Increasing the correlation between field and remotely sensed tree carbon estimates from 0.91 to 0.95, however, changed the set of Pareto optima from four to six field-remote sampling intensity combinations: (210, 480), (200, 1380), (190, 2280), (180, 3280), (170, 4180), and (160, 5080). Increasing the correlation between field and remotely sensed soil carbon estimates to 0.4 removed one efficient point (210 field plots and 480 lidar plots), while increasing all three correlations to 0.9 led to a unique efficient solution with 180 field plots and 3280 lidar plots. While we do not necessarily consider either of these scenarios to be realistic, the exercise does show that other Pareto optimal solutions can exist. Lastly, we note that the estimates of soil carbon always have the greatest associated standard errors, since correlation of soil carbon with lidar data is small.

4. Discussion and Conclusions

In this study, we illustrated how multi-objective mathematical programming can be used to optimize forest carbon sampling with a combination of field and lidar measurements in order to minimize the standard error of estimates of three specific carbon pools: tree carbon, down woody material carbon, and soil carbon pools. These results help illustrate the utility of lidar observations in supporting data collection in the face of limited time and budget constraints.
Using the Tanana Valley in Alaska, United States, as an illustrative case study, we identified four Pareto optimal combinations of field and lidar sampling intensities (Figure 2).
We also found, via sensitivity analyses, that these four sampling intensity combinations were rather robust to changing parameter values. Since none of these four sampling strategies consist solely of using either the field or the lidar method, we conclude that combining the two is the optimal course of action. This result matches what Stephen et al. [21] found for New Zealand forest inventories. Given the budget constraints, field sampling alone was not feasible. We also found that spending the entire budget on field sampling produces standard errors that can be reduced for all three carbon stocks simultaneously (unless the correlation between field and lidar observations is less than zero). In decision scientific terminology, this means that sampling strategies that comprise solely field measurements are dominated. The same is true for sampling schemes with very low field plot intensities, unless field and lidar observations are nearly perfectly correlated across carbon stocks. The location of the four Pareto optimal sampling strategies in 3-dimensional standard error space (Figure 2) suggests that choosing (210, 480) for field and lidar plot sampling intensities would likely be advantageous for the decision maker, as this option would lead to minimum standard errors for both down woody material carbon and soil carbon. Nonetheless, it is noteworthy that there is a lot of standard error space left between the (200, 1380) and the (210, 480) sampling intensity combinations. Alternatively, if the decision maker was particularly interested in minimizing standard error in tree C, given that this carbon pool is most influenced by dominant disturbance processes in interior Alaska (fires and insects) and management activities (fire suppression, bioenergy production, and timber harvesting), this could motivate them to choose the alternative with smaller standard error in tree C estimation (200, 1380) while still maintaining reasonable standard error for the DWM and soil C pools. If the decision makers were interested in exploring further options between these alternatives, then our model could be run for all nine field sampling intensities (201, 202 … 209) to see if any one of them would lead to Pareto optimality. We also note that the standard error ranges for the three carbon pools we examined are relatively narrow across the four Pareto optimal sampling intensity combinations. In other words, the tradeoffs with regard to our ability to accurately estimate carbon in these three pools are not that great. This finding speaks to the power of multi-objective optimization. No matter which of the four strategies the analyst chooses, the expected standard errors will be reasonable regardless of what the underlying carbon variances or correlations with lidar are, within extremely broad ranges.
As to which of the Pareto optimal sampling intensity combinations should be chosen, this is ultimately a decision that the inventory program (e.g., the USFS FIA) and their data user groups will have to make. While the national FIA program has established precision guidelines for area and volume estimates in the lower 48 states (Bechtold and Patterson, 2005) [7], these precision guidelines are not applicable to other inventory attributes (e.g., soil carbon) and in cases where the standard FIA sampling design is not logistically or economically feasible (e.g., interior Alaska). In addition, the Global Forest Observations Initiative has published methods and guidance to inform the development of monitoring, reporting, and verifying (MRV) systems in the context of REDD+ and national greenhouse gas reporting, including consideration of costs, sampling design, bias, and precision (GFOI, 2020) [22]. While these precision guidelines are useful points of reference, in the context of a multi-resource inventory (such as the FIA program) the decision maker must weigh the value of precision for a number of inventory attributes against the costs in arriving at an acceptable solution, taking into account the needs of a wide variety of stakeholders. Our contribution with the present work is to provide a mathematical framework to quantify and visualize the Pareto optimal tradeoffs that come with any one of those decisions. Unlike other multi-objective sampling optimization models that were documented in the literature (e.g., Lark 2006) [23], ours maps the tradeoffs associated with the expected standard errors across inventory attributes such as tree carbon, down woody material carbon, and soil carbon.
We acknowledge that when weighing which sampling scheme to employ, knowing which carbon stock is the most important in minimizing standard errors would be helpful, as standard errors do change under different assumptions of carbon variability and field-lidar correlation. To this end, we created an optimization framework that is flexible and can be updated as improved estimates of carbon stock variability and correlations between field and remotely sensed data become available. In addition, the proposed multi-objective mathematical programming method can optimize sampling for other metrics as well, such as invasive species cover or mean forest stand age. The model can also be used to choose optimal intensities for any two sampling strategies with known estimates of standard error rates and correlations. Thus, it would also be applicable to other remote sensing technologies that might become available in the future, as long as information about their error and correlation profiles is provided. More broadly, we suggest that the multi-objective linear programming framework we established here is even more generalizable to other spatial sampling applications where the logistical and budgetary constraints are entirely different. This is due to the general mathematical structure of linear programs, comprising a set of linear functions that can represent various sampling objectives accompanied by a set of linear inequalities that can represent various logistical and budgetary constraints. Another example of the flexibility of this mathematical structure is that the model can easily be modified to minimize the sum total of standard errors (Equation (8)) or the sum of weighted standard errors (Equation (9)), or even to minimize the maximum standard error over the three pools (also known as the minimax formulation; see Equation (10):
min ( i , j ε i j t x i j + i , j ε i j w x i j + i , j ε i j s x i j )
min ( w a i , j ε i j t x i j + w b i , j ε i j w x i j + w s i , j ε i j s x i j )  
where w a , w b , and w s denote the relative weights assigned to estimating tree carbon, DWM carbon, and soil carbon, and
m i n ( max ( i , j ε i j t x i j , i , j ε i j w x i j , i , j ε i j s x i j ) )
Whether any of these make sense will depend on the objectives and the unique features of the survey at hand. Although standard errors do propagate in an additive fashion across the pools, and thus the estimates of total carbon can be optimized by minimizing the total of standard errors, inventory analysts might still want to optimize the individual standard errors associated with each pool, thereby developing a more holistic understanding of the tradeoffs in the totality of the expected standard error space.
Lastly, we would like to emphasize that multi-objective linear programming is not the only optimization method that can be used to aid forest carbon sampling design. We would argue, however, that it is likely the best choice today for the following reasons. Ad hoc heuristics, such as Monte Carlo simulation, simulated annealing, tabu search, and genetic algorithms, cannot guarantee optimality, which is a problem when one seeks to make the best possible use of limited budgets for a specific purpose. Unlike heuristics, LPs are exact mathematical models that can find proven optimal solutions given enough computing time. Time can be a barrier to finding true optima when some of the variables in an LP are restricted to be discrete (also known as integer programs, or IPs) and the size of the problem is large in terms of the number of variables. However, even under these circumstances, IPs provide an upper bound on the gap that remains between the true optimum and the best incumbent found in the allotted time. Thus, one would know how good the solution is in terms of optimality. Heuristics cannot provide such information. What they can do, though, is find good solutions quickly for some hard problems that might otherwise be intractable by IPs (some combinatorial problems fall in this category, where worst-case IP solve times are more than a polynomial function of problem size; these problems are called NP-Hard). However, the forest carbon inventory problem we examined in this paper is computationally trivial; we were able to find the true optima in a fraction of a second under all parameter value combinations. Apart from being able to guarantee optimality or an upper bound on the objective value of the incumbent, perhaps the biggest advantage of LPs over heuristics is that encoding and solving LP models do not require custom application development. There are many off-the-shelf solvers that are easy to access and free to use (e.g., Excel Solver). Model representation, parameterization, the addition of constraints and variables, and sensitivity analyses are trivial and streamlined for the average user. With heuristics, on the other hand, not only does the developer have to build the model using custom application, but the developer must also develop and tune the solution algorithm itself. While the general algorithmic design is widely available for most popular heuristics (e.g., simulated annealing, or SA, and evolutionary algorithms), coding and executing these algorithms is not trivial. Even if the developer succeeds in finely tuning the parameters (e.g., setting the temperature parameter in SA), any addition of new constraints or variables to the model will require going back to the drawing board to adjust the code and re-tune the parameters. This is a big issue for sampling optimization, because logistical and budgetary constraints must often be added, modified, or dropped on a case-by-case basis. Unlike the situation with LP solvers, this can only be done by the developer. Last, the computational performance of software and hardware resources have been increasing at a breakneck speed for the last 30 years, rapidly making ad hoc custom heuristics less and less advantageous when compared with exact mathematics in optimization.
In sum, the principal value of our work is a novel application of multi-objective mathematical programming that not only allows forest inventory analysts to determine what combination of field and lidar plot intensities should be used for sampling design to minimize expected standard errors in estimating boreal forest carbon pools, but also to streamline what prior data must be collected or estimated, and with what accuracy, in order to find optimal sampling designs. As an example, our demonstrative case study in the Tanana Valley, Alaska, United States revealed that the accuracy in estimating variances in the amount of carbon in the various pools via pilot surveys is much less important for finding optimal sampling designs than having more accurate estimates of the budget that is available for the surveys.

Author Contributions

S.F.T., K.L.O. and C.C.S. designed the study. H.-E.A. acquired and pre-processed the Alaska forest dataset. K.L.O. and C.C.S. built the program to solve the proposed multi-objective model. K.L.O. and C.C.S. wrote the initial draft. S.F.T. implemented additional model runs and revised the manuscript. H.-E.A. interpreted and contextualized the results. He also provided substantial revisions. All authors have read and agreed to the published version of the manuscript.

Funding

The United States Forest Service, Pacific Northwest Research Station, provided financial support for this study. One of the co-authors, H.-E.A., is an employee of the funding body. H.-E.A. provided some of the input data for the study, participated in analyzing the results, and helped revise the manuscript for submission.

Data Availability Statement

The datasets used in the present study are available from the corresponding author on reasonable request.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Keyser, A.R.; Kimball, J.S.; Nemani, R.R.; Running, S.W. Simulating the effects of climate change on the carbon balance of North American high-latitude forests. Glob. Chang. Biol. 2000, 6, 185–195. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Tans, P.P.; Fung, I.Y.; Takahashi, T. Observation constraints on the global atmospheric CO2 budget. Science 1990, 247, 1431–1438. [Google Scholar] [CrossRef] [PubMed]
  3. Watson, R.T. Land Use, Land-Use Change, and Forestry; Cambridge University Press: Cambridge, UK, 2000. [Google Scholar]
  4. Barrett, T.; Christensen, G. Forests of Southeast and South-Central Alaska, 2004–2008: Five-Year Forest Inventory and Analysis Report; General Technical Report; Pacific Northwest Research Station, USDA Forest Service: Portland, OR, USA, 2011.
  5. Zech, R.; Huang, Y.; Zech, M.; Tarozo, R.; Zech, W. High carbon sequestration in Siberian permafrost loess-paleosols during glacials. Clim. Past 2011, 7, 501–509. [Google Scholar] [CrossRef] [Green Version]
  6. Intergovernmental Panel on Climate Change. Climate Change 2007: Mitigation of Climate Change: Contribution of Working Group III to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, UK, 2007. [Google Scholar]
  7. Bechtold, W.; Patterson, P. The Enhanced Forest Inventory and Analysis Program: National Sampling Design and Estimation Procedures; USDA Forest Service, Southern Research Station: Asheville, NC, USA, 2005; Volume 80.
  8. Smith, W.B.; Miles, P.D.; Perry, C.H.; Pugh, S.A. Forest Resources of the United States, 2007; General Technical Report WO-78; U.S. Department of Agriculture, Forest Service, Washington Office: Washington, DC, USA, 2009; 336p.
  9. USDA. Forest Inventory and Analysis Fiscal Year 2018 Business Report. 2018. Available online: https://www.fia.fs.fed.us/library/bus-org-documents/docs/17973%20FS%20FIA%20Fiscal%20Year%202018%20Business%20Reportv3%20508.pdf (accessed on 23 September 2020).
  10. UNFCCC. Warsaw Framework for REDD+ Decision 11/CP.19. 2013. Available online: http://unfccc.int/resource/docs/2013/cop19/eng/10a01.pdf#page=31 (accessed on 17 December 2020).
  11. FAO. Global Forest Resources Assessment 2020–Key Findings; FAO: Rome, Italy, 2020. [Google Scholar] [CrossRef]
  12. UNFCCC. Adoption of the Paris Agreement. In Proceedings of the 21st Conference of the Parties, United Nations, Paris, France, 12 December 2015. [Google Scholar]
  13. Barrett, T.; Andersen, H.-E.; Winterberger, K. Integrating field and lidar data to monitor Alaska’s boreal forest. In Proceedings of the IUFRO Division 4 Conference: “Extending Forest Inventory and Monitoring over Space and Time”, Quebec City, QC, Canada, 19–22 May 2009. [Google Scholar]
  14. Köhl, M.; Lister, A.; Scott, C.T.; Baldauf, T.; Plugge, D. Implications of sampling design and sample size for national carbon accounting systems. Carbon Balance Manag. 2011, 6, 10. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Cochran, W. Sampling Techniques, 3rd ed.; Wiley: New York, NY, USA, 1977. [Google Scholar]
  16. Pattison, R.; Andersen, H.-E.; Gray, A.; Schulz, B.; Smith, R.J.; Jovan, S. Forests of the Tanana Valley State Forest and Tetlin National Wildlife Refuge, Alaska: Results of the 2014 Pilot Inventory; General Technical Report PNW-GTR-967; U.S. Department of Agriculture, Forest Service, Pacific Northwest Research Station: Portland, OR, USA, 2018; 80p.
  17. Vincent, L.; Soille, P. Watersheds in Digital Spaces: An Efficient Algorithm Based on Immersion Simulations. IEEE Trans. Pattern Anal. Mach. Intell. 1991, 13, 583–598. [Google Scholar] [CrossRef] [Green Version]
  18. Næsset, E.; Gobakken, T. Estimation of above- and below-ground biomass across regions of the boreal forest zone using airborne laser. Lidar Environ. 2008, 112, 3079–3090. [Google Scholar] [CrossRef]
  19. Tóth, S.F.; McDill, M.E.; Rebain, S.A. Finding the Efficient Frontier of a Bi-Criteria, Spatially Explicit, Harvest Scheduling Model. For. Sci. 2006, 52, 93–107. [Google Scholar]
  20. Tóth, S.F.; McDill, M.E. Finding Efficient Harvest Schedules under Three Conflicting Objectives. Forest Sci. 2009, 55, 117–131. [Google Scholar]
  21. Stephens, P.R.; Kimberley, M.O.; Beets, P.N.; Paul, T.S.; Searles, N.; Bell, A.; Brack, C.; Broadley, J. Airborne scanning lidar in a double sampling forest carbon inventory. Remote Sens. Environ. 2012, 117, 348–357. [Google Scholar] [CrossRef]
  22. GFOI. Integration of Remote-Sensing and Ground-Based Observations for Estimation of Emissions and Removals of Greenhouse Gases in Forests: Methods and Guidance from the Global Forest Observations Initiative; Group on Earth Observations: Rome, Italy, 2020; Available online: https://www.reddcompass.org/mgd/resources/GFOI-MGD-3.1-en.pdf (accessed on 3 June 2022).
  23. Lark, R.M. Multi-objective optimization of spatial sampling. Spat. Stat. 2016, 18, 412–430. [Google Scholar] [CrossRef] [Green Version]
Figure 1. The Tanana Unit study area in Alaska, United States. Approximate extent of forest cover shown in green.
Figure 1. The Tanana Unit study area in Alaska, United States. Approximate extent of forest cover shown in green.
Forests 13 00972 g001
Figure 2. Pareto optimal sampling intensity combinations for field and lidar plots with projected standard errors in tree carbon, down woody material carbon, and soil carbon.
Figure 2. Pareto optimal sampling intensity combinations for field and lidar plots with projected standard errors in tree carbon, down woody material carbon, and soil carbon.
Forests 13 00972 g002
Table 1. Pareto optimal sampling schemes and their associated standard errors using parameter estimates based on observed data with a budget of USD 2,000,000.
Table 1. Pareto optimal sampling schemes and their associated standard errors using parameter estimates based on observed data with a budget of USD 2,000,000.
BUDGET: USD 2,000,000
Standard Error (tC/ha)Intensity (Number of Plots)Total Sampling Cost (USD)
Tree CarbonDown Woody Material CarbonSoil CarbonFieldLidar
1.1.466790.386903.49962210480$1,997,850
21.111060.393293.584732001380$1,995,000
3.1.035550.402813.677571902280$1,992,150
4.1.010550.413523.778211803280$1,999,300
Table 2. Pareto optimal sampling schemes and their associated standard errors using parameter estimates based on observed data with a budget of USD 2.25 M, USD 2.5 M, and USD 2.75 M.
Table 2. Pareto optimal sampling schemes and their associated standard errors using parameter estimates based on observed data with a budget of USD 2.25 M, USD 2.5 M, and USD 2.75 M.
Standard Error (tC/ha)Intensity (Number of Plots)
Tree CarbonDown Woody Material CarbonSoil CarbonFieldLidar
BUDGET: USD 2,250,000
1.2.120590.370653.27726240180
21.131640.367433.343062301080
3.1.007250.374633.417762201980
4.0.962980.383013.498012102980
5.0.952590.392263.584312003880
BUDGET: USD 2,500,000
1.0.987420.351713.206262501780
20.931180.358463.272162402680
3.0.909790.365923.342442303580
4.0.903970.373983.417492204480
BUDGET: USD 2,750,000
1.0.863220.350923.205932504280
20.860790.358033.271992405180
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Tóth, S.F.; Oken, K.L.; Stawitz, C.C.; Andersen, H.-E. Optimal Survey Design for Forest Carbon Monitoring in Remote Regions Using Multi-Objective Mathematical Programming. Forests 2022, 13, 972. https://doi.org/10.3390/f13070972

AMA Style

Tóth SF, Oken KL, Stawitz CC, Andersen H-E. Optimal Survey Design for Forest Carbon Monitoring in Remote Regions Using Multi-Objective Mathematical Programming. Forests. 2022; 13(7):972. https://doi.org/10.3390/f13070972

Chicago/Turabian Style

Tóth, Sándor F., Kiva L. Oken, Christine C. Stawitz, and Hans-Erik Andersen. 2022. "Optimal Survey Design for Forest Carbon Monitoring in Remote Regions Using Multi-Objective Mathematical Programming" Forests 13, no. 7: 972. https://doi.org/10.3390/f13070972

APA Style

Tóth, S. F., Oken, K. L., Stawitz, C. C., & Andersen, H. -E. (2022). Optimal Survey Design for Forest Carbon Monitoring in Remote Regions Using Multi-Objective Mathematical Programming. Forests, 13(7), 972. https://doi.org/10.3390/f13070972

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop