LiDAR-Assisted Multi-Source Program (LAMP) for Measuring Above Ground Biomass and Forest Carbon

Abstract

Forest measurement for purposes like harvesting planning, biomass estimation and mitigating climate change through carbon capture by forests call for increasingly frequent forest measurement campaigns that need to balance cost with accuracy and precision. Often this implies the use of remote sensing based measurement methods. For any remote-sensing based methods to be accurate, they must be validated against field data. We present a method that combines field measurements with two layers of remote sensing data: sampling of forests by airborne laser scanning (LiDAR) and Landsat imagery. The Bayesian model-based framework presented here is called Lidar-Assisted Multi-source Programme—or LAMP—for Above Ground Biomass estimation. The method has two variants: LAMP2 which splits the biomass estimation task into two separate stages: forest type stratification from Landsat imagery and mean biomass density estimation of each forest type by LiDAR models calibrated on field plots. LAMP3, on the other hand, estimates first the biomass on a LiDAR sample using models calibrated with field plots and then uses these LiDAR-based models to generate biomass density estimates on thousands of surrogate plots, with which a satellite image based model is calibrated and subsequently used to estimate biomass density on the entire forest area. Both LAMP methods have been applied to a 2 million hectare area in Southern Nepal, the Terai Arc Landscape or TAL to calculate the emission Reference Levels (RLs) that are required for the UN REDD+ program that was accepted as part of the Paris Climate Agreement. The uncertainty of these estimates is studied with error variance estimation, cross-validation and Monte Carlo simulation. The relative accuracy of activity data at pixel level was found to be 14 per cent at 95 per cent confidence level and the root mean squared error of biomass estimates to be between 35 and 39 per cent at 1 ha resolution.

1. Introduction

Traditionally large area forest inventory has been carried out by systematically designed field campaigns [1,2]. They result in district level and national forest inventories (NFI). Recently, Light Detection and Ranging (LiDAR) data is often used to improve the accuracy of the forest inventory parameter estimation in such campaigns. Estimation approaches vary from design-based to model assisted and model-based estimates. Design-based estimates establish a probabilistic sampling design that can be proven to be asymptotically un-biased when the number of plots increases. Typical such designs are systematic, simply random and clustered sampling designs. Model assisted sampling designs also establish an unbiased estimator but instead of direct probability sampling, a model can be used to create a sample with uneven selection probabilities. Popular model assisted sampling designs are two-phased sampling designs and designs that rely on a linear model based on some auxiliary variables. In all these cases, mean and variance estimators can be explicitly constructed [3]. Model-assisted estimation is a further extension of such methodologies where multiple different sources, such as field plots and LiDAR sampling, can be combined for producing unbiased estimates of forest parameters, as in [4].

In model-based estimation, the viewpoint changes from using models to assist sampling based estimation into using a model that is based on some auxiliary variables, such as LiDAR metrics derived from the point cloud, as the primary source of information. Field plots no longer provide the data for the estimates but instead they become a teaching set for the model so that model parameters can be estimated. This latter paradigm is at the heart of Bayesian estimation theory, where the model is used to construct a prior estimator, to which field samples are coupled via a likelihood function to produce the posterior estimate that is a weighted average of the model forecast and the field plots, with weights inversely proportional to a prior estimate of error in the model and in the measurements, respectively. As soon as the model chosen becomes non-linear, it is very challenging to provide explicit mean and variance estimators. Instead, the burden of proof of the unbiased nature of the estimator is shifted to the calibration plot set that can be sampled based on a design or with a model assisted design. If this is the case, then asymptotic unbiasedness is tested by cross-validation methods, such as Leave-One-Out (LOO). Some recent Bayesian studies of forest inventory that have adopted this approach include [5,6,7].

Another important factor in forest inventory has been the goal of the inventory. In the National Forest Inventory context, the goal is accurate and precise estimates of total tree population statistics for large areas. By implication these estimates should then have small systematic error and small standard deviation. For operational forest planning, the primary need is accurate and precise estimates on much smaller areas, such as forest stands. But it is also possible to bridge the two goals and produce accurate and precise estimates along a range of scales from forest stands to entire countries or provinces. In this case model-based estimation is in practice the only choice, because the cost of dense field surveys is prohibitive.

In [4] the authors studied model-assisted approach to estimation and inference when using LiDAR as a tool to inventory Above Ground Biomass (AGB). The article features a methodological presentation of the estimators of total biomass and biomass per hectare as well as variance estimators. In [8] the authors concluded that the most urgent problem facing LiDAR-assisted estimation based on systematic sampling is the very large overestimation of estimator variance by assuming a simple random sampling design. In [9] the authors observed that standard errors were consistently lower for the LiDAR-assisted designs than those based purely on sample plots. Results shown in [10] indicate that an ALS-based survey produces valid inference under design-based and model-based frameworks. The authors of [11] determined whether profiling LiDAR data can be used to improve the precision of an existing ground survey. The preliminary results show that use of model-dependent, two-phase approach can improve the precision of the AGB estimate when considering smaller geographical areas, i.e., smaller political units (∼1000 km${}^2$). In another study, it was discovered that the model precision was improved when using a pixel-based regression model estimated from scanner data rather than using the models estimated from profiler data [12].

However, if we hope to measure forest biomass for the purpose of measuring the carbon captured in tropical forests and for their potential in mitigating climate change, a much finer spatial resolution is needed to be able to address and quantify forest degradation at a local level. A few attempts have been made to produce global or national biomass maps based on satellite data alone [13,14,15,16] but use of satellite data without local ground truth they may contain large systematic errors [17]. For this reason, combining several data sets that feature both field plots, satellite imagery and recently also Airborne Laser Scanning (ALS or LiDAR) have become attractive. The articles [18,19,20,21] summarize many recent efforts in this direction. Even national scale biomass maps have been produced [22,23,24,25,26]. Other studies on using LiDAR sampling to assist in tropical biomass estimation include [27,28]. Also other remote sensing sources have been proposed for biomass estimates, such as SAR interferometry [29,30]. At the same time, several new guidelines that cover several methodologies have been assembled [31]. However, in most cases of such estimation processes that involve at least three spatial scales and two hierarchical steps of inference, comprehensive and reliable uncertainty analyses have been difficult to create, even if several authors have made significant steps in this direction [13,32,33]. A remarkably comprehensive study of estimation errors of LiDAR-based estimation of biomass, applied to both tropical and boreal forests, and across multiple spatial scales, is presented in articles [34,35]. The analysis presented in them of LiDAR-based model uncertainty starts from tree-level biomass estimation errors and scales up to millions of hectares. Some of the principal findings of these two studies are that tree-level errors only dominate the error budget of the estimation at the very local level of trees. From plot level upwards, the error budget is first dominated by model residual error up to a scale of 380 m, or some 15 ha, and from that scale onwards by model parameter estimation errors [35].

A recent study produced theoretical variance and mean error estimators for many model-based forest measurement approaches that involve hierarchical model building in cases where both field samples, samples of LiDAR-scanned areas and Landsat imagery are present [36]. This paper features a Bayesian framework for forest estimation, as does the current article, but the formulas estimating prediction variance are classical ones. When studying the uncertainty of model-based forest maps at different scales, one faces two principal approaches: validation by independent samples that may be probabilistic, or error propagation. The latter can be done in cases where models are linear orcan be linearized by classical variance estimation, whereas in the case of more complex and non-linear models the prevailing approach is Monte Carlo simulation.

In [34], where the study site was a tropical forest, relative prediction uncertainty was relatively high, staying at 20% level at 1 ha resolution. While this is partly due to the small size of validation plots, it is also to some extent caused by the simple two-parameter LiDAR-based model used in the study. Non-linear models with many more covariates, such as the ones discussed in this article, have mutually dependent covariate error distributions. If a linear variance analysis is properly applied to them, e.g., by linearizing model nonlinearities first by Taylor series, they typically yield very wide confidence intervals that do not reflect the true uncertainty of the estimates. Monte Carlo simulation and cross-validation therefore are often the most appropriate tools available to test the uncertainty of nonlinear multivariate models.

As discussed in [37], rational decision related to the maintenance and enhancement of the multiple functions, such as carbon storage, provided by the forests needs to be based on objective, reliable information. Integration of mapping and inventory provides an effective framework for the support of forest management [38]. The research reported in this article is an attempt to combine high spatial resolution with good overall accuracy for the purpose of mitigating climate change under the auspices of the United Nations programme for reducing emissions from deforestation and forest degradation, the so-called REDD+ scheme that is part of the Paris Climate Agreement and a number of earlier climate summits. REDD+ explicitly requires that the uncertainty of of forest inventory results is estimated and this is also the goal in the current research effort. The family of methods proposed here combines several data sources, namely a field plot campaign, sampling by LiDAR and optical satellite imagery.

Principles of REDD+ and Monitoring, Reporting and Verification (MRV)

A combination of remote sensing and ground-based forest carbon inventory approaches is necessary to design and implement an efficient REDD+ MRV system. There are two principal methods for accounting annual GHG emissions and removals in forest areas. One of them is the gain-loss method relying on activity data monitoring and emission/removal factors. Activity data monitoring is based on land use and land cover change analysis by activities over a reference period. Emission factors refer to emissions of GHGs per unit of activity, and removal factors to emission removals of GHGs per unit of activity. The other method is the stock change method which is used to estimate net annual emissions or removals from the difference in total carbon stocks at two points in time. National Forest Inventory (NFI) permanent sample plot measurements including the relevant carbon pools provide a basis for these kinds of difference analysis.

IPCC describes methods at three levels of detail, called tiers [39]. Tier 1 is considered as the default method, and are often based on global datasets to acquire activity data and emission and removal factors by broad forest types. Tier 2 usually uses the same mathematical structure as Tier 1, but the countries apply data specific to their national circumstances. These methods require conducting field inventory campaigns to estimate the values required if they are not available from historical archives. Tier 3 methods are generally more complex, normally involving modelling and higher resolution land use and land-use change data. IPCC expects that higher Tier (Tier 2 or Tier 3) methods are targeted unless the demanded efforts exceed the overall benefits that can be achieved.

Measuring, Reporting and Verification (MRV) of carbon capture in forests—whether tropical, sub-tropical, temperate or boreal—calls for forest biomass estimation methods that need to be accurate, precise and affordable for REDD+ to become an effective vehicle in helping humankind to mitigate climate change. High accuracy means minimal error in estimates, so that tons of carbon can be reliably converted to financial rewards. Moreover, it is highly desirable to have high spatial resolution that is needed to target climate mitigation efforts so that performance can be rewarded even at a local level. Affordability allows frequent, for example biennial, updating of change in biomass and correspondingly frequent compensation for performance. A family of methods with the generic title of LiDAR-Assisted Multi-source Program (LAMP), have been proposed that fulfil these requirements. LAMP methods combine systematic field measurement campaigns, airborne LiDAR sampling and geographic extension of both of these by satellite imagery.

2. Materials and Methods

In this study we generated a Reference Level (RL) and a biomass map at 1 hectare-scale for the sub-national region Terai Arc Landscape (TAL) of Nepal. This section introduces the study area and the data sources used, in particular field sample plots, LiDAR scanning data and Landsat satellite imagery. Two methods (LAMP2 and LAMP3) are described that were used to calculate the RL and the biomass map, respectively.

2.1. LiDAR-Assisted Multisource Programme

The family of methods introduced here, called LAMP for LiDAR-Assisted Multi-source Program for MRV, a term coined by Eric Dinerstein and George Powell, then with WWF US, seeks to fulfil the three requirements of efficiency, effectiveness and affordability [40]. It achieves this by combining an initial, relatively small field campaign with also an initial affordable Airborne Laser Scanning (ALS or LiDAR) sampling campaign with spatial extension and follow-up measurement based on satellite imagery. Also calculating the past time series of Reference Levels (hereafter RL’s) will be based on satellite imagery, typically Landsat imagery.

In order to ensure that the three different kinds of materials used to produce RL estimates keep their accuracy, they are incorporated in a unified statistical framework. This framework shares some features of both design-based, model-assisted and model-based spatial statistical paradigms [2,3,41,42], but on the bottom line it is model-based.

Measurements used in LAMP are first initiated with a statistical sampling of blocks or strips to be flown with LiDAR, typically covering between one and five percent of the total forest area. The sampling may be weighted to ensure that all significant forest types are covered by LiDAR. After LiDAR blocks have been assigned, an either systematic, simple random or clustered random sampling of field plots is conducted on the blocks. These sample plots are used to calibrate a Bayesian model of AGB that will be used to estimate AGB from LiDAR metrics.

After AGB estimates are produced on blocks of high spatial resolution (e.g., 20 m by 20 m), LiDAR blocks will further be sampled in simple random sampling for 1 ha sample areas or as surrogate plots that will be used to contribute to the estimation of RL’s from satellite data as the teaching set to replace real field plots. The 1 ha size of these sample areas was chosen as a compromise between observed LAMP model noise and desired spatial resolution that targets decision making at areas that can be regarded as forest stands. In other conditions, such as in Finnish boreal forests, forest stands are typically of 1 hectar size. A similar choice of finest scale was also made in a recent study of biomass mapping across multiple scales in [35]. AGB estimates on these sample areas are aggregated from the AGB values on the LiDAR estimation grid. The size of the sample areas can be chosen to correspond to the size of the 1 ha primary estimation units of the AGB map that are calculated on a 100 m by 100 m grid.

At this point LAMP methods branch at a junction into LAMP2, for Tier 2 level RLs, and LAMP3 for Tier 3 level RLs with a very fine spatial resolution. In LAMP2 methods, analysis of satellite images is used to stratify the forest into forest types and two forest condition classes for each forest type: intact and degraded. The areas covered by each stratum are used as Activity Data, and Emission Factors are calculated as AGB averages of sample areas that fall into each stratum. Other carbon pools are modelled by empirical functions from AGB values. With LAMP2, RLs can be calculated at district level. A diagrammatic overview of the LAMP2 algorithm is presented in Appendix A.

In LAMP3, both Activity Data and Emission Factors are deduced directly from AGB estimates on the chosen estimation units for the satellite model that is a Bayesian model of the same kind as that used for LiDAR based estimates. Circular 1 ha surrogate plots are used as calibration plots for this model. In LAMP3, Emission Factors are assigned at the level of primary estimation units, so they can vary even at 1 hectare spatial scale. A particular challenge with LAMP3 is the complex dependence of satellite image band values of optical satellites on imaging conditions, season and time. For this reason, satellite images must be carefully mosaicked before LAMP3 can be applied.

In the study reported in the remainder of this article, both kinds of LAMP methods are described in the section Materials and Methods, along with the description of their application in the TAL jurisdiction in Southern Nepal.

Table 1. Overview of LiDAR-Assisted Multi-source Program (LAMP)2 and LAMP3 algorithms.

Material:		Section
LiDAR:	LiDAR data (5% coverage)	2.3
Field data:	Vegetation plots (738 plots of 500 m${}^2$ within LiDAR coverage)	2.4
Satellite:	Satellite data (medium resolution such as Landsat, 100% coverage)	2.5.1
Step	Contents
1.	Stratify the forest of the study area into the main forest types and forest condition	2.3, 2.5.1,
	classes using satellite data (= forest strata map). (Satellite)	2.5.4
2.	Sampling of locations for LiDAR data acquisition and field plot collection. Weighted	2.3
	random sampling by incorporating the forest strata map, covering all important
	forest types.
3.	Calibrate LiDAR-to-AGB model with field based AGB. (LiDAR and Field data)	2.5.4
4. LAMP2	Randomly select 1000 circular LiDAR sample areas of 1 ha size for each forest	2.5.4
	strata within the LiDAR-area. Purpose: They will be used for calculating a mean
	biomass value for each stratum (forest type and condition class).
4. LAMP3	Select 10,000 circular LiDAR sample areas (surrogate plots) using a weighted	2.6.2
	random sampling within the LiDAR area. Weights should be the inverse of LiDAR
	block sampling. Purpose: To be used as training data (surrogate field data) in
	satellite-based model.
5.	Use LiDAR-to-AGB model to estimate AGB for the LiDAR sample areas (LAMP2)	2.5.4,
	or surrogate plots (LAMP3)	2.6.2
6. LAMP2	Calculate a mean AGB value for each stratum from the LiDAR-model estimates	2.5.4
	on LiDAR sample areas. To be used for calculating Emission Factors. Combine
	these forest class-specific mean AGB values with the forest strata map of
	the entire area.
6. LAMP3	1. Extract satellite-based features (band values, vegetation indices) from mosaicked	2.6.1
	satellite-imagery of the entire area. (Satellite)
	2. Calibrate Satellite-to-AGB model with the surrogate plot AGB estimates.	2.6.3
	3. Estimate AGB for each satellite image pixel with the Satellite-to-AGB model.
	4. Post-process the satellite data based AGB estimates with histogram matching	2.6.4
	method to avoid saturation effect.
7.	The previous steps result in mapped AGB for the entire area, at strata level (LAMP2)	2.5.4, 2.6.3,
	or at 1 ha level (LAMP3), respectively.	2.6.4
8. LAMP2	Time-series analysis of satellite data to generate Activity Data for Reference Level,	2.5.5–2.5.8
	using stratified satellite imagery of two successive time instances T1 and T2.
8. LAMP3	Time-series analysis based on AGB value differences at 1 ha grid level, estimated	2.6
	with the Satellite-to-AGB model from mosaicked satellite-imagery of the entire area
	over the whole time period.

presents a step-by-step overview of both methods, indicating also the sections where the two methods differ, with reference to the section below where each step is described in more detail.

LAMP2 was used for comprehensive analysis and calculation of RLs for the whole TAL region over a twelve year period. LAMP3 was used only for producing high resolution AGB maps for a single year. A series of subsequent analyses, including several validation field campaigns, cross-validation studies and Monte Carlo simulation of error accumulation were conducted to verify the accuracy of the estimation processes and they are reported in the last section of Results to provide the uncertainty analysis of especially the LAMP2 method. A summary of this process and the resulting Reference Level analysis has been reported in [40] where much technical detail was omitted because of page restrictions. An application of LAMP towards analyzing the dependence of AGB estimates on field plot distance from roads in TAL have been reported in [43].

2.2. Study Site

Terai Arc Landscape (TAL) is situated along the foothills of the Himalayas in the southernmost part of Nepal, ranging from the lowlands of Terai region up to the southern slopes of the Himalayas in Churia hills. The average altitude varies from less than 100 m to 2200 m [44]. The TAL area is influenced by tropical and subtropical climate. About half of the TAL is covered by subtropical, mainly deciduous forests. The dominating forest types are sal (Shorea robusta) terai mixed hardwood, khair-sisau (Acacia catechu/Dalbergia sissoo) and chir pine (Pinus roxburghii).

According to [45], about 1.18 million ha (51.5%) of the total land area was under forest cover in TAL in 2001. About 79% (0.9 million ha) of the forest is located outside of protected areas and 21% (0.3 million ha) is within protected areas [46]. In 2013, about 241,484 ha of forest were under the community forest management regime (i.e., 20.5% of the total forest area) and about 45,154 ha of forest were under the collaborative forest management regime (i.e., 3.8% of the total forest area) [47]. The remaining forests are mostly government-managed forest. Sal (Shorea robusta) is the dominant species found in the TAL. In the recent forest resource assessment (FRA) project, the Terai forest was classified into four major types: Sal Forest, Terai Mixed Hardwood Forest, Sal Mixed with Terai Hardwood Forest, and Khair-Sissoo.

The TAL is linked with eleven trans-boundary protected areas across Nepal and India. TAL is home to flagship species like tigers, rhinos, Asiatic wild elephants, and many other endangered species. This landscape has the second largest population of rhinos and one of the highest densities of tiger populations in the world. TAL plays an important role in maintaining linkage among these eleven protected areas in Nepal and India. Along the corridor, the forests connecting these protected areas vary from dense intact forests to degraded forest patches. Due to human pressure on forest resources, the forest cover of low land Terai and Churia has decreased during the last three decades. Connectivity among protected areas is crucial for effective and sustainable landscape level conservation. Hence, TAL program was started to create a single landscape level functioning unit by connecting 11 protected areas in Nepal and India through restoration of degraded forest corridors. This would give rise to trans-boundary dispersal corridors and migration paths for tigers, rhinos, elephants and many other species, which are crucial for maintaining biological diversity and gene flow. Conservation of the Churia forests is crucial for preventing soil erosion, flash floods, and recharging the water table of the Terai, the most productive land in the country. Therefore, sustainable management of TAL will help maintain biological diversity and also meet national demand of forest products and food supply for its rapidly growing human population.

2.3. Conducting the LiDAR Campaign

The National Forest Inventory (NFI), 1994 defined "forests" as having a crown cover >$10\%$ and an area >1 hectare. The forest in the study area was stratified using forest type map of TAL based on LANDSAT 7 [48], with an overall accuracy of 84.5% (Kappa = 0.75) at 30 m resolution, to produce a LiDAR sample that reflects the full range of variation in biomass over the study area. Different weights were assigned to the grid cells based on importance of forest types and amount of remaining forest in each type. The weights were scaled to sum up to one. To give higher weight to rare vegetation types, the initial expert weights were then scaled and normalized by the inverse of area fraction of the vegetation types from total area. Probability proportional-to-size sampling [3] based on the forest type map was used to select 20 sampling areas (5 km × 10 km blocks) covering 5% of the study area for LiDAR acquisition.

Prior to the field campaign to measure calibration plots for the LiDAR model, the location of sample plots was designed using systematic cluster sampling within rectangular sample areas, please see Figure 1. The airborne laser scanning (ALS) campaign was carried in March/April 2011. All blocks were scanned in full coverage from 2200 m average height above ground using a local helicopter equipped with a Leica ALS50-II LiDAR-scanner device. Resulting nominal outgoing pulse density at ground level was in average 0.8 points per square meter. The collected LiDAR data was evaluated after each flight, and supporting scans were conducted if data gaps or other problems occurred.

Each designed LiDAR block contained six clusters of eight sample plots each (Figure 2). The distance between cluster centres was 3333 m in West-East and 2500 m in North-South direction. Within the clusters, the sample plots were aligned in two parallel columns in North-South direction, with four plots per column. The distance between plots was 300 m in West-East direction, and 300 m and 150 m in North-South direction in Terai and Siwaliks, respectively. The smaller North-South distance for Siwaliks was chosen because of the large variations in altitude in the mountainous region. The plots were of fixed circular shape with 12.62 m radius, equivalent to an area of 500 m${}^2$.

2.4. Field Campaigns

The field data was collected in March/April 2011 in Terai and April/May 2011 in Siwaliks. Plot centre coordinates were recorded using differential GPS with ProMark 3 and MobileMapper CX devices, and corrected in post-processing mode (GNSS Solutions software (version 3.10.13) and MobileMapper Office software (version 3.40a)). 792 sample plots that were located in forest with at least 10% canopy cover were measured in the field. The measurements at tree-level included all living trees and shrubs above 5 cm diameter within the plot area.

The tree-level data was divided into 23 different tree species groups. For each field sample plot the following attributes were derived from the tree-level measurements, by species group and totals: Stem count (1/ha), mean diameter at breast height weighted by basal area (cm), basal area (m${}^2$/ha), mean tree height weighted by basal area (m), stem volume (m${}^3$/ha), and AGB (tons/ha). While mean diameter at breast height (dbh) and mean basal area were a direct output of the field measurements, mean tree height, mean volume and mean biomass were estimated using species group-specific functions and coefficients. In the following we explain in more detail how these estimates were derived.

Tree height was measured in the field only for a subset of trees per plot. Mean tree height per plot was then calculated using species group-specific height-diameter relationships with non-linear mixed-effect models. Mixed-effects models are an appropriate tool for modelling the relationship between tree height and field-measured tree diameter because the explanatory variables are clustered and thus spatially correlated (compare [49,50,51,52,53]). For the mixed-effect modelling we utilised height-diameter relationships based on power, Korf and Näslund functions, depending on the species group.

Tree level volume was calculated from tree height and diameter at breast height, based on species group-specific volume equations published by [54] for the 23 species groups applied in this project. Stem volume was converted to stem biomass by applying wood density coefficients documented for 41 species groups in the Master Plan for the Forestry Sector Nepal (Ministry of Forests and Soil Conservation of Nepal, 1989). Stem biomass was expanded to foliage and branches based on species group-specific equations for different diameter classes, taken from [54]. Above-ground tree biomass was calculated by summing up biomass of stem, foliage and branches.

Finally, the field plot data was screened for outliers. A plot is considered an outlier if it has plausible measurement or positioning errors, or if we suspect a disturbance of the area during the time between field campaign and LiDAR campaign. To detect outliers the field data was checked against LiDAR observations. If the difference between mean tree height calculated from field measurements and the 90th percentile of the first-pulse returns that was used as a surrogate for mean canopy height, was more than 10 m, the field plot was removed from the dataset. A rough estimate of AGB per plot was modelled from LiDAR data by regressing LiDAR-derived vegetation height and field-measured biomass. If the modelled biomass value differed more than 500 tonnes per hectare from the biomass that was estimated from field measurements, the plot was removed. The removal of outliers from the model parameter calibration step is essential in Bayesian modelling to regularize model formation: outliers will cause more harm than benefit to model construction. “Viable” outliers should be kept in the validation set, though. In total, 54 field plots were excluded from modelling, so that the remaining number of plots used in our study was 738. Statistics of the field data are shown in

Table 2. Statistics of the field data (712 plots).

Variable Name (Unit)	Min	Max	Mean	StD
Mean diameter weighted by basal area (cm)	5.9	127.9	34.2	17.0
Mean tree height weighted by basal area (m)	2.9	36.0	15.8	6.1
Basal area (m${}^2$/ha)	0.1	53.4	18.4	10.6
Number of trees (1/ha)	20	2219	679.3	450.1
Stem volume (m${}^3$/ha)	0.3	680.9	149.8	114.0
Total above-ground biomass (tons/ha)	0.4	829.1	189.1	142.6

.

2.5. LAMP2 with Stratification for Reference Level Generation

The Reference Level is calculated by multiplying Activity Data (forest area changes) with Emission Factors. Therefore it is essential to know (1) how much land changed within each forest type from one structural class to another in a given time period (Activity Data); and (2) how much carbon will be emitted when a forest class changes to another class (Emission Factors). To derive such Emission Factors, we need to know how much biomass and carbon is contained in each forest type and structural class. In the following we describe how Activity Data and Emission Factors were derived.

2.5.1. Satellite Data Acquisition and Processing

The best available Landsat 5 and Landsat 7 data, based on minimizing cloud cover from 1999, 2002, 2006, 2009 and 2011, were used as the raw data for generating activity data. ImgTools software (version 2.2) was used to conduct Spectral Mixture Analysis (SMA) of Landsat satellite imagery, Figure 3, into fractions with natural break points, known as endmembers [48]. SMA uses these endmembers to develop generic spectral libraries for green vegetation (GV), non-photosynthetic vegetation (NPV), and bare soil. The software has algorithms to generate water mask as well as shadow mask which are used to generate a normalized difference factional index (NDFI) and the shade-normalized green vegetation (GVs) [55,56].

A decision tree, built in the software was adjusted for the TAL based on the spectral curves of SMA components, to classify images into forest, non-forest, water bodies using fractional cover and GVs. The forest was further classified into intact and degraded forest using NDFI values. In order to avoid spectral confusion in areas previously deforested or degraded, this historical contextual information was used in combination with spectral curves to delineate areas of regrowth.

2.5.2. Image Processing

Image processing was done using different modules in ImgTools which are described below (Figure 3).

• Spectral Mixture Analysis (SMA): ImgTools was used to carry out spectral mixture analysis for each Landsat scene. The SMA module of ImgTools decomposes the spectral mixture, commonly found in the pixel reflectance values of remotely sensed data, into fractions with natural break points, known as endmembers. SMA module uses these endmembers to develop generic spectral libraries for green vegetation (GV), non-photosynthetic vegetation (NPV), bare soil and clouds [55,56].
• Water Mask: This module creates a water mask as a layer using fractional image.
• Cloud and Shade Mask: This module creates a cloud and shade mask layer that is used in deriving NDFI.
• Normalized Difference Factional Index (NDFI): In this module, the fractions developed from the SMA analysis: GV, NPV, Soil are processed to quantify the percentage of pixels lying outside the range of zero to 100% and to evaluate fraction value consistency for pixels over time (i.e., that pixels with intact forest values were similar over time). Only pixels with at least 98% of the values within zero to 100% and those that showed mean fraction value consistency over time were used by the software algorithm for computing Normalized Difference Fraction Index [55].

  $$\mathrm{NDFI}=\frac{{\mathrm{GV}}_{\mathrm{shade}}-(\mathrm{NPV}+\mathrm{Soil})}{{\mathrm{GV}}_{\mathrm{shade}}+(\mathrm{NPV}+\mathrm{Soil})}$$

  where ${\mathrm{GV}}_{\mathrm{shade}}$ (or $\mathrm{GVs}$) is the shade-normalized GV fraction given by ${\mathrm{GV}}_{\mathrm{shade}}=\frac{\mathrm{GV}}{100-\mathrm{shade}}$.

2.5.3. Image Classification

A decision tree to provide the supervised classification of forest structure built in ImgTools was adjusted for the TAL (Figure 4) based on the spectral curves of SMA components, to classify images into forest, non-forest, water bodies using fractional cover and GVs. The forest was further classified into intact and degraded forest using NDFI values. In order to avoid spectral confusion in areas previously deforested or degraded, this historical contextual information was used in combination with spectral curves to delineate areas of regrowth.

• Non-Forest—An area is classified as non-forest when it meets one of following criteria:

  • $53<\mathrm{GVs}<65$
  • $\mathrm{GVs}>65$ and $\mathrm{GV}>68$
  • $\mathrm{GVs}<52$ but $\mathrm{soil}+\mathrm{NPV}>14$
• Water: $\mathrm{GVs}<52$ but $\mathrm{soil}+\mathrm{NPV}<15$
• Forest: $\mathrm{GVs}\ge 66$ and $\mathrm{GV}<69$ (Justification here is forest will have shade from tall trees but the grassland will have virtually no shade)

  • Intact forest: $66<\mathrm{GV}<69$ and $\mathrm{NDFI}>168$
  • Degraded forest: $66<\mathrm{GV}<69$ and $\mathrm{NDFI}<168$
• Regeneration

  • Classified as intact forest in step 3 above and classified in previous time period as non-forest or degraded
  • Classified as degraded forest in step 3 above and classified in previous time period as deforested

The decision tree classification was then used to classify each satellite image into 5 classes: intact (undisturbed) forest, degraded forest, non-forest, water and cloud-shadow classes to produce a forest type map of TAL [48] with four major forest types: (1) Sal forest; (2) Sal dominated mixed forest (here after “Sal mixed forest”); (3) other than Sal dominated forest (here after “other mixed forest”) and (4) Riverine. These four forest types were overlaid on the forest structural map to generate forest type and condition maps for each time period. The study assumed forest types do not change from one type to another type (i.e., from Sal forest to mixed forest or riverine forest or vice versa) in 10–20 years.

2.5.4. Generation of Emission Factors Using Tier-2 LiDAR-Assisted Multi-Source Programme (LAMP2)

The Tier–2 version of LAMP is based on using stratification to increase accuracy of biomass estimates. The change in biomass density of forest per hectare over time is calculated as the difference in biomass content of three forest conditions: deforested, degraded and intact. For deforested area carbon density is always taken to be zero. Degraded and intact forest in each forest type—these are the strata of the LAMP Tier-2 stratification—are assumed to have a spatially and temporally constant biomass density. This constant biomass density for each stratum is estimated by random sampling thousands of grid cells per stratum in each forest type over biomass density map, generated in LiDAR blocks using Sparse Bayesian method. The aboveground carbon density values (tons/ha) are calculated by using carbon fraction 0.47 of AGB [57].

• LAMP2 step 1: Stratifying of forest on the study area using satellite data
  In the first step of the LAMP2 approach, the forest extent over the entire study area was stratified based on the forest types into Sal, Sal mixed, other mixed and riverine [48]. These strata were further divided into two conditions, intact and degraded, resulting in a total of eight forest classes.
• LAMP2 step 2: Estimating forest parameters for LiDAR blocks
  In the second step of the LAMP approach, a regression model was generated based on the relationship between LiDAR metrics (height and density distribution) and field based biomass data. It has been shown that Sparse Bayesian methods offer a flexible and robust tool for regressing LiDAR pulse histograms with forest parameters. While performing comparably to traditional regression methods, they are computationally more efficient and allow better flexibility than step-wise regression [7,58]. To correspond to the field plot size of 500 m${}^2$, the modelling of forest parameters was carried out at 22.4 m × 22.4 m grid-cell level. By using this grid size we also reduce the impact of potential Lidar-DEM errors introduced by steep terrain that will have a more pronounced effect on smaller grid cells. The Lidar metrics selected by the model for estimating above-ground biomass are described in [40]. The model was validated against an independent sample of 46 plots.
• LAMP2 step 3: Deriving forest class-specific mean biomass values
  In the third step of the LAMP2 approach, LiDAR model estimates are generated for a random sample of locations within the LiDAR blocks. These estimates are combined with the forest strata map to calculate mean biomass for each forest class. The procedure of this calculation is described in more detail in the paragraph below.

Within each of the eight forest type and condition classes 1000 circular LiDAR sample areas of 1 ha size were randomly allocated. To do so, the forest types and condition map of 2011 (LAMP2 step 1) was overlaid with the LiDAR data coverage. In this case a random sampling without weights was sufficient since the sampling was done separately for each forest class. If less than 50% of a LiDAR sample area belonged to same class as its center point, it was removed. If a LiDAR sample area had 50% or more of its area outside the forest mask or outside the LiDAR blocks, it was also removed. The remaining 7710 LiDAR sample areas were sub-divided into 500 m${}^2$ rectangular cells (estimation units) for the LiDAR model output. To allow forest class-specific biomass estimations, only those cells were used that were falling with their center into the same class as the original center point of the LiDAR sample area (intended forest class). The regression model based on LiDAR features (LAMP2 step 2) was applied to predict AGB for the cells. As the LiDAR sample areas contained a varying number of cells, the final results were aggregated as area-weighted mean for each forest class. Statistics of the results by forest class are shown in

Table 3. Statistics of the forest class- specific estimations for above-ground biomass, mean carbon density and carbon dioxide equivalent CO2e (t/ha).

Class	Nr of Plots	Above-Ground Biomass (t/ha)				C and CO2e Values
		Mean	Min	Max	StD	C (t/ha)	CO2e (t/ha)
1-Sal intact	988	235.6	20.4	509.5	84.1	110.7	406.0
2-Sal degraded	969	173.2	0.0	425.3	72.9	81.4	298.5
3-Salmix intact	966	183.2	0.0	556.9	84.7	86.1	315.7
4-Salmix degraded	946	146.4	0.0	539.6	106.2	68.8	252.3
5-Othermix intact	985	186.1	5.5	479.5	94.0	87.4	320.7
6-Othermix degraded	943	143.2	0.4	461.6	86.8	67.3	246.8
7-Riverine intact	934	171.1	0.0	405.5	46.8	80.4	294.9
8-Riverine degraded	979	99.4	0.0	505.6	57.9	46.7	171.3

. The mean biomass-per-hectare values calculated from LiDAR features were then applied to the corresponding forest classes to create a stratified biomass map.

The carbon stock is calculated as 47% of the AGB consistent with IPCC GPG (Chapter 4, Table 4.4) [39]. Therefore, the emission factors for each forest type and condition were calculated by multiplying the AGB by 0.47 (Table 3). When the forest changes from intact or degraded forest to deforestation all carbon was assumed to be released. But when forest goes from intact to degraded the difference in the mean carbon contents between intact and degraded forest is assumed to be emitted, for example when intact Sal forest changes to degraded Sal forest, 29.3 tC/ha or 107.5 tCO2/ha are emitted. The emission factors for regeneration forest changing to deforestation or degradation, and for sequestrations due to regeneration are calculated with the IPCC default value of 2.8 tC/ha/year or 10.3 tCO2/ha/year. Emission factors were derived by calculating the difference between the carbon and carbon dioxide equivalent, CO2e, values in Table 3 to reflect the loss of carbon or amount of emissions, calculated in tons of carbon dioxide CO2e when land area containing various forest types transitions from one structure to another.

2.5.5. Calculation of Emissions from Below-Ground Biomass

Based on IPCC GPG we used 20% of above-ground CO${}_2$ emissions as the below-ground emissions. Below-ground biomass was assumed to result in emissions at the time of mortality.

2.5.6. Time-Series Analysis of Satellite Data to Generate Activity Data

To delineate areas of deforestation, degradation and regeneration, we completed a time-series analysis of forest change for the TAL for four time periods, 1999–2002, 2002–2006, 2006–2009 and 2009–2011, using the classified images. A pair of classified images for the same satellite scene was run through a change detection algorithm in the ERDAS Imagine software (version 9.3) [55], to produce a change matrix at pixel level. This resulted in a 25–class matrix for the first set of image pairs (1999 and 2002). Any forested area under the cloud and cloud shadow (cloud-shadow class) was considered as unchanged within each period (1999–2002) for the purpose of this study. Likewise, areas remaining in same classes within each period were also considered unchanged. The change classes derived from the change matrix are listed below (

Table 4. New classes derived from the change matrix.

Change Matrix	Change Class
Intact forest to non-forest	Deforestation 1
Intact forest to degraded forest	Degradation
Degraded forest to non-forest	Deforestation 2
Non-forest to dense regenerating forest	Regeneration 1
Non-forest to sparse regenerating forest	Regeneration 2
Degraded forest to regenerating forest	Regeneration 3
Regeneration forest to non-forest	Deforestation 3

) as Deforestation 1–3, Degradation, and Regeneration 1–3.

For the subsequent time-series analysis the base classified image for that series (older of 2 images) was adjusted to reflect changes in the previous time period; for example, change classes derived in Table 4 as a change between 1999 and 2002 were delineated and re-coded in the 2002 scene, before comparing change between 2002 and 2006. All three types of deforestation were merged into one deforestation class because they represent areas going from forest to non-forest. Therefore, each base image potentially has nine classes: Intact Forest, Degraded Forest, Non-forest, Water, Cloud/Shadow, Deforestation, Regeneration 1, Regeneration 2, and Regeneration 3. The change analysis between 2002 and 2006 resulted in a 45-class change matrix with nine classes (described above) representing actual change in forest conditions. These nine change classes were adjusted in the base image (2006) for analyzing time series between 2006 and 2009. The same process was repeated for 2009 and 2011 series. The areas under each activity (Deforestation 1–3, Degradation, and Regeneration 1–3) for each time series analysis were used to generate activity data (

Table 5. Activity data for different forest types between 1999 and 2011.

Forest Type	Transition	Activity	Activity Data (ha)
			1999–2002	2002–2006	2006–2009	2009–2011	12-Year Total
Sal Forest	Intact to Deforested	Deforestation 1	11,583	2085	9488	17,914	41,070
	Degraded to Deforested	Deforestation 2	4322	679	615	1651	7268
	Regenerated to Deforested	Deforestation 3		905	2117	6655	9677
	Intact to Degraded	Degradation	10,831	1342	3141	17,488	32,803
	Deforested to regrowth	Regeneration	24,635	35,951	6313	10,008	76,907
Sal Mixed	Intact to Deforested	Deforestation 1	8487	2291	10,588	20,332	41,697
	Degraded to Deforested	Deforestation 2	7632	1395	964	1927	11,918
	Regenerated to Deforested	Deforestation 3		1996	3405	12,821	18,222
	Intact to Degraded	Degradation	10,186	1661	10,003	10,375	32,225
	Deforested to regrowth	Regeneration	32,597	40,999	4995	11,886	90,477
Other Mixed	Intact to Deforested	Deforestation 1	2029	273	2661	3308	8271
	Degraded to Deforested	Deforestation 2	674	175	514	284	1647
	Regenerated to Deforested	Deforestation 3		174	870	1536	2580
	Intact to Degraded	Degradation	1570	216	380	1250	3417
	Deforested to regrowth	Regeneration	2483	5239	1251	3461	12,434
Riverine	Intact to Deforested	Deforestation 1	918	160	255	1663	2995
	Degraded to Deforested	Deforestation 2	458	59	39	163	719
	Regenerated to Deforested	Deforestation 3		76	147	752	974
	Intact to Degraded	Degradation	697	81	225	877	1881
	Deforested to regrowth	Regeneration	2202	3306	510	244	6262

). Activities Regeneration 1–3 were combined to a single Regeneration activity because all these activities were differentiated only based on activities in the previous time period that resulted in regeneration in the current period, thus their growth rates and mean carbon content are assume to be same.

2.5.7. Generating Reference Level (RL)

The RL is generated by multiplying areas changed under each activity by the appropriate emission factor, i.e., mean carbon stocks in each forest type to calculate amount of CO${}_2$ emission due to that particular activity.

$$\mathrm{RL}=\mathrm{Activity}\mathrm{data}\times \mathrm{Emission}\mathrm{factors}$$

(1)

The amount of CO${}_2$ released due to loss of forest carbon resulting from deforestation and degradation is termed as gross emissions while intake of CO${}_2$ by growing plants during forest regeneration is called sequestration. Therefore, net carbon loss is equal to gross emissions minus sequestrations. The reference emissions level (RL) for TAL is based on net carbon accounting process.

2.5.8. Calculating Net Emissions Level

Following formula was used to calculate RL for TAL:

$$\mathrm{Reference}\mathrm{level}=\frac{\sum {\mathrm{Em}}_{\mathrm{def}}1+\sum {\mathrm{Em}}_{\mathrm{def}}2+\sum {\mathrm{Em}}_{\mathrm{def}}3+\sum {\mathrm{Em}}_{\mathrm{deg}}-\sum {\mathrm{Seq}}_{\mathrm{reg}}}{y},$$

(2)

where $\sum {\mathrm{Em}}_{\mathrm{def}}1$ is the sum of emissions from deforestation of intact forest over y years, $\sum {\mathrm{Em}}_{\mathrm{def}}2$ is the sum of emissions from deforestation of degraded forest over y years, $\sum {\mathrm{Em}}_{\mathrm{def}}3$ is the sum of emissions from deforestation of regenerated forest over y years, $\sum {\mathrm{Em}}_{\mathrm{deg}}$ is the sum of emissions from degradation over y years, and $\sum {\mathrm{Seq}}_{\mathrm{reg}}$ is the sum of sequestrations from regeneration over y years.

2.6. LAMP3 with Estimation of Above-Ground Biomass at 1 ha-Scale

The LAMP3 method is based on a two-stage regression modelling by incorporating field data (sample plots), LiDAR data (sample blocks) and full-coverage satellite imagery. It results in estimates of AGB on 1-hectare resolution.

In addition to the collected field plot data and LiDAR sample data that have been described in Section 2.3, five medium-resolution Landsat 5 TM scenes of processing Level 1T from years 2009/2010 were acquired to cover the entire study area. The data were ortho-rectified and corrected for atmospheric and radiometric effects. After masking out cloud and snow areas, pixel values from overlapping areas between adjacent images were extracted. To obtain a homogeneous mosaic data set, the five Landsat scenes were spectrally normalised relative to each other.

2.6.1. Variance-Preserving Landsat Image Mosaicking

Due to variation in acquisition conditions (e.g., solar illumination, atmospheric scattering and atmospheric absorption) the same ground object on two overlapping images can result in different spectral values [59]. Because of this, radiometrically uniform mosaics using multitemporal scenes should be created before employing satellite imagery into carbon assessment. To overcome radiometric differences between multitemporal scenes, relative normalization is performed under an assumption of a linear relationship between overlapping regions of multitemporal images [60,61,62,63,64,65], where the linear transformation of a pixel intensity for image $P_k$ is defined by two coefficients: shift $a_k$ and scale factor $b_k$, $k=1,...,N$, here N is the number of images in the set. In [66], we proposed a variance-preserving mosaic (VPM) algorithm that minimizes overall error of the normalization and aims to preserve average variance of normalized images. We introduced the corresponding additive penalty function and used the Limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm to find solution to the optimization problem. Then, to avoid down-scaling of the output image intensity with $b_k<1$, we calculate the final scale factors for all images: $b_k^{*}=b_k/min{b}_k\ge 1$, and recalculate the corresponding shift coefficients.

With the normalisation algorithm in use, the differences between overlapping pixels is minimized while at the same time the original variance in the images is preserved. For each image, optimised normalisation coefficients (shift and scale factor) were found and applied to adjust the image data (

Table 6. Normalisation coefficients a (shift) and b (scale factor) applied to the Landsat scenes during mosaicking.

Path / Row	Band	a	b
141/41	1	$-561.12$	2.12
	2	$-365.57$	1.62
	3	$-135.51$	1.28
	4	0	0.84
	5	$-140.97$	1.02
	7	$-180.50$	1.01
142/41	1	$-276.10$	1
	2	$-123.43$	1
	3	$-73.20$	1
	4	$-170.56$	1
	5	$-79.19$	1
	7	$-98.55$	1
143/40	1	$-1146.87$	2.16
	2	$-618.48$	1.55
	3	$-326.08$	1.23
	4	$-709.97$	1.43
	5	$-258.92$	1.24
	7	$-154.35$	1.17
143/41	1	$-267.44$	1.26
	2	$-111.18$	1.18
	3	$-23.62$	1.12
	4	$-152.39$	1.06
	5	$-178.36$	1.15
	7	$-79.64$	1.02
144/40	1	$-603.32$	1.44
	2	$-383.87$	1.28
	3	$-229.26$	1.16
	4	$-387.02$	1.22
	5	$-1.37$	1.08
	7	23.84	1.01

). The normalisation was performed separately for each spectral band.

2.6.2. Applying the LiDAR Model to Calculate AGB Estimates on Surrogate Plots

In the first stage of modelling, the LiDAR-based regression model as described in Section 2.5.4 was applied to predict forest characteristics for a set of 10,000 circular-shaped surrogate plots (simulated field plots) of 1-hectare size within the forested area of the LiDAR-scanned blocks. The locations of the surrogate plots were selected through weighted random sampling using the inverse of the block weights applied in LiDAR block sampling (see Section 2.5.4). The chosen surrogate plot size of one hectare is on the one hand large enough to decrease edge effects and to compensate for geometrical discrepancies between LiDAR and satellite data, and on the other hand it is consistent with the desired output resolution of the final LAMP3 biomass estimates. To predict forest characteristics for the surrogate plots, each 1-hectare surrogate plot was first divided using a grid with cell size of 500 m${}^2$. Estimates were calculated for each grid cell within the surrogate plot, so that the estimation area unit was the same as used in model formation to prevent a possible bias due to the scale-dependent LiDAR features. Finally, grid cell estimates were aggregated over the 1-hectare surrogate plot area.

2.6.3. LAMP3 Model Construction

In the second modelling stage, the AGB values that were estimated for the surrogate plots from LiDAR data are applied as simulated ground-truth to generate a regression model between above-ground biomass and features derived from the Landsat mosaic in order to derive biomass values for all forests in the project area. LiDAR data acquired through airborne laser scanning over forested areas provides a sufficiently accurate reference when estimating AGB from satellite imagery. Therefore, LiDAR-derived estimates of AGB can be used as a substitute for real field measurements because they represent each average AGB values over a much larger area than that of a traditional field plot. Average AGB values over 1 ha are more stable than those over 500 m${}^2$. Unlike field plots, surrogate plots also include forest on difficult and inaccessible terrain in mountainous regions while still reaching an accuracy that is comparable or close to those of field measurements (described in Section 2.5.1 and in Section 3.4.4).

For each surrogate plot area, the spectral mean values were extracted from Landsat 5 TM imagery acquired for a similar time as the LiDAR data. In addition to the values in the visible and infrared image bands, two vegetation indices were calculated from the Landsat data to facilitate the biomass modelling, namely the normalized difference vegetation index (NDVI) and the atmospherically resistant vegetation index (ARVI). The ARVI is an enhancement of the NDVI that is relatively resistant to atmospheric factors like aerosol. It uses reflectance in the blue band to correct red reflectance for atmospheric scattering [67]. Sparse-Bayesian method was used to regress the satellite-derived variables with the LiDAR-derived forest characteristics for the locations of the surrogate plots. The regression model was validated using k-fold cross-validation against the LiDAR-based estimates on all 1-hectare surrogate plots. The model was applied to estimate basal area, volume and above-ground biomass from the satellite-based variables for the entire study area using an output grid of 100 m × 100 m (1 ha) cell size. Non-forested areas were clipped off from the resulting map using the forest mask published by the FRA Nepal project [68].

2.6.4. Variance-Preserving Histogram Matching

Optical remote-sensing does not capture information from below dense canopy covers as is the case for sub-tropical forests in TAL. This results in a “saturation effect” present in the biomass estimates which leads to underestimation of biomass in areas with high AGB concentrations [69,70]. This problem does not concern LiDAR which can penetrate even through a closed canopy cover to return information from ground-level. To compensate for the saturation effect, the satellite-based biomass estimates were post-processed by applying a histogram matching method using the biomass distribution in the LiDAR-based surrogate plots as a reference. A similar approach of calibrating satellite-derived predictions to reference data has been used e.g., by [25] for the creation of a pan-tropical biomass map.

In our case, the reference biomass values were sorted by their value and binned into equally sized quantiles. 1000 of such quantiles consisting each of 10 surrogate plot values have been used. For each quantile, the mean biomass value was calculated from the surrogate plots of that quantile. In the same way also the predicted biomass values of the Landsat-based model applied to the entire study area were binned into quantiles, using the same number of quantiles as for the reference data. The mean value of each quantile of the estimated data was replaced by the mean value of the corresponding quantile of the reference data. This way the value range of the satellite-based estimates was matched with the value range of the reference data (surrogate plots), preserving the original standard deviation in the reference data and improving the regression trend line.

By applying this method we assume that the surrogate plots are a representative and unbiased sample of the biomass distribution over the entire TAL area. We base our assumption on the fact that the locations of surrogates were indeed generated as a random sample by applying the inverse weights of those used in the LiDAR block sampling (compare Section 2.6.2). However, it must be acknowledged that a deviation of the sampled distribution from the distribution of the entire population may lead to a bias in the final biomass results.

3. Results

3.1. Reference Emissions Level (RL) Estimation

The RL analysis shows that during the 12-year period between 1999 and 2011 total of 52,245,991 tons CO${}_2$ (tCO2e) was emitted from the forest sector in the TAL, an average emission of 4,353,833 tons CO2e per year (

Table 7. Forest-related CO${}_2$ emissions in Terai Arc Landscape (TAL) between 1999 and 2011.

Period	CO${}_2$ Emissions (tCO2e)		Total
	Above-Ground	Below-Ground
1999–2002	13,136,430	2,627,286	15,763,716
2002–2006	1,736,537	347,307	2,083,845
2006–2009	9,644,698	1,928,940	11,573,637
2009–2011	19,020,661	3,804,132	22,824,793
Total 12-year	43,538,325	8,707,665	52,245,991
Average annual	3,628,193.79	725,639	4,353,833

). In the period 2006–2011, emissions averaged 6,879,686 tCO2e per year, 58% higher than the 12–year average, and in the period 2009–2011, emissions averaged 11,412,396 tCO2e per year or 162% higher than the 12–year average.

3.2. Reference Level at District Level

TAL falls under 12 districts or administrative units so district-level analysis was conducted to better understand geographic trends. District-level RL analysis is presented in

Table 8. Total CO${}_2$ emission (tCO2e) by districts for 4 time intervals.

	1999–2002	2002–2006	2006–2009	2009–2011	12-Year Emissions
Kahchanpur	1,326,570	120,105	296,008	3,499,486	5,242,169
Kailali	3,736,460	93,151	911,511	7,891,560	12,632,682
Bardia	425,756	151,066	312,516	3,116,150	4,005,488
Banke	1,227,909	304,491	2,515,125	567,689	4,615,215
Dang	2,600,210	582,332	4,759,420	892,183	8,834,146
Kapilbastu	1,594,386	113,716	1,025,029	380,993	3,114,124
Rupandehi	597,963	(24,121)	72,593	224,251	870,686
Nawalparasi	1,869,896	171,651	758,771	456,103	3,256,421
Chitwan	1,388,989	267,881	250,988	1,315,372	3,223,230
Parsa	189,225	76,152	142,864	872,272	1,280,513
Bara	395,579	96,825	207,383	1,615,801	2,315,588
Rautahat	410,772	130,596	321,429	1,992,933	2,855,730

. In addition to the significant differences in rates of deforestation and degradation for the various time intervals, there are also significant geographic variations in the distribution of forest-related emissions. Three of the 12 districts—Kailali, Kachnapur and Dang—accounted for 51% of the carbon loss of the TAL during the RL period.

3.3. High-Resolution AGB Maps Calculated in TAL with LAMP3

Estimation of AGB Change with LAMP3

The biomass estimates predicted at 1-hectare with the LAMP3 method were compared against the highly accurate LiDAR predictions for the locations of the surrogate plots described in Section 2.6.2. The initial predictions of total AGB significantly underestimate areas rich in biomass (saturation effect) and overestimate areas of low biomass. The histogram matching step (compare Section 2.6.4) overcomes this problem but causes an increase of the root mean square error (RMSE).

After histogram matching, the final estimates of total AGB with LAMP3 achieve an RMSE of 39%. Mean tree height and basal area are the most accurately estimated variable with 26% RMSE. For all predicted variables the mean values and standard deviations of the estimates match with the mean values and standard deviations of the reference data. The validation statistics before and after histogram matching are shown in

Table 9. Error statistics for the LAMP3 estimates before histogram matching, when compared to LiDAR predictions, on 9805 surrogate plots of 1-hectare size each.

Variable	Estimates		Surrogate Plots		RMSE	RMSE (%)	Bias	Bias (%)
	Mean	Std	Mean	Std
AGB, t/ha	198.2	53.1	198.1	88.1	69.8	35.2	0.12	0.06
Volume, m${}^3$/ha	157.1	42.5	157.0	71.0	56.5	36.0	0.10	0.06
Basal area, m${}^2$/ha	19.2	4.1	19.2	6.1	4.6	23.8	0.05	0.00
Diameter, cm	35.0	4.3	35.0	10.2	9.3	26.6	0.02	0.00
Height, cm	16.2	2.1	16.2	4.2	3.7	22.6	0.00	0.00

and

Table 10. Error statistics for the LAMP3 estimates after histogram matching, when compared to LiDAR predictions, on 9805 surrogate plots of 1-hectare size each.

Variable	Estimates		Surrogate Plots		RMSE	RMSE (%)	Bias	Bias (%)
	Mean	StD	Mean	StD
AGB, t/ha	198.0	88.2	198.1	88.1	77.5	39.1	$-0.04$	$-0.02$
Volume, m${}^3$/ha	157.0	71.1	157.0	71.0	63.0	40.1	$-0.03$	$-0.02$
Basal area, m${}^2$/ha	19.2	6.1	19.2	6.1	5.0	25.9	$-0.00$	$-0.00$
Diameter, cm	35.0	10.2	35.0	10.2	11.2	31.9	$-0.00$	$-0.00$
Height, cm	16.2	4.2	16.2	4.2	4.3	26.4	$-0.00$	$-0.00$

. Figure 5 illustrates the comparison between the LAMP3 results and LiDAR predictions, before and after histogram matching.

To compare the results obtained with LAMP2 and LAMP3, a high resolution AGB map of TAL was produced with LAMP3. In the Figure 6 below, a part of a district in Terai is displayed, illustrating the high spatial resolution obtainable with direct AGB estimates with LAMP3. The areas identified with both methods as deforested or degraded coincide reasonably well by visual analysis, but a comprehensive comparative analysis has not yet been possible to carry out. It can be seen from the images that inhomogeneity of satellite images causes some artefacts to appear into the biomass maps, but the level of these artificial features stays within the standard deviation of the model based estimator. Further research into the calibration of satellite images will be conducted in order to reduce the impact of image inhomogeneity.

3.4. Uncertainty Assessment

Uncertainty assessment of forest inventory often has two somewhat divergent goals. On the one hand, we wish to bind population totals of large forest areas between tight confidence intervals. On the other hand, medium and high resolution forest parameter maps, such as biomass or carbon maps, should also be as accurate and precise as possible. The latter purpose serves the needs of operational planning of forest interventions and in that case large-scale forest statistics are often not very relevant.

LAMP methods aspire to serve both purposes as well as possible, although some degree of trade-off is still necessary. The principal goal of LAMP processes is to produce accurate and precise high-resolution forest biomass maps. But a secondary goal is to construct that map in a manner that allows high-resolution results to be aggregated to large area averages without engendering significant systematic error.

There are two main principles by which this is achieved. Firstly, a field campaign to collect calibration plots for the model-based estimation must be design-based [3]. Secondly, the regression method applied in model building must be an unbiased estimator, such as a linear regression model. An underlying assumption in this latter process is that a linear correlation exists between biomass and the co-variates used to estimate it in the model. The first condition is fulfilled for both LAMP model trials in TAL, where as the latter condition is fulfilled for an individual LAMP3 model, such as the one built for AGB estimation in TAL, but not for the methodology as a whole, because covariates are picked up and dropped by the Sparse Bayesian algorithm in a non-linear fashion, and also because it often turns out that the relationship between most, if not all, satellite-based covariates and AGB is quite non-linear.

To address the fullest range of these divergent goals, the uncertainty assessment of the biomass estimates described in the following subsections has been conducted from several different angles. Firstly, classical variance estimation was carried out on the two-level regression of LAMP3 to assess the impact of model residuals. Secondly, several different field plot assessments of both LAMP estimates were conducted that feature analysing the impact of field plot size on error variance and the accuracy of the LiDAR model at high spatial resolution, so as to validate the usability of sample areas for Emission Factor estimation in LAMP2 and that of surrogate plots as a teaching set in LAMP3. Furthermore, a separate accuracy assessment was conducted on Activity Data in LAMP2 using a confusion matrix and validation on high-resolution satellite imagery. And finally, a cross-validation analysis of both the LiDAR model and the LAMP models was conducted against two different validation sets: the calibration plot set of 738 field plots and a set of ten thousand surrogate plots interpreted through the LiDAR model.

3.4.1. Variance Estimation of Two-Level Regression Models

Recently, the authors in [36] have produced variance estimators for homogeneous and inhomogeneous two-layer estimation processes based on linear models. Since the models used in the current LAMP3 are linear, a study was conducted to apply the methods introduced in [36] to these estimates. Following this reference, the tables included in Appendix B indicates the Standard Deviations at different scales based on a Monte Carlo study.

A heteroscedastic model construction was assumed and covariance matrices of both the LiDAR model and the satellte model were constructed. Based on these matrices, a variance estimator for the model predictions was constructed and applied at several different forest biomass map resolutions, ranging from 500 ha to 10,000 ha. The resulting absolute and relative standard deviations are presented in

Table A1. Average estimated surrogate plot aboveground biomass (AGB-LAMP), and variance and standard deviation (SD) of OLSLAMP and LAMP3 methods calculated using four different sets of field and surrogate plots. Field plots (m): 50, 100, 500. Surrogate plots (M): 500, 1000, 5000. The last set has all field plots (738) and surrogate plots (9805). All sets, excluding the full set, have been ran 100 times and an average of sample results have been added to table.

		OLSLAMP	OLSLAMP	LAMP3	OLSLAMP	LAMP3	OLSLAMP	LAMP3
Data Size	Area, ha	Estimates,	Variance,	Variance,	SD,	SD,	Relative SD,	Relative SD
		tons/ha	(tons/ha)${}^2$	(tons/ha)${}^2$	tons/ha	tons/ha	%	%
$m=50$	2.5 ha
$M=500$	500 ha	198.2	32972.2	33,985.4	181.6	184.4	91.6	93.0
$m=100$	5 ha
$M=1000$	1000 ha	198.2	6211.1	6453.1	78.8	80.3	39.8	40.5
$m=500$	25 ha
$M=5000$	5000 ha	198.2	260.9	280.4	16.2	16.7	8.2	8.5
$m=738$	36.9 ha
$M=9805$	9805 ha	198.2	141.7	149.5	11.9	12.2	6.0	6.2

.

3.4.2. Validation of Activity Data through Additional Field Verification

A weighted random stratified sampling design was used to select 200 field plots of 1-ha (100 m × 100 m) covering intact (no change), deforested, degraded, and regenerated areas based on time series analysis. The goal was to cover about 5% area of each change category. However, after field visits, the team found 110 (>50%) sampling plots were concentrated in 4 eastern districts. Therefore, to maintain consistent sampling across the study area, only 50 plots were selected randomly from these 4 district resulting in 140 potential sampling plots across the TAL. Among these plots, the field team was able to measure only 103 plots, other plots were inaccessible. Using GPS, the field crew navigated to the center of 1-ha plot to collect information on forest condition types (intact, degraded, deforested and regeneration). The field crew also estimated crown closure, ground cover based on visual observation. A relascope was used to estimate basal areas of tress in each 1-ha plot. The plots were categorized based on these information as intact, degraded, deforested and regeneration. These plots were then overlaid over the forest change map from time series analysis between 2009 and 2011, to generate an error-matrix and 95% confidence interval for accuracy assessment. The field plots were overlaid on the changed map resulted from a time series analysis of 2011 and 2014 satellite data, for an accuracy assessment of the activity data. The activity data viz. intact, deforested, degraded, and regenerations derived from the time series data and data observed/measured in the field were tallied. Tallied numbers were then multiplied by the proportions of area in each activity, based on the change map derived from time series analysis of 2011 and 2014 data, to generate an error-matrix and 95% confidence intervals for each activity (Table 4), following the process used by [71]. Overall accuracy of activity data was 85% ± 14% at 95% confidence interval. The producer’s accuracy and user’s accuracy for each activity with 95% confidence intervals are presented in the

Table 11. An error matrix showing accuracy of forest change between 2009 and 2011 with 95% confidence intervals.

Activity	Intact	Deforestation	Degradation	Regeneration	Total	Mapped Area	Proportion Wi
						(ha)	(ha)
Intact	0.704	0.016	0.008	0.142	0.871	858,910	0.871
Deforestation	0.008	0.063	0.001	0.002	0.074	72,700	0.074
Degraded	0.003	0.005	0.024	0.000	0.032	31,398	0.032
Regeneration	0.001	0.003	0.001	0.020	0.024	23,623	0.024
Total	0.716	0.086	0.034	0.164	1.000	986,631	1.000
Overall accuracy	0.81 ± 0.09
Producer’s accuracy	0.98 ± 0.065	0.73 ± 0.024	0.72 ± 0.017	0.87 ± 0.061
User’s	0.81 ± 0.092	0.86 ± 0.007	0.76 ± 0.009	0.82 ± 0.004

.

3.4.3. Impact of Field Plot Size

The impact of validation field plot size on the reported root mean square error (RMSE) in forest inventory has been studied in detail [72]. In general, larger field plot size yields smaller RMSE when the estimates do not contain significant systematic error. In general, field validation exercises should be conducted with field plots that have a size that matches the size of the primary estimation unit. If the latter is rather large, as is the case in the current study with 1 ha unit size, comprehensive field campaigns that calculate and measure every tree on a statistically adequate number of large field plots are not feasible. We therefore must resort to an analysis of the impact of field plot size to the precision of forest parameter estimates.

To study the relative impact of both variance between plots and that of the edge effect on different plot sizes, additional field plots were collected from LiDAR-covered area in Terai Arc Landscape (TAL), Nepal, during March 2013. The main objective was to analyse the impact of field plot size on LiDAR model accuracy. Field measurements on larger sized plots of 30 m radius were carried out within the area of two LiDAR blocks, in order to compare results and accuracy between LiDAR predictions derived from field plots of different size (small plots versus large plots). The aim was to see if it is possible to reduce costs by collecting few larger plots instead of many smaller plots. In earlier studies [73,74,75] researchers have demonstrated that the larger plot minimizes edge effects, increases sample variance, and maintains a greater amount of spatial overlap between ground-reference and LiDAR. The effect of the GPS error on model accuracy is smaller in larger plots because of relatively larger overlap. Two LiDAR blocks were selected for extra field plot collection: one representing typical Terai Sal and associated forests, and another one representing Siwaliks Sal and dry deciduous forest types. The plot design was based on weighted random sampling using LiDAR canopy-height information from existing LiDAR data to capture the full heterogeneity of the forest. Field measurements were taken from 48 plots, while 2 plots were inaccessible. For 38 plots, a fixed circular plot with outer radius of 30 m was used, equivalent to an area of 2826 m${}^2$. For 10 plots, an outer radius of 40 m was used, equivalent to an area of 5026.5 m${}^2$. For every tree on these large plots, the distance to the plot centre was also recorded, so as to facilitate simulations of validation results also on smaller circular plots by restricting the tree lists by this distance.

The overall AGB variance within the plot sample strongly decreases with increasing plot size from 5 m to ∼10 m. However, it was found that beyond a plot radius of ∼15 m, increasing plot size does not decrease between-plot variance of AGB anymore, see Figure 7a.

To study the edge effect, the amount of edge trees within a plot was calculated for varying plot sizes. An edge tree was defined as a tree which is within 1 m distance from plot edge on either side. It was found that the portion of edge trees in relation to all trees on the plot greatly decreases as the plot size increases. Beyond a plot radius of ∼15 m, increasing plot size does not significantly decrease the edge effect anymore, see Figure 7b. After the change in edge effect stabilizes, the difference between the plot sizes is more or less random in nature.

3.4.4. LiDAR Model Errors on Different Plot Sizes

For each plot set (5 m plots, 6 m plots, …, 40 m plots), a linear regression model was created between AGB and LiDAR variables for the plot locations. The Sparse Bayesian method was applied to find the best predictors out of the 46 LiDAR variables to estimate biomass, for each plot set separately. This means, for each plot set a different linear regression model was used to estimate biomass. Leave-one-out cross validation technique was then used for predicting biomass for each plot within a plot set. The biomass estimates were validated against 30 m plot data, for all plot sets. This was done in order to make the validation comparable between the different plots sets: When validated against the same field data, differences between plot sets resulting from measurement errors could be excluded.

The results show that model errors and bias significantly decrease up to a radius of about 15 m, while there are only small improvements beyond that (Figure 8a–c). The Coefficient of determination improves slightly longer, but does also reach its apex after some 20 m (Figure 8d). Since there have been only 10 plots available with radius above 30 m, the model errors for those plot sets are higher because 10 plots are not sufficient to create a good linear regression model.

3.4.5. Validation of Results by a Separate Field Campaign

The LiDAR model was calibrated with 738 field plots of 12.6-m radius (500 m${}^2$) and validated against plots of 30-m radius. To calculate biomass estimates for the 30-m plots, each validation plot (with area of 2826 m${}^2$) was split into 500 m${}^2$ estimation units (cells). This was done because the optimal estimation area for a model is equal to the size of the calibration plots, i.e., in this case 500 m${}^2$. For each estimation unit, a separate estimate was produced by the model. Finally, an area-weighted average of all the biomass estimates within a validation plot was derived from the cell-level estimates of that plot. The estimates of the 30-m plots were compared to the field based values, and error statistics were calculated. The relative RMSE was 17%, and the achieved R2 0.92. No significant bias was present (relative bias 1.3%). Validation results are shown in

Table 12. Error statistics for the LiDAR model built from 738,12.6 m-plots. The model predictions of AGB were validated against field data from 30 m plots.

	Estimates		Ref. Plots		Error
Variable	Mean	Std	Mean	Std	RMSE	Rel. RMSE (%)	Bias	Rel. Bias (%)
AGB in Tonnes/hectare	182.8	104.2	180.4	108.5	30.8	17.1	2.4	1.3

and Figure 9. As comparison, the figure includes also the validation results for the model with the original calibration field plots using leave-one-out cross-validation. This comparison demonstrates the error dependency on scale. The validation results directly depend on the validation plot size: The LiDAR model achieves significantly better validation results when being validated against the larger (30 m) plots, even though the model itself is the same.

4. Discussion and Conclusions

In order for MRV methods to succeed in mitigating climate change, they have to be designed so that results are efficient, effective and verifiable. Efficiency means that monitoring can be conducted frequently, ideally on a biennial basis, i.e., every second year, to match the agreed reporting frequency of the Paris Climate Agreement. On the other hand efficiency cannot sacrifice verifiability. MRV processes must produce reliable data on carbon captured in forests that donors can trust. Therefore MRV results must be accurate and contain no systematic errors. Finally, in order for REDD+ to be also effective, it must reward mitigation efforts that have a true impact. Often this means addressing sustainable forest use practices at a very local level, so that rewards can be distributed on a small scale to benefit the people that in practice tend or log in forests.

To achieve these goals, MRV cannot constantly depend solely on comprehensive field campaigns across a country or a jurisdiction. Such large-scale measurement efforts can, when well designed, be accurate and deliver on verifiability ([1,2,41,76,77] and references therein), but will fall short either on effectiveness, if plots are laid out too sparsely, or affordability, if plots are assigned to, say, every hectare. Methods based purely on remote sensing will, on the other hand, fall short on accuracy since without ground truth there is no way to guarantee that they contain no systematic error.

A combination of field measurements and remote sensing based MRV processes is therefore called for for simultaneously achieving the triple goals of accuracy, precision and affordability. Different remote sensing modalities have different properties. Optical satellite images provide frequent and comprehensive coverage of almost any forest area, but there are many difficult issues in associating band values with biomass. Airborne Laser Scanning (ALS or LiDAR) on the other hand provides high quality height information of trees but is expensive to collect over large areas.

The results obtained with LiDAR-Assisted Multi-source Program are encouraging. Nepal was accredited with 14 million tons’ worth of CO${}_2$ emission rights by the Forest Carbon Partnership Facility (FCPF) of the World Bank, based on the Reference Levels at district level calculated with LAMP2. These calculations were incorporated in the required Emission Reduction Project Idea Note (ERPIN) in 2014. Currently the Government of Nepal is pursuing the next stage of that process and producing an Emission Reductions Program Document, ERPD. This work will involve projecting future emissions and their reduction on the same area: the Terai Arc Landscape.

As seen from the additional field validation of Activity Data in Table 12, there is some amount, roughly 85%, of mis-classification between forest classes. These can be caused e.g., by optical effects, such as mountain shadows, or other differences caused e.g., by different times of day, weather conditions and seasonal variation. Mis-classification will possibly result in artificial carbon dioxide emissions and/or carbon sinks when they occur between different assessment time intervals. Such optical artefacts are common in all land use class based biomass assessment approaches and they are a prime motivator to introducing LiDAR wall-to-wall surveys or LiDAR sampling, as in the LAMP3 method, into REDD+ MRV processes.

As explicitly calculated and reported in [35], the source of biomass estimation errors across many spatial scales in LiDAR based forest estimation changes as a function of scale. At the smallest, individual tree scale allometric model errors are prevalent, as the authors found out from a comprehensive destructively tested tree database. But already at typical field plot levels, LiDAR model residuals become the dominant source of estimation error. According to these authors, the next change in error domination occurs at roughly 380 m resolution that corresponds to 15 ha resolution, after which LiDAR model parameter estimation errors become dominant. Their model uses only two LiDAR metrics as covariates whereas the ones adopted in the current research typically have between ten and twenty LiDAR or satellite covariates in the model and a relatively comprehensive calibration set of 738 plots of 500 m${}^2$ size or 9908 surrogate plots of 1 ha size. The exact error budget might therefore shift from the conclusions of [35], but its overall pattern will very likely remain quite similar.

The process and cost of acquiring the necessary LiDAR scans and the corresponding field plots has been analysed in [78]. In that study, LAMP was seen as a viable alternative to traditional comprehensive field campaigns and to possess the additional potential of producing high-resolution carbon density maps also in the follow-up stages. In that study it was estimated that the total cost of LAMP will remain lower than continuous monitoring of field plots already from the second monitoring stage onwards.

While LAMP2 has been vindicated so far by a donor process, LAMP3 still needs further development steps. Inhomogeneous satellite imagery remains a formidable issue for direct estimation of AGB. In that respect, new remote sensing modalities provide reason for hope. Radar interferometry is capable of providing accurate height information. Also new developments in LiDAR technology, such as so-called single-photon or Geiger-mode LiDARs, have the potential of dramatically reducing LiDAR acquisition costs over large areas in the near future [79].

As can be seen from Table A1, the classical variance estimator produces a higher standard deviation of estimation error on small estimation units when compared to a LOO validation, such as the one shown in Table 10. This is to be expected, since the posterior error variance in a Bayesian estimator is typically smaller than in the prior estimator [6], and the variance estimators presented in [36] correspond to prior estimators since they assume additive variance accumulation in both stages of the estimation process which ignores the fact that both models have been calibrated based on the same field plot set. But if the field plot set were not a probabilistic sample of the forest area, then LOO variance estimates could be much too low.