S. Magnussen, F. Gougeon, D. Leckie and M. Wulder

Natural Resources Canada, Canadian Forestry Service, Pacific Forestry Center, 506 W. Burnside Road, Victoria, British Columbia, Canada, V8Z 1M5.

Email: {smagnuss| fgougeon| dleckie| mwulder @pfc.forestry.ca}

Tree heights can be derived from LIDAR (Light Detection And Ranging) data of canopy heights when an estimate of the expected number (N) of ‘tree top’ hits exists. The probability of a LIDAR pulse hitting a tree top is governed by the number and size of ‘tree tops’. Counts of trees with a sunlit crown were derived from Compact Airborne Spectrographic Imager images with 55 x 55 cm pixels. Pairing this count with an assumed ‘tree top’ size of 3 pixels suggested an average of 6 ‘tree top’ hits in a 400 m² plot. The topmost N_hit canopy heights can be considered as proxies of tree height. This simple concept is demonstrated in a thinning trial with 39-year-old Douglas-fir. Extreme value distributions fitted to expected ‘tree-top’ hits generated various statistics of predicted dominant and codominant tree heights. Our approach lowered the average discrepancy between mean canopy height and Lorey’s mean tree height from 20% to an insignificant 3% (P=0.09).

Ground-based tree-counts improved the results slightly.

Accurate canopy height profiling of forest stands is now a reality due to improvements of differential global positioning systems, laser scanner technology, and inertial navigation systems (Flood and Gutelius 1997, Ritchie 1995, Blair et al. 1994, Bufton 1989). Canopy height data provide valuable information about the forest stand, canopy structure and possibly crown closure (Nelson 1997, Ritchie et al. 1993). However, to estimate mean stem volume and site productivity of a stand a measure of stand height is needed. Canopy height is the vertical distance from an arbitrary point on the forest floor to the topmost part of the canopy at the chosen location. Tree height is the vertical distance from the base of a tree to its topmost terminal bud (Avery and Burkhart 1983). The relationship between the two can be very complex. For dominant and codominant trees the mean canopy height is bound to be lower than the mean tree height (Avery and Burkhart 1983). Crown shape, stem density and spatial distribution of the trees will all influence the discrepancy.

Laser scanning of forest stands provides the mean canopy height of dominant and codominant trees. The chance that a laser pulse hits an intermediate or even a suppressed tree is low (Magnussen and Boudewyn 1998). Converting canopy heights to tree heights can be done by either a model-based approach (Magnussen and Boudewyn 1998) or through heuristic methods consisting of selecting the largest canopy height within spatially binned areas (Nśsset 1997, Nilsson 1996). Extreme values of canopy heights are more likely to be associated with a tree top than less extreme values. Both approaches have demonstrated their validity. The advantages of a model-based approach include a reduced reliance on auxiliary information from ground-based surveys. Model-based approaches, on the other hand, tend to be complex and computationally demanding (Magnussen et al. 1999). This study demonstrates how tree counts and an estimate of the combined area occupied by tree "tops” (target area) leads to improved estimates of tree heights from laser canopy height data.

Data were collected from one 50-ha stand of 49-year-old second-growth Douglas-fir near Shawnigan Lake ( 48^o 38' N, 123^o 43^' W) on Vancouver Island of British Columbia. The stand was fairly homogenous and on a gentle slope. The stand contained a thinning and fertilization trial (Brix 1993) with 42 square plots of 400 m² each. Geographic coordinates for each plot were obtained by bearing and distance from a known control point.

All trees in each plot were measured in the winter of 1994-1995 (age 49) for diameter at breast height (DBH cm, 1.3 m above ground). At that time the plots contained between 32 and 158 trees per plot (median 64, mean 66, standard deviation 31.9). Fifteen trees per plot representing the plot range of DBH values were also measured for height (HT). HT of the remaining trees were predicted from plot-specific nonlinear height diameter relationships (McWilliams and Therien 1997). DBH and HT values of all trees were prorated to their expected 1995-1996 status by means of individual and trait-specific relative growth rates (RGR). The RGR values used were those registered for the period between the last two measurements (1989/90 and 1994/95). A control sample of 46 trees was established to check the predicted height growth. A t-test of predicted versus measured heights did not reject the hypothesis of no difference (P=0.29). Heights of dominant and codominant trees were computed as Lorey’s height (Löetch and Haller 1964), i.e. the basal-area weighted tree height.

Crown shape data were obtained from a 1992 biomass study of 120 trees from 12 plots. The study furnished 2214 paired measurements of branch length (BL) and the distance (DT) between the stem location of the branch and the top of the tree (Brix 1993). Branch length (BL) was found to be a power function of distance to the top (DT). The average relationship was BL = DT^0,421 (root mean square error = 0.19, r²=0.91). This functional relationship was used to compute the expected total crown area at various canopy depths.

2.3 CASI data

The Compact Airborne Spectrographic Imager (CASI) is an airborne sensor capable of acquiring multispectral imagery in the visible and near infrared part of the spectrum (Anger et al. 1994). CASI was flown to acquire a high spatial resolution multispectral image (nine spectral bands, mostly ~14 nm wide) over the study site. Using an onboard inertial navigation system, differential GPS technology and proprietary algorithms, the supplier provided the image already geo-referenced to UTM coordinates and at a nominal spatial resolution of 55 cm/pixel.

2.4 Digital stem counts on CASI images

Stem counts were obtained from the CASI image using a technique originally developed to detect and classify the species of mature trees in medium resolution airborne sensor images (1-2 m/pixel) (Gougeon and Moore 1989). The technique consists of detecting local maxima in a smoothed image from the most appropriate spectral band (typically the near infrared). In medium to dense coniferous stands, where individual tree crowns appear separated by areas of shade, the algorithm generally isolates a single pixel per tree, usually corresponding closely to a well illuminated tree top. An appropriately chosen threshold suppresses local maxima in the shaded areas. Recent applications of this technique includes regeneration assessment in dense, moderate, and even open stands using 30 cm/pixel MEIS images (Gougeon, 1997). Similar techniques have been used by Eldrige and Edwards (1993), and Dralle and Rudemo (1996). In this study we improved the stem count technique by "doubling the spatial resolution” (to 27.5 cm/pixel) before smoothing the image.

Laser scanner data were acquired in December 1996 with the Optech ALTM 1020 laser scanning system supplied with a Trimble 4000 SSE GPS receiver. Two flight lines were each flown in a forward direction then backward recording distances to the first and the last laser pulse, respectively. Flying altitude was about 600 m with a ground speed of 40 mˇs^-1. Laser pulses were emitted at a wavelength of 1047 nm with a frequency of 2 to 8 kHz. Scan angles were less than ą12^o.

Two data files were generated from the flights and some proprietary post-flight processing: one containing 25835 triplets of latitude (X), longitude (Y), and elevation (Z) of first returned pulses and the second 17413 triplets (X, Y, Z) of last returned pulses. The former is assumed to constitute ‘canopy hits’ and the latter ‘ground hits’. We refer to these data as (X_c , Y_c , Z_c) for the imputed canopy hits, and (X_g , Y_g , Z_g) for the imputed ground hits. Positional errors across the scanned area are 1 to 3 m for x and y and about 20 cm for z (Magnussen and Boudewyn 1998). That is, all hits from a single flight line may have a constant location error of 1 to 3 m in either horizontal direction. Relative positions of points gathered in a single over-flight are accurate to within about 20 cm. Six of the 42 plots were only partially scanned by the over-flight, leaving 36 plots for the analyses. Further details about the data capture are available (Magnussen and Boudewyn 1998).

With a beam divergence of 0.25 mrad, an average flight height of 600 m, and an average canopy height of 18 m the laser footprint of canopy hits is a circle with radius about 0.14 cm.

Canopy heights were computed as
where is a prediction of the ground elevation at the location (x_c , y_c) obtained from a digital elevation model (DEM) fitted to the ground triplets (Magnussen and Boudewyn 1998). DEM residuals had a standard error of about 20 cm and appeared to be normally distributed.

Scatter plots of versus indicated a large number of near ground canopy heights. Thus we conjectured that the canopy hits were a mixture of ‘near ground’ hits and actual canopy hits. Application of an Estimation-Maximization (EM) algorithm (Dempster et al. 1977) based on two distinct linear relationships between and in this mixture enabled us to mark 14.6% of the imputed canopy hits as ‘near ground’ hits.The ‘near ground’ data were discarded from the analyses; Magnussen and Boudewyn (1998) provide further details. Canopy heights (CH) were then equated to for the remaining 22061 records. When allocated to the 36 plots each field plot contained between 29 and 71 estimates of canopy height (median 44, mean 48, standard deviation 11.4).

2.7 Expected number of laser pulses hitting a tree "top”

An estimate of the number of laser pulses hitting a tree "top”

is obtained from:

(1)

where

denote the number of laser pulses hitting the canopy within an area A,

an estimate of the number of visible tree tops within area A, and the average horizontally exposed area of a tree top. Exposed in this context simply means that the crown area is visible from above. Implicit in (1) is the assumption that laser pulses and trees are randomly distributed within A (Cressie 1991). Note in (1) includes the convolution of the actual target area (tree top) with the spatial extension of the laser pulses (Stoyan et al. 1987).

CASI-based estimation of and

For the CASI-based estimates of derived from the automated stem detection, a target area of 3 pixels was assumed. While the tree top proper surely occupies less than one pixel the convolution with the laser footprint and allowance for hits ‘near’ (within the topmost annual growth whorl) the top makes the assumption reasonable.

Ground-based estimation of and

Ground-based tree-counts () while readily available need to be converted before they enter (1). In general, as some trees are obscured from the view of the laser pulses by adjacent trees. Also, the average target area

must be determined from considerations based on the average geometry of tree crowns.

, in this context, is defined as the horizontal projection of the space defined by the uppermost whorl and the terminal bud of a tree (a cone). The projected space is assumed to be circular with a radius equal to
where is the average current height increment of the trees (Magnussen and Boudewyn 1998). For a sub-sample of 46 randomly selected trees we found which gives a target radius of about 0.6 m. The effective target radius is slightly larger if we allow half a laser footprint inside a target area to trigger a hit. From these considerations we arrived at

1.7 m².

Conversion of to

relied on considerations of canopy geometry. We essentially computed the probability that a tree would have an exposed tree top (P_ET ) and then derived .

An estimate of P_ET was obtained from simple geometric considerations of the vertical distribution of target areas (tree tops) and crown areas. First we assumed that target and crown areas were distributed between a height of 10.2 and 26.2 m (range = 16.0 m). These limits are the observed minimum and maximum laser canopy height in the 36 field plots. Actual limits obtained from ground data would be 5.5 to 31.7 m. Next, we posited a symmetrical bell-shaped distribution of target and canopy areas. Magnussen and Boudewyn (1998) found support for a beta distribution with parameters ą (= 2.1) and ˛ (= 2.1). Finally, we assumed that the projected crown profiles were circular with a radius (CR) determined by the distance to the tree top (DT) as in CR = DT^0.421 (Magnussen and Boudewyn 1998). By computing the ratio of the expected target area to the expected crown area in each of a series of canopy segments we essentially estimate the probability that a target area is exposed within a segment. Canopy segments were defined as a horizontal canopy slice with a width of 0.29 m equal to the average height increment. Note the independence of view angle: the area multipliers for target and crown areas due to off-nadir view angles are identical and cancel when the taking ratios. Summing these probabilities for all canopy segments yields the final probability. In short:

(2)

where the symbols denotes the largest integer smaller than k and the smallest integer larger than k. pdf(Beta( ą ,˛),(x-10.2)/14.0) is the probability density function of a beta distribution with parameters ą ( = 2.1) and ˛ ( = 2.1) evaluated at an interior point x of the interval [10.2; 26.2]. The last term in (2) is a division by the total probability of a target area.

2.10 Extracting the topmost canopy heights (CH)

To ensure a spatially balanced sample of the topmost canopy heights a double stratified sampling scheme was applied. First, each plot was subdivided into

sections or strips along an East-West direction, and then sectioned along a North-South direction. Subsequently, for each of the two sectioning directions the topmost laser canopy height in each section was selected as a proxy height of dominant and codominant tree heights. Thus two samples of the topmost

canopy heights were generated from each plot. The two samples, although not independent of each other, are treated as simple replicates in the statistical analyses.

Statistical inference for the selected topmost canopy heights is limited by small sample sizes (here about 5 per plot and direction of sectioning), the data (unequal number of canopy laser heights in each section), and the selection process (the topmost). More robust and extended inferences are possible if we can obtain the expected distribution of the selected canopy heights. According to extreme value theory (Smith 1986), when a sample of n canopy heights are independent and identically distributed and there exist and such that the largest canopy height in the sample (CH_1:n ) when scaled to converges in distribution to a random variable with a non-degenerate cdf, and that cdf is one of three types (Frechet, Weibull, or Gumbel). The convergence occurs as . In other words, if we can derive estimates of the two parameters (a,b) we would have a proxy distribution of TH. For convenience and ease of estimation the Gumbel model was chosen to quantify the distribution of extreme canopy heights.

The maximum canopy height in a sample of n canopy heights (CH_1:n ) is assumed to have the following Gumbel (G) probability density distribution (pdf) with parameters a and b (Smith 1986):

From this proposition it follows that plot specific estimates of the Gumbel parameters a and b based on observations of the largest canopy height within each of sections in a plot allow us to predict the pdf of heights of dominant and codominant trees. The mean tree height and other descriptive statistics of dominant and codominant trees are derived from this pdf by conventional calculus (Johnson et al. 1994).

Plot specific maximum likelihood estimates of a and b were obtained by maximizing the following likelihood (Strand and Boes 1998):

where subscript refers to plots () and

is the total number of included laser canopy heights in plot and is the number of included laser canopy heights in the kth section of plot. As a rule we excluded the lower half of the laser canopy heights from each plot section. This exclusion of the lower half had virtually no impact on the estimated parameters but improved the computational aspects considerably, especially convergence. Whenever fell below three, the number of plot sections were lowered by one, and so on until all sections contained at least three canopy heights. For these plots n in (3) was lowered accordingly. Maximization of (4) was done with a Newton-Raphson algorithm (McCulloch 1997).

Table 1 summarizes the ground-based maximum likelihood results. CASI-based results were similar (not shown). All estimates of a and b were more than 3 times as large as their estimated standard deviations. The correlation between estimates of a and b was strong and negative (-0.80). Estimates from the East-West and the North-South sectioning were, in general, within 10% of each other. Occassionally (four plots), when the spatial distribution of canopy heights were very uneven, larger discrepancies emerged. For each plot we used the simple average of the East-West and the North-South estimates as the parameters of a plot specific Gumbel pdf of presumed height of dominant and codominant trees.

Figure 1 shows the laser canopy heights in a randomly chosen plot (no. 10) by East-West and North-South partitioning of the plot. In this particular example only the topmost 4 canopy heights in each section were considered for fitting. Figure 2 provides further details on the estimated Gumbel pdf and the selected topmost laser canopy heights.

Given estimates of a_i and b_i, the predicted mean height of dominant and codominant trees in the ith plot was:

(5)

where is mean number of canopy heights used for estimation in plot . The divisor in (5) is the total probability integral of the Gumbel pdf. Integration limits were arbitrarily set to 0 and 40 m. Variance estimates of were computed by replacing the integrand h by .

Ground- and CASI-based estimates of the number of ‘visible’ trees (

) showed considerable divergence. The minimum, mean, and median CASI-based counts were about 50% above the estimates derived from purely geometric considerations. This is hardly surprising because the geometric approach treated the crowns and tree tops as solids and implicitly assumed trees to be randomly distributed when in reality tree crowns are less than opaque (Carlson and Ripley 1997), and trees in this trial tend to be spaced more regularly than random (, based on a Chi-Square goodness-of-fit statistics for the quadrat method of testing of complete spatial randomness, Cressie 1991). Number of quadrats per plot was 20. Magnussen and Boudewyn (1998) argued for a crown transparency of 50% (that is only half the crown projection area is occupied by foliage, branches and stems) which would help narrow the gap between the two estimates of . Less divergence emerged in the expected number of tree top hits . Here the ground-based estimates surpassed the CASI-based estimates by about 20%. On average we expect 4.7 tree top hits given the ground-based approach and 5.8 given the CASI-based counts. Given intrinsic differences in spatial resolution and hence definitions of target area, the concordance is fortuitous.

Estimated average dominant and codominant tree height derived from the expected number of tree tops hit and plot specific Gumbel pdf’s were generally within one meter of estimated Lorey’s heights (Figure 3). Simple means of the topmost canopy height from each plot were equally close to their ground counterparts (not shown). Ground-based estimates were slightly closer to Lorey’s heights than the CASI-based estimates. The average difference of 13 cm was, however, not statistically significant (

). In six plots the discrepancy between estimated and actual height exceeds 3 m. Previous analyses (Magnussen and Boudewyn 1998) made it clear that the six plots could reasonably be considered as outliers. The number of laser canopy hits were in all six instances well below average and spatially clustered towards one end of the respective plots.

On average, Lorey’s mean height of 22.4 m exceeded the ground-based estimate by 60 cm. A simple two-sided t-test did not reject the hypothesis of no difference (. The corresponding mean absolute deviation (MAD) was 1.34 m. Parallel CASI-based results were 73 cm () with MAD at 1.39 m and a correlation of 0.62. Standard deviations of the recovered heights, as obtained from the fitted Gumbel pdf’s varied from 0.7 to 1.6 m (mean 1.0m). Thus most (30 of 36) estimated plot means of dominant and codominant tree height had an estimated confidence interval that included Lorey’s height. In contrast, mean laser canopy heights, were, on average 4.3 m shorter than Lorey’s height ( and MAD averaged 4.3m

A one-tree difference in the mean number of expected hits between the ground- and CASI-based method explicates the consistent off-set of about 20 cm in the point estimates. The off-set is consistent with the differential selection intensity (Magnussen 1999). In the CASI-based approach 1 laser canopy height were selected out of an average of 8 while 1 of 10 were selected in the ground-based approach. Assuming that heights of dominant and codominant trees are nearly normally distributed the stated difference in selection intensity alone accounts for 85% of the offset (see Falconer 1981 for details on the effect of selection differential on expected means of selected items).

Part of the persistent downward bias in our recovered tree heights may actually be explicated by the selection procedure of the topmost laser canopy heights in each plot. Our partitioning of each plot and selection of only the topmost in each section inadvertently yields lower heights than a selection of the topmost canopy heights within the entire plot. Given the estimated standard deviations of recovered height we expect a priori a downward bias of about 20 cm. However, we rationalize that the spatial partitioning is consistent with the averaging over all trees in a plot.

For laser scanning of forest canopies to become a realistic forest inventory alternative to a direct (traditional) measurement of tree heights it is key that easy-to-use robust methods for extracting tree heights from canopy heights become available (Magnussen and Boudewyn 1998, Nśsset 1997, Nilsson 1996). The typical significant discrepancy between canopy heights and tree heights prevents the direct use of canopy heights as a proxy for tree height. Site index and volume estimates based on mean canopy heights would be seriously downward biased. Results from this study confirmed that easy-to-use methods transform laser canopy heights into good approximations of desired tree heights.

Magnussen et al. (1999) showed that data processing based on double-blind deconvolution is capable of removing the average bias between laser-based and ground-based estimates of tree heights, albeit at a higher price. Additional costs are to be paid for analysis and computing.

Reasonable approximations of tree height have been obtained from laser canopy data by selecting only a few topmost canopy heights from within a plot (Tickle et al. 1998, Nćsset 1997, Nilsson 1996). Local experiences have guided the establishment of sampling guidelines and the results indicate that estimated tree heights are typically within 0.5 to 2 m of the ground-based results when trees are between 12 and 20 m tall. Magnussen and Boudewyn (1998) found support for this practical approach in a thinning trial with Douglas-fir by demonstrating that the vertical distribution of canopy heights obtained from airborne laser scanning closely resembled the vertical distribution of foliage surface. Thus, selecting a quantile of the canopy heights that matches in probability the quantile of foliage surface area corresponding to a desired tree height would produce unbiased estimates of tree height.

A model-based approach to estimate tree heights directly from canopy heights has the advantage of establishing a framework for optimization and statistical inference. Derivation of tree heights from canopy heights without auxiliary information is feasible (Magnussen et al. 1999) but the necessary deconvolution algorithms are computationally intensive.

Auxiliary ground information or information from other remotely sensed images such as CASI can greatly reduce the computational complexities of recovering tree heights from laser canopy heights. A pivotal parameter to be gained from the auxiliary information is the number of tree tops and size of tree crowns. Pragmatism must prevail in the definition of ‘tree top’. A reasonable definition should be possible in most circumstances.

Once an estimate of this number has been obtained, it is straightforward to either use this number of topmost canopy heights from each plot as proxies for tree heights or, as demonstrated, as input to fit a distribution model of extreme values. The difference is bound to be inconsequential. Fitting a distribution (say Gumbel) offers the additional advantage of sufficient statistics of the estimated quantities.

5. Acknowledgments

The CASI imagery was acquired under a Forest Renewal of British Columbia supported project entitled "Development of Certified Forestry Applications Using Compact Airborne Spectrographic Imager (CASI) Data", a co-operation between MacMillan Bloedel Ltd., ITRES Research Ltd., and the Canadian Forest Service (Pacific Forestry Centre).

Estimated number of visible tree tops per plot (, see (2))