Remote Sensing-Based Biomass Estimation

Abstract

Over the past two decades, one of the research topics in which many works have been done is spatial modeling of biomass through synergies between remote sensing, forestry, and ecology. In order to identify satellite-derived indices that have correlation with forest structural parameters that are related with carbon storage inventories and forest monitoring, topics that are useful as environmental tools of public policies to focus areas with high environmental value. In this chapter, we present a review of different models of spatial distribution of biomass and resources based on remote sensing that are widely used. We present a case study that explores the capability of canopy fraction cover and digital canopy height model (DCHM) for modeling the spatial distribution of the aboveground biomass of two forests, dominated by Abies Religiosa and Pinus spp., located in Central Mexico. It also presents a comparison of different spatial models and products, in order to know the methods that achieved the highest accuracy through root-mean-square error. Lastly, this chapter provides concluding remarks on the case study and its perspectives in remote sensing-based biomass estimation.

1. Introduction

Forest ecosystems are about 31% of the total land cover of the earth [1], being one of the most important ecosystems due to economic goods and environmental services they provide. One of these services is as an environmental regulator, reducing the concentration of carbon dioxide (greenhouse gas) from the atmosphere and transforming it into oxygen and biomass through photosynthesis, thereby playing an important role in the global carbon cycle [2–4].

Biomass is defined as the dry weight of both aboveground biomass (AGB) and belowground biomass (BGB) living mass of vegetation, such as wood, bark, branches, twigs, stumps, or roots as well as dead mass of litter associated with the soil [4–6]. According to this, it can be considered as a measure of forest structure and function. Thus, by knowing the spatial distribution of biomass, it is possible to calculate the net flow of terrestrial carbon, nutrient cycling, forest productivity, biomass energy, and carbon storage and sequestration by the forest, reducing the uncertainty of carbon emission and sequestration measures to support climate change modeling studies [5–8].

Since calculating field measures of BGB is difficult, most studies have focused on calculating AGB. The most accurate way to obtain AGB data is by using field measurements and allometric equations for individual trees; however, these techniques are difficult to implement because they are time consuming and labor intensive. Furthermore, forests are a complex and widely distributed ecosystem, which makes these techniques expensive to apply in large areas; therefore, they cannot provide the spatial distribution of biomass [4,6].

An alternative form to map and monitor spatial distribution of AGB is through the use of remote sensing-based techniques, because through them it is possible to obtain a continuous and repetitive collection of digital data from the same area with different spatial resolutions, covering large areas and reducing processing time and costs [5,6]. Due to the increasing availability of satellite imagery, several researches have been developed to prove the effectiveness of both imagery data provided by different sensors and diverse modeling approaches [4,6,9]. In order to estimate AGB, two main kinds of models have been used: direct and indirect. The direct models measure biomass throughout the relationship between spectral data response and biomass field measurements. For the indirect models, biomass is estimated from biophysical parameters or forest structural metrics [10].

This chapter reviews the main models used for estimating biomass and key resources used in remote sensing (Sections 1 and 2). The case study integrates some models and resources applied in forests dominated by Abies Religiosa and Pinus spp. located in Central Mexico. In the last section, concluding remarks are provided about the best biomass estimation models as well as their limitations.

2. Biomass modeling

Several factors can affect the remote sensing-based AGB estimation, such as insufficient sample data, atmospheric conditions, complex biophysical environments, scale of the study area, availability of software, spatial resolution of remotely sensed data, or mixed pixels, among others [6, 10]. In order to introduce different approaches that have been developed to reduce the uncertainties produced by these factors on the estimation of the spatial distribution of AGB, the most commonly used models will be described in this section.

2.1. Field measurements and allometric equations

All spatial distribution AGB estimation models need high quality and representative field data to be implemented. Therefore, forest inventories are the most common approach to obtain detailed and periodic data. They are held for monitoring, modeling, and predicting several biophysical processes, such as stocking levels, harvests, diseases, and pests, among others. Therefore, they are generally implemented at several scales to obtain different structural parameters, either through the data aggregated from stand-level management inventories or by plots established through a statistical sampling design [3,4,6,11].

Since the grouping of data from stand-level management inventories tend to underestimate the forested areas and stock volume calculation, presently in most of the countries a statistical sampling design is used. In this procedure, the sampling plots are randomly selected from a population where each one of them has the same positive known probability to be chosen. They can be selected by different methods. One of them is a systematic sampling procedure based on the use of grids of randomly selected points in two dimensions, with a 0.5–20-km separation range. In other cases, plot locations are randomly selected by regular polygons created by a tessellation of large areas, which can be stratified when different sampling intensities are required, for example, when different land uses and covers exist [11].

In order to make a more efficient sampling, a plot generally conforms clusters, commonly of four plots but, in some cases, they can be as large as 18, with a huge variety of shapes and sizes. Circular plots are commonly used in Boreal and Temperate forests, whereas square and rectangular ones are used in rainforests. Plots can be combined with other sampling methods like transect for measuring deadwood or soil pits to assess soil carbon [11].

Commonly, the information gathered about a plot is location, number of trees, species, health, and site description, among others. In addition, individual tree dendrometric variables are considered, such as diameter at breast height (DBH), tree height, crown size, and canopy cover, DBH and tree height being the most commonly used parameters to derivate AGB through allometric equations [3,4,10].

The allometric biomass equation is the most common and accurate method to translate forest inventory reports of individual tree data to tree and stand biomass. It is a mathematical relation between total tree biomass (stem, branch, and foliage) and its DBH or both DBH and height, applying a least-squares regression of logarithmic equation. It can be both species- or site specific (e.g., Pinus montezume or Abies religiosa) or more generic (e.g., pine gender or tropical hardwoods); however, it has been observed that biomass equations at the site-specific level produce better results than generic equations [3,8,12].

The most accurate method to obtain allometric equations is a destructive process in which individual trees from a wide range of DBH, distributed in a local forest, must be felled and separated into boles, branches, and leaves. Then boles are cut into sections and weighted in the field as leaves and branches. After that, a thick disc sample must be cut from the base of each bole section and a subsample is extracted for the other side of the disk, both being dried at 105°C until a constant weight is reached to obtain the dry mass, to estimate the moisture content and the wood density. The total biomass therefore is the sum of the dry mass of the branches, the leaves, and the various sections of the stem [3,13].

Since it has been observed that growing plants maintain the weight proportion between different parts, it is possible to build a nonlinear mathematical model to relate biomass with DBH using these parameters [12]. One of the most common models used is the logarithmic equation (1):

$\text{Ln}\left(B\right)=\text{Ln}\left(a\right)+b\text{Ln}\left(D\right)$E1

where B is the biomass, a and b are scalar coefficients estimated by least-squares linear regression, and D is the DBH.

Using allometric equations, it is possible to compute the total AGB for a given area using biomass expansion factors or conversion tables; however, these approaches do not provide the spatial distribution of AGB [3,10]. Later in this chapter, different models created to obtain spatial distribution of AGB are explained.

2.2. Regression models

One of the most common methods to estimate biomass is the regression analysis, which is a statistical technique to investigate and model the relationship between variables. Traditionally, in the remote sensing approach, the regression analysis techniques applied to AGB estimation are based on the quantitative relationship between ground-based data and satellite information, such as spectral reflectance, radar, or light detection and ranging (LiDAR) data [6,14,15]. Models based on regression analysis are considered to be relatively easy to implement and can provide accurate results through their application at all spatial scales [16]. Generally, these methodologies consist of three major steps: biomass estimation based on fieldwork, establishment of regression model between field biomass and satellite information of corresponding pixels, and the use of regression models to generate a biomass image with the spatial prediction.

These remote sensing-based biomass estimation methods assume that the forest information, obtained by the sensors, is highly correlated with AGB. According to this, the keys for biomass estimation are the use of appropriate variables and the development of suitable estimation models for sufficient sample plots, using regression methods that aim for an efficient integration of multisource data, necessary to get better biomass estimation. Spectral data, radar, and LiDAR have their own positive and negative characteristics and proper integration of them can improve biomass estimation accuracy [17].

After data integration, the correct use of regression methods for establishing biomass estimation models is also important. Many models have been developed based on multiple combinations of in situ tree parameters calculated through linear regression (LR) or nonlinear regression (NLR) models [18–20]. Multiple regression analysis may be the most frequently used approach for developing biomass estimation models [21–23]. In both cases, these parametric algorithms assume that the relationships between dependent (e.g., biomass) and independent (derived from remote sensing data) variables have explicit model structures that can be specified a priori by parameters [15]. Generally, the independent variables can be spectral bands, vegetation indices, textural images, LiDAR height, and synthetic aperture radar (SAR) backscatter; in some cases, the use of subpixel information offered better estimation results than per-pixel-based spectral signatures [24].

However, the linear regression approach has been known to mislead the prediction of the studied variable at values beyond a saturation point of the canopy reflectance [25]. Since biomass is usually nonlinearly related to remote sensing variables, nonlinear models such as power models [26], logistic regression models [27], and geographically weighted regression models [16] were often used to estimate biomass with more accuracy. Nonetheless, some estimation methods have been established as a nonparametric alternative to the use of regression approaches for biomass modeling: k-nearest neighbor (k-NN), artificial neural network (ANN), regression tree, random forest, support vector machine (SVM), and maximum entropy (MaxEnt) [15].

2.3. Geostatistical models

Some studies have used geostatistics as the main approach for biomass forest estimation in order to predict variables related to forest structure (e.g., tree height and volume) and aboveground biomass and carbon measurements in unsampled sites based on known values of adjacent spatial data as forest inventory sites [28,29]. Other recent works have explored the synergy between geostatistical models with remotely sensed data to improve estimations using remote sensing indices as spatial secondary variables [10, 16, 30–32].

Geostatistics was defined in the 1960s by Georges Matheron, who generalized a set of techniques developed by Krige (1951) in order to exploit the spatial correlation to make predictions in evaluating reserves of gold mines in South Africa. This generalization is detailed in his regionalized variable theory in 1970 [33]. The purpose of geostatistics is the estimation, prediction, and simulation of the values of a variable that is distributed through space [34]. This theory assumes that a variable measured in a spatial domain corresponds to a random variable z(x), assuming that the structure of the phenomenon having spatial correlation is considered a regionalized variable; therefore, a set of spatially distributed random variables will be a random function Z(x). This provides the theoretical basis for establishing the spatial structural characteristics of natural phenomena. Moreover, it can be used as a tool for calculating the value of a variable in a certain position in space, knowing the values of that variable among adjacent positions in space, which is known as interpolation [35].

There are diverse methods of interpolations, which can be classified into two main groups: deterministic and geostatistical. The deterministic techniques are based directly on some properties of similarity of adjacent measured values (e.g., distance), which establish a set of mathematical formulas that determine the smoothness of the resultant surface interpolated. Examples of these are the inverse distance weighting (IDW), nearest neighbors, splines, and triangular irregular network (TIN) [35]. Geostatistical techniques studied spatial autocorrelation of the variables in order to fit a spatial dependence model to a set of random variables. This approach produces predictions and also generates an error surface concerning the uncertainty-associated analyzed model [33].

The spatial dependence is the spatial behavior of a phenomenon, derived from spatial patterns in terms of distance and similarity and/or contrast of a spatial unit or relative spatial location to other spatial units [36]. It has its basis in Tobler’s first law of geography proposed in 1970 that says, "everything is related to everything else, but near things are more related than distant things" [37].

The kriging algorithms are one such example, which is mostly used in geosciences, ecology, and geomatics [35]. Kriging is a generic name for a family of generalized least-squares regression techniques, where the spatial structural characteristics are accomplished by the semivariogram function as a metric of the spatial autocorrelation [33].

All kriging estimators are variants of the following basic equation (2):

$\widehat{Z}\left(x_0\right)-\mu ={\displaystyle \sum }_{i=1}^n{\lambda }_i\left[Z\left(x_i\right)-\mu \left(x_0\right)\right]$E2

where µ is a known stationary mean, assumed to be constant over the whole domain and calculated as the average of the data. The parameter λ is the kriging weight; n is the number of sampled data points used to make the estimation, and µ(x₀) is the mean of the samples within the search window [33].

Different studies have applied univariate and bivariate geostatistical interpolations in order to calculate the forest volume [28], aboveground biomass [10,16,29–31], and carbon in the aboveground biomass [32]. The most commonly used technique for univariate-based modeling is kriging [28], whereas in bivariate-based modeling regression-kriging [10,30–32], cokriging [31], kriging with external drift [29], cokriging regression [31], and geographically weight regression (GWR) are used [16].

2.4. Nonparametric models

Similar to regression models, nonparametric algorithms are based on the use of different sensor data, for example, spectral, radar, and LiDAR [21,38], using many of these models in the forest attributes estimation [39–41]. They are a framework for creating complex nonlinear biomass models based on the use of remote sensing variables and as alternatives for the parametric approaches. Common nonparametric algorithms include k-nearest neighbor (k-NN), artificial neural network (ANN), random forest, support vector machine (SVM), and maximum entropy (Max Ent).

One of the most applied nonparametric methods is the nearest neighbor approach (NN). In the context of forest attribute estimation, the NN methods have been first introduced in the late 1980s [42]. In the NN methods, the value of the target variable at a certain location is predicted as a weighted average of the values of neighboring observations, with k-nearest neighbors spectrally using a weighting method [43]. Several methods have been offered to measure the distance from the target unit to the neighbors. In the NN approach, the choice of the k value, type of distance measure including Euclidean and Mahalanobis, along with weighed function are the critical factors influencing the estimation accuracy [44,45].

The ANN provides a more robust solution for complicated and nonlinear problems due to its properties [46]. The network commonly consists of one input layer, one or more hidden layers, and one output layer. Since it does not require the assumption that data have normal distribution and linear relationships between biomass and independent variables, the ANN can deal with different data through approximation, using various complex mathematical functions, with independent variables from different data sources such as remote sensing and ancillary data [15]. A detailed overview of ANN approach is provided in Ref. 47.

Regression tree and random forest are a family of tree-based models; in the first one, data are stratified into homogeneous subsets by decreasing the within-class entropy, whereas in the second one, a large number of regression trees are constructed by selecting random bootstrap samples from the discrete or continuous dataset. In fact, the random forest algorithm is now widely used for biomass estimation [48,49].

SVM is an important method to estimate forest biophysical parameters using remote sensing data [50, 51]. It is a statistical learning algorithm with the ability to use small training sample data to produce relatively higher estimation accuracy than other approaches like ANN [15]. Ref. 51 provides a detailed overview of the SVM approach used in remote sensing. The Max Ent approach is a general-purpose machine-learning method for predicting or inferring target probability distribution from incomplete information [52]. These kinds of nonparametric algorithms have become popular in biomass modeling when large representative field datasets exist for calibration [53,54].

3. Remote sensing products for biomass modeling

3.1. From optical sources

Optical sensors are those that detect electromagnetic radiation emitted or reflected from the earth, the main source of light being the sun. Among the passive sensors are photographic and optical-electronic sensors that combine similar photographic optics and electronic detection system (detectors and push scanning) and image spectrometers [55,56]. Optical remote sensing refers to methods and technologies that acquire information from the visible, near-infrared, shortwave infrared, and thermal infrared regions of the electromagnetic spectrum. They are called optical, because energy is directed through optical components such as lenses and mirrors.

3.2. Spectral indices

Remotely sensed spectral vegetation indices represent an integrative measure of both vegetation photosynthetic activity and canopy structural variation that are widely used and have benefited numerous disciplines interested in the assessment of biomass estimation [57,58]. Different kinds of vegetation indices have been in use for a long time in AGB estimation and the number of publications is immense [6].

The key to the development of vegetation indices is the ability of the canopy of green vegetation to interact differently with certain portions of the electromagnetic spectrum. Since this contrast is particularly strong between red and near-infrared regions (NIR), it has been the focus of a large variety of attempts to build up quantitative indices of vegetation condition using remotely sensed imagery [59]. Theoretically, the ideal vegetation indices should be particularly sensitive to different vegetation covers, insensitive to soil brightness and color, and little affected by atmospheric effects [60]. However, in reality, different factors affect the reflectance of vegetation and consequently the vegetation index (e.g., atmospheric correction is essential when biomass is extracted from the vegetation indices as a final product).

According to Ref. 58, vegetation indices can be classified into (1) slope-based, (2) distance-based, and (3) transformation indices. The slope-based indices are simple arithmetic combinations that focus on the contrast between the spectral responses patterns of vegetation in the red and near-infrared portions of the electromagnetic spectrum. The most known of them are the ratio vegetation index (RVI) proposed by Birth and Mc Vey [61], normalized difference vegetation index (NDVI) introduced by Rouse et al. [62], and soil-adjusted vegetation index (SAVI) developed by Huete [63]. The vegetation indices show better sensitivity than individual spectral bands for the detection of biomass [64].

Distance-based indices measure the degree of vegetation from the soil background (known as a soil line) to the pixel with the highest content of vegetation in a perpendicular incremental distance. In this group of indices, the slope and intercept of the soil line have to be defined for each particular image; perpendicular vegetation index (PVI) introduced by Richardson and Weigand [65] cancels the effect of soil brightness in cases where vegetation is sparse and the pixels contain a mixture of green vegetation and soil background. The effect of the background soil is a major limiting factor in certain statistical analyses geared toward the quantitative assessment of AGB [59].

Transformation indices are transformations of the available spectral bands to form a new set of uncorrelated bands within which a vegetation index band can be defined. Tasseled cap (TC) is one of the most widely used indexes of this type and may apply to various remote sensing images with multiple resolutions [66–72].

In the specialized literature, many vegetation indices have been proposed, and depending on the complexity of the forest stand structure, indices vary in their relationships with biomass [23,60,73]. For forest sites with complex stand structures, vegetation indices including near-infrared wavelength have weaker relationships with biomass than those including shortwave infrared wavelength. In contrast, for areas with poor soil conditions and relatively simple forest stand structure, near-infrared vegetation indices had a strong relationship with biomass, and finally the results of transformation indices showed stronger relationships with biomass independent of different biophysical conditions [15].

3.3. Textural indices

Some studies have used regression analyses between remote sensing textural indices and biomass data from sampling sites [30,32,74–77]. This framework has been applied to different optical and synthetic aperture radar-derived indices in order to use the textural parameters as continuous spatial variables to improve biomass estimations.

In 1973, Haralick proposed a statistical analysis based on a set of parameters according to spatial dependence of gray tones in an image, defined as second-order statistics [78]. Textures are intrinsic properties of surfaces and its importance lies in image–objects segmentation, because they are related to structural arrangements of the land surface and the connections between neighboring spatial objects [79]. The most common mathematical method used to measure texture parameters is the co-occurrence matrix (gray-level co-occurrence matrix) [78]. It describes the frequency of a gray level displayed, in a specific spatial relationship, to another gray value within a neighborhood represented by a mobile window or kernel [79]. Its construction is based on four steps: (1) Window or kernel size definition, (2) band selection input, (3) texture parameters selection, and (4) spatial dependency criteria.

Textures applied parameters are described below:

Contrast; ${\displaystyle \sum }_{i,j=0}^{N-1}{P}_{i,j}{\left(i-j\right)}^2$	1. Dissimilarity; ${\displaystyle \sum }_{i,j=0}^{N-1}{P}_{i,j}\left|i-j\right|$
2. Homogeneity; ${\displaystyle \sum }_{\text{i},\text{j}=0}^{\text{N}-1}\frac{{\text{P}}_{\text{i},\text{j}}}{1+{\left(\text{i}-\text{j}\right)}^2}$	3. Angular Second Moment; ${\displaystyle \sum }_{\text{i},\text{j}=0}^{\text{N}-1}{\text{P}}_{\text{i},{\text{j}}^2}$
4. Maximum Probability; Largest Pij value found; within the window	5. Entropy; ${\displaystyle \sum }_{i,j=0}^{N-1}{P}_{i,j}\left(-ln{P}_{i,j}\right)$
6. Mean; ${\mu }_i={\displaystyle \sum }_{i,j=0}^{N-1}i\left(P_{i,j}\right)$; ${\mu }_j={\displaystyle \sum }_{i,j=0}^{N-1}j\left(P_{i,j}\right)$	7. Correlation; ${\displaystyle \sum }_{i,j=0}^{N-1}{P}_{i,j}\left[\frac{\left(i-{\mu }_i\right)\left(j-{\mu }_j\right)}{\sqrt{\left({\sigma }_i^2\right)\left({\sigma }_j^2\right)}}\right]$
8. Variance; ${\sigma }_i^2={\displaystyle \sum }_{i,j=0}^{N-1}{P}_{i,j}{\left(i-{\mu }_j\right)}^2$; ${\sigma }_j^2={\displaystyle \sum }_{i,j=0}^{N-1}{P}_{i,j}{\left(j-{\mu }_j\right)}^2$

Table 1.

Textures

Pij = Pr (Is = i ∩ It = j) = pixel probability of s is i and t is j, for separated pixels through one pixel distance in a relative displacement vector between neighboring pixels.

These texture parameters and the indices derived (e.g., ratios) have been applied in different satellite inputs, for example, on high spatial resolution optical-infrared images as SPOT-V images [30,32] and ALOS AVNIR-2 [74], medium spatial resolution resources as Landsat TM [76] and Landsat OLI [76], and SAR images as Jers-1 [76] and ALOS Palsar [32,75].

3.4. Biophysical variables

A forest biophysical parameter is a measure that simplifies the aboveground organization of plant materials [4]. In this context, several studies in remote sensing field have focused on determining canopy structural parameters such as leaf area index (LAI), canopy height or canopy fraction [80–86], the second one being the most commonly used for biomass estimation. Because canopy height is most commonly derived from active sensors [87], this part of the chapter is focused on canopy fraction.

The reflectance of forest canopy cover recorded by the instantaneous field of view (IFOV) of the sensor is a spectral mixture obtained from the interaction between electromagnetic radiation and both canopy and forest elements such as nonphotosynthetic vegetation (NPV; such as branches, stem, and litter), photosynthetic vegetation (PV; such as leaves), and others, such as bare soil and shadow. Therefore, image pixels are generally composed of more than one element, making the image interpretation difficult, which can result in a poor relation between AGB and spectral bands [88–91].

One of the most widely used remote sensing approaches to derive and extract fraction covers from mixed pixels is the spectral mixture analysis (SMA) [92]. It decomposes the mixed pixel using a collection of constituent spectra (end members) to obtain their areal proportions or abundances in a pixel, and therefore unmixing a multispectral image into fraction images of end members [88,93]; it can be linear or nonlinear. The linear model assumes a single interaction between each incident photon and the surface object, and therefore the mixed pixel is a linear combination of pure spectral signatures (end members) of the surface materials weighted by their area covered. The nonlinear model is the opposite of linear model, since electromagnetic radiation can intercept more than one element of surface, with mixed pixels resulting from a multiple-scattered signal [94].

As green leaves scatter radiation at NIR spectra, vertical structure of vegetation commonly produces a multiple-scattered signal; however, nonlinear spectral mixture approaches are barely used because they require more specific information, such as scattering properties of end members, the illumination of the sensor, and certain geometrical parameters of the scene [91,93]. For these reasons, linear approaches have been largely implemented [92].

The linear mixture model is expressed in matrix form in the following equation (3):

$p=Mf+\epsilon$E3

where p = [p1...pn]T is the mixed pixel, f = [f1...fm]T is the fractional end-member abundance, M is an $n\times m$ matrix with n end-member spectra as column vectors, and ε is the residual not explained by the model. If p and M are known, it is possible to estimate it from ordinary least-squares procedure [93,95]), with two common constraints: the full additive condition, which determines that the fraction must sum to one, and the nonnegative condition, which makes all abundances nonnegative, thus making the model physically meaningful [92].

SMA has been successfully applied in vegetation studies to estimate the land cover fractions of PV, NPV, bare soil, and shade [96–98] or mapping the fractional cover of coniferous species [99–101]. In biomass studies, it has been proved that using ASTER fraction images (green vegetation, soil, and shade) in regression models improve the AGB estimation in Mediterranean forests, being better than NDVI or tasseled cap components [88]. Ref. [15] uses the same Landsat Thematic Mapper (TM) fractions and TM spectral signatures to relate with AGB, finding that fractional images perform better for successional forest biomass estimation than TM spectral signatures, but not for primary biomass estimation. SMA has also been applied to remove subpixel atmospheric and soil reflectance contamination in order to improve dry biomass estimations, showing that unmixed vegetation indices are better [102] than those which are not unmixed. Through geometric-optical reflectance models, Peddle et al. also estimated areal fractions of sunlit canopy, sunlit background, and shadow at subpixel scales showing higher accuracy than NDVI [103]. The other case in which sunlit crowns, background, and shadow fractions were compared with seven different vegetation indices (NDVI, SR, MSR, RDVI, WDVI, GEMI, and NLI) and three different soil-adjusted vegetation indices (SAVI, SAVI-1, and SAVI-2) to estimate biomass, LAI, net primary productivity (NPP), DBH, stem density, and basal area was the study conducted by Peddle et al. In this study, the authors concluded that the SMA shadow fraction improves the results by about 20% compared to vegetation indices [104].

3.5. From active sources

Active sensors are those that provide its own energy source in order to control the double operation of signal emission and reception of known characteristics. These sensors have the advantage of an operational capacity to produce information both at night and in the day, in addition to working in a region of the electromagnetic spectrum that makes them less sensitive to atmospheric conditions. Of these, radar and LiDAR systems [55,56] are the most known. In this section, we briefly describe each of them, pointing out the relevant examples of their application in biomass modeling.

3.6. Radio detection and ranging

RAdio detection and ranging (RADAR) is the system name of active sensors that work in microwave region of the electromagnetic spectrum. Their mechanism is performed through signal transmission–reception of a portion of energy that interacts with the surface, which is referred as backscattered, being a measure of strength and time delay of the returned signals [105].

Such energy is considered consistent or coherent (illumination beam has same wavelength and phase), which makes it possible to use different polarization schemes (orientation of emitted and detected electromagnetic fields) to generate images [106]. The spatial resolution of radar images is strongly dependent on the antenna length (aperture) of the receptor and sensor inclination angle.

Synthetic aperture radar is widely used in forest monitoring through remote sensing, which is able to generate high-resolution imagery by taking advantage of the movement of the aircraft or satellite platform. SAR simulates a long virtual antenna that comprises long coherent successive radar signals, transmitted and received by a small antenna, which simultaneously moves along a given flight path [105,106].

Since microwave energy can penetrate forest canopies, the backscattered energy of SAR systems is modulated or influenced by the structural parameters of trees (e.g., branches, leaves, and stems), which in turn depends on different ecological variables [6,107,108]. Analyses of these data have been used to determine forest state [109], forest types [110], biomass density [32,111], forest canopy height [112], forest fire degradation [49,69], deforestation [113], and forest soil moisture [114].

The sensor sensitivity to forest parameters is a function of the wavelength, for instance, bands X and C with wavelengths 2.4 and 7.5 cm, respectively, are more sensitive to backscatter of leaves, whereas bands L and P with wavelengths 15 and 100 cm, respectively, are associated with backscatter of branches and stems [55,106,115].

Three of the main approaches from SAR systems that are widely used include (1) SAR backscattering coefficient [108,111,116], (2) interferometric SAR data [32,84], and (3) polarimetric SAR data [49,117].

3.7. SAR backscattering coefficient

Forest components scatter energy transmitted by the SAR systems in all directions. A portion of energy recorded by radar is translated to a proportional ratio between density of energy scattered and density of energy transmitted from the surface targets per unit area. Backscatter coefficient (σ°) or sigma nought is the amount of radar cross section [106,108]. Generally, this magnitude is expressed as a logarithm through decibel units as the following linear-form equation (4):

$\sigma {°}_{\left(db\right)}=10{\text{*log}}_{10}\sigma °$E4

The backscatter coefficient value is related to two variables of sensor and target parameters. Sensor characteristics are a function of wavelength, polarization, and incidence angle, whereas target characteristics are associated with roughness, geometry, and dielectric properties [106,108].

Biomass modeling through this approach has been usually applied from simple regression models under the assumption of correlation between backscatter coefficient and aboveground biomass/carbon [32,108,111,116,118]. The results are different because they rely on saturation of the signal, which is a function of wavelength, polarization, and the characteristics of the vegetation cover as well as of the difficulties caused by the specific properties of the ground as slope and aspect.

Recently, some of them have combined spatial models with remotely sensed data to improve geostatistical estimations using backscatter coefficient as spatial secondary variables [32].

3.8. Interferometry SAR data

Interferometric synthetic aperture radar (Figure 1) is a framework containing diverse methods or techniques that use phase information derived by recording phase difference or state of vibration of the wave at the instant that is received by the radar between two SAR images (known as master and slave) acquired from different sensor positions [106], called the interferometric phase (ΔΦInt). The interferometric phase can be written as equation (5):

Figure 1.

Interferometric synthetic aperture radar.

$\text{Δ}\mathrm{\Phi }\text{INT}=\mathrm{\Phi }S–\mathrm{\Phi }M=\mathrm{\Phi }\text{Topo }+\mathrm{\Phi }\text{Mov }+\mathrm{\Phi }\text{Atm }+\mathrm{\Phi }\text{Noise}$E5

where ΦS and ΦM are the phase observations of slave image and phase observations of master image, respectively, ΦTopo is the topographic component, ΦMov is the shift component, ΦAtm is the atmospheric component, and ΦNoise is the phase noise.

One parameter obtained by this approach is the coherence or correlation image as an indicator of the phase stability [105]. Its mathematical expression is (6)

$\gamma =\frac{\left(S_1{S}_2^{\text{*}}\right)}{\sqrt{\left(S_1{S}_1^{\text{*}}\right)\left(S_2{S}_2^{\text{*}}\right)}}$E6

where $\gamma$ is the coherence, the brackets () denote the ensemble average, and * denotes the complex conjugate; S₁ denotes the complex value of the master single-look complex (SLC) for pixels x,y and $S_2^{*}$ indicates the complex conjugate value for pixels x,y in the second image known as slave. Coherence values are included between 0 and 1 and a coherence image quantifies the decorrelation between two SLC images as the loss of coherence [106]. Decorrelation is the combination of impacts in the radar phase: (1) the baseline decorrelation due to changes in the acquisition geometry of the images (which increases as the distance between the satellite orbits at the moment of acquisition increases) and (2) the temporal decorrelation due to variations in reflectivity in the zone, which can be caused by rain, phonological changes, and agricultural factors [105].

The interferometric coherence in biomass modeling is used under the assumption that for forested areas, coherence diminishes with the increment in vegetation density, as the volumetric scattering increases with movement (wind) and forest growth. Biomass modeling through this approach has been typically used from simple regression models under the assumption of correlation between interferometric coherence and aboveground biomass/carbon [32,84,116], another approach is by combining methods such as regression-kriging [32] and classification algorithms [119]. In this case, results are related to baseline and temporal decorrelations, forest type, polarization, and sensor wavelength.

3.9. Polarimetric SAR data

Antennas of radar systems can be configured to transmit and receive electromagnetic radiation polarized either horizontally or vertically. The two most common bases of polarizations are horizontal linear or H, and vertical linear or V. When the energy transmitted is polarized in the same direction as the received, it is called as like-polarized and when the transmitted energy depolarizes in a direction orthogonal to the received system, it is known as cross-polarization [120]. The polarization schemes are HH (for horizontal transmit and horizontal receive), VV (for vertical transmit and vertical receive), HV (for horizontal transmit and vertical receive), and VH (for vertical transmit and horizontal receive).

A radar system can be configured in different levels of polarization complexity:

• Single polarized – HH, VV, HV, or VH

• Dual polarized – HH and HV, VV and VH, or HH and VV

• Quad polarized – HH, VV, HV, and VH

A quadrature polarization or polarimetric radar uses these four polarizations in order to measure the magnitudes and relative phase difference between the polarization schemes or channels through an ellipse shape [106,120]. These kinds of radar systems promoted a new framework called polarimetry of synthetic aperture radar, which describes the surface through different combinations of polarization under the assumption that the interaction of electromagnetic energy and elements of the land surface can change the polarization of a portion of the wavelength transmitted by the sensor, and therefore receive information of amplitude and relative phase of the same target in four channels of information, which is considered as a basis for description of scattering polarimetric of surface targets [106,120]. It is mathematically simplified in the so-called scattering matrix (equation (7)):

$S=\left[\begin{array}{cc}{S}_{hh}& S_{hv}\\ S_{vh}& S_{vv}\end{array}\right]$E7

which describes different forms of the polarized electric fields between incident wave and the scatter wave in order to be the basis for diverse ways to analyze the scattering properties of a target (e.g., the covariance and coherency matrices) and diverse transformations as polarimetric decompositions and through the synergy between polarimetry and interferometry (polarimetric interferometric coherence) [120].

The use of polarimetry in biomass modeling is under the assumption of correlation between forest structural properties and polarimetric behavior. It may be construed through scattering mechanism analysis. This has been addressed mainly by polarimetric decompositions, such as Freeman Durden [49,121], eigenvector–eigenvalue [49], and Cloude and Pottier [31,121]. Biomass modeling through this framework has been usually performed from simple and multiple regression models [31,121] and nonparametric model random forest regression [49]. In this framework, results are related to forest type, spatial resolution, and sensor wavelength. 3.10. Light detection and ranging

Light detection and ranging is an active laser sensor, which emits pulses of polarized light or pulse echo, which can be calibrated within a narrow range of wavelength. The most commonly used wavelength is 1,064 nm (near-infrared), although it can range from ultraviolet to infrared range of electromagnetic spectrum (500–1500 nm) [122–124].

These laser scanners consist of a range finder, global positioning system (GPS), inertial measuring unit (IMU), and a clock capable of recording travel times to within 0.2 of a nanosecond. The integration of these systems produce accurate measurements of the position and orientation of objects registered. These technologies allow us to measure elapsed time of pulse echo between laser transmitter and objects on the surface. The energy that interacts with surfaces is backscattered over different times exhibiting multiple laser pulse returns associated with distinct surface layers toward a laser scanner that can be mounted on an aerial or terrestrial platform [123,125,126].

LiDAR information is essentially a three-dimensional point cloud composed of simple derivative returns and multiple laser pulses, this type of LiDAR data is known as discrete LiDAR returns. In addition to three-dimensional information, most LiDAR systems record the intensity as a fraction of pulse energy reflected at that location [127].

The use of LiDAR in forest areas is mainly to analyze forest vertical structural metrics under the assumption that laser can be sensed from the top of the canopies, elements of different canopies, or even to the ground, which will be reflected in the number of returns. The depth of laser penetration depends on the density of canopies and density of point clouds, which vary from less than one point per square meter to several dozens, with vertical accuracies around 12.5 cm [123,127].

One of the most widely used products in forest analysis is the result of the processing of three-dimensional point cloud, the canopy height model (CHM) (Figure 2). It is derived from the difference or subtraction between digital elevation models and digital terrain models, both datasets generally are a result of different interpolation methods, such as nearest neighbor, splines, inverse distance weighting, and kriging [127]. The first one is associated with first returns and the second one is related to the last returns. Other forest structural measurements are the fractional crown cover, crown area, crown diameter, basal area [38,125], and canopy volume [126], which are of key interest to the managers and represent information that is expensive and time consuming to collect in the field. When small-footprint LiDAR data are acquired at very high enough densities, individual tree crowns can readily be observed in the point clouds, processing algorithm for automated measuring and modeling of vegetation at individual tree crown segmentation (e.g., watershed segmentation) [95,128,129].

Figure 2.

Canopy height model.

Biomass modeling through this approach has been usually used in simple and multiple regression models [38,125,126,130]. Other works have explored the use of learning machines [131,132]. Other approach is through combining methods that integrate LiDAR information with other sensors [22,31,124,133].

4. Case study

Mexico City is a continuous truss of multiple ecosystems, which is administratively divided into two large areas: urban land (41%) governed by the Urban Development Programs and Environmental Conservation Zone (ECZ; 61%) steered by the General Ecological Planning Program (Figure 3). The ECZ provides Mexico City with environmental services such as carbon capture, aquifer recharge, biodiversity, and scenic beauty. The zone is under anthropogenic pressure, including human settlements, land use changes, and extraction of natural resources, and therefore immediate action for conservation and appropriate resource management is necessary. This has led to deforestation, degradation, development of pest infestations, fires, and erosion. Models of the spatial variability of forest density are required in order to obtain an inventory of carbon storage, useful for public policies in areas with high environmental value in order to facilitate decision-making by reducing the complexity involved in integrating and interpreting values at a pixel level.

Figure 3.

Environmental Conservation Zone.

The study area lies within the ECZ of Mexico City (882 km²) and is covered by sacred fir or oyamel (Abies religiosa) and Mexican mountain pine (Pinus hartwegii) forests. Fir forests generally grow at 2,700–3,500 m above sea level. They are evergreen forests with heights of 20–40 m and their understory is densely shaded. This type of forest contains important elements, including alders (Alnusfir mifolia), white cedar (Cupressus lindleyi), oak (Quercus laurina), Mexican Douglas-fir (Pseudotsuga macrolepis), willows (Salix oxylepis), and black cherry (Prunus serotina). Pine forests (Pinus hartwegii) grow above 3,000 m, this species being most tolerant to the extreme environmental conditions of the high mountains. This pine develops along with Festuca and Muhlenbergia grasses [134].

The present case study is a comparative analysis between regression-cokriging and multiple regression approaches using satellite-derived indices for modeling the aboveground biomass of forests in the Environmental Conservation Zone of Mexico City. The objectives are:

1. to identify satellite-derived indices that are better associated with aboveground biomass, either from LiDAR or fraction cover (SPOT-5) imagery

2. to quantify spatial patterns of the residuals derived from simple regression between satellite indices and carbon values, using spatial autocorrelation

3. to determine whether spatial statistical methods improve the estimates of aboveground biomass carbon pools over nonspatial conventional regression methods

In order to achieve these objectives, a correlation analysis was performed between digital canopy height model (LiDAR data) and vegetation fraction cover (SPOT-5 data) and, on the other hand, ground biomass estimates at forest inventory sites. Then, the spatial autocorrelation was calculated for residuals in order to define the variables to be used in multiple models and regression-cokriging methods. Once models were obtained, the root mean square error (RMSE) was computed for each approach.

4.1. Methods

4.1.1. Models of aboveground biomass carbon used three sources of data: Forest inventory data from in situ measurements

Since aboveground carbon is the amount of carbon stored in aboveground biomass, comprising all living plant material above the soil (e.g., trunks, branches, and leaves) [3], the calculation of carbon stock from biomass consists of multiplying the total biomass by a conversion factor that represents the average carbon content in biomass. A common assumption is that biomass is around 50% carbon expressed in tons of dry matter per unit area [135]. Typically, the terms of measurement are density of biomass expressed as mass per unit area (e.g., t/ha). Here, forest inventory data were obtained from Mexico City Environmental and Land Planning Authority (PAOT, because of its name in Spanish) and were derived from sampling 283 plots during 2008–2010. Their sampling is based on the design of the National Forest and Soil Inventory of the National Forest Commission (CONAFOR, because of its name in Spanish). In it, each sampling conglomerate is composed of four circular secondary sampling plots in an inverted “Y” shape, each of which covers an area of 400 m² and peripheral plots are at 45.15 m from the center of the conglomerate (Figure 4). Of the 283 plots, 155 were among pines, 86 in fir forest, 30 in mixed forest, 10 in scrub, and 2 in planted forests. Per-tree carbon was estimated from allometric carbon equations developed by Acosta-Mireles et al., Jiménez, and Avendaño-Hernández et al. for the species of the region [136–138]. Conversion of biomass carbon from conglomerate to hectares [139] used the “ratio of means” as shown in equation (8):

$\stackrel{⌢}{R}=\frac{\overline{Y}}{\overline{X}}=\frac{{{\displaystyle \sum }}_{i=1}^n{Y}_i}{{{\displaystyle \sum }}_{i=1}^n{X}_i}$E8

where Y_i is the total aboveground carbon in all plots of 400 m² and X_i is the total area sampled in i plots.

Figure 4.

Sampling plots and conglomerate form.

4.1.1.1. SPOT image

Four multispectral SPOT-5 HRG images were used: two from February 25, 2010 (zenith 51.72° and 52.03°, azimuth 136.96° and 136.43°) and two from March 28, 2010 (zenith 62.67° and 62.89°, azimuth 125.66° and 124.75°). These were radiometric and atmospherically corrected with the second simulation of a satellite signal in the solar spectrum (S-6) code through CLASlite software and orthorectified with the polynomial coefficients and geoid information based on Geocover 2000 of Landsat as reference images.

4.1.2. LiDAR data

LiDAR data used were acquired from the ALS50-II sensor flown by the National Institute for Statistics and Geography (INEGI, because of its name in Spanish) between November and December of 2007 over the entire Mexico Valley. The data had an average horizontal distance of 2.0 m, minimum point density of 0.433 points/m², and vertical root mean square error of 7.3 cm. These points are used as the basis for the generation of digital terrain model (DTM) and digital surface model (DSM) with a resolution of 5 m.

4.1.2.1. Photosynthetic vegetation fraction cover

Photosynthetic vegetation fraction cover was estimated throughout the Automated Monte Carlo Unmixing (AutoMCU) model. This model integrates spectral mixture analysis and spectral end-member libraries resulting from fieldwork (ground spectrometer) and high-resolution hyperspectral information of Hyperion Sensor, in order to separate photosynthetic vegetation, non-photosynthetic vegetation, and bare substrate. The photosynthetic vegetation fraction cover was calculated by CLASlite v3.2 software (Figure 5) [140].

Figure 5.

Photosynthetic vegetation fraction cover.

4.1.2.2. Canopy height model

The calculation of canopy height model used altitude values of different digital terrain model (DTM) and digital surface model (DSM) in order to extract differences between both models (Figure 6).

Figure 6.

Digital canopy height model.

4.1.2.3. Correlation and autocorrelation coefficients

Correlation analysis based on multiple regressions explored statistical relationships between aboveground biomass carbon and satellite indices. The sampling points were randomly divided into 50% for model calibration and 50% for model verification. Residuals from regressions were retained and their spatial autocorrelation was quantified [141]. Moran’s I index was used to identify the type and intensity of spatial pattern, measuring the degree of autocorrelation or dependence of a distribution. Moran’s I index can be written as in the equation (9):

$I=\frac{n}{{{\displaystyle \sum }}_i{{\displaystyle \sum }}_j{\varpi }_{ij}}\text{*}\frac{{{\displaystyle \sum }}_i{{\displaystyle \sum }}_j{\varpi }_{ij}\left(x_i-\overline{x}\right)\left(x_j-\overline{x}\right)}{{{\displaystyle \sum }}_i{\left(x_i-\overline{x}\right)}^2}$E9

where n is the number of spatial units indexed by i and j, x is the variable of interest, and ${\varpi }_{ij}$ is the spatial weight matrix. Moran’s I values near to 1 indicate clustering negative values near to −1 represent spatial dispersion, and a value of zero indicates randomness. Statistical significance was expressed in terms of the Z descriptor and confidence level 1−α. Construction of the spatial weight matrix was distance based, since the spatial representation units are points. This distance or spatial lag includes at least 12 samples as recommended by Isaaks and Srivsatava [33].

4.1.3. Regression models

In addition to multiple regression, the present study compared models derived from regression-cokriging. The regression-cokriging was calculated through the estimation of a simple linear regression approach between aboveground carbon and canopy height model (equation (10)) and the addition of interpolated layer via ordinary cokriging integrated by regression residuals and the secondary variable (photosynthetic vegetation fraction cover) [31, 142]. The predictions were carried out separately for drift and residuals, and were added together later as in the following equation (11):

$\widehat{z}RCoK\left(S_0\right)=\widehat{m}\left(S_0\right)+\widehat{e}\left(S_0\right)$E10

$\widehat{z}RCoK\left(S_0\right)={\displaystyle \sum }_{k=0}^p{\widehat{\beta }}_k\text{*}{q}_k\left(S_0\right)+\left({\displaystyle \sum }_{i=1}^n{\varpi }_{i}e\left(So\right)+{\displaystyle \sum }_{j=1}^n{\varpi }_{2j}{Z}_2\left(S_i\right)\right)$E11

where βk are the coefficients of the drift model, qk is the number of auxiliary variables, ϖi(S₀) and ${\varpi }_{2j}$ are the weights determined by covariogram for regression residuals, and secondary variable and e are the regression residuals [31].

4.1.3.1. Model accuracy

The regression models were validated with data from the field sampling [143]. The RMSE criterion was used to determine which regression models have more precision in the estimation of stored carbon in the area, the RMSE can be written as equation (12):

$\text{RMSE}=\sqrt{{\displaystyle \sum }_{i=1}^n\frac{{\left(Z_{\left(si\right)}-z_{\left(si\right)}\right)}^2}{n}}$E12

where Z(si) is the reference value, z(si) is the estimated value, and n is the total number of samples.

4.2. Results

4.2.1. Correlation

The degree of association between carbon stored and each index derived from remote sensing and multiple associations (Table 2) was the synergy between canopy height model and photosynthetic vegetation fraction cover with an r coefficient of 0.88. All these correlations were positive, indicating that, as stored carbon increases, vegetation indices increase.

Satellite indices	r	Error
Canopy Hieght Model	0.85	24.17
Photosynthetic Vegetation Fraction Cover	0.52	39.29
Multiple regression	0.88	22.26

Table 2.

Correlation coefficients above-ground carbon and remote sensing indices.

4.2.2. Moran’s I index

Once the satellite-derived index, that is the most associated with carbon storage, was identified, Moran’s I was calculated from regression residuals according to sampling sites in order to identify spatial autocorrelation and, hence, the information could be included as auxiliary information in the regression-cokriging model. A value of 0.31 (z = 2.96, α = 0.01) was obtained for the regression residuals between canopy height model and aboveground carbon. In this case, a low-positive spatial autocorrelation was present, with statistical significance, indicating that the neighboring spatial units presented near or close values and a slight trend toward plots.

4.2.3. Spatial distribution

To identify the spatial distribution of stored carbon, models were developed with the application of equations resulting from multiple regression and regression-cokriging spatial methods. By obtaining the covariograms for the theoretical fit, no significant variation was found between anisotropic and isotropic covariograms, therefore, the isotropic model was settled.

As the empirical covariogram showed strong spatial dependence, it did not present constant semivariances as a function of distance. In this case, the adjusted theoretical models ranged from 10,000 to 15,000 m, distances at which the observations were independent.

We used the best fit (according to root mean squares in prediction errors) for the ordinary cokriging interpolation to evaluate which is more effective in the predictions throughout the study area. Table 3 shows the parameters obtained.

Remote Sensing indices	model	Sill	Range	Nugget
Photosynthetic Vegetation Fraction Cover and regression residuals	Exponential	0.89	12,007	0

Table 3.

Covariogram indices

4.2.4. Model accuracy

Comparison of the models with the set of verification sites produced the RMSE in tC/ha (Table 4). The models based on regression-cokriging presented the least error. Figure 5 shows the spatial distribution obtained by regression-cokriging and multiple regressions for stored carbon. The delineation of forest types (fir and pine) was based on the map of land use and vegetation of PAOT [144].

Models	RMSE (tC/ha)
Multiple regression	34.1
Regression-Cokriging	28.6

Table 4.

RMSE

Figure 7.

Spatial distribution of the aboveground carbon storage in Environmental Conservation Zone, estimated by two regression models.

The results of this research indicated that consideration of spatial autocorrelation can improve estimates of carbon content in aboveground biomass, specifically using the regression cokriging method. This could be due to its sensitivity to local variations [142], since it is particularly developed to consider the adjustment of the spatial variance model in order to improve predictions obtained from global models. The present work uses spatial modeling of indices related to carbon for the purpose of exploring the spatial autocorrelation of auxiliary variables, which is useful for representing the way in which a phenomenon radiates through spatial units. Although the primary method used for the estimation of carbon stocks in Mexico is the stratify-and-multiply approach, which assigns a single value or a range of values to each vegetation type and then multiplies these values by the areas covered by the vegetation to estimate total carbon stock values, this investigation has demonstrated a more accurate, spatial-explicit, repeatable alternative.

According to Ref. [145], autocorrelation is “perhaps, after the average and variance, the most important property of any geographic variable and, unlike these, it is explicitly linked with spatial patterns.” The present comparative analysis demonstrates the importance of the use of spatial methods to model carbon stored in the aboveground biomass, since these methods consider the spatial pattern of the data. The hypothesis of the homogeneity of the relationships between stored carbon and remote sensing indices sometimes does not consider the spatial heterogeneity of many factors affecting this relationship, such as geographic differences in orientation and climatic and soil conditions [28].

This analysis provided a synoptic mapping of aboveground biomass as a potentially valuable tool for environmental protection policies in the ECZ of Mexico City, one of the most important ecological reserves for the inhabitants of the Mexico Valley in the economic, cultural, and social sense, as well as for the volume and quality of the environmental services it provides.

5. Conclusion

Remote sensing-derived indices play a major role in forest monitoring, because traditional methodologies derive their estimates of carbon content in the biomass through forest inventories and, for its implementation, they require much time and money and are generally limited to 10-year intervals. The information resulting from them is designed to present average timber volumes linked to administrative regions, which not represent the spatial variability and therefore it generate a bias in carbon measures..

This has led to great interest to estimate, map, and monitor the carbon stored in forests more precisely, enhancing the recognition of their role in the global carbon cycle, particularly in the mitigation of greenhouse gases. Through the estimation of carbon content, a base line for calculating the dynamics of this gas as a mitigation strategy can be established.

An overview of modeling options and remote sensing resources that have been used for monitoring and researching the forest biomass was presented. Techniques ranging from collecting georeferenced data in the field to the information extraction methods from satellite images and synergies between remote sensing and geostatistics were described. A case study was selected to illustrate the application of some of these techniques in the modeling of spatial distribution of aboveground carbon in Mexican coniferous forests.

Based on the study, the following conclusions can be drawn:

1. According to these results, the synergy between remote sensing and geostatistics has the potential to estimate forest biomass to improve estimations using remote sensing indices as spatial secondary variables.

2. Geospatial methods have a better modeling adjustment (e.g., RMSE) than conventional statistical methods as multiple regressions, because geospatial methods considerer local spatial variations.

3. Best Pearson coefficient between the two variables tested in the study is the digital canopy model, which resulted from LiDAR data. This kind of information is very expensive, so integrating multispectral information can be a way to capitalize on multitemporal study of biomass.