Use of Stereology in Forest Inventories—A Brief History and Prospects for the Future

Abstract

Several forest inventory techniques utilize approaches that are similar to stereological approaches often applied in microscopy and other fields. Stereology is characterized by the description and estimation of properties of objects based on samples of lower dimension than the object, e.g., 2-dimensional slices from 3-dimensional objects, 1-dimensional probes from 3-dimensional or 2-dimensional objects and dimensionless points from higher dimensional objects. The stereological character of many forest inventory methods was historically developed independently of recognition of a relationship with stereology. Strip sampling of forests, common in the late 19th and early 20th century, can be considered as a sterelogical approach if the strip centerline is viewed as a 1-dimensional probe of tree inclusion zones on a land area. The stereological character of plot sampling and Bitterlich sampling becomes evident if one views these methods as samples of 1-dimensional probes for volume within tree inclusion zones, or dimensionless points sampling for basal area in inclusion zones. Traditional methods of estimation of tree stem volume include samples of 2-dimensional cross-sectional area at fixed points along the tree stem to estimate 3-dimensional volume. Though these traditional methods usually use a shape assumption (e.g., parabolic frustum) for short stem segments, we show how a random-systematic start estimator of stem cross-sections can provide a design-unbiased estimate of stem volume without using any stem shape assumptions. Monte Carlo integration estimators of tree volume such as importance sampling that are designed to depend on only a few (usually one or two) tree upper-stem height or cross-sectional samples can also be viewed as stereological methods. Several forest inventory methods such as Matern’s individual tree basal area estimator and sector sampling can be viewed as local stereology, in which sample lines or slices pass through a central point. Finally, we suggest potential applications of stereological principles in the emerging “big data” era characterized by lidar and other remote sensing data and the assemblage of large tree and stand datasets. We suggest a new stem volume estimator which may have potential for future use with terrestrial lidar.

1. Introduction

Many forest measurement and inventory techniques are related to stereological methods because they frequently use samples of lower dimension, such as points, lines, or tree cross-sectional areas, to estimate quantities of higher dimension such as 3-dimensional tree or forest cubic volume. However, for most of their histories, the sciences of stereology and forest measurements have developed independently so the relationship between them has only been recognized in recent years (e.g., [1] (pp. 195–201), [2]) but to our knowledge this relationship is still rarely discussed in the forest mensurational literature.

Historically, many stereological methods have arisen from the microscopic study of sample slices. These methods are often characterized by the use of 2-dimensional slices of 3-dimensional objects and 1-dimensional probes to make inferences (such as volume estimates) for 3-dimensional objects. The word “stereos” in the Greek language means “solid” and is the source of the name “stereology.” Many sterelogical methods can also be used to make inferences for 2-dimensional regions based on 1-dimensional probes or dimensionless points. Currently, sterelogical inferences are also used in imaging technology in medical science, geology, and materials science.

According to Baddeley and Jensen [3] (p. 3) “Modern stereology may be defined as ’sampling inference for geometrical objects.”’ Stereology is “...especially useful when a sample has lower spatial dimensions than the original material. Hence, stereology is often defined as the science of estimating higher dimensional information from lower dimensional samples” (International Society for Stereology 2017). Many forest inventory methods have these same characteristics. For example, Bitterlich sampling uses randomly-located dimensionless points on land areas to estimate forest volume. Traditionally, 2-dimensional stem cross-sectional area samples have been used to develop individual tree stem volume estimates for individual trees.

Since the late 1920s and 1930s, there has been an active relationship between the statistical sciences and forest mensuration. Influential in initiating this relationship were seminars and training sessions arranged by F. X. Schumacher (then a United States Forest Service scientist) and colleagues featuring R. A. Fisher, one of the primary developers of widely used statistical techniques in the 20th century [4] (pp. 319–321). However traditional statistical survey sampling as classically developed, for example, by Cochran [5] (pp. 4–8) posits a finite population composed of a finite number of units from which sample units are selected to make statistical inferences. This theoretical framework did not formally include the possibility that sample units could be of lower dimension than the estimated quantities resulting from statistical inference, and so differs formally from the stereological paradigm. Nevertheless, it was often possible to formulate forest sampling problems in ways that made the traditional statistical survey sampling approaches useful (e.g., [6]). Although many traditional expositions of statistical methods are based on infinite population models, often based on the normal distribution, they usually have not explicitly recognized the possibility of making inferences about 3-dimensional objects based on lower dimensional samples.

During the 1980s, applications of Monte Carlo integration to forest mensurational problems began to appear, such as importance sampling (IS) for tree bole volume estimation [7]. Hammersley and Handsomb [8] (pp. 6–8) state that the term “Monte Carlo” method along with its development date from about 1944 and the origination of these methods stems from the work on the atomic bomb during the second world war. Although, of course, Monte Carlo methods are quite broad, Monte Carlo integration is one application. It is difficult to integrate functions that are not regular, especially higher dimensional integrals, using traditional numerical methods such as Simpson’s rule and the trapezoidal rule, so Monte Carlo integration techniques have been developed to accomplish this [9] (pp. 114–115). Monte Carlo integration essentially views the integral as an expected value and samples values of the integrand within the domain of integration in order to estimate the value of definite integrals. Valentine et al. [10] show how approaches to forest sampling using lines or points can often be interpreted in the context of Monte Carlo integration.

Our objectives will include demonstration of the relationship between stereology and several forest inventory methods by using examples at the individual tree level, the forest stand level, methods that can be viewed as local stereology, and estimation of coarse woody debris. We also propose two new design-unbiased systematic random-start stereological estimators for tree stem volume that have connections to traditional model-based estimators. These will be tested using simulation methods. Finally, we discuss prospects for future applications of stereological methods to forest sampling problems.

2. Methods

Estimators with stereological properties have been developed at both the individual tree level and at the stand level. Some of the these methods such as strip cruising and fixed-size plot sampling have histories going back to at least the 19th century, while other techniques, such as Monte Carlo integration methods for individual tree volume and surface area estimation have been developed in the late 20th and early 21st centuries.

We will see that several estimators that have been developed for forest sampling can be viewed stereologically as using samples of dimension lower that the dimension of the quantity to be estimated. In forest sampling, 0-dimensional points and 1-dimensional probes have often been used to estimate 2-dimensional quantities such as area and 3-dimensional quantities such as cubic volume. The estimation of 3-dimensional volumes has also been accomplished by using 2-dimensional sample “slices” in the context of forest sampling.

Foresters have long used 1-dimensional measurements to estimate 2-dimensional quantities such as basal area and other stem cross-sectional areas, as well as 3-dimensional stem volume. Most commonly, a model has been used to convert 1-dimensional measurements to 2-dimensional areas and 3-dimensional volumes. The circle has been the most common model for the estimation of tree cross-sectional areas based on diameter measurements, however, its deficiencies have been recognized by writers such as Matérn [11] and more recently by Pulkkinen [12]. Dot grids and line transects were often used to compute areas on maps and aerial photos [13] (p. 84), practices which have analogs in stereology. These methods are less widely used currently due to the availability of Geographic Information Systems (GIS) and related software. The polar planimeter was once widely used to compute maps and other 2-dimensional areas based on 1-dimensional curve length [13] (pp. 84–86). The method was based on the principle of a line integral.

2.1. Stereology Preliminaries

A fundamental stereological formula for volume estimation is known as Crofton’s formula [3] (p. 72). The following version of Crofton’s formula relates to the estimation of 3-dimensional volume by sampling 2-dimensional slices, which is frequently the paradigm for estimation of individual tree stem volume:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{V}\left(Y\right)={\int }_U\mathcal{A}(Y\cap T_u)du\end{array}\end{array}$$

(1)

where $\mathcal{A}$ is area, Y is the 3-dimensional area, and $T_u$ is a plane section through Y at $u\in U$. A special case of the use of Crofton’s formula would be the familiar “slices” solid of rotation integral to obtain stem volume from an individual tree taper curve (e.g., [14] (pp. 180–181)). Generalizations of Crofton’s formula can handle other situations such as surface area. Generalized forms of Crofton’s formula may be found in [3] (p. 75) and [15] (pp. 72–75).

Because the areal “slice” in Crofton’s formula can itself be conceived as an integral over one of the dimensions in which the solid is embedded, the formula fundamentally depends on Fubini’s theorem (e.g., [16] (pp. 269–270)) which indicates that double integrals can be computed as iterative “single” integrals, and that the order of integration for double integrals may be interchanged. Let us consider a tree stem composed of cross-sections $c\in \mathcal{C}$ at various stem heights $h\in \mathcal{H}$, and $I(c,h)=1$ if the stem contains a cross-section c at height h, otherwise zero, for a tree stem with a volume measure of $V_{\mathcal{S}}$. Fubini’s theorem allows us to view the stem volume in the following two ways for sampling purposes:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill V_{\mathcal{S}}={\int }_{\mathcal{C}\times \mathcal{H}}I(c,h)d(c,h)={\int }_0^H{\int }_0^{C}I(c,h)dcdh={\int }_0^{H}c\left(h\right)dh\end{array}\end{array}$$

(2)

where $c\left(h\right)$ is the cross-sectional stem area at h, $\mathcal{C}\times \mathcal{H}$ is the product of sets, $d(c,h)$ may be considered an infinitesimal volume and H is total height. Alternatively, the Fubini theorem allows us to interchange the order of integration and use the iterated integral:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill V_{\mathcal{S}}={\int }_{\mathcal{C}\times \mathcal{H}}I(c,h)d(c,h)={\int }_0^C{\int }_0^{H}I(c,h)dhdc={\int }_0^{C}h\left(c\right)dc\end{array}\end{array}$$

(3)

where $h\left(c\right)$ is the upper-stem height to stem cross-section c and C is the cross-section at the stem base. However, most forest measurement applications work from Equation (2); as will be seen below, critical height sampling and certain other estimators are based on Equation (3).

Cavalieri’s estimator of volume is named after the Italian mathematician who noted that two solid objects which have equal cross-sectional areas on every plane have equal volume [3] (p. 15). The estimator is a systematic random-start estimator based on a series of equally-spaced parallel section planes [3] (pp. 65–66). Following [3] (pp. 65–66), a stack of parallel planes with spacing s can be denoted as:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathcal{T}=\{T_{U+js}:j\in Z\}\end{array}\end{array}$$

(4)

where $T_u$ is a horizontal plane at height u and the $j^{th}$ plane of the stack is $T_{\left(j\right)}=T_{U+js}$ with the index j a signed integer. Although the stack contains infinitely many planes, a particular estimation requires only a finite number of them. The stack starting position is defined by U randomly and uniformly chosen over the interval $[0,s)$. The Cavalieri estimate of volume is then [3] (p. 66):

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{V}\left(Y\right)=s\sum _j\mathcal{A}(Y\cap T_{\left(j\right)})\end{array}\end{array}$$

(5)

This can be viewed as a finite-sum approximation to the integral in Crofton’s formula Equation (1) and hence to the volume of the solid.

2.2. Individual Tree Estimators

Foresters have usually obtained the cubic volumes of tree sections by using model-based estimators such as Huber’s, Newton’s or Smalian’s formula [13] (pp. 97–104). The equation below is an example of a traditional approach to tree stem volume estimation with n stem diameter measurements at fixed intervals of length L meters using Huber’s formula [13] (p. 103)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{V}=k\times L\left(d^2(L/2)+d^2\left((3/2)L\right)+\cdots +d^2\left((n-1/2)L\right)\right)\end{array}\end{array}$$

(6)

where $\widehat{V}$ is the estimated stem cubic meter volume, $k=\pi /\mathrm{40,000}$ and $d\left(\right(i-1/2\left)L\right)$ is $i^{th}$ sectional stem diameter measurement in centimeters at stem height $(i-1/2)L$. Huber’s formula is based on a diameter measurement at the midpoint of a stem section and is exactly correct for the frustum of a paraboloid [13] (pp. 99–101). For other stem section shapes, there will be a bias, although the bias may be modest for short stem sections.

2.2.1. Systematic Random-Start Estimation of Tree Bole Volume

A stereological design-unbiased estimator can be constructed that is superficially similar to Equation (6) above. Let the first diameter measurement $d(u\times L)$ be made at a point $u\times L$ from the base of the stem where u is randomly-selected from a uniform distribution on the interval $[0,1)$. Then let the $i^{th}$ diameter $d\left(\right(u+i-1\left)L\right)$ be measured at a point $(u+i-1)L$ from the base of the stem:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{V}=k\times L\left(d^2(u\times L)+d^2\left((u+1)L\right)+\cdots +d^2\left((u+(n-1))L\right)\right)\end{array}\end{array}$$

(7)

where $\widehat{V}$ is the estimated stem cubic meter volume. This is a systematic random-start sample of upper-stem cross-sectional areas that is design-unbiased for actual cubic meter stem volume V. It is clearly an example of the use of the Cavalieri volume estimator Equation (5) based on the integral volume representation (2). The Estimator (7) is essentially the same as a stereological systematic random-start sample for solid volume based on Cavalieri’s principle presented by Baddeley and Jensen [3] (pp. 65–68). We tested the properties of this estimator using simulation to assess variability for an example tree stem, $h=2$ m, and stem shapes including neiloid, cone, and paraboloid as well as a shape intermediate between a cone and a paraboloid known as a paracone [17]. A tree stem 30 m long with basal diameter 37.5 cm was selected. Both Estimators (6) and (7) could be used on standing trees with upper-stem dendrometer measurements or on felled tree stems. This method is similar to the stereological method of volume estimation by local vertical planes in the “vertical rotator” [3] (pp. 200–202).

Another option for systematic random-start sampling would be to systematically sample n heights $h_c\left(c\right)$ (the distances from stump to stem cross-sections c) at constant cross-sectional area increments ${\Delta }_c$ with an initial random decrease from stump diameter $\delta$ randomly selected on the uniform interval $(0,{\Delta }_c)$.

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{V}={\Delta }_c(h_c(\delta \times {\Delta }_c)+h_c((\delta +1)\times {\Delta }_c)+\cdots +h_c((\delta +(n-1))\times {\Delta }_c))\end{array}\end{array}$$

(8)

This estimator is based on the integral volume representation of Equation (3).

We also used simulations to test Estimator (8) on the same example tree sizes discussed above. The motivation for this is similar to Grosenbaugh’s proposal for height accumulation ([13] (pp. 103–104), [18]), that the method might better sample against variation in taper. However, the height accumulation procedure uses constant taper steps rather than cross-sectional area steps. In order to compare the method of constant cross-sectional area steps to the use of constant length increments in Equation (7), the stem cross-sectional area increment ${\Delta }_c$ was determined so that the number of stem measurements was the same as when $h=2$ m was used in Equation (7) for a stem having basal diameter 37.5 cm and length 30 m. This resulted in ${\Delta }_c=0.04$ square meters.

2.2.2. Stereology and Crude Monte Carlo Integration for Tree Bole Volume

It is often possible to formulate the forest and tree attributes which forest mensurationists want to estimate as integrals, e.g., tree bole volume has often been represented as an integral based on a solid of rotation [13] (pp. 90–91):

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{V}=\frac{\pi }{\mathrm{40,000}}{\int }_0^H{d}^2\left(h\right)dh\end{array}\end{array}$$

(9)

where $\widehat{V}$ is tree bole volume in cubic meters, $d\left(h\right)$ is diameter in centimeters at height h on the tree bole and H is total height in meters. Clearly, the above is an example of volume integral representation (2). Crude Monte Carlo (CMC) [8] (pp. 50–52) integration (with no variance reduction techniques) could estimate the volume by selecting n random samples of the integrand with uniform probability over the interval $(0,H)$:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{V}=\frac{kH}{n}\sum _{i=1}^n{d}^2\left(h_i\right)\end{array}\end{array}$$

(10)

where $k=\frac{\pi }{\mathrm{40,000}}$ and $d\left(h_i\right)$ is the ith sample diameter at height $h_i$. Figure 1a,b illustrates the concept of Crude Monte Carlo integration for individual tree stem cubic volume where the dark lines represent sample slices of the tree stem.

This unbiased estimator is sometimes called the sample mean Monte Carlo estimator [9] (p. 119). Note that the sample values of the integrand $d^2\left(h\right)$ are 2-dimensional quantities that are being used in the estimator Equation (10) above to estimate the 3-dimensional cubic volume V. Therefore, this Monte Carlo integration estimator also falls within the definition of a stereological estimator which uses lower-dimensional samples to estimate higher-dimensional volumes. Estimator (10) is essentially the same as stereological estimators based on Cavlieri’s principle [3] (p. 15) which states that two solid objects having the same cross-sectional area on every horizontal plane will consequently have the same volume. In many instances, quantities such as the cubic volume of solids or 2-dimensional areas can be formulated as integrals in which the integrand is of lower dimension. Thus, it would seem that many Monte Carlo integration estimators could also be considered to be stereological in nature because the value of a definite integral often has higher dimension than its integrand. Of course Monte Carlo integration can be used to estimate the value of definite integrals that cannot be interpreted as areas or volumes and in those instances Monte Carlo integration estimators might not be closely related to stereological methods.

Several Monte Carlo integration estimators for tree stem volume have been proposed. As indicated above, there is a strong relationship between stereological estimation methods and Monte Carlo integration estimators when Monte Carlo estimation is applied to a 3-dimensional object such as a tree stem and the integrand can be computed from a measured tree dimension such as stem diameter.

2.2.3. Stereology and Variance Reduction Techniques from Monte Carlo Integration

The Crude Monte Carlo estimator Equation (10) is not practical to apply directly to tree stems because it would require an excessive number of upper-stem diameter measurements in the field. However, variance reduction techniques can be applied to substantially reduce the variance of the Monte Carlo integration stem volume estimate so that a reasonable estimate can be made with only one or two upper-stem diameter measurements. Variance reduction techniques work well for Monte Carlo integration estimators for tree volume because we have a great deal of prior knowledge concerning tree stem shape once a DBH and total tree height have been measured. Perhaps the earliest proposal for variance reduction in Monte Carlo integration estimation of tree stem volume was Importance Sampling [7]. Importance sampling uses a proxy taper function to sample upper-stem cross-sectional areas with probability density $p\left(h\right)$ proportional to volume $V_p$ as computed by the proxy taper function, where $d_p\left(h\right)$ is upper-stem diameter at stem height h predicted by the proxy taper function:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill p\left(h\right)=\frac{k{d}_p^2\left(h\right)}{V_p}\end{array}\end{array}$$

(11)

The importance sampling stem volume estimator can be expressed as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\widehat{V}}_{HIS}=\frac{1}{n}\sum _{i=1}^n\frac{k{d}^2\left(h_i\right)}{p\left(h_i\right)}=\frac{1}{n}\sum _{i=1}^n\frac{d^2\left(h_i\right)}{d_p^2\left(h_i\right)}{V}_p\end{array}\end{array}$$

(12)

where ${\widehat{V}}_{HIS}$ is the IS estimator of stem cubic volume. The importance sampling estimator is design-unbiased for stem volume [7]. Like the sample mean Monte Carlo estimator above, it is similar to stereological estimators based on Cavalieri’s principle [3] (p. 15). Because samples are selected with probability density proportional to proxy volume, sample upper stem diameters are more likely to occur on the lower portion of the tree stem than with the sample mean Monte Carlo Estimator (10) and the variance of the importance sampling estimator is substantially lower than the variance of the sample mean Monte Carlo estimator for the same number n of sample upper stem diameters [7]. Figure 1c,d shows the placement of stem sample slices (dark lines) with IS tending to be lower on the stem than the CMC sample Figure 1a,b.

Control variate estimators provide another approach to variance reduction in Monte Carlo integration estimation of tree stem volume. The control variate estimator ${\widehat{V}}_{HCS}$ of tree stem volume for n randomly-located upper-stem sample diameter measurements is:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\widehat{V}}_{HCS}=H\times \frac{k}{n}\sum _{i=1}^n{d}^2\left(h_i\right)-\beta [H\times \frac{k}{n}\sum _{i=1}^n{d}_p^2\left(h_i\right)-V_p]\end{array}\end{array}$$

(13)

where $\beta$ is a constant. Rubenstein [9] (p. 127) gives the equation to compute the optimal value of $\beta$ but the information needed to compute this value is often not known prior to sampling, so a value of $\beta =1$ has sometimes been used. Valentine et al. [19] used control sampling with $\beta =1$ to compare vertical and horizontal modes of importance and control sampling for bole estimation. Lynch [20,21] used control sampling with $\beta =1$ in an assessment of the effects of upper-stem measurement error on bole volume estimation using Monte Carlo integration techniques including importance sampling and control variate estimation. If the optimal value of $\beta$ is estimated, then the control variate estimator is essentially the same as the classical regression estimator of statistical survey sampling [5] (pp. 189–190). Like the importance sampling Equation (12) stem volume estimator, the control variate estimator is based on an estimator that is similar to stereolgical volume estimators that trace their roots to Cavalieri’s principle [3] (p. 15). Applications of control variates to tree stem volume have been discussed by Van Deusen [22] and Gregoire and Valentine [23].

Antithetic variates are another variance reduction technique that has been applied to individual tree stem volume estimation. Antithetic pairs of upper-stem diameter measurements can be obtained by selected the measurements at heights $h=uH$ and $h_a=(1-u)H$ where u is a uniform random variate selected from the interval $[0,1)$. Due to typical decrease in tree stem diameter as stem heights associated with these diameter measurements increase, these pairs of antithetic diameter measurements have a strong tendency to be negatively correlated, so that the variance of the sum of these two antithetic diameter measurements will be less than the variance of two independent stem diameter measurements. An estimator based on the sample mean Monte Carlo Estimator (10) is:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{V}=\frac{kH}{2n_a}\sum _{i=1}^{n_a}[d^2\left(h_i\right)+d^2\left(h_{ai}\right)]\end{array}\end{array}$$

(14)

where $n_a$ is the number of antithetic pairs, $h_i=u_{i}H$, $h_{ai}=(1-u_i)H$ and $u_i$ is the $i^{th}$ uniform random variate drawn from the interval $[0,1)$, $i=1,\dots ,n_a$. Thus within a pair, a larger value of $d\left(h_i\right)$ tends to be associated with a smaller value of $d\left(h_{ai}\right)$, inducing a negative correlation among antithetic pairs.

The antithetic variates technique can be combined with other variance reduction techniques such as importance sampling. Van Deusen and Lynch [24] found that a combination of antithetic variates and importance sampling reduced variance in stem volume substantially and was better than the average of two independent importance samples. A disadvantage of antithetic variates for standing tree stem volume estimation is that one must make at least two time-consuming individual tree upper-stem diameter measurements, although it is likely that one of them may be on the lower bole and within reach of sampling personnel.

Although the examples of Monte Carlo integration individual tree volume estimators presented above are based on selection of stem cross-sectional areas at randomly-selected stem heights, one may also develop stem volume estimators that select sample upper-stem heights based on randomly-selected upper-stem cross-sectional areas. Fundamentally, this is due to the fact that Fubini’s theorem (e.g., [16] (pp. 269–270)) allows one to interchange the order of integration and then to apply Crofton’s formula Equation (1) on the resulting integral. As a result, it is possible to sample on the basis of volume integral Equation (3). Individual tree volume estimators of this kind (based on Equation (3)) were developed by Lynch et al. [25] by applying importance sampling and antithetic variates. They found this approach worked well if one applies both importance sampling and antithetic variates, however, as indicated above, antithetic variates requires the measurement of two upper-stem sample heights.

2.3. Forest-Level Estimators

For many years, foresters have employed stand-level estimators of per-hectare forest attributes such as volume, number of trees and basal area per hectare which could be interpreted as stereological estimators although this has not been widely noted in the past. Going back at least to the 19th and early 20th century, strip sampling and plot sampling methods have been used to estimate per-hectare quantities of forest attributes.

Strip cruising is a forest sampling technique that samples all trees within a fixed distance of a centerline, traditionally marked by a forester’s chain, where the centerline crosses a forested property at fixed intervals (e.g., [26] (pp. 276–289)). The strip cruising method can be related to stereology if we associate with each tree a “bracket” having width equal to the strip width and perpendicular to the direction of the centerline (typically strip cruising employs equally-spaced centerlines that are parallel to each other). Then, we may view the centerline as a stereological “probe” that samples trees if the probe crosses the “bracket” associated with the tree. Alternatively, we may posit circular inclusion zones about each tree in the population and sample trees when the strip cruising centerline crosses into a tree inclusion zone and extends at least as far enough into the inclusion zone to cross a radius perpendicular to the centerline. Lynch [27] described a scheme similar to this in the context of Monte Carlo integration for sampling riparian areas, although in that case the width of the “bracket” varied proportionally to tree dbh.

Plot sampling samples all trees within a fixed-sized plot of land. Circular plots are probably the most popular in forestry although square and rectangular plots have also been used. Fixed-radius circular plots sample all trees within a fixed-length radius of a randomly or systematically located plot center within the forest population being sampled. In the first quarter of the 20th century in the USA, Chapman [26] (p. 297) mentions plot sampling for forest inventory but seems to assume that strip crusing would be most widely used at that time. Chapman [26] (pp. 312–314) believed that permanently established plots should be square or rectangular.

To view plot sampling from a stereological perspective, it is useful conceptually to shift the plot center to the tree, and to consider that a tree is selected when a randomly or systematically-located plot in the field falls into the plot associated with that tree. This was termed the “tree-centered” view of plot sampling by Husch et al. [13] (p. 225) and has come to be more commonly termed the “inclusion zone” in the forest sampling literature (e.g., [23] (pp. 223–224)). Under this concept, the fixed-sized plot method of forest sampling uses a 0-dimensional point to sample 3-dimensional quantities such as cubic meter volume or 2-dimensional quantities such as basal area per hectare. As indicated above, using lower-dimensional points or probes to estimate higher-dimensional quantities is an inherently stereological procedure.

2.3.1. Plot Sampling

A fruitful way to conceive of circular plot sampling stereologically is to associate a cylinder centered on each tree in the forest population. Let the base of the cylinder have area ${\mathcal{A}}_p$ (the area of the fixed-radius sample plot) and the height of the cylinder equal to $V/{\mathcal{A}}_p$ where V is individual tree cubic volume. Then, we may envision the sampling problem as sampling the total volume of tree cylinders in the population by using a 1-dimensional probe. Where tree cylinders have bases that overlap, they are to be “stacked.” The average volume per hectare of the forest population is then found by the average depth of n randomly or systematically located sample probes multiplied by the sample plot area ${\mathcal{A}}_p$. This concept is illustrated by Figure 2 below.

2.3.2. Bitterlich Sampling

The concept of point sampling, often termed “Bitterlich sampling” after its originator, Walter Bitterlich [28] (pp. 1–6) depends on selecting trees whose cross-sections at breast height subtend a fixed angle which has a vertex at a randomly or systematically located point in the field [13] (pp. 220–223). It is well-known that forest stand basal area per hectare can be estimated by multiplying a count of sample trees from point sampling [13] (p. 220) by the Basal Area Factor ($BAF$):

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{B}=\frac{BAF}{n}\sum _{i=1}^n{m}_i\end{array}\end{array}$$

(15)

where $\widehat{B}$ is estimated basal area per acre, n is the number of sample points and $m_i$ is the number of trees subtended by the point sampling angle at point i.

It can be shown that individual tree inclusion zones in point sampling are circular and proportional in size to individual tree basal area. One may view the basal area per hectare estimator as a process of counting the number of randomly or systematically located points that fall into tree inclusion zones. This is equivalent to “hit or miss” Monte Carlo integration [9] (pp. 115–119) where “hits” happen when sample points fall into tree inclusion zones. Point counts are also widely used in stereology. Point sampling for estimation of basal area per hectare with systematically located field points is very similar to stereological “planar point grid” methods [3] (pp. 158–159). The same principle is employed by the “dot grid” Avery and Burkhart [14] (pp. 79–80) once commonly used by foresters for area determination on maps and aerial photographs, though now largely superseded by GIS technology.

Bitterlich or point sampling can be used to estimate forest volume by summing the sample tree volumes weighted by tree basal area and multiplying by the basal area factor:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\widehat{V}}_{PS}=\frac{BAF}{n}\sum _{i=1}^n\sum _{j=1}^{m_i}\frac{v_{ij}}{g_{ij}}\end{array}\end{array}$$

(16)

where ${\widehat{V}}_{PS}$ is point sampling estimated stand volume per hectare, $v_{ij}$ is the volume of sample tree j at sample point i, and $g_{ij}$ i is the basal area of tree i at sample point j. As in the plot sampling illustration of Figure 2, the method could be viewed stereologically as measurement of the height of a probe because the inclusion zone for each tree in point sampling is proportional to tree basal area g. Figure 3 shows this conceptually using the volume sampling surfaces [29,30] for each tree as generated by the R [31] package ‘sampSurf’ [32].

2.3.3. Critical Height Sampling

Critical height sampling measures the height at which sample tree cross-sectional area fits exactly into the point sampling angle projected from the sample point in the field. Again, this can be viewed as the use of 1-dimensional vertical stereological probes placed randomly or systematically in a forested area. The average probe depth multiplied by forest land area would be a stereological estimate of cubic wood volume. To increase the chance of interception by a sample probe, the point sampling angle is used to expand the stem at every height, leading to the concept of the expanded tree stem in critical height sampling as explained by Iles [33] and illustrated in Figure 12.15 of Iles [34] (p. 578). This concept is illustrated via a sampling surface approach in Figure 3b.

The critical height tree cubic volume estimator originally developed by Kitamura [35,36] is the product of the $BAF$ and the sum of critical heights:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\widehat{V}}_c=\frac{BAF}{n}\sum _{i=1}^n\sum _{j=1}^{m_i}{h}_{c_{ij}}\end{array}\end{array}$$

(17)

where ${\widehat{V}}_c$ is the critical height stand volume estimator and $h_{c_{ij}}$ is the critical height of sample tree j at sample point i. Lynch [37] and Van Deusen and Meerschaert [38] showed that critical height sampling is based on computation of stem volume from the “cylindrical shells” integral by integrating with respect to tree radius or cross-sectional area rather than height. Recall that the Fubini theorem allows for interchanging the order of integration in Equation (9) and then the Crofton Formula (1) can be applied. Note that cross-sectional area is a transformation of stem radius r by the formula $c=\pi r^2$. Thus, critical height sampling is implicitly based on integral Equation (3).

Lynch and Gove [39] used an antithetic variate with importance sampling to address the problem that critical heights occur very high on the stem for trees close to the sample point. By using the antithetic variate associated with critical height, positions are reversed so that antithetic critical height occurs low on the stem near the sample point and high on the stems of trees distant from the sample point.

2.3.4. Ueno’s Method

Ueno [40,41] developed an estimator based on critical height sampling in which the cubic volume is estimated on the basis of counting the number points randomly generated in 3-dimensional space which fall within the expanded tree stems of trees in the forest of interest. This method is similar to stereological methods of counts using 3-dimensional grids as discussed in Baddeley and Jensen [3] (pp. 159–160). One locates n sample points in the field randomly or systematically. At each sample point, a point sampling angle gauge is used to tally sample trees. For each sample tree, a random height between 0 and a pre-specified maximum $H_m$ is compared to the critical height $h_c$ and the tree is selected if the random height is below critical height. Savings of time in the field come from the fact that one does not have to measure the critical heights of all sample trees because in many cases the randomly-generated height will be obviously greater or less than critical height. Ueno’s forest volume estimator can be expressed as [42]:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\widehat{V}}_u=\frac{H_{m}BAF}{n}\sum _{i=1}^N\sum _{j=1}^n{Y}_{ij}\end{array}\end{array}$$

(18)

where ${\widehat{V}}_u$ is Ueno’s estimate of forest volume per hectare, N is the number of trees in the forest of interest and $Y_{ij}=1$ if tree i is selected at field sample location j, $Y_{ij}=0$ if tree i is not selected at field sample location j. Ueno’s method is illustrated in Figure 4a.

Lynch [42] showed how the Monte Carlo integration variance reduction techniques of importance sampling, control variates, and antithetic variates could be applied to Ueno’s method. The cubic volume estimator for Ueno’s method with antithetic variates is:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\widehat{V}}_{ua}=\frac{H_{m}BAF}{2n}\sum _{i=1}^N\sum _{j=1}^n(Y_{ij}+Y_{ija})\end{array}\end{array}$$

(19)

where ${\widehat{V}}_{ua}$ is the cubic volume estimate from Ueno’s method with antithetic variates, $Y_{ij}=1$ if tree i is sampled at point j and the random height $u{H}_m$ is below the critical height, otherwise zero, and $Y_{ija}=1$ if the antithetic random height $(1-u)H_m$ is less than critical height, otherwise zero.

Simulations using published diameter and height measurement times indicated that antithetic variates reduced the variance in Ueno’s method substantially, making it more efficient [42]. Simulations by Sterba [43] (also see [28] (pp. 139–141), [42]) indicated that Ueno’s method can be competitive with point sampling in which sample tree heights are measured, which might make Ueno’s method attractive for use where appropriate individual tree volume equations are not available. Figure 4b illustrates the sample selection process in Ueno’s method with antithetic variates.

2.3.5. Cylinder Sampling

Lynch [44] developed a method of sampling termed “cylinder sampling” in Lynch [42] which is also similar to the use of 3-dimensional grids in stereology. In cylinder sampling, a combination of 3-dimensional cylinders is constructed that are proportional to a combined variable volume equation [45] (p. 24), which is a linear equation in the variable $x=D^{2}H$. A cylinder with a cross-sectional area equal to each tree’s point sampling inclusion zone, and height equal to total tree height is constructed for each tree in the forest population. A smaller cylinder proportional to the additive constant term in the combined variable volume equation is either added or subtracted from each tree cylinder, depending on whether the constant is positive or negative. Sample points are located in 3-dimensions in a way similar to Ueno’s method above, by generating a random height between zero and maximum expected height $H_m$ for each sample tree. Trees are selected when that sample point falls within the sampling cylinder. The cylinder sampling cubic volume estimator for a combined variable equation with a negative intercept ($v=b_0+b_1{D}^{2}H$ where v is tree volume and $b_1$ and $b_0<0$ are constants) is:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\widehat{V}}_{cyl}=\frac{b_1{H}_{m}BAF}{n}\sum _{i=1}^N\sum _{j=1}^n{Y}_{ij}\end{array}\end{array}$$

(20)

where ${\widehat{V}}_{cyl}$ is the estimate of cubic volume using cylinder sampling, $b_1$ is the coefficient of $D^{2}H$ in the combined variable equation and $Y_{ij}=1$ if tree i is selected at sample point j, otherwise $Y_{ij}=0$. Similar methods are available where $b_0>0$ in the combined variable volume equation [44].

In simulations, Lynch [42] found that cylinder sampling was more efficient than ordinary point sampling in which the heights of all sample trees are measured, or Ueno’s method. However, Ueno’s method avoids possible bias in volume equations, and can be used for situations in which volume tables do not exist. Antithetic variates reduced the standard error of cylinder sampling but not sufficiently to improve the efficiency of cylinder sampling due to the extra time needed for another measurement. Figure 5 illustrates the sample selection process in cylinder sampling.

2.3.6. Stereology from the Perspective of Modifying Sampling Designs

From sampling theory, it is known that estimators of population parameters of interest, such as the population total or mean, must properly reflect the sampling design. In developing sampling schemes that support stereological principles, a reverse manner of thinking can be applied, i.e., the design can sometimes be modified so that the estimator becomes very straightforward. The reason may be to replace impractical or costly measurements on sampled elements with simple counts or measures such as lengths. Mandallaz [1] (p. 203) discusses “pseudo Horvitz–Thompson estimators” for probability proportional to length line intercept estimation from a stereological perspective. In the following, we demonstrate the principle by using the Horvitz–Thompson (HT) estimator (e.g., [23]) as a starting point.

With the HT estimator, a population total, $\tau$, can be estimated under any probability sampling design as

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{\tau }=\sum _{iϵS}\frac{y_i}{{\pi }_i}\end{array}\end{array}$$

(21)

where $y_i$ is the quantity of interest for the $i^{th}$ sampled element and ${\pi }_i$ its probability of being included in the sample; the sum extended over all elements of the sample S. For example, if the total volume of all trees in a forest area is estimated under a sampling design that selects individual trees into the sample, $y_i$ would be the volume of a sampled tree and ${\pi }_i$ the probability by which the particular tree was included in the sample. Normally, we would need to measure or accurately model, based on measurements of the predictor variables used in the model, the volume of each sampled tree and divide it by the tree’s inclusion probability. The value of the inclusion probability follows from the design, and in some cases additional measurements are required to calculate it. A well-known case is Bitterlich sampling [28] (pp. 1–8) for estimating the total cross-sectional area of trees in a forest stand at $1.3$ m stem height (known as the basal area). With an angle count device, a relascope, trees are selected into the sample with probabilities proportional to their cross-sectional areas. Thus, the cross-sectional area of a tree appears in both the numerator and the denominator of Equation (21) and since these terms will always cancel in the estimator, measurements become unnecessary and can be substituted by simple tree counts. Thus, a general principle for developing a design based on which stereological principles can be applied is to select elements into the sample in a manner so that the inclusion probability of an element is either fully proportional to the study variable or at least proportional to some dimension of the study variable. In the first case, the estimator would simply be based on element counts; in the latter case, it would be based on simplified measurements, such as measuring length when the target quantity is volume.

The principle embodied in Equation (21) is demonstrated in the development count-based estimators for basal area per hectare (Bitterlich sampling Equation (15)) as well as cubic volume per hectare (Ueno’s method Equation (19) and Cylinder sampling Equation (20)) as developed above. In each case, the probability of sample selection was proportional to the quantity to be estimated on sample trees. A good example of this principle applied to the problem of estimating volume of downed coarse woody debris is perpendicular distance sampling [46] in which a count-based volume estimator is developed. The importance sampling individual tree volume estimator reviewed above (Equation (12)) similarly selects sample tree measurements (upper-stem diameters) with probability density proportional to a value of the measurement “predicted” by the proxy taper function.

2.4. Local Stereology and Forest Sampling

Local stereology is the quantitative study of spatial structures which can be viewed as neighborhoods of reference points. According to Baddeley and Jensen [3] (pp. 192–207), local stereology frequently uses 2-dimensional lines or 2-dimensional slices through a common reference point to infer attributes of 3-dimensional objects (also see [15] (pp. 27–28)). An example could be the use of 1-dimensional probes through the nucleus of a cell used to determine cell volume. For example, the stereological estimator known as the nucleator ([3] (pp. 198–199), [15] (pp. 179–180)) is based on two perpendicular intersecting isotropic lines.

One of the earliest examples of a forest sampling application that could be considered local stereology is the Matérn [11] (p. 17) basal area (or more generally cross-sectional area) estimator which is based on measuring lengths of radial sample semi-lines from a common vertex within the tree stem where all lines are in the same horizontal plane. Let n sample radians ${\theta }_i$ be selected from a uniform distribution on $(0,2\pi )$ and measure $r\left({\theta }_i\right)$ the length of the radial semi-line from a fixed vertex interior to the tree stem to the exterior of the stem, where all samples are taken from the same vertex and lie in a common horizontal plane then:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{G}=\frac{\pi }{n\mathrm{10,000}}\sum _{i=1}^n{r}^2\left({\theta }_i\right)\end{array}\end{array}$$

(22)

where $r\left({\theta }_i\right)$ is the length of random sample semi-line at sample angle ${\theta }_i$ in centimeters and $\widehat{G}$ is the cross-sectional area estimate in square meters. This form of the estimator implicitly assumes that the boundary of the stem slice in which the lines are sampled forms a “star-shaped” set with respect to the common vertex of the sample semi-lines such that for any line emanating from the vertex and extending to the outermost boundary point on that semi-line, there will be no points exterior to the set on that line. Baddeley and Jensen [3] (pp. 194–195) have provided a formal definition of a star-shaped set and have indicated the importance of such sets for local stereology. It is reasonable to assume that most, but not all, horizontal stem slices can be used to construct star-shaped sets from some interior point which could serve as a common vertex for the estimator $\widehat{G}$. Of course, Matérn [11] (p. 17) indicates that the estimator is not currently practical for standing trees due to the fact that the interior vertex point is not accessible for such trees. However, perhaps that barrier can be overcome in the future by imaging devices which could possibly use the method to estimate cross-sectional area with a very large number of radial sample lines. Gregoire and Valentine [47] have applied Estimator (22) above to the problem of area determination for irregularly-shaped planar areas more generally including land area. They also suggested variance reduction techniques and methods for perimeter length estimation.

Another forest sampling technique that can be considered to be related to local stereology is the sector sampling method proposed by Iles and Smith [48]. The method samples one or more sectors having a common vertex and the same projected angle. When applied to forest populations, all trees within a sample sector area are selected as sample trees. Just as one may view point or plot sampling as a process of sampling tree-centered inclusion zones with random sample points, one may associate a “tree-centered” sector with each tree in the population. Trees selected by sector sampling are those for which a randomly located radius is included in the sector having angle $2\alpha$ radians centered on that tree. Then, the sector sampling estimator can be formulated as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\widehat{Y}}_A=\frac{\pi }{\alpha \times n}\sum _{j=1}^n\sum _{i=1}^N{A}_i{I}_{ij}\left({\theta }_j\right)\end{array}\end{array}$$

(23)

where ${\widehat{Y}}_A$ is the estimated total amount of tree attribute A, $A_i$ is the amount of tree attribute A associated with tree i of a population of N total trees sampled by random azimuth ${\theta }_j$ with ${\theta }_j$ selected from a uniform distribution on $[0,2\pi )$, $I_{ij}\left({\theta }_j\right)=1$ if $({\lambda }_i-\alpha )<{\theta }_j<({\lambda }_i+\alpha )$, ${\lambda }_i$ the azimuth of tree i as measured from the vertex (sometimes referred to as “pivot point”) otherwise $I_{ij}\left({\theta }_j\right)=0$ and n is the number of random azimuths ${\theta }_j$ selected for the sample. This process is equivalent to sampling all trees located in a randomly chosen sector having a vertex angle of $2\alpha$. The method could be suitable for sampling smaller forest areas or it could be used to select sub-samples from fixed-radius plots or perhaps other types of plots. Lynch [49] indicated how to apply variance reduction techniques including importance sampling, antithetic variates and control variates. Use of squared length of the radius $r^2\left({\theta }_j\right)$ between the sample vertex and the tract border might be advantageous as a control variate because long sectors might be expected to contain more sample trees. It would be desirable to place the sample vertex so that the tract area is a star-shaped set as described by Baddeley and Jensen [3] (pp. 194–195) and indicated above, although if the tract is not star-shaped, the sector sampling process could still work if field crews are able to exclude non-tract areas within sample sectors. The estimation process is illustrated in Figure 6 below.

Sector sampling can be classified as local stereology because all sample sectors are referenced to a common point. The method could be viewed as the use of a radial probe to select sample trees whose tree-centered sectors contain the probe.

Estimation of vegetative cover using radial line sampling is also a method of forest sampling that can be classified as local stereology. With respect to a single vertex point within the area to be sampled, a random azimuth ${\mathsf{\Theta }}_j$ is sampled from the interval $[0,2\pi )$. A radial sample semi-line emanating from the vertex is established. At the entry point $C_j\left({\mathsf{\Theta }}_j\right)$ of each vegetative clump i encountered on the line, the distance $c_i\left({\mathsf{\Theta }}_j\right)$ is measured from the vertex to that point. At the exit point $D_i\left({\mathsf{\Theta }}_j\right)$ of each corresponding vegetative clump i, the distance $d_i\left({\mathsf{\Theta }}_j\right)$ is measured. The estimate of vegetative cover $\widehat{C_j}$ for each sample line j according to Lynch [2] is:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\widehat{C}}_j=\pi \sum _{i=1}^{m\left({\mathsf{\Theta }}_j\right)}(d_i^2\left({\mathsf{\Theta }}_j\right)-c_i^2\left({\mathsf{\Theta }}_j\right))\end{array}\end{array}$$

(24)

where $m\left({\mathsf{\Theta }}_j\right)$ is the number of vegetative clumps crossed by the radial semi-line associated with random azimuth ${\mathsf{\Theta }}_j$. Note that the radial nature of the sampling necessitates computing the difference between squared distances between entry and exit points rather than the simple linear distance as would be required for ordinary line intercept sampling. The coverage estimate for n randomly chosen angles ${\mathsf{\Theta }}_j$ is then:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{C}=\frac{1}{n}\sum _{j=1}^n\widehat{C_j}\end{array}\end{array}$$

(25)

The radial sampling process is illustrated for one sample line in Figure 7 below.

Lynch [2] shows how importance sampling, antithetic variates, or control variates may be used to reduce the variance in this estimator under appropriate conditions. Like sector sampling, the radial line coverage estimator would probably be most appropriate for small areas or for sub-sampling for coverage area on fixed-sized plots. Lynch [50] discusses relationships between three radially-oriented forest sampling methods, sector sampling, radial line sampling coverage estimation, and radial line sampling for estimation of planar areas, and also the effects of randomly- versus purposely-located pivot points (pivot point being one of the names of the common vertex for radial sampling within the tract area).

2.5. Line Sampling

As indicated above, some of the traditional methods of forest sampling could be considered as forms of line intercept sampling and are similar to stereological sampling with 1-dimensional probes. “Strip cruising” could be considered as line intercept sampling of fixed radius circular inclusion zones surrounding potential sample trees. When lines are similarly used to intercept tree inclusion zones arising from point sampling (Bitterlich sampling), the technique known as horizontal line sampling results. Horizontal line sampling was proposed by Strand [51] and has been discussed by Husch et al. [13] (pp. 224–231). Mandallaz [1] (pp. 195–200) devotes Chapter 12 of his forest sampling book to a detailed development of forestry applications of line intercept sampling from a stereological viewpoint.

Strand [51] also proposed the method of vertical line sampling (also see [13] (pp. 224–231)), which could be interpreted as line intercept sampling of circular tree inclusion zones with radii proportional to tree height. Vertical line sampling can be used to construct a sample tree count estimator of cubic volume by constructing 3-dimensional inclusion zones based on vertical line sampling [52]. By measuring tree stem diameters where they intersect a vertical angle projected from a sample line, Minowa [53,54,55], [28] (pp. 46–48) constructed a forest cubic volume estimator that, like the critical height sampling estimator, does not require individual tree volume equations.

Ståhl [56] has used a line sampling concept to develop transect relascope sampling for quantification of coarse woody debris. Addition discussion of sampling for coarse woody debris is given in the subsection below. Foresters have applied line intercept sampling to the similar problem of sampling logging slash (e.g., [57]).

Sampling with 1-dimensional lines or probes can be done in stereological applications to estimate the lengths of features in a plane [3] (p. 163). The forester and statistician Matérn [58] proposed a similar approach to the estimation of forest road length. In stereological applications, systematically spaced parallel lines are used where the system has a random orientation. The “bed of nails,” a 3-dimensional system of parallel lines, can be used for stereological estimation of surface area [3] (pp. 166–168).

There is a large literature relating to line intercept sampling in forestry and ecology. Chapter 9 in Gregoire and Valentine [23] (pp. 279–325) provides more information regarding some of these applications. Gregoire and Valentine [23] (p. 336) show how line intercept sampling can be viewed in the context of Monte Carlo integration. Barabesi and Fattorini [59] developed a stereological approach to line intercept sampling that could be used in natural resource applications.

2.6. Coarse Woody Debris

Much recent work has gone into the development of sampling methods for down coarse woody debris (hereafter simply ‘logs’), building upon older work that concentrated largely on line intersect sampling methods. Here, we describe perpendicular distance sampling. Many other methods for sampling logs that can be cast in a stereological framework have also been devised. Two recent sources providing in-depth reviews of many of these sampling methods are found in Gregoire and Valentine [23] (Chapters 9–10) and Kershaw, Jr. et al. [60] (Chapter 12).

Perpendicular Distance Sampling

Perpendicular distance sampling [46] estimates each log’s contribution to the total volume on a tract by a simple count when a sample point falls within the log’s inclusion zone. This is an ideal situation, which is similar to the estimation of basal area from the expanded tree circle under Bitterlich sampling (Section 2.3.2). In addition, both critical height sampling (Section 2.3.3) and Ueno’s method (Section 2.3.4) are similar to perpendicular distance sampling [46]. However, critical height sampling requires the measurement of a sampled tree’s critical height (which is tied to the critical cross-section for a given basal area factor) in order to estimate volume through its spatially referenced (or constrained) [39] implementation of crude Monte Carlo from a vertical probe (at the sample point); this extra information allows the jump from point (line probe) to volume over neighboring dimension area which is implicit in the critical height. Interestingly, if a tree were laid down on the forest floor with its accompanying three-dimensional critical height sampling surface, which was then projected (collapsed) onto the forest floor in a two-dimensional outline (or coverage area) to form the inclusion zone, Ueno’s method and perpendicular distance sampling become conceptually quite similar. Like Ueno’s method, where a simple count suffices to estimate each sampled tree’s contribution to volume, in perpendicular distance sampling a simple count is also sufficient to represent a sampled log’s contribution to the total (or per unit area) volume on the tract; this contribution is the same for all logs under a constant design parameter.

The development of perpendicular distance sampling on a downed log in a stereological sense would derive from either Fubini’s theorem (2) or Crofton’s Formula (1) in terms of cross-sectional area: The volume of a log is the integral of its cross-sectional area along its full length. For any downed log, a vertical probe hitting it would sample the log. However, to make sampling the log more likely, imagine simply widening the area on the ground that can be sampled on one side of the log with a vertical probe (sample point) by multiplying the cross sectional area along the full length of the log by an expansion factor, termed $K_{\mathrm{PDS}}$ (the design parameter). For example, $K_{\mathrm{PDS}}=2$ would give a one-sided expanded sample area (inclusion zone) that is twice as wide as the cross-section at any point along the log. Double this ground area to include the same amount of area on the other side of the log to complete the inclusion zone (Figure 8a). Thus, if V is the true log volume, then the inclusion zone area is $2K_{\mathrm{PDS}}V$ (via (9)). When a random sample point falls within the inclusion zone (i.e., if the distance from the point to the log is less than $K_{\mathrm{PDS}}$ times the perpendicular—point to log—cross sectional area of the log), then the log is included in the sample. Its contribution to the total volume on a tract of area $\left|A\right|$ m${}^2$ is

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\tau }_{\mathrm{vol}}=\left|A\right|\frac{1}{2K_{\mathrm{PDS}}}\end{array}\end{array}$$

(26)

An estimator for the total volume over n sample points would simply replace the numerator of (26) with the count of all sampled logs over all points. A constant volume factor associated with each log sampled can be easily determined based on the expansion factor as $V_{\mathrm{PDS}}=\raisebox{1ex}{$\mathrm{10,000}$}\!\left/ \!\raisebox{-1ex}{$2K_{\mathrm{PDS}}$}\right.$. Each log sampled represents $V_{\mathrm{PDS}}$ m${}^3$ha${}^{-1}$, where $K_{\mathrm{PDS}}$ has units m${}^{-1}$. Thus, the volume sampling surface for an individual log within its inclusion zone is not only flat, it is also the same height for each log in the population for a given value of $K_{\mathrm{PDS}}$ (Figure 8).

Williams and Gove [46] provide a simple proof demonstrating that the above Estimator (26), is design-unbiased. The connection to volume in perpendicular distance sampling comes from the implicit use of cross-sectional area in the development of (26) through the $K_{\mathrm{PDS}}$ factor. The $K_{\mathrm{PDS}}$ expansion factor allows one to ignore the measurement of cross-sectional area in most cases as it will be obvious whether a log is sampled or not based on comparison of the diameter (often an ocular estimate is good enough) to a limiting distance table for the corresponding $K_{\mathrm{PDS}}$ factor [46]. Again, this is very similar to Ueno’s method, where, in most cases, one can judge whether the random height at the sample point is less than the critical height, making measurement unnecessary. An overview of other forms of PDS, including the most recent developments, is given in [61] and [60] (p. 414).

2.7. Prospects for Future Applications of Stereology in Forest Sampling

We believe that there are a number of opportunities to apply stereological thinking to forest inventory problems that are emerging in the twenty-first century and have the potential to become important applications in the future.

2.7.1. Lidar and Remote Sensing

Stereology may also inform certain aspects of remote sensing. In particular, lidar might be viewed as the use of 1-dimension laser probes in a forested population. In a typical lidar survey, there are numerous laser “hits” on forested vegetation including trees stems and branches. Just as stereological probes are used to estimate contents, curve lengths, and other attributes in the usual stereological applications, perhaps laser “probes” can be used in similar ways to estimate quantities of interest in forest populations such as tree stem contents, which are critical not only for forest products but for estimates of carbon contents of global forests. Lidar has been used to characterize tree heights and forest height distributions (e.g., [62,63]). Lidar is being integrated into forest inventory projects (e.g., [64]).

In addition to airborne lidar, a significant research effort in estimation of forest attributes using terrestrial lidar is underway (e.g., Ducey et al. [65]). These systems use ground-based rather than aerial based laser scanners.

Consider the problem of estimating individual tree stem volume from terrestrial lidar. To discuss the possible estimation of tree stem contents from lidar “hits” stereologically, we will call the expression below the “bottle brush” estimator for potential use with lidar systems, because radii from laser “hits” on the stem to a central axis within the stem look like spines on a “bottle brush”:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{V}=L\times \frac{\pi }{n}\sum _{i=1}^n{r}_i^2\end{array}\end{array}$$

(27)

where $\widehat{V}$ is the estimated volume of a tree stem segment of length L, $r_i$ is the radial distance between the ith laser hit on the stem segment and a central axis positioned within the stem segment, and n the number of laser “hits” on the stem segment.

This estimator is like Matérn’s basal area estimator [11] (p. 17) but expanded to 3-dimensions. Like Matern’s individual tree basal area estimator, it is probably not practical to apply at the current time but may be in the future. It would be a design-unbiased estimator if we could say points were randomly-generated from a process that selects height at random and then a random rotation to place the radius. However, there are probably some deviations from this because the random process is coming from “outside” and not from a random rotation about the central axis of the stem.

The mathematical framework for this estimator could be cylindrical coordinates, with a vertical coordinate approximating the center of the tree stem and polar coordinates emanating from this central axis. Because tree stems are often not perfectly straight, the central axis could be a curve or a “piece-wise linear” segmented series of straight lines. An example of a stereological estimator based on cylindrical coordinates is the “vertical rotator” [3] (pp. 199–202). The tree stem Estimator (27) above is similar to estimation of volume using local vertical planes in the stereological “vertical rotator” [3] (pp. 200–202).

Because tree “circles” with larger circumferences will be more likely to receive laser “hits,” other things being equal, it might be advisable to weight the radial samples in the estimator above inversely by the radius, resulting in:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \widehat{V}=L\times \pi \sum _{i=1}^n\frac{r_i^2/r_i}{{\sum }_{i=1}^{n}1/r_i}=L\times \pi \frac{{\sum }_{i=1}^n{r}_i}{{\sum }_{i=1}^{n}1/r_i}\end{array}\end{array}$$

(28)

This estimator must be regarded as an “ad hoc” procedure that is probably somewhat biased but it is possible that the bias will be small in practical applications. Tests could be performed in the field on the two estimators above with conic frustums of known cubic volume to evaluate their performance.

To account for the possibility that the probability of a stem hit might vary with height on the stem, a stratified sampling estimator could be used to stratify the sample stem hit by elevation above ground. Within each stratum, Estimator (27) or (28) could be applied to find the stem volume for that stratum. The addition of estimates for all strata could be used to find the total volume.

Basal area or other cross-sectional areas might be estimated with a sufficiently dense point cloud by using Estimator (27) with $L=1$ at a stem height of $1.3$ m or other height that may be desired.

2.7.2. Large Databases and “Big Data” Applications

The “big data” era which is evolving is resulting in the assemblage of large databases of individual tree measurements from many historical studies. Radtke et al. [66] have assembled a large database of individual tree taper, volume, weight, biomass and physical properties (e.g., specific gravity) for a number of tree species in the USA (tree data available for download at legacytreedata.org). This might make possible the use of stereological estimators based on nearest-neighbor individual tree measurements in cases where model-based estimates have been used previously.

Individual tree volume estimates may be made by using the average of k nearest-neighbors from a large database, where volumes are estimated to merchantable limits, if needed, using upper-stem diameter measurements from large databases of tree measurements. In some situations, this could replace the now more commonly used estimates based on regression models of tree content. If large databases were available from systematic-random start upper-stem measurements, Estimator (7) could be used to obtain a stereological type of individual tree stem volume estimate. Perhaps unfortunately, few if any historical databases were obtained in this way. However, close approximations might be made using interpolation methods, which possibly could be more sophisticated than the linear interpolation techniques often used by foresters in the past. Unbiased or nearly-unbiased estimates obtained by design-based stereological estimators might have appeal for tree volume and weight estimation. Frequently, hundreds or thousands of individual tree volume estimates may be made in large forest inventories. In this situation, the modest individual tree volume variance of stereological estimates might be tolerable if the reduction of estimator bias is obtained compared to the usual model-based regression estimates, because bias cannot be reduced even by very large sample sizes.

3. Results and Discussion

3.1. Systematic Random-Start Sampling for Tree Bole Volume

A design-unbiased method of estimating tree stem volume by systematic random start sampling was proposed above in the form of Estimator (7). Similar model-based systematic estimators have been used for many years but do not to our knowledge employ a random starting point. We tested this estimator using several standard tree stem shapes and two tree sizes using simulation methods. The results are indicated in

Table 1. Systematic-random start bole volume estimation with cross-sections $h=2$ m apart and rotated power function having power of height $r/2$, tree bole basal diameter $=37.5$ cm, length $=30$ m, based on 50,493 random trials.

Shape	r	Volume $\left({\mathbf{m}}^3\right)$	Std. Error $\left({\mathbf{m}}^3\right)$	RMSE% ${}^1$
Paraboloid	1	8.84	0.34	3.85%
Paracone	$\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.$	7.57	0.34	4.49%
Cone	2	5.89	0.34	5.78%
Neiloid	3	4.42	0.34	7.68%

${}^1$ RMSE is Root Mean Square Error.

below.

These results indicate that the standard error of the method is the same for all the tree shapes considered. We chose these shapes because a tree stem is sometimes considered to be similar to the frustum of a neiloid in its base, similar to the frustim of a paraboloid in its midsection, and similar to a cone in the upper section leading to the tip of the tree [13] (p. 100). The stem shape corresponding to $r=4/3$ is Forslund’s paracone [17], which is midway between a cone and a paraboloid and proposed as a simple shape that can be used to approximate tree stem cubic volume. The RMSE% values increase from 3.85% to 7.68% as stem content decreases for the paraboloid, paracone, cone and neiloid respectively. This is simply due to the decrease in cubic meter stem content associated with the respective shapes, and because the height and diameter of each shape was the same, and the convexity of the stem profile decreases for these shapes in the order given.

We also tested the systematic random-start Estimator (8) where cross-sectional area is selected according to a constant decrease of ${\Delta }_c=0.04$ square meters with the result that the linear distance between cross-sections can be variable. For purposes of comparison, the same size tree bole was used in these simulations as in the simulations of Table 1. The results of these simulations are given in

Table 2. Systematic-random start bole volume estimation with cross-sections decreasing by steps of ${\Delta }_c=0.04$ square meters consequently with variable linear distances apart and rotated power function having power of height $r/2$, tree bole basal diameter $=37.5$ cm, length $=30$ m, based on 10,493 random trials.

Shape	r	Volume $\left({\mathbf{m}}^3\right)$	Std. Error $\left({\mathbf{m}}^3\right)$	RMSE% ${}^1$
Paraboloid	1	8.84	0.34	3.84%
Paracone	$\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.$	7.57	0.32	4.28%
Cone	2	5.89	0.31	5.24%
Neiloid	3	4.42	0.29	6.57%

${}^1$ RMSE is Root Mean Square Error.

below.

It should not be surprising that the results in Table 2 for the paraboloid bole shape are the same as those for Table 2, because for a paraboloid height is linearly related to cross-sectional area, so the methods turn out to be equivalent in that case. In a similar way, the height accumulation method of Grosenbaugh [18] is essentially equivalent to selecting cross-sections at equally spaced bole heights if the tree bole is shaped exactly like a simple cone. The slight difference in RMSE% is due to decimal places not reported for standard error, which in turn are likely due to variation in simulation results. Of more interest are the results for other non-parabolic shapes. Table 2 indicates that standard errors and RMSE% for the paracone, cone, and neiloid are smaller than those associated with Table 1, but the differences are slight.

These results do not suggest that the method of systematic random-start cross-sections with constant decrease in cross-sectional area ${\Delta }_c$ is substantially better from a practical point of view than the method of systematic random-start cross-sections with a constant stem interval h, however it is possible that the differences could be somewhat greater on actual tree stems which have more complex bole shapes. The standard errors in Table 2 associated with the method of constant cross-sectional area decrease as bole shape departs more substantially from parabolic and are slightly better than those of Table 1 for the nonparbolic shapes. The percentage improvement for the neiloidal shape is about 15%, improvement for the cone, 9% and for the paracone 6%. Of these, the paracone is likely the closest of these very simple shapes to a typical tree bole shape. These results do appear to support a slight improvement in precision from taking measurements at fixed intervals of cross-sectional area decrease ${\Delta }_c$.

3.2. Potential for Application of Stereological Variance Estimators

Systematic sampling has been extensively used in stereology as well as in forest sampling. Many systematic samples in forestry have used variance formulas developed under random sampling assumptions [6] (pp. 60–61). An early investigation of the problem of variance estimation in forestry applications has been provided by Matérn [67]. Variance reductions on the order of $O\left(n^{-1}\right)$ have been indicated for systematic sampling in some stereological applications [3] (pp. 302–314). Chapter 13 of Baddeley and Jensen [3] (pp. 295–328) discusses the variance of stereological estimators. Several approaches to quantifying the variance of stereological estimators involve modeling the covariogram because of the covariance among systematic samples. Baddeley and Jensen [3] (p. 329) indicate that most past approaches have either modeled covariances directly or have modeled the spatial structure of the population and then derived model covariances. Gunderson et al. [68] discuss estimation of the variance associated with systematic sampling in stereology. They discuss the case of estimating a parameter that can be expressed as a 1-dimensional integral. An example of such a parameter in forestry applications would be individual tree stem volume as discussed in examples above. In some situations, line samples emanating from a baseline might be characterized as samples from an integral over the baseline, for example, the proposal of Lynch [27] for horizontal line sampling to sample riparian areas.

According to Baddeley and Jensen [3] (p. 328), stereological approaches to variance estimation have been informed by previous work in geostatistics by Matheron [69]. There have been several applications of geostatistics to forest sampling and inventory. Köhl and Gertner [70] discussed the application of geostatistics to forest damage surveys as well as describing several previous applications of geostatistics to forest sampling.

It may be that improved variance estimation for some forestry sampling designs, especially those which employ systematic sampling, could eventually result from the study of variance estimation techniques for certain stereological estimators. Precise field locations for sample points, plots, and lines are now commonly available from GPS technology, perhaps facilitating the computation of spatial covariances. Most forest inventory estimates are currently being provided by using specialized software, making the use of complex variance estimators potentially more feasible.

Barabesi and Marcheselli [71] obtained significant variance reductions on the order of $O\left(n^{-1}\right)$ for systematic non-aligned line-intercept sampling. Lynch [2] indicated that systematic nonaligned sampling could be applied to radial line sampling by dividing the interval $[0,2\pi )$ into n sub-intervals with a radial line sample randomly located within each sub-interval. Systematic non-aligned sampling might be useful for other forest sampling applications. Currently, GIS software often allows one to assign sampling locations randomly with restrictions on the minimum distances between sample locations. Foresters then navigate to these sample locations using GPS. This type of procedure could be extended to perform systematic non-aligned sampling for inventories of forested areas. It may be that the variance estimators proposed by Barabesi and Marcheselli [71] could be adapted for this situation.

3.3. Monte Carlo Integration and Stereology

We have noted above that several applications of Monte Carlo integration for estimation of tree stem volume could be interpreted in stereological terms. However, we have not found recognition of the relationship between Monte Carlo integration and stereological methods in either the literature of Monte Carlo integration or stereolgy (although, we do not claim to have exhaustively studied the vast literature in these areas). Probably, this is due to the fact that most applications of Monte Carlo integration outside the literature of forest mensuration have not been based on measurements of physical dimensions of objects, but rather methods to solve intractable definite integrals. Conversely, most stereological applications are based on physical measurements and do not attempt to solve intractable integrals.

As indicated above, foresters have found that variance reduction strategies used in Monte Carlo integration including importance sampling, control variates and antithetic variates are advantageous for estimation of tree stem volume. This is largely due to the fact that foresters have a lot of prior knowledge concerning the shapes of tree stems. This may not as frequently be the case for physical objects commonly encountered in stereological problems, such as the internal organs of mammalian species. As well, importance sampling may be considered related to estimation using inverse selection probabilities which has been applied in stereology. Likewise, control variates is essentially similar to regression estimation as applied in survey sampling.

4. Conclusions

Several forest sampling techniques are similar to methods that have been used in stereology. Many forest sampling methods can be viewed as samples of 2-dimensional or 3-dimensional quantities by using lower dimensional samples such as 2-dimensional tree stem cross-sections, 1-dimensional lines or counts of 0-dimensional points. Basal area per hectare, a 2-dimensional quantity, can be estimated using counts of 0-dimensional points by the Bitterlich point sampling technique. Also, systematic sampling is common in stereology and in many forest sampling designs.

Many forest estimation problems such as estimation of tree stem content from samples of tree cross-sections more closely resemble the stereology paradigm than the classical survey statistics paradigm in which sample units are obtained from a finite population of units. Thus, development of forest sampling theory and practice might be positively impacted by the study of the large literature of stereology.