Evaluation of Image-Assisted Forest Monitoring: A Simulation

Abstract

Fiscal uncertainties can sometimes affect national continuous forest monitoring efforts. One solution of interest is to lengthen the time it takes to collect a “full set” of plot data from five to 10 years in order to reduce costs. Here, we investigate using ancillary information to partially offset this proposed solution’s negative effects. We focus our discussion on the corresponding number of years between measurements of each plot while we investigate how thoroughly the detrimental effects of the reduced sampling effort can be ameliorated with change estimates obtained from temporally-dense remotely-sensed images. We simulate measured plot data under four sampling error structures, and we simulate remotely-sensed change estimates under three reliability assumptions, integrated with assumptions about the additional unobserved growth resulting from the lengthened observation window. We investigate a number of estimation systems with respect to their ability to provide compatible annual estimates of the components of change during years spanned by at least half of the full set of plot observations. We show that auxiliary data with shorter observation intervals can contribute to a significant improvement in estimation.

1. Introduction

Fiscal uncertainties can sometimes affect national continuous forest monitoring efforts. The sample designs for these large long-term efforts are often inflexible and can be adversely affected by budget shortfalls. When a budget shortfall affects a national forest inventory (NFI) program to the extent that it reduces the interpretive ability of the effort, the effects can be wide-ranging and profound. For instance, in Europe, there has been a concerted effort to interpret multiple national forest inventories, in unison, in order to understand continent-wide (and world-wide) forest population change by harmonizing monitoring information at the supra-national level (e.g., Tomppo and Schadauer [1]). This harmonization process, as well as monitoring forest population change in general, can be hampered severely by national budget shortfalls, which are occurring with increasing frequency. A further complication lies in the fact that limitations of capital expenses can be very different in different countries. This has led to many different solutions for the problem. In some cases, the scope and objectives of the monitoring have been limited by strategies involving a reduction of the number of sample plots or an increase in the number of years between measurements of the sample plots. Decreasing the number of the plots might be effected randomly over the entire area covered by a particular monitoring effort or it might be restricted to areas considered to be of minor economic or managerial significance. Reducing the scope or objectives of the monitoring effort can result in the collection of a reduced set of variables or in fewer measurements of a particular variable. For instance, because the measurement of tree height can be very time consuming, limiting the number of tree height measurements is often one of the effected solutions. Increasing the number of years between measurements of sample plots can involve all of the plots or only the plots in particular subpopulations. In addition, a staggered measurement schedule might be incorporated, such as measuring 50% of the plots in one measurement year and the remaining 50% of the plots in another one. By this scheme, the number of years between measurements of each set of plots might remain unchanged. An extension of this idea, which was intended to result in annually-consistent budgets, led to the United States Department of Agriculture (USDA) Forest Service’s panelized annual inventory system conducted by the Forest Inventory and Analysis (FIA) Program. Even with this level of budgetary planning, further shortfalls can occur. Here, we investigate the effects of one proposed solution to a potential budget shortfall for the national monitoring effort conducted by FIA. The solution of interest is to lengthen the “cycle”, that is the time it takes to collect a “full set” of plot data from five years to 10 years. This solution can have manifold implications, a few of which have been discussed in Roesch [2] and Roesch and Van Deusen [3].

We use a simulation to investigate how thoroughly the detrimental effects of the reduced sampling effort can be ameliorated with change estimates obtained from temporally-dense remotely-sensed images. Within the simulation, we use four sampling error structures for the measured plot data and three assumptions of applicability of the remotely-sensed image change estimates (ICE).

Specifically, our objective is to obtain annual estimates of the components of change, given multi-year data collected by temporally-overlapping panels. Van Deusen and Roesch [4] and Roesch and Van Deusen [5] explored estimation of the change in forest land classification, in the context of the national forest inventory in the United States, while Roesch [6] concentrated on estimation of the change in the tree population on land that remains forested throughout the period of interest. The current work focuses on land that is forested at some point during a period of interest. Eriksson [7] presented a set of definitions for the continuous components of change. In this paper, we subscribe to the definitions of Roesch [8], which presented a discrete version of the Eriksson definitions. We define entry as the cubic meter volume (or value) of trees as they attain the entry criterion; live growth as the annual growth in volume that occurs on trees after a defined entry criterion has been achieved; mortality as the volume of trees as they die; and harvest as the volume of trees as they are harvested.

These population components of change are compatible, that is:

$$Y_{t+1}=Y_t+L_t+E_t-M_t-H_t$$

(1)

where:

• Y_t = the value of interest at the beginning of year t;
• L_t = growth in the value of interest on live trees during year t;
• E_t = the value of interest on live trees as they enter the population during year t;
• M_t = the value of interest on trees as they die during year t and;
• H_t = the value of interest on trees as they are harvested during year t.

It has often been argued that because these components are compatible, the estimators of these components should also be compatible. The recursive estimation systems that we explore here are compatible and are compared to non-recursive estimation systems. For clarity and without loss of generality, assume that the estimands of interest are the cubic meter volume per hectare of all live trees in a fixed area in each component of change category during each year of a multi-decadal period. The plot sampling design considered here consists of a number (g) of mutually-exclusive, spatially-disjoint and temporally-successive panels. The first panel to be measured is selected by locating a random point within the population area. One panel per year is then measured, in a predetermined order of succession, for g years. After the entire set of panels has been measured, the sequence reinitiates. Bechtold and Patterson [9] discuss an example of this type of design. Various statistical aspects of this type of design have been discussed by McCollum [10], Van Deusen [11], Roesch et al. [12] and Van Deusen [13], among others. This panel design provides remeasurement observations that are spatially disjoint and temporally overlapping. Usually, analysts are interested in a temporal scale that is finer than the scale of observation, which has resulted in a number of suggestions for obtaining those finer-scale estimates. Roesch [8] argued that the average annual growth within each individual panel is best applied to the center of the measurement interval, which is analogous to an assumption of linear change between observations. This was thought to be a reasonable first approximation, in lieu of contradictory evidence. All analytical methods proposed to date for this class of forest monitoring sample designs have been predicated on this or similar assumptions. For the interested reader, a related discussion can be found in Westfall et al. [14].

Roesch and Van Deusen [3] discussed the effects of ignoring differences in the temporal aspects of a realized sample from the intended design. Specifically, they showed the effects of two (usually tacit) assumptions in NFI designs. The first was that variation in the time of observation for an individual areal sample is ignorable. The second was that variation in the remeasurement period lengths between individual plots in successive areal samples was ignorable. They explored the effects of these assumptions and discussed how inference can be improved by a judicious accounting of these sampling disparities. They showed the remeasurement period assumption to be especially problematic. That is, plots in NFI systems are never remeasured on exact temporal intervals, and large biases can be introduced when there is little effort made to restrict the distribution of temporal interval lengths. They concluded that further research was needed to determine what restrictions should be placed on the distribution of temporal intervals to achieve specific objectives. This paper constitutes one specific extension of that research in that we explore some of the effects of an unusually-long survey cycle.

2. Methods

2.1. Simulated Population

To construct the simulated population, we used methods similar to those used to construct one of the populations in Roesch [6], specifically Population 1 in the referenced work. We started with the FIA plot data measured at least twice by the annual sample design in the state of Georgia (USA) between 1995 and 2012. This resulted in 7330 plots, most of which had 3 measurement times (i.e., two observed growth intervals for each component). Note that although the FIA sample design includes a plot remeasurement interval of 5 years, logistical adjustments, especially in the early years of usage of the design, led to the actual remeasurement interval varying quite widely around the 5-year target for these plots. To create a plausible population, we first created a seed population (Population 0) utilizing the theory and methods described in Roesch [6] for deriving annual values from the multi-year observations. The seed population allowed us to simulate a population (Population 1), from which we might assume the observed sample data could have been drawn. Because exact harvest times were unknown, harvested volume was randomly allocated to a year within each observation interval. Linear interpolation and extrapolation were used to obtain an initial value for the live growth, entry and mortality change components for each year from 1995 through 2012, as well as a starting cubic-meter per hectare value in the beginning of 1995, with temporal adjustments made as necessary for high levels of harvest and mortality (note that, as in Roesch [6], we are not attempting to reconstruct the plot, but rather a reasonable facsimile to the forested condition from which it could have been drawn). Construction of Population 0 then proceeded with 500 variance-interjected copies of Set 1, resulting in 3,665,000 hectares, 2,360,411 of which were forested at some time during the period of interest. The population is represented by 5 matrices, each with 3,665,000 rows and 18 columns (one column for each year), one matrix for each change component and one for initial annual volume. Variance was interjected at two levels. In Level 1, to maintain the trend while adding variance to the seed, all values for each component on each hectare were multiplied by a unique random variate, drawn from an N (1, 0.025) distribution. The second level of variance was introduced temporally by multiplying the result of Step 1 for each annual value for each component on each hectare by a unique random variate drawn from an N (1, 0.0025) distribution, thereby completing the construction of Population 0. For Population 1, a mild (latent) non-linear trend was introduced into each of the components of Population 0, $T{1}_i=\left[0.95+\left(0.05ln\left(i-1997\right)\right)\right], i=1998 \text{to} 2011$, where each value in each year i is multiplied by T1_i. Table 1 gives the Population 1 distribution statistics for 1998 to 2011. Population 1 is available from the first author upon request.

2.2. Sampling Simulations

Each simulation consisted of 1000 iterations of 1000 plots each (without replacement) from each population, under each of the four sampling error structures. In the current investigation, we use a method similar to that described in Roesch [6], to consider 4 sampling error structures. The sampling errors used here are elevated relative to those described in Roesch [6] in order to account for additional errors resulting from the longer cycle length or observation window. Note that the use of the sampling error structures to account for unobserved live growth should be greater for a 10-year observation window than for a 5-year observation window. This affects all live growth on mortality and cut trees between the last observation and the time of death or harvest. Therefore, the components of live growth, mortality, and harvest are more heavily affected by this additional error than entry. The entry component is affected by the cohort of trees that entered the population and died or were harvested prior to the next observation.

Each of the four sampling error structures was introduced by multiplying a unique random normal deviate of mean b and standard deviation d from Table 2 by each sampled observation of initial volume, entry, live growth, mortality and harvest.

Table 1. Distribution statistics for Population 1.

Component	Statistic	Year
		1998	1999	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011
Volume	Minimum	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	1st Quartile	0.00	0.00	0.42	1.89	4.33	7.11	9.86	12.56	14.84	16.47	17.58	18.59	19.17	19.59
	Median	9.91	14.93	19.98	25.31	30.71	35.84	41.21	45.78	49.51	53.05	56.53	59.07	61.52	63.46
	Mean	50.53	54.21	57.79	61.45	65.60	69.50	73.21	76.16	78.78	81.30	83.62	85.70	87.46	89.00
	3rd Quartile	77.83	83.76	89.56	95.08	101.26	106.63	110.80	113.99	116.78	120.38	123.21	126.65	130.25	132.89
	Maximum	813.34	813.36	814.69	815.08	816.37	818.58	821.72	825.90	831.01	837.24	844.48	852.75	862.01	872.21
Live Growth	Minimum	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	1st Quartile	0.00	0.00	0.00	0.00	0.29	0.66	0.73	0.66	0.50	0.22	0.00	0.00	0.00	0.00
	Median	1.19	1.49	2.02	2.49	3.04	3.36	3.43	3.40	3.36	3.20	2.35	0.73	0.00	0.00
	Mean	3.57	3.59	3.79	4.01	4.27	4.43	4.49	4.47	4.48	4.53	4.13	3.51	2.75	1.71
	3rd Quartile	5.30	5.39	5.70	6.03	6.35	6.50	6.65	6.68	6.72	6.94	6.56	5.65	3.94	0.00
	Maximum	52.52	50.58	66.35	58.15	48.79	39.41	33.73	35.04	39.47	43.90	44.57	41.38	39.08	46.09
Entry	Minimum	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	1st Quartile	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	Median	0.00	0.04	0.11	0.17	0.22	0.25	0.26	0.24	0.22	0.18	0.08	0.00	0.00	0.00
	Mean	0.45	0.47	0.52	0.57	0.62	0.64	0.65	0.64	0.63	0.64	0.57	0.47	0.35	0.22
	3rd Quartile	0.46	0.50	0.58	0.64	0.72	0.76	0.78	0.75	0.72	0.67	0.56	0.38	0.16	0.00
	Maximum	21.21	21.60	29.03	27.63	17.04	12.20	9.91	11.56	13.28	15.23	16.71	15.15	15.01	14.90
Mortality	Minimum	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	1st Quartile	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	Median	0.00	0.00	0.00	0.00	0.00	0.03	0.06	0.06	0.02	0.00	0.00	0.00	0.00	0.00
	Mean	0.74	0.74	0.76	0.78	0.82	0.86	0.88	0.91	0.94	0.97	0.90	0.77	0.60	0.38
	3rd Quartile	0.21	0.28	0.39	0.46	0.55	0.65	0.71	0.73	0.72	0.69	0.48	0.20	0.00	0.00
	Maximum	108.46	93.27	78.00	77.97	78.57	78.77	79.11	79.44	79.54	80.31	80.39	80.51	80.83	80.94
Harvest	Minimum	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	1st Quartile	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	Median	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	Mean	0.63	0.75	0.84	0.99	1.21	1.45	1.78	1.99	2.16	2.20	1.90	1.55	1.05	0.55
	3rd Quartile	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	Maximum	341.31	373.88	352.31	358.61	447.97	466.30	466.29	462.06	465.60	483.98	485.19	498.90	488.34	510.38

Table 1. Distribution statistics for Population 1.

Component	Statistic	Year
		1998	1999	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011
Volume	Minimum	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	1st Quartile	0.00	0.00	0.42	1.89	4.33	7.11	9.86	12.56	14.84	16.47	17.58	18.59	19.17	19.59
	Median	9.91	14.93	19.98	25.31	30.71	35.84	41.21	45.78	49.51	53.05	56.53	59.07	61.52	63.46
	Mean	50.53	54.21	57.79	61.45	65.60	69.50	73.21	76.16	78.78	81.30	83.62	85.70	87.46	89.00
	3rd Quartile	77.83	83.76	89.56	95.08	101.26	106.63	110.80	113.99	116.78	120.38	123.21	126.65	130.25	132.89
	Maximum	813.34	813.36	814.69	815.08	816.37	818.58	821.72	825.90	831.01	837.24	844.48	852.75	862.01	872.21
Live Growth	Minimum	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	1st Quartile	0.00	0.00	0.00	0.00	0.29	0.66	0.73	0.66	0.50	0.22	0.00	0.00	0.00	0.00
	Median	1.19	1.49	2.02	2.49	3.04	3.36	3.43	3.40	3.36	3.20	2.35	0.73	0.00	0.00
	Mean	3.57	3.59	3.79	4.01	4.27	4.43	4.49	4.47	4.48	4.53	4.13	3.51	2.75	1.71
	3rd Quartile	5.30	5.39	5.70	6.03	6.35	6.50	6.65	6.68	6.72	6.94	6.56	5.65	3.94	0.00
	Maximum	52.52	50.58	66.35	58.15	48.79	39.41	33.73	35.04	39.47	43.90	44.57	41.38	39.08	46.09
Entry	Minimum	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	1st Quartile	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	Median	0.00	0.04	0.11	0.17	0.22	0.25	0.26	0.24	0.22	0.18	0.08	0.00	0.00	0.00
	Mean	0.45	0.47	0.52	0.57	0.62	0.64	0.65	0.64	0.63	0.64	0.57	0.47	0.35	0.22
	3rd Quartile	0.46	0.50	0.58	0.64	0.72	0.76	0.78	0.75	0.72	0.67	0.56	0.38	0.16	0.00
	Maximum	21.21	21.60	29.03	27.63	17.04	12.20	9.91	11.56	13.28	15.23	16.71	15.15	15.01	14.90
Mortality	Minimum	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	1st Quartile	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	Median	0.00	0.00	0.00	0.00	0.00	0.03	0.06	0.06	0.02	0.00	0.00	0.00	0.00	0.00
	Mean	0.74	0.74	0.76	0.78	0.82	0.86	0.88	0.91	0.94	0.97	0.90	0.77	0.60	0.38
	3rd Quartile	0.21	0.28	0.39	0.46	0.55	0.65	0.71	0.73	0.72	0.69	0.48	0.20	0.00	0.00
	Maximum	108.46	93.27	78.00	77.97	78.57	78.77	79.11	79.44	79.54	80.31	80.39	80.51	80.83	80.94
Harvest	Minimum	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	1st Quartile	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	Median	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	Mean	0.63	0.75	0.84	0.99	1.21	1.45	1.78	1.99	2.16	2.20	1.90	1.55	1.05	0.55
	3rd Quartile	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	Maximum	341.31	373.88	352.31	358.61	447.97	466.30	466.29	462.06	465.60	483.98	485.19	498.90	488.34	510.38

Table 2. Mean (b) and standard deviation (d) of the normal random variates used to simulate each sample observation.

Component	Error Structure
	1		2		3		4
	b	d	b	d	b	d	b	d
Initial Volume	1.01	0.10	1.00	0.03	1.00	0.05	0.99	0.10
Entry	1.01	0.10	1.00	0.03	1.00	0.05	0.99	0.10
Live Growth	1.01	0.10	0.99	0.03	0.98	0.05	0.95	0.10
Mortality	1.01	0.10	0.99	0.03	0.98	0.05	0.95	0.10
Harvest	1.01	0.10	0.99	0.03	0.98	0.05	0.95	0.10

Table 2. Mean (b) and standard deviation (d) of the normal random variates used to simulate each sample observation.

Component	Error Structure
	1		2		3		4
	b	d	b	d	b	d	b	d
Initial Volume	1.01	0.10	1.00	0.03	1.00	0.05	0.99	0.10
Entry	1.01	0.10	1.00	0.03	1.00	0.05	0.99	0.10
Live Growth	1.01	0.10	0.99	0.03	0.98	0.05	0.95	0.10
Mortality	1.01	0.10	0.99	0.03	0.98	0.05	0.95	0.10
Harvest	1.01	0.10	0.99	0.03	0.98	0.05	0.95	0.10

2.3. Moving-Window Mean of Ratios Estimator

Our estimation approach uses the (semi-centralized) moving-window mean of ratios estimator (MWMOR) in Roesch and Van Deusen [3] and Roesch [6], an estimator that arose from a different perspective than previously-developed estimators. The idea was simple. One stacks the observations on a temporal scale (or a function of the temporal scale) and then slices through the stack (say to create annual segments) to determine what proportion of each observation contributes to the estimate for each year spanned by the observation. For this problem, as in Roesch and Van Deusen [3], we use the general three-dimensional selection model given in Roesch [2] with the exception that time will be re-scaled relative to the proportion of the growing season elapsed within each year. Assign to each observation of variable x labels for plot i and a superscript representing the beginning value and ending value as x_i^b and x_i^e, respectively. Because there are no observations between x_i^b and x_i^e, the distribution of the volume growth between the two observations must be modeled. In Roesch and Van Deusen [3], two simplifying assumptions were made. In the first, it was assumed that the growing season begins and ends on the same dates for each year in the area of interest. In the second, it was assumed that growth for each plot is uniform throughout the growing season. Given these assumptions, we can temporally order each observation by the year of observation plus the proportion (p) of the growing season that has elapsed. We use s_i to represent the difference in these values to obtain the growing season-adjusted temporal span between the beginning and ending observations. We then allocate the proportion of growth observed over s_i to the proportion of each year spanned by s_i, (thereby accounting for the marginal probability of the time dimension). Modeling growth between observations allows us to allocate growth within components to the years the growth occurred. A simple time-adjusted estimator for annual volume growth (within the growth component) is the moving-window mean of ratios estimator (MWMOR) for component C in year t:

$$C_t^M=\frac{1}{n_t}{\displaystyle \sum }_{i=1}^{n_t}\frac{c_{i,t}}{p_{i,t}}$$

(2)

where:

• $n_t=$ the number of plots observing growth in year t;
• $p_{i,t}=$ the product of portion of year t growing season observed by plot i and the portion of plot i area within the area of interest and;
• $c_{i,t}=$ the value of component C observed on plot i, assignable to year t.

The general statistical properties of ratio estimators are well known and can be found in very early works, such as Raj [15], Walton and DeMars [16], Cassel et al. [17] and Cochran [18].

In some of the estimators below, we use a trend distributed version of $C_t^M$ in which $c_{i,t}$ in Equation (2) is replaced by $c_{i,t}^T$. The value of $c_{i,t}^T$ is determined by the value of the trend in initial cubic meter volume at year t. A proportion is allocated to each year observed in a plot remeasurement interval, $c_{i,t}^T={\pi }_y{x}_i$, where x_i is the total value of a change component (such as live growth) observed over an interval of years y = 1 to Y, and ${\pi }_y=v_y/{\displaystyle \sum }_{j=1}^Y{v}_j$, where v_j is the sample mean initial volume for year j. Note that we use the variable t to refer to the general year within the population or estimation interval and the variable y to refer to the year within a particular sample plot’s remeasurement interval. Because plots are measured in panels, y selects a panel-specific subset of t. We use this change of variable below, where it is appropriate. The underlying assumption of $c_{i,t}^T$ is that the level of a change component during a particular year will be proportional to the standing volume at the beginning of the year. The trend distributed estimators then take the form:

$$C_t^T=\frac{1}{n_t}{\displaystyle \sum }_{i=1}^{n_t}\frac{c_{i,t}^T}{p_{i,t}}$$

(3)

2.4. Incorporating Image Change Estimates

Roesch [6] demonstrated the utility of a weighting method that can be used in conjunction with the MWMOR estimators above. We use a variant of that approach here in order to use the remotely-sensed information to better approximate the timing of harvests and the corresponding allocation of the components of change within a remeasurement window.

We simulated the incorporation of ICE observations from remotely-sensed images every two years under three different assumptions of reliability. In Assumption 1 (A1), a clearcut, defined as the harvest of at least 95% of the cubic-meter volume, is successfully identified in the image 95% of the time, and there are no false positives. Under this assumption, when a clearcut is identified, it is known to have occurred since the last plot visit, but may have been missed during intervening ICE observations. Under Assumption 2 (A2), all clearcut harvests are correctly identified 100% of the time and assigned to the correct year. Again, there are no false positives. Assumption 3 (A3) reflects the highest level of reliability in which every harvest (whether a partial harvest or a full harvest) is identified and assigned to the correct year; we can estimate the proportion harvested to within ±5% (truncated at 0% and 100%), and there are no false positives.

Under all three of these reliability assumptions, ICE data collected every two years will be in sync with about half of the plot observations collected under the panel design and out of sync with the other half. That is, one panel or 1/10 of the plots is observed every year, and the ICE observations are made every year. For half of the plots, the first ICE observation made after the plot observation will be in one year, and for the other half of the plots, the first ICE observation after the plot observation will be in two years.

There are a number of different approaches that could be taken to determine the weights used to allocate the change components to the intervening years, depending on the analyst’s prior knowledge with respect to the reliability of the ICE data and the availability of reasonable growth models. For instance, by the clearcut definition above, up to 5% of the cubic-meter volume could remain on live trees, and new trees could develop and grow subsequent to the harvest and prior to the next sample plot observation. For partial harvests, the possibility for subsequent growth is much higher. In large countries with a very diverse range of species and forest types, there are often not reasonable growth models available for many of those types. Additionally, there will usually be unobserved growth on trees that die or are harvested during the observation interval. In the case of a clearcut, this unobserved growth occurred on the prior stand, so it is reasonable to assume that it will usually far outweigh the observed growth, because growth is observed only on trees present at the next plot observation. Therefore, under Assumptions 1 and 2, we simplify the annual allocation problem by assuming that all observed growth occurred before the ICE clearcut observation. Under Assumption 3, for all harvests, we use an allocation that separates pre-harvest and post-harvest change based on an estimate of the proportion of volume harvested.

2.4.1. Estimation under A1

Under Assumption A1, when a clearcut is identified, it is known to have occurred since the last plot observation, but may have been missed during intervening ICE observations. To establish weights, we could either use an estimate of the probability of having missed the clearcut at previous ICE observations or we could assume that the probability is unknown. Because A1 is intended to be our worst case assumption in the simulation, we assumed that this probability was unknown. The allocation of partial harvest volumes is unaffected by ICE observations under this assumption.

We define $y_{i,I}$ as the number of years of the first ICE observation of clearcut for plot i since the previous plot measurement. Although it would be a rare event, we also define $y_{i,II}\ge y_{i,II}+2$ as the year of the second ICE observation of clearcut for plot i during the remeasurement interval, in order to estimate the weight for harvest. Additionally, when there are two ICE-identified clearcuts on plot i, we define r_i,I as the proportion of total clearcut volume in the first clearcut and r_i,II = 1 − r_i,I as the proportion of total clearcut volume in the second clearcut.

$${\pi }_{i,y}^{ThA1}=\{\begin{array}{cc}{\pi }_y& \text{if no ICE-identified clearcut}\hfill \\ {y_{i,I}}^{-1}& \text{if one ICE-identified clearcut and }y\le y_{i,I}\hfill \\ r_{i,I}& \text{if ICE is out-of-sync with two identified clearcuts and }y=y_{i,I}\hfill \\ .5r_{i,I}& \text{if ICE is in-sync with two identified clearcuts and }y=y_{i,I}-1\hfill \\ .5r_{i,I}& \text{if ICE is in-sync with two identified clearcuts and }y=y_{i,I}\hfill \\ .5r_{i,II}& \text{if two ICE-identified clearcuts and }y=y_{i,II}-1\hfill \\ .5r_{i,II}& \text{if two ICE-identified clearcuts and }y=y_{i,II}\hfill \\ 0& \text{otherwise}\hfill \end{array}$$

(4)

The proportion r_i,I would have to be estimated. An estimate might be obtained from the ICE data, but it is not directly available from the sample plot data. Under A1, we do not assume that we can confidently estimate r_i,I from the images, so we set it equal to 0.5. The allocation for the harvest component on plot i for year t under A1 is:

$$h_{i,t}^{ThA1}={\pi }_{i,y}^{ThA1}{h}_i$$

(5)

where h_i is the observation of cubic meter harvested volume on plot i during the remeasurement interval. The estimator for the harvest component for year t under A1 is:

$$C_t^{hA1}=\frac{1}{n_t}{\displaystyle \sum }_{i=1}^{n_t}\frac{h_{i,t}^{ThA1}}{p_{i,t}}$$

(6)

For the other components of change (live growth, entry and mortality), the allocation to each year between plot observation and ICE clearcut call:

$${\pi }_{i,y}^{ToA1}=\{\begin{array}{c}\begin{array}{cc}{\pi }_y& \text{if no ICE-identified clearcut }\\ {y_{i,I}}^{-1}& \text{if one ICE-identified clearcut and }y\le y_{i,I} \end{array}\\ \begin{array}{cc}0 & \text{if no ICE-identified clearcut and }y>y_{i,I }\\ {\pi }_y & \text{if two or more ICE-identified clearcuts }\end{array}\end{array}$$

(7)

The allocation for the other components of change on plot i for year t under A1 is represented as $o_{i,t}^{ToA1}={\pi }_{i,y}^{ToA1}{o}_i$, where o_i is the observation of a specific component o for plot i during the interval. The estimator for the other components of change for year t under A1 is:

$$C_t^{oA1}=\frac{1}{n_t}{\displaystyle \sum }_{i=1}^{n_t}\frac{o_{i,t}^{ToA1}}{p_{i,t}}$$

(8)

2.4.2. Estimation under A2

Under Assumption 2 (A2), we know if a clear-cut has occurred within the past two years, and we can assign it to the correct year. As under A1, the allocation of partial harvest volumes is unaffected by ICE observations. The weight for the harvest component is defined as:

$${\pi }_{i,y}^{ThA2}=\{\begin{array}{cc}{\pi }_y& \text{if no ICE-identified clearcut}\hfill \\ 1& \text{if one ICE-identified clearcut and }y=y_{i,I}\hfill \\ r_{i,I}& \text{if two ICE-identified clearcuts and }y=y_{i,I}\hfill \\ r_{i,II}& \text{if two ICE-identified clearcuts and }y=y_{i,II}\hfill \\ 0& \text{otherwise}\hfill \end{array}$$

(9)

As we mentioned above, the proportion r_i,I would have to be estimated, and it is not directly available from the sample plot data. However, by definition, a clearcut consists of the harvest of 95% of the volume present during the year of harvest, so r_i,I could be approximately estimated by:

$${\widehat{r}}_{i,I}=\frac{v_I}{v_I+v_{II}}$$

(10)

where $v_I$ is the sample mean initial volume for year for the year of the first ICE-identified clearcut and $v_{II}$ is the sample mean initial volume for the year of the second ICE-identified clearcut. An estimator for $r_{i,II}$ would then be ${\widehat{r}}_{i,II}=1-{\widehat{r}}_{i,I}$.

The allocation for the harvest component on plot i for year t under A2 is:

$$h_{i,t}^{ThA2}={\pi }_{i,y}^{ThA2}{h}_i$$

(11)

This leads to the estimator for the harvest component for year t under A2:

$$C_t^{hA2}=\frac{1}{n_t}{\displaystyle \sum }_{i=1}^{n_t}\frac{h_{i,t}^{ThA2}}{p_{i,t}}$$

(12)

For the other components of change, if there were a single clearcut, we assign all observed growth to the years prior to the ICE clearcut observation, but otherwise, we use proportional allocation as described above. We define the value of the trend by a proportion that is allocated to each year y in an interval beginning with the year of the sample plot observation for plot i and ending with the year of the first ICE observation of clearcut ($y_{i,I}$):

$${\pi }_{i,y}^I=\raisebox{1ex}{$v_y$}\!\left/ \!\raisebox{-1ex}{${{\displaystyle \sum }}_{j=1}^{y_I}{v}_j$}\right.$$

(13)

Then, the weight for the other components is:

$${\pi }_{i,y}^{ToA2}=\{\begin{array}{cc}{\pi }_y& \text{if no ICE-identified clearcut}\hfill \\ {\pi }_{i,y}^I& \text{if one ICE-identified clearcut and }y\le y_{i,I}\hfill \\ {\pi }_y& \text{if two ICE-identified clearcuts}\hfill \\ 0& \text{otherwise}\hfill \end{array}$$

(14)

The allocation for the other components on plot i for year t under A2 is then:

$$o_{i,t}^{ToA2}={\pi }_{i,y}^{ToA2}{o}_i$$

(15)

The estimator for the other components of change for year t under A2 is:

$$C_t^{oA2}=\frac{1}{n_t}{\displaystyle \sum }_{i=1}^{n_t}\frac{o_{i,t}^{ToA2}}{p_{i,t}}$$

(16)

2.4.3. Estimation under A3

Under Assumption 3 (A3), we know if any harvest has occurred within the past two years and which year the harvest occurred. We will assume that proportional growth occurred before and after ICE harvest observations and that we can estimate the proportion harvested to within ±5% (truncated at 0% and 100%). Let ${\pi }_{i,y}^h$ be the ICE-estimated proportion of volume harvested in year y on plot i. Then, the weight for the harvest component under A3 is:

$${\pi }_{i,y}^{ThA3}=\{\begin{array}{cc}{\pi }_y& \text{if no ICE-identified harvest in observation interval}\hfill \\ {\pi }_{i,y}^h& \text{if ICE-identified harvest in year }y\hfill \\ 0& \text{otherwise}\hfill \end{array}$$

(17)

The harvest allocation on plot i for year t under A3 is:

$$h_{i,t}^{ThA3}={\pi }_{i,y}^{ThA3}{h}_i$$

(18)

This results in the estimator for the harvest component for year t under A3:

$$C_t^{hA3}=\frac{1}{n_t}{\displaystyle \sum }_{i=1}^{n_t}\frac{h_{i,t}^{ThA3}}{p_{i,t}}$$

(19)

To determine the allocation for the other components of change, let ${\pi }_{i,F}^h$ be the first (in year F, y_F) ICE-estimated proportion of volume harvested on plot i and ${\pi }_{i,S}^h$ be the second (in year S, y_S) ICE-estimated harvest proportion on plot i. We set ${\pi }_{i,F}^h=0$ if there is not an ICE-identified harvest during the interval, and we set ${\pi }_{i,S}^h=0$ if there are not two ICE-identified harvests during the interval. The proportion remaining following the first harvest on plot i is ${\pi }_{i,F}^R=1-{\pi }_{i,F}^h$. Likewise, the proportion remaining following the second harvest on plot i is ${\pi }_{i,S}^R=1-{\pi }_{i,S}^h$. We calculate the annual allocation in three parts:

$$a_y=\{\begin{array}{cc}\left({\pi }_{i,F}^h{\pi }_y\right)/{\displaystyle \sum }_{j=1}^{y_F}{\pi }_j& \text{if one or two ICE-identified harvests and }y\le y_F\\ 0& \text{otherwise }\end{array}$$

(20)

$$b_y=\{\begin{array}{cc}\left({\pi }_{i,F}^R{\pi }_{i,S}^h{\pi }_y\right)/{\displaystyle \sum }_{j=1}^{y_S}{\pi }_j& \text{if two ICE-identified harvests and }y\le y_S\\ 0& \text{otherwise }\end{array}$$

(21)

and:

$$c_y=\left({\pi }_{i,F}^R{\pi }_{i,S}^R{\pi }_y\right)$$

(22)

We then sum the parts to obtain the annual proportion:

$${\pi }_{i,y}^{ToA3}=a_y+b_y+c_y$$

(23)

The allocation under A3 then becomes:

$$o_{i,t}^{ThA3}={\pi }_{i,y}^{ToA3}{o}_i$$

(24)

Finally, the estimator for each of the other components under A3 is:

$$C_t^{oA3}=\frac{1}{n_t}{\displaystyle \sum }_{i=1}^{n_t}\frac{o_{i,t}^{ToA3}}{p_{i,t}}$$

(25)

2.5. Compatibility and the Estimation of Initial Volume

In the previous subsections, we have described the estimators for the components of change. A complete estimation system also requires an estimator for initial volume. If the estimation system is also required to be compatible, as defined in Section 1, then some means of ensuring compatibility is required. We caution the reader that to ensure compatibility in an equation system is to enforce a constraint, and every constraint leads to a sub-optimization of one or more of the system’s estimators.

The plot sample provides single-panel direct annual estimates of standing volume, while we desire estimates of standing volume at the beginning of the year or prior to the growing season. We assume that each panel’s estimate of mean standing volume is a mid-year estimate, because sample plots are measured throughout the year. There are a number of approaches that could be taken, but the simplest approach would be to take the mean of successive panel means to obtain the initial annual estimates. For i = 1998 to 2006:

$${\widehat{Y}}_{i+1}^P=.5\left({\overline{v}}_i^P+{\overline{v}}_{i+1}^P\right)$$

(26)

where ${\overline{v}}_j^P$ is the (assumed mid-season) standing volume estimate for panel j. We denote this series of estimates as ${\widehat{Y}}^P$.

Alternatively, because ${\widehat{Y}}^P$ does not ensure compatibility of the estimation system, we can use recursive estimation to estimate initial annual volume. That is, we can start with an initial volume estimate for a particular year and successively apply the annual estimators of the components of change to obtain successive estimates of initial annual volume. In the simulation, we started with the initial volume estimate for the year 2003, (${\widehat{Y}}_{2003}^R={\widehat{Y}}_{2003}^P$), which is in the center of the estimation interval. To obtain the recursive series of estimates ( ${\widehat{Y}}^R$) for 1999 to 2007, we used the algorithm:

For i = 2003 to 2006:

$${\widehat{Y}}_{i+1}^R={\widehat{Y}}_i^R+{\widehat{L}}_i+{\widehat{E}}_i-{\widehat{M}}_i-{\widehat{H}}_i$$

(27)

For earlier years we use, for i = 2003 to 2000:

$${\widehat{Y}}_{i-1}^R={\widehat{Y}}_i^R-{\widehat{L}}_{i-1}-{\widehat{E}}_{i-1}+{\widehat{M}}_{i-1}+{\widehat{H}}_{i-1}$$

(28)

2.6. Estimation Systems

The estimators described in the preceding sections can be combined in various ways to define an estimation system. Furthermore, variance reduction can be achieved in the estimators by combining successive estimates with a moving window estimator. For any initial estimator for time t (${\epsilon }_t$), a moving window estimator of size s (s is an odd positive integer) is:

$${\epsilon }_i^s={\displaystyle \sum }_{i=t-\left(\left(s-1\right)/2\right)}^{t+\left(\left(s-1\right)/2\right)}{\epsilon }_i$$

(29)

The estimation systems that we tested in the simulations will be easiest to follow if we define an estimation system nomenclature. In the nomenclature, the character “E” is followed by a 4-character code (e.g. E1234). The 1st character is for harvest; the 2nd character is for the other growth components; the 3rd character is for initial volume; and the 4th character is for the moving window size:

• ECCP1 = C^M for all change components and ${\widehat{Y}}^P$ for initial annual volume,
• ECTR1 = C^M for harvest, C^T for the other change components, ${\widehat{Y}}^R$ for annual volume,
• E11R1 = C^hA1 for harvest, C^oA1 for the other change components, ${\widehat{Y}}^R$ for annual volume,
• E22R1 = C^hA2 for harvest, C^oA2 for the other change components, ${\widehat{Y}}^R$ for annual volume,
• E33R1 = C^hA3 for harvest, C^oA3 for the other change components, ${\widehat{Y}}^R$ for annual volume,
• ECCP3 = 3-year moving window on ECCP1,
• ECCP9 = 9-year moving window on ECCP1,
• ECTR9 = 9-year moving window on ECTR1,
• E11R9 = 9-year moving window on E11R1,
• E22R9 = 9-year moving window on E22R1,
• E33R3 = 3-year moving window on E33R1,
• E33R5 = 5-year moving window on E33R1 and,
• E33R9 = 9-year moving window on E33R1.

2.7. Estimator Evaluation

For each iterate, for each year, we calculated the empirical bias (EB) and the empirical mean squared error (MSE), over the 1000 iterations, between each estimator and the true population values under each of the four error structures.

That is:

$$E{B}_{PES}=\frac{1}{1000}{\displaystyle \sum }_{i=1}^{1000}\left({\widehat{x}}_{PE{S}_i}-X_P\right)$$

(30)

where ${\widehat{x}}_{PE{S}_i}$ is the sample estimate of any variable, X, in population P for estimator E, under error structure S for iterate i. Likewise:

$$MS{E}_{PES}=\frac{1}{1000}{\displaystyle \sum }_{i=1}^{1000}{\left({\widehat{x}}_{PE{S}_i}-X_P\right)}^2$$

(31)

Here, we use the empirical mean squared error as the overriding criterion for judging the effectiveness of the external information while being cognizant of any introduced bias. Estimator robustness was tested in a simulation by sampling the population under four different assumptions of sampling error structure for the plot data and three different assumptions of the accuracy of the remotely-sensed data. Each simulation consisted of 1000 iterations of 1000 plots each (without replacement) from the population.

3. Results

Figure 1 gives the empirical means of the components of change and initial annual volume plotted with the true population mean on the left-hand side and mean squared errors on the right-hand side of the initial annual estimators (ECTR1, ECCP1, E11R1, E22R1 and E33R1) for estimation years 1999 through 2007. The results are given subsequent to 1000 iterations of 1000 samples each under Sample Error Structure 1. The corresponding results for these components of change for Error Structures 2 through 4 are available from the first author upon request. They are not given here, in the interest of brevity, because they are very similar to the results for Error Structure 1.

The differences between each annual mean and the true population mean on the left-hand side of Figure 1 is the empirical bias (EB_PES) described above, for each estimator. These estimation years are shown because they are the years that are the center of at least 1 observation window (or panel), for a sample drawn under this design and sample error structure from a population spanning 1995 to 2012. From Figure 1, it is obvious that the ICE observations usually served to lower the bias for estimates of annual harvest. This usually, though not exclusively, also resulted in lowered empirical mean squared errors for the annual harvest estimates. Figure 1 also shows that the general trend model had a positive effect on mortality estimates, while using the trend model for entry and live growth showed mixed results.

This work was predicated on the assumption that the improvement in harvest estimates obtainable through the use of the ICE observations will also result in better estimates of annual volume. The bottom row of Figure 1 investigates that assumption. Of the single-panel estimators for initial (1 January) annual volume in the bottom row of Figure 1, ECCP1 can be seen to be the most variable.

Figure 1. The empirical means of the initial annual estimators (ECTR1, ECCP1, E11R1, E22R1 and E33R1) plotted with the true population mean for (a) live growth; (b) entry; (c) mortality; (d) harvest; and (e) initial volume; and the corresponding mean squared errors for (f) live growth; (g) entry; (h) mortality; (i) harvest and (j) initial volume, from the simulation of 1000 iterations of 1000 samples each under Sample Error Structure 1.

Figure 1. The empirical means of the initial annual estimators (ECTR1, ECCP1, E11R1, E22R1 and E33R1) plotted with the true population mean for (a) live growth; (b) entry; (c) mortality; (d) harvest; and (e) initial volume; and the corresponding mean squared errors for (f) live growth; (g) entry; (h) mortality; (i) harvest and (j) initial volume, from the simulation of 1000 iterations of 1000 samples each under Sample Error Structure 1.

Figure 2 gives the empirical means plotted with the true population mean on the left-hand side and mean squared errors on the right-hand side of the moving-window annual estimators (ECTR9, ECCP9, ECCP3, E11R9, E22R9, E33R3, E33R5 and E33R9) over 1000 iterations of 1000 samples, each under Sample Error Structure 1. As expected, the moving window estimators “smooth” the mean annual estimates and usually serve to lower the annual empirical mean squared errors for the components relative to the corresponding results in Figure 1. The bottom row of Figure 2 gives the results for the moving window annual estimators of annual volume. As in Figure 1, this row shows some advantage to combining panels. For instance, the MSEs are slightly less variable in the early years for ECCP3, relative to the results for ECCP1 in the top row of the figure. ECCP9 shows an even further reduction in MSE for most years, relative to ECCP3.

Figure 2. The empirical means of the moving window annual estimators (ECTR9, ECCP9, ECCP3, E11R9, E22R9, E33R3, E33R5 and E33R9) plotted with the true population mean for (a) live growth; (b) entry; (c) mortality; (d) harvest; and (e) initial volume; and the corresponding mean squared errors for (f) live growth; (g) entry; (h) mortality; (i) harvest and (j) initial volume, from the simulation of 1000 iterations of 1000 samples each under Sample Error Structure 1.

Figure 2. The empirical means of the moving window annual estimators (ECTR9, ECCP9, ECCP3, E11R9, E22R9, E33R3, E33R5 and E33R9) plotted with the true population mean for (a) live growth; (b) entry; (c) mortality; (d) harvest; and (e) initial volume; and the corresponding mean squared errors for (f) live growth; (g) entry; (h) mortality; (i) harvest and (j) initial volume, from the simulation of 1000 iterations of 1000 samples each under Sample Error Structure 1.

Although we are only showing the results for Error Structure 1, we note that the ranking of the mean for each annual estimator remained constant through the four error structures, although the position of each annual group of estimators did change slightly with error structure.

4. Discussion and Conclusions

Roesch [6] explored some special problems that arise in estimation of the components of change when the temporal scale of the population estimand of interest is finer than the scale of observation under both biased and unbiased sampling error structures. Here, we go further and attempt to ameliorate the effects of an unplanned increase in the temporal interval between observations. In the example simulations, the temporal scale of observation was increased from five to 10 years, while the temporal dimension of the population of interest and the estimands of interest remained fixed at one year.

As discussed in Roesch [6], by definition, the MW estimator is non-centralized and uses the information for all plots whose intervals span a particular year. This can allow the use of more data in the years at the extremes of a period of interest, but will lead to a reduced ability to detect trend changes. This effect is intensified for longer cycle lengths. That is, when the design observes five-year windows, an estimate of average annual change “smooths” the actual annual change, and the estimators, while drawing strength from overlapping panels, provide further smoothing. Both of these effects are substantially increased when the cycle length is changed to 10 years. The simulations showed the variance/bias trade-off encountered when the moving-window mean of ratios estimator was used in the extremity years of observation. Although the MWMOR estimator is sometimes biased in the presence of the trend in the extremity years, the empirical mean-squared error was often much lower when the extraneous information was used. These simulations have accentuated the difficulty that exists when attempting to compensate for an inadequately-informed sample. Large observation intervals frustrate the estimation of even the simplest of trends. We investigated the estimator performance of compatible annual estimators of the components of change during years spanned by at least half of the full set of plot observations, without attempting to address the additional problem arising from the known loss of adequate information in the extreme years of the monitoring effort, which resulted from the lengthened cycles. Never-the-less, we have shown that highly-informative auxiliary observations with shorter observation intervals (as was available under ICE Assumption 3) can contribute to significant improvement in estimation.

Auxiliary information, such as simulated under ICE, may arise from several remote sensing efforts. For example, the National Agriculture Imagery Program (NAIP) collects fine-scale imagery across the United States on a two- to three-year rotating schedule. The manual interpretation of these data for land use, land cover and change information would provide auxiliary observations similar to the ICE scenarios presented in these analyses. Data from other remote sensing efforts may also be appropriate. For example, Li et al. [19] used a time series of Landsat TM imagery to classify changes in forest vegetation and the year that change occurred. Model output arising from Li’s et al. [19] approach could also be informative in our application if modeling error in the classification is in line with the ICE scenarios.

Although our focus in this research was on estimating components of change, ICE-type information, such as described by Webb et al. [20], could aid in estimating several other parameters of interest. For example, forestland provides a suite of ecosystem services, and changes in the use of the forestland base influence the availability of those services (Coulston et al. [21]). Van Deusen and Roesch [4] and Coulston et al. [21] provide temporally-specific approaches to estimate changes in the use of forestland, based on remeasured plot data alone. Including auxiliary data on land use and/or land cover change may improve the precision of land change estimates. However, further development and testing of appropriate temporally-specific estimators is required to take full advantage of the auxiliary information. Temporally-specific estimators serve to increase the relevance and timeliness of NFI data. Using temporally-dense remote sensing data to assist in the estimation may provide an opportunity for national forest inventories whether or not fiscal climates are uncertain. However, the inclusion of auxiliary data that is more temporally dense than the plot data has the potential to provide improvement only if time is specifically accounted for in the estimation process, as is the case here.