Aboveground Biomass Changes in Tropical Montane Forest of Northern Borneo Estimated Using Spaceborne and Airborne Digital Elevation Data

: Monitoring anthropogenic disturbances on aboveground biomass (AGB) of tropical montane forests is crucial, but challenging, due to a lack of historical AGB information. We examined the use of spaceborne (Shuttle Radar Topographic Mission Digital Elevation Model (SRTM) digital surface model (DSM)) and airborne (Light Detection and Ranging (LiDAR)) digital elevation data to estimate tropical montane forest AGB changes in northern Borneo between 2000 and 2012. LiDAR canopy height model (CHM) mean values were used to calibrate SRTM CHM in di ﬀ erent pixel resolutions (1, 5, 10, and 30 m). Regression analyses between ﬁeld AGB of 2012 and LiDAR CHM means at di ﬀ erent resolutions identiﬁed the LiDAR CHM mean at 1 m resolution as the best model (modeling e ﬃ ciency = 0.798; relative root mean square error = 25.81%). Using the multitemporal AGB maps, the overall mean AGB decrease was estimated at 390.50 Mg / ha, but AGB removal up to 673.30 Mg / ha was estimated in the managed forests due to timber extraction. Over the 12 years, the AGB accumulated at a rate of 10.44 Mg / ha / yr, which was attributed to natural regeneration. The annual rate in the village area was 8.31 Mg / ha / yr, which was almost 20% lower than in the managed forests (10.21 Mg / ha / yr). This study identiﬁed forestry land use, especially commercial logging, as the main driver for the AGB changes in the montane forest. As SRTM DSM data are freely available, this approach can be used to estimate baseline historical AGB information for monitoring forest AGB changes in other tropical regions. (CHM) for developing an SRTM-to-LiDAR canopy height calibration model and to determine the best spatial resolution for developing an AGB estimation model using ﬁeld and LiDAR data that can be applied to the SRTM canopy height values. Using the developed model, we examined the AGB changes between 2000 and 2012 in the study area.


Introduction
Emissions from land-use changes, such as deforestation, logging, and intensive cultivation of cropland, are the second-largest anthropogenic emissions source after fossil fuel emissions [1]. Annually, 15-25% of global greenhouse gas emissions are produced by the loss of tropical rainforests due to human activities [2]. Although deforestation contributes to carbon emissions, forest degradation is the result of human-induced activities that lead to a long-term reduction in forest carbon stocks. In Borneo, most lowland primary forest has disappeared as a result of deforestation and forest degradation over the past decades [3,4]. The remaining upland rainforests are severely threatened by increasing anthropogenic activities, particularly in the mountains of the Malaysian Borneo adjacent to the Indonesian border

Study Area
The study area, Ulu Padas, is located at the northern part of Borneo Island near the border between Sabah and Sarawak of Malaysia and Kalimantan of Indonesia (4 • 23 -27 N, 115 • 42 -47 E; Figure 1). The topography of Ulu Padas varies from undulating to hilly. Muruk Miau (2049 m) is the highest peak. This area receives rainfall between 2000 and 3500 mm annually (Ministry of Tourism, Sabah, unpublished). Dominant forest types in the study area include hill dipterocarp, stunted montane mossy, swamp, and oak chestnut forests. Many endemic plant species are found in Virgin Jungle Reserves located in the northern part of the Ulu Padas [44]. developing an SRTM-to-LiDAR canopy height calibration model and to determine the best spatial resolution for developing an AGB estimation model using field and LiDAR data that can be applied to the SRTM canopy height values. Using the developed model, we examined the AGB changes between 2000 and 2012 in the study area.

Study Area
The study area, Ulu Padas, is located at the northern part of Borneo Island near the border between Sabah and Sarawak of Malaysia and Kalimantan of Indonesia (4°23′-27′ N, 115°42′-47′ E; Figure 1). The topography of Ulu Padas varies from undulating to hilly. Muruk Miau (2049 m) is the highest peak. This area receives rainfall between 2000 and 3500 mm annually (Ministry of Tourism, Sabah, unpublished). Dominant forest types in the study area include hill dipterocarp, stunted montane mossy, swamp, and oak chestnut forests. Many endemic plant species are found in Virgin Jungle Reserves located in the northern part of the Ulu Padas [44].

Field Measurement
Field data were collected at both sites from October 2011 to October 2012. Within the two sites, a total of 64 plots were established with 49 and 15 plots at sites 1 and 2, respectively, with different plot sizes depending on DBH: (1) 59 plots of 30 × 30 m, (2) four plots of 20 × 20 m, and (3) one plot of 50 × 50 m. Coordinates of the plot centers were determined using the differential Global Navigation Satellite System (DGNSS) from Ashtech ProMark 100, Spectra Precision, and Triumph-1, Javad GNSS. Then, positions of the plot corners were determined using compass and laser rangefinder from the plot center. In total, 4611 trees were measured for DBH and tree height. Trees with a minimum DBH of 10 cm (5 cm for 20 × 20 m plots) were measured for tree height and DBH. Plant specimens were collected for species identification at the herbarium of Forest Research Center (FRC), Sandakan, Sabah. Tree species were identified to the genus or species level to obtain wood density values for AGB estimation. More than 400 species were found in the study area, including Macaranga spp. from the Euphorbiaceae family, Litsea spp. from the Lauraceae family, and Syzygium spp. from the Myrtaceae family. AGB was calculated using the allometric equation of Chave et al. (2014) as shown in Equation (1). Allometry has been widely used across the Afro-tropical region for both old-growth and secondary forests. AGB = 0.0673 × (ρD 2 H) 0.976 (1) in which, ρ = wood specific gravity (g per cm 3 ), D =DBH (cm)H = tree height (m) The calculated AGB of all trees in a plot were then summed and converted to Mg per hectare (Mg/ha).

SRTM Data
The National Aeronautics and Space Administration (NASA) successfully implemented SRTM to collect a set of interferometric C-band SAR data in February 2000. The DEM produced by the mission had a ±16 m absolute vertical accuracy at a 90% confidence level. Since interferometric SAR signals transmitted from SRTM have a weak penetration rate, the DEM elevation values of vegetation areas likely represent the canopy level or a digital surface model (DSM). We, therefore, considered the DEM as DSM in this study. The SRTM DSM (30 m resolution) data were downloaded from USGS's Earth Explorer website (http://earthexplorer.usgs.gov/). The downloaded SRTM DSM has a geo-reference of the World Geodetic system 1984 (WGS84) and EGM96 geoid height for horizontal and vertical datum, respectively.

LiDAR Data
LiDAR data acquisition was carried out from 5th to 8th of October in 2012 using Riegl LMS-Q560 laser measurement system (Riegl LMS GmbH, Horn, Austria) mounted on a helicopter (Bell 206B3) that flew approximately 400 m above ground level at a speed of 50 to 60 m/s with a 45 • field of view (±22.5 • scanning angle). The laser measurement system was operated using a near-infrared laser with a pulse repetition frequency of 240 kHz and less than 0.5 mrad beam divergence. The resulting point clouds had an average density of 14.9 shots/m 2 at site 1 and 16.2 shots/m 2 at site 2. Overall, 15 lines were scanned with 30% to 50% side overlap in site 1 and 20 lines with about 50% side overlap in site 2. Processed point clouds had a vertical accuracy of 0.15 and 0.25 m for clear ground and vegetated area, respectively. Ground return points were first classified while the remaining of the returns as vegetation points using TerraScan software. To compare the LiDAR data to the SRTM data, the elevation values of these datasets need to be assigned to the same horizontal and vertical datum. Hence, the vertical datum of LiDAR data was converted from WGS84 ellipsoid to Earth Gravitational Model (EGM96) geoid to match the vertical datum of SRTM DSM [45,46] using the ArcGIS Data Management Tools.

Remotely Sensed Digital Elevation Data Processing
For LiDAR data, the point clouds were classified as ground and vegetation points using the TerraScan software. The ground points were interpolated to a grid with 1 m spacing as a digital terrain model (LiDAR DTM), while a digital surface model (DSM) was generated using the return with the maximum altitude, also at a 1 m pixel resolution (LiDAR DSM). The LiDAR DSM and DTM models were resampled to 5, 10, and 30 m spatial resolutions using the nearest neighbor resampling method. A LiDAR canopy height model (LiDAR CHM) was calculated by subtracting the LiDAR DSM from the LiDAR DTM for each resolution.
We analyzed the Pearson correlation (r, Equation (2)) and compared the z or elevation value differences by calculating the root mean square error (RMSE, Equation (3)) of z, the bias of z (biasz, Equation (4)), and the relative bias of z (rbiasz, Equation (5)) of bare ground surfaces, such as bare land and roads, between LiDAR and SRTM DSM data. More than 200 points of bare land and roads were interpreted and digitized on sites 1 and 2. Pearson's r showed that the SRTM DEM and LiDAR DTM had an almost one-to-one relationship ( Figure 2). There were some expected differences in the z or elevation values between the SRTM DSM and LiDAR DTM as these were products of different remote sensing technologies on different platforms (Table 1).
where n = sample size; x = elevation value of the SRTM DSM; y = elevation value of the LiDAR DTM; y i = mean elevation value of the LiDAR DTM for plot I;ŷ i = mean elevation value of the SRTM DSM for plot I; y = mean elevation value of the LiDAR DTM.

Remotely Sensed Digital Elevation Data Processing
For LiDAR data, the point clouds were classified as ground and vegetation points using the TerraScan software. The ground points were interpolated to a grid with 1 m spacing as a digital terrain model (LiDAR DTM), while a digital surface model (DSM) was generated using the return with the maximum altitude, also at a 1 m pixel resolution (LiDAR DSM). The LiDAR DSM and DTM models were resampled to 5, 10, and 30 m spatial resolutions using the nearest neighbor resampling method. A LiDAR canopy height model (LiDAR CHM) was calculated by subtracting the LiDAR DSM from the LiDAR DTM for each resolution.
We analyzed the Pearson correlation (r, Equation (2)) and compared the z or elevation value differences by calculating the root mean square error (RMSE, Equation (3)) of z, the bias of z (biasz, Equation (4)), and the relative bias of z (rbiasz, Equation (5)) of bare ground surfaces, such as bare land and roads, between LiDAR and SRTM DSM data. More than 200 points of bare land and roads were interpreted and digitized on sites 1 and 2. Pearson's r showed that the SRTM DEM and LiDAR DTM had an almost one-to-one relationship ( Figure 2). There were some expected differences in the z or elevation values between the SRTM DSM and LiDAR DTM as these were products of different remote sensing technologies on different platforms (Table 1).    As our study emphasized forest AGB estimation, we focused on the calibration of canopy height values of the SRTM data using the LiDAR CHM. We resampled the SRTM DSM data (30 m) at 1, 5, and 10 m pixel resolutions to match the LiDAR data. A CHM for SRTM was derived by subtracting the LiDAR DTM from the SRTM DEM (SRTM CHM). As historical field data were unavailable, we used LiDAR CHM values of the old-growth forest plots for correlation analysis and linear regression analysis to derive an equation to calibrate the SRTM CHM values. The old-growth forest locations were checked in the field and confirmed with the local villagers to make sure no disturbance activities had occurred prior to our study. Table 2 shows correlations between LiDAR CHM mean and both SRTM CHM and mean height using the old-growth forest plots. As expected, the resampling of SRTM CHM to finer resolution almost did not affect the correlation with mean height. Correlations between LiDAR CHM mean and both mean height and SRTM CHM remained strong (r = 0.860-0.918) up to 10 m resolution, but resampling to 30 m resolution led to a great reduction in the Pearson's r values. Table 3 shows the calibration equations of SRTM CHM for all resolutions. All calibration models had fair to strong goodness-of-fit values (R 2 = 0.740-0.843) except for 30 m resolution (R 2 = 0.467). The calibrated SRTM CHM had an improved regression fit with R 2 > 0.950 and was well distributed along the x = y line in all spatial resolutions (Figure 3).

AGB Estimation Models
Least-square regression analysis is widely used to derive a model for estimating the AGB of forests. In this study, we employed the regression analysis to develop an estimation model (Equation (6)) using LiDAR CHM mean as the predictor to estimate the AGB of 2012 (AGB 2012). The LiDAR CHM mean values of the plots were extracted for each pixel resolution. For each pixel resolution, the LiDAR CHM mean was regressed against field AGB using the following model: where AGB est is the estimated above-ground biomass in Mg/ha and X 1 is the SRTM CHM cal in meters. β 0 , β 1, and ε are intercept, regression coefficient, and additive error term, respectively. These independent and dependent variables were also natural-log transformed because height is known to have a nonlinear relationship with AGB [47]. We randomly selected 50 plots for the model development, and the remaining 14 plots were used to calculate RMSE and relative RMSE (RMSE divided by the mean AGB and expressed in percent) as independent validation (RMSE Spliting ). Splitting the field data into model and validation data is often done in cross-validation that evaluates goodness-of-fit of regression models [48]. Leave-one-out cross-validation (LOOCV) was also applied to the modeling plots to calculate the RMSE (RMSE LOOCV ). The best model was selected mainly based on modeling efficiency (EF) (Equation (7)) and RMSE of the modeling plots (RMSE Model ). The EF, analogous to R 2 , allows direct comparison of predicted values and observed values [49]. The RMSE LOOCV and RMSE Spliting were calculated for the best model to check for consistency. To account for the bias due to the back-transformation, a correction factor (Equation (8)) was applied when calculating the RMSE values [50].
where n = sample size; y = mean field measured AGB;ŷ i = estimated AGB derived from the linear regression model; y i = field measured AGB i. The best model was then applied to the calibrated SRTM CHM (SRTM CHM Cal ) to estimate the AGB of 2000 (AGB 2000). An AGB change map was generated by subtracting the AGB 2000 from AGB 2012. Table 4 shows the summary statistics of DBH, tree H, wood density, and AGB. Although AGB can be calculated using DBH and/or tree height, an allometric equation that includes the wood density value of each tree is more accurate. Field AGB values calculated using Chave et al. (2014) in site 1 ranged between 41.26 and 901.29 Mg/ha. The lowest field AGB value at site 2 was 94.52 Mg/ha, almost twice the minimum value at site 1. However, the maximum field AGB at site 2 was only half the maximum value at site 1, which was 468.04 Mg/ha.

AGB Estimation
The regression analysis results showed that the field AGB could be estimated with LiDAR CHM mean as the predictor for all spatial resolutions except 30 m, which had a very low EF (Table 5). Overall, the model with LiDAR CHM mean at 10 m resolution predicted the field AGB with the lowest RMSE Model value at 74.04 Mg/ha (Relative RMSE Model = 25.79%), but the model with LiDAR CHM mean at 1 m resolution had the highest EF (EF = 0.798). The 1 m resolution model was selected as the best model because it performed better than the other resolution models in explaining the AGB variance, while its relative RMSE Model (25.81%) was only 0.02% higher than the 10 m resolution model. Validation of the best model using the leave-one-out cross-validation (LOOCV) method and data splitting yielded a slightly higher RMSE value of 77.10 (Relative RMSE LOOCV = 26.86%) and 83.50 Mg/ha (RMSE Splitting = 32.41%), respectively (Figure 4). The best model was applied to the SRTM CHM to estimate the AGB 2000.

AGB Estimation
The regression analysis results showed that the field AGB could be estimated with LiDAR CHM mean as the predictor for all spatial resolutions except 30 m, which had a very low EF (Table 5). Overall, the model with LiDAR CHM mean at 10 m resolution predicted the field AGB with the lowest RMSEModel value at 74.04 Mg/ha (Relative RMSEModel = 25.79%), but the model with LiDAR CHM mean at 1 m resolution had the highest EF (EF = 0.798). The 1 m resolution model was selected as the best model because it performed better than the other resolution models in explaining the AGB variance, while its relative RMSEModel (25.81%) was only 0.02% higher than the 10 m resolution model. Validation of the best model using the leave-one-out cross-validation (LOOCV) method and data splitting yielded a slightly higher RMSE value of 77.10 (Relative RMSELOOCV = 26.86%) and 83.50 Mg/ha (RMSESplitting = 32.41%), respectively (Figure 4). The best model was applied to the SRTM CHM to estimate the AGB 2000. Table 5. Aboveground biomass (AGB) estimation models using LiDAR CHM of different spatial resolutions.

Above-Ground Biomass Changes
AGB change maps were derived from the subtraction of AGB 2012 from AGB 2000 for sites 1 and 2 ( Figures 5 and 6, respectively) to assess the AGB changes over the 12 years. Both an increase and a decrease in AGB could be observed in the study area over the 12 years. AGB decrease (gradation of green color) ranged from 41.26 to 673.27 Mg/ha with a mean decrease of 390.49 Mg/ha (±133.53 Mg/ha). A decrease in AGB occurred in both the managed forests and village areas, but the highest decrease was found in the managed forests in which commercial logging is conducted. On the other hand, the AGB increase during the study period ranged between 21.60 and 901.30 Mg/ha. The mean AGB accumulation was 125.29 Mg/ha (±39.82 Mg/ha) over the years or an annual rate of 10.44 Mg/ha/yr (Table 6). In the managed forests, the mean AGB increase was 122.57 Mg/ha or an annual rate of 10.21 Mg/ha/yr. In contrast, the AGB increase in the village areas was almost 20% lower than

Above-Ground Biomass Changes
AGB change maps were derived from the subtraction of AGB 2012 from AGB 2000 for sites 1 and 2 ( Figures 5 and 6, respectively) to assess the AGB changes over the 12 years. Both an increase and a decrease in AGB could be observed in the study area over the 12 years. AGB decrease (gradation of green color) ranged from 41.26 to 673.27 Mg/ha with a mean decrease of 390.49 Mg/ha (±133.53 Mg/ha). A decrease in AGB occurred in both the managed forests and village areas, but the highest decrease was found in the managed forests in which commercial logging is conducted. On the other hand, the AGB increase during the study period ranged between 21.60 and 901.30 Mg/ha. The mean AGB accumulation was 125.29 Mg/ha (±39.82 Mg/ha) over the years or an annual rate of 10.44 Mg/ha/yr (Table 6). In the managed forests, the mean AGB increase was 122.57 Mg/ha or an annual rate of 10.21 Mg/ha/yr. In contrast, the AGB increase in the village areas was almost 20% lower than that in the managed forests (99.70 Mg/ha or 8.31 Mg/ha/yr).   The time series AGB estimates allowed us to examine the pattern and driver of AGB changes in the study area. Most mountain communities in Sabah depend on slash-and-burn cultivation for their livelihoods, and this transforms the vegetation into a secondary forest landscape [55]. Nevertheless, AGB removal by the local villagers is likely to be small-scale and limited to the immediately accessible areas. The main driver of AGB removal identified in this study was timber harvesting in the managed forests. The mean AGB decrease was 158.97 up to 673.30 Mg/ha, suggesting that timber extraction occurred in both logged-over and primary forests. The SFI confirmed that all compartments except compartment P50 at sites 1 and 2 were commercially logged. These logging activities started in 1997 and continued until the early 2000s when the state's forestry policy changed from conventional logging to sustainable forest management [56].
Overall, there was a significant gain of AGB in the previously disturbed areas as a result of natural regeneration, especially in the managed forests. In a regenerating forest, small trees with DBH less than 10 cm can store up to 50% of the total AGB [57]. The mean AGB accumulation rate estimated in this study was 10.44 Mg/ha/yr over the 12 years. Our estimates match those reported elsewhere in the tropics. In a secondary tropical montane forest in southern Ecuador, the annual rate of AGB increase was estimated at 10 Mg/ha/yr [58]. In the montane forest of Mount Kinabalu, Sabah, the mean annual rates of aboveground net primary productivity range widely from 1.60 to 9.44 Mg C/ha/yr with a mean of 4.30 Mg C/ha/yr [59]. Figures 5 and 6 showed that some areas had AGB increases of more than three standard deviations (>20.4 Mg/ha/yr), which were likely to be

Discussion
This study investigated the use of LiDAR CHM and the canopy height information in SRTM data to estimate the AGB changes. Remotely sensed digital elevation data from different sensors, LiDAR and SRTM CHMs, correlate well with the mean canopy height of the forest [23,27]. In our study, we found a strong correlation between the LiDAR CHM mean and the mean tree height. LiDAR CHM mean has been widely used to estimate AGB in tropical forests [36,51,52]. This suggests that SRTM CHM can be calibrated using LiDAR CHM for AGB estimation if these CHMs are strongly correlated. In this study, the LiDAR CHM mean was strongly correlated with SRTM CHM except for 30 m resolution, thus allowing the calibration of SRTM CHM for spatial resolutions of 1, 5, and 10 m using the LiDAR CHM mean. While we found a high goodness-of-fit for the calibration models of the CHMs (R 2 = 0.740-0.843) in the tropical montane forest, past studies indicate that the regression fit might depend on forest type. Notable regression fits were reported for the CHMs in mangrove and boreal pine forests [27,53] but not for the CHMs in boreal hardwood forests [53]. In this study, the calibrated SRTM CHMs showed an improved regression fit when plotted against the LiDAR CHMs (R 2 ≥ 0.959) compared to the uncalibrated ones and became well distributed along the x = y line (Figure 3).
In the tropical montane forest, the LiDAR-CHM mean was examined as the estimator of field AGB for different spatial resolutions. The AGB estimation model using 1 m resolution LiDAR CHM had the best fit and explained the AGB variance better than other spatial resolutions in the regression analysis. The model EF decreased with decreasing the spatial resolution of the LiDAR CHM. When the LiDAR CHM resolution became coarser, the extracted information became more uncertain, thus leading to a decrease in spatial resolution. When the pixel size increased, the small forest canopy gaps could be overlooked, which could have led to an overestimation of the forest canopy height [54]. The best model using LiDAR CHM mean at 1 m resolution had a relatively low estimation error (RMSE LOOCV = 77.10 Mg/ha, corresponding to about 26.86% of the average AGB). By randomly repeating the sample partitioning, we found that the relative RMSE of the independent validation ranged between 25.89% and 34.54%, with an average of 29.84%. Overall, these RMSE values were in a similar range to the RMSE of a multivariate AGB estimation model described by Ioki et al. (2014).
The time series AGB estimates allowed us to examine the pattern and driver of AGB changes in the study area. Most mountain communities in Sabah depend on slash-and-burn cultivation for their livelihoods, and this transforms the vegetation into a secondary forest landscape [55]. Nevertheless, AGB removal by the local villagers is likely to be small-scale and limited to the immediately accessible areas. The main driver of AGB removal identified in this study was timber harvesting in the managed forests. The mean AGB decrease was 158.97 up to 673.30 Mg/ha, suggesting that timber extraction occurred in both logged-over and primary forests. The SFI confirmed that all compartments except compartment P50 at sites 1 and 2 were commercially logged. These logging activities started in 1997 and continued until the early 2000s when the state's forestry policy changed from conventional logging to sustainable forest management [56].
Overall, there was a significant gain of AGB in the previously disturbed areas as a result of natural regeneration, especially in the managed forests. In a regenerating forest, small trees with DBH less than 10 cm can store up to 50% of the total AGB [57]. The mean AGB accumulation rate estimated in this study was 10.44 Mg/ha/yr over the 12 years. Our estimates match those reported elsewhere in the tropics. In a secondary tropical montane forest in southern Ecuador, the annual rate of AGB increase was estimated at 10 Mg/ha/yr [58]. In the montane forest of Mount Kinabalu, Sabah, the mean annual rates of aboveground net primary productivity range widely from 1.60 to 9.44 Mg C/ha/yr with a mean of 4.30 Mg C/ha/yr [59]. Figures 5 and 6 showed that some areas had AGB increases of more than three standard deviations (>20.4 Mg/ha/yr), which were likely to be overestimated. As LiDAR is currently the most accurate remote sensing method for estimating AGB, the overestimation was probably due to underestimation of the AGB 2000 values, especially in rugged terrain areas, because the accuracy of SRTM DSM of steep-sloped areas can be significantly affected by slope and aspect characteristics in relation to the incidence angle [60].

Conclusions
In the tropics, monitoring of AGB changes has been challenging due to a lack of historical AGB information. We examined the use of spaceborne (SRTM) and airborne (LiDAR) digital elevation data for estimating AGB changes in a tropical montane forest in northern Borneo between 2000 and 2012. Coupled with field data, an AGB estimation model was developed using the 1 m resolution LiDAR data. Not only was the AGB of 2012 estimated accurately, but the model can also be applied to estimate the historical AGB of 2000 by resampling and calibrating the SRTM DSM data. With the AGB estimates, it was possible to examine the AGB changes in two different land uses over the 12 years. As SRTM DSM data are freely available, this data can be used to estimate historical forest AGB as baseline information for monitoring forest changes in other tropical regions and to determine the main driver of AGB changes.