The Interpolation Method for Estimating the Above-Ground Biomass Using Terrestrial-Based Inventory

This paper examined several methods for interpolating biomass on logged-over dry land forest using terrestrialbased forest inventory in Labanan, East Kalimantan and Lamandau, Kota Waringing Barat, Central Kalimantan. The plot-distances examined was 1,000−1,050 m for Labanan and 1,000−899 m for Lamanda. The main objective of this study was to obtain the best interpolation method having the most accurate prediction on spatial distribution of forest biomass for dry land forest. Two main interpolation methods were examined: (1) deterministic approach using the IDW method and (2) geo-statistics approach using Kriging with spherical, circular, linear, exponential, and Gaussian models. The study results at both sites consistently showed that the IDW method was better than the Kriging method for estimating the spatial distribution of biomass. The validation results using chi-square test showed that the IDW interpolation provided accurate biomass estimation. Using the percentage of mean deviation value (MD)(%), it was also recognized that the IDWs with power parameter (p) of 2 provided relatively low value , i.e., only 15% for Labanan, East Kalimantan Province and 17% for Lamandau, Kota Waringing Barat Central Kalimantan Province. In general, IDW interpolation method provided better results than the Kriging, where the Kriging method provided MD of about 27% and 21% for Lamandau and Labanan sites, respectively.


Introduction
Since clean development mechanism (CDM) and reducing emissions from deforestation and forest degradation (REDD +) were launched, issue on methods for vegetation biomass estimation is becoming increasingly important.Having the second largest tropical forests in the world, Indonesian forestry sector receives much international attention, particularly related to deforestation and forest degradation.In Asian countries, Indonesia contributes the biggest carbon stocks of tropical forests.Carbon densities and pools in vegetation and soil varied widely by eco floristic zone and country.
In Indonesia, in line with the framework of REDD+ mechanism as mentioned above and the "one map policy" that has been declared by the President of the Republic of Indonesia in 2010, estimation and mapping of vegetation biomass as well as their carbon stocks for each forest ecosystem has been a fundamental task.Map of biomass is required to estimate the total carbon stock at any given time and to quantify the amount of emission reduction that has been achieved through the implementation of REDD+ activities.Such map can be used to support the government's commitment in reducing green house gasses (GHG).The Government of Indonesia had committed to voluntary reduce GHG emissions up to 26% with its own fund or 41% with international assistance by 2020 (Perpres Number 61 of 2011).
Various techniques are readily available that can be used to estimate forest biomass, either using terrestrial or remote sensing techniques.Currently, Indonesia possesses various data sources that are considerably potential to be used in estimating forest biomass.One such data sources can be derived from the periodically performed forest inventory, often referred to as inventarisasi hutan menyeluruh berkala-comprehensive and periodic forest inventory (IHMB) performed by each forest concession holders (IUPHHK) both in natural forests and plantation forests.Although the ultimate goal for data collection of IHMB is to generate data on standing stocks as the basis to formulate the concession's work plan, the available data has the potential use to estimate the spatial distribution of forest biomass and carbon content.In this study, the use of IHMB data is studied to produce a map of forest biomass.
Presenting the spatial information of forest biomass either in the form of a digital map or a paper map is an important task to support sustainable forest management, particularly for the management of forest carbon "payment/offset".However, to obtain the appropriate method with accurate estimation, a spatial interpolation technique development in line with the ecological characteristics of each forest ecosystem is needed.To date, there are few researches were related to biomass estimation using terrestrial data, and almost no research related to biomass estimation using IHMB data in Indonesia.Tojo et al. (2008) examined the powerful use of interpolation technique for estimating the biomass distribution in the Pacific Coast.The main questions of this study were (1) what are the characteristics of spatial distribution within each forest ecosystem type, for example, within the dry forest, swamp forest, peat swamp forest or mangrove forest?; (2) is the spatial distribution pattern of old forest biomass similar to the newly logged over forests or disturbed logged-over forest or degraded forest due to fires?
In addition, geostatistics and deterministic methods are useful tools to describe and map the spatial variability and estimation of forest variables.Therefore, this research was conducted to investigate which method is suitable to be applied in ecosystem with high spatial variability and secondary forest with high dynamic changes in East and Central Kalimantan.Currently, geostatistical methods have found their applications in forestry.Biomass mapping of mangrove forest has been explored by Simard et al. (2006), while Caspian forest has been examined by Akhawan et al. (2010).Du et al. (2010) use geostatistics theory to examine the spatial heterogeneity of the above-ground biomass (AGB) of bamboo and apply a point Kriging interpolation method to estimate and map its spatial distribution.
At present, there are several approaches that can be used to estimate the distribution of AGB level.With terrestrial approach, the AGB distribution can be estimated using the results from terrestrial measurements, i.e., census and sampling methods, or by using passive or active remote sensing technology.Previous researchers (Mulyanto & Jaya 2004;Septriana et al. 2004;Jaya 2005;Jaya & Abe 2006;Gunawan et al. 2010) present a statistical approach that integrates satellite imagery, and ground inventory data to map the spatial distribution of land cover as well as forest cover changes.In this study, the focus will be on the utilization of IHMB data collected by each IUPHHK to estimate the standing biomass distribution, especially in natural forests.
The main objective of this study was to obtain the most accurate spatial interpolation method to estimate the spatial distribution of forest biomass in the dry forest.Additional objective was to determine the spatial distribution patterns of dry land forest stand biomass.

Methods
The study covered 2 sites, site 1 is within PT.Inhutani I concession area located in Labanan District, Berau Regency, East Kalimantan Province and site 2 is within PT TSI concession area located Lamandau District, Kota Waringin Barat Regency, Central Kalimantan Province.These 2 study sites are mainly covered by logged-over dry land forest.
The important step on measuring above-ground biomass is to select sample plots followed by tree measurement, living or dead, necromass, litter and ground-cover.In this study the sample plots were selected using systematic sampling with random start, where distance between plots was 1,000 × 1,050 m for Labanan site and 1,000 × 899 m for Lamandau site (Figure 1 and Figure 2).The sample plot for large tree (dbh > 30 cm) took a rectangular form with dimension of 20 × 125 m.Small tree (dbh between 20 and < 30 cm), pole (dbh between 10 and < 20 cm), and sapling (vegetation height > 1.5−dbh < 10 cm) were measured within 20 × 20 m, 10 × 10 m, and circle with a radius of 2.82 m sub-plot, respectively.
Seedling, ground-cover and litter were measured within 1 × 1 m sub-plot.The main data used in this study were derived from the IHMB data of dry land forest that were collected by the concession holders of PT.Inhutani I, in Labanan, Berau Regency, East Kalimantan and PT.TSI in Lamandau Regency, Central Kalimantan.In this study, the number of sample plots available was divided into 2 parts, one part was used for model development, and the other was used for model validation.The plots for the model and the validation were taken alternately from IHMB data, so that the distance between the path and the distance between the plots on track would be 2 times as long as the distance between plots of IHMB.The illustration of this data sampling is shown in Figure 1 for Labanan and Figure 2 for Lamandau.Since data for either the model or for verification were selected alternately, therefore the distance between plots examined would be 1,450 × 1,450 m for Labanan site and 1,345 × 1,345 m for Lamandau site.
The main issue to be addressed in this study was the AGB of forest.The AGB is defined as all living biomass above the soil including stem, stump, branches, bark, seeds and foliage (IPCC 2006).The biomass content of live standing trees were computed using the equation developed by Kettering et al. (2001), as shown in Equation [ -1 note: B = biomass content (kg tree ), is specific gravity for each tree species and D = tree diameter at breast height.In this study, since all dead materials such as trunk, branches, stump, and litter were included, then total aboveground biomass (AGB ) for each plot was derived by All IHMB The spatial interpolation was carried out with 2 approaches, namely deterministic (non-geostatistical) and stochastic (geostatistical) methods.In this study, 2 interpolation methods, i.e., inverse distance weighted (IDW) and Kriging methods that fall into the first and second categories were examined.

Inverse distance weighting
The first method examined was IDW method that estimated the values of an attribute at non sampled points using a linear combination of values at sampled points weighted by an inverse function of the distance from the point of interest to the sampled points.The assumption is that the values of sampled points closer to the non sampled point were more similar to it than those further away.The weights could be expressed as shown in Equation  The main factor affecting the accuracy of IDW was the value of the power parameter, the smaller the value, the sharper interpolation shape would be obtained.Weights diminished as the distance increased, especially when the value of the power parameter increased, hence nearby samples had a heavier weight and more influence on the estimation, and the resultant spatial interpolation is local.The choice of power parameter and neighbourhood size was arbitrary (Webster & Oliver 2001).The most popular choice of p was 2 and the resulting method was often called the inverse square distance or inverse distance squared (IDS).This method was widely used for interpolating the standing stocks of timber in tropical forest of Indonesia.

Geostatistics
The second method examined was geostatistics.This method was often believed to have originated from the work in geology and mining by Krige (1951), but it could be traced back to the early 1910s in agronomy and 1930s in meteorology (Webster & Oliver 2001).Geostatistics includes several methods that use Kriging algorithms for estimating continuous attributes.Kriging is a generic name for a family of generalised leastsquares regression algorithms, used in recognition of the pioneering work of Danie Krige. Semivariance (γ) of Z between 2 data points is an important concept in geostatistics and is defined as shown in Equation [4]. [4] note: x = the prediction location 0 x = measurement at the ith location i h = the distance between point x −x and γ(h) is the i 0 semivariogram (commonly referred to as variogram) (Webster & Oliver 2001).A plot of γˆ(h) against h is known as the experimental variogram (Figure 3), which displays several important features (Burrough & McDonnell 1998).The first is the "nugget", a positive value of γˆ(h) at h close to 0, which is the residual reflecting the variance of sampling errors and the spatial variance at shorter distance than the minimum sample spacing.The "range" is the value of distance at where the "sill" is reached  where χ < χ at 95% confidence level.
cal table [6] note: The model is acceptable when the AD value lies between the range of -1 and 1. [7] [8] If the value of bias was positive, then the estimation results tent to be overestimated (likely exaggerated estimated value).On the contrary, if the value was negative, then the estimated value or predictions tent to be smaller (underestimated).A better model tent to have bias values close to zero. [9] The formula for RMSE in percentage, RMSE(%) is written as follows: [10] The generic name for measuring the average of absolute difference between estimated and actual value is the mean absolute error (MAE), formulated as follows: Modification of the MAE is percentage of mean deviation (MD%), formulated as follows: [12]    In general, the MD% is categorized as good when the value is less than or equal to 10%.

Results and Discussion
This study analyzed the above-groung biomass stand variables that were converted from the standing stock derived from IHMB data using the equation model produced by Jaya et al. (2012).
Labanan study site Spatial interpolation performed for data of study site 1 showed that the simple IDW method with power value p = 2 gave better estimation than using Kriging method.This finding coincided with the research of Akhawan et al. (2010).The forest cover and stock in logged over forest of Labanan was quite varying, where many gaps due to logging activity, shifting cultivation and fire were found.Therefore, Kriging could not produce accurate estimations due to the high spatial variability of forest growing stocks related variables in this heterogeneous and uneven-aged forest.These results were different from the study performed by Due et al. (2010), which show that Kriging interpolation based on geostatistical theory provides a promising method for estimating AGB.Their results showed that the exponential semi-variance model demostrated the best performance for AGB estimation and for examining its spatial heterogeneity.
As summarized in Table 1, interpolation method using IDW with power p = 0.5 provided a MD of 19.2%.The MD then decreased as the power increased to p =2. Furthermore, the MD would increase gradually when the power p was increased.In general, the MD of IDW showed much better results than those obtained from Kriging's.Figure 4 showed the trends of MD, RMSE, BIAS, Chi-sq.and AD at several values of p and Kriging models.Low errors were consistently shown by IDW method at p = 2.The Kriging method consistently provided more errors than the IDW.
On the basis of score for each interpolation method as tabulated in Table 2, the p = 2 consistently provided very low score, having a total score of 1.3.This may conclude that the best performance was provided by IDW method with p = 2.The highest score was given by Gaussian model (GS) expressing that model was not implementable.Graphically, the total score for each estimation model were shown in Figure 5.
For the geostatistics approach, the effect of "nugget" was quite typical in forest inventory data (Akhawan et al. 2010).The large nugget in the variogram could be explained by 3 causes (Chiles & Delfiner 1999): (a) estimation error of random effects or stand density, (b) Micro-structure that is the component of a range shorter than the sampling plot, and (c) structures with ranges shorter than the smallest interpoint distance.In this study, for practical reason, the plot distance examined was larger than their original value.For Labanan area (site 1), the original distance was 1,000 × 1050 m, while the analyzed plot distance was 1,450 × 1,450 m (note the data were alternately sampled); for the Lamandau (site 2) the original distance was 1,000 × 899 m, while the evaluated plot distance became 1,245 × 1,345 m.Normally, nugget effect should decline as the sampling distance decreases (Akhawan et al. 2010).Inversely, increasing the sampling plot distance would augment the nugget effect.
In the secondary forest, forest growing stocks were mainly scattered randomly over the study area, with almost no significant continuity.This condition may be due to the selective cutting system applied in the study area.Besides the randomly distribution of standing stock (trees) as the main cause, the second dominating cause of this abrupt spatial variations were due to human intervention such as shifting cultivation, road construction for forest harvesting (skidding and logging roads), arson, slash and burn, land encroachment, and natural causes such as wind throw, insects outbreaks, diseases and forest fire.For shorter plot distance, physiographic and topographic agents were the major causes.In fact, the site was quite heterogeneous in nature.Stand with high heterogenity and weak spatial structure caused anisotropy, which is expected in such naturally heterogeneous area, vanished and could not be seen in experimental variograms.
Results of Kriging showed that due to smoothing effect of Kriging, the variance of estimated values was much lower than the measured data, while the estimated mean was close to the data mean.However, the cross-validation results indicated a significant bias existed between kriged and measured data.Therefore, the biomass (even their original data such as tree density and basal area) was not a good candidate for Kriging.Our study results were in line with the result of Akhawan et al. (2010).
Lamandau study site Similar pattern was also provided during interpolating biomass distribution at the study site, where IDW method with power value p = 2 showed much better estimation than those obtained using Kriging method.This suggested that forest cover and stocks in logged over forest of Labanan were varied, with many gaps found due to logging activity, shifting cultivation and fire.This finding is in line with the research of Akhavan et al. (2010).This research finding also noted that the use of Kriging could not be recommended to estimate the spatial distribution of vegetation biomass when the stand condition varies widely.
As depicted in Table 3, the best interpolation method was IDW with power p = 0.5 which gave the mean deviation (MD) value of 17.48%.In this evaluation, the MD was also decreased as the power increased to p = 2.After p = 2, the MD increased gradually when the power p was increased.Similar patterns indicated that MD of IDW were much better than those obtained from Kriging. Figure 6 shows the trends of MD, RMSE, BIAS, Chi-sq, and AD at several values of p and Kriging models.The IDW method with p = 2 was consistently provided better estimation than Kriging method.
In general, the p = 2 consistently provides very low score, having a total score of 1.3.This also concluded that the best performance was provided by IDW method with p = 2.The highest score was given by Gaussian model (GS) which suggested that model was not suitable for estimating biomass distribution.Graphically, the score for each validation is shown in Figure 7.
Examples of the spatial distribution maps provided by IDW and Kriging for Labanan and Lamandau sites were shown in Figure 8 and Figure 9 respectively.It is clearly shown that most gaps existed in the secondary forest or disturbed forest could be accurately estimated by using IDW method.Clustered island-like gap (dark red) demonstrated in Figure 8a and Figure 9a would be more suitable using the IDW approach.The biomass accumulation or biomass scarcity that hampered in more massive area would be better suited using Kriging method (Figure 8b and Figure 9b).
As shown in Figures 8 and Figure 9, the spatial distribution of dry forest wass strongly influenced by forest disturbances.The ground investigations and satellite imagery interpretation showed that secondary dry land forest in either Lamandau or Labanan sites had been under heavy pressures from illegal logging, encroachment, illegal mining, and forest fires.The small openings caused by these disturbances led to the formation of many small gaps.Furthermore, these small gaps tent to lead to a reduction or even loss of autocorrelation with increasing distance between the plots (more than 1,000 m).This provided the rationale as to why the Kriging interpolation method did not provide good biomass estimation.In many forest regions in Kalimantan, it is well known by foresters that the timber standing stock variation was quite high, with coeficient of variance (CV) value of approximately 65%.This means that there were wide ranges of timber or biomass volume.This condition might also reduce the successful use of Kriging interpolation method.As noted by Wutzler et al. (2011), at stand scale, uncertainties of tree biomass carbon stocks in temperate forest were also high, having coefficient of variation ranging from 30 to 50%.The timber volume or biomass derived from tree density or basal area would not behave as a regionalized variable in the forest.Consequently, spatial distribution of the biomass was not auto-correlated over distance.As a result of these discontinuities, Kriging might not be a suitable alternative for estimating biomass of natural forest stock.Therefore, the best estimator would be using deterministic (the simple mean) such as spline and IDW.This contradited to the results of Montes et al. (2005) who used ordinary Kriging to estimate the cork oak production in Spain.
In this study we found that human interventions produced abrupt changes in the forest, whereas geostatistical methods were best suited for data, in which the value of the measured data attributes changed slowly.Based on the study results, we proposed to apply geostatistical approach for intact natural forests, plantation forest such as forest reserves where lack of abrupt changes available.This research was also intended to obtain insights of natural processes which weres necessary for close to nature forestry or in plantation forests.

Conclusion
From the foregoing discussion, it concluded that in logged-over dry land forest such as in Lamandau and Labanan sites, where forest received much human disturbances such as illegal logging, encroachment, illegal mining and forest fires, the deterministic spatialinterpolation method, i.e. the IDW, consistently provided better estimations than the geostatistical spatialinterpolation, i.e. the Kriging.Under such forest conditions, it could be concluded that the shorter sample-plot distance, namely 1,345 × 1,345 m in Lamandau and 1,450 × 1,450 m in Labanan, did not present better estimation results.For Labanan site, the IDW method with p = 2 gave estimation of 15% mean deviation, better than the Kriging method with 21% mean deviation; whereas for Lamandau site, an IDW with p = 2 generated a 17% mean deviation smaller than the Kriging with 27% mean deviation.
Figure 1 (a) Spatial distribution of IHMB data of Labanan site used for model development (green circle) and for verification (red dot); (b) distance between plots.(a) (b) Allconverting the AGB derived from IHMB-data (AGB ) IHMB using the equation developed byJaya et al. (2012), as shown in Equation [2].AGB = 7.355 + 1.003*AGB [2] unknown weight for the measured value at the ith i location of the point of interest d = the distance between x and x , i 0 i p = a power parameter, n represents the number of sampled points used for the estimation, x = the prediction location, 0 x = measurement at the ith location. i Figure 2 (a) Spatial distribution of IHMB data of Lamandau site used for model development (green circle) and for verification (red dot); (b) distance between plots.(a) (b) Model evaluationTo compare between the estimation results and field observation data collected separately from the model, model verifications were performed using 5 measures, i.e, Chi-square (Chi-sq), aggregative deviation (AD), bias error (BIAS), root mean squared error (RMSE) and the mean deviation (MD), with the formulas shown inEquation [5]  until Equation[12], respectively.
= the number of samples used Vti = estimated AGB value from the model Vai = actual value of AGB from the measurement The model is categorized as good when no significant difference between the estimated and actual value existed, 2 2

Figure 9
Figure 9 Spatial distribution of biomass stand at Lamandau site using (a) IDW with p =2 and (b) Ordinary Kriging with exponential model.

Table 1
Model verification for Labanan site, Berau, East Kalimantan •Chi-test result shows that all interpolation models are acceptable with no significant differences between the actual biomass and the predicted biomass (plot distance 1,000 × 1,050 m) •The mean deviation evaluation values indicate that the IDW method with several values of p provides better estimations than the Figure5Score performance of each interpolation method for Labanan site, Berau, East Kalimantan (the lower the better).

Table 2
Score for each interpolation method at Labanan site, Berau, East KalimantanThe study found that IDW with p = 2 provided the lowest score, indicating the lowest estimation error.Figure4Validation value of each interpolation method for Labanan site, Berau, East Kalimantan (the lower the better).