**2. Materials and methods**

298 Remote Sensing of Biomass – Principles and Applications

improve the health of forests (Eriksson et al., 2002; Raison, 2002). Environmental benefits are also obtained in the application phase. Considering only the combustion process, the energy produced is almost neutral in terms of carbon dioxide, as the amount of CO2 emitted to the atmosphere due to FRB combustion is the same than the fixed in the FRB formation (net primary production). In addition, sulphur and chlorine emissions are very low. Consequently, taking into account only the combustion process, the generation of energy from FRB does not contribute to an increase in greenhouse gasses (Hakkila & Parikka, 2002; IDAE, 2005a). The social and economic benefits can be separated into two categories: (i) those that occur at the national scale, associated with a reduction in the use of nonrenewable fossil fuels (Domínguez, 2002; Eriksson et al., 2002; IDAE, 2005a), and (ii) those that occur at the local scale, as the utilization of a waste product leads to increasing harvests, transportation, and utilization of forest residues in power stations, which also leads to an increase in rural employment. These considerations are important for rural areas in which the level of unemployment and depopulation is a public policy issue, and where increased employment can help to support a population that can maintain the natural environment

Despite the benefits outlined above, previous studies conducted before the subsequent revision of PPRE, the Plan to Renewable Energies 2005–2010 (PRE), showed that the use of FRB had a long way to go in achieving the anticipated objectives. In the PRE analysis, the lack of specific methodologies to assess the regional-scale quantity of FRB was identified as one of the main problems (IDAE, 2005a). This is a fundamental concern because power stations that use FRB require knowledge of the amount of resource available (Esteban et al., 2004; IDAE, 2005b, 2007; Velázquez, 2006). To overcome this problem, a methodology is needed that quantifies the potential of this forest resource in a given area. This is an important matter considering that the new Plan to Renewable Energies in Spain 2011-2020 (currently being written using the provisions of Directive 2009/28/EC) set a final energy consumption target of 20% derived from renewable energy, doubling the previous mark in

Several studies have demonstrated the effectiveness of satellite images in estimating forest variables, including biomass. These studies were carried out using both passive (optical) and active (radar and lidar) sensors. Focusing on optical sensors, different experiments have been conducted using multispectral or hyperspectral coarse spatial resolution data (i.e. Anaya et al., 2009; Muukkonen & Heiskanen, 2007) and fine spatial resolution data (i.e. Gonzalez et al., 2010; Proisy et al., 2007); although the most frequent experiments have been performed using medium spatial resolution data, mainly using Landsat Thematic Mapper (TM) or Enhanced Thematic Mapper Plus (ETM+) data (Lu, 2006; Powell et al., 2010). Radar data has also been used extensively, establishing significant correlations between biomass and radar backscatter at C-band (Kurvonen et al., 2002) and L-band (Austin et al., 2003; Lucas et al., 2010), with the latter being more consistent and robust. Recently, forest biomass assessment research has focused on using polarimetric synthetic aperture radar interferometry (PolInSAR) and laser scanning systems (lidar) data, mostly by means of airborne sensors. Koch (2010) offers a

Images derived from remote sensing register continuous and complete information across a landscape and such images can be obtained at frequent intervals. These characteristics help to overcome some of the problems associated with inventory methods exclusively based on field work, interpolation techniques and GIS (Franklin, 2001; Lu, 2006; Salvador & Pons, 1998a).

complete review of the recent advances and future developments in this issue.

(Borsboom et al., 2002; Domínguez, 2002; IDAE, 2005a, 2007).

the 2010 Plan (MITYC, 2010).

The methodology applied in this study was developed in the framework of geographic information technologies (Remote Sensing and Geographical Information Systems, GIS) as information sources and tools for forest management. The software Erdas Imagine was used to process the Landsat images, ARC-GIS was used for the management of ancillary information, and SPSS was employed for statistical analyses.

The methodology proposed and developed in this study was divided into five phases or steps. Firstly, FRB data was obtained by means of field work and pre-existing forest data. Secondly, suitable satellite images were selected and different treatments were applied in order to convert the data to a suitable format for quantitative analysis and to guarantee the validity of the results. Thirdly, three different methods were used to relate field data and radiometric information in order to overcome the difficulties involved in estimating forest parameters in heterogeneous Mediterranean forests. The fourth phase focused on developing and running regression models to map FRB in the study area, and to evaluate and compare the results at pixel level obtained using the applied extraction methods. Finally, the best model was selected for application to a recent image to quantify the amount of FRB at present time in the study area.

#### **2.1 Study area**

The study area is the Province of Teruel, which is located in the northeastern Iberian Peninsula (Figure 1). This province was selected in the context of the LIGNOSTRUM project

Using Remote Sensing to Estimate a Renewable Resource: Forest Residual Biomass 301

At the start of the LIGNOSTRUM project in 2002 the most recent, accurate, and complete information on the dimensional parameters of trees in the forests of the study area was found in the Second National Forest Inventory (NFI-2), which was completed in 1994 (MMA, 1996). The Spanish NFI is carried out every 10 years as a systematic field sampling of permanent plots. These plots are located in the corners of the UTM 1x1 km grid of the 1:50,000 National Topographic Cartography that are within forested areas. The placement of plots in the field was performed using georeferenced 1:30,000 aerial photographs and topographic cartography. Plots have a circular shape with radii ranging from 5 m to 25 m depending on the DBH of the trees; only trees with a DBH greater than 7.5 cm were considered. The NFI-2 fieldwork in Teruel Province was undertaken from March to August 1994, sampling a total of 2083 plots. The DBH (two perpendicular diameters measured at a height of 1.30 m) and total height (accurate to within half a meter) were measured for all trees. A detailed description of NFI-2 is provided in the summary report of the inventory

Taking into account these characteristics of NFI-2, we devised a species-based sampling strategy to obtain FBR regressions. We located all the sample areas in forests managed by the Teruel Environment Service owing to the destructive nature of this sampling. Fieldwork was carried out from November 2003 to June 2004. For each species, the distribution of sampling points was proportional to the number of trees, the volume of wood, and the stand basal area documented in NFI-2 data. We considered a normal diameter range from 7.5 to 40 cm because these are the minimum diameters recorded in NFI-2 data and they are

In the field, we used a scale with an accuracy of 250 grams to obtain the wet weight of the residual biomass. In addition, we took two size measurements for each sampled tree: DBH (accuracy of 5 mm) and height (centimeter accuracy). Finally, to avoid the influence of the variations in the degree of wetness in the samples, we collected samples of leaves and branches in the different trees in order to calculate the dry weight by applying the method

The total number of sampled trees was 186 (30 of *Pinus sylvestris*, 59 of *Pinus halepensis*, 57 of *Pinus nigra*, and 40 of *Pinus pinaster*). We obtained a regression equation for each species, but two in the case of *Pinus pinaster* regarding the tree origin (natural or reforested) (Table 1).

**Models** *R2*

*P. sylvestris 30 FRB=0.064·DBH3.3/H1.5* 0.974 12.298 - *P. nigra 57 FRB=338.416·e-35.116/DBH* 0.910 18.836 1/*DBH P. halepensis 59 FRB=0.067 DBH3/H* 0.969 13.063 1/(*DBH*\*(1+ORI))

*natural 28 FRB=1.101·10-3·DBH4/H* 0.973 6.004 -

height, H is tree height, and ORI is the tree origin (0=reforested; 1=natural)

Table 1. Regression equations obtained by species, R2 adjusted, standard error associated and weight expression to accomplish the assumption of residual homocedasticity. Source: Alonso et al., (2005), where FRB is forest residual biomass, DBH is tree diameter at breast

*a*

*4·DBH3.823·H0.337* 0.974 12.175 -

**Std. Desv. (kg)** 

**Weight (X-k)** 

the most frequently observed in the tables of production of the study area.

(MMA, 1996).

described by Joosten et al. (2004).

**Nº of sampled trees** 

*reforested 12 FRB=1.97·10-*

**Species** 

*P. pinaster* 

*P. pinaster* 

(2002-2005), of which the present research was a component. The aim of the Spanish government-financed LIGNOSTRUM project was to increase the use of agricultural and forest residues as energy resources. There were two main reasons for selecting this territory as a setting to measure the benefits of energy use of forest residues: (i) the presence of large agricultural and forest areas; and (ii) the existence of economically disadvantaged rural areas.

Fig. 1. Location of the study area

There are three characteristics of this province that mean difficulties to achieve the objective of estimate FRB using satellite images: (i) the size of the province (14,804 km2), (ii) the concentration of forest resources in the mountain areas (as the topography modifies the reflectance values of equal ground covers), and (iii) the Mediterranean characteristics of these forests (presence of multiple species with a complex spatial structure resulting from multiple interrelationships between different ecological factors and a history of human activity extending back for many centuries). These peculiarities determined the types of images that were used, the treatments applied to the images and the methods applied to extract the radiometric data. Among all types of forests present in the study area, only pine forests were selected for study since they represent slightly over 71% of the total forest area of the province and they have the greatest potential for FRB generation from harvest (IDAE, 2007; MMA, 1996).

#### **2.2 FRB data**

Aboveground biomass (AGB) is normally estimated using allometric equations that relate it to tree data, usually diameter at breast height (DBH) and height (Ketterings et al., 2001; Montero et al., 2005; Pilli et al., 2006; Wang, 2006; Zianis & Mencuccini, 2004). These equations show a large variation related to tree species, geographic location (Ketterings et al., 2001; Montero et al., 2005), and other factors such as climate and tree population age (Zianis & Mencuccini, 2004). Few biomass equations have been developed for Spain (Montero et al., 2005), with even fewer focusing on residual biomass (in some cases, crown biomass has been assessed). Consequently, it was necessary to develop specific FRB regressions for each pine species in the study area for application to pre-existing tree dimensional data.

(2002-2005), of which the present research was a component. The aim of the Spanish government-financed LIGNOSTRUM project was to increase the use of agricultural and forest residues as energy resources. There were two main reasons for selecting this territory as a setting to measure the benefits of energy use of forest residues: (i) the presence of large agricultural and forest areas; and (ii) the existence of economically disadvantaged rural

There are three characteristics of this province that mean difficulties to achieve the objective of estimate FRB using satellite images: (i) the size of the province (14,804 km2), (ii) the concentration of forest resources in the mountain areas (as the topography modifies the reflectance values of equal ground covers), and (iii) the Mediterranean characteristics of these forests (presence of multiple species with a complex spatial structure resulting from multiple interrelationships between different ecological factors and a history of human activity extending back for many centuries). These peculiarities determined the types of images that were used, the treatments applied to the images and the methods applied to extract the radiometric data. Among all types of forests present in the study area, only pine forests were selected for study since they represent slightly over 71% of the total forest area of the province and they have the greatest potential for FRB generation from harvest (IDAE,

Aboveground biomass (AGB) is normally estimated using allometric equations that relate it to tree data, usually diameter at breast height (DBH) and height (Ketterings et al., 2001; Montero et al., 2005; Pilli et al., 2006; Wang, 2006; Zianis & Mencuccini, 2004). These equations show a large variation related to tree species, geographic location (Ketterings et al., 2001; Montero et al., 2005), and other factors such as climate and tree population age (Zianis & Mencuccini, 2004). Few biomass equations have been developed for Spain (Montero et al., 2005), with even fewer focusing on residual biomass (in some cases, crown biomass has been assessed). Consequently, it was necessary to develop specific FRB regressions for each pine species in the study area for application to pre-existing tree

areas.

Fig. 1. Location of the study area

2007; MMA, 1996).

dimensional data.

**2.2 FRB data** 

At the start of the LIGNOSTRUM project in 2002 the most recent, accurate, and complete information on the dimensional parameters of trees in the forests of the study area was found in the Second National Forest Inventory (NFI-2), which was completed in 1994 (MMA, 1996). The Spanish NFI is carried out every 10 years as a systematic field sampling of permanent plots. These plots are located in the corners of the UTM 1x1 km grid of the 1:50,000 National Topographic Cartography that are within forested areas. The placement of plots in the field was performed using georeferenced 1:30,000 aerial photographs and topographic cartography. Plots have a circular shape with radii ranging from 5 m to 25 m depending on the DBH of the trees; only trees with a DBH greater than 7.5 cm were considered. The NFI-2 fieldwork in Teruel Province was undertaken from March to August 1994, sampling a total of 2083 plots. The DBH (two perpendicular diameters measured at a height of 1.30 m) and total height (accurate to within half a meter) were measured for all trees. A detailed description of NFI-2 is provided in the summary report of the inventory (MMA, 1996).

Taking into account these characteristics of NFI-2, we devised a species-based sampling strategy to obtain FBR regressions. We located all the sample areas in forests managed by the Teruel Environment Service owing to the destructive nature of this sampling. Fieldwork was carried out from November 2003 to June 2004. For each species, the distribution of sampling points was proportional to the number of trees, the volume of wood, and the stand basal area documented in NFI-2 data. We considered a normal diameter range from 7.5 to 40 cm because these are the minimum diameters recorded in NFI-2 data and they are the most frequently observed in the tables of production of the study area.

In the field, we used a scale with an accuracy of 250 grams to obtain the wet weight of the residual biomass. In addition, we took two size measurements for each sampled tree: DBH (accuracy of 5 mm) and height (centimeter accuracy). Finally, to avoid the influence of the variations in the degree of wetness in the samples, we collected samples of leaves and branches in the different trees in order to calculate the dry weight by applying the method described by Joosten et al. (2004).

The total number of sampled trees was 186 (30 of *Pinus sylvestris*, 59 of *Pinus halepensis*, 57 of *Pinus nigra*, and 40 of *Pinus pinaster*). We obtained a regression equation for each species, but two in the case of *Pinus pinaster* regarding the tree origin (natural or reforested) (Table 1).


Table 1. Regression equations obtained by species, R2 adjusted, standard error associated and weight expression to accomplish the assumption of residual homocedasticity. Source: Alonso et al., (2005), where FRB is forest residual biomass, DBH is tree diameter at breast height, H is tree height, and ORI is the tree origin (0=reforested; 1=natural)

Using Remote Sensing to Estimate a Renewable Resource: Forest Residual Biomass 303

effects in the surface reflectivity registered in both Landsat images. In this process, the relative data on surface reflectivity given in digital numbers (DN) is transformed into an absolute form in reflectance values. In areas of rough terrain, as in the present study area, a good topographic normalization is required to compensate for the varying solar illumination associated with the irregular shape of the terrain. To accomplish this, first, the default transmittance method, proposed in Chavez (1996), was applied to eliminate the atmospheric effects present in both Landsat images. Second, conversion from DN values to reflectance measurements was carried out using the Minnaert Correction Method (Colby, 1991). Different transformations (determine in previous work focusing on forest parameter estimation) were applied to the 1994 image to increase the amount of spectral information, including different vegetation indices. Figure 2 shows the location of the 482 NFI-2 plots

Fig. 2. Location of the 482 NFI-2 plots over a false colour band combination 7/4/3 (RGB) of

As mentioned previously, heterogeneity is an inherent characteristic of Mediterranean forests and is one of the main factors that complicates forest parameter estimation in this area using remote sensing images. This is mainly because spatial heterogeneity is registered in the satellite image, resulting in a large spectral variability in the spectral response of the areas covered by these forests. This introduces difficulties in building accurate predictive models. As a result, two NFI-2 plots with equal amounts of FRB can show different reflectance values because of the presence of other landscape elements (i.e. firebreaks, trails), the position of the FRB plot on the border of another land cover type (i.e. scrublands, farmlands) or because of high variability within the forested area (presence of different tree species, ages). In addition, other problems detected in previous research in Mediterranean

**2.4 Linking FRB ground data to 1994 radiometric data** 

the study area

over a false colour band combination 7/4/3 (RGB) of the study area.

In applying the obtained FRB equations to the NFI-2 plots, we only chose those in which pines were among the range used in the equations, thus obtaining 617 plots. These residual biomass data were linked to a table with all of the other NFI-2 variables, including their locations in UTM coordinates. To avoid complexity associated with the mixture of tree species in the spectral data of the Landsat image, we selected only monospecies pine plots (Salvador & Pons, 1998b). In addition, we removed from the dataset plots located in an area affected by a wildfire that occurred in 1994 and those affected by cloud shadows in the southeast corner of the Landsat image. Following this selection process, a total of 482 plots were finally obtained. We created a point map with these plots to read the Landsat TM spectral information. The estimated residual biomass in selected plots ranged from 0.107 to 64.720 tons/ha. This large range is related to the high degree of heterogeneity of Mediterranean forests.

#### **2.3 Image pre-processing and development of new spectral indices**

Among the different types of remote sensing data available to achieve the objective of this research, Landsat images were selected because they are one of the most common in forestry-related applications and estimates of aboveground biomass (AGB) at regional-local scales (i.e. Fazakas et al., 1999; Foody et al., 2001, 2003; Gasparri et al., 2010; Hall et al., 2006; Labrecque et al., 2003, 2006; Lu, 2005; Lu & Batistiella, 2005; Lu et al., 2004; Mäkelä & Pekkarinen, 2004; Mallinis et al., 2004; Meng et al., 2009; Powell et al., 2010; Roy & Ravan, 1996; Salvador & Pons, 1998a, 1998b; Steininger, 2000; Tangki & Chappell, 2008; Wulder et al., 2008; Zheng et al., 2004, 2007). In addition, taking into account the research objectives, there were two other important reasons for using Landsat images. Firstly, essentially all of the study area was represented in one scene, an advantage over other approaches in which errors are introduced owing to inconsistencies in radiometric corrections, necessary when multiple scenes are used (Chuvieco, 2002; Foody et al., 2003). Secondly, the methodology would be useful and replicable in other areas with similar characteristics because of the systematic and frequent coverage of Landsat images over time as well as their ease of distribution at no cost to users via the Internet (Chander et al., 2009).

The study used two Landsat 5 TM images recorded on 10 June 1994 and 5 July 2008 (Path 199/Row 32). The earlier image was selected on the basis of temporal coincidence with NFI-2 field work, so it was useful in performing the methodology, whereas the later image was chosen for a recent FRB estimation and inventory in the study area.

We applied three pre-processing techniques to both images to obtain radiometric variables consistent and suitable to be related to field data, and to guarantee the applicability of the estimation model carried out with the 1994 image to the 2008 image: (i) geometric correction, (ii) radiometric correction, and (iii) transformations and ratios. The importance of accurate geometric rectification in the present research arises from the fact that the image data were linked with ground data. Both Landsat TM images were geometrically rectified into a local UTM projection using a second-order polynomial model. The rectification model incorporated a Digital Elevation Model (DEM) with 25 m spatial resolution. Ground Control Points (GCPs) were taken from high-resolution ortho-photographs (pixel size = 1 x 1 m). A total of 125 GCPs were used to re-project each image with an estimated Root Mean Squared Error (RMSE) of 0.60 and 0.48 pixels, respectively. A nearest-neighbor resampling technique was used to minimize changes in the radiometric values of the ground data, with the pixel re-projected to 25 m. The radiometric correction addresses atmospheric and topographic

In applying the obtained FRB equations to the NFI-2 plots, we only chose those in which pines were among the range used in the equations, thus obtaining 617 plots. These residual biomass data were linked to a table with all of the other NFI-2 variables, including their locations in UTM coordinates. To avoid complexity associated with the mixture of tree species in the spectral data of the Landsat image, we selected only monospecies pine plots (Salvador & Pons, 1998b). In addition, we removed from the dataset plots located in an area affected by a wildfire that occurred in 1994 and those affected by cloud shadows in the southeast corner of the Landsat image. Following this selection process, a total of 482 plots were finally obtained. We created a point map with these plots to read the Landsat TM spectral information. The estimated residual biomass in selected plots ranged from 0.107 to 64.720 tons/ha. This large range is related to the high degree of heterogeneity of

Among the different types of remote sensing data available to achieve the objective of this research, Landsat images were selected because they are one of the most common in forestry-related applications and estimates of aboveground biomass (AGB) at regional-local scales (i.e. Fazakas et al., 1999; Foody et al., 2001, 2003; Gasparri et al., 2010; Hall et al., 2006; Labrecque et al., 2003, 2006; Lu, 2005; Lu & Batistiella, 2005; Lu et al., 2004; Mäkelä & Pekkarinen, 2004; Mallinis et al., 2004; Meng et al., 2009; Powell et al., 2010; Roy & Ravan, 1996; Salvador & Pons, 1998a, 1998b; Steininger, 2000; Tangki & Chappell, 2008; Wulder et al., 2008; Zheng et al., 2004, 2007). In addition, taking into account the research objectives, there were two other important reasons for using Landsat images. Firstly, essentially all of the study area was represented in one scene, an advantage over other approaches in which errors are introduced owing to inconsistencies in radiometric corrections, necessary when multiple scenes are used (Chuvieco, 2002; Foody et al., 2003). Secondly, the methodology would be useful and replicable in other areas with similar characteristics because of the systematic and frequent coverage of Landsat images over time as well as their ease of

The study used two Landsat 5 TM images recorded on 10 June 1994 and 5 July 2008 (Path 199/Row 32). The earlier image was selected on the basis of temporal coincidence with NFI-2 field work, so it was useful in performing the methodology, whereas the later image was

We applied three pre-processing techniques to both images to obtain radiometric variables consistent and suitable to be related to field data, and to guarantee the applicability of the estimation model carried out with the 1994 image to the 2008 image: (i) geometric correction, (ii) radiometric correction, and (iii) transformations and ratios. The importance of accurate geometric rectification in the present research arises from the fact that the image data were linked with ground data. Both Landsat TM images were geometrically rectified into a local UTM projection using a second-order polynomial model. The rectification model incorporated a Digital Elevation Model (DEM) with 25 m spatial resolution. Ground Control Points (GCPs) were taken from high-resolution ortho-photographs (pixel size = 1 x 1 m). A total of 125 GCPs were used to re-project each image with an estimated Root Mean Squared Error (RMSE) of 0.60 and 0.48 pixels, respectively. A nearest-neighbor resampling technique was used to minimize changes in the radiometric values of the ground data, with the pixel re-projected to 25 m. The radiometric correction addresses atmospheric and topographic

**2.3 Image pre-processing and development of new spectral indices** 

distribution at no cost to users via the Internet (Chander et al., 2009).

chosen for a recent FRB estimation and inventory in the study area.

Mediterranean forests.

effects in the surface reflectivity registered in both Landsat images. In this process, the relative data on surface reflectivity given in digital numbers (DN) is transformed into an absolute form in reflectance values. In areas of rough terrain, as in the present study area, a good topographic normalization is required to compensate for the varying solar illumination associated with the irregular shape of the terrain. To accomplish this, first, the default transmittance method, proposed in Chavez (1996), was applied to eliminate the atmospheric effects present in both Landsat images. Second, conversion from DN values to reflectance measurements was carried out using the Minnaert Correction Method (Colby, 1991). Different transformations (determine in previous work focusing on forest parameter estimation) were applied to the 1994 image to increase the amount of spectral information, including different vegetation indices. Figure 2 shows the location of the 482 NFI-2 plots over a false colour band combination 7/4/3 (RGB) of the study area.

Fig. 2. Location of the 482 NFI-2 plots over a false colour band combination 7/4/3 (RGB) of the study area

#### **2.4 Linking FRB ground data to 1994 radiometric data**

As mentioned previously, heterogeneity is an inherent characteristic of Mediterranean forests and is one of the main factors that complicates forest parameter estimation in this area using remote sensing images. This is mainly because spatial heterogeneity is registered in the satellite image, resulting in a large spectral variability in the spectral response of the areas covered by these forests. This introduces difficulties in building accurate predictive models. As a result, two NFI-2 plots with equal amounts of FRB can show different reflectance values because of the presence of other landscape elements (i.e. firebreaks, trails), the position of the FRB plot on the border of another land cover type (i.e. scrublands, farmlands) or because of high variability within the forested area (presence of different tree species, ages). In addition, other problems detected in previous research in Mediterranean

Using Remote Sensing to Estimate a Renewable Resource: Forest Residual Biomass 305

method avoids some problems related to the use of fixed windows, as they do not completely eliminate errors in image registration and location of the sample plots, may intersect several stands with different spectral and forest characteristics, and not fit well with the limited spatial and spectral resolution of Landsat data (Mäkelä & Pekkarinen, 2001; Mäkelä & Pekkarinen,

FRB cartography with stands larger than NFI-2 plots did not exist in the study area; however, the interpretation of high-resolution aerial photographs has shown to be useful in forestinventory applications such as stratification, volume estimation, and measurements of forest characteristics (Lu, 2006; Mäkelä & Pekkarinen, 2004). In this context, the high-resolution ortho-photographs used in geometric rectification were applied to extend the plot areas to larger sizes with visually similar composition and forest structure. This was performed using a "heads-up," on-screen digitizing technique within a GIS application. The selected NFI-2 plots, with their respective radii size, were displayed over the composite digital aerial photograph with 1 m spatial resolution. Where possible, larger homogeneous areas containing the in situ plots were then defined. Aragon 1:50,000 digital forest cartography was used to guarantee the pure composition of the new areas. To avoid errors in the results, only homogeneous areas with a high degree of similarity between observations from the aerial photographs of the NFI-2 plots were selected. As a result, a total of 131 areas were selected. In addition, when extracting mean values, only those pixels located in the core areas of the delimited stand were selected. Pixels located on the borders were avoided, as these can reflect properties of

landscape elements in the immediate vicinity of the delimited stands (Figure 4).

Fig. 3. (a) 3 x 3 pixel window centred on a plot with high heterogeneity (high CV); (b) 3 x 3

pixel window centred on a plot with low heterogeneity (low CV)

2004; Muukkonen & Heiskanen, 2005; Pekkarinen, 2002).

forests are inaccuracies in the localization of inventory field plots, small plot sizes, and the small number of plots used in the analysis (Mallinis et al., 2004; Maselli & Chiesi, 2006; Salvador & Pons, 1998a, 1998b; Shoshany, 2000; Vázquez de la Cueva, 2005).

The present study tested three different methods to extract the radiometric data in order to overcome the problems outlined above and achieve accurate FRB regression models: (i) fixed pixel windows or kernels, (ii) visual analysis, and (iii) spectral segmentation.

#### **2.4.1 Fixed pixel windows**

The use of kernels or fixed pixel windows greater than one pixel is the most common method to extract radiometric information in works focused on forest parameter estimation by means of satellite images using field plots data (i.e. Foody et al., 2001; Hall et al., 2006; Labrecque et al., 2003; Lu et al., 2004; Lu, 2005; and Lu & Batistiella, 2005 Roy & Ravan, 1996). In our case, taking into account the high variability of Mediterranean forests, a 3 x 3 pixel window centered on each plot was used to extract the mean value and standard deviation for all radiometric variables derived from the 1994 image. Following the methodology used in Labrecque et al. (2003), the variability in each plot was calculated over the six spectral reflectance bands using Pearson's coefficient of variation (CV) (1):

$$\text{CV} = \frac{\text{S}}{\left| \overline{\overline{\text{X}}} \right|} \tag{1}$$

where S is the standard deviation calculated for the 3 x 3 pixel window centered on each plot and X is the average calculated for the same window.

Thus, information is obtained regarding the degree of spectral heterogeneity within the immediate surroundings of each plot: plots with high CV have a high degree of spatial heterogeneity, while plots with low CV are the most homogeneous. With the aim of determining the influence of spectral heterogeneity on estimates of FRB, the CV values were used to divide the sample into 10 homogeneity levels. Initially, we calculated the thresholds needed to delimit the 10 CV percentiles in each reflectance band and, after that, these thresholds were used to create 10 groups of plots as follows: the first group contained all of the plots, the second contained only those plots whose CV values were lower than that of percentile 9 in all TM bands, the third contained only those plots whose CV values were lower than that of percentile 8 in all of the TM bands, and so on, until the 10 groups were delineated. As a result, the first of the delineated groups contains the plots with the highest degree of spectral heterogeneity, while the last contains the plots that are the most homogeneous. For the rest of this chapter, these plot groups are named using the CV percentile employed in creating them. Figure 3 shows the differences regarding the degree of spectral heterogeneity within the immediate vicinity among one plot classified in group CV9 and another classified in group CV4.

#### **2.4.2 Visual analysis**

In the second method to extract the radiometric data from Landsat image, we use visual image analysis to delimit homogeneous FRB units larger than the NFI-2 plots. Different studies have demonstrated that biomass can be estimated more accurately if the ground data is referenced over areas larger than the pixel level using similar units in forest and spectral characteristics (Hyyppä & Hyyppä 2001; Muukkonen & Heiskanen, 2005; Zheng, 2004). In addition, this

forests are inaccuracies in the localization of inventory field plots, small plot sizes, and the small number of plots used in the analysis (Mallinis et al., 2004; Maselli & Chiesi, 2006;

The present study tested three different methods to extract the radiometric data in order to overcome the problems outlined above and achieve accurate FRB regression models: (i)

The use of kernels or fixed pixel windows greater than one pixel is the most common method to extract radiometric information in works focused on forest parameter estimation by means of satellite images using field plots data (i.e. Foody et al., 2001; Hall et al., 2006; Labrecque et al., 2003; Lu et al., 2004; Lu, 2005; and Lu & Batistiella, 2005 Roy & Ravan, 1996). In our case, taking into account the high variability of Mediterranean forests, a 3 x 3 pixel window centered on each plot was used to extract the mean value and standard deviation for all radiometric variables derived from the 1994 image. Following the methodology used in Labrecque et al. (2003), the variability in each plot was calculated over

Salvador & Pons, 1998a, 1998b; Shoshany, 2000; Vázquez de la Cueva, 2005).

**2.4.1 Fixed pixel windows** 

fixed pixel windows or kernels, (ii) visual analysis, and (iii) spectral segmentation.

the six spectral reflectance bands using Pearson's coefficient of variation (CV) (1):

plot and X is the average calculated for the same window.

CV9 and another classified in group CV4.

**2.4.2 Visual analysis** 

<sup>S</sup> CV

where S is the standard deviation calculated for the 3 x 3 pixel window centered on each

Thus, information is obtained regarding the degree of spectral heterogeneity within the immediate surroundings of each plot: plots with high CV have a high degree of spatial heterogeneity, while plots with low CV are the most homogeneous. With the aim of determining the influence of spectral heterogeneity on estimates of FRB, the CV values were used to divide the sample into 10 homogeneity levels. Initially, we calculated the thresholds needed to delimit the 10 CV percentiles in each reflectance band and, after that, these thresholds were used to create 10 groups of plots as follows: the first group contained all of the plots, the second contained only those plots whose CV values were lower than that of percentile 9 in all TM bands, the third contained only those plots whose CV values were lower than that of percentile 8 in all of the TM bands, and so on, until the 10 groups were delineated. As a result, the first of the delineated groups contains the plots with the highest degree of spectral heterogeneity, while the last contains the plots that are the most homogeneous. For the rest of this chapter, these plot groups are named using the CV percentile employed in creating them. Figure 3 shows the differences regarding the degree of spectral heterogeneity within the immediate vicinity among one plot classified in group

In the second method to extract the radiometric data from Landsat image, we use visual image analysis to delimit homogeneous FRB units larger than the NFI-2 plots. Different studies have demonstrated that biomass can be estimated more accurately if the ground data is referenced over areas larger than the pixel level using similar units in forest and spectral characteristics (Hyyppä & Hyyppä 2001; Muukkonen & Heiskanen, 2005; Zheng, 2004). In addition, this

<sup>X</sup> (1)

method avoids some problems related to the use of fixed windows, as they do not completely eliminate errors in image registration and location of the sample plots, may intersect several stands with different spectral and forest characteristics, and not fit well with the limited spatial and spectral resolution of Landsat data (Mäkelä & Pekkarinen, 2001; Mäkelä & Pekkarinen, 2004; Muukkonen & Heiskanen, 2005; Pekkarinen, 2002).

FRB cartography with stands larger than NFI-2 plots did not exist in the study area; however, the interpretation of high-resolution aerial photographs has shown to be useful in forestinventory applications such as stratification, volume estimation, and measurements of forest characteristics (Lu, 2006; Mäkelä & Pekkarinen, 2004). In this context, the high-resolution ortho-photographs used in geometric rectification were applied to extend the plot areas to larger sizes with visually similar composition and forest structure. This was performed using a "heads-up," on-screen digitizing technique within a GIS application. The selected NFI-2 plots, with their respective radii size, were displayed over the composite digital aerial photograph with 1 m spatial resolution. Where possible, larger homogeneous areas containing the in situ plots were then defined. Aragon 1:50,000 digital forest cartography was used to guarantee the pure composition of the new areas. To avoid errors in the results, only homogeneous areas with a high degree of similarity between observations from the aerial photographs of the NFI-2 plots were selected. As a result, a total of 131 areas were selected. In addition, when extracting mean values, only those pixels located in the core areas of the delimited stand were selected. Pixels located on the borders were avoided, as these can reflect properties of landscape elements in the immediate vicinity of the delimited stands (Figure 4).

Fig. 3. (a) 3 x 3 pixel window centred on a plot with high heterogeneity (high CV); (b) 3 x 3 pixel window centred on a plot with low heterogeneity (low CV)

Using Remote Sensing to Estimate a Renewable Resource: Forest Residual Biomass 307

Fig. 5. Procedure to extract radiometric data combining the use of segmentation and 3 x 3

We used Pearson's coefficient to explore the feasibility of building accurate and representative predictive models using the plot groups derived from the three extraction methods. After this selection, each of the three groups was divided into two samples: 80% of the sample was used to carry out the predictive model and the remaining 20% was used to validate it. This sample division in each group was made randomly to guarantee the execution of the estimation equations and the validation processes. Moreover, to assess the robustness of the models, the sample division was done five times, completing the

The multiple linear regression model was performed, using the option "stepwise" to include only the significant variables (p < 0.05). In addition, performance was verified for all of the principles assumed for this type of regression at the model and variable level. To evaluate the performance of models, the coefficient of determination (R2) was used, and the Root Mean Square Error (RMSE) and the relative RMSE (RMSEr) were calculated using 20% of the sample previously reserved for the validation. Finally, the best model conducted with each of the extraction methods was applied to the 1994 Landsat data in order to obtain the FRB estimation cartography. The three derived cartographies were validated using plots not included in the groups considered in the regression models, and the RMSE and RMSEr were

The results obtained with the June 1994 Landsat image (selected on the basis of its temporal coincidence with NFI-2 fieldwork) at the regression model level and in the cartography validation (in terms of R2 and RMSE and RMSEr, respectively) were analyzed. Then, the most suitable estimation model was selected for application to the July 2008 Landsat image. As a result, current information was obtained about the potential quantity of this energy

**3.1 Models run using fixed pixel windows to extract the radiometric information**  Table 3 shows the correlation coefficients tendency of the original bands and some of the variables derived from them correlated to FRB in the first nine groups delimited using the CV. As expected, higher correlations were obtained with increasing spectral homogeneity

pixel windows with restrictions

calculated to evaluate the results.

**3. Results** 

**2.6 Model application and inventory** 

resource and its spatial distribution in the study area.

**2.5 FRB mapping and validation method** 

respective estimation model and its validation each time.

Fig. 4. (a) Homogeneous area delimited using high-resolution aerial photographs; (b) Core pixels of the homogeneous area used to extract spectral data from the Landsat TM image
