MODELLING THE CARBON STOCKS ESTIMATION OF THE TROPICAL LOWLAND DIPTEROCARP FOREST USING LIDAR AND REMOTELY SENSED DATA

Tropical forest embraces a large stock of carbon in the global carbon cycle and contributes to the enormous amount of above and below ground biomass. The carbon kept in the aboveground living biomass of trees is typically the largest pool and the most directly impacted by the anthropogenic factor such as deforestation and forest degradation. However, fewer studies had been proposed to model the carbon for tropical rain forest and the quantification still remain uncertainties. A multiple linear regression (MLR) is one of the methods to define the relationship between the field inventory measurements and the statistical extracted from the remotely sensed data which is LiDAR and WorldView-3 imagery (WV-3). This paper highlight the model development from fusion of multispectral WV-3 with the LIDAR metrics to model the carbon estimation of the tropical lowland Dipterocarp forest of the study area. The result shown the over segmentation and under segmentation value for this output is 0.19 and 0.11 respecti vely, thus D-value for the classification is 0.19 which is 81%. Overall, this study produce a significant correlation coefficient (r) between Crown projection area (CPA) and Carbon stocks (CS); height from LiDAR (H_LDR) and Carbon stocks (CS); and Crown projection area (CPA) and height from LiDAR (H_LDR) were shown 0.671, 0.709 and 0.549 respectively. The CPA of the segmentation found to be representative spatially with higher correlation of relationship between diameter at the breast height (DBH) and carbon stocks which is Pearson Correlation p = 0.000 (p < 0.01) with correlation coefficient (r) is 0.909 which shown that there a good relationship between carbon and DBH predictors to improve the inventory estimates of carbon using multiple linear regression method. The study concluded that the integration of WV-3 imagery with the CHM raster based LiDAR were useful in order to quantify the AGB and carbon stocks for a larger sample area of the Lowland Dipterocarp forest. * Corresponding author


INTRODUCTION 1.1 Aboveground Biomass and carbon stocks in relation to Climate Change
Forests are crucial for human life.They provide important resources for medicine, foods, habitat for animal and also play an essential role in carbon sequestration (Abd Latif et al., 2011).In the past few centuries, anthropogenic activities and natural consequences has increased the concentration of carbon dioxide including greenhouse gases (GHG) to the atmosphere (FAO, 2015).In year 1992, the United Nation Framework Convention on Climate Change (UNFCC) formed a framework to decrease the GHG conveyed at the United Nations Conference on Environment and Development (UNCED) held in Rio de Janeiro (UNFCC, 2015).In 2007, the Bali Climate Change Conference under United Nations had adopted Kyoto protocol by setting rightfully binding obligation for reduction of GHG gases emission (UNFCC, 1998).During this conference, an important agreement was reached for the developing countries to initiate actions in reducing emissions from deforestation and forest degradation (REDD) (Pelletier et al., 2012).It sets a collective global target of reducing GHG emissions by about 5% of the 1990 levels by the first commitment period from 2008 to 2012.An amendment made to the protocol in Doha, Qatar, on 8 December 2012, known as the "Doha Amendment to the Kyoto Protocol", aims to reduce GHG emissions by at least 18% below 1990 levels from 2013 to 2020 during the second commitment period (UNFCC, 2015).In continuation to the persistent efforts to combat global warming and restrain associated global risks, developing a method is increasingly important for monitoring, reporting and verification (MRV) mechanism of Reducing Emissions from Deforestation and Forest Degradation (REDD+).
In South East Asia generally, due to the economic factors, it has been historically proven that a lot of primary forest had been cleared for timber extraction and extensive conversion to oil palm plantation, and rubber plantation (Hansen et al., 2013).It is confirmed in a recent forest cover research that shows the main cause of forest loss in South East Asia is the conversion of the forest to cash crop plantation (Stibig et al., 2014).Extensive deforestation and conversion have significant impacts from local to regional scales on wide ranges of ecosystem offered by forest such as hydrological and habitat provision (Blackburn et al., 2014).Consequently, many habitats of plants and animal are lost and it impacts the biodiversity and the richness of species contained in the forests.

Aboveground Biomass for tropical rain forest.
Estimating the aboveground biomass of the forest, and investigating the carbon stock for the tropical rainforest are not easy tasks due to the complex architecture of trees and plants that vary in sizes (Mohd Zaki et al., 2015;Mohd Zaki & Latif, 2016).In order to predict the tree structure parameter using integration of field data and remote sensing technology, models are normally developed by previous researchers (Babcock et al., 2015;Gonzalez et al., 2010;Li et al., 2015;Saeidi et al., 2014).With the present complex structure of tropical forest, there is an enormous uncertainty for carbon stocks estimation.(Figure 1).In order to measure the aboveground biomass and carbon stocks of an area, the fusion of small footprint discrete Light Detection and Ranging (LiDAR) and Very high resolution (VHR) WV-3 (30cm spatial resolution) data had been used to precisely determine the AGB/C stocks for the dominant tree of the study area.Primarily, there are six dominant family species which are Dipterocarpaceae, Euphorbiaceaea, Sapotaceae, Burseraceae, Moraceae, and Lamiaceae.Among them, one of the species is categorized as Critically Endangered (CR) by (IUCN, 2014) which is Hopea Sulcata or Merawan Meranti.(Turner, 1996) The model conveyed in this study aims at filling the knowledge gap by assess modelling the carbon stocks of lowland Dipterocarp forest for Malaysia cases.The models for aboveground biomass and carbon stocks needed for forest biomass with the remote sensing environment, and the option that geospatial technology, can be provided in order to support effective forest management.Therefore, the objective of this study includes (i) investigating the relationship between carbon stocks and tree parameter (DBH, tree height, and crown diameter), (ii) developing a model of carbon stocks estimation using Crown Projection Area (CPA) of the emergent and canopy layer of the tree crown and height derived from LIDAR, (iii) producing the map of carbon stocks generated from the model and (iv) predicting the carbon estimation for the broader area of Ayer Hitam Forest Reserve.

Study Area
The study area is located at Latitude 3°00'24.19"N,Longitude 101°38'25.24"E in the Ayer Hitam Forest Reserve (AHFR), Selangor State, Malaysia (Figure 2).The data collected from this lowland Dipterocarp forest of Ayer Hitam Forest Reserve chosen for field sampling are owned by University Putra Malaysia (UPM) under custodian of Faculty of Forestry.This secondary forest comprises of various species that are dominant by family tree of Dipterocarpacaea.The altitude that comprises in this lowland forest varies from 15 m to 233 m height, and the terrain slope undulating up to 34º (Shida et al., 2014).The average annual rainfall is 2178 mm while the average temperature annually is 25.3ºC with maximum 27.7º C and minimum 22.9º C (Shida et al., 2014).

Airborne Laser Scanning (ALS) Data and WV-3 data
The ALS survey had been flown using a Eurocopter 120 on August 2013.The LiteMapper-Q560 that consist of RIEGL LMS-Q560 laser scanner for LiDAR scanner was mounted on the aircraft, along with the Hassleblad digital camera 39 megapixel with 50mm focal length.Table 1 summarizes the flight parameter used for capturing LiDAR to achieve point density of 11.98 per square units.The projection used is WGS1984 UTM Zone 47N.The ALS survey of raw files was pre-processed using AeroOffice IGI to analyse the quality of IMU and GPS data.Terrasolid tools were used to process the ALS dataset (Bentley Systems, 2015) into LIDAR 3D point clouds.Super-spectral high resolution of WV-3 (8 bands) 25 km 2 with the projection of UTM47 N and WGS84 datum, have been used in this study.The data was obtained on 9th December 2014 that consists of 31cm panchromatic band, 4 bands VNIR colours (blue, green, red and NIR-1) and 4 added VNIR colours which are coastal, yellow, red edge and near-IR2.The spatial resolution was 0.30 m (panchromatic) and 1.2 m (multispectral).The classification was performed on a pansharpen WV-3 image using hyperspherical colour sharpening (HCS) techniques for better visualization.

Field sampling plot
Field sampling had been collected using stratification number of sampling, estimated using equation ( 1).Tree stand sampling and forest mensuration data was collected on May 2015 at 2 ha rectangular that consists 50 subplots included tree id, species, tree ages, stem diameter at breast height (1.3) for live trees which have 10 cm and above, sample of tree height (living tree), crown base height (the height for 1 st branches up to the crown), leaf area index and also crown diameter of the canopies.
(1) Nplot = Minimum number of sampling plots t 2 = Value of tree distribution of Nplot CV = Coefficient of Variance (%) E = Allowable error Field plot and tree coordinates were located by traversing using total station to the plot edge coordinate.Each of coordinate were tied by Global Positioning System (GPS) observation of Topcon GTS using static observation L1/L2 and L5 bands.The positional accuracy of the plot coordinates was generally within ± 0.5 m of the acceptable tolerance.Single tree locations that consist of coordinate, bearing and slope distance were measured for each individual tree in relation with the centre plots for the future use.In all plots, the tree height of forest stands were measured using digital DISTO D5 laser ranger while the DBH (1.3 m above ground height) was measured using a DBH tape.

Figure 3. Methodology workflow
This study was performed according to three major phases.This comprises of the data collection of the forest stands, data acquisition, pre-processing of the data, object based image analysis (OBIA), model development, validation of the data and finally the statistical analysis.Figure 3 illustrate the process involved in this study.

Canopy Height Model generation:
Canopy height model is commonly calculated in the form of raster representation, generate from the highest echo of first return laser pulse known that create the digital surface model (DSM) subtract with the lowest echo of last return which is known as digital terrain model (DTM) (Popescu et al., 2003).
Normalization is applied to the CHM which gives the absolute point of height of tree height on Z-axis, and according to the dataset, the maximum tree was 44m.Thus, we drop all the LiDAR points that are higher than the value of maximum tree, as a noise.The CHM was created using LAStools software in ArcGIS plug-in.Figure 4 show the CHM of the tree height in 3Dimension (3D) view.On the basis, carbon were converted by applying the conversion factor of 0.47 which represents 47% of the dry biomass assumed to be carbon for all part of the tree as the default value that had been recommend by IPCC (IPCC, 2006).

Crown Projection Area (CPA) delineation and validation:
There are many approaches in calculating the accuracy assessment of segmentation based on literature but the most common method used is by visual (Mohd Zaki et al., 2015)Previous researchers have used different type of segmentation validation by overlapping the area that intersect between output segmented area with reference area (Möller et al., 2007).Another researcher used the distance between two centroids to assess the segmentation accuracy (Ke, Quackenbush, & Im, 2010).On the other hand , Gougeon (1995) had developed the segmentation validation using 1:1 spatial corresponding based on goodness of fit (D).However, Clinton et al., (2010) had improve the segmentation validation method developed by Ke et al., (2010) and Möller et al., (2007) by modifying the relative area metrics by calculating Equation ( 3) and ( 4) and measuring the goodness of fit by using Equation (5).Segmentation accuracy is measured based on distance index (D) ranging from 0 to 1, where 0 shows that the classification is an ideal match between xi and yi, while 1 is the minimum mismatch.In order to evaluate the multi-resolution segmentation, this study has applied evaluation approach by Clinton et al., (2010) where it measures the segmentation goodness measures.

Statistical analysis:
The correlation between CS stocks, height, DBH and CPA were analysed using Pearson's correlation as well as multiple linear regression model.Bootstrap sampling was executed in order to account for the small value of tree stands.The statistical analysis is proceeded by performed by calculating statistical parameter such as residual R 2 , adjusted R 2 , Root Mean Square Error (RMSE) (Equation 7), segmentation validation and height validation.The multi-collinearity test also had been done in order to avoid any collinearity problem between variables, with a tolerance cut of value of variance inflation factor (VIF) should be less than 10 (O'Brien et al., 2007) (Equation 6).

….. (6)
Where is the squared multiple correlation of the variables of jth explanatory variable in the regression model (Equation 6).During the classification of the segmentation of the individual tree, one to one matching of reference and segmentation tree crowns was used for developing the model.Then, the model were validated with 30% of field sampling data out of 911 number of tree stands that had been measured.

…….. (7)
Where n is the number of samples, yi is the calculated value of carbon, is the predicted carbon by model of the response variables.

ANALYSIS AND DISCUSSION
The analysis of this research comprises of the description of studies, normality, and correlation analysis.The following paragraphs will explain the results of these analyses.

Descriptive statistics
Overall, the amount of tree sampling of the whole study area at Ayer Hitam Forest Reserve was 911 trees collected from 50 subplots within 2 ha plots.The accuracy of GPS observation at the base station was 3 mm for horizontal (x and y) and 5mm for vertical (z) plot.In total, 911 of tree had been measured from the observation where six dominant genuses had been identified: Dipterocarpaceae, Euphorbiaceaea, Burseraceae, Sapotaceae, Moraceae, Rhizophoraceae and Rubiaceae.Dominant species nominated by Endospermun diadenum, Hopea Sulcata, Shorea Macroptera and Dipterocarpus Verucosus.In order to ensure representative, 183 of tree were randomly selected in each subplot, while the other 62 trees (30%) were used as a reference (Figure 5).crown projection area derived from the fusion of LiDAR and multispectral WV-3.Based on the finding, this shows the output of the Kolmogorov-Smirnov and the ShapiroWilk tests for the distribution of DBH, height (field), height (LiDAR), and CPA.For both test, p-value was nearly normal distributed (P > 0.01) (Table 2).The log transformation had been applied to make sure that the data is normally distributed for the all dependent variables thus to make a perfectly linear regression (Piaw, 2014).

Validation of LiDAR derived height
With regards to the analysis, the highest positive r value was recorded for height LiDAR versus height field which is 0.995 (Figure 6) which shows that there is a strong correlation between height LiDAR and height field.Besides, there is no significant difference (p > 0.01) between height measured from field and height derived from LiDAR.However, there is a strong correlation between height measured from field and height derived from LiDAR.Table 3 and Table 4 summarize the output result of the compare means between heights LiDAR versus height field of 245 sample of tree.
Figure 6.The relationship between tree height from field (Ht_m) and tree height derived from LiDAR (Ht_LIDAR).

Accuracy assessment of segmentation output
Accuracy assessment of the tree crown segmentation had been applied in two techniques, which is i.)By using techniques 1:1 matching of the references and segmentation polygon, ii.) Quantify the goodness of fit (D) which represents the ideal segmentation.Based on Table 5, it is shown that the average accuracy of the segmentation is 79.5%.Out of 288 manually delineated of the crown, only 266 were matched 1:1 (78% match).The over segmentation and under segmentation value for this output is 0.19 and 0.11 respectively ,thus D-value for the classification is 0.19 which is 81% .The result shows that the segmentation output are acceptable for further analysis.(Figure 7).

Correlation analysis -Pearson correlation
In an effort to see the relationship between these variables, for example, between DBH and height LiDAR, height field and CPA, Pearson correlation had been done to measure the correlation between independent variables of carbon estimation.From the analysis, the highest positive r value was recorded for H_LDR versus H_Field which is 0.995 which shows that there is a strong correlation between height LiDAR and height field (P < 0.01).As a result, there is strong evidence that suggests that this correlation is significant.The positive correlation between Ht_LiDAR and DBH showed in Table 6, which depicted that the high correlation between this two variables (0.763) with (P < 0.01).The results show that there is a high correlation between height LiDAR and DBH.

Variable
Pearson correlation analysis of the independent variables is used to calculate the carbon stocks where carbon as the dependent variable.High positive relationship between DBH and carbon stocks was found with (P < 0.01) (Table 6).The data yielded by this study provides strong evidence.Multiple regression model was developed using LiDAR derived height and CPA as independent variable and carbon as a dependent predictors.(See Table 7, Table 8 and Table 9)

The multiple linear regression analysis
Based on the regression, there was a negative linear relationship between predictor and output variables, β = -2971.429,ρ = 0.000 (p < 0.01).Multi-collinearity between the two variables had been checked by applying the variance inflation factor (VIF) test where the result shows that there was no sign of multi-collinearity found between these two variables (VIF= 1.363).By using multiple regression model to predict the carbon stocks estimation over the height derived from LiDAR and crown projection area was a negative predictors, β = -2971.429,p = 0.000 (p < 0.01) disputed that every unit increase of carbon stocks, the height from LiDAR and CPA were decreased by -2971.429m.Thus, the relationship of carbon stocks can be written as the following equations: … (7) .......( 8) The coefficient demonstrated in Table 9 represents the effect size, beta = 0.726.Moreover, the result shows that the Durbin-Watson value is 1.502 which exemplifies that autocorrelation is within the tolerance where the Analysis of variance (ANOVA) test indicates that the overall model are significant, p = 0.000 (p < 0.01).

Validation of multiple linear regression model
Pearson's productmodel of tree crown was validated using a random selection of 30% of the independent dataset which is in this case 62 trees from various types of tree height and species.The resulting coefficient of determination between carbon predicted and the carbon observed was 0.6737 (R 2 = 0.6737) (Figure 8).

Carbon Stocks Map
In brief, the range of the predicted carbon is around 800 kg/trees to 7634 kg/tree where it can be concluded that the higher the CPA and height value, the greater the amount of carbon it would make.Figure 9 shows that the darker colour of the CPA which indicates the highest value of carbon, while the white colour of the background of the map is the gap of each individual tree.The carbon stocks map of the study area is shown in (Figure 9).
Figure 9. Carbon stocks of the study area (Karna (2013) found out that a significant correlation coefficient (r) between (CPA -Carbon), (height -Carbon), and (CPA -Height) is 0.73, 0.76 and 0.63 respectively.Similar to this relationship, this study produced 0.671, 0.709 and 0.549 for the type of Lowland Dipterocarp forest which mainly focuses on tropical rain forest.However, the highest correlation exemplified by the relationship between DBH and carbon stocks which is Pearson Correlation (0.909) shows that there is a good relationship between carbon and DBH predictors to improve the inventory estimation of carbon using multiple linear regression method.

CONCLUSION AND RECOMMENDATION
This paper presents a model development from fusion of CHM derived from LiDAR and WorldView-3 imagery using the MLR-based methodology to estimates the forest biomass from remote sensing technology.The combination of LiDAR and Very High Resolution Multispectral Imagery, WV-3 has demonstrated a promising capability to model the aboveground biomass and carbon stocks needed for forest biomass estimation for lowland Dipterocarp forest.Results indicate that the relationship between carbon stocks with LiDAR and CPA obtained in this study are similar in terms of correlation produce and the levels of variances with other studies that have been done previously (Karna et al., 2013).The output MLR shown that there is non-linear equation of the predicted model where the natural logarithmic had been used and get the better prediction.Future works should address the gaps that did not cover in this paper.The primary limitation of the methods is the robust of the complexity intermingles of the tropical forest itself and the improvement of the CPA delineation should be address.Finally, the synergistic use of the non-parametric test comparing with this traditional regression analysis, MLR will be the future works in order to estimate the biomass and carbon stocks of the tropical rainforest.

Figure 2 .
Figure 2. A map of the location of the study area.Figure 2a.Shows the location of Selangor District, at Peninsular Malaysia.Fig. 2b shows a location of Ayer Hitam Forest Reserve Lowland Dipterocarp forest covered by the WV-3 satellite image data.
Figure 2. A map of the location of the study area.Figure 2a.Shows the location of Selangor District, at Peninsular Malaysia.Fig. 2b shows a location of Ayer Hitam Forest Reserve Lowland Dipterocarp forest covered by the WV-3 satellite image data.

Figure 4 .
Figure 4. Canopy height model.Figure 4a.The top view of the CHM of the study area.Fig. 4b.The 3D view from the top of CHM, Fig. 4c.The interpolation of height of the CHM and finally Fig. 4d.The perspective view of the CHM of the study area.2.4.2Estimation Carbon stocks using allometric equation: AGB estimation were calculated using allometric equations (Chave et al., 2014) (Equation2) that uses wood density, height and DBH as predictors.Wood density is an important predictor in the calculation of AGB based on species despite of the height and DBH(Fayolle et al., 2013).The details of each tree were recorded as: total height ht, bole height hb, crown diameter ci, diameter at 1.30m height di and species.

Figure 5 .
Figure 5. Box plot of data collection from field data (DBH, height, crown diameter and bole height).Fig. 5a.Show the boxplot of DBH of the tree stands, Fig. 5b.The height of the tree, Fig. 5c.The crown diameter of the tree stands and Fig. 4d.The bole height which is the height from of the first branch.3.2Normality test for model developmentThe selection predictor variables had analysed the normality test used for the model development.The result represents the distribution of DBH, height from field, height from LiDAR and over segmentation, US = under segmentation and D = Goodness of fit) . a. Dependent Variables: Carbon Stocks b.Predictors: (Constant), H_LDR and CPA

Figure 8 .
Figure 8. Scatter plot of multiple regression model validation

Table 2 .
Normality test for model development.

Table 3 .
Summary of statistical -Tree height from field and LiDAR data

Table 4 .
Overall summary of fit

Table 5 .
Accuracy assessment of OBIA output Figure 7.The segmented (Red colour) and Reference (Yellow) polygon

Table 6 .
Pearson correlation of the DBH, H_LDR and CPA.

Table 7 .
Model summary

Table 9 .
Coefficients a