**Remote Sensing Support for the Gain-Loss Approach for Greenhouse Gas Inventories**

## **Ronald E. McRoberts 1,2,\*, Erik Næsset 3, Christophe Sannier 4, Stephen V. Stehman <sup>5</sup> and Erkki O. Tomppo 6,7**


Received: 21 April 2020; Accepted: 8 June 2020; Published: 11 June 2020

**Abstract:** For tropical countries that do not have extensive ground sampling programs such as national forest inventories, the gain-loss approach for greenhouse gas inventories is often used. With the gain-loss approach, emissions and removals are estimated as the product of activity data defined as the areas of human-caused emissions and removals and emissions factors defined as the per unit area responses of carbon stocks for those activities. Remotely sensed imagery and remote sensing-based land use and land use change maps have emerged as crucial information sources for facilitating the statistically rigorous estimation of activity data. Similarly, remote sensing-based biomass maps have been used as sources of auxiliary data for enhancing estimates of emissions and removals factors and as sources of biomass data for remote and inaccessible regions. The current status of statistically rigorous methods for combining ground and remotely sensed data that comply with the good practice guidelines for greenhouse gas inventories of the Intergovernmental Panel on Climate Change is reviewed.

**Keywords:** statistical estimator; IPCC good practice guidelines; activity data; emissions factor; removals factor

#### **1. Introduction**

The Intergovernmental Panel on Climate Change (IPCC) [1] (p.17) identifies five carbon pools to be monitored for forest-related carbon emissions and removals: aboveground biomass, below-ground biomass, litter, dead wood and soil organic carbon. Within the agriculture, forestry and other land use sector, greenhouse gas (GHG) inventories for all these pools are typically conducted using either the stock difference approach or the gain-loss approach [1] (p. 22), [2] (Vol. 4, Chap 2, p. 2.10). With the stock-difference approach, mean annual carbon emissions or removals are estimated as the ratio of the difference in carbon stock estimates at two points in time and the number of intervening years [1] (Chap. 3), [2] (Vol. 4, Chap. 3). For countries with established and comprehensive forest sampling programs, such as national forest inventories (NFI), the stock difference approach is fairly easy to implement. However, for tropical countries that do not have NFIs or that only have a single set of NFI measurements, the gain-loss approach is used more often. With the gain-loss approach, emissions are estimated as differences between additions to and removals from carbon pools. Specifically, emissions are estimated as the product of activity data defined as the areas of "human activity causing emissions and removals" and emissions factors defined as the per unit area responses of carbon stocks for those activities [1] (pp. xvii, 22], [2] (Vol. 1, Chap. 1, Sect. 1.2).

The IPCC characterizes GHG methods for estimating emissions and removals with respect to three tiers or levels of data detail and analytical complexity [1] (p. 19), [3]. Tier 1 is the default method and permits use of default, national, ecological zone, or global emissions factors based on canopy cover reductions, thus making it applicable for any country. However, emissions factors from these sources are subject to considerable uncertainty. Tier 2 is based on the same conceptual structure as Tier 1, but with the expectation that activity data aggregate land use changes between categories and emissions factors are based on national level data. Tier 3 uses models, data from repeated inventories, and fine resolution land use and land use change activity data to produce spatially continuous, sub-national estimates. Gain-loss methods can be implemented at all three tier levels.

Within tiers, the IPCC further documents three approaches for estimating activity data, all involving land use class areas [1] (p. 19). Approach 1 uses aggregations of land use class areas such as would be reported by an NFI, but without regard to the specific geographic locations of those areas of those land use classes; Approach 2 differentiates among land use change classes, but like Approach 1, does not identify the specific locations of individual segments of the classes; and Approach 3 tracks individual land parcels over time.

Regardless of the approach and method, the IPCC specifies two good practice guidelines for GHG inventories: (i) "neither over- nor underestimates so far as can be judged," and (ii) "uncertainties are reduced as far as practicable" [1] (p. 15), [2] (Vol 1. Chap. 1, Sect. 1.2). From a statistical perspective, the satisfaction of these criteria requires use of unbiased estimators and, at minimum, rigorous estimation of uncertainty. In particular, there can be little assurance that uncertainties are reduced until they are first rigorously estimated [4,5].

Considerable recent attention has been devoted to developing and illustrating methods for implementing the gain-loss approach, particularly for estimating activity data, but less so for estimating emissions and removals factors for the aboveground biomass pool [1] (p. 17). The objective of the review is to document and summarize the current role of remote sensing for gain-loss methods for Tier 1 and Tier 2 GHG inventories, with particular emphasis on the two IPCC good practice guidelines. A subordinate objective is to submit statistical guidance developed for the Methods and Guidance Document (MGD) of the Global Forest Observations Initiative to journal-level peer review [1]. As per the MGD, the review focuses on forest land use, particularly forest remaining forest and conversions from and to forest land use. The structure of the review closely follows the structure of the MGD, with an initial focus on estimation of activity data followed by estimation of emissions and removals factors. Within these two categories, the focus is on the statistical estimators, as is the MGD, with the discussion logically proceeding from less complex to more complex.

#### **2. Estimating Total Emissions**

The ultimate objective is an inference in the form of a confidence interval,

$$\text{H}^{\text{tot}} \neq \text{t} \cdot \sqrt{\text{Var}\left(\text{E}^{\text{tot}}\right)} \tag{1}$$

where Eˆtot is the estimate of total emissions and removals, Varˆ Eˆtot is an estimate of the corresponding variance, and t is a percentile from Student's t-distribution corresponding to the desired confidence and degrees of freedom determined by the sample size and sampling design. With the gain-loss approach, total emissions and removals are estimated as the sum over activities of the products of estimates of activity class areas and corresponding estimates of activity class emissions and removals factors,

$$\mathbf{\hat{E}}^{\text{tot}} = \sum\_{\mathbf{c}=1}^{\mathbf{C}} \mathbf{\hat{A}}^{\mathbf{c}} \cdot \mathbf{\hat{E}} \mathbf{\hat{F}}^{\mathbf{c}} \tag{2}$$

where c=1, ... , C indexes activity classes, Aˆ <sup>c</sup> is the estimate of the area for activity class c, and EFˆ <sup>c</sup> is the estimate of the emission factor for activity class c [1] (pp. xvii, 22), [2] (Vol.1, Chap.1, Sect. 1.2). An estimator of the variance of Eˆtot can be formulated using a first-order Taylor series as,

$$\mathrm{Var}\left(\hat{\mathbf{E}}^{\mathrm{tot}}\right) = \sum\_{\mathbf{c}=1}^{\mathsf{C}} \left[ \left(\hat{\mathbf{A}}^{\mathrm{c}}\right)^{2} \cdot \mathrm{Var}\left(\hat{\mathbf{E}}\mathbf{\hat{F}}^{\mathrm{c}}\right) + 2 \cdot \hat{\mathbf{A}}^{\mathrm{c}} \cdot \hat{\mathbf{E}}\mathbf{\hat{F}}^{\mathrm{c}} \cdot \mathrm{Cov}\left(\hat{\mathbf{A}}^{\mathrm{c}}, \hat{\mathbf{E}}\mathbf{\hat{F}}^{\mathrm{c}}\right) + \left(\hat{\mathbf{E}}\mathbf{\hat{F}}^{\mathrm{c}}\right)^{2} \cdot \mathrm{Var}\left(\hat{\mathbf{A}}^{\mathrm{c}}\right) \right] \tag{3}$$

Reference [6] (Appendix). If activity areas and emission factors are estimated independently, then Covˆ Aˆ c , EFˆ <sup>c</sup> = 0 may be assumed. The technical challenge with the gain-loss method is to formulate the activity data estimators, Aˆ <sup>c</sup> and Varˆ Aˆ <sup>c</sup> , and the estimators for the emissions and removals factors, EFˆ <sup>c</sup> and Varˆ EFˆ <sup>c</sup> . Multiple country-level examples of the use of the gain-loss method for estimating total emissions and removals can be found in the country Fact Sheets at the REDD+ Web Platform [7], although rigorous uncertainty estimation is not always documented. Illustrations that include uncertainty assessments and that document estimation details can be found in [1] (p. 158) and [6].

#### **3. Activity Data**

The estimation of activity data typically entails three components: reference data, auxiliary data, and statistical estimators. When estimating activity data, auxiliary data are typically acquired from external sources, so that the primary purpose of reference data is to serve as the basis for estimates. However, when estimating emissions and removals factors (Section 4), auxiliary data are often developed by the user, so that reference data may have multiple purposes, including as training data for constructing remote sensing-based maps, as map validation data, and as a source of data for calculating estimates. For estimating activity data, reference data can be obtained from both ground and remote sensing sources. The primary purpose of auxiliary data is to enhance or improve estimates based on reference data, most often to increase precision. Auxiliary data, like reference data, can be obtained from both ground and remote sensing sources.

#### *3.1. Data*

#### 3.1.1. Reference Data

For the estimation of activity data, the reference data are the most accurate information on land use change classes available and are typically acquired as measurements of ground plots, or more commonly as image interpretations. To facilitate design-based inference, reference data are most often acquired using a probability sampling design [8,9] (Section 3.2).

Regardless of the estimator used with reference data (Section 3.2), observation errors and other uncertainties in reference data tend to induce bias into estimators [10–13]. Sun et al. [14] reported that manual interpretations of Google Earth and other fine resolution imagery by three trained interpreters were not as reliable as ground observations for seven land cover classes in Central Asia. Thus, for most applications, ground reference data are considered to be preferable, although Foody [10,11] argues that even ground reference data may be subject to error. However, the acquisition of probability samples of ground reference data for remote and inaccessible regions may be prohibitively expensive, if not logistically infeasible. The common alternative is to acquire reference data in the form of visual image interpretations, albeit with the caveat that such reference data are of greater quality than

any auxiliary map data (Section 3.1.2), with respect to factors such as resolution and accuracy [1] (pp. 125, 139) [8,15,16].

Mowrer and Congalton [17] characterize reference data with non-negligible uncertainty as imperfect reference data. Næsset [18,19] reported that interpretations of crown coverage for structurally homogenous Norwegian boreal forests differed statistically significantly among interpreters and among different times of year for the same interpreter. Further, interpretations of broad tree species groups by 12 professional interpreters using stereo aerial photography produced only 31–79% agreement with field reference data. These results were consistent with the results of a review of more than 10 studies from the Nordic countries from the 1970s and 1980s. An additional finding was, as might be expected, that interpretations were less accurate for more complex forests. For five trained interpreters of videography, Powell et al. [20] reported interpreter disagreement of almost 30% for five land cover classes in the Brazilian Amazon. Pengra et al. [21] examined duplicate interpretations of land-cover reference class labels obtained by well-trained interpreters from a probability subsample of 2900 pixels (30 m × 30 m) of the United States of America. For the 8-class land cover legend, the overall agreement between interpreters was 88%, but agreement by class ranged from 46% for the disturbed class to 94% for water, with more prevalent classes having greater agreement than rare classes. Thus, reference data in the form of visual interpretations of remotely sensed data, even by well-trained professional interpreters, could be subject to substantial interpreter disagreement and error.

McRoberts et al. [13] used a combination of photo interpretations by multiple professional interpreters and simulation data to assess the effects of imperfect reference data on the bias and precision of estimators of land cover class proportions, to characterize conditions that affect the magnitudes of bias and precision, and to develop a variance estimator that incorporates the effects of interpreter error. Several relevant conclusions were drawn from the study. First, estimator bias resulting from interpreter error is greater for greater inequality in areas of land cover classes, greater for smaller map and interpreter accuracies, greater for fewer interpreters, and greater for greater correlations among interpretations for different interpreters. For some scenarios, seven or more interpreters were necessary to reduce biases to acceptable levels. Second, failure to incorporate the effects of interpreter error into variance estimators led to underestimates of standard errors by factors as great as 2.0.

The important lesson from these studies is that imperfect reference data induce bias into statistical estimators, sometimes substantial bias despite only small errors [10–13], and lead to non-compliance with the first IPCC good practice guideline. In addition, failure to incorporate the uncertainty associated with imperfect reference data into statistical variance estimators leads to under-estimation of variances and non-compliance with the second IPCC good practice guideline.

Multiple strategies for dealing with imperfect reference data may be considered. First, greater numbers of interpreters reduce the effects of imperfect reference data [13]. Second, prior to operational interpretation, interpreters can calibrate their interpretations, with respect to known field conditions and/or to each other [22] (p. 82). Third, in the absence of unanimous interpretations, interpreters may discuss the specific sample units and agree on a consensus interpretation. Finally, instead of using majority interpretations leading to categorical reference observations (e.g., 0 for non-forest, 1 for forest), continuous reference observations in the form of the proportions of forest interpretations among interpreters for the same sample unit are possible.

#### 3.1.2. Auxiliary Data

For purposes of estimating activity data, auxiliary data are secondary, i.e., the reference data are the primary source of information on which estimates are based, whereas the auxiliary data are used only to improve the estimation process by increasing the precision of estimates. Because activity data are estimates of change in areas of land use classes, both the reference and auxiliary data are closely related to change. In particular, auxiliary data are often in the form of a map or spatial product whose classes reflect or can be aggregated to reflect land use change classes related to human activities, such as deforestation, reforestation, afforestation, and forest remaining forest, including the special

case of forest degradation, which entails the substantial loss of biomass, but no conversion from forest to other land use.

Although the scope of this review does not include techniques and methods for constructing remote sensing-based maps, regional and global land use and land use change maps can be used as auxiliary data when estimating activity data. The most widely known example is the Global Forest Change dataset [23], a percentage tree cover map for the year 2000, in which trees are defined as vegetation taller than 5m in height at a 30-m pixel size. Annual forest gain and loss data from 2001 to 2019 are also included. Forest loss is a binary layer (1: loss, 0: no loss), and is understood as complete or the comprehensive removal of forest cover and "defined as a stand-replacement disturbance, or a change from a forest to non-forest state." Forest gain is the complete or comprehensive recovery of forest cover and "defined as the inverse of loss, or a non-forest to forest change entirely within the study period". Examples of uses of the Global Forest Change dataset for estimating land cover class areas include Sannier et al. [24], Næsset et al. [25], and McRoberts et al. [26]. In addition, the Global Forest Change dataset can be adapted to match national forest definitions [24,25,27]. Although few examples appear in the literature, maps for this purpose can be produced locally, and are typically more accurate than global maps [22,28]. McRoberts et al. [26] demonstrate how data from multiple maps can be combined to produce a new, more accurate map.

A potentially important issue is that while activity data are defined in terms of land use change, auxiliary data in the form of maps based on satellite spectral data inevitably depict land cover change. Further, land cover change does not necessarily correspond to land use change. For example, forest land that has been completely harvested loses its forest cover, but retains its forest land use. For use with the stratified (Section 3.2.3) and model-assisted (Section 3.2.4) estimators, auxiliary data in the form of land cover change maps in lieu of land use change maps do not induce bias into the estimators, but rather just reduce precision.

#### *3.2. Statistical Estimators*

#### 3.2.1. Design-Based Inference

Design-based inference, also characterized as probability-based inference, is based on three assumptions. First, a probability sample incorporating some form of randomization is used and constitutes the basis for validity. Second, apart from negligible observation or measurement error, each population unit is assumed to have one, and only one, possible value. Third, the selection of population units into the sample is based on positive and known probabilities of selection. Familiar sampling designs include simple random (SRS), systematic (SYS), stratified (STR), multi-phase and multi-stage sampling designs. The design-based simple expansion (Section 3.2.2), stratified (Section 3.2.3), post-stratified (Section 3.2.4) and model-assisted (Section 3.2.5) estimators herein considered for the estimation of activity data are either unbiased or asymptotically unbiased. The uncertainty estimation for these estimators entails comparing observations to their corresponding means or model predictions. Of importance, for statistical estimators to be unbiased, they must be consistent with the probability sampling design used to collect the reference data, i.e., the known probabilities of selection must be incorporated into the estimators. For example, the simple expansion estimator, Equation (6), is generally biased if used with reference data collected with a stratified sampling design.

If the reference data are categorical, then for land use change class, c, a new variable is defined, such that for each reference sample unit, y<sup>c</sup> <sup>i</sup> <sup>=</sup> 1 if class c is observed, or yc <sup>i</sup> = 0 if any class other than c is observed. If the reference data are continuous, such as differences in proportion forest cover at two times, then for each land use change class, c, the new variable takes on a value in the interval [0,1] corresponding to the difference in proportions of the reference sample unit in class c. Letting μ<sup>c</sup> be the proportion of the area of interest in class c, the estimator of the area of class c is,

$$\mathbf{A}^{\mathbf{c}} = \mathbf{A}^{\text{tot}} \cdot \mathbf{\hat{\mu}}^{\mathbf{c}} \tag{4}$$

with

$$\text{SE}(\hat{\mathbf{A}}^{\text{c}}) = \mathbf{A}^{\text{tot}} \cdot \text{SE}(\hat{\mathbf{\mu}}^{\text{c}}) \tag{5}$$

where Atot is the total area which is assumed to be without error. The activity data challenge, then, is to formulate the estimators, ˆμ<sup>c</sup> and SE(μˆ <sup>c</sup>).

#### 3.2.2. Simple Expansion Estimator

The simple expansion (EXP) estimator of μ<sup>c</sup> is,

$$\hat{\mathbf{n}}\_{\text{EXP}}^{c} = \frac{1}{\mathbf{n}} \sum\_{i=1}^{n} \mathbf{y}\_{i}^{c} \tag{6}$$

with

$$\text{SE}\{\hat{\mu}\_{\text{EXP}}^{c}\} = \sqrt{\frac{1}{\text{n} \cdot (\text{n} - 1)} \sum\_{i=1}^{n} \left(\mathbf{y}\_{i}^{c} - \hat{\mu}\_{\text{EXP}}^{c}\right)^{2}} \tag{7}$$

where n is the sample size. Of importance, the EXP estimators use no auxiliary data, which further establishes that the reference data are the primary source of information on which estimates are based.

#### 3.2.3. Stratified Estimators

Stratified (STR) estimation uses a map or similar spatial product to increase the precision of estimates. If the map classes are the same as the land use change activity classes, an intuitive estimator of μ<sup>c</sup> is simply the sum of the areas of all map units classified as activity class c. However, this estimator, characterized as pixel-counting, is biased, because it does not account for classification errors and, therefore, does not comply with the first IPCC good practice guideline pertaining to avoiding under- and over-estimation; this approach should be avoided. The reason for the bias is that due to classification error, the map class of interest may include some units that are not actually of the activity class of interest and, similarly, other map classes may include some map units that are of the activity class of interest.

For the purposes of estimating activity data, STR estimators assume that reference data are acquired using a STR sampling design, guided by an activity class map, for which the classes serve as strata. An advantage of the combination of stratified sampling and estimation is that within-stratum sample sizes can be controlled and selected for the primary purpose of more precisely estimating the areas of key land use change classes such as deforestation. In particular, when the class of interest is small, SRS and SYS sampling designs generally do not produce sample sizes sufficiently large enough to satisfy precision requirements. STR designs facilitate the allocation of greater sample sizes for smaller strata corresponding to land use change classes, and smaller sample sizes for larger strata corresponding to stable land use classes.

Cochran [29] (p. 134) recommends no more than 6–8 strata. Särndal et al. [30] (pp. 267, 407) recommend minimum within-stratum sample sizes of 10–20; Cochran [29] (p. 134) recommends minimum within-stratum sample sizes of 20; and for temperate forest inventories, Westfall et al. [31] recommend within-stratum sample sizes of at least 20.

The STR estimator of the proportion of the total area in class c is,

$$\mathfrak{\mu}\_{\rm STR}^{\rm c} = \sum\_{\mathbf{h}=1}^{\rm H} \mathbf{w}\_{\mathbf{h}} \cdot \mathfrak{\mu}\_{\mathbf{h}}^{\rm c} \tag{8}$$

with

$$\text{SE}\{\hat{\mu}\_{\text{STR}}^{\text{c}}\} = \sqrt{\sum\_{\mathbf{h}=1}^{\text{H}} \mathbf{w}\_{\text{h}}^{2} \cdot \frac{\hat{\sigma}\_{\text{h}}^{2}}{\mathbf{n}\_{\text{h}}}} \tag{9}$$

where h = 1, ... , H indexes the strata, wh is the stratum weight calculated as the proportion of the population in the *h*th stratum, and nh is the sample size for the *h*th stratum,

$$\mathfrak{h}\_{\mathrm{h}}^{\mathrm{c}} = \frac{1}{\mathrm{n}\_{\mathrm{h}}} \sum\_{i=1}^{n\_{\mathrm{h}}} \mathbf{y}\_{\mathrm{hi}'}^{\mathrm{c}} \tag{10}$$

$$\boldsymbol{\sigma}\_{\rm h}^{2} = \frac{1}{\mathbf{n}\_{\rm h} - 1} \sum\_{i=1}^{n\_{\rm h}} \left( \mathbf{y}\_{\rm hi}^{\rm c} - \boldsymbol{\mu}\_{\rm h}^{\rm c} \right)^{2},\tag{11}$$

and yc hi is the observation for the *i*th sample reference unit in the *h*th stratum. An additional advantage of the combination of stratified sampling and estimation is that the SEs are typically smaller than for the simple expansion estimators.

For many forests, both natural change and human-induced disturbance occur most frequently in the vicinity of the forest/non-forest boundary. In addition, classification uncertainty for most maps is greatest at boundaries between map classes. Finally, because of both classification error and geo-reference errors, the erroneous assignment of ground plots to strata inconsistent with their observations is often most severe near map class boundaries. Some of the adverse effects of these phenomena on precision can be alleviated by constructing an additional small stratum with corresponding small stratum weight along class boundaries. Thus, heterogeneous ground plots found near the boundaries can be confined to small strata with small stratum weights, thereby having less detrimental effects on precision. For a stratification based on a forest/non-forest map, McRoberts et al. [32] constructed a forest edge stratum consisting of a multi-pixel buffer on the forest side of the forest/non-forest map boundary and a non-forest edge stratum consisting of a multi-pixel buffer on the non-forest side of the boundary. The effect of the two additional strata was to decrease the SEs of the estimated forest area by as much as 12%.

These buffer strata have an additional beneficial effect. Activity data are often required at relatively short intervals, perhaps as frequently as annually or biennially. For such short reporting intervals, area change may be small. For the purposes of precisely estimating area change, STR designs based on change maps commonly entail using a small change stratum with a small stratum weight, but a disproportionately large sample size; similarly, no-change strata would have large stratum weights, but with disproportionately small sample sizes. A risk associated with this approach is that only a few sample units with change observations erroneously assigned to a no-change stratum can both greatly inflate standard errors and induce a much greater range in estimates of the mean, even though the estimator remains unbiased. Buffer strata, as previously described, can be used to alleviate at least some of these adverse effects [33].

#### 3.2.4. Post-Stratified Estimators

Even though reference data may have been acquired using an SRS or SYS sampling design, some of the benefits of stratified estimation can still be realized using the post-stratified (PSTR) estimators. The primary difference between the combination of STR sampling and stratified estimation and the combination of SRS or SYS sampling and post-stratified estimation is that with the former combination, the map-based stratification guides the sampling, whereas with the latter combination, the map has no influence on the sampling. Further, whereas within-stratum sample sizes are considered fixed with the STR estimators, they are considered random with the PSTR estimators. The PSTR estimator of the class proportion is the same as for the STR estimator of Equation (8), i.e., μˆ <sup>c</sup> PSTR <sup>=</sup> <sup>μ</sup><sup>ˆ</sup> <sup>c</sup> STR. However, the PSTR variance estimator has a slightly different form to accommodate the random within-stratum sample sizes,

$$\text{SE}\{\hat{\mu}\_{\text{PSTR}}^{c}\} = \sqrt{\sum\_{\mathbf{h}=1}^{\text{h}} \left[ \mathbf{w}\_{\text{h}} \cdot \frac{\hat{\sigma}\_{\text{h}}^{2}}{\mathbf{n}} + (1 - \mathbf{w}\_{\text{h}}) \cdot \frac{\hat{\sigma}\_{\text{h}}^{2}}{\mathbf{n}^{2}} \right]^{2}} \tag{12}$$

Some researchers consider the variance and SE to be conditional on the sample, in which case the estimator of Equation (9) can be used instead of Equation (12) [34] (pp. 152–156). The same recommendations regarding the maximum number of strata and the minimum within-stratum sample sizes for stratified estimation pertain for post-stratified estimation.

#### 3.2.5. Model-Assisted Estimators

The model-assisted (MA) estimators can be used with any probability sampling designs. The first component of model-assisted estimators of the proportion, μc, is formulated as the synthetic (SYN) estimator,

$$\hat{\mu}\_{\text{SYN}}^{c} = \frac{1}{\text{N}} \sum\_{i=1}^{N} \hat{\mathbb{y}}\_{i\prime}^{c} \tag{13}$$

where N is the total number of population units and yˆ <sup>c</sup> <sup>i</sup> is a prediction for the *i*th population unit. However, the SYN estimator may be biased because of prediction error. For SRS and SYS designs, the bias can be estimated as,

$$\text{B}\hat{\text{ias}}\{\hat{\mu}\_{\text{SYN}}^{\text{c}}\} = \frac{1}{\text{n}}\sum\_{i=1}^{n} \left(\hat{\mathbf{y}}\_{i}^{\text{c}} - \mathbf{y}\_{i}^{\text{c}}\right). \tag{14}$$

The asymptotically unbiased model-assisted estimator of the proportion of the total area is then,

$$\begin{split} \boldsymbol{\upmu}\_{\text{MA}}^{\text{c}} &= \boldsymbol{\upmu}\_{\text{SYN}}^{\text{c}} - \text{Bias} \left( \boldsymbol{\upmu}\_{\text{SYN}}^{\text{c}} \right) \\ &= \frac{1}{N} \sum\_{i=1}^{N} \boldsymbol{\upmu}\_{\text{i}}^{\text{c}} - \frac{1}{n} \sum\_{i=1}^{n} \left( \mathbf{\uphat{y}}\_{\text{i}}^{\text{c}} - \mathbf{y}\_{\text{i}}^{\text{c}} \right) \end{split} \tag{15}$$

with standard error,

$$\text{SE}(\hat{\mu}\_{\text{MA}}^c) = \sqrt{\frac{1}{\text{n} \cdot (\text{n} - 1)} \cdot \sum\_{i=1}^{n} \left(\varepsilon\_i - \overline{\varepsilon}\right)^2} \tag{16}$$

where

$$
\varepsilon\_{\mathbf{i}} = \mathbf{y}\_{\mathbf{i}}^{\circledast} - \mathbf{y}\_{\mathbf{i}}^{\circledast} \text{ and } \ \overline{\varepsilon} = \frac{1}{\mathbf{n}} \sum\_{\mathbf{i}=1}^{n} \varepsilon\_{\mathbf{i}}.\tag{17}
$$

When the prediction technique used to obtain yˆ <sup>c</sup> <sup>i</sup> is formulated and calibrated using data external to the area of interest, or otherwise does not use observations of *y* from the sample, the estimator is characterized as the difference (DIF) estimator, whereas if the prediction technique is calibrated using data internal to the area of the interest, the estimator is characterized as the generalized regression (GREG) estimator [35]. For sampling designs other than SRS and SYS, Equations (14)–(16) must be revised to accommodate the features of the designs (e.g., [25]).

Representative examples, as opposed to an exhaustive list of tropical applications which are sparse, of uses of the STR, PSTR, DIF and GREG estimators for estimating land cover and/or land use area and area change are reported in Table 1. For categorical response variables, the STR and PSTR estimators are generally expected to produce greater precision than the model-assisted estimators [26,36], although Stehman [36] notes three exceptions: (i) large overall accuracies, (ii) small true proportions, and (iii) small reference sets. For continuous response variables, the model-assisted estimators are generally more precise than the STR and PSTR estimators. Exceptions are when the PSTR estimator is formulated as a model-assisted estimator [16], and possibly when the map resolution is much coarser than the resolution of the reference data [37].


**Table 1.** Representative examples of land cover class area and area change estimation.

\* STR: Stratified, PSTR: Post-stratified, DIF: Difference, GREG: Generalized regression.

#### *3.3. Time Series Analyses*

With the increasing availability of free Landsat and Sentinel 2 spectral data at regular temporal intervals over many years, opportunities for tracking land cover and land use changes over time are greatly enhanced [3]. Two applications are attracting attention. First, trajectory analyses track metrics for individual pixels over time for purposes of characterizing phenomena, such as gradual forest degradation, abrupt disturbances such as harvest followed by recovery, and dates of change [53–60]. The second application entails estimating differences or trends in forest area or forest area change at multi-year intervals, or perhaps even annually. With this approach, a remote sensing-assisted estimate of area or area change and an associated SE are calculated, as per Section 3.2, for each temporal point of interest.

For both applications, compliance with the second IPCC good practice guideline related to rigorous uncertainty assessment presents challenges. For the first application, for which dates of individual pixel-level changes are estimated, very large numbers of tests of significance will inevitably lead to large numbers of false positives and false negatives. For the second application which focuses on multi-year trends of annual differences, an intuitive approach would be to use a regression or similar prediction technique for estimating trends or an ANOVA technique followed by a simultaneous multiple inference method for determining which, if any, temporal estimates differ statistically significantly from other temporal estimates. A complicating factor, however, is that both regression and ANOVA techniques assume that the response variable, in this case forest area change, is observed with at most negligible uncertainty. For this kind of trend analysis, this assumption is not satisfied, because the standard errors associated with the individual change estimates are likely non-negligible. Special techniques similar to hybrid inference would be necessary to accommodate and account for the uncertainty in the temporal area or area change estimates [35,61–63].

#### **4. Emissions and Removals Factors**

Emissions and removals factors, by definition, represent carbon change per unit area. Estimates of these factors can be obtained from three primary sources: published default values [64], ground biomass observations that are converted to carbon, and remote sensing-based biomass maps from which estimates can be converted to carbon. For this review, the focus is remote sensing-assisted estimation of biomass change per unit area, particularly aboveground, live tree biomass change per unit area.

The state-of-the-science for estimating biomass change as a precursor to estimating emissions and removals factors is considerably less mature than for estimating activity data. To estimate biomass change, reference data for two dates are required, although an exception may be for complete deforestation, for which there is no remaining biomass at the second time. Two approaches for estimating change can be used, the indirect approach and the direct approach. With the indirect approach, total or mean biomass per unit area is estimated for an activity as the difference in two biomass estimates, one based on a sample acquired at time 1 and the second based on a sample acquired at time 2, usually for a different set of sample unit locations, because otherwise, the direct approach would be used. With the direct approach, biomass change is estimated directly, using biomass change observations as reference data. Of necessity, the direct approach requires ground level observations of biomass for the same locations at the two times, a requirement that is currently difficult to satisfy for many tropical applications. Thus, at the current time, the indirect approach is more commonly used.

For estimating either biomass or biomass change, remote sensing-based maps are used for three purposes: first, as auxiliary data to enhance design-based estimates using probability samples of ground reference data; second, for direct estimation when constructed using either probability or non-probability samples of ground reference data; and, third, as sources of reference data. Multiple biomass maps are available for the first and third purposes: a circa year 2001, 250-m, MODIS-based map for the USA [65]; a global, circa year 2005, 1-ha GLAS, ALOS, and Landsat-based map constructed by NASA's Jet Propulsion Laboratory [66]; the year 2010, 1-ha, lidar and SAR-based GlobBiomass map [67]; a circa 2007-2008, 500-m, GLAS and MODIS-based, pan-tropical map [68]; and the European Space Agency's year 2017 global datasets of AGB [69]. Global, regional, or even large area biomass change maps are not commonly available, at least partially as a result of the lack of ground biomass change observations. The emphasis of the sections that follow relates to the use of maps to estimate biomass, although the estimation of biomass change is completely analogous.

#### *4.1. Probability Samples of Ground Reference Data*

Remote sensing-based forest attribute maps, not just biomass or biomass change maps, can be used as sources of auxiliary data for increasing the precision of ground-based estimates of both biomass and biomass change. If observations of biomass or biomass change are obtained using a probability sampling design, any of the design-based STR, PSTR, or model-assisted estimators described in Section 3.2 can be used, albeit with μ<sup>c</sup> denoting mean biomass or mean biomass change rather than an area proportion.

The statistical forms of the STR and PSTR estimators are the same as described for estimating activity data. For both estimators, the primary purpose of the maps is to serve as a basis for stratifications, which, in turn, are used to increase precision. For both these estimators, the map values are aggregated into a small number of contiguous classes that serve as strata [6,35,70]. The map-based stratifications can also be used to guide sample allocations to strata.

For use with the model-assisted estimators, the primary purpose of the maps is to serve as a source of predictions. The forms of both model-assisted estimators are the same as described for estimating activity data in Equations (13)–(16). Although many of the early applications used linear regression models, additional prediction techniques have been used, including nonlinear models [36,71,72]; k-nearest neighbors [73,74], and random forests [73,75].

Representative examples, as opposed to an exhaustive list of tropical applications which are sparse, of the uses of maps as auxiliary data for the estimation of biomass and biomass change, are reported in Table 2.


**Table 2.** Representative examples of biomass and biomass change estimation.

\* STR: Stratified, PSTR: Post-stratified, DIF: Difference, GREG: Generalized regression.

#### *4.2. Non-Probability Samples of Ground Reference Data*

Although probability samples of ground data may not be available to support design-based inferential methods, ground observations and measurements may still be available from other sources, such as long-term research plots, local forest management plots, and pre-harvest plots. If conditions associated with these plots conform to the features of activities of interest such as forest-remaining-forest, pre-harvest as a form of deforestation, or thinning treatments as a form of degradation, they can be used for estimating emissions and removals factors. Two key challenges are associated with the use of data from these kinds of plots. First, differences in data features acquired from these kinds of plots must be reconciled before the ground data can be combined with remotely sensed data, and second, the lack of a consistent probability sampling design means that less familiar and more complex model-based inferential methods must be used instead of design-based inferential methods.

#### 4.2.1. Data-Related Challenges

The data-related challenges result from the different plot configurations and different measurement protocols inevitably used for the different sources of reference data. In general, smaller plots tend to have more extreme per unit area observations than larger plots that are at least marginally more representative of the entire population. Thus, combining data for mixed size plots runs the risk of skewing relationships with remotely sensed data toward the data from the smaller, less representative plots. For modeling applications, multiple studies have shown the advantages of larger plots with smaller area to perimeter ratios that minimize edge effects [78,80,81]. Although no studies evaluating the effects of constructing models using data for mixtures of small and large plots are known, the effects are expected to be a form of heteroscedasticity and possibly model predictions skewed toward the data for the smaller plots. Differences in measurement protocols include, but are not limited to, differences in minimum diameters of trees to be measured and minimum height at maturity. For otherwise comparable plot configurations, plot-level biomass would be larger for smaller minimum diameters. Moreover, plots with trees that only marginally satisfy a small minimum height criterion might not even be measured if the criterion was larger.

If data for plots with different configurations and measurement protocols are to be combined, some form of data harmonization is necessary. For example, the smallest among multiple plot radii can be selected, and data for the now smaller plot can be recalculated. Similarly, the largest among multiple minimum diameters can be selected. The issue of harmonization of national forest inventories in Europe has received considerable attention, including the development of useful harmonization

methods. Although developed for temperate forests, these methods are likely also applicable for tropical forests [82,83].

#### 4.2.2. Model-Based Inference

Model-based inference, also characterized as model-dependent inference, relies on a quite different set of underlying assumptions than the more familiar design-based inference [72]. First, validity is based on correct model specification rather than a probability sample. Second, an entire distribution of possible values is assumed for each population unit, rather than just a single value. Third, randomization is via the realized observations from the distributions characterizing the population units selected for the sample rather than the sampling design. An important consequence of the first assumption is that model-based inference does not require probability samples of reference data for constructing the model. Although probability samples may be used and, in fact, may be preferable, purposive and other non-probability samples may also produce entirely valid model-based inferences [30] (p. 534). The absence of a requirement for a probability sample means that model-based inference can be used for applications for which design-based inference is not possible, such as when a non-probability sample is used to construct the model.

The model-based estimator of the mean is simply the synthetic estimator of Equation (13). As noted in Section 3.2.5, the synthetic estimator may be biased. The model-based SE of the estimate requires three sources of uncertainty information: (1) variances and covariances among population unit predictions due to sampling variability associated with the training data used to construct the map, (2) residual uncertainty in the form of differences between the map's training data and the corresponding map unit values, and (3) spatial correlation among the residuals. Of these three sources, the effects of sampling variability are often the greatest. Examples of model-based inference relevant for greenhouse gas inventories include McRoberts et al. [72] and Saarela et al. [84].

#### *4.3. Remote Sensing-Based Maps as Sources of Reference Data*

In the absence of ground data, large area, regional, or global biomass maps can be used as sources of reference data for estimating biomass. However, because the maps consist of map unit predictions subject to uncertainty, compliance with the IPCC good practice guidelines requires special considerations. In particular, the map must be validated, if not in its entirety, then preferably for a validation unit coincident with the area of interest and/or perhaps for a sample of validation units. Validation consists of a statistically rigorous test of the hypothesis of no difference between the map-based estimate and an estimate obtained using independent (ind) data of the form,

$$\mathbf{t} = \frac{\boldsymbol{\upmu}\_{\text{map}} - \boldsymbol{\upmu}\_{\text{ind}}}{\sqrt{\text{Var}\left(\boldsymbol{\upmu}\_{\text{map}}\right) + \text{Var}\left(\boldsymbol{\upmu}\_{\text{ind}}\right)}}.\tag{18}$$

where t follows Student's t-distribution with degrees of freedom determined by the sum of the sample sizes for the two variance estimates. Although these independent data serve as reference data for this analysis, for this discussion, they continue to be characterized only as independent data to avoid confusion with the reference data to be acquired from the map. Four estimates are therefore required: (1) the estimate based on the independent data, μˆind; (2) the variance of the independent data estimate, Varˆ (μˆind); (3) the map-based estimate for the validation unit, μˆ map; and (4) the variance of the map-based estimate, Varˆ <sup>μ</sup><sup>ˆ</sup> map .

#### 4.3.1. Reference Data Estimates

McRoberts et al. [37] demonstrate two inferential approaches for acquiring the independent data necessary for the comparison of Equation (18). The first approach uses a sample of the independent data and auxiliary data to calculate an estimate for the validation unit that can be compared to the map estimate. If the independent data are obtained from a probability sample, one of the design-based estimators of Section 3.2 can be used, for which the map that is subject to validation can be used as auxiliary data. If the independent data are not obtained using a probability sampling design, the model-based estimators described in Section 4.2.2. and by McRoberts et al. [37] can be used.

The second approach uses a local map of greater quality than the biomass map as the source of independent data. Greater quality is with respect to factors such as finer resolution, use of additional auxiliary data, and use of local training data. The local map is sampled to obtain independent data using any convenient sampling intensity and sampling design, although probability sampling designs greatly simplify estimation. Because the local map also consists of predictions rather than observations, the uncertainty in the local map must be incorporated into the estimation of overall uncertainty using hybrid inference [35,61–63].

#### 4.3.2. Map-Based Estimates

The map-based mean for Equation (18) is simply the mean over all map units obtained using the SYN estimator of Equation (13). Because the data used to construct the map are typically not available, estimation of the SE for the map-based estimate can be particularly challenging and typically entails model-based inferential methods, as described in Section 4.2.2. However, unlike Section 4.2.2, where the remote sensing-based map is constructed by the user, global maps used as sources of reference data are generally not constructed by the user. Thus, information on the sources of uncertainty necessary for calculation of the SE must be provided by the map authors. Unfortunately, however, they rarely provide such information, presumably because they do not know it is necessary or they do not know how to produce it. McRoberts et al. [37] describe some approximations and bounds on the uncertainty, but the validation process would be greatly facilitated if map authors were cognizant of these requirements.

#### **5. Summary**

Remotely sensed imagery and land use change maps are crucial for the statistically rigorous estimation of activity data. The post-stratified estimators are appropriate and effective for use with land use change maps and data obtained from an existing probability ground sampling program. When such ground data are unavailable or cannot be readily obtained, the combination of land use change maps, land use change interpretations of stratified samples from imagery, and the stratified estimators have been demonstrated to be particularly efficient. When sample unit data are continuous, such as plot- or pixel-level proportions of land use change, the model-assisted estimators can be used. In general, the design-based stratified, post-stratified and model-assisted estimators are more frequently used and are easier to use than model-based estimators.

Methods for estimating emissions and removals factors are much less mature than methods for estimating activity data, likely because the former rely much more heavily on ground data which may not be available. When probability-based samples of ground data are available or can be acquired, remote sensing-based biomass maps can be used as sources of auxiliary data for increasing the precision of estimates of emissions and removals factors obtained using any of the design-based stratified, post-stratified, and model-assisted estimators.

When some ground data are available, but were not acquired using probability sampling designs and/or not with similar plot configurations and measurement protocols, model-based inferential methods can be used to estimate biomass or biomass change. If the ground data are acquired from multiple sources, a degree of data harmonization will inevitably be necessary.

When sufficient ground data are not available or cannot be acquired, biomass maps can be used in lieu of ground data. However, the required statistical methods are more complex, because the uncertainty in the maps must be accommodated to comply with the IPCC good practice guidelines for greenhouse gas inventories. These methods typically require model-based inferential methods that are difficult to implement, apart from substantial meta-data regarding sources and effects of uncertainty (Section 4.2.2) that must be provided by the map authors.

Additional research should focus in several areas: (1) the documentation of actual tropical applications, (2) the development of rigorous uncertainty assessment methods for use with time series methods, and (3) greater attention by map authors on providing meta-data for global maps to facilitate uncertainty estimation.

**Author Contributions:** Conceptualization: R.E.M., E.N., C.S., S.V.S.; writing—original draft preparation, R.E.M.; writing—review and editing: R.E.M., E.N., C.S., S.V.S., E.O.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors thank two anonymous reviewers for their constructive comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Remote Sensing* Editorial Office E-mail: remotesensing@mdpi.com www.mdpi.com/journal/remotesensing

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18