Statistical properties of hybrid estimators proposed for GEDI—NASA's global ecosystem dynamics investigation

The GEDI instrument will produce ~25 m (diameter) footprint observations along tracks paralleling the ISS obit. The intersection of a track with a 1 km cell produces a collection of GEDI footprints, to be denoted here as a cluster. We propose to view each 1 km grid cell as tessellated into equal-area non-overlapping population elements, each representing a potential GEDI footprint. GEDI footprints are population elements sampled by the GEDI instrument. Figure 2 illustrates this conception, adapted to the spatial dimensions of lidar and field data available for this study. Footprint (20 m) and grid cell dimensions (approximately 800 m) used here (illustrated in figure 2) were considered sufficiently similar to GEDI to support applied testing of the estimators described in the following. For simplicity, we consider tracks as having either a 45° or 135° degree inclination for descending or ascending orbits, respectively, although ISS inclination patterns will vary across latitudes. In this idealization, we assert that there are two potential disjoint clusters for the intersection of a track with a grid cell—one starting on the border footprint, and one starting with the subsequent along-track footprint. The sample will be the set of clusters from the tracks that intersect the 1 km grid cell; due to variability of the ISS orbit this sample can be considered a simple random sample of clusters of GEDI footprints. In this exercise, we do not consider missing data issues that will arise because of clouds.

Zoom In Zoom Out Reset image size

The population parameter of interest is the average of the true AGBD over the N population elements (potential GEDI footprints), where N is the total number of population elements in a grid cell. For a 1 km grid cell, N is equal to 2500 20 × 20 m population elements, where the footprints in this realization would have a 20 m diameter. There is a parametric model, g(x, α), of AGBD in Mg ha⁻¹, where x are characteristics of the GEDI waveform for the GEDI footprint, and α is a vector of parameters. The true AGBD for a footprint is assumed to deviate from the g(x, α) by a random quantity, ε, with an expected sum of zero. Under the assumption that the model is an unbiased predictor of the true AGBD, and ignoring model misspecification error, the modeled value g(x, α) will be used in place of the true AGBD. The mean can be expressed in terms of GEDI ground clusters. There are M non-overlapping clusters distributed across all possible tracks. There are T_i footprints in the ith cluster. The reader may have noted that each population element is in one ascending cluster and one descending cluster. We can express the population attribute, the average expected predicted AGBD over the N population elements, in terms of the clusters; that is:

$$\begin{eqnarray}&&{\mu }_{Y}=\displaystyle \frac{\displaystyle {\sum }_{i=1}^{M}\displaystyle {\sum }_{t=1}^{{T}_{i}}g\left({{\boldsymbol{x}}}_{it},{\boldsymbol{\alpha }}\right)}{2N},\end{eqnarray} \tag{ 1 }$$

where the modeled AGBD for each population element occurs in the numerator twice.

Equation (1) will be expressed in an equivalent form, for which there is an estimator and a method of calculating an approximate variance. Let ${G}_{i}\,=\displaystyle {\sum }_{t=1}^{{T}_{i}}g\left({{\boldsymbol{x}}}_{it},{\boldsymbol{\alpha }}\right),$ equal the cluster total of the predicted AGBD per hectare for the footprints in the ith cluster. Equation (1) can be expressed as the ratio of the mean (or sum of) AGBD per cluster and the mean (or sum of) number of elements per cluster:

$$\begin{eqnarray}&&{\mu }_{Y}=\displaystyle \frac{\displaystyle {\sum }_{i}^{M}{{G}}_{i}}{\displaystyle {\sum }_{i}^{M}{T}_{i}}=\displaystyle \frac{{M}^{-1}\displaystyle {\sum }_{i}^{M}{{G}}_{i}}{{M}^{-1}\displaystyle {\sum }_{i=1}^{M}{T}_{i}}.\end{eqnarray} \tag{ 2 }$$

Parameter estimates, $\hat{{\boldsymbol{\alpha }}},$ for the model g(x, α) are based on a separate sample, S_m, that is independent of the 1 km grid cell. In equations (1) and (2) the parameters, α, are replaced by the parameter estimates, $\hat{{\boldsymbol{\alpha }}},$ and the population number of clusters, M, by the sample number of clusters, m, to get the estimator ${\hat{\mu }}_{Y},$ i.e. G_i with ${\hat{G}}_{i}=\displaystyle {\sum }_{t=1}^{{T}_{i}}g\left({{\boldsymbol{x}}}_{it},\hat{{\boldsymbol{\alpha }}}\right)$

$$\begin{eqnarray}&&{\hat{\mu }}_{Y}=\displaystyle \frac{\displaystyle {\sum }_{i}^{m}\displaystyle {\sum }_{t=1}^{{T}_{i}}g\left({{\boldsymbol{x}}}_{it},\hat{{\boldsymbol{\alpha }}}\right)}{\displaystyle {\sum }_{i=1}^{m}{T}_{i}}=\displaystyle \frac{\displaystyle {\sum }_{i=1}^{m}{\hat{G}}_{i}}{\displaystyle {\sum }_{i=1}^{m}{T}_{i}}.\end{eqnarray} \tag{ 3 }$$

The second component of the equation expresses the estimator ${\hat{\mu }}_{Y},$ the sample mean of all predicted AGBD (per sample element) which is an estimator of the population average expected AGBD (per element). The population true average AGBD differs randomly from the expected value ${\mu }_{Y}$ by the population average $\bar{\varepsilon }$ of the N deviations ε_it. For large N the contribution of $\bar{\varepsilon }$ to the total random error can be assumed negligible compared to the other two sources of uncertainty (the sampling error and the variability due to the model, see below), even if spatial autocorrelation is present. Hence, it is sufficient to assume the size of the grid cell (and thus N) is large enough to imply that the average $\bar{\varepsilon }$ will be close to zero. This assumption, along with the assumptions related to model fit will be covered in the discussion section.

As noted in the Introduction, the GEDI sample will be treated as a random sample of clusters of GEDI footprints. The estimator in equation (3) can be viewed as estimating the mean of cluster totals and mean number of footprints separately, and then combining as a ratio estimator:

$$\begin{eqnarray}&&{\hat{\mu }}_{Y}=\displaystyle \frac{\displaystyle {\sum }_{i}^{m}{\hat{G}}_{i}}{\displaystyle {\sum }_{i=1}^{m}{T}_{i}}=\displaystyle \frac{{m}^{-1}\displaystyle {\sum }_{i}^{m}{\hat{G}}_{i}}{{m}^{-1}\displaystyle {\sum }_{i=1}^{m}{T}_{i}}=\displaystyle \frac{\bar{\hat{G}}}{\bar{T}}.\end{eqnarray} \tag{ 4 }$$

This estimation strategy combines a probability sample of auxiliary information (as opposed to wall-to-wall auxiliary information) with a prediction for the population element value, AGBD, instead of directly observing AGBD. Ståhl et al (2011) proposed an estimator of the form of equation (4) and derived an estimator of the approximate variance. Due to the small sample size used in that study in relation to the population size, Ståhl et al (2011) did not include a finite population correction. A finite population correction may need to be taken into account in GEDI's case. Expression (5) below is the variance estimator proposed in Ståhl et al (2011), with the addition of the finite population correction:

$$\begin{eqnarray}\begin{array}{lll}\hat{V}\left({\hat{\mu }}_{Y}\right) & = & \displaystyle \frac{1}{{\bar{T}}^{2}}\left(1-\displaystyle \frac{m}{M}\right)\displaystyle \frac{\displaystyle {\sum }_{i=1}^{m}{\left({\hat{G}}_{i}-{\hat{\mu }}_{Y}{T}_{i}\right)}^{2}}{m\left(m-1\right)}\\ & & +\,\displaystyle \frac{1}{{\bar{T}}^{2}}\displaystyle \sum _{j=1}^{p}\displaystyle \sum _{k=1}^{p}{\widehat{{\rm{C}}{\rm{o}}{\rm{v}}}}_{{S}_{m}}\left({\hat{\alpha }}_{j},{\hat{\alpha }}_{k}\right){\widehat{\bar{G}}}_{j}^{{\prime} }{\widehat{\bar{G}}}_{k}^{{\prime} }\end{array}\end{eqnarray} \tag{ 5 }$$

The first term is due to the sampling error and the second term is due to effects of uncertainty of the $\hat{{\boldsymbol{\alpha }}}$ estimates (Ståhl et al 2011). The $\widehat{{\rm{C}}{\rm{o}}{\rm{v}}}\left(\hat{{\boldsymbol{\alpha }}}\right)$ is the estimated covariance matrix of the p-parameter estimates, where the estimate is based on the separate set of data, S_m, that is independent of the 1 km grid cell. For simplicity the first term is referred as the sample component and second term the model component. If we assume a linear model, that is $g\left({{\boldsymbol{x}}}_{it},\hat{{\boldsymbol{\alpha }}}\right)\,={\hat{{\boldsymbol{\alpha }}}}_{1}+\displaystyle {\sum }_{j=2}^{p}{x}_{itj}{\hat{{\boldsymbol{\alpha }}}}_{j},$ where x_itj is the jth component of x_it, then, ${\widehat{\bar{G}}}_{j}^{{\prime} }=\displaystyle \frac{1}{m}\displaystyle {\sum }_{i=1}^{m}\displaystyle {\sum }_{t=1}^{{T}_{i}}{x}_{itj}.$ Let $\widehat{\bar{{\boldsymbol{x}}}}$ be the vector of the means of the cluster totals of the predictor variables, i.e. ${\widehat{\bar{{\boldsymbol{x}}}}}_{j}=\displaystyle \frac{1}{m}\displaystyle {\sum }_{i=1}^{m}\displaystyle {\sum }_{t=1}^{{T}_{i}}{x}_{itj},$ for j = 1, ..., p. Then the second component on the right side of equation (5) can be expressed in matrix notation:

$$\begin{eqnarray}&&\displaystyle \sum _{j=1}^{p}\displaystyle \sum _{k=1}^{p}{\widehat{{\rm{C}}{\rm{o}}{\rm{v}}}}_{{s}_{m}}\left({\hat{\alpha }}_{j},{\hat{\alpha }}_{k}\right){\widehat{\bar{G}}}_{j}^{{\prime} }{\widehat{\bar{G}}}_{k}^{{\prime} }={\hat{\bar{{\boldsymbol{x}}}}}^{T}\widehat{{\rm{C}}{\rm{o}}{\rm{v}}}\left(\hat{{\boldsymbol{\alpha }}}\right)\widehat{\bar{{\boldsymbol{x}}}}.\end{eqnarray} \tag{ 6 }$$

Note that ${\widehat{\bar{{\boldsymbol{x}}}}}_{1}$ is equal to $\bar{T}.$ In the construction of equation (5), a Taylor linearization process is used twice; first with regards to the model parameters of $g\left({\boldsymbol{x}},\hat{{\boldsymbol{\alpha }}}\right)$ (Ståhl et al 2011) and second to estimate variance of a ratio estimator (Särndal et al 1992).

The empirical data described in section 2.2 were used to calibrate the simulation study described in section 2.3 to evaluate the statistical properties of these estimators in the GEDI setting.

Lidar and field data collected over six diverse areas in the United States supported this study. These areas were defined by the non-overlapping window (Thiessen Scene Area) of a local Landsat scene. The six scenes (given by numeric WRS-2 Path/Row) included one each in: northern Maine ('ME': 12/28, excluding the Canadian portion); eastern Pennsylvania and central New Jersey ('PA/NJ': 14/32); coastal South Carolina ('SC': 16/37); northern Minnesota ('MN': 27/27); northwestern Colorado ('CO': 35/32); and western Oregon ('OR': 45/30) (figure 3). The selected areas represented a wide range of forest ecosystems and disturbance processes, as described by Cohen et al (2017) and Healey et al (2018). Small-footprint lidar was collected at each site in the pattern displayed in figure 3, and a waveform simulator developed by the GEDI Science Definition Team was used to simulate GEDI waveforms from the airborne data. Data collection and waveform simulation is described in appendix S1.

Zoom In Zoom Out Reset image size

Simulations were conducted to examine the properties of the above estimators of mean AGBD and variance of the estimated mean. Modeled AGBD values treated as truth for each population element were developed for sixty square grid cells (ten randomly chosen grid cells at each of the six test areas in figure 3) of dimensions approximately similar to GEDI's product specifications. These populations were created by applying a linear model to the GEDI metrics simulated from actual small-footprint lidar acquisitions (figure 2). Appendix S2 details the process of both establishing the true population and simulating different potential cluster patterns and model parameter combinations to test the proposed estimators under a range of conditions. The advantage of using real data to develop the 'true' population, as opposed to simulating arbitrary AGBD surfaces, was that it ensured realistic covariance among lidar metrics. To the degree that lidar metrics are correlated with AGBD (see appendix S3), this approach also created a realistic spatial representation of biomass and of residual prediction error when different realizations of model parameters were applied to the lidar metrics. The swath width of the available lidar data limited the size of these test cells to between 680 × 680 m and 800 × 800 m (instead of 1 km); as mentioned earlier, this was considered adequate for the tests described below.

Earlier work with NASA's Geoscience Laser Altimeter System instrument pioneered the use of spaceborne lidar to make forest AGBD maps, often employing models that related field biomass to lidar and then, in a separate stage, related lidar-modeled biomass to synoptic data from sources like the MODIS satellites. While some work has attempted to impose a statistical framework on the use of such maps (e.g. Nelson et al 2017), most have employed ad hoc approaches to uncertainty. Mitchard et al (2014) pointed out two such ad hoc error budgeting approaches that yielded estimates across the Amazon that did not overlap with each other's relatively narrow confidence intervals, suggesting unacknowledged uncertainty in one or both estimates.

Saarela et al (2016) concluded that ignoring error occurring in the field-to-lidar step, as is often done, can lead to underestimation of variance by a factor of three. The spatial mismatch inherent in modeling AGBD in large pixels from constituent smaller-footprint measurements can also introduce uncertainties not clearly addressed in variance estimates (Réjou-Méchain et al 2014). The proposed hybrid approach does explicitly account for field-to-lidar model error, and it also accounts for sample error as footprints are combined to infer mean 1 km AGBD.

However, the hybrid approach summarized in figure 1 is subject to at least three important assumptions. The first is that the GEDI waveform simulator (Hancock et al 2019; described in appendix S2) correctly translates pulse counts to waveform data while representing the important sources of error. The benefit of using simulated GEDI data at precisely known locations to build footprint-level models is that it obviates potential error related to the positional uncertainty of GEDI footprints (currently estimated to be on the order of 10 m, post-processed). The risk is that any systematic simulator error may be propagated through footprint-level models to bias the 1 km mean AGBD estimates. However, the simulator is based on well-validated work by Hofton and Blair and Hofton (1999) and it is unlikely that violations of this assumption are substantial.

A second assumption is that one square kilometer is a large enough domain, given spatial autocorrelation of the population, that residual model error may be considered negligible. Mean residual error in the footprint AGBD model is presumed to tend toward zero over a large number of predictions. In an area that is small, particularly if spatial autocorrelation of residual error is high, residual errors may not sum to zero, introducing additional error. While GEDI's 1 km spatial domain is smaller than other applications of hybrid inference (e.g. Ståhl et al 2011), initial exploration using the data from this study suggests that spatial autocorrelation at GEDI's 25 m grain size is relatively low, as is the effect of omitting residual error (also called individual error in appendix 3) from the proposed hybrid variance estimator (appendix 3). In some cases, however, the impact of residual error may reach 15% or more of the variance, so further work is needed to develop ways to include it in the hybrid approach.

Finally, it is assumed that the footprint-level AGBD model is correctly specified and 'applies' to the population in the cell. This notion implies that the number and order of terms are correctly specified, and that the parameter covariance matrix used in the variance calculation (equations (5) and (6)) adequately reflects uncertainty in the relationship between lidar and AGBD within the cell. Since even small model difference can generate strongly divergent population estimates when multiplied over large areas, this is a critical assumption. Although the GEDI team is collating a comprehensive global database of training data, the reality is that ground data are sparse in some areas. In such areas, model misspecification is an important risk. In this respect, hybrid inference is no different from any other remote sensing approach that relies upon field data to calibrate remotely sensed observations.

Properties of the hybrid estimators proposed for the GEDI mission were evaluated here using simulations in which thousands of potential GEDI cluster patterns were tested in the context of model covariance across forests in 60 diverse grid cells. Bearing in mind the above assumptions associated with the approach outlined in figures 1, 4 and 5 suggest that the proposed estimator for mean AGBD exhibits negligible bias across cells. Our tests and the results in figures 4 and 5 presumed proper model specification, an assumption of model-based inference, although residual plots (figure S4(1)) suggested that trends in model residual error may exist in at least one study region. In such cases, corrective measures such as variable transformations may be required.

The variance estimator, on the other hand, appeared to be only asymptotically unbiased (figure 6); that is, bias approached zero as larger numbers of cluster samples were acquired. The Taylor linearization used to compute the sample component of the variance (as opposed to the model component) likely contributed to this underestimation. While Taylor linearization simplifies calculation of the variance, it is known to lead to variance underestimation for small sample sizes (Särndal et al 1992, p 176). One area of future work may center around methods proposed to adjust for under-estimation (Cochran 1977, p 156). Linearization was also used in deriving the model component of the variance estimator (equation (5)), but the linear model used in this study for footprint-level prediction likely eliminated the impact of linearization in the model component.

Negative bias in the variance estimates at low numbers of clusters should be considered when hybrid methods are used with GEDI data, and information about the number of clusters will accompany GEDI hybrid estimates for each cell. Low sample numbers will occur both at the beginning of the mission and in equatorial regions where the ISS orbit provides fewer 'looks' and where persistent cloud cover may frequently obscure the surface. An advantage of hybrid estimation is that the same methodology applied at the 1 km scale may also be applied to much larger, irregularly shaped areas such as entire countries; instead of having to combine grid cell estimates, all clusters intersecting a given country may be used in the context of a single hybrid estimate. GEDI country-scale sample size will be very large, leading to greatly diminished design-based contribution to variance estimates. The apparent dominance of the sample component in the variance estimator (figure 8) may be relevant to GEDI operational planning. There are signal-to-noise parameters that may allow GEDI to consider a larger sample of noisier footprints (and models) or a smaller sample of footprints collected under ideal atmospheric conditions. It should be recognized that the relative contribution of the sampling and model components of the variance estimator may differ in different settings as properties of the model change due, for example, to the size of S_m or complexity of the underlying lidar/forest structure relationship. Across the diverse ecosystems we studied, though, results suggest that including a larger sample will have a bigger impact on the magnitude and unbiasedness of estimated variance than a marginal improvement in model fit.