Estimating forest above-ground carbon using object-based analysis of very high spatial resolution satellite images

The potentials and limitations of very high spatial resolution images for aboveground carbon (AGC) estimation are unknown and the methods are not developed. This research was designed to develop a method that predicts AGC at the individual tree level using object-based analysis of very high resolution QuickBird satellite images and in situ diameter at breast height (DBH) measurements. This study was based on the fact that, crown projected areas (CPA) are strongly correlated to DBH. Assuming that CPAs are delineated with higher accuracy, a spatial model that predicts AGC can be developed using in situ DBH measurements and allometric equations. The DBH (1.3 m) of sample coniferous and broadleaf trees was measured and converted to above ground biomass and then to carbon. The panchromatic and pan-sharpened QuickBird satellite images were processed through object-based analysis to derive the CPA of coniferous and broadleaf trees and then followed by an accuracy assessment. The developed model predicted AGC stock and linearly explained about 58 and 55% of the variances for coniferous and broadleaf trees, respectively. Errors of CPAs resulted from over- and under-segmentation, sampling errors, and allometric errors as well as uncertainties in the developed model caused by structural errors.  
 
   
 
 Key words: Very high spatial resolution (VHSR), above-ground carbon (AGC), diameter at breast height (DBH), crown projected areas (CPA).


INTRODUCTION
Forests are major sources and sinks of CO 2 emissions and play a crucial role in the global carbon cycle. As global efforts are made to mitigate climate change through the reduction of emissions from deforestation and forest degradation, the need to accurately estimate the carbon stock of forest is growing. This demand is particularly associated with the introduction of carbon trading schemes in the Kyoto Protocol and through United Nations Framework Climate Change Convention (UNFCCC) commitments (Watson, 2010). Such E-mail: tenever4@gmail.com.
Author(s) agree that this article remain permanently open access under the terms of the Creative Commons Attribution License 4.0 International License assessments are useful in the creation of baseline data for measuring future changes of above-ground carbon (AGC) stock. Thus, forest AGC estimations are required at various spatial scales, with methods of quantification largely reliant on remote sensing techniques (Brown, 2001).
In combination with in situ biophysical sample tree measurements and information acquired from synthetic aperture radar (SAR), light detection and ranging (LIDAR) and optical and multi-sensor synergy measurements are widely used for above ground biomass (AGB) mapping (Goetz et al., 2009). AGB and carbon stock can be directly or indirectly estimated using remotely sensed data. Direct approaches extend satellite measurements into maps by calibrating them to field measurements of above ground biomass using a number of statistical techniques such as multiple regression analysis, knearest neighbour classification, and neural networks (Roy and Ravan, 1996;Nelson et al., 2000;Steininger, 2000;Foody et al., 2003;Zheng et al., 2004); indirect methods rely on estimated canopy structural parameters, such as crown diameter, which are first derived from remotely sensed data using multiple regression analysis or different canopy reflectance models (Wu and Strahler, 1994). In the category of medium spatial resolution images, Landsat time series images have been widely used for AGB estimation. K-nearest neighbour and nonlinear regression (Tomppo et al. 2002) and linear regression techniques (Hame et al. 1997) are used to estimate AGB based on spectral signatures, vegetation indices and image textures. However, the saturation of canopy reflectance and vegetation indices over time are problems associated with these techniques (Lu, 2006).
Advanced very high resolution radiometer (AVHRR) data is one of the most extensively used datasets for studies of vegetation dynamics on a continental scale. By using AVHRR data, AGB can be estimated based on the relationship between middle infrared reflectance and AGB (Boyd et al., 1999). However, the use of coarse spatial-resolution data for AGB estimation is limited because of the occurrence of mixed pixels as well as huge differences between the size of field measurement data and the size of pixels in the image, making the integration of sample data and remote sensing-derived variables more difficult (Lu, 2006). Hence, these methods are entrenched with considerable uncertainty and more robust methods are needed (Kohl et al., 2009). Nonetheless, the applicability of very high spatial resolution images for AGC estimation is unknown and methods are not yet developed, nor have their limitations and potentials been assessed.
With the advent of commercially very high spatial resolution images and developments in image analysis, conventional methods of pixel-based analysis are being replaced by object-based analysis. Such transformations in image analysis have created opportunities for detailed forest surveys such as the automation of individual tree crown outlines (Gougeon and Leckie, 2006), which are equivalent to crown projection areas (CPAs). Tree crown delineation algorithms such as the watershed transformation algorithm and the region growing and valley-following approaches are now widely applicable for the automation of tree CPAs. The success of these algorithms is largely dependent on the forest type, tree species and age composition, the image used and the image pre-processing such as the smoothing window size employed (Ke, 2008). For example, these algorithms delineate tree crowns of coniferous trees more accurately than broadleaf trees (Gougeon and Leckie, 2006;Ke, 2008). Delineating broadleaf trees is often problematic due to the presence of image noise, which could result in irregular local maxima within singing tree crown (Gougeon and Leckie, 2006). Image noise is often reduced by applying a certain smoothing window size (Derivaux et al., 2010). However, selecting the appropriate window size is difficult when the forest type is mixed in terms of both species and age since different trees require specific optimum smoothing window sizes (Wang, 2010;Ke, 2008;Li et al., 2008). Multiple local maxima limit the success of marker-controlled watershed transformation (Wang, 2010) and region growing, resulting in severe over-segmentation of single tree crowns and consequent commission errors. Overlapping tree crowns resulting in the clustering of tree crowns, in which one or more tree crowns are delineated as a single object causing omission errors, particularly limits the applicability of the valley-following approach (Gougeon and Leckie, 2006). Due to limitations of each algorithm under certain conditions, Ke (2008) indicated that developing a more robust algorithm with broad applicability would require taking advantage of the characteristics of multiple algorithms.
Studies have found a strong relationship between CPA and tree diameter at breast height (DBH) (Shimano, 1997). DBH is often used to estimate tree above ground biomass with the help of allometric equations (Muukkonen, 2007). This research was designed with the aim of developing a remote sensing-based method which first derives CPAs from very high spatial resolution QuickBird images through object-based analysis, then makes use of in situ field measurements of DBH and allometric biomass equations before finally predicting AGC at the individual tree level. However, before using CPAs as a proxy estimator of other parameters, a reasonable level of accuracy of tree crown delineation must be ensured.

The study area
The study area is 7 km from the city of Enschede in the Eastern part of the Netherlands. The forest lies between

METHODOLOGY
The study employed cloud-free QuickBird images, both panchromatic and multispectral bands, acquired in 2006 with less than 10° off-nadir angles (Table 1). These high spatial resolution QuickBird images can be used to retrieve forest canopy structure parameters, including CPA. However, the mutual shadows between tree canopies create difficulties in accurately retrieving these parameters. Off-nadir angles can also lead to geometric distortion of the CPA. Therefore, it is important to acknowledge that the image acquisition characteristics can affect the accuracy of automated CPAs. A QuickBird MSS image (2.4 m) was pan-sharpened using the panchromatic image (spatial resolution 0.61 m) to obtain a multispectral image with 0.61 m spatial resolution using IHS image fusion technique. As a result, trees with smaller tree crown sizes were made to be represented by many pixels, which were expected to improve the accuracy of their segmentation.

Sample data collection
Forest sampling, which addresses land cover variability results with greater precision than simple random sampling (Bertram et al., 2003), requires the creation of homogenous strata through image classification. The pan-sharpened QuickBird MSS image was classified into coniferous and broadleaf trees using a supervised classification technique with a maximum likelihood classifier. Therefore, 60 sample plots were distributed randomly in proportion to the areal extent of the coniferous and broadleaf forests ( Figure  2). A circular plot with a radius of 12.62 m (500 m 2 ) was chosen as a sampling unit. To establish this, a 12.62 m radius buffer was created from sample plot centres. The shape-file of sample plots was overlaid on the pan-sharpened image, and printouts of the sample plots were prepared for annotation and measurement of biophysical characteristics in the field. The images were saved with the sample plot shape-file in the Ipaq global positioning system (GPS) to facilitate navigation to sample plots. The GPS was mainly used to aid navigation to the neighbourhood of sample plots since it had errors reaching up to 10 m but in combination with the printout image, the sample plots were located on the ground. The DBH of all trees within the sample plot having DBH ≥ 10 cm were identified in the field corresponding to counterparts on the printout, measured using measuring tape, and recorded. Figure 2 shows forest types (pixel-based classification) and sample plot centres.

Allometric biomass regression analysis of sample trees
Allometric equations were used to estimate above ground biomass from the volumetric and structural dimension of the trees. Although DBH and tree height are usually used to compute the biomass of the trees (IPCC, 2003), tree height was not considered in the biomass equation since tree height was not measured in the field. Due to the scarcity of generalised allometric equations developed for Dutch and European forests, general allometric equations developed for Canadian hardwood and softwood forests were adopted for coniferous and broadleaf trees in the study area, respectively as most of the tree species in the study area are present in the Canadian forests. The allometric biomass equation takes the form of the power function as follows (Lambert et al., 2005).
Softwood AGB (kg) = Stem wood biomass + Branch biomass + Foliage biomass + Bark Biomass (5) Hardwood AGB (kg) = Stem wood biomass + Branch biomass + Foliage biomass + Bark biomass (10) Using these allometric biomass equations and sample field measurements of DBH, the dry above-ground biomass of respective sample trees were computed and converted to AGB and carbon, which constitute 50% of the dry AGB (Solicha, 2007).

Individual tree crown delineation
Tree crown delineations are one of the major factors affecting the accuracy of CPAs. In this study, a combination of valley following approaches and watershed transformations were used in ITC and eCognition software. Tree crown delineation in ITC software follows a valley-following approach performed in two major steps: shadow masking and rule-based tree crown delineation. These two processes are applicable to vegetated areas only, hence nonvegetated areas were masked prior to these processes. In the first step of the valley-following process, valleys of shade are created pixel-by-pixel by following small areas of shade, areas devoid of significant vegetation, and local minimas. A valley pixel is defined as a pixel continuing along a valley (minimum 8-connected pixels) with a radiance value lower than the pixels on either side when going in the valley direction (Gougeon and Leckie, 2006). As a result, good but often incomplete separations between tree crowns are obtained. This is followed by a rule-based crown delineation process attempting to finish the separation of tree crowns and produce tree crown outlines (Gougeon and Leckie, 2006). This result in a bitmap of objects referred to as isols (isolations) that typically represent individual tree crowns or, under certain circumstances and/or poorer spatial resolution, can be tree clusters (Gougeon and Leckie, 2006). To separate the clustered into individual tree crowns which result partly due to the spatial resolution of 60 cm and the presence of overlapping tree crowns, the clustered trees are transformed into a marker-free watershed segmentation after selecting them based on a threshold of a CPA greater than 100 pixels; equivalent to the maximum field-observed CPA. Watershed transformation considers the image to be processed as a topographic surface under three basic notions: local maxima, catchment basins and watershed lines (Chen et al., 2004). It can be performed with either marker controlled or marker free techniques. Markers are local maxima which usually lie at the centre of tree crowns (Wang et al., 2004). If a greyscale image is inverted, the local maxima becomes the local minima (Zhao and Popescu, 2007). The area between the local maxima and minima are the catchment basins that correspond to the tree crowns. Watershed lines are the local maxima of an inverted image. When the watershed is flooded from its minima continually, the water rises and the watershed lines form dams that prevent water entering from the adjacent catchment (Ke, 2008). These catchments correspond to the grey area of the tree crowns. Markers help avoid over-segmentation of a tree crown by watershed transformation (Zhao and Popescu, 2007).

Accuracy assessment of tree crown delineation
Following tree crown delineation, the accuracy of the resultant image objects must be assessed with respect to reference objects. The accuracy assessment technique must be appropriate for the specific purpose of tree crown delineation and must be free of bias. In this research, two accuracy assessment techniques were applied: area-based (Clinton et al., 2008) and count-based, which considers the degree of correspondence between reference objects and delineated objects (Ke, 2008). To avoid biases in the assessment of the accuracy of CPAs delineated through semiautomated techniques, clearly visible and well-distributed trees (70 coniferous and 30 broadleaf trees) were manually digitised from the pan-sharpened image prior to tree crown delineation. In the areabased accuracy assessment techniques, manually digitised CPAs were overlaid on the semi-automated CPAs in the JAVA environment using the JTS Topology Suite and GeoTools. The measure of accuracy was made in terms of the degree of oversegmentation, under-segmentation and the goodness of fit (Clinton et al., 2008).
In a supervised interpretation of the segmentation result, let X = {xi: i=1…n} be the set of n training objects, or assumed polygons, relative to which the segmentation is to be judged. Let Y = {yj: j=1…m} be the set of all segments in the segmentation. Let Ỹi be a subset of Y such that: For convenience, let area (xi ∩ yj) = the area of the geographic intersection of training object xi and segment yj. For each training object xi, the following subsets of Y exist: Workie 591 The union of these subsets is the subset Yi* = Ya ∪ Yb ∪ Yc ∪ Yd where Yi* is assumed to be the subset of segments that are relevant to training object xi. Processing over Y* is designed to minimize, if not eliminate, the effects of spurious intersections with very small parts of very large segments. Define # (Yi*) = pi and ∑ Thus, for each training object, there are pi segments deemed relevant. Define the following properties of the segments in Yi*: To supplement the area based accuracy assessment technique, count-based accuracy assessment techniques (with a degree of 1:1 correspondence between reference and delineated CPA) were also used. This accuracy assessment technique measures whether one reference tree corresponds to one or more trees when segmented. A single tree might be segmented into more than one tree due to the performance of the segmentation approach.

Object-based Isol classifications and accuracy assessment
Object-based classification separates objects created through image segmentation into similar categories based on their spectral, textural, shape and contextual properties (Mathieu and Aryal, 2005). Classification of different tree species was difficult due to the spectral similarities between many species in the area. As a result, all species were broadly categorised as either coniferous or broadleaf. First, the isols obtained from image segmentation were overlaid on the pan-sharpened QuickBird image and subjected to nearest neighbour classification. Feature space optimisation, which can give a better separation between classes, was used to select variables related to spectral, textural and geometrical characteristics. Feature space optimisation was made using variables such as the mean near infrared (NIR) and red band, maximum, minimum, and standard deviation values of the NIR band and area of the objects. Half of the field-collected data was used to train the nearest neighbour classifier (training data) while the remaining half was used for accuracy assessment data (ground truth).

Above ground carbon stock modeling and mapping
Regression modeling was used to quantify the relationship between dependent and independent variables. Regression equations are common statistical techniques used in the process of above ground biomass estimation (Lu, 2005). The relationship between semiautomated CPA and AGC stock were investigated through regression equations. To develop the model, sample field data (168 coniferous and 90 broadleaf trees) were divided systematically into two parts after identifying sample data, which are eligible for model development. Eligible sample CPAs have a 1:1 correspondence with automated CPAs and those that are not misclassified. Therefore, the first half was used as training data for developing the regression model. This was followed by model validation using the remaining data. This process aims to systematically establish a level of confidence of the model (Buranathiti et al., 2006). For validating the developed model, mean square error (MSE) techniques are applied (Oredein et al., 2011). MSE is a frequently used measure of the difference between values predicted by a model and the values actually observed. The developed model predicts total AGC using CPA as independent variable. Using the coefficients of the developed model, the CPAs of the validation data was used as independent variable and the model was run to generate the predicted total AGC (Ymodel

Above ground carbon stock of individual sample trees
The AGC stock which is about 50% of the dry AGB of the different components of a tree (steam, foliage, bark and branch) takes the form of power sigmoid. The rate and amount of AGC stock increases with the increment of DBH (Figure 3).

Tree crown delineation, accuracy assessment and classification
The success of tree crown delineation in ITC software depends on the existence of shade in between tree crowns. Due to the high basal density in some parts of the forest and the existence of trees with overlapping crowns, the isol delineation gave rise to clustered tree crowns (Figure 4b). These trees were separated from neighbouring tree crowns by transforming them into a marker-free watershed segmentation after selecting clustered trees based on a threshold of a CPA greater than 100 pixels equivalent to the maximum field-observed CPA ( Figure 5). Unlike marker-controlled transformation, marker-free watershed transformation does not need a marker; however, it may result in severe oversegmentation unless the image is well smoothed (Vincent and Soille, 1991). The problem of over-segmentation was disregarded due to the fact that high basal density and overlapping tree crowns resulted in a canopy that appeared relatively smooth. Figure 4 indicate the major steps and intermediary results involved before generating in CPAs of trees. Tree crown delineation was generally more successful for coniferous trees than broadleaf trees both using the valley-following approach and after combining that with MFWT ( Table 2). As with area-based accuracy assessment techniques, coniferous trees had a higher percentage of 1:1 correspondence than broadleaf trees but with a lower percentage of accuracy. Broadleaf trees severely over-segment to the extent that a single tree crown is sometimes segmented into more than four objects (Table 3).
The classification of coniferous trees was more reliable than broadleaf trees; hence, the likelihood for the user to find the same coniferous tree on the ground is higher ( Table 4). The producer's accuracy was slightly higher for broadleaf trees, which indicates that representation of broadleaf trees on the map was mostly accurate with minimal omission error. In general, classification errors arise due to: (i) the spectral similarity of species, (ii) effects of shadows due to low sun elevation angle (Table  1) and variability in trees age composition, and (iii) errors in tree crown delineation. Figure 6 indicates the CPA classification results.

Above ground carbon stock modeling and mapping
Linear models predicted AGC more accurately than nonlinear ones for both broadleaf and coniferous trees. The models explained more than 50% of the variance in AGC stock. In line with the accuracy level of the CPAs, the model for coniferous trees has lower residual variance than the model for deciduous trees which is about 42 and 45%, respectively.
Both models indicate that the carbon stock continues to increase linearly with the increase of the CPA. However, Shimano (1997) indicated the best fit curve between DBH (proxy indicator of carbon stock) and CPA relationships to be a power sigmoid for coniferous and deciduous trees. The power sigmoid curve shows that the rate of CPA increment decreases over time due to the competition of nearby trees for sunlight, nutrients and space. The relationship between CPA and the competition between trees for sunlight, nutrient and space also appears to be dependent on the DBH ranges. A study by Hemery et al. (2005)  crown diameter. Strong linear relationships exist for trees with a DBH range of 20 to 50 cm in particular, and slight reduction in the rate of CPA growth afterwards appears due to the effect of senility. Moreover, basal density andforest management appears to have its influence on the CPA-DBH relationship. In less dense forests, the CPA can continue to grow without the effect of competition for space from the nearby trees. Thinning and selective harvesting, which are common management practices in the study area, can also result in distance between trees and lead to a linear relationship. Hence, the linear relationship obtained in this research can partly indicate the minimal and insignificant effect of competition on the CPA of trees. Generally, despite minimal residual variance for coniferous trees in comparison to broadleaf trees, a considerable amount of the variance of AGC remains unexplained by the two models (Figure 7). The model validation indicates that on average the predicted total AGB deviates from the observed total AGC by about 173 and 392 kg of AGC for coniferous and deciduous trees in each observation, respectively (Table  5).
AGC was mapped (Figure 8) using the linear regression equation developed for coniferous and broadleaf trees (Figure 7). In general, broadleaf trees had a mean carbon stock higher than coniferous trees which is about 390 and 170 kg, respectively. This is partly due to the larger relative size of broadleaf trees. Of the total carbon stock in the forest, the highest proportion was found in the coniferous trees. Moreover, the area of private forest, which constituted the highest number of         (Nabuurs et al., 2000) of about 59 Mg C/ha (Nabuurs and Mohren, 1993). Moreover, remote sensing-based estimation of the carbon stock in woody biomass ranges from 25 to 60 Mg C/ha for different countries, the average being 42.91 Mg C/ha (USDA, 2003). However, remote sensing estimates, which are often made from very coarse satellites, have larger uncertainties particularly due to the occurrence of mixed pixels and considerable differences between the size of field-measurement data and pixel size (Lu, 2006).

Uncertainties and sources of errors in CPA-based AGC modelling
Object-based tree crown delineation using very high spatial resolution QuickBird images enables us to obtain CPA and estimate AGC at the individual tree level. The developed models were validated and used to estimate AGC linearly. Root mean square errors in the models were in the order of 0.173 and 0.392 Mg C per tree for coniferous and broadleaf trees, respectively. The carbon estimation for coniferous and broadleaf trees laid ± 0.17 and ± 0.38 Mg C with a 95% confidence interval, respectively. The uncertainty in coniferous trees is relatively small in comparison to broadleaf trees. This is in line with the accuracy level of tree crown delineation, partly revealing the interdependency between the accuracy of CPAs and uncertainties in model prediction. General sources of error may include sampling errors, systematic and random errors during DBH measurements (Zhang et al., 2010), and errors in the allometric biomass equations and tree crown delineation. Each tree species has a specific allometric pattern. However, due to the difficulty of classifying the CPA for each tree species, general allometric equations were enforced to a composite of different coniferous and sample Scots Pines (P. sylvestris) were involved in the model development, which affects the accuracy of the model. Accurate tree crown delineation is also a key factor for estimating other variables based on CPAs derived from image segmentation (Ke, 2008) as a slight error in over segmentation and under segmentation can affect the relationship between AGC and CPA. As indicated in Table 2, the goodness of fit (D) reveals about 20 and 34% mismatch between delineated and reference crowns due to over-and under-segmentation for coniferous and broadleaf trees, respectively.

SUMMARY, CONCLUSION AND RECOMMENDATIONS
CPAs automated from very high spatial resolution QuickBird images showed the potential for estimating AGC at the individual tree level despite some prevailing uncertainties arising due to sampling errors, errors in tree crown delineation and in the general allometric biomass equations, and the propagation of these errors through the model. Increasing the sample size is one method of reducing sampling errors. Taking more representative samples through the creation of more detailed strata by classifying all the tree species would be advantageous. However, this may require intensive field sampling and images with more spectral bands such as WorldView-2 satellite images. Through object-based analysis of very high spatial resolution QuickBird images, it was possible to automate CPAs with considerable accuracy. The level of accuracy was a function of the tree crown delineation algorithm, the tree type, and the basal density. The success of the algorithms in turn is dependent on the tree type, image pre-processing (smoothing filter window size) and the basal density. In mixed-species forests with a large variability in species, tree crown size and basal density, applying a single algorithm and image preprocessing technique proved less successful. Thus, a method which addresses individual tree requirements for an appropriate algorithm and includes pre-processing that takes advantage of different tree crown delineation algorithms is fundamental. Errors introduced through general allometric biomass equations mainly arise from structural errors in the model and when equations are forced to fit tree species growing in a different biogeography if not entirely different tree species. As these errors affect the overall accuracy, applying speciesspecific allometric equations and classifying all tree species is crucial. Quantifying the error propagation could also improve the model, but requires further research.