Reducing Uncertainty in Mapping of Mangrove Aboveground Biomass Using Airborne Discrete Return Lidar Data

: Remote sensing techniques offer useful tools for estimating forest biomass to large extent, thereby contributing to the monitoring of land use and landcover dynamics and the effectiveness of environmental policies. The main goal of this study was to investigate the potential use of discrete return light detection and ranging (lidar) data to produce accurate aboveground biomass (AGB) maps of mangrove forests. AGB was estimated in 34 small plots scatted over a 50 km 2 mangrove forest in Rio de Janeiro, Brazil. Plot AGB was computed using either species-speciﬁc or non-species-speciﬁc allometric models. A total of 26 descriptive lidar metrics were extracted from the normalized height of the lidar point cloud data, and various model forms (random forest and partial least squares regression with backward selection of predictors (Auto-PLS)) were tested to predict the recorded AGB. The models developed using species-speciﬁc allometric models were distinctly more accurate (R 2 (calibration) = 0.89, R 2 (validation) = 0.80, root-mean-square error (RMSE, calibration) = 11.20 t · ha − 1 , and RMSE(validation) = 14.80 t · ha − 1 ). The use of non-species-speciﬁc allometric models yielded large errors on a landscape scale (+14% or − 18% bias depending on the allometry considered), indicating that using poor quality training data not only results in low precision but inaccuracy at all scales. It was concluded that under suitable sampling pattern and provided that accurate ﬁeld data are used, discrete return lidar can accurately estimate and map the AGB in mangrove forests. Conversely this study underlines the potential bias affecting the estimates of AGB in other forested landscapes where only non-species-speciﬁc allometric equations are available.


Introduction
Forests play a crucial role in the global carbon cycle (C), capturing CO 2 from the atmosphere and storing large quantities of organic matter. Forest carbon reservoirs include aboveground biomass (AGB) and below-ground biomass (BGB), both living and dead [1][2][3]. Susceptible to changes through natural processes and human impacts, biomass monitoring can provide an indication of the sequestration, storage, and/or emission of carbon in the atmosphere [4,5]. Tropical forests' C stocks were estimated of the mangrove canopy height [15]. These spaceborne sensors cover a considerable geographic area, providing publicly available data, but at the cost of reduced accuracy of height and biomass estimates.
Given airborne lidar scanning accuracy, precision, and capacity of repeated acquisitions, it seems particularly well-adapted for characterizing mangrove AGB over large areas. However, some particularities of mangrove ecosystems make them specifically challenging. First, mangroves are periodically flooded, and because conventional lidar does not penetrate water, obtaining a correct digital terrain model can be a problem. The obvious strategy is to collect lidar data during low tide, but this can become a strong constraint in cases where large areas must be surveyed. In the present case, a flight duration of less than two hours and a carefully chosen flight time helped circumvent the problem. Another specificity of the mangrove ecosystem relates to the vegetation structure. The prop roots of the mangrove forest may account for 15-17% of the AGB in mature stands [44]. Because the aerial prop roots are submerged in water for several hours or on a permanent basis, they have been in some cases included in the below-ground biomass category [44] or more commonly included in the aboveground category [45]. This lack of standard procedure affects the consistency between studies. Discounting aerial roots from AGB also creates a discrepancy between the vegetation probed by lidar at low tide and the vegetation contributing to the AGB measured in the field. The estimation of mangrove biomass using lidar data has so far been restricted to a handful of studies [36,37,42,46,47], most of which used satellite lidar only with the notable exception of the study by Fatoyinbo et al. [48]. In the latter study, a coarse lidar metric (highest point per 10 × 10 m grid cell) was used as the sole predictor in the AGB model. The choice of a simple height metric was made in order to allow for metrics from current spaceborne sensors and digital elevation modeling techniques, to be used rather than relatively costlier aerial lidar.
The demonstration of the capacity of aerial lidar scanning to deliver accurate, high-resolution biomass maps over tens of square kilometers in a mangrove forest still remains to be made. The present work investigates the potential of lidar point cloud data to map the aboveground biomass of a mangrove forest with different degrees of disturbance. Our specific research objectives are: to compare different regression methods (partial least squares regression with backward selection of predictors (Auto-PLS) and random forest) for modeling AGB with lidar; to examine how sampling characteristics affect the process of upscaling from plot estimates to landscape estimates of AGB; and to assess the loss of accuracy entailed when using models relying on non-species-specific allometric equations on local and landscape scales.

Materials and Methods
The main steps of the global workflow to produce the AGB map were the following: (i) individual tree AGB was estimated from field measurements, and the values were added for each plot; (ii) lidar metrics were extracted from the cloud of height-normalized points; (iii) the predictive models of plot AGB were adjusted and compared for their performance; (iv) the biomass map was generated based on the best predictive model; and (v) uncertainty of pixel-level and landscape-level predictions was analyzed.

Study Area
The study area is located in the region of Guanabara Bay, Rio de Janeiro State, Brazil ( Figure 1) and represents the last conserved remnant of mangrove forests in a landscape subjected to strong anthropogenic pressure. The climate is hot and humid tropical Atlantic or Aw, according to the Köppen classification, with the rainy season in the austral summer and a drier winter. The region has a mean, minimum, and maximum annual rainfall of 1709 mm, 1155 mm, and 2396 mm, respectively, and a mean annual temperature of 23 • C [49]. The driest month is June with rainfall of 109 mm, and the most humid month is December with a rainfall of 249 mm [50]. Tidal amplitude is less than 2 m (microtidal).
The mangrove forests of the Guanabara Bay (GB) region show a high structural diversity as a result of direct and indirect human action, presenting different degrees of disturbance and Remote Sens. 2018, 10, 637 4 of 21 regeneration stages [51]. The mangrove system studied in this work is located in two protected areas: the Guapimirim Environmental Protection Area (APA Guapimirim) and the Guanabara Ecological Station (ESEC Guanabara). The mangrove forests that belong to these two protected areas show a higher structural development (mean height, live basal area, and mean diameter at breast height) than those in the adjacent regions of the protected areas [52]. Three typical mangrove tree species are found in the region-Avicennia schaueriana Stapf. & Leechman: Acanthacea, Laguncularia racemosa (L.) C.F. Gaertn.: Combretaceae, and Rhizophora mangle L.: Rhizophoraceae-as well as associated species, such as Acrostichum aureum L.: Pteridaceae (a fern) and Hibiscus pernambucensis Arruda: Malvaceae (a shrub).

Ground Data
Vegetation structure information was obtained from 34 plots scattered across the mangrove area inside the APA Guapimirim and ESEC Guanabara, labeled C01, C02, ..., C34 (Figure 1). The plots are distributed along four main rivers (Guapimirim, Guaraí, Caceribu, and Guaxindiba). Their locations were chosen so as to cover a wide range of vegetation biomass but to exclude areas dominated by

Ground Data
Vegetation structure information was obtained from 34 plots scattered across the mangrove area inside the APA Guapimirim and ESEC Guanabara, labeled C01, C02, ..., C34 (Figure 1). The plots are Remote Sens. 2018, 10, 637 5 of 21 distributed along four main rivers (Guapimirim, Guaraí, Caceribu, and Guaxindiba). Their locations were chosen so as to cover a wide range of vegetation biomass but to exclude areas dominated by associated (non-mangrove) species. The size of each plot was adjusted to the structural characteristics of each mangrove stand. Plot sizes varied between 121 and 560 m 2 , depending on the density of individuals and their homogeneity in terms of structural features and species composition. The following tree measures were obtained during field surveys in 2010-2011, according to the methods proposed by [45,53] and described in [54,55]: (i) individual tree height; (ii) stem diameter at breast height (DBH); (iii) identification of the plant at the species level; (iv) living status of the stems (alive or dead). DBH of all stems >1 m tall were obtained using a measuring tape. The measurement of total tree height from the base of the tree to the top of the canopy was obtained with an optical rangefinder or a measuring pole. Mangrove biomass in the study area was calculated using two approaches: (a) using species and local-specific allometric equations (referred to as "species-specific" throughout the text) developed for Rio de Janeiro mangrove forests [45,54], and (b) using pantropical mangrove equations provided by Komiyama et al. [44] and Chave et al. [56] ( Table 1). Note that for all the equations in Table 1, we used the species-specific wood density value from [45,54] (i.e, A. schaueriana = 0.6751 g·cm −3 , L. racemosa = 0.5292 g·cm −3 , and R. mangle = 0.6868 g·cm −3 ). We added the AGB of all trees to derive the total AGB at the plot level. Live trees and standing dead trees were included in the AGB estimation.

Geolocation
The geographical coordinates of the corners of each plot were registered using a GPS Receptor Sokkia model Stratus L1, in a static relative mode [57]. In this method of survey, two GPS Receptors are used, and the coordinates of the points of interest are determined in relation to a referential point with precisely known coordinates, achieving a mean precision of approximately RMS = 0.08 m of the final positioning.

Lidar Data
Lidar data were acquired with an airborne Riegl LMSQ 560 laser scanner and pre-processed by Hansa Geophysics and Aerial Survey. The lidar point cloud covered the entire APA Guapimirim and ESEC Guanabara. The main specifications of the lidar data used in this work are shown in Table 2. The lidar point cloud (in LAS format) was filtered and classified by the data provider into three main classes: (i) ground, (ii) vegetation, and (iii) noise.

Lidar Processing
Lidar data were processed using the software package LAStools 2.1 (http://www.lastools.org/). The main processing steps are summarized as follows: (i) extraction of normalized point cloud height (also referred as height above ground) and (ii) derivation of lidar metrics of each corresponding plot. In the present work, the normalized point cloud heights were derived for the whole study area from the classified pre-processed data set, using the ground data as the bottom reference.
According to previous studies [22,26,[58][59][60], one source of error that may influence the relationship between the lidar-and the field-estimated AGB at different spatial resolutions is the disagreement between the lidar and field plot measurements over which trees or parts of trees are inside the calibration plots. In the lidar measurements, tree crowns are bisected exactly at the plot edge, while in the field plot measurements, an individual is included in the plot if its basal area contributes to the basal area of the plot; in other words, if the DBH trunk point of measure is inside the plot. In our work, for each plot area, two polygons were considered. One was the original polygon area of each plot, and the other included a 5-m buffer extension around the edges of the original plot perimeter. Adding a buffer around the ground plots increases the total area of the lidar data, hopefully capturing more of the plot-based data by systematically including tree crowns located at the edge of some plots. While the area of the 34 original plots ranged from 121 to 560 m 2 (mean = 270 m 2 ), the area of the extended plots ranged from 361 to 1053 m 2 (mean = 614 m 2 ) after adding the 5-m buffer.
The widely used area-based lidar metrics for biomass estimation are height metrics [61]. Height metrics are usually calculated from vegetation returns, which are typically defined as returns with a certain height, such as 0.5 m [62] above the ground surface, thus excluding the contribution of the understory.

Biomass Models: Predicting Biomass on a Plot Level from Lidar Statistics
Different models were tested to predict biomass on a plot level from the lidar statistics (see Table 4): (i) random forest regression (RF) [63] (randomForest Package version 4.6-12 [64] of R software version 3.4.0 [65]), and (ii) partial least squares (PLS) regression with backward selection of predictors (Auto-PLS) [66] (autopls Package version 1.3 [66] of R software version 3.4.0 [65]). RF is a classifier and a regression method based on the generation of a large number of tree-structured classifiers. The RF method is a robust and nonlinear multiple regression formed by growing trees and is a non-parametric statistical method that optimizes predictive accuracy by fitting an ensemble of trees to stabilize model estimates [63]. The analyses of RF results are based on mean square errors (MSE), calculated as the sum of squared residuals divided by number of observation. The importance of the predictor variable for RF can be estimated by a "variable of importance" measurement, calculated by the mean decrease in accuracy, which is based on MSE, and the mean decrease in node impurity, calculated by the residual sum of squares [63]. The RF method uses out-of-bag estimates to monitor error, strength, and correlation. The out-of-bag estimates are based on combining one-third of a training set that removes the need for a separate test set. The out-of-bag estimates can be computed in the same run that constructs the bagged predictor. The result can be compared to cross-validating bagged predictors, and these estimates are close to optimal. The PLS regression is particularly suited for cases when the matrix of predictors has more variables than observations and when there is multicollinearity among predictor values. PLS is a method for relating two data matrices, X and Y, in a multivariate model, but decomposing X and Y into orthogonal scores, finding a new x variable and y variable, which are estimates of the latent variables or their rotations [67,68]. In the PLS method, cross-validation is used to estimate the optimal rank, analyzing the minimum predicted residual sum of squares. Auto-PLS is a wrapper for PLS [63], which incorporates backward variable selection into standard PLS regression. Leave-one-out cross-validation (LOO) was performed for the Auto-PLS methods.
In each case, the predicted model was tested with the lidar metrics derived from the original polygon of each plot (without buffer, named Model 1sp) for species-specific AGB estimate, and the extended polygons (with the 5-m buffer) for all types of AGB estimates: species-specific AGB estimate [45,54] (Model 2sp), Komiyama et al. [44] (Model 3K), and Chave et al. [56] (Model 4C). The accuracy of the different models was evaluated in terms of root-mean-square error (RMSE) and coefficient of determination between prediction and observation (R 2 ). Both quantities are reported, for the calibration (CAL) set (Auto-PLS), as well as for the validation set (all models). Validation statistics (R 2 _val and RMSE_val) use the LOO cross-validation procedure for Auto-PLS and the bootstrapped out-of-bag error for RF. In addition, the observed (Y-axis) versus predicted (X-axis) AGB regression line was plotted, and its departure from the 1:1 line was taken as an indicator of model bias [69]. Table 4. Designation of the regression models tested for predicting biomass on a plot level from the lidar statistics. Model 1 species-specific (M1sp), Model 2 species-specific (M2sp), Model 3 pantropical_K (M3K), and Model 4 pantropical_C (M4C) using 34 plots for (a) Auto-PLS, (b) and random forest (RF).

Model
Auto

Mapping Biomass on a Landscape Level from Lidar Statistics
The biomass mangrove map of APA Guapimirim and ESEC Guanabara was obtained through the application of the best predictive model using only the lidar metrics calculated for the landscape level in a 25 m × 25 m cell grid. The polygon and delineation of the mangrove area over which the model was applied was taken from Arasato et al. [70]. In their study [70], spectral information (from a high-resolution WorldView2 image) and lidar data were combined to map land use, landcover, and mangrove types in APA Guapimirim and ESEC Guanabara on a landscape level. The authors used an object-based image analysis (OBIA) approach using spectral, geometric, and textural attributes derived from lidar and optical imagery data. In addition, they visited 150 geo-located points of mangrove to help with the visual interpretation and the validation of the maps. Based on their mangrove-type map, we excluded non-mangrove areas (water bodies, open areas, and large areas of non-mangrove vegetation) from the final biomass mangrove map.

Sample Plot Coverage Assessment
We evaluated how well the 34 sample plots captured the mangrove structural variability in the study area using a principal component analysis (PCA) approach, similar to [71,72]. We performed PCA on the lidar statistical metrics dataset for the whole mangrove study area as a dimension reduction method by finding the principal components of the input data. The lidar metric values were used to predict the field plot coordinates that were overlaid on the space plane of principal component 1 (PC1) and principle component 2 (PC2). The result is the field plot distribution on a space plane of PC1 versus PC2, obtained with a set of lidar statistical metrics for mangrove forest area. This type of analysis can be helpful in identifying pixels in the predictor space that are poorly represented in the training set and consequently, may be poorly predicted by the model.

Evaluating Uncertainty on the Landscape Level
A "design-based model assisted" inference framework was used to estimate the uncertainty of the biomass on the landscape level [29]. The lidar data serve as auxiliary variables via the best regression model trained on the sample plots. The standard general regression estimator (GREG) of biomass and its variance are then used to build a population (landscape) level estimate of biomass and its uncertainty [73].
The general regression estimators of population (landscape) biomass mean and its variance are given by Equations (1) and (2) [29].B Remote Sens. 2018, 10, 637 where N is the population size (N = 84,752, the number of 25 m × 25 m pixels grid), and n is the number of sampling units in sampling set S (either 34 or 17). ∑ k∈S is the estimated error of prediction at plot scale (MSE of LOO or two-fold cross-validation). GREG estimates of landscape-level mean and variance were computed for the AGB predictions from regression models M2sp, M3K, and M4C.
We also analyzed the prediction discrepancy between models per area dominated by each mangrove species in the APA Guapimirim and ESEC Guanabara. The mangrove species mask for that analysis was obtained from [70]. We report the difference in the prediction of AGB maps obtained with pantropical (M3K and M4C) and the AGB map obtained with M2sp.
In addition to the standard procedure outlined above, an alternative resampling approach was used to estimate the uncertainty of the predictions on a landscape scale. The 34 sampling plots were randomly divided in two subsets with each one covering most of the range of field-estimated AGB by stratified random sampling. Plots were first ordered by increasing values of AGB. They were then split into 17 strata of two consecutive observations. One out of two observations per stratum was then randomly allocated either to sample 1 or sample 2, thus creating two subsets of sample plots of equal size and balanced with regard to AGB. For each random draw, each one of the two subsets of field plots was used to predict landscape-scale biomass (i.e., fitting the prediction model and applying the model on a landscape level). The difference in average AGB per hectare at landscape level, predicted by the paired models, was recorded. The entire process (random splitting, model adjustment, and recording of the difference in landscape-level prediction) was repeated 100 times. The distribution of the mean square difference between the paired predictions reflects the uncertainty of prediction on a landscape level.

Results
The mean AGB value of all 34 field plots using species-specific equations was 125.71 t·ha −1 , ranging from a minimum of 43.56 t·ha −1 to a maximum of 186.10 t·ha −1 (Figure 2a The AGB estimates obtained with pantropical allometries (Chave et al. [56] and Komiyama et al. [44]) essentially differ from one another by a scaling factor (Table 1). However, the plot AGB predictions of both models differ markedly from the predictions obtained using species-specific allometries (see Figure 3). Figure 4 shows the correlation between structural parameters, AGB species-specific metrics, and lidar metrics for 34 plots. Figure 5   The AGB estimates obtained with pantropical allometries (Chave et al. [56] and Komiyama et al. [44]) essentially differ from one another by a scaling factor (Table 1). However, the plot AGB predictions of both models differ markedly from the predictions obtained using species-specific allometries (see Figure 3). Figure 4 shows the correlation between structural parameters, AGB species-specific metrics, and lidar metrics for 34 plots. Figure 5       The results show that for all the regression models tested for predicting biomass, the best models (considered here with the highest R 2 and the lowest RMSE) were obtained using the AGB of speciesspecific equations, such as Model 1sp and Model 2sp. Model 3K and Model 4C were generated with AGB calculated using pantropical equations and presented much lower accuracy (see Table 5) with roughly doubled RMSE and halved R 2 .
The inclusion of the 5-m buffer around the polygons of the field plots improved model performance (see Table 5). The best model (highest R 2 and lowest RMSE) was the M2sp.autopls (Auto-PLS with extended polygons) with a R 2 (LOO) of 0.803, R 2 (CAL) of 0.89, RMSE(LOO) of 14.8 t·ha −1 , and RMSE(CAL) of 11.5 t·ha −1 . The plot of the 1:1 line in Figure 6 indicates that model M2sp.autopls has the lowest bias. The M2sp.rf model (RF with extended polygons) was relatively less successful with a R 2 = 0.71 and an out-of-bag error of 17.8 t·ha −1 . The eight variables included in the best performing M2sp.autopls model were as follows: avg, min, max, d02, d03, d04, d05, and d08. Figure 4 shows the correlation of these variables. The ground AGB values were most strongly positively correlated with avg and max and negatively with d02, while all the other variables showed an absolute correlation coefficient lower than 0.8. The eight most important variables for the M2sp.rf model were p99, p50, max, d06, p95, d08, avg, and qav, in that order. Including species dominance information as an additional variable did not improve the model prediction capability (not shown); although, wood density among mangrove species varies from 0.53 g·cm −3 (Laguncularia racemosa) to 0.68 g·cm −3 (Rhizophora mangle) [45,54]. The mangrove aboveground biomass map (Figure 7) was generated using the best performing M2sp.autopls model. The most impacted mangrove areas with lower AGB values were observed near urban and degraded sites (the right sector in the map). Table 5. Results of the regression models tested for predicting biomass. Model 1 species-specific (M1sp), Model 2 species-specific (M2sp), Model 3 pantropical_K (M3K), and Model 4 pantropical_C (M4C) using 34 plots for (a) Auto-PLS and (b) random forest. (Calibration (CAL), Leave-one-out crossvalidation (LOO), root-mean-square error (RMSE) and coefficient of determination (R 2 )).  The results show that for all the regression models tested for predicting biomass, the best models (considered here with the highest R 2 and the lowest RMSE) were obtained using the AGB of species-specific equations, such as Model 1sp and Model 2sp. Model 3K and Model 4C were generated with AGB calculated using pantropical equations and presented much lower accuracy (see Table 5) with roughly doubled RMSE and halved R 2 .

Model Auto-PLS (a) Random Forest (b)
The inclusion of the 5-m buffer around the polygons of the field plots improved model performance (see Table 5). The best model (highest R 2 and lowest RMSE) was the M2sp.autopls (Auto-PLS with extended polygons) with a R 2 (LOO) of 0.803, R 2 (CAL) of 0.89, RMSE(LOO) of 14.8 t·ha −1 , and RMSE(CAL) of 11.5 t·ha −1 . The plot of the 1:1 line in Figure 6 indicates that model M2sp.autopls has the lowest bias. The M2sp.rf model (RF with extended polygons) was relatively less successful with a R 2 = 0.71 and an out-of-bag error of 17.8 t·ha −1 . The eight variables included in the best performing M2sp.autopls model were as follows: avg, min, max, d02, d03, d04, d05, and d08. Figure 4 shows the correlation of these variables. The ground AGB values were most strongly positively correlated with avg and max and negatively with d02, while all the other variables showed an absolute correlation coefficient lower than 0.8. The eight most important variables for the M2sp.rf model were p99, p50, max, d06, p95, d08, avg, and qav, in that order. Including species dominance information as an additional variable did not improve the model prediction capability (not shown); although, wood density among mangrove species varies from 0.53 g·cm −3 (Laguncularia racemosa) to 0.68 g·cm −3 (Rhizophora mangle) [45,54]. The mangrove aboveground biomass map (Figure 7) was generated using the best performing M2sp.autopls model. The most impacted mangrove areas with lower AGB values were observed near urban and degraded sites (the right sector in the map).

Sample Coverage Assessment
Field plot representativeness was evaluated graphically (Figure 8). The full range of structural variability in the mangrove forest area is represented by grey circles (10,000 grid cells randomly selected from the entire mangrove area and by black circles the 1000 cells randomly selected from the non-mangrove area), while the field plots are represented by red circles and labels. The first two principal components jointly explained 64% of the total variance. Field plots covered almost all of the range of structural variability of the mangrove forest in the study area (Figure 8). The 95% confidence ellipse (solid line) covers 88% of the mangrove structural variability, and the 98% confidence ellipse (dashed line) covers 96% of all the field plots. The area outside the confidence ellipse ( Figure 8) appears to be mostly dominated by the non-mangrove forest, such as associated species (black points). Figure 7 shows the mangrove area and also the non-mangrove area (in black).

Sample Coverage Assessment
Field plot representativeness was evaluated graphically (Figure 8). The full range of structural variability in the mangrove forest area is represented by grey circles (10,000 grid cells randomly selected from the entire mangrove area and by black circles the 1000 cells randomly selected from the non-mangrove area), while the field plots are represented by red circles and labels. The first two principal components jointly explained 64% of the total variance. Field plots covered almost all of the range of structural variability of the mangrove forest in the study area (Figure 8). The 95% confidence ellipse (solid line) covers 88% of the mangrove structural variability, and the 98% confidence ellipse (dashed line) covers 96% of all the field plots. The area outside the confidence ellipse ( Figure 8) appears to be mostly dominated by the non-mangrove forest, such as associated species (black points). Figure 7 shows the mangrove area and also the non-mangrove area (in black). The plot vegetation average return height (avg, 34 plots with 5-m buffer) was further compared to landscape-level average return height (mangrove area covered by lidar with spatial resolution of 25 m, see Figure 9), suggesting that low mangrove was slightly under-sampled.  The plot vegetation average return height (avg, 34 plots with 5-m buffer) was further compared to landscape-level average return height (mangrove area covered by lidar with spatial resolution of 25 m, see Figure 9), suggesting that low mangrove was slightly under-sampled.

Sample Coverage Assessment
Field plot representativeness was evaluated graphically (Figure 8). The full range of structural variability in the mangrove forest area is represented by grey circles (10,000 grid cells randomly selected from the entire mangrove area and by black circles the 1000 cells randomly selected from the non-mangrove area), while the field plots are represented by red circles and labels. The first two principal components jointly explained 64% of the total variance. Field plots covered almost all of the range of structural variability of the mangrove forest in the study area (Figure 8). The 95% confidence ellipse (solid line) covers 88% of the mangrove structural variability, and the 98% confidence ellipse (dashed line) covers 96% of all the field plots. The area outside the confidence ellipse ( Figure 8) appears to be mostly dominated by the non-mangrove forest, such as associated species (black points). Figure 7 shows the mangrove area and also the non-mangrove area (in black). The plot vegetation average return height (avg, 34 plots with 5-m buffer) was further compared to landscape-level average return height (mangrove area covered by lidar with spatial resolution of 25 m, see Figure 9), suggesting that low mangrove was slightly under-sampled.

Uncertainty at the Landscape Level
The general regression estimators of population (landscape) biomass mean and its variance are given by Equations (1) and (2) in Section 2. In Table 6, we report the AGB estimates and uncertainty in terms of mean per hectare over landscape (rather than the cumulated AGB) (i.e.,B/N and standard error V (B) N 2 for M2sp.autopls, M3K.autopls, and M4C.autopls). For the best model, M2sp.autopls, the estimate of the mean AGB from resampling is higher than the GREG estimate (see Table 6). However, confidence intervals of both estimates largely overlap. The prediction and discrepancy between the models per area dominated by each mangrove species in the APA Guapimirim and ESEC Guanabara are shown in Table 7. Table 7 reports the difference in the prediction of AGB maps obtained with pantropical equations (M3K.autopls and M4C.autopls) and the AGB map obtained with M2sp.autopls.

Discussion
The models trained with AGB species-specific equations (M1 and M2) vastly out-performed the others models trained with AGB non-species-specific equations (M3 and M4). This indicates that these species-specific equations better explained the structure and AGB of the mangrove of the study area. Models trained with AGB estimates derived from pantropical allometric equations (M3 and M4) appear to be quite severely biased ( Figure 6).
Model M2sp.autopls presented the best performance with R 2 = 0.803 and RMSE(LOO) = 14.7 t·ha −1 RMSE(CAL) = 11.2 t·ha −1 and was used to map the mangrove AGB of APA Guapimirim and ESEC Guanabara (Figure 7). In general, the lowest biomass values were observed upstream in the rivers near the continental borders of the mangrove forests, on more degraded areas, and with a higher occurrence of invasive and associated species. Higher AGB values are observed in the central area of the mangrove, which is also the best preserved area included in the ESEC Guanabara (Figure 7). Those distribution patterns are in agreement with a previous study [52]. It is interesting to note very low AGB values in the mangrove forests located along the coast of some rivers. This pattern is the result of a pest of Lepdoptera caterpillars that occurred in 2009 and severely impacted the A. schaueriana population, which was dominant in this area of the estuary.
In the present study, the maximum height (max) and average height (avg) are important variables for mangrove biomass prediction with the M2sp.autopls and M2sp.rf models. Point density metrics (d02, d03, d04, d05, and d08) were also influential variables for M2a. Point density metrics were also used in the final models in [60] for predicting the vertical distribution of heights and penetration depths in an area-based approach. The best independent variables to estimate AGB for a boreal forest showed by [23] were one variable related to canopy height and another one related to canopy density (similar to point density here). A simple canopy height metric was also used to predict AGB in East African mangroves [48].
An analysis of reported biomass accuracy estimates of terrestrial vegetation from more than 70 referenced articles using different remote sensing techniques was conducted in [22]. Studies based on discrete return lidar (DRL) and full return lidar (FRL) had a mean R 2 of 0.76 and 0.80, respectively, with mean residual standard error (RSE) of 39.4 t·ha −1 for DRL and 50.2 t·ha −1 for FRL. The average relative residual standard error (RSE%), calculated as the RSE standardized by mean AGB from field measurements, was 27.0% and 22.3% for DRL and FRL, respectively [22]. For our M2sp.autopls model, the RSE% was 8.9%. The RSE refers to the absolute model error, while RSE% expresses model performance relative to the mean biomass from field measurements. The M2sp.autopls model showed a good estimation of AGB, with a lower RMSE (11.17) than the mean RSE (39.4) of all the DRL studies, and a RSE% of 8.9% lower than 20% (or ±20 t·ha −1 , the greater of the two), which is considered the accuracy requirement of a global forest biomass mapping mission for at least 80% of grid cells [22]. Compared to other studies [22], the M2sp.autopls model showed better results than those found for the following: coniferous forest (R 2 = 0. 64 [77]; R 2 = 0.78 [78]; and R 2 = 0.72, RMSE = 40.2 t·ha −1 [33]). Other studies reported higher accuracies, such as in rainforest (R 2 = 0.90, RMSE = 38.3 t·ha −1 ) [79] and temperate forest (R 2 = 0.89, RMSE = 50.2 t·ha −1 [80] and R 2 = 0.93, RMSE = 33.9 t·ha −1 [81]), among other results reported in the literature. In comparable settings, the study of [48] reports relative RMSE errors of 23-33% to be compared with the calibration RMSE(LOO) of 12%.
Spatial mismatch between inventory data and associated canopy area may be particularly large with small plots due both to more severe border effects and the higher impact of positioning errors. The models used here that considered the 5-m buffer extension (M2sp.autopls and M2sp.rf) showed better results than the models that used the original plot perimeters (M1sp.autopls and M1sp.rf). The increase in plot size enables the inclusion of tree crowns with trunks within the plot but with part of the crown outside, while the exact boundary of the plot on the lidar data cannot include all the tree crowns of the plot. Conversely, extending the point cloud beyond the field plot limits led to the inclusion of trees and crowns located outside the plot. This suggests that model improvement cannot be ascribed to a better match between the trees sampled within the field plot and the area covered by the lidar data. Rather, the benefit of extending the lidar area (and the point cloud size) probably arose from improving the lidar metric estimates. On average, the point cloud size was doubled, and it was even tripled for the smallest plots. As long as the surrounding areas remained similar to the core plot area (in terms of structure), including a buffer around the plot tended to decrease the uncertainty of the lidar metrics and thus improve the model by the sheer effect of increasing the lidar point sample size. Plot sampling should be designed to capture the entire range of possible structural variation encountered by the lidar data. Our sample coverage assessment can be considered as satisfactory (see Figures 8 and 9).
The estimated landscape-level error in the AGB prediction using the GREG estimator is 2.5% (Table 6) (i.e., about one-fifth of the plot-level error of 12% (see Table 5)). However, the GREG estimated uncertainty might be slightly underestimated. Indeed, Figure 6a suggests that the plot-level model may be slightly biased. Any plot-level model bias will translate into an error on the landscape level, which is not accounted for by the GREG estimator. The resampling approach provides an estimate of the mean biomass on a landscape scale slightly higher than the one estimated via GREG (106.7 vs. 105.0 t·ha −1 ) with a higher standard deviation (6.7% vs. 2.5%). A single plot (C025) may be responsible for the bulk of the difference between the two estimates. Plot C025, with the lowest measured biomass (see Figure 5), is responsible for the slight apparent bias in the plot-level model (see Figure 6a), as AGB for this plot is overestimated by all the models, albeit only slightly by model M2sp.autopls. The resampling technique, which compares predictions made from subsets of training plots either including or not this outlier observation, will be sensitive to the high leverage of this particular observation. As a matter of fact, the mean prediction for the 100 random subsets, which include the C025 plot, is 103.1 t·ha −1 , as opposed to a mean of 110.1 t·ha −1 for the 100 random subsets without C025.
The large bias detected in Model 3K and Model 4C (see Figure 6) indicates that the assumptions required for calculating the GREG estimator are not met, and indeed, the landscape uncertainty estimations are irrelevant (and inconsistent with the observed prediction discrepancies between the various models). The landscape errors in M3K and M4C appear to be on the order of +16 t·ha −1 or +14% (M3K) and −22 t·ha −1 or −18% (M4C). The error can be even larger in areas dominated by a particular species (see Table 7). Hence, the inaccuracy in the field data (resulting from the use of poorly adapted allometries) yields bias in the prediction models, which translates into large errors when such models are applied on a landscape scale. The large prediction errors reported in [48], and the bias detected with height class (underestimation of low AGB and overestimation of high AGB values) may also be a consequence of not using species-specific allometric equations.
Despite a reasonably good coverage (see Figure 8), the unbalanced sampling across vegetation height classes and the notable underrepresentation of low mangrove (see Figure 6a) have probably impaired the landscape-level estimate accuracy. A better coverage of low mangrove (predominantly juvenile stages) could have been ascertained using a canopy height map when the field work was conducted, thereby potentially reducing the uncertainty of AGB on a landscape scale.
While under sampling of low canopy mangrove seems to have impaired the model performance, the small size of the plots seems not to have done so. This was not entirely expected as previous studies have shown an inverse relationship between model errors and plot size [22,27]. The use of large plots reduces errors in terms of RMSE [26]. Extending the area from which the lidar statistics were derived slightly outside the plot borders improved prediction. This suggests that smaller plots should probably be avoided. Furthermore, such small plots were adapted to the low canopy (<18 m) characteristic of this mangrove area, which represents forests with low structural development and may not be recommended for higher canopy forests (structurally well-developed forests).

Conclusions
The present study underscores the effectiveness of lidar as a means to estimate mangrove forest AGB in areas with different degrees of disturbance. The model (M2sp.autopls) used to map the mangrove AGB of APA Guapimirim and ESEC Guanabara presented the best performance with R 2 = 0.80% and 8.9% (RMSE CAL %) for plot-level error. The final uncertainty (RMSE LOO %) of plot-level prediction was~12%, whereas landscape uncertainty was reduced to 2.5-6% (depending on the method used). This best model (M2sp.autopls) was trained with AGB obtained with species and local-specific equations and extended plot perimeter (5-m buffer extension).
The most important factor found to affect the quality of the predictions is the quality of the field data. Inaccuracy in the field data yields bias in the prediction models, which translates into large errors when such models are applied on a landscape scale. In this study, the models trained with AGB species-specific equations (M1 and M2) presented more accurate results than models trained with pantropical specific equations. Species-specific allometric scaling equations seem to be required to build accurate (unbiased) models. The combination of the strong degree of spatial segregation of tree species common in mangrove forests and the small size of the sampling plots may have exacerbated this sensitivity to inaccurate (non-species-specific) allometric scaling equations. Nonetheless, because the uneven distribution of tree species in a terrestrial forest ecosystem is also the rule rather than the exception, similar problems-in nature if not in amplitude-are most likely to affect AGB estimates obtained for non-mangrove forests for which species-specific allometry are not available.
Lidar can provide accurate maps of mangrove AGB and can be a useful tool for mangrove monitoring, carbon stock assessments, and coastal management in general. Lidar is highly accurate but relatively expensive. It could be used to efficiently "train" more cost-effective solutions, such as textural analysis of high-resolution optical imagery [26]. Such imagery is available at a much lower cost. Textural features may be further combined with multi-spectral information to produce robust landcover maps. In this sense, future work could compare lidar-derived biomass estimates with other remote sensing data and techniques, such as SAR and optical imagery. Lidar data can also be merged with other remote sensing data to predict biomass and to quantify mangrove changes in terms of AGB and carbon stock.