Analysis of a lidar voxel-derived vertical profile at the plot and individual tree scales for the estimation of forest canopy layer characteristics

ABSTRACT The goal of the current study was to develop methods of estimating the height of vertical components within plantation coniferous forest using airborne discrete multiple return lidar. In the summer of 2008, airborne lidar and field data were acquired for Loblolly pine forest locations in North Carolina and Virginia, USA, which comprised a variety of stand conditions (e.g. stand age, nutrient regime, and stem density). The methods here implement both field plot-scale analysis and an automated approach for the delineation of individual tree crown (ITC) locations and horizontal extents through a marker-based region growing process applied to a lidar derived canopy height model. The estimation of vertical features was accomplished through aggregating lidar return height measurements into vertical height bins, of a given horizontal extent (plot or ITC), creating a vertical ‘stack’ of bins describing the frequency of returns by height. Once height bins were created the resulting vertical distributions were smoothed with a regression curve-line function and canopy layers were identified through the detection of local maxima and minima. Estimates from Lorey’s mean canopy height was estimated from plot-level curve-fitting with an overall accuracy of 5.9% coefficient of variation (CV) and the coefficient of determination (R2) value of 0.93. Estimates of height to the living canopy produced an overall R2 value of 0.91 (11.0% CV). The presence of vertical features within the sub-canopy component of the fitted vertical function also corresponded to areas of known understory presence and absence. Estimates from ITC data were averaged to the plot level. Estimates of field Lorey’s mean canopy top height from average ITC data produced an R2 value of 0.96 (7.9% CV). Average ITC estimates of height to the living canopy produced the closest correspondence to the field data, producing an R2 value of 0.97 (6.2% CV). These results were similar to estimates produced by a statistical regression method, where R2 values were 0.99 (2.4% CV) and 0.98 (4.9% CV) for plot average top canopy height and height to the living canopy, respectively. These results indicate that the characteristics of the dominant canopy can be estimated accurately using airborne lidar without the development of regression models, in a variety of intensively managed coniferous stand conditions.


Introduction
A forest stand can be described in terms of its structural, compositional and functional properties, which are determined by site factors and the current forestry management objective. A great deal of spatial data is required for these tasks, such as estimates of stem density, canopy height and canopy base height (Shugart, Saatchi, and Hall 2010). The assessment of forest vertical structures is of importance in assessing regeneration success, competition, and biodiversity. Forest structure is determined by several factors, including species composition and the three-dimensional (3D) distribution of leaves/ needles and woody biomass. Many anthropogenic and non-anthropogenic processes alter forest structure including forest management and natural disturbances. The dimensions of the forest canopy is related to tree vigour and therefore to the timing and response to thinning (Smith 1986). Individual crown dimensions are correlated with stem size, and therefore timber volume (Avery and Burkhart 1994). This information can be of use to habitat modelling applications (e.g. Humphrey et al. 1999). Fire risk and damages due to wind and pests are dependent to a degree on the distribution of aboveground vegetation, formed by the collection of plant crowns comprising the canopy. Stand vulnerability can depend on the presence/absence of different vertical layers, such as the dominant, co-dominant, sub-dominant, and understorey (Sandberg, Ottmar, and Cushon 2001).
The measurement of detailed forest stand structural parameters is typically accomplished through the use of manual field-based assessment. Examples of such measures include tree size and height measurements, in addition to spatial distribution. This approach is often constrained by site accessibility, the availability of objective and effective measurement techniques, and the cost of personnel and equipment. These restraints may preclude the establishment of a sufficient number of field plots to sample the diversity of forest structural parameters throughout the environment. Such an approach is generally considered impractical for anything other than local-scale studies (Aplin 2005;Vierling et al. 2008). The application of remote-sensing technologies is often used as a means of extrapolating plot-based data to the landscape scale (Cohen and Spies 1992). The use of passive optical sensors is limited to recording horizontal landcover information observable from above, and thus do not provide direct information on the vertical distribution of canopy elements, particularly in uneven aged mixed species forest or where multiple canopy layers exist.
Airborne lidar remote sensing can characterize the 3D organization of vegetation (or physical structure) within forested environments that are of interest to forest managers. The use of lidar has rapidly come into prominence in estimating forest structural characteristics, such as canopy height, basal area, timber volume, biomass ), and leaf area index (Sumnall et al. 2016). This is frequently accomplished through the use of regression methods using lidar-derived metrics related to the vertical distribution and intensity of lidar returns over a range of scales and locations (Anderson et al. 2008;Lefsky et al. 2002;Lim et al. 2003aLim et al. , 2003bMaltamo et al. 2005;Naesset 2002). The derivation of individual tree metrics is possible using small-footprint laser data with a sufficient point density (Kaartinen et al. 2012).
High correlations have been found to exist between lidar data and some fieldworkderived forest stand structural attributes. However, there are a number of potential issues that limit the transferability of such models to other sites, and that relate to the range of forest types previously studied in order to derive relationships, with many developed for a specific region and limited number of species types. In addition, considerations over differences in the acquisition parameters of the lidar sensors between studies pose problems (Hopkinson 2006;Naesset 2009), and the time difference between field and remote-sensing data capture, for example models calibrated using leaf-on field and lidar data, can produce erroneous results when applied to leaf-off lidar data (Villikka, Packalén, and Maltamo 2012). The usage of lidar intensity remains a contested issue due to the proprietary methods that commercial sensors use to report return intensity, which can change in flight, making it impossible to directly compare two returns (Lim et al. 2003a).
Given airborne lidar's ability to partially penetrate through the forest canopy to the lower vegetated elements and ground, there is the potential to generate discrete statistics by classifying the lidar data point-cloud into vertical layers to determine statistics for the dominant canopy and understory layers using a process similar to the vertically recorded waveforms of large-footprint lidar systems (e.g. SLICER) (Lefsky et al. 1999). A number of studies utilize an approach where lidar returns are binned at set heights, also known as volumetric pixels (voxels). Each voxel requires the horizontal and vertical extent to be defined by the user, for example 1 m × 1 m × 1 m, allowing vertical stacks of bins to be generated for each 1 m × 1 m of horizontal extent. Each voxel element will contain a single value for the lidar returns located within the set threedimensional extent, such as the sum of returns, or the average intensity value of those returns (e.g. Chasmer, Hopkinson, and Treitz 2004;Lee, Lucas, and Brack 2004;Popescu and Zhao 2008;Wang, Weinacker, and Koch 2008). The resulting vertical distribution of the sequential voxel values within that column is referred to as the vertical profile.
A curve function can be fit to the values of each of the vertically sequential voxels in a vertical column (e.g. the frequency of returns by height), which allows vertical vegetation layers to be characterized and undesirable features (e.g. noise) to be suppressed. For example, Muss, Mladenoff, and Townsend (2011) apply a spline function to the vertical profile to smooth the distribution and detect layers and compute metrics based on these distributions for the estimation of forestry metrics (e.g. mean stem diameter). Dean et al. (2009) fitted truncated Weibull functions to vertical profiles, of return frequencies, to identify the dominant canopy layer and estimate canopy crown base and median heights, which were then input into regression models to predict average stem diameter. Unfortunately, at this stage there is little consensus on what the optimal curve-fitting procedure is.
The highest proportion of lidar returns intercept with the top of the forest canopy, with a reducing number of returns from the lower strata as permeability decreases (Chevalier, Steinvall, and Larsson 2007). Therefore, the use of smaller sample extents, or lower pulse densities, increases the potential for less-accurate estimates of forest vertical layer heights and sub-canopy features going undetected (e.g. Chasmer, Hopkinson, and Treitz 2004;Dean et al. 2009). Lindberg et al. (2012) presented an approach to compensate for this effect based upon the Beer-Lambert law for small-footprint waveform lidar, which provided somewhat more accurate estimates of vegetation volume when compared with discrete return data when analysing the vertical profile.
A number of studies have been implemented to estimate the presence of understory features, for example Hill and Broughton (2009) and Martinuzzi et al. (2009) demonstrate the potential to map the presence and absence of such features. Other studies go beyond this to estimate sub-canopy layer heights, for example Wang, Weinacker, and Koch (2008) and Jaskierniak et al. (2011), which implemented a plot-level height-bin analysis, where Gaussian curve-fitting or mixture modelling, respectively, were used to determine dominant and sub-dominant forest layers characteristics.
A number of methods have been developed to automatically identify tree crowns automatically from airborne lidar data (e.g. Kaartinen et al. 2012). The variability of canopy 3D characteristics can be assessed in a neighbourhood. For example, Holmgren and Persson (2004) utilized a regression approach to estimate crown base height and discriminate between broad species types. A vertical profile voxel-based approach was employed by Popescu and Zhao (2008). The authors initially delineated ITCs using a local-maxima-based detection approach applied to a canopy height model (CHM). Vertical profiles were created for each ITC object intersecting the lidar data, whereby a Fourier filter was applied to estimate crown base height.
The voxel-based approaches have commonly been applied to characterize the dominant canopy layers at plot-level scales. The separation of different vegetation layers holds the potential for mapping the vertical characteristics of vegetation, which would assist with assessments of productivity, fire fuel assessment, or evaluating competition. Therefore, the overall goal of this study was to develop an approach to analyse airborne lidar data to detect and characterize the individual vertical layers, with particular focus upon the dominant canopy, for multiple Loblolly pine dominated sites in the southeast of the USA at the field-plot and individual tree scale in stands of varying ages, structural types, and stem densities. More specific objectives were to: (1) Perform a regression analysis to statistically derive models for the estimation of average canopy top height and HTLC for comparison with layer finding methods; (2) Develop an automated method to identify the dominant canopy layer and estimate its height and depth from the vertical profile, and to determine the presence of sub-dominant layers; (3) Develop and apply an automated individual tree detection approach and then to test the automated canopy layer detection method at the individual tree-scale; (4) Compare plot-level area-based and averaged-ITC estimates for the height and crown depth of the dominant layer.

Study sites
Four study sites located within North Carolina and Virginia, USA, were used for the current research project. The sites were established and maintained in support of research studies investigating the role of intensive management in optimizing Loblolly pine (Pinus taeda L.) production. These studies were a joint effort between the Forest Productivity Cooperative (FPC) (http://forestproductivitycoop.net/), academic institutions, the USDA Forest service, the Virginia Department of Forestry and private industry. The United States Department of Agriculture (2008) defines typical Loblolly pine normal yield of natural, even-aged Loblolly pine stands on averages sites, such as those of height 27 m, ranges from approximately 133 m 3 ha −1 in trees 9 cm in diameter at breast height (DBH), to 29.4 m 3 ha −1 in trees 24 cm DBH at age 20 years, to 427 m 3 ha −1 DBH at age 60 years. Stands are commonly rotated at either 15 or 30 years. It should be stated, however, that many of the stands sampled in the current study utilize an experimental design and may differ from the above relationships.
The first of the study sites was the RW195501 trial (RW19), which is part of a study examining the effects of fertilization and thinning in mid-rotation stands. This trial is located in the Piedmont of Virginia in Appomattox County at 37°26ʹ 32ʺ N and 78°39ʹ 43ʺ W. A total of 32 plots were installed in a 13-year-old stand. The plots vary in size from approximately 400 to 1280 m 2 , and could vary in shape from square to rectangular. At the time of the lidar acquisition, no fertilizer had been added; each plot had been cleared of understory vegetation, however.
The second site in Virginia, the RW180601 (RW18), was part of study designed with the objective of understanding the optimal rates and frequencies of nutrient additions for rapid growth in young stands. The trial was located in the Piedmont site of Virginia, Brunswick County, at 36°40ʹ 51ʺ N and 77°59ʹ 13ʺ W. A total of 40 plots were installed in 1999 in a 6-year-old planted stand, and again could be of a square or rectangular shape. A total of 32 plots (out of 40) were subject to thinning. Four control plots were installed in thinned locations, which had no nutrient additions but were weed controlled. The remaining plots had complete weed control and nutrient treatments. Plot size varied from approximately 400 to 470 m 2 .
The third site was the Southeast Tree Research and Education Site (SETRES), and was located in sand hills of North Carolina, in Scotland County, at 34°54ʹ 17ʺ N and 79°29ʹ W. This trial was established in 1992 in an 8-year-old plantation. The goal of the study was to quantify the effects of nutrient and water availability on above-and belowground productivity and growth efficiency in Loblolly pine. Treatments consisted of four different nutrient additions types and irrigation (Albaugh et al. 1998) and weed control. Plot size was 900 m 2 (30 m × 30 m), four blocks and four plots per block, for a total of 16 plots. One block had no fertilizer treatment and was used as control.
The final site was the Henderson Long Term Site Productivity Study (Henderson) located in Vance County, NC, at 36°26ʹ 52ʺ N and 78°28ʹ 23ʺ W. It was established in 1982 for monitoring the effects of soil management practices on soil structure, organic matter, and nutrient contents and pine growth. Treatments consisted of two levels of biomass harvest, stem wood only or whole tree removals; two site preparation methods, chop and burn, or shear, pile, and disc. There are three groups, with eight plots per group (four pairs), totalling 24 plots in the study. Each pair had the same treatment type; however, vegetation control for the first five years was implemented in only one of the paired plots. Plot measurement size is 450 m 2 (15 × 30 m). For more information, see Vitousek and Matson (1985).

Field data
The field data utilized here was collected as part of Peduzzi et al. (2012). All study sites were measured during the 2008 dormant season and located using a handheld GPS. Total number of tree stems with their DBH greater than 10 cm (measured at 1.3 m above ground) was recorded within each plot in addition to DBH. Individual stem locations within the plot were not recorded. Total tree height and height to live crown (HTLC), measured up to the height of the lowest living branch, and was recorded for each tree using the Haglöf Vertex hypsometer for all stems within the plot extent for locations at Henderson and SETRES. Subsampling of stems was implemented in RW18 and RW19, where between 30 and 50 tree heights and HTLCs were recorded for each plot depending on plot size. And finally the presence or absence of understory vegetation was recorded through visual assessment. Lorey's mean height was calculated by multiplying the individual tree heights by their basal area, and then dividing the sum of this calculation by the total stand basal area. Plot-level summary values were then computed for HTLC, which consisted of the mean, standard deviation, and coefficient of variation (CV). The CV is a standardized measure of dispersion of a frequency distribution. Table 1 contains the summarized field measurements.

Lidar data
2.3.1. Pre-processing Small-footprint, discrete-return, lidar data were acquired for all study sites in late August 2008 to the same lidar acquisition specifications. The system was an Optech ATLM 3100. The system could record up to four returns per pulse, with an average sampling density of five pulses per square metre. The scan angle was below 20°. The instrument's vertical accuracy over bare ground was 15 cm, and horizontal accuracy was 0.5 m. The flying altitude was approximately 1200 m with a pulse repetition frequency of 70 kHz. The laser operated at 1064 nm with a beam divergence of approximately 0.3 mrad.
A number of pre-processing steps were required before metrics could be derived from the lidar data for analysis. All of these steps were performed using the RSC LAS Tools software (Armston 2014). The classification of ground returns utilized a progressive morphological filter as outlined in Zhang et al. (2003). Above-ground heights were calculated by subtracting the corresponding ground point heights, interpolated into a surface (nearest neighbour method), from the original (unclassified) data set. The aboveground height lidar metrics were then calculated. All lidar returns that were intersected within field plot horizontal extents were clipped from the data sets, and metrics were generated utilizing the R software (version 3.1.1.) (http://www.r-project.org/). Summary metrics were calculated for the heights of all returns, the heights of vegetation returns (all returns with height >0.2 m), the intensity of all first returns, and the intensity of vegetation first returns (returns with height >0.2 m). The summary metrics of mean, median, standard deviation, variance, coefficient of variation, skewness, kurtosis, maximum, and minimum and percentile heights, at 5% intervals, were computed.

Individual tree crown delineation
Individual trees were delineated using the lidar-derived CHM by an automated processing approach developed in the R software, utilizing an approach similar to that defined in Chen et al. (2006), a brief description of which follows.
A high-resolution CHM (0.5 m × 0.5 m) was created for each of the lidar data sets, and an initial circular weighted matrix focal (moving window) analysis was implemented, using a search radius of 1 m, in order to smooth the height features within the original layer. The resulting 'smoothed' raster layer provides a three-dimensional surface of the highest vegetation heights above ground using a linear interpolation method. A local maxima search, using a 3 × 3 cell search window, was applied to the CHM in order to identify local maxima, or peaks, within the surface, which should correspond to individual tree-top locations. A proximity function was then applied to identify the clustering of maxima points, and to remove or merge maxima points defined as too close to one another. For each maxima location, the distance was calculated to each of its immediate neighbours in a search window of radius 5 m. The proximity limit here was defined as less than 1 m in order to sample each of the eight neighbouring 0.5 m × 0.5 m cells in the raster layer. If one or more maxima locations was within range, the maxima IDs were added to a unique group. Groups were further merged if the IDs were duplicated. If these maxima locations within a group were not of equal height, then the tallest maxima was preserved, otherwise the most central point was kept.
Once maxima points had been identified, an iterative region growing approach was applied to the raster cells of the CHM in proximity to the maxima point location, whereby the nearby cells were searched (in eight directions from the focal point), and were added to the 'crown object' if the cells occupied a lower height value to that of the maximum initial height, and neighbouring cells classed as 'crown' in subsequent iterations. Region growing was prohibited to expand into a raster cell of 0 m (i.e. ground). Crown radius was estimated from the average distance from the maxima point in the four cardinal directions to the perimeter of the crown object (north, south, east, and west), and was represented as a generalized horizontal circle centred on the maxima location. Crown area was estimated by applying the formula for calculating the area of a circle, using the estimated crown radius. In order to illustrate the position and extent of the delineated ITC objects, a shapefile polygon of the estimated circular horizontal area was created for each. Polygon overlap is present for multiple tree crowns, due to the method of creating circular polygons, this means one ITC circle may include elements from one or more trees.

Creation of height bins and vertical profiles
The basic concept behind estimating the heights of vegetation layers from the frequency of returns within height bins in a vertical profile is to detect the position where the frequency of returns drops abruptly below larger number of returns. Such a technique is dependent on the 'shape' or distribution of return frequencies stratified by height. According to Magnussen and Boudewyn (1998), lidar returns within the tree crown mimic sampling with probability proportional to the projected tree crown sizes. Within their study, an estimated 50% of lidar returns hitting the canopy in a plot are presumably at or above a height that coincides with the height of mean leaf area index. What follows is a description of how the drop in return frequency of the vertical profile is detected.
A method was developed to automatically predict the vertical position of canopy layer features, a summary of which follows. All points that intersected with the field plot extents were clipped from the lidar data sets, and height bin metrics were generated, where for each 3D bin the number of lidar returns within were summed. Height bins were generated for each plot at 0.2 m intervals, which produced 200 layers (0-40 m) per plot. Figure 1 shows a cross-sectional illustration of the height bins generated at the plot level and how this can be visualized as the number of returns stratified by height. In a similar approach to Holmgren and Persson (2004), a noise filtering step was implemented where a value of zero was assigned to each bin, which contained less than 0.5% of the total returns.
An automated method of stratifying layer features from the vertical profile for each plot was developed using the R software. A Nadaraya-Watson kernel regression was fit to the distribution of lidar returns within the vertical profile above 0.2 m. A threshold of 0.2 m was selected so the high number ground returns did not influence the line fitting and the detection of features near the ground level. This approach was selected as it can be fit to multi-modal distributions rather than other smoothing approaches such as Gaussian and Weibull functions, which are restricted to unimodal distributions.
The kernel regression approach functions by forming a weighted average of all the y-values corresponding to points whose x-values are close to the x-value of a point being plotted. The function that defines the weights is known as a kernel, and the number of points involved in the weighted average is based on a parameter known as the bandwidth. Due to the variation in average tree canopy heights between plots, and the difference in the total number of bins containing data, a single bandwidth parameter for all sites proved to be insufficient. A subset of plot locations was selected, which covered a variety of stand heights and lidar subsets produced, where a curvefitting procedure was implemented using a range of bandwidth values in order to test their ability to delineate vertical layers. It was found that higher plot heights were better represented by larger bandwidth values, and too smaller value resulted in the detection of false positives. Therefore, an equation was developed to automatically select an appropriate value relative to the plot maximum height (as defined by the maximum height of lidar returns within the current horizontal extent) for the curve-fitting procedure. The formula is expressed as: where H plot is the maximum lidar above-ground height within the plot. Note that the equation was developed for the current field plots where the canopy height ranged from 12 to 27 m. Vertical layers were then identified through the detection of local maxima and minima from the fitted regression line. This is illustrated in Figure 2. If bins neighbouring a maxima or minima location were of the same value, the middle height value is reassigned to that of the maxima or minima. Locations along the fitted line function representing bins of zero returns were all classified as minima locations. Bins in which no returns were recorded were considered minima locations in the top half of the vertical profile only. The presence of no return bins in the lower half may be potentially caused by the lower proportion of returns beneath the forest canopy. The minima detected above and below a maxima formed the basis of extracting layer height information. Layer top height, bottom height, and maxima height were recorded for each detected layer. The dominant canopy layer was classified as a layer above one-third height of the maximum plot height with the largest maxima (i.e. the largest number of lidar returns at the layer 'middle' point). For all field plot locations this was always the highest layer detected. Therefore, for each field plot extent, lidar-derived estimates of average canopy to height and HTLC were calculated.

Height to the living canopy estimate correction
Due to the distribution of returns, a number of vertical profiles exhibited large separation between maxima and minima below causing an underestimation of height to the living canopy (e.g. as in Figure 2). This applied to the dominant canopy layer only in taller stands (above 15 m). To correct this, an additional processing step was applied when the separation between maxima and minima below it was beyond an arbitrary threshold of 6 m, but only when canopy top height was greater than 15 m. The former threshold, of maxima-minima separation, was determined through the trends observed in the field data values of crown length (HTLC subtracted from top height), where the average value varied from 6.19 m (Henderson) to 7.36 m (RW19). Thus, the largest separation between maxima (already beneath the top) and minima could not be more than 6 m. The current threshold should be applicable for other similar canopy architectures; however, as crown length increases, the threshold should be increased. If this criteria was met, a test was applied to the section of the vertical profile between the maxima and minima in question to assess the possibility of multiple layers existing within this vertical extent. A function was developed to locate the point at which the smoothed line function between the maxima and minima location was at maximum curvature, for example where there was a change in the rate at which the number of returns decreased with height (see Figure 3). This was accomplished by computing the difference between each bin frequency value and a straight line vector between the start and end points within this section of the smoothed regression line, then identifying the locations of the difference between largest values. A new layer was inferred if the largest difference was larger than 1.5 times the median difference value. If sufficient change was detected, the height value at this point was selected as the bottom height of the canopy layer above it. A new layer was then created between the new bottom height and the minima beneath it. It is not recommended to apply this method to all vertical profiles, as the shape of the curve and position of the minima below can vary between sample locations, and cause the above approach to overestimates layer heights. Therefore, each field plot extent could have an extra estimate of HTLC if the above criteria were met.

Individual tree crown vertical profiles
As noted previously, the circular ITC shapefiles may overlap one another to some extent. Therefore, larger ITC circles with a radius larger than 2 m were reduced to 66% of its original size, to increase the probability that the lidar returns sampled were from the correct tree. As with plot-level metrics, lidar returns intersecting with the horizontal extent of tree crowns were extracted and utilized to produce 1 m vertical height bins and create a vertical profile for each ITC object. The selection of 1 m height bins was to account for the smaller number of returns incident within each delineated ITC object. The maximum bin height considered for each ITC vertical profile was defined by the maxima height derived from the CHM. Again, Nadaraya-Watson kernel regression lines were fit to smooth the vertical profiles, and local maxima and minima defined. For the ITC profile only the dominant crown characteristics were estimated (see Figure 4). A number of the ITC vertical profiles exhibited multiple maxima within the tree crown, as in Figure 4. A simple function was created to define the crown vertical extent; this can be summarized as follows: (i) maxima were interrogated from highest to lowest; (ii) maxima must be located above one-third of the tree height (defined by CHM maxima); (iii) the selected maxima must be sequential; (iv) maxima cannot be separated by a minima with zero returns; (v) maxima must not be separated by more than 3 m.
If these criteria were met, the group of maxima were considered part of the living tree crown.
As with the plot scale analysis, layer top, middle (largest or highest maxima), and bottom statistics were calculated for each ITC object. No HTLC corrections (as in Section 2.3.4) were applied at the ITC level. Those ITC objects intersecting with the field plot locations were extracted. These intersecting ITC objects were then used to calculate averaged (mean) top height and HTLC values for each plot in order to make them comparable to those metrics calculated at the plot scale. In addition, the number of ITC objects at each plot level was extracted. Estimates of HTLC variability at the plot scale were also computed.

Statistical analysis
In order to assess if there were any relationships between the field measured average canopy top height and HTLC, and to compare the relative accuracy of each of the three estimation methods from the lidar data to the reference field measurements, three statistical measures were calculated. The first of these was to compare the correlation between the field plot data and predicted using univariate regression analysis, implemented in R. The coefficient of determination (R 2 ) was calculated for each comparison, and is considered to indicate the proportional amount of variation in the response variable explained by the independent variable (Field 2013). root mean square error (RMSE) was also calculated as a means of assessing performance: ( 2) where n is the number of observations, i denotes individual plots, y i is the field plot value, and ỹ i is the corresponding lidar-derived prediction. The coefficient of variation of the RMSE (CV RMSE ) was also calculated by where y mean is the mean of all field observations.

Relationships between field variables
The field-measured plot-averaged values for average canopy height and HTLC were input into a linear regression analysis. An R 2 value of 0.96 was produced, indicating a high level of correlation between the two field metrics. The relationship is expressed as HTLC ¼ 1:082 Â H canopy À 8:3839; where H canopy is the average canopy top height.

Plot-level-based regression analysis
The correlations of the regression analysis for the two plot-level summary metrics are presented in Table 2. Both the models developed for plot average top height and HTLC produced high correlations (R 2 > 0.9), were statistically significant according to Student t-test (p < 0.05), and only one covariate was used in the model. Overall RMSE was 0.42 m (2.44% CV RMSE ) for average canopy top height and 0.48 m (4.91% CV RMSE ) for average HTLC in a plot. A study site summary is presented in Table 3.

Plot-scale voxel estimates of canopy characteristics
The estimation of dominant canopy characteristics using the plot-level voxel-based analysis, average canopy top height (Lorey's mean height), average HTLC, and average HTLC with correction from the fitting of functions to the bin lidar returns is summarized in Table 4. Lorey's mean canopy height was estimated with an overall accuracy of 1.0 m RMSE (5.9% CV RMSE ) and overall R 2 value of 0.93. Broken down to the study site level, RMSE varies from 0.7 m to 1.5 m, whereas CV RMSE varies from 4% to 9%. The CV RMSE value corresponds to the field data summary presented in Table 1, where the smallest CV RMSE is for RW18, and the largest was the SETRES site.
Estimates of HTLC, without the correction, produced an overall R 2 value of 0.78 and an RMSE of 1.6 m (17% CV RMSE ). Whereas estimates of HTLC with the correction produced an overall R 2 value of 0.91, and RMSE of 1.1 m (11.0% CV RMSE ). The correction was applied for 73 out of a total of 112 plots. Uncorrected HTLC error was larger than corrected HTLC for all four sites, and the corresponding values of the RMSE were between 0.5 and 1.3 m larger. The uncorrected estimates of HTLC of RW18 have a larger RMSE, but a smaller CV RMSE value than for the corrected, however. The 11 plots (9 were thinned plots) for which the correction was applied, resulted in HLTC overestimates. It should be noted that all of these 11 plots were located at the overlap between two flight line swaths. The RW19 site exhibited the lowest RMSE for estimates of HTLC, where the correction was applied for all but one of the field plot extents. The plot-level results using the HTLC with the correction are presented in Figure 5a. The RMSE values for each study site are also similar to the average field data plot standard deviations for each study site (Table 1) for the HTLC estimates with the correction.  Understory vegetation was removed in all the field plots assessed as part of this project, apart from those located at Henderson. The field plot design at Henderson incorporated four different initial site treatment types and understory control, resulting in clusters of eight plots. Figures 6 and 7 illustrate the vertical profile generation for one such cluster, allowing the user to compare the vertical distribution of canopy elements, such as differences in dominant canopy height and the detection of sub-dominant vegetation layers and the presence of understory. From examination of the vertical profile, where there is a presence or absence of understory layers is clear, corresponding to no vegetation control and vegetation control site management types, respectively. As illustrated in Figures 6 and 7, the vertical position of the dominant canopy can alter between sites, where those subject to herbicide use were typically taller.

Individual tree crown metrics
There were a number of difficulties in validating the results of the ITC delineation approach that relate to determining an objective correspondence between fieldrecorded and lidar-determined positions, as stem positions were not recorded in the field. The current study therefore attempts to estimate individual tree crown parameters and report results by aggregating estimates to the plot level.
The number of stems recorded in the field for each individual plot was compared with the number of ITC objects recorded within the field plot horizontal extent. The relationship produced an R 2 value of 0.85 and an RMSE of 23.9 (39% CV RMSE ); however, as seen in Figure 8(a), as the stem populations in each of the field plots increased, the number of ITC objects created was not as many, resulting in an underestimation. A number of plots within the Henderson site were the exception to this. The regression equation (with a zero intercept) indicated there were approximately 1.35 times more stems recorded in the field data. The largest RMSE for the number of ITC objects relative to the number of stems recorded in the plot was for the RW19 and SETRES sites, i.e. where higher RMSE values are present in locations where plot sizes were larger, and contained more stems, for example RW19 encompassed plots with a minimum of 41 to a maximum of 184 stems. The CV RMSE value is smallest for the Henderson site (18%). Whilst the other sites have higher CV RMSE values (>30%), the largest is for RW19, which may be due to variable plot size.
In order to test if the ITC delineation was detecting the co-dominant and dominant trees only, the number of trees in the sub-dominant category were removed from the field plot totals. Percentage quantile heights were calculated, and stems that were below 30% were removed. For the subsampled RW18 and RW19, the number of stems was considered representative and the relative proportion subtracted from the total. The new field plot totals were then compared with the number of intersecting ITC objects and produced an R 2 value of 0.83, with an RMSE of 12.4 (25% CV RMSE ) (see Figure 8(b)). The regression equation (with zero intercept) produced a coefficient value of 1.056. As Figure 6. Vertical profiles created for four of eight field plots in one cluster located at the Henderson site, with initial treatment type of stem removal. In this example it is possible to compare the distribution of vertical canopy components (dashed-red lines), such as differences in dominant canopy height, in addition to detecting sub-dominant vegetation layers, and understory presence (dashed-blue lines). The green triangle represents the height of a vegetation layer intersecting the ground. (SPD, shear pile and disc; CD, chopped and burned). before, the plots located at Henderson created more ITC objects, whilst the overall correspondence was improved.
Plot-level summaries for the estimates of top heights and HTLC for ITC objects intersecting with the plot horizontal extent produced slightly higher R 2 values in comparison with the plot-level layer fitting equivalent. Overall, the R 2 values for the mean ITC top height and field-calculated mean canopy height was 0.96, RMSE was 1.4 m (7.9% CV RMSE ). The overall estimates of ITC mean HTLC produced an R 2 value of 0.97, with a lower overall RMSE of 0.61 m (6.2% CV RMSE ). Table 5 shows the results calculated for each study site. The correspondence between field-recorded aggregated plot-level values of HTLC is presented in Figure 5(b).
ITC plot summary metrics for average canopy top height were slightly higher than those for the plot-level layer fitting approach (0.3-0.8 m). The RMSE for the SETRES site was lower, however (by 0.8 m). The SETRES site has the highest CV of plot heights, according to Table 1. Plot-averaged ITC HTLC RMSE values were lower than those produced from the plot-level layer corrected fitting approach (see Table 4); the differences ranged from 0.9 m to 1.9 m, whilst the estimates for RW19 were similar. There was multiple ITC-level estimates of top height and HTLC for each plot, and there was some within-plot variability in the field-recorded values. Figures 9 and 10 illustrate the variability of HTLC within individual field plots within the context of the RW18 and SETRES sites, respectively, for field data (a) and lidar-derived ITC estimates (b). Note that there are similarities between the inter-plot variability, such as between thinned and unthinned plots, grey and white, respectively in the case of the RW18 site. RW18 ITC HTLC estimates typically have a larger range of values, where the medians are often skewed when comparing them to the field equivalent. The largest variability is typically encountered for those plots that have been thinned (grey). In order to exemplify the causes of this issue, two plots will be explored in greater detail. Plot number 7 for example has a large range of ITC HTLC values (4 m to 11 m), in contrast to field measures (~8 m to 10 m). This indicates that non-canopy returns have been included or canopy returns were missing, resulting in an inaccurate representation of HTLC range. A number of plots within this example exhibit crown centres that are very close together  and radii larger than 2 m; for example, plot number 30 has eight ITC centre points less than 2 m apart. Therefore, returns from neighbouring overlapping tree crowns may cause erroneous estimates. The large variations encountered for RW18 are not observed to the same extent in the SETRES comparison, indeed which corresponds closely to the patterns observed in the field data. Overall, the ITC values of HTLC have larger ranges than their field-measured counterparts.

Discussion
The current study has presented one statistical and two methods of fitting a curve to lidar-derived voxel data at both the plot and ITC scales for estimating canopy attributes Figure 9. The plot-level variability of HTLC for each field plot is presented for field data (a) and lidarderived ITC estimate data (b) for all plot locations at the RW18 site. The boxplots are colour coded to correspond to the broad management type, where red represents thinned control, grey is thinned and fertilized, and white is unthinned plots. Figure 10. The plot-level variability of HTLC for each field plot is presented for field data (a) and lidar-derived ITC estimate data (b) for all plot locations at the SETRES site. The boxplots are colour coded to correspond to the broad management type, where grey represents fertilized and irrigated, and red represents control plots (irrigated).
in Loblolly pine plantation forests, in addition to demonstrating the potential ability to detect the presence or absence of understory layers within this context.

Plot-level canopy vertical layer analysis
The results presented here for plot-level lidar estimates of average canopy top height (Lorey's mean height) from the plot-level statistical, plot-level voxel-based curve-fitting and averaged ITC voxel-based curve-fitting analysisall correlated well with fieldmeasured values with R 2 values of 0.99, 0.93, and 0.96, respectively. All three methods produced similar RMSE values of 0.4, 1.0, and 1.4 m, respectively. Canopy height estimation accuracy in coniferous stands are consistent with those in other studies, such as in Gaveau and Hill (2003) where estimations of tree canopy height differed by an average of ±2.12 m when compared to CHM samples points for forest locations with canopy height ranging between 10 and 18 m. Naesset (2002) reports an accuracy of 0.61-1.17 m (R 2 = 0.80-0.95), and for Andersen, Mcgaughey, and Reutebuch (2005) 1.3-1.5 m (R 2 = 0.98) for statistically estimated average plot canopy height. Other studies that employ methods of layer detection, for example Muss, Mladenoff, and Townsend (2011), used such metrics into statistical models and estimated canopy top height to ±2.09 m (R 2 = 0.82), whereas Morsdorf et al. (2010) used 2D cluster analysis to estimate canopy top height with a bias of 1.1 m (±0.5 standard deviation) in pine stands with maximum canopy heights between 14.8 and 17.8 m.
Likewise, plot-level lidar-derived estimates of HTLC corresponded well with field measurements. This section will discuss the statistical and plot-level voxel-based curve-fitting approaches. The estimates derived from ITC analysis will be reviewed in Section 4.2. Statistical estimates produced a model with an R 2 value of 0.98 (RMSE 0.48 m), which used one predictor variable, the 90th height percentile. This was very similar to the regression model produced for average canopy top height ( Table 2). The analysis of the two field-recorded metrics with one another indicated they were highly correlated (R 2 = 0.96), it is therefore likely that the lidar predictive model parallels this relationship. These characteristics are likely typical for the mid-rotation closed-crown stands used for this study. Its transferability to other sites and to data sets using different acquisition parameters is likely an issue (Naesset 2009).
The two estimates of plot-averaged HTLC derived from the plot-scale voxel-based analysis provided a differing correspondence to field data, where estimates not including the correction produced an overall smaller R 2 value (0.78) when compared to estimates including the correction (0.91), and reduced the largest RMSE by 1 m in the case of the RW19 site. The RMSE of estimated values also closely corresponds to the standard deviation of values recorded in the field data. Within the context of the current study, uncorrected estimates of HTLC produced lower values than those recorded in the field. The segregation of the vertical profile into layers is dependent on the fitted curve function and the presence or absence of returns within height strata. Explanation for the phenomena of inaccurate HTLC may be caused where vegetation layers overlap or the presence of features such as dead-branches below the living crown, thus altering the shape of the fitted function. Predicted HTLC (including the correction) RMSE and CV RMSE were similar for all four study sites.
HTLC estimates for RW18 increased in terms of CV RMSE when using the correction method, indicating a possible difference in vertical physical structure, which caused problems. A total of 9 of the 11 plots that employed the correction were from locations that had been thinned 9 years before. This may account for the differences in vertical structure for the thinned plots, but not for the other two that were unthinned. This implies that the arbitrary threshold between dominant canopy maxima and minima below of 6 m is not applicable in the case of the RW18 site. As observed in Figure 9(a), there is a degree of heterogeneity of HTLC within each field plot extent. Crowns of different sizes, spatial distributions and 3D positions will undoubtedly influence the proportion of lidar returns from various vertical strata and thus the shape of the curve function being fit to the data. Alternatively, for those plots with a higher canopy density could exhibit a less open structure for laser pulses to penetrate through, resulting in fewer returns from within the canopy and the lower strata (e.g. Chasmer, Hopkinson, and Treitz 2004;Dean et al. 2009), and therefore larger estimation errors as the shape of the curve reflects this. The majority of these plots locations were covered by two flight lines that incorporated larger scan angles (≥12°). It is possible that the vertical distribution of returns was influenced by this (Holmgren, Nilsson, and Olsson 2003). It should be noted that the various metrics at the plot level were estimated for plot sizes and shapes, which varied considerably over the four study site extents, and may introduce more heterogeneity in canopy heights, in addition to the possibility of edge effects. Future work would be needed to further develop the rules for applying the aforementioned correction.
The results from this study compare well with other studies for estimating plot-level average HTLC either statistically or using the estimates derived from voxel-based curvefitting analysis. For example, Andersen, Mcgaughey, and Reutebuch (2005) produced a linear regression model that produced an R 2 value of 0.84 and RMSE values between 3.9 and 4.1 m. The following studies utilize layer detection techniques. Muss, Mladenoff, and Townsend (2011) were able to estimate HTLC with an R 2 of 0.71 and RMSE of 1.52 m by fitting a frequency-based curve. In contrast, Morsdorf et al. (2010) estimated the top coniferous layer base height from sub-canopy broadleaved species layers with a bias of 0.07 m (±1 m standard deviation) by separating clusters of returns based upon height and intensity values.
While a quantitative validation for the sub-dominant canopy layers detected within the course of this project was unavailable, and was not present in the majority of the plantation study site areas due to management regime. The Henderson site data set, however, contained a mix of sub-canopy layers due to their respective initial management types. Through the inspection of the vertical profile, it was possible to determine the presence and potential height of the sub-canopy and understory layers, which corresponded to known field-assessed presence and absence. Plot locations that received the initial five years of vegetation control were all noticeably different from the plots that received no herbicide use. Given the nature of the experiments at Henderson, the pairing of vegetation control and no-control plots for each of the initial site treatment types allowed an easy comparison of vertical profiles between plots. As illustrated in Figures 6 and 7, the vertical position of the dominant canopy can alter between sites, where those subject to herbicide use were typically taller. Such a technique that is able to compare locations may be of benefit to forest managers, for example in assessing the influence of understory competition on tree growth. Further work needs to be implemented to validate the quantification of the heights of understory components.

Individual tree crown analysis
The results of the ITC delineation approach used within this study cannot be validated for actual tree-specific characteristics due to insufficient field data. Instead a number of plot-level summary values were calculated. These were the number of ITC objects, average ITC top height, and average ITC HTLC within the plot extent.
The number of Loblolly pine stems recorded for each field plot extent was compared to the number of ITC objects generated within the same polygon extent. A relatively high R 2 value was obtained (0.85); however, the number of ITC objects was almost always lower than the field equivalent, with larger underestimation in locations with higher number of stems or larger plot sizes (e.g. RW19 and SETRES). One site, Henderson, had a number of plot locations with a larger number of ITC objects than field-recorded stems. The level of correspondence increased when comparing against dominant and co-dominant stems only; however, the relationship was not exact. As stated in Lee, Lucas, and Brack (2004) and Kaartinen et al. (2012), many of the image-based methods, which attempt to define ITC objects from a raster CHM, have a number of challenges relating to how well ITC objects can be delineated due to the complexity and unpredictability of the shape and geometry for groups of tree-crowns, and no consideration of features obscured by the overstorey (especially so when tree crowns in the dominant canopy merge).
The estimation of plot-averaged canopy height from ITC analysis was slightly inferior to that of the plot-level curve-fitting analysis in terms of RMSE, even though the R 2 was high (0.96). Presumably this is due to the uncertainty inherent in the delineation of ITC objects and the disparity with the number of stems, which introduced a bias. This is also a cause for concern for the estimates of average HTLC; however, additional concerns are raised when considering the sampling of lidar points for each ITC circular object, where returns may be located in overlapping, or merged, tree crowns and structures beneath the canopy. The results for estimating HTLC produced, however, higher correlations (R 2 = 0.96) as is evident in both Figure 6, and in the smaller RMSE, which was 0.4 m smaller, when compared to the plot-level curve-fitting approach. While not a direct comparison, due to the spatially explicit nature of ITC-level analysis, these results in the current study correspond favourably to the results reported by Holmgren and Persson (2004). The aforementioned study utilized a linear discriminant function to estimate HTLC for a coniferousdominated forest area, producing an R 2 value of 0.84 and RMSE value of 2.82 m. Similarly to the current project, Popescu and Zhao (2008) were able to estimate crown base height for individual tress through a combination of analysing vertical profiles and regression models, producing an R 2 value of 0.8 and RMSE of 2.03 m.
There is still uncertainty in estimating the actual heterogeneity of HTLC at the plot level through ITC analysis, as evidenced by Figure 10 (SETRES) and especially by Figure 9 (RW18). It should be noted that the fieldwork for RW18 and RW19 sites employed some degree of subsampling for each plot. This introduces uncertainty concerning the actual heterogeneity of stems encountered. Indeed, the larger variability in lidar-derived estimates of crown length highlights the uncertainties that lidar returns are actually incident on features of interest at an appropriate scale, i.e. a tree crown, rather than from the stem or competing vegetation, or indeed this method's ability to successfully classify the position of the crown. With reference to RW18 (Figure 9), far more variability in lidarestimated HTLC is observed. The most variability in HTLC value was encountered in thinned plots for RW18, both above and below the expected range of values. Some of the difference in accuracy is no doubt due to the resolution of the height bin size utilized (1.0 m) in addition to the low proportion of number of lidar returns available for analysis at the ITC scale and their likely biased distribution of returns towards the top of the canopy. However, there are challenges relating to the difficulties in delineating meaningful 3D objects and separating them from one another in the point cloud and excluding unwanted features (e.g. returns from the tree stem) in complex environments. RW18 represented the most extreme example, whereas the SETRES site ( Figure 10) more closely reflected the pattern of field values caused by management type, i.e. for the four groups of four plots with different treatment types of irrigation and fertilizer additions.
It can be inferred that the ITC-level analysis can somewhat account for the heterogeneity encountered at the plot level caused by individual stem HTLC variability over the age range of 13-26 years and multiple management types in order to provide an accurate plot-level average estimate in variable plot sizes. Therefore, at the plot level there seems to be a benefit in estimating HTLC by aggregating ITC predictions, or by further subdividing the plot. Future work should aim to validate the estimates of HTLC at the individual tree level, as with studies such as Popescu and Zhao (2008), and to attempt to improve the 3D delineations.

Considerations and future work
Within the context of relatively homogenous planation forest the statistical and plotlevel curve-fitting approaches appear to be sufficient; however, more heterogeneous environments may fail to be accurately characterized. Further sub-dividing the plot extent in some manner and attempting to estimate canopy characteristics could potentially provide more accurate results, for HTLC at least, given an appropriate data set.
Forest management type can have an impact upon the forest's vertical structure and by extension the distribution of lidar returns within that extent, at least for plantation Loblolly pine forest. This is indicated by the vertical profiles fit to the Henderson plotlevel data and derived ITC HTLC variability, where canopy and sub-canopy vertical features varied by management type. Whilst there are a number of issues to overcome, such as at smaller scales (ITC), the methods in the current project have provided the potential extraction of features pertinent to assessing the response to management over variable plot sizes.
This voxel-based curve-fitting method employs the caveat that the lidar data must penetrate the vegetated canopy sufficiently in order to be able to detect the required features. Small-footprint full-waveform lidar systems, which offer a much higher potential for detecting returns beneath the canopy (Wagner et al. 2006), may offer a number of advantages to this type of method in the future. Lindberg et al. (2012) demonstrated that incorporating a correction based upon the Beer-Lambert law applied to the waveform intensity could reduce the effects of decreasing returns from beneath the canopy. Lindberg et al. (2012) were able to produce vertical profiles from both waveform and discrete return lidar, in a similar approach to that presented here, in order to estimate vertical vegetation structure for one or more layers (e.g. the canopy and shrub layers). Estimates from the processed waveforms were slightly more accurate, implying that further processing of the waveform may improve prediction accuracies.
Plot-level curve-fitting-based estimates of HTLC could be both above and below field measurements. The depth of laser pulse penetration will vary with canopy structural characteristics and lidar device configuration (Gaveau and Hill 2003). Therefore, the estimation error could be a consequence of there not being enough laser pulse energy returning to the sensor sufficient for recording as a return (Harding et al. 2001), or indeed an issue prevalent in discrete return lidar, where there is a 'blind spot' following each detected return (up to 1.2-5 m), during which nothing can be detected (Reitberger, Krzystek, and Stilla 2008). Alternatively, the occlusion of lower canopy components as reported in Chasmer, Hopkinson, and Treitz (2004) could be considered. One potential explanation for the underestimation is the inclusion of returns from features such as dead-branches or from the tree stems beneath the living crown altering the curve shape minima and maxima locations. Alternately, the occurrence of HTLC heterogeneity within the plot could also impact the comparison as noted previously. Future work incorporating different forest structural arrangements and species, and including explicit spatial information from individual trees would be necessary in order to explore this source of uncertainty.
Validation of the sub-dominant vertical layer heights is needed beyond a mere presence and absence in future work, for example in order to determine fire fuel load (Sandberg, Ottmar, and Cushon 2001), and the understory's connection with nutrient cycling (Noss 1990). Areas which are heterogeneous in structure and composition, or contain vegetation layers with less distinct vertical boundaries, also impose a number of challenges for future work. Other studies have made use of lidar return intensity values to stratify the returns beneath the forest canopy, in order to separate the vertical components of the coniferous overstorey, broad-leaved understory, and shrub layer (Morsdorf et al. 2010). The study outlined in Wing et al. (2012) determined understory presence by stratifying returns using field-derived height information, allowing the mapping of understory attributes at the plot-level.
Capturing variation in these vertical attributes at smaller scales, such as that of ITC, can produce accurate results for dominant canopy layer characteristics; however, subdominant layers may not be consistent in height or horizontal coverage at the plot level, so far little attempt has been made to map of information from 'hidden' vertical layers (e.g. understory) at relatively small scales.

Conclusions
The results of this study show that airborne lidar data can be used to estimate the average canopy to height and HTLC at the field plot level. The approach was intended to characterize the vertical structure through fitting a Nadaraya-Watson kernel regression curve to the vertically binned frequencies of returns. While not the only potential technique available, the method described in the current study appears to be effective for delineating the dominant canopy vertical layer at the field plot scale.
The results of the automated curve-fitting approach at the resolution of the field-plot have produced high correlations with fieldwork. Higher relationships were achieved for the majority of cases when a correction was applied in order to alter where the automated procedure considered the dominant layer vertical to be, which depended on the form of the section of curve corresponding to the dominant canopy layer. In addition, the method including the correction produced similar level of accuracy when compared to conventional regression-derived estimates for both average canopy height and HTLC (e.g. Naesset 2002Naesset , 2009).
Many of the study site plots were established in single-storey forest locations; however, the Henderson site contained locations with variable levels of understory due to differing site management. For each of these eight management types, a vertical profile and curve were created, where there were more features detected from the subdominant portion of the profile in areas that had no herbicide used in the first five years of growing. Likewise few or very small features were detected for locations where herbicide was initially used. This implies that statistics beneath that of the dominant canopy may be derived given an appropriate method.
An approach for automatically delineating ITCs from a raster CHM was developed. lidar returns within a circular radius of each ITC centre were extracted and Nadaraya-Watson kernel regression curve was fit to the vertically binned-frequencies of returns. Those ITC objects intersecting with each field plot extend were averaged to produced plot-level estimates for comparison. Plot-level comparisons were implemented as individual tree comparisons could not be explored. The number of delineated ITC objects was lower than the number of recorded Loblolly pine stems in the majority of plot extents. Whilst estimates of average canopy top height were comparable for the plotlevel and averaged ITC curve-fitting methods, the latter produced a higher overall correlation with HTLC.
Within the context of relatively homogenous planation forest the statistical and plotlevel curve-fitting approaches appear to be sufficient; however, more heterogeneous environments may fail to be accurately characterized. Even within the context of this study, many of the plots exhibited some degree of HTLC variation within. While a number of limitations in the current ITC approach were highlighted, such as heightbin resolution, sampling overlap, under-and over-estimation of HTLC effects at the ITC scale that increased uncertainty, high correlations were produced. This implies that subdividing a forest plot-location with heterogeneous HTLC and aggregating analysis at a smaller scale could improve estimates.
The approach is dependent on the capacity of the lidar pulses penetrating the vegetated canopy and with a density of returns beneath the canopy sufficient to identify features, where dense vegetation may impede its overall accuracy. With a potentially higher return density than conventional discreet returns systems, small-footprint fullwaveform lidar data may be of great benefit for future work. In addition lidar return intensity could potentially assist in identifying compositional differences between layers (as in Morsdorf et al. 2010).