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

Abstract

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² mangrove forest in Rio de Janeiro, Brazil. Plot AGB was computed using either species-specific or non-species-specific 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-specific allometric models were distinctly more accurate (R²(calibration) = 0.89, R²(validation) = 0.80, root-mean-square error (RMSE, calibration) = 11.20 t·ha⁻¹, and RMSE(validation) = 14.80 t·ha⁻¹). The use of non-species-specific 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 field 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-specific allometric equations are available.

1. Introduction

Forests play a crucial role in the global carbon cycle (C), capturing CO₂ 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 to be ~471 ± 93 PgC—that is 55% of the total C of the world’s forests [6]. This estimation included C in live biomass (above and below ground), soil carbon, deadwood, and litter.

Mangroves are a typical coastal ecosystem occurring at the interface of land and sea in tropical and subtropical regions; they are also highly productive and present complex physiological and growth processes [7,8]. Mangrove communities and processes can be influenced by salinity, insolation, nutrients, fresh water, flooding, soil anoxia, and tidal action. Mangrove forests store more carbon than other ecosystems on an area-specific basis, with a mean carbon stock estimate of 956 t·C·.ha⁻¹ compared with 241 t·C·ha⁻¹ for rainforests [1,9]. Soil stocks are typically much larger than in most other forests and are estimated to account for 49–98% of the total carbon storage [1]. Deforestation and aquaculture are major threats to mangrove forests [10]. Mangrove areas are declining and are today one of the most threatened environments on Earth [11]. All human threats on mangroves can impact the mangrove carbon cycle.

Aboveground biomass, defined as all living biomass above the ground, including stem, stump, branches, bark, seeds, and foliage [2], is used to estimate terrestrial carbon pools and carbon sequestration from the atmosphere [4,5]. Plot AGB estimates are derived from individual tree measurements and allometric equations derivation and extrapolated on a larger scale based on remote sensing techniques, such as optical imagery [12], synthetic aperture radar (SAR) [13], or light detection and ranging (lidar) surveys [14,15,16]. Lidar is an active system that has been used to estimate and map forest structure, aboveground biomass, and forest carbon [15,17,18,19,20,21,22,23,24]. When developing a spatially explicit model of biomass to be applied on a landscape scale, there are many sources of errors to consider. The first and possibly foremost source of uncertainty stems from the limited accuracy of the field-based estimates of AGB [25]. This is notably due to the lack of species-specific allometric equations for the vast majority of species found in natural tropical forests. However, mangrove forests, which have many fewer species compared to other tropical forests, provide an interesting opportunity to explore the uncertainty accruing with the use of biome-specific rather than species-specific allometric scaling equations. Many other sources of uncertainty can be listed. Spatial and temporal inconsistency between lidar and field observations may stem from inaccuracies in global positioning system (GPS) records, temporal discrepancy between different data acquisitions, disagreement between lidar and field plots edge measurements, and calibration plots that include or exclude parts of the trees [19,26]. Furthermore, errors in biomass models can be influenced by plot sizes [26,27], forest types, and sensor types [22]. However, the use of lidar systems for estimating forest aboveground biomass tends to show better results than the use of optical or SAR data [22]; although, optical and SAR data provide a greater coverage area as an advantage. Due to the numerous sources of errors that may affect the modeling process, the evaluation of a model’s uncertainty at various spatial scales is an inherent part of model development [19,26,27,28]. In particular, the estimation of the uncertainty of predictions made over large areas from sample plots can be derived from classical sampling theory and the use of the general regression estimator, with lidar metrics serving as auxiliary variables [29].

Lidar has already proven its ability to map AGB and to extract biophysical parameters in many forest types, such as deciduous [24,30], coniferous [31,32], and tropical forests [17,20,33]. In the specific case of mangrove forests, lidar has previously been used for mapping the extension of the mangrove, detecting gaps [34], and estimating and mapping structural parameters and AGB from the integration of lidar and other remote sensing data [15,35,36,37,38]. Lidar has been also used to detect individual mangrove trees and estimate their structure [39,40,41]. One of the advantages of using airborne lidar data is the possibility of extracting surface and terrain models and, consequently, obtaining the canopy height model in a given forest area. Hence, it is possible to obtain a three-dimensional estimate of the forest structure. The use of lidar and interferometric synthetic aperture radar (InSAR) derived Shuttle Radar Topography Mission (SRTM) data for biomass estimation in large mangrove areas has also been reported [15,37,42,43]. Ice, Cloud, and Land Elevation Satellite/Geoscience Laser Altimeter System (ICEsat/GLAS) data were used to calibrate the SRTM elevation data and to produce a map 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.

2. 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.

2.1. 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 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).

2.2. 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 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², 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. Mangrove aboveground biomass (AGB) equations, using diameter at breast height (DBH) and wood density (WD).

Equations
Species-Specific-AGB for Live Tree
A. schaueriana	Ln (AGB total) = 4.8017 + 2.5282 × Ln (DBH)	Estrada et al. (2014)	[54]
L. racemosa	Ln (AGB total) = 5.2394 + 2.2792 × Ln (DBH)	Soares et al. (2005)	[45]
R. mangle	Ln (AGB total) = 5.2985 + 2.4810 × Ln (DBH)	Soares et al. (2005)	[45]
Species-Specific-AGB for Dead Tree
A. schaueriana	Ln (AGB total) = 4.4117 + 2.5578 × Ln (DBH)	Estrada et al. (2014)	[54]
L. racemosa	Ln (AGB total) = 4.9308 + 2.2951 × Ln (DBH)	Soares et al. (2005)	[45]
R. mangle	Ln (AGB total) = 4.9851 + 2.5142 × Ln (DBH)	Soares et al. (2005)	[45]
Pantropical-AGB Total
AGB pantropical	AGB total = 251 × WD × (DBH 2.46)	Komiyama et al. (2008)	[44]
AGB pantropical	AGB total = 167.6 × WD × (DBH 2.47)	Chave et al. (2005)	[56]

). 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⁻³, L. racemosa = 0.5292 g·cm⁻³, and R. mangle = 0.6868 g·cm⁻³). 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.

2.3. 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.

2.4. 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. Light detection and ranging (Lidar) data main specifications.

Products	Main Specifications
Laser	Riegl LMS Q560
Lidar point cloud	LAS format
Total surveyed area	90 km2
Geodesic reference system	WGS84
Projection system	Universal Transverse Mercator (UTM)
Point density	5 pulse m−2 (6 points m−2)
Altimetry precision	15 cm
Planimetric accuracy	50 cm
Swath angle	60° (±30°)
Day and time of data acquisition Coordinated universal time (UTC)	12 November 2012-17:20–19:02 UTC 12 December 2012-15:29–18:10 UTC
Tide height (tide gauge in GB)	12 November 2012-17:20 UTC~0.86 m, 19:02 UTC~0.34 m 12 December 2012-15:29 UTC~0.89 m, 13:53 h~1.1 m, 18:10 UTC~0.9 m

. 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.

2.5. 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² (mean = 270 m²), the area of the extended plots ranged from 361 to 1053 m² (mean = 614 m²) 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.

Based on canopy height, the following lidar metrics (

Table 3. Lidar metrics calculated from point clouds and used on predictive models.

Lidar Metrics		Lidar Metrics
Avg	Mean height	std	Standard deviation
d00	Density points of 0.5–2 m	min	Minimum height
d01	Density points of 2–4 m	p01	Height percentile of 1%
d02	Density points of 4–6 m	p05	Height percentile of 5%
d03	Density points of 6–8 m	p10	Height percentile of 10%
d04	Density points of 8–10 m	p25	Height percentile of 25%
d05	Density points of 10–12 m	p50	Height percentile of 50%
d06	Density points of 12–14 m	p75	Height percentile of 75%
d07	Density points of 14–16 m	p90	Height percentile of 90%
d08	Density points of 16–18 m	p95	Height percentile of 95%
dns_gap	Density gap	p99	Height percentile of 99%
kur	Kurtosis	qav	Quadratic mean
max	Maximum height	ske	Skewness

) were calculated for the point clouds (with or without the buffer) corresponding to each plot: height percentiles (01, 05, 10, 25, 50, 75, 90, 95, 99 percentile, denoted as P01, P05, …, P99, respectively); average (avg), minimum (min) and maximum (max) height; several statistics describing the height distribution, such as coefficient of skewness (ske), standard deviation (std), and coefficient of kurtosis (kur); quadratic mean height (qav); density gap (dns-gap), defined as the number of points below the cover cutoff (50 cm) divided by the number of all returns; and density metrics (relative height density calculated by dividing counts per the total number of points and scaling to a percentage) of d00 (0.5–2 m), d01 (2–4 m), d02 (4–6 m), d03 (6–8 m), d04 (8–10 m), d05 (10–12 m), d06 (12–14 m), d07 (14–16 m), and d08 (16–18 m). The same metrics were calculated on a landscape level for the entire mangrove study area on a 25 m × 25 m cell grid.

2.6. 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. 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-PLS (a)	RF (b)
M1sp species-specific	34 plots (without buffer)	M1sp.autopls	M1sp.rf
M2sp species-specific	34 plots with 5 m buffer (extended polygons)	M2sp.autopls	M2sp.rf
M3K pantropical_K		M3K.autopls	M3K.rf
M4C pantropical_C		M4C.autopls	M4C.rf

): (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²). Both quantities are reported, for the calibration (CAL) set (Auto-PLS), as well as for the validation set (all models). Validation statistics (R²_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].

2.7. 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.

2.8. 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.

2.9. 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].

$$\widehat{B}={\displaystyle \sum }_{k=1}^N{\widehat{y}}_k+\frac{N}{n}{\displaystyle \sum }_{kϵS}{y}_k-{\widehat{y}}_k$$

(1)

$$\widehat{V}\left(\widehat{B}\right)=N^2\left(\frac{1}{n}-\frac{1}{N}\right){\displaystyle \sum }_{k\in S}\frac{{\left(y_k-{\widehat{y}}_k\right)}^2}{n-1}$$

(2)

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).${{\displaystyle \sum }}_{k\in S}\frac{{\left(y_k-{\widehat{y}}_k\right)}^2}{n-1}$ 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.

3. Results

The mean AGB value of all 34 field plots using species-specific equations was 125.71 t·ha⁻¹, ranging from a minimum of 43.56 t·ha⁻¹ to a maximum of 186.10 t·ha⁻¹ (Figure 2a). The median value was 121.92 t·ha⁻¹, and the first and third quantiles were 99.71 t·ha⁻¹ and 150.29 t·ha⁻¹, respectively.

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 shows the lidar point cloud for three mangrove field plots with contrasting AGB values: plot 25 that presents the lower AGB_ssp value with 45.65 t·ha⁻¹, plot 09 with AGB_ssp value of 126.37 t·ha⁻¹, and plot 20 the highest AGB_ssp with 186.10 t·ha⁻¹.

The results show that for all the regression models tested for predicting biomass, the best models (considered here with the highest R² 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. 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 cross-validation (LOO), root-mean-square error (RMSE) and coefficient of determination (R²)).

Model	Auto-PLS (a)	Random Forest (b)
M1sp 34 plots without buffer species-specific	RMSE(CAL) = 14.70 RMSE(LOO) = 17.30 RMSE% =11.69% R2(CAL) = 0.80 R2(LOO) = 0.73	RMSE = 18.60 RMSE% =14.79% R2 = 0.68
M2sp 34 plots with 5-m buffer (extended polygons) species-specific	RMSE(CAL) = 11.17 RMSE(LOO) = 14.80 RMSE% = 8.88% R2(CAL) = 0.89 R2(LOO) = 0.80	RMSE = 17.86 RMSE% = 14.20% R2 = 0.71
M3K pantropical_K	RMSE(CAL) = 22.50 RMSE(LOO) = 24.90 RMSE% = 18.21% R2(CAL) = 0.39 R2(LOO) = 0.25	RMSE = 25.99 RMSE% =21.04% R2 = 0.18
M4C pantropical_C	RMSE(CAL) = 15.50 RMSE(LOO) = 17.10 RMSE% = 18.31% R2(CAL) = 0.39 R2(LOO) = 0.25	RMSE = 17.60 RMSE% =20.79% R2 = 0.20

) with roughly doubled RMSE and halved R².

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² and lowest RMSE) was the M2sp.autopls (Auto-PLS with extended polygons) with a R²(LOO) of 0.803, R²(CAL) of 0.89, RMSE(LOO) of 14.8 t·ha⁻¹, and RMSE(CAL) of 11.5 t·ha⁻¹. 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² = 0.71 and an out-of-bag error of 17.8 t·ha⁻¹. 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⁻³ (Laguncularia racemosa) to 0.68 g·cm⁻³ (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).

3.1. Landscape-Level Accuracy

3.1.1. 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.

3.1.2. 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. Mean AGB per hectare (t·ha⁻¹) and standard error (SE) on a landscape scale using either a generalized regression estimator or a resampling technique. The resampling technique repeatedly fits the model to a stratified subsample of 17 plots only. (Standard general regression estimator (GREG)).

Model	Landscape Mean (GREG)	Landscape SE (GREG)	Landscape Mean (Resampling)	Landscape SE (Resampling)
M2sp.autopls	105.04	2.54	106.72	6.67
M3K.autopls	121.60	4.27	122.20	11.68
M4C.autopls	83.29	2.93	83.78	6.98

, we report the AGB estimates and uncertainty in terms of mean per hectare over landscape (rather than the cumulated AGB) (i.e., $\hat{\mathrm{B}}/\mathrm{N}$ and standard error $\sqrt{\frac{\hat{\mathrm{V}}\left(\hat{\mathrm{B}}\right)}{{\mathrm{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. AGB map of error for analysis of each class of mangrove species.

Species	Map Area	Mean AGB M2sp.autopls	Mean AGB M3K.autopls	Mean AGB M4C.autopls	RMSD K	RMSD C	Mean Error K	Mean Error C
Av	156,800	93.27	112.97	77.35	30.08	28.57	19.70	-15.92
Lg	606,775	93.041	118.78	81.34	33.43	25.74	26.26	-11.70
Rh	1026,500	125.13	127.94	87.67	23.75	44.408	2.81	-37.47

. 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.

4. 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² = 0.803 and RMSE(LOO) = 14.7 t·ha⁻¹ RMSE(CAL) = 11.2 t·ha⁻¹ 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² of 0.76 and 0.80, respectively, with mean residual standard error (RSE) of 39.4 t·ha⁻¹ for DRL and 50.2 t·ha⁻¹ 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⁻¹, 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² = 0.64, RMSE = 28.5 t·ha⁻¹ [74] and R² = 0.75, RMSE = 45.6 t·ha⁻¹ [75]); deciduous forest (R² = 0.71, RMSE = 39.3 t·ha⁻¹) [75]; taiga forest (R² = 0.72, RMSE = 14.2 t·ha⁻¹) [76]; and tropical forest (R² = 0.36, RMSE = 22.8 t·ha⁻¹ [77]; R² = 0.78 [78]; and R² = 0.72, RMSE = 40.2 t·ha⁻¹ [33]). Other studies reported higher accuracies, such as in rainforest (R² = 0.90, RMSE = 38.3 t·ha⁻¹) [79] and temperate forest (R² = 0.89, RMSE = 50.2 t·ha⁻¹ [80] and R² = 0.93, RMSE = 33.9 t·ha⁻¹ [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 Figure 8 and Figure 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⁻¹) 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⁻¹, as opposed to a mean of 110.1 t·ha⁻¹ 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⁻¹ or +14% (M3K) and −22 t·ha⁻¹ 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).

5. 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² = 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.