Estimation of Voxel-Based Above-Ground Biomass Using Airborne LiDAR Data in an Intact Tropical Rain Forest, Brunei

Abstract

The advancement of LiDAR technology has enabled more detailed evaluations of forest structures. The so-called “Volumetric pixel (voxel)” has emerged as a new comprehensive approach. The purpose of this study was to estimate plot-level above-ground biomass (AGB) in different plot sizes of 20 m × 20 m and 30 m × 30 m, and to develop a regression model for AGB prediction. Both point cloud-based (PCB) and voxel-based (VB) metrics were used to maximize the efficiency of low-density LiDAR data within a dense forest. Multiple regression model AGB prediction performance was found to be greatest in the 30 m × 30 m plots, with R², adjusted R², and standard deviation values of 0.92, 0.87, and 35.13 Mg∙ha⁻¹, respectively. Five out of the eight selected independent variables were derived from VB metrics and the other three were derived from PCB metrics. Validation of accuracy yielded RMSE and NRMSE values of 27.8 Mg∙ha⁻¹ and 9.2%, respectively, which is a reasonable estimate for this structurally complex intact forest that has shown high NRMSE values in previous studies. This voxel-based approach enables a greater understanding of complex forest structure and is expected to contribute to the advancement of forest carbon quantification techniques.

1. Introduction

In recent years, awareness of forests under climate change has reached unprecedented levels [1]. The United Nations Framework Convention on Climate Change (UNFCCC) has introduced the REDD+ program, which provides incentives for developing countries to reduce carbon emissions caused by deforestation and forest degradation, and thus going beyond the role of conservation, and sustainable management of forests to enhance forest carbon [1,2]. There have been notable advances in remote sensing technologies in recent decades that have sufficiently enabled the monitoring, reporting, and verification (MRV) systems required by REDD+ [1,2]. The Fifth Assessment Report (AR5) of the United Nations Intergovernmental Panel on Climate Change (IPCC) showed that there are technical opportunities to use remote sensing data toward REDD+ initiatives and that the practicality of the application has been verified through field research in high-biomass forests [2,3]. Nevertheless, most studies have focused on estimating above-ground biomass (AGB) in the relatively low-biomass temperate and boreal forests, while few studies have examined high-biomass tropical forests using Light Detection and Ranging (LiDAR) [4,5]. Tropical forests cover 2.5 billion ha globally, which is 48% of the total global forested area [1]. They store approximately 470 billion tons of carbon in biomass and soil and are responsible for one-third of global terrestrial primary productivity [6]. The importance of tropical forests has been extensively discussed in international discourse and previous studies have investigated whether tropical forests have become more vulnerable to climate change than other forest types in different climate zones. [2,6,7,8].

In particular, completely intact primary tropical rain forests have become extremely rare due to indiscriminate logging and deforestation [4,9], an issue that is expected to contribute to the global greenhouse gas production [4]. In light of this, the current study aims to efficiently estimate plot-level above-ground biomass (AGB) in a government-protected intact tropical rain forest of Brunei using remote sensing data.

Optical imagery has been widely used to estimate AGB in previous studies [10]. However, since this method provides only two-dimensional (2-D) spectral information of the upper canopy, it is ineffective for quantifying AGB and forest carbon stock, especially within dense forests [4,10,11]. The three-dimensional (3-D) point cloud data of LiDAR, however, is a key way to represent and understand structural characteristics of forests [12]. A LiDAR point cloud consists of a number of points that penetrate through the vegetation to the ground, namely from the top canopy layer to the lower canopy layer. It thus provides a more detailed evaluation of the dense and complex forest structure [4,13].

A design of optimal plot size is important in forest inventory because variation in this parameter can reduce or inflate the perceived impact of edge effects when estimating forest stand information [14,15]. LiDAR analysis can be fine-tuned by setting the optimal plot size. For example, a tree with its center located in the boundary of Plot A will be logged in field inventory as belonging to Plot A. However, its crown might broadly extend outside of plot A. To reduce these edge issues, many studies have proposed the use of larger plots [5,16,17], which have smaller perimeter-to-area ratios and are thus less impacted by edge effects. Larger plots, however, require more complicated field measurement [5,18]. Therefore, this study aims to determine the optimal plot size for the study area by analyzing field data collected in both 20 m × 20 m and 30 × 30 m plots.

The use of LiDAR data can be divided into two methods: point cloud-based (PCB) analysis and canopy height model (CHM)-based analysis. Most studies have been performed using CHM, which converts raw LiDAR point cloud data into regular grid format in order to extract individual trees and crown shape, and to estimate biomass [19]. The CHM-based approach has the advantages of relative simplicity and of applicability to various image segmentation methods, such as local minimum and maximum filtering. However, details contained in the original dataset are lost in the process of resampling from 3-D point LiDAR data into the 2-D grid [20,21]. The method therefore does not take advantage of LiDAR’s primary virtue of 3-D point cloud data [22]. PCB approaches that make use of the entire LiDAR point cloud, such as LiDAR height and intensity distributional metrics, have been used in a number of studies. Kwak et al. [23] estimated AGB using height and intensity percentiles, skewness, and kurtosis, etc., and Yoon et al. [24] used these parameters to verify the optimal plot size in the same study area. Additionally, a ‘height bins’ approach has been implemented recently in a number of cases. This approach separates vertical space at various height intervals (e.g., 5 m intervals) and then examines the occurrence and total number of points within each height bin [25,26,27].

The so-called “volumetric pixel (volume pixel or voxel)” emerged as a comprehensive approach to illustrate both vertical and horizontal slices of LiDAR point clouds [26,27,28]. A voxel is defined as an array of 3-D cubes, in the form of either a regular or an irregular grid, separated by user-defined horizontal and vertical resolutions, which are associated with the concept of the horizontal slice and vertical slice (or height bin), respectively [28,29]. In this study, voxel and sub-voxel (SV) arrays were created with user-defined horizontal and vertical resolutions. LiDAR point cloud data were processed to fall within each voxel space and these were used to build voxel-based (VB) metrics. VB metrics enable stratum-specific analyses to be conducted because, unlike PCB metrics, they produce variables using points derived from specific height layers. These height-specific metrics are expected to be useful in tropical rain forests where multi-layer structured vegetation is prevalent. Thus, this study aims to use both PCB and VB metrics to maximize the efficiency of low-density airborne LiDAR data with a 3-D point cloud format and to resolve the optimal plot size for the metrics by comparing the accuracy. Such new voxel-based improvements to LiDAR analysis are expected to produce enhanced methodologies and to contribute to the advancement of plot-level forest AGB and forest carbon quantification techniques.

2. Materials and Methods

2.1. Study Area

The study area was located in the Kuala Belalong Field Studies Centre (KBFSC; 04°32′50.1″ N, 115°09′29.6″ E) in the Ulu Temburong National Park, Brunei Darussalam. The KBFSC is 25 ha of intact tropical rain forest, which consists of lowland mixed dipterocarp forest (MDF), at an elevation range of 49 m to 307 m above sea level. Meteorological data for the study area were obtained from the Semabat area of Semabat Agricultural Station, which is 9 km west-northwest of the KBFSC. These data were analyzed by collecting climate data from 2007 to 2010 and rainfall data from 2006 to 2010 [9,13,24]. The daily maximum temperature ranged from 30.58 °C to 35.08 °C and the daily minimum temperature ranged from 25.6 °C to 26.2 °C. The average annual rainfall during the study period was 4582 mm, with 240 rain days per year. The spatial pattern of the rainfall was heavily concentrated in a small area over a short period of time [30]. While small and young understory trees had a monopodial crown shape, mature trees had a sympodial crown shape in which the crown width is broader than its height [10,31].

Researchers at the Center for Tropical Forest Science (CTFS) and University of Brunei Darussalam (UBD) established 25 ha of forest dynamics plots for the field study of carbon stock estimation and carbon sequestration, and have monitored the area since 2011 (Figure 1a). Field inventory data for six 60 m × 60 m plots, located in the northeastern part of the forest dynamics plot, were provided by UBD at a horizontal collection interval of 10 m. We subdivided the plots further into 20 m × 20 m and 30 m × 30 m for the estimation of AGB (Figure 1b), for a total of fifty-four 20 m × 20 m plots and twenty-four 30 m × 30 m plots (Figure 1b).

Field inventory data were obtained in 2011 from the CTFS. Data on diameter at breast height (DBH) over 1 cm and the relative position of individual trees in each plot were collected. DBH of individual trees was less than 10 cm in 91% of the study area. Above-ground biomass (AGB) was calculated from trees whose DBH was ≥1 according to the allometric equation, and plot-level AGB of the study area was estimated by summing the number of individual trees located in each plot. The allometric equation was developed in lowland dipterocarp forests by Basuki et al. ([32], Equation (1)).

$$\ln \left(\mathrm{TAGB}\right)=-1.232+2.178\times \ln \left(\mathrm{DBH}\right)$$

(1)

where TAGB is the total above-ground biomass measured in kg/tree and DBH is measured in cm.

The number of individual trees with DBH < 10 cm occupied 91% of the total area and yet their AGB contributed only 8.47% of the total biomass in the 2.16 ha area. In contrast, trees with DBH > 10 cm occupied only 9% of the study site and yet their AGB corresponded to 91.53% of the total biomass in the area (

Table 1. Descriptive statistics of the field inventory data.

Estimator	Statistics	20 m × 20 m Plot	30 m × 30 m Plot
DBH (cm)	avg.	4.6
	max.	161.5
	min.	1.0
	standard deviaion	7.8
Above-ground Biomass (DBH ≥ 1 cm) (Mg∙ha−1)	avg.	313.75	302.71
	max.	904.63	585.88
	min.	77.39	154.12
	standard deviation	150.81	98.51

). To reduce the influence of edge effects, plot-level AGB for 20 m × 20 m and 30 m × 30 m plots was calculated. The edge effect determines that as the size of plot becomes smaller, big trees have greater influence and result in a higher standard deviation (SD). Supporting this theoretical knowledge, the average AGB ± SD of the 20 m and 30 m plots were 313.75 ± 150.81 Mg∙ha⁻¹ and 302.71 ± 98.51 Mg∙ha⁻¹, respectively. Additionally, 20 m plots were characterized by a significantly greater variation in each plot compared with the 30 m plots.

The number of individual trees with DBH < 10 cm occupied 91% of the total area and yet their AGB contributed only 8.47% of the total biomass in the 2.16 ha area. In contrast, trees with DBH > 10 cm occupied only 9% of the study site and yet their AGB corresponded to 91.53% of the total biomass in the area (Table 1). To reduce the influence of edge effects, plot-level AGB for 20 m × 20 m and 30 m × 30 m plots was calculated. The edge effect determines that as the size of plot becomes smaller, big trees have greater influence and result in a higher standard deviation (SD). Supporting this theoretical knowledge, the average AGB ± SD of the 20 m and 30 m plots were 313.75 ± 150.81 Mg∙ha⁻¹ and 302.71 ± 98.51 Mg∙ha⁻¹, respectively. Additionally, 20 m plots were characterized by a significantly greater variation in each plot compared with the 30 m plots.

2.2. Airborne LiDAR Data

Airborne LiDAR data were acquired at an average of 1400 m above ground level (AGL) in April 2009 with a sampling density of 1.13 pt∙m⁻² (Figure 2). LiDAR acquisition was conducted with a laser footprint on the ground of 34 cm and a scan angle of ± 20° at the rate of 35 Hz. The LiDAR system recorded four returns from a single laser pulse; 1st return, intermediate return 1 (2nd), intermediate return 2 (3rd), and last return (4th). Each return contained X, Y, Z, and intensity values.

For the normalization of the raw LiDAR point cloud data, the normalized height, referring to the Z value of terrain, was calculated (Figure 2). First, a digital terrain model (DTM) was generated from last returns using a triangulated irregular network (TIN) interpolation method at a spatial resolution of 1 m [19]. Second, the normalized height values were determined as the height difference between each raw LiDAR point and its corresponding DTM value [10,33,34]. Last, abnormal values such as NULL values and negative height values were removed.

2.3. LiDAR Metrics

Full usage of LiDAR point cloud data allows a more sophisticated and detailed approach, unlike CHM based analysis. Typical examples of full usage are LiDAR height and intensity distributional metrics and statistics such as percentile, maximum height, minimum height, and intensity. Additionally, the volumetric pixel (voxel), implemented in this study, is a comprehensive approach that illustrates both vertical and horizontal slices of LiDAR point cloud data. This study used both point cloud-based (PCB) and voxel-based (VB) metrics to gain a better understanding of the 3-D forest structure in the multi-layered forest of the KBFSC study area (

Table 2. Description of LiDAR point cloud based (PCB) metrics and voxel-based (VB) metrics.

LiDAR Metrics		Description
Point cloud based (PCB) metrics
Height	HP10, HP20, …, HP100	Percentile height of total returns at 10 percentile intervals
	Hmax, Hmin, Havg, Hmed	Maximum, minimum, average, and median height of total returns
	Hrange, Hstd	Difference of the maximum and minimum height of total returns, Standard deviation of total returns
Intensity	Imax, Imin, Iavg, Imed	Maximum, minimum, average, and median intensity of total returns
	Irange, Istd	Difference of the maximum and minimum intensity of total returns, Standard deviation of total returns
Ratio	PT	LiDAR points of total returns
	IT	Summed intensity of total returns
Voxel-based (VB) metrics (i = 1, 2, …, 8)
SVi	P_SVi/P_SViF/P_SViFS	LiDAR points of total returns/1st returns/the sum of 1st and 2nd returns within each sub-voxel (SVi) a
	FR_SVi/FR_SViF/FR_SViFS	Frequency ratio of the SVi returns/1st returns/the sum of 1st and 2nd returns over all returns
	I_SVimed/I_SVimed_F/I_SVimed_FS	Median intensity of total returns/1st returns/the sum of 1st and 2nd returns within each SVi
Di	P_Di/P_DiF/P_DiFS	LiDAR points of total returns/1st returns/the sum of 1st and 2nd returns above their corresponding sub-voxel (Di) b
	FR_Di/FR_DiF/FR_Di FS	Frequency ratio of the Di returns/1st returns/the sum of 1st and 2nd returns over all returns c
	I_Dimed/I_Dimed_F/I_Dimed_FS	Median intensity of total returns/1streturns/the sum of 1st and 2nd returns above corresponding sub-voxel (Di) b
SVM	H_SVMmed/H_SVMmed_F//H_SVMmed _FS	Median height of total returns/1st returns/the sum of 1st and 2nd returns within SVM
	H_SVM_Dmed/H_SVM_Dmed_F/H_SVM_Dmed _FS	Median height of total returns/1st returns/the sum of 1st and 2nd returns above SVM

^a ‘P_SV1’ was removed to avoid duplication; ^b Both ‘P_D1’ and ‘P_D8’ were removed to avoid duplication; ^c Both ‘FR_D1’ and ‘FR_D8’ were removed to avoid duplication.

). In total, 118 LiDAR metrics were derived from the 20 m × 20 m plots and 30 m × 30 m plots.

2.4. Point Cloud Based (PCB) LiDAR Metrics

LiDAR height and intensity distributional parameters were built using the normalized LiDAR point cloud height and intensity values of each point. These were treated as independent variables, which are the explanatory variables for estimating plot-level AGB. A total of 29 variables were classified into three data sets for interpreting the vertical height profile of each plot: (1) height data set; (2) intensity data set; and (3) ratio metrics data set [23,33,35]. These variables were recorded at 10 percentile intervals from 0 to 100 percentile height (H_P10, H_P20, …, H_P100), for the purpose of defining the relationship between AGB and height [36]. In addition, maximum, minimum, average, and median statistics of height and intensity values were also examined. A more detailed description of the variables is presented in Table 2.

2.5. Voxel-Based (VB) Metrics

2.5.1. Extraction of LiDAR Point Cloud Data Based on the Voxel Approach

Unless using CHM analysis, full usage of LiDAR points is difficult to achieve due to the massive quantity of points. To solve this issue, we developed an algorithm in this study to conduct a time-efficient voxel-based approach analysis in RStudio (R version 3.2.3). Voxels in the form of a regular 3-D grid (similar to a 3-D cube) were created at a pre-defined horizontal and vertical resolution. In accordance with a pre-defined plot size (60 m × 60 m) for field measurements, the X, Y values were set as 60 m × 60 m and the Z value was set to the maximum LiDAR point height, 80 m (Figure 3a). X, Y values were then subdivided into both 20 m × 10 m and 30 m × 30 m intervals (Figure 3b,d) and the Z value was subdivided into user-defined 10 m vertical intervals (Figure 3c,e).

Finally, the voxel-based LiDAR metrics were computed based on the LiDAR points within each sub-voxel (SV). From this point forward, each SV that was divided into 10-m vertical intervals is referred to as SV1 (0–10 m), …, SV7 (60–70 m), SV8 (70–80 m) (Figure 3c,e).

2.5.2. Concepts of Voxel-Based (VB) Metrics

The voxel-based LiDAR metrics employed in this study can be divided into three basic variable concepts as follows; Variable Sub-Voxel i (SVi), Variable Density i (Di), Variable Sub-Voxel Maximum (SVM) (Figure 4a–c). Sheridan et al. [27] proposed the concepts of Variable SVi and Variable Di, and a new variable approach of Variable SVM was added in this study. Using these three concepts for building VB metrics, a number of points, frequency ratios, median intensities, and median height values were calculated (Table 2). To find out how dense the corresponding voxel was compared to all returns, the frequency ratio variable was calculated as the fraction of total LiDAR points versus the points that existed within or above SVi. Variable SVi accounted for the variables derived from the points in their corresponding SV (Figure 4a). For example, P_SV5 represented the total LiDAR returns within SV5, ranging from 40 m to 50 m. In contrast, Variable Di described the cumulative points located above the corresponding SV (Figure 4b). For example, I_D5_med represented the median intensity of accumulative LiDAR points above SV5, ranging from 40 m to 80 m. Variable SVM was the location of the greatest point density, where the highest frequency of LiDAR returns occurred within the entire range of SV1 to SV8 (Figure 4c,d). After locating SVM, the points within (SVi) or above (Di) SVM were extracted and used for building further variables. For example, H_SVM_med represented the median height of SV, which exhibited the highest point density throughout SV1 to SV8.

Intermediate return is expected to have high applicability in dealing with multi-layered forests, such as tropical rain forests, because these returns have the ability to reflect mid-canopy forest structure [12]. However, only a few studies have attempted to incorporate intermediate returns (2nd, 3rd returns) on a deeper level [37]. In this study, 2nd returns were added to the database of the 1st and last returns in hopes of increasing the accuracy of the analysis and obtaining a better understanding of the multi-layered forest. For various purposes, all variables were extracted using either total LiDAR returns, first returns only (e.g., for P_Di_F), and the sum of 1st and 2nd returns (e.g., for P_Di_FS).

2.6. Statistical Analysis

To predict plot-level AGB using LiDAR data, a multiple linear regression analysis (Equation (2); [23]) was performed with a stepwise selection method at a significance level of 0.05 with ten or less predictors as a result.

$$\mathrm{y}=\mathsf{\alpha }+{\mathsf{\beta }}_1{x}_1+{\mathsf{\beta }}_2{\mathrm{x}}_2+\cdots +{\mathsf{\beta }}_n{\mathrm{x}}_{\mathrm{n}}+\mathsf{\epsilon }$$

(2)

where AGB is the plot-level AGB in Mg∙ha⁻¹; x₁, x₂, …, x_n are the selected LiDAR variables as predictors; α is the intercept of the regression; β₁, β₂, …, β_n are the slope coefficients of the regression; and ε is the residual.

Multiple regression analysis was applied separately to the 20 m × 20 m plot and the 30 m × 30 m plot. To determine the optimal model for each by plot size, the statistics coefficient of determination (R²), adjusted coefficient of determination (Adj. R²) and standard error (SE) were calculated for each model using RStudio (R version 3.2.3, RStudio, Inc., Boston, MA, USA) and compared. Furthermore, to avoid multicollinearity among independent variables, a variance inflation factor (VIF) of less than 10 was required for selection as an independent variable [23,27]. To avoid inaccuracies from an excessive number of the predictors, the total number of selected independent variables was limited to 10. Predicted models from multiple regression analyses were validated using root-mean-square error (RMSE) and normalized root-mean-square error in percentage (NRMSE%). RMSE indicates the absolute value of error, whereas NRMSE% is the ratio of RMSE versus the average value of the dependent variables ($\overline{y}$) and thus represents the relative value of the error with respect to the average plot-level AGB (Equations (3) and (4)).

$$\mathrm{RMSE}\left(\mathrm{Mg}·\text{h}{\text{a}}^{-1}\right)=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^n{\left(y_{observed}-y_{predicted}\right)}^2}{n}}$$

(3)

$$\mathrm{NRMSE}(\%)=\frac{RMSE}{\overline{y}}\times 100$$

(4)

where y_observed is the plot-level field AGB in Mg∙ha⁻¹, y_predicted is the LiDAR-estimated plot-level AGB in Mg∙ha⁻¹, and $\overline{y}$ is the average value of plot-level field AGB in Mg∙ha⁻¹.

3. Results and Discussion

3.1. Development of Models for AGB Estimation

Multiple linear regression analysis for each plot size was performed to predict plot-level AGB at a significance level of 0.05. In total, 118 LiDAR metrics combining point cloud based (PCB) metrics and voxel-based (VB) metrics were used for developing and assessing the accuracy of the regression models The R², Adj. R², and SE values of the 20 m × 20 m plot model (n = 54) were 0.72, 0.69, and 107.56 Mg∙ha⁻¹, respectively (

Table 3. Results of multiple regression analysis on 20 m × 20 m and 30 m × 30 m plot size.

Plot Size	Variables		Estimate	Standard Error	t Value	Pr (>|t|)	VIF
20 m × 20 m (0.04 ha)	Intercept		710.10	187.41	3.79	0.00
	VB	P_SV1	−1.74	0.56	−3.12	0.00	4.89
		I_SV2med_FS	−60.19	27.15	−2.22	0.03	1.86
		H_SVMmed	9.77	1.55	6.31	0.00	1.36
	PCB	HP20	−25.98	5.98	−4.34	0.00	7.31
30 m × 30 m (0.09 ha)	Intercept		757.02	186.19	4.07	0.00	−
	VB	FR_SV1	−3048.11	806.61	−3.78	0.00	4.41
		I_SV3med	110.27	27.34	4.03	0.00	4.36
		I_SV3med_FS	52.26	13.25	3.95	0.00	1.54
		I_D8med	94.03	19.18	4.90	0.00	2.36
		H_SVMmed	3.06	0.79	3.87	0.00	1.70
	PCB	HP10	−12.34	4.32	−2.86	0.01	6.40
		Irange	−149.38	41.78	−3.58	0.00	6.66
		Istd	232.04	87.58	2.65	0.02	2.37

). Four independent variables were selected, including one variable from PCB metrics and three variables from VB metrics. The R², Adj. R², and SE values of the 30 m × 30 m plots model (n = 24) were 0.92, 0.87, and 35.13 Mg∙ha⁻¹, respectively. The model defined for both scenarios selected more independent variables from VB metrics than from PCB metrics. It is postulated that VB metrics held more power than PCB metrics in explaining AGB in stratum-specific analysis because its variables were produced from points at various height layers, whereas the PCB metrics did not implement height stratification (Table 3).

The coefficients of all selected independent variables were interpreted. Since the coefficients for variables P_SV1, H_P10, H_P20, and FR_SV1 were negative, it is understood that AGB tends to increase with the number of points in the understory layer of 0 to 20 m. This indicates that in plots where the emergent and canopy layer are dominant, the growth of understory vegetation is hindered because less than 2%–15% of sunlight can reach the ground. On the other hand, plots with less understory vegetation can also be interpreted as having a prevalence of large trees. In this respect, it is determined that the variables representing understory layer have a negative relationship with plot AGB.

Variable “H_SVM_med” was notable as it was selected in both 20 m × 20 m and 30 m × 30 m plot models. This variable is comprised of median height values from SVM, which had the highest point density among SV1 to SV8 (0–80 m). This indicates the greater the elevation of the densest forest layer, the higher its explanatory power of AGB.

Intensity values have often been used for species classification and are one of the key factors used in understanding forest biophysical parameters [23,38,39,40,41,42]. Six intensity related variables were also adopted in this study. VB intensity variables I_SV2_{med_FS}, I_SV3_med, I_SV3_{med_FS}, and I_D8_med and PCB intensity variables I_Range and I_Std were selected as the predictors of the regression models. The mean intensity value of points was calculated at 2-m intervals between the heights of 0 m and 80 m, and it was observed that intensity value tended to increase with canopy height (Figure 5a) with a high correlation between height group and mean intensity (Figure 5b). This finding indicates that intensity-related variables and height values of LiDAR points can be representative descriptors of plot AGB.

3.2. Statistical Analysis for Model Validation

AGB values predicted using the regression model were assessed for accuracy using field-measured AGB (observed AGB) data (

Table 4. Summary of field-measured (observed) above-ground biomass (AGB) and LiDAR-predicted (predicted) AGB (Mg∙ha⁻¹) by plot size.

Plot size	Estimator	Average	Standard Error	Min	Max	n
20 m × 20 m (0.04 ha)	Observed AGB	313.75	150.81	77.39	904.63	54
	Predicted AGB	290.85	98.16	135.23	565.9	54
30 m × 30 m (0.09 ha)	Observed AGB	302.71	98.51	154.12	585.88	24
	Predicted AGB	302.71	94.33	147.59	591.17	24

; Figure 6). To validate the model accuracy, values of RMSE (Mg∙ha⁻¹) and NRMSE% were examined. Model validation for the 20 × 20 m plots yielded RMSE and NRMSE values of 77.57 Mg∙ha⁻¹ and 24.72%, respectively. Observed AGB was 313.75 ± 150.81 Mg∙ha⁻¹, while the LiDAR model under-estimated AGB at 290.85 ± 98.16 Mg∙ha⁻¹. Model validation for the 30 × 30 m plot yielded RMSE and NRMSE% values of 27.77 Mg∙ha⁻¹ and 9.17%, respectively. Observed and predicted AGB in the 30-m plots showed a high level of agreement, as demonstrated by the same average AGB value of 302.71 Mg∙ha⁻¹. The discrepancy in accuracy between the 20-m and 30-m plot AGB predictions may have been caused by the differing sample size available (20-m plots n = 54 and 30-m plots n = 24). However, since the field inventory data indicated that the SE value of the 20-m plots (150.81 Mg∙ha⁻¹) was greater than the 30-m plots (94.33 Mg∙ha⁻¹), it was determined that the smaller plot size is more highly influenced by edge effects. Accordingly, the 30 m × 30 m plots were determined to be more reliable for research purposes than 20 m × 20 m plots in this study area.

Very few studies have focused on high-biomass forests with average AGB over 200 Mg∙ha⁻¹ [3,4,43,44,45]. High-biomass forests such as those found in the tropics could have a significant bias in AGB estimation owing to greater absolute error caused by their more complicated structure relative to other forest types (biomes). To verify how well our high-biomass AGB prediction models performed, we compared our results with those of previous airborne LiDAR studies that were conducted in high-biomass forests (>200 Mg∙ha⁻¹) on 30 m × 30 m plots (

Table 5. Comparison of preceding plot-level AGB estimation studies.

Study	Forest Biomes	Nation	Method	Average AGB a (Mg∙ha−1)	LiDAR Point Density (pts∙m−2)	R2/Adj. R2	RMSE (Mg∙ha−1)	NRMSE%
This Study (2016)	Intact tropical rain forest	Brunei	VB + PCB	302.7 b,c (DBH ≥ 1)	1.1	0.92/0.87	27.8	9.2
Phua et al. (2016)	Tropical rain forest	Malaysia	PCB	358.6 (DBH ≥ 10)	25.0	0.67/- c	80.0	22.3
Ioki et al. (2014)	Tropical rain forest	Malaysia	PCM + CHMB	235.5 (DBH ≥ 10)	14.9	- c/0.81	61.26	26.0
Mauya et al. (2015)	Tropical rain forest	Tanzania	PCB	359.4 (DBH ≥ 5)	10.6	- d/0.74	- c	29.2
Hansen et al. (2015)	Tropical submontane rain forest	Tanzania	PCB	461.9 (DBH ≥ 10)	10.0	0.71/- d	158.0	34.4
Lu et al. (2012)	Temperate rain forest	USA	PCM + CHMB	239.1 (DBH ≥ 5)	2.0	0.76/- d	82.1	34.3

VB, Voxel-based; PCB, Point cloud based; CHMB, Canopy Height Model based; ^a Field measured AGB; ^b the reference to 30 m × 30 m plots; ^c 302.5 Mg∙ha⁻¹ (DBH ≥ 5 cm) and 302.3 Mg∙ha⁻¹ (DBH ≥ 10 cm); ^d Corresponding information is not given.

). For quantitative evaluation, average AGB (Mg∙ha⁻¹), LiDAR point density (pts∙m⁻²), RMSE (Mg∙ha⁻¹), and NRMSE% were reviewed. Overall LiDAR point density of the comparison studies was 2–25 times higher than that used in this study and their average field AGB was 235.5 to 461.9 Mg∙ha⁻¹.

Phua et al. [4] performed a PCB approach on a logged tropical rain forest in Malaysia, close to our study area in Borneo. Their study area was also lowland MDF, akin to our study. They created height percentile variables as predictors and derived R² and NRMSE values of 0.67% and 22.3%, respectively. Thus the current study yielded higher confidence performance despite a 25-fold lesser LiDAR point density. Ioki et al. [43] also estimated plot-level AGB in Borneo, Malaysia, in the transition zone between lowland and mountain forests through the combined use of mean CHM value and height percentiles derived from the point cloud basis. Their average AGB was 235.5 Mg∙ha⁻¹, which is a relatively lower level than that found in this study. The values of Adj. R² and NRMSE were 0.81% and 26.0%, respectively. Mauya et al. [3] used the PCB approach for their analysis with optimal plot size in the tropical rain forest of Tanzania. This study concluded with an NRMSE value of 29.2% at the plot size of 2000 m², with a high LiDAR point density of 10.6 pts∙m⁻². Additionally, the study area of Hansen et al. [44] had a considerably high AGB of 461.9 Mg∙ha⁻¹, determined using a PCB approach to describe height and density of the forest canopy for R² and NRMSE values of 0.71% and 34.4%, respectively. Lastly, Lu et al. [45] predicted an average field AGB of 239.1 Mg∙ha⁻¹ in a temperate rain forest within the USA. The R² value was 0.76, but NRMSE was relatively high at 34.4%, and LiDAR point density was a low 2 pts∙m⁻².

The current study is notable in that it efficiently analyzed the complicated intact tropical rain forest with a high accuracy using the new approach. The lowland MDF of this study area contains a morphologically different crown shape according to the age of each tree. Small and young understory trees have a monopodial crown shape, while mature trees have a sympodial crown shape with a broader crown width than crown height [10,31]. Furthermore, the characteristics of the tropical rain forest, especially ‘intact’ rainforest are more complex and multi-layered than other forests, with different structural objects and a high stand density. The current study implemented a very low LiDAR density (1 pt∙m⁻²) compared to previous research, which ranged from 2 to 25 pts∙m⁻². Nevertheless, the AGB prediction models developed in this study showed more accurate results than previous studies.

4. Conclusions

The present study aimed to estimate plot-level AGB and to develop a regression model through a voxel-based approach using airborne LiDAR data from an intact tropical rain forest. The comprehensive volumetric pixel (voxel) approach illustrated both vertical and horizontal slices of LiDAR point cloud data in the form of a 3-D cube, and was expected to sufficiently represent the structural complexity of the forest.

By using both point cloud-based (PCB) metrics and voxel-based (VB) metrics, we were able to promote the advantages of LiDAR data, which are able to describe 3-D forest structures with 3-D point cloud data. After calculating plot-level field AGB from DBH data of individual trees at two different plot sizes of 20 m × 20 m and 30 m × 30 m, multiple linear regression analyses were performed to derive the best models for estimating AGB at each plot size. The results of regression analysis indicated that the 30 m × 30 m plot model was more accurate in predicting AGB than the 20 m × 20 m plot model, with R², adj. R², and SE values of 0.92, 0.87, and 35.13 Mg∙ha⁻¹, respectively. Both 20 m × 20 m and 30 m × 30 m models selected a greater number of independent variables from VB metrics than from PCB metrics. It is postulated that VB metrics were selected because their stratum-specific AGB quantification of AGB provides greater power in AGB prediction. PCB metrics do not have this advantage because they are not derived from points located in specific height layers. We therefore expect the VB metrics to be more applicable in tropical rain forests where multi-layer structured vegetation is prevalent.

Comparing predicted and observed AGB for model validation, the value of RMSE was 27.8 Mg∙ha⁻¹. To relate the value of RMSE to the average field value of AGB, NRMSE% was calculated as 9.2%. This study provides a high confidence estimate for complicated intact tropical rain forest, relative to those found in previous research (NRMSE% from 22.3% to 34.4%), even though the density of our LiDAR returns was low. This new voxel-based approach is expected to contribute to the advancement of forest carbon quantification techniques. However, this study is limited in that it is constrained to an intact MDF forest in a tropical region within a relatively small number of 24 sample plots. This research can be furthered by applying the same techniques in other forest climatic zones (e.g., the temperate forests of Korea) to examine their effectiveness more broadly.