Impacts of Airborne Lidar Pulse Density on Estimating Biomass Stocks and Changes in a Selectively Logged Tropical Forest

Airborne lidar is a technology well-suited for mapping many forest attributes, including aboveground biomass (AGB) stocks and changes in selective logging in tropical forests. However, trade-offs still exist between lidar pulse density and accuracy of AGB estimates. We assessed the impacts of lidar pulse density on the estimation of AGB stocks and changes using airborne lidar and field plot data in a selectively logged tropical forest located near Paragominas, Pará, Brazil. Field-derived AGB was computed at 85 square 50 × 50 m plots in 2014. Lidar data were acquired in 2012 and 2014, and for each dataset the pulse density was subsampled from its original density of 13.8 and 37.5 pulses·m−2 to lower densities of 12, 10, 8, 6, 4, 2, 0.8, 0.6, 0.4 and 0.2 pulses·m−2. For each pulse density dataset, a power-law model was developed to estimate AGB stocks from lidar-derived mean height and corresponding changes between the years 2012 and 2014. We found that AGB change estimates at the plot level were only slightly affected by pulse density. However, at the landscape level we observed differences in estimated AGB change of >20 Mg·ha−1 when pulse density decreased from 12 to 0.2 pulses·m−2. The effects of pulse density were more pronounced in areas of steep slope, especially when the digital terrain models (DTMs) used in the lidar derived forest height were created from reduced pulse density data. In particular, when the DTM from high pulse density in 2014 was used to derive the forest height from both years, the effects on forest height and the estimated AGB stock and changes did not exceed 20 Mg·ha−1. The results suggest that AGB change can be monitored in selective logging in tropical forests with reasonable accuracy and low cost with low pulse density lidar surveys if a baseline high-quality DTM is available from at least one lidar survey. We recommend the results of this study to be considered in developing projects and national level MRV systems for REDD+ emission reduction programs for tropical forests.


Introduction
The Amazon is the largest remaining tropical forest in the world, however, its original extent has been steadily reduced due to deforestation and forest degradation, although deforestation rates in Brazil have decreased by 70% since 2004 [1,2]. In recent decades, selective logging of valuable tree species has been an important land use of tropical forest in the Brazilian Amazon [3,4]. Selective logging timber extraction removes only the most valuable tree species from the forest [1]. It contributes substantially to gross carbon fluxes from the Brazilian Amazon and in other tropical regions as well [5]. Selective logging has continued apace with degradation from forest fires and forest fragmentation, and may also degrade the Amazon forest through long term changes in structure, loss of forest carbon and species diversity [6]. Characterizing the spatial distribution of forest structure, aboveground biomass (AGB), and AGB changes are important prerequisites for understanding carbon cycle dynamics and for monitoring the impact of selective logging in tropical forests over time [7]. Accurate, landscape-wide estimates of AGB stocks and changes from selective logging in tropical forest are also desired for ongoing climate mitigation efforts to Reduce Emissions from Deforestation and Forest Degradation (REDD+) and for Measuring, Report and Verification (MRV) systems [7,8].
Airborne lidar is a technology well-suited to measure forest structure and estimate AGB stocks and changes in tropical forests (e.g., [8][9][10][11][12][13]). Lidar can provide high resolution, three-dimensional information on forest structure and the underlying topography [14,15]. Recently, d'Oliveira et al. [7] and Andersen et al. [8] have used airborne lidar for detecting selective logging activities and mapping AGB stocks and changes in the Brazilian Amazon forest. To monitor selective logging impacts on forest structure and AGB changes over large areas, multiple lidar data acquisitions must be acquired that increase the cost of data collection and processing over time. Many factors influence the cost of lidar data, including variables such as project size, location and deliverables, as well as market variables, such as competition amongst lidar vendors. A major variable that affects the cost of acquisition of lidar data is the pulse density [16], defined as the number of pulses emitted by the sensor per m 2 (pulses·m −2 ) (Evans et al., 2009). As pulse density increases, so does acquisition cost, due to the direct link between pulse density, aircraft altitude and flight time [17,18].
While airborne lidar can facilitate timely and accurate estimates of forest structure in tropical forest, trade-offs still exist between lidar pulse density and accuracy. For instance, it is unclear how much the lidar pulse density can be reduced and still maintain an adequate level of accuracy for AGB change estimation in tropical forests. Leitold et al. [19] and Ota et al. [20] have carried out studies to examine the relationship between lidar pulse density and AGB stock estimation accuracy in tropical forests. Even though they found that AGB can be accurately estimated from lidar using low-pulse density, neither of these authors assessed the impact of pulse density on estimating and mapping AGB stocks and AGB change at landscape level, and in the context of selectively logged tropical forest. Here, we focus on the impacts of pulse density in estimating AGB change in tropical forests and provide recommendations for specification of lidar data acquisitions for forest monitoring, REDD+ projects and MRV systems. We work with data collected in an eastern selectively logged Amazonian forest and lidar data with high pulse density acquired in 2012 and 2014. The study quantifies how reduced pulse density reduces the accuracy of AGB stock and AGB change estimation at plot and landscape levels. We evaluate accuracy based on lidar data acquired over dense tropical forests with variations of terrain characteristics and topography. The results of our study are then discussed in the context of implementation of airborne lidar systems in monitoring forest AGB change for REDD+ and emission reduction programs.

Study Area
The study was conducted at the Fazenda Cauaxi in the Paragominas Municipality of Pará State, Brazil, in the eastern Amazon ( Figure 1). The climate of the Cauaxi region is typically humid, with an Remote Sens. 2017, 9, 1068 3 of 19 average annual temperature of about 25 • C and an average annual precipitation of 2200 mm, which primarily falls between the months of January and June [21]. The forest is predominantly classified as tropical dense moist forest (IBGE, 2004). The terrain ranges from flat to steep slopes and the soils within the region are classified predominately as dystrophic yellow latosols following the Brazilian classification system [22]. The study area is divided into 12 logging units, where ten of them have been logged through the reduced-impact logging (RIL) since 2007 and the remaining two are still unlogged ( Figure 1D). Remote Sens. 2017, 8, 1068 3 of 19 primarily falls between the months of January and June [21]. The forest is predominantly classified as tropical dense moist forest (IBGE, 2004). The terrain ranges from flat to steep slopes and the soils within the region are classified predominately as dystrophic yellow latosols following the Brazilian classification system [22]. The study area is divided into 12 logging units, where ten of them have been logged through the reduced-impact logging (RIL) since 2007 and the remaining two are still unlogged ( Figure 1D).

Field Data Collection
A total of 22 field transects of 20 × 500 m were stratified randomly across the study area in 2012, and 88 plots of 50 × 50 m (0.25 ha) were spaced at 100 m intervals along the transects in 2014. In the field, plot corners were registered using differential GNSS (GeoXH6000, Trimble Navigation, Ltd., Dayton, OH, USA). At each plot, a sub-plot (strip) along one side of the plot with dimensions of 5 × 50 m (250 m 2 ) was also demarcated. Because three of the plots were not covered by the lidar data, we selected 85 plots for further analysis. For each plot, all living trees with diameter at breast height (dbh) ≥ 35 cm were identified by parataxonomists familiar with the flora of the region and their dbh measured. In the sub-plots, all trees with dbh ≥ 10 cm were measured. Dbh was measured at 1.3 m above ground or above buttresses. A total of 1757 living trees were measured. The AGB (kg) of each tree was estimated using the Chave et al.

Field Data Collection
A total of 22 field transects of 20 × 500 m were stratified randomly across the study area in 2012, and 88 plots of 50 × 50 m (0.25 ha) were spaced at 100 m intervals along the transects in 2014. In the field, plot corners were registered using differential GNSS (GeoXH6000, Trimble Navigation, Ltd., Dayton, OH, USA). At each plot, a sub-plot (strip) along one side of the plot with dimensions of 5 × 50 m (250 m 2 ) was also demarcated. Because three of the plots were not covered by the lidar data, we selected 85 plots for further analysis. For each plot, all living trees with diameter at breast height (dbh) ≥ 35 cm were identified by parataxonomists familiar with the flora of the region and their dbh measured. In the sub-plots, all trees with dbh ≥ 10 cm were measured. Dbh was measured at 1.3 m above ground or above buttresses. A total of 1757 living trees were measured. The AGB (kg) of each tree was estimated using the Chave et al. [23]  where AGB (kg) is the live tree aboveground biomass in Kg; dbh is the diameter at breast height (1.30 m); ρ is the wood density and E is a measure of environmental stress. In this study area location E = −0.103815. The total live aboveground biomass (AGB) of the plots and sub-plots was obtained by aggregating the individual tree biomass values and converting to Mg·ha −1 . The summary of the dbh and AGB measurements in 2014 at the sample plots is presented in Table 1.

Lidar Data Acquisition and Processing
Airborne lidar data were collected as part of Sustainable Landscapes Brazil, a joint project of the Brazilian Corporation of Agricultural Research Corporation (EMBRAPA) and the United States Forest Service. The first lidar acquisition occurred in July of 2012 with a pulse density of 13.8 pulses·m −2 , with the second lidar collection in December of 2014 with a pulse density of 37.5 pulses·m −2 . The total area covered by the lidar survey for the multitemporal analysis was 1200 ha. The data attributes of the lidar sensor and flight characteristics are listed in Table 2. The lidar data processing can be summarized in three steps: (i) Lidar data thinning: Lidar data from both 2012 and 2014 were thinned from the original pulse density of 13.8 and 37.5 pulses·m −2 to lower densities of 12, 10, 8, 6, 4, 2, 0.8, 0.6, 0.4 and 0.2 pulses·m −2 . An example of the pulse density reduction is shown in Figure 2. The reduction of pulse density was executed using the algorithm implemented in the ThinData utility of the FUSION toolkit [24]. The algorithm first identifies all the returns that belong to a pulse, and randomly reduces the number of pulses until achieving the desired pulse density within a certain grid cell size (e.g., in this case 50 m). To evaluate the uncertainty in the thinning process, we generated 30 random replicates for each lidar dataset and target density. (ii) Digital terrain models and lidar data above ground height normalization: Ground returns were classified using the Progressive Triangulated Irregular Network (TIN) densification algorithm implemented in lasground [25] (settings: step is 10 m, bulge is 0.5 m, spike is 1 m, offset is 0.05 m), and 1 m DTMs were created for each of the reduced pulse density datasets using the blast2dem utility in Lastools [25]. Afterwards, the lidar datasets were normalized to height aboveground by subtraction of the DTM elevation from the Z coordinate of each return projected above the ground. In particular, the impacts of pulse density were assessed under two DTM scenarios. DTM scenario DS2, where to simulate the impact of pulse density on subsequent acquisition, the DTMs generated at the highest pulse density (37.5 pulse·m −2 ) from 2014 were used to normalize the lidar datasets from both 2012 and 2014. (iii) Lidar-derived Mean Height (HMEAN): In this study, HMEAN, the mean height of all returns above 1.3 m in height, was computed at plot and landscape levels. Herein, even though the plot corners were geolocated with differential GNSS to within <1 m in most cases (Longo et al., 2016), we chose to optimize plot location to reflect canopy conditions. It is common to find trees with crown diameters larger than 30 m in the study area. When large tree stems are found outside of a plot, substantial proportions of their crowns may fall inside a plot thereby influencing lidar metrics such as HMEAN. To avoid these effects, we iteratively shifted each plot within a 25 m square neighbourhood on the lidar canopy height model (CHM) for improved co-registration, and consequently to achieve a better correlation between AGB and HMEAN. After this procedure, we observed that the plot centers were moved an average of 18.38 m (sd ± 6.32 m) from the initial plot locations. Finally, at the landscape level, HMEAN was computed in grids with a cell size of 50 m × 50 m. (iii) Lidar-derived Mean Height (HMEAN): In this study, HMEAN, the mean height of all returns above 1.3 m in height, was computed at plot and landscape levels. Herein, even though the plot corners were geolocated with differential GNSS to within < 1 m in most cases (Longo et al., 2016), we chose to optimize plot location to reflect canopy conditions. It is common to find trees with crown diameters larger than 30 m in the study area. When large tree stems are found outside of a plot, substantial proportions of their crowns may fall inside a plot thereby influencing lidar metrics such as HMEAN. To avoid these effects, we iteratively shifted each plot within a 25 m square neighbourhood on the lidar canopy height model (CHM) for improved co-registration, and consequently to achieve a better correlation between AGB and HMEAN. After this procedure, we observed that the plot centers were moved an average of 18.38 m (sd ± 6.32 m) from the initial plot locations. Finally, at the landscape level, HMEAN was computed in grids with a cell size of 50 m × 50 m.

Aboveground Biomass Change Estimation and Mapping
For both DTM scenarios described in Section 2.3, we used the nls function in R [26] to calibrate the relationship between plot level AGB measurements and lidar-derived HMEAN. Non-linear leastsquares regression models (the power-law models) were used to model AGB across 30 replicates of each reduced target density in 2014, and used to predict and map AGB stocks in 2012 (AGB ) and 2014 (AGB ) at plot and landscape levels. The AGB change (∆ ( ) ) estimation was then computed as the difference in AGB prediction from 2012 to 2014.
where and are the estimates' parameters of the power-law models in 2014. Leave-oneout cross-validation (LOOCV) was developed (e.g., [14,15]), and the prediction precision of the LOOCV models was evaluated in terms of the coefficient of determination (R 2 ), absolute and relative Root Mean Square Error (RMSE), and absolute and relative bias from the linear relationship between observed and LOOCV predicted AGB values:

Aboveground Biomass Change Estimation and Mapping
For both DTM scenarios described in Section 2.3, we used the nls function in R [26] to calibrate the relationship between plot level AGB measurements and lidar-derived HMEAN. Non-linear least-squares regression models (the power-law models) were used to model AGB across Leave-one-out cross-validation (LOOCV) was developed (e.g., [14,15]), and the prediction precision of the LOOCV models was evaluated in terms of the coefficient of determination (R 2 ), absolute and relative Root Mean Square Error (RMSE), and absolute and relative bias from the linear relationship between observed and LOOCV predicted AGB values: where n is the number of plots, y i is the observed value for plot i, and y i is the predicted value for plot i. Moreover, relative RMSE and bias were calculated by dividing absolute RMSE and bias (Equations (5) and (6)) by observed AGB mean. In order to have prediction accuracy equal to or higher than a conventional forest inventory in tropical forest, we defined accepted model accuracy as relative RMSE and bias of ≤20%. At the landscape level, maps representing the mean and standard deviation of AGB stocks and AGB changes from the 30 replications were calculated and used as final estimates to assess the impact of pulse density. The mean of AGB stocks in 2012 and 2014 and AGB changes from the final maps were computed as follows: where An uncertainty analysis of the AGB 2012 and AGB 2014 at landscape level for each pulse density target and DTM scenario was also performed by integrating the pixel level errors and accounting for spatial autocorrelation of the errors as follows [27][28][29]: where σ 2 AGB is the variance of the estimator for the mean AGB stock for the entire study area; m is the number of pixels in the area; ρ(d) is the spatial autocorrelation function of the distance, d, based on an exponential semi-variogram model; and σ i is the estimated standard error of AGB stock values at the i-th pixel.
The variance of the estimator of the mean σ 2 ∆AGB at landscape level was computed as: where the variances ( σ 2 AGB(2012) and σ 2 AGB(2014) ) were computed as described in Equation (10) and the cross-time covariance of the AGB was computed according to McRoberts et al. [30]: Remote Sens. 2017, 9, 1068 7 of 19

Assessing Effects of Pulse Density on Lidar-Derived Mean Canopy Height
As the HMEAN values vary from one replication to another for a given plot and target density, we calculated the mean and standard deviation of HMEAN at the plot level across the 30 repetitions. The impacts of pulse density on HMEAN at plot level was then evaluated by the reliability ratio, which is the ratio of the variance of HMEAN among sample plots, to the total variance of the HMEAN across 30 repetitions (Fuller, 1987): where S 2 u is the estimated among-plot variance of HMEAN and S 2 w is the estimated average within-plot variance. Reliability ratio ranges from 0 (no reliability) to 1 (complete reliability), and large replication variance makes HMEAN a low reliability predictor of AGB.

Assessing Impacts of Pulse Density on the Aboveground Biomass Stocks and Change Estimation at the Plot and Landscape Levels
The impacts of the lidar pulse density on AGB stocks and AGB change estimations were assessed at the plot level by comparing the R 2 , relative and absolute RMSE and bias across pulse densities and DTM scenarios. Boxplots were created to compare the variability of LOOCV AGB stock estimates, RMSE and bias at plot level. Because we built a model for each repeated dataset and target density, we also calculated the mean and standard deviation of the a 2014 and b 2014 model parameters. The two-sided Kolmogorov-Smirnov (KS) [26] test was used to test if the distributions of the HMEAN, observed and LOOCV AGB stock estimates in 2014 differed significantly, with significance level of 0.05, across pulse densities and DTM scenarios. Moreover, besides the maps of AGB stock and changes, we created a map of the slope over the study area, and the impacts of the lidar pulse density on AGB stocks and change estimations were evaluated by the difference on the AGB estimates from 12 pulses·m −2 to lower pulse densities, in each DTM scenario and across slope gradients of 0 to 12%; 12 to 24% and 24 to 35%. An overview of the methodology is outlined in Figure 3. Reliability ratio (RR) = + where is the estimated among-plot variance of HMEAN and is the estimated average withinplot variance. Reliability ratio ranges from 0 (no reliability) to 1 (complete reliability), and large replication variance makes HMEAN a low reliability predictor of AGB.

Assessing Impacts of Pulse Density on the Aboveground Biomass Stocks and Change Estimation at the Plot and Landscape Levels
The impacts of the lidar pulse density on AGB stocks and AGB change estimations were assessed at the plot level by comparing the R 2 , relative and absolute RMSE and bias across pulse densities and DTM scenarios. Boxplots were created to compare the variability of LOOCV AGB stock estimates, RMSE and bias at plot level. Because we built a model for each repeated dataset and target density, we also calculated the mean and standard deviation of the and model parameters. The two-sided Kolmogorov-Smirnov (KS) [26] test was used to test if the distributions of the HMEAN, observed and LOOCV AGB stock estimates in 2014 differed significantly, with significance level of 0.05, across pulse densities and DTM scenarios. Moreover, besides the maps of AGB stock and changes, we created a map of the slope over the study area, and the impacts of the lidar pulse density on AGB stocks and change estimations were evaluated by the difference on the AGB estimates from 12 pulses·m −2 to lower pulse densities, in each DTM scenario and across slope gradients of 0 to 12%; 12 to 24% and 24 to 35%. An overview of the methodology is outlined in Figure 3.

Impacts of Lidar Pulse Density on Mean Canopy Height (HMEAN)
The impacts of the pulse density on the lidar derived HMEAN at plot level are shown in

Impacts of Lidar Pulse Density on Mean Canopy Height (HMEAN)
The impacts of the pulse density on the lidar derived HMEAN at plot level are shown in  (Figure 4(a1,a2)). The variability of HMEAN represented by the standard deviation (Figure 4(b1,b2)) within plots for the 30 replications was larger at lower than at higher pulse densities for both DTM scenarios and years. Reduced pulse density resulted in a decreased reliability ratio (RR); however, RR still showed high stability of HMEAN across DTM scenarios and years with RR values higher than 0.96 (Figure 4(c1,c2)).  (Figure 4(a1,a2)). The variability of HMEAN represented by the standard deviation ( Figure 4(b1,b2)) within plots for the 30 replications was larger at lower than at higher pulse densities for both DTM scenarios and years. Reduced pulse density resulted in a decreased reliability ratio (RR); however, RR still showed high stability of HMEAN across DTM scenarios and years with RR values higher than 0.96 (Figure 4(c1,c2)).

Effects of Lidar Pulse Density on AGB Modelling
The parameters and of the models adjusted in 2014 to predict AGB from lidarderived HMEAN for 2014 across scenarios and pulse densities are presented in Table 3. Reduced pulse density resulted in increased variation of these parameter values in both scenarios, where DS2 showed less variation then DS1. For DS1 the mean of model parameters and showed significant differences (KS: D ≥ 0.40, p-values ≤ 0.015) at pulse densities ranging from 0.2 to 4 pulses·m −2 compared to the value at 12 pulses·m −2 ; for DS2, significant differences from 12 pulses·m −2 (KS: D ≥ 0.46, p-values < 0.002) were found only at pulse densities ranging from 0.2 to 2 pulse·m −2 . Both parameters and showed significant differences (KS: D ≥ 0.47, p-values ≤ 0.002) at pulse densities ranging from 0.2 to 2 pulse·m −2 when compared at the same pulse density, but across scenarios.

Effects of Lidar Pulse Density on AGB Modelling
The parameters a 2014 and b 2014 of the models adjusted in 2014 to predict AGB from lidar-derived HMEAN for 2014 across scenarios and pulse densities are presented in Table 3. Reduced pulse density resulted in increased variation of these parameter values in both scenarios, where DS2 showed less variation then DS1. For DS1 the mean of model parameters a 2014 and b 2014 showed significant differences (KS: D ≥ 0.40, p-values ≤ 0.015) at pulse densities ranging from 0.2 to 4 pulses·m −2 compared to the value at 12 pulses·m −2 ; for DS2, significant differences from 12 pulses·m −2 (KS: D ≥ 0.46, p-values < 0.002) were found only at pulse densities ranging from 0.2 to 2 pulse·m −2 . Both parameters a 2014 and b 2014 showed significant differences (KS: D ≥ 0.47, p-values ≤ 0.002) at pulse densities ranging from 0.2 to 2 pulse·m −2 when compared at the same pulse density, but across scenarios. The performance of the models to estimate AGB stock in 2014 was further assessed by leave-one-out cross-validation. The HMEAN was an important lidar metric to explain AGB stock variation. Reduced pulse density resulted in decreased R 2 and increased relative and absolute RMSE and bias for both DTM scenarios, particularly for DS1 ( Figure 5). Mean R 2 values ranged from 0.60 (±0.09) to 0.73 (±0.00) and 0.71 (±0.02) to 0.73 (±0.00) across pulse densities in DS1 and DS2, respectively, and showed significant differences (KS: D ≥ 0.37, p-value ≤ 0.034) at pulse densities ranging from 0.2 to 6 pulse·m −2 and from 0.2 to 2 pulse·m −2 compared to the results at 12 pulses·m −2 in DS1 and DS2, respectively (Figure 5a). Mean relative RMSE values ranged from 18.81 (±0.06) to 22.80 (±0.270) % and 18.81 (±0.06) to 19.50 (±0.69) %, across pulse densities in DS1 and DS2, respectively, and also showed significant differences (KS: D ≥ 0.37, p-values < 0.035) at pulse densities ranging from 0.2 to 6 pulse·m −2 and from 0.2 to 2 pulse·m −2 compared to the value at 12 pulses·m −2 in DS1 and DS2, respectively ( Figure 5(b1)). Bias was less affected by pulse density in DS2 than in DS1 ( Figure 5(c1)). While mean bias across pulse densities ranged from −0.48 (±0.38) to −0.04 (±0.01) % in DS1, and showed significant differences (KS: D ≥ 0.57, p-values < 0.001) at pulse densities ranging from 0.2 to 2 pulse·m −2 compared to at 12 pulses·m −2 , mean bias remained constant around 0.04% in DS2, and showed significant differences (KS: D ≥ 0.43, p-values < 0.006) only at pulse densities ranging from 0.2 to 0.6 pulse·m −2 when compared to 12 pulses·m −2 .
The estimated AGB change among plots ranged from −163.845 to 354.29 Mg·ha −1 and from −159.79 to 155.46 Mg·ha −1 across pulse densities for DS1 and DS2, respectively. The estimated mean

Impacts of Lidar Pulse Density on AGB Stocks and AGB Change Estimation at the Plot and Landscape Levels
The AGB stocks estimates from the leave-one-out cross-validation at plot level are presented in Figure 6. In 2012, mean predicted AGB ranged from 226.81 (±76.29) to 230.35 (±75.98) Mg·ha −1 and from 226.51 (±75.78) to 227.02 (±74.46) Mg·ha −1 across pulse densities for DS1 and DS2, respectively, and significant differences (KS: D ≥ 0.035, p-value ≤ 0.03) were found at pulse densities ranging from 0.2 to 0.4 pulse·m −2 compared to the value at 12 pulses·m −2 , but only in DS1. In 2014, mean predicted AGB ranged from 237.83 (±74.37) to 238.89 (±70.21) Mg·ha −1 and 273.83 (±74.28) to 273.84 (±73.43) Mg·ha −1 across pulse densities for DS1 and DS2, respectively, and significant differences (KS: D ≥ 0.04, p-value ≤ 0.015) were found only at pulse densities ranging from 0.2 to 0.4 pulse· m −2 compared to the value at 12 pulses·m −2 , but only in DS1 as well (Figure 6(a1,a2)).  11.01 (±38.68) Mg·ha −1 across pulse densities for DS1 and DS2, respectively, and significant differences (KS: D ≥ 0.38, p-values ≤ 0.042) were found at pulse densities ranging from 0.2 to 0.8 pulse· m −2 and from 0.2 to 0.6 pulse·m −2 compared to 12 pulses·m −2 in DS1 and DS2, respectively. The reduced pulse density increased the variance of the AGB estimates within plots across the 30 replicates, and showed significant differences (KS: D ≥ 0.59, p-values < 0.001) in the standard deviation of AGB stocks and changes, from 0.2 to 0.8 pulse·m −2 compared to the value at 12 pulses·m −2 . DS2 shows slightly less variation than DS1, but after pulse density reaches values higher than 0.8 pulses·m −2 , both scenarios show very low and not significant differences in standard deviation of AGB stocks and changes within plots across replicates.  11.01 (±38.68) Mg·ha −1 across pulse densities for DS1 and DS2, respectively, and significant differences (KS: D ≥ 0.38, p-values ≤ 0.042) were found at pulse densities ranging from 0.2 to 0.8 pulse·m −2 and from 0.2 to 0.6 pulse·m −2 compared to 12 pulses·m −2 in DS1 and DS2, respectively. The reduced pulse density increased the variance of the AGB estimates within plots across the 30 replicates, and showed significant differences (KS: D ≥ 0.59, p-values < 0.001) in the standard deviation of AGB stocks and changes, from 0.2 to 0.8 pulse·m −2 compared to the value at 12 pulses·m −2 . DS2 shows slightly less variation than DS1, but after pulse density reaches values higher than 0.8 pulses·m −2 , both scenarios show very low and not significant differences in standard deviation of AGB stocks and changes within plots across replicates.  Table 4.  Table 4.    Landscape-wide standard deviation of estimated AGB change was also mapped at a 50 m × 50 m grid cell resolution (Figure 8c,d and Figure S4); reduced pulse density increased variation in estimated AGB change within replicates, and showed large standard deviation of AGB at 0.2 pulses·m −2 , for both DS1 and DS2 (Figure 8(c1,d1)). In general, the standard deviation of estimated AGB change was higher in DS1 than DS2. Landscape-wide elevation and slope were also mapped (Figure 8a,b), and for both DS1 and DS2 the large variability in AGB change occurred in high slope areas, reaching up to 33 Mg· ha −1 in areas with slopes ranging from 24 to 36% (Figure 8(c1)). Estimated AGB change at landscape level at pulse densities ranging from 0.2 to 10 pulses·m −2 were compared with those estimates at 12 pulses·m 2 across slopes and DTM scenarios (Figure 9a-c). In general, when compared with the AGB change estimates from 12 pulses·m −2 , reduced pulse density underestimated AGB change at landscape level in areas with slopes higher than 12% and showed significant differences (KS: D ≥ 0.45, p-value > 0.001) at pulse densities ranging from 0.2 to 0.8 pulses·m −2 , but only for DS1 at slopes of 12-24% (Figure 9b) and 24-35% (Figure 9c). At slopes of 0-12% (Figure 9a), reduced pulse density increased the difference in AGB change when compared with those estimates at 12 pulses·m 2 in both DS1 and DS2, but their differences were no higher than 20 Mg· ha −1 (Figure 9a). Landscape-wide standard deviation of estimated AGB change was also mapped at a 50 m × 50 m grid cell resolution (Figure 8c,d and Figure S4); reduced pulse density increased variation in estimated AGB change within replicates, and showed large standard deviation of AGB at 0.2 pulses·m −2 , for both DS1 and DS2 (Figure 8(c1,d1)). In general, the standard deviation of estimated AGB change was higher in DS1 than DS2. Landscape-wide elevation and slope were also mapped (Figure 8a,b), and for both DS1 and DS2 the large variability in AGB change occurred in high slope areas, reaching up to 33 Mg·ha −1 in areas with slopes ranging from 24 to 36% (Figure 8(c1)). Estimated AGB change at landscape level at pulse densities ranging from 0.2 to 10 pulses·m −2 were compared with those estimates at 12 pulses·m 2 across slopes and DTM scenarios (Figure 9a-c). In general, when compared with the AGB change estimates from 12 pulses·m −2 , reduced pulse density underestimated AGB change at landscape level in areas with slopes higher than 12% and showed significant differences (KS: D ≥ 0.45, p-value > 0.001) at pulse densities ranging from 0.2 to 0.8 pulses·m −2 , but only for DS1 at slopes of 12-24% (Figure 9b) and 24-35% (Figure 9c). At slopes of 0-12% (Figure 9a), reduced pulse density increased the difference in AGB change when compared with those estimates at 12 pulses·m 2 in both DS1 and DS2, but their differences were no higher than 20 Mg·ha −1 (Figure 9a).

Discussion
This research assessed the impact of lidar pulse density on AGB stocks and changes estimations in an Amazon tropical forest. Previous studies have evaluated the impact of lidar pulse density on forest attribute estimation from lidar data (e.g., [31][32][33]), yet few studies have evaluated the impacts of lidar pulse density on forest AGB stock estimates in tropical forest [19,20,34]. To our knowledge, this is the first study to assess the impact of airborne lidar pulse density on AGB stocks and AGB change estimations in tropical forest, and in the context of using an airborne lidar system in selective logging for monitoring forest AGB change for REDD+ and emission reduction programs.
Many lidar-derived metrics have been used for modelling forest attributes [7,[35][36][37][38][39][40][41][42]. Hansen et al. [31] evaluated the effects of lidar pulse density on DTM and canopy structure metrics in a tropical forest, and showed also that HMEAN was one of the most stable predictor variables for modelling forest attributes using airborne lidar data. In our study, under both DTM scenarios, reduced pulse density did not significantly affect the variability of HMEAN among plots. Magnusson et al. [32] recommended a calibration of lidar models when the reliability ratio (RR) of one or more predictors is below 0.9. In our study, even though reduced pulse density increased the standard deviation of HMEAN within plots, the RR of HMEAN across all pulse densities and both DTM scenarios remained very high (RR > 0.9). Therefore, further calibration of models was not necessary. This is not surprising given that HMEAN, unlike other metrics such as the top mean canopy height (MCH), considers all returns above a certain height threshold (e.g., 1.37 m) to compute a vertical mean height, and therefore uses more information to describe the canopy structure. Unlike MCH, HMEAN is computed directly from the point cloud, and does not have any influence of CHM interpolation methods or grid cell size. Garcia et al. [43] assessed the impact of lidar point density on the prediction of AGB across different forest ecosystems, and found that predictions were more affected when using CHM-derived metrics then those computed directly from the 3D point cloud, even if the point density was as low as 1 point m −2 . Herein, we simulated low pulse density lidar datasets by removing pulses randomly, however, other approaches of removing lidar pulses have also been found in the literature [16][17][18][19]32,44] and may lead to different outcomes when considering the covarying effect of pulse density when survey parameters are changed. Nevertheless, our results on HMEAN variation patterns agree with a previous study [31], and independent of the approach used, a realistic thinning approach on real lidar data is always extremely challenging [45].
In addition to the high stability, HMEAN correlated well with AGB, and we obtained linear relationships between observed and predicted AGB via LOOCV that explained at least 40% of the variation at the lowest pulse density (0.2 pulses·m −2 ) in DS1, for example. For comparison purposes, the R 2 values were substantially greater than adjusted R 2 = 0.43 obtained by Leitold et al. [19] using a regression model for predicting AGB stocks in the Brazilian Atlantic forest, but similar to those

Discussion
This research assessed the impact of lidar pulse density on AGB stocks and changes estimations in an Amazon tropical forest. Previous studies have evaluated the impact of lidar pulse density on forest attribute estimation from lidar data (e.g., [31][32][33]), yet few studies have evaluated the impacts of lidar pulse density on forest AGB stock estimates in tropical forest [19,20,34]. To our knowledge, this is the first study to assess the impact of airborne lidar pulse density on AGB stocks and AGB change estimations in tropical forest, and in the context of using an airborne lidar system in selective logging for monitoring forest AGB change for REDD+ and emission reduction programs.
Many lidar-derived metrics have been used for modelling forest attributes [7,[35][36][37][38][39][40][41][42]. Hansen et al. [31] evaluated the effects of lidar pulse density on DTM and canopy structure metrics in a tropical forest, and showed also that HMEAN was one of the most stable predictor variables for modelling forest attributes using airborne lidar data. In our study, under both DTM scenarios, reduced pulse density did not significantly affect the variability of HMEAN among plots. Magnusson et al. [32] recommended a calibration of lidar models when the reliability ratio (RR) of one or more predictors is below 0.9. In our study, even though reduced pulse density increased the standard deviation of HMEAN within plots, the RR of HMEAN across all pulse densities and both DTM scenarios remained very high (RR > 0.9). Therefore, further calibration of models was not necessary. This is not surprising given that HMEAN, unlike other metrics such as the top mean canopy height (MCH), considers all returns above a certain height threshold (e.g., 1.37 m) to compute a vertical mean height, and therefore uses more information to describe the canopy structure. Unlike MCH, HMEAN is computed directly from the point cloud, and does not have any influence of CHM interpolation methods or grid cell size. Garcia et al. [43] assessed the impact of lidar point density on the prediction of AGB across different forest ecosystems, and found that predictions were more affected when using CHM-derived metrics then those computed directly from the 3D point cloud, even if the point density was as low as 1 point m −2 . Herein, we simulated low pulse density lidar datasets by removing pulses randomly, however, other approaches of removing lidar pulses have also been found in the literature [16][17][18][19]32,44] and may lead to different outcomes when considering the covarying effect of pulse density when survey parameters are changed. Nevertheless, our results on HMEAN variation patterns agree with a previous study [31], and independent of the approach used, a realistic thinning approach on real lidar data is always extremely challenging [45].
In addition to the high stability, HMEAN correlated well with AGB, and we obtained linear relationships between observed and predicted AGB via LOOCV that explained at least 40% of the variation at the lowest pulse density (0.2 pulses·m −2 ) in DS1, for example. For comparison purposes, the R 2 values were substantially greater than adjusted R 2 = 0.43 obtained by Leitold et al. [19] using a regression model for predicting AGB stocks in the Brazilian Atlantic forest, but similar to those reported by d'Oliveira et al. [7] and Andersen et al. [8] using lidar for detecting and quantifying AGB changes in selective logging in western Amazonia, respectively. Herein, a substantial decrease in R 2 and increase in RMSE and bias occurred when the pulse density was reduced from 12 pulses·m −2 to pulse densities lower than 2 pulses·m −2 . For instance, AGB models under DS1 were much more affected at pulse densities ≤2 pulses·m −2 . Magnusson et al. [32] and Watt et al. [46] in mixed conifer forests also evaluated the effect of lidar pulse density on the prediction accuracy of forest attributes under DS1 and DS2, and found increased RMSE at relatively low pulse densities. However, similar to findings by those authors, our predictions were more affected in DS1 than DS2. Nevertheless, in both scenarios, AGB stocks were underestimated. An underestimation of AGB stock with reduced pulse density was also found in Leitold et al. [22], but as AGB stocks at low pulse densities were estimated from a single model adjusted at high pulse density, the author attributed the underestimation to a systematic error in the DTM propagated to the canopy. Herein, the underestimation in AGB stock in DS1 is attributed to the deterioration of the DTM and HMEAN quality due to combined random effects derived from pulse density reduction, while in DS2, the underestimation in AGB stock is only attributed to the poor HMEAN quality. Moreover, as we kept a constant interpolation method and grid size for DTM modeling in DS1 and DS2 as in Magnusson et al. [32] and Watt et al. [44], the pulse density is the only factor affecting the DTM and HMEAN quality. Therefore, the differences in R 2 , RMSE and bias between DS1 and DS2 were only attributable to the deteriorating DTM and HMEAN quality.
In this study, we mapped AGB change at the landscape level at a spatial resolution of 50 m. Therefore, unlike most previous studies (e.g., [31][32][33]), we not only evaluated the impact of pulse density on AGB stocks estimation at plot the level, but also at the landscape scale. Andersen et al. [8] used repeated lidar flights to monitor selective logging in western Amazonia, and AGB stock and changes maps were accurately derived from lidar data acquired in 2010 and 2011 with pulse densities of 25 pulses·m −2 and 14 pulses·m −2 , respectively. In the study, the authors found that multi-temporal lidar data can be used to detect and quantify AGB changes due to selective logging activities, even when the level of AGB change is low (10-20 Mg·ha −1 ). While the pulse density was appropriate for the study, they are not economically feasible for large-area acquisitions and for monitoring selective logging over time. Hansen et al. [31] suggest that canopy metrics derived from sparse pulse density ALS data can be used for AGB estimation in a tropical forest; however, the authors either estimated AGB or expanded the analysis to landscape level in their studies. Wilkes et al. [47] found that structural metrics (canopy height, canopy cover and vertical canopy structure) derived from pulse densities < 0.5 pulses·m −2 returned larger differences, particularly for tropical forest. Herein, while our uncertainty analysis showed that reduced pulse density did not affect the accuracy of mean estimated AGB change at the landscape level in both DTM scenarios, reduced pulse density did significantly affect the standard deviation of estimated AGB change at the local scale, especially in areas of steeper slopes. However, this effect was only observed for the DS1 scenario at pulse densities ≤ 0.8 pulses·m −2 . Variation in AGB estimates increased with decreasing pulse density; this is illustrated by mapped standard deviation at landscape level (Figure 8).
From a carbon monitoring perspective, our results show that it is not necessary to acquire lidar data with pulse density higher than 0.2 pulses·m −2 for mapping AGB stocks and AGB change in selective logging areas in the Amazonia when a DTM is already available from previously flown high pulse density airborne lidar. In cases where no such DTM already exists, our results suggest that a lidar data acquisition with a minimum pulse density of 2 pulses·m −2 is necessary. Furthermore, we have demonstrated that it is possible to capture AGB stocks and AGB change variability at the stand level even in steep terrain with low pulse density lidar under DS2 and with pulse density equal to or higher than 2 pulses·m −2 in DS1. Despite the change in AGB stocks, selective logging can also substantially alter forest structure and affect tree survival, growth, and recruitment rates for up to a decade post-harvest [48]. In this study, we did not evaluate the impact of pulse density for detecting forest impacts associated with selective logging, nor the combined effect of pulse density and plot size on AGB change estimation. Therefore, we suggest caution when acquiring new lidar datasets, because the accuracy of the AGB change estimates in selective logging areas may depend also on other factors, such as plot and grid cell sizes for sampling and mapping, which were not evaluated in this study. In some small, randomly distributed areas, AGB change was highly overestimated (AGB change ≥ 100 Mg ha −1 ), which is not biologically possible in only two years. These overestimates could result from subtracting AGB stocks predicted from models calibrated with only 2014 data, and may not have resulted had independent AGB models been calibrated with data from both years (e.g., 41]. Also, these overestimates could be attributable to small co-registration errors (<0.5 pixel) between the two lidar surveys, but this would lead to a comparable number of randomly distributed underestimates. In summary, it is unlikely that these errors would alter the sensitivity analysis to pulse density.
Although a quantitative evaluation of lidar data acquisition cost was not a central objective in this study, it is nonetheless an important factor to consider because it can be a primary factor driving choices made about forest and AGB monitoring across a wide range of spatial scales. Because pulse density has a strong influence on the acquisition cost of lidar data, and even though the cost for using lidar with high or low pulse density for AGB in tropical forest might be lower than the cost of a conventional inventory [49], airborne lidar can be cost prohibitive for forest carbon monitoring, in selected logged areas for REDD+ at large spatial extents. Although field-based AGB estimations remain necessary for these purposes, integrating lidar remote sensing into AGB inventory schemes allows recovery of carbon content and spatially-explicit estimates across landscapes, while reducing the total costs and need for extensive field-based sampling.

Conclusions
We evaluated the impacts of airborne lidar pulse density on AGB stocks and AGB change estimation in a selectively logged Amazon tropical forest. First, we confirmed that HMEAN is a stable lidar-derived metric for estimating AGB stock in selective logging. Second, we found that the accuracy of AGB stocks and AGB change estimates decreased as the pulse density decreased, but it remained relatively high except at low pulse densities of 0.8 and 0.2 pulses·m −2 for the DS1 and DS2 scenarios, respectively. Furthermore, AGB stocks estimations at the landscape level were strongly underestimated at pulse densities lower than 0.8 pulses·m −2 in areas with steep slopes, but only in DS1, where the lidar datasets from both 2012 and 2014 were height normalized using the DTMs created from their respective thinned dataset. Therefore, these results demonstrate that high lidar pulse density is not necessary to estimate and map AGB stock and changes in selective logging in tropical forest, especially when there is already an accurate DTM derived from high pulse density lidar. Third, we showed that low pulse density lidar data (~2 pulses·m −2 ) has the ability to map ground topography, allowing accurate estimation of canopy height even over rough terrain and as a baseline for subsequent low density lidar acquisition for AGB change studies. Lower point densities can cover larger areas at reasonable cost and be used to complement satellite remote sensing measurements, e.g., NISAR-National Aeronautics and Space Administration-Indian Space Research Organization Synthetic Aperture Radar (http://nisar.jpl.nasa.gov) and GEDI-Global Ecosystem Dynamics Investigation (http://science.nasa.gov/missions/gedi/), that may have limitations in estimating tree height in areas with complex topography. Finally, although we focus on AGB stocks and AGB change estimation in a selectively logged tropical moist forest in Brazil, our methodology may also be applicable for inventorying and monitoring AGB changes to support REDD+ monitoring efforts in selective logging elsewhere across the tropics.